고체 로켓 노즐의 경계층 해석과 유한차분법을 이용한 탄소/페놀릭의 열반응 해석 연구

2.1 로켓 노즐 코어

로켓 노즐의 유동은 크게 코어유동과 경계층유동으로 나눌 수 있다. 그 중 로켓 노즐의 코어유동은 Euler 식을 통한 수치계산으로 구할 수 있다. 또는 이론적 접근을 통한 대수 해를 이용하여 계산하여도 복잡하지 않은 노즐의 경우 잘 맞는 것으로 알려져 있다. 이를 열량적 완전기체의 1차원 등엔트로피 유동으로 가정하여 얻는 Eq. 1과 같다[4].

Eq. 1에서 A는 임의 위치에서의 노즐 단면적, A_th는 노즐목 단면적, M은 마하수, γ는 비열비를 나타낸다. 여기서 마하수는 2분법(Bisection Method)을 통하여 구할 수 있다[5].

Fig. 1은 고체 로켓 모터 노즐 내면의 형상을 나타낸 그림이며, r_o은 임의위치에서의 노즐 반경, r_th은 노즐목 반경을 나타낸다. 주어진 형상을 통해 각 위치에서 연소가스의 마하수를 Eq. 1을 통해 구할 수 있다. 각 위치에서 얻어진 마하수와 Eq. 2 ~ Eq. 4를 통해 온도 및 연소가스 속도를 구할 수 있다.

Fig. 1

Coordinate of rocket nozzle.

Eq. 2는 에너지방정식 및 음속을 이용하여 유도된 식이며, 마하수를 이용한 노즐 온도분포식이다. Eq. 3은 노즐목에서의 속도를 나타내며, Eq. 4는 위치별 전엔탈피 값이 동일하다는 가정으로부터 유도된 속도분포를 나타내는 식이다. 여기서 T는 온도, u는 속도를 나타내고, th는 노즐목의 값, o은 연소관을 의미한다. 노즐 코어유동에서의 속도분포 및 온도분포는 경계층 해석에 필요한 자료가 된다.

2.2 로켓 노즐 경계층유동

노즐 경계층유동 해석에서 좌표계는 축대칭 좌표계를 사용하면 편리하다. Fig. 2는 노즐에서의 좌표계를 나타내며, x는 노즐표면의 접선방향, y는 수직방향, z는 축방향, r은 반경방향, r_o은 노즐표면까지의 반경길이를 나타낸다.

Fig. 2

Coordinate of rocket nozzle.

로켓 노즐 유동은 압축성 난류유동이며, 이에 대해 마하수 5 이하에 적용할 수 있는 Morkovin의 가설을 적용할 수 있다[6]. Morkovin 가설 적용 및 경계층유동 가정을 적용한 연속 방정식, 운동량 방정식 및 에너지 방정식은 Eq. 5 ~ Eq. 7로 정리된다. 축대칭 좌표계를 이용할 경우 물체의 길이 및 반경에 비하여 경계층 두께가 매우 작다는 가정이 적용된 것이며, 이 가정은 노즐 유동에서 타당하다. 여기서 ρ는 연소가스 밀도, u는 접선방향 속도, v는 표면 수직방향 속도, p는 압력, τ_eff는 유효전단응력, J는 전엔탈피, q_eff는 유효열유속을 나타내며, ^-과 ′은 각각 시간평균과 교란값을 의미한다.

축대칭 경계층 미분방정식에 코어유동 조건을 적용하여 노즐 표면으로부터 경계층(y=0~max(δ,δ_T))까지 적분하여 정리가 가능하다. 연속방정식 적분 후 발생하는 변수를 이용하여 운동량방정식와 에너지방정식을 적분하면 각각 Eq. 8과 Eq. 9로 정리된다. 편의상 ^-을 생략하여 정리하였다. 여기서 ρ는 연소가스 밀도, u와 v는 각각 x방향과 y방향 속도, θ는 운동량 두께, H는 형상계수, C_f는 표면마찰계수, B는 운동량 경계층 적분방정식에서의 분출계수, J는 전엔탈피, θ_T는 엔탈피 두께, St는 Stanton수, B_T는 에너지 경계층 적분방정식에서 분출계수, 하첨자 e는 코어유동, 하첨자 B는 분출가스를 나타낸다.

여기서

$$ \begin{array}{c}{\delta }^{*}={\int }_0^{{\left(\delta ,{\delta }_T\right)}_{\text{max}}}\left(1-\frac{u}{u_e}\right)dy,\hfill \\ \theta ={\int }_0^{{\left(\delta ,{\delta }_T\right)}_{\text{max}}}\frac{\rho u}{{\rho }_e{u}_e}\left(1-\frac{u}{u_e}\right)dy,\hfill \\ \begin{array}{c}{\theta }_T={\int }_0^{{\left(\delta ,{\delta }_T\right)}_{\text{max}}}\frac{\rho u}{{\rho }_e{u}_e}\left(1-\frac{J-J_w}{J_e-J_w}\right)dy,\hfill \\ H=\frac{{\delta }^{*}}{\theta },C_f=\frac{{\tau }_w}{{\rho }_e{u}_e^2}, St=-\frac{q_w}{{\rho }_e{u}_e^2\left(J_{aw}-J_w\right)},\hfill \\ B=-\frac{1}{C_f}\frac{{\rho }_B{v}_B}{\rho u}, B_T=-\frac{1}{St}\frac{{\rho }_B{v}_B}{\rho u}\hfill \end{array}\hfill \end{array}$$

결과적으로 압축성 난류유동의 경계층 적분방정식은 운동량 두께(θ)와 엔탈피 두께(θ_T)에 대한 미분방정식으로 표현된다. 난류유동에서 비압축성 유동에 대해 잘 정리되어 있지만, 압축성 유동에 대해서는 비압축성 유동에 비해 아직 충분히 완전한 연구가 진행되지 않았다. 본 논문에서는 비압축성 난류유동으로 정리된 표면마찰계수(C_f)와 Stanton수(St)의 반경험식을 기준온도방법을 이용하여 압축성 난류유동에 적용하여 계산을 수행하였다. 경계층 적분방정식을 통하여 얻은 운동량 두께 및 엔탈피 두께를 이용하여 표면마찰계수 및 Stanton수를 구할 수 있다. 이로부터 대류열전달계수를 구할 수 있다.

2.3 경계층 적분방정식 수치해석

우선 경계층 해석에 있어 필요한 데이터는 노즐 내면 형상이며, 축방향(z) 및 노즐 반경길이(r_o)로 표현이 가능하다. Fig. 3은 노즐 내면형상을 축방향 및 노즐 반경길이 격자점으로 표현한 것이다.

Fig. 3

Grids for inner contour of rocket nozzle.

Eq. 10은 축방향(z)과 노즐 반경길이(r_o)로 표면의 접선방향(x)을 나타낸 식이다. Eq. 11은 Eq. 10을 적분하여 접선길이를 근사화한 것이다. 본 논문에서의 미분은 Eq. 12와 같이 근사화 하였고, 적분은 Simpson’s rule을 이용하였다[5].

압축성 난류 운동량 경계층 적분방정식인 Eq. 8과 에너지 적분방정식인 Eq. 9를 통하여 운동량 두께와 엔탈피 두께를 구한 후 표면에서 대류열전달계수를 얻을 수 있다. 이때 형상계수(H), 표면마찰계수(C_f) 및 Stanton수(St)는 운동량 두께 및 엔탈피 두께의 함수이지만, 초기값 설정으로 상수로 간주하면 미분방정식을 통해 운동량 두께와 엔탈피 두께를 계산할 수 있다. 구해진 운동량 두께와 엔탈피 두께를 이용하여 H, C_f 및 St을 갱신(update)할 수 있다. 갱신된 변수를 이용하여 다시 운동량 두께 및 엔탈피 두께를 재계산할 수 있고, 이와 같은 방법을 반복하여 운동량 두께 및 엔탈피 두께를 계산하였다.

Fig. 4는 전술한 계산 흐름도를 나타낸다. 여기서 i는 노즐 내면 격자번호, k는 운동량 두께 및 엔탈피 두께 반복계산 횟수를 의미한다. 여기서 형상계수의 초기조건은 전형적인 난류유동에서의 값으로 설정하였다[6]. Eq. 13은 Stanton수와 대류열전달계수의 관계를 나타내며, Fig. 4를 통해 구해진 Stanton수를 Eq. 13에 적용하여 대류열전달계수를 얻을 수 있다. 여기서 α는 온도기준 대류열전달계수이다. 이는 노즐 내열재료의 열전도 해석의 경계조건으로 이용된다.

Fig. 4

Process for calculation of momentum thickness and enthalpy thickness.

3.1 탄소/페놀릭 열반응

고체 추진기관에서 탄소섬유에 페놀 수지가 함침된 탄소/페놀릭 복합재료를 노즐 내열재료로 사용한다. 페놀 수지를 이용하여 탄소섬유의 높은 열전도도를 낮춤(단열성능)으로 단열재료로도 사용이 가능하다. 고온의 연소가스로 인해 노즐 재료는 열반응(thermal response)이 발생하므로 이에 대한 예측이 필요하다. 로켓 추진기관 연소 시 발생하는 연소가스와 탄소/페놀릭의 열반응은 (1) 373 K에 재료에 포함된 수분 증발이 발생하며, (2) 570 K~600 K에서 페놀 수지가 열분해(thermal decomposition)되어 숯(char) 생성 및 열분해가스(Pyrolysis Gas)를 발생시키고, (3) 900 K ∼ 1000 K에서 열분해가 완료되어 표면에 숯이 형성된 후 숯 및 탄소섬유가 연소가스에 의해 표면 삭마가 발생한다. Fig. 5는 탄소/페놀릭의 열반응 과정을 나타낸다. 열반응이 발생한 탄소/페놀릭 재료는 표면으로부터 숯층(Char Layer), 열분해층(Decomposition Zone), 처녀층(Virgin Layer)으로 이루어진다.

Fig. 5

Thermal response of Carbon/Phenolic.

열분해를 고려한 탄소/페놀릭 복합재료의 축대칭 열전도 지배방정식은 Eq. 14와 같다[8]. 여기서 c는 비열, λ는 열전도계수, Q는 흡열량, χ는 바인더의 분해비율, 하첨자 m은 탄소/페놀릭 재료, 하첨자 d는 열분해층, 하첨자 G는 열분해가스(Gas), 하첨자 B는 페놀 수지 바인더(Binder), 하첨자 dA와 dB는 각각 열분해층 전면과 후면, 상첨자 A는 겉보기 밀도(Bulk density)가 아닌 열분해가스 자체의 밀도를 의미한다. Eq. 15는 열분해가스 질량 유속을 나타내는 식이며, 계산을 위해 적용된 물성값 및 열분해 가스 특성 관련 값은 Table 1과 같다. 탄소/페놀릭 재료의 랩핑(warpping)에 따른 열전도도 변화는 래핑 각도에 따른 변환식을 적용하였다[9]. 여기서 탄소/페놀릭 처녀밀도는 1,340 (kg/m³)이었다.

Table 1.

Material properties of Carbon / Phenolic and pyrolysis gas.

탄소/페놀릭 재료가 연소가스에 노출 시 발생하는 열반응은 열분해뿐만 아니라 표면 삭마도 발생한다. 탄소/페놀릭의 표면 삭마는 크게 화학적 삭마, 기계적 삭마, 물리적 삭마로 구분 할 수 있다. 화학적 삭마는 연소가스와 숯의 산화 반응에 의한 삭마를 의미하고, 기계적 삭마는 내열재가 용융되어 유동장에 의한 전단, 입자 충돌 및 침식에 의한 삭마, 물리적 삭마는 용융, 증발 또는 승화에 의한 표면 소실을 의미한다. 탄소/페놀릭 재료는 화학적 삭마가 지배적인 것으로 알려져 있으며, Eq. 16과 Eq. 17은 연소가스 중 화학적 삭마의 주요 성분인 H₂O, CO₂와 탄소(C)의 반응식이다[10]. 여기서 상첨자 S는 고체, 상첨자 G는 기체를 의미한다.

탄소/페놀릭 재료의 화학적 삭마속도는 화학양론 관계식과 질량전달 관계식의 등치를 통하여 구할 수 있다. 탄소/페놀릭 열분해가스의 영향을 고려하지 않은 표면 삭마속도는 Eq. 18로 귀결된다[11]. 열분해가스가 삭마에 미치는 영향은 엔탈피 기준 대류열전달계수 대비 열분해가스의 질량유속이 작으면 안전하게 무시하여 사용할 수 있다[12,13]. 여기서 Y_C,ω는 내열재 표면층의 탄소의 질량분율, M_C는 탄소의 분자량, M_mix,e은 연소가스 혼합물의 평균분자량, X_H2O,_e와 X_CO2,_e는 각각 H₂O와 CO₂의 질량분율, B_m은 산화포텐셜, $$ \overline{R}$$은 기체상수, T_ω는 노즐벽 온도, P_ω는 압력, E_C와 K_C는 각각 표면산화반응에서의 활성화에너지와 반응상수, 하첨자 w는 벽근처, e는 코어유동층을 나타낸다. 사용되는 산화물질의 질량분율은 CEA 프로그램을 이용하여 계산하였다[14].

여기서

$$ \text{Ω}=\frac{\alpha }{c_p}\frac{\overline{R}{T}_w}{P_w}\frac{e^{\frac{E_c}{\overline{R}{T}_w}}}{K_C{M}_C\text{'}} B_m=\left(X_{H_{2}O}+X_{C{O}_2}\right)\frac{M_C}{M_{mix,e}}$$

3.2 탄소/페놀릭의 열전도 유한차분법

탄소/페놀릭의 열분해를 고려한 비정상 열전도 지배방정식인 Eq. 14를 유한차분법을 이용한 해석을 수행하기 위해 이산화(discretization)할 수 있다. 시간은 후진차분, 공간은 중앙차분을 이용한 BTCS 방법을 기반으로 이산화하면, Eq. 19로 정리된다. 여기서 i는 유한차분법 내부격자 절점번호, n은 시간단계, r_s은 내열재 표면에서의 반경, ∆r은 내열재료의 격자크기, a는 i번째 절점에서의 반경, ∆t은 시간단계크기, ε은 공극률이다. Fig. 6은 노즐에서 탄소/페놀릭 열전도 유한차분법 격자를 나타낸다.

여기서

$$ \begin{array}{c}a=r_s+\left(i-1\right)∆r,\hfill \\ A_i^n=∆t{\lambda }_{i+\frac{1}{2}}^n\left(1+\frac{\mathrm{\Delta }r}{2a}\right)+0.5c_{G,d}{G}_{di}^n\mathrm{\Delta }t\mathrm{\Delta }r,\hfill \\ \begin{array}{c}{C}_i^n=\mathrm{\Delta }{r}^2{\rho }_0\left(1-{\epsilon }_{mi}^n\right)c_{mi}^n+\mathrm{\Delta }t\mathrm{\Delta }{r}^2\chi \left(1-\frac{{\rho }_{g,d}^A}{{\rho }_0}\right)c_{G,d}{D}_{mi}^n+∆t{\lambda }_{i+\frac{1}{2}}^n\left(1+\frac{\mathrm{\Delta }r}{2a}\right)+∆t{\lambda }_{i-\frac{1}{2}}^n\left(1-\frac{\mathrm{\Delta }r}{2a}\right)\hfill \\ B_i^n=∆t{\lambda }_{i-\frac{1}{2}}^n\left(1-\frac{\mathrm{\Delta }r}{2a}\right)-0.5c_{G,d}{G}_{di}^n\mathrm{\Delta }t\mathrm{\Delta }r,\hfill \\ F_i^n=\mathrm{\Delta }{r}^2{\rho }_0\left(1-{\epsilon }_{mi}^n\right)c_{mi}^n{T}_i^n-\mathrm{\Delta }t\mathrm{\Delta }{r}^2{D}_{mi}^n{Q}_{B,d}\hfill \end{array}\hfill \end{array}$$

Fig. 6

Grids in solid rocket nozzle.

유한차분법을 이용하여 온도분포를 구한 다음 시간단계 내에서 발생한 삭마두께를 구해야한다. 해당 시간단계에서 발생한 삭마두께는 Eq. 20과 같이 근사화하여 계산을 수행하였다. 발생한 삭마두께를 고려하여 다음 시간단계에서 격자크기를 재구성하여 계산을 수행하였다.