初めに

お味噌汁が好きな理由は様々ある．

美味しさは勿論のこと，統計学（cf. みどりの窓口の枕），熱対流を感じられることなど．

Methods

支配方程式

熱と流れをそれぞれ独立に考えたことはあるが（e.g., heat conduction in a rod, shear-driven flow in a square / cubic cavity），同時に考えたことはない．

先ず，それぞれの方程式は以下．

\begin{align}
    \partial_t \rho
    + \nabla \cdot \left( \rho \boldsymbol{u} \right)
    &= 0
    \\
    \rho
    \left(
        \partial_t \boldsymbol{u}
        + \left( \boldsymbol{u} \cdot \nabla \right) \boldsymbol{u}
    \right)
    &= - \nabla p
    + \mu \nabla^2 \boldsymbol{u}
    + \rho \boldsymbol{f}
    \\
    \rho c_p \left(
        \partial_t T
        + \left( \boldsymbol{u} \cdot \nabla \right) T
    \right)
    &= \lambda \nabla^2 T
    + \rho c_p s
\end{align}

ここで，$\boldsymbol{u}$は流速$[ \mathrm{m / s} ]$，$p$は圧力$[ \mathrm{Pa} ]$，$\boldsymbol{f}$は物体力$[ \mathrm{m / s^2} ]$であり，これには重力を与える．また，運動に関するパラメータとして，$\rho$は密度$[ \mathrm{kg / m^3} ]$，$\mu$は粘度$[ \mathrm{Pa \cdot s} ]$である．

さらに，$T$は温度$[ \mathrm{K} ]$，$s$は熱源$[ \mathrm{K / s} ]$であり，ゼロとする．なお，熱に関するパラメータとして，$c_p$は比熱容量$[ \mathrm{J / (kg \cdot K)} ]$，$\lambda$は熱伝導率$[ \mathrm{W / (m \cdot K)} ]$である¹．

いま，密度$\rho$を基準点$\rho_0$周りでTaylor展開する．

\begin{align}
    \rho
    &= \rho_0
        + \left. \frac{\partial \rho}{\partial T} \right|_{p} ( T - T_0 )
        + \left. \frac{\partial \rho}{\partial p} \right|_{T} ( p - p_0 )
        + (\text{HOT})
    \\
\end{align}

但し，密度が圧力に依存して変化することはないとし，温度のみの関数として取り扱う（Boussinesq近似，$\rho \left( T, p, ... \right) \simeq \rho \left( T \right)$²）．したがって

\begin{align}
    \rho
    &\simeq \rho_0
        + \left. \frac{\partial \rho}{\partial T} \right|_{p} ( T - T_0 )
    \\
    &= \rho_0 + \rho^{\prime}
        \quad \left(
            \rho^{\prime}
            := \left. \frac{\partial \rho}{\partial T} \right|_{p} ( T - T_0 )
        \right)
\end{align}

ここで，$\rho^{\prime}$は基準点からの変動であり，$\rho^{\prime} \ll \rho_0$．熱膨張率$\beta \ [ \mathrm{1 / K} ]$を

\begin{align}
    \beta
    &= - \frac{1}{\rho_0} 
        \left. \frac{\partial \rho}{\partial T} \right|_{p}
\end{align}

と定めると，

\begin{align}
    \rho^{\prime}
    &= \left. \frac{\partial \rho}{\partial T} \right|_{p} \left( T - T_0 \right)
    \\
    &= - \rho_0 \beta \left( T - T_0 \right)
    \\
    \therefore
    \frac{\rho^{\prime}}{\rho_0}
    &= - \beta \left( T - T_0 \right)
\end{align}

である．いま，外力$\boldsymbol{f}$に重力$\boldsymbol{g}$を与えたことを思い出せば，

\begin{align}
    \frac{\rho}{\rho_0} \boldsymbol{f}
    &= \left( 1 + \frac{\rho^{\prime}}{\rho_0} \right) \boldsymbol{g}
    \\
    &= \left( 1 - \beta \left( T - T_0 \right) \right) \boldsymbol{g}
\end{align}

であるから，

\begin{align}
    \nabla \cdot \boldsymbol{u}
    &= 0
    \\
    \partial_t \boldsymbol{u}
    + \left( \boldsymbol{u} \cdot \nabla \right) \boldsymbol{u} 
    &= - \frac{1}{\rho_0} \nabla p 
    + \frac{\mu}{\rho_0} \nabla^2 \boldsymbol{u} 
    + \left( 1 - \beta \left( T - T_0 \right) \right) \boldsymbol{g}
    \\
    \partial_t T    
    + \left( \boldsymbol{u} \cdot \nabla \right) T
    &= \kappa \nabla^2 T
\end{align}

を解くこととする．但し，$\kappa$は熱拡散率$[ \mathrm{m^2 / s} ]$である（$\kappa = \lambda / \rho c_p$）．

Results

空間微分をArakawa-B grid上で3次+2次精度で，時間微分を1次精度で近似した．離散化パラメータは$h = 2\times10^{-2} \text{ [m]}$と$\tau = 2\times10^{-3} \text{ [s]}$．

Shear-driven flow

先ずは，ただのshear-driven flowを長方形領域上で計算する．速度ノルム$| \boldsymbol{u} |$と，動圧で正規化した圧力$q = p / (1/2 \rho | \boldsymbol{u} |^2)$を示す．

$\text{Re}$	$\text{Results}$
$\sim 10^2$
$\sim 10^3$
$\sim 10^4$

squareの自然な拡張となっており良し．

Buoyancy-driven flow

深さ$1 \ [ \mathrm{m} ]$の鍋を水で満たし，底面を沸点$100 \ [ ^{\circ}\mathrm{C} ] = 373.15 \ [ \mathrm{K} ]$で温め，表面を室温$20 \ [ ^{\circ}\mathrm{C} ] = 293.15 \ [ \mathrm{K} ]$に解放する．その他，パラメータはこちらから借用する．

密度　　　: $\rho_0 = 1.00 \times 10^{3} \ [ \mathrm{kg / m^3} ]$
粘性率　　: $\mu = 1.00 \times 10^{-3} \ [ \mathrm{Pa \cdot s} ]$
動粘性率　: $\nu = \mu / \rho_0 = 1.00 \times 10^{-6} \ [ \mathrm{m^2 / s} ]$
熱伝導率　: $\lambda = 0.60 \ [ \mathrm{W / (m \cdot K)} ]$
比熱容量　: $c_p = 4.18 \times 10^{3} \ [ \mathrm{J / (kg \cdot K)} ]$
熱拡散率　: $\kappa = \lambda / \rho c_p = 1.44 \times 10^{-7} \ [ \mathrm{m^2 / s} ]$
体積膨張率: $\beta = 2.10 \times 10^{-4} \ [ \mathrm{1 / K} ]$
重力加速度: $g = 9.81 \ [ \mathrm{m / s^2} ]$
Prandtl数 : $\mathrm{Pr} = \nu / \kappa = 6.97$
Grashof数 : $\mathrm{Gr} = g \beta L^3 \Delta T / \nu^2 = 1.65 \times 10^{11}$
Rayleigh数: $\mathrm{Ra} = \mathrm{Pr} \cdot \mathrm{Gr} = 1.15 \times 10^{12}$

先ず，以下は，

流速: 全ての壁面で非滑り
圧力: 全ての壁面で斉次Neumann
温度: 上下の壁面に温度差，その他全ての壁面で室温

としたときの結果．

$\text{Dim}$	$\text{Results}$
$\text{2D}$
$\text{3D}$

次いで，以下は，

流速: 全ての壁面で非滑り
圧力: 全ての壁面で斉次Neumann
温度: 上下の壁面に温度差，その他全ての壁面で斉次Neumann

としたときの結果．

$\text{Dim}$	$\text{Results}$
$\text{2D}$
$\text{3D}$

最後に，以下は，

流速: 上下の壁面で非滑り，その他全ての壁面で周期的
圧力: 上下の壁面で斉次Neumann，その他全ての壁面で周期的
温度: 上下の壁面に温度差，その他全ての壁面で周期的

としたときの結果．

$\text{Dim}$	$\text{Results}$
$\text{2D}$
$\text{3D}$

終わりに

お盆休み，僅かばかりの遊び心．

$h$	$1 \times 10^{-2} \text{ [m]}$	$5 \times 10^{-3} \text{ [m]}$	$2 \times 10^{-3} \text{ [m]}$	$1 \times 10^{-3} \text{ [m]}$
-
-

鶏冠を立ち上げたニワトリか，湯煎されたチョコレートのよう．その昔，高校の地理で，ペルー・チリ沖で起きるという湧昇流なるものを習ったことを思い出す．また，高解像度の解を眺めていると，木星の大赤斑が思い浮かぶ．

そう言えば，密度差を陽な浮力としているから，Rayleigh-Taylor instabilityのような構造も当然生まれる（というか，普通に体積力の存在）．

1	2	3	4

Appendix

無次元量

幾つかの無次元量を導入する．

Prandtl数

動粘性と熱拡散の比．Ludwig Prandtlに由来．

\begin{align}
    \mathrm{Pr}
    &= \frac{\nu}{\kappa}
     = \frac{\mu / \rho}{\lambda / \rho c_p}
     = \frac{\mu c_p}{\lambda}
\end{align}

ここで，

$\mu$: 粘性率$[ \mathrm{Pa \cdot s} ]$
$\nu$: 動粘性率$[ \mathrm{m^2 / s} ]$
$\rho$: 密度$[ \mathrm{kg / m^3} ]$
$\kappa$: 熱拡散率$[ \mathrm{m^2 / s} ]$
$\lambda$: 熱伝導率$[ \mathrm{W / (m \cdot K)} ]$
$c_p$: 比熱容量$[ \mathrm{J / (kg \cdot K)} ]$

である．

Grashof数

浮力と粘性力の比．Franz Grashofに由来．

\begin{align}
    \mathrm{Gr}
    &= \frac{g \beta L^3 \Delta T}{\nu^2}
\end{align}

ここで，

$g$: 重力加速度$[ \mathrm{m / s^2} ]$
$\beta$: 体積膨張率$[ \mathrm{1 / K} ]$
$\nu$: 動粘性率$[ \mathrm{m^2 / s} ]$
$L$: 特徴的な長さ$[ \mathrm{m} ]$
$\Delta T$: 特徴的な温度差$[ \mathrm{K} ]$

である．

Rayleigh数

浮力と熱拡散の比．John William Strutt, 3rd Baron Rayleighに由来．

\begin{align}
    \mathrm{Ra}
    &= \frac{g \beta L^3 \Delta T}{\kappa \nu}
     = \mathrm{Pr} \cdot \mathrm{Gr}
\end{align}

ここで，

$g$: 重力加速度$[ \mathrm{m / s^2} ]$
$\beta$: 体積膨張率$[ \mathrm{1 / K} ]$
$\kappa$: 熱拡散率$[ \mathrm{m^2 / s} ]$
$\nu$: 動粘性率$[ \mathrm{m^2 / s} ]$
$L$: 特徴的な長さ$[ \mathrm{m} ]$
$\Delta T$: 特徴的な温度差$[ \mathrm{K} ]$

である．

「熱拡散率（thermal diffusivity）」と「熱伝導率（thermal conductivity）」と「熱伝達率（heat transfer coefficient）」は，全て異なる．熱拡散率$\kappa \ [ \mathrm{m^2 / s} ]$は，熱伝導率$\lambda \ [ \mathrm{W / (m \cdot K)} ]$，比熱容量$c_p \ [ \mathrm{J / (kg \cdot K)} ]$，密度$\rho \ [ \mathrm{kg / m^3} ]$を用いて，$\kappa = \lambda / \rho c_p$で与えられる．熱伝達率$\eta \ [ \mathrm{W / (m^2 \cdot K)} ]$は，壁面と流体などの2種類の物質間での熱エネルギーの伝わり易さである．熱移動量$Q \ [ \mathrm{W} ]$，熱流束密度$J \ [ \mathrm{W / (m^2)} ]$，電熱面積$A \ [ \mathrm{m^2} ]$，壁面の温度$T_0 \ [ \mathrm{K} ]$，媒質の温度$T \ [ \mathrm{K} ]$を用いて，$\eta = Q / A (T_0 - T) = J / (T_0 - T)$で与えられる．但し，$T_0 \gt T$． ↩
非圧縮性近似とは，密度$\rho$が圧力$p$に依存して変化することはない，と仮定するものである．直ちに$\rho = (\text{const.})$を言うものではない（参考）． ↩