Ausgangsfragestellung

The conventional strategies for mesh refinement in finite element methods are mostly based on a posteriori error estimates in the global energy norm. Such estimates reflect the approximation properties of the trial functions by local interpolation constants while the stability property of the continuous model enters through a global coercivity constant. However, meshes generated on the basis of such global error estimates are not appropriate in computing local quantities as, e.g., point values or contour integrals, and in the case of nonlinear material behavior. Recently a new technique for a posteriori error control and adaptive mesh design for finite element methods has been proposed in [25], [26], and [31]. Following the general approach of C. Johnson and his co-workers, [11] and [4], residual-based a posteriori error estimates are derived via duality arguments. Here, we carry this idea further into a feed-back method in which the dual solutions are approximated on the current meshes and used as local weights in the a posteriori error estimates multiplying the residuals of the computed solution. In this way local information about the mechanism of error propagation is captured which results in most economical meshes as well as accurate error bounds in the course of a dynamic mesh adaptation process. The additional work required by the computation of the weights is acceptable as it usually amounts to less than 30% of the total cost. This approach has been presented in [31] for primal as well as dual-mixed finite element methods in linear elasticity. Here, it is developed further for solving finite element models of Hencky- and Prandtl-Reuss-type in linear-elastic perfect plasticity. We concentrate on the primal-mixed formulation as in this case the yield-constraint can be incorporated into a nonlinear material law which reduces the linear variational inequality to a nonlinear variational equation and allows for the use of efficient solution methods. In the dual-mixed formulation the pointwise constraints have to be explicitly satisfied by the stresses which makes the solution of the resulting algebraic problems much harder.

We recall the basic concept of our method for a posteriori error estimation at a simple model problem due to Strang [18]. The physical problem is that of an infinitely long straight pipe, with quadratic cross-section  $\Omega$ $\subset$ $\mathbb {R}$² , filled with plastic material adherent to the walls and subjected to a volume force f acting in vertical direction (see Figure 1). The mathematical model seeks a scalar displacement  u  in the vertical direction and a stress vector $\sigma$ = ($\sigma_{1}^{}$,$\sigma_{2}^{}$) as functions on  $\Omega$ . The plastic behaviour of the material is taken into account by the nonlinear restriction  |$\sigma$| $\leq$ 1 . This results in the system

$$\sigma \sigma \Pi \nabla \Omega \partial \Omega$$ (1)

where $\Pi$ denotes the pointwise projection onto the unit circle. For the moment, let us assume that the stresses are small, so that the side condition can be neglected leading to essentially elastic behaviour. Then, (1) reduces to the linear boundary value problem

$$\Delta \Omega \partial \Omega$$ (2)

Abbildung 1: Geometry sketch of the Strang example

We discretise this problem by a conforming finite element method using piecewise linear or (isoparametricly) bilinear shape functions on triangular or quadrilateral meshes $\mathbb {T}$_h = {K}, respectively, satisfying the usual condition of shape regularity (cf. [2]). For ease of mesh refinement and coarsening hanging nodes are allowed in our implementation. The width of the mesh $\mathbb {T}$_h is characterised in terms of the piecewise constant mesh size function h = h(x),  0 < h$\le$1, where h_K : = h_{| K} = diam(K) and h_max = $\max_{K\in\mathbb T_h}^{}$ h_K. Using the notation V_h $\subset$ V = H¹₀($\Omega$) for the corresponding finite element subspaces, the approximate solution u_h $\in$ V_h is determined by the discrete equation

$$\nabla \nabla \varphi \varphi \forall \varphi \in$$ (3)

Here and below, ( ^. , ^. ) denotes the inner product of L²($\Omega$) and | ^. | the corresponding norm, while ($\nabla$ ^. ,$\nabla$ ^. ) is the natural energy inner product of problem (2). Further, H^m($\Omega$), for m $\in$ $\bf N$, denotes the usual m-th order Sobolev space with norm | ^. |_m, and H¹₀($\Omega$)  the subset of  H¹($\Omega$) of functions vanishing on $\partial$$\Omega$. These are all spaces of either scalar or vector-valued functions and no distinction will be made in the notation of the corresponding inner products and norms.

In [26] a new concept has been proposed for estimating the error in the scheme (3) for general error measures given in terms of linear functionals J( ^. ) defined on the space V, or on a suitable subspace containing the finite element space V_h and the exact solution u. Relevant examples are torsion moments, stress values, or the mean surface tension,

$$\psi \int_{\Omega}^{} \psi \sigma_{ij}^{} \Gamma \int_{\Gamma}^{} \sigma_{nn}^{}$$

Then, following a common strategy for residual-based a posteriori error estimation (see [11], [21], [4], and [26]), we utilize the solution  z $\in$ V  of the corresponding dual problem

$$\nabla \varphi \nabla \varphi \forall \varphi \in$$ (4)

to derive a representation of the error  e = u - u_h,

$$\nabla \nabla$$ (5)

Further, using the Galerkin orthogonality relation  ($\nabla$e,$\nabla$$\varphi$_h) = 0 , $\varphi$ $\in$ V_h, and element-wise integration by parts, it follows that

$$\sum_{K\in\bf {T}_h}^{} \left\{\vphantom{(f+\Delta u_h,z-z_h)_K - \tfrac{1}{2} ([\partial_nu_h],z-z_h)_{\partial K} }\right. \Delta {\textstyle\tfrac{1}{2}} \partial_{n}^{} \partial \left.\vphantom{(f+\Delta u_h,z-z_h)_K - \tfrac{1}{2} ([\partial_nu_h],z-z_h)_{\partial K} }\right\}$$ (6)

where  [$\partial_{n}^{}$u_h]  denotes the jump of the normal derivative  $\partial_{n}^{}$u_h  accross the interelement boundaries. From (6), we deduce the a posteriori error bound

$$\le \sum_{K\in\mathbb T_h}^{} \omega_{K}^{} \rho_{K}^{}$$ (7)

with the local residuals $\rho_{K}^{}$ and weights $\omega_{K}^{}$ defined by

$$\rho_{K}^{} \Delta {\frac{1}{2}} \nabla \partial$$

$$\omega_{K}^{} \Big\{ \partial \Big\}$$

Weighted a posteriori error estimates in global norms as the energy norm or the L² norm can obtained within this framework by taking the special error functionals

$$\varphi \nabla \int_{\Omega}^{} \nabla \varphi \nabla \varphi \int_{\Omega}^{} \varphi$$ (8)

which depend on the error  e  itself. In this case, the evaluation of the weight factors  $\omega_{K}^{}$ , i.e., the computation of the dual solution  z , requires to approximate the functional  J( ^. ) . This may be achieved by replacing the unknown error  e  by some approximation  $\tilde{e}$ = $\tilde{u}_{h}^{}$ - u_h  obtained by extrapolation from two consecutively coarser meshes.

In general, the weights  $\omega_{K}^{}$  cannot be determined analytically, but have to be computed by solving the dual problem numerically on the available mesh. To this end one uses the interpolation estimate

$$\omega_{K}^{} \le \nabla^{2}_{}$$ (9)

for z $\in$ H²(K), and evaluates the right hand side by simply taking second order difference quotients of the approximate dual solution $\tilde{z}_{h}^{}$ $\in$ V_h,

$$\omega_{K}^{} \approx \tilde{\omega}_{K}^{} \nabla^{2}_{h} \tilde{z}$$ (10)

where  x_K  is the mid-point of element  K. This results in approximate a posteriori error bounds

$$\approx \eta_{weight}^{} \sum_{K\in\mathbb T_h}^{} \tilde{\omega}_{K}^{} \rho_{K}^{}$$ (11)

It has been demonstrated in [26] that this approximation has only minor effects on the quality of the resulting meshes. The interpolation constant  C_i  may be set equal to one. In our test calculations, we actually set it to  C_i = 0.2 , in order to compensate for the overestimation effect in the resulting approximate a posteriori estimate.

The a posteriori error estimate (11) provides the basis of a feed-back process by which optimally refined meshes can be generated. On the current mesh the dual solution is computed numerically yielding approximate weights. Using the resulting a posteriori error estimate the mesh is refined according, for example, to the strategy described below. This process is repeated yielding more and more accurate weights, i.e., a posteriori error estimates, until the prescribed stopping criterion is fulfilled. This approach allows us to construct almost optimal meshes for various kinds of error measures, where ,,optimal'' can mean ,,most economical for achieving a prescribed accuracy  TOL'' or ,,most accurate for a given maximum number N_max of mesh points''.

The outlined approach to error control for finite element Galerkin schemes by weighted a posteriori error estimates is rather universal as it can, in principle, be applied to almost any problem, linear or nonlinear, which is posed in a variational setting. This has already been demonstrated for simple linear model cases in [26], for the nonlinear Navier-Stokes equations in fluid mechanics in [25], for ordinary differential equations in [27], for one-dimensional combustion problems in [28], for integro-differential equations modelling radiative transfer in astrophysics in [30], and for problemes in linear elasticity in [31].