Airfoil shape optimization using adjoint method and automatic - - PowerPoint PPT Presentation

Airfoil shape optimization using adjoint method and automatic differentiation

• Praveen. C

praveen@math.tifrbng.res.in

Tata Institute of Fundamental Research Center for Applicable Mathematics Bangalore 560065 http://math.tifrbng.res.in

AeSI CFD Symposium 11-12 August, 2009

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 1 / 25

Objectives and controls

• Objective function I(β) = I(β, Q)

mathematical representation of system performance

• Control variables β

◮ Parametric controls β ∈ Rn ◮ Infinite dimensional controls β : X → Y ◮ Shape β ∈ set of admissible shapes

• State variable Q: solution of an ODE or PDE

R(β, Q) = 0 = ⇒ Q = Q(β)

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 2 / 25

Mathematical formulation

• Constrained minimization problem

min

β

I(β, Q) subject to R(β, Q) = 0

• Find δβ such that δI < 0 (to first order)

δI = ∂I ∂β δβ + ∂I ∂QδQ = ∂I ∂β + ∂I ∂Q ∂Q ∂β

• G

δβ

• Steepest descent

δβ = −ǫG⊤, ǫ > 0 δI = −ǫGG⊤ = −ǫG2 < 0 How to compute gradient G cheaply and accurately ?

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 3 / 25

Elements of shape optimization

1 Shape parameterization 2 Surface grid generation/deformation 3 Domain grid generation/deformation 4 Flow solution (Euler/Navier-Stokes solver) 5 Adjoint flow solution 6 Optimization method

Shape parameters β Surface grid Xs Volume grid X CFD solution Q I

dI dβ = dI dX dX dXs dXs dβ

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 4 / 25

Adjoint approach

• For shape optimization: I = I(X, Q)

dI dX = ∂I ∂X + ∂I ∂Q ∂Q ∂X

• Flow sensitivity ∂Q

∂X ; costly to evaluate

• Differentiate state equation R(X, Q) = 0

∂R ∂X + ∂R ∂Q ∂Q ∂X = 0

• Introducing an adjoint variable Ψ, we can write

dI dX = ∂I ∂X + ∂I ∂Q ∂Q ∂X + Ψ⊤ ∂R ∂X + ∂R ∂Q ∂Q ∂X

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 5 / 25

Adjoint approach

• Collect terms involving the flow sensitivity

dI dX = ∂I ∂X + Ψ⊤ ∂R ∂X + ∂I ∂Q + Ψ⊤ ∂R ∂Q ∂Q ∂X

• Choose Ψ so that flow sensitivity vanishes

∂I ∂Q + Ψ⊤ ∂R ∂Q = 0

• r

∂R ∂Q ⊤ Ψ + ∂I ∂Q ⊤ = 0

• Gradient

dI dX = ∂I ∂X + Ψ⊤ ∂R ∂X

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 6 / 25

Optimization steps

• β =

⇒ Xs = ⇒ X

• Solve the flow equations to steady-state

dQ dt + R(X, Q) = 0 = ⇒ Q, I(X, Q)

• Solve adjoint equations to steady-state

dΨ dt + ∂R ∂Q ⊤ Ψ + ∂I ∂Q ⊤ = 0 = ⇒ Ψ

• Compute gradient wrt grid X

dI dX = ∂I ∂X + Ψ⊤ ∂R ∂X dI dβ = dI dX dX dXs dXs dβ = ⇒ β ← − β − ǫdI dβ

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 7 / 25

Continuous and discrete approaches

• Continuous approach (differentiate and discretize)

PDE − → Adjoint PDE − → Discrete adjoint

• Discrete approach (discretize and differentiate)

PDE − → Discrete PDE − → Discrete adjoint

• We use the discrete approach

◮

R(X, Q) = 0 represent the finite volume equations which are algebraic equations

◮ Use ordinary calculus to differentiate ◮ Need to compute

∂I ∂Q, ∂I ∂X , ∂R ∂Q ⊤ Ψ, ∂R ∂X ⊤ Ψ

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 8 / 25

Automatic differentiation

• Computer code available to compute

I(X, Q), R(X, Q)

• Code is made of composition of elementary functions

T0 = X, Tr = Fr(Tr−1) Y = F(X) = Fp ◦ Fp−1 ◦ . . . ◦ F1(T0)

◮ Use differentiation by parts formula

˙ Y = F ′(X) ˙ X = F ′

p(Tp−1)F ′ p−1(Tp−2) . . . F ′ 1(T0) ˙

X

◮ Automated using AD tools

Computer code P Automatic Differentiation New code ˙ P

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 9 / 25

Reverse differentiation

• Reverse mode computes transpose: (X, ¯

Y ) − → ¯ X ¯ X = [F ′(X)]⊤ ¯ Y = [F ′

1(T0)]⊤[F ′ 2(T1)]⊤ . . . [F ′ p(Tp−1)]⊤ ¯

Y

• Forward sweep and then reverse sweep

Func: T0 − → F1

T1

− → F2

T2

− → . . .

Tp−1

− → Fp Grad: Tp−1, ¯ Y − → [F ′

p]T ¯ Tp−2

− → [F ′

p−1]T ¯ Tp−3

− → . . .

¯ T0

− → [F ′

1]T

Forward variables Tj required in reverse order: store or recompute

• Reverse mode useful to compute

∂I ∂Q ⊤ , ∂I ∂X ⊤ , ∂R ∂Q ⊤ Ψ, ∂R ∂X ⊤ Ψ

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 10 / 25

Differentiation: Example

• A simple example

f = (xy + sin x + 4)(3y2 + 6)

• Computer code, f = t10

t1 = x t2 = y t3 = t1t2 t4 = sin t1 t5 = t3 + t4 t6 = t5 + 4 t7 = t2

2

t8 = 3t7 t9 = t8 + 6 t10 = t6t9

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 11 / 25

F77 code: costfunc.f

subroutine costfunc (x , y , f ) t1 = x t2 = y t3 = t1 ∗ t2 t4 = sin ( t1 ) t5 = t3 + t4 t6 = t5 + 4 t7 = t2 ∗∗2 t8 = 3.0∗ t7 t9 = t8 + 6.0 t10 = t6 ∗ t9 f = t10 end

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 12 / 25

Automatic Differentiation: Reverse mode

SUBROUTINE COSTFUNC_B(x, xb, y, yb, f, fb) t1 = x t2 = y t3 = t1*t2 t4 = SIN(t1) t5 = t3 + t4 t6 = t5 + 4 t7 = t2**2 t8 = 3.0*t7 t9 = t8 + 6.0 t10b = fb t6b = t9*t10b t9b = t6*t10b t8b = t9b t7b = 3.0*t8b t5b = t6b t3b = t5b t2b = t1*t3b + 2*t2*t7b t4b = t5b t1b = t2*t3b + COS(t1)*t4b yb = t2b xb = t1b fb = 0.0 END

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 13 / 25

Implementation of AD

• CFD code written with many subroutines
• Subroutines differentiated individually
• Then assembled together to form adjoint solver
• Only non-linear portions differentiated with AD

◮ numerical flux (Roe) ◮ Limiters

• Linear portions differentiated manually
• Leads to an efficient code with less memory requirements
• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 14 / 25

Shape parameterization

• Parameterize the

deformations xs ys

• =
• x(0)

s

y(0)

s

• +

nx ny

• h(ξ)

h(ξ) =

m

• k=1

βkBk(ξ)

• Hicks-Henne bump functions

Bk(ξ) = sinp(πξqk), qk = log(0.5) log(ξk)

• Move points along normal to

reference line AB

A B

• n

ξ

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ξ h(ξ)

Exact derivatives dXs

dβ can be

computed

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 15 / 25

Grid deformation

• Interpolate displacement of

surface points to interior points using RBF ˜ f(x, y) = a0 + a1x + a2y +

N

• j=1

bj| r − rj|2 log | r − rj| where

• r = (x, y)
• Results in smooth grids
• Exact derivatives

dX dXs can be

computed Initial grid Deformed grid

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 16 / 25

NUWTUN flow solver

Based on the ISAAC code of Joseph Morrison http://isaac-cfd.sourceforge.net

• Finite volume scheme
• Structured, multi-block grids
• Roe flux
• MUSCL reconstruction with Koren limiter
• Implicit scheme

Source code of NUWTUN available online http://nuwtun.berlios.de

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 17 / 25

Convergence tests

500 1000 1500 2000 2500

Number of iterations

1x10-15 1x10-14 1x10-13 1x10-12 1x10-11 1x10-10 1x10-9 1x10-8 1x10-7 1x10-6 0.00001 0.0001 0.001 0.01 0.1 1 10 100

Residue

Flow residual Adjoint residual

500 1000 1500 2000 2500 3000

Number of iterations

0.2 0.4 0.6 0.8

Cl

0.02 0.04 0.06 0.08 0.1

Cd

Cl Cd

Convergence characteristics for the flow and adjoint solutions, and convergence of lift and drag coefficients, for RAE2822 airfoil at M∞ = 0.73

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 18 / 25

Validation of adjoint gradients

• Dot-product test
• Limiters can cause

non-differentiability

• Koren limiter: dependance
• n parameter
• Check adjoint derivatives

against finite difference

5 10 15 20

Hicks-Henne parameter

• 40
• 20

20 40 60 80

Gradient

AD FD 5 10 15 20

Hicks-Henne parameter

• 100
• 50

50 100 150

Gradient

AD FD

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 19 / 25

Test cases

• NACA0012: M∞ = 0.8, α = 1.25o

I = Cd Cl Cl0 Cd0

• RAE2822: M∞ = 0.729, α = 2.31o

◮ Penalty approach

I = Cd Cd0 +

• 1 − Cl

Cl0

• ◮ Constrained minimization

min I = Cd Cd0 s.t. Cl = Cl0

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 20 / 25

NACA0012: Maximize L/D

0.2 0.4 0.6 0.8 1

x/c

• 1
• 0.5

0.5 1

• Cp

Initial conmin_frcg

• ptpp_q_newton

steep 0.2 0.4 0.6 0.8 1

• 0.06
• 0.04
• 0.02

0.02 0.04 0.06 Initial conmin_frcg

• ptpp_q_newton

steep

Method I 100Cd Cl Nfun Ngrad Initial 1.000 2.072 0.295

• conmin frcg

0.170 0.389 0.325 50 13

• ptpp q newton

0.142 0.329 0.329 50 39 steep 0.166 0.335 0.287 50 45

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 21 / 25

RAE2822: Drag minimization, penalty approach

0.2 0.4 0.6 0.8 1

x

• 1
• 0.5

0.5 1 1.5

• Cp

Initial conmin_frcg

• ptpp_q_newton

steep 0.2 0.4 0.6 0.8 1

x

• 0.06
• 0.04
• 0.02

0.02 0.04 0.06

y

Initial conmin_frcg

• ptpp_q_newton

steep

Method I 100Cd Cl Nfun Ngrad Initial 1.000 1.150 0.887

• conmin frcg

0.355 0.405 0.890 50 13

• ptpp q newton

0.351 0.400 0.884 50 51 steep 0.341 0.388 0.884 50 47

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 22 / 25

RAE2822: Lift-constrained drag minimization

0.2 0.4 0.6 0.8 1

x/c

• 1
• 0.5

0.5 1 1.5

• Cp

Initial conmin_mfd fsqp ipopt 0.2 0.4 0.6 0.8 1

• 0.06
• 0.04
• 0.02

0.02 0.04 0.06 Initial conmin_mfd fsqp ipopt

Method I 100Cd Cl Nfun Ngrad Initial 1.000 1.150 0.887

• conmin mfd

0.342 0.394 0.881 50 22 fsqp 0.325 0.374 0.887 50 32 ipopt 0.335 0.385 0.887 45 62

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 23 / 25

RAE2822: Lift-constrained drag minimization

0.2 0.4 0.6 0.8 1

x/c

• 1
• 0.5

0.5 1 1.5

• Cp

ipopt Initial 0.2 0.4 0.6 0.8 1

• 0.06
• 0.04
• 0.02

0.02 0.04 0.06 ipopt Initial

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 24 / 25

Sensitivity to perturbations

0.7 0.71 0.72 0.73 0.74 0.75 0.76 Mach number 0.01 0.02 0.03 0.04 Drag coefficient Initial Optimized 0.7 0.71 0.72 0.73 0.74 0.75 0.76 Mach number 50 100 150 200 Lift/Drag Initial Optimized

(a) (b) Variation of (a) drag coefficient and (b) L/D with Mach number for RAE2822 airfoil and optimized airfoil Need for robust aerodynamic optimization

• Praveen. C

(TIFR-CAM) Shape Optimization AeSI, Aug 2009 25 / 25