Flurrer Analysis of Aircrafr Wing Using Equivalenr-Plare Models wirh Orrhogonal Polynomials

Tang Jian（唐?。?，Xie Changchuan（謝長川），Yang Chao（楊超）

School of Aeronautic Science and Engineering，Beihang University，Beijing 100191，P.R.China

（Received 5 March 2013；revised 30 May 2014；accepted 21 June 2014）

Flurrer Analysis of Aircrafr Wing Using Equivalenr-Plare Models wirh Orrhogonal Polynomials

Tang Jian（唐健），Xie Changchuan（謝長川）*，Yang Chao（楊超）

School of Aeronautic Science and Engineering，Beihang University，Beijing 100191，P.R.China

（Received 5 March 2013；revised 30 May 2014；accepted 21 June 2014）

Absrracr：A method of equivalent simplification，using equivalent-plate models（EPMs），is developed.It is to achieve goals of rapid modeling and effective analysis in structural dynamics and flutter analysis of complex wing structures.It is on the assumption that the wing structures discussed are composed of skin，beams and ribs，and the different plate units（such as skin，beam web，rib web）are not distinguished in modeling，which is to avoid the complex pre-processing and make it more generalized.Taking the effect of transverse shear deformation into consideration，the equivalence is based on the first-order shear deformation theory，and it can import the model files of MSC/NASTRAN and process the information to accomplish the equivalent modeling.The Ritz method is applied with the Legendre polynomials，which is used to define the geometry，structure and displacements of the wing. Particularly，the selection of Legendre polynomials as trial functions brings good accuracy to the modeling and can avoid the ill-conditions.This is in contrast to the EPM method based on the classical plate theory.Through vibration and flutter analysis，the results obtained by using EPM agree well with those obtained by the finite element method，which indicates the accuracy and effectiveness in vibration and flutter analysis of the EPM method.

equivalent-plate model（EPM）；flutter analysis；vibration analysis；Ritz method；orthogonal polynomials

0 Inrroducrion

In structural analysis process，the finite element method is well known and has wide applications.There are many finite element method codes，such as MSC/NASTRAN，ANSYS，ABAQUS，and they are widely used in solving engineering problems.While in the preliminary design of aircraft，the structural parameters change quite often，which brings many obstacles to modify the finite element models（FEMs）. The preparation time for FEM data may be prohibitive and for complex structures a detailed finite element analysis（FEA）needs huge amount of CPU time and computation capacity，which diminishes the utility of the finite element method and makes the cost soar［1］.

Flutter problem has been taken into consideration in preliminary design of modern aircraft［2］.To avoid the inconvenience of complex experiments，detailed FEA and complicated calculations in preliminary design，it is necessary to develop simplified equivalent models to simulate complex structures and obtain global solutions. Degrees of freedom are reduced and only necessary modes for flutter analysis are retained in these equivalent models.The key point to the equivalent models is rapid modeling，accuracy and effective analysis.

The equivalent-plate models（EPMs）were first proposed by Giles［1］from NASA Langley Research Center.This method rapidly became a useful tool for structure analysis in designing. Based on the classical plate theory（CPT），Giles［3-4］developed EPM in forms of trapezoidal plates by using the Ritz method.Tizzi［5-6］presen-ted a method which was similar to that of Giles in many aspects.In his work，several trapezoidal segments in different plane（Such as winglets，T-tails，etc）were considered，but the internal parts of wing structures（Spars，ribs，etc）were neglected for simplification.Livne［7］developed his EPM method by using the first-order shear deformation theory（FSDT），taking simple polynomial as the trial function to define the displacements of structures.This method was more accurate than that of Giles.Later，Livne and Navarro et al.［8-9］further developed the EPM method and applied it to geometry nonlinear research of wing box structures.Due to the simple-polynomial trial functions they selected，it was easy to face numerical ill-conditioning problems and the order of the polynomial is limited.Kapania et al.［10］used the Rayleigh-Ritz method with the Chebyshev polynomials as the trial functions to obtain EPM. Chebyshev polynomials were orthogonal so that ill-conditioning problems could be avoided.But in calculation，numerical integration must be used so that CPU time was much longer.Later，Kapania and Liu［11］developed a similar method by using Legendre polynomials as the trial functions. Ranjan Vepa［12］formulated the Reissner-Mindlin first-order shear deformation model based on governing equation for a nonisotropic plate and made the aeroelastic analysis of a wing modeled as EPM.Recently，Krishnamurthy et al.［13-14］developed EPM that could be used to design a wind tunnel model to match the stiffness characteristics of a full-scale aircraft wing box while satisfying strength-based requirements.EPM analysis technique is also presented to predict the static and dynamic response of a wing with damage.Na and Shin［15］developed a method similar to that of Livne to analyze a low-aspect-ratio composite wing with a control surface.

The past EPM method has been widely used in preliminary design phase of aircraft structures in recent years.However，in these studies，the different kinds of plate units（Such as skin，beam web，rib web）of the wing composed of skin，beams and ribs，are distinguished in modeling. That is，the integrals of the methods are computed in different parts，respectively.It brings complex pre-processing work.In contrast，the present method will not distinguish different plate units in modeling by using a different integration method，which is to avoid the complex pre-processing work and make it more generalized.Also，a comparison is made between different plate deformation theories to select a more appropriate one for the present study.

In the present study，the effect of transverse shear deformation is also taken into consideration.This method is based on FSDT，and it can import the model files of MSC/NASTRAN and process the information to accomplish the equivalent modeling.Legendre polynomials are used as the trial functions to define the displacements of the wing structures，and the boundary conditions are approximately dealt with artificial springs，whose stiffness is rather large.Then，the Ritz method is used to solve the problem in order to get equivalent stiffness matrix and equivalent mass matrix.The natural frequencies and vibration modes of the wing structure are easily obtained by solving eigenvalue problem.Particularly，the selection of Legendre polynomials as trial functions can bring good accuracy to the modeling and avoid the ill-conditions.Results of vibration and flutter analysis can be easily obtained by using the modes of EPM.

The main work of this paper aims at derivation of the equivalent stiffness and equivalent mass formula of equivalent plate model，which is the key to the equivalent plate modeling method and forming a complete set of modeling program. Programming code is developed for vibration and flutter analysis with the EPM method，through which we can easily get various vibrational mode shapes and input files for flutter analysis using software ZAERO.As two examples to verify the present method，a rectangular plate and a wing with medium aspect ratio are analyzed，respectively，and the results are compared with those obtained from FEA by using MSC/NASTRAN. Also，to verify the accuracy of the present meth-od，comparison is made between the EPM method based on CPT and that based on FSDT.In addition，we also establish the interface between equivalent plate model and MSC/NASTRAN so that the model files of this commercial software can be imported，and element information can be used for equivalent modeling and further analysis. The present method emphasizes accuracy，fast modeling and effective analysis，in order to get a fast and reasonable estimate of structural and aerodynamic design.

1 Marhemarical Model and Formularions

1.1 Assumprions and Rirz base

To get the equivalent stiffness matrix and mass matrix of structures，the displacements and strains should first be defined.Using the Reissner-Mindlin method，taking transverse shear effects into consideration，F(xiàn)SDT［16-17］is used to define the displacements field of the plate as follows

where u0，v0，w0are the displacements in x，y，z directions，respectively andφx，φythe rotations about y and x axes，respectively.

Then we can obtain

where，if CPT，instead of FSDT is used for the solution of plates，γxz，γyzwill not be considered. Then the strain matrix and the constitutive matrix based on CPT are less dimensional and simpler than those based on FSDT.Therefore，the solution of the plate using the former theory is less accurate.Especially for wing with large thickness-chord ratio，CPT is not reasonable.

According to the Ritz method，u0，v0，w0，φx，φycan be expressed as following forms

Define q=｛qˉU，qˉV，qˉW，qˉX，qˉY｝Tas the displacement description using generalized coordinates，one can write

Then，u0，v0，w0，φx，φyin Eq.（3）can be written as

In matrix form

The Legendre polynomials as trial functions to describe the Ritz base，we have

Due to the orthogonality of the Legendre polynomials，matrix B in Eq.（6）and matrixin Eq.（7）are symmetric and have amounts of zero elements，which makes the solution of the matrix easier and more accurate.If simple-polynomial trail functions are used here，the ill-conditions are prone to occur.

1.2 Isoparamerric rransformarion and numerical inregrarion

The wing we want to analyze is assumed to behave as a plate.For the convenience of calculation，a transformation of the coordinate system established on the finite element is performed from（ x，y，z） to（ξ，η，ζ）.

As shown in Fig.1，the transformation can be formulated as［18］

where（ξi，ηi）is the location and V3ithe vector in the normal direction of the middle plane.

Fig.1 Coordinate system and its transformation

Let integrands be f（x，y，z），the transformation can be written using the chain rule of derivative as

where J is the Jacobian matrix.

1.3 Equivalenr sriffness marrix

The expression of potential energy of the wing structure is

Using Eq.（5），the expression of strain of the wing structure in the global coordinate system is

Substitute Eqs.（14—17）into Eq.（13），and one has

The equivalent stiffness matrix of the structure in terms of q is

The equivalent stiffness matrix of the whole wing includes the matrix of all structures of the wing，such as skin，beam，rib and that of the artificial spring of boundary conditions［7］.

1.4 Equivalenr mass marrix

The expression of kinetic energy of the wing structure is

where V is the velocity vector，and it can be written as

Substitute Eqs.（21，22）into Eq.（20），and one has

Then the equivalent mass matrix of the structure in terms of q is

1.5 Numerical inregrarion of sriffness and mass marrices

In the present method，in order to avoid the complex pre-processing and make EPM more generalized，different plate units（Such as skin，beam web，rib web）are not distinguished in modeling.For a specific structure，numerical integration is used for calculation of Eqs.（19，24）to achieve this goal.

By Gaussian integration，one has ξi，ηj，ζkare the sampling points of Gaussian integration and gi，gj，gkthe weights.

Then，using Gaussian integration，the equivalent stiffness and mass matrices of the structure in Eqs.（19，24）can be written as follows

2 Vibrarion and Flurrer Analysis

2.1 Vibrarion analysis

For free vibration，according to Lagrange equations，natural frequencies and mode shapes of EPM can be easily obtained by solving the following eigenvalue problem

whereλdenotes the matrix of eigenvalues，λ= diag（λ1，λ2，…，λi），λi=is an eigenvalue of the system equations，withe corresponding frequency，and x the corresponding eigenvectors.K and M are the whole equivalent stiffness matrix and equivalent mass matrix of the wing.

Normalized mass matrix and stiffness matrix can be obtained as

2.2 Flurrer analysis

Flutter is a dangerous phenomenon encountered in flexible structures subjected to aerodynamic forces.Due to the interactions between aerodynamics，stiffness and inertial forces on the structure，flutter occurs as the wind speed increases to a point at which the structural damping is insufficient to damp out the motions［2］.

Flutter equation of motion is

Unsteady aerodynamic force is calculated by software ZAERO using the ZONA6/ZONA7 method，based on which the results of flutter analysis can be easily obtained by the use of g method［2］.

3 Resulrs and Discussion

3.1 EPM of recrangular plare

A rectangular plate is selected to verify the present method.Its geometric and material parameters are given with the span of 2.0 m，the width of 0.7 m，the thickness of 0.02 m，the mass density of 2.7×103kg/m3，and E= 200 GPa，μ=0.3.The plate is clamped at the root and the finite model is shown in Fig.2.

Fig.2 Finite model of rectangular plate

Modal analysis based on FSDT and CPT using EPM are conducted.The comparison between the natural frequencies obtained by the present methods and those by the FEA software MSC/NASTRAN is shown in Table 1.The comparisons of the mode shapes obtained by FEM and EPM are shown in Fig.3.

Table 1 Narural frequencies of recrangular plare

Fig.3 Mode shapes of rectangular plate

Separate flutter computation of the rectangular plate using finite element model and equivalent plate model.Comparisons are made between the flutter speed as well as the flutter frequencies as obtained by the two methods in Table 2.Comparisons between V-g curves and V-f curves are also shown in Fig.4，respectively.

Table 2 Flurrer compurarion resulrs of recrangular plare

For the rectangular plate，the relative difference of natural frequencies between the two EPM methods and the FEM method is within 5%，while the EPM method based on FSDT is more accurate than that based on CPT.Therefore，for the EPM calculations，F(xiàn)SDT brings smaller errors.And the higher the order of the natural frequencies is，the more advantages of FSDT there are.Also，mode shapes obtained by the present method and the FEM method are almost the same.According to Table 2，flutter frequencies and flutter speed obtained by the EPM method agree well with those obtained by the FEM method.It shows good accuracy and efficiency of the present method.

3.2 EPM of wing

Fig.4 V-g curves and V-f curves of rectangular plate

An example of a wing with medium aspect ratio is selected to verify the present method.The wing is made of metals and it has a two-spar structure and 11 ribs with NACA0012 airfoil. The root of the wing is fixed in the analysis process.Its specific geometric shape and the distribution of ribs and beams are shown in Fig.5.

Fig.5 Finite element model of wing

The comparisons of natural frequencies between EPM and FEM is shown in Table 3.Comparisons are also made in Fig.6 between the vibrational mode shapes as obtained by the present method and those by the FEA software MSC/ NASTRAN employing hundreds of shell elements（CQUAD4）.The natural frequencies，as obtained using EPM，agree well with those obtained using FEM.Also，the relative differences of the natural frequencies between the two methods for the first three modes are within 3.67%.In addition，the vibrational mode shapes of EPM and FEM are also in good agreement.All results above show that the results obtained using the present method are rather accurate for complex aircraft wing structures.

Table 3 Narural frequencies of wing

Fig.6 Vibrational mode shapes of wing

After separating flutter computations of the wing structures using FEM and EPM，flutter characteristics of the two different methods are easily obtained.Comparisons between the flutter speed as well as the flutter frequencies as obtained by the two methods are shown in Table 4.Comparisons between V-g curves and V-f curves are also shown in Fig.7，respectively.

According to Table 4 and Fig.7，it is clearly shown that the relative differences between flutter frequencies obtained using FEM and those obtained using EPM is 7.22%，while the relative differences of flutter speed is 7.04%，both of which are acceptable.Both V-g curves and V-f curves shown in Fig.7，as obtained using EPM，agree well with those obtained using FEM.Also，the computation time of the FEM method using MSC/NASTRAN with a computer（i5CPU，4 GB system memory）is about 6 s，while that of the EPM method using MATLAB code is less than2 s.The reason of the decrease of the computation is that the degrees of freedom of the model are largely reduced.Especially for the complex FEM with more tens of thousands of nodes，the computation advantage of EPM is more obviously.In general，the results obtained in flutter analysis using the present EPM method are accurate for complex aircraft wing structures，and the computational efficiency is better.

Table 4 Flurrer compurarion resulrs of wing

Fig.7 V-g curves and V-f curves of wing

4 Conclusions

A method of equivalent-plate models simplification has been developed，on the assumption that the wing behaves like a plate whose deformation can be modeled using FSDT of Reissner-Mindlin.By this，vibration and flutter analysis of aircraft wing structures composed of skins，spars，and ribs using the present method become available.The suggested method can be applied to metallic wings at present，and it also has the promising applications in composite wings.According to the results of the examples using the two methods，we can conclude：

（1）It turns out that the EPM method with the Legendre polynomials as the trial functions has achieved the goal of accuracy，fast modeling and effective flutter analysis，on condition that the relative thickness of the wing structures is relative small.And for the EPM calculations and analysis，F(xiàn)SDT is more accurate than CPT. Therefore，the EPM method can play an important role in preliminary design phases of aircraft design.

（2）Vibration and flutter analysis using equivalent plate model is available for complex wing structures.Comparisons for a series of examples with commercial FEA calculations have shown the accuracy and reliability of the method for design purposes.The element information of MSC/NASTRAN can be imported and processed to accomplish the equivalent modeling by the present method，which shows the interface established between EPM and FEM works well.

（3）When building an equivalent model of a simple structure，we will get a rather high accuracy.However，as the complexity increases，the accuracy of the analysis results will reduce. Therefore，analysis time of EPM is also an important factor in the modeling，which means that the EPM method should take both accuracy and efficiency into account.

［1］ Giles G L.Equivalent plate model for conceptual design of aircraft wing structures［R］.AIAA 95-3945，1995.

［2］ Chen Guibin，Zou Congqing，Yang Chao.Elements of aeroelastic dynamics（Version 2）［M］.Beijing：Beihang University Press，2010.（in Chinese）

［3］ Giles G L.Further generalization of an equivalent plate representation for aircraft structural analysis［J］.Journal of Aircraft，1989，26（1）：67-74.

［4］ Giles G L.Equivalent plate analysis of aircraft wing box structures with general planform geometry［J］. Journal of Airxraft，1986，23（11）：859-864.

［5］ Tizzi S.Numerical procedure for the dynamic analysis of three-dimensional aeronautical structures［J］. Journal of Aircraft，1997，34（1）：120-130.

［6］ Tizzi S.Improvement of a numerical procedure for the dynamic analysis of aircraft structures［J］.Journal of Aircraft，2000，37（1）：144-154.

［7］ Livne E.Equivalent plate structural modeling for wing shape optimization including transverse shear［J］.AIAA Journal，1994，32（6）：1278-1288.

［8］ Livne E，Navarro I.Nonlinear equivalent plate modeling of wing-box structures［J］.Journal of Aircraft，1999，36（5）：851-865.

［9］ Demasi L，Livne E.Structural Ritz-based simplepolynomial nonlinear equivalent plate approach：An assessment［J］.Journal of Aircraft，2006，43（6）：1685-1697.

［10］Kapania R K，Lovejoy A E.Free vibration of thick generally laminated cantilever quadrilateral Plates［J］.AIAA Journal，1996，34（7）：1474-1486.

［11］Kapania R K，Liu Youhua.Static and vibration analyses of general wing structures using equivalent-plate models［J］.AIAA Journal，2000，38（7）：1269-1277.

［12］Ranjan Vepa.Aeroelastic analysis of wing structures using equivalent plate models［J］.AIAA Journal，2008，46（5）：1216-1224.

［13］Krishnamurthy T.Frequencies and flutter speed estimation for damaged aircraft wing using scaled equivalent plate analysis［R］.AIAA 2010-2769，2010.

［14］Krishnamurthy T，Tsai F J.Static and dynamic structural response of an aircraft wing with damage using equivalent plate analysis［R］.AIAA 2008-1967，2008.

［15］Na Youngho，Shin Sangjoon.Equivalent-plate analysis for a composite wing with a control surface［J］. Journal of Aircraft，2013，50（3）：853-862.

［16］Mindlin R D.Influence of rotatory interia and shear on flexural motions of isotropic，elastic plates［J］. Journal of Applied Machanics，1951（18）：31-38.

［17］Reissner E.The effect of transverse shear deformation on the bending of elastic plates［J］.Journal of Applied Mechanics，1945（12）：69-77.

［18］Wang Xucheng.Fundamentals and numerical methods（Version 2）［M］.Beijing：Tsinghua University Press，1996.（in Chinese）

（Executive Editor：Xu Chengting）

V211.47Documenr code：AArricle ID：1005-1120（2015）05-0508-09

*Corresponding aurhor：Xie Changchuan，Lecturer，E-mail：xiechangc@163.com.

How ro cire rhis arricle：Tang Jian，Xie Changchuan，Yang Chao.Flutter analysis of aircraft wing using equivalent-plate models with orthogonal polynomials［J］.Trans.Nanjing U.Aero.Astro.，2015，32（5）：508-516.

http：//dx.doi.org/10.16356/j.1005-1120.2015.05.508