Geometrical nonlinear aeroelastic stability analysis of a composite high-aspect-ratio wing 1

A composite high-aspect-ratio wing of a high-altitude long-endurance (HALE) aircraft was modeled with FEM by MSC/NASTRAN, and the nonlinear static equilibrium state is calculated under design load with follower force effect, but without load redistribution. Assuming the little vibration amplitude of the wing around the static equilibrium state, the system is linearized and the natural frequencies and mode shapes of the deformed structure are obtained. Planar doublet lattice method is used to calculate unsteady aerodynamics in frequency domain ignoring the bending effect of the deflected wing. And then, the aeroelastic stability analysis of the system under a given load condition is successively carried out. Comparing with the linear results, the nonlinear displacement of the wing tip is higher. The results indicate that the critical nonlinear flutter is of the flap/chordwise bending type because of the chordwise bending having quite a large torsion component, with low critical speed and slowly growing damping, which dose not appear in the linear analysis. Furthermore, it is shown that the variation of the nonlinear flutter speed depends on the scale of the load and on the chordwise bending frequency. The research work indicates that, for the very flexible HALE aircraft, the nonlinear aeroelastic stability is very important, and should be considered in the design progress. Using present FEM software as the structure solver (e.g. MSC/NASTRAN), and the unsteady aerodynamic code, the nonlinear aeroelastic stability margin of a complex system other than a simple beam model can be determined.


Introduction
In recent decades, the High-Altitude Long-Endurance (HALE) aircraft has been given great interest, and the aeroelastic stability of high-aspect-ratio wings has been a key subject for many people [1][2][3].Dowell et al. [4], pointed out clearly that very high aspect ratio wings with both structural and aerodynamic nonlinearities is one of the important research subjects."A statically nonlinear, but dynamic linear, . . .Then the question of the dynamic stability of the statically nonlinear fluid-structural (aeroelastic) system may be addressed by a linear static dynamic perturbation analysis about this nonlinear static equilibrium."For a complex structure, the Finite Element Method (FEM) gives out so many degrees of freedom that post instability motions, such as Limit Cycle Oscillation (LCO), are hard and inefficient to calculate.For this reason, the dynamic linearization is the practicable engineering approach to determine the critical aeroelastic instability conditions of a wing or a complete aircraft.
In the last ten years, the beamlike wing -as the simplest model of the high-aspect-ratio wing -has been used to study the mechanism of its aeroelastic behavior.Patil and Hodges [5] studied the geometrical nonlinearity effects on the static and dynamic aeroelastic behavior using the dynamic linearization method.In their model, the chordwise bending frequency was higher than the torsion frequency, the geometric nonlinearities induced the coupling of these two modes where torsion frequency was decreasing, and the corresponding flutter speed decreased as the load magnitude was growing.And the non-planar effects of aerodynamics were negligible even for large deformation.Tang and Dowell [6] tested a beamlike wing in the wind tunnel, and compared with the theoretical results, in which the dynamic linearized calculation gave out a quite precise prediction of the flutter speed.Recently, the dynamic linearization method started to be applied to the complex structure in engineering.Strong et al. [7], contractedly studied a complex pre-stressed wing to demonstrate the differences comparing with linear analysis.
The composite material can be beneficial to the structure, but the stiffness characteristics of a composite structure are different from the metallic ones and the geometrical nonlinear effects are different too.For the traditional metallic high-aspect-ratio wing, the bending and torsion are elastically uncoupled, and the stiffness of chordwise bending is inherently quite large.So the chordwise bending mode has a little contribution to the linear flutter.As for the deformed wing, the geometrical nonlinear effect is mainly showed in the change of the torsion mode leading to the changing of the flutter speed.For composite structures, the bending and torsion stiffness coupling is usually introduced to satisfy the linear aeroelastic constrains and light weight, and then the chordwise bending stiffness would sometimes be much smaller than that of torsion.However, for the high-aspect-ratio wing with large deformation during flight, the chordwise bending and torsion coupling is strengthened due to the geometrical nonlinearities.The low frequency chordwise bending mode with large torsion component can lead to severe nonlinear aeroelastic problems.
In the present work, the aeroelastic stability of a real imbedded composite high-aspect-ratio wing with large symmetrical deformation is studied.The static nonlinear equilibrium is calculated firstly, and then the system is dynamic linearized to determine the aeroelastic unstable condition.The flutter of the wing is caused by flap bending and chordwise bending which has a large torsion component due to geometrical nonlinearities.

Structural and aeroelastic theory
For the high-aspect-ratio wing of HALE aircraft, the structure has quite large deformation in the climbing maneuver process, even during the straight flight.The static and dynamic analysis of the structure should consider the deformed configuration of the aircraft with the specific flight conditions, where the flight load and geometric nonlinearity are important.Consequently, the aeroelastic stability will be changed dramatically comparing with the linear analysis.The complete nonlinear aeroelastic problem is quite complex, but for the engineering calculation of subsonic aeroelastic stability, many factors can be trimmed off.The assumptions can be made as follow, 1) the structure vibrates slightly about the state of nonlinear equilibrium state; 2) the aerodynamic nonlinearity due to large attack angle can be ignored for the design limitation of flight load and stall, 3) the non-planar effect of the aerodynamics can be ignored.Under these assumptions, it is adequate that the frequency domain method of the traditional aeroelastic stability analysis can be used with little modification.After calculation of nonlinear static equilibrium, the normal modes can be carried out, and the unsteady aerodynamics in the frequency domain can be calculated by the doublet lattice method.Thus the p-k method [8] can be used to give out the critical unstable speed of the system.

Geometrical nonlinear elasticity
Flexible long slender wing standing large aerodynamic forces has finite bending and torsion deflection, so the infinitesimal deformation condition is disobeyed, while the material is thought to be not beyond the elastic limitation for the little strain.This means that the nonlinear geometric equation should include the quadric term of the displacement differential, and that the nonlinear force equilibrium equation should be defined with reference to the deformed state of the structure.The incremental finite element method is commonly used to solve structural geometric nonlinear problems, which has two formulas called Total Lagrange Formulation (TLF) and Updated Lagrange Formulation (ULF) [9].The ULF is presented and used in the current work.
The relationship between nonlinear Lagrange/Green strain and displacement is where, t u i,j means the partial derivative of displacement component u i to the coordinate x j at time t.The conjugate Kirchhoff stress tensor S ji at time t satisfies, where, t n j is direction cosine of small area element ds at time t, dT j is the corresponding surface force in which the follower force effect is considered.The linear elastic constitutive relation is as follows, D ijkl is the elastic tensor, which has different form for isotropic or anisotropic material.

Finite element method
FEM based on energy principles is an effective approach to solve structural problems.For the geometric nonlinear problems with follower forces, the incremental FEM is used.The strain ε ij can be decomposed into linear part e ij and nonlinear part The stress is decomposed by increments, where t S ij represents the equilibrium stress at time t, and t τ ij represents the incremental stress to be calculated at each time step. t+∆t The integral equation is established by linearization in each incremental step, where, t+∆t Q is the incremental outer force including the aerodynamic force, engine thrust and gravity, etc., at the new time step.Considering a number of shape functions, the relationship between strain and deformation is presented as Substituting them into Eq.( 6), the element governing equation results in The corresponding dynamic equation is,

Dynamic linearization
The assumption of small amplitude around the static equilibrium state is suitable in many dynamic problems, including dynamic stability.That is where, ū is the large deflect equilibrium deformation from Eq. ( 8).x is a small vibration deformation.According to equation ( 9) and static equilibrium condition, the vibration equation of the system under steady forces reduces to where M T is the inertial matrix of the structure at static equilibrium configuration and K T is the corresponding stiffness matrix.From Eq. ( 11), the mode shapes and frequencies are deduced.

Unsteady aerodynamics
The subsonic doublet lattice method in frequency domain is used to calculate unsteady aerodynamics, ignoring the non-planar effect of aero surface, and that the local attack angle is thought to stay below the stall limit.The unsteady aerodynamic pressure at the aero element is where, ρ is the air density at a given altitude; V is the air speed of inflow; φ H is the component of the mode shape normal to the aero surface at the control point of the aero element; φ H is the mode shape slope along the inflow direction at the control point of the aero element; q is the general coordinate of the structural mode; k = bω/V is the reduced frequency, ω is the radius frequency, and b is the reference chord length; C is the aerodynamic coefficient matrix.The details can be found in [8,10].

Flutter analysis
Selecting a number of structural modes, the general aeroelastic equation in the frequency domain is established.
where, M is the general inertial matrix; K is the general stiffness matrix; A is the general aerodynamic coefficient matrix, which is a complex function of the reduced frequency.Equation ( 13) is usually solved by the p-k method in the same way as in traditional flutter analysis.Transforming the equation as follows, then it turns into an eigenvalue problem: Where, A R and A I are the real and imaginary parts of A respectively; p = ω(γ ± i) is the complex eigenvalue of the system at the given value of V .The structural damping coefficient is g = 2γ, and the frequency is f = ω/2π.The critical flutter speed can be identified through the V-g locus (where g goes through zero to positive values), and the flutter frequency from the V-f locus [8].

Analysis procedure
Because the structural deformation and stress affect the system dynamic characteristics, the aeroelastic stability of the high-aspect-ratio wing should be a unified problem of statics and dynamics.So, the finite element model of a real system should give out exact information of the structural deformation and stress distribution, as well as the dynamic characteristics at the same time.
The stability analysis would be carried out as follows: firstly, the nonlinear static problem is solved at a given flight condition.The air load redistribution and the trim state of static aeroelasticity are calculated to get the exact flight condition, and rigorously, the object should be a complete airplane.Secondly, calculate the vibration modes and frequencies at the state obtained from the first step, and the corresponding unsteady aerodynamics as well.Therefore, the stability of the system at this state can be confirmed with traditional flutter analysis methods.Finally, iterating the flight speed and repeat the first two steps, the instability margin can be determined.But in some cases, the load distribution patterns of different flight conditions, such as load factors, or even maneuver types, would change the aeroelastic characters.
In this paper, the aim is not to detail the nonlinear air load redistribution and trim analysis.The dynamic characters and aeroelastic dynamic stability of the high-aspect-ratio wing with pre-stress and finite deformation is the main focus of the work.So, the mode shapes and frequencies of the wing under design load are calculated, and the modes coupling mechanism of the flutter is analyzed.Furthermore, the effect of the system parameters, such as load scale, frequency of the key mode, is studied.

Example study
A composite imbedded high-aspect-ratio wing with large deflection is selected to illustrate the structural geometrical nonlinearity effects on the aeroelastic stability.The full span of the wing is approximately 19.0 meters and the aspect ratio is about 18.0.The average swept angle is about 5 degree backwards.The material fiber directions are designed to make the wing have outer-wash, which means the upward bending results in the nose-down twist [11].
The wing is analyzed at the nominal flight condition, where the aerodynamic force is call designed load in this paper.The dynamic pressure of the designed load is q D , and the equivalent flight speed at sea level is V D .In the static analysis, the steady aerodynamic force is always normal to the aero surface to consider the "follower force" effect.

Static and normal modes analysis
Figure 1 shows the finite element model of the un-deflected wing which is elastically constrained at the wing/fuselage connect points, and the nonlinear static deflection under design load with the follower force effect.Table 1 shows the ratio of vertical displacement of the wing tip to semi-span.With the same load amplitude, the nonlinear result is greater than the linear one.
Table 2 gives out the unloaded and loaded structural mode frequencies divided by a reference frequency f * within the concerning range, as well as the percent of the frequency differences relative to the unloaded frequencies.Most of the structural mode frequencies at the nonlinear equilibrium state are less than the unloaded linear ones, except for the aileron deflection and torsion modes.The primary effect of the geometrical nonlinearity is the stiffness coupling between the chordwise bending and torsion, which causes the frequency of the chordwise bending to decrease and the frequency of the torsion to increase.At the same time, the pre-stressed chordwise bending shapes manifest much more coupling with wing torsion than the unloaded ones, as they have quite a large torsion component.
The greater deflection and the lower frequencies of the loaded nonlinear wing mean that the stiffness of the structure is weakened by the geometrical nonlinearities, such as the displacement and stress status.Because of the "follower force" effect of the steady aerodynamics, the wing is compressed due to the force component along the spanwise.So the nonlinear structure stiffness is overall reduced when compared to the linear one.

Aeroelastic stability
Based on the calculated mode shapes at the nonlinear static equilibrium state, the unsteady aerodynamics in the frequency domain is calculated by the double lattice method, and then, aeroelastic stability analysis is done at the given load condition without mode damping.Table 3 lists the stability and the predicted critical speed divided by the reference speed V * using linear and nonlinear pre-stressed approaches.Because the load and the deformation are all symmetrical, there is no aeroelastic coupling between the symmetric and anti-symmetric modes.The V-g locus

third flap bending
Anti-sym.first chordwise bending of nonlinear pre-stressed wing is showed in Fig. 2, in which the flight speed V D , i.e. the flight speed of the designed load, is presented by the dashed-line.Only the first six anti-symmetric modes are included in the Fig. 2.
In the linear analysis, the system is stable under the flight load, and the predicted flutter is the aileron deflection with the wing bending type.The predicted flutter speed V F /V *, about 98.5, is greater than the flight speed V D /V *, about 66.2.That means that the linear wing is aeroelastically stable below the flight speed V D .Meanwhile, the nonlinear analysis gives an unstable result under the flight load, and the flutter type is changed to the coupling of the first flap and chordwise bending, leading to a very low unstable speed.It indicates that the nonlinear wing is unstable at the flight speed V D .The dramatic change of aeroelastic stability is mainly induced by the nonlinear effects and the switch of the flutter mode.The strengthen fiber orientation of the wing is designed to make it have outer-wash effect under flight load, so that the flap and chordwise bending are coupled with torsion even in linear analysis.When the geometric nonlinearity of large deformation is considered, this coupling is strengthened.Consequently, the aeroelastic stability becomes worse because the chordwise bending modes of the deformed wing have quite large torsion component and lower frequencies.

Load scale investigation
At the given load condition, the flutter mode of the wing is changed.The variation of the new flutter depending on the load scale is studied, in which the load distribution is considered to remain unchanged.Figure 3 shows the trend of the predicted flutter speed caused by the first flap and chordwise bending modes, represented by the symbol-line.The dashed-line represents the equivalent air speed of the load condition at the sea level.The calculated points in the upper area of the dashed-line mean that the system is stable at the given load.The intersection implies that the system is critically stable when the load scale is about 0.29, where the ratio of flight speed divided by reference speed V * is about 36.0.
When the load scale is equal to zero, the nonlinear calculation is equivalent to the linear analysis.In the linear and the small load scales nonlinear analysis, the first two modes are weakly unstable, the unstable speed of which is not the critical speed of the system.As the load scale gets larger, the nonlinear effect is strengthened and the chordwise bending mode has much more torsion component.So the unstable coupling of the two modes gets stronger and the flutter speed drops dramatically.

Chordwise bending frequency investigation
From the load scale analysis, it can be concluded that the coupling of the first flap and chordwise bending modes brings aeroelastic instability of the pre-stressed wing.And the chordwise bending frequency and mode shape are affected significantly by the nonlinearities.Here, the scale of the first chordwise bending frequency is studied to illustrate the variation of the critical speed of instability, while the mode shape remains unchanged.Figure 4 shows variation of the predicted flutter speed V F /V * with different chordwise bending frequency scale, where the origin chordwise bending frequency f /f * 4.03 is noted as unit.The flutter speed gets larger quickly as the frequency scale increases, because the frequency interval between the first bending and the first chordwise bending is getting larger, reducing the aeroelastic coupling between the two modes.That means that strengthening the chordwise bending stiffness can improve the aeroelastic stability of the wing efficiently.

Conclusion
The geometric nonlinear equilibrium state of a composite high-aspect-ratio wing is calculated, and the nonlinear stiffness of the pre-stressed wing is considered to carry out the aeroelastic stability analysis with a small perturbation assumption.Furthermore, the parameter analysis illustrates that the load condition changes the dynamics character of the system dramatically, leading to the decrease of flutter speed for the discussed model.Since nonlinearity induces the coupling of chordwise bending and torsion, which makes the chordwise bending to have a quite large torsion component, the chordwise and flap bending modes manifest themselves at very low flutter speed.Strengthening the chordwise bending stiffness can improve the aeroelastic stability of the wing.
Otherwise, the mode damping is omitted in the stability analysis.Especially for the chordwise bending motion of the composite structure, the mode damping and the aerodynamic drag would play essential roles of promoting the stability.However the structural damping and the aerodynamic drag values are hard to obtain accurately, the analysis without these factors usually gives out a conservative prediction of stability.For the high-aspect-ratio wing, there will be the LCO when the flutter has happened at the flight speed larger than the critical speed.The LCO of the simply beamlike wing had been researched by theoretical and experimental methods [6], but there is little literature to report the LCO of a complex wing or a complete airplane.They should have some common phenomenon and mechanism between simple and complex structures.
The calculated results have shown that for the High-Altitude Long-Endurance aircraft with high-aspect-ratio wing, the linear flutter analysis would give not only an imprecise result, but a possible dramatic mistake.The geometric nonlinearity of the pre-stressed structure is important and should be considered rigorously in the aeroelastic stability analysis.

Table 1
Vertical displacement of the wing tip Fig.1.Finite element mesh of the wing and its nonlinear deflection.

Table 2
Mode frequencies of the wing divided by f *