Aeroelastic method to investigate nonlinear elastic wing structures

Stiffness directions of wing structures are already part of the optimisation in aircraft design. Aircraft like the A350 XWB and the Boeing 787 mainly consist of such composite material, whose stiffness directions can be optimised. To proceed with this stiffness optimisation, the aim of this work is to modify and optimise also the linear stress-strain relation. On that account, the Hooke’s law is exchanged by a multi-linear formulation to analyse any nonlinear elastic structural technology on wing structures. The wing structures, which are used to investigate the nonlinear behaviour, are deduced from a mid-range and a long-range aircraft configuration. These wings are analysed with an extended beam method and coupled with a VLM solution to calculate the aeroelastical loading. The proposed beam method is capable of analysing any multi-linear wing structure technology. A degressive structural behaviour shows up a good potential to reduce the bending moment which is one of the main drivers of the structural weight.


Introduction
Besides the level flight loadcase, an aircraft carries out manoeuvres and encounters gusts. These two types of loadcases define the structural dimensions of the wing primarily.
To reduce structural weight, load alleviation in passive and active variations are prevalent in aircraft design for those loadcases.
On the passive side, carbon fibre reinforced plastics are the key drivers. Aeroelastic tailoring describes the optimisation of those reinforced plastics, set up in laminates, to get the optimal anisotropic stiffness behaviour. Additionally, well designed stacking sequences can activate the 1 3 bending-torsion coupling to influence the load distribution [1,2].
On the active side, several control surfaces can affect the lift distribution to create a beneficial loading (e.g. [3,4]). The advantage is that it can be activated and adjusted in the moment when the high loading occurs. Disadvantageous, the weight for the mechanical set up and the energy consumption increase.
According to statistical load data, an imbalance exists in the occurring between the design loadcases and the level flight loadcases [5]. The design loadcases are determining the structural weight, but occur seldom. This imbalance can be taken into account by active load alleviation (due to their on/off-states) but not by passive concepts. In this thought, it would be helpful, when the structure could act automatically and passively on its own and without further motorised mechanical system.
Thus, the next step of aeroelastic tailoring is to modify also the linear stiffness relation given by the Hooke's law ( = E ⋅ ). To counter these high loading cases passively, an optimal nonlinear elastic stiffness has to be designed. During the level flight the wing shall behave linearly and more rigid in order to stay close to the optimal flight shape. During manoeuvres and gusts the wing shall behave more flexible to use the bending-torsion coupling as load alleviation. The upward bending of an backward-swept wing, decreases the incident angle at the outer part of the wing which results in load shift towards the root. Exemplarily for a 2.5 g manoeuvre, when the loading and deformation increase, the wing becomes softer and bends up more compared to a linear elastic wing. After the manoeuvre the wing returns to its initial stiffness state. In doing so, the high load can be reduced while the level flight can be strengthened. This idea is comparable to a palm tree which bends downwards in stormy weather and stays straight upwards during sunny conditions. Natural fibres like sisal or coir show up such hyperelastic behaviours [6]. Some alloys (e.g. NiTi-alloy) have such a behaviour in a specific temperature range [7]. Elastomers are also part of this kind of materials. In the experiment of Brojan et al. such an elastomer is tested as a material for a beam structure [8]. Brojan et al. developed also a comprehensive method to analyse such stiffness relation in beams. Albeit, this method is connected with to much computational effort for an aeroelastic coupled analysis at this point.
One existing concept for nonlinear wing structures is the semi aeroelastic hinged wing which has been tested experimentally with the AlbatrossOne by Airbus [9]. The outer part of the wing is hinged nonlinearly to reduce the root bending moment for the design loadcases.
The goal of this paper is to propose a method which is capable to investigate any structural nonlinearities in the aeroelastic regime, and to reach a better understanding of nonlinear elastic stiffness for wing structures.
The following Sect. 2 introduces the proposed nonlinear beam method. The various models for the aeroelastic analysis are presented in Sect. 3. The validations against numerical and experimental cantilever beams are performed in Sect. 4. The capabilities of this approach, detailed insights into nonlinear elastic design, and a realistic application are provided in Sect. 5. The last Sect. 6 lists and discusses the key messages of the new method for nonlinear elastic stiffnesses.

Nonlinear elastic beam method
The following method is valid for a bending beam calculation with nonlinear elastic stiffnesses. More details of the bending theory are given in Gross et al. [10] In the classical beam theory the bending moment is equal to the integrated normal strains of the cutting surface: That is the static equilibrium which has to be fulfilled. As far as the normal stresses are linearly dependent on the z-coordinate (following the Hooke's law = E ⋅ and the kinematic relation = ⋅ z ) the equilibrium equation becomes: with the moment of inertia This Eq. 2 has to be inverted to get the curvature of the beam. The curvature is approximately equal to the second derivative of the bending line w" for small slopes.
With this Eq. 4, the differential equation of the deflection curve, the elastic line w can be derived. Now, for nonlinear elastic behaviour Hooke's law will be modified. The relation between the stresses and the strains , in Fig. 1, is defined stepwise with the following sets of endpoints. So, any multi-linear stiffness relation can be addressed.
With the elastic modulus E can be calculated for every linear step. Because E i is valid for i ⩽ < i−1 the initial stiffness E 0 can be defined as equal to E 1 . Therefore, the given elastic moduli can be written as: The stresses in the static equilibrium (Eq. 1) are then substituted to the following equation with the intermediate steps n. The number of steps n are defined by the maximum strain value max at the upper bound ( z max = b∕2 ) with ( n−1 < max ⩽ n ).
To solve the integral, the strains have to be substituted with the kinematic relations in Eq. 10. The left part is equal to the linear beam theory.
With integration by parts (Eq. 12), and expanding the binomial formulas with 2 and 3 (Eq. 13), the relation between bending moment and curvature can be derived. (Note: , and M b,y are variables while i , i and M b,y,i are constants.) The curvature endpoints, which are needed for the resubstiuation between Eqs. 12 and 13 can be defined from the strain endpoints and the half of the height of the beam b/2.
The static equilibrium in Eq. 9 is then solved with the rearranged formulation in Eq. 16, to highlight the nonlinear relation between the bending moment and the curvature.
which is valid for ( n−1 < ⩽ n ) with three stiffness constants A, B and C for each step.
This results in the following sets of stiffness constants, in which A 1 and C 1 are both equal to 0.
To obtain the unknown curvature , the inverse function of Eq. 16 is needed (compare Eq. 4). In this case three solutions can be found for each step. The following equations gives us the real part of the solution: The bending moment related endpoints for the sections of Eq. 23 can be calculated with Eq. 16.
The relation between curvature and bending moment (Eq. 23) is important for realising any desired nonlinear behaviour (see Sect. 5.2). The linear elastic solution is still included in the Eqs. 16 and 23 for n = 1.
Consequently, the bending line w can be calculated with two integration steps for a constant bending moment, because of w �� ≈ . To be capable of arbitrary bending moment distributions, the beam in Fig. 2 is analysed stepwise with the stress resultants M b,y , F z and M T . Shear correction and torsion are calculated linearly with the linear shear modulus G.
To create a different behaviour for upward and downward bending, the stiffness relation (Eq. 5 and 6 ) can be extended by sets of endpoints with negative values. To find the real part of the solution (Eq. 23), the absolute of the negative bending moment must be used. Also, the curvature endpoint set (Eq. 15) has to be used with absolute values. In the following, the positive curvature solution has to be considered as negative, to get the correct downward bending line w.
Otherwise, Eq. 16 is a complex solution of the inverse function of Eq. 23. Figure 3 shows the distribution of the strains and stresses for two different bending moments. The strains are assumed to be distribute linearly. Following this, the linear stresses lin act dependently also with a linear distribution. In contrast, the nonlinear stress distribution of nlin shows up the predefined knee when the maximum strain max is greater than the given strain endpoint 1 The z-coordinate of this knee is not fixed and moves towards the neutral axis in case of an increasing bending moment.

Simulation models
In order to analyse nonlinear elastic effects on wing structures, two backward swept wings are chosen: (a) in a typical size of a mid-range aircraft wing and (b) in a typical size of a long-range aircraft wing. And (c) a simplified rectangular swept wing is set-up for the validation, which is also based on the long-range aircraft properties (see Fig. 4). The midrange aircraft is based on the SE2A-MR configuration [11] and is comparable to an A320, while the long-range aircraft is of a similar size as the XRF1 [12] or an A350. The dimensions are given in Table 1.
The backward sweep ensures that the aeroelastic bending torsion coupling between the aerodynamics and the structures is affected. In addition, the influences of engines, fuel tanks, and structural kinks are omitted to show the pure effect of the nonlinear elastic behaviour.

Aerodynamic model
To create aerodynamic forces, the vortex lattice method (VLM) from the in-house tool LoadsKernel is chosen. Load-sKernel [13] is a comprehensive load analysis tool and in use for several configurations at the DLR-Institute of Aeroelasticity. In this work the VLM-part of the tool is utilised. This VLM follows the descriptions of Katz and Plotkin [14] closely. The lifting surfaces of all three variants are represented by a grid of 8 × 30 elements on each side of the wing. The angle of attack can be varied and a predefined twist distribution can be applied by an additional downwash.

Structure model
From the wingboxes of the given configurations, SE2A-MR and XRF1, Timoshenko beams are derived (see Fig. 5). The beam represents the wingbox stiffnesses with a shear factor of 0.83. The necessary moments of inertia I y , I z and I T are calculated from the derived wing box dimensions, which are given in Table 2. Consequently, the dimensions of the beams are dependent on the moments of inertia. To assure the same torsional stiffness, the correction factor c in Eq. 32 is varied.
The nonlinear elasticity is realised with stepwise variations of the E-modulus. (Note: This is done only for bending not for torsion.) The beam is clamped at its origin. While the mass distribution of the wing is not considered, the mass of the aircraft configuration (see Table 1) is represented by a point mass at the origin. For this reason, the corresponding MTOW-masses are fitted to match the maximum root bending moment and the wing tip deflection of more detailed models.

Coupling model
The calculated loads are attached to the structure by rigid body connections. To build up these connections, for each aerodynamic grid point the nearest grid of the structural model is searched for. The structural deformations are then returned to the aerodynamic model by affecting the downwash. This spline can be expressed as the transformation matrix T kf which relates the structural grid deformations with the aerodynamic grid deformations. Figure 4 displays (32) I T = c ⋅ ab 3 Fig. 4 Coupling models of the three configurations: red dots for the beam grids, blue dots for VLM-knots

Load analysis
To simulate different loadcases, the aircraft mass, altitude, speed, and load factor can be varied. The complete loads analysis is done by the loops in Fig. 6. An initial angle of attack init and the predefined load factor n z are given as input for the trim solver, which searches for the angle of attack needed to create the required lift of each loadcase. Included in this iterative process, an aeroelastic solver is composed of the VLM and the proposed nonlinear beam method (NLBM). This loop is run until the convergence in the wing tip deflection is constant ( w tip → 0 ). The outer loop is conducted until the predefined load factor n z is achieved ( n z → 0 ). To be more precise, the weight force at the centre of the wing multiplied by the load factor must be equal to the total lift produced by the angle of attack . For comparison of two different stiffness variants, the twist distributions of the cruise loadcase are adjusted by an additional downwash to create an equal lift distribution for both cruise flights.

Validation
In this section the proposed method is validated with numerical results of MSC.Nastran, with experimental results of a rubber-like material, and with a more comprehensive numerical solution of the experiment.
For the first validation with the beam analysis of MSC. Nastran, the rectangular wing (c) with a linear material and two multi-linear stiffness relations is considered. The stressstrain relations are visualised in Fig. 7. Figure 8 shows the resulting wing deflection for the three variations and for three loadcases (−1.0 g, 1.0 g, 2.5 g). For the linear case both methods obtain the exact same results. However, the nonlinear behaviour differs from each other. The degressive wing is less easing and the progressive wing is less stiff than calculated with the proposed method.
According to the Reference Manual, just the first and the last eighth longwise of each beam element are affected of the stiffness variation. Also, multi-linear stiffness relations are not recommended by MSC.Nastran in the nonlinear beam analysis with the element CBEAM [15]. Thus, the blue stiffness variation is overreacting and the red stiffness variation is underreacting. The two arrows show the changing direction and delta which can be calculated with the proposed method.   In a second and a third validation step the nonlinear part is validated with an experiment and a more sophisticated beam analysis. The experiment was conducted by Brojan et al. in 2012 [8]. They performed a cantilever beam experiment with a rectangular beam made of rubber-like material and loaded it with different weights at the free end. The results are shown for four of their ten different loading conditions in Fig. 9 with cross marks.
The more comprehensive beam method of Brojan et al. considers also the differences between tension and compression parts of the bending beam and its influence to the neutral axis. The dotted lines show the results of their numerical solution for such beams composed with linear (black) and nonlinear elastic material (green). In comparison, the results of the proposed method are given with the solid lines. In order to get reasonable results with this simpler method, a combined approximation of the stressstrain curves for tension and compression is utilised (see Fig. 10).manuscript.v7_proofread_VH3 The comparison of the beam deflection shows overall good results. The analytical results of both methods are matching up well. The nonlinearity of the material has to be considered to match the experimental data. The experiment and numeric results with deflections above 50% show up also geometric nonlinearities which are not considered in the proposed method, for both the linear and the nonlinear results. Slight differences between the deflections of the least loaded nonlinear beams can be seen, which can be explained with the approximated stress-strain curve which is linear up to strains of = 0.003 , due to its stepwise formulation.
Conclusively, the proposed method is valid for beams with nonlinear and linear stiffnesses up to 15%. Nevertheless, the results of the proposed method are still usable to get a good preview for deflections above 15%.

Results
The following section carries out the capabilities of the new approach, the discussion on the relation visualisation between curvature and bending moment, and the results for two wings for a simple use case of nonlinear stiffness design. The plot's x-axis, labelled with "y in m", is always the global y-axis of the wings.

Capabilities
With the new method any stiffness relation is possible, also multi-nonlinear elastic variants. To show up the possibilities of the new approach to its full extend, the stress-strain curves in Fig. 12 are adjusted to generate a specific behaviour around the cruise loadcases (1.0 g, dotted lines) and a differing behaviour for the manoeuvring loadcases (−1.0 g, dashed lines / 2.5 g, permanent lines). The results of these Fig. 9 Deflection of a rubber material beam in comparison between experiment [8], numerical results [8], and the proposed beam method Fig. 10 Stress-strain curves of a rubber material [8] and the approximated for the proposed beam method for Fig. 9 multi-nonlinear elastic stiffnesses are shown in Fig. 11. The deflection of the linear variant (black) are calculated for an equidistant set of load factors. Thus, these lines create an isometric visualisation (like in a contour map). In Fig. 11 it can be seen, that the blue variant is concentrating around the cruise loadcase and the red variant is clinched towards the manoeuvring loadcases.
The slight influence on the bending moment distribution is just visible at the intermediate loadcases (dash-dotted lines), because the deflection is designed to be the same for the loadcases 1.0 g and 2.5 g for all variants by the predefined stiffness relations.

Curvature-bending relation
In order to archive such relations and to design any required deflection of a wing structure, the view on the relation between and M b is more important than the relation between and . For optimised nonlinear behaviours this relation can show more accurately on which loading point M b the compared bending lines will cross each other.
For example, in Fig. 12 a fourth variant (light blue) is introduced to point out the difference between the stressstrain visualisation and the curvature-bending visualisation. In the stress-strain relation the blue and lightblue variants seem to lead to an equal deflection behaviour of the wing.
However, in the curvature-bending plot (in Fig. 12) the lightblue variant shows that its curvature is more alike the linear stiffness. For bending moments between 10.0 MNm and 17.5 MNm the linear and the lightblue variant would exhibit almost the same curvature. This has to be considered while optimising for a multi-nonlinear elastic structural behaviour which is best suitable for a configuration. The concentration of the wing deflection in Fig. 11 is consequently just achievable with such exaggerative variant in blue.

Softening
In the following, a simpler approach of a nonlinear elastic stiffness is considered. Now, for both wings a degressive behaviour is analysed. This leads to a softening of the wing structure for loadcases with a load factor above 1.0 g. Thus, the stiffness relations in Fig. 13 are split into two linear areas. Its E-modulus is reduced for strains above > 1200 m m for the mid-range wing and for strains > 800 m m for the long-range wing.

Mid-range wing (a)
The mid-range wing perfectly demonstrates the principles of the softening concept. The reduced stiffness leads to more deflection for higher loadcases (Fig. 16A). This deflection reduces the angle of attack at the wing tip due to the bending-torsion coupling (Fig. 16B). Consequently, this change of the twist distribution results in an inboard shift of the lift distribution (Fig. 16C), which is responsible for the root bending moment reduction (Fig. 16D).
The wing tip deflection of this 20% stiffness approach is 7.58 m (34.11%) for the 2.5 g manoeuvre. So, geometrical nonlinearity effects already have an influence. But these preview results still reveal the helpful tendency in the twist distribution with a lower incidence angle at the wing tip and a higher incidence angle at root, which is necessary for the required inboard shift. The lift rise at root is 2.8 kN higher resulting in a root bending moment reduction of -4.47%.

Long-range wing (b)
The long-range wing results are shown in Fig. 16E-H. Again, the principles of the easing wing above the cruise loadcase leads to a lift distribution shift and a root bending moment reduction. For the wing tip deflection of 5.89 m (19.6%) the geometrical nonlinearities can almost be neglected. So, the 7.6 kN lift rise at root yields a root bending moment reduction of -4.21%.
The twist distribution shift is a little bit stronger than the one of the mid-range wing, since the sweep angle is bigger. Thus, a lower wing deflection change is as effective as the wing deflection of the less swept wing, while the relative reduction is quite similar. This softening approach is conducted for a parameter analysis of sweep angle and stiffness reduction above cruise. This is presented in Figs. 14 and 15 . The most obvious result is that the root bending moment reduction is stronger for a higher stiffness reduction. It seems to be a quadratic relation for both wings. The change of the sweep angle shows more efficiency for a higher sweep angle, which is connected with the addressed bending-torsion coupling. (Fig. 16)

Conclusion and outlook
A new analytical method for nonlinear elastic beams is presented. It is implemented in a simple aeroelastic trim calculation and is tested and validated against numerical and experimental results. This new approach is capable to investigate any nonlinear elastic stiffness behaviour in the scope of aeroelasticity.
The view on the curvature-bending relation gives a better insight to these nonlinearities. Since this view is the integrated solution it is closer to its final result and creates a better preview for this kind of usage.
The root bending moment can be reduced for both wings with a degressive stress-strain relation. This softening approach has a good potential to reduce the loading   Fig. 16 Wing characteristics of the mid range wing and the long range wing with linear and nonlinear elastic stiffness for 11 loadcases between 2.5 g and −1.0 g with 20% of the initial stiffness above the 1g-cruise flight of design loadcases while the cruise flight shape is not affected.
Geometrical nonlinearities are not considered at this stage. In the experimental validation (Fig. 9) the results are overestimated for the deflections above 15% because of these geometrical effects. But a statement considering the bending moment reduction above 15% is not possible at this point, due to the contradictory effects of shortening and follower forces.
In connection with active load alleviation, it is difficult to say if both alleviation concepts can be superpositioned to its full extend, because both concepts deal with modifying the twist distribution. The advantages of this new passive concept are its automatic reaction to any loadcase and the possibility to strengthen the performance of cruise conditions. The disadvantages are the small experiences with this technology and the few possibilities of realisation for such stiffness behaviours.
In future, the process shall be more detailed and shall take the following aspects into account: wing mass distributions, kinks, and initial twist distributions. In order to also analyse the flight mechanical stability of such configurations, the dihedral of the wing and an empennage have to be modelled. Additionally, the flight performance has to be considered to better evaluate any kind of stiffness.
Eventually, the proposed method helps to get an impression of stiffness related nonlinearities and can be utilised to optimise the stiffness relation for aircraft structures.