[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 VII International Conference on Coupled Problems in Science and Engineering Rhodes Island, Greece, June 12-14, 2017

† Luis W. Alvarez Postdoctoral Fellow

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

SLIDE 2

Optimization of unsteady fluid-structure systems

Discover energetically optimal flapping motions
understand biological systems, design Micro Aerial Vehicles (MAVs)
Design minimal weight structures and vehicles that will not flutter
Design energy harvesting mechanisms

c θ(µ, t) h(us)

2 / 23

SLIDE 3

Unsteady single physics 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

3 / 23

SLIDE 4

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

4 / 23

SLIDE 5

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

5 / 23

SLIDE 6

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

6 / 23

SLIDE 7

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

7 / 23

SLIDE 8

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]

8 / 23

SLIDE 9

Energetically optimal flapping vs. required thrust

Energy = 0.21935 Thrust = 0.0000 Energy = 3.00404 Thrust = 1.5000 Energy = 6.2869 Thrust = 2.5000 Optimal Tx = 0 Optimal Tx = 1.5 Optimal Tx = 2.5

9 / 23

SLIDE 10

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) 10 / 23

SLIDE 11

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

11 / 23

SLIDE 12

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]

12 / 23

SLIDE 13

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)
Solve: (1) implicit structure, (2) implicit fluid, (3) explicit structure
Due to choice of IMEX partition: no explicit fluid stages

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

12 / 23

SLIDE 14

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 13 / 23

SLIDE 15

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

14 / 23

SLIDE 16

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

15 / 23

SLIDE 17

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

16 / 23

SLIDE 18

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)

17 / 23

SLIDE 19

Flow around Membrane, 3-D

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

18 / 23

SLIDE 20

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

19 / 23

SLIDE 21

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

∂µ

20 / 23

SLIDE 22

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)

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

1 ≈ 45◦ 21 / 23

SLIDE 23

Summary and future work

Summary

Extended standard fully discrete adjoint framework to partitioned, high-order

multiphysics setting

Demonstrated on simple optimal energy harvesting model problem

Future work

Extend structure model to fully deformable model
High-order, energy conserving load transfer from fluid to structure
Handle discontinuities between fluid elements that arise from DG discretization
Study optimal 3D flapping with deformable wing

22 / 23

SLIDE 24

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.

23 / 23

SLIDE 25

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.

24 / 23