A1.ナビエ・ストークス方程式

物体の運動方程式は\(ma=F\)であった。 流体の運動方程式は、ナビエ・ストークス方程式と呼ばれ、以下の形を持つ。

\[\frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \cdot \nabla)\boldsymbol{v} =　- \frac{1}{\rho}\nabla p + \frac{\mu}{\rho}\nabla^2\boldsymbol{v} +\boldsymbol{F}\]

ここでは、この方程式のかんたんな導出と、円筒座標系への変換について書いておく。

流体素片の運動方程式の導出

いま、とある流体素片があったとして、この流体素片の運動方程式を導出してみよう。 流体素片の体積を\(\delta V\)としておく。すなわち、この流体素片の質量は

\[(流体素片の質量)=\rho \delta V\]

となっている。この流体素片の加速度はラグランジュ微分を用いて\(\frac{d\boldsymbol{v}}{dt}\) と与えられそうなので、運動方程式ma=Fは、

\[\rho \delta V \frac{d\boldsymbol{v}}{dt} = \delta \boldsymbol{F}_{\rm body}+\delta \boldsymbol{F}_{\rm surface}\]

と書ける。ここで、右辺のちからを2つに分けた。右辺第一項は、流体素片の全体に掛かる力である（体積力）。例えば重力など。そのため、この力はかんたんに

\[\delta \boldsymbol{F}_{\rm body}=\rho \delta V \boldsymbol{F}\]

と書けるであろう。これは質点の運動方程式と変わらない。 問題になるのは第二項、表面力である（連続体の力学では面積力とか言うかも）。この力は、流体素片の表面にかかる力だとしよう。

ここで、流体素片の面積要素をベクトルdS(太字)とおくと、おそらく表面力は面積要素の大きさに比例すると思われるので、二階のテンソルを用いて

\[(\delta \boldsymbol{F}_{\rm surface})_{i}= -P_{ij} d\boldsymbol{S}_{ij}\]

とかけそうな気がしてくる。これを流体素片の表面全体で積分するとき、体積積分に直すときは一回微分しておけばよかったので

\[(\delta \boldsymbol{F}_{\rm surface})_{i}= -\int \frac{\partial P_{ij}}{\partial x_j} dV\]

これを元の式に戻すと流体素片の運動方程式は、i成分を取り出して書くと

\[\rho \frac{dv_{i}}{dt} = \rho F_{i} - \frac{\partial P_{ij}}{\partial x_j}\]

という形に書ける。

表面力の分解 - 圧力勾配力と粘性項

ここで、この右辺第二項を、圧力項と粘性項に分けて書いてみよう。クロネッカーのデルタを思い出して、対角成分だけ別に書くと

\[P_{ij}=p\delta_{ij}+\pi_{ij}\]

と分けることができる。

もし、第二項\(\pi_{ij}\)が0であれば、この場合は、表面力は全てその法線方向に働く。これを理想流体とよんでいる。 だが我々は、日常的な経験の中で、流体の抵抗力を知っている。例えば「ねばりけ」というのは、とある流体素片が動いているときに、そのそばにいるが、速度ベクトルの先にはいない別の流体素片に対してなにかしら力が作用している、ということである。

…ちょっと例を出そう。流れ場にシアがあった場合、\(\pi_{ij}\)はどんな形を持つだろうか。

(途中。P79あたり。5Viscous Flows)

圧力項

どこでもドアを月につないだらどうなるだろう？たぶん地球の空気が凄まじい勢いで月に噴出するであろう。これが圧力勾配力である。 何も難しいことは言っておらず、もし圧力に空間分布があり、勾配があれば、すなわち圧力pのグラジエントに比例して力がかかる効果である。 圧力は

\[p=\rho c_s^2\]

だったことを思い出すと、原始惑星系円盤では基本的に圧力は動径方向には内側ほど高く、鉛直方向には中心面ほど高い。 すなわち、圧力勾配力は外側に押し出すように働いている。

粘性項

(まだ書いてない)

円筒座標系への変換

まず、円筒座標系における発散等をおさらいしよう。 \(\psi\)をスカラー場、\(\mathbf{A}=A_r \hat{\mathbf{e}}_r+A_\theta\hat{\mathbf{e}}_\theta+A_z\hat{\mathbf{e}}_z\)をベクトル場とすると、下記が成り立つ。

\[ \begin{aligned} & \nabla \psi=\frac{\partial \psi}{\partial r} \hat{\mathbf{e}}_r+\frac{1}{r} \frac{\partial \psi}{\partial \theta} \hat{\mathbf{e}}_\theta+\frac{\partial \psi}{\partial z} \hat{\mathbf{e}}_z, \\ & \nabla \cdot \mathbf{A}=\frac{1}{r} \frac{\partial}{\partial r}\left(r A_r\right)+\frac{1}{r} \frac{\partial A_\theta}{\partial \theta}+\frac{\partial A_z}{\partial z}, \\ & \nabla \times \mathbf{A}=\left(\frac{1}{r} \frac{\partial A_z}{\partial \theta}-\frac{\partial A_\theta}{\partial z}\right) \hat{\mathbf{e}}_r+\left(\frac{\partial A_r}{\partial z}-\frac{\partial A_z}{\partial r}\right) \hat{\mathbf{e}}_\theta+\frac{1}{r}\left[\frac{\partial}{\partial r}\left(r A_\theta\right)-\frac{\partial A_r}{\partial \theta}\right] \hat{\mathbf{e}}_z, \\ & \nabla^2 \psi=\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial \psi}{\partial r}\right)+\frac{1}{r^2} \frac{\partial^2 \psi}{\partial \theta^2}+\frac{\partial^2 \psi}{\partial z^2} . \end{aligned} \]

(人生で一度くらいは導いてみてもいいかもしれない)。

Chouduri, Appendix C1, P393

これをナビエ・ストークス方程式に当てはめると、円筒座標系におけるナビエ・ストークス方程式が得られる。

\[ \begin{aligned} & \frac{\partial v_r}{\partial t}+v_r \frac{\partial v_r}{\partial r}+\frac{v_\theta}{r} \frac{\partial v_r}{\partial \theta}+v_z \frac{\partial v_r}{\partial z}-\frac{v_\theta^2}{r}\\ =-\frac{1}{\rho} \frac{\partial p}{\partial r} &+\nu\left(\frac{\partial^2 v_r}{\partial r^2}+\frac{1}{r^2} \frac{\partial^2 v_r}{\partial \theta^2}+\frac{\partial^2 v_r}{\partial z^2}+\frac{1}{r} \frac{\partial v_r}{\partial r}-\frac{2}{r^2} \frac{\partial v_\theta}{\partial \theta}-\frac{v_r}{r^2}\right)+F_r, \\ & \frac{\partial v_\theta}{\partial t}+v_r \frac{\partial v_\theta}{\partial r}+\frac{v_\theta}{r} \frac{\partial v_\theta}{\partial \theta}+v_z \frac{\partial v_\theta}{\partial z}+\frac{v_r v_\theta}{r}=-\frac{1}{\rho r} \frac{\partial p}{\partial \theta} \\ &+\nu\left(\frac{\partial^2 v_\theta}{\partial r^2}+\frac{1}{r^2} \frac{\partial^2 v_\theta}{\partial \theta^2}+\frac{\partial^2 v_\theta}{\partial z^2}+\frac{1}{r} \frac{\partial v_\theta}{\partial r}+\frac{2}{r^2} \frac{\partial v_r}{\partial \theta}-\frac{v_\theta}{r^2}\right)+F_\theta, \\ &\frac{\partial v_z}{\partial t} +v_r \frac{\partial v_z}{\partial r}+\frac{v_\theta}{r} \frac{\partial v_z}{\partial \theta}+v_z \frac{\partial v_z}{\partial z} \\ & =-\frac{1}{\rho} \frac{\partial p}{\partial z}+\nu\left(\frac{\partial^2 v_z}{\partial r^2}+\frac{1}{r^2} \frac{\partial^2 v_z}{\partial \theta^2}+\frac{\partial^2 v_z}{\partial z^2}+\frac{1}{r} \frac{\partial v_z}{\partial r}\right)+F_z . \end{aligned} \]

Chouduri, Appendix C2, P393

円筒座標系での連続の式

ついでだし、円筒座標系における連続の式もメモしておこう。 ベクトルの形で書いた連続の式

\[\frac{\partial \rho}{\partial t}+\nabla \cdot( \rho v)=0\]

に対して、先程の円筒座標系における\(\nabla\)をただ代入すると

\[\frac{\partial \rho}{\partial t}+\frac{1}{r}\frac{\partial}{\partial r}(r \rho v_r) + \frac{1}{r}\frac{\partial}{\partial \theta}(\rho v_\theta) + \frac{\partial}{\partial z} (\rho v_z)=0\]

となる