Numerical shape optimization for compressible flows Praveen. C - PowerPoint PPT Presentation

Numerical shape optimization for compressible flows Praveen. C praveen@math.tifrbng.res.in Tata Institute of Fundamental Research Center for Applicable Mathematics Bangalore 560065 http://math.tifrbng.res.in ICMPDE TIFR-CAM, Bangalore 13–17 August, 2010 Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 1 / 63

Shape optimization framework • Shape is parameterized in terms of x ∈ D ⊂ R d • PDE-constrained minimization min x ∈ D J ( x, u ) s.t. R ( x, u ) = 0 • Solving R = 0 is computationally expensive • d can be large - Curse of dimensionality Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 2 / 63

Navier-Stokes Equations I • In d-dimensions d d ∂U ∂F i ∂G i � � ∂t + = ∂x i ∂x i i =1 i =1 • Conserved quantities and fluxes       ρ ρu i 0 ρu 1 pδ i 1 + ρu 1 u i τ i 1             U = ρu 2 , F i = pδ i 2 + ρu 2 u i , G i = τ i 2             ρu 3 pδ i 3 + ρu 3 u i τ i 3       E ( E + p ) u i τ ij u j − q i ρ = Density ( u 1 , u 2 , u 3 ) = Velocity p = Pressure E = Energy per unit volume τ ij = Viscous stress tensor q i = Heat flux Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 3 / 63

Navier-Stokes Equations II • Ideal gas equation of state � � E − 1 2 ρ | u | 2 p = ( γ − 1) • Constitutive law � ∂u i � + ∂u j − 2 τ ij = ( µ + µ t ) 3( ∇ · u ) δ ij ∂x j ∂x i � µ � ∂T Pr + µ t q i = − Pr t ∂x i • Additional equations, Turbulence models , to determine µ t Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 4 / 63

Quantities of interest • Forces on a solid body � F i = ( − pn i + τ ij n j )d S S • Lift and drag L = F · V ⊥ D = F · V ∞ ∞ • Optimization problem min D s.t. L = W, etc. Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 5 / 63

Finite volume method ∂U ∂t + ∇ · H ( U ) = 0 • Divide Ω into non-overlapping, polygonal, finite volumes Ω i Ω = ∪ i Ω i • Conservation principle on each finite volume � � d U d x + H · n d s = 0 d t Ω i ∂ Ω i • Semi-discrete scheme | Ω i | d U i � d t + H ( U i , U j , n ij ) = 0 j ∈ N ( i ) Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 6 / 63

Classical approach • PDE-constrained problem min x J ( u, x ) s.t. R ( u, x ) = 0 or L ( x, u, v ) = J ( x, u ) + ( v, R ( u, x )) • Need to develop complex adjoint solvers • PDE models not well motivated, e.g., turbulence models min x J h ( u h , x ) s.t. R h ( u h , x ) = 0 • Discrete approach: Automatic Differentiation • Noisy objective functions • Only local optimum Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 7 / 63

Discrete adjoint in presence of shocks 0.8 Automatic Differentiation 0.11 Divided Differences 0.7 0.1 0.6 0.09 0.08 0.5 0.07 d c d /d Mach drag coefficient 0.4 0.06 0.3 0.05 0.2 0.04 0.03 0.1 0.02 0 0.01 -0.1 0 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Mach-0.83 Mach-0.83 Figure 7: Drag derivatives with respect to Mach number in Figure 6: Drag with respect to Mach number in transonic transonic regime. regime. (Martinelli et al.) Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 8 / 63

Discrete adjoint in presence of shocks 0.045 Naca 256x64 0.038 Naca 128x32 0.04 Naca 64x16 0.036 0.035 0.034 C-lift C-lift 0.03 0.032 0.03 0.025 0.028 0.02 0.026 0.7 0.75 0.8 0.85 0.8 0.81 0.82 0.83 0.84 0.85 Mach Mach 3 3 2 2 d C-lift / d Mach d C-lift / d Mach 1 1 0 0 -1 -1 -2 -2 0.7 0.75 0.8 0.85 0.8 0.81 0.82 0.83 0.84 0.85 Mach Mach Figure 1. Lift coefficient and derivatives against Mach number for a transonic NACA0012 aerofoil on three grids. The bars in the top plots display the derivatives computed using the discrete adjoint. (Dwight et al.) Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 9 / 63

Discrete adjoint: Frozen turbulence 0.2 C-lift 0.14 1 C-lift C-drag C-drag 0.18 Grad - exact Grad - frozen turb 1 0.12 0.16 0.14 0.1 0.5 0.12 0.5 C-drag C-drag 0.08 C-lift C-lift 0.1 0.06 0.08 0 0.06 0 0.04 0.04 0.02 0.02 0 -4 -2 0 2 4 6 8 10 0 5 10 Alpha Alpha Exact adjoint Frozen turb. 0.04 1 d C-drag / d alpha 0.02 C-lift 0.8 0 0.6 -0.02 -0.04 0.4 -4 -2 0 2 4 6 8 10 1 2 3 4 5 6 Alpha Alpha Figure 2. Lift and drag against angle-of-attack for an RAE 2822 aerofoil. Line segments represent gradients computed with a discrete adjoint code, with a full linearization of the turbulence model (black), and a frozen- turbulence approximation (red). (Dwight et al.) Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 10 / 63

Derivative-free methods • Using function values only • Global, stochastic search ◮ Genetic algorithm ◮ Particle swarm method D ⊂ R d min x ∈ D J ( x ) , • Collection of N p solutions at any iteration n P n = { x n 1 , x n 2 , . . . , x n N p } ⊂ D • Solutions evolve according to some rules P n +1 = E ( P n ) Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 11 / 63

Particle swarm optimization • Kennedy and Eberhart (1995) • Modeled on behaviour of animal swarms: ants, bees, birds • Cooperative behaviour of large number of individuals through simple rules • Emergence of swarm intelligence Optimization problem D ⊂ R 2 min x ∈ D J ( x ) , Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 12 / 63

Particle swarm optimization • Kennedy and Eberhart (1995) • Modeled on behaviour of animal swarms: ants, bees, birds • Cooperative behaviour of large number of individuals through simple rules • Emergence of swarm intelligence Optimization problem D ⊂ R 2 min x ∈ D J ( x ) , Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 12 / 63

Particle swarm optimization Particles distributed in design space x i ∈ D, i = 1 , ..., N p X2 X1 Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 13 / 63

Particle swarm optimization Each particle has a velocity v i ∈ R d , i = 1 , ..., N p X2 X1 Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 14 / 63

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 TIFR-CAM, Aug 2010 15 / 63

Surrogate Models Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 16 / 63

Metamodels • Replace min x ∈ D J ( x, u ) s.t. R ( x, u ) = 0 with ˜ min J ( x ) x ∈ D • Model error estimate/indicator σ ( x ) x ∈ D J ρ ( x ) := ˜ min J ( x ) − ρσ ( x ) , ρ ≥ 0 • Local and global search x 0 = argmin J 0 ( x ) x 3 = argmin J 3 ( x ) x ∈ D x ∈ D ⇒ Update ˜ J ( x 0 ) , J ( x 3 ) = J Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 17 / 63

Kriging I Unknown function f : D ⊂ R d → R Given the data as F N = { f 1 , f 2 , . . . , f N } ⊂ R sampled at X N = { x 1 , x 2 , . . . , x N } ⊂ D , infer the function value at a new point x N +1 ∈ D . Treat result of a computer simulation as a fictional gaussian process F N is assumed to be one sample of a multivariate Gaussian process with joint probability density � � − 1 N C − 1 2 F ⊤ p ( F N ) = exp N F N (1) � (2 π ) N det( C N ) where C N is the N × N covariance matrix. Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 18 / 63

Kriging II When adding a new point x N +1 , the resulting vector of function values F N +1 is assumed to be a realization of the ( N + 1)-variable Gaussian process with joint probability density � � − 1 N +1 C − 1 2 F ⊤ p ( F N +1 ) = exp N +1 F N +1 (2) � (2 π ) N +1 det( C N +1 ) Using Baye’s rule we can write the probability density for the unknown function value f N +1 , given the data ( X N , F N ) as � � − ( f N +1 − ˆ f N +1 ) 2 p ( f N +1 | F N ) = p ( F N +1 ) = 1 Z exp 2 σ 2 p ( F N ) f N +1 where ˆ f N +1 = k ⊤ C − 1 σ 2 f N +1 = κ − k ⊤ C − 1 N F N , N k (3) � �� � � �� � Inference Error indicator Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 19 / 63

Kriging III Covariance matrix: Given in terms of a correlation function, C N = [ C mn ], C mn = corr( f m , f n ) = c ( x m , x n ) � � d ( x i − y i ) 2 − 1 � c ( x, y ) = θ 1 exp + θ 2 r i 2 2 i =1 Parameters Θ = ( θ 1 , θ 2 , r 1 , r 2 , . . . , r d ) determined to maximize the likelihood of known data max log( p ( F N )) Θ Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 20 / 63

Kriging: Illustration 12 10 DACE predictor 8 6 12 4 standard error 10 2 of the predictor 8 0 0 2 4 6 8 10 12 6 4 2 0 0 2 4 6 8 10 12 1 1.5 1.5 2 2.5 2.5 3 3.5 3.5 4 4.5 Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 21 / 63

Minimization of 2-D Branin function: Initial database 15 10 5 0 −5 0 5 10 Praveen. C (TIFR-CAM) Shape optimization TIFR-CAM, Aug 2010 22 / 63