Username(ID) Password Login

Forgot
my username Forgot
my password Register

Sorry.

You are not permitted to access the full text of articles.

If you have any questions about permissions,

please contact the Society.

죄송합니다.

회원님은 논문 이용 권한이 없습니다.

권한 관련 문의는 학회로 부탁 드립니다.

발간 현황

Journal of the Korean Society of Propulsion Engineers - Vol. 25 , No. 1

[Paper List] [Go to Volume List]


[ Research Paper ]
Journal of the Korean Society of Propulsion Engineers - Vol. 25, No. 1, pp. 1-11
Abbreviation: KSPE
ISSN: 1226-6027 (Print) 2288-4548 (Online)
Print publication date 28 Feb 2021
Received 16 Oct 2020 Revised 08 Dec 2020 Accepted 13 Dec 2020
DOI: https://doi.org/10.6108/KSPE.2021.25.1.001
다목적함수 최적화 기법을 이용한 우주발사체의 포고억제기 설계
윤남경^a ; 유정욱^b ; 박국진^b ; 신상준^a^{, *}


Pogo Suppressor Design of a Space Launch Vehicle using Multiple-Objective Optimization Approach
NamKyung Yoon^a ; JeongUk Yoo^b ; KookJin Park^b ; SangJoon Shin^a^{, *}
aInstitute of Advanced Aerospace Technology, Seoul National University, Korea
bDepartment of Mechanical and Aerospace Engineering, Seoul National University, Korea
Correspondence to : ^* E-mail: ssjoon@snu.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.


Funding Information ▼ Space Propulsion Research Center Ministry of Science, ICT and Future Planning National Research Foundation of Korea NRF-2013R1A5A1073861

초록

포고 현상은 액체추진 로켓에서 발생하는 축방향의 동적 불안정 진동이다. 동체의 고유진동수와 추진제 공급계의 주파수가 가까와 지면 전체 시스템이 불안정 현상을 보인다. 포고 현상을 예측하기 위해 1단의 추진제 (산화제 및 연료) 탱크는 쉘 요소로, 나머지 구성 요소인 엔진 및 상단은 mass-spring으로 모델링하여 구조해석을 수행하였다. 추진제 공급계의 압력 및 유량 섭동예측에는 transmission line model이 사용되었다. 본 논문에서는 이와 같이 수행된 구조 및 유체 모델링을 통합하여 폐루프 전달함수를 구성하였다. 포고 억제기는 수동적인 방법으로 압력 섭동을 흡수하는 분기관 및 accumulator로 구성되며 추진제 공급계 중간에 위치한다. 발사체의 비행과정 동안 포고현상을 억제하는 설계 최적화를 위한 설계변수로는 분기관 및 accumulator의 직경 및 길이로 설정하였다. 목적함수로는 포고 억제기의 질량, 그리고 추진제 질량에 따른 폐루프 전달함수의 에너지 최소화로 설정하여 다목적함수 최적화를 수행하였다.

Abstract

POGO is a dynamic axial instability phenomenon that occurs in liquid-propelled rockets. As the natural frequencies of the fuselage and those of the propellant supply system become closer, the entire system will become unstable. To predict POGO, the propellant (oxidant and fuel) tank in the first stage is modeled as a shell element, and the remaining components, the engine and the upper part, are modeled as mass-spring, and structural analysis is performed. The transmission line model is used to predict the pressure and flow perturbation of the propellant supply system. In this paper, the closed-loop transfer function is constructed by integrating the fuselage structure and fluid modeling as described above. The pogo suppressor consists of a branch pipe and an accumulator that absorbs pressure fluctuations in a passive manner and is located in the middle of the propellant supply system. The design parameters for its design optimization to suppress the decay phenomenon are set as the diameter, length of the branch pipe, and accumulator. Multiple-objective function optimization is performed by setting the energy minimization of the closed loop transfer function in terms of to the mass of the pogo suppressor and that of the propellant as the objective function.


Keywords: POGO, Liquid Rocket Engine, Propellant Feedline, Longitudinal Stability, Multiple-Objective Optimization 키워드: 포고, 액체로켓엔진, 추진제 공급관, 축방향 안정성, 다목적함수 최적화

1. 서 론

우주발사체에서는 다분야 간의 연성에 의해 불안정 현상이 발생할 수 있다. 그 중에서 포고 현상은 액체 추진엔진을 사용하는 발사체에서 구조물, 추진제 공급시스템 및 추력시스템의 상호작용에 의해 발생하는 축 방향 불안정성 현상이다. 구조물의 진동은 추진제 공급라인에 영향을 주고, 추진제의 유량 및 압력 섭동은 추력에 섭동을 일으킨다. 추력 섭동은 다시 구조물의 진동을 발생(self-excited vibration)시키는 폐루프가 형성되어 포고 현상을 발생시킨다. 이 현상에 대한 최초 발견은 1960년대에 이루어진 유인 우주 계획에서 보고되었다[1,2]. Titan II Gemini 프로그램에서 발사 90초후에 약 30초 동안 10-13Hz로 축방향 진동이 발생하였다. 11Hz에서 최대 2.5g 가속도로 진동이 발생하였으며, NASA에서 0.25g 이하로 낮추는 것으로 요구되었다[2]. Saturn V 발사체에서 1단 추력 비행구간 105-140초 사이에 5Hz 대역의 최대 가속도 0.6g의 축방향 진동이 포고불안정 현상이 발생되었다. 우주 왕복선에는 이러한 불안정 현상을 막기 위해 주 엔진(Space Shuttle Main Engine)에 수동 억제기를 장착하였다. 포고 현상은 비행 초기 단계에서도 발생할 수 있으며 사고 발생시 인명 및 재산 피해를 초래할 수 있다. 따라서 포고 현상의 정확한 예측이 요구되며 발사체의 축 방향 진동에 대한 해석과 추진제 공급계 및 추진기관 시스템 전체에 대한 정밀한 모델링이 요구된다. 포고 현상을 예측하기 위해서는 Fig. 1과 같이 동체의 구조 모드, 연료공급계 및 추력시스템으로 구성된 폐루프 해석이 요구된다.

Fig. 1
POGO analysis procedure of a launch vehicle.

본 논문에서는 분기관 형식(branch type)의 수동 포고억제기를 설계하고자 한다. 우주발사체는 연료가 전체 중량의 80∼90%를 차지하므로[3] 비행중에 추진제가 소모되면서 고유진동수가 변하게 된다. 구조모드 해석 시 잔류 추진제량을 고려하여 구조 전달함수를 해석하고 여기에 추진제 공급계의 전달함수와 함께 폐루프를 구성하였다. 기체 구조와 추진공급계, 두 시스템 간의 상호작용에 의한 에너지 전달을 최소화하도록 목적함수를 설정하였다. 전 비행구간동안 포고현상을 억제하는 포고억제기 설계를 위해 질량특성을 100%에서 0%까지 6단계로 구분하여 구조모드 해석을 수행하였다. 비행동안 산화제 공급계의 전달함수는 일정하다고 가정하여 한 개의 전달함수로 구성하였다. 구조전달 함수 6개와 공급계 전달함수 1개로 총 6개의 폐루프 전달함수를 구성하였으며 폐루프 전달함수는 최적화 수행 시 목적함수 계산에 사용된다. 여기에 분기관과 억제기의 부피를 목적함수로 추가하여 총 7개의 목적함수로 구성된 다목적 함수 최적화 설계로 구성하였다. 포고억제기의 형상에 의해 공급계의 고유 진동수가 낮아져서 폐루프 시스템의 진폭이 감소 한다. 따라서 설계변수로 분기관과 억제기의 직경과 길이로 설정하였다. 최적화를 통해 도출된 분기관의 직경과 길에 해당하는 포고억제기의 형상을 반영한 폐루프에 sine sweep 함수로 가진하여 포고 억제기의 효과를 살펴보고자 한다.

2. 본 론

2.1 개요

포고해석을 위해 필요한 다분야 해석, 즉 동체 구조 해석과 추진제 공급 해석을 진행하였다. 산화제 및 연료가 소모되면서 발사체의 중량이 감소하고 이에 따라 비행시간별로 구조 모드 고유진동수가 증가하게 된다. 기체의 축 방향 구조 모드의 고유진동수와 공급계-추력시스템의 압력 주파수가 가까와지면 시스템이 불안정해진다. 포고억제기를 설치 함으로써 공급계-추력시스템의 고유진동수와 구조 모드의 고유진동수를 분리하여 포고 불안정 현상을 방지 할 수 있다[4]. 따라서 이 섹션에서는 분기관 형식의 포고억제기 설계를 위한 다목적함수 최적화 문제를 정의하고자 한다. 다음으로 폐루프 전달함수를 구성하기위한 구조 모델링, 산화제 공급계 모델링 및 추력시스템 모델링을 살펴본다. 폐루프 전달함수 구성 및 최적화 과정이 Fig. 2에 묘사되어 있다.

Fig. 2
Flowchart of the optimization and construction of the closed-loop transfer function.

Fig. 3
1D-3D integrated structural modeling [10].

2.2 최적화 문제 정의

포고 억제기의 설계는 메타 모델링과 다목적 최적화 과정으로 이루어진다. 목적함수는 이전 절에서 소개한 연료소모를 고려한 6개의 폐루프 전달함수의 파워 스펙트럼 밀도(power spectral density, PSD)과 분기관을 포함함 포고억제기의 부피로 구성된다. 이러한 목적함수를 3개 이상 갖는 다목적 함수는 일반적인 단일 또는 두개의 목적함수에 비하여 최적화 결과 선택의 어려움이 존재한다. 단일 목적함수의 경우 목적함수가 최대 또는 최소가 되는 설계변수를 선택하거나 두개의 목적함수를 갖는 경우 Pareto front에 기반하여 최적해 선택 가이드라인을 제시할 수 있으나, 목적함수가 3개 이상이면 선택조건을 제시하기 어려운 점이 있다. 따라서 도메인 지식을 활용하여 선택하여야 한다. 본 논문에서는 폐루프 시스템의 파워 스텍트럼 밀도 6개와 포고억제기 부피로 이루어진 목적함수 중에서 가장 큰 목적함수의 값을 최소화 하는 설계변수를 최적해로 선택 하였다. Eq. 1에 정의된 목적함수는 연료 잔량 상태 6개에 대한 구조 전달함수와 공급계의 전달함수가 결합된 폐루프 전달함수의 파워 스펙트럼 밀도와 포고억제기 부피의 최소화를 의미한다. 포고 억제를 위한 분기관 및 억제기는 일정한 원형 단면적을 가지는 실린더로 가정하여 설계변수를 각각의 직경 및 길이로 설정하였다. 설계변수는 Eq. 2-4에 제시된 상한과 하한 사이에서 Latin hypercube 샘플링을 통하여 각각 100,000개로 구성된다. 샘플링된 설계변수 값을 사용하여 폐루프 전달함수의 고유 응답을 6×10⁵개 구성하였다. 샘플링된 설계변수를 입력으로 설정하고 폐루프 전달함수의 응답을 출력으로 설정하여 2차 함수로 근사된 반응 표면을 생성 하였다. 반응표면 구성은 MATLAB의 내장함수인 rstool을 사용하였다[5]. 다목적최적화는 유전자 알고리즘에 기반을 둔 kRVEA 알고리즘을 사용하여 수행되었다[6].

f=mini,j⁡PSDGisHjs,V,i=1:6, j=1:105

(1)

0.010≤xi,x2≤0.015 m

(2)

0.005≤x3≤0.010 m

(3)

0.010≤x4≤0.020 m

(4)

여기서 G_i(S)는 잔류 추진제량 상태 6개([100, 80, 60, 40, 20, 0]%)에 해당하는 구조 전달함수, H_i(S)는 설계변수(x₁, x₂, x₃, x₄)_j 샘플링에 대한 산화제 공급계의 전달함수, V는 분기관과 포고억제기의 부피, x₁는 accumulator 길이, x₂는 분기관 길이, x₃는 분기관 반지름, 그리고 x₄는 accumulator의 반지름을 의미한다.

Fig. 4에 잔류 추진제량의 변화에 따른 구조 모드의 전달함수 크기변화와 Table 1에 고유진동수의 변화가 제시되어 있다. 분기관 및 포고억제의 부피는 샘플링된 설계변수로부터 계산된다 (V=πx42x1+πx32x2). 섹션 2.3-2.4에서 Eq. 1에서 사용되는 폐루프 전달함수를 구성하는데 필요한 구조 및 산화제 공급계의 전달함수 생성과정을 살펴본다.

Fig. 4
Transfer fuction of the fuselage structure in terms of fuel consumption [9].

Table 1.
Changes of the fuselage natural frequency in terms of the fuel consumption.

Fuel mass	100% j=1	80% j=2	60% j=3	40% j=4	20% j=5	0% j=6
ω_f,j [Hz]	9.3	10.1	10.5	11.3	13.3	16.2

2.3 구조 모드 해석

발사체의 축 방향 구조모드 해석 즉, 전달함수 G(s) 해석을 위해서는 추진체 탱크와 공급관의 유탄성(hydroelasticity) 효과를 반영하는 것이 필수적이다. 하지만, 3차원 동체 구조 모델의 유한요소 해석은 많은 자유도 수와 계산시간을 요구하기 때문에, 축 방향 모드의 고유진동수를 직관적으로 예측하기 쉽지 않다. 그러므로 추진제를 포함한 탱크를 3차원 유탄성 해석으로 수행하고, 유상하중, bulkhead, 어댑터 등 구성품들은 1차원 질량-스프링 요소로 구성하였다. 특히, 구조 이득(structural gain), G_e을 구하기 위해, 엔진 시스템이 1차원 질량-스프링 요소로 표현되어 하나의 모드 변위를 가지는 축 방향 움직임을 한다고 가정한다. 이러한 1차원 및 3차원 요소의 결합 해석 기법은 발사체의 축 방향 진동 모드를 정확히 예측할 수 있도록 하고, 포고 해석의 메커니즘을 단순화시킬 수 있다.

동체 구조 해석에서는 동체가 축 방향으로 받는 추력에 대한 응답을 계산하며, Eq. 5와 같이 구성된다.

Gs=Ges2s2+2ζωns+ωn2

(5)

이때 구조 이득(gain)을 결정짓는 요소는 Ge=φeφTM이며, 이는 진동 해석 결과로부터 얻어지는 축 방향 고유벡터 각 점에서의 모드 변위(modal displacement) 및 일반화된 질량(generalized mass)으로부터 얻어진다.

모드 변위는 각 구성품이 연결되어 있는 점에서의 변위이며, 포고 시스템에서는 각 엔진점에서의 변위, 산화제 공급 관과 탱크, 혹은 펌프와의 접점 부분 변위 등에 해당된다. 이에 해당하는 각 점에서의 변위와 n차 모드에 대한 지배방정식은 Eq. 6-7과 같다.

r¯it=∑n=1NMr¯nit=∑n=1NMϕ¯niqnt

(6)

Mnq¨nt+2ζnωnq˙nt+ωn2qnt=Qnt

(7)

축방향 진동 특성 G(s)를 구하기 위해서, 추진제를 포함한 탱크를 3차원으로, 나머지 구성품들을 1차원 요소로 모델링하였다. 해석의 주요 관심사는 엔진 모드 변위 φ_e와 발사체의 축방향 모드 즉, ω_n, ζ_n, M_n이다

액체 추진제를 사용하는 우주발사체에서 탱크와 추진제는 발사체 중량의 대부분을 차지한다. 또한 탱크의 경우 그 크기에 비해 두께가 매우 얇아 유연한 구조를 가진다. 따라서 탱크가 액체 추진제에 미치는 연성 효과가 크고, 추진제가 탱크 구조의 동적 응답 특성에 중요한 영향을 끼치기 때문에 구조, 유체 간 상호작용인 유탄성 효과를 고려한 해석이 필요하다.

Fig. 3에 본 논문의 해석 대상물인 Atlas/Centaur/Surveyor 발사체의 1-3차원 모델링을 표시하였다. Atlas 1단에서 포고의 발생가능성이 높으므로 추진제의 영향을 고려한 유탄성 해석을 위하여 쉘로 모델링하고 나머지 요소들을 1차원 질량-스프링 요소로 모델링하였다. 엔진을 질량으로 표현했기 때문에 모드 변위, φ_e를 구할 수 있다. 추진제의 영향을 고려한 유탄성 해석을 위해 가상 질량 기법이 적용된 NASTRAN의 MFLUID 요소를 사용하여 부가 질량(added mass) 행렬을 구하였다. NASTRAN에서는 간단한 해를 제공하는 소스(source)를 외부 경계(outer boundary)에 분포 시켜 라플라스 방정식을 도출하고 이방정식을 Helmholtz 방법으로 푼다. 알고 있는 경계 움직임과 소스의 움직임을 일치시킴으로써 선형 행렬 연산을 통해 소스의 크기를 구할 수 있다. 소스의 크기로부터 압력이 결정되며 이를 통해 격자(grid) 점에 작용하는 하중을 구할 수 있다[11-13]. 만약 유체의 포인트 소스 σ_j(단위면적당 체적 유량율)가 r_j에 위치하고 면적 A_j에 작용한다고 가정하면, 위치 r_i에서 속도벡터 u˙i는 다음과 같다.

u˙i=∑j∫Ajσjeijri-rj2dAj

(8)

여기서 e_ij는 j에서 i로의 단위 벡터이다. Eq. 8에서 u˙i 의 gradient는 라플라스 방정식을 만족하는 포텐셜 함수가 됨을 알 수 있다. 위치 r_i에서 압력 p_i는 아래와 같다.

pi=∑j∫Ajρσ˙jeijri-rjdAj

(9)

Eq. 8-9를 유한 요소의 표면에 대한 적분을 수행하여 행렬로 표현하면 다음과 같다. Eq. 9의 압력을 힘으로 변환하기 위해 NASTRAN에서 면적분(area integration)이 수행된다.

u˙i=χσ

(10)

F=Λσ˙

(11)

Eq. 10을 Eq. 11에 대입하여 정리하면 다음과 같다.

F=Λχ-1u¨=Mfu¨

(12)

결과적으로, 유체의 포인트 소스 σ_j에 대해 압력 및 속도 벡터를 정식화 하면 탱크 내 추진제를 단순한 질량 행렬, [M^f]로 표현할 수 있다. 따라서 연료 탱크 구조의 질량 행렬에 추진제의 질량이 더해져 축 방향 진동모드를 구한다.

2.4 공급관 및 추력해석 모델링

발사체의 경우 공급관을 흐르는 유체의 압력 및 유량 섭동 해석을 위해 Michalopoulos[7], Rubin[4], Oppenheim[8] 등에 의해 유도된 transmission line modeling이 사용되었다. 압력이 있는 유체가 흐르는 관의 유동 해석을 위해 Fig. 5와 같은 형상을 가진 직선관에 대하여 N개의 요소로 분할하여 포고억제기를 포함한 공급관의 압력 및 유량 섭동을 계산하게 된다. 압력과 유량을 평균 유동(mean flow)에 대한 미소-섭동을 가정하여 선형 방정식을 유도 한다. 원형 단면을 가진 공급관을 흐르는 유체의 점성과 압축성을 포함한 1차원 모멘텀 방정식(convective term은 무시)과 질량보존법칙(단열흐름을 가정) 을 나타내면 Eq. 13-14와 같다[7].

IiQ˙i=Pi-Pi+1-Rf,iQi0Qi

(13)

Cf,iP˙i=Qi-1-Qi

(14)

Fig. 5
Straight feeline along with a POGO suppressor (branch + accumulator) [9].

여기서 하첨자 i는 요소번호(i=1:N), Q_i_-1는 i요소로 들어가는(양의부호) 유량, Q_i는 i요소에서 유출되는(음의 부호) 유량, f는 마찰계수, P_i는 i번째 요소의 압력을 의미한다. Eq. 13-14에서 정의된 inertance, resistance, compliance는 Eq. 15와 같다.

Ii=ρLOXLiAi,Cf,i=1Ks,LOX=1ρLOXcLOX2, Rf,i=fρLOX2Ai2Li+LeD

(15)

여기서 L_i는 요소길이, A_i는 요소 단면적, D는 관(pipe)의 내경, c_LOX는 음향속도(acoustic speed), L_e는 손실에 대한 등가길이를 의미한다. Eq. 13-14에서 감쇠를 무시하면, 즉 R_f = 0 으로 가정하면 Eq. 16과 같이 정리할 수 있다.

AQ˙+BP=F1CP˙-BTQ=F2

(16)

여기서 Q = [Q₁,Q₂,....,Q_N]^T, P = [P₁,P₂,....,P_N]^T, F₁ = [0,0,....,-P_N+1]^T, F₂ = [Q₀,0,....,0]^T, 그리고 행렬 B는 단일 직선관의 경우 Eq. 17과 같고 Fig. 5와 같이 분기관을 갖는 경우 N×(N+1) 크기가 되며 Eq. 18과 같이 표현된다.

Bi,i=-1,Bi,i+1=1,Bij=0otherwise

(17)

Bi,i=-1; i=1,N; i≠kBk,j=-1;Bi,i+1=-1; i=1,N

(18)

Eq. 19-20과 같이 미소 섭동을 가정하면 비선형인 Eq. 13-14을 선형화할 수 있다. 이를 행렬로 나타내면 Eq. 21-22와 같다.

δQLOX=QLOX,0+δqLOX,δqLOX≪QLOX,0

(19)

δPLOX=PLOX,0+δpLOX,δpLOX≪PLOX,0

(20)

Aq¨+Dq˙+EQq=g1

(21)

Cp¨+Hp˙+Epp=g2

(22)

여기서 행렬 [A], [D] 그리고 [C]는 Eq. 15에 정의된 inertance, resistance 및 compliance를 대각요소로 갖는 행렬로서 Eq. 23-25와 같고, 행렬 H,E_Q,E_p,g₁,g₂가 Eq. 26-28에 정의되어 있다.

C=C10000C20000⋱0000CN=CiI

(23)

D=2Rf1Q1000002Rf2Q200000⋱00002RfNQN0=2RfiQi0I

(24)

A=I10000I20000⋱0000IN=IiI

(25)

H=BTA-1DB-TC

(26)

EQ=BC-1BT, EP=BT A-1B

(27)

g1=f˙1-BC-1f2, g2=f˙2+BTA-1f1+DB-Tf2

(28)

Eq. 18-19에서 입구와 출구에서의 압력과 유량 섭동을 일반화된 좌표계와 고유값, qnpωnp,∑i=1Npxi,np 및 qnqωnq,∑i=1Nqxi,nq로 표시 할 수 있고 이를 이용하여 전달함수로 나타내면 아래와 같다.

AnQnQs2+DnQnQs+ωnQ2qnQ=G1s

(29)

CnPnPs2+HnPnPs+ωnP2qnP=G2s

(30)

여기서 G₁(s) = L[g₁(t)], G₂(s) = L[g₂(t)], L[•]는 Laplace 변환을 의미한다. Eq. 29-30를 비감쇠 시스템으로 가정하여 추진제 공급계의 고유진동수와 모드를 고유치 해석으로부터 구한다[9]. 섹션 2.3-2.4를 통해 Fig. 2에 제시된 폐루프 전달함수 생성을 위한 구조 해석 및 공급관 해석을 살펴보았다. Fig. 2에 제시된 캐비테이션 및 추력의 전달함수는 [10]의 결과가 사용되었다.

2.5 수치해석 결과

Eq. 1에서 정의된 최적화 문제의 해를 구할 때, 매 최적화마다 폐루프 전달함수를 모든 시간에 대하여 계산하는 방식에는 많은 계산 시간이 요구된다. 따라서 파워스펙트럼밀도 함수를 계산하여 설계변수의 관계를 2차 함수로 반응표면을 구하여 최적화를 수행하였다. 파워스펙트럼밀도 함수는 전달 함수에 비해 공진을 일으키는 폐루프 전달함수의 주파수를 파악하기 용이하며 반응표면을 구성하기 용이하다. 폐루프 전달함수로부터 계산된 파워스펙트럼밀도 계산 결과가 Fig. 6과 같으며, 연료 0%에 해당하는 파워스펙트럼밀도와 설계변수의 관계의 반응표면은 Fig. 7과 같다. kRVEA알고리즘[6]을 적용하여 섹션 2.2에서 정의한 7개의 목적함수를 갖는 다목적함수 최적화를 수행하였으며 결과가 Fig. 8에 제시되어 있다. Fig. 8에서 볼 수 있듯이 목적함수 별로 최적의 해가 다르기 때문에, 예를 들어 2번째 목적함수에서 최소값을 갖는 설계변수가 목적함수 3번에서 최소값을 가지지 않기 때문에 최적해를 선택하는 데 어려움이 있다. POGO 현상이 주로 1단 연소시 발생함을 고려하여 2번째 목적함수와 3번째 목적함수, 즉 잔류 추진제량 80%일 때와 60% 일 때 목적함수가 최소가 되는 설계변수를 비교하여 선택할 수 있다. Eq. 2-4에 제시된 설계변수의 상한과 하한 구속조건이 적용된 최적 설계변수 중 하나가 Table 2에 제시되어 있으며 분기관과 accumulator의 길이는 소수점 셋째자리 버림 된 값이 계산에 사용되었다.

Fig. 6
Power spectral density of the closed-loop transfer function (G(s)H(s)_ij, i = 1, j =1:10).

Fig. 7
Response surface between the design variables and closed-loop transfer fuction (fuel=100%).

Fig. 8
Resuls of the multi-objective optimization.

Table 2.
Optimized POGO supperssor configuration.

Design variables	x1	x2	x3	x4
Optimized results [m]	0.1406	0.1139	0.075	0.1973

추진제 공급관은 내경 0.15m, 길이 10m인 직선관의 형상에 대하여 요소개수 1000개를 사용하여 Eq. 29-30을 해석하여 Table 3에 1차부터 4차 모드까지 제시하였다. Table 3에서 볼 수 있듯이 포고 억제기를 설치하기 이전에는 Table 1에 제시된 구조 모드의 저차 고유 진동수와 추진제 공급시스템의 고유 진동수가 가깝지만 최적화된 포고 억제기를 설치함으로써 두 시스템 간 고유 진동수가 분리됨을 볼 수 있다. Table 2의 형상에 대해 잔류 추진제량 0%에 해당하는 폐루프 전달함수에 0∼60Hz의 sine 함수를 가진 후 이의 응답을 Fig. 9에 제시 하였다. Fig. 9의 왼쪽은 포고억제기를 설치하기 이전의 폐루프 응답이며 오른쪽은 포고억제기 형상을 최적화하였을 때의 폐루프 응답이다. Fig. 9에서 볼 수 있듯이 포고억제기의 형상 최적화를 통한 추진제 공급 시스템의 고유진동수 변화가 폐루프 응답의 크기를 감소시는 데에 효과적임을 알 수 있다.

Table 3.
Changes of the feedline natural frequency by the accumulator.

Mode number	j=1	j=2	j=3	j=4
ω_f,j [Hz] without accumulator	12.51	37.52	62.57	87.60
ω_f,j [Hz] with accumulator	23.71	25.04	71.19	75.12

Fig. 9
Time response of the closed-loop plant (fuel=0%, 0∼60 Hz sine sweep).

3. 결 론

본 논문에서는 우주발체의 포고 억제기 설계를 위한 다목적 최적화 과정을 제시하였다. 또한 목적함수에 사용되는 폐루프 전달함수 구성을 위해 구조 및 추진제 공급시스템 해석을 수행하는 다분야 해석을 수행하였다. 다목적 함수 설계는 목적함수가 3개 이상이므로 이러한 목적에 맞게 개발된 알고리즘의 사용이 요구된다. 또한 단일 또는 2중 목적함수와 다르게 다중의 목적함수에 대한 여러 개의 최적 값 중 어느 것이 좋고 나쁨을 평가에 대한 척도를 제시하기 어렵기 때문에 도메인 지식을 활용하여 설계변수를 선택해야 한다. 해석결과를 요약하면 다음과 같다.

1. NASTRAN의 MFLUID요소[11]를 사용하여 연료 소모에 따라 동체 모드의 전달함수가 변함을 Fig. 4에 구조 모드 전달함수의 크기로 제시하였다.
2. Table 1과 Table 3으로부터 최적화된 형상의 포고억제기를 설치함으로써 1차 구조 모드와 추력-공급계 고유 진동수가 분리됨을 확인하였다.
3. 구조 모드와 추진제 공급계의 고유진동수가 분리되어 Fig. 9처럼 포고현상 발생 시 진폭이 감소함을 확인하였다.
4. Fig. 8을 통해 다목적 최적화 결과 2번째 목적함수에서 에너지 전달이 가장 크기 때문에 해당 목적함수를 최소화하는 설계변수를 선택하는 것이 1단에서 발생하는 포고현상을 억제하는 데에 효과적임을 예측할 수 있다.

본 논문에서 사용된 구조 형상은 Atlas 발사체로 질량-스프링과 shell로 1단의 탱크를 모델링하는 방법을 참조하였다. 그러나 이러한 방법은 격자로 구성된 벽면, 어댑터 링 등의 질량-스프링 요소의 조합하는 과정에 대한 시행착오가 요구된다. 따라서 3차원 구조물을 1차원과 3차원 요소로 축소 모델링하는 방법에 대한 연구가 이후에 필요하다. 그리고 폐루프 전달함수를 구성하는 과정에서 연료의 양을 기준으로 구조 전달 함수를 모델링하였으나 비행 이벤트(비행 시간)가 고려되지 않았다. 따라서 궤도 해석과 통합하여 관심 이벤트 구간에 대한(예를 들어 킥턴, 최대동압 구간) 연료량을 고려하여 모델링함으로써 관심있는 비행시간에서 포고해석이 가능하도록 확장이 요구된다. 추가적으로 현재 재사용 발사체에 대한 관심이 증가로 추진제 공급라인에 밸브를 추가하여 재연소할 수 있는 추진제 공급계에 대한 연구를 필요로 한다. 따라서 밸브의 조작에 의한 water-hammer 현상 등을 검토할 수 있도록, 추진제 공급계 해석 시 과도응답의 해석 능력이 또한 요구된다.

Acknowledgments

본 연구는 서울대학교 차세대 우주추진 연구센터와 연계된 미래창조과학부의 재원으로 한국연구재단의 지원을 받아 수행한 선도연구센터지원사업(NRF-2013R1A5A1073861)의 연구 결과입니다.

References


1.	Larsen, C.E., NASA Experience with Pogo in Human Spaceflight Vehicles, NATO RTO Symposium ATV-152 on Limit-Cycle Oscillations and Other Amplitude-Limited, Self-Excited Vibrations, May 5, 2008.
2.	Fenwick, J., “POGO,” Threshold: Pratt & Whitney Rocketdyne’s Engineering, Journal of Power Technology, 1992.
3.	Hopkins Jr., J., Hopkins, J. and Isakowitz, S., International Reference Guide to Space Launch Systems, Fourth Edition, American Institute of Aeronautics and Astronautics, Inc. 2004.
4.	Rubin, S., “Longitudinal Instability of Liquid Rockets Due to Propulsion Feedback (POGO),” J. Spacecraft and Rockets, Vol. 3, No. 8, pp. 1188-1195, 1966.
5.	Anonymous, “rstool,” https://kr.mathworks.com/help/stats/rstool.html.
6.	Tian, Y., Cheng R., Zhang X. and Jin, Y., “PlatEMO: A MATLAB Platform for Evolutionary Multi-Objective Optimization [Educational Forum],” IEEE Computational Intelligence Magazine, Vol. 12, No. 4, pp. 73-87, 2017.
7.	Michalopoulos, C.D., Clark, R.W. and Doiron, H.H., “Fourth Annual Thermal and Fluids Analysis Workshop : Acoustic Modes in Fluid Networks,“ NASA Conference, Cleveland, Ohio, US., 1992.
8.	Oppenheim, B.W. and Rubin, S., “Advanced Pogo Stability Analysis for Liquid Rockets,” Journal of Spacecraft and Rockets, Vol. 30, No. 3, pp. 360-373, 1993.
9.	Park, K.J., Yoo, J.U., Lee, S.H., Nam, J.H., Kim, H.J., Lee, J.Y., Roh, T.S., Yoh, J.J., Kim, C.A. and Shin, S.J., “Pogo Accumulator Optimization Based on Multiphysics of Liquid Rockets and Neural Networks,“ Journal of Spacecraft and Rockets 2020, pp. 809-822.
*10.*	Park, K.J., Lee, S.H., Lee, S.G. and Shin, S.J., “Longitudinal Characteristics Analysis of a Space Launch Vehicle using One- and Three-Dimensional Combined Modeling for Pogo Prediction,“ 17-19 September 2018, Orlando, FL, 2018 AIAA SPACE and Astronautics Forum and Exposition.
*11.*	Anonymous, “Advanced Dynamic Analysis User’s Guide,” DOC9180, MSC.Software Corporation 2014.
*12.*	Lee, D,Y. and Choi, H.S., “A Study on the sloshing of Cargo Tanks Including Hydroelastic Effects,” Journal of the Society of Naval Architecture of Korea, Vol. 35, No. 4, pp. 27-31, 1998.
*13.*	Byeon, J.H., Cho, H.J., Baek, S.J., Prabowo, A.R., Bae, D.M. and Sohn, J.M., “Numerical approaches in idealizing mass of fluid in tank for ship vibration analysis,“ MATEC Web of Conferences, Vol. 159, 2018.