📂流体力学

圧縮性ナビエ・ストークス方程式の導出

定理

$$ \mathbf{u} = \mathbf{u} \left( t ; \mathbf{x} \right) = \left( u_{1} \left( t ; \mathbf{x} \right) , u_{2} \left( t ; \mathbf{x} \right) , u_{3} \left( t ; \mathbf{x} \right) \right) $$ 特に、3次元空間で時刻 $t$ と空間座標 $\mathbf{x} = \left( x_{1} , x_{2} , x_{3} \right)$ の速度場を上のような速度ベクトルで表すとする。それと同様に、$p : \mathbb{R}^{3} \to \mathbb{R}$ は各座標で加えられる圧力 $p = p \left( \mathbf{x} \right)$ を表す。$\mathbf{u}$ がニュートン流体の流速であるならば、次の支配方程式に従う。 $$ {\frac{ \partial \mathbf{u} }{ \partial t }} + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} = - \nabla w + \nu \nabla^{2} \mathbf{u} + \left( {\frac{ 1 }{ 3 }} \nu + \xi \right) \nabla \left( \nabla \cdot \mathbf{u} \right) + \mathbf{g} $$ ここで $\nabla \cdot$ は発散、$\nu = \mu / \rho$ は動粘性係数 $\mu$ と密度 $\rho$ の比、$\nabla w = \nabla p / \rho$ は熱力学的仕事^{thermodynamic work}、$\xi = \zeta / \rho$ は体積粘性 $\zeta$ と密度の比である体積動粘性係数^{bulk kinematic viscosity}、$\mathbf{g}$ は重力加速度である。

説明

第一印象では式の長さだけでもおそろしいが、オイラー方程式と非圧縮性ナビエ–ストークス方程式の導出過程に従ってきたならば、項が一つずつ増えているに過ぎず、意外と簡単な一般化を経ていることが分かる。 $$ {\frac{ \partial \mathbf{u} }{ \partial t }} + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} = - \nabla w + \mathbf{g} $$ オイラー方程式は上のように圧力と重力だけを考える。 $$ {\frac{ \partial \mathbf{u} }{ \partial t }} + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} = - \nabla w + \nu \nabla^{2} \mathbf{u} + \mathbf{g} $$ ナビエ–ストークス方程式には粘性項 $\nu \nabla^{2} \mathbf{u}$ が追加される。 $$ {\frac{ \partial \mathbf{u} }{ \partial t }} + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} = - \nabla w + \nu \nabla^{2} \mathbf{u} + \left( {\frac{ 1 }{ 3 }} \nu + \xi \right) \nabla \left( \nabla \cdot \mathbf{u} \right) + \mathbf{g} $$ 圧縮性ナビエ–ストークス方程式は非圧縮性と同等の条件である $\nabla \cdot \mathbf{u} = 0$ を仮定できないため、右辺の三番目の項が残るに過ぎない。

圧縮性流体といえば通常気体を想定するので、非圧縮性ナビエ–ストークス方程式は液体を扱い、圧縮性ナビエ–ストークス方程式は気体を扱うと考えて差し支えない。

導出

発散定理: 3次元ベクトル関数 $\mathbf{F}$ に対して次が成り立つ。 $$ \int_{\mathcal{V}} \nabla \cdot \mathbf{F} dV = \oint_{\mathcal{S}} \mathbf{F} \cdot d \mathbf{S} $$ ここで $\nabla \cdot \mathbf{F}$ は ダイバージェンス、$\int_{\mathcal{V}}$ は体積積分、$\oint_{\mathcal{S}}$ は閉曲面積分である。

基本的にナビエ–ストークス方程式の導出過程と同様にコーシー応力テンソルに発散定理を適用するつもりだ。突然 $\nabla \cdot \sigma$ が出てきて理解できないなら、ナビエ–ストークス方程式の導出過程でより丁寧に説明してあるので該当ポストを熟読することを勧める。

$$ \begin{align*} \mathbf{f} =& \oint_{\partial V} \sigma \cdot d S \\ =& \int_{V} \nabla \cdot \sigma d V \\ \implies {\frac{ d \mathbf{f} }{ d V }} =& \nabla \cdot \sigma \end{align*} $$

発散定理により、$\rho D \mathbf{u} / D t$ の右辺に $d \mathbf{f} / d V = \nabla \cdot \sigma$ を加えれば我々の求める支配方程式が得られる。

コーシー応力テンソルとラメパラメータ: 特に、3次元空間で時刻 $t$ と空間座標 $\mathbf{x} = \left( x_{1} , x_{2} , x_{3} \right)$ の速度場を上のような速度ベクトルで表すとする。この流体が粘性と圧縮性を持つニュートン流体であるとする。等方性を仮定するコーシー応力テンソル $\sigma$ は対称化された勾配 $\varepsilon$ に対して次のように表せる。 $$ \sigma = - p I + 2 \mu \varepsilon + \lambda \tr \left( \varepsilon \right) I $$ ここで対称化された勾配 $\varepsilon$ は次のように定義される。 $$ \varepsilon (\mathbf{u}) = {\frac{ 1 }{ 2 }} \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} \right) $$

粘性と圧縮性の両方が考慮される場合のコーシー応力テンソルは圧力を含めて三つの項から成る。これらを $\mathbf{u}$ に対して表すために $\tr \left( \varepsilon \right)$ から計算してみよう。 $$ \tr \left( \varepsilon \right) = {\frac{ \partial u_{1} }{ \partial x_{1} }} + {\frac{ \partial u_{2} }{ \partial x_{2} }} + {\frac{ \partial u_{3} }{ \partial x_{3} }} = \nabla \cdot \mathbf{u} $$ $\varepsilon$ は定義そのものを代入すればよいので、$\sigma \left( \mathbf{u} \right)$ のベクトル形は体積粘性 $\zeta = \lambda + {\frac{ 2 }{ 3 }} \mu$ に対して次のようになる。 $$ \begin{align*} & \sigma \left( \mathbf{u} \right) \\ =& - p I + \mu \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} \right) + \lambda \left( \nabla \cdot \mathbf{u} \right) I \\ =& - p I + \lambda \left( \nabla \cdot \mathbf{u} \right) I + {\frac{ 2 }{ 3 }} \mu \left( \nabla \cdot \mathbf{u} \right) I - {\frac{ 2 }{ 3 }} \mu \left( \nabla \cdot \mathbf{u} \right) I + \mu \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} \right) \\ =& - p I + \zeta \left( \nabla \cdot \mathbf{u} \right) I + \mu \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} - {\frac{ 2 }{ 3 }} \left( \nabla \cdot \mathbf{u} \right) I \right) \end{align*} $$

$\nabla \mathbf{u}$ のダイバージェンス: $$ \nabla \cdot \nabla \mathbf{u} = \nabla^{2} \mathbf{u} $$ $\left( \nabla \mathbf{u} \right)^{T}$ のダイバージェンス: $$ \nabla \cdot \left( \nabla \mathbf{u} \right)^{T} = \nabla \left( \nabla \cdot \mathbf{u} \right) $$

さて $\nabla \cdot \sigma$ を直接展開してみよう。 $$ \begin{align*} & \rho {\frac{ D \mathbf{u} }{ D t }} - \rho \mathbf{g} \\ =& \nabla \cdot \sigma \\ =& \nabla \cdot \left( - p I + \zeta \left( \nabla \cdot \mathbf{u} \right) I + \mu \left[ \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} - {\frac{ 2 }{ 3 }} \left( \nabla \cdot \mathbf{u} \right) I \right] \right) \\ =& - \nabla p + \zeta \nabla \left( \nabla \cdot \mathbf{u} \right) + \mu \nabla \cdot \left[ \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} - {\frac{ 2 }{ 3 }} \left( \nabla \cdot \mathbf{u} \right) I \right] \\ =& - \nabla p + \zeta \nabla \left( \nabla \cdot \mathbf{u} \right) + \mu \left[ \nabla^{2} \mathbf{u} + \nabla \left( \nabla \cdot \mathbf{u} \right) - {\frac{ 2 }{ 3 }} \nabla \left( \nabla \cdot \mathbf{u} \right) \right] \\ =& - \nabla p + \zeta \nabla \left( \nabla \cdot \mathbf{u} \right) + \mu \left[ \nabla^{2} \mathbf{u} + {\frac{ 1 }{ 3 }} \nabla \left( \nabla \cdot \mathbf{u} \right) \right] \\ =& - \nabla p + \mu \nabla^{2} \mathbf{u} + {\frac{ 1 }{ 3 }} \mu \nabla \left( \nabla \cdot \mathbf{u} \right) + \zeta \nabla \left( \nabla \cdot \mathbf{u} \right) \\ =& - \nabla p + \mu \nabla^{2} \mathbf{u} + \left( {\frac{ 1 }{ 3 }} \mu + \zeta \right) \nabla \left( \nabla \cdot \mathbf{u} \right) \end{align*} $$

最後に $\rho \mathbf{g}$ を移項して両辺を $\rho$ で割れば、我々が求めていた方程式が得られる。

■