1. 서 론
최근 나로호가 성공적으로 발사되고, 정보산업의 발달로 인해 인공위성의 수요가 증가하고 있다. 이에 따라 지구 궤도로 올리기 위한 로켓 추진기관의 수요도 증가하고 있다. 이 중, 고체연료를 사용하는 로켓의 경우, 다른 추진기관에 비해 구조가 단순하고, 높은 추력을 발생시키며, 사용 및 취급이 쉽다는 장점이 있다. 그러나 고온·고압의 극한 환경에서 작동하며, 연소 중단·제어가 어려운 단점도 존재한다. 이 때문에, 로켓의 작동 환경이 변하거나 비교적 사소한 결함일지라도 발생 시 여러 가지 고장모드가 발생하게 되고, 로켓 전체의 소실이 발생할 수 있다. 이러한 고장이 발생하는 경우, 재개발 및 기타 비용으로 인해 많은 손실이 발생하기 때문에 고체로켓 모터에 대한 신뢰도예측은 필수적이다.
이 때문에 전 세계적으로 2~30년 전부터 여러 가지 기법을 적용하여 고체로켓의 신뢰도예측 연구를 진행하였다. Cruse 등(1990)은 우주 추진체 부품(터보펌프, 터빈의 회전날개)의 응력파괴에 관한 확률적 구조해석을 위해서 mean value first-order second moment(MVFOSM)과 advanced- MVFOSM을 통해 신뢰도지수(reliability index)를 탐색하여 신뢰도를 예측하는 기법을 제시하였다. Zheng(1993)에서는 로켓 모터 케이스의 응력파괴, 균열파괴에 대해 역시 MVFOSM 기법을 적용하여 신뢰도예측을 수행하였다. 그러나 MVFOSM은 평균값을 이용하여 신뢰도지수를 계산하기 때문에 고장확률을 계산할 때 정확하지 못한 단점이 있다.
Bozkaya 등(2005)에서는 고체로켓 모터의 신뢰도예측을 위해 여러 가지 고장함수(응력파괴, 균열파괴, total impulse 미달 등)에 대해 반응표면을 구축하고, Monte Carlo Simulation(MCS) 방법으로 고장확률을 구하였다. MCS를 이용하여 고장확률을 계산하는 경우, 정확한 결과를 얻을 수 있다는 장점이 있지만, 계산을 할 때마다 결과값이 달라진다. 또한, 신뢰도 높은 결과를 얻기 위해서는 105~106회 이상의 반복 수행해야 하기 때문에 고장함수가 복잡하거나 유한요소해석과 같은 해석이 포함된 경우, 많은 계산 시간이 소요되어 현실적으로 불가능하다는 단점이 있다.
이러한 문제를 해결하기 위해서 최근에는 여러 가지 기법을 적용하는 연구가 많이 진행되고 있는데, 신뢰도예측에서 널리 사용되는 first-order reliability method(FORM)을 이용한 연구가 많이 진행되었다. Yang(1997)은 고체로켓 모터의 gas generator 시스템의 고장모드(점화 실패)의 신뢰도예측을 수행하였다. 이를 위해, Generator 내부의 압력 및 온도 등의 물성치 변화를 예측하는 모델링을 수행하고, 이를 통해 한계상태함수를 도출하였다. 그리고 FORM을 적용하여 신뢰도지수를 도출하고, 이를 통해 gas generator의 신뢰도예측을 수행하였다. Kuran(2007)은 3가지 고장모드(total impulse 미달, 로켓의 탈출 속도 미달, 최대 가속도 미달)에 대해 고장확률을 구하였다. 이를 위해 고장모드에 FORM, 반응표면기법을 하나씩 적용하여 각 기법별 고장확률을 계산하고 MCS와 비교하였다. 그러나 FORM은 비선형성이 강한 문제나 고장확률이 매우 낮은 사례에는 사용하기 어렵다는 단점이 존재한다.
또한, 이외의 다른 기법들을 사용하여 신뢰도 예측을 수행한 사례도 존재한다. Morio(2011)는 중요도 추출법의 일종인 non-parametric adaptive importance sampling(NAIS)을 이용하여 고체로켓 부스터의 신뢰도예측을 수행하였다. 이를 위해 launch impact scattering에 대해 한계상태방정식을 정의하고, NAIS를 이용하여 impact scattering의 발생확률을 계산하였다. 또한, 이 결과를 MCS의 결과와 비교하여 보았다. Balesdent(2013)는 위의 Morio(2011)와 동일하게 launch impact scattering에 대해서 kriging surrogate model을 기반으로 한 중요도 추출법을 적용하여 신뢰도예측을 수행하였다. 이러한 중요도 추출법은 고장이 발생하는 영역을 중점으로 하여 샘플링을 진행하는 기법이기 때문에 MCS에 비해 적은 수의 계산으로도 정확도가 높은 결과를 얻을 수 있다. 그러나 중요도 추출법을 적용하는 영역 선정 시, 잘못된 영역을 선택하게 되면 전혀 다른 계산 결과를 얻을 수도 있다.
위에서 제시된 연구의 단점을 보완하고자 본 연구에서는 고체로켓이 연소하면서 발생하는 내압으로 인한 고장확률을 계산하기 위한 새로운 방법을 제시하였다. 이를 위해 내압으로 인한 3가지 고장모드(응력파괴, 균열파괴, 볼트 체결 부 파손)를 정의하고, 확률기반 연소해석을 통해 maximum expected operating pressure(MEOP)의 확률분포를 구하여 케이스의 고장모드 계산에 적용하였다. 또한 신뢰도예측 정확도 향상을 위해 FORM과 most probable failure point(MPP)를 중심으로 한 반응표면기법, 그리고 중요도 추출법을 적용하여 보다 정확하고 효율적인 신뢰도예측을 수행하는 새로운 방법을 제시하였다.
2. 본 론
Fig. 1은 본 연구에서 진행된 신뢰도예측 Process를 나타낸다. 첫째로 고체로켓 모터 내에서 발생하는 연소과정을 모델링하여 MEOP의 확률분포를 예측하고, 이를 응력해석과 신뢰도예측에 입력변수로 사용한다. 다음 단계에서는 고체로켓 모터의 유한요소 모델링을 통해 응력해석을 수행하며, FORM을 이용하여 각각의 입력변수들의 MPP를 계산한다. 마지막으로 이 점을 중심으로 반응표면모델을 구축하고, 중요도 추출법을 통해 고체로켓 모터의 고장확률을 계산한다. 그리고 본 연구에서 제시한 기법과 FORM, crude MCS의 결과와 비교 분석한다.
2.1 고체로켓 모터의 MEOP 예측
이 장에서는 앞서 언급한 신뢰도예측 절차 중 고장모드에 가장 영향력이 큰 MEOP를 계산하고, 확률분포를 예측하는 과정에 대해 상세히 설명한다.
2.1.1 고체로켓 모터 연소해석
MEOP는 고체로켓의 추진 시, 시간에 따라 발생하는 고체로켓 연소관의 압력 중 최대값을 의미하며 이를 계산하기 위해서는 아래 지배방정식으로 주어지는 고체로켓의 연소해석을 수행해야 한다.
여기서 Pc는 연소관 압력, 
는 추진제 질량 유량, 
는 후퇴율(regression rate), ρp는 추진제 밀도, Ab는 추진제 연소 면적, CD는 항력계수, At는 노즐 목 면적, Tc는 연소관 내부 온도, M은 추진제 몰 질량, Vc는 연소가스의 속도를 의미한다. 식 (2)에서 a는 후퇴율 상수(burning rate constant), n은 후퇴율 지수(burning rate exponent)를 의미한다. 본 연구에서 선정한 고체로켓은 단일 port 추진제를 사용하는 것으로 가정하였으며, 이상 로켓 조건으로 가정하여 연소해석을 진행하였다. 이를 통해 위의 지배방정식을 이용하여 시간에 따른 연소관 내부의 압력을 계산하고, MEOP를 도출한다. 자세한 내용은 Sutton 등(2010)을 참고하기 바란다.
2.1.2 확률변수 선정 및 MEOP 확률분포 계산
MEOP를 계산하기 위해서는 식 (1)의 해를 구해야 하는데 이 과정에서 많은 변수들이 사용된다. 본 연구에서는 Table 1과 같이 8개의 입력변수를 선정하고 해석을 수행하였다. 이때 변수들은 계산의 편의를 위해 모두 정규분포로 가정하였으며, Bozkaya 등(2005)과 Kim 등(2012)에서 제시된 값에 따라 COV를 표와 같이 정하였다.
위의 변수를 통해 연소해석을 한 결과를 Fig. 2에 나타내었다. 그래프에서 압력이 최대가 될 때의 값이 MEOP이며 그 값은 약 12MPa이다.
Table 1
Input parameter of MEOP calculation
Table 1에서 선정된 확률변수를 이용하여 MCS, 즉 많은 수의 입력변수 샘플을 생성하고 각 케이스마다 MEOP를 구하면 확률분포를 얻을 수 있다. 그러나 고체로켓의 경우, 높은 신뢰도(99%이상의 신뢰수준)의 결과를 얻기 위해서는 105회 이상의 반복 계산이 필요하다. 본 연구에서는 이러한 계산비용을 줄이기 위해 민감도 해석을 수행하여 영향력이 큰 입력변수를 선별한 후, 이들에 관한 반응표면을 다항식으로 구축하여 압력을 계산하였다.
입력변수의 평균과 COV가 정해지면 민감도 해석을 수행하고 이를 이용하여 각각 민감도 분석을 수행하여 MEOP 결과에 영향력이 큰 변수를 선별한다. 민감도 해석을 위해 각 변수를 +σ만큼 변화시켜서 민감도(입력변수에 따라 결과 값이 변화하는 정도)를 구하고, 그 결과가 5% 이상인 변수를 선별한 결과 Fig. 3에서 보는 바와 같이 3개의 변수(burning rate constant, burning rate exponent, throat radius)가 선별되었다.
선별된 변수를 이용하여 반응표면을 구축하기 위해 이들의 평균과 COV 값을 이용하여 latin hypercube sampling (LHS)을 이용하여 각 변수의 실험점을 생성하였다. 이 점은 변수의 평균을 기준으로 ±3σ의 범위 내에서 생성하며, 샘플 개수는 변수 3개의 2차 다항식의 계수 개수와 동일하게 10개를 생성한다. 샘플 생성이 완료되면 MEOP에 대해 2차 다항식으로 근사화하여 반응표면을 구축하고, 반응표면의 검증을 위해 최소자승법(least square method)을 통해 결정계수 R2을 계산한다. 계산 시, 반응표면의 계수의 수와 동일하게 10개의 LHS 실험점이 적용되었다. Fig. 4에 반응표면을 이용하여 구한 MEOP의 확률분포를 표시하였으며, 그 결과 평균, 표준편차, 왜도 및 첨도가 각각 12.57, 3.21, 0.6538 및 3.5934로 나타나 대칭적이지 않은 분포임을 알 수 있었다.
2.2 고체로켓 모터 케이스의 응력해석
MEOP의 확률 분포가 계산되면, 이를 이용하여 로켓 모터 케이스에 대해서 응력해석을 수행한다. 이를 위한 고체로켓의 형상 변수는 Table 2와 같다.
Table 2
Structural parameter
이를 이용하여 연소관, 점화기, 토출관, 노즐, 체결 볼트를 모델링하였고, 재질은 AISI 4130 Steel 재질로 설정하였다. 입력변수는 계산의 편의를 위해 모든 변수를 정규분포로 가정하였다(Bozkaya et al, 2005; Kim et al. 2012). 여기서 볼트체결을 위한 조임 길이 δ는 케이스 체결을 위해 볼트에 예력을 부여하기 위해 가하는 조임 길이를 의미한다. 볼트의 예력을 입력변수로 정하는 경우에는 이를 조임 길이로 변환해야 하며, 이는 Budynas 등(2011)을 참고한다. 한편 문헌 Choi 등(2010)에 의하면 로켓 모터 케이스는 2D 축대칭 구조물로 모델링하는 것이 가능함을 입증하여 본 연구에서도 2D로 가정하였다.
유한요소해석 시 연소관-볼트 체결 부는 겹쳐있는 node를 하나로 통합하여 bonding하는 tied 접촉조건을, 연소관-점화기, 볼트-점화기, 노즐-볼트, 노즐-연소관 연결 부위에는 접촉조건을 부여하였고, 이를 Fig. 5에 나타내었다. 붉은 선(굵은선)으로 표시된 부분은 접촉조건을, 노란 선으로 표시된 부분은 tied 접촉이 적용된 곳을 의미한다.
Fig. 6는 각 변수의 평균을 이용하여 로켓 모터 케이스의 von-Mises 응력을 구한 결과이다. 유한요소해석이 수행되면 케이스와 노즐, 볼트에 가해지는 von-Mises 응력과 주응력의 최대값을 계산하는데 von-Mises 응력은 케이스의 응력파괴와 볼트파손의 고장확률을, 주응력은 케이스 균열파괴의 고장확률을 계산하는데 사용된다. 고체로켓 모터 케이스의 최대 응력은 내부의 fillet에서 발생한다. 이 곳에서의 응력의 크기는 최대 482.1MPa로 이는 벽면에 가해지는 hoop stress(σ=Pr/t=153.4MPa)의 약 3배 정도의 크기를 갖는 것을 알 수 있다.
2.3 로켓 모터의 신뢰도 해석
2.3.1 고장모드 정의
유한요소해석을 수행한 후, 각각의 고장모드에 대해 신뢰도예측을 수행한다. 신뢰도예측에 앞서 3가지 고장모드(응력파괴, 균열파괴, 볼트파손)에 대해 한계상태방정식을 정의한다.
식 (3), (4), (5)는 각각 응력파괴, 균열파괴, 볼트파손의 한계상태방정식을 나타낸다. 식 (3)에서 P는 MEOP, E는 케이스 재질의 탄성계수, S는 Table 2의 형상변수를 의미한다. 식 (4)에서 KIC는 파괴 인성치, KI는 응력확대계수, α는 균열크기, β는 균열형상에 따라 정해지는 계수를 의미한다. 식 (5)에서 Sp는 볼트의 보증강도, σb는 볼트에 가해지는 응력, Fi는 체결 토크로 인한 볼트의 예력을 나타낸다. 이 중 균열파괴의 경우, 주응력 방향으로 반월형 균열이 발생한다고 가정하여 한계상태방정식을 정의한다. 균열이 반원형일 때, 응력확대계수 KI는 다음과 같이 계산한다(Yang et al, 1993).
여기서, t,R은 케이스의 두께 및 반경을 나타낸다. 각각의 고장함수 값이 0보다 작은 경우 고장이 발생한다.
2.3.2 유한요소해석을 적용한 FORM 수행
각 고장모드의 한계상태방정식이 정의되면 이를 이용하여 FORM을 이용하여 MPP를 탐색하고 신뢰도지수를 계산한다. FORM을 이용한 MPP 탐색 개념도를 Fig. 7에 보였다. 자세한 수행 과정은 Mahadevan 등(2000)을 참조한다.
여기서 MPP는 고장함수가 0이 되는 곡선 중, 확률 분포의 평균과 가장 가까운 지점의 변수 값을 의미한다. 신뢰도지수 β는 정규분포로 근사화된 고장확률을 계산할 수 있는데 일반적으로 비선형성이 높지 않거나 COV가 작은 경우에는 신뢰도지수만으로도 어느 정도 정확한 고장확률을 예측할 수 있다. 그러나 매우 높은 신뢰도를 요구하는 경우에는 정확도를 높이기 위해 추가 과정이 필요하다. 이를 위해서 반응표면기법과 중요도 추출법을 적용하였는데, 이는 다음 장에서 설명하기로 한다.
한편 FORM을 적용할 때, MEOP와 같이 정규분포가 아닌 입력변수에 대해서는 아래의 식을 통해 등가정규분포로 변환하여 계산을 수행해야 한다.
식 (7)에서 
, 
은 각각 등가정규분포의 평균과 표준편차, 
는 i번째 iteration에서의 변수 값을 의미한다.
2.3.3 반응표면 구축
MPP가 계산되면 이 점을 중심으로 반응표면을 구축한다. 반응표면 구축 시, MEOP의 반응표면을 생성하였던 과정과 동일한 과정을 수행한다. 각각의 변수에 대해서 평균을 기준으로 하여 ±3σ의 범위에서 LHS를 수행하여 샘플을 생성하고, 이를 통해 반응표면을 구축한다. 반응표면은 2차 다항식으로 구성한다. 이 때, 앞의 과정과는 다르게 고장모드별로 입력변수의 민감도가 달라질 수 있기 때문에 민감도 분석 통한 변수선별은 하지 않는다. 반응표면을 구축하기 위해 각각의 고장모드별로 총 10, 15, 15개의 LHS 실험점이 적용되었다.
2.3.4 중요도 추출법 및 고장확률 계산
반응표면이 구축되면 중요도 추출법에 의해 고장확률을 계산한다. 중요도 추출법의 개념도는 Fig. 8과 같다. 여기서 p(X)는 원래의 확률 분포, q(X)는 새롭게 가정된 proposal PDF를 나타내며 이는 MPP를 중심으로 주어진다. q(X)를 이용하여 샘플을 추출하면 다음과 같은 식에 의해 고장확률을 계산할 수 있다.
이런 방법으로 고장확률을 계산하게 되면 원래의 확률분포로 샘플 생성할 경우 보다 훨씬 적은 개수로 계산이 가능하여 계산비용을 크게 절감할 수 있다.
2.4 고장확률 계산 결과
앞서 제시한 방법을 통해 로켓 모터의 고장확률을 계산하였다. 이렇게 계산된 고장확률의 정확도를 FORM에서 계산된 신뢰도지수 기반 고장확률과 반응표면 및 중요도 추출법 없이 매우 많은 계산횟수의 crude MCS를 통해 구한 고장확률과 비교하였다. Crude MCS의 경우 100,000개의 샘플을 이용했는데 이때 구해지는 고장확률의 COV는 근사적으로 식 (9)와 같이 주어진다.
Table 3
Failure Probability
Table 3은 고장확률 계산결과를 나타내며, 괄호 안의 %값은 MCS 대비 상대 정확도를 나타낸다. 여기서 식 (9)를 통해 계산된 MCS의 COV값이 각각 약 0.215, 0.033, 0.25가 계산되어 응력파괴 및 볼트파괴에 대한 MCS 신뢰도는 낮았다. 즉 이 두 고장모드에 대해서는 MCS 결과가 큰 변동성이 있다는 것이다. 이를 향상시키기 위해서는 보다 많은 샘플링이 필요하나, 본 연구에서는 그 대신 균열파괴 하나에 대해서만 MCS 결과와 비교하였다. 그 결과 본 방법의 정확도가 MCS 대비 상대 정확도가 103%로 나타난 반면, FORM의 경우 130%로 나타났다. 다른 고장모드에 대해서도 MCS를 정확하다고 간주하면 유사한 정확도 차이를 보이고 있어 본 방법이 FORM 대비 상대적으로 더 우수함을 나타내고 있다. 한편 계산시간을 보면 MCS 소요시간이 약 3주인데 반해 FORM은 11분, 본 방법은 27분의 시간이 소요되어 시간적으로 매우 유리함을 알 수 있다.
3. 결 론
본 연구에서는 반응표면기법과 중요도 추출법을 통해 고장확률을 계산하는 로켓 모터의 내압으로 인한 신뢰도예측 절차를 제안하였다. 내압 및 응력해석 과정에 있어서 LHS와 반응표면기법을 이용하여 MEOP의 확률분포를 계산하였으며, 2D 축대칭모델의 유한요소해석을 통해 MEOP으로 인한 로켓 모터 내부의 von-Mises 응력과 주응력을 계산하였다. 또한, 계산된 결과와 앞서 정의된 고장함수를 이용하여 FORM을 수행하여 MPP를 탐색하고, 이 점을 중심으로 반응표면기법과 중요도 추출법을 적용하여 각각의 고장모드의 고장확률을 계산하였다.
본 논문에서 제시된 기법이 현재 범용적으로 사용되는 FORM 방법과 비교하여 고체로켓 모터의 고장확률을 정확하고 효율적으로 예측할 수 있음을 알 수 있었다. 또한, 본 연구에서 제시된 기법은 근사식으로도 비교적 정확한 결과를 얻을 수 있으며, 계산 과정에서 시간적 · 비용적으로 매우 큰 이점을 갖는 것을 알 수 있었다.
