Flight Loads Prediction of High Aspect Ratio Wing Aircraft Using Multibody Dynamics

Abstract

A framework based on multibody dynamics has been developed for the static and dynamic aeroelastic analyses of flexible high aspect ratio wing aircraft subject to structural geometric nonlinearities. Multibody dynamics allows kinematic nonlinearities and nonlinear relationships in the forces definition and is an efficient and promising methodology to model high aspect ratio wings, which are known to be prone to structural nonlinear effects because of the high deflections in flight. The multibody dynamics framework developed employs quasi-steady aerodynamics strip theory and discretizes the wing as a series of rigid bodies interconnected by beam elements, representative of the stiffness distribution, which can undergo arbitrarily large displacements and rotations. The method is applied to a flexible high aspect ratio wing commercial aircraft and both trim and gust response analyses are performed in order to calculate flight loads. These results are then compared to those obtained with the standard linear aeroelastic approach provided by the Finite Element Solver Nastran. Nonlinear effects come into play mainly because of the need of taking into account the large deflections of the wing for flight loads computation and of considering the aerodynamic forces as follower forces.

[ProQuest: [...] denotes non US-ASCII text; see PDF]

Michele Castellani 1 and Jonathan E. Cooper 1 and Yves Lemmens 2

Academic Editor:Kenneth M. Sobel

1, Department of Aerospace Engineering, Faculty of Engineering, University of Bristol, Bristol BS8 1TR, UK

2, Aerospace Centre of Competence, Siemens PLM Software, 3001 Leuven, Belgium

Received 10 August 2016; Revised 4 November 2016; Accepted 29 November 2016; 27 December 2016

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In recent years there has been a strong push in the aviation world towards the reduction of fuel consumption and the design of ecoefficient aircraft. Many research initiatives are currently addressed to investigate and develop design solutions that would lead to achieve these goals. The improvement of aerodynamic performance is at the forefront of these efforts and one of the most promising concepts being sought is the design of high aspect ratio wings. High aspect ratio wings can lead to significant fuel savings due to the reduction in induced drag. For future designs, a number of high aspect ratio wing configurations are currently being considered and both Airbus [1] and Boeing [2] have published their own concepts.

High aspect ratio wings nevertheless suffer from certain structural drawbacks. Due to the large span, the bending moment increases, resulting in higher structural weight. In order to achieve an effective performance benefit, a lightweight wing design is needed, which in turn leads to very flexible structures, where geometric nonlinearities due to large displacements cannot be neglected anymore. The greater flexibility and lower structural natural frequencies could also result in a strong coupling between structural dynamics and rigid body (flight mechanics) modes leading to undesirable effects on the handling qualities.

The move away from a linear behavior means that a nonconventional approach needs to be taken for the loads and aeroelastic analysis, in order to deal with geometric nonlinearities, and also the nonlinear aerodynamics and flight mechanics characteristics [3]. The ability to predict accurately limit loads, including these nonlinear effects, from the conceptual design phase onwards is paramount in achieving an optimized structural sizing and eventually reaching success with these configurations.

A great deal of work has considered the aeroelasticity of very flexible aircraft [4-11]. Most approaches have used nonlinear beam models coupled to aerodynamic models ranging from strip theory to unsteady vortex lattice method and CFD. However, less focus has instead been devoted to the use of multibody simulation for the modelling of high aspect ratio wings, the two most relevant pieces of work being those presented by Krüger [7] and Zhao and Ren [9]. Recently, Castellani et al. [12] developed two nonlinear methodologies, based, respectively, on nonlinear Finite Element Method (FEM) and multibody dynamics, for the static aeroelastic trim analyses including structural nonlinearities and applied these to a very flexible High-Altitude Long Endurance Unmanned Aerial Vehicle test case.

In this work, a framework based upon multibody dynamics is developed for the static and dynamic aeroelastic analyses of high aspect ratio wing aircraft including structural nonlinearities. The nonlinearities considered are the so-called geometric nonlinearities, arising because of the large deflections that a flexible high aspect ratio wing undergoes when loaded. Following this assumption, a further source of nonlinearity that must be introduced is the follower nature of the aerodynamic forces.

The studies performed are limited to structures undergoing large displacements, but small strains, so that the material constitutive law is still linear, and to attached subsonic flow, so that transonic and stall effects can be neglected.

The focus of this paper is on static and dynamic flight loads prediction, in accordance with the loads requirements set by airworthiness regulations (EASA CS-25 and FAR-25). Most of the research efforts dealing with structural nonlinearities in aeroelasticity have focused on the prediction of aeroelastic and flight dynamics instabilities; less focus has been instead devoted to the impact of geometric nonlinearities on flight loads and studies on this topic have been performed, for example, by Garcia [6] and De Breuker et al. [13]. There is therefore a need in the industry to develop tools and methodologies able to take into account these effects and assess their importance in the design of future high aspect ratio wing aircraft.

2. Aeroelastic Modelling in Multibody Dynamics

Multibody dynamics simulation is a convenient tool capable of simulating multiphysics systems with arbitrary types of nonlinearities and both rigid and flexible components [14]. In the fixed-wing aeroelasticity field, it has been employed for the trim and simulation of manoeuvring flexible aircraft coupled with aerodynamic methods of various levels of fidelity [7, 15].

For the nonlinear aeroelasticity of very flexible aircraft, there have been applications of multibody simulation by Krüger [7] and Zhao and Ren [9], respectively, for the study of the flight mechanics stability of a HALE configuration and for the aeroelastic stability analysis and flight control in manoeuvres of a UAV-like flexible aircraft.

Multibody dynamics allows for arbitrary large displacements and rotations, generic force definition (follower and nonfollower) and inherent coupling between large rigid body motion, linked to flight mechanics, and elastic deformation, without the need of developing dedicated formulations. These are distinct advantages that make multibody dynamics attractive for the analysis of high aspect ratio wings including structurally nonlinear effects.

The multibody software employed for this work is LMS Virtual.Lab Motion v.13.1, a Commercial Off The Shelf (COTS) software developed by Siemens PLM [16].

In the following the equations of motion of a multibody system are briefly outlined (for more details see Shabana [14]). Each body is described by a set of Cartesian coordinates, identifying the location of its centre of gravity in the global reference frame. The vector of the generalized coordinates of the ith body is thus [figure omitted; refer to PDF] where x, y, and z are the Cartesian coordinates and _e0 , _e1 , _e2 , and _e3 the (redundant) Euler parameters used to describe the orientation of the body and to avoid the singularity occurring with other representations, for example, Euler angles.

The bodies in the system are connected together by joints and kinematic relationships, which are expressed as general nonlinear algebraic constraint equations [figure omitted; refer to PDF] Differentiating these equations twice with respect to time t, one obtains the kinematic acceleration equations [figure omitted; refer to PDF] where _Qd =-_Ctt -_{(Cq q )q} q -2_Cqt q . The dynamic equations of motion, for example, derived from Lagrange method, are, for the ith body, written as [figure omitted; refer to PDF] with _Mi mass matrix, _λi vector of Lagrange multipliers, _Qe,i vector of generalized applied forces, and _Qv,i vector of velocity dependent terms. Adding the kinematic relationships to the equations of motion, a system of nonlinear Differential Algebraic Equations (DAE) describing the kinematics and dynamics of a multibody system is obtained [figure omitted; refer to PDF] These equations are nonlinear, as the matrices are a function of the vector of generalized coordinates itself, and are solved using a Backward Differentiation Formula integrator.

The bodies can be considered either as rigid or as flexible. The most common approach to model flexibility is a modal representation based on Component Mode Synthesis from FEM [14], which adds to the generalized coordinates the modal participation factors of each mode used to represent a body's flexibility. This approach however limits the applicability to linear structures with small elastic displacements. Formulations based on nonlinear FE beams [17] and generic nonlinear FEM elements [18] have been also proposed to this purpose.

The work presented herein employs a simpler, yet efficient, approach to model a flexible wing with arbitrary large elastic displacements. It is based on the discretization of the wing by a series of rigid bodies, to which inertial properties are assigned, interconnected by beam force elements, representing the stiffness distribution. The CG of each body can have any arbitrary offset with respect to the elastic axis chordwise location. In the literature this modelling technique has been referred to as the Finite Segment approach [19] and has been successfully used for very flexible aircraft [7, 9]. Since the multibody formulation allows arbitrarily large rigid body motion, each wing section can undergo large displacements and rotations, and the ensuing internal forces are determined based on this displacement field. Each multibody beam element connects two consecutive rigid bodies and has a stiffness matrix derived from FE linear 6-degree-of-freedom (DOF) beam theory and the usual cross-sectional properties (EA, EI, and GJ) are assigned to it. The relative forces and moments _Fel exchanged between two connected bodies are calculated as [figure omitted; refer to PDF] where x and x are the relative displacements and velocities K and D are the linear stiffness and damping matrices. The stiffness matrix is a 6 × 6 symmetric matrix given by [figure omitted; refer to PDF] and the damping is taken as being proportional to the diagonal of the stiffness matrix by a damping factor ξ; that is, D=ξ·diag[...](K).

The aerodynamic model is based on quasi-steady strip theory. Though more simplistic than higher-fidelity methods, this approach is suitable and still accurate for high aspect ratio wings. Besides, the assumption of quasi-steady aerodynamics is deemed acceptable because the first natural frequencies of a flexible high aspect ratio wing aircraft are generally low (refer to Table 2 for the natural frequencies of the test case considered in this work) and, considering the speeds forming the typical flight envelope of a commercial transport aircraft, the resulting reduced frequencies are also low.

To further support this choice, strip theory can be straightforwardly integrated with the wing Finite Segment representation because no interpolation process is required between the aerodynamic and structural meshes, and the aerodynamic forces and moment are in fact applied at the aerodynamic centre of each rigid body, which represents a strip.

The aerodynamic forces on each strip are given by [figure omitted; refer to PDF] where j represents drag _CD or lift _CL and the aerodynamic pitching moment by [figure omitted; refer to PDF] Using u, v, and w to indicate the relative airflow velocities in body axes for each strip, the local angle of attack α is calculated as [figure omitted; refer to PDF] and includes all the contributions due to the aircraft states (aircraft angle of attack, sideslip, and angular rates) and to the elastic deformation of each section.

The quasi-steady aerodynamics stems from two contributions: the first being the inclusion in the sectional α of the kinematic boundary conditions due to the heave and pitch motion of the wing section and the second being the terms proportional to the angle of attack time derivative. As pointed out by Dowell [20], there is ambiguity in the definition of quasi-steady approximation; in this work, it is assumed that the quasi-steady approximation is an expansion in reduced frequency of the unsteady aerodynamics for sinusoidal motion truncated to the first power of frequency, which in time domain corresponds to the first-time derivative, represented by the term proportional to α .

In order to compare the results of the multibody nonlinear approach to the linear FEM, which employs linear DLM aerodynamics, limiting the sources of discrepancies between the methodologies, equivalent strip theory coefficients are derived from the DLM aerodynamic matrix.

In the light of the quasi-steady approximation, an expansion, truncated to the first derivative, of the DLM unsteady aerodynamic matrix about zero reduced frequency k is performed, such that [figure omitted; refer to PDF] where p=g+ik=s_la /_V∞ is the complex reduced frequency and _Qjα indicates the complex matrix, tabulated versus a set of reduced frequencies, relating aerodynamic forces on aerodynamic panels to a change in the local downwash, that is, local angle of attack. The expansion of _Qjα about k=0 delivers real matrices which, in the time domain, relates the aerodynamic force on each panel, acting along the panel normal as per the DLM assumption, to the local angle of attack, the matrix _Qjα (0), and to the time derivative of the local angle of attack, the matrix ^{Qjα[variant prime]} (0). This latter term is computed using finite differences as [figure omitted; refer to PDF] where k- is a value of reduced frequency sufficiently close to zero. In deriving (12), the following properties of _Qjα have been used:

(i) The matrix _Qjα (p) is assumed to be analytic and, as a result, satisfies the Cauchy-Riemann equations so that (∂/∂p)[_Qjα ]=(∂/∂(ik))[_Qjα ].

(ii) The real part of _Qjα is an even function of k so that (∂/∂k)Re[...][_Qjα (0)]=0.

(iii): The imaginary part of _Qjα is an odd function of k so that Im[...](_Qjα (k))=-Im[...](_Qjα (-k)).

Equivalent sectional lift coefficients _CL,α and _CL,α are derived from this expansion by summing the matrix terms corresponding to a strip of chordwise panels along the wing span.

By computing the coefficients from a 3D aerodynamic method such as DLM it is possible to correct strip theory for the sweep angle and tip loss effects. In addition to the lift coefficients, a constant drag coefficient _CD0 is assigned to each strip, representing the airfoil viscous drag.

3. High Aspect Ratio Wing Aircraft Model

The multibody framework presented in this paper has been applied to a high aspect ratio wing aircraft representative of a future concept of a narrow-body commercial transport aircraft, similar to the Boeing Sugar Volt configuration [2] and depicted in Figure 1. It features a high wing with moderate sweep angle, two wing-mounted engines, and a conventional aluminium construction. The main geometric and inertial properties of the test case aircraft are given in Table 1.

Table 1: Main properties of the high aspect ratio wing aircraft.

Property	Value
MTOW	75012 kg
Wingspan	50 m
Mean aerodynamic chord	3.174 m
Wing area	138.42 m2
AR	18.1
Length	42.2 m
Quarter chord sweep	18.1°
Taper ratio	0.201

Table 2: Lowest natural frequencies of free-free aircraft.

Mode description	Frequency [Hz]
1st symm. bending	0.836
1st antisymm. bending	1.238
1st antisymm. bending + VT bending	1.792
2nd symm. bending	2.059
1st antisymm. bending + VT torsion/HT bending	2.235
2nd antisymm. bending	2.757
1st antisymm. mixed bending + VT torsion/HT bending	3.281
1st symm. mixed bending + fuselage/HT bending	3.640
3rd symm. bending + torsion	3.710
1st in-plane bending + fuselage/HT bending	4.019
3rd antisymm. bending + torsion	4.848

Figure 1: High aspect ratio wing aircraft.

[figure omitted; refer to PDF]

In the first step, a Finite Element (FE) model of this aircraft, based on the FE package NX Nastran, has been created. This model forms the basis of the multibody aeroelastic model derivation and of the comparison between linear FEM aeroelastic analyses and nonlinear aeroelastic analyses by multibody dynamics.

The model includes both structural and aerodynamic meshes and has been created with the free software NeoCASS (for more details refer to [21]).

The structural model is a hybrid stick-shell model (shown in Figure 2), where the fuselage, horizontal tail, and vertical tail are represented by beam elements and the wing box is instead a 3D model, with shell elements for the skin, ribs, and spar webs and beam elements for the stringers and spar caps. The structural mass is directly represented by the density on the finite elements, whereas engines, landing gears, systems, payload, and fuel are introduced as concentrated masses. The lowest natural frequencies of the free-free aircraft are presented in Table 2.

Figure 2: 3D FE structural mesh of the high aspect ratio wing aircraft.

[figure omitted; refer to PDF]

Regarding the aerodynamic model, a flat plate mesh of the lifting surfaces, such as that required by vortex lattice and Doublet Lattice Method (DLM), has been created. DLM has been employed to generate the steady and unsteady aerodynamics matrices needed to build the multibody strip theory aerodynamic model, as described previously.

4. Stick Model Development from 3D FE Model

As previously stated, the multibody structural model of the wing consists of rigid bodies interconnected by beam elements, which are defined by a 6-DOF stiffness matrix. Starting from the 3D FE wing box model shown in Figure 2, an equivalent stick model in the multibody environment is then generated.

The reduction of a 3D FE model (3D FEM) to a stick model has been the subject of various investigations in the past. Since the beam stiffness matrix is fully defined by cross-sectional properties, these can be calculated analytically from a built-up 3D FE or CAD model of a wing box, as the geometric (wing box height and width, skin and spar webs thickness, stringers, and spar caps area) and material properties are known. Bindolino et al. [22] applied cross-sectional analysis to estimate the cross-sectional properties of a wing box in the framework of a multilevel structural optimization. For more complex composite sections, where coupling terms between all the deformation components become important, specific cross-sectional analyses tools have been developed by Giavotto et al. [23], ANBA, and by Cesnik and Hodges [24], VABS, and mainly applied to helicopter blades and wind turbine blades.

A second approach consists of identifying the classical cross-sectional stiffness (bending stiffness EI and torsional stiffness GJ) by loading the wing, assuming cantilever boundary condition, with unitary load cases and working out the stiffness, at each section of interest along the span, from the relative displacements and rotations. Singh and Nichols [25] proposed a procedure to derive the elastic axis and the equivalent stiffness of a beam model from a built-up wing box Nastran model. Their procedure consists in applying unit moments at the free end of the cantilevered wing box structure and estimating the bending stiffness EI and torsional stiffness GJ from the relative rotations between reference points along the wing box axis. Recently, Jones and Cesnik [26] applied this technique to develop a nonlinear beam model of the X-56A Multi-Utility Demonstrator from a Nastran FE model to perform aeroelastic analysis. Similarly, Malcolm and Laird [27] extracted equivalent beam properties from an ANSYS FE model of a wind turbine blade by applying unit loads a the tip and processing the nodal displacements in order to obtain a 6 × 6 stiffness matrix at each blade section; these have been subsequently used to generate a multibody model in MSC Adams of a wind turbine. Elsayed et al. [28] reviewed the most common methodologies employed in the industry to generate stick models from 3D ones and proposed an improved procedure based on applying unit tip moments, deriving bending rotations from displacements and eventually bending and torsional stiffness.

Another category of methodologies includes mathematical reduction techniques, such as Guyan reduction [29] and Improved Reduced System [30]. In both methods, a set of master nodes of the FE model is selected and the mass and stiffness matrices reduced to this. Guyan reduction, also known as static condensation, is well established in the aerospace industry and the reduced equations are developed using only the stiffness matrix, leading to an exact reduction of the stiffness, but only an approximate reduction of the mass matrix. IRS is an extension of the former methodology that includes mass effects in the development of the system reduction transformation matrix. Wang et al. [31] proposed a procedure to identify a geometric nonlinear modal-based 1D intrinsic beam model by the application of Guyan reduction to a 3D FEM onto a small set of nodes along the axis of slender beam-like structures. A potential drawback of the reduction methodologies, compared to those previously described, is that the reduced mass and stiffness matrices are fully populated and lose the link to the physical distribution of the usual cross-sectional stiffness and mass.

In the present work, three approaches are investigated to build a multibody equivalent stick model from the 3D FEM of the high aspect ratio wing aircraft. These are as follows:

(A) Cross-sectional analysis.

(B) Stiffness identification by unitary loadings.

(C) Guyan reduction.

4.1. Method A: Cross-Sectional Analysis

The first method, cross-sectional analysis, applies the classical formula from thin walled structural analysis to calculate the sectional area A and moments of area _I1 , _I2 , and J. Knowing the geometry, thickness, and material properties of the structural elements and under the assumption of material isotropy, valid since the wing considered has a metallic construction, this procedure is straightforward and delivers the axial, bending, and torsional stiffness characterizing the beam at each span section. A simplified cross-section representation of the wing box is assumed, shown in Figure 3, made up of upper and lower skin, N stringers, and spar caps and spar webs. With reference to Figure 3, the formulas employed are presented hereafter. The sectional properties obtained are then directly input into the multibody beam stiffness matrices (7). [figure omitted; refer to PDF]

Figure 3: Wing box cross-section.

[figure omitted; refer to PDF]

4.2. Method B: Stiffness Identification by Unitary Loadings

The second method follows the procedure outlined by Elsayed et al. [28]. The wing box is clamped at the root and unit tip moments along the wing reference x-, y-, and z-axes are applied independently and linear static analyses are performed. The stiffness properties are extracted at each wing box section corresponding to the locations of the multibody beam elements, which have been described previously, from the relative rotations of each cross-section. In order to retrieve these, interpolation elements (RBE3) are introduced at the multibody beam locations and attached to the surrounding nodes lying on a cross-section. The interpolation element provides displacements and rotations of a dependent node by averaging the degrees of freedom to which it is connected.

Taking the well-known relationships between load and displacement/rotation from the beam theory [32], for each section, the bending and torsional stiffness properties can be easily calculated as [figure omitted; refer to PDF] where l is the distance along the beam axis between two consecutive sections, [straight phi] represents the rate of twist of the section considered, and _θ1,l and _θ2,l represent, respectively, the relative rotations (difference between the rotations of two consecutive sections) about the cross-section axes _x1 -_x2 (Figure 3). The subscript l is either 1 or 2 and refers to the two unit tip moment load cases, respectively, _{[M1 ,M2 ]1} =(1,0) Nm and _{[M1 ,M2 ]2} =(0,1) Nm, used to identify the stiffness distributions. The above equations assume that the vertical and in-plane bending are coupled whereas the torsion (i.e., rotation about the beam reference axis) is independent.

The sectional properties obtained are then directly input into the multibody beam stiffness matrices (7). The sectional area is instead directly calculated from the cross-section geometry (13).

4.3. Method C: Guyan Reduction

Guyan reduction is chosen as the technique to derive the reduced stiffness matrix from the 3D FEM. Regarding the mass distribution of the multibody model, both the structural (distributed) and concentrated masses are first discretized along the span, considering wing segments between each reference point of the equivalent multibody stick model, and subsequently lumped to these locations, taking into account proper CG offsets and moments and products of inertia. In this way, there is no need to use the mass matrix obtained by Guyan reduction, which is known to be inaccurate [30].

To reduce the size of the 3D FEM, the degrees of freedom are first divided into two subsets, those retained in the reduced model (a set) and those omitted (o set), such that [figure omitted; refer to PDF] Then perform the following transformation on the stiffness matrix: [figure omitted; refer to PDF] where the transformation matrix T is given by [figure omitted; refer to PDF] As previously mentioned, this methodology is purely a mathematical reduction and the reduced stiffness matrix does not provide information about the physical distribution of the cross-sectional properties along the wingspan. However, recalling the way a beam is represented in the multibody environment, (7), it is possible to assign to each of the beam elements along the wing a full 6 × 6 stiffness matrix. The coefficients of these matrices are indeed taken from the reduced stiffness matrix _Ka of the whole model by extracting the 6 × 6 block diagonal submatrices.

4.4. Comparison of the Methodologies

Three stick multibody models of the high aspect ratio wing aircraft have been generated from the 3D FEM employing the three presented methodologies. Nonlinear static analyses have been carried out in Nastran and in the multibody environment to validate the structural modelling. For this validation, only the RHS wing has been considered, clamping it at the root and applying two different loading conditions: a 50000 N tip force and 2.5 g trim loads, both aerodynamic and inertial, along the wing; these have been obtained from a linear aeroelastic trim analysis performed in Nastran using the 3D FEM. Figures 4 and 5 present the wing deflected shape for these two load cases. For reference, the linear solutions obtained with the 3D FEM in Nastran are also reported. For the tip force load case, little difference is evident between the linear and nonlinear 3D FEM results; in the nonlinear solution there is, as expected, a wing lateral shortening. Regarding the multibody results, the model obtained by the stiffness identification technique is the one delivering the most accurate results with respect to the reference 3D FEM. It is interesting to note that the multibody model obtained by cross-sectional analysis, despite being nonlinear, is still overpredicting the vertical displacement. This result suggests that the cross-sectional analysis approach, which is based on approximate formulas, is not accurate for complex geometries and structural layouts.

Figure 4: Wing deflected shape, tip force 50000 N, and comparison of models.

[figure omitted; refer to PDF]

Figure 5: Wing deflected shape, 2.5 g linear aeroelastic trim loads, and comparison of models.

[figure omitted; refer to PDF]

The wing deflected shape for the 2.5 g trim load condition, in Figure 5, shows a deviation between the linear and nonlinear 3D FEM results but still confirms the considerations presented above. The stick model generated by cross-sectional analysis overestimates the displacements and bending curvature. Similarly, the model generated by Guyan reduction shows higher displacements than the reference solution for both load cases. The best agreement with the reference nonlinear 3D FEM is again obtained with the multibody model generated by the stiffness identification procedure.

Following the static analyses, a normal modes analysis of the 3D FEM RHS wing with free-free boundary conditions has been carried out and the natural frequencies compared to those obtained with the three multibody models. In the multibody environment natural frequencies and mode shapes are computed by performing a linearization of the equations, through the finite differences method, about the undeformed configuration.

Table 3 reports the first five natural frequencies and the relative errors between the 3D FEM and the three multibody models. These results confirm the trend shown by the static analyses. The cross-sectional analysis method exhibits the largest differences, underpredicting both the bending and torsional natural frequencies, which confirms that the stiffness is underestimated. Likewise, Guyan reduction, though more accurate, leads to lower natural frequencies whereas the model created by stiffness identification shows the best agreement. Anyhow, for all the three methods the greatest discrepancy is in the 1st torsional frequency, which is to be expected since the torsional behavior of a complex built-up structure is difficult to predict accurately with a stick model.

Table 3: Comparison of the first five natural frequencies of the free-free RHS wing.

Mode description	3D FEM	Method A	Method B	Method C
Frequency [Hz]	Frequency [Hz]	Error [%]	Frequency [Hz]	Error [%]	Frequency [Hz]	Error [%]
1st bending	2.421	2.363	-2.40%	2.412	-0.37%	2.405	-0.66%
2nd bending	5.593	5.405	-3.36%	5.553	-0.72%	5.527	-1.18%
1st in-plane bending	7.694	7.417	-3.60%	7.628	-0.86%	7.599	-1.23%
3rd bending	10.204	9.779	-4.17%	10.109	-0.93%	10.064	-1.37%
1st torsion	15.350	13.062	-14.91%	14.536	-5.30%	13.609	-11.34%

In the light of the structural validation of the three multibody models, the equivalent stick model generated by stiffness identification (Method B) is deemed the most accurate and it is the one selected for the nonlinear static and dynamic aeroelastic analyses presented in the following. The model is shown in Figure 6.

Figure 6: Multibody model of the high aspect ratio wing aircraft.

[figure omitted; refer to PDF]

As a further validation of the multibody structural and inertial modelling selected (Method B), the first prestressed symmetric natural frequencies of the free-free aircraft under increasing trim loads (undeformed/0 g and from 1 g to 2.5 g) are compared in Table 4 to the 3D FEM results. These results confirm the good agreement of the multibody modelling selected and also illustrate that a stiffening effect occurs on the bending frequencies when the wing is loaded and deformed. For instance, in the 2.5 g deformed configuration, the frequency 1st symmetric bending mode increases by 9%.

Table 4: Natural frequencies of the free-free aircraft under trim loads: linear FEM versus multibody (Method B).

Mode description	3D FEM - frequency [Hz]	Multibody - frequency [Hz]
0 g	1 g	1.5 g	2 g	2.5 g	0 g	1 g	1.5 g	2 g	2.5 g
1st symm. bending	0.836	0.851	0.867	0.887	0.910	0.833	0.848	0.864	0.884	0.907
2nd symm. bending	2.059	2.076	2.095	2.118	2.144	2.045	2.059	2.080	2.097	2.123
1st symm. mixed bending + VT torsion/HT bending	3.640	3.739	3.754	3.766	3.775	3.607	3.646	3.649	3.657	3.673
3rd symm. bending + torsion	3.710	3.364	3.204	3.134	3.106	3.625	3.327	3.169	3.091	3.066
1st in-plane bending + fuselage/HT bending	4.019	4.197	4.253	4.292	4.322	3.971	4.147	4.193	4.232	4.257

5. Nonlinear Aeroelastic Trim

Nonlinear aeroelastic trim analyses have been performed with the multibody stick model of the high aspect ratio wing aircraft. The results of such analyses are compared to those obtained by standard linear trim analyses carried out in Nastran (SOL144) using the 3D FEM, with the purpose of highlighting the effects of structural nonlinearities on flight loads prediction.

In the multibody approach, the trim solution is sought by performing a dynamic settling simulation with the implementation of controllers in order to achieve a steady trimmed state. Details of the trim methodology developed are provided by Castellani et al. [12].

Regardless of the structural method employed (FEM, multibody dynamics, Ritz-Raleigh method, etc.) and within the framework of linear aerodynamics, the main differences between a standard linear aeroelastic trim solution and an aeroelastic trim procedure including geometric nonlinearities, such as the one proposed, are as follows:

(i) Large displacements and rotations: high aspect ratio flexible wings undergo large displacements and rotations that cannot be neglected and second-order effects, such as the wing end shortening, affect wing deformation and eventually flight loads.

(ii) Follower force effect: aerodynamic forces, arising from pressure distributions, are inherently follower forces and it is paramount to take into account this change of aerodynamic force orientation for high aspect ratio flexible wings undergoing large displacements and rotations.

(iii): Computation of wing integrated loads based on the deformed shape: since the hypothesis of small displacements is not valid anymore for high aspect ratio flexible wings, the actual deformed shape must be considered to compute the integrated loads along the wing, which are expressed in the local reference frame of each displaced wing section.

The aircraft is trimmed at two load factors, 1 g and 2.5 g, for a flight condition with Mach number 0.60 at 25000 ft. The 2.5 g load factor corresponds to the maximum positive load factor for a large commercial aircraft (as per EASA CS-25) and typically forms part of the critical loads envelope of the inboard and midboard wing [33].

Table 5 reports the trim angle of attack and the computational time resulting from the linear FEM and the multibody trim analyses. This latter predicts slightly higher trim angles of attack, the reason being that, in the nonlinear approach, the follower force effect of the lift is accounted for and thus, as the wing bends upwards, the lift is progressively tilted inboard and its vertical component, the one balancing the weight, is reduced. As a result the angle of attack required to balance the aircraft must be increased compared to a linear solution and the greater is the load factor, since the bending on the wing increases with it.

Table 5: Trim angle of attack and CPU time: linear FEM versus multibody results.

Load factor [g]	Linear trim FEM	Trim multibody
α [deg]	CPU time [s]	α [deg]	CPU time [s]
1	4.91	<1	5.04	103.0
2.5	12.28	<1	13.40	70.3

In order to gain further insights about the effects of structural geometric nonlinearities, the lift distribution is plotted in Figure 7 versus the undeformed (for linear analysis) and the deformed (for nonlinear analysis) wing y coordinate, together with its lateral, _Fy,aero , and vertical, _Fz,aero , components in aircraft body axes. Due to the significant wing bending at 2.5 g, the lift is tilted inboard and generates a lateral force component, whose maximum magnitude reaches 35% of the total lift at the corresponding span station (51% span). As previously mentioned, this force is neglected by the linear analysis. Furthermore, the lift in the nonlinear solution is shifted inboard because of the wing tip shortening, a second-order effect not captured by a linear structural formulation.

Figure 7: Lift distribution (a) and lateral and vertical lift components (b) at 2.5 g trim, linear FEM versus multibody.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

The wing integrated loads at 2.5 g are presented next in Figure 8, showing the forces, and in Figure 9, showing the moments, and are both plotted along the undeformed wing y-axis. The nonlinear analysis predicts a lower vertical bending moment _MX , -10.8% at the root, and this can be explained by noting that the lift in the nonlinear solution (Figure 7) is shifted inboard and acts through a smaller moment arm because of the wing tip shortening. The lateral force component arising from the follower force effect is contributing as well to the bending moment; nevertheless this contribution is not high enough to compensate for the moment arm shortening and lift redistribution. This result is the opposite of what has been obtained by Castellani et al. [12] on a very flexible high aspect ratio unswept wing, where, due to the extreme bending, the lateral lift component acting out-of-plane overcomes the bending moment reduction caused by the wing tip shortening and leads to higher bending moment predicted by the nonlinear analysis.

Figure 8: Wing integrated loads-forces, 2.5 g trim linear FEM versus multibody.

[figure omitted; refer to PDF]

Figure 9: Wing integrated loads-moments, 2.5 g trim linear FEM versus multibody.

[figure omitted; refer to PDF]

The main differences between the linear and nonlinear results occur however in the in-plane loads, shear _FX and moment _MZ , and axial force _FY , as also pointed out by Castellani et al. [12] on a very flexible high aspect ratio unswept wing. The sources of these differences are the aforementioned lateral lift component and the longitudinal forces arising from the rotation of the lift vector, perpendicular to the airspeed, from wind to body axes, and the rotation of the gravity vector from global (Earth-fixed) to body axes through the trim aircraft pitch attitude.

5.1. Flexible Stability Derivatives

In addition to computing flight loads, nonlinear aeroelastic trim analyses at increasing load factor have been performed in order to estimate the stability derivatives of the flexible aircraft taking into account geometric structural nonlinearities and follower force effects.

Stability derivatives are an important measure of the flight dynamics characteristics and handling qualities of an aircraft and can be highly affected by aeroelastic effects; for instance, the most common issue for aircraft with sweptback wings is the loss of aileron effectiveness at high speed that can lead to aileron reversal. Traditionally, aeroelastic trim analyses are performed to compute the stability derivatives including aeroelastic effects; these derivatives depend both on the Mach number (aerodynamic effect) and on the dynamic pressure (aeroelastic effect), as the aerodynamic forces due to structural displacements grow with increasing dynamic pressure. However, for very flexible wings undergoing large deformations, an additional parameter dependence for the stability derivatives comes into play: the load factor or current loaded wing deformed shape. Because of the large displacements, the current deformed wing shape and the current CG location of the deformed aircraft cannot be approximated anymore with the undeformed configuration and this has an impact on the aerodynamics of the wing and the CG location, leading thus to a change in stability derivatives too.

For the present work, multibody trim analyses at increasing load factors (1 g to 2.5 g) and Mach number 0.60 at 25000 ft have been carried out in order to compute the longitudinal stability derivatives _CL,α and _CM,α and the ailerons control derivative _CL,δ , with _CL being the roll moment coefficient and δ the ailerons deflection. The derivatives are calculated by finite differences, perturbing either the angle of attack α or the ailerons deflection δ by a small quantity. The results are presented in Figure 10, normalised to the stability derivatives values obtained by a linear trim analysis in Nastran, for which stability derivatives, in the light of the linear approach, do not depend on the load factor.

Figure 10: Stability derivatives versus load factor computed by nonlinear trim analysis.

[figure omitted; refer to PDF]

It can be noted how the lift curve slope _CL,α decreases with increasing load factor. This is to be expected as the wing progressively bends and a significant component of the total lift is tilted inward, leading to a loss of vertical lift effectiveness. However, _CM,α shows an opposite trend, that is, an increase (in magnitude) versus load factor. The cause of this is the aft shift of the wing aerodynamic centre due to the large bending deformation coupled with the sweptback wing configuration; because of the sweep angle, the outboard wing sections undergo a significant aft motion and their moment arms with respect to the CG increase, leading thus to an increase in _CM,α . This effect overcompensates the aforementioned loss of vertical lift. Regarding the ailerons control derivative, there is a decrease with load factor, but not as sensitive as the two longitudinal derivatives. The ratio at 1 g is already low, about 70%, but this can be attributed to the different aerodynamic modelling of the ailerons between Nastran and the multibody model.

6. Gust Response

Gust response analyses have been performed with the multibody model of the high aspect ratio wing aircraft and compared to the linear gust response of the 3D FEM, performed with the standard Nastran dynamic aeroelastic solution (SOL146) which employs the modal approach. Following a convergence study, modes up to 30 Hz, including the rigid body ones, are retained in the linear gust analysis of the 3D FEM.

Certification requirements specify the discrete gust load cases considering the aircraft in level flight subject to symmetrical vertical gusts with a "1-cosine" velocity profile, having gust gradient H (half of the gust wavelength) and asking for several gust gradients between 30 ft and 350 ft to be investigated in order to determine the critical conditions [34].

The flight condition assumed for the gust analyses is a Mach number of 0.60 at 25000 ft. In the linear approach, the superposition principle is applied and thus the gust response is performed assuming zero initial conditions for all the modal coordinates and modal velocities; that is, the aircraft is in the undeformed configuration. However, in a nonlinear approach the superposition principle is not valid anymore and the actual initial conditions must be considered. Therefore, in the multibody solution, the aircraft is first trimmed at 1 g and then flown into the gust field.

The time history of the incremental load factor for three gust gradients (90 ft, 220 ft, and 350 ft) and upward gust is shown in Figure 11. There is a close agreement between the linear and nonlinear peak load factors, with the multibody peaks being higher for the shorter gusts (+10%). The incremental root bending moment response, presented in Figure 12, is driven by the load factor. The peak load prediction of the linear FEM and multibody approach is similar, except for the shortest gust gradient, where the root bending moment resulting from the multibody analysis is +33% higher. This discrepancy on the shortest gust can be also caused by the adoption of a quasi-steady aerodynamic approximation in the nonlinear analysis, while the linear one is based on the unsteady DLM. The second (negative) peak shows instead more differences, with the multibody values being considerably lower. The negative bending moment peak occurs when the gust excitation has already subsided and it is mainly driven by the free response of the elastic structure, as the wing springs back after having reached the maximum upward bending. Table 4 shows that the wing bending modes experience a stiffening effect when under loads. The linear FEM gust response is based on a fixed modal basis, made up of the normal modes in the undeformed configuration, and therefore misses this effect, which is instead captured by the multibody method, which does not employ the modal approach and does not have such approximation. This stiffening effect could be the cause of a lower negative peak, as the overswing of the wing is limited by an increased stiffness.

Figure 11: Gust incremental load factor, linear FEM versus multibody.

[figure omitted; refer to PDF]

Figure 12: Gust wing root bending moment, linear FEM versus multibody.

[figure omitted; refer to PDF]

The maximum bending moment throughout the time history along the wing is extracted and plotted in Figure 13. Comparing the linear FEM and multibody loads predicted, there is a cross-over point at 50% span, the latter being lower in the inboard wing (between -6% and -9%) and higher in the outboard sections.

Figure 13: Maximum gust bending moment versus span, linear FEM versus multibody.

[figure omitted; refer to PDF]

Finally, the wing root correlated loads plot, bending versus torsion, is shown in Figure 14. It is generated by performing both upward and downward gust analyses and taking into account the initial 1 g loads (in the linear case by superimposition, in the nonlinear case as the actual gust initial condition). The nonlinear (multibody) plot mostly bounds the linear one, except for the maximum (upward) bending moment points, which are higher for the linear analysis.

Figure 14: Correlated loads plot wing root bending versus torsion, linear FEM versus multibody.

[figure omitted; refer to PDF]

Following the results presented, it can be noted that, for the gust response, the effects of the structural nonlinearities and follower forces considered are less predictable and definite than in static aeroelastic analyses. In addition, for both static and dynamic loads, no general conclusion can be drawn about whether a nonlinear approach can be beneficial in terms of delivering less conservative loads. This outcome and the quantitative differences are distinctly configuration dependant; in particular wing geometry (span, sweep angle, dihedral angle, and planform), stiffness, and mass distributions are the most important factors.

7. Conclusions

Multibody dynamics provides a powerful framework to model flexible high aspect ratio wing aircraft and perform aeroelastic analyses including structural nonlinearities.

The multibody approach presented in this paper has been applied to a commercial aircraft featuring an aspect ratio of 18 with the aim of predicting flight loads, both static and dynamic. From a 3D FE representation of the wing box, a multibody equivalent stick model has been generated by identifying the stiffness distribution through the application of unitary loads at the wing tip. This technique has been shown to deliver an excellent agreement with the reference 3D FEM, with the advantage of a significant computational saving due to the reduced size of the resulting stick equivalent model.

A nonlinear trim procedure has been implemented in the multibody environment through PID controllers and static flight loads have been computed and compared to the linear results. This comparison highlights that there are significant differences between a linear and a nonlinear static aeroelastic approach and the reasons have been identified in the need to consider large displacements and rotations of the wing under loads and from the follower force effects of the aerodynamic forces, factors that both are neglected in the typical linear aeroelastic analyses carried out nowadays in the industry. These effects have an impact on the flexible stability derivatives too, which become a function not only of the Mach number and dynamic pressure, but also of the load factor, indicating that flight dynamics characteristics and handling qualities change when the aircraft is manoeuvring.

Finally, gust responses have been performed and the comparison between linear and nonlinear loads predictions has been shown, including correlated loads plot.

The results presented demonstrated that flight loads prediction including structural nonlinearities can deliver significantly different results than the usual linear approach, confirming the need to develop reliable methodologies to take into account these effects. The trend and the quantitative differences between loads predicted with a linear and a nonlinear method are highly dependent on each aircraft configuration considered.

Acknowledgments

This work is supported by the European Commission (EC FP7) under the Marie Curie European Industrial Doctorate Training Network ALPES (Aircraft Loads Prediction using Enhanced Simulation, Grant Agreement no. 607911) and the Royal Academy of Engineering.

References

[1] "The Airbus Concept Plane," http://www.airbus.com/innovation/future-by-airbus/the-concept-plane/the-airbus-concept-plane/

[2] N. Kehayas, "Propulsion system of jet-flapped subsonic civil transport aircraft design,", Journal of Aircraft , vol. 48, no. 2, pp. 697-702, 2011.

[3] C. E. S. Cesnik, R. Palacios, E. Y. Reichenbach, "Reexamined structural design procedures for very flexible aircraft,", Journal of Aircraft , vol. 51, no. 5, pp. 1580-1591, 2014.

[4] D. Tang, E. H. Dowell, "Experimental and theoretical study on aeroelastic response of high-aspect-ratio wings,", AIAA Journal , vol. 39, no. 8, pp. 1430-1441, 2001.

[5] M. J. Patil, D. H. Hodges, C. E. S. Cesnik, "Nonlinear aeroelasticity and flight dynamics of high-altitude long-endurance aircraft,", Journal of Aircraft , vol. 38, no. 1, pp. 88-94, 2001.

[6] J. A. Garcia, "Numerical investigation of nonlinear aeroelastic effects on flexible high-aspect-ratio wings,", Journal of Aircraft , vol. 42, no. 4, pp. 1025-1036, 2005.

[7] W. R. Krüger, "Multibody dynamics for the coupling of aeroelasticity and flight mechanics of highly flexible structures," in Proceedings of the International Forum on Elasticity and Structural Dynamics (IFASD '07), Stockholm, Sweden, June 2007.

[8] R. Palacios, J. Murua, R. Cook, "Structural and aerodynamic models in nonlinear flight dynamics of very flexible aircraft,", AIAA Journal , vol. 48, no. 11, pp. 2648-2659, 2010.

[9] Z. Zhao, G. Ren, "Multibody dynamic approach of flight dynamics and nonlinear aeroelasticity of flexible aircraft,", AIAA Journal , vol. 49, no. 1, pp. 41-54, 2011.

[10] H. Hesse, R. Palacios, "Consistent structural linearisation in flexible-body dynamics with large rigid-body motion,", Computers and Structures , vol. 110-111, pp. 1-14, 2012.

[11] A. Arena, W. Lacarbonara, P. Marzocca, "Nonlinear aeroelastic formulation and postflutter analysis of flexible high-aspect-ratio wings,", Journal of Aircraft , vol. 50, no. 6, pp. 1748-1764, 2013.

[12] M. Castellani, J. E. Cooper, Y. Lemmens, "Nonlinear static aeroelasticity of high-aspect-ratio-wing aircraft by finite element and multibody methods,", Journal of Aircraft , 2016.

[13] R. De Breuker, M. M. Abdalla, Z. Gürdal, "On the effect of geometric nonlinearities on static load alleviation," in Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Palm Springs, Calif, USA, May 2009.

[14] A. Shabana, Dynamics of Multibody Systems , Cambridge University Press, Cambridge, UK, 2010. 3rd.

[15] L. Cavagna, P. Masarati, G. Quaranta, "Coupled multibody/computational fluid dynamics simulation of maneuvering flexible aircraft,", Journal of Aircraft , vol. 48, no. 1, pp. 92-106, 2011.

[16] LMS Virtual. Lab 13 Online Help, Rev 13, June 2014

[17] G. L. Ghiringhelli, P. Masarati, P. Mantegazza, "Multibody implementation of finite volume C beams,", AIAA Journal , vol. 38, no. 1, pp. 131-138, 2000.

[18] O. A. Bauchau, C. L. Bottasso, Y. G. Nikishkov, "Modeling rotorcraft dynamics with finite element multibody procedures,", Mathematical and Computer Modelling , vol. 33, no. 10-11, pp. 1113-1137, 2001.

[19] J. D. Connelly, R. L. Huston, "The dynamics of flexible multibody systems: a finite segment approach--I. Theoretical aspects,", Computers & Structures , vol. 50, no. 2, pp. 255-258, 1994.

[20] E. H. Dowell, A modern course in aeroelasticity , vol. 217, of Solid Mechanics and its Applications, Springer, Switzerland, 2015. 5th.

[21] L. Cavagna, S. Ricci, L. Travaglini, "NeoCASS: an integrated tool for structural sizing, aeroelastic analysis and MDO at conceptual design level," in Proceedings of the AIAA Atmospheric Flight Mechanics Conference, Toronto, Canada, August 2010.

[22] G. Bindolino, G. Ghiringhelli, S. Ricci, M. Terraneo, "Multilevel structural optimization for preliminary wing-box weight estimation,", Journal of Aircraft , vol. 47, no. 2, pp. 475-489, 2010.

[23] V. Giavotto, M. Borri, P. Mantegazza, G. Ghiringhelli, V. Carmaschi, G. C. Maffioli, F. Mussi, "Anisotropic beam theory and applications,", Computers and Structures , vol. 16, no. 1-4, pp. 403-413, 1983.

[24] C. E. S. Cesnik, D. H. Hodges, "VABS: a new concept for composite rotor blade cross-sectional modeling,", Journal of the American Helicopter Society , vol. 42, no. 1, pp. 27-38, 1997.

[25] A. K. Singh, C. W. Nichols, "Derivation of an equivalent beam model from a structural finite element model," in Proceedings of the MSC 1988 World Users Conference, Los Angeles, Calif, USA, 1988.

[26] J. R. Jones, C. E. S. Cesnik, "Nonlinear aeroelastic analysis of the X-56A multi-utility aeroelastic demonstrator," in Proceedings of the AIAA SciTech-15th Dynamics Specialist Conference, San Diego, Calif, USA, January 2016.

[27] D. J. Malcolm, D. L. Laird, "Extraction of equivalent beam properties from blade models,", Wind Energy , vol. 10, no. 2, pp. 135-157, 2007.

[28] M. S. A. Elsayed, R. Sedaghati, M. Abdo, "Accurate stick model development for static analysis of complex aircraft wing-box structures,", AIAA Journal , vol. 47, no. 9, pp. 2063-2075, 2009.

[29] R. J. Guyan, "Reduction of stiffness and mass matrices,", AIAA Journal , vol. 3, no. 2, pp. 380, 1965.

[30] J. O'Callahan, "A procedure for an Improved Reduced System (IRS) model," in Proceedings of the 7th International Modal Analysis Conference, pp. 17-21, Las Vegas, Nev, USA, 1989.

[31] Y. Wang, R. Palacios, A. Wynn, "A method for normal-mode-based model reduction in nonlinear dynamics of slender structures,", Computers and Structures , vol. 159, pp. 26-40, 2015.

[32] T. H. G. Megson, Aircraft Structures for Engineering Students , Butterworth-Heinemann, Oxford, UK, 1999. 3rd.

[33] T. D. Lomax, Structural Loads Analysis for Commercial Transport Aircraft: Theory and Practice , AIAA Education Series, Reston, Va, USA, 1996. 2nd.

[34] J. R. Wright, J. E. Cooper, Introduction to Aircraft Aeroelasticity and Loads , John Wiley & Sons Ltd, Chichester, UK, 2015. 2nd.