Comparison of structural model reduction methods applied to a large-deformation wing box

Abstract

In this paper, the accuracy and practical capabilities of three different reduced-order models (ROMs) are explored: an enhanced implicit condensation and expansion (EnICE) model, a finite element beam model, and a finite volume beam model are compared for their capability to accurately predict the nonlinear structural response of geometrically nonlinear built-up wing structures. This work briefly outlines the different order reduction methods, highlighting the associated assumptions and computational effort. The ROMs are then used to calculate the wing deflection for different representative load cases and these results are compared with the global finite element model (GFEM) predictions when possible. Overall, the ROMs are found to be able to capture the nonlinear GFEM behaviour accurately, but differences are noticed at very large displacements and rotations due to local geometrical effects.

NOMENCLATURE

f

column matrix of generalised loads

${f_{NL}}$

column matrix of nonlinear modal forces

n

number of linear modes

q

column matrix of amplitudes of the linear modes

r

column matrix of amplitudes of the dual modes

s

number of dual modes

C

flexibility matrix

${I_1}$ , ${I_2}$ , ${I_{12}}$

bending stiffness constants

${I_n}\, n \times n$

identity matrix

$\overline{J}$

torsion stiffness constant

K

stiffness matrix of the structure

M

mass matrix of the structure

R

matrix of moment arms of internal loads in the finite volume formulation

T

kinetic energy of the system

U

potential energy of the system

Greek characters:

$\vartheta $

column matrix of internal loads in the finite volume formulation

$\xi $

arc-length along a beam

$\Lambda $

linear modal stiffness matrix

$\Phi $

matrix of linear modes along its columns

$\Psi $

matrix of dual modes along its columns

1.0 INTRODUCTION

Much effort is currently being devoted to reducing the environmental effect of commercial jet aircraft and future designs are heading towards more flexible higher-aspect-ratio wing (HARW) configurations to take advantage of the reduction of induced drag⁽¹⁾. Previous experiences with the design and test of high-aspect-ratio aircraft have indicated the need for accurate modelling in the time domain of nonlinear structural and aerodynamic effects, cited as a key reason for the crash of the NASA Helios experimental aircraft⁽²⁾. Both the structural response and flight dynamics can be significantly affected by large deformations with tip displacements likely to be in the order of 30% or more of the semi-span^(3,4) . The review of Afonso et al. ⁽⁵⁾ compiled different nonlinearities related to high-aspect-ratio configurations and suggested that aeroelastic concerns may be hindering the development of regional jets with longer wings.

In order to support the development of accurate and efficient nonlinear aeroelastic tools, capable of performing time-domain simulations with numerous test cases, the development of robust structural reduced-order models (ROMs) is important. Wing structure global finite element models (GFEMs) with ${10^5}$ to ${10^6}$ degrees-of-freedom are frequently used in industry to accurately calculate the internal load distribution for stress analysis purposes; however, their computational cost limits the amount of design exploration that can be done, a limitation on exploring new configurations such as HARW aircraft. For aeroelastic loads and stability analysis, reduced-order models of GFEMs can be used to cut the computational costs, while still accurately capturing the overall behaviour of the structure. Nonlinear structural ROMs may be used to efficiently analyse the dynamic aeroelastic behaviour of high-aspect-ratio wings undergoing large deflections⁽⁶⁾. ROMs could also be integrated into aeroelastic optimisation frameworks to reduce run times and allow wider design space explorations. Another application of ROMs is the design of modern control strategies, such as model predictive control (MPC)⁽⁷⁾ and potentially, those ROMs can enable real-time simulation of aeroelastic systems⁽⁸⁾ by including nonlinear solutions.

Different strategies may be employed for the simplification of a GFEM into a nonlinear structural ROM. Two possible approaches are the equivalent beam modeling and the nonlinear modal solution. The variational asymptotic beam sectional analysis (VABS)⁽⁹⁾ is one of the methods for equivalent beam modeling of anisotropic composite high-aspect-ratio structures. It uses information on the section’s material, geometry and layup to calculate beam properties. This method has been developed for decades and often reduces in two or three orders of magnitude the computational costs, compared to finite element solutions⁽¹⁰⁾. VABS is able to quickly calculate displacements and recover 3D stress distributions. Refinements on VABS have been introduced to include modeling of initially twisted and curved beams⁽¹¹⁾ and distributed loads⁽¹²⁾ and more recently for initially curved and twisted smart beams⁽¹³⁾ and oblique cross-sections⁽¹⁴⁾. Since VABS is based on cross-section analysis, novel designs with intricate geometries may be harder to model, and methods to smear the stiffness and mass contributions of complex components have been suggested⁽¹⁰⁾, as well as novel methodologies with variational asymptotic approach that can be generalised to structures other than beams, such as in the Mechanics of Structure Genome solution⁽¹⁵⁾. However, all asymptotic solutions rely on approximations of small parameters and one of the requirements is the cross-section dimension to be much smaller than the deformation wavelength⁽¹⁰⁾. Usually, this assumption is not valid when there is local buckling between ribs, for example. Structure geometry and layup also affects the accuracy of VABS models. According to comparisons carried out with 3D GFEM solutions, the solution from VABS shows higher errors for wings with bend-twist coupling and sweep⁽¹⁶⁾. In such a scenario, methodologies for model reduction without the need of cross-section analysis or asymptotic behaviour may be useful, and these solutions are the focus of this paper.

Equivalent beam modeling may be performed using GFEM results at particular reference points along the wing. Different approaches consider particular beam models or beam cross-sections to match the GFEM results. Elsayed et al.⁽¹⁷⁾ proposed a methodology to fit properties for an Euler-Bernoulli beam and compared the results against other techniques such as empirical calculations. Bindolino et al.⁽¹⁸⁾ demonstrated static solutions approximating a GFEM with a wing of rectangular cross-section containing stiffeners. An optimisation is employed in this case to find suitable properties for the wing cross-section along the reference axis. More recently, Stodieck et al.⁽¹⁹⁾ proposed a method that calculates an isotropic Timoshenko beam based on 13 properties extracted from GFEM linear results. This last method is able to capture the modal frequencies and eigenvectors with high accuracy under linear assumptions. The equivalent beam can then be solved with a geometrically nonlinear solution like the intrinsic beam theory, displacement-based or strain-based nonlinear FE analysis or even finite volume (FV). This method was selected for inclusion in this paper, due to its automated procedure and potential for application in optimisation.

Another approach for nonlinear structural ROMs is to increment the classical modal solutions with nonlinear terms accounting for stiffness and displacement changes in nonlinear conditions of high displacements^{(6, 20)}. There are different methods following this nonlinear modal approach. An excellent review of those methods is provided by Mignolet et al.⁽²¹⁾. In general, nonlinear static FE solutions are used for identification. The nonlinear modal solutions were initially applied to problems of shell structures constrained at all sides, where the nonlinear displacements are in the order of magnitude of the thickness⁽²¹⁾. Applying this method to cantilevered structures was challenging, due to the identification of many nonlinear terms, but good results were achieved for both beams and wing boxes^(22–25).

This work compares different ROMs extracted from a wing GFEM at high displacements conditions. Extraction of ROMs for linear analysis is a common practice, but it is important to evaluate the accuracy of equivalent models at nonlinear conditions. It is known that phenomena like local buckling may affect the general structural behaviour in bending and torsion, and these phenomena may be predicted even for beam solutions⁽²⁶⁾, but the behaviour of a beam may differ from the built-up structure due to nonlinearities such as the Brazier effect⁽²⁷⁾. Usually, high-aspect-ratio wings are not designed to have instabilities, but it is important to check if the different ROM strategies that have been employed are yielding accurate predictions for high displacements. The main contribution is to explore the suitability of different ROMs when modeling large wing displacements, up to conditions where the cross-section deforms significantly.

Three different GFEM ROMs are explored: an enhanced implicit condensation and expansion (EnICE) model^{(6, 20)}, a finite element beam model, and a finite volume beam model calculated from the same geometrically nonlinear built-up wing structure. An FV beam is used in addition to the FE one in order to rule out the possibility of shear locking affecting the solution and to provide an additional check. These three models were chosen due their previously encouraging performance and also their availability for benchmark comparisons in the research groups involved in this work.

The University of Bristol Ultra-Green (BUG) wing GFEM is used as a reference model which was derived from the truss-braced high-aspect-ratio wing of the SUGAR Volt aircraft configuration⁽²⁸⁾. Bend-twist coupling was introduced in the wing structure by defining anisotropic material properties in the skins. This model is representative of a realistic high-aspect-ratio wing.

Following a brief description of the formulations of the different ROMs used for the comparisons, the reference model is detailed and the accuracy of the equivalent beam obtained checked in terms of frequencies and mode shapes. Then, the nonlinear behaviour of the BUG wing is investigated when large displacements are achieved. Finally, the static and dynamic tip responses are presented and the accuracy and usage of the different ROMs is discussed.

2.0 STRUCTURAL ROM FORMULATIONS

Three different models were used to compare performance and accuracy in this paper: an FE beam, an FV beam and a modal structural ROM using the EnICE method. This section describes briefly the estimation of the equivalent beam model from the GFEM and the formulations of the FV and the EnICE methods.

2.1 Equivalent beam modeling

As discussed in the introduction, multiple approaches may be followed to compute a representative beam for a built-up wing. Depending on the kind of beam solution sought, different parameters may be calculated. In this work, an automated method that calculates an equivalent Timoshenko beam from the $6 \times 6$ flexibility matrix was used⁽¹⁹⁾.

The extraction of the flexibility matrix follows the Beam Property Extraction procedure outlined by Malcolm and Laird⁽²⁹⁾. A reference line approximately aligned to the wing is proposed and reference points are selected close to the rib locations. The reference points are linked to the nodes of associated ribs using rigid elements that transform the motion of a set of pre-defined nodes into an average motion of the reference point. These rigid elements avoid local stiffening. The reference points also define the end nodes of the beam elements. Then, a series of linear FE solutions are calculated, applying forces and moments in three directions, individually over each element. Results of displacements and rotations are used to construct the flexibility matrix C.

Once the flexibility matrix is determined, the following Timoshenko beam properties are calculated analytically from its entries:

• section offsets of the shear center with respect to the element end nodes

• torsion and bending stiffness constants $\overline{J}$ , ${I_1}$ , ${I_2}$ and ${I_{12}}$

• section centroid offsets from the elastic axis

• section area

• Timoshenko shear coefficients

The offsets of the shear center at the ends of each element are calculated from rotation angles needed to cancel the bend-twist coupling terms of the flexibility matrix, namely ${C_{45}}$ and ${C_{46}}$ . Bending and torsional stiffness are then calculated along the elastic axis. The neutral axis is calculated from the beam kinematics, using the already calculated shear axis. From the results of displacements due to axial loads, section areas are calculated. Finally, shear coefficients are obtained from the rotation results, subtracting the bending contributions from them.

For details about the properties calculations, refer to Stodieck et al.⁽¹⁹⁾. Using the properties of a Timoshenko beam, the parameters of the beam elements are defined for the FE analysis. An advantage is that there is no need to define a particular cross-section for the beam definition. Based on the beam properties calculated, any FE code can be used to define the Timoshenko beam. For this work, the NASTRAN SOL 400⁽³⁰⁾ sequence (version 2014) was selected, which is an implicit nonlinear solver, but any other commercial code could have been used.

2.2 Finite volume beam

The FV beam was included in the comparisons to provide an additional check for the beam solution and to avoid the problem of shear locking that is common in slender structures modeled with first-order elements. The formulation of the finite volume beam was introduced by Ghiringhelli et al.⁽³¹⁾ and recently applied to nonlinear aeroelastic analysis^(1,32) . It considers the equilibrium equation, i.e.

(1) \begin{align}{\vartheta ^{\prime}} - {R^T}\vartheta + f = 0,\end{align}

where $\vartheta $ is the column matrix of internal loads, R represents the moment arm of the internal loads when analysed in an infinitesimal element and f is the column matrix of external loads. The derivative of $\vartheta $ is performed with respect to the scalar quantity $\xi $ that represents the arc-length along the beam.

For the FV formulation, the Equation (1) is written in a weak form, integrated over the length of each element. Different from variational formulations, no energy quantities are considered. The integral of the equilibrium equation is given by

(2) \begin{align}\int_a^b \left( {{\vartheta ^{\prime}} - {R^T}\vartheta + f} \right)d\xi = 0,\end{align}

with a and b representing the end-points of the segment where the integration is performed. After an integration of Equation (2) by parts, a relation is obtained between the internal loads $\vartheta $ at the ends a and b and the integrated external loads along the beam. The complete formulation depends on the constitutive relations that are needed to relate internal loads and generalised strains, as described in Ghiringhelli et al.⁽³¹⁾.

The integral of Equation (2) is performed for every element of the beam. Each element has three nodes and two evaluation points, as illustrated in Fig. 1. There may be offsets between the each node and the reference line for the beam. The displacements and rotations are interpolated as quadratic functions from the reference values at the three nodes. The evaluation points are chosen so that an exact solution is obtained for a tip load case.

Figure 1. Element of the FV beam formulation.

The framework with the FV formulation is suitable for multibody dynamics simulations, incorporating large rotations by means of modified Gibbs-Rodriguez parameters. Constitutive relations are considered in the material reference frame.

Different from an FE formulation, the FV approach leads directly to equilibrium conditions at collocated points, instead of energy integrals over the elements. Fundamentally, the disadvantage is a loss of symmetry of the linear and linearised beam matrices, but the lack of energy relations involving shear will automatically remove the problem of shear locking. For more details, the reader is referred to Ghiringhelli et al.⁽³¹⁾.

A set of nonlinear differential algebraic equations are obtained for the FV beam in which the equilibrium equation is augmented with the inertial load and the enforced motions are imposed through the use of holonomic constraints. An implicit scheme is finally used for time-integration.

2.3 Enhanced implicit condensation and expansion

The enhanced implicit condensation and expansion approach was developed initially to incorporate nonlinear stiffness into a modal solution, identifying additional terms with nonlinear static FE solutions⁽³³⁾. Later, the approach was expanded by including nonlinear displacements⁽³⁴⁾. Recently, the nonlinear displacements were used to correct the inertia forces in large displacements conditions^(6,20) . The method used in this work is the same of the cited references, where it was presented as the implicit condensation and expansion method. The name of the approach has been updated with the word enhanced to differentiate it from the method introduced by Hollkamp and Gordon⁽³⁴⁾. With the enhancement, the kinetic energy related to the nonlinear displacement is included in the dynamic equations, and this effect is visible for tip out-of-plane displacements around 30% of the semi-span, depending on the loading applied⁽³⁵⁾.

For the nonlinear modal approach, an additional basis of shape functions enhance the set of linear modes to express the nonlinear displacements, like the shortening in the span-wise direction. Using the FE notation and discretisation, the column matrix of displacements and rotations is expressed in terms of modal amplitudes, i.e.

(3) \begin{align}u = \Phi q + \Psi r\left( q \right),\end{align}

where u is the column of displacements/rotations, $\Phi $ is the matrix of n different linear modes along its columns and $\Psi $ is the matrix of s different dual modes, which is the additional basis of shape functions and not properly a set of modes. The amplitudes of the linear modes are in the column matrix q and the amplitudes of the dual modes are in the column matrix r. One of the key assumptions of this method is that the amplitudes of the dual modes r are functions of the amplitudes of the linear modes q.

Using the above representation of displacements, it is possible to express the kinetic energy as a function of the rates $\dot q$ and the mass matrix M of the structure, as well as the derivatives of the dual modes amplitudes, such as

(4) \begin{align}T = \frac{1}{2}{\dot q^T}\underbrace {\left( {{I_n} + 2{M_X}\frac{{\partial r}}{{\partial q}} + {{\frac{{\partial r}}{{\partial q}}}^T}{M_\Psi }\frac{{\partial r}}{{\partial q}}} \right)}_{M'}\dot q,\end{align}

where ${I_n}$ is an identity matrix of size $n \times n$ and the terms ${M_X}$ and ${M_\Psi }$ are defined as

(5) \begin{align}{M_X} = {\Phi ^T}M\Psi ,\end{align}

(6) \begin{align}{M_\Psi } = {\Psi ^T}M\Psi .\end{align}

The derivative of modal amplitudes $\frac{{\partial r\left( q \right)}}{{\partial q}}$ used in Equation (4) is a Jacobian matrix where each column corresponds to a derivative of r with respect to a different component of q. It is important to notice that the contribution of nonlinear displacements is equivalent to a modified mass matrix $M'$ in the expression for the kinetic energy.

From the kinetic energy, applying Euler-Lagrange equations is a straightforward path to describe the dynamics of the system as

(7) \begin{align}\frac{d}{{dt}}\left( {\frac{{\partial T} }{ {\partial {{\dot{q}}_i}}}} \right) - \frac{{\partial T} }{ {\partial {q_i}}} +\frac{{\partial U} }{ {\partial {q_i}}} = {f^T}\frac{{\partial u}}{ {\partial {q_i}}}\,\,\,\,\,\,for\,\,\,\,\,i = 1, \ldots ,n,\end{align}

where the the potential energy U is the elastic strain energy of the system and f is the column matrix of applied loads, using the FE notation, where these loads may be forces or moments, according to the degree-of-freedom involved. The number of nonlinear equations corresponds to the number of linear modes included in the displacements description.

The second key idea of the EnICE method is that the elastic modal forces $\frac{{\partial U}}{{\partial {q_i}}}$ in Equation (7) are a function of the linear modal amplitudes q. In order to simplify the fitting, the elastic forces are split into linear and nonlinear components, with the linear term calculated using the linear modal stiffness matrix $\Lambda $ . The nonlinear modal force ${f_{NL}}(q)$ is calculated from the reference solutions isolating the nonlinear residue from the elastic force, i.e.

(8) \begin{align}{f_{NL}}\left( q \right) = \frac{{\partial U}}{{\partial q}} - \Lambda q,\end{align}

where $\Lambda $ is the linear modal stiffness matrix. In the reference cases, the elastic force $\frac{{\partial U}}{{\partial q}}$ is known a priori, because it corresponds to the external applied load, and the amplitudes of the linear modes q are calculated by projecting the obtained displacements onto the linear modes of the structure.

In the end, two fittings are used to represent the system: a fitting for the amplitudes of dual modes $r(q)$ and a fitting of nonlinear elastic modal forces ${f_{NL}}(q)$ .

The final equations of motion are developed from the Euler-Lagrange set (Equation (7)) as

(9) $$\eqalign{ & \left( {{I_n} + {M_X}{{\partial r} \over {\partial q}} + {{{{\partial r} \over {\partial q}}}^T}M_X^T + {{{{\partial r} \over {\partial q}}}^T}{M_\Psi }{{\partial r} \over {\partial q}}} \right)\ddot q \cr & \,\,\,\,\,\,\,\, + [{M_X}{d \over {dt}}{{\partial r} \over {\partial q}} + {d \over {dt}}{{{\partial r} \over {\partial q}}^T}M_X^T + {d \over {dt}}{{{\partial r} \over {\partial q}}^T}{M_\Psi }{{\partial r} \over {\partial q}} + {{{\partial r} \over {\partial q}}^T}{M_\Psi }{d \over {dt}}{{\partial r} \over {\partial q}} \cr & \,\,\,\,\,\,\,\, - \left[ {{{\dot q}^T}\left( {{M_X} + {{{{\partial r} \over {\partial q}}}^T}{M_\Psi }} \right){{\left( {{{{\partial ^2}r} \over {\partial {q^2}}}} \right)}_{s \times {n^2}}}} \right]_{n \times n}^T]\dot q \cr & = {\left( {\Phi + \Psi {{\partial r} \over {\partial q}}} \right)^T}F - \Lambda q - {f_{NL}}\left( q \right), \cr} $$

where the sub-indexes indicate a particular shape for the matrices. For example, the second derivatives of the s dual modes with respect to the n linear modes in Equation (9) is expressed as

(10) \begin{align}{\left( {\frac{{{\partial ^2}r}}{{\partial {q^2}}}} \right)_{s \times {n^2}}} = \left[ {\frac{{{\partial ^2}r}}{{\partial {q_1}\partial q}}\quad \frac{{{\partial ^2}r}}{{\partial {q_2}\partial q}} \ldots \frac{{{\partial ^2}r}}{{\partial {q_n}\partial q}}} \right],\end{align}

where the n different derivatives of the $s \times n$ Jacobian $\frac{{\partial r}}{{\partial q}}$ are stacked side by side. Reshaping is also used in Equation (9) to represent the time derivative of the Jacobian, i.e.

(11) \begin{align}\frac{d}{{dt}}\frac{{\partial r}}{{\partial q}} = {\left[ {{{\left( {\frac{{{\partial ^2}r}}{{\partial {q^2}}}} \right)\!\!}_{sn \times n}}\dot{q}} \right]_{s \times n}}{\kern 1pt} \;\textrm{where}\;\;\;\;{\left( {\frac{{{\partial ^2}r}}{{\partial {q^2}}}} \right)_{\!\!sn \times n}} = \left[ \begin{array}{l}\dfrac{{{\partial ^2}r}}{{\partial {q_1}\partial q}}\\ \\[-7pt] \dfrac{{{\partial ^2}r}}{{\partial {q_2}\partial q}}\\ \\[-7pt] \vdots \\\\[-7pt] \dfrac{{{\partial ^2}r}}{{\partial {q_n}\partial q}}\end{array} \right].\end{align}

Equation (9) defines a set of n second-order ordinary differential equations, and the system can be integrated in time to solve for the amplitudes of the linear modes, allowing to recover the magnitudes of dual modes and the complete nonlinear displacements.

3.0 MODEL DESCRIPTION

3.1 Wing box model

The structure used for all comparisons is a wing box built out of graphite/epoxy composites modeled in MSC NASTRAN with 48 750 degrees-of-freedom, using CQUAD4 (quadrilateral), CTRIA3 (triangular) and CBEAM (beam) elements, besides concentrated masses and rigid connectors. Due to large displacements conditions, nonlinear properties were introduced for the shell elements using the PSHLN1 card, a piece of the input defining properties like integration scheme and classification of the shell elements.

The wing box corresponds to a high-winged design, the University of Bristol Ultra-Green (BUG) wing, which is derived from the truss-braced high-aspect-ratio wing of the SUGAR Volt aircraft⁽²⁸⁾. This wing has negative dihedral, and the fiber orientation is designed to allow bending-torsion coupling. The semi-span of the wing box is $25.9$ m. A general view of the model is presented in Fig. 2.

Figure 2. BUG wing model, with detail for wing ribs and engine in the bottom view.

The boundary conditions are compatible with a high-wing design. The root nodes of the wing are constrained with symmetry conditions: no span-wise translations (y-direction) and no rotations about the x and z axes. Besides that, the nodes on the section connected to the fuselage are constrained in translation along the x and z directions.

The quality of reduced-order models for small displacements depends on how well the first linear modes are reproduced. With this set of constraints, the first free-vibration mode has a frequency of $1.30$ Hz. The first four modes are shown in Fig. 3, but engine displacement is not displayed along with the wing box. The first mode is a bending one, with a relatively low frequency, characteristic of high-aspect-ratio wings. The second mode is an in-plane mode, with displacements along the x axis, while the third one shows some torsion towards the wing tip. The fourth mode is is a bending one, similar to the first mode, but has a higher curvature close to the tip.

Figure 3. First three free-vibration modes of the BUG wing box.

The mass of the wing is mostly lumped. Points connected to neighbouring nodes on the structure via rigid connections of the kind RBE3 (Fig. 4) are used to transfer inertia forces. The rigid connections RBE3 do not constrain the relative motion of the nodes on the wing, avoiding unrealistic stiffening. For the RBE3 connections, the motion of a reference point is the uniformly weighted average of the motion of chosen neighbouring nodes. These same reference points with lumped masses are also used to apply external loads and obtain a reference solution for the dynamic load case investigated later.

Figure 4. Point mass connected to neighbouring nodes in a section of the wing.

3.2 Beam models

In order to evaluate the beam model that reproduces the wing box, it is necessary to first compare frequencies and mode shapes of the linear free-vibration problem. This discussion is valid for both the FE and the FV beam models. A reference line was chosen to represent the wing box, composed of 30 distributed nodes. The reference line and points are shown in Fig. 5. It was chosen to pass approximately through the center of the wing box. The method used to reduce the box to a beam allows for disparities between the elastic axis and the reference line. This approach allows for a better representation of the element stiffness properties. The elastic axes of the individual elements are represented in red in Fig. 5. The mass related to the set of elements close to each reference point is concentrated at a single point representing their center of gravity, which may have an offset with respect to the corresponding reference point. This offset is especially significant at the node closest to the engine mass, as shown in Fig. 6. Rotational inertias are accounted for using the parallel axis theorem.

Figure 5. Reference line and points selected to represent the wing box structure.

Figure 6. Point masses and its offsets relative to the reference points.

From the beam reduction process, an equivalent beam model was generated with match in frequencies with the complete wing box model within 1% up to the seventh linear mode. Table 1 compares the frequencies of the two models. In this case, the errors are defined relative to the wing box frequencies, which are the reference values. In general, the matching is better for the lowest modes, up to mode 7. For the higher modes, some discrepancy happens due to local effects, but the low-frequency response will still be dominated by the first modes.

Table 1 Comparison of modal frequencies for the wing box and the equivalent beam model

In terms of mode shapes, the Modal Assurance Criterion (MAC) was used to compare displacements for the points of the equivalent beam and the corresponding points connected to the wing box’s neighbouring nodes via RBE3 connections. Table 2 compares the beam modes and the wing box modes. The MAC is an approximate comparison between two mode shapes. It disregards the mass in the process, but it will be $1.0$ if two mode shapes are identical. For this table, the degrees-of-freedom considered were the three translations and the rotation along the y-direction. A comparison was made between the mode shapes from the beam model and the mode shapes from the GFEM, taking results from the reference nodes only. The formula of the MAC is

(12) \begin{align}MAC = \frac{{{{\left( {\Phi _1^T{\Phi _2}} \right)}^2}}}{{\left( {\Phi _1^T{\Phi _1}} \right)\left( {\Phi _2^T{\Phi _2}} \right)}},\end{align}

where ${\Phi _1}$ and ${\Phi _2}$ are the two modal shapes from the different models compared, the GFEM and the equivalent beam. For the example shown in Table 2, all the MAC values in the diagonal are 1 until mode 7, while there are discrepancies for mode 8 onwards, just like in the frequencies comparison. Even though the diagonal terms in the MAC matrix are not all ones, it can be said that the agreement in low-frequency response is good. Obtaining a beam with that level of agreement is only possible because the offset between the shear center and reference axis is included in the equivalent beam model.

Table 2 Modal Assurance Criterion comparison of modal shapes for the equivalent beam along rows and the wing box along columns

It is interesting to note that mode 8 on the GFEM corresponds to a local box cross section deformation mode. That occurs at the location where the engine inertia forces are being transferred into the wing through four discrete nodes on the wing section. Clearly, this local mode cannot be captured by the beam model. The mode mismatch also appears in the results of Table 2, where the MAC between the beam mode 8 and the GFEM mode 9 is equal to 1, indicating that the GFEM mode 8 was skipped. The GFEM in Stodieck et al.⁽¹⁹⁾ was modified compared to the model used in this study, by connecting the discrete engine mass to a larger number of nodes on the wing cross section. Using this updated GFEM improved results were achieved and it was found that the beam and GFEM normal modes up to 20Hz were closely correlated (2% maximum frequency error).

3.3 EnICE model

The EnICE model of the BUG wing is based on 18 degrees-of-freedom, corresponding to the amplitudes of the first 18 linear modes. The nonlinear displacements and elastic forces are fitted with artificial neural networks composed of four neurons in the hidden layer. Feed-forward neural networks with a single hidden layer are used. The fitting error is minimised by adjusting the weights and bias of the neural network. In the hidden layer, a tangent sigmoid activation function is employed, and the training follows a Bayesian approach. Since the neural networks are mathematical functions with many local minimum points, the optimisation is performed a few times, using different random initial conditions. Each trained network is then evaluated on a separate subset of reference data that was not used for training and the best one is selected.

For the training of the EnICE model, a total of 21 560 nonlinear static solutions were obtained. Loads in the shape of linear modes excited large amplitudes for extraction of nonlinear displacements and nonlinear forces. A separate code for calculation of reference results was used. This routine called the finite element solutions successively, storing results of displacements and rotations for all the degrees-of-freedom of the model. The input load for each case is given by

(13) \begin{align}F = K\Phi {q^*},\end{align}

where ${q^*}$ is the column matrix with of amplitudes of linear modes that would be obtained if the response was linear. These amplitudes define the magnitude of the loads to be used for the reference solutions. Initially, the ranges of these magnitudes are explored for each one of the components. These loads should be high enough to capture the displacements of interest. However, excessively high loads may lead to convergence problems, and exploratory solutions are needed to establish a range of values for ${q^*}$ . The particular values of $q^{*}$ selected for the reference solutions are random, inside the range determined from exploratory analyses.

The nonlinear displacements were represented as combinations of 15 dual modes.

4.0 COMPARISON OF RESULTS

Different load cases were used to analyse the accuracy of the reduced-order models compared. The benchmark comparisons were performed to (1) check each ROM with different load excitations and conditions, and (2) compare the accuracy of each ROM against the other ones and the full model whenever possible.

This section compares dynamic and static results obtained for different load cases. Tip displacements were used to monitor the differences between the investigated models.

One aspect of practical utility regarding GFEM solution in large displacements conditions is that the presence of rigid connections may lead to convergence problems. Taking the BUG wing as an example and using the MSC NASTRAN SOL 400, dynamic solutions cannot be obtained easily for tip displacements beyond 16% of the semi-span unless sufficiently small time steps are adopted. If the rigid connections are removed, however, then larger displacements can be calculated.

4.1 Load cases for the benchmark comparisons

The load cases were selected to explore different input profiles, as well as the effects of inclusion of follower forces. Solutions were obtained for the ROMs and the full model. For higher loads, the dynamic nonlinear solution with the wing FE model was not obtained anymore, despite the numerous attempts to achieve convergence with different parameters sets. In order to explore larger displacements, only static loads were used.

While Case 1 is the application of a step load input at two different points along the wing, Case 2 explores a sinusoidal load with frequency of $0.125$ Hz. For both dynamic cases, the damping is a stiffness-proportional Rayleigh damping, and the values listed are the constants of proportionality for the stiffness term of the Rayleigh damping. Case 3 is a comparison of static results considering larger deflections. Table 3 summarises the three load cases compared.

Table 3 Table of load cases for benchmark comparisons

4.2 Responses with different methods

This subsection details the results obtained from the load cases. Time-domain tip responses were selected to compare the behaviour of the different ROMs. The out-of-plane displacement is particularly important to indicate the degree of nonlinearity. Tip displacements around 15% of the semi-span will generally be associated to a mild nonlinearity.

4.2.1 Cases 1 and 2: Dynamic analyses

The first two cases represent dynamic loads. The amplitudes were chosen with smaller values to make the comparison with GFEM results possible. For Case 1, the out-of-plane tip displacement is shown in Fig. 7. Since the displacements are kept below 15% for this analysis, the nonlinearity is not strong, and the differences between the full solution and the linear one are relatively small. A zoom (Fig. 8) reveals the differences, showing that the linear solution has the highest error among the ROMs, and the beam-based approaches are clustered, but preserving relatively small errors. The EnICE results achieved a good accuracy for this level of displacements, because the linear behaviour matches the first 12 modes and the light nonlinearity is well-captured by the fitted functions. The span-wise displacements are compared in Fig. 9. As expected, the linear solution does not present an accurate span-wise displacement, because the shortening along this direction is a typical nonlinear effect. The wing still has a small linear response due to its negative dihedral, but the nonlinearity dominates this direction, and all nonlinear ROM responses were able to capture it.

Figure 7. Out-of-plane displacements for Case 1.

Figure 8. Zoom at the out-of-plane displacements.

Figure 9. Tip spanwise displacements for Case 1.

In Case 2, the excitation has the same amplitude as Case 1, but it is sinusoidal with a low frequency: $0.125$ Hz. This case turned out to have small amplitude response, as shown in Fig. 10. In a long simulation of 100s, all the solutions remained accurate, and there was no phase error accumulated. The only differences are visible with a zoom into the first peak, as shown in Fig. 11, where it can be noticed a small error of the linear solution, while the nonlinear solutions remained closer to the GFEM reference. For the subsequent peaks, the differences are even smaller. This difference happens because the transient response is more affected by the presence of nonlinear displacement components in the solution. For harmonic excitation, it was already expected that no phase error accumulation would happen, because the response is driven by the applied load, in the long runs.

Figure 10. Out-of-plane displacements for Case 2.

Figure 11. Zoom at the out-of-plane displacements.

4.2.2 Case 3: Static analyses

As explained in the beginning of this section, dynamic results could not be obtained for large displacements, due to lack of convergence with the rigid elements included in the model. However, it is possible to delete the rigid elements together with the concentrated masses and obtain static results up to very large displacements.

For this comparison, a tip vertical load was applied, and the results for tip displacements predicted by the the different ROM methods were compared to the GFEM reference up to a tip vertical force of 70kN. Since this loading corresponds to a tip force, a high stress is generated at the outboard region of the wing. That contributes to local buckling and deformation of the cross section.

Figures 12 and 13 show the out-of-plane and span-wise tip displacements as the load is increased. In the span-wise direction, the linear model shows an error relative to the nonlinear ones, even for small loads. Until 60% of the maximum load, the agreement between all models is really good in the vertical direction. When the tip vertical displacement achieves values higher than 15% of the semi-span, however, a significant difference is observed between the GFEM solution and the beam-based models, while the EnICE predictions remain close to the reference. This can be observed in the detail of out-of-plane displacements shown in the Fig. 14.

Figure 12. Out-of-plane displacements for static loading.

Figure 13. Span-wise displacements for static loading.

Figure 14. Out-of-plane displacements for static loading, with detail for large displacements.

Finally, another interesting aspect of the ROM solution is the possibility to capture local effects with an enrichment of the basis of dual modes considered in the nonlinear displacements. These are purely nonlinear features accessible to modal ROM predictions. In order to analyse this capability, the cross section of the minimum principal strain was selected for this study. When a tip load is applied, this cross section is the one indicated in Fig. 15.

Figure 15. Position of cross section selected for analysis.

In order to capture the local effects of the cross section deformation, it was necessary to introduce additional dual modes specific for the nodes on the selected cross section. With an additional set of 15 dual modes, a result was obtained for the cross-section shape that indicates both the local buckling of the upper surface and the crushing of the cross section due to the Brazier effect. The comparison of cross-section geometries for the reference and the ROM solutions is presented in Fig. 16, along with the undeformed cross section. Even though the accuracy of the EnICE solution may still be increased, it indicates that local effects can be predicted with ROM models enriched with nonlinear features on the spot of interest.

Figure 16. Cross-section deformation after tip load.

4.3 Discussion of results

From the tip responses of dynamic load cases, a good agreement was observed between the solutions of three different nonlinear ROMs considered. The step response pointed out to accumulation of phase errors between the beam solutions and the modal approach. However, for the harmonic excitation, excellent matching was observed. Due to the large deformations, the linearised GFEM failed to accurately capture the correct response for span-wise components.

When higher tip deflections are achieved, the beam solutions are not able to represent phenomena like local buckling and Brazier effect, as expected. Therefore, softening due to the deformation of the cross section is not properly modeled using these approaches. The load cases investigated tip loads, which are associated to a high curvature in the outboard region of the wing. Since this region is not designed for such high loads, local buckling is anticipated with tip loads. Figure 17 shows a qualitative comparison of minimum principal strain distributions when the load applied is a tip one and when it is a distributed load in the shape of the first bending mode. In the first case, the loads are concentrated at the tip, and there is buckling. In the second case, the tip displacement is similar, but the strain is more uniformly distributed and there is no buckling.

Figure 17. Comparison of minimum principal strain distributions for two different loading conditions.

The EnICE is able to capture effects of cross section deformation, which explains the matching between tip displacements obtained with the EnICE and the GFEM solutions, for static solutions. In terms of torsion, the EnICE prediction is also matching the reference, if the translational degrees-of-freedom are observed. Figure 18 shows the geometry of the tip section of the deformed wing under the tip vertical load of 70kN, comparing the reference solution with the EnICE prediction. This view is interesting to show the ROM errors relative to the cross-section dimensions and to indicate the torsion discrepancies. The torsion angle is $ - 21.8$ degrees for the GFEM reference and $ - 22.4$ degrees for the ROM prediction. Although no torsion comparisons were presented previously, this result gives an indication that the modal ROM is able to capture torsion correctly, if the translational degrees-of-freedom are used to find the cross-section shape.

Figure 18. Tip cross section of deformed wing.

It was not possible to obtain the GFEM reference for dynamic cases in large displacements, but from the results of lower amplitudes, it is expected that the ROM predictions would be accurate for large amplitude dynamic solutions. In this case, the ROM in dynamic conditions has the versatility of predicting a response when the GFEM solution shows convergence difficulties.

In terms of computational efficiency and model preparation, it can be said that the modal solution represented by the EnICE approach is less automated in the off-line phase, requiring considerations of reasonable training loads for the identification of nonlinear stiffness and displacements. The calculation of the high-fidelity solutions took approximately 37 hours from a Xeon E5 processor. However, the reference solutions may be calculated in parallel. After getting the reference solutions, the training of the neural networks involves additional time, depending on the number of degrees-of-freedom considered and the number of trainings required to achieve the desired error levels. With the current level of automation, the preparation of an EnICE model takes a few days. The extraction of equivalent beams, however, is based on linear solutions for the extraction of the flexibility matrix. These reference data are much simpler to obtain, compared to the nonlinear static solutions used for the EnICE model. The process of generating an equivalent beam model is thus considerably faster than building a nonlinear modal ROM. The number of degrees-of-freedom considered in the ROMs is also different. The EnICE approach deals with fewer degrees-of-freedom due to its modal nature, while the beam solutions will usually include more than 100 degrees-of-freedom. Even though the modal solution has less dependent variables involved, the computational efficiency is not as high as the classical modal approach due to the complexity of the nonlinear functions and its derivatives, which are neural networks, for now. Table 4 presents an approximate comparison of computational times required from each method to simulate 1s of dynamics.

Table 4 Computational time comparisons for 1s of physical time simulation

5.0 FINAL REMARKS

This work compared the accuracy of three different ROMs representative of a built-up wing structure of a high-aspect-ratio aircraft. Two of the models were beam representations while the third one was a nonlinear modal method. The objective was to check the different models for different load cases related to high-amplitude dynamic responses and static simulations. The ROM results were compared to a nonlinear NASTRAN SOL 400 reference solution.

In general, all three ROMs showed good accuracy up to tip vertical displacements around 15% of the wing semi-span. The dynamic solutions had a good agreement. Only the step response showed a mild accumulation of phase error between the different ROMs over time. As the displacement amplitudes were increased, however, nonlinear deformations resulted in a disparity between the modal ROM and the beam-based ones, as noticed in the static analysis.

An interesting point observed is that the reference NASTRAN SOL 400 did not converge for the high-amplitude dynamic responses (> 16% of the wing semi-span), for the conditions simulated. However, the ROM solutions were obtained easily for such large displacements. This observation indicates that ROMs may be a more practical approach for cases where the convergence is harder. Further investigations are needed to verify this issue with the FEM solution and to check the accuracy of ROM predictions in such conditions.

Each ROM may be applicable to different conditions. The equivalent beam models are cheaper to obtain, but those are applicable when the cross-section deformation does not affect the stiffness significantly. For conventional aircraft, that condition is usually met for all structural members. The nonlinear modal ROMs represented by the EnICE method in this work should be useful for non-conventional designs with higher cross-section deformation or for development of control techniques when the number of structural degrees-of-freedom is limited but the nonlinear displacements cannot be neglected.

Both the beam solutions and the modal approach revealed the capability of greatly reducing the computational costs of GFEMs while retaining a good accuracy in the nonlinear range of displacements. It is expected that further research on ROM methodologies will decrease even more the simulation costs and ultimately allow real-time simulation frameworks to use the structural ROMs.