Adjoint Optimization of a Wing Using the CSRT
Method
Michiel H. Straathof
∗, Michel J.L. van Tooren
† Delft University of Technology, Delft, 2629 HS, NetherlandsThis paper will demonstrate the potential of the Class-Shape-Refinement-Transformation (CSRT) method for aerodynamically optimizing three-dimensional surfaces. The CSRT method was coupled to an in-house Euler solver and this combination was used in an optimization framework to optimize the ONERA M6 wing in transonic conditions. The gradients of the flow variables with respect to the design parameters were computed using an adjoint solver integrated in the Euler code. The optimization was performed by a trust region reflective algorithm. A two-step approach was used to optimize the wing. First, a “general” optimization was done using the Bernstein coefficients of the shape function. Second, a “regional” refinement was performed using the B-spline coefficients of the refine-ment function. It was shown that using this strategy a considerable improverefine-ment of the lift-to-drag ratio of 22 % could be achieved. This result proves that an adjoint optimization using the CSRT method is an effective way to improve transonic wing performance.
Nomenclature
B Bernstein basis function
c Chord, m C Class function CD Drag coefficient CL Lift coefficient L/D Lift-to-drag ratio M Mach number
n Number of B-spline control points + 1
N B-spline basis function
N 1 First exponent of class function
N 2 Second exponent of class function
p Order of B-spline Bezier curves + 1
p′ Order of Bernstein polynomials ¯
Pi Vector of Bernstein (curve) coefficients
¯
Pij Matrix of Bernstein (surface) coefficients
¯
Pi B-spline (curve) control polygon
¯
Pij B-spline (surface) control polyhedron R Refinement function
S Shape function
u, w Parametric variables x Chordwise position, m
y Spanwise position, m
z Distance from airfoil centerline, m
α Angle of attack, deg
ζ Relative chordwise position
∗Ph.D. Researcher, Faculty of Aerospace Engineering, Kluyverweg 1, 2629 HS, Delft, Member, e-mail:
M.H.Straathof@tudelft.nl, phone: +31 15 2789008.
†Professor, Faculty of Aerospace Engineering, Kluyverweg 1, 2629 HS, Delft, AIAA Senior Member, e-mail:
M.J.L.vanTooren@tudelft.nl, phone: +31 15 2784794. 29th AIAA Applied Aerodynamics Conference
27 - 30 June 2011, Honolulu, Hawaii
AIAA 2011-3804
η Relative spanwise position
ψ Relative distance from airfoil centerline
I.
Introduction
“Shape parameterization” or “geometric representation” is an essential part of any shape optimization problem. It allows a shape to be expressed into a finite number of variables, preferably as little as possible to minimize computation cost. According to Kulfan1 a good parameterization technique should possess certain features, like smoothness, mathematical efficiency, appropriate design space, and ability to handle constraints. She presented a method that fulfills all these criteria: the Class-Shape-Transformation (CST) method, which she later expanded to more complex and three-dimensional shapes.2 This method combined a class function, representing a specific class of shapes, and a shape function, which defined the deviation from the class function. The shape function was based on Bernstein polynomials and their coefficients formed the design variables. Several other parameterization techniques for use in the conceptual design phase have been proposed, for example by Rodriguez and Sturdza.3 After evaluating a number of different geometric representation techniques, Samareh4concluded that polynomial and spline representations are the best approach in terms of:
1. Ability to handle surfaces 2. Consistency across disciplines
3. Ability to handle large geometry changes 4. Availability of sensitivity derivatives
Straathof et al.5 expanded the CST method with a refinement function based on B-splines, allowing for local modifications of the shape. B-splines made it possible to locally modify the shape of a curve or surface without having to increase the order of the Bernstein polynomials of the shape function. An additional advantage of having local control of a shape is that it provides a larger range of possibilities for handling geometric constraints, which was explained by Straathof and Van Tooren.6 In the same paper, it was demonstrated that coupling the CSRT method to an Euler solver and using this in an optimization framework led to significant (wave) drag reduction on an airfoil.
This paper will show that similar results can be achieved for a three-dimensional case. The CSRT method was again used, but this time the class, shape, and refinement functions represented surfaces instead of curves. This is explained in detail in Section II. The flow solution was again calculated by the Euler solver developed by Carpentieri et al.7 The gradients of the flow variables with respect to the design parameters were computed using an adjoint solver integrated in the Euler code.8 The integration of the CSRT method into the Euler and adjoint codes, as well as a validation, are presented in Section III. In Section IV an optimization test case is presented of the ONERA M6 wing9in transonic conditions. Finally, the conclusions and recommendations are presented in Section V.
II.
Class-Shape-Refinement-Transformation method
The basis of the parameterization method discussed in this paper is the Class-Shape-Transformation method, developed by Kulfan.2 Section II.A will start by introducing this method and explaining the way it has been extended using B-splines. A small discussion will be presented in Section II.B.
II.A. Definition
The Class-Shape-Transformation (CST) method, as developed by Kulfan, combines an analytical function (the class function) with a parametric curve (the shape function) and is defined as follows in two dimensions:
ζ(ψ) = CN 2N 1(ψ)· S(ψ) (1)
with ζ = z/c, ψ = x/c, and c the chord length of the wing. CN 2N 1(ψ) and S(ψ) represent the class and shape function respectively. The class function is defined as:
CN 2N 1(ψ) = (ψ) N 1
(1− ψ)N 2 (2)
If N 1 and N 2 are chosen equal to 0.5 and 1.0 respectively, then the class function describes a basic airfoil shape. The shape function then represents the deviation from this basic airfoil. Kulfan used Bernstein polynomials as a shape function. The advantage of this is that any choice of Bernstein coefficients leads to a realistic, smooth airfoil geometry. A limitation, however, is that any change in the coefficient vector leads to a global change in shape, therefore not allowing for local modifications. The use of B-splines remedies this problem. In order to make use of the advantages of B-splines, an extra function is added, called the “refinement function”. For two dimensions, it is defined in terms of the parametric variable u as follows:
R(u) = n ∑ i=0 ¯ PiNi,p(u) (3)
Where ¯Pi describes a control polygon and Ni,p are the B-spline basis functions. The definition of the
B-spline basis functions can be found in literature, for example in Mortenson.10 For the three-dimensional case, the refinement function is defined as a B-spline surface:
R(u, w) = m ∑ i=0 n ∑ j=0 ¯
Pi,jNi,puu(u)N w
j,pw(w) (4)
The matrix ¯Pij describes a control polyhedron. The B-spline extension to the CST method will be called
the Class-Shape-Refinement-Transformation (CSRT) method and is defined as follows for curves:
ζ = CN 2N 1(ψ)· S(ψ) · R(u) (5) or: ζ = (ψ)N 1(1− ψ)N 2· p′ ∑ i=0 ¯ PiBi,p′(ψ)· n ∑ i=0 ¯ PiNi,p(u) (6)
Notice the use of the symbol p′instead of p in the shape function. This is done to distinguish between the order of the Bernstein polynomials (p′) and the order of the B-spline Bezier curves (p−1). A problem arising in multiplying the refinement function with the class and shape functions is the fact that the former is defined in terms of the parametric variable u, while the latter two are defined as a function of the Cartesian coordinate
ψ. There are two ways of solving this problem. The first way is to find the values of ψ corresponding to
the values of u and then calculate the Bernstein polynomials using these values of ψ. The second way is to compute the shape and refinement functions separately and then interpolate either of them afterwards. The first method is simple and efficient, but it removes the possibility of using different resolutions for the shape and refinement functions. For this reason it is more convenient to use the second method.
For surfaces Eq. (5) and (6) change to:
ζ = CN 2N 1(ψ)· S(ψ, η) · R(u, w) (7) or: ζ = (ψ)N 1(1− ψ)N 2· p′ψ ∑ i=0 p′η ∑ j=0 ¯ Pi,j· B ψ i,p′ψ(ψ)· B η j,p′η(η)· n ∑ i=0 m ∑ j=0 ¯
Pi,jNi,puu(u)N w
j,pw(w) (8)
with η = y/c. In order to align ψ and η with u and w, the same techniques can be applied as for two dimensions.
II.B. Discussion
Naturally, adding an extra function to the CST method requires a higher number of design variables compared to the original CST method, which contributes to the complexity of the design problem. However, it is believed that the benefits outweigh this disadvantage.
One of the main benefits of the CSRT-method is the possibility to perform a shape optimization in a two-step approach. In the first “general” step, the Bernstein coefficients of the shape function are varied. In the second “regional” step, the B-spline coefficients of the refinement function are varied to refine the shape. This technique is applied in Section IV.
III.
In-house Euler code
For the test case presented in Section IV, use was made of an Euler code developed in-house by Carpentieri
et al.7 The body of this code consists of a flow solver and an adjoint solver. These two components will be
discussed in Sections III.A and III.B respectively.
III.A. Flow solver
Euler solvers neglect viscosity, but they do include rotational flow. Consequently, they are capable of modeling shock waves and hence transonic flow. Carpentieri et al.8 developed an Euler code based on an unstructured finite-volume formulation that discretizes the Euler equations on a median-dual mesh. It was shown by Carpentieri11 that this solver is capable of accurately predicting wave drag. In three-dimensional cases, also induced drag can be computed. The code is written in Fortran and consists of 5 main modules, as can be seen in Fig. 1: a mesh reading module (READMESH), a preprocessor (PREPROC), a deformation module (DEFORM), the flow solver (FLOW), and a postprocessor (POSTPROC). An effective way of implementing the CSRT method is to create a mesh of the class function and subsequently deform the mesh using the shape and refinement functions, which were programmed into the deformation module.
READMESH PREPROC DEFORM FLOW POSTPROC
Mesh CSRT coefficients Optimization Forces Aerodynamic coefficients Convergence? no yes Euler code
Figure 1. Schematic of the flow solver
To validate the implementation of the three-dimensional CSRT method into the Euler code, the ONERA M6 wing was modeled and analyzed for two different flow conditions. Figure 2 (a) shows the pressure distribution on the ONERA M6 wing at a Mach number of 0.84 and an angle of attack of 3.06 degrees, computed using the CSRT method coupled to the in-house Euler solver. Figure 2 (b) shows the pressure distribution on the same wing in the same conditions, as calculated by Sung and Kwon.12 The pressure distributions are qualitatively very similar. Both show one shock wave along the leading edge, and another parallel to the trailing edge. The shape of the rest of the isobars also seems to correspond. There appears to be a difference in the position of the aft shock wave, but that is because both figures have been plotted in different ways. The pressure distribution in Fig. 2 (b) was plotted on the actual geometry, while the one in Fig. 2 (a) was plotted on a flat surface. This makes the shock wave in the former figure look more aft. In fact, both shock waves are located at around 68 % chord.
For a quantitative validation, the results of the ONERA M6 model were compared to the results calculated by Carpentieri11 and Viviand.13 Because both references did not document details about the meshes used, a mesh sensitivity study was performed to determine which mesh density would be sufficient. The results of this study can be found below. Figure 3 shows the lift coefficient as a function of the number of mesh elements, as well as the reference values. A logarithmic trend line is also included in the figure. The value of the lift coefficient seems to converge to the value as mentioned by Carpentieri. Figures 4 and 5 show the drag coefficient as a function of mesh size, for two different flow conditions. For the first case (M = 0.84), the
a) b)
Figure 2. Pressure distributions on the ONERA M6 wing
drag coefficient converges to a value in between those mentioned by Carpentieri and Viviand. For the second case (M = 0.92) the convergence value of the drag coefficient lies just above that found by Carpentieri. From this study it was concluded that a mesh consisting of about 200000 elements would be sufficient.
0.260 0.265 0.270 0.275 0.280 0.285 0.290 0 100000 200000 300000 400000
Number of mesh elements, [~]
L if t c o e ff ic ie n t, Cl [~ ] Carpentieri Viviand CSRT CSRT trend
Figure 3. Lift coefficient as a function of mesh size (M = 0.84, α = 3.060)
The mesh used for the validation consisted of 37448 nodes and 184683 elements. Table 1 shows the results for two different flow conditions. It also lists the reference values.
Table 1. ONERA M6 results compared to Carpentieri et al.8 and Viviand13
M∞[∼] α [deg] CL[∼] CD[∼]
CSRT Carpentieri8 Viviand13 CSRT Carpentieri8 Viviand13
0.84 3.06 0.279 0.281 0.284 0.021 0.020 0.017
0.92 0.00 0.000 0.000 0.000 0.026 0.025 0.022
III.B. Adjoint solver
An additional feature of the Euler code is the ability to calculate the flow derivatives with respect to the design variables, using an adjoint solver. In this case, two extra modules are required: the adjoint module (ADJOINT) and a module to calculate the gradients (GRADIENT). Figure 6 shows the program flow with the additional two modules included.
Before the adjoint solver can be used in an optimization test case, it has to be validated to show that it works in combination with the CSRT method. In order to do this, the gradients of the aerodynamic coefficients with respect to the design variables were calculated using both (forward) finite differences (FD) and the adjoint solver. Two validation cases were run: one two-dimensional case using the NACA 0012
0.015 0.018 0.020 0.023 0.025 0.028 0.030 0 100000 200000 300000 400000
Number of mesh elements, [~]
D rag c o ef fi c ien t, C d [~ ] Carpentieri Viviand CSRT CSRT trend
Figure 4. Drag coefficient as a function of mesh size
(M = 0.84, α = 3.060) 0.020 0.022 0.024 0.026 0.028 0.030 0 100000 200000 300000 400000
Number of mesh elements, [~]
D rag c o ef fi c ien t, C d [~ ] Carpentieri Viviand CSRT CSRT trend
Figure 5. Drag coefficient as a function of mesh size
(M = 0.92, α = 0.000)
READMESH PREPROC DEFORM FLOW POSTPROC
Mesh CSRT coefficients Optimization Forces Aerodynamic coefficients Convergence? no yes Euler code
ADJOINT GRADIENT Gradients
Figure 6. Schematic of flow solver, including the adjoint solver
airfoil at M = 0.80 and α = 1.00o and one three-dimensional case using the ONERA M6 wing at M = 0.84
and α = 3.06o. The gradients with respect to both the Bernstein and B-spline coefficients were computed.
For the two-dimensional case, a step size of 0.01 was used for both the FD and adjoint calculations. For the three-dimensional case, a step size of 0.1 was used.
For both cases, the adjoint method proved to be faster than using FD. There was however a significant difference between the two- and three-dimensional cases. For the airfoil, FD took about 2 hours to compute the 20 gradients (of the upper and lower sides), while the adjoint solver needed about 90 minutes. Even though this difference would save time during an optimization process, it is not very significant. However, looking at the three-dimensional case, the strength of the adjoint method becomes apparent. It required approximately 10 hours to compute the 200 gradients using the adjoint, versus over 70 hours for FD. This means that an optimization using 10 design iterations could be performed in approximately 4 days using the adjoint code, while it would take about a month using finite differences.
The results of the validation for the two-dimensional case are shown in Fig. 7 and 8. They show the
L/D-gradients of both the Bernstein and the B-spline coefficients for the upper side of the NACA 0012 airfoil.
The graphs show very good agreement between the FD and adjoint gradients. For the three-dimensional case, the gradients are shown in Fig. 9 to 12. The FD and adjoint gradients are presented separately as contour plots. Although it would be more accurate to show the gradient values at discrete points, like in Fig. 7 and 8, this allows for a better visual comparison. There again seems to be good agreement between the FD and adjoint results.
IV.
Test case: optimizing the ONERA M6 wing
This section will describe a three-dimensional test case showing the full potential of the CSRT method. A definition of the test case will be given in Section IV.A. The results will be presented in Section IV.B.
-150 -100 -50 0 50 100 1 2 3 4 5 6 7 8 9 10 Coefficient number, [~] L /D g ra d ie n t, [ ~ ] FD Adjoint
Figure 7. FD and adjoint gradients of L/D w.r.t. Bern-stein coefficients (upper side)
-60 -40 -20 0 20 40 1 2 3 4 5 6 7 8 9 10 Coefficient number, [~] L /D g rad ien t, [~ ] FD Adjoint
Figure 8. FD and adjoint gradients of L/D w.r.t.
B-spline coefficients (upper side)
Coefficient number, [~] (in x−direction) Coefficient number, [~] (in y−direction) 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 −7 −6 −5 −4 −3 −2 −1 0 1 2 3 L/D gradient
Figure 9. FD gradients of L/D w.r.t. Bernstein coef-ficients (upper side)
Coefficient number, [~] (in x−direction) Coefficient number, [~] (in y−direction) 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 −7 −6 −5 −4 −3 −2 −1 0 1 2 3 L/D gradient
Figure 10. Adjoint gradients of L/D w.r.t. Bernstein coefficients (upper side)
Coefficient number, [~] (in x−direction) Co eff icie nt n um b er , [~] (i n y − d ir e ct io n ) 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 −2 −1.5 −1 −0.5 0 0.5 1 L/D gradient
Figure 11. FD gradients of L/D w.r.t. B-spline coef-ficients (upper side)
Coefficient number, [~] (in x−direction) Co eff icie nt n um b er , [~] (i n y − d ir e ct io n ) 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 −2 −1.5 −1 −0.5 0 0.5 1 L/D gradient
Figure 12. Adjoint gradients of L/D w.r.t. B-spline
coefficients (upper side)
IV.A. Test case definition
The ONERA M6 wing, which was used to validate the implementation of the CSRT method into the Euler code in Section III.A, was chosen as the initial solution. The optimization strategy was the same as explained in Section II.A: first a “general” optimization was performed using the shape function, followed by a “regional” refinement using the refinement function. For both steps the angle of attack was kept constant. However, if the value of the lift coefficient resulting from the “general” optimization differed significantly from the initial value, the angle of attack was changed to give approximately the same value of CL. This
new value of the angle of attack was then used for the refinement step.
For the two-dimensional test case presented by Straathof and Van Tooren,6a SQP method was employed, using finite differencing to find the gradients of the shape variables. However, for this three-dimensional case, the sharp increase in number of shape variables and time required to solve the flow field renders this approach unfeasible. Instead, adjoint sensitivity analysis was used in combination with a trust region reflective algorithm.
IV.B. Test case results
This section will present the results of the ONERA M6 test case. The pressure distributions for each optimization step will be shown, as well as the corresponding wing shapes. Additionally, the shape and refinement functions required to reach these wing shapes will be shown. The initial solution will be presented in Section IV.B.1. The results of the “general” optimization step will be shown in Section IV.B.2. The results of the “regional” refinement step are described in Section IV.B.3. Finally, Section IV.B.4 will present the results of an optimization using only the B-spline coefficients.
IV.B.1. Initial solution: the ONERA M6 wing
Figures 13 and 14 show the pressure distribution on the ONERA M6 wing for an angle of attack of 3.06 degrees and a Mach number of 0.84. This is the same pressure distribution that was used for validating the Euler code in Section III.A. It shows supersonic regions and two strong shock waves: one along the leading edge and one parallel to the trailing edge. The shock waves seem to merge near the wing tip. For this condition, the lift and drag coefficients are equal to 0.276 and 0.0237 respectively, resulting in a lift-to-drag ratio of 11.7. The reason why the coefficients differ slightly from those presented in Section III.A is that for this test case a courser grid was used (20568 nodes, 100069 cells). The pressure distribution shown in the figures formed the basis for the optimization.
Spanwise position [m] C h o rd w is e p o si ti o n [ m] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.2 1 0.8 0.6 0.4 0.2 −1 −0.8 −0.6 −0.4 −0.2 0 Pressure coefficient [~]
Figure 13. ONERA M6 pressure contours Figure 14. ONERA M6 pressure distribution
Figure 15 shows the shape function that was used to generate the mesh of the ONERA M6 wing. The latter is shown in Fig. 16. Considering that the undeformed mesh had the shape of the class function, the shape function represents the ratio between the ONERA M6 wing and the class function.
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.1 0.15 0.2 0.25 Spanwise direction Chordwise direction S h a p e f u n ct io n 0.1 0.12 0.14 0.16 0.18 0.2
Figure 15. ONERA M6 shape function Figure 16. ONERA M6 mesh
IV.B.2. “General” optimization step
For this step, the thickness of the ONERA M6 wing was taken as a minimum thickness constraint. In other words, the Bernstein coefficients of the shape function were only allowed to increase. Figures 17 and 18 show the pressure distribution after the first optimization step. It is clear from both figures that the shock wave along the leading edge has been greatly reduced. Additionally, the pressure field on the rest of the wing has flattened considerably, resulting in a more even distribution of the lift. This results in a lift coefficient of 0.381 and a drag coefficient of 0.0276 (L/D = 13.8). In order to make a fair comparison with the initial solution, the lift coefficient should be approximately equal. This was achieved by lowering the angle of attack from 3.06 to 2.00 degrees, leading to a lift coefficient of 0.279 and a drag coefficient of 0.0200. The resulting lift-to-drag ratio of 13.9 means an increase of 19 % compared to the ONERA M6 wing at the same lift coefficient. Spanwise position [m] C h o rd w is e p o si ti o n [ m] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.2 1 0.8 0.6 0.4 0.2 −1 −0.8 −0.6 −0.4 −0.2 0 Pressure coefficient [~]
Figure 17. “General” optimization pressure contours Figure 18. “General” optimization pressure distribution
The shape function belonging to the pressure distributions of Fig. 17 and 18 is shown in Fig. 19. Comparing this shape function to the one in Fig. 15, it can be seen that the ‘leading edge’ and ‘trailing edge’ of the shape function have not changed position. This is correct, since the angle of attack is adjusted through the flow solver and not through the geometry. It is also clear from the figure that the rear of the wing has generally become thicker, especially near the wing tip. It is this thickening that led to the flattening of the pressure field apparent from Fig. 18. This also explains why the top of the shape function lies near the wing tip, since this is where the shock wave was strongest. The resulting wing mesh is shown in Fig. 20.
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.1 0.15 0.2 0.25 Spanwise direction Chordwise direction S h a p e f u n ct io n 0.1 0.12 0.14 0.16 0.18 0.2
Figure 19. “General” optimization shape function Figure 20. “General” optimization mesh
IV.B.3. “Regional” refinement step
In a first attempt, the same bounds were used as for the “general” optimization step, meaning that the B-spline coefficients were only allowed to increase. This, however, did not lead to an improvement of the lift-to-drag ratio, as can be seen in Table 2. After some consideration, it was concluded that only allowing the B-spline coefficients to increase did not leave enough freedom for the refinement function to really implement small local changes. Furthermore, Fig. 21 illustrates that it is possible for some of the B-spline coefficients to be less than 1, without the B-spline itself becoming smaller than 1 at any point. For this reason, a new constraint strategy was conceived that allowed the refinement function to do its job while still maintaining the overall thickness of the wing. This time, the bounds on the B-spline coefficients were set to [0.8, 2.0], but with the added constraint that the mean value of all the coefficients had to be greater than or equal to 1. This resulted in the pressure distributions of Fig. 22 and 23. A further reduction of the shock waves and flattening of the pressure field can be distinguished, leading to a lift coefficient of 0.281 and a drag coefficient of 0.0197. The resulting lift-to-drag ratio of 14.3 means an increase of 3 % from the first optimization step and 22 % compared to the initial solution.
0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 Control point Control polygon B−spline
Figure 21. B-spline example
Figure 24 shows the refinement function resulting from the second optimization step and Fig. 25 shows the corresponding mesh. It is clear from these figures that the change in shape is very subtle. In fact, all B-spline coefficients lay between 0.9875 and 1.0625.
IV.B.4. B-spline only optimization
As a final check, it was investigated what would happen if the B-spline optimization was performed directly on the ONERA M6 wing. The new set of constraints that was used for the refinement step was also used for this case. The result of this optimization is shown in Fig. 26 and 27. The corresponding lift and drag coefficients (at α = 2.75o) are 0.283 and 0.0224 respectively. With a lift-to-drag ratio of 12.6, only a 8 % increase in L/D has been achieved, compared to a 22 % increase for the two-step approach. More importantly though, the pressure distribution on this wing is very bumpy and would probably cause the boundary layer to separate early, rendering this solution unusable. Note that this result is very similar to
Spanwise position [m] C h o rd w is e p o si ti o n [ m] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.2 1 0.8 0.6 0.4 0.2 −1 −0.8 −0.6 −0.4 −0.2 0 Pressure coefficient [~]
Figure 22. “Regional” refinement pressure contours, with new constraints
Figure 23. “Regional” refinement pressure distribution,
with new constraints
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.9 1 1.1 1.2 1.3 Spanwise direction Chordwise direction R e f n e m e n t fu n ct io n 0.9 0.95 1 1.05 1.1 1.15
Figure 24. “Regional” refinement function, with new con-straints
Figure 25. “Regional” refinement mesh, with new
con-straints
that of the two-dimensional case presented by Straathof and Van Tooren.6 Spanwise position [m] C h o rd w is e p o si ti o n [ m] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.2 1 0.8 0.6 0.4 0.2 −1 −0.8 −0.6 −0.4 −0.2 0 Pressure coefficient [~]
Figure 26. B-spline only pressure contours Figure 27. B-spline only pressure distribution
The bumpy nature of the solution is also very apparent from the shape of the refinement function and the wing mesh, as illustrated in Fig. 28 and 29 respectively. The B-spline coefficients varied significantly more than during the refinement step, their extreme values being 0.8656 and 1.3522.
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.9 1 1.1 1.2 1.3 Spanwise direction Chordwise direction R e f n e m e n t fu n ct io n 0.9 0.95 1 1.05 1.1 1.15
Figure 28. B-spline only refinement function Figure 29. B-spline only pressure distribution
IV.B.5. Summary of results
The results of the ONERA M6 test case are summarized in Table 2.
V.
Conclusions and recommendations
This paper presented a three-dimensional application of the Class-Shape-Refinement-Transformation (CSRT) technique. As a test case, the lift-to-drag ratio of the ONERA M6 wing was optimized using an in-house Euler solver. The sensitivity analysis was performed using an adjoint code integrated in the flow solver. The implementation of the CSRT method in both the flow solver and the adjoint code was validated. A trust region reflective method was chosen as the optimization algorithm. The results of the test case verified the results of earlier work performed on the NACA 0012 airfoil. Using the Bernstein shape function, a significant improvement in lift-to-drag ratio of 19 % was achieved. An extra optimization step using the B-spline refinement function led to an additional improvement of 3 %. With this two-step approach, the lift-to-drag ratio was improved from 11.7 to 14.3. It was also shown that an optimization using only the B-spline coefficients as shape variables can give convergence problems and lead to an unrealistic solution.
Table 2. ONERA M6 test case results for similar values of CL(M = 0.84)
Test case α[deg] CL[∼] CD[∼] L/D[∼] ∆L/D[%]
ONERA M6 3.06 0.276 0.0237 11.7
-“General” optimization 2.00 0.279 0.0200 13.9 19
“Regional” refinement 1.95 0.281 0.0202 13.9 19
“Regional” refinement, new constraints 2.00 0.281 0.0197 14.3 22 B-spline only optimization 2.75 0.283 0.0224 12.6 7
The work presented in this paper proves that the CSRT method is a very intuitive and effective way of parameterizating aircraft shapes, both in two as well as in three dimensions. The method allows for a two-step approach which has the potential to significantly increase the lift-to-drag ratio of various aircraft shapes. It was also shown that using an adjoint algorithm provides the computational efficiency necessary to perform true three-dimensional shape optimization.
A recommendation for future research is to investigate the applicability of the CSRT to more diverse aircraft shapes. A first step in that direction has already been taken by optimizing a blended-wing-body aircraft. This work will be presented in a forthcoming paper.
References
1Kulfan, B.M., ““Fundamental” Parametric Geometry Representation for Aircraft Component Shapes”, Proceedings of the 11thAIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Portsmouth, Virginia, U.S.A., 6-8 September 2006
2Kulfan, B.M., “Universal Parametric Geometry Representation Method”, Journal of Aircraft, Vol. 45, No. 1, pp. 142-158, 2008
3Rodriguez, D.L. and Sturdza, P., “A Rapid Geometry Engine for Preliminary Aircraft Design”, Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, AIAA 2006-929, Reno, Nevada, U.S.A., 9-12 January 2006
4Samareh, J.A., “Survey of Shape Parameterization Techniques for High-Fidelity Multidisciplinary Shape Optimization”, AIAA Journal, Vol. 39, No. 5, pp. 877-884, 2001
5Straathof, M.H., Tooren, M.J.L. van, Voskuijl, M. and Koren, B., “Aerodynamic Shape Parameterisation and
Optimisation of Novel Configurations”, Proceedings of the RAeS Aerodynamic Shape Parameterisation and Optimisation of Novel Configurations Conference, London, U.K., 27-28 October 2008
6Straathof, M.H. and Tooren, M.J.L. van, “Extension to the Class-Shape-Transformation Method Based on B-Spline”, AIAA Journal, Vol. 49, No. 4, pp. 780-790, 2011
7Carpentieri, G., Koren, B. and Tooren, M.J.L. van, “Adjoint-Based Aerodynamic Shape Optimization on Unstructured Meshes”, Journal of Computational Physics, Vol. 224, No. 1, pp. 267-287, 2007
8Carpentieri, G., Koren, B. and Tooren, M.J.L. van, “Development of the Discrete Adjoint for a Three-Dimensional Unstructured Euler Solver”, Journal of Aircraft, Vol. 45, No. 1, pp. 237-245, 2008
9AGARD 138
10Mortenson, M.E., “Geometric Modeling”, Second Edition, John Wiley & Sons, 1997
11Carpentieri, G., “An Adjoint-Based Shape-Optimization Method for Aerodynamic Design”, Delft University of
Technology, Delft, The Netherlands, 2009
12Sung, C. and Kwon, J.H., “Efficient Aerodynamic Design Method Using a Tightly Coupled Algorithm”, AIAA Journal,
Vol. 40, No. 9, pp. 1839-1845, 2002
13Viviand, H., “Numerical Solutions of Two-Dimensional Reference Test Cases”, Test Cases for Inviscid Flow Field Methods, AGARD-AR-211, May 1985
