Adaptive Finite Elements with Thin Obstacles

This work describes concepts for a posteriori error estimation and adaptive mesh design for dual-mixed finite element models where the solution is subjected to inequality constraints. These methods are developed here especially for variational problems with a so-called thin obstacle, i.e. the restrictions are imposed on the boundary of the domain under consideration. The studies presented in this work are partly motivated by our memberships in research groups , where one focus is on the numerical FE-analysis of 3D-contact problems. Eventually, numerical tests confirm our theoretical results.


Introduction
The studies presented in this work are partly motivated by our memberships in some research groups supported by the Deutsche Forschungsgemeinschaft (DFG).In collaboration with scientists from mechanical engineering, we provide the numerical analysis for problems arising in the field of continuum mechanics.
In this context (c.f.Suttmeier [15]), systems of partial differential equations, where the solution is subjected to inequality constraints, have to be treated.We employ the finite element Galerkin (FE) method to obtain approximate solutions of such systems, which for instance typically arise in the field of contact mechanics.Examples are plastic materials, where certain norms of the stresses are bounded or unilateral problems, where the displacement is restricted by a rigid obstacle.The basis for applying an FE discretisation is a suitable mathematical setting, which in the topics under consideration takes the form of variational inequalities (VI).For illustration, the situation of a workpiece pressed onto a grinding disk (Figures 1,2) is approximated by a FE-scheme (Figures 3,4) to approximate the resulting surface forces.
The underlying, mathematical model situation for contact problems in elasticity is Signorini's problem, which in classical notation takes the form (cf. Kikuchi and Oden [8]) This idealised model describes the deformation of an elastic body occupying the domain Ω ⊂ R 3 , which is unilaterally supported by a frictionless rigid foundation.The boundary ∂Ω is assumed to be divided into non-empty subsets Γ N ,Γ D ,Γ C , which have no points in common.The displacement u and the corresponding stress tensor σ are caused by a body force f and a surface traction t along Γ N .Along the portion Γ D of the boundary the body is fixed and Γ C ⊂ ∂Ω denotes the part which is a candidate contact surface.We use the notation where n is the outward normal of ∂Ω, and G denotes the gap between Γ C and the foundation.
Further, the deformation is supposed to be small, so that the strain tensor can be written as ε(u) = 1 2 (∇u + ∇u T ).The compliance tensor A is assumed to be symmetric and positive definite.We intend to derive a posteriori error control techniques for this equation.In order to demonstrate the concept for our method for a posteriori error estimation, we consider the simplified case  Problem (2) is to be solved by the finite element Galerkin method on adaptively optimized meshes.
The basis for applying the finite element method to (2) is the formulation as a variational inequality.Here as first option a solution u ∈ K is sought fulfilling where we set Here and in what follows, (., .)represents the L 2 inner product of a bounded domain Ω in R 2 and .the corresponding norm.Furthermore H m = H m (Ω) denotes the standard Sobolev space of L 2 -functions with derivatives in L 2 (Ω) up to the order m, and H 1 0 ⊂ H 1 is the subspace of H 1 -functions vanishing on ∂Ω.Existence of a unique solution can be established by using the theoretical framework presented e.g. in Glowinski [7].
We will apply the finite element method on decompositions satisfying the usual condition of shape regularity.For ease of mesh refinement and coarsening, hanging nodes are allowed in our implementation.The width of the mesh T h is characterised in terms of a piecewise constant mesh size function h = h(x) > 0, where h T := h |T = diam(T ).Eventually, following Blum & Suttmeier [2], the finite element approximation u h of u in (3) is determined by with K h = V h ∩ K, where V h by standard linear finite elements and for simplicity g is assumed to be polygonal.This finite dimensional problem can be shown to be uniquely solvable following the same line of arguments as in the continuous case.Optimal order a priori error estimates in the energy norm have been given for example in Falk [6] and Brezzi et al. [3].Dobrowolski and Staib [5] show O(h)-convergence in the energy norm without additional assumptions on the structure of the free boundary.Error estimates with respect to the L ∞ -norm have been obtained, e.g., by Nitsche [9] based on a discrete maximum principle.

A mixed setting
In general (c.f.Suttmeier [14]), for applying a finite element method to a given problem, it has to be written in a variational setting, which is accessible for a Galerkin scheme.Additional constraints, as for example the incompressibility condition in the Stokes problem, have to be considered as restrictions.These can be treated by the Lagrangian formalism yielding saddle point problems.
One important application of such mixed systems and the corresponding finite element schemes is the following: For higher order problems one defines auxiliary variables for the derivatives.The equations describing these new quantities are handled as restrictions in the reformulated problem.In this way certain derivatives represented by auxiliary variables are approximated with higher accuracy.Additionally the new formulation often allows to weaken the regularity assumptions on the primary solution.As one example we mention problems in perfect plasticity.Here, discontinuities may occur in the displacement field u due to slip lines in the micro-structure, whereas the stresses σ, determined by linear combinations of derivatives of u, have a smoother behaviour (see Suquet [12], Seregin [11]).From a practical point of view, the stresses are required with high accuracy.Therefore mixed methods treating the stresses directly are often more adequate for solving problems in continuum mechanics.
In our case, introducing the stress vector σ = (σ 1 , σ 2 ), the problem can be stated in form of a first order system, The dual-mixed variational formulation of (5) reads, see e.g.Wang & Yang [16], with Exploiting − div σ = f so called least-squares concepts may be applied to obtain with a(., .):= (., .)+ δ(div ., div .)and suitable δ ≥ 0. The discrete analogue seeks for a pair {σ h , u h } fulfilling The finite dimensional space H h × W h ⊂ H × L 2 for a discretisation for this saddle point problem is determined by approximating each component of the stresses by the standard linear shape functions for H h .W h is constructed by elementwise constant functions.Above the brackets [.] denote the jump across the element boundaries ∂T and δ h is an appropriatly chosen mesh-dependent parameter.For details to this notation we refer to Becker [1].In computations below we choose δ h = O(h)

A posteriori error analysis
In what follows, we present the error analysis for the discretisation error (e, e p ) = (u − u h , p − p h ) between the continuous problem ( 7) and the discrete system (8).Below, the subscript i denotes the standard interpolant of its argument.Additionally defining τ 2 a := τ 2 + δ div τ 2 one obtains Here the last term may be controlled as follows.We start with taking the difference between ( 7) and ( 8) tested with pairs (0, Now locally, using the standard trace theorem, one proceeds by where u l i and u r i denote the values of u i on the two triangles having an edge of ∂T in common.Furthermore we set Next, introducing u L ∈ V h as a suitable approximation of u h with the additional property u L − g ≥ 0 on Γ C , we get Remark: Since u h is a cell-wise constant function, u L may easily be obtained by local averaging techniques (Zienkiewicz & Zhu [17]) and a simple postprocess ensuring u L − g ≥ 0 on Γ C .
Recalling the inequality σn ≥ 0 one obtains Using interpolation gives us with the definition Consequently, observing σ − σ h div = σ − σ h + f + div σ h , we can state Theorem 3.1.For the discretisation error σ − σ h measured in the natural H div -norm there holds the "a posteriori estimate" with notations introduced above.

Numerical test
The numerical results presented throughout this work are obtained by FE-implementations based on the DEAL-library [4,13].
The a posteriori mesh design is organised as follows.Let an error tolerance T OL or a maximal number of cells N max be given.Starting from some initial coarse mesh the refinement criteria are chosen in terms of the local error indicators η T .Then for the mesh refinement, we use the following fixed fraction strategy: In each refinement cycle, the elements are ordered according to the size of η T and then a fixed portion of about 30% of the elements with largest η T is refined resulting in approximately a doubling of the number N of cells.This process is repeated until the stopping criterion η ≤ T OL is fulfilled or N max is exceeded.A more detailed discussion about mesh refinement strategies can be found, e.g., in [10].As a test example we consider the system (5) on the unit square.The external load is f = 0. Zero boundary conditions are imposed on the vertical lines, whereas on the horizontal lines the thin obstacle is defined by g = 100 * (x * x * x − 1.5 * x * x + 0.5 * x) (see Figure 7 (right)).Now based on the error estimate derived above, an adative finite element algorithm is applied to the test problem producing results on a sequence of locally refined triangulations.
Assuming η(σ h , u h , f, g) to behave like O(h α ) we determine an approximation α on each level, by comparing the error bounds with the estimator on the previous level.The results are listed in Table 1.One observes α tending towards 1, indicating the components of the estimate to behave like O(h), as one would anticipate.Furthermore in order to illustrate that in addition the complementarity condition (u − g)σ • n = 0 is well approximated u L − g and σ h • n along one contact boundary are depicted in Figure 7 (left).Eventually a sequence of adaptively refined grids, generated on the basis of the weighted estimate, is

Figure 1 :
Figure 1: Snap shot of a grinding process.

Figure 2 :
Figure 2: Sketch of the model situation.

Figure 6 :
Figure 6: Sequence of locally refined grids produced by our numerical simulation, showing that especially the critical zone between contact zone and free boundary is well resolved.

Figure 7 :
Figure 7: Snapshots along contact line: plot of u L − g and σ h • n (left) demonstrating that the complementarity condition (u − g)σ • n = 0 is well approximated.Plot of the obstacle −g along upper boundary (right).

Table 1 :
Computational results for the test example.Estimated exponents describing rates of convergence O(h α ) of the proposed error estimate.