COMPUTATIONAL ANALYSIS OF AERODYNAMICS AND STATIC STABILITY OF RE-ENTRY VEHICLE AT LOW ALTITUDES

B.G. Park; S.J. Jo; S.M. Choi; S.J. Hwang; R.S. Myong

doi:10.6112/kscfe.2024.29.4.010

Preview

Original Article

Journal of Computational Fluids Engineering. 31 December 2024. 10-23
https://doi.org/10.6112/kscfe.2024.29.4.010

COMPUTATIONAL ANALYSIS OF AERODYNAMICS AND STATIC STABILITY OF RE-ENTRY VEHICLE AT LOW ALTITUDES

재진입 비행체의 저고도에서의 공력 및 정적 안정성에 관한 전산해석

B.G. Park¹

S.J. Jo¹

S.M. Choi¹

S.J. Hwang¹

R.S. Myong¹²^†

박 보근¹

조 성진¹

최 수민¹

황 선재¹

명 노신¹²^†

¹School of Mechanical and Aerospace Engineering, Gyeongsang National University

²Research Centrer for Aircraft Core Technology, Gyeongsang National University

¹경상국립대학교 기계항공우주공학부

²경상국립대학교 항공핵심기술선도연구센터

^{*Corresponding Author}

License (open-access, http://creativecommons.org/licenses/by-nc/4.0):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0) which permits unrestricted noncommercial use, distribution, and reproduction in anymedium, provided the original work is properly cited.

ABSTRACT

Re-entry vehicles may experience unstable flow during descent to Earth at low altitudes and the flow can negatively affect the stability of the vehicle. In this study, we investigate the characteristics of rear wake flow and frontal shock waves of re-entry vehicles for the angle of attack of 0º to -60º in low altitude transitional regime using a Navier-Stokes code. Five species and the transition SST turbulence model are used to describe the low-altitude flows at subsonic to hypersonic speed. It was shown that recirculation rings are getting bigger at lower altitudes and smaller with increasing angle of attack. By calculating the pitching moment at the center of gravity, the vehicle was shown to be statically stable at a low angle of attack. The angle of attack of neutral static stability decreases as the Mach number increases. Additionally, the aerodynamic center of the vehicle is calculated at each altitude and moves near the vehicle as the Mach number increases.

Keywords

Re-entry vehicle

Static stability

Transitional regime

High-speed flow

Computational fluid dynamics

키워드

재진입 비행체

정적 안정성

천이 영역

고속 유동

전산유체역학

MAIN

1. 서 론
2. 지배 방정식
3. 수치해석 기법
3.1 해석 기법
3.2 해석 형상 및 경계조건
3.3 해석 조건 및 화학종 모델
3.4 해석 기법 검증
4. 공력해석 결과
4.1 고도 별 받음각에 따른 후류 특성
4.2 Pitching Moment Coefficient를 통한 정적 안정성 계산
4.3 KRC 공력중심 위치
5. 결 론

1. 서 론

극초음속 비행과 우주 재진입은 최근 몇십 년 동안 꾸준히 연구되어 왔으며, 현재는 전 세계에서 우주 재진입 기술의 상용화를 위한 기술이 개발되고 있다[1]. 미국, 유럽 등 전 세계에서 지구 재진입 및 화성 대기 진입 비행체를 연구 및 개발 중에 있으며, 대표적으로 미국 나사의 Apollo, 유럽의 IXV, Space Rider, 일본의 OREX, Hayabusa, 중국의 Shenlong 등이 있다. 최근 몇 년 동안에는 미국의 민간 기업 SpaceX가 우주 비행체의 재사용을 위해 우주 재진입 기술을 개발하여 Dragon, Starship 등을 개발하여 재진입 기술의 상용화가 많이 이루어지고 있다. 최근 국내에서도 누리호의 발사가 성공적으로 이루어짐에 따라, 우주 발사체 뿐만 아니라 우주 재진입 기술을 개발하여 발사체에 탑재하여 기술의 검증하고자 하는 연구가 증가하는 추세이다. 대표적으로 한국항공우주연구원에서 재진입 비행체에 대한 연구[2]를 진행하였고, KAIST가 CubeSat의 재진입을 통해 재진입 기술을 위한 데이터를 수집하는 연구도 진행하였다. 이러한 재진입 기술에 대한 연구를 위해서는 실제 비행 시험 또는 풍동실험을 통해 수행할 수 있지만 이는 시험체 제작, 우주로 발사하는 시간적, 금전적 비용이 매우 크다. 또한 비행체에 장착 가능한 측정 장비의 한정되어 있어 데이터의 수급에 한계가 있다. 따라서 비용 및 공간적인 제약이 비교적 적고, 선행 연구들로 인해 정확성이 증명된 CFD(Computational Fluid Dynamics)를 많이 사용하고 있으며 이는 우주 재진입 비행체의 성능 해석에 핵심적인 역할을 하고 있다[3].

재진입 비행체는 행성 대기권 밖으로 발사되어 고도 100 km 이상에서 임무를 수행한 뒤 대기로 재진입 비행을 하게 된다. Fig. 1에서 나타나는 것처럼 비행체는 여러 대기층을 지나게 되며, 다양한 형태의 혼합 기체를 지나게 된다[4, 5]. 재진입 비행체는 초기에 마하수 약 40의 매우 빠른 속도로 진입을 시작하며 점점 속도가 감소한다. 고도 50~80 km의 고고도에서는 높은 마하수와 함께 충격파 이후 발생하는 다양한 이온화 반응으로 인한 진동 에너지와 병진 및 회전 에너지가 다른 비평형 상태와 매우 높은 열이 지배적인 환경이며[6, 7], 8], 9], 이 구간에서는 비행체가 열로부터 스스로를 보호할 수 있도록 열보호 시스템을 설계하여야 한다[10]. 또한 고도 30~40 km의 저고도에서는 5 이하의 비교적 낮은 마하수에 함께 비교적 밀도가 높아진 대기와의 마찰과 후류의 불안정한 유동으로 인해 비행체가 안정적인 비행이 어려워진다. 이 구간에서 비행체의 생존성을 확인하기 위하여 재진입 비행체 주변 공력계산 및 비행체의 안정성 분석이 필수적으로 이루어져야 한다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F1.jpg

Fig. 1.

Trajectories of various re-entry vehicles[4, 5]

극초음속 재진입 비행 및 재진입 비행체 에 대한 연구는 꾸준하게 이어져 왔다. 극초음속 및 초음속 단계에서 재진입 비행체에 대한 CFD 공력해석 연구는 전 세계에서 다양하게 진행되었다. Lee[11]는 2-온도 모델을 적용하여 뭉뚝한 형상에 대하여 화학 반응을 고려한 공력 가열 예측 연구를 하였고, Liever 등[12]은 CFD를 사용하여 Beagle 2 화성탐사 재진입 비행체에 대한 공력해석을 진행하였다. Pezzella 등[13]은 Stardust, Genesis, Hayabusa, Marco polo 재진입 비행체의 비행에서 발생하는 열의 차폐 설계에 대한 연구를 진행하였고, Desai 등[14]은 그 중 Genesis 태양풍 탐사 비행체에 대하여 지구 귀환 시나리오에 대한 연구를 진행하였다. 국내에서는 Mankodi 등[15]이 Stardust 재진입 비행체 주위 비평형 유동에 대하여 2-온도 모델과 NETT 모델의 비교 연구를 진행하였고, Kang[16] 등은 대기 재진입 우주비행체의 고온반응기체 효과에 대하여 전산해석을 통해 외부 열유동장을 비교하는 연구를 진행하였다. Lee 등[17]은 뭉뚝한 물체의 선두부 반경이 비평형 상태 흐름에서의 고온 특성에 관한 연구를 진행하였고, Chae 등[18]은 OREX, Stardust 재진입 비행체의 초기 지구 대기 진입 속도와 열적 특성과의 연관성에 관한 연구를 진행하였다.

한편, 저고도에서는 극초음속 대기 조건과는 다르게 높은 열 환경이 재진입 비행체에 미치는 영향성은 적어지며 비행체의 후류 등 주변 유동이 비행체에 미치는 영향성이 높아진다. 이러한 유동이 재진입 비행체의 안정성 및 궤도에도 영향을 미치게 되어, 결국 비행체의 생존성이 불안정해지므로 이 구간에 대한 분석이 이루어져야 한다. Mehta[19]는 5 이하의 낮은 마하수에서 ARD, OREX, Soyuz 등 여러 재진입 비행체의 후방 유동장을 분석하고 표면 압력계수를 통한 비행체 안정성 연구를 하였다. Komurasaki 등[20]은 RANS(Reynolds Average Navier-Stokes) 방정식을 통한 소성 표면을 가지는 재진입 비행체에 대한 저고도에서 발생하는 저 레이놀즈수 효과로 인한 Transition 현상의 공력해석을 진행하는 연구를 하였다. Shiva 등[21]은 FIRE2와 OREX 재진입 비행체의 저 마하수 조건에서 다양한 받음각에서의 유동 특성을 분석하는 연구를 하였다. Kang 등[22]은 IXV 우주 재진입 비행체에 대하여 수직, 수평방향 힘 계수와 피칭 모멘트 계수를 계산함으로써 안정성 및 조종성에 대한 연구를 진행하였다.

재진입 비행체에 대한 안정성을 확인하는 연구는 아직 미비하다. 특히 재진입 비행체의 공력중심을 계산하는 연구는 매우 적다. 따라서 본 연구에서는 새로운 재진입 비행체 형상의 재진입 궤도 중 저고도의 후방 유동 등 주변 유동 특성을 분석하고 고도에 따른 유동의 특성을 비교하였다. 또한 각 고도에서 받음각 별 피칭 모멘트 계수를 확인하고, 나아가 고도에 따른 공력중심 위치를 계산하여 비행체의 정적 안정성에 대하여 연구하였다.

2. 지배 방정식

고도 100 km 이상의 고고도에서는 DSMC와 같은 Mesoscale 동역학 기법을 사용하여 희박기체 효과를 고려한 열공력 해석을 수행한다. 이 고도에서는 기체의 이온화 및 해리 현상으로 인하여 연속체 이론을 바탕으로 하는 Euler 및 Navier-Stocks 방정식은 맞지 않다. 그러나 DSMC는 30~40 km의 비교적 저고도에서는 계산 비용이 매우 크게 요구되고, 해당 고도에서 Navier-Stokes 방정식을 근사적으로 적용 가능함을 이전 연구에서 확인하였다[2]. 연속체 이론을 판단하는 무차원 수인 Knudsen Number는 고도의 상승에 따라 선형적으로 증가하는데[23], Euler 나 Navier-Stokes 방정식은 Knudsen Number가 0.1 이하여야 기체의 연속성을 만족한다. 고도 100 km에서 Knudsen Number가 0.01로 나타나 고도 100 km 이하의 고도에서는 Navier-Stocks 방정식이 근사적으로 적용 가능함에 대한 타당성을 확보하였다.

이를 근거로 본 연구에서는 극초음속으로 재진입하는 비행체의 저고도에서 유동장 해석을 위하여 지배방정식으로 3차원 압축성 Naiver-Stokes 방정식으로 하였다. 이는 대기권의 공기 흐름을 모사하고, 점성과 압축성 효과를 고려할 수 있는 특징이 있다.

극초음속 유동에서 단일기체에 대한 압축성 Navier-Stokes 보존 방정식을 아래와 같이 표현한다.

(1)

\frac{\partial U}{\partial t} + \frac{\partial (E_{i} - E_{v})}{\partial x} + \frac{\partial (F_{i} - F_{v})}{\partial y} + \frac{\partial (G_{t} - G_{v})}{\partial z} = 0

위 식의 $U$ 는 보존 변수 벡터, $E$ , $F$ , $G$ 는 각각 x, y, z 방향 유동량(flux) 벡터, $i$ 는 비점성 유동량 벡터, $v$ 는 점성 유동량 벡터를 나타낸다. 보존 변수 벡터는 다음과 같이 표현한다.

(2)

U = [\begin{array}{l} ρ \\ ρ u \\ ρ v \\ ρ w \\ ρ e_{t} \end{array}]

위 벡터의 𝜌는 밀도, $u$ , $v$ , $w$ 는 각각 x, y, z 방향 속도 성분, $e_{t}$ 는 단위 질량당 총 에너지를 나타낸다. 비점성 및 점성 유동량 벡터의 x, y, z 방향 벡터는 다음과 같이 표현한다.

(3)

E_{i} = [\begin{array}{l} ρ u \\ ρ u^{2} + p \\ ρ u v \\ ρ u w \\ (ρ e_{t} + p) u \end{array}] G_{i} = [\begin{array}{l} ρ u \\ ρ u^{2} + p \\ ρ u v \\ ρ u w \\ (ρ e_{t} + p) u \end{array}] G_{i} = [\begin{array}{l} ρ u \\ ρ u^{2} + p \\ ρ u v \\ ρ u w \\ (ρ e_{t} + p) u \end{array}]

(4)

E_{v} = [\begin{array}{l} 0 \\ τ_{x x} \\ τ_{x y} \\ τ_{x z} \\ u τ_{x x} + v τ_{x y} + w τ_{x z} + q_{x} \end{array}] G_{v} = [\begin{array}{l} 0 \\ τ_{y x} \\ τ_{y y} \\ τ_{y z} \\ u τ_{y x} + v τ_{y y} + w τ_{y z} + q_{y} \end{array}] G_{v} = [\begin{array}{l} 0 \\ τ_{z x} \\ τ_{z y} \\ τ_{z z} \\ u τ_{z x} + v τ_{z y} + w τ_{z z} + q_{z} \end{array}]

여기서 $p$ 는 압력, $τ_{i j}$ 는 i,j 방향 점성 응력 성분, $q$ 는 각 방향으로 Fourier의 열전도 법칙을 사용하는 열전달 flux 성분을 나타낸다. 이러한 Navier-Stokes 보존 방정식은 x, y, z 방향으로의 유동량을 고려하여 유체의 상태를 시간에 따라 기술한다. Navier-Stokes 방정식을 사용하는 CFD에서는 이러한 형태의 방정식을 사용하여 복잡한 유체 흐름을 수치적으로 계산한다.

3. 수치해석 기법

3.1 해석 기법

본 연구에서는 초음속 영역인 마하수 1에서 5까지 변화시키며 진행하였다. 상용 CFD인 Ansys Fluent를 사용하여 전산해석을 수행하였다. 주요 수치 기법으로는 밀도 기반(Density-based), 시간 차분으로 내재적 시간 전진법(Implicit), 공간 차분을 위한 유한 체적법(Finite Volume Method)을 사용하였다. 수치해석 기법으로는 초음속 영역에서 충격파의 불연속성을 정확하게 포착하기 위하여 AUSM(Advection Upstream Splitting Method)을 사용하였다. Turbulence Closure를 위해 난류 경계층에서 Transition 영역을 고려할 수 있는 γ-Re_θ 모델로 알려진 Fluent에서 제공하는 Transition SST 모델을 사용하였다. Transition SST 모델은 2-방정식의 k-ω SST 난류모델을 바탕으로 하여 층류-난류 천이에 대한 정확한 예측을 위해 천이 레이놀즈 수 Re_θ와 간접 천이 변수 γ에 관련한 2개의 방정식을 추가적으로 계산하는데, Laminar 모델에 비해 재순환 유동 영역을 잘 묘사하는 것으로 알려져 있다. Michna[24]의 연구와 Aliaga[25]의 연구를 통해 Transition SST 난류모델의 타당성을 증명하였다.

3.2 해석 형상 및 경계조건

본 연구에서 사용한 해석 형상은 KRC(Korean Reentry Capsule)로, Fig. 2에 나타내었으며 JAXA의 Hayabusa Reentry Capsule을 기반으로 변형시켜 새롭게 설계한 모델이다. 비행체 지름은 400 mm, 높이는 200 mm, Fore-body Shell Angle은 34.6˚로 설계하였다. 비행체의 축 방향 길이를 D로 설정하여 격자의 경계는 비행체 전방 5D, 후방 15D, 상하 경계는 10D로 설정하였다. Fig. 3과 같이 Fluent를 사용하여 Polyhedral 격자로 약 110만 개의 격자를 생성하였다. Polyhedral 격자는 Tetrahedral 격자에 비해 격자 수를 최대 절반까지 절감하지만 높은 정확성을 보이고 유동의 세부적인 변화를 잘 포착하여 경계층의 해석에서 우수한 성능을 보여 재진입 비행체의 충격파 및 후방 유동을 해석하는데 적합하다. 해석의 정확성을 향상시키기 위하여 비행체 주변 유동 영역에 대하여 Body of Influence(BOI)를 설정하여 격자의 밀집도를 증가하였다. 비행체의 표면은 Adiabatic, Non-catalytic, Non-slip, Emissivity를 0.8의 고정된 값을 가지는 Wall로 하였다. 비행체 전방은 Velocity Inlet, 비행체 상, 하단, 후방은 Pressure Farfield의 경계조건을 설정하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F2.jpg

Fig. 2.

KRC configuration

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F3.jpg

Fig. 3.

Computational grid and boundary conditions

3.3 해석 조건 및 화학종 모델

해석에 사용된 비행 데이터는 Table 1에 나타내었다. 안정성을 확인하기 위해 받음각을 5º ~ -15º, -30º ~ -50º 범위에서 5º 간격으로 해석을 진행하였다. 비행체가 재진입 하는 저고도의 궤도에서 극초음속부터 아음속까지 마하수가 나타나는 고도를 선정하여 이 구간의 안정성을 확인하였다. 지구 대기를 모사하기 위하여 5가지의 화학종(N₂, O₂, NO, N, O)으로 구성된 혼합기체를 사용하였다. 고고도에서는 높은 온도로 인해 화학종끼리 반응하여 이온화되는 이온화 반응도 공력해석에 중요한 역할을 하지만, 비행체의 정체점에서의 온도가 9,000 K 이하에서는 이온화 반응에 의한 효과는 무시할 수 있는 수준으로 작은 것으로 나타났다[26]. 따라서 본 연구에서는 화학 반응은 고려하지 않고, 화학 종만 고려하여 진행하였다. 화학 종들의 질량 분율(Mass Fraction)은 Table 2에 나타냈다.

Table 1.

Flight conditions

Altitude(km)	Velocity(m/s)	Mach number	Temperature(K)	Pressure(Pa)	Reynolds number
34	1797.63	5.86	233.7435	663.41	234,580
32	987.49	3.26	228.4897	889.06	180,230
30	444.96	1.47	226.5091	1197.0	110,690
28	154.03	0.51	224.5272	1616.2	52,460

Table 2.

Species mass fraction

Species	N₂	O₂	NO	N	O
Mass fraction	0.747	0.25	0.001	0.001	0.001

3.4 해석 기법 검증

수치해석 기법에 대하여 검증하기 위해 극초음속 유동에서 지구 재진입 비행체 IXV 형상에 대한 전산해석을 진행하였다. 충격파 포착에 우수한 적응 격자를 고려하였으나, 비적응 일반격자로 구한 유동의 최대 마하수 및 온도, 벽면의 최대 압력 등 물리량이 적응 격자 결과와 3% 이내의 차이를 보여, 계산시간이 적게 소요되는 일반적인 격자를 사용하였다. 받음각 45도에서 마하수 3~11 에 대하여 유동 속도를 비교하였으며, 3D 형상으로 Pressure Farfield, Wall, Symmetry 경계조건을 사용하여 효과적으로 검증을 수행하였다. Normal Force Coefficient( $C_{N}$ ), Axial Force Coefficient( $C_{A}$ )의 전산해석 결과와 기존 형상에 대한 풍동실험 결과[27]와 비교하여 유사한 결과를 Fig. 4를 통해 확인하였다. 풍동 실험과 CFD 해석 결과 값의 최대 오차는 약 2% 이며 이를 통해 해석 기법에 대한 신뢰성을 확보하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F4.jpg

Fig. 4.

Aerodynamic coefficient of IXV

4. 공력해석 결과

4.1 고도 별 받음각에 따른 후류 특성

재진입 비행에서 비행체 주변 저고도에서의 유동은 충격파와 함께 후류가 복잡하게 발달한다. Fig. 5에 나타나듯이 비행체 전방에서 유동은 충격파를 통해 높은 압력과 밀도, 온도를 가지며 속도가 느려진다. 느려진 유동은 비행체 벽면을 따라 이동해, 비행체의 어깨에서 형상의 급격한 변화로 인해 Prandtl-Meyer 확장이 발생하여 유속이 가속된다[28]. 아음속, 천음속 단계의 고도 28, 30 km에서는 유동이 비행체의 어깨 부분을 지나며 자유류 유동의 속도보다 마하수 약 0.12, 0.44가 가속된 반면 초음속, 극초음속 단계의 고도 32, 34 km에서는 유동이 비행체의 어깨를 지나며 가속되지만 자유류 유동의 속도보다는 가속되지 않았다. 이는 높은 마하수에서 발생한 충격파로 인해 발생한 에너지가 벽면에 더 많이 영향을 주어 Prandtl-Meyer 확장에서 자유류 유동의 속도만큼 충분히 가속을 하지 못한 것으로 판단된다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F5.jpg

Fig. 5.

Mach number contours and streamlines at AOA 0°

가속된 유동은 경계층을 따라 흘러 재압축되며 분리 충격파를 발생시킨다. 이 때 경계층의 운동에너지가 물체에 발생하는 역압력구배를 극복하지 못하고 역류가 발생하면서 유동박리가 발생하여 비행체 후방에 재순환 영역이 발생하게 된다. 재순환 유동은 비행체 후방 부분을 따라 흐르며 다시 경계층에 영향을 미치게 되며 불안정한 흐름을 유발한다. 이 과정을 반복하게 되면서 비행체 후방 재순환 고리가 형성된다[29]. 재순환 고리의 크기는 비행체 주변 유동의 속도에 영향을 받으며, 속도가 빠를수록 고리가 작고 비행체에 가깝게 형성되고 속도가 느릴수록 고리가 크고 비행체와 멀리 형성되는 것을 확인하였다.

충격파가 나타나지 않은 고도 28 km에서는 유동이 비행체 전방 벽면과의 상호작용으로 유속이 느려졌다가 비행체 Shoulder에서 가속되며 Boundary Layer를 따라 흐르다 분리 충격파를 발생시키지는 않지만 Boundary Layer의 에너지가 역압력구배를 극복하지는 못하여 재순환 유동 및 고리가 발생하는 것을 확인하였다.

Fig. 6에는 받음각 -30º의 다양한 고도에서 유동을 나타냈다. 비행체 전방 충격파는 받음각에 따라 정체점에서 어깨 방향으로 이동하여 발생하였고, 이로 인해 비행체 양 어깨를 지나는 유속이 달라지게 된다. 유동이 흐르는 과정은 받음각 0º과 같으며 다만 비행체 양 어깨 주변 유동의 속도와 압력의 차이로 비행체 후방의 재순환 고리가 받음각 0º에 비해 더 작게 나타나는 것을 확인하였다. 또한 비행체 후방 분리 충격파 이후의 유동도 복잡한 형태를 띄게 되어 불안정한 유동의 형태를 확인하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F6.jpg

Fig. 6.

Mach number contours and streamlines at AOA -30°

Fig. 7에는 다양한 고도에서의 받음각 -60º의 유동을 나타냈다. 높은 절댓값의 받음각으로 인해 유동이 비행체에 가하는 압력이 어깨 부분에 집중되어 비행체가 회전하려는 힘이 평형 상태를 유지하려는 힘보다 강해져 비행체가 회전하게 된다. 이는 받음각에 따른 피칭 모멘트 계수 계산을 통해 확인하였다. 고도 28, 30 km에서는 상대적으로 낮은 속도로 인해 비행체 후방에 복잡한 재순환 유동이 발생한 반면, 고도 32, 34 km에서는 빠른 속도와 높은 받음각으로 인해 비행체 후방에 재순환 유동이 충분하게 발생하지 못한 것을 확인하였다. 이로 인해 비행체 전,후방의 압력 및 유동 균형이 이루어지지 않아 비행체가 매우 불안정하게 나타났다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F7.jpg

Fig. 7.

Mach number contours and streamlines at AOA -60°

Fig. 8은 고도 32 km의 받음각 0º, -30º, -60º 의 비행체 벽면 압력 분포를 나타낸 것이다. 압력 분포를 통해 비행체 전방의 충격파로 인한 압력 변화를 확인할 수 있다. 받음각 0º, -30º에 비해 –60º에서의 압력 분포가 비행체의 어깨에 집중되어 비행체가 회전하려는 힘이 증가하여 안정성에 영향을 미치게 된다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F8.jpg

Fig. 8.

Pressure distribution of KRC wall surface; (a) 0° (b) -30° (c) -60°

4.2 Pitching Moment Coefficient를 통한 정적 안정성 계산

재진입 비행체의 정적 안정성을 확인하기 위해서는 재진입 비행체의 중심축으로부터 수직 방향을 중심축으로 하는 모멘트인 피칭 모멘트를 계산하여야 한다. 정적 안정성이란 비행체가 힘의 균형을 유지하고 있는 조건에서 외부 교란을 받았을 때 원래의 평형 상태로 돌아가려는 복원력을 갖는 초기 경향성을 뜻하며, 이는 비행체의 무게중심에서 받음각에 따라 피칭 모멘트를 계산하여 받음각이 커질수록 피칭 모멘트가 감소하면 정적으로 안정하다고 판단한다. 받음각에 따른 피칭 모멘트 계수를 계산하는 방법은 다음과 같은 식으로 나타낼 수 있다.

(5)

L_{r e f} \frac{\partial C_{m}}{\partial α} = - C_{Y} \frac{\partial x_{c p}}{\partial α} + C_{X} \frac{\partial y_{c p}}{\partial α} + \frac{\partial C_{Y}}{\partial α} (x_{cog} - x_{c p}) - \frac{\partial C_{X}}{\partial α} (y_{cog} - y_{c p})

위의 식에서 $C_{m}$ 은 피칭 모멘트 계수, $x_{c p}, y_{c p}$ 는 비행체 압력중심의 x, y좌표, $x_{c o g}, y_{c o g}$ 는 비행체 무게중심의 x, y좌표, $C_{X}$ 는 Axial Force Coefficient, $C_{Y}$ 는 Normal Force Coefficient를 나타낸다[30]. 두 방향 힘의 계수는 다음과 같이 나타낸다.

(6)

C_{X} = C_{D} \cos α - C_{L} \sin α

(7)

C_{Y} = C_{L} \cos α + C_{D} \sin α

위의 식에서 $C_{L}$ 은 양력계수, $C_{D}$ 는 항력계수를 나타낸다. 위의 식들을 통하여 $\partial C_{m} / \partial α$ 가 음수일 때 비행체는 정적으로 안정하다. KRC 의 무게중심의 위치는 정체점으로부터 중심축을 따라 113.042 mm를 이동한 지점에 위치하며, 무게중심에서의 피칭 모멘트 계수를 계산한 결과를 Fig. 9에 나타내었다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F9.jpg

Fig. 9.

Pitching moment coefficient at AOA; (a) 0°~-15° (b) –30°~-60°

KRC 형상이 축대칭 형상이므로 받음각의 절댓값이 증가할수록 피칭 모멘트 계수가 감소하면 정적으로 안정하다고 판단하였다. Fig. 9에 따르면 받음각이 0º 에서 –15º 에서는 모든 고도에서 $\partial C_{m} / \partial α$ 가 음수로 나타나므로 정적으로 안정함을 확인하였다. 또한 아음속 조건인 고도 28 km를 제외하고 고도가 높을수록 $\partial C_{m} / \partial α$ 의 감소율이 더 큰 것으로 나타났다. 추가적으로 더 낮은 받음각까지 계산을 진행하여 비행체의 정적 안정성이 불안정해지는 받음각을 확인하기 위해 해석을 진행하여 Fig. 9에 나타내었다.

받음각 –30º 에서 –60º 의 피칭 모멘트 계수를 확인하였을 때 특정 받음각에서 $\partial C_{m} / \partial α$ 값이 음수에서 양수로 바뀌는 것을 확인할 수 있다. 이 상태를 중립 정적 안정성이라고 하여, $\partial C_{m} / \partial α$ 값이 양수가 되면서 정적으로 불안정해진다. 고도 28 km에서는 받음각 약 –50º에서 정적 중립안정성이 나타나며, 고도가 상승할수록 정적 중립안정성이 나타나는 받음각이 더 높아지며 고도 34 km에서는 약 –40º에서 정적 중립안정성이 나타남을 확인하였다. 이는 높은 고도일수록 자유류의 속도가 빠를 경우 항력이 증가하므로 Axial Force Coefficient와 Normal Force Coefficient가 증가하여 $\partial C_{m} / \partial α$ 값이 커지기 때문으로 판단하였다. 이를 통해 고도가 높을수록 정적 안정성을 만족하는 받음각의 범위가 감소하는 것으로 나타났다.

4.3 KRC 공력중심 위치

피칭 모멘트 계수 계산을 통해 재진입 비행체의 공력중심 위치를 확인하였다. 공력중심이란 받음각이 변하여도 피칭 모멘트의 값이 변하지 않는 지점을 뜻한다. 이를 $\partial C_{M a c} / \partial α = 0$ 로 나타낼 수 있다. 비행체의 안정성을 받음각에 따른 피칭 모멘트 계수의 변화율 뿐만 아니라 공력중심의 위치로도 확인할 수 있다. 비행체의 무게중심 위치가 공력중심의 위치보다 앞에 위치할 경우 안정하고, 뒤에 위치하면 불안정하다고 정의한다. KRC의 공력중심 위치를 Table 3과 Fig. 10을 통해 다음과 같이 나타냈다.

Table 3.

Position of aerodynamic center at various altitudes

Altitude(km)	28	30	32	34
Position(mm) (from stagnation point)	395.14	307.55	268.89	250.41

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290402/images/jkscfe_2024_294_010_F10.jpg

Fig. 10.

Aerodynamic center of KRC at various altitudes

일반적인 초음속 비행기 에어포일의 경우, 코드 길이가 에어포일의 두께보다 10배 이상 긴 형태를 띄므로 코드 길이의 1/2 위치에 공력중심이 형성된다. 그러나 본 캡슐형 재진입 비행체의 경우, 비행체의 길이보다 높이가 2배 짧은 형상을 가지기 때문에 비행체 후방 바깥에 공력중심이 위치한다. KRC의 공력중심 위치는 고도 28 km에서 정체점으로부터 250.41 mm, 고도 34 km에서 정체점으로부터 395.14 mm에 위치하며 약 144.73 mm 차이가 나는 것으로 나타났다. 고도가 높아질수록, 속도가 빨라질수록 비행체에 가까운 후방에 위치하는 것을 확인하였다. 이는 속도가 증가할수록 비행체 전방의 충격파로 인한 압력이 강해지며 비행체 후방 압력이 상대적으로 낮아지고, 이로 인해 재순환 영역이 비행체에 가까이 위치하며 영역의 중심점인 저압 영역이 앞으로 이동하며 공기역학적 힘의 중심점을 비행체에 가깝게 위치시키는 것으로 판단하였다. 공력중심과 무게중심의 위치의 차이를 정적 여유라고 하며, 이 정적 여유가 클수록 안정적이다고 판단한다. 이에 따라 고도가 낮을수록, 속도가 느릴수록 정적 여유가 크게 되므로 더 안정적이라고 판단하였다. 다만 모든 고도에서 무게중심의 위치보다 공력중심의 위치가 뒤에 위치하므로 저고도에서 비행체는 안정하다고 해석할 수 있다.

5. 결 론

저고도에서 새롭게 설계한 재진입 비행체 KRC의 정적 안정성을 확인하기 위하여 아음속부터 극초음속까지의 유동에 대하여 KRC 주변 유동장에 대한 전산해석을 수행하였다. 고도 별 여러 받음각에 따른 마하수 유동장 및 유선을 비교하여 비행체 주변 유동 흐름의 특성을 분석하고, 재진입 비행체의 무게중심에 가해지는 피칭 모멘트 계수를 계산하여 정적 안정성을 분석하였다. 또한 받음각에 따라 피칭 모멘트 계수가 변하지 않는 공력중심 지점의 계산을 통해 비행체의 안정성을 추가로 분석하였다.

저고도에서 발생하는 Transitional Regime을 계산하기 위하여 Transition SST 난류 모델을 사용하였고, 저고도에서 화학반응이 거의 일어나지 않는 것으로 판단하여 5개의 화학 종만 고려하여 해석을 진행하였다. 극초음속 영역인 고도 34 km, 초음속 영역인 고도 32, 30 km, 아음속 영역인 고도 28 km에 대하여 받음각이 변화할 때 유동 특성을 분석하였다. 고도가 낮아질수록 충격파가 재진입 비행체 멀리 발생하며 아음속 영역에서는 발생하지 않았다. 또한 고도가 낮아질수록 비행체 후방에서 발생하는 재순환 고리가 크게 형성되는 것을 확인하였고, 받음각의 절댓값이 높아질수록 재순환 고리가 작아지며 비행체 전방에서의 압력이 어깨 부분에 집중되며 비행체가 회전하려는 힘이 강해지는 것을 확인하였다. 이를 받음각에 따른 피칭 모멘트로 확인하였을 때, 낮은 받음각에서는 모든 고도에서 받음각의 절댓값이 커질수록 $\partial C_{m} / \partial α < 0$ 으로 나타나 KRC 형상이 정적으로 안정함을 확인하였다. 또한 정적으로 안정적이지 않아 지는 받음각인 정적 중립안정성을 나타내는 받음각의 절댓값이 고도가 높을수록, 속도가 빠를수록 더 낮아지는 것을 확인하였다. 추가적으로 공력중심의 위치를 계산하였을 때, 고도가 높아질수록 공력중심의 위치가 비행체에 더 가깝게 나타나 정적 여유가 더 감소하는 것을 확인하였지만, 무게중심의 위치보다 후방에 위치하여 비행체가 안정함을 확인하였다. 본 연구를 통해 저고도에서의 공력중심 및 재진입 비행체 정적 안정성을 확인하였으며, 시간에 따른 CFD 계산을 통해 동적 안정성을 계산하는 연구에 활용될 것으로 기대한다.

Acknowledgements

본 논문은 2024년 정부(과학기술정보통신부)의 재원으로 한국연구재단 스페이스챌린지사업(NRF-2022M1A3B 8075129)의 지원을 받아 수행된 연구입니다.

References

2005, Lips, T. and Fritsche, B., "A Comparison of Commonly Used Re-entry Analysis Tools," Acta Astronaut., Vol.57, No.2-8, pp.312-323.

10.1016/j.actaastro.2005.03.010

2008, Kim, C.W., Lee, Y.G. and Lee, D.S., "Aerodynamic Analysis of Sub-orbital Re-entry Vehicle," J. Comput. Fluids Eng., Vol.13, No.2, pp.1-7.

2020, Viviani, A., Aprovitola, A., Iuspa, L. and Pezzela, G., "Aeroshape Design of Reusable Re-entry Vehicles by Multidisciplinary Optimization and Computational Fluid Dynamics," Aerosp. Sci. Technol., Vol.105, 106029.

10.1016/j.ast.2020.106029

2014, Colombo, C., Letizia, F., Alessi, E.M. and Landgraf M., "End-of-life Earth Re-entry for Highly Elliptical Orbits : The Integral Mission," The 24th AAS/AIAA Space Flight Mechanics Meeting.

2010, Savino, R., Paterna, D., Fumo D.S. and D'Elia, M., "Plasma-Radiofrequency Interactions Around Atmospheric Re-Entry Vehicles : Modelling and Arc-Jet Simulation," The Open Aerospace Engineering Journal, Vol.3, pp.76-85.

10.2174/1874146001003010076

2020, Chae, J.H., Mankodi, T.K., Choi, S.M. and Myong, R.S., "Combined Effects of Thermal Non-equilibrium and Chemical Reactions on Hypersonic Air Flows Around An Orbital Reentry Vehicle," Int. J. Aeronaut. Space Sci., Vol.21, pp.612-626.

10.1007/s42405-019-00243-9

2020, Ejtehadi, O. and Myong, R.S., "A Modal Discontinuous Galerkin Method for Simulating Dusty and Granular Gas Flows in Thermal Non-equilibrium in the Eulerian Framework," J. Comput. Phys., pp.1-43.

2020, Mankodi, T.K., Bhandarkar, U.V. and Myong, R.S., "Collision Corss Sections and Nonequilibrium Viscosity Coefficients of N₂ and O₂ based on Molecular Dynamics," Phys. Fluids, Vol.32, 036102.

10.1063/1.5141148

2019, Mankodi, T.K. and Myong, R.S., "Quasi-classical Trajectory-based Non-equilibrium Chemical Reaction Models for Hypersonic Air Flows," Phys. Fluids, Vol.31, 106102.

10.1063/1.5119147

2023, Won, B., Park, Y., Han, N., Lee, E. and Shin, E., "A Study on Characteristics of Needle-Punching C/C in High-Enthalpy Plasma and High-Velocity Flame Environments," J. Korean Soc. Aeronaut. Space Sci., Vol.51, No.4, pp.235-242.

10.5139/JKSAS.2023.51.4.235

1998, Lee, M.G., "Aerodynamic Heating Prediction of Reacting Blunt Body Flow with an Impinging Shock Wave," International Heat Transfer Conference, Vol.4, pp.207-212.

10.1615/IHTC11.3690PMC124835

2003, Liever, P.A., Habchi, S.D., Burnell, S.I. and Lingard, J.S., "Computational Fluid Dynamics Predistion of the Beagle 2 Aerodynamic Database," J. Spacecr. Rockets, Vol.40, No.5, pp.632-638.

10.2514/2.6911

2012, Pezzella, G. and Trifoni, E., "Aerothermodynamic Design of an Ablative Heat Shield Design for a Sample Return Vehicle," 4th Sino-Italian Conference on Space Aerothermodynamics and Hot Structures, pp.52-64.

2001, Desai, P.N. and Cheatwood, F.M., "Entry Dispersion Analysis for the Genesis Sample Return Capsule," J. Spacecr. Rockets, Vol.38, No.3, pp.345-350.

10.2514/2.3707

2019, Mankodi, T.K. and Myong, R.S., "Investigation of Non-Equilibrium Reacting Air Flow Around Stardust Re-entry Vehicle," J. Korean Soc. Aeronaut. Space Sci., pp.20-21.

2014, Kang, E.J., Kim, J.Y., Park, J.H. and Myong, R.S., "Computational Investigation of the High Temperature Reacting Gas Effects on Re-entry Vehicle Flowfields," J. Comput. Fluids Eng., Vol.19, No.1, pp.7-14.

10.6112/kscfe.2014.19.1.007

2003, Lee, C.H. and Park, S.O., "Effects of Nose Radius of Blunt Body on Aerodynamic Heating in Thermochemical Nonequilibrium Flow," J. Comput. Fluids Eng., Vol.8, No.4, pp.34-40.

2019, Chae, J.H., Lee, J.W. and Myong, R.S., "Computational Thermal Analysis of Earth Reentry Vehicles for Different Flight Missions," J. Comput. Fluids Eng., Vol.2019, No.10, pp.34.

2017, Mehta, R.C., "Effect of Geometrical Parameters of Reentry Capusle over Flowfield at High Speed Flow," Adv. Aircr. Spacecr. Sci., Vol.4, No.4, pp.487-501.

2000, Komurasaki, K. and Candler, G.V., "Laminar-to-Turbulent Transitions over an Ablating Reentry Capsule," Acta Astronaut., Vol.47, No.10, pp.745-751.

10.1016/S0094-5765(00)00117-X

2012, Shiva, P.U. and Srinivas, G., "Flow Simulation over Re-Entry Bodies at Supersonic & Hypersonic Speeds," Int. J. Eng. Res. Dev., Vol.2, No.4, pp.29-34.

2014, Kang, E.J., Kim, J.Y., Chun, S.H., Myong R.S. and Cho T.H., "Investigation of the Stability and Controllability of a Re-entry Vehicle Using Compuational Aerodynamic Analysis," J. Comput. Fluids Eng., Vol.2012, No.5, pp.247-250.

2018, Fadgyas, M.C., Pricop, M.V., Niculescu, M.L., Cojocaru, M.G. and Dumitrache, A., "Semi-Empirical Relations for Pressure Distributions in Hypersonic Regime," International Conference of Numerical Analysis and Applied Mathematics 2017, 370007.

10.1063/1.5043975

2021, Michna, J., Rogowski, K., Bangga, G. and Hansen, M.O.L., "Accuracy of the Gamma Re-Theta Transition Model for Simulating the DU-91-W2-250 Airfoil at High Reynolds Numbers," Energies, Vol.14, No.24, p.8224.

10.3390/en14248224

2020, Aliaga, C., Guan, K., Selvanayagam, J., Stokes, J., Viti, V. and Menter F., "Hypersonic Applications of the Laminar-Turbulent Transition SST Model in ANSYS Fluent," AIAA Aviation 2020 Forum, p.3290.

10.2514/6.2020-3290

1959, Hansen, C.F., "Approximations for the Thermodynamic and Transport Properties of High-temperature Air," NASA Technical Note 4150.

2007, Paolo, B., Sylvian G., Patrice, P. and Jacques M., "Pre-X Atmospheric Re-entry Experimental Lifting Body: Program Status and System Synthesis," Acta Astronaut., Vol.61, pp.459-474.

10.1016/j.actaastro.2007.01.053

2004, Kim, H.D., Matsuo, S. and Setoguchi, T., "Passive Control Applied to Prandtl-Meyer Expansion Flow with Homogeneous Condensation," Aerosp. Sci. Technol., Vol.8, pp.369-376.

10.1016/j.ast.2004.02.002

2020, Mehta, R.C., "Computation of Base Pressure Based on to Fill-up the Growing Space Applied to Reentry Capsules," AIAA Aviation 2020 Forum.

10.2514/6.2020-2711

2009, Hirschel, E.H. and Weiland C., Selected Aerothermodynamic Design Problems of Hypersonic Flight Vehicles, American Institute of Aeronautics and Astronautics, Springer, pp.357-369.

Journal of Computational Fluids Engineering ISSN:1598-6071(Print) 3022-6252(Online) 한국전산유체공학회지

Preview

COMPUTATIONAL ANALYSIS OF AERODYNAMICS AND STATIC STABILITY OF RE-ENTRY VEHICLE AT LOW ALTITUDES

ABSTRACT

MAIN

Fig. 1.

Trajectories of various re-entry vehicles[4, 5]

(1)

(2)

(3)

(4)

Fig. 2.

KRC configuration

Fig. 3.

Computational grid and boundary conditions

Table 1.

Flight conditions

Table 2.

Species mass fraction

Fig. 4.

Aerodynamic coefficient of IXV

Fig. 5.

Mach number contours and streamlines at AOA 0°

Fig. 6.

Mach number contours and streamlines at AOA -30°

Fig. 7.

Mach number contours and streamlines at AOA -60°

Fig. 8.

Pressure distribution of KRC wall surface; (a) 0° (b) -30° (c) -60°

(5)

(6)

(7)

Fig. 9.

Pitching moment coefficient at AOA; (a) 0°~-15° (b) –30°~-60°

Table 3.

Position of aerodynamic center at various altitudes

Fig. 10.

Aerodynamic center of KRC at various altitudes

Acknowledgements

References