A High Resolution Total Variation Diminishing Scheme for Hyperbolic Conservation Law and Related Problems M. K. Kadalbajoo Department of Mathematics, Indian Institute of Technology, Kanpur, India. email: [email protected]

Ritesh Kumar Department of Mathematics, Indian Institute of Technology, Kanpur, India. email: [email protected]

Abstract A high resolution, explicit second order accurate, Total Variation Diminishing (TVD) difference scheme using flux limiter approach for scalar hyperbolic conservation law and closely related convectively dominated diffusion problem is presented. Bounds and TVD region is given for these limiters such that the resulting scheme is TVD. A class of such limiters is given which gives second order accurate TVD scheme.Numerical results for one dimensional scalar conservation law and convectively dominated diffusion problem are presented. Comparison of numerical results with some of the classical difference schemes are given. Key words: Conservation Law, Convection, Finite Difference Schemes, High Resolution Schemes, Total Variation Diminishing Schemes.

1

Introduction

For simplicity we consider the scalar conservation law ∂u ∂f (u) + =0 ∂t ∂x ∂u ∂u + a(u) =0 ∂t ∂x

(1)

and closely related one dimensional convectively dominated diffusion equation Preprint submitted to Elsevier Science

29 November 2006

∂u ∂f (u) ∂ ∂u + = (Γ ) ∂t ∂x ∂x ∂x ∂u ∂ ∂u ∂u + a(u) = (Γ ) ∂t ∂x ∂x ∂x

(2)

together with appropriate initial and boundary conditions. Here a(u) = ∂f∂u(u) is the characteristic (convective) speed depending on the variable u, Γ may be a turbulent diffusion coefficient and u may be the concentration of a passive tracer, the temperature, salinity or the velocity component. It is easy to see that depending on the values of a and Γ, Eq.(2) changes its character. For a = 0 and Γ 6= 0 it is parabolic; it is hyperbolic when a 6= 0 and Γ = 0. If a, Γ are functions of x and t, the character changes locally with time. The numerical treatment of such problems involves specific difficulties, which mainly originate from the different scales of convective and turbulent motion. It has been seen that for parabolic (pure diffusion) problems, Euler forward in time and central in space difference scheme is suitable, but for convectively dominated problems, the problems associated with hyperbolic conservation laws comes into the play which makes difference schemes for convective terms quite sensitive to stability and accuracy. Although the classical first order finite difference methods are monotonic and stable, but they suffer from numerical dissipation, and made the solution of the conservation law to become smeared out and are often grossly inaccurate. On the other hand, classical higher order finite difference methods (e. g. Lax-Wendroff, Warming-Beam, QUICK ([1]),([2])) are less dissipative in nature but are susceptible to numerical instabilities and exhibit spurious oscillations around the point of discontinuity, i. e. at turning points.Infact it happens that turning points and steep gradients tend to appear together, and more repidly the gradient changes, the more exaggerated the oscillations. This is why some authers write this phenomena occurs in the region of large gradients.([3]). For more details one can look into ([3],[4]). Apart from the numerical instabilities, fundamental mechanical or thermo-dynamical principles can be violated. Fronts will not be resolved due to numerical diffusion. It has been obsevered that, even if the initial condition for conservation law is very smooth, the wiggles get introduced into the solution by the difference scheme. One can avoid this by requiring that the total variation of the solution does not increase from time step to time step and try to find a scheme that satisfy this requirement. In [5] it has been shown that non-conservative schemes converge to wrong solution. In fact using conservative schemes to model physical problems help to avoid violation of findamental mechanical or thermo-dynamical principles. So an obvious choice for numerical treatement 2

of conservation laws is to use a TVD, conservative difference scheme. In recent years, efforts has been made into obtaining a class of schemes for conservation laws which are conservative, generally (i.e. at most places) be second order accurate (or higher) on smooth section of the solution and non oscillatory in nature and capable of resolving discontinuities in the solution. This class of schemes are known as high resolution schemes. In this paper we are proposing a high resolution TVD scheme for scalar linear conservation law and convectively dominated diffusion equation. In section 2, we are giving some useful existing results for system of hyperbolic conservation laws. we have used them to prove some results in our work as they are true for scalar case too. In section 3, we have drive the numerical scheme and proved some results for our scheme. We have given some numerical results in 4.

2

Conservative difference scheme

Integrating Eq.(2) over the rectangle [xj− 1 , xj+ 1 ] × [tn , tn+1 ] and introducing 2 2 the definitions of the spatial and temporal cell averages xk+ 1

1 unj = ∆x x hk+ 1 2

Dk+ 1 2

Z

2

u(x, tn )dx,

k− 1 2

1 1 = h(U; k + ) = 2 ∆t 1 1 = D(U; k + ) = 2 ∆t

n+1 tZ

f (u(xk+ 1 , t))dt, 2

(3)

tn

tZn+1

(Γ

tn

∂u )(u(xk+ 1 , t))dt 2 ∂x

one can obtain a difference scheme of the form un+1 = unk − k

o o ∆t n ∆t n hk+ 1 − hk− 1 + Dk+ 1 − Dk− 1 2 2 2 2 ∆x ∆x

(4)

where the subscript k and superscript n denote the k th space step and nth time step, hk± 1 and Dk± 1 are convection fluxes and diffusion fluxes on the 2 2 cell boundaries at xk± 1 respectively. 2

Difference schemes of the form given by eq.(3),and eq. (4) are called the conservative (flux form) difference schemes. The main advantage of using conservative form to model the transport of the physical variable is that with conservative form it is simpler to assure that the total physical variable is 3

conserved. The other advantage is that it becomes easier to avoid non-linear instabilities. It has been pointed out [3] that for most flows of practical situations, the classical, second order central scheme for diffusion terms is entirely adequate. Therefore, without any loss of generality, In this work we are considering only the convection term of the Eq. (1) and (2). 2.1 Some results and definitions [11] Consider hyperbolic conservation law given by eq. (1) ∂u ∂f (u) + =0 ∂t ∂x ∂u ∂u + a(u) =0 ∂t ∂x

(5)

where u is a K-vector and f(u) is a function that maps RK into RK . The corresponding conservative difference scheme is given by Ukn+1 = Ukn −

o ∆t n hk+ 1 − hk− 1 2 2 ∆x

(6)

where ∆t hk± 1 approximate the flux of material across the sides x = xk± 1 . 2

2

The approximate fluxes are written as hnk+ 1 = h(unk−p , .........., unk+q ) 2

hnk− 1 = h(unk−p−1, .........., unk+q−1)

(7)

2

where numerical flux fuction hk± 1 depend on u at p + q + 1 points. 2

Proposition 1 If h(u, u, .....u) = f (u) then difference scheme (6) with numerical flux function Eq.(7) will be consistent with conservation law Eq.(5) (cf. [11], page 157). Proposition 2 The truncation error of a conservative scheme of the form un+1 = Q(unk−(p+1) , ....... , unk+q ) k

(8)

is given by τkn = −∆t[q(u)ux ]x u=un + O(∆t2 + ∆x2 ), k

4

(9)

where u = u(x, t) is a solution to conservation law eq.(5) and 

1 1 q(u) =  2 2 R

q X

j=−(p+1)

′



j 2 Qj (u, ....... , u) − [f (u)]2  .

(10)

where R = ∆t/∆x and Qj represents the partial derivative of Q with respect to the (p + 2 + j)-th variable. (cf. [11], page 184). Definition 3 A difference scheme is said to be total variation diminishing scheme (TVD) if the numerical solution obtained by the scheme satisfies ∞ X

|δ+ un+1 k | ≤

k=−∞

3

∞ X

|δ+ unk |

for all n ≥ 0

(11)

k=−∞

High resolution scheme

In this section, we are deriving a higher resolution TVD scheme using flux limiters in much the same way as described by Sweby [6]. The approach he followed was similar to the Flux Corrected Transport (FCT) of Boris and Book [7], the application of a low order scheme along with an additional “limited” (or Corrected) flux. This flux is a difference between the flux of a high order scheme and that of the low order scheme which has been “corrected” in such a way that the resulting scheme is TVD. There are two main differences between the approach adopted by Sweby and of Boris and Book [7]. Firstly the FCT algorithm was essentially a two-step procedure, whereas Sweby [6] adopted a single-step procedure; and secondly the FCT limiter was constricted by unity whilst the Sweby allowed a more generous upper limit. We are giving a single-step procedure for the hyperbolic scalar conservation law similar to Sweby [6]. Though the approach is similar but the our proposed TVD scheme, TVD region and the class of limiters do not match with Sweby’s work. Infact first order upstream and second order upstream scheme, both obey the hyperbolicity condition. Hence using their linear combination seems to give better accuracy. The numerical results of proposed scheme also support this fact. For clarity of approach we consider the linear scalar equation ut + aux = 0,

a>0

(12) 5

The approach that we are using to construct high resolution scheme as follows. We define the numerical flux function of the new scheme in terms of numerical flux function of the first order upstream scheme and second order upstream scheme as hnk+ 1 = hnL

1 k+ 2

2

+ φnk [hnH

1 k+ 2

− hnL

1 k+ 2

]

(13)

where hL and hH are the numerical flux function of first order and second order upstream schemes respectively, and are given by hnL

k+ 1 2

hnH

k+ 1 2

= aunk

(14)

a = (3unk − unk−1) 2

(15)

and φ is the limiter which is defined in such a way that the resulting numerical flux function has the smoothing capabilities of the lower order scheme when it is necessary (i.e near discontinuities) and the accuracy of the higher order scheme (i.e. in the smooth sections of solution) when it is possible. Note that, since we formulate the scheme by defining the numerical flux function, the resulting scheme will automatically be conservative. Since the flux given by hL is considered to be too diffusive, the correction term φnk [hnH 1 − hnL 1 ] is k+

2

k+

2

referred to as the anti-diffusive flux. We try to add as much as anti-diffusive flux as possible without increasing the variation of the solution. The resulting conservative difference scheme is un+1 k

=

unk

−c



(unk

−

unk−1)

o 1n n n φk (uk − unk−1) − φnk−1(unk−1 − unk−2) (16) + 2 

△t where c = a △x is the courant number, and φ is taken to be nonnegative so as to maintain the anti-diffusive flux.

Like Roe [8], Vanleer [9] and Sweby [6] and many others we take the limiter φ to be a function of consecutive gradients (in the case of fixed mesh size it turns out to be function of consecutive differences), i.e. φnk = φ(θkn ). The smoothness parameter θ is given by θkn

=

∆unk− 1

(17)

2

∆unk+ 1

2

where ∆− uk+1 = ∆+ uk = ∆uk+ 1 = uk+1 − uk . 2

6

We now seek to choose the function φ(θ) in such a way that the limited anti-diffusive flux is maximized in amplitude subject to the constraint of the resulting scheme being TVD. Theorem 4 The scheme given by Eq.(16 ) is TVD under the constraint 0 ≤ φ ≤ 2, 0 ≤ φ ≤ 2θ . Proof Let us rewrite the given scheme as un+1 k

=

unk

φn − c 1+ k 2 "

!

∆− unk−1 c − φnk−1 ∆− unk n 2 ∆− uk #

(18)

Comparing scheme with the incremental form (I-form) uk+1 = unk + Ck+ 1 ∆+ unk − Dk− 1 ∆− unk k 2

(19)

2

where Ck+ 1 and Dk− 1 depends on unk and its neighboring values, so we have 2

2

Ck+ 1 = 0, Dk− 1 2

2

φn = c 1+ k 2 "

!

c ∆− unk−1 − φnk−1 2 ∆− unk

#

(20)

Harten [10] has proved that a sufficient condition for any scheme of the form Eq.(19) for solving Eq.(1) to be TVD is that the coefficients Ck+ 1 , Dk+ 1 2 2 satisfy Ck+ 1 ≥ 0; Dk+ 1 ≥ 0; 0 ≤ Ck+ 1 + Dk+ 1 ≤ 1 2

2

2

2

(21)

Hence scheme (18) will be TVD if Ck+ 1 and Dk+ 1 given by Eq.(20), satisfy 2 2 Conditions (21), i.e. if 0 ≤ Dk− 1 ≤ 1 2

i.e., φn 0≤ c 1+ k 2 "

!

∆− unk−1 c − φnk−1 ≤1 2 ∆− unk #

(22)

using the definition of smootheness parameter θ, we see that Inequality (22) 7

satisfies if (dropping out the suffixes and prefixes) −2 ≤ φ − φθ ≤ 2

(1 − c) c

(23)

under the CFL condition 0 ≤ c ≤ 12 , (which is the stability condition for second order conservative scheme using numerical flux given by Eq.(15) ). We see that Inequality Eq.(23) satisfies for 0 ≤ |φ − φθ| ≤ 2

(24)

In addition to requiring φ(θ) to be nonnegative we must have φ(θ) = 0, θ < 0. It is easily seen that Inequality (24) along with non-negativity of φ are satisfied for the bounds 0 ≤ (φ, φθ) ≤ 2

(25)

Hence for the scheme Eq.(16) to be TVD the limter function must satisfy bounds given by (25). Geomatrically limiter function φ(θ) must lie in the shaded region of Fig. 1(a). Proposed TVD region is compared with TVD region proposed by Sweby in fig 1. 3

3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

1

2

3

4

5

6

7

8

Proposed TVD region 1(a)

0

0

1

2

3

4

5

Sweby’s TVD region 1(b). Fig. 1.

We need to maximize the limiter φ in order to maximize the antidiffusive flux that we add to the first order scheme, subject to the TVD constraints. An obvious choice is φ(θ) = min(2, 2/θ),

θ>0

which is the upper boundary of the region shown in Fig. 1. The extra condition imposed on φ in order to make resulting scheme Eq.(16) be second order accurate whenever posssible is that φ(1) = 1. This is also required for Lipschitz continuity of φ(θ). Note that the measure of smoothness θ breaks down near extreme points of solution, where the denominator may be close to 8

zero and θ become arbitrary large, or negative even if the solution is smooth. Hence maintaining the second order accuracy at extreme points is impossible with TVD methods. Following theorem gives conditions on φ which garantees second order accuracy away from extrema of solution. Theorem 5 If the flux limiter φ is bounded, then the difference scheme (16) is consistent with the linear scalar equation (12). If φ(1) = 1 and φ is Lipschitz continuous at θ = 1, then the difference scheme (16) is second order accurate on smooth solutions with ux bounded away from zero. Proof Consider the associated numerical flux function with difference scheme (16) i a h hk+ 1 = h(unk+1 , unk , unk−1) = aunk − φnk unk − unk−1 2 2

Obviously, since φ is bounded, for linear scalar equation (12) h(u, u, u) = f (u) First claim follows from Proposition 1. We prove a slight variation of the second claim. So that we can apply Proposition 2, we prove the second claim under the assumption that φ is differentiable at θ = 1. The proof using the assumption that φ is only Lipschitz continous at θ = 1 is a slight, technical variation of the proof given. The difference scheme (16) can be written as, un+1 = Q(unk−2, ....., unk+1) k = unk

−c



(unk

−

unk−1)

o 1n n n + φk (uk − unk−1) − φnk−1(unk−1 − unk−2 ) (26) 2 

where p = q = 1, and note that φnk , φnk−1depend on unk+1, unk , unk−1, unk , unk−1, unk−2 respectively. An easy calculation shows that 1 X

j 2 Qj (u, u, u, u) = 0

(27)

j=−2

by Eq. (10) q(u) =

−a2 2

(28)

9

hence by Eq. (9) τkn

= −∆t

(

a2 − ux 2

! ) x

+ O(∆t2 + ∆x2 ),

(29)

u=un k

Eq.(29 shows that truncation error is of O(∆t, ∆x2 ), i.e. the difference scheme is atleast second order accurate. In the light of above results we are proposing a new class of limiters which are Lipschitz continuous and lie within the proposed TVD region such that φ(1) = 1. This class of limiters is given by

φ(θ) =

  

  min

0, h

θ≤0

i

1+δ , θ > 0, for any fixed δ ∈ [0, ∞) 2, θ2 , δ+θ

We see that for every value θ > 0 and any fixed δ ∈ [0, ∞), the limiters are bounded, Lipschitz continuous and satisfies the condition φ(1) = 1, hence gives second order accuracy in the smooth region of solution. Also since this class of limiters lies in the proposed TVD region hence yields TVD scheme.

4

Numerical Results

In all the Numerical results we have used the limiters corresponding to δ = 1 and δ = 9 for proposed scheme. We have found camparable results for limiters corresponding to some of the other test values of δ ∈ [0, ∞).

4.1 Pure Convection problem with discontinous Initial Condition [3]

One of the simplest model problems is the one-dimensional convection at constant velocity of an initial data with step function in the variable u without physical diffusion. ut + aux = 0,

(30) 10

The analytical solution of the above equation is u = f (x − at), where f (x) is the initial distribution of u given by    1,

f (x) = 

 0,

x ≤ 0.1

(31)

x > 0.1

The solution describes a wave propagating in the positive x-direction (if a > 0) with velocity a. Comparison of proposed scheme with some of the classical schemes and existing high resolution schemes are presented. Numerical results for various scheme are shown in fig. 2 and fig. 3 for dimensionless convection velocity a = 0.01, grid size ∆x = 0.005 for the domain x ∈ [0, 1], Courant number c = 0.01, and c = 0.5 at dimensionless time t = 50, at this time the jump is moving through x = 0.6 from its initial position x = 0.1. We are using the Courant number c = 0.01 for numerical comparison because our main interest is in various spatial differnce schemes for convection . Since all the campared schemes are not second order accurate in time, so in order to avoid numerical error due to discretization in time as far as possible, we have chosen such a small Courant number. For Courant Number c = 0.5 comperison are given only with Sweby’s scheme in fig.4 as most of the higher order schemes become unstable and give too oscillatory results . The classical first order Upstream scheme is diffusive and solution smeared out because of the inherent numerical diffusion of these schemes. For first order schemes reducing the grid size ∆x (increasing grid number N) reduces numerical diffusion, but it makes using these schemes computationally expenssive. On the other hand Standard second or higher order difference methods, e.g., Central, Lax-Wendroff, Second order upstream, Beam-Warming, Fromm, QUICK schemes eliminate such diffusion but introduce dispersive effects that lead to unphysical oscillations in the numerical solution. The central and Lax-wendroff difference schemes have leading odd order error term which corresponds to propagating numerical dispersion terms which corrupt large region of flow with unphysical oscillations. Also, these are behind the advancing front and are damped with the distance to the front. The second order upstream and Beam-Warming schemes eliminate artificial diffusion, while minimizing numerical dispersion. Their leading truncation errors (potentially oscillatory) third order derivative terms. The damped oscillations before the advancing fronts are typical of these second order upstream difference methods. However, the fourth order derivative (corresponds to dissipative error) is large enough to dampen short-wavelength components of the dispersion to some extent. QUICK scheme has a leading fourth order derivative truncation error term which is dissipative, but higher order dispersion term still cause over11

shoots and a few oscillations when excited by nearly discontinuous behavior of the advection variable. But, they are considerably smaller than the other second order schemes see fig. 2. Numerical results for Courant number C = 0.5 are compared in Fig.4.

12

1.4

1.4 ’exsol’

1.2 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

0

0.2

0.4

0.6

0.8

’exsol.txt’

1.2

1

−0.4

0

I Order Upstream

0.6

0.8

1

1.4 ’exsol.txt’

1.2

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2 0

0.2

0.4

0.6

0.8

’exsol.txt’

1.2

1

1

−0.4

0

0.2

Fromm

0.4

0.6

0.8

1

Central

1.4

1.4 ’exsol’

1.2 1

1 0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2 0

0.2

0.4

0.6

0.8

’exsol.txt’

1.2

0.8

−0.4

0.4

II-Order Upstream

1.4

−0.4

0.2

1

−0.4

0

Lax-Wendroff

0.2

0.4

0.6

0.8

1

Beam-Warming

1.4 1.2

exsol−−−−−

1 0.8 0.6 0.4 0.2 0 −0.2 −0.4

0

0.2

0.4

0.6

0.8

1

QUICK scheme Fig. 2. Numerical Results of some of the classical schemes for travelling shock problem corresponding to Courant number C=0.01.

Numerical results of proposed scheme with limiter function corresponding to δ = 1 and δ = 9, are compared with Sweby’s [6] scheme un+1 k

=

unk

−c



∆− unk

1 + (1 − c)∆− (φnk ∆+ unk ) 2 13



(32)

using Superbee, Vanleer and BW-LW(minmod) limiters. Superbee Vanleer

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

BW-LW

|θ| + θ 1 + |θ|

(34)

φ(θ) = max [0, min(θ, 1)]

1.4

1

1 0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2 0

0.2

0.4

0.6

0.8

’exsol’

1.2

0.8

1

−0.4

0

Proposed Scheme for δ = 9

0.2

0.4

0.6

0.8

1

Proposed scheme for δ = 1

1.4

1.4 ’exsol’

1.2

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2 0

0.2

0.4

0.6

0.8

’exsol’

1.2

1

−0.4

(35)

1.4 ’exsol’

1.2

−0.4

(33)

1

−0.4

0

0.2

0.4

0.6

0.8

1

Sweby’s scheme (superbee limiter) Sweby’s scheme (vanleer limiter) 1.4

exact sol .........

1.2 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4

0

0.2

0.4

0.6

0.8

1

Sweby’s scheme (BW-LW limiter) Fig. 3. Numerical Results of the High Resolution schemes corresponding to Courant number C=0.01.

14

1.4

1.4 ’exsol’

1.2 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

0

0.2

0.4

0.6

0.8

’exsol’

1.2

1

−0.4

0

0.2

Lax-Wendroff

0.6

0.8

1

Beam-Warming 1.4

1.4 ’exsol’

1.2

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2 0

0.2

0.4

0.6

0.8

’exsol.txt’ ’exsol’

1.2

1

−0.4

0.4

1

−0.4

0

0.2

0.4

0.6

0.8

1

Sweby’s scheme (superbee limiter) Sweby’s scheme (vanleer limiter) 1.4

1.4 ’exsol’

1.2

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

0

0.2

0.4

0.6

0.8

’exsol’

1.2

1

1

−0.4

Sweby’s scheme (BW-LW)

0

0.2

0.4

0.6

0.8

1

Proposed scheme δ = 1

1.4 ’exsol’

1.2 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4

0

0.2

0.4

0.6

0.8

1

Proposed scheme δ = 9 Fig. 4. Numerical Results of some of the classical and High resolution schemes for corresponding to Courant number C=0.5.

Numerical results shows that proposed sceme gives remarkably sharp profile and better resolution near points of discontinuities, Fig. 3. and Fig. 4.

15

4.2 Pure Convection problem with smooth initial condition [14] Consider the Eq.(30) along with the following initial condition f (x) =

   sin(10πx)

0

 

0 ≥ x ≤ 0.1

(36)

0.1 ≥ x ≤ 1.0

with boundary conditions u(0, t) = 0 and u(1, t) = 0 The solution describes a sine wave propagating in the positive x-direction (if a > 0) with velocity a. The exact solution, upto t = 0.9/u; is

u (x, t)) =

      

0

0 ≤ x ≤ ut

sin[10π(x − ut)]       0

ut < x ≤ ut + 0.1

(37)

ut + 0.1 < x ≤ 1.0

The sine wave propagates with no reduction in amplitude at a speed a = 0.1. The exact solution and approximate solution with C = 0.4, ∆t = 0.1, ∆x = 0.025, and convection velocity a = 0.1 at time t = 8, are shown in fig. 5.

16

1

1 ’exact sol’

’exact sol’

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

0

0.2

0.4

0.6

0.8

1

1.2

−0.4

0

0.2

Lax-Wendroff

0.4

0.6

0.8

1

I order Upwind

1

1 ’exact sol’

’exact sol’

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

1.2

0

0.2

0.4

0.6

0.8

1

1.2

−0.4

0

0.2

0.4

0.6

0.8

1

1.2

Sweby’s scheme (superbee limiter) Sweby’s scheme (vanleer limiter) 1

1 ’exact sol

’exact sol’

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

0

0.2

0.4

0.6

0.8

1

1.2

−0.4

0

Proposed scheme for δ = 1

0.2

0.4

0.6

0.8

1

1.2

Proposed scheme for δ = 9

Fig. 5. Comparision of numerical results travelling sine wave problem corresponding to Courant number C=0.4.

4.3 Propagating Temperature Front [14] The simplest model to describe propagating temperature front is given by one dimensional linear convection dominated diffusion equation, ∂T ∂T ∂T ∂T +a = (Γ ) ∂t ∂x ∂x ∂x

(38)

At t = 0 a sharp front is located at x = 0. For subsequent time the front convect to the right with speed a and its profile loses its sharpness under the influence of thermal diffusivity Γ. Consequently, for a given time t, the larger the value of cell Reynolds number Rcell = a.∆x/Γ, the sharper the profile of the front. 17

The initial condition for this problem is shown in the Fig. 5. 3

2.5

2

1.5

1

0.5

0 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Fig. 6. Initial condition for propagating temperature front

with the following boundary conditions. T (−2, t) = 1.0, T (2, t) = 0.0

(39)

The series solution of this problem (cf.[[14], page 306]), by the separation of variables technique, is given by N 2X π(x − at) exp [−Γ(2k − 1)2 πt/L2 ] T (x, T ) = 0.5 − sin (2k − 1) π k=1 L 2k − 1

#

"

Numerical results for different schemes are given and compared with the exact solution. It has been seen that for low cell Reynolds number, most of the classical schemes give acceptable numerical solution Fig. 6. but for higher cell Reynold number they exhibit oscillatory behaviour, and give unphysical results. Numerical results have shown that for high Reynolds number proposed scheme give comparable accuracy with Sweby’s scheme as shown in Fig. 7.

18

1.2

1.2 ’exsol.txt’

’exsol.txt’

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2 −2

−0.2 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

−1.5

upstream

−1

−0.5

0

0.5

1

1.5

2

Lax-Wendroff

1.2

1.2 ’exsol.txt’

’exsol.txt’ 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2 −2

−0.2 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Sweby’s scheme (Superbee limiter)

−1.5

−1

−0.5

0

0.5

1

1.5

2

Sweby’s scheme (Vanleer limiter)

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

’exsol’

0

0

−0.2 −2

−0.2 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Proposed scheme for δ = 1

−1.5

−1

−0.5

0

0.5

1

1.5

2

Proposed scheme for δ = 9

Fig. 7. Numerical results when Rcell =3.33, ∆x = 0.1,c = 0.5, at time t = 1.0, a = 1.0, Γ = 0.03003

19

1.2

1.2 ’exsol.txt’

’exsol.txt’

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2 −2

−0.2 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

−1.5

−1

Upstream

−0.5

0

0.5

1

1.5

2

Lax-Wendroff

1.2

1.2 ’exsol.txt’

’exsol.txt’

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

−0.2 −2

−0.2 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

−1.5

Sweby’s scheme (Superbee limiter) 1.2

−1

−0.5

0

0.5

1

1.5

2

Sweby’s scheme (Vanleer limiter)

1.2 ’exsol’

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 −0.2 −2

0

−1.5

−1

−0.5

0

0.5

1

1.5

Proposed scheme for δ = 1

−0.2 2 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Proposed scheme for δ = 9

Fig. 8. Numerical results when Rcell =100, ∆x = 0.1,c = 0.5, at time t = 1.0, a = 1.0, Γ = 0.001

5

Concluding Remarks

The derivation of high resolution second order accurate scheme is investigated by means of adding anti-diffusive flux to a general first order monotone scheme. Constraints on the limiters, as function of gradient ratios, have been obtained so that the resulting scheme is TVD, and a class of limiters proposed which satisfy these constraints and yields a second order accurate TVD scheme for the smooth section of solution. Numerical results for one dimensional linear case are compared with some of the classical schemes along with the Sweby’s scheme using Superbee, Vanleer and BW-LW limiters. As the proposed scheme obeys the hyperbolicity condition hence for Hyperbolic cases it gives better accuracy. Numerical results have shown that in case of shocks and travelling sine wave, proposed scheme give better resolution near points of discontinuity. In 20

convectively dominated diffusion problems proposed scheme give comparable results with other high resolution schemes.

References [1] B. P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comp. Mathods in App. Mech. and Eng., 19 (1979), 59-98. [2] B. P. Leonard, Order of accuracy of QUICK and related convection-diffusion schemes, Applied Mathematics Modelling, 1995, 19(11), 640-653. [3] Yongqi Wang Kolumban Hutter, Comparisons of numerical methods with respect to convectively dominated problems, Int. J. Numer. Meth. Fluids, 37 (2001), 721-745. [4] Richard E. Ewing and Hong Wang, A summary of numerical methods for timedepandent advection-dominated partial differential equations, J. of Comp. and App. Math. 128(2001), 423-445. [5] Thomas Y. Hou and Philippe G. Lefloch, Why Non conservative schemes converge to wrong solutions: Error analysis, Math. Comp. Vol 62, No. 206, (1994), 497-530. [6] P. K. Sweby, High resolution schemes using flux limiters for Hyperbolic conservation Laws, SIAM J. Numer. Anal, Vol 21, No. 5, 1984. [7] J. P. Boris and D. L. Book, Flux corrected transport, I, SHASTA, A fluid transport algorithm that works, J. Comp. Phys., 11 (1973), 38-69. [8] P.L. Roe, Numerical algorithms for the linear wave equation, Royal Aircraft Establishment Technical Report 81047,1981. [9] B. Vanleer, Towards the ultimate conservative scheme, II. Monotonicity and conservation combined in a second order scheme, J. Comp. Phys., 14 (1974), 361-370. [10] A. Harten, High resolution schemes for conservation law, J. Comp. Phys., 49 (1983), 357-393. [11] J. W. Thomas, Numerical Partial Differential Equations, conservation laws and elliptic equations, Text in Applied mathematics 33, Springer, 1999. [12] R. J. Levaque, Numerical methods for conservation Laws, Lectures in Mathematics, ETH Zurich, Birkhauser, 1999. [13] E. F. Toro, Riemann Solvers and Numerical methods for Fluid Dynamics, A practical introduction, 2nd edition, Springer, 1999. [14] C. A. J. Fletcher, Computational Techniques for Fluid Dynamics, Part I, Springer Series in Computational Physics, 1988.

21

[15] Serge Piperno and Sophie Depeyre, Criteria for the design of limiters yielding efficient high resolution TVD schemes, Comp. and Fluids,vol. 27, No. 2, 1998, 183-197. [16] B. Jonathan and Randall J. Levaque, A Geometric Approach to High Resolution Schemes, Siam J. Numer. Anal., vol. 25 No. 2, 1998, 268-284. [17] C.W. Shu, An overview on high order numerical methods for convection dominated PDEs, in Hyperbolic Problems: Theory, Numerics, Applications, T.Y. Hou and E. Tadmor, editors, Springer-Verlag, Berlin, 2003, pp.79-88.

22