LINEBURG

ńņšąķčöą 13 |

including a computer graphics simulation of water, and ļ¬nally a chapter on

fully two-phase incompressible ļ¬‚ow. We wrap up this incompressible ļ¬‚ow

material with a chapter on incompressible ļ¬‚ames, including applications in

Part IV computer graphics, and a chapter on techniques for coupling a compressible

ļ¬‚uid to an incompressible ļ¬‚uid. Finally, we turn our attention to heat ļ¬‚ow

Computational Physics with a chapter discussing the heat equation and a chapter discussing level

set techniques for solving Stefan problems.

While the ļ¬eld of computational ļ¬‚uid dynamics is quite broad, a large por-

tion of it is dedicated to computations of compressible ļ¬‚ow, incompressible

ļ¬‚ow, and heat ļ¬‚ow. In fact, these three classes of problems can be thought

of as the basic model problems for hyperbolic, elliptic, and parabolic par-

tial diļ¬erential equations. Volumes have been ļ¬lled with both broad and

detailed work dedicated to these important ļ¬‚ow ļ¬elds; see, for example,

[86], [87], [7], [8], and the references therein. One might assume that there

is little left to add to the understanding of these problems. In fact, many

papers are now concerned with smaller details, e.g., the number of grid

points in a shock wave or contact discontinuity. Other papers are devoted

to rarely occurring pathologies, e.g., slow-moving shock waves and shock

overheating at solid wall boundaries.

One area where signiļ¬cant new ideas are still needed is multicomponent

ļ¬‚ow, for example, multicomponent compressible ļ¬‚ow where diļ¬erent ļ¬‚u-

ids have diļ¬erent equations of state, multicomponent incompressible ļ¬‚ow

where diļ¬erent ļ¬‚uids have diļ¬erent densities and viscosities with surface

tension forces at the interface, and Stefan-type problems where the individ-

ual materials have diļ¬erent thermal conductivities. These problems have

interfaces separating the diļ¬erent materials, and special numerical tech-

niques are required to treat the interface. The most commonly used are

front tracking, volume of ļ¬‚uid, and level set methods. In the next three

chapters we discuss the use of level set methods for the canonical equations

of computational ļ¬‚uid dynamics.

The ļ¬rst chapter discusses basic one-phase compressible ļ¬‚ow, and the

subsequent chapter shows how the level set method can be used for two-

phase compressible ļ¬‚ow problems where the equations of state diļ¬er across

the interface. Then, the use of level set techniques for deļ¬‚agration and det-

onation discontinuities is discussed in the third chapter. The fourth chapter

introduces techniques for coupling an Eulerian grid to a Lagrangian grid.

This is useful, for example, in compressible solid/ļ¬‚uid structure problems.

After this, we turn our attention to incompressible ļ¬‚ow with a chapter

focused on the basic one-phase equations including a computer graphics

simulation of smoke, a chapter on level set techniques for free surface ļ¬‚ows

150 14. Hyperbolic Conservation Laws and Compressible Flow

system of conservation laws. These also form the basis for their numerical

modeling.

A conserved quantity, such as mass, can be transported by convective

14 or diļ¬usive ļ¬‚uxes. The distinction is that diļ¬usive ļ¬‚uxes are driven by

gradients, while convective ļ¬‚uxes persist even in the absence of gradi-

ents. For most ļ¬‚ows where compressibility is important, e.g., ļ¬‚ows with

Hyperbolic Conservation Laws and shock waves, one needs to model only the convective transport and can ig-

nore diļ¬usion (mass diļ¬usion, viscosity and thermal conductivity) as well

Compressible Flow as the source terms (such as chemical reactions, atomic excitations, and

ionization processes). Moreover, convective transport requires specialized

numerical treatment, while diļ¬usive and reactive eļ¬ects can be treated with

standard numerical methods, such as simple central diļ¬erencing, that are

independent of those for the convective terms. Stiļ¬ reactions, however, can

present numerical diļ¬culties; see, for example, Colella et al. [50]. Conserva-

tion laws with only convective ļ¬‚uxes are known as hyperbolic conservation

laws. A vast array of physical phenomena are modeled by such systems,

e.g., explosives and high-speed aircraft.

The important physical phenomena exhibited by hyperbolic conserva-

tion laws are bulk convection, waves, contact discontinuities, shocks, and

rarefactions. We brieļ¬‚y describe the physical features and mathematical

model equations for each eļ¬ect, and most importantly note the implica-

We begin this chapter by addressing general systems of hyperbolic conserva- tions they have on the design of numerical methods. For more details on

tion laws including numerical techniques for computing accurate solutions numerical methods for conservation laws, see, e.g., LeVeque [105] and Toro

to them. Then we discuss the equations for one-phase compressible ļ¬‚ow as [164].

an example of a system of hyperbolic conservation laws.

14.1.1 Bulk Convection and Waves

14.1 Hyperbolic Conservation Laws Bulk convection is simply the bulk movement of matter, carrying it from

one spot to another, like water streaming from a hose. Waves are small-

A continuum physical system is described by the laws of conservation of amplitude smooth disturbances that are transmitted through the system

mass, momentum, and energy. That is, for each conserved quantity, the rate without any bulk transport like ripples on a water surface or sound waves

of change of the total amount in some region is given by its ļ¬‚ux (convective through air. Whereas convective transport occurs at the gross velocity

or diļ¬usive) through the region boundary, plus whatever internal sources of the material, waves propagate at the āspeed of soundā in the system

exist. The integral form of this conservation law is (relative to the bulk convective motion of the system). Waves interact by

superposition, so that they can either cancel out (interfere) or enhance each

d

f (u) Ā· dA =

u dV + s(u) dV (14.1) other.

dt ā„¦ ā‚ā„¦ ā„¦

The simplest model equation that describes bulk convective transport is

where u is the density of the conserved quantity, f (u) is the ļ¬‚ux, and s(u) the linear convection equation

is the source rate. By taking ā„¦ to be an inļ¬nitesimal volume and applying

ut + v Ā· u = 0, (14.3)

the divergence theorem, we get the diļ¬erential form of the conservation

law, where v is a constant equal to the convection velocity. The solution to this

is simply that u translates at the constant speed v. This same equation can

ā‚u

+ Ā· f (u) = s(u), (14.2) also be taken as a simple model of wave motion if u is a sine wave and v is

ā‚t

interpreted as the speed of sound. The linear convection equation is also an

which is the basis for the numerical modeling of continuum systems. A important model for understanding smooth transport in any conservation

physical system can be described by a system of such equations, i.e., a

14.1. Hyperbolic Conservation Laws 151 152 14. Hyperbolic Conservation Laws and Compressible Flow

14.1.3 Shock Waves

law. As long as f is smooth and u has no jumps in it, the general scalar

conservation law

A shock is a spatial jump in material properties, like pressure and tem-

perature, that develops spontaneously from smooth distributions and

Ā· f (u) = 0

ut + (14.4)

then persists. The shock jump is self-forming and also self-maintaining.

This is unlike a contact discontinuity, which must be put in the sys-

can be rewritten as

tem initially and will not resharpen itself if it is smeared out by some

ut + f (u) Ā· u = 0, (14.5) other process. Shocks develop through a feedback mechanism in which

strong impulses move faster than weak ones, and thus tend to steepen

where f (u) acts as a convective velocity. That is, locally in smooth parts themselves up into a āstepā proļ¬le as they travel through the system. Fa-

of the ļ¬‚ow, a conservation law behaves like bulk convection with velocity miliar examples are the āsonic boomā of a jet aircraft and the ābangā

f (u). This is called the local characteristic velocity of the ļ¬‚ow. from a gun. These sounds are our perceptions of a sudden jump in air

Bulk convection and waves are important because they imply that sig- pressure.

nals propagate in deļ¬nite directions at deļ¬nite speeds. This is in contrast The simplest model equation that describes shock formation is the one-

to diļ¬usion, which propagates signals in all directions at arbitrarily large dimensional Burgersā™ equation

speeds depending on the severity of the driving gradients. Thus we antici-

u2

pate that suitable numerical methods for hyperbolic systems will also have

ut + = 0, (14.6)

directional biases in space, which leads to the idea of upwind diļ¬erencing 2 x

and a deļ¬nite relation between the space and time steps (discrete propa-

which looks like the convection equation (14.3) with a nonconstant convec-

gation speed), which will roughly be that the discrete propagation speed

tive speed of u, i.e., ut + uux = 0. Thus larger u values move faster, and

āx/āt must be at least as large as the physical propagation speeds (char-

they will overtake smaller values. This ultimately results in the develop-

acteristic speeds) in the problem. The general form of this relation is called

ment of, for example, a right-going shock if the initial data for u constitute

the Courant-Friedrichs-Lewy (CFL) restriction.

any positive, decreasing function.

Wave motion and bulk convection do not create any new sharp features

Shocks move at a speed that is not simply related to the bulk ļ¬‚ow speed

in the ļ¬‚ow. The other remaining phenomena are all special because they

or characteristic speed, and they are not immediately evident from exam-

involve discontinuous jumps in the transported quantities. Because smooth

ining the ļ¬‚ux, in contrast to contacts. Shock speed is controlled by the

features can be accurately represented by a polynomial interpolation, we

diļ¬erence between inļ¬‚ux and outļ¬‚ux of conserved quantity into the region.

expect to be able to develop numerical methods of extremely high accuracy

Speciļ¬cally, suppose a conserved quantity u with conservation law

for the wave and convective eļ¬ects. Conversely, since jump functions are

poorly represented by polynomials, we expect little accuracy and perhaps

ut + f (u)x = 0 (14.7)

great diļ¬culty in numerically approximating the discontinuous phenomena.

has a step function proļ¬le with constant values extending both to the left,

uL , and to the right, uR , with a single shock jump transition in between

14.1.2 Contact Discontinuities moving with speed s. Then the integral form of the conservation law (14.1),

applied to any interval containing the shock, gives the relation

A contact discontinuity is a persistent discontinuous jump in mass density

moving by bulk convection through a system. Since there is negligible mass s(uR ā’ uL ) = f (uR ) ā’ f (uL ), (14.8)

diļ¬usion, such a jump persists. These jumps usually appear at the point

which is just another statement that the rate at which u appears, s(uR ā’

of contact of diļ¬erent materials; for example, a contact discontinuity sep-

arates oil from water. Contacts move at the local bulk convection speed, uL ), in the interval of interest is given by the diļ¬erence in ļ¬‚uxes across

or more generally the characteristic speed, and can be modeled by using the interval. Thus we see that the proper speed of the shock is directly

step-function initial data in the bulk convection equation (14.3). Since con- determined by conservation of u via the ļ¬‚ux f . This has an important

tacts are simply a bulk convection eļ¬ect, they retain any perturbations implication for numerical method design; namely, a numerical method will

they receive. Thus we expect contacts to be especially sensitive to numeri- ācaptureā the correct shock speeds only if it has āconservation form,ā i.e.,

cal methods; i.e., any spurious alteration of the contact will tend to persist if the rate of change of u at some node is the diļ¬erence of ļ¬‚uxes that are

and accumulate. accurate approximations of the real ļ¬‚ux f .

14.1. Hyperbolic Conservation Laws 153 154 14. Hyperbolic Conservation Laws and Compressible Flow

14.2 Discrete Conservation Form

The self-sharpening feature of shocks has two implications for numerical

methods. First, it means that even if the initial data are smooth, steep

gradients and jumps will form spontaneously. Thus, our numerical method To ensure that shocks and other steep gradients are captured by the scheme,

must be prepared to deal with shocks even if none are present in the ini- i.e., that they move at the right speed even if they are unresolved, we must

tial data. Second, there is a beneļ¬cial eļ¬ect from self-sharpening, because write the equation in a discrete conservation form. That is, a form in which

modest numerical errors introduced near a shock (smearing or small os- the rate of change of conserved quantities is equal to a diļ¬erence of ļ¬‚uxes.

cillations) will tend to be eliminated, and will not accumulate. The shock This form guarantees that we discretely conserve the total amount of the

is naturally driven toward its proper shape. Because of this, computing states u (e.g., mass, momentum, and energy) present, in analogy with the

strong shocks is mostly a matter of having a conservative scheme in order integral form given by equation (14.1). More important, this can be shown

to get their speed correct. to imply that steep gradients or jumps in the discrete proļ¬les propagate at

the physically correct speeds; see, for example, LeVeque [105].

Usually, conservation form is derived for control volume methods, that

is methods that evolve cell average values in time rather than nodal values.

14.1.4 Rarefaction Waves In this approach, a grid node xi is assumed to be the center of a grid cell

A rarefaction is a discontinuous jump or steep gradient in properties that (xiā’1/2 , xi+1/2 ), and we integrate the conservation law (14.7) across this

dissipates as a smooth expansion. A common example is the jump in air control volume to obtain

pressure from outside to inside a balloon, which dissipates as soon as the

ut + f (ui+1/2 ) ā’ f (uiā’1/2 ) = 0,

ĀÆ (14.9)

balloon is burst and the high-pressure gas inside is allowed to expand.

Such an expansion also occurs when the piston in an engine is rapidly

where u is the integral of u over the cell, and uiĀ±1/2 are the (unknown)

ĀÆ

pulled outward from the cylinder. The expansion (density drop) associated

values of u at the cell walls. This has the desired conservation form in that

with a rarefaction propagates outward at the sound speed of the system,

the rate of change of the cell average is a diļ¬erence of ļ¬‚uxes. The diļ¬culty

relative to the underlying bulk convection speed. A rarefaction can be

with this formulation is that it requires transforming between cell averages

modeled by Burgersā™ equation (14.6) with initial data that start out as

of u (which are directly evolved in time by the scheme) and cell wall values

a steep increasing step. This step will broaden and smooth out during the

of u (which must be reconstructed) to evaluate the needed ļ¬‚uxes. While this

evolution.

is manageable in one spatial dimension, in higher-dimensional problems the

A rarefaction tends to smooth out local features, which is generally ben-

series of transformations necessary to convert the cell averages to cell wall

eļ¬cial for numerical modeling. It tends to diminish numerical errors over

quantities becomes increasingly complicated. The distinction between cell

time and make the solution easier to represent by polynomials, which form

average and midpoint values can be ignored for schemes whose accuracy

the basis for our numerical representation. However, a rarefaction often

is no higher than second order (e.g., TVD schemes), since the cell average

connects to a smooth (e.g., constant) solution region and this results in

and the midpoint value diļ¬er by only O(āx2 ).

a ācorner,ā which is notoriously diļ¬cult to capture accurately. The main

Shu and Osher [150, 151] proposed a fully conservative ļ¬nite diļ¬erence

numerical problem posed by rarefactions is that of initiating the expansion.

scheme on uniform grids that directly evolves nodal values (as opposed

If the initial data is are perfect, symmetrical step, such as u(x) = sign(x),

to the cell average values) forward in time. They deļ¬ned a numerical ļ¬‚ux

it may be āstuckā in this form, since the steady-state Burgersā™ equation is

function F by the property that the real ļ¬‚ux divergence is a ļ¬nite diļ¬erence

satisļ¬ed identically (i.e., the ļ¬‚ux u2 /2 is constant everywhere, and similarly

of numerical ļ¬‚uxes

in any numerical discretization). However, local analysis can identify this

F (x + āx/2) ā’ F (x ā’ āx/2)

stuck expansion, because the characteristic speed u on either side points

f (u)x = (14.10)

away from the jump, suggesting its potential to expand. In order to get āx

the initial data unstuck, a small amount of smoothing must be applied to

at every point x. We call F the numerical ļ¬‚ux, since we require it in our

introduce some intermediate-state values that have a nonconstant ļ¬‚ux to

numerical scheme, and also to distinguish it from the closely related āphys-

drive expansion. In numerical methods this smoothing applied at a jump

ical ļ¬‚uxā f (u). It is not obvious that the numerical ļ¬‚ux function exists, but

where the eļ¬ective local velocity indicates expansion should occur is called

from relationship (14.10) one can solve for its Taylor expansion to obtain

an āentropy ļ¬x,ā since it allows the system to evolve from the artiļ¬cial

low entropy initial state to the proper increased entropy state of a free (āx)2 7(āx)4

F = f (u) ā’ f (u)xxxx ā’ Ā· Ā· Ā· ,

f (u)xx + (14.11)

expansion.

24 5760

14.3. ENO for Conservation Laws 155 156 14. Hyperbolic Conservation Laws and Compressible Flow

which shows that the physical and numerical ļ¬‚ux functions are the same only nodal values of the conserved variables. Their method is faster and eas-

to second-order accuracy in āx. Thus, a ļ¬nite diļ¬erence discretization can ier to implement than the cell-averaged formulation. In addition, the ļ¬nite

be based on diļ¬erence ENO method extends to higher dimensions in a dimension-by-

dimension fashion, so that the one-dimensional method applies unchanged

F (x + āx/2) ā’ F (x ā’ āx/2)

ut + =0 (14.12) to higher-dimensional problems. We emphasize that this is not dimensional

āx

splitting in time, which has accuracy limitations unlike the dimension-by-

to evolve point values of u forward in time using numerical ļ¬‚ux functions F dimension approach. Shu and Osher also use the method of lines for time

at the cell walls. integration, decoupling the time and space discretizations.

We consider the treatment of a one-dimensional contact discontinuity

to illustrate how the method works. Assuming that the time evolution

takes place exactly, each time step āt should rigidly translate the spatial

14.3 ENO for Conservation Laws

proļ¬le by the amount vāt as governed by equation (14.3). Spatially, the

contact is initially represented by a discrete step function, i.e., nodal values

14.3.1 Motivation that are constant at one value uL on nodes x1 , . . . , xJ , and constant at a

Essentially nonoscillatory (ENO) methods were developed to address the diļ¬erent value uR on all remaining nodes xJ+1 , . . . , xN . To update the

special diļ¬culties that arise in the numerical solution of systems of non- value ui in time at a given node xi , we ļ¬rst reconstruct the graph of a

linear conservation laws. Numerical methods for these problems must be function u(x) near xi by interpolating nearby nodal u values, shift that

able to handle steep gradients, e.g., shocks and contact discontinuities, u(x) graph spatially by vāt (the exact time evolution), and then reevaluate

that may develop spontaneously and then persist in these ļ¬‚ows. Classical it at the node xi to obtain the updated ui . We require our local interpolant

numerical schemes had a tendency either to produce large spurious oscil- be smooth at the point xi , since in actual practice we are going to use it to

lations near steep gradients or to greatly smear out both these gradients evaluate the derivative term ux there. The simplest symmetric approach to

and the ļ¬ne details of the ļ¬‚ow. The primary goal of the ENO eļ¬ort was smooth interpolation near a node xi is to run a parabola through the nodal

to develop a general-purpose numerical method for systems of conservation data at xiā’1 , xi , and xi+1 . This interpolation is an accurate reconstruction

laws that has high accuracy (e.g., third order) in smooth regions and cap- of u(x) in smooth regions, where it works well. However, near the jump

tures the motion of unresolved steep gradients without creating spurious between xJ and xJ+1 the parabola will signiļ¬cantly overshoot the nodal u

data by an amount comparable to the jump uL ā’ uR , and this overshoot

oscillations and without the use of problem-dependent ļ¬xes or tunable pa-

rameters. The philosophy underlying the ENO methods is simple: When will show up in the nodal values once the shift is performed. Successive

reconstructing a proļ¬le for use in a convective ļ¬‚ux term, one should not time steps will further enhance these spurious oscillations. This approach

use high-order polynomial interpolation across a steep gradient in the data. corresponds to standard central diļ¬erencing.

Such an interpolant would be highly oscillatory and ultimately corrupt the To avoid the oscillations from parabolic interpolation, we could instead

computed solution. ENO methods use an adaptive polynomial interpola- use a smooth linear interpolation near xi , noting that there are two linear

tion constructed to avoid steep gradients in the data. The polynomial is also interpolants to choose from, namely the line through the data at nodes xi

biased to use data from the direction of information propagation (upwind) and xiā’1 , and the line through the data at xi and xi+1 . The direction of in-

for physical consistency and stability. formation propagation determines which should be used. If the convection

The original ENO schemes developed by Harten et al. [81] were based speed v is positive, the data are moving from left to right, and we use xi

on the conservative control volume discretization of the equations, which and xiā’1 . This linear interpolation based on upwind nodes will not intro-

yields discrete evolution equations for grid cell averages of the conserved duce any new extrema in u as long as the shift vāt is less than the width

of the interval āx = xi ā’ xiā’1 , which is exactly the CFL restriction on the

quantities, e.g. mass, momentum, and energy. This formulation has the

disadvantage of requiring complicated transfers between cell averages and time step. The main problem with the linear upwind biased interpolant is

cell center nodal values in the algorithm. In particular, the transfer pro- that it has low accuracy smearing out the jump over more and more nodes.

cess becomes progressively more complicated in one, two, and three spatial If we naively go to higher accuracy by using a higher-order upwind biased

dimensions. The formulation also results in space and time discretizations interpolant, such as running a parabola through xi , xiā’1 , and xiā’2 to ad-

that are coupled in a way that becomes complicated for higher-order accu- vance ui , we run into the spurious oscillation problem again. In particular,

rate versions. To eliminate these complications, Shu and Osher [150, 151] at nodes xJ+1 and xJ+2 , this upwind parabola will interpolate across the

developed a conservative ļ¬nite diļ¬erence form of the ENO method that uses jump and thus have large overshoots. By forcing the parabola to cross a

14.3. ENO for Conservation Laws 157 158 14. Hyperbolic Conservation Laws and Compressible Flow

jump, it no longer reļ¬‚ects the data on the interval that will be arriving at where the second equality sign comes from

xJ+1 (or xJ+2 ) during the next time step. i

xi+1/2 xj+1/2

The motivation for ENO is that we must use a higher-degree polynomial H(xi+1/2 ) = h(y) dy = h(y) dy (14.18)

interpolant to achieve more accuracy, and it must involve the immediate up- xā’1/2 xjā’1/2

j=0

wind node to properly represent the propagation of data. But we must also i

avoid polluting this upwind data with spurious oscillations that come from = x f (u(xj )), (14.19)

interpolating across jumps. Thus, the remaining interpolation nodes (after j=0

the ļ¬rst upwind point) are chosen based on smoothness considerations. In

and the higher divided diļ¬erences are

particular, this approach will, if at all possible, not run an interpolant across

f (u(xi+1 )) ā’ f (u(xi )) 11

a jump in the data. However, very small interpolation overshoots do occur 2

Di+1/2 H = = Di+1/2 f, (14.20)

near extrema in the nodal data, as they must, since any smooth function 2x 2

will slightly overshoot its values as sampled at discrete points near extrema. 12

3

Di H = Di f, (14.21)

This is the sense in which the method is only essentially nonoscillatory. 3

continuing in that manner.

According to the rules of polynomial interpolation, we can take any path

14.3.2 Constructing the Numerical Flux Function

along the divided diļ¬erence table to construct H, although not all paths

We deļ¬ne the numerical ļ¬‚ux function through the relation give good results. ENO reconstruction consists of two important features.

1

First, choose Di H in the upwind direction. Second, choose higher-order

Fi+1/2 ā’ Fiā’1/2

f (ui )x = (14.13) divided diļ¬erences by taking the smaller in absolute value of the two pos-

āx

sible choices. Once we construct H(x), we evaluate H (xi+1/2 ) to get the

as in equation (14.12). To obtain a convenient algorithm for computing this

numerical ļ¬‚ux Fi+1/2 .

numerical ļ¬‚ux function, we deļ¬ne h(x) implicitly through the equation

x+ x/2

1

14.3.3 ENO-Roe Discretization (Third-Order Accurate)

f (u(x)) = h(y) dy, (14.14)

x xā’ x/2

For a speciļ¬c cell wall, located at xi0 +1/2 , we ļ¬nd the associated numerical

and note that taking a derivative on both sides of this equation yields

ļ¬‚ux function Fi0 +1/2 as follows. First, we deļ¬ne a characteristic speed

x/2) ā’ h(x ā’

h(x + x/2) Ī»i0 +1/2 = f (ui0 +1/2 ), where ui0 +1/2 = (ui0 + ui0 +1 )/2 is deļ¬ned using a

f (u(x))x = , (14.15)

x standard linear average. Then, if Ī»i0 +1/2 > 0, set k = i0 . Otherwise, set

k = i0 + 1. Deļ¬ne

which shows that h is identical to the numerical ļ¬‚ux function at the cell

walls. That is, FiĀ±1/2 = h(xiĀ±1/2 ) for all i. We calculate h by ļ¬nding its 1

Q1 (x) = (Dk H)(x ā’ xi0 +1/2 ). (14.22)

primitive 2 2 2

If |Dkā’1/2 H| ā¤ |Dk+1/2 H|, then c = Dkā’1/2 H and k = k ā’ 1. Otherwise,

x

2

c = Dk+1/2 H and k = k. Deļ¬ne

H(x) = h(y) dy (14.16)

xā’1/2

Q2 (x) = c(x ā’ xkā’1/2 )(x ā’ xk+1/2 ). (14.23)

using polynomial interpolation, and then take a derivative to get h. Note

3 3 3 3

If |Dk H| ā¤ |Dk +1 H|, then c = Dk H. Otherwise, c = Dk +1 H. Deļ¬ne

that we do not need the zeroth-order divided diļ¬erences of H that vanish

with the derivative. Q3 (x) = c (x ā’ xk ā’ xk ā’ xk

ā’1/2 )(x +1/2 )(x +3/2 ). (14.24)

0

The zeroth order divided diļ¬erences, Di+1/2 and all higher-order even

Then

divided diļ¬erences of H exist at the cell walls and have the subscript iĀ±1/2.

1

The ļ¬rst order divided diļ¬erences Di and all higher-order odd divided Fi0 +1/2 = H (xi0 +1/2 ) = Q1 (xi0 +1/2 )+Q2 (xi0 +1/2 )+Q3 (xi0 +1/2 ), (14.25)

diļ¬erences of H exist at the grid points and will have the subscript i. The

which simpliļ¬es to

ļ¬rst-order divided diļ¬erences of H are

1

3(i0 ā’ k )2 ā’ 1 ( x)2 .

Fi0 +1/2 = Dk H + c (2(i0 ā’ k) + 1) x+c

H(xi+1/2 ) ā’ H(xiā’1/2 )

1

Di H = = f (u(xi )), (14.17) (14.26)

x

14.3. ENO for Conservation Laws 159 160 14. Hyperbolic Conservation Laws and Compressible Flow

14.4 Multiple Spatial Dimensions

14.3.4 ENO-LLF Discretization (and the Entropy Fix)

The ENO-Roe discretization can admit entropy-violating expansion shocks

In multiple spatial dimensions, the ENO discretization is applied inde-

near sonic points. That is, at a place where a characteristic velocity changes

pendently using a dimension-by-dimension discretization. For example,

sign (a sonic point) it is possible to have a stationary expansion shock so-

consider a two-dimensional conservation law

lution with a discontinuous jump in value. If this jump were smoothed out

even slightly, it would break up into an expansion fan (i.e., rarefaction) ut + f (u)x + g(u)y = 0 (14.30)

and dissipate, which is the desired physical solution. For a speciļ¬c cell wall

on a rectangular 2D grid. Here, we sweep through the grid from bottom

xi0 +1/2 , if there are no nearby sonic points, then we use the ENO-Roe dis-

to top performing ENO on 1D horizontal rows of grid points to evaluate

cretization. Otherwise, we add high-order dissipation to our calculation of

the f (u)x term. The g(u)y term is evaluated in a similar manner, sweeping

Fi0 +1/2 to break up any entropy-violating expansion shocks. We call this

through the grid from left to right performing ENO on 1D vertical rows

entropy-ļ¬xed version of the ENO-Roe discretization ENO-Roe ļ¬x or just

of grid points. Once we have a numerical approximation to each of the

ńņšąķčöą 13 |

Copyright © Design by: Sunlight webdesign