고두노프 도식

1. 개요[편집]

고두노프 도식 Godunov's scheme
고안	S. K. Godunov (1959), Mat. Sbornik 47
틀	유한체적 + 셀 경계마다 리만 문제
정확도	공간·시간 모두 1차
안정 조건	CFL ≤ 1 (특성속도 max(|u|+a) 기준)
강점	엔트로피 해로 수렴, 인공점성 없이 충격 포착
대가	수치점성 ≈ |a|Δx/2 → 접촉면이 번진다

셀 하나하나를 상수로 뭉개 놓고 나면, 격자는 저절로 수백만 개의 충격파 튜브가 된다. 그걸 전부 정직하게 풀자는 게 이 도식이다.

고두노프 도식(Godunov’s scheme)은 유한체적법에서 셀 안의 해를 셀 평균 상수로 근사하고, 그 결과 모든 셀 경계에 저절로 생기는 국소 리만 문제를 매 스텝 풀어 경계 플럭스를 얻는 1차 보존형 도식이다. 세르게이 고두노프(Sergei K. Godunov)가 1959년 논문에서 제시했고, 오늘날 압축성 전산유체역학에서 쓰이는 충격 포착 도식의 절대다수가 이 골격 위에 서 있다.¹

이 문서는 도식을 다룬다. 리만 문제 자체의 파동 구조와 솔버 계보의 개관은 리만 솔버, 셀 재구성으로 차수를 올리는 이야기는 MUSCL·WENO 도식, 왜 비선형 제한자가 불가피한지는 전변분 감소가 소스 오브 트루스다. 여기서는 “경계마다 리만 문제를 푼다”는 결정이 무엇을 사고 무엇을 팔았는가만 판다.

2. REA 사이클 — 재구성, 발전, 평균[편집]

고두노프 도식은 세 동작의 반복으로 읽는 게 가장 깔끔하다. 판 레이어와 르베크가 정리한 REA(Reconstruct–Evolve–Average) 틀이다.

1. 재구성(Reconstruct): 셀 평균 $\bar{\mathbf{U}}_i^n$ 로부터 셀 내부 해를 복원한다. 고두노프는 여기서 가장 게으른 선택을 한다 — 셀마다 상수. 그 결과 해는 계단 함수가 되고, 모든 셀 경계 $x_{i+1/2}$ 는 좌 $\bar{\mathbf{U}}_i$, 우 $\bar{\mathbf{U}}_{i+1}$ 인 리만 문제가 된다.
2. 발전(Evolve): 이 계단 초기조건을 $\Delta t$ 동안 정확히 시간 전진시킨다. 각 경계에서 부챗살처럼 퍼지는 파동들이 그 정확해다.
3. 평균(Average): 발전된 해를 다시 셀 위에서 평균해 $\bar{\mathbf{U}}_i^{n+1}$ 을 얻는다.

3단계를 보존형으로 정리하면 결국 익숙한 플럭스 갱신식이 된다.

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

여기서 수치 플럭스는 리만 해를 $x/t = 0$ 선 위에서 평가한 물리 플럭스다.

$$\mathbf{F}_{i+1/2} = \mathbf{F}\!\left(\mathbf{U}^{\mathrm{RP}}\!\left(0;\ \bar{\mathbf{U}}_i,\ \bar{\mathbf{U}}_{i+1}\right)\right)$$

경계에서 자기 위치의 값만 필요하다는 것이 핵심이다. 리만 부채 전체를 알 필요 없이 원점 방향의 한 상태만 뽑으면 된다. 시간 적분이 정확한 이유도 여기 있다 — 리만 해는 자기유사(self-similar)라 $x/t=0$ 에서의 값이 $\Delta t$ 동안 변하지 않는다.

안정 조건은 이웃 경계에서 출발한 파동이 한 스텝 안에 서로의 경계를 침범하지 않을 것이라는 요구에서 나온다. 그것이 곧 CFL 조건

$$\frac{\Delta t}{\Delta x}\max_i\left(|u_i| + a_i\right) \le 1$$

이고, 이걸 넘기면 “정확히 발전시킨다”는 2단계의 전제가 무너져 도식이 통째로 근거를 잃는다.

3. 정확 리만 해법기[편집]

오일러 방정식의 리만 문제는 정확히 풀린다. 좌우 음향파(각각 충격 또는 팽창)와 가운데 접촉면 사이에 낀 별 영역(star region)의 압력 $p_*$ 하나만 알아내면 나머지는 대수적으로 따라 나오기 때문이다. 그 $p_*$ 는 다음 비선형 대수방정식의 근이다.

$$f(p) \;=\; f_L(p;\ \mathbf{U}_L) \;+\; f_R(p;\ \mathbf{U}_R) \;+\; (u_R - u_L) \;=\; 0$$

각 항은 그쪽 파가 충격이냐 팽창이냐에 따라 갈라진다($K = L, R$).

$$f_K(p) = \begin{cases} \left(p - p_K\right)\sqrt{\dfrac{A_K}{p + B_K}}, & p > p_K \quad (\text{충격파, 랭킨-위고니오}) \\[2.2ex] \dfrac{2 a_K}{\gamma - 1}\left[\left(\dfrac{p}{p_K}\right)^{\frac{\gamma-1}{2\gamma}} - 1\right], & p \le p_K \quad (\text{팽창파, 등엔트로피 리만 불변량}) \end{cases}$$ $$A_K = \frac{2}{(\gamma+1)\rho_K}, \qquad B_K = \frac{\gamma-1}{\gamma+1}\,p_K$$

$f$ 는 $p>0$ 에서 단조 증가하고 위로 볼록하지 않은(오목한) 함수라 근이 유일하고, 뉴턴-랩슨법이 안정적으로 수렴한다. 초기 추정을 두 팽창파 근사(TRRS)나 선형화 근사(PVRS)로 잡으면 상대오차 $10^{-6}$ 까지 보통 3~5회 반복이면 끝난다. 근을 얻은 뒤 접촉면 속도는

$$u_* = \frac{1}{2}\left(u_L + u_R\right) + \frac{1}{2}\left[f_R(p_*) - f_L(p_*)\right]$$

이다. 해가 존재할 조건, 즉 진공이 생기지 않을 조건은

$$\frac{2}{\gamma-1}\left(a_L + a_R\right) > u_R - u_L$$

이며, 양쪽이 너무 세게 서로 멀어지면 이 부등식이 깨지고 가운데에 진공이 열린다. 이때는 $f$ 에 근이 없으므로 별도의 진공 분기를 코딩해 두어야 한다. 이 예외 처리를 빠뜨린 코드는 강한 팽창 문제에서 조용히 NaN을 뱉는다.

문제는 값이 아니라 비용이다. 셀 경계 하나마다 초월함수와 거듭제곱이 섞인 반복법을 도는데, 3차원 100만 셀이면 스텝당 경계가 300만 개다. 그래서 정확 해법기는 실무 솔버의 기본값이 아니라 다른 해법기를 검증할 기준해로 주로 쓰인다.²

4. 고두노프 정리 — 1차라는 천장[편집]

고두노프는 도식만 준 게 아니라 그 도식의 한계를 증명한 정리도 같은 시기에 남겼다.

단조성을 보존하는 선형 도식은 정확도가 1차를 넘을 수 없다.

즉 계수를 해와 무관하게 고정한 채로는 “진동 없음”과 “2차 이상”을 동시에 가질 수 없다. 원래 고두노프 도식이 딱 1차인 것은 게을러서가 아니라 이 천장에 붙어 있기 때문이다.

탈출구는 하나뿐이다 — 도식을 비선형으로 만들 것. 해의 국소 거동을 보고 계수를 스스로 바꾸는 장치, 즉 제한자와 비선형 가중이 그것이다. 정리의 증명과 TVD 설계 규격은 전변분 감소가, 조각선형 재구성으로 2차를 얻는 구체적 방법은 MUSCL이, 국소 극값에서도 차수를 지키려는 시도는 ENO 도식·WENO 도식이 각각 다룬다. 여기서 반복하지 않는다. 중요한 것은 그 모든 고해상도 도식이 여전히 셀 경계에서 리만 문제를 푼다는 사실이다. 재구성이 바뀌었을 뿐 고두노프의 골격은 그대로다.

5. 근사 리만 해법기[편집]

반복법을 없애고 대수식 한 번으로 플럭스를 만드는 것이 근사 리만 해법기다. 어떤 파를 살리고 어떤 파를 뭉갤지가 설계의 전부다.

해법기	명시적으로 세는 파	접촉면·전단층	견고성	성격
Rusanov (국소 락스-프리드리히)	없음 (최대 속도 하나)	심하게 번짐	매우 높음	가장 싸고 가장 무딤
HLL	좌·우 2개	번짐	높음	중간을 상수로 뭉갬
HLLE	좌·우 2개 (로 평균 기반 속도 추정)	번짐	양의 보존 보장	밀도·압력 음수 방지
HLLC	좌·접촉·우 3개	유지	높음	사실상 업계 표준
Roe	계의 모든 특성파	가장 선명	엔트로피 픽스 필요	선형화 후 파동 분해

HLL(Harten–Lax–van Leer, 1983)은 리만 부채를 좌우 두 속도 $S_L, S_R$ 로 감싸고 그 사이를 상수 하나로 대체한다.

$$\mathbf{F}^{\mathrm{HLL}} = \frac{S_R \mathbf{F}_L - S_L \mathbf{F}_R + S_L S_R\left(\mathbf{U}_R - \mathbf{U}_L\right)}{S_R - S_L}$$

(단 $S_L \le 0 \le S_R$ 인 경우이고, 부채가 한쪽으로 쏠리면 그냥 $\mathbf{F}_L$ 또는 $\mathbf{F}_R$ 이다.) 세 줄이면 구현되고 어지간해선 안 터지지만, 가운데 접촉파를 아예 세지 않으므로 물질 경계와 전단층이 뭉개진다. HLLE(Einfeldt, 1988)는 여기에 로 평균 기반의 파속 추정을 못 박아 밀도와 내부에너지가 음수가 되지 않음을 증명한 판본이다. 강한 팽창이나 저밀도 영역에서 코드가 살아남게 해 주는 안전벨트다.

HLLC(Toro–Spruce–Speares, 1994)는 잃어버린 접촉파(C = Contact)를 되살려 중간 영역을 둘로 나눈다. 견고성은 HLL 수준으로 유지하면서 접촉면과 전단층을 살리기 때문에, 오늘날 압축성 솔버의 기본값 자리를 차지했다.

Roe(1981)는 접근이 다르다. 비선형 계를 좌우 상태에 의존하는 상수 행렬 $\tilde{\mathbf{A}}(\mathbf{U}_L,\mathbf{U}_R)$ 로 선형화하되, 그 행렬이 정확한 플럭스 차를 재현($\mathbf{F}_R - \mathbf{F}_L = \tilde{\mathbf{A}}(\mathbf{U}_R - \mathbf{U}_L)$)하고 실고유값·완전 고유벡터를 갖도록 요구한다. 오일러 방정식에서 이를 만족하는 것이 그 유명한 로 평균이다.

$$\tilde{u} = \frac{\sqrt{\rho_L}\,u_L + \sqrt{\rho_R}\,u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}, \qquad \tilde{a}^2 = (\gamma - 1)\left(\tilde{H} - \tfrac{1}{2}\tilde{u}^2\right)$$

($\tilde{H}$ 도 같은 $\sqrt{\rho}$ 가중 평균이다.) 플럭스는 파동별 기여의 합으로 쓰인다.

$$\mathbf{F}^{\mathrm{Roe}} = \frac{1}{2}\left(\mathbf{F}_L + \mathbf{F}_R\right) - \frac{1}{2}\sum_k \left|\tilde{\lambda}_k\right| \tilde{\alpha}_k \tilde{\mathbf{K}}_k$$

밀도 제곱근 가중이라는 기묘한 평균이 나오는 이유는 순전히 “플럭스 차를 정확히 재현하라”는 조건을 풀었더니 그렇게 나왔기 때문이다.

6. 엔트로피 조건과 팽창 충격[편집]

정확 해법기를 쓴 고두노프 도식은 스칼라 보존법칙에서 단조 도식이고, 단조 도식은 엔트로피 해로 수렴한다는 정리가 있다. 인공점성을 손으로 얹지 않아도 물리적으로 옳은 불연속만 골라 잡는 것이 이 도식의 최대 자산이다.

그런데 근사 해법기는 이 보증을 함께 사 오지 못한다. 위 Roe 플럭스에서 수치 소산은 $|\tilde{\lambda}_k|$ 에 실려 있는데, 음속점($\tilde{\lambda}_k = 0$, 즉 팽창파가 $u = \pm a$ 를 가로지르는 자리)에서 이 계수가 0이 된다. 소산이 사라진 자리에서 도식은 정상 팽창파 대신 팽창 충격(expansion shock)이라는 불연속을 붙잡아 버린다. 엔트로피가 감소하는 해라 물리적으로 존재할 수 없는데도 이산 방정식은 만족시키기 때문에, 잔차는 얌전히 떨어지고 결과만 틀린다.

표준 처방이 하튼의 엔트로피 픽스다. 작은 고유값을 인위적으로 바닥에서 들어 올린다.

$$|\lambda| \;\longrightarrow\; \begin{cases} |\lambda|, & |\lambda| \ge \delta \\[1ex] \dfrac{\lambda^2 + \delta^2}{2\delta}, & |\lambda| < \delta \end{cases}$$

$\delta$ 는 국소 파속에 비례하는 작은 값으로 잡는다. 이 한 줄이 없으면 라발 노즐의 목처럼 음속점이 격자에 딱 걸리는 문제에서 해가 비물리적으로 갈라진다. 반대로 $\delta$ 를 너무 키우면 접촉면 해상도가 HLL 수준으로 주저앉는다 — 이 저울질에 걸린 마법 상수 하나가 여전히 실무의 튜닝 손잡이로 남아 있다.

7. 수치점성 — 무엇이 얼마나 번지나[편집]

선형 이류 $u_t + a u_x = 0$ 에 고두노프 도식을 적용하면 정확히 1차 상류 차분이 된다. 그 수정 방정식(modified equation)은

$$u_t + a u_x = \underbrace{\frac{|a|\,\Delta x}{2}\left(1 - |\nu|\right)}_{\text{수치점성 } \varepsilon}\,u_{xx} + \mathcal{O}(\Delta x^2), \qquad \nu = \frac{a\,\Delta t}{\Delta x}$$

이다. 즉 도식이 주문하지 않은 확산항을 공짜로 끼워 판다. 계수의 규모가 $|a|\Delta x/2$ 라는 게 요점이고, 괄호 안의 $(1-|\nu|)$ 때문에 쿠랑수를 1에 가깝게 밀수록 소산이 줄어든다(스칼라 1차원에서 $\nu = 1$ 이면 정확해다). 실무에서 CFL을 0.2가 아니라 0.8~0.9로 잡는 습관에는 이 이유가 절반쯤 깔려 있다.³

결과의 비대칭이 재미있다.

• 충격파는 안 번진다. 비선형 압축이 특성선을 스스로 모아 수치확산과 균형을 이루므로, 1차 고두노프도 충격을 보통 2~3셀 폭에 붙잡아 둔다.
• 접촉 불연속과 전단층은 계속 번진다. 여기엔 압축 메커니즘이 없어서 순수 확산이 이긴다. 폭은 $\sqrt{\nu(1-\nu)\,n}$ 꼴, 즉 스텝 수의 제곱근으로 자란다. 1만 스텝이면 100셀 규모로 퍼진다는 뜻이다.
• 매끈한 극값도 깎인다. 사인파를 오래 이류시키면 진폭이 지수적으로 줄어든다.

다상 계면이나 화염면, 물질 경계를 오래 끌고 가야 하는 계산에서 1차 고두노프가 못 쓰이는 이유가 두 번째 항목이다. 자세한 위상·진폭 오차 분해는 수치 소산과 분산 참고.

8. 카벙클 — 너무 잘 듣는 해법기의 병[편집]

접촉면을 살리는 해법기(정확·Roe·HLLC)를 격자선과 나란한 강한 충격파에 쓰면, 정상 충격 앞면이 국소적으로 튀어나오며 종기 같은 돌기가 자란다. 1988년 피어리와 임레이가 이름 붙인 카벙클 현상(carbuncle phenomenon)이다. 극초음속 뭉툭한 물체의 이탈 충격파 정체선에서 가장 잘 재현되며, 결과적으로 정체점 압력과 공력가열이 통째로 틀린다.

퀴크(1994)의 분류가 유명하다. 접촉면을 잘 해상하는 해법기일수록 카벙클에 취약하고, 접촉·전단파를 뭉개는 해법기(HLL, HLLE, Rusanov)는 면역이다. 원인은 충격면 횡방향(전단 방향) 수치소산의 부족이다. 충격 내부의 몇 셀이 물리적으로 정의되지 않은 중간 상태를 갖는데, 그 상태의 횡방향 섭동을 눌러 줄 소산이 없으면 섭동이 자라 충격면을 왜곡한다. 격자 종횡비가 클수록, 충격이 격자선과 정확히 나란할수록 잘 터진다.

처방은 전부 “필요한 곳에만 소산을 더한다”는 형태다.

• 혼합 스위치: 충격 감지기로 강한 충격 근처에서만 HLL/Rusanov로 갈아탄다.
• H-보정(Quirk 1994): 이웃 경계의 파속을 참조해 횡방향 소산 하한을 준다.
• 회전 리만 해법기: 플럭스를 충격 법선·접선 방향으로 분해해 법선 쪽엔 정확한 해법기를, 접선 쪽엔 소산 큰 해법기를 쓴다.
• 격자 정렬 회피: 충격면과 격자선이 나란해지지 않게 만든다. 적응 격자 세분화로 충격을 잡을 때 오히려 정렬이 심해지는 역설이 있으니 주의.

완전한 해결책은 아직 없다. “견고함과 정확도는 함께 살 수 없다”는 이 바닥의 오래된 저울질이 가장 노골적으로 드러나는 자리다.⁴

9. 확장과 사촌들[편집]

고두노프 골격은 오일러 방정식 밖으로도 잘 퍼졌다. 다차원은 방향 분할이나 코렐라의 CTU(corner transport upwind)로, 자기유체역학(MHD)은 $\nabla\cdot\mathbf{B}=0$ 제약을 붙인 8파 리만 해법기로, 다상유동은 계면을 가로지르는 리만 문제로 각각 확장됐다. 아예 리만 문제를 안 푸는 계보도 있다 — 플럭스를 상·하류 성분으로 쪼개는 플럭스 벡터 분할(Steger–Warming, van Leer)과, 대류항과 압력항을 분리해 다루는 AUSM 도식(Liou)이 대표적이며, 저마하수 거동과 견고함에서 좋은 평판을 얻었다.

10. 관련 문서[편집]

• 리만 솔버 · 유한체적법
• MUSCL · ENO 도식 · WENO 도식 · 전변분 감소
• 압축성 유동 · 충격파 · 오일러 방정식 (비점성 유동)
• CFL 조건 · 수치 소산과 분산 · 차분 도식
• 초음속 유동 · 극초음속 유동 · 라발 노즐
• 뉴턴-랩슨법 · 적응 격자 세분화
• 전산유체역학 · 검증 및 확인

11. Footnotes[편집]

1. 고두노프(1929~2023)는 이 도식을 소련 핵무기 개발과 얽힌 계산 업무 중에 만들었다. 본인 회고에 따르면 정작 “선형 단조 도식은 1차가 한계”라는 정리 쪽이 더 마음에 들었다는데, 세상은 정리보다 도식에 그의 이름을 붙였다. 부정적 결과보다 쓸 수 있는 물건이 이기는 건 언제나 그렇다. ↩

2. 그래서 정확 해법기의 진짜 용처는 검증 및 확인이다. Sod 충격파 튜브 같은 표준 문제의 해석해를 뽑아 근사 해법기의 오차를 재는 자, 즉 눈금자 역할. 프로덕션 코드에서 정확 해법기를 켜 놓고 돌리는 사람은 대개 뭔가 잘못됐거나, 논문 마감이 아직 멀었거나 둘 중 하나다. ↩

3. 물론 CFL을 1에 붙이면 좋다는 건 1차원 스칼라 얘기다. 다차원 계에서는 파속 추정이 보수적이지 않으면 CFL 0.95에서 조용히 터지고, 그때 초보자는 “CFL을 낮췄더니 안정해졌다”며 0.2로 내려 버린다. 안정해진 게 아니라 소산으로 덮은 것이다. ↩

4. 카벙클을 처음 본 사람은 대개 격자를 의심하고, 그다음엔 경계조건을 의심하고, 마지막에야 플럭스 함수를 의심한다. 순서가 거꾸로인 게 문제인데 — 격자를 촘촘히 할수록 카벙클은 더 선명해지므로, 통상적인 격자 수렴성 검사가 오히려 “격자가 부족하다”는 오진을 강화한다. ↩