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 β ∈ R n ◮ 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 ∂β δβ + ∂ I ∂β + ∂ I ∂Q ∂QδQ = δβ ∂Q ∂β � �� � G • Steepest descent δβ = − ǫG ⊤ , ǫ > 0 δ I = − ǫGG ⊤ = − ǫ � G � 2 < 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 Volume grid CFD solution I Q β X s X d β = d I d I d X d X s d X d X s d β Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 4 / 25

Adjoint approach • For shape optimization: I = I ( X, Q ) d X = ∂ I d I ∂X + ∂ I ∂Q ∂Q ∂X • Flow sensitivity ∂Q ∂X ; costly to evaluate • Differentiate state equation R ( X, Q ) = 0 ∂X + ∂R ∂R ∂Q ∂X = 0 ∂Q • Introducing an adjoint variable Ψ, we can write � ∂R � d X = ∂ I d I ∂X + ∂ I ∂Q ∂X + ∂R ∂Q ∂X + Ψ ⊤ ∂Q ∂Q ∂X Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 5 / 25

Adjoint approach • Collect terms involving the flow sensitivity � ∂ I � ∂Q d X = ∂ I d I ∂X + Ψ ⊤ ∂R ∂Q + Ψ ⊤ ∂R ∂X + ∂Q ∂X • Choose Ψ so that flow sensitivity vanishes � ∂ I � ∂R � ⊤ � ⊤ ∂Q + Ψ ⊤ ∂R ∂ I ∂Q = 0 or Ψ + = 0 ∂Q ∂Q • Gradient d X = ∂ I d I ∂X + Ψ ⊤ ∂R ∂X Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 6 / 25

Optimization steps • β = ⇒ X s = ⇒ X • Solve the flow equations to steady-state d Q d t + R ( X, Q ) = 0 = ⇒ Q, I ( X, Q ) • Solve adjoint equations to steady-state � ∂ I � ∂R � ⊤ � ⊤ dΨ d t + Ψ + = 0 = ⇒ Ψ ∂Q ∂Q • Compute gradient wrt grid X d X = ∂ I d I ∂X + Ψ ⊤ ∂R ∂X d β = d I d I d X d X s − β − ǫ d I = ⇒ β ← d X d X s d β 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 � ∂R � ∂R � ⊤ � ⊤ ∂ I ∂ I ∂Q, ∂X , Ψ , Ψ ∂Q ∂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 T 0 = X, T r = F r ( T r − 1 ) Y = F ( X ) = F p ◦ F p − 1 ◦ . . . ◦ F 1 ( T 0 ) ◮ Use differentiation by parts formula Y = F ′ ( X ) ˙ ˙ 1 ( T 0 ) ˙ X = F ′ p ( T p − 1 ) F ′ p − 1 ( T p − 2 ) . . . F ′ X ◮ Automated using AD tools Computer code New code Automatic ˙ Differentiation P P Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 9 / 25

Reverse differentiation • Reverse mode computes transpose: ( X, ¯ → ¯ Y ) − X X = [ F ′ ( X )] ⊤ ¯ 2 ( T 1 )] ⊤ . . . [ F ′ p ( T p − 1 )] ⊤ ¯ ¯ Y = [ F ′ 1 ( T 0 )] ⊤ [ F ′ Y • Forward sweep and then reverse sweep T p − 1 T 1 T 2 Func: T 0 − → F 1 − → F 2 − → . . . − → F p ¯ ¯ ¯ T p − 2 T p − 3 T 0 T p − 1 , ¯ [ F ′ p ] T [ F ′ p − 1 ] T [ F ′ 1 ] T Grad: Y − → − → − → . . . − → Forward variables T j required in reverse order: store or recompute • Reverse mode useful to compute � ∂ I � ∂ I � ∂R � ⊤ � ⊤ � ∂R � ⊤ � ⊤ , , Ψ , Ψ ∂Q ∂X ∂Q ∂X Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 10 / 25

Differentiation: Example • A simple example f = ( xy + sin x + 4)(3 y 2 + 6) • Computer code, f = t 10 t 1 = x t 2 = y t 3 = t 1 t 2 t 4 = sin t 1 t 5 = t 3 + t 4 t 6 = t 5 + 4 t 2 t 7 = 2 t 8 = 3 t 7 t 9 = t 8 + 6 t 10 = t 6 t 9 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 n � deformations � � � x s � � n x � A x (0) s = + h ( ξ ) ξ y (0) y s n y B s 1 0.9 m � 0.8 h ( ξ ) = β k B k ( ξ ) 0.7 k =1 0.6 h( ξ ) 0.5 • Hicks-Henne bump functions 0.4 0.3 q k = log(0 . 5) B k ( ξ ) = sin p ( πξ q k ) , 0.2 log( ξ k ) 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ξ • Move points along normal to Exact derivatives d X s d β can be reference line AB computed Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 15 / 25

Grid deformation Initial grid • Interpolate displacement of surface points to interior points using RBF ˜ f ( x, y ) = a 0 + a 1 x + a 2 y + N � r j | 2 log | � b j | � r − � r − � r j | j =1 Deformed grid where � r = ( x, y ) • Results in smooth grids d X • Exact derivatives d X s can be computed 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 100 0.1 10 1 0.1 0.8 0.08 0.01 0.001 0.0001 Adjoint residual 0.06 0.00001 Cl 0.6 Residue 1 x 10 -6 Cd Cd Cl 1 x 10 -7 0.04 1 x 10 -8 Flow residual 1 x 10 -9 0.4 1 x 10 -10 0.02 1 x 10 -11 1 x 10 -12 1 x 10 -13 0.2 0 1 x 10 -14 1 x 10 -15 0 500 1000 1500 2000 2500 3000 0 500 1000 1500 2000 2500 Number of iterations Number of iterations 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 80 AD FD 60 40 Gradient 20 0 • Dot-product test -20 • Limiters can cause -40 non-differentiability 5 10 15 20 Hicks-Henne parameter • Koren limiter: dependance 150 on parameter AD FD 100 • Check adjoint derivatives 50 Gradient against finite difference 0 -50 -100 5 10 15 20 Hicks-Henne parameter Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 19 / 25

Test cases • NACA0012: M ∞ = 0 . 8, α = 1 . 25 o I = C d C l 0 C l C d 0 • RAE2822: M ∞ = 0 . 729, α = 2 . 31 o ◮ Penalty approach � � I = C d � 1 − C l � � + � � C d 0 C l 0 � ◮ Constrained minimization min I = C d s.t. C l = C l 0 C d 0 Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 20 / 25

NACA0012: Maximize L/D 0.06 1 0.04 0.5 0.02 Initial -Cp conmin_frcg 0 0 optpp_q_newton steep Initial -0.02 conmin_frcg -0.5 optpp_q_newton steep -0.04 -1 -0.06 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x/c Method 100 C d C l N fun N grad I Initial 1.000 2.072 0.295 - - conmin frcg 0.170 0.389 0.325 50 13 0.142 0.329 0.329 50 39 optpp q newton steep 0.166 0.335 0.287 50 45 Praveen. C (TIFR-CAM) Shape Optimization AeSI, Aug 2009 21 / 25