よくわかるナビエ・ストークス方程式ーver. 2。

ここ3週間くらい研究が進んでいなかったんですが，ついさっき良い進捗（半群の exponential decay estimate）を得たので，暇つぶしに記事を更新しようと思います．先輩にアウトリーチ活動は大切と言われたので，しっかり実行しようと思います．さすがに土日に論文を書きたくないのです...

さて，弊ブログの記事のアクセスランキングを見ると第2位がナビエ・ストークス方程式になっています．閲覧いただきありがとうございます．
watanabeckeiich.hatenablog.com

以前サークルの後輩に「ナビエ・ストークス方程式について調べたらべっくさんのブログが出てきましたよ」と言われたことがある．この記事を書いた当初は，数学専攻の読者を想定して書いたが，どうやら工学系の人も見ているみたいだ．そこで，今日は，前よりも物理的な説明を加え，実際にナビエ・ストークス方程式を（簡単に）導出しようと思う．

【目次】

1. 応力
2. 平衡方程式
3. 体積モーメント
4. 主応力
5. 等方応力・偏差応力
6. 変形速度テンソル
7. ナビエ・ストークス方程式の導出
8. 参考書

1. 応力

まずは，応力について考えよう．固体（例えば氷など）に外力が作用すると，固体は変形し，その内部にない力が発生する．ここでは，固体の変形が小さいと仮定しよう．いま，内力について考察したいので，固体内部の図形の面積・方向などの変化は小さいとして無視しよう．すなわち，物体の点は移動しないとして内力を考察する．

固体への外力の作用は

① 圧力や支持力のように，固体の境界面を通じて固体の各部分に間接的に作用する力

② 重力のように境界面を通じてではなく，固体の各部分に直接的に作用する力

の2通りが考えられる．① の力を表面力 (surface force) といい，これは境界面の単位面積当たりの量として定義されるベクトル量である．一方，② の力を体積力 (volume force) といい，単位体積当たりの量として定義されるベクトル量である．

表面力を $\mathbf{T}(\mathbf{n}) dS$ としたとき， $\mathbf{T}(\mathbf{n})$ を応力テンソルという．ここで，境界面を $S$ ，単位法線ベクトルを $\mathbf{n}$ ，単位面積を $dS$ とした．応力テンソルは，境界面上の点 $P$ と単位法線ベクトル $\mathbf{n}$ を指定すると決定される．また，作用・反作用の法則より $\mathbf{T}(-\mathbf{n})= - \mathbf{T}(\mathbf{n})$ が成り立つ．

2. 平衡方程式

ここでは，固体が体積力 $\mathbf{f}$ の作用のもとで運動しているする．また，固体の密度を $\rho$ とし，固体内に任意の領域 $V$ をとる．領域 $V$ の境界面 $\partial V$ の単位法線ベクトル $\mathbf{n}$ を $V$ の外側に向けてとるものとする．

領域 $V$ が受ける力は

・表面力： $\displaystyle \int_{\partial V} \mathbf{T}(\mathbf{n})\, dS,$

・体積力： $\displaystyle \int_V \mathbf{f} \, dV$

の2つである．また，領域 $V$ の慣性力は $\displaystyle \int_V \rho \mathbf{a} \, dV$ である．ここで， $\mathbf{a}$ は固体の各部分の加速度を表す．したがって，領域 $V$ における運動方程式（運動量保存則）は

$\displaystyle \int_{\partial V} \mathbf{T}(\mathbf{n})\, dS + \int_V \mathbf{f} \, dV = \int_V \rho \mathbf{a} \, dV$

と書ける．ガウスの発散定理と $\mathbf{f}= f_i e_i$ , $\mathbf{a}=a_i e_i$ , $\mathbf{T}(\mathbf{n})=(\sigma_{ij} n_j)e_i$ より

$\nabla_j \sigma_{ij} + f_i - \rho a_i=0$

が成り立つ．これを平衡方程式 (equation of equilibrium) という．特に，固体が静止している場合は， $a_i=0$ より

$\nabla_j \sigma_{ij} + f_i=0$

が成り立つ．

3. 体積モーメント

電場内で誘電体が影響を受ける場合，あるいは磁場内で磁性体が影響を受ける場合などでは，固体の各部分がモーメント（= 物体を回転させる力）を受けることがある．このとき，単位体積当たりに働くモーメントはベクトル量である．このベクトル量を体積モーメントといい $\mathbf{M}$ で表す．

上で説明した平衡方程式と角運動方程式を用いると

$\varepsilon_{ijk} \sigma_{kj} + M_i = 0$

を得る．ここで， $M_i$ は体積モーメント， $\sigma_{kj}$ は応力テンソル， $\varepsilon_{ijk}$ は $\varepsilon$ テンソル，すなわち，

$\varepsilon_{ijk} = \begin{cases} 1 & \text{($(i,j,k)$ が偶順列),}\\ -1 & \text{($(i,j,k)$ が奇順列),}\\ 0 & \text{($i,\,j,\,k$ のうち等しいものがある)}\\ \end{cases}$

である．特に，体積モーメント $\mathbf{M}$ が 0 であれば，応力テンソルは対称テンソル，すなわち $\sigma_{ij} = \sigma_{ji}$ である．

4. 主応力

固体内の点を通る微小面分 $dS$ をとり， $\mathbf{n}$ を微小面分 $dS$ の単位法線ベクトルとする．この微小面分を通じて作用する応力ベクトル $\mathbf{T}(\mathbf{n})$ の $\mathbf{n}$ 方向への成分： $\widetilde{\mathbf{T}}(\mathbf{n}) = \mathbf{n} \cdot \mathbf{T} (\mathbf{n})$ を $\mathbf{n}$ 方向の垂直応力という．

応力テンソルを $\sigma_{ij}$ とすると $\mathbf{T}(\mathbf{n}) = \sigma_{ij} n_j e_i$ であるから， $\mathbf{n}$ の第 $i$ 成分を $n_i$ とすれば，

$\widetilde{\mathbf{T}}(\mathbf{n}) = \sigma_{ij} n_i n_j,$

$n_i n_j =1$

である．微小面分 $dS$ 上の点 $P$ を固定し，単位ベクトル $\mathbf{n}$ を変化させたときの $\widetilde{\mathbf{T}}(\mathbf{n})$ の極値を点 $P$ での主応力といい， $\widetilde{\mathbf{T}}(\mathbf{n})$ が主応力となるような $\mathbf{n}$ の向きを点 $P$ における主方向という．

主応力は応力テンソル $\sigma_{ij}$ （ $\sigma_{ij}$ は対称テンソル）の固有値であり，主方向は各々の主応力（固有値）に対する $\sigma_{ij}$ の固有ベクトルの向きを表す．よって，ある点 $P$ における応力テンソルの互いに垂直な3つの主方向を座標軸とする直交座標系に関して，応力テンソルの成分は

$\displaystyle (\sigma_{ij}) = \begin{pmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_3 \end{pmatrix}$

と書ける．ここで， $\sigma_1,\, \sigma_2,\, \sigma_3$ は主応力である．

5. 等方応力・偏差応力

応力テンソル $\sigma_{ij}$ のトレース $\sigma_{ii}$ は不変量である．主応力を $\sigma_1, \, \sigma_2, \, \sigma_3$ とすれば， $\sigma_{ii} =\sigma_1 +\sigma_2 +\sigma_3$ である．そこで，

$p = \displaystyle \frac{1}{3}\sigma_{ii} = \frac{1}{3} (\sigma_1 +\sigma_2 +\sigma_3)$

を等方テンソル（hydrostatic stress）という．さらに， $\delta_{ij}$ をクロネッカーのデルタとして， $S_{ij}$ を

$S_{ij} = - p \delta_{ij} + \sigma_{ij}$

と定義する．このとき， $S_{ij}$ を偏差応力（deriatoric stress）という．ここで，応力テンソル $\sigma_{ij}$ と偏差応力 $S_{ij}$ の主方向は一致している．

6. 変形速度テンソル

一般に，流体の応力 $\mathbf{T}(\mathbf{n})$ は流体の速度の空間的な変化 $\displaystyle \frac{\partial u_i}{\partial x_j}$ によって引き起こされると考えられている．すなわち， $\sigma_{ij}$ は $\displaystyle \frac{\partial u_i}{\partial x_j}$ の関数である．ここで，ベクトル値関数 $\mathbf{u}=(u_1,u_2,u_3)$ の微分 $\nabla \mathbf{u}$ を

$\nabla \mathbf{u} = \begin{pmatrix} \frac{\partial u_1}{\partial x} & \frac{\partial u_1}{\partial y} & \frac{\partial u_1}{\partial z} \\\frac{\partial u_2}{\partial x} & \frac{\partial u_2}{\partial y} & \frac{\partial u_2}{\partial z}\\ \frac{\partial u_3}{\partial x} & \frac{\partial u_3}{\partial y} & \frac{\partial u_3}{\partial z}\end{pmatrix}$

と定義する．このとき， $\nabla \mathbf{u}$ の反対称部分・対称部分をそれぞれ

$\displaystyle \mathbf{G}(\mathbf{u}) =\frac{1}{2}(\nabla \mathbf{u} - {}^\top\!(\nabla \mathbf{u})),$
$\displaystyle \mathbf{D}(\mathbf{u}) =\frac{1}{2}(\nabla \mathbf{u} + {}^\top\!(\nabla \mathbf{u}))$

とおくと， $\nabla \mathbf{u} = \mathbf{G}(\mathbf{u}) + \mathbf{D}(\mathbf{u})$ と書くことができる．ここで， $\mathbf{D}(\mathbf{u})$ を変形速度テンソルといい，対称行列である．応力 $\mathbf{T}(\mathbf{n})$ は変形速度テンソル $\mathbf{D}(\mathbf{u})$ からのみ引き起こされることが知られている．

7. ナビエ・ストークス方程式の導出

さて，以上の準備を経て，ようやくナビエ・ストークス方程式の導出を説明することができる．以下， $\rho$ を流体の密度， $\mathbf{u}$ を流体の速度場， $\mu$ をずり粘性係数， $\zeta$ を体積粘性係数とする．ただし， $\rho , \, \mu, \, \zeta > 0$ を仮定する．

ここでは導出しないが，流体の質量保存則に当たる連続の式が成り立つことが知られている：

$\displaystyle \frac{\partial}{\partial t} \rho + \mathrm{div}\, (\rho \mathbf{u}) =0 .$

また，流体はニュートン流体であると仮定する．このとき，応力テンソル $\sigma_{ij}$ は

$\sigma_{ij} = 2 \mu D_{ij}(\mathbf{u}) + \zeta (\mathrm{div}\,\mathbf{u}) \delta_{ij}$

と書けることが知られている．これは，応力 $\sigma_{ij}$ が体積保存の変形運動に起因する成分 $2 \mu D_{ij}(\mathbf{u})$ と体積変化に起因する等方成分 $\zeta (\mathrm{div}\,\mathbf{u}) \delta_{ij}$ の線形結合で表されることを示している．ここで， $(D_{ij}(\mathbf{u}))_{ij} = \mathbf{D}(\mathbf{u})$ は変形速度テンソルである．これより，運動量流束テンソル $P_{ij}$ を

$P_{ij} = \rho u_i u_j + p \delta_{ij} - \sigma_{ij} = \rho u_i u_j + S_{ij}$

と定義すると，次の運動量保存則を得る：

$\displaystyle \frac{\partial}{\partial t} (\rho u_i) + \frac{\partial}{\partial x_j} ( \rho u_i u_j + p \delta_{ij} - \sigma_{ij}) = \rho f_i$

なぜ，運動量流束テンソルを上のように定義したかについては後述の流体力学のテキストを参照のこと．連続の式を使って式を整理すると

$\displaystyle \rho \bigg( \frac{\partial u_i}{\partial t}+ u_j \frac{\partial u_i}{\partial x_j } \bigg)= -\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} \sigma_{ij} + \rho f_i$

を得る．

さて，ずり粘性係数 $\mu$ および体積粘性係数 $\zeta$ は一般には圧力 $p$ や温度 $\theta$ の関数であり，定数でない．しかし，多くの場合に変化が少なく，定数と見做せる場合が多い．係数 $\mu,\, \zeta$ が定数のときは，

$\begin{align*}\displaystyle \frac{\partial}{\partial x_j} \sigma_{ij} &= \mu \bigg(\frac{\partial^2}{\partial x_j^2} u_i + \frac{\partial}{\partial x_i} \frac{\partial u_j}{\partial x_j} -\frac{2}{3}\frac{\partial}{\partial x_i} (\mathrm{div}\,\mathbf{u}) \bigg)+ \zeta \frac{\partial}{\partial x_i}D\\ &= \mu \Delta u_i + \bigg(\frac{1}{3} \mu + \zeta \bigg)\frac{\partial}{\partial x_i}(\mathrm{div}\,\mathbf{u})\end{align*}$

と計算できる．これを運動量保存の式に代入すれば，

$\displaystyle \rho \bigg( \frac{\partial u_i}{\partial t}+ u_j \frac{\partial u_i}{\partial x_j } \bigg)= -\frac{\partial p}{\partial x_i} + \mu \Delta u_i + \bigg(\frac{1}{3} \mu + \zeta \bigg)\frac{\partial}{\partial x_i}(\mathrm{div}\,\mathbf{u}) + \rho f_i$

を得る．これをベクトル表記に直すと

$\displaystyle \rho \bigg( \frac{\partial \mathbf{u}}{\partial t}+ (\mathbf{u} \cdot\nabla) \mathbf{u} \bigg)= - \nabla p + \mu \Delta \mathbf{u} + \bigg(\frac{1}{3} \mu + \zeta \bigg)\nabla(\mathrm{div}\,\mathbf{u}) + \rho \mathbf{f}$

を得る．これと連続の式をあわせた

$\begin{align*} \frac{\partial}{\partial t} \rho + \mathrm{div}\, (\rho \mathbf{u}) &=0,\\ \displaystyle \rho \bigg( \frac{\partial \mathbf{u}}{\partial t}+ (\mathbf{u} \cdot\nabla) \mathbf{u} \bigg) &= - \nabla p + \mu \Delta \mathbf{u} + \bigg(\frac{1}{3} \mu + \zeta \bigg)\nabla(\mathrm{div}\,\mathbf{u}) + \rho \mathbf{f} \end{align*}$

をナビエ・ストークス方程式（Navier-Stokes equations）という．

特に，流体が非圧縮，一様密度の場合，すなわち， $\mathrm{div}\,\mathbf{u} =0$ のときは式がもう少し簡単になる．実際に，密度 $\rho$ は正定数 $\rho_0$ になるので，方程式は

$\begin{align*} \mathrm{div}\, \mathbf{u} &=0,\\ \displaystyle \frac{\partial \mathbf{u}}{\partial t}+ (\mathbf{u} \cdot\nabla) \mathbf{u} &= - \nabla \bigg( \frac{p}{\rho_0} \bigg) + \nu \Delta \mathbf{u} + \rho \mathbf{f} \end{align*}$

となる．通常これをナビエ・ストークス方程式ということが多い．ここで， $\displaystyle\nu = \frac{\mu}{\rho_0}$ は動粘性係数である．これに対してさっきの粘性係数 $\mu$ を力学的粘性係数という．これは， $\mu$ が粘性摩擦力の係数となるからである．

8. 参考書

私は流体力学の授業を受けたことがないので，独学で流体力学の知識を得た（つもりです）．そのときに参考になった参考書を挙げておく．

・神部勉著「流体力学」
知りたかったことがコンパクトにまとまっていて，大変参考になった．いまでもたまに参考にする．数学専攻の人はこの本で十分だと思う．

流体力学

作者: 神部勉,石井克哉
出版社/メーカー: 裳華房
発売日: 1995/09/01
メディア: 単行本
クリック: 4回
この商品を含むブログを見る

・巽友正著「流体力学」
これも丁寧でわかりやすかった．ただ，神部先生の本の方が薄かったので，この本はあまり読んでない．神部先生の本でわからなかったことを調べるのに使用した．

流体力学 (新物理学シリーズ)

作者: 巽友正
出版社/メーカー: 培風館
発売日: 1995/09/01
メディア: 単行本
購入: 1人クリック: 1回
この商品を含むブログ (1件) を見る

・今井功著「流体力学（前編）」
いろいろ詳しかった．等角写像を使った流体の解析が詳しく載っていた．ただ，等角写像を使った解析はどうやら時代遅れらしいので，あんまり一生懸命勉強する必要はないのかもしれない（あんまり勉強してないのでえらそうなことはいえませんが...）．前編があるなら後編もあるはずだと思って，一生懸命探したけども，後編を執筆中に今井先生が亡くなったので，後編は幻の本です．