Abstract
Aerodynamic optimization based on continuous adjoint method for a flexible wing is developed using FORTRAN 90 in the present work. Aerostructural analysis is performed on the basis of high-fidelity models with Euler equations on the aerodynamic side and a linear quadrilateral shell element model on the structure side. This shell element can deal with both thin and thick shell problems with intersections, so this shell element is suitable for the wing structural model which consists of two spars, 20 ribs, and skin. The continuous adjoint formulations based on Euler equations and unstructured mesh are derived and used in the work. Sequential quadratic programming method is adopted to search for the optimal solution using the gradients from continuous adjoint method. The flow charts of rigid and flexible optimization are presented and compared. The objective is to minimize drag coefficient meanwhile maintaining lift coefficient for a rigid and flexible wing. A comparison between the results from aerostructural analysis of rigid optimization and flexible optimization is shown here to demonstrate that it is necessary to include the effect of aeroelasticity in the optimization design of a wing.
1. Introduction
Thanks to the rapid development of computer performance and Computational Fluid Dynamics (CFD), CFD methods have been increasingly used in aerodynamic analysis of aircrafts. In the meantime, aerodynamic shape optimization based on CFD technology and optimization strategies has been more and more interesting to aircraft designers. Continuous adjoint formulation method was developed by Jameson to design aircrafts in transonic speeds, which can significantly reduce computational cost of aerodynamic shape optimization [1–5]. Subsequently, continuous adjoint method was utilized by Anderson and Venkatakrishnan on unstructured grid [6]. And recently Economon et al. in Stanford University have developed an open-source CFD code named SU2 [7–12], which is an integrated computational environment consisting of multi-physics simulation and design. The main module of SU2 is aerodynamic shape optimization on the basis of continuous adjoint method and Reynolds-averaged Navier-Stokes (RANS) equations. Now SU2 has also developed aerodynamic shape optimization based on discrete adjoint method [13, 14]. Discrete adjoint method was derived from the discrete flow equations by Baysal et al. [15, 16]. Subsequently discrete adjoint method has been developed rapidly [17–22]. Discrete adjoint equations which come from discrete flow equations can provide a little more accurate gradients than continuous adjoint equations, but discrete adjoint equations are more difficult to be derived and implemented numerically. Continuous adjoint equations are directly derived from the governing equations, which can be solved with similar methods of the governing equations and implemented more easily. However, since the coupling between aerodynamics and structures is so tight that these two disciplines cannot be isolated from each other especially for flexible wings, aerodynamic shape optimization alone may lead to the performance that is not appropriate in realistic flight conditions. But successive discipline optimization may lead to suboptimal results of coupled systems, so MultiDisciplinary Optimization (MDO) is essential to aircraft design. Relatively low fidelity models have been used in relevant disciplines in conceptual design, but these models cannot accurately capture nonlinear phenomena such as wave drag. So MDO of CFD and CSM (Computational Structural Mechanics) has been more and more attractive to aircraft researchers. MultiDisciplinary design optimization based on Kriging and response surface methodology has been a UAV wing and supersonic fighter wing [23, 24]. Martins et al. [25–27] have developed systematic aerostructural optimization design codes, which use coupled high-fidelity sensitivity analysis based on discrete adjoint method, and have been applied to design full aircraft configurations with coupled adjoint method. Mavriplis et al. [28–30] have developed time-dependent aeroelastic shape optimization based on discrete adjoint method.
However, aerostructural optimization that has been performed so far mostly is grounded on discrete adjoint method. In present work, aerodynamic optimization on the basis of continuous adjoint method including the effect of aeroelasticity is developed using FORTRAN 90, which is for the purpose of the implementation of coupled aerostructural adjoint optimization in the future. The structural model of the wing is created consisting of two spars, 20 ribs, and skin. Aerostructural analysis is performed with high-fidelity models, with Euler equations on the aerodynamic side, and with a linear shell element model on the structure side. This shell element is suitable for thin and thick shell problems and is easy to implement, and in particular this shell element can deal with the shell problems with intersections which are common in the structural model of aircrafts. Thin-Plate Spline (TPS) method is applied to ensure the interpolation between aerodynamics and structures consistent and conservative [31]. sequential quadratic programming (SQP) method, a gradient-based optimizer, is used to search for the optimal point in a coupled system consisting of aerodynamics and structures.
The paper begins with a description of aerodynamic analysis in Section 2, where the flow governing equations and the solution method are described. Section 3 presents structural analysis method with the theory of the shell element and structural design of the wing. The interpolation method between aerodynamic and structural disciplines is introduced in Section 4. The derivation of the continuous adjoint method is shown in Section 5. Design strategies containing parameterization methods, mesh deformation, and optimization algorithm are presented in Section 6. In the end, Section 7 presents verification of sensitivities of ONERA M6 and the optimization results for rigid and flexible DPW-W1.
2. Aerodynamic Analysis
The governing equations of 3D steady inviscid compressible flow arewhere is conservative variable vector and is convective flux vector. They are defined aswhere is density, , , and are Cartesian velocity components, is total energy per unit mass, , , and are Cartesian components of unit normal vector, is static pressure, is total enthalpy, and is contravariant velocity, which is represented as
The finite-volume method of HLLC (Harten, Lax, van Leer, Contact) scheme is applied to solve Euler equations, and in order to reduce the computational time, an efficient implicit parallel hybrid LU-SGS (Lower-Upper Symmetric Gauss-Seidel) method [32] is used here.
3. Structural Analysis
3.1. Shell Element
In the structural analysis of the present work, a linear quadrilateral shell element [33] is applied to model the wing structure. This four-node quadrilateral shell element with 5/6 nodal degrees of freedom is constructed based on Hellinger-Reissner variational principle with independent displacements and stress resultants. The element stiffness matrix is integrated analytically with the advantage of fast stiffness computation. The shell element can get rid of locking which many shell theories may encounter and it can deal with both thin and thick shell problems. And in essence this shell element can be implemented easily for the linear isotropic aeroelasticity.
With the stationary condition, the variational formulation based on Hellinger-Reissner principle iswhere and are displacement vector and virtual displacement vector, respectively, and are independent stress resultants and virtual stress resultants, respectively, is a vector containing membrane strains, curvatures, and shear strains, is corresponding virtual vector, is the constitutive matrix, is translational displacement vector, is volume load in , and is the boundary load on a part of boundary . The membrane strains , curvatures , and shear strains are evaluated, respectively, as follows:where where is the nodal coordinate vector, is a unit vector normal to the shell midsurface, is the deviation of , and are isoparametric coordinates with the value , is the Jacobian matrix with local coordinates to isoparametric coordinates and local coordinates are evaluated with the center of the elements, and , , , and are the middle nodes of relevant edges of the elements [33]. From (5) it is shown that the shear strains are evaluated with independent interpolation functions, while the others are interpolated from displacements. These interpolation functions are adopted to fulfill the bending patch test. The interpolation of the stress resultants is defined as follows:where and are the second and third order unit matrices, respectively, the vector contains 14 parameters in total and 8 for the constant part and the others for the nonconstant part of the stress field, respectively, and are the coordinates of the center of gravity of the element, and is the parameter of the Jacobian matrix with and [33]. After a series of transformation, (4) can be written aswhere nele denotes the number of elements, is the element displacement vector, and is the element load vector which comes from the external virtual work. The matrices and are, respectively,
In view of (9), we know that the matrix only contains polynomials of and , so it can be evaluated analytically. While the matrix can only be carried out analytically in the projected flat surface [33]. As the matrix is evaluated with the flat surface, then the corresponding transformation should be adopted to get the sound stiffness matrix.
In order to analyze the shell with intersections, a second transformation should be utilized [33]. The drilling degree of freedom is not available for the nodes on intersections. So the nodes on intersections have six degrees of freedom and the other nodes have five. The intersections can usually be seen in the internal structure of the wing, such as the intersection between the spars, the ribs, and the skin. Eventually the stiffness matrix of the shell element takes the form as follows [33]:where denotes block parameters of the element stiffness matrix .
In the analysis of static aeroelasticity, structural equilibrium equation should be solved to obtain the structural deformation:where is the global stiffness matrix which is assembled with the element stiffness matrix, is the displacement vector of nodes, and is nodal force vector which is interpolated from the aerodynamic forces. Equation (11) is a classical linear equation set, and here it is solved using GMRES [34] method to get the values of the nodal displacements. The global stiffness matrix is sparse so Compressed Sparse Row (CSR) which is popular for storing general sparse matrix is used here to store it.
3.2. Structural Design
The wing designed here is the wing-alone geometry called DPW-W1 which was developed for the Third AIAA Drag Prediction Workshop [35]. Figure 1 shows the overview of DPW-W1 (dimensions in m) and the geometric quantities are located in Table 1. The structural model is not necessary to represent exactly the internal components but is sufficient to capture the features of the wing. Furthermore since the realistic and precise structural model may have too many degrees of freedom, it is appropriate to simplify the model with the identical attributes. As the structural model of the DPW-W1 wing is not delivered, we create a structural model conforming to the wing which is a typical structure for a modern wing. The model used here consists of front and rear spars, ribs, and skin. There are 20 ribs distributed evenly along the span of the wing, not including the rib at the root for wing-body attachment. The front spar is located at 10% and the rear spar at 60% of the chord from the leading edge. In order to investigate appropriate effect of aerostructural deformation, the thickness of the two spars is adjusted to make sure that the displacement of the wing tip is about 3% of the semispan. Figure 2 shows the structural model of the wing. The skin is mirrored to the left to have a clear vision of the interior structure. The structural model of the wing is constructed with 1452 shell elements in total. The material used here for the structural model is assumed to be 7000 series of aluminium alloy with the properties shown in Table 2.
4. Aerostructural Analysis
In the aerostructural analysis, the aerodynamics provides loads to the structure and in turn the deflections are transferred from the structure to the CFD mesh. In general, the wall surface of the CFD mesh and the outer surface of the CSM mesh do not conform with each other. Thus a consistent and conservative data transfer method needs to be used. Here in the present work, TPS method which can deal with 3D problems with large deformation is applied.
The relational expression of displacements between the aerodynamics and structure isAnd the relational expression of forces iswhere the matrices and are defined as follows:where and are aerodynamic displacements and forces, respectively, and are structural displacements and forces, respectively, , , and are Cartesian coordinates of CFD mesh points, and , , and are Cartesian coordinates of CSM mesh nodes.
5. Sensitivity Analysis with Continuous Adjoint Method
Most cost functions in aerodynamic shape optimization can be represented as integrals of a function of flow variables over the solid surface: where V denotes flow variables.
As the wall surface deforms, the change of the cost function comes from the variation of the shape and the alteration of the solution due to the deformation of shape. So the variation of the cost function contains two parts as follows:wherewhere denotes the change of location of points on the wall surface and stands for the variation of unit normal vector of boundary. The cost function in the work is drag coefficient and subject to that lift coefficient is constant. As a result,where is angle of attack, is angle of sideslip, , and is reference area. After a series of transformations [36], (16) turns into
Consider the governing equations of 3D steady inviscid compressible flow as follows:where is convective flux vector and , , and are Cartesian components of . In order to circumvent the computation of in (19), the adjoint variables are introduced. Firstly linearizing Euler equations, calculating the inner product with the adjoint variables, and integrating the resulting equations over the entire domain by parts, then we obtainwhere is convective flux Jacobian matrix and denotes conservative variables. As to the third part in (21) on the far-field boundary, we can getwhere is the diagonal matrix of eigenvalues and , are right and left eigenvectors, respectively. For shape optimization, the variation of characteristic variables of the flow is , so the part on the far-field boundary of (21) can be eliminated with appropriate far-field boundary conditions for the adjoint variables [1, 6]. Applying wall boundary conditions of Euler equations, the integration over wall boundary turns into
Using linearized nontranspiration wall boundary conditions , (23) giveswhere . Substituting (24) into (21) and subtracting (21) from (19), we can get the variation of the cost function as follows:
For eliminating and in (25), set the first two terms in (25) to zero, and then we can get the adjoint equations,The wall boundary conditions of the adjoint equations are as follows:In the end, the variation of the cost function is
It can be seen that the right hand side terms of (28) contain flow variables (such as gradient of pressure and normal derivative of velocity to the surface ) and adjoint variables (such as ), so the geometry sensitivities can be acquired only after the governing and adjoint equations are solved to get flow and adjoint variables.
From (26), it can be known that the Jacobian matrix of the adjoint equations is transposed of the flow convective flux Jacobian matrix, so the adjoint equations have the same form with Euler equations. Then the adjoint equations can be solved using a finite-volume type of method similar to that used for Euler equations. As the cost of solving the adjoint equations is approximately equivalent to that of solving the Euler equations and the computational cost of the gradients does not depend on the number of design variables, the gradients with respect to all design variables for each objective function or constraint can be obtained with the cost roughly twice that of solving Euler equations. So the adjoint method is quite suitable for aerodynamic shape optimization with a large number of design variables.
6. Design Strategy
A Free-Form Deformation (FFD) method is used for shape parameterization. In FFD, the design variables are the control points of the control volume and as the design variables change, then the surface mesh inside the control volume makes corresponding changes. Then a simple classical edge spring method is adopted to deform the CFD mesh to conform to the new designed configuration. The optimization algorithm used here to drive the design to the optimal point is sequential quadratic programming (SQP) method which is gradient-based and can deal with linear or nonlinear problems with or without constraints.
The flow chart of aerodynamic optimization is shown as Figure 3(a), and Figure 3(b) is the flow chart of aerodynamic optimization including the effect of aeroelasticity. Comparing with Figure 3(a), the significant difference of Figure 3(b) is that in every design cycle aerostructural analysis is included. Because of the presence of the deflection of the static aeroelasticity, there are two shapes of the wing before and after aerostructural analysis, respectively. In the design, the deflection of the optimization is performed on the former, while the flow variables, adjoint variables, and gradients are evaluated on the latter.
(a) Aerodynamic optimization
(b) Aerodynamic optimization including aeroelasticity
7. Test Cases
The airfoil of DPW-W1 is representative of a supercritical section which is usually adopted by most modern commercial aircrafts. Here the wing has been modified with a sharp trailing edge in order to be compatible with inviscid flow solution. Firstly verification of sensitivities and aerodynamic analysis and aerodynamic optimization results for ONERA M6 wing are presented in Section 7.1. Then aerodynamic optimization for rigid DPW-W1 wing is shown in Section 7.2. As the aspect ratio of this wing is relatively large and as it deforms in realistic flight conditions, then the optimization results including the effect of static aeroelasticity are presented in Section 7.3. A comparison between the results from aerostructural analysis of rigid optimization and flexible optimization is shown. In both cases, drag coefficient is minimized with lift coefficient and the wing volume being fixed.
7.1. Verification
In this section, all the verifications are implemented using ONERA M6 wing. The flow field is computed using Euler equations with Mach 0.8395 and angle of attack on unstructured mesh with 310000 points. The surface mesh is shown in Figure 4.
The verification of the inviscid solver for ONERA M6 is presented. Figure 5 shows the comparison of pressure coefficients at different sections between aerodynamic analysis and experimental results. It can be seen that the pressure coefficient distributions match each other well.
(a)
(b)
The sensitivities of drag coefficient are verified. The sensitivities of the continuous adjoint method are compared with those of the finite difference method. The design variables are the control points on the upper surface of FFD control volume (shown as Figure 4), 5 of which are used to validate the sensitivities of the objective function . The comparison of sensitivities of drag coefficient with respect to the design variables is shown in Table 3. It is seen that the sensitivities are accurate enough to be used for aerodynamic optimization.
The pressure coefficient of ONERA M6 before and after aerodynamic optimization is shown in Figure 6, respectively. Figure 7(a) is the pressure coefficient distributions at the section of 20% semispan while Figure 7(b) shows the pressure coefficient distributions at the section of 80% semispan. The optimization results of ONERA M6 are shown in Table 4. Figures 6 and 7 show that the shock wave on the upper surface has been eliminated. So the optimization results can show that the optimization method developed in the present work can be applied to aerodynamic optimization for a flexible wing.
(a) Before design
(b) After design
(a)
(b)
7.2. Aerodynamic Optimization for Rigid DPW-W1 Wing
In this case, drag coefficient is minimized with lift coefficient fixed for rigid DPW-W1 wing. The design here is performed in typical transonic flow conditions with Mach 0.76 and angle of attack . Aerodynamic analysis is implemented on unstructured mesh which has 270000 points. The surface meshes of DPW-W1 with the FFD control volume before and after the design are shown as Figure 8. There are 208 control points in total in the FFD control volume, while 91 of them are taken as design variables which are all located on the upper surface. Figure 9(a) is the original pressure coefficient contours, while Figure 9(b) is the pressure coefficient contours after optimization. Comparing the two figures, we can see that after optimization the shock wave on the upper surface is significantly weakened. Figure 10(a) is the pressure coefficient distributions at the section of 20% semispan while Figure 10(b) shows the pressure coefficient distributions at the section of 80% semispan. It is seen directly that the shock wave along the span is eliminated. The shapes of corresponding sections are shown in Figure 11, and it is clear that the incidences of the optimized sections are decreased, equivalently with decreased local angle of attack to reduce drag. Figure 12 shows the convergence history. In this case, the drag coefficient reduces from 251 counts to 242 counts, namely, a reduction of 9 counts. While the lift coefficient varies from 0.6211 to 0.6210 which almost keeps constant. The optimization results of rigid wing are presented in Table 5.
(a) Before design
(b) After design
(a) Before design
(b) After design
(a)
(b)
(a)
(b)
7.3. Aerodynamic Optimization for DPW-W1 Wing including Aerostructural Deflections
As shown in Figure 13, the aeroelasticity has apparent influences on the performance of the wing. Figure 13(a) shows the comparison of the rigid shape (black) and the deflected one (red), and in Figure 13(b) the section shape at 80% semispan is presented. It is seen that the wing has a deflection of negative twist angle. Due to the aerostructural deflection, drag and lift are affected by about 10% and 6%, respectively, with respect to the rigid wing. So it is necessary to include the effect of aeroelasticity during the optimization design.
(a) Deflection of the wing
(b) Section shape at
Here DPW-W1 wing is optimized to minimize drag coefficient with constant lift coefficient including aerostructural deflections. The flow conditions are the same as the rigid optimization case. As seen in Figures 14 and 15, after the optimization design of the flexible wing, the shock wave on the upper wing surface is eliminated. Since the effect of aeroelasticity is introduced, it is obvious that the pressure coefficient contours before optimization in Figure 14 are significantly different from those in Figure 9, as well as the pressure coefficient contours after optimization. The comparison of different section shapes before and after optimization is shown in Figure 16 and likewise the incidence of the wing tip is reduced after design. As shown in Figure 16(b) because of aeroelasticity the section at 80% semispan is deflected upward apparently. Convergence history for flexible wing design is shown in Figure 17. In the aerostructural analysis of the original DPW-W1 wing, the resulting values of drag coefficient and lift coefficient are 223.7 counts and 0.5815, while after the optimization introducing the effect of aeroelasticity, drag coefficient decreases to 218.6 counts with the reduction of 5.1 counts and lift coefficient is 0.5807 which nearly keeps constant. The optimization results of flexible wing are presented in Table 6.
(a) Before design
(b) After design
(a)
(b)
(a)
(b)
A comparison of wing pressure distributions resulting from aerostructural analysis of rigid optimization and flexible optimization is shown in Figure 18, as well as the comparison of pressure coefficient distribution of the section at 80% semispan. In this comparison it is seen that the results from the two optimization cases are not the same, especially at the sections near the wing tip where the aerostructural deflections are large. We can see that the wave on the upper surface of aerostructural analysis of rigid optimization is weakened less completely than that of flexible optimization. So it is proved that it is necessary to include the effect of aeroelasticity in the wing optimization design. Figure 19 shows the comparison of the corresponding resulting section shapes of the two optimization cases.
(a)
(b)
(a)
(b)
8. Conclusions
The main feature of this work is to develop an optimization strategy for a flexible wing based on the continuous adjoint method. The aerostructural analysis is implemented based on high-fidelity models, namely, the Euler equations on unstructured mesh and a linear quadrilateral shell element model which is suitable for both thin and thick shells but is relatively easily implemented. Sequential quadratic programming method is adopted to search for the optimal solution using the gradients from continuous adjoint method. The flow charts of rigid and flexible optimization are presented and compared. The continuous adjoint method is used to design the rigid wing and flexible wing and these two optimization cases are compared with each other to demonstrate that it is necessary to consider the effect of aeroelasticity in the wing optimization design.
The future work is to implement the continuous adjoint methodology based on Navier-Stokes equations and to develop the optimization design system on the basis of aerostructural continuous adjoint method.
Nomenclature
| Convective flux Jacobian matrix | |
| , , , | Middle nodes of relevant edges of the elements |
| Constitutive matrix | |
| Unit vector normal to the shell midsurface | |
| Aerodynamic displacements | |
| Structural displacements | |
| Total energy per unit mass | |
| Convective fluxes vector | |
| Nodal force vector | |
| Aerodynamic forces | |
| Structural forces | |
| , , | Cartesian components of convective fluxes vector |
| Total enthalpy | |
| Coordinate transformation Jacobian matrix | |
| Objective or constraint function | |
| Global stiffness matrix | |
| Element stiffness matrix | |
| , , | Cartesian components of unit normal vector |
| Static pressure | |
| Element load vector | |
| Displacement vector of nodes | |
| Solid wall domain boundary | |
| Reference area | |
| Time | |
| Right eigenvectors | |
| Left eigenvectors | |
| Conservative variables vector | |
| Nodal translational displacement vector | |
| , , | Cartesian velocity components |
| Flow variables | |
| Structural nodal displacement vector | |
| Element displacement vector | |
| Contravariant velocity | |
| Nodal coordinate vector | |
| , , | Cartesian coordinates of CFD mesh points |
| , , | Cartesian coordinates of CSM mesh nodes |
| Diagonal matrix of eigenvalues | |
| Fluid density | |
| Stress resultants | |
| Vector containing membrane strains, curvatures, and shear strains | |
| Adjoint variable vector | |
| Angle of attack | |
| Angle of sideslip | |
| Boundary load vector on a part of boundary | |
| Flow field computational domain | |
| Structural nodal virtual displacement vector | |
| Virtual stress resultants | |
| Virtual vector containing membrane strains, curvatures, and shear strains | |
| Volume load vector in | |
| Changes of location of points on the wall surface | |
| Deviation of | |
| Structural boundary | |
| Structural computational domain | |
| Variation of unit normal vector of boundary | |
| , | Coordinates of the center of gravity of the element |
| , | Isoparametric coordinates. |
Competing Interests
The authors declare no competing interests regarding the publication of this article.
Acknowledgments
This work is supported by a project funded by the Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions.