ANALYSIS OF A SPACE-TIME DISCRETIZATION FOR DYNAMIC ELASTICITY PROBLEMS BASED ON MASS-FREE SURFACE ELEMENTS

. In this paper, a new space-time discretization is proposed which is based on a modiﬁed mass matrix. The mass associated with a surface layer of elements is redistributed such that the inertia at the boundary is removed. This approach is motivated by the observation that standard space-time discretization schemes applied to dynamic contact problems yield spurious oscillations. A widely used approach for the numerical simulation of these problems is based on Lagrange multipliers which represent the contact stresses. But the algebraic contact conditions in combination with the inertia volume terms often yield nonphysical results for the contact stresses, and the stability of the algorithm can be lost. Our modiﬁed matrix is calculated via nonstandard quadrature formulas that require no extra computational eﬀort. In addition, the conservation properties of the underlying algorithm are carried over to the modiﬁed method, and the standard optimal a priori estimates are still satisﬁed. Numerical examples conﬁrm the optimality of the approach and its stabilization eﬀect applied to contact problems.

1. Introduction and motivation.The numerical simulation of contact problems in elasticity plays an important role in many applications in mechanics.The interest in this type of problem has led to extensive research activities from both the numerical and the theoretical points of view (see [7,15,17,20,22] and the references therein for an overview of the topic).A common approach for solving frictional contact problems is the use of Lagrange multipliers [21] to ensure the weak fulfillment of the contact and friction conditions.Mathematically, the Lagrange multipliers represent the dual variable at the contact boundary, whereas physically, they can be identified with the surface traction.Unfortunately, continuum mechanic models for the contact conditions do not capture all features of the dynamic contact properly.As a consequence, numerical simulation results often show spurious oscillations in the Lagrange multiplier (see, e.g., [3,10]).An example is shown in Figure 1.1, where the calculated values for the normal component of the Lagrange multiplier are displayed, computed with a standard finite element discretization scheme in space and a midpoint scheme in time.
Recently, different approaches have discussed how to overcome this numerical artifact [9,11,12,13,14,16].In [12,13,14], a space-time discretization based on a modified mass matrix for the inertia term has been proposed.The mass associated with the elements adjacent to the contact boundary is shifted to the inner elements such that a mass-free layer is obtained.By this, the original PDE decouples into an algebraic equation in time for the contact nodes and a PDE in time for the other nodes, which leads to a higher regularity in time for the contact stresses, as shown in [11].Thereby no spurious oscillations in the numerical computation of the Lagrange multipliers occur.But in the above references, the new mass matrix is defined as the solution of a constrained minimization problem, which makes its computation quite expensive.This drawback is overcome in [9], where an appropriate mass matrix is constructed by local modifications of the original one.In particular, the changes can be assembled in terms of special quadrature formulas, and no extra computational cost comes into play by the modification.
The aim of the current paper is twofold.First, the implementation of the new discretization scheme is studied.Here we propose and discuss two different prototypes of suitable nonstandard quadrature rules.The first formula, already presented in [9], is based on a macroelement decomposition, whereas the second variant can easily be applied to completely unstructured shape-regular triangulations in two and three dimensions.Second, we analyze the error of the new discretization scheme and provide optimal a priori estimates under suitable regularity assumptions by applying techniques from the theory of variational crimes.We emphasize that, in contrast to standard mass lumping techniques, the new method results in a singular mass matrix.Hence, norm equivalences based on the regularity of the mass matrix no longer apply.
The work is organized as follows: In section 2, we state the setting of linear elasticity considered here and introduce our notation.In section 3, we describe two representative examples of quadrature rules that can be used to define the modified mass matrix.We focus on the key requirements that these formulas satisfy.A different mathematical interpretation of these modified quadrature rules as a combination of stable interpolation operators and standard quadrature is presented in section 4. In the next two sections, we prove a priori estimates for the semidiscrete system in section 5 and for the fully discrete system in section 6, after stating the regularity assumptions.Section 7 contains numerical examples confirming the error reduction properties proved in the previous sections.Moreover, we demonstrate the stabilization effect of the method by the numerical solution of a dynamic two-body contact problem.
The boundary ∂Ω is partitioned into two disjoint subsets, ∂Ω= ΓD ∪ ΓN ,w i t h meas (Γ D ) > 0. As we can always reduce nonhomogeneous Dirichlet boundary conditions to the homogeneous case, the value of u on the Dirichlet boundary Γ D is set to zero, whereas the boundary traction on Γ N is given by g N .The outward normal on the boundary ∂Ω is denoted by n,a n dl stands for the volume forces acting on Ω.
The stress tensor σ(u) is given by the constitutive equations of linear elasticity where we assume constant material parameters ̺, μ, λ for ease of notation.
Next we discretize (2.2) in space using a shape-regular triangulation T h of Ω into quadrilaterals/hexahedrals or simplices.The discretization is such that the Dirichlet boundary ΓD can be written as K∈T D K| ∂Ω for a subset T D ⊂T h .
For ω ⊂ R d and k ∈ N 0 , we define P k (ω) as the space of polynomials with the basis functions x α with x ∈ ω, α ∈ N d 0 ,a n d|α|≤k.F u r t h e r m o r e ,Q k (ω)i st h e polynomial space spanned by the functions x α with max 1≤i≤d (α i ) ≤ k.L e tX h be the space of lowest order (k = 1) conforming finite element functions on the triangulation T h ,a n ds e tV h := X h ∩ V.The former space is spanned by the nodal finite element basis functions {φ p } p∈N h ,w h e r eN h is the set of all vertices of the triangulation T h .
With this, we obtain the following semidiscrete problem: Find u h ∈ L 2 ((0,T), V h ) with uh , üh as before and The boundary values u 0h and v 0h are suitable approximations of the continuous initial data which we are going to specify later.Writing the first equation of (2.3) in matrix notation gives We hereby use the notation u h for the function as well as its vector representation.
In the rest of this work, c and C denote generic constants independent of the mesh size h := max K∈T h h K or time variables.Further, we use the abbreviations for the static and • l,k,ω := • H l ((0,T ),H k (ω)) for the timedependent norms.In addition, we set 3. Construction of M H .As motivated in the introduction, we are now going to describe the construction of a new mass matrix M H as a local modification of the standard mass matrix M h such that the degrees of freedom of the potential contact nodes do not contribute to the inertia term.In order to formulate the construction of M H , we need some preliminary definitions.
We assume that contact can occur only on a connected subset Γ C of Γ N satisfying ΓC ∩ ΓD = ∅.Thus, the modification of the mass matrix involves the nodes within the set C h defined by (3.1) As illustrated in Figure 3.1, we denote the union of all elements along Γ C and its complement by The interface ∂Ω C1 ∩ ∂Ω H is called Γ H . Further, we introduce layer-like domains For discrete functions v h , w h ∈ X h , we now define a modified bilinear form with a suitable quadrature rule Q i that is composed of local quadrature formulas Q i K .Two representative examples follow.
Example 3.1 (quadrature rule Q 0 ).On each element K ∈T h , K ⊂ Ω H ,w e take any standard quadrature formula of sufficient accuracy.The only modification is on elements K ⊂ Ω C1 which satisfy by definition K ∩ Γ C = ∅.Depending on the q 1 q 2 q 3 q 1 q 2 q 3 q 4 q 5 q 6 ).On each element K ∈T h , K ⊂ Ω\Ω C2 , we again take any standard quadrature formula of sufficient accuracy.Modifications occur on elements K ⊂ Ω C2 .In order to describe the construction of the quadrature formula, we need the notion of so-called macroelements (see [9] for a more detailed explanation).
We assume that there exists a second triangulation T 1 h possibly with hanging nodes such that each element of T 1 h can be written as the union of elements in T h .Moreover, if K ∈T 1 h with K ⊂ Ω H ,t h e nK ∈T h .E a c hK ∈T 1 h \T h contains at least one element of T h being in Ω C1 (see the shaded parts in Figure 3.4) and exactly one element K 1 ∈T h with K 1 ⊂ Ω H .We note that such triangulation T 1 h exists for any T h as no further conditions on the shape of the macroelements are imposed, but it is not uniquely defined.If T h is obtained from T 2h by uniform refinement, T 1 h can easily be constructed.
In Figure 3.4, representative reference macroelements K in two dimensions are depicted, each with one example of a suitable quadrature formula.The rule for the quadrilateral element (on the left side of Figure 3.4) exhibits a tensor product structure with Gauß nodes in the direction tangential to the contact part of ∂ K and equidistant nodes in the remaining direction; thus, the formula can easily be generalized to the three-dimensional (3D) case.
Using m i H (•, •), i =0 , 1, for the definition of the mass matrix (now denoted by M i H ), we obtain the modified version of (2.4):The solutions of (2.4) and (3.3) are not the same, but from now on u h refers only to the solution of (3.3), also neglecting the index i in order to keep the notation simple.
Several properties are worth noticing about the quadrature formulas presented in Examples 3.1 and 3.2: • They are based on some standard quadrature rule, where the modifications are only locally performed near Γ C .• No quadrature node is located inside the support of basis functions φ p for p ∈C h .T h ee n t r i e sm i H (•, e j φ p )=m i H (e j φ p , •) of the mass matrix M i H with the jth unit vector e j ∈ R d vanish, which leads to zero rows and columns for the degrees of freedom on the boundary Γ C .Hence, the resulting matrix M i H is singular.
• Some of the formulas contain negative weights which in principle could cause stability problems.But as the local mass matrix for K is positive semidefinite for all quadrature rules shown above, the stability is not affected.• Both quadrature formulas are exact for functions in P 0 ( K).In addition, it can easily be checked that Q 1 K is exact for functions in P 2 ( K)f o rt h e simplicial case and for Q 3 ( K) in the quadrilateral case.
• As the definition of Q 1 needs the construction of the macroelement triangulation T 1 h , this rule is more suited for structured meshes which often consist of quadrilateral/hexahedral elements.In contrast, Q 0 can easily be applied even to unstructured simplicial meshes.Examples 3.1 and 3.2 can be considered as prototypes for two different ways of constructing the modified mass matrix M i H ; the former is associated with the triangulation T 0 h := T h , and the latter is based on the macrotriangulation T 1 h .Written concisely, we can define the formula Q i , i ∈{ 0, 1}, locally with respect to the triangulation T i h as follows: Below, we summarize the main features of Examples 3.1 and 3.2 as requirements which a suitable quadrature formula Q i has to satisfy: is exact for all functions in P 2i (K).Clearly, a quadrature formula of type Q 0 satisfying condition Q1 cannot be locally exact for all functions in P 1 (K)w i t hK ∈T 0 h , K ⊂ Ω C1 .But a formula of type Q 1 K can exhibit a higher degree of exactness than P 0 (K), as recorded in condition Q2.
In fact, as stated before, the examples given in Figure 3.4 on the reference macroelement K are exact for all functions in Q 3 ( K) in the quadrilateral/hexahedral case and for P 2 ( K) in the simplicial case.Using the transformation F K : K → K to obtain the nodes and weights on the actual macroelement K ∈T 1 h ,w eh a v et ot a k e into account the determinant det(F −1 K ) which is piecewise constant for simplicial elements but can be piecewise in Q 1 (K) for quadrilateral/hexahedral elements.Hence both quadrature rules given in Figure 3.4 are locally exact for all functions in P 2 (K).
The following lemma is a direct consequence of condition Q2.Lemma 3.3.If the quadrature formula Q i , i ∈{ 0, 1}, is chosen according to condition Q2, the new mass matrix M 0 H in (3.3) conserves the total mass; i.e., . Moreover, the mass matrix M 1 H conserves the zeroth, first, and second order moments of the original system (2.4) (i.e., the total mass, the center of gravity, and the moments of inertia).

Different interpretation of m i H (•, •).
In this section, we will look at the modified bilinear form 2) from a different mathematical point of view.Let us assume that there exists an interpolation operator I i h , i ∈{ 0, 1},o n X h satisfying the three conditions P 1.
Remark 4.1.Condition P 1 is a common assumption for a well-defined interpolation operator, whereas P 2 is motivated by the observation that r e m e n tP 3 reflects the fact that Q i is exact for all functions in P 2i (K), K ∈T i h , i ∈{0, 1}.Remark 4.2.Conditions Q1, P 1, P 2 and the relation (4.1) imply the following inequality (as shown in [9]): Hence we can state that |•| H,i := m i H (•, •)i se q u i v a l e n tt ot h eL 2 (Ω H )-norm and thus is a seminorm on X h .Condition P 1 and (4.1) imply that By condition Q1, all quadrature points of Q i are placed in ΩH , and thus we get by (3.2), (3.4), and P 2 From (4.1), we find (4.3) This formula holds if the quadrature formula h for the product of the functions 4.1.Interpolation operator I 0 h .The action of the linear operator I 0 h can be characterized by the images of the basis functions φ p , p ∈N h .Hence, we compute the right-hand side of (4.3) for functions χ h , η h ∈{e j φ p , 1 ≤ j ≤ d, p ⊂ K ∈T 0 h }.Equation (4.3) and condition Q1implythatI 0 h e j φ p = 0 has to hold for any p ∈C h defined in (3.1).For nodes p ∈N H ,w ek n o wf r o mQ1 that for any modified basis function φ 0 p constructed according to the specification Hence we are looking for scalar values β pq such that and we define (4.6) The first condition of (4.5) is valid for any β pq ∈ R, whereas the second relation depends on the quadrature formula Q 0 K and is fulfilled only for certain choices of β pq .
One can easily verify that the operator I 0 h given by (4.6) satisfies conditions P 1-P 3 and by construction (4.1).

Interpolation operator I 1
h .We proceed analogously to the previous subsection and look for modified local basis functions φ 1 p , p ∈N H , defined according to (4.4) such that A possible construction of the functions φ 1 p , p ∈N H (which are shown in Figure 4.3 for the 2D case), can be done as follows: For each K ∈T 1 h ,l e tK 1 ∈T h be the unique element satisfying K 1 ⊂ K ∩ Ω H (see Figure 4.3).For each vertex p ⊂ K1 , we define φ 1 p as the polynomial extension of h .Now, we can define the interpolation operator I 1 h and the modified finite element space X 1 h as in (4.6) and (4.7), respectively, with the upper index 0 replaced by 1.The conditions P 1-P 3 are easy to verify.
Remark 4.3.We note that the modified basis functions φ i p , i ∈{0, 1},forp ⊂ Γ H are in general not globally continuous because of their elementwise definition, which leads to X i h ⊂ X h .If one wants to avoid such a nonconforming situation or if the mesh is strongly anisotropic, the nodal values of the modified basis functions can be adapted to be continuous or to respect the anisotropy.However, this results in supp φ i p = supp φ p in general.Hence an elementwise definition of the interpolation operator and the corresponding quadrature formula according to (4.3) would not be possible.

5.
A priori error estimate for the semidiscrete system.In order to gain an a priori error estimate for the solution of (3.3), we proceed similarly to the proof of the lumped mass matrix method in [19] which can be interpreted as a variational crime introduced by a quadrature formula (for details see [1]).
In the following, we will need a projection operator onto the discrete space X h which allows for local estimates.An example is the Scott-Zhang operator Z h [18] which preserves the homogeneous Dirichlet boundary conditions on Γ D and satisfies the following for each element K ∈T h : j ∈{0, 1},s∈{1, 2} : where ωK := { T : T ∈T h , T ∩ K = ∅} is the patch of elements around K.
We start with a lemma about the approximation properties of the interpolation operator I i h , i ∈{0, 1}.
Another tool we are going to use in the following is the Sobolev embedding theorem [5] which states that, for an interval A corollary of this theorem is the following lemma given in [6].
Lemma 5.2.Let u ∈ L 2 ((0,T),X) and u ∈ L 2 ((0,T),X ′ ) for some Banach space X.Then we have u ∈ C 0 ([0,T],X) and We define the quadrature error of the modified bilinear form using (4.1): The next lemma gives a bound on this bilinear form.Lemma 5.3.Let u ∈ X, χ h ∈ X h and assume that I i h meets conditions P 1-P 3. Then the following estimate holds: For u, v ∈ H 2 (Ω),w eo b t a i n Proof.Recalling condition P 2 which implies that (η h − I i h η h )v anisheson ΩH for η h ∈ X h , we consider (5.6): w h e r ew em a d eu s eo fP 1, (5.2), (5.3), and (5.4).Now we turn to (5.7).Due to the Sobolev embedding theorem and standard arguments, we get for v ∈ X v 2 0,Ω Cj ≤ Ch v 2 1,Ω ,j ∈{1, 2, 3}.
Here we employed the fact that the width of Ω Cj is bounded in terms of a constant multiple of h.Similarly, we obtain v 2 1,Ω C1 ≤ Ch v 2 2,Ω for v ∈ H 2 (Ω), which leads to We further state another auxiliary result about the time derivative of the discrete operators Z h which holds by definition (see [18]).
Lemma 5.4.For u ∈ X we have Z h u = (Z h u).Now we are able to derive a priori bounds on the difference u h − u.We start with an H 1 (Ω)-estimate which can be shown for I i h satisfying conditions P 1-P 3 without additional assumptions on Ω or the boundary conditions.
To this end, we introduce w h ∈ L 2 ((0,T), V h ) as a suitable elliptic projection of u onto V h : and a V-coercive bilinear form a(•, •) this equation uniquely defines w h .Now we bound the difference between u and w h .
Lemma 5.5.Choose I i h such that conditions P 1-P 3 hold.If the solution u of (2.2) meets the requirements for some k ∈{0, 1, 2}, then the kth time derivative of the solution w h of (5.8) satisfies (5.9) Proof.We give the proof only for k = 0, as the other estimates follow with Lemma 5.4 and differentiation with respect to time.
We start with (2.2) and (5.8) for a test function χ h ∈ L 2 ((0,T), V h ): We cho ose χ h = w h − Z h u; then we get with the V-coercivity of a(•, •), the continuity of m(•, •), integration from 0 to T , and the division of both sides by With Lemma 5.3 and the Cauchy-Schwarz inequality, we obtain In terms of the approximation properties (5.1) and the above estimates, we find With these lemmas, we are able to prove the following result for the H 1 (Ω)-norm of the error.
Theorem 5.6.Assume that the solution u of (2.2) meets the smoothness requirements u, u, ü ∈ L 2 ((0,T), H 2 (Ω)) and ∂ 3 u ∂t 3 , ∂ 4 u ∂t 4 ∈ L 2 ((0,T), X).L e tI i h satisfy P 1-P 3. Take the initial conditions u h (0) = Z h u 0 and uh (0) = Z h v 0 .Then the following estimate holds: Proof.L e tu h , w h be defined by the weak form of (3.3) and (5.8), respectively.Setting θ h := u h − w h gives for any test function We cho ose χ h = θh in (5.12) and thus obtain 1 2 Using the Cauchy-Schwarz and Young inequalities and adding the positive term a(θ h , θ h ) on the right-hand side, we arrive at From (5.13) we conclude with Gronwall's lemma, the seminorm equivalence (4.2), and the coercivity of a(•, •)w i t hac o n s t a n tC(T ) that may depend on T : This inequality holds for all t ∈ (0,T] and hence also for the supremum.Now we have to estimate the terms on the right-hand side.First, we can state , where we made use of (5.1) and Lemma 5.5 for k = 2. Second, we get where we applied Lemmas 5.2 and 5.5.Finally, we estimate similarly If we insert all these results into (5.14),we conclude that (5.17) θh Due to Lemmas 5.2 and 5.5, we can estimate u − w h ∞,1,Ω and u − ẇh ∞,0,Ω H such that we arrive at (5.11).
Next, we turn to the estimate in the L 2 (Ω)-norm, where we need a further assumption, namely, the H 2 (Ω)-regularity of the boundary value problem (2.2).
If we choose v = w h − u in (5.20) and use (5.10), we get for ).Now we take χ h = Z h ϕ and integrate the above equality from 0 to T : For the first term on the right-hand side, we get with the continuity of a(•, •), the approximation property (5.1) of the operator Z h , and Lemma 5.5 The second and third terms are estimated by means of (5.1), (5.3), and Lemma 5.3: Using (5.19), (5.21), and these estimates, we can now conclude that This lemma implies the following estimate.Theorem 5.8.Let the problem (2.2) be H 2 (Ω)-regular and assume that ∂ s u ∂t s ∈ L 2 ((0,T), H 2 (Ω)) for all s ∈{0,...,3}, ∂ 4 u ∂t 4 ∈ L 2 ((0,T), X).Choose I i h according to P 1-P 3 and take the initial conditions u h (0) = Z h u 0 and uh (0) = Z h v 0 .Then the following inequality holds: Proof.Following the proof of [1, Theorem 4.1] and using the seminorm equivalence (4.2) and Lemma 5.7, we find (5.24) Equation ( 5.24) and Theorem 5.6 then conclude the proof.Remark 5.9.The regularity assumptions that we pose on the exact solution u in Theorem 5.8 are a little more strict than those from [1, Theorem 4.1].For Theorem 5.6, we need less regularity for the higher time derivatives of u than [1, Theorem 4.2], but we can bound only the L 2 (Ω H )-norm of the error ( u − uh )i n s t e a d of the norm on the full domain Ω.Still, Theorem 5.8 provides an error estimate of O(h 2 ) with respect to the L 2 (Ω)-norm.
6. Error estimate for the fully discrete system.In this section, we look at the time-discretized version of (3.3):Let Jτ := T for some integer J and the stepwidth τ .S e tw n := w(t n )f o rn ∈ N 0 , w ∈ C( Ω × [0,T]), t n := nτ , and define The fully discrete problem we now consider reads as follows: Find sequences (u n h ) J n=0 , (v n h ) J n=0 in V h , corresponding to the displacement and the velocity, such that The main results are the following two theorems with the constant C possibly depending on T ; the proofs are given in Appendix A for the sake of completeness.
Theorem 6.1.Under the assumptions of Theorem 5.6, we have Let the assumptions of Theorem 5.8 be satisfied as well as ∂ 4 u ∂t 4 ∈ L 2 ((0,T), H 2 (Ω)).Then the following estimate holds: The L 2 (Ω)-error can be bounded as follows.Corollary 6.3.Under the assumptions of Theorems 6.1 and 6.2,w eo b t a i n Proof.The estimate follows from the Poincaré-type estimate (5.25), the triangle inequality, and Theorems 6.1 and 6.2.
7. Numerical results.For the numerical tests, we use the finite element toolbox UG [2].All problems are discretized in time using the midpoint rule as investigated in section 6 which can be regarded as a Hilber-Hughes-Taylor scheme with the parameters 2α =2β = γ = 1 (see [17]).
7.1.Linear beam in two dimensions.First, we validate our error estimates proved in sections 5 and 6 via numerical computation.To this end, we consider the cross-section of an elastic beam of length 10 and height 1.
As material parameters for the linear elasticity law, we use Young's modulus E = 200, Poisson's ratio ν =0 .3, and mass density ̺ =1 .0.We choose the righthand side f and Neumann boundary conditions on ∂Ω such that the exact solution is given by u(x,t)= 0, 0.02 This leads to homogeneous Neumann boundary conditions on the bottom boundary given by x 2 = 0, which we define as Γ C .In the spatial domain, we use a uniform quadrilateral grid with mesh widths from h = 1 4 to h = 1 32 , and for the time discretization we compute 250 time steps with τ =10 −5 .T h el e f ts i d eo fF i g u r e7 .1p r e s e n t st h ee x a c te r r o ra tt i m et 30 with respect to the mesh size h.As the computation with M 1 H and M 0 H leads to almost exactly the same error, only one value is plotted.We obtain an error reduction of O(h 2 )i nt h e L 2 -norm and of O(h)i nt h eH 1 -norm which corresponds to the theoretical results of sections 5 and 6.On the right side of Figure 7.1, the absolute value of the difference between the errors for the modified and the standard computation is plotted.Here we observe numerically that the modified error converges to the standard error with a convergence order of 4 in the H 1 -norm and 5 in the L 2 -norm.Again, there is no visible difference between the results for M 0 H and M 1 H .In Figure 7.2, the vertical displacement and velocity of the node (5.0, 1.0) are depicted.We see that the behavior of the beam is accurately resolved.7.2.Frictionless contact in two dimensions.Further, we investigate the stabilization property of the new discretization scheme applied to contact problems.For this, we simulate the contact of two elastic discs, each with the same radius R =8 and the data E = 100, ν =0 .3,̺ =2 .0• 10 −8 .The initial distance is 0.1, and both circles move toward each other at a speed of 800.We consider only the setting without friction as the frictional work is very small for this example.The grid is refined near the potential contact boundary Γ C as shown in Figure 7 is τ =10 −6 .The numerical treatment of the contact conditions is explained in more detail in [9].
We assemble the mass matrix with two distinct grids: The standard matrix M h and the mass matrix M 0 H are computed for a triangular as well as a quadrilateral mesh; the matrix M 1 H is constructed for the quadrilateral grid only.We use the quadrature formulas described in Examples 3.1 and 3.2 on the appropriate reference elements and transform them into the actual elements.The computed energies for M 0 H (const) and M h (stand) are displayed in Figure 7.4, showing the discrete kinetic energy at time t n given by E n kin := We obtain the same total energy for each computation except for M 0 H on the quadrilateral grid.This is due to the linear determinant of the transformation that is not necessarily integrated exactly (see the discussion before Lemma 3.3).Hence we rather use the quadrature formula Q 1 if we are dealing with a quadrilateral mesh and curvilinear boundaries.The energy results for M 1 H are exactly the same as for the computation with M h , and thus only the standard value is depicted in Figure 7.4..5 shows the improvement in the contact stresses due to the modified mass matrix computed for the simplicial as well as the quadrilateral grid.We see that the standard method exhibits unphysical oscillations, whereas the modification of the mass matrix leads to smoother results.We observe that the beginning and the end of the contact period are the same for all our calculations.
We conclude that solving dynamic contact problems with the modification of the mass matrix gives stable and smooth results without increasing the computational work.
Appendix A. Proofs of section 6.We set θ h := Z h u − w h and further write the following estimate holds with a constant C(T ) depending on T : Proof.Using (4.2) instead of the norm equivalence with respect to • 0,Ω ,t h e proof is done as in [1, Lemma 5.1].
In order to prove the discrete equivalent of Theorem 5.6, we introduce another sequence (ν n h ) J−2 n=0 by Lemma A.2.The following estimate holds: Proof.Using the definitions, we obtain (for details see [1, Proof of Lemma 5.1]) We choose the test function χ h = ∂ t ξ n+1/2 h ∈ V h in (A.10) and arrive at Summing from 0 to (l − 1) for 1 ≤ l ≤ (J − 1) and the use of Schwarz's inequality give for any α>0 Choosing α =4T in (A.11) and using (A.12), we get This holds for each l with 1 ≤ l ≤ (J − 1); hence we can take the maximum of the left-hand side.Subtracting 1 2 A 2 o nb o t hs i d e sa sw e l la su s i n gt h ec o e r c i v i t y ,t h e continuity of a(•, •), and the seminorm equivalence (4.2), we finally arrive at (A.9).This completes the proof.
The estimation of J−1 n=0 ν n h 2 0,Ω H is carried out in the next lemma.Lemma A.3.Let the solution u of (2.2) satisfy ü ∈ L 2 ((0,T), H 2 (Ω)), ∂ 4 u ∂t 4 ∈ L 2 ((0,T), X),a ndletI i h fulfill conditions P 1-P 3. Then the following inequality holds: Proof.By the definition of ν n h (A.8), we get Using Taylor expansion, we obtain for the first term on the right-hand side (see [1, Proof of Lemma 5.2] for details) Hence we arrive at Recalling the definition of θ h , we conclude with Lemma 5.5 and (5.1) that The second term on the right-hand side of (A.14) gives with (A.3) and (A.4) Taylor expansion at t =(n +1)τ , summation over n,a n d( 5 .Combined with (A.18), this yields the desired inequality (6.2).For the estimates in the L 2 (Ω)-norm we again need to impose restrictions as in Theorem 5.8.We begin with bounding the right-hand side of (A.7).
Lemma A.4. Assume that the problem (2.2) is H 2 (Ω)-regular and that ∂ s u ∂t s ∈ L 2 ((0,T), H 2 (Ω)) holds for all s ∈{ 1,...,4}.L e tI i h satisfy conditions P 1-P 3,a n d let u be the solution of (2.2) and (u n h ) J n=0 the solution of (6.1).Then we have the following estimate for the sequence (ǫ n h ) J N =0 defined by (A.Proof.We proceed along the lines of [1, Lemma 5.2], using Lemma 5.7 and the norm equivalence (4.2).
We can now prove Theorem 6.2.

Fig. 1 . 1 .
Fig. 1.1.Normal Lagrange multiplier for a two-body contact problem given in section 7.2 with standard mass matrix M h .

2 .
Problem setting.Let Ω ⊂ R d , d ∈{2, 3}, be an open, bounded domain with piecewise smooth boundary ∂Ω.We consider the following dynamic linear elasticity problem: Find the displacement vector u

Fig. 3 . 4 . 1 K
Fig. 3.4.Nodes and weights for the quadrature rule Q 1 K on macroelements K in two dimensions.

Figures 4. 1
and 4.2 display several representatives of the functions φ 0 p , p ∈N H ,w h i c h correspond to the quadrature rules in Figures 3.2 and 3.3.They span a modified finite element space X 0 h : (4.7)

Fig. 4 . 3 .
Fig. 4.3.Nodal values of the modified basis functions φ 1 p on K for K ⊂ Ω C2 for the 2Dc a s e .
and the total energy E n tot := E n kin + E n pot .The picture at the right of Figure 7.3 visualizes the effective stress σ

Figure 7
Figure 7.5 shows the improvement in the contact stresses due to the modified mass matrix computed for the simplicial as well as the quadrilateral grid.We see that the standard method exhibits unphysical oscillations, whereas the modification of the mass matrix leads to smoother results.We observe that the beginning and the end of the contact period are the same for all our calculations.We conclude that solving dynamic contact problems with the modification of the

..∂ t ξ 0 h 2 0
The second and third terms of (A.19) are estimated by means of (5.3), Lemma 5.2, and Taylor expansion: All together, the initial terms are bounded by