A weighted mass explicit scheme for convection-diffusion equations

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: [e-mail] / [e-mail] / [e-mail]

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 u_t 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 u⁰∈ 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 h_min 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 N_h be the number of nodes of _h, denoted by P_j, j = 1, 2,..., N_h. We assume that these are numbered in such a manner that the first I_h nodes are located in the interior of Ω and the remaining N_h - I_h nodes are located on ∂Ω. Now for every K ∈ _h we denote by P₁(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 C⁰(∂Ω) the following manifold of V_h:

Now let be the field of satisfying (P_j) = u⁰(P_j) for every vertex P_j of _h, and Δt > 0 be a given time step. Defining gⁿ on ∂Ω by gⁿ(·) = g(·, nΔt), fⁿ in Ω by fⁿ(·) = f(·, nΔt) and aⁿ in Ω by aⁿ(·) = 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 V_h associated with the j-th node of _h, i.e. P_j, ∈ ℜ being the value of at P_j. We denote by S_j the support of φ_j and by Π_j its measure.

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

where, setting := aⁿ(P_i) and := fⁿ(P_i), 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 M_i be the number of nodes different from P_i lying in the closure of S_i, i.e. _i, and be such nodes for j = 1, 2,..., M_i with 1 < k_j < N_h. Let also be the measure fractions associated with given by:

and be corresponding strictly positive weights satisfying,

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

Now selecting the nodes in _i different from P_i, 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 P_j does not lie in _i.

Typically we may choose = for every j and for every node P_i, 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 = | a_i(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 W^m,∞(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 P_i be a node of _h, for i ∈ {1, 2,..., I_h}, and

be the vector leading from P_i to its neighbor , that is, the j-th node belonging to _i, j = 1, 2,..., M_i. 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,..., I_h} and ∀j ∈ {1, 2,..., M_i}, satisfy (12)-(5). Assume that for a given finite time T > 0 both the solution u of (1) and u_t belong to W^{2, ∞}(Ω). Assume also that ∀t ∈ [0, T] , a(·, t) ∈ [W^{1, ∞}(Ω)] ^N and f(·, t) ∈ W^{1, ∞}(Ω), and (u_t)_t belongs to L^∞(Ω). Finally let a strictly positive integer k_T 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 M_i strictly positive weights, for each mesh inner node P_i, 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 P_i is uniquely defined by the equations:

where and are the lengths of the intervals of _h having P_i 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 M_i, i.e. M_i = 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 E_i = ∪ . Setting K_i = card(E_i), we clearly have K_i > M_i. We number the M_i vectors of E_i ∩ in the same manner as the vectors of .

Let us first consider the case where = = E_i. Then for every j₁∈ {1,..., M_i} there is necessarily another j₂∈ {1,..., M_i} such that + = 0. This implies in particular that M_i 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 E_i in such a manner that the first M_i/2 ones form a subset of E_i whose vectors do not have any vector opposite to it in this subset, and from M_i/2 + 1 up to M_i the vectors in the complementary subset consisting of corresponding opposite vectors. In so doing the weights satisfy:

for all pairs (j₁, j₂) ∈ {1,..., M_i/2} ×{M_i/2 + 1,..., M_i}, 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, E_i≠ ). In this case we have K_i > M_i, and we number the K_i - M_i vectors in E_i that do not belong to from M_i + 1 up to K_i in an arbitrary order, say , k = M_i + 1,..., K_i. By assumption the intersection of _i with the half straight line with origin at P_i and oriented in the sense and direction of for k > M_i, cannot be an edge P_i with j < M_i. Letting the segment P_ibe such intersection, where is necessarily a point of the boundary of S_i, it follows that it is contained in either a single N-simplex of S_i 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 P_i to , for k = M_i + 1,..., K_i, still denoted by , is a non trivial convex combination of either exactly N edge vectors with 1 < j < M_i pertaining to the same N-simplex of S_i, or of two such edge vectors pertaining to a common face of two neighboring tetrahedra of S_i. 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 < m_l < M_i, be the vertices of S_i whose convex combination yields for k = M_i + 1,..., K_i. Let us denote the modulus of by , for k = M_i + 1,..., K_i 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 = M_i + 1,..., K_i, with

Now we momentarily assign to each node of S_i different from P_i, and to each point , j = M_i + 1,..., K_i, 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 = M_i + 1,..., K_i, 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) C_i is a normalizing constant allowing (5) to hold. Notice that, provided the weight increments are all non-negative, the value of C_i is strictly positive. The values of in turn, simply account for the sum of the contributions of for a given j < M_i, to the vectors for k > M_i, expressed by (16), respectively multiplied by , and the corresponding convex combination coefficient. More specifically we have,

where = 0 if does not belong to S_i ∩ and = for the pertaining convex combination coefficients in (16), that is, all the 's such that m_l = 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 K_i = M_i is trivial and hence we give a detailed proof only for the case K_i > M_i.

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 c₀ we have C_i > c₀∀i. Indeed K_i < 2M_i and for no mesh of the quasi-uniform family of meshes {_h }_h under consideration, M_i 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 c_N such that

Remark 4.2. (24) clearly indicates that c_N 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 2M_i vectors not necessarily distinct, say := , where the first M_i vectors are those of and the last M_i vectors are given by , for j = 1,..., M_i. Now let the segment P_ibe the intersection of _i with the half straight line with origin at P_i and oriented in the sense and direction of for k > M_i, being necessarily a point of the boundary of S_i. It follows that P_ieither coincides with an edge P_i, where is a vertex of S_i, or is contained in either a single N-simplex of S_i 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 P_i to , for k = M_i + 1,..., 2M_i, still denoted by , is a convex combination of at most N edge vectors with 1 < j < M_i pertaining to the same N-simplex of S_i. Let J be the number of such edge vectors, i.e., J = 1 if coincides with a given vertex of S_i, J = 2 for N = 3 only if is contained in a common face of two neighboring tetrahedra of S_i, or J = N for N = 2 or N = 3 otherwise. Let , for l = 1,..., J with 1 < m_l < M_i, be the vertices of S_i whose convex combination yields for k = M_i + 1,..., 2M_i. Here again we denote the modulus of by , for k = M_i + 1,..., 2M_i 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 = M_i + 1,..., 2M_i, with

Now we assign to each point , k = M_i + 1,..., 2M_i the measure fraction , and the weight = , where are provisional weights respectively associated with the vertices of S_i different from P_i 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 = M_i + 1,..., 2M_i, 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) C_i is a normalizing constant allowing (5) to hold. Notice that, provided the weight increments are all non-negative, the value of C_i is strictly positive. The values of in turn, simply account for the sum of the contributions of for a given j < M_i, to the vectors for k > M_i, expressed by (25), respectively multiplied by , and the corresponding convex combination coefficient. More specifically we have,

where = 0 if does not belong to S_i ∩ and = for the pertaining convex combination coefficients in (25), that is, all the 's such that m_l = 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 L_i 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 x₁ = x₂. 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 10^k/π. An exact solution is considered to be a function with a double boundary layer in the neighborhood of the edges given by x₁ = 1 and x₂ = 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(x₁, x₂, t) = e^-t [s(x₁) cos(πx₂) + s(x₂) cos(πx₁) ] , and the prescribed initial and boundary values are those of u.

In self explanatory Tables 1 through ⁴ we summarize the results generated by both methods in the solution of the above test problem, by showing both L² 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 10². 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 L²-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 10⁵, as clearly indicated in Tables 3 and ⁴, 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.

Received: 22/IX/12.

Accepted: 28/IX/12.

#CAM-698/12.

• [1] R.A. Adams, Sobolev Spaces, Academic Press, N.Y., (1975).
• [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.
• [3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, Noth-Holland,Amsterdam (1978).
• [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.
• [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.
• [6] V. Ruas and A.P. Brasil Jr., A stable explicit method for time-dependent convection-diffusion equations, Proc. ICNAAM, Greece, (2007), 480-483.
• [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.
• [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).
• [9] G. Strang and G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs (1973).
• [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.

*

Supported by CNPq grant 307996/2008-5.

Publication Dates

• Publication in this collection
  28 Nov 2012
• Date of issue
  2012

History

• Received
  22 Sept 2012
• Accepted
  28 Sept 2012