발간 현황

연소표면적을 계산하기 위한 방법에는 크게 해석적 방법과 수치적 방법이 있다. 해석적 방법은 그레인 형상을 사면체, 육면체, 구 등의 기하학적 형상으로 나누어 계산하는 것이다. 이러한 방법은 정확한 해를 도출하지만 복잡한 3D 그레인 형상에 적용하는 것에는 한계를 가진다.

수치적 방법은 그레인을 유한요소로 나누어 계산하는 것으로 자유 경계면을 정의하는 방법에 따라 라그랑지안 법(Lagrangian method)과 오일러 법(Euler method)으로 나뉜다[4]. 라그랑지안 법은 유체입자 자체의 움직임을 따라 경계면을 계산하는 방법으로써 대표적으로 마커-스트링(Marker-string)법이 포함된다. 라그랑지안 법은 자유 경계면에 있는 유체 입자를 따라 계산하기 때문에 경계면 이동 중 입자가 사라지거나(Dissipation) 여러 입자가 겹치는(Shock) 등의 계산상 단점을 가진다[4].

이와 달리 오일러 법은 유체입자의 움직임을 따르지 않고 계산 영역 내 셀(cell)의 유체 특성 값을 계산하는 방식으로, 라그랑지안 법이 가지는 단점을 보완할 수 있다.

VOF와 레벨셋(level set) 등이 대표적인 오일러 법으로, 이 중 VOF는 셀 내에 유체가 차지하는 부피를 이용하여 자유 경계면을 정의하는 방법이며, Fig. 1은 VOF 기법의 개념을 보여준다[5-7]. 일반적으로 VOF는 레벨셋 방법보다 곡률체적 및 상 질량 보존에 장점을 가지는 것으로 알려져 있다[10]. 따라서 본 연구에서는 VOF를 이용하여 연소면을 추적하고, 연소속도에 따라 후퇴하는 연소표면적을 계산하였다.

Fig. 1 
A Concept of VOF method[7].

연소 현상을 고려하기 위하여 반응진행변수(C)를 고려한 C-방정식 모델을 사용하였다[8]. 연소 물질 내에서 연소 유무와 관련된 반응진행변수의 전달해석은 Eq. (1)과 같은 연소 스칼라 전달방정식으로 표현된다.

여기서, c¯는 평균반응진행 변수, S는 Schmidt 수, S_C는 반응진행 소스항을 나타낸다. Schmidt 수는 동적점성도(ν)에 대한 질량확산도(D)의 비로 나타내며 Eq. (2)와 같고, 반응진행 소스항은 Eq. (3)과 같다.

여기서 u는 화염속도를 의미한다. 추진제를 단일 물질로 가정하고, 고체 추진제의 경우 대류에 의한 확산을 고려하지 않으면 Eq. (1)은 Eq. (4)와 같이 정리된다.

c는 0과 1사이의 값을 가지며, 완전 연소 시 c=1, 비 연소 시에는 c=0의 값을 갖는다. 연소 경계면은 c=0.5일 때로 정의하였다.

본 연구에서 개발한 BAAP의 신뢰성을 검증하기 위해 NAWC no.6 모터를 대상으로 계산을 수행하였으며 연소표면적 계산은 다음과 같은 과정으로 수행하였다.

2.3.1 3D 모델 생성

본 연구에서는 3D 모델 생성을 위해 상용 소프트웨어인 Solidworks 2016을 사용하였다. 생성된 NAWC no.6 형상은 Fig. 2와 같다[2,4]. 전체길이는 약 1,690 mm, 외경은 63.5 mm이며 축 방향을 따라 그레인의 단면이 연속적으로 변화하는 형상을 갖는다[3,5]. 대상 모델은 축 대칭 형상을 가지므로 계산시간의 단축을 위해 1/6 모델을 사용하였다.

Fig. 2 
A schematic picture of NAWC no.6 motor grain.

2.3.3 격자 생성(Meshing)

BAAP에서는 사용자가 최대 격자크기(∆x_max) 를 설정할 수 있으며, 설정된 격자크기에 따라 데카르트 격자(cartesian mesh)를 생성한다. NAWC no. 6 모터 1/6 모델의 경우 최대 격자크기가 2 mm일 때 약 32만개, 최대 격자크기가 1 mm일 때는 약 280만개의 셀이 생성되었다.

Fig. 3은 데카르트 격자를 사용하여 각각 최대 크기 2 mm와 1 mm인 격자를 생성한 결과이다. 격자 크기가 충분히 작지 않으면 그레인 표면 형상이 격자를 따라 고르지 못한 것을 확인할 수 있었으며, 격자 크기에 따른 해석 결과의 차이는 2.4절에서 비교하였다.

Fig. 3 
Captured pictures of grain mesh.

2.3.4 연소속도 및 단위 계산시간 설정

2.2절의 Eq. (3)으로부터 난류화염속도(u_t)가 반응진행 소스항을 결정하는 것을 알 수 있다. 난류화염속도가 클수록 반응진행 소스항이 커지고, 연소면이 빠르게 확산되어간다. BAAP에서는 사용자가 난류화염속도를 설정할 수 있으며, 본 연구에서는 1~30 mm/s 범위 내에서 계산을 수행하였다. 또한 수치해석계산을 위해 단위 계산시간(∆_t)은 0.01~0.05 s 범위에서 설정하였다.

그레인 웹 두께(d_web)와 난류화염속도를 통해 결정된 계산 종료시간(t_end)을 이용하여 BAAP에서는 t=0초부터 t=t_end초까지 C-방정식 모델을 계산하고 계산 시간별로 결과를 저장한다. Fig. 4는 BAAP 결과를 나타낸 예로 연소면 후퇴현상을 확인할 수 있다.

Fig. 4 
An example of graphical burnback result.

BAAP에서 사용자가 변경가능 한 매개변수는 격자 크기, 난류화염속도, 단위 계산시간이다. 이 변수들은 CFL (Courant- Friedrichs- Lewy condition)수를 정의하는 매개변수이며, 그 관계는 Eq. (5)와 같다.

일반적으로 수치해의 수렴을 위해 CFL 수는 1보다 작아야한다[1]. 본 연구에서는 매개변수에 따른 연소표면적 결과의 의존성을 평가하기 위해 Table 1과 같은 조건에서 계산을 수행하였다.

계산 결과는 3D 모델링 소프트웨어(Solidworks 2016)를 이용한 연소표면적 산출 결과(Table 1의 Ref.)와 비교하여 신뢰성을 검증하였다. 3D 모델링을 이용한 면적 계산법은 초기 추진제 그레인의 형상을 모델링하고 초기 연소면이 표면에 수직방향으로 일정 길이만큼 후퇴하는 것을 반복적으로 모델링 한 뒤, 각 모델에 대한 연소 표면적 정보를 3D 모델링 소프트웨어를 통해 얻는 방법이다. 이러한 방법은 많은 노력과 시간을 요하는 단점을 가지나 기하학적 형상 변화를 실제적으로 반영하여 매우 단순하면서도 정확한 값을 제공하는 장점이 있어 본 계산 결과와 비교하였다.

2.4.1 동일 격자크기 내에서 CFL 수의 영향

Table 1의 Test 1~3은 동일한 격자와 CFL 수를 가지며, 그 결과는 Fig. 5와 같다. 동일 격자크기 내에서 CFL 수가 같으면 연소표면적 계산 결과는 동일하였다. 격자크기가 동일할 때 CFL 수는 난류화염속도와 단위 계산시간에 의해 결정된다. 따라서 Fig. 5와 같은 결과는 난류화염속도와 단위 계산시간에 관계없이 격자크기와 CFL수가 같으면 연소표면적 결과도 동일함을 의미한다. 격자크기는 연소면 모사 정확도에 영향을 주고, CFL수에 의해 수치적 수렴성이 결정되기 때문에 CFL수와 격자크기가 동일하면 연소표면적 계산결과도 동일하게 나타나는 것으로 사료된다.

Fig. 5 
The results of test 1, 2 and 3.

이러한 결과를 바탕으로 동일 격자크기 내에서 CFL 수의 영향을 확인하기 위해 Test 3~7의 조건을 비교하였으며 그 결과를 Fig. 6에 나타내었다. CFL 수가 커질수록 웹 두께가 길어지고, 전체적인 경향이 Ref.와 멀어지는 것을 확인하였다. 이는 CFL 수 조건이 수치해석 수렴의 범위를 넘어 과도하게 주어졌기 때문으로 판단된다.

Fig. 6 
The effect of CFL number with identical maximum mesh size.

Test 3, 4의 경우 웹 두께 및 연소표면적 증감의 경향은 Ref.와 가장 유사하나 웹 두께 5~30 mm에서 연소표면적이 과대추정 되었다. 이는 연소면이 원주방향에 접근할수록 격자의 왜곡이 증대되어 발생한 결과로 사료된다.

2.4.2 격자 크기의 영향

Fig. 7은 CFL 수가 동일할 때 최대 격자 크기에 따른 영향을 비교한 결과이다. CFL 수가 0.05, 0.1인 경우 격자 크기가 커지면 연소표면적의 과대추정 정도는 다소 감소하는 경향을 보이나 웹 두께는 증가하는 결과가 나타났다.

Fig. 7 
The effect of maximum mesh size.

해의 정확도를 확보하기 위해 격자의 크기는 추진제 형상의 가장 작은 곡률을 기준으로 설정되어야 한다. NAWC no.6 모터의 경우 약 4.8 mm이므로 본 연구에서는 격자 크기 1 mm를 최소 기준으로 분석하였다.

내탄도해석 시 추진제 질량은 총역적을 결정하는 중요한 변수이다. Test No. 1~10의 연소면적 결과를 바탕으로 부피를 역 추산한 결과는 Table 2와 같다. 3D 모델을 바탕으로 계산한 결과 대비 최소 14.7%의 차이를 보이며, 이는 전술한 바와 같이 격자의 크기와 형태, 경계면 정의에 의해 나타나는 것으로 판단된다.

Table 2. 
Volume and error based on the ref.

Test No.	Volume Based on the Burning Area (m3)	Error Based on the Ref. (%)
Ref.	0.0181	-
1	0.0208	14.8
2	0.0207	14.7
3	0.0208	14.8
4	0.0209	15.4
5	0.0224	12.4
6	0.0267	14.8
7	0.0348	19.3
8	0.0212	17.1
9	0.0213	17.8
10	0.0212	12.2

이러한 오차는 본 프로그램을 실제 내탄도 해석에 적용하기에는 다소 높은 것으로 추후 개선될 필요성이 있다.

추진제는 참고문헌에서 제시된 저연 추진제 조성을 사용하였다[5]. 추진제 주 조성 및 연소 속도 관련 변수는 Table 3과 같다. 변수 a, n, rb는 각각 연소속도 압력상수, 연소속도 압력지수, 연소속도를 나타내며 MPa와 mm/s를 기준으로 한다. 연소속도는 1000 psia에서 6.75 mm/s이다.

Table 4의 추진제 조성을 기초로 NASA CEA 코드를 사용하여 화학평형계산을 수행하였으며[9], 결과는 Table 4에 나타내었다. 본 해석에 적용한 추진제는 연소온도 2974.3 K, 특성속도 1527.5 m/s, 비추력 247.3 s의 특성을 가진다.

본 내탄도해석에서는 Table 1, Test 4의 연소표면적 결과를 기초로 하였다. 단열상태를 가정하여 0-D에 대한 계산을 수행하였고, 추진제 연소속도는 r_b=αPⁿ에 의하여 압력에 따라 변화한다. 시간에 따른 연소표면적 및 연소관 내부 압력은 Fig. 8과 같이 나타났다. 본 연구에서는 웹 연소시간 4.98초, 웹 평균압력 5.3 MPa, 최대 압력 7.8 MPa로 나타났으며, 기존의 연구결과인 웹 연소시간 4.70초, 웹 평균압력 4.4 MPa, 최대압력 6.0 MPa 비교하여 평균압력 기준 약 20%의 차이를 보였다.

Fig. 8 
Burning area and pressure with time.