[PPT] - Adjoint-Based Optimization of Time-Dependent Fluid-Structure Systems PowerPoint Presentation

SLIDE 1

Adjoint-Based Optimization of Time-Dependent Fluid-Structure Systems using a High-Order Discontinuous Galerkin Discretization

Matthew J. Zahr† and Per-Olof Persson

ECCOMAS High-Order Nonlinear Numerical Methods (HONOM) University of Stuttgart, Stuttgart, Germany March 27-31, 2017

† Luis W. Alvarez Postdoctoral Fellow

Department of Mathematics Lawrence Berkeley National Laboratory University of California, Berkeley 1 / 37

SLIDE 2

PDE optimization is ubiquitous in science and engineering

Design: Find system that optimizes performance metric, satisfies constraints Aerodynamic shape design of automobile Optimal flapping motion of micro aerial vehicle

2 / 37

SLIDE 3

PDE optimization is ubiquitous in science and engineering

Control: Drive system to a desired state Boundary flow control Metamaterial cloaking – electromagnetic invisibility

3 / 37

SLIDE 4

PDE optimization is ubiquitous in science and engineering

Inverse problems: Infer the problem setup given solution observations Left: Material inversion – find inclusions from acoustic, structural measurements Right: Source inversion – find source of airborne contaminant from downstream measurements Full waveform inversion – estimate subsurface of Earth’s crust from acoustic measurements

4 / 37

SLIDE 5

Time-dependent PDE-constrained optimization

Introduction of fully discrete adjoint method

emanating from high-order discretization of governing equations

Time-periodicity constraints
Extension to high-order partitioned solver for

fluid-structure interaction

Applications: flapping flight, energy harvesting

Micro Aerial Vehicle Vertical Windmill Volkswagen Passat LES Flow past Airfoil

5 / 37

SLIDE 6

Unsteady PDE-constrained optimization formulation

Goal: Find the solution of the unsteady PDE-constrained optimization problem minimize

U, µ

J (U, µ) subject to C(U, µ) ≤ 0 ∂U ∂t + ∇ · F (U, ∇U) = 0 in v(µ, t) where

U(x, t)

PDE solution

µ

design/control parameters

J (U, µ) =

Tf

T0

Γ

j(U, µ, t) dS dt

bjective function
C(U, µ) =

Tf

T0

Γ

c(U, µ, t) dS dt constraints

6 / 37

SLIDE 7

Nested approach to PDE-constrained optimization

PDE optimization requires repeated queries to primal and dual PDE

Primal PDE Dual PDE Optimizer 7 / 37

SLIDE 8

Nested approach to PDE-constrained optimization

PDE optimization requires repeated queries to primal and dual PDE

Primal PDE Dual PDE Optimizer µ 7 / 37

SLIDE 9

Nested approach to PDE-constrained optimization

PDE optimization requires repeated queries to primal and dual PDE

Primal PDE Dual PDE Optimizer J (U, µ) 7 / 37

SLIDE 10

Nested approach to PDE-constrained optimization

PDE optimization requires repeated queries to primal and dual PDE

Primal PDE Dual PDE Optimizer J (U, µ) µ U 7 / 37

SLIDE 11

Nested approach to PDE-constrained optimization

PDE optimization requires repeated queries to primal and dual PDE

Primal PDE Dual PDE Optimizer J (U, µ)

dJ dµ (U, µ)

7 / 37

SLIDE 12

High-order discretization of PDE-constrained optimization

Continuous PDE-constrained optimization problem

minimize

U, µ

J (U, µ) subject to C(U, µ) ≤ 0 ∂U ∂t + ∇ · F (U, ∇U) = 0 in v(µ, t)

Fully discrete PDE-constrained optimization problem

minimize

u0, ..., uNt∈RNu, k1,1, ..., kNt,s∈RNu, µ∈Rnµ

J(u0, . . . , uNt, k1,1, . . . , kNt,s, µ) subject to C(u0, . . . , uNt, k1,1, . . . , kNt,s, µ) ≤ 0 u0 − g(µ) = 0 un − un−1 −

s

i=1

bikn,i = 0 Mkn,i − ∆tnr (un,i, µ, tn,i) = 0

8 / 37

SLIDE 13

Highlights of globally high-order discretization

Arbitrary Lagrangian-Eulerian formulation:

Map, G(·, µ, t), from physical v(µ, t) to reference V ∂U X ∂t

X

+ ∇X · F X(U X, ∇XU X) = 0

Space discretization: discontinuous Galerkin

M ∂u ∂t = r(u, µ, t)

Time discretization: diagonally implicit RK

un = un−1 +

s

i=1

bikn,i Mkn,i = ∆tnr (un,i, µ, tn,i)

Quantity of interest: solver-consistency

F(u0, . . . , uNt, k1,1, . . . , kNt,s)

X1 X2 NdA V x1 x2 nda v G, g, vX

Mapping-Based ALE DG Discretization c1 a11 c2 a21 a22 . . . . . . . . . ... cs as1 as2 · · · ass b1 b2 · · · bs Butcher Tableau for DIRK

9 / 37

SLIDE 14

Adjoint method to efficiently compute gradients of QoI

Consider the fully discrete output functional F(un, kn,i, µ)
Represents either the objective function or a constraint
The total derivative with respect to the parameters µ, required in the context
f gradient-based optimization, takes the form

dF dµ = ∂F ∂µ +

Nt

n=0

∂F ∂un ∂un ∂µ +

Nt

n=1

s

i=1

∂F ∂kn,i ∂kn,i ∂µ

The sensitivities, ∂un

∂µ and ∂kn,i ∂µ , are expensive to compute, requiring the solution of nµ linear evolution equations

Adjoint method: alternative method for computing dF

dµ that require one linear evolution equation for each quantity of interest, F

10 / 37

SLIDE 15

Adjoint equation derivation: outline

Define auxiliary PDE-constrained optimization problem

minimize

u0, ..., uNt∈RNu, k1,1, ..., kNt,s∈RNu

F(u0, . . . , uNt, k1,1, . . . , kNt,s, µ) subject to R0 = u0 − g(µ) = 0 Rn = un − un−1 −

s

i=1

bikn,i = 0 Rn,i = Mkn,i − ∆tnr (un,i, µ, tn,i) = 0

Define Lagrangian

L(un, kn,i, λn, κn,i) = F − λ0

T R0 − Nt

n=1

λn

T Rn − Nt

n=1

s

i=1

κn,i

T Rn,i

The solution of the optimization problem is given by the

Karush-Kuhn-Tucker (KKT) sytem ∂L ∂un = 0, ∂L ∂kn,i = 0, ∂L ∂λn = 0, ∂L ∂κn,i = 0

11 / 37

SLIDE 16

Dissection of fully discrete adjoint equations

Linear evolution equations solved backward in time
Primal state/stage, un,i required at each state/stage of dual problem
Heavily dependent on chosen ouput

λNt = ∂F ∂uNt

T

λn−1 = λn + ∂F ∂un−1

T

+

s

i=1

∆tn ∂r ∂u (un,i, µ, tn−1 + ci∆tn)T κn,i M T κn,i = ∂F ∂uNt

T

+ biλn +

s

j=i

aji∆tn ∂r ∂u (un,j, µ, tn−1 + cj∆tn)T κn,j

Gradient reconstruction via dual variables

dF dµ = ∂F ∂µ + λ0

T ∂g

∂µ(µ) +

Nt

n=1

∆tn

s

i=1

κn,i

T ∂r

∂µ(un,i, µ, tn,i) [Zahr and Persson, 2016]

12 / 37

SLIDE 17

Energetically optimal flapping under x-impulse constraint

minimize

µ

− 3T

2T

Γ

f · v dS dt subject to 3T

2T

Γ

f · e1 dS dt = q U(x, 0) = g(x) ∂U ∂t + ∇ · F (U, ∇U) = 0

Isentropic, compressible,

Navier-Stokes

Re = 1000, M = 0.2
y(t), θ(t), c(t) parametrized via

periodic cubic splines

Black-box optimizer: SNOPT

y(t) θ(t) l l/3 c(t)

Airfoil schematic, kinematic description

13 / 37

SLIDE 18

Optimal control - fixed shape

Fixed shape, optimal Rigid Body Motion (RBM), varied x-impulse Energy = 9.4096 x-impulse = -0.1766 Energy = 0.45695 x-impulse = 0.000 Energy = 4.9475 x-impulse = -2.500 Initial Guess Optimal RBM Jx = 0.0 Optimal RBM Jx = −2.5

14 / 37

SLIDE 19

Optimal control, time-morphed geometry

Optimal Rigid Body Motion (RBM) and Time-Morphed Geometry (TMG), varied x-impulse Energy = 9.4096 x-impulse = -0.1766 Energy = 0.45027 x-impulse = 0.000 Energy = 4.6182 x-impulse = -2.500 Initial Guess Optimal RBM/TMG Jx = 0.0 Optimal RBM/TMG Jx = −2.5

15 / 37

SLIDE 20

Optimal control, time-morphed geometry

Optimal Rigid Body Motion (RBM) and Time-Morphed Geometry (TMG), x-impulse = −2.5 Energy = 9.4096 x-impulse = -0.1766 Energy = 4.9476 x-impulse = -2.500 Energy = 4.6182 x-impulse = -2.500 Initial Guess Optimal RBM Jx = −2.5 Optimal RBM/TMG Jx = −2.5

16 / 37

SLIDE 21

Energetically optimal flapping in three-dimensions

Energy = 1.4459e-01 Thrust = -1.1192e-01 Energy = 3.1378e-01 Thrust = 0.0000e+00

17 / 37

SLIDE 22

Time-periodic solutions desired when optimizing cyclic motion

To properly optimize a cyclic, or periodic problem, need to simulate a

representative period

Necessary to avoid transients that will impact quantity of interest and may

cause simulation to crash

Task: Find initial condition, u0, such that flow is periodic, i.e. uNt = u0

18 / 37

SLIDE 23

Time-periodic solutions desired when optimizing cyclic motion

Vorticity around airfoil with flow initialized from steady-state (left) and time-periodic flow (right) 2 4 −60 −40 −20 time power 2 4 −4 −2 time power Time history of power on airfoil of flow initialized from steady-state ( ) and from a time-periodic solution ( )

19 / 37

SLIDE 24

Time-periodicity constraints in PDE-constrained optimization

Slight modification leads to fully discrete periodic PDE-constrained optimization minimize

u0, ..., uNt∈RNu, k1,1, ..., kNt,s∈RNu, µ∈Rnµ

J(u0, . . . , uNt, k1,1, . . . , kNt,s, µ) subject to C(u0, . . . , uNt, k1,1, . . . , kNt,s, µ) ≤ 0 u0 − uNt = 0 un − un−1 +

s

i=1

bikn,i = 0 Mkn,i − ∆tnr (un,i, µ, tn,i) = 0

20 / 37

SLIDE 25

Adjoint method for periodic PDE-constrained optimization

Following identical procedure as for non-periodic case, the adjoint equations

corresponding to the periodic conservation law are λNt = λ0 + ∂F ∂uNt

T

λn−1 = λn + ∂F ∂un−1

T

+

s

i=1

∆tn ∂r ∂u (un,i, µ, tn−1 + ci∆tn)T κn,i M T κn,i = ∂F ∂uNt

T

+ biλn +

s

j=i

aji∆tn ∂r ∂u (un,j, µ, tn−1 + cj∆tn)T κn,j

Dual problem is also periodic
Solve linear, periodic problem using Krylov shooting method

[Zahr et al., 2016]

21 / 37

SLIDE 26

Energetically optimal flapping: x-impulse, time-periodic

minimize

µ

− T

Γ

f · v dS dt subject to T

Γ

f · e1 dS dt = q U(x, 0) = U(x, T) ∂U ∂t + ∇ · F (U, ∇U) = 0

Isentropic, compressible,

Navier-Stokes

Re = 1000, M = 0.2
y(t), θ(t), c(t) parametrized via

periodic cubic splines

Black-box optimizer: SNOPT

y(t) θ(t) l l/3

Airfoil schematic, kinematic description

22 / 37

SLIDE 27

Solution of time-periodic, energetically optimal flapping

Energy = 9.4096 x-impulse = -0.1766 Energy = 0.45695 x-impulse = 0.000

23 / 37

SLIDE 28

Structure: semi-discretization, first-order form

M s ∂us

∂t = rs(us; t) = rss(us) + rsf · t

Semidiscretization (CG-FEM) of continuum (hyperelasticity)

∂p ∂t − ∇ · P (G) = b in Ω0 P (G) · N = t

n ΓN

x = xD

n ΓD
Force balance on rigid body

M ∂2q ∂t2 + C ∂q ∂t + Kq = t

c θ(µ, t) h(us) 24 / 37

SLIDE 29

Coupled fluid-structure formulation

Write discretized fluid and structure equations as ODEs

M f ˙ uf = rf(uf; x) M s ˙ us = rs(us; t) = rss(us) + rsf · t in the fluid uf and structure us variables

Apply couplings
Structure-to-fluid: deform fluid domain x = x(us)
Fluid-to-structure: apply boundary traction t = t(uf)
Write coupled system as M ˙

u = r(u) u =

uf

us

r(u) =
rf(uf; x(us))

rs(us; t(uf))

M =
M f

M s

Structure of linearized residual

∂r ∂u(u) =     ∂rf ∂uf ∂rf ∂x ∂x ∂us ∂rsf ∂t ∂t ∂uf ∂rs ∂us    

25 / 37

SLIDE 30

High-order partitioned FSI solver: IMEX Runge-Kutta1

Exploit linear dependence of structure residual (rs) on traction (t)

r(u) =

rf(uf; x(us)

rs(us; t(uf))

=
rsf · (t(uf) − ˜

t)

f(u)

+

rf(uf; x(us))

rs(us; ˜ t)

g(u)
Apply high-order implicit-explicit Runge-Kutta scheme to discretize

M ∂u ∂t = r(u) = f(u)

explicit

+ g(u)

implicit
Explicit Runge-Kutta scheme ˆ

c, ˆ A, ˆ b for f(u)

Diagonally implicit scheme c, A, b for g(u)

un = un−1 +

s

i=1

ˆ biˆ kn,i +

s

i=1

bikn,i Mkn,i = ∆tng

un−1 +

i−1

j=1

ˆ aij ˆ kn,j +

i

j=1

aijkn,j

M ˆ

kn,i = ∆tnf

un−1

i−1

j=1

ˆ aij ˆ kn,j +

i

j=1

aijkn,j

1[van Zuijlen and Bijl, 2005, Froehle and Persson, 2014]

26 / 37

SLIDE 31

High-order partitioned FSI solver: IMEX Runge-Kutta1

Exploit linear dependence of structure residual (rs) on traction (t)

r(u) =

rf(uf; x(us)

rs(us; t(uf))

=
rsf · (t(uf) − ˜

t)

f(u)

+

rf(uf; x(us))

rs(us; ˜ t)

g(u)
Structure of linearized implicit residual

∂g ∂u(u) =     ∂rf ∂uf (uf; x(us)) ∂rf ∂x (uf; x(us)) ∂x ∂us (us) ∂rs ∂us (us; ˜ t)    

Solve: (1) implicit structure, (2) implicit fluid, (3) explicit structure
Sequence of solves can be performed at linearized or nonlinear level
Due to choice of IMEX partition: no explicit fluid stages

1[van Zuijlen and Bijl, 2005, Froehle and Persson, 2014]

26 / 37

SLIDE 32

High-order partitioned FSI solver: IMEX Runge-Kutta

Define stage solutions

us

n,i = us n−1 + i

j=1

aijks

n,j + i−1

j=1

ˆ aij ˆ k

s n,j

uf

n,i = uf n−1 + i

j=1

aijkf

n,j

Define traction predictor as true traction at previous stage

˜ tn,i = t(un,i−1)

Solve for stage velocities (i = 1, . . . , s)

M sks

n,i = ∆tnrs(us n,i; ˜

tn,i) M fkf

n,i = ∆tnrf(uf n,i; x(us n,i))

M sˆ k

s n,i = ∆tnrsf(t(uf n,i) − ˜

tn,i)

Update state solution at new time

uf

n = uf n−1 + s

j=1

bjkf

n,j,

us

n = us n−1 + s

j=1

bjks

n,j + s

j=1

ˆ bj ˆ k

s n,j 27 / 37

SLIDE 33

Validation: benchmark pitching airfoil system

Simple FSI benchmark problem for studying the high-order accuracy of the

IMEX scheme

Rigid pitching/heaving NACA 0012 airfoil, torsional spring
Smooth heaving step y(t) prescribed, angle θ(t) measured

Setup Mach number

28 / 37

SLIDE 34

Validation: benchmark pitching airfoil system

Up to 5th order of convergence in time.
Similar accuracy as solving fully coupled system

0.0 0.5 1.0 1.5 2.0 T 0.00 0.01 0.02 0.03 0.04 0.05 0.06 theta nofluid 2e-1 1e-1 4e-2 2e-2 2e-3 2e-4

Angle θ(t) vs time t Entropy

10-2 10-1 100 Time step ∆t 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 Relative error in θ(t) 1 1 1 3 1 4 1 5

Weak Coupling ARK3 FC-ARK3 ARK4 FC-ARK4 ARK5 FC-ARK5

29 / 37

SLIDE 35

Validation: cantilever system

Standard FSI benchmark problem.
Elastic cantilever behind a square bluff body in incompressible flow.

5.5 14.0 12.0 1.0 1.0 4.0 0.06

Cantilever:

ρs = 100 kg/m3, νs = 0.35, E = 2.5 × 105 Pa.

Fluid & Flow:

ρf = 1.18 kg/m3, νf = 1.54 × 10−5 m2/s, vf = 0.513 m/s, Re = 333, Ma = 0.2.

Vortex shedding frequency: ∼ 6.3 Hz
Cantilever first mode: 3.03 Hz

30 / 37

SLIDE 36

Validation: cantilever system

Entropy

5 10 15 20 Time (s) 1.5 1.0 0.5 0.0 0.5 1.0 1.5 Vertical tip displacement (cm)

Tip displacement

Tip frequency:

f = 3.14 Hz (Literature: 2.98 – 3.25 Hz)

Tip displacement:

dmax = 1.09 cm (Literature: 0.95 – 1.25 cm)

31 / 37

SLIDE 37

Flow around Membrane, 3-D

Angle of attack 22.6◦, Reynolds number 2000.
Flexible structure prevents leading edge separation.

32 / 37

SLIDE 38

Adjoint equations for high-order partitioned IMEX FSI solver

Define

rf

n,i = rf(uf n,i; x(us n,i))

rs

n,i = rs(us n,i; ˜

tn,i)

Final condition for state Lagrange multipliers (F is quantity of interest)

λf

Nt =

∂F ∂uf

Nt T

, λs

Nt =

∂F ∂us

Nt T

Solve for stage Lagrange multipliers (j = s, . . . , 1)
Explicit structure stage

M sT ˆ κs

n,j =

∂F ∂ˆ k

s n,j T

+ˆ bjλs

n + ∆tn s

i=j+1

ˆ aij ∂rf

n,i

∂us

T

κf

n,i + ∆tn s

i=j+1

ˆ aij ∂rs

n,i

∂us

T

κs

n,i

Implicit fluid stage

M f T κf

n,j =

∂F ∂kf

n,j T

+ bjλf

n + ∆tn s

i=j

aij ∂rf

n,i

∂uf

T

κf

n,i + ∆tn s

i=j+1

aij ∂˜ tn,i ∂uf

T

rsf T κs

n,i

− ∆tn

s

i=j

aij ∂tn,i ∂uf

T

rsf T ˆ κs

n,i + ∆tn s

i=j+1

aij ∂˜ tn,i ∂uf

T

rsf T ˆ κs

n,i

Implicit structure stage

M sT κs

n,j =

∂F ∂ks

n,j T

+ bjλs

n + ∆tn s

i=j

aij ∂rf

n,i

∂us

T

κf

n,i + ∆tn s

i=j

aij ∂rs

n,i

∂us

T

κs

n,i

33 / 37

SLIDE 39

Adjoint equations for high-order partitioned IMEX FSI solver

Update state Lagrange multipliers at new time

λf

n−1 = λf n +

∂F ∂uf

n−1 T

+ ∆tn

s

i=1

∂rf

n,i

∂uf

T

κf

n,i + ∆tn s

i=1

∂˜ tn,i ∂uf

T

rsf

n,i T κs n,i

+ ∆tn

s

i=1

∂˜ tn,i ∂uf − ∂tn,i ∂uf T rsf

n,i T ˆ

κs

n,i

λs

n−1 = λs n +

∂F ∂us

n−1 T

+ ∆tn

s

i=1

∂rf

n,i

∂us

T

κf

n,i + ∆tn s

i=1

∂rs

n,i

∂us

T

κs

n,i

Reconstruct total derivative of quantity of interest F as

dF dµ = ∂F ∂µ + λf

T ∂¯

uf ∂µ + λs

T ∂¯

us ∂µ −

Nt

n=0

∆tn

s

i=1

κf

n,i T ∂rf n,i

∂µ −

Nt

n=0

∆tn

s

i=1

κs

n,i T ∂rs n,i

∂µ −

Nt

n=0

∆tn

s

i=1

ˆ κs,T

n,i

∂rsf

n,i

∂µ

34 / 37

SLIDE 40

Optimal energy harvesting from foil-damper system

Goal: Maximize energy harvested from foil-damper system maximize

µ

1 T T (c˙ h2(us) − Mz(uf) ˙ θ(µ, t)) dt

Fluid: Isentropic Navier-Stokes on deforming domain (ALE)
Structure: Force balance in y-direction between foil and damper
Motion driven by imposed θ(µ, t) = µ1 cos(2πft); µ1 ∈ (−45◦, 45◦)

c θ(µ, t) h(us) µ∗

1 = 45◦ 35 / 37

SLIDE 41

References I

Alexander, R. (1977). Diagonally implicit Runge-Kutta methods for stiff o.d.e.’s. SIAM J. Numer. Anal., 14(6):1006–1021. Froehle, B. and Persson, P.-O. (2014). A high-order discontinuous galerkin method for fluid–structure interaction with efficient implicit–explicit time stepping. Journal of Computational Physics, 272:455–470. Peraire, J. and Persson, P.-O. (2008). The Compact Discontinuous Galerkin (CDG) method for elliptic problems. SIAM Journal on Scientific Computing, 30(4):1806–1824. Persson, P.-O., Bonet, J., and Peraire, J. (2009). Discontinuous galerkin solution of the navier–stokes equations on deformable domains. Computer Methods in Applied Mechanics and Engineering, 198(17):1585–1595.

36 / 37

SLIDE 42

References II

van Zuijlen, A. H. and Bijl, H. (2005). Implicit and explicit higher order time integration schemes for structural dynamics and fluid-structure interaction computations. Computers & structures, 83(2):93–105. Zahr, M. J. and Persson, P.-O. (2016). An adjoint method for a high-order discretization of deforming domain conservation laws for optimization of flow problems. Journal of Computational Physics. Zahr, M. J., Persson, P.-O., and Wilkening, J. (2016). A fully discrete adjoint method for optimization of flow problems on deforming domains with time-periodicity constraints. Computers & Fluids.

37 / 37