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^{(Reference Stodieck, Cooper, Neild, Lowenberg and Iorga19)}.

The extraction of the flexibility matrix follows the Beam Property Extraction procedure outlined by Malcolm and Laird^{(Reference Malcolm and Laird29)}. 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.^{(Reference Stodieck, Cooper, Neild, Lowenberg and Iorga19)}. 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.^{(Reference Ghiringhelli, Masarati and Mantegazza31)} and recently applied to nonlinear aeroelastic analysis^{(Reference Calderon, Cooper, Lowenberg, Neild and Coetzee1,Reference Calderon, Cooper, Lowenberg, Neild and Coetzee32)} . 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.^{(Reference Ghiringhelli, Masarati and Mantegazza31)}.

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.^{(Reference Ghiringhelli, Masarati and Mantegazza31)}.

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^{(Reference McEwan, Wright, Cooper and Leung33)}. Later, the approach was expanded by including nonlinear displacements^{(Reference Hollkamp and Gordon34)}. Recently, the nonlinear displacements were used to correct the inertia forces in large displacements conditions^{(Reference Medeiros, Cesnik and Coetzee6,Reference Medeiros, Cesnik and Coetzee20)} . 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^{(Reference Hollkamp and Gordon34)}. 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^{(Reference Medeiros35)}.

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.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^{(Reference Bradley, Allen and Droney28)}. 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.^{(Reference Stodieck, Cooper, Neild, Lowenberg and Iorga19)} 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.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