Other Examples for φ(θ)

1. φ(θ) = 1. This is the Lax-Wendroff scheme.

2. φ(θ) = θ. This is Beam-Warming.

3. Anyφ between φ_B−W andφ_L−W with0 ≤ φ ≤ 2, 0 ≤ ^φ(θ)_θ ≤ 2 is second order.

4. Van Leer’s minmod

φ(θ) = θ + |θ|

1 + |θ|. It is a smooth limiter withφ(1) = 1.

5. Roe’s superbee

φ(θ) = max(0, min(1, 2θ), min(θ, 2))

7.1. FLUX SPLITTING METHODS 89

There are two kinds of extensions. One is thea < 0 case, and the other is the linear system case.

Fora < 0, we let

In the linear system case, our equation is

u_t+ Au_x= 0. (1.3)

We can decompose A so that A = RΛR⁻¹ with Λ = diag(λ₁, · · · , λn) constituting by A’s

7.2 High Order Godunov Methods

Algorithm

1. Reconstruction: start from cell averages{Ujⁿ}, we reconstruct a piecewise polynomial func-tionu(x, t˜ n).

If 2. is an exact solver, using Z _tn+1

7.2. HIGH ORDER GODUNOV METHODS 91 3’. U_jⁿ⁺¹= U_jⁿ+_∆x^∆t( ˜f_j−1

2 − ˜f_j+1 2)

1. Reconstruction: We want to construct a polynomial in each cell. The main criterions are (1) high order in regions whereu is smooth and u_x 6= 0

(2) total variation no increasing.

In other words, suppose we are given a functionu(x), let U_j = 1

∆x Z x_{j+ 1}

2

x_{j− 1}

2

u(x) dx

From {Uj}, we can use some reconstruct algorithm to construct a function ˜u(x). We want the reconstruction algorithm to satisfy

(1) |˜u(x) − u(x)| = O(∆x)^r, whereu is smooth in I_j = (x_j−1 2, x_j+1

2) and u_x 6= 0 near Ij. (2) T.V.˜u(x) ≤ T.V.u(x)(1 + O(∆x))

7.2.1 Piecewise-constant reconstruction

Our equation is

u_t+ f (u)_x = 0 (2.4)

Following the algorithm, we have

(1) approximate u(t, x) by piecewise constant function, i.e., {Ujⁿ} represents the cell average of u(x, t_n) over (x_j−1

2, x_j+1 2).

x_j−1

2 x_j+1

2

∆x

(2) solve Riemann problem (u_j, u_j+1) on the edge x_j+1

2, its solutionu(x˜ _j+1

2, t),t_n < t < t_n+1can be found, which is a constant.

(3) integrate the equation (2.4) over(x_j−1

Riemann problem givesu(x, t) =˜

( uj ifx − xj+¹₂ < at, tn < t < tn+1 This is precisely the upwind scheme.

Example 2 Linear system

7.2. HIGH ORDER GODUNOV METHODS 93

Therefore the structure of the solution of Riemann problem is composed ofn waves ℓ1(Uj+1− U_j)r₁, · · · , ℓn(U_j+1− Uj)r_nwith left stateU_j and right stateU_j+1. Each wave propagate at

7.2.2 piecewise-linear reconstruction

a) high order approximation in smooth regions.

b) TVD or TVB or ENO Whenu has discontinuities or u_xchanges sign, we need to put a “limiter” to avoid oscillation ofu.˜

Example of limiters

7.2. HIGH ORDER GODUNOV METHODS 95 (2) Exact solver for small time step

Consider the linear advection equation

u_t+ au_x = 0.

Its graph is shown in Fig.(7.3).

Ifa < 0, then

0

Figure 7.3: The limiter of second order Godunov method

For System Case

u_t+ Au_x= 0 (2.5)

(1) Reconstruction

Constructu(x, t˜ _n) to be a piecewise linear function.

˜

We trace back along the characteristic curve to getu in half time step.

uⁿ⁺

7.3. MULTIDIMENSION 97

In another viewpoint, letuⁿ⁺¹²

j+¹₂,Lbe the solution of (2.5) with initial data = Then we solve (2.5) with(uⁿ⁺

1

There are two kinds of methods.

1. Splitting method.

2. Unsplitting method.

We consider two-dimensional case.

7.3.1 Splitting Method

We start from

u_t+ Au_x+ Bu_y = 0. (3.6)

This equation can be viewed as

u_t= (−A∂x− B∂y)u.

Then the solution operator is:

e^{−t(A∂x+B∂y)},

which can be approximate bye^−tA∂xe^−tB∂y for smallt. Let A = −A∂x, B = −B∂y, we have u = e^t(A+B)u0.

Considere^t(A+B),

e^t(A+B) = 1 + t(A + B) +t²

2(A²+ B²+ AB + BA) + · · · e^tB· e^tA = (1 + tB +t²

2B²+ · · · )(1 + tA + t²

2A²+ · · · )

= 1 + t(A + B) +t²

2(A²+ B²) + t²BA + · · · .^·.e^t(A+B)− e^tB· e^tA= t²(AB − BA

2 ) + O(t³).

Now we can design splitting method as:

Given{Ui,jⁿ},

1. For eachj, solve u_t+ Au_x= 0 with data {Ujⁿ} for ∆t step. This gives ¯U_i,jⁿ. U¯_i,jⁿ = U_i,jⁿ + ∆t

∆x(F (U_i−1,jⁿ , U_i,jⁿ) − F (Ui,jⁿ, U_i+1,jⁿ )) whereF (U, V ) is the numerical flux for u_t+ Au_x = 0.

2. For eachi, solve u_t+ Bu_y = 0 for ∆t step with data { ¯U_i,jⁿ}. This gives Ui,jⁿ⁺¹. U_i,jⁿ⁺¹= ¯U_i,jⁿ + ∆t

∆y(G( ¯U_i,j−1ⁿ , ¯U_i,jⁿ) − G( ¯U_i,jⁿ, ¯U_i,j+1ⁿ )) The error is first order in timen(∆t)² = O(∆t).

To reach higher order time splitting, we may approximatee^t(A+B)by polynomialsP (e^tA, e^tB) or rationalsR(e^tA, e^tB). For example, the Trotter product (or strang splitting) is given by

e^t(A+B)= e^12tAe^tBe^12tA+ O(t³).

Fort = n∆t,

e^t(A+B)u₀ = (e^12∆tAe^∆tBe^12∆tA) · · · (e^12∆tAe^∆tBe^12∆tA)(e^12∆tAe^∆tBe^12∆tA)u₀

= e^12∆tAe^∆tBe^∆tAe^∆tBe^∆tA· · · e^∆tAe^∆tBe^12∆tAu₀ Trotter product is second order.

7.3. MULTIDIMENSION 99 7.3.2 Unsplitting Methods

The PDE is

u_t+ f (u)_x+ g(u)_y = 0 (3.7)

Integrate this equation over(x_i−1 2, x_i+1 We consider Godunov type method.

1. Reconstruction

Finally, solve Riemann problemu_t+ Au_x= 0 with data





 Uⁿ⁺

1 2

i+¹₂,L,j

Uⁿ⁺

1 2

i+¹₂,R,j

.^·.Uⁿ⁺¹²

i+¹₂,j = Uⁿ⁺¹²

i+¹₂,L,j+ X

λ^x_k≥0

ℓ_k· δU_i+¹

2,jr_k

Chapter 8

Systems of Hyperbolic Conservation Laws

8.1 General Theory

We consider

u_t+ f (u)_x = 0, u =





 u₁ u₂ . ..

u_n





 f :Rⁿ→ Rⁿthe flux (1.1)

The system (1.1) is called hyperbolic if ∀u, the n × n matrix f^′(u) is diagonalizable with real eigenvaluesλ₁(u) ≤ λ2(u) ≤ · · · ≤ λn(u). Let us denote its left/right eigenvectors by ℓ_i(u)/r_i(u), respectively.

It is important to notice that the system is Galilean invariant, that is , the equation is unchanged under the transform:

t −→ λt, x −→ λx, ∀λ > 0.

This suggests we can look for special solution of the formu(^x_t).

We plugu(^x_t) into (1.1) to yield u^′· (−x

t²) + f^′(u)u^′·1 t = 0

=⇒ f^′(u)u^′ = x tu^′

This implies that there existsi such that u^′ = ri(u) and ^x_t = λi(u(^x_t)). To find such a solution, we first construct the integral curve ofr_i(u): u^′ = r_i(u). Let R_i(u₀, s) be the integral curve of r_i(u) passing throughu₀, and parameterized by its arclength. AlongR_i, the speedλ_ihas the variation:

d

dsλ_i(R_i(u₀, s)) = ∇λi· R^′i = ∇λi· ri. 101

We have the following definition.

Definition 1.15. Thei-th characteristic field is called 1. genuinely nonlinear if∇λⁱ(u) · rⁱ(u) 6= 0∀u.

2. linearly degenerate if∇λi(u) · ri(u) ≡ 0

3. nongenuinely nonlinear if∇λi(u) · ri(u) = 0 on isolated hypersurface inRⁿ.

For scalar equation, the genuine nonlinearity is equivalent to the convexity( or concavity) of the fluxf , linear degeneracy is f (u) = au, while nongenuine nonlinearity is nonconvexity of f . 8.1.1 Rarefaction Waves

When thei-th field is genuiely nonlinear, we define

R⁺_i (u₀) = {u ∈ Ri(u₀)|λi(u) ≥ λi(u₀)}.

Now supposeu₁ ∈ R⁺_i (u₀), we construct the centered rarefaction wave, denoted by (u₀, u₁):

(u₀, u₁)(x t) =





u₀ if ^x_t ≤ λi(u₀) u₁ if ^x_t ≥ λi(u₁)

u ifλ_i(u₀) ≤ ^x_t ≤ λi(u₁)andλ_i(u) = ^x_t It is easy to check this is a solution. We call(u₀, u₁) an i-rarefaction wave.

t

x λ_i(u0) λ_i(u₁)

u₀ u₁

λ_i(u) = ^x_t

λ_i

u₀

λ_i(u1) u₁

λ_i(u₀)

Figure 8.1: The integral curve ofu^′ = r_i(u) and the rarefaction wave.

8.1.2 Shock Waves

The shock wave is expressed by:

u(x t) =

 u₀ for ^x_t < σ u₁ for ^x_t > σ Then(u₀, u₁, σ) need to satisfy the jump condition:

f (u₁) − f (u0) = σ(u₁− u0). (1.2)

8.1. GENERAL THEORY 103 Lemma 8.3. (Local structure of shock waves)

1. The solution of (1.2) for (u, σ) consists of n algebraic curves passing through u0 locally, named them byS_i(u₀), i = 1, · · · , n.

2. S_i(u₀) is tangent to R_i(u₀) up to second order. i.e., S_i^(k)(u₀) = R^(k)_i (u₀), k = 0, 1, 2, here the derivatives are arclength derivatives.

3. σ_i(u₀, u) −→ λi(u₀) as u −→ u0, andσ_i^′(u₀, u₀) = ¹₂λ^′_i(u₀)

Proof. 1. LetS(u₀) = {u|f (u)−f (u0) = σ(u −u0) for some σ ∈ R}. We claim that S(u0) = Sn

i=1

S_i(u₀), where S_i(u₀) is a smooth curve passing through u₀ with tangent r_i(u₀) at u₀. Whenu is on S(u₀), rewrite the jump condition as

f (u) − f (u0) = [ Z ₁

0

f^′(u₀+ t(u − u0) dt](u − u0)

= A(u˜ ₀, u)(u − u0)

= σ(u − u0)

.^·. u ∈ S(u0) ⇐⇒ (u − u0) is an eigenvector of ˜A(u₀, u).

AssumeA(u) = f^′(u) has real and distinct eigenvalues λ₁(u) < · · · λn(u), ˜A(u₀, u) also has real and distinct eigenvalues ˜λ1(u0, u) < · · · < ˜λⁿ(u0, u), with left/right eigenvectors ℓ˜_i(u₀, u) and ˜r_i(u₀, u), respectively, and they converge to λ_i(u₀), ℓ_i(u₀), r_i(u₀) as u → u0

respectively. Normalize the eigenvectors: k˜rik = 1, ˜ℓ_ir˜_j = δ_ij. The vector which is parallel tor_ican be determined by

ℓ˜_k(u0, u)(u − u⁰) = 0 for k 6= i, k = 1, · · · , n.

Now we define

S_i(u₀) = {u|˜ℓk(u₀, u)(u − u0) = 0, k 6= i, k = 1, · · · , n}

We claim this is a smooth curve passing throughu₀. Choose coordinate systemr₁(u₀), · · · , rn(u₀).

Differential this equation ˜ℓ_k(u0, u)(u − u⁰) = 0 in rj(u0).

∂

∂rj

 

u=u0

(˜ℓ_k(u₀, u)(u − u0)) = ˜ℓ_k(u₀, u₀) · rj(u₀) = δ_jk,

which is a full rank matrix. By implicit function theorem, there exists unique free parameter smooth curveS_i(u₀) passing through u₀. ThereforeS(u₀) = Sⁿ

i=1

S_i(u₀).

2,3. R_i(u₀) = u₀ = S_i(u₀)

f (u) − f (u0) = σ_i(u₀, u)(u − u0) ∀u ∈ Si(u₀) Take arclength derivative alongS_i(u₀)

f^′(u)u^′ = σ_i^′(u − u0) + σ_iu^′andu^′ = S_i^′. Whenu −→ u0

f^′(u₀)S_i^′(u₀) = σ_i(u₀, u₀)S_i^′(u₀)

=⇒ S_i^′(u₀) = r_i(u₀) and σ_i(u₀, u₀) = λ_i(u₀).

Consider the second derivative.

(f^′′(u)u^′, u^′) + f^′(u)u^′′ = σ_i^′′(u − u0) + 2σ^′_i· u^′+ σ_iu^′′

Atu = u₀,u^′ = S_i^′(u₀) = R^′_i(u₀) = r_i(u₀) and u^′′= S_i^′′(u₀),

=⇒ (f^′′r_i, r_i) + f^′S_i^′′= 2σ^′_ir_i+ σ_iS_i^′′

On the other hand, we take derivative off^′(u)r_i(u) = λ_i(u)r_i(u) along R_i(u₀), then evaluate atu = u₀.

(f^′′r_i, r_i) + f^′(∇ri· ri) = λ^′_ir_i+ λ_i∇ri· ri, where∇ri· ri = R^′′_i.

=⇒ (f − λi)(S_i^′′− R^′′i) = (2σ_i^′− λ^′i)r_i

=⇒ 2σ_i^′ = λ^′_i

LetS^′′_i − Ri^′′=P

k

α_kr_k(u₀) = α_ir_i P

k6=i(λ_k− λi)α_kr_k = 0

=⇒ α_k= 0 ∀k 6= i andλ^′_i = 2σ^′_iatu₀

·.^· (R^′_i, R^′_i) = 1 (S_i^′, S_i^′) = 1 and (R^′′_i, R_i^′) = 0 (S_i^′′, S_i^′) = 0

.^·. (R^′′_i − Si^′′)⊥ri

.^·. α_i= 0 HenceR^′′_i = S_i^′′atu₀.

8.1. GENERAL THEORY 105

Let S_i⁻(u0) = {u ∈ Sⁱ(u0)|λⁱ(u) ≤ λ_i(u₀)}.

Ifu₁ ∈ S −i(u₀), define (u0, u1) =

 u₀ for ^x_t < σ_i(u₀, u₁) u₁ for ^x_t > σ_i(u₀, u₁) (u₀, u₁) is a weak solution.

u₀

S_i⁻

R_i⁺

R_i S_i

We propose the following entropy condition: (Lax entropy condition)

λ_i(u₀) > σ_i(u₀, u₁) > λ_i(u₁) (1.3) If thei-th characteristic field is genuinely nonlinear, then for u1 ∈ Si⁻(u0), and u1 ∼ u⁰, (1.3) is always valid. This follows easily fromλ_i = 2σ^′_i andσ_i(u₀, u₀) = λ_i(u₀). For u₁ ∈ Si⁻(u₀), we call the solution(u₀, u₁) i-shock or Lax-shock.

8.1.3 Contact Discontinuity (Linear Wave)

If∇λi(u) · ri(u) ≡ 0, we call the i-th characteristic field linearly degenerate (ℓ. dg.). In the case of scalar equation, this correspondf^′′= 0. We claim

Ri(u0) = Si(u0) and σi(u0, u) = λi(u0) for u ∈ Sⁱ(u0) or Ri(u0).

Indeed, alongRi(u0), we have

f^′(u)u^′ = λ_i(u)u^′.

andλ_i(u) is a constant λ_i(u₀) from the linear degeneracy. We integrate the above equation from u₀ tou along R_i(u₀), we get

f (u) − f (u⁰) = λi(u0)(u − u⁰).

This gives the shock condition. Thus,S_i(u₀) ≡ Ri(u₀) and σ(u, u₀) ≡ λi(u₀).

Homeworks.

(u₀, u₁) =

 u₀ ^x_t < σ_i(u₀, u₁) u₁ ^x_t > σ_i(u₀, u₁)

LetT_i(u₀) = R⁺_i (u₀) ∪ Si⁻(u₀) be called the i-th wave curve. For u₁ ∈ Ti(u₀), (u₀, u₁) is either a rarefaction wave, a shock, or a contact discontinuity.

Theorem 8.8. (Lax) For strictly hyperbolic system (1.1), if each field is either genuinely nonlin-ear or linnonlin-ear degenerate, then for u_L ∼ uR, the Riemann problem with two end states (u_l, u_R) has unique self-similar solution which consists ofn elementary waves. Namely, there exist u₀ = u_L, · · · , un= u_Rsuch that(u_i−1, u_i) is an i-wave.

Proof. Given (α₁, · · · , αn) ∈ Rⁿ, define u_i inductively u_i ∈ Ti(u_i−1), and the arclength of (u_i−1, u_i) on T_i= α_i.

u_i = f (u₀, α₁, · · · , αi) We want to findα₁, · · · , αnsuch that

uR= f (uL, α1, · · · , αⁿ).

u₀

u₁

T₁ T₂(u1) α₂

α₁

FirstuL= f (uL, 0, · · · , 0), as u^R= uL, (α1, · · · , αⁿ) = (0, · · · , 0) is a solution. When u^R ∼ u^L and{ri(u₀)} are independent,

∂

∂α_i  

α=0

f (u_L, 0, · · · , 0) = ri(u₀)andf ∈ C²

By Inverse function theorem, foru_R∼ uL, there exists uniqueα such that u_R= f (u_L, α). Unique-ness leaves as an exercise.

8.2 Physical Examples

8.2.1 Gas dynamics

The equations of gas dynamics are derived based on conservation of mass, momentum and energy.

Before we derive these equations, let us review some thermodynamics. First, the basic thermo variables are pressure (p), specific volume (τ ), called state variables. The internal energy (e) is a function ofp and τ . Such a relation is called a constitutive equation. The basic assumption are

∂e

∂p  

τ

> 0, ∂e

∂τ  

p

> 0 Sometimes, it is convinient to expressp as a function of (τ, e).

1-wave

2-wave

n-wave t

x

u₀ = uL u_n= u_R

u₁ u₂ u_n−1

8.2. PHYSICAL EXAMPLES 107 In an adiabetic process (no heat enters or losses), the first law of thermodynamics (conservation of energy) reads

de + pdτ = 0. (2.4)

This is called a Pfaffian equation mathematically. A functionσ(e, τ ) is called an integral of (2.4) if there exists a functionµ(e, τ ) such that

dσ = µ · (de + pdτ).

Thus, σ = constant represents a specific adiabetic process. For Pfaffian equation with only two independent variables, one can always find its integral. First, one can derive equation forµ: from

σ_e= µ and σ_τ = µp and usingσ_eτ = σ_{τ e}, we obtain the equation forµ:

µτ = (µp)e.

This is a linear first-order equation for µ. It can be solved by the method of characteristics in the region τ > 0 and e > 0. The solutions of µ and σ are not unique. If σ is a solution, so doesσ with d¯¯ σ = ν(σ)dσ for any function ν(σ). We can choose µ such that if two systems are in thermo-equilibrium, then they have the same value µ. In other words, µ is only a function of emperical temperature. We shall denote it by1/T . Such T is called the absolute temperature. The correspondingσ is called the physical entropy S. The relation dσ = µ(de + pdτ ) is re-expressed as

de = T dS − pdτ. (2.5)

For ideal gas, which satisfies the laws of Boyle and Gay-Lussac:

pτ = RT, (2.6)

whereR is the universal gas constant. From this and (2.5), treating S and τ as independent variables, one obtains

Re_S(S, τ ) + τ e_τ(S, τ ) = 0.

We can solve this linear first-order equation by the method of characteristics. We rewrite this equa-tion as a direcequa-tional differentiaequa-tion:

 R ∂

∂S + τ ∂

∂τ

 e = 0.

This means thate is constant along the characteristic curves Rdτ

dS = τ.

These characteristics can be integrated as

τ e^−S/R= φ.

Hereφ is a positive constant. The energy e(τ, S) is constant when τ e^−S/R is a constant. That is, e = h(φ) for some function h. We notice that h^′ < 0 because p = −(_∂τ^∂e)_S= −e^−S/Rh^′(τ H) > 0.

FromT = (_∂S^∂e)_τ = −_R¹h^′(φ) · φ, we see that T is a function of φ. In most cases, T is a decreasing function ofφ. We shall make this as an assumption. With this, we can invert the relation between T and φ and treat φ as a decreasing function of T . Thus, we can also view e as a function of T , say

Hence, we can choose two of as independent thermo variables and treat the rest three as dependent variables.

For instance,e is a linear function of T , i.e. e = cvT , where cv is a constant called specfic heat at constant volume. Such a gas is called polytropic gas. We can obtain

pτ = RT and e = c_vT = pτ

γ − 1 (2.7)

or in terms of entropy,

p = A(S)τ^−γ

8.2. PHYSICAL EXAMPLES 109 In general,c_p > c_v. Becausec_pis the amount of heat added to a system per unit mass at constant pressure. In order to maintain constant pressure, the volume has to expand (otherwise, pressure will increase), the extra amount of work due to expansion is supplied by the extra amount of heatc_p−cv. Next, we derive the equation of gas dynamics. Let us consider an arbitrary domain Ω ⊂ R³. The mass flux from outside to inside per unit time per unit areadS is −ρv·, where n is the outer normal ofΩ. Thus, the conservation of mass can be read as

d This holds for arbitraryΩ, hence we have

ρt+ div(ρ v) = 0. (2.8)

This is called the continuity equation.

Now, we derive momentum equation. Let us suppose the only surface force is from pressure (no viscous force). Then the momentum change inΩ is due to (i) the momentum carried in through boundary, (ii) the pressure force exerted on the surface, (iii) the body force. The first term is−ρvv·n, the second term is−pn. Thus, we have

d

Here, the notation∇ · ρv ⊗ v stands for a vector whoes ith component isP

j∂_j(ρvⁱv^j). The energy per unit volume isE = ¹₂ρ v² + ρe. The energy change in Ω per unit time is due to (i) the energy carried in through boundary (ii) the work done by the pressure from boundary, and (iii) the work done by the body force. The first term is−Ev · n. The second term is −pv · n. The third term is F · v. The conservation of energy can be read as

d By applying divergence theorem, we obtain the energy equation:

E_t+ div[(E + p)v] = ρF · v. (2.10)

In one dimension, the equations are (without body force)

ρ_t+ (ρu)_x = 0 (ρu)_t+ (ρu²+ p)_x = 0 (1

2ρu²+ e)_t+ [(1

2ρu²+ e + p)u]_x = 0.

Here, the unknowns are two thermo variable ρ and e, and one kinetic variable u. Other thermo variablep is given by the constitutive equation p(ρ, e).

8.2.2 Riemann Problem of Gas Dynamics We use(ρ, u, S) as our variables.

S. The eigenvalues and corresponding eigenvectors are

8.2. PHYSICAL EXAMPLES 111

Figure 8.2: The integral curve of the first and the third field on the(u, P ) phase plane.

OnR₂, which is a contact discontinuity,du = 0, dP = 0. Therefore u = u₀, P = P₀.

Let

m = ρ₀v₀= ρv

which is from the first jump condition. The second jump condition says that ρ₀v₀²+ P₀ = ρv²+ P

mv₀+ P₀ = mv + P

m = −P − P0

v − v0

= − P − P0

mτ − mτ0

whereτ = _ρ¹ is the specific volume.

.^·. m² = −^{P −P}_{τ −τ0}⁰ v − v0 = −^{P −P}_m ⁰

(u − u0)² = (v − v0)² = −(P − P0)(τ − τ0) The third one is

(1

2ρ₀v²₀+ ρ₀e₀+ P₀)v₀ = (1

2ρv²+ ρe + P )v

=⇒ 1

2v₀²+ e₀+ P₀τ₀ = 1

2v²+ e + P τ Byv²₀ = m²τ₀², v² = m²τ², m²= −^{P −P}_{τ −τ0}⁰,

=⇒ H(P, τ ) = e − e0+ P + P₀

2 (τ − τ0) = 0 Recalle = _γ−1^{P τ} . FromH(P, τ ) = 0,

P τ

γ − 1− P₀τ₀

γ − 1+ (P + P₀

2 )(τ − τ0) = 0.

Solve fotτ in terms of P, P₀, τ₀, then plug into

(u − u0)² = −(P − P0)(τ − τ0)

Setφ(P ) = (P − P0) vu

ut ^γ+12 τ₀ P +^γ−1_γ+1P0

Then

S₁ : u = u₀− φ0(P ) S₃ : u = u₀+ φ₀(P )

8.2. PHYSICAL EXAMPLES 113 Therefore,

T₁^(ℓ): u =

 u0− ψ⁰(P ) P < P0

u₀− φ0(P ) P ≥ P0

T₃^(ℓ): u =

 u₀+ ψ₀(P ) P > P₀ u0+ φ0(P ) P ≤ P⁰

T₁^(r): u =

 u₀− ψ0(P ) P > P₀ u₀− φ0(P ) P ≤ P0

T₃^(r): u =

 u0+ ψ0(P ) P < P0

u₀+ φ₀(P ) P ≥ P0

P

u (ℓ)

S₁⁻

R⁺₃

R⁺₁ S₃⁻

P

u S₁⁺

S₃⁺ R⁻₁

R⁻₃ (r)

Figure 8.3: The rarefaction waves ans shocks of 1,3 field on(u, P ) phase plane at left/right state.

Now we are ready to solve Riemann Problem with initial states(ρ_L, P_L, u_L) and (ρ_R, P_R, u_R).

Recall that in the second field,[P ] = [u] = 0.

ρ_I, P_I, u_I ρ_I, P_II, u_II

ρ_L, P_L, u_L ρ_R, P_R, u_R

PI = PII = P_∗ u_I = u_II = u_∗

u S₁

S₃

ℓ

P R₁

R₃ r

S₃ S₁

R₁ R₃

ρ_L, P_L, u_L ρ_R, P_R, u_R ρ_I ρ_II

P_∗, u_∗

The vaccum state P

u ℓ

r Solution must satisfyP > 0. If u_ℓ+ ℓ_ℓ is less thanu_r−ℓr, there is no solution.

For numerical, Godunov gives a procedure, ”Godunove iteration”. The algorithm to findP_∗: u_ℓ− fℓ(P ) = uI = uII = ur+ fr(P )

f₀(P ) =

 ψ₀(P ) P < P₀ φ₀(P ) P ≥ P0

We can solve this equation.

Godunov iteration is 

Z_R(u_∗− uR) = P_∗− PR

−ZL(u_∗− uL) = P_∗− PL. Where

Z_R =

rP_R τ_RΦ(P_∗

P_R) Z_L =

rP_L τ_LΦ(P_∗

P_L)

8.2. PHYSICAL EXAMPLES 115 and

Φ(w) =



 qγ+1

2 w +^γ−1₂ w > 1

2√γγ−1 1−w

1−w^γ−12γ w ≤ 1 This can be solved by Newton’s method.

Approximate Riemann SolverOur equation isut+f (u)x = 0 with Riemann data (uL, uR).

We look for middle state. Supposeu_L∼ uR, the original equation can be replaced by u_t+ f^′(u)u_x = 0

ut+ A(u)ux = 0 Chooseu =¯ u_L+ u_R

2 , solveu_t+ A(¯u)u_x = 0 with Riemann data u(x, 0) =

 u_L x < 0 u_R x > 0 Let λ_i, ℓ_i, r_ibe eigenvalues and eigenvectors ofA. Then

u(x

t) = u_L+ X

λi<^x_t

(ℓ_i· (uR− uL)) · ri.

One severe error in this approximate Riemann solver is that rarefaction waves are approximated by discontinuities. This will produce nonentropy shocks. To cure this problem, we expand such a linear discontinuity by a linear fan. Precisely, supposeλ_i(u_i−1) < 0, λ_i(u_i) > 0, this suggests that there exists rarefaction fan crossing ^x_t = 0. We then expand this discontinuity by a linear fan. At x/t = 0, we thus choose

um = (1 − α)ui−1+ αui, α = −λi(u_i−1)

λ_i(u_i) − λi(u_i−1).

Chapter 9

Kinetic Theory and Kinetic Schemes

9.1 Kinetic Theory of Gases 9.2 Kinetic scheme

Assume the equilibrium distribution isg₀(ξ). It should satisfy (i) momentum conditions, (ii) equa-tion of states (or flux condiequa-tion), (iii) positivity. That is, the moment condiequa-tion:

Z

g₀ψ_α(ξ) dξ = U_α

and flux condition: Z

g₀ξψ_α(ξ) dξ = F_α(U ).

For non-convex casef ,

117

在文檔中 Finite Difference Methods for Solving Differential Equations (頁 92-0)