복사자기유체역학(Radiation Magnetohydrodynamics) 코드 개발 I
복사 자기유체역학(Rad-MHD) 방정식의 수치적 해법: FLD 근사와 M1 폐쇄 사용
서론
동기부여와 물리적 배경
우주의 많은 천체물리 현상은 복사, 자기장, 그리고 플라즈마가 복잡하게 얽혀 있는 환경에서 일어납니다. 몇 가지 대표적인 예를 들면:
- 강착원반: 블랙홀이나 중성자별 주위에서 물질이 강하게 끌려들어가는 과정에서 방출되는 복사가 에딩턴 한계에 가까워지면, 복사 압력이 물질의 낙하를 저해할 수 있고, 자기장은 강력한 제트를 형성하는 데 중요한 역할을 합니다.
- 별 형성 영역: 밀도가 높은 분자운이 자체 중력으로 붕괴할 때, 복사는 온도를 조절하며, 자기장은 각운동량 수송을 통해 붕괴 속도에 영향을 미칩니다.
- 초신성 잔해: 폭발한 항성의 충격파가 주변 성간 물질과 부딧히면서 강한 복사가 방출되고 자기장이 증폭됩니다.
이러한 현상들을 정확히 이해하고 예측하려면 복사, 자기장, 그리고 플라즈마 운동을 동시에 다루는 이론적 틀이 필요합니다. 이것이 바로 복사 자기유체역학(Radiation Magnetohydrodynamics, Rad-MHD)입니다.
계산적 도전과 근사법의 필요성
그런데 복사 전달 방정식을 엄밀히 풀려면 어마어마한 계산량이 필요합니다. 복사는 각 위치(3차원)에서 모든 방향(2차원)과 모든 주파수(1차원)에 대해 고려해야 합니다.. + 시간(1차원) = 7차원
이런 계산 부담을 현실적인 수준으로 줄이기 위해 다양한 근사법이 개발되었습니다:
- 플럭스 제한 확산(Flux-Limited Diffusion, FLD): 복사를 확산 현상으로 간주하여 계산을 단순화하면서도, 비물리적인 결과를 방지하기 위해 플럭스에 제한을 둡니다.
- M1 모멘트 폐쇄(M1 moment closure): 복사의 각도 분포를 저차 모멘트로 근사하여 계산 효율성을 높입니다.
이 글에서는 이 두 방법의 물리적 배경, 수학적 유도 과정, 그리고 실제 코드에서 어떻게 구현하는지 살펴보겠습니다.
Rad-MHD의 지배 방정식
Rad-MHD를 수학적으로 기술하려면 일반적인 MHD 방정식에 복사와 관련된 항들을 추가해야 합니다. 각 방정식을 차례로 살펴보겠습니다:
질량 보존 방정식: \(\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0\)
가장 기본적인 질량 보존 법칙으로, 어떤 체적 내에서 질량이 생성되거나 사라지지 않는다는 것을 나타냅니다.
운동량 보존 방정식: \(\frac{\partial (\rho \mathbf{v})}{\partial t} + \nabla \cdot \left(\rho \mathbf{v} \otimes \mathbf{v} + \left(P + \frac{B^2}{2} \right) \mathbf{I} - \mathbf{B} \otimes \mathbf{B} \right) = \mathbf{G}_{rad}\)
복잡해 보이지만, 이것은 단순히 뉴턴의 제2법칙이 MHD와 복사로 확장된 것입니다. 우변의 $\mathbf{G}$는 복사가 물질에 가하는 힘을 나타내며, 광자가 물질에 흡수되거나 방출될 때 운동량을 주고받는 것을 표현합니다.
에너지 보존 방정식: \(\frac{\partial E}{\partial t} + \nabla \cdot \left[(E + P + \frac{B^2}{2})\mathbf{v} - \mathbf{B}(\mathbf{v} \cdot \mathbf{B})\right] = cG^{0}\)
에너지 방정식에서 핵심은 우변의 소스 항입니다. 물질과 복사 사이의 에너지 교환과 복사 압력이 물질에 하는 일을 나타냅니다.
자기 유도 방정식: \(\frac{\partial \mathbf{B}}{\partial t} + \nabla \times (\mathbf{B} \times \mathbf{v}) = \nabla \times (\eta \nabla \times \mathbf{B})\)
패러데이의 유도 법칙을 나타내는 이 방정식은 플라즈마 운동에 따른 자기장의 변화를 기술합니다. 우변의 $\eta$ 항은 실제 플라즈마에서 나타나는 저항성 효과를 표현하며, 이상적인 경우에는 무시할 수 있습니다. 중요한 것은 자기장에는 단극(모노폴)이 없다는 제약조건입니다: \(\nabla \cdot \mathbf{B} = 0\)
복사 에너지 방정식: \(\frac{1}{\hat{c}}\frac{\partial E_r}{\partial t} + \nabla \cdot \mathbf{F}_r = -G^0\)
복사 에너지의 보존을 나타내는 이 방정식은 앞서 나온 에너지 방정식과 짝을 이룹니다. 여기서, $\hat{c}$ 는 reduced speed of light 입니다. 비상대론적인 경우를 다룰 때, 오버헤드(overhead)를 줄이기 위해 사용된 수치기법 입니다.
복사 운동량 방정식 (M1 폐쇄의 경우): \(\frac{1}{\hat{c}}\frac{\partial \mathbf{F}_r}{\partial t} + \nabla \cdot \mathbb{P}_r = -\mathbf{G}\)
M1 폐쇄를 사용할 때만 필요한 방정식으로, 복사 플럭스의 진화를 기술합니다. $\mathbb{P}_r$은 복사 압력 텐서로, M1 폐쇄에서는 이것을 $E_r$과 $\mathbf{F}_r$의 함수로 표현합니다.
상태 방정식: \(P = (\gamma - 1)\rho e\) \(T = \frac{\mu m_H P}{\rho k_B}\)
마지막으로, 이 방정식들은 밀도, 압력, 온도, 내부 에너지 사이의 관계를 연결해줍니다. 일반적으로 이상 기체 가정을 사용합니다.
플럭스 제한 확산(FLD) 근사법
FLD 근사의 기본 아이디어는 간단합니다. 복사를 열처럼 확산하는 현상으로 보는 것입니다. 물질이 충분히 불투명하고 복사가 모든 방향에서 거의 균등하게 들어오는 경우(광학적으로 두꺼운 영역), 이 근사는 매우 잘 작동합니다.
문제는 광학적으로 얇은 영역에서 발생합니다. 순수한 확산 방정식을 사용하면 복사가 빛의 속도보다 빠르게 퍼질 수 있다는 비물리적인 결과가 나옵니다. 이는 복사 플럭스가 물리적 한계인 $|\mathbf{F}_r| \le cE_r$를 초과할 수 있기 때문입니다.
이 문제를 해결하기 위해 플럭스 제한자(flux limiter)라는 영리한 장치를 도입합니다. 이것은 마치 자동차의 속도 제한 장치처럼, 복사 플럭스가 물리적 한계를 넘지 않도록 조절합니다. FLD에서 복사 플럭스는 다음과 같이 표현됩니다:
\[\mathbf{F}_r = -D_{rad} \nabla E_r = -\frac{c\lambda(R)}{\kappa_R\rho}\nabla E_r\]여기서 $\kappa_R$은 로셀란드 평균 불투명도로, 주파수에 대한 특별한 평균 방법입니다. $\lambda(R)$는 플럭스 제한자로, 다음과 같이 정의되는 무차원 매개변수 $R$의 함수입니다:
\[R = \frac{|\nabla E_r|}{\kappa_R\rho E_r}\]$R$은 복사 에너지가 얼마나 급격하게 변하는지를 나타내는 지표입니다. 이 매개변수에 따라 플럭스 제한자의 동작이 결정됩니다:
- 광학적으로 두꺼운 영역($R \ll 1$): $\lambda(R) \to 1/3$로 접근하여 표준 확산 근사를 회복합니다.
- 광학적으로 얇은 영역($R \gg 1$): $\lambda(R) \approx 1/R$로 감소하여 플럭스가 $cE_r$을 초과하지 않도록 합니다.
대표적인 플럭스 제한자
여러 플럭스 제한자가 제안되었는데, 가장 널리 사용되는 것들은:
레버모어-폼레이닝(Levermore-Pomraning) 제한자: \(\lambda(R) = \frac{2+R}{6+3R+R^2}\)
미너보(Minerbo) 제한자: \(\lambda(R) = \frac{2}{3+\sqrt{9+12R^2}}\)
두 제한자 모두 필요한 극한 조건을 만족하며, 실제 사용에서는 큰 차이가 없습니다.
FLD 근사를 복사 에너지 방정식에 대입하면 다음과 같은 비선형 확산 방정식을 얻습니다:
\[\frac{\partial E_r}{\partial t} - \nabla \cdot \left(\frac{c\lambda(R)}{\kappa_R\rho}\nabla E_r\right) = c\kappa_E\rho(aT^4 - E_r)\]이 방정식은 확산 계수가 $\lambda(R)$에 의해 비선형적으로 변하기 때문에 특별한 수치적 처리가 필요합니다. 또한 시간 스케일의 강성성 때문에 대부분 암시적 방법을 사용합니다.
FLD의 장단점
장점:
- 계산 비용이 낮고 구현이 비교적 간단합니다
- 광학적으로 두꺼운 영역에서 매우 정확합니다
- 안정적인 수치 해법이 잘 정립되어 있습니다
한계:
- 복사의 방향성을 제대로 표현하지 못합니다 (그림자, 빔 현상 등)
- 광학적으로 얇은 영역과 두꺼운 영역의 경계에서 부정확할 수 있습니다
- 플럭스 제한자 선택에 따라 결과가 달라질 수 있습니다
M1 폐쇄 모델
M1 폐쇄는 FLD보다 정교한 방법으로, 복사의 방향성을 어느 정도 표현할 수 있습니다. 기본 개념은 복사 전달 방정식을 각도에 대해 적분하여 모멘트 방정식을 유도하는 것입니다:
- 0차 모멘트: 복사 에너지 밀도 $E_r$ (모든 방향의 복사를 합친 것)
- 1차 모멘트: 복사 플럭스 $\mathbf{F}_r$ (알짜 흐름의 방향과 크기)
- 2차 모멘트: 복사 압력 텐서 $\mathbb{P}_r$ (복사의 방향성 정보)
문제는 이 모멘트들이 무한히 계속된다는 것입니다. 어딘가에서 끊어야 하는데, M1 폐쇄는 2차 모멘트에서 끊습니다.
M1 폐쇄의 핵심은 복사 압력 텐서 $\mathbb{P}_r$을 이미 알고 있는 $E_r$과 $\mathbf{F}_r$로 표현하는 것입니다:
\[\mathbb{P}_r = \mathbb{D}(\mathbf{f})E_r\]여기서 $\mathbf{f} = \mathbf{F}_r/(cE_r)$는 축소된 플럭스로, 복사가 얼마나 방향성을 가지는지를 나타내는 지표입니다. 물리적으로 $|\mathbf{f}| \le 1$이어야 하는데, 이는 복사 플럭스가 복사 에너지에 빛의 속도를 곱한 것보다 클 수 없다는 의미입니다.
$\mathbb{D}(\mathbf{f})$는 에딩턴 텐서로, Levermore(1984)가 제안한 형태는:
\[\mathbb{D}(\mathbf{f}) = \frac{1-\chi(f)}{2}\mathbb{I} + \frac{3\chi(f)-1}{2}\mathbf{n}\otimes\mathbf{n}\]$\mathbf{n} = \mathbf{f}/f$는 플럭스 방향의 단위 벡터이고, $\chi(f)$는 에딩턴 인자로:
\[\chi(f) = \frac{3+4f^2}{5+2\sqrt{4-3f^2}}\]이 함수는 두 가지 극한 상황을 자연스럽게 연결합니다:
- 등방성 확산 ($f \to 0$): 복사가 모든 방향에서 균등하면 $\chi \to 1/3$
- 자유 흐름 ($f \to 1$): 복사가 한 방향으로만 흐르면 $\chi \to 1$
M1 폐쇄를 사용하면 복사 에너지와 복사 플럭스에 대한 쌍곡선형 편미분 방정식을 얻게 됩니다. 이것은 FLD의 포물선형 확산 방정식과는 달리 파동 형태로 전파하는 해를 가지므로, 복사의 전파 특성을 더 잘 표현할 수 있습니다.
M1 폐쇄의 장단점
장점:
- FLD보다 복사의 방향성을 훨씬 잘 표현합니다
- 그림자나 빔 현상을 어느 정도 모사할 수 있습니다
- 확산 한계와 자유 흐름 한계를 자연스럽게 연결합니다
한계:
- FLD보다 계산 비용이 높습니다 (방정식이 더 많고 복잡함)
- 여러 개의 복사 빔이 교차하는 상황은 제대로 처리하지 못합니다
- $f > 1$이 되는 비물리적 상황을 방지하기 위한 추가 처리가 필요할 수 있습니다
FLD와 M1 선택 가이드
언제 어떤 방법을 선택해야 하는가?
적절한 방법 선택은 문제의 물리적 특성과 계산 자원에 따라 달라집니다:
FLD를 선택해야 하는 경우:
- 광학적으로 두꺼운 매질이 지배적인 경우 (예: 항성 내부, 광학적으로 두꺼운 강착원반)
- 복사장이 거의 등방성에 가까운 경우
- 계산 자원이 제한적이고 빠른 계산이 필요한 경우
- 장시간 진화를 추적해야 하는 경우
- 복사의 정확한 방향성이 중요하지 않은 경우
M1을 선택해야 하는 경우:
- 광학적으로 얇은 영역과 두꺼운 영역이 공존하는 경우
- 복사 빔이나 그림자 효과가 중요한 경우
- 복사압에 의한 유체 가속이 중요한 경우 (예: 복사 구동 바람)
- 충격파 근처의 복사 전달을 정확히 모사해야 하는 경우
- 계산 자원이 충분하고 높은 정확도가 필요한 경우
하이브리드 접근법
일부 코드에서는 영역에 따라 다른 방법을 사용하는 하이브리드 접근법도 고려됩니다:
- 광학적으로 두꺼운 영역: FLD 사용
- 광학적으로 얇은 영역: M1 또는 더 정교한 방법 사용
- 전이 영역: 두 방법을 부드럽게 연결
💡 실용적 조언: 처음에는 FLD로 시작하여 문제를 이해하고, 필요시 M1으로 전환하는 것이 효율적입니다.
수치 해법 접근법
Rad-MHD 방정식을 FLD 또는 M1 폐쇄와 함께 푸는 것은 복잡한 다중물리 문제입니다. 일반적으로 다음과 같은 수치적 요소들을 고려해야 합니다:
1. 공간 이산화 (Spatial Discretization)
유한 체적법(Finite Volume Method, FVM)이 보존 법칙을 정확히 만족시키기 때문에 널리 사용됩니다. MHD 부분에는 고차 정확도 기법(예: Piecewise Parabolic Method, PPM; Weighted Essentially Non-Oscillatory, WENO; Monotonic Upstream-centered Scheme for Conservation Laws, MUSCL-Hancock)과 함께 근사 리만 솔버(Approximate Riemann Solver, 예: HLL, HLLC, HLLD)가 주로 사용됩니다.
$\nabla \cdot \mathbf{B} = 0$ 조건을 유지하기 위해 Constrained Transport (CT) 기법이나 발산 세정(divergence cleaning) 기법이 필요합니다.
2. 시간 적분 (Time Integration)
Rad-MHD 방정식은 서로 다른 시간 스케일을 갖는 항들(예: 유체 역학적 시간, 자기장 확산 시간, 복사 전파 시간, 방사-물질 결합 시간)을 포함하므로 매우 강성적(stiff)일 수 있습니다.
- 명시적 방법 (Explicit methods): CFL 조건에 의해 시간 단계 $\Delta t$가 제한되지만, 구현이 비교적 간단합니다. MHD의 쌍곡선 항에 주로 사용됩니다.
- 암시적 방법 (Implicit methods): 시간 단계 제한이 덜 엄격하여 강성적인 항(예: 확산 항, 방사-물질 결합 항) 처리에 유리하지만, 각 시간 단계마다 비선형 연립방정식을 풀어야 합니다.
- IMEX (Implicit-Explicit) 방법: 일부 항은 명시적으로, 다른 항(주로 강성적인 항)은 암시적으로 처리하여 효율성과 안정성을 동시에 추구합니다.
3. 연산자 분할 (Operator Splitting)
Rad-MHD 방정식의 복잡성을 줄이기 위해 연산자 분할 기법이 자주 사용됩니다. 예를 들어, 한 시간 단계 $\Delta t$ 동안 MHD 부분과 복사 전달 부분을 분리하여 순차적으로 풀 수 있습니다. 스트랭 분할(Strang splitting)과 같은 2차 정확도 분할 방식이 선호됩니다:
이러한 분할은 각 물리 과정에 최적화된 수치 기법을 독립적으로 적용할 수 있게 하지만, 분할 오차(splitting error)를 유발할 수 있으므로 주의해야 합니다.
FLD를 위한 수치 방법
1. 암시적 확산 해법
FLD 근사에서 얻어지는 비선형 확산 방정식 (복사 에너지 방정식)은 일반적으로 그 강성성 때문에 암시적으로 풀어야 합니다. 시간 $n$에서 $n+1$로 전진하는 이산화된 방정식은 다음과 같은 형태를 가집니다:
\[\frac{E_r^{n+1} - E_r^n}{\Delta t} - \nabla \cdot \left(D(E_r^{n+1}) \nabla E_r^{n+1}\right) = S(T^{n+1}, E_r^{n+1})\]여기서 $D(E_r) = \frac{c\lambda(R(E_r))}{\kappa_R\rho}$는 비선형 확산 계수이고, $S = c\kappa_E\rho(aT^4 - E_r)$는 소스 항입니다. 일반적으로 국소 열역학적 평형(LTE) 조건에서는 에너지 평균, Planck 평균, 흡수 불투명도가 모두 같지만($\kappa_E = \kappa_P = \kappa_A$), 비LTE 상황에서는 다를 수 있습니다. $T^{n+1}$ 또한 $E_r^{n+1}$ (그리고 다른 유체 변수)과 결합될 수 있어 문제는 더욱 비선형적입니다.
2. 확산 방정식의 이산화 및 선형화
유한 체적법을 사용하여 셀 $i$에 대해 이산화하면, 각 면(face)을 통한 플럭스와 셀 내부의 소스 항을 포함하는 대수 방정식을 얻습니다. 확산 계수 $D$와 소스 항 $S$의 비선형성 때문에, 각 시간 단계에서 반복적인 해법이 필요합니다. 일반적으로 이전 반복 단계 $k$의 값으로 $D(E_r^k)$와 $S(T^k, E_r^k)$를 계산하고, $E_r^{k+1}$에 대해 선형 방정식을 풉니다 (Picard 반복). 또는 Newton-Krylov 방법과 같은 더 정교한 비선형 솔버를 사용할 수 있습니다.
\[(A_{diff} + A_{source}) E_r^{n+1, k+1} = \mathbf{b}\]여기서 $A_{diff}$는 확산 연산자, $A_{source}$는 소스 항 중 $E_r^{n+1}$에 비례하는 부분 (예: $-c\kappa_A\rho E_r^{n+1}$), $\mathbf{b}$는 나머지 항들을 포함합니다.
3. 행렬 해법
이산화된 방정식은 거대한 희소 행렬(sparse matrix) 시스템을 형성합니다. 이 시스템은 GMRES, BiCGSTAB과 같은 Krylov 부공간 반복법(Krylov subspace iterative methods)과 함께 적절한 예비 조건자(preconditioner, 예: AMG, ILU)를 사용하여 효율적으로 풀 수 있습니다.
M1 폐쇄를 위한 수치 방법
1. 쌍곡선 시스템 해법
M1 모델 방정식은 쌍곡선 시스템을 형성하며, 이는 적절한 리만 해법(Riemann solver)을 사용한 고두노프(Godunov) 유형 방법으로 해결할 수 있습니다:
2. 주요 단계
M1 모델의 수치 해법은 일반적으로 다음 단계들을 포함합니다:
- 재구성 (Reconstruction): 셀 평균값으로부터 셀 경계에서의 좌우 상태($W_L, W_R$)를 고차 정확도로 재구성합니다. (예: MUSCL-Hancock)
- 리만 문제 해결 (Approximate Riemann Solver): $W_L, W_R$을 초기 조건으로 하는 리만 문제를 풀어 셀 경계에서의 수치 플럭스 $\mathcal{F}_{face}$를 계산합니다.
- HLL 플럭스의 경우, 가장 빠른 신호 속도 $S_L, S_R$을 추정하여 플럭스를 계산합니다.
- $\mathcal{F}_{HLL} = \frac{S_R \mathbf{F}(\mathbf{U}_L) - S_L \mathbf{F}(\mathbf{U}_R) + S_L S_R (\mathbf{U}_R - \mathbf{U}_L)}{S_R - S_L}$
- 여기서 $\mathbf{U} = [E_r, F_{r,x}/c, F_{r,y}/c, F_{r,z}/c]^T$ 이고 $\mathbf{F}(\mathbf{U})$는 해당 방향의 물리적 플럭스 벡터입니다. $\mathbb{P}_r$ 계산 시 $\chi(f)$가 필요합니다.
- 상태 업데이트: 계산된 수치 플럭스를 사용하여 셀 평균값을 업데이트합니다. \(\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{V_i} \sum_{faces} \mathcal{F}_{face} A_{face} + \Delta t \mathbf{S}_i^n\) 여기서 $\mathbf{S}_i$는 소스 항입니다.
- 소스 항 처리: 방사-물질 결합 소스 항($S_{rad}$와 $-\mathbf{G}_{rad}$)은 종종 매우 강성적이므로, 암시적 또는 반-암시적(semi-implicit)으로 처리하여 시간 단계 제한을 완화할 수 있습니다.
3. Realizability
M1 폐쇄에서 중요한 점은 축소된 플럭스 $f=|\mathbf{F}_r|/(cE_r)$가 물리적 범위 $[0, 1]$ 내에 있도록 보장하는 것입니다. 수치적 오류로 인해 이 범위가 위반될 수 있으며, 이를 방지하기 위해 realizability limiter가 필요합니다. 이는 일반적으로 에딩턴 인자 $\chi(f)$를 수정하거나 플럭스를 직접 제한하는 방식으로 적용됩니다. 또한 $E_r \ge 0$ 조건도 유지되어야 합니다.
결론 및 향후 과제
본 문서는 Rad-MHD 방정식의 수치 해법, 특히 FLD와 M1 폐쇄 근사를 중심으로 심층적인 이론적 배경과 수치적 고려 사항을 논의했습니다. FLD는 광학적으로 두꺼운 영역에서 효율적이고 안정적인 해를 제공하는 반면, M1 폐쇄는 복사의 방향성을 고려하여 더 넓은 범위의 광학적 깊이에 적용 가능하지만 계산 비용과 복잡성이 증가합니다.
이러한 방법론들은 천체물리학적 시뮬레이션에서 중요한 도구로 사용되지만, 여전히 많은 도전 과제가 남아있습니다. 예를 들어,
- 효율적인 강성 소스 항 처리: 방사-물질 결합 항은 매우 강성적이어서 전체 계산 시간을 지배할 수 있습니다. 더 효율적인 암시적/IMEX 적분법 개발이 필요합니다.
- 다차원에서의 정확도 및 안정성: 고차 정확도 공간 이산화 기법과 강인한 리만 솔버, 그리고 $\nabla \cdot \mathbf{B}=0$ 조건 유지는 여전히 연구 주제입니다.
- Realizability 문제: M1 및 다른 모멘트 폐쇄 방법에서 물리적 해를 보장하는 것은 중요합니다.
- 불투명도 데이터: 정확하고 광범위한 온도/밀도 범위에 대한 불투명도 데이터베이스의 구축 및 효율적인 사용이 중요합니다.
- 더 발전된 복사 전달 모델: FLD와 M1보다 더 정확하지만 계산적으로 다루기 용이한 새로운 근사법 (예: Variable Eddington Factor, VEF; 혹은 선택된 방향에 대한 이산 종좌표법, Discrete Ordinates $S_N$의 축소 버전)에 대한 연구가 지속되고 있습니다.
참고문헌
- Levermore, C. D. (1984). “Relating Eddington factors to flux limiters.” Journal of Quantitative Spectroscopy and Radiative Transfer, 31(2), 149-160.
- Mihalas, D., & Mihalas, B. W. (1984). Foundations of Radiation Hydrodynamics. Oxford University Press.
- Turner, N. J., & Stone, J. M. (2001). “A module for radiation hydrodynamic calculations with ZEUS-2D using flux-limited diffusion.” The Astrophysical Journal Supplement Series, 135(1), 95.
- González, M., Audit, E., & Huynh, P. (2007). “HERACLES: a three-dimensional radiation hydrodynamics code.” Astronomy & Astrophysics, 464(2), 429-435.