A truly third order finite volume scheme on the quadrilateral mesh Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 15, № 4, 2010

A truly third order finite volume scheme on the quadrilateral mesh

D. Petkovic, M. Petrovic Pristina University, Kosovska Mitrovica, Serbia e-mail: [email protected]

This article addresses construction of the finite volume scheme of the third-order accuracy using approximation of flows by some known function. The problem is solved by means of procedures of the multidimensional numerical integration.

Keywords: conservation law, finite volume method, multidimensional numerical integration.

1. Finite volume methods

The aim of this work is to obtain the solution of the conservation laws whose PDE formulation is

qt + f (q)x = 0. (1)

The methods that actually give the mentioned solution are finite volume methods. The building element of the structure of these methods is grid cell or grid volume (see Figure 1). The i-th grid cell is the basic element of the subdivision of the spatial domain and it is noted as:

Ci = > xi+\

Our task is to approximate the average value over the i-th interval at time tn:

xi+1

^ ~ ~Kx / q<ytn'x")dx = ¿| JQ(tn,x)dx, (2)

x. i Ci

where Ax =

the length of the cell [1]. By integrating the integral form of the conservation law and dividing by Ax afterwords, we update the cell average of q in one step time:

' in+i in+i

f (q(t,xi+i^dt- j f (q(t,Xi_i^dt

tn

The next task is to construct approximations for the integrals (fluxes) from (3)

1

(3)

tn+1

/ /(?(*.*;Li))dt. (4)

At J

tn

tn+1

T n Fi-1/2 n Fi+1/2

Ci-1 Ci Ci+1

Q-1 Qn QN+1

Fig. 1

The formulation above is an approximation to the average flux at the edge x

Xi_ i. ,2

The flux at the point x — Cctll be obtained respectively. Although we cannot know the exact value of Q™ (since q(xi_i,t) varies with the time along each edge of the cell), we can

at least suppose how the numerical method should look like:

At

Q

n+1

Ax

Q'I-^AF^-fz,

The methods of this form are so-called 3-point stencil, meaning that the value Q™+1 depends on three values: Q™-1, Q™, Qn+1 at the previous time level.

These kind of methods have the conservation property as next: let L = {C1, C2,..., CJ} be a set of cells. Then,

j

Ax Q

n+1 / y Q i

i=1

Ax

i= 1

JQi — At ( F'j+i — Fj_i

As we can see, every value of the flux function cancels out except at the extreme edges. In other words, over the entire domain we have exact conversation.

2

2. Local double logarithmic reconstruction in 1D

In order to derive a truly third order formula, a starting point is probably finding a "sufficiently" good approximation of the unknown function. Some recommendations were based on development of piecewise polynomial reconstructions with higher degree polynomials such as MUSCLE, PPM, ENO, WENO. Yet, the fact is that the polynomial of higher degree can produce oscillations. Considering the nature of the solution of conservation law which are sensitive to artificial oscillations, this can be interpreted as a disadvantage. So, polynomials are not always good enough as model functions. There are several reconstruction techniques based on the hyperbolas, such as LHR (Local Hyperbolic Reconstruction) and LHHR (Local Harmonic Hyperbolic Reconstruction). Since hyperbolas are monotonic functions, the order of a hyperbolic reconstruction will be reduced at the local extremes. In this case using some sort of limiters is necessary. To avoid working with limiters and yet keeping the desired order even at the points of discontinuity, Local Double Logarithmic Reconstruction was developed. Based on two logarithmic functions, this technique improves its forerunners LLR (Local Logarithmic Reconstruction) and DLR (Double Logarithmic Reconstruction) in the following properties:

— compared to LLR, the LDLR holds third order of accuracy at the local extremes as well and it is also local variation bounded;

— LDLR upgrades DLR because it exists for all input data.

2.1. Defining a LDLR Reconstruction Function

The aim is to achieve the full third order reconstruction scheme. The basic fact which provides wanted order of accuracy is the next lemma: Lemma. Let f,g G C3[a,b], h = b — a, and

f (x)dx = g(x)dx,

(f — g)'(a) = O(h2) = (f — g)'(b) then (f — g)(x) = O(h3) for all x G [a, b].

We are looking for well defined piecewise smooth function which would exist in each cell. This function is LOCAL, because it depends only on information from the central and nearest neighbouring cells. This is so from the reason that the central difference approximation is obtained right here. Therewith, it is composed of two logarithmic functions and it has a DOUBLE LOGARITHMIC form.

Let xj = ih, where i is an integer and h is a positive number, be the considered computational grid. Since we are dealing with a piecewise function defined in each cell Cj = [ i, xi+i], it is enough to analyse the cell Co only. As described, the reconstructing function in cell Co has the following form:

ro(x) « Ci + C2 log(x + C3) + C4 log(x + C5). (5)

Let

^o(x) = —

Ch

a

hi 2

X - x0 - o--1

2a

dh -ylog

h (2

x — x0 — — 7 — 1

2b

(6)

be some function defined for x G C0, where a,b,c,d are unknown parameters that have to be found. Now, define our reconstruction function as:

r0(x) = v0 + (p0(x) -Jlf (7)

Co

where vq = — [ ro(x)dx, is the cell average of unknown function in Co-h

Co

The reconstructing task is to determine parameters a, b, C, d according to some conditions and properties of the low variation bounded condition (for more information see [2, 3]). After deriving of needed parameters, we are able to write down a LDLR scheme. This scheme is not given in the direct form namely because, the reconstruction function ro is not used directly. Instead, the lateral derivatives are much more suitable for direct computation:

r0 ^'Xq ± —^ = v0 + ch^^a) + dhu^ib),

li+{t) = - log-

t2

(¿-l)log(l-t)-t

¡1 (t) =---. (8)

b

b

The functions ^+ and ^ are related in the next way:

¡jj^ ( b =- ) = (a — 1 )/i±(a).

V a — V

Knowing this relation the computation can be realized more easily.

3. Second order approximations to first derivatives on quadrilateral cell

When the problem we are dealing with is posed on the irregular quadrilateral mesh then we are actually dealing with quadrilateral cells. Each of these cells is unique has generally different side lengths. In short, we are dealing with specific local geometry [4]. With the given known components our problem can be posed as follows:

— function f e C3(tt), tt c R;

1 f

— /q = —(n / fix)y)dxdy, are given cell averages (i G Z);

i vol(Ci) J

Ci

— the goal: find the second order accurate approximations to the lateral derivatives of function f in a given point M' = (x',y').

The first step in realization of stated intention is to expand f in a certain point M = (x, y) e tt, which lies in the neighbourhood of the known point M'. By Taylor we get:

f{x,y) = f(M') + fx(M')(x - x') + fy{M'){y - y') + l-fxx{M'){x - x')2+

+fxy(M')(x - x')(y - y') + l-fyy{M'){y - y'f + 0(\\M - M'||3). (9)

After integrating the derived expression we are getting formulas for the averages:

JCi = f(M') + ^^ JJ(x - x')dxdy + JJ(y- y')dxdy+

Ci Ci

Ci Ci

i Ci

where Ai = diam(Ci) = max ||p — q||. Conventionally, we can rewrite the previous state-

p,q&Ci

ment as:

fa = f(M') + aJx(M') + A/»(M') + \%fxx(M') + 6tfxy(M') + ^/„(M') + 0(A3), (11)

where

a* = ^cJI{x~x')dxdy y')dxdy>

Ci Ci

Yi

volCj

(x — x')2 dxdy Si

volCi

(x — x')(y — y')dxdy,

Ci

Ci

£i =

volCi

(y — y')2dxdy.

(12)

Ci

The desired approximation of the gradient of f is described by next expression:

/am') ■

(13)

where j = 0, 5 is a local index. Therewith, the following conditions have to hold for Kj, j = 0, 5:

= 0, constant terms ajKj = 1, ^^ftjKj = 0, j j j

first order integral terms

Yl Yj Kj = Y15j Kj = Y1 £j Kj = ^ (14)

second order integral terms. These expressions form a linear system of six equations with six unknown variables Kj, j = 0,5 which matrix equation is

where

A - k = 12

1 1 11 1 1 " 0 " ko

ao a\ a2 a3 a4 a5 1 ki

A= Po Pi P2 P3 P4 P5 , 12 = 0 , k = k2

Yo Yi Y2 Y3 Y4 Y5 0 k3

So Si S2 S3 S4 S5 0 k4

. £o £1 £2 £3 £4 £5 _ 0 k5

Analogously, the formula for the y-derivatives states:

J2pjJct = fy(M') + j2pMA33)

(15)

(16)

(17)

and the system which determines pj coefficients is:

Ap =13 = [ 0 0 1 0 0 0 ] .

Let' concentrate on the x-derivative only (the y-derivative, of course, can be obtained analogously).

3.1. Deriving the Coefficients a, 7, 5 and e

As analysed previously, the solution of the system A • k = 12 will give us Kj. Afterwards we are able to derive the approximation of fx from expression (13). But the question is how to derive the elements of the matrix A from (16).

1

1

1

If we observe the expressions (12) we can notice their dependence of the mesh and of the relevant point M' = (x',y'). How to derive the integrals (12) is actually the crucial question in this work. By using the basic computations and the property of additivity we are able to simplify these integrals. For example:

volC,

(x — x')(y — y')dxdy

volC,

(xy — xy' — x'y + x'y')dxdy

Ci

Ci

volC,

IJ xydxdy — J J xy'dxdy — J J x'ydxdy + J J x'y'dxdy

Ci Ci Ci Ci

So, the functions we are dealing with are constant, linear and quadratic functions (or the product of two linear functions).

The first try in obtaining the coefficients a, ft, y and e was TRIANGULATION in R2. This triangulation formula works perfectly accurate with linear and constant functions. Unfortunately, when the integrand is quadratic function the triangulation formula does not give satisfying results.

Another try was by multidimensional numerical integration [5]. The concept of this integration is also settled on a triangular decomposition of a certain polygonal domain. Actually, we are looking for two-dimensional composite quadratures. Briefly, the formula

I (f )

m

3

Ef (j ),

(18)

where aj, j = 1,3 are the midpoints of the edges of a generic triangle T G Th and \T\ is the area of T. Formula (18) absolutely satisfies our demands for desired degree of exactness to 2 and it is suitable for the efficient computations of the integrals (12). In the forthcoming part we will introduce complete deriving of the coefficient aj. Deriving of the coefficients fy, Yj, Sj and £j will be written in shorter form.

Consider a certain quadrilateral cell C (for the sake of simplicity we omit the index j, for a moment). As already mentioned that cell can be divided into two triangles: T and T2 (see Figure 2).

The vertices of the triangle T are: (xo,yo), (xi,yi), (x3,y3). Similary, the vertices of the triangle T2 are (x1,y1), (x2,y2) and (x3,y3); also the midpoints are:

xo + xi yo + yi

2

2

xi + x3 yi + y3

2

2

,T1

xo + x3 yo + y3

2

2

1

1

1

3

1

2

3

and

,T2

xi + x2 yi + y2

,T2

x2 + x3 y2 + y3

2

2

2

2

-T2

xi + x3 yi + y3

T1

2

2

i

2

3

2

To derive the integral a we will do the integrations along the the border of the both triangles Ti and T2. That is how ai and a2 ( the integrals over the border of Ti and T2 respectively)

( *3, y3)

Fig. 2

will be obtained. Considering M'(x',y') as a given point and M(x,y) as the unknown point around M' in which we expand the function, we can start with calculations:

3

1 1 lTll

j=1

(#0 + + — 3x') = + + x3) — x',

1 |TiL , _ , _ „ ,, 1

|Ti| 3 v ' 3

3

j=i

(#1 + X-2 + X3 — 3x') = -(^1 + + X3) — x'.

1 TL _ , _ 9 ,, 1

|T2| 3 v ' 3

Finally:

a = Q'i + Q'2 = + x\ + x3) - x' + ^(rri + x2 + x3) - x' = ^(rro + 2xi + x2 + 2x3) - 2x'. 3 3 3

The coefficients ß, 7, 5 and e can be derived in a similar way. Before we write down all expressions for the derived coefficients we must bring back the omited index j. In this way the coefficients will be specified according to the relevant cell

Q'i = + 2xji + Xj2 + 2^3) - 2./;' . ßj = ^•//,,, + 2//? 1 + yj2 + 2yj3) - 2//',

12222

Tj = g y^'jO '^Xji + Xj2 + 2^3 + XjoXji + Xj\Xj2 + 2Xji'Xj3 + XjoXjs + Xj2Xjs) —

2

+ 2^-1 + Xj2 + 2^3) + 2^2, 1

6' r'"r' ' ' ' ' 12 '

^ = 7: 5 •''/".'//•■ + 2./-? I //? I + Xj2Vj2 + 2xj3yj3) + —•;.'>,//, I + I//,,, + /1 ///:: + 2.r?;;//?| •

y ■

+xm0 + Xjoyji + Xj\]jj2 + ^ + Xj2Vj3 + ^-2) " f (**> + + ^2 + 2^)-

Tp//./•• + 2//, I + ///2 + —.'//.; i + 2.r'//',

12222

= + —///1 + //,2 + 2//?;; + ///"//.; I + //,l//,2 + 2//? I//,;■, + //,,,//,;■, + //,2//,:'. Ï

2 2

rXi'UIr- + 2yji + yi2 + 2yi3) + 2yf.

3.2. Conclusion

Now we are able to obtain the coefficients Kj from the (14). However, in this general case it is worthless to derive explicit formulas for Kj. So, we present these variables as a functions of derived coefficients aj, ftj, Yj, and £j or said in another way, as a functions of x',y' and Xi,iji, i = 0, 3 and this clearly describes the local geometry dependence of the method. In practice, these coefficients are computed and stored in each relevant cell before main simulation

Ko = <Maj , Yj ,Sj ,£j) = po(ih(xi ,yi,x',y ))

Kl = Maj , Yj ,àj ,£j) = <pi(^i(xi ,yi,x',y ))

K2 = P2(aj , Yj ,àj ,£j) = P2(ih(xi ,yi,x',y ))

K3 = P3(aj , Yj ,àj ,£j) = p3(ih(xi ,yi,x',y ))

K4 = P4(aj , Yj ,àj ,£j) ,yi,x',y ))

k5 = = f5(ip5(xi,yi,xf,y')), i = 0,3.

Finally, we can write the x-derivative of the function f at the point M', fx(M') is completely determined by following statement:

E '<•,/< • = fAM') + E '

j j

5 5

U(M') = £ «;C1c, - 0(AJ)) = o M7c, - 0(A*)). j=0 j=0

The corresponding formula for the y-derivative, fy(M') is:

5

j=0

4. The finite volume higher order scheme on the quadrilateral cells

Knowing that the LDLR was defined as: (7), we can extend previous expression in two space dimension on a non-uniform quadrilateral mesh. Then we get:

rij Qij

vol(Cij )

Ci.

Vij (x,y) -

Axij

<Pij(Ç,9)dÇd9 -

ci.

Ay

pij(n, ^)dr/d^

ij

Ax

1 r vy&dWM- 1

dxdy + pij (x,y)-Pij(n, ^)drjd^,

ij

Ci'

Ayij

Ci.

which presents LDLR in two space dimension.

1

1

1

4.1. Flux Balancing and Applying Gaussian Quadrature

Let's now evaluate the conservation law on each quadrilateral cell:

jt J J q(x, y, t)dxdy + JJ f(q(x, y, t)) • n dS = 0. (20)

Cij dCij

Our final task is to find a flux balance along the cell boundary which is needed for describing the dynamics of cell average:

W(t) = _L_ J J q{x,y,t)dxdy. (21)

The flux balance is described by next expression:

i f (q(x,y,t)) • ndS. (22)

dqij(t)

dt vol (Cij)

dCij

How can we compute the right hand side of (22)? One of the possibility is to evaluate the integral in (22) by the fourth order Gauss quadrature formula. To obtain this approximation we need to evaluate the LDL-reconstructed flux function over all four edges of the cell boundary. The sum of the quadratures over all four edges gives the right hand side of equation (22) and at the same time our third order finite volume scheme on a quadrilateral mash. In Figure 3 the Gaussian nodal points are presented and marked as Aki, k,l = 0,3, k = l. For fourth order Gauss quadrature formula on each of the the four edges there are two Gauss nodal points:

Am = Mkl - ii.o.///• - (xk, ///,;•;•. BkI = Mki + ■;■;.'•/. yi) - (xk, ///,;•).

where Mki is the midpoint of a relevant edge whose vertices are (xk,yk) and (xi,yi),

Mki = (xk, Vk) + (xi, Vi)

(Xj3, yj3) ^ V B23 (xJ2' YJ2)

((X2' Y2)- (Xl> Yl))

12 = M12+^12((x2> Y2)- (X1. Yl))

(j o) A01

Fig. 3

A third order flux-balance scheme is presented by the next formula:

1 1 f f(q(x,y,t)) • ndS = M01rl3(A01) + M01rl3(B01) + Ml2rl3{Al2) +

vol(Cij) J

dCij

+Ml2rj(Bl2) + +M23fj(A23) + M23f'ij (B23) + M30rj(A30) + M3on3(B3o) =

3

= Mki [rij (Aki) + rij (Bki)].

k,l=0,k=l

Finally, desired finite volume scheme, third order accurate, on the quadrilateral mesh is achieved by integration of the previous expression.

References

[1] Randall J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge Univ. Press, 2002.

[2] Artebrant R., Sohroll H.J. Limiter-free third order logarithmic reconstruction //J. Sci. Comp. 2006. Vol. 28(1). P. 359-381.

[3] Artebrant R., Schroll H.J. Conservative Logarithmic reconstruction and Finite Volume Methods // Ibid. 2005. Vol. 27(1). P. 294-314.

[4] Schroll H.J., Svensson F. A bi-hyperbolic finite volume method on quadrilateral meshes // Ibid. 2006. Vol. 26(2). P 237-260.

[5] Artebrant R. Third order accurate non-polynomial reconstruction on rectangular and triangular meshes // Ibid. 2006. Vol. 30(2). P. 193-221.

Received for publication 14 April 2010