The complete flux scheme for conservation laws in curvilinear coordinates

• A submitted manuscript is the author's version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


Introduction
Conservation laws are ubiquitous in continuum physics, they occur in disciplines like fluid mechanics, combustion theory, plasma physics, etc.These conservation laws are often of advection-diffusion-reaction type, describing the interplay between different processes such as advection or drift, diffusion or conduction and chemical reactions or ionization.Examples are the conservation equations for reacting flow [8] or plasmas [9].Sometimes, these conservation laws hold in spherical or cylindrical geometries, and in such cases it is convenient to reformulate the conservation laws in the corresponding coordinate system.In combustion theory, for example, the study of spherical and cylindrical flames is useful for finding parameters such as burning velocity or flame curvature [1].
For space discretization of these conservation laws we consider the finite volume method in combination with the complete flux scheme to approximate the fluxes at the cell interfaces.The complete flux scheme for Cartesian coordinates is introduced in [13].The purpose of this contribution is to generalize the complete flux scheme to conservation laws in spherical or cylindrical coordinates.
The development of the complete flux scheme is inspired by papers by Thiart [10,11].The basic idea of the complete flux scheme is to compute the numerical flux at a cell interface from a local (one-dimensional) boundary value problem for the entire equation, including the source term.As such, the scheme is a generalization of the exponential scheme, where the flux is determined from a local, constant coefficient, homogeneous equation [4,6].Our approach is to first derive an integral representation for the flux, and subsequently apply suitable quadrature rules to obtain the numerical flux.As a consequence, the numerical flux is the superposition of a homogeneous and inhomogeneous flux, corresponding to the advection-diffusion operator and the source term, respectively.The resulting discretization has a three-point coupling in each spatial direction, shows uniform second order convergence and virtually never generates spurious oscillations [13].The purpose of this chapter is to extend this approach to conservation laws, where the advection-diffusion operation is formulated in spherical or cylindrical coordinates.Another important issue is the extension to time-dependent problems.The key idea is then to consider the time derivative as a source term, and to include it in the inhomogeneous flux.The resulting implicit ODE system often has small dissipation and dispersion errors [15].
We have organised our paper as follows.The finite volume method for conservation laws in spherical and cylindrical coordinates is outlined in Section 2. In Section 3 we briefly repeat the complete flux scheme for stationary, one-dimensional conservation laws in Cartesian coordinates.The extension to spherical coordinates is presented in Section 4, and the next logical extension to cylindrical coordinates, is discussed in Section 5. How to deal with time dependent conservation laws is demonstrated in Section 6 for spherical coordinates.As an example, we present in Section 7 the numerical solution of a premixed, spherical flame, and finally in Section 8, we give a summary and formulate conclusions.

Finite volume discretization
In this section we outline the finite volume method (FVM) for a generic conservation law of advection-diffusion-reaction type, defined on a domain in R d (d = 1, 2, 3).Therefore, consider the following model equation where u is a mass flux or (drift) velocity, ε ≥ ε min > 0 a diffusion coefficient, and s a source term describing, e.g., chemical reactions or ionization.The unknown ϕ is then the mass fraction of one of the constituent species in a chemically reacting flow or a plasma.The parameters ε and s are usually (complicated) functions of ϕ whereas the vector field u has to be computed from (flow) equations corresponding to (2.1).However, for the sake of discretization, we will consider these parameters as given functions of the spatial coordinates x and the time t.Moreover, in the derivation of the numerical flux, we assume that the vector field u is incompressible, i.e., ∇•u = 0.
(2.2) Equation (2.1) is a prototype of a conservation law for a mixture, defining the mass balance for ϕ, and equation (2.2) is a simplified version of the corresponding continuity equation, describing conservation of mass or charge in the mixture.
Associated with equation (2.1) we introduce the flux vector f , defined by Consequently, equation (2.1) can be concisely written as ∂ϕ/∂t + ∇•f = s.Integrating this equation over a fixed domain Ω ⊂ R d and applying Gauss' theorem we obtain the integral form of the conservation law, i.e., where n is the outward unit normal on the boundary Γ = ∂Ω.This equation is the basic conservation law, which reduces to (2.1) provided ϕ is smooth enough.
In the FVM we cover the domain with a finite number of disjunct control volumes or cells Ω j and impose the integral form (2.4) for each of these cells.The index j is an index vector for multi-dimensional problems.We restrict ourselves to uniform tensor product grids for an orthogonal, curvilinear coordinate system ξ = ξ 1 , ξ 2 , ξ 3 and adopt the vertex-centred approach [16], i.e., we first choose the grid points ξ j = ξ 1 j , ξ 2 k , ξ 3 l with j =( j, k, l), where the unknown ϕ has to be approximated and subsequently choose the control volume The boundary Γ j = ∂Ω j is then the union of six interface surfaces Γ j,j±e i (i = 1, 2, 3) where, e.g., Γ j,j+e l and perpendicular to the line segment connecting ξ j and ξ j+e 1 .The (integral) conservation law for such a control volume reads where N (j)={j ± e i | i = 1, 2, 3} is the index set of neighbouring grid points of ξ j and where Γ j,k is the face of the boundary Γ j connecting the adjacent cells Ω j and Ω k .The unit normal n on Γ j,k is directed from ξ j to ξ k .Obviously, the volume element dV and the surface elements dS have to be expressed in terms of the curvilinear coordinates ξ.Approximating the volume and surface integrals in (2.5) by the midpoint rule, we obtain the following semi-discrete conservation law for ϕ j (t) ≈ ϕ(ξ j , t), i.e., where V j is the volume of Ω j , A j,k the area of Γ j,k , φj (t) ≈ ∂ϕ/∂t(ξ j , t) and s j (t)=s(ξ j , t).Furthermore, (F •n) j,k is the normal component on Γ j,k , at the interface point ξ j,k := 1 2 ξ j + ξ k of the numerical flux vector F , approximating (f •n) ξ j,k , t .Obviously, for stationary problems the time derivatives in (2.5) and (2.6) can be discarded.
In this paper we consider the formulation of the conservation law (2.1) in terms of the spherical coordinates (r, φ, θ) and the cylindrical coordinates (r, θ, z).In the first case, we assume spherical symmetry, i.e., ϕ = ϕ(r, t) and f = f (r, t)e r .As a typical example we mention a spherical flame; see Section 7. A control volume is then given by the spherical shell Ω j =[r j−1/2 , r j+1/2 ] × [0, π] × [0, 2π) and the surface integral over Γ j = ∂Ω j can be written as where we used the shorthand notation r = r j+1/2 to denote the sphere {r j+1/2 }× 0, π × 0, 2π .Note that this expression for the surface integral of the flux is exact and replaces the second term in (2.5).For the approximation of the volume integrals in (2.5) we apply the midpoint rule, so we find where F j+1/2 (t) denotes the numerical flux approximating f (r j+1/2 , t) etc..
The FVM has to be completed with expressions for the numerical flux.We require that (F • n) j,k depends on ϕ and a modified source term s in the neighbouring grid points x j and x k , i.e., we are looking for an expression of the form where d j,k := |x j − x k |.The variable s includes the source term and an additional terms like the cross flux or time derivative, when appropriate.The derivation of expressions for the numerical flux is detailed in the next sections.

Numerical flux for Cartesian coordinates
In this section we outline the derivation of the complete flux scheme for the steady, one-dimensional conservation laws in Cartesian coordinates, which is based on the integral representation of the flux.The derivation is a summary of the theory in [3,13].The conservation law can be written as d f /dx = s with f = uϕ − ε dϕ/dx.The integral representation of the flux f j+1/2 := f (x j+1/2 ) at the cell edge x j+1/2 located between the grid points x j and x j+1 is based on the following model boundary value problem (BVP) for the variable ϕ: In accordance with (2.13), we derive an expression for the flux f j+1/2 corresponding to the solution of the inhomogeneous BVP (3.1), implying that f j+1/2 not only depends on u and ε, but also on the source term s.It is convenient to introduce the variables a, P, A and S for x ∈ (x j , x j+1 ) by Here, P and A are the Peclet function and Peclet integral, respectively, generalizing the well-known (numerical) Peclet number.Integrating the differential equation d f /dx = s from x j+1/2 to x ∈ (x j , x j+1 ) we get the integral balance f (x) − f j+1/2 = S(x).Using the definition of A in (3.2), it is clear that the flux can be rewritten as f = −ε e A d e −A ϕ /dx.Substituting this into the integral balance, isolating the derivative d e −A ϕ /dx, and integrating from x j to x j+1 we obtain the following expressions for the flux: where f h j+1/2 and f i j+1/2 are the homogeneous and inhomogeneous part, corresponding to the homogeneous and particular solution of (3.1), respectively.
In the following we assume that u and ε are constant; extension to variable coefficients is discussed in [3,13].In this case we can determine all integrals involved.Moreover,

87
The Complete Flux Scheme for Conservation Laws in Curvilinear Coordinates www.intechopen.comsubstituting the expression for S(x) in (3.3c) and changing the order of integration, we can derive an alternative expression for the inhomogeneous flux.This way we obtain see Figure 2. Note that G relates the flux to the source term and is different from the usual Green's function, which relates the solution to the source term.G is a function of the normalized coordinate σ =(x − x j )/Δx (0 ≤ σ ≤ 1) between x j and x j+1 and depends on the parameters P and σ b , the σ-coordinate of the cell boundary.Obviously, σ j+1/2 = σ x j+1/2 = 1 2 .For the constant coefficient homogeneous flux we introduce the function to denote the dependence of f h j+1/2 on the parameter values ε/Δx and P and on the function values ϕ j and ϕ j+1 ; cf.(2.13).The homogeneous flux (3.6) is the well-known exponential flux [7].
Next, we give the numerical flux F j+1/2 .For the homogeneous component F h j+1/2 we obviously take (3.6), i.e., F h j+1/2 = F h ε/Δx, P; ϕ j , ϕ j+1 .The approximation of the inhomogeneous component f i j+1/2 depends on P. For dominant diffusion (|P|≪1) the average value of G(σ; P) is small, which implies that the inhomogeneous flux is of little importance.On the contrary, for dominant advection (|P|≫1), the average value of G(σ; P)  on the half interval upwind of σ = 1 2 , i.e., the interval [0, 1 2 ] for u > 0 and [ 1 2 ,1] for u < 0, is much larger than the average value on the downwind half.This means that for dominant advection the upwind value of s is the relevant one, and therefore we replace s(x(σ)) in (3.4b) by its upwind value s u,j+1/2 , i.e., s u,j+1/2 = s j if u ≥ 0 and s u,j+1/2 = s j+1 if u < 0, and evaluate the resulting integral exactly.This way we obtain where W(z) := e z − 1 − z / z e z − 1 ; see Figure 1.From this expression it is once more clear that the inhomogeneous component is only of importance for dominant advection.We refer to (3.7) as the complete flux (CF) scheme, as opposed to the homogeneous flux (HF) scheme for which we omit the inhomogeneous component.

Numerical flux for spherical coordinates
Our objective in this section is to extend the derivation in the previous section to spherical coordinates, assuming spherical symmetry.
The stationary conservation law can be written as d r 2 f /dr = r 2 s with f = uϕ − εdϕ/dr.The expression for the flux f j+1/2 := f (r j+1/2 ) at the cell boundary r j+1/2 is based on the following model BVP for the unknown ϕ: ) where ε and s are sufficiently smooth functions of r.Moreover, we assume that u > 0 and, in view of (2.2), u satisfies the relation U := r 2 u = Const for r ∈ (r j , r j+1 ).where r 2 f h j+1/2 and r 2 f i j+1/2 are the homogeneous and inhomogeneous part of r 2 f j+1/2 , corresponding to the homogeneous and particular solution of (4.1), respectively; cf.(3.3).
We consider next the expression for the inhomogeneous flux, and first take ε(r)=Const on (r j , r j+1 ).Substituting the expression for S(r) in (4.6c) and changing the order of integration, we can derive the representation with Pj+1/2 defined in (4.9) and with G the Green's function for the flux defined in (3.5), provided the normalized coordinate σ(r) and the coordinate of the cell boundary σ j+1/2 are chosen as For arbitrary ε we can generalize (4.11) as follows where the correction factor D(r(σ))/ Dj+1/2 in (4.12a) is a consequence of the relation Dj+1/2 dr = ΔrD(r(σ))dσ.Note that a(r) > 0 implies that σ(r) defined in (4.12b) is monotonically increasing from 0 to 1 on the interval (r j , r j+1 ).Summarizing, the flux f j+1/2 is the superposition of the homogeneous and inhomogeneous flux, defined in (4.10) and (4.12), respectively.
Then, applying all approximations mentioned above, we obtain the numerical flux We refer to (4.13) as the complete flux (CF) scheme for spherical coordinates, with as special case the homogeneous flux (HF) scheme (4.13b).

Numerical flux for cylindrical coordinates
In this section we present the complete flux scheme for conservation laws in cylindrical coordinates, assuming rotational symmetry about the z-axis.Consequently, the problem does not depend on the azimuthal coordinate θ.We proceed in two steps.First, we derive the r-component of the flux in polar coordinates, so we solve an essentially one-dimensional problem, and second, we extend the scheme by including the z-component, to derive the full two-dimensional scheme.
The stationary, rotationally symmetric conservation law in polar coordinates reads d(rf)/dr = rs with f = uϕ − εdϕ/dr.We give a very concise derivation of the CF scheme, since it is quite similar to the CF scheme in spherical coordinates.To determine the integral relation for the flux f j+1/2 := f (r j+1/2 ), we consider the one-dimensional model BVP: We can essentially repeat the derivation in the previous section: integrate the conservation law from the cell boundary r j+1/2 to r ∈ (r j , r j+1 ), rewrite the flux in terms of its integrating factor, substitute the flux in the integral relation and subsequently integrate over the interval (r j , r j+1 ), to arrive at the following expressions: Next, we have to elaborate (5.3b) and (5.3c).Evaluating all integrals involved, we recover relation (4.10) for the homogeneous flux.Substituting the expression for S in (5.3c) and changing the order of integration, we obtain the expression where the normalized coordinate σ is defined in (4.12b).Finally, to derive expressions for the numerical flux, we need approximations for Dj+1/2 , a,1 , σ j+1/2 and for the integral in the right hand side of (5.4).For the latter, we replace the term rs(r) in the integrand by its upwind value rs(r) u,j+1/2 , i.e., rs(r) u,j+1/2 = r j s j if ūj+1/2 ≥ 0 and rs(r) u,j+1/2 = r j+1 s j+1 if ūj+1/2 < 0. Approximating ε by its average εj+1/2 , we obtain similar results as in Section 4, except that the harmonic average Dj+1/2 is now approximated as Dj+1/2 ≈ εj+1/2 rj+1/2 , rj+1/2 = r j+1 − r j ln r j+1 /r j .
The flux corresponding to equation (2.1) is given by We outline the derivation of the r-component of the numerical flux F r,j+1/2,l at the eastern edge of the control volume Ω j,l ; see Figure 3.The derivation of the z-component F z,j,l+1/2 of the numerical flux at the northern edge is completely analogous and is therefore omitted.The key idea is to include the cross flux term ∂ f z /∂z in the evaluation of the flux.Therefore we determine the numerical flux F r,j+1/2,l from the quasi-one-dimensional BVP: where the modified source term s r is defined by The derivation of the expression for the numerical flux is essentially the same as for (5.5), the main difference being the inclusion of the cross flux term ∂ f z /∂z in the source term.In the computation of s r we replace ∂ f z /∂z by its central difference approximation and for f z we take the homogeneous numerical flux.A similar procedure applies to the z-component of the flux, which is actually the Cartesian flux from Section 3, albeit with nonconstant coefficients, [13].
Putting everything together, we obtain the following two-dimensional complete flux scheme.

two-dimensional CF scheme
Peclet function source term with cross wind diffusion

Aspects of time integration
Next, we extend the derivation to time-dependent conservation laws, restricting ourselves to spherically symmetric conservation laws; for Cartesian coordinates see [13,15].
The semidiscrete conservation law for ϕ j (t) ≈ ϕ(r j , t) can be written as where φj (t) ≈ ∂ϕ/∂t(r j , t) and s j (t)=s(r j , t).In the following we shall omit the explicit dependence on the variable t.

Numerical example
In this section we apply the complete flux scheme to a model problem, describing a premixed spherical flame stabilized by a point source of combustible mixture.
A point source at the origin issues a mass flux of 4πU of combustible mixture.After ignition, a stable spherical flame is formed, provided the value of U is in the proper range.The governing equations for this system are given by [2, 12]: ) with β the dimensionless activation energy.In the unburnt gas mixture, far ahead of the flame front, there is no combustion product and the temperature equals the temperature of the unburnt gas.In the burnt gas, beyond the flame, we assume that the reaction is completed, and consequently the combustion product is the only species and the temperature is equal to the adiabatic temperature of the burnt gas mixture.These conditions lead to the following boundary conditions As initial conditions, we take the linear profiles C(r,0)=r/r max and T(r,0)=r/r max on the truncated domain (0, r max ) and let the solution evolve to its steady state.We take r max = 120.

97
The Complete Flux Scheme for Conservation Laws in Curvilinear Coordinates www.intechopen.comFor space discretisation of (7.1) we employ the TCF scheme (6.4a) in combination with the θ-method for time integration [5].The resulting nonlinear system at each time step is solved applying one Newton iteration step.Moreover, to enhance the robustness of the method, we bound the numerical solutions according to 0 ≤ C j , T j ≤ 1, followed by a smoothing step as follows: C j := 1 4 C j−1 + 2C j + C j+1 ), and likewise for T j .
As an example, the numerical solutions at t = 100 for U = 1.0475 × 10 4 ,Le= 1 and β = 10, 8 are shown in Figure 4, computed with grid size Δr = 0.4 and time step Δt = 0.25.The solutions exhibit a steep interior layer, the so-called flame front, connecting the (virtually) constant unburnt and burnt states.Since Le = 1, the numerical solutions for C and T are identical.The solution for β = 10 is very close to the asymptotic solution [12] C(r,0)= and likewise e T .The time histories of e C and e T corresponding to the numerical solutions in Figure 4 are shown in Figure 5.We observe a regular convergence to the steady state.Finally, in order to study the effect of preferential diffusion, the numerical simulations are repeated for Le = 0.3, and the results are shown in Figure 6.As expected, the interior layer for C is slightly wider than for T.

Conclusions and future research
In this paper we have derived complete flux schemes for spherically or cylindrically symmetric conservation laws of advection-diffusion-reaction type.An integral relation for the flux is derived from a local one-dimensional BVP for the entire equation, including the source term.Applying suitable quadrature rules, we derived expressions for the numerical flux.As a result of this procedure, we obtained a numerical flux that is the superposition of a homogeneous flux, corresponding to the advection-diffusion operator, and an inhomogeneous flux, corresponding to the reaction term.For time-dependent conservation laws, we included the time derivative in the inhomogeneous flux, resulting in an implicit ODE system.The CF-scheme has been applied to a thermo-diffusive model for a spherical flame.
Possible directions of future research are the following: first, a rigourous convergence analysis of the (stationary) CF-schemes for spherical and cylindrical coordinates, and second, a dispersion analysis of time-dependent CF scheme; cf.[15] where such analysis is presented for Cartesian coordinates.Finally, from a more fundamental point of view, it would be very interesting to base the derivation of the time-dependent CF scheme on a local initial-boundary-value problem for a truly time-dependent equation, rather than computing the flux from a quasi-stationary BVP.

99
The Complete Flux Scheme for Conservation Laws in Curvilinear Coordinates www.intechopen.com

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

3 )
flux in Cartesian coordinates, we derive an integral relation for the flux that is the superposition of the homogeneous flux, depending on the advection-diffusion operator, and the inhomogeneous flux, taking into account the effect of the source term s.Approximating all integrals involved gives us the expression for the numerical flux F j+1/2 .Analogous to (3.2) we introduce the variables D, a, P, A and S, defined by D := r 2 ε, a := We refer to P and A as the Peclet function and Peclet integral, respectively.Integrating the conservation law from r j+1/2 to r ∈ (r j , r j+1 ), we obtain the relation r 2 f (r) − r 2 f (r j+1/2 )=S(r).
Scheme for Conservation Laws in Curvilinear Coordinates www.intechopen.comUsing the definitions of D and A in (4.3), it is clear that the expression for the flux can be rewritten as r 2 f = U ϕ − D dϕ dr = −D e A d dr ϕ e −A .(4.5) Inserting this expression in (4.4), isolating the derivative d ϕ e −A /dr, and integrating the resulting equation from r j to r j+1 we obtain the following expressions for the flux: )b(r) dr. ( .1b) where ε and s are sufficiently smooth functions of r and where, because of (2.2), u satisfies the relation U := ru = Const.The definitions of the variables a, P and A in (4.3) still hold, whereas the definitions of D and S change to D := εr, S(r) := r r j+1/2 ηs(η) dη.(5.2) 92 Finite Volume Method -Powerful Means of Engineering Design www.intechopen.com .3a) rf h j+1/2 = e −A j ϕ j − e −A j+1 ϕ e −A dr, (5.3c) thus, as anticipated, the flux f j+1/2 is the superposition of the homogeneous flux f h j+1/2 and the inhomogeneous flux f i j+1/2 ; cf.(4.6).
u,e = s r,C if ūr,e ≥ 0 s r,E if ūr,e < 0 s z,u,n = s z,C if ūz,n ≥ 0 s z,N if ūz,n < 0 inhomogeneous fluxrF i r e = Δr σ e − W(P r,e ) rs r u,e F i z,n = Δz 1 2 − W( Pz,n ) s z,u,n σ e = ln r e /r C ln r E /r C complete flux rF r e = rF h r e + rF i r e F z,n = F h z,n + F i z,n .The stencil of the flux approximation for F r,e is depicted in Figure3.Assume first that ūr,e > 0.Then F r,e depends on ϕ in the grid points x C and x E ,o ns in the central point x C and on the homogeneous fluxes F h z,n and F h z,s and through these fluxes again on ϕ in x N and x S .For ūr,e < 0, F r,e again depends on ϕ C and ϕ E , but this time on the source term s E and the homogeneous fluxes F h z,En and F h z,Es , inducing a further dependency on ϕ NS and ϕ SE .Thus, in95The Complete Flux Scheme for Conservation Laws in Curvilinear Coordinates www.intechopen.comaddition to the direct neighbours, F r,e depends on a few other values of ϕ, determined by the local upwind direction.
1i f r ≥ r f , , with r f = 93.4 the radius of the flame.For decreasing β the reaction slows down, resulting in a slightly wider flame front and a location of the flame front closer to the source.We define e C := C n+1 − C n /Δt 1 /N with C n = C n j T and N the number of grid points,