# Investigation of Aeroelasticity Effect on Missile Separation from the Internal Bay

Authors: Shuling Tian; Rongjie Li; Ke Xu
Journal: International Journal of Aerospace Engineering
 (2023)

---

## Abstract

There is a strong aerodynamic interference when launching the missile in the embedded mode. During the separation process, the carrier aircraft safety may be threatened due to large slenderness ratio, low structural stiffness, and aeroelasticity effects of the missile. The present study simulates missile separation in the presence of the aeroelasticity effects based on the computational fluid dynamics (CFD), rigid body dynamics (RBD), and computational structure dynamics (CSD) coupling method. A hybrid dynamic grid method consisting of the mixed overset unstructured grid and deformation grid is utilized. In order to verify the accuracy of the coupled numerical method, store separation from a wing and AGARD 445.6 wing flutter are first simulated as two standard test cases. The verification results imply that the present coupled numerical method is reliable and capable in simulation of the aeroelastic effect in missile separation. The influence of aeroelasticity on the separation trajectory of a missile from the internal bay is systematically studied at different states. Numerical results show that aeroelasticity substantially affects the missile angular displacement, while it has a slight impact on the linear displacement of the center of mass. Mach number and flight altitude are two important flight parameters that characterize the aeroelasticity effect on missile separation from the internal bay.


---

## Body

## 1. Introduction

Relative motion between multiple bodies is a common problem in aviation and aerospace technology which causes unsteady flow due to moving boundaries. Separations of launch vehicles, aircraft mounts, and ejection seats are some examples in this regard. The multibody separation problem is significantly characterized by the coupling between flow field and motion. This leads to strong aerodynamic interference between moving objects, which is often nonlinear and unsteady [1–3]. Although the embedded missile internal bay can ensure the stealth and mobility of the fighter, it causes new problems in releasing the missile. After opening the internal bay door, complex unsteady flow phenomena such as boundary layer separation and reattachment, shock wave, and boundary layer interference are likely to appear in the high-speed air flow through the missile internal bay. They result in an unstable state when the missile and fighter are separated [4, 5]. However, air-to-air missiles generally involve large slenderness ratios, such that their structural stiffness is small and the elastic effect becomes prominent especially with the wide application of new composite materials. The complicated unsteady flow near the missile internal bay and the interactions due to the elastic effects deteriorate the launch performance of the missile and even threaten the carrier safety [6]. Therefore, it is necessary to study the influence of aeroelasticity on the separation trajectory of missile when it is released from the embedded missile bay.Multibody separation has been successfully investigated in numerous studies by means of flight test, wind tunnel test, and numerical simulation [7–10]. However, it is difficult to consider the aeroelastic effect in the multibody separation in the wind tunnel test, and limited research has been conducted by this approach. Over the past several decades, computational fluid dynamics has considerably developed, and thanks to its low cost and high accuracy, it has been utilized as a crucial technology in the development of advanced aerospace vehicles. It has been used in the simulation of unsteady process of multibody separation and aeroelasticity. Izhak et al. [11] studied the influence of elastic wing on the release of external stores. They indicated that the rolling motion of external stores caused by the bending deformation of the wing affects the separation safety. Zhang et al. [12] investigated the aeroelastic characteristics of the wing during the separation of external stores. Wang et al. [13] used the loose coupling method to analyze the influence of elastic and rigid wings on the separation trajectory of external stores and found quite different trajectories for these two cases. Yang et al. [14] simulated the release of external stores for elastic wings based on the modal superposition method. They revealed that the elastic deformation of wings at the initial stage of separation greatly affects the trajectory of external stores. However, these works have only considered the wing elasticity, not the missile elasticity. Hua et al. [15] studied the influence of aeroelasticity on external missile separation trajectory based on the overset unstructured grid. It was observed that elastic deformation has a great impact on the missile attitude angle, and elastic deformation can reduce the attitude angle oscillation. Zhu et al. [16] numerically simulated launching of airborne missile with either elastic or rigid bodies. They indicated that the elastic deformation of the missile slightly affects the displacement while it greatly characterizes the attitude angle. Although the aeroelasticity of the released missile has been considered in the simulation of separation in the abovementioned studies, they have concerned the external missiles. However, there is a stronger aerodynamic interference in the separation of missile from the internal bay. In addition, it involves diverse characteristics at different Mach numbers and flight altitudes. Therefore, the aeroelastic effect on the separation of missile from the internal bay should also be thoroughly studied.The present study utilizes a hybrid dynamic grid to accurately study the impact of aeroelasticity on missile separation from the internal bay. In this hybrid approach, the overset unstructured grid is utilized to simulate the relative motion between separated bodies and the deformation grid is adopted to handle the aeroelastic deformation of the structure. Based on the hybrid dynamic grid, the unsteady flow equations are solved by coupling the six-degree-of-freedom (6-DoF) motion equations and structural dynamics equations.The remainder of this paper is organized as follows. The second section thoroughly describes the numerical methodology for the problem of multibody separation with consideration of aeroelastic effect. It introduces unsteady CFD solver, 6-DoF motion model, structural dynamics solver, hybrid dynamic grid method, and the simulation flowchart. Section3 presents the test results for verification of the numerical method. The subsequent sections investigate aeroelastic effects on the separation process of missile from the internal bay and then present the concluding remarks.


## 2. Numerical Method of Multibody Separation with Aeroelasticity

The aeroelasticity in the separation process of missile from the internal bay is a coupled multiphysics problem. It can be simulated by solving the coupled equations of unsteady flow and motion of free elastic body. The motion of an elastic body under large rigid motion can be decomposed into a large rigid body displacement and a small elastic displacement with respect to some inertial frame [17]. Using mean axis as the reference frame can minimize dynamic coupling between the elastic and reference frames [18]. Under the assumption of small deformation, the mean axis can be approximately regarded as the principal axis of the moving deformable body and is considered to remain unchanged. The information of center of mass and angular orientation of this reference frame for the deformable body can be directly determined from the solution of the equations of motion [19]. In solving the structure dynamics equation in the mean axis reference frame, the inertia forces of each finite element node due to the large rigid motion are obtained from their acceleration which can be determined from the solution of the equation of rigid body motion. The unsteady flow solver, rigid body motion model, structural solver, and the dynamic grid method used to treat the moving boundary are presented in this section.


### 2.1. Unsteady Flow Solver

The conservative integral form of unsteady compressible Reynolds-averaged Navier-Stokes (RANS) equations in three-dimensional rectangular coordinate system can be written as [20]
(1)∂∂t∫ΩW⟶dΩ+∮∂ΩFc⟶−Fv⟶dS=0,where Ω is the control volume domain and ∂Ω represents the control surface domain. W⟶, Fc⟶, and Fv⟶ represent the conserved variable vector, the convective flux vector, and the viscous flux, respectively, which are defined as follows [21]:
(2)W⟶=ρρuρvρwρE,Fc⟶=ρVρuV+nxpρvV+nypρwV+nzpρHV+pVb,Fv⟶=0nxτxx+nyτxy+nzτxznxτyx+nyτyy+nzτyznxτzx+nyτzy+nzτzznxΘx+nyΘy+nzΘz,where ρ represents fluid density, p is static pressure, and E and H indicate total energy and total enthalpy per unit mass, respectively. u, v, and w notify velocity components and nx, ny, and nz represent the unit normal vector in the outward direction of the boundary in the x, y, and z directions, respectively. V indicates the contravariant velocity, which is defined as the scalar product of velocity vector and the unit normal vector:
(3)V≡v⟶−v⟶g·n⟶=u−ugnx+v−vgny+w−wgnz,where v⟶g=ugi⟶+vgj⟶+wgk⟶ is the grid motion velocity. The finite volume method based on unstructured grid is used to discretize the three-dimensional unsteady compressible RANS equations. For an arbitrary control volume Ωm, Equation (1) can be discretized as follows:
(4)dW⟶mdt=−1Ωm∑n=1NFFc⟶n∆Sn=−1ΩmR⟶m,where NF represents the number of faces over the control body surface, ∆Sn indicates the area of the nth face, and Rm notifies the residual. The convective flux vector Fc⟶ is computed based on the Riemann scheme of HLLC/E [21].The dual time-stepping scheme proposed by Jameson [22] is adopted to march in the physical time. Accordingly, Equation (4) can be expressed as follows:
(5)∂ΩIn+1WI∗∂τ=3WI∗ΩIn+1−4WInΩIn+WIn−1ΩIn−12∆t+RIW∗=−RI∗W∗,whereτ is a pseudo-time step which is marched based on the implicit LU-SGS [23] scheme.


### 2.2. Rigid Body Dynamics Motion Model

The trajectory of a separation body is determined by the force and moment imposed on that body. Therefore, the trajectory and attitude angle of the body in the process of multibody separation are determined by solving the six-degree-of-freedom equation of the rigid body. The motion of a rigid body can be decomposed into translation of its center of mass and rotation around the center of mass [24]. The six-degree-of-freedom motion equation of the rigid body can be obtained based on this decomposition. According to Newton’s law, the motion equation of the center of mass in the inertial coordinate system can be expressed as follows:
(6)md2xcmidt2=Fxi,(7)md2ycmidt2=Fyi,(8)md2zcmidt2=Fzi,where xcmi, ycmi, and zcmi represent the coordinates of the center of mass and Fxi, Fyi, and Fzi indicate the force components in the x, y, and zdirections, respectively. The superscript i indicates that the equation is established on the inertial coordinate system.The equation of rigid body rotation around the center of mass is expressed as follows:(9)Ixxbdωxbdt=Mxb+Iyyb−Izzbωybωzb,(10)Iyybdωybdt=Myb+Izzb−Ixxbωzbωxb,(11)Izzbdωzbdt=Mzb+Ixxb−Iyybωxbωyb,where subscript b represents the body frame, Mxb, Myb, and Mzb are the moments relative to the center of mass, ωxb, ωyb, and ωzb indicate the angular velocity, and Ixxb, Iyyb, and Izzb are the moments of inertia of the rigid body.The quaternion method is used to solve the equation of motion, which is expressed as follows:(12)P=e0e1e2e3.The transformation matrix composed of directional cosines can be expressed as(13)A=2e02+e12−0.5e1e2−e0e3e1e3+e0e2e1e2+e0e3e02+e22−0.5e2e3−e0e1e1e3−e0e2e2e3+e0e1e02+e32−0.5.The direction change caused by angular velocity can be calculated as follows:(14)dPdt=12Lωb,where
(15)L=−e1e0e3−e2−e2−e3e0e1−e3e2−e1e0.The force in Equation (7) can be decomposed into three parts, namely, the aerodynamic force, the gravity force, and the external force. The same decomposition can be performed for the moment in Equation (10). The aerodynamic force and moment are calculated based on CFD. The position and direction of the rigid body are updated by numerical integration of Equations (7), (10), and (14) using the predictor-corrector scheme [25]. The new position of the rigid body can be obtained as follows:
(16)rin+1=rc.m.in+1+An+1rb−rc.m.b,where rin+1 and rc.m.in+1 represent the coordinates of the rigid body and the center of mass in the inertial frame at the time step n+1, respectively. rb and rc.m.b indicate the coordinate of the rigid body and the center of mass in the body frame.


### 2.3. Structural Dynamics Solver

The structural dynamics system can be expressed as a system of second-order ordinary differential equations [26] as follows:
(17)Mu¨t+Cu̇t+Kut=Ft,where M, C, and K are the mass, damping, and stiffness matrices, respectively. u¨t, u̇t, and ut denote structural acceleration, velocity, and displacement vectors, respectively. F is the difference of the aerodynamic loads and inertial forces which are determined from the solution of equations of rigid body motion.In the present study, the Newmark method [27] is used to discretize the structural dynamics equation. Equation (17) can be written as follows:
(18)Mu¨n+1+Cu̇n+1+Kun+1=Fn+1.In addition to Equation (17), the Newmark method updates the displacement and velocity vectors as follows:
(19)u̇n+1=u̇n+1−δu¨n+δu¨n+1∆t,(20)un+1=un+u̇n∆t+12−αu¨n+αu¨n+1∆t2,whereδ and α are the Newmark integral parameters. The Newmark time integration algorithm can be determined based on the Newmark integration parameters.From Equations (18), (19), and (20), the one-step integration method can be expressed as follows:
(21)a0M+a1C+Kun+1=Fn+1+Ma0un+a2u̇n+a3u¨n+Ca1un+a4u̇n+a5u¨n,where
(22)a0=1α∆t2,a1=δα∆t,a2=1α∆t,a3=12α−1,a4=δα−1,a5=∆t2δα−2.un+1 can be calculated from Equation (21), and u̇n+1 and u¨n+1 can be obtained as follows:
(23)u̇n+1=a1un+1−un−a4u̇n−a5u¨n+1,u¨n+1=a0un+1−un−a2u̇n−a3u¨n+1.According to the stability analysis, the method is unconditionally stable when the Newmark parameters satisfy the following conditions:(24)δ≥12,α≥1412+δ2.


### 2.4. Hybrid Dynamic Grid Method

In the problem of missile release from the internal bay with consideration of aeroelastic effect, both the structure deformation and the overall motion relative to the aircraft are involved. The dynamic grid method can simultaneously treat both of them and is practical to simulate such complex process. A hybrid dynamic grid method is developed, in which the deformation grid deals with the aeroelastic deformation of structure and the overset unstructured grid handles the relative motion between separation bodies.


#### 2.4.1. Deformation Grid Method

When considering the aeroelasticity effect, the structural deformation should be concerned in the process of release and separation, and this is simulated by the structure dynamics solver. In order to capture the flow after the structure deformation, the CFD near-wall grid should adapt to the deformation. Therefore, the deformable grid method is adopted in the current study, because it can deal with the deformable wall surface. The basic idea of the deformable mesh method [28] is to keep the mesh topology unchanged and absorb the movement of the solid wall boundary by adjusting the nodes position in the field grid. Although this method is only applicable for small-scale motion, it has the advantages of preserving the correlation information between meshes and avoiding interpolation at each step. The spring-tension-based deformation method is utilized to adjust the field grid node position.The spring-tension method models the connection between each two grid nodes by a spring. The displacement of a given boundary node produces a force, which is directly proportional to the displacement of all springs connected to that node [29]. According to Hooke’s law, the force imposed on the mesh node can be expressed as follows:
(25)F⟶i=∑jnikij∆x⟶j−∆x⟶i,where ∆x⟶j and ∆x⟶i represent the displacement of node i and its adjacent node j and ni denotes the number of nodes nearby node i. kij indicates the elasticity coefficient between nodes i and j, which is defined as follows:
(26)kij=kfacx⟶i−x⟶j,where kfac is the spring constant. In the equilibrium condition, the resultant force exerted by all springs connected to the central node is zero, leading to the following iterative equation:
(27)∆x⟶im+1=∑jnikij∆x⟶jm∑jnikij,where m represents the number of iteration steps.After running the structure dynamics solver, the wall boundary displacement of the CFD grid is obtained by interpolating the structural grid. Once the boundary deformation is known, Equation (27) can be solved by the Jacobian scanning for all the field grid nodes.


#### 2.4.2. Dynamic Overset Unstructured Grid Method

In the missile release process, the missile moves relative to the aircraft and the grid must be adjusted to follow its movement. In the present study, a dynamic overset unstructured grid method is utilized to handle the moving boundary. This method was first introduced by Nakahashi [30], and it has the advantages of both unstructured grid and Chimera grid methods. The dynamic overset unstructured grid method is quite suitable to deal with unsteady flow past complex geometries and/or multiple bodies with relative movement. In the simulation of missile separation, subgrids are first generated around the missile and the carrier aircraft. The subgrids need to be generated only once, and the subgrids around the missile just move rigidly with its body during simulation of the whole separation process.In the present work, a wall distance-based overset unstructured grid assembly method is adopted. In order to obtain the hole region, this method applies a wall distance criterion to classify the grid nodes into three groups, i.e., active nodes, inactive nodes, and intergrid boundary nodes. A grid node is active when the distance between this node and its body surface is less than the distances to all other bodies. In contrast, a grid node is inactive if the distance to its own body is larger than the distances to the other bodies. The inactive nodes which are immediately next to the active ones can be further classified as the intergrid boundary nodes. The active nodes are used in the computations while the inactive nodes generate the hole region which is skipped in the computations. The intergrid boundary nodes, also referred to as interpolation nodes, are used to transfer the solution between the subgrids. The donor cells of the interpolation nodes can be efficiently found based on a reliable neighbor-to-neighbor search algorithm coupled with the alternating digital tree data structure [31, 32]. A general interpolation method [31] which is compatible with any element type is utilized for information transfer between the subgrids.In the unsteady simulation of missile separation, the overset unstructured grid assembly is applied with the new grid position before solving the flow equations. This is determined by superimposing the rigid and deformation motions.


### 2.5. Simulation Process of Multibody Separation with Aeroelasticity

By providing the flow solver, the rigid body dynamics model, the structure solver, and the dynamic grid, numerical simulation of missile separation in the presence of aeroelastic effects is achievable by coupling CFD/RBD/CSD on the dynamic grid. Under the assumption of small deformation, the influence of deformation on the mass distribution and rotational inertia of the elastomer is ignored. This implies that the rigid body motion and structure deformation is a one-way coupling. Figure1 presents the flowchart of numerical simulation scheme for this problem. The simulation procedure can be described as follows:
(1)
Initialization of the flow field for unsteady calculations based on the steady flow solution(2)
Calculation of the aerodynamic force and moment at the current physical time step by integrating the distributed force on the wall surface(3)
Transfer of the aerodynamic force and moment to the RBD model and, at the same time, interpolation of the distributed aerodynamic force (imposed on the wall boundary of CFD grid) to the CSD grid(4)
Solving the RBD equation and execution of the structure dynamics solver to obtain the whole motion and structure deformation(5)
Evaluate whether the solution is converged: if the number of iterations is insufficient, update the CFD grid by the overset unstructured grid and deformation grid methods; otherwise, terminate the simulationFigure 1
Simulation procedure for the problem of multibody separation in the presence of aeroelastic effect.


## 2.1. Unsteady Flow Solver

The conservative integral form of unsteady compressible Reynolds-averaged Navier-Stokes (RANS) equations in three-dimensional rectangular coordinate system can be written as [20]
(1)∂∂t∫ΩW⟶dΩ+∮∂ΩFc⟶−Fv⟶dS=0,where Ω is the control volume domain and ∂Ω represents the control surface domain. W⟶, Fc⟶, and Fv⟶ represent the conserved variable vector, the convective flux vector, and the viscous flux, respectively, which are defined as follows [21]:
(2)W⟶=ρρuρvρwρE,Fc⟶=ρVρuV+nxpρvV+nypρwV+nzpρHV+pVb,Fv⟶=0nxτxx+nyτxy+nzτxznxτyx+nyτyy+nzτyznxτzx+nyτzy+nzτzznxΘx+nyΘy+nzΘz,where ρ represents fluid density, p is static pressure, and E and H indicate total energy and total enthalpy per unit mass, respectively. u, v, and w notify velocity components and nx, ny, and nz represent the unit normal vector in the outward direction of the boundary in the x, y, and z directions, respectively. V indicates the contravariant velocity, which is defined as the scalar product of velocity vector and the unit normal vector:
(3)V≡v⟶−v⟶g·n⟶=u−ugnx+v−vgny+w−wgnz,where v⟶g=ugi⟶+vgj⟶+wgk⟶ is the grid motion velocity. The finite volume method based on unstructured grid is used to discretize the three-dimensional unsteady compressible RANS equations. For an arbitrary control volume Ωm, Equation (1) can be discretized as follows:
(4)dW⟶mdt=−1Ωm∑n=1NFFc⟶n∆Sn=−1ΩmR⟶m,where NF represents the number of faces over the control body surface, ∆Sn indicates the area of the nth face, and Rm notifies the residual. The convective flux vector Fc⟶ is computed based on the Riemann scheme of HLLC/E [21].The dual time-stepping scheme proposed by Jameson [22] is adopted to march in the physical time. Accordingly, Equation (4) can be expressed as follows:
(5)∂ΩIn+1WI∗∂τ=3WI∗ΩIn+1−4WInΩIn+WIn−1ΩIn−12∆t+RIW∗=−RI∗W∗,whereτ is a pseudo-time step which is marched based on the implicit LU-SGS [23] scheme.


## 2.2. Rigid Body Dynamics Motion Model

The trajectory of a separation body is determined by the force and moment imposed on that body. Therefore, the trajectory and attitude angle of the body in the process of multibody separation are determined by solving the six-degree-of-freedom equation of the rigid body. The motion of a rigid body can be decomposed into translation of its center of mass and rotation around the center of mass [24]. The six-degree-of-freedom motion equation of the rigid body can be obtained based on this decomposition. According to Newton’s law, the motion equation of the center of mass in the inertial coordinate system can be expressed as follows:
(6)md2xcmidt2=Fxi,(7)md2ycmidt2=Fyi,(8)md2zcmidt2=Fzi,where xcmi, ycmi, and zcmi represent the coordinates of the center of mass and Fxi, Fyi, and Fzi indicate the force components in the x, y, and zdirections, respectively. The superscript i indicates that the equation is established on the inertial coordinate system.The equation of rigid body rotation around the center of mass is expressed as follows:(9)Ixxbdωxbdt=Mxb+Iyyb−Izzbωybωzb,(10)Iyybdωybdt=Myb+Izzb−Ixxbωzbωxb,(11)Izzbdωzbdt=Mzb+Ixxb−Iyybωxbωyb,where subscript b represents the body frame, Mxb, Myb, and Mzb are the moments relative to the center of mass, ωxb, ωyb, and ωzb indicate the angular velocity, and Ixxb, Iyyb, and Izzb are the moments of inertia of the rigid body.The quaternion method is used to solve the equation of motion, which is expressed as follows:(12)P=e0e1e2e3.The transformation matrix composed of directional cosines can be expressed as(13)A=2e02+e12−0.5e1e2−e0e3e1e3+e0e2e1e2+e0e3e02+e22−0.5e2e3−e0e1e1e3−e0e2e2e3+e0e1e02+e32−0.5.The direction change caused by angular velocity can be calculated as follows:(14)dPdt=12Lωb,where
(15)L=−e1e0e3−e2−e2−e3e0e1−e3e2−e1e0.The force in Equation (7) can be decomposed into three parts, namely, the aerodynamic force, the gravity force, and the external force. The same decomposition can be performed for the moment in Equation (10). The aerodynamic force and moment are calculated based on CFD. The position and direction of the rigid body are updated by numerical integration of Equations (7), (10), and (14) using the predictor-corrector scheme [25]. The new position of the rigid body can be obtained as follows:
(16)rin+1=rc.m.in+1+An+1rb−rc.m.b,where rin+1 and rc.m.in+1 represent the coordinates of the rigid body and the center of mass in the inertial frame at the time step n+1, respectively. rb and rc.m.b indicate the coordinate of the rigid body and the center of mass in the body frame.


## 2.3. Structural Dynamics Solver

The structural dynamics system can be expressed as a system of second-order ordinary differential equations [26] as follows:
(17)Mu¨t+Cu̇t+Kut=Ft,where M, C, and K are the mass, damping, and stiffness matrices, respectively. u¨t, u̇t, and ut denote structural acceleration, velocity, and displacement vectors, respectively. F is the difference of the aerodynamic loads and inertial forces which are determined from the solution of equations of rigid body motion.In the present study, the Newmark method [27] is used to discretize the structural dynamics equation. Equation (17) can be written as follows:
(18)Mu¨n+1+Cu̇n+1+Kun+1=Fn+1.In addition to Equation (17), the Newmark method updates the displacement and velocity vectors as follows:
(19)u̇n+1=u̇n+1−δu¨n+δu¨n+1∆t,(20)un+1=un+u̇n∆t+12−αu¨n+αu¨n+1∆t2,whereδ and α are the Newmark integral parameters. The Newmark time integration algorithm can be determined based on the Newmark integration parameters.From Equations (18), (19), and (20), the one-step integration method can be expressed as follows:
(21)a0M+a1C+Kun+1=Fn+1+Ma0un+a2u̇n+a3u¨n+Ca1un+a4u̇n+a5u¨n,where
(22)a0=1α∆t2,a1=δα∆t,a2=1α∆t,a3=12α−1,a4=δα−1,a5=∆t2δα−2.un+1 can be calculated from Equation (21), and u̇n+1 and u¨n+1 can be obtained as follows:
(23)u̇n+1=a1un+1−un−a4u̇n−a5u¨n+1,u¨n+1=a0un+1−un−a2u̇n−a3u¨n+1.According to the stability analysis, the method is unconditionally stable when the Newmark parameters satisfy the following conditions:(24)δ≥12,α≥1412+δ2.


## 2.4. Hybrid Dynamic Grid Method

In the problem of missile release from the internal bay with consideration of aeroelastic effect, both the structure deformation and the overall motion relative to the aircraft are involved. The dynamic grid method can simultaneously treat both of them and is practical to simulate such complex process. A hybrid dynamic grid method is developed, in which the deformation grid deals with the aeroelastic deformation of structure and the overset unstructured grid handles the relative motion between separation bodies.


### 2.4.1. Deformation Grid Method

When considering the aeroelasticity effect, the structural deformation should be concerned in the process of release and separation, and this is simulated by the structure dynamics solver. In order to capture the flow after the structure deformation, the CFD near-wall grid should adapt to the deformation. Therefore, the deformable grid method is adopted in the current study, because it can deal with the deformable wall surface. The basic idea of the deformable mesh method [28] is to keep the mesh topology unchanged and absorb the movement of the solid wall boundary by adjusting the nodes position in the field grid. Although this method is only applicable for small-scale motion, it has the advantages of preserving the correlation information between meshes and avoiding interpolation at each step. The spring-tension-based deformation method is utilized to adjust the field grid node position.The spring-tension method models the connection between each two grid nodes by a spring. The displacement of a given boundary node produces a force, which is directly proportional to the displacement of all springs connected to that node [29]. According to Hooke’s law, the force imposed on the mesh node can be expressed as follows:
(25)F⟶i=∑jnikij∆x⟶j−∆x⟶i,where ∆x⟶j and ∆x⟶i represent the displacement of node i and its adjacent node j and ni denotes the number of nodes nearby node i. kij indicates the elasticity coefficient between nodes i and j, which is defined as follows:
(26)kij=kfacx⟶i−x⟶j,where kfac is the spring constant. In the equilibrium condition, the resultant force exerted by all springs connected to the central node is zero, leading to the following iterative equation:
(27)∆x⟶im+1=∑jnikij∆x⟶jm∑jnikij,where m represents the number of iteration steps.After running the structure dynamics solver, the wall boundary displacement of the CFD grid is obtained by interpolating the structural grid. Once the boundary deformation is known, Equation (27) can be solved by the Jacobian scanning for all the field grid nodes.


### 2.4.2. Dynamic Overset Unstructured Grid Method

In the missile release process, the missile moves relative to the aircraft and the grid must be adjusted to follow its movement. In the present study, a dynamic overset unstructured grid method is utilized to handle the moving boundary. This method was first introduced by Nakahashi [30], and it has the advantages of both unstructured grid and Chimera grid methods. The dynamic overset unstructured grid method is quite suitable to deal with unsteady flow past complex geometries and/or multiple bodies with relative movement. In the simulation of missile separation, subgrids are first generated around the missile and the carrier aircraft. The subgrids need to be generated only once, and the subgrids around the missile just move rigidly with its body during simulation of the whole separation process.In the present work, a wall distance-based overset unstructured grid assembly method is adopted. In order to obtain the hole region, this method applies a wall distance criterion to classify the grid nodes into three groups, i.e., active nodes, inactive nodes, and intergrid boundary nodes. A grid node is active when the distance between this node and its body surface is less than the distances to all other bodies. In contrast, a grid node is inactive if the distance to its own body is larger than the distances to the other bodies. The inactive nodes which are immediately next to the active ones can be further classified as the intergrid boundary nodes. The active nodes are used in the computations while the inactive nodes generate the hole region which is skipped in the computations. The intergrid boundary nodes, also referred to as interpolation nodes, are used to transfer the solution between the subgrids. The donor cells of the interpolation nodes can be efficiently found based on a reliable neighbor-to-neighbor search algorithm coupled with the alternating digital tree data structure [31, 32]. A general interpolation method [31] which is compatible with any element type is utilized for information transfer between the subgrids.In the unsteady simulation of missile separation, the overset unstructured grid assembly is applied with the new grid position before solving the flow equations. This is determined by superimposing the rigid and deformation motions.


## 2.4.1. Deformation Grid Method

When considering the aeroelasticity effect, the structural deformation should be concerned in the process of release and separation, and this is simulated by the structure dynamics solver. In order to capture the flow after the structure deformation, the CFD near-wall grid should adapt to the deformation. Therefore, the deformable grid method is adopted in the current study, because it can deal with the deformable wall surface. The basic idea of the deformable mesh method [28] is to keep the mesh topology unchanged and absorb the movement of the solid wall boundary by adjusting the nodes position in the field grid. Although this method is only applicable for small-scale motion, it has the advantages of preserving the correlation information between meshes and avoiding interpolation at each step. The spring-tension-based deformation method is utilized to adjust the field grid node position.The spring-tension method models the connection between each two grid nodes by a spring. The displacement of a given boundary node produces a force, which is directly proportional to the displacement of all springs connected to that node [29]. According to Hooke’s law, the force imposed on the mesh node can be expressed as follows:
(25)F⟶i=∑jnikij∆x⟶j−∆x⟶i,where ∆x⟶j and ∆x⟶i represent the displacement of node i and its adjacent node j and ni denotes the number of nodes nearby node i. kij indicates the elasticity coefficient between nodes i and j, which is defined as follows:
(26)kij=kfacx⟶i−x⟶j,where kfac is the spring constant. In the equilibrium condition, the resultant force exerted by all springs connected to the central node is zero, leading to the following iterative equation:
(27)∆x⟶im+1=∑jnikij∆x⟶jm∑jnikij,where m represents the number of iteration steps.After running the structure dynamics solver, the wall boundary displacement of the CFD grid is obtained by interpolating the structural grid. Once the boundary deformation is known, Equation (27) can be solved by the Jacobian scanning for all the field grid nodes.


## 2.4.2. Dynamic Overset Unstructured Grid Method

In the missile release process, the missile moves relative to the aircraft and the grid must be adjusted to follow its movement. In the present study, a dynamic overset unstructured grid method is utilized to handle the moving boundary. This method was first introduced by Nakahashi [30], and it has the advantages of both unstructured grid and Chimera grid methods. The dynamic overset unstructured grid method is quite suitable to deal with unsteady flow past complex geometries and/or multiple bodies with relative movement. In the simulation of missile separation, subgrids are first generated around the missile and the carrier aircraft. The subgrids need to be generated only once, and the subgrids around the missile just move rigidly with its body during simulation of the whole separation process.In the present work, a wall distance-based overset unstructured grid assembly method is adopted. In order to obtain the hole region, this method applies a wall distance criterion to classify the grid nodes into three groups, i.e., active nodes, inactive nodes, and intergrid boundary nodes. A grid node is active when the distance between this node and its body surface is less than the distances to all other bodies. In contrast, a grid node is inactive if the distance to its own body is larger than the distances to the other bodies. The inactive nodes which are immediately next to the active ones can be further classified as the intergrid boundary nodes. The active nodes are used in the computations while the inactive nodes generate the hole region which is skipped in the computations. The intergrid boundary nodes, also referred to as interpolation nodes, are used to transfer the solution between the subgrids. The donor cells of the interpolation nodes can be efficiently found based on a reliable neighbor-to-neighbor search algorithm coupled with the alternating digital tree data structure [31, 32]. A general interpolation method [31] which is compatible with any element type is utilized for information transfer between the subgrids.In the unsteady simulation of missile separation, the overset unstructured grid assembly is applied with the new grid position before solving the flow equations. This is determined by superimposing the rigid and deformation motions.


## 2.5. Simulation Process of Multibody Separation with Aeroelasticity

By providing the flow solver, the rigid body dynamics model, the structure solver, and the dynamic grid, numerical simulation of missile separation in the presence of aeroelastic effects is achievable by coupling CFD/RBD/CSD on the dynamic grid. Under the assumption of small deformation, the influence of deformation on the mass distribution and rotational inertia of the elastomer is ignored. This implies that the rigid body motion and structure deformation is a one-way coupling. Figure1 presents the flowchart of numerical simulation scheme for this problem. The simulation procedure can be described as follows:
(1)
Initialization of the flow field for unsteady calculations based on the steady flow solution(2)
Calculation of the aerodynamic force and moment at the current physical time step by integrating the distributed force on the wall surface(3)
Transfer of the aerodynamic force and moment to the RBD model and, at the same time, interpolation of the distributed aerodynamic force (imposed on the wall boundary of CFD grid) to the CSD grid(4)
Solving the RBD equation and execution of the structure dynamics solver to obtain the whole motion and structure deformation(5)
Evaluate whether the solution is converged: if the number of iterations is insufficient, update the CFD grid by the overset unstructured grid and deformation grid methods; otherwise, terminate the simulationFigure 1
Simulation procedure for the problem of multibody separation in the presence of aeroelastic effect.


## 3. Numerical Method Verification

In order to validate the accuracy and robustness of the coupled CFD/RBD/CSD solver on the hybrid dynamic grid, two test cases are considered in this section. One is the store separation from a wing, and the other is the flutter boundary simulation of AGARD445.6 wing.


### 3.1. Simulation of Store Separation from a Wing

The separation of a finned store from a delta wing is simulated to verify the accuracy of the coupled CFD/RBD method and the robustness of the dynamic overset unstructured grid method. This is a standard test case for multibody separation problem, including detailed model and wind tunnel test data [33–35]. The model is shown in Figure 2. A trimmed delta wing is considered with spanwise section of NACA64a010 airfoil, leading edge sweep angle of 45°, wing root chord length of 7.62 m, and half-span length of 6.6 m. The middle section of the missile model has a cylindrical geometry with diameter of 0.5 m. It has a length of about 3 m with pointed arch head and tail. Four tail wings with NACA0008 airfoil and 45° sweep angle are also included. In order to ensure a safe separation, the ejection forces are initially applied at the front and rear positions of the store centroid. When the missile movement is greater than a given threshold, the ejection force disappears. Table 1 presents the parameters of the store and the ejection force [36].Figure 2
Wing-store model geometry: (a) the overall 3D view of the computational configuration and the coordinate system; (b) the geometry of store and ejection model.
(a)(b)Table 1
Mass attributes and ejection force parameters of the missile.
Mass,m (kg)907.184Center of mass,XCG (m)1.417Moment of inertia,Ixx (kg·m3)27.116Moment of inertia,Iyy=Izz (kg·m3)488.094The position of ejection force on the leading edge (m)1.2375 lThe ejection force applied at the leading edge (N)10668.48The position of ejection force on the trailing edge (m)1.7465The ejection force applied at the trailing edge (N)42672.0Operating distance (m)0.1The flow condition of simulation is considered based on the wind tunnel test condition, such that the free stream Mach number is 1.2, the angle of attack is 0°, and the flight altitude is 11.6 km. The Spalart–Allmaras (S-A) one equation turbulence model is employed in the simulation. The dynamic overset unstructured grid includes the subgrids of the store and the wing. In unsteady condition, the wing subgrid remains static while the store subgrid moves based on the output of the RBD mode. The physical time step adopted for this test case is 0.01 s.Figure3 illustrates the missile position at four different times during the separation process. As can be seen, the missile moves downward, backward, and outward over time because of the gravitational effect and ejection and aerodynamic forces. Figure 4 compares the simulation and wind tunnel test results of the missile trajectory and attitude angle. As can be seen, the simulated displacement and velocity in the three directions are in good agreement with the experimental results. According to Figure 4(a), the slope of Y-component velocity sharply changes around t=0.06s, and the acceleration decreases significantly. This is due to the presence of an ejection force causing a considerable acceleration in the early stage of separation, while the ejection force disappears after 0.06 s leading to significant reduction of the acceleration. As can be seen from Figure 4(c), the general trend of all components of the simulated angular velocity is consistent with the experimental results, though some deviations can be observed. The noticeable difference between numerical and experimental values of the rolling angle speed is due to the small value of moment of inertia Ixx of the store. Even if there is a small difference in the rolling moment, it will cause a large change in the rolling angle speed. The simulated yaw angular velocity is in good agreement with the experimental results before 0.3 s, while their deviation gradually increases afterwards (with the same general trend). The pitch angular velocity increases linearly in the initial 0.06 s period due to the ejection force, whereas it decreases gradually after removing this force. According to Figure 4(d), in contrast to good agreement between numerical and experimental results of the yaw and roll angles, the deviation gradually increases for the pitch angle after 0.2 s. It can be observed from Figure 4(c) that after disappearance of the ejection force, there is always a certain difference between numerical and experimental values of the pitch angular velocity. This leads to a steady increase of the error obtained for the pitch angle.Figure 3
The store position at four different times during the separation process: (a) the side view; (b) the front view.
(a)(b)Figure 4
Temporal variation of various characteristics of the store during the separation process: (a) linear velocity; (b) linear displacement; (c) angular velocity; (d) angular displacement.
(a)(b)(c)(d)In the simulation of store separation, the dynamic overset unstructured grid method can capture the store motion automatically. Numerical results show that the CFD/RBD coupled method is robust and reliable in simulating the multibody separation problem.


### 3.2. Simulation of AGARD445.6 Wing Flutter

The AGARD445.6 wing has been frequently used to verify the reliability of the proposed method for aeroelastic problems. This model was used by Langley Research Center to study the wing flutter characteristics at transonic condition based on public wind tunnel test data [37]. The spanwise section of the wing is based on NACA65A004 airfoil with aspect ratio of 1.6440, root tip ratio of 0.659, and 1/4 chord sweep angle of 45°. The wing model geometry is depicted in Figure 5.Figure 5
The AGARD445.6 wing model geometry.Vibration of this wing is mainly characterized by the first four modes, namely, first-order bending, first-order torsion, second-order bending, and second-order torsion. To show the mode analysis of the structure model of this wing, the first four modes and their respective frequencies are presented in Figure6.Figure 6
Mode shape and the corresponding frequency of the AGARD445.6 wing at the first four principal modes: (a) first mode at the frequency of 9.73 Hz; (b) second mode at the frequency of 37.54 Hz; (c) third mode at the frequency of 49.62 Hz; (d) fourth mode at the frequency of 92.37 Hz.
(a)(b)(c)(d)The simulation Mach numbers are 0.499, 0.678, 0.901, 0.96, 1.072, and 1.141. The angle of attack of the incoming flow is 0°, and the critical flutter velocity at each Mach number is determined by changing the incoming flow pressure. The dimensionless flutter velocity is expressed as follows:(28)V∗=Vbsωαμ¯,where V represents the incoming flow velocity, bs is the half chord length of the wing root, and ωα denotes the first-order uncoupled torsional frequency.The simulation results of the present study are compared with both experimental and Liu’s computational results of Liu [38] and FUN3D [39], as illustrated in Figure 7. The transonic “pit” phenomenon is clearly captured. In the subsonic stage, relatively small differences can be observed between simulated and experimental values of the critical flutter velocity. As the fine CFD grid with 3.2 million elements is used in the unsteady flow simulation, in the transonic region, the simulated critical flutter velocity in the present paper has a better agreement with the experimental results. This implies reliability of the implemented numerical approach for aeroelasticity simulation based on the CFD/CSD coupled method.Figure 7
The flutter boundary of the AGARD445.6 wing.


## 3.1. Simulation of Store Separation from a Wing

The separation of a finned store from a delta wing is simulated to verify the accuracy of the coupled CFD/RBD method and the robustness of the dynamic overset unstructured grid method. This is a standard test case for multibody separation problem, including detailed model and wind tunnel test data [33–35]. The model is shown in Figure 2. A trimmed delta wing is considered with spanwise section of NACA64a010 airfoil, leading edge sweep angle of 45°, wing root chord length of 7.62 m, and half-span length of 6.6 m. The middle section of the missile model has a cylindrical geometry with diameter of 0.5 m. It has a length of about 3 m with pointed arch head and tail. Four tail wings with NACA0008 airfoil and 45° sweep angle are also included. In order to ensure a safe separation, the ejection forces are initially applied at the front and rear positions of the store centroid. When the missile movement is greater than a given threshold, the ejection force disappears. Table 1 presents the parameters of the store and the ejection force [36].Figure 2
Wing-store model geometry: (a) the overall 3D view of the computational configuration and the coordinate system; (b) the geometry of store and ejection model.
(a)(b)Table 1
Mass attributes and ejection force parameters of the missile.
Mass,m (kg)907.184Center of mass,XCG (m)1.417Moment of inertia,Ixx (kg·m3)27.116Moment of inertia,Iyy=Izz (kg·m3)488.094The position of ejection force on the leading edge (m)1.2375 lThe ejection force applied at the leading edge (N)10668.48The position of ejection force on the trailing edge (m)1.7465The ejection force applied at the trailing edge (N)42672.0Operating distance (m)0.1The flow condition of simulation is considered based on the wind tunnel test condition, such that the free stream Mach number is 1.2, the angle of attack is 0°, and the flight altitude is 11.6 km. The Spalart–Allmaras (S-A) one equation turbulence model is employed in the simulation. The dynamic overset unstructured grid includes the subgrids of the store and the wing. In unsteady condition, the wing subgrid remains static while the store subgrid moves based on the output of the RBD mode. The physical time step adopted for this test case is 0.01 s.Figure3 illustrates the missile position at four different times during the separation process. As can be seen, the missile moves downward, backward, and outward over time because of the gravitational effect and ejection and aerodynamic forces. Figure 4 compares the simulation and wind tunnel test results of the missile trajectory and attitude angle. As can be seen, the simulated displacement and velocity in the three directions are in good agreement with the experimental results. According to Figure 4(a), the slope of Y-component velocity sharply changes around t=0.06s, and the acceleration decreases significantly. This is due to the presence of an ejection force causing a considerable acceleration in the early stage of separation, while the ejection force disappears after 0.06 s leading to significant reduction of the acceleration. As can be seen from Figure 4(c), the general trend of all components of the simulated angular velocity is consistent with the experimental results, though some deviations can be observed. The noticeable difference between numerical and experimental values of the rolling angle speed is due to the small value of moment of inertia Ixx of the store. Even if there is a small difference in the rolling moment, it will cause a large change in the rolling angle speed. The simulated yaw angular velocity is in good agreement with the experimental results before 0.3 s, while their deviation gradually increases afterwards (with the same general trend). The pitch angular velocity increases linearly in the initial 0.06 s period due to the ejection force, whereas it decreases gradually after removing this force. According to Figure 4(d), in contrast to good agreement between numerical and experimental results of the yaw and roll angles, the deviation gradually increases for the pitch angle after 0.2 s. It can be observed from Figure 4(c) that after disappearance of the ejection force, there is always a certain difference between numerical and experimental values of the pitch angular velocity. This leads to a steady increase of the error obtained for the pitch angle.Figure 3
The store position at four different times during the separation process: (a) the side view; (b) the front view.
(a)(b)Figure 4
Temporal variation of various characteristics of the store during the separation process: (a) linear velocity; (b) linear displacement; (c) angular velocity; (d) angular displacement.
(a)(b)(c)(d)In the simulation of store separation, the dynamic overset unstructured grid method can capture the store motion automatically. Numerical results show that the CFD/RBD coupled method is robust and reliable in simulating the multibody separation problem.


## 3.2. Simulation of AGARD445.6 Wing Flutter

The AGARD445.6 wing has been frequently used to verify the reliability of the proposed method for aeroelastic problems. This model was used by Langley Research Center to study the wing flutter characteristics at transonic condition based on public wind tunnel test data [37]. The spanwise section of the wing is based on NACA65A004 airfoil with aspect ratio of 1.6440, root tip ratio of 0.659, and 1/4 chord sweep angle of 45°. The wing model geometry is depicted in Figure 5.Figure 5
The AGARD445.6 wing model geometry.Vibration of this wing is mainly characterized by the first four modes, namely, first-order bending, first-order torsion, second-order bending, and second-order torsion. To show the mode analysis of the structure model of this wing, the first four modes and their respective frequencies are presented in Figure6.Figure 6
Mode shape and the corresponding frequency of the AGARD445.6 wing at the first four principal modes: (a) first mode at the frequency of 9.73 Hz; (b) second mode at the frequency of 37.54 Hz; (c) third mode at the frequency of 49.62 Hz; (d) fourth mode at the frequency of 92.37 Hz.
(a)(b)(c)(d)The simulation Mach numbers are 0.499, 0.678, 0.901, 0.96, 1.072, and 1.141. The angle of attack of the incoming flow is 0°, and the critical flutter velocity at each Mach number is determined by changing the incoming flow pressure. The dimensionless flutter velocity is expressed as follows:(28)V∗=Vbsωαμ¯,where V represents the incoming flow velocity, bs is the half chord length of the wing root, and ωα denotes the first-order uncoupled torsional frequency.The simulation results of the present study are compared with both experimental and Liu’s computational results of Liu [38] and FUN3D [39], as illustrated in Figure 7. The transonic “pit” phenomenon is clearly captured. In the subsonic stage, relatively small differences can be observed between simulated and experimental values of the critical flutter velocity. As the fine CFD grid with 3.2 million elements is used in the unsteady flow simulation, in the transonic region, the simulated critical flutter velocity in the present paper has a better agreement with the experimental results. This implies reliability of the implemented numerical approach for aeroelasticity simulation based on the CFD/CSD coupled method.Figure 7
The flutter boundary of the AGARD445.6 wing.


## 4. Results and Discussion

Because of the large slenderness ratio and low structural stiffness of the air-to-air missile, aeroelasticity plays an important role in the separation process and its characteristics. In the severe circumstances, it can even threaten the safety of the carrier aircraft. This section examines the impact of aeroelasticity in multibody separation process based on the CFD/RBD/CSD coupled method.


### 4.1. Computational Model

Figure8 illustrates the aircraft carrier and missile model used in this paper. Two embedded bomb pods are mounted in the middle of the aircraft carrier, which are disconnected from each other internally. The hatch opening angle is 150°, and the size of the internal bomb compartment is 4.0m×2.0m×0.5m. The air-to-air missile has total length of 2.87 m, diameter of 0.137 m, wingspan of 0.63 m, and launch mass of 87 kg. The missile is situated in the center of the internal bay. The front rudder and the tail wing of the missile have X-shaped configurations.Figure 8
The wall surface grid topology in the study of aeroelastic effect in the missile separation process: (a) the carrier aircraft; (b) the missile head section; (c) the missile tail section.
(a)(b)(c)In order to utilize the overset unstructured grid, the subgrids around the carrier aircraft and the missile are initially generated. Before the unsteady simulations of the missile separation, the convergence of grid has been studied by steady computation. We studied the influence of the number of grid size on the calculation results. In addition, to improve the grid resolution in dynamic computation, the subgrid of carrier aircraft is refined in the internal bay and the region where the missile may pass and the local mesh spacing is set similar to that of the outer boundary elements in the subgrid of the missile. The carrier subgrid with 5 million elements and the missile grid pertaining to about 4 million elements are utilized in the following unsteady simulation tests. Figure8 demonstrates the surface grids of the carrier and missile.In order to study the aeroelastic effect of the missile on the separation trajectory, two structural models are considered, namely, the rigid and the elastic models. Both of these models have a density of 2150 kg/m3 and Poisson’s ratio of 0.3. Young’s modulus of the rigid and elastic models is 2.0×1011 Pa and 4.0×109 Pa, respectively. The ejection force is applied in the form of an initial velocity which has the magnitude of 8.2 m/s and a downward direction (towards the negative y-axis) in the same direction of the gravitational acceleration. The first four modes of the elastic model and their respective frequencies are presented in Figure 9 to show the main vibration characteristic.Figure 9
Mode shapes of the elastic model at the first four vibration modes: (a) the first mode at the frequency of 19.01 Hz; (b) the second mode at the frequency of 19.03 Hz; (c) the third mode at the frequency of 54.86 Hz; (d) the fourth mode at the frequency of 54.91 Hz.
(a)(b)(c)(d)


### 4.2. Aeroelasticity Effect on the Missile Separation for Different Mach Numbers

The incoming Mach numbers of 0.8, 1.2, and 1.5 are considered to study the aeroelastic effect on the missile separation trajectory. The time step convergence is studied on the verified grid. It is found that when the time step is not larger than 0.5 ms, the unsteady separation results of the missile are consistent. Figure10 illustrates the contour plots of pressure field for different Mach numbers at various times. It is evident that the flow field structure is very different for the Mach numbers of 0.8 and 1.2. In the supersonic condition, an expansion wave existed near the missile internal bay and an oblique shock wave is generated by the forebody far from the internal bay. Therefore, shortly after releasing from the internal bay, the missile passes through the expansion zone which can generate a strong detached shock wave ahead the missile, and then, it enters into the oblique shock wave. However, at the Mach number of 0.8, although aerodynamic interference is available in the missile separation process, no obvious expansion wave and shock wave can be observed under the internal bay. Therefore, the existence of a stronger aerodynamic interference in the supersonic flight condition intensifies the aeroelasticity impact on the missile separation. This interference can be reflected in the history of separation.Figure 10
Contour plots of pressure field for the Mach numbers of 0.8 (left frames, I) and 1.2 (right frames, II) at four different times: (a) 0.0 s; (b) 0.05 s; (c) 0.1 s; (d) 0.2 s.Figure11 presents the temporal variation of displacement and attitude angle of the rigid and elastic models for different Mach numbers. As can be seen, the X-component displacement difference between rigid and elastic models is very small for various Mach numbers, while the Y-component displacement difference grows gradually by increase of the Mach number. When the incoming Mach number is 1.5, the displacement trajectories of the two models have an obvious change of trend in the vertical direction. At Mach numbers of 0.8 and 1.2, the Z-component displacement difference between rigid and elastic models is small, whereas this difference is relatively large at Mach number of 1.5.Figure 11
Temporal variation of the linear and angular displacements for the rigid and elastic models at different Mach numbers: (a) X-component displacement; (b) Y-component (vertical) displacement; (c) Z-component (lateral) displacement; (d) pitch angle displacement; (e) yaw angle displacement; (f) roll angle displacement.
(a)(b)(c)(d)(e)(f)Parts (d) to (f) of Figure11 indicate that there are considerable differences in the attitude angles of rigid and elastic models. Temporal variation of the pitch angle at Ma=0.8 is different from that experienced at Mach numbers of 1.2 and 1.5. At the Mach number of 0.8, the pitch angles of both models gradually increase, while the difference between the two models begins to increase after 0.1 s. More specifically, the pitch angle of the elastic model increases by 160% compared to that of the rigid model at 0.2 s. At Mach numbers of 1.2 and 1.5, large amplitude oscillations can be observed in the pitch angle of both models, though the oscillation amplitude and frequency are lower in the case of Ma=1.2. In addition, the rigid model has a smaller oscillation frequency than the flexible model. Figure 11(e) indicates a slight difference of the yaw angle between rigid and elastic models at the initial stage for all three Mach numbers. This difference becomes apparent after about 0.1 s when the missile enters the mainstream from the internal bay. During the simulation, the yaw angle displacement of the flexible model is larger than that of the rigid model as the yaw angle oscillation period of the flexible model is improved. Figure 11(f) shows that aeroelasticity has the most significant effect on the temporal variation of the roll angle, since difference of the roll angle between the rigid and elastic models is quite noticeable. Although the roll angle of the rigid model changes very slightly and it stabilizes quickly, the roll angle of the elastic model oscillates with an elevated period and amplitude. Furthermore, in the case of the flexible model, the oscillation frequency of the roll angle displacement increases as the Mach number rises.In order to further analyze the impact of aeroelasticity in the missile separation process, Figure12 compares structural deformations of the rigid and elastic models at different times for the Mach number of 0.8. Based on temporal variation of separation characteristics, it was previously discussed that the difference between these characteristics for the two models is very small at the initial stage, while downward movement of the missile gradually increases the difference of angular displacement between the rigid and elastic models. In the process of missile separation, with gradual increase of the pitch angle, the elastic deformation produced by the elastic model enhances the head down moment. This further increases the pitch angle velocity and consequently leads to higher deformation of the missile head. This coupled process makes the pitch angle of the elastic model much greater than that of the rigid model. The same trend can be observed for the coupling of yaw angle displacement and structural deformation.Figure 12
Comparison of the rigid model deformation (blue) with the elastic model deformation (red) at three different times: (a) 0.035 s; (b) 0.055 s; (c) 0.075 s. The left (I) and right (II) plots represent the side and the top views, respectively.


### 4.3. Aeroelasticity Effect on Missile Separation at Different Flight Altitudes

According to the results of Section4.2, the most dominant aeroelastic effect on the missile separation was at the Mach number of 1.5. Therefore, this Mach number is considered in the present section to study the aeroelastic effects with the change of flight altitude. Three flight altitudes of 0 km, 5 km, and 8.8 km are adopted. Temporal variations of the linear displacement and attitude angle displacement for the rigid and elastic models are shown in Figure 13.Figure 13
Temporal variation of the linear and angular displacements for the rigid and elastic models for the Mach number of 1.5 at different flight altitudes: (a) X-component displacement; (b) vertical displacement; (c) lateral displacement; (d) pitch angle displacement; (e) yaw angle displacement; (f) roll angle displacement.
(a)(b)(c)(d)(e)(f)Similar to the results obtained at different Mach numbers, Figures13(a)–13(c) indicate that except for the lateral displacement, the impact of aeroelasticity on the displacement of the center of mass changes slightly with flight height. Although the difference of the lateral displacement between the rigid and elastic models at the last stage of separation seems to be considerable, its magnitude is small.Figures13(d)–13(f) reveal that there are obvious differences in the attitude angle displacement between rigid and flexible models. At the initial separation stage, as shown in Figures 13(d) and 13(e), aeroelasticity has a minimal effect on the pitch and yaw angle displacements, while after 0.11 s, the missile enters the mainstream causing enhancement of the aeroelastic effects. Furthermore, the flight altitude at a fixed Mach number can contribute to the aeroelastic effects. Such effect is more obvious at a lower flight altitude that corresponds to a higher dynamic pressure. In addition, during the missile separation process, aeroelasticity increases the oscillation period of angular displacement. The same trend is true for the roll angle displacement, and aeroelasticity has an even greater influence on the roll angle. As shown in Figure 13(f), the roll angle of the rigid model changes negligibly and it stabilizes quickly. However, the roll angle of the flexible model varies periodically with amplitude which is similar to the pitch angle magnitude. This can easily lead to the vertical and horizontal coupled movements of the missile, which threaten the separation security.


## 4.1. Computational Model

Figure8 illustrates the aircraft carrier and missile model used in this paper. Two embedded bomb pods are mounted in the middle of the aircraft carrier, which are disconnected from each other internally. The hatch opening angle is 150°, and the size of the internal bomb compartment is 4.0m×2.0m×0.5m. The air-to-air missile has total length of 2.87 m, diameter of 0.137 m, wingspan of 0.63 m, and launch mass of 87 kg. The missile is situated in the center of the internal bay. The front rudder and the tail wing of the missile have X-shaped configurations.Figure 8
The wall surface grid topology in the study of aeroelastic effect in the missile separation process: (a) the carrier aircraft; (b) the missile head section; (c) the missile tail section.
(a)(b)(c)In order to utilize the overset unstructured grid, the subgrids around the carrier aircraft and the missile are initially generated. Before the unsteady simulations of the missile separation, the convergence of grid has been studied by steady computation. We studied the influence of the number of grid size on the calculation results. In addition, to improve the grid resolution in dynamic computation, the subgrid of carrier aircraft is refined in the internal bay and the region where the missile may pass and the local mesh spacing is set similar to that of the outer boundary elements in the subgrid of the missile. The carrier subgrid with 5 million elements and the missile grid pertaining to about 4 million elements are utilized in the following unsteady simulation tests. Figure8 demonstrates the surface grids of the carrier and missile.In order to study the aeroelastic effect of the missile on the separation trajectory, two structural models are considered, namely, the rigid and the elastic models. Both of these models have a density of 2150 kg/m3 and Poisson’s ratio of 0.3. Young’s modulus of the rigid and elastic models is 2.0×1011 Pa and 4.0×109 Pa, respectively. The ejection force is applied in the form of an initial velocity which has the magnitude of 8.2 m/s and a downward direction (towards the negative y-axis) in the same direction of the gravitational acceleration. The first four modes of the elastic model and their respective frequencies are presented in Figure 9 to show the main vibration characteristic.Figure 9
Mode shapes of the elastic model at the first four vibration modes: (a) the first mode at the frequency of 19.01 Hz; (b) the second mode at the frequency of 19.03 Hz; (c) the third mode at the frequency of 54.86 Hz; (d) the fourth mode at the frequency of 54.91 Hz.
(a)(b)(c)(d)


## 4.2. Aeroelasticity Effect on the Missile Separation for Different Mach Numbers

The incoming Mach numbers of 0.8, 1.2, and 1.5 are considered to study the aeroelastic effect on the missile separation trajectory. The time step convergence is studied on the verified grid. It is found that when the time step is not larger than 0.5 ms, the unsteady separation results of the missile are consistent. Figure10 illustrates the contour plots of pressure field for different Mach numbers at various times. It is evident that the flow field structure is very different for the Mach numbers of 0.8 and 1.2. In the supersonic condition, an expansion wave existed near the missile internal bay and an oblique shock wave is generated by the forebody far from the internal bay. Therefore, shortly after releasing from the internal bay, the missile passes through the expansion zone which can generate a strong detached shock wave ahead the missile, and then, it enters into the oblique shock wave. However, at the Mach number of 0.8, although aerodynamic interference is available in the missile separation process, no obvious expansion wave and shock wave can be observed under the internal bay. Therefore, the existence of a stronger aerodynamic interference in the supersonic flight condition intensifies the aeroelasticity impact on the missile separation. This interference can be reflected in the history of separation.Figure 10
Contour plots of pressure field for the Mach numbers of 0.8 (left frames, I) and 1.2 (right frames, II) at four different times: (a) 0.0 s; (b) 0.05 s; (c) 0.1 s; (d) 0.2 s.Figure11 presents the temporal variation of displacement and attitude angle of the rigid and elastic models for different Mach numbers. As can be seen, the X-component displacement difference between rigid and elastic models is very small for various Mach numbers, while the Y-component displacement difference grows gradually by increase of the Mach number. When the incoming Mach number is 1.5, the displacement trajectories of the two models have an obvious change of trend in the vertical direction. At Mach numbers of 0.8 and 1.2, the Z-component displacement difference between rigid and elastic models is small, whereas this difference is relatively large at Mach number of 1.5.Figure 11
Temporal variation of the linear and angular displacements for the rigid and elastic models at different Mach numbers: (a) X-component displacement; (b) Y-component (vertical) displacement; (c) Z-component (lateral) displacement; (d) pitch angle displacement; (e) yaw angle displacement; (f) roll angle displacement.
(a)(b)(c)(d)(e)(f)Parts (d) to (f) of Figure11 indicate that there are considerable differences in the attitude angles of rigid and elastic models. Temporal variation of the pitch angle at Ma=0.8 is different from that experienced at Mach numbers of 1.2 and 1.5. At the Mach number of 0.8, the pitch angles of both models gradually increase, while the difference between the two models begins to increase after 0.1 s. More specifically, the pitch angle of the elastic model increases by 160% compared to that of the rigid model at 0.2 s. At Mach numbers of 1.2 and 1.5, large amplitude oscillations can be observed in the pitch angle of both models, though the oscillation amplitude and frequency are lower in the case of Ma=1.2. In addition, the rigid model has a smaller oscillation frequency than the flexible model. Figure 11(e) indicates a slight difference of the yaw angle between rigid and elastic models at the initial stage for all three Mach numbers. This difference becomes apparent after about 0.1 s when the missile enters the mainstream from the internal bay. During the simulation, the yaw angle displacement of the flexible model is larger than that of the rigid model as the yaw angle oscillation period of the flexible model is improved. Figure 11(f) shows that aeroelasticity has the most significant effect on the temporal variation of the roll angle, since difference of the roll angle between the rigid and elastic models is quite noticeable. Although the roll angle of the rigid model changes very slightly and it stabilizes quickly, the roll angle of the elastic model oscillates with an elevated period and amplitude. Furthermore, in the case of the flexible model, the oscillation frequency of the roll angle displacement increases as the Mach number rises.In order to further analyze the impact of aeroelasticity in the missile separation process, Figure12 compares structural deformations of the rigid and elastic models at different times for the Mach number of 0.8. Based on temporal variation of separation characteristics, it was previously discussed that the difference between these characteristics for the two models is very small at the initial stage, while downward movement of the missile gradually increases the difference of angular displacement between the rigid and elastic models. In the process of missile separation, with gradual increase of the pitch angle, the elastic deformation produced by the elastic model enhances the head down moment. This further increases the pitch angle velocity and consequently leads to higher deformation of the missile head. This coupled process makes the pitch angle of the elastic model much greater than that of the rigid model. The same trend can be observed for the coupling of yaw angle displacement and structural deformation.Figure 12
Comparison of the rigid model deformation (blue) with the elastic model deformation (red) at three different times: (a) 0.035 s; (b) 0.055 s; (c) 0.075 s. The left (I) and right (II) plots represent the side and the top views, respectively.


## 4.3. Aeroelasticity Effect on Missile Separation at Different Flight Altitudes

According to the results of Section4.2, the most dominant aeroelastic effect on the missile separation was at the Mach number of 1.5. Therefore, this Mach number is considered in the present section to study the aeroelastic effects with the change of flight altitude. Three flight altitudes of 0 km, 5 km, and 8.8 km are adopted. Temporal variations of the linear displacement and attitude angle displacement for the rigid and elastic models are shown in Figure 13.Figure 13
Temporal variation of the linear and angular displacements for the rigid and elastic models for the Mach number of 1.5 at different flight altitudes: (a) X-component displacement; (b) vertical displacement; (c) lateral displacement; (d) pitch angle displacement; (e) yaw angle displacement; (f) roll angle displacement.
(a)(b)(c)(d)(e)(f)Similar to the results obtained at different Mach numbers, Figures13(a)–13(c) indicate that except for the lateral displacement, the impact of aeroelasticity on the displacement of the center of mass changes slightly with flight height. Although the difference of the lateral displacement between the rigid and elastic models at the last stage of separation seems to be considerable, its magnitude is small.Figures13(d)–13(f) reveal that there are obvious differences in the attitude angle displacement between rigid and flexible models. At the initial separation stage, as shown in Figures 13(d) and 13(e), aeroelasticity has a minimal effect on the pitch and yaw angle displacements, while after 0.11 s, the missile enters the mainstream causing enhancement of the aeroelastic effects. Furthermore, the flight altitude at a fixed Mach number can contribute to the aeroelastic effects. Such effect is more obvious at a lower flight altitude that corresponds to a higher dynamic pressure. In addition, during the missile separation process, aeroelasticity increases the oscillation period of angular displacement. The same trend is true for the roll angle displacement, and aeroelasticity has an even greater influence on the roll angle. As shown in Figure 13(f), the roll angle of the rigid model changes negligibly and it stabilizes quickly. However, the roll angle of the flexible model varies periodically with amplitude which is similar to the pitch angle magnitude. This can easily lead to the vertical and horizontal coupled movements of the missile, which threaten the separation security.


## 5. Conclusion

The aeroelasticity effect in the multibody separation process is analyzed in the present study by solving the flow equations coupled with the rigid body dynamics model and the structural dynamics solver. A hybrid dynamic grid method is applied based on the overset unstructured grid and deformation grid. After verification of the robustness and reliability of the numerical method, the aeroelasticity effect on the missile separation from an internal bay is investigated at different states. Based on a given initial velocity, numerical results indicate that aeroelasticity has a little influence on the linear displacement of the missile, while it considerably affects temporal variation of the angular displacement during the separation process. It can increase the oscillation amplitude and period of the angular displacement. The flight Mach number and altitude have also important contributions to the aeroelasticity effect on the missile separation from the internal bay. A greater Mach number or a lower flight altitude can increase the impact of aeroelasticity on the angular displacement. Dependence of the missile attitude angle to the aeroelasticity effects is because of production of bending deformation of the missile in the separation process. This generates an additional moment which can further affect the angular displacement.



---

*Source: 9875622-2023-02-16.xml*