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

SLIDE 1

Numerical Shape Optimization

Praveen. C

praveen@math.tifrbng.res.in

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

IMI Workshop on Control and Inverse Problems IISc, Bangalore 1-15 December, 2009

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 1 / 114

SLIDE 2

Aerodynamic shape optimization

http://www.aviation-history.com

Lift = Weight Drag = Thrust min

β

Drag(β), such that Lift(β) = W

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 2 / 114

SLIDE 3

Effect of shape on flow

http://www.aerospaceweb.org

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 3 / 114

SLIDE 4

Wing and Airfoils

http://www.centennialofflight.gov

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 4 / 114

SLIDE 5

Compressible flow and shocks

Mach number, M = speed of air speed of sound Range, R = M a cT CL CD log Wi Wf

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 5 / 114

SLIDE 6

Objectives and controls

Objective function I(β)

mathematical representation of system performance I(β) = I(β, Q(β))

Control variables β

β represents the shape We want to minimize/maximize I wrt β

State variable Q: solution of an ODE or PDE

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

Basic problem

How to find the derivative of the objective function I wrt the design variables (shape) β ?

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 6 / 114

SLIDE 7

Gradient-based minimization: blackbox/FD approach

min

β∈RN I(β, Q(β))

Initialize β0, n = 0
For n = 0, . . . , Niter

◮ Solve R(βn, Qn) = 0 ◮ For j = 1, . . . , N ⋆ βn (j) = [βn 1 , . . . , βn j + ∆βj, . . . , βn N]⊤ ⋆ Solve R(βn (j), Qn (j)) = 0 ⋆ dI dβj ≈ I(βn

(j),Qn (j))−I(βn,Qn)

∆βj ◮ Steepest descent step

βn+1 = βn − sn dI dβ (βn)

Cost of FD-based steepest-descent

Cost = O(N + 1)Niter = O(N + 1)O(N) = O(N2)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 7 / 114

SLIDE 8

Accuracy of FD: Choice of step size

d dxf(x0) = f(x0 + δ) − f(x0) δ + O(δ) In principle, if we choose a small δ, we reduce the error. But computers have finite precision. Instead of f(x0) the computers gives f(x0) + O(ǫ) where ǫ = machine precision. [f(x0 + δ) + O(ǫ)] − [f(x0) + O(ǫ)] δ = f(x0 + δ) − f(x0) δ + C1 ǫ δ = d dxf(x0) + C2δ + C1 ǫ δ

Total error

, C1, C2 depend on f, x0, precision Total error is least when δ = δopt = C3 √ǫ, C3 =

C1/C2
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 8 / 114

SLIDE 9

Accuracy of FD: Choice of step size

Error δ δopt C1δ C2ǫ/δ

Precision ǫ δopt Single 10−8 10−4 Double 10−16 10−8 See: Brian J. McCartin, Seven Deadly Sins of Numerical Computation The American Mathematical Monthly, Vol. 105, No. 10 (Dec., 1998), pp. 929-941

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 9 / 114

SLIDE 10

Drag gradient using FD (Samareh)

Errors in Finite Difference Approximations

Scaled Step Size % Error

10

9

10

8

10

7

10

6

10

5

10

4

10

3

10

2

10

1

10

6

10

5

10

4

10

3

10

2

10

1

10 10

1

Mid chord LE Sweep2 Tip chord LE sweep3 Twist (root) Twist (mid) Twist (tip)

Tip chord Mid chord LE Sweep2 LE Sweep3 Twist (root) Twist (mid) Twist (tip)

Roundoff Truncation

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 10 / 114

SLIDE 11

Iterative problems

I(β, Q), where R(β, Q) = 0

Q is implicitly defined, require an iterative solution method.
Assume a Q0 and iterate Qn −

→ Qn+1 until ||R(β, Qn)|| ≤ TOL

If TOL is too small, need too many iterations
Many problems, we cannot reduce ||R(β, Qn)|| to small values
This means that numerical value of I will be noisy

Finite difference will contain too much error, and is useless RAE5243 airfoil, Mach=0.68, Re=19 million, AOA=2.5 deg. iter Lift Drag 41496 0.824485788042416 1.627593747613790E-002 41497 0.824485782714867 1.627593516695762E-002 41498 0.824485777387834 1.627593285794193E-002 41499 0.824485772061306 1.627593054909022E-002 41500 0.824485766735297 1.627592824040324E-002

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 11 / 114

SLIDE 12

Complex variable method

f(x0 + iδ) = f(x0) + iδf′(x0) + O(δ2) + iO(δ3) f′(x0) = 1 δ imag[f(x0 + iδ)] + O(δ2)

No roundoff error
We can take δ to be very small, δ = 10−20
Can be easily implemented

◮ fortran: redeclare real variables as complex ◮ matlab: no change

Iterative problems: β −

→ β + i∆β

◮ Obtain ˜

Q = Q(β + i∆β) by solving R(β + i∆β, ˜ Q) = 0

◮ Then gradient

I′(β) ≈ 1 ∆β imag[I(β + i∆β, Q(β + i∆β)]

Computational cost is O(N2) or higher (due to complex

arithmetic)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 12 / 114

SLIDE 13

Sensitivity equation approach

PDE constrained optimization problem

min

β

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

Derivative of J(β) = J(β, Q(β))

J′(β)˜ β = lim

ǫ→0

J(β + ǫ˜ β) − J(β) ǫ = Jβ ˜ β + JQ ˜ Q

Differentiating the state equation

Rβ ˜ β + RQ ˜ Q = 0

Suppose β ∈ RN. To compute dJ

dβi

◮ Set ˜

β(i) = [0, . . . , 0, 1, 0, . . . , 0]

◮ Solve for ˜

Q(i) from Rβ ˜ β(i) + RQ ˜ Q(i) = 0

◮ Compute

dJ dβi = Jβ ˜

β(i) + JQ ˜ Q(i)

◮ No problem of roundoff

error

◮ Cost of gradient is O(N 2)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 13 / 114

SLIDE 14

Adjoint equation approach

J′(β)˜ β = Jβ ˜ β + JQ ˜ Q, Rβ ˜ β + RQ ˜ Q = 0 Introduce the adjoint variable Ψ J′(β)˜ β = Jβ ˜ β + JQ ˜ Q + (Rβ ˜ β + RQ ˜ Q, Ψ) = (Jβ ˜ β, 1) + (JQ ˜ Q, 1) + (Rβ ˜ β + RQ ˜ Q, Ψ) = (˜ β, J∗

β) + ( ˜

Q, J∗

Q) + (˜

β, R∗

βΨ) + ( ˜

Q, R∗

QΨ)

= (˜ β, J∗

β + R∗ βΨ) + ( ˜

Q, J∗

Q + R∗ QΨ)

Since Ψ was arbitrary, and we want to eliminate ˜ Q, we can set J∗

Q + R∗ QΨ = 0

Then the derivative is given by J′(β)˜ β = (˜ β, J∗

β + R∗ βΨ)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 14 / 114

SLIDE 15

Adjoint approach

The adjoint equation R∗

QΨ + J∗ Q = 0

does not depend on the design variables β. Irrespective of the number

f design variables, we have to solve just one adjoint equation.

The gradient can be computed cheaply from J′(β)˜ β = (J∗

β + R∗ βΨ, ˜

β) which only requires scalar products. Thus the cost of adjoint approach to compute gradient is O(1), i.e., it is independent of the number of design variables.

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 15 / 114

SLIDE 16

One-shot versus iterative methods

One-shot method: solve for (Q, Ψ, β) simultaneously1

R(β, Q) = R∗

QΨ + J∗ Q

= J∗

β + R∗ βΨ

=

Iterative method: steepest descent

1 Initial guess for β 2 Solve for Q from R(β, Q) = 0 3 Solve for Ψ from R∗

QΨ + J∗ Q = 0

4 Compute gradient J′(β) = J∗

β + R∗ βΨ

5 Update β ←

− β − εJ′(β)

6 If not converged, goto step 2

1Hazra and Schulz, Simultaneous pseudo-timestepping for aerodynamic shape

ptimization
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 16 / 114

SLIDE 17

Adjoint: Two approaches

Continuous or Discrete or differentiate-discretize discretize-differentiate

PDE Adjoint PDE Discrete adjoint PDE Discrete PDE Discrete adjoint

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 17 / 114

SLIDE 18

Shape parameterization for gradient computation

Mapping2 a reference domain ˆ

Ω ⊂ Rd T : Rd → Rd, T ∈ Hm(Rd) Admissible shapes = {Ω : Ω = T(ˆ Ω)}

Normal displacements

◮ reference domain ˆ

Ω

◮ Γ ⊂ ∂ ˆ

Ω to be optimized

◮ perturb Γ by normal displacements

Γǫ = {x + ǫα(x)nΓ(x) : x ∈ Γ}, α ∈ Hm(Γ)

2Delfour & Zolesio, Shapes and geometries: Analysis, Differential Calculus and

Optimization, SIAM

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 18 / 114

SLIDE 19

Derivative of volume functional

J(Ω, u) =

Ω

F(u)dx, u = u(Ω) If A, B are two sets then A = (A − B) ∪ (A ∩ B) J(Ωǫ, uǫ)−J(Ω, u) =

Ωǫ−Ω

F(uǫ)dx−

Ω−Ωǫ

F(u)dx+

Ω∩Ωǫ

[F(uǫ)−F(u)]dx

n Γ+ Γ− Γ+ Γ−

Define the following boundary parts Γ+ = {s ∈ Γ : α(s) > 0} Γ− = {s ∈ Γ : α(s) ≤ 0}

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 19 / 114

SLIDE 20

Then we approximate the integrals as follows:

Ωǫ−Ω

F(uǫ)dx = +ǫ

Γ+ α(s)F(u)ds + O(ǫ2)
Ω−Ωǫ

F(u)dx = −ǫ

Γ− α(s)F(u)ds + O(ǫ2)

Define the sensitivity of u ˜ u = lim

ǫ→0

1 ǫ (uǫ − u)

r

uǫ − u = ǫ˜ u + O(ǫ2) Then

Ω∩Ωǫ

[F(uǫ) − F(u)]dx = ǫ

Ω

F ′(u)˜ udx + O(ǫ2) Shape derivative is 1 ǫ [J(Ωǫ) − J(Ω)] =

Γ

α(s)F(u)ds +

Ω

F ′(u)˜ udx + O(ǫ)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 20 / 114

SLIDE 21

Derivative of boundary functional

J(Γ, u) =

Γ

f(u)ds We perturb Γ by normal displacements Γǫ = {x + ǫα(s(x))n : x ∈ Γ} We must consider change in area element ds due to this perturbation. First look at a line element ds with radius of curvature R, ds = Rdθ Due to boundary perturbation, line element changes to dsǫ =

1 − ǫ α

R

ds

Same formula holds in 3-D but R is mean curvature defined using two

rthogonal tangential coordinates

1 R = 1 R1 + 1 R2

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 21 / 114

SLIDE 22

First express

Γǫ f(uǫ)ds in terms of integral and quantities on Γ

Consider x ∈ Γ and corresponding xǫ ∈ Γǫ given by xǫ = x + ǫαn u(xǫ) = u(x) + ǫα∂u ∂n(x) + O(ǫ2) uǫ(xǫ) = u(xǫ) + ǫ˜ u(xǫ) + O(ǫ2) = u(x) + ǫα∂u ∂n(x) + ǫ˜ u(x) + O(ǫ2) f(uǫ(xǫ)) = f

u(x) + ǫα∂u

∂n(x) + ǫ˜ u(x)

+ O(ǫ2)

= f(u) + ǫ

α∂u

∂n + ˜ u

f′(u) + O(ǫ2)

f(uǫ)dsǫ = f(u)ds + ǫ

αf′(u)∂u

∂n − α Rf(u)

ds + ǫf′(u)˜

u + O(ǫ2)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 22 / 114

SLIDE 23

Finally the sensitivity is given by 1 ǫ [J(Γǫ) − J(Γ)] =

Γ

α

f′(u)∂u

∂n − f(u) R

ds +
Γ

f′(u)˜ uds + O(ǫ)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 23 / 114

SLIDE 24

Euler equations: compressible, inviscid flow

ρ = density u = (u1, u2, u3) = velocity p = pressure Total energy per unit volume E = ρe + 1 2ρ|u|2 e = internal energy per unit mass For an ideal gas, e = p (γ − 1)ρ = ⇒ E = p γ − 1 + 1 2ρ|u|2 γ = Cp Cv For air γ = 1.4

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 24 / 114

SLIDE 25

Conservation of mass

Rate of change of mass inside volume Ω = – flux of mass going out through ∂Ω d dt

Ω

ρdx = −

∂Ω

ρ(u · n)ds Using divergence theorem d dt

Ω

ρdx = −

Ω

div(ρu)dx Since this should hold for any control volume Ω ∂ρ ∂t + div(ρu) = 0

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 25 / 114

SLIDE 26

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 dt

Ω(t)

ρuidx =

∂Ω(t)

Rids +

Ω(t)

ρfidx For an inviscid fluid, R = −pn. The left hand side is d dt

Ω(t)

ρuidx =

Ω(t)

∂ ∂t(ρui)dx +

∂Ω(t)

(ρui)(u · n)ds =

Ω(t)

∂ ∂t(ρui)dx +

Ω(t)

div(ρuiu)dx

Ω(t)

∂ ∂t(ρui)dx +

Ω(t)

div(ρuiu)dx = −

Ω

∂p ∂xi dx +

Ω

ρfidx ∂ ∂t(ρui) + div(ρuiu) + ∂p ∂xi = ρfi

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 26 / 114

SLIDE 27

Conservation of energy

Rate of change of energy = rate of work done by all forces d dt

Ω(t)

Edx =

∂Ω(t)

(R · u)ds +

Ω(t)

ρ(f · u)dx Using divergence theorem

Ω

∂E ∂t dx +

Ω

div(Eu)dx = −

Ω

div(pu)dx +

Ω(t)

ρ(f · u)dx ∂E ∂t + div[(E + p)u] = ρf · u

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 27 / 114

SLIDE 28

Euler equations

Mass conservation

∂ρ ∂t + ∂ ∂x1 (ρu1) + ∂ ∂x2 (ρu2) + ∂ ∂x3 (ρu3) = 0

Momentum conservation

∂ ∂t(ρu1) + ∂ ∂x1 (p + ρu2

1) +

∂ ∂x2 (ρu1u2) + ∂ ∂x3 (ρu1u3) = ∂ ∂t(ρu2) + ∂ ∂x1 (ρu1u2) + ∂ ∂x2 (p + ρu2

2) +

∂ ∂x3 (ρu2u3) = ∂ ∂t(ρu3) + ∂ ∂x1 (ρu1u3) + ∂ ∂x2 (ρu2u3) + ∂ ∂x3 (p + ρu2

3)

=

Energy conservation

∂E ∂t + ∂ ∂x1 [(E + p)u1] + ∂ ∂x2 [(E + p)u2] + ∂ ∂x3 [(E + p)u3] = 0

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 28 / 114

SLIDE 29

Euler equations: vector conservation form

∂U ∂t + ∂F1 ∂x1 + ∂F2 ∂x2 + ∂F3 ∂x3 = 0 where U is the conserved variables U =       ρ ρu1 ρu2 ρu3 E       and (F1, F2, F3) are the flux vectors F1 =       ρu1 p + ρu2

1

ρu1u2 ρu1u3 (E + p)u1       , F2 =       ρu2 ρu1u2 p + ρu2

2

ρu2u3 (E + p)u2       , F3 =       ρu3 ρu1u3 ρu2u3 p + ρu2

3

(E + p)u3      

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 29 / 114

SLIDE 30

Euler equations: Hyperbolic property

Flux jacobians

A = ∂F1 ∂U , B = ∂F2 ∂U , C = ∂F3 ∂U

Quasi-linear form

∂U ∂t + A ∂U ∂x1 + B ∂U ∂x2 + C ∂U ∂x3 = 0

Hyperbolic property: A(U, ω) = ω1A + ω2B + ω3C,

ω ∈ R3

◮ All eigenvalues are real ◮ Eigenvectors are linearly independent

Eigenvalues of A(U, ω) with |ω| = 1 are

u · ω − a, u · ω, u · ω, u · ω, u · ω + a a = speed of sound = γp ρ

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 30 / 114

SLIDE 31

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., min

Γ

1 2

Γ

(p − p∗)2ds The pressure p = p(U) is obtained from the solution of the steady Euler equations Fx + Gy + Hz = 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

SLIDE 32

Due to change in shape Γ, the solution changes from U to U + ǫ ˜ U F(U + ǫ ˜ U) = F(U) + ǫ∂F ∂U (U) ˜ U + O(ǫ2) = 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 ˜ U)x = −

Ω

Vx · A ˜ U +

∂Ω

V · A ˜ Unx = −

Ω

AtVx · ˜ U +

∂Ω

V · A ˜ Unx Then

Ω

V · [(A ˜ U)x + (B ˜ U)y + (C ˜ U)z] = 0

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 32 / 114

SLIDE 33

becomes −

Ω

(AtVx + BtVy + CtVz) · ˜ U +

∂Ω

V · (Anx + Bny + Cnz) ˜ U = 0 Define D = Anx + Bny + Cnz we note that DU =   ρ(u · n) pn + ρu(u · n) (E + p)(u · n)   = normal flux and on solid wall Γ, u · n = 0 DU|Γ =   pn  

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 33 / 114

SLIDE 34

Since DU = Fnx + Gny + Hnz

DU

= ˜ Fnx + ˜ Gny + ˜ Hnz = A ˜ Unx + B ˜ Uny + C ˜ Unz = D ˜ U Also

pn = ˜

pn + p˜ n, ˜ n = −

d−1

j=1

∂α ∂tj tj Denoting V = [v1, v, v5] with v = (v2, v3, v4)

Γ

V · D ˜ U =

Γ

V · DU =

Γ

V · [0

pn

0] =

Γ

v · pn =

Γ

˜ p(v · n) −

Γ

p

d−1

j=1

∂α ∂tj (v · tj)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 34 / 114

SLIDE 35

First variation in objective function 1 ǫ ˜ J =

Γ
(p − p∗)˜

p + α(p − p∗) ∂p ∂n − α(p − p∗)2 2R

To obtain gradient of J wrt α we must eliminate ˜
p. We add the first
rder perturbation equation to the above equation

1 ǫ ˜ J = RHS −

Ω

(AtVx + BtVy + CtVz) · ˜ U +

∂Ω

V · D ˜ U To eliminate ˜ U from the volume integral, we choose V to satisfy AtVx + BtVy + CtVz = 0, in Ω

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 35 / 114

SLIDE 36

1 ǫ ˜ J = RHS +

Γ

V · D ˜ U +

Γo

V · D ˜ U =

Γ

(p − p∗ + v · n)˜ p +

Γ

α

(p − p∗) ∂p

∂n − (p − p∗)2 2R

−
Γ

p

d−1

j=1

∂α ∂tj (v · tj) +

Γo

DtV · ˜ U To eliminate ˜ p we choose the adjoint boundary conditions p − p∗ + v · n = 0,

n

Γ

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 36 / 114

SLIDE 37

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:
Γo

DtV · ˜ U =

Γo

DtV · T −1T ˜ U =

Γo

T −tDtV · T ˜ U We have (T ˜ U)1,2,3,4 = 0 while (T ˜ U)5 is arbitrary, so we choose (T −tDtV )5 = 0

Subsonic outflow: (T ˜

U)1 is specified while (T ˜ U)2,3,4,5 is arbitrary. Hence we choose (T −tDtV )2,3,4,5 = 0

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 37 / 114

SLIDE 38

Adjoint Euler equations

AtVx + BtVy + CtVz = 0, in Ω with boundary conditions

supersonic inflow

none

supersonic outflow

V = 0

subsonic inflow

(T −tDtV )5 = 0

subsonic outflow

(T −tDtV )2,3,4,5 = 0

solid wall

p − p∗ + v · n = 0

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 38 / 114

SLIDE 39

Shape derivative

1 ǫ ˜ J =

Γ

α

(p − p∗) ∂p

∂n − (p − p∗)2 2R

−
Γ

p

d−1

j=1

∂α ∂tj (v · tj) =

Γ

α

(p − p∗) ∂p

∂n − (p − p∗)2 2R − divΓ(pv)

Hence

∇αJ = (p − p∗) ∂p ∂n − (p − p∗)2 2R − divΓ(pv) Steepest descent update αn+1 = αn − sn∇αJ(αn)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 39 / 114

SLIDE 40

Discrete adjoint approach

Approximate PDE using finite volume method

R(X, Q) = 0 where X = grid coordinates, X ∈ Rm Q = solution vector, Q ∈ Rn R : Rm × Rn → Rn

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

SLIDE 41

Elements of practical shape optimization

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

SLIDE 42

Adjoint approach:

dI dX

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

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

Differentiate state equation R(X, Q) = 0

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

Flow sensitivity ∂Q

∂X ; costly to evaluate

Introducing an adjoint variable Ψ ∈ Rn, we can write

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

+ Ψt

∂R ∂X + ∂R ∂Q ∂Q ∂X

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 42 / 114

SLIDE 43

Adjoint approach

Collect terms involving the flow sensitivity

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

+

∂I ∂Q + Ψt ∂R ∂Q ∂Q ∂X

Choose Ψ so that flow sensitivity vanishes

∂I ∂Q + Ψt ∂R ∂Q = 0 = ⇒ ∂R ∂Q t Ψ + ∂I ∂Q t = 0

Gradient

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

Cost of Adjoint-based steepest-descent

Cost = O(1)Niter = O(N)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 43 / 114

SLIDE 44

Optimization steps

β =

⇒ Xs = ⇒ X

Solve the flow (primal) equations to steady-state

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

Solve adjoint equations to steady-state

X, Q = ⇒ dΨ dt + ∂R ∂Q t Ψ + ∂I ∂Q t = 0 = ⇒ Ψ

Compute gradient wrt grid X

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

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 44 / 114

SLIDE 45

Automatic Differentiation

Discrete adjoint method requires

∂J ∂X , ∂J ∂Q, ∂R ∂X t Ψ, ∂R ∂Q t Ψ

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

Code for J(X, Q) AD tool Code for

∂J ∂X

X

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 45 / 114

SLIDE 46

Automatic Differentiation: forward and backward mode

X ∈ Rm, Y ∈ Rn, Y = F(X)

Forward mode

X ∈ Rm, ˙ X ∈ Rm → ∂F ∂X

˙

X ∈ Rn

Backward mode

X ∈ Rm, ¯ Y ∈ Rn → ∂F ∂X t ¯ Y ∈ Rm

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 46 / 114

SLIDE 47

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 Generates new code residue b.f : X, Q, Ψ →

∂R

∂Q

t Ψ 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)

Ψ

∂R ∂Q ⎛ ⎝ ⎜ ⎞ ⎠ ⎟

T

Ψ

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 47 / 114

SLIDE 48

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. t1 = L1(T0) t2 = L2(T1) . . . = . . . tp = Lp(Tp−1) f = Lp+1(Tp) x ∈ Rn, T0 = x tr ∈ R Tr = Tr−1 tr

The function f as coded, in general, depends on all of these variables

x, t1, t2, . . . , tp

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 48 / 114

SLIDE 49

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

SLIDE 50

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

SLIDE 51

Examples of shape optimization

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 51 / 114

SLIDE 52

Quasi 1-D flow

Inflow Outflow h(x) x

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 52 / 114

SLIDE 53

Quasi 1-D flow

Quasi 1-D flow in a duct

∂ ∂t(hU) + ∂ ∂x(hf) = dh dxP, x ∈ (a, b) t > 0 U =   ρ ρu E   , f =   ρu p + ρu2 (E + p)u   , P =   p   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

J = b

a

(p − p∗)2dx

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 53 / 114

SLIDE 54

Quasi 1-D flow

1 i i − 1 / 2 i + 1 / 2 N Inflow face Outflow face

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 54 / 114

SLIDE 55

Quasi 1-D flow

Finite volume scheme

hi dUi dt + hi+1/2Fi+1/2 − hi−1/2Fi−1/2 ∆x = (hi+1/2 − hi−1/2) ∆x Pi

Discrete cost function

J =

N

i=1

(pi − p∗

i )2

Control variables

h1/2, h1+1/2, . . . , hi+1/2, . . . , hN+1/2

N = 100
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 55 / 114

SLIDE 56

Duct shape

2 4 6 8 10 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 x h Duct shape Target shape Current shape

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 56 / 114

SLIDE 57

Target pressure distribution p∗

2 4 6 8 10 0.5 1 1.5 2 2.5 x Pressure Target pressure AUSMDV KFVS LF

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 57 / 114

SLIDE 58

Current pressure distribution

2 4 6 8 10 0.5 1 1.5 2 2.5 x Pressure Starting pressure AUSMDV KFVS LF

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 58 / 114

SLIDE 59

Adjoint density

2 4 6 8 10 −35 −30 −25 −20 −15 −10 −5 5 x Adjoint density AUSMDV KFVS LF

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 59 / 114

SLIDE 60

Convergence history: Explicit Euler

1000 2000 3000 4000 5000 6000 7000 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 No of iterations Residue Convergence history with AUSMDV flux Flow Adjoint

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 60 / 114

SLIDE 61

Shape gradient

2 4 6 8 10 5 10 15 20 25 30 35 x Gradient Gradient AUSMDV KFVS LF

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 61 / 114

SLIDE 62

Shape gradient

2 2.5 3 3.5 4 4.5 5 5.5 6 −5 5 10 15 20 25 30 35 40 x Gradient Gradient AUSMDV KFVS LF

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 62 / 114

SLIDE 63

Validation of Shape gradient

2 4 6 8 10 5 10 15 20 25 30 35 x Gradient Gradient with AUSMDV flux A B C

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 63 / 114

SLIDE 64

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

SLIDE 65

Gradient smoothing

Non-smooth gradients G especially in the presence of shocks
Smooth using an elliptic equation
1 − ǫ(x) d2

dx2

¯

G(x) = G(x) ǫi = {|Gi+1 − Gi| + |Gi − Gi−1|}Li Li = |Gi+1 − 2Gi + Gi−1| max(|Gi+1 − Gi| + |Gi − Gi−1|, TOL)

Finite difference with Jacobi iterations
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 65 / 114

SLIDE 66

Gradient smoothing

2 4 6 8 10 5 10 15 20 25 30 35 x Gradient Gradient using AUSMDV flux Original gradient Smoothed gradient

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 66 / 114

SLIDE 67

Quasi 1-D optimization: Shape

1 2 3 4 5 6 7 8 9 10 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 x h Target Initial

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 67 / 114

SLIDE 68

Quasi 1-D optimization: Final shape

1 2 3 4 5 6 7 8 9 10 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 x h Target Initial Final

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 68 / 114

SLIDE 69

Quasi 1-D optimization: Pressure

1 2 3 4 5 6 7 8 9 10 0.5 1 1.5 2 2.5 x Pressure Target Initial Final

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 69 / 114

SLIDE 70

Quasi 1-D optimization: Convergence

5 10 15 20 25 30 35 40 45 50 10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10

No. of iterations

Cost/Gradient Cost Gradient

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 70 / 114

SLIDE 71

Quasi 1-D optimization: Adjoint density

2 4 6 8 10 −25 −20 −15 −10 −5 x Adjoint density Initial Final

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 71 / 114

SLIDE 72

Airfoil shape optimization

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 72 / 114

SLIDE 73

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

SLIDE 74

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

SLIDE 75

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

SLIDE 76

Validation of adjoint gradients

Dot-product test

˙ R := R′(Q) ˙ Q ≈ R(Q + ǫ ˙ Q) − R(Q) ǫ ¯ Q := [R′(Q)]⊤ ¯ R ˙ R⊤ ¯ R

?

= ˙ Q⊤ ¯ Q

Upwind schemes and limiters

can cause problems due to non-differentiability

Check adjoint derivatives

against finite difference

◮ NACA0012: Cd/Cl ◮ RAE2822: Cd

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

SLIDE 77

Convergence and limiter: RAE2822 airfoil

L(a, b) = a(b2 + 2ǫ) + b(2a2 + ǫ) 2a2 − ab + 2b2 + 3ǫ , 0 < ǫ ≪ 1 ǫ = 10−8 ǫ = 10−4

1e-05 0.0001 0.001 0.01 0.1 1 200 400 600 800 1000 1200 1400 1600 Residue Number of iterations Density x-momen y-momen z-momen Energy 1e-10 1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 200 400 600 800 1000 1200 1400 1600 Residue Number of iterations Density x-momen y-momen z-momen Energy

Adjoint iterations blow-up No blow-up

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 77 / 114

SLIDE 78

Iterative convergence tests

Primal residual Adjoint residual Rn = ||R(Qn)|| Rn = ||[R′(Q∞)]⊤Ψn + I′(Q∞)||

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

SLIDE 79

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

SLIDE 80

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

SLIDE 81

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

SLIDE 82

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

SLIDE 83

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

SLIDE 84

RANS computation

0.2 0.4 0.6 0.8 1 x/c

1
0.5

0.5 1 1.5

Cp

RAE2822 Optimized M=0.729, Re=6.5 million, Cl=0.88, k-w turbulence model

Shape Inviscid RANS AOA Cl 100Cd AOA Cl 100Cd 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

SLIDE 85

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

SLIDE 86

Adjoints with shocks: Scalar conservation law

Consider the functional J = +1

−1

G(u(x, T))dx where u is the solution of the conservation law ut + f(u)x = 0, x ∈ (0, 1), t ∈ (0, T) u(x, 0) = u0(x) To compute the derivative of J wrt u0 we make a small perturbation u0 − → u0 + ˜

u0. Then solution changes to u + ˜

u and ˜ J =

g(u(x, T))˜

udx, g(u) = dG du The equation for ˜ u is obtained by linearizing the conservation law (˜ u)t + (f′(u)˜ u)x = ˜ u(x, 0) = ˜ u0

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 86 / 114

SLIDE 87

To write ˜ J in terms of ˜ u0, we introduce the adjoint variable v(x, t) ˜ J =

g(u(x, T))˜

u(x, T)dx − T

v
(˜

u)t + (f′(u)˜ u)x

dxdt

Integrating by parts, we get ˜ J =

v(x, 0)˜

u0(x)dx +

(g − v)˜

u(x, T)dx − T vt + f′(u)vx

˜

u(x, t)dxdt To eliminate ˜ u(x, t) and ˜ u(x, T) from above equation, choose vt + f′(u)vx = v(x, T) = g(u(x, T)) The derivative is given by ˜ J =

v(x, 0)˜

u0(x)dx

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 87 / 114

SLIDE 88

Note that the adjoint equation must be solved backward in time, starting from t = T to t = 0.

Example: Inviscid Burger’s equation

ut + (u2/2)x = 0, x ∈ R, t ∈ (0, T) u(x, 0) =

ul

x < 0 ur x > 0 Exact solution: u(x, t) =

ul

x < St ur x > St S = 1 2(ul + ur) Consider the objective functional J = 1 2

u2(x, T)dx
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 88 / 114

SLIDE 89

Then the corresponding adjoint equation is vt + uvx = v(x, T) = u(x, T) =

ul

x < ST ur x > ST

x t t = T

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

SLIDE 90

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 = xs, 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 = dxs

dt is the shock speed.

The linearization of the conservation law must also include the linearization of the jump conditions. xs → xs + ˜ xs, dxs dt → dxs dt + d˜ xs dt d˜ xs dt [u] + dxs dt [˜ u] + dxs dt ˜ xs ∂u ∂x

=

df du ˜ u

+ ˜

xs df du ∂u ∂x

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 90 / 114

SLIDE 91

J = xs(T)

−1

G(u(x, T))dx + +1

xs(T)

G(u(x, T))dx Perturbation in J is ˜ J =

xs(T)

−1

g(u(x, T))˜ u(x, T)dx+

+1

xs(T)

g(u(x, T))˜ u(x, T)dx−[G]xs(T)˜ xs(T) −

T

xs(t)
−1

v

˜

ut + (f′(u)˜ u)x

dxdt −

T

+1
xs(t)

v

˜

ut + (f′(u)˜ u)x

dxdt

−

T

vs(t)

d˜ xs dt [u] + dxs dt [˜ u] + dxs dt ˜ xs ∂u ∂x

−

df du ˜ u

− ˜

xs df du ∂u ∂x

dt
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 91 / 114

SLIDE 92

After some work we get ˜ J =

+1

−1

v(x, 0)˜ u0(x)dx provided the adjoint variables v(x, t) and vs(t) satisfy the following equations vt + f′(u)vx = v(x, T) = g(u(x, T)) vs(t) = v(x+

s (t), t) = v(x− s (t), t)

d dtvs(t) = vs(T) = [G] [u] Adjoint variable requires an interior boundary condition along the shock (xs(t), t). This uniquely determines the adjoint solution.

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 92 / 114

SLIDE 93

Steepest descent: limitations, problems

Multiple optima are common
Convergence towards local optimum:

dependance on starting guess So

Shape gradient: difficult to compute, or impossible
Noisy cost function J
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 93 / 114

SLIDE 94

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 = {xn

1, xn 2, . . . , xn M} ⊂ Rd

Solutions evolve according to some rules

P n+1 = E(P n)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 94 / 114

SLIDE 95

Particle swarm optimization

Kennedy and Eberhart (1995)
Modeled on behaviour of animal swarms: ants, bees, birds
Swarm intelligence

Optimization problem

min

x∈D J(x),

D ⊂ Rd

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 95 / 114

SLIDE 96

Particle swarm optimization

Kennedy and Eberhart (1995)
Modeled on behaviour of animal swarms: ants, bees, birds
Swarm intelligence

Optimization problem

min

x∈D J(x),

D ⊂ Rd

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 95 / 114

SLIDE 97

Particle swarm optimization

Particles distributed in design space xi ∈ D, i = 1, ..., Np

X2 X1

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 96 / 114

SLIDE 98

Particle swarm optimization

Each particle has a velocity vi ∈ Rd, i = 1, ..., Np

X2 X1

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 97 / 114

SLIDE 99

Particle swarm optimization

Particles have memory (t = iteration number)

Local memory : pt

i = argmin 0≤s≤t

J(xs

i)

Global memory : pt = argmin

i

J(pt

i)

Velocity update

vt+1

i

= ωvt

i + c1rt 1(pt i − xt i)

Local

+ c2rt

2(pt − xt i)

Global
Position update

xt+1

i

= xt

i + vt+1 i

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 98 / 114

SLIDE 100

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

SLIDE 101

PSO: Parallelizability

xn

1

xn

2

. . . xn

Np

↓ ↓ ↓ J(xn

1)

J(xn

2)

. . . J(xn

Np)

↓ ↓ ↓ vn+1

1

vn+1

2

. . . vn+1

Np

↓ ↓ ↓ xn+1

1

xn+1

2

. . . xn+1

Np

Parallel evaluation of cost functions using MPI: Message Passing Interface

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 100 / 114

SLIDE 102

Free Form Deformation

Originated in computer graphics field
Embed the object inside a box and deform the box
Independent of the representation of the object
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 101 / 114

SLIDE 103

Free Form Deformation: Example

(Duvigneau, OPALE, INRIA)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 102 / 114

SLIDE 104

Free Form Deformation

Xo(P) = coordinate of point P wrt reference shape
Movement of point P under the deformation

X(P) = Xo(P) +

ni

i=0

nj

j=0

nk

k=0

Bni

i (ξP )Bnj j (ηP )Bnk k (ζP )Yijk

Bernstein polynomials

Bn

p (t) = Cn p tp(1 − t)n−p

Design variables

{Yijk}, 0 ≤ i ≤ ni, 0 ≤ j ≤ nj, 0 ≤ k ≤ nk

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 103 / 114

SLIDE 105

Test case: Wing shape optimization3

Minimize drag under lift

constraint min Cd Cdo s.t. Cl Clo ≥ 0.999

FFD parameterization, n = 20

design variables

Particle swarm optimization:

120 particles M∞ = 0.83, α = 2o (Piaggio Aero. Ind.) Grid: 31124 nodes

Cost function

J = Cd Cdo + 104 max

0, 0.999 − Cl

Clo

3Joint work with R. Duvigneau, OPALE, INRIA
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 104 / 114

SLIDE 106

Test case: Wing shape optimization3

Minimize drag under lift

constraint min Cd Cdo s.t. Cl Clo ≥ 0.999

FFD parameterization, n = 20

design variables

Particle swarm optimization:

120 particles M∞ = 0.83, α = 2o (Piaggio Aero. Ind.) Grid: 31124 nodes

Cost function

J = Cd Cdo + 104 max

0, 0.999 − Cl

Clo

3Joint work with R. Duvigneau, OPALE, INRIA
Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 104 / 114

SLIDE 107

Wing optimization

Initial shape

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 105 / 114

SLIDE 108

Wing optimization

Optimized shape

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 106 / 114

SLIDE 109

PSO computational cost

Slow convergence: a few hundred iterations
Require large swarm size: O(100)
CFD is expensive: few minutes to hours
Example: Wing optimization (coarse CFD grid)

(10 minutes/CFD) (120 CFD/iteration) (200 iterations) = 4000 hours

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 107 / 114

SLIDE 110

Robust design

Variability of the fitness J(x, A) due to uncertain parameters A (mach number, angle of attack, shape, etc.)

Example : effect of Mach number fluctuations

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 108 / 114

SLIDE 111

Optimization problem

Optimization

min

x∈Rd J (x, ao)

C(x, ao) ≤ 0

Robust optimization

min

x∈Rd

µJ(x) =

Ω(A)

J (x, a) ρA(a) da σ2

J(x) =

Ω(A)

[J (x, a) − µJ]2 ρA(a) da Prob[C(x, A) ≤ 0] ≥ p

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 109 / 114

SLIDE 112

Optimization problem

Optimization

min

x∈Rd J (x, ao)

C(x, ao) ≤ 0

Robust optimization

min

x∈Rd

µJ(x) =

Ω(A)

J (x, a) ρA(a) da σ2

J(x) =

Ω(A)

[J (x, a) − µJ]2 ρA(a) da Prob[C(x, A) ≤ 0] ≥ p

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 109 / 114

SLIDE 113

Monte-Carlo estimation using meta-models

For each design x: compute

J(x, a1), J(x, a2), . . . , J(x, aP ) ⇓ meta-model ˜ Jx(a)

Monte-Carlo estimation of the statistics: Ns ≫ 1

MJ(x) = 1 Ns

Ns

i=1

˜ Jx(ai) S2

J(x) =

1 Ns − 1

Ns

i=1

[ ˜ Jx(ai) − MJ]2

Robust optimization problem: 0 ≤ ω ≤ 1

min

x∈Rd

ωMJ + (1 − ω)SJ + constraint penalty

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 110 / 114

SLIDE 114

Monte-Carlo estimation using meta-models

For each design x: compute

J(x, a1), J(x, a2), . . . , J(x, aP ) ⇓ meta-model ˜ Jx(a)

Monte-Carlo estimation of the statistics: Ns ≫ 1

MJ(x) = 1 Ns

Ns

i=1

˜ Jx(ai) S2

J(x) =

1 Ns − 1

Ns

i=1

[ ˜ Jx(ai) − MJ]2

Robust optimization problem: 0 ≤ ω ≤ 1

min

x∈Rd

ωMJ + (1 − ω)SJ + constraint penalty

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 110 / 114

SLIDE 115

Test case: Wing shape optimization

Minimize the drag mean and

the drag variance

Probabilistic lift constraint

p = 0.95

FFD parameterization, n = 32

design variables

Particle swarm optimization:

32 particles

Uncertain Mach number

M∞ = 0.83, α = 2o (Piaggio Aero. Ind.) Grid: 31124 nodes Number of CFD computations/iteration = 32 × 4 = 128

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 111 / 114

SLIDE 116

PDF of Mach number

Mean = 0.83, std. dev. = 0.0166

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 112 / 114

SLIDE 117

Wing optimization: result

Drag vs Mach number PDF of drag

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 113 / 114

SLIDE 118

Wing optimization: flow solutions

Optimal design for M∞ = 0.83 M∞ = 0.81 M∞ = 0.83 M∞ = 0.85 Robust design, weights (0.5, 0.5)

Praveen. C

(TIFR-CAM) Shape Optimization IISc, Dec 2009 114 / 114