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.
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 th body is thuswhere , , and are the Cartesian coordinates and , , , and 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 equationsDifferentiating these equations twice with respect to time , one obtains the kinematic acceleration equationswhere . The dynamic equations of motion, for example, derived from Lagrange method, are, for the th body, written aswith mass matrix, vector of Lagrange multipliers, vector of generalized applied forces, and 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 obtainedThese 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 (, , and ) are assigned to it. The relative forces and moments exchanged between two connected bodies are calculated aswhere and are the relative displacements and velocities and are the linear stiffness and damping matrices. The stiffness matrix is a 6 × 6 symmetric matrix given byand the damping is taken as being proportional to the diagonal of the stiffness matrix by a damping factor ; that is, .
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 bywhere represents drag or lift and the aerodynamic pitching moment byUsing , , and to indicate the relative airflow velocities in body axes for each strip, the local angle of attack is calculated asand 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 is performed, such thatwhere is the complex reduced frequency and 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 about 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 , and to the time derivative of the local angle of attack, the matrix . This latter term is computed using finite differences aswhere is a value of reduced frequency sufficiently close to zero. In deriving (12), the following properties of have been used:(i)The matrix is assumed to be analytic and, as a result, satisfies the Cauchy-Riemann equations so that .(ii)The real part of is an even function of so that .(iii)The imaginary part of is an odd function of so that .
Equivalent sectional lift coefficients and 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 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.
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.
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 and torsional stiffness ) 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 and torsional stiffness 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 and moments of area , , and . 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, 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).
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 -, -, and -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 aswhere is the distance along the beam axis between two consecutive sections, represents the rate of twist of the section considered, and and represent, respectively, the relative rotations (difference between the rotations of two consecutive sections) about the cross-section axes (Figure 3). The subscript is either 1 or 2 and refers to the two unit tip moment load cases, respectively, Nm and 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 ( set) and those omitted ( set), such that Then perform the following transformation on the stiffness matrix:where the transformation matrix is given byAs 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 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.
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.
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.
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%.
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.
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 coordinate, together with its lateral, , and vertical, , 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.
(a)
(b)
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 -axis. The nonlinear analysis predicts a lower vertical bending moment , −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.
The main differences between the linear and nonlinear results occur however in the in-plane loads, shear and moment , and axial force , 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 and and the ailerons control derivative , with 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.
It can be noted how the lift curve slope 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, 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 . 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 (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.
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.
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.
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.
Competing Interests
The authors declare that they have no competing interests.
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.