1. 서 론

고체 추진기관은 연소 중 임의로 추력을 조절하기 어렵다. 이 때문에 다양한 임무에서 요구되는 추력 조건을 만족하기 위해 사전에 추진제를 적절하게 성형하여야 한다. 이를 위해서는 추진제의 초기 형상에 따른 추력 거동 예측이 필수적이다. 추진기관의 추력은 연소실 내로 유입되는 유량에 의존한다. 고체 추진기관의 경우, 추진제의 표면에서 연소가 진행되므로 연소면적이 유입 유량과 추력 거동에 큰 영향을 미친다[1]. 이 때문에 고체 추진제의 초기 형상에 따라 나타나는 연소면적 변화 예측에 관한 다양한 연구들이 진행됐다.

연소면적의 변화 문제는 임의의 경계면의 각 요소가 방향과 속도를 갖고 움직이는 이동경계면 문제라고 볼 수 있다[2]. 이러한 문제에 대한 기존 기법에는 해석적 기법, Drafting 기법이 있다. 해석적 기법은 기하학적 수식을 이용하여 그레인의 형상 변화를 나타내는 방법으로, 가장 정확도가 높다[3,4]. 그러나 그레인의 형상이 복잡하거나, 연소가 진행되며 기하적 위상 변화(topology change)가 나타나는 경우 이를 모사하기 위한 여러 복잡한 수식이 필요하다는 단점이 있다. 또한, 같은 이유로 3D 문제에 적용이 어렵다. Drafting 기법은 형상 설계 소프트웨어를 이용하여 연소면적 변화를 예측하는 방법이다[5,6]. 이 기법 역시 정확도가 높지만, 설계자가 직접 소프트웨어를 조작해야 하므로 시간이 오래 걸리고 절차가 복잡하여 인적 요소에 의한 오차 위험이 있다.

기존 기법을 보완하기 위해 많은 선행연구에서 연소면적 해석(burn-back analysis)을 더 간단하고, 효율적으로 수행하기 위한 수치적 접근 방식을 연구하였다. 이동경계면 문제에 대한 수치적 접근 방식은 라그랑지(Lagrangian) 접근 방식과 오일러(Eulerian) 접근 방식이 있다. 라그랑지 접근 방식의 경우 경계면의 형상을 여러 개의 점으로 직접 표현하며, 경계면의 이동을 점들의 움직임으로 기술한다. 이에 대하여 Marker Particle Method[7], Partial Interface Tracking Method[8], Surface Tracking Method[9]와 같은 다양한 기법들이 연구되었다. 이러한 기법들은 정확도가 높고 계산시간이 짧다는 장점이 있다. 그러나, 경계면의 확장, 수축, 충돌 등으로 인한 기하적 위상 변화가 나타나는 경우, 복잡한 예외 처리가 필요하다는 단점이 있다.

반면에 오일러 접근 방식은 경계면의 형상을 주변의 격자점과 경계면 사이의 상대 정보를 통해 간접적으로 표현하는 방식이다. 대표적인 기법으로는 Level Set Method (LSM)[10,11], Volume of Fluid Method(VOF)[12]가 있다. 이러한 기법들은 지배방정식을 푸는 과정에서 기하적 위상 변화가 자연스럽게 반영되기 때문에 많은 예외 처리가 필요하지 않고, 복잡한 형상에 대해서 강건하다는 장점이 있다. 게다가 경계면의 이동 속도가 경계면 위치, 방향에 의존하여도 적용할 수 있어 이질성, 이방성을 갖는 실제 추진제에 대한 적용도 가능하다[2,13]. 그러나 경계면의 표현 정확도를 높이기 위해서 많은 격자점이 필요하기에 계산량이 크다는 단점이 있다.

한편, Fast Marching Method(FMM)는 오일러 접근 방식의 기법으로, 보편적으로 알려진 LSM의 계산량 문제를 개선하기 위해 Sethian 등이 제안한 기법이다. 이들의 연구에서는 FMM을 적용할 경우, 계산 복잡도가 O(N2logN) 임으로 O(N2tfinal/△t)의 계산 복잡도를 갖는 LSM이나 VOF 기법보다 계산시간을 개선할 수가 있음을 주장하였다[14,15]. FMM의 경우에는 기존 기법과 달리 이질성, 이방성을 갖는 고체 추진제의 연소면적 해석에 적용된 사례를 찾기는 어렵다. 하지만 균일, 등방성 가정이 가능한 고체 추진제의 연소면적 해석 문제에 FMM을 적용한다면, 여러 복잡한 형상에 대해서도 강건한 결과를 얻음과 동시에 계산시간에 대한 개선도 가능해 보인다. 그러나 FMM을 연소면적 해석 문제에 적용하고 그 결과를 기존 기법의 결과와 정량적인 기준을 통해 비교 분석한 문헌은 찾기 어렵다. 따라서, FMM이 연소면적 해석 문제에서 적용되었을 때 실질적으로 계산시간의 개선이 가능한지는 확인이 필요하다.

본 연구의 목적은 균일, 등방성 가정이 가능한 추진제에 대하여 FMM을 적용한 연소면적 해석 프로그램을 개발하고 계산시간의 개선 정도를 보편적으로 사용되는 LSM 기반 프로그램과 비교하는 것이다. 이를 위해 Matlab 소프트웨어를 이용한 환경에서 해석 프로그램을 개발하였다. 또한, 두 기법이 동등한 계산 정확도를 갖는 조건을 확인하고, 여러 그레인 형상에 대한 해석을 수행하여 FMM과 LSM 기반 프로그램의 계산시간을 비교하였다.

2. 이동경계면 해석

연소면적 해석을 위해서는 이동경계면의 해석과 후처리 과정이 필요하다. 이 장에서는 이 논문에서 사용된 기법인 LSM, FMM에 대해서 다룬다. 또한, 두 기법으로부터 얻은 격자에서 연소관의 형상을 반영하고, 연소면적을 계산하기 위한 후처리 과정을 요약하였다.

한편, 균일, 등방성 추진제라면, 이동경계면의 모든 점에서 경계면 이동 속도가 일정하며, 표면에서 수직 방향으로 연소가 진행된다는 가정이 유효하다[1,2,3,4,5]. 이와 같은 가정을 통해 이동경계면 해석을 진행하면 일정 시간 간격을 갖는 연소면적 데이터를 얻을 수 있다. 이때, 경계면의 이동 속도가 일정하므로 이를 경계면 이동 거리 혹은 연소거리(burn distance)에 따른 테이블로 나타낼 수 있다. 내탄도 해석을 수행할 때는 매 순간의 연소 압력을 통해 얻은 연소속도를 적분하여 연소거리를 구하고, 연소면적 테이블을 참조하여 순간 연소면적을 구할 수 있다. 본 논문에서도 연소면적 해석 결과를 연소거리에 대해 나타내었다.

2.1 Level Set Method

LSM에서는 Signed Distance Function(SDF)을 이용하여 경계면을 표현한다. SDF는 각각의 격자점에서 경계면까지의 최소거리와 경계면 이동 방향에 대한 상대위치에 의해 결정되는 부호로 나타내어진다. Fig. 1에 1차원 SDF의 예시가 나타내져 있다. 그림에서 SDF 값이 0인 부분이 경계면의 현재 위치에 해당하며, 경계면 이동 방향의 반대 방향에 있는 격자점들의 값은 음의 부호로, 경계면 이동 방향의 격자점들의 값은 양의 부호를 갖는다. 2차원, 3차원의 경우, 격자점들의 값이 0인 부분을 모두 찾아낸다면 경계면의 모양을 알 수 있다.

Fig. 1.

1-D Signed distance function example.

추진제의 모든 면에서 순간 후퇴율이 일정하며, 표면의 수직 방향으로 연소가 진행된다고 가정하면, LSM에서는 Eq. 1을 이용하여 경계면의 이동을 표현할 수 있다[10].

Eq. 1에서 𝜙는 SDF 값, t는 시간, V는 경계면의 이동 속도를 나타낸다. 모든 면에서 순간 후퇴율이 일정하므로 V는 상수다. 본 논문에서는 Eq. 1을 다음과 같이 1차 정확도로 이산화하고, 수치 적분하여 계산하였다.

Eq. 2에서 n은 n번째 시간 간격, i,j는 각각 x 및 y 방향 격자점 인덱스, △t는 시간 간격을 나타낸다. 차분 연산자 D는 각각 다음과 같다.

Eq. 3, 4, 5, 6에서 △x와 △y는 각각 x 방향, y 방향 격자 크기를 나타낸다. 이때, 해의 안정성을 보장하기 위해서 CFL값이 특정 값을 만족할 수 있도록 속도, 시간 간격, 격자 크기를 조정하여야 한다. CFL값 및 변수들의 설정에 관련된 내용은 본 논문의 3.1절에 설명되어있다. CFL값은 Eq. 7과 같이 정의된다.

이처럼 구현된 LSM의 경우, 이동경계면 문제를 N×N의 2차원 직교 격자를 이용하여 해석할 때 N2tfinal /Δt의 계산 복잡도를 갖게 된다. 이때, tfinal은 최종 시간을 의미한다.

2.2 Fast Marching Method

LSM이 경계면의 위치 표현을 위해 SDF를 이용하는 것과 다르게 FMM은 경계면이 각각의 격자점에 도달하는 도착시간(arrival time)을 이용하여 경계면의 위치를 표현한다. 또한, LSM에서는 n번째 시간 간격에서의 경계면 위치를 알기 위해서 Eq. 2를 n번 계산하여 n번째 격자 정보를 얻어야 하지만, FMM에서는 도착시간 정보가 한 격자 안에 모두 포함되어 있다는 차이점이 있다. Fig. 2에 FMM을 통해 얻은 1차원 격자의 예시가 나타내져 있다. 그림에서 경계면의 초기위치가 0으로 표현되어 있으며, 경계면의 이동 방향 쪽으로 각각의 격자점이 경계면의 도착시간 값을 갖는다. 2차원 이상의 해석을 수행한 경우, 격자점들의 값이 같은 모든 점을 찾아낸다면 임의의 원하는 시간 t에서의 경계면의 위치와 모양을 알 수 있다.

Fig. 2.

1-D Arrival time example.

Eq. 1과 마찬가지로, 추진제의 모든 면에서 순간 후퇴율이 일정하며, 표면의 수직 방향으로 연소가 진행된다고 가정한다면, FMM에서는 다음의 식을 이용하여 경계면의 이동을 표현할 수 있다[14,15].

Eq. 8에서 T는 도착시간을 나타내고, V는 경계면의 이동 속도를 나타낸다. Eq. 1과 마찬가지로 V는 상수이다. 본 논문에서는 Eq. 8을 다음과 같이 1차 정확도로 이산화했다.

Eq. 9에서 차분 연산자 D는 Eq. 3, 4, 5, 6과 같은 정의를 갖는다. Eq. 9에서 구하고자 하는 격자점에서의 도착시간 Ti,j가 우변에 포함되어 있으므로, Eq. 10과 같이 정리하고 Eq. 11과 같이 근의 공식을 활용하여 값을 얻을 수 있다.

Eq. 10, 11에서 a,b,c는 모두 상수다.

LSM의 경우 n+1번째 시간 간격의 격자 정보를 얻기 위해 n번째 격자의 정보가 사용되므로 각각의 격자점을 순서에 상관없이 계산하면 다음 시간 간격의 정보를 얻을 수 있다. 그러나, FMM의 경우에는 한 격자 안에 모든 도착시간 정보를 포함하므로 아래의 알고리즘에 따라 계산을 수행한다.

Fig. 3에 FMM 알고리즘의 설명이 그려져 있다. 계산이 진행되는 동안 격자점들은 계산이 끝난 상태(Accepted), 계산 중인 상태(In-band), 모르는 상태(Unknown)의 세 가지 상태를 갖게 된다. FFM은 이동경계면의 초기형상을 통해 구한 Accepted 상태의 격자점 값들에 기반하여 Unknown 상태의 격자점에서의 도착시간을 계산한다. 먼저, 초기 경계면을 이용하여 격자점들을 초기화한다. 이때, 경계면이 안쪽에서 바깥쪽으로 이동한다면 경계면의 바깥쪽 격자점은 큰 값으로(∞), 경계면의 안쪽의 격자점은 0으로 초기화 한다. 경계면 안쪽의 격자점들은 이미 계산이 끝난 것으로 간주하여 Accepted 상태를 갖고, 바깥쪽 격자점은 Unknown 상태를 갖는다. 다음으로는 0으로 초기화된 격자점(Accepted)에 인접한 격자점들의 도착시간을 Eq. 6을 이용하여 계산한다. 이때, 인근의 Accepted 상태를 갖는 격자점 정보만 이용한다. 계산된 격자점들은 In-band 상태를 갖는다.

Fig. 3.

Fast marching algorithm.

초기화 과정이 끝나면 In-band 상태를 갖는 격자점에서부터 Iteration을 시작한다. Iteration 과정은 다음과 같다. 첫째, In-band 상태의 격자점들에서 최솟값을 찾고, 해당 격자점의 상태를 Accepted로 갱신한다. 둘째, 해당 격자점의 인근에서 Unknown 상태의 격자점들의 값을 Eq. 10, 11을 이용하여 계산한다. 이때, 계산되는 격자점 인근의 In-band나 Accepted 상태의 격자점 값들만을 이용한다. 셋째, 계산된 격자점들의 상태를 In-band로 갱신한다. 이와 같은 Iteration을 모든 격자점이 Accepted 상태가 될 때까지 반복한다.

이를 Matlab 환경에서 구현하였다. 이같이 구현된 FMM의 경우, 이동경계면 문제를 N×N의 2차원 직교 격자를 이용하여 해석할 때 N2logN의 계산 복잡도를 갖게 된다.

2.3 연소관 형상을 반영한 연소면적 계산

2.1과 2.2절에서는 각각 LSM과 FMM을 적용하여 이동경계면을 해석하고, 시간에 따른 형태와 위치를 표현하는 격자를 얻어내는 방법에 대해 다루었다. 2.3절에서는 격자를 후처리하여 경계면의 형상으로부터 직접 연소면적을 계산하는 방법을 요약하였다.

연소면적계산을 위해서는 경계면을 표현하는 격자점만을 모두 찾아서 적분해주면 된다. FFM의 경우에는 원하는 시간 값을 갖는 격자점을 찾고, LSM의 경우에는 SDF가 0 값을 갖는 점들을 찾으면 된다. 본 논문에서는 디락 델타(dirac delta) 함수를 통해서 원하는 영역의 격자점만을 선택하였다[16], 적용한 함수는 다음식과 같다.

Eq. 12는 FMM의 경우에 적용된 디락 델타 함수를 나타내고, Eq. 13은 LSM에 적용된 디락 델타 함수를 나타낸다. Eq. 12에서 T는 도착시간을 의미하고, t는 경계면의 형태와 면적을 알고 싶은 특정 시점의 시간을 의미한다. Eq. 13에서 𝜙는 SDF 값을 나타낸다. 두 식에서 𝜎는 사용자 정의 변수로, 본 논문에서는 LSM과 FMM에서 각각 △x, △x/V값으로 설정하였다. Fig. 4a에 적용된 디락 델타 함수의 예시가 그려져 있다.

Fig. 4.

Example of Dirac delta function(a) and Heavy-sided function(b).

한편, 연소가 진행되는 중에는 웹 두께가 얇은 부분의 연소 경계면이 다른 부분보다 먼저 연소관 벽까지 진행될 수 있다. 이 부분에서는 웹 두께가 0이 되기 때문에 유량이 발생하지 않는다. 따라서 연소관 벽 바깥까지 진행된 이동경계면의 부분은 연소면적 계산에서 제외해야 한다. 본 논문에서는 다음 식과 같이 헤비사이드(heavy-sided) 함수를 이용하여 연소관 벽 외부의 격자점들을 계산에서 제외하였다.

Eq. 14에서 x,y는 각각 격자점들의 x 방향 및 y 방향 좌표를 의미하고, casing은 연소관 내부의 격자점 집합을 의미한다. Fig. 4b에 적용된 헤비사이드 함수의 예시가 그려져 있다.

최종적으로 연소면적 계산은 다음과 같이 가능하다.

Eq. 15와 16은 각각 FMM에 적용할 경우와 LSM에 적용할 경우의 연소면적 계산식이다.

3. 다양한 형상에 대한 연소면적 해석

본 논문의 목적은 FMM을 적용한 연소면적 해석 프로그램을 개발하고 계산시간의 개선 정도를 LSM 기법 기반 프로그램과 비교하는 것이다. LSM과 FMM의 계산시간은 공통으로 격자 크기에 의존한다. 그러나 LSM의 경우, 격자 크기만이 아니라 시간 간격에도 의존한다. 이 때문에 단순히 동일 격자 조건을 설정한 후 계산시간을 비교할 수 없다. 따라서 이 장에서는 격자 해상도와 시간 간격에 따른 해석 결과의 정확도를 확인하고, 두 기법이 동등한 수준의 계산 정확도를 갖는 격자 크기와 시간 간격 조건을 설정한 후 계산시간을 비교하고자 하였다.

계산 조건 설정 후에는 여러 형상에 대하여 프로그램의 정확도를 검증하고, 계산시간을 비교하기 위해 성형(star) 그레인과 슬롯형(slot) 그레인에 대한 해석을 진행하고 결과를 비교하였다. 두 그레인 형상은 이미 많은 연구가 진행된 표준 형상이다. 또한, 연소가 진행됨에 따라 기하적 위상 변화가 여러 단계에 걸쳐 나타남에도 연소면적을 간단한 수식을 통해 계산할 수 있어 프로그램의 검증에 적합하다[3,5]. 이 장에서 사용된 그레인 형상은 모두 임의의 숫자를 기반으로 정의된 것이며, 단위는 없다.

3.1 오차 비교를 통한 계산 조건 설정

계산 정확도는 격자의 해상도, 시간 간격만이 아니라 이동경계면의 초기 형상에 따라서도 달라질 수 있다. 이 때문에 여러 선행연구에서는 복잡한 형상에 해석 기법을 적용하기 전에 점이나 원, 사각형과 같은 간단한 기준형상을 이용하여 기법의 정확도를 분석하였다[2,14]. 특히, 참고문헌 [2]에서는 경계면이 확장하며 곡률이 커지는 곡률 효과(rounding effect)가 나타나는 부분에서 오차가 커짐을 확인한 바 있다. 본 논문에서도 곡률 효과를 관찰 가능하며, 구현이 간단한 사각형 형상을 기준형상으로 선정했다. 그리고 여기서 나타나는 오차를 기준으로 계산 조건을 설정하였다. Fig. 5는 같은 격자 조건에서 LSM과 FMM을 사각형 형상에 적용하여 계산한 결과를 나타낸다. 그림에서 두 기법을 통해 얻은 수치해가 모두 해석해와 유사함을 확인할 수 있다. 그러나 선행연구와 마찬가지로 두 기법 모두 곡률 효과가 나타나는 부분에서 오차가 나타난다. 또한, 해당 부분에서 두 기법으로 얻어진 경계면 형상 사이에도 차이가 있음을 확인할 수 있다.

Fig. 5.

Comparison of rectangle shaped front propagation results.

한편, LSM과 FMM을 이용하면 시간에 따른 면적 변화 값을 얻을 수 있고, 경계면의 모든 점에서 속도 값이 일정하므로 이를 이동 거리 혹은 연소거리에 따른 값으로 나타낼 수 있다. Fig. 6은 연소거리에 따른 면적과 오차를 나타낸다. 그림으로부터 같은 격자 조건에서도 두 기법의 연소거리 별 오차 양상이 서로 다름을 확인할 수 있다. LSM의 경우, 연소거리가 늘어남에 따라 오차가 0 근처에서 시작하여 완만히 증가한다. 반면, FMM의 경우에는 연소거리 0 근처에서 큰 오차가 나타난 이후 작아졌다가 다시 완만히 증가함을 알 수 있다. 이처럼 서로 다른 오차 양상에는 크게 두 가지 원인이 있다. 첫째로는 앞서 설명한 바와 같이 곡률 효과가 나는 부분에서 두 기법을 통해 계산된 경계면의 형상이 서로 다르다는 점이다. 또한, 두 기법의 격자점 초기화 방식의 차이도 영향을 미친다. 이에 관한 내용은 3.2절에 자세하게 설명되어있다. 본 논문에서는 이처럼 서로 다른 양상을 보이는 오차를 모두 반영할 수 있도록 적분 오차(Integral of Absolute Error, IAE) 값을 기준으로 두 기법의 정확도를 비교하였다. 오차 적분 값은 다음과 같이 계산하였다.

Eq. 17에서 d는 경계면 이동 거리, Aanalytic, An은 각각 해석적 방식으로 구한 면적과 수치적 방식으로 구한 면적을 의미한다.

Fig. 6.

Calculated burn area & error.

앞서 설명한 것처럼 시간에 따른 면적계산 결과는 연소거리에 따른 함수로 활용하게 되므로, 경계면의 속도는 임의로 설정할 수 있다. 따라서 본 논문에서는 경계면 속도를 0.1로 고정하고, 여러 격자 크기와 시간 간격 조건에서 계산을 수행했다. Fig. 7은 여러 CFL 값에서 계산한 오차 적분 값을 비교한 것이다. LSM의 결과는 점선으로 표기하였으며, FMM의 경우에는 실선으로 표시했다. 서로 다른 격자 크기 조건은 색을 통해 구분하였다. FMM의 경우, 격자 크기에 의해서만 오차가 결정되기 때문에, CFL값에 따른 오차 적분 값 그래프가 상수로 나타나게 된다. 반면, LSM의 경우, 모든 격자 크기 조건에서 CFL값이 커짐에 따라 약 0.7까지는 오차 적분 값이 감소하다가, 0.7을 초과하면 경계면이 이동함에 따라 계산 결과가 불안정해지면서 급격하게 증가하는 경향을 보였다. 한편, 두 기법 모두 격자 크기가 줄어듦에 따라서 오차 적분 값이 감소하는 경향을 보였다. 또한, CFL값 0.7 근처에서 두 기법의 오차 적분 값이 비슷한 것을 확인할 수 있다. 따라서 본 논문에서는 CFL값이 0.7이 되도록 속도, 시간 간격, 격자점 간격을 설정했다.

Fig. 7.

IAE under various calculation conditions.

한편, 오차 적분 값을 통해서는 여러 계산 조건에 따른 정확도 양상을 확인할 수 있지만, 수치해가 해석해에 얼마나 근접한 결과를 제시하는지는 확인할 수 없다. 따라서 적절한 계산 조건 설정을 위해서는 추가적인 기준이 필요하다. 고체 추진기관의 설계에서 주요하게 고려되는 변수 중 하나는 추력을 연소시간에 따라 적분한 총 역적 값이다[1]. 고체 추진기관의 순간 추력과 연소면적은 비례관계에 있으므로 연소면적의 적분 값 오차는 총 역적 값에 직접적으로 영향을 미친다. 따라서 본 논문에서는 두 번째 기준으로 오차 적분 값을 해석적 면적의 적분 값으로 나눈 오차 비율(IAE ratio)을 이용하였다. 오차 비율은 Eq. 18과 같이 계산할 수 있다.

Eq. 18에서 RIAE는 오차 비율을 나타낸다. Table 1에 CFL값 0.7을 만족하는 계산 조건과 각각의 기법의 결과가 요약되어있다. 본 논문에서는 CFL값이 0.7이고, 두 기법 모두 오차 비율 0.5% 미만을 만족하는 경계면 이동 속도 0.1, 격자 크기 0.005, 시간 간격 0.035초를 계산 조건으로 설정하였다.

3.2 연소면적 해석 결과

Fig. 8과 Fig. 9에 각각 본 연구에 이용된 성형 그레인과 슬롯형 그레인의 초기 형상 정의가 나타나 있다. Table 2에는 본 연구에서 사용한 형상 매개변수가 요약되어있다. 사용된 숫자는 선행연구[3,5]에서 제시한 값을 참고하여 임의로 설정하였으며, 단위는 없다. 슬롯형 그레인의 경우 초기 형상 변수에 따라서 12가지의 연소 시나리오가 나타난다. 본 논문에서는 이중 총 네 단계에 걸쳐서 연소가 진행되는 A형과 B형 슬롯 형상을 참고했다[3].

Fig. 8.

Star grain configuration.

Fig. 9.

Slot grain configuration.

성형 그레인의 경우 연소면적이 일정하게 유지되는 1단계, 점진적으로 증가하는 2단계, 급격하게 작아지는 3단계에 걸쳐서 연소면적이 변화한다. 특히 2단계에서 3단계로 진행될 때는 가장 얇은 웹의 두께가 0이 되면서 불연속적인 연소면적 변화가 나타난다. Fig. 10은 성형 그레인에 LSM과 FMM을 적용한 연소면적 해석 결과를 나타낸다. 두 기법 모두 해석해와 유사한 결과를 얻을 수 있음을 확인할 수 있다.

Fig. 10.

Star grain front propagation result.

Fig. 11은 연소거리에 따른 연소면적 값을 나타낸다. 전체적으로 수치적 연소면적 계산 값이 해석적 계산 값에 근접함을 알 수 있다. 그러나, 오차가 상대적으로 크게 나타나는 두 구간을 확인할 수 있다. Fig. 12는 두 구간에서의 연소면적 그래프를 확대한 것이다. 왼쪽 그림은 연소거리 0 근처에서 두 기법을 이용한 계산 결과를 나타낸다. FMM을 적용한 연소면적 계산 결과에서만 약 15 이상의 오차가 발생함을 확인할 수 있다. 오른쪽 그림은 연소면적 변화 2단계에서 3단계로 진행되는 연소거리 1.9 근처를 확대한 그림이다. 두 기법 모두 최대 약 3.5 정도의 오차를 확인할 수 있다.

Fig. 11.

Star grain burn area calculation.

Fig. 12.

Error in star grain burn area calculation.

해당 구간에서 오차가 커진 원인은 각 기법의 격자점 초기화 방식과 후처리 과정에서 기인한다[16]. Fig. 13은 LSM을 1-D 이동경계면 문제에 적용할 경우, 경계면의 초기위치(Fig. 13a)와 경계면이 연소관 위치까지 이동한 말기 위치(Fig. 13b)에서의 연소면적 계산예시를 나타낸다. Fig. 13a에는 이동경계면의 초기위치와 면적계산을 위해 디락 델타 함수가 적용된 모습이 그려져 있다. 연소관의 형상을 반영하기 위해 모든 격자점에 헤비사이드 함수가 적용되어 있으므로 연소관 벽 바깥쪽의 격자점들의 값은 0이다. 또한, 그림에서 디락 델타 함수나 헤비사이드 함수에 의해 0이 되지 않고 면적계산에 참여하는 격자점들은 파란색으로, 계산에 참여하지 않는 격자점들은 하얀색 점으로 표시되어있다. Fig. 13a에서 경계면의 이동 초기에는 디락 델타 함수의 영역 안에 있는 모든 격자점이 연소면적 계산에 참여하는 것을 확인할 수 있다.

Fig. 13.

Example of LSM applied area calculation at initial(a) and final(b) front position.

반면, Fig. 13b에서 확인할 수 있는 것처럼 경계면이 연소관에 접촉하는 연소 말기에는 연소관 바깥쪽 일부 격자점들이 헤비사이드 함수에 의해 0이 되므로 연소면적 계산에 참여하지 못하게 된다. 이 때문에 LSM을 통한 연소면적 계산에서는 경계면이 연소관에 접촉하기 시작하는 3단계에서 오차가 발생한 것으로 보인다.

한편, Fig. 14는 FMM을 1-D 이동경계면 문제에 적용할 경우, 경계면의 초기위치(Fig. 14a)와 말기 위치(Fig. 14b)에서의 연소면적 계산예시를 나타낸다. Fig. 13과 마찬가지로 연소면적 계산에 참여하지 못하는 격자점들은 하얀색으로, 참여하는 격자점들은 초록색으로 표시되어있다. LSM의 경우에는 각 격자점이 연속적인 SDF 값을 갖도록 초기화되기 때문에 디락 델타 함수의 영역에 있는 모든 격자점 연소면적 계산에 참여한다(Fig. 13a). 반면, FMM의 경우에는 경계면 안쪽의 격자점들이 모두 0으로 초기화되므로 초기 경계면의 연소면적 계산에 참여하지 못하는 격자점들이 있음을 알 수 있다(Fig. 14a). 또한, 연소관 형상 반영을 위해 헤비사이드 함수가 적용되어 있으므로 연소관 바깥쪽 점들은 연소 말기의 연소면적 계산에 참여하지 못하게 된다(Fig. 14b). 따라서 FMM의 경우에는 연소거리 0 근처, 그리고 경계면이 연소관에 접촉하는 4단계에서 큰 오차가 발생하였다.

Fig. 14.

Example of FMM applied area calculation at initial(a) and final(b) front position.

Fig. 15와 Fig. 16은 각각 A형 슬롯 그레인의 연소면적 해석 결과와 연소면적 계산 결과를 나타낸다. 또한, Fig. 17과 Fig. 18은 각각 B형 슬롯 그레인을 해석한 결과와 연소면적 계산 결과를 나타낸다. 해당 그림들에서 두 기법을 통해 얻은 수치해가 해석해에 근접함을 확인할 수 있다. 한편, 슬롯형 그레인의 연소면적 계산 결과에서도 경계면이 연소관에 접촉하는 시점부터 오차가 상대적으로 커짐을 확인할 수 있다. A형 슬롯 그레인의 경우 3단계에서, B형 슬롯 그레인의 경우에는 4단계에서부터 접촉이 시작된다. 따라서 Fig. 16과 Fig. 18에서 해당 구간에서 오차가 커짐을 확인할 수 있다. 또한, FMM의 경우에는 연소거리 0 근방에서 오차가 발생함을 확인할 수 있는데, 이러한 오차는 앞서 설명한 성형 그레인에서 나타난 오차와 같은 원인을 갖는 것으로 보인다. Table 3에 그레인의 종류와 기법별 계산 결과가 요약되어있다. 오차 비율이 모든 결과에서 1% 이하임을 확인할 수 있다.

Fig. 15.

Slot type A grain front propagation result.

Fig. 16.

Slot type A grain burn area calculation.

Fig. 17.

Slot type B grain front propagation result.

Fig. 18.

Slot type B grain burn area calculation.

3.3 계산시간 비교

Fig. 19는 LSM과 FMM을 이용한 프로그램의 흐름도를 나타낸다. 두 프로그램 모두 크게 두 단계로 계산을 수행한다. 초기화 단계에서는 초기 경계면 형상을 반영하여 격자점을 초기화한다. 연소면적 계산단계에서는 경계면의 이동과 후처리를 진행한다. 오일러 접근법을 이용한 연소면적 해석 프로그램의 경우, 많은 격자점 개수로 인하여 초기화 단계에 많은 시간이 소요된다. 따라서 본 논문에서는 총 계산시간 이외에도 초기화 단계와 연소면적 계산단계에서의 계산시간을 구분하여 측정하였다. LSM의 경우에는 초기 SDF 값으로 격자점을 초기화한다. 따라서 모든 격자점이 경계면의 안쪽이나 바깥에 위치하는지 검사하는 과정과 해당 격자점에서 경계면까지의 최소거리를 구하는 과정을 거친다. 격자가 초기화되면 Eq. 2에 따라 경계면의 이동을 계산하며, 계산이 완료될 때까지 모든 시간 간격에서 경계면 이동 계산, 후처리 과정을 실행한다. 반면, FMM의 경우에는 격자점이 초기 경계면의 안쪽이나 바깥에 위치하는지만 검사하면 초기화가 완료된다. 초기화 이후에는 2.2절의 알고리즘에 따라 계산이 진행되며, 모든 격자점의 계산이 완료되면 사용자가 설정한 시간 간격에 따라 후처리를 진행한다. 이때, 두 프로그램이 같은 시간 간격의 데이터를 출력할 수 있도록 설정하여 같은 횟수의 후처리 과정을 수행하였다.

Fig. 19.

Flow chart of LSM and FMM applied burn-back analysis program.

Table 3에 각 그레인 형상에 대한 계산시간이 요약되어있다. 모든 경우에서 FMM의 초기화 시간이 LSM의 30~55% 수준임을 확인할 수 있다. 이는 FMM의 초기화 과정이 LSM보다 더 간단하기 때문이다. 최대 연소거리가 각각 0.9, 1.4인 슬롯 A, B형 그레인의 경우에 두 기법이 서로 유사한 연소면적 계산시간을 보였다. 그러나 최대 연소거리가 2 이상이었던 성형 그레인의 경우에서는 FMM의 계산시간이 LSM에 비해 약 70% 정도 수준임을 확인할 수 있다. 또한, 모든 경우에서 FMM의 총 계산시간이 LSM에 비해 50~65% 수준임을 확인하였다.

4. 결 론

본 논문의 목적은 균일, 등방성 가정이 가능한 추진제에 대하여 Fast Marching Method(FMM)를 적용한 연소면적 해석 프로그램을 개발하고 계산시간의 개선 정도를 Level Set Method(LSM) 기반 프로그램과 비교하는 것이다. 동등한 조건에서의 계산시간 비교를 위해 사각형 형상의 기준형상을 도입하고, 적분 오차(IAE)와 오차 비율(IAE ratio)를 기준으로 동등한 수준의 정확도를 갖는 계산 조건을 설정하였다. 이후 다양한 그레인 형상에 대해 해석을 수행한 결과 모든 그레인에서 두 기법의 오차 비율이 1% 이하인 것을 확인하였다. 또한, FMM을 이용한 프로그램의 총 계산시간이 LSM 대비 50~65% 수준임을 확인하였다. 본 연구에서 개발한 FMM 기반의 연소면적 해석 프로그램은 균일, 등방성 가정이 가능한 추진제에 대하여 유효하며, 이를 이질성, 이방성 등이 나타나는 추진제에 적용하기 위해서는 추가 연구가 필요하다.