Effect of elastic deformation on flight dynamics of projectiles with large slenderness ratio

Aerospace Science and Technology 71 (2017) 347–359

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Effect of elastic deformation on flight dynamics of projectiles with large slenderness ratio Hua RuHao ∗ , Ye ZhengYin, Wu Jie School of Aeronautics, Northwestern Polytechnical University, Xi’an, 710072, China

a r t i c l e

i n f o

Article history: Received 12 September 2016 Received in revised form 23 July 2017 Accepted 18 September 2017 Available online 21 September 2017 Keywords: Elastic deformation Flight dynamics Six degree of freedom motion Large slenderness ratio

a b s t r a c t The elastic deformation of modern projectiles with large slenderness ratio cannot be ignored with the increasing of flight speed and maneuverability. Unsteady Reynolds-averaged Navier–Stokes (URANS) Equations are solved through CFD technique in this paper. Based on the frame of unstructured mesh, techniques of rigid-motion mesh and inverse-distance-weighted (IDW) morphing mesh are adopted to treat the rigid motion caused by flight dynamics and flexible structure deformation due to aeroelasticity, respectively. Moreover, the six degree of freedom (SDOF) dynamic equations and static aeroelastic equation are solved through the aerodynamic coupling. Numerical results of both free flight case and aeroelastic case calculated by the in-house code agree well with the experimental data, validating the numerical method. A projectile model with X–X configuration is constructed to investigate the effect of elastic deformation on the flight dynamics. Comparison results show that the longitudinal oscillation is more affected by the elastic deformation than the centroid motion, and the oscillation cycle of the orientation angle increases. Furthermore, the trajectories of rigid models with various centroid locations are simulated, illustrating that the elastic deformation could move the aerodynamic center forward and weaken the margin of the static stability margin. In the end, detailed analysis and comparison of the pressure distribution indicates the mechanism by which the elastic deformation leads to the movement of the aerodynamic center and changes the flight dynamic characteristics of the flexible projectile. © 2017 Elsevier Masson SAS. All rights reserved.

1. Introduction High thrust–weight ratio and large length-diameter (slenderness) ratio are needed for the design of projectiles to reduce the aerodynamic drag, launch heavier payloads and increase the range. As the slenderness ratio increases, the structure of the projectile becomes more and more flexible, and the natural frequencies of longitudinal and transverse vibration decrease. Moreover, in order to improve the maneuverability and agility, stability margin of modern aircrafts is limited or even unstable configuration is adopted. These requirements lead to more attention on the research of dynamic response and vibration characteristics of flexible flight vehicles. The early missiles and projectiles with smaller slenderness ratio and denser structure can be reasonably regarded as rigid models for dynamic analysis. Unfortunately, the elastic deformation of the new practical missile cannot be ignored, which might be beyond the range of common tolerance introduced in the process of manufacture [1–5]. Besides, fairings of large carrier rockets are more prone to aeroelastic problems such as breathing

*

Corresponding author. E-mail address: [email protected] (R.H. Hua).

https://doi.org/10.1016/j.ast.2017.09.029 1270-9638/© 2017 Elsevier Masson SAS. All rights reserved.

deformation due to thin-wall structures, which makes it more difficult to predict the trajectories of the separated components suffering from serious aerodynamic interference [6]. Since the structural rigidity is lower, it is easier for the air-to-air missiles and rocket missiles with large slenderness ratio to deform during the process of maneuvering rapidly, which not only has adverse effect on the trajectory, but also causes uncertain threat to the carrier. Accordingly, it is of great importance to regard the objects as flexible models to ensure the flight safety and firing accuracy. Since flight dynamics of flexible flight vehicles has been involved in severely mutual interaction among aerodynamics, flight mechanics, elastic forces and control systems [7–15], previous studies mainly handled the problem with simplified models which had various assumptions. Reis [7] addressed the problem of aeroelastic bending of a spinning rocket using a simple two-rigid-body model and pointed out that the transverse bending of the free– flight rocket could lead to obvious lateral angular velocity. Womack [8] developed a new “integrated method” based on normal modes to model the flexibility of sounding rockets and investigated the coupling of the vibratory and roll motions. Based on the Euler– Bernoulli beam model and Hamilton principle, Wang [9] undertook the trajectory simulation of a rocket and found that the range

348

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

of the rocket could be increased while the flight level became lower compared to the results of rigid-body models. Tuzcu [10] investigated the stability of flexible aircraft with various dynamic models and assumptions. Oliveira [11] treated the sounding rocket as a mass-varying Euler–Bernoulli free–free beam and presented a mathematical model to calculate the trajectory, flight dynamic parameters and the elastic deformations. Li [12] addressed the rigid-elastic coupled equations of motion for trajectory simulation of a spinning symmetric vehicle with non-uniform free–free Euler– Bernoulli beam model and quasi-steady theory. Two typical related studies provided by Waszak [13] and Buttrill [14], developed the equations of motion for elastic aircraft to meet the compromise between the accuracy and the simplification in dynamic modeling. The emphasis of the first study, which could be referred to as the Waszak’s study, was the assembly of a mathematical model which integrated “rigid-body” and “elastic” degrees of freedom with particular emphasis on the assumptions made at the various stages in the development and on obtaining a set of equations that constituted a literal model. The second study, which could be referred to as the Buttrill’s study, focused on including effects of nonlinear inertial coupling between rigid-body angular rates and structural deformations and rates which were usually ignored in conventional aircraft modeling. Although a great number of researches have been expanded to investigate the effects of flexibility on the stability of projectiles and launch vehicles, the aerodynamic loads which dominate the trajectory were normally predicted with engineering methods such as lifting-line method and perturbation method. Furthermore, most researches involved with the coupling of elastic deformation and flight dynamic characteristics focus on the coupling phenomena itself, while the investigation for the mechanism by which the elastic deformation affects flight dynamics is very limited. Since the accuracy of the engineering method is limited and will lead to more cumulative errors, it is necessary to develop more exact methods to predict the characteristics of the flight dynamics for the flexible projectiles. With the rapid development of numerical methods and computer technology, computational fluid dynamics (CFD) has been widely used to investigate the aerodynamic performance of flight vehicles [15–21]. Schütte [17] at DLR simulated the unsteady aerodynamics of a free–flight aeroelastic combat aircraft by using the flow-solver TAU which coupled aerodynamics, flight mechanics, and aeroelastic computations. The CREATE-AV project established by the DoD (the United States Department of Defense) aimed to develop, deploy and support a set of multidisciplinary, physics-based simulation software products for the engineering workforces to support air vehicle acquisition programs, and the project can handle the problems related to flexible aircraft with SDOF motion and control surfaces deflection [18–21]. Abbas [22] developed a numerical method for the deformation calculation of a rocket and its influence on aerodynamic parameters was studied based on static analyzes using CFD/CSD workflow. In this paper, based on the frame of unstructured hybrid mesh, the rigid-motion mesh technique and the mesh deformation method are applied to treat the large displacement due to SDOF motion and the elastic deformation due to aeroelasticity, respectively. URANS equations are solved to obtain the instant aerodynamic forces/moments at different time. Generalized flight dynamic equation composed of SDOF motion and static aeroelastic equation are constructed and solved to investigate the effect of elastic deformation on the flight dynamic characteristics.

Fig. 1. Flow chart of the solving methodology.

grid-treating module, CFD solver module and fluid–solid interaction module. The grid-treating module is capable of adapting elastic deformation of the meshes caused by aeroelasticity and adjusting the new computational domains according to the large displacement due to SDOF motion. The CFD solver module is used to obtain the numerical solution of the fluid governing equations at each time step. With the solution obtained in CFD solver module, the aerodynamic forces/moments on the projectile are computed by integrating the pressure and skin friction over the body. Meanwhile, the generalized forces needed for the solution of the aeroelastic equation can be obtained. The rigid motion and the elastic deformation can be computed by solving generalized flight dynamic equations for flexible aircrafts in the fluid–solid interaction module. The methodology in this paper is a loosely coupled system similar with Waszak’s study [13], and the effect of the aeroelasticity on flight dynamics is embodied in the new location of the mesh nodes. The combined mesh node displacement is the summation of the rigid displacement and the elastic deformation from the fluid–solid interaction module. Correspondingly, the flow field is composed of synthetic effect of both flight dynamics and aeroelasticity, and the focus of this paper is mainly on the aerodynamic coupling between rigid-body motion and structural deformation, while the effects of nonlinear inertial coupling between rigid-body angular rates and structural deformations and rates are not included in the flexible aircraft modeling. Various components of this framework are described as follows.

2. Overview of the framework for flight dynamics coupled with aeroelasticity

2.1. Dynamic mesh technique

The frame of the solving methodology presented in this paper is shown in Fig. 1, which mainly consists of three modules:

With the rapid progress of computer science and technology, more detailed research of the flow field around complicated con-

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

figurations can be treated with CFD technique. It is widely acknowledged that unstructured meshes are more convenient and flexible than structured grids to handle complex models. Hybrid unstructured meshes are selected for the discretization of the flow field in this paper. The rigid-motion mesh technique is adopted to treat the large rigid-body displacement due to the SDOF motion, while the flexible structure deformation caused by aeroelasticity is dealt with by the mesh morphing technique. 2.1.1. Rigid-motion mesh technique The basic idea of the rigid-motion mesh technique is to update the location of all nodes of the computational mesh according to the translational displacement of the centroid of the flight vehicle and the rotational displacement around the centroid. Here, the displacement of the centroid and the orientation angles could be obtained through solving the flight dynamic equation in each time step. The distinct advantage of the rigid-motion mesh technique is that the computational mesh can be obtained directly without being reconstructed during the numerical simulation process, which not only reduces the computational cost dramatically, but also keeps the topological relationship and the quality of the computational mesh. 

Here, r 0 is defined as the coordinate of the node of the computational mesh in the inertial system at the initial time, and the 

new location r of the corresponding node during the simulation process can be expressed as 







r =r 0 + r cg + r rot

(1)



Where  r cg is the displacement of the centroid of the flight vehicle, and it has no distinction for all the nodes of the computational 

mesh.  r rot is the displacement due to the rotation around the centroid, which can be obtained through the following formula 





⎡

T B−I =

ψ T B−I T B−I T B−I θ

1 0 = ⎣ 0 cos φ 0 sin φ

⎤⎡

⎡

cos θ ×⎣ 0 − sin θ

⎤

0 − sin φ ⎦ cos φ

0 sin θ cos ψ 1 0 ⎦ ⎣ sin ψ 0 0 cos θ

⎤ − sin ψ 0 cos ψ 0⎦ 0

− − →

wi =

b∈Ξ (B) [(



−−→ xib )a w b ]

b∈Ξ (B) (



∂ ∂t

˚

¨

Q dV + Ω

xib )a

(4)

¨

F V (Q) · n dS

F(Q, Vgrid ) · n dS = ∂Ω

(5)

∂Ω

where Q = [ρ , ρ u , ρ v , ρ w , e ] T is the vector of conservative variables; ρ is the density; u , v , w are the velocity components in the x, y and z directions respectively; e is the total energy. F(Q , Vgrid ) and F V (Q ) are the inviscid flux vector and the viscous flux vector, respectively; Vgrid is the velocity vector of a control volume due to SDOF motion and elastic deformation. Ω represents a control volume with the boundary expressed as ∂Ω and n is the normal vector of the boundary face of a control volume. A cell-centered finite volume method is employed based on unstructured hybrid meshes which are composed of hexahedrons, prisms, tetrahedrons and pyramids. Discretize Eq. (5) at a control volume and it can be expressed as:

dt

  = − R(Q n , Vgrid ) − R V (Q n )

(6)

where V n is the volume of the control volume signed as n, R(Q n , Vgrid ) and R V (Q n ) are the residual values of the inviscid flux and the viscous flux, respectively. The time derivative term dQ n V n represents the unsteady characteristic of the control voldt ume n. While the pseudo-time variable τ is introduced, Eq. (6) becomes:

(3)





Written in the time-dependent integral form of a moving and/or deforming control volume Ω with a surface element dS, the Navier–Stokes equations are as follows:

dτ

where



=−

3Q nk+1 V nk+1 − 4Q nk V nk + Q nk−1 V nk−1

+1

− R¯ Q nk+1 , Vkgrid

1

2.1.2. Mesh deformation method Elastic deformation caused by external forces exists besides the large displacement due to rigid-body motion during the simulation of the free flight process of the flexible projectiles. Each nodal position on the boundaries should be determined by the summation of the position of the rigid motion and the elastic deformation through solving the generalized flight dynamic equations. In this paper an inverse distance-weighted (IDW) method [23] is employed for the establishment of a special weighted function to diffuse the displacements of wall surfaces to the interior mesh nodes of the computational zone. The interpolation function is expressed as:



2.2. Governing equations of fluid flow and numerical method

d(Q nk+1 V nk+1 )

where φ , θ and ψ are the rotation angle, the pitching angle and the yaw angle, respectively.





solution of the aeroelastic equation and  = xib = r 0,b − r 0,i represents the vector between node b and node i. Exponent a decides the propagation speed of the deformation. In order to ensure the quality of the small-scale meshes in the boundary layer for the simulation of viscous flow, a can be set within the range of 3–5. Moreover, the method based on the distributed KD-tree structure developed by Luke [24] is also applied to improve computational efficiency for the case with numerous boundary nodes.

(2)

Where r 0,cg the coordinate of the centroid at the initial time, and the transformation matrix in Eq. (2) can be expressed as φ

− − →

where Ξ (B) represents the boundary point set, w i is the displace− −− → ment of the interior mesh node i, w b is the displacement of the boundary mesh node b which can be exactly obtained through the

d(Q n V n )



 r rot = T B − I • (r 0 − r 0,cg )

349







2t (7)



¯ Q nk+1 , Vk+1 = R Q nk+1 , V k+1 − R V Q nk+1 R grid grid



The superscript k represents the number of the time step and the conservative variables of the current time step Q n can be solved with the pseudo-time iterative method. Moreover, AUSM+ scheme is used for spatial discretization with an implicit LU-SGS scheme for temporal integration. The oneequation Spalart–Allmaras turbulence model is implemented to treat the turbulent boundary layers for computations of the viscous flow. More details about the CFD solver can be found in Ref. [25]. 2.3. Solution of flight dynamics equations and aeroelastic equations As illustrated in Fig. 1, the loosely coupled method is used for the research on flight dynamics of flexible aircrafts. When finishing the solution of the unsteady flow at a real time step, we obtain the aerodynamic parameters all over the flow field.

350

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

With the aerodynamics forces and moments over the wall surfaces, the movement of the objects can be calculated by numerically integrating the flight dynamic equation of SDOF motion, while the structural deformation can be obtained by solving the aeroelastic equation. The nodal positions of the new meshes can be regarded as the summation of the rigid-body movement caused by SDOF motion and the elastic deformation due to aeroelasticity. According to Waszak’s study, it is found that nonlinear inertial coupling can be ignored if the following characteristics are reflected in the flexible vehicle [13]: 1) The magnitude of aerodynamic load is not smaller than that of the gravity, 2) The order of the expected rotational rates (i.e. the shortperiod frequencies of flight dynamics) is smaller than that of the elastic frequencies, 3) The model geometry is not very complex and transverse deflections won’t result in notable changes in mass distribution. When the above restricted conditions are met, the additional effect accounting for variations in the inertia terms (including the location of the centroid and the moments of inertia) with structural deformation (deflection induced changes in mass distribution) and the coupling associated with the product between structural displacement and rate are both zero. Furthermore, when these terms are neglected, the dynamic modeling in Buttrill’s study is completely analogous to that of Waszak’s study, where the rigidbody SDOF and elastic degrees of freedom of flexible vehicles are decoupled. 2.3.1. Rigid body motion equations of SDOF The rigid body motion is disassembled into the centroid translation and the rotational motion around the centroid. The centroid motion can be calculated by Newton’s Law as follows I mx¨ cg = FaI + FeI + F Ig

(8)

I where m is the mass of the object, x¨ cg is the acceleration vector

of the centroid, and FaI , FeI , F Ig are the aerodynamic force, external forces (such as the thrust of the engine) and the forces of the gravity, respectively. Moreover, the superscript I represents the inertial coordinates. The dynamic equation of the rigid body rotation in the body coordinate can be readily written as

dH dt

=M

(9)

where H is the vector of the moment of momentum around the rigid centroid and M is the resultant moment vector caused by the aerodynamic moments and the other external forces. In addition, H could also be expressed in the body coordinate as Eq. (10).

H = Iω

(10)

where ω = [ωx , ω y , ωz ] is the angular velocity vector in the body coordinate and the inertia matrix I can be expressed as T

⎡

Ix I = ⎣ − I xy − I xz

− I xy Iy

− I yz

⎤ − I xz − I yz ⎦ Iz

where I x , I y , I z are the rotational inertia about the x, y and z axes in the body coordinate respectively, and I xy , I yz , I xz are the corresponding products of inertia respectively. When the principal axes of inertia are set as the initial axes of the body coordinate, the inertia matrix I becomes diagonal. While Eq. (9) is transformed into the formula under the body coordinates, it can be written as

δH +ω×H=M δt

(11)

Compared to the translation of the centroid, the rotational motion is more readily computed in the body coordinate to keep inertia properties time-invariant. Thus, the angular velocities can be calculated via

ω˙ = −I−1 ˙Iω − I−1 ω ss Iω + I−1 M

(12)

ss

where ω is the skew symmetric matrix of pression as

⎡

0

ω ss = ⎣ ωz −ω y

−ωz 0

ωx

ω , and it has the ex-

⎤

ωy −ωx ⎦ 0

It is also worth noting that inertia matrix I in the body coordinate is approximately constant based on the reasonable assumption that the effect of the elastic deformation on the mass distribution and the moments of inertia can be neglected. Then Euler motion equation is added to obtain the rate of the orientation angles.

ϕ˙ = Tω

(13)

where ϕ = ( ϕ θ ψ ) T and T is the transformation matrix which can be expressed as

⎡

1 sin ϕ tan θ cos ϕ T=⎣0 0 sin ϕ / cos θ

⎤

cos ϕ tan θ − sin ϕ ⎦ cos ϕ / cos θ

(14)

Combining Eq. (8), (12) and (13), we obtain the system of Eq. (15) for SDOF motion regardless of time derivatives of inertia matrix I

 ⎧ I ⎨ mx¨ cg = F ϕ˙ = Tω ⎩ ω˙ = −I−1 ω ss Iω + I−1 M

(15)

A hybrid linear multi-step scheme is applied to solve the firstorder ordinary differential equation expressed in Eq. (15). Secondorder semi-implicit linear multi-step scheme is employed for the time marching of flight dynamic equations [26]. 2.3.2. Static aeroelastic equation based on the mode-summation method Structural deformation needs to be obtained by means of computational structural dynamics (CSD) technique coupled with CFD approach when the load-carrying structure is of high nonlinearity or when large geometric deformation exists. Unfortunately, much computational time should be paid on the iterative process of CSD. Nevertheless, while the elastic deformation is much smaller than the characteristic length of the flight vehicle and the structure is of good linear characteristics, the structural deformation could be introduced from a modal decomposition of the discrete finite element model. Thus it leads to an approximately linear model based on reduced basis of modal coordinates, which could be titled as the modal approach [27,28]. This method has been used as one of the most practical approach during the analysis of structural dynamics of linear systems for its higher efficiency compared to CSD technique. Since the scale of the elastic deformation in this paper is much smaller than the characteristic length of the rigid projectile model, the modal approach is adopted to obtain the structural deformation. While N order structural modes are used to describe the structural deformation in the modal approach, the elastic physical deformation of the wall is subsequently determined by

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

w(t ) =

N 

Φ i qi (t )

351

(16)

i =1

where w is the displacement due to the structural deformation, Φ i is the mass-orthogonally scaled modal basis which could be obtained by FEM or modal tests and it should be interpolated to the nodes on the wall boundary of the fluid meshes only once at the pre-processing step. q i is the generalized coordinate of the ith-order mode. The final aeroelastic equations based on the modal approach could be written as

Mq¨ + Cq˙ + Kq = Q

Fig. 2. Geometric parameter and inertia properties of spinning projectile.

(17)

where M, K and C are the generalized mass matrix, the generalized stiffness matrix and generalized structural damping matrix, respectively. The generalized force Q can be calculated by

¨

Qi =

P Φ i · dS,

i ∈ [1,

···,

N]

s

where P (x, y , z, t ) is the pressure distribution on the wall surface. Noticeably, Eq. (15) is constructed in the body coordinates, which could keep the structural modes invariant during the free flight process. Furthermore, when static aeroelasticity is dealt with, all terms of Mq¨ and Cq˙ should be set to 0. Hence, Eq. (17) is simplified and we obtain an algebraic equation expressed as Eq. (18) which can be solved straightly.

Kq = Q

(18)

Combining Eq. (15) and (18), we can obtain the system of Eq. (19) which is needed in the fluid–solid interaction module in Fig. 1.

⎧ I  x¨ cg = F/m ⎪ ⎪ ⎪ ⎨ ϕ˙ = Tω ⎪ ω˙ = −I−1 ω ss Iω + I−1 M ⎪ ⎪ ⎩ q = K− 1 Q

(19)

Here, the reason for replacing the dynamic aeroelastic equation (17) with the static aeroelastic equation (18) is that the natural frequency of the structure is much larger than the frequency of the short-period mode. Previous study [29] indicated that quasistatic assumption can be adopted to simply the dynamic equation of flexible flight vehicles. That is to say, the structural deformation immediately converges to the steady value without dynamic vibration. Since the time step adopted for time marching of Eq. (19) depends on the frequency of the short-period mode but not the natural frequencies of the structure, this quasi-static assumption can dramatically increase the time step and computational efficiency. During time marching for the simulation of the free flight process, the system of Eq. (19) is solved individually after the computation of the unsteady flow field. The interference between the elastic deformation and the displacement of the SDOF motion is coupled through their contribution to the new positions of the mesh nodes in the moving and deformation zone. 3. Computational cases and analysis Because of the limitation of experimental facility and the space of the test section in wind tunnel, it is difficult to simulate the free flight process on the long-time history, especially for the flexible vehicles. Very few cases involved in the aerodynamics/structure/SDOF-motion coupled simulation can be found in the open

Fig. 3. Computational mesh of spinning projectile.

literatures. Accordingly, the coupling issue related to aerodynamics, structural mechanics and flight dynamics is decomposed into the two issues related to aerodynamics/structure coupling and aerodynamics/SDOF-motion coupling, respectively. 3.1. Verification case of spinning projectile An ARL spinning projectile configuration is selected as the computational model because of its detailed experimental data [30,31]. The geometric parameters can be seen in Fig. 2 where the unit is millimeter. The supersonic projectile model in this case is an ogivecylinder-finned configuration. More geometric parameters can be found in Ref. [31]. Fig. 2 also presents the inertia properties of the finned projectile. The computational mesh of the projectile is illustrated in Fig. 3. Table 1 presents the initial conditions of the trajectory simulation. The equivalent free flight Mach number is 3. The Reynolds number is 7.083 × 107 /m. Fig. 4 shows the variation of the attitude angle with the distance traveled. Both the amplitude and the frequency

352

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

coordinate system is also shown in Fig. 6 and it also agrees well with the result of Ref. [24] during the simulation process.

Table 1 Initial conditions of finned projectile. Direction

Distance of C.G. (meter)

Speed of C.G. (m/s)

Euler angle (rad)

Angular velocity (rad/s)

X Y Z

4.593 0.200 0.159

1030.81 22.064 86.278

2.051 0.088 −0.023

2518.39 −52.802 22.233

agree well with the experimental data. The pitching angle is of cyclical fluctuation. The amplitude decreases with the distance traveled. The yaw angle also changes periodically. The amplitude decreases in the beginning and then goes into stable oscillation. Fig. 5 shows the computed y and z distances as a function of the range. The change trend is basically the same. Compared to the lateral displacement in the y direction, the calculated result in the z direction is closer to the experimental data. Moreover, the oscillation of the displacement in the z direction is larger than that of the y direction. The reason is that the amplitude of the pitching angle is greater than that of the yaw angle during the flight process. It leads to the aerodynamic force in the z direction changing more dramatically than that of the y direction. In addition, the orientation angles decrease due to air resistance and the oscillation amplitudes of the aerodynamic forces/moments decrease accordingly. The y and z displacement curves also reflect the attenuation trend of the oscillation. Fig. 6 shows the comparison of the predicted aerodynamic forces in three directions with the data extracted from Ref. [24]. The x direction aerodynamic force agrees well with the reference data when x < 30 m. When x > 30 m, the amplitude and the period gradually decline and the calculated value is a little smaller than the reference results. The aerodynamic forces in the y and z directions both fit well with the reference results and they are very similar to the trend of the yaw angle and the pitching angle respectively, which means that the aerodynamic forces are mainly due to the orientation angles. The pitching moment in the inertial

3.2. Verification case of static aeroelasticity The aeroelastic study of the HIRENASD wing-body [28] is selected to evaluate the accuracy of aeroelastic simulation of the solver presented in this paper. The configuration and the hybrid meshes (see Fig. 7) are provided by the AePW Conference sponsors. This simulation is performed at a Mach number of 0.80 and a Reynolds number of 1.4 × 10e7 , since interesting transonic flow phenomena and significant wing deformations were found in the experiment [28,32]. The level of the dimensionless dynamic pressure q/E is 0.47 × 10e−6 , where q is the dynamic pressure of the free stream and E is the elastic modulus of the structure. Fully turbulent solution is adopted for viscous flow with SA model and the IDW method is applied for mesh deformation. Fig. 8 plots the comparison between the calculated elastic deformations and the experimental data at the leading edge and the trailing edge. There are significant deformations at the wingtip and the calculated deformations agree well with the experimental data, which gives a reliable impression of the numerical simulation accuracy concerning the flow field and the interaction with the elastic structure. 3.3. Free flight case of flexible projectile In order to improve the speed and the range of modern projectiles, engineers tend to increase the slenderness ratio and to decrease the structural weight ratio for more fuel space. However, these structures are prone to elastic deformation due to the low rigidity of the structure. Accordingly, the trajectory of the flexible projectiles cannot be predicted accurately regardless of the effect of elastic deformation. Here, a projectile model similar to AIM-7 is constructed to study the influence of aeroelasticity on flight dynamics.

Fig. 4. Attitude angle vs. range.

Fig. 5. Trajectory simulation vs. range.

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

353

Fig. 6. Aerodynamic forces/moments vs. range.

3.3.1. Computational model The computational geometry and corresponding parameters are illustrated in Fig. 9. The mesh for CFD solver is shown in Fig. 10 and Fig. 11. Three meshes were generated for this study, and the grid statistics are summarized in Table 2. The surface grid resolution near the region of the wing for each of these meshes is shown in Fig. 11. As mentioned in Section 3.1, Gambit and TGrid are used to generate the surface meshes and the volume meshes respectively. Although the densities of surface mesh in the computational domains are different as given in Fig. 11, the ratios of the volume mesh keep the same. The distance between the first layer and the wall in the perpendicular direction is 1.4 × 10−5 relative to the projectile diameter, to ensure that the dimensionless length y + is around 1. There are 30 layers of prism cells in the boundary layer region with a growth rate of 1.2. The free flight event is simulated at a Mach number of 1.5 at an altitude of 8000 meters and the Reynolds number is about 1.51 × 107 based on the reference length of the projectile diameter. Here, according to the mass distribution of the structure, the inertial properties of the projectile are computed with the software Patran. Since the displacement due to the elastic deformation is much smaller than the characteristic length of the projectile, the centroid position and the moments of inertia are assumed to be constant during the flight process, which will be validated in the following section. It is necessary to get the generalized mass, natural frequencies, and the elastic mode shapes of the model when using the mode-summation method. MSC Nastran is adopted here to undertake modal analysis and obtain the modes of the first sixteen orders including the main dynamic character of the body, the wings and the control surfaces. Fig. 12 presents the first six order mode shapes. For the projectile model constructed with X–X configuration and axial symmetry, only the three degrees of freedom of the longitudinal motion is researched here. Fig. 13 and Fig. 14 depict the centroid displacements and the orientation angles of different meshes as a function of time. Noticeably, the displacement in the x direction represents the distance between the free flight model

and the reference vehicle traveling at a constant speed of the initial condition. It is apparent that the displacements in z (vertical) direction match very closely for all three grid refinements. This is mainly because the normal aerodynamic force is less sensitive to mesh resolution compared to the drag. In the x direction the agreement is not very good. The discrepancy in the x direction is expected because the drag of viscous effect is overestimated with the coarse mesh. Although some deviation exists in the amplitudes of the orientation angles of different meshes, the oscillation periods match well. Comparison data indicates that as the density of the mesh is refined, the spatial resolution of the solver is enhanced, and qualitative consistency of the numerical results performs well. Therefore, the medium grid is adopted in this research considering the balance between efficiency and accuracy. In addition, as shown in Fig. 14, it is noticeable that frequency of the longitudinal oscillation is about 1 hertz, which is much smaller than the natural frequencies of the bending modes (17.35 hertz). Accordingly, the second assumption of Waszak’s study is reasonable when applied in this paper. 3.3.2. Results and analysis The centroid position and the velocity of the rigid projectile are compared with the results of the elastic model in Fig. 15 and Fig. 16. The trajectories of the two models are qualitatively the same. The difference increases as the free flight proceeds. Noticeably, since the elastic deformation change the streamlined shape of the original rigid model, the drag is increased which leads to more reduction of the flight speed. The pitching angle and the angular velocities around the centroid are depicted in Fig. 18 and Fig. 19, respectively. Although there is obvious difference during the response process of the two models, the final steady state of the orientation angles and the angular velocities are almost the same. The oscillation of the pitching angle attenuates, since the configuration of the projectile is stable. However, the oscillatory cycles increase when elastic deformation is taken into consideration, which is consistent with the conclusion of Ref. [33]. The reason for this phenomenon is that the aero-

354

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

Fig. 7. Boundary mesh and clip of HIRENASD model.

Fig. 10. Computational mesh of the projectile. Table 2 Grid statistics of different refinements.

Fig. 8. Comparison of the deformation between the calculated and experimental results (α = 3◦ ).

Fig. 9. Geometric model and parameter of the projectile with large slenderness ratio.

dynamic center of the flexible has moved because of the elastic deformation, resulting in the change of the stability margin. Responses of aerodynamic forces/moments are illustrated in Fig. 17 and Fig. 20, respectively. Aerodynamic oscillatory cycle of the elastic model is longer than that of the rigid model, which is consistent with the results of the velocities (see Fig. 16) and the orientation angles (see Fig. 18). Although the response histories are different, the steady states of the elastic and the rigid models agree well with each other. It’s worth noting that the oscillation of the pitching angle of the elastic model is weakened compared to that of the rigid model due to the vibration-absorbing effect of the elastic deformation. Fig. 21 presents the generalized coordinates of the first six structural modes. It is obvious that the first two bending modes of the body play a dominant role in the structural deformation. The generalized coordinates of the third and fourth order modes are much smaller than that of the first two modes. The amplitudes

Mesh refinement

Node

Cell

Coarse Medium Fine

116152 268446 445480

612818 1294038 2368426

of the generalized coordinates oscillate at the beginning and converge to the steady states as the flight proceeds. Here, with the lump mass model to describe the mass distribution along the axial direction of the projectile (as shown in Fig. 22), posterior estimate is adopted to compute the variation of the flexible projectile’s inertial parameters with the deformation in Fig. 21 by means of Eq. (20) and Eq. (21)

⎧  N  N Lm  ⎪   ⎪ ⎪ ⎪ xcg = mj q i φi j , x m ⎪ ⎪ ⎨ j =1 i =1  N  N Lm ⎪    ⎪ ⎪ ⎪ ⎪  z = m q φ m cg j i i j , z ⎪ ⎩ j =1 i =1 ⎧ N Lm ⎪ 

⎪ ⎪ ⎪ m j r 2j ,x − r02, j ,x  I xx = ⎪ ⎪ ⎨ j =1

N Lm ⎪  ⎪

⎪ ⎪ ⎪ I = m j r 2j , y − r02, j , y  y y ⎪ ⎩

(20)

(21)

j =1

Where m j is the mass of the j-th lump mass element, N Lm is the number of the lump mass elements, φi j ,x is the x component of i-th order structural mode at the j-th lump mass element. r j ,x is the distance between the j-th lump mass element and y axis of

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

355

Fig. 11. Surface meshes for coarse, medium, and fine resolution.

flexible projectile, while r0, j ,x is the corresponding distance of the rigid projectile. Similarly, r j , y and r0, j , y can be obtained. Fig. 23 and Fig. 24 present the variation of the centroid location and the moments of inertia of the flexible projectile. Due to the lower bending modes, the maximum variation of the centroid in z direction is about 3 millimeters, i.e., smaller than 0.1% of the length of the projectile, while the change of the centroid in x direction is very close to zero. The variation of I xx is larger than that of I y y , which is mainly caused by the misalignment value  zcg due to the bending deformation. Overall, the relative variations z x of the inertial parameters are quite small ( Lcg < 0.1%, Lcg <

Fig. 12. First six order modes and natural frequencies.

I 0.01%, I I xx < 0.3%, I y y < 0.02%), which can be neglected during xx

yy

the flight process. In contrast, the variation of the inertial parameters in Ref. [34] is much larger (with the magnitude up to 1% or even larger). Thus, the third assumption of Waszak’s study is also reasonable when applied in this paper. In order to understand the mechanism by which the effect of the elastic deformation affects the flight dynamics of flexible projectiles, the trajectory of another rigid model with the centroid moving backward is investigated. The position of the centroid is 2.00 meters away from the nose of the body (i.e., 0.10 meter backward compared to the model in Fig. 9). Computational results of the pitching angle shown in Fig. 25 indicate that the response history of the rigid model with the centroid moving backward (marked with “_rigid_bw”) is much closer to the result of the flex-

Fig. 13. Centroid displacements of different mesh resolution.

ible model with the original centroid, despite some differences in the amplitudes. Fig. 26 presents the aerodynamic center relative to the centroid during the flight process. It can be found that all the three models are static stable during the flight process, (i.e.,

356

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

Fig. 14. Pitching angle of different mesh resolution,

Fig. 18. History of orientation angles.

Fig. 15. History of the position of the centroid.

Fig. 19. History of angular velocities.

Fig. 16. History of centroid velocities.

Fig. 17. History of aerodynamic forces.

Fig. 20. History of aerodynamic moments.

Fig. 21. General displacement of flexible model.

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

357

Fig. 22. Lump mass model of the projectile.

Fig. 25. History of orientation angles of various models.

Fig. 23. Variation of centroid of the flexible projectile.

Fig. 26. Aerodynamic center relative to center of gravity of various models.

Fig. 24. Variation of moments of inertia of the flexible projectile.

Fig. 27. Elastic deformation sketch of the flexible projectile model.

xac − xcg > 0). At the beginning of the flight process, the relative position of the aerodynamic center of the flexible model is much closer to that of the rigid model with the centroid moving backward. In comparison, the average relative position of the aerodynamic center of the original rigid model is much larger than the above two models (see Fig. 26). Compared to the oscillation of the aerodynamic center mainly caused by the instant pitching angle of the rigid model, the fluctuation of the aerodynamic center of the flexible model is mainly due to the effect of the elastic deformation. As the flight process goes on, the oscillation of the pitching angle converges to a stable value. Both the aerodynamic centers of the flexible model and the original rigid model converge to the same stable position. And the dispersion of the aerodynamic centers of the flexible model and the one with the centroid moving backward is equal to the variation of the centroids. Moreover, compared to the variation of the centroid in Fig. 23, the variation of the aerodynamic center due to the elastic deformation is much larger. Thus, the effect of the structural deformation on the stability margin (i.e., the value of (xac − xcg )) mainly depends on the change of the aerodynamic center but not the location of the centroid. Comparison results could provide the explanation that the elastic deformation moves the aerodynamic center forward. From Ref. [35], the angular velocity of the oscillation in the longitudinal direction could be derived in Eq. (22)

ω2Lg =

1 Iy





q∞ S W C AW C L α (xac − xcg ) − C mq /2μ)

(22)

Where ω Lg is the angular velocity of the longitudinal oscillation, q∞ is the dynamic pressure of the free stream, S W and C AW are the wing area and the aerodynamic chord of the wing respectively, C L α is the slope of the lift coefficient relative to the angle of attack, xac and xcg are the aerodynamic center and the centroid position respectively, C mq is the slope of the longitudinal moment relative to the pitching rate, μ is the dimensionless density. For the projectile model researched in this paper, q∞ , S W , C AW and μ are unrelated to the rigidity of the models, while C L α and C mq are less sensitive to the deformation. Thus, ω Lg is mainly determined by the distance between the aerodynamic center and the centroid position. Since the centroid of the original rigid model is in front of the aerodynamic center, i.e. and the configuration is stable, the elastic deformation decreases the distance between the centroid and the aerodynamic center, which could decrease ω Lg and increase the oscillation cycle correspondingly. As illustrated in Fig. 21, the elastic deformation is mainly caused by the first two bending modes presented in Fig. 12. Detailed research is focused on understanding the mechanism of the elastic deformation on the movement of the aerodynamic center. Fig. 27 depicts the deformation sketch of the flexible projectile due to aerodynamic loads, where α0 represents the angle of attack relative to the rigid model, τ (x) is the torsion angle due to the elastic deformation and τt is the torsion angle at the tail. The local angle of attack and the deflection angle of the control surface can be approximately defined as

358

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

Fig. 28. Pressure coefficient distribution on the symmetry clip (dashed line—rigid, solid line—flexible, y = 0 and α0 = 6◦ ).

⎧ N  ⎪ ∂φi ,z ⎪ ⎪ α ( x ) = α − τ ( x ) = α − qi ⎪ 0 0 ⎪ ⎨ ∂x i =1

N ⎪  ⎪ ∂φi ,z ⎪ ⎪ θ = θ − qi |x=xt ⎪ cs,0 ⎩ cs ∂x

(23)

i =1

Where θcs,0 is the deflection angle of the control surface, θcs is the effectual deflection angle when the torsion angle is taken into consideration. Bounded by the position of the centroid, the torsion angle of the fore body of the flexible projectile is negative, while the one of the aft body is positive. Fig. 28 presents the comparison of the pressure coefficient distribution on the symmetry clip. Here, three segments are analyzed for the effect of the elastic deformation on the pressure coefficient distribution: 1) The pressure on the lower surface of the fore body increases for the increment of the local attack angle due to the torsion angle. And the increment of the pressure becomes larger as the position is closer to the nose. This effect will increase the lift and the upward moment. 2) As to the aft body relative to the centroid, the pressure on the lower surface of the aft body decreases for the reduction of the local attack angle due to the torsion angle. And the reduction of the pressure becomes more remarkable as the position is closer to the tail. Compared to the rigid model, this effect will reduce the lift and generate more upward moment. 3) For the elevator at the tail of the body, the torsion angle is more prominent. While added to the deflection angle, the effectual deflection angle will produce more upward moment than that of the rigid control surface. Thus, a smaller deflection angle is needed to trim the flexible projectile. All the three factors above will generate extra upward momentum compared to the rigid model. Since the increment of the lift provided by the fore body makes up for the reduction of the aft body, the lift of the whole body will not change so markedly. Consequently, the aerodynamic center of the flexible projectile will move forward compared to that of the rigid model, resulting in low stability margin. Furthermore, if the stability margin is much smaller or the structural stiffness is much lower, the aerodynamic center of the flexible model might be moved to the front of the centroid, which could even make the aerodynamic configuration unstable. And attention should be paid to the design of the control system for the effect of the elastic deformation accordingly. 4. Conclusion In this paper the rigid-motion mesh technique and an IDW mesh morphing method are adopted to treat the flight dynamics of a flexible projectile. URANS equations are solved with CFD technique to obtain the aerodynamic forces/moments. Flight dy-

namic equation coupled with static aeroelastic equation is solved. According to the comparison results of the projectile models of X–X configuration, some conclusions are obtained as follows: 1) The deformation of the model is dominated by the lower order bending modes of the body and could converge to the steady values as the flight proceeds. 2) The variations of the inertial parameters due to elastic deformation are much smaller compared to the length of the projectile. And the effect of the structural deformation on the stability margin mainly depends on the change of the aerodynamic center but not the location of the centroid. 3) The elastic deformation will not qualitatively change the trajectory of the centroid. Because of the change the streamlined shape of the original rigid model due to elastic deformation, the drag is increased and leads to more reduction of the flight speed. 4) The effect of the elastic deformations increases the oscillation cycles and decrease the static stability margin. And the vibration-absorbing effect of the elastic model can weaken the oscillation amplitudes of the pitching angle. 5) As the elastic deformation changes the local attack angle of the projectile, the summation of the extra upward momentum of the body and the elevator makes the aerodynamic center move forward, which decreases the static stability margin. Thus, attention should be paid to the design of the control system of the flexible projectiles. In addition, the authors will extend the model and methodology and use more accurate structural model and rigid– elastic coupling equations for dynamic analysis in future work. Conflict of interest statement None declared. Acknowledgements The work presented in this paper is supported by the National Natural Science Foundation of China (Grant No. 11272262). References [1] J.E. Cochran, D.E. Christensen, Free-flight rocket attitude motion due to transverse vibration, J. Spacecr. Rockets 17 (5) (1980) 425–431. [2] T.R. Beal, Dynamic stability of flexible missile under constant and pulsating thrusts, AIAA J. 3 (3) (1965) 486–494. [3] L. Meirovitch, D.A. Wesley, On the dynamic characteristics of a variable-mass slender, elastic body under high accelerations, AIAA J. 5 (8) (1967) 1439–1447. [4] T.C. Zang, H.X. Hu, A review of great slenderness ratio projectile elastic effect research, J. Ballist. 11 (3) (1999) 89–93. [5] D.H. Latus, Aeroelastic stability of slender, spinning missile, J. Guid. Control Dyn. 15 (1) (1992) 144–191. [6] Z.H. Ma, Analysis of large-scale flexible fairing separation, Sci. China, Technol. Sci. 39 (2009) 482–489. [7] G.E. Reis, W.D. Sundery, Calculated aeroelastic bending of a sounding rocket based on flight data, J. Spacecr. Rockets 4 (11) (1967) 1489–1494. [8] W.C. Womack, C.W. Bert, F.J. Perdreauville, Dynamics of sounding rockets at burnout, J. Spacecr. Rockets 11 (10) (1973) 716–720. [9] L.M. Wang, An analysis on the flexibility in flight of projectiles or rockets having high L/D ratios, Acta Armamentarii 21 (2) (2000) 108–111. [10] I. Tuzcu, On the stability of flexible aircraft, Aerosp. Sci. Technol. 12 (2008) 376–384. [11] E. Oliveira, P. Gasbarri, I.M. Fonseca, Flight dynamics numerical computation of a sounding rocket including elastic deformation model, in: AIAA Atmospheric Flight Mechanics Conference, Dallas, 2011, AIAA Paper 2015-2710. [12] M.J. Li, X.T. Rui, L.K. Abbas, Elastic dynamic effects on the trajectory of a flexible launch vehicle, J. Spacecr. Rockets 52 (6) (2015) 1586–1602. [13] M.R. WaszakI.M., C.S. Buttrill, D.K. Schmidt, Modeling and Model Simplification of Aeroelastic Vehicles: An Overview, NASA Technical Memorandum 107691, 1992. [14] C.S. Buttrill, T.A. Zeiler, P.D. Arbuckle, Nonlinear simulation of a flexible aircraft in maneuvering flight, in: Proceedings of the AIAA Flight Simulation Technologies Conference, Monterey, CA, 1987, AIAA Paper 87-2501. [15] L. Yang, Z.Y. Ye, The interference aerodynamics caused by the wing elasticity during store separation, Acta Astronaut. 121 (2016) 116–129.

Hua R.H. et al. / Aerospace Science and Technology 71 (2017) 347–359

[16] L. Yang, Z.Y. Ye, J. Wu, The interference of the elastic vibration of the carrier to the aerodynamics of the external store in air-launch-orbit process, Acta Astronaut. 128 (2016) 440–454. [17] A. Schütte, G. Einarsson, A. Raichle, Numerical simulation of maneuvering aircraft by aerodynamic, flight-mechanics, and structural-mechanics coupling, J. Aircr. 46 (1) (2009) 53–64. [18] R.L. Meakin, C.A. Atwood, N. Hariharan, Development, deployment and support of a set of multi-disciplinary, physics-based simulation software products, in: 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Orlando, 2011, AIAA Paper 2011-1104. [19] S.A. Morton, D.R. McDaniel, D.R. Sears, Kestrel v2.0-6DOF and control surface additions to a CREATE simulation tool, in: 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Orlando, 2011, AIAA Paper 2011-1104. [20] S.A. Morton, D.R. McDaniel, HPCMP CREA-TETM-AV kestrel current capabilities and future direction for fixed-wing aircraft simulations, in: 53rd AIAA Aerospace Sciences Meeting, Kissimmee, 2015, AIAA Paper 2015-0039. [21] S.I. Silton, J. Sahu, F. Fresconi, Comparison of uncoupled and coupled CFDbased simulation techniques for the prediction of the aerodynamic behavior of a complex projectile, in: 34th AIAA Applied Aerodynamics Conference, Washington, D.C., 2011, AIAA Paper 2011-1104. [22] L.K. Abbas, D. Chen, X. Rui, Numerical calculation of effect of elastic deformation on aerodynamic characteristics of a rocket, Int. J. Aerosp. Eng. 3–4 (2014) 1–11, http://dx.doi.org/10.1155/2014/478534. [23] J.A. Witteveen, H. Bijl, Explicit mesh deformation using inverse distance weighting interpolation, in: 19th AIAA Computational Fluid Dynamics Conference, AIAA, San Antonio, Texas, 2009, AIAA Paper 2009-3996. [24] E. Luke, E. Collins, E. Blades, A fast mesh deformation method using explicit interpolation, J. Comput. Phys. 231 (2) (2012) 586–601.

359

[25] G. Wang, New Type of Grid Generation Technique Together with the High Efficiency and High Accuracy Scheme Researches for Complex Flow Simulation, Northwestern Polytechnical University, Xi’an, 2005. [26] W.W. Zhang, Y.W. Jiang, Z.Y. Ye, Two better loosely coupled solution algorithms of CFD based aeroelastic simulation, Eng. Appl. Comput. Fluid Mech. 1 (4) (2007) 253–262. [27] W.T. Thomson, M.D. Dahleh, Theory of Vibration with Application, Prentice-Hall Inc, New Jersy, 1998. [28] J. Neumann, M. Ritter, Steady and unsteady aeroelastic simulations of the HIRENASD and wind tunnel experiment, in: IFASD 2009-132, 2009. [29] Z.D. He, H. Gao, Advanced Flight Dynamics, Northwestern Polytechnical University Press, Xi’an, 1990. [30] B.S. Davis, B.J. Guidos, T.E. Harkins, Complementary Roles of Spark Range and Onboard Free-Flight Measurements for Projectile Development. US Army Research Lab, Aberdeen Proving Ground, MD, Weapons and Materials Research Directorate, 2009. [31] J. Sahu, Time-accurate numerical prediction of free-flight aerodynamics of a finned projectile, J. Spacecr. Rockets 45 (5) (2008) 946–954. [32] https://c3.nasa.gov/dashlink/static/media/other/HIRENASD_base.htm. [33] R.H. Hua, C.X. Zhao, Z.Y. Ye, et al., Effect of elastic deformation on the trajectory of aerial separation, Aerosp. Sci. Technol. 45 (9) (2015) 128–139. [34] P. Gasbarri, R. Monti, C.D. Angelis, M. Sabatini, Effects of uncertainties and flexible dynamic contributions on the control of a spacecraft full-coupled model, Acta Astronaut. 94 (2014) 515–526. [35] Z.P. Fang, Flight Dynamics of Airplane, Beijing University of Aeronautics and Astronautics Press, Beijing, 2005.