Accurate Methods for Computing High Mach Number Reactive Flow - - PowerPoint PPT Presentation

▶
accurate methods for computing high mach number reactive
SMART_READER_LITE
LIVE PREVIEW

Accurate Methods for Computing High Mach Number Reactive Flow - - PowerPoint PPT Presentation

Mar 02, 2023 302 likes •1.43k views

Accurate Methods for Computing High Mach Number Reactive Flow Matthew J. Zahr and Joseph M. Powers Department of Aerospace and Mechanical Engineering University of Notre Dame, USA Lawrence Livermore National Laboratory Livermore, California

slide-1
SLIDE 1

Accurate Methods for Computing High Mach Number Reactive Flow

Matthew J. Zahr and Joseph M. Powers

Department of Aerospace and Mechanical Engineering University of Notre Dame, USA

Lawrence Livermore National Laboratory Livermore, California 24 June 2020

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 1 / 62

slide-2
SLIDE 2

Motivation

We investigate numerical methods to accurately simulate high Mach number flow with shocks and exothermic reaction: Shock-capturing: (WENO, Roe, etc), less than first order convergence, Shock-fitting: high accuracy and potential high order convergence; limited to simple topologies and algorithmically complicated, Wavelet adaptive methods: high accuracy, appropriate for multiscale, algorithmically complicated, Implicit shock-tracking: high accuracy and potential high order convergence, finite element based. We find, relative to WENO, implicit shock tracking yields remarkable improvement in accuracy on a standard verification problem in oblique detonation.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 2 / 62

slide-3
SLIDE 3

WENO shock-capturing: looks good, but converges at < O(∆x)

1,986 journal articles on WENO (Shu, et al, 1996), cited 35,489 times! WENO shock-capturing “looks” better than Lax-Wendroff, etc. shock-capturing in the “picture norm”. But, WENO converges at O(∆x1−1/(r+1)), where r is the rate of convergence for smooth problems. Our r = 5 WENO converges at O(∆x5/6) for shocks. Henrick, Aslam, Powers, Journal of Computational Physics, 2005.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 3 / 62

slide-4
SLIDE 4

WENO with shock fitting: looks good and converges at O(∆x5)

Algorithmically complicated shock-fitting requires exact solution of shock jump conditions and mapping of the domain to shock-attached coordinates. It enables high accuracy and convergence rates. Our r = 5 WENO and shock-fitting converges at O(∆x4.907) for a problem with shock and reaction. Henrick, Aslam, Powers, Journal of Computational Physics, 2006.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 4 / 62

slide-5
SLIDE 5

Pseudo-spectral with shock fitting: looks good and converges spectrally!

Shock-fitting may be coupled with a pseudo-spectral method for very high accuracy and convergence rates. Solution here for an inert blunt body re-entry problem at M = 3.5. Algorithmically complicated! Karhunen-Loéve-based reduced-order model for shape

  • ptimization
  • Brooks and Powers, Journal of Com-

putational Physics, 2004.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 5 / 62

slide-6
SLIDE 6

Physical diffusion can eliminate carbuncle instability, but must resolve thin zones

Roe-based shock-capturing may induce non-physical carbuncle instabilities for Euler equation solution. Inclusion of physical viscosity in Navier-Stokes solution stabilizes the flow. Solution here from OpenFOAM for an inert blunt body re-entry problem at M = 5.73. Resolution of thin zones is expensive!

Pressure (Pa) 250 230 210 190 170 150 130 110 90 70 50 a = 1.5 x 10 -4 m 5.377 x 10 -4 m Pressure (Pa) 250 230 210 190 170 150 130 110 90 70 50 a = 1.5 x 10 -4 m 5.377 x 10 -4 m

Powers, Bruns, Jemcov, AIAA- 2015-0579, 2015.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 6 / 62

slide-7
SLIDE 7

Adaptive wavelet method can capture viscous shock and detailed kinetics reaction

Wavelet basis can efficiently capture multiscale signals. Solutions must be continuous, so shocks require viscosity. Paolucci, Powers, et al., 2001, initiated this method for realistic combustion flows with detailed kinetics. Nonlinear chaotic dynamics resolved. Multidimensional viscous detonations in H2 − air resolved. Algorithmically complicated! Romick, Aslam, Powers, Jour- nal of Fluid Mechanics; Romick, Ph.D. Dissertation, 2015.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 7 / 62

slide-8
SLIDE 8

Implicit shock tracking: a relatively new method

The new method of implicit shock tracking is a competitive and

  • ften highly advantageous method, as will be described shortly.

To quantify its advantage, we will test it against the widely used WENO shock-capturing method on a benchmark verification problem in two-dimensional steady detonation. Powers, J. M., and Aslam, T. D., 2006, “Exact Solution for Multidimensional Compressible Reactive Flow for Verifying Numerical Algorithms,” AIAA Journal, 44(2): 337-344. We develop this useful solution next.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 8 / 62

slide-9
SLIDE 9

Partial review of oblique detonation

Samaras, D. G., “Gas Dynamic Treatment of Exothermic and Endothermic Discontinuities,” Canadian Journal of Research A, Vol. 26, No. 1, 1948, pp. 1-21. Gross, R. A., “Oblique Detonation Waves,” AIAA Journal, Vol. 1, No. 5, 1963,

  • pp. 1225-1227.

Pratt, D. T., Humphrey, J. W., and Glenn, D. E., “Morphology of Standing Oblique Detonation Waves,” Journal of Propulsion and Power, Vol. 7, No. 5, 1991,

  • pp. 837-645.

Lee, R. S., “A Unified Analysis of Supersonic Nonequilibrium Flow over a Wedge: I. Vibrational Nonequilibrium,” AIAA Journal, Vol. 4, No. 1, 1966, pp. 30-37. Powers, J. M., and Stewart, D. S., “Approximate Solutions for Oblique Detonations in the Hypersonic Limit,” AIAA Journal, Vol. 30, No. 3, 1992, pp. 726-736. Grismer, M. J., and Powers, J. M., “Comparisons of Numerical Oblique Detonation Solutions with an Asymptotic Benchmark,” AIAA Journal, Vol. 30, No. 12, 1992,

  • pp. 2985-2987.

Powers, J. M., and Gonthier, K. A., “Reaction Zone Structures for Strong, Weak Overdriven, and Weak Underdriven Oblique Detonations,” Physics of Fluids A,

  • Vol. 4, No. 9, 1992, pp. 2082-2089.

Grismer, M. J., and Powers, J. M., “Numerical Predictions of Oblique Detonation Stability Boundaries,” Shock Waves, Vol. 6, No. 3, 1996, pp. 147-156.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 9 / 62

slide-10
SLIDE 10

Oblique detonation schematic

Straight shock. Curved wedge. Orthogonal coordinate system aligned with shock.

x y

X Y

β s t r a i g h t s h

  • c

k curved wedge

u1 u cos β

1

u sin β

1

p = p ρ = ρ T = T λ = 0

1 1 1 λ = 1

streamline reaction zone ~

Powers and Aslam, AIAA Journal, 2006.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 10 / 62

slide-11
SLIDE 11

Model: reactive Euler equations

two-dimensional, steady, inviscid, irrotational,

  • ne-step kinetics with zero activation energy,

calorically perfect ideal gases with identical molecular masses and specific heats

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 11 / 62

slide-12
SLIDE 12

Model: reactive Euler PDEs

∂ ∂X (ρU) + ∂ ∂Y (ρV ) = 0, ∂ ∂X

  • ρU 2 + p
  • + ∂

∂Y (ρUV ) = 0, ∂ ∂X (ρUV ) + ∂ ∂Y

  • ρV 2 + p
  • =

0, ∂ ∂X

  • ρU
  • e + 1

2(U 2 + V 2) + p ρ

  • + ∂

∂Y

  • ρV
  • e + 1

2(U 2 + V 2) + p ρ

  • =

0, ∂ ∂X (ρUλ) + ∂ ∂Y (ρV λ) = αρ(1 − λ)H(T − Ti), e = 1 γ − 1 p ρ − λq, p = ρRT.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 12 / 62

slide-13
SLIDE 13

Model reductions: PDEs→ODEs

Assume no Y variation, so d dX (ρU) = 0, d dX

  • ρU 2 + p
  • =

0, d dX (ρUV ) = 0, d dX

  • ρU
  • e + 1

2(U 2 + V 2) + p ρ

  • =

0, d dX (ρUλ) = αρ(1 − λ)H(T − Ti).

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 13 / 62

slide-14
SLIDE 14

Model reductions: ODEs→DAEs

ρU = ρ1u1 sin β, ρU 2 + p = ρ1u2

1 sin2 β + p1,

V = u1 cos β, γ γ − 1 p ρ − λq + 1 2

  • U 2 + u2

1 cos2 β

  • =

γ γ − 1 p1 ρ1 + 1 2u2

1,

dλ dX = α1 − λ U H(T − Ti).

ZND reaction zone structure ODE supplemented with ex- tended Rankine-Hugoniot algebraic conditions.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 14 / 62

slide-15
SLIDE 15

Model reductions: inversion of algebraic relations

with M1 ≡ M1 sin β,

ρ(λ) = ρ1(γ + 1)M2

1

1 + γM2

1 ±

  • (1 + γM2

1)2 − (γ + 1)M2 1

  • 2 + γ−1

γ 2λq RT1 + (γ − 1)M2 1

, U(λ) = ρ1u1 sin β ρ(λ) , p(λ) = p1 + ρ2

1u2 1 sin2 β

1 ρ1 − 1 ρ(λ)

  • ,

T(λ) = p1 ρ(λ)R + ρ2

1u2 1 sin2 β

ρ(λ)R 1 ρ1 − 1 ρ(λ)

  • ,

q ≤ γRT1(M2

1 − 1)2

2(γ2 − 1)M2

1

,

CJ limitation. + shocked; − unshocked. Take the shocked branch.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 15 / 62

slide-16
SLIDE 16

Reaction zone structure solution

dλ dX = α ρ1u1 sin β ρ(λ)(1 − λ), λ(0) = 0,

X(λ) = a1       2a3

  • 1 − a4λ − 1
  • + ln

     

  • 1

1 − λ a2    

  • 1 −
  • 1−a4λ

1−a4

1 +

  • 1

1−a4

  • 1 +
  • 1−a4λ

1−a4

1 −

  • 1

1−a4

  

a3

√1−a4             ,

a1 = 1 (γ + 1)M1 √γRT1 α , a2 = 1 + γM2

1,

a3 = M2

1 − 1,

a4 = 2 M2

1

  • M2

1 − 1

2 γ2 − 1 γ q RT1 .

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 16 / 62

slide-17
SLIDE 17

Parametric values

Independent Parameter Units Value R J/kg/K 287 α 1/s 1000 β rad π/4 γ

  • 6/5

T1 K 300 M1

  • 3

ρ1 kg/m3 1 q J/kg 300000 Ti K 131300/363

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 17 / 62

slide-18
SLIDE 18

Reaction zone structure normal to shock

  • 1 0

1 2 3 4 X (m) 1 2 3 4 ρ (kg/ m )

3

  • 1 0

1 2 3 4 X (m) 300 400 500 600 T(K)

  • 1 0

1 2 3 4 X (m) 1 2 3 MX

  • 1 0

1 2 3 4 X (m) 0.2 0.4 0.6 0.8 1 λ

a) b) c) d)

200 0.0 LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 18 / 62

slide-19
SLIDE 19

Exact solution: streamlines

Curved streamlines identical to wedge contour. Streamline curvature approaches zero as reaction completes.

2 4 1 2 3 4 1 3 x (m) y (m) λ ~ . 9 straight oblique shock, β = π / 4 curved wedge

Powers and Aslam, AIAA Journal, 2006.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 19 / 62

slide-20
SLIDE 20

Reaction zone structure normal to shock

  • 1 0

1 2 3 4 X (m) 1 2 3 4 ρ (kg/ m )

3

  • 1 0

1 2 3 4 X (m) 300 400 500 600 T(K)

  • 1 0

1 2 3 4 X (m) 1 2 3 MX

  • 1 0

1 2 3 4 X (m) 0.2 0.4 0.6 0.8 1 λ

a) b) c) d)

200 0.0 LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 20 / 62

slide-21
SLIDE 21

Verification of WENO shock capturing

Algorithm of Xu, Aslam, and Stewart, 1997, CTM. Uniform Cartesian grid. Embedded internal boundary with level set representation. Nominally fifth order weighted essentially non-oscillatory (WENO) discretization. Non-decomposition based Lax-Friedrichs solver. Third order Runge-Kutta time integration.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 21 / 62

slide-22
SLIDE 22

Exact versus WENO solution

256 × 256 uniform numerical grid. good agreement in picture norm. numerical solution stable.

0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 2.0

x (m) y (m) exact numerical

ρ = 2.0 kg/m3 ρ = 2.9 kg/m3 ρ = 2.6 kg/m3 ρ = 2.3 kg/m3

Powers and Aslam, AIAA Journal, 2006.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 22 / 62

slide-23
SLIDE 23

Iterative convergence to steady state: various grids

Coarse grids relax quickly; fine grids relax slowly. All grids iteratively converge to steady state. Iterative convergence is distinct from grid convergence.

0.01 0.1 1 10 0.002 0.004 0.006 0.008 0.01 64 x 64 128 x 128 256 x 256 512 x 512 1024 x 1024 t (s) L (kg/m)

1

Powers and Aslam, AIAA Journal, 2006.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 23 / 62

slide-24
SLIDE 24

Grid convergence

Convergence rate: O

  • ∆x0.779

. Both shock capturing and embedded boundary induce the low convergence rate.

0.01 0.1 1 0.001 0.01 0.1 ∆x (m) L (kg/m)

1

Powers and Aslam, AIAA Journal, 2006.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 24 / 62

slide-25
SLIDE 25

Conclusions

Modeling flows with thin zones and surfaces of discontinuity is challenging. Most common shock-capturing methods converge at less than first

  • rder.

Heroic measures are often needed to achieve high order convergence Implicit shock tracking......

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 25 / 62

slide-26
SLIDE 26

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-27
SLIDE 27

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-28
SLIDE 28

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-29
SLIDE 29

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-30
SLIDE 30

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-31
SLIDE 31

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-32
SLIDE 32

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-33
SLIDE 33

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-34
SLIDE 34

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-35
SLIDE 35

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-36
SLIDE 36

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-37
SLIDE 37

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-38
SLIDE 38

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-39
SLIDE 39

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-40
SLIDE 40

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-41
SLIDE 41

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-42
SLIDE 42

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-43
SLIDE 43

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-44
SLIDE 44

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-45
SLIDE 45

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-46
SLIDE 46

High-order implicit tracking: reacting inviscid flow

Original mesh: 100 triangular elements (unstructured)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 26 / 62

slide-47
SLIDE 47

Reacting inviscid flow: density (ρ)

p = q = 1 p = q = 2 p = q = 3 p: polynomial degree of solution q: polynomial degree of mesh

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 27 / 62

slide-48
SLIDE 48

Reacting inviscid flow: density (ρ)

p = q = 1 p = q = 2 p = q = 3 p: polynomial degree of solution q: polynomial degree of mesh

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 27 / 62

slide-49
SLIDE 49

Reacting inviscid flow: reaction progress (λ)

p = q = 1 p = q = 2 p = q = 3 p: polynomial degree of solution q: polynomial degree of mesh

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 28 / 62

slide-50
SLIDE 50

Isolines of density indicate excellent agreement with analytical solution on coarse (100 element) mesh for p > 1

Analytical ( ) p = 1 ( ) p = 2 ( ) p = 3 ( ) −0.20 0.5 1 1.5 1.9 0.5 1 1.5 2 x1 x2 ρ = 2 ρ = 2.3 ρ = 2.6 ρ = 2.9 shock

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 29 / 62

slide-51
SLIDE 51

Implicit shock tracking recovers optimal O(hp+1) convergence rates; compares favorably to WENO

10−2 10−1 10−4 10−3 10−2 10−1 mesh size parameter (h) L1 error (ρ) WENO ( ), p = 1 ( ), p = 2 ( ), p = 3 ( )

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 30 / 62

slide-52
SLIDE 52

Implicit shock tracking

Goal: Align element faces with (unknown) discontinuities to perfectly capture them and approximate smooth regions to high-order DG discretization: inter-element jumps, high-order Discontinuity-aligned mesh is the solution of an optimization problem constrained by the discrete PDE = ⇒ implicit tracking Simultaneously converge solution and mesh to ensure solution of PDE never required on non-aligned mesh

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 31 / 62

slide-53
SLIDE 53

Implicit shock tracking

Goal: Align element faces with (unknown) discontinuities to perfectly capture them and approximate smooth regions to high-order DG discretization: inter-element jumps, high-order Discontinuity-aligned mesh is the solution of an optimization problem constrained by the discrete PDE = ⇒ implicit tracking Simultaneously converge solution and mesh to ensure solution of PDE never required on non-aligned mesh

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 31 / 62

slide-54
SLIDE 54

Discontinuous Galerkin discretization of conservation law

Inviscid conservation law: ∇ · F(U) = 0 in Ω Test Vh,p′ and trial Vh,p spaces, where h is the mesh size and p/p′ is the polynomial degree, to define the finite-dimensional DG residual: rK

h,p′(Uh,p) :=

  • ∂K

ψ+

h,p′ · H(U + h,p, U− h,p, n) dS −

  • K

F(Uh,p) : ∇ψh,p′ dV Introduce basis for polynomial spaces to obtain discrete residuals r(u, x) (p′ = p), R(u, x) (p′ = p + 1), where u is the discrete state vector and x are the coordinates of the mesh

  • nodes. The “standard” DG equations are: r(u, x) = 0.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 32 / 62

slide-55
SLIDE 55

Discontinuous Galerkin discretization of conservation law

Inviscid conservation law: ∇ · F(U) = 0 in Ω Test Vh,p′ and trial Vh,p spaces, where h is the mesh size and p/p′ is the polynomial degree, to define the finite-dimensional DG residual: rK

h,p′(Uh,p) :=

  • ∂K

ψ+

h,p′ · H(U + h,p, U− h,p, n) dS −

  • K

F(Uh,p) : ∇ψh,p′ dV Introduce basis for polynomial spaces to obtain discrete residuals r(u, x) (p′ = p), R(u, x) (p′ = p + 1), where u is the discrete state vector and x are the coordinates of the mesh

  • nodes. The “standard” DG equations are: r(u, x) = 0.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 32 / 62

slide-56
SLIDE 56

Discontinuous Galerkin discretization of conservation law

Inviscid conservation law: ∇ · F(U) = 0 in Ω Test Vh,p′ and trial Vh,p spaces, where h is the mesh size and p/p′ is the polynomial degree, to define the finite-dimensional DG residual: rK

h,p′(Uh,p) :=

  • ∂K

ψ+

h,p′ · H(U + h,p, U− h,p, n) dS −

  • K

F(Uh,p) : ∇ψh,p′ dV Introduce basis for polynomial spaces to obtain discrete residuals r(u, x) (p′ = p), R(u, x) (p′ = p + 1), where u is the discrete state vector and x are the coordinates of the mesh

  • nodes. The “standard” DG equations are: r(u, x) = 0.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 32 / 62

slide-57
SLIDE 57

Implicit shock tracking: constrained optimization

We formulate the problem of tracking discontinuities with the mesh as the solution of an optimization problem constrained by the discrete PDE (DG discretization) minimize

u,x

f(u, x) subject to r(u, x) = 0. We propose an objective function that balances the tracking objective with maintaining a high-quality mesh f(u, x) = 1 2R(u, x)T R(u, x) + κ2 2 Rmsh(x)T Rmsh(x).

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 33 / 62

slide-58
SLIDE 58

Implicit shock tracking: constrained optimization

We formulate the problem of tracking discontinuities with the mesh as the solution of an optimization problem constrained by the discrete PDE (DG discretization) minimize

u,x

f(u, x) subject to r(u, x) = 0. We propose an objective function that balances the tracking objective with maintaining a high-quality mesh f(u, x) = 1 2R(u, x)T R(u, x) + κ2 2 Rmsh(x)T Rmsh(x).

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 33 / 62

slide-59
SLIDE 59

Implicit shock tracking: SQP solver

Define z = (u, x) and use interchangeably. To solve the optimization problem, we define a sequence {zk} updated as zk+1 = zk + αk∆zk. The step direction ∆zk is defined as the solution of the quadratic program (QP) approximation of the tracking problem centered at zk minimize

∆z∈RNz

gz(zk)T ∆z + 1 2∆zT B(zk, ˆ λ(zk))∆z subject to r(zk) + Jz(zk)∆z = 0, where gz(z) = ∂f ∂z (z)T , Jz(z) = ∂r ∂z (z), Bz(z, λ) ≈ ∂2L ∂z∂z (z, λ), L(z, λ) = f(z) − λT r(z), ˆ λ(z) = ∂r ∂u(z)−T ∂f ∂u(z)T .

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 34 / 62

slide-60
SLIDE 60

Implicit shock tracking: SQP solver

Define z = (u, x) and use interchangeably. To solve the optimization problem, we define a sequence {zk} updated as zk+1 = zk + αk∆zk. The step direction ∆zk is defined as the solution of the quadratic program (QP) approximation of the tracking problem centered at zk minimize

∆z∈RNz

gz(zk)T ∆z + 1 2∆zT B(zk, ˆ λ(zk))∆z subject to r(zk) + Jz(zk)∆z = 0, where gz(z) = ∂f ∂z (z)T , Jz(z) = ∂r ∂z (z), Bz(z, λ) ≈ ∂2L ∂z∂z (z, λ), L(z, λ) = f(z) − λT r(z), ˆ λ(z) = ∂r ∂u(z)−T ∂f ∂u(z)T .

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 34 / 62

slide-61
SLIDE 61

Practical considerations: element collapse

Despite measures to keep mesh well-conditioned, best option may be to remove element from the mesh: tag elements for removal based on volume, collapse shortest edge, trivial DG solution transfer Only meshing operation required by implicit shock tracking Simplices: arbitrary dimensions, arbitrary polynomial degree

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 35 / 62

slide-62
SLIDE 62

Practical considerations: element collapse

Despite measures to keep mesh well-conditioned, best option may be to remove element from the mesh: tag elements for removal based on volume, collapse shortest edge, trivial DG solution transfer Only meshing operation required by implicit shock tracking Simplices: arbitrary dimensions, arbitrary polynomial degree

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 35 / 62

slide-63
SLIDE 63

Linear advection, straight shock

p = 0 space for solution, q = 1 space for mesh L1 error = 3.84 × 10−11

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 36 / 62

slide-64
SLIDE 64

Newton-like convergence when solution lies in DG space

Linear advection with straight shock 2 4 6 8 10 10−16 10−12 10−8 10−4 100 Iteration Solver convergence r(z) ( ), R(z) ( ),

  • ∇xL(z, ˆ

λ(z))

  • (

)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 37 / 62

slide-65
SLIDE 65

Linear advection, trigonometric shock

p = 0 space for solution, q = 2 space for mesh L1 error = 1.15 × 10−3

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 38 / 62

slide-66
SLIDE 66

Linear advection, trigonometric shock, 3D

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 39 / 62

slide-67
SLIDE 67

Linear advection, trigonometric shock, 3D

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 39 / 62

slide-68
SLIDE 68

Inviscid Burgers’ equation: space-time formulation

p = 0 space for solution, q = 3 space for mesh

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 40 / 62

slide-69
SLIDE 69

Temporal snapshots show discontinuity perfectly represented, smooth solution well-approximated elsewhere

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 x U(x, t) Temporal snapshots of inviscid Burgers’ equation from p = 4 implicit tracking simulation

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 41 / 62

slide-70
SLIDE 70

Supersonic flow past NACA0012 airfoil (M = 1.5)

Initialization p = 1 tracking eH = 1.30 × 10−3 p = 2 tracking eH = 6.73 × 10−5

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 42 / 62

slide-71
SLIDE 71

Supersonic flow past NACA0012 airfoil (M = 1.5)

Initialization p = 1 tracking eH = 1.30 × 10−3 p = 2 tracking eH = 6.73 × 10−5

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 42 / 62

slide-72
SLIDE 72

Spatial slices show both discontinuities are tracked and solution well-approximated elsewhere

−0.5 0 0.5 1 1.5 1 1.5 2 2.5 x ρ −0.5 0 0.5 1 1.5 0.8 1 1.2 1.4 1.6 x M −0.5 0 0.5 1 1.5 1 1.5 2 x P Spatial slice of supersonic flow past NACA0012 airfoil (M = 1.5) from p = 3 implicit tracking simulation

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 43 / 62

slide-73
SLIDE 73

High-order, implicit shock tracking

Primary benefit: highly accurate solutions on coarse meshes, recover optimal convergence rates Traditional barrier to tracking (explicitly meshing unknown discontinuity surface) replaced with solving constrained

  • ptimization problem

Future extensions: 3D, viscous, time-dependent, model reduction, hypersonics Zahr, Shi, Persson, Journal of Computational Physics, 2020.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 44 / 62

slide-74
SLIDE 74

Implicit shock tracking for hypersonics applications

Strong local features: shock waves and contact discontinuities Improved solver robustness (time-stepping, vanishing viscosity) Viscous effects: thin boundary layer, critical to predict heating Extension to viscous problems straightforward Expect r-adaptivity near boundary layer to improve approximation Interacting features: shock/shock and shock/boundary layer Handled by optimization formulation and element collapse Unsteady: prediction of combustion stability Slab-based space-time formulation (moving features) Method-of-lines-based formulation (stationary features) Multi-scale/physics: turbulence, ablation, conjugate heat transfer Clip objective function to avoid tracking all turbulent features Integrate in partitioned multiphysics setting

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 45 / 62

slide-75
SLIDE 75

Reduction of parametrized discontinuities

Fundamental issue: linear subspace approximation ill-suited for advection- dominated features (slowly decay Kolmogorov n-width) apply parameter-dependent domain mapping to align features use linear subspace in reference domain to reduce dimension push forward to physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 46 / 62

slide-76
SLIDE 76

Reduction of parametrized discontinuities

Fundamental issue: linear subspace approximation ill-suited for advection- dominated features (slowly decay Kolmogorov n-width) apply parameter-dependent domain mapping to align features use linear subspace in reference domain to reduce dimension push forward to physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 46 / 62

slide-77
SLIDE 77

Reduction of parametrized discontinuities

Fundamental issue: linear subspace approximation ill-suited for advection- dominated features (slowly decay Kolmogorov n-width) apply parameter-dependent domain mapping to align features use linear subspace in reference domain to reduce dimension push forward to physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 46 / 62

slide-78
SLIDE 78

Reduction of parametrized discontinuities

Fundamental issue: linear subspace approximation ill-suited for advection- dominated features (slowly decay Kolmogorov n-width) apply parameter-dependent domain mapping to align features use linear subspace in reference domain to reduce dimension push forward to physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 46 / 62

slide-79
SLIDE 79

Reduction of parametrized discontinuities

Fundamental issue: linear subspace approximation ill-suited for advection- dominated features (slowly decay Kolmogorov n-width) Proposed solution apply parameter-dependent domain mapping to align features use linear subspace in reference domain to reduce dimension push forward to physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 46 / 62

slide-80
SLIDE 80

Reduction of parametrized discontinuities

Fundamental issue: linear subspace approximation ill-suited for advection- dominated features (slowly decay Kolmogorov n-width) Proposed solution apply parameter-dependent domain mapping to align features use linear subspace in reference domain to reduce dimension push forward to physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 46 / 62

slide-81
SLIDE 81

Reduction of parametrized discontinuities

Fundamental issue: linear subspace approximation ill-suited for advection- dominated features (slowly decay Kolmogorov n-width) Proposed solution apply parameter-dependent domain mapping to align features use linear subspace in reference domain to reduce dimension push forward to physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 46 / 62

slide-82
SLIDE 82

Numerical experiment: parametrized linear advection

The governing parametrized conservation law is ∇ · (βµU) + σµU = fµ in Ω := (0, 1)2 U = ˆ Uµ

  • n Γi := {x ∈ ∂Ω | βµ · n < 0}

µ1 : shock angle (βµ) µ2 : magnitude of source term (σµ) µ3 : magnitude of boundary condition ( ˆ Uµ)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 47 / 62

slide-83
SLIDE 83

Optimization-based alignment of discontinuities in reference domain (similar to implicit shock tracking)

Snapshot 1 Reference domain Physical domain

The blue line indicates the true orientation of the discontinuity in the physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 48 / 62

slide-84
SLIDE 84

Optimization-based alignment of discontinuities in reference domain (similar to implicit shock tracking)

Snapshot 5 Reference domain Physical domain

The blue line indicates the true orientation of the discontinuity in the physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 48 / 62

slide-85
SLIDE 85

Optimization-based alignment of discontinuities in reference domain (similar to implicit shock tracking)

Snapshot 6 Reference domain Physical domain

The blue line indicates the true orientation of the discontinuity in the physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 48 / 62

slide-86
SLIDE 86

With feature alignment, POD modes do not have to resolve moving discontinuity since it is (mostly) fixed

Mode 1 Mode 2 Mode 3

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 49 / 62

slide-87
SLIDE 87

Without feature alignment, POD modes must resolve moving discontinuity, which leads to slowly decaying Kolmogorov n-width

Mode 1 Mode 2 Mode 3

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 50 / 62

slide-88
SLIDE 88

Alignment framework leads to reduction in error for discontinuities not encountered during training

Test point 96 (discontinuity not in training) HDM ROM (aligned) ROM (non-aligned)

The blue line indicates the true orientation of the discontinuity in the physical domain

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 51 / 62

slide-89
SLIDE 89

Alignment framework leads to reduction in error for discontinuities not encountered during training

Testing set Maximum error Mean error discontinuity in training 0.92% 0.45% discontinuity not in training 26.5% 19.5% Classical minimum-residual model reduction Testing set Maximum error Mean error discontinuity in training 2.61% 1.22% discontinuity not in training 4.18% 3.31% Minimum-residual model reduction with feature alignment Parameters: 12 training, 125 testing

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 52 / 62

slide-90
SLIDE 90

Appendix

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 53 / 62

slide-91
SLIDE 91

Implicit shock tracking: SQP solver

The solution of the QP leads to the following linear system   Buu(zk, ˆ λ(zk)) Bux(zk, ˆ λ(zk)) Ju(zk)T Bux(zk, ˆ λ(zk))T Bxx(zk, ˆ λ(zk)) Jx(zk)T Ju(zk) Jx(zk)     ∆uk ∆xk ηk   = −   gu(zk) gx(zk) r(zk)   , where g(z) = ∂f ∂(z)T , J(z) = ∂r ∂(z), the approximate Hessian of the Lagrangian is partitioned as B△(z, λ) ≈ ∂2L ∂∂△(z, λ), and ηk are the Lagrange multipliers of the QP.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 54 / 62

slide-92
SLIDE 92

SQP solver: Levenberg-Marquardt Hessian approximation

The proposed objective function takes the form of a residual norm f(z) = 1 2 F (z)2

2 ,

F (z) =

  • R(z)

κRmsh(x)

  • and therefore the Hessian of the Lagrangian has the special structure

H(z, λ) = ∂F ∂z (z)T ∂F ∂z (z) + F i(z) ∂2F i ∂z∂z (z) − λi ∂2ri ∂z∂z (z) The last two terms are dropped: they are difficult to compute and neg- ligible if the F and λ are small (Gauss-Newton approximation) H(z, λ) ≈ ∂F ∂z (z)T ∂F ∂z (z) To guard against ill-conditioning, a regularization matrix is added Buu = ∂F ∂u

T ∂F

∂u , Bux = ∂F ∂u

T ∂F

∂x , Bxx = ∂F ∂x

T ∂F

∂x + γD.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 55 / 62

slide-93
SLIDE 93

SQP solver: Levenberg-Marquardt Hessian approximation

The proposed objective function takes the form of a residual norm f(z) = 1 2 F (z)2

2 ,

F (z) =

  • R(z)

κRmsh(x)

  • and therefore the Hessian of the Lagrangian has the special structure

H(z, λ) = ∂F ∂z (z)T ∂F ∂z (z) + F i(z) ∂2F i ∂z∂z (z) − λi ∂2ri ∂z∂z (z) The last two terms are dropped: they are difficult to compute and neg- ligible if the F and λ are small (Gauss-Newton approximation) H(z, λ) ≈ ∂F ∂z (z)T ∂F ∂z (z) To guard against ill-conditioning, a regularization matrix is added Buu = ∂F ∂u

T ∂F

∂u , Bux = ∂F ∂u

T ∂F

∂x , Bxx = ∂F ∂x

T ∂F

∂x + γD.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 55 / 62

slide-94
SLIDE 94

SQP solver: Levenberg-Marquardt Hessian approximation

The proposed objective function takes the form of a residual norm f(z) = 1 2 F (z)2

2 ,

F (z) =

  • R(z)

κRmsh(x)

  • and therefore the Hessian of the Lagrangian has the special structure

H(z, λ) = ∂F ∂z (z)T ∂F ∂z (z) + F i(z) ∂2F i ∂z∂z (z) − λi ∂2ri ∂z∂z (z) The last two terms are dropped: they are difficult to compute and neg- ligible if the F and λ are small (Gauss-Newton approximation) H(z, λ) ≈ ∂F ∂z (z)T ∂F ∂z (z) To guard against ill-conditioning, a regularization matrix is added Buu = ∂F ∂u

T ∂F

∂u , Bux = ∂F ∂u

T ∂F

∂x , Bxx = ∂F ∂x

T ∂F

∂x + γD.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 55 / 62

slide-95
SLIDE 95

Practical considerations: regularization matrix D

The mesh regularization matrix D is taken as the stiffness matrix of the linear elliptic PDE ∇ · (k∇vi) = 0 in Ω for i = 1, . . . , d. The coefficient is constant over each element and in- versely proportional to the element volume k(x) = min

K′∈Eh,q

|K′| |K| , x ∈ K for each element K in the mesh T : critical to maintain well-conditioned search directions for meshes where element size varies significantly.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 56 / 62

slide-96
SLIDE 96

Practical considerations: regularization matrix D

The mesh regularization matrix D is taken as the stiffness matrix of the linear elliptic PDE ∇ · (k∇vi) = 0 in Ω for i = 1, . . . , d. The coefficient is constant over each element and in- versely proportional to the element volume k(x) = min

K′∈Eh,q

|K′| |K| , x ∈ K for each element K in the mesh T : critical to maintain well-conditioned search directions for meshes where element size varies significantly.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 56 / 62

slide-97
SLIDE 97

SQP solver: step length (αk)

The step length, αk ∈ (0, 1], is selected using a backtracking line search to ensure sufficient decrease of a merit function ϕk : R → R ϕk(αk) ≤ ϕk(0) + cαkϕ′

k(0),

c ∈ (0, 1). We use the ℓ1 merit function ϕk(α) := f(zk + α∆zk) + µ r(zk + α∆zk)1 where µ >

  • ˆ

λ(zk)

  • ∞ because it is “exact”, i.e., any minimizer of the
  • riginal optimization problem is a minimizer of ϕk.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 57 / 62

slide-98
SLIDE 98

SQP solver: step length (αk)

The step length, αk ∈ (0, 1], is selected using a backtracking line search to ensure sufficient decrease of a merit function ϕk : R → R ϕk(αk) ≤ ϕk(0) + cαkϕ′

k(0),

c ∈ (0, 1). We use the ℓ1 merit function ϕk(α) := f(zk + α∆zk) + µ r(zk + α∆zk)1 where µ >

  • ˆ

λ(zk)

  • ∞ because it is “exact”, i.e., any minimizer of the
  • riginal optimization problem is a minimizer of ϕk.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 57 / 62

slide-99
SLIDE 99

Practical considerations: termination criteria

The termination criteria for the solver is based on the Karush-Kuhn- Tucker (KKT) conditions: z⋆ is a solution if there exist Lagrange multi- pliers λ⋆ such that ∇uL(z⋆, λ⋆) = 0, ∇xL(z⋆, λ⋆) = 0, r(z⋆) = 0 Our choice for the Lagrange multiplier estimate ˆ λ(z) ensures ∇uL(z, ˆ λ(z)) = 0 and therefore termination is based on the remaining KKT conditions

  • ∇xL(z, ˆ

λ(z))

  • < ǫ1,

r(z) < ǫ2, where ǫ1, ǫ2 > 0 are convergence tolerances.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 58 / 62

slide-100
SLIDE 100

Practical considerations: termination criteria

The termination criteria for the solver is based on the Karush-Kuhn- Tucker (KKT) conditions: z⋆ is a solution if there exist Lagrange multi- pliers λ⋆ such that ∇uL(z⋆, λ⋆) = 0, ∇xL(z⋆, λ⋆) = 0, r(z⋆) = 0 Our choice for the Lagrange multiplier estimate ˆ λ(z) ensures ∇uL(z, ˆ λ(z)) = 0 and therefore termination is based on the remaining KKT conditions

  • ∇xL(z, ˆ

λ(z))

  • < ǫ1,

r(z) < ǫ2, where ǫ1, ǫ2 > 0 are convergence tolerances.

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 58 / 62

slide-101
SLIDE 101

Practical considerations: initialization

Since using gradient-based optimization, good initialization of u and x is crucial – For p = 1 simulations, x0 comes from mesh generation without knowledge of the shock location and u0 is the DG(p = 0) solution – For p > 1 simulations, x0 and u0 are the p = 1 tracking solution Reference mesh, p = 0 DG solution

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 59 / 62

slide-102
SLIDE 102

Practical considerations: initialization

Since using gradient-based optimization, good initialization of u and x is crucial – For p = 1 simulations, x0 comes from mesh generation without knowledge of the shock location and u0 is the DG(p = 0) solution – For p > 1 simulations, x0 and u0 are the p = 1 tracking solution p = 1 tracking solution

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 59 / 62

slide-103
SLIDE 103

Practical considerations: initialization

Since using gradient-based optimization, good initialization of u and x is crucial – For p = 1 simulations, x0 comes from mesh generation without knowledge of the shock location and u0 is the DG(p = 0) solution – For p > 1 simulations, x0 and u0 are the p = 1 tracking solution p = 2 tracking solution

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 59 / 62

slide-104
SLIDE 104

Choice of numerical flux function for implicit tracking

Shock tracking places strict requirements on numerical flux function since inter-element jumps will not tend to zero on shock surface

i

consistent [DG stability]: H(U, U, n) = F(U)n

ii

conservative [DG stability]: H(U, U′, n) = −H(U ′, U, −n)

iii

preservation of Rankine-Hugoniot [conservation at discontinuities] F(U +)n = F(U −)n = ⇒ H(U +, U−, n) = F(U +)n = F(U −)n

iv

smoothness w.r.t. n [gradient-based optimization]

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 60 / 62

slide-105
SLIDE 105

Choice of numerical flux function for implicit tracking

Shock tracking places strict requirements on numerical flux function since inter-element jumps will not tend to zero on shock surface

i

consistent [DG stability]: H(U, U, n) = F(U)n

ii

conservative [DG stability]: H(U, U′, n) = −H(U ′, U, −n)

iii

preservation of Rankine-Hugoniot [conservation at discontinuities] F(U +)n = F(U −)n = ⇒ H(U +, U−, n) = F(U +)n = F(U −)n

iv

smoothness w.r.t. n [gradient-based optimization]

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 60 / 62

slide-106
SLIDE 106

Choice of numerical flux function for implicit tracking

Shock tracking places strict requirements on numerical flux function since inter-element jumps will not tend to zero on shock surface

i

consistent [DG stability]: H(U, U, n) = F(U)n

ii

conservative [DG stability]: H(U, U′, n) = −H(U ′, U, −n)

iii

preservation of Rankine-Hugoniot [conservation at discontinuities] F(U +)n = F(U −)n = ⇒ H(U +, U−, n) = F(U +)n = F(U −)n

iv

smoothness w.r.t. n [gradient-based optimization]

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 60 / 62

slide-107
SLIDE 107

Choice of numerical flux function for implicit tracking

Shock tracking places strict requirements on numerical flux function since inter-element jumps will not tend to zero on shock surface

i

consistent [DG stability]: H(U, U, n) = F(U)n

ii

conservative [DG stability]: H(U, U′, n) = −H(U ′, U, −n)

iii

preservation of Rankine-Hugoniot [conservation at discontinuities] F(U +)n = F(U −)n = ⇒ H(U +, U−, n) = F(U +)n = F(U −)n

iv

smoothness w.r.t. n [gradient-based optimization]

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 60 / 62

slide-108
SLIDE 108

Choice of numerical flux function for implicit tracking

Shock tracking places strict requirements on numerical flux function since inter-element jumps will not tend to zero on shock surface

i

consistent [DG stability]: H(U, U, n) = F(U)n

ii

conservative [DG stability]: H(U, U′, n) = −H(U ′, U, −n)

iii

preservation of Rankine-Hugoniot [conservation at discontinuities] F(U +)n = F(U −)n = ⇒ H(U +, U−, n) = F(U +)n = F(U −)n

iv

smoothness w.r.t. n [gradient-based optimization] Remark: Difficult to satisfy (i)-(iv); select standard numerical flux that satisfies (i)-(iii) and smooth any nonsmooth/discontinuous features

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 60 / 62

slide-109
SLIDE 109

Example: Upwind numerical flux for linear advection

The upwind numerical flux for linear advection (Fadv(U) = UβT ) is: Hup(U +, U−, n) = (β · n)

  • U + · H(β · n) + U − · (1 − H(β · n))
  • ,

which satisfies (i)-(iii) but is not smooth where β · n = 0 ((iv) fails). To recover smoothness, we replace the Heaviside function H(x) with the smoothed heaviside function Ha(x) :=

  • 1 + e−2ax−1: a = 5 (

), a = 10 ( ), a = 30 ( ), a = ∞ ( ). −1 −0.5 0.5 1 0.5 1 x Ha(x)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 61 / 62

slide-110
SLIDE 110

Example: Upwind numerical flux for linear advection

The upwind numerical flux for linear advection (Fadv(U) = UβT ) is: Hup(U +, U−, n) = (β · n)

  • U + · H(β · n) + U − · (1 − H(β · n))
  • ,

which satisfies (i)-(iii) but is not smooth where β · n = 0 ((iv) fails). To recover smoothness, we replace the Heaviside function H(x) with the smoothed heaviside function Ha(x) :=

  • 1 + e−2ax−1: a = 5 (

), a = 10 ( ), a = 30 ( ), a = ∞ ( ). −1 −0.5 0.5 1 0.5 1 x Ha(x)

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 61 / 62

slide-111
SLIDE 111

Non-smooth numerical flux can cause SQP solver to fail

Linear advection with trigonometric shock (p = 2) 20 40 60 80 10−14 10−8 10−2 Iteration (k) r(zk) 20 40 60 80 10−9 10−5 10−1 Iteration (k)

  • ∇xL(zk, ˆ

λ(zk))

  • Convergence of solver with standard (non-smooth) upwind flux (

) and smoothed upwind flux (a = 10) ( )

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 62 / 62

slide-112
SLIDE 112

Non-smooth numerical flux can cause SQP solver to fail

Linear advection with trigonometric shock (p = 2) 20 40 60 80 10−14 10−8 10−2 Iteration (k) r(zk) 20 40 60 80 10−9 10−5 10−1 Iteration (k)

  • ∇xL(zk, ˆ

λ(zk))

  • Convergence of solver with standard (non-smooth) upwind flux (

) and smoothed upwind flux (a = 10) ( )

LLNL Virtual Seminar Accurate CFD for Shocks 24 June 2020 62 / 62