Numerical Shape Optimization Praveen. C praveen@math.tifrbng.res.in - PowerPoint PPT Presentation

Derivative of boundary functional � J (Γ , u ) = f ( u )d s Γ We perturb Γ by normal displacements Γ ǫ = { x + ǫα ( s ( x )) n : x ∈ Γ } We must consider change in area element d s due to this perturbation. First look at a line element d s with radius of curvature R , d s = R d θ Due to boundary perturbation, line element changes to � � 1 − ǫ α d s ǫ = d s R Same formula holds in 3-D but R is mean curvature defined using two orthogonal tangential coordinates R = 1 1 + 1 R 1 R 2 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 21 / 114

� Γ ǫ f ( u ǫ )d s in terms of integral and quantities on Γ First express Consider x ∈ Γ and corresponding x ǫ ∈ Γ ǫ given by x ǫ = x + ǫαn u ( x ) + ǫα∂u u ( x ǫ ) ∂n ( x ) + O ( ǫ 2 ) = u ǫ ( x ǫ ) u ( x ǫ ) + ǫ ˜ u ( x ǫ ) + O ( ǫ 2 ) = u ( x ) + ǫα∂u u ( x ) + O ( ǫ 2 ) = ∂n ( x ) + ǫ ˜ � � u ( x ) + ǫα∂u f ( u ǫ ( x ǫ )) + O ( ǫ 2 ) = f ∂n ( x ) + ǫ ˜ u ( x ) � � α∂u f ′ ( u ) + O ( ǫ 2 ) = f ( u ) + ǫ ∂n + ˜ u � � αf ′ ( u ) ∂u ∂n − α f ( u ǫ )d s ǫ u + O ( ǫ 2 ) d s + ǫf ′ ( u )˜ = f ( u )d s + ǫ Rf ( u ) Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 22 / 114

Finally the sensitivity is given by � � � � 1 f ′ ( u ) ∂u ∂n − f ( u ) f ′ ( u )˜ ǫ [ J (Γ ǫ ) − J (Γ)] = α d s + u d s + O ( ǫ ) R Γ Γ Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 23 / 114

Euler equations: compressible, inviscid flow ρ = density u = ( u 1 , u 2 , u 3 ) = velocity p = pressure Total energy per unit volume E = ρe + 1 2 ρ | u | 2 e = internal energy per unit mass For an ideal gas, p γ − 1 + 1 p 2 ρ | u | 2 e = = ⇒ E = ( γ − 1) ρ γ = C p For air γ = 1 . 4 C v Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 24 / 114

Conservation of mass Rate of change of mass inside volume Ω = – flux of mass going out through ∂ Ω � � d ρ d x = − ρ ( u · n )d s d t Ω ∂ Ω Using divergence theorem � � d ρ d x = − div( ρu )d x d t Ω Ω Since this should hold for any control volume Ω ∂ρ ∂t + div( ρu ) = 0 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 25 / 114

Conservation of momentum: Newton’s law We consider a material volume Ω( t ), i.e, Ω( t ) moves with the fluid Rate of change of momentum = Net force acting on the body � � � d ρu i d x = R i d s + ρf i d x d t Ω( t ) ∂ Ω( t ) Ω( t ) For an inviscid fluid, R = − pn . The left hand side is � � � d ∂ ρu i d x = ∂t ( ρu i )d x + ( ρu i )( u · n )d s d t Ω( t ) Ω( t ) ∂ Ω( t ) � � ∂ = ∂t ( ρu i )d x + div( ρu i u )d x Ω( t ) Ω( t ) � � � � ∂ ∂p ∂t ( ρu i )d x + div( ρu i u )d x = − d x + ρf i d x ∂x i Ω( t ) Ω( t ) Ω Ω ∂t ( ρu i ) + div( ρu i u ) + ∂p ∂ = ρf i ∂x i Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 26 / 114

Conservation of energy Rate of change of energy = rate of work done by all forces � � � d E d x = ( R · u )d s + ρ ( f · u )d x d t Ω( t ) ∂ Ω( t ) Ω( t ) Using divergence theorem � � � � ∂E ∂t d x + div( Eu )d x = − div( pu )d x + ρ ( f · u )d x Ω Ω Ω Ω( t ) ∂E ∂t + div[( E + p ) u ] = ρf · u Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 27 / 114

Euler equations • Mass conservation ∂ρ ∂ ∂ ∂ ∂t + ( ρu 1 ) + ( ρu 2 ) + ( ρu 3 ) = 0 ∂x 1 ∂x 2 ∂x 3 • Momentum conservation ∂ ∂ ∂ ∂ ( p + ρu 2 ∂t ( ρu 1 ) + 1 ) + ( ρu 1 u 2 ) + ( ρu 1 u 3 ) = 0 ∂x 1 ∂x 2 ∂x 3 ∂ ∂ ∂ ∂ ( p + ρu 2 ∂t ( ρu 2 ) + ( ρu 1 u 2 ) + 2 ) + ( ρu 2 u 3 ) = 0 ∂x 1 ∂x 2 ∂x 3 ∂ ∂ ∂ ∂ ( p + ρu 2 ∂t ( ρu 3 ) + ( ρu 1 u 3 ) + ( ρu 2 u 3 ) + 3 ) = 0 ∂x 1 ∂x 2 ∂x 3 • Energy conservation ∂E ∂ ∂ ∂ ∂t + [( E + p ) u 1 ] + [( E + p ) u 2 ] + [( E + p ) u 3 ] = 0 ∂x 1 ∂x 2 ∂x 3 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 28 / 114

Euler equations: vector conservation form ∂U ∂t + ∂F 1 + ∂F 2 + ∂F 3 = 0 ∂x 1 ∂x 2 ∂x 3 where U is the conserved variables   ρ   ρu 1     U = ρu 2     ρu 3 E and ( F 1 , F 2 , F 3 ) are the flux vectors       ρu 1 ρu 2 ρu 3  p + ρu 2      ρu 1 u 2 ρu 1 u 3  1            p + ρu 2 F 1 = ρu 1 u 2 , F 2 = , F 3 = ρu 2 u 3       2       p + ρu 2 ρu 1 u 3 ρu 2 u 3 3 ( E + p ) u 1 ( E + p ) u 2 ( E + p ) u 3 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 29 / 114

Euler equations: Hyperbolic property • Flux jacobians A = ∂F 1 B = ∂F 2 C = ∂F 3 ∂U , ∂U , ∂U • Quasi-linear form ∂U ∂t + A ∂U + B ∂U + C ∂U = 0 ∂x 1 ∂x 2 ∂x 3 ω ∈ R 3 • Hyperbolic property: A ( U, ω ) = ω 1 A + ω 2 B + ω 3 C, ◮ All eigenvalues are real ◮ Eigenvectors are linearly independent • Eigenvalues of A ( U, ω ) with | ω | = 1 are u · ω − a, u · ω, u · ω, u · ω, u · ω + a � γp a = speed of sound = ρ Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 30 / 114

Shape design using Euler equations Pressure matching problem Find the shape of the airfoil Γ such that the pressure on the airfoil is close to a specified pressure p ∗ , i.e., � 1 ( p − p ∗ ) 2 d s min 2 Γ Γ The pressure p = p ( U ) is obtained from the solution of the steady Euler equations F x + G y + H z = 0 The fluxes also satisfy the homogeneity property F ( U ) = A ( U ) U, G ( U ) = B ( U ) U, H ( U ) = C ( U ) U Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 31 / 114

Due to change in shape Γ, the solution changes from U to U + ǫ ˜ U F ( U ) + ǫ∂F F ( U + ǫ ˜ ∂U ( U ) ˜ U + O ( ǫ 2 ) U ) = F ( U ) + ǫA ˜ U + O ( ǫ 2 ) = Equation governing first order perturbation ( A ˜ U ) x + ( B ˜ U ) y + ( C ˜ U ) z = 0 For an arbitrary function V � � � V · ( A ˜ V x · A ˜ V · A ˜ U ) x = − U + Un x Ω Ω ∂ Ω � � A t V x · ˜ V · A ˜ = − U + Un x Ω ∂ Ω Then � V · [( A ˜ U ) x + ( B ˜ U ) y + ( C ˜ U ) z ] = 0 Ω Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 32 / 114

becomes � ( A t V x + B t V y + C t V z ) · ˜ − U Ω � V · ( An x + Bn y + Cn z ) ˜ + U = 0 ∂ Ω Define D = An x + Bn y + Cn z we note that   ρ ( u · n )   = normal flux DU = pn + ρu ( u · n ) ( E + p )( u · n ) and on solid wall Γ, u · n = 0   0   DU | Γ = pn 0 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 33 / 114

Since DU = Fn x + Gn y + Hn z � Fn x + ˜ ˜ Gn y + ˜ DU = Hn z A ˜ Un x + B ˜ Un y + C ˜ = Un z D ˜ = U Also d − 1 � ∂α pn = ˜ � pn + p ˜ n, ˜ n = − t j ∂t j j =1 Denoting V = [ v 1 , v, v 5 ] with v = ( v 2 , v 3 , v 4 ) � � V · D ˜ V · � U = DU Γ Γ � � = V · [0 pn � 0] = v · � pn Γ Γ � � d − 1 � ∂α = p ( v · n ) − ˜ p ( v · t j ) ∂t j Γ Γ j =1 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 34 / 114

First variation in objective function � � � ∂n − α ( p − p ∗ ) 2 1 p + α ( p − p ∗ ) ∂p ˜ ( p − p ∗ )˜ J = ǫ 2 R Γ To obtain gradient of J wrt α we must eliminate ˜ p . We add the first order perturbation equation to the above equation � 1 ˜ ( A t V x + B t V y + C t V z ) · ˜ J = RHS − U ǫ Ω � V · D ˜ + U ∂ Ω To eliminate ˜ U from the volume integral, we choose V to satisfy A t V x + B t V y + C t V z = 0 , in Ω Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 35 / 114

� � 1 ˜ V · D ˜ V · D ˜ J = RHS + U + U ǫ Γ Γ o � ( p − p ∗ + v · n )˜ = p Γ � � � ∂n − ( p − p ∗ ) 2 ( p − p ∗ ) ∂p + α 2 R Γ � d − 1 � ∂α − p ( v · t j ) ∂t j Γ j =1 � D t V · ˜ + U Γ o To eliminate ˜ p we choose the adjoint boundary conditions p − p ∗ + v · n = 0 , on Γ Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 36 / 114

We choose V on Γ o to eliminate the integral over Γ o • Supersonic inflow: U is specified so ˜ U = 0. No b.c. are imposed on V • Supersonic outflow: No b.c. for U , hence ˜ U � = 0. Choose V = 0 • Subsonic inflow: � � � D t V · ˜ D t V · T − 1 T ˜ T − t D t V · T ˜ U = U = U Γ o Γ o Γ o We have ( T ˜ U ) 1 , 2 , 3 , 4 = 0 while ( T ˜ U ) 5 is arbitrary, so we choose ( T − t D t V ) 5 = 0 • Subsonic outflow: ( T ˜ U ) 1 is specified while ( T ˜ U ) 2 , 3 , 4 , 5 is arbitrary. Hence we choose ( T − t D t V ) 2 , 3 , 4 , 5 = 0 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 37 / 114

Adjoint Euler equations A t V x + B t V y + C t V z = 0 , in Ω with boundary conditions • supersonic inflow none • supersonic outflow V = 0 • subsonic inflow ( T − t D t V ) 5 = 0 • subsonic outflow ( T − t D t V ) 2 , 3 , 4 , 5 = 0 • solid wall p − p ∗ + v · n = 0 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 38 / 114

Shape derivative � � � � d − 1 � ∂n − ( p − p ∗ ) 2 1 ( p − p ∗ ) ∂p ∂α ˜ J = α − p ( v · t j ) ǫ 2 R ∂t j Γ Γ j =1 � � � ∂n − ( p − p ∗ ) 2 ( p − p ∗ ) ∂p = α − div Γ ( pv ) 2 R Γ Hence ∂n − ( p − p ∗ ) 2 ∇ α J = ( p − p ∗ ) ∂p − div Γ ( pv ) 2 R Steepest descent update α n +1 = α n − s n ∇ α J ( α n ) Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 39 / 114

Discrete adjoint approach • Approximate PDE using finite volume method R ( X, Q ) = 0 where X ∈ R m X = grid coordinates , Q ∈ R n Q = solution vector , R : R m × R n → R n • Q depends implicitly on X through R ( X, Q ) = 0 • Discrete approximation to objective function I ( X, Q ) Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 40 / 114

Elements of practical shape optimization 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 IISc, Dec 2009 41 / 114

d I Adjoint approach: d X • For shape optimization: I = I ( X, Q ) d X = ∂ I d I ∂X + ∂ I ∂Q ∂Q ∂X • Differentiate state equation R ( X, Q ) = 0 ∂X + ∂R ∂R ∂Q ∂X = 0 ∂Q • Flow sensitivity ∂Q ∂X ; costly to evaluate • Introducing an adjoint variable Ψ ∈ R n , we can write � ∂ I � � ∂R � d I ∂X + ∂ I ∂Q ∂X + ∂R ∂Q + Ψ t d X = ∂Q ∂X ∂Q ∂X Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 42 / 114

Adjoint approach • Collect terms involving the flow sensitivity � ∂ I � � ∂ I � ∂Q d I ∂X + Ψ t ∂R ∂Q + Ψ t ∂R d X = + ∂X ∂Q ∂X • Choose Ψ so that flow sensitivity vanishes � ∂R � t � ∂ I � t ∂Q + Ψ t ∂R ∂ I ∂Q = 0 = ⇒ Ψ + = 0 ∂Q ∂Q • Gradient d X = ∂ I d I ∂X + Ψ t ∂R ∂X Cost of Adjoint-based steepest-descent Cost = O (1) N iter = O ( N ) Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 43 / 114

Optimization steps • β = ⇒ X s = ⇒ X • Solve the flow (primal) equations to steady-state d Q X = ⇒ d t + R ( X, Q ) = 0 = ⇒ Q, I • Solve adjoint equations to steady-state � ∂R � t � ∂ I � t dΨ X, Q = ⇒ d t + Ψ + = 0 = ⇒ Ψ ∂Q ∂Q • Compute gradient wrt grid X d X = ∂ I d I ∂X + Ψ t ∂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 IISc, Dec 2009 44 / 114

Automatic Differentiation • Discrete adjoint method requires � ∂R � t � ∂R � t ∂J ∂J ∂X , ∂Q, Ψ , Ψ ∂X ∂Q • Differentiate J and R by hand and write new code to compute above derivatives • J and R are already available as a computer code ◮ elementary operations: ∗ , /, + , − ◮ elementary functions: sin, exp, log, etc. ◮ chain rule of differentiation ◮ use Automatic Differentiation • AD tool: automates the application of chain rule of differentiation X Code for Code for AD tool ∂J J ( X, Q ) ∂X Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 45 / 114

Automatic Differentiation: forward and backward mode X ∈ R m , Y ∈ R n , Y = F ( X ) • Forward mode � ∂F � X ∈ R m , X ∈ R m ˙ X ∈ R n ˙ → ∂X • Backward mode � ∂F � t X ∈ R m , Y ∈ R n ¯ Y ∈ R m ¯ → ∂X Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 46 / 114

Automatic Differentiation: TAPENADE Code residue.f to compute X, Q → R ( X, Q ) = RES REAL X(M), Q(N), RES(N) CALL RESIDUE(X, Q, RES) Differentiate RES wrt Q tapenade -backward -vars Q -outvars RES residue.f � � t ∂R Generates new code residue b.f : X, Q, Ψ → Ψ ∂Q REAL X(M), Q(N), RES(N) REAL QB(N), RESB(N) CALL RESIDUE_B(X, Q, QB, RES, RESB) Ψ SUBROUTINE RESIDUE_B(X, Q, QB, RES, RESB) ⎛ ⎞ T ∂ R Ψ ⎜ ⎟ ⎝ ∂ Q ⎠ Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 47 / 114

Reverse mode Consider a subroutine which takes in a vector x ( n ) and returns a scalar function f = f ( x ). A computer program to compute f contains a sequence of computations. Each line takes some set of previously computed values and returns a new value. t 1 = L 1 ( T 0 ) x ∈ R n , T 0 = x t 2 = L 2 ( T 1 ) . . t r ∈ R . . . = . � T r − 1 � t p = L p ( T p − 1 ) T r = t r f = L p +1 ( T p ) The function f as coded, in general, depends on all of these variables x, t 1 , t 2 , . . . , t p Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 48 / 114

AD tools • Source transformation and operator overloading • See http://www.autodiff.org Tool Modes Method Lang ADIFOR F ST Fortran ADOLC F/B OO C TAPENADE F/B ST Fortran/C TAF/TAMC F/B ST Fortran/C MAD F OO Matlab Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 49 / 114

Important criteria for MDO applications of complex, 3D models Consistent: Is the parameterization consistent across multiple disciplines? Airplane shape design variables: Are the design variables directly related to the airplane shape design variables such as camber, thickness, twist, shear, and planform? Compact: Does the parameterization provide a compact set of design variables? 10s vs 1000s Smooth: Does the shape perturbation maintain a smooth geometry? Local control: Is there any local control on shape changes? Analytical sensitivity: Is it feasible to calculate the sensitivity analytically? Grid deformation: Does the parameterization allow the grid to be deformed? Setup time: Can a shape optimization application be set up quickly? hours, days, weeks, months? Existing grid: Does the parameterization allow the existing grid to be reused? Does it require to reverse engineer the baseline design parameters? CAD: Is there a direct connection to the CAD system? Samareh: Survey of shape parameterization techniques Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 50 / 114

Examples of shape optimization Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 51 / 114

Quasi 1-D flow h ( x ) Inflow Outflow x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 52 / 114

Quasi 1-D flow • Quasi 1-D flow in a duct ∂t ( hU ) + ∂ ∂ ∂x ( hf ) = d h d xP, x ∈ ( a, b ) t > 0       ρ ρu 0   ,   ,   p + ρu 2 U = ρu f = P = p E ( E + p ) u 0 h ( x ) = cross-section height of duct • Inverse design: find shape h to get pressure distribution p ∗ • Optimization problem: find the shape h which minimizes � b ( p − p ∗ ) 2 d x J = a Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 53 / 114

Quasi 1-D flow Inflow face 1 i N 2 / 1 2 − / 1 i + i Outflow face Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 54 / 114

Quasi 1-D flow • Finite volume scheme d t + h i +1 / 2 F i +1 / 2 − h i − 1 / 2 F i − 1 / 2 = ( h i +1 / 2 − h i − 1 / 2 ) d U i h i P i ∆ x ∆ x • Discrete cost function N � ( p i − p ∗ i ) 2 J = i =1 • Control variables h 1 / 2 , h 1+1 / 2 , . . . , h i +1 / 2 , . . . , h N +1 / 2 • N = 100 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 55 / 114

Duct shape Duct shape 1.9 Target shape Current shape 1.8 1.7 1.6 1.5 h 1.4 1.3 1.2 1.1 1 0 2 4 6 8 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 56 / 114

Target pressure distribution p ∗ Target pressure 2.5 AUSMDV KFVS LF 2 Pressure 1.5 1 0.5 0 2 4 6 8 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 57 / 114

Current pressure distribution Starting pressure 2.5 AUSMDV KFVS LF 2 Pressure 1.5 1 0.5 0 2 4 6 8 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 58 / 114

Adjoint density 5 AUSMDV KFVS 0 LF −5 −10 Adjoint density −15 −20 −25 −30 −35 0 2 4 6 8 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 59 / 114

Convergence history: Explicit Euler Convergence history with AUSMDV flux 0 Flow Adjoint −1 −2 −3 −4 Residue −5 −6 −7 −8 −9 −10 1000 2000 3000 4000 5000 6000 7000 No of iterations Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 60 / 114

Shape gradient Gradient 35 AUSMDV KFVS LF 30 25 20 Gradient 15 10 5 0 0 2 4 6 8 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 61 / 114

Shape gradient Gradient 40 AUSMDV KFVS 35 LF 30 25 Gradient 20 15 10 5 0 −5 2 2.5 3 3.5 4 4.5 5 5.5 6 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 62 / 114

Validation of Shape gradient Gradient with AUSMDV flux 35 30 B 25 20 Gradient 15 10 C A 5 0 0 2 4 6 8 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 63 / 114

Validation of shape gradient ∂J ∂h ≈ J ( h + ∆ h ) − J ( h − ∆ h ) 2∆ h ∆ h A B C 0.01 0.4191069499 35.18452823 2.545316345 0.001 0.4231223499 36.10982621 2.556461900 0.0001 0.4231624999 36.11933154 2.556573499 0.00001 0.4231599998 36.11942125 2.556575000 0.000001 0.4231999994 36.11942305 2.556550001 0.0000001 0.4229999817 36.11942329 2.556499971 AD 0.4231628330 36.11941951 2.556574450 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 64 / 114

Gradient smoothing • Non-smooth gradients G especially in the presence of shocks • Smooth using an elliptic equation � � 1 − ǫ ( x ) d 2 ¯ G ( x ) = G ( x ) d x 2 ǫ i = {| G i +1 − G i | + | G i − G i − 1 |} L i | G i +1 − 2 G i + G i − 1 | L i = max( | G i +1 − G i | + | G i − G i − 1 | , TOL) • Finite difference with Jacobi iterations Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 65 / 114

Gradient smoothing Gradient using AUSMDV flux 35 Original gradient Smoothed gradient 30 25 20 Gradient 15 10 5 0 0 2 4 6 8 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 66 / 114

Quasi 1-D optimization: Shape 1.9 Target Initial 1.8 1.7 1.6 1.5 h 1.4 1.3 1.2 1.1 1 0 1 2 3 4 5 6 7 8 9 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 67 / 114

Quasi 1-D optimization: Final shape 1.9 1.8 1.7 1.6 1.5 h 1.4 Target Initial 1.3 Final 1.2 1.1 1 0 1 2 3 4 5 6 7 8 9 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 68 / 114

Quasi 1-D optimization: Pressure 2.5 2 Target Pressure Initial 1.5 Final 1 0.5 0 1 2 3 4 5 6 7 8 9 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 69 / 114

Quasi 1-D optimization: Convergence 0 10 Cost Gradient −1 10 −2 10 Cost/Gradient −3 10 −4 10 −5 10 −6 10 0 5 10 15 20 25 30 35 40 45 50 No. of iterations Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 70 / 114

Quasi 1-D optimization: Adjoint density 0 −5 Initial Final Adjoint density −10 −15 −20 −25 0 2 4 6 8 10 x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 71 / 114

Airfoil shape optimization Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 72 / 114

Shape parameterization • Parameterize the n � deformations � � � x s � � n x � x (0) s A B = + h ( ξ ) ξ y (0) y s n y 1 s 0.9 0.8 m � 0.7 h ( ξ ) = β k B k ( ξ ) 0.6 k =1 h( ξ ) 0.5 0.4 • Hicks-Henne bump functions 0.3 0.2 q k = log(0 . 5) B k ( ξ ) = sin p ( πξ q k ) , 0.1 log( ξ k ) 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ξ Exact derivatives d X s d β can be • Move points along normal to computed reference line AB Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 73 / 114

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 IISc, Dec 2009 74 / 114

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 Hemker-Koren limiter • Implicit scheme Source code of NUWTUN available online http://nuwtun.berlios.de Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 75 / 114

Validation of adjoint gradients 80 AD • Dot-product test FD 60 ˙ R ′ ( Q ) ˙ 40 R := Q Gradient 20 R ( Q + ǫ ˙ Q ) − R ( Q ) 0 ≈ ǫ -20 [ R ′ ( Q )] ⊤ ¯ ¯ Q := R -40 R ⊤ ¯ ? Q ⊤ ¯ 5 10 15 20 ˙ ˙ R = Q Hicks-Henne parameter 150 • Upwind schemes and limiters AD FD 100 can cause problems due to 50 non-differentiability Gradient 0 • Check adjoint derivatives -50 against finite difference -100 ◮ NACA0012: C d /C l ◮ RAE2822: C d 5 10 15 20 Hicks-Henne parameter Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 76 / 114

Convergence and limiter: RAE2822 airfoil L ( a, b ) = a ( b 2 + 2 ǫ ) + b (2 a 2 + ǫ ) , 0 < ǫ ≪ 1 2 a 2 − ab + 2 b 2 + 3 ǫ ǫ = 10 − 8 ǫ = 10 − 4 1 1 Density Density x-momen x-momen y-momen 0.1 y-momen z-momen z-momen Energy Energy 0.1 0.01 0.001 0.01 0.0001 Residue Residue 1e-05 0.001 1e-06 1e-07 0.0001 1e-08 1e-09 1e-05 1e-10 0 200 400 600 800 1000 1200 1400 1600 0 200 400 600 800 1000 1200 1400 1600 Number of iterations Number of iterations Adjoint iterations blow-up No blow-up Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 77 / 114

Iterative convergence tests Primal residual Adjoint residual R n = || R ( Q n ) || R n = || [ R ′ ( Q ∞ )] ⊤ Ψ n + I ′ ( Q ∞ ) || 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 Cl Cd 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 IISc, Dec 2009 78 / 114

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 IISc, Dec 2009 79 / 114

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 I 100 C d C l N fun N grad Initial 1.000 2.072 0.295 - - 0.170 0.389 0.325 50 13 conmin frcg optpp q newton 0.142 0.329 0.329 50 39 0.166 0.335 0.287 50 45 steep Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 80 / 114

RAE2822: Drag minimization, penalty approach 1.5 0.06 1 0.04 0.02 0.5 Initial -Cp conmin_frcg y 0 optpp_q_newton 0 steep -0.02 Initial -0.5 conmin_frcg optpp_q_newton -0.04 steep -1 -0.06 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x Method I 100 C d C l N fun N grad Initial 1.000 1.150 0.887 - - 0.355 0.405 0.890 50 13 conmin frcg optpp q newton 0.351 0.400 0.884 50 51 0.341 0.388 0.884 50 47 steep Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 81 / 114

RAE2822: Lift-constrained drag minimization 1.5 0.06 1 0.04 0.5 0.02 Initial conmin_mfd -Cp fsqp 0 ipopt 0 Initial -0.02 conmin_mfd -0.5 fsqp ipopt -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 I 100 C d C l N fun N grad Initial 1.000 1.150 0.887 - - 0.342 0.394 0.881 50 22 conmin mfd fsqp 0.325 0.374 0.887 50 32 0.335 0.385 0.887 45 62 ipopt Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 82 / 114

RAE2822: Lift-constrained drag minimization 1.5 0.06 1 0.04 0.5 0.02 ipopt Initial -Cp 0 0 ipopt -0.02 Initial -0.5 -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 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 83 / 114

RANS computation M=0.729, Re=6.5 million, Cl=0.88, k-w turbulence model 1.5 RAE2822 Optimized 1 0.5 -Cp 0 -0.5 -1 0 0.2 0.4 0.6 0.8 1 x/c Shape Inviscid RANS AOA C l 100 C d AOA C l 100 C d RAE2822 2.31 0.887 1.150 3.28 0.887 2.558 Optimized 2.31 0.887 0.386 3.22 0.885 1.692 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 84 / 114

Sensitivity to perturbations 0.04 Initial 200 Initial Optimized Optimized 0.03 150 Drag coefficient Lift/Drag 0.02 100 0.01 50 0 0 0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.7 0.71 0.72 0.73 0.74 0.75 0.76 Mach number Mach number (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 IISc, Dec 2009 85 / 114

Adjoints with shocks: Scalar conservation law Consider the functional � +1 J = G ( u ( x, T ))d x − 1 where u is the solution of the conservation law u t + f ( u ) x = 0 , x ∈ (0 , 1) , t ∈ (0 , T ) u ( x, 0) = u 0 ( x ) To compute the derivative of J wrt u 0 we make a small perturbation u 0 − → u 0 + ˜ u 0 . Then solution changes to u + ˜ u and � g ( u ) = d G ˜ J = g ( u ( x, T ))˜ u d x, d u The equation for ˜ u is obtained by linearizing the conservation law u ) t + ( f ′ ( u )˜ (˜ u ) x = 0 ˜ u ( x, 0) = u 0 ˜ Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 86 / 114

To write ˜ J in terms of ˜ u 0 , we introduce the adjoint variable v ( x, t ) � � T � � � ˜ u ) t + ( f ′ ( u )˜ J = g ( u ( x, T ))˜ u ( x, T )d x − v (˜ u ) x d x d t 0 Integrating by parts, we get � � ˜ J = v ( x, 0)˜ u 0 ( x )d x + ( g − v )˜ u ( x, T )d x � T � � � v t + f ′ ( u ) v x − u ( x, t )d x d t ˜ 0 To eliminate ˜ u ( x, t ) and ˜ u ( x, T ) from above equation, choose v t + f ′ ( u ) v x = 0 v ( x, T ) = g ( u ( x, T )) The derivative is given by � ˜ J = v ( x, 0)˜ u 0 ( x )d x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 87 / 114

Note that the adjoint equation must be solved backward in time, starting from t = T to t = 0. Example: Inviscid Burger’s equation u t + ( u 2 / 2) x = 0 , x ∈ R , t ∈ (0 , T ) � u l x < 0 u ( x, 0) = u r x > 0 Exact solution: � u l x < St S = 1 u ( x, t ) = 2( u l + u r ) u r x > St Consider the objective functional � J = 1 u 2 ( x, T )d x 2 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 88 / 114

Then the corresponding adjoint equation is v t + uv x = 0 � u l x < ST v ( x, T ) = u ( x, T ) = u r x > ST t t = T x Adjoint solution cannot be determined in the shaded region . In other regions, solution can be determined by going forward in time along the characteristic until you hit t = T line. Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 89 / 114

Back to scalar conservation law Conservation laws can have discontinuous solutions, but the discontinuity must satisfy the Rankine-Hugoniot jump conditions. If there is a discontinuity at x = x s , then it must satisfy f ( u ( x + s , t )) − f ( u ( x − s , t )) = S · ( u ( x + s , t ) − u ( x − s , t )) where S = d x s d t is the shock speed. The linearization of the conservation law must also include the linearization of the jump conditions. d x s d t → d x s d t + d˜ x s x s → x s + ˜ x s , d t � ∂u � � d f � � d f � d˜ d t [ u ] + d x s x s u ] + d x s ∂u d t [˜ d t ˜ x s = d u ˜ u + ˜ x s ∂x d u ∂x Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 90 / 114

� x s ( T ) � +1 J = G ( u ( x, T ))d x + G ( u ( x, T ))d x − 1 x s ( T ) Perturbation in J is x s ( T ) � � +1 ˜ J = g ( u ( x, T ))˜ u ( x, T )d x + g ( u ( x, T ))˜ u ( x, T )d x − [ G ] x s ( T ) ˜ x s ( T ) − 1 x s ( T ) x s ( t ) � T � � T � +1 � � � � u t + ( f ′ ( u )˜ u t + ( f ′ ( u )˜ − v ˜ u ) x d x d t − v ˜ u ) x d x d t 0 − 1 0 x s ( t ) � d˜ � ∂u � � d f � � d f �� � T d t [ u ] + d x s x s u ] + d x s ∂u − v s ( t ) d t [˜ d t ˜ x s − d u ˜ u − ˜ x s d t ∂x d u ∂x 0 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 91 / 114

After some work we get � +1 ˜ J = v ( x, 0)˜ u 0 ( x )d x − 1 provided the adjoint variables v ( x, t ) and v s ( t ) satisfy the following equations v t + f ′ ( u ) v x = 0 v ( x, T ) = g ( u ( x, T )) v ( x + s ( t ) , t ) = v ( x − v s ( t ) = s ( t ) , t ) d d tv s ( t ) = 0 [ G ] v s ( T ) = [ u ] Adjoint variable requires an interior boundary condition along the shock ( x s ( t ) , t ). This uniquely determines the adjoint solution. Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 92 / 114

Steepest descent: limitations, problems • Multiple optima are common • Convergence towards local optimum: dependance on starting guess S o • Shape gradient: difficult to compute, or impossible • Noisy cost function J Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 93 / 114

Gradient-free methods • Simplex method (Nelder-Mead, Torczon) • Global search methods ◮ Genetic algorithm ◮ Particle swarm optimization ◮ Ant colony optimization ◮ ... • Do not require gradient • Collection of M solutions at any iteration n P n = { x n 1 , x n 2 , . . . , x n M } ⊂ R d • Solutions evolve according to some rules P n +1 = E ( P n ) Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 94 / 114

Particle swarm optimization • Kennedy and Eberhart (1995) • Modeled on behaviour of animal swarms: ants, bees, birds • Swarm intelligence Optimization problem D ⊂ R d min x ∈ D J ( x ) , Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 95 / 114

Particle swarm optimization • Kennedy and Eberhart (1995) • Modeled on behaviour of animal swarms: ants, bees, birds • Swarm intelligence Optimization problem D ⊂ R d min x ∈ D J ( x ) , Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 95 / 114

Particle swarm optimization Particles distributed in design space x i ∈ D, i = 1 , ..., N p X2 X1 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 96 / 114

Particle swarm optimization Each particle has a velocity v i ∈ R d , i = 1 , ..., N p X2 X1 Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 97 / 114

Particle swarm optimization • Particles have memory ( t = iteration number) Local memory : p t J ( x s i = argmin i ) 0 ≤ s ≤ t Global memory : p t = argmin J ( p t i ) i • Velocity update 2 ( p t − x t v t +1 = ωv t i + c 1 r t 1 ( p t i − x t + c 2 r t i ) i ) i � �� � � �� � Local Global • Position update x t +1 = x t i + v t +1 i i Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 98 / 114

Basic PSO algorithm !"# $%&!&'(&)*+ ,-.&!&-%/+0*(-1&!2 3-4,5!*+1-.!+ 65%1!&-% 7,8'!*+(-1'(+'%8+ 9(-:'(+4*4-;2 7,8'!*+0*(-1&!2+ '%8+,-.&!&-% 3-%0*;9*%1*+< ?- !"!=> A*. @!-, Praveen. C (TIFR-CAM) Shape Optimization IISc, Dec 2009 99 / 114