우리는 3. Rocket Nozzle Theory 에서, 같은 로켓엔진이더라도 대기압에 따라 그 성능이 변화한다는 것을 알았다. 진공에서 주로 사용하는 2~3단의 경우에는 사실 큰 문제가 없다. 하지만 1단 엔진의 경우 지표면에서부터 비행기가 날아다니는 고도보다 훨씬 높은 곳까지 그 작동 환경이 크게 변화하므로, 대기압(고도)에 따른 성능 변화를 고려해야 한다.
로켓엔진의 추력에 대한 식을 다시 보자.
$ F = \dot m V_e+A_e(P_e-P_{atm}) $
배기가스 질량유량 $ \dot m $ 의 경우, 엔진 설계 시 상수값으로 특정된다. 노즐 출구 면적 $ A_e $ 또한 상수값을 가진다. 나머지의 변수는 로켓엔진의 작동 환경에 따라 달라지게 된다.
그런데, 3장에서 본 노즐 내부에서 압력에 대한 그래프를 다시 살펴보자.
노즐 확산부에서 노즐 면적이 커질수록 압력이 감소하는데, 노즐이 충분히 확산되지 않아(부족 팽창, $ P_{atm}<P_e $) 노즐 출구 배기가스의 압력이 대기압보다 높다면, 위 식의 $ A_e(P_e-P_{atm}) $ 에 의한 추력이 존재할 것이다. 이는 대기압이 낮은 고고도의 상황을 의미한다.
노즐이 충분히 확산되고 노즐 출구의 압력이 대기압과 일치한다면, 이 상태는 노즐이 해당 고도(압력)에서 최적화되었다고 할 수 있다. 이 상태에서는 고압의 에너지를 최대한으로 가스의 속도에너지로 변환시켜 최고의 추력과 효율을 가진다.
하지만 대기압이 높은 저고도의 상황에서, 계산된 노즐 출구의 압력이 대기압보다 낮다면 어떻게 될까?
노즐 확산부가 너무 크게 확산되어, 노즐 출구의 압력이 대기압보다 낮게 된다면(초과 팽창, $ P_{atm}>P_e $) 효율에 손해를 보는 것은 물론, 대기압에 의해 배기가스가 압축되며 노즐 출구 근방에서 유동 박리가 나타나고 큰 충격파를 유발시킨다. 이는 상단 그림의 C에 해당한다.
최악의 경우, 상단 그림의 A, B와 같이 노즐 내부에서 발생한 충격파로 인해 배기가스가 아음속까지 감속되고 추력을 상실할 수 있다. 이는 엔진은 물론 발사체 시스템 전체에 문제를 야기시킬 수 있기 때문에 가급적이면 이를 회피하도록 엔진이 설계된다.
해당 구간에서는 노즐 확산부의 끝까지 배기가스가 확산되지 못하고, 대기압과 같은 수준까지만 팽창이 되는 것으로 해석할 수 있다. 이 경우, $ A_e(P_e-P_{atm}) $ 에 의한 추력을 상실하고 $ F = \dot m V_e $ 의 추력을 가지는 엔진으로 모델링할 수 있다.
여기서 주목해야 할 것은 출구 배기가스 속도 $ V_e $ 는 출구압력에 대한 함수로 나타내어진다는 점이다. 따라서 우리는 노즐 내의 특정 압력 지점에서의 배기가스의 속도를 알아내면, 아래와 같이 고도에 따른 로켓 엔진의 전 작동 과정에서의 추력을 모델링할 수 있을 것이다.
"Rocket Prupulsion Element", George P. Sutton, WILEY
노즐 내부의 유동이 이상기체 등엔트로피 유동이라고 한다면, 노즐 내의 임의의 위치에서 배기가스의 온도 $ T $ 와 등압비열 $ c_p $ 를 이용해 아래와 같이 나타낼 수 있다.
$ T_c = T+\frac{V}{2c_p} $
여기서 $ T_c $ 는 노즐을 지나기 전 연소실에서의 연소 가스의 온도이다. 여기에 $ c_p=\frac{kR}{k-1} $, $ c^2=kRT $, $ Ma=\frac{V}{c} $ 의 관계식을 대입하면 아래와 같다.
$ \frac{V^2}{2c_pT}=\frac{V^2}{2\left[\frac{kR}{k-1}\right]T}=\left(\frac{k-1}{2}\right)\frac{V^2}{c^2}=\left(\frac{k-1}{2}\right)Ma^2 $
이를 대입하면 아래와 같은 식을 얻을 수 있다.
$ \frac{T_c}{T}=1+\left(\frac{k-1}{2}\right)Ma^2 $
연소실에서의 연소 가스의 온도는 상수값으로 간주할 수 있다. 따라서 위 식을 이용해 노즐 내부의 특정 위치에서의 온도를 계산할 수 있다.
또한 이상 기체는 아래의 식을 따른다.
$ TP^\left(-\frac{k-1}{k}\right)=const. $
이를 위 식에 대입하면, 노즐 내부의 특정 위치에서의 압력을 계산할 수 있다.
$ \frac{P_c}{P}=\left[1+\left(\frac{k-1}{2}\right)Ma^2\right]^{\frac{k}{k-1}} $
$ T_c $ 와 마찬가지로, 연소실에서의 압력 $ P_c $ 또한 상수값으로 간주할 수 있다. KSLV-II 누리 1단에 사용된 KRE-075 엔진의 경우, $T_c=3500$°C, $P_c=6$MPa 의 값을 가진다.
위 식을 이용해 특정 출구 압력에서의 배기가스의 마하수 $ Ma $ 를 계산할 수 있다. 또한 해당 마하수를 바탕으로 온도, 즉 음속을 구할 수 있다. 따라서 위 과정으로 임의의 고도에서의 로켓 엔진의 추력을 계산할 수 있다.
위 방법을 이용한 엔진 모델을 간단히 정리하면 아래와 같다.
1. 노즐이 최적화되었을때의 노즐 출구 압력 $P_{e, opt}$ , 배기가스 속도 $V_{e, opt}$ 를 계산
2. $P_{e, opt} > P_{atm} $ ?
2.1. Yes : $ F = \dot m V_{e, opt} + A_e(P_{e, opt}-P_{atm}) $
2.2. No :
2.2.1. $ \frac{P_c}{P_{atm}}=\left[1+\left(\frac{k-1}{2}\right)Ma^2\right]^{\frac{k}{k-1}} $ 를 만족하는 양수의 $ Ma $ 계산
2.2.2. 위에서 구한 $ Ma $ 를 이용해, $ \frac{T_c}{T_e}=1+\left(\frac{k-1}{2}\right)Ma^2 $ 를 만족하는 $ T_e $ 계산
2.2.3. 위에서 구한 $ Ma $ 와 $ T_e $ 를 이용해 배기가스 출구속도 $ V_e $ 계산
2.2.4. $ F = \dot m V_{e} $
3. End
위에서 구한 방법을 바탕으로, KSLV-II 누리호 1단 엔진에 위 식을 적용해 보자. 공개된 1단의 제원은 아래와 같다.
"75톤급 액체로켓엔진 연소기 기본설계", 한영민 외 6인, 한국추진공학회 2009년도 추계학술대회 논문집 pp.125~129
위 값을 바탕으로, 상술한 엔진 모델을 MATLAB 코드로 작성하면 아래와 같다(첨부파일 참조).
MATLAB Code
1. 노즐이 최적화되었을때의 노즐 출구 압력 $P_{e, opt}$ , 배기가스 속도 $V_{e, opt}$ 를 계산
2. $P_{e, opt} > P_{atm} $ ?
2.1. Yes : $ F = \dot m V_{e, opt} + A_e(P_{e, opt}-P_{atm}) $
2.2. No :
2.2.1. $ \frac{P_c}{P_{atm}}=\left[1+\left(\frac{k-1}{2}\right)Ma^2\right]^{\frac{k}{k-1}} $ 를 만족하는 양수의 $ Ma $ 계산
2.2.2. 위에서 구한 $ Ma $ 를 이용해, $ \frac{T_c}{T_e}=1+\left(\frac{k-1}{2}\right)Ma^2 $ 를 만족하는 $ T_e $ 계산
2.2.3. 위에서 구한 $ Ma $ 와 $ T_e $ 를 이용해 배기가스 출구속도 $ V_e $ 계산
2.2.4. $ F = \dot m V_{e} $
3. End
노즐이 최적화되었을때의 출구 압력은 맨 처음에만 계산하면 되므로, persistent 변수를 이용해 그 값을 저장한다. isp는 엔진의 효율이라고 볼 수 있는 값이며, 5장에서 더 자세히 설명하겠다.
1. 노즐이 최적화되었을때의 노즐 출구 압력 $P_{e, opt}$ , 배기가스 속도 $V_{e, opt}$ 를 계산
2. $P_{e, opt} > P_{atm} $ ?
2.1. Yes : $ F = \dot m V_{e, opt} + A_e(P_{e, opt}-P_{atm}) $
2.2. No :
2.2.1. $ \frac{P_c}{P_{atm}}=\left[1+\left(\frac{k-1}{2}\right)Ma^2\right]^{\frac{k}{k-1}} $ 를 만족하는 양수의 $ Ma $ 계산
2.2.2. 위에서 구한 $ Ma $ 를 이용해, $ \frac{T_c}{T_e}=1+\left(\frac{k-1}{2}\right)Ma^2 $ 를 만족하는 $ T_e $ 계산
2.2.3. 위에서 구한 $ Ma $ 와 $ T_e $ 를 이용해 배기가스 출구속도 $ V_e $ 계산
2.2.4. $ F = \dot m V_{e} $
3. End
Atmosphere_model 함수는 현재 고도[m] 를 입력받아 대기압[kPa] 를 출력하는 함수이다. 첨부파일 참조.
1. 노즐이 최적화되었을때의 노즐 출구 압력 $P_{e, opt}$ , 배기가스 속도 $V_{e, opt}$ 를 계산
2. $P_{e, opt} > P_{atm} $ ?
2.1. Yes : $ F = \dot m V_{e, opt} + A_e(P_{e, opt}-P_{atm}) $
2.2. No :
2.2.1. $ \frac{P_c}{P_{atm}}=\left[1+\left(\frac{k-1}{2}\right)Ma^2\right]^{\frac{k}{k-1}} $ 를 만족하는 양수의 $ Ma $ 계산
2.2.2. 위에서 구한 $ Ma $ 를 이용해, $ \frac{T_c}{T_e}=1+\left(\frac{k-1}{2}\right)Ma^2 $ 를 만족하는 $ T_e $ 계산
2.2.3. 위에서 구한 $ Ma $ 와 $ T_e $ 를 이용해 배기가스 출구속도 $ V_e $ 계산
2.2.4. $ F = \dot m V_{e} $
3. End
부족 팽창의 경우는 식을 간단하게 나타낼 수 있다.
1. 노즐이 최적화되었을때의 노즐 출구 압력 $P_{e, opt}$ , 배기가스 속도 $V_{e, opt}$ 를 계산
2. $P_{e, opt} > P_{atm} $ ?
2.1. Yes : $ F = \dot m V_{e, opt} + A_e(P_{e, opt}-P_{atm}) $
2.2. No :
2.2.1. $ \frac{P_c}{P_{atm}}=\left[1+\left(\frac{k-1}{2}\right)Ma^2\right]^{\frac{k}{k-1}} $ 를 만족하는 양수의 $ Ma $ 계산
2.2.2. 위에서 구한 $ Ma $ 를 이용해, $ \frac{T_c}{T_e}=1+\left(\frac{k-1}{2}\right)Ma^2 $ 를 만족하는 $ T_e $ 계산
2.2.3. 위에서 구한 $ Ma $ 와 $ T_e $ 를 이용해 배기가스 출구속도 $ V_e $ 계산
2.2.4. $ F = \dot m V_{e} $
3. End
초과 팽창의 경우에는, 위 알고리즘을 바탕으로 비선형 방정식을 풀어 그 해를 찾아내야 한다. 이를 위해 MATLAB의 syms 와 수치해석적 solver 함수인 vpasolve 를 사용하였다.
위 식에서 $ Ma^2 $ 이 사용되므로 해당 식을 그대로 사용할 경우 vpasolve 로 구한 해는 양수와 음수가 존재한다. 두 가지 모두 수식적으로는 맞는 해 이지만, 물리적으로는 노즐 내의 속도가 음수가 될 수는 없으므로 음수의 경우를 제외해야 한다. 이를 위해 $ Ma^2 $ 의 값을 Ma_sq 의 변수로 지정하고 실제 $ Ma $ 값은 Ma_sq 변수에 제곱근을 취하여 그 값을 계산하였다.
KSLV-II 누리호 1단의 경우, KRE-075 엔진 4기를 클러스터링 하였으므로 최종적으로 구한 추력에 4배를 적용하여, Rocket Dynamics 에 적용하였다.
동역학 모델은 오일러 1차 근사를 통한 적분을 사용하였다.
누리호 1단의 경우 127초 동안 작동하므로 같은 시간 동안 시뮬레이션을 진행하였으며 그 결과는 아래와 같다.
처음에는 고도가 낮기에 추력과 isp 가 낮은 값을 가지지만, 엔진을 사용하고 점차 고도가 올라가며 추력과 비추력이 커지는 모습을 볼 수 있다.
본 시뮬레이션은 로켓이 계속 수직으로 날아가는 것을 상정하였으나, 실제 누리호의 궤적은 비행 중 남쪽을 향해 비행방향을 꺾는다.
그렇기에 실제 누리호 1단 종료 시점에서의 비행 고도는 약 60 km 이지만, 본 시뮬레이션에서는 그보다 높은 약 80km 의
고도를 가진다.
그러나 로켓이 받는 가속도는 실제 누리호 비행 궤적과 거의 흡사한 형태를 가진다. 누리호의 1단 종료 시점 최대 가속도는 3.3G로, 본 시뮬레이션으로 계산된 3.2G와 큰 차이를 가지지 않는다.
본 글을 통해 작동 환경, 특히 고도에 따른 로켓 엔진의 출력 변화를 이론적으로 유도하고, 해당 모델을 MATLAB 을 통해 구현하여 실제 로켓과 유사한 형태를 가짐을 확인하였다. 이제 해당 모델을 사용하여 로켓의 동역학 모델링을 구현할 수 있을 것이다.
Comments
Post a Comment