파이썬으로 구현하는 초음속 노즐 CFD 시뮬레이터

들어가며

로켓 엔진이나 초음속 항공기의 추진 노즐을 설계할 때 핵심 과제는 하나다. 연소실에서 발생한 고온 고압의 가스를 최대한 효율적으로 가속시켜 추력을 극대화하는 것이다. 이 글에서는 MATLAB으로 작성된 초음속 노즐 CFD 코드를 파이썬으로 재구현하는 전 과정을 소개한다. 노즐 형상 설계부터 격자 생성, 그리고 MacCormack 유한체적법을 이용한 2차원 오일러 방정식 수치 해석까지, 각 단계의 이론적 배경과 구현 방법을 함께 설명한다.

 

2-D Nozzle Design

 

---

1. 노즐 형상 설계 — 특성곡선법

노즐 설계의 출발점은 원하는 출구 마하수를 달성하는 최소 길이 노즐(Minimum Length Nozzle)의 형상을 구하는 것이다. 이를 위해 특성곡선법(Method of Characteristics, MOC)을 사용한다. 특성곡선법은 초음속 영역에서 편미분 방정식을 상미분 방정식으로 변환해 유선을 따라 적분하는 방법으로, 2차원 초음속 유동에서 파동 전파 특성을 정확하게 반영한다.

 

구체적으로는 Prandtl-Meyer 팽창 이론을 기반으로, 목(throat)에서 출발하는 특성선들이 대칭축과 노즐 벽면 사이를 반사하며 교차하는 점들의 좌표를 순차적으로 계산한다. 각 교점에서는 리만 불변량(Riemann invariants) $K^+ = \theta + \nu$, $K^- = \theta - \nu$가 보존되며, 이를 통해 유동각 $\theta$와 Prandtl-Meyer 함수 $\nu$를 결정한다. 최종적으로 각 특성선을 벽면까지 연장함으로써 노즐 벽면의 좌표를 얻는다.

 

파이썬 구현에서는 nozzle_shape.py의 generate_bell_moc() 함수가 이 역할을 담당한다. Newton-Raphson 반복법으로 Prandtl-Meyer 함수를 역산하여 각 점의 마하수를 구하고, 2×2 선형 연립방정식의 해를 구해 특성선의 교점 좌표를 계산한다. 특성선 수 num을 늘릴수록 설계 정밀도가 높아지며, 기본값 15를 사용할 때 16개의 벽면 좌표점이 생성된다.

---

2. 격자 생성 — 구조 격자

노즐 형상이 결정되면 CFD 해석을 위한 계산 격자를 만들어야 한다. noz_mesh.py의 build_mesh() 함수는 체형(body-fitted) 구조 격자를 생성한다.

 

특성곡선법으로 얻은 벽면 좌표는 불균일한 간격의 16개 점이기 때문에, 먼저 scipy의 CubicSpline으로 300점으로 재샘플링하여 연속적인 벽면 곡선을 만든다. 이후 x 방향은 균일한 간격으로 분할하고, 각 x 위치에서 y 방향은 대칭축(y=0)부터 벽면까지 선형으로 배분한다. 격자 해상도는 mesh_factor 파라미터로 제어하며, 0.05(Coarse)부터 0.20(Fine)까지 설정할 수 있다. Medium 해상도(mesh_factor=0.10)에서는 약 91×72의 격자가 생성된다.

---

3. 지배 방정식 — 2차원 오일러 방정식

이 시뮬레이터는 점성 효과를 무시한 2차원 오일러 방정식을 지배 방정식으로 사용한다. 보존형으로 표현하면 다음과 같다.

$$\frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} = 0$$

여기서 보존 변수 벡터 $\mathbf{Q}$, x 방향 유속 벡터 $\mathbf{F}$, y 방향 유속 벡터 $\mathbf{G}$는 각각 다음과 같다.

$$\mathbf{Q} = \begin{pmatrix} \rho \ \rho u \ \rho v \ e \end{pmatrix}, \quad \mathbf{F} = \begin{pmatrix} \rho u \ \rho u^2 + p \ \rho u v \ (e+p)u \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} \rho v \ \rho u v \ \rho v^2 + p \ (e+p)v \end{pmatrix}$$

완전 기체 상태 방정식 $p = (\gamma - 1)\left(e - \frac{1}{2}\rho(u^2+v^2)\right)$으로 압력을 구하며, 비열비 $\gamma = 1.25$는 연소 가스의 특성을 반영한다.

---

4. 수치 기법 — MacCormack 유한체적법

오일러 방정식의 시간 적분에는 MacCormack 기법을 사용한다. 이 방법은 2단계 예측-수정(predictor-corrector) 구조로, 공간과 시간 방향 모두에서 2차 정확도를 가진다.

 

예측 단계(Predictor, dd=-1)에서는 전진 차분으로 임시 해 $\bar{\mathbf{Q}}$를 계산한다.

$$\bar{\mathbf{Q}}^n = \mathbf{Q}^n - \Delta t \cdot \mathcal{L}^n(\mathbf{Q}^n)$$

 

수정 단계(Corrector, dd=0)에서는 후진 차분을 적용해 최종 해를 얻는다.

$$\mathbf{Q}^{n+1} = \frac{1}{2}\left(\mathbf{Q}^n + \bar{\mathbf{Q}}^n - \Delta t \cdot \mathcal{L}(\bar{\mathbf{Q}}^n)\right)$$

 

유한체적법을 적용하므로, 각 셀의 플럭스는 네 면(right, left, top, bottom)에서의 기여를 합산해 계산한다. 각 면의 법선 벡터는 격자 좌표로부터 기하학적으로 계산하며, 이를 통해 비직교 격자에서도 정확한 플럭스 계산이 가능하다.

 

파이썬 구현에서 가장 중요한 부분은 MATLAB(1-based)과 파이썬(0-based) 인덱스 변환이다. MATLAB의 반복문 변수 $I=2\ldots n_x$(1-based)는 파이썬에서 $i=1\ldots n_x-1$(0-based)로 대응되며, flux 벡터 참조 인덱스는 아래 표와 같이 변환된다.

MATLAB (1-based)	Python (0-based)	dd=-1	dd=0
F(I+1+dd, J)	F[i+dd, j]	F[i-1,j]	F[i+1,j]
F(I+dd, J)	F[i+dd-1, j]	F[i, j]	F[i, j]
F(I, J+1+dd)	F[i, j+dd]	F[i, j-1]	F[i, j+1]
F(I, J+dd)	F[i, j+dd-1]	F[i, j]	F[i, j]

 

초기 구현은 MATLAB 원문에 충실한 이중 for 루프로 작성했으나, NumPy 슬라이싱으로 벡터화한 후 약 35배의 속도 향상을 달성했다. 예를 들어 예측 단계의 플럭스 계산은 다음과 같이 간결해진다.

# dd=-1 (forward predictor) — 벡터화
Fc = F[1:nx, 1:ny, :]    # F[i, j]   — right/top face
Fl = F[0:nx-1, 1:ny, :]  # F[i-1, j] — left face
Fb = F[1:nx, 0:ny-1, :]  # F[i, j-1] — bottom face

div = (Fc*b(sfpx) + Gc*b(sfpy) +   # right
       Fl*b(sfmx) + Gl*b(sfmy) +   # left
       Fc*b(sgpx) + Gc*b(sgpy) +   # top  (Fc = F[i,j])
       Fb*b(sgmx) + Gb*b(sgmy))    # bottom

---

5. 경계조건

올바른 경계조건 설정은 수치 안정성에 직결된다. 각 경계에 적용되는 조건은 다음과 같다.

• 입구(inflow, i=0): 노즐 목에서 음속 유입 조건을 적용한다. 축 방향 속도를 음속 $u = \sqrt{\gamma R T_c}$으로 고정하고, 횡 방향 속도는 인접 내부 셀에서 외삽한다.
• 출구(outflow, i=nx): 완전 초음속 출구이므로 모든 유동 변수를 내부 셀에서 0차 외삽한다.
• 대칭축(symmetry, j=0): 압력, 밀도, 축 방향 속도는 내부 셀과 동일하게, 횡 방향 속도는 부호를 반전시켜 반사 경계를 구현한다.
• 벽면(wall, j=ny): 오일러 방정식의 slip wall 조건으로, 압력과 밀도는 내부 셀에서 외삽하고 양쪽 속도 성분을 0으로 설정한다.

---

6. 시간 적분과 안정성

시간 스텝 $\Delta t$는 CFL(Courant-Friedrichs-Lewy) 조건에 따라 자동으로 결정된다.

$$\Delta t = \frac{\text{CFL}}{\max_{i,j}\left(\frac{a_{ij} + b_{ij} + c_{ij}}{V_{ij}}\right)}$$

여기서 $a_{ij}$, $b_{ij}$는 각 방향 속도의 기여이며, $c_{ij}$는 음속에 의한 항이다. $V_{ij}$는 각 셀의 면적으로, shoelace 공식으로 계산한다. 기본값 CFL=0.8을 사용하며, 이는 매 스텝마다 파동이 한 셀을 넘지 않도록 보장한다.

---

7. GUI 구현

사용자 인터페이스는 파이썬 기본 내장 라이브러리인 Tkinter와 Matplotlib을 조합해 구현했다. PyQt5 대신 Tkinter를 선택한 이유는 별도 설치 없이 어느 환경에서도 바로 실행할 수 있기 때문이다.

 

GUI는 왼쪽의 입력 패널과 오른쪽의 결과 탭으로 구성된다. 입력 패널은 스크롤 가능한 영역 안에 노즐 형상 선택(Cone/Bell), 기하학적 파라미터, 유동 조건, CFD 설정을 담고, 항상 보이는 하단 고정 영역에 Preview, Run CFD 버튼을 배치했다. CFD 계산은 threading.Thread를 이용해 백그라운드에서 실행하므로 계산 중에도 GUI가 응답을 유지한다.

 

결과는 5개 탭으로 표시된다. Mesh Preview는 생성된 구조 격자와 노즐 형상을 보여주고, Mach Contour는 RdYlGn 컬러맵으로 2차원 마하수 분포를, Axial Profiles는 대칭축을 따른 마하수와 압력 분포를 나타낸다. 추가된 Characteristics 탭은 Bell 노즐 설계에 사용된 특성선망과 벽면 교점을 시각화하며, Summary 탭은 노즐 기하 정보, CFD 결과, 이상 이센트로픽 예측값을 한 번에 비교할 수 있는 요약표를 제공한다.

---

8. 결과 검증

MATLAB 원문과 동일한 조건(h_th=0.025m, P_c=1.2MPa, T_c=2000K, γ=1.25, M_e=2.26)으로 계산한 결과, 마하수 범위가 약 1.2~2.4로 MATLAB의 1.0~2.45와 잘 일치한다. 입구 근방의 마하수 최솟값이 정확히 1.0이 아닌 이유는 입구 경계조건에서 속도만 음속으로 고정하고 압력과 밀도는 내부에서 외삽하는 MATLAB 원문의 방식을 그대로 따랐기 때문이다.

 

성능 측면에서는 for 루프 버전 대비 NumPy 벡터화 후 약 35배 속도가 향상되어, 91×72 격자에서 500 스텝 계산이 약 2초 만에 완료된다.

---

마치며

이 프로젝트는 MATLAB 원문의 물리적·수치적 구조를 최대한 보존하면서 파이썬의 장점을 활용하는 방향으로 진행됐다. 특성곡선법, MacCormack 기법, 유한체적법이라는 세 가지 핵심 기법이 유기적으로 연결되어 초음속 노즐 유동을 시뮬레이션한다. 향후에는 인공 점성(artificial viscosity) 추가로 수치 진동을 억제하고, 입구 경계조건을 이센트로픽 조건으로 개선하면 더욱 정확한 결과를 얻을 수 있을 것이다.