流体の熱エネルギー保存則(1) ブジネスク近似、熱流体の取り扱い

非圧縮性流体

熱流体の扱い方

前回でやっと、「非圧縮等温流体」の運動量保存則が完成しました。一応書いておきます。

美しいですね。流体の挙動を表す式が、こんなシンプルな方程式にまとまるなんて。これを作ったNavierさんとStocksさんは天才ですね。ちなみに、Navierがほとんどを作っていたのですが、(分子)粘性項のところがうまくいかず、最後にそこを完成させたのがStocksのようです。

$$
\begin{align}
&\frac{ \partial u_i }{ \partial t} + \frac{\partial u_i u_j}{\partial x_j} = \ – \frac{1}{\rho} \frac{\partial P}{\partial x_i} + \nu \frac{\partial^2 u_i}{\partial x_j^2}
\tag{1} \label{eq:scale-factor-1} 
\end{align}
$$

さて、ここからはこの非圧縮の方程式に「熱」を組み込み、熱流体にする方法を考えていきます。

先に全体像を示すために、熱流体を扱う方法は3つです。

  1. 非圧縮のNavier-Stocks方程式に、そのまま熱エネルギーの保存則を加えて解く
  2. 非圧縮のNavier-Stocks方程式に、熱エネルギーの保存則を加え、熱による影響をNavier-Stocks方程式に浮力項としてフィードバックさせ解く
  3. 圧縮性の方程式で解く

まず、1.はよく機械設計で使われます。強制対流と呼ばれ、機械の力で強制的に流体を動かす時などに使います。機械で強制的に流体を動かすので、熱による浮力の影響など無視されるから、そのまま並列して解けばいいよね、という考え方です。

2.は、室内気流などのあまり動かない流体に対して、非圧縮のまま浮力を考慮するために使われます。この方法をBoussinesq近似(ブジネスク近似)と言います。理解するためのいいトレーニングになるので、ブジネスク近似を主に勉強していきます。

3.はそのままです。基本的に圧縮性の方程式は、圧縮性用の質量保存則(連続式)、運動量保存則に加え、熱力学ベースの、エネルギー保存則と、状態方程式を並列させて解きます。圧縮性では、圧力や熱(内部エネルギー)に対して密度が変わることを考慮して解いていきます。考え方としてはこちらの方がシンプルではあるのですが、ケースバイケースでいろんな処理や方程式を変える必要があり、難しさや計算コストは大きくなります。

ただ、あらゆる場面に対応できるので、いずれ圧縮性のところでやっていきます。


Boussinesq近似

Boussinesq近似とは何かと説明されれば、以下の2点です。

  • 非圧縮性を仮定
  • 密度変化は温度変化として重力項のみで考慮

まずは適用範囲に付いて考えていきます。密度が温度によって一定とみなせる範囲は、30℃差以内と言われています。つまり、Boussinesq近似の適用範囲は、温度差30℃以内です。それ以上だと圧縮性流体で解きます。

また、流速についてもついでに考えます。これはそもそもBoussinesq近似するかどうかではなくて、非圧縮性か圧縮性かなのですが、基本的に マッハ数 = 流速 / 音速 <= 0.3 くらいが非圧縮で解いていい範囲と言われています。音速はだいたい30℃くらいで340m/s くらいなので、 だいたい流速が100m/sくらいからでしょうか。

と言った感じで、自分が解く問題に対して方程式を変えて解く必要があるわけです。それでも非圧縮性はほぼ同じ方程式で解けるので、便利なものです。

 


温度と密度のモデリング

まず、温度と密度の関係式を作ってモデリングすることを考えます。

物理現象として、密度と温度の関係式と言えばボイル・シャルルとか状態方程式ですね。状態方程式を考えてみます。

$$
\begin{align}
&P = \rho r T
\tag{2} \label{eq:scale-factor-2} 
\end{align}
$$

すでに知っているのと少し違うかもしれませんが、こっちが流体でよく使う状態方程式です。
一応PV = nRTからの変形までやっておきます。

$$
\begin{align}
& n = \frac{w}{M}より、 \\\\
P V &= nRT \\\\
& = \frac{w}{M} RT \\\\
& = \frac{w}{V} \cdot \frac{1}{M} RT \\\\ 
& = \rho \frac{R}{M} T \\\\
ここで、&\frac{R}{M} = rとして、\\\\
P &= \rho r T
\tag{3} \label{eq:scale-factor-3} 
\end{align}
$$

とまあこんな感じのを圧縮性でよく使うので覚えておいてください。ここでの気体定数rは、Mが入っているので、気体の種類によって変わります。

 

この関係を踏まえ、モデリングしていきます。まず、足し算で簡単にモデリングして関係式を作れないかと考えてみます。

$$
\begin{align}
& \rho – \rho_0 = \hspace{15pt} (T – T_0) + \hspace{15pt}(P-P_0)
\tag{4} \label{eq:scale-factor-4} 
\end{align}
$$

しかし、このままだと単位が違うので、同じ物理量で割って無次元化を行って作ります。

$$
\begin{align}
& \frac{1}{\rho_0}(\rho – \rho_0) = \frac{1}{T_0} (T – T_0) + \frac{1}{P_0}(P-P_0)
\tag{5} \label{eq:scale-factor-5} 
\end{align}
$$

こうすれば、単位もオーダーも揃います。しかし、状態方程式を考えると\(\rho\)とTは反比例で、Pとは比例関係なのでこのままではいけません。

$$
\begin{align}
& \frac{1}{\rho_0}(\rho – \rho_0) = \ – \ \frac{1}{T_0} (T – T_0) + \frac{1}{P_0}(P-P_0)
\tag{6} \label{eq:scale-factor-6} 
\end{align}
$$

そして、圧力はここでは考慮しないので無視します。つまり、

$$
\begin{align}
& \frac{1}{\rho_0}(\rho – \rho_0) = \ – \ \frac{1}{T_0} (T – T_0) = \ – \ \beta (T – T_0)
\tag{7} \label{eq:scale-factor-7} 
\end{align}
$$

ここで、\(\beta = \frac{1}{T_0}\)を体膨張率と言います。

こんな簡単な形ですが、密度と温度のモデリング完了です。

 


運動量保存則に組み込み

ここから、この影響を運動量保存則に組み込んでいきます。

まず、圧力項を重力項を吸収する前の形の運動量保存則に戻します。

$$
\begin{align}
&\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{8} \label{eq:scale-factor-8} 
\end{align}
$$

ここで基準密度を\(\rho_0\)として、重力項を変形します。

$$
\begin{align}
& \rho g \delta_{i3} = ( \rho – \rho_0) g \delta_{i3} + \rho_0 g \delta_{i3}
\tag{9} \label{eq:scale-factor-9} 
\end{align}
$$

ここで、密度変化をしないという近似を入れて、各項の\(\rho\)を\(\rho_0\)にします。

$$
\begin{align}
&\frac{ \partial \rho_0 u_i }{ \partial t} + \frac{\partial \rho_0 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 – \rho_0) g \delta_{i3} \ – \rho_0 g \delta_{i3}
\tag{10} \label{eq:scale-factor-10} 
\end{align}
$$

そして、また\( – \rho_0 g \delta_{i3}\)を圧力項に戻します。

$$
\begin{align}
&\frac{\partial P}{\partial x_i}  = \frac{\partial P’}{\partial x_i} \ – \rho_0 g \delta_{i3}
\tag{11} \label{eq:scale-factor-11} 
\end{align}
$$

そして、\eqref{eq:scale-factor-7} 式を使い、

$$
\begin{align}
&( \rho – \rho_0) \ g \ \delta_{i3} = \ – \rho_0 \  \beta (T – T_0) g \ \delta_{i3}
\tag{12} \label{eq:scale-factor-12} 
\end{align}
$$

と変形します。そして、\(\rho_0\)で全体を割り、\(\rho_0  \rightarrow  \rho\) 、\(P’  \rightarrow  P\)、\((T – T_0)  \rightarrow \Delta T \)と置き換えます。

$$
\begin{align}
&\frac{ \partial u_i }{ \partial t} + \frac{ \partial u_i u_j }{ \partial x_j} = -\frac{1}{\rho} \frac{\partial P}{\partial x_i} + \nu \frac{\partial^2 u_i}{\partial x_j^2} + g\beta \Delta T \delta_{i 3} \tag{13} \label{eq:scale-factor-13}
\end{align}
$$

これで、運動量保存則の修正が終わりました。

ここでの、\(\rho_0\)と\(T_0\)は\(\rho_0\)を考えた時のものであるべきです。これらはどーいった値にすればいいかと言うと、なんでもいいです。まじでなんでもいいです。

代表値をとればいいです。この時、\(T_0\)を決めれば\(\rho_0\)も決まるので、\(T_0\)を考えて決めます。

これで、非圧縮性非等温流体の方程式が完成です。最後の項を浮力項と言います。

ここでの浮力項の中に出てくるgは、\(g = +9.8 m/s\)を使っています。

ベクトル表記するとこうなります。

$$
\begin{align}
&\frac{ \partial \mathbf{V}}{ \partial t} + (\mathbf{V} \cdot \nabla) \mathbf{V} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{V} –  \beta \Delta T \mathbf{g} \tag{14} \label{eq:scale-factor-14}
\end{align}
$$

最後の浮力項の正負が変わっていることに注意してください。ベクトル表記の場合、重力ベクトルを使って表記するので、重力の方向と逆の方向に浮力は働きますので、マイナスが付いています。



[
前の記事へ]  [非圧縮性流体の目次へ]  [次の記事へ]