ナビエストークス方程式をひたすら作る

このスクラップでは、物理の概念と計算の流れに重点を置いて説明をしています。そのため、説明に文量が必要な数学的内容は別ページになっています。（もしくはないです）

参考図書は

「入門連続体の力学」　半揚 捻雄　日本評論社　と

「連続体力学（物理学基礎シリーズ）」　佐野 理　朝倉書店　です。

できるだけ正しく書いているつもりですが、間違っているところもあるかもしれません。この記事はイメージをつかむためにして、上記のような教科書でしっかり勉強することをお勧めします。（特に調べ学習の時）

０．どのような話の流れで作るのか、概略を説明します。

まず、このスクラップでは連続体力学において重要な式のうち２つを作ります。

違っている点としては、普通の運動方程式は一つの質点を考えていましたが、連続体の運動方程式は大きさのある変形する物体を考えます。

それぞれの物体はどのように変形するのか、物体特有の性質を表すには次の構成方程式が必要です。

のように書かれます。(現時点では式の意味を理解しなくて大丈夫です。)

そして、取り扱いたい物体に合わせて構成方程式を選び、まっさらな運動方程式に適用することで特定の物体の運動方程式が得られます。

特に今回の場合はニュートン流体の構成方程式を使うことで、ニュートン流体の運動方程式であるナビエストークス方程式が得られます。

1．まずは取り扱うものとはどんなものか考える。

まずは流体が空間に余すところなく流れている状態を考えます。ここから、流体の運動を特徴づける物理量を取り出します。

動力学を考えるなると、運動方程式の考えから一番に考えられるものは、位置だと思います。しかし、流体は考えたい空間の外から中に、中から外にとめどなく流れています。これでは、考えている流体部分が流れて行ってしまったら位置を把握することができません。

では、何をもって流体の運動を考えるかですが、これはそれぞれの時間、それぞれの位置の流れの速度を考えます。数学的に言うと、時間と三次元位置変数のベクトル場を考えます。

https://gyazo.com/8730514eb75eb6c2b0b1cfb75f558a6f

図１：流れ場

２．流体の運動方程式を作る

２－１．左辺、質量と加速度

https://gyazo.com/9caeb5c075dc44a205c57b75bf6880e2

図２：流れ場中の微小体積

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

２－２．右辺、力

続いて、物体にかかる力について考えます。

２－２－１．外力

２－２－２．応力

この章からは力(force)と応力(stress)を明確に区別して話を進めます。

なぜ応力は方向だけでなく、面の情報を持っているかというと、方向だけを持っている力では物体の変形が表現できないからです。

そして、面の方向の情報を持った応力は、物体の変形と適切な関係を持った力として使われます。

図３，４：xの面のx軸方向変形：xの面のy軸方向変形

応力を式に導入するためには一般的な形で応力を記述する必要があるので、以下で応力の記法を導き出します。

これを考えるときには、図のような四面体を考えます。

https://gyazo.com/bf73be7e41c08e0bd5398327bd26fe1d

図５：体積力と応力が釣り合っている四面体

$ \rho dV\mathbf{K}+\mathbf{p}_{\mathbf{n}}dS+ \mathbf{p}_{-x}dS_x + \mathbf{p}_{-y}dS_y + \mathbf{p}_{-z}dS_z=\mathbf{0}

として力の釣り合いの式が書けます。

私達は応力について考えたいのですが、体積力が邪魔です。そこでまずは体積力を消すことにします。

となり、邪魔な体積力を消すことができました。

邪魔な体積力がなくなったので、それぞれの面の面積と応力について具体的に議論を進めます。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

こうすると、

$ \mathbf{p}_{\mathbf{n}}dS+ \mathbf{p}_{-x}dS_x + \mathbf{p}_{-y}dS_y + \mathbf{p}_{-z}dS_z=\mathbf{p}_{\mathbf{n}}dS+\mathbf{p}_{-x}ldS+\mathbf{p}_{-x}mdS+\mathbf{p}_{-z}ndS=\mathbf{0}

としてより具体的に書けます。

https://gyazo.com/a9c9bc6be7529211dc5bce172120c71b

図６：応力を内向きにする

すると、式は

$ \mathbf{p}_{\mathbf{n}}=\begin{bmatrix} p_{xx}&p_{xy}&p_{xz} \\ p_{yx}&p_{yy}&p_{yz} \\ p_{zx}&p_{zy}&p_{zz} \end{bmatrix} \begin{bmatrix} l \\ m \\ n \end {bmatrix}

とできます。

$ \mathbf{p}_{\mathbf{n}}=P\mathbf{n}

として、法線ベクトルとの積で表すことができます。この法線ベクトルとの積での表現はとても便利で、連続体の運動方程式を微分形式で書く時に効果的です。

さて、応力の表現法がわかったのですが、実際に流体の部分に応力が働く時はどうなるのでしょうか。

https://gyazo.com/4efceff258c4322d1ecc56801ba9ce22

図７：体積Vの流体領域の状態

２－３．連続体の運動方程式の形にする

そして、質量×加速度＝力として運動方程式を書きます。

$ \rho \frac{D \mathbf{v}}{D t}dV= \rho \mathbf{K}dV+ \mathbf{p}_{\mathbf{n}}dS

まずは体積についての積分を考えます。

続いて面積についての積分は、面は裏表があるので、法線ベクトルを考える必要がありますが、ほとんど同じです。

さて、体積、面積について積分することで、微分形式の運動方程式は

$ \iiint_V \rho \frac{D \mathbf{v}}{D t}dV=\iiint_V \rho \mathbf{K}dV+\iint_S \mathbf{p}_{\mathbf{n}}d\mathbf{S}

となります。

ここで、右辺第2項については

$ \iint_S P\mathbf{n}d\mathbf{S}

と書き換えられて、ガウスの発散定理から

$ \iint_SP\mathbf{n}d\mathbf{S}=\iiint_V divPdV

として書き換えられて、体積積分で式を統一的に表現できます。

そして、式全体を体積について微分すると、

$ \rho \frac{D\mathbf{v}}{Dt}=\rho\mathbf{K}+divP

となり、これが連続体の運動方程式(微分形)になります。

３．構成方程式

以上の運動方程式をニュートン流体のものにするために、ニュートン流体における応力と変形の関係を特徴付ける構成方程式を導きます。

３－１－１．弾性体の変形テンソル・変形の記法

構成方程式を作るには、応力と変形の情報が必要です。応力は上記からよくわかっているので、物体の変形はどのようにして記述されるのかを決めなければなりません。以下ではまず弾性体である金属で変形の表現方法を確立し、それを流体に応用します。

まず簡単に、金属の棒が伸縮する場合を考えます。

$ p=E\frac{\Delta l}{l}

と書かれます。

https://gyazo.com/45289ee68cfe700c79a4267c6d0f1a44

図８：金属棒の伸び変形

(応力は外力でないという話はありますが、簡単のため省略します)

このひずみをさらに一般的に表現します。

https://gyazo.com/816b2ed81f513d284efa5d4d895b5d0b

図９：金属棒中の近接した点の変形

では、同じような話を3次元でしていきます。

$ \delta\mathbf{u}(x,y,z)=\mathbf{u}(x+dx,y+dy,z+dz)-\mathbf{u}(x,y,z)

と書きます。

$ \mathbf{u}(x+dx,y+dy,z+dz)=\mathbf{u}(x,y,z)+\textbf{a}

変形量の3次元での表現はわかったので、それを使って、1次元のときと同じようにひずみを微分で表現します。

$ \delta\mathbf{u}=\begin{bmatrix} \delta u_x \\ \delta u_y \\ \delta u_z \end{bmatrix}

$ = \begin{bmatrix} u_x(x+dx,y+dy,z+dz)-u_x(x,y,z)\\ u_y(x+dx,y+dy,z+dz)-u_y(x,y,z)\\ u_z(x+dx,y+dy,z+dz)-u_z(x,y,z)\end{bmatrix}

$ =\begin{bmatrix} \frac{\partial u_x}{\partial x} & \frac{\partial u_x}{\partial y} & \frac{\partial u_x}{\partial z} \\ \frac{\partial u_y}{\partial x} & \frac{\partial u_y}{\partial y} & \frac{\partial u_y}{\partial z} \\ \frac{\partial u_z}{\partial x} & \frac{\partial u_z}{\partial y} & \frac{\partial u_z}{\partial z} \end{bmatrix} \begin{bmatrix} dx \\ dy \\ dz \end{bmatrix}

変形テンソルの成分は1次元のときと同じように、ひずみを表しています。

さらに、3次元になったことで物体のどの面がどの方向にひずんでいるかまで考えられています。

この変形テンソルの中身を詳細に吟味していきます。

対角成分は変形量の方向と、基準の長さの方向が同じです。

https://gyazo.com/6574d87f193d362711e090bc1e852ce6

図１０：対角成分の伸縮ひずみ

そして、全体的には微小面がずれるようにひずみます。形は平行四辺形になります。

同じような話が他の対角でない成分にもできて、対角でない成分は物体のせん断ひずみを表現しているということになります。

https://gyazo.com/fcb5a512b57354ecb854491cdb55f208

図１１：対角でない成分のせん断ひずみ

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

３－１－２．弾性体のひずみテンソル

$ E= \begin{bmatrix} \frac{\partial u_x}{\partial x} &\frac{1}{2}(\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}) & \frac{1}{2}(\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x}) \\ \frac{1}{2}(\frac{\partial u_y}{\partial x}+\frac{\partial u_x}{\partial y}) & \frac{\partial u_y}{\partial y} & \frac{1}{2}(\frac{\partial u_y}{\partial z}+\frac{\partial u_z}{\partial y}) \\ \frac{1}{2}(\frac{\partial u_z}{\partial x}+\frac{\partial u_x}{\partial z}) & \frac{1}{2}(\frac{\partial u_z}{\partial y}+\frac{\partial u_y}{\partial z}) & \frac{\partial u_z}{\partial z} \end {bmatrix}

$ E = \begin{bmatrix} e_{xx} & e_{xy} & e_{xz} \\ e_{yx} & e_{yy} & e_{yz} \\ e_{zx} & e_{zy} & e_{zz} \end{bmatrix}

変形テンソルと同様に、ひずみテンソルも成分について詳細に見ていきます。

変形テンソルでは対角でない成分は、1つにつき1種のせん断ひずみを表していました。

https://gyazo.com/92aa970276a8c238da181c00606b46fa

図１２：ひずみテンソルによるせん断ひずみ

ここまで読むと「なんで変形に関する同じようなテンソルを2個も作ったんだ」と思われると思います。

それは、変形テンソルよりもひずみテンソルの方がよりよくひずみを表現できるからです。

図１３,１４：変形テンソルとひずみテンソルのせん断ひずみの違い

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

ちなみに、変形テンソルの回転の要素は変形テンソルから作られる反変テンソル

$ \Omega =\frac{1}{2}(D-D^T)

として出てきます。回転テンソルと言います。気象科学で大気にかかるコリオリ力を考えるときに使います。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

３－２－３．流体のひずみテンソル

ひずみ速度テンソルの成分を書き下すと以下のようになります。

$ \dot{E}= \begin{bmatrix} \frac{\partial u}{\partial x} &\frac{1}{2}(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}) & \frac{1}{2}(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}) \\ \frac{1}{2}(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}) & \frac{\partial v}{\partial y} & \frac{1}{2}(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}) \\ \frac{1}{2}(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}) & \frac{1}{2}(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}) & \frac{\partial w}{\partial z} \end {bmatrix}

さて、ようやく変形の記述方法を確立するという大仕事(主に字数的に)が終わりました。ここで一度応力とひずみについてまとめておきます。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

また、上三角の成分と下三角の成分が等しく対応するので、ひずみの種類は実質6種類になります。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

ここの変形やひずみの話は材料力学の教科書にも載っていますので、持っていたら確認してみると理解が深まると思います。

３－３－１．構成方程式の具体的な形

まず実は、先ほど出てきた金属の応力とひずみの関係式は金属の構成方程式の一部です。

これを原子的な視点から、具体的・物理的に導出してみます。

またこの議論は、このあとの抽象的な構成方程式の議論で出てくるマクローリン展開の物理的イメージを与えてくれます。

まず、力に垂直な物体の変形を考えたいとき、その物体の変形を相当にミクロな視点で考えると、その変形は物体を構成する原子の間の距離の変動にの総和と見れます。

まずは、一つの原子間の距離の変動を考えてみます。

図１５：原子のポテンシャル

ここで、このポテンシャルのグラフを、大胆に近似します。

具体的には、相手の原子が居るポテンシャルの谷を頂点とした下に凸の二次関数でポテンシャルを近似します。

そうして、近似がだめにならない範囲で原子間の距離が変動するとしたら、変動に必要な力は近似二次関数の微分で、

以上の情報から、ひずみと力の関係式を立てると

$ \frac{F}{\sigma A} = f = \frac{NE'\delta}{Nd}

$ p = E\frac{\Delta l}{l}

となって、ミクロな世界のポテンシャルと原子間距離の変動から、物理的な議論で構成方程式の一端が得られました。

続いて、一般的に、抽象的に構成方程式を考えていきます。

$ \begin{matrix} p_{xx}=f_{xx}(e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \\ p_{xy}=f_{xy}(e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \\ p_{xz}=f_{xz}(e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \\ p_{yx}=f_{yx} (e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \\ p_{yy}=f_{yy}(e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \\ p_{yz}=f_{yz}(e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \\ p_{zx}=f_{zx}(e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \\ p_{zy}=f_{zy}(e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \\ p_{zz}=f_{zz}(e_{xx},e_{xy},e_{xz},e_{yx},e_{yy},e_{yz},e_{zx},e_{zy},e_{zz}) \end{matrix}

これ全部をまとめて構成方程式です。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

となり、応力とひずみ速度の線形関係が導かれます。

とも表現できます。

では、ここからはこの構成方程式をニュートン流体のものにしていきます。

３－３－２．ニュートン流体の構成方程式

ということで、添字が同じ文字の時に1となり、異なるとき0と定義した

と表現できます。

また、流体はどの方向から力を加えても、同じように変形します。(実験的にはそのはずです)

木のように力をかける方向によって変形の仕方が異なるのは異方性という性質です。

とします。これは、座標系の変換によって値の変わらない４階のテンソルを考えると得られます。

座標系によって値が変わらないというところが、等方性を表現しています。

それがこちらです

３－３－３．ニュートン流体の構成方程式の特徴

急にこうなりますと言われても意味不明なので、ニュートン流体の構成方程式それぞれについて説明します。

$ p_{xx}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{xx}

$ p_{yy}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{yy}

$ p_{zz}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{zz}

つまり流体がどの位置でも同じ速度のとき、もしくは静止流体のとき

$ p_{xx}=p_{yy}=p_{zz}=-p

しかし、流体が動いている場合はこの限りではありません。

ここでは、考えやすくするために、非圧縮性を仮定します。

連続の方程式（流体が変形しても質量が一定であることを表す）より、

垂直方向の構成方程式：

$ p_{xx}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{xx}

$ p_{yy}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{yy}

$ p_{zz}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{zz}

から

$ p_{xx}=-p+2 \mu e_{xx}

$ p_{yy}=-p+2 \mu e_{yy}

$ p_{zz}=-p + 2\mu e_{zz}

一方向にしか速度勾配がないとき、等方性からその他の方向の応力は等しいと考えられます。

そこで、式を変形して

$ p_{xx}=-p+2 \mu e_{xx}

$ p_{yy}=-p-2 \mu (e_{zz} + e_{xx})

$ p_{yy} = p_{zz}

$ p_{xx}=-p+2 \mu e_{xx}

$ p_{yy} + p_{zz} = -2p - 2\mu(e_{xx} + e_{yy} + e_{zz} +e_{xx}) = -2p - 2\mu e_{xx}

となり、「一方に速く流れている時にそれ以外の方向の圧力は小さくなる」というベルヌーイの法則のお気持ち的なところがわかります。

ーーーー

ニュートン流体は伸び変形すると伸び方向の圧力は大きくなり、それ以外の方向の圧力は小さくなる性質を持っていると構成方程式から見ることができます。

（粘土の塊を伸ばしてみると、力を掛けた方向に伸びて細くなっていくのと同じようなものだと思います。）

https://gyazo.com/6f458277a4bf335d76c140091d17da38

図１６：面に垂直な応力による変形パターン

ーーーー

続いて、構成方程式の残りの六個分についてです。

書き下します。

$ p_{xy}=2\mu e_{xy}

$ p_{yx}=2\mu e_{yx}

$ p_{yz}=2\mu e_{yz}

$ p_{zy}=2\mu e_{zy}

$ p_{zx}=2\mu e_{zx}

$ p_{xz}=2\mu e_{xz}

となります。

$ p_{xy}=p_{yx}=2\mu e_{xy}

$ p_{yz}=p_{zy}=2\mu e_{yz}

$ p_{zx}=p_{xz}=2\mu e_{zx}

ともできます。一種類のせん断ひずみについて二つの応力が必ず対応するということを表しています。式の意味としてはニュートンの粘性法則を表しています。

ひずみがあるとき、面に対して水平な応力を流体はかけられていることになります。特に顕著な例が翼面での粘性摩擦抵抗力です。

粘性流体は物体と触れているところは速度ゼロと考えます。

実際に飛行機が飛行しているときも翼面では空気の速度はゼロで、その面から少し離れたところでは対気速度に近い速さになっています。

ここで、翼面が流れる空気の面に水平な応力をかけることで、流れを引き止めていると考えると、翼面は空気にせん断応力をかけているということになります。翼面が力をかけているということは、その反作用を翼面は受けるわけで、これが摩擦抵抗力です。

https://gyazo.com/26aa83b858d6cdb0ee54af592ebae226

図１７：せん断ひずみと応力の関係図

４．ナビエストークス方程式

さて、連続体の運動方程式とニュートン流体の構成方程式を作れたので、ナビエストークス方程式を作ります。

まず材料を並べます。

$ \rho \frac{D\mathbf{v}}{Dt}=\rho\mathbf{K}+divP

$ p_{xx}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{xx}

$ p_{yy}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{yy}

$ p_{zz}=-p+\lambda(e_{xx}+e_{yy}+e_{zz})+2 \mu e_{zz}

$ p_{xy}=2\mu e_{xy}

$ p_{yx}=2\mu e_{yx}

$ p_{yz}=2\mu e_{yz}

$ p_{zy}=2\mu e_{zy}

$ p_{zx}=2\mu e_{zx}

$ p_{xz}=2\mu e_{xz}

以上です。代入してゴイゴイ計算したらおわりです。

$ divP = div\begin{bmatrix} p_{xx}&p_{xy}&p_{xz} \\ p_{yx}&p_{yy}&p_{yz} \\ p_{zx}&p_{zy}&p_{zz} \end{bmatrix}=\begin{bmatrix}\frac{\partial p_{xx}}{\partial x} + \frac{\partial p_{xy}}{\partial y} + \frac{\partial p_{xz}}{\partial z} \\ \frac{\partial p_{yx}}{\partial x} + \frac{\partial p_{yy}}{\partial y} + \frac{\partial p_{yz}}{\partial z} \\ \frac{\partial p_{zx}}{\partial x} + \frac{\partial p_{zy}}{\partial y} + \frac{\partial p_{zz}}{\partial z} \end{bmatrix}

$ \begin{bmatrix} p_{xx}&p_{xy}&p_{xz} \\ p_{yx}&p_{yy}&p_{yz} \\ p_{zx}&p_{zy}&p_{zz} \end{bmatrix}=\begin{bmatrix}-p & 0 & 0 \\ 0 & -p & 0 \\ 0 & 0 & -p \end{bmatrix}+ \lambda \begin{bmatrix}e_{xx}+e_{yy}+e_{zz} & 0 & 0 \\ 0 & e_{xx}+ e_{yy}+e_{zz} & 0 \\ 0 & 0 & e_{xx}+ e_{yy}+ e_{zz} \end{bmatrix}+ 2\mu \begin{bmatrix} e_{xx} & e_{xy} & e_{xz} \\ e_{yx} & e_{yy} & e_{yz} \\ e_{zx} & e_{zy} & e_{zz} \end{bmatrix}

$ div\begin{bmatrix} p_{xx}&p_{xy}&p_{xz} \\ p_{yx}&p_{yy}&p_{yz} \\ p_{zx}&p_{zy}&p_{zz} \end{bmatrix}=-\nabla p + \lambda \nabla(div\mathbf{v})+\mu \nabla(div \mathbf{v}) + \mu \nabla^2 \mathbf{v}

$ = -\nabla p + (\lambda + \mu)\nabla(div\mathbf{v})+\mu\nabla^2\mathbf{v}

なので式は

$ \rho \frac{D\mathbf{v}}{Dt}=\rho\mathbf{K}-\nabla p + (\lambda + \mu)\nabla(div\mathbf{v})+\mu\nabla^2\mathbf{v}

となり、ナビエストークス方程式ができました！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！

（おわり）