SciELO - Scientific Electronic Library Online

 
vol.31 issue3Numerical analysis of the nonlinear subgrid scale methodHow good are MatLab, Octave and Scilab for computational modelling? author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

Share


Computational & Applied Mathematics

On-line version ISSN 1807-0302

Comput. Appl. Math. vol.31 no.3 São Carlos  2012

http://dx.doi.org/10.1590/S1807-03022012000300004 

A weighted mass explicit scheme for convection-diffusion equations

 

 

Vitoriano Ruas*

Research associate, Departamento de Engenharia Mecânica, PUC-Rio - Catholic University of Rio de Janeiro, Brazil. Research collaborator, Institut Jean le Rond d'Alembert, UPMC - Université Pierre et Marie Curie, Paris 6, France. Visiting professor, Instituto de Ciências Matemáticas e Computação, USP - University of São Paulo, São Carlos, Brazil. E-mails: vitoriano.ruas@pq.cnpq.br / vitoriano.ruas@upmc.fr / vitoriano@icmc.usp.br

 

 


ABSTRACT

An explicit scheme based on a weighted mass matrix, for solving time-dependent convection-diffusion problems was recently proposed by the author and collaborators. Convenient bounds for the time step, in terms of both the method's weights and the mesh step size, ensure its stability in space and time, for piecewise linear finite element discretisations in any space dimension. In this work we study some techniques for choosing the weights that guarantee the convergence of the scheme with optimal order in the space-time maximum norm, as both discretisation parameters tend to zero.

Mathematical subject classification: Primary: 65M60; Secondary: 76Rxx.

Key words: Convection-diffusion, convergence, explicit scheme, finite elements, piecewise linear, time-dependent, weighted mass.


 

 

1 Introduction

This work deals with an explicit scheme introduced in [6] , for the numerical time integration of the convection-diffusion equations, discretised in space by techniques based on variational formulations such as the finite element method. In this framework, since the mid-eighties, the most widespread manner to deal with dominant convection has been the use of stabilizing procedures based on the space mesh parameter, among which the streamline upwind Petrov-Galerkin (SUPG) technique introduced by Hughes & Brooks (cf. [2] ) is one of the most popular.

The author and collaborators studied in [7] a contribution in this direction, based on a standard Galerkin approach, and a space discretisation of theconvection-diffusion equations with piecewise linear finite elements, combined with a non standard explicit forward Euler scheme for the time integration. The main theoretical result in that work, states that the numerical solution is stable in the maximum norm in both space and time (and even convergent with order h|lnh| if the mesh is of the acute type [10] ), provided that roughly the time step is bounded by the space mesh parameter h multiplied by a mesh-independent constant, for a high Péclet number. As it should be clarified, the scheme under consideration follows similar principles to the one long exploited by Kawahara and collaborators, for simulating convection dominated phenomena (see e.g. [4] , among several other papers published by them before and later on). The originality of our contribution relies on the fact that it not only introduces a reliable scheme for any space dimension, but also exhibits rigorous conditions for it to provide converging sequences of approximations in the sense of the space-time maximum norm.

The main purpose of this work, is to specify procedures for the determination of sets of weights that characterize our explicit scheme. As for the theoretical contribution, we prove that for such choices of the weights the method converges with optimal order in the maximum space-time norm.

An outine of the paper is as follows: In Section 2 we recall the problem to solve together with the type of discretisation corresponding to our explicit scheme; more particularly the weighted manner to deal with the mass matrices on both sides of the discrete equations is described. In Section 3 we further recall the stability results that hold for the method being considered, in the sense of the space and time maximum norm, together with the conditions to be fulfilled by the sequences of meshes and time steps in order to ensure convergence. We proceed in Section 4 by studying in detail some particular choices of the weights associated with the scheme, that satisfy the conditions leading to both stability and convergence. Finally in Section 5 we consider some implementation aspects of the method and give corresponding numerical results.

 

2 The problem to solve and its discretisation

Let us consider a time-dependent convection-diffusion problem described as follows:

Find a scalar valued function u(x, t) defined in × [0, ), Ω being a bounded open subset of N with boundary ∂Ω, N = 1, 2 or 3, such that,

where ut represents the first order derivative of u with respect to t, ν is a positive constant and a is a given solenoidal convective velocity at every time t, assumed to be uniformly bounded in Ω × (0, ). The data f and g are respectively, a given forcing function belonging to L[Ω × (0, )] , and a prescribed value in L[∂Ω × (0, )] . We further assume that u0 L(Ω) and that for every x ∈ Ω g(x, ·) is of bounded variation in (0, ). In (1) ν represents the inverse of the Péclet number.

Without loss of essential aspects, in all the sequel we assume that Ω is an interval if N = 1, a polygon if N = 2 or a polyhedron if N = 3. In so doing we consider a partition h of Ω into N- simplices, with maximum edge length equal to h. We assume that h satisfies the usual compatibility conditions for finite element meshes, and that it belongs to a quasi-uniform family of partitions. We further define a second mesh parameter hmin as the minimum height of all the elements of h if N = 2 or 3 and the minimum length of K h if N = 1.

Let Nh be the number of nodes of h, denoted by Pj, j = 1, 2,..., Nh. We assume that these are numbered in such a manner that the first Ih nodes are located in the interior of Ω and the remaining Nh - Ih nodes are located on ∂Ω. Now for every K h we denote by P1(K) the space of polynomials of degree less than or equal to one defined in K. In so doing we introduce the following spaces or manifolds associated with h:

We further introduce for any function ϕ defined in C0(∂Ω) the following manifold of Vh:

Now let be the field of satisfying (Pj) = u0(Pj) for every vertex Pj of h, and Δt > 0 be a given time step. Defining gn on ∂Ω by gn(·) = g(·, nΔt), fn in Ω by fn(·) = f(·, nΔt) and an in Ω by an(·) = a(·, nΔt), for n = 1, 2,..., idealistically we wish to determine approximations (·) of u(·, nΔt) for n *, by solving the following finite element discrete problem described below, corresponding to a modification of the first order forward Euler scheme.

For n successively equal to 1, 2,..., we wish to determine of the form

where φj is the canonical basis function of Vh associated with the j-th node of h, i.e. Pj, ∈ ℜ being the value of at Pj. We denote by Sj the support of φj and by Πj its measure.

The unknowns for n = 1, 2,..., are recursively determined by the following expressions:

where, setting := an(Pi) and := fn(Pi), the coefficients and are given by:

Coefficient is the well-known lumped mass diagonal matrix (cf. [9] ) given by . The mass matrix coefficients on the right hand side of (2) in turn are defined by a weighted quadrature formula described as follows.

Let Mi be the number of nodes different from Pi lying in the closure of Si, i.e. i, and be such nodes for j = 1, 2,..., Mi with 1 < kj < Nh. Let also be the measure fractions associated with given by:

and be corresponding strictly positive weights satisfying,

Notice that since each N-simplex in Si appears in exactly N measure fractions , we necessarily have:

Now selecting the nodes in i different from Pi, we define,

together with

where is the i-th diagonal coefficient of the standard consistent massmatrix ([9] ) given by

Naturally enough, by definition, = 0 if Pj does not lie in i.

Typically we may choose = for every j and for every node Pi, thereby generating a weighted combination of the lumped mass and the consistent mass matrix (cf. [9] ) on the right hand side of (2), with weights equal to

respectively. However, except for the case of uniform meshes, in principle this is not the choice to make, if one wishes to reach the best results in terms of accuracy.

 

3 Stability and convergence of the scheme

In this Section we recall the stability and convergence results proven in [7] that hold for the above defined weighted mass scheme. In short they state that, provided Δt is chosen conveniently small with respect to the spatial mesh parameter, the scheme (2) is stable in the sense of the maximum norm. Moreover under a suitable condition on the mesh it converges in both space and time in the same sense, as h and Δt go to zero.

First we have to define the following quantities:

• A = | ai(t)|;

ω =

The following theorem proved in [7] states the stability result that holds for scheme (2).

Theorem 1. If Δt fulfills the condition

then the finite element solution sequence given by generated by (2) for n = 1, 2,... satisfies the following stability result for every m , whereby ||F||0, , D denotes the L-norm of a function F defined in an open set D of N, and BV[G] represents the standard norm of a function G(t) having bounded variation for t (0, ):

The above stability result can be refined as follows, in the particular casewhere the partition h is of the acute type (see e.g. [10] ).

Theorem 2. (cf. [7] ). Assume that the partition h is of the acute type (cf. [10] ). Then if Δt satisfies the condition

the finite element solution sequence given by generated by (2) for n = 1, 2,... satisfies the stability condtion (10).

In [7] we give error estimates for the approximations of the solution of (1) generated by (2) under condition (9). In particular we recall here that, provided the weights are suitably chosen, this scheme provides convergent approximations in the maximum norm, as both h and Δt go to zero, under the assumption (11) of Theorem 2. For this purpose we use Sobolev spaces Wm,(D) equipped with the standand norm and seminorm denoted respectively by ||·||m, , D and |·|m, , D, where m is a non negative integer and D is a subset of N (cf [1] ).

As usual a suitable consistency result is needed, which together with the stability results given in this Section leads to convergence. Actually the consistency of our scheme is a consequence of the following lemma:

Lemma 1. (cf. [7] ). Let Pi be a node of h, for i {1, 2,..., Ih}, and be the vector leading from Pi to its neighbor , that is, the j-th node belonging to i, j = 1, 2,..., Mi. Then there exists strictly positive weights satisfying (5) such that

Then we can establish the validity of the following convergence result for scheme (2):

Theorem 2. (cf. [7] ). Let the strictly positive weights , i {1, 2,..., Ih} and j {1, 2,..., Mi}, satisfy (12)-(5). Assume that for a given finite time T > 0 both the solution u of (1) and ut belong to W2, (Ω). Assume also that t [0, T] , a(·, t) [W1, (Ω)] N and f(·, t) W1, (Ω), and (ut)t belongs to L(Ω). Finally let a strictly positive integer kT be defined as the minimum of all integers k such that the quantity Δt : = fulfills the condition (11). Let us further assume that, besides belonging to a quasiuniform family, all the partitions h in use are of the acute type, and the quantity ω is bounded below away from zero independently of h. Then there exists a constant C independent of u, h and Δt, such that the following estimate applies:

 

4 Consistent choices of the weights

In this section we describe two coherent strategies to determine a set of Mi strictly positive weights, for each mesh inner node Pi, that are proven to be bounded below away from zero, independently of the mesh parameter h.

Let us first consider the one-dimensional case. From (12) and (5) it is trivially seen that the pair of weights (, ) associated with inner node Pi is uniquely defined by the equations:

where and are the lengths of the intervals of h having Pi as the right and the left end, respectively. This yields = and = .

Next we switch to the case N > 1. In principle for N = 2 or N = 3 there are infinitely many solutions, except for the case of the least possible value of Mi, i.e. Mi = N + 1, in which the solution is necessarily unique. The two constructions described below allow for the unique determination of a set of weights satisfying (12) and (5), and incidentally they apply even to the particular case where this set is unique.

4.1 A first set of weights

Let := be the set of unit vectors corresponding to the vectors , that is, = /, where is the modulus of . Let also := , and Ei = . Setting Ki = card(Ei), we clearly have Ki > Mi. We number the Mi vectors of Ei in the same manner as the vectors of .

Let us first consider the case where = = Ei. Then for every j1 {1,..., Mi} there is necessarily another j2 {1,..., Mi} such that + = 0. This implies in particular that Mi must be an even number. Hence we may choose the weights in pairs, say (, ), in the same way as in the one-dimensional case (cf. (14)). More specifically, we number the vectors in Ei in such a manner that the first Mi/2 ones form a subset of Ei whose vectors do not have any vector opposite to it in this subset, and from Mi/2 + 1 up to Mi the vectors in the complementary subset consisting of corresponding opposite vectors. In so doing the weights satisfy:

for all pairs (j1, j2) {1,..., Mi/2} ×{Mi/2 + 1,..., Mi}, such that + = 0. Notice that the set of weights determined by solving (15) trivially satisfy both (12) and (5).

Next we assume that (or equivalently, Ei ). In this case we have Ki > Mi, and we number the Ki - Mi vectors in Ei that do not belong to from Mi + 1 up to Ki in an arbitrary order, say , k = Mi + 1,..., Ki. By assumption the intersection of i with the half straight line with origin at Pi and oriented in the sense and direction of for k > Mi, cannot be an edge Pi with j < Mi. Letting the segment Pibe such intersection, where is necessarily a point of the boundary of Si, it follows that it is contained in either a single N-simplex of Si for N = 2 or N = 3, or in a common face of exactly two neighboring tetrahedra for N = 3. In any case, the vector leading from Pi to , for k = Mi + 1,..., Ki, still denoted by , is a non trivial convex combination of either exactly N edge vectors with 1 < j < Mi pertaining to the same N-simplex of Si, or of two such edge vectors pertaining to a common face of two neighboring tetrahedra of Si. Let J be the number of such edge vectors, i.e., either J = N for N = 2 or N = 3, or J = 2 for N = 3 only, and , for l = 1,..., J with 1 < ml < Mi, be the vertices of Si whose convex combination yields for k = Mi + 1,..., Ki. Let us denote the modulus of by , for k = Mi + 1,..., Ki too. In so doing we have:

where the 's for l = 1,..., J, are coefficients of a non trivial convexcombination, that is, 0 < < 1, l = 1,..., J and k = Mi + 1,..., Ki, with

Now we momentarily assign to each node of Si different from Pi, and to each point , j = Mi + 1,..., Ki, the same measure fraction, say

and the weight

Thanks to the fact that by construction, we necessarily have:

together with

as one can easily check.

Now we replace in (17) the 's for k = Mi + 1,..., Ki, by the expressiongiven by (16). Then rearranging the terms in the resulting expression, we establish that relation (12) holds for weights defined in the following manner:

In (19) Ci is a normalizing constant allowing (5) to hold. Notice that, provided the weight increments are all non-negative, the value of Ci is strictly positive. The values of in turn, simply account for the sum of the contributions of for a given j < Mi, to the vectors for k > Mi, expressed by (16), respectively multiplied by , and the corresponding convex combination coefficient. More specifically we have,

where = 0 if does not belong to Si and = for the pertaining convex combination coefficients in (16), that is, all the 's such that ml = j. Notice that by construction (5) holds for the weights defined by (19)-(20).

An important result that holds in connection with the above construction is

Theorem 4.1. The weights generated in the way prescribed in this sub-section are all bounded below by a strictly positive constant ω independent of the mesh step size, whatever the mesh one might consider in the quasi-uniform family of meshes in use.

Proof. The case where Ki = Mi is trivial and hence we give a detailed proof only for the case Ki > Mi.

First of all we note that

Next from (20) and (5) we have,

Since

and from (18) it holds that , taking into account (21) we obtain,

It follows from (22) that for a suitable mesh independent constant c0 we have Ci > c0 i. Indeed Ki < 2Mi and for no mesh of the quasi-uniform family of meshes {h }h under consideration, Mi exceeds the value c- N, where c is a mesh independent constant such that ρ > ch for every h, ρ being the minimum over the elements of h of the radii of the largest inscribed balls in the elements of h for a given h (cf. [3] ).

Finally, since > 0 for all i, j we have,

It immediately follows from (23) that there exists a suitable mesh independent strictly positive constant cN such that

Remark 4.2. (24) clearly indicates that cN may play the role of the parameter ω in the relations (9) and (11). However for very distorted meshes such a value of ω may be largely under-evaluated, and for this reason in practical computations it is advisable to determine this parameter simply as the minimum of all the 's for a given mesh, keeping in mind that, according to (24) such a value is necessarily bounded below away from zero independently of the mesh size.

4.2 A second set of weights

Here we consider again the sets and . The case where = is to be treated in the same manner as above. Hence we may assume that . Distinguishing eventually coincident vectors in and , we put together all the vectors in both sets. Otherwise stated we are dealing with a set of exactly 2Mi vectors not necessarily distinct, say := , where the first Mi vectors are those of and the last Mi vectors are given by , for j = 1,..., Mi. Now let the segment Pibe the intersection of i with the half straight line with origin at Pi and oriented in the sense and direction of for k > Mi, being necessarily a point of the boundary of Si. It follows that Pieither coincides with an edge Pi, where is a vertex of Si, or is contained in either a single N-simplex of Si for N = 2 or N = 3, or in a common face of exactly two neighboring tetrahedra for N = 3. In any case, the vector leading from Pi to , for k = Mi + 1,..., 2Mi, still denoted by , is a convex combination of at most N edge vectors with 1 < j < Mi pertaining to the same N-simplex of Si. Let J be the number of such edge vectors, i.e., J = 1 if coincides with a given vertex of Si, J = 2 for N = 3 only if is contained in a common face of two neighboring tetrahedra of Si, or J = N for N = 2 or N = 3 otherwise. Let , for l = 1,..., J with 1 < ml < Mi, be the vertices of Si whose convex combination yields for k = Mi + 1,..., 2Mi. Here again we denote the modulus of by , for k = Mi + 1,..., 2Mi too. In so doing we have:

where the 's for l = 1,..., J, are coefficients of a convex combination, that is, 0 < < 1, l = 1,..., J and k = Mi + 1,..., 2Mi, with

Now we assign to each point , k = Mi + 1,..., 2Mi the measure fraction , and the weight = , where are provisional weights respectively associated with the vertices of Si different from Pi satisfying = . Then similarly to the case of (15) we necessarily have:

together with

as one can easily check.

Now we replace in (26) the 's for k = Mi + 1,..., 2Mi, by the expression given by (25). Then rearranging the terms in the resulting expression, we establish that relation (12) holds for weights defined in the following manner:

In (28) Ci is a normalizing constant allowing (5) to hold. Notice that, provided the weight increments are all non-negative, the value of Ci is strictly positive. The values of in turn, simply account for the sum of the contributions of for a given j < Mi, to the vectors for k > Mi, expressed by (25), respectively multiplied by , and the corresponding convex combination coefficient. More specifically we have,

where = 0 if does not belong to Si and = for the pertaining convex combination coefficients in (25), that is, all the 's such that ml = j. Notice that by construction (5) holds for the weights defined by (28)-(29).

By arguments in all similar to those in the proof of Theorem 4.1 we can prove,

Theirem 4.3. The weights generated in the way prescribed in this sub-section are all bounded below by a strictly positive constant ω independent of the mesh step size, whatever the mesh one might consider in the quasi-uniform family of meshes in use.

Remark 4.4. The pairs of expressions (19)-(20) and (28)-(29) defining two a priori different sets of weights attempt to take into account three main factors playing a role in the influence of the vector in equation (12). First of all its modulus, since the smaller it is, the larger must be the corresponding weight. Next the associated measure fraction , which goes in the same sense as . Finally the increments , which introduce the necessary adjustment in the weight values, in order to take into account the eventual existence of vectors in Li roughly opposite to , that is, making angles with it close to π. In this sense we may assert that the above described procedure for determining the 's, can be viewed as both close to optimal choices of the weights of our method, and straightforward manners to compute them.

 

5 Implementation and Numerical Aspects

It is not difficult to figure out from the construction of the sets of weights described in the last Section, that both choices give rise to a straightforward implementation of the method. Nevertheless we would like to point out that it is wiser to assemble node by node both the weighted mass matrix and the matrix generated by the discretisation of the convection-diffusion operator, instead of the usual element-by-element procedure. All that is needed for this purpose is a table of node numbers and an associated integer pointer vector having as many components as there are nodes. In the table the series of numbers of the neighbors of each inner node are successively stored, while in the pointer vector the i - th component contains the position in the table corresponding to the first neighbor of the i - th inner node, thereby allowing to uniquely identify all the series of neighbors in the table.

Next we check the performance of the method studied in this paper as compared to a well-established technique to deal with dominant convection in convection-diffusion problems, namely, the least squares formulation. Here the latter is implemented in connection with the Crank-Nicholson scheme for the time-integration, as described in [5] . This comparative study is illustrated by means of some results extracted from [8] , in the framework of a test-problem described below, where uniform meshes are used, thereby allowing the use of weights all equal to 1/(N + 2).

Take Ω to be the unit square (0, 1) ×(0, 1) and a space discretisation based on a uniform L ×L mesh in which every square cell is subdivided into two triangles by means of the diagonal parallel to the line x1 = x2. We take T = 0.1 and ν = 10-k, for k = 2 and k = 5, and a = (1/π; 1/π). For this choice the Péclet number equals 10k/π. An exact solution is considered to be a function with a double boundary layer in the neighborhood of the edges given by x1 = 1 and x2 = 1. More specifically we take

where for z [0, 1] ,

Equation (1) has effectively the above solution, provided the right hand side is given by f(x1, x2, t) = e-t [s(x1) cos(πx2) + s(x2) cos(πx1) ] , and the prescribed initial and boundary values are those of u.

In self explanatory Tables 1 through 4 we summarize the results generated by both methods in the solution of the above test problem, by showing both L2 and L relative errors.

 

 

 

 

 

 

 

 

Observation of Tables 1 and 2 leads to the conclusion that both methods being compared simulate correctly thicker boundary layers, that is, those corresponding to a moderate Péclet number close to 102. Moreover, as one can infer from the above results, our scheme is much less accurate than the least-squares formulation in this case. This is quite natural for the latter is a second order method in the L2-norm, and empirically in the L-norm too, whereas the former is a quasi first order method in the L-norm (cf. [7] ). Nevertheless the least-squares approach failed completely in the case of a thin boundary layer corresponding to Pé roughly equal to 105, as clearly indicated in Tables 3 and 4, while our explicit scheme showed a much better behavior. However even so the errors in the maximum norm are not so small. In this respect it is worthwhile commenting that such results are to be expected. Indeed the maximum error values tend to occur precisely in the interior of the thin boundary layer, which cannot be reached by any numerical method, at least not for the degree of mesh refinement used in the above test.

Acknowledgement. The author is thankful to Carlos Tomei from the Department of Mathematics of PUC-Rio for fruitful discussions.

 

References

[1] R.A. Adams, Sobolev Spaces, Academic Press, N.Y., (1975).         [ Links ]

[2] A.N. Brooks and T.J.R. Hughes, The streamline upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis on the Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering, 32 (1982), 199-259.         [ Links ]

[3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, Noth-Holland,Amsterdam (1978).         [ Links ]

[4] M. Kawahara, N. Takeuchi and Y. Yoshida, Two step finite element methods for Tsunami wave propagation analysis, International Journal of Numerical Methods in Engineering, 12 (1978), 331-351.         [ Links ]

[5] R. Leal Toledo and V. Ruas, Numerical analysis of a least-squares finite element method for the time-dependent advection-diffusion equation, Journal of Computational and Applied Mathematics, 235 (2011), 3615-3631.         [ Links ]

[6] V. Ruas and A.P. Brasil Jr., A stable explicit method for time-dependent convection-diffusion equations, Proc. ICNAAM, Greece, (2007), 480-483.         [ Links ]

[7] V. Ruas, A.P. Brasil Jr. and P. Trales, An explicit scheme for solving convection-diffusion equations, Japan Journal of Industrial and Applied Mathematics, 26 (2009), 65-91.         [ Links ]

[8] V. Ruas, M. Kischinhevsky and R. Leal Toledo, Elementos finitos em formulação mista de mínimos quadrados para a simulação da convecção-difusão em regime transiente, to appear in Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería (January 2013).         [ Links ]

[9] G. Strang and G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs (1973).         [ Links ]

[10] M. Tabata, Uniform convergence of the upwind finite element approximation for semilinear parabolic problems, Journal of Mathematics of the Kyoto University, 18-2 (1978), 327-351.         [ Links ]

 

 

Received: 22/IX/12.
Accepted: 28/IX/12.

 

 

#CAM-698/12.
*Supported by CNPq grant 307996/2008-5.

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License