[論文レビュー] 30 K でブシネスクが嘘をつくとき — Demou(2018) NOB自然対流と圧力分割

Gray–Giorginiの水の限界1.25 Kに対しΔT=30 Kは24倍で、密度補正は1%台に残ります。

Gray と Giorgini は 1976 年に数字を二つ残しました。常温の空気で Oberbeck–Boussinesq が 10% 誤差の内にいる最大温度差は 28.6 K です。水は 1.25 K です。太陽熱タワーの空気は数百 K を行き来し、PV/T の水は 50 K 温まります。Demou・Frantzis・Grigoriadis(2018)はその仮定を捨てます。物性は温度に従い、圧力 Poisson だけが定係数のままです。

hot 壁と cold 壁を動かすと、同じ ΔT で密度の曲線はほとんど寝て、粘性の曲線だけが大きく傾きます。

T_hot: 333 K (60 °C)T_cold: 293 K (20 °C)

Dashed line marks the reference value at T₀ = 313 K. Δ% is the relative change between hot and cold. For ΔT ≈ 50 K, density moves about 2%, but viscosity changes by tens of percent — the OB approximation hides that.

50 K で密度は約 2% です。粘性は同じ区間で約 30% 動き、その差が非対称な境界層を作ります。

論文メモ

著者: A. D. Demou, C. Frantzis, D. G. E. Grigoriadis
誌名: International Journal of Heat and Mass Transfer 121 (2018) 1148–1162
DOI: 10.1016/j.ijheatmasstransfer.2018.04.144
寄与: 非 OB(NOB)自然対流で物性を温度依存のまま保ち、可変係数の圧力 Poisson を Dodd–Ferrante 分割で定係数 Poisson と明示的補正に変えます。FFT 直接ソルバがそのままで、報告された加速は 2.5–5 倍です。

Boussinesqが破れる場所#

標準 Boussinesq は三つを同時に置きます。密度は重力項だけで変わります。他の物性は定数です。粘性散逸は無視します。ΔT が大きくなると粘性 $\mu(T)$ 、熱伝導率 $k(T)$ 、比熱 $C_p(T)$ が一緒に動きます。液体は密度の変化が小さくても、粘性の変化は大きいです。

物性をすべて残すと、運動量方程式に可変係数が入ります。

\frac{\partial u_i}{\partial t} + \frac{\partial(u_i u_j)}{\partial x_j} = -\frac{1}{\hat\rho}\frac{\partial P}{\partial x_i} + \frac{1}{\hat\rho}\frac{\Pr}{\sqrt{\mathrm{Ra}}}\frac{\partial}{\partial x_j}\!\left[\hat\mu\!\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)\right] + \frac{1}{\mathrm{Fr}^2}\delta_{i3}

$\hat\rho, \hat\mu, \hat C_p, \hat k$ は基準温度の値で正規化した無次元物性です。同じ可変性はエネルギー式にも入ります。

可変係数Poissonの罠#

分数ステップ法をそのまま使うと、圧力補正は発散ゼロを次の形で要求します。この補正は Chorin が圧力を運動量から消したあとに来る段です。

\frac{\partial}{\partial x_i}\!\left(\frac{1}{\hat\rho^{n+1}}\frac{\partial P^{n+1}}{\partial x_i}\right) = \frac{1}{\Delta t}\frac{\partial u_i^*}{\partial x_i}.

左辺に $1/\hat\rho^{n+1}$ が掛かっています。係数が空間と時間で変わると、行列も毎ステップ変わります。FFT 直接ソルバは定係数ラプラスしか受けません。結局、多重格子か反復ソルバに落ちます。NS ソルバの時間の 60–80% が圧力なら、その落ち方が請求書です。

Dodd–Ferrante分割 — 一行の置換#

借りてきた工夫は二流体(Dodd & Ferrante 2014)からです。圧力勾配を二つに割ります。

\frac{1}{\hat\rho^{n+1}}\frac{\partial P^{n+1}}{\partial x_i} \longrightarrow \frac{1}{\hat\rho_0}\frac{\partial P^{n+1}}{\partial x_i} + \left(\frac{1}{\hat\rho^{n+1}} - \frac{1}{\hat\rho_0}\right)\frac{\partial \widehat P}{\partial x_i}

第一項は定数 $1/\hat\rho_0$ が圧力に掛かります。第二項は予測圧力 $\widehat P = 2P^n - P^{n-1}$ に掛ける補正です。発散を取ると Poisson は定係数になります。

\frac{\partial^2 P^{n+1}}{\partial x_i \partial x_i} = \frac{\hat\rho_0}{\Delta t}\frac{\partial u_i^*}{\partial x_i} + \frac{\partial}{\partial x_i}\!\left[\left(1 - \frac{\hat\rho_0}{\hat\rho^{n+1}}\right)\frac{\partial \widehat P}{\partial x_i}\right]

左辺のラプラスは毎ステップ同じです。二方向の FFT のあと、z 方向の Helmholtz を直接解きます。右辺の補正だけを陽的に更新します。

大きさを測る量は $(1 - \hat\rho_0/\hat\rho)$ です。二流体では 1 に近く、NOB の水対流ではとても小さいです。

NumPyで一サイクル#

論文 §2.2–2.3 を 1D の RB スラブ一行に縮めました。2D DNS はなく、圧力補正を一歩だけ見ます。

import numpy as np
from numpy.fft import fft, ifft
 
# 格子: x周期、z壁境界 (簡略化のため1D Helmholtz一束のみ)
Nz = 64
Lz = 1.0
dz = Lz / Nz
z = (np.arange(Nz) + 0.5) * dz
 
# 温度依存密度 (水polynomial簡略化)
def rho_hat(T_field, T0=313.0):
    dT = T_field - T0
    return 1.0 * (1 - 3.736e-4 * dT - 3.98e-6 * dT**2)
 
# 初期温度場: 線形 + 小さな摂動
T0 = 313.0
dT = 30.0
T_field = T0 - dT/2 + dT * z + 0.1 * np.sin(2*np.pi*z)
rho = rho_hat(T_field)
 
# 基準密度: 最も軽い (最も熱い) 側
rho_ref = rho.min()
 
# 仮の発散場 (実際は分数ステップ後の∇·u*から)
div_u_star = 0.01 * np.sin(2*np.pi*z)
dt = 1e-3
 
# 予測圧力 (論文 Eq. 15)
P_n_minus_1 = np.zeros(Nz)
P_n = 0.01 * np.cos(2*np.pi*z)
P_hat = 2 * P_n - P_n_minus_1
dP_hat = np.gradient(P_hat, dz)
 
# 補正項 (Eq. 16 右辺第二)
correction = (1 - rho_ref / rho) * dP_hat
rhs = rho_ref / dt * div_u_star + np.gradient(correction, dz)
 
# 定係数Laplace解 (FFT、周期境界仮定)
k = 2*np.pi * np.fft.fftfreq(Nz, dz)
laplacian_eigenvalues = -(k**2)
laplacian_eigenvalues[0] = 1.0  # null空間除去
P_new = np.real(ifft(fft(rhs) / laplacian_eigenvalues))
P_new -= P_new.mean()
 
# 補正項のサイズ確認
correction_magnitude = np.max(np.abs(1 - rho_ref / rho))
print(f"max |1 - rho_0/rho|  = {correction_magnitude:.4f}")
print(f"P_new range           = [{P_new.min():.3e}, {P_new.max():.3e}]")

出力は次です。

max |1 - rho_0/rho|  = 0.0117
P_new range           = [-1.572e-06, 1.572e-06]

ρ_min / ρ_ref: 0.978ρ_max / ρ_ref: 1.000x-range:

The splitting scheme rewrites the variable Poisson coefficient as constant ρ₀ plus a residual proportional to |1 − ρ₀/ρ|. For water at ΔT = 50 K the correction stays under 3% — that is why a constant-coefficient FFT solver can replace an iterative variable solver without losing accuracy.

既定のスライダーは水の 50 K です。補正は 1% の内側です。二流体モードに切り替えると、同じ項が 0.5 に近づきます。コードの ΔT = 30 K の出力 0.0117 は、その 50 K の窓より狭いです。

Gray–Giorgini(1976)の水の限界は 1.25 K です。コードの ΔT = 30 K は $30/1.25 = 24$ 倍です。空気の 10% 線は 28.6 K なので、30 K はその線もわずかに越えます。PV/T の 50 K は水の限界の 40 倍です。その 50 K で密度は約 2%、粘性は約 30% です。

密度多項式の線形係数は $3.736\times 10^{-4}$ です。30 K を掛けると $3.736\times 10^{-4}\times 30 = 0.0112$ です。出力の $\max|1-\rho_0/\rho| = 0.0117$ と並びます。二次の係数 $3.98\times 10^{-6}$ に $(30/2)^2 = 225$ を掛けると約 $9.0\times 10^{-4}$ です。0.0117 の大部分は線形項です。 $T_0 = 313\,\mathrm{K}$ 、 $dt = 10^{-3}$ 、発散の振幅 0.01 で、 $P$ の範囲は $\pm 1.572\times 10^{-6}$ です。

二流体の例は $\rho_{\min} = 1$ 、 $\rho_{\max} = 1000$ です。 $1 - 1/1000 = 0.999$ です。水の 0.0117 とは桁が二つ違います。圧力求解が一ステップの 60–80% のとき、この分割が FFT を残すと、論文の 2.5–5 倍に会います。 $\mathrm{CFL}_{\max} = 0.2$ だけを見ると、 $\mathrm{VSL} = \nu\Delta t/\Delta x^2$ が先に 0.03 を越えます。0.3 までは理論上安定でも、高次統計は 3.5% ずれます。

ΔT = 30 K の水で補正が 1.2% に収まるという文と、出力 0.0117 は同じ数字を指します。同じコードを $\rho_{\min} = 1$ 、 $\rho_{\max} = 1000$ で回すと、その値は 0.999 まで大きくなります。二流体ソルバが同じ分割で小さな $\Delta t$ を要求する理由です。

格子は $N_z = 64$ 、 $L_z = 1$ なので $dz = 1/64$ です。予測圧力は $\widehat P = 2P^n - P^{n-1}$ で、コードの $P^n$ の振幅は 0.01、 $P^{n-1}$ は 0 です。だから $\widehat P$ の振幅は 0.02 です。出力 0.0117 はその勾配に掛かる係数です。陽的補正は予測勾配の 1.17% です。Table 2 の $Ra$ は $10^6$ から $10^8$ まで二桁です。その区間で VSL の上限 0.03 を 0.3 まで緩めると、高次統計が 3.5% ずれます。

検証時に出会う三つの落とし穴

論文の Table 2(空気 DHC、 $Ra=10^6\!-\!10^8$ )を再現すると、三つの場所で足を取られます。

$\hat\rho_0$ の選択 — Dodd–Ferrante の原形は最小密度を使います。 $(1-\hat\rho_0/\hat\rho)$ を 0 以上 1 未満に閉じるためです。NOB で基準温度の密度を使っても結果は近いと論文は書いていますが、大きな ΔT ではその差が安定余裕を削ります。
CFL と VSL の両方 — 陽的 Adams–Bashforth では $\mathrm{CFL}_{\max} = 0.2$ だけでは粘性限界が先に破れます。 $\mathrm{VSL} = \nu\Delta t/\Delta x^2 < 0.03$ が実用上の上限です。0.3 までは理論上安定でも、高次統計は 3.5% ずれます。
Arakawa 型の対流項 — 運動量と運動エネルギーを離散レベルで同時に保存する必要があります。単純な中心差分はエネルギーを漏らし、大きな Ra で時間平均が静かに歪みます。

ρ₀を最小密度に固定するか

Dodd–Ferrante は最小密度を $\hat\rho_0$ に取ります。 $(1-\hat\rho_0/\hat\rho)$ を 0 以上 1 未満に置くためです。論文は NOB で基準温度の密度を使っても結果が近いと書いています。ΔT が 30 K を越えると、その差が余裕を削ります。出力 0.0117 は最小密度を使った 1D スラブの数字です。2D DNS の請求書ではありません。

2.5–5 倍は、圧力時間が 60–80% のときの話です。FFT が食べる周期方向があるときだけ、その倍率が出ます。境界層の非対称を密度 2% で説明すると、粘性 30% を落とします。その仮定が閉じる温度は Rayleigh 数 1708 の自然対流の OB の窓とは別に、水では 1.25 K で既に終わります。