솔리드 모델링과 레벨셋 기법을 이용한 고체 로켓 추진기관 성능해석

2.1 Burn-Back 해석

그레인 Burn-Back 해석은 고체 로켓 모터가 작동하는 동안 그레인 연소 표면의 진행을 계산하는 것이다. 추진제 연소 과정에서 그레인 연소 표면은 후퇴 진행하여 내부로부터 외부로 수직 전파된다. 이러한 방식으로 고체 상태에 있는 그레인 연소 표면층 아래의 추진제는 고온 가스가 되어 연소면으로부터 방출되어 연소실을 가로질러 유동하며, 최종적으로 로켓 노즐을 통해 빠져나간다.

추진제 그레인 연소 표면의 진행은 고체 로켓 모터의 성능에 결정적인 영향을 미친다. 고체 로켓 추진기관의 성능은 그레인 Burn-Back 해석과 연계된 내탄도 시뮬레이션으로 도출된다. 내탄도 시뮬레이션에는 연소 표면적, 포트 면적, 연소실 부피와 같은 그레인 Burn-Back 해석에서 계산된 입력 자료들을 요구한다.

Burn-Back 해석은 해석적 방법과 수치적 방법이 있다. 해석적 방법은 추진제 그레인의 기하학적 관계로 연소 면적 분석을 수행하는 정확한 해석이다. 그러나 그레인 형상이 복잡해지면 수학적 관계식으로 구성하기 어렵고, 특히 3차원으로 그레인이 구성될 때 거의 해석이 불가능하다. 수치적 방법은 3차원 그레인 형상을 포함한 다양하고 복잡한 그레인에 적용이 가능하나, 수치적 오차가 동반된다. 수치적 방법은 초기 경계치 문제로 알려진 Hamilton-Jacobi 미분 방정식의 해를 구함으로써 경계면인 연소면을 해석할 수 있으며, 이를 레벨셋 기법이라고 한다.

2.2 레벨셋 기법

레벨셋 기법은 Osher와 Sethian의 1988년 논문⁷ 이후 결정 성장, 화염 전파, 그레인 Burn-Back 해석,^8,9 최단 경로 계산, 이미지 분할과 내탄도 같은 많은 분야에 적용되었다.

LSM은 전면 전파와 관련된 문제를 해결하는데 매우 적합한 수치 알고리즘이다. Fig. 1과 같이 이 방법은 전파 속도장이 주어지면 경계 또는 인터페이스의 움직임을 예측할 수 있는데, LSM의 가장 큰 장점은 인터페이스의 토폴로지 변화뿐만 아니라 첨예한 기울기와 복잡한 모양을 쉽게 처리할 수 있다는 것이다.

Fig. 1

Transformation of front motion into initial value problem

LSM은 다음의 Hamilton-Jacobi 미분 방정식의 해를 구하는 방법으로 전개된다.

여기서 Ø는 레벨셋 함수이고, F는 인터페이스 전파속도이다. Ø 함수는 Inviscid Burgers 방정식과 같이 파(Wave)의 전파를 나타내는 함수로 생성된다. Fig. 1에서처럼 레벨 설정 방식에서 전파면 인터페이스는 함수의 0(Zero) 레벨 세트를 사용하여 음적(Implicit)으로 다음 식(2)와 같이 표현한다.

식(2)는 인터페이스 포인트가 함수의 0 레벨셋에 의해 정의되며, 경계를 포함하여 모든 x에 대해 정의된다. 여기서 0 레벨 셋은 계의 경계면인 추진제 그레인면으로 고려할 수 있다.

2.3 레벨셋 기법 적용 방법

추진제 그레인의 형상이 복잡하고 3차원적으로 구성되어 있다면, 압력과 추력에 상당한 영향을 미칠 수 있는 추진제 표면적은 수학적으로나 기하학적인 관계로 계산하기 매우 어렵다.

과거 10여년 전부터 LSM에 대하여 최근에 여러 연구가 진행되고 적용되어 왔다. 추력 요구 조건에 따라 복잡해지는 그레인 형상 해석을 위해 3차원적 면적 분석이 필요하다. 고체 로켓은 대부분 원통형으로 구성되어 있어 원통형 좌표를 사용할 수 있지만, 격자의 종횡비나 비틀림이 없는 이상적인 격자 사용을 위해 데카르트 형식의 계산 격자를 적용하였다. 계산 격자는 Fig. 2처럼 구성이 될 경우 LSM을 적용시킬 수 있다. 이에 따른 3차원 격자구조를 갖는 형상에 따라 식(1)을 이산화하면 다음 식(3)과 같다.

Fig. 2

Values of the discretized Ø function at the grid points surrounding an i,j,k-th point

여기서 F_i,j,k는 경계 전파속도(Interface Propagation Speed)이며, 내탄도 해석에서 사용된 추진제 연속 속도에 해당된다. $$ {\text{∇}}_{i,j,k}^{+}$$와 $$ {\text{∇}}_{i,j,k}^{-}$$는 각각 다음 식(4)와 식(5)와 같다.

그리고, $$ D_{i,j,k}^{-x}$$와 $$ D_{i,j,k}^{+x}$$는 각각 x 방향의 역방향 및 순방향차분 연산자이다.

마찬가지로, $$ D_{i,j,k}^{-y}$$와 $$ D_{i,j,k}^{+y}$$ 그리고 $$ D_{i,j,k}^{-z}$$와 $$ D_{i,j,k}^{+z}$$는 각각 y 방향 및 z의 역방향 및 순방향 차분 연산자이다.

LSM 함수 계산을 위해 추진제 그레인 형상에 따른 격자점을 알아야 한다. 그 격자점을 별도로 형상에 따라 구하는 것도 복잡하고, 어려운 절차로 수행되어야 한다. 하지만 모델링 작업을 용이하게 하는 Solidworks^®나 CATIA^®를 사용하면 어렵지 않게 작업할 수 있으며, 모델링 작업 후 일반 구조해석에서 사용되는 STL (Streolithography) 파일로 저장하면 그 격자점을 사용할 수 있다. STL 파일로 저장하게 되면 삼각형 격자가 자동 생성되고, 면(Face), 모서리(Edge)와 꼭짓점(Vertex)들로 구성되는데, 본 연구에서 이 꼭짓점은 LSM에서 격자점으로 사용하였다. 특히 고체 추진기관은 요구된 연소면을 설정하기 위해 충전치구나 충전용 맨드렐을 이용하여 충전한다.

충전 코어 모델의 표면이 초기 연소면이 되고, 초기 연소면은 법선 방향으로 전파되기 때문에 추진제 연소면을 모델링하고, Matlab^®의 PDE 툴 박스에 포함되어있는 “Createpde”를 사용하여 자동 격자를 구성하면 Fig. 3과 같다. 레벨셋 기법은 추진제 초기 표면으로부터 Burn-Back 되어가는 형상을 지속해서 계산할 수 있는데, 연소관을 벗어나는 계산 표면은 연소면적 계산에서 제외되어야 한다. 따라서 3차원 계산을 위해서 연소관의 라이너와 같은 연소가 일어나지 않는 외곽 형상에 대한 정보를 요구한다. 연소관의 형상은 2차원 AutoCAD 도면에서 작성하면 좌푯값을 생성할 수 있는데, 좌푯값 생성 방법은 명령어 Divide로 직선이나 곡선부를 점으로 분할하고, 분할된 좌푯값을 “Data-Extraction” 명령어를 사용하여 asci 파일로 저장하였다. 연소관은 축 대칭 형상이므로 2차원 형상 위치를 코드 내에서 3차원으로 변환하여 사용하였다.

Fig. 3

Mesh generation of the STL file for the solid propellant grain model

본 논문에서 적용한 연소관 내부의 추진제가 아닌 연소되지 않는 부분 형상은 Fig. 4와 같으며 이를 해석하기 위한 격자점을 갖는 두 번째 입력 파일로 사용하였다.

Fig. 4

Geometry example of combustion case configuration by Solidworks®

2.4 연소면적 계산

LSM 함수 Ø=0인 경우의 경계면이 추진제 그레인 연소면이므로 격자점으로부터 연소 면적을 계산하였다.

Fig. 5는 경계를 나타내는 Ø=0인 경우가 연소면일 경우 삼각형 메쉬의 격자점과 모서리에서 Vector를 보여주고 있다. 3차원 공간상에 있는 Face ①, ②, ③에서 삼각형을 이루는 Edge Vector가 $$ \overrightarrow{a},\overrightarrow{b}$$라 하면 Face의 면적은 다음 식(12)와 같이 계산된다.

Fig. 5

Grid point and mesh edge vector for the calculation of the burning area

마찬가지로 각 Face의 면적을 고려하면 전체 연소 면적을 식(13)과 같이 쓸 수 있다.

2.5 성능해석

고체 로켓 추진기관의 내탄도 해석은 각 시간단계에서 계산된 식(13)의 연소 면적을 계산한 후 식(14)부터 식(17)에 따라 압력을 Runge-Kutta 4th Order를 사용하여 연소실 압력을 구하고, 이에 상응한 추력을 계산하였다.

여기서 r_b는 추진제의 연소속도로 $$ r_b=a{p}_c^n$$ 이다.

연소실 가스가 노즐로 배출될 때, 배출 질량 유량은 노즐 목에서 쵸킹될 경우와 쵸킹이 일어나지 않는 경우에 대하여 다음 식(18)과 식(19)와 같이 사용된다.

(1) Non Chocking

(2) Chocking

본 연구에서 성능해석 코드는 Matlab^® 2020a를 이용하여 Fig. 6과 같은 절차로 개발하였다. 개발된 코드는 Solidworks^®를 이용하여 모델링 작업 후 형성된 메쉬에서 각 노드를 구하여 LSM을 해석하기 위한 격자점으로 사용하였다. LSM에 따라 식(13)의 연소 면적을 계산한 후, 0D 준정상상태(Quasi-Steady State) 조건의 내탄도 해석을 수행하였다.

Fig. 6

Flow chart for the performance analysis by using Matlab® code

3.1 NAWC 추진기관 성능해석 적용

본 연구에서 성능해석의 검증을 위해 NAWC Motor No. 13을 사용하여 해석하였다. 추진기관 형상은 Fig. 7과 같이 단순한 원통형의 형상이며, 노즐 목 직경은 26.42 mm이고, 추진기관 그레인 외경은 121.92 mm이며, 길이는 860 mm인 로켓 모터이다.

Fig. 7

NAWC motor No. 13 propellant grain geometry

적용 추진제의 연소속도는 6.9MPa (@1,000 psia)에서 5.41 mm/s이고, 압력 지수는 n = 0.46이다.

계산을 위해 3D 모델링 작업은 Solidworks^®로 작업을 수행하였는데, Fig. 8과 같이 LSM을 해석과 내탄도해석을 위해서 STL 파일로 저장하였다.

Fig. 8

3D Solid model for the propellant grain

Fig. 9는 추진제 그레인의 시간에 따른 형상 변화를 보여준다. 그레인의 전방과 후방은 인히비트(Inhibit)되어 원통 부분만 연소되어 진행됨을 보여준다. 외곽 음영 처리된 부분은 연소관 부분이고, 내부 형상은 시간에 따라 연소면의 변화를 나타내는데, 시간에 따라 반경 방향으로 전파되어 내탄도에서 사용되는 연소 면적이 증가됨을 보여주고 있다.

Fig. 9

Evaluation of the grain burning surface of NAWC No. 13

성능해석 결과는 Fig. 10과 같이 압력-시간, 추력-시간 성능 선도로 나타내었다. 성능 선도는 증가하는 연소 면적에 따라 압력의 증가 양상을 나타내고있다.

Fig. 10

NAWC motor No. 13 ground test and numerical 0D quasi steady simulation with other simulated results

지상 연소 시험과 비교하면 연소 시험 결과와는 제작상의 오차, 추진제 특성치의 편차 및 추진제 초기 표면 온도 등의 원인으로 차이가 발생되지만, 연소 초기에 압력의 튐이 있는 점을 제외하고, Wilcox, et al.²과 Cavallini, et al.³ 결과들과 0D Quasi-Steady 해석 결과와 유사함을 보여주고 있어서, 개발된 코드가 성능해석에 있어서 적합하다고 판단한다.

3.2 6인치 시험용 추진기관 성능해석 적용

성능해석 검증을 위하여 연소 시험 결과가 존재하는 6인치 시험용 추진기관에 적용하였다. 추진기관은 Fig. 11과 같이 추진제 성능 시험과 노즐 내삽 효과를 위하여 제작된 시험용 추진기관이다.¹⁰

Fig. 11

Schematic of an experimental rocket motor (Upper: Test motor, Lower left: Charged propellant of the motor, Lower right: Feature of ground test)

추진기관은 추진제 길이가 264이고, 외경은 67 mm인 중공형 실린더 형상이다. 노즐 목 직경은 25이고, 노즐 외경은 70 mm이다. 적용된 추진제 특성으로 연소속도는 1,000 psia에서 11.1 mm/s이며, 압력 지수가 n = 0.43인 HTPB 추진제이다. 성능해석을 위해 추진제 그레인 연소면은 Fig. 12와 같이 실린더로 구현되고, 양 끝단에서 축 방향으로 연소가 되기 때문에 디스크 모양의 연장된 부분을 추가시켰다.

Fig. 12

3D Solid model for the test motor propellant grain

Fig. 13은 6인치 추진기관의 성능해석 결과를 보여준다. 추진제 연소면은 원통 내면과 양 끝단에 존재하여 시간이 지날수록 원통부는 넓어지고, 양 끝단 연소면은 중앙 방향으로 서로 가까워지면서 연소된다. 2.3 sec에서는 연소가 종료되어 연소관 형상만 존재하고 있음을 보여준다.

Fig. 13

Grain burn-back evaluations of 6 inch test motor

추진기관의 성능을 나타내는 압력-시간, 추력-시간은 Fig. 14로 실제 지상 연소 시험 결과와 비교하였다. 추진제 연소실 형상이 비교적 단순하고, 길이에 대한 직경비가 L/D-2.2로서, 예측된 해석값이 시험값에 잘 일치되고 있음을 보여준다.

Fig. 14

Comparison of the analysis (Predict) and experimental results (Upper: Pressure time history, Lower: thrust time history)

3.3 2.75인치 추진기관 성능해석 적용

세 번째로는 Fig. 15처럼 기존 2.75인치 추진기관을 대상으로 당사에서 개발된 추진기관에 적용하였다. 사용된 추진제는 바인더로 HTPB의 질량 분율이 8.5%, 산화제인 AP는 87%, 그리고 가소제와 경화제등 첨가제 4.5%로 구성된 HTPB 추진제로 연소 속도는 1,000 psia에서 9.58 mm/s이고, 압력 지수 n = 0.39이다.

Fig. 15

LIG 2.75 inch rocket motor

마찬가지로 해석에 앞서 솔리드 모델링을 통한 초기 연소면은 Fig. 16과 같다. 형상은 꼭짓점이 6개인 Star 그레인이고, 노즐 쪽은 연소관과 추진제 형상에 따라 덮개 모양으로 모델링하였다.

Fig. 16

Initial grain surface for LIG 2.75 inch rocket motor

Fig. 17은 시간 간격 0.25 sec마다 추진제 그레인의 Burn-Back 모습을 나타내었다. 1.75 sec부근에서 별(Star) 꼭짓점이 연소관에 도달하면 연소면이 사라지면서, 2 sec 이후 모든 추진제 그레인이 소멸하여 연소가 종료되고 있음을 가시적으로 보여주고 있다.

Fig. 17

Grain burning surface evolution of LIG 2.75 inch rocket motor

Fig. 18은 2.75인치 로켓 모터의 시험대 장착 모습과 시험 장면을 나타내었다. 지상 연소 시험에서 계측한 성능 데이터와 본 연구에서 개발된 코드로부터 결과를 비교한 성능곡선은 Fig. 19에 나타내었다. 연소 초기에 압력과 추력의 예측값과 시험값은 많은 차이가 발생하였다.

Fig. 18

Ground test for the performance of 2.75 inch rocket motor

Fig. 19

Histories of pressure and thrust of 2.75 inch rocket motor

세장비(L/D)가 9.16인 NAWC No. 13 추진기관 해석에서도 연소 초기에 약간 나타난 동일한 현상으로, 당사에서 개발한 2.75인치 추진기관은 매우 큰 L/D-12.3인 추진기관으로, 침식연소(Erosive Burning) 효과가 연소 초기에 상당한 비중이 있다고 분석된다. 연소 초기에 추진제 그레인 포트의 직경이 작아서 연소면을 따른 유속이 연소 표면으로 전달되는 열 유속을 증가시켜 연소속도가 증가되기 때문이다. 침식연소에 따른 추진제 질량 유속 증가는 연소실 압력의 상승을 유발시켜, 예측된 설계값과의 오차가 발생하고 있다고 분석된다. 추진제 연소가 진행됨에 따라 포트의 직경이 증가하면서 유동에 따른 질량 유속이 작아져서 침식연소 효과가 소멸되고, 식(17)의 정상연소속도의 질량 발생만 존재하여 예측값으로 접근한다. 추력-시간 성능 선도이하 면적을 나타내는 총 역적(Total Impulse)이 동일함을 고려하면 초기 상승된 압력에 따라 연소 시간은 짧아진다. 본 연구에서는 대부분 추진기관(L/D < 8)에서 사용하는 완전 해석적 방법(0차 해석)을 적용하였는데, 추진기관 내부 유동에 의한 연소속도 증가를 모사하기에는 한계가 있다고 생각된다. 소형 추진기관에 따른 설계에 대하여 추진제 충전에서 발생하는 제작 오차도 실험값과의 편차에 기여함도 무시할 수 없다. 그러나 제작상의 오차를 제외하고, 정확한 해석을 위해서 추진기관 내부 연소실의 유동을 동반한 1차원 해석이 수행된다면 보다 정확한 결과를 도출할 수 있다고 판단되고, 추가적인 연구와 코드의 보완을 진행하고 있다.

3.4 Finocyl 추진기관 성능해석

다음의 해석 대상은 향후 개발 예정인 Fig. 20과 같은 직경 160 mm급 Finocyl (= Fin+Cylinder) 추진제 그레인 형상을 갖는 추진기관이다.

Fig. 20

160 mm rocket motor configuration with a propellant grain core

추진기관의 길이는 1,550 mm이고, Finocyl 추진제 그레인은 길이가 1,240 mm이며, 원통과 4개의 슬롯으로 구성된 대표적인 3차원 형상이다. 이는 연소 면적의 Burn-Back 해석에서 수학적으로 해석하기 어렵다.

본 연구에서 LSM과 연계된 수치적 Burn-Back 해석을 위해 추진기관 연소실 형상의 맨드렐을 Fig. 21처럼 Solidworks^®로 모델링하여 삼각형 격자를 갖는 STL 파일로 구현하였다. 추진제 그레인은 슬롯과 원통부를 포함하였고, 슬롯 하단 단면도 전방으로 진행하는 연소면이 존재하기 때문에 슬롯 후방부에 메쉬를 추가하여 모델링하였다. 연소가 진행되면 바로 후방의 원판부는 사라지고, 슬롯 단면만 남는 그레인이 Burn-Back할 수 있도록 형상화하였다. 또한 충전 코어의 이형을 위해서 전방원통부의 직경과 슬롯 전단 직경을 각각 50과 54 mm의 형상으로 경사지게 취하였다. Fig. 21은 LSM을 계산하기 위해 Matlab^® 코드에서 사용하는 pdeplot3D 모듈로 그려진 그레인의 메쉬 형상이다.

Fig. 21

Triangle mesh for a propellant grain core mandrel

격자들은 면과 모서리 그리고 꼭짓점으로 구성된 삼각형 메쉬로 자동 생성된 형상이다. 이 형상과 연소관의 형상 격자점 파일을 입력으로 하여 LSM과 내탄도 해석을 수행하였다. 적용된 추진제 특성과 노즐 형상 데이터는 Fig. 22에 나타내었다.

Fig. 22

Input data of propellant characteristics and nozzle data

동일한 절차에 따라 해석된 결과로 추진제 그레인 연소면의 형상 변화는 Fig. 23에 나타내었다. 추진제는 시간에 따라 그레인이 Burn-Back되는 형상을 보여주고 있으며, 슬롯 부분은 3 sec를 지나면서 연소관 벽면에서 소멸되어 진행되고, 5 sec 경과 이후 원통형만 남아 소멸(Tail-Off)된다. 이는 Fig. 24에서 산출된 성능곡선으로 확인된다. 추진제의 온도민감계수가 0.18인 경우를 적용하여, 초기 표면 온도가 각각 상온(20^oC), 저온(-32^oC), 그리고 고온(43^oC)에서 성능곡선을 도출하였다. 상온을 기준으로 저온에서 낮은 압력/추력과 긴 연소 시간을, 고온에서 높은 압력/추력과 짧은 연소 시간이 얻어지고 있음을 보여준다. 표준 상태인 상온에서 해석된 대표적인 성능 수치는 다음과 같이 얻었다.

Fig. 23

Change of the grain burning surface with each time for the finocyl grain rocket motor

Fig. 24

Performance curves for 160 mm rocket motor