Exactly equilibrated fields, can they be efficiently used for a posteriori error estimation? Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Том 148, кн. 4 Физико-математические пауки 2006

UDK 519.63

EXACTLY EQUILIBRATED FIELDS, CAN THEY BE EFFICIENTLY USED FOR A POSTERIORI ERROR ESTIMATION?

I.E. Anufriev, V.G. Korneev, V.S. Kostyle.v

Abstract

The answer given in the paper to the question in the title: yes, they can. We advocate the approach to the a posteriori error estimation, which can be called “classical”, and for the theory elasticity problems stems from the Lagrange and Castigliano variational principles. In it, the energy of the error of an approximate solution, satisfying geometrical restrictions, is estimated by the energy of the difference of the stress tensor corresponding to the approximate solution and any stress tensor, satisfying the equations of equilibrium. Notwithstanding a popular point of view that the construction of equilibrated stress fields requires considerable computational effort, we show that it can be practically always done for the number of arithmetic operations, which is asymptotically optimal. We derive also new general a posteriori estimates, in which equilibrated fields are replaced by arbitrary fields of fiuxes/stresses. Numerical experiments show that our a posteriori error estimators provide very good coefficients of effectiveness, which in many cases can be convergent to the unity. At the same time they have linear complexity and are robust.

1. Introduction

Publications on a posteriori error estimates for approximate, e.g., numerical solutions of partial differential equations are numerous. The earliest a posteriori error estimates were apparently known in mechanics from the time when the Lagrange and Castigliano principles, which, from the mathematical point of view, provide the primal and dual formulations in theory elasticity, obtained a mature form. Such estimates are deduced from the fact that approximate solutions obtained on the basis of these principles approach the exact solution in the energy sense from the opposite directions, and. namely, from above and from below, respectively. Let u be the vector of exact displacements of a linearly elastic body, a = a (u) is the corresponding stress tensor, | • |U and | • \a are the potential energy norms expressed in displacements and stresses. If u is an arbitrary displacement vector of finite energy, satisfying geometric boundary conditions, and t is an arbitrary stress tensor of finite energy, satisfying the equilibrium equations (including the boundary conditions in stresses), then the classical a posteriori estimate is

|u - U|u < | a (U) - t \a. (1.1)

and may be found, e.g., in the Mikhlin’s book [33]. In spite of its simple form and enormous amount of publications on a posteriori estimates, the authors were unable to find references where it had been directly used in practical FE applications.

During several last decades, a few groups of a posteriori error estimation techniques have been developed and, first of all, so called residual-based techniques found in Babuska I and Reinbolt [5, 6], Verfiirth [51] Stewart and Hudges [48] and more recent publications. Among them there are distinguished the explicit residual method,

see Bernardi and Girault [9]. Carstensen [11] and the related paper by Clement [12] on special interpolation, and the implicit and equilibrated residual methods, for which we refer to Babuska I and Reinbolt [5], Kelly [18], Ladevese and Leguillon [31], and Ainsworth, Demkowicz and Kim [1]. Many papers are dedicated to obtaining of indicators of the error, which do not necessarily bound it, but are approximately proportional to the error, and, therefore, can be used for the mesh refinement in adaptive computations. A pioneering paper Zienkiewicz and Zhu [55], see also Ainsworth and Oden [2], commenced the group of such techniques widely used in applications and employing superconvergence properties of finite element solutions. Many contributions have been related to a posteriori error estimation for specific problems of mathematical physics. The development of a posteriori estimation techniques in recent decades, as well as

ii

Oden [2] Babuska and Strobolis [7] and Neittaanmaki and Repin [35].

The main idea of some of the mentioned approaches to the error estimation is to use fields of stresses, which can be obtained from the FE solution and at the same time are most close to the exact equilibrated fields (i.e., representing the exact solution of the problem). An example is the equilibrated residual method, which now gains more attention as a method which allow to obtain reliable bounds, often without solving some global systems of algebraic equations, see, e.g., Ainsworth, Demkowicz and Kim [1], Luce and Wohlniuth [32], Vejchodsky [50] and Braess and Schoberl [10]. However, it is also true that the purpose of most authors is to outflank construction of exactly equilibrated fields at all1. One way of obtaining equilibrated fields2, which approximate equilibrated fields of the exact solution of the primal problem, is approximate solution of the dual problem, which in the theory of elasticity is expressed by the Castigliano principle of virtual equilibrated states. As a rule a motivation for avoiding the use of equilibrated fields is that the solution of the dual problem or other ways of finding such fields are computationally too expensive.

The purpose of this paper is to illustrate that in many cases the estimate (1.1) can be directly used as an efficient and cheap error estimator. This is for the reason that indeed equilibrated fields are not difficult to find in a variety of ways. As one of the options, as we will see not the most efficient in many cases, the numerical solution of the dual problem can be considered. For advocating this option, the following fact is important: numerical solutions of the discretizations of the primal and dual problems, having the same (in the order) accuracy in the energy norm, can be found for the same (in the order) computational cost. More over, under some conditions, the discretization of the dual problem may be obtained in such a way that its matrix will coincide with the FE matrix for the primal problem up to the boundary conditions. Therefore, practically the same solver can be used for solution of the both discrete problems.

From the above discussion, one concludes that the option of solving the dual problem for evaluating the equilibrated fields deserves examination. Suppose that discretizations of the same order of accuracy are used for the primal and dual problems. In general, one can expect that the efficiency coefficient will converge to the unity at h ^ 0, if the error estimator is super-convergent. In the practice, such convergence was observed for a number of alternative a posteriori error estimators, considered in this paper and papers of other authors, see, e.g., Luce and Wohlniuth [32]. In our numerical experiments with the use of classical error estimate of the type (1.1) and of the equilibrated fields, obtained by solution of the dual problem with the same accuracy, the observed effectiveness coefficient remained close to 1.25. Further improvement of the efficiency coefficient is

1In general discussions we use the theory elasticity problem for a model without special remarks.

2In the paper the both terms equilibrated and exactly equilibrated field imply that a field satisfies equilibrium equations exactly in classical or generalized sense.

possible, e.g., if we solve the dual problem on a denser, then for the primal problem, mesh, which results in a greater cost of an error estimator, than the cost of the solution to be validated.

At least not less promising approach can be referred as direct evaluation of the equilibrated fields. It is based on the fact that the equilibrium equations (in stresses) are under-determined. For instance, in the theory of elasticity, the symmetric stress tensor

& = iaki}k i=1 with six stresses for the entries satisfies three equilibrium equations. Therefore, in order to satisfy the equilibrium equations it is sufficient to perform two steps:

1) to specify three stresses, say shear stresses

&ki, k = l, (1.2)

by arbitrary sufficiently smooth functions and

2) to find the rest stresses from the equilibrium, equations by evaluating 1-d integrals. The presence of the boundary conditions for stresses does not make this procedure significantly more difficult. When this procedure is used in the a posteriori estimator. e.g., for the FE solution, the stresses (1.2) are found from the FE solution with the accurate use of its superconvergence properties.

The approach of direct evaluation of the equilibrated fields allows us not only to design cheap algorithms for evaluating a posteriori error estimators, often producing convergent the efficiency coefficients, but to come to a new type of general a posteriori estimates. For simplicity, let us consider the Dirichlet problem for the Poisson equation in the rectangle Q = (a1,b1) x (a2, b2)

-Au = f (x), x = (xi,x2) G Q, u|dQ = 0. (1-3)

and any function v from H1(Q) = {w G H 1(Q) : u|dfi = 0}, which is considered as an “approximation” for u. Then for any e > 0 the error estimate can be written in the form

l|V(v - u)||0 , n < (1+ e)||Vv - y||0 ,n+

xk

+ + J ak(f ~ V-y)(?7fc,x3_fc)d?7fc

ak (x3- k)

(1.4)

y = (y1, y2)T is any sufficiently smooth vector-function, functions ak = ak(x) satisfy a1 + a2 = 1, and % convention V • y = divy. If f G L2(Q) and ak G C(Q) it is natural to consider y G H(Q, div) = {z : V-z G L2(Q)}. If v = ufem is the FE solution, then y may be obtained % the averaging of the derivatives dufem/dxk at the nodes and interpolation. In our numerical experiments, algorithms of that type almost always produced good and convergent effectiveness coefficients. In the algorithms presented in the paper we implement a variety of techniques of direct evaluation of the exactly-equilibrated fields for their use in the a posteriori estimators.

Clearly, the estimate (1.4) fetches additional opportunities in comparison with the known estimates of a close appearance, in which the second term in the right part is. e-9--, If - V • y||-1in- Numerical evaluation of this negative norm is not at all easy, whereas its replacement in (1.4) can be often computed for a number of arithmetic operations proportional to the number of unknowns. In particular, this is true for the FE discretizations by means of the orthogonal grids. Another example is Galerkin methods with the coordinate functions specified by analytical expressions in the whole domain. In

the paper we discuss also the procedures, which can provide the optimal computational cost in more general situations.

The equilibrium equations in the theories of thin shells and shells of moderate thickness. see. e.g., Gol’denweizer [17] Novozhilov [36. 37] and Reissner [44]. are written in terms of internal forces, i.g., shear forces and bending and twisting moments. These equations are more complicated, than the equilibrium equations in the theory of elasticity. However, a quite similar approach can be implemented for obtaining the exactly-equilibrated functions of internal forces. This approach originates from the papers devoted to numerical methods for solving the bending thin plate and shell problems on the basis of the Castigliano principle and the method of splitting of the thin plate and shell partial differential equations, studied by Rozin [45. 46]. Korneev/Rozin [29. 30]. Korneev [19 21. 25. 27].

Considerable part of the algorithms contained in the paper were tested numerically. We present the graphs and tables of numbers illustrating the dependence of effectiveness indices and arithmetical work on numbers of unknowns. Obviously, all tested algorithms are optimal in the computational work. Additionally to this, the practical computational costs of the a posteriori error estimators and of the optimal multigrid solvers for the primal problem were compared. As a rule the latter exceeded the former in two times at the least.

The paper is arranged as follows. In Section 2. we consider a posteriori estimators for the Poisson equation in the unite square and arbitrary Lipshits continuous domain with different boundary conditions. The special case of the differential operator with the discontinuous coefficient in the main term is treated in Subsection 2.3. It is shown that the algorithms of a posteriori estimators are easily adjusted to this case. The results of numerical tests for this case, presented in Subsection 5.2. show that discontinuity practically does not affect the efficiency of the a posteriori error estimator. Section 3 deals with the a posteriori error estimators for the plane elasticity problem. All a posteriori error estimators of Sections 2 and 3 are based on the direct evaluation of the balanced fluxes and equilibrated stress tensors, i.e., without solving any systems of algebraic equations. In Subsection 4.1. we consider an alternative approach based on solution of the dual problem equivalent to the Castigliano principle of virtual complementary-work. We show that it is possible to chose a basis in the space of the self-equilibrated tensors of stresses in such a way that the system of algebraic equations will possess practically the same properties as the FE system for the primal problem. Throughout the paper alongside with the general algorithms, we present the algorithms, which we tested numerically. Results of numerical experiments are discussed in Section 5.

Throughout the paper, we use the notations listed below.

P^d Q p are the spaces of polynomials of the total order p and of the order p in each variable, e1 = (1,0), e2 = (0,1), d is the dimension.

L2(Q) is the space [L2(Q)]d with the norm || • || = || • ||2iq and the same notation is used for the norm in L2(Q),

| • |k,n > II • ||k,n stand for the semi-norm and the norm in the Sobolev space Hk(Q),

Dqx v := d|x|v/dxf dxf ,...,dxf, q = (q1,q2, . .. ,qd), qk > 0, |q| = 91+92+------+qd,

k

where

H/1(Q) := {v G H 1(Q) : v|dQ = 0} is the subspace of functions from H 1(Q) vanishing on the boundary dQ.

Wo use also the abbreviations: a.o. arithmetic operations. FE finite element. In the vectors of the space variables x = (x1; x2) or x = (x1? x2, x3), sometimes we interchange the positions of variables and write x = (xk,x3-k) and x = (xk, xk+1,xk+2), assuming in the latter case that indices k + l are taking values modulo 3.

2. Poisson equation

2.1. An outline of the approach. In this section, we illustrate basics of the approach of the direct evaluation of equilibrated fields on a simple model problem. Let us consider a boundary value problem for the Poisson equation in the unite square Q = (0,1) x (0,1) with the mixed boundary conditions

where

—Am = f (x), x = (x1, x2) G Q, dQ = Td U Tn, w|td = 0, dw/dv|rN = 0.

rd = {x| x1 G (0,1] , x2 = 1} u {(x1,x2)| x1 = 1 , x2 G (0,1]},

(2.1)

(2.2)

rN = {x| x1 G [0,1) , x2 = 0} U {(x1, x2)| x1 = 0 , x2 G [0, 1)} ,

and v is the distance from the boundary along the outward normal. The generalized formulation of this boundary value problem reads

a(u,v) = (f,v)n , Vv G V(Q), (2.3)

where V(Q) = {v G H 1(Q) : v|rD = 0}

»(„,») = / V'v -Vwdx, = j«»*.

n n

Let V (Q) be the finite element space of the piece wise bilinear functions on the uniform square grid of size h = 1/n, n > 1, with the nodes x(i) = h(i1,i2), ik = = 0,1,... ,n, and Vo(Q) is the subspace of functions from V(Q), vanishing on rD. By ufem is denoted the finite element solution belonging to V0(Q) and satisfying the identity

a(wfem,v) = (f,^)n , Vw GVo(Q). (2.4)

In order to be able to efficiently implement the a posteriori estimate (1.1). it is necessary

with the use of the obtained FE solution ufem to construct the vector valued function

t = (t1, t2)T , which obeys the two conditions:

a) it satisfies the balance differential equation

-V-1 = f, (2.5)

and the boundary conditions

V • t|rN =0, (2.6)

where v is a unite vector normal to dQ, and

ft) is as much close as possible to the gradient Vm of the exact solution.

We shall use notations Qf, Q0 for the sets of functions satisfying (2.5), (2.6) with the given f and f = 0, respectively, from which Q0 is clearly a liner space. The balance law

(2.5) models equilibrium equations in the case of the theory elasticity boundary value Qf

Q0 - self-balanced or self-equilibrated. The sets Qf, Q0 can be defined constructively

by moans of tho splitting technique, which was introduced by Korneev and Rozin [27. 29. 45] at developing numerical methods for solving problems of the theory of elastic thin plates and shells in solid mechanics on the basis of the Castigliano principle. For the problem under consideration constructive definition of the set Qf of the “equilibrated” vectors is quite simple. If q is an arbitrary sufficiently smooth function, then the vector t = (ti,t2)T with the components

satisfies equation (2.5) and boundary conditions (2.6). Clearly, Qf = Qo + tf, where tf any element of Qf, and Q0 is defined % (2.7) with f = 0. For all q G L2(Q) one comes to the space Q0 with the norm |t|CT = ||t||, where || • || stands for the [L2(^)j2 norm.

The estimate (1.1) takes the form

Obviously, a better approximation of du/dx1 by t1 (e.g., with the use of values of the gradient of the finite element solution at superconvergence points) will result in a better

q

the second derivative d2u/dx1.

Taking for sufficiently smooth functions satisfying a1 + a2 = 1, one can use more “symmetric” formulas instead of (2.7):

They can provide more accurate a posteriori estimates, especially with a good choice of functions but require more a.o.

If the approach, presented above, is used for a posteriori estimation, then different boundary conditions should be given attention. Suppose, w is an approximate twice differentiable solution of (2.1), (2.2), e.g., obtained by the Galorkin method or any other function from V(Q). For instance, we can set q = d2w/dx1 and come to the expressions

(2.7)

0

0

||V(u - Wfem)|| < ||Vufem - t||, V t G Qf.

(2.8)

(2.9)

Xi

0

(2.10)

Similarly, we can proceed from setting q = d2u/dx2, coming to the a posteriori estimates

Xk

||V(m-m)||“ </| J Q^^k,x3-k)d^k

X 3-k

f d2u

:+ J (f+Q^)(xk,&-k)dZ3-

dx, k = 1, 2. (2.11)

Taking into account boundary conditions at rN and triangular inequality, one comes from (2.11) to

l|V(u -SOIK ]T k = 1,2

du i

Xk

I Xk=0

+

(0,1)

(f + Au)(x1,£2) d&

(2.12)

t Qf t

Xi

t1 = J (a1d2u/dx1 — a2(f + d2w/dx2)) (C1;x2) d^, 0

X2

t2 = J («2d2M/dx2 - «1(f + 52w/dx2)) (x1,^2) d^2.

(2.13)

The corresponding a posteriori error estimates are ||V(u - S)|| <

Xk

I f du f ( d2M d2M A

“I VkM~ ^Lk)

1/2

dx V <

<

k=1,2

du

— +

xk Xk=0 (0,1)

Xk

J «3-fc (f + Aw)(£fc,x3_fc) d^fc

(2.14)

which for ak = 0.5 is invariant with respect to xk, k =1,2. It easy to see, that adding and subtracting «3_kd2M/dx| inside round brackets, taking into account boundary conditions at rN and triangular inequality, we obtain the same estimate (2.12).

In the case of the Dirichlot boundary value problem

a(u,v) = (f,v)n u, Vv G V(fi) = id1(0), (2.15)

the estimates (2.12),(2.14) take especially simple forms. Instead of the latter we have

II V(u - S)|| < ]T

k=1,2

Xk

/“3-k(f+Au)(£k 'xs-k) d£t

(2.16)

2

Let us consider for simplicity the case a2 = 0. Since in this case no boundary conditions are imposed on the equilibrated fluxes, we can set

ti = du/dxi,

X2

t2 = dw/dx2(xi, 0) - J(f + d2w/dx2)(xi, ^2)d^2,

0

and. therefore, du

~-----#1=0,

dxi

X2 _ X2

du du du f ( d2u \ , /■

l^2~t2 = ~^2~ ~dx2 J \ + J (Xl’^2) = J ( ~~ Au)(xi’&)^2’

0 1 0

completing the proof. In the case of an arbitrary sufficiently smooth domain, the proof is similar.

Lemma 2.1. Let Q be Lipshits continuous domain, f G L2(Q), u - solution of (2.15), and u be any function in H2(Q) satisfying boundary condition u|dfi = 0. Then u - uu

2.2. Examples of algorithms for finite element solutions of Poisson equation. Solutions obtained by FE methods compatible in C, which are primarily used in practice for second order elliptic equations, do not have second derivatives. Basically, three ways to outflank this obstacle can be distinguished. All of them start from the procedure of constructing some smooth approximation of the second, or first derivatives of FE solution, or the FE solution itself. In what follows, this procedure is termed smoothing procedure. After smoothing procedure has been applied, we proceed in one of the ways described in Subsection 2.1. The distinctions between three types of a posteriori error estimation algorithms for our model problem can be illustrated on example of

f

of only one of the fluxes. Since this flux is uniquely defined by the balance equation and the boundary condition on rN, it is sufficient to point out the way of evaluation of one flux, which is calculated first. Briefly, three types of such a posteriori error estimation algorithms are the following:

a) Calculate second derivative of the FE solution along one of the axes approximately with the use of finite differences at some set of discrete points (e.g., FE nodes). Define q (e.g., as a function of the FE space V (Q)) by interpolation of the calculated approximate values of the second derivative. Define the corresponding flux by 1-d into-

q

flux, given by the boundary conditions on rN.

b) Calculate the first derivative of the FE solution in one of the directions at the nodal points, e.g. by averaging. Define the tentative flux in the chosen direction as the FE function, which belongs to V (Q) and takes at the nodes calculated values. Define the flux by adjustment of the tentative flux to the boundary condition on rN.

c) Construct twice differentiable approximation of the FE solution. Define the tentative flux in one of the axes as the first derivative of the smoothed FE solution. Define the flux by adjustment of the tentative flux to the boundary condition on rN.

Indeed each of a), b) and c) allows to define two equilibrated fluxes t(k), k = 1,2, corresponding to the direction , the flux along which is defined first. The flux for the a posteriori estimator can be defined as t = a1t(1) + a2t(2).

Apart from “symmetric” versions, other variations of the outlined algorithms are numerous. For instance, in the smoothing procedures of a) and b) some other finite dimensional functional spaces W (Q) can be used instead of the basic FE space V (Q).

qQ subdomains Qj and in each define q by the list squares method with the use of polynomials of some specific order pj. In the vicinity of singularities in the exact solution

q

solution can be implemented. In the latter case it is more appropriate to use the term preprocessed FE solution. It is also worth adding that in different subdomains one can use algorithms of different types, i.e., a), b) or c). for obtaining equilibrated fluxes.

In this subsection, we present a few simple examples of the outlined algorithms for the Poisson equation in the unit square and arbitrary sufficiently smooth domain with different boundary conditions. For the problems in the unit square, the FE space

V (Q) is the space of continuous piece wise bilinear functions on the square mesh of size h = 1/n. For the nodes of this mesh we use notation x(i) = hi, i = (i1, i2). We start from the algorithm of the type c) for the problem (2.1) (2.2).

Algorithm 2.1.

Step 1. For each node g QQ calculate the value of the mesh function = = (v2i1)"1ij2=cM which is the ffnite-difference approximation of the second derivative d2u/dx12(x(i)). For internal nodes of horizontal mesh lines x2 = hi2,

i2 =0,1,..., n, use

(i) ufem(h(i1 + 1^ hi2) 2ufem(hi1, hi2) + ufem(h(i1 1^ hi2) ■

*’2,1 = -----------------------------p-------------------------------> *1 = 1,..., n-l.

For the nodes (0, i2) on the axis x1 = 0 set

v2i1 = d2uo/dx12(0, hi2),

where u0 is the 3-d order interpolation polynomial of x1 over the values ufem(x(i)) for i1 = 0,1,2, Mid dufem/dx1(0, hi2). For the nodes (n, i2) calculate

v2i)1 = d2u1/dx12(1, hi2),

where u1 is the 3-d order Lagrange interpolation polynomial of x1 over the values ufem(x(i)) for i1 = n — 3, n — 2, n — 1, n

Step 2. Define I2,1(x) c V (Q) as the piece wise bilinear interpolation of t

Us 1

t1(x) = y ^2,1 (^1, X2 ) d^1,

0 X2 (2-18)

#2(x) = — J(f ^2) + ^2,1 (x1, ^2)) d^2.

0

Step 4. Evaluate the estimator n := ||Vufem — t||2.

Remark 2.1. Formulas (2.18) correspond to (2.7), (2.10). Since I2,2 ~ d2u/dx12 can be calculated in a similar way, a more ‘'symmetric” formulas, corresponding to (2.9),

(2.13). can bo used:

Xi Xi

tl = J [^1/2,1 — a2(f + /2,2)](Clj X2) = J [/2,1 — «2(7 + /2,1 + I"2,2)](£l 5X2)

0 0 X2 X2

t2 = J\<.»h, — M/ + d6 = /ft.2 — »1(/ + /2, + /2,2)](X1.^2) '«2-

00

(2.19)

In Step 1. we used finite-difference approximations of the second derivative d2u/dx2, which provide the same order h2 of approximation for a sufficiently smooth u. In Step 3, the way of the evaluation of integrals of /^d / may be different. In particular, in many cases the analytical integration can be performed. In general, procedures, used in Steps 3, depend on the way of evaluating the norm || Vufem — 1||2 in Step 4. For instance, for each finite element := h(i1 — 1, i1) x h(i2 — 1, i2), we can use quadratures exact for polynomials of some order p* > 1. Then it is necessary to evaluate tk only at the quadrature nodes, and for doing this other quadratures can be used. The type of the quadratures, used element wise for evaluating t2 and the norm ||Vufem — 11|2 may depend on the local smoothness of /, if the integral of / is not evaluated exactly. Bounds for quadrature errors will enter the resulting error bound for

ufem /

not elaborate on these subjects in the present paper aimed to illustrate main features of the approach.

In the case rN = 0, i.e. when only the Dirichlet boundary conditions are imposed

1

balance equation (2.5) and not subjected any boundary conditions. But when boundary conditions are different from the ones considered above in (2.1) and rN = 0 then some

1

rD = {(X1,X2)| X1 e [0, 1], X2 = 1} , Tn = Tn,1 U Tn,2 U rN,3.

where

rN,1 = {(X1,X2)| X2 e [0, 1], X1 =0},

rN,2 = {(X1,X2)| X1 e [0, 1], X2 =0},

rN,3 = {(X1,X2)| X2 e [0, 1], X1 = 1} ,

and the boundary conditions are

ulr.

0,

u|rf

(2.20)

where g G C(rn) and

| g1 (x2 ): on Tn,1

g = < g2 (x1 ): on r N,2

[ g3(x2): on r N,3

(2.21)

1

and the boundary condition (2.6). The steps 1 and 2 may be the same as in the algorithm adduced above, but they produce only a tentative flux, which should be adjusted to the

boundary condition (2.6). Lot

Xi

t1(x) = J(/2,1 + r)(^1,X2) d^1 + g1 (x2),

0 X2 (2-22)

t2(x) = — y (/ + /2,1 + r)(x1 ,£2) d^2 + g2 (x1)

where r(x) is chosen with the purpose to fulfill the boundary condition (2.6) on rN,3. It is sufficient to take r depending only on x2

1

r(x1,x2)= g3(x2) — g1 (x2) — y /2,1 (£1, X2) d£1 . (2.23)

Hence

X1 1

t1(x)^y* /2,1 (£1, X2) d^1 — X1 J /2,1 (£1, X2) d^1 + (1 — x1)g1(x2)+ X1g3(x2),

00

X2 X2 1

t2(x) = — J(/ + /2,1)(x1,6) d£2 +J J /2,1(^1,C2)) d^1d^2+ (2.24)

0 0 0 X2

/ (g1 — g3)(£2) d£2 + g2(x1).

In the case of the homogeneous Neumann boundary condition the expressions (2.24) on rN takes a simpler form

Xi 1

t1 = / /2,1 (£1, X2) d^1 — X1 J /2,1 (£1, X2) d^1,

00

X2 X2 1

^ = —J(/ + /2,1)(x1,£2) d£2 +J J /2,1(£1,£2)) d£1d£2.

0 0 0

The described approach is easily realized in a much more general situation. Suppose. Q is the domain occupied % the arbitrary triangulation Sh with the triangles Tr, r = 1, 2,..., R, V (Q) is the space of the continuous piece wise linear functions and V (Q) = {v e V (Q) : v|dfi=0}. We will turn to the algorithm of the type b) for the problem with the Dirichlet boundary condition. Namely, we consider the problem

a(u, v) = (/, v)n u, V v e V(Q) = V1(Q), (2.26)

ufem

a(ufem, v) = (/,v)n ufem, Vv e V (Q). (2.27)

Let us assume for simplicity that each line xk = const crosses Q not more than in two points, rk,_ is the part of the boundary containing the points of such pairs, having

lesser coordinate xk, xk = ak(x3_k) is the equation of rk,_, x3_k = a3_k,63_k are

the coordinates of the ends of I\_. We also use notation Bt = {r : rr G Q, x^ G rr} for the set of the numbers of finite elements, having x(i) for a vertex, and notation

V (rk,_) for the space of the traces of functions from V (Q) on rk,_. Among simplest is the algorithm, which is not invariant with respect to xk, k = 1,2, and is based on the averaging and the procedure reflected in Remark 2.1.

Algorithm 2.2.

Step 1. For each node calculate the average

rEBi

(r)

where u(erm = ufem

is the restriction of the FE solution to Tr.

Step 2. Define the interpolation /(wfem,i) G V(Q) satisfying

I{U fem,i)(xW) =4L,1 Vi«Gn.

Step 3. For each node g p2,- calculate the average

-« _ y- dutl( (i))

fem,2 - 2^ dxo ^% h rEBi 2

and for x G T2j_ define the piece wise linear continuous interpolation irfc,_ («fem,2) satisfying

/(«fem,2 ) (*(i) ) = «11,2 V £ Fa, - •

Step 4. For components of an equilibrated vector t = (t1,t2)T set

#1 = /(Wfem,l), ^2 | p2 _ = <f>2{xi)) := Irk'- («fem,2)|r2

and evaluate

X2

h(x) = h(xi, ch(xi)) ~ J (xi>6)^2-

«2

Step 5. Evaluate the bound ||Vufem — t||2.

The algorithm, based on averaging, can be made invariant with respect to xk, k = = 1, 2

Algorithm 2.3.

Step 1. For each node and k = 1,2 calculate the averages

o (r)

-« _ V- dufem( (ih

fem,fc - 2^ dx ^ ^

rEBi k

Step 2. Define the interpolations /(wfem,fc) G V (Q) satisfying

I(Ufem,k)(x{i)) =v£m,k VxW Gfi.

Clearly, the vector t = i,^2)T, = I(ufem,k)i approximates the equilibrated

vector.

Step 3. For k = 1, 2 calculate

Qfem:k = :k) / &Xk: Sf = f ^fem,l ^fem^*

Step 4. For ai(x) + a2(x) = 1 and 0fc(x3_fc) = I(uiemik)(ak(x3-k), x3_fc) fc =

1, 2,

Step 5. Evaluate the bound ||Vufem — t||2.

Functions 0k(x3_k) specify boundary conditions for t. They can be defined differently from Step 4. and for their evaluation by means of the FE solution special more accurate procedures can be used. The simplest choice for a is a = 1/2, however a number of more sophisticated procedures for the evaluation of this function can be

a

ysis. Obviously, the procedures of Algorithm 2.2 are a part of Algorithm 2.3. The whole process of the a posteriori estimation can be arranged in the following way. One uses Algorithm 2.2. If the a posteriori estimate is unsatisfactory then additional calculations are performed according Algorithm 2.3 with some choice of ak. Further perfection of the a posteriori estimator is possible in a variety of ways. For instance, it is not necessary

t

t

which, e.g., is orthogonal inside the domain and provide some hp FE interpolation for fluxes, determined by the FE solution ufem.

Remark 2.2. Evaluation of a posteriori error bounds according Algorithm 2.1 Algorithm 2.3 involves only three operations

numerical differentiation with the use of finite differences, interpolation, and evaluating of 1-d integrals.

For this reason, these algorithms are obviously optimal in the arithmetic operations count, if the mesh is orthogonal. In the case of an arbitrary quasiuniform triangulation, the evaluation of integrals may be often arranged by layers of elements. From layer to layer the number of points, at which we need to evaluate an equilibrated flux may in general double. Therefore, the computational cost of the third among listed operations is estimated as O (nk n|_k), where nk is the maximal number of nodes in one layer and n3_k is the number of layers. In this paper, we concentrate on basic facts of a posteriori estimation, but several recipes can be immediately suggested for the reduction of the computational work, even in the case of nonuniform unstructured meshes. For instance, we can cover the computational domain by the nonuniform orthogonal mesh with the hanging nodes, matching in density the FE grid. Then we calculate one or both fluxes at the nodes of this mesh by means of averaging and interpolation. After that with the use of the introduced orthogonal mesh, we obtain corrections, which are necessary in order to make fluxes equilibrate.

tk (x) $k (x3_k) (qfem,k + ak ) (Ck, x3_k ) dCk

Remark 2.3. There are known a posteriori estimates

|| V'Ufem - V'u||njf2 < (1 + e)|| V'Ufem “ y||o,0 + ^ llV ' Y “ /11-1,0,

||Vufem - Vw||0 n < (1 + e)y Vufem - y||0,o+ (2-28)

+ cn(l + ^||V-y-/||gin ve>0,

where y is an arbitrary vector and cq is the constant from the Friedrich’s inequality, see Repin/Frolov [41] and Neittaanniakki/Repin [35]. One of our estimates can be written in the form

||Vufem - Vu||2 о < ||Vufem - t||2 о < (1 + e)||Vufem - У||0 П"

+ (1 + _) ak(f -V ■y)(i]k,x3-k)di]k

' e' k=1,2 ,■ \

Ofc(®3-fc)

(2.29)

о,о

In particular, in Algorithm 2.3 we used у = (/(wfem,i), Д«*т,2))Т • The right part / enters the last a posteriori estimate in a more adequate way, than in the second estimate (2.28). At the same time, it is easily computable, whereas the negative norm, entering the right part of the first bound (2.28), makes the bound difficult for the use in practice.

2.3. Heat conduction problem with discontinuous coefficient. The extension of the advocated approach and in particular of the algorithms of the previous section to the elliptic equations with the discontinuous coefficients is straightforward. Let us consider as an example the boundary value problem

-V • (p(x)Vu) = f (x)

du | 8v'Tn

e fi = (Q, 1) x (0,1),

ulrB =

7- L =0,

(2.30)

with rD, rN defined as in (2.2), p(x) > 0 and

pi = const for x e 01 := {x e 0 : 0 < x1 < 0.5},

P(x) =

P2 = const for x 6 O2 := \ ^l}-

For simplicity it is assumed that the boundary conditions are consistent and there exits such u0 e H2(Q) that = #. We define the approximate solution ufem of this

problem as the function belonging to the set L (Q) = V0(Q) + u0 and satisfying the identity

ap(ufem,«) = (/,v)n Vv e Vo(Q), (2.31)

where

ap(v,w) = y pVv -Vwdxidx2

and V0(Q) is the FE space defined in Subsection 2.1. If to introduce the norms

1/2

Mv,v))

1/2

p t • t dx

(2.32)

v

p

and imply by Qf the set of the equilibrated vectors satisfying (2.5), (2.6), then the a posteriori error estimate (2.8) takes the form

H V(u — Wfem)||p < UpVufem — t||p-i . (2.33)

Algorithms 2.1 2.3 are easily adapted to the problem under consideration. For instance, the first one is written as follows.

Algorithm 2.4.

Step 1. For each node g QQ calculate the value of the mesh function Vh = (v2,i(*))"1 j2=0> which is the finite-difference approximation of d(pdw/dxi)dxi(x(i)). For internal nodes of horizontal mesh lines x2 = hi2, i2 =0,1,..., n, use

vi,l(i) = p((il — 0.5)h, *2)[ufem(hi) — Ufem(h(il — 1), hi2)],

V2,l(i) = J^[vi,l(il + 1, *2) - *’l,l(*)] , il = 1, .. . ,11 - 1.

For the nodes (0, i2) on the axis x1 = 0 set

V2,i(i) = pid2Mo/dxi2(0, h*2),

where w0 is the 3-d order interpolation polynomial of xi over the values ufem(x(i)) for ii = 0,1,2, Mid dufem/dxi(0, hi2). For the nodes (n, i2) calculate

V2,i(i) = P2d2«i/dxi2(1, h*2),

where wi is the 3-d order Lagrange interpolation polynomial of xi over the values ufem(x(i)) for ii = n — 3, n — 2, n — 1,n.

Step 2. Define I2ji(x) c V (Q) as the piece wise bilinear interpolation of Step 3. Define t by evaluating the integrals

xi

ti(x) = J l2,i(Ci,X2) d£i,

0 x2 (2.34)

t2(x) = — J(f (xi, £2) + I2,i (xi ? ^2 )) d£2.

Step 4. Evaluate ||pVufem — t|p-1.

3. Linear elasticity problems

3.1. A posteriori estimation for plane problems. Let E be the Young’s modulus, v - the Poisson’s ratio, I - the unit tensor and tr( k ) = k : I - the trace of a tensor k . The linearly elastic plain strain problem in some domain Q is formulated in terms of the displacement vector u(x) = (ui(x), u2(x))T and symmetric strain and stress tensors

^ii <7i2\ _ / eii ei2

u

, CT2i ^22 j ' V£2i £22

related by the system of equations

div u + f = 0 , (3.1)

e(u) = (eki(u))kl=12 , £ Id = ^(duk/dxi + dui/dxk), (3.2)

t ' E <T(£

tr(e)I + є

1 - 2v

(3.3)

supplied by boundary conditions. The mixed homogeneous boundary conditions may have the form

u|rc =0, u n = t „|rN = 0, (3.4)

where un(x) and t n(x) are stresses normal and tangential to the boundary at a point x G Tn-

We turn to the case of the first boundary value problem when rN = 0 and for the boundary conditions we have

u|dn = 0. (3.5)

We assume that for positive constants cE,... ,cu

cE < E(x) < eg, cu < v(x) <cv< 0.5,

and introduce the space V = [i/i(Q)]2. The generalized solution u of the problem

(3.1) (3.3), (3.5) formulated in respect to displacements satisfies

I u (u) : e (v) dx = j f • v dx,, V v G V. (3-6)

n n

For simplicity it is assumed that the domain of the FE assemblage coincides with Q. For approximate solution of (3.6), the subspace V = [V (Q)]2 = [V (Q)]2 n V, where

V (Q) is the space of FE scalar functions. The FE solution is found from the identity

J u (ufem) : e (v) dx = y f • v dx, V v G V(Q). (3.7)

nn

The Hook’s law (3.3) can be written in the inversed form

e(cr) = +z/ti-(cr)I], (3.8)

ue

i/2 I \ i/2

Ieje = \ I u (e) : e dx I , | u |CT = I u : e (u) dx I . (3.9)

If у rnd є (3.8), then clearly і є|є = і у |ct , if additionally у = у (u) and

є = є (u) are related by equations (3.2), (3.3), then the energy norm for displacements

is defined according to the expression

|u|y = [у у (u) : є (u) dx| . (3.10)

The set Qf of equilibrated stress tensors is specified as Qf = т + Q0, where т is

any tensor satisfying the equilibrium equations

div т + f = О

and Qo is the linear space of self-equilibrated functions, i.e., satisfying equilibrium equations with f = 0. The latter space can be considered as a Hilbert space with the scalar product

[a, a ']CT = j a : e (a') dx n

and the norm | • ||CT. For the exact and the FE solutions u = (w1;w2)T, ufem = = (ufemj1, ufem,2)T, respectively, of the plain strain problem (3.6), we have the a

posteriori estimate

|u - Ufemlu <11 a (ufem) - T |CT V T G Qf . (3.12)

For obtaining a good tensor t , we use the same approach as before. The system of the equilibrium equations are sub-defined and specify the set Qf of equilibrated vectors Qo

f=0

Qf

Tkl by a sufficiently smooth arbitrary function and find other components from the equilibrium equations. Where no special assumptions on f G L2(Q) are made, we always assume that f G L2(Q). In algorithms of a posteriori estimators, it is sufficient to point out the way of definition of sufficiently smooth tensors t g Qf.

Algorithm A.

1. We specify an arbitrary function ^12 G L^(Q), arbitrary functions ^kk,rfc,_ (x3-k) G L^ [03-k, 63—k], and set Ti2 = ^12 •

2. Find q1 = d^12/dx2 and

Xi

T11 = ^11,ri,_ (x2) - J (/1 + ?1)(^1,x2) d^1. (3.13)

oi(x2)

3. Find q2 = d^12/dx1 and

X2

T22 = ^22,r2,_ (X1) - J (/2 + q2)(x1 ,£2) d^2. (3.14)

02(1x1)

As well one can also start from specifying

two functions ^fcfc,rfcj_ (X3_k) G LTO [a3_fc, 63—k], function ^12,rfcj_ (x3_k) G LTO [03—k, 63—k], and function q(x).

Then the tangential stress is defined by the integral

T12 = ^12,rfc,_ (x3—k) - J q(£k ,x3_k) d^k

°fc(xX3_k )

for one of k = 1,2, and other stresses are defined according to the above algorithm.

If steps 1-3 are used for obtaining an a posteriori estimate, functions ^12 and ^kk,rfc,_ are calculated by means of the obtained FE solution. This should be done in the most accurate way (e.g., with the use of superconvergence properties of the FE solution), since the closeness of these functions to the true stresses a12(u) on Q and to

akk (u) on the part of the boundary rk,_ is crucial for the accuracy of the a posteriori estimate. In particular, ^12 can be specified as an element of the FE space V (Q) with the nodal values obtained by the procedure of averaging similar to the used in Steps

2, 3 of Algorithm 2.2. In the case of a rectangular domain and orthogonal mesh, this procedure is especially simple

^12(x(l)) = 0"fem, 12 (i 1 h - 0,i2h - 0) + <7fem,12(*1h + 0, *2^ - 0)+ (3 15)

+ ^"fem, 12 (* 1 h - 0,*2h + 0)+ <7fem,12(*1h + 0,*2h + 0).

Here <Tfem,12 = <r12(ufem) is the stress defined by the FE solution and it is assumed that ^fem,12(*1h ± 0,i2h ± 0) = 0, if the corresponding element Tj, j = (j ± 1, j2 ± 1) does Q

In the case of a more general domain, rk,_ c rN and the homogeneous boundary condition (3.4) on rN it is adopted ^kk,rfc,_ (x3—k) = 0. If the second boundary condition is nonhomegeneous, i.e.,

a n|rN = t, (3-16)

then we set ^kk,rfcj_ (x3—k) = tkk or ^kk,rfcj_ (x3—k) = ?kk, where tkk is a more convenient for the use in (3.13), (3.14) but sufficiently accurate approximation of t. The estimate of the error tkk - tkk of approximation in some norm will enter the right part of the a posteriori estimate.

The approach under consideration may be realized in a number of ways. For the first step one can specify Tkk = ^kk and then find t12,t3—k,3—k from the equilibrium equa-

xk , k = 1 , 2

process of obtaining such equilibrated symmetric stress tensor, this function is differentiated twice in xk and integrated twice along x3—k. The latter requires more smoothness from ^kk at least in xk. In general an additional differentiation will certainly result in cruder a posteriori error estimates, if it is not compensated by the integration along the same direction, and this was observed in our numerical experiments. However, more complicated algorithms, but having the same asymptotical computational complexity,

xk , k = 1 , 2

time provide continuous equilibrated stress tensors. An example of such algorithms for the case of Q = (0,1) x (0,1) is Algorithm B.

Algorithm B.

1. Using the finite element solution, define a continuous piece wise bilinear function

I (k h2) G V, which “approximates” the second mixed derivative

7Y h ) ^ dqi2 _ E ( d3t<-i d3u2 \

12 dxidxo 2(1 + 7/) \dx\dx2 dx^dxo J

of the stress <r12 = <r12(u). Here, kh2 = (x12(i)j is the mesh function and, e.g., I(kh2)(x(i)) = K12(i). Inverted commas stand for the reason that indeed da12(ufem)/dx1dx2 may not be defined even on finite elements, and, therefore, some special ways of evaluation of k ^2 should be implemented. They should allow us to expect approximation in some sense of <r12(u) by T12(ufem), evaluated in a posteriori estimator by means of k h2. For instance, below t12 is defined by the backward double integration of I (k h2) in such a way that under some conditions one can expect the same accuracy from t12 as from <r12(ufem). If the square bilinear or higher order elements are used, one can evaluate da12(ufem)(x(i))/dx1dx2 for each element at its nodes, than for each node of FE assemblage calculate K12(i) as the average of da12(ufem)(x(i))/dx1dx2 for each element.

2. Evaluate

X2 Xi Xi X2

T12 =J J I (k ^2) dx1dx2 + Co + J C1(x1) dx1 +J C2(x2) dx2,

oo

Xi

T11 = J /1 - JI (k 5*2) dx1

0 L 0

X2 f X2

T22 = J /2 - JI (k 12) dx2

dx1 + C3 (x2) - C2 (x2) x1,

dx2 + C4 (x1) - C1 (x1) x2,

(3.17)

where

Co — ^12 C3(x2) — <711

Xi = 0,X2=0’ C1(x1) — d712/dx1 |X2=0, c2 (x2) — d712/dx2

xi = 0:

xi=0‘

C4(x1) — 722

(3.18)

X2 = 0"

and 7kl may be understood as exact or FE values of the stresses.

There are many other ways of finding the appropriate stresses rkl by means of the FE solution, which can be used as starting ones for evaluation of the equilibrated stress

T

In order to illustrate the essential difference from the approaches used by other authors, we formulate below basic a posteriori estimates in a form, in which the error in the smoothed FE stresses and the residual are separated.

Lemma 3.1. Let v be arbitrary vector in V = [V1(Q)]2, a(v) be the stress tensor satisfying (3.2), (3.3) for u = v and, y = {ykl}| l=1 be an arbitrary symmetric tensor with the components in H 1(Q). Then for u - v either of the estimates

|u - ufemlu < |a(ufem) - T||CT,

|u - ufemlu < |a(ufem) - y|<r + |St|CT,

|u - ufem|u < |a(ufem) - y|<r +

E

k=1,2

1 -V1

E

2 \ 1/2

ak(x3_fc)

, dyk,k dy1,2 w , ,

fk - —— 7^---------) (m,x3-k)dm

dx

k

dx

3—k

(3.19)

T

xk

Tl2=yi2, Tkk=ykk(ak(x3-fc))+ J fk ~ Q^’ ^j (Vk, %3-k) dl]k

ak (x3_k)

and

S t

ST11 0

0 ST22 1 '

St

kk

ak (x3_k)

r dyk,k dy1,2 w ,

fk-----—------~------ ) №,x’3-fc)c%-

dxk dx3—k

T

k = 1 k = 2

Xk

Tkk=Vkk, T12 = 2/12(a.*(x3_fc)) + J ^ fk - g^fc’fc ^ x’3-fc) di]k,

0k(X3_k)

Xi

0

T3-k,3-k = y-3-k;3-k(a-3-k(xk)) + J ^f3-k,3-k---------------------& ^ (xk, Vs-k) d??3-fc-

«3-fc(xfc)

Proof. Since all tensors t , appearing in Lemma, satisfy the equilibrium equations, the first estimate (3.19) clearly holds. Tensor t = y + St also satisfies the equilibrium equations. For this reason, the proof of the rest estimates requires only the inversion of the stress-strain relations and application of the Conchy and triangular inequalities. □

y

and made matching the equilibrium equations in a weak sense (remind that in general y

3.2. Linear elasticity and more general problems of solid mechanics in

3-d. In the case of 3-d elasticity problem, the stress tensor

/711 712 713

a = 1 721 722 723

\731 732 733

satisfying the equilibrium equations

d(Jk l do-fc 2 dfffc 3

dxi dx2

= /fc, k = 1, 2, 3. (3.20)

can be obtained in a similar to the described above ways. For, instance we specify the shear stresses 712, <713, 723 by some sufficiently smooth functions, approximating the stresses, specified by the FE solution. Then the rest stresses are found from the equilibrium equations and their boundary values, specified either by the boundary conditions in stresses or by the approximate values, found by means of the FE solution. We will not describe these obvious algorithms and restrict ourselves to the formulation of a statement similar to Lemma 3.1.

For definiteness of the norms | • | u, | • |CT, | • |e one can assume Hooke’s law for the homogeneous linearly elastic body

E E

7fcfc = ^ _ 9^ [(1 ~ v)£kk + ^(gfc+l.fc+l + £fc+2,fc+2)], Cfci = Y^\~^£kh

However, according to the above discussion, the first two estimates (3.21), given below, hold in a much more general situation under assumption of the proper definition of the norms | • |u, | • |<t.

Lemma 3.2. Let v be arbitrary vector in V = [ij1^)]3, a (v) be the stress tensor satisfying (3.20), (3.2) for u = v and, y = }|,l=1 be an arbitrary tensor with the

components in H 1(Q). Then for u — v either of the estimates

|u — Ufem|u < |a(Ufem) — T|ct,

|U — Ufem|u < |a(Ufem) — y|<7 + |St|ct ,

|U — Ufem|u < |a(Ufem) — y|<7 + ^ ^

3

1 f ft dyfc,fc+1 dyfc,fc+2A , S J

-7= [Ik - -7;------------^------------^------ {i]k,xk+1,xk+2)di]k

ak

+ E

k=1

holds, where ak = ak(xk+1,xk+2), t is stress tensor with the components

Tk1 = 7 k ^

Xk / \

Tfcfc ykk{,Qjk'> *^fe+2) H~~ I ( fk —o —o ) %k-\-2) djTjk

J \ dxk+i dxk+2 /

ak

and

/Jin 0 0 \

6t = I 0 Jt22 0 I ,

\ 0 0 6T33/

xk

x f ft dyfc,fc dyfc,fc+i dyfe,k+2A , s ,

v7~kk I I Jk r\ r\ o ) 7 7 *^fc+2 j dTjk •

7 V dxk dxk+i dxk+2 /

ak

Proof. The proof is similar to the proof of Lemma 3.1. □

There are several other sequences of constructing symmetric equilibrated stress tensors. We can start from setting Tk1 = yk1 with arbitrary functions yk1 G H 1(0) for three components Tk1. Besides (1,2), (1,3), (2,3), there are other admissible combinations of M: (11), (22), (12); (22), (33), (23); (11),(33),(31); (12), (13), (22); (31), (32), (22). The rest stresses Tmp are found from the equilibrium equations by 1-d integration. Boundary values, entering these integrals, are specified by arbitrary sufficiently smooth functions

Tmp I r Tmp (ap(xp+1 ,xp+2),xp+1 ,xp+2) G H ^ (rp,-)*

' 1 P,—

Let us underline that the equilibrium conditions do not depend on the type of the Hooke’s law, e.g., for orthotropic, transversally isotropic or other types of elastic bodies. As well they are not changed for a wide range of physically and geometrically nonlinear solid bodies. Therefore, the ways of obtaining of equilibrated and self-equilibrated stress tensors, introduced in this paper, are applicable to a wide range of problems in solid mechanics. All mentioned factors influence only techniques of evaluation of the (smoothed, if necessary) stress tensor, corresponding to the approximate solution, being subjected to a posteriori error estimation, and the specific energy norms, in which error estimation is produced.

Let us turn to a general case of nonlinear problems of solid mechanics, for which the approximate solutions obtained by means of the Lagrange principle of virtual work and Castigliano principle of complementary work provide an upper and a lower bounds for the true potential energy of the body. As it is well known, in this case the a posteriori estimate can be written in the form

L(v) — L(u) <L(v) — C (t), (3.22)

and under some conditions

^||u — v||v <L(v) — C(t), 0 = const, (3.23)

where L ( v)

v

uL

C ( t ) is the functional of the complementary work on the stress tensor t , satisfying the equilibrium conditions, and

|| • ||v is a norm satisfying the inequality

^||u — v||v <L(v) — L(u). (3.24)

The estimate (3.22) expresses basic properties of the Lagrange and Castigliano principles see. e.g., de Vebeke [15]. Arthurs [4] Mosolov/Myasnikov [34]. Washizu [54]. Berdichevskii [8]. where references on the earlier publications can be found. The estimate (3.23) follows from (3.22) under the condition that (3.24) is fulfilled. In a mathematical setting the basic facts for validity of the bounds (3.22). (3.23) may be found. e.g., in Ekeland/Temam [14] and Duvaut/Lions [13]. Glowinski [16]. They are found in the duality theory of the variational calculus and the theory of monotone/coercive operators. The latter allows to formulate conditions on the smoothness of data and the type of nonlinearity under which (3.24) holds.

The primal problem to be solved may be formulated in the following way: find u G U such that

L (u) = inf L (v), U G V, (3.25)

LV space with the norm || • ||v and U is a closed convex subset of V. The variable t and the functional C ( t ) of the complementary work are the dual variable and the dual functional with respect to the primal variational problem (3.25). At that t belongs to the set Q/ = Qo + t / of tensors satisfying equilibrium equations, e.g., (3.20), with Q0 being the space of the self-equilibrated tensors. The problem of finding the stresses bymeans of the Castigliano principle is the dual problem: find a G Q = Q0 + t„, such that

C (a )= sup C (t ), U G V. (3.26)

T £Q

Assume that ( —C (t0)) is also a proper convex, lower simicontinuous functional, which

Q0

C ( t ) < C (a ) = L (u) < L (v) V v G V, V t g Qo,

see.fi.g., Ekeland/Temam [14].

Many authors consider the use of the a posteriori estimates (3.22), (3.23) computationally very costly for the two reasons. One is that it is allegedly impossible to find the equilibrated tensor t close to the exact stress tensor a in a direct and sufficiently cheap way. Another reason is based on the conviction that numerical solution of the dual problem for finding t is much more difficult that the numerical solution of the primal problem. Such reasons are soundly pronounced in some contemporary publications. However, in the preceding sections we have demonstrated that the first reason in a delusion. In what follows we consider some additional ways (including practical algorithms optimal in the arithmetic operations count) for finding equilibrated stress tensors close to the solution. Apart from that we will show that it is with no doubt feasible to develop numerical techniques for solving dual problems, which are comparable with the most efficient numerical techniques for solving primal problems in respect of the computational cost.

3.3. Examples of algorithms for numerical testing. Below we illustrate the described approach by the two algorithms for obtaining a posteriori estimates in the case of the linear plain strain elasticity problem in the square Q = (0,1) x (0,1). We use notations e1 = (1,0), e2 = (0,1), whereas afem stands for the stress tensor, determined by the FE solution.

Algorithm 3.1.

Step 1. For each point = (h*i+h/2, hi2+h/2), ik =0,1,..., n—1, calculate 12(s<.>,= + using the finite element solution Ufem •

Step 2. Calculate approximate values of dafem,12/dx2 at the middle points = (hi1 + h/2, hi2), of the horizontal mesh intervals:

l(i) ^fem, 12 (?/^) (7fem,12(l/^ 62'>) • n 1 1 ' 1 o 1

' =---------------------------------, *1 = 0,1,..., n - 1, *9 = 1,2,..., n - 1.

h

Step 3. Evaluate approximate values of dafem,12/dx2 at the internal nodes x(i):

#° = ^(^ei)+^W), ifc = 1,2,1.

Step 4. Calculate values ^(n,i2) for i2 = 1, 2, ...,n — 1 by means of

linear extrapolation from the interior of ^ over the two nearest values on the mesh line x = i2h. In the same way calculate the values ^(i1,0) and for

ii = 1, 2,... ,n — 1.

Step 5. Determine values of for i1, i2 = 0, n corresponding the vertices, for instance, as the mean value of the two linear extrapolations along the two edges.

Step 6. Determine the piece wise bilinear continuous interpolation I(^h) e G V (0) of the mesh function ^h = (^(i) )"1,j2=0-

Step 7. Evaluate components of the stress tensor t satisfying the equilibrium

equation (3.11).

Step 7.1. For xk e (0,1), define as piece linear continuous functions c12(x1) —

— ^fem,12(x1, 0), Cn(x2) — CTfem,11(0, X2^, and C22(x1) — CTfem,22(x1, 0) . For c12

C12(*1h) = CTfem,12(*1h, 0), for i1 = 0, n,

C12(i1h) = <7fem,12(*1 — 0, 0) + CTfem,12(*1 + 0, 0), for i1 = 1, 2, . . . , n — 1,

and similar procedures are used for ckk.

Step 7.2. Determine components

X2

T12 = C12(x1)+y I(^h) dx2, T21 = T12, (3.27)

0

Xi

T11 = C11 (x2) — J (/1 + I(^h)) dx1, (3.28)

0

X2

t22 = c22(xi) - y ^/2 + -7^7^ dx2 (3.29)

0

of the equilibrated stress tensor t.

Step 8. Calculate |a (ufem) — t |ct for the a posteriori estimate.

Let us denote the stress tensor obtained by Algorithm 3.1 by t(1). If to change

x1, x2 x2, x1

equilibrated stress tensor t(2). Сіевхіу, the tensor t — «1 t(D + a2t(2) also belongs to Qf and

|u - Ufem|u <| CT (Ufem) - («it(1) + «2t(2))|CT, V «1 + «2 = 1. (3.30)

Another invariant to variables ж ^ ж2 procedure for finding a tensor t, satisfying the equilibrium equations (3.11). follows Algorithm В and is presented in the algorithm 3.2.

Algorithm 3.2.

Ufem —

— (ufem,b ufem,2) at the centers of the mesh cells y(i) — (h(i1 + 0.5), h(i2 + 0.5)), ifc — 0,1,..., n - 1,

(»i + 1/2,»2 + 1/2) _ ^'»-fem,l, s (ii + l/2,»2 + l/2) _ ^'»-fem,2 , s \

fem, 1,2 — qXo “fem,2,1 — Qx У'ьЧ,і2) {0.01 j

for i1, i2 — 0, 1, . . . , П — 1. Step 2. Define values

u(i1 + 1/2>i2 + 1/2) _ u(i1 + 1/2,i2 —1/2) -(ii+l/2,i2) _ “fem,1,2 “fem,1,2

*1,22 - i ■.

(3.32)

-.(ii + V2,i2 + 1/2) „((ii —1/2,i2 + 1/2)

~(»i ,i-, + l/2) _ “fem,2,1 “fem,2,1

«2,11 - h

for *i, i2 =0,1,..., n — 1.

Step 3. Define values

u(ii + 1/2,i2) (jl-1/2,j2) u(ii,i2 + 1/2) u(ii,i2-1/2)

~(*i,*2> U1,22 u1,22 ~(ii,i2) _ U2,11 u2,11

ul,221 — j? J w2,112 — j? (3.33)

for i1, i2 = 1, 2,..., n — 1.

Step 4. At the mesh nodes (hi1, 0) of the boundary for i1 = 1, 2,..., n — 1 . calculate "1*221 with the use of the linear extrapolation by the two nearest values u1i2211) and u1i12’221) • Calculate

u(ii,n) u(°,i2) u(n,i2) u(ii,°) u(ii,n) u(°,i2) U(n,i2)

“1,221, “1,221, “1,221, “2,112, “2,112, “2,112, “2,112

similarly. For the corner point (0,0), determine ■“ 1(°221 > e-9-> as the arithmetic mean of the two values obtained by linear extrapolations along two axes x1, x2 with the use of the two nearest values. Determine

~(°,n) ~(n,°) ~(n,n) ~(°,°) ~(°,n) ~(n,°) ~(n,n)

u1,221, u1,221, u1,221, u2,112, u2,112, "“2,112, "“2,112

similarly.

Step 5. Calculate

^_____/Vі) ,

12,21 — 9(1 -|- г/) v ’221 2>112J '

Step 6. Evaluate components of the stress tensor t satisfying the equilibrium equation (3.11).

Step 6.1. Determine

X2 Xi Xi X2

T12 = J 1I(cri2,2i) + Co + / Ci(xi) dxi + J 02(^2) dx2,

0 0 0 0

T2i = ti2, where I(<Ti2j2i) is the bilinear interpolation for <t(2 2i on the finite element mesh and

CO - 0fem,12(0,0), Cl(.X’l) ~ ^em’12(xi,0), C2(x2) ~ ^12(0,X3).

dxi dx2

Step 6.1?. For c3(x2) ~ fffem,ii(0, x2), determine

xi / xi \

Tii = I fi - I(<7i2,2i) dxi I dxi + C3(x2) - xiC2(x2). (3.34)

oo

Step 6.3. Define

X2 j X2 \

T22 = J (/2 " /1 (^2i) dx2) dx2 + C4(xi) " x2Ci<xi)’

0 \ 0 /

where C4(xi) ~ CTfem,22(xi, 0).

Step 7. Calculate the a posteriori estimator |a (ufem) — t | a.

We do not give formulas for evaluation of functions Ck, implying, however, averaging and interpolation procedures similar to those used in algorithms for Poisson equation. They provide the accuracy O (h2) for stresses, if understood as direct approximations of stresses corresponding to smooth displacements. Let us emphasize that less accurate approximation of the boundary stresses, than inside of the domain, can damage the accuracy of the a posteriori estimator. The choice of Ck can be optimized on purpose to minimize the posteriori estimator. The system of algebraic equations for finding such Ck has by the order of h smaller dimension. This allows to arrange computations in such a way that the optimization will not compromise the optimality of the a posteriori estimator in the computational cost.3

4. Equilibrated fluxes/stresses obtained by means of Castigliano principle

To some authors, dual formulations of the boundary value problems, expressing in mechanics of solid bodies the Castigliano principle, seem difficult for numerical solution. By this reason dual formulations are often discarded from consideration as a tool in the process of the a posteriori estimation. However, we will illustrate that numerical solution of dual problems may be as simple as of primal problems.

4.1. Poisson equation. Solution of the problem (2.15) minimizes the functional

J(v) = - [f,v)n Vv G V =

3See, V.S. Kostylev. A posteriori estimates optimal in the computational cost . Master thesis. Chair of Applied Mathematics. St. Petersburg State Polytechnical University. St. Petersburg, Russia, ‘2006 (in Russian).

If to define Qf as the set of functions satisfying the equilibrium equations in the generalized sense

The solution of the dual problem can be represented in the form z = z0 + tf, where tf any vector from Qf and the vector z0 G Q0 satisfies the equation

For deriving a discrete approximation of the dual problem, we use the same as before square mesh of size h = 1/n, n > 1, and define a subset Qh C Qf, which is represented

tf

For the purpose of discretization of the dual problem the simplest one can be used, e.g., as in (2.7) at q = 0. At the same time it is worth emphasizing that a good choice of tf

Remark 4.1. Instead of the vector t f, one can use to approximation f. Then the bound for the error of approximation ||tf — thy = O(hY)N(f) with some y and some norm N(f) of the function f will appear in the right part of the a posteriori estimate. The approximation f may be chosen on purpose, e.g., to simplify integration (since f can be obtained % approximation of f).

The self-equilibrated vectors 0ik = 1,2,..., n, should be chosen in a way which lead to the system (4.5) with good computational properties. This can be anticipated,

(4.1)

then the dual formulation of the problem is: find such z G Qf that

J*(z) = min J*(t), J*(t) = - [ t t clx.

y ’ tea, y h y ’ 2 J

(4.2)

(4.3)

n

Qh = tf + Qh,

where Qh is a ^^^^e dimensional sub space of Q0 and tf is any fixed vector from Q f. The approximate solution zh G Qh satisfies the equation

(4.4)

n

Qh = span [0 (i)(x) = (^^(x), ^2*^ (x))T]ie/h, where Ih is the appropriate set of indices i, then (4.4) is reduced to solving the system

Cw = f,

(4.5)

with the matrix C and the vector f defined as

C = {c*,j , w = , f = {f ,

(4.6)

n

n

if they have localized supports. Each 0(j) introduced below has for the support the set k* = w* n Q, where w* = {x : h(ik — 1) < xk < h(ik + 1), k = 1, 2}. Before doing this, we

remind the notation t* = {x : h(ik — 1) < xk < hik} for the square nests of the mesh.

First we define subsidiary local functions q(j) which at the definition of self-equilibrated

q

{1, x G t* U Ti1+i,i2+i.

— 1, x G t*i+ii*2 U t* 1 2+i, (4-7)

0, Q\Ui.

For the master linearly independent vectors, denoted as ^(j) = (yu^,^^)Tj and the

tf

Xi X2

M(ij)(x) = / q(j) (n,x2) dn, ^2j)(x) = — f q(i)(xi,n) dn,

0 X2 0 (4-8)

tf (x) = (0, t/,2)T, tf,2(x)^y f (xi,n) dn-

0

Clearly, the supports of these vectors are the sets w* , and instead of (4.8) one can use

Xi X2

Mij)(x) = J q(j)(n,x2) dn, ^2j)(x) = — J q(i)(xi,n) dn. (4.9)

(j i —i)h (i 2 i) h

The vectors 0(j ) are defined as restrict ions to k(* ) = w(j) n Q of the vect ors ^(j ), determined in (4.9), whereas Ih = {i : 0 < ik < n}.

Now we compare the system (4.5), generated with the use of the coordinate vectors 0(i)

KDuD,fem = fD!fem, KNuN,fem = fN!fem, (^-10)

generated by the spaces V (Q^d V (Q) for the Poisson equation in the unite square with the homogeneous Dirichlet and Neumann boundary conditions on its boundary, respectively. Remind that V (Q) is the space of the continuous piece wise bilinear functions and V (Q) is its subspace of functions van ishing on dQ .Let p(j) (x) be the standard piece wise bilinear continuous function, satisfying conditions p(j) (x (j)) = £*,3 with dij being the Kronecker’s delta, and p(j) (x) be its restriction to Q. It is easy to conclude that

= *(.) = fMl (4 H,

1 dx,dx,' V \dx2 dx, j ' ' '

whence it immediately follows that C = KN.

The system (4.10) with the matrix KD defines the FE solution of the primal problem. In order to obtain the equilibrated vector valid for the a posteriori estimation of

C = KN

Clearly, the both can be solved very efficiently by many fast solvers, developed for FE methods.

Remark 4.2. It is easy to note that the introduced localized fluxes satisfy the equality

^0(j) =0 Vx G Q, (4.12)

ie/h

C = KN

y0 = 1 = {y0j) = 1}*£/h with the unity for all entries. This eigenvector corresponds to zero eigenvalue. According to (4.12). the right part in (4.5) satisfies the solvability condition, and if w is the solution, then w + cl is also solution, where c is an arbitrary number.

K C

coincides with the FE matrix for the problem

—Am =f (x), x G Q, dQ = Td U Tn,

I 0 d I 0 (413)

m|tn = 0, dM/dn|rD = 0,

with rD, rN defined as in (2.2).

C K

grids by higher order rectangular finite elements. The situation is more complicated for triangulations by triangular elements compatible in C, because the second mixed derivative do not exist. However, if any elements of the class Ci are used, then the analogy is

2

retained. If the differential operator has the form Lm = ^ d(dakl/dxk)/dx;. Coeffi-

fc,i=i

cients in the dual formulation will be entries of the matrix A-1, where A = {akl }| l=i

C

induced by the same coordinate functions as in the FE method for the primal formula-C

the matrix A-1.

Below we obtain another characterization of the of the space Qh and the set of vectors {0(j) }ie/h, used for generating system (4.5) with the matrix C = KN. This characterization will illuminate what kind of approximation is used for the dual problem. Yet another characterization by means of the integral equation with respect to the q

We return again to the Dirichlet boundary value problem (2.26) in the unite square.

Let

g(j) = ^ f0r x G Tj n Q, (4.14)

0 for x G Q \ t* ,

and

—*3 — fc / \ j 1 for x3 — k G h(i3—k ^-, i3 — k^ i3 — k 1, 2,...,n, / a -t

Mfc0 (x3—k) = <n . .. ! ■ \ (4.15)

^0 I°r x3 —k G h(i3—k 1,i3 —k).

The set {g(j)}ni i2 = i is the basis in the space Lh(Q), which is a discrete approximation

for L2(Q) and contains piece wise constant functions on the uniform square mesh of

h

g = £ a(j)g(j)( *1,*2=i

x)

play the same role as q in (2.7). In ton, fanctions {/^0 fc (x3—k)}n3—k=i serve as the basis for the approximation of the value of the boundary flux at xk = 0. In

accordance with such understanding, we define the finite dimensional space of the selfequilibrated fluxes as = sPan [t(i)]ie/h where = {i = (ik, i3—k), ik =0,1,..., n, i3—k = 1,2,..., n} and

Xi X2

t(1i)(x)^y g(i)(Ci,X2) d£i, 4*}(x) = g(i)(xi,6)) d^2 (4.16)

0 0

40’i3-fc)(x) = ^ (X3_k), 40—Tfc )(x) = 0. (4.17)

It is easy to see that functions q(i), 1 < ik < n, are linear combinations of g(i) for the same i. Therefore, for these i fluxes 0(i) are linear combinations of fluxes t(i). Fluxes 0(0>*3-fc) are obtained by means of q(i), ik =0, 1 < i3—k < n and boundary fluxes k (xk )• Therefore, we have proved the following Lemma.

Lemma 4.1. Tfte space Qq(Q) = span [t(i)]ie/h, spanned over fluxes (4.16), (4.17), and the space Qq(Q) = span [{ 0(i) (x)]ie/h, spanned, over localized fluxes (4.9), coincide.

Remark 4.4. Solution of the system generated with the use of the coordinate functions {t(i)}ie/h may be unstable. The reason is not that this system has a bad condition number, which indeed is O (h_2). The system is equivalent to the discretized integral equation of the first kind - see next subsection - and at h ^ 0 it smallest eigenvalue tends to zero.

Remark 4.5. Suppose that we have to solve the Poisson equation with the Neumann boundary condition

z„|dfi = t„, (4.18)

which for the dual formulation is an essential one. In this case, the approximate solution of the dual problem is represented as

n— 1

Z° = Zq + tf,N = Z° + tN + tf, Z0 = £ w(i) 0i,

il,i2 = 1

where the vector tN found from the Neumann boundary condition and the vector tD,f is defined as in the preceding case of the Dirichlet boundary condition, see (4.8). The vector tn may be specified by the vector, defined on Q and, therefore, one simply can set tN = tn. Coefficients, w(i) for i G i0 := {i : 1 < ik < n — 1, k = 1, 2,} are found from the system

Cw = f, (4.19)

with the matrix C and the vector f defined as

C = K,}n—i=0, f = {/(->}n—1=0, Cj = / 0«'0«dx,

Q X2 (4.20)

/(i) = J 0(i) ■ tf dx, tf (x) = (0,t/i2)T, t/,2(x)^y /(xi,n) dn.

Q 0

In the particular case of the homogeneous boundary condition du/dv|Q = 0, one has tN = 0, and (4.20) uniquely define (4.19), which in turn has the unique solution. There is no difficulties in defining appropriate vector tN satisfying the nonhomoge-neous boundary condition (4.18). Particularly, tN can be defined in such a way that

h

approximation for tN by some vector tN- We can set

tN = £ c* 0, i G idQ = {i : x(i) G dQ},

with the coefficients c« chosen in such a way that the trace tN,dQ = tN1N is the

tn

Zo

||tn — tN dQ||2 , will appear in its right part.

Clearly, the case of an arbitrary sufficiently smooth domain does can be treated similarly

4.2. Some remarks and generalizations. Several obvious, but having important consequences, remarks can be made.

a) Since the Dirichlet boundary condition is natural for the dual formulation, we can use orthogonal grid for obtaining equilibrated fluxes in the case of arbitrary sufficiently smooth domain Q.

3) For solving the dual problem it is not necessary to use the same mesh, which was used for FE discretization. More over, since equilibrated fluxes from Q0, see (4.1), are not supposed to satisfy compatibility conditions, we can easily add any equilibrated coordinate vectors to the basis {0(i) }n i2 =0 and enrich the space Q°. For instance, additional coordinate vectors can be defined with the use of local finer mesh, arbitrarily oriented with respect to the mesh used for the definition of the basis { 0(i) }n i2=0- We can add also coordinate vectors with specific properties admitting a better approximation of concentration of fluxes or their singularities.

Arbitrary domain. Since the Dirichlet boundary condition is natural for the dual formulation, we can use orthogonal grid for obtaining equilibrated fluxes in the case of

Q

Suppose, we would like to obtain the equilibrated fluxes for the Dirichlet boundary value problem (2.26). Formally, the formulation of the dual problem is not changed, and again we have to solve integral identity (4.3) with the use of Q0 defined for a Q

h

= {i : mes [ w(i) n Q = 0]} and Q0 = spanieJ^ [0(i)], where each 0(i) is defined as the restriction of ^(i) to Q. It is necessary to underline that since mes x(i) can be small for some i G the condition of matrix C of the system (4.5) for the problem under consideration can be bad. However, due to the discussed at the end of the preceding subsection analogy with the FE systems, several simple remedies for improving the condition can be used. We refer in this relation to Korneev [26] and Oganesian/Ruhovets [38].

Densening of the mesh. For solving the dual problem it is convenient to use the same mesh, which was used for FE discretization (e.g., for evaluation of the norms entering a posteriori error estimate), but not necessary. For instance, since equilibrated fluxes from Q0, see (4.1), are not supposed to satisfy compatibility conditions, we can easily add any equilibrated coordinate vectors to the set of such vectors, spanning the space Q0 > see, e.g., (4.16), (4.17), and enrich this space upto some space Q0, * ^ Qo- For instance, additional coordinate vectors can be defined with the use of local finer mesh, generating Qq . We can add also coordinate vectors with specific properties admitting a better approximation of singularities in fluxes. A good source of functions, which can serve for generating the localized equilibrated functions are coordinate functions used

in meshlcss methods, see Oden/Duarte/Zienkiewcz [40] and Strobolis/Babuska/Copps [47].

We consider only a simple example of densening the mesh. Suppose that in a strictly internal subdomain D c Q it is necessary to use more accurate approximation. Since the case of an arbitrary sufficiently smooth domain was discussed above, we restrict

Q

way. Let Dh is the least “mesh domain” covering D , i.e.,

Vh = UieJ,7t IhD = {i: T? nv^0},

tq is the nest of the mesh of size h. Each square nest tq , i G ID, is subdivided in four squares, defining the mesh of size ft = h/2 on Dh. Retaining old in dices i = (ii, i2) for old nodes, we add indices i = (i1 ± 1/2, i2 ± 1/2) for new nodes. We use the notation td

i

of td , and introduce sets

ID = {i : tr nDo = 0}, IQ = {i : tQ n Q = 0},

^ = {i : t/0 n Q = 0}, /M = IQ\Dh U ^-

It is convenient to use the common notation t« for t0, i G d h anfl f°r td ,

i G /d. Now we can directly use (4.14), (4.15) and (4.16), (4.17) for defining the coordinate vectors of equilibrated fluxes, spanning the space of equilibrated fluxes, which we denote Qh’D. The dimension of Q°’d is card[/h,D] + 2n.

The above consideration shows that in essence the densening is quite simple. However, the basis vectors (4.16), (4.17) are not localized and, therefore, the matrix of the corresponding system will have considerable fill in. Apart from that, solving this system may be unstable for the reason pointed out in Remark 4.4. The instability can be removed by the transformation of the introduced coordinate self-equilibrated fluxes to the localized ones. Instead of one type, two types of the localized self-equilibrated fluxes are used. If to use the notation 0 = 0(i) for the fluxes, introduced in (4.7)-(4.9) on

the mesh of size h, the second type fluxes 0 Q/ are similarly defined on the mesh of h/2

Interpretation as solution of integral equation. We turn to the Dirichlet problem (2.26) in an arbitrary sufficiently smooth domain and consider the equivalent system of two equations

d2u(1) (i)i

—j- = axf-q, U^|an = 0,

d2u(2) , , (2), „ (4.21)

-^- = «2 f + q, «W|an = 0,

u(1) = u(2) Vx £ Q.

In the mechanical sense, this system describes the two systems of strings stretched along axes xk with each string of one direction fastened to the strings of other direction at the cross points. Function q is the internal force, acting between the two systems of strings. For simplicity we assume again that each line xk = const crosses Q not more than in two points, rk,_ and rk,+ are the parts of the boundary containing the points of such pairs, having lesser and larger coordinates xk, respectively. We write the equations defining the curves rk,_ and rk,+ as xk = ak(x3_k), and xk = (x3-k)

for a3-k < x3-k < b3-k .Let G k (xk , x3-k , yk) be the Grin’s functions for the ordinary differential operators in (4.21). so that

bfc(x3-fc)

u(k)(x) = J Gk(xk,x3-k,yk)(«k/ - q)(yk,x3-k) dyk-

ak (x3-k)

Satisfying the equality u(1) = u(2), one comes to the integral equation

bk(x3-k)

^ ^ I Gk (xk, x3-k, yk)q(y^ x3-^ dyk

k=1,2 , v

ak(x3-k)

bk(x3-k)

= 53 / Gk(xk,x3-k,yk)ak/(yk,x3-k) dyk- (4.22)

k=1,2 , s

ak(x3-k)

In Rozin [45]. Korneev/vRozin [29. 30]. the class of integral equation of a more general but similar to (4.22) type was termed integral equation of the method of splitting.

For discretization of the integral equation (4.22), one can use the space Gh(Q) of the piece wise constant functions with the basis {g(i)}ie/!h, /q = {i : t0 n Q = 0}. As we will show below, in this way we come to the system equivalent to (4.5) up to the choice of the basis functions and ak. However, from (4.22) it becomes clear that this system is not good for the numerical solution. Since (4.22) is an integral equation of the 1-st kind, this basis will lead to the unstable system of algebraic equations, the matrix of which has considerable fill in. The use of another set of coordinate fluxes {q(i)}, which produce local self-equilibrated fluxes, results (as in Subsection 4.1) in the system, which computational properties are the same as of the FE system for the Poisson equation with the Neumann boundary condition. In order to make more clear the interrelation between solving procedures of the integral equation and the dual problem (4.2), we note that the generalized formulation of (4.22) is: find q G L2(Q) such that for any q G L2(Q) we have

bk(x3-k)

/ q(x){ ^2 / Gk (xk ,x3-k,yk)q(yk ,x3-k) dyk} dx =

bk(x3-k)

q(x){ ^ / Gk(xk,x3-k,yk)ak/(yk,x3-k) dyk} dx. (4.23)

k=1,2 {° , ak (x3-k)

Taking into account the equalities

bk(x3-k)

/ G.(x.,x3-t,K)q(yk,x3-k) dyk =0 for = ak(x3-k>, (x3-t).

ak(x3-k)

and integrating by parts, we obtain

bk(x3-k) [ bk(x3-k)

J q(x) [ J Gk(xk,x3-k,yk)q(yk,x3-k) dyk [ dxk

a k ( x 3 -k) ^a k ( x 3 - k )

bk(x3-k) [ Xk

J < ^(x3-k)+ J q(nk,x3-k) dnkX

ak(x3-k) ^ ak(x3-k) J

(bk(x3-k)

J Gk{xk,x3-k, yk)q(yk,x3-k) d,yk > cixk,

ak(x3-k)

where ^(x3_k) is an arbitrary sufficiently smooth function. By the definition of the Gk

bk(x3-k) Xk

J G k{xk,x3-k, yk)p(yk,x3-k) dyk = -t0,k{x3-k) - J p{m,x3-k) diik,

ak(x3-k) ak(x3-k)

in which —q0,k(x3-k) = q0,k(p(x), ak(x3-k), x3-k) is the boundary value for the derivative in the left part and, therefore, is uniquely defined by the function p. Suppose that for discretizing the problem we use the basis {g(i)}ie/fc and

bk(x3-k)

+(*)

toi}(x) = (4*i >4*2) , 4*k = ~^Tk J Gk(xk,x3-k,yk)q(i:)(yk,x3-k)dyk,

ak(x3-k)

bk(x3-k)

^/(*^) J ~3x~ j" 3'3—ki Vk^&kf (jjki 3'3—k') d/ljki

ak(x3-k)

Qq = span [t^i)]ie/^ . Then (4.23) can be reformulated: find z = zQ + tf, where tf was defined above and the vector zQ G Q0 satisfies the equation

J(z„ — tf )t0j) dx = 0 V i G lO - (4.25)

Let us underline that the coordinate functions g of the internal forces are not self-equilibrati functions.

equilibrated, but the fluxes t0J) are due to their definition by means of the green’s

Gk

o,i 0,i / \ f f,,

k = xk (x3-k), xk = xk

ating vectors t0l),tf. If the points x.’* = xk’*(x3-k), x. = xk’k(x3-k) are such that

bk(x3-k) bk(x3-k)

J (xk — xk’*)q (i) (x) dxk = 0, J (xk — xk)a/(x) dxk = 0,

ak(x3-k) ak(x3-k)

then

t0i,)k(x) = t0i,)k’-(x3-k) + (—1)(3 k) J q(i)(yk,x3-k) dУk,

ak(x3-k)

(xk

tf,k(x)=tf,k’-(x3-k)+ / ak/(yk,x3-k)dyk,

(4.26)

ak(x3-k)

where

t0I-(x3-k ) =

bk-x°k bk ak

^f’k,- (x3— k)

7 0,i bk — xk’

bk(x3-k)

/

ak(x3-k)

bk(x3-k)

q(i) (x) dxk,

bk ak

ak/(xk) dxk-

ak(x3-k)

The formulation (4.25) differs from (4.3) not only by the choice of the basis, but also in the spaces of coordinate functions. Indeed, we have Q0 c Q0 and instead of (4.8) the relationships (4.26) are used. That means that vectors t0*kj f,k satisfy additional conditions in comparison with t0« k ? f,k entering (4.3). If one defines

,,(<)

(x)

t0«k(yk,x3-k) dyk, Wk(x)

f,k (yk, x3— k ) dyk,

ak(x3-k)

ak(x3-k)

then

(i) I

vk ,wk|ao

0-

according to (4.24), (4.26) and definitions of Green’s functions.

Remark 4.7. In this paper, we do not discuss convergence of the described Riesz Galerkin methods for the minimization of the functionals of the complementary work or corresponding integral equations of the splitting method. The analysis of the convergence does not meet difficulties. For the techniques which can be applied and results we refer to Korneev/Rozin [29, 30] and Korneev [19 21, 25, 27]. In these works, simpler for the realization, but more complicated for the analysis discretizations were studied. For instance, in Korneev [20, 21] 1-d integrals in (4.22) were approximated by the trapezium quadratures on a rectangular grid and then the collocation method was applied for obtaining the system of algebraic equations. Namely, it was required that (4.22) holds for the quadrature nodes. Let us also note, that integral equations obtained by the method of splitting for the bending problem of thin plates and cylindrical shells were studied by Korneev/Rozin [29, 30] and Korneev [22] [25]. Construction of localized

equilibrated functions of internal forces for thin shells and shells of moderate thickness was completed in Korneev [27] under rather general assumptions on the configuration of the middle surface. Also in [27] the analysis of the convergence may be found for numerical methods, based on the use of these equilibrated functions for the minimization of the functional of the complementary work for shells with arbitrary sufficiently smooth middle surfaces.

Remark 4.8. For the boundary value problems of mechanics of solid bodies two level of splitting are distinguished: i) splitting of the equilibrium equations and ii) splitting of partial differential equations of a boundary value problem in respect to displacements. The latter means that we are able to split the equilibrium, stress-strain. strain-displacements relations and boundary conditions. When this is possible, one can obtain integral equations of the method of splitting, which were introduced in mechanics of solid bodies by Rozin [45. 46]. From his and other works, mentioned in Remark 4.7. it follows that in rather general case ii) can be accomplished, if the Poisson ratio v is zero. However, for an efficient use of the a posteriori error estimation algorithms of the specific class under consideration one needs only to split equilibrium equations, which is always possible. The latter is true for the possibility of obtaining discrete dual formulations based on the Castigliano principle, which are comparable in the computational cost of their solution with FE methods for primal formulations in respect to displacements.

Remark 4.9. Let ufem be the function from the FE space interpolating the exact solution u. We rewrite (2.8) in the form

||Vu — Vufemll < ||V(«fem — Ufem)|| + ||VUfem — z|| V Z G Qf, (4.27)

and assume that the FE space is the space of the continuous piece wise bilinear functions. z

Note that all functions in the right part are from the finite dimensional spaces. For convenience, let us call by the gauge order the order of convergence of the norm in the left, given by the a priori estimate. Due to the superconvergence property, the first term in the right can be estimated with an additional with respect to the gauge order, see, e.g., Oganesian/Ruhovets [38], Korneev [28] and Whalbin [53], Babuska/'Strobolis [7].

q

estimate of the second norm with the gauge order. Therefore, the accuracy of the a posteriori error estimate is at list the same in the order as of the a priori estimate.

Comparison of discrete primal and dual formulations. The equality C/ = = K takes place in a much more general case. To illustrate this we turn to the Dirichlet problem (2.26) in an arbitrary sufficiently smooth domain and its FE discretization (2.27). We can assume that the finite elements of the FE assemblage are arbitrary which can provide that V (Q) G C(Q) n H1(Q). In other words, the finite elements are allowed to be be curvilinear and associated with the triangular or rectangular reference element with any compatible in C(Q) n H1(Q) shape functions. The Hermite finite elements are not excluded, but we number the FE Galerkin basis functions of the space V (Q) consecutively with the use of the number l = 1, 2,..., L without making difference between basis functions, corresponding to the values of FE functions or their derivatives at the nodes. Therefore, L is the total number of the FE Galerkin basis functions, for which we use now the notation p[1] (x). The number of the internal basis functions is denoted by L/ so that and V (Q) = span|p[i1]i=1 . The finite element solution satisfies the identity (2.27). The basis self-equilibrated vectors 001 = (^00 ]i, ^0 ’2)T in the space Q0 = span[00 1];-Li can kg defined by means of the FE basis functions according to

'0 Ji=1

i-fc

x3-k

^]fc = (-l)1-fc7p—, k = 1,2, / = 1, 2,..., £/. (4.28)

0x3-k

[l

It is clear that vectors ^0 1 k satisfy the equilibrium equation

df}_ df}_= = 9xi <9x2 dxidxn dxodxi

in classical sense only when the finite elements are compatible in C1. In general, for elements compatible in C, second derivatives in (4.29) are Dirack’s deltas on the borders of finite elements, and therefore (4.29) involves equalities for the Dirack’s deltas corresponding to /dxkdx3-k for k = 1, 2. However, in the weak sense, see (4.1),

the equilibrium equations are satisfied. Now we see that K/ = C/, where

K = {kl,m};,me£J , C = ,

k;,m = J VpM • Vpm dx, c;,m = J 0(1) • 0(m) dx. (4-30)

n n

Clearly, in the case of more general equations of second order, e.g., when the coefficients

K/ C/

turn to the equation V- p Vu = f where p is a diagonal 2 x 2 matrix p = diag[p1, p2], then K/ is the same as C/ for the similar equation with p = diagjp-1, p-1]. Indeed, the basis in Qf is the same, and the coefficients of the stiffness and deflection matrices are the integrals

k;,m = I VpM • p Vpm dx, c;,m = I 0(1) • p-10(m) dx.

n n

From the above considerations, one can conclude that at least for regular elliptic problems computational properties of the discretizations of primal and dual problems are in essential the same and solutions of these diecretizations can be obtained by fast solvers of the same types.

Remark 4.10. For 3-d Dirichlet problem (2.15) we have for the flux t = (t1; t2, t3)T the balance equation

!^+!^+!^+/=o. (4.3D

dx1 dx2 dx3

Therefore, in order to obtain an equilibrated flux, two of the components can be specified by arbitrary functions and only third found from (4.31). Again, sufficiently smooth local functions can be used for generating self-balanced fluxes. Suppose that ^>(x) has local support S and d3^/dx1dx2dx3 is bounded in the vicinity of S. Then for sufficiently smooth functions k =1, 2,3, components of a self-balanced flux are defined by

d2t

tk = o.k-------£----, At = 1,2, 3. (4.32)

dxfc+1dxfc+2

Suppose, V (Q), Q = (0,1)3, is the space of the continuous piece wise linear function on the FE cubic mesh of size ^d ^(i) (x) satisfies ^(i) (x(j)) = Sjj, where i = (i1, i2, i3), j = j j2, j3), x(i) = hi, x(j) = fe^d h = 1/n. Then substituting ^(x) = ^(i)(x) in

(4.32), one obtains local self-balanced fluxes t(i) = (t^, 4*), 4J))T •

4.3. Linear elasticity problems. State of plain stress. In this section we will use the matrix-vector form of the stress strain relations for the state of plain stress

E h -v v 0 De, D = ----------—--------- v l—i/ 0

(1+ v)(1 - 2v) 1 0 0 (1 - 2v)/2y

where

= (a11, a22, ^12)T, £ = (£11, £22, Y12)T, 712 =2e12,

and wo turn to the problem (3.6). The solution of the dual problem is the stress vector z = zo + tf, where t f any vector from Qf and the vector zo G Qo satisfies the equation

J(zo - tf )D-1to dx = 0 Vto G Qo. (4.33)

n

In order to discretize this integral identity, one can proceed along the lines of Algorithms A or B. Suppose for simplicity that the domain is covered by the square mesh of size h Qh is some mesh domain containing Q and defined below, and

V(nh) = {v : G C(jf),t>l G Qp, p > 1},

where Q p is the space of polynomials of the the order p in each variable. For defining *h(

a discrete subspace Qh(Q) C Qo(Q), at first we define the space of stresses <r12 as the

restriction V (Q) of the space V (Qh) to Q. The stresses are evaluated according

(3.13), (3.14) with = 0 and functions ^kkirkj_ from appropriate finite dimensional spaces of traces. The coefficients before the basis functions of these spaces are clearly additional unknowns in the discrete dual formulation. It is possible to avoid special description of functions ^kfc,rfcj_ % choosing a proper basis in V(Qh). At the same time, it is possible to define the basis functions in such a way that their supports will be localized on the squares, containing 9 mesh sells. Since the generalization to any p > 1 is obvious, we will describe one of the bases for the case p = 1.

Let p(i)(x), p(i)( x(j)) = be the usual nodal coordinate function in the space of

continuous piece wise bilinear functions, and

^(i)(x) = p(i)(x) — p(i1-1,i2)(x) — p(i1,i2-1)(x) + p(i1-1,i2-1)(x), (4.34)

so that supp (x)] = w(l)^d = {x : (ik — 2)h < xk < (ik + 1)h}. First, we

define the master basis vector ^by the equalities

Xk

A‘12 = 4>{t)(x), Mfcfc = - J ^(m,x3-k)di]k.

(ik-2)h

Elements of the basis {0o’i}ih, Ih = {* : K(i) := ^(i) n Q = 0} in Qh(Q) are the

restrictions of ^ (i^o Q. Note that for Qh one can take the domain with the closure

Cl = Uj-hZjt1'1. Discrete formulation of (4.33) is the following one: find such vector z = zo,h + tf with z° G Qh(Q) that

/(zh — tf )D-10h,i dx = 0, V 0h,i G Qh. (4.35)

3-d elasticity. As it is seen from Section 3, the construction of the space of the self-equilibrated stresses for 3-d is similar to 2-d case. The same is true about the construction of the basis functions in the space Qh, which have local supports. In the matrix-voctor form the stress strain relations are

a = De.

(1 + v)(1 — 2v)

/1 — v v v 0 0 0

1 — v v 0 0 0

1 — v 0 0 0

(1 — 2v)/2 0 0

SYM (1 — 2v)/2 0

(1 — 2v)/2

(4.36)

where

a = (^11, ^22, ^33, ^12, ^23, ^13)T, e = (£11, £22, £33, Y12, Y23, Y12)T, Ykl = 2£fc;, k = l.

As for the 2-d problem, we cover the domain by the square mesh of size h, consider some domain Clh D Cl containing Cl, and the space V(Clh) = {u : v G C(Cl ), u|r G G Qp, p > 1}, where Qp is the space of polynomials of the the order p in each variable. For Qh we take the mesh domain containing all nests ri, i = (i1, i2, i3) involved in the

Q oh

Let p(i)(x), p(i)( x(j)) = Ji,j be the usual nodal coordinate function in the space of continuous piece wise bilinear functions, and

#(x) = p(i) (x) — p(il-1’i2’i3)(x) — p(il’i2-1’i3)(x) + p(il-1’i2-1’i3)(x) —

— p(ii,i2,i3-1)(x) + p(ii-1,i2,i3-1)(x) + p(ii,i2-M3-1)(x) — p(il-1’i2-1’i3-1)(x), (4.37)

so that supp [^(i) (x)] = w(i)^d w(i) = {x : (ik — 2)h < xk < (ik + 1)h, k = 1, 2, 3}. We define the master basis vector ^(i) by the equalities

Xk

Llki = 4>{t)(x), A4fc = “ J (ydxk+1 + dxk+2 ) ^k,xk+i,xk+2)di]k, k ± I,

(ik-2)h

and the elements of the basis { 0 h,i}j h, Ih = {i : K(i) := w(i) n Q = 0^^n Qh(Q) as the restrictions of ^(i) to Q.

The discrete formulation of the dual problem has the same form (4.35) and requires solution of the system of linear algebraic equations with the banded matrix.

5. Numerical results

In this section, we discuss the results of numerical experiments with the equilibrium based a posteriori error estimates described in previous sections. The purpose of our experiments is to demonstrate that our algorithms are able to produce such estimates with the very good effectiveness index and for the optimal in the order number of arithmetic operations. For model problems, linear and nonlinear second order elliptic equations in the unite square were used, including the equation with jumping coefficients and the plain strain linear elasticity problem. Main conclusions made from numerical results are that our a posteriori estimates

are asymptotically exact, and, more over, in many cases convergence of the effectiveness index to the unity was observed at h ^ 0,

are asymptotically optimal in the computational cost, the number of arithmetic operations was always proportional to the number of unknowns,

can be easily made robust in respect to coefficients jumps after necessary modifications of the algorithms.

We tested also the a posteriori estimators in which the equilibrated fields were obtained by solving the dual problem, expressing the Castigliano principle.

Fig. 1. Energy norm of FE error and a posteriori estimator against number of unknowns for the problem (5.1)

5.1. The Poisson equation.

5.1.1. Direct evaluation algorithms. Consider the model problem

. 5n2 3n n

— Au = —— cos ~^xi cos ~^x2 , \Xi,X2) G £2 = (0,1) x (0,1),

u\rD = 0, ulr n = 0, (5.1)

rD = {(xi,X2)\ xi G [0,1] ,X2 = 1}U{(xi,X2)\ Xi = 1 , X2 G [0,1]} ,

rN = {(xi,X2)\ Xi G [0, 1] ,X2 =0}U{(xi,X2)\ Xi =0 ,X2 G [0, 1]} .

with the exact solution

/ \ 3n n . .

U{Xi, X2) = cos — X\ COS — X2 . (5-2)

The FE space of the piece wise bilinear functions for this problem V 0 (^) was defined in Section 2. For the FE solution ufe^ ^esh of size h, we used Algorithm 2.1 in

order to calculate the vector-valued function t(x) = (ti(x),t2(x))T , which satisfies the balance equation (2.5) and the boundary conditions (2.6). Then we calculated the energy norm of the error e = ||V(u — ufem)^^ ^ ^^^riori estimator n = ||Vufem — t||,

i.e., the left and right sides correspondingly in (2.8).

Fig. 1 shows the dependence of the energy norm of the FE error and of the a posteriori estimator on the number of unknowns N. It demonstrates the same asymptotic behavior of the both values. The number of unknowns in this experiment exceeded

4 • 10^ ^ ^mctically coincide for N greater than 104. The a posteriori es-

timator n stays greater than the energy norm of the error e (see Table 1). This validates that the equilibrium based a posteriori estimate guarantees the upper asymptotically exact estimate.

Fig. 2 shows the dependence of (/eff — 1 on N), where

Jeff = - (5.3)

e

is the effectiveness index of the a posteriori estimate. We see that Ieff converges to 1 rather fast and always stays grater then 1, see also Table 1.

Table 1

N e r] Jeff

16 8.40422 10-1 2.87729 10-1 3.42362

64 4.08785 10_1 8.38575 10_1 2.05138

256 2.02318 10_1 2.67228 10_1 1.32083

1024 1.00877 10_1 1.09368 10_1 1.08417

4096 5.04023 10“2 5.14654 10-2 1.02109

16384 2.51966 10-2 2.53287 10-2 1.00525

65536 1.25977 10“2 1.26141 10“2 1.00131

262144 6.29880 10“3 6.30085 10“3 1.00033

1048576 3.14939 10“3 3.14965 10“3 1.00008

4194304 1.57469 10“3 1.57473 10“3 1.00002

left - 1

101 102 103 104 105 106 107 Ar

Fig. 2. Dependence Ieff — 1 on the number of unknowns N for the problem (5.1)

We also compared the computational costs of the a posteriori estimator and optimal multigrid solver for the problem (5.1). These results are shown on the Fig. 3 and demonstrate that the equilibrium based a posteriori estimator is optimal with respect to the number of arithmetic operations. Moreover, computation of the a posteriori estimator is about twice cheaper than solving the finite element system (this results were obtained on AMD Athlon 64 3200+ 2.01 GHz with 2 Gb of RAM).

5.1.2. Algorithms based on the Castigliano principle. In the unite square, we considered the Dirichlet problem having for the exact solution

u = sin(2nx1) sin(nx2) + (x1 + 1)(2x2 + 1)

The corresponding right part is f = 5n2 sin(2nx1) sin(nx2) and nonhomogeneous Dirichlet boundary condition is u\dQ = (x1 + 1)(2x2 + 1).

The FE solution ufem was obtained by means of the space of the continuous piece wise bilinear functions. For obtaining the approximate solution zh of the discretized dual problem (4.4), we used the subspace Q with the local flux basis vectors 0, described in Subsection 4.1. Two choices for the vector tf and respectively for the set Qh were implemented: one according to second line of (4.8) and another according to

Fig. 3. Computational costs (in ms) of multigrid solver and n against the number of unknowns for the problem (5.1)

Table 2

N e Jeff Jeff T2 1eS

16 1.59912 2.93605 2.05747 1.20610

64 7.60498 • 10_1 1.95927 2.16500 1.23228

256 3.70781 • 10_1 1.26742 2.21586 1.24493

1024 1.84023 • 10_1 1.07113 2.23084 1.24869

4096 9.18344 • 10“2 1.01826 2.23475 1.24967

16384 4.58949 • 10“2 1.00460 2.23574 1.24992

65536 2.29446 • 10-2 1.00115 2.23599 1.24998

262144 1.14720 • 10“2 1.00029 2.23605 1.24999

1048576 5.73594 • 10“3 1.00007 2.23606 1.25000

the formulas

Xk

tf(x) = tfAx) = \ J fiVk, xs-k) drjk-

0

For these two choices of fluxes, satisfying the balance equation, we introduce notations tf = tf1^, tf2). These fluxes result in the two a posteriori estimators

nc,j = ||VUfem - ||, Zh,j = ZQ,j + tfj}, j = 1, 2,

with the effectiveness indices denoted Iff. The direct evaluation algorithm was also applied to the problem. In essence it is Algorithm 2.1, but adapted to the case, when one has the Dirichlet boundary condition at rN as it is defined in (2.2). The adaptation was performed with the use of special algorithm, which optimizes the a posteriori estimator among different boundary fluxes in the appropriate finite dimensional space. The effectiveness index for the the estimator produced by the adapted Algorithm 2.1 is denoted Ieff.

The dependence of the energy norm for the error and of the described effectiveness indices on N is presented in Table 2.

Other problems were also solved by FE method on purpose of testing the Casigliano principle based a posteriori estimators. The numerical results showed similar behavior of the estimators of this class and allow us to come to the following conclusions.

a) The equilibrated fluxes obtained by means of the Castigliano principle provide good error estimators with good effectiveness indices, one of which stays below 1.3. ft) However, both in dices Iff do not converge at h ^ 0, whereas Ieff does.

Y) The computational time for solving the system of algebraic equations (4.5), resulting from the Castigliano principle, and for evaluating the a posteriori estimators is proportional to the number of unknowns N. Therefore, the algorithms are asymptotically optimal in the computational cost. However, the computer time is greater approximately in 1.5 times, than the computational time for the direct evaluation algorithm with the effectiveness index Ieff.

Item ft) can be explained by the fact that the FE method for solving the primal and dual problems have the same rate of convergence in the energy norms. Namely,

Apparently, the convergence of the effectiveness index to the unity for our algorithms of

of the FE solution ufem to the continuous piece wise bilinear interpolation of the exact solution u.

Item y ) completely approves conclusions made in Subsection 4.2 in the part titled Comparison of discrete primal and dual formulations.

The set Q f does not depend от the choice of the vector t f, which enter the definition of this set. However, the sets Q0’j := Qo + tfj) depend, and, according to the numerical results, some allow to approximate true fluxes better. Besides, the error estimator nc,2 is more symmetric with respect to the axes xk, Лап ПсД • Probably, these factors caused the difference of the effectiveness indexes, reflected in a).

5.2. Second order elliptic equation with discontinuous coefficient. We

tested also our a posteriori estimator as applied to the problem:

with the same boundary conditions as in (5.1) and the piece-wise constant coefficient, which has a jump across the common boundary y for the two parts of Q:

The right-hand side f as well as the mixed boundary conditions corresponded to the exact solution

||V(u - ufem)|| x || Vu - z|| = O (h).

the direct evaluation of the balanced fluxes is related to the fact of the superconvergence

V • (p(x)Vu) = f (x), x € П = (0,1) x (0,1),

(5.4)

= {x| xi € (0, 0.5), x2 € (0, 1)}, ^2 = {x| xi € (0.5, 1), x2 € (0, 1)} .

(5.5)

(5.6)

For the function p, we used

x € ^i, x € ^2-

(5.7)

(x - l)2£l + 0.25— + 1.25, ж > 0.5.

V 'P2 P2 ’

x < 0.5,

For obtaining FE solution, we used the space V0(Q) of the continuous piece-wise bilinear functions, satisfying the Dirichlet boundary condition on rD.

Fig. 4. Energy norm of the finite element error and the a posteriori estimator against the number of unknowns for the problem (5.4)

Fig. 5. Dependence of Ieff — 1 on the number of unknowns N for the problem (5.4)

The vector t(x) = (t(x)1, t(x)2)T, satisfying the balance equation (2.5) and the Neumann boundary conditions (5.1) on rN, was calculated according to Algorithm 2.4 in Subsection 2.3. In turn, this vector allowed to evaluate the a posteriori estimator

Vp = UpVufem — t||p-i, which enters the estimate (2.33). Fig. 4 shows the dependence of the energy norm, of the error e = ||V(u — ufem)||P and of nP on the number of unknowns N ( in this numerical experiment N also exceeded 1 • 106, but for N > 104 the lines on the graph coincide). The effectiveness index tends to 1 rather fast, as it is illustrated by Fig. 5, in which the value of Ieff — 1 is plotted against N. At the same time the effectiveness index is always grater then 1 (see also Table 2) that validates that the a posteriori estimate (2.33) is a guaranteed upper asymptotically exact bound.

The comparison of the computational costs of the a posteriori estimator and the optimal multigrid solver for the FE system of linear algebraic equations is presented in Fig. 6. These results demonstrate the optimality of the a posteriori estimator. Note, that the a posteriori estimator is more than twice cheaper, than solving the FE system by the multigrid method.

Fig. 6. Computational costs of the a posteriori estimator and the multigrid solver against the number of unknowns for the problem (5.4)

Fig. 7. Energy norm of FE error and a posteriori estimator n(1) against the number of unknowns for the problem (5.8)

5.3. Linear elasticity problem. Algorithms 3.1 and 3.2 were applied to several linear elasticity problems with nonhomogeneous Dirichlet boundary conditions. Since the results reflect similarly the difference of these algorithms and the level of their efficiency, we present them only for one problem (3.1)-(3.3) in the unite square Q = [0,1] x [0,1]. The vector f and the Dirichlet boundary conditions correspond to the exact solution u = (u1, u2)T ,

ui (x) = sin(nx1) sin(2nx2) + xi + x2,

U2{x) = sin(27TXi) sin(7nr2) + j(xi + l)(x2 + 1).

(5.8)

Algorithms 3.1 and 3.2 produce the a posteriori estimators n(1) and n(2), respectively, according to their description in Subsection 3.3.

Fig. 7 and 8 demonstrate the behavior of n(1) • The obtained for n(2) are

presented in Fig. 9-11, see also Table 3 for numbers. The numerical results show, that

Fig. 8. Computational cost of multigrid solver and of the a posteriori estimator n(1) against the number of unknowns for the problem (5.8)

Fig. 9. Energy norm of FE error and the a posteriori estimator nagainst the number of unknowns for the problem (5.8)

Table 3

N e r] Jeff

64 256 1024 4096 16384 65536 262144 1048576 9.14308 4.44483 2.20395 1.09958 5.49486 • 10_1 2.74705 • 10_1 1.37348 • 10_1 6.86734 • 10-2 1.02954 • 101 4.77467 2.25318 1.10628 5.50364 • 10_1 2.74818 • 10_1 1.37362 • 10_1 6.86751 • 10-2 1.12604 1.07421 1.02234 1.00609 1.00160 1.00041 1.00010 1.00003

the a posteriori estimator n(2) outperforms n(1) • The effectiveness index /ffj of n(2) tends to 1 staying grater then 1, whereas the effectiveness index of n(1) does not converge and stays slightly greater than 2. Both a posteriori estimators are optimal in the computational cost. Fig. 8 and 11 present the comparison of their computational cost

10'4-------------------...........................’-'-l..............1.......------------------------------------------------------------------------------------‘...............................................................-

10' 102 103 10* 1(f 106

Fig. 10. Dependence of Ieff — 1 for the aposteriori estimator n(2) on the number of unknowns N for the problem (5.8)

Fig. 11. Computational costs of the multigrid solver and the a posteriori estimator n(2) against the number of unknowns for the problem (5.8)

Table 4

N e I(1) 7eff v(2) r(2) 7eff

16 3.29126 • 10_1 2.11430 6.42399 2.45022 7.44462

64 1.63481 • 10_1 5.54014 • 10_1 3.38885 3.54162 • 10_1 2.16638

256 8.17444 • 10“2 1.94481 • 10_1 2.37914 1.46349 • 10_1 1.79033

1024 4.08766 • 10“2 8.62833 • 10“2 2.11083 5.45946 • 10“2 1.33560

4096 2.04389 • 10-2 4.18226 • 10-2 2.04622 2.56000 • 10“2 1.10387

16384 1.02196 • 10“2 2.07433 • 10“2 2.02976 1.05054 • 10-2 1.02797

65536 5.10979 • 10“3 1.03501 • 10“2 2.02554 5.14644 • 10“3 1.00717

262144 2.55490 • 10-3 5.17229 • 10-3 2.02446 2.55953 • 10“3 1,00181

1048576 1.27745 • 10-3 2.58580 • 10-3 2.02419 1.27804 • 10“3 1.00046

and the computational cost of the multigrid solver for solving the problem (3.1)-(3.3) with the Dirichlet boundary conditions.

The difference in the behavior of the a posteriori estimators n(1) and (see Table 4) shows that the quality of the error estimators strongly depend on how accurately the properties of the superconvergence of FE solutions are taken into account at their evaluation. Clearly, the convergence of the effectiveness index can take place only if an a posteriori error estimator is superconvergent. Algorithm 3.1 is simpler, but it is not invariant in respect to axes xk and the accuracy of approximation of the boundary values for stresses provides only an approximation of the order of h.

Research is supported by the grant from the Russian Fund of Basic Research No 05-01-00779-a. The second author has partially been supported by the Austrian Science Fund (FWF) through the special research programme (SFB) F013. project 16.

Резюме

И.Е. Ануфриев, В.Г. Корнеев, B.C. Косты,лев. Точно уравновешенные поля: могут ли они быть эффективно использованы для получения апостериорных оценок погрешности.

Показывается, что конструктивное определение линеалов тензоров напряжений, удовлетворяющих уравнениям равновесия, является простой для численной реализации операцией. Это позволяет эффективно применять классические апостериорные оценки погрешности приближенных решений краевых задач, проистекающие из двух взаимно дополнительных принципов принципа Лагранжа минимума энергии деформации и принципа Кастильяпо минимума дополнительной работы, являющегося двойственным по отношению к первому. Применительно к задачам линейной теории упругости в таких оценках энергия погрешности приближенного решения, удовлетворяющего всем геометрическим условиям, оценивается энергией, отвечающей разности тензора напряжений приближенного решения и любого тензора напряжений, удовлетворяющего уравнениям равновесия. Вопреки распространенному мнению о большой вычислительной трудоемкости построения уравновешенных тензоров, близких получаемым посредством МКЭ (метода конечных элементов), мы показываем, что во многих случаях это может быть сделано за оптимальное число арифметических действий. Доказываются также новые апостериорные оценки посредством неуравновешенных тензоров напряжений. По сравнению с известными оценками, содержащими, например, в случае уравнения Пуассона норму невязки (в уравнении баланса для используемого вектора потока) в пространстве H-1, они вычисляемы и более точны. Приводятся ряд алгоритмов вычисления апостериорных оценок для уравнения Пуассона и системы уравнений теории упругости и результаты численных экспериментов, подтверждающих весьма высокую эффективность алгоритмов и их робастность.

Literature

1. Ainsworth M., Demkowicz L., Kim C.-W. Analysis of the equilibrated residual method for a posteriori estimation on meshes with hanging nodes. - Submitted for publication.

2. Ainsworth M., Oden J.T. A posteriori estimation in finite element analysis. - N. Y.: John Wiley Sons, Inc., 2000. - xx+240 p.

3. Aubin J.-P. Approximation of elliptic boundary-value problems. - N. Y.-London-Sydney-Toronto: Wiley Interscience, a division of John Wiley & Sons, Inc., 1972.

4. Arthurs A.M. Complementary variational principles. - Oxford: Calderon Press, 1980.

5. Babus ka I., Reinbolt W.C. Error estimates for adaptive finite element computations //

SIAM J. Num. Anal. - 1978. - V. 15, No 4. - P. 736-754.

6. Babus ka I., Reinbolt W.C. A posteriori error analysis for finite element methods for one

dimensional problems // SIAM J. Num. Anal. - 1981. - V. 13, No 3. - P. 565-589.

7. Babus ka I., Strobolis T. Finite element method and its reliability. - N. Y.: Oxford University Press, 2001. - xi+802 p.

8. Berdichevskii V.L. Variatsionnye principy mahaniki sploshnoi sredy // Variational principles of mechanics of conytinuous media. - M.: Nauka, 1983 (in Russian).

9. Bernardi Ch., Girault V. A local regularization operator for triangular and quadrilateral finite elements // SIAM J. Num. Anal. - 1998. - V. 35, No 5. - P. 1893-1916.

10. Braess D., Scho berl J. Equilibrated residual error estimator for Maxswell’s equations. -Submitted for publication.

11. Carstensen C. Residual-based a posteriori error estimate for a nonconforming Reissner-Mindlin plate finite element // SIAM J. Num. Anal. - 2002. - V. 39, No 6. - P. 2034-2044.

12. Cle ment Ph. Approximation by finite element functions using local regularization // Rev. Francaise Automat. Informat. Recherche Operationelle Ser. Rouge Anal. Numer. -1925. - T. 9, F. 2. - P. 72-84.

13. Duvaut G., Lions I.-I. Inequalities in mechanics and physics. - Berlin; Hedelberg; N. Y.: Springer-Verlag, 1976.

14. Ekeland I., Temam R. Convex analysis of variational problems. - North-Holland; Amsterdam; N. Y.; Oxford, 1976.

15. Fraeijs de Vebeke. Displacement and equilibrium models in the finite element method // Stress Analysis / Eds. O.C. Zienkiewicz, G.S. Holister - London-N. Y.: Wiley, 1965. -P. 145-197.

16. Glowinski R. Numerical methods for nonlinear variational problems. - N. Y.-Berlin-Heidelberg-Tokio: Springer-Verlag, 1984.

17. Gol’denweizer A.L. Theory of elastic thin shells (Translation of original Russian edition of 1953). - N. Y.: Pergamon Press, 1961.

18. Kelly D. The self-equilibration of residuals and complementary a posteriori error estimates in the finite element method // Internat. J. Num. Methods Engrg. - 1984. - V. 20. -P. 1491-1506.

19. Korneev V. G. Primenenie metoda mehanicheskih kvadratur dlia postroenia raznostnoi shemy dlia ellipticheskogo uravnenia chetvertogo poriadka (Implementation of quadrature method for costruction of difference method for elliptic equation of 4-th order) // Metody vychislenii. - Leningrad: Isdatel’stvo Leningradskogo universiteta, 1971.- No 7. - P. 46-55 (in Russian).

20. Korneev V. G. O sviazi raznostnyh shem dlia ellipticheskih uravnenii vtorogo poriadka so shemami kvadratur dlia integral’nyh uravnenii metoda raschlenenia (On interrelation between difference schemes for elliptic equations of second order and quadrature method for integral equations of splitting method) // Vestnik Leningradskogo universiteta. -1973. - No 1. - P. 18-27 (in Russian).

21. Korneev V.G. Raznostnaya shema poriadka tochnosti h3/2 v norme || • ||h Wi dlia ellipticheskogo uravnenia vtorogo poriadka v proizvol’noi oblasti (Difference scheme of accuracy h3/2 in the norm || • || h Wj for elliptic equation of second order in arbitrary domain) // Vestnik Leningradskogo universiteta. - 1973. - No 7. - P. 25-33 (in Russian).

22. Korneev V. G. Svedenie zadachi rascheta arochnoi plotiny k zadache na minimum raboty uprugih sil nepreryvnoi sterzhnevoi sistemy (Reduction of the problem to minimization of work of elastic forces work for continuous rod system ) // Trudy Gidroproekta / Pod red. D.M. Yurinova. - M., 1973. - No 22. - P. 4-12 (in Russian).

23. Korneev V. G. Diskretizatsia raboty uprugih sil nepreryvnoi sterzhnevoi sistemy (Discretization of functional of elastic forces work for continuous rod system) // Trudy Gidroproekta / Pod red. D.M. Yurinova. - M., 1973. - No 22. - P. 12-32 (in Russian).

24. Korneev V. G. O ratsional’nom vybore koordinatnyv funktsii pri minimizatsii diskretnogo funktsionala raboty uprugih sil (On rational choice of the coordinate functions at minimization of the discrete functional of the work of elastic forces) // Trudy Gidroproekta / Pod red. D.M. Yurinova. - M., 1973. - No 22. - P. 86-102 (in Russian).

25. Korneev V. G. Issledovanie shodimosti priblizhennyh reshenii (Analysis of convergence of approximate solutions) // Trudy Gidroproekta / Pod red. D.M. Yurinova. - M., 1973. -No 22. - P. 141-182 (in Russian).

26. Korneev V. G. Shemy metoda konechnyh elementov poriadkov tochnosti (The finite elel-ment methods of high oreders of accuracy). - Leningrad: Izdatel’stvo Leningradskogo universiteta, 1977 (in Russian).

27. Korneev V. G. On numerical solution of the shell theory problems in stresses on skew-angular meshes // Comput. Maths and Math. Physics. - 1981. - V. 21, No 2. - P. 441-450.

28. Korneev V. G. Superconvergence of the finite element solutions in the mesh norms // Comput. Maths and Math. Physics. - 1982. - V. 22, No 5. - P. 1133-1148.

29. Korneev V.G., Rozin L.A. Algorithmy rascheta obolochek po sterzhnevoi sheme i ih matematicheskoe issledovanie (Algorithms for numerical solution of shell problems and their mathematical analysis) // Materialy 6-i vsesoyuznoi konferetsii po teorii obolochek i plastin, Baku, 1966, - M.: Nauka, 1966. - P. 555-560 (in Russian).

30. Korneev V.G., Rozin L.A. Chislennaia relizatsia metoda raschlenenia na primere zadachi izgiba plastinki (Numerical realization of the splitting method for the plate bending problem) // Izvestia Akademii Nauk SSSR, Mehanika tverdogo tela. - 1967. - No 2 (in Russian).

31. Ladevese P., Leguillon D. Error estimate procedure in the finite element method and applications // SIAM J. Num. Anal. - 1983. - V. 20. - P. 485-509.

32. Luce R., Wohlmuth B. A local a posteriori error estimator based on equilibrated fluxes // SIAM J. Num. Anal. - 2004. - V. 42. - P. 1394-1414.

33. MikhMn S. G. Variational methods in mathematical physics. - Oxford: Pergamon, 1964.

34. Mosolov P.P., Myasnikov V.P. Mechanics of rigid plastic bodies. - M.: Nauka, 1981 (in Russian).

35. Neittaanmaki P., Repin S.I. Reliable methods for computer simulation Error control and a posteriori estimates. - N. Y.: Elsevier, 2004. - 305 p.

36. Novozhilov V. V. Teoria tonkih obolochek (Theory of thin shells). Gosudarstvennoe soyuz-noe izdatel’stvo sudostroitel’noi promyshlennosti. - Leningrad, 1962 (in Russian).

37. Novozhilov V. V. Theory of thin shells. - Dordrecht: Noordoff, 1959.

38. Oganesian L.A., Ruhovets L.A. Variatsionno-raznostnye metody reshenia ellipticheskih uravnenii. - Erevan: Izdatel’stvo Armianskoi akademii nauk, 1979.

39. Oden J.T., Demkowicz L., Rachowicz W., Westermann T.A. Towards a universal h-p adaptive finite element strategy. Part 2, A posteriori error estimation // Comput. Methods Appl. Mech. Engrg. - 1989. - V. 77. - P. 113-180.

40. Oden T., Duarte C., Zienkiewcz O. A new cloud based finite element method // Comput. Methods Appl. Mech. Engrg. - 1998. - V. 153. - P. 117-126.

41. Repin S., Frolov M. Ob a posteriornyh otsenkah tochnosti priblizhennyh reshenii kraevyx zadach dlia ellipticheskih uravnenii. (On a posteriori estimates of approximate solutions of boundary value problems for elliptic equations) // Zhurnal vychislit. matem. i matem. phiziki. - 2002. - V. 42, No 12. - P. 1774-1787.

42. Repin S., Neittaanmakki P., Frolov M. On computational properties of a posteriori estimates based upom the method of duality error majorants // Numerical Mathematics and Advanced Applications (ENUMATH-2003). - Berlin-Heidelberg: Springer-Verlag, 2004. - P. 346-357.

43. Repin S.I., Xanthis L.S. A posteriori error estimation for elasto-plastic problems based on duality theory // Comput. Methods Appl. Mech. Engrg. - 1996.- V. 138, No 1. - P. 317339.

44. Reissner E. On the derivation of the theory of thin elastic shells // J. Math. Phys. -1963. - V. 42, No 4. - P. 263-277.

45. Rozin L.A. Metod raschlenenia v teorii obolochek (The splitting method in the shell theory) // Prikladnaia matematika i mehanika. - 1961. - V. 25, No 5 (in Russian).

46. Rozin L.A. Hekotorye voprosy raschlenenia i diskretizatsii uravnenii teorii obolochek (Some questions of splitting and discretization of equations of the shell theory) // Issle-dovania po teorii uprugosti i plastichnosti. - Leningrad: Izdatel’stvo Leningradskogo universiteta, 1965. - No 4.

47. Strobolis T., Babuska I., Copps K. The design and analysis of the generalized finite

element method // Comput. Methods Appl. Mech. Engrg. - 2000. - V. 181. - P. 43-69.

48. Stewart J.R., Hughes T.J.R. A tutorial in elementary finite element error analysis. A sys-

tematic presentation of a priori and a posteriori error estimates // Comput. Methods Appl. Mech. Engrg. - 1998. - V. 158, No 1-2. - P. 1-22.

49. Temam R. Problemes Mathematique en Plasticite. - Paris: Bordas, 1983.

50. Vejchodsky T. Local a posteriori error estimator based on the hypercircle method // Proc. of the European Congress on Computational Method in Applied Mechanics and Engrg. (ECCOMUAS 2004). Yavaskyla, Finland, 2004.

51. Verfu rth R.A. A posteriori error estimates for nonlinear problems. Finite element discretizations of elliptic equations // Math. of Comput. - 1994. - V. 62, No 206. - P. 445-465.

52. Verfu rth R.A. A review of a posteriori error estimation and adaptive mesh refinement

techniques. - BG Teubner: John Wiley & Sons, 1996. - vi+127 p.

53. Walbin L.B. Superconvergence in Galerkin finite element methods. Lecture notes in Mathematics. - Berlin: Springer, 1995.

54. Washizu K. Variational principles in elasticity and plasticity. - N. Y.: Pergamon Press, 1982.

55. Zienkiewicz O.C., Zhu J.Z. The superconvergence path recovery (SPR) and adaptive finite element refinement // Comput. Methods Appl. Mech. Engrg. - 1992. - V. 29,

No 1. - P. 78-88.

Поступила в редакцию 13.12.06

Ануфриев Игорь Евгеньевич кандидат физико-математических паук, доцепт кафедры «Прикладная математика» Санкт-Петербургского государственного политехнического университета.

E-mail: AnufrievQamd.stи.neva.ru

Корнеев Вадим Глебович доктор физико-математических паук, профессор кафедры «Прикладная математика» Санкт-Петербургского государственного политехнического университета.

E-mail: Korneev Qamd. stu. пелт. ru

Костылев Владимир Сергеевич инженер Всероссийского паучпо-исследователь-ского института гидротехники, г. Санкт-Петербург.

E-mail: Vul KosQmail.ru