Coupled adjoint aerostructural wing optimization using quasi-three-dimensional aerodynamic analysis

This paper presents a method for wing aerostruc-tural analysis and optimization, which needs much lower computational costs, while computes the wing drag and structural deformation with a level of accuracy comparable to the higher fidelity CFD and FEM tools. A quasi-three-dimensional aerodynamic solver is developed and connected to a finite beam element model for wing aerostruc-tural optimization. In a quasi-three-dimensional approach an inviscid incompressible vortex lattice method is coupled with a viscous compressible airfoil analysis code for drag prediction of a three dimensional wing. The accuracy of the proposed method for wing drag prediction is validated by comparing its results with the results of a higher fidelity CFD analysis. The wing structural deformation as well as the stress distribution in the wingbox structure is computed using a finite beam element model. The Newton method is used to solve the coupled system. The sensitivities of the outputs, for example the wing drag, with respect to the inputs, for example the wing geometry, is computed by a


Introduction
Selection of the fidelity of analysis in a complex multidisciplinary design optimization (MDO) such as wing aerostructural optimization is a challenge.Using high fidelity models makes the MDO problem computationally intensive and in many cases impossible to solve without the use of parallel high performance computational resources (Kenway and Martins 2014;Martins et al. 2004).This is a serious barrier against using high fidelity optimizations in early design stages, where many different designs have to be evaluated and optimized.On the other hand, using lower fidelity analysis has its own drawbacks.Lower fidelity methods sacrifice the level of accuracy and design sensitivity to achieve lower computational cost.For example empirical wing weight estimation methods usually compute the wing weight just based on the wing maximum thickness and do not take into account the effect of the whole airfoil shape on the wing weight.Besides, low fidelity methods cannot capture unconventional designs and their simplifying physical assumptions may break down during optimization.
According to industry criteria (van Dam 2003) a drag prediction with accuracy of one drag count (one ten thousandth of the drag coefficient) is required.The need for such a high level of accuracy is confirmed by Meredith (1993), where he showed that one drag count is equal to the weight of one passenger for a long-haul aircraft.Although achieving such a high level of accuracy for drag estimation in the conceptual design phase seems impossible (since many key elements of drag such as interference drag, power plant installation drag or excrescence drag are still missing), the first step toward reaching this goal is to replace the semiempirical methods such as ESDU and DATCOM with more physics based analysis methods as early as possible in a design process.
Aerodynamic solvers with different levels of fidelity have been used for wing aerodynamic and aerostructural optimization.Kennedy and Martins (2010) used a panel code for aerodynamic analysis in a wing aerostructural optimization.Since panel codes are not able to predict viscous and wave drag, they applied the aerostructural optimization to only minimize the induced drag of a subsonic passenger aircraft.The same panel code was used by Liem et al. (2013) for optimization of a transonic wing.A combined use of panel code, for induced drag estimation, and semiempirical methods, for viscous and wave drag estimation, was used by Kennedy andMartins (2014, 2012).Piperni et al. (2007) used a three-dimensional transonic small disturbance (TSD) code coupled to boundary layer calculation for wing aerostructural optimization in the transonic regime.In general, TSD codes are reliable for drag estimation at transonic conditions with relatively weak shock waves and attached flow.A higher fidelity Euler code is used by Maute et al. (2001) for wing aerostructural optimization.However they did not implement any method (semi-empirical or boundary layer method) for viscous drag prediction.Barcelos and Maute (2008) compared the results of the optimization using Euler and Navier-Stokes flow models to investigate the importance of accounting for viscous effects.They concluded that a general idea about the overall layout of the optimum wing can be achieved by an optimization using an inviscid flow model, however the viscous effects need to be taken into account for fine-tuning the design and for obtaining reliable optimization results.It should be mentioned that they only used aerodynamic lift and drag in the objective function.If the objective function is constructed based on both aerodynamic (drag for example) and structural properties (weight for example) an inviscid formulation cannot be used, because it does not provide correct value for the total drag, that negatively affects the compromise between the drag and the weight.Finally (Kenway et al. 2014) performed aerostructural optimization using RANS code for minimizing aircraft fuel burn during cruise flight.That can be counted as the highest level of fidelity applied so far for aerostructural optimization.
In all the methods reviewed above, a three-dimensional wing was analyzed for drag prediction.Generating and deforming a three-dimensional mesh, as well as solving a three-dimensional domain using CFD is extremely time consuming and requires high performance, parallel computing systems.In contrast, two-dimensional airfoil analysis can be executed much faster and cheaper.An interesting way to compute wing drag with sufficient level of accuracy and low computational cost is to combine twodimensional viscous airfoil data with an inviscid threedimensional wing lift calculation.This approach is named quasi-three-dimensional (Q3D) analysis.Examples of Q3D wing aerodynamic analysis and optimization can be found in the works of Drela (2010a), van Dam et al. (2001), Elham and van Tooren (2014a, b), Mariens et al. (2014), Willcox and Wakayama (2003) and Jansen et al. (2010).
The Q3D approach for drag estimation, which can be counted as medium level of fidelity, is a very useful technique for aircraft design and optimization in early design stages.Using this approach the aerodynamic characteristics of an aircraft can be estimated with much higher accuracy than semi-empirical methods, while the computational time is much lower than high fidelity three-dimensional analysis.
In this paper a Q3D aerodynamic solver is connected to a finite element based structural solver for wing aerostructural analysis and optimization.The coupled tool is able to compute the derivative of outputs with respect to the inputs using analytical methods.That ability makes the tool suitable for optimization using gradient based algorithms.This tool can be integrated with other aircraft design disciplines such as flight dynamics and performance for aircraft optimization in conceptual and preliminary design steps.

Aerodynamic analysis
In general there are two ways to compute drag of a body using CFD analysis.The first way is called near field analysis (van Dam 1999), in which the drag is computed by integrating the pressure and the friction forces around the body.In that way the drag includes two components: the pressure drag and the friction drag.The second way for drag computation is called far field analysis (Meheut and Bailly 2008;Gariepy et al. 2013).In a far field analysis drag of a body is computed by analyzing the inflow and outflow of a control volume around the body.Viscous drag, vortex drag and wave drag are the drag components that can be computed by a far field analysis.
As mentioned earlier a Q3D aerodynamic analysis is used for wing drag computation.The concept of the Q3D aerodynamic analysis presented in this paper is based on the Q3D analysis method presented by the author (Elham 2015), with some additions and modifications.The modifications are mainly applied to connect the aerodynamic solver to a finite element analysis tool for aerostructural optimization.In this Q3D wing aerodynamic analysis the lift distribution on the wing is computed using a Vortex Lattice Method (VLM), and then the strip theory (Flandro et al. 2012) is applied to compute the wing viscous drag in several spanwise positions.
Figure 1 shows the force and angles at a typical wing 2D cross section.From this figure one can observe that the total drag (parallel to the free stream velocity) is a function of both the effective lift and the effective drag as well as the downwash angle (α i ).The effective drag of the 2D section is the sum of the pressure and friction drag of the section at the effective angle of attack: The first component of the drag is the pressure drag caused by the shape of the airfoils.This drag component is also known as the form drag (Torenbeek 2013).The second component is the friction drag.The sum of the form drag and the friction drag is called the parasite drag (Torenbeek 2013).The third component of the drag in ( 1) is in fact the drag caused by tilting the lift vector due to the downwash angle resulted from the wing tip vortex.This drag component is known as the induced drag.Therefore based on (1) the drag of a wing in a Q3D analysis can be decomposed into the form drag, the friction drag and the induced drag.It should be noted that in a 3D wing drag computation using near field analysis the induced drag is included in the pressure drag, so the wing drag consists of the pressure drag and the friction drag.However in a Q3D drag analysis, the section pressure drag does not include the induced drag, so the wing induced drag is counted as a separate (third) component of the total drag.Therefore the total wing drag consists of the parasite drag and the induced drag.The Q3D approach used to compute the wing total drag is shown in Fig. 2. Each steps shown in this figure is explained in more details in the followings.
In the first step the lift distribution on a wing is computed using a VLM.A VLM code has been developed based on the method presented by Katz and Plotkin (2001).In such a method a ring vortex is placed around each collection point.The collection points are placed at the center of three-quarter chord lines of the panels and the leading segment of the ring vortices are placed at the quarter chord lines of the panels.The wing is followed by free wake vortices starting from the trailing edge.In order to take into account the effect of airfoil shape on the wing loading, the boundary conditions are applied on the wing camber line.The aileron deflection is simulated by rotating the vortex and the collection points that are placed on the aileron.That is used to compute the aileron effectiveness, see Section 3.
Using the wing geometry and the angle of attack, the Aerodynamic Influence Coefficients (AIC) matrix and the Right Hand Side (RHS) vector are computed.Then the strengths of the vortex rings (Γ ) are computed using the following equation: The wing lift distribution is calculated by using the Kutta-Joukowski theorem based on Γ .The results are corrected for compressibility effects at high Mach numbers using Prandtl-Glauert compressibility correction.The wing induced drag is also computed using the Trefftz plane analysis.
In the second step of Q3D drag computation, the wing is divided into several sections for 2D analysis.The flow properties at each section can be determined from three dimensional flow properties using two steps of transformation.The first step of transformation is performed based on the sweep theory (Holt 1990) to find the airfoil geometry (y coordinate for normalized x between 0 and 1) and the flow characteristics perpendicular to the sweep line: where Λ is the sweep angle and in (6) is the wing local geometrical twist angle.The value of C l in ( 7) is determined by interpolating the spanwise lift distribution from the VLM analysis for the given spanwise position.In the subsonic regime the wing quarter chord sweep angle can be used in Fig. 2 Steps of the quasi-three-dimensional approach for computing wing total drag (4) to (7) (Obert 2009).However for the transonic regime a sweep line that coincides with the shock wave should be used (Drela 2010a;b) because the pressure drag acts perpendicular to the isobars (or shock wave line).Therefore in transonic regime half-chord sweep angle is a better choice than the quarter chord sweep angle.The second step of transformation is performed to determine the airfoil effective angle of attack, Mach number and Reynolds number at each strip.Those data are required to perform 2D airfoil analysis.
It should be noted that in order to compute the effective flow properties the wing angle of attack as well as the downwash angles at each strip are required.The method used to compute them is explained in Section 4. Using the effective properties, as well as the airfoil geometry, the airfoil effective pressure drag, friction drag and lift (see Fig. 1) can be computed using an airfoil analysis tool such as MSES (Drela 2007): MSES is an interactive viscous/inviscid Euler method that features the design and analysis of single and multi-element airfoils at low Reynolds numbers and transonic Mach numbers.In addition, MSES can also predict flows with transitional separation bubbles, shock waves, trailing edge, and shock-induced separations (Drela and Giles 1987).Different parameterization methods are implemented in MSES for airfoil geometry modeling.In this study the Chebychev polynomials are used to parameterize the airfoil geometry.
Using that method the airfoil geometry perturbation normal to its current surface (Δn) is determined based on the basis functions g j , which are the Chebychev polynomials, and the mode amplitudes G j that are defined as design variables: where s is the fractional arc length of each side of the airfoil.MSES uses analytical methods to compute the derivatives of the outputs (lift, drag etc.) with respect to the inputs, including airfoil geometry, angle of attack, Mach number and Reynolds number.That ability of MSES is used to compute the sensitivity of the wing drag with respect to the wing geometry, see Section 5.
Using (1) the airfoil total drag parallel to V ⊥ (normal to the sweep line) can be calculated based on the airfoil effective drag: The term cos 2 α i in the denominator is due to the transformation of the force into the coefficient.The last step is to use the sweep theory again to calculate the wing drag parallel to the free stream velocity using the drag coefficient perpendicular to the sweep line.According to Drela (2010a, b) it is reasonable to assume that the friction drag scales with the free stream dynamic pressure and acts mostly along the free stream flow direction.On the other hand, the pressure drag from the shock and viscous displacement is assumed to scale with the wing normal dynamic pressure and to act normal to the shock wave line (or sweep line for subsonic cases).Therefore the drag parallel to the free stream velocity is calculated as follows: As mentioned earlier the third term in the drag equation is the induced drag.However the induced drag is already computed using Trefftz plane analysis in the VLM code.Trefftz plane analysis is a type of far field analysis method for drag computation and that is more accurate than a near field analysis.Therefore in the proposed Q3D analysis the wing parasite drag is computed by integrating the parasite drag coefficient of the 2D sections over the span and the induced drag is computed using the Trefftz plane analysis:

Structural analysis
The FEMWET tool (Elham and van Tooren 2016) is used as the core of the structural analysis in this research.However some modifications have been applied to the tool to couple it with the Q3D aerodynamic solver.In FEMWET the wingbox structure is simulated using equivalent panels.In such a way the upper skin, stringers and spars caps are modeled as the equivalent upper panel, the lower skin, stringers and spars caps are modeled as the lower equivalent panel, and the spars webs are modeled by two vertical panels.The wing aeroelastic deformation can be computed using a shell element as well as a beam element model.Dorbath et al. (2010) compared the results of a shell element and a beam element model of a wing and showed that the difference between different outputs of those two models (such as wing deflection) is about 5 %.Therefore a beam model was used in FEMWET to increase the computational speed.A finite beam element model is used to compute the wing deformation.The beam is placed at the wing box elastic axis (assumed to be the same as the sections shear centers).The positions of the shear centers are computed using the wing geometry and the thickness of the four equivalent panels.The consistent shape functions for a 3D 2-node Timoshenko beam element (Luo 2008) are used to construct the stiffness, mass and force matrices of the beam.The wing box properties such as EA, EI, GJ, etc., that are required to construct the beam stiffness matrix, are computed at each node based on the geometry, material and structural properties of the real wing box (not the beam model).For more detailed information about the finite element analysis see Elham and van Tooren (2016).As soon as the stiffness and force matrices are constructed the displacement vector, U, can be computed by solving the following equation: Using the displacement vector U, the stress distribution in the wingbox structure can be computed.Using the stress distribution both the failure criteria due to material yield and the failure criteria due to structural buckling are calculated for each element.The upper and lower equivalent panels are in fact stiffened panels, therefore the stiffened panel efficiency method (Niu 1997) is used to compute the buckling load for those panels.For the spars webs, the shear buckling is used as another failure criteria.The method presented in Niu (1997) is used to determine the shear buckling load as a function of the wingbox geometry and the material properties of the spars webs.Details of calculation of the material allowables and the buckling stresses for wing box panels can be found in the previous publication of the authors (Elham and van Tooren 2014c).
In order to compute the wing total weight an empirical equation is used.Based on the following equation from Kennedy and Martins (2014) the total wing weight is computed as a function of the optimum wingbox weight (computed based on finite element analysis) and the wing area: The factor 1.5 in (18) counts the weights that are not modeled in FEM.The second term represents the secondary weights such as leading edge, trailing edge, flaps, slats etc.
In that equation the wing area is in square meters and the wing weight is in kilogram.

Solving the coupled system
The Q3D aerodynamic solver is integrated with the FEMWET for an aerostructural analysis and optimization.Figure 3 shows an example of VLM and finite element mesh.The wing drag is not considered for computing the wing structural deformation, since its order of magnitude is negligible in comparison with the wing lift and pitching moment and the wing stiffness I yy is several orders higher than I zz .However the wing structural deformation is considered for wing drag computation.The aileron effectiveness is a very important constraint in a wing aerostructural optimization.The aileron effectiveness is defined as the ratio of elastic to rigid roll moment of the wing due to an aileron deflection: The roll moment due to an aileron deflection (L δ = dL/dδ) of an elastic wing can be computed by coupling the aerodynamic solver to the structural solver, as explained further in this section.The general aerostructural system has the following four governing equations: Fig. 3 Example for an aerodynamic and structure mesh The first and the second equations are the governing equations of the VLM and the finite element method respectively.The third equation indicates that the total lift should be equal to the design weight multiplied by the design load factor.The fourth equation is related to the Q3D analysis.It indicates that the sectional lift of 2D airfoils calculated by the airfoil analysis tool (MSES in this case) should be the same as the lift computed using the VLM (corrected for the sweep, see ( 7)).This equation needs to be satisfied to make sure that the lift distribution calculated using the strip theory is the same as the lift distribution calculated using VLM.In fact the values of downwash angles are determined based on this governing equation.C l 2d in ( 23) is determined as follows (see Figs. 1 and 4): In order to solve the coupled system the Newton method for iteration is used.Using the Newton method, the updates on the state variables (for fixed X) are found using the following equation: It should be noted that when the wing drag is not required, for example for computing stress distribution in  wingbox during a pull up maneuver, the system of governing equations can be reduced to three equations.In such cases ( 23) and all the related terms in ( 25) can be excluded.

Sensitivity analysis
In order to perform Newton iteration the partial derivatives of the governing equations with respect to the state variables (matrix J in ( 25)) are required.In addition to that when a gradient based optimization algorithm is used for wing aerostructural optimization, the sensitivities of any function of interest, for example wing aerodynamic drag or structural failure criteria, with respect to the design variables, for example wing outer shape or thickness of the equivalent panels, are needed.The presented tool computes all the required derivatives by a combined use of coupled-adjoint method, Automatic Differentiation (AD) and chain rule of differentiation.The method implemented for sensitivity analysis is explained in this section.
The partial derivatives of the governing equations R 1 to R 4 with respect to the state variables are summarized in Table 1.Starting from the first row, the partial derivative of R 1 with respect to Γ is simply the AI C matrix.To compute the partial derivatives of R 1 with respect to U and α the partial derivatives of AI C and RH S with respect to U and α are required.They are computed using AD.The Matlab AD toolbox Intlab (Rump 1999) is used for that purpose.In the second row, ∂F ∂Γ is computed using AD.The partial derivative of R 2 with respect to U is simply the stiffness matrix K.The partial derivative of lift (L) with respect to Γ in the third row is also computed using AD.
Computing the partial derivatives in the forth row of Table 1 is more challenging.The partial derivative of R 4 with respect to Γ is equal to the partial derivative of the lift coefficient perpendicular to the sweep line (C l ⊥ ) computed by VLM with respect to Γ .This term can be computed using AD easily.However in order to compute the partial derivatives of C l 2d with respect to U , α and α i the adjoint sensitivity analysis and the chain rule of differentiation are needed in addition to AD.From ( 24) one can observe that C l 2d is a function of C l eff , C d p eff , C d f eff as well as α i and θ , which is a component of U .The airfoil effective lift, pressure and friction drag are computed using MSES software as functions of airfoil geometry, effective angle of attack, Mach number and Reynolds number, see (11).MSES computes the sensitivity of the outputs with respect to the inputs using the adjoint method.Therefore using the chain rule of  The derivatives of effective lift and drag with respect to effective angle of attack, Mach number and Reynolds number are computed by MSES.The derivatives of effective angle of attack, Mach number and Reynolds number with respect to α can be computed using AD or analytically based on the method presented in Section 2. The same approach can be used to compute the sensitivity of C l 2d with respect to U and α i .
As mentioned before, in order to facilitate a gradient based optimization, the coupled adjoint sensitivity analysis method (Kenway et al. 2014) is implemented in the tool.Using that method the total derivative of a function of interest I with respect to a design variable x is computed as:  where λ = [λ 1 λ 2 λ 3 λ 4 ] T is the Adjoint vector and computed using the following equation: From ( 28) one can observe that the matrix of partial derivatives of the residuals with respect to the state variables is the same as the matrix J required for Newton iteration.
Therefore as soon as that matrix is generated during the Newton iteration, it can be used also to compute the adjoint vector.However in addition to the matrix J , the partial derivatives of the function of interest I , with respect to the state variables, the partial derivatives of residuals (R 1 to R 4 ) with respect to the design variable x and the partial derivative of I with respect to x are required to compute the total derivative of I with respect to x.All those partial derivatives are computed by a combined use of analytical methods and AD.
As mentioned earlier, in a wing aerostructural optimization the aileron effectiveness is an important constraint.If η a is defined as a constraint, its derivatives with respect to the design variables are required.In (19) the term L δ is the derivative of the rolling moment L, with respect to the aileron deflection δ.Therefore the derivative of the aileron effectiveness with respect to any design variable x, includes the second derivative of L. Although computing the second derivative using adjoint method is possible, it is computationally very expensive.Therefore a semi-analytical method is used to compute the derivative of L δ with respect to the design variables.In such an approach dL δ /dx at point (x 0 , δ 0 )is defined as follows: In ( 29) the first derivative of L with respect to x is calculated twice using coupled adjoint method, one for the aileron deflection δ 0 and one for δ 0 + Δδ.Although using this approach the coupled system needs to be solved one more time, it is still computationally more efficient than using a fully analytical approach for second derivative computation.The accuracy of dL δ /dx is strongly affected by the value of Δδ.The optimum value of Δδ, that minimizes the error of derivatives, was determined to be equal to 10 −6 .

Aircraft performance analysis
The aircraft mission fuel weight is computed using the method presented by Roskam (1986).Using that method the required fuel for the cruise is calculated using the Bréguet range equation, while some statistical factors are used to estimate the fuel weight of the other segments of the flight mission see Table 2.Each fuel weight fraction M ff i indicates the ratio of the total aircraft weight at the end of the flight segment to the total aircraft weight at the beginning of the segment.The total fuel weight fraction (M ff ) indicates the consumed fuel as a ratio of the total aircraft weight at the end of the flight mission to the total aircraft weight at the beginning: Hence, the fuel weight can be determined including a 5 % of the total fuel weight as reserve fuel using the following equation: To compute the aircraft lift over drag ration the aircraft total drag is assumed to be the sum of the aircraft wing drag and the drag of the rest of aircraft.The wing drag is computed as a function of design variables, while the drag of the rest of aircraft is assumed to be constant.The aircraft range, cruise Mach number and altitude and the engine parameters are assumed to be constant, see Section 8.The aircraft  MTOW assumed to be equal to the aircraft fuel weight, aircraft wing weight and the weight of the rest of aircraft.
The third component is also assumed to be constant.

Validation
The proposed tool has been validated in three different area: the accuracy of the tool in estimation of the wing drag, the accuracy of the tool in estimation of the wing deformation and verifying the coupled adjoint sensitivity analysis method.For validating the accuracy of the Q3D method for drag prediction, a higher fidelity CFD tool, named MATRICS-V code (van der Wees et al. 1993), is used.The MATRICS-V flow solver is based on fully conservative full potential outer flow in quasi-simultaneous interaction with an integral boundary layer method on the wing.The code uses a far field analysis method for drag prediction in transonic regime (van der Vooren and van der Wees 1991).The MATRICS-V tool was developed by NLR and has been validated using wind tunnel test as well as the flight test results for Fokker 100 aircraft, see Figs. 5 and 6.Therefore in order to validate the Q3D solver, different drag components of the Fokker 100 wing drag in cruise conditions (1g loaded wing in Mach number of 0.75) are computed by both the Q3D solver and the MATRICS-V codes.The results are shown in Fig. 7 and summarized in Table 3.The results shows a high accuracy of Q3D solver for drag prediction.
In order to validate the accuracy of the tool for computing the wing stiffness and deformation, the wing twist of the A320 aircraft under 1g load is used.Reference (Obert 2009) presents the actual wing jig twist and the wing twist under 1g load for A320-200 aircraft.In order to predict the wing twist of A320 using FEMWET, an aeroelastic   18) and the wingbox weight resulting from the optimization, the total wing weight is computed equal to 8791kg.Comparing to the actual wing weight of A320-200 (Obert 2009), which is equal to 8801kg, the error of weight estimation is -0.12 %.Of course one case is not enough to validate a tool.This case can be counted more as verification than validation.Figure 8 shows the A320 wing twist distribution under 1g load computed by FEMWET (for the optimum wing structure resulted from the optimization) and the actual twist distribution.The maximum error in wing structural deformation prediction is 8.5 % at wing tip.
Eventually in order to verify the sensitivity analysis method implemented in the proposed tool, the derivatives of different functions of interest with respect to different design variables are computed using both the coupled adjoint method and finite differencing.The results are shown in Table 4.

Test case application
As a test case, aerostructural optimization of an A320 like aircraft wing is considered.The geometry of the wing is shown in Fig. 9.The characteristics of the test case aircraft  5.In order to initialize the wingbox structure, an aeroelastic optimization was performed to the thickness distribution of the four equivalent panels from wing root to wing tip.The optimization is formulated in a way to minimize the wing structural weight subject to constraints on wing failure criteria (material tensile, compressive and structural buckling), structural fatigue and aileron effectiveness.The method suggested by Hurlimann ( 2010) is used to simulate the effect of fatigue on the wingbox structural weight.Using that method the stress in the wing box lower panel is limited to 42 % of the maximum allowable stress of the material in a 1.3g gust load case.As mentioned earlier aileron effectiveness is an important constraint in wing aerostruictural design.Elham and van Tooren (2016) showed that the wing structural weight increases quadraically by increasing the aileron effectiveness.On the other hand reducing the aileron effectiveness may results in the aircraft being unable to satisfy the rool requirements, or in worse case aileron reversal may happen.A constraint is defined to keep the aileron effectiveness in the critical roll case higher or equal to 0.52.This number is selected based on data published by BDM (1989).
Five different load cases are considered for calculation of the failure criteria.Two 2.5g pull up maneuver cases, a -1g push over maneuver, a 1.3g gust load to simulate the fatigue of the wing lower panel and a roll maneuver to calculate the aileron effectiveness.The flight condition of those mentioned load cases as well as the cruise condition are determined based on the load diagram of a similar aircraft (Dillinger et al. 2013) and listed in Table 6.This table also shows the aircraft weights used in each load case, where MTOW is the aircraft maximum take-off weight, ZFW is the aircraft zero fuel weight and W des is the aircraft design weight equal to the aircraft mid cruise weight.To compute the ultimate loads, a safety factor of 1.5 is used.The effects of the aircraft tail and the location of the center of gravity are ignored for load estimation.
The SNOPT optimization algorithm (Gill et al. 2005) is used as the optimizer.The results of the optimization (the thickness of the equivalent panels) are used for validation of the tool (see Section 7) and also as the initial wingbox structure for aerostructural optimization.
In the second step a full aerostructural optimization is formulated.The aircraft fuel weight is defined as the objective function.The aircraft design weight (see ( 22)) is defined as a function of the aircraft MTOW and aircraft fuel weight: The design vector consists of four groups of design variables.The design variables of the first group are the thicknesses of the upper equivalent panel t u , the lower equivalent panel t l , the front spar t f s and the rear spar t rs .Those thicknesses are defined at 10 spanwise positions from root to tip, so in total 40 design variables are used to optimize the wingbox structure.The second group includes the design variables defining the wing planform geometry.The wing planform geometry is parametrized using 6 design variables: root chord C r , span b, taper ration λ, leading edge sweep angle Λ, twist angle at kink kink and twist angle at tip tip .The location of the wing kink is kept constant at 37 % of the wing semi-span the same as the original wing.The third group of design variables is used to define the wing airfoils shapes.The airfoils shapes at 8 spanwise position are parameterized using Chebychev polynomials and defined as design variables.10 modes are used for parameterizing each airfoil surface, so 20 per section.As mentioned 8 sections are used for optimizing the airfoil shape, so in total 160 design variables are used for wing outer shape optimization.
The fourth group includes two surrogate variables for the aircraft fuel weight and the aircraft MTOW, that are used to avoid iterations for aeroelastic analysis.The sensitivity of ( 22) with respect to the design vector was modified accordingly.
The aerostructural optimization is subject to several constraints.The first group of the constraints includes the constraint on the structure failure.The same load cases as the initial aeroelastic optimization are used to compute  Fig. 13 Lift distribution on the initial wing and the optimized wing for minimum fuel weight the failure criteria including tensile, compressive, buckling and fatigue failure.In order to reduce the number of constraints on structural failure, the Kresselmeier-Steinhauser (KS) function (Kreisselmeier and Steinhauser 1980) is used.All the 960 original constraints on structural failure were aggregated into 22 constraints using KS function.Selection of a proper KS parameter is a challenge.Using a low value for the KS parameter results in a conservative optimization, while a large value may cause convergence difficulties for the optimization (Poon and Martins 2007), therefore a compromise is required to select the value of the KS parameter.In this research the value of the KS parameter was set to 50 as suggested by Poon and Martins (2007) as a reasonable value.
As mentioned earlier the aileron effectiveness is the ratio of L δ of the elastic wing to L δ of the rigid wing.In the initial aeroelastic optimization a constraint on the aileron effectiveness (η a ) was used.However in an aerostructural optimization this constraint is not enough to satisfy the requirements on the aircraft roll performance.The roll performance is a function of the aircraft roll moment as well as the aircraft moment of inertia.When the wing planform geometry changes both these variables change.Changing the planform geometry changes the aileron area as well as the aileron arm.Ailerons with the same η a but with different geometries may result in different L δ .Therefore to better present the effect of aerostructural optimization on the aircraft roll moment the absolute value of L δ is used as a constraint instead of η a .Computing the aircrat total moment of inertia is beyond the scope of this research and needs detailed data about the whole aircraft geometry and mass distribution.So the effect of wing geometry and mass on aircraft moment of inertia is ignored.
Another constraint is defined to keep the wing loading (aircraft MTOW divided by the wing area) lower or equal to the initial value of the wing loading.This constraint is required to make sure the aircraft can satisfy the take-off and landing requirements.
The Multidisciplinary Feasible (MDF) strategy is used to solve this MDO problem (Martins and Lambe 2013).As mentioned earlier the aircraft fuel burn is used as the objective function.Therefore aircraft performance analysis is also required to compute the value of the objective function.This will be the third discipline in addition to aerodynamics and structure.Using the MDF strategy this discipline should also be integrated with the other two.However since the focus of this research is to develop a stand-alone aerostructural analysis and optimization method two surrogate variables are used to avoid iteration between the performance analysis and the aerostructural analysis.In such a way the aerostructural analysis is decoupled from the performance analysis.The extended design structure matrix of such a problem is shown in Fig. 10.The mathematical formulation of the optimization is as follows: The history of the optimization is shown in Fig. 11.Fig. 12 shows the planform of the optimized wing.From this figure one can observe that the optimized wing has a higher aspect ratio, 13.36 for the optimized versus 9.26 for the initial wing, and lower leading edge sweep, 17.4 • for the optimized versus 27.5 • for the initial wing.The optimized wing resulted in more than 10 % lower fuel weight, more than 2 % reduction in the aircraft MTOW and more than 2 % increase in wing structural weight, see Table 7.The total drag of the optimized wing is about 28 % lower than the total drag of the initial wing.From the wing planform geometry one can find that the optimized wing has lower induced drag.The optimized wing has a higher aspect ratio and the lift distribution on the optimized wing is closer to the elliptical load distribution comparing to the initial wing, see Fig. 13.Therefore the induced drag of the optimized wing is 48 % lower than the induced drag of the initial wing.However from Fig. 13b one can observe that the outer part of the optimized wing works under larger values of lift coefficient.Besides, the optimized wing has a lower sweep angle that resulted in higher normal Mach number.Therefore although the optimizer managed to optimize the airfoil shapes in such a way to minimize the wave drag by removing shock waves (or weakening them) on the airfoils (see Figs. 14 and 15), the total pressure drag of the optimized wing is just about 3 % lower than the pressure drag of the initial wing.The  Fig. 15 The shape of the sections perpendicular to the sweep line in different wing spanwise positions friction drag coefficint is the same for both wings since a forced transition at the leading edge is used for boundary layer analysis in both wings.The drag breakdown of the wing airfoils used for Q3D analysis is shown in Table 8.The value of C d in that table is in fact the local value of parasite drag before applying the sweep correction.As explained before the optimizer tried to minimize the induced drag by increasing the wing span and pushing the load distribution on the wing toward elliptical lift distribution.Those changes are not preferred from structural point of view, and they results in heavier structure.In order to compensate those effects on the wing structural weight, the optimizer has tried to minimize the structural weight penalty in several ways.The wing sweep was reduced from 27.5 • to 17.4 • that resulted in a huge amount of structural weight reduction.The new wing is also more flexible.The initial wing tip vertical and twist deformation under 1g load are 0.57m and -1.4 degree respectively.Those values for a 2.5g load are 1.48m and -3.8 degree.However for the optimized wing the tip vertical and twist deformation under 1g load are 1.42m and -2.3 degree respectively.For a 2.5g load the wing tip is deformed 3.09m vertically and twisted by -5.6 degree, see Fig. 16.Reducing the wing stiffness also helped for more structural weight reduction.Another important factor that affects the wing structural weight is the aileron effectiveness.The aileron effectiveness is usually an active constraint and the wing weight increases quadratically with the value of it (see Elham and van Tooren (2016)).The optimized wing has a larger span that resulted in a larger aileron surface and a larger aileron arm.Therefore the same amount of L δ has been achieved with a lower aileron effectiveness.The amount of L δ for the initial and the optimize wings are the same and equal to 4.1177 × 10 6 .In fact the constraint on the wing rolling moment is an active constraint.However the aileron effectiveness (η a ) for the initial and optimized wings are 0.52 and 0.42 respectively.The lower value of the aileron effectiveness allowed the optimizer for moving toward a more flexible wing and more reduction in the structural weight was achieved.
In order to investigate whether the optimization process has exploited any fidelity issues related to the Q3D solver or not, the 1g deformed shape of the optimized wing was analyzed using the MATRICS-V code.Table 9 compares the drag reduction achieved for the optimized wing computed by the MATRICS-V code and the Q3D solver.From this table one can conclude that the optimization process has not exploited any fidelity issues related to the Q3D solver.

Conclusion
A quasi-three-dimensional method for wing drag prediction was presented.The wing total drag was decomposed into the induced drag, that was computed using a vortex lattice code by Trefftz plane analysis, and the parasite drag, that was computed by analysis of several wing sections from wing root to wing tip.The comparison of the results of the Q3D solver with a high fidelity CFD tool showed that different drag components of a transonic wing can be predicted using a VLM and a 2D airfoil analysis code, if they are connected rationally.The wingbox structure is modeled using four equivalent panels for the structural analysis.A finite beam element model was presented to predict the wing deformation under the aerodynamic and inertia loads as well as the stress distribution in the wingbox structure.The validation of that method also showed a good level of accuracy for wing weight structural estimation and prediction of the wing structural deformation.
The Q3D aerodynamic solver was connected to the finite element model for wing aerostructural analysis and optimization.The Newton method was used to solve the coupled system.In order to facilitate the optimization the gradient of the outputs, such as wing drag, weight and structural failure, with respect to the inputs, such as the wing geometry and the internal structure, were computed using analytical method.The analytical methods and the automatic differentiation were combined to compute the required gradient.
The proposed tool was used for a wing aerostructural optimization.The wing planform and airfoil geometry as well as the wingbox structure were used as design variables.The optimization is defined to minimize the aircraft fuel weight, while satisfying several constraints on the wing structural failure, wing loading and the roll requirement.The optimization result showed more than 10 % reduction in the aircraft fuel weight.That amount of fuel weight reduction was achieved mainly by reducing the wing induced drag.The optimum wing has a larger span, higher aspect ratio and lower sweep angle.The optimizer could satisfy the roll requirement by increasing the aileron surface as well as the aileron arm, while reducing the aileron effectiveness.Lower aileron effectiveness allowed a more flexible wing, that resulted in a structural weight reduction.Although due to the larger span and aspect ratio and almost an elliptical lift distribution the optimized wing is about 2 % heavier that the initial wing.
The Q3D analysis introduced in this paper can be performed on an ordinary computer in a few minutes depending on the wing geometry, the Mach number and the lift coefficient.However the same analysis using a RANS or Euler code may take a couple of hours.It should be noted that the same analysis using semi-empirical methods (for example combining a VLM with empirical profile/wave drag) can be performed in a coupled of seconds.So the computational cost of the Q3D analysis is close to the semi-empirical methods, however its accuracy is much higher than those methods and closer to the high fidelity RANS or Euler methods.
The optimization performed in this paper is a single point optimization with some important constraints neglected.Therefore a huge amount of reduction in fuel weight was achieved.A more realistic design can be obtained using a multi-point optimization including off-design points, constraints on flutter, buffet, maximum cruise speed, and geometrical constraint on wing span and aspect ratio.Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/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.

Fig. 1
Fig.1Force and angles at a typical wing spanwise 2D section

Fig. 4
Fig. 4 Wing sections perpendicular to the elastic axis before and after wing deformation

Fig. 7
Fig. 7 Comparison of computed drag by the MATRICS-V and Q3D solvers for cruise condition (1g loaded wing and M = 0.75)

Fig. 10
Fig. 10 Extended design structure matrix for wing aerostrcutural optimization

Fig. 14
Fig. 14 Pressure distribution on sections perpendicular to the sweep line in different wing spanwise positions

Fig. 16
Fig. 16 Wing jig shape and deformed shape under 2.5g pull up load

Table 1
Partial derivatives of the governing equations R 1 to R 2 with respect to the state variables

Table 2
Fuel fraction for each segment in a typical flight mission for a passenger aircraft differentiation the sensitivity of C l 2d with respect to e.g.α can be computed as follows:

Table 3
Comparison of drag prediction by the MATRICS-V and Q3D solvers

Table 4
Verification of the sensitivities computed using coupled adjoint method

Table 5
Characteristics of the test case aircraft

Table 6
Load cases for wing aeroelastic optimization

Table 7
Characteristics

Table 9
Drag count reduction of the optimized wing in cruise condition predicted by MATRICS-V and Q3D