[ Research Paper ]

Journal of the Korean Society of Propulsion Engineers - Vol. 26, No. 5, pp.35-43

ISSN: 1226-6027 (Print) 2288-4548 (Online)

Print publication date 31 Oct 2022

Received 01 Aug 2022 Revised 08 Oct 2022 Accepted 15 Oct 2022

DOI: https://doi.org/10.6108/KSPE.2022.26.5.035

VOF와 FSI 방법을 적용한 CBT 연소 정밀 모델링 및 해석

강정석^a ; 박종근^a ; 성홍계^b^{, *}

CBT Combustion Precise Modeling and Analysis Using VOF and FSI Methods

Jeongseok Kang^a ; Jonggeun Park^a ; Hong-Gye Sung^b^{, *}

aDepartment of Aerospace and Mechanical Engineering, Korea Aerospace University, Korea
bSchool of Aerospace and Mechanical Engineering, Smart Air Mobility Engineering, Korea Aerospace University, Korea

Correspondence to: ^* E-mail: hgsung@kau.ac.kr

Copyright Ⓒ The Korean Society of Propulsion Engineers
This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License(http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

초록

본 연구에서는 고체추진제를 사용하는 closed bomb test(CBT)의 연소에 대한 정밀 모델링 및 해석을 수행하였다. 기상과 고상을 동시에 해석하기 위해 fluid structure interaction(FSI) 기법을 사용하였으며 기체상과 그레인의 연소해석은 Eulerian 방법을, 그레인의 이동은 Lagrangian 방법을 적용하였다. 고체상의 그레인과 연소가스의 상호 작용은 소스텀을 통해 완전 결합(fully coupled) 되도록 하였다. 그레인의 연소거리와 연소면의 이동을 모사하기 위하여 volume of fluid(VOF) 방법을 사용하였고, 그레인에 작용하는 힘은 그레인 연소면에 작용하는 압력과 중력을 고려하고, VOF의 속도항에 그레인 연소속도와 그레인 이동속도를 고려하였다. 개발한 수치모델을 바탕으로 1개와 3개 그레인에 대한 연소해석을 수행하여 실험결과와 비교 검증하였다. 연소시에 나타나는 압력 섭동에 대한 음향장을 분석하였다.

Abstract

Precise modeling and analysis of closed bomb test(CBT) combustion using solid propellants was performed. The fluid structure interaction(FSI) method was implemented to analyze the gas and solid phases at the same time. The Eulerian analysis method was applied for the gas phase and grain combustion, and the Lagrangian analysis method was implemented for the grain movement. The interaction between the solid phase grains and the combustion gas was fully coupled through the source term. The volume of fluid(VOF) method was used to simulate the burning distance of the grain and the movement of the combustion surface. The force acting on the grain was comprised of the pressure and gravity acting on the grain burning surface, and the grain burning rate and grain movement speed were considered in the velocity term of the VOF. The combustion analysis was performed for both one and three grains, and fairly compared with the experiments. The acoustic field during grain combustion due to pressure fluctuations was also analyzed.

Keywords:

Solid Propellant, Closed Bomb Test, Volume of Fluid, Fluid-Structure Interaction

키워드:

고체추진제, 밀폐용기실험, 유체체적법, 유체-고체 연성

1. 서 론

고체 추진제는 화포, 로켓, 총, 폭약 및 에어백 등 다양한 분야에서 사용되고 있다[1]. 특히 고체 추진제는 에너지 밀도가 높고 폭발위험이 적으며 탑재체로 장기간 보관이 용이하다. 이러한 고체 추진제는 사용 목적에 따라 다양한 구성성분을 조합해 연소속도를 조절하여 사용한다[2]. 고체추진제의 연소속도는 drop tube, strand burner, closed bomb test(CBT)와 같은 실험장치를 통해 측정하며 이는 고체추진기관의 성능에 가장 중요한 인자 중 하나이다[3,4]. 고체추진제의 연소속도는 곧 고체추진제 연소면의 이동속도를 의미한다. 따라서 고체추진기관의 성능을 예측하기 위해 고체추진제 그레인의 형상에 따른 연소면을 추적하는 연구들이 수행되고 있다.

그레인의 연소면 해석을 위해 초기에는 해석적 방법인 generalized coordinate(GC) 기법이 연구되었다[5]. 하지만 GC 기법은 해석형상이 복잡해지면 오차가 발생하는 한계점이 있다. 이후 Eulerian 해석방법인 level-set 해석 기법[6]과 volume of fluid(VOF) 해석기법[7]을 이용해 연구가 진행되고 있다. 이외에도 연소면과 격자점 사이의 최소 거리를 이용하는 minimum distance function(MDF) 기법도 연구되고 있다[8]. Eulerian 해석방법 외에도 Lagrangian 해석방법이 있으며 대표적으로 face offsetting method(FOM) 기법이 있다[9]. FOM 기법은 Lagrangian 요소를 적용하여 Eulerian 해석방법보다 강건하고 정확한 해석이 가능하다는 장점이 있다.

그레인의 연소면 해석을 위해 위와 같이 활발한 연구가 진행되고 있지만, 모두 고체로켓모터에 장착된 그레인 해석을 가정하였기 때문에 그레인의 움직임은 고려하지 않았다. 또한 그레인의 연소면 해석과 연소에 따른 고체추진기관 성능해석을 별개로 진행하여 연소시에 발생하는 기상과 고상 간의 양방향 물리적 현상이 결여되어 있다.

따라서 본 연구에서는 위와 같은 한계점들을 극복하고자 그레인의 burn-back뿐 아니라 그레인의 이동을 함께 모사하고 기상과 고상 간의 양방향 연성해석 수치 모델링을 수행하였다. 기상과 고상이 완전 결합(fully coupled)되어 있기 때문에 보다 물리적인 연소해석이 가능하다. 기상과 고상의 연성해석을 위해 fluid structure interaction(FSI) 방법을 적용하였으며 유체 거동과 그레인의 연소면은 Eulerian 해석기법으로, 그레인의 이동은 Lagrangian 해석기법을 적용하였다. 그레인 연소면과 연소면 이동은 VOF방법을 적용하였으며 CBT 내부의 난류 유동은 동적 Smagorinsky LES 난류모델을 이용하였다. 구성한 수치 모델을 바탕으로 단공 원통 그레인을 사용하는 CBT해석을 수행하여 실험결과와 비교 검증하였다.

2. 이 론

2.1 수치해석 흐름도

Fig. 1은 본 연구에서 구성한 수치해석의 전체 흐름도이다. 기상과 고상을 각각 Eulerian, Lagrangian 해석방법을 사용해 독립적으로 해석한다. Eulerian 해석방법에서 기상의 거동은 Navier-Stokes 방정식에 의해 해석하며 그레인의 연소면 거동은 VOF 방정식을 통해 계산한다. 이를 통해 그레인의 연소질량과 에너지를 계산하고 이를 소스항에 저장하여 시간전진 기법에 따라 변수값을 업데이트한다. 이후 Lagrangian 해석에서는 그레인의 이동을 해석한다. 계산의 효율성을 위해 병렬연산을 수행하고, 1개 그레인 해석은 100개 cpu를, 3개 그레인 해석은 200개 cpu를 사용하였다. 각 cpu에서 계산된 그레인에 작용하는 압력힘과 중력힘을 합산하여 그레인의 이동 속도를 계산한다. 각각의 그레인은 독립된 속도를 가지며 그레인 속도값을 VOF의 속도항에 대입하여 그레인의 다음 위치를 업데이트 한다. 이 과정은 그레인이 모두 연소될 때까지 반복하여 계산을 수행한다.

Fig. 1

Numerical analysis flow chart.

2.2 Eulerian 해석방법(기체상)

∂ ρ - ∂ t + ∂ ρ - u ~ i ∂ x i = m ˙ g r a i n V o l

(1)

∂ ρ - u ~ j ∂ t + ∂ ρ - u ~ i u ~ j + p - δ i j ∂ x i = ∂ τ ~ i j - τ ~ i j s g s ∂ x i

(2)

∂ ρ - E ~ ∂ t + ∂ ρ - E ~ + p u ~ i ∂ x i = ∂ ∂ x i u ~ j τ ~ i j + λ ∂ T ~ ∂ x i - H i s g s + σ j s g s + m ˙ g r a i n e g r a i n V o l

(3)

∂ ρ - Y ~ k ∂ t + ∂ ρ - Y ~ k u ~ i ∂ x i = ∂ ∂ x i - ρ - V ~ k, i Y ~ k - Φ k, i s g s - Θ k, i s g s + m ˙ g r a i n Y k V o l, k = 1, …, N

(4)

CBT 내부에서 발생하는 복잡한 난류유동을 모사하기 위해 큰 와류는 직접 풀고, 격자보다 작은 와류는 모델링하는 LES 난류 모델을 적용하였다. Favre-averaged된 질량, 운동량, 에너지, 화학종 보존방정식이 사용되었다. 미완결항인 SGS항은 다음과 같이 정의된다.

τ i j s g s = ρ - u i u j - u i ~ u j ~

(5)

H i s g s = ρ - E u i - E u i ~ + ρ u i - - p - u i ~

(6)

σ i s g s = u j τ i j - - u j ~ τ i j ~

(7)

Φ k, i s g s = ρ Y k u i - - ρ - Y k ~ u i ~

(8)

Θ i s g s = ρ Y k V i - - ρ - Y k ~ V k, i ~

(9)

미완결항은 난류유동을 정확히 모사하기 위해 동적 Smagorinsky 모델을 적용하였다[10]. 그리고 Eq. 1-4의 우측항에 있는 그레인 연소에 의한 소스항은 2.3절에 설명되어 있다.

2.3 그레인 연소면 해석방법

∂ α ∂ t + V → r e g r e s s i o n ⋅ ∇ α = 0

(10)

V → r e g r e s s i o n = r ˙ ⋅ n →

(11)

r ˙ = a P n

(12)

그레인의 연소면을 해석하기 위해 Eq. 10과 같이 VOF 방정식을 해석한다. 여기서 체적분율은 그레인의 체적분율을 의미한다. 그레인 이동거리 해석을 위해 그레인 표면 위치별로 압력에 따른 연소거리를 계산하고 해당 속도를 VOF 방정식의 속도항에 대입한다. 그레인 표면 속도벡터는 Eq. 11과 같이 각 격자 위치에서 그레인의 연소율과 표면 기울기 벡터와의 곱으로 표현된다. 기울기는 least square 방법을 이용해 계산한다. 그레인의 연소율은 압력에 비례하며 Eq. 12와 같이 표현된다[11]. 여기서 a와 n은 실험적으로 구해지는 상수이다.

m ˙ g r a i n = △ α ⋅ V o l ⋅ ρ g r a i n △ t

(13)

Eq. 10을 통해 계산된 체적분율의 변화량을 Eq. 13에 대입하여 시간당 그레인의 연소질량을 구할 수 있으며 이 값은 Eq. 1-4의 소스항에 사용된다.

2.4 Lagrangian 해석방법(고체상)

∂ α ∂ t + V → g r a i n ⋅ ∇ α = 0

(14)

∑ F = F g r a v i t y + F p r e s s u r e = m g r a i n a g r a i n

(15)

V g r a i n n + 1 = V g r a i n n + a g r a i n ⋅ △ t

(16)

본 연구에서는 기상과 고상을 함께 해석하기 위해 FSI 기법을 적용하였다. 기상과 고상의 경계는 그레인 체적분율 0.5를 기준으로 하였으며 그레인의 이동만을 고려하기 위해 고상은 강체로 가정하였다. 경계에서 그레인은 벽과 동일한 조건으로 설정하여 기상유동은 그레인을 통과하지 못한다. 기상유동은 2.2절과 같은 방법으로 해석되며 그레인의 이동을 해석하기 위해 각 그레인을 Lagrangian 방법을 통해 개별 해석한다.

그레인의 이동을 해석하기 위해서는 그레인에 작용하는 중력힘과 압력힘을 합산해야 한다. 압력힘을 구하기 위해 인접한 격자에서의 압력정보를 사용한다. Eq. 15-16을 통해 그레인의 가속도 및 속도를 업데이트하고 Eq. 14의 속도항에 대입한다. 각 그레인은 개별적인 속도를 가지며 Eq. 14를 통해 이동한다. 만약 그레인의 이동방향에 벽이나 다른 그레인이 존재하면 해당 그레인은 더 이상 이동하지 않는다.

2.5 수치기법

기상과 고상의 지배방정식은 유한체적법을 적용해 차분화 하였으며 시간차분화는 3차 Runge-Kutta 기법을 적용하였다[12]. 기상해석에서 대류 플럭스 기법은 HR-SLAU2[13]를 적용하였으며 점성 플럭스는 중심 차분을 사용하였다. 상태방정식은 연소실 내부의 고온, 고압을 정확히 모사하기 위해 Noble-Abel 상태방정식을 사용하였다[14]. 고상의 경우, 경계면에서 발생하는 수치적 확산을 막고 경계면을 정확히 예측할 수 있는 Tangent of Hyperbola for Interface Capturing(THINC) 재건설(Reconstruction)기법을 적용하였다[15].

3. 결 과

3.1 해석모델

본 연구에서는 단공 원통 그레인을 사용하는 CBT에 대한 수치해석을 수행하여 실험결과와 비교 검증하였다. 실험에서는 100개의 그레인을 사용하지만 효율적인 계산을 수행하기 위해 1개, 3개 그레인에 대한 해석을 수행하였다. 이 과정에서 챔버의 체적은 그레인 양의 비율과 같은 1/100, 3/100배로 조절하여 상사성을 유지하였다. 해석형상은 Fig. 2와 같으며 CBT 상단 중앙에 압력센서가 장착되어 있는 형태이다. 경계조건은 모두 벽면조건으로 설정하였다. 실험형상과 각 그레인별 챔버 크기 및 그레인 크기는 Table 1에 요약하였다. 해석조건과 그레인 연소가스 조성은 Table 2와 같다. 그레인 연소가스의 조성은 그레인 추진제의 조성을 기반으로 CEA 해석을 통해 추출하였다[16]. 그레인 연소속도에 대한 상수 a와 n값은 기존에 수행한 연소속도 추정 연구를 기반으로 실험데이터를 통해 구하였으며 Table 3에 요약하였다[14].

Fig. 2

Computational domain.

Fig. 3

Grain shape.

Table 1.

Computational geometry.

Table 2.

Operating conditions.

Table 3.

Burning rate constants.

3.2 해석결과

3.2.1 1개 그레인

Fig. 4는 1개 그레인의 시간에 대한 해석결과로 온도장과 그레인의 변화를 나타낸 그림이다. 그레인이 연소되면서 모든 공간으로 연소가스와 에너지가 방출되어 온도가 상승하는 것을 볼 수 있다. 특히 그레인이 단공형이기 때문에 그레인 중앙에서는 단공 공간을 따라 그레인 외부(CBT 공간)로 연소가스가 빠르게 분사된다(Fig. 4(a)). 그레인의 연소로 인해 CBT 내부 압력이 상승하고 상승한 압력은 그레인 표면의 연소속도를 가속한다. 따라서 그레인은 시간이 지남에 따라 더욱 빠르게 연소가 진행된다(Fig. 4).

Fig. 4

Temporal evolution of temperature field and grain surface for single-grain burning.

Fig. 5는 VOF 기법으로 계산된 시간에 따른 그레인의 질량변화를 해석적 결과와 비교한 그림이다. 해석적 방법은 기하학적 해석방법을 사용하였으며 그레인 연소에 따른 내경, 외경, 높이의 변화를 통해 그레인 질량을 도출하였다. 전반적으로 서로 잘 일치하며, 연소후반인 17 msec 이후에 다소 차이를 보인다. 이는 Eulerian 해석방법의 특성으로 그레인이 작아짐에 따라 격자가 부족하여 나타나는 현상이다. 격자를 작게 할수록 오차를 더 줄어들게 할 수 있으나 계산시간은 크게 증가한다. 현재의 격자를 사용한 계산 결과의 압력기울기 및 최대압력이 실험결과와 유사하므로 이러한 오차는 무시할 수 있다(Fig. 6).

Fig. 5

Grain mass vs. time history.

Fig. 6

Pressure time history with 2.5-4.5 msec zoom.

Fig. 6은 시간에 따른 CBT 내부압력을 실험결과와 함께 나타낸 그림이다. 실험결과와 수치해석결과는 같은 압력기울기를 보이며 최종압력은 0.3% 내의 차이를 보인다. 확대한 그래프를 보면압력진동이 관찰되며 이는 약 12.3 kHz이다. CBT는 길이방향에 대해 음향학적으로 닫혀진 계이므로 다음과 같이 해석적으로 주파수를 계산할 수 있다.

f = c 2 π λ m, n 2 R c 2 + l 2 π 2 L c 2

(17)

여기서 λ_m,n는 횡방향 모드의 고유값이다[17]. 음속은 CBT 내부의 평균 음속으로 2.5-5 msec에서는 650 m/sec이며 이를 통해 계산된 길이방향 모드(l=1, m=0, n=0)의 주파수는 11.7 kHz이다. 따라서 연소에서 나타나는 압력진동은 길이방향 1 L 모드(Longitudinal half wave)임을 알 수 있다.

Fig. 7은 1개 단공 원통 그레인의 중앙을 통과하는 중심선을 따라 나타낸 정규화된 압력을 도시한 그래프이다. 약 0-10 mm 위치는 그레인의 단공 내부이며 이후는 그레인 외부이다. 원통 내부에서 축방향 1 L 모드(Longitudinal half wave)를 형성함을 알 수 있다. 이는 앞에서 서술한 주파수 분석 결과와 일치한다.

Fig. 7

Normalized pressure along the center line.

3.2.2. 3개 그레인

Fig. 8은 3개 그레인의 해석결과로 온도장과 그레인 변화를 나타낸 그림이다. 1개 그레인은 처음부터 벽면과 접촉상태로 해석을 진행하였지만 3개 그레인의 경우에는 그레인의 이동모사를 확인하기 위해 Fig. 8(a)와 같이 그레인을 허공에 띄운 상태로 해석을 수행하였다. 시간이 지남에 따라 그레인이 중력에 의해 아래 방향으로 이동하는 것을 볼 수 있다. 가장 아래 그레인의 경우에는 그레인이 연소되면서 지속적으로 아랫방향으로 이동하기 때문에 바닥면에 붙어있는 것처럼 보인다. 1개 그레인과 마찬가지로 연소가 진행됨에 따라 그레인 표면의 연소속도가 가속되고 CBT 내부는 연소가스로 채워지며 온도가 상승한다.

Fig. 8

Temporal evolution of temperature field and grain surface for three-grains burning.

Fig. 9는 3개 그레인의 시간에 따른 CBT 내부압력을 실험결과와 함께 나타낸 그림이다. 1개 그레인과 마찬가지로 실험결과와 수치해석 결과는 같은 압력기울기를 보이며 최종압력은 0.5%내의 차이를 보인다. 반면 확대한 그래프를 보면 CBT의 길이가 길어졌기 때문에 압력진동의 주파수가 6.5 kHz로 1개 그레인의 주파수보다 낮게 나타난다. Eq. 17을 통해 계산된 길이방향 모드(l=1, m=0, n=0)의 주파수는 5.85 kHz이므로 해당 압력진동은 길이방향의 1 L 모드임을 알 수 있다.

Fig. 9

Pressure time history with 2.0-3.0 msec zoom.

Fig. 10은 단공 원통 그레인의 중앙을 통과하는 중심선을 따라 나타낸 정규화된 압력을 도시한 그래프이다. 약 0-30 mm 위치는 3개 그레인의 단공 내부이며 이후는 그레인 외부이다. 1개 그레인과 마찬가지로 각 그레인 단공 내부에서 축방향 1 L모드가 형성된다.

Fig. 10

Normalized pressure along the center line.

4. 결 론

본 연구에서는 고체추진제를 사용하는 CBT 연소 정밀 모델링과 수치해석을 수행하여 그레인의 연소 과정을 관찰하고, 실험결과와 비교 검증하였다. 그레인의 표면연소와 연소면 이동을 모사하기 위해 VOF방법을 사용하였으며 기상과 고상을 함께 해석하기 위해 FSI 기법을 적용하였다. 1개와 3개 그레인 해석을 위해 각각 100개 cpu와 200개 cpu를 사용하여 연소해석을 수행하였다. 연소가 진행됨에 따라 그레인의 연소면 변화와 중력에 의한 이동이 모사됨을 확인하였다. 그레인 내부에서 압력이 상승하고 내부에서 외부로 분출되는 화염 구조가 관찰되었다. 압력-시간 선도를 실험결과와 비교한 결과 수치해석 결과는 실험결과와 동일한 압력기울기를 보였으며 최대 압력은 0.3%의 오차가 발생하였다. 시간에 따른 그레인의 질량변화는 연소 후반부를 제외하고는 해석해와 매우 유사하였다. 1개 그레인에 대한 해석결과로 CBT 중앙에서의 정규화된 압력을 보면 단공 그레인 내부에서는 지속적인 연소가스와 에너지의 생성으로 축방향 정상파(standing wave)가 관찰되었다. 압력섭동의 주파수 분석(FFT 결과)은 12.3 kHz이며 CBT 축방향 이론적 주파수는 11.7 kHz로 1 L 모드가 발생하였다. 3개 그레인의 해석결과, 연소 초반의 압력섭동의 주파수는 6.5 kHz이며 CBT 1 L 모드의 이론적 주파수는 5.85 kHz로 연소 초반에 나타나는 불안정은 1 L 모드이다. CBT 중앙에서의 정규화된 압력을 보면 각 그레인마다 축방향 1 L 모드가 관찰되었으며 단공 내부 압력이 그레인 외부압력에 비해 약 5%높게 형성되었다. 따라서 단공 내부가 외부에 비해 비교적 빠르게 연소된다. 추후 단공형 이외에도 다양한 형상을 고려하여 그레인 설계시 활용할 수 있을 것으로 예측된다.

Nomenclature


c :	Speed of sound
D_c :	Chamber diameter
D_g,i :	Grain inner diameter
D_g,o :	Grain outer diameter
e_grain :	Grain combustion energy
E :	Total specific energy
f :	Frequency
F :	Force
H :	Total specific enthalpy
l :	Longitudinal mode number
L_c :	Chamber length
L_g :	Grain length
m :	Tangential mode number
$m ˙ g r a i n$ :	Grain mass flow rate
n :	Radial mode number
p :	Pressure
$r ˙$ :	Regression rate
R_c :	Chamber radius
t :	Time
T :	Temperature
u_i :	I direction velocity
V_c :	Chamber volume
V_k,i :	Diffusion velocity
Vol :	Volume
Y_k :	K^th species mass fraction
α :	Grain volume fraction
δ_i,j :	Kronecker delta
τ_i,j :	Viscous stress tensor
λ :	Thermal conductivity
ρ :	Density
σ :	Viscous work
$Φ k, i s g s$ :	SGS species flux
$Θ k, i s g s$ :	SGS species diffusive flux

Acknowledgments

본 연구는 방위산업기술지원센터의 지원(사업명: 화포 추진제 및 추진성능 예측 SW기술, 과제번호: UC190003D)하에 수행되었습니다.

References

Lysien, K., Stolarczyk, A. and Jarosz, Tomasz., “Solid Propellant Formulations: A Review of Recent Progress,” Materials, Vol. 14, No. 21, pp. 1-17, 2021. [https://doi.org/10.3390/ma14216657]
Chaturvedi, S. and Dave, P.N., “Solid Propellants: AP/HTPB Composite Propellants,” Arabian Journal of Chemistry, Vol. 12, No. 8, pp. 2061-2068, 2019. [https://doi.org/10.1016/j.arabjc.2014.12.033]
Han, D.H., Yoo, Y.L. and Sung, H.G., “An Analysis of the Different Flow Characteristics of a Closed Bomb Test in Cylindrical and Spherical Closed Vessels,” International Journal of Aeronautical and Space Sciences, Vol. 20, pp. 150-156, 2018. [https://doi.org/10.1007/s42405-018-0116-1]
Han, D.H., Sung, H.G. and Ryu, B.T., “Numerical Simulation for the Combustion of a Zirconium Potassium Perchlorate Explosive inside a Closed Vessel,” Propellants, Explosives, Pyrotechnics, Vol. 42, pp. 1168-1178, 2017. [https://doi.org/10.1002/prep.201700068]
Dunn, S.S. and Coats, D.E., “3-D Grain Design and Ballistic Analysis Using SPP97 Code,” AIAA 1997-3340, 1997. [https://doi.org/10.2514/6.1997-3340]
Wei, R., Bao, F., Liu, Y. and Hui, W., “Combined Acceleration Methods for Solid Rocket Motor Grain Burnback Simulation Based on the Level Set Method,” International Journal of Aerospace Engineering, Vol. 2018, pp. 1-12, 2018. [https://doi.org/10.1155/2018/4827810]
Kim, S. and Kim, S., “A Study on Internal Ballistic Analysis of Solid Rocket Motor Using VOF Method,” Journal of Korean Society of Propulsion Engineers, Vol. 21, No. 3, pp. 61-67, 2017. [https://doi.org/10.6108/KSPE.2017.21.3.061]
Willcox, M.A., Brewster, M.Q. and Tang, K.C., “Solid Propellant Grain Design and Burnback Simulation Using a Minimum Distance Function,” Journal of Propulsion and Power, Vol. 23, No. 2, pp. 465-475. 2007. [https://doi.org/10.2514/1.22937]
Oh, S.H. Lee, S.B., Kim, Y.C., Cha, S.W., Kim, K.R., Kim, D.M., Lee, H.J. and Roh, T.S., “Study on Solid Propellant Grain Burn-back Analysis Applying Face Offsetting Method,” Journal of Korean Society of Propulsion Engineers, Vol. 23, No. 4, pp. 81-91, 2019. [https://doi.org/10.6108/KSPE.2019.23.4.081]
Germano, M., Piomelli, U., Moin, P. and Cabot, W.H., “A Dynamic Subgrid-scale Eddy Viscosity Model,” Physics of Fluids, Vol. 3, pp. 1760-1765, 1991. [https://doi.org/10.1063/1.857955]
Kubota, N., Propellants and Explosives: Thermochemical Aspects of Combustion, 3rd ed., John Wiley & Sons Inc., New York, N.Y., U.S.A., 2015. [https://doi.org/10.1002/9783527693481]
Yoo, Y.L. and Sung, H.G., “Numerical Investigation of an Interaction Between Shock Waves and Bubble in a Compressible Multiphase Flow Using a Diffuse Interface Method,” International Journal of Heat and mass Transfer, Vol. 127, No. 1, pp. 210-221, 2018. [https://doi.org/10.1016/j.ijheatmasstransfer.2018.08.012]
Kitamura, K. and Hashimoto, A., “Reduced Dissipation AUSM-family Fluxes: HR-SLAU2 and HR-AUSM+-up for High Resolution Unsteady Flow Simulations,” Computers and Fluids, Vol. 126, No. 1, pp. 41-57, 2016. [https://doi.org/10.1016/j.compfluid.2015.11.014]
Choi, H.Y., Lee, D.S., Sung, H.G., Lee, W.M. and Kim, E.M., “Burning Rate Estimate Method of Solid Propellants at High Pressure Condition,” Journal of Korean Society of Propulsion Engineers, Vol. 26, No. 1, pp. 28-37, 2021. [https://doi.org/10.6108/KSPE.2022.26.1.028]
Yoo, Y.L., Kim, J.C. and Sung, H.G., “Homogeneous Mixture Model Simulation of Compressible Multi-phase Flows at All Mach Number,” International Journal of Multiphase Flows, Vol. 143, pp. 1-17, 2021. [https://doi.org/10.1016/j.ijmultiphaseflow.2021.103745]
Gordon, S. and McBride, B.J., “Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications: I. Analysis,” NASA, Cleveland, O.H., U.S.A., NASA RP-1311, 1994.
Laudien, E., Pongratz, R. and Preclik, D., “Experimental Procedures Aiding the Design of Acoustic Cavities,” Progress in Astronautics and Aeronautics, Vol. 169, pp. 377-399, 1995.

	unit	experiment	1 grain	3 grains
L_c	mm	110.5	27.63	55.26
D_c	mm	48.0	9.6	11.76
V_c	cc	200.0	2.0	6.0
L_g	mm	9.68	9.68	9.68
D_g,i	mm	1.221	1.221	1.221
D_g,o	mm	4.148	4.148	4.148

	unit	1 grain	3 grains
Grain mass	g	0.2	0.6
Initial pressure	MPa	0.7	0.7
Y_CO		0.4822	0.4822
Y_CO₂		0.1986	0.1986
Y_H₂O		0.1759	0.1759
Y_N₂		0.1334	0.1334
Y_H₂		0.0099	0.0099
e_grain	MJ/kg	2.7120	2.7120

Pressure range	a	n
0-10%	0.5655	0.5131
10-25%	0.6665	0.4750
25-100%	0.1979	0.8440