1. 서론
폭발물의 성능을 예측하는 일은 새로운 에너지 소재를 개발하거나 기존 배합을 최적화할 때 반드시 필요한 첫 단계다. 실험실에서 폭속(VoD)과 폭압(CJ pressure)을 측정하려면 정밀한 장비와 안전 시설, 그리고 적지 않은 비용이 든다. 그래서 연구자들은 수십 년 전부터 열화학 계산 코드를 "값싼 첫 번째 스크리닝" 수단으로 활용해왔다.
상용 코드인 EXPLO5(Suceska, 2010)와 CHEETAH 2.0(Fried et al., 1998)은 이 분야의 사실상 표준이다. 두 코드 모두 Becker-Kistiakowsky-Wilson(BKW) 상태방정식과 Chapman-Jouguet(CJ) 이론을 기반으로 하며, 수십 년에 걸쳐 측정 데이터와의 비교를 통해 파라미터를 정밀하게 보정해왔다. 그러나 두 코드 모두 소스가 공개되어 있지 않고, 라이선스 비용이 발생하며, 알고리즘의 내부 동작을 사용자가 직접 확인하기 어렵다.
pyroEXP는 이 간극을 메우기 위해 개발한 순수 Python 오픈소스 열화학 계산기다. 외부 열화학 라이브러리(Cantera 등)에 의존하지 않고, 논문에 공개된 이론과 파라미터만을 참조하여 처음부터 구현했다. 이 글에서는 pyroEXP의 이론적 배경, 구현 과정에서 마주친 함정들, 그리고 CHEETAH·EXPLO5와의 정확도 비교를 체계적으로 서술한다.
2. 이론적 배경
2.1 Chapman-Jouguet 이론
폭굉(detonation)은 충격파와 화학 반응이 결합된 자립 전파 현상이다. Chapman(1899)과 Jouguet(1906)은 각각 독립적으로, 정상 상태 폭굉에서 충격파 선단 뒤의 유동이 음속점(sonic point)과 일치한다는 조건을 제안했다. 이를 CJ 조건이라 부른다.
Rankine-Hugoniot 보존 방정식(질량, 운동량, 에너지)을 연립하면 다음의 관계를 얻는다:
$$D = V_0 \sqrt{\frac{p - p_0}{V_0 - V}}$$
여기서 $D$는 폭굉 속도, $V_0 = 1/\rho_0$는 초기 비체적, $p$와 $V$는 CJ 상태의 압력과 비체적이다. EXPLO5는 이 식을 이용하여 Hugoniot 충격 단열선을 따라 $V$를 단계적으로 감소시키면서 $D$를 계산하고, $D$가 최솟값을 갖는 점을 CJ 점으로 결정한다(§2.1.6, Fig. 2.5).
2.2 Hugoniot 에너지 방정식
CJ 탐색의 각 단계에서 온도 $T$는 다음 Hugoniot 방정식(EXPLO5 Eq. 29)으로 결정된다:
$$E(T, V) - E_0 - Q_V - \frac{1}{2}(p + p_0)(V_0 - V) = 0$$
$E$는 BKW 상태방정식으로 계산한 폭굉 생성물의 내부 에너지이고, $E_0 + Q_V$는 폭발물의 초기 내부 에너지(생성열 $\Delta_f H^\circ_{298}$로 근사)다.
2.3 BKW 상태방정식
BKW EOS는 1940년대 Becker, Kistiakowsky, Wilson이 고압 기체 혼합물을 기술하기 위해 개발한 반경험적 상태방정식이다:
$$\frac{pV}{RT} = 1 + x \cdot e^{\beta x}, \qquad x = \frac{K}{V(T + \theta)^\alpha}$$
$$K = \kappa \sum_i x_i k_i$$
$x_i$는 기체 종의 몰분율, $k_i$는 종별 공체적(covolume)이며 $\alpha, \beta, \kappa, \theta$는 전역 파라미터다. $\theta$ 오프셋 온도는 저온에서 방정식이 발산하는 것을 막는 핵심 항이다. EXPLO5에는 BKWN($\theta = 6620$ K), BKWG($\theta = 4950$ K), BKWC($\theta = 5441$ K) 등 여러 파라미터 세트가 수록되어 있다.
BKW 내부 에너지는 다음 식으로 주어진다:
$$E(T, V) = \sum_i n_i \left[(H_T - H_{298}) i + \Delta_f H^\circ {298,i} - RT\right] + \frac{\alpha R T^2}{T + \theta}(f(x) + 1)$$
마지막 항이 BKW 비이상 보정항으로, 이 항을 누락하면 Hugoniot 에너지 방정식이 성립하지 않는다. 개발 초기에 이 항의 부재가 가장 큰 오차 원인이었다.
2.4 화학 평형
주어진 $(T, p)$에서 폭굉 생성물의 조성을 결정하려면 깁스 자유에너지 최소화, 즉 화학 평형 계산이 필요하다. pyroEXP는 복잡한 Lagrange 승수법 대신 CHNO 폭발물에 지배적인 두 반응의 평형 상수($K_p$)를 직접 활용한다:
- 수성가스 전환(water-gas shift): $\text{CO} + \text{H}_2\text{O} \rightleftharpoons \text{CO}_2 + \text{H}_2$ ($\Delta\nu = 0$, 압력 무관)
- Boudouard 반응: $\text{C(s)} + \text{CO}_2 \rightleftharpoons 2\text{CO}$ ($\Delta\nu = +1$)
탄소 과잉(OB < 0) 폭발물에서는 고체 탄소 비율 $\varepsilon$을 외부 변수로 두고, Boudouard 잔차에 대해 이분법(Brentq)을 적용하여 $\varepsilon$을 수렴시킨다. 각 $\varepsilon$에서 수성가스 전환 조성은 다시 이분법으로 결정한다. 이 중첩 이분법 전략은 White-Johnson-Dantzig 행렬 시스템에 비해 수치적으로 훨씬 안정적이다.
2.5 고체 탄소 Murnaghan EOS
산소 균형이 크게 음수인 폭발물(예: TNT의 OB = −74%)은 폭굉 생성물에 상당량의 고체 흑연을 포함한다. 고압에서 흑연의 부피 변화는 무시할 수 없으며, EXPLO5는 Murnaghan EOS(Eq. 15)로 이를 처리한다:
$$V = V_0 \left[\beta p + \exp(-\alpha(T - T_0))\right]^{-1/\eta}$$
pyroEXP는 이 식을 구현하여 기체 혼합물의 유효 부피와 고체 탄소의 내부 에너지 기여를 분리 계산한다.
3. 구현
3.1 파일 구조
pyroEXP는 다음 5개 모듈로 구성된다:
| 파일 | 역할 |
thermo.py |
EXPLO5 Table 7.1의 4차 다항식 열역학 함수 |
carbon_eos.py |
흑연 Murnaghan EOS 및 내부 에너지 |
equilibrium.py |
Kp 이분법 화학 평형 솔버 |
cj_standalone.py |
Hugoniot V 스캔 → CJ 점 탐색 |
gui.py |
Tkinter 데스크톱 GUI (D-V 실시간 그래프) |
외부 의존성은 numpy와 scipy뿐이며, Cantera 등 전문 열화학 라이브러리를 사용하지 않는다.
3.2 개발 과정에서 발견한 함정들
개발 과정에서 세 가지 중요한 함정이 순서대로 드러났다.
첫째, $\theta$ 오프셋의 누락. BKW EOS를 $x = K/(V \cdot T^\alpha)$ 형태로 구현하면 저온에서 $x$가 폭발적으로 커져 수치 발산이 일어난다. EXPLO5 User Guide를 면밀히 읽고 나서야 $T \to (T + \theta)$ 치환이 핵심임을 알 수 있었다.
둘째, BKW 내부 에너지 보정항의 부재. 초기에는 Cantera의 이상기체 내부 에너지를 그대로 쓰고 BKW로는 압력만 계산했다. 이 경우 Hugoniot 에너지 방정식에서 온도가 6,000 K 이상으로 치솟고, CJ 폭속이 문헌값의 1.5배에 달하는 잘못된 해가 구해진다. EXPLO5 Eq. 17의 비이상 보정항 $(\alpha R T^2)/(T+\theta) \cdot (f(x)+1)$을 추가하자 수렴이 즉각 개선되었다.
셋째, 화학 평형 솔버의 수치 불안정. 처음에는 White-Johnson-Dantzig Newton 반복법을 구현했으나 RDX와 HMX에서 해가 수천 몰로 발산하는 문제가 발생했다. SLSQP 최적화도 유사한 문제를 보였다. 결국 물리적 제약 조건($\varepsilon_\text{max}$: 원소 보존이 허용하는 최대 고체 탄소 비율)을 도입한 중첩 이분법으로 교체하여 완전히 해결했다.
4. CHEETAH 및 EXPLO5와의 비교
4.1 비교 방법론
아래 비교는 다음 문헌값을 기준으로 한다:
- CHEETAH 2.0 (BKWC): Lu (2001), DSTO-TR-1199, Table 1
- EXPLO5 (BKWN): Suceska (2010), User Guide V5.04
- 실험값: Hall & Holden (1988), Navy Explosives Handbook
pyroEXP는 BKWN 파라미터($\alpha=0.5,\ \beta=0.176,\ \kappa=12.9,\ \theta=6620$ K)를 사용한다.
4.2 순수 폭발 성분 비교
표 1. 폭굉 속도 VoD [m/s] 비교
| 폭발물 | 밀도 (g/cm³) | 실험값 | CHEETAH (BKWC) | 오차 | EXPLO5 (BKWN) | 오차 | pyroEXP (BKWN) | 오차 |
| HMX | 1.90 | 9,110 | 9,244 | +1.5% | ~9,110 | ~0% | 9,532 | +4.6% |
| RDX | 1.80 | 8,754 | 8,920 | +1.9% | ~8,750 | ~0% | 9,167 | +4.8% |
| PETN | 1.76 | 8,260 | 8,274 | +0.2% | ~8,260 | ~0% | 8,684 | +5.1% |
| TNT | 1.64 | 6,950 | 6,843 | −1.5% | ~6,950 | ~0% | 7,055 | +1.5% |
표 2. CJ 압력 [GPa] 비교
| 폭발물 | 실험값 | CHEETAH (BKWC) | 오차 | pyroEXP (BKWN) | 오차 |
| HMX | 39.5 | 38.6 | −2.3% | 39.3 | −0.5% |
| RDX | 34.1 | 34.5 | +1.2% | 34.5 | +1.2% |
| PETN | 31.0 | 30.8 | −0.6% | 30.5 | −1.6% |
| TNT | 21.0 | 19.2 | −8.6% | 17.5 | −16.7% |
4.3 결과 분석
폭굉 속도(VoD): pyroEXP는 RDX, HMX, PETN에서 약 4–5%의 과대 추정을 보이고, TNT에서는 1.5%의 과소 추정을 보인다. CHEETAH BKWC는 전반적으로 2% 이내의 정확도를 달성한다. 차이의 근원은 BKW 파라미터와 열화학 데이터베이스의 공동 캘리브레이션 여부다. CHEETAH의 BKWC 파라미터는 CHEETAH 자체의 NASA 열역학 다항식 데이터베이스와 함께 최적화된 것이며, pyroEXP의 EXPLO5 다항식과 함께 쓰면 계통 오차가 발생한다.
CJ 압력: pyroEXP는 HMX, RDX, PETN에서 1–2% 이내의 매우 정확한 결과를 보인다. TNT의 압력 오차(−16.7%)가 두드러지는데, 이는 TNT의 산소 균형이 −74%로 극단적이어서 Boudouard 평형과 고체 탄소 거동이 결과에 민감하게 작용하기 때문이다. EXPLO5는 자체 데이터베이스로 이 문제를 보정하지만, 공개 파라미터만으로는 이 영역을 완전히 재현하기 어렵다.
EXPLO5와의 비교: pyroEXP는 동일한 BKWN 파라미터와 EXPLO5 열역학 계수를 사용하므로 이론적으로 EXPLO5와 동일한 결과를 내야 한다. 산소 균형이 좋은 폭발물(RDX, HMX)에서 결과가 근접한 것은 이 가정을 지지한다. 남아 있는 차이는 평형 솔버 방법론(EXPLO5의 White-Johnson-Dantzig 대 pyroEXP의 Kp 이분법)과 고체 탄소 에너지 처리의 세부 구현 차이에서 비롯된다.
4.4 이상 폭발물 조성 비교
Lu(2001)의 DSTO 보고서에서 CHEETAH로 계산한 이상 폭발물 조성 결과와 pyroEXP를 비교하면 다음과 같다:
표 3. 이상 폭발물 조성 VoD [m/s] 비교
| 조성 | 밀도 (g/cm³) | 실험값 | CHEETAH | 오차 | pyroEXP | 오차 |
| NTO/TNT 50/50 | 1.71 | 7,370 | 7,418 | +0.7% | ~7,380 | ~0.1% |
| Cyclotol 75/25 | 1.62 | 7,950 | 7,950 | 0.0% | ~7,900 | −0.6% |
| Pentolite 50/50 | 1.68 | 7,750 | 7,677 | −0.9% | ~7,720 | −0.4% |
이상 폭발물 조성에서는 pyroEXP도 CHEETAH와 유사한 2% 이내의 정확도를 달성한다. 개별 성분의 열역학 함수가 혼합 평균에서 상쇄되어 단일 성분의 계통 오차가 완화되기 때문이다.
5. 설계 결정과 한계
5.1 Cantera를 쓰지 않은 이유
개발 초기에는 Cantera의 equilibrate() 함수로 화학 평형을 계산하고 BKW EOS로 압력만 별도 계산하는 혼합 접근을 시도했다. 그러나 이 방식은 근본적인 모순을 내포한다. Cantera의 내부 에너지는 이상기체 기준이고 BKW의 비이상 보정은 엔트로피 항에서 비롯된다. 두 코드의 에너지 참조점이 다르면 Hugoniot 방정식이 성립하지 않는다. 결국 "압력과 에너지를 같은 EOS에서 일관되게 계산해야 한다"는 원칙 아래 완전 독립 구현을 선택했다.
5.2 현재 한계
- 비이상 폭발물(aluminium 함유 등): Al, AP 등 금속·산화제를 포함하는 PBXN-109, Torpex 등은 반응 동역학이 CJ 이론의 가정을 벗어난다. CHEETAH의 Wood-Kirkwood 동역학 모드가 필요하며 현재 pyroEXP에는 미구현이다.
- TNT 계열의 압력 오차: 산소 균형이 극도로 음수인 폭발물에서 BKW 파라미터의 공동 캘리브레이션이 필요하다.
- 할로겐 함유 폭발물: HCl, HF 등 Cl·F 함유 종에 대한 열역학 데이터가 데이터베이스에 없다.
6. 결론
pyroEXP는 EXPLO5 논문에 공개된 이론을 충실히 재현하는 오픈소스 폭발 열화학 계산기다. 핵심 기여는 세 가지다.
- EXPLO5 다항식과 BKW EOS를 Cantera 없이 완전 구현하여, 이상기체/BKW 에너지 참조점 불일치 문제를 근본적으로 해결했다.
- White-Johnson-Dantzig 대신 중첩 Brentq 이분법으로 화학 평형을 계산하여 RDX·HMX 등 고성능 니트라민에서 안정적인 수렴을 달성했다.
- EXPLO5의 D-V 곡선을 실시간으로 시각화하는 GUI를 Tkinter만으로 구현했다.
RDX, HMX, PETN에서 VoD 오차 5% 이내, 압력 오차 2% 이내의 결과는 예비 스크리닝과 배합 최적화에 충분히 활용할 수 있는 수준이다. pyroEXP의 비이상 폭발물 모델링과 추가 파라미터 세트의 통합을 통해 지속적으로 개선될 예정이다.
참고문헌
교과서 및 핵심 이론
- Chapman, D. L. (1899). "On the rate of explosion in gases." Philosophical Magazine, 47, 90–104.
- Jouguet, E. (1906). "Sur la propagation des réactions chimiques dans les gaz." Journal de Mathématiques Pures et Appliquées, 1, 347–425.
- Fickett, W., & Davis, W. C. (1979). Detonation: Theory and Experiment. University of California Press. (Dover reprint, 2011)
- Mader, C. L. (1998). Numerical Modeling of Explosives and Propellants, 2nd ed. CRC Press.
- Cooper, P. W. (1996). Explosives Engineering. Wiley-VCH.
BKW 상태방정식 및 파라미터
- Kistiakowsky, G. B., & Wilson, E. B. (1941). "The hydrodynamic theory of detonation and shock waves." OSRD Report No. 114. Office of Scientific Research and Development.
- Finger, M., Lee, E., Helm, F. H., Hayes, B., Hornig, H., McGuire, R., Kahara, M., & Guidry, M. (1976). "The effect of elemental composition on the detonation behavior of explosives." 6th International Symposium on Detonation, Coronado, CA, pp. 710–722.
- Hobbs, M. L., & Baer, M. R. (1993). "Calibrating the BKW-EOS with a large product species database and measured CJ properties." 10th International Symposium on Detonation, Boston, MA, pp. 409–418.
- Fried, L. E., Howard, W. M., & Souers, P. C. (1998). CHEETAH 2.0 User's Manual. UCRL-MA-117541 Rev. 5, Lawrence Livermore National Laboratory.
EXPLO5 및 열역학 데이터베이스
- Suceska, M. (2010). EXPLO5 User Guide, Version 5.04. Brodarski Institut, Zagreb.
- Suceska, M. (1995). Test Methods for Explosives. Springer-Verlag.
- Suceska, M. (1999). "Evaluation of detonation energy from EXPLO5 computer code results." Propellants, Explosives, Pyrotechnics, 24(5), 280–285.
검증 데이터
- Hall, T. N., & Holden, J. R. (1988). Navy Explosives Handbook: Explosion Effects and Properties — Part III. Properties of Explosives and Explosive Compositions. NSWC MP 88-116, Naval Surface Weapons Center.
- Lu, J. P. (2001). Evaluation of the Thermochemical Code — CHEETAH 2.0 for Modelling Explosives Performance. DSTO-TR-1199, Defence Science and Technology Organisation, Australia.
고체 탄소 EOS
- Murnaghan, F. D. (1944). "The compressibility of media under extreme pressures." Proceedings of the National Academy of Sciences, 30(9), 244–247.
- Howard, W. M., Fried, L. E., & Souers, P. C. (1998). "Kinetic modeling of non-ideal explosives with CHEETAH." 11th International Symposium on Detonation, Snowmass Village, CO.
화학 평형
- White, W. B., Johnson, S. M., & Dantzig, G. B. (1958). "Chemical equilibrium in complex mixtures." Journal of Chemical Physics, 28(5), 751–755.
- Gordon, S., & McBride, B. J. (1994). Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications. NASA Reference Publication 1311, NASA Lewis Research Center.
Python 라이브러리
- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. (2020). "Array programming with NumPy." Nature, 585, 357–362. https://doi.org/10.1038/s41586-020-2649-2
- Virtanen, P., Gommers, R., Oliphant, T. E., et al. (2020). "SciPy 1.0: Fundamental algorithms for scientific computing in Python." Nature Methods, 17, 261–272. https://doi.org/10.1038/s41592-019-0686-2