Static aeroelastic stiffness optimization of a forward swept composite wing with CFD-corrected aero loads

This paper presents an aero load correction strategy applicable to the static aeroelastic optimization of composite wings. The optimization framework consists of a successive convex subproblem iteration procedure, in which a gradient-based optimizer consecutively solves a local approximation problem. Responses are approximated as a linear and/or reciprocal function of the laminate membrane and bending stiffness matrices. Together with the laminate thicknesses h, they constitute the design variables in the optimization process. Internally, the design space is transformed from stiffness matrices to lamination parameters, resulting in a continuous and convex optimization problem. Structural responses considered in the stiffness optimization are strength, local buckling and mass; aileron effectiveness, divergence, and twist constitute the aeroelastic responses. Steady aeroelastic loads are calculated with a doublet lattice method (DLM) embedded in the applied finite element solver, allowing for the generation of response sensitivities that incorporate the effects of displacement-dependent aeroelastic loads. To incorporate flow phenomena that cannot be reproduced with DLM, a higher order aerodynamic method is considered. The developed correction methods and their application are presented in this paper. The correction is twofold, first, aiming at a correction of DLM by means of camber and twist modifications applied directly to the doublet lattice mesh and second, by employing the capabilities of a higher order computational fluid dynamics (CFD) solver, like the DLR-based TAU code. To this end, DLM loads transferred to the structure are rectified by means of the supposedly superior CFD results. The aero load correction method is applied in the stiffness optimization of a forward swept wing. First, a trim application without structural optimization is discussed, to demonstrate the convergence behavior of the correction forces. The results of a wing skin mass minimization with balanced and unbalanced laminates are presented. In particular, the differences between optimizations with and without aero correction are highlighted. Eventually, a stacking sequence optimization based on the continuous optimization results is demonstrated.


Introduction
The application of composite material in load carrying structural components of an aircraft is rapidly gaining momentum. Three principal reasons for this can be identified. One, the continuous enhancements in the area of automated production technologies. Two, the progress in the development of computational methods to analyze and optimize composite structures to fully exploit their possible advantages. Three, the increasing confidence of designers in composite materials as a result of growing experience.
The identification of potential benefits to be achieved with composites, accordingly has entailed a considerable amount of research work. [39], for example, apply an approximationbased optimization to minimize the weight of a composite wing structure subjected to different combinations of buckling, strength, displacement and twist constraints. Although not optimizing, analyzing the effect of bending-twist coupling on the aeroelastic behavior of composites plates with various symmetric, unbalanced stacking sequences [17], demonstrate the effects induced by the coupling terms in 1 3 the bending stiffness matrix. [14] investigates the influence of non-symmetric laminates on the aeroelastic behavior of high aspect ratio wings by applying an integrating matrix method to derive divergence and flutter speeds of various stacking sequences. A valuable survey on aeroelastic tailoring effects as a result of the directional stiffness in orthotropic composites is provided in [37]. They investigate possible influences of rotated fiber angles on minimum weight, twist, normal modes, flutter and aerodynamic performance for various configurations, ranging from fighter aircraft to forward swept configurations. A general overview on optimization technologies ranging from optimizations on panel level to the aeroelastic optimization of composites in aircraft wings is provided in [42]. [10] investigate the influence of layup orientation in a straight fiber design on the optimized mass of a low aspect ratio wing, constrained by strength, roll-reversal and flutter velocity.
In a more recent work [27], maximize the flutter eigenfrequency and speed of a composite plate wing with ply angles as design variables, considering symmetric and non-symmetric stackings. A related research, but for a tow-steered composite plate featuring symmetric layups is described in [40]. The effect of varying tow angles on eigenmodes, elastic axis, as well as divergence and flutter speed is investigated. [15] present a two-stage procedure, where in a first step, layer thicknesses and angles serve as the design variables in minimizing mass, subject to strength and damage tolerance constraints. The second stage focuses on a reduction of the wing's response to a 1-cos-gust, measured by its transverse tip deflection.
Composite optimization with a focus on forward swept wings were first performed by [26], who demonstrated a considerable weight saving potential with the application of tailored stiffness distributions. [43,44] performed detailed investigations into the effect of the spanwise stiffness distribution and bending-torsion coupling on divergence velocity, aileron effectiveness and spanwise center of pressure of forward swept wings.
Aeroelastic tailoring studies based on lamination parameters have been attempted previously. [24] using a composite plate wing demonstrated the influence of lamination parameters on the flutter and divergence characteristics. They consider a parametrization of the bending stiffness matrix only, while implying symmetric laminates. [16] demonstrate a two-step optimization scheme to minimize mass of an aeroelastically loaded wing. In a first step, lamination parameters describing symmetric but potentially unbalanced laminates in the skins and spars are used. The second step comprises a genetic algorithm-based stacking sequence optimization. Minimization of compliance of a variable stiffness slender wing represented by beam elements is demonstrated in [1]. The cross-sectional properties of the beam are parametrized using lamination parameters that define the membrane stiffnesses of the box cross section. Another two-level optimization strategy has recently been proposed by [29]. Allowing for symmetric and unbalanced laminates, a lamination parameter-based weight minimization is performed before applying genetic algorithms for the derivation of stacking sequences.
In [8], the authors describe a stiffness optimization strategy based on lamination parameters and the consideration of static aeroelastic constraints. The strategy is adopted in the stiffness optimization of a parametrically defined set of forward swept wings, [7], featuring equivalent wing area and span and variable leading edge sweep. The influence of several static aeroelastic constraints on the minimized wing skin mass is investigated, in particular focusing on potential benefits of unbalanced over balanced laminates. The present paper describes an advancement of the optimization framework by means of a correction of the aero loads adopted in the optimization, as well as a stacking sequence optimization following the continuous stiffness optimization. Correction methods, with a strong focus on doublet lattice, have received a lot of attention. [13] suggest methods that require either a correcting of DLM pressures or modifications to the downwash based on the aerodynamic influence coefficients. [32] provided a survey of the available correction techniques, focusing on unsteady aerodynamic forces. An industrial application of DLM correction is presented in [4], who suggest a least square correction method that makes use of unsteady CFD results, and a twodimensional correction method. Both methods are aimed at correcting the aerodynamic influence coefficient matrix. In a more recent work, [9] compare gust responses obtained from DLM, CFD, and DLM that are corrected using quasi-steady CFD data.
The outline of the paper is as follows. Section 2 provides an overview of the continuous stiffness optimization framework. The aero load correction methods are introduced in Sect. 3. They are twofold, one, an adaptation of the camber line represented by the doublet lattice method that is integrated in the applied finite element solver, Sect. 3.1, and two, a correction of aero loads by means of a higher order CFD method, Sect. 3.2. An introduction to the stacking sequence optimization is provided in Sect. 4. Section 4.1 describes the methodology of stacking sequence tables, Sect. 4.2 introduces an interpolation strategy for response approximations. Results for the optimization of a forward swept wing are summarized in Sect. 5, starting with a model description, Sect. 5.1, a basic trim analysis, Sect. 5.2, a mass minimization, Sect. 5.3, and a stacking sequence optimization, Sect. 5.4.

Continuous stiffness optimization
For the optimization of composite structures, the problem is split into a three-step process, as detailed in [18], allowing for efficient use of suitable methods in each step. The first 1 3 gradient-based design step targets optimal stiffness distribution using laminate stiffness matrices through lamination parameters as the design variables.
Step one, which constitutes the central topic in the previous works on aeroelastic tailoring of composite wings, [7,8], will briefly be introduced in this section. The second step targets an optimal stacking sequence design, with the optimization seeded from the design obtained in the first step. Details will be provided in Sect. 4. The third step deals with retrieving optimal fiber paths for manufacturing and is not part of this research.
An overview of the step one optimization process is depicted in Fig. 1. The boxes featuring a gray background constitute the original process without aero correction. It corresponds to a successive convex subproblem iteration procedure, in which a gradient-based optimizer consecutively solves a local approximation problem.
The process is geared towards the optimization of the load carrying, shell-like structural components in a wing box, the properties of which can be represented as membrane A and bending stiffness matrices D, and shell thickness h. Considering only symmetric laminates, the coupling stiffness matrix B is zero and does not have to be considered.
The optimization is based on a finite element (FE) model of the wing structure that serves as an analysis model for the desired responses f, and for the evaluation of sensitivities f of the responses with respect to the design variables x. In this research, the NastraN finite element solver is applied. The DLR-Institute of Aeroelasticity in-house tool ModGeN [25] is used as a pre-processor to parametrically define and generate the FE model of the load carrying wing box, the aerodynamic and coupling model, as well as fuel and mass models. The responses and sensitivities serve as input for the derivation of an analogous analysis model that describes the behavior of each response in the surrounding of the analyzed design. For this purpose, each response is approximated as a function of potentially each design variable, while satisfying the essential properties of convexity, separability and conservativeness, as shown in the below equation: The sensitivities generated with NastraN are converted to linear, , and reciprocal, , sensitivities with respect to the stiffness matrices, where superscripts [⋅] m and [⋅] b denote sensitivities with respect to membrane and bending stiffness, respectively. j is the sensitivity with respect to the thickness design variable. The approximation model replaces the actual analysis model in the search of the optimizer for a minimum of the objective function f 0 , greatly accelerating the function evaluations required during the optimization. Within the optimizer, stiffness matrices A and D are represented as functions of lamination parameters, along with consideration of their feasible region as described for example in [36].
Each optimization step results in a new set of design variables that represent the global optimum of the convex approximation subproblem. Convergence is monitored by determining the change in the objective function for subsequent feasible iterations, Eq. (2). In case the change drops below a prescribed value stop , the optimization is stopped, otherwise the process continuous with the derivation of new responses and sensitivities: A detailed description of the optimization process and the components involved is provided in [6]. Aside from regular structural responses strength, buckling and mass, the aeroelastic responses aileron effectiveness, twist and divergence are considered in the optimization process. While their approximations are summarized in Table 1, the reader is referred to [6][7][8] for details on their derivation.
The process depicted in Fig. 1 is based on aeroelastic loads generated with the NastraN internal doublet lattice method. Methods for correcting those loads will be discussed in the following section.

Aero correction methods
The reason to perform a correction of the aeroelastic loads is the potential differences in wing surface pressures between the doublet lattice method developed for subsonic flows and presumably more correct higher order aerodynamic methods, usually designated as computational fluid dynamics (CFD). The latter types of methods allow for the consideration of flow phenomena that cannot be reproduced with DLM. Among the flow phenomena that should be considered the most important are potential flow singularities, distributed on a flat mesh, are superimposed with the undisturbed free stream. The strength of the singularities as a result of a linear system of equation has to be such, that the flow tangency at each box making up the mesh is fulfilled. This finally leads to a so-called aerodynamic influence coefficient matrix AIC that relates the pressure p in a box to the local downwash angle w of each box: Along with a surface integration matrix S, the aeroelastic loads in each DLM box a can be related to the pressures, Eq. (4): and ultimately to the aeroelastic loads in the structural coupling nodes DLM by means of a coupling matrix H, Eq. (5). Aero correction methods modifying either Eqs. (3) or (5) will be introduced in the following Sects. 3.1 and 3.2.
It should be noted that this paper does not necessarily focus on the most correct aerodynamic modeling strategy rather than on the provision of a general methodology which will allow for the consideration of all kinds of CFD-based aerodynamic fidelity levels in the aeroelastic optimization. To demonstrate the functionality of the method while maintaining reasonable modeling and computational effort, an Euler solver is applied in the present work, noting its shortcomings especially with respect to point 3 in the above listing of flow phenomena to be addressed. More details on this will be given in Sect. 3.2.

Camber correction
It is possible in NastraN to influence the geometric downwash w in Eq. (3) by a so-called W2GJ correction matrix. A detailed investigation of this contribution is provided in [23]. To this end, an individual downwash can be defined for each DLM box, corresponding to a local angle of attack change, as illustrated in Fig. 2. The required box rotations (4) a = aj j , Table 1 Response for the emulation of a camber line are shown in Fig. 2a. The chordwise constant rotation of each DLM box as shown in Fig. 2b is used to emulate a twist of the wing section. Both downwash types described in Fig. 2 can be varied in spanwise direction, allowing for the simulation of cambered airfoil blending and a geometric twist distribution. The parametric FE model generator ModGeN by default provides three W2GJ correction matrices for camber, twist, and the combination of both, which are generated based on the wing planform and the airfoil data provided in the ModGeN input file. Airfoil and planform data suffice to span the underlying aerodynamic surface.
In Sect. 5, it will be demonstrated that by a sole correction for airfoil camber and twist, an adequate agreement of pressure distributions computed with DLM and higher order aerodynamics can be achieved. This, however, only holds for lower, recompression shock-free Mach numbers. In the presence of shocks, the necessity for a correction by means of CFD is inevitable when aiming at an improvement of aeroelastic loads. Such a correction will be discussed in the following section.

CFD correction
Computational fluid dynamics, CFD, denotes methods used to solve the governing equations of a fluid flow. While the term is usually applied in conjunction with higher order volume mesh methods like Navier-Stokes solvers, strictly speaking it also applies to lower order methods. In the scientific community, as in this work, its meaning is dedicated to higher order methods exclusively.
Decision was made in favor of the DLR-German Aerospace Center unstructured CFD solver tau, [11,35]. It can be applied either for solving the full Navier-Stokes equations with a dedicated boundary layer and turbulence modeling, or in a simplified version, in which the viscous terms in the flow equations are neglected, hence resulting in Euler equations. Despite its inability to account for skin friction and flow separation-except at sharp corners such as trailing edges-as well as a biased shock prediction for strong shocks, [22,45], the Euler solver is applied in the present work. The main reasons are the considerable time saving compared to a full Navier-Stokes calculation, along with faster mesh generation by avoiding the necessity for a dedicated boundary layer mesh. With respect to the consideration of recompression shocks, Euler constitutes a reasonable compromise between computational effort and the gain in accuracy of the aerodynamic loading, well suited for the application in the structural optimization process. It should be noted, however, that capturing viscous effects is still out of scope and requires the application of higher order methods.
The basic idea behind the entitled CFD correction consists of rectifying aerodynamic loads obtained using the doublet lattice method by means of the supposably superior CFD results. To this end, the relevant sizing load cases are analyzed with tau, considering an appropriate volume mesh deformation that resembles the structural displacements. With the doublet lattice forces concentrated onto the coupling nodes, compare Eq. (5), the same nodes are selected for splining the surface forces obtained from the Euler calculation. The difference between the CFD force vector and the DLM force vector at each coupling node is applied as a static amendment to the respective load case. The force vectors on a deformed wing structure are shown in Fig. 3.

Process
The process of correcting the aerodynamic loads obtained with DLM constitutes an addition to the existing optimization framework [8]. In Fig. 1, the CFD correction part is already incorporated. It is defined as a stand-alone module that collaborates with the optimization via well-defined interfaces.
For each load case that is to be corrected using CFD aero loads, the necessary input is generated along with the Nas-traN runs required for sensitivity derivation and passed to the CFD correction module. The input consists of structural displacements at the coupling nodes, utilized for the CFD mesh deformation, and the flow parameters Mach number, stagnation pressure, density and the angle of attack. Two important notes have to be addressed at this point. One, the Euler calculation is performed for the same angle of attack as identified in the NastraN DLM trim maneuver analysis. Two, the displacements passed to the load correction process and used in the mesh deformation routines in the beginning are not necessarily converged displacements. Only throughout the structural optimization process, the loads and displacements eventually will converge, as will be shown below and be supported by computational results in Sect. 5.3.
Only symmetric maneuver and cruise load cases are designated for CFD correction, the first due to their potential to drive the structural sizing and the latter due to the evaluation and possible constraining of wing twist. In both cases, the corrected aero loading can have a significant influence on the results. The reason to consider only symmetric load cases is twofold: one, the currently limited analysis capabilities confined to symmetric steady cases and two, in lack of further information available, the effort to capture the largest loads by means of only a few load cases, see for example [6].
Once the CFD forces are generated and condensed to the coupling nodes, denoted as in Fig. 1, they are subtracted from the appropriate DLM forces DLM , yielding the correction forces Δ : where superscript k indicates the structural iteration steps for which tau corrections are generated.
The correction forces are grouped by load case and saved in the appropriate NastraN FORCE card format. A new optimization process is initiated, this time including the correction forces during the responses and sensitivity generation. They will have a direct influence on responses such as strains and displacements, and only an indirect influence on sensitivities via the forces influence on the aeroelastic behavior. The loop of optimizing and computing new correction forces is continued until an overall convergence is achieved.
Investigations have shown that a recalculation of the correction forces is not necessarily required during each, but only every nth iteration step. Details will be provided in Sect. 5. Usually, an aeroelastic trim calculation in NastraN ensures that the lift forces generated by the doublet lattice model exactly balance the weight vector multiplied with the vector defining the load factor. Performing a trim calculation along with the additional correction forces implies a force distribution at the coupling nodes in the trimmed solution, which exactly matches the CFD results. Since this only is true if the structural properties and, therefore, the stiffness remains unchanged, this statement strictly speaking does not apply during the optimization, where stiffness properties deliberately are modified. The reason is that during an optimization sequence, the correction forces remain constant and will only be updated after n iterations performed in the optimization block. Details on the influence on the load displacement behavior will be addressed in more detail in Sects. 5.2 and 5.3.
In the following, an examination of the convergence behavior will be given to demonstrate the acceptance of this approach. Setting up the static equilibrium equation for the (k + n)th iteration leads to where A is the aerodynamic stiffness matrix, which links structural displacements to the aerodynamic forces generated by DLM according to Superscript k denotes the last structural iteration step for which tau correction forces were generated, and accordingly superscript n is the n th subsequent structural iteration step. This notation implies that tau corrections not necessarily have to be performed for each structural iteration. Vector ie represents inertial and external forces. Inserting Eqs. (6) and (8) in (7) yields where Equation (10) states that in a converged solution where the residual force Δ res vanishes, the static equilibrium is determined entirely by means of tau aerodynamic forces , keeping in mind that according to Eq. (6) the aerodynamic force vector is a combination of doublet lattice forces and correction forces, = DLM + Δ . At this point, it should be stressed that the converged aeroelastic deformation complies with the tau aerodynamic forces, and likewise also accounts for the displacement-dependent doublet lattice forces according to Eq. (8).

Implementation
The CFD correction module depicted in Fig. 1 takes over the task of computing Euler CFD forces at the coupling nodes, based on the input that is required to perform a tau calculation on a deformed mesh. The essential steps successively performed in the module are illustrated in Fig. 4. Equivalent to the coupling matrix applied in Eq. (5), a coupling matrix H to map FE node deformations onto the CFD surface mesh, Eq. (12), needs to be constructed. The transpose of H can also be used to achieve the second goal of transforming CFD forces onto the structure, Eq. (13): In this research, the coupling method by [31] is applied. The volume mesh deformation following the surface mesh deformation is part of the tau software suite. Details are described in [6].
The CFD meshes employed in this work were generated using the surface modeler and unstructured volume mesher suMo, [41]. Depending on the intended surface pressure resolution, result accuracy, in connection with, for example, the computed shock position, and the CFD convergence and computing time, the generated Euler meshes usually k+n s exhibit ≈ 200,000 to 300,000 triangular surface elements and tetrahedra elements ranging from ≈ 0.7 × 10 6 to 3 × 10 6 . suMo allows for a mesh export directly in the appropriate tau mesh format.

Stacking sequence optimization
In this section, the second stacking sequence design step is described using the method developed in [30]. It constitutes the next step subsequent to the continuous optimization process depicted in Sect. 2, which in turn provides the required stiffness matrix input and sensitivities. The stacking sequence design tool [30] consists of a stacking sequence table (SST)-based genetic algorithm (GA) for the design of blended structures, [21], combined with a modified Shepard's method for improving approximation accuracy in the GA, [20].

Stacking sequence table
In an efficient variable stiffness design, usually a finite number of panels each having constant stacking sequence are considered. Laminate blending [46] is a technique that assures continuity of fiber and material content between such adjacent panels having different ply layups, thus enhancing structural integrity and manufacturability. [19,38] present some of the successful attempts at achieving fully blended optimal designs. A stacking sequence table is an intuitive way of representing the ply layup of panels in a blended composite structure. Each column in an SST essentially corresponds to the stacking sequence of a certain number of plies. Starting with the stack featuring the smallest thickness, plies are added successively, following several design guidelines, up to the maximum number of plies. This ensures that the plies in the thinnest stack are carried on to the thicker ones, thus inherently satisfying the requirement of blending, without enforcing it as additional constraints.
In comparison to previous works on blended composite design, [2,19,28], the SST-based optimizer [21] offers the following advantages: -implementation of several industry-standard laminate and ply-drop design guidelines; -explicit information of ply-drop sequence among adjacent panels for manufacturability; -fully blended designs according to the generalized blending definition [5].
An advantage inherent to this approach is the compact and efficient form of encoding of the entire SST using only two chromosomes. SST lam defines the stacking sequence of the thickest laminate and SST ins defines the order in which plies are inserted from the thinnest to thickest plies. By suitably optimizing for these two chromosomes and the ply numbers in the panels, an optimal blended structure can be obtained. An example of an SST having ply numbers between six and ten is shown in Fig. 5.

Successive approximations
The concept of successive approximations as introduced in Sect. 2 helps to reduce computational costs related to expensive FE analyses. This is realized by carrying out the optimization on approximations of responses, Eq. (1), rather than expensive FE responses themselves. The optimization is iteratively followed by a design update and construction of new approximations for the next optimization loop. A similar successive approximation framework is used in this stacking sequence step also. The aeroelastic response approximations listed in Table 1 are applied also in the GA to evaluate the optimization objective and constraints. The stacking sequence optimizer presented in [30] addresses strength and buckling constraints using a modified approach. Rather than using direct approximations of strain and buckling as in [8], panel loads are directly approximated instead: where m and b are panel load sensitivities evaluated from NastraN. The panel loads are then used to evaluate strain and buckling responses using suitable analytical tools. The benefit of such an approach is that using a direct analysis of approximated loads compensates for the increase in computational time with a higher accuracy in the buckling and strain responses themselves. The approximations mentioned above are single-point approximations constructed at a design point and are accurate only at that particular point. A modified Shepard's method [20] constructs multi-point approximations by interpolating previously analyzed design points, thus improving the quality of the approximations. For instance, the aileron effectiveness response is evaluated using such a multi-point approximation as where ̃a il i is a single-point approximation of the aileron effectiveness response at the ith point as in Table 1, w i ( ) is a weighting function defined in [20] and n is the total number of approximation points. An improvement in the quality of the approximations naturally results in a faster convergence.

Optimization steps
1. At the optimal design obtained from the continuous stiffness optimization step, response approximations are evaluated. This forms the first design point. 2. The optimal solution to this approximate subproblem is then found using the GA for SST, wherein response approximations are used to evaluate the objective and constraints. In this step, approximations for mass, twist, divergence and aileron effectiveness are evaluated to obtain the respective approximate responses. Panel strength and buckling are estimated analytically from the panel loads, which are approximated from the load sensitivities. 3. An FE analysis at the optimal design point obtained in step (2) yields new response approximations at this point.

A multi-point approximation of the responses (as in
Eq. 15) is formulated using all single-point approximations obtained till the present iteration. Steps (2)-(4) are then sequentially repeated till convergence.
Two possible objective functions for the GA optimization can be utilized. In the first case, the objective function of the GA is directly chosen as the particular response to be minimized, for instance structural mass. In the second case, the objective function aims to match the stiffness of the GA Fig. 5 Example of SST with six to ten number of plies design to the optimal stiffness obtained in the previous continuous optimization step. The optimizer described here and in the ensuing results utilizes the former approach.

Results
A stiffness optimization with corrected aero loads, followed by a stacking sequence optimization, likewise featuring aero load correction, is demonstrated in this section.

Model description
With this work following up on the research presented in [7], the same reference model for demonstration of the developed processes was adopted. An overview on the analysis model is presented in Fig. 6. In Fig. 6a, the dimensions of wing and wing box are depicted. The corresponding FE model is shown in Fig. 6b, along with nonstructural masses in front and aft of the wing box and clamping at the root rib. The fuel model, featuring one concentrated mass per rib bay is shown in Fig. 6c. The DLM and aero-structural coupling model are depicted in Fig. 6d. The coupling nodes are based on the so-called load reference axis (LRA), which consists of a virtual axis in spanwise direction that is marked by grid points in each rib plane. The grid points are attached with RBE3 multipoint connections to the corresponding circumferential rib nodes. Each LRA grid point comprises two RBE2 rigid elements, one extending to the leading edge of the underlying planform and one to the trailing edge.
All the above-mentioned models were parametrically generated with the ModGeN pre-processor. A detailed description of the model generation process and the applied mass and material properties are given in [6].
Load cases applied in the optimization are listed in Table 2. They can be grouped into maneuver load cases used for sizing of the structure (1-4), a cruise load case (7), load cases for the determination of aileron effectiveness (12)(13)(14)(15), and a divergence load case (16). The load cases are combined with five mass cases resulting from variations in fuel and passenger mass. Masses for half the fuselage and tails were attached to the clamping node on the centerline. Presumably, the most unfavorable combination of empty wing tanks and maximum passenger load was considered for sizing load cases (1)(2)(3)(4). The cruise load case (7) was investigated for a wing fuel loading approximately corresponding to begin/mid, and end of cruise flight, and for maximum and half passenger loading, totaling four more load cases. The aileron effectiveness at a constant roll rate and divergence load case were independent of the mass distribution and could, therefore, be computed along with one of the depicted mass cases (further details on the responses can be found in [6]   The optimization model comprised in total 68 design fields, Fig. 7, where a design field specifies a group of elements that all refer to the same, unique set of stiffness matrices and shell thickness. Accordingly, the design variables are made up of 68 independent (A, D, h)-sets.
Approximations for strain and buckling failure, aileron effectiveness, divergence and twist were implemented according to the derivations summarized in Table 1. This included the convexification of aeroelastic responses, described in [6]. The strain allowables required for the failure envelope construction were set to [ t , c , xy ] = [0.5%, −0.4%, 0.4%].
The CFD mesh was generated with suMo and based on the same planform and airfoil coordinate input defined for the generation of the finite element model. To determine the appropriate mesh resolution that would be required to achieve convergence of the important aerodynamic quantities such as lift and moment, a convergence study was performed. The values of lift, moment and drag coefficient, C L , C m and C D , for six meshes featuring different tetrahedra numbers are presented in Fig. 8a. As stated by the figure, the coefficients converged for meshes with element numbers above ≈ 1.5 × 10 6 .
The mesh that was finally selected for the aero correction process consisted of ≈ 2.7 × 10 6 tetrahedra and is marked by a blue star in Fig. 8a. It comprised ≈ 315, 000 surface triangles, a spanwise section of which is shown in Fig. 8b, to give an impression on the element density and distribution as it was achieved using diverse suMo meshing parameters.
In response to the vast amount of significant results generated by trim calculation and stiffness optimization using aero correction forces, only the most interesting results and load cases are selected for presentation. They are chosen depending on the particular structural or aerodynamic focus. The reader is referred to [6] for an in-depth discussion of the findings.

Trim application
To separate the effects of a newly applied CFD aero correction from the effects induced by a stiffness optimization, the general functionality of the CFD correction module was demonstrated using a classical static aeroelastic trim application. In doing so, the most prominent differences between the results of the applied aerodynamic methods can be highlighted.
The iterative procedure to determine the aeroelastic deformation considering tau correction forces can be expressed as with where superscript k denotes the iteration step and ie corresponds to the summation of constant inertial and external forces. An iterative procedure which performs deformation and aerodynamic analysis in a sequential order is referred to as a weakly coupled system, whereas in a closely coupled system the aerodynamic forces can be expressed directly as a function of displacement, thus allowing for a direct solution of the static equilibrium equation.
According to Eq. (16), convergence can be tested by monitoring either a characteristic node deflection or the iterative behavior of the aerodynamic forces. This is shown in Fig. 9,   Table 2, where the last two digits identify the load case and the first digit the mass case.
The deflection shown in Fig. 9a is monitored at the outermost spanwise load reference axis grid point, representing the tip of the wing. In the finite element solution of iteration k = 1 , tau forces 0 were not yet included. They are generated only afterwards, for a CFD mesh deformation which is based on the displacement results of the FE analysis. Accordingly, only DLM forces affected the deflection for iteration k = 1 . From the second iteration onwards, tau forces based on the previous FE analysis were included. Figure 9a shows that once the correction forces are considered from iteration two on, convergence was achieved after three to four iteration steps.
The same holds for the monitored aero loads, Fig. 9b. Depicted in the upper axis is the force development for DLM and tau, respectively. | aero | denotes the magnitude of the sum over all aerodynamic forces in the xz-plane. Aside from quick convergence, a second very important finding is the convergence of the resulting tau force towards exactly the same value as the DLM force in the first iteration step k = 1 , which did not, as yet, comprise the correction forces. This implies that the combination of DLM and correction Δ in the converged solution exactly reflects , as it was stated in Eq. (17), = DLM + Δ . The corresponding force magnitude reflects the constant lift force required to balance aircraft weight multiplied by the load factor.
The lower axes in Fig. 9b depict the relative differences between the DLM and tau force magnitudes. It was mentioned previously that the tau Euler calculation is performed for the same angle of attack as that resulting from the Nas-traN trim analysis. Thus, the converged difference between | DLM | and | | is based only on aerodynamic discrepancies among the two discretization and analysis methods.
To demonstrate the impact of the two aero load correction methods presented in Sect. 3, a comparison of pressure difference ΔC p between upper and lower surfaces for DLM and tau is plotted in Fig. 10.
A comparison for LC 1002, hence a sizing load case with n z = + 2.5 g and a rather low Mach number of M = 0.597 is shown in Fig. 10a. Evidently, the camber-corrected DLM was in good agreement with the higher order aerodynamic CFD method tau. In all the depicted sections, DLM reproduced the main trends, although slightly downstream of the tau results. With the DLM mesh featuring a considerably coarser discretization, no data were available at the trailing edge; the largest deviations, therefore, occurred at these positions. An entirely different behavior is revealed when looking at the second n z = + g sizing load case 1004, Fig. 10b. The free stream Mach number was considerably higher, resulting in a strong shock on the upper and on the lower wing surface and thus a considerable influence on the pressure difference ΔC p . As a result of the large negative pressure difference in the back, and a decrease in the frontal part with respect to the doublet lattice results, the tau distribution was expected to cause a larger nose-down, thus negative twisting moment.

Numerical results: mass minimization
Having demonstrated the functionality of the aero correction process in a pure trim analysis, the application and interaction within the stiffness optimization process were tested. While the wing skin mass was defined as the objective to be minimized, the structural constraints again comprised strain and buckling failure for all shell elements belonging to the specified design fields. In terms of aeroelastic responses, only aileron effectiveness and divergence were constrained, while no constraint was defined for twist. Aileron effectiveness for the four load cases listed in Table 2 was limited by a lower bound to ail ≥ 0 , thus no aileron reversal. Divergence was constrained to a lower limit of q div ≥ 35, 000 Pa at a Mach number of M = 0.87 , ensuring a reasonable clearance with respect to the aeroelastic stability margin (1.15 times the dive Mach number/velocity). To confirm the mass minimization results, several variations of the starting design, comprising modifications on thickness and stacking sequence, were optimized in parallel. Except for some exceptions, all the starting designs led to the same optimum in terms of minimum mass and optimized thickness and stiffness.
The minimized masses of the investigated combinations with and without aero correction as well as balanced and unbalanced laminates are listed in Table 3.
As illustrated by the table, consideration of the aero correction did not imply a fixed impact on optimized mass, given by the fact that the mass increased for balanced, and decreased for unbalanced laminates. Anyway, the intention of the aero correction lies in the improvement of sizing loads and the enhanced determination of the aeroelastic responses. Owing to the considerable mass savings of unbalanced over balanced laminates, the sections focus will be placed on the results attained with unbalanced laminates.
In a first step, the aero force development during the stiffness optimization was reviewed to ensure convergence of the structural properties, and of the correction forces. The results for two sizing load cases are depicted in Fig. 11. In the optimization process depicted in Fig. 1, a tau correction run was requested every five structural optimization steps. The dashed blue lines in Fig. 11 indicate a tau correction run at the 4th , 9th , 14th … iteration step. Accordingly, the new correction forces were only available for these iterations, while DLM forces were generated during each iteration step. The graphs state that the tau and DLM forces for all sizing load cases converged, while in parallel, the optimization process minimized the mass objective by modifying the stiffness properties and thus the aeroelastic behavior. The trim application already proved a fast convergence with constant structural properties, and the stiffness changes during the optimization did not considerably deteriorate the convergence behavior. Nevertheless, Fig. 11 reveals a more gradual correction force change during the optimization, compared to the already good agreement seen with respect to the final state for the second tau correction step in the trim application.
The optimized thickness and stiffness distribution are plotted in Fig. 12, where Ê 11 ( ) = 1∕Â −1 11 ( ) is the polar thickness normalized engineering modulus of elasticity. It allows for a visual assessment of the directional membrane stiffness distribution. In search of a weight optimal solution, bending-torsion coupling was introduced by tilting the major stiffness direction from inner to outer wing gradually forward, thus leading to a negative twisting tendency when bending the wing up. Thereby the center of lift could be shifted inward, eventually reducing the root bending moment.
The optimization with aero correction led to an active divergence pressure constraint, while it was not active when optimizing with DLM only. Divergence pressure as developing throughout the iteration steps is shown in Fig. 13a. The aileron effectiveness remained clearly in the feasible domain; the lowest, still inactive response is shown in Fig. 13b.
More prominent differences when applying aero correction were revealed when looking at the load cases that accounted for the highest strength or buckling failure index in each design field, Fig. 14. Plots in the left column depict upper and lower skins of the optimized model including aero correction (denoted "DLM & tau "), the equivalent, but without aero correction (denoted "DLM"), is shown in the right column. While wing skin sizing for the optimization without aero correction was clearly dominated by n z = + 2.5 g LC 1002 and n z = − 1.0 g LC 1003 and only a few design fields being sized by LC 1004 and none by LC 1001, the allowance for aero correction led to a perceptibly different distribution. Distinct spanwise and chordwise regions developed, comprising all sizing load cases considered in the optimization. It should be noted that the distributions shown in Fig. 14 are only an indicator of the dominant load cases. Even a slightly lower failure index in an element resulting from another load case is not represented in this plot.
The differences in twist for an optimization with aero correction compared to an optimization without is depicted in Fig. 15, in which the spanwise twist distribution for two high Mach number sizing load cases 1003 and 1004, and for two representative cruise load cases 2007 and 5007, featuring the  largest mass variation among the investigated cruise conditions, are shown. Unexpectedly, the n z = + 2.5 g pull-up maneuver LC 1004 in both optimizations showed a negative tip twist, differing by ≈ 2 • . The reason for this could be found in the bending-torsion coupling evoked by the variable stiffness orientation, Fig. 12b. The twist distribution promoted compliance with the divergence pressure constraint and helped to alleviate loads in the outer wing, thus supporting mass minimization. The twist being more negative when considering aero correction could be attributed to the different structural designs and to the more negative aerodynamic airfoil moment. As a result of the superposition of geometric coupling of the forward swept wing, and the negative aerodynamic twisting moment, and despite the lower inertial and, therefore, aerodynamic loading in case of the n z = − 1.0 g push-down maneuver LC 1003, the wing twisted considerably more negatively for LC 1003 compared to LC 1004. Again, the difference between the aero-corrected design and non-corrected design could mainly be attributed to the different aerodynamic moment distribution.

Numerical results: stacking sequence optimization
The results from the stacking sequence optimization as described in Sect. 4 are presented here. The optimization hence constitutes the second step in the three-step design process. The results for the forward swept model are introduced above, subjected to the same load cases, objectives and constraints. The continuous stiffness optimization presented in the previous section showed an optimal wing box mass of 403.9 kg with unbalanced laminates. This mass represents the theoretical upper bound in performance that can be achieved since a lamination parameter-based continuous optimization does not impose any restriction on ply angles, ply thickness or laminate continuity. This continuous design constituted the first design point in the stacking sequence optimization, helping to efficiently find the region of interest for the GA to search. For this step, the wing was modeled to comprise six independently blended regions-hence, six independent SSTs. This is more appealing from an industrial perspective, wherein a large structure such as an aircraft wing is manufactured as segments before being joined together. Enforcing blending over the complete wing hence only constrains the design space by combining structurally separate entities. Shown in Fig. 16 is the distribution of the independently blended regions over the upper wing surface. A similar distribution is also used for the lower surface.
The GA parameters and the design guidelines enforced are presented in Table 4. The optimal stacking sequence design was found to have a wing mass of 485.9 kg, or ≈20% higher mass than the stiffness-optimal design. Importantly, the optimal stacking design was feasible, i.e., having a safety factor for the design constraints of exactly 1.0. The increase in mass over the stiffness-based design can be attributed mainly to the following reasons: allowable thickness only in discrete steps of the ply thickness, discrete set of allowable The optimized thickness and stiffness distribution are plotted in Fig. 17. A direct comparison of Fig. 17b with the continuous optimization counterpart, Fig. 12b, identifies an equivalent tendency of the wing to have a forward-tilted stiffness, indicating the dominant direction of the ply angles. The discrete design, however, shows a rather smeared distribution of the stiffness with a less articulate directional alignment.
A comparison of the thickness in Fig. 17a with Fig. 12a shows a more gradual decrease in the thickness of the plies from root to tip in the stacking sequence design. Some of the elements, in fact, have a lower thickness than in the theoretical continuous optimum itself. The above observations are a consequence of the blending requirement-each panel in a blended design is influenced by the ply layup in all the other panels in that blended region. This results in a smooth and 'smeared' distribution of the thickness and ply angles, resultantly in the stiffness itself.

Conclusion
A detailed insight into the application and implications of a doublet lattice force correction using a higher order CFD method in an aeroelastic stiffness optimization was presented.
The trim application demonstrated the general functionality and convergence of the correction procedure, and highlighted the differences to be expected when using the two aerodynamic methods tau and DLM. The improvements made in the doublet lattice method due to using a W2GJ camber correction and the limitations of this correction with the emergence of recompression shocks were illustrated. It was shown that the consideration of camber correction did greatly improve DLM quality when applied in shock-free conditions. The effects of aero correction consideration on mass minimization, emphasizing the differences with respect to optimizations that did not feature an aero correction, were highlighted.
The applied Euler method proved to converge reliably with the mesh generated in suMo. Nevertheless, it should be mentioned that the Euler method has a limited application range. With increasing angle of attack and thus lift coefficient, Euler predicts the shock to move more and more downstream, increasing in strength. Flow separation, other than at the sharp trailing edge, cannot be modeled.   Accordingly, in the case of severe aerodynamic load conditions, the Euler results will start to deviate from what can be expected in reality. A possible solution to this problem is to increase CFD fidelity further, and thus consideration of the full Navier-Stokes equations along with viscous boundary layers and turbulence modeling. Apart from the need for a new CFD mesh topology including a prismatic sub-layer, the required changes to the developed CFD correction module are minimal.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.