# Analysis of a Space-Time Discretization for Dynamic Elasticity Problems Based on Mass-Free Surface Elements

C. Hager, B. Wohlmuth

► To cite this version:

C. Hager, B. Wohlmuth. Analysis of a Space-Time Discretization for Dynamic Elasticity Problems Based on Mass-Free Surface Elements. SIAM Journal on Numerical Analysis, 2009, 47 (3), pp.1863-1885. 10.1137/080715627. hal-03349902

HAL Id: hal-03349902

https://hal.science/hal-03349902

Submitted on 20 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ANALYSIS OF A SPACE-TIME DISCRETIZATION FOR DYNAMIC
ELASTICITY PROBLEMS BASED ON MASS-FREE
SURFACE ELEMENTS

C. HAGER† AND B. I. WOHLMUTH†

**Abstract.** In this paper, a new space-time discretization is proposed which is based on a modified 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 modified matrix is calculated via nonstandard quadrature formulas that require no extra computational effort. In addition, the conservation properties of the underlying algorithm are carried over to the modified method, and the standard optimal a priori estimates are still satisfied. Numerical examples confirm the optimality of the approach and its stabilization effect applied to contact problems.

**Key words.** dynamical contact problem, mass lumping, quadrature formula, moment conservation, Lagrange multiplier, a priori error estimate

**AMS subject classifications.** 65D32, 65M60, 70E55, 74B05, 74M15, 74S05

**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 [6, 13, 19, 29] 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

<footnote>†Institute of Applied Analysis and Numerical Simulation (IANS), Universität Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany (hager@ians.uni-stuttgart.de, wohlmuth@ians.uni-stuttgart.de).</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1.1. Normal Lagrange multiplier for a two-body contact problem given in section 7.2 with standard mass matrix $M_h$.</figcaption></figure>

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.

**2. Problem setting.** Let $\Omega \subset \mathbb{R}^d$, $d \in \{2, 3\}$, be an open, bounded domain with piecewise smooth boundary $\partial\Omega$. We consider the following dynamic linear elasticity
problem: Find the displacement vector $\mathbf{u} = \mathbf{u}(t, \mathbf{x}) \subset \mathbb{R}^d$ satisfying

$$ (2.1) \qquad \begin{array}{ll} \varrho \ddot{\mathbf{u}} - \operatorname{div} \sigma(\mathbf{u}) = \mathbf{l} & \text{in } (0, T] \times \Omega, \\ \mathbf{u} = \mathbf{0} & \text{in } (0, T] \times \Gamma_D, \\ \sigma(\mathbf{u})\mathbf{n} = \mathbf{g}_N & \text{in } (0, T] \times \Gamma_N, \\ \mathbf{u}(0, \cdot) = \mathbf{u}_0 & \text{in } \Omega, \\ \dot{\mathbf{u}}(0, \cdot) = \mathbf{v}_0 & \text{in } \Omega. \end{array} $$

The boundary $\partial\Omega$ is partitioned into two disjoint subsets, $\partial\Omega = \bar{\Gamma}_D \cup \bar{\Gamma}_N$, with meas $(\Gamma_D) > 0$. As we can always reduce nonhomogeneous Dirichlet boundary conditions to the homogeneous case, the value of $\mathbf{u}$ on the Dirichlet boundary $\Gamma_D$ is set to zero, whereas the boundary traction on $\Gamma_N$ is given by $\mathbf{g}_N$. The outward normal on the boundary $\partial\Omega$ is denoted by $\mathbf{n}$, and $\mathbf{l}$ stands for the volume forces acting on $\Omega$. The stress tensor $\sigma(\mathbf{u})$ is given by the constitutive equations of linear elasticity

$$ \varepsilon(\mathbf{u}) := \frac{1}{2} (\nabla \mathbf{u} + (\nabla \mathbf{u})^T), \quad \sigma(\mathbf{u}) := \lambda \operatorname{tr}(\varepsilon(\mathbf{u}))\mathrm{Id} + 2\mu \varepsilon(\mathbf{u}), $$

where we assume constant material parameters $\varrho, \mu, \lambda$ for ease of notation.

From now on, we denote vector-valued spaces by a bold letter, e.g., $\mathbf{L}_2(\Omega) := [L_2(\Omega)]^d$. To obtain the weak formulation of (2.1), we define $\mathbf{X} := \mathbf{H}^1(\Omega)$ and take $\mathbf{V} := \{\chi \in \mathbf{X} : \chi|_{\Gamma_D} = \mathbf{0}\}$ as the test space for the displacements. Thus, we arrive at the following weak formulation: Find $\mathbf{u} \in L^2((0, T), \mathbf{V})$ such that $\dot{\mathbf{u}} \in L^2((0, T), \mathbf{L}^2(\Omega))$, $\ddot{\mathbf{u}} \in L^2((0, T), \mathbf{V}')$, and

$$ (2.2) \qquad \begin{array}{ll} m(\ddot{\mathbf{u}}, \chi) + a(\mathbf{u}, \chi) = \mathbf{f}(\chi), & \chi \in \mathbf{V}, \quad t \in (0, T], \\ (\mathbf{u}(0, \cdot), \chi) = (\mathbf{u}_0, \chi), & \chi \in \mathbf{V}, \\ (\dot{\mathbf{u}}(0, \cdot), \chi) = (\mathbf{v}_0, \chi), & \chi \in \mathbf{V}, \end{array} $$

with the notation $m(\mathbf{u}, \chi) := \int_{\Omega} \varrho \mathbf{u} \cdot \chi d\mathbf{x}$, $a(\mathbf{u}, \chi) := \int_{\Omega} \sigma(\mathbf{u}) : \varepsilon(\chi) d\mathbf{x}$, $\mathbf{f}(\chi) := \int_{\Omega} \mathbf{l} \cdot \chi d\mathbf{x} + \int_{\Gamma_N} \mathbf{g}_N \cdot \chi ds$, and $(\mathbf{u}, \chi) = \int_{\Omega} \mathbf{u} \cdot \chi d\mathbf{x}$. We remark that by these definitions, $m(\cdot, \cdot)$ and $a(\cdot, \cdot)$ are symmetric and continuous bilinear forms [5]. Further, meas $(\Gamma_D) > 0$ implies that $a(\cdot, \cdot)$ is $\mathbf{V}$-elliptic. The definition of the time-dependent spaces and norms used here can be found in [8].

Next we discretize (2.2) in space using a shape-regular triangulation $\mathcal{T}_h$ of $\Omega$ into quadrilaterals/hexahedrals or simplices. The discretization is such that the Dirichlet boundary $\bar{\Gamma}_D$ can be written as $\bigcup_{K \in \mathcal{T}_D} \bar{K}|_{\partial\Omega}$ for a subset $\mathcal{T}_D \subset \mathcal{T}_h$.

For $\omega \subset \mathbb{R}^d$ and $k \in \mathbb{N}_0$, we define $\mathcal{P}_k(\omega)$ as the space of polynomials with the basis functions $\mathbf{x}^\alpha$ with $\mathbf{x} \in \omega$, $\alpha \in \mathbb{N}_0^d$, and $|\alpha| \le k$. Furthermore, $\mathcal{Q}_k(\omega)$ is the polynomial space spanned by the functions $\mathbf{x}^\alpha$ with $\max_{1 \le i \le d}(\alpha_i) \le k$. Let $\mathbf{X}_h$ be the space of lowest order ($k=1$) conforming finite element functions on the triangulation $\mathcal{T}_h$, and set $\mathbf{V}_h := \mathbf{X}_h \cap \mathbf{V}$. The former space is spanned by the nodal finite element basis functions $\{\phi_p\}_{p \in \mathcal{N}_h}$, where $\mathcal{N}_h$ is the set of all vertices of the triangulation $\mathcal{T}_h$.

With this, we obtain the following semidiscrete problem: Find $\mathbf{u}_h \in L^2((0, T), \mathbf{V}_h)$ with $\dot{\mathbf{u}}_h, \ddot{\mathbf{u}}_h$ as before and

$$ (2.3) \qquad \begin{array}{ll} m(\ddot{\mathbf{u}}_h, \chi_h) + a(\mathbf{u}_h, \chi_h) = \mathbf{f}(\chi_h), & \chi_h \in \mathbf{V}_h, \quad t \in (0, T], \\ (\mathbf{u}_h(0), \chi_h) = (\mathbf{u}_{0h}, \chi_h), & \chi_h \in \mathbf{V}_h, \\ (\dot{\mathbf{u}}_h(0), \chi_h) = (\mathbf{v}_{0h}, \chi_h), & \chi_h \in \mathbf{V}_h. \end{array} $$
The boundary values $\mathbf{u}_{0h}$ and $\mathbf{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

$$ (2.4) \qquad M_h \ddot{\mathbf{u}}_h + A_h \mathbf{u}_h = \mathbf{f}_h. $$

We hereby use the notation $\mathbf{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 \in \mathcal{T}_h} h_K$ or time variables. Further, we use the abbreviations $\| \cdot \|_{k,\omega} := \| \cdot \|_{H^k(\omega)}$ for the static and $\| \cdot \|_{l,k,\omega} := \| \cdot \|_{H^l((0,T),H^k(\omega))}$ for the time-dependent norms. In addition, we set $\| \cdot \|_{\infty,k,\omega} := \| \cdot \|_{L^\infty((0,T),H^k(\omega))}$.

**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 $\Gamma_C$ of $\Gamma_N$ satisfying $\bar{\Gamma}_C \cap \bar{\Gamma}_D = \emptyset$. Thus, the modification of the mass matrix involves the nodes within the set $C_h$ defined by

$$ (3.1) \qquad C_h := \{p \in \mathcal{N}_h : p \in \bar{\Gamma}_C\}, \quad \mathcal{N}_H := \mathcal{N}_h \setminus C_h. $$

As illustrated in Figure 3.1, we denote the union of all elements along $\Gamma_C$ and its complement by

$$ \bar{\Omega}_{C1} := \bigcup_{p \in C_h} (\text{supp } \phi_p), \quad \Omega_H := \Omega \setminus \bar{\Omega}_{C1}. $$

The interface $\partial\Omega_{C1} \cap \partial\Omega_H$ is called $\Gamma_H$. Further, we introduce layer-like domains

$$ \bar{\Omega}_{C2} := \bigcup_{\substack{p \subset \bar{\Omega}_{C1}, \\ p \in \mathcal{N}_h}} (\text{supp } \phi_p), \quad \bar{\Omega}_{C3} := \bigcup_{\substack{p \subset \bar{\Omega}_{C2}, \\ p \in \mathcal{N}_h}} (\text{supp } \phi_p). $$

For discrete functions $\mathbf{v}_h, \mathbf{w}_h \in \mathbf{X}_h$, we now define a modified bilinear form

$$ (3.2) \qquad m_H^i(\mathbf{v}_h, \mathbf{w}_h) := Q^i(\varrho \mathbf{v}_h \cdot \mathbf{w}_h) $$

with a suitable quadrature rule $Q^i$ that is composed of local quadrature formulas $Q_K^i$. Two representative examples follow.

**Example 3.1** (quadrature rule $Q^0$). On each element $K \in \mathcal{T}_h$, $K \subset \Omega_H$, we take any standard quadrature formula of sufficient accuracy. The only modification is on elements $K \subset \Omega_{C1}$ which satisfy by definition $\bar{K} \cap \Gamma_C \neq \emptyset$. Depending on the

<figure><img src="image_3.png" /><figcaption>FIG. 3.1. Examples for $\Omega_{C1}$, $\Omega_{C2}$, $\Omega_{C3}$, and $\Omega_H$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3.2. Nodes and weights of the quadrature rule $Q_{\hat{K}}^0$ in two dimensions for $K \subset \Omega_{C1}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3.3. Nodes and weights of the quadrature rule $Q_{\hat{K}}^0$ in three dimensions for $K \subset \Omega_{C1}$.</figcaption></figure>

shape of the element and its position relative to $\Gamma_C$, different situations can occur. In Figures 3.2 and 3.3, several cases are depicted; for each situation, the corresponding quadrature nodes and weights with respect to the reference element $\hat{K}$ ($[0, 1]^d$ or the unit simplex, respectively) are given. The intersection of $\partial \hat{K}$ with the boundary $\Gamma_C$ is on the bottom side of each element and marked in gray, whereas the quadrature nodes are indicated by filled black bullets.

*Example 3.2 (quadrature rule $Q^1$).* On each element $K \in \mathcal{T}_h$, $K \subset \Omega \setminus \Omega_{C2}$, we again take any standard quadrature formula of sufficient accuracy. Modifications occur on elements $K \subset \Omega_{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 $\mathcal{T}_h^1$ possibly with hanging nodes such that each element of $\mathcal{T}_h^1$ can be written as the union of elements in $\mathcal{T}_h$. Moreover, if $K \in \mathcal{T}_h^1$ with $K \subset \Omega_H$, then $K \in \mathcal{T}_h$. Each $K \in \mathcal{T}_h^1 \setminus \mathcal{T}_h$ contains at least one element of $\mathcal{T}_h$ being in $\Omega_{C1}$ (see the shaded parts in Figure 3.4) and exactly one element $K_1 \in \mathcal{T}_h$ with $K_1 \subset \Omega_H$. We note that such triangulation $\mathcal{T}_h^1$ exists for any $\mathcal{T}_h$ as no further conditions on the shape of the macroelements are imposed, but it is not uniquely defined. If $\mathcal{T}_h$ is obtained from $\mathcal{T}_{2h}$ by uniform refinement, $\mathcal{T}_h^1$ can easily be constructed.

In Figure 3.4, representative reference macroelements $\hat{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 $\partial \hat{K}$ and equidistant nodes in the remaining direction; thus, the formula can easily be generalized to the three-dimensional (3D) case.

Using $m_H^i(\cdot, \cdot)$, $i = 0, 1$, for the definition of the mass matrix (now denoted by $M_H^i$), we obtain the modified version of (2.4):

$$ (3.3) \qquad M_H^i \ddot{u}_h + A_h u_h = f_h. $$
<figure><img src="image_6.png" /></figure>

$$
\begin{align*}
w_1 &= w_2 = 2|\hat{K}|, \\
w_3 &= w_4 = -\frac{15}{4}|\hat{K}|, \\
w_5 &= w_6 = 3|\hat{K}|, \\
w_7 &= w_8 = -\frac{3}{4}|\hat{K}|.
\end{align*}
$$

<figure><img src="image_7.png" /></figure>

$$
\begin{align*}
w_1 &= \frac{1}{9}|\hat{K}|, & w_2 &= \frac{11}{9}|\hat{K}|, \\
w_3 &= \frac{8}{9}|\hat{K}|, & w_4 &= -\frac{5}{9}|\hat{K}|, \\
w_5 &= -\frac{13}{9}|\hat{K}|, & w_6 &= \frac{7}{9}|\hat{K}|.
\end{align*}
$$

<caption>FIG. 3.4. Nodes and weights for the quadrature rule $Q_{\hat{K}}^1$ on macroelements $\hat{K}$ in two dimensions.</caption>

The solutions of (2.4) and (3.3) are not the same, but from now on $\mathbf{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 $\Gamma_C$.

* No quadrature node is located inside the support of basis functions $\phi_p$ for $p \in C_h$. The entries $m_H^i(\cdot, \mathbf{e}_j\phi_p) = m_H^i(\mathbf{e}_j\phi_p, \cdot)$ of the mass matrix $M_H^i$ with the $j$th unit vector $\mathbf{e}_j \in \mathbb{R}^d$ vanish, which leads to zero rows and columns for the degrees of freedom on the boundary $\Gamma_C$. Hence, the resulting matrix $M_H^i$ is singular.

* Some of the formulas contain negative weights which in principle could cause stability problems. But as the local mass matrix for $\hat{K}$ is positive semidefinite for all quadrature rules shown above, the stability is not affected.

* Both quadrature formulas are exact for functions in $\mathcal{P}_0(\hat{K})$. In addition, it can easily be checked that $Q_{\hat{K}}^1$ is exact for functions in $\mathcal{P}_2(\hat{K})$ for the simplicial case and for $\mathcal{Q}_3(\hat{K})$ in the quadrilateral case.

* As the definition of $Q^1$ needs the construction of the macroelement triangulation $\mathcal{T}_h^1$, 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_H^i$; the former is associated with the trian-
gulation $\mathcal{T}_h^0 := \mathcal{T}_h$, and the latter is based on the macrotriangulation $\mathcal{T}_h^1$. Written
concisely, we can define the formula $Q^i$, $i \in \{0,1\}$, locally with respect to the trian-
gulation $\mathcal{T}_h^i$ as follows:

$$
(3.4) \qquad Q^i(\mathbf{v}_h \cdot \mathbf{w}_h) := \sum_{K \in \mathcal{T}_h^i} Q_K^i(\mathbf{v}_h \cdot \mathbf{w}_h), \quad \mathbf{v}_h, \mathbf{w}_h \in \mathbf{X}_h.
$$

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:

Q1. No quadrature point of $Q^i$ is placed in $\bar{\Omega}_{C1} \setminus \bar{\Gamma}_H$.

Q2. For $K \in \mathcal{T}_h^i$, the local formula $Q_K^i$ is exact for all functions in $\mathcal{P}_{2i}(K)$.

Clearly, a quadrature formula of type $Q^0$ satisfying condition $Q1$ cannot be locally
exact for all functions in $\mathcal{P}_1(K)$ with $K \in \mathcal{T}_h^0$, $K \subset \Omega_{C1}$. But a formula of type $Q_K^1$
can exhibit a higher degree of exactness than $\mathcal{P}_0(K)$, as recorded in condition $Q2$.

In fact, as stated before, the examples given in Figure 3.4 on the reference
macroelement $\hat{K}$ are exact for all functions in $\mathcal{Q}_3(\hat{K})$ in the quadrilateral/hexahedral
case and for $\mathcal{P}_2(\hat{K})$ in the simplicial case. Using the transformation $F_K: \hat{K} \to K$ to
obtain the nodes and weights on the actual macroelement $K \in T_h^1$, we have to take into account the determinant det($F_K^{-1}$) 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 \in \{0,1\}$, is chosen according to condition Q2, the new mass matrix $M_H^0$ in (3.3) conserves the total mass; i.e., $1^T M_H^0 1 = 1^T M_h 1$ holds with the vector $1 = (1, \dots, 1)^T \in \mathbb{R}^{\dim(V_h)}$. Moreover, the mass matrix $M_H^1$ 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).

*Proof.* For the proof, see [9]. □

4. Different interpretation of $m_H^i(\cdot, \cdot)$. In this section, we will look at the modified bilinear form $m_H^i(\cdot, \cdot)$ defined in (3.2) from a different mathematical point of view. Let us assume that there exists an interpolation operator $I_h^i$, $i \in \{0,1\}$, on $X_h$ satisfying the three conditions

P1. $I_h^i$ is $L^2(\Omega_{C2})$-stable; i.e., $\|I_h^i \chi_h\|_{0,\Omega_{C2}} \le c\|\chi_h\|_{0,\Omega_{C2}}$, $\chi_h \in X_h$,

P2. $I_h^i \chi_h|_{\bar{\Omega}_H} = \chi_h|_{\bar{\Omega}_H}$, $\chi_h \in X_h$,

P3. $I_h^i \chi_h|_K = \chi_h|_K$, $\chi_h|_K \in P_i(K)$, $K \in T_h^i$,

and

$$ (4.1) \qquad m_H^i(\chi_h, \eta_h) = m(I_h^i \chi_h, I_h^i \eta_h), \quad \chi_h, \eta_h \in X_h. $$

*Remark 4.1.* Condition P1 is a common assumption for a well-defined interpolation operator, whereas P2 is motivated by the observation that $m_H^i(\chi_h, \eta_h) = m(\chi_h, \eta_h)$ for $\chi_h, \eta_h \in X_h$ with supp $\chi_h \subset \bar{\Omega}_H$ or supp $\eta_h \subset \bar{\Omega}_H$. Requirement P3 reflects the fact that $Q^i$ is exact for all functions in $P_{2i}(K)$, $K \in T_h^i$, $i \in \{0,1\}$.

*Remark 4.2.* Conditions Q1, P1, P2 and the relation (4.1) imply the following inequality (as shown in [9]):

$$ (4.2) \qquad c\|\chi_h\|_{0,\Omega_H}^2 \le m_H^i(\chi_h, \chi_h) \le C\|\chi_h\|_{0,\Omega_H}^2, \quad \chi_h \in X_h. $$

Hence we can state that $|\cdot|_{H,i} := \sqrt{m_H^i(\cdot, \cdot)}$ is equivalent to the $L^2(\Omega_H)$-norm and thus is a seminorm on $X_h$. Condition P1 and (4.1) imply that $m_H^i(\cdot, \cdot)$ is continuous on $X_h$.

By condition Q1, all quadrature points of $Q^i$ are placed in $\bar{\Omega}_H$, and thus we get by (3.2), (3.4), and P2

$$ m_H^i(\chi_h, \eta_h) = \sum_{K \in \mathcal{T}_h^i} Q_K^i(\varrho \chi_h \cdot \eta_h) = \sum_{K \in \mathcal{T}_h^i} Q_K^i(\varrho I_h^i \chi_h \cdot I_h^i \eta_h) = m_H^i(I_h^i \chi_h, I_h^i \eta_h). $$

From (4.1), we find

$$ (4.3) \qquad \sum_{K \in \mathcal{T}_h^i} Q_K^i(\varrho I_h^i \chi_h \cdot I_h^i \eta_h) = \sum_{K \in \mathcal{T}_h^i} \int_K \varrho I_h^i \chi_h \cdot I_h^i \eta_h \,dx. $$

This formula holds if the quadrature formula $Q_K^i$ is exact on each element $K \in \mathcal{T}_h^i$ for the product of the functions $I_h^i \chi_h \cdot I_h^i \eta_h$, $\chi_h, \eta_h \in X_h$. Thus, for a given quadrature formula $Q^i$, we can define a corresponding interpolation operator $I_h^i$ via the relation (4.3).

In the following, we return to Examples 3.1 and 3.2 and construct appropriate operators $I_h^0$ and $I_h^1$.
<figure><img src="image_8.png" /><figcaption>FIG. 4.1. Nodal values of the modified basis functions $\phi_p^0$ on $\hat{K}$ for $K \subset \Omega_{C1}$ for the two-dimensional (2D) case.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 4.2. Nodal values of the modified basis functions $\phi_p^0$ on $\hat{K}$ for $K \subset \Omega_{C1}$ for the 3D case.</figcaption></figure>

**4.1. Interpolation operator $I_h^0$.** The action of the linear operator $I_h^0$ can be characterized by the images of the basis functions $\phi_p$, $p \in \mathcal{N}_h$. Hence, we compute the right-hand side of (4.3) for functions $\chi_h, \eta_h \in \{\mathbf{e}_j\phi_p, 1 \le j \le d, p \subset K \in \mathcal{T}_h^0\}$.

Equation (4.3) and condition Q1 imply that $I_h^0 \mathbf{e}_j \phi_p = 0$ has to hold for any $p \in \mathcal{C}_h$ defined in (3.1). For nodes $p \in \mathcal{N}_H$, we know from Q1 that for any modified basis function $\phi_p^0$ constructed according to the specification

$$ (4.4) \qquad \phi_p^0 := \begin{cases} \phi_p & \text{on } K \subset \Omega_H, \\ \phi_p + \sum_{q \in \mathcal{C}_h, q \subset K} \beta_{pq} \phi_q & \text{on } K \subset \Omega_{C1}, \end{cases} $$

with $\beta_{pq} \in \mathbb{R}$, we have

$$ Q_K^0(\mathbf{e}_j \phi_p \cdot \mathbf{e}_k \phi_q) = \delta_{jk} Q_K^0(\phi_p \phi_q) = \delta_{jk} Q_K^0(\phi_p^0 \phi_q^0). $$

Hence we are looking for scalar values $\beta_{pq}$ such that

$$ (4.5) \qquad Q_K^0(\phi_p^0 \phi_q^0) = \begin{cases} 0, & p \in \mathcal{C}_h \text{ or } q \in \mathcal{C}_h, \\ \int_K \phi_p^0 \phi_q^0 d\mathbf{x} & \text{else,} \end{cases} $$

and we define

$$ (4.6) \qquad I_h^0 \mathbf{e}_j \phi_p := \mathbf{e}_j \phi_p^0, \quad 1 \le j \le d, p \in \mathcal{N}_h. $$

The first condition of (4.5) is valid for any $\beta_{pq} \in \mathbb{R}$, whereas the second relation depends on the quadrature formula $Q_K^0$ and is fulfilled only for certain choices of $\beta_{pq}$. Figures 4.1 and 4.2 display several representatives of the functions $\phi_p^0$, $p \in \mathcal{N}_H$, which correspond to the quadrature rules in Figures 3.2 and 3.3. They span a modified finite element space $\mathbf{X}_h^0$:

$$ (4.7) \qquad \mathbf{X}_h^0 := \left\{ \sum_{p \in \mathcal{N}_H} \alpha_p \phi_p^0 : \alpha_p \in \mathbb{R}^d \right\}. $$

One can easily verify that the operator $I_h^0$ given by (4.6) satisfies conditions P1-P3 and by construction (4.1).
<figure><img src="image_10.png" /><figcaption>FIG. 4.3. Nodal values of the modified basis functions $\phi_p^1$ on $\hat{K}$ for $K \subset \Omega_{C2}$ for the 2D case.</figcaption></figure>

**4.2. Interpolation operator $I_h^1$.** We proceed analogously to the previous subsection and look for modified local basis functions $\phi_p^1$, $p \in \mathcal{N}_H$, defined according to (4.4) such that

$$Q_K^1(\mathbf{e}_j\phi_p \cdot \mathbf{e}_k\phi_q) = \delta_{jk}Q_K^1(\phi_p^1\phi_q^1) = \begin{cases} 0, & p \in \mathcal{C}_h \text{ or } q \in \mathcal{C}_h, \\ \delta_{jk} \int_K \phi_p^1\phi_q^1 d\mathbf{x} & \text{else.} \end{cases}$$

A possible construction of the functions $\phi_p^1$, $p \in \mathcal{N}_H$ (which are shown in Figure 4.3 for the 2D case), can be done as follows: For each $K \in \mathcal{T}_h^1$, let $K_1 \in \mathcal{T}_h$ be the unique element satisfying $K_1 \subset K \cap \Omega_H$ (see Figure 4.3). For each vertex $p \subset \bar{K}_1$, we define $\phi_p^1$ as the polynomial extension of $\phi_p|_{K_1}$ onto $K$. This implies $\phi_p^1|_{K'} = \phi_p|_{K'}$ for $K' \in \mathcal{T}_h \cap \mathcal{T}_h^1$.

Now, we can define the interpolation operator $I_h^1$ and the modified finite element space $X_h^1$ as in (4.6) and (4.7), respectively, with the upper index 0 replaced by 1. The conditions P1–P3 are easy to verify.

**Remark 4.3.** We note that the modified basis functions $\phi_p^i$, $i \in \{0, 1\}$, for $p \subset \Gamma_H$ are in general not globally continuous because of their elementwise definition, which leads to $X_h^i \not\subset 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 $\text{supp } \phi_p^i \neq \text{supp } \phi_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 $\Gamma_D$ and satisfies the following for each element $K \in \mathcal{T}_h$:

$$(5.1) \quad j \in \{0,1\}, s \in \{1,2\} : \quad \|Z_h \mathbf{v} - \mathbf{v}\|_{j,K} \le C h^{s-j} |\mathbf{v}|_{s,\omega_K}, \quad \mathbf{v} \in \mathbf{H}^s(\Omega),$$

$$(5.2) \qquad |Z_h \mathbf{v}|_{1,K} \le C |\mathbf{v}|_{1,\omega_K}, \qquad \mathbf{v} \in \mathbf{H}^1(\Omega),$$

$$(5.3) \qquad \|Z_h \mathbf{v}\|_{0,K} \le C \|\mathbf{v}\|_{0,\omega_K}, \qquad \mathbf{v} \in \mathbf{L}^2(\Omega),$$

where $\bar{\omega}_K := \bigcup \{\bar{T} : T \in \mathcal{T}_h, \bar{T} \cap \bar{K} \neq \emptyset\}$ is the patch of elements around $K$.

We start with a lemma about the approximation properties of the interpolation operator $I_h^i$, $i \in \{0, 1\}$.
LEMMA 5.1. Let $I_h^i$ satisfy conditions P1-P3 (given in section 4). Then the following inequalities hold:

$$ (5.4) \qquad \|I_h^i \chi_h - \chi_h\|_{0,\Omega} \le Ch \|\chi_h\|_{1,\Omega_{C2}}, \quad \chi_h \in \mathbf{X}_h, $$

$$ (5.5) \qquad \|I_h^i Z_h \mathbf{v} - \mathbf{v}\|_{0,\Omega_{C1}} \le Ch \|\mathbf{v}\|_{1,\Omega_{C3}}, \quad \mathbf{v} \in \mathbf{X}. $$

*Proof.* The estimate (5.4) follows directly from P2 and P3, whereas (5.5) results from the triangle inequality, (5.1), (5.2), and (5.4). □

Another tool we are going to use in the following is the Sobolev embedding theorem [5] which states that, for an interval $[a, b] \subset \mathbb{R}$, $L^\infty([a, b])$ is continuously embedded in $H^1([a, b])$. A corollary of this theorem is the following lemma given in [6].

LEMMA 5.2. Let $\mathbf{u} \in L^2((0,T), X)$ and $\dot{\mathbf{u}} \in L^2((0,T), X')$ for some Banach space $X$. Then we have $\mathbf{u} \in C^0([0,T], X)$ and

$$ \sup_{t \in [0,T]} \|\mathbf{u}(t)\|_{X} \le C \left( \|\mathbf{u}\|_{L^2((0,T),X)} + \|\dot{\mathbf{u}}\|_{L^2((0,T),X')} \right). $$

We define the quadrature error of the modified bilinear form using (4.1):

$$ \varepsilon_H^i(\chi_h, \eta_h) := m_H^i(\chi_h, \eta_h) - m(\chi_h, \eta_h) = m(I_h^i \chi_h, I_h^i \eta_h) - m(\chi_h, \eta_h) $$

for functions $\chi_h, \eta_h \in \mathbf{X}_h$. The next lemma gives a bound on this bilinear form.

LEMMA 5.3. Let $\mathbf{u} \in \mathbf{X}$, $\chi_h \in \mathbf{X}_h$ and assume that $I_h^i$ meets conditions P1-P3. Then the following estimate holds:

$$ (5.6) \qquad |\varepsilon_H^i(Z_h \mathbf{u}, \chi_h)| \le Ch (\|\mathbf{u}\|_{1,\Omega_{C3}} \|\chi_h\|_{0,\Omega_{C2}} + \|\mathbf{u}\|_{0,\Omega_{C2}} \|\chi_h\|_{1,\Omega_{C2}}). $$

For $\mathbf{u}, \mathbf{v} \in \mathbf{H}^2(\Omega)$, we obtain

$$ (5.7) \qquad |\varepsilon_H^i(Z_h \mathbf{u}, Z_h \mathbf{v})| \le Ch^2 (\|\mathbf{u}\|_{2,\Omega} \|\mathbf{v}\|_{1,\Omega} + \|\mathbf{u}\|_{1,\Omega} \|\mathbf{v}\|_{2,\Omega}). $$

*Proof.* Recalling condition P2 which implies that $(\eta_h - I_h^i \eta_h)$ vanishes on $\bar{\Omega}_H$ for $\eta_h \in \mathbf{X}_h$, we consider (5.6):

$$ \begin{aligned} |\varepsilon_H^i(Z_h \mathbf{u}, \chi_h)| &= |m(I_h^i Z_h \mathbf{u}, I_h^i \chi_h) - m(Z_h \mathbf{u}, \chi_h)| \\ &= |m((I_h^i - \mathrm{Id})Z_h \mathbf{u}, I_h^i \chi_h) + m(Z_h \mathbf{u}, (I_h^i - \mathrm{Id})\chi_h)| \\ &\le C(\|I_h^i - \mathrm{Id}\|Z_h \mathbf{u}\|_{0,\Omega_{C1}}\|I_h^i \chi_h\|_{0,\Omega_{C1}} + \|Z_h \mathbf{u}\|_{0,\Omega_{C1}}\|(I_h^i - \mathrm{Id})\chi_h\|_{0,\Omega_{C1}}) \\ &\le Ch (\|Z_h \mathbf{u}\|_{1,\Omega_{C2}}\|\chi_h\|_{0,\Omega_{C2}} + \|\mathbf{u}\|_{0,\Omega_{C2}}\|\chi_h\|_{1,\Omega_{C2}}) \\ &\le Ch (\|\mathbf{u}\|_{1,\Omega_{C3}}\|\chi_h\|_{0,\Omega_{C2}} + \|\mathbf{u}\|_{0,\Omega_{C2}}\|\chi_h\|_{1,\Omega_{C2}}), \end{aligned} $$

where we made use of P1, (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 $\mathbf{v} \in \mathbf{X}$

$$ \|\mathbf{v}\|_{0,\Omega_{Cj}}^2 \le Ch \|\mathbf{v}\|_{1,\Omega}^2, \quad j \in \{1, 2, 3\}. $$

Here we employed the fact that the width of $\Omega_{Cj}$ is bounded in terms of a constant
multiple of $h$. Similarly, we obtain $\|\mathbf{v}\|_{1,\Omega_{C1}}^2 \le Ch\|\mathbf{v}\|_{2,\Omega}^2$ for $\mathbf{v} \in \mathbf{H}^2(\Omega)$, which leads to

$$
\begin{aligned}
|\varepsilon_H^i(Z_h \mathbf{u}, Z_h \mathbf{v})| & \le \|(I_h^i - \mathrm{Id})Z_h \mathbf{u}\|_{0,\Omega_{C1}} \|(I_h^i Z_h \mathbf{v})\|_{0,\Omega_{C1}} \\
& \quad + \|Z_h \mathbf{u}\|_{0,\Omega_{C1}} \|(I_h^i - \mathrm{Id})Z_h \mathbf{v}\|_{0,\Omega_{C1}} \\
& \le Ch (\|Z_h \mathbf{u}\|_{1,\Omega_{C2}} \|Z_h \mathbf{v}\|_{0,\Omega_{C2}} + \|Z_h \mathbf{u}\|_{0,\Omega_{C1}} \|Z_h \mathbf{v}\|_{1,\Omega_{C2}}) \\
& \le Ch^{3/2} (\|\mathbf{u}\|_{1,\Omega_{C3}} \|Z_h \mathbf{v}\|_{1,\Omega} + \|Z_h \mathbf{u}\|_{1,\Omega} \|\mathbf{v}\|_{1,\Omega_{C3}}) \\
& \le Ch^2 (\|\mathbf{u}\|_{2,\Omega} \|\mathbf{v}\|_{1,\Omega} + \|\mathbf{u}\|_{1,\Omega} \|\mathbf{v}\|_{2,\Omega}). \quad \square
\end{aligned}
$$

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 $\mathbf{u} \in \mathbf{X}$ we have $Z_h \dot{\mathbf{u}} = (Z_h \mathbf{u})$.

Now we are able to derive a priori bounds on the difference $\mathbf{u}_h - \mathbf{u}$. We start with an $\mathbf{H}^1(\Omega)$-estimate which can be shown for $I_h^i$ satisfying conditions P1–P3 without additional assumptions on $\Omega$ or the boundary conditions.

To this end, we introduce $\mathbf{w}_h \in L^2((0, T), \mathbf{V}_h)$ as a suitable elliptic projection of $\mathbf{u}$ onto $\mathbf{V}_h$:

$$ (5.8) \qquad a(\mathbf{w}_h, \chi_h) = (\mathbf{f}, \chi_h) - m_H^i(Z_h \ddot{\mathbf{u}}, \chi_h), \quad \chi_h \in \mathbf{V}_h, \quad 0 \le t \le T. $$

For $\mathbf{V}_h \subset \mathbf{V}$ and a $\mathbf{V}$-coercive bilinear form $a(\cdot, \cdot)$ this equation uniquely defines $\mathbf{w}_h$. Now we bound the difference between $\mathbf{u}$ and $\mathbf{w}_h$.

**LEMMA 5.5.** Choose $I_h^i$ such that conditions P1–P3 hold. If the solution $\mathbf{u}$ of (2.2) meets the requirements $\frac{\partial^k \mathbf{u}}{\partial t^k} \in L^2((0, T), \mathbf{H}^2(\Omega))$, $\frac{\partial^{k+2} \mathbf{u}}{\partial t^{k+2}} \in L^2((0, T), \mathbf{X})$ for some $k \in \{0, 1, 2\}$, then the $k$th time derivative of the solution $\mathbf{w}_h$ of (5.8) satisfies

$$ (5.9) \qquad \left\| \frac{\partial^k \mathbf{w}_h}{\partial t^k} - \frac{\partial^k \mathbf{u}}{\partial t^k} \right\|_{0,1,\Omega} \le Ch \left( \left\| \frac{\partial^k \mathbf{u}}{\partial t^k} \right\|_{0,2,\Omega} + \left\| \frac{\partial^{k+2} \mathbf{u}}{\partial t^{k+2}} \right\|_{0,1,\Omega} \right). $$

*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 $\chi_h \in L^2((0, T), \mathbf{V}_h)$:

$$ (5.10) \qquad \begin{aligned} a(\mathbf{w}_h - Z_h \mathbf{u}, \chi_h) = & a(\mathbf{u} - Z_h \mathbf{u}, \chi_h) + m(\ddot{\mathbf{u}} - Z_h \ddot{\mathbf{u}}, \chi_h) \\ & + m(Z_h \ddot{\mathbf{u}}, \chi_h) - m_H^i(Z_h \ddot{\mathbf{u}}, \chi_h). \end{aligned} $$

We choose $\chi_h = \mathbf{w}_h - Z_h \mathbf{u}$; then we get with the $\mathbf{V}$-coercivity of $a(\cdot, \cdot)$, the continuity of $m(\cdot, \cdot)$, integration from 0 to $T$, and the division of both sides by $\|\mathbf{w}_h - Z_h \mathbf{u}\|_{0,1,\Omega}$

$$ \begin{aligned} \|\mathbf{w}_h - Z_h \mathbf{u}\|_{0,1,\Omega} \le C & \left( \|\ddot{\mathbf{u}} - Z_h \ddot{\mathbf{u}}\|_{0,0,\Omega} + \|\mathbf{u} - Z_h \mathbf{u}\|_{0,1,\Omega} \right) \\ & + \sup_{\substack{\chi_h \in L^2((0,T), \mathbf{V}_h) \\ \chi_h \ne 0}} \frac{C}{\|\chi_h\|_{0,1,\Omega}} \int_0^T |\varepsilon_H(Z_h \ddot{\mathbf{u}}, \chi_h)| dt. \end{aligned} $$

With Lemma 5.3 and the Cauchy–Schwarz inequality, we obtain

$$ \int_0^T |\varepsilon_H(Z_h \ddot{\mathbf{u}}, \chi_h)| dt \le Ch \|\ddot{\mathbf{u}}\|_{0,1,\Omega} \|\chi_h\|_{0,1,\Omega}. $$
In terms of the approximation properties (5.1) and the above estimates, we find

$$ \|\mathbf{w}_h - \mathbf{u}\|_{0,1,\Omega} \le \|\mathbf{w}_h - Z_h \mathbf{u}\|_{0,1,\Omega} + \|Z_h \mathbf{u} - \mathbf{u}\|_{0,1,\Omega} \le Ch (\|\mathbf{u}\|_{0,2,\Omega} + \|\ddot{\mathbf{u}}\|_{0,1,\Omega}). \quad \square $$

With these lemmas, we are able to prove the following result for the $\mathbf{H}^1(\Omega)$-norm of the error.

**THEOREM 5.6.** Assume that the solution $\mathbf{u}$ of (2.2) meets the smoothness requirements $\mathbf{u}, \dot{\mathbf{u}}, \ddot{\mathbf{u}} \in L^2((0,T), \mathbf{H}^2(\Omega))$ and $\frac{\partial^3 \mathbf{u}}{\partial t^3}, \frac{\partial^4 \mathbf{u}}{\partial t^4} \in L^2((0,T), \mathbf{X})$. Let $I_h^i$ satisfy P1-P3. Take the initial conditions $\mathbf{u}_h(0) = Z_h \mathbf{u}_0$ and $\dot{\mathbf{u}}_h(0) = Z_h \mathbf{v}_0$. Then the following estimate holds:

$$ (5.11) \quad \|\dot{\mathbf{u}}_h - \dot{\mathbf{u}}\|_{\infty,0,\Omega_H} + \|\mathbf{u}_h - \mathbf{u}\|_{\infty,1,\Omega} \le Ch \left( \sum_{s=0}^2 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,1,\Omega} \right). $$

*Proof.* Let $\mathbf{u}_h, \mathbf{w}_h$ be defined by the weak form of (3.3) and (5.8), respectively. Setting $\theta_h := \mathbf{u}_h - \mathbf{w}_h$ gives for any test function $\chi_h \in L^2((0,T), \mathbf{V}_h)$

$$ (5.12) \quad m_H^i(\ddot{\theta}_h, \chi_h) + a(\theta_h, \chi_h) = m_H^i(Z_h \ddot{\mathbf{u}}, \chi_h) - m_H^i(\ddot{\mathbf{w}}_h, \chi_h). $$

We choose $\chi_h = \dot{\theta}_h$ in (5.12) and thus obtain

$$ \frac{1}{2} \frac{d}{dt} \left( |\dot{\theta}_h|^2_{H,i} + a(\theta_h, \theta_h) \right) = m_H^i(Z_h \ddot{\mathbf{u}} - \ddot{\mathbf{w}}_h, \dot{\theta}_h). $$

Using the Cauchy-Schwarz and Young inequalities and adding the positive term $a(\theta_h, \theta_h)$ on the right-hand side, we arrive at

$$ (5.13) \quad \frac{d}{dt} \left( |\dot{\theta}_h|^2_{H,i} + a(\theta_h, \theta_h) \right) \le C |Z_h \ddot{\mathbf{u}} - \ddot{\mathbf{w}}_h|^2_{H,i} + \left( |\dot{\theta}_h|^2_{H,i} + a(\theta_h, \theta_h) \right). $$

From (5.13) we conclude with Gronwall's lemma, the seminorm equivalence (4.2), and the coercivity of $a(\cdot, \cdot)$ with a constant $C(T)$ that may depend on $T$:

$$ (5.14) \quad \|\dot{\theta}_h\|_{0,\Omega_H}^2 + \|\theta_h\|_{1,\Omega}^2 \le C \left( |\dot{\theta}_h|^2_{H,i} + a(\theta_h, \theta_h) \right) \\ \le C(T) \left( \|\dot{\theta}_h(0)\|_{0,\Omega_H}^2 + \|\theta_h(0)\|_{1,\Omega}^2 + \|Z_h \ddot{\mathbf{u}} - \ddot{\mathbf{w}}_h\|_{0,0,\Omega_H}^2 \right). $$

This inequality holds for all $t \in (0, T]$ and hence also for the supremum. Now we have to estimate the terms on the right-hand side. First, we can state

$$ \|Z_h \ddot{\mathbf{u}} - \ddot{\mathbf{w}}_h\|_{0,0,\Omega_H} \le \|Z_h \ddot{\mathbf{u}} - \ddot{\mathbf{u}}\|_{0,0,\Omega_H} + \|\ddot{\mathbf{u}} - \ddot{\mathbf{w}}_h\|_{0,0,\Omega_H} \\ \le Ch \left( \|\ddot{\mathbf{u}}\|_{0,2,\Omega} + \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,1,\Omega} \right), $$

where we made use of (5.1) and Lemma 5.5 for $k=2$. Second, we get

$$ \begin{align} \|\dot{\theta}_h(0)\|_{0,\Omega_H} &\le \|Z_h \mathbf{v}_0 - \mathbf{v}_0\|_{0,\Omega_H} + \|\mathbf{v}_0 - \dot{\mathbf{w}}_h(0)\|_{0,\Omega_H} \nonumber \\ &\le C (h \|\mathbf{v}_0\|_{1,\Omega} + \|\dot{\mathbf{u}} - \dot{\mathbf{w}}_h\|_{0,0,\Omega_H} + \|\ddot{\mathbf{u}} - \ddot{\mathbf{w}}_h\|_{0,0,\Omega_H}) \nonumber \\ (5.15) \quad &\le Ch \left( \sum_{s=1}^2 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,1,\Omega} \right), \end{align} $$
where we applied Lemmas 5.2 and 5.5. Finally, we estimate similarly

$$
\begin{align}
\|\boldsymbol{\theta}_h(0)\|_{1,\Omega} &\le \|Z_h \mathbf{u}_0 - \mathbf{u}_0\|_{1,\Omega} + \|\mathbf{u}_0 - \mathbf{w}_h(0)\|_{1,\Omega} \nonumber \\
&\le C (h \|\mathbf{u}_0\|_{2,\Omega} + \|\mathbf{u} - \mathbf{w}_h\|_{0,1,\Omega} + \|\dot{\mathbf{u}} - \dot{\mathbf{w}}_h\|_{0,1,\Omega}) \nonumber \\
(5.16) \qquad &\le Ch \left( \sum_{s=0}^{1} \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + \sum_{s=2}^{3} \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,1,\Omega} \right). \nonumber
\end{align}
$$

If we insert all these results into (5.14), we conclude that

$$
(5.17) \quad \| \dot{\theta}_h \|_{\infty,0,\Omega_H} + \| \theta_h \|_{\infty,1,\Omega} \le C(T)h \left( \sum_{s=0}^{2} \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + \sum_{s=3}^{4} \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,1,\Omega} \right).
$$

Now we recall $\mathbf{u}_h - \mathbf{u} = \boldsymbol{\theta}_h + \mathbf{w}_h - \mathbf{u}$ and get

$$
\| \dot{\mathbf{u}}_h - \dot{\mathbf{u}} \|_{0,\Omega_H} + \| \mathbf{u}_h - \mathbf{u} \|_{1,\Omega} \le \| \dot{\theta}_h \|_{0,\Omega_H} + \| \boldsymbol{\theta}_h \|_{1,\Omega} + \| \dot{\mathbf{u}} - \dot{\mathbf{w}}_h \|_{0,\Omega_H} + \| \mathbf{u} - \mathbf{w}_h \|_{1,\Omega}.
$$

Due to Lemmas 5.2 and 5.5, we can estimate $\|\mathbf{u} - \mathbf{w}_h\|_{\infty,1,\Omega}$ and $\|\dot{\mathbf{u}} - \dot{\mathbf{w}}_h\|_{\infty,0,\Omega_H}$ such that we arrive at (5.11). $\square$

Next, we turn to the estimate in the $L^2(\Omega)$-norm, where we need a further assumption, namely, the $\mathbf{H}^2(\Omega)$-regularity of the boundary value problem (2.2).

LEMMA 5.7. Let problem (2.2) be $\mathbf{H}^2(\Omega)$-regular and assume that $\frac{\partial^k \mathbf{u}}{\partial t^k}, \frac{\partial^{k+2} \mathbf{u}}{\partial t^{k+2}} \in L^2((0,T), \mathbf{H}^2(\Omega))$ for some $k \in \{0,1\}$. Let the interpolation operator $I_h^i$ satisfy conditions P1-P3. Then the following estimate holds:

$$
(5.18) \quad \left\| \frac{\partial^k \mathbf{w}_h}{\partial t^k} - \frac{\partial^k \mathbf{u}}{\partial t^k} \right\|_{0,0,\Omega} \le C h^2 \left( \left\| \frac{\partial^k \mathbf{u}}{\partial t^k} \right\|_{0,2,\Omega} + \left\| \frac{\partial^{k+2} \mathbf{u}}{\partial t^{k+2}} \right\|_{0,2,\Omega} \right).
$$

Proof. Again, it suffices to consider the case $k=0$. We start with the observation that for $\psi \in L^2((0,T), \mathbf{L}^2(\Omega))$, we get

$$
(5.19) \quad \|\psi\|_{0,0,\Omega} = \sup_{\mathbf{g} \in L^2((0,T), \mathbf{L}^2(\Omega)) \setminus \{0\}} \frac{|\int_0^T (\psi, \mathbf{g}) dt|}{\|\mathbf{g}\|_{0,0,\Omega}},
$$

as $L^2((0,T), \mathbf{L}^2(\Omega))$ is its own dual space. Now we choose $\varphi \in L^2((0,T), \mathbf{V})$ as the solution of the following auxiliary problem for any fixed $\mathbf{g} \in L^2((0,T), \mathbf{L}^2(\Omega))$:

$$
(5.20) \qquad a(\varphi, \mathbf{v}) = (\mathbf{g}, \mathbf{v}), \quad \mathbf{v} \in \mathbf{V}, \quad 0 \le t \le T.
$$

By our assumptions on the regularity of the boundary value problem (5.20), we obtain
$\varphi \in L^2((0,T), \mathbf{V} \cap \mathbf{H}^2(\Omega))$ (see [4]) and

$$
(5.21) \qquad \|\varphi\|_{0,2,\Omega} \le C \|\mathbf{g}\|_{0,0,\Omega}.
$$

If we choose $\mathbf{v} = \mathbf{w}_h - \mathbf{u}$ in (5.20) and use (5.10), we get for $\chi_h \in L^2((0,T), \mathbf{V}_h)$

$$
(\mathbf{w}_h - \mathbf{u}, \mathbf{g}) = a(\mathbf{w}_h - \mathbf{u}, \varphi - \chi_h) + m(\ddot{\mathbf{u}} - Z_h \ddot{\mathbf{u}}, \chi_h) + \varepsilon_H^i (Z_h \ddot{\mathbf{u}}, \chi_h).
$$

Now we take $\chi_h = Z_h \varphi$ and integrate the above equality from 0 to T:

$$
\begin{equation} \tag{5.22}
\begin{aligned}
\left| \int_0^T (\mathbf{w}_h - \mathbf{u}, \mathbf{g}) dt \right| &\le \int_0^T |a(\mathbf{w}_h - \mathbf{u}, \varphi - Z_h \varphi)| dt + \int_0^T |m(\ddot{\mathbf{u}} - Z_h \ddot{\mathbf{u}}, Z_h \varphi)| dt \\
&\quad + \int_0^T |\varepsilon_H^i (Z_h \ddot{\mathbf{u}}, Z_h \varphi)| dt.
\end{aligned}
\end{equation}
$$
For the first term on the right-hand side, we get with the continuity of $a(\cdot, \cdot)$, the approximation property (5.1) of the operator $Z_h$, and Lemma 5.5

$$ \int_0^T |a(\mathbf{w}_h - \mathbf{u}, \varphi - Z_h \varphi)|dt \leq Ch^2 (\|\mathbf{u}\|_{0,2,\Omega} + \|\ddot{\mathbf{u}}\|_{0,1,\Omega}) \|\varphi\|_{0,2,\Omega}. $$

The second and third terms are estimated by means of (5.1), (5.3), and Lemma 5.3:

$$ \int_0^T |m(\ddot{\mathbf{u}} - Z_h \ddot{\mathbf{u}}, Z_h \varphi)|dt + \int_0^T |\varepsilon_H^i(Z_h \ddot{\mathbf{u}}, Z_h \varphi)|dt \leq Ch^2 \|\ddot{\mathbf{u}}\|_{0,2,\Omega} \|\varphi\|_{0,2,\Omega}. $$

Using (5.19), (5.21), and these estimates, we can now conclude that

$$ \begin{aligned} \|\mathbf{w}_h - \mathbf{u}\|_{0,0,\Omega} &= \sup_{\mathbf{g} \in L^2((0,T), \mathbf{L}^2(\Omega)) \setminus \{0\}} \frac{|\int_0^T (\mathbf{w}_h - \mathbf{u}, \mathbf{g}) dt|}{\|\mathbf{g}\|_{0,0,\Omega}} \\ &\leq Ch^2 (\|\mathbf{u}\|_{0,2,\Omega} + \|\ddot{\mathbf{u}}\|_{0,2,\Omega}). \quad \square \end{aligned} $$

This lemma implies the following estimate.

**THEOREM 5.8.** Let the problem (2.2) be $\mathbf{H}^2(\Omega)$-regular and assume that $\frac{\partial^s \mathbf{u}}{\partial t^s} \in L^2((0,T), \mathbf{H}^2(\Omega))$ for all $s \in \{0, \dots, 3\}$, $\frac{\partial^4 \mathbf{u}}{\partial t^4} \in L^2((0,T), \mathbf{X})$. Choose $I_h^i$ according to P1-P3 and take the initial conditions $\mathbf{u}_h(0) = Z_h \mathbf{u}_0$ and $\dot{\mathbf{u}}_h(0) = Z_h \mathbf{v}_0$. Then the following inequality holds:

$$ (5.23) \qquad \|\mathbf{u}_h - \mathbf{u}\|_{\infty,0,\Omega} \leq Ch^2 \left( \sum_{s=0}^3 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,1,\Omega} \right). $$

*Proof.* Following the proof of [1, Theorem 4.1] and using the seminorm equivalence (4.2) and Lemma 5.7, we find

$$ (5.24) \qquad \|\mathbf{u}_h - \mathbf{u}\|_{\infty,0,\Omega_H} \leq Ch^2 \sum_{s=0}^3 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega}. $$

A Poincaré-type estimate yields for $\mathbf{v} \in \mathbf{X}$

$$ (5.25) \qquad \|\mathbf{v}\|_{0,\Omega_{C1}} \leq C(\|\mathbf{v}\|_{0,\Omega_H} + h\|\mathbf{v}\|_{1,\Omega}). $$

Equation (5.24) and Theorem 5.6 then conclude the proof. $\square$

**Remark 5.9.** The regularity assumptions that we pose on the exact solution $\mathbf{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 $\mathbf{u}$ than [1, Theorem 4.2], but we can bound only the $\mathbf{L}^2(\Omega_H)$-norm of the error $(\ddot{\mathbf{u}} - \dot{\mathbf{u}}_h)$ instead of the norm on the full domain $\Omega$. Still, Theorem 5.8 provides an error estimate of $O(h^2)$ with respect to the $L_2(\Omega)$-norm.

**6. Error estimate for the fully discrete system.** In this section, we look at the time-discretized version of (3.3): Let $J\tau := T$ for some integer $J$ and the stepwidth $\tau$. Set $\mathbf{w}^n := \mathbf{w}(t_n)$ for $n \in \mathbb{N}_0$, $\mathbf{w} \in C(\bar{\Omega} \times [0, T])$, $t_n := n\tau$, and define

$$ \partial_t \mathbf{w}^n := \frac{1}{\tau}(\mathbf{w}^{n+1} - \mathbf{w}^n), \qquad \mathbf{w}^{n+1/2} := \frac{1}{2}(\mathbf{w}^{n+1} + \mathbf{w}^n). $$
The fully discrete problem we now consider reads as follows: Find sequences $(\mathbf{u}_h^n)_{n=0}^J$, $(\mathbf{v}_h^n)_{n=0}^J$ in $\mathbf{V}_h$, corresponding to the displacement and the velocity, such that

$$ 
\begin{array}{l}
\mathbf{u}_h^0 = Z_h \mathbf{u}_0, \\
(6.1) \qquad \mathbf{v}_h^0 = Z_h \mathbf{v}_0, \\
\qquad m_H^i(\partial_t \mathbf{v}_h^n, \mathbf{\chi}_h) + a(\mathbf{u}_h^{n+1/2}, \mathbf{\chi}_h) = (\mathbf{f}^{n+1/2}, \mathbf{\chi}_h), \quad \mathbf{\chi}_h \in \mathbf{V}_h, \quad 0 \le n \le J-1, \\
\qquad \partial_t \mathbf{u}_h^n = \mathbf{v}_h^{n+1/2}.
\end{array}
$$

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

$$ 
\begin{array}{l}
(6.2) \qquad \max_{0 \le n \le J-1} \left( \|\partial_t (\mathbf{u}_h^n - \mathbf{u}^n)\|_{0,\Omega_H} + \|\mathbf{u}_h^{n+1/2} - \mathbf{u}^{n+1/2}\|_{1,\Omega} \right) \\
\qquad \le C \left( h \left( \sum_{s=0}^2 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,1,\Omega} \right) + \tau^2 \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,0,\Omega} \right).
\end{array}
$$

**THEOREM 6.2.** Let the assumptions of Theorem 5.8 be satisfied as well as $\frac{\partial^4 \mathbf{u}}{\partial t^4} \in L^2((0,T), \mathbf{H}^2(\Omega))$. Then the following estimate holds:

$$ (6.3) \qquad \max_{0 \le n \le J} \|\mathbf{u}_h^n - \mathbf{u}^n\|_{0,\Omega_H} \le C \left( h^2 \sum_{s=0}^3 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + h^2 \tau^2 \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,2,\Omega} + \tau^2 \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,0,\Omega} \right). $$

The $\mathbf{L}_2(\Omega)$-error can be bounded as follows.

**COROLLARY 6.3.** Under the assumptions of Theorems 6.1 and 6.2, we obtain

$$ (6.4) \qquad \max_{0 \le n \le J-1} \|\mathbf{u}_h^{n+1/2} - \mathbf{u}^{n+1/2}\|_{0,\Omega} \\
\le C \left( h^2 \sum_{s=0}^3 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + h^2 \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,1,\Omega} + h^2 \tau^2 \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,2,\Omega} + \tau^2 \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,0,\Omega} \right). $$

*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 tool-box 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\alpha = 2\beta = \gamma = 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 $\nu = 0.3$, and mass density $\rho = 1.0$. We choose the right-hand side $\mathbf{f}$ and Neumann boundary conditions on $\partial\Omega$ such that the exact solution is given by

$$ \mathbf{u}(\mathbf{x}, t) = (0, \ 0.02 x_1 \cdot (10 - x_1) \cdot x_2^2 \cdot \sin(2\pi t))^T. $$
This leads to homogeneous Neumann boundary conditions on the bottom boundary given by $x_2 = 0$, which we define as $\Gamma_C$. In the spatial domain, we use a uniform quadrilateral grid with mesh widths from $h = \frac{1}{4}$ to $h = \frac{1}{32}$, and for the time discretization we compute 250 time steps with $\tau = 10^{-5}$.

<figure><img src="image_11.png" /><figcaption>FIG. 7.1. Error reduction with respect to the mesh size $h$.</figcaption></figure>

The left side of Figure 7.1 presents the exact error at time $t_{30}$ with respect to the mesh size $h$. As the computation with $M_H^1$ and $M_H^0$ leads to almost exactly the same error, only one value is plotted. We obtain an error reduction of $O(h^2)$ in the $L^2$-norm and of $O(h)$ in the $H^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_H^0$ and $M_H^1$.

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.

<figure><img src="image_12.png" /><figcaption>FIG. 7.2. Vertical displacement and velocity of the node (5.0, 1.0) for $h = \frac{1}{8}$ and $M_h$, $M_H^0$.</figcaption></figure>

**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$, $\nu = 0.3$, $\rho = 2.0 \cdot 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 $\Gamma_C$ as shown in Figure 7.3; the size of the time step
<figure><img src="image_13.png" /><figcaption>FIG. 7.3. Initial grids and effective stresses for the two-circle contact problem without friction.</figcaption></figure>

is $\tau = 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_H^0$ are computed for a triangular as well as a quadrilateral mesh; the matrix $M_H^1$ 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_H^0$ (const) and $M_h$ (stand) are displayed in Figure 7.4, showing the discrete kinetic energy at time $t_n$ given by $\mathcal{E}_{\text{kin}}^n := \frac{1}{2}(\mathbf{v}_h^n)^T M \mathbf{v}_h^n$, $M \in \{M_h, M_H^0\}$, the potential energy $\mathcal{E}_{\text{pot}}^n := \frac{1}{2}(\mathbf{u}_h^n)^T A_h \mathbf{u}_h^n$, and the total energy $\mathcal{E}_{\text{tot}}^n := \mathcal{E}_{\text{kin}}^n + \mathcal{E}_{\text{pot}}^n$. The picture at the right of Figure 7.3 visualizes the effective stress $\sigma_{\text{eff}} := \sum_{i,j=1}^d |\sigma_{ij} - \delta_{ij} \frac{1}{d} \text{tr}(\sigma)|^2$.

We obtain the same total energy for each computation except for $M_H^0$ 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_H^1$ are exactly the same as for the computation with $M_h$, and thus only the standard value is depicted in Figure 7.4.

<figure><img src="image_14.png" /><figcaption>FIG. 7.4. Energy results for the two-circle contact problem without friction.</figcaption></figure>

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
<figure><img src="image_15.png" /><figcaption>FIG. 7.5. LM for simplicial and quadrilateral grid at the bottom slave node.</figcaption></figure>

mass matrix gives stable and smooth results without increasing the computational work.

**Appendix A. Proofs of section 6.** We set $\theta_h := Z_h u - w_h$ and further write

$$ (A.1) \qquad \xi_h^n := u_h^n - w_h^n, \quad 0 \le n \le J, $$

$$ (A.2) \qquad \mathbf{p}_h^n := \mathbf{v}_h^n - \dot{\mathbf{v}}_h^n, \quad 0 \le n \le J. $$

**LEMMA A.1.** If we define sequences $(\rho_h^n)_{n=0}^{J-1}$, $(\pi_h^n)_{n=0}^{J-1}$, and $(\epsilon_h^n)_{n=0}^{J-1}$ in $V_h$ by

$$ (A.3) \qquad \rho_h^n := Z_h \ddot{u}^{n+1/2} - \partial_t (Z_h \dot{u}^n), $$

$$ (A.4) \qquad \pi_h^n := \partial_t (Z_h u^n) - Z_h \dot{u}^{n+1/2}, $$

$$ (A.5) \qquad \epsilon_h^0 := \partial_t \theta_h^0 - \pi_h^0 + \frac{\tau}{2} \rho_h^0, $$

$$ (A.6) \qquad \epsilon_h^n := \partial_t \theta_h^n - \pi_h^n + \frac{\tau}{2} \left( \sum_{k=0}^{n} \rho_h^k + \sum_{k=0}^{n-1} \rho_h^k \right), \quad 1 \le n \le J-1, $$

the following estimate holds with a constant $C(T)$ depending on $T$:

$$ (A.7) \qquad \max_{0 \le n \le J} \|\xi_h^n\|_{0, \Omega_H}^2 \le C(T) \left( \tau \sum_{n=0}^{J-1} \|\epsilon_h^n\|_{0, \Omega_H}^2 + \|\xi_h^0\|_{0, \Omega_H}^2 \right). $$

*Proof.* Using (4.2) instead of the norm equivalence with respect to $\| \cdot \|_{0,\Omega}$, the proof is done as in [1, Lemma 5.1]. □

In order to prove the discrete equivalent of Theorem 5.6, we introduce another sequence $(\nu_h^n)_{n=0}^{J-2}$ by

$$ (A.8) \qquad \nu_h^n := \partial_t (\partial_t \theta_h^n) - \partial_t \pi_h^n + \rho_h^{n+1/2}. $$

**LEMMA A.2.** *The following estimate holds:*

$$ (A.9) \qquad \begin{aligned} & \max_{0 \le n \le (J-1)} \left( \|\partial_t \xi_h^n\|_{0, \Omega_H}^2 + \|\xi_h^{n+1/2}\|_{1, \Omega}^2 \right) \\ & \le C(T) \left( \tau \sum_{n=0}^{J-2} \|\nu_h^n\|_{0, \Omega_H}^2 + \|\partial_t \xi_h^0\|_{0, \Omega_H}^2 + \|\xi_h^{1/2}\|_{1, \Omega}^2 \right). \end{aligned} $$
Proof. Using the definitions, we obtain (for details see [1, Proof of Lemma 5.1])

$$ (A.10) \qquad m_H^i (\partial_t (\partial_t \xi_h^n), \chi_h) + \frac{1}{2} a \left( \xi_h^{n+1/2} + \xi_h^{n+3/2}, \chi_h \right) = m_H^i (\nu_h^n, \chi_h). $$

We choose the test function $\chi_h = \partial_t \xi_h^{n+1/2} \in V_h$ in (A.10) and arrive at

$$ \begin{aligned} & |\partial_t \xi_h^{n+1}|_{H,i}^2 - |\partial_t \xi_h^n|_{H,i}^2 + a(\xi_h^{n+3/2}, \xi_h^{n+3/2}) - a(\xi_h^{n+1/2}, \xi_h^{n+1/2}) \\ & \le 2\tau m_H^i (\nu_h^n, \partial_t \xi_h^{n+1/2}). \end{aligned} $$

Summing from 0 to $(l-1)$ for $1 \le l \le (J-1)$ and the use of Schwarz's inequality give for any $\alpha > 0$

$$ (A.11) \qquad |\partial_t \xi_h^l|_{H,i}^2 - |\partial_t \xi_h^0|_{H,i}^2 + a(\xi_h^{l+1/2}, \xi_h^{l+1/2}) - a(\xi_h^{1/2}, \xi_h^{1/2}) \le 2\tau \left( \alpha \sum_{n=0}^{l-1} |\nu_h^n|_{H,i}^2 + \frac{1}{\alpha} \sum_{n=0}^{l-1} |\partial_t \xi_h^{n+1/2}|_{H,i}^2 \right). $$

We set $A := \max_{0 \le n \le (J-1)} |\partial_t \xi_h^n|_{H,i}$, which leads to

$$ (A.12) \qquad \sum_{n=0}^{l-1} |\partial_t \xi_h^{n+1/2}|_{H,i}^2 \le A^2 J. $$

Choosing $\alpha = 4T$ in (A.11) and using (A.12), we get

$$ |\partial_t \xi_h^l|_{H,i}^2 + a(\xi_h^{l+1/2}, \xi_h^{l+1/2}) \le 8\tau T \sum_{n=0}^{J-2} |\nu_h^n|_{H,i}^2 + \frac{1}{2}A^2 + |\partial_t \xi_h^0|_{H,i}^2 + a(\xi_h^{1/2}, \xi_h^{1/2}). $$

This holds for each $l$ with $1 \le l \le (J-1)$; hence we can take the maximum of the left-hand side. Subtracting $\frac{1}{2}A^2$ on both sides as well as using the coercivity, the continuity of $a(\cdot, \cdot)$, and the seminorm equivalence (4.2), we finally arrive at (A.9). This completes the proof. $\square$

The estimation of $\sum_{n=0}^{J-1} \|\nu_h^n\|_{0,\Omega_H}^2$ is carried out in the next lemma.

LEMMA A.3. Let the solution $\mathbf{u}$ of (2.2) satisfy $\ddot{\mathbf{u}} \in L^2((0,T), \mathbf{H}^2(\Omega))$, $\frac{\partial^4 \mathbf{u}}{\partial t^4} \in L^2((0,T), \mathbf{X})$, and let $I_h^i$ fulfill conditions P1-P3. Then the following inequality holds:

$$ (A.13) \qquad \tau \sum_{n=0}^{J-2} \|\nu_h^n\|_{0,\Omega_H}^2 \le C \left( h^2 \left( \|\ddot{\mathbf{u}}\|_{0,2,\Omega}^2 + \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,1,\Omega}^2 \right) + \tau^4 \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,0,\Omega}^2 \right). $$

Proof. By the definition of $\nu_h^n$ (A.8), we get

$$ (A.14) \qquad \|\nu_h^n\|_{0,\Omega_H}^2 \le \|\partial_t(\partial_t \theta_h^n)\|_{0,\Omega_H}^2 + \|\rho_h^{n+1/2} - \partial_t \pi_h^n\|_{0,\Omega_H}^2. $$

Using Taylor expansion, we obtain for the first term on the right-hand side (see [1, Proof of Lemma 5.2] for details)

$$ \begin{aligned} \partial_t (\partial_t \theta_h^n) &= \frac{1}{\tau^2} (\theta_h^{n+2} - 2\theta_h^{n+1} + \theta_h^n) \\ &= \frac{1}{\tau^2} \left( \int_{(n+1)\tau}^{(n+2)\tau} ((n+2)\tau - t) \ddot{\theta}_h(t) dt + \int_{n\tau}^{(n+1)\tau} (t - n\tau) \ddot{\theta}_h(t) dt \right). \end{aligned} $$
Hence we arrive at

$$ \|\partial_t(\partial_t \boldsymbol{\theta}_h^n)\|_{0,\Omega_H}^2 \le C \frac{1}{\tau} \int_{n\tau}^{(n+2)\tau} \|\ddot{\boldsymbol{\theta}}_h(t)\|_{0,\Omega_H}^2 dt. $$

Recalling the definition of $\boldsymbol{\theta}_h$, we conclude with Lemma 5.5 and (5.1) that

$$ (A.15) \qquad \tau \sum_{n=0}^{J-2} \|\partial_t(\partial_t \boldsymbol{\theta}_h^n)\|_{0,\Omega_H}^2 \le Ch^2 \left( \|\ddot{\boldsymbol{u}}\|_{0,2,\Omega}^2 + \left\| \frac{\partial^4 \boldsymbol{u}}{\partial t^4} \right\|_{0,1,\Omega}^2 \right). $$

The second term on the right-hand side of (A.14) gives with (A.3) and (A.4)

$$ \boldsymbol{\rho}_h^{n+1/2} - \partial_t \boldsymbol{\pi}_h^n = Z_h \left( \frac{1}{4} \ddot{\boldsymbol{u}}^{n+2} + \frac{1}{2} \ddot{\boldsymbol{u}}^{n+1} + \frac{1}{4} \ddot{\boldsymbol{u}}^n - \frac{1}{\tau^2} \boldsymbol{u}^{n+2} + \frac{2}{\tau^2} \boldsymbol{u}^{n+1} - \frac{1}{\tau^2} \boldsymbol{u}^n \right). $$

Taylor expansion at $t = (n + 1)\tau$, summation over $n$, and (5.3) lead to

$$ (A.16) \qquad \tau \sum_{n=0}^{J-2} \| \boldsymbol{\rho}_h^{n+1/2} - \partial_t \boldsymbol{\pi}_h^n \|_{0,\Omega_H}^2 \le C \tau^4 \left\| \frac{\partial^4 \boldsymbol{u}}{\partial t^4} \right\|_{0,0,\Omega}^2. $$

Inequalities (A.15) and (A.16) conclude the proof. □

Now we are able to prove Theorem 6.1.

*Proof.* With $\boldsymbol{u}_h^n - \boldsymbol{u}^n = \boldsymbol{\xi}_h^n + (\boldsymbol{w}_h^n - \boldsymbol{u}^n)$ we can estimate

$$ (A.17) \qquad \begin{aligned} & \max_{0 \le n \le J-1} \left( \|\partial_t (\boldsymbol{u}_h^n - \boldsymbol{u}^n)\|_{0,\Omega_H} + \|\boldsymbol{u}_h^{n+1/2} - \boldsymbol{u}^{n+1/2}\|_{1,\Omega} \right) \\ & \le \max_{0 \le n \le J-1} \left( \|\partial_t (\boldsymbol{w}_h^n - \boldsymbol{u}^n)\|_{0,\Omega_H} + \|\boldsymbol{w}_h^{n+1/2} - \boldsymbol{u}^{n+1/2}\|_{1,\Omega} \right. \\ & \qquad \left. + \|\partial_t \boldsymbol{\xi}_h^n\|_{0,\Omega_H} + \|\boldsymbol{\xi}_h^{n+1/2}\|_{1,\Omega} \right). \end{aligned} $$

Using Taylor expansion as well as Lemmas 5.2 and 5.5, the first term gives

$$ \begin{aligned} \|\partial_t (\boldsymbol{w}_h^n - \boldsymbol{u}^n)\|_{0,\Omega_H}^2 &\le \|\partial_t (\boldsymbol{w}_h^n - \boldsymbol{u}^n) - (\dot{\boldsymbol{w}}_h - \dot{\boldsymbol{u}})^{n+1/2}\|_{0,\Omega_H}^2 + \|(\dot{\boldsymbol{w}}_h - \dot{\boldsymbol{u}})^{n+1/2}\|_{0,\Omega_H}^2 \\ &\le \tau \int_0^\tau \|\ddot{\boldsymbol{w}}_h(t) - \ddot{\boldsymbol{u}}(t)\|_{0,\Omega_H}^2 dt + \|(\dot{\boldsymbol{w}}_h - \dot{\boldsymbol{u}})^{n+1/2}\|_{0,\Omega_H}^2 \\ &\le C \left( \|\ddot{\boldsymbol{w}}_h - \ddot{\boldsymbol{u}}\|_{0,0,\Omega_H}^2 + \|\dot{\boldsymbol{w}}_h - \dot{\boldsymbol{u}}\|_{0,0,\Omega_H}^2 \right) \\ &\le Ch^2 \left( \sum_{s=0}^1 \left\| \frac{\partial^s \boldsymbol{u}}{\partial t^s} \right\|_{0,2,\Omega}^2 + \sum_{s=2}^3 \left\| \frac{\partial^s \boldsymbol{u}}{\partial t^s} \right\|_{0,1,\Omega}^2 \right). \end{aligned} $$

The second term yields the same upper bound, due to Lemmas 5.2 and 5.5.

For the last two terms of (A.17) we employ Lemmas A.2 and A.3 to obtain

$$ (A.18) \qquad \begin{aligned} & \max_{0 \le n \le (J-1)} \left( \|\partial_t \boldsymbol{\xi}_h^n\|_{0,\Omega_H}^2 + \|\boldsymbol{\xi}_h^{n+1/2}\|_{1,\Omega}^2 \right) \le C \left( \|\partial_t \boldsymbol{\xi}_h^0\|_{0,\Omega_H}^2 + \|\boldsymbol{\xi}_h^{1/2}\|_{1,\Omega}^2 \right) \\ & \qquad + C \left( h^2 \left( \|\ddot{\boldsymbol{u}}\|_{0,2,\Omega}^2 + \left\| \frac{\partial^4 \boldsymbol{u}}{\partial t^4} \right\|_{0,1,\Omega}^2 \right) + \tau^4 \left\| \frac{\partial^4 \boldsymbol{u}}{\partial t^4} \right\|_{0,0,\Omega}^2 \right). \end{aligned} $$
Only the initial terms remain to be bounded. By the definitions (A.1) and (A.5), we get

$$m_H^i(\partial_t \xi_h^0, \chi_h) + \frac{\tau}{2} a(\xi_h^{1/2}, \chi_h) = m_H^i(\epsilon_h^0, \chi_h).$$

Hence, by choosing the test function $\chi_h = \partial_t \xi_h^0$, we obtain with Young's inequality

$$|\partial_t \xi_h^0|_{H,i}^2 + \frac{1}{4} (a(\xi_h^1, \xi_h^1) - a(\xi_h^0, \xi_h^0)) = m_H^i(\epsilon_h^0, \partial_t \xi_h^0) \le C |\epsilon_h^0|_{H,i}^2 + \frac{1}{2} |\partial_t \xi_h^0|_{H,i}^2.$$

Using the seminorm equivalence (4.2), this leads to

$$\|\partial_t \xi_h^0\|_{0,\Omega_H}^2 + \|\xi_h^1\|_{1,\Omega}^2 \le C (\|\epsilon_h^0\|_{0,\Omega_H}^2 + \|\xi_h^0\|_{1,\Omega}^2).$$

As we have $\xi_h^0 = \theta_h(0)$, the last term is treated as in (5.16). The first term yields by definition (A.5)

$$(A.19) \qquad \|\epsilon_h^0\|_{0,\Omega_H}^2 \le 3 \left( \|\partial_t \theta_h^0\|_{0,\Omega_H}^2 + \|\pi_h^0\|_{0,\Omega_H}^2 + \frac{\tau^2}{4} \|\rho_h^0\|_{0,\Omega_H}^2 \right).$$

Further, we have with Lemma 5.5

$$\begin{aligned} \|\partial_t \theta_h^0\|_{0,\Omega_H}^2 &= \frac{1}{\tau} \left\| \int_0^\tau \dot{\theta}_h(t) dt \right\|_{0,\Omega_H}^2 \le \int_0^\tau \|\dot{\theta}_h(t)\|_{0,\Omega_H}^2 dt \\ &\le \|Z_h \dot{\mathbf{u}} - \dot{\mathbf{u}}\|_{0,0,\Omega_H}^2 + \|\dot{\mathbf{u}} - \dot{\mathbf{w}}_h\|_{0,0,\Omega_H}^2 \\ &\le C h^2 \left( \|\dot{\mathbf{u}}\|_{0,2,\Omega}^2 + \left\| \frac{\partial^3 \mathbf{u}}{\partial t^3} \right\|_{0,1,\Omega}^2 \right). \end{aligned}$$

The second and third terms of (A.19) are estimated by means of (5.3), Lemma 5.2, and Taylor expansion:

$$\|\pi_h^0\|_{0,\Omega_H}^2 + \frac{\tau^2}{4} \|\rho_h^0\|_{0,\Omega_H}^2 \le C \tau^4 \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,0,\Omega}^2.$$

All together, the initial terms are bounded by

$$\begin{aligned} & \|\partial_t \xi_h^0\|_{0,\Omega_H}^2 + \|\xi_h^1\|_{1,\Omega}^2 \\ & \le C \left( h^2 \left( \sum_{s=0}^1 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega}^2 + \sum_{s=2}^3 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,1,\Omega}^2 \right) + \tau^4 \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,0,\Omega}^2 \right). \end{aligned}$$

Combined with (A.18), this yields the desired inequality (6.2). □

For the estimates in the $L^2(\Omega)$-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(\Omega)$-regular and that $\frac{\partial^s \mathbf{u}}{\partial t^s} \in L^2((0, T), H^2(\Omega))$ holds for all $s \in \{1, \dots, 4\}$. Let $I_h^i$ satisfy conditions P1-P3, and let $\mathbf{u}$ be the solution of (2.2) and $(\mathbf{u}_h^n)_{n=0}^J$ the solution of (6.1). Then we have the following estimate for the sequence $(\epsilon_h^n)_{N=0}^J$ defined by (A.6):

$$(A.20) \qquad \tau \sum_{n=0}^{J-1} \|\epsilon_h^n\|_{0,\Omega_H}^2 \le C \left( h^2 \sum_{s=1}^3 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + h^2 \tau^2 \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,2,\Omega} + \tau^2 \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,0,\Omega} \right)^2.$$
*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.

*Proof.* With the given initial conditions, Lemma 5.2, (A.1), and the triangle inequality, we have

$$ \|\xi_h^0\|_{0,\Omega_H} \le \|Z_h u_0 - u_0\|_{0,\Omega_H} + \|u - w_h\|_{0,0,\Omega_H} + \|\dot{u} - \dot{w}_h\|_{0,0,\Omega_H}. $$

Using (5.1) as well as Lemma 5.7, we get

$$ (A.21) \qquad \|\xi_h^0\|_{0,\Omega_H} \le C h^2 \sum_{s=0}^{3} \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega}. $$

Combining the results of Lemmas A.1 and A.4 and (A.21), we obtain

$$ \max_{0 \le n \le J} \|\xi_h^n\|_{0,\Omega_H} \le C \left( h^2 \sum_{s=0}^3 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,2,\Omega} + h^2 \tau^2 \left\| \frac{\partial^4 \mathbf{u}}{\partial t^4} \right\|_{0,2,\Omega} + \tau^2 \sum_{s=3}^4 \left\| \frac{\partial^s \mathbf{u}}{\partial t^s} \right\|_{0,0,\Omega} \right). $$

With Lemma 5.2 we finally get

$$ \max_{0 \le n \le J} \|\mathbf{u}_h^n - \mathbf{u}^n\|_{0,\Omega_H} \le \max_{0 \le n \le J} \|\xi_h^n\|_{0,\Omega_H} + \|\mathbf{u} - \mathbf{w}_h\|_{0,0,\Omega_H} + \|\dot{\mathbf{u}} - \dot{\mathbf{w}}_h\|_{0,0,\Omega_H}, $$

which concludes the proof using Lemma 5.7. □

## REFERENCES

[1] G. A. BAKER AND V. A. DOUGALIS, *The effect of quadrature errors on finite element approximations for second order hyperbolic equations*, SIAM J. Numer. Anal., 13 (1976), pp. 577–598.

[2] P. BASTIAN, K. BIRKEN, K. JOHANNSEN, S. LANG, N. NEUSS, H. RENTZ–REICHERT, AND C. WIENERS, *UG—a flexible software toolbox for solving partial differential equations*, Comput. Vis. Sci., 1 (1997), pp. 27–40.

[3] M. BORRI, C. BOTTASSO, AND L. TRAINELLI, *Integration of elastic multibody systems by invariant conserving/dissipating algorithms. II. Numerical schemes and applications*, Comput. Methods Appl. Mech. Engrg., 190 (2001), pp. 3701–3733.

[4] D. BRAESS, *Finite Elements*, Cambridge University Press, Cambridge, UK, 1997.

[5] P. CIARLET, *Mathematical Elasticity*, Vol. I, North–Holland, Amsterdam, 1998.

[6] R. DAUTRAY AND J.-L. LIONS, *Mathematical Analysis and Numerical Methods for Science and Technology. Evolution Problems*, Vol. 5, Springer, Berlin, 1992.

[7] C. ECK, J. JARUŠEK, AND M. KRBEC, *Unilateral Contact Problems. Variational Methods and Existence Theorems*, Chapman & Hall/CRC, Boca Raton, FL, 2005.

[8] L. EVANS, *Partial Differential Equations*, Grad. Stud. Math. 19, AMS, Providence, RI, 1998.

[9] C. HAGER, S. HÜEBER, AND B. WOHLMUTH, *A stable energy conserving approach for frictional contact problems based on quadrature formulas*, Internat. J. Numer. Methods Engrg., 73 (2008), pp. 205–225.

[10] P. HAURET AND P. LE TALLEC, *Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact*, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 4890–4916.

[11] H. KHENOUS, *Problèmes de contact unilatéral avec frottement de Coulomb en élastostatique et élastodynamique. Etude mathématique et résolution numérique.*, Ph.D. thesis, INSA de Toulouse, Toulouse, France, 2005.

[12] H. KHENOUS, P. LABORDE, AND Y. RENARD, *Comparison of two approaches for the discretization of elastodynamic contact problems*, C. R. Math. Acad. Sci. Paris, 342 (2006), pp. 791–796.

[13] H. KHENOUS, P. LABORDE, AND Y. RENARD, *On the discretization of contact problems in elastodynamics*, in Analysis and Simulation of Contact Problems, Lect. Notes Appl. Comput. Mech. 27, P. Wriggers and U. Nackenhorst, eds., Springer, Berlin, Heidelberg, 2006, pp. 31–38.
[14] H. KHENOUS, P. LABORDE, AND Y. RENARD, *Mass redistribution method for finite element contact problems in elastodynamics*, Eur. J. Mech. A Solids, 27 (2008), pp. 918–932.

[15] N. KIKUCHI AND J. T. ODEN, *Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods*, SIAM Stud. Appl. Math. 8, SIAM, Philadelphia, 1988.

[16] R. KORNHUBER, R. KRAUSE, O. SANDER, P. DEUFLHARD, AND S. ERTEL, *A monotone multigrid solver for two body contact problems in biomechanics*, Comput. Vis. Sci., 11 (2008), pp. 3–15.

[17] T. LAURSEN, *Computational Contact and Impact Mechanics*, Springer, Berlin, 2002.

[18] L. R. SCOTT AND S. ZHANG, *Finite element interpolation of nonsmooth functions satisfying boundary conditions*, Math. Comp., 54 (1990), pp. 483–493.

[19] V. THOMÉE, *Galerkin Finite Element Methods for Parabolic Problems*, Springer, Berlin, 1997.

[20] K. WILLNER, *Kontinuums- und Kontaktmechanik*, Springer, Berlin, Heidelberg, 2003.

[21] B. I. WOHLMUTH, *A mortar finite element method using dual spaces for the Lagrange multiplier*, SIAM J. Numer. Anal., 38 (2000), pp. 989–1012.

[22] P. WRIGGERS, *Computational Contact Mechanics*, J. Wiley & Sons, New York, 2002.