ホースで庭に水をまいているとき、ホースの先端を指で潰すと水の勢いが急に強くなります。あるいは、川幅が狭くなっている場所では、流れが速くなっていることに気づいた経験はないでしょうか。これらの身近な現象の裏にある物理法則が、本記事で扱う連続の方程式（equation of continuity）です。
連続の方程式は「流体の質量は途中で消えたり、どこからか湧き出したりしない」という当たり前の事実を数式で表現したものです。一見すると単純なこの法則ですが、流体力学のあらゆる問題を解く出発点になっています。
連続の方程式を理解すると、以下のような応用に直結します。
- 配管設計: 管路の断面積変化による流速・圧力の予測。工場のプラントや水道システム設計の基本です
- 気象シミュレーション: 大気の流れを数値的に解くとき、質量保存は必ず満たすべき拘束条件になります
- 航空工学: 翼まわりの流れ解析で、連続の方程式とナビエ・ストークス方程式を連立して揚力を計算します
- 血流のモデリング: 血管の狭窄部での流速増加を定量評価し、動脈硬化のリスク評価に活用されています
このように、連続の方程式は流体力学の根幹をなす方程式であり、ベルヌーイの定理やナビエ・ストークス方程式と並んで、流体現象を理解するための三本柱の一つです。
本記事の内容
- 質量保存則の直感的な理解
- 積分形の連続の方程式の導出（検査体積の考え方から丁寧に）
- 微分形への変換（発散定理を使ったステップごとの変換）
- 非圧縮性流体への適用と $\nabla \cdot \bm{v} = 0$ の物理的意味
- 1次元管路への適用と具体的な計算例
- Pythonによる流れ場の可視化と流線の描画
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ベクトル場の発散（div） — 発散定理を使った積分形から微分形への変換で必要です
- 流体の性質 — 密度、粘性、圧縮性など流体の基本的な性質を確認できます
- ベクトル解析 — スカラー場・ベクトル場の基本と勾配・発散・回転の定義を扱っています
質量保存則とは — 「水は消えない、湧き出さない」
連続の方程式の出発点は、日常の直感そのものです。蛇口をひねって水を出しているとき、ホースの途中で水が消えてなくなることはありません。逆に、ホースの途中から勝手に水が増えることもありません。これは「質量は保存される」という自然界の根本原理で、物理学の歴史の中で古くから認識されてきました。
流体力学における質量保存則を、もう少し厳密に言い換えると次のようになります。
ある閉じた領域に注目したとき、その領域内の質量の時間変化は、領域の境界を通って出入りする流体の質量流量だけで決まる。
つまり、領域の中で質量が勝手に生成されたり消滅したりすることはなく、質量の増減は必ず境界を通じた出入りで説明できるということです。
たとえば、バスタブにお湯をためている場面を想像してください。蛇口から入ってくるお湯の量と、排水口から出ていくお湯の量の差が、バスタブ内のお湯の増減量に等しくなります。蛇口を閉じて排水口も塞げば、バスタブの中のお湯は一定量を保ちます。この直感をそのまま数式にしたものが、連続の方程式の積分形です。
ここで大切なのは、この原理があらゆるスケールで成り立つことです。バスタブのような大きなスケールでも、ミクロな流体要素のスケールでも、質量保存は破れません。この「あらゆるスケールで成立する」という性質こそが、積分形と微分形という2つの表現を生み出す原動力になっています。
それでは、この直感を数式として定式化していきましょう。まずは「閉じた領域」を数学的に定義するところから始めます。
検査体積の考え方
連続の方程式を導出するには、流体の中に仮想的な領域を設定して、そこへの出入りを追跡する必要があります。この仮想的な領域を検査体積（control volume）と呼びます。
検査体積は、空間に固定された任意の閉じた領域 $V$ で、その表面（検査面）を $S$ で表します。大切なのは、検査体積は「物理的な壁」ではなく、あくまで分析のために頭の中で設定した領域だということです。流体は検査面を自由に通過できます。
検査体積の設定は、流体力学における解析の最も基本的なアプローチです。これは電磁気学でガウスの法則を適用するときに閉曲面を設定することと同じ発想です。ある空間領域を「切り取って」注目することで、複雑な流れの中に明確な収支関係を見出すことができます。
検査体積には2種類の量が関係します。
- 検査体積内の質量 — 領域 $V$ の中に存在する流体の質量
- 検査面を通過する質量流量 — 単位時間に表面 $S$ を通って出入りする流体の質量
この2つの量の間に成り立つ関係式が、まさに連続の方程式です。次のセクションでは、これらを数式で表現し、積分形の連続の方程式を導出します。
積分形の導出
検査体積内の質量
空間に固定された検査体積 $V$ を考えます。体積内の各点における流体の密度を $\rho(\bm{x}, t)$ とすると、検査体積内に含まれる流体の全質量 $M$ は体積積分で表されます。
$$ \begin{equation} M = \int_V \rho \, dV \end{equation} $$
この式は「密度 $\rho$ を体積全体にわたって足し合わせれば全質量が求まる」という直感そのままの表現です。密度が一様でない流れ（たとえば温度差がある水の流れ）では、場所によって $\rho$ の値が異なるため、積分が必要になります。
単位時間あたりの質量変化
検査体積内の質量が時間とともに変化する速さは、上の式を時間で微分することで得られます。
$$ \begin{equation} \frac{dM}{dt} = \frac{\partial}{\partial t}\int_V \rho \, dV \end{equation} $$
ここで時間微分を偏微分 $\partial / \partial t$ で書いているのは、検査体積 $V$ が空間に固定されている（時間によらない）ためです。この偏微分は積分の中に入れることができ、次のように書けます。
$$ \frac{dM}{dt} = \int_V \frac{\partial \rho}{\partial t} \, dV $$
この式は「検査体積内の全質量の時間変化は、各点での密度の時間変化を足し合わせたもの」であることを表しています。
表面を通る質量流量
次に、検査面 $S$ を通って流出する質量を考えます。流体の速度ベクトルを $\bm{v}$ 、表面の外向き法線ベクトルを $\bm{n}$ とします。
面の微小要素 $dS$ を単位時間に通過する流体の体積は $\bm{v} \cdot \bm{n} \, dS$ です。この量は、速度ベクトルの法線成分 $v_n = \bm{v} \cdot \bm{n}$ に微小面積をかけたもので、正ならば流体が外に出ていく（流出）、負ならば中に入ってくる（流入）ことを意味します。
体積流量に密度をかければ質量流量になるので、表面全体を通して流出する質量の総量は次のように表されます。
$$ \begin{equation} \dot{m}_{\text{out}} = \oint_S \rho \bm{v} \cdot \bm{n} \, dS \end{equation} $$
ここで $\oint$ は閉曲面全体にわたる積分を意味します。
質量保存の定式化
質量保存則は「検査体積内の質量の時間変化は、表面からの正味の流出量の符号を反転したものに等しい」と述べています。なぜ符号を反転するかというと、流出が正のとき内部の質量は減少するからです。すなわち、
$$ \frac{dM}{dt} = -\dot{m}_{\text{out}} $$
これを先ほどの式で書き換えると、
$$ \int_V \frac{\partial \rho}{\partial t} \, dV = -\oint_S \rho \bm{v} \cdot \bm{n} \, dS $$
右辺を左辺に移項して整理すると、連続の方程式の積分形が得られます。
$$ \begin{equation} \frac{\partial}{\partial t}\int_V \rho \, dV + \oint_S \rho \bm{v} \cdot \bm{n} \, dS = 0 \end{equation} $$
この式の物理的な意味を言葉で述べると、「検査体積内の質量の増加率と、表面からの質量の正味の流出率を合計するとゼロになる」ということです。先ほどのバスタブの例に立ち返ると、バスタブ内のお湯の増加量は、蛇口から入る量から排水口から出る量を引いたものに等しい、という直感と完全に一致しています。
積分形は、有限の大きさを持つ領域全体での質量収支を表しています。しかし、流体の各点での局所的な条件を知りたい場合には、微分形の方が便利です。次のセクションでは、積分形から微分形への変換を行います。
微分形の導出 — 発散定理による変換
積分形の連続の方程式を微分形に変換するには、ベクトル場の発散（div）の記事で扱ったガウスの発散定理を使います。
ガウスの発散定理の復習
ガウスの発散定理は、ベクトル場 $\bm{F}$ の閉曲面 $S$ を通る流束（面積分）を、体積 $V$ 内での発散の体積積分に変換する定理です。
$$ \oint_S \bm{F} \cdot \bm{n} \, dS = \int_V \nabla \cdot \bm{F} \, dV $$
直感的には、「表面全体から出ていく量の合計は、内部の各点で湧き出す量の合計に等しい」という意味です。ここで、$\nabla \cdot \bm{F}$ はベクトル場 $\bm{F}$ の発散で、各点における「湧き出しの強さ」を表すスカラー場です。
積分形への発散定理の適用
連続の方程式の積分形に現れる面積分に、$\bm{F} = \rho \bm{v}$ として発散定理を適用します。ここで $\rho \bm{v}$ は質量流束ベクトル（mass flux vector）と呼ばれ、単位面積・単位時間あたりに通過する質量を表すベクトルです。
$$ \oint_S \rho \bm{v} \cdot \bm{n} \, dS = \int_V \nabla \cdot (\rho \bm{v}) \, dV $$
この結果を、積分形の連続の方程式に代入します。
$$ \int_V \frac{\partial \rho}{\partial t} \, dV + \int_V \nabla \cdot (\rho \bm{v}) \, dV = 0 $$
左辺の2つの体積積分はどちらも同じ検査体積 $V$ 上の積分なので、1つにまとめることができます。
$$ \begin{equation} \int_V \left[\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v})\right] dV = 0 \end{equation} $$
被積分関数をゼロにする論法
ここで決定的な論法を使います。上の等式は任意の検査体積 $V$ に対して成り立たなければなりません。なぜなら、検査体積は私たちが自由に設定できる仮想的な領域であり、どこにどんな大きさ・形で設定しても質量保存は成立するからです。
もし被積分関数 $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v})$ がある点で正の値を取ったとすると、その点を含む十分小さな体積 $V$ を選べば、積分値は正になってしまい等式に矛盾します。同様に、負の値を取る点があっても矛盾します。
したがって、被積分関数はすべての点で恒等的にゼロでなければなりません。
$$ \begin{equation} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v}) = 0 \end{equation} $$
これが連続の方程式の微分形です。この式は流体の各点（局所的）で成り立つ条件であり、偏微分方程式の形をしています。
第1項 $\frac{\partial \rho}{\partial t}$ は「その場所での密度の時間変化率」を表し、第2項 $\nabla \cdot (\rho \bm{v})$ は「質量流束の発散 = その場所からの質量の正味の湧き出し率」を表しています。両者の和がゼロということは、密度が増加している場所には正味の質量流入があり、密度が減少している場所には正味の質量流出がある、ということです。
微分形は偏微分方程式なので、数値シミュレーションや解析的な問題解法に直接使うことができます。次に、この式を直交座標系で成分ごとに書き下してみましょう。
直交座標系での展開
微分形の連続の方程式を、直交座標系 $(x, y, z)$ において速度ベクトル $\bm{v} = (u, v, w)$ を使って成分ごとに書き下します。ここで $u$, $v$, $w$ はそれぞれ $x$, $y$, $z$ 方向の速度成分です。
発散 $\nabla \cdot (\rho \bm{v})$ を展開すると、各方向の質量流束の偏微分の和になります。
$$ \begin{equation} \frac{\partial \rho}{\partial t} + \frac{\partial(\rho u)}{\partial x} + \frac{\partial(\rho v)}{\partial y} + \frac{\partial(\rho w)}{\partial z} = 0 \end{equation} $$
この各項の意味を確認しましょう。
- $\frac{\partial \rho}{\partial t}$ : その点での密度の時間変化率
- $\frac{\partial(\rho u)}{\partial x}$ : $x$ 方向の質量流束の空間変化率
- $\frac{\partial(\rho v)}{\partial y}$ : $y$ 方向の質量流束の空間変化率
- $\frac{\partial(\rho w)}{\partial z}$ : $z$ 方向の質量流束の空間変化率
物質微分を使った表現
上の式の空間微分項を、積の微分の公式を使って展開してみましょう。たとえば $x$ 成分は $\frac{\partial(\rho u)}{\partial x} = u \frac{\partial \rho}{\partial x} + \rho \frac{\partial u}{\partial x}$ となります。他の成分も同様に展開すると、
$$ \frac{\partial \rho}{\partial t} + u\frac{\partial \rho}{\partial x} + v\frac{\partial \rho}{\partial y} + w\frac{\partial \rho}{\partial z} + \rho\left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}\right) = 0 $$
ここで、左辺の前半4項に注目します。これは密度 $\rho$ の物質微分（material derivative）$\frac{D\rho}{Dt}$ そのものです。物質微分とは、流体要素に乗って移動しながら観測したときの時間変化率で、次のように定義されます。
$$ \frac{D\rho}{Dt} = \frac{\partial \rho}{\partial t} + u\frac{\partial \rho}{\partial x} + v\frac{\partial \rho}{\partial y} + w\frac{\partial \rho}{\partial z} = \frac{\partial \rho}{\partial t} + \bm{v} \cdot \nabla \rho $$
また、後半の括弧内は速度場の発散 $\nabla \cdot \bm{v}$ です。したがって、連続の方程式は物質微分を使って次のようにコンパクトに書けます。
$$ \begin{equation} \frac{D\rho}{Dt} + \rho \nabla \cdot \bm{v} = 0 \end{equation} $$
この形は、流体要素に乗った視点（ラグランジュ的視点）での質量保存を表しています。「流体要素を追いかけたときの密度の変化率は、速度場の発散（= 流体の膨張率）に比例する」と読むことができます。流体が膨張すれば（$\nabla \cdot \bm{v} > 0$）密度は減少し、収縮すれば（$\nabla \cdot \bm{v} < 0$）密度は増加する — これは直感にもよく合います。
この物質微分を使った表現は、非圧縮性流体の条件を導く際に特に便利です。次のセクションで見ていきましょう。
非圧縮性流体の場合 — $\nabla \cdot \bm{v} = 0$ の物理的意味
非圧縮性とは何か
日常で扱う液体（水、油など）の密度は、圧力が変わってもほとんど変化しません。このような流体を非圧縮性流体と呼びます。厳密には「流体要素を追いかけたときに密度が変わらない」、すなわち、
$$ \frac{D\rho}{Dt} = 0 $$
が成り立つことを意味します。なお、気体であっても流速がマッハ数 $0.3$ 以下（目安）であれば、密度変化は無視できるため非圧縮の仮定が使えます。レイノルズ数や流れの速度スケールによって、圧縮性の効果を考慮すべきかどうかが決まります。
非圧縮連続の方程式の導出
物質微分の形の連続の方程式に $D\rho/Dt = 0$ を代入すると、即座に次の結果が得られます。
$$ \begin{equation} \nabla \cdot \bm{v} = 0 \end{equation} $$
直交座標系で書き下すと、
$$ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0 $$
物理的な意味
$\nabla \cdot \bm{v} = 0$ は「速度場の発散がゼロ」を意味します。発散とは、ある点から流体がどれだけ「湧き出して」いるかを表す量です。これがゼロということは、流体のどの微小領域をとっても、流入量と流出量が常に等しい — つまり、流体要素が膨張も収縮もしない、ということを意味します。
この条件は、ポテンシャル流れの解析や、渦度と循環の議論において基本的な仮定として頻繁に登場します。また、数値流体力学（CFD）の計算では、非圧縮性流体のシミュレーションにおいて $\nabla \cdot \bm{v} = 0$ を満たすように圧力場を調整する「圧力ポアソン方程式」を解く手法が広く使われています。
2次元の場合、$\nabla \cdot \bm{v} = 0$ はさらに直感的です。$x$ 方向に流速が増加する場所（$\partial u / \partial x > 0$）では、連続の方程式を満たすために必ず $y$ 方向に流速が減少（$\partial v / \partial y < 0$）しなければなりません。あたかも、一方に押し出された流体の分だけ別の方向から引き込まれるかのように振る舞います。
非圧縮性の仮定は多くの実用的な場面で優れた近似を与えてくれます。しかし、最もシンプルで直感的な連続の方程式の適用例を見るには、さらに次元を下げて1次元の管路の流れを考えるのが最適です。
1次元定常流への適用 — $A_1 v_1 = A_2 v_2$
定常流の仮定
まず、定常流（steady flow）を考えます。定常流とは、流れの様子が時間的に変化しない流れのことで、$\partial \rho / \partial t = 0$ が成り立ちます。蛇口を開けて一定の水量を流し続けているような状態がこれに相当します。
管路での質量保存
1次元の管路（パイプ）を考え、断面1（断面積 $A_1$、流速 $v_1$、密度 $\rho_1$）と断面2（断面積 $A_2$、流速 $v_2$、密度 $\rho_2$）の間で質量保存を適用します。
定常流なので、断面1から流入する質量流量と断面2から流出する質量流量は等しくなければなりません。単位時間に断面を通過する質量は「密度 $\times$ 断面積 $\times$ 流速」ですから、
$$ \begin{equation} \rho_1 A_1 v_1 = \rho_2 A_2 v_2 \end{equation} $$
この式は積分形の連続の方程式を1次元管路に適用した結果に他なりません。
非圧縮性流体の場合
水のような非圧縮性流体では $\rho_1 = \rho_2$ なので、密度が消えて、
$$ \begin{equation} A_1 v_1 = A_2 v_2 \end{equation} $$
すなわち、断面積と流速の積（= 体積流量）が一定であることを意味します。管の断面積が小さくなれば流速は増加し、断面積が大きくなれば流速は減少します。
これこそが、冒頭で述べたホースの現象の数学的な説明です。ホースの先端を指で潰して断面積を小さくすると、同じ体積流量を維持するために流速が増加し、水が勢いよく飛び出すのです。
この $A_1 v_1 = A_2 v_2$ の関係は、ベルヌーイの定理と組み合わせると非常に強力です。ベルヌーイの定理は流速と圧力の関係を表しますが、流速の変化自体は連続の方程式から決まります。つまり、連続の方程式が流速を決め、ベルヌーイの定理が圧力を決める、という役割分担になっています。
具体的な数値を入れて、連続の方程式の使い方を確認してみましょう。
具体例 — 管路の絞りにおける流速変化
問題設定
直径 $D_1 = 10$ cm のパイプが、途中で直径 $D_2 = 5$ cm に絞られています。入口（断面1）での流速が $v_1 = 2$ m/s のとき、絞り部（断面2）での流速 $v_2$ はいくらになるでしょうか。流体は水（非圧縮性）とします。
計算
円管の断面積は $A = \frac{\pi D^2}{4}$ なので、非圧縮連続の方程式 $A_1 v_1 = A_2 v_2$ を使います。
まず、$v_2$ について解きます。
$$ v_2 = \frac{A_1}{A_2} v_1 $$
断面積の比を直径で表すと、
$$ \frac{A_1}{A_2} = \frac{\pi D_1^2 / 4}{\pi D_2^2 / 4} = \left(\frac{D_1}{D_2}\right)^2 $$
数値を代入します。
$$ \begin{align} v_2 &= \left(\frac{D_1}{D_2}\right)^2 v_1 \\ &= \left(\frac{10}{5}\right)^2 \times 2 \\ &= 4 \times 2 \\ &= 8 \text{ m/s} \end{align} $$
直径が半分（断面積が $1/4$）になったため、流速は4倍になりました。この結果は直感にもよく合います。同じ量の水が狭い場所を通るためには、それだけ速く流れなければなりません。
体積流量の確認
念のため、両断面での体積流量 $Q = Av$ が一致することを確認しましょう。
$$ \begin{align} Q_1 &= A_1 v_1 = \frac{\pi \times 0.10^2}{4} \times 2 = \frac{\pi \times 0.01}{4} \times 2 \approx 0.01571 \text{ m}^3\text{/s} \\ Q_2 &= A_2 v_2 = \frac{\pi \times 0.05^2}{4} \times 8 = \frac{\pi \times 0.0025}{4} \times 8 \approx 0.01571 \text{ m}^3\text{/s} \end{align} $$
確かに $Q_1 = Q_2$ となり、質量保存（非圧縮なので体積保存）が成り立っています。
ここまでの理論を、Pythonを使って可視化してみましょう。まずは管路の絞りにおける流速変化をグラフで確認し、続いて2次元の流れ場を可視化します。
Pythonで連続の方程式を可視化する
管路の絞りによる流速変化
最初に、管路の断面積が変化するときに連続の方程式に従って流速がどのように変わるかをPythonで計算・可視化します。管路の形状をガウス関数で滑らかに絞る形にモデル化し、各位置での流速を $A_1 v_1 = A(z) v(z)$ から求めます。
import numpy as np
import matplotlib.pyplot as plt
# --- 管路の形状を定義 ---
z = np.linspace(0, 2, 500) # 管路の長さ方向 [m]
R_inlet = 0.05 # 入口の管半径 [m]（直径10cm）
R_min = 0.025 # 絞り部の管半径 [m]（直径5cm）
# ガウス関数で滑らかに絞る
R = R_inlet - (R_inlet - R_min) * np.exp(-0.5 * ((z - 1.0) / 0.15) ** 2)
A = np.pi * R ** 2 # 断面積
# --- 連続の方程式で流速を計算 ---
v_inlet = 2.0 # 入口流速 [m/s]
A_inlet = np.pi * R_inlet ** 2
v = A_inlet * v_inlet / A # A1*v1 = A(z)*v(z)
# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
# (a) 管路の形状
axes[0].fill_between(z, R * 100, 8, alpha=0.2, color='gray', label='Wall')
axes[0].fill_between(z, -R * 100, -8, alpha=0.2, color='gray')
axes[0].plot(z, R * 100, 'k-', linewidth=2)
axes[0].plot(z, -R * 100, 'k-', linewidth=2)
axes[0].set_ylabel('r [cm]')
axes[0].set_title('(a) Pipe geometry')
axes[0].set_ylim(-8, 8)
axes[0].grid(True, alpha=0.3)
axes[0].legend()
# (b) 断面積の変化
axes[1].plot(z, A * 1e4, 'g-', linewidth=2)
axes[1].set_ylabel('A [cm²]')
axes[1].set_title('(b) Cross-sectional area')
axes[1].grid(True, alpha=0.3)
# (c) 流速の変化
axes[2].plot(z, v, 'b-', linewidth=2, label='v(z)')
axes[2].axhline(y=v_inlet, color='r', linestyle='--', alpha=0.5,
label=f'v_inlet = {v_inlet} m/s')
axes[2].axhline(y=v_inlet * (R_inlet / R_min) ** 2, color='orange',
linestyle='--', alpha=0.5,
label=f'v_max = {v_inlet * (R_inlet / R_min) ** 2:.1f} m/s')
axes[2].set_xlabel('z [m]')
axes[2].set_ylabel('v [m/s]')
axes[2].set_title('(c) Flow velocity (continuity equation)')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("pipe_constriction.png", dpi=150, bbox_inches='tight')
plt.show()
このグラフから、以下のことが読み取れます。
- 管路が絞られる位置（$z = 1.0$ m 付近）で流速が急激に増加している — 断面積が最も小さくなる場所で流速がピーク（8.0 m/s）に達しており、先ほど手計算で求めた結果と一致しています
- 断面積と流速は逆の関係にある — 断面積のグラフ (b) と流速のグラフ (c) を比較すると、ちょうど鏡像のように対称的な振る舞いをしています。これは $A \times v = \text{const}$ という連続の方程式の直接的な帰結です
- 絞りの前後で流速は入口の値に戻る — 管径が元に戻れば流速も元に戻ります。これは質量保存が常に成り立っていることの現れです
2次元非圧縮流れ場の可視化
次に、2次元の非圧縮性流れ場を可視化します。ここでは、一様流とソース（湧き出し）の重ね合わせによる流れを扱います。この流れはポテンシャル流れの基本的な例で、半無限物体まわりの流れのモデルとして知られています。
import numpy as np
import matplotlib.pyplot as plt
# --- 一様流 + ソースの重ね合わせ ---
N = 30
x = np.linspace(-3, 3, N)
y = np.linspace(-3, 3, N)
X, Y = np.meshgrid(x, y)
# 一様流の速度
U_inf = 1.0
# ソース（原点に湧き出し）の強さ
Q_source = 5.0
# 原点でのゼロ割を防ぐ
r = np.sqrt(X ** 2 + Y ** 2)
r = np.where(r < 0.1, 0.1, r)
# 速度成分の計算
U = U_inf + Q_source / (2 * np.pi) * X / r ** 2
V = Q_source / (2 * np.pi) * Y / r ** 2
# 速さ（色付けに使用）
speed = np.sqrt(U ** 2 + V ** 2)
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
# (a) 速度ベクトル場
axes[0].quiver(X, Y, U, V, speed, cmap='coolwarm', alpha=0.8)
axes[0].set_title('(a) Velocity vector field (uniform + source)', fontsize=12)
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_aspect('equal')
axes[0].grid(True, alpha=0.3)
axes[0].plot(0, 0, 'ko', markersize=8, label='Source')
axes[0].legend()
# (b) 流線
axes[1].streamplot(X, Y, U, V, color=speed, cmap='coolwarm',
density=1.8, linewidth=1.5, arrowsize=1.2)
axes[1].set_title('(b) Streamlines (uniform + source)', fontsize=12)
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_aspect('equal')
axes[1].grid(True, alpha=0.3)
axes[1].plot(0, 0, 'ko', markersize=8, label='Source')
axes[1].legend()
plt.tight_layout()
plt.savefig("2d_flow_field.png", dpi=150, bbox_inches='tight')
plt.show()
このプロットから、以下のことが読み取れます。
- ソース（原点）から放射状に流体が湧き出している — 左図のベクトル場を見ると、原点付近では全方向に向かう速度成分が存在します。これはソースの特徴で、$\nabla \cdot \bm{v} > 0$ となる点です（ソースのある点だけは質量の湧き出しが許容されている特異点です）
- 一様流との重ね合わせにより、非対称な流れパターンが生じる — 右図の流線を見ると、ソースの左側（上流側）では流線が分岐し、ソースを避けるように流れています。この分岐する流線は半無限物体の表面とみなすことができます
- 流線の間隔が狭い場所ほど流速が大きい — 流線が密集している領域は赤色で表示されており、高速であることが色からも確認できます。これは連続の方程式の2次元版が表す現象です
発散の数値検証 — 非圧縮条件の確認
理論で学んだ $\nabla \cdot \bm{v} = 0$ を、数値的に確かめてみましょう。一様流単体、渦流れ、そしてソースを含む流れのそれぞれについて発散を計算し、非圧縮条件が満たされているかどうかを可視化します。
import numpy as np
import matplotlib.pyplot as plt
# --- 高解像度のグリッドを作成 ---
N = 100
x = np.linspace(-3, 3, N)
y = np.linspace(-3, 3, N)
X, Y = np.meshgrid(x, y)
dx = x[1] - x[0]
dy = y[1] - y[0]
r = np.sqrt(X ** 2 + Y ** 2)
r = np.where(r < 0.2, 0.2, r)
# --- 3種類の流れ場を定義 ---
# (a) 一様流 + 自由渦（非圧縮: div v = 0）
Gamma = 5.0 # 循環
U_a = 1.0 - Gamma / (2 * np.pi) * (-Y) / r ** 2
V_a = Gamma / (2 * np.pi) * X / r ** 2
# (b) 二重渦 (非圧縮: div v = 0)
U_b = -Y / r ** 2
V_b = X / r ** 2
# (c) ソース流 (原点で div v ≠ 0)
Q_s = 5.0
U_c = Q_s / (2 * np.pi) * X / r ** 2
V_c = Q_s / (2 * np.pi) * Y / r ** 2
# --- 発散を中心差分で数値計算する関数 ---
def compute_divergence(U, V, dx, dy):
"""速度場の発散を中心差分で計算"""
dU_dx = np.gradient(U, dx, axis=1)
dV_dy = np.gradient(V, dy, axis=0)
return dU_dx + dV_dy
div_a = compute_divergence(U_a, V_a, dx, dy)
div_b = compute_divergence(U_b, V_b, dx, dy)
div_c = compute_divergence(U_c, V_c, dx, dy)
# --- 可視化 ---
fig, axes = plt.subplots(2, 3, figsize=(18, 11))
titles = [
'(a) Uniform + vortex',
'(b) Pure vortex',
'(c) Source'
]
U_list = [U_a, U_b, U_c]
V_list = [V_a, V_b, V_c]
div_list = [div_a, div_b, div_c]
for i in range(3):
speed = np.sqrt(U_list[i] ** 2 + V_list[i] ** 2)
# 上段: 流線
axes[0, i].streamplot(X, Y, U_list[i], V_list[i],
color=speed, cmap='coolwarm',
density=1.5, linewidth=1)
axes[0, i].set_title(f'{titles[i]} — Streamlines', fontsize=11)
axes[0, i].set_aspect('equal')
axes[0, i].grid(True, alpha=0.3)
axes[0, i].set_xlim(-3, 3)
axes[0, i].set_ylim(-3, 3)
# 下段: 発散のカラーマップ
vmax = max(abs(div_list[i]).max(), 0.1)
im = axes[1, i].pcolormesh(X, Y, div_list[i],
cmap='RdBu_r', vmin=-vmax, vmax=vmax,
shading='auto')
axes[1, i].set_title(f'{titles[i]} — div v', fontsize=11)
axes[1, i].set_aspect('equal')
axes[1, i].set_xlim(-3, 3)
axes[1, i].set_ylim(-3, 3)
plt.colorbar(im, ax=axes[1, i], shrink=0.8)
plt.tight_layout()
plt.savefig("divergence_check.png", dpi=150, bbox_inches='tight')
plt.show()
# --- 数値で確認 ---
# 原点付近を避けた領域で発散の平均値を表示
mask = r > 1.0 # 原点付近の特異点を除外
print(f"一様流+渦: div v の平均 = {np.mean(np.abs(div_a[mask])):.6f}")
print(f"純粋渦: div v の平均 = {np.mean(np.abs(div_b[mask])):.6f}")
print(f"ソース: div v の平均 = {np.mean(np.abs(div_c[mask])):.6f}")
このコードの実行結果から、以下のことが読み取れます。
- (a) 一様流 + 渦、(b) 純粋渦は、原点付近を除いて発散がほぼゼロ — 下段のカラーマップがほぼ白色（ゼロ）になっており、$\nabla \cdot \bm{v} = 0$ が数値的にも確認できます。わずかに値が出るのは数値微分の離散化誤差によるもので、グリッドを細かくすれば小さくなります
- (c) ソース流は原点付近で発散が大きな正の値を持つ — 赤色の領域が原点近くに現れており、質量が湧き出していることを示しています。原点から離れると発散はゼロに近づきます
- 非圧縮流れでは、流線パターンが複雑でも発散はゼロを保つ — 渦の流線は回転するような複雑なパターンですが、膨張・収縮がないため発散はゼロです。発散と回転（渦度）は独立した概念であることが、この可視化から確認できます
管路流速の計算ツール
最後に、実用的な計算ツールとして、さまざまな管径比に対する流速変化を一覧表示するコードを紹介します。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ ---
v_inlet = 2.0 # 入口流速 [m/s]
D_inlet = 0.10 # 入口直径 [m]
# 絞り比（D2/D1）を変えて流速を計算
ratios = np.linspace(0.2, 1.0, 50)
D_outlet = D_inlet * ratios
# 連続の方程式: A1*v1 = A2*v2 → v2 = (D1/D2)^2 * v1
v_outlet = v_inlet / ratios ** 2
# --- 可視化 ---
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(ratios, v_outlet, 'b-', linewidth=2.5, label='Outlet velocity $v_2$')
ax.axhline(y=v_inlet, color='r', linestyle='--', alpha=0.6,
label=f'Inlet velocity $v_1$ = {v_inlet} m/s')
# 具体的な値をマーク
specific_ratios = [0.25, 0.5, 0.75, 1.0]
for ratio in specific_ratios:
v_val = v_inlet / ratio ** 2
ax.plot(ratio, v_val, 'ko', markersize=8)
ax.annotate(f'D₂/D₁={ratio}\nv₂={v_val:.1f} m/s',
xy=(ratio, v_val),
xytext=(ratio + 0.05, v_val + 2),
fontsize=9, ha='left',
arrowprops=dict(arrowstyle='->', color='gray'))
ax.set_xlabel('Diameter ratio $D_2 / D_1$', fontsize=13)
ax.set_ylabel('Outlet velocity $v_2$ [m/s]', fontsize=13)
ax.set_title('Continuity equation: velocity vs. diameter ratio', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_xlim(0.15, 1.05)
ax.set_ylim(0, 40)
plt.tight_layout()
plt.savefig("velocity_vs_diameter_ratio.png", dpi=150, bbox_inches='tight')
plt.show()
# --- 数値テーブル ---
print("=" * 50)
print(f"{'D2/D1':>8} {'D2 [cm]':>10} {'v2 [m/s]':>10} {'v2/v1':>8}")
print("=" * 50)
for ratio in [1.0, 0.75, 0.5, 0.4, 0.3, 0.25, 0.2]:
v2 = v_inlet / ratio ** 2
print(f"{ratio:>8.2f} {D_inlet * ratio * 100:>10.1f} {v2:>10.1f} {v2 / v_inlet:>8.1f}")
このコードの実行結果から、以下のことがわかります。
- 流速は直径比の2乗に反比例して急激に増加する — グラフは双曲線的な形をしており、管径が少し小さくなるだけで流速は急激に上がります。たとえば $D_2/D_1 = 0.5$（直径半分）で流速は4倍、$D_2/D_1 = 0.25$（直径1/4）では流速は16倍になります
- 直径比が小さくなるにつれ、流速変化はますます敏感になる — グラフの傾きが左に行くほど急になっていることから、細い管ではわずかな管径変化が大きな流速変化を引き起こすことがわかります。これは工業的な配管設計において、絞り部の寸法精度が重要であることを示唆しています
- 直径比が1のとき（管径一定）、流速は変化しない — これは当然ですが、連続の方程式の帰結として自然に導かれます
圧縮性流体への拡張
ここまでは主に非圧縮性流体を中心に議論してきましたが、連続の方程式はもちろん圧縮性流体にもそのまま適用できます。圧縮性が重要になるのは、音速に近い（あるいは超える）高速の流れです。
圧縮性1次元定常流
圧縮性流体の1次元定常流では、
$$ \rho_1 A_1 v_1 = \rho_2 A_2 v_2 $$
が成り立ちます。非圧縮の場合と異なり、密度 $\rho$ も変化するため、断面積と流速だけでは流れの状態が決まりません。
たとえば超音速ノズル（ラバルノズル）では、断面積を一度絞ってから再び広げることで、亜音速の流れを超音速まで加速できます。これは一見すると「断面積が大きくなったのに流速が上がる」という直感に反する現象ですが、密度の急激な減少（膨張）が断面積の増加を上回るために生じます。連続の方程式 $\rho A v = \text{const}$ は、このような複雑な現象も正しく記述します。
微分形での圧縮性
微分形の連続の方程式
$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v}) = 0 $$
は、圧縮性・非圧縮性の区別なく一般に成り立つ形です。圧縮性流体の場合、この式に加えてエネルギー方程式と状態方程式が必要になり、密度・速度・圧力・温度を連立して解くことになります。これに運動量の保存則（ナビエ・ストークス方程式）を加えた方程式系が、圧縮性流体力学の基礎方程式です。
非圧縮性の仮定が使えるかどうかの判断基準は、マッハ数 $\text{Ma} = v / c$（流速と音速の比）です。$\text{Ma} < 0.3$ 程度であれば密度変化は約5%以下に収まるため、非圧縮近似が妥当とされています。
圧縮性の話題はより発展的ですが、連続の方程式が圧縮性・非圧縮性の両方の場合を統一的に記述できる強力な方程式であることを認識しておくことが大切です。
連続の方程式の他の座標系での表現
実際の問題では、流れの形状に合わせて直交座標系以外の座標系を使う方が便利なことがよくあります。ここでは、よく使われる円筒座標系と球座標系での表現を紹介します。
円筒座標系 $(r, \theta, z)$
パイプ内の流れや回転する流体を扱う際に便利な座標系です。速度成分を $(v_r, v_\theta, v_z)$ として、連続の方程式は次のようになります。
$$ \begin{equation} \frac{\partial \rho}{\partial t} + \frac{1}{r}\frac{\partial(r \rho v_r)}{\partial r} + \frac{1}{r}\frac{\partial(\rho v_\theta)}{\partial \theta} + \frac{\partial(\rho v_z)}{\partial z} = 0 \end{equation} $$
$r$ に関する微分に $1/r$ の因子と $r$ の積がかかっているのは、円筒座標系の体積要素が $r \, dr \, d\theta \, dz$ であることに起因します。直感的には、原点から離れるほど「円周」が長くなるため、動径方向の流束は $r$ を掛けた形で表す必要があるのです。
非圧縮性流体の場合は、
$$ \frac{1}{r}\frac{\partial(r v_r)}{\partial r} + \frac{1}{r}\frac{\partial v_\theta}{\partial \theta} + \frac{\partial v_z}{\partial z} = 0 $$
となります。
球座標系 $(r, \theta, \phi)$
地球大気の流れや球対称な流れを扱う際に用いられます。速度成分を $(v_r, v_\theta, v_\phi)$ として、
$$ \begin{equation} \frac{\partial \rho}{\partial t} + \frac{1}{r^2}\frac{\partial(r^2 \rho v_r)}{\partial r} + \frac{1}{r\sin\theta}\frac{\partial(\rho v_\theta \sin\theta)}{\partial \theta} + \frac{1}{r\sin\theta}\frac{\partial(\rho v_\phi)}{\partial \phi} = 0 \end{equation} $$
この式は一見複雑ですが、各項の構造は直交座標系の場合と同じです。各方向の質量流束の空間変化率を足し合わせてゼロとおく、という質量保存の原理がそのまま反映されています。$r^2$ や $\sin\theta$ の因子は、球座標系の体積要素 $r^2 \sin\theta \, dr \, d\theta \, d\phi$ から来ています。
座標系を変えても物理法則そのものは変わりません。変わるのは数式の見た目だけです。問題の対称性に合った座標系を選ぶことで、式が大幅に簡単になることがあります。たとえば球対称な流れでは $v_\theta = v_\phi = 0$ かつ $\theta, \phi$ 依存性がないため、連続の方程式は $\frac{1}{r^2}\frac{\partial(r^2 \rho v_r)}{\partial r} = 0$ と単純化されます。
連続の方程式と他の基本方程式との関係
連続の方程式は、流体力学における支配方程式の一つですが、単独では流れの完全な解を得ることはできません。ここでは、連続の方程式が他の基本方程式とどのように関係しているかを整理します。
ナビエ・ストークス方程式との関係
ナビエ・ストークス方程式は運動量の保存則を表す方程式です。連続の方程式（質量保存）とナビエ・ストークス方程式（運動量保存）を連立することで、初めて流れの速度場と圧力場を決定できます。
非圧縮性流体の場合、
- 連続の方程式: $\nabla \cdot \bm{v} = 0$ （拘束条件）
- ナビエ・ストークス方程式: $\rho \left(\frac{\partial \bm{v}}{\partial t} + \bm{v} \cdot \nabla \bm{v}\right) = -\nabla p + \mu \nabla^2 \bm{v} + \rho \bm{g}$
ここで、連続の方程式は速度場が満たすべき拘束条件として機能します。ナビエ・ストークス方程式を解いた速度場が連続の方程式を自動的に満たすように、圧力場が調整されます。
ベルヌーイの定理との関係
ベルヌーイの定理は、定常・非粘性・非圧縮の流れにおいて流線に沿ったエネルギー保存を表します。
$$ p + \frac{1}{2}\rho v^2 + \rho g h = \text{const} $$
この定理を使って管路の圧力変化を求めるとき、流速 $v$ の値は連続の方程式から決まります。つまり、連続の方程式はベルヌーイの定理に流速の情報を供給する役割を担っています。管路の絞りで流速が増加すると、ベルヌーイの定理から圧力が低下する — この因果関係を正しく理解するには、両方の方程式が必要です。
渦度方程式との関係
非圧縮性流体の連続の方程式 $\nabla \cdot \bm{v} = 0$ は、渦度と循環の解析でも重要な役割を果たします。ナビエ・ストークス方程式の回転（curl）を取って渦度方程式を導く際、$\nabla \cdot \bm{v} = 0$ の条件が式の簡略化に使われます。
このように、連続の方程式は流体力学のあらゆる場面で他の方程式と連携して働く、まさに「基盤」と呼べる方程式です。
まとめ
本記事では、連続の方程式（質量保存則）について、直感的な理解から数学的な導出、そしてPythonによる可視化まで一貫して解説しました。
- 質量保存の本質: 流体は途中で消えたり湧き出したりしない — この直感を数式にしたものが連続の方程式です
- 積分形: $\frac{\partial}{\partial t}\int_V \rho \, dV + \oint_S \rho \bm{v} \cdot \bm{n} \, dS = 0$ — 検査体積全体での質量収支を表します
- 微分形: $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v}) = 0$ — 各点での局所的な質量保存を表します。ガウスの発散定理を使って積分形から導出しました
- 物質微分の形: $\frac{D\rho}{Dt} + \rho \nabla \cdot \bm{v} = 0$ — 流体要素に乗った視点での表現です
- 非圧縮性流体: $\nabla \cdot \bm{v} = 0$（速度場の発散がゼロ）— 流体要素が膨張も収縮もしないことを意味します
- 1次元管路: $A_1 v_1 = A_2 v_2$（断面積と流速の積が一定）— ホースの現象を定量的に説明します
- Python可視化: 管路の流速変化、2次元流れ場の流線、発散の数値検証を通じて理論を確認しました
連続の方程式は、流体力学の最も基本的な方程式でありながら、配管設計から気象予報、航空工学まで幅広い応用を持つ強力なツールです。この方程式を出発点として、運動量保存（ナビエ・ストークス方程式）やエネルギー保存（ベルヌーイの定理）と組み合わせることで、さまざまな流体現象を解析できるようになります。
次のステップとして、以下の記事も参考にしてください。
- ベルヌーイの定理 — 連続の方程式と合わせて使うことで、管路の圧力変化を予測できます
- ナビエ・ストークス方程式 — 質量保存に加えて運動量保存を考慮した、より一般的な流体の支配方程式です
- ベクトル場の発散（div） — 連続の方程式の核心である発散の概念をさらに深く学べます
- ポテンシャル流れ — 非圧縮・非粘性の仮定のもとでの流れの解析手法を扱っています