Technologies of finite volume-finite element method for the solution of convection-diffusion problems on unstructured grids Текст научной статьи по специальности «Физика»

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

Том 7, № 3, 2002

TECHNOLOGIES OF FINITE VOLUME-FINITE ELEMENT METHOD FOR THE SOLUTION OF CONVECTION-DIFFUSION PROBLEMS ON UNSTRUCTURED GRIDS *

E. P. Shurina Novosibirsk State Technical University, Russia O.P. SOLONENKO, T.V. VOITOVICH Institute of Theoretical and Applied Mechanics SB RAS

Novosibirsk, Russia e-mail: [email protected], [email protected], [email protected]

Работа посвящена построению и исследованию схем методов конечных объемов/ конечных элементов на симплициальных сетках для решения задач конвекции —диффузии. Рассматривается проблема точного интегрирования интерполяционных полиномов, построения противопотоковых схем повышенного порядка точности на компактных шаблонах триангуляций. Предлагаемые технологии используются для моделирования вязких течений несжимаемого газа.

Introduction

Finite volume methods on unstructured grids have become widespread tool for the solution of computational fluid dynamics problems. They combine such attractive features of finite element unstructured approaches as local grid refinement realization, adaptation techniques, unstructured multigridding and advantages of classical finite volume difference method (FVDM) [1, 2], namely discrete local conservation, simplicity of approximation and high accuracy of one-dimensional upwind schemes. At the same time some points — development of high-order-accurate genuinely multidimensional upwind schemes, asymmetric nature of the finite volume-finite element approximations, existence of p-version of the finite volume element method — are still open.

Unstructured technologies are rooted in the finite element context, with the advantages of the method being, in particular, inherence of the symmetrical nature of the self-adjoint part of differential operators in their discrete analogues and the possibility of increment of interpolation polynomial degree (p-version of the finite element method [3]). At present, technologies of the finite element methods on simplicial grids are well-settled and include interpolation procedures

*This research has been partially supported under the Grants No. 00-01-00-899, 98-0217810 of the Russian Foundation of Basic Research and also by the Grant 2000-1 of Integration program of SB RAS.

© E. P. Shurina, O.P. Solonenko, T.V. Voitovich, 2002.

and technologies of the inner product approximation with the use of piecewise-polynomial representation of the solution and coefficients of differential equations. Special classes of integration formulas for the exact integration of barycentric coordinates monomials obtained by Eisenberg and Malvern [4] are used as the basis of the finite element approximations.

Attempts at systematization of FVDM approximation techniques and extension of the method to irregular geometries lead to partial combination of the finite element method technologies and discretization principles of the finite volume method. The first class of methods combining ideas from both the finite element and finite volume techniques has been proposed by B. R. Baliga, S. Ramadhyani, C. Prakash and S. V. Patankar in the mid-1970s for the solution of convection-diffusion problems and referred to as control-volume-based finite element methods (CVFEMs) [1, 5]. The control-volume-based finite element methods inherit such principles of the finite element method as specification of element based interpolation polynomials for the solution, source terms and coefficients of the boundary value problems and element-by-element assembly of discrete analogue matrices and vectors.

Both the finite volume element and the control-volume-based finite element method use the family of conforming spaces of piecewise-linear functions over triangles and belong to the class of the cell-vertex finite volume approximations. We shall refer to these classes of discretization techniques as the finite volume-finite elements methods (FVFEMs).

At present the finite volume-finite element methods have no unified technologies of exact integration of element-based interpolation polynomials and their combinations for the following particular reasons. Unlike the finite element method, finite volume-finite element approximations are restricted by the discrete spaces of low degrees: with the use of several types of dual meshes corresponding to interpolation polynomials of high and low degrees, some desirable local conservation properties are violated [5]. The finite volume methods are diverse as characterized by variety of possibilities for defining the weighting functions that are related to both the arrangement of computational points (cell-centered, edge-based, cell-vertex arrangements) and construction of the control volumes (with the use of circumcenter, ortocenter, incenters, or barycenters).

The only paper concerned with the problem is the work of Y. Liu and M. Vinokur [6]. It is important to note that the approach presented uses the exact integration of basis functions to obtain high-order accurate FV schemes; formulas for the exact integration of generic monomials in Cartesian coordinated are used as the base of the approximation. The present paper develops alternative technologies of the exact integration of polynomials in the finite volume-finite element discretization processes on unstructured simplicial grids, we restrict the consideration to the barycentric control volumes.

The second problem we are interested in is the development of high-order accurate upwind schemes in the cell-vertex framework on the compact triangulation stencils. In fact, the comparison of various finite volume and the finite element methods for the convection-dominated boundary value problems comes to the comparison of corresponding upwind schemes.

Finite volume methods due to their ability of sharp capturing shocks and boundary layers are used mainly for the numerical solution of hyperbolic conservation laws. The unstructured finite volume formulations of Jameson and Mavriplis use some form of artificial dissipation [7]. Currently finite volume upwind schemes are dominated by cell-centered schemes based on the solution of local Riemann problem for each cell face [8, 9], the disadvantages being one-dimensional nature and high numerical viscosity of solvers [10]. It should be noted that for the construction of second-order accurate FV schemes both in the upwind biasing and artificial diffusion frameworks essential widening of the stencil support is necessary [7, 8].

Upwind residual distribution schemes proposed by P.L. Roe [11] following the concept of fluctuation also combines principles of both the finite element and finite volume methods. Unlike Riemann-solver-based upwind schemes the residual distribution schemes adequately capture the multidimensional flow structure and allow to obtain second-order approximations on the compact stencils corresponding to conforming linear finite element spaces. These schemes have been validated in both compressible and incompressible cases [10].

Applications of the finite volume schemes to the incompressible flow modeling are relatively rare [5, 12, 13] and corresponding upwind schemes are currently not well-developed and systematized. Some of the (cell-vertex) upwind finite volume formulations actually come to the use of one upstream nodal value of the solution for the approximation of local convection fluxes. Early control volume based finite element formulations (CVFEMs) were based on flow-oriented upwind schemes [5] that successfully reduced the false diffusion but suffered from appearance of negative coefficients in the discrete analogue matrices. Development of unconditionally monotonic mass-weighted upwind schemes in the unstructured CVFEM framework is relatively recent and is an extension of the skewed mass-weighted upwind schemes of Schneider and Raw [13] to the simplicial grids.

Development of high-order accurate finite volume upwind schemes using the same compact stencils as RDS does is of high interest. In the present paper a new principle of weighting the local mass fluxes with the use of some asymmetric profile for development of upwind schemes is proposed.

1. FVFEM discretization of convection-diffusion-reaction equations

1.1. The model problem

Let us consider the following convection-diffusion boundary value problem

dM

dt

+ V (pu0 - A^V0) + y0 = S^ in Qt ,

(1)

0 (■, 0) = 0o (■) in Q,

0 = gD on rD x (0,T),

A^V0n = gN on rN x (0,T),

where Q is a bounded domain in the plane, with the boundary r = rD U rN, QT is the space time cylinder, p is the fluid density, u is the velocity vector, A^ is the diffusion coefficient, S^ is the volumetric source of 0; gD, gN are given real-valued functions.

(2)

(3)

(4)

Q x (0,T)

1.2. Domain discretization

Let Th be a triangulation of Q: Uenerh en = Q, such as the intersection of two distinct triangles is either empty or consists of one common vertex or edge. Assume that computational points and finite element vertexes coincide ("the cell-vertex" arrangement). In this paper we restrict ourselves to dual meshes constructed on centroids of triangles and mid-points of the edges (the Donald meshes), such that each node has the corresponding complete or "incomplete" finite volume bounded by median segments (Fig. 1, b).

Fig. 1. Fragment of primary (FE) and secondary (FV) grids (a); notation for the dual mesh segments and finite element subregions (b), outward unit normals (c).

Let Si denote the boundary of Qj. Assume that all elements and nodes are ordered and a local numbering of nodes is fixed. Let us assume fixed global and local numbering on the set of triangulation nodes and on each finite element.

1.3. Integral conservation form

The integral conservation form of equation (1) is

J dQ + J (pu0 - A^V0) ndS + J y0dQ = J S^dQ

or

j dQ + j J ■ ndS + j Y^dQ = j S^dQ (5)

where J = Jc + Jd is the combined convection-diffusion flux of 0,

Jc = pu0, Jd =

n is the outward unit normal to Si.

Consider a finite element en. Let SHb denote a segment of the dual mesh lying on a median, which emanates from the triangle node nb, and Jvnb, v, nb = 1, 2, 3 is a value of the combined convection-diffusion flux through the median segment Snnb in the direction of the outward normal nvnb to the boundary of the finite volume Qv (Fig. 1, b, c):

Jv nb I Jdnb + J v nbnvnbdS-

S\

nb

It will be sufficient to approximate three of the six introduced fluxes for the element of interest, because of the conservativity property, J21 = — J31, J13 = — J23, J32 = — J12. We shall refer to the fluxes J21, J13, J32 as the determining fluxes of the element of interest.

1.4. Technology of the finite volume-finite element approximations 1.4.1. Stiffness matrix

Let us demonstrate features of the proposed technology on a problem of approximation of the diffusion flux with variable A^.

Using Green's formula, the diffusion flux through a median segment S^ can be represented as follows

Td i ^ a

J

v nb

I ^"dx^y - ^"dydX'

Sn,

nb

Consider the usual finite element space of continuous piecewise polynomial functions

Sh = {v e C°(H): v|en e Pi Vera eTh} ,

(6)

(7)

where P1 is the space of polynomials in two variables of first degree. The basis {^i(x)}N1 of the discrete space Sh consists of continuous piecewise linear functions ^(x) over triangulation Th, being equal to unity at the vertex xi and 0 elsewhere: ^(xj) = 5ij, i,j = 1,N.

While approximating diffusion fluxes we shall assume that FVFEM solution 0h belongs to the finite element space SE, that is

nh SE

{v e C0(H): v|rD = gD} ;

let us use the same discrete space (7) for the local presentation of the diffusion coefficient. Thus, the solution and interpolation function for the transfer coefficient À^ G Sh can be locally represented as follows:

3

0h(x) = ^ (f)m"tjjm(x) = a<p + b^x + c^y,

Aj(x) = \m-tpm(x), x e

i=1

i=1

where ^(x) are local linear basis functions, Am are the local nodal values of À^, m = 1, 2, 3, and coefficients a^, b^, c<p can be expressed in terms of the local discrete unknowns <f1, 02, 03 via the inverse of the Vandermonde matrix containing coordinates of the triangular nodes:

K0.

a0 1 X1 y1 -1 "01 "

= 1 X2 y2 02

_ 1 X3 ya _ . & _

Thus, approximation of the diffusion flux (6) becomes

J

v nb

í l}mdy > k2l 0l + < í ^mdx > foli

i=1 Snb J i=1 J i=1 Snnb J i=1

(8)

It is clear from (8) that it is necessary to introduce special formulas for integration of the local basis functions along the median segments to complete discretization of the diffusion flux. Note that in the case of continuous triangular based linear interpolation, local basis functions coincide with simplex barycentric coordinates Lj, ^(x) = Lj, i = 1, 2, 3.

For every local basis function consider a pair of trapezoids and a triangle, which are constructed using finite element median segments and their images on the surface of the basis

Fig. 2.

function (Fig. 2). Projection of the trapezoids and the triangle onto the coordinate planes gives the following formulas for the exact integration of barycentric coordinates over median segments:

1

Lv dx —

Lv dy

Snb

nb

1l

12 vx)

6 *

5

12

vy,

v — nb, v — nb;

v — nb,

lvy, v — nb.

(9)

Here lvx and lvy are the lengths of the projections of the median segment Sn onto the x- and y-axis, respectively, with the signs corresponding to the anticlockwise integration; v,nb £ {1, 2, 3}. Taking into consideration formula (9), we will have for three determining diffusion fluxes of the element:

J D

J31

l

ly

J

l

2y

1A 5 - 5 -6 Al + 12 A2 + 12 A3

5 A 1. 5 A 12 Al + 6-2 + 12 A3

3 f 1 - 5 - 5 -

y^ k2lil + llx\ 6-1 + 12-2 + 12 A3

i=1

y^ hi il + l

i=1

5 A 1A 5 A 12 Al + 6 A2 + 12 A3

3

E

i=1

k3i ii,

i=1

JD

J

23

5 A 5 A 1 A

-l3y\ ^Al + ^ + 6A3

£k

i=1

5 A 5 A 1 A

2lil + l3x \ 12A1 + 12 + 6

Y,k3iii- (10)

i=1

Finally, to determine coefficients in the first, second and third row of the local diffusion

matrix of the element of interest, the combinations Jf2 — J23

J2c

3 Jc i and Jc i Ji <

are

23 ^ " 12

considered in terms of conservation unknowns. Coefficients at (f)1, (f)2, (f>3 in the combinations give three elements in the current row d™, i = 1, 2, 3 of the local diffusion matrix An

a d

dn d3

1.4.2. Source terms

The volumetric source S<p contribution to the ith component of RHS global vector can be represented as f S^dQ, gathering the contributions from all triangles that share the

n ninen

Snb

nb

3

3

dn

2

vertex i. Let us consider the case of continuous piecewise-linear interpolation of volumetric source function, therefore, at each simplex en £ Th we have:

3

(x) = y1 x £ e„,

m=1

where are values of the source function at three vertexes of the element en,m = 1, 2, 3. Let denote an intersection of the finite element en of interest and the finite volume of the node with local number nb, nb = 1, 2, 3. Let us introduce the formulas for the exact integration of the local form functions over barycentric subregions Q(Fig. 1, b):

— meas en, v = nb,

L dft = < 57 (11)

„ n , —meas en, v = nb.

nn I 108 n

Hence, the contribution of the element en to the global RHS vector can be represented by the following three values:

'11 sVi + -!_ + _!

54 ^ 108 108

= ( + 77^+ TT^) meas en

= (158 +14 +108S<^3) meas

= 1158 +158 + 5iS<^3) meas

being the element s of the local RHS vector f'. 1.4.3. Mass matrix

It is possible to use formulas (11) of the exact integration of barycentric coordinates over barycentric subdomains for the calculation of the mass matrix. Assume that the FVFEM solution is a continuous piecewise linear function over the triangulation Th, £ SE; let us use an averaged value of 7 at each simplex. The approximation becomes

I 70hdH « Yn /" ^ ^m(x)dQ, x £ en, nb = 1, 2, 3. (12)

J J ¿=1

n™ nn

nb nb

Taking into consideration (11) we will have for the local mass matrix:

11n

—7' meas en, i = j,

A?) 54 (13) 7 ij ' n • / • -Y' meas en, i = j,

I 108 ' n

where 7" is the averaged value of 7 on the element en, i, j = 1, 2, 3.

Thus, the implicit time discretization scheme leads to the following local mass matrix:

i11 Pi .

—- ——meas en, i = j,

54 AP; . =. (14)

--:—meas en, i = j

108 At, n' ^

n n

n n

n n

(a non-lumped scheme); corresponding RHS vector contributions are calculated as A^(fs-1, where pi are local values of density, Ats — a time step, (fs-1 — local vector of the conservative variable from the previous time step. The lumped discretization of the non-steady terms will be characterized by the following diagonal mass matrix

AL = diag

A L) 1 , (A L) (A L ,

(15)

1.5. Multidimensional upwind schemes on unstructured meshes 1.5.1. Concept of local weights matrices

In this section we will be interested in the development of high-order accurate upwind schemes on the compact triangulation stencils corresponding to the conforming linear finite element spaces [14, 15]. Let us introduce the notion of a cell-vertex upwind scheme internal with respect to simplexes. An upwind scheme issaid to be internal with respect to simplexes if for the approximation of three determining convection fluxes of an element only the local nodal values of the scalar variable at the element are used. From the realization point of view this means the possibility of the usage of the most simple grid data structures. Thus, we will be interested in the construction of internal with respect to simplexes upwind schemes in the cell-vertex FVFEM framework.

One of the determining convection fluxes at the element of interest will be approximated as follows:

Jnb = Jp(u ■ nvnbdS « („¿Jpu ■ nvnbdS, (16)

Sn, Sn,

nb nb

where {vnb} £ {{12}, {23}, {31}}. In fact, for every median segment Snb we have to express a value (f)nb, that is the value of the scalar quantity used for the approximation of the determining convection flux through the segment (Fig. 3) in terms of the local discrete unknowns as their linear combination

4>nb = Y, ßnbk (17)

i= 1

3

h

Fig. 3. Notation used for the local conservative unknowns vectors, (f and (.

Upwind character. While approximating the convection flux through median segments of the element of interest, let downstream nodal values of (h be used in the combination (17) with

less weights Thus, some desirable upwind bias exists, resulting in some false diffusion level. It should be noted that the above general upwind character principle is less categorical than one usually introduced as an upwind character design principle for Residual Distribution Schemes by Roe, Deconinck and co-authors [11]. Remind that the upwind character principle of RDS prescribes for the distribibution coefficients of downstream nodes being equal to exactly zero. Upon our experience, the schemes biasing the "information center" of the simplex of interest but using all nodal values 0j for the approximation of convection fluxes have much less false diffusion than schemes projecting the information center onto an upstream edge (or replacing the center by an upstream node of the simplex) [12].

Let us introduce the concept of weights local matrix of an upwind scheme. B matrix is called a weights local matrix of the upwind scheme on the element en if its first, second and third rows contain the coefficients fi^ corresponding to the median segments S^, S2?, S3 , respectively. For example, the following weights local matrix

" 10 0 B = 1 0 0 100

prescribes for all determining convection fluxes of the element to be approximated with the use of the upstream node with local number 1, while approximation scheme with the matrix

111

B

3 1

3 1

3

3 1

3 1

3

3 1

3 1

3 J

will have the stability properties of the central scheme.

Therefore, development of an upwind scheme comes to the determination of matrix relations of the following form

0 = B0 (18)

at each simplex en £ Th, where 0 = (01 02 03 ) is the local vector of the conservative

unknown values which are used for the approximation of the determining convection fluxes through the median segments S^, S2?, S3*"; 0 = (01 02 03 ) is the local vector of conservative

unknowns. Once the local weight matrices are determined at each simplex, the corresponding upwind scheme is completely specified, not depending on the accuracy of the approximation of the local mass fluxes in (16). We shall refer to this principle as to the design principle of separate approximations of the local mass fluxes and the local values of the conservative unknown at the dual mesh segments.

Let us use the following notations for the local mass fluxes:

m 1 = / p0u ■ n31dS, m2 = / p0u ■ n12dS, m3 = / p0u ■ n23dS.

(19)

It is obvious that the local mass fluxes (19) have to be calculated as exact as possible. Assume the components of vector p(x)u(x) varying linearly at each simplex:

p(x)u(x) = p(xi)u(X¿)-0m(x)

(20)

i=1

Sn

cn S2

cn S3

where p(xi) are nodal values of density, «(x^) = (u(xi),v(xi)) are nodal velocity vectors and -0m(x) are local linear form functions. Using the integration formula (9) introduced above one can obtain the following approximation of the mass fluxes (19):

(1 5 5 \

m 1 ~ ( 6P(x 1)u(x 1) + 12P(x2)u(x2) + 12P(x3)u(xsH hy —

/1 5 5 \

— ( 6P(xi)v(xi) + 12P(x2)v(x2) + 12P(x3)v(xsn hx,

5 5 1

«2 ~ ( 12P(X1)u(Xl) + 12p(x2)u(x2) + 6P(Xa)u(x3^ l2y -

5 5 1

12P(x 1)v(X 1) + 12P(X2)v(X2) + 6P(X3)v(X3^ l2x,

5 1 5

«3 ~ I 12P(X 1)u(X 1) + 6P(X2)u(X2) + 12P(X3)u(X3^ l3y -

(5 1 5

V 12P(X 1)V(X 1) + 6P(X2)v(X2) + 12P(X3)v(X3^ l3x,

where lvx and lvy are the lengths of the projections of the median segment S^ onto the x- and y-axis, respectively, with the signs corresponding to the anticlockwise integration; v, nb G {1, 2, 3}.

For the discretization of the momentum equations in the compressible conditions it can appear essential for shock capturing to have a possibility of linear representation of both velocity components and density at each simplex:

P(x) = P(Xj)V>m(X), U(X) = ^ u(Xj)V>m(X),

i=1

i=1

instead of (20). Special formulas for the exact integration of barycentric coordinates monomials of the second order, which will necessary in this case are introduced in Section 2.5.3.

Finally, expressing the combinations J^2 — J£3, J£3 — JSq and Ji^ — in terms of local conservation unknowns, we will have for the local convection matrix of the element:

A c

ß21«2 - ß1 «3 ß3« 3 - ß1 m 1 ß^m 1 - ß2«-2

ß|«2 - ß|«3 ß|«3 - ßi2« 1 ß2« 1 - ß|«2

ß3 m« 2 - ß3« 3 ßfm3 - ßfm 1 ß^rn 1 - ß|m2

(21)

where weights Plnb are the elements of the local weights matrix of the upwind scheme. The weights local matrices will be calculated with the use of various upwind principles, determining a scheme of approximation of convection fluxes.

1.5.2. Mass Weighted Upwind scheme

The MAW scheme defines a mass-weighted average of 0 at each of the median segments of the element of interest following [12], thus projecting the information center onto an upstream

edge (or replacing the center by an upstream node of the simplex). The treatment of the local mass fluxes is performed depending on the local mass fluxes values:

01

02

03

fp02 + ( - fp)03, where fp = min [max(m2/rn 1, 0)] if mn 1 > 0,

fp03 + ( - fp ) 02, where fp = min [max(mh,3/fn 1, 0)] if mn 1 < 0,

fp03 + ( fp)01, where fp = min [max(mh,3/rh2, 0)] if mn 2 > 0,

fp01 + ( where fp = min [max(m 1/rh2, 0)] if mn 2 < 0,

fp01 + ( - fp ) 02, where fp = min [max(mh, 1/rh3, 0)] if m 3 > 0,

fp02 + ( fp)01, where fp = min [max(mh,2/rh3, 0)] if mm 3 < 0.

(22)

Let us obtain the expression for the corresponding local weights matrix. Let us write equations (22) in the matrix form 0 = C0 + C0, where all non-zero elements of matrices C, C are placed out of its diagonals depending on the local mass flow rates m 1; m2, m2 signs. Therefore, the local weights matrix of the MAW scheme can be found as

B

E - C

-i

C,

thus, for realization of the MAW upwind scheme a 3 x 3 matrix have to be inverted on each element.

1.5.3. Modification of the Flow Oriented Upwind scheme

The interpolation function used in the Flow-Oriented Upwind scheme [5] is defined using a local flow-oriented coordinate system (X, Y), with the origin in the barycenter of an element and X axis oriented along the element-averaged velocity Uav (Fig. 4). The convection fluxes are approximated by assuming that over any element 0h varies exponentially in the direction of the average velocity vector and linearly in the normal direction:

0h = A + + CY, where

£ =

A,

pUa

exp

- 1

PeA = pUa

XX

Xmax Xmin

A,

Xn

PeA (X — Xmax) (Xmax Xmin)

max(Xi, X2, X3), Ye

(23)

(24)

max(Yi, Y2, Y3),

Fig. 4. Auxiliary coordinate system used in FLO upwind scheme.

and Xi, Yi, i = 1, 2, 3 are coordinates of the simplex nodes in the new coordinate system. Introduction of auxiliary variable £ results in the modification of the linear shape function (23) and appearance of numerical diffusion only in the averaged flow direction.

Let us use the profile (23) for the local weighting of the mass fluxes mi, m2, m3, thus taking into consideration all vertexes of the element, unlike the MAW scheme. The new upwind scheme will be called the Modified Flow-Oriented Upwind scheme (FLOM). While constructing the local weights matrix for the FLOM scheme, we shall represent the combinations 02, 03 as images of the median segments Sn, Sn, Sn midpoints on the surface (23). Thus, the combinations are close to the nodal values of the conservative unknown in the upstream nodes (Fig. 5). From the interpolation condition we have for the coefficients A, B, C:

(25)

" A " " 1 6 YT -1

B = 1 6 Y2 02

C 1 Ca Y3 . 03 .

Fig. 5. FLOM: weighting of the local mass fluxes with the use of the asymmetric profile.

On the other hand, components of the local vector 0 = 02 03^ being the images of the median midpoints, can be found as:

= A + B£i + CYi, 02 = A + ££2 + CY2, 03 = A + ££3 + CY3, (26)

where £i, Y, i = 1, 2, 3 are coordinates of the median segment middles in the new coordinate system. Substitution of (25) to (26) allows to express 0 in terms of the local unknowns:

0nb = (kii + k2i£nb + &3iY^ 01 + (ki2 + k22£nb + fc^Y«^ 02 + (fci3 + k^fnb + fe^Y^ 03 .

Thus, the component (nb, i), nb, i = 1, 2, 3 of the FLOM local weights matrix becomes

ffnh = k1i + k2i£nb + k3iYnb. (27)

Resulting convection local matrix of the element is calculated using (21) with the weights (27). One of the advantages of the FLO profiles lies in its flexible dependence on local Peclet number Pe^, so the shape can vary from nearly linear at the rim of the flow to boundary layer form at its kernel.

1.5.4. Analogue of the FLOM upwind scheme for the cell-centered arrangement

Let us construct a FLOM type scheme for the cell-centered finite volume discretization of the conservation equation (1). Consider a control volume consisting of an elementary triangle Tj with boundary Sj. Integration of the conservation equation (1) over Tj gives

j dQ + y J ■ ndS + y = j S^dQ.

Ti Si Ti Ti For the approximation of convection flux some average 0 has to be expressed in terms of values of , , Consider the triangle (l j k) constructed with the use of the barycenters of the neighboring triangles. Let a, b, c be intersections of lines connecting node i with centers of the neighboring cells and edges of the control triangle ej (Fig. 6). The intersection points seemed to lie on the edges of the control simplex for quite general triangulation. The exponential profile

Fig. 6. Analogue of the FLOM scheme for the cell-centered arrangement.

(23) can be constructed with the use of three centers of adjacent triangles and velocity vector corresponding to the node i (the velocity vector will play the same role as the averaged velocity vector at an element in vertex-based approaches), as shown in Fig. 6. The average values

which are necessary for the approximation of the convection fluxes, can be picked as images of points a, b, c at the profile, which is exponential in the direction of Uj and linear in the direction normal to Uj. The scheme combines features of both central and upwind differencing. On the one hand, the upstream nodes are used with more significant weights. On the other hand, the control node i is used only in the sense of the velocity vector Uj analyzed, hence, the average , are expressed in terms of , and the scheme characterized by zero

diagonal contribution in the steady-state case.

In the light of this fact the scheme proposed seems to be suitable for the construction of some composite schemes as follows. In finite element literature several examples are reported [16] concerning the development of a composite upwind scheme by means of calculation of local convection contribution as linear combination of the form

A = ( 1 — positive + less diffusive,

where Apositive is the local matrix corresponding to some positivity preserving but crude

upwind scheme and Aiess diffusive is the local matrix calculated according to some less diffusive but oscillating scheme. Following [16], the weighting coefficient a is specified such that the resulting matrix is maintained as an M-matrix.

By analogy with the FEM approach, one can calculate contribution of the control triangle Ti as combination of local contribution vectors corresponding to different schemes:

q = (1 - a)qpositive + aqFLOM' where a is chosen in such a way to ensure non-positivity of non-diagonal convection contributions to the ith row of the global matrix.

1.6. Exact integration of barycentric coordinates monomials of second order over median segments

Consider the control element (Fig. 7). With the use of barycentic coordinates of the centroid,

(o.o,i)

Fig. 7. Barycentric coordinates of the simplex points used for the exact integration of barycentric coordinates monomials over dual mesh lines.

midpoints of the edges and midpoints of the median segments one can obtain the following integration formulas:

27/fci '

v

Lv LMd£ = <

7

Ï08 19

Ï08

kÇ,

kÇ,

ß = k, v = ß, k = ß or k = v = ß, k = ß, k = v,

(28)

where is the length of the projections of the median segment S% onto the £ axis, with the signs corresponding to the anticlockwise integration; £ G {x,y}, k = 1, 2, 3. While performing discretization it is more convenient to use the following tables, which are equivalent to (28). For the integration over and SJ respectively, we have:

L1 L2 L3 L1 L2 L3 L1 L2 L3

L1 1 7 7 L1 19 7 19 L1 19 9 7

27 108 108 108 108 108 108 108 108

L2 7 19 19 L2 7 1 7 L2 19 19 7

108 108 108 108 27 108 108 27 108

L3 7 19 19 L3 19 7 19 L3 7 7 1

108 108 108 108 108 108 108 108 21

v

2. Numerical experiments

2.1. Stability analysis on the solutions of boundary layer type

In this section the results of numerical experiments are presented for CVFEM P1/FLOM and CVFEM P1/MAW techniques (piecewise linear interpolation for 0 in diffusion terms and FLOM (MAW) upwind schemes for the discretisation of convection fluxes, respectively) from its comparison and accuracy estimation point of view.

Test problem 1. The first problem considered in this section is widely used for stability analysis of discretization techniques of convection dominated problems:

dX0+0 - ui- i) * <»' i>-

0(0, y) = sin ny, 0(1, y) = 0(x, 0) = 0(x, 1) = 0. The exact solution has the boundary-layer form in the neighborhood of x = 1:

0(x y) = ea-yb («*+" - e^"") ,

2 \ 1/2 / 2 \ 1/2 U I u2 2\ u / u2 2X

" 2a + + n) , b =2^ + ^

Two types of techniques - CVFEM P1/MAW, CVFEM P1/FLOM were applied to this problem for different local Peclet number, with results depicted in Fig. 8, a, b. While the MAW

1

0.9 0.8 0.7 0.6

io

0 5

0.4 0.3 0.2 0.1 0i

Pe„=1/6

o CVFEM P1/MAW x CVFEM P1/FLOM — exact solution

0.5 x

a

1

0.9 0.8 0.7 0.6

io

0 5

u

0.4 0.3 0.2 0.1 0

o CVFEM P1/MAW x CVFEM P1/FLOM — exact solution

b

Fig. 8. Computed solutions of 0(x, 0.5) for different upwind schemes: different local Peclet numbers, 21 x 21 uniform triangulation (a), 51 x 51 uniform triangulation (b).

scheme is positivity preserving by construction (at the level of element contribution [12, 13]), the FLOM scheme does not meet the requirement of positivity. Existence of some positive non-diagonal convective contribution which is characteristic of FLOM more probably results in local extrema if simplexes with significantly varying volumes share the same node (Fig. 9, a). In Table 1 the computed global L2 error norms are presented together with experimental rates of convergence for different local Peclet numbers. The FLOM and the MAW schemes are seen to have second and first order behavior, respectively, but every time with the growth of the

0.25

0.5

0.75

1X-0.9 0.8 0.7 0.6

5

. 0.5

!

0.4 0.3 0.2 0.1

00

o CVFEM P1 / MAW x CVFEM P1/FLOM — exact solution

1X-0.9 0.8 0.7 )0.6

7 0.5

!

0.4 0.3 0.2 0.1

00

o CVFEM P1 /MAW x CVFEM P1 /FLOM - exact solution

1,

a b

Fig. 9. Appearance of the local maximum in the neighborhood of boundary layer region while using FLOM scheme and unconditionally monotonie solution given by MAW scheme, 21 x 21 non-uniform triangulation (a); more accurate capturing of the sharp profile given by the MAW scheme in conjunction with a uniform grid refinement (b).

Table 1

Discrete L2 Error norm for the solutions computed with the use of the CVFEM P1/FLOM and CVFEM P1/MAW schemes for different local peclet numbers and estimated rate of convergence a

1/h CVFEM P1/MAW a CVFEM P1/FLOM a

PeA = 1/6

10 6.812e-3 — 9.934e-4 —

20 3.241e-3 1.07 2.632e-4 1.916

40 1.570e-3 1.045 6.714e-5 1.971

80 7.708e-4 1.026 1.692e-5 1.988

160 3.816e-4 1.014 4.249e-6 1.994

320 1.898e-4 1.007 1.069e-6 1.990

PeA = 5/6

10 2.277e-2 — 6.093e-3 —

20 1.642e-2 0.472 1.799e-3 1.760

40 9.730e-3 0.755 4.689e-4 1.940

80 5.279e-3 0.882 1.192e-4 1.976

160 2.747e-3 0.942 3.009e-5 1.986

320 1.401e-3 0.971 7.567e-6 1.992

PeA = 25/3

10 7.491e-3 — 8.485e-3 —

20 1.052e-2 -0.490 9.929e-3 -0.227

40 1.371e-2 -0.382 8.445e-3 0.234

80 1.378e-2 -0.017 4.154e-3 1.0236

160 1.021e-2 0.442 1.340e-3 1.632

320 6.178e-3 0.724 3.588e- 4 1.901

local Peclet number the approximation properties of both schemes are rapidly getting worse at a fixed reasonable grid resolution. Experimental rate of convergence a is computed as

a = log2

-

||0 - 0h/2||

X

X

X

X

X

0.25

0.5

0.75

0.25

0.5

0.75

with the use of approximate solutions 0h, 0h/2, corresponding to a triangulation Th and Th/2 twice as dense; the discrete norm is defined as follows

1/2

-

£ I (0 - 0h)2 dxdy

where the integral over the nth element is approximated by the 1-point Gauss barycenters based quadrature.

While analyzing influence of grid refinement on the accuracy of the schemes investigated, it can be pointed out that Mass Weighted scheme in combination with grid refinement gives better results than Modified Flow Oriented scheme (Fig. 9, b). Taking into consideration the possibility of appearance of unphysical maxima in FLOM solutions, the MAW with an adaptive strategy and a local grid refinement algorithm should be recommended for real problems.

It is well known that expected false diffusion contribution often depends on the main flow direction with respect of the grid lines. In this work two patterns of discretization were selected

Fig. 10. The two types of triangulation patterns used for numerical experiments.

(Fig. 10). The results obtained by both CVFEM P1/MAW and CVFEM P1/FLOM were found to be absolutely independent of the triangulation pattern used.

Test problem 2. The second problem has been considered to make a comparison between CVFEM and some of quadrilateral-based FEM high-resolution upwind techniques [16]:

50 50

U— + V —

dx dx

'd20 d 20 dx2 dy2

(x,y) G Q = (0,1) x (0,1),

0(0, y)

The analytical solution:

sinny sinnys

exp

L2

(y - ys) , 0(1, y) = 0(x, 0) = 0(x, 1) = 0.

0(x, y)

sin(ny) exp(r2x) exp [(v/2)(y - ys)] {1 - exp [(ri - r2)(x - 1)]}

{1 - exp(ri - r2)} sin(nys)

where h is the grid size and

U

4hcot0, v = 4h 0 = tan-1(U), ys =

, 1 if

1--tan-M — ,

n V V /

r1

U

+ Vu2 + v2 + 4n2

r2

U

- Vu2 + v2 + 4n2

The local Peclet numbers are h-independent:

1

Pex

cot0, Pey

1

4'

For small local Peclet number the quadrilateral based FEM composite scheme proposed in [16] was found to give a little better approximation than the CVFEM upwind schemes, but CVFEM P1/MAW and CVFEM P1/FLOM give essentially the best results for moderate and large Peclet numbers (Tables 2-4). It is also clear from Figs 11, a and 11, b that the CVFEM schemes allow for more exact treatment of sharp profiles in the solutions.

v

--------------------FEM, Legendre polynomial

------------FEM, characteristic Galerkin

---FEM, composite scheme

ooooo CVFEMP1/MAW xxxxx CVFEM P1/FLOM --exact solution

j. 1.2

/! ; ! 1.1

"x X 0.9

• \ i \ 0.8

\ ■ „0.7

V!

il! °0.6

il! il ^0.5

il! il! 0.4

t 0.3

i : 0.2

0.1

—1—1— » 0

--------------------FEM, Legendre polynomial

------------FEM, characteristic Galerkin

---FEM, composite scheme

ooooo CVFEM P1/MAW xxxxx CVFEM P1/FLOM --exact solution

0.25

0.5

0.75

0

0.25

0.5

0.75

b

a

Fig. 11. Comparison of CVFEM P1/FLOM, CVFEM P1/MAW, finite element model of Rice and Schnipke, Legendre polynomial finite element model and high-resolution composite scheme of Sheu. d = 2° (a), d = 1° (b).

Table 2

Discrete L2 error norm for the solutions computed with the use of the composite scheme [16], for different local peclet numbers

1/h 9 = 15°, Pex = 0.033 9 = 2°, Pex = 7.159 9 = 1°, Pex = 14.322

10 2.125e-3 7.522e-2 5.800e-2

20 1.660e-3 5.132e-2 3.977e-2

30 8.976e-4 4.694e-2 3.703e-2

40 4.757e-4 4.507e-2 3.589e-2

50 3.036e-4 4.332e-2 3.465e-2

Table 3

Discrete L2 error norm for the solutions computed with the use of the CVFEM P1/MAW scheme, for different local peclet numbers

1/h 9 = 15°, Pex = 0.033 9 = 2°, Pex = 7.159 9 = 1°, Pex = 14.322

10 1.780e-2 2.532e-2 1.475e-2

20 1.884e-2 1.772e-2 1.003e-2

30 1.593e-2 1.375e-2 7.801e-3

40 1.346e-2 1.137e-2 6.522e-3

50 1.157e-2 9.850e-3 5.733e-3

Table 4

Discrete L2 error norm for the solution computed with the use of the CVFEM P1/FLOM scheme, for different local peclet numbers

1/h 9 = 15°, Pex = 0.033 9 = 2°, Pex = 7.159 9 = 1°, Pex = 14.322

10 3.028e-3 1.914e-2 1.487e-2

20 2.030e-3 1.284e-2 9.411e-3

30 1.476e-3 9.609e-3 6.988e-3

40 1.103e-3 7.562e-3 5.514e-3

50 8.473e-4 6.159e-3 4.519e-3

2.2. Laminar flow over a backward-facing step 2.2.1. Comparison of MAW and FLOM upwind schemes

The MAW and FLOM schemes give the following results for the laminar flow over backward-facing step flow benchmark. Experimental data of Armaly et al. [17] are used for the comparison. The case considered is Re = 389 — the maximal value, for which the experimental data are reported with the planar laminar regime of the flow. Figure 12 shows the geometry and a sketch of the flow field. In accordance with [17], upstream channel height h = 5.2 mm, step height S = 4.9 mm. The outlet boundary is at a distance of 32S downstream of the step. Components

45--325

Fig. 12. Geometry and schematic of the laminar flow over a backward facing step benchmark.

x/S = 0. 2.55 3.06 3.57 4.18 5.41 7.14 7.76 8.52 9.74 13.57

Fig. 13. Laminar flow over a backward-facing step, Re = 389. Experimental and calculated velocity profiles for Re=389 at different locations; "o" — experiment Armaly et al., "—" — the CVFEM P1/MAW solution (a), CVFEM P1/FLOM solution (b).

of viscous stress tensor were imposed to vanish at the outlet boundary, thus, for the edges

approximating outlet boundary, only corresponding convection fluxes are approximated. Two types of the inlet boundary placement are considered: four times upstream of the step and at the line of sudden expansion. The parabolic velocity profile approximating experimental one was imposed at the inlet boundary.

The first series of numerical experiments has been performed in the computational domain limited to the region downstream of the expansion, with the CVFEM P1/FLOM results much more closer to the experimental one in comparison with CVFEM P1/MAW. Longitudinal velocity profiles calculated with the use of different upwind schemes at various sections behind the step have been plotted in Figs 13 a and b. The reattachment point computed is 7.14S for the CVFEM P1/MAW and 8.22S for the CVFEM P1/FLOM technique. The result is quite close to the experimental one 7.81S. The velocity profiles at two axial stations, the first at the rim of the main vortex and the second in the neighborhood of the reattachment point computed with the use of different upwind schemes are shown in Figs 14 and 15. Convergence of predicted velocity profiles to the experimental data with the mesh refinement is presented. The use of the modified flow-oriented upwind scheme clearly results in much better prediction of the experimental values than the use of first-order accurate MAW scheme.

Fig. 14. Velocity profiles at x/S = 4.18.

2.2.2. On the mutual compensation of the entrance effect and the false diffusion effect

The so-called entrance effect for the laminar flows over backward facing step has been demonstrated first by Barton [18] and consists in the difference between the finite difference solutions obtained with the different placement of the inlet boundary with respect to sudden expansion. It has been shown by the author of [18], that for the calculations where the domain of interest is limited by the wide channel, the main vortex appears to be longer than for the same calculations taking into consideration the small channel. On the other hand, remind that the use of highly diffusive upwind schemes amounts to the introduction of additional viscosity in the flow, thus we will have the smaller reattachment lengths.

The essence of our contribution lies in the demonstration of the interaction (mutual compensation) of the entrance effect, demonstrated by I. Barton, and numerical viscosity effect introduced by unstructured upwinding. Several series of computational experiments have been

Fig. 15. Velocity profiles at x/S = 7.14.

performed for the solution of the laminar flow over backward-facing step and it has been shown that the following is true. If the computational domain is limited to the region downstream of the expansion, then erroneous arguments of high accuracy of reattachment point prediction can be obtained as a result of the mutual compensation of the entrance effect and numerical viscosity.

Thus, the numerical modeling of the flows over a backward-facing step is the case, then experimental data are in use for the validation of codes, but the physics of the problem is not well understood at the moment.

x/S

Fig. 16. Laminar flow over a backward-facing step. Pressure distribution along the centerline of the small channel.

The reattachment points computed with respect to the small channel are 6.79S and 7.57S for the CVFEM P1/MAW and CVFEM P1/FLOM techniques respectively (triangulation for the wide channel is the same as in the previous section, and a matching grid is used for the small one). Thus, the CVFEM P1/MAW scheme, which had resulted in the underestimation

of the reattachment length (7.14S), gives even worse result. At the same time results of the CVFEM P1 /FLOM scheme, which had lead to the exceeding values of the reattachment length (8.22S), have become closer to the experiment since the calculations are performed in the correct manner.

The pressure fields calculated with the use of different upwind schemes and four pressure distributions along the centerline of the small channel calculated with and without consideration of the small channel with the use of the MAW and FLOM schemes are presented in Fig. 16. It is clear from the figures that the pressure field downstream of the step essentially depends on the upwind scheme used, especially in the neighborhood of the reattachment point; the entrance effect and the false diffusion effect are comparable with each other.

3. Conclusions and future works

With the unsteady incompressible flows modeling in mind, we consider FVFEM technologies for approximation of convection-diffusion problems on unstructured simplicial grids. The paper is concerned, in particular, with the design and investigation of higher-order-accurate upwind schemes for complex incompressible flow modeling by means of cell-vertex finite volume method based on simplicial grids. New finite-element-like technology of exact integration of polynomials in finite volumes methods proposed by the authors and based on the exact integration of simplex barycentric coordinates monomials is used as the basis of the accurate approximations. Concept of weights local matrix of an upwind scheme and the principle of weighting of local mass fluxes with the use of asymmetric profile have been introduced first in cell-vertex finite volume context. It has been shown experimentally for the solutions of boundary layer type that with the use of asymmetric profiles of the Flow-Oriented Upwind scheme, exponential in the direction of the flow and linear in the normal direction, proposed schemes are second-order accurate while Mass-Weighted scheme (MAW) has first order behavior. A comparison of MAW and proposed Modified Flow-Oriented upwind scheme (FLOM) is presented for the complex incompressible flow modeling, with the results that are much more accurate for the last one. Excellent agreement of FVFEM P1/FLOM solutions with experimental data of Armaly et al. for laminar flow over backward-facing step is demonstrated. For the benchmark — laminar flow over backward-facing step problem — an interaction (mutual compensation) of the entrance effect demonstrated by I. Barton, and numerical viscosity effect introduced by unstructured upwinding has been demonstrated first.

Acknowledgements. The authors would like to acknowledge Prof. Dr. V. P. Il'in, Prof. Dr. Yu. M. Laevskii, and Prof. Dr. D. Kroner for inspiring discussions.

References

[1] PATANKAR S. V. Numerical Heat Transfer and Fluid Flow. Washington: Hemisphere, 1980.

[2] Samarskii A. A., Vabishevich P.N. Numerical Solution of Convection-Diffusion Problems. M.: Editorial URSS, 1999.

[3] BABUSKA I., SzABO B. A., Katz I.N. The p-version of the finite element method // SIAM J. Numer. Anal. 1981. Vol. 18, No. 516.

[4] Eisenberg M.A., Malvern L. E. On finite element integration in natural coordinates // Intern. J. Numer. Methods Eng. 1973. Vol. 7, No. 574.

[5] PRAKASH C., PATANKAR S. V. A control volume-based finite-element method for solving the Navier — Stokes equations using equal-order velocity-pressure interpolation // Numer. Heat Transfer. 1985. Vol. 8, No. 259.

[6] Liu Y., VlNOKUR M. Exact integration of polynomials and symmetric quadrature formulas over arbitrary polyhedral grids //J. Comput. Phys. 1998. Vol. 140, No. 122.

[7] JAMESON A., MAVRIPLIS D. Finite volume solution of the two-dimensional Euler equations on a regular triangular mesh // AIAA J. 1987. Vol. 26, No. 824.

[8] BARTH T.J., JESPERSEN D. C. The design and application of upwind schemes on unstructured meshes. AIAA Paper 89-0336.

[9] DELANAYE M., ESSERS J. A. Quadratic-reconstruction finite volume scheme for compressible flows on unstructured adaptive grids // AIAA J. 1997. Vol. 35, No. 631.

[10] DECONINCK H., DEGREZ G. Multidimensional upwind residual distribution schemes and applications // Proc. of Second Intern. Symp. on Finite Volumes for Complex Appl., July 19-22, 1999, Duisburg, Germany. P. 27-40.

[11] ROE P. L. Linear advection schemes on triangular meshes. CoA Report 8720, Cranfield Inst. of Tech., 1987.

[12] RiDA S., McKENTY F., MENG F. L., REGGIO M. A staggered control volume scheme for unstructured triangular grids // Intern. J. Numer. Methods Fluids. 1995. Vol. 25, No. 697.

[13] SCHNEIDER G.E., RAW M.J. A skewed positive influence coeffitient procedure for control-volume-based finite-element convection-diffusion computation // Numer. Heat Transfer. 1986. Vol. 9, No. 1.

[14] VOITOVICH T. V., SOLONENKO O.P., SHURINA E. P. On the construction of high-order-accurate upwind schemes on the compact triangulation stencils // Proc. of Fourth Summer Conf. "Numerical Modelling in Continuum Mechanics", July 31 - August 3, 2000, Prague, Czech Republic. P. 315-326.

[15] Shurina E.P., VOITOVICH T.V. New classes of integration formulas for CVFEM discretization of convection-diffusion problems // Proc. of Second Intern. Symp. on Finite Volumes for Complex Appl., July 19-22, 1999, Duisburg, Germany. P. 377-384.

[16] Sheu T. W. H., Wang S. K., TSAI S. F. Development of a high-resolution scheme for a multi-dimensional advection-diffusion equation //J. Comput. Phys. 1998. Vol. 144, No. 1.

[17] Armaly B., Durst F., Pereira J. C. F., SCHOENUNG B. Experimental and theoretical investigation of backward-facing step flow //J. Fluid Mech. 1983. Vol. 127, No. 473.

[18] BARTON I. E. The entrance effect of laminar flow over a backward-facing step geometry // Intern. J. Numer. Methods Fluids. 1995. Vol. 25, No. 633.

Received for publication January 31, 2002