自然対流計算におけるブーシネスク近似の拡張

By M.Sato

June 1, 2024

概要

今回は矩形の閉鎖領域内における自然対流の計算を紹介する。

対象とする流体は水である。水密度は温度が4[℃]付近で最大となるが、水温が4[℃]をまたがって変動する場合を計算した。

熱膨張係数の符号が変化することになるため、ここではブーシネスク近似モデルを拡張し、浮力効果を体積力として考慮することで流れ場、温度場を計算するようにした。

また、ユーザ定義関数を使わずに名前付き式（expression language）を使って計算する方法を示した。

解析対象

解析対象は $38mm \times 38mm$の矩形領域である。

左側が高温(8℃)壁、右側が低温(0℃)壁となっている。４側面は全て速度ゼロとした。

解析領域の中心位置にモニター点を設置した。

geom

解析モデル

メッシュ分割

格子状にメッシュを作成した。

mesh

物性値

流体は、低温の水である。

左右の壁面温度が0℃、8℃であることから平均を取り、4[℃]を基準として物性値を設定した。

なお、今回の計算は浮力を考慮するため密度については温度依存性を考慮する必要がある。

ここではNIST(米国立標準技術研究所)のデータベースに基づき、密度データを取得した。

浮力の実装方法については後述する。

物性値	値（４℃を想定）
密度（基準値）	$999.974 [kg/m^3]$
粘度	$1.567\times 10^{-3} [Pa.s]$
比熱	$4.207\times 10^{3} [J/kg/K]$
熱伝導率	$5.655\times 10^{-1} [W/m/K]$

計算に使用する密度は下図に示すような温度関数とした。温度が4℃付近で最も密度が大きくなっている。

なお計算での扱いを容易にするために、密度の温度依存性を３次式で回帰するようにした。

温度の単位は絶対温度[K]である（fluentではソルバーが使用する単位系は全てSI単位になることに注意）。

回帰式

$$ \rho(T) = 8.8751\times 10^{-5} T^3 -8.1777\times 10^{-2} T^2 +2.4877\times 10^1 T -1.5026\times 10^3 \quad [kg/m^3] $$

density-vs-temperature

NISTからのデータ取得方法

以下に、NISTから熱物性データを取得する手順を示す。

NISTホームページにアクセスする

NIST

流体名、等温/等圧条件の指定、単位系の指定などを行う

想定圧力、温度範囲/刻みの指定、表示桁数の指定などを行う。

物性値データの表示

物性値データのダウンロード

ダウンロードしたデータから必要な情報を取り出す

浮力のモデル化

重力の影響下にある流体の場合、流体密度が大きい箇所では下降する流れ、流体密度が小さい箇所では上昇する流れが発生する。

この流れは、運動方程式に浮力の効果を考慮することで再現することができる。

一般的に密度変化を伴う流れの基礎方程式は次式で表される。

$$ \begin{align} &\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x_i}(\rho u_i) = 0 \\ &\frac{\partial \rho u_i}{\partial t} + \frac{\partial}{\partial x_j}(\rho u_i)=0 \\ &\frac{\partial \rho c_p T}{\partial t}+ \frac{\partial}{\partial x_j}( \rho c_p u_j T)= k \frac{\partial^2 T}{\partial x_j \partial x_j} \end{align} $$

運動方程式(2)式の右辺の圧力勾配項と重力項は以下のように表すことができる。

$$ \begin{align} -\frac{\partial p}{\partial x_i}+ \rho g_i &= -\frac{\partial p}{\partial x_i}+ (\rho_0 + (\rho-\rho_0)) g_i \\ &=-\frac{\partial (p-\rho_0 g_i x_i)}{\partial x_i} + (\rho-\rho_0)g_i \end{align} $$

密度 $\rho(T)$をテイラー展開すると次式が得られる。

$$ \rho(T) = \rho(T_0) + \frac{\partial \rho}{\partial T}(T-T_0)+\frac{1}{2}\frac{\partial^2 \rho}{\partial T^2}(T-T_0)^2+… $$

$\rho_0 = \rho(T_0)$で両辺を割り、$\Delta T=T-T_0$と書いて整理すると

$$ \frac{\rho(T)-\rho_0}{ \rho_0} = \frac{1}{\rho_0}\frac{\partial \rho}{\partial T}\Delta T +\frac{1}{2\rho_0}\frac{\partial^2 \rho} {\partial T^2} \Delta T^2+… \simeq - \beta \Delta T $$

ここで $\beta=-\frac{1}{\rho_0}\frac{\partial \rho}{\partial T}$は熱膨張係数である。

$\beta$を一定値として密度を温度の線形関数として近似すればブーシネスク近似となる。

この場合、運動方程式は次のようになる。

$$ \frac{\partial \rho_0 u_i}{\partial t} + \frac{\partial}{\partial x_j}\left( \rho_0 u_i u_j \right)= -\frac{\partial \hat{p}}{\partial x_i} -\rho_0 g_i \beta (T-T_0) +\mu \frac{\partial^2 u_i}{\partial x_j \partial x_j} \tag{4} $$

ここで $\hat{p}=p-\rho_0 g_i x_i$である。

今回の解析対象は密度が温度により大きく増減(熱膨張係数$\beta$の符号が逆転)するため、浮力項を単純な温度の線形関数として扱うのは不適切である。

したがって、$\beta$を使わずに直接、浮力項 $(\rho-\rho_0)g_i$ を体積力として近似する。浮力項以外の密度は一定値とする。

以下に今回の計算で使用する基礎方程式を示す。

$$ \begin{align} &\frac{\partial u_i}{\partial x_i}=0 \tag{5} \\ &\rho_0\left( \frac{\partial u_i}{\partial t} + \frac{\partial u_i u_j}{\partial x_j}\right)= -\frac{\partial \hat{p}}{\partial x_i}+ (\rho-\rho_0) g_i +\mu \frac{\partial^2 u_i}{\partial x_j \partial x_j} \tag{6} \\ &\rho_0 c_p\left( \frac{\partial T}{\partial t} + \frac{\partial u_j T}{\partial x_j} \right)= k \frac{\partial^2 T}{\partial x_j \partial x_j} \tag{7} \end{align} $$

無次元パラメータの評価

浮力流れの特性は$R_a$（レーリー）数で評価される。

平均の熱膨張係数（$0$～$4$℃）を以下の式で近似する。

$$ \beta =\frac{1}{\rho_0}\frac{\Delta \rho}{\Delta T} = \frac{1}{999.974}\frac{999.974-999.843}{4-0}=3.275\times 10^{-5}[K^{-1}] $$

代表温度差$\Delta T=4$℃、代表長さ$L=0.038$m とすると、動粘性係数 $\nu=\mu/\rho_0=1.567\times 10^{-3}/999.974=1.567\times 10^{-6}[m^2/s]$、熱拡散係数 $\alpha=\mu c_p/k=1.567\times 10^{-6}\times4207/0.5655=11.66[m^2/s]$であることから

$$ R_a = \frac{ g \beta \Delta T L^3}{\nu \alpha} = \frac{ 9.81 \times 3.275\times 10^{-5} \times 4 \times 0.038^3}{1.567\times 10^{-6}\times 11.66} =3.35\times 10^{5} $$

$R_a$数は$10^9$付近で乱流に遷移すると思われるため、ここでは層流として計算する。