8 Multidimensional Fluid Solvers and Higher-order Schemes

Higher-order Schemes

8.1 Short repeat: Godunov’s method and Riemann solvers

It is useful to introduce another interpretation of common finite-volume discretiza- tions of fluid dynamics, so-called Reconstruct-Evolve-Average (REA) schemes. We use this here for a short repeat of Godunov’s important method, and the way Rie- mann solvers come into play in it.

An REA update scheme of a hydrodynamical system discretized on a mesh can be viewed as a sequence of three steps:

1. Reconstruct: Using the cell-averaged quantities, this defines the run of these quantities everywhere in the cell. In the sketch, a piece-wise constant recon- struction is assumed, which is the simplest procedure one can use and leads to 1st order accuracy.

2. Evolve: The reconstructed state is then evolved forward in time by ∆t. In

Godunov’s approach, this is done by treating each cell interface as a piece-wise

constant initial value problem which is solved with the Riemann solver exactly

or approximately. This solution is valid for arbitrary times, but in practice one needs to limit the timestep ∆t such that the waves emanating from the opposite sides of a cell do not yet start to interact (at this point, the solution constructed for an isolated interface would become invalid).

3. Average: The evolved solution of the previous step is averaged at time ∆t to compute new average states U ⁿ⁺¹ in each cell in a conservative fashion. Then the whole cycle repeats again.

Fortunately, the averaging step does not need to be done explicitly. Instead it can be carried out by accounting for the fluxes that enter or leave the control volume of the cell. Consider the following integral over the region [x _i−

2

, x _i+

2

] × [t _n , t _n+1 ], corresponding to cell i and the course of the timestep ∆t = t n+1 − t n . We have

Z x

i+ 12

x

dx Z t

t

dt (∂ _t U + ∂ _x F) = 0 (8.1)

because the solution we evolve in Godunov’s approach is actually an exact solution of the hyperbolic equations. This then yields

Z x

i+ 12

x

U(x, t _n+1 )dx−

Z x

i+ 12

x

U(x, t _n )dx+

Z t

t

F(x _i+

2

, t)dt−

Z t

t

F(x _i−

2

, t)dt = 0.

(8.2) Notice that the first term is simply the average of U at the new time over the whole cell of length ∆x = x _i+

2

−x _i+

2

, hence this is the cell average U _n+1 we are looking for, multiplied by ∆x. It may seem as if the fluxes at the cell interfaces are needed as a function of time during the timestep in order to do the time integration of the fluxes.

However, in Godunov’s method we estimate the fluxes through Riemann problems with piece-wise constant initial state. Here the solutions are scale-invariant in the variable x/t, and the fluxes at the initial interfaces are in fact constant in time. We can hence write

F(x _i+

2

, t) = F ^? _i+

, (8.3)

where the asterisk is meant to indicate that the corresponding flux is the appropriate flux from the Riemann solver at the interface x _i+

2

. Equation (8.2) therefore gives:

U _n+1 = U _n + ∆t

∆x

F ^? _i−

− F ^? _i+

(8.4)

which is the update formula for the state of a cell in the REA schemes. What is

hence needed is a prescription to either exactly or approximately solve the Riemann

problem for a piece-wise linear left and right state that are brought into contact at

time t = t _n .

Formally, this can be written as

F ^? = F _Riemann (U _L , U _R ). (8.5) In practice, a variety of approximate Riemann solvers F _Riemann are commonly used in the literature, as discussed earlier in the lecture. For the ideal gas and for isothermal gas, it is also possible to solve the Riemann problem exactly, but not in closed form (i.e. the solution involves an iterative root finding of a non-linear equation).

There are now two main issues left:

• How can this be extended to multiple spatial dimensions?

• How can it be extended such that a higher order integration accuracy both in space and time is reached?

We’ll discuss these issues next.

8.2 Extending the 1D problem to multiple dimensions

So far, we have considered one-dimensional hyperbolic conservation laws of the form

∂ _t U + ∂ _x F(U) = 0. (8.6) For example, for isothermal gas with soundspeed c _s , the state vector U and flux vector F(U) are given as

U = ρ

ρu

, F =

ρu ρu ² + ρc ² _s

, (8.7)

where u is the velocity in the x-direction.

In three dimensions, the PDEs describing a fluid become considerably more in- volved. For example, the Euler equations for an ideal gas are given in explicit form as

∂ _t











 ρ ρu ρv ρw ρe











 + ∂ _x













ρu ρu ² + P

ρuv ρuw ρu(ρe + P )











 + ∂ _y













ρu ρuv ρv ² + P

ρvw ρv(ρe + P )











 + ∂ _z













ρw ρuw ρvw ρw ² + P ρw(ρe + P )













= 0

(8.8)

These equations are often written in the following notation

∂ _t U + ∂ _x F + ∂ _y G + ∂ _z H = 0. (8.9) Here the functions F(U), G(U) and H(U) give the flux vectors in the x-, y- and z-direction, respectively.

8.2.1 Dimensional splitting

Let us know consider the three dimensionally split problems derived from equation (8.9):

(1) ∂ _t U + ∂ _x F = 0, (8.10)

(2) ∂ _t U + ∂ _y G = 0, (8.11)

(3) ∂ _t U + ∂ _z H = 0, (8.12)

Note that the vectors appearing here have still the same dimensionality as in the full equations. They are augmented one-dimensional problems, i.e. the transverse variables still appear but spatial differentiation happens only in one direction. Be- cause of this, these additional transverse variables do not make the 1D problem more difficult compared to the ‘pure’ 1D problem considered earlier, but the fluxes appearing in them still need to be included.

Now let us assume that he have a method to solve/advance each of these one- dimensional problems. We can for example express this formally through time- evolution operators X (∆t), Y(∆t), and Z (∆t), which advance the solution by a timestep ∆t. Then the full time advance of the system can for example be approxi- mated by

U ⁿ⁺¹ ' Z (∆t)Y (∆t)X (∆t)U ⁿ (8.13)

This is one possible dimensionally split update scheme. In fact, this is the exact solution if (8.10)-(8.11) are the linear advection problem, but for more general non- linear equations it only provides a first order approximation. However, higher-order dimensionally split update schemes can also be easily constructed. For example, in two-dimensions

U ⁿ⁺¹ = 1

2 [X (∆t)Y (∆t) + Y(∆t)X (∆t)U ⁿ (8.14) and

U ⁿ⁺¹ = X (∆t/2)Y(∆t)X (∆t/2)U ⁿ (8.15) are second-order accurate. Similarly, for three dimensions the schemes

U ⁿ⁺² = X (∆t)Y (∆t)Z(∆t)Z(∆t)Y(∆t)X (∆t)U ⁿ (8.16)

U ⁿ⁺¹ = X (∆t/2)Y (∆t/2)Z(∆t)Y (∆t/2)X (∆t/2)U ⁿ (8.17)

are second-order accurate. As a general rule of thumb, the time evolution opera-

tors have to be applied alternatingly in reverse order to reach second-order accu-

racy. We see that the dimensionless splitting hence reduces the problem effectively

to a sequence of one-dimensional solution operations which are applied to multi- dimensional domains. Note that each one-dimensional operator leads to an update of U, and is a complete step for the corresponding augmented one-dimensional prob- lem. Gradients, etc., needed for the next step have then to be recomputed before the next time-evolution operator is applied. In practical applications of mesh codes, these one-dimensional solves often called sweeps.

8.2.2 Unsplit schemes

In an unsplit approach, all flux updates of a cell are applied simultaneously to a cell, not sequentially. This is for example illustrated in this 2D situation.

The unsplit update of cell i, j is then given by U _i,j ⁿ⁺¹ = U _i,j ⁿ + ∆t

∆x

F _i−

2

,j − F _i+

2

,j

+ ∆t

∆y

G _i,j−

2

− G _i,j+

2

(8.18)

Unsplit approaches can also be used for irregular shaped cells like those appearing in unstructured meshes. For example, integrating over a cell of volume V and denoting with U the cell average, we can write the cell update with the divergence theorem as

U ⁿ⁺¹ = U ⁿ − ∆t V

Z

F · dS, (8.19)

where the integration is over the whole cell surface, with outwards pointing face area vectors dS.

8.3 Higher order schemes

We should first clarify what we mean with higher order schemes. Loosely speaking, this refers to the convergence rate of a scheme in smooth regions of a flow. For example, if we know the analytic solution ρ(x) for some problem, and then obtain a numerical result ρ _i at a set of N points at locations x _i , we can ask what the typical error of the solution is. One possibility to quantify this would be through a L1 error norm. We can define

L1 = 1 N

X

i

|ρ _i − ρ(x _i )|. (8.20)

If we now measure this error quantitatively for different resolutions of the discretiza- tion, we want to find that L1 decreases with increasing N . In such a case our nu- merical scheme is converging, and provided we use sufficient numerical resources, we have a chance to get below any desired absolute error level. But the rate of convergence can be very different between different numerical schemes when applied to the same problem. If a method shows a L1 ∝ N ⁻¹ scaling, it is said to be first- order accurate; a doubling of the number of cells will then cut the error in half. A second-order method has L1 ∝ N ⁻² , meaning that a doubling of the number of cells can actually reduce the error by a factor of 4. This much better convergence rate is of course highly desirable. It is also possible to construct schemes with still higher convergence rates, but they quickly become much more complex and computation- ally involved, so that one one often reaches a point of diminishing return quickly.

But the extra effort one needs to make to go from first to second-order is often very

small, sometimes trivially small, so that one basically should always strive to do at

least this.

A first step in constructing a 2nd order extension of Godunov’s method is to replace the piece-wise constant with a piece-wise linear reconstruction. This involves that one first estimates gradients for each cell (usually by a simple finite difference formula). These are then slope-limited if needed such that the linear extrapolations of the cell states to the cell interfaces do not introduce new extrema. This slope- limiting procedure is quite important; it needs to be done to avoid that real fluid discontinuities introduce large spurious oscillations into the fluid.

Given slope limited gradients, for example ∇ρ for the density, one can then esti- mate the left and right state adjacent to an interface x _i+

2

by spatial extrapolation from the centers of the cells left and right from the interface:

ρ ^L _i+

= ρ i + (∇ρ) i

∆x

2 , (8.21)

ρ ^R _i+

= ρ _i+1 − (∇ρ) _i+1 ∆x

2 . (8.22)

The next step would in principle be to use these states in the Riemann solver. In doing this we will ignore the fact that our reconstruction has now a gradient over the cell; instead we still pretend that left and right of the interface, the fluid state can be taken as piece-wise constant as far as the Riemann solver is concerned. However, it turns out that the spatial extrapolation needs to be augmented also with a temporal extrapolation one half timestep into the future, such that the flux estimate is now effectively done in the middle of the timestep. This is necessary both to reach second-order accuracy in time and also for stability reasons.

Hence we really need to use ρ ^L _i+

2

= ρ _i + (∇ρ) _i ∆x 2 +

∂ρ

∂t

i

∆t

2 (8.23)

ρ ^R _i+

= ρ _i+1 − (∇ρ) _i+1 ∆x 2 +

∂ρ

∂t

i+1

∆t

2 (8.24)

for extrapolating to the interfaces. More generally, this has to be done for the whole state vector of the system, i.e.

U ^L _i+

= U i + (∂ x U) i

∆x

2 + (∂ t U) i

∆t

2 (8.25)

U ^R _i+

= U _i+1 − (∂ _x U) _i+1 ∆x

2 + (∂ _t U) _i+1 ∆t

2 . (8.26)

Note that here the quantity (∂ _x U) _i is a (slope-limited) estimate of the gradient in cell i, based on finite-differences plus a slope limiting procedure. Similarly, we somehow need to estimate the time derivative encoded in (∂ _x U) _i . How can this be done? One smart way to do this is to exploit the Jacobian matrix of the Euler equations. We can write the Euler equations as

∂ t U = −∂ x F(U) = − ∂F

∂U ∂ x U = −A(U) ∂ x U, (8.27) where A(U) is the Jacobian matrix. Using this, we can simply estimate the required time-derivative based on the spatial derivatives:

(∂ _t U) _i = −A(U _i ) (∂ _x U) _i (8.28) Hence the extrapolation can be done as

U ^L _i+

= U i + ∆x

2 − ∆t

2 A(U i )

(∂ x U) i (8.29)

U ^R _i+

= U i+1 +

− ∆x 2 − ∆t

2 A(U i+1 )

(∂ x U) i+1 (8.30) This procedure defines the so-called MUSCL-Hancock scheme, which is a 2nd-order accurate extension of Godunov’s method.

Higher-order extensions such as the piece-wise parabolic method (PPM) start out

with a higher order polynomial reconstruction. In the case of PPM, parabolic shapes

are assumed in each cell instead of piece-wise linear states. The reconstruction is

still guaranteed to be conservative, however, i.e. the integral underneath the recon-

struction recovers the total values of the conserved variables individually in each

cell. So-called ENO and WENO schemes use yet higher-order polynomials to recon-

struct the state in a conservative fashion. Here many more cells in the environment

need to be considered (i.e. the stencil of these methods is much larger) to robustly

determine the coefficients of the reconstruction. This can for example involve a

least-square fitting procedure.