Chao AN, Chao YANG, Changchuan XIE, Lan YANG

School of Aeronautic Science and Engineering, Beihang University, Beijing 100083, China

KEYWORDS

Abstract In this paper,a framework is established for nonlinear flutter and gust response analyses based on an efficient Reduced Order Model(ROM).The proposed method can be used to solve the aeroelastic response problems of wings containing geometric nonlinearities. A structural modeling approach presented herein describes the stiffness nonlinearities with a modal formulation. Two orthogonal spanwise modes describe the foreshortening effects of the wing. Dynamic linearization of the ROM under nonlinear equilibrium states is applied to a nonlinear flutter analysis, and the fully nonlinear ROM coupled with the non-planar Unsteady Vortex Lattice Method (UVLM) is applied to gust response analysis.Furthermore,extended Precise Integration Method(PIM)ensures accuracy of the dynamic equation solutions. To demonstrate applicability and accuracy of the method presented, a wind tunnel test is conducted and good agreements between theoretical and test results of nonlinear flutter speed and gust response deflection are reached. The method described in this paper is suitable for predicting the nonlinear flutter speed and calculating the gust responses of a large-aspect-ratio wing in time domain.Meanwhile,the results derived highlight the effects of geometric nonlinearities obviously.

1. Introduction

High-Altitude Long-Endurance (HALE) aircraft, which are representative of very flexible airplanes,have recently attracted extensive attention. Because of their lightweight design and substantial flexibility, the wings of HALE aircraft produce large deformation during flight; however, such large deformation of the wing will lead to vital changes in the aerodynamic configuration and stiffness characteristics, and therefore constitute the special aeroelastic problem of large flexible aircraft.Geometric nonlinearity is a very important factor that affects the aeroelastic stability and response. To satisfy the design requirements of HALE aircraft,in 1999,Hodges and Patil proposed a geometric, nonlinear aeroelastic problem for fixedwing aircraft and emphasized the importance of corresponding research.1,2Subsequently, many studies have been conducted across the world on large flexible HALE aircraft to investigate their static deflections,Limit Cycle Oscillation(LCO)behavior and responses. Tang and Dowell studied the aeroelastic response of large-aspect-ratio wing with nonlinear beam theory and ONERA dynamic stall model.3Frulla et al. investigated the dynamic characteristics of flexible wing.4They pointed out that chordwise stiffness and elastic axis eccentricity parameters play an important role in nonlinear aeroelastic behavior. Abbas performed response analysis of HALE aircraft considering the imbalanced effects originated from noncoincidence of center of mass with the elastic axis.5Kim and Strganac studied aeroelastic response of a large flexible wing including the store-induced effect.6Zhang and Xiang have studied the LCO, velocity and drag effect with structural geometric nonlinearities on the static and dynamic aeroelastic characteristics.7Eskandary et al. presented analysis results of large-aspect-ratio wing considering various parameters such as stiffness ratio, mass ratio and distance from elastic axis to mass center.8Two famous frameworks, namely, the UM/NAST toolbox and the DLR solver, are employed for these types of aeroelastic problems associated with very flexible aircraft.9Nevertheless,all studies have shown that their aeroelastic behavior is very different when considering nonlinear flexibility.

Because of variations in the aerodynamic configuration and structural stiffness with deformation, the modeling of large flexible structures is vital for geometric aeroelastic analysis.The nonlinear Finite Element Method (FEM), which is a traditional modeling method, has been widely used in geometric nonlinear analysis through an iterative process in which the aerodynamic mesh is updated in every loop and interpolated from nonlinear deformations.10Xie et al. used nonlinear FEM and strip theory to simulate aeroelastic behavior of a metal wing under large deformation.11As another example,Schetz et al. used quasi-nonlinear aeroelastic analysis to predict the flutter speed based on the static equilibrium state.12Furthermore, the intrinsic beam model by Hodges is an exact model that describes the beam dynamics.13Patil et al. applied the exact intrinsic beam model in the nonlinear aeroelastic analysis with linear state-space aerodynamics.14Palacious and Cesnik applied this model coupled with finite-state inflow theory to analyze nonlinear dynamic behavior.15Wang et al.used the exact intrinsic beam model for discrete gust response simulation of very flexible aircraft with unsteady vortex lattice method.16Multi-body dynamic simulations constitute another tool for analyzing aeroelastic systems with arbitrary types of nonlinearities. For instance, Kruger et al. used this technique to research the stability and control of HALE aircraft in which different sections are connected together.17,18The advantage of multi-body simulation is the convenience to consider the interactions between rigid body motion and aeroelastic deformation.19Other structural modeling methods, such as co-rotational formulation and strain-based FEM, have also been applied in geometric nonlinear aeroelastic analysis.20,21However, beam modeling may be insufficient for capturing the complex structural details of aircraft wings, and nonlinear FEM may be excessively costly for performing dynamic simulations.

Modal analysis is the standard method for linear aeroelastic analysis. However, large deformations change the stiffness characteristics, and such deformation, especially including the foreshortening effect, cannot be represented by the linear modal approach; moreover, geometric nonlinearities cannot be captured in a linear manner. A Reduced Order Model(ROM) based on the modal form was developed to resolve these problems. Muravyov et al. used the surrogate model to satisfy the structural nonlinear dynamics equation.22In addition, Mignolet et al. reported that the nonlinear stiffness can be described using an equation with the quadratic and cubic terms of the basis function.23Based on that formulation,McEwan et al. implemented the Modal/FE (MFE) method to conduct analysis with a number of serial static load cases; the MFE approach is useful for performing dynamic simulations of shells under intense acoustic excitations.24Hollkmap and Gordon improved this method to recover the in-plane deformation with the outer-plane deformation.25Furthermore,Harmin and Cooper implemented the MFE method to model the geometric nonlinearities of large flexible wing in aeroelastic analysis.26However,although the vertical deflection of a wing can be solved exactly,the foreshortening effect cannot be modeled. Consequently, Kim et al. used Proper Orthogonal Decomposition (POD) to capture the foreshortening effect;the spanwise deflection was solved, but the foreshortening effect could not be reflected directly in the nonlinear stiffness matrix.27Cestino et al. reduced the partial differential equations of nonlinear beam model to a dimensionless form in terms of three ordinary differential equations with Galerkin’s method. Nonlinear flutter analysis is presented based on nonlinear equilibrium condition.28Additionally, previous studies rarely considered the follower force effect.

This paper is committed for aeroelastic problems associated with flutter and gust responses including geometric nonlinearities based on modified structural ROM and non-planar aerodynamics. First, the structural ROM builds on Ref.23with modification. Nonlinear stiffness coefficients of structure are described in a specified polynomial formulation. Except linear modes, additional orthogonal spanwise modes are used to describe the foreshortening effects of the wing. Displacements under the follower forces as load cases are solved by nonlinear FEM. Then unknown nonlinear stiffness coefficients are obtained from a curve-fitting procedure, and useless terms are rejected by stepwise regression. The structural equations are obtained for geometrical nonlinearity problem.

Second, the stability problem can be solved through the dynamic linearization of the structural ROM under static equilibrium states.Meanwhile,the gust response analysis based on the structural ROM and the non-planar Unsteady Vortex Lattice Method (UVLM) can be solved in time domain. Application of extended Precise Integration Method (PIM) guarantee the accuracy and stability in numerical calculation. Accordingly, a wind tunnel test is performed to validate the analysis framework. Comparison of nonlinear flutter speed and gust response deflection in theoretic calculation and wind tunnel test measurement is presented, applying on a large-aspectratio wing model. The theoretic and test results show that the method presented in this paper is both effective and accurate.

2. Formulation

2.1. Structural ROM

2.1.1. Structure equation

The development of the structural ROM is based on equations derived from a Galerkin approach to solve geometric nonlinear dynamics in a weak form.23The equation of motion of the structure may often be given in dynamic equations as

where the tensorSis the second Piola-Kirchoff stress tensor,the tensorFis the deformation gradient tensor, ρ0is the mass density of the structure,andb0is the vector of the body force.Xdenotes the position vector of the structure in the reference configuration, andxdenotes the deformed position vector insomuch that the displacement vector is given asu=x-X.The deformed gradient tensor can be defined as

and

To approximate the solution with the Galerkin approach,a set of basis functions that describes the displacement can be written as

where Φij（X）satisfies all of the boundary conditions and represents the componentiof thejth basis function,qjrepresents the time dependent generalized coordinates, anduirepresents the components of the displacement of pointXin the reference configuration. In this research, a truncated basis of linear modes are chosen as the basis functions with no consideration of foreshortening effects first,and the truncated basis of linear modes includes the first two bending modes, the first two torsion modes and the first chordwise bending mode.

Applying Eq. (5) to Eq. (1), which should be expressed in weak form and proceeding with a Galerkin approach as Mignolet shows,23a third-degree polynomial form describes the nonlinear relationship and the dynamic equation in terms of the generalized coordinates can be expressed as

whereMirepresents the modal mass term of theith basis function,andEiis the modal stiffness term of theith basis function.The formulation of the nonlinear dynamic equations corresponding to theith basis function can be written as

That is to say, the structure dynamic equations in modal space have been obtained and generalized coordinatesqiare modal coordinates namely.

The stiffness is related to the deformation of the structure considering geometrical nonlinearities. The tangent stiffness equation for the incremental forces and displacements which are expressed with generalized coordinates in modal space can be written as

whereF= [F1，F(xiàn)2，···，F(xiàn)n]Trepresents the vector of the modal force andq= [q1，q2，···，qn]Tis the vector of the generalized coordinates.Considering Eq. (9),the termsKTijof the tangent stiffness matrixKTunder an equilibrium deformationq*can be described as

The mass matrixMis assembled as Eq. (9) including the modal mass term:

This is generally a good approximation, even for large structural deformations.

As a consequence, by introducing the harmonic oscillating assumption, the eigenvalue problem of a structure under a stable deformation can be expressed as

where ω is the vibration circular frequency under equilibrium deformation. The eigenvectors, which are also called quasimodes,can also be solved;they are useful in the stability analysis of geometric nonlinear aeroelastic problems.

2.1.2. Foreshortening effects under large deformations

Foreshortening effects under large deformations have a substantial influence on geometric nonlinear problems, especially aeroelastic problems. The shortening effect on the spanwise projection of a wing will change the aerodynamic distribution considerably. A truncated basis of the linear modes is chosen as the basis function in the ROM model presented herein.However, a few linear modes truncated will not show the shortening effect of the wing in ROM analysis. Therefore,two orthogonal spanwise modes are taken into the established ROM to describe the foreshortening effects of the wing. A combination of truncated linear modes and orthogonal spanwise modes is generated as a basis function in the nonlinear structural ROM described above.

2.1.3. Regression analysis for solving nonlinear stiffness coefficient

Evidently,if there are a set of specified static loads and corresponding structural deformations, the unknown nonlinear stiffness terms related to the applied loads and the structural displacement resultant can be solved by using regression analysis. The set of specified loads and corresponding structural deformations can be denoted as the static test load case and calculated by a commercial FEM software package.

Considering that there are NT sets of static test load cases,NT sets of corresponding displacements can be obtained after performing static FEM analysis of NT sets of test loads on the model in the commercial finite element software. The loads and deformations are translated into the modal space.To find the unknown nonlinear modal stiffness terms, the left side of Eq. (14) can be fitted with a curve in a regression progress.The regression problem can be presented as

Einstein summation convention is applied in Eqs. (14) and(15). It should be noted that the superscripts without bracket denote the serial number of test cases from 1 to NT instead of power value.

As described before, two bending modes, two twist modes,one edgewise bending mode and two spanwise modes are used to establish the nonlinear ROM. The number of nonlinear stiffness coefficients will be large;as a consequence,overfitting may occur. Therefore, a stepwise regression method will be used to decrease the number of nonlinear stiffness coefficients.

2.1.4. Strategy for generating test load cases

Through the abovementioned analysis, regression analysis is performed using the commercial software package on the actual deformation and load testing after FEM analysis, and thus the accuracy of the nonlinear stiffness coefficients directly depends on the rationality of the selected test load case,which is related to the success of recovery of the nonlinear structure equation. The selection of test load cases must emphasize that the aerodynamic force on the wing should be a follower force,which more closely resembles the actual characteristics of the aerodynamic force.That is,taking an oriented load as the load test case cannot satisfy the requirements. In this paper, the aerodynamic force under the deformation combined bending and torsion modes is chosen as the test load case. The formulation of the wing deformation, which generates aerodynamic forces, should be

Here, ｛φi｝bendand ｛φi｝torsionrepresent the investigated bending and torsion modes, respectively, andai，jrepresents the scalar mode weight factors,through which the selected test cases contain the nonlinear characteristics of the structure investigated in our research.

2.1.5. Selection of nonlinear stiffness coefficients

The selection of the nonlinear stiffness coefficients constitutes a main problem in establishing a structural ROM. In general,to ensure high accuracy for the ROM, the nonlinear stiffness coefficients should have sufficient flexibility to approximate the training points. However, if the number of coefficients is excessively large,overfitting may occur.The number of nonlinear stiffness coefficients utilizingNmodesNCis

where

With 7 modes used in our research, the number of nonlinear stiffness coefficients is 112.Therefore,the overfitting problem cannot be neglected in the nonlinear structural ROM method.

Stepwise regression is a statistical procedure used to select variables in a regression space.This approach is computationally efficient because it provides intermediate statistical information at each stage of the calculation; this information is then used to select the most appropriate coefficients to be added into the model. The procedure of stepwise regression can be summarized as follows:

(1) The regression analysis problem corresponding to Eq.

(2) Choose the term in unintroduced variables which has the maximum value of partial regression quadratic sum as the introduced term for regression equation and do the significant test. If it passes the test, the term can be reserved in the equation.Otherwise,it will not be considered.

(3) Choose the term in regression equation which has the minimal value of partial regression quadratic sum and do the significant test. If it passes the test, the term can be reserved in the equation. Otherwise, it will be removed.

(4) Repeat the above Process (2)–(3), until no term can be introduced or removed.

By stepwise regression procedure, the best nonlinear ROM model will be determined. The joint hypotheses test (F-test) is used for the specified criterion in significance test and is expressed as29

2.2. Non-planar unsteady aerodynamics method

The UVLM, which constitutes a time domain aerodynamics computation, is an efficient method employed to calculate the aerodynamic loads for various aircraft. This technique can be coupled with structure dynamic computations, such as the structural ROM, to obtain the response results for aeroelastic systems. Additionally, the exact boundary condition at the actual wing surface is satisfied,and thus the UVLM can be conveniently applied to very flexible wings whose aerodynamic surfaces are subjected to large deformations. The UVLM is based on full potential equations without linearization,and it can effectively reflect the unsteadiness effects of the 3-D low-speed flow around a flexible lifting surface.30

Vortex ring elements are used to discretize the boundaries of the aerodynamic surface in the UVLM for both the wing and the wake in Fig.1.U（t），V（t），W（t） are the components velocity of air flow in aerodynamic coordinates system withx，y，zaxes, Δbij，Δcij，nij，Γijrepresent the span length, chord length, normal vector and vortex strength ofijth lattice,t0，Δtare the initial time and time step. The leaving segment of the vortex ring is placed on the panel’s quarter-chord line,and the collocation point is located at the center of the three-quarter-chord line. The whole flow domain is represented by vortex rings, and the aerodynamic influence coefficient can be obtained via Biot-Savart’s law.

Specific vortex distribution on the wing is determined by the geometric exact boundary condition expressed as

where Ψ is the full potential function about vortex,vis the velocity of the local boundary movement, andnis the local normal vector of the bound vortex. Formulation above can be developed as

whereAis the normal aerodynamic influence coefficient of the bound vortex, Γ is a vector consisting of the panel vortex circulation on the wing,V（t）is the vector of time dependent kinematic velocity due to the unsteady coming flow with gust excitation discussed in this paper,andVwis the vector of velocity induced by the wake vortices.

At each discrete time step in the computation, the wing is moved along its flight path,and each trailing edge vortex panel sheds a wake panel with a vortex strength equal to its circulation in the previous time step.For a gust response problem,the gust load represents the main factor influencing the structure response aerodynamically, and thus the fixed-wake model(shown in Fig.2), which is much more computationally efficient than the free-wake model,is deemed sufficient for a single flexible wing response problem.

2.3. Surface spline interpolation

Fig.1 Non-planar UVLM model.

Fig.2 Fixed-wake model.

Surface spline interpolation is applied to couple the aerodynamics and structure. In most situations, the configuration of the structure is embedded into a 3-D space; however, while the configuration of the deformation is usually 3-D, the undeformed configuration might be 1-D, 2-D or 3-D.

When structural grids and the corresponding deformation vectorUSare confirmed,the deformation vectorUAof aerodynamic grids can be obtained by displacement interpolation:

Here,Gdenotes the spline matrix. The structural equivalence between the structural force system and aerodynamic force system is satisfied when the virtual work corresponding to the virtual deflections generated by the aerodynamic loadsFAand their equivalent structural forcesFSis equal in the force interpolation:

2.4. Response analysis methodology

An analysis of the static aeroelastic and dynamic stability characteristics is necessary to perform before the response analysis.An aircraft structure is subjected to steady aerodynamic loads before encountering a gust regardless of whether it is under real flight conditions or wind tunnel test conditions.Therefore,the dynamic characteristics will change with variations in large static deformation. All of these factors will have an enormous influence on the response analysis.

Static aeroelastic analysis can be implemented with an iterative computation including the non-planar Vortex Lattice Method (VLM) under steady conditions and the structural ROM/FEM.31Surface spline interpolation is also used for exchanging information during the coupling process. With the results from the nonlinear static aeroelastic analysis,stability analysis can be conducted through quasi-mode analysis, which considers the nonlinear stiffness and stress effects under large deformations. Coupled with the non-planar Doublet Lattice Method(DLM),the flutter results under a nonlinear large deformation equilibrium state can be obtained. This iterative procedure is illustrated in Fig.3. A convergent result can be obtained through nonlinear static aeroelastic analysis with an inner iterative process under given flight velocity conditions. Then, the flutter equation can be solved by coupling the non-planar DLM computation with the quasi-modes at the nonlinear static aeroelastic equilibrium state. The flutter speed can be precisely determined once the flutter analysis and the flight velocity in the current iterative process reach the same flight speed; otherwise, the analysis varies the flight speed, after which the calculation is repeated.30A dichotomy is commonly used to search for the exact flutter speed. If the predicted flutter speed is higher than the current computed flight speed, an intermediate speed can be selected as the next computed flight speed.

Fig.3 Nonlinear flutter analysis flowchart.

The gust response analysis is implemented in the discrete time domain. At the beginning of each computational time step, the unsteady aerodynamic load is computed by the non-planar UVLM. The structural deformation and velocity at the end of the preceding computational time step will be treated as the initial condition in the next structural transient dynamic analysis, thereby guaranteeing the continuity of the structural response analysis.Each structural transient dynamic analysis is carried forward over a time step, during which the unsteady load is kept unchanged. The structural transient dynamic analysis is calculated by the nonlinear integration method with the nonlinear ROM described above, the results of which will be accurate if the time step is sufficiently small.The resultant structural displacement and velocity are used to update the aerodynamic surfaces and the exact geometric boundary conditions for the next step of the aerodynamic computation. Furthermore, although the computation is implemented in the discrete time domain, the structural displacement and velocity are continuous, and the updated aerodynamic computation helps ensure that the unsteady aerodynamic computation is both accurate and practical.The analysis procedure is illustrated in Fig.4.

The nonlinear dynamic solution is solved by the following equation:

whereFi（t）=Φif（t） is the generalized unsteady aerodynamic load in modal space under an atmospheric gust andf（t） is the unsteady aerodynamic load solved by UVLM in timedomain. ‘1-cos’ gust model is chosen as gust model in this paper.Eq.(25)constitutes a set of Ordinary Differential Equations (ODEs). Thus, numerical integration represents an important component of the ROM solutions for a set of ODEs. The Finite Difference Method (FDM) is usually applied to solutions that could introduce errors and numerical difficulties. In particular, in this research, the number of nonlinear terms in the equation will grow rapidly with the addition of modes chosen for the structure ROM. Consequently, an efficient and accurate numerical integration method is necessary.

A PIM32is proposed for the numerical integration of the ROM.The PIM will be extended to the time-variant nonlinear ODEs system and applied to the ROM solutions.Accordingly,the ROM can be expressed in matrix/vector form as

wherev（t） is ann-dimensional vector function to be determined,Ais a givenn×nconstant matrix, andf（t） is a givenn-dimensional vector that includes the nonlinear stiffness terms and external forces.The response can be computed as follows:

Fig.4 Time-domain nonlinear gust response analysis flowchart.

whereT=exp（Aη） is the exponential matrix, η is the time step. If an approximation method with linear interpolation is used in the intervaltk～tk+1:

wherer0is constant vector andr1is linear approximation vector in linear interpolation of nonlinear vectorf（tk+ξ）.

Eq. (28) can be integrated as follows:

The exponential function exp（Aη） is solved by the PIM,providing a highly precise numerical result that approaches the full computational precision.

3. Wing model and wind tunnel test system

3.1. Wing model

To validate the accuracy of the ROM and the complete aeroelastic response analysis including geometric nonlinearities, a wing model is constructed.The parameters of the wing,which is characterized by a large aspect ratio and a considerable flexibility, are shown in Table 1.

A rectangular steel ruler was chosen for the beam simulation of the wing. The beam is located at the 50% chord line of the wing, and the dimensions of the rectangular section are 35 mm×1.5 mm. The density of the material is 7.75×103kg/m3, and the elastic modulus is 219 GPa. The wing frame includes 8 light wood boxes and each box attaches to the wing beam at a single point. A sufficient amount of clearance is left between each section to ensure that no stiffness will be added to the wing beam by the external shell. Wingtip stores are used to regulate the flutter characteristics.The wingtip store is 200 mm long, and the weight is 62 g. The CATIA model, FEM model and wind tunnel test model are shown in Fig.5.

The main linear modes of the wing model are presented in Table 2.The frequencies of the first two modes are low,which means that the flexibility of the model is very high.

3.2. Design of wind tunnel test subsystem

A wind tunnel test was performed for the static aeroelastic test,nonlinear flutter test and gust response test of a very flexible wing. All tests were performed in the FD-09 low-speed wind tunnel at the China Academy of Aerospace Aerodynamics.The Reynolds number as maximum test velocity is about 1.4 million and the turbulence coefficient of wind tunnel is lower than 0.1%.

Table 1 Design parameters of wing model.

Fig.5 Very flexible models.

Table 2 Modes of wing model.

3.2.1. Support system

The wing model is vertically fixed to a support system in a wind tunnel with a lamped root. The root of the wing model can be changed from an angle of attack α=0° to 5° which are related to the undeformed wing model. At the joint between the support system and the wing model, a force balance measures the forces and moments in six directions. Furthermore,the support system is designed to be absolutely rigid.

3.2.2. Gust generator

A gust generator is applied to create the expected gust condition during the wind tunnel test. Two rectangular blades deflect sinusoidally and synchronously. The span length of each blade is 2000 mm, and the chord length of each blade is 300 mm. The distance between the two blades is 600 mm, and the distance from the model to the generator is 2000 mm.NACA0015 airfoils are used for both blades. The deflecting angle of each blade ranges from-7°to 7°.The guide generator in the wind tunnel is shown in Fig.6.

When the blades are deflected sinusoidally at a certain frequency,an approximately sinusoidal lateral gust can be generated in the test field, and the gust velocity can be written as

whereamis the amplitude of the blade deflecting angle andAgis the gust disturbance coefficient, which is relevant to the velocityVand the blade deflecting frequencyf.

3.2.3. Fiber optic sensor and accelerometer

In the wind tunnel test,Fiber Bragg Grating(FBG)sensors are embedded in the structure of the wing to measure the strain,see Fig.7.The FBG sensors are created by exposing an optical fiber to an ultraviolet interference pattern, which produces a periodic change in the core index of refraction.

It is often difficult to measure large deformations during nonlinear aeroelastic wind tunnel tests, as there are various limitations.The FBG sensors are arranged and installed at different locations so that strain values can be obtained at these stations. After strain values are obtained, the rotation angles and the displacements can be calculated by numerical integration with a Strain-Displacement Transformation (SDT).33

An ADXL345 light accelerometer is chosen to measure the vertical acceleration, and it is located at a distance of approximately 33% of the span on the wing beam from the wing tip.

Fig.6 Gust generator in wind tunnel.

Fig.7 GBT sensor and accelerometer on wing beam.

4. Theoretical analysis and wind tunnel test results

4.1. Validation of nonlinear ROM

Fig.8 Validation procedure.

Fig.9 Displacement of wing tip under validation loads.

A nonlinear ROM has been obtained by the abovementioned regression analysis of test load cases and the corresponding deformations. Subsequently, the nonlinear ROM calculations from the structural equations must be validated. Then, the ROM can be reasonably applied to structural and aeroelastic response analyses. The validation procedure is illustrated in Fig.8. The nonlinear FEM analysis is performed in MSC.NASTRAN, the module SOL 106 is used for static nonlinear FEM analysis, and the module SOL 400 is used for dynamic nonlinear FEM analysis.

4.1.1. Static validation with nonlinear FEM

For a static validation of the nonlinear ROM with data from the nonlinear FEM,50 sets of growing validation loads, all of which have distributions of actual aerodynamic loads, are employed as examples. The displacements of the wing tip under these sets of validation loads are presented in Fig.9(a). The selected validation loads result in vertical displacements of the wing tip within the range from 10% to 25%,which satisfies the requirements of nonlinear analysis. The spanwise displacement of the wing tip is also considered.Consequently, it is important in geometric nonlinear aeroelastic analysis which will change aerodynamic distribution.

The calculation results of the 50 cases are shown in Fig.9(b) with a comparison between the ROM solutions and FEM solutions. Here, the black lines with circles represent the relative deviation in the vertical displacement of the wing tip,and the red lines with rectangles represent the relative deviation in the spanwise displacement of the wing tip.The relative deviation between the ROM solutions and nonlinear FEM solutions under these validation loads does not exceed 1%.Therefore, the ROM solutions demonstrate great agreement with the nonlinear FEM solutions.

Fig.10 Response of wing model under load case 1.

4.1.2. Dynamic response validation with nonlinear FEM

Two load cases are established for a comparison for the dynamic response validation of the nonlinear ROM with data from the nonlinear FEM:

(1) Load case 1: 0.5 N follower force applied at the tip of the beam in the vertical direction as a step input at the initial time.

(2) Load case 2:a set of follower forces has a distribution of the actual aerodynamic load as a step input at the initial time, and the load can impose a vertical deflection of nearly 20% the length of the span.

Comparisons of the FEM solution results with the ROM solution results are shown in Figs. 10 and 11. The nonlinear ROM solution results for both the vertical deflection and the spanwise deflection under load cases 1 and 2 coincide with the nonlinear FEM solution results. Therefore, the structural ROM is reliable.

4.2. Nonlinear flutter analysis and test results

Fig.11 Response of wing model under load case 2.

The static aeroelastic stable deformation and flutter speed are calculated prior to performing the response analysis and wind tunnel test. The nonlinear analysis uses the non-planar VLM and non-planar DLM coupled with the structural ROM, and the geometric nonlinearities of the structure and aerodynamics are fully considered. Meanwhile, the linear analysis employs the planar VLM and planar DLM that can be easily solved in MSC.NASTRAN.

Fig.12 Static aeroelastic deformation results.

Fig.12 shows the static aeroelastic deformation solutions under the 3° and 4° conditions with varying velocities from 10 m/s to 18 m/s. The tip displacement of the wing changes with an increase in the velocity. The solutions derived from the nonlinear method proposed in this paper are in good agreement with the wind tunnel test results. In contrast, the linear analysis solutions exhibit a large difference from the nonlinear solutions and test results. The deformation of linear analysis solution is bigger and the differences are observed increasing with velocity growing up. A nonlinear analysis method is therefore deemed necessary for geometric nonlinear aeroelastic problems.

As described above in the context of geometric nonlinear analysis,the stiffness of the structure is a function of the deformation. Dynamic linearization can provide useful dynamic information concerning the static equilibrium. According to Eqs. (11)–(13), a solution of quasi-modes can be obtained.The numerical example of the α=3°and 18 m/s velocity conditions is taken as an example.Figs.13 and 14 show the deformation of the wing model and the distribution of the lift,respectively, while Fig.15 and Table 3 show the quasi-mode solutions. The gray zone denotes the initial configuration of the wing model, the purple zone represents the deformation under the aerodynamic loads at static equilibrium,and the yellow zone represents the quasi-mode solutions.

In the nonlinear analysis,the frequencies of the 1st and 2nd bending modes are less than the unloaded frequencies. The effect of geometric nonlinearity is mainly the existence of coupling between the torsion and chordwise bending;this coupling causes a substantial decrease in the chordwise bending frequencies and an increase in the torsion frequency.Meanwhile,the 3rd quasi-mode as the chordwise mode has a large torsion component;this will have a large influence on the dynamic stability analysis. The larger deformation and the lower frequencies of the loaded wing structure indicate that the stiffness changes considerably due to the presence of geometric nonlinearities. Furthermore, because of the effects of the follower force on the aerodynamics,loads on the wing have a compressive effect compared with aerodynamic analyses that ignore follower force effects.The nonlinear structural stiffness is overall reduced in comparison with the linear stiffness and the nonlinear stiffness without follower force effect.

Fig.13 Deflection under α=3° and 18 m/s velocity conditions.

Fig.14 Aerodynamic distribution under α=3° and 18 m/s velocity conditions.

Fig.15 Quasi-mode solutions.

Table 3 Mode frequencies and descriptions.

The linear flutter analysis results are shown in Fig.16.TheV-gcurves andV-fcurves with the same legend show that the results exhibit a typical bend/twist coupling flutter,and the critical flutter speed is 29 m/s which are irrelative with velocity and α.TheV-fcurves show that the chordwise bending mode does not contribute to the flutter. In contrast, the nonlinear flutter analysis under α=3° performed with the abovementioned iterative procedure reveals different results from the linear analysis.Fig.17 shows the nonlinear flutter analysis results under a static equilibrium of 14 m/s;the flutter speed is 24 m/s,which is lower than the linear analysis results.The iterative procedure is used to narrow this disagreement and search for the exact nonlinear flutter speed. Fig.18 shows the analysis results under a static equilibrium of 19 m/s; the predicted flutter speed is 19 m/s, which represents the critical nonlinear flutter speed under α=3°. Moreover, the flutter speed decreases considerably and the chordwise bending mode contributes to the flutter under large deformation. TheV-fcurves andV-gcurves also have the same legend in Figs. 17 and 18.

Fig.16 Linear flutter analysis results.

Fig.17 Predicted flutter analysis results at 14 m/s velocity.

Fig.18 Nonlinear flutter analysis results at 19 m/s velocity.

Fig.19 Wind tunnel test results of vertical acceleration.

Fig.20 Wind tunnel test results of chordwise acceleration.

The wind tunnel test results are shown in Figs. 19 and 20.Nonlinear flutter occurs when the speed is close to 18.5 m/s.This result agrees with our nonlinear analysis method. Meanwhile, the chordwise bending mode obviously contributes to the nonlinear flutter.

4.3. Gust response analysis and test results

In this section,the time-domain response analysis results of the vertical deflection are compared with the wind tunnel test results. The test deflections can be obtained from numerical integration with an SDT based on measured strain data. In the gust response test,α of wing model is 3°,and the gust generator oscillates to produce a continuous sinusoidal gust.Three cases are taken for comparison between the simulation and test.The values of the velocity and gust frequency in each case are shown in Table 4.

The results of the wind tunnel test and ROM simulation are presented in Fig.21. The choice of the initial value willintroduce some error into the numerical integration with an SDT.Despite the small deviation in the numerical integration,the stable response amplitudes from the gust response analysis solutions based on the nonlinear ROM and the wind tunnel test results are quite similar, proving that the established procedure is accurate and reasonable. Furthermore, because the gust frequency in the test is lower than the frequency of the 1st chordwise bending mode, the response amplitude is nearly unchanged with variation in the gust frequency at a fixed velocity.

Table 4 Comparison cases.

Fig.21 Response results in time domain.

5. Conclusions

(1) In this paper,an effective method for analyzing the flutter and gust responses of a large flexible wing including geometric nonlinearities is developed. The research herein mainly establishes an ROM for geometric nonlinear structural modeling.Moreover,the proposed framework for nonlinear flutter and gust response analyses based on the structural ROM and non-planar aerodynamics is proven effective compared with wind tunnel tests.

(2) The structural modeling method presented herein describes the geometric nonlinear stiffness coefficients in an equation with the quadratic and cubic terms of a linear modal function. The follower force effect is also considered in the ROM. Additional orthogonal spanwise modes are introduced to describe the foreshortening effects. Furthermore, the spanwise deflection can be solved, and the influences of foreshortening effects can be directly included into non-planar aerodynamics calculation.The ROM solutions of both static deflection and dynamic response are in good agreement with the nonlinear FEM solutions.

(3) The dynamic linearization of the ROM under static equilibrium conditions can be used to achieve a solution for the stability. Quasi-modes obtained from this linearization are coupled with the non-planar DLM to predict the nonlinear flutter speed using an iterative process.The flutter speed decreases under large deformations relative to the linear flutter speed.

(4) The fully nonlinear ROM is coupled with the nonplanar UVLM to calculate the gust response in the time domain, and the aerodynamic loads are taken as follower forces corresponding to the actual conditions.An extended PIM is used to solve the dynamic equations.

(5) A wind tunnel test is performed with a large flexible wing to validate the proposed method.The support system, gust generator and measurement system are described accordingly. The wing deflection can be measured by FBG sensors through numerical integration.The simulation results are consistent with the wind tunnel test results, and the proposed method for nonlinear aeroelastic analysis is proven to be effective and accurate.

Acknowledgement

The study was supported by the National Key Research and Development Program of China (No. 2016YFB 0200703).