난류 점성 모델 — Boussinesq, 혼합 길이, k-ε, k-ω

5장에서 닫히지 않았던 레이놀즈 응력을 “거대한 점성”으로 모델링하는 공학적 트릭. 혼합 길이부터 산업 표준 두 방정식 모델까지의 계보를 정리한다.

들어가며

5장에서 우리는 RANS 방정식을 유도하다가 막다른 길에 도달했다 — 레이놀즈 응력 $-\overline{u'_i u'_j}$ 이라는 새로운 미지수 여섯 개가 등장했지만, 식은 그대로 네 개뿐이었다(연속방정식 1 + 운동량방정식 3). 이것이 RANS의 닫힘 문제(closure problem) 다. 이 장에서는 한 세기 동안 공학자들이 사용해 온 가장 단순한 해법, “점성처럼 다루자”를 따라간다. 이 장을 마치면 OpenFOAM 케이스 파일에서 kEpsilon, kOmegaSST 같은 이름을 보았을 때 그것이 어떤 가정을 하고 어떤 한계를 갖는지 한 문장으로 답할 수 있다.

본론 1 — Boussinesq 가설: 난류를 점성처럼 다루기

부시네스크(Boussinesq, 부시네스크)는 1877년, 레이놀즈 응력이 평균 변형률과 비례한다라는 단순한 가정을 제안했다. 분자 점성에서 점성응력이 변형률에 비례하듯이( $\tau = \mu \, \partial u / \partial y$ ), 난류에서도 같은 형태의 관계가 성립한다고 가정하는 것이다:

-\overline{u'_i u'_j} = \nu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) - \frac{2}{3} k \delta_{ij}

여기서 $\nu_t$ (nu-t, 와점성 또는 난류 운동학적 점성)는 새로 도입된 계수이고, $k = \tfrac{1}{2}\overline{u'_i u'_i}$ 는 단위 질량당 난류 운동에너지(turbulent kinetic energy), $\delta_{ij}$ 는 크로네커 델타(Kronecker delta) — $i=j$ 일 때 1, 그 외 0인 기호다. 마지막 항 $\tfrac{2}{3} k \delta_{ij}$ 는 대각합(trace) 보정인데, 처음 읽을 때는 “ $i = j$ 일 때 식이 깨지지 않도록 더해 둔 항” 정도로만 이해해도 된다.

여기서 결정적인 차이점에 주의해야 한다 — 분자 점성 $\nu$ 는 유체의 물성(공기는 약 $1.5 \times 10^{-5} \text{ m}^2/\text{s}$ 로 고정)이지만, 와점성 $\nu_t$ 는 흐름의 성질이다. 같은 공기라도 정지 상태에서는 $\nu_t = 0$ , 강한 난류 흐름에서는 $\nu_t$ 가 $\nu$ 의 수백~수천 배가 될 수 있다. 즉 $\nu_t$ 는 미지수이고, 이걸 어떻게 결정하느냐가 이후 모든 모델의 출발점이다.

이 가정의 이점은 명확하다 — 레이놀즈 응력 여섯 개 미지수가 스칼라 한 개 $\nu_t$ 로 줄어든다. 이제 닫힘 문제는 “ $\nu_t$ 를 어떻게 줄 것인가”라는 한 문제로 압축된다. 한계도 명확하다 — 실제 난류에서 응력과 변형률이 정확히 비례하지 않는 흐름(곡률이 강한 흐름, 회전, 분리 등)에서는 이 가정 자체가 깨진다. 그러나 산업용 CFD의 8할 이상은 여전히 Boussinesq에 기반한다.

본론 2 — Prandtl의 혼합 길이 모델

프란틀(Prandtl, 프란틀)이 1925년 제안한 혼합 길이 모델(mixing length model) 은 가장 단순한 $\nu_t$ 결정 방식이다:

\nu_t = \ell_m^2 \left| \frac{\partial \bar{u}}{\partial y} \right|

여기서 $\ell_m$ (ell-m, 혼합 길이)은 “유체 덩어리가 다른 덩어리와 섞이기 전에 평균적으로 이동하는 거리”라는 직관적 의미를 갖는다. 기체 분자가 충돌 전까지 평균적으로 이동하는 자유 행로(mean free path)에서 영감을 얻은 개념이다.

직관은 다음과 같다. 어떤 위치 $y$ 에서 유체 덩어리가 거리 $\ell_m$ 만큼 위로 이동했다면, 그 덩어리는 원래 있던 곳의 속도를 그대로 갖고 와 주변과 충돌한다. 두 위치의 평균 속도 차이는 약 $\ell_m \cdot \partial \bar{u}/\partial y$ 이고, 이것이 곧 속도 변동 $u'$ 의 전형적 크기다. 이 추론을 응력 정의에 대입하면 위 공식이 나온다.

핵심은 $\ell_m$ 이 위치에 따라 변한다는 점이다. 벽면 근처에서는 유체 덩어리가 벽 너머로 갈 수 없으므로 $\ell_m \to 0$ 이고, 따라서 $\nu_t \to 0$ 이다. 자유 흐름에서는 $\ell_m$ 이 흐름 폭에 비례한다. 벽면 근처의 가장 간단한 가정은 카르만(Kármán, 카르만) 상수 $\kappa \approx 0.41$ 을 써서:

\ell_m(y) = \kappa \cdot \min(y, 2h - y)

이 식은 양쪽 벽에서 모두 거리 0이고, 채널 중앙에서 최대가 된다. 혼합 길이 모델은 단순하고, 미지수가 없고, 벽 근처 흐름을 잘 맞추기 때문에 지금도 경계층 코드와 교과서 예제에서 살아있다. 한계는 — 흐름마다 $\ell_m$ 을 사람이 미리 정해 줘야 한다는 것. 분리, 재부착, 충돌 제트처럼 기하가 복잡한 흐름에서는 어떤 $\ell_m$ 이 옳은지 미리 알 수 없다.

본론 3 — k-ε 두 방정식 모델

혼합 길이의 한계를 풀려면 $\nu_t$ 를 사람이 주는 대신 방정식이 풀게 만들어야 한다. k-ε 모델(Launder & Spalding, 1974)은 두 개의 추가 수송 방정식을 도입한다:

$k$ — 난류 운동에너지, 단위 질량당 변동 운동에너지 $\tfrac{1}{2} \overline{u'_i u'_i}$ . 단위는 m²/s²이다. 난류의 “세기” 를 나타낸다.
$\varepsilon$ (epsilon, 입실론) — 난류 소산율(dissipation rate). 단위 시간당 단위 질량에서 점성에 의해 열로 흩어지는 운동에너지로, 단위는 m²/s³이다. 난류가 작은 와로 부서져 사라지는 속도 를 나타낸다.

두 양으로 차원 해석을 하면 점성과 같은 차원 m²/s를 갖는 조합은 $k^2/\varepsilon$ 하나뿐이고, 여기에 보정 상수 $C_\mu \approx 0.09$ 를 붙인 형태가 정의식이다:

\nu_t = C_\mu \frac{k^2}{\varepsilon}

$k$ 와 $\varepsilon$ 각각은 자체 수송 방정식(생성 = 평균 흐름이 난류를 만들어 내는 양, 소산 = $\varepsilon$ 자체, 확산 = $\nu_t$ 로 퍼지는 양)을 풀어서 매 격자점마다 값을 얻는다. 본 책에서는 그 방정식을 유도하지 않고, “이 두 양을 풀면 $\nu_t$ 가 자동으로 정해진다”는 사실만 짚고 간다.

k-ε는 자유 흐름(제트, 후류, 대기 경계층 등)에서 안정적이고 정확하다. 반면 벽면 근처의 점성 부층(viscous sublayer) 에서는 가정이 깨져서 보정을 위한 “벽 함수(wall function)“가 필요하다. 산업적으로는 1990년대까지 가장 널리 쓰였고 지금도 외부 공기역학, 환기, 연소 등에서 표준이다.

본론 4 — k-ω 모델과 SST

윌콕스(Wilcox, 윌콕스)의 k-ω 모델(1988)은 두 번째 방정식을 $\varepsilon$ 대신 비소산율 $\omega$ (omega, 오메가)에 대해 푼다. $\omega \approx \varepsilon / (C_\mu k)$ 의 관계로 정의되며 단위는 1/s, 즉 주파수 차원이다. 와점성은:

\nu_t = \frac{k}{\omega}

이 단순한 변경이 가져다주는 가장 큰 이점은 벽 근처 거동이다. $\omega$ 방정식은 벽면에서 자연스럽게 발산하는 경계 조건을 가져서 벽 함수 없이도 점성 부층까지 적분할 수 있다. 반면 한계는 — 자유 흐름의 $\omega$ 경계값에 결과가 민감하다는 점이다. 입구 조건을 조금만 바꿔도 후류 위치가 크게 달라지는 사례가 보고되어 왔다.

멘터(Menter, 멘터)는 1994년 두 모델의 장점만 취한 SST k-ω(Shear-Stress Transport k-ω)를 제안했다. 핵심 아이디어는 단순하다 — 벽 근처에서는 k-ω로, 자유 흐름에서는 k-ε으로 자동 전환하는 블렌딩 함수를 도입한다. SST k-ω는 현재 항공우주·자동차·터보기계 산업에서 사실상 기본값이며, 12장에서 산업 CFD의 모델 선택을 다룰 때 다시 등장한다.

모델	강점	약점
혼합 길이	단순, 미지수 0개	$\ell_m$ 을 사람이 결정
k-ε	자유 흐름에서 안정	벽 근처에서 벽 함수 필요
k-ω	벽 근처 정확	자유 흐름 경계값에 민감
SST k-ω	두 영역 모두 안정	계산비용 약간 증가

파이썬으로 확인

혼합 길이 모델이 채널 단면에서 어떤 $\nu_t$ 분포를 만들어 내는지 그려 본다. 채널 반높이 $h=1$ , 평균 속도는 포물선 $\bar{u}(y) = 1 - ((y-h)/h)^2$ 로 가정한다(중앙에서 최대 1, 양쪽 벽에서 0). 그러면 $\partial \bar{u}/\partial y = -2(y-h)/h^2$ 이고, 혼합 길이는 $\ell_m(y) = 0.41 \cdot \min(y, 2h - y)$ . 결과 곡선은 양쪽 벽에서 0이고 채널 중앙 근처에서 최대가 되어, “벽 근처에서 난류 확산이 약하고 중앙에서 강하다”는 1장의 직관과 일치한다.

import numpy as np
import matplotlib.pyplot as plt

# 채널 형상: 반높이 h=1, 전체 높이 2h=2
h = 1.0
y = np.linspace(0.0, 2.0 * h, 401)

# 평균 속도 프로파일 (포물선): 중앙에서 최대, 벽에서 0
u_bar = 1.0 - ((y - h) / h) ** 2

# 속도 기울기 (해석적): du/dy = -2(y-h)/h^2
du_dy = -2.0 * (y - h) / (h ** 2)

# 혼합 길이: l_m = κ · min(y, 2h - y), κ = 0.41 (카르만 상수)
kappa = 0.41
l_m = kappa * np.minimum(y, 2.0 * h - y)

# Prandtl 혼합 길이 와점성: ν_t = l_m^2 · |du/dy|
nu_t = (l_m ** 2) * np.abs(du_dy)

# 시각화
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(y, nu_t, label=r'$\nu_t = \ell_m^2\,|\partial \bar{u}/\partial y|$')
ax.set_xlabel('y (벽으로부터의 거리)')
ax.set_ylabel(r'$\nu_t$ (와점성)')
ax.set_title('혼합 길이 모델: 채널 단면의 와점성')
ax.axvline(h, color='gray', linestyle='--', alpha=0.5, label='채널 중앙')
ax.legend()
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig('nu_t_mixing_length.png', dpi=120)
print(f"최대 ν_t = {nu_t.max():.4f} (y = {y[nu_t.argmax()]:.2f})")

곡선을 보면 $\nu_t$ 가 양쪽 벽에서 0으로 떨어지고( $\ell_m \to 0$ 때문), 중앙 근처에서 부드러운 봉우리를 이룬다. 정확히 $y = h$ 에서는 $\partial \bar{u}/\partial y = 0$ 이라 작은 골이 패고, 봉우리는 그 양옆에 두 개로 나뉜다. 이 작은 디테일이 바로 혼합 길이 모델이 “대칭축에서 $\nu_t$ 가 한 점에서 0이 되는 비물리적 거동을 보인다”라는 잘 알려진 한계의 시각적 증거다. k-ε 같은 두 방정식 모델은 이 결함을 자동으로 보정한다.

다음 장으로

지금까지는 흐름 전체에 같은 모델을 균일하게 적용했지만, 실제 흐름의 가장 흥미로운 부분은 벽 근처에 집중되어 있다. 7장: 경계층 이론에서는 점성 부층, 버퍼층, 로그 법칙층이라는 세 부분의 구조를 정리하고, 6장의 모델들이 각 층에서 어떻게(혹은 왜 잘 맞지 않는지) 동작하는지 본다. 7장이 끝나면 CFD 결과를 볼 때 “벽 근처가 맞는가”를 직접 판정할 수 있다.