[논문 리뷰] 30 K가 넘으면 Boussinesq가 거짓말한다 — Demou(2018) NOB 자연대류와 압력 분할

Gray와 Giorgini는 1976년에 한 가지 숫자를 계산했다. 상온 공기에서 Oberbeck–Boussinesq 근사가 10% 오차 안에 머무는 최대 온도차는 28.6 K. 물에서는 단 1.25 K. 솔라 타워 안의 공기가 수백 K를 오가고, PV/T 시스템의 물이 50 K씩 끓어오르는 현실에서 이 가정이 얼마나 빨리 꺾이는지를 보여주는 숫자다. Demou·Frantzis·Grigoriadis(2018)는 그 가정을 버린 자연대류 솔버를 제안한다 — 모든 식에서 물성이 온도에 따라 변하고, 그런데도 압력 Poisson은 여전히 상수계수로 풀린다. 그 한 줄짜리 트릭을 따라가 본다.

한 줄 정리#

저자: 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배 빨라진 NOB DNS.

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$ 는 모두 기준 온도에서의 값으로 정규화된 무차원 물성이다. 변동성은 운동량 식뿐 아니라 에너지 식에도 같이 들어간다.

물성을 직접 보자#

물의 물성을 온도 함수로 그렸다. 슬라이더로 hot 벽과 cold 벽 온도를 바꾸면, 네 패널이 각각 얼마나 다르게 변하는지 보인다.

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% 정도밖에 안 움직인다. 그런데 점성은 같은 ΔT에서 30%가 흔들린다. 평균 밀도만 보고 OB 근사를 쓰면, 점성 변화가 만든 비대칭 경계층을 통째로 놓친다. Demou의 본문 §3.2는 물 RB 대류에서 정확히 이 비대칭이 Nusselt 수 분포를 비뚤게 만드는 과정을 보여준다.

가변계수 Poisson의 함정#

분수단계법(fractional-step)을 그대로 쓰면 압력 보정 단계에서 발산-자유 조건이 다음을 요구한다.

\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 기반 직접 솔버(FDS)는 상수계수 라플라스만 받아들이므로 못 쓰고, 결국 다중격자나 반복 솔버로 빠진다. 60–80% 시간이 압력 풀이에 들어가는 NS 솔버에서 이 차이는 결정적이다.

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 물 대류에서는 매우 작다. 아래 그래프는 그 차이가 얼마나 다른 차원인지 보여준다.

ρ_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에 근접한다 — 분할 스킴이 두-유체에서는 작은 $\Delta t$ 를 요구하지만, 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]

ΔT = 30 K의 물에서 보정 항이 1.2%로 묶인다. 같은 코드를 두-유체(예: $\rho_{\min} = 1, \rho_{\max} = 1000$ )에 돌리면 그 값이 0.999로 폭발한다. 두-유체 솔버가 같은 분할 스킴에서 작은 $\Delta t$ 를 요구하는 이유다.

검증 시 만난 함정#

논문의 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에서 시간평균 통계가 슬며시 비뚤어진다.

마지막 정리#

NOB 자연대류는 가변 점성과 가변 전도율이 만든 비대칭 경계층을 그대로 풀어야 의미가 있다. Demou·Frantzis·Grigoriadis는 그 가변성을 운동량·에너지 식에는 살려두고, 압력 Poisson에서만 분할 트릭으로 빼냈다. 결과는 한 줄짜리 치환과 한 항짜리 명시적 보정이지만, 이 한 줄이 FFT 직접 솔버를 그대로 쓸 수 있게 해주고, 그 차이가 3D DNS를 가능하게 만든다.