Riemann 문제: CFD의 심장

유한체적법(Finite Volume Method)에서 셀 경계의 **수치 플럭스(numerical flux)**를 어떻게 계산하느냐가 해의 정확도와 안정성을 좌우합니다. 이 수치 플럭스를 결정하는 열쇠가 바로 Riemann 문제입니다.

Riemann 문제란?

초기 조건이 하나의 불연속면을 사이에 두고 좌우로 다른 상수 상태를 갖는 초기값 문제입니다:

$$\mathbf{U}(x, 0) = \begin{cases} \mathbf{U}_L & \text{if } x < 0 \\ \mathbf{U}_R & \text{if } x > 0 \end{cases}$$

1차원 Euler 방정식에 대해 이 문제의 해석해는 세 개의 파(wave)로 구성됩니다:

1. 좌향 파 (rarefaction 또는 shock)
2. 접촉 불연속 (contact discontinuity)
3. 우향 파 (rarefaction 또는 shock)

Godunov의 아이디어 (1959)

Godunov의 핵심 통찰은 간단하지만 강력합니다:

셀 경계에서의 플럭스를 구하려면, 인접한 두 셀의 평균값을 좌우 상태로 놓고 Riemann 문제를 풀면 된다.

유한체적법의 업데이트 공식:

$$\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x}\left(\hat{\mathbf{F}}_{i+1/2} - \hat{\mathbf{F}}_{i-1/2}\right)$$

여기서 $\hat{\mathbf{F}}_{i+1/2}$가 바로 Riemann 문제의 해로부터 결정되는 수치 플럭스입니다.

Exact Riemann Solver

정확한 Riemann 풀이는 Newton-Raphson 반복법으로 star region의 압력 $p^*$를 구합니다. 좌우 파의 종류(shock/rarefaction)에 따라 Rankine-Hugoniot 관계식 또는 등엔트로피 관계식을 적용합니다.

Shock wave 관계식 ($p^* > p_K$, $K = L$ or $R$):

$$f_K(p^*) = (p^* - p_K)\left[\frac{A_K}{p^* + B_K}\right]^{1/2}$$

여기서 $A_K = \frac{2}{(\gamma+1)\rho_K}$, $B_K = \frac{\gamma-1}{\gamma+1}p_K$ 입니다.

Rarefaction wave 관계식 ($p^* \leq p_K$):

$$f_K(p^*) = \frac{2c_K}{\gamma - 1}\left[\left(\frac{p^*}{p_K}\right)^{\frac{\gamma-1}{2\gamma}} - 1\right]$$

$p^*$는 다음 조건을 만족합니다:

$$f_L(p^*) + f_R(p^*) + (u_R - u_L) = 0$$

Approximate Riemann Solver: 왜 필요한가?

Exact solver는 반복 계산이 필요하므로 비쌉니다. 실제 CFD 코드에서는 대부분 **근사 Riemann 풀이(approximate Riemann solver)**를 사용합니다.

Roe Solver (1981)

Roe는 비선형 Riemann 문제를 **선형화(linearization)**하여 행렬 $\hat{\mathbf{A}}$에 대한 선형 Riemann 문제로 바꿉니다:

$$\hat{\mathbf{F}}_{i+1/2} = \frac{1}{2}(\mathbf{F}_L + \mathbf{F}_R) - \frac{1}{2}\sum_{k=1}^{3}|\hat{\lambda}_k|\hat{\alpha}_k\hat{\mathbf{r}}_k$$

Roe 평균(Roe average)을 사용하여 $\hat{\mathbf{A}}$를 구성합니다:

$$\hat{u} = \frac{\sqrt{\rho_L}\,u_L + \sqrt{\rho_R}\,u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}, \quad \hat{H} = \frac{\sqrt{\rho_L}\,H_L + \sqrt{\rho_R}\,H_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}$$

장점: 접촉 불연속을 정확히 포착. 단점: 엔트로피 조건 위배 가능 (entropy fix 필요).

HLLC Solver (Toro, 1994)

HLL 풀이를 개선하여 접촉 불연속(Contact)을 복원한 것이 HLLC입니다. 세 개의 파속(wave speed) $S_L$, $S^*$, $S_R$을 추정합니다:

$$\hat{\mathbf{F}}_{i+1/2} = \begin{cases} \mathbf{F}_L & \text{if } S_L > 0 \\ \mathbf{F}^*_L & \text{if } S_L \leq 0 < S^* \\ \mathbf{F}^*_R & \text{if } S^* \leq 0 \leq S_R \\ \mathbf{F}_R & \text{if } S_R < 0 \end{cases}$$

star region 플럭스:

$$\mathbf{F}^*_K = \mathbf{F}_K + S_K(\mathbf{U}^*_K - \mathbf{U}_K)$$

HLLC는 구현이 간단하면서도 robust하여 압축성 다상유동 코드에서 가장 널리 사용됩니다.

AUSM+ (Liou, 1996)

AUSM 계열은 플럭스를 대류항(convective)과 압력항(pressure)으로 분리합니다:

$$\hat{\mathbf{F}}_{i+1/2} = \dot{m}_{1/2}\boldsymbol{\Psi}_{1/2} + \hat{p}_{1/2}\mathbf{D}$$

여기서 $\dot{m}_{1/2}$은 셀 경계 질량 유량, $\boldsymbol{\Psi}$는 대류되는 양입니다. 저마하 유동에서도 안정적이어서 all-speed scheme으로 확장하기 좋습니다.

다상유동으로의 확장

다상유동에서는 계면 양쪽의 EOS가 다르기 때문에, Riemann solver를 multi-material 버전으로 일반화해야 합니다.

대표적인 접근:

• Ghost Fluid Method: 계면 근처에서 상대편 유체의 ghost cell을 구성하고, 각 상에서 독립적인 single-material Riemann 문제를 품
• HLLC for multi-material: star region에서 양쪽 EOS를 각각 적용하여 $p^*$, $u^*$를 결정

다음 글에서는 **계면 포착 기법(VOF, Level Set, Diffuse Interface)**을 비교합니다.