[논문 리뷰] 30 K가 넘으면 Boussinesq가 거짓말한다 — 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의 함정#

분수단계법(fractional-step)을 그대로 쓰면 압력 보정 단계가 발산-자유를 이렇게 요구한다. 이 보정은 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)
 
# 상수계수 라플라스 풀이 (FFT, 주기경계 가정)
k = 2*np.pi * np.fft.fftfreq(Nz, dz)
laplacian_eigenvalues = -(k**2)
laplacian_eigenvalues[0] = 1.0  # null space 제거
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$ 과 나란히 선다. 2차 계수 $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에서 이미 끝난다.