Supercritical wing design based on airfoil optimization and 2.75D transformation

Accepted Manuscript Supercritical wing design based on airfoil optimization and 2.75D transformation

Tong Zhao, Yufei Zhang, Haixin Chen, Yingchun Chen, Miao Zhang

PII: DOI: Reference:

S1270-9638(16)30307-8 http://dx.doi.org/10.1016/j.ast.2016.07.010 AESCTE 3722

To appear in:

Aerospace Science and Technology

Received date: Revised date: Accepted date:

18 February 2016 7 July 2016 18 July 2016

Please cite this article in press as: T. Zhao et al., Supercritical wing design based on airfoil optimization and 2.75D transformation, Aerosp. Sci. Technol. (2016), http://dx.doi.org/10.1016/j.ast.2016.07.010

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Supercritical Wing Design Based on Airfoil Optimization and 2.75D Transformation*

1 2

Tong Zhao1†, Yufei Zhang1‡, Haixin Chen1§, Yingchun Chen2**, Miao Zhang3††

3 4 5 6 7 8

1. Tsinghua University, Beijing, 100084, China 2. Commercial Aircraft Corporation of China, Ltd, Shanghai 200120, China 3. Shanghai Aircraft Design and Research Institute, Shanghai 201210, China

Abstract

9

An optimization design method for a supercritical wing is developed in this paper

10

based on a pressure distribution transformation between the airfoil and wing. This

11

transformation takes the stream-wise variation of the local sweep angle and the local

12

curvature of a tapered swept wing into account. With this method, the desired pressure

13

distribution of the wing, which is critical to its supercritical characteristics, as well as

14

the constraints and objectives, can be mapped onto a 2D airfoil. When the airfoil is

15

optimized, the desired performances and the pressure distributions can be well realized

16

on a 3D wing. Test cases based on the RAE2822 airfoil and on the common research

17

model are used to validate the CFD tool and the transformation. A supercritical wing is

18

then designed using this method. Comparisons among the results from the proposed

19

method, an optimization based on the traditional transformation and a fully 3D wing

20

optimization are conducted. These comparisons show that the proposed method, by

21

essentially performing 2D optimization, greatly improves the design efficiency.

22

Moreover, the more accurate transformation helps make the optimization more reliable

23

and engineering applicable.

24

Keywords: Optimization Design, Supercritical Wing, 2D/3D Transformation, Local

25

Sweep Angle, Curvature.

*

Supported by the National Key Basic Research Program of China (2014CB744806) and the National Natural Science Foundation of China (11372160). Ph. D candidate, School of Aerospace Engineering. ‡ Assistant professor, School of Aerospace Engineering. § Professor, School of Aerospace Engineering. Corresponding author, email: [email protected] ** Research professor, General Design and Aerodynamic Department. †† Research professor, General Design and Aerodynamic Department. †

1

Introduction

1 2

The supercritical wing design is the most important issue in the aerodynamic

3

design of modern civil jets. It is crucial to an airplane’s safety and efficiency [1-4]. A

4

well-designed wing must be evaluated using many design objectives and constraints,

5

such as lift/drag ratio, stall performance, buffet margin, fuel tank capacity, and so forth.

6

Currently, the aerodynamic design is taking advantage of the rapid development of

7

modern optimization methods [5-7].

8

The genetic algorithm [8] is a popular optimization method for engineering

9

applications. This algorithm is theoretically a global optimum searching method. It has

10

good flexibility with regard to setting design constraints and objectives. The genetic

11

algorithm also has independence from the aerodynamic equations and can avoid

12

complicated formulations. However, using fully 3D wing optimization based on the

13

genetic algorithm and Reynolds-averaged Navier-Stokes equations (RANS)

14

computation to achieve a real globally optimal solution is still difficult because of the

15

considerable computational cost.

16

The pressure distribution of a supercritical wing is decisive to its aerodynamic

17

performance. However, the patterns of a "supercritical"-type pressure distribution have

18

rarely been quantitatively defined. The suction peak, shock wave location, and aft

19

loading intensity are important elements. These elements are sensitively connected to

20

the geometry due to the transonic non-linearity, which greatly increases the design

21

difficulty. Li et al. [9] and Vavalle et al. [10] achieved a feasible supercritical pressure

22

distribution curve through optimization. Zhang et al. [11] summarized the

23

characteristics of the typical types of supercritical pressure distributions. These authors

24

then developed a method that optimizes the supercritical wing with constraints on the

25

pressure distribution and successfully applied it in the design of a civil airplane. 2

1

Considering that most of the computational cost arises from the 3D CFD

2

computation, the optimization efficiency could be greatly improved if the CFD analysis

3

could be even partially simplified into a 2D computation. Moreover, for the sake of the

4

cross flow, as well as manufacturing, airfoils along the span-wise direction should have

5

smooth variations in thickness, twist angle and geometry. Vassberg et al. [2] provided

6

the geometry of the common research model (CRM) with a well-optimized supercritical

7

wing. In this design, the airfoils and the pressure distributions exhibit very good span-

8

wise similarity. An immature 3D “optimized” wing might possess poor span-wise

9

smoothness, leading to the need for meticulous manual post-adjustments. In this sense,

10

a wing “stretched” from a 2D airfoil could help guarantee a sufficiently smooth

11

geometry and pressure distribution in the span-wise direction.

12

To obtain an optimized 3D wing by simply performing optimization on an airfoil,

13

a conversion method for the pressure distributions and the aerodynamic performances

14

between the airfoil and wing has to be established. A simple traditional model called

15

the cosine rule [1] has been widely used to transform the aerodynamic parameters from

16

an airfoil to an infinite swept wing. The core principle of this rule is that the Mach

17

number can be decomposed into the components perpendicular and parallel to the

18

sweep line. Such a conversion that establishes a bridge between a 2D airfoil and a 3D

19

swept wing is therefore called a 2.5D transformation in this paper.

20

A realistic wing is always tapered to approach an elliptic load distribution and to

21

relieve the bending moment at the wing root. For a tapered wing, different chord lines

22

correspond to different sweep angles. The tapering greatly reduces the accuracy of the

23

above 2.5D transformation, which simply uses a unified sweep angle. Lock [12]

24

proposed that the pressure coefficient would be identical only if the Mach number

25

component perpendicular to the local sweep angle is equal. Based on this assumption, 3

1

he presented a pressure transformation equation for a tapered swept wing by using the

2

local sweep angle. Van der Velden and Kroo [13] improved Lock’s equation and

3

proposed a design method for the wing pressure distribution and geometry.

4

Although Lock’s method could reflect some tendency of the taper effects on the

5

pressure variation, the aerodynamic drag integrated by the “corrected” pressure has a

6

considerable error because the drag is quite sensitive. Petruzzelli and Keane [14]

7

presented an airfoil-to-wing transformation method for the wave drag considering the

8

taper effect. Atkin and Gowree [15] further proposed the Callisto-Viscous-Garabedian-

9

Korn (CVGK) method. This method revealed that the proportions of the friction drag

10

to the pressure drag are highly different for the airfoil and wing, which would strongly

11

affect the optimization results when minimizing the total drag is set as an objective.

12

For the majority of these transformations, a nominal sweep angle is required.

13

Aircraft engineers generally tend to use the sweep angle of the 1/4 chord line. However,

14

Streit et al. [16] suggested that the sweep angle near the shock wave location is better

15

for the transonic wing.

16

3D wing optimization of a civil aircraft is prohibitively expensive, whereas the

17

traditional 2.5D method is not sufficiently accurate because of the absence of the taper

18

effect. A method with a low computational cost and high reliability is always desired

19

for engineering applications. The objective of this paper is to develop a design method

20

for supercritical wings. The proposed method is based on a 2D/3D transformation

21

considering the effects of both the sweep angle and the taper ratio. With the so-called

22

2.75D transformation, the optimization can be conducted upon the airfoil but targeted

23

on the wing’s performance. An improved selection of the nominal sweep angle for a

24

more accurate transformation will be proposed. With a reallocation of its components

25

during the drag’s transformation from the airfoil to the wing, the reliability of the multi4

1

point design can be improved. By setting constraints on an airfoil’s pressure distribution,

2

the optimization can help a wing achieve the desired pressure distribution.

3

This paper first introduces the computational fluid dynamics (CFD) and

4

optimization tools used for the research. Then, the process and formulae of the 2.75D

5

transformation are introduced and numerically validated through test cases. Finally, the

6

optimization process, including the geometry parameterization and the objective and

7

constraint settings, is presented for the wing design of a realistic wide-body civil aircraft.

8

The optimizations based on the traditional 2.5D transformation, the proposed 2.75D

9

method and the fully 3D wing optimization are compared and analyzed.

10 11 12

1. Numerical Method and Optimization Method

13

1.1 CFD Solver

14

The CFD analysis in this paper is based on an in-house multi-block structured code

15

called NSAWET [17,18], which has been validated through many test cases and applied

16

in many design activities [19-21]. Based on a cell-centered finite volume architecture,

17

the code has multiple types of numerical schemes and turbulence models for solving

18

the compressible Navier-Stokes equations. For the research in the present paper, it

19

employs the MUSCL scheme for the flux reconstruction, Roe's scheme for the spatial

20

discretization, the lower-upper symmetric Gauss-Seidel method for the time advancing,

21

and the shear stress transport (SST) model for the turbulence modeling. Because the

22

present paper will use 2D computations for the optimization and 3D computations for

23

the validation and discussion of the results, the solver is validated by both 2D and 3D

24

test cases in this section.

25 5

1 2

1.1.1 2D Validation

3

Computation on case 9 of the ONERA RAE2822 airfoil (M=0.734, Re=6.5 million,

4

and AoA=2.54°) [22] is illustrated in Fig. 1. Four C-type grids with 183× 69 ,

5

257× 97 , 363×137 and 513× 193 points in the circumferential direction and the

6

normal direction, respectively, are used to test the grid convergence. Fig. 1(a) shows

7

the grid with 257× 97 points. Normalized by the chord length, the heights of the first

8

cell in the normal direction are all set to be 2.0 × 10 −6 , which guarantees that the Y+

9

values are less than 1.

0.0164

-1 0.0162

-0.5

Cd

Cp

0.016

0

0.0158

Exp ExtraCoarse Coarse Medium Fine

0.5

1

10 11 12

0

0.2

0.4

0.6

0.8

0.0156

0.0154

1

X

(a) C-Type Grid

(b) Pressure Coefficient Results

2.5E+04

5.0E+04

7.5E+04

1.0E+05

Number of Grids

(c) Grid Convergence

Fig. 1. A Grid for the RAE2822 Airfoil and the Cp Results by Different Grids

13 14

Fig. 1(b) shows the pressure distribution results computed using the four grids.

15

These results are compared to the experimental data. The numerical results are almost

16

identical, and they are in good agreement with the experimental data. Table 1 provides

17

the drag coefficients with different grids computed for the same lift coefficient. In Fig.

18

1(c), from the extra coarse grids to the fine grids, there is a clear trend of grid

19

convergence. As shown in Table 1, the drag difference between the coarse grids to the

20

fine grids is 0.00026, which is less than 2 percent. These results demonstrate the good

21

accuracy of the solver and the adequate resolution of the grid. Throughout the remainder

22

of this paper, the coarse grid ( 257× 97 points) is used for the computation in the airfoil 6

1

optimization to reduce the computational cost.

2

Table 1. Comparison of the CD Results by Different Grids (M=0.734, CL=0.8, and Re=6.5 million) Extra Coarse

183× 69 0.01623

Grid CD

Coarse 257 × 97 0.01589

Medium 363× 137 0.01571

Fine 513×193 0.01563

3 4

1.1.2 3D Validation

5

The CRM is used to validate the NSAWET’s accuracy of 3D computation. A 30

6

million mesh is used to analyze the wing-body configuration at Re=5.0 million and

7

M=0.85. The computed polar curve is shown in Fig. 2 and is compared to the

8

experimental result from the National Transonic Facility (NTF) [23]. These two curves

9

clearly have very good agreement. The pressure distribution results are also compared

10

to the experimental data at CL=0.521 in Fig. 3. The agreements are satisfactory at all

11

four sections, which are at 13.1%, 28.3%, 50.2% and 72.7% of the half-span. At this

12

point, the CFD-predicted drag coefficient is 0.02578, whereas the experimental result

13

is 0.02591.

0.6

Cl

0.5 0.4

NSAWET-SST EXP-NTF

0.3 0.2

14 15

0.02

0.03

0.04

0.05

Cd

Fig. 2 Computed Polar Curve of the CRM (Re=5.0 million, M=0.85 and CL=0.521)

16

7

13.1%

28.3%

-1

-0.5

Cp

Cp

-0.5

0

0

NSAWET-SST Exp_NTF

0.5

NSAWET-SST Exp_NTF 0.5

0

0.2

0.4

1

X

0.6

0.8

1

0

0.2

0.4

X

0.6

50.2%

0.8

1

72.7% -0.5

Cp

Cp

-0.5

0

0

NSAWET-SST Exp_NTF

0

2 3

0.2

0.4

X

0.6

NSAWET-SST Exp_NTF

0.5

0.5 0.8

1

0

0.2

0.4

X

0.6

0.8

1

Fig. 3 Pressure Distributions of the CRM (Re=5.0 million, M=0.85 and CL=0.521)

4 5

1.2 Optimization Method

6

According to the authors’ previous work [7,11,24], the genetic algorithm (GA) is

7

used as the optimization method in the present paper because of its easy implementation,

8

simple handling of constraints and excellent multi-objective performances. One of the

9

GA’s variants, NSGA-II (non-dominated sorting genetic algorithm II) [8] is chosen,

10

which features a fast sorting procedure and a more effective elite strategy. This strategy

11

guarantees that each generation will not retrograde and ensures a faster optimum

12

searching. More details of NSGA-II can be found in reference [8].

13 14 15

2. Transformation Methods

16

2.1 2.5D Transformation between Airfoil and Infinite-Span Swept Wing 8

1

As described in Fig. 4, the 2.5D transformation or the cosine rule [1] between an

2

airfoil and an infinite-span swept wing can be written as Eq. (1), where Λ is the sweep

3

angle. Denoted with superscript “Ԣ”, C p′ , C L′ , C D′ , M ∞′ and

4

of the corresponding airfoil. C p , C L , C D , M ∞ and

5

wing.

t c

′

are the parameters

t are the parameters of the c

6

Eq. (1) is widely used for determining the lift coefficient of the airfoil for the wing

7

design. However, because of the root and tip effects of a realistic wing, the right side of

8

Eq. (1) always needs to be magnified, as C ′L = (1.1 ~ 1.2) ⋅ C L / cos 2 Λ 1/4C [25,26].

9

Although this approximation is simple to use, it has too much uncertainty on the airfoil

10

optimization.

Cp′ = Cp / cos 2 Λ CL′ = CL / cos 2 Λ 11

12 13

CD′ = CD / cos 2 Λ M ∞′ = M ∞ ⋅ cos Λ ′ t t = / cos Λ c c

Fig. 4. The Baseline Airfoil (Black) and the Corresponding Airfoil (Red) of a Swept Wing

14 15

2.2 2.75D Transformation between Airfoil and Tapered Swept Wing 9

(1)

1

As shown in Fig. 5, let x , y and z be the coordinates of the chord direction, span-

2

wise direction and normal direction, respectively. Assume that the extension lines of a

3

tapered swept wing’s leading edge and trailing edge intersect at point S. Λ  is the

4

sweep angle of the leading edge. r is a vector from point S to any point on the surface

5

of the wing. Let θ be the spatial angle formed by r and the leading edge of the wing

6

when r sweeps along the surface of the wing. r and θ form a curvilinear

7

coordinate system. Let r = r .

→

→

8 9

→

→

→

Fig. 5. Definitions of the Parameters on a Tapered Swept Wing

10 11

For a point on the surface, Lock [12] defined the local normal Mach number Mθ ,

12

which is the component normal to the local contour line of pressure or velocity.

13

Considering an ideally designed wing and using the conical flow model hypothesis, as

14

proposed in reference [16], these contour lines should be parallel to the local sweep line

15

everywhere. According to Lock’s hypothesis, an airfoil and a wing are connected if they 10

1

have the same Mθ distributions. This is the core of the 2D/3D transformation.

2

Different approaches have different methods for computing Mθ .

3 4 5 6 7

2.2.1 2.75D Local Sweep Angle Pressure Coefficient Transformation →

Neglecting the thickness of the wing, the projection of r on the x-y plane is treated as the local sweep line, with a local sweep angle Λ . With the assumption of isentropic flow, Lock’s equation [12] establishes a

8

relationship between the wing’s M θ and C p , as shown in Eq. (2).

9

1 + 12 (γ − 1) M θ2 =

10 11

1 + 12 (γ − 1) M ∞2 ⋅ cos 2 Λ

(2)

(1 + 12 γM ∞2 C p ) (γ −1) / γ

For the 2D pressure distribution on the corresponding airfoil, the M distribution should be identical to that of M θ . It has:

12

~ 1 + 12 (γ − 1) M ∞2 1 + (γ − 1) M θ′ = ~ (1 + 12 γM ∞2 C p′ ) (γ −1) / γ

13

~ where C ′p is the pressure coefficient computed on the corresponding 2D airfoil. M ∞

14

~ ~ is the far-field Mach number of the airfoil, computed as M ∞ = M ∞ ⋅ cos Λ in Eq. (1).

15

~ For a tapered swept wing, a nominal sweep angle Λ should be assigned.

16

17

18 19 20

2

1 2

(3)

From Eq. (2) and Eq. (3): ~ ~ § 1 + 12 γM ∞2 C p′ 1 + 12 (γ − 1) M ∞2 ¨ = 1 + 12 (γ − 1) M ∞2 cos 2 Λ ¨© 1 + 12 γM ∞2 C p

· ¸ ¸ ¹

( γ −1) / γ

(4)

A simple equation for the pressure distribution transformation from airfoil to wing was then deduced by Lock, as in Eq. (5).

Cp =

1 2

f −1 ~ + fC p′ ⋅ cos 2 Λ 2 γM ∞ 11

(5)

1 2

f is the correlation coefficient between the local sweep angle and the nominal

sweep angle, as in Eq. (6). γ 2 2 ~ γ −1 1 § · 1 + 2 (γ − 1) M ∞ ⋅ cos Λ ~ ¸ f ( M ∞ , Λ, Λ ) = ¨¨ 2 2 1 ¸ 1 ( γ 1 ) M cos + − ⋅ Λ ∞ 2 © ¹

3

(6)

4

Eqs. (5)~(6) form the transformation between the pressure distributions of a 2D

5

airfoil and a 3D tapered swept wing. This method considers only the local sweep angles,

6

and it eliminates the need for explicitly computing M θ . This method is denoted as

7

"2.75D_localsweep" in the following discussion.

8

The above deduction is based on the assumption of isentropic flow. However, for

9

a transonic wing, the flow may have a shock wave. This issue is discussed in this paper.

10

On both sides of the shock wave, the flow can still be treated as isentropic. However,

11

the isentropic relation for M θ in Eq. (2) and Eq. (3) needs to be solved separately

12

before and after the shock wave. The normal shock wave relations need to be used for

13

determining the total pressure loss caused by the shock wave.

14

M θ′ _ shock and Mθ _shock are defined as the normal Mach number immediately

15

before the shock wave of the airfoil and the wing, respectively. The subscripts “1” and

16

“2” denote the points before and after the shock wave, respectively. The total pressure

17

ratio is

′ p01 p for the airfoil and 01 for the wing, as shown in Eq. (7). ′ p02 p02

18

′ p 01 = k 0 ( M θ′_shock ), ′ p 02

19

where

ª γ +1 M 2 º γ −1 k 0 ( M ) = « 2γ −1 2 » ¬1 + 2 M ¼

20

At a point on the airfoil after the shock wave, if M θ′ 2 is the Mach number component

21

′ can be formulated as Eq. normal to the sweep line, then the relation of M θ′ 2 and C p2

p 01 = k 0 ( M θ_shock ) p02

(7)

γ

12

[

2γ γ +1

M 2 − γγ +−11

]

−1

γ −1

(8)

1

(9) by introducing the shock relation: ~ ~ ~ 1 + 12 (γ − 1) M ∞2 1 + 12 (γ − 1) M ∞2 1 + 12 (γ − 1) M ∞2 = = p′ p′ p′ p′ p′ ( 2 ) (γ −1) / γ ( 2 ⋅ ∞1 ) (γ −1) / γ ( 2 ⋅ 01 ) (γ −1) / γ ′ p ∞′ 2 p ∞′ 1 p ∞′ 2 p ∞′ 1 p 02 ~ 1 + 12 (γ − 1) M ∞2 = ~ ′ ) (γ −1) / γ ⋅ (k 0 ( M θ′_shock )) (γ −1) / γ (1 + 12 γM ∞2 C p2

1 + 12 (γ − 1) M θ′ 22 =

2

(9)

3

where p′2 is the local static pressure, p∞′ 2 is the static pressure with the Mach

4

number isentropically varied to

5

flow prior to the shock wave.

6

Similarly, the equation of M θ 2 and C p2 for a point after the shock wave on the

7

3D wing is written as Eq. (10):

8

1 + 12 (γ − 1) M θ22 =

9

1 + 12 (γ − 1) ⋅ M ∞2 cos 2 Λ (1 + 12 γM ∞2 C p2 ) (γ −1) / γ ⋅ ( k 0 ( M θ_shock )) (γ −1) / γ

(10)

Therefore, Eq. (4) becomes: (γ −1) / γ ~ § 1 + 12 γM ∞2 C p2 ′ · 1 + 12 (γ − 1) M ∞2 ¨ ¸ = 1 + 12 (γ − 1) M ∞2 cos 2 Λ ¨© 1 + 12 γM ∞2 C p2 ¸¹

10

11

, and p∞′ 1 is the static pressure of the incoming

§ k 0 ( M θ′_shock ) · ¸ ⋅¨ ¨ k (M ¸ ) © 0 θ_shock ¹

(γ −1) / γ

(11)

Eq. (5) becomes:

C p2 =

12

f ⋅ k −1 ~ ′ ⋅ cos 2 Λ + f ⋅ kC p2 2 1 2 γM ∞

k ( M θ′_shock , M θ_shock ) =

k 0 ( M θ′_shock )

(12)

(13)

13

where

14

Because the M θ distributions of the airfoil and the wing should be identical, namely,

15

M θ′_shock = M θ_shock , k = 1 . Eq. (12) is therefore identical to Eq. (5). The existence of a

16

shock wave does not change the formulation of the transformation.

k 0 ( M θ_shock )

17 18

2.2.2 2.75D Local Curvature Pressure Coefficient Transformation 13

1

The effects of the airfoil’s curvature are significant, particularly at locations such

2

as the leading edge. Lock [12] further investigated the flow path along the wing surface

3

and proposed a correction to the transformation based on the local curvature. For a wing, we can construct the geometry with the baseline airfoil by defining

4 5

z = g ( x ) , x ∈ [ 0,1]

6

as:

at the wing root (y=0); then, the equation of the wing surface is written

z = (1 −

7 8

where ξ = ( xy0 − x0 y ) /( y 0 − y )

9

normalized by the local chord.

10 11

12 13 14 15 16 17 18 19 20

y ) g (ξ ) y0

(14)

is the chord-wise coordinate of a surface point

→

Therefore, the vector r for a surface point (x,y) is: ­r = ( x − x ) 2 + ( y − y ) 2 + z 2 0 0 °° ® ξ 1 + [ g ′(ξ )] 2 dξ °θ = y 0 ³ 0 ( x − ξ ) 2 + y 2 + g 2 (ξ ) °¯ 0 0

(15)

→

uθ and ur are defined as the flow velocities perpendicular and parallel to r , respectively. The total velocity is: q 2 ( r , θ ) = u r2 + uθ2

(16)

As suggested by Lock [12], on a wing, with a good approximation, q will be obtained by the isentropic relations from the pressure coefficient: 1 + 12 (γ − 1) M ∞2 ⋅ (1 −

q2 ) = (1 + 12 γM ∞2 C p ) ( γ −1) / γ U ∞2

(17)

The local sonic speed can also be estimated by the isentropic relations: 2 1 clocal q2 1 = + ( γ − 1 )( 1 − ) 2 U ∞2 M ∞2 U ∞2

The local normal Mach number Mθ can be obtained by: 14

(18)

Mθ =

1 2

uθ clocal

With Eq. (17) and Eq. (18), Eq. (19) becomes 1 + 12 (γ − 1) M ∞2 ⋅ (1 − 1 + 12 (γ − 1) M θ2 =

3

4

(19)

u r2 ) U ∞2

(1 + 12 γM ∞2 C p ) (γ −1) / γ

(20)

With Eq. (3), the relation of C p and C p′ is:

~ 1 + 12 (γ − 1) M ∞2

§ 1 + 12 γM ∞2 C p′ ¨ = ¨ 1 + 1 γM 2 C u r2 2 ∞ p 2 1 1 + 2 (γ − 1) M ∞ ⋅ (1 − 2 ) © U∞

5

· ¸ ¸ ¹

(γ −1) / γ

(21)

6

In the curvilinear coordinate system of Fig. 5, let hr and hθ be the metrical

7

coefficients. The vorticity in the direction normal to the surface on the curvilinear

8

coordinate system is as shown in Eq. (22):

9 10

Ωn =

∂ ∂ ( hr u r ) = ( hθ uθ ) ∂θ ∂r

15 16

(23)

Because hr = 1 and hθ = r , Eq. (23) becomes: 1 ∂u r uθ ∂uθ = + r ∂θ r ∂r

13 14

(22)

In a well-designed wing flow, such a vorticity should be negligible, which leads to

11 12

º 1 ª ∂ ∂ « ( hr u r ) − ∂ ( hθ uθ ) » hθ hr ¬ ∂θ r ¼

(24)

Again, by the conical flow model, uθ is independent of r ; then, Eq. (24) becomes: ∂u r = uθ = q 2 − u r2 ∂θ

(25)

17

Through a numerical integration from the leading edge to the trailing edge, the ur

18

and then uθ are calculated from Eq. (25). The boundary condition is set on the leading

19

edge of the wing, where uθ = 0 , θ = 0 , and ur = −q = −U∞ ⋅ sin Λ 0 . 15

1

After obtaining u r , the Cp distribution of the wing can be obtained using Eq. (21)

2

with the pressure coefficient C′p of the airfoil. This transformation method is referred

3

to as the 2.75D local curvature method in this paper.

4

For the issue of the shock wave, the discussion for the 2.75D local sweep angle

5

method is still applicable. Because the distribution of M θ should be identical for the

6

wing and the corresponding airfoil, Eq. (21) remains unchanged with the introduction

7

of the shock wave effects.

8 9

2.2.3 Drag Coefficient Transformation

10

Drag reduction is one of the most important objectives for wing design. The

11

computation of an airfoil can provide the pressure drag and the friction drag. However,

12

the ratio between the two components needs to be adjusted when the total airfoil drag

13

is transformed into the wing drag. The induced drag cannot be calculated from the 2D

14

airfoil and transformed onto the 3D wing. However, because the present optimization

15

does not change the planform and the twist angle distribution along the span-wise

16

direction of the wing and because the span-wise load distribution of the wing is

17

controlled by the pressure distribution constraints, the change in the induced drag of the

18

wing during the optimization is negligible.

19

In previous “airfoil optimization for 3D wing” studies [27-29], many researchers

20

chose to minimize the total drag of the corresponding airfoil as the optimization

21

objective, which implies that the weights of the pressure drag and the friction drag are

22

kept unchanged when being transformed from the 2D airfoil to the 3D wing. The drag

23

components might have a negligible influence on the result of a single-point

24

optimization. However, it may have a considerable effect on a multi-point optimization.

25

For example, when the free-stream Mach number is increased from the cruise condition 16

1

to the drag divergence condition, the friction drag changes only little. It is the pressure

2

drag that dominates the increase in the total drag. Hence, for different Mach numbers

3

or lift coefficients, if the drag transformation holds a constant ratio of the

4

friction/pressure drag, then the variation of the total drag will not be correct and the

5

optimization will be misled. Because multiple design points with different Mach

6

numbers and CL need to be set for the optimization for the supercritical wing design,

7

the transformation of the drag should be able to correctly reflect the change in the drag

8

components for these design points.

9

As shown in the introduction, Atkin and Gowree [15] introduced a relation of drag

10

~ ~ components, called the CVGK method, which proposed cos 0.2 Λ and cos3 Λ as the

11

correction factors for friction drag and pressure drag, respectively. The total drag of the

12

wing can be written as

13

~ ~ ′ ⋅ cos 3 Λ + C Df ′ ⋅ cos 0.2 Λ C D = C Dp

14

′ are the pressure drag and the friction drag of the corresponding ′ and CDf where C Dp

15

airfoil, respectively.

(26)

16

Such a reallocation leads to a dramatic reduction in the weights of the form drag

17

and the wave drag computed for a wing. In the present 2.75D methods, this drag

18

transformation is adopted.

19 20

2.2.4 Validation by a Simple Wing Model.

21

A simple tapered swept wing is constructed as the test case, as shown in Fig. 6(a).

22

The taper ratio is 0.33, and the aspect ratio is 9.38. The sweep angles of the leading

23

edge, the 1/4 chord line and the 1/2 chord line are 30 degrees, 27.6 degrees and 25.2 17

1

degrees, respectively. The RAE 2822 airfoil is used as the corresponding 2D airfoil.

2

Because of the cross flow effect of the swept wing, the lift can hardly be maintained

3

constant in a span-wise manner. A linear twisting of -4 degrees is found by trial-and-

4

error and employed to obtain similar lift coefficients along the span-wise direction to

5

satisfy the conical wing hypothesis. Extrapolation boundary conditions are adopted in

6

the span-wise direction to reduce the tip and root effects. The wing is computed by

7

NSAWET to obtain a 3D RANS result. The flow conditions are M=0.812, Re=7.3

8

million, and CL=0.655. Fig. 6 (a) shows the contours of the pressure coefficient. The

9

straight contour lines indicate that the tip and root effects are essentially negligible. The

10

pressure distributions at 35%, 50% and 65% of the half-span are compared in Fig. 6(b),

11

which shows that the pressure distributions are quite close in the span-wise direction.

12

Hence, the pressure distribution at 50% of the half-span is chosen to validate the

13

transformation methods.

14

-1

Cp

-0.5

0

Cp_35% Cp_50% Cp_65%

0.5

0

15 16 17

0.2

0.4

0.6

0.8

1

X

(a) pressure contour

(b) Cp distributions

Fig. 6. Pressure Contour (Left) and Cp Distributions (Right) for Tapered Swept Wing

18 19

A nominal sweep angle is required for the 2.75D methods. Normally, the sweep

20

angle of the 1/4 chord line is often chosen. In this paper, three sweep angles of the

21

leading edge, the 1/4 chord line and the half chord line are tested in the 2.75D local 18

1

curvature method. The transformed results from the RAE2822 airfoil computation

2

(M=0.734 and CL=0.82) for different sweep angles are shown in Fig. 7. As shown, the

3

transformation results with the leading edge sweep angle and the 1/4 chord sweep angle

4

have large deviations from the 3D results on the upper surface. In contrast, the result

5

based on the 1/2 chord sweep angle matches the 3D result well. The selection of the

6

nominal sweep angle has a significant influence on the transformation. Streit et al. [16]

7

showed that for the transonic cases, the sweep angle of the shock wave location line

8

performs better than that of the 1/4 chord. The shock wave location is approximately

9

50~60% chord in the present case. The 1/2 chord line is to some extent consistent with

10

Streit's selection. Most modern civil aircraft employ supercritical wings. In fact, their

11

locations of shock waves are normally designed between 40% and 60% chord.

12

Apparently, the sweep angle of the 1/2 chord line is considerably easier for engineering

13

implementation. It is therefore selected as the nominal sweep angle in the present paper.

-1 -0.8 -0.6 -0.4

Cp

-0.2 0 0.2 2.75D_50%CordLine 2.75D_LeadingEdge 2.75D_25%CordLine 3D_Result

0.4 0.6 0.8 0

14 15

0.2

0.4

0.6

0.8

1

X

Fig. 7. Comparison of 2.75D_Curvature Results Using Different Nominal Sweep Angles

16 17

The 2.75D local sweep angle method and the 2.75D local curvature method are

18

compared in Fig. 8, in which the “2.5D” is the simple cosine rule and “3D” is the 3D

19

CFD result of the wing. The 2.5D method fails to predict the suction peak of the wing,

20

as shown in Fig. 8(b). Both of the 2.75D methods have better agreement with the 3D 19

1

results on both the suction peak in Fig. 8(b) and the front-loading zone in Fig. 8(c). In

2

Fig. 8(a), the blue dashed line quantifies the pressure coefficient deviation between the

3

2.5D and the 2.75D local curvature results. The improvement is small but solid. In Fig.

4

8(b), the 2.75D local curvature method is more accurate than the 2.75D local sweep

5

angle method at the leading edge. With the error being accumulated in a stream-wise

6

manner, the advantage of the 2.75D local curvature method becomes more evident on

7

both the upper and lower surfaces of the trailing edge, as shown in Fig. 8(d).

-1

-1 0.04

-0.8

-0.95

-0.6 0.02

Cp

-0.2 0

0

-0.9

ΔCp Cp

-0.4

-0.85

0.2 0.4 0.6 0.8 0

-0.8

-0.02

2.5D_Result 2.75D_LocalSweep 2.5D_Curvature 3D_Result 2.75DCur_Delta

-0.04

0.5

2.5D_Result 2.75D_LocalSweep 2.75D_Curvature 3D_Result

-0.75

1

0

0.05

X

8 9

0.1

X

(a) Overall

(b) Suction Peak

-0.2 0.2

Cp

Cp

0

0.4

0.2 2.5D_Result 2.75D_LocalSweep 2.75D_Curvature 3D_Result

0

10 11 12

0.05

0.1

2.5D_Result 2.75D_LocalSweep 2.75D_Curvature 3D_Result

0.4 0.15

0.65

X

0.7

0.75

0.8

0.85

0.9

0.95

1

X

(c) Front Part of Lower Surface

(d) Trailing Edge

Fig. 8. Comparison of Cp among the Transformations for Tapered Swept Wing

13

Although the 2D/3D transformation is based on inviscid flow assumptions, it is

14

well applied for the above transformation between the 2D and 3D viscous results. For

15

the high Reynolds number viscous flow, which is mainly considered in this paper, the

16

boundary layer is relatively thin. The existence of the non-isentropic boundary layer’s 20

1

effect on the transformation is negligible. Therefore, the transformation’s accuracy is

2

satisfactory.

3

Fig. 9 presents a comparison of the drag coefficients obtained by the transformation

4

from the airfoil results by Eq. (1) (2.5D) and Eq. (26) (2.75D) and by the 3D CFD

5

results. The drag increments with the Mach number and the lift coefficient are plotted.

6

The vertical axis is for the lift coefficient in Fig. 9(a) and for the Mach number in Fig.

7

9(b). The horizontal axis is for the drag increment ǻC D . In Fig. 9(a), the Mach number

8

is fixed at 0.812, and the wing’s lift coefficient is increased from 0.49 to 0.655; the drag

9

coefficient of the CL=0.49 case is set as the baseline value for the calculation of ǻC D .

10

In Fig. 9(b), the lift coefficient is fixed at 0.655, and the Mach number is from 0.775 to

11

0.80. The drag coefficient of the M=0.775 case is set as the baseline. It shows that the

12

discrepancy of the 2.5D method increases with both the Mach number and the lift

13

coefficient. However, the 2.75D results are quite close to the 3D results. These

14

comparisons demonstrate that Atkin's relation can well capture the drag’s variation

15

tendencies with the Mach number and the lift coefficient. Such a property is important

16

for the multi-point optimization.

17

0.8 0.65 0.795

0.6

Cl

Ma

0.79

0.785 0.55

Cd_2.5D Cd_2.75D Cd_3D

Cd_2.5D Cd_2.75D Cd_3D

0.78 0.5 0

18 19 20

0.001

0.002

0.003

0.004

0

ΔCd

0.0005

0.001

0.0015

0.002

0.0025

0.003

ΔCd

(a) Vs. lift coefficient

(b)Vs. Mach Number

Fig. 9. Comparison of Wing CD Variations Computed by 2.5D Transformation, 2.75D 21

1

Transformation and 3D Computation

2 3 4

2.2.5 Validation by the CRM Model

5

The wing-body configuration of the CRM [30] is also used to validate the present

6

transformation methods. The configuration was calculated at a cruise condition of

7

M=0.85, CL=0.5 and Re=5.0 million. Fig. 10(a) shows the pressure contour of the CRM

8

model. Because the pressure distribution is largely affected by the root-tip effect and

9

the kink, only the mid-outer part of the wing (from 45% to 85% of the half-span) can

10

be regarded as a quasi-conical flow region. The airfoils within this range have very little

11

difference, as shown in Fig. 10(b). Therefore, the transformation methods can be

12

applied there. The sweep angle of the 1/2 chord line is employed as the nominal sweep

13

angle for the 2.75D transformation, which is 32.6°.

14

52% of Half-Span 60% of Half-Span 68% of Half-Span 76% of Half-Span 84% of Half-Span

Y

0.05

0

-0.05

15 16 17 18

(a) Pressure Contour

0

0.2

0.4

X

0.6

0.8

1

(b) 5 Airfoils at Outer Part

Fig. 10. Pressure Contour and Airfoils of the CRM model

19

Streit et al. [16] used the pressure distribution at 60% of the half-span wing to

20

validate their method. In the present paper, the airfoil of this section is also used as the

21

baseline airfoil to test the transformations. Fig. 11 shows the pressure comparisons 22

1

among the 2.5D, 2.75D local sweep angle, 2.75D local curvature and 3D results. It is

2

clear that the results of the 2.75D local curvature method match the 3D results well. In

3

contrast, the results of the 2.5D method and the 2.75D local sweep angle method have

4

large deviations on both the lower surface and the suction peak.

5 -0.8

-0.85

-0.6

-0.8

-0.4

-0.75 -0.7

Cp

Cp

-0.2 0

-0.65 -0.6

0.2 2.5D_Method 2.75D_LocalSweep 2.75D_Curvature 3D_Method

0.4 0.6 0.8

2.5D_Method 2.75D_LocalSweep 2.75D_Curvature 3D_Method

-0.55 -0.5 -0.45

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

0

0.1

X

6 7

(a) Overall

(b) Suction Peak

-0.1

-0.4

0

-0.2

Cp

Cp

0.2

X

0.1

2.5D_Method 2.75D_LocalSweep 2.75D_Curvature 3D_Method

0

0.2 0.2

2.5D_Method 2.75D_LocalSweep 2.75D_Curvature 3D_Method

0.3

0.4 0.1

8 9 10 11 12

0.2

0.3

0.5

X

0.6

0.7

0.8

0.9

1

X

(c) Front Part of Lower Surface

(d) Aft Part

Fig. 11. Comparison of Cp Distributions Computed by the 2.5D Method, 2.75D Methods and 3D Computation

13

The key elements of a pressure distribution, such as the suction peak, the front and

14

aft loading, and the pressure recovery, are critical to the aerodynamic efficiency and

15

off-design performances of a supercritical wing. In this sense, the 2.75D_curvature

16

method, which can better reproduce the 3D pressure distribution, is more applicable for

17

the wing design. In the remainder of this paper, the 2.75D local curvature method will

18

be designated as the 2.75D method unless otherwise indicated. 23

1 2 3

3. Wing Optimization of a Wide-Body Airplane

4

In this section, the 2.75D local curvature transformation is used to optimize an

5

airfoil geared to the wing of a wide-body airplane, whose cruise condition is M=0.85,

6

CL=0.5, and Re=40.0 million. As discussed above, the “3D wing design by airfoil

7

optimization” is focused only on the outer part of the wing, which can be considered to

8

be a quasi-conical region. In the aerodynamic design of a supercritical wing, this region

9

is expected to maintain a similar supercritical pressure distribution along the span-wise

10

direction. The sweep angles of the 1/4 chord line and the 1/2 chord line are 32.52° and

11

29.41°, respectively.

12 13 14 15

3.1 Geometry Parameterization The wing shape is constructed based on 7 airfoils distributed span-wise as shown in Fig. 12. The range covered by profiles No. 4-7 is treated as the “outer part”.

16 17

Fig. 12. 3D Wing Shaped by CST Method

18 24

1

The CST (class function/shape function transformation) [24,31] technique is used

2

for the geometric parameterization. Airfoil geometries and other types of geometries

3

can be exactly represented by well-behaved and simple functions. For example, the

4

upper surface of a round nose/sharp tail 2D airfoil can be described by Eq. (27). (27)

N1 ζ Upper (ψ ) = C N2 (ψ ) ⋅ S Upper (ψ ) + ψ ⋅ ζ T / 2

5

ψ=

x ∈ (0,1) c

6

where

7

S Upper (ψ ) = ¦ Aui ⋅ Cni ⋅ψ i (1 − ψ ) n−i

n

is

the

normalized

chord-wise

coordinate,

N1 is the shape function, C N2 (ψ ) = ψ N1 ⋅ (1 − ψ ) N2 is the

i =0

8

class function, N1 = 0.5 and N 2 = 1.0 are the control coefficients, and ζ T is the

9

trailing edge thickness.

10

In this paper, a 6th-order Bernstein polynomial is used as the shape function

11

according to reference [7]. Both the upper surface and the lower surface have 7 design

12

variables. The ranges of the variables are set to be the same as those in reference [7].

13

Fig. 13 shows the initial airfoils defined by the CST method. Compared to the outer

14

part airfoils, the 3 inner airfoils have quite different shapes and relatively large

15

thicknesses.

16

To design an entire wing with good performance, iterations between the inner part

17

and the outer part of the wing are certainly needed. Because the present paper is focused

18

on the outer wing, a simplified design loop is specially constructed. It is assumed that

19

there is already a good design for the inner part. When the outer part is optimized,

20

further modifications can be applied to the inner part to eliminate the mismatch.

25

0.2

InnerPart_Section1 InnerPart_Section2 InnerPart_Section3 OuterPart_Section4567

0.15

Y

0.1 0.05 0 -0.05

1 2

0

0.2

0.4

X

0.6

0.8

1

Fig. 13. Airfoils at 7 Different Span-wise Positions

3 4

3.2 Optimization Objective and Procedures

5

Multi-point optimization is conducted to achieve a well-balanced optimal design.

6

Optimization designs based on the 2.5D and 2.75D transformations and on the 3D

7

CFD computation are conducted and compared.

8

(I) 2.5D method: The multi-point optimization is conducted on the 2D airfoil.

9

Minimizing the drag coefficients of the airfoil on multiple design points is set as the

10

objective. The pressure distribution of the airfoil is transformed onto the infinite-span

11

swept wing with the nominal sweep angle. The constraints on the pressure distribution

12

of the wing (will be discussed in section 3.3) are then checked. The CFD analysis of

13

one individual design costs approximately 2 minutes on a CPU core. The number of

14

total design variables is 7 × 2 = 14 .

26

1 2

Fig. 14. Flow Chart of the 2.75D Optimization Method

3 4

(II) 2.75D method: The pressure distribution and the drag coefficient are transformed

5

onto the tapered swept wing by the 2.75D transformation. The multi-point optimization

6

objective is set for the performance of the wing rather than the airfoil. Pressure

7

constraints are also applied on the transformed wing pressure distribution. Because the

8

CFD computation and geometry definition are still on the airfoil, the CPU cost is almost

9

the same as that of the 2.5D method. The optimization procedures are shown in Fig. 14.

10 11

Fig. 15. The Surface Mesh of the Wing-body Model and the Wing

12 27

1

(III) 3D method: Fully 3D optimization is performed directly as in reference [11].

2

A wing/body mesh of 6.5M grid points is used for the CFD analysis, as shown in Fig.

3

15. Only the geometry of the outer part wing is modified in the optimization. The

4

number of design variables is 7 × 2 × 4 = 56 . Restricted by the computational costs,

5

the fully 3D optimization performed only a single-point optimization. The objective is

6

to minimize the cruise drag coefficient. The CFD cost for the computation of a single

7

individual is 2.5 hours with 4 CPU cores.

8

The sweep angle of the 1/2 chord line is used as the nominal sweep angle for the

9

2.5D and 2.75D transformations. The span-wise lift distribution of the initial wing is

10

used to determine the lift coefficient of the cruise condition of the airfoil. Table 2 shows

11

the design conditions for the wing and the corresponding airfoil. The pressure

12

distributions of the cruise condition (a) and the lower lift condition (d) are constrained

13

during the optimization. Minimizing the drag coefficients of the (a), (b) and (c) design

14

points are set as the design objectives.

15

Table 2 Multi-point Design Conditions of the Wing and the Corresponding Airfoil Airfoil design conditions Wing flight conditions

(for 2.5D and 2.75D)

(a) Cruise condition

M=0.85᧨CL=0.5

M=0.741᧨CL=0.82

(b) Lower Mach number

M=0.83᧨CL=0.5

M=0.723᧨CL=0.82

(c) Higher Mach number

M=0.87᧨CL=0.5

M=0.758᧨CL=0.82

(d) Lower lift coefficient

M=0.85᧨CL=0.45

M=0.741᧨CL=0.74

16 17 18 19

3.3 Constraint Settings According to reference [11], the optimization constraints consist of pressure distribution constraints and geometry constraints, as denoted in Fig. 16.

28

1 Fig. 16. Pressure Distribution Constraints and Geometry Constraints

2 3 4

(I) Pressure distribution constraints:

5

(1) Suction peak: the Cp of the suction peak should not be lower than -1.3 to avoid

6

excessive acceleration at the leading edge and to prevent the airfoil from becoming a

7

“peaky” type airfoil;

8

(2) Suction platform: a low pressure gradient (-0.2 < dCp/dx < 0.5) should be

9

maintained within the range of 20 % < x < 50% to obtain sufficient lift and to maintain

10

good stability of the shock wave;

11

(3) Shock wave location: the shock wave should be located between 45% and 55%

12

of the chord to achieve an adequate buffet margin. This constraint is applied only on

13

the cruise condition;

14 15

(4) Shock wave strength: the Cp increase across the shock wave should be less than 0.55 to limit the shock wave drag;

16

(5) Pressure recovery: the pressure gradient should be less than 3 (dCp/dx < 3) at the

17

rear part of the upper surface (70% < x < 90%) to restrict the trailing edge flow

18

separation; 29

1 2

(6) Aft-loading: the difference in Cp between the upper surface and lower surface at 70% of the chord length should be less than 0.7 to restrict the pitching moment;

3

(7) Maximum speed point of the lower surface: the suction peak of the lower surface

4

should be higher than -0.4 to avoid the supersonic zone appearing on the lower surface

5

at higher Mach numbers.

6

(II) Geometry constraints

7 8 9

(8) Leading edge radius: the leading edge radius should be no smaller than the original value to maintain good low-speed performance. (9) Thickness: the thickness of the corresponding airfoil is fixed to 10.9%.

10 11

3.4 Optimization Algorithm Settings

12

The NSGA-II is used as the optimization method. The parameters of NSGA-II are

13

set as follows: the crossover probability is 0.8, and the mutation probability is 1/n

14

(where n is the number of design variables). For the 2.5D and 2.75D optimizations, all

15

the objectives have the same weight to guarantee a well-balanced result. Each

16

generation of the genetic algorithm has 32 individuals.

17 18

3.4.1 Multi-Point Optimization History

19

The 2.5D and 2.75D optimizations are conducted on a 32-core workstation. The

20

airfoil of the initial wing and the RAE 2822 airfoil, together with other random

21

individuals, form the first generation’s population. With the same initial population, the

22

2.5D optimization and the 2.75D optimization both iterate 100 generations and cost

23

approximately 15 hours. Fig. 17 shows the 2.75D method’s convergence histories of

24

the drag coefficients. The drag coefficients at all 4 design points are obviously

25

decreased. 30

1 M = 0.723 CL = 0.82

M = 0.741 CL = 0.74

0 .0 1 1 5 0 .0 0 9 6

'Cd'

'Cd'

0 .0 1 1 0 .0 0 9 4

0 .0 1 0 5 0 .0 0 9 2

0 .0 1 0 .0 0 9 20

40

'G e n e r a t i o n '

60

80

20

40

'G e n e r a t io n '

60

M = 0.741 CL = 0.82

0 .0 1 3 5

80

M = 0.758 CL = 0.82 0 .0 1 6

0 .0 1 3 0 .0 1 2 5

'Cd'

'Cd'

0 .0 1 2

0 .0 1 4

0 .0 1 1 5 0 .0 1 1 0 .0 1 0 5

2 3

0 .0 1 2 20

40

'G e n e r a t i o n '

60

80

20

40

'G e n e r a t io n '

60

80

Fig. 17. Convergence Histories of CD for the 4 Design Points

4 5

3.4.2 Pareto Front Analysis

6

Fig. 18 shows the Pareto fronts of the 2.75D method (black dots) and the 2.5D

7

method (green tri-angles) of the CD at the design point (M=0.758 and CL=0.82) and the

8

lower Mach number point (M=0.741 and CL=0.82). The drag of the 2.5D method results

9

is re-calculated by Eq. (26) in Fig. 18 for comparison purposes.

10 11

Fig. 18. Pareto Fronts of CD for Drag Divergence Point and Cruise Point

12

31

-1

Ma = 0.723 Cl = 0.85

-0.8

-0.6

-0.6

-0.4

-0.4

-0.2

Cp

Cp

-0.2 0

0.4

0.4

Initial_Airfoil 2.5D_Method 2.75D_Method

0.6

Initial_Airfoil 2.5D_Method 2.75D_Method

0.6 0.8

0.8 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

0

X

1

1

Ma = 0.758 Cl = 0.85

-0.8 -0.6

-0.4

-0.4

-0.2

-0.2

Cp

-0.6

Cp

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

X

Ma = 0.741 Cl = 0.85

-0.8

0

0 0.2

0.2 0.4

0.4

Initial_Airfoil 2.5D_Method 2.75D_Method

0.6 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Initial_Airfoil 2.5D_Method 2.75D_Method

0.6 0.8

0.8

3

0 0.2

0.2

2

Ma = 0.741 Cl = 0.77

-0.8

1

0

X

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

X

Fig. 19. Optimized Cp Distributions at 4 Design Points by the 2.5D Method and the 2.75D Method

4 5

The Pareto fronts can be clearly observed in Fig. 18. Because of the weights

6

changing the pressure drag and friction drag, the importance of the higher Mach number

7

or the drag divergence (M = 0.758 and CL=0.82) point is lowered in the 2.75D

8

optimization. Moreover, the importance of the cruise condition is increased. The

9

optimization tends to better improve the cruise lift/drag ratio while appropriately

10

sacrificing the lift/drag ratio at the higher Mach number point. The Pareto front of the

11

2.75D optimization clearly tends to approach the left in Fig. 18. Therefore, the designs

12

from the 2.75D method tend to show better performance for the cruise condition.

13

Two solutions on their Pareto fronts are selected and tested as the optimal designs.

14

They are highlighted by red and dark edges, respectively. Fig. 19 shows the pressure

15

distribution comparison among the two selected optimal airfoils and the initial airfoil. 32

1

The results show that with the pressure distribution constraints, both of the methods can

2

obtain improved designs with desired supercritical pressure distributions.

3

Table 3. Performances of Airfoils Compared between 2 Optimization Methods

Initial Airfoil 2.5D Method 2.75D Method

Ma

CL_2D

CDp_2D

CDf_2D

CD_2D

CL/CD

CD_2.5D

CD_2.75D

0.723

0.820

0.00803

0.00501

0.01304

65.2

0.00990

0.01018 0.01061

0.741

0.820

0.00872

0.00498

0.01370

62.0

0.01039

0.758

0.820

0.01139

0.00502

0.01641

51.8

0.01246

0.01242

0.723

0.820

0.00615

0.00467

0.01082

78.6

0.00821

0.00861 0.00879

0.741

0.820

0.00643

0.00467

0.01110

75.2

0.00842

0.758

0.820

0.00877

0.00463

0.01340

63.4

0.01017

0.01030

0.723

0.820

0.00552

0.00471

0.01023

83.1

0.00776

0.00823 0.00856 0.01035

0.741

0.820

0.00604

0.00469

0.01073

79.2

0.00814

0.758

0.820

0.00883

0.00463

0.01346

63.2

0.01021

4 5

Table 3 presents the aerodynamic coefficients of the three airfoils. CL_2D is the

6

lift coefficient of the 2D airfoil. CDp_2D, CDf_2D and CD_2D are the pressure drag, the

7

friction drag and the total drag coefficients of the airfoil, respectively. CD_2.75D is the

8

drag coefficient transformed by Eq. (26). The optimizations based on the 2.5D and

9

2.75D transformations both achieve improvements on the initial design. At the cruise

10

point, the lift/drag ratios are increased from 62.0 to 75.2 for the 2.5D result and to 79.2

11

for the 2.75D result.

12 13

3.5 Integrated Verification of the Designed Wings

14

The optimized 2.5D and 2.75D airfoils are installed onto the outer part of the wing

15

and verified by 3D RANS computation. The pressure distributions at 60% of the half-

16

span location are compared with the transformations from the corresponding airfoils.

17

Fig. 20 presents the comparisons. The results of the 2.5D optimized wing show large

18

discrepancies on both the upper and lower surfaces. There is a tendency of forming a

19

secondary shock wave in the 3D results. Apparently, the 2.5D transformation from the

20

airfoil fails in predicting this tendency. In contrast, the results of the 2.75D-optimized

21

wing show considerably better agreement. There is only a slight difference after the 33

1

shock wave. This difference might be caused by the effects of the boundary layer. The

2

2.75D method is able to better pursue the performance objectives and better follow the

3

constraints on the pressure distributions.

4

The wings obtained by the 2.5D, the 2.75D and the 3D methods are numerically

5

analyzed by 3D RANS computation and compared at the cruise point with a fine grid

6

of 20 million cells. The pressure coefficients of the sections at 40%, 60% and 85% of

7

the half-span are compared in Fig. 21.

8

2.5D

-0.8

-0.4

-0.2

-0.2

Cp

-0.4

Cp

-0.6

0 0.2

0 0.2

Result_Airfoil Result_Wing

0.4

Result_airfoil Result_Wing

0.4

0.6

0.6 0

9 10 11

2.75D

-0.8

-0.6

0.2

0.4

0.6

0.8

1

0

0.2

0.4

X

(a) Wing-body Configuration

0.6

0.8

1

X

(b) 2.5D Result

(c) 2.75D Result

Fig. 20. Comparison of the Pressure Distribution Results by the 2.5D and 2.75D Methods

12 40%

60%

-0.6

-0.2

Cp

Cp

-0.4

0

0.4

13 14 15

0.2

0.4

0.6

-0.6

-0.6

-0.4

-0.4

-0.2

-0.2

0

2.5D 2.75D 3D

0.4

0.8

X

(a) 40% of the Half-span

1

85%

0

2.5D 2.75D 3D

0.2 0.4 0.6

0.6

0.6 0

-0.8

0.2

2.5D 2.75D 3D

0.2

-0.8

Cp

-1 -0.8

0

0.2

0.4

0.6

0.8

X

(b) 60% of the Half-span

1

0

0.2

0.4

0.6

0.8

1

X

(c) 85% of the Half-span

Fig. 21. Comparison of the Pressure Distributions among the 2.5D, 2.75D and 3D Methods

16 17

In the 2.5D method results, an excessive acceleration region after the shock wave

18

appears at 60% of the half-span location. The shock position at 85% of the half-span

19

section is too close to the trailing edge, which might lead to unexpected flow separation 34

1

at an off-design condition. Because of the error in the 2.5D transformation, the pressure

2

distribution constraints in the optimization failed to detect and suppress these.

3

The pressure distributions of the 3D optimized wing do not show these types of

4

acceleration regions. However, the pressure distributions of the three sections on the

5

lower surface are not consistent, which indicates that the span-wise smoothness of the

6

wing might not be good enough and will cause problems in the cross flow and in

7

manufacturing.

8 9

Comparatively, the pressure distributions from the 2.75D result are in good shape and demonstrate better span-wise consistency.

10

Table 4 shows the drag coefficients of the wing-body configuration with the

11

different wings. The wings optimized by the 2.75D and 3D methods have similar

12

performances and are both better than the 2.5D result. The drag difference between the

13

2.75D and 3D wings is less than 1 drag count. However, the CPU time costs of the two

14

methods are quite different. With the same 32-core workstation, both run for 100

15

generations: the 2.75D method spends 15 hours for a multi-point optimization, while

16

the 3D method costs more than 30 days for a single-point optimization. The 2.75D

17

transformation dramatically improves the efficiency of the optimization.

18

Table 4. Comparison of Drag Coefficients among 3 Optimization Methods

Initial design

M=0.83 CL=0.5

M=0.85 CL=0.45

M=0.85 CL=0.5

M=0.87 CL=0.5

0.02278

0.02106

0.02354

0.02703

2.5D Method

0.02119

0.01973

0.02143

0.02352

2.75D Method

0.02066

0.01917

0.02083

0.02284

3D Method

0.02059

0.01922

0.02077

0.02294

19 20

Fig. 22 presents the pressure contours of the 4 wings. This figure shows that the

21

2.75D-optimized wing has the smoothest span-wise variation of the pressure

22

distribution. Moreover, Fig. 23 compares the buffet and drag divergence performances 35

1

of the 4 wings. With satisfactory application of the well-designed pressure distribution

2

constraints, the 2.75D-optimized wing achieved good performance at both the design

3

and off-design points.

4

5 6

7 8 9

(a) Initial Wing

(b) 2.5D Optimized Result

(c) 2.75D Optimized Result

(d) 3D Optimized Result

Fig. 22. Comparison of the Pressure Contours of Wings Optimized by Different Methods

10

36

0.027

Initial 2.5D 2.75D 3D

0.0265 0.7

0.026 0.0255 0.025

0.6

0.024

0.5

Cd

Cl

0.0245

0.0235

Initial 2.5D 2.75D 3D

0.4

0.3

0.023 0.0225 0.022 0.0215 0.021

0.2 0

1 2 3

1

2

Alfa

3

4

5

0.83

(a) Buffet Performance

0.835

0.84

0.845

0.85

Ma

0.855

0.86

0.865

0.87

(b) Drag Divergence Performance

Fig. 23. Comparison of the Off-design Performances of Wings Optimized by Different Methods

4 5 6

4. Conclusion

7

This paper introduces a design method for supercritical wings. By restricting the

8

CFD analysis on only the airfoil, superior optimization efficiency is achieved. With a

9

more accurate 2D/3D transformation and appropriate application of constraints on the

10

pressure distribution, satisfactory design-point and off-design performances of

11

supercritical wing can be achieved. The contributions of this paper can be summarized

12

as follows:

13

(1) As the core of the design method, a 2.75D transformation is introduced for both the

14

pressure distribution and the lift/drag coefficients. The effects of the taper ratio, the

15

sweep angle and the surface curvature are better taken into account. The sweep angle

16

of the 1/2 chord line is suggested to be the nominal sweep angle. Different correction

17

factors are proposed for the pressure drag and friction drag. Test cases verified the

18

transformation’s satisfactory accuracy and better applicability for the multi-point

19

optimization of the supercritical wing.

20

(2) Constraints are applied on the pressure distribution during the optimization. With 37

1

the 2.75D transformation, these constraints can well control the pressure distribution of

2

the wing and make it capable of achieving satisfactory "supercritical" properties.

3

Although the method is only applied to the outer part, the designed wing can have great

4

improvements in both design and off-design performances.

5

(3) For the large aspect ratio wing design, such a philosophy of “2D Optimization for

6

3D Wing” can improve not only the design’s efficiency and quality but also the

7

consistency of the geometry. This will considerably help in forming the geometry and

8

in manufacturing.

9 10

38

1

References

2

[1]

Aircraft, IOS Press, Amsterdam, 2009, pp. 135-197.

3 4

E. Obert, Swept wing, in: E. Obert (Ed.), Aerodynamic Design of Transport

[2]

J.C. Vassberg, M.A. DeHaan, S.M. Rivers, R.A. Wahls, Development of a

5

common research model for applied CFD validation studies, in: 26th AIAA

6

Applied Aerodynamics Conference, 2008, AIAA 2008-6919.

7

[3]

application of CFD at Boeing commercial airplanes, Seattle. AIAA 2003-3439.

8 9

[4]

[5]

[6]

A. Jameson, S. Shankaran, L. Martinelli, Continuous adjoint method for unstructured grids, AIAA J. 46 (2008) 1226-1239.

14 15

A. Jameson, S. Kim, Reduction of the adjoint gradient formula for aerodynamic shape optimization problems, AIAA J. 41 (2003) 2114-2129.

12 13

P. Edi, J.P. Fielding, Civil-transport wing design concept exploiting new technologies, J. Aircraft 43 (2006) 932-940.

10 11

F.T. Johnson, E.N. Tinoco, N.J. Yu, Thirty years of development and

[7]

Y. Zhang, X. Fang, H. Chen, S. Fu, Z. Duan, Y. Zhang, Supercritical natural

16

laminar flow airfoil optimization for regional aircraft wing design, Aerosp. Sci.

17

Technol. 43 (2015) 152-164.

18

[8]

genetic algorithm: NSGA-II, IEEE Trans. Evol. Computat. 6 (2002) 182-197.

19 20 21

K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective

[9]

W. Li, S. Krist, R. Campbell, Transonic airfoil shape optimization in preliminary design environment, J. Aircraft 43 (2006) 639-651. 39

1

[10]

transonic airfoil design, J. Aircraft 44 (2007) 365-376.

2 3

A. Vavalle, N. Qin, Iterative response surface based optimization scheme for

[11]

Y.F. Zhang, H.X. Chen, M. Zhang, M.H. Zhang, T.J. Liu, W.S. Zhang, S. Fu,

4

Supercritical wing design and optimization for transonic civil airplane, in: 49th

5

AIAA Aerospace Sciences Meeting including the New Horizons Forum and

6

Aerospace Exposition, 2011, AIAA 2011-27.

7

[12]

distributions, RAE R.&M 3346 (1962).

8 9

R.C. Lock, An equivalence law relating three- and two-dimensional pressure

[13]

A.J.M. Van der Velden, I. Kroo, Numerical method for relating two- and three-

10

dimensional pressure distributions on transonic wings, in: Aircraft Design,

11

Systems and Operations Conference, AIAA-90-3211-CP.

12

[14]

Aircr. 38 (2001) 778-782.

13 14

N. Petruzzelli, A.J. Keane, Wave drag estimation for use with panel codes. J.

[15]

C.J. Atkin, E.R. Gowree, Recent developments to the viscous Garabedian and

15

Korn method, in: 28th International Congress of the Aeronautical Sciences,

16

2012.

17

[16]

T. Streit, G. Wichmann, F. Knoblauch zu Hatzbach, R.L. Campbell,

18

Implications of conical flow for laminar wing design and analysis, in: 29th

19

AIAA Applied Aerodynamics Conference, 2011, AIAA 2011-3808.

20 21

[17]

Y.F. Zhang, H.X. Chen, S. Fu, Improvement to patched grid technique with high-order conservative remapping method, J. Aircraft 48 (2011) 884-893. 40

1

[18]

Large Eddy Simulation, Comput. Fluids 117 (2015) 233–246.

2 3

Z. Li, Y.F. Zhang, H.X. Chen, A low dissipation numerical scheme for implicit

[19]

H.X. Chen, X. Huang, K. Shi, S. Fu, M. Ross, M.A. Bennington, J.D. Cameron,

4

S.C. Morris, S. McNulty A. Wadia, A computational fluid dynamics study of

5

circumferential groove casing treatment in a transonic axial compressor, J.

6

Turbomachinery 136 (2014) 031003.

7

[20]

Y. Zhang, H. Chen, M. Zhang, M. Zhang, Z. Li, S. Fu, Performance prediction

8

of conical nozzle using Navier–Stokes computation, J. Propul. Power 31 (2015)

9

192-203.

10

[21]

dimensional Euler solutions, J. Aircraft 53 (2016) 294–299.

11 12

Z. Li, H. Chen, Y. Zhang, S. Fu. (2016), Grid-convergence study of two-

[22]

P.H. Cook, M.A. Mcdonald, M.C.P. Pirmin, Aerofoil RAE 2822-pressure

13

distributions, and boundary layer and wake measurements. Experimental Data

14

Base for Computer Program Assessment. AGARD Advisory Report No.138.

15

pp. A6-1~76.

16

[23]

M.B. Rivers, A. Dittberner, Experimental investigation of NASA common

17

research model (Inveted), in: 28th AIAA Applied Aerodynamics Conference,

18

2010, AIAA 2010-4218.

19

[24]

T. Zhao, Y.F. Zhang, H.X. Chen, Multi-objective aerodynamic optimization of

20

supercritical wing with substantial pressure constraints, in: 53rd AIAA

21

Aerospace Sciences Meeting, 2015, AIAA 2015-0763. 41

1

[25]

University Press, 2007, Section 3, Chapter 6, pp. 490-492. (In Chinese).

2 3

F.W. Li, Airfoil design, in: Aerodynamic Theory. Northeast Poly-Technical

[26]

P.M. Sforza, Wing design, in: P.M. Sforza (Ed.),Commercial Airplane Design Principles, vol. 4, Elsevier, Amsterdam, 2014, pp. 119-210.

4 5

[27]

C.D. Harris, NASA supercritical airfoils, NASA Technical Paper 2969, 1990.

6

[28]

S. Painchaud-Ouellet, C. Tribes, J.Y. Trepanier, D. Pelletier, Airfoil shape

7

optimization using a non-uniform rational B-splines parameterization under

8

thickness constraint. AIAA J. 44 (2006) 2170-2178.

9

[29]

optimization, J. Aircraft 50 (2013) 1046–1059.

10 11

A. Elham, M.J.L. van Tooren, Weight indexing for airfoil multi-objective

[30]

A.J. Sclafani, J.C. Vassberg, C. Winkler, A.J. Dorgan, M. Mani, M.E. Olsen,

12

J.G. Coder, DPW-5 analysis of the CRM in a wing-body configuration using

13

structured and unstructured meshes, in: 51st AIAA Aerospace Sciences Meeting

14

including the New Horizons Forum and Aerospace Exposition, 2013, AIAA-

15

2013-0048.

16 17

[31]

B.M. Kulfan, Universal parametric geometry representation method, J. Aircraft 45 (2008) 142-158.

18 19

42