[PPT] - PSLV-C23 launch, 30 June 2014 @ 9:52am IST PSLV-C23 will carry a 714 PowerPoint Presentation

SLIDE 1

PSLV-C23 launch, 30 June 2014 @ 9:52am IST Satish Dhawan Space Centre SHAR, Sriharikota PSLV-C23 will carry a 714 kg French Earth Observation Satellite SPOT-7 as the main

payload. Also, the 14 kg AISAT
f Germany, NLS7.1 (CAN-X4)

& NLS7.2 (CAN-X5) of Canada each weighing 15 kg, and the 7 kg VELOX-1 of Singapore

1 / 58

SLIDE 2

Aryabhatta, India’s first satellite Launched on 19 April 1975, by KOSMOS-3M (Soviet-Union)

2 / 58

SLIDE 3

DG-KFVS schemes for convection-diffusion equations

Praveen. C

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

Dept. of Mathematics

University of W¨ urzburg 30 June, 2014

3 / 58

SLIDE 4

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

4 / 58

SLIDE 5

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

5 / 58

SLIDE 6

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

Pressure

5

Kinetic meshless method

6 / 58

SLIDE 7

Flow through scramjet intake

Density Kinetic meshless method

7 / 58

SLIDE 8

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 = ρ.

8 / 58

SLIDE 9

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.

9 / 58

SLIDE 10

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ψg ∂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

10 / 58

SLIDE 11

Kinetic model for convection equation

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

β

π exp

−β(v − c)2

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

β

π exp

−β(v − c)2

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

11 / 58

SLIDE 12

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

12 / 58

SLIDE 13

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) = ∞ vgdv, F −(u) =

−∞

vgdv, F ±(u) = 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.

13 / 58

SLIDE 14

Kinetic numerical flux

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.

Central flux with dissipation

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

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

F(u+, u−) =

cu+,

c > 0 cu−, c < 0

14 / 58

SLIDE 15

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 (1) 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 ) 15 / 58

SLIDE 16

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

(5) 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.

16 / 58

SLIDE 17

Kinetic model for convection diffusion equation

Consider the linear convection-diffusion equation ut + cux = µuxx

r

ut + (cu − µux)x = 0, µ > 0 (6) Kinetic formulation using the Boltzmann equation with BGK model ft + vfx = g − f τR (7) τR is relaxation time scale towards local equilibrium distribution g. Chapman-Enskog expansion f = f0 + τRf1 + τ 2

Rf2 + . . .

(8) Substituting this in equation (7) 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

(9)

17 / 58

SLIDE 18

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 (6), 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 (7). We can derive an equation for f ce by differentiating equation (9) 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 (6).

18 / 58

SLIDE 19

Flux vector splitting

F = cux − µux Split fluxes due to positive and negative velocity particles F ± = v ± |v| 2 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 (4), while the dissipative flux is given by Fd(u+

x , u− x ) = −µu+ x A+ − µu− x A− =: F + d (u+ x ) + F − d (u− x )

19 / 58

SLIDE 20

DG scheme for convection-diffusion equation

The DG discretization of equation (6) 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

(10) In the above equation, the inter-element flux is approximated by the kinetic 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.

20 / 58

SLIDE 21

Test Case

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

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

21 / 58

SLIDE 22

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 (10) and P1 basis functions. (a) Test case 1, (b) Test case 2

22 / 58

SLIDE 23

DG scheme with stabilization

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

(13) 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

(14) 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 (15) which shows that the scheme is stable.

23 / 58

SLIDE 24

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 (14): (a) P1 and (b) P2 basis functions

24 / 58

SLIDE 25

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 (16) 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.

25 / 58

SLIDE 26

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

(17) 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.

26 / 58

SLIDE 27

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 + Cip

µ h

j

uhj+ 1

2 φhj+ 1 2

1 If ǫ = −1 we obtain the scheme of equation (14) 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.

27 / 58

SLIDE 28

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

28 / 58

SLIDE 29

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

29 / 58

SLIDE 30

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

30 / 58

SLIDE 31

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

31 / 58

SLIDE 32

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

32 / 58

SLIDE 33

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

33 / 58

SLIDE 34

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.

34 / 58

SLIDE 35

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

35 / 58

SLIDE 36

Euler equations

∂tU + ∂iFi = 0 (18) 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

36 / 58

SLIDE 37

Euler equations

η(U) is a strictly convex function. Define the entropy variables by V ⊤ = η′(U) (19) 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 is positive-definite.

37 / 58

SLIDE 38

Euler equations

Making a change of variables from U to V in equation (18), 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

38 / 58

SLIDE 39

DG scheme for Euler equations

39 / 58

SLIDE 40

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 (20)

Numerical flux Hc(V +

h , V − h , n) =

(v · n)+ψg(V +

h )

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

h )

40 / 58

SLIDE 41

DG scheme for Euler equations

41 / 58

SLIDE 42

Entropy stability of DG scheme

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

K

∂tηdx +

∂K

θinidσ ≤ 0 substitute Wh = Vh in equation (20) 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)

42 / 58

SLIDE 43

Chapman-Enskog distribution function

Consider the Boltzmann equation with the BGK model for the collision term ∂tf + vi∂if = g − f τR (21) If we write the solution as a series in the relaxation time τR as in equation (8), 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

(22)

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

43 / 58

SLIDE 44

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 (21); 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

44 / 58

SLIDE 45

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 (23) 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

45 / 58

SLIDE 46

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)

46 / 58

SLIDE 47

DG scheme for NS equations

Here Hc is the KFVS numerical flux for Euler equationsas given by equation (39) 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σ

47 / 58

SLIDE 48

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.

48 / 58

SLIDE 49

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.

49 / 58

SLIDE 50

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.

50 / 58

SLIDE 51

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

51 / 58

SLIDE 52

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

52 / 58

SLIDE 53

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

53 / 58

SLIDE 54

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.

54 / 58

SLIDE 55

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

55 / 58

SLIDE 56

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

56 / 58

SLIDE 57

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

57 / 58

SLIDE 58

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 inequality
Holds for arbitrary degree of polynomials
Applications: coupling continuum and molecular models

58 / 58