流体の運動量保存則(4) 面積力のモデリング
モデリング
引き続き、流体の運動量保存則とはどのようなものになるかを考えていきます。
前回の復習で、面積力を\(f_{s_i} \)、体積力を\(f_{v_i} \)として、運動量保存則は以下のように書けるのでした。
$$
\begin{align}
&\frac{ \partial \rho u_i \Delta x_1 \Delta x_2 \Delta x_3}{ \partial t} =- \frac{\partial \rho u_i u_j}{\partial x_j} \Delta x_1 \Delta x_2 \Delta x_3 +f_{v_i} + f_{s_i} \\\\
\tag{1} \label{eq:scale-factor-1}
\end{align}
$$
そして、これを、
$$
\begin{align}
&\tau_{ij} = \tau_{ji}
\tag{2} \label{eq:scale-factor-2}
\end{align}
$$
を使いながらtaylor展開で面積力と体積力を具体化すると、
$$
\begin{align}
&\frac{ \partial \rho u_i \Delta x_1 \Delta x_2 \Delta x_3}{ \partial t} =- \frac{\partial \rho u_i u_j}{\partial x_j} \Delta x_1 \Delta x_2 \Delta x_3 + \frac{\partial \tau_{ij}}{\partial x_j} \Delta x_1 \Delta x_2 \Delta x_3 – \rho g \delta_{i3} \Delta x_1 \Delta x_2 \Delta x_3 \\\\
&両辺\Delta x_1 \Delta x_2 \Delta x_3で割って、 \\\\
&\frac{ \partial \rho u_i }{ \partial t} =- \frac{\partial \rho u_i u_j}{\partial x_j} + \frac{\partial \tau_{ij}}{\partial x_j} – \rho g \delta_{i3}
\tag{3} \label{eq:scale-factor-3}
\end{align}
$$
この式こそが流体の運動量保存則です。
しかし、具体的ではないので、もう少し工夫する必要があります。
また、ここで、Taylor展開で高次元無視をしていることと、\(\tau_{ij} = \tau_{ji}\)をしていることには留意する必要があります。
ここで、大胆なことを考えます。
$$
\begin{align}
&\tau_{ij} = \tau_{ij}^{p} + \tau_{ij}^{v}
\tag{4} \label{eq:scale-factor-4}
\end{align}
$$
\(\tau_{ij}^{p} \) : 圧力起因 、 \( \tau_{ij}^{v}\) : 速度起因
なんてことをしてしまいます。つまり、面積力である、\(\tau_{ij}\)を圧力から受ける成分と速度から受ける成分に分解します。しかし、こんなことしていいのでしょうか？だって、速度変化によって圧力が生じたりするので、綺麗に分解など不可能なはずです。また、なんで足し算？
これは、モデリングです。そう、ここでモデリングが入ります。本当はそんな綺麗にはなってないだろうけど、こういう風に表してみましょうということです。
圧力起因
ここで、圧力起因の、\(\tau_{ij}^{p} \) は簡単です。圧力から受ける力は、圧力をPとすれば、Pですから。しかし、テンソル量なので、工夫が入ります。また、圧力は勾配の向きと力の向きが逆なので、マイナスがつきます。そうすると〜〜〜
$$
\begin{align}
&\tau_{ij}^{p} = – P \delta_{ij}
\tag{5} \label{eq:scale-factor-5}
\end{align}
$$
分かりましたか？Pをテンソル量に合わせるため、クロネッカーのデルタで処理しています。11,22,33の成分しか残らないので、法線応力だけが残り、うまく表せました。Pにマイナスがつくのは、この運動量保存則の第一成分を考えてみると、
$$
\begin{align}
&\frac{\partial \rho u_i}{\partial t} = \cdots – \frac{\partial P}{\partial x_1}
\tag{6} \label{eq:scale-factor-6}
\end{align}
$$
となっていれば、右側がx_1の正の向きとして、\(\frac{\partial P}{\partial x_1}\)が負の時、つまり圧力が左側が高くて右側が低い時に、右への速度は速くなるはずなので、マイナスをつけて\(\frac{\partial P}{\partial x_1}\)が正になれば、\(\frac{\partial \rho u_i}{\partial t}\)が大きくなり、つまり次の瞬間に速度が速くなります。
というより、基本的にスカラー値の微分と力の方向は逆になります。後々使いますが、例えばフーリエの法則とか。
これで、圧力の方はできました。
ニュートンの粘性法則
そして、速度起因を考えねばならないのですが、その前にニュートンの粘性法則を知っておく必要があります。なぜなら、これを使ってモデリングするからです。
下図のような地面と上の板に挟まれた水で、上の板が動くシチュエーションを考えてみます。教科書で良くみるやつですね。層流になるくらいの狭めのものと考えてください。
意外と奥深いのです。
まず、流体に粘性があると考えます。その時、上の板が一定の速度で動くことで、それに引きずられてすぐ下の流体は引きずられて速度が同じになろうとします。しかし、下は固体面との摩擦で止まろうとします。その結果、流体には上の図のような速度勾配ができます。
流体には下の固体の摩擦面から受けるせん断応力が流体に伝わり、どこでも同じせん断応力がかかります。そして、この時のせん断応力が以下のような式になります。
$$
\begin{align}
&\tau = \mu \frac{d u_i}{d x_3}
\tag{7} \label{eq:scale-factor-7}
\end{align}
$$
この時の、\(\mu\)を、分子粘性係数[Pa・s]または、分子動粘性係数[m2/s]と言います。これらは単に単位の違いで、分子動粘性係数は、分子粘性係数を密度で割ったものです。また、わざわざ分子とつけているのは、分子の粘性によるものだと示すためです。のちに乱流による同じような効果が出てくるので、それと区別する意図です。
そして、この力と速度の関係性をニュートンの粘性法則といい、この法則に従う流体をニュートン流体と言います。言葉にすると、流れのせん断応力（接線応力）と流れの速度勾配（ずり速度、せん断速度）が比例した粘性の性質を持つ流体のことです。本来は全方向に効いてくるのでテンソルで考えますが、ややこしいので置いときます。
あと、実はこの速度勾配は層流出ないとありえません。乱流が起きていたらこのような速度勾配になりません。
余談ですが、分子粘性係数\(\mu\)と分子動粘性係数\(\nu\)は面白くて、水と空気だと、
空気 : \(\mu = 1.822 \times 10^{-5}\)[Pa・s]、\(\nu = 1.512 \times 10^{-5}\)[m2/s]
水 : \(\mu = 1.004 \times 10^{-3}\)[Pa・s]、\(\nu = 1.004 \times 10^{-6}\)[m2/s]
という感じで、分子粘性係数は空気の方が小さいのですが、密度で割った分子動粘性係数は水の方が小さくなってます。動粘性係数は[m2/s]で、動きやすさを表し、つまりは水の方が動きにくいということです。でも、密度を考慮しなければ空気の方が動きにくいということです。
比熱で熱容量の関係に似ています。比熱が小さくても密度が大きいため、熱容量は大きいコンクリートのような。
層流じゃないとありえないとしましたが、実際はこんな感じになります。
同じように、この場合でもニュートンの粘性法則はありますが、乱流の影響もあるのでニュートンの粘性法則から出される摩擦と、実際の摩擦とは違います。粘性底層では同じです。
粘性底層は、本当にわずかなエリアで、例えば普通の外の建物くらいのスケールだと、数mmも行かないそんなレベルの狭い領域です。この範囲は層流のようになり、分子粘性が支配的です。
また、上図のような見方をすることもできます。上の黒い矢印のように、運動量が供給されているという見方です。まず、上の板が止まれば、この流体は全体的に止まります。そのように、上から運動量が供給されている、つまり流体から見た摩擦が運動量を供給しているものだという見方です。
実際、摩擦は応力なので、[N/m2]で、これは[kg・m/s / m・s]つまり、運動量フラックスになるのでした。
また、上図の右のように、せん断力の摩擦と、運動量の供給量が同じであり、\(\tau_{13} = \tau_{31}\)であり、対称テンソルとなっています。これは、対称テンソルになることの別の観点から見た説明であり、摩擦は運動量の供給と同じあるとも言えます。
速度起因
さて、\eqref{eq:scale-factor-4} 式で分けた、速度起因の部分をどう表現するべきか考えていきましょう。先ほど考えたモデルであれば、
$$
\begin{align}
&\tau_{13} = \mu \frac{\partial u_1}{\partial x_3}
\tag{8} \label{eq:scale-factor-8}
\end{align}
$$
こうなります。
全ての方向をまとめてかくと、
$$
\begin{align}
&\tau_{ij} = \mu \frac{\partial u_i}{\partial x_j}
\tag{9} \label{eq:scale-factor-9}
\end{align}
$$
これでどうでしょうか？ニュートンの粘性法則を基にした、速度起因の面積力を表現できたと思います。
しかし、本当にこれでいいのでしょうか？対称テンソルの議論で、
$$
\begin{align}
&\tau_{ij} = \tau_{ji}
\tag{10} \label{eq:scale-factor-10}
\end{align}
$$
でなければならなかったのに対し、このままでは、
$$
\begin{align}
&\frac{\partial u_i}{\partial x_j} \neq \frac{\partial u_i}{\partial x_j}
\tag{11} \label{eq:scale-factor-11}
\end{align}
$$
となってしまっています。たまたま一致することはあるかもしれませんが、常に一致しなければなりません。
そこで、こんなことをしてみます。
$$
\begin{align}
&\tau_{ij} = \mu ( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} )
\tag{12} \label{eq:scale-factor-12}
\end{align}
$$
ニュートンの粘性法則は物理的なモデリングですが、これは完全なモデリングです。
とりあえず、簡単だしこれを使ってみたら、実際の速度と合っているのでみんなこれを使っています。
同じモデリングなら、例えば以下の式でもいいんじゃないでしょうか？対称テンソルを満たしていますよね？
$$
\begin{align}
&\tau_{ij} = \mu ( \frac{\partial u_i}{\partial x_j} \cdot \frac{\partial u_j}{\partial x_i} )^{\frac{1}{2}}
\tag{13} \label{eq:scale-factor-13}
\end{align}
$$
しかし、掛け算だとダメな点があります。
上図のような速度勾配を考えると、
$$
\begin{align}
&\frac{\partial u_1}{\partial x_3} \gg \frac{\partial u_3}{\partial x_1}
\tag{14} \label{eq:scale-factor-14}
\end{align}
$$
となっていると思います。この条件を基に、\eqref{eq:scale-factor-12}を考えると、
$$
\begin{align}
&\tau_{31} = \tau_{13} = \mu \frac{\partial u_1}{\partial x_3}
\tag{15} \label{eq:scale-factor-15}
\end{align}
$$
となり、だいたいニュートンの粘性法則と一致するのに対し、\eqref{eq:scale-factor-13}式を使うと、ほぼ0になってしまい、実際との乖離が大きくなります。なので、足し算系で行かないとダメそうです。
そして、その足し算で一番シンプルな\eqref{eq:scale-factor-12}式を使っています。モデリングは計算の面でも理解の面でも出来るだけシンプルなのがいいのです。
これにて、面積力が完成し、Navier-Stocks方程式はほぼ完成します。
ここまでのモデリングを元の式にぶちこみます。
$$
\begin{align}
&\frac{ \partial \rho u_i }{ \partial t} =- \frac{\partial \rho u_i u_j}{\partial x_j} + \frac{\partial }{\partial x_j}(- P \delta_{ij}) + \frac{\partial} {\partial x_j}\{ \mu (\frac{\partial u_i}{\partial x_j} +\frac{\partial u_j}{\partial x_i}) \} – \rho g \delta_{i3} \\\\
&\frac{ \partial \rho u_i }{ \partial t} =- \frac{\partial \rho u_i u_j}{\partial x_j} \ – \frac{\partial P}{\partial x_i} + \frac{\partial} {\partial x_j}\{ \mu (\frac{\partial u_i}{\partial x_j} +\frac{\partial u_j}{\partial x_i}) \} – \rho g \delta_{i3}
\tag{16} \label{eq:scale-factor-16}
\end{align}
$$
となります。一応解説しておくと、一番最初の記事くらいに書いてある、クロネッカーのデルタの置き換えテンソルの性質を使いました。
$$
\begin{align}
&\frac{\partial }{\partial x_j}(- P \delta_{ij}) = – \frac{\partial P}{\partial x_i}
\tag{17} \label{eq:scale-factor-17}
\end{align}
$$
で、クロネッカーのデルタで縮約をとってやると結局添字が変わるというものです。忘れた人は一番最初の流体の前に基礎知識という記事をみてみてください。
$$
\begin{align}
&\frac{\partial }{\partial x_j}(- P \delta_{ij}) = – \frac{\partial P}{\partial x_i}
\tag{18} \label{eq:scale-factor-18}
\end{align}
$$
$$
\begin{align}
&a_i \delta_{ij} = a_j
\tag{19} \label{eq:scale-factor-19}
\end{align}
$$
というやつです。これでやっと終わりかーと思いきやもう少し続きます。長くなったので、次の記事に回します。