A method for nonlinear aeroelasticity trim and stability analysis of very flexible aircraft based on co-rotational theory

Journal of Fluids and Structures 62 (2016) 209–229

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

A method for nonlinear aeroelasticity trim and stability analysis of very flexible aircraft based on co-rotational theory Wei Wang a,1, Xiaoping Zhu b,2, Zhou Zhou a,n,3, Jingbo Duan c,4 a b c

College of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China Science and Technology of UAV Laboratory, Northwestern Polytechnical University, Xi’an 710065, China Ordance Engineering College, Shijiazhuang 050003, China

a r t i c l e i n f o

abstract

Article history: Received 15 June 2015 Accepted 21 January 2016

Very flexible aircraft with high aspect ratio wings subjected to aerodynamic loads undergoes large deformation, which will lead to distinct changes on the mass distributions, stiffness characteristics and aerodynamic characteristics of the complete aircraft. The aeroelasticity and flight dynamics of such aircrafts are nonlinear and the linear elastic theory model cannot be used. A new method is developed for the analysis of nonlinear aeroelasticity and flight dynamics of very flexible aircraft through combining the corotational beam theory with the modified ONERA dynamic stall model. Based on a form of co-rotational technique which is external to the element, a spatial two-node beam element, which depicts the geometrically nonlinear dynamic characteristics of the flexible wing, is developed. Both tangential stiffness matrix and mass matrix of the beam element are formulated to establish the nonlinear dynamic equations. In addition, the modified ONERA dynamic stall model is adapted to evaluate the unsteady nonlinear aerodynamic loading of the very flexible wing. Using the present method, the nonlinear aeroelastic response, trim and stability characteristics of a very flexible aircraft are predicted in this paper. The obtained results show a good agreement to the literature, which indicates that the present method is accurate and efficient. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Very flexible aircraft Nonlinear aeroelasticity Trim Flight dynamics Co-rotational approach Geometrically nonlinearity

1. Introduction HIgh-altitude and ultra-long-endurance solar-powered unmanned aerial vehicles (UAVs) which can fill the gap between low-orbit satellites and conventional fuel-powered UAVs, are designed to perform a very large spread of low cost continuous coverage of a geographic area missions, such as meteorological monitoring, coastal surveillance, tele-communication relay, and so on (Flittie and Curtin, 1998; Cestino, 2006; Romeo et al., 2004; Noll et al., 2004; Min et al., 2010; Wei et al., 2014). As the energy available from solar cells is very limited, the flexible structural design strategy is usually introduced to minimize the structural weight and to maximize aerodynamic efficiency of high aspect ratio wings. For example, the conventionally n

Corresponding author. E-mail addresses: [email protected] (W. Wang), [email protected] (X. Zhu), [email protected] (Z. Zhou), [email protected] (J. Duan). 1 Department of Aeronautics A 212 Hang Kong Building. 2 UAV Research Institute, A 106 Hang Kong Building. 3 Department of Aeronautics A 211 Hang Kong Building. 4 A 212 Hang Kong Building. http://dx.doi.org/10.1016/j.jfluidstructs.2016.01.009 0889-9746/& 2016 Elsevier Ltd. All rights reserved.

210

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

Nomenclature Aρ Aa D e fl f f mass g I Ixyz Jρ J trim k L M Kl KT KTσ N pl p R Sk;T S

mass per unit length angular acceleration in body attached frame drag unit vector of coordinate system local nodal generalized forces global nodal generalized forces inertial force vector gravitational acceleration unit matrix global reference coordinate system mass moment of inertial tensor criterion function linear momentum per unit length lift of aerodynamics moment of aerodynamics linear stiffness matrix tangential stiffness matrix geometric stiffness matrix interpolation function local nodal displacements global nodal displacements rotation matrix Jacobian matrix skew-symmetric matrix

Sm T u ul u_ V δW w π γ Δ ϑ_ ϑ€ α,β,r Γ ρ

first moment area of a section transformation matrix aircraft body attached reference coordinate system axial deformation linear velocity flow velocity virtual work angular velocity in body attached frame angular momentum per unit of length pseudo-vector increment angular velocity in inertial frame angular acceleration in inertial frame integration parameters vortex strength density of air

Subscript i,j m e l

nodal number of beam element ‘mean configuration’ coordinate system element coordinate system local

powered high altitude long endurance (HALE) UAV of Global Hawk has an aspect ratio of 25 and a surface density of the wing structure of 53 kg/m2. While, the aspect ratio of the solar-powered UAV of Helios is 31 and the surface density of the wing structure is only 3.2 kg/m2. High aspect ratio wings with extremely low weight will undergo large deformation during normal flight, which will lead to distinct changes on the structural stiffness characteristics, the mass and the aerodynamic loading distribution of the aircraft. In that case, nonlinearity coupling arising from large wing deflection may play an important role on the vehicle aeroelasticity and flight dynamics, as it was sadly illustrated by the mishap of NASA’s Helios. Due to the complexity of the problem, these nonlinear problems are still hotspots of research. Therefore, precise analysis methods for nonlinear aeroelasticity and flight dynamics of high aspect ratio flexible aircraft need urgently to be developed to deal with its nonlinear aeroelastic, trim, flight stability and other characteristics (Noll et al., 2004). The elastic deformation of a very flexible wing under aerodynamic loads belongs to typical geometrically nonlinear problems with large displacement but small strains. Currently, there are three main kinds of approaches to analyze the flexible structure with the geometric nonlinearity, which are all based on finite element method (FEM) but with every one choosing a different selection of independent degrees of freedom including displacements, strains, and combined velocities and internal forces of the intrinsic formulation (Palacios et al., 2010). The intrinsic beam model was used to describe the geometrically large deformation of flexible wing (Hodges, 2003). Combining intrinsic beam equations with finite-state induced flow model (Peters, 1985), Patil and Dowell (2001), Patil and Hodges (2006), Chang (2008), and Patil et al. (1999) developed a nonlinear aeroelastic trim and stability analysis method that was suitable to HALE aircraft, and their computation codes are referred as NATASHA. In literature (Zhang and Xiang, 2008), Zhang researched the nonlinear aeroelastic response of a very flexible wing based on intrinsic beam element and ONERA dynamic stall model (McAlister and Lambert, Petot). Further, in the development of nonlinear aeroelastic and flight dynamics analysis tools, Cesnik and Su (2005), Su and Cesnik (2006), Shearer and Cesnik (2007) Shearer and Cesnik (2008), and Su (2008) developed the University of Michigan’s Nonlinear Aeroelastic Simulation Toolbox (UM/NAST) which combines nonlinear strain-based beam model, unsteady aerodynamics with simplified stall model and a six degree of freedom flight dynamic formulation. Other research works include: Shams et al. (2012) and Zhao and Ren (2011) established a geometrically nonlinear model of flexible wing structure through using nonlinear general flexible Euler–Bernoulli beam equations and multi rigid body theory, separately. The advantage of relatively low order of intrinsic beam theory or strain-based theory can provide significant reduction in computation cost of solving dynamics problems. However, its independent variants are not straightforward and lead to this method being less versatile than the displacement based element method (Patil and Dowell, 2001). Presently, there are three main methods for the solving of displacement based beam element, which are Total Lagrangian (TL), Update Lagrangian (UL) and Co-rotational (CR) formulation, respectively. Based on the linearization of UL methods, Yang et al. (2012) (Chuan and Chao, 2011) studied the nonlinear aeroelastic and flight load characteristics of a very flexible

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

211

Fig. 1. Sketch of motion and coordinate systems of a UAV.

aircraft. In literature (Wei et al., 2015), a linearization of CR methods was used to explore the nonlinear aeroelastic stabilities of a flexible aircraft with the layout similar to Helios, and some measures of improving adverse effects caused by large deformation were also proposed. Summarily, the classical TL and UL method usually improve the calculation accuracy through increasing load steps, thus leading to significant reduction in computational efficiency. Unlike the TL or UL method which derives tangential stiffness matrix from direct nonlinear relations between displacement and stain, the tangential stiffness matrix of CR method is obtained by decomposing the structural motions into rotation and translation of the rigid body and purely elastic deformation of the element in local coordinate system. So CR method has a preferably computational accuracy and efficiency, and is widely used in the solving of geometrically nonlinear problems of large displacements but small strains (Belytschko and Schwer, 1977; Crisfield et al., 1997; Le et al., 2012). In the present paper, a new nonlinear aeroelastic trim and flight stability analysis method is established and its iterative solution process is presented by coupling a CR finite element model and the unsteady nonlinear aerodynamic strip theory based on the ONERA dynamic stall model. Then, taking a flexible wing as an example, the dynamic response of the nonlinear aeroelasticity is studied by taking into account the large deformation and dynamic stall nonlinear effects. Secondly, the longitudinal nonlinear trim and flight load characteristics of a very flexible aircraft with the layout similar to Helios under different payloads are studied by setting constant trim altitude and velocity as prerequisites. The convergence of the proposed nonlinear trim algorithm is also discussed in this paper. Finally, the flight dynamics of the previously mentioned aircraft is explored with the nonlinear aeroelastic effects being considered.

2. Geometrically nonlinear structural modeling 2.1. Static equations based on CR theory in view of geometrically nonlinearity High-aspect-ratio flexible wing is generally described as a geometrically nonlinear beam. The equilibrium equation is firstly established in the local coordinate system when CR beam elements are chosen for the flexible structural model. The local coordinate system can be obtained through using the global coordinate system Ixyz, the nodal coordinate system ui and the ‘mean configuration’ coordinate system um, as shown in Fig. 1. Based on the configuration after the rigid motion and elastic deformation of a wing structure, the pseudo-vector that transforms global coordinate system to node coordinate system of the element and the current position of the node described in the global coordinate system can be obtained. We can obtain the nodal coordinate system by the following formulation: ui ¼ Ri ðθi ÞI xyz :

ð1Þ

For a two node spatial beam element with each node having six degrees of freedom, the nodal coordinate system ui and uj of the element can be obtained by Eq. (1). The pseudo-vectors (θix, θiy, θiz)T and (θjx, θjy, θjz)T corresponding to rotation matrix Ri and Rj respectively, are evaluated from the initial conditions or from the current solution of the problem. Let γ be the pseudo-vector associate with the relation from Ri to Rj and the ‘mean configuration’ coordinate system um can be defined with the relation:
 ð2Þ um ¼ R γ=2 ui ; 2

where Rðγ=2Þ ¼ I þ

sin ðλÞ 1  cos ðλÞ Sðγ=2ÞSðγ=2Þ, λ Sðγ=2Þ þ λ2

T 1=2

λ ¼ ððγ=2Þðγ=2Þ Þ

0 6 , Sðγ=2Þ ¼ 4 γ 3 =2  γ 2 =2

 γ 3 =2 0 γ 1 =2

γ 2 =2

3

γ 1 =2 7 5, and I repre0

sents the unit matrix with three orders. The element coordinate system ue is defined in such a way that its first vector ee1

212

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

lies along the element between nodes i and nodes j. It is possible to define ue as 8 e ¼ ðd2  d1 Þ=ln > < e1 ee2 ¼ em2  ððeTm2 ee1 Þ=2Þðee1 þem1 Þ ; > : e ¼ e  ððeT e Þ=2Þðe þe Þ e3

m3

m3 e1

e1

ð3Þ

m1

where d1 and d2 are the current positions of nodes i and j separately, and ln ¼((d2 d1)T( d2 d1))1/2. All of the nodal coordinate system, ‘mean configuration’ coordinate system and element coordinate system are functions of the structural deformation which will be updated during the process of solving. The element coordinate system should be distinguished between the deformed structure and the undeformed one for considering the deflection effects. Therefore, CR method can fully indicate the large displacement but small strain characteristics that is different to the conventional linear finite element algorithm. Incremental method is usually selected for the CR finite element analysis. The tangential stiffness matrix described in the global coordinate system can be obtained through the rotation matrix and its differential formulation from element coordinate system to global coordinate system, which is a sticking point of establishing the incremental static equilibrium equations. Combining Eqs. (1) and (3), the local displacements are obtained as 8 ul ¼ ln l0 > > > > > > > > > > 2 sin θl1;1 ¼  eTi3 ee2 þeTi2 ee3 > > > > > > > > > > > 2 sin θl1;2 ¼  eTi1 ee3 þeTi3 ee1 > > > > > > > < 2 sin θl1;3 ¼  eTi2 ee1 þeTi1 ee2 ; ð4Þ > > > > > > > > 2 sin θl2;1 ¼  eTj3 ee2 þeTj2 ee3 > > > > > > > > T T > > > 2 sin θl2;2 ¼  ej1 ee3 þej3 ee1 > > > > > > > > T T > : 2 sin θl2;3 ¼  ej2 ue1 þ ej1 ee2 where ln is the current length and l0 is the initial length of the element, and eij, j¼1, 2, 3 are the components of the rotation matrix Ri. The local nodal displacements vector pl can be written as  
 T T pTl ¼ dl1 ; θTl1 ; dl2 ; θTl2 ¼ dl1;1 ; dl1;2 ; dl1;3 ; θl1;1 ; θl1;2 ; θl1;3 ; dl2;1 ; dl2;2 ; dl2;3 ; θl2;1 ; θl2;2 ; θl2;3 ; The global nodal displacements vector p can be written as  
 T T pT ¼ d1 ; θT1 ; d2 ; θT2 ¼ d1;1 ; d1;2 ; d1;3 ; θ1;1 ; θ1;2 ; θ1;3 ; d2;1 ; d2;2 ; d2;3 ; θ2;1 ; θ2;2 ; θ2;3 ; where dl1 is the local linear displacements, θl1 is the local angular displacements. The subscript e indicates the element, m indicates the ‘mean configuration’, (1,1) indicates the components of displacements in the global coordinates system, (l1,1) indicates the components of displacements in the local coordinates system. Since the origin of element coordinate system always consists of node i and the vector ee1 of ue pass through node j, then dl1;1 ¼ dl1;2 ¼ dl1;3 ¼ dl2;2 ¼ dl2;3 ¼ 0: The incremental transformation matrix T is obtained by differentiating Eq. (4) as δpl ¼ Tδp;

ð5Þ




 where ðδpl ÞT ¼ 0; 0; 0; δθl1;1 ; δθl1;2 ; δθl1;3 ; δul ; 0; 0; δθl2;1 ; δθl2;2 ; δθl2;3 . The local nodal generalized forces fl can be written as h f l ¼ 0 0 0 M l1;1 M l1;2 M l1;3 Nl 0 0 M l2;1

M l2;2

M l2;3

The global nodal generalized forces f can be written as  f ¼ F 11 F 12 F 13 M11 M 12 M 13 F 21 F 22 F 23

M 21

M 22

iT

:

M 23

T

:

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

213

By equating the virtual work in the local and global systems, then T

T

f l δpl ¼ f δp:

ð6Þ

With the assumption that the deformation in local system still satisfies the linear displacement–strain relationships, then f ¼ T T f l ¼ T T K l pl :

ð7Þ

The tangential stiffness matrix can be expressed by differentiation of Eq. (7) δf ¼ T T δf l þ δT T f l ¼ T T K l Tδp þK Tσ ðf l ÞδP ¼ K T δp:

ð8Þ

According to the definition of the tangential stiffness matrix in the global frame, the formulation can be written as K T ¼ T T K l T þ K Tσ ðf l Þ;

ð9Þ

T

where T K l T is generally called material stiffness matrix and K Tσ is the geometric stiffness matrix is a function of ql . Because of the employment of element coordinate system and nodal coordinate system, the large displacement matrix that must be solved in the traditional TL method is not needed for the CR method. Thus, the calculation is extremely simplified. Based on the derived structural tangential stiffness matrix in Eq. (9), the global tangential stiffness matrix can be easily assembled. Then, the incremental equilibrium equation for structural statics analysis is given by K T Δp ¼ ΔF:

ð10Þ

The above equation is usually solved by Newton–Raphson or modified Newton–Raphson iterative algorithm with the incremental load or displacement control method, and more details can be found in Crisfield et al. (1997). 2.2. Dynamics equations Though the deformation of a flexible wing is large, it is reasonable that the area and shape of the wing cross sections may remain unchanged and the angular velocity for any point in that section is the same. Then the linear momentum and angular momentum per unit length are defined as Z k ¼ ðu_ þ ϑ_  ξÞdm ZS ð11Þ π ¼ ξ  ðu_ þ ϑ_  ξÞdm; S

where Z ZS ZS ZS S

dm ¼ Aρ ξdm ¼ S m _ ðξ  uÞdm ¼ S m  u_ _ ξ  ϑ_  ξdm ¼ I ρ ϑ;

ð12Þ

where Aρ is the mass per unit length and Sm is the first moment area of the cross section. The vector of u_ ¼ ½u_ 1 ; u_ 2 ; u_ 3 T is the  T linear velocity at the shear center of the cross-section and ϑ_ ¼ ϑ_ 1 ; ϑ_ 2 ; ϑ_ 3 is the angular velocity. The mass moment of 2 3 J1 þ J2 0 0 6 J1 0 7 inertial tensor is J ρ ¼ ρ4 0 5 and I ρ is defined as 0 0 J2 T ð13Þ I ρ ¼ UJ ρ U ; where U is the rotations from local coordinates system to global coordinates system. Then, the inertial forces and moments of any cross section of an element can be written as ! ! Aρ u€ þ ϑ€  US m d k s ¼ f mass ¼ ; dt π US m  u€ þ UJ ρ U T ϑ€ þ ϑ_  UJ ρ U T ϑ_

ð14Þ

where ϑ_  UJ ρ U T ϑ_ is gyroscopic term which does not appear in 2-D problems. Shear center accelerations of the cross section and the section angular velocities and angular accelerations can be obtained by interpolating the linear accelerations, angular velocities and angular accelerations of the node of the element. We take into consideration beam elements with two nodes and therefore we can obtain u€ ¼ N1 d€ i þ N2 d€ j ϑ_ ¼ N1 ϑ_ i þ N2 ϑ_ j

214

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

ϑ€ ¼ N1 ϑ€ i þ N2 ϑ€ j :

ð15Þ

Similarly, the virtual nodal displacement of a section can be obtained by interpolation and the virtual work done by inertial force on the virtual displacement can be expressed as Z l s ðδuT ; δϑT Þf mass dx: ð16Þ δW ¼ 0

Inserting Eqs. (14) and (15) into Eq. (16) and integrating through the entire beam element, one can obtain 0 1 N1 I 0 B C Z l N1 I C s B 0 T T Cf ðδdi ; δϑi T ; δdj ; δϑj T ÞB δW ¼ mass dx: B N2 I 0 C 0 @ A 0 N2 I

ð17Þ

According to the principle of virtual work, the virtual work done by the equivalent inertia forces is obtained as T

T

T

δW ¼ f mass ðδdi ; δϑi T ; δdj ; δϑj T ÞT ; 0

ð18Þ

1

N1 I Z lB 0 B B f mass ¼ B 0 @ N2 I 0

0 C N1 I C s Cf mass dx: 0 C A N2 I

ð19Þ

After inserting Eq. (14) into Eq. (19) and dropping the nodal virtual displacements, the inertial forces f mass can be derived as

0

1

N1 Aρ u€  N 1 USðS m ÞU T ϑ€

B Z l B N1 USðS m ÞU T u€ þ N 1 UðJ ρ Aa þ SðwÞJ ρ wÞ B B f mass ¼ B 0 B N2 Aρ u€  N 2 USðS m ÞU T ϑ€ @ N2 USðS m ÞU T u€ þ N 2 UðJ ρ Aa þ SðwÞJ ρ wÞ 2

_ where w ¼ U ϑ, T

¼ Aa ¼ dw dt

0 € SðwÞ ¼ 6 U ϑ, 4 w3 w2

ð20Þ

2 0 7 w1 5, SðSm Þ ¼ 6 4 Sm1 0  Sm2

w3

T

C C C Cdx; C C A w2

0 w1

3

 Sm3 0 Sm1

Sm2

3

Sm1 7 5. Aa is the angular accelera0

tions and w is the angular velocity in the element attached frame, respectively. By introducing the assumption that U is constant along the whole element, U can be replaced by Ue. So, the inertial forces vector can be written as f mass ¼ ½U e f mass;in þ f mass;o ; 0

I B0 B where ½U e  ¼ B @0

0 Ue 0

I

0

0

0

0

0 0

ð21Þ 0

1

1 0 0 B 3Sðw ÞJ w þSðw ÞJ w þ Sðw ÞJ w þSðw ÞJ w C 0 C B 1 ρ 1 1 ρ 2 2 ρ 1 2 ρ 2C C l B C ¼ M mass p€~ þ 12 C, f B C 0 0 A mass;in @ A Sðw1 ÞJ ρ w1 þ Sðw1 ÞJ ρ w2 þ Sðw2 ÞJ ρ w1 þ 3Sðw2 ÞJ ρ w2 Ue

Ue lB B 0 f mass;o ¼ B 6@ 0

0 Ue

0 0

0

Ue

0

0

0

10 0 0 B 2SðS m Þ 0 C CB CB 0 AB @ 0 SðS m Þ Ue

2SðS m Þ

0

0

SðS m Þ

 SðS m Þ

0

0

2SðS m Þ

 SðS m Þ

10

UT CB e 0 CB 0 CB B  2SðS m Þ C A@ 0 0 0

0

0

I

0

0

U Te

0

0

0

1

C 0C Cp~€ ; 0C A I

where f mass;o is the inertial force vector that represents the coupling of bending and torsion. It will become zero if it is assumed that the mass center is chosen to be located at the shear center of the cross-section. 0 1 Aρ I 0 2Aρ I 0 B 0 Iρ C 2I ρ C lB 0 C; M mass ¼ B B 0 2Aρ I 0 C 6@ Aρ I A 0 2I ρ 0 Iρ  T   T p~€ T ¼ d€ 1 ; ATa;1 ; d€ 2 ; ATa;2 Þ ¼ d€ 1;1 ; d€ 1;2 ; d€ 1;3 ; Aa;11 ; Aa;12 ; Aa;13 ; d€ 2;1 ; d€ 2;2 ; d€ 2;3 ; Aa;21 ; Aa;22 ; Aa;23 : It should be noted that the linear velocities are defined in the global coordinate system whereas the angular ones are given in the element coordinate frame.

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

215

According to the inertial formulation described in Eq. (21), the dynamic equilibrium equations with the damping ignored can be written as F e f mass f i ¼ 0;

ð22Þ

where F e is the external load vector at time t, f mass is the inertial forces vector and f i is the internal force vector. Obviously, Eq. (22) includes the inertial forces that related to accelerations and velocities compared with the pure statics equilibrium equations, and more details can be found Crisfield et al. (1997). The CR theory has the advantage of directly providing physical degrees of freedom for nonlinear aeroelastic problems. It should be noted that the mass matrix and tangential stiffness matrix need to be regularly updated in the solution.

3. Nonlinear unsteady aerodynamic modeling Dynamic stall refers to a strong nonlinear unsteady stall phenomenon that the flow field around the airfoil surface separated and then reattached after the effective angle of attack of oscillating lift surface exceeds its critical angle of attack. The local angle of attack of a very flexible aircraft may be very large during normal flight, which will lead to separation between the airflow and suction surface. This is the typical dynamic stall phenomenon that extensively influences the lift and moment of the airfoil section. In the present paper, the modified ONEAR (EDLin) dynamic stall model is used (Laxman and Venkatesan, 2007). The formulation of lift and moment that acts at the quarter-chord point is given as i 1 h _ ~ _ L ¼ ρS sbW ð23Þ 0 þ kbW 1 þ VΓ 1 þ VΓ 2 ; 2   2 V _ V Γ€ 1 þB1 Γ 1 þ B2 Γ 1 ¼ Γ 1;r ; b b "  #   2  V _ V V 2 V _ W0 ; Γ 2 þr Γ€ 2 þa Γ2 ¼  r VΔC L jW 0 =V þ E b b b b

ð24Þ

ð25Þ

n h io
2
2

 π _ 0 þ A2 V σ W € 0 þA1 σ W € 1 . s ¼ π þ 5π ð1  M 2 Þ0:285  1 _ 1 þ A1 ∂C zL W Γ 1;r ¼ A3 Vb ∂C∂αzL W 0 þ A3 σ Vb W 1 þ A2 Vb ∂C∂αzL W 1 180, ∂α b  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi π π 2π ffi π , λ ¼ 0:17  0:13M 1 , α~ ¼ 0:53 þ0:25ð 1 M 21 1Þ, d ¼ d1 jΔC L j, k~ ¼ 1  M 21 1 180 , σ ¼ pffiffiffiffiffiffiffiffiffiffiffi 2 þ 1:96π 1  M 21 180 pffiffiffi a ¼ a0 þa1 ΔC 2L , r ¼ r 0 þr 1 ΔC 2L , E ¼ E1 ΔC 2L . Γ 1 is the aerodynamic state in unstalled region in lift equation. Γ 2 is the aerodynamic state in stalled region in lift equation. ρ is the density of air and V is the velocity of inflow. M 1 is the mach number of inflow. S is the area of airfoil section and b is the semi-chord. Besides, A1 ¼0.5, A2 ¼0.393, A3 ¼0.0439425, B1 ¼0.551 and B2 ¼0.0439075 (Laxman and

_ þαÞ and W ¼ bα. Venkatesan, 2007). W 0 ¼ Vðh=V _ h_ represent the heaving velocity at elastic axis. ΔC L W 0 =V is the difference 1 between the linear static aerodynamic coefficient extrapolated to the stalled region to actual static aerodynamic coefficient _ þ α. of lift, measured at an effective angle of attack W 0 =V ¼ h=V Eqs. (24) and (25) can be simplified as # "
  "
2 #" _ #  Γ1 1 0 Γ€ 1 Γ 1;r B2 Vb B1 Vb and þ ¼ Γ1 0 0 1 Γ_ 1 1 0 where



1 0

0 1

" € # "
V  Γ2 a b þ Γ_ 2 1

#
 2 #" _ # "
2
 Γ2 r Vb ½r Vb VΔC L jW 0 =V þ E Vb W_ 0  ¼ : Γ2 0 0

The unsteady moment on the airfoil is given as i

1 ~ h 2 _ 0 þ σ m VW 1 þ sm bW _ 1 þ VΓ m2 ; V C mL W 0 =V þðσ~ m þ dm ÞbW M ¼ ρS2b 2 " 

  2 

V _ V V 2 V _ Γ m2 þ r m W0 : Γ€ m2 þ am Γ m2 ¼  r m VΔC m

W 0 =V þ Em b b b b Eq. (27) is written in a compact matrix form as i# # "

  "
2 #" _ # " h
V 2 _ 0 Γ m2 1 0 Γ€ m2  r m b VΔC m jW 0 =V þ Em Vb W am Vb r m Vb ; þ ¼ Γ m2 0 1 Γ_ m2 1 0 0 σ~ m ¼ 

i π  π πh 3π  1 þ 1:4M 21 ; dm ¼ σ 1m jΔC z j; sm ¼   1:26 1:53arctan½15ðM1  0:7Þ ; 4 180 16 180

ð26Þ ð27Þ

216

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

Fig. 2. Sketch of the aerodynamic loads acting at the airfoil.

 σ 0m ¼

i π  πh pffiffiffiffiffiffi 3π  1:26 1:53arctan½15ðM 1  0:7Þ  1 þ 1:4M 21 ; am ¼ am0 þ am1 ΔC 2z ; r m ¼ r m0 þ r m1 ΔC 2z ; 16 2 180

Em ¼ Em1 ΔC 2z ; σ m ¼ σ 0m þσ 1m jΔC z jo ;

ΔC m W 0 =V is the difference between the linear static aerodynamic coefficient extrapolated to the stalled region to actual _ þ α. C j static aerodynamic coefficient of moment, measured at an effective angle of attack W 0 =V ¼ h=V mL W 0 =V is the static _ moment coefficient in linear region measured at an effective angle of attack W 0 =V ¼ h=V þα. The differential equation of the modified ONERA dynamic stall model is solved by the fourth-order four-steps Runge– Kutta formulation as 8 Γ ¼ Γ tn þ Δt > > 6 ðK 1 þ 2K 2 þ 2K 3 þK 4 Þ > tn þ 1 > > > > > > > > K ¼ f ðt n ; Γ tn Þ > > 1 > > > > < K 2 ¼ f ðt n þ 12Δt; Γ tn þ 12ΔtK 1 Þ ð28Þ : > > > > > > > > > K 3 ¼ f ðt n þ 12Δt; Γ tn þ 12ΔtK 2 Þ > > > > > > > > : K 4 ¼ f ðt n þ Δt; Γ t þ ΔtK 3 Þ n

The effective angle of attack for the calculation of aerodynamic loading of a trim case can be indicated by the angle between the chord line and local air flow direction as shown in Fig. 2. The induced flow velocity V T and the effective angle of attack α are obtained according to the velocity projection in the element coordinate system: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V T ¼ V 2ey þ V 2ez ; ð29Þ
 α ¼ arctan V ey =V ez : Then, the projection of the aerodynamic loading in the element coordinate system can be written as 2 3 2 32 3 aero 0 0 0 0 M þf y U ea 6 7 6 6 7 aero cos α sin α 7 ¼4 fe ¼41 M aero 0 5: 54 L 5; e 0  sin α cos α D 0 When control surfaces are arranged on the wing sections, the lift L, drag D and moment M can be expressed as

L ¼ L δ ¼ 0 þ 12ρV 2 SC L;δ δ

1 D ¼ D

δ ¼ 0 þ ρV 2 SC D;δ δ 2

1

M ¼ M δ ¼ 0 þ ρV 2 ScC M;δ δ: 2

ð30Þ

ð31Þ

ð32Þ

4. Gravitational model The force due to gravitational effect is described in the inertial coordinate system as G

f g ¼ mgI z ; where g is the gravitational acceleration.

ð33Þ

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

217

Then, the distributed gravity model in the element coordinate system can be obtained by transforming the inertial coordinate system: G

f g ¼ ReI mgI z ;

ð34Þ

where ReI is the rotation matrix between the inertial coordinate system and the element coordinate system. Generally, the mass center of the wing section does not coincide with its elastic axis, which will lead to a moment caused by gravity. We can obtain the moment through projecting the gravity in inertial coordinate system to the element frame. G

M Gg ¼ mgSðξe ÞReI f g ; 2

0 6 where Sðξe Þ ¼ 4 ξez  ξey

 ξez 0 ξex

ð35Þ ξey

3

 ξex 7 5. 0

5. Propulsive model The propulsive systems are generally installed before the leading edge of the wing, as shown in Fig. 3. yp is the distance between the thrust line and the elastic axis, and φp is the angle between the thrust line and the chord line (i.e., the installation angle). The expression of thrust and its counter-torque in the element frame is given as 2

p fe

3 0 6 T p sin φ 7 ¼4 p 5; T p cos φp

2

3 T p yp 6 7 Mpe ¼ 4 M p sin φp 5:  Mp cos φp

ð36Þ

6. Nonlinear aeroelastic equations Combination of the previous structural and aerodynamics models provide a full description of the nonlinear aeroelastics of flexible wing. By assembling the nodal aerodynamic loadings, the unsteady aerodynamics vector Q e in global coordinate system can be obtained, which replaces the external force vector F e . Then, the formulation of the nonlinear aeroelastic dynamic response can be rewritten as Q e  f mass  f i ¼ 0:

ð37Þ

The solving of the nonlinear aeroelastic dynamic response by Eq. (37) is different to the pure structural dynamics analysis. The reason is that the external loads are known at tn þ 1 for the structural dynamics analysis, however, the unsteady aerodynamics Q e is related to the structural dynamic response which need to be computed previously. By estimating the unsteady aerodynamics at tn þ 1, we can obtain the updated aerodynamics according to the structural dynamic response results. Based on the equidistant Newton forward interpolation algorithm, we can predict the unsteady aerodynamics Q e;t n þ 1 by the extrapolation polynomial: Q e;t n þ 1 ¼ 3Q e;tn  3Q e;tn  1 þQ e;tn  2 :

ð38Þ

The loosely coupling calculation algorithm is used for solving the developed nonlinear aeroelastic system, and the procedure can be expressed as Fig. 4.

Fig. 3. Propulsive system.

218

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

Fig. 4. Procedures of solving nonlinear aeroelastic system by loosely coupled algorithms.

7. Calculation of nonlinear trim Choosing the ground axis I as the inertial reference coordinate system and the mean axis with the origin at the center of gravity as the body attached reference axial system, we can obtain the balance equation of the deformed aircraft as F x mGg g sin θ ¼ 0 F y þ mGg g sin φ cos θ ¼ 0 ; F z þmGg g cos φ cos θ ¼ 0

Mφ ¼ 0 Mθ ¼ 0 :

Mψ ¼ 0

ð39Þ

The nonlinear trim analysis of very flexible aircraft considering the geometrically nonlinear static aeroelastic effects needs the iterative formulations, since the structural stiffness, aerodynamic loading and its direction, the mass center of the whole aircraft and the space location of distributed propellers are associated with the wing deformation and must be updated along with the structural motions. The imbalance relations are formulated as 2 3 F x  mGg g sin θ 6 7 6 F y þmG g sin φ cos θ 7 g 6 7 6 7 6 F z þ mGg g cos φ cos θ 7 7: F trim ¼ 6 ð40Þ 6 7 Mφ 6 7 6 7 6 7 Mθ 4 5 Mψ The criterion of the above is defined as J trim ¼ F Ttrim F trim :

ð41Þ

For trimming performed of 1 g level flight with constant altitude and velocity, a simple numerical Newton–Raphson method is used to evaluate the Jacobian tangent matrix. The trim parameters of thrust Tp, elevator deflection angle δe and angle of attack α are written as a vector S¼(Tp, δe, α)T and the Jacobian tangent matrix can be obtained as  ∂F : ð42Þ Sk;T ¼ ∂S

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

219

The Jacobian tangent matrix in Eq. (42) can be calculated numerically through finite difference algorithms based on finite increment around the current trim cases. Having known F trim;k that represents the trim residual error of the k-th iteration, we can obtain the incremental trim parameters by the following formulation: ΔS ¼  ðSk;T Þ  1 F trim;k :

ð43Þ

Then, Sk þ 1 ¼ Sk þ ΔS:

ð44Þ

In summary, the nonlinear trim analyze system includes two modules: (1) static aeroelastic solving module that considers the geometrically nonlinearity; (2) trim iterative module that is similar to that of a rigid aircraft. The entire procedure is outlined in Fig. 5.

8. Numerical examples 8.1. Deformation of a very flexible cantilever beam under concentrated moment at the tip In order to validate the method suggested here, a very flexible cantilever beam with a concentrated moment at the tip is presented, as shown in Fig. 6. The length of beam is l ¼20 m. The physical properties of the beam are EIz ¼10 N m2 and EA¼10 N. For comparison, the analytical solutions to this problem with the load conditions M¼0.2π, 0.4π, 0.6π, 0.8π and π is an arc with a length of 20 m and central angle θ¼0.4π, 0.8π, 1.2π, 1.6π and 2π, respectively. With the present method, 20 elements are used for the beam modeling. The deformations obtained by linear and nonlinear structural model are shown in Figs. 7 and 8, respectively. It can be seen that, when the linear structural model is used, there is no displacement in x direction and the vertical displacements in y direction increase linearly which cannot describe the rigid rotational component of geometrically nonlinear deflection. So the result obtained from the linear structural model is obviously different from the nonlinear analytical solutions. When the nonlinear model is used, very good

Fig. 5. Procedures of nonlinear trim.

y M

x 20m Fig. 6. Cantilever beam with a concentrated tip moment.

220

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

Fig. 7. Linear deformation of a cantilever beam under concentrate tip moment.

Fig. 8. Nonlinear deformation of a cantilever beam under concentrate tip moment.

Table 1 Flexible wing model data (Patil and Dowell, 2001). Parameter

Value

Parameter

Value

Half-span (m) Chord (m) Mass per unit length (kg/m) Moment of inertial (kg m) Span-wise elastic axis

16.0 1.0 0.75 0.1 50% chord

Center of gravity Bending rigidity (N m2) Torsional rigidity (N m2) Chordwise bending rigidity (N m2) Density of air (kg/m3)

50% chord 2  104 1  104 4  106 0.0889

agreement between the theory results and present nonlinear results is observed with the maximum error no more than 3%. Besides, the sub-iteration will converge after 3 or 4 iterations. In conclusion, the above results show that the beam model established by CR theory is accurate and efficient. 8.2. Nonlinear aeroelastic dynamic response analysis of a very flexible wing In this section, the nonlinear aeroelastic dynamic response of a very flexible wing is investigated to validate the effectiveness of the developed method. The involved geometrical parameters, structural stiffness parameters, mass parameters and atmosphere parameters are listed in Table 1. The time responses of the wing tip obtained by the present method are compared with the existing results in literatures Patil and Dowell (2001) and Shams et al. (2012) (as shown in Fig. 9). Obviously, the results show a very close correlation among each other. It should be noted that the results in literature (Shams et al., 2012) is gained by the Wagner function in the Duhamel integration which is a linear aerodynamic model and cannot consider the dynamic stall effects. The

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

y/(m)

4.0 3.8 3.6 3.4 3.2 3.0 2.8 2.6 2.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0

221

Reference 10 Reference 21 Present work

0

2

4

6

8

10

12

14

16

18

20

t/(s) Fig. 9. Nonlinear wingtip vertical displacement vs. time in V ¼28 m/s, y0 ¼4.0 m.

0.4

Tip vertical velocity /(m/s)

5

0.0

-0.2

3 2 1 0 -1 -2 -3 -4 -5

-0.4 -8

-6

-4

-2

0

2

4

6

0.8

8

1.2

1.6

Tip angular velocity/(rad/s)

2.0

2.4

2.8

0.3 0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 0.5

3.2

3.6

Tip vertical displacement /(m)

0.4

Tip twist /(rad)

Tip twist /(rad)

0.2

4

1.0

1.5

2.0

2.5

3.0

3.5

Tip vertical displacement /(m) Fig. 10. Various phase-plane plots in V ¼ 28 m/s, y0 ¼4.0 m.

4.0

4.5

4.0

4.4

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

y/(m)

222

6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5

V=29m/s V=30m/s V=31m/s V=32m/s

0

5

10

15

20

25

t/(s) Fig. 11. Tip vertical displacement responses at various flow velocity in y0 ¼ 4.0 m.

2.0

Effective angle of attack/(rad)

1.8 1.6 1.4 1.2 1.0 0.8

V=29m/s V=30m/s V=31m/s V=32m/s

0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

t/(s) Fig. 12. Effective angle of attack of wingtip at various flow velocity in y0 ¼ 4.0 m.

phase-plane of the limit cycle oscillation at a speed of 28.0 m/s is plotted in Fig. 10. It can be observed that the proposed method can satisfactorily predict the limit cycle oscillation of a very flexible wing. Further, the nonlinear dynamic response characteristics of the very flexible wing at the velocity greater than the critical speed are studied with the initial wing tip displacement y0 ¼4.0 m, as shown in Figs. 11 and 12. It can be seen that, for the flow velocity slightly greater than 28 m/s, the wing tip displacement response of the very flexible wing is still periodic. However, when the flow velocity reaches 32 m/s, the oscillations of wing tip displacement lose periodicity and become chaotic, which is similar to the phenomenon obtained in literature (Patil and Dowell, 2001). What’s more, the amplitude of wing tip vertical displacement and twist angle get increased as the flow speed increases and the limit cycle oscillations become more complex. Besides, at the flow velocity v¼28.0 m/s, the oscillation magnitude of each period is the same. But, as the flow velocity increases, the oscillation will come back to the initial state through a large magnitude period and a small magnitude one. Also, the dynamic stall effects become more significant as the flow velocity increases. The nonlinear dynamic response characteristics of the very flexible wing with the velocity lower than the critical speed are studied still with the initial the wing tip displacement y0 ¼4.0 m, as shown in Figs. 13 and 14. It is viewed that, 25 s is needed for the wing moving to its equilibrium state at v¼27 m/s and 15 s is needed at v ¼26 m/s. It is indicated that the oscillation magnitude of the wing becomes smaller with the decrease of the flow velocity and the same to the oscillation convergence time of the flexible wing. 8.3. Nonlinear trim and flight dynamics analysis of a very flexible aircraft As the third example, the nonlinear trim and flight dynamics for a very flexible flying-wing similar to Helios are researched. The configuration and geometry of the flexible aircraft referred in the literature (Patil and Hodges, 2006) are

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

223

4.5 4.0

V=28m/s V=27m/s V=26m/s

3.5 3.0

y/m

2.5 2.0 1.5 1.0 0.5 0.0 -0.5 -5

0

5

10

15

20 t/s

25

30

35

40

Fig. 13. Vertical displacement responses history of the wing tip.

V=28m/s V=27m/s V=26m/s

0.3 0.2

z/m

0.1 0.0 -0.1 -0.2 -0.3 -5

0

5

10

15

20 t/s

25

30

35

40

Fig. 14. Chord-wise displacement responses history of the wing tip.

Node1

Node11

Pro_1

A

Node21

Node31

Node41

Pro_2

Pro_3

Pro_4

C

Node51

Node61

Pro_5

B

Fig. 15. Geometry of the aircraft (Patil and Hodges, 2006).

shown in Fig. 15. The wingspan is 72.8 m, the constant chord length is 2.44 m and the reference area of the wing is 177.6 m2. A dihedral angle of 10° is on the five-sixth of the wingspan and the flap spreading all over the wing is located on the wing trailing edge. There are five propulsive units and three pods denoted by A, B and C arranged in the straight segment (see Fig. 15). The aerodynamic properties and structural parameters for the cross section of the wing and the pods parameters are

224

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

given in Table 2. The payload is assumed to be a concentrate mass hanging below the elastic axis of the wing with a distance of 0.914 m. Both Pod A and B have a constant mass of 22.68 kg, and the mass of Pod C is assumed to range from 27.216 kg to 254 kg.

Table 2 Parameters of the aircraft (Patil and Hodges, 2006). Parameter

Value

Structural and aerodynamic parameters of the wing cross section Elastic axis 25%chord Bending rigidity (N m2) 1.032  106 Torsional rigidity (N m2) 1.65  105 Mass moment of inertial about y axis (kg m) 0.69 Mass moment of inertial about x axis (kg m) 4.146 Cmα (1/rad) 0.0 CD0 0.01 Cmδ (1/rad)  0.25 Pod parameters CLα 1/rad

5.0

Parameter

Value

Center of gravity Chord wise bending rigidity (N m2) Mass per unit length (kg/m) Mass moment of inertial about z axis (kg m) CLα (1/rad) CLδ (1/rad) Cm0

25%chord 1.24  107 8.93 3.455 2π 1 0.025

CD0

0.02

50

30

20

Rigid (present result) Nonlinear (present result) Linear (present result) Nonlinear(reference11)

10

0 0

20

40

60

80

100

120

140

160

180

200

220

240

Payload /(kg) Fig. 16. Trim thrust variation with payload mass.

Rigid (present result) Nonlinear (present result) Linear (present result) Nonlinear(reference 11)

6

Flap deflection in degrees

Thrust /(N)

40

5 4 3 2 1 0 -1 0

20 40 60 80 100 120 140 160 180 200 220 240

Payload /(kg) Fig. 17. Flap deflection variation with payload mass.

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

5.2

Flight angle of attack in degrees

5.0 Rigid ( present result ) Nonlinear ( present result ) Linear (present result) Nonlinear(reference 11)

4.8 4.6 4.4 4.2 4.0 3.8 3.6 3.4 3.2 3.0 0

20

40

60

80 100 120 140 160 180

200 220 240

Payload /(kg) Fig. 18. Flight angle of attack variation with payload mass.

37.470

Nonlinear Linear

Thrust /(N)

37.469 37.468 37.467 37.466 37.465 37.464 0

10

20

30

40

Iteration times Fig. 19. The convergence history of thrusts.

Flap deflection in degrees

4 Nonlinear Linear

3 2 1 0 -1 0

10

20

30

Iteration times Fig. 20. The convergence history of flap deflection.

40

225

226

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

Flight angle of attack in degrees

5.3 5.2 5.1 5.0 4.9 4.8 4.7 4.6 4.5 4.4 4.3 4.2 4.1 4.0

Nonlinear Linear

0

5

10

15

20

25

30

35

40

45

Iteration times Fig. 21. The convergence history of flight angle of attack.

Nonlinear trim Inertia Structural dynamics analysis

Aerodynamic derivatives

Root locus calculation Fig. 22. Procedures of nonlinear flight dynamics of very flexible aircraft.

0.6 0.5 0.4

reference 17 reference 11 present

0.3

imaginary

0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.15

-0.10

-0.05

0.00

0.05

0.10

0.15

Real Fig. 23. Phugoid modes of the flexible aircraft. Note: The arrows indicate the payload increasing direction.

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

227

8.3.1. Nonlinear trim analysis of the aircraft With the assumption that the deflected angle of all flap-like control surfaces and the thrust of all the propulsive units are the same, the trim of the aircraft is studied. Figs. 16–18 show the variation of trim parameters, including the thrust, the flaplike control surface along the trailing edge and the flight angle of attack, along with the payload at the center pod at the velocity of 12.19 m/s on the sea-level. It can be seen that the trimmed thrusts obtained respectively by rigid, linear and nonlinear model are close to each other. That is mainly attributed to the fact that a constant CD0 is assumed for the 2-D aerodynamic model, which does not changes with the local angle of attack. The trimmed deflecting angles of the control surfaces are researched with a series of payload weights. Fig. 17 plotted the results derived by rigid, linear and nonlinear model, respectively. As expected, there is no significant difference for the trimmed deflecting angles when the payload mass is light. However, the differences of the trimmed deflecting angles obtained by three models gradually become significant with the increasing of payload mass. In particular, the different results illustrate that both the rigid and linear models are invalid in solving the trim properties of a very flexible aircraft with large deformation. This can be attributed to the decrease in flap-like effectiveness of the outboard section of the wing with increasing bending deformation. Besides, the fact that the present nonlinear results agreed well with the literature (Patil and Hodges, 2006) show that the present method is accurate. The trimmed angles of attack for the aircraft are investigated with a series of payload weights, as shown in Fig. 18. It can be seen that, with the increase of the payload weight, the trimmed angle of attack with linear structural model is less than that of rigid model among all the payload regions. The reason is that the twist deformation arose by the nose-up moment caused by aerodynamic loading. Secondly, the trimmed angle of attack with linear structural model is almost the same to that with the nonlinear model under a light payload. When the weight of the payload is heavier than 180.0 kg, the trimmed angle of attack with linear structural model will be exceeded by that with the nonlinear model. The reason is that the invalidity of linear model in describing the rotation effects will lead to obvious error in predicting of wing bending deformation, which dominantly reduces the lift coefficient of the wing through changing the direction of local lift of the outboard wing section. Further, the convergence of the algorithm is discussed to evaluate the accuracy and stability of the proposed method, as shown in Figs. 19–21. Take the full load condition for example and select the convergence residuals error J trim o 10  8 for each iteration.

Fig. 24. Phugoid modes of the rigid aircraft. Note: The arrows indicate the payload increasing direction. Table 3 Short-period mode roots variation payload mass. Payload (kg)

Real

Imaginary

Real

Imaginary

27.216 47.216 67.216 87.216 107.216 127.216 147.216 167.216 187.216 207.216 227.216 254

 10.98034144  10.63985627  10.32080248  10.00631419  9.692862168  9.380468612  9.066068735  8.745925462  8.414171136  8.061555099  7.671265787  7.015088307

0 0 0 0 0 0 0 0 0 0 0 0

 4.4735089  3.5269651  2.7329262  2.1521724  1.7827723  1.5221019  1.3382898  1.2068161  1.1098216  1.0347906  0.972584  0.8600181

0 0 0 0 0 0 0 0 0 0 0 0

228

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

It can be observed that, 12 iterations will be generally required for a nonlinear trimming, the calculation cost of which is less than that of the CFD technology. It should be pointed that the iteration times of nonlinear model is almost the same to the linear one and no more iteration will be introduced by the geometrical nonlinearity. Though the iteration times of nonlinear model is almost the same to the linear one, the computation cost may be higher for its introducing of subiteration. In conclusion, the trim results of thrust, flap deflection and flight angle of attack show good agreements with that in literatures, which indicates that the proposed algorithm is accurate and valid and can depict the geometric nonlinear effects in the large deformation of a very flexible aircraft. 8.3.2. Flight dynamics of the very flexible aircraft Finally, the flight dynamics of the very flexible aircraft previously given in Section 8.3.1 are studied. Based on the trim state, the coupled framework to describe the flight dynamic model is given in Fig. 22, and the flight stabilities are plotted in Figs. 23 and 24. In the figures, there is some difference of roots between the present results and the reference ones when the payload is light. It may be brought by the calculation error of aerodynamic derivatives or inertial properties. But, both of the trends of the results are the same as the payload increases. Particularly, phugoid roots have good agreement with the literature when the payload is heavy. So, the proposed flight dynamics analysis algorithms have a good accuracy and can be used as a tool to investigate the stability of flight dynamics of a very flexible aircraft. In Figs. 23 and 24, with the increasing of payload, the phugoid mode of the elastic aircraft will be unstable while the mode of rigid aircraft is stable even in heavy payload state. So the results obtained with the rigid aircraft assumption for the very flexible aircraft may be unsuitable. The short-period mode roots with the increasing of payload mass are shown in Table 3. With the increasing of payload mass, the short-period mode of the aircraft is always stable. So, it can be concluded that the short-period mode is not very sensitive to the changing of inertial properties and aerodynamics due to elastic deformation.

9. Conclusion This paper provides a method for the analysis of nonlinear aeroelastic response, trim and stability characteristics of a very flexible aircraft. The whole calculation model includes geometrically nonlinear structural model, nonlinear unsteady aerodynamic model, gravitational model, propulsive model, nonlinear aeroelastic equations and calculation of nonlinear trim. The geometrically nonlinear structural model is based on the CR beam theory and is composed with static equations and dynamics equations which are solved with the HHT α algorithm. The nonlinear unsteady aerodynamic model is based on the ONERA dynamic stall model to considering the typical dynamic stall phenomenon that extensively influences the lift and moment of the airfoil section. The loosely coupling calculation algorithm is used for solving the nonlinear aeroelastic system and the iterative formulations for the nonlinear trim analysis of the nonlinear aeroelastic system, respectively. With the present method, the nonlinear deformation and aeroelastic response of the very flexible wing are studied firstly. Then, the trim for a representative example case of a Helios-like aircraft are predicted and the flight dynamic stability of the aircraft was analyzed at the trimmed condition with different payloads. The good agreements between the present solutions and the literature ones indicate that the present method is accurate and efficient.

References Belytschko, T., Schwer, L., 1977. Large displacement, transient analysis of space frames. International Journal for Numerical Methods in Engineering 11, 65–84. Cesnik C.E.S., Su W., 2005. Nonlinear Aeroelastic Modeling and Analysis of Fully Flexible Aircraft, AIAA-2005-2169. Cestino, E., 2006. Design of Very Long-Endurance Solar Powered UAV (Ph.D. dissertation). Department of Aerospace Engineering, Politecnico di Torino. Chang, C.S., Hodges, D.H., Patil, M.J., 2008. Flight dynamics of highly flexible aircraft. Journal of Aircraft 45 (2), 538–545. Changchuan, Xie, Chao, Yang, 2011. Linearization methods of nonlinear aeroelastic stability for complete aircraft with high-aspect-ratio wings. Science China Technological Sciences 54, 403–411. Crisfield, M.A., Galvanetto, U., Jelenic,́ G., 1997. Dynamics of 3-D co-rotational beams. Computational Mechanics 20, 507–519. Flittie, Kirk, Curtin, Bob, 1998. Pathfinder Solar-Powered Aircraft Flight Performance, AIAA Paper No.1998-4446. Hodges, D.H., 2003. Geometrically exact, intrinsic theory for dynamics of curved and twisted anisotropic beams. AIAA Journal 41 (6), 1131–1137. Laxman, V., Venkatesan, C., 2007. Chaotic response of an airfoil due to aeroelastic coupling and dynamic stall. AIAA Journal 45 (1), 271–280. Le, T.N., Battini, J.M., Hjiaj, M., 2012. Dynamics of 3d beam elements in a corotational context: a comparative study of established and new formulation. Finite Elements in Analysis and Design 61, 97–111. McAlister K.W., Lambert O., Petot D., 1984. Application of the ONERA Model of Dynamic. NASA, Technical Paper, No. 2399. Min, Chang, Zhou, Zhou, Zhicheng, Zheng, 2010. Flight principles of solar-powered airplane and sensitivity analysis of its conceptual parameters. Journal of Northwestern Polytechnical University 28 (5), 792–796. Noll, T.E., Brown, J.M., Perez-Davis, M.E., Ishmael, S.D., 2004. Investigation of the Helios Prototype Aircraft Mishap. Technical Report. NASA. Palacios, Rafael, Murua, Joseba, Cook, Robert, 2010. Structural and aerodynamic models in nonlinear flight dynamics of very flexible aircraft. AIAA Journal 48 (11), 2648–2659. Patil, M.J., Dowell, E.H., 2001. Limit cycle oscillation in high-aspect-ratio wings. Journal of Fluids and Structures 15, 107–132. Patil, M.J., Hodges, D.H., 2006. Flight dynamics of highly flexible wings. Journal of Aircraft 43 (6), 1790–1799.

W. Wang et al. / Journal of Fluids and Structures 62 (2016) 209–229

229

Patil, M.J., Hodges, D.H., Cesnik, C.E.S., 1999. Nonlinear Aeroelasticity and Flight Dynamics of High-Altitude Long-Endurance Aircraft, AIAA-1999-1470. Peters, D.A., 1985. Toward a unified lift model for use in rotor blade stability analysis. Journal of American Helicopter Society 30 (3), 32–42. Romeo, G., Frulla, G., Cestino, E., Corsino, G., 2004. HELIPLAT: design, aerodynamic, structural analysis of long-endurance solar-powered stratospheric platform. Journal of Aircraft 41 (6), 1505–1520. Shams, S., Sadr, M.H., Haddadpour, H., 2012. An efficient method for nonlinear aeroelasticy of slender wings. Nonlinear Dynamics 67, 659–681. Shearer, C.M., Cesnik, C.E.S., 2007. Nonlinear flight dynamics of very flexible aircraft. Journal of Aircraft 44 (5), 1528–1545. Shearer, C.M., Cesnik, C.E.S., 2008. Trajectory control for very flexible aircraft. Journal of Guidance, Control and Dynamics 31 (2), 340–357. Su, W., 2008. Coupled Nonlinear Aeroelasticity and Flight Dynamics of Fully Flexible Aircraft (Ph.D. dissertation). The university of Michigan, Ann Arbor, MI. Su, W., Cesnik, C.E.S., 2006. Dynamic Response of Highly Flexible Flying Wings, AIAA-2006-1636. Wei, Wang, Zhou, Zhou, Xiaoping, Zhu, Rui, Wang, 2014. Static aeroelastic characteristics analysis of a very flexible solar powered UAV with geometrical nonlinear effect considered. Journal of Northwestern Polytechnical University 32 (4), 499–504. Wei, Wang, Zhou, Zhou, Xiaoping, Zhu, Rui, Wang, 2015. Exploring aeroelastic stability of very flexible solar powered UAV with geometrically large deformation. Journal of Northwestern Polytechnical University 33 (1), 1–8. Yang, C., Wang, L.B., Xie, C.C., et al., 2012. Aeroelastic trim and flight load analysis of flexible aircraft with large deformations. Science China Technological Sciences 55, 2700–2711. Zhang, J., Xiang, J.W., 2008. Nonlinear aeroelastic response of high-aspect-ratio flexible wings. Chinese Journal of Aeronautics 22 (4), 355–363. Zhao, Z.J., Ren, G.X., 2011. Multibody dynamic approach of flight dynamics and nonlinear aeroelasticity of flexible aircraft. AIAA Journal 49 (1), 41–54.