Effects of rotational motion on dynamic aeroelasticity of flexible spinning missile with large slenderness ratio

JID:AESCTE AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.1 (1-14)

Aerospace Science and Technology ••• (••••) ••••••

1

Contents lists available at ScienceDirect

67 68

2 3

Aerospace Science and Technology

4

69 70 71

5

72

6

www.elsevier.com/locate/aescte

7

73

8

74

9

75

10

76

11 12 13

Effects of rotational motion on dynamic aeroelasticity of flexible spinning missile with large slenderness ratio

16 17

81

School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China

83

82 84 85

19

a r t i c l e

i n f o

a b s t r a c t

86 87

21 22 23 24 25

Article history: Received 23 January 2019 Received in revised form 19 June 2019 Accepted 3 September 2019 Available online xxxx

26 27 28 29 30 31

79

Li Heng ∗ , Ye ZhengYin

18 20

78 80

14 15

77

Keywords: Dynamic response Spinning missile Aeroelasticity Structural deformation and rate Mesh deformation

32 33 34 35 36 37 38 39

The structural rigidity of a spinning missile with large slenderness ratio is usually small, and the structural deformation and rate should not be ignored. Furthermore, rotational motion makes the aeroelasticity more complicated. Therefore, unsteady Euler equations and generalized dynamic aeroelastic equations are coupled simultaneously to simulate the dynamic aeroelastic response of a spinning missile with large slenderness ratio using rigid-motion mesh and radial-basis-function (RBF) morphing mesh techniques. The unsteady Euler equations are solved by computational fluid dynamics (CFD) technique by the in-house code. The Coriolis term and centrifugal loading term due to rotational motion are both considered in the generalized dynamic aeroelastic equations. The rigid-motion mesh and RBF morphing mesh techniques are both based on unstructured mesh, and the rigid-motion mesh is adopted to treat the rigid motion due to rotational motion, while the RBF morphing mesh is employed for flexible structural deformation caused by aeroelasticity. Numerical results of aeroelastic case are well agreed with the experimental results, which validates the numerical method. A missile model with X-X configuration is constructed to investigate the effects of rotational motion on dynamic aeroelasticity. The dynamic aeroelastic responses of the missile with and without rotational motion are simulated, respectively. Comparison results show that the lateral modes and longitudinal modes are coupled together because of rotational motion. In addition, the structural natural frequencies are changed due to rotational motion. In the end, detailed numerical analysis of the generalized dynamic aeroelastic equations used in this paper indicates the mechanism by which the rotational motion leads to the coupling of lateral modes and longitudinal modes and changes the structural natural frequencies. © 2019 Elsevier Masson SAS. All rights reserved.

88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105

40

106

41

107

42

108 109

43 44

1. Introduction

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

The net force of the missile in flight is not always located at the center of mass, therefore, the eccentric moment will try to flip the missile in air. The spinning flight is widely used by modern missiles with fixed tails because of its advantages as follows. Firstly, it can overcome or reduce ballistic dispersion caused by thrust eccentricity, mass eccentricity and aerodynamic eccentricity, and consequently it will enhance the attack accuracy. Secondly, it is possible to simplify control system composition for the missile with low rotational speed, because both pitching and yawing channels can be controlled by one control surface, thus related control equipment can be reduced. Thirdly, it is helpful for a reentry vehicle to avoid asymmetric ablation and improve the penetration ability. In recent years, greater missile payload and range are increasingly required with the development of missile technology,

60 61 62 63 64 65 66

*

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

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

therefore, the lighter structure and larger slenderness ratio of missile are more urgently needed than before. But unfortunately, the lighter structure and larger slenderness ratio consequently result in lower structural stiffness, smaller natural frequency and larger structural deformation. As for the missile with small slenderness ratio, aeroelasticity is not very important, and the missile could be treated as a rigid body. However, as for the missile with large slenderness ratio, aeroelasticity should be considered carefully. As a flexible body, a missile without rigid-body motion is a pure aeroelastic problem which could be solved by coupling aerodynamics equations and structural dynamics equations [1–8], because the process is only associated with aerodynamics and structural dynamics. However, as for the missile with rigid-body motion, especially the spinning missile with large rotational angular velocity, it is more complicated and more disciplines should be considered in the dynamic process. Since the problem is involved in several interactions among aerodynamics, flight mechanics and aeroelasticity, previous studies mainly handled the problem with simplified models which had various assumptions [9–22]. Reis et al. [9] investigated aeroelastic bending of the sounding rocket consisted of

110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.2 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

two rigid bodies hinged at their interface and pointed out that the transverse bending of the free-flight rocket could lead to obvious lateral angular velocity. Buttrill et al. [11] proposed a mathematical model which integrates nonlinear rigid-body flight mechanics and linear aeroelastic dynamics, and the model considers aerodynamic coupling and inertial coupling, and can be simplified to Waszak’s model [10] when meeting some extra assumptions, in addition, Buttrill’s model was further improved by utilizing nonlinear straindisplacement relations [12]. Waszak, Buttrill and Schmidt [13] further discussed the method to simplify these theories. Platus et al. [14] derived the complete coupled equations of motion that govern the static and dynamic aeroelastic stability of flexible spinning missiles. Furthermore, a lot of mathematical models [15–22] were also presented and simplified based on various assumptions to investigate these questions. Over the years, with the development of computer science and numerical algorithm, CFD has been introduced to investigate these problems [6,23–33]. Blades et al. [23], Abbas et al. [6], Yin et al. [28,29] and Yang et al. [30,31] concentrated on aerodynamics and investigated the effects of aeroelasticity or given structural deformation on the aerodynamic performance of the aircraft or the missile. Schütte et al. [24] coupled aerodynamic, flight mechanics and aeroelastic computations to simulate the unsteady aerodynamics of a free-flight aeroelastic combat aircraft. Hua et al. [25,26] coupled aerodynamics equations, rigid body dynamics equations and static aeroelastic equation based on modal superposition method to investigate the effects of elastic deformation on mass location, inertial moment and missile trajectory, furthermore, flight dynamics based on closed-loop feedback control was also studied [33]. Li et al. [27] employed quasi-steady CFD method to construct an aerodynamic model, and coupled the model with flexible dynamics equation derived by Lagrange equation to investigate the effects of elasticity on missile trajectory and vibration. Lu et al. [32] studied the aerodynamic coupling effects of spinning and coning motions for a finned vehicle. The aerodynamic coupling effects are mainly produced by the shockwaves and expansion waves around fins, which are significant at low Mach numbers and high coning and spin rates. Although a lot of researches have been conducted to investigate the aeroelasticity of the flexible missile and aircraft, there are still some limitations. Firstly, the aerodynamic loads which dominate the structural deformation and structural vibration were normally predicted with engineering methods such as lifting-line method, but more exact methods such as CFD could be used to investigate these problems, and there were some successful applications of CFD on the aeroelasticity [2,3,5–8,23–33]. Secondly, many researches involved with coupling of elastic deformation and rigidbody motion focused on the change of aerodynamic characteristics or flight trajectory due to structural static deformation or given elastic deformation, but the dynamic response of structural vibration were rarely considered. Thirdly, the aeroelastic investigations involved with rigid-body motion mainly focused on translation of rigid body, that is to say, the angular velocity is usually negligible, and the effects of rotational motion on structural vibration are less investigated. Due to the limitations mentioned above, in this paper, based on the frame of unstructured hybrid mesh, the rigid-motion mesh and RBF mesh deformation techniques are adopted to treat the rigid motion due to rotational motion and flexible structural deformation caused by dynamic aeroelasticity, respectively. Unsteady Euler equations are solved to obtain the instant aerodynamic forces at different physical time. Generalized dynamic aeroelastic equations composed of rotational motion and dynamic aeroelastic equation are constructed and solved to simulate the dynamic aeroelastic response of the spinning missile with large slenderness ratio. The investigation is focused on the effects of the rotational motion on

dynamic aeroelasticity, and the dynamic aeroelastic responses of the flexible missile with and without rotational motion are simulated, respectively.

67 68 69 70

2. Numerical method

71 72

2.1. Coupling strategy and framework of the aeroelasticity of flexible spinning missile

73 74 75

The flowchart to solve the aeroelastic problem of the flexible spinning missile is shown in Fig. 1. When solving this problem, three modules are mainly employed: grid-treating module, CFD solver module and aeroelastic module of flexible spinning missile. The grid-treating module contains two submodules: rigidmotion mesh submodule and mesh deformation submodule. The rigid-motion mesh submodule is able to adjust the new computational domains according to the large deformation due to translation and rotation. The mesh deformation submodule is capable of adapting elastic deformation of the meshes caused by aeroelasticity. The CFD solver module is used to obtain the numerical solution of the fluid governing equations at every physical time step. In this paper, inviscid Euler equations is solved based on unsteady Reynolds-averaged Navier-Stokes (URANS) equations. With the solution obtained in CFD solver, the aerodynamic forces on the spinning missile are computed by integrating the pressure on the surface of the missile. At the same time, generalized forces needed by the solution of the dynamic aeroelastic equations could be obtained. The aeroelastic module of flexible spinning missile is mainly used to solve generalized aeroelastic equations, and the module is loosely coupled with the CFD solver module. The aeroelastic module of flexible spinning missile contains two submodules: rigidbody motion submodule and structural dynamics submodule. The result of this module can be divided into two parts: rigid-body motion result due to large rigid-body motion and structural vibration result caused by aeroelasticity. The rigid-body motion result contains rigid displacement, which combined with rigid-mesh motion method to adjust the computational domains. Structural vibration result contains boundary deformation, which combined with mesh deformation method to adapt volume mesh. The displacement of any mesh point includes two parts: the deformation resulting from rigid-motion mesh and the deformation resulting from mesh deformation. Correspondingly, the flow field is composed of synthetic effects of both rigid-body motion and aeroelasticity. Various components of this framework are described as follows.

76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113

2.2. Grid-treating module

114 115

With the rapid progress of computer science and numerical algorithm, more detailed research of the flow field around complicated configurations can be treated with CFD technique. It is well known to all that unstructured meshes are more flexible and convenient than structured grids to handle models. Hybrid unstructured meshes are used for spatial discretization of the flow field in this paper. The rigid-motion mesh method is employed to treat the large rigid-body displacement caused by translation and rotation, while the mesh morphing method is used to treat the flexible structural deformation due to aeroelasticity.

116 117 118 119 120 121 122 123 124 125 126

2.2.1. Rigid-motion mesh method The basic idea of the rigid-motion mesh method is to update the location of all nodes of the computational mesh according to the translational displacement of the center of mass and the rotational displacement around the center of mass. The distinct advantage of the rigid-motion mesh method is that the computational

127 128 129 130 131 132

JID:AESCTE AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.3 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

3

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

Fig. 1. Flowchart of the solving methodology.

33

⎡

34

1 0 TB-I = ⎣ 0 cos φ 0 sin φ

35 36

⎡

37

cos ψ × ⎣ sin ψ 0

38 39 40

⎤⎡

0 cos θ − sin φ ⎦ ⎣ 0 cos φ − sin θ

⎤ − sin ψ 0 cos ψ 0⎦ 0

99

⎤

100

0 sin θ 1 0 ⎦ 0 cos θ

101 102 103 104

(3)

1

Where φ , θ and ψ are rolling angle, pitching angle and yawing angle, respectively. In this paper, only the rotational term is used and the translational term is zero.

42 43 44

Fig. 2. Vector definition.

46 47

50 51 52 53 54 55

mesh can be obtained directly without being reconstructed during the numerical simulation process, and what we need to do is just some transformation on coordinate position. The method is very efficient, and is able to keep the topological relationship and the quality of the computational mesh. As shown in Fig. 2, in inertial frame, define r0 = [x0 , y 0 , z0 ]T as initial position of mesh vertex, and new position at any time r = [x, y , z]T can be defined as

56 57 58 59 60 61 62 63 64 65 66

108 109 110 111

45

49

106 107

41

48

105

r = r0 + rcg + rrot

(1)

Where, rcg is the translation of the center of mass, and is same for all vertices in mesh. rrot is the displacement owing to rotation around the center of mass, which is different at different position, and can be written as

rrot = TB-I (r0 − rcg,0 ) − (r0 − rcg,0 )

(2)

Here, rcg,0 is the position of the center of mass, and the transition matrix from body frame to inertial frame TB-I is

2.2.2. Mesh deformation method The main purpose of the mesh deformation method is to diffuse the displacements of wall surfaces to the interior mesh nodes of the computational zone. The displacement of wall surfaces can be obtained from aeroelastic module. In this paper, mesh deformation method is based on RBF mesh deformation method [34,35]. In general, the number of nodes in boundary is large, therefore there will be a large amount of computational cost in mesh deformation. In order to improve the mesh deformation efficiency, a “double-edge” greedy supporting point selection algorithm [36] using a multi-level subspace method is adopted to reduce data. In the condition of meeting the requirement of the grid precision, this method can reduce computational cost. The RBF mesh deformation method with greedy point selection algorithm is robust and efficient, and the mesh quality of volume element in computational domains can be well kept.

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128

2.3. CFD solver module

129 130

The flow governing equation used in this paper is three dimensional unsteady Euler equations which is based on three dimen-

131 132

JID:AESCTE

AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.4 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

sional URANS equations. In order to consider the rigid-body motion and structural elastic deformation, the flow control equations described by arbitrary Lagrange-Euler method (ALE) are employed. For a grid cell with arbitrary motion and deformation, the integral form of URANS equations in Cartesian coordinate system is given as

∂ ∂t







∂



V

F Q, V grid · ndS =

Q dV + 



 F (Q) · ndS

(4)

∂

T

Where, Q = ρ ρ u ρ v ρ w e is conservation variable; ρ is density; u, v, w are velocity components in three different directions; e is total energy; F(Q, V grid ) is inviscid flux and FV (Q) is viscous flux; V grid is mesh velocity of a control volume due to rigid-body motion and structural deformation;  represents a control volume and ∂  represents the boundary of control volume ; n is outer normal direction of the boundary face ∂  of the control volume . In this paper, the finite volume method based on cell-center method are employed on hybrid unstructured mesh. Physical time is marched by dual time stepping method and pseudo time is marched by Gauss-Seidel iteration. Since the flow field of computational domains has been obtained, the aerodynamic forces can be obtained by integrating the pressure and skin friction over the body, and meanwhile, the generalized forces needed by aeroelastic module can be obtained too.

29 30 31

2.4. Aeroelastic module of flexible spinning missile

32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62

Aeroelastic module is mainly used to solve the generalized dynamic aeroelastic equations of flexible spinning missile. The input data is generalized forces, which can be obtained from CFD solver module, and the output data is rigid-body displacement and structural flexible deformation, which is needed by grid-treating module. The generalized dynamic aeroelastic equations of flexible missile are derived from Buttrill’s equations [11] The Buttrill’s equations are full governing equations of flight dynamics of flexible spinning missile, which are based on the following six assumptions: 1. The aircraft is idealized as a collection of lumped-mass elements, each being a finite rigid body, and each having an associated mass and moments of inertia. 2. The elastic restoring force resulting from displacement of any mass element is linear and proportional to that displacement. 3. The total rotational displacement of any lumped mass with respect to its undeformed orientation is small. That is to say, the change in orientation of a lumped mass element due to elastic deformation is assumed to be negligible, however, the change in orientation of a lumped mass element due to rigid-body motion is included. Therefore, the assumption is not only applicable to the aircraft, but also to the spinning missile with large rotational angular velocity. 4. Deformation is described by a linear sum of mode shapes multiplied by their time-dependent participation coefficients. 5. Gravity is constant over the aircraft. 6. The sea-level local Earth frame is assumed to be an inertial frame.

63 64 65 66

The six assumptions are easy to be satisfied and are appropriate for most of questions. And the equations in practical mean-axes reference frame [37] can be written as:

˙ = F + mg − m(ω × V) mV

67

˙ + h jk η j η¨ k = L − ω × Jω − ˙Jω − h jk η˙ j η˙ k − ω × h jk η j η˙ k Jω

68

(5)

˙ · h jk ηk = Q j − K jk ηk + 2ω · h jk η˙ k M jk η¨ k − ω 1





72

2

73

Here, the term such as x jk yk follows Einstein summation convention, which means indicial sum over k, that is to say

x jk yk = x j1 y 1 + x j1 y 1 + ... + x jN y N

(6) dx˙ dt

x˙ = is time rate of change of x and x¨ = is time rate of

change of x˙ ; m = i dmi is total mass of the aircraft; i indicates sum over the lumped mass; the subscript i indicates lumped mass id i; dmi is mass of lumped mass id i; V is velocity in body ˙ is acceleration in body frame; F is aerodynamic force; g frame; V is acceleration of gravity; ω is angular velocity in body frame; J is real-time two-order inertia tensor, and can be expressed as:

J = J0 + J j η j +

1 2

2 J jk η j ηk

J0 =

{(ri · ri )I − ri rTi }dmi

74 75 76 77 78 79 80 81 82 83 84 85 86

(7)

87 88

Here, the subscripts j and k indicate elastic modes j and k, respectively; J0 is two-order inertia tensor in undeformed condition, and is defined as



70 71

+ ωT J j + 2 J jk ηk ω

dx dt

69

89 90 91 92

(8)

i

93 94

Where, ri is position vector locating the center of the lumped mass i respect to the body frame when the body is in the undeformed reference condition; I is two-order unit tensor. J j is first partial derivative of the inertia tensor with respect to elastic mode j, and it is a two-order tensor defined as:

 J j = {(2ri · φi j )I − φi j rTi − ri φTi j }dmi

(9)

i

95 96 97 98 99 100 101 102 103

2 J jk is second partial derivative of the inertia tensor with respect

104

to elastic modes j and k, and it is a two-order tensor defined as:

105 106

2 J jk = A jk + Akj

(10)

2

 A jk = {(φi j · φik )I − φi j φTik }dmi

109

(11)

i



110 111 112

T Here A jk = Akj . η j is generalized displacement of mode j; η˙ j is generalized velocity of mode j; η¨ j is generalized acceleration of ˙ is angular acmode j; φi j is mode j shape at lumped mass i; ω celeration in body frame; h jk indicates energy coupling between rigid-body angular velocity and modal velocity in mode k due to deformation in mode j, and it is a vector defined as:

h jk =

107 108

Here  J jk =  Jkj . A jk is a two-order tensor defined as: 2

113 114 115 116 117 118 119

φi j × φik dmi

(12)

120 121

i

Here h jk = −hkj . L is total applied moment on the aircraft; ˙J is

time derivative of inertia tensor J; M jk = i φi j · φik dmi is generalized mass coupling between elastic mode j and mode k; K jk is generalized stiffness coupling between elastic mode j and mode k; Q j is generalized aerodynamic force of mode j. What should be mentioned here is mode shapes are all orthogonal to each other, and the matrices M and K are diagonal, and the relationship between the terms in the two matrices is

122 123 124 125 126 127 128 129 130 131

K jj = M jjω jω j

j = 1, 2, ..., N

(13)

132

JID:AESCTE AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.5 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

Here, ω j is undamped structural natural frequency of mode j without rotational motion, and N is the number of modes. Eq. (5) is flight dynamics equations of flexible aircraft, which not only consider the aerodynamic coupling effects between rigid body six degree of freedom motion and structural dynamics, but also consider the variation of inertia tensor and extra moment terms caused by structural vibration, and even more importantly, the extra force term caused by angular acceleration of rigid body, Coriolis term and centrifugal loading term caused by angular velocity of rigid body are considered, too. And it’s easy to find that there is no inertial coupling between translational equation and elastic equation because of the linearized mean-axes constraints [11], but there is inertial coupling between angular equation and elastic equation. In this paper, the investigation is focused on the aeroelasticity of flexible missile with rotational angular velocity, therefore, the Eq. (5) could be simplified to

18 19

˙ · h jk ηk = Q j − K jk ηk + 2ω · h jk η˙ k M jk η¨ k − ω

20

1

25 26 27

(14)

2

22 24



+ ωT J j + 2 J jk ηk ω

21 23



Eq. (14) is generalized dynamic aeroelastic equations for flexi˙ is zero, that is ble spinning missile. If the angular acceleration ω to say, if angular velocity is invariable, Eq. (14) could be further simplified to





29

1 M jk η¨ k = Q j − K jk ηk + 2ω · h jk η˙ k + ωT J j + 2 J jk ηk ω (15) 2

30

Eq. (14) could be also written as

28

˙ · h jk − M jk η¨ k − 2ω · h jk η˙ k + ( K jk − ω

33 34 35

˜ jk = M jk M

40 41 42 43 44 45

(16)

2

37 39

2

ωT 2 J jk ω)ηk

1

Now, let

38

1

= Q j + ωT J j ω

36

G˜ jk = −2ω · h jk

˙ · h jk − K˜ jk = K jk − ω Q˜ j = Q j +

1 2

1 2

(17)

ωT 2 J jk ω

T

ω J j ω

Thus, Eq. (16) could be written as

46 47 48 49 50 51 52 53 54 55 56 57 58

˜ jk η¨ k + G˜ jk η˙ k + K˜ jk ηk = Q˜ j M

(18)

Remove ∼ and Eq. (18) could be written in matrix form:

Mη¨ + Gη˙ + Kη = Q

(19)

Where, M, K, G are the equivalent generalized mass matrix, the equivalent generalized stiffness matrix and the equivalent generalized damping matrix respectively and Q is equivalent generalized force. Define state variable x = [η1 , η2 , · · · , η N , η˙ 1 , η˙ 2 , · · · , η˙ N ]T , and therefore the second-order ordinary differential equation could be written as

x˙ = F(x, t ) = Ax + BQ(x, t )

63 64 65 66

69 70 71 72 73 74 75 76 77 78 79 80

Fig. 3. Generalized displacements history at V cf = 0.330.

81

A hybrid linear multi-step (HLM) scheme [38] is applied to solve Eq. (20) in time domain. For the same missile, the structural natural frequencies are different between the cases with and without rotational motion, that is to say, the structural natural frequencies are changed due to rotational motion. The structural natural frequencies could be determined from the imaginary part of eigenvalues of the matrix A by numerical method in physical time marching process. Here, the natural frequencies represent the frequencies of structural free vibration viewed in body frame (not in inertial frame).

V cf =

bs ωα



Eq. (20) is generalized dynamic aeroelastic equation which is a first-order ordinary differential equation described in state space. Where



A=

0 I −M−1 K −M−1 G





, B=

0 M−1



83 84 85 86 87 88 89 90 91 92 94 95 96

The dynamic aeroelastic response of AGARD 445.6 weakened modal wing case [39] is chosen to verify the aeroelastic module. The incoming free flow Mach number is 0.901 and the angle of attack is 0 degree, and there is no sideslip angle. Euler equation is solved, and dual time stepping method is used to march physical time, and pseudo time is marched by Gauss-Seidel iteration, and spatial discretization method is Roe scheme, and the first four modes are considered during the simulation. Physical time step is 1/40 of the natural period of the highest mode, and the sub-iteration number is 200 in order to make sure the residual convergent, therefore the simulation precision can be ensured. Generalized displacement responses at different nondimensional velocity are depicted from Fig. 3 to Fig. 5. Here, the non-dimensional velocity V cf is defined as:

V

82

93

3. Verification case of numerical method

97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112

μ¯

(21)

Where V is stream velocity, bs is reference length, ωα is natural ¯ is mass ratio, circular frequency of first uncoupled torsion mode, μ and more details can be found in reference [39]. As shown in Fig. 4, when non-dimensional velocity V cf is about 0.342, generalized displacement responses are in equal amplitude oscillation, thus the non-dimensional flutter critical speed is 0.342, which agrees well with experimental result 0.37 [39]. Therefore, the code can give a reliable result for the fluid-structure coupling simulation, and can be used in the simulations related to aeroelasticity.

113 114 115 116 117 118 119 120 121 122 123 124 125 126

(20)

61 62

68

4. Computational results and analysis

59 60

67

3.1. AGARD 445.6 dynamic aeroelastic case

31 32

5

4.1. Geometrical model and simulation states

127 128

The spinning missile is shown as Fig. 6, the diameter of the missile body is 0.2 m, and the length of the missile is 4.74 m, and the slenderness ratio is 23.7. The sweep angle of head fins are 45 degree, and the sweep angle of tail fins are 35 degree. The

129 130 131 132

JID:AESCTE

6

1 2 3 4 5 6 7 8 9 10 11

AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.6 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

far boundary mesh, wall boundary mesh and space mesh slice are shown in Fig. 7. The number of the volume elements of the mesh is about 425 thousand, while the number of the volume nodes of the mesh is about 79 thousand. The incoming free flow Mach number is 2.5, and flight altitude is 8 km. Euler equations are chosen as flow governing equations, and AUSM+ scheme is used as spatial discretization method. Dual time stepping method is employed to march physical time, and pseudo time is marched by GaussSeidel iteration. The generalized dynamic aeroelastic equations of the flexible spinning missile is solved by second order precision of HLM scheme [38].

The structural dynamics analysis software ANSYS is used for modal analysis, and the first eight modes are chosen for further calculation and analysis. Fig. 8 shows structural natural frequencies

67 68 69 70 71 72 73 74 75 76 77

12

78

13

79

14

80

15

81

16

82

17

83 84

18

Fig. 5. Generalized displacements history at V cf = 0.350.

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92 93

27 28

Fig. 4. Generalized displacements history at V cf = 0.342.

Fig. 6. Geometric parameter.

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130 131

65 66

Fig. 7. Computational mesh.

132

JID:AESCTE AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.7 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

7

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112 113

47 48

Fig. 8. The First eight modes of structure.

50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

114 115

49

and mode shapes for the missile without rigid-body motion, and here the generalized displacements are all 1.0. Structural modes are calculated only for the missile body, and the head and tail fins are supposed as rigid bodies, and the mode shapes of the fins are interpolated from the missile body. Viewed in body frame, the first, third, fifth and seventh modes mainly vibrate in lateral direction, while the second, fourth, sixth, and eighth modes mainly vibrate in longitudinal direction. In fact, these modes can be grouped into four bending modes. The first and second modes are the first bending mode, the third and fourth modes are the second bending mode, and it is similar for the fifth-sixth and seventh-eighth modes. Initial conditions are very important when simulating dynamic aeroelastic response of flexible spinning missile. Initial flight angle of attack is 2.07 degree without fin deflection angle. The first initial generalized displacement is 1.0 (physical deformation is shown in Fig. 8a), and all other initial structural generalized displacements

are 0.0. The initial structural generalized velocities are all 0.0. The rotational speed is 31.3 rad/s (4.98 Hz), which is 1/4 of first structural natural frequency, and it is constant in the whole simulation process. Physical time step should be chosen according to flow field, rotational motion and structural natural frequencies. In this paper, physical time step is chosen as 2.0e-4 second, which is about a 1/29 of the structural natural period of the highest mode, and the spinning missile will rotate about 0.36 degree in each physical time step, therefore the simulation precision can be ensured. The aeroelastic simulation starts from a converged flow field calculated by the steady method, and the pseudo iteration number of CFD solver module in the aeroelastic simulation is 100 per physical time step in order to make sure pseudo iterations convergent. When solving the generalized dynamic aeroelastic equations depicted by Eq. (16), there are three terms h, J and 2 J that need to be evaluated beforehand. The terms h and J are both

116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.8 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

8

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79 80

14 15

Fig. 9. Generalized displacement for the first mode with rotational motion.

Fig. 11. Generalized displacement for the third mode with rotational motion.

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95 96

30 31

Fig. 12. Generalized displacement for the fourth mode with rotational motion.

32 33 34 35 36 37 38 39 40 41 42

97

Fig. 10. Generalized displacement for the second mode with rotational motion.

100 101 102 103 104 105 106 107 108 109

43 44

110

4.2. Results and analysis

111

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

98 99

three-order tensors, while the term 2 J is a four-order tensor. The three terms are very complicated and are related to mass distribution and mode shapes. In this paper, the three terms are evaluated by numerical method. The mesh used by the software ANSYS in modal analysis is used to discretize the missile, then a collection of lumped-mass elements are obtained, and the mode shapes at lumped-mass center can be obtained by interpolation method. Now, since the location and mode shapes of each lumped mass have already been obtained, the three terms can be evaluated according to the definitions.

The generalized displacements for all eight modes with rotational motion are indicated from Fig. 9 to Fig. 16. Overall, the generalized displacements of first two modes are far larger than those of others, which indicates the first two modes play a major role among all eight modes. The amplitude and period of the generalized displacements of the first and second modes are almost same. The phenomenon also happens on the third and fourth modes, fifth and sixth modes, and seventh and eighth modes. In initial stage, the amplitudes of the all eight modes are large, however, in final stage, the amplitudes are small. The amplitudes decrease due to aerodynamic damping. Fig. 9 and Fig. 10 indicate the generalized displacements of the first and second modes, respectively. The initial generalized displacement of the first mode (lateral direction) is 1.0, while the initial generalized displacement of the second mode (longitudinal direction) is 0.0. Although the initial generalized displacement of the two modes are different, the response histories are similar, which indicates the two modes are coupled with each other. The reason can be explained from the structural vibration equation depicted by Eq. (19). When there is no rotational angular velocity, the matrix M and K are diagonal, and the matrix G is zero,

112 113 114

Fig. 13. Generalized displacement for the fifth mode with rotational motion.

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131

Fig. 14. Generalized displacement for the sixth mode with rotational motion.

132

JID:AESCTE AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.9 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

9

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79 80

14 15

Fig. 15. Generalized displacement for the seventh mode with rotational motion.

16

Fig. 17. Shape of generalized damping matrix G (the symbol × denotes large values, while the symbol • denotes negligible values).

81 82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28 29 30 31

Fig. 16. Generalized displacement for the eighth mode with rotational motion.

32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

therefore, the vibration equation of each mode can be separately solved. Except for aerodynamic coupling among modes, there is no other coupling between modes. However, when rotational angular velocity exists, the equivalent generalized damping matrix G( G˜ jk = −2ω · h jk ) depicted by Fig. 17 is not zero, and it is a

antisymmetric matrix ( G˜ jk = −G˜ kj ), and this is to say, except for aerodynamic coupling among modes, there is other coupling between different modes, and the coupling between different modes is due to damping term caused by rotational motion. Numerical results show that: the terms G˜ 1,2 , G˜ 2,1 , G˜ 3,4 , G˜ 4,3 , G˜ 5,6 , G˜ 6,5 , G˜ 7,8 and G˜ 8,7 in matrix G are large, while other terms are negligible. Therefore, the coupling between the first and second modes are strong, and it’s similar for the third and fourth, fifth and sixth, and seventh and eighth modes. The generalized displacements of the first and second modes contain vibrations with different frequencies, and it could be seen more clearly in the detail figure of the final stage. In final stage, the period of the generalized displacements of the two modes is about 0.2 s, namely, the frequency of the two modes is about 5 Hz, which is almost equal to the frequency of the rotational motion (4.98 Hz). It indicates the vibration of the two modes is a forced vibration, and the exciting force is caused by rotational motion. Fig. 11 and Fig. 12 indicate the generalized displacements of the third and fourth modes, respectively. The initial generalized displacements of the two modes are both 0. However, the amplitudes of the two modes in initial stage are large, and then they gradually decrease, which is similar to the first and second modes. The reason why they are similar can be explained as follows. As Fig. 18 shown, except for the coupling effects between the third and fourth modes due to rotational motion, there are coupling effects between the first and third modes, and coupling effects between the second and fourth modes due to aerodynamic forces. In a word, rotational motion and aerodynamic forces couple the four modes together. Therefore, the third and fourth modes are affected

Fig. 18. Coupling effects between modes (solid line represents coupling effects due to rotational motion, dashed line represents coupling effects due to aerodynamic forces).

94 95 96 97

by the first and second mode, and they are similar. In final stage, the period of the generalized displacement of the third and fourth modes is about 0.2 s, which is almost the same as that of the first and second modes, and is almost equal to the frequency of the rotational motion. Fig. 13 and Fig. 14 indicate the generalized displacements of the fifth and sixth modes, respectively. And Fig. 15 and Fig. 16 indicate the generalized displacements of the seventh and eighth modes, respectively. These four modes are similar to the third and fourth modes, and therefore details are omitted here. The generalized displacements for all eight modes without rotational motion are indicated from Fig. 19 to Fig. 26. Overall, the first mode play a major role in structural vibration, while the generalized displacement of the second mode is small, which is different from the case with rotational motion. This is mainly because there is no rotational motion, thus the coupling effects between the first and second mode are little. The other lateral vibrations including the third, fifth and seventh modes are affected by the first mode due to aerodynamic coupling. There is longitudinal static aeroelastic bending due to the attack of angle, therefore, in final stage, the average amplitude of longitudinal vibrations including the second, fourth, sixth and eighth mode are not zero. Fourier analysis are made for the generalized displacements of the eight modes, and the results are indicated from Fig. 27 to Fig. 34. Fig. 27 indicates the distribution of amplitude of generalized displacement for first and second modes with rotational motion. There are three peak frequencies in Fig. 27, and the values are 4.98 Hz, 15.14 Hz and 25.02 Hz, respectively. The peak frequency 4.98 Hz is equal to the frequency of the rotational motion, and this peak frequency indicates the structural vibration is directly affected by rotational motion. The peak frequencies 15.14 Hz and 25.02 Hz are close to the structural natural frequency of the first and second modes (19.93 Hz), but not equal. However, as for the case without rotational motion depicted by Fig. 28, there is only one peak frequency and its value is equal to the structural natural frequency. Compared with the case without rotational mo-

98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

10

AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.10 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79 80

14 15

Fig. 19. Generalized displacement for the first mode without rotational motion.

Fig. 23. Generalized displacement for the fifth mode without rotational motion.

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96 97

31 32

Fig. 20. Generalized displacement for the second mode without rotational motion.

Fig. 24. Generalized displacement for the sixth mode without rotational motion.

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113 114

48 49

Fig. 21. Generalized displacement for the third mode without rotational motion.

Fig. 25. Generalized displacement for the seventh mode without rotational motion.

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130 131

65 66

Fig. 22. Generalized displacement for the fourth mode without rotational motion.

Fig. 26. Generalized displacement for the eighth mode without rotational motion.

132

JID:AESCTE AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.11 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

11

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79 80

14 15 16

Fig. 27. Distribution of amplitude of generalized displacement for first and second modes with rotational motion.

Fig. 29. Distribution of amplitude of generalized displacement for third and fourth modes with rotational motion.

81 82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32 33

Fig. 28. Distribution of amplitude of generalized displacement for first and second modes without rotational motion.

Fig. 30. Distribution of amplitude of generalized displacement for third and fourth modes without rotational motion.

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

99 100

34 35

98

tion, the peak frequencies 15.14 Hz and 25.02 Hz in the case with rotational motion may be the structural natural frequencies with rotational motion, and actually they are, and the reason will be explained later. Fig. 29 indicates the distribution of amplitude of generalized displacement of third and fourth modes with rotational motion. There are five peak frequencies, one is the frequency of rotational motion, and another two are the frequencies close to the structural natural frequency of the first and second modes, while the last two are the frequencies close to the structural natural frequency of the third and fourth modes. The five peak frequencies indicate there are coupling effects between first-second mode and third-fourth mode. However, as for the case without rotational motion depicted by Fig. 30, there are only two peak frequencies, and one is the frequency of the structural natural frequency of the first and second modes, and another is the frequency of the structural natural frequency of the third and fourth modes. The two peak frequencies indicate there are coupling effects between first mode and third mode, and coupling effects between second mode and fourth mode. The reason for coupling effects of the case with/without rotational motion are different. For the case without rotational motion, the coupling effects are due to aerodynamic forces. However, for the case with rotational motion, except for the coupling effects due to aerodynamic coupling, there are coupling effects resulting from rotational motion. Fig. 31 and Fig. 32 indicate the distribution of amplitude of generalized displacement of fifth and sixth modes with/without rotational motion, respectively. Fig. 33 and Fig. 34 indicate the distribution of amplitude of generalized displacement of seventh and eighth modes with/without rotational motion, respectively. These four modes are similar to the third and fourth modes, therefore more details are omitted here.

101 102 103 104 105 106 107 108 109 110 111 112 113 114

Fig. 31. ARABIC Distribution of amplitude of generalized displacement for fifth and sixth modes with rotational motion.

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130

Fig. 32. Distribution of amplitude of generalized displacement for fifth and sixth modes without rotational motion.

131 132

JID:AESCTE

12

AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.12 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

1

Table 1 BIC theoretical natural frequency of eight modes when rotational frequency is 4.98 Hz.

2 3 4

67 68 69 70

5

Theoretical natural frequency (Hz)

71

6

14.99 24.87 51.22 61.00 102.09 111.80 166.52 176.20

72

7 8 9 10 11 12

15

74 75 76 77 78 79

13 14

73

80

Fig. 33. Distribution of amplitude of generalized displacement for seventh and eighth modes with rotational motion.

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

Fig. 36. Theoretical structural natural frequency history for all eight modes without rotational motion.

28 30 32

Fig. 34. Distribution of amplitude of generalized displacement for seventh and eighth modes without rotational motion.

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Fig. 35. Theoretical structural natural frequency history for all eight modes with rotational motion.

49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

94 95

29 31

93

In order to figure out what the peak frequencies in Fourier analysis (from Fig. 27 to Fig. 34) mean, the structural natural frequencies for the case with and without rotational motion are calculated by numerical method. What should be mentioned here is that for the elastic vibration with N modes, the size of the matrices M, G and K are all N × N, while the size of the matrix A is 2N × 2N, and there are 2N eigenvalues for the matrix A in the form of N complex conjugate pairs. In fact, there should be no difference in structural natural frequencies between a longitudinal mode and a lateral mode that consist of a bending mode, therefore, it would be easy to deduce that each pair of complex conjugate eigenvalues represents two structural natural frequencies, and one is for the longitudinal mode, and another is for the lateral mode. Therefore, there are 2N structural natural frequencies for the vibration with N modes, that is to say, for each mode, there are two structural natural frequencies. For simplicity, only one structural natural frequency for each pair of complex

conjugate eigenvalues is shown here. Fig. 35 indicates theoretical structural natural frequency history with rotational motion for all eight modes, and they are changeless in the whole process, and more detail values are shown in Table 1, and Fig. 36 indicates those without rotational motion. Except for the frequency of rotational motion, the peak frequencies of the Fourier analysis (depicted in Fig. 27, Fig. 29, Fig. 31 and Fig. 33) are well agreed with structural natural frequencies, namely, the peak frequencies reflect the structural natural frequency. Furthermore, the reason why structural natural frequencies are changed due to rotational motion could be explained as follows. The structural natural frequency with rotational motion for each mode is different from that without rotational motion, because rotational motion will introduce extra terms in structural vibration equation depicted by Eq. (16). Compared with the structural vibration equation without rotational motion, there are four extra terms in structural vibration ˙ · h jk ηk , equation depicted by Eq. (16), and they are −2ω · h jk η˙ k , −ω − 12 ωT 2 J jk ωηk and 12 ωT J j ω, respectively. The term −2ω · h jk η˙ k is a Coriolis term, which represents the effects of rotational angular velocity on structural vibration, and it acts as a damping term ˙ · h jk ηk represents the in structural vibration equation. The term −ω effects of rotational angular acceleration on structural vibration, which acts as a stiffness term in structural vibration equation, and it is zero when rotational angular velocity is invariable. The term − 12 ωT 2 J jk ωηk is centrifugal loading, which represents the effects of rotational angular velocity on structural vibration, and it acts as a stiffness term in structural vibration equation. The last term 1 T ω J j ω is centrifugal loading too, however, it is a generalized 2 force term in structural vibration equation. The first three terms have effects on structural natural frequency. For rotational motion with invariable rotational angular velocity, in other word, when ˙ · h jk ηk is zero, Fig. 37 indicates theoretical structural the term −ω natural frequencies with rotational frequency. When there is no rotational motion, for each mode, the two structural natural frequencies are equal. However, when there is rotational motion, the two structural natural frequencies are not equal any more, and that is to say, there are two different structural natural frequencies for

96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.13 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Fig. 37. Theoretical structural natural frequency with rotational frequency.

15 16 17 18 19

each mode. The structural natural frequencies are changed due to rotational motion, and when rotational angular velocity increases, the difference between the two structural natural frequencies for each mode will increase too.

20 21

5. Concluding remarks

22 23 24 25 26 27 28 29 30 31

Based on rigid-motion mesh and RBF mesh morphing method, unsteady Euler equations and generalized dynamic aeroelastic equations of flexible spinning missile derived from Buttrill’ equation are coupled simultaneously to simulate the dynamic aeroelastic response history of the spinning missile with large slenderness ratio. The investigation is focused on the effects of rotational motion on dynamic aeroelasticity, and the structural dynamic responses for the missile with and without rotational motion are simulated, respectively. The following conclusions can be drawn:

32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

1. Rotational motion will directly change the generalized damping matrix in structural vibration equation which will couple the lateral modes and longitudinal modes together. When there is no rotational motion, the equivalent damping term is zero, therefore, the coupling are negligible between the lateral modes and longitudinal modes. However, when there is rotational motion, the equivalent damping term is not negligible, and the lateral modes and longitudinal modes are coupled. 2. Rotational motion will also change the stiffness term in structural vibration equation which has effects on structural natural frequencies. Under the effects of rotational motion, the generalized stiffness matrix and generalized damping matrix are both changed, and therefore the structural natural frequencies are changed, and for each mode there will be two different structural natural frequencies whose values are not equal. When rotational angular velocity increases, the difference between the two values of structural natural frequencies for each mode will increase too.

51 52

Declaration of competing interest

53 54

There is no conflict of interest.

55 56

Acknowledgement

57 58 59 60 61 62 63 64 65 66

This work could not be accomplished successfully without the financial support from National Natural Science Foundation of China under the grant No. 11732013. References [1] P. Reisenthel, D. Lesieutre, D. Nixon, Prediction of aeroelastic effects for seaskimming missiles with flow separation, in: 32nd Structures, Structural Dynamics, and Materials Conference, 1991.

13

[2] M. Karpel, S. Yaniv, D.S. Livshits, Integrated solution for computational static aeroelasticity of rockets, J. Spacecr. Rockets 35 (1998) 612–618. [3] D. Raveh, M. Karpel, S. Yaniv, Nonlinear design loads for maneuvering elastic aircraft, J. Aircr. 37 (2000) 313–318. [4] S. Chae, D. Hodges, Dynamics and aeroelastic analysis of missiles, in: 44th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2003. [5] E. Baskut, ¸ A. Akgül, Development of a coupling procedure for static aeroelastic analyses, Sci. Tech. Rev. 61 (2011) 39–48. [6] L.K. Abbas, D. Chen, X. Rui, Numerical calculation of effect of elastic deformation on aerodynamic characteristics of a rocket, Int. J. Aerosp. Eng. 2014 (2014). [7] D. Chen, L.K. Abbas, X. Rui, G. Wang, Aerodynamic and static aeroelastic computations of a slender rocket with all-movable canard surface, Proc. Inst. Mech. Eng., G J. Aerosp. Eng. 232 (2018) 1103–1119. [8] D. Tihomirov, D.E. Raveh, Nonlinear aerodynamic effects on static aeroelasticity of flexible missiles, in: AIAA Scitech 2019 Forum, 2019. [9] G. Reis, W. Sundberg, Calculated aeroelastic bending of a sounding rocket based on flight data, J. Spacecr. Rockets 4 (1967) 1489–1494. [10] M.R. Waszak, D.K. Schmidt, Flight dynamics of aeroelastic vehicles, J. Aircr. 25 (1988) 563–571. [11] C. Buttrill, P. Arbuckle, T. Zeiler, Nonlinear simulation of a flexible aircraft in maneuvering flight, in: Flight Simulation Technologies Conference, Monterey, CA, U.S.A., 1987. [12] T. Zeiler, C. Buttrill, Dynamic analysis of an unrestrained, rotating structure through nonlinear simulation, in: 29th Structures, Structural Dynamics and Materials Conference, 1988. [13] M.R. Waszak, C.S. Buttrill, D.K. Schmidt, Modeling and Model Simplification of Aeroelastic Vehicles: An Overview, 1992. [14] D. Platus, Aeroelastic stability of slender, spinning missiles, J. Guid. Control Dyn. 15 (1992) 144–151. [15] D.K. Schmidt, D.L. Raney, Modeling and simulation of flexible flight vehicles, J. Guid. Control Dyn. 24 (2001) 539–546. [16] N. Nguyen, I. Tuzcu, Flight dynamics of flexible aircraft with aeroelastic and inertial force interactions, in: AIAA Atmospheric Flight Mechanics Conference, 2009. [17] E.J. Oliveira, P. Gasbarri, I. Milagre da Fonseca, Flight dynamics numerical computation of a sounding rocket including elastic deformation model, in: AIAA Atmospheric Flight Mechanics Conference, Dallas, TX, USA, 2015. [18] Z. Shi, L. Zhao, Effects of aeroelasticity on coning motion of a spinning missile, Proc. Inst. Mech. Eng., G J. Aerosp. Eng. 230 (2016) 2581–2595. [19] T. Xu, J. Rong, D. Xiang, C. Pan, X. Yin, Dynamic modeling and stability analysis of a flexible spinning missile under thrust, Int. J. Mech. Sci. 119 (2016) 144–154. [20] D. Zhang, S. Tang, Q.-j. Zhu, R.-g. Wang, Analysis of dynamic characteristics of the rigid body/elastic body coupling of air-breathing hypersonic vehicles, Aerosp. Sci. Technol. 48 (2016) 328–341. [21] C. Xie, L. Yang, Y. Liu, C. Yang, Stability of very flexible aircraft with coupled nonlinear aeroelasticity and flight dynamics, Aircr. Eng. 55 (2017) 862–874. [22] S.A. Keyes, P. Seiler, D.K. Schmidt, Newtonian development of the mean-axis reference frame for flexible aircraft, Aircr. Eng. (2018) 1–6. [23] E. Blades, J. Newman, Aeroelastic effects of spinning missiles, in: 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2007. [24] A. Schütte, G. Einarsson, A. Raichle, B. Schöning, W. Mönnich, M. Orlt, J. Neumann, J. Arnold, T. Forkert, Numerical simulation of maneuvering aircraft by aerodynamic, flight mechanics and structural mechanics coupling, Aircr. Eng. 46 (2009) 53–64. [25] R.H. Hua, C.X. Zhao, Z.Y. Ye, Y.W. Jiang, Effect of elastic deformation on the trajectory of aerial separation, Aerosp. Sci. Technol. 45 (2015) 128–139. [26] R.H. Hua, Z.Y. Ye, J. Wu, Effect of elastic deformation on flight dynamics of projectiles with large slenderness ratio, Aerosp. Sci. Technol. 71 (2017) 347–359. [27] M.J. Li, X.T. Rui, L.K. Abbas, Elastic dynamic effects on the trajectory of a flexible launch vehicle, J. Spacecr. Rockets 52 (2015) 1–17. [28] J. Yin, J. Lei, X. Wu, T. Lu, Effect of elastic deformation on the aerodynamic characteristics of a high-speed spinning projectile, Aerosp. Sci. Technol. 45 (2015) 254–264. [29] J. Yin, J. Lei, X. Wu, T. Lu, Aerodynamic characteristics of a spinning projectile with elastic deformation, Aerosp. Sci. Technol. 51 (2016) 181–191. [30] L. Yang, Z.Y. Ye, The interference aerodynamics caused by the wing elasticity during store separation, Acta Astronaut. 121 (2016) 116–129. [31] L. Yang, Z.Y. Ye, J. Wu, The influence of the elastic vibration of the carrier to the aerodynamics of the external store in air-launch-to-orbit process, Acta Astronaut. 128 (2016) 440–454. [32] T. Lu, X. Wu, J. Lei, J. Yin, Numerical study on the aerodynamic coupling effects of spinning and coning motions for a finned vehicle, Aerosp. Sci. Technol. 77 (2018) 399–408. [33] R. Hua, X. Yuan, Z. Tang, Z. Ye, Study on flight dynamics of flexible projectiles based on closed-loop feedback control, Aerosp. Sci. Technol. (2019). [34] T.C. Rendall, C.B. Allen, Efficient mesh motion using radial basis functions with data reduction algorithms, J. Comput. Phys. 228 (2009) 6231–6249.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

14

1 2 3 4 5 6 7

AID:105384 /FLA

[m5G; v1.261; Prn:9/09/2019; 11:51] P.14 (1-14)

H. Li, Z.Y. Ye / Aerospace Science and Technology ••• (••••) ••••••

[35] T. Rendall, C. Allen, Parallel efficient mesh motion using radial basis functions with application to multi-bladed rotors, Int. J. Numer. Methods Eng. 81 (2010) 89–105. [36] G. Wang, H.H. Mian, Z.-Y. Ye, J.-D. Lee, Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions, AIAA J. 53 (2014) 1016–1025. [37] M.R. Waszak, D.K. Schmidt, On the flight dynamics of aeroelastic vehicles, in: Astrodynamics Conference, 1986.

[38] 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 (2007) 253–262. [39] E.C. Yates Jr, AGARD Standard Aeroelastic Configurations for Dynamic Response. Candidate Configuration I.-wing, 1987.

67 68 69 70 71 72 73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65

131

66

132