Journal of the Korean Society for Aviation and Aeronautics

The Korean Society for Aviation and Aeronautics

J. Korean Soc. Aviat. Aeronaut. 2019; 27(4):9-20

pISSN: 1225-9705, eISSN: 2466-1791

DOI: https://doi.org/10.12985/ksaa.2019.27.4.009

Original Article

이동변형격자 기법을 활용한 외부장착물 분리운동 해석

안병희^*, 김동현^**

External Store Separation Analysis Using Moving and Deforming Mesh Method

Byeong Hui Ahn^*, Dong-Hyun Kim^**

^*경상대학교 기계항공정보융합공학부 대학원생

^**경상대학교 기계항공정보융합공학부 및 항공기부품 기술연구소 교수

연락저자 E-mail : dhk@gnu.ac.kr 연락저자 주소 : 경남 진주시 가좌동 900번지 공과 대학 407동 710호

© Copyright 2019 The Korean Society for Aviation and Aeronautics. This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received: Oct 20, 2019; Revised: Dec 26, 2019; Accepted: Dec 30, 2019

Published Online: Dec 31, 2019

ABSTRACT

A military aircraft generally includes external stores such as fuel tanks or external arming, depending on the purpose of the operation. When a store is dropped from a military aircraft at high subsonic, transonic, or supersonic speeds, the aerodynamic forces and moments acting on the store can be sufficient to send the store back into contact with the aircraft. This can cause damage to the aircraft and endanger the life of the crew. In this study, time accurate computational fluid dynamics (CFD) with dynamic moving grid (moving and deformable mesh, MDM) technique has been used to accurately calculate store trajectories. For the verification of the present numerical approach, a wind tunnel test model for the wing-pylon-finned store configuration has been considered and analyzed. The comparison results for the ejected store trajectories between the present numerical analysis and the wind tunnel test data at the Mach number of 0.95 and 1.2 are presented. It is also importantly shown that the numerical parameter of MDM technique gives significant effect for the calculated store trajectory in the low-supersonic flow such as Mach 1.2.

Keywords: External Store Separation(외부장착물 분리); CFD(전산유체역학); Moving Grid(이동격자); Deforming Grid(변형격자); Unstructured Mesh(비정렬격자)

Ⅰ. 서 론

실제 운용되고 있는 다양한 종류의 군용 항공기들은 연료탱크 또는 무장 등과 같은 외부장착물이 장착된 상태로 운용되고 있다. 외부장착물(external store)은 종류가 많지만, 유도제어 개념이 아닌 경우는 공중에서 분리 투하되어 외부유동 하중 및 동역학적 평형 관계에 의한 운동이 결정되게 된다. 항공기에서 분리투하된 외부장착물은 분리 후에 항공기와 충돌 위험성이 없어야 하고, 설계된 궤적으로 정확한 운동이 요구되기 때문에 동적 안전성을 사전에 충분히 검토 및 검증할 필요성이 있다. 외부장착물 투하 안전성 검증은 풍동시험을 수행하거나, 실제 비행 투하 비행시험(flight test)을 하는 방법과 전산유체역학(computational fluid dynamics, CFD) 해석 기법을 활용하는 방법이 있다(Fig. 1).

Fig. 1. Relation between CFD, wind tunnel and flight tests

Download Original Figure

외부장착물이 있는 항공기 날개는 복잡한 형상으로 인하여 천음속 및 초음속 영역에서 강한 충격파 간섭현상이 유발되게 된다. 이로 인하여 장착물 투하 시 발생 가능한 동적 불안정성을 실험적 혹은 해석적 방법을 통하여 사전에 예측 및 분석하기 위한 연구들이 수행되어 왔다. Carman[1]에 의해 미국 공군 Arnold Engineering Development Center(AEDC)에서 무장분리에 관한 풍동시험 연구가 수행되었는데, captive trajectory support (CTS)를 활용한 시험기법을 제시하였다. Heim[2]은 AEDC의 풍동시험 설비를 이용하여 천음속 영역에서 3차원 날개에 파일런 및 외부장착물이 부착된 형상에 대해 날개 및 외부장착물 표면의 압력 분포를 계측하였다. 또한 파일런에서 외부장착물이 분리된 후, 시간에 대한 위치와 피치, 롤, 요 각도 등을 계측하고, 상세한 풍동시험 데이터를 제시하였다.

전산유체역학(CFD) 기법을 활용한 장착물 분리운동 해석은 정렬중첩격자(structured overset mesh)와 동적비정렬격자(dynamic unstructured mesh)에 기반을 둔 연구가 주로 수행되어 왔다[3-12]. 정렬중첩격자 방식은 수치 안정성은 우수한 편이나 물체의 이동에 따른 물체주변 이동격자와 배경격자의 중첩영역에서 데이터 교환과정에서 수치보간(numerical interpolation) 오차 발생 가능성이 잠재되어 있다. 반면에 동적비정렬격자 방식은 수치보간 오차 발생은 없으나, 물체의 움직임을 고려한 유동해석 과정에서 이동 및 변형격자 생성에서 격자 꼬임현상으로 인한 수치에러 발생 가능성이 있다.

기존 국내외 연구사례에서 천음속 장착물 분리운동 CFD 해석기법의 정량적인 검증을 위해 참고문헌[2]의 풍동시험 결과가 대표적으로 비교되고 있다. 하지만 동적비정렬격자에 기반한 대부분의 연구는 마하 0.95에서만 CFD 해석 및 풍동시험 결과를 비교 제시하고 있으며, 마하 1.2 조건에 대해서는 풍동시험과 잘 일치하는 비교결과를 찾아보기 어렵다.

본 연구에서는 동적비정렬격자 기반의 이동변형격자(moving and deforming mesh, MDM) CFD기법을 활용하여 마하 0.95 및 1.2의 천음속(transonic) 영역에서 외부장착물 분리운동 해석을 수행하고, 기존 풍동시험과 비교결과를 제시하였다. 특히 MDM 기법 적용 시 필수적으로 설정해야 하는 spring constant factor (SCF)가 외부장착물 분리운동 수치해석 결과에 미치는 영향에 대해 고찰하였다. 본 연구는 실제적으로 외부장착물 분리할 경우, 반발력에 따른 탄성체 날개의 공력탄성학적 변형 및 진동 효과의 정확한 해석을 위한 선행연구로 수행되었다. 외부장착물 공력간섭 효과를 고려한 2-way coupled CFD-CSD 비선형 전산공력탄성학해석[13-15]은 고난이도의 수치해석 기법이 요구되기 때문에 동적비정렬격자 기법이 우선적으로 활용될 필요성이 크다.

Ⅱ. 수치해석 기법

2.1 비정상 유동해석

압축성 유동에 대한 비정상(unsteady) 오일러(Euler) 방정식은 다음과 같이 나타낼 수 있다.

∂ ∂ t ∫ Ω Q d V + ∮ ∂ Ω F (Q) ⋅ n → d S = 0

(1)

여기서, Q와 F는 각각 보존변수 및 비점성 플럭스를 의미하며, 아래와 같이 표현된다.

Q = [ρ, ρ u, ρ v, ρ w, e o] T

(2)

F (Q) ⋅ n → = (ρ u ¯ ρ u u ¯ + p n x ρ v u ¯ + p n y ρ w u ¯ + p n z e o u ¯ + p u n)

(3)

u ¯ = n → ⋅ (u − u g)

(4)

u n = n → ⋅ u

(5)

위 식에서 Q, u_g, e_o, $n →$ 은 각각 보존변수, 격자속도, 단위체적당 전에너지(total energy) 및 물체 표면에 수직한 단위벡터를 의미한다. 또한 이상기체 가정을 적용하면 다음과 같은 압력관계식을 얻을 수 있으며, 공기의 경우 비열비 γ는 1.4이다.

p = (γ − 1) [e − 12 (u 2 + v 2 + w 2)]

(6)

식 (1)~(5)는 전산유체해석을 위해서 공간이산화와 시간이산화 과정을 거쳐 연립대수방정식으로 변환되고, 이 식을 풀어 이산화 지배방정식에 대한 해를 구하게 된다. 본 연구에서는 유한체적법((finite volume method) 기반의 CFD 유동해석 프로그램인 Fluent(Ver.18.1)과 장착물 분리운동 모사를 위한 C언어 기반의 User-Defined Function (UDF) 코드를 별도로 작성하여 수행하였다. 각각의 셀 표면을 통한 플럭스는 Roe의 flux-difference splitting (FDS) 기법을 적용하였으며, 공간차분화 gradient는 least squares cell based 및 2차 정확도 풍상차분법(second-order upwind scheme)을 적용하였다.

2.2 동적격자 생성기법

CFD 유동해석 영역 내에서 물체가 운동하는 문제나 유동과 구조진동 효과를 동시에 고려한 전산공력탄성학 해석[ ]을 수행하는 경우, 동적격자 재생성 과정이 필요하다. 물체의 운동 변위 또는 회전이 큰 경우는 격자의 volume 또는 skewness 조건에 근거하여 낮은 품질의 격자 셀(cell)들은 필요한 경우, 국부적으로 격자재생성(remesh)을 하게 된다. 반면에 물체의 움직임이 작은 경우는 국부적인 smoothing 기법이 사용된다. 이 경우 격자점(node)들은 cell의 품질을 향상시키기 위해 이동하게 되나, connectivity는 유지하게 된다. Spring-based smoothing 기법은 동적격자 생성을 위해 cell의 모서리(node와 node 사이)를 서로 연결된 spring 조합으로 고려하는 개념이다.

물체 초기 형상을 고려하여 생성된 격자계에 대해 물체가 이동 또는 변형되면 우선 표면에 위치한 격자점들을 물체 면을 따라 강제로 이동시키게 된다. 이후 이 이동변위를 spring force로 등가 고려하여 구조 평형상태에서 모든 spring 하중의 합은 0이 되어야 하는 조건을 반영하여 아래와 같이 반복계산(iterative) 방정식이 얻어질 수 있다. 또한 임의의 i 격자점에 대해 물체 이동 및 회전 변형에 의한 모든 격자점들의 신규 변위(displacement) 벡터는 다음과 같이 나타낼 수 있다.

δ → i n + 1 = ∑ j = 1 n i k i j ⋅ δ n → ∑ j = 1 n i k i j

(7)

여기서, $δ → i n + 1 = Δ x n + 1 i → + Δ y n + 1 j → + Δ z n k →$ 이고, k_ij는 각각 j 격자점들에 연결되어 있는 i 격자점의 신규 변위벡터와 스프링 강성을 의미한다.

2.3 동적격자 보존방정식

동적격자 재생성 기법을 적용하는 경우, 초기 구성된 CFD 격자 cell 단위에서 격자 경계면의 움직임으로 인해 상대속도가 유발되어 보존방정식(conservation equation)의 보정이 필요하다. 동적격자와 관련하여 경계가 움직이는 임의의 제어체적(V)에서 일반적인 스칼라(φ)에 대한 보존방정식은 다음과 같다.

d d t ∫ V ρ ϕ d V + ∫ ∂ V ρ ϕ (u → − u g →) ⋅ d A → = ∫ ∂ V Γ ∇ ϕ ⋅ d A → + ∫ V S ϕ d V

(8)

여기서, ρ는 유체의 밀도, $u →$ 는 유체 속도벡터, $u g →$ 는 격자 속도벡터, Γ는 디퓨전(diffusion) 계수, S_φ는 ϕ의 source 항을, ∂V는 제어체적(V)의 경계를 의미한다.

1차 후방차분(backward difference) 공식을 사용하여 시간미분 항은 다음과 같이 나타낼 수 있다.

d d t ∫ V ρ ϕ d V = (ρ ϕ V) n + 1 − (ρ ϕ V) n Δ t

(9)

여기서 n과 n+1은 현재 및 다음 시간스텝을 의미하며, n+1 시간스텝에서 체적 Vⁿ⁺¹은 다음과 같이 계산된다.

V n + 1 = V n + d V d t Δ t

(10)

여기서 dV/dt는 제어체적의 시간미분을 의미하며, 격자 보존법칙을 만족시키기 위하여 다음식과 같이 계산된다.

d V d t = ∫ ∂ V u g → ⋅ d A → = ∑ ​ f n f u g j → ⋅ A j →

(11)

위 식에서 n_f는 제어체적의 면(face) 수를 나타내고, $A j →$ 는 j번째 면의 면적벡터를 의미한다. 각각의 제어체적 면에 대해 식 (8)의 우변 항은 다음 식으로 계산된다.

u g j → ⋅ A j → = δ V j Δ t

(12)

여기서 δV_f는 모든 시간스텝에서 모든 제어체적 면 j에 대해 계산되게 된다.

2차 정확도(second order)의 후방차분 공식을 사용하여 식 (8)의 시간미분 항은 다음과 같이 나타낼 수 있다.

d d t ∫ V ρ ϕ d V = 3 (ρ ϕ V) n + 1 − 4 (ρ ϕ V) n + (ρ ϕ V) n − 1 2 Δ t

(13)

위와 같이 2차 정확도의 시간차분법을 적용하는 경우 제어체적의 부피 시간미분은 각 제어체적 면에 대해 다음과 같이 계산된다.

(u g j → ⋅ A j →) n + 1 = 32 (u g j → ⋅ A j →) n − 12 (u g j → ⋅ A j →) n − 1 = 32 (δ V j δ t) n − 12 (δ V j δ t) n − 1

(14)

여기서, (δV_j)ⁿ과 (δV_j)^n-1은 모든 시간스텝에서 모든 제어체적 면 j에 대해 현재 시간스텝 및 이전 시간스텝에 대해 계산되게 된다.

2.4 이동물체 6-자유도계 운동해석

외부장착물 투하 후 운동궤적을 예측하기 위해서는 CFD 유동해석 시 투하된 장착물의 6-자유도계 동역학 운동방정식을 연성한 수치해석이 수행되어야 한다. 운동하는 장착물에 대해 임의의 시간스텝에서 계산된 공력 하중과 모멘트 및 중력 하중을 6-자유도계 동역학 운동방정식에 대입하여 다음 시간스텝에서 위치 및 자세를 결정하게 된다. 이후 격자 재생성을 수행하고, CFD 해석을 통해 다시 공력하중을 구하는 과정을 반복하면서 해석을 수행하게 된다.

장착물 무게 중심의 병진 운동을 해석하기 위한 동역학 지배 방정식은 다음과 같다.

v ˙ G → = 1 m Σ f G →

(15)

여기서, $υ ˙ G →$ 는 물체 무게 중심의 병진운동 속도벡터를 나타내고, m은 질량, $f G →$ 는 물체 무게 중심에 작용하는 적분된 공력하중 및 중력하중 벡터를 의미한다.

물체의 회전 동역학 운동은 국부 좌표계에 대해 동역학적 모멘트 평형방정식을 이용하여 다음과 같이 계산할 수 있다.

ω ˙ B → = 1 I G (Σ M B → − ω B → × I G ω B →)

(16)

여기서 I_G은 질량관성 텐서, $M B →$ 는 물체에 작용하는 모멘트 벡터, $ω B →$ 는 각속도 벡터를 의미한다. 모멘트는 관성 좌표계에서 물체 좌표계로 다음과 같이 변환될 수 있다.

M B → = R M G →

(17)

여기서, R은 변환행렬로 다음과 같다.

[C θ C ψ C θ S ψ − S θ S φ S θ C ψ − C φ S ψ S φ S θ S ψ + C φ C ψ S φ C θ C φ S θ C ψ + S φ S ψ C φ S θ S ψ − S φ S ψ C φ C θ]

(18)

위 식에서 C_χ=cos(χ) and S_χ=sin(χ)를 의미한다. 또한 φ는 z축에 대한 회전각(yaw), θ는 y축에 대한 회전각(pitch), ψ는 x축에 대한 회전각(roll)으로 Euler angle을 나타낸다. 식 (7)과 식 (9)에 대한 수치적분을 통해 병진운동 속도 및 회전속도 벡터를 구할 수 있으며, 이는 다음 시간스텝의 동적격자(dynamic mesh) 재생성 과정에 사용되게 된다.

Ⅲ. 해석결과 및 검토

3.1 해석 모델

본 연구에서 수치실험을 위해 Arnold Engineering Development Center[2]에서 수행한 1/20 축소모델 풍동시험 모델 형상을 대상으로 하였다. 또한 참고문헌에서 제시한 외부장착물의 무게와 ejector의 사출 하중을 모사하기 위하여 실제 크기로 형상 모델링을 수행하였다(Fig. 2). 날개의 익형은 NACA 64A010이며, 앞전 후퇴각은 45°, 날개 시위길이(chord length)는 7,620 mm이다. 파일런의 위치는 날개 스팬 길이의 중간 지점에 위치하고 있다. 4개의 외부장착물 핀의 익형은 NACA 0008이다. Table 1에는 본 해석에 필수적으로 필요한 외부장착물의 질량, 질량중심, 질량관성모멘트, 파일런에 위치한 2개의 store ejector 위치 및 ejector force 등에 대한 상세 수치를 제시하였다.

Fig. 2. Geometrical shape of the generic wing-pylon-store model (unit: mm)[]

Download Original Figure

Table 1. Dynamic properties of the store and ejector

Parameter↩	Value↩
Weight↩	907 kg↩
Center of Gravity↩	1417 mm (aft of the store nose)↩
Roll Inertia↩	27.116 kg · m²↩
Pitch Inertia↩	488.094 kg · m²↩
Yaw Inertia↩	488.094 kg · m²↩
Roll Damping Coefficient↩	-4/rad↩
Pitch Damping Coefficient↩	-40/rad↩
Yaw Damping Coefficient↩	-40/rad↩
Forward Ejector Location↩	1237.49 mm↩
Forward Ejector Force↩	10675.7 N↩
Aft Ejector Location↩	1746.50 mm↩
Aft Ejector Force↩	42702.7 N↩
Ejector Stroke Length↩	100.584 mm↩

Download Excel Table

항공기 날개에서 외부장착물의 시간에 대한 분리운동해석을 위해서는 전산유체역학(CFD) 및 다물체동역학(Multi-body Dynamics, MBD) 연성해석 기법 적용이 요구된다.

Fig. 3은 CFD-MBD 연성해석을 위해 생성된 3차원 CAD 모델 형상을 보여주고 있다. 좌표계정의는 유동 방향이 x축, 날개 스팬 방향이 y축 그리고 수직 방향이 z축이며, 원점은 외부장착물의 질량중심에 위치시켰다.

Fig. 3. Generated CAD geometry for the generic wing-pylon-store model

Download Original Figure

Fig. 4는 Fig. 3의 형상에 대한 3차원 CFD 유동해석 격자를 보여주고 있다. 외부 유동장의 크기는 유동 방향으로 날개 뿌리 시위길이(root chord length)의 25배, 높이 방향으로 20배, 스팬 방향으로 5배를 설정하였다. CFD 격자 유형은 사면체 tetrahedron grid이며, 총 cell 수는 3,858,000개이다. 본 모델의 경우, 외부장착물의 분리운동이 고려되므로 파일런과 외부장착물 사이에 미소 공간(gap)이 존재하게 된다. 장착물 분리 운동 시 수치안정성 확보를 위해 미소 공간에 수직방향으로 최소 4개 이상의 격자가 분포되도록 격자 분포를 조절하여 생성하였다. 이를 위해 별도의 body of influence (BOI) 영역을 생성하였으며, 이 영역에서는 격자의 밀집도 함수를 공간격자와 개별적으로 조정할 수 있도록 하였다.

Fig. 4. Tetrahedral CFD grid domain of the generic wing-pylon-store (Total cells = 3,858,000)

Download Original Figure

Fig. 5는 설정된 유동해석 경계조건을 보여주고 있다. 날개-파일런-외부장착물 표면은 wall 경계조건을, 날개 뿌리에서 x-z 면은 symmetry경계조건을, 외곽경계는 far field 경계조건을 설정하였다.

Fig. 5. CFD boundary conditions for the generic wing-pylon-store model

Download Original Figure

CFD-MBD 연성해석 과정에서 2개의 사출기 반발력을 받아 외부장착물이 분리운동을 시작하게 되면 CFD 유동영역 내에서 장착물의 위치와 자세가 실시간으로 변화하게 된다. 장착물의 운동은 장착물 표면에 작용하는 외부압력 하중을 실시간으로 고려한 6-자유도계(6-DOF) 동역학 평형방정식을 풀어서 결정하게 된다. 이 과정에서 비정상(unsteady) 유동해석을 정교하게 수행하기 위해서는 매 시간 스텝마다 장착물의 신규 위치 및 자세를 반영한 3차원 유동해석 격자의 재생성 과정이 필수적으로 필요하다. CFD 유동해석에서 장착물의 운동을 비교적 용이하게 고려할 수 있는 중첩격자법(overset mesh method)이 있지만, 이 방법은 중첩격자 영역에서 매 시간스텝 보간 오차(interpolation error)가 누적될 수 있는 단점이 있다. 본 논문에서는 중첩격자법을 적용하지 않고 외부장착물의 운동을 해석할 수 있도록 spring analogy[13] 기반의 이동변형격자(moving and deforming mesh, MDM) 기법을 활용하였다. MDM 공간격자 재생성 기법은 중첩격자법의 단점인 수치보간 오차는 배제할 수 있지만, 격자재생성 과정에서 격자의 꼬임 등에 의한 수치불안정성이 발생할 수 있는 단점이 있다. MDM 기법 적용 시 관계되는 주요 입력변수는 spring constant factor(S.C.F.), convergence tolerance (C.T.) 및 max iteration이다. 본 논문에서는 마하 0.95 및 1.2 운항조건에서 이러한 수치 매개변수가 외부장착물 분리운동 해석결과에 미치는 영향을 고찰하였다.

3.2 Mach 0.95 운항조건 외부장착물 분리운동 해석

천음속 영역인 마하수 0.95에서 외부장착물 분리운동 해석을 위해 고려한 운항조건은 고도 26,000ft, 대기압력 35977.12 Pa 및 대기온도 236.54 K 이다. MDM 공간격자 재생성 기법의 초기 매개변수 설정은 S.C.F. 0.5 및 C.T. 0.001이며, 해석결과, 개선을 위해 재설정한 매걔변수 값은 S.C.F. 0.1, C.T. 0.0005이다.

Fig. 6은 운항 마하수 0.95에서 외부장착물 분리운동해석을 수행한 결과로, 시간에 대한 외부장착물의 위치 및 자세와 표면 마하수 분포 가시화 결과를 보여주고 있다. 2개의 Ejector 반발력을 작용한 이후 장착물의 6-자유도계 운동이 안정적으로 수치해석되었음을 확인할 수 있다. CFD-MBD 연성해석을 위한 시간간격(time-step)은 0.001초이며, 장착물의 실시간 운동이 반영된 비정상(unsteady) CFD 유동해석 과정에서 시간정확도 향상을 위한 sub-iteration 회수는 총 20회를 설정하였다.

Fig. 6. Store motion and surface Mach number contour during separation motion (M = 0.95, h=26,000 ft)

Download Original Figure

Fig. 7은 마하 0.95 운항조건에서 MDM 기법 매개변수 설정을 S.C.F. 0.5 및 C.T. 0.001로 한 경우 분리된 외부장착물의 이동변위 및 자세를 풍동시험 결과와 비교한 결과이다. CFD-MDB 연성해석 결과와 시간에 대한 3축 방향 질량중심 이동 변위는 모두 풍동시험 결과와 잘 일치하고 있음을 확인할 수 있다. 장착물 분리운동 시 회전각 자세의 경우는 피치(pitch)와 요(yaw) 방향으로는 기존 풍동시험 결과와 잘 일치하는 결과를 나타내고 있으나, 롤링(rolling) 각도의 경우만 다소 차이를 나타내고 있다.

Fig. 7. Comparison of store C.G. position and rotation angle with wind-tunnel test data (M=0.95, h=26,000 ft, MDM Initial parameter)

Download Original Figure

Fig. 8은 마하수 0.95 조건에서 분리된 장착물이 운동하는 동안 0초, 0.17초 및 0.32초 순간에서 외부장착물 표면에서 몸체 길이 방향으로 압력계수 분포를 비교한 결과이다. 본 연구에서의 수치해석 결과가 기존 풍동시험 결과와 잘 일치하는 것을 확인할 수 있다.

Fig. 8. Comparison of the pressure coefficients distribution along circumferential angle 5° at different time steps

Download Original Figure

Fig. 9는 Fig. 7과 동일한 조건인 마하 0.95에서 MDM 매개변수 설정을 S.C.F. 0.5 및 C.T. 0.001에서 S.C.F. 0.1 및 C.T. 0.0005로 변경하여 해석한 경우의 결과를 보여주고 있다. 장착물 무게중심 운동변위의 경우, 변경 전 결과와 거의 유사한 경향을 나타내고 있으며, 회전변위 운동응답의 경우 pitch 응답이 실험결과에 약간 더 근접한 결과를 나타내었다.

Fig. 9. Comparison of store CG position and rotation angle at M=0.95 for the different spring control factors

Download Original Figure

3.3 Mach 1.2 초음속 운항조건 외부장착물 분리운동 해석

초음속 영역인 마하수 1.2에서 외부장착물 분리운동 해석을 위해 고려한 운항조건은 고도 38,000ft, 대기압력 20,589 Pa 및 대기온도 216.54 K이다. Fig. 10은 운항 마하수 1.2에서 외부장착물 분리운동해석을 수행한 결과로 시간에 대한 외부장착물의 위치 및 자세와 표면 마하수 분포 가시화 결과를 보여주고 있다. 마하 0.95의 경우와 마찬가지로 파일런에서 2개의 ejector 반발력을 가한 후, 장착물의 6-자유도계 운동이 안정적으로 해석되었음을 확인할 수 있다. CFD-MBD 연성해석을 위한 시간간격(time-step)은 0.001초이며, 비정상 유동해석 과정에서 시간정확도 향상을 위한 sub-iteration 회수는 총 20회를 설정하였다. 장착물 투하 조건에 대한 총 해석시간은 0.8초이다. 마하 0.95의 경우와 비교해 보면 비행 마하수 증가로 인한 외부장착물 주변의 유동환경 변화로 동일한 ejector 반발력 조건에서 운동궤적이 변화됨을 확인할 수 있다.

Fig. 10. Store motion and surface Mach number contour during separation motion (M=1.20 and h=38,000 ft)

Download Original Figure

Fig. 11에서는 마하 1.2 운항조건에 대해 MDM 초기 매개변수를 S.C.F. 0.5 및 C.T. 0.001로 설정하여 수치해석을 수행한 경우, 장착물의 운동변위 및 회전운동 궤적을 풍동시험 결과와 비교하였다. 초음속 유동의 경우 z축 방향으로의 운동변위는 풍동시험 결과를 잘 예측하는 반면에 나머지 자유도 운동응답은 큰 차이를 보이고 있다.

Fig. 11. Comparison of store C.G. position and rotation angle with wind-tunnel test data (M=1.2, h=38,000 ft, MDM Initial parameter)

Download Original Figure

Fig. 12는 외부장착물 분리 전 동일한 CFD 격자에 대해 분리 후 0.3초가 지난 시점에서 초기 MDM 변수설정과 개선된 MDM 변수설정 값에 따라 다시 생성된 격자분포를 보여주고 있다. 그림에서는 큰 차이가 없는 것처럼 보이나, 투하 운동하는 외부장착물 주변으로 비정상 유동해석 시간 스텝마다 재생성되는 격자분포 및 품질이 MDM 변수설정 값에 따라 차이가 발생하게 된다. 본 연구에서는 다양한 MDM 변수 값에 대한 수치실험을 수행하였으나, 지면 관계상 모든 결과를 제시하지는 못하였다.

Fig. 12. Comparison of regenerated grids using different MDM spring parameters (cutting plane view, M=1.20)

Download Original Figure

Fig. 13은 Fig. 11과 동일한 해석조건의 경우로 마하 1.2에서 개선된 MDM 매개변수를 적용한 경우, 외부장착물 운동궤적 해석결과를 풍동시험 결과와 비교한 경우이다. 개선된 MDM 매개변수를 적용한 경우 더욱 잘 일치하는 결과가 도출됨을 확인할 수 있다. 게다가 마하 1.2의 비행조건에서는 MDM 변수설정에 따라 재생성되는 주변 격자의 영향이 마하 0.95의 경우보다 민감하게 작용하게 됨을 파악할 수 있다. 마하 0.95와 마하 1.2 조건에서 발생하는 충격파의 상호작용 특성이 다르고, 재생성되는 공간격자 분포가 장착물의 운동을 결정하는 공력하중 결과에 지속적인 되먹임(feedback) 작용으로 영향을 미치면서 누적되기 때문이다. 이러한 수치실험 결과는 기존 연구에서 보고된 사례가 없었다. 동적이동격자 CFD 해석기법으로 외부장착물 분리운동을 해석하는 경우, 초음속 유동에서 MDM 변수설정 조건에 따라 더욱 민감하게 영향을 미칠 수 있음이 파악되었다. 본 논문의 결과는 향후 관련된 수치해석 연구의 정확도 향상을 위해 유용한 활용이 기대된다.

Fig. 13. Comparison of store CG position and rotation angle at M=1.2 for the different spring control factors

Download Original Figure

IV. 결 론

본 논문에서는 충격파의 상호작용이 존재하는 천음속 및 초음속 운항조건에서 항공기 외부장착물 CFD 분리운동 해석을 MDM 기법으로 수행 및 비교하였다. 마하 0.95 및 1.2의 비행속도 조건에 대해 기존 풍동시험 결과와 잘 일치하는 결과를 얻을 수 있었다. 천음속 영역보다 초음속 영역에서 MDM 기법의 매개변수가 장착물 운동궤적 해석 결과에 더욱 민감한 영향을 미칠 수 있음을 파악하였으며, 이의 개선을 통해 더욱 정확한 해석 결과를 도출할 수 있음을 보였다. 본 연구의 결과는 향후 다양한 종류의 외부장착물 투하해석을 위한 MDM 기법 적용 시 정확한 투하궤적 결과를 도출하는 데 유용하게 활용될 수 있을 것으로 기대된다.

후 기

이 성과는 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 연구임(No. NRF-2017R1A2B4010797).

References

[1].

Carman, J. B. Jr., Hill D. W. Jr., and Christopher, J. P., Store Separation Testing Techniques at the Arnold Engineering Development Center, Volume II, Description of Captive Trajectory Store Separation Testing in the Aerodynamic Wind Tunnel (4T), AEDC-TR-79-1, June 1980.

[2].

Heim. E. R.,“CFD wing/pylon/finned store mutual interference wind tunnel experiment,” Arnold Engineering Development Center, Arnold Air Force Base, Tennessee, AEDC-TSR-91-P4, February 1991.

[3].

Lijewski, L. E. and Suhs, N. E., “Time-accurate computational fluid dynamics approach to transonic store separation trajectory prediction,” Journal of Aircraft, (31)4, July-August 1994, pp. 886-891.

[4].

Prewitt, N. C., Belk, D. M., and Maple, R. C., “Multiple body trajectory calculations using the begger code,” Journal of Aircraft, (36)5, September-October 1999.

[5].

Cline, D. M., Riner, W. C., Jolly B. A., and Lawrence, W. R., “Calculation of generic store separations from an F-16 aircraft,” AIAA-96-2455-CP, 1996.

[6].

Coleman, L., Jolly, B., Chesser, B. L. Jr., and Brock, J. M. Jr., “Numerical simulation of a store separation event from an F-15E aircraft,” AIAA Paper 96-3385, July 1996.

[7].

Tomaro, R. F., Witzeman, F. C., and Strang, W. Z., “A solution on the F-18C for store separation simulation using cobalt60,” AIAA Paper 99-0122, 1999.

[8].

Lee, S., Park, M., Cho, K. W., and Kwon, J. H., “New fully automated procedure for the prediction of store trajectory,” Journal of Aircraft, 37(6), Novermber-December 2000.

[9].

Panagiotopoulos, E. E., and Kyparissis S. D., “CFD transonic store separation trajectory predictions with comparison to wind tunnel investigations,” International Journal of Engineering (IJE), (3)6, January 2010, pp. 538-553.

[10].

Jafari, M., Toloei, A., Ghasemlu, S., and Parhizkar, H., “Simulation of store separation using low-cost CFD with dynamic meshing,” International Journal of Engineering, (27)5, May 2014, pp. 775-784.

[11].

Zheng, J., Chen, J., Zheng, Y., Yao, Y., Li, S. and Xiao, Z., “An improved local remeshing algorithm for moving boundary problems”, Engineering Applications of Computational Fluid Mechanics, (10)1, 2016, pp. 403-426.

[12].

Huang, H., Young, M., Giacobello, M., and Blyth, R., “A comparison of approaches to multi-body relative motion using the Kestrel CFD solver,” 21st Australasian Fluid Mechanics Conference, Adelaide, Australia, 10-13 December 2018.

[13].

Kim, D. H., Park, Y. M., Lee I., and Kwon O. J., “Nonlinear aeroelastic computation of a wing/pylon/finned-store using parallel computing,” AIAA Journal, (43)1, January 2005, pp.53-62.

[14].

Han, L., Kim, D. H., Kim, Y. S., Kim, Y. H., and Kim, S. S., “Transonic flutter characteristics of supercritical airfoils considering shockwave and flow separation effects”, Journal of the Korean Society for Aviation and Aeronautics, (17)2, June 2009, pp. 8-17.

[15].

Kim, Y. H., and Kim, D. H., “Transonic flutter characteristics of the AGARD 445.6 Wing Considering DES Turbulent Model and Different Angle-of-Attack”, Journal of the Korean Society for Aviation and Aeronautics, (18)1, March 2010, pp. 27-32.