DG-KFVS schemes for convection-diffusion equations Praveen. C - - PowerPoint PPT Presentation

DG-KFVS schemes for convection-diffusion equations

• Praveen. C

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

• Dept. of Mathematics, IISER Pune

25 March, 2014

1 / 73

Compressible Navier-Stokes model

• Mass conservation equation

∂ρ ∂t + ∇ · (ρu) = 0

• Momentum equation

∂ ∂t(ρu) + ∇ · (pI + ρuu) = ∇ · τ

• Energy conservation

∂E ∂t + ∇ · (E + p)u = ∇ · (u · τ) + ∇ · q

2 / 73

Compressible Navier-Stokes model

• Equation of state: ideal gas

E = p γ − 1 + 1 2ρ|u|2, γ > 1

• Constitutive laws

τij = µ(∂iuj + ∂jui) − 2 3µ(∂kuk)δij, qi = −κ∂iT

• System of conservation laws

∂U ∂t + ∂ ∂xi Fi(U) + ∂ ∂xi Gi(U, ∇U) = 0

3 / 73

Flow over NACA0012 airfoil: M = 0.85, α = 1

Pressure

5

4 / 73

Flow through scramjet intake

Density

5 / 73

Model problems

• Linear convection equation

∂ρ ∂t + u∂ρ ∂x = 0

• Inviscid Burger’s equation

∂u ∂t + u∂u ∂x = 0

• Heat equation

∂u ∂t = µ∂2u ∂x2

6 / 73

Heat equation: Finite difference method

∂u ∂t = µ∂2u ∂x2 , x ∈ (0, 1) u(x, 0) = f(x), u(0, t) = u(1, t) = 0 Partition domain (0, 1) with grid points xj = jh, j = 0, 1, . . . , N + 1, uj(t) = u(xj, t) Approximate second derivative by finite difference uj−1 − 2uj + uj+1 h2 = ∂2u ∂x2 (xj) + O(h2) Semi-discrete scheme (system of ODE) duj dt = µuj−1 − 2uj + uj+1 h2 , j = 1, 2, . . . , N Central difference scheme, stencil (j − 1, j, j + 1) is symmetric

7 / 73

Heat equation: Finite element method

Find u(t) ∈ V = H1

0(0, 1) such that for all φ ∈ V

1 ∂u ∂t φdx + µ 1 ∂u ∂x ∂φ ∂xdx = 0 Main idea: Approximate V by a finite dimensional space Vh. Finite element method: Vh is made of piecewise polynomial functions. Partition (0, 1) into disjoint elements Ij = (xj− 1

2 , xj+ 1 2 ),

(0, 1) = ∪jIj Vh = V k

h = {φ ∈ C(0, 1) : φ(0) = φ(1) = 0, φ|Ij ∈ Pk(Ij) ∀j} ⊂ H1

Galerkin method: Find uh(t) ∈ Vh such that for all φh ∈ Vh 1 ∂uh ∂t φhdx + µ 1 ∂uh ∂x ∂φh ∂x dx = 0

8 / 73

Linear convection equation

Linear, scalar convection/advection equation (Initial value problem) ut + aux =0 x ∈ R, t > 0 u(x, 0) =f(x) x ∈ R (1) Exact solution u(x, t) = f(x − at) Initial condition is convected with speed a without change of form.

x f(x) u(x, t) at

Hence the extrema of the solution do not change with time. Also the L2-norm of the solution does not change with time. If E(t) is the solution operator u(x, t) = E(t)f(x) = ⇒ E(t)u = u in both sup-norm and L2-norm.

9 / 73

Forward time, backward space (FTBS) un+1

j

− un

j

∆t + aun

j − un j−1

h = 0 ⇒ un+1

j

= (1 − aλ)un

j + aλun j−1

Forward time, central space (FTCS) un+1

j

− un

j

∆t + aun

j+1 − un j−1

2h = 0 ⇒ un+1

j

= un

j + aλ

2 (un

j−1 − un j+1)

Forward time, forward space (FTFS) un+1

j

− un

j

∆t + aun

j+1 − un j

h = 0 ⇒ un+1

j

= (1 + aλ)un

j − aλun j+1

λ = ∆t h

10 / 73

x u x u x u

Backward Forward Central Upwind scheme: switch between backward and forward difference un+1

j

− un

j

∆t + a+ un

j − un j−1

h + a− un

j+1 − un j

h = 0

11 / 73

Stable in maximum norm for any a provided |a|λ ≤ 1 is satisfied.

Courant-Friedrichs-Levy (CFL) number, CFL condition

CFL = |a|λ = |a|∆t h , CFL ≤ 1, ∆t = O(h)

12 / 73

Hyperbolic conservation law

∂u ∂t + ∂ ∂xf(u) = 0 e.g., Burger’s equation ∂u ∂t + ∂ ∂x u2 2 = 0 d dtu(x(t), t) = 0 along d dtx(t) = u(x(t), t)

0.5 1.0 1.5 2.0 2.5 0.5 1.0 1.5 2.0

Solution at times t0, ttc and ttc

1.5 2.0 2.5 0.5 1.0 1.5 2.0

13 / 73

Hyperbolic conservation law

Definition: Weak solution

A function u : R × R+ → R is a weak solution of the IVP ut + f(u)x = 0, (x, t) ∈ R × R+, u(x, 0) = u0(x) together with locally integrable initial data u0 if u is locally integrable and satisfies

∞

• ∞
• −∞

(uφt + f(u)φx)dxdt +

∞

• −∞

u0(x)φ(x, 0)dx = 0, ∀φ ∈ C1

0(R × R+)

Rankine-Hugoniot condition: f(u+) − f(u−) = s(u+ − u−), s = speed of discontinuity

14 / 73

Entropy function

Consider a convex scalar conservation law ut + f(u)x = 0 Assume that there exists a convex function η(u) and another function θ(u) such that η′(u)f ′(u) = θ′(u) Such a pair (η, θ) is called an entropy-entropy flux pair. For Burgers equation, we can choose η(u) = u2, θ(u) = 2 3u3 For smooth solutions ut + f ′(u)ux = 0, η′(u)ut + η′(u)f ′(u)

• θ′(u)

ux = 0, leads to another conservation law ηt + θx = 0

15 / 73

Entropy function

In reality, the conservation law includes some dissipation ut + fx = ǫuxx = ⇒ η′(u)ut + η′(u)f ′(u)ux = ǫη′(u)uxx leads to the entropy equation ηt + θx = ǫ(η(u)ux)x − ǫη′′(u)u2

x ≤ ǫ(η(u)ux)x

since η′′(u) > 0 In the limit of ǫ → 0, we get ηt + θx ≤ 0 This condition must be satisfied in weak sense for all φ ∈ C1

0(R × R+), φ ≥ 0

∞

• R

(η(u)φt + θ(u)φx)dxdt +

• R

η(u0(x))φ(x, 0)dt ≥ 0

16 / 73

Kruzkov’s result

The scalar Cauchy problem ut + f(u)x = 0, f ∈ C1(R) with initial condition u(0, x) = u0(x), u0 ∈ L∞(R) has a unique entropy solution u ∈ L∞(R+ × R) which fulfills (important for numerics)

1 Stability: ||u(t, ·)||L∞ ≤ ||u0||L∞, a.e. in t ∈ R+ 2 Monotone: if u0 ≥ v0 a.e. in R, then

u(t, ·) ≥ v(t, ·) a.e. in R, a.e. in t ∈ R+

17 / 73

Kruzkov’s result

3 TV-diminishing: if u0 ∈ BV (R) then

u(t, ·) ∈ BV (R) and TV (u(t, ·)) ≤ TV (u0)

4 Conservation: if u0 ∈ L1(R) then

• R

u(t, x)dx =

• R

u0(x)dx, a.e. in t ∈ R+

5 Finite domain of dependence: if u, v are two entropy solutions corresponding

to u0, v0 ∈ L∞ and M = max

φ {|f ′(φ)| : |φ| ≤ max(||u0||L∞, ||v0||L∞)}

then

• |x|≤R

|u(t, x) − v(t, x)|dx ≤

• |x|≤R+Mt

|u0(x) − v0(x)|dx

18 / 73

Finite volume method

Divide space domain into finite volumes Ω =

• j

(xj− 1

2 , xj+ 1 2 ),

hj = xj+ 1

2 − xj− 1 2 ,

xj = 1 2(xj− 1

2 + xj+ 1 2 )

Integrate conservation law over finite volume (xj− 1

2 , xj+ 1 2 ) and time slab

(tn, tn+1) tn+1

tn

xj+ 1

2

xj− 1

2

∂u ∂t + ∂f ∂x

• dxdt = 0

Cell average value uj(t) = 1 hj xj+ 1

2

xj− 1

2

u(x, t)dx gives conservation law (exact) (un+1

j

− un

j )hj +

tn+1

tn

[f(xj+ 1

2 , t) − f(xj− 1 2 , t)]dt = 0 19 / 73

Finite volume method

Approximate time integral of flux using solution at tn (explicit scheme) tn+1

tn

f(xj+ 1

2 , t)dt ≈ f(xj+ 1 2 , tn)∆t

leads to finite volume method vn+1

j

− vn

j

∆t + f n

j+ 1

2 − f n

j− 1

2

hj = 0 Cell average values are the unknowns in the finite volume method. vn

j ≈ uj(tn) = 1

hj xj+ 1

2

xj− 1

2

u(x, tn)dx

20 / 73

Godunov’s idea: Riemann problem

Riemann problem defined at each cell face xj+ 1

2

w(x, tn) =

• un

j

x < xj+ 1

2

un

j+1

x > xj+ 1

2

Find exact solution w(x, t) = wR((x − xj+ 1

2 )/(t − tn); un

j , un j+1)

Compute flux f n

j+ 1

2 = f(wR(0; un

j , un j+1))

21 / 73

Linear convection eqn: upwind scheme

ut + aux = 0 Upwind scheme un+1

j

− un

j

∆t + a+ un

j − un j−1

h + a− un

j+1 − un j

h = 0 can be written as finite volume scheme un+1

j

− un

j

∆t + f n

j+ 1

2 − f n

j− 1

2

h = 0 where fj+ 1

2 =

• auj

a > 0 auj+1 a < 0 = 1 2(auj + auj+1) − 1 2|a|(uj+1 − uj)

22 / 73

Kinetic description of gases

• Gas composed of many molecules
• Describe state of gas by a velocity distribution function

f(x, v, t), x ∈ R3, v ∈ R3

• Number density of gas molecules

n(x, t) =

• R3 f(x, v, t)dv =: f
• Mass density

ρ(x, t) = mn(x, t), m = mass of one molecule We will assume that f is rescaled to give mass density, f = ρ.

23 / 73

Kinetic description of gases

• Basic conserved quantities

ρ = f ρu = vf ρe =

• (v2/2)f
• U = ψf ,

ψ =   1 v

1 2|v|2

 

• Evolution of f governed by Boltzmann equation

∂f ∂t + v · ∇f = J(f, f)

• Collisions do not change mass, momentum, energy

ψJ(f, f) = 0 ψ are the only collisional invariants.

24 / 73

Kinetic description of gases

• Equilibrium distribution: Maxwell-Boltzmann distribution function

J(f, f) = 0 = ⇒ f = g = exp(a + b · v + c|v|2)

• r in the more familiar form

g(v) = ρ β π 3/2 exp

• −β|v − u|2

, β = 1 2RT

• Moments of g lead to Euler equations

U = ψg , Fi = viψf ∂g ∂t + vi ∂g ∂xi = 0

• =

⇒ ∂U ∂t + ∂Fi ∂xi = 0 Construct numerical scheme for Boltzmann eqn (linear convection eqn), average over all velocities to obtain scheme for Euler equation

25 / 73

Kinetic model for convection equation

The scalar convection equation ut + cux = 0 (2) can be obtained from a “Maxwell-Boltzmann” distribution g(v, u) = u

• β

π exp

• −β(v − c)2

(3) since g = u and vg = cu The ”Boltzmann equation” for the distribution function is obtained by differentiating equation (3) gt + vgx = (v − c)ux

• β

π exp

• −β(v − c)2

=: Q (4) Note that the moment of the collision term Q is zero, Q = 0 so that taking moments of equation (4) gives us the convection equation (2).

26 / 73

Kinetic numerical flux

Discontinuous solution in a finite volume

• r finite element method

u±(x) = lim

ǫց0 u(x ∓ ǫ) j − 1/2 j + 1/2 j + 3/2 Ij Ij+1 u+

j+1/2

u−

j+1/2

In the kinetic model, molecules carry information due to their streaming motion. This allows us to use the upwind principle to approximate the velocity distribution function at xj+ 1

2 based on the sign of the molecular velocity as

gj+ 1

2 =

• g(v, u+

j+ 1

2 )

v > 0 g(v, u−

j+ 1

2 )

v < 0

27 / 73

Kinetic numerical flux

Numerical flux function Fj+ 1

2 = vgj+ 1 2 = F +(u+

j+ 1

2 ) + F −(u−

j+ 1

2 ) =: F(u+

j+ 1

2 , u−

j+ 1

2 )

The split fluxes are obtained by integrating over half velocity spaces, and are given by F ±(u) =

• v±g
• = cuA± + uB±

where A± = 1 2[1 ± erf(s)], B± = ± 1 2√πβ exp(−s2), s = c

• β
• F + is due to all molecules moving to the right (v > 0) while the flux F − is

due to all molecules moving to the left (v < 0)

• F(u−, u+) is consistent, i.e., F(u, u) = cu and is obviously a smooth

function of u.

• The split flux F +(u) is an increasing function of u while F −(u) is a

decreasing function. Hence the numerical flux F(u−, u+) is a monotone flux.

28 / 73

Kinetic numerical flux

• Central flux with dissipation

F(u+, u−) = 1 2(cu+ + cu−) + 1 2D(u+ − u−) (5) where D = c erf(s) + 1 √πβ exp(−s2) > 0

• In the limit of β → ∞, D → |c|

F(u+, u−) =

• cu+,

c > 0 cu−, c < 0

29 / 73

DG scheme for convection equation

Ω =

N

• j=1

Ij, Ij = [xj− 1

2 , xj+ 1 2 ],

j = 1, . . . , N The DG method uses the broken space of polynomials of degree k V k

h =

• φ ∈ L2(Ω) : φ|Ij ∈ Pk(Ij)
• Multiply equation (2) by a test function φh ∈ V k

h and integrate over element Ij

• Ij

φh∂tuhdx −

• Ij

cuh∂xφhdx + F(uh)j+ 1

2 φh(x+

j+ 1

2 ) − F(uh)j− 1 2 φh(x−

j− 1

2 ) = 0

Use numerical flux at the element interface F(uh)j+ 1

2 = F(u+

j+ 1

2 , u−

j+ 1

2 ) 30 / 73

DG scheme for convection equation

If we take φh = uh and sum over all elements, we obtain the energy equation 1 2 d dtuh2 −

N

• j=1
• Ij

cuh∂xuhdx +

N

• j=1

F(uh)j+ 1

2 uhj+ 1 2 = 0

(6) where we have defined the inter-element jump in the solution by uhj+ 1

2 := uh(x+

j+ 1

2 ) − uh(x−

j+ 1

2 )

After some simple manipulations, the energy equation can be written as 1 2 d dtuh2 + 1 2D

N

• j=1

uh2

j+ 1

2 = 0

and hence the DG scheme is stable in energy norm for any degree of basis functions.

31 / 73

Kinetic model for convection diffusion equation

Consider the linear convection-diffusion equation ut + cux = µuxx, µ > 0 (7) Kinetic formulation using the Boltzmann equation with BGK model ft + vfx = g − f τR (8) τR is relaxation time scale towards local equilibrium distribution g. Chapman-Enskog expansion f = f0 + τRf1 + τ 2

Rf2 + . . .

(9) Substituting this in equation (8) and collecting terms with the same coefficient of τR, we get f0 = g, f1 = −(gt + vgx) The first two terms constitute the Chapman-Enskog distribution function f ce = f0 + τRf1 = [u − τR(v − c)ux]

• β

π exp

• −β(v − c)2

(10)

32 / 73

Kinetic model for convection diffusion equation

Upon taking moments of this distribution function, we obtain f ce = u, vf ce = cu − τR 2β ux In order to recover the correct flux of equation (7), we have to choose τR = 2βµ = ⇒ vf ce = cu − µux Since f ce is a truncated solution in the Chapman-Enskog expansion it does not satisfy equation (8). We can derive an equation for f ce by differentiating equation (10) f ce

t + vf ce x =

• (v − c)ux + (µ − τR(v − c)2)uxx − τRµ(v − c)uxxx

β π exp

• −β(v − c)2

=: Q Upon taking moments of this equation, the right hand side vanishes, i.e., Q = 0 and we obtain the convection-diffusion equation (7).

33 / 73

Flux vector splitting

F = cux − µux Split fluxes due to positive and negative velocity particles F ± =

• v±f ce

= (cu − µux)A± + uB± Numerical flux function F(u+, u−) = F +(u+) + F −(u−) = Fc(u+, u−) + Fd(u+

x , u− x )

where the convective flux Fc is identical to equation (5), while the dissipative flux is given by Fd(u+

x , u− x ) = −µu+ x A+ − µu− x A−

Note that the diffusive flux can also be written in split form as Fd(u+

x , u− x ) = F + d (u+ x ) + F − d (u− x ),

F ±

d (ux) = −µuxA±

34 / 73

DG scheme for convection-diffusion equation

The DG discretization of equation (7) is obtained by multiplying it by a test function φh and integrating over any element Ij

• Ij

φh∂tuhdx −

• Ij

cuh∂xφhdx + µ

• Ij

(∂xuh)(∂xφh)dx + F(uh)j+ 1

2 φh(x+

j+ 1

2 ) − F(uh)j− 1 2 φh(x−

j− 1

2 ) = 0

(11) In the above equation, the inter-element flux is approximated by the numerical flux which is defined as F(uh) = Fc(u+

h , u− h ) + Fd(∂xu+ h , ∂xu− h )

and it includes both convective and diffusive fluxes. Substituting φh = uh and summing over all elements we obtain the following energy equation 1 2 d dtuh2+1 2D

N

• j=1

uh2

j+ 1

2 +µ

N

• j=1
• Ij

(∂xuh)2dx+

N

• j=1

Fd(∂xuh)j+ 1

2 uhj+ 1 2 = 0

(12)

35 / 73

Test Case

We will consider the convection-diffusion equation (7) in the domain (−1, +1) together with the initial condition u(x, 0) = − sin(πx) (13) and periodic boundary conditions. The exact solution to this problem is given by u(x, t) = − exp(−µπ2t) sin(π(x − ct)) (14) Test Case 1 2 µ 0.001 1 tf 30 0.5

Table : Parameters for the test case based on convection-diffusion equation (7)

36 / 73

Naive DG scheme for convection-diffusion equation

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 Exact KFVS

• 1
• 0.5

0.5 1

• 0.5

0.5 Exact KFVS

(a) (b)

Figure : Convection-diffusion problem using equation (11) and P1 basis functions. (a) Test case 1, (b) Test case 2

37 / 73

DG scheme with stabilization

Motivated by the energy analysis, we propose the following stabilized DG scheme

• Ij

φh∂tuhdx −

• Ij

cuh∂xφhdx + µ

• Ij

∂xuh∂xφhdx + F(uh)j+1/2φh(x+

j+ 1

2 ) − F(uh)j− 1 2 φh(x−

j− 1

2 )

− F +

d (∂xφ+ h )j+ 1

2 uhj+ 1 2 − F −

d (∂xφ− h )j− 1

2 uhj− 1 2 = 0

(15) Choosing φh = uh and summing up over all elements gives us 1 2 d dtuh2 + 1 2D

N

• j=1

uh2

j+ 1

2 + µ

N

• j=1
• Ij

(∂xuh)2dx = 0 (16) which shows that the scheme is stable.

38 / 73

Stabilized DG scheme for convection-diffusion equation

• 1
• 0.5

0.5 1

• 0.4
• 0.2

0.2 0.4 Exact KFVS

• 1
• 0.5

0.5 1

• 0.4
• 0.2

0.2 0.4 Exact KFVS

(a) (b)

Figure : Results for Test case 2 using equation (15): (a) P1 and (b) P2 basis functions

39 / 73

Cell energy balance

Multiply the convection-diffusion equation by u and integrate over element Ij d dtu2

Ij +

1 2cu2 j+ 1

2

j− 1

2

+ [−µuux]

j+ 1

2

j− 1

2 = −µ

• Ij

(ux)2dx ≤ 0 (17) Set φh = uh in the stabilized DG scheme to obtain d dtu2

Ij+

¯ Fc j+ 1

2

j− 1

2 +

¯ Fd j+ 1

2

j− 1

2 = −µ

• IJ

(∂xuh)2dx−1 4D uh2

j− 1

2 −1

4D uh2

j+ 1

2 ≤ 0

where ¯ Fc(uh)j+ 1

2 = Fc(uh)j+ 1 2 {

{uh} }j+ 1

2 − c

2

• u2

h

• j+ 1

2

and ¯ Fd(uh)j+ 1

2 = F +

d (∂xu+ h )j+ 1

2 uh(x−

j+ 1

2 ) + F −

d (∂xu− h )j+ 1

2 uh(x+

j+ 1

2 )

which are consistent numerical fluxes for 1

2cu2 and −µuux respectively.

40 / 73

Global variational form

Find uh(t) ∈ V k

h such that for all φh ∈ V k h , the following equation is satisfied

• Ω

φh∂tuhdx −

• Ω

cuh∂xφhdx + µ

• Ω

(∂xuh)(∂xφh)dx +

• j

Fc(uh)j+ 1

2 φhj+ 1 2

+

• j

Fd(∂xuh)j+ 1

2 φhj+ 1 2 + ǫ

• j

Fd(∂xφh)j+ 1

2 uhj+ 1 2

+

• j

δh(uh)j+ 1

2 φhj+ 1 2 = 0

(18) where the last term is an additional term that has been added, called the interior penalty term, and is of the form δh(u) = Cip µ h u where Cip is a non-dimensional positive number.

41 / 73

Global variational form

Define bilinear form due to viscosity terms aǫ(uh, φh) =µ

• Ω

(∂xuh)(∂xφh)dx +

• j

Fd(∂xuh)j+ 1

2 φhj+ 1 2

+ǫ

• j

Fd(∂xφh)j+ 1

2 uhj+ 1 2 +

• j

δh(uh)j+ 1

2 φhj+ 1 2

1 If ǫ = −1 we obtain the scheme of equation (15) without the interior penalty

• term. The scheme is non-symmetric since

a−1(uh, φh) = a−1(φh, uh) Hence this is refered to as the NIPG scheme. As we have seen, this scheme satisfies a cell energy inequality.

2 If ǫ = +1, the discretization of the diffusive terms becomes symmetric

a+1(uh, φh) = a+1(φh, uh) and the scheme is known as SIPG scheme. However, the SIPG scheme does not satisfy a cell energy inequality. This scheme is stable only with the addition of a penalty term δh, with a sufficiently large value of Cip.

42 / 73

Test Case 1 with IP scheme: 20 cells, Cip = 10

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 Exact KFVS

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 Exact KFVS

NIPG, P1 NIPG, P2

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 Exact KFVS

• 1
• 0.5

0.5 1

• 1
• 0.5

0.5 1 Exact KFVS

SIPG, P1 SIPG, P2

43 / 73

Test Case 2 with IP scheme: 20 cells, Cip = 10

• 1
• 0.5

0.5 1

• 0.4
• 0.2

0.2 0.4 Exact KFVS

• 1
• 0.5

0.5 1

• 0.4
• 0.2

0.2 0.4 Exact KFVS

NIPG, P1 NIPG, P2

• 1
• 0.5

0.5 1

• 0.4
• 0.2

0.2 0.4 Exact KFVS

• 1
• 0.5

0.5 1

• 0.4
• 0.2

0.2 0.4 Exact KFVS

SIPG, P1 SIPG, P2

44 / 73

cells dofs L2 H1 20 40 1.259e-02

• 2.757e-01
• 40

80 3.535e-03 1.83 1.356e-01 1.02 80 160 1.245e-03 1.51 6.639e-02 1.03 160 320 3.511e-04 1.83 3.257e-02 1.03 320 640 3.986e-05 3.14 1.614e-02 1.01

Table : Test case 1: NIPG, P1 basis

cells dofs L2 H1 20 60 3.376e-04

• 1.416e-02
• 40

120 7.361e-05 2.20 3.612e-03 1.97 80 240 1.660e-05 2.15 9.186e-04 1.98 160 480 3.433e-06 2.27 2.243e-04 2.03 320 960 6.892e-07 2.32 5.103e-05 2.14

Table : Test case 1: NIPG, P2 basis

45 / 73

cells dofs L2 H1 20 40 1.245e-02

• 2.866e-01
• 40

80 1.843e-03 2.76 1.403e-01 1.03 80 160 3.170e-04 2.54 6.807e-02 1.04 160 320 6.181e-05 2.36 3.293e-02 1.05 320 640 1.351e-05 2.19 1.615e-02 1.03

Table : Test case 1: SIPG, P1 basis

cells dofs L2 H1 20 60 1.354e-04

• 1.367e-02
• 40

120 1.622e-05 3.06 3.296e-03 2.05 80 240 1.900e-06 3.09 7.692e-04 2.10 160 480 2.153e-07 3.14 1.748e-04 2.14 320 960 2.487e-08 3.11 4.324e-05 2.02

Table : Test case 1: SIPG, P2 basis

46 / 73

cells dofs L2 H1 20 40 1.592e-02

• 1.088e-01
• 40

80 8.018e-03 0.99 5.451e-02 1.00 80 160 4.031e-03 0.99 2.727e-02 1.00 160 320 2.022e-03 1.00 1.364e-02 1.00 320 640 1.013e-03 1.00 6.820e-03 1.00

Table : Test case 2: NIPG, P1 basis

cells dofs L2 H1 20 60 7.467e-04

• 4.778e-03
• 40

120 1.860e-04 2.01 1.195e-03 2.00 80 240 4.648e-05 2.00 2.989e-04 2.00 160 480 1.162e-05 2.00 7.472e-05 2.00 320 960 2.906e-06 2.00 1.868e-05 2.00

Table : Test case 2: NIPG, P2 basis

47 / 73

cells dofs L2 H1 20 40 2.255e-03

• 1.124e-01
• 40

80 5.653e-04 2.00 5.638e-02 1.00 80 160 1.415e-04 2.00 2.823e-02 1.00 160 320 3.538e-05 2.00 1.412e-02 1.00 320 640 8.846e-06 2.00 7.063e-03 1.00

Table : Test case 2: SIPG, P1 basis

cells dofs L2 H1 20 60 1.355e-04

• 9.936e-03
• 40

120 1.713e-05 2.98 2.497e-03 1.99 80 240 2.158e-06 2.99 6.268e-04 1.99 160 480 2.709e-07 2.99 1.571e-04 2.00 320 960 3.394e-08 3.00 3.932e-05 2.00

Table : Test case 2: SIPG, P2 basis

48 / 73

Comparison with BGK finite volume

Torrillhon and Xu u(x, 0) = 4+ 8 π sin(πx/2)+ 16 3π sin(3πx/2) BGK finite volume scheme for scalar convection-diffusion equation exhibits

• nly first order accuracy in the

asymptotic limit, though on coarse grids it shows third order accuracy. They find a non-monotone convergence of the error with grid size.

49 / 73

DG KFVS

β = 1, µ = 0.005, c = 1 The grid is refined starting with 50 cells and increased by 50 cells in each refinement step. The initial condition is taken to be same as in Torrillhon-Xu.

10

2

10

3

10

−4

10

−3

10

−2

log(N) ||u−uh||

Slope = 1 Slope = 2

NIPG SIPG 10

2

10

3

10

−7

10

−6

10

−5

10

−4

log(N) ||u−uh||

Slope = 2 Slope = 3

NIPG SIPG

P1 P2

50 / 73

Euler equations

∂tU + ∂iFi = 0 (19) Hyperbolicity: for any vector n = (n1, n2, n3), the matrices A(U, n) = Ai(U)ni, Ai(U) = F ′

i(U)

have all real eigenvalues and a complete set of eigenvectors. The Euler equations admit an entropy function η(U) with entropy fluxes θi(U), i = 1, 2, 3 such that θ′

i(U) = η′(U)F ′ i(U)

This implies that for smooth solutions, an additional conservation law is satisfied ∂tη + ∂iθi = 0 while for discontinuous solutions ∂tη + ∂iθi ≤ 0

51 / 73

Euler equations

η(U) is a strictly convex function. Define the entropy variables by V ⊤ = η′(U) (20) and the mapping U → V is one-to-one. Define the conjugate functions η∗(V ) = V · U(V ) − η(U(V )), θ∗

i (V ) = V · Fi(U(V )) − θi(U(V ))

Differentiating these expressions gives η∗′(V ) = U ⊤ θ∗′

i (V ) = Fi(U(V ))⊤

It then follows that the matrices U ′(V ) = η∗′′(V ) and F ′

i(U(V ))U ′(V ) = θ∗′′ i (V )

are symmetric. Moreover the matrix U ′(V ) = η′′(U(V ))−1

52 / 73

Euler equations

is positive-definite. Making a change of variables from U to V in equation (19), we get ˜ A0∂tV + ˜ Ai∂iV = 0 where we have defined ˜ A0 = U ′(V ) = η∗′′(V ), ˜ Ai = Ai ˜ A0 = θ∗′′

i (V )

the matrix ˜ A0 is symmetric positive definite, while the matrices ˜ Ai are symmetric. Existence of convex entropy ⇐ ⇒ symmetrizability

53 / 73

DG scheme for Euler equations

54 / 73

DG scheme for Euler equations

n K K′ V +

h

V −

h

Triangulation Th = {K} of the domain Ω by polygonal, non-overlapping

• elements. We assume that the edges of

the triangulation are oriented, so that each edge has a unit normal vector n with trace values of the discontinuous functions on the edge. The DG scheme for an interior element K is given by multiplying the Euler equation by a test function Wh and doing some integration by parts to obtain

• K

Wh·∂tU(Vh)dx−

• K

Fi(Vh)·∂iWhdx+

• ∂K

Hc(V +

h , V − h , n)·W − h dσ = 0 (21)

Numerical flux Hc(V +

h , V − h , n) =

• (v · n)+ψg(V +

h )

• +
• (v · n)−ψg(V −

h )

• 55 / 73

DG scheme for Euler equations

56 / 73

Entropy stability of DG scheme

∂tη + ∂iθi ≤ 0 = ⇒

• K

∂tηdx +

• ∂K

θinidσ ≤ 0 substitute Wh = Vh in equation (21) to get

• K

∂tη(Vh)dx = −

• ∂K
• Θ(V +

h , V − h , n) + D(V + h , V − h , n)

• dσ

where we have defined Θ(V +, V −, n) = { {V } }+

− · Hc(V +, V −, n) − {

{θ∗ · n} }+

−

which is a consistent and conservative numerical flux for the entropy flux θ · n, i.e., Θ(V, V, n) = V · (Fini) − θ∗

i ni = θini

and D(V +, V −, n) = 1 2

• V +

− · Hc(V +, V −, n) − θ∗ · n+ −

• ≥ 0

? Yes (Barth)

57 / 73

Chapman-Enskog distribution function

Consider the Boltzmann equation with the BGK model for the collision term ∂tf + vi∂if = g − f τR (22) If we write the solution as a series in the relaxation time τR as in equation (9), we

• btain the following first two terms of the solution

f0 = g, f1 = −(∂tg + vi∂ig) The distribution function g depends on the macroscopic variables ρ, u, T; we substitute the inviscid equations governing these quantities into the expression for f1 to obtain the following Chapman-Enskog velocity distribution function f ce = f0 + τRf1 = g

• 1 − ρ

p2 τijcicj − ρqici p2

• 1 − 2

5β|c|2

• (23)

τR = µ/p and ci = vi − ui

58 / 73

Chapman-Enskog distribution function

U = ψf ce and Fi + Gi = viψf ce , i.e., the moments of the Chapman-Enskog distribution yield the conserved variables and fluxes for the Navier-Stokes equations. Since f ce is a truncated solution in the Chapman-Enskog expansion, it does not satisfy the Boltzmann equation (22); we can derive an equation for f ce as ∂tf ce + vi∂if ce = Q and it can be shown through explicit but lengthy computations that ψQ = 0

59 / 73

DG scheme for NS equations

In terms of the entropy variables, the Navier-Stokes equations can be written as ˜ A0∂tV + ˜ Ai∂iV − ∂i(Kij∂jV ) = 0, Gi = −Kij∂jV and the matrix K = [Kij] is symmetric positive semi-definite (Hughes, Dutt). Taking dot product with entropy variables, we obtain ∂tη + ∂iθi + ∂i(V · Gi) = −(∂iV )⊤Kij∂jV ≤ 0 (24) which is essentially the entropy condition. ΓI = set of all interior edges/faces of the finite element grid and let Γ = set of all boundary edges/faces

60 / 73

DG scheme for NS equations

The time continuous DG scheme for Navier-Stokes equations is given by

• Ω

Wh · ∂tU(Vh)dx −

• Ω

[Fi(Vh) + Gi(Vh, ∇Vh)] · ∂iWhdx +

• ΓI

H(V +

h , ∇V + h , V − h , ∇V − h , n) · Wh dσ

+ ǫ

• ΓI

Hd(V +

h , ∇W + h , V − h , ∇W − h , n) · Vh dσ

+

• ΓI

δh(Vh) · Wh dσ + NΓ(Vh, Wh) = 0 Choosing ǫ = −1 yields the NIPG scheme while ǫ = +1 yields the SIPG scheme. The numerical flux function is given by kinetic splitting as H(V +, ∇V +, V −, ∇V −, n) = Hc(V +, V −, n) + Hd(V +, ∇V +, V −, ∇V −, n)

61 / 73

DG scheme for NS equations

Here Hc is the KFVS numerical flux for Euler equationsas given by equation (54) while Hd is the KFVS numerical flux for the viscous fluxes obtained from Chapman-Enskog distribution. This flux can be written as Hd(V +, ∇V +, V −, ∇V −, n) = H+

d (V +, ∇V +, n) + H− d (V −, ∇V −, n)

The interior penalty term is of a similar form as in the scalar problem and given by δh(Vh) = Cip k2ν h U(Vh) The boundary conditions are imposed through the boundary terms NΓ in a weak form and given by NΓ(Vh, Wh) =

• Γ

H(ΓV +

h , Γ∇V + h , ΓV − h , Γ∇V − h , n) · W + h dσ

+

• Γ

Hd(ΓV +

h , Γ∇W + h , ΓV − h , Γ∇W − h , n) · [V + h − ΓV + h ]dσ

+

• Γ

δΓ(Vh) · W +

h dσ

62 / 73

DG scheme for NS equations

The operator Γ modifies the boundary trace of the solution to apply the boundary conditions; for example, in the case of a no-slip boundary ΓV +

h = ΓV − h =

      V +

1

V +

5

      , Γ∇V +

h = Γ∇V − h = ∇V + h

For an isothermal wall, V +

5

is modified by using the specified wall temperature. For an adiabatic wall, ∇V +

h is modified so that the heat flux is zero.

63 / 73

Entropy stability of NIPG scheme

• K

∂tηdx +

• ∂K

[θini + V · (Gini)]dσ = −

• K

(∂iV )⊤Kij∂jV dx ≤ 0 For an interior element K and substitute Wh = Vh to obtain

• K

∂tη(Vh)dx +

• ∂K
• Θ(V +

h , V − h , n) + D(V + h , V − h , n)

• dσ

+

• ∂K
• H+

d (V + h , ∇V + h , n) · V − h + H− d (V − h , ∇V − h , n) · V + h

• dσ

−

• K

Gi(Vh, ∇Vh) · ∂iV dx = 0 The first two integrals are common with the Euler equations, while the last two terms are the contributions from NS equations. We observe that H+

d (V + h , ∇V + h , n) · V − h + H− d (V − h , ∇V − h , n) · V + h

is a consistent and conservative numerical flux for the viscous entropy flux V · (Gini), while from the Navier-Stokes equations, we have Gi(Vh, ∇Vh) · ∂iVh = −(∂iVh)⊤Kij∂jVh ≤ 0 which leads to a cell entropy inequality for the DG scheme.

64 / 73

Convergence rate for NS equations

We assume the following form for the exact solution ρ(x) = 1+ 1 2 cos(2πx), u(x) = 10x2(1−x)2 sin(2πx), T(x) = 1+2x2(1−x)2 which solves the following stationary 1-D NS equations with source terms (ρu)x = f1, (p + ρu2)x − τx = f2, ((ρe + p)u)x − (τu)x + qx = f3 together with the boundary conditions u(0) = u(1) = 0, q(0) = q(1) = 0 The source terms f1, f2, f3 are computed by substituting the exact solutions into the NS equations.

65 / 73

Convergence rate for NS equations

N u error Rate ux error Rate T error Rate Tx error Rate 40 5.210e-04

• 1.406e-01
• 8.919e-04
• 1.550e-02
• 80

1.752e-04 1.57 7.029e-02 1.00 4.427e-04 1.01 8.036e-03 0.94 160 7.539e-05 1.21 3.514e-02 1.00 2.416e-04 0.87 4.165e-03 0.94 320 3.657e-05 1.04 1.756e-02 1.00 1.285e-04 0.90 2.130e-03 0.96 Table : Convergence rate of NIPG scheme for 1-D NS equation using P1 basis functions N u error Rate ux error Rate T error Rate Tx error Rate 40 2.027e-05

• 6.637e-03
• 2.237e-05
• 7.755e-04
• 80

2.431e-06 3.06 1.500e-03 2.14 5.263e-06 2.08 2.181e-04 1.83 160 3.569e-07 2.76 3.504e-04 2.09 1.403e-06 1.90 5.842e-05 1.90 320 6.960e-08 2.35 8.537e-05 2.03 3.739e-07 1.90 1.506e-05 1.95 Table : Convergence rate of NIPG scheme for 1-D NS equation using P2 basis functions

66 / 73

Convergence rate for NS equations

N u error Rate ux error Rate T error Rate Tx error Rate 40 3.128e-07

• 1.630e-04
• 4.935e-07
• 8.941e-06
• 80

2.400e-08 3.70 2.000e-05 3.02 7.072e-08 2.80 1.059e-06 3.07 160 2.147e-09 3.48 2.465e-06 3.02 8.772e-09 3.01 1.255e-07 3.07 320 1.885e-10 3.50 3.033e-07 3.02 9.027e-10 3.28 1.419e-08 3.14 Table : Convergence rate of NIPG scheme for 1-D NS equation using P3 basis functions N u error Rate ux error Rate T error Rate Tx error Rate 40 4.702e-04

• 1.404e-01
• 4.295e-04
• 1.435e-02
• 80

1.175e-04 2.00 7.025e-02 0.99 1.038e-04 2.04 8.765e-03 0.71 160 2.940e-05 1.99 3.513e-02 0.99 2.576e-05 2.01 4.946e-03 0.82 320 7.355e-06 1.99 1.757e-02 0.99 6.437e-06 2.00 2.634e-03 0.90 Table : Convergence rate of SIPG scheme for 1-D NS equation using P1 basis functions

67 / 73

Convergence rate for NS equations

N u error Rate ux error Rate T error Rate Tx error Rate 40 1.849e-05

• 6.253e-03
• 2.914e-06
• 8.815e-04
• 80

2.108e-06 3.13 1.453e-03 2.10 3.823e-07 2.93 2.340e-04 1.91 160 2.477e-07 3.08 3.473e-04 2.06 4.967e-08 2.94 6.134e-05 1.93 320 3.022e-08 3.03 8.541e-05 2.02 6.317e-09 2.97 1.569e-05 1.96 Table : Convergence rate of SIPG scheme for 1-D NS equation using P2 basis functions N u error Rate ux error Rate T error Rate Tx error Rate 40 2.851e-07

• 1.629e-04
• 4.876e-08
• 6.248e-06
• 80

1.840e-08 3.95 2.002e-05 3.02 2.968e-08 4.03 8.944e-07 2.80 160 1.190e-09 3.95 2.467e-06 3.02 1.570e-10 4.24 1.209e-07 2.88 320 7.722e-11 3.94 3.035e-07 3.02 5.401e-12 4.86 1.561e-08 2.95 Table : Convergence rate of SIPG scheme for 1-D NS equation using P3 basis functions

68 / 73

NS shock structure

Mach number ahead of the shock is M1 = 1.5, γ = 5/3, Pr = 2/3 viscosity law is given by µ = µ1(T/T1)0.8 where the subscript “1” denotes pre-shock conditions and µ1 = 0.0005.

69 / 73

N=100, P1 N=100, P2

0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Density x NIPG SIPG Exact 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Density x NIPG SIPG Exact

N=200, P1 N=200, P2

0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Density x NIPG SIPG Exact 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Density x NIPG SIPG Exact

Figure : Density for shock structure problem using entropy variables

70 / 73

N=100, P1 N=100, P2

• 0.045
• 0.04
• 0.035
• 0.03
• 0.025
• 0.02
• 0.015
• 0.01
• 0.005

0.005

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Shear stress x NIPG SIPG Exact

• 0.045
• 0.04
• 0.035
• 0.03
• 0.025
• 0.02
• 0.015
• 0.01
• 0.005

0.005

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Shear stress x NIPG SIPG Exact

N=200, P1 N=200, P2

• 0.045
• 0.04
• 0.035
• 0.03
• 0.025
• 0.02
• 0.015
• 0.01
• 0.005

0.005

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Shear stress x NIPG SIPG Exact

• 0.045
• 0.04
• 0.035
• 0.03
• 0.025
• 0.02
• 0.015
• 0.01
• 0.005

0.005

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Shear stress x NIPG SIPG Exact

Figure : Shear stress for shock structure problem using entropy variables

71 / 73

N=100, P1 N=100, P2

• 0.04
• 0.035
• 0.03
• 0.025
• 0.02
• 0.015
• 0.01
• 0.005

0.005

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Heat flux x NIPG SIPG Exact

• 0.04
• 0.035
• 0.03
• 0.025
• 0.02
• 0.015
• 0.01
• 0.005

0.005

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Heat flux x NIPG SIPG Exact

N=200, P1 N=200, P2

• 0.04
• 0.035
• 0.03
• 0.025
• 0.02
• 0.015
• 0.01
• 0.005

0.005

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Heat flux x NIPG SIPG Exact

• 0.04
• 0.035
• 0.03
• 0.025
• 0.02
• 0.015
• 0.01
• 0.005

0.005

• 0.14
• 0.135
• 0.13
• 0.125
• 0.12
• 0.115
• 0.11

Heat flux x NIPG SIPG Exact

Figure : Heat flux for shock structure problem using using entropy variables

72 / 73

Summary

• KFVS for diffusion - unstable scheme
• Stabilize using KFVS fluxes - energy and entropy stability
• Non-symmetric scheme - cell-wise stability, but sub-optimal accuracy
• Symmetric scheme necessary for optimal accuracy, globally stable
• Holds for arbitrary degree of polynomials
• Applications: coupling continuum and molecular models

73 / 73