Numerical Schemes for Hyperbolic Balance Laws – Applications to Fluid Flow Problems

Balance laws arise from many areas of engineering practice specifically from the fluid mechanics. Many numerical methods for the solution of these balanced laws were developed in recent decades. The numerical methods are based on two views: solving hyperbolic PDE with a nonzero source term (the obvious description of the central and central-upwind schemes; (Kurganov & Levy, 2002; LeVeque, 2004)) or solving the augmented quasilinear nonconservative formulation (Gosse, 2001; Le Floch & Tzavaras, 1999; Parès, 2006). Furthermore, the methods can be interpreted using flux-difference splitting (or flux-vector splitting), or by selecting adaptive intervals and the transformation to the semidiscrete form (for example (Kurganov & Petrova, 2000)). We prefer the augmented quasilinear nonconservative formulation solved by the flux-difference splitting in our text. We try to formulate the methods in the most general form. The range of this text does not give the complete overview of currently used methods.


Introduction
Balance laws arise from many areas of engineering practice specifically from the fluid mechanics.Many numerical methods for the solution of these balanced laws were developed in recent decades.The numerical methods are based on two views: solving hyperbolic PDE with a nonzero source term (the obvious description of the central and central-upwind schemes; (Kurganov & Levy, 2002;LeVeque, 2004)) or solving the augmented quasilinear nonconservative formulation (Gosse, 2001;Le Floch & Tzavaras, 1999;Parès, 2006).Furthermore, the methods can be interpreted using flux-difference splitting (or flux-vector splitting), or by selecting adaptive intervals and the transformation to the semidiscrete form (for example (Kurganov & Petrova, 2000)).We prefer the augmented quasilinear nonconservative formulation solved by the flux-difference splitting in our text.We try to formulate the methods in the most general form.The range of this text does not give the complete overview of currently used methods.

Mathematical models
In this section we describe the specific mathematical models based on hyperbolic balanced laws.There are many models that describe fluid flow phenomena but we are interested in the two type of them: models described open channel flow and urethra flow.

Shallow water equations
We are interested in solving the problem related to the fluid flow through the channel with the general cross-section area described by a t + q x = 0, (1) where a = a(x, t) is the unknown cross-section area, q = q(x, t) is the unknown discharge, b = b(x) is given function of elevation of the bottom, g is the gravitational constant and here η is the depth integration variable, h(x, t) is the water depth and σ(x, η) is the width of the cross-section at the depth η.The derivation can be found in e.g.(Cunge at al., 1980).
The first special case are the equations reflecting the fluid flow through the spatially varying rectangular channel a t + q x = 0, (4) with l = l(x) being the function describing the width of the channel, and second one the system for the constant rectangular channel In the above equation, h(x, t) is the water depth and u(x, t) is the horizontal velocity.It also possible to add some friction term to the system described above.For example, we can write where M is Manning's coefficient.
All of the presented systems can be written in the compact matrix form q t +[f(q, x)] x = ψ(q, x), with q(x, t) being the vector of conserved quantities, f(q, x) the flux function and ψ(q, x) the source term.This relation represents the balance laws.
It is possible to use any augmented formulation, which is suitable for rewriting the system to the quasilinear homogeneous one.For example, we can obtain Numerical Schemes for Hyperbolic Balance Laws.Applications to Fluid Flow Problems 3

Urethra flow
We now briefly introduce a problem describing fluid flow through the elastic tube represented by hyperbolic partial differential equations with the source term.In the case of the male urethra, the system based on model in (Stergiopulos at al., 1993) has the following form where a = a(x, t) is the unknown cross-section area, q = q(x, t) is the unknown flow rate (we also denote v = v(x, t) as the fluid velocity, v = q a ), ρ is the fluid density, a 0 = a 0 (x) is the cross-section of the tube under no pressure, β = β(x, t) is the coefficient describing tube compliance and λ(Re) is the Mooney-Darcy friction factor (λ(Re)=64/Re for laminar flow).Re is the Reynolds number defined by where μ is fluid viscosity.This model contains constitutive relation between the pressure and the cross section of the tube where p e is surrounding pressure.

Conservative problems and numerical schemes
We consider the conservation law in the conservative form q(x,0)=q 0 (x), x ∈ R, The numerical scheme based on finite volume discretization in the conservation form can be written as follows, We use also the semidiscrete version of ( 14) The relation ( 14) can be derived as the approximation of the integral conservation law at the interval x j−1/2 , x j+1/2 from time level t n to t n+1

37
Numerical Schemes for Hyperbolic Balance Laws -Applications to Fluid Flow Problems

www.intechopen.com
The previous relations lead to the following approximations of integral averages Numerical fluxes F n j+1/2 are usually defined by the approximate solution of the Riemann problem between states Q n j+1 and Q n j (this technique is called the flux difference splitting; it will be described in the following parts) or by the Boltzmann approach (flux vector splitting; it will be described later).In what follows, we use the notation Q + j+1/2 and Q − j+1/2 for the reconstructed values of unknown function.Reconstructed values represent the approximations of limit values at the points x j+1/2 .The most common reconstructions are based on the minmod function (see for example (Kurganov & Tadmor, 2000)) or ENO and WENO techniques (C ˇrnjaric ˇ-Zic ˇat al., 2004).
If the exact solution of the problem has a compact support in the interval 0, T then it is possible to show, that the scheme (15) is conservative, i.e.
The uniqueness of discontinuous solutions to the conservation laws is not guaranteed.Therefore the additional conditions, based on physical considerations, are required to isolate the physically relevant solution.The most common condition is called entropy condition.The unique entropy satisfying weak solution q holds for the convex entropy functions η(q) and corresponding entropy fluxes ϕ(q) (for example see (LeVeque, 2004)).
If the conservation law ( 13) is solved by the consistent method in conservation form (15) the Lax-Wendroff theorem is valid (LeVeque, 2004): if the approximate function Q(x, t) converged to the function q(x, t) for the ∆x, ∆t → 0, the function q(x, t) is the weak solution of the problem (13).
Many theoretical results in the field of hyperbolic PDEs can be found in the literature.For example, for scalar problems with the convex flux function, the convergence of the some method to the entropy-satisfying weak solutions, is proven.It means that if the solution q ǫ of the problem exists and the limit q * = lim ǫ→0 + q ǫ exists too than q * is the entropy-satisfying weak solutions of the problem (20).
It is possible to define the appropriate properties of the methods.For example, in the case of scalar problem, TVD (Total Variation Diminishing) property, which ensures that the total 38 Finite Volume Method -Powerful Means of Engineering Design www.intechopen.comvariation of the solution is non-increasing, i.e.
for all time layers t n .This property is important for limitation of the oscillations in the solution.

Nonconservative problems
We consider the nonlinear hyperbolic problem in nonconservative form q(x,0)=q 0 (x), x ∈ R.
The numerical schemes for solving problems ( 22) can be written in fluctuation form where ) are so called fluctuations.They can be defined by the sum of waves moving to the right or to the left.The directions are dependent on the signs of the speeds of these waves, which are related to the eigenvalues of matrix A(q).When the problem ( 22) is derived from the conservation form (13), i.e. f ′ (q)=A(q) is the Jacobi matrix of the system, fluctuations can be defined as follows (24)

Riemann problem for conservative systems
The Riemann problem is the special problem based on finite volume discretization with the discontinuous initial condition.In the nonlinear case it has the form We solve the transitions between two states, but this solution could not exists in general.These transitions can be rarefaction waves, shock waves or contact discontinuities.Rarefaction wave is case of a continuous solution when the following equality holds where q′ (ξ)=α(ξ)r p (ξ) (27)

39
Numerical Schemes for Hyperbolic Balance Laws -Applications to Fluid Flow Problems www.intechopen.comfor any function ξ(x, t), where α(ξ) is a coefficient dependent on the function ξ and R p (ξ) is the corresponding p-th eigenvector of the Jacobi matrix f ′ (u).
Shock waves and contact discontinuities are special cases of discontinuous solutions.The requirement that this solution should be (see (LeVeque, 2004)) a weak solution of the problem (25) leads to the following relation where s is the speed of the propagation of the discontinuities.The relation ( 28) is known as the Rankine-Hugoniot jump condition.In some cases it is possible to construct the solution of Riemann problem as the sequence of the transitions between the discontinuous states (in the case of strong nonlinearity the discontinuous solution consists of the shock waves, in the case of linear degeneration solution consists of the contact discontinuities).
In the case of a linear problem f(q)=Aq it is known that the initial state of each characteristic variable w p is moving at a speed that corresponds to the eigenvalue λ p of the matrix A.
The solution is a system of constant states separated by discontinuities that move at speeds correspondent to eigenvalues.Therefore, the jump Q j+1 − Q j over p-th discontinuity can be expressed as, (w where r p is the p-th eigenvector of matrix A. The relation ( 29) represents the initial jump in the characteristic variable w p and at the same time q = Rw.Therefore, the solution of linear Riemann problem can be defined by the decomposition of the initial jump of the unknown function to the eigenvectors r p of the Jacobi matrix A = f ′ (q) The discontinuities W p = α p r p are called waves and they are propagated by the speeds λ p .For details see (LeVeque, 2004).

Riemann problem for nonconservative systems
In this section we are interested in nonlinear systems in nonconservative form where q ∈ R m , q → A(q) is smooth locally bounded matrix-valued map, matrix A(q) is strictly hyperbolic (diagonalizable, with real and different eigenvalues).We suppose that A(q) is not Jacobi matrix so it is not possible to rewrite nonconservative system in conservative form (13). Here, BV[(R)] m is function space contains functions with bounded total variation.
Above mentioned system makes sense only if q is differentiable.In the case when q admits discontinuities at a point, A(q) may admits discontinuities as well and q x contains delta function with singularity at this point.Then A(q)q x is product of Heaviside function with delta function and in general is not unique.There is possibility to smooth out" of this function over width ǫ, for example by adding viscosity or diffusion.Then we get well defined product 40 Finite Volume Method -Powerful Means of Engineering Design www.intechopen.com of continuous functions.The limiting behavior for ǫ → 0 is strongly depend on "smoothing out." Under a special assumption, the nonconservative product A(q)q x can be understood as a Borel measure.In the following we introduce basic theorems and definitions, for details see (Gosse, 2001;Le Floch, 1989;Le Floch & Tzavaras, 1999;Parès, 2006).
Definition 4.1.A path φ in Ω ∈ R m is a family of smooth maps 0, 1 ×Ω × Ω → Ω satisfying: • φ(0; q l , q r )=q l and φ(1; q l , q r )=q r , ∀q l , q r ∈ Ω, ∀s ∈ 0, 1 , • for each bounded set O∈Ω, there exists a constant k > 0 such that ∂φ ∂s (s; q l , q r ) ≤ k |q r − q l | for any q l , q r ∈Oand almost all s ∈ 0, 1 , • for each bounded set O∈Ω, there exists a constant K > 0 such that ∂φ ∂s (s; for any q 1 l , q 1 r , q 2 l , q 2 r ∈Oand almost all s ∈ 0, 1 .Theorem 4.1 (Dal Maso, Le Floch, Murat).Let q : (a, b) → R m be a function with bounded variation and A : R m → R m×m a locally bounded function.There exists a unique signed Borel measure μ on (a, b) characterized by following properties: where we denote q(x − 0 )= lim Remark 4.1.Borel measure μ is called nonconservative product and is usually written [A(q)q x ] φ .
Remark 4.2.In the case where A(q)=f ′ (q) then Borel measure, [A(q)q x ] φ = f(q) x is independent of the path φ.
Definition 4.2 (Weak solution).Let φ be a family of paths in the sense of definition 4.
] m is a weak solution of system (31), if it satisfies as a bounded Borel measure on R × R + .
A weak solution is said to be entropic if it satisfies the inequality in the sense of distribution.
Function q(x, t) is weak solution if and only if, across a discontinuity with speed ζ, it satisfies generalized Rankine-Hugoniot condition, Now we define the Riemann problem for nonconservative strictly hyperbolic system: with an initial condition q(x,0)=q 0 (x)= q l pro x < 0, where q l , q r ∈ R m are vectors of constants.
Theorem 4.2.Let φ be a family of path in the sense of definition 4.1.Assume that system (36) is strictly hyperbolic with genuinely nonlinear or linearly degenerate characteristic field and the family of path φ satisfies ∂φ ∂q 1 (1; q 0 , q 0 ) − ∂φ ∂q 1 (0; q 0 , q 0 )=I, Then, for |q r − q l | small enough, the Riemann problem ( 36) and (37) has a solution q(x, t) with bounded variation, which depends only on x t and has Lax's structure.That is q(x, t) consist of m + 1 constant states separated by shock waves, rarefaction waves or contact discontinuities.
The solution of the Riemann problem for nonconservative system is related to the solution of the Riemann problem for conservative system.The only difference is in the case of shock wave, precisely in Rankine-Hugoniot condition formulation.
Before we define a class of useful numerical methods, we introduce a brief motivation, which can be in details found in (Gosse, 2001).Suppose nonconservative system q t + A(q)q x = 0, x ∈ R, t > 0 (39)

42
Finite Volume Method -Powerful Means of Engineering Design www.intechopen.comand suppose family of path in the sense of definition 4.1.If q(x, t) is piecewice regular weak solution then for given time t, the Borel measure can be written in following form where index k represents number of discontinuities in solution, x k (t) are discontinuity points in the time t > 0 and The amount of the quantity on the cell boundaries x j+1/2 can be splitted into contribution to cell I j+1 and the contribution to the cell I j .In other words, we split it into two terms A ± j+1/2 .

3.
A(q,...,q)= x j+1/2 x j−1/2 A(q)q x dx. (45) Remark 4.3.If A(q)=f ′ (q) then path-conservation scheme is consistent and conservative (that is if the nonconservative system can be written as conservative one, path-conservative method become classical conservative and consistent).Notice here that A ± are fluctuations, see 3.2.

43
Numerical Schemes for Hyperbolic Balance Laws -Applications to Fluid Flow Problems www.intechopen.com

The properties of exact solution
Consider a system in conservative form with special right hand side (this system agrees with shallow water system (5)) q t + f(q) x = S(q)σ x .
(46) This system can be rewritten in a following homogeneous nonconservative form where where ∂f ∂u (w) denotes Jacobian matrix of flux f.In order to define weak solutions of the system (47) we have to choose family of paths (see definition 4.1 and theorem 4.2).And following requirement become natural: if u =[q, σ] T is a weak solution of the nonconservative system (47) and σ is constant, then q must be weak solution of the homogeneous system of conservation laws, i.e. ( 46) with zero right hand side.By this requirement we impose addition condition of family of paths: if u l and u r are such that σ l = σ r = σ, then For these systems, there are no difficulties with convergence.Moreover, the shock waves propagating in regions where σ is continuous are correctly captured independently of the choice of path.For details see (Le Floch at al., 2008).

Approximate Riemann solvers
There are many numerical schemes for solving (7) with different properties and possibilities of failing.The main types of the finite volume schemes are the central, upwind and central-upwind schemes.All these schemes relate together in variety of ways.We can interpret them as schemes with different deep knowledge about structure of the solution of Riemann problem.For example, the same scheme can be interpreted as the HLL solver (an approximate Riemann solver) or as the central-upwind scheme.
The main requirements on the numerical schemes are the consistency (in the finite volume sense, i.e. consistency with the flux function), the conservativity (if there is possibility to rewrite the problem to the conservative form it is required to have conservative numerical scheme), positive semidefiniteness (the scheme preserves nonnegativity of some quantities, which are essentially nonnegative from their physical fundamental) and the well-balancing (the schemes maintain some or all steady states which can occur).The next properties are the order of the schemes and stability.From the computational point of view there are other properties such as robustness, simplicity and computational efficiency.The second and third of these properties are typical for the central methods, but they are not common for scheme based on Roe's solver.
The steady states mean that the unknown quantities do not change in the time, i.e. q t = 0 in (7), and the flux function must balance the right hand side, i.e. [f(q)] x = ψ(q, x).Some schemes are constructed to preserve some special steady states for example called "lake at rest" in open channel flow problems, where there is no motion and the free surface is constant.

44
Finite Volume Method -Powerful Means of Engineering Design www.intechopen.com

Steady states
Roe's solver is based on upwind technique.For preserving steady states, it is possible to upwind the source terms too.See (Bermudez & Vasquez, 1994) for details.It is based on approximate Jacobi matrix constructed by the following relation (for simplicity the lower indexes are neglected) where λ p and r p are the eigenvalues and eigenvectors of matrix A. Since the matrix A is diagonalizable we have where Λ is the diagonal matrix of the eigenvalues of A and the columns of matrix R are the corresponding eigenvectors of A. Then we can derive the decomposition of the source term where Λ + and Λ − are the positive and negative parts of Λ so that Λ + + Λ − = Λ.The jump of unknown function is decomposed to the eigenvectors of matrix A Because and we get The decomposition ( 89) is used for construction of fluctuation by the same way as in the following part 5.4 in the case of flux-difference splitting method.It can be shown (Bermudez & Vasquez, 1994), that such approximation preserves steady state "lake at rest" for the appropriate approximation of the source terms.

Positive semidefiniteness
Roe's solver is not positive semidefinite in general.The main reason is in using the linearization to construct the approximate Jacobi matrix.

50
Finite Volume Method -Powerful Means of Engineering Design www.intechopen.comone.In the case of urethra flow we obtain the system of four equations, where the augmented vector of unknown functions is w =[ a, q, a 0 β , β] T .Furthermore we formally augment this system by adding components of the flux function f(u) to the vector of the unknown functions.We multiply balance law (7) by Jacobian matrix f ′ (u) and obtain following relation Because of f ′ (q)q t =[f(q)] t we obtain hyperbolic system for the flux function In the case of the urethra fluid flow modelling we add only one equation for the second component of the flux function i.e. φ = av 2 + a 2 2ρβ (the first component q is unknown function of the original balance law), which has the form Finally augmented system can be written in the nonconservative form briefly w t + B(w)w x = 0, where matrix B(w) has following eigenvalues and corresponding eigenvectors We have five linearly independent eigenvectors.The approximation is chosen to be able to prove the consistency and provide the stability of the algorithm.In some special cases this scheme is conservative and we can guarantee the positive semidefiniteness, but only under the additional assumptions (see (Brandner at al., 2009)).
The fluctuations are then defined by where ) is a suitable approximation of the source term and r p j+1/2 are suitable approximations of the eigenvectors (98).

Steady states
The steady state for the augmented system means B(w)w x = 0, therefore w x is a linear combination of the eigenvectors corresponding to the zero eigenvalues.The discrete form of the vector ∆w corresponds to the certain approximation of these eigenvectors.It can be shown (Brandner at al., 2009) that (100) Therefore we use vectors on the RHS of (100) as approximations of the fourth and fifth eigenvectors of the matrix B(w) to preserves general steady state.

Positive semidefiniteness
Positive semidefiniteness of this scheme is shown in (George, 2008) for the case of shallow water equation.It is based on a special choice of approximations of the eigenvectors (98).This, in the case of urethra flow, is more complicated because of structure of the eigenvectors.Some necessary conditions for approximation of these eigenvectors are presented in (Brandner at al., 2009).

The other methods
Many other methods exist that are suitable for solving nonhomogeneous hyperbolic PDEs.These methods are often derived from the ideas described above.Some approaches are very different.However, we describe at least briefly the two of them.

ADER schemes
ADER scheme is an approach for constructing conservative nonlinear finite volume type methods of arbitrary accuracy in space and time.The first step in ADER algorithm is the reconstruction of point-wise values of solution from cell averages at time t n via high-order polynomials.To design non-oscillatory schemes we can use for example WENO reconstruction.After this step we solve High-order Riemann problem (Derivative or Generalized Riemann problem).This problem is defined as a clasical Riemann problem with polynomial-wise initial condition (polynomials arise in reconstruction step).The solution of High-order Riemann problem is used for numerical fluxes evaluation.For details see for example (Toro & Titarev, 2006) and (Toro & Titarev, 2002).

Flux-vector splitting
This approach is based on solving Riemann problem for the augmented formulation of the system (for example (8)) and nonconservative reformulation of the zero-order terms of the 53 Numerical Schemes for Hyperbolic Balance Laws -Applications to Fluid Flow Problems www.intechopen.comright-hand-side of the equations in the form (46).This solution is decomposed to the stationary contact discontinuity described by the suitable Rankine-Hogoniot condition and the remaining part of the solution solved by the flux splitting for conservation systems.For details see for example (Gosse, 2001).

Numerical experiments
Now we present several numerical experiments of the urethra flow based on mathematical model ( 10) and shallow water flow based on model ( 5).The experiments are presented only for illustration of properties of some described methods.Detailed comparison of all presented methods can be found in using literature.

Steady urethra flow
This experiment simulates the urethra flow with no viscosity.The parameter p e = 1200 and parameter β is illustrated at the figure 2. The initial condition is defined by the following p(x,0)= 3000 Pa, for x = 0 500 Pa, else, , v(x,0)= 1 m/s, for x = 0 0 m/s, else.
It is used the constant input discharge at the point x = 0 during the whole simulation.The functions described discharge and cross section of the tube are extrapolated from the boundary st the outflow (x = 0.19).At the last figures 2 and 3 we can see the difference between the methods.While the method based on decomposition of augmented system (Sec.5.6) preserve steady state (q ≡ const.) the central-upwind method (Sec.5.2) do not preserve such general steady state.

Shallow water flow
These experiments compare solution computed by the method based on decomposition of augmented system (see section 5.6) using piecewise constant reconstruction (first order method) and the piecewise linear reconstruction (second order method).The preservation of general steady state and positive semidefiniteness of the scheme is shown.
In all experiments we used 200 grid points and the following bottom topography b )

Small perturbations from the steady state
This experiment simulates propagation of the small perturbation from the steady state "lake at rest".The initial condition is given by Water height and discharge are extrapolated at the boundaries.The propagation of the perturbation is illustrated at the Fig. 4 54  This flow demonstrates also the positive semidefinitness of the scheme which insures the water depth remains non-negative.

Conclusion
In this chapter we have tried to summarize the basic approaches for solving hyperbolic equations with non-zero right hand side.We described all methods using the flux-difference or flux-vector splitting.We paid special attention to approximations that maintain steady states.
Other desired property of the proposed schemes is positive semidefiniteness.It should be noted that it is very difficult to construct schemes that meet multiple properties simultaneously.
For example, the central and central-upwind schemes are positive semidefinite, but maintain only the special steady state lake at rest.Similarly, we can explore the links between the just mentioned two properties and the order of schemes, the amount of numerical diffusion, etc.It can not be unambiguously decided which method is the best one.Before we choose the method that we have to decide which properties are essential for the concrete simulation.We hope that our overview provides to readers at least partial guidance on how to choose the most appropriate numerical approach to solve nonhomogeneous hyperbolic partial differential equations by finite volume methods.

) 52 Finite
Volume Method -Powerful Means of Engineering Design www.intechopen.com

Fig. 2 .
Fig. 2. Urethra flow simulation solved by the augmented system decomposition6.2.2 Drainage of the reservoirHere the initial condition is defined byh + b = 0.8, q = 0. (104)This simulates reservoir initially at rest, draining onto a dry bed through its boundaries.So we use open boundary conditions at the right boundary, and reflecting boundary conditions at the left boundary.More specifically it is implemented by zero discharge Q(0, t)=0 during the whole simulation, while the water height is extrapolated from the domain.At the outflow (right boundary), the boundary conditions are implemented as follows: if the flow is supercritical, H + B and Q are extrapolated from the interior of the domain, while if the flow is subcritical, the water height H is prescribed H(1, t)=10 −16 and Q is extrapolated.Solution is depicted in Fig.5.