발간 현황

Fig. 1은 고온의 연소가스에 노출되었을 때 탄소/탄소와 같이 숯이 생성되지 않는 재료와 탄소/페놀릭과 같이 숯이 생성되는 재료에 대한 열반응 메카니즘을 보여준다[3]. 삭마 현상은 연소가스와 내열재 표면 탄소(carbon) 간의 산화 반응으로 인해 반응질량이 표면에서 소모되면서 유입되는 많은 외부 열에너지를 소멸시키는 열과 물질의 전달과정이다. 고체 로켓 성능을 높이기 위해 추진제에 알루미늄이 포함된 경우에는 연소과정에 4Al+3O₂→2Al₂O₃ 반응으로 인해 연소가스 내에 산화알루미늄(Al₂O₃) 입자가 존재하게 되고 Al₂O₃ 입자가 노즐 표면에 충돌하여 국부적으로 침식(erosion)되는 현상이 발생하며, 침식의 정도는 Al₂O₃의 입자 크기, 함량 및 분포도에 영향을 받는다[4].

Fig. 1 
Charring and non-charring ablation mechanisms [3].

산화재 분말(AP 등)을 고분자질 바인더(HTPB 등)와 혼합하여 제조한 복합형 추진제(composite propellant)가 고체 추진기관에 주로 사용되며, 추진제 연소가스 중에서 OH, O, O₂는 탄소와 아주 민감하게 반응하지만 농도가 매우 낮고, H₂와 N₂는 3,200K 이상의 높은 온도에서 탄소와 반응하므로 삭마를 야기 시키는 산화 반응에는 H₂O와 CO₂가 주된 반응물로 작용한다[5]. 내열재 표면에서 Eq. 1의 비균질 화학반응(heterogeneous chemical reaction)을 통해 삭마가 진행된다.

Eq. 1의 비균질 화학반응을 자세히 관찰하면 다음의 단계로 진행된다[5]. 즉, ① 단계 : 연소가스 유동장에 있는 H₂O와 CO₂가 내열재 표면으로 확산된다. ② 단계 : 내열재 표면의 온도가 상승하여 탄소 입자가 활성화되고 유동장으로부터 확산되어 들어 온 H₂O와 CO₂가 탄소 입자와 Eq. 1의 반응을 일으킨다. ③ 단계 : 생성된 CO 기체가 내열재 표면을 이탈할 상태가 된다. ④ 단계 : CO와 H₂는 연소가스 유동장으로 확산되고 또한 연소가스 유동장에 있던 H₂O와 CO₂가 내열재 표면으로 계속 확산된다.

위의 4가지 단계 중에서 ① 단계와 ④ 단계가 가장 느리게 진행되면 (즉, ② 단계와 ③ 단계가 더 빠르게 진행되면) 이 반응은 확산이 지배하고 그 반대이면 화학반응이 전체 반응을 지배한다. 후자의 경우를 화학반응 지배(kinetic-controlled) 영역이라 하고, 전자의 경우를 확산지배(diffusion-controlled) 영역이라 하며, 그 천이온도는 약 1,500K 이다[6]. 화학반응 지배영역은 내열재 표면 온도가 낮을수록, 연소가스 유동장의 압력이 낮을수록 또 연소가스 유속이 낮을수록 현저히 나타나고 이 경우 삭마는 내열재 표면 온도에 크게 영향을 받으며, 화학 반응률은 재료 표면에서 화학반응을 나타내는 Arrhenius 모델에 의해 결정된다. 확산지배 영역에서 화학 반응은 경계층으로부터 재료 표면으로의 산화물질 확산율에 의해 결정되고 질량 손실률(삭마율)은 연소가스내 H₂O, CO₂의 몰분율에 비례하며, 내열재 표면 온도에는 거의 무관하다[7,8].

Fig. 1(b)는 페놀릭 수지(binder)의 열분해로 인해 숯층이 생성되는 탄소/페놀릭과 같은 내열재 단면에서의 표면과 내부 열반응 개략도를 보여준다. 탄소, 수소, 산소의 고분자 화합물(C_lH_mO_n)인 페놀릭 수지의 열분해로 열분해층(pyrolysis zone)과 숯층(char zone)이 존재한다. 탄소/페놀릭은 표면온도가 373K가 되면 내열재에 함유된 수분이 증발하고, 온도가 573~600K로 더 상승하면 탄소, 수소, 산소의 고분자 화학물인 페놀릭 바인더가 열분해 반응을 시작하면서 약 40~45%는 열분해 가스(H₂ 27%, CO 6%, H₂O 6%, 기타)로 변하고 나머지는 탄소 성분의 숯층을 형성하며 삭마과정이 진행함에 따라 숯층의 두께가 증가하게 된다. 1,000K에 도달하면 페놀 수지의 열분해 반응은 거의 완료되고 복합재 표면에 노출된 숯 잔류물인 탄소는 노즐 유동장 내의 H₂O와 CO₂ 성분과 산화반응으로 삭마가 진행되기 시작한다. 따라서 페놀 수지는 열분해로 숯 잔여물과 열분해 가스로 변하며, 열분해가 완료되면 페놀릭 질량의 55~60% 정도가 숯으로 남는다. 열분해는 흡열 반응으로 내열재로 전달되는 열의 상당 부분을 흡수하고 처녀(virgin) 층으로 전달되는 열을 감소시킨다. 또한 열분해 가스는 표면의 경계층으로 유출되면서 내열재 표면으로의 열전달을 감소시키는 역할을 한다.

내열재 밀도는 열분해가 진행되면서 감소하며, 밀도 변화율은 Eq. 2와 같이 Arrhenius 식으로 표시된다[9].

여기서, k_m은 열분해 반응 비례상수, E_m은 반응 활성화 에너지, R¯는 일반기체상수, T는 현재 온도, n은 반응차수, ρ는 내열재 현재 밀도, ρ_c는 숯층 밀도, ρ₀는 내열재 처녀층 밀도이다.

Fig. 2는 실리카(S_iO₂) 섬유에 페놀릭(C_lH_mO_n) 수지를 함침하여 경화시킨 실리카/페놀릭 내열재 단면에서의 내부 열반응 개략도를 보여주며, 탄소/페놀릭과는 달리 표면에 실리카 용융층(melting layer)이 추가되어 4개 층이 존재한다[10]. 표면에 가까울수록 온도(T)와 열분해 가스 질량유속(G)은 증가하고, 밀도(ρ_m)는 감소하며, 공극(porosity)은 많아진다. 용융층에는 공극이 없어지고 밀도가 불연속적으로 증가한다.

Fig. 2. 
Schematic of in-depth thermal response of silica/phenolic composite [10].

실리카/페놀릭의 경우 상온에서 1,000K 까지는 내부 열반응이 탄소/페놀릭과 동일하지만, 온도가 1,450~1,500K로 상승하면 실리카 섬유가 연화(softening)되면서 실리카와 탄소(숯)가 흡열반응(C+S_iO₂→S_iO+CO)을 시작하여 CO와 SiO의 반응 가스를 생성한다. 1,900~2,000K로 온도가 더 상승하면 실리카가 완전히 용융되고 용융 상태에서 탄소-실리카 반응은 계속된다. 용융 잠열로 인하여 많은 열을 소모하며 탄화 및 탄소-실리카 반응에 의해서 형성된 공극들은 용융 실리카로 채워져 용융층의 체적이 수축되고 이에 따라 표면이 후퇴하게 되며, 용융 실리카는 연소가스 유동에 의해 전단(shear)이 되기도 하며 일부는 숯층의 탄소 입자를 덮어 탄소성분이 산화되는 것을 방지해 주는 작용도 한다. 온도가 더욱 상승하여 2,400~2,500K 정도가 되면 용융 실리카는 증발을 시작한다[11].

노즐 내부의 압축성, 난류 열·유동장 지배방정식을 유한체적법 등의 CFD 프로그램으로 풀어서 노즐 벽면의 대류 열전달계수를 구할 수 있다. 또한 반 실험식을 이용하는 방법으로는 원관내 완전히 발달된 비압축성 난류유동에서 압축성 효과 보정을 위해 기준온도(벽온도와 경계층 온도의 평균치)을 적용하고, 경계층을 가로지르는 물성치 변화 등을 고려한 Eq. 3의 Bartz 식으로부터 대류 열전달계수를 구할 수 있다[12].

여기서, D_t는 노즐목 직경, Pr은 Prandtl 수, p_c는 연소실 압력, c^*는 특성속도, r_c는 노즐목 곡률 반경, A_t는 노즐목 단면적, A는 노즐 임의 지점 단면적, σ는 경계층을 가로지르는 물성치 변화에 대한 보정계수를 나타낸다. Eq. 3에서 [ ]내 수식은 노즐 목에서 대류 열전달계수이며, 이는 연소실 압력(p_c)의 0.8승에 비례하고, 노즐목 직경(D_t)의 0.2승에 반비례함을 알 수 있다.

노즐 벽면의 대류 열전달계수를 구하는 또 다른 방법으로는 경계층 내에서 밀도가 변하고 경계층 두께가 아주 얇으며 벽면을 통한 물질전달이 있는 축대칭 난류경계층 유동에 대한 운동량 적분 방정식과 벽면이 가열되는 경우의 에너지 적분 방정식을 푸는 것이 있다. 이때 압축성 효과를 고려한 표면마찰계수(skin-friction coefficient) 및 Stanton 수 뿐 아니라 난류 경계층 내에서의 속도분포, 온도분포 및 점성계수에 대한 관계식을 추가로 정의해 주어야 한다[13]. 운동량과 에너지 적분 방정식을 풀어서 노즐 벽면의 대류 열전달계수를 구하는 프로그램으로는 BLIMP(Boundary Layer Integral Matrix Procedure), BLIMP-J, ARGEIBL(Aerotherm Real Gas Energy Integral Boundary Layer) 등이 있다[14].

내열재 열반응에 대한 에너지 방정식은 Eq. 4와 같이 표시되며, 에너지 보존 법칙에 따라 재료내부에 축적되는 열량(좌측항)은 열전도에 의해 재료에 유입되는 열량(우측 첫째항), 수지가 분해되어 분출되는 가스에 전달되는 열량(우측 두번째항) 및 재료내 수지가 반응하는데 소모되는 열량(우측 세번째항)의 합으로 표시된다[15].

여기서, m˙g는 분해가스의 질량유속, Δh_g는 생성된 가스의 열분해 열량이다.

고속 유동인 노즐 내부는 열역학적 상태의 조건이 변함에도 불구하고 화학반응이 진행되지 않고 화학조성이 일정하게 유지되는 동결 유동(frozen flow)으로, 속도가 낮은 연소실 내부는 평형 유동(equilibrium flow)으로 보통 고려한다.

Fig. 3의 노즐 내면에서 에너지 평형식은 Eq. 5와 같으며, 우측 두번째 항은 벽면에서 화학반응으로 인해 흡수된 열량을, 우측 세번째와 마지막 항은 복사로 인한 열량 및 Al₂O₃ 입자 충돌 등으로 인해 흡수하는 열량을 나타낸다. 그리고 노즐 외면에서의 에너지 평형식은 Eq. 6과 같으며, 우측 마지막 항은 외부의 열(노즐 연소가스 등)로 인해 흡수하는 열량 등을 나타낸다[16].

Fig. 3. 
Convective and radiant fluxes on rocket nozzle surfaces.

여기서, V_n는 삭마율, Q_s는 연소가스와 벽면의 화학반응에 의한 열적 효과(thermal effect)로써 온도의 함수이며, 온도가 증가할수록 선형적으로 감소하는 경향을 나타낸다.

내열재의 내부 열반응 모델링은 Table 1과 같이 3가지로 나눌 수 있다[17]. Level 1에서는 내열재의 열분해 반응을 무시하고, 열분해로 인해 흡수된 열량은 상변화 문제에서 잠열 효과와 비슷하게 비열 값을 조정함으로써 고려한다. 계산되는 변수는 온도뿐이다. Level 2에서는 열분해로 인한 질량 손실, 이로 인한 열흡수 그리고 표면으로 유출되는 분해가스의 대류냉각 효과를 고려한다. 그러나 생성된 분해가스는 재료 내에 저장되지 않고 순식간에 가열면을 향해 유출된다고 가정한다. 또한 분해가스의 밀도와 속도를 결합하여 질량유속(mass flow rate)이라는 한 개의 변수로 고려한다. 계산되는 변수는 내열재의 내부 밀도와 온도 및 분해가스의 질량유속이다. Level 3에서는 숯(기공)을 통과하는 분해가스 유동과 관련되는 국소 밀도, 속도, 압력, 기공 등의 변수를 고려한다. 분해가스 속도와 압력의 상관 관계식(Darcy 법칙) 및 분해가스의 질량보존 법칙을 추가로 적용한다.

내열재 해석 In-house 프로그램은 국외 반출 통제로 입수가 어렵지만 학회 논문이나 보고서를 통해 제한적으로(적용 수식, 해석 과정/결과 등) 접근이 가능하며, Table 2는 미국의 내열재 열반응 해석용 In-house 프로그램을 보여준다. ACE[18]와 CMA[19]는 1960년대에 극초음속 재진입 비행체의 첨두부(nose cone) 삭마현상 해석을 위해 Aerotherm 사에서 개발한 프로그램으로, 2개 프로그램은 서로 연동이 되어 운용된다. ACE 프로그램은 CMA의 표면 에너지 평형식을 계산하는데 필요한 물성치 즉, 표면과 경계층 가장자리에서의 열·화학적 상태를 결정하며, CMA 프로그램은 표면에서 삭마율과 표면으로부터 내부의 시간에 따른 온도분포 등을 계산한다. FIAT 프로그램[20]은 CMA 프로그램의 문제점(에너지 방정식이 열분해와 삭마 방정식과 explicitly하게 연결되어 있어 사용자가 정의하는 시간 간격이나 격자수에 민감하고, 분해가스 질량유속 또는 삭마율이 아주 큰 경우에 계산 수렴 어려움)을 개선한 프로그램으로 NASA Ames 연구센터에서 개발하였다. CMA와 FIAT 프로그램은 1차원 비정상 열전달 방정식을 푸는 프로그램으로, 재진입 비행체의 대기권 첨두부와 같이 2차원 특성이 크게 나타나는 곳에서는 적용할 수 없어 Aerotherm사에서 1980년대에 ASCC 프로그램[21]을 개발하였다. ASCC 프로그램은 삭마재와 주변 기체를 연계하여 2차원(축대칭 또는 평면)으로 해석을 수행하지만 삭마재 열분해로 인한 분해가스 영향은 고려하지 않았다. TITAN [22]은 NASA Ames 연구센터에서 2000년경에 개발한 프로그램으로, 분해가스의 영향을 고려하면서 2차원으로 해석하며, 삭마로 인한 격자 이동을 고려하기 위해 경계밀착좌표계(boundary- fitted coordinate)를 사용하였고 유한체적법으로 열분해가 고려된 에너지 방정식을 차분화하였다. 3dFIAT 프로그램[23]은 비행체 첨두부의 다양한 재진입 각도에 따른 영향 해석이 가능 하도록 FIAT 프로그램을 3차원으로 확장한 프로그램으로, 삭마와 열분해 해석이 가능하며 복잡한 형상이나 여러 삭마재로 구성된 구조물의 해석이 가능하도록 다중 블록 이동 격자계(multiple-block moving grid system)를 적용하였다.

Hero[24]는 미국 ATK사에서 개발한 2차원(축대칭) 노즐 해석 프로그램으로, 표면 열화학 반응(삭마), 내부 열분해 반응(숯) 및 열전달을 유한요소법으로 해석한다. 열해석과 동일한 격자로 구조해석을 수행할 수 있으며, parallel processing 적용으로 계산시간을 크게 향상시키었다. CHAR 프로그램[25]은 NASA Lyndon B. Johnson 우주센터에서 개발한 프로그램으로, 1, 2, 3차원 열전달, 삭마 및 열분해 반응을 해석하기 위해 비정렬 격자 유한요소법을 사용하였다. 온도나 표면 하중으로부터 야기되는 내부 응력은 선형 열탄성 해석으로 수행한다.

SEP(Société Européenne de Propulsion)사는 1992년경에 Marc 상용코드를 이용하여 Ariane 5 로켓 부스터 노즐의 열·구조해석을 2차원 축대칭 모델링으로 수행하였으며, 노즐 전방영역의 유한요소 격자 형상은 Fig. 5와 같다[26]. 노즐벽면의 대류열전달계수는 2차원 압축성 경계층 유동에 대해 유한체적법으로 Navier-Stokes 방정식을 푸는 ONERA-CERT에서 개발한 프로그램으로 구하였고, 노즐벽면 탄소와 연소가스내 3개 화학종(H₂O, CO₂, H₂) 간의 화학반응을 고려한 체 Arrhenius 식으로 노즐 삭마율을 계산하였다. 해석의 단순화를 위해 내열재의 내부 열반응을 무시한 열전도 방정식으로 노즐 벽면 온도를 계산하였다. 이는 노즐 목 부근의 탄소/탄소 내열재에는 합당하지만 페놀릭 수지의 열분해 반응이 동반되는 탄소/페놀릭이나 실리카/페놀릭 해석에는 부적절하다. 삭마로 인한 이동면을 고려하기 위해 표면에 인접한 영역에는 이동 격자(moving node)를, 내부 영역에는 고정 격자를 사용하였고, 이동 격자점에서의 온도는 이전 시간 단계에서 계산된 온도를 이용하여 보간법으로 구하였다. 구조해석에서 직교 탄성 비선형 거동 모델을 적용했고, 응력-변형률 관계는 온도의 함수로 고려했으며, 벽면 온도구배나 노즐내부 압력으로 인해 야기되는 내열재 변형률(strain)을 계산하였다.

Fig. 5 
FEM mesh of Ariane 5 rocket booster nozzle [26].

2002년경에 Snecma사(구 SEP사)에서는 탄소/페놀릭과 같이 내부 열반응이 존재하는 내열재의 열해석이 가능하도록 Marc 코드를 개발한 MSC사와 공동으로 기존 Marc 코드를 upgrade한 MSC.Marc-ATAS 상용코드를 개발하였다[17]. Fig. 6은 표면에서의 열평형 관계를 보여주며, 대류/복사 열전달, 물질전달(diffusion), Al₂O₃ 입자 충돌, 내열재 열분해, 분해가스 영향 등이 고려되는 것을 볼 수 있다. Fig. 7은 Snecma사의 in-house 코드와 MSC.Marc-ATAS 상용코드로 각각 계산한 내열재 밀도변화를 보여주며, 연소시간(1,2, ... 40초) 별로 유사한 결과를 얻었다.

Fig. 6 
Surface energy balance for MSC.Marc-ATAS [17].

Fig. 7 
Comparison of density field between Snecma in-house code and MSC.Marc-ATAS [17].

2008년경에 Snecma사와 MSC사는 공동으로 MSC.Marc-ATAS 코드에 열응력 해석 모듈을 추가하여 열전달과 열응력 해석을 연동으로 수행하였다[27]. 모든 물성치는 온도의 함수로 고려했으며, 구조 해석에서 삭마재 열분해에 따른 물성 변화를 고려하기 위해 열분해나 숯의 정도가 반영된 유효 열전도도(effective conductivity), 유효 Young‘s modulus, 유효 열팽창계수 등을 사용하였다.

ATK사는 2012년경에 Fig. 8에서 보듯이 Fluent 상용코드로 추진기관 CFD 해석을 통해 노즐벽면의 경계조건(대류열전달계수 등)을 구한 후 FEM Builder을 사용하여 내열재 열반응(삭마) 해석용 Hero In-house 코드로 데이터를 넘기며, Hero 코드에서는 연소시간 15초 동안 내열재 해석을 수행하여 노즐벽면의 삭마형상을 구한 후 이를 다시 Fluent 상용코드로 넘겨서 추진기관 CFD 해석을 수행한다. 이렇게 매 15초 마다 Fluent와 Hero 코드 간에 FEM Builder를 통해 데이터를 이동하면서 155초까지 추진기관 내부 유동과 내열재 열반응 해석을 연동으로 수행하였다[24]. Fluent 상용코드를 이용한 CFD 해석은 2차원 축대칭 2상 유동 비반응 정상상태로, k-Ω SST 난류모델 (y⁺= 30~100) 적용, Al₂O₃ 입자는 34.6wt% (24%는 40~300μm의 큰 응축입자, 76%는 미세 입자) 적용, ATK의 입자 breakup 모델을 적용하여 수행하였다. CFD 해석에서는 추진제 그레인 연소로 인한 연소 영역 변화로 야기되는 유동 특성을 고려하였으며, 내열재 삭마와 열분해 반응을 고려한 2차원 열해석을 수행하였다. 또한 유동(CFD)-열(Hero)-구조해석의 상호연동(coupling) 해석 과정을 제시하였다.

Fig. 8 
Rocket nozzle analysis using Fluent, Hero and FEM builder [24].

ATA Engineering 및 CUBRI사는 NASA 지원 하에 Fig. 9에서 보듯이 유동장 해석은 Loci/CHEM In-house 코드로, 삭마와 열전달 해석은 CHAR In-house 코드로, 구조해석은 ABAQUS 상용코드로 서로 연동하여 흑연(graphite) 노즐에 대해 유동-열반응-구조해석을 수행한 결과를 2017년에 발표하였다[28]. Loci/CHEM 코드로 구한 유동해석 결과인 압력, 열유속, 연소가스 화학종(H₂O, CO₂)의 질량 분율은 CHAR 코드로 넘기고, CHAR 코드의 해석 결과인 내열재의 벽면 온도, 생성된 화학종(CO, H₂) 질량유속, 삭마율은 Loci/CHEM 코드로, 내열재의 온도, 삭마율, 숯 분율은 ABAQUS 코드로 넘기는 것을 볼 수 있다. CFD 해석은 매 계산시간 단계(time step)마다 CHAR에서 계산된 삭마율을 반영한 표면 이동을 고려하여 CFD 격자를 재구성하였다. Fig. 10은 노즐 영역의 격자 형상을, Fig. 11은 노즐 내부 유동장의 마하수와 내열재의 온도분포를 보여준다. CHAR와 Loci/CHEM 코드 간에 노즐 벽면에서 격자점(node)들은 서로 일치하지 않으며, 노즐 유동장의 CFD 해석은 육각형(hexahedral) 격자를 사용했고, 벽면 인근 경계층 영역의 격자는 상당히 조밀하게 구성된 것을 볼 수 있다.

Fig. 9 
Information exchanged between Loci/CHEM, CHAR, and Abaqus [28].

Fig. 10 
Computational grid of nozzle region [28].

Fig. 11 
Mach number (fluid) and temperature (solid) contours at t=0.5 sec [28].

Fig. 12는 1차원 내열재 열반응(삭마/숯) 해석이 가능하도록 ABAQUS 상용코드에 DFLUX, UMATHT 등 5개의 서브루틴과 ALE remesh 알고리즘이 추가된 해석 흐름도를 보여준다[29]. 기존의 유한차분법(FDM) 또는 유한체적법(FVM) 기반의 CMA나 FIAT in-house 코드에서 사용되는 것과 동등한 열반응 수식(Table 2의 Level 2)이 유한요소법(FEM) 기반의 ABAQUS에 적용되었다. Fig. 13에서 보듯이 내열재 표면/내부 고정 지점에서 시간에 따른 온도변화는 ABAQUS와 FIAT 코드 간에 매우 유사하였으며, 이를 통해 내열재 열반응 해석용으로 ABAQUS 상용코드의 적용 가능함을 확인하였다. 추후 2차원, 3차원 열분해 삭마 모델을 ABAQUS 코드와 연계할 연구 계획이 있으며, 이것이 완성되면 상용코드의 활용도는 더욱 증대될 것으로 사료된다.

Fig. 12 
Schematic of ABAQUS analysis procedure for solving pyrolyzing ablation problems [29].

Fig. 13 
Comparison of temperature histories between ABAQUS and FIAT [29].

중국 Harbin 공대의 Meng 등[30]은 Fig. 14에서 보듯이 Fluent와 ABAQUS 상용코드를 연동하여 비행체 첨두부(nose cone)의 열반응(삭마/온도) 해석을 수행하였다. Fluent를 이용한 CFD 해석에서는 O, N, NO, O₂, N₂ 등의 5개 화학종 kinetic 모델을 사용했으며, ABAQUS를 이용한 열반응 해석에서는 첨두부 소재가 탄소/탄소라서 내부 열반응이 없는 열전도 방정식을 적용하였다. 삭마 표면 열반응에 대해서는 기존 연구자들의 연구[31]와 유사하게 표면의 탄소와 고속·고온의 주변 기체인 O, O₂, N 사이에 3개 화학반응(C+O→CO, 2C+O₂→2CO, C+N→CN)을 적용하였다. 삭마로 인한 이동 격자를 고려하기 위해 ABAQUS 코드에 UMESHMOTION 서브루틴을 적용하였다. Fig. 15는 쇄기(wedge) 모양의 첨두부 형상과 격자를 보여주며, 첨두부 전방 고온 기체 영역에서 충격파와 경계층으로 인한 복잡한 열·유동 특성을 잘 묘사하기 위해 벽면 부근의 기체 영역 격자는 조밀하게 구성된 것을 볼 수 있다.

Fig. 14 
Coupling scheme between Fluent and ABAQUS [30].

Fig. 15 
Computational grid of nose cone region [30].

국방과학연구소(ADD)는 연세대학교와 공동으로 Fluent 상용코드로 노즐 유동과 내열재 열반응 해석을 연동하여 수행할 수 있도록 Fig. 16과 같이 Fluent-2C 해석 모듈을 개발하였다[33,34]. ADD In-house 코드인 DAT-1C[32]에서 사용하던 열반응 수식(프로그램)을 Fluent 코드용 UDF(User Defined Function)로 만들었으며, 내열재 에너지 방정식인 Eq. 4의 우측 두 번째(분해가스에 전달 열량)와 세 번째항(수지의 흡열반응 열량)은 Fluent 코드의 표준 에너지 방정식에서 소스 항으로 처리하였고, 삭마로 인한 격자 이동을 모든 방정식에서 고려하였다. 내열재의 밀도 변화는 Eq. 2의 밀도변화율에 시간증분을 곱한 후 이전단계 밀도에서 빼서 구했고, 이때 기공(porosity)과 반응율 등 사용되는 값은 UDF 기능을 이용하여 각 격자마다 저장해 두었다가 사용하였다. 분해가스 질량유량은 표면으로부터 열분해 층이 위치하는 거리까지의 밀도 변화율을 적분하여 구하였다. 2차원 형태의 분해가스 거동을 묘사하기 위하여 다공성 매체내의 유체 속도가 압력 구배에 비례한다는 Darcy’s law를 적용하여 분해가스 보존방정식을 풀었다.

Fig. 16 
Schematic of Fluent-2C analysis procedure for solving pyrolyzing ablation problems [34].

Fig. 16에서 보듯이 BIM, S.Energy 등 10개의 UDF가 사용되었고, 각 계산 모듈 간에 다양한 변수들이 이동하며, 특수 기호의 의미는 다음과 같다. Ma는 마하수, k_c는 처녀층 밀도(ρ₀)에 대한 숯층 밀도(ρ_c)의 비율, m˙p는 분해가스의 질량유속, Q_d는 수지의 열분해 열량, G는 분해가스의 질량유량, V_s는 삭마속도, χ_H2O는 H₂O의 몰분율, k_s와 E_s는 표면 열반응의 Arrhenius 식에서 비례상수와 활성화 에너지를, 또한 하첨자 s와 m은 각각 내열재 표면과 내부를 나타낸다.

Fig. 17은 Fluent-2C를 이용한 내열재 열반응 해석 과정을 보여준다. 먼저 노즐 형상에 대해 수치해석 Geometry를 만든 후 격자를 생성하고, 노즐 내부 유동장 해석을 수행하여 노즐 벽면의 경계조건(대류열전달계수, 압력 등)을 구한다. 그 다음에 열반응 해석에 필요한 UDF를 정의하고 계산시간을 증가시키면서 비정상(unsteady) 내열재 열반응 해석을 수행한다.

Fig. 17 
Nozzle analysis using Fluent-2C [34].

Fig. 18은 토출관(blast tube)이 있는 노즐[8]에 대해 토출관 축 방향 가운데 위치에서 표면으로부터 깊이 방향으로 온도와 밀도 변화를 연소시간(3초, 6초, 10초)별로 Fluent-2C를 이용하여 해석한 결과와 DAT-1C 해석 결과를 보여준다. 그림에서 보듯이 온도와 밀도 변화는 2개 해석 결과가 거의 유사하였으며, 이를 통해 Fluent-2C의 열반응 해석 타당성을 입증할 수 있었다.

Fig. 18 
Comparison of temperature and density profiles between DAT-1C and Fluent-2C [34].