# Time domain boundary elements for elastodynamic contact

Alessandra Aimi Department of Mathematical, Physical and Computer Sciences,
University of Parma, Parco Area delle Scienze, 53/A, 43124, Parma, Italy,
email<EMAIL_ADDRESS><EMAIL_ADDRESS>
Members of the INDAM-GNCS Research Group, Italy. Giulia Di Credico∗
Engineering Mathematics, University of Innsbruck, Innsbruck, Austria, email:
<EMAIL_ADDRESS>Heiko Gimperlein†

###### Abstract

This article proposes a boundary element method for the dynamic contact
between a linearly elastic body and a rigid obstacle. The Signorini contact
problem is formulated as a variational inequality for the Poincaré-Steklov
operator for the elastodynamic equations on the boundary, which is solved in a
mixed formulation using boundary elements in the time domain. We obtain an a
priori estimate for the resulting Galerkin approximations. Numerical
experiments confirm the stability and convergence of the proposed method for
the contact problem in flat and curved two-dimensional geometries, as well as
for moving obstacles.

Key words: boundary element methods; space-time methods; Signorini contact
problem; variational inequality; elastodynamics.

##  1 Introduction

Contact problems between elastic bodies play a crucial role in applications
ranging from fracture dynamics to rolling car tires [44]. As the contact takes
place at the interface of two materials, boundary elements as well as the
coupling of finite and boundary elements lead to efficient and much-studied
tools for numerical simulations for time-independent situations as well as for
strongly dissipative materials [30, 43]. Mathematically, the analysis of the
resulting variational inequalities has been well understood for elliptic and
parabolic problems.
In spite of the practical relevance of dynamic contact and the challenges of
simulations based on finite elements, rigorous boundary element methods are
only starting to be developed. We refer to [10, 11, 12, 16, 27, 28, 32, 35,
37, 39] as examples from the extensive mathematical finite element literature.
Indicating the challenges, without dissipation only stability or energy
conservation, but not the convergence of finite element methods, is known.
The computational difficulties relate to the mathematical challenges of the
time-dependent contact problem [17]. Even the existence of weak solutions is
known only for viscoelastic materials or modified contact conditions, such as
in [13]. Recent progress towards the numerical analysis problems with
unilateral constraints and without dissipative terms was made in [22], for the
simplified problem of a scalar wave equation which arises in the physical
limit when transversal stresses can be neglected [38]. Unlike for elastic
contact, in this scalar problem refined information about the Poincaré-Steklov
operator is available for the analysis [14, 36].
Building on [22] and [7], in this work we propose a time domain boundary
element method for the elastodynamic Signorini contact problem corresponding
to contact with an impenetrable obstacle.

To be specific, we consider the dynamics of a linear elastic body with
homogeneous mass density $\varrho$ in a bounded domain
$\Omega\subset\mathbb{R}^{d}$, $d=2,3$. Its dynamics is described by the
Navier-Lamé equations

$\nabla\cdot\sigma(\textbf{u})-\varrho\ddot{\textbf{u}}=\textbf{0}$ (1)

for times $t\in[0,T]$. Here u is the unknown displacement vector and the upper
dots denote time derivatives. The Cauchy stress tensor $\sigma(\textbf{u})$ is
defined as follows (see e.g. [18])

$\sigma(\textbf{u})=\varrho(c_{\mathtt{P}}^{2}-2c_{\mathtt{S}}^{2})\,(\nabla\cdot\textbf{u})I+\varrho
c_{\mathtt{S}}^{2}\,(\nabla\textbf{u}+\nabla\textbf{u}^{\top})\,,$

with pressure and shear velocities $c_{\mathtt{P}}$, respectively
$c_{\mathtt{S}}$. From $\sigma(\textbf{u})$ we obtain the elastic traction

$\textbf{p}=\sigma\left(\textbf{u}\right)_{|_{\Gamma}}\bf n\,$ (2)

on the boundary $\Gamma$ of $\Omega$. Here n denotes the outward-pointing unit
normal vector to $\Gamma$.

Indicating with the subscript $\perp$ the normal component of a vector, we
consider the non-penetration of an obstacle using the following (non-linear)
Signorini boundary conditions for ${u}_{\perp},p_{\perp}$ on a given contact
boundary $\Gamma_{C}\subset\Gamma$:

$\begin{cases}{u}_{\perp}\geq g\ ,\,{p}_{\perp}\geq{f}_{\perp}\ ,\\\
{u}_{\perp}>g\ \Longrightarrow{p}_{\perp}={f}_{\perp}\ .\end{cases}$ (3)

The mechanical interpretation of the gap function $g$ and the prescribed force
${\bf f}$ are discussed in Section 2.1 below. As described there, the
Signorini conditions are complemented by boundary conditions prescribing the
displacement, respectively the traction, on the remainder of $\Gamma$.
In this article we reduce the contact problem to a variational inequality on
the boundary $\Gamma$, using the elastodynamic Poincaré-Steklov operator.
Because the contact area and forces are often relevant in applications, we
reformulate this variational inequality in an equivalent mixed formulation,
which we discretize with a time domain Galerkin boundary element method. The
mixed formulation is solved with an Uzawa algorithm to obtain the displacement
and the contact forces.
The detailed algebraic formulation and implementation of an energetic space-
time boundary element method are presented in this article. As a main
theoretical result we obtain an a priori estimate for the numerical error in
Theorem 7, as well as the convergence of the proposed Uzawa algorithm. Given
the analytical challenges described above, the theoretical results are subject
to an assumption known for simplified situations [22].
Numerical results confirm the stability and convergence of the proposed method
in two dimensions. Both flat and curved contact boundaries are considered, as
is the contact with time-dependent, moving obstacles.
Let us finally note that our approach relies on the recent advances in time
domain boundary elements and coupled finite element / boundary element
procedures for interface problems, where both space-time Galerkin and
convolution quadrature methods have been of interest [1, 3, 6, 19, 21, 33, 34,
40, 41, 42]. More specifically, we refer to [2, 4, 8, 9] for works on the
theoretical and numerical analysis of the elastodynamic boundary integral
operators.
This article is structured as follows: in Section 2 the differential model
problem is set up together with its reformulation in terms of space-time
boundary integral operators; in Section 3 the adopted discretization is
introduced, while in Section 4 the theoretical analysis is conducted.
Algorithmic details of the overall implementation are described in Section 5
and extensive numerical simulations are presented and discussed in Section 6.
Conclusions are briefly drawn in Section 7.

##  2 Dynamic contact and boundary integral formulations

###  2.1 Differential problem formulation

Figure 1: Schematic depiction of contact problem.

In this subsection we introduce the equations which govern the contact
problems considered in this article. As illustrated in Figure 1 for $d=2$, we
consider the time dependent deformation of a linearly elastic body, described
by the Navier-Lamé equations (1) for the displacement u in a bounded domain
$\Omega\subset\mathbb{R}^{d}$, $d=2,3$. Starting from the reference
configuration, $\textbf{u}=\textbf{0}$ for times $t\leq 0$, the dynamics of
the body is due to surface forces applied on a subset $\Gamma_{N}$ of its
boundary $\Gamma$, where the traction $\textbf{p}=\mathbf{f}$ is prescribed.

The unilateral Signorini contact problem here considered describes the
physical impossibility of the body to penetrate an adjacent rigid and
frictionless surface. In Figure 1 the impenetrable obstacle is given by the
positive $x$\- and $y$-axes. Nonpenetration of the obstacle leads to the
contact condition (3) on the part $\Gamma_{C}\subset\Gamma$ of the boundary
where contact may occur. More precisely, from (3) contact takes place when the
normal displacement satisfies $u_{\perp}=g$, where $g$ describes the gap
between the reference configuration and the obstacle. Contact is avoided when
$u_{\perp}>g$, and then only the applied surface forces
${p}_{\perp}={f}_{\perp}$ act on $\Gamma_{C}$.

It will be convenient to denote by $\bar{\Gamma}_{\Sigma}$ the union
$\bar{\Gamma}_{N}\cup\bar{\Gamma}_{C}$ of the traction and contact boundaries.
To complete the description of the problem, the body is fixed on the remaining
subset $\Gamma_{D}=\Gamma\setminus\bar{\Gamma}_{\Sigma}$ of the boundary. This
results in the boundary condition $\textbf{u}=\textbf{0}$ for the displacement
on $\Gamma_{D}$.
Indicating with the subscript $\parallel$ the component of a vector tangential
to $\Gamma$, we can summarize the governing equations as follows:

$\displaystyle\nabla\cdot\sigma(\textbf{u})-\varrho\ddot{\textbf{u}}$
$\displaystyle=\textbf{0}$ $\displaystyle\text{in }[0,T]\times\Omega$ (4a) u
$\displaystyle=\textbf{0}$ $\displaystyle\text{on }[0,T]\times\Gamma_{D}$ (4b)
p $\displaystyle=\mathbf{f}$ $\displaystyle\text{on }[0,T]\times\Gamma_{N}$
(4c) $\displaystyle\textbf{p}_{\parallel}$
$\displaystyle=\textbf{f}_{\parallel}$ $\displaystyle\text{on
}[0,T]\times\Gamma_{C}$ (4d)

together with the contact conditions (3) on $[0,T]\times\Gamma_{C}$ and
homogeneous condition $\mathbf{u}\equiv\mathbf{0}$ in $\Omega$ for $t\leq 0$.

###  2.2 Boundary integral formulations: variational inequality and mixed
formulations

Boundary integral formulations are well-known to lead to efficient numerical
methods for time dependent contact problems (see [22] in the context of
acoustic wave equation), as the dynamics in $\Omega$ can be fully described by
an integral equation on $[0,T]\times\Gamma$.
For the remainder of this article we assume that $\Omega$ is a bounded
polygonal or polyhedral Lipschitz domain. To derive the corresponding
formulations for the elastodynamic equations (1), we introduce the time-
dependent single layer and double layer potential operators $V$ and $K$, that
act onto the generic fields $\boldsymbol{\phi}$ and $\boldsymbol{\psi}$, for
$(t,{\bf x})\in[0,T]\times\Omega$, as follows:

$[V\boldsymbol{\phi}](\mathbb{x},t)=\int_{0}^{t}\int_{\Gamma}G(t,\tau;\mathbb{x},\mathbb{y})\>{\boldsymbol{\phi}}(\tau,\mathbb{y})d\Gamma_{\mathbb{y}}d\tau,\>\>[K\boldsymbol{\psi}](\mathbb{x},t)=\int_{0}^{t}\int_{\Gamma}\left[\sigma_{\mathbb{y}}\left(G\right)^{\top}(t,\tau;\mathbb{x},\mathbb{y})\mathbb{n}_{\mathbb{y}}\right]\boldsymbol{\psi}(\mathbb{y},\tau)d\Gamma_{\mathbb{y}}d\tau,$

where $G$ is the fundamental solution to (1). The subscript $\mathbb{y}$
applied to the stress tensor $\sigma$ denotes the variable for the application
of the spatial derivative, while for the vector $\mathbb{n}$ it declares the
point of $\Gamma$ where we are considering the normal direction. For $d=2$ the
fundamental solution reads as

$\displaystyle G_{ij}(\textbf{x},\textbf{y};t,\tau):=$
$\displaystyle\dfrac{H[c_{\mathtt{P}}(t-\tau)-r]}{2\pi\varrho
c_{\mathtt{P}}}\left\\{\dfrac{r_{i}r_{j}}{r^{4}}\dfrac{2c^{2}_{\mathtt{P}}(t-\tau)^{2}-r^{2}}{\sqrt{c_{\mathtt{P}}^{2}(t-\tau)^{2}-r^{2}}}-\dfrac{\delta_{ij}}{r^{2}}\sqrt{c^{2}_{\mathtt{P}}(t-\tau)^{2}-r^{2}}\right\\}$
(5) $\displaystyle-$
$\displaystyle\dfrac{H[c_{\mathtt{S}}(t-\tau)-r]}{2\pi\varrho
c_{\mathtt{S}}}\left\\{\dfrac{r_{i}r_{j}}{r^{4}}\dfrac{2c^{2}_{\mathtt{S}}(t-\tau)^{2}-r^{2}}{\sqrt{c_{\mathtt{S}}^{2}(t-\tau)^{2}-r^{2}}}-\dfrac{\delta_{ij}}{r^{2}}\dfrac{c^{2}_{\mathtt{S}}(t-\tau)^{2}}{\sqrt{c^{2}_{\mathtt{S}}(t-\tau)^{2}-r^{2}}}\right\\},\quad
i,j=1,2,$

being $H[\cdot]$ the Heaviside step function, $r_{j}$ and $r$ the $j$-th
cartesian component and the euclidean norm of the vector
$\textbf{r}=\textbf{x}-\textbf{y}$, respectively, and $\delta_{ij}$ the
Kronecker delta. For $d=3$ the fundamental solution reads instead as

$\displaystyle G_{ij}(\textbf{x},\textbf{y};t,\tau):=$
$\displaystyle\frac{t-\tau}{4\pi\varrho
r^{2}}\left(\frac{r_{i}r_{j}}{r^{3}}-\frac{\delta_{ij}}{r}\right)(H[c_{\mathtt{P}}(t-\tau)-r]-H[c_{\mathtt{S}}(t-\tau)-r])$
$\displaystyle+\frac{r_{i}r_{j}}{4\pi\varrho
r^{3}}\left(c_{\mathtt{P}}^{-2}\delta(c_{\mathtt{P}}(t-\tau)-r)-c_{\mathtt{S}}^{-2}\delta(c_{\mathtt{S}}(t-\tau)-r)\right)$
$\displaystyle+\frac{\delta_{ij}}{4\pi\varrho
rc_{\mathtt{S}}^{2}}\delta(c_{\mathtt{S}}(t-\tau)-r),\qquad i,j=1,2,3,$

where $\delta(\cdot)$ is the Dirac delta distribution.

The unknown displacement u can be expressed by the representation formula

$\textbf{u}=V\textbf{p}-K\textbf{u},\quad\textrm{in}\;[0,T]\times\Omega.$ (6)

Letting $\textbf{x}\in\Omega\rightarrow\textbf{x}\in\Gamma$ in (6), we deduce
the boundary integral equation in the standard notation

$\frac{1}{2}\textbf{u}=\mathcal{V}\textbf{p}-\mathcal{K}\textbf{u},\quad\textrm{in}\;[0,T]\times\Gamma,$
(7)

taking into account the free term $\frac{1}{2}\textbf{u}$ generated by the
operator $K$ pursuant the limiting process. Moreover, we need also to
introduce the adjoint double layer operator $\mathcal{K}^{\star}$ and the
hypersingular integral operator $\mathcal{W}$ for
$(\textbf{x},t)\in\Gamma\times[0,T]$:

$\displaystyle[\mathcal{K}^{\star}\boldsymbol{\phi}](\mathbb{x},t)$
$\displaystyle=\int_{0}^{t}\int_{\Gamma}\left[\sigma_{\mathbb{x}}\left(G\right)(t,\tau;\mathbb{x},\mathbb{y})\mathbb{n}_{\mathbb{x}}\right]{\boldsymbol{\phi}}(\tau,\mathbb{y})d\Gamma_{\mathbb{y}}\
d\tau,$ $\displaystyle[\mathcal{W}\boldsymbol{\psi}](\mathbb{x},t)$
$\displaystyle=\int_{0}^{t}\int_{\Gamma}\left[\sigma_{\textbf{x}}\left(\sigma_{\mathbb{y}}\left(G\right)^{\top}(t,\tau;\mathbb{x},\mathbb{y})\mathbb{n}_{\mathbb{y}}\right)\textbf{n}_{\textbf{x}}\right]\boldsymbol{\psi}(\mathbb{y},\tau)d\Gamma_{\mathbb{y}}\
d\tau\,,$

involved in the classical boundary integral equation

$\frac{1}{2}\textbf{p}=\mathcal{K}^{*}\textbf{p}-\mathcal{W}\textbf{u},\quad\textrm{in}\;[0,T]\times\Gamma.$
(8)

The above recalled integral operators will be used, as detailed in Section
5.1, in the boundary integral problem we are going to solve, which is based on
the Poincaré-Steklov operator $\mathcal{S}$, defined by

$\mathcal{S}\left({\bf
u}_{|_{\Gamma}}\right):=\sigma\left(\textbf{u}\right)_{|_{\Gamma}}{\bf
n}=\mathbf{p},$ (9)

where ${\bf u}$ is a solution to the elastodynamic equations (1) in
$[0,T]\times\Omega$, given a Dirichlet datum ${\bf u}_{|_{\Gamma}}$, and the
last equality is due to the definition of traction in (2).
Basic properties of $\mathcal{S}$ are summarized in Theorem 5 below. Moreover,
let us remark that this operator can be expressed by two equivalent forms,
matching properly the previously introduced integral operators, namely: the
non-symmetric formulation

$\mathcal{S}=\mathcal{V}^{-1}\left(\mathcal{K}+\frac{1}{2}\right)$ (10)

and the symmetric formulation

$\mathcal{S}=\left(\mathcal{K}^{*}+\frac{1}{2}\right)\mathcal{V}^{-1}\left(\mathcal{K}+\frac{1}{2}\right)-\mathcal{W}\,.$
(11)

In the following, to simplify the notation the subscript $|_{\Gamma}$ in the
argument of the operator $\cal S$ will be omitted, whenever clear from the
context.

Now, let us denote $L^{2}$ space-time scalar products useful in the sequel by

$\displaystyle\langle\textbf{u},\textbf{v}\rangle_{0,\Gamma,(0,T]}:=\int_{0}^{T}\int_{\Gamma}\textbf{u}(\textbf{x},t)\cdot\textbf{v}(\textbf{x},t)\>d\Gamma_{\textbf{x}}\>dt,$
(12)
$\displaystyle\langle\textbf{u},\textbf{v}\rangle_{\sigma,\Gamma,\mathbb{R}^{+}}:=\int_{0}^{\infty}e^{-2\sigma
t}\int_{\Gamma}\textbf{u}(\textbf{x},t)\cdot\textbf{v}(\textbf{x},t)\>d\Gamma_{\textbf{x}}\>dt,$
(13)

where $\sigma>0$, changing the subscript $\Gamma$ whenever dealing with
boundary subsets and including in the above definitions the simpler case of
scalar functions. For precise statements and the error analysis we also
require space-time Sobolev spaces
$H^{r}(\mathcal{I},\tilde{H}^{s}({\Gamma^{\prime}}))$ on subsets
${\Gamma^{\prime}}\subset\Gamma$ and for the time intervals
$\mathcal{I}=[0,T],\mathbb{R}^{+}$, which are introduced in Appendix A.

Given $g\in H^{1/2}([0,T],H^{1/2}(\Gamma_{C}))$ and $\textbf{f}\in
H^{1/2}([0,T],\tilde{H}^{-1/2}(\Gamma_{\Sigma}))^{d}$, the precise functional
analytic formulation of the contact problem (3)-(4) as a variational
inequality in terms of $\mathcal{S}$ reads:

find $\textbf{u}\in\mathcal{C}:=\left\\{\textbf{v}\in
H^{1/2}([0,T],\tilde{H}^{1/2}({\color[rgb]{0,0,0}\Gamma_{\Sigma}}))^{d}:v_{\perp}\geq
g\text{\leavevmode\nobreak\ a.e.\leavevmode\nobreak\ on\leavevmode\nobreak\
}[0,T]\times\Gamma_{C}\right\\}$ such that

$\displaystyle\langle{\mathcal{S}{\bf u}},{\bf v}-{\bf
u}\rangle_{0,\Gamma_{\Sigma},(0,T]}\geq\langle{\bf f},{\bf v}-{\bf
u}\rangle_{0,\Gamma_{\Sigma},(0,T]}\qquad\quad\forall{\bf v}\in\mathcal{C}.$
(14)

For the following result, we recall that, as typical for direct boundary
element methods, the displacement $\mathbf{u}$ in the domain
$[0,T]\times\Omega$ is recovered from its boundary trace
$\mathbf{u}|_{\Gamma}$ on $[0,T]\times\Gamma$ using the representation formula
(6).

######  Proposition 1.

The variational inequality (14) for the displacement $\mathbf{u}|_{\Gamma}$ on
$[0,T]\times\Gamma$ is equivalent to the contact problem (3)-(4) for the
solution $\mathbf{u}$ in $[0,T]\times\Omega$.

###### Proof.

First, we show that the boundary trace $\mathbf{u}|_{\Gamma}$ of the solution
$\mathbf{u}$ to (3)-(4) satisfies (14). Note that the Signorini boundary
condition (3) on $\Gamma_{C}$ can be expressed in terms of $\mathcal{S}$ using
(9):

$\begin{cases}{u}_{\perp}\geq g\
,\,\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}\geq{f}_{\perp}\ ,\\\
{u}_{\perp}>g\
\Longrightarrow\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}={f}_{\perp}\
.\end{cases}$ (15)

Equivalently, it is written as

${u}_{\perp}\geq g\ ,\
\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}\geq{f}_{\perp}\ ,\
({u}_{\perp}-g)(\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}-f_{\perp})=0\ .$
(16)

Due to (4d), the parallel component on $\Gamma_{C}$ satisfies

$\mathcal{S}(\mathbf{u}|_{\Gamma})_{\parallel}=f_{\parallel}.$

Further, from (4c), the Neumann boundary conditions on $\Gamma_{N}$ are
equivalent to

$\mathcal{S}(\mathbf{u}|_{\Gamma})=\mathbf{f}.$

We therefore find for any $\mathbf{v}\in\mathcal{C}$

$\displaystyle\langle\mathcal{S}(\mathbf{u}|_{\Gamma})-\mathbf{f},\mathbf{v}-\mathbf{u}|_{\Gamma}\rangle_{0,\Gamma_{\Sigma},(0,T]}$
$\displaystyle=\langle\mathcal{S}(\mathbf{u}|_{\Gamma})_{\parallel}-f_{\parallel},(\mathbf{v}-\mathbf{u}|_{\Gamma})_{\parallel}\rangle_{0,\Gamma_{\Sigma},(0,T]}+\langle\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}-f_{\perp},(\mathbf{v}-\mathbf{u}|_{\Gamma})_{\perp}\rangle_{0,\Gamma_{\Sigma},(0,T]}$
$\displaystyle=0+\langle\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}-f_{\perp},v_{\perp}-g\rangle_{0,\Gamma_{\Sigma},(0,T]}-\langle\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}-{f}_{\perp},u_{\perp}-g\rangle_{0,\Gamma_{\Sigma},(0,T]}$
$\displaystyle=\langle\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}-f_{\perp},v_{\perp}-g\rangle_{0,\Gamma_{\Sigma},(0,T]}+0$
$\displaystyle\geq 0.$

Equivalently,

$\langle\mathcal{S}(\mathbf{u}|_{\Gamma}),\mathbf{v}-\mathbf{u}\rangle_{0,\Gamma_{\Sigma},(0,T]}\geq\langle\mathbf{f},\mathbf{v}-\mathbf{u}\rangle_{0,\Gamma_{\Sigma},(0,T]},$

and (14) follows.

To show that, conversely, a sufficiently smooth solution ${\bf
u}_{|_{\Gamma}}$ of (14) leads to a solution ${\bf u}$ of the contact problem
(3)-(4), we first recall that the representation formula (6) recovers the
solution $\mathbf{u}$ to (4a) in the domain $[0,T]\times\Omega$ from its
restriction to $[0,T]\times\Gamma$. It remains to show that $\mathbf{u}$
satisfies the boundary conditions on $\Gamma$. Because
$\textbf{u}\in\mathcal{C}\subset
H^{1/2}([0,T],\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}$, the homogeneous
Dirichlet condition is satisfied on $\Gamma_{D}$.

We now choose $\mathbf{v}=\mathbf{u}|_{\Gamma}+\mathbf{w}\in\mathcal{C}$, for
$\mathbf{w}\in H^{1/2}([0,T],\tilde{H}^{1/2}({\Gamma_{\Sigma}}))^{d}$ with
$w_{\perp}=0$ a.e. on $[0,T]\times\Gamma_{C}$. Then (14) implies

$\langle\mathcal{S}(\mathbf{u}|_{\Gamma}),\mathbf{w}\rangle_{0,\Gamma_{\Sigma},(0,T]}\geq\langle\mathbf{f},\mathbf{w}\rangle_{0,\Gamma_{\Sigma},(0,T]}.$

Similarly, $\mathbf{v}=\mathbf{u}|_{\Gamma}-\mathbf{w}\in\mathcal{C}$ and
therefore also

$\langle\mathcal{S}(\mathbf{u}|_{\Gamma}),\mathbf{w}\rangle_{0,\Gamma_{\Sigma},(0,T]}\leq\langle\mathbf{f},\mathbf{w}\rangle_{0,\Gamma_{\Sigma},(0,T]}.$

We conclude

$\langle\mathcal{S}(\mathbf{u}|_{\Gamma}),\mathbf{w}\rangle_{0,\Gamma_{\Sigma},(0,T]}=\langle\mathbf{f},\mathbf{w}\rangle_{0,\Gamma_{\Sigma},(0,T]}.$

As this holds for all $\mathbf{w}\in
H^{1/2}([0,T],\tilde{H}^{1/2}({\Gamma_{\Sigma}}))^{d}$ with $w_{\perp}=0$ a.e.
on $[0,T]\times\Gamma_{C}$, we conclude
$\mathbf{p}=\mathcal{S}(\mathbf{u}|_{\Gamma})=\mathbf{f}$ on $\Gamma_{N}$ and
$p_{\parallel}=\mathcal{S}(\mathbf{u}|_{\Gamma})_{\parallel}=f_{\parallel}$ on
$\Gamma_{C}$.

It remains to show the Signorini conditions on $\Gamma_{C}$. For this we note
that $u_{\perp}\geq g$ on $\Gamma_{C}$ because
$\mathbf{u}|_{\Gamma}\in\mathcal{C}$. We now choose
$\mathbf{v}=\mathbf{u}|_{\Gamma}+\mathbf{w}\in\mathcal{C}$ with $w_{\perp}\geq
0$ on $\Gamma_{C}$. As we have already concluded
$\mathcal{S}(\mathbf{u}|_{\Gamma})_{\parallel}=f_{\parallel}$, we obtain

$0\geq\langle-\mathcal{S}(\mathbf{u}|_{\Gamma})+\mathbf{f},\mathbf{w}\rangle_{0,\Gamma_{\Sigma},(0,T]}=\langle-\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}+f_{\perp},w_{\perp}\rangle_{0,\Gamma_{\Sigma},(0,T]}.$

Therefore the boundary condition
$\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}\geq f_{\perp}$ is satisfied almost
everywhere on $\Gamma_{C}$.

Finally, for every $\mathbf{v}\in\mathcal{C}$ with
$\mathbf{v}|_{\Gamma_{C}}=g\mathbf{n}$, we similarly note

$0\geq\langle-\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}+f_{\perp},g-u_{\perp}\rangle_{0,\Gamma_{\Sigma},(0,T]}.$
(17)

In the previous step we saw that
$0\geq-\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}+f_{\perp}$ almost everywhere,
while $0\geq g-u_{\perp}$ for $\mathbf{u}|_{\Gamma}\in\mathcal{C}$. Hence also
$0\leq(-\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}+f_{\perp})(g-u_{\perp})$
almost everywhere. Together with (17) we conclude that

$0=(-\mathcal{S}(\mathbf{u}|_{\Gamma})_{\perp}+f_{\perp})(g-u_{\perp})$

almost everywhere on $\Gamma_{C}$, i.e. the equality in (16).

Therefore, $\mathbf{u}$ satisfies the boundary conditions of the contact
problem (3)-(4). ∎

In spite of the interest in the dynamic contact problem (3)-(4), from a
rigorous mathematical perspective the existence of solutions is only known for
certain materials with dissipation or for not perfectly rigid, but dissipative
obstacles, see e.g. [13]. Without dissipation the existence of a solution has
been proven for simplified problems involving the scalar wave equation in
special geometries [14, 36]. In these cases refined information about the
coercivity of the Poincaré-Steklov operator is available, as discussed in
[22].

For the numerical approximation we will consider a mixed formulation. Let

$\displaystyle M^{+}:=\left\\{{\mu}\in
H^{1/2}([0,T],\tilde{H}^{-1/2}(\Gamma_{C})):\left\langle{\mu},{w}\right\rangle_{0,\Gamma_{C},(0,T]}\geq
0,\ \forall\ 0\leq{w}\in
H^{-\frac{1}{2}}([0,T],\tilde{H}^{1/2}(\Gamma_{\Sigma}))\right\\}$ (18)

be the set of admissible Lagrange multipliers, in which the representative

${\lambda}=(\mathcal{S}{\bf u}-{\bf f})_{\perp}\in M^{+}$ (19)

is sought. Note that by definition ${\lambda}\geq 0$ and $\lambda=0$ outside
$\Gamma_{C}$. Then, the mixed formulation reads:

find $(\textbf{u},{\lambda})\in
H^{1/2}([0,T],\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}\times M^{+}$ such that

$\displaystyle\left\langle\mathcal{S}\textbf{u},\textbf{v}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}-\left\langle\lambda,v_{\perp}\right\rangle_{0,\Gamma_{C},(0,T]}$
$\displaystyle=\left\langle\textbf{f},\textbf{v}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}$
$\displaystyle\forall\textbf{v}\in
H^{1/2}([0,T],\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}$ (20a)
$\displaystyle\left\langle{\mu}-{\lambda},u_{\perp}\right\rangle_{0,\Gamma_{C},(0,T]}$
$\displaystyle\geq\left\langle
g,\mu-\lambda\right\rangle_{0,\Gamma_{C},(0,T]}$ $\displaystyle\forall{\mu}\in
M^{+}.$ (20b)

######  Theorem 2.

The mixed formulation (20) is equivalent to the variational inequality (14).

###### Proof.

The proof is identical to the proof of Theorem 14 in [22] for the wave
equation. ∎

##  3 Discretization

To solve the mixed formulation (20), in a discretized form, we consider a
uniform decomposition of the time interval $[0,T]$ with time step $\Delta
t=\frac{T}{N_{\Delta t}}$, $N_{\Delta t}\in\mathbb{N}^{+}$, generated by the
time instants $t_{\ell}=\ell\Delta t$, $\ell=0,\ldots,N_{\Delta t}$. We define
the corresponding spaces

$\begin{array}[]{ll}V^{-1}_{\Delta t}=&\left\\{v_{\Delta t}\in
L^{2}([0,T])\>:\>{v_{\Delta
t}}_{|_{[t_{\ell},t_{\ell+1}]}}\in\mathcal{P}_{0},\>\forall\ell=0,...,N_{\Delta
t}-1\right\\},\\\\[4.0pt] V^{0}_{\Delta t}=&\left\\{r_{\Delta t}\in
C^{0}([0,T])\>:\>{r_{\Delta
t}}_{|_{[t_{\ell},t_{\ell+1}]}}\in\mathcal{P}_{1},\>\forall\ell=0,...,N_{\Delta
t}-1,\>v(0)=0\right\\},\end{array}$ (21)

where $\mathcal{P}_{s}$, $s\geq 0$, is the space of the algebraic polynomials
of degree s. For the space discretization with $d=2$, we introduce a boundary
mesh constituted by a set of straight line segments
$\mathcal{T}=\left\\{e_{1},...,e_{M}\right\\}$ such that
$h_{i}:=length(e_{i})\leqslant h$, $e_{i}\cap e_{j}=\emptyset$ if $i\neq j$
and $\cup_{i=1}^{M}\overline{e}_{i}=\overline{\Gamma}$ if $\Gamma$ is
polygonal, or a suitably fine approximation of $\Gamma$ otherwise. For $d=3$,
we assume that $\Gamma$ is triangulated by
$\mathcal{T}=\\{e_{1},\cdots,e_{M}\\}$, with $h_{i}:=diam(e_{i})\leqslant h$,
$e_{i}\cap e_{j}=\emptyset$ if $i\neq j$ and, if
$\overline{e_{i}}\cap\overline{e_{j}}\neq\emptyset$, the intersection is
either an edge or a vertex of both triangles. We denote by $h$ the maximum of
the $h_{i}$. On $\mathcal{T}$ we consider the spaces of piecewise polynomial
functions

$X^{-1}_{h,\Gamma}=\left\\{w_{h}\in
L^{2}(\Gamma)\>:\>w_{h}|_{e_{i}}\in\mathcal{P}_{s},\>e_{i}\in\mathcal{T}\right\\}\subset
H^{-1/2}(\Gamma),\quad
X^{-1}_{h,\Gamma^{\prime}}=X^{-1}_{h,\Gamma}\cap\widetilde{H}^{-1/2}(\Gamma^{\prime}),$
(22) $X^{0}_{h,\Gamma}=\left\\{w_{h}\in
C^{0}(\Gamma)\>:\>w_{h}|_{e_{i}}\in\mathcal{P}_{s},\>e_{i}\in\mathcal{T}\right\\}\subset
H^{1/2}(\Gamma),\quad
X^{0}_{h,\Gamma^{\prime}}=X^{0}_{h,\Gamma}\cap\widetilde{H}^{1/2}(\Gamma^{\prime}),$
(23)

where $\Gamma^{\prime}\subset\Gamma$.
The full discretization of (20) involves the subspace

$\displaystyle M^{+}_{H,\Delta T}:=\left\\{{\mu}_{H,\Delta T}\in
X^{-1}_{H,\Gamma_{C}}\otimes V^{-1}_{\Delta T}:\mu_{H,\Delta T}\geq 0\text{ on
}\Gamma_{C}\times[0,T]\right\\}\subset M^{+}$ (24)

and a discretized version ${\cal S}_{h,\Delta t}$ of the Poincaré-Steklov
operator. Two specific discretizations ${\cal S}_{h,\Delta t}$ will be
discussed in detail in Section 5.1, corresponding to the symmetric and the
non-symmetric formulations (10), (11) of the operator ${\cal S}$. The full
discretization of (20) then reads:

find $(\textbf{u}_{h,\Delta t},{\lambda}_{H,\Delta
T})\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t})^{d}\times
M^{+}_{H,\Delta T}$ such that

$\displaystyle\left\langle\mathcal{S}_{h,\Delta t}\textbf{u}_{h,\Delta
t},\textbf{v}_{h,\Delta
t}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}-\left\langle{\lambda}_{H,\Delta
T},{v}_{\perp,h,\Delta
t}\right\rangle_{0,\Gamma_{C},(0,T]}=\left\langle\textbf{f},\textbf{v}_{h,\Delta
t}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}\quad\forall\textbf{v}_{h,\Delta
t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t})^{d}$ (25a)
$\displaystyle\left\langle{\mu}_{H,\Delta T}-{\lambda}_{H,\Delta
T},{u}_{\perp,h,\Delta t}\right\rangle_{0,\Gamma_{C},(0,T]}\geq\left\langle
g,\mu_{H,\Delta T}-\lambda_{H,\Delta
T}\right\rangle_{0,\Gamma_{C},(0,T]}\quad\quad\quad\quad\quad\quad\>\forall{\mu}_{H,\Delta
T}\in M^{+}_{H,\Delta T}.$ (25b)

######  Theorem 3.

The discretized mixed formulation (25) admits a unique solution provided that
${\cal S}_{h,\Delta t}$ is positive definite.

###### Proof.

As ${\cal S}_{h,\Delta t}$ is positive definite, the existence of a unique
solution follows from standard results for saddle point problems. ∎

Theoretically, it is not generally known that ${\cal S}_{h,\Delta t}$ is
positive definite, a main difficulty in the analysis of contact problems as
discussed in [22]. Empirically, we obtain a well-conditioned, discretized
mixed formulation with, in particular, a unique solution.

A standard solver for the discrete formulation (25) is given by the Uzawa
algorithm. This is described in the following and it involves the
$L^{2}$-projection $\text{Pr}_{C}:L^{2}((0,T];L^{2}(\Gamma))\to
M^{+}_{H,\Delta T}$:

Algorithm 1 (Uzawa algorithm)

Fix $\rho>0$.

$k=0$, $\lambda_{H,\Delta T}^{(0)}=0$

while stopping criterion not satisfied do

solve  equation (25a) for ${\bf u}_{h,\Delta t}^{(k)}$

compute  ${\lambda}^{(k+1)}_{H,\Delta T}=\text{Pr}_{C}({\lambda}_{H,\Delta
T}^{(k)}-\rho({u}_{\perp,h,\Delta t}^{(k)}-g))$

$k\leftarrow k+1$

end while

######  Lemma 4.

Assume that ${\cal S}_{h,\Delta t}$ is positive definite. The space-time Uzawa
algorithm converges, provided that $0<\rho<2C$, with $C$ the smallest
eigenvalue of ${\cal S}_{h,\Delta t}$.

###### Proof.

Note that $\lambda_{H,\Delta T}=\text{Pr}_{C}\ {\lambda}_{H,\Delta
T}=\text{Pr}_{C}\ ({\lambda}_{H,\Delta T}-\rho({u}_{\perp,h,\Delta t}-g))$,
because $\lambda_{H,\Delta T}\in M^{+}_{H,\Delta T}$ and
$\rho({u}_{\perp,h,\Delta t}-g)\geq 0$. From the algorithm and the contraction
property of the projection $\textrm{Pr}_{C}$, one has

$\displaystyle\|{\lambda}_{H,\Delta T}^{(k+1)}-{\lambda}_{H,\Delta
T}\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}^{2}$ $\displaystyle=\|\text{Pr}_{C}\
({\lambda}_{H,\Delta T}^{(k)}-\rho({u}_{\perp,h,\Delta
t}^{(k)}-g))-\text{Pr}_{C}\ ({\lambda}_{H,\Delta T}-\rho({u}_{\perp,h,\Delta
t}-g))\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}^{2}$
$\displaystyle\leq\|{\lambda}_{H,\Delta T}^{(k)}-{\lambda}_{H,\Delta
T}-\rho({u}_{\perp,h,\Delta t}^{(k)}-{u}_{\perp,h,\Delta
t})\|_{L^{2}([0,T],L^{2}(\Gamma_{\Sigma}))}^{2}$
$\displaystyle=\|{\lambda}_{H,\Delta T}^{(k)}-{\lambda}_{H,\Delta
T}\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}^{2}+\rho^{2}\|{u}_{\perp,h,\Delta
t}^{(k)}-{u}_{\perp,h,\Delta t}\|^{2}_{L^{2}([0,T],L^{2}(\Gamma_{\Sigma}))}$
$\displaystyle\qquad-2\rho\langle{\lambda}_{H,\Delta
T}^{(k)}-{\lambda}_{H,\Delta T},{u}_{\perp,h,\Delta
t}^{(k)}-{u}_{\perp,h,\Delta t}\rangle_{0,\Gamma_{C},(0,T]}\ .$

Using $\|{u}_{\perp,h,\Delta t}^{(k)}-{u}_{\perp,h,\Delta
t}\|_{L^{2}([0,T],L^{2}(\Gamma_{\Sigma}))}\leq\|{\bf{u}}_{h,\Delta
t}^{(k)}-{\bf{u}}_{h,\Delta t}\|_{L^{2}([0,T],L^{2}(\Gamma_{\Sigma}))^{d}}$,
we find that

$\displaystyle\|{\lambda}_{H,\Delta T}^{(k)}-{\lambda}_{H,\Delta
T}\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}^{2}-\|{\lambda}_{H,\Delta
T}^{(k+1)}-{\lambda}_{H,\Delta T}\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}^{2}$
$\displaystyle\geq 2\rho\langle{\lambda}_{H,\Delta
T}^{(k)}-{\lambda}_{H,\Delta T},{u}_{\perp,h,\Delta
t}^{(k)}-{u}_{\perp,h,\Delta
t}\rangle_{0,\Gamma_{C},(0,T]}-\rho^{2}\|{\bf{u}}_{h,\Delta
t}^{(k)}-{\bf{u}}_{h,\Delta
t}\|_{L^{2}([0,T],L^{2}(\Gamma_{\Sigma}))^{d}}^{2}\ .$

Further note

$\displaystyle\langle{\lambda}_{H,\Delta T}^{(k)}-{\lambda}_{H,\Delta
T},{u}_{\perp,h,\Delta t}^{(k)}-{u}_{\perp,h,\Delta
t}\rangle_{0,\Gamma_{C},(0,T]}$ $\displaystyle=\langle\mathcal{S}_{h,\Delta
t}({\bf{u}}_{h,\Delta t}^{(k)}-{\bf{u}}_{h,\Delta t}),{\bf{u}}_{h,\Delta
t}^{(k)}-{\bf{u}}_{h,\Delta t}\rangle_{0,\Gamma_{\Sigma},(0,T]}$
$\displaystyle\geq C\|{\bf{u}}_{h,\Delta t}^{(k)}-{\bf{u}}_{h,\Delta
t}\|_{L^{2}([0,T],L^{2}(\Gamma_{\Sigma}))^{d}}^{2}\ ,$

so that

$\|{\lambda}_{H,\Delta T}^{(k)}-{\lambda}_{H,\Delta
T}\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}^{2}-\|{\lambda}_{H,\Delta
T}^{(k+1)}-{\lambda}_{H,\Delta
T}\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}^{2}\geq(2\rho\,C-\rho^{2})\|{\bf{u}}_{h,\Delta
t}^{(k)}-{\bf{u}}_{h,\Delta
t}\|^{2}_{L^{2}([0,T],L^{2}(\Gamma_{\Sigma}))^{d}}\ .$

If $0<\rho<2C$ If $0<\rho<2C$, the right hand side is non-negative, and
$\|{\lambda}_{H,\Delta T}^{(k)}-{\lambda}_{H,\Delta
T}\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}$ is a decreasing sequence. As
$\|{\lambda}_{H,\Delta T}^{(k)}-{\lambda}_{H,\Delta
T}\|_{L^{2}([0,T],L^{2}(\Gamma_{C}))}\geq 0$, it converges, and we conclude
$\|{\bf{u}}_{h,\Delta t}^{(k)}-{\bf{u}}_{h,\Delta
t}\|_{L^{2}([0,T],L^{2}(\Gamma_{\Sigma}))^{d}}\to 0$. ∎

In addition to the Uzawa algorithm for the space-time problem, as above, also
Uzawa algorithms in each time-step were studied in [22]. They provide an
alternative solver for problem (25).

##  4 Error estimates

In this section we show an a priori estimate for the mixed formulation of the
contact problem (20) and its discretization (25). It builds on the
corresponding analysis for the scalar wave equation in [22] where, however,
refined information about the Poincaré-Steklov operator is available. Building
on the classical theory for time-domain boundary integral equations [31], the
analysis requires the formulation of (20) presented for times
$t\in\mathbb{R}^{+}$, using the generalized inner product (13) for $\sigma>0$
and considering

$\displaystyle{\widetilde{M}}^{+}:=\left\\{{\mu}\in
H^{1/2}_{\sigma}(\mathbb{R}^{+},\tilde{H}^{-1/2}(\Gamma_{C})):\left\langle{\mu},{w}\right\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}\geq
0,\ \forall\ 0\leq{w}\in
H^{-1/2}_{\sigma}(\mathbb{R}^{+},\tilde{H}^{1/2}(\Gamma_{\Sigma}))\right\\}\,.$
(26)

For $\textbf{f}\in
H^{1/2}_{\sigma}\left(\mathbb{R}^{+},H^{-1/2}(\Gamma_{\Sigma})\right)^{d}$ and
$g\in H^{1/2}_{\sigma}(\mathbb{R}^{+},H^{1/2}(\Gamma_{C}))$ the formulation is
given by:

Find $(\textbf{u},{\lambda})\in
H^{1/2}_{\sigma}(\mathbb{R}^{+},\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}\times{\widetilde{M}}^{+}$
such that

$\displaystyle\left\langle\mathcal{S}\textbf{u},\textbf{v}\right\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}-\left\langle\lambda,v_{\perp}\right\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$
$\displaystyle=\left\langle\textbf{f},\textbf{v}\right\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}$
$\displaystyle\forall\textbf{v}\in
H^{1/2}_{\sigma}(\mathbb{R}^{+},\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}$ (27a)
$\displaystyle\left\langle{\mu}-{\lambda},u_{\perp}\right\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$
$\displaystyle\geq\left\langle
g,\mu-\lambda\right\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$
$\displaystyle\forall{\mu}\in{\widetilde{M}}^{+}.$ (27b)

The corresponding discretization, generalizing (25), is given by:

Find $(\textbf{u}_{h,\Delta t},{\lambda}_{H,\Delta
T})\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t})^{d}\times
M^{+}_{H,\Delta T}$ such that

$\displaystyle\left\langle\mathcal{S}_{h,\Delta t}\textbf{u}_{h,\Delta
t},\textbf{v}_{h,\Delta
t}\right\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}-\left\langle{\lambda}_{H,\Delta
T},{v}_{\perp,h,\Delta
t}\right\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}=\left\langle\textbf{f},\textbf{v}_{h,\Delta
t}\right\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}\quad\forall\textbf{v}_{h,\Delta
t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t})^{d}$ (28a)
$\displaystyle\left\langle{\mu}_{H,\Delta T}-{\lambda}_{H,\Delta
T},{u}_{\perp,h,\Delta
t}\right\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}\geq\left\langle
g,\mu_{H,\Delta T}-\lambda_{H,\Delta
T}\right\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}\quad\quad\quad\quad\,\,\>\forall{\mu}_{H,\Delta
T}\in M^{+}_{H,\Delta T}\ ,$ (28b)

where the involved discrete functional spaces contain functions supported in a
finite number of time steps.
For the proof, we recall basic mapping properties of the Poincaré-Steklov
operator $\mathcal{S}$. They are summarized in the following Theorem 5. See
[14, 22] for corresponding results for the wave equation, with analogous
proofs.

######  Theorem 5.

For $r\in\mathbb{R}$ and $\sigma>0$,
$\mathcal{S}:H^{r}_{\sigma}(\mathbb{R}^{+},{H}^{\frac{1}{2}}(\Gamma))^{d}\to
H^{r}_{\sigma}(\mathbb{R}^{+},{H}^{-\frac{1}{2}}(\Gamma))^{d}$ continuously.

Here and below we write $A\lesssim B$ provided there exists a constant $C$
such that $A\leq CB$. If the constant $C$ is allowed to depend on a parameter
$\sigma$, we write $A\lesssim_{\sigma}B$.

Like in [22] an inf-sup estimate is a crucial ingredient to estimate the
error:

######  Theorem 6.

Let $C>0$ sufficiently small and $\frac{\max\\{h,\Delta t\\}}{\min\\{H,\Delta
T\\}}<C$. Then there exists $\alpha>0$ such that $\forall\,\lambda_{H,\Delta
T}\in X^{-1}_{H,\Gamma_{C}}\otimes V^{-1}_{\Delta T}$:

$\sup_{\mathbf{v}_{h,\Delta t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes
V^{0}_{\Delta t})^{d}}\frac{\langle v_{\perp,h,\Delta t},\lambda_{H,\Delta
T}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}}{\|\mathbf{v}_{h,\Delta
t}\|_{0,\frac{1}{2},\sigma,\ast}}\geq\alpha\|\lambda_{H,\Delta
T}\|_{0,-\frac{1}{2},\sigma}\ .$

The proof is found in [22, Theorem 15]. Now, we can prove the following
result:

######  Theorem 7.

Let $(\textbf{u},{\lambda})\in
H^{1/2}_{\sigma}(\mathbb{R}^{+},\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}\times{\widetilde{M}}^{+}$
be a solution to the mixed problem (27) and $(\textbf{u}_{h,\Delta
t},{\lambda}_{H,\Delta T})\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta
t})^{d}\times M^{+}_{H,\Delta T}$ a solution of the discretized mixed problem
(28). Assume that $\mathcal{S}$ is coercive. Then for a sufficiently small
constant $C>0$ and $\frac{\max\\{h,\Delta t\\}}{\min\\{H,\Delta T\\}}<C$, the
following a priori estimates hold:

$\displaystyle\|\lambda-\lambda_{H,\Delta T}\|_{0,-\frac{1}{2},\sigma}$
$\displaystyle\lesssim_{\sigma}\inf\limits_{\tilde{\lambda}_{H,\Delta T}\in
M^{+}_{H,\Delta T}}\|\lambda-\tilde{\lambda}_{H,\Delta
T}\|_{0,-\frac{1}{2},\sigma}+(\Delta t)^{-\frac{1}{2}}\|{\bf u}-{\bf
u}_{h,\Delta t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}\ ,$ (29)
$\displaystyle\|{\bf u}-{\bf u}_{h,\Delta
t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}$
$\displaystyle\lesssim_{\sigma}\inf\limits_{{\bf v}_{h,\Delta
t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t})^{d}}\|{\bf u}-{\bf
v}_{h,\Delta t}\|_{\frac{1}{2},\frac{1}{2},\sigma,\ast}$
$\displaystyle\qquad+\inf\limits_{\tilde{\lambda}_{H,\Delta T}\in
M^{+}_{H,\Delta T}}\left\\{\|\tilde{\lambda}_{H,\Delta
T}-\lambda\|_{\frac{1}{2},-\frac{1}{2},\sigma}+\|\tilde{\lambda}_{H,\Delta
T}-{\lambda}_{H,\Delta T}\|_{\frac{1}{2},-\frac{1}{2},\sigma}\right\\}\ .$
(30)

###### Proof.

Using equations (27) and (28), we note that for arbitrary
$\tilde{\lambda}_{H,\Delta T}\in M^{+}_{H,\Delta T}$ the following identity
holds:

$\displaystyle\langle\lambda_{H,\Delta T}-\tilde{\lambda}_{H,\Delta
T},v_{\perp,h,\Delta t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$
$\displaystyle=\langle\mathcal{S}\mathbf{u}_{h,\Delta t},\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}-\langle\mathbf{f},\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}-\langle\tilde{\lambda}_{H,\Delta
T},v_{\perp,h,\Delta t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$
$\displaystyle=\langle\mathcal{S}\mathbf{u}_{h,\Delta t},\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}-\langle\mathcal{S}\mathbf{u},\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}+\langle\lambda,v_{\perp,h,\Delta
t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$
$\displaystyle\qquad-\langle\tilde{\lambda}_{H,\Delta T},v_{\perp,h,\Delta
t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$
$\displaystyle=\langle\mathcal{S}(\mathbf{u}_{h,\Delta
t}-\mathbf{u}),\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}+\langle\lambda-\tilde{\lambda}_{H,\Delta
T},v_{\perp,h,\Delta t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}.$ (31)

From (4) and the inf-sup condition in Theorem 6, we obtain:

$\displaystyle\alpha\|\lambda_{H,\Delta T}-\tilde{\lambda}_{H,\Delta
T}\|_{0,-\frac{1}{2},\sigma}$ $\displaystyle\leq\sup_{\mathbf{v}_{h,\Delta
t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta
t})^{d}}\frac{\langle\lambda_{H,\Delta T}-\tilde{\lambda}_{H,\Delta
T},v_{\perp,h,\Delta
t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}}{\|\mathbf{v}_{h,\Delta
t}\|_{0,\frac{1}{2},\sigma,\ast}}$ $\displaystyle=\sup_{\mathbf{v}_{h,\Delta
t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta
t})^{d}}\frac{\langle\mathcal{S}(\mathbf{u}_{h,\Delta
t}-\mathbf{u}),\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}+\langle\lambda-\tilde{\lambda}_{H,\Delta
T},v_{\perp,h,\Delta
t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}}{\|\mathbf{v}_{h,\Delta
t}\|_{0,\frac{1}{2},\sigma,\ast}}\ .$

We separately estimate the two terms of the numerator. Using the continuity of
the duality pairing and an inverse inequality in time [23], the first term is
estimated by

$\displaystyle|\langle\mathcal{S}(\mathbf{u}_{h,\Delta
t}-u),\mathbf{v}_{h,\Delta t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}|$
$\displaystyle\leq\|\mathcal{S}(\mathbf{u}_{h,\Delta
t}-\mathbf{u})\|_{-\frac{1}{2},-\frac{1}{2},\sigma}\|{\mathbf{v}_{h,\Delta
t}}\|_{\frac{1}{2},\frac{1}{2},\sigma,\ast}$
$\displaystyle\lesssim\,\|\mathbf{u}_{h,\Delta
t}-\mathbf{u}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}(\Delta
t)^{-\frac{1}{2}}\|\mathbf{v}_{h,\Delta t}\|_{0,\frac{1}{2},\sigma,\ast}\ .$

For the second term, a similar argument yields:

$\displaystyle|\langle\lambda-\tilde{\lambda}_{H,\Delta T},v_{\perp,h,\Delta
t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}|$
$\displaystyle\leq\|\lambda-\tilde{\lambda}_{H,\Delta
T}\|_{0,-\frac{1}{2},\sigma}\|\mathbf{v}_{h,\Delta
t}\|_{0,\frac{1}{2},\sigma,\ast}\ .$

We conclude the a priori estimate as in (29), i.e.

$\displaystyle\|\lambda-\lambda_{H,\Delta T}\|_{0,-\frac{1}{2},\sigma}$
$\displaystyle\leq\inf\limits_{\tilde{\lambda}_{H,\Delta
T}}\left(\|\lambda-\tilde{\lambda}_{H,\Delta
T}\|_{0,-\frac{1}{2},\sigma}+\|{\lambda}_{H,\Delta
T}-\tilde{\lambda}_{H,\Delta T}\|_{0,-\frac{1}{2},\sigma}\right)$
$\displaystyle\lesssim_{\sigma}\inf\limits_{\tilde{\lambda}_{H,\Delta
T}}\|\lambda-\tilde{\lambda}_{H,\Delta T}\|_{0,-\frac{1}{2},\sigma}+(\Delta
t)^{-\frac{1}{2}}\|\mathbf{u}_{h,\Delta
t}-\mathbf{u}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}\ .$

Next we combine the Galerkin orthogonality

$\langle\mathcal{S}(\mathbf{u}-\mathbf{u}_{h,\Delta t}),\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}=\langle\lambda-{\lambda}_{H,\Delta
T},v_{\perp,h,\Delta t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$

with the coercivity of the Poincaré-Steklov operator, leading to

$\displaystyle\|\mathbf{u}_{h,\Delta t}-\mathbf{v}_{h,\Delta
t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}^{2}$
$\displaystyle\lesssim_{\sigma}\langle\mathcal{S}(\mathbf{u}_{h,\Delta
t}-\mathbf{v}_{h,\Delta t}),\mathbf{u}_{h,\Delta t}-\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}$
$\displaystyle=\langle\mathcal{S}(\mathbf{u}-\mathbf{v}_{h,\Delta
t}),\mathbf{u}_{h,\Delta t}-\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}+\langle\mathcal{S}(\mathbf{u}_{h,\Delta
t}-\mathbf{u}),\mathbf{u}_{h,\Delta t}-\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}$
$\displaystyle=\langle\mathcal{S}(\mathbf{u}-\mathbf{v}_{h,\Delta
t}),\mathbf{u}_{h,\Delta t}-\mathbf{v}_{h,\Delta
t}\rangle_{\sigma,\Gamma_{\Sigma},\mathbb{R}^{+}}$
$\displaystyle\qquad+\langle\tilde{\lambda}_{H,\Delta
T}-\lambda+{\lambda}_{H,\Delta T}-\tilde{\lambda}_{H,\Delta
T},u_{\perp,h,\Delta t}-v_{\perp,h,\Delta
t}\rangle_{\sigma,\Gamma_{C},\mathbb{R}^{+}}$

for all $\mathbf{v}_{h,\Delta t}$ and $\tilde{\lambda}_{H,\Delta T}$. The
mapping properties of $\mathcal{S}$ and the continuity of the duality pairing
then show

$\displaystyle\|\mathbf{u}_{h,\Delta t}-\mathbf{v}_{h,\Delta
t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}^{2}$
$\displaystyle\lesssim\|\mathbf{u}-\mathbf{v}_{h,\Delta
t}\|_{\frac{1}{2},\frac{1}{2},\sigma,\ast}\|\mathbf{u}_{h,\Delta
t}-\mathbf{v}_{h,\Delta t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}$
$\displaystyle\quad+\|\tilde{\lambda}_{H,\Delta
T}-\lambda\|_{\frac{1}{2},-\frac{1}{2},\sigma}\|\mathbf{u}_{h,\Delta
t}-\mathbf{v}_{h,\Delta t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}$
$\displaystyle\quad+\|\tilde{\lambda}_{H,\Delta T}-{\lambda}_{H,\Delta
T}\|_{\frac{1}{2},-\frac{1}{2},\sigma}\|\mathbf{u}_{h,\Delta
t}-\mathbf{v}_{h,\Delta t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}\ .$

Dividing by $\|\mathbf{u}_{h,\Delta t}-\mathbf{v}_{h,\Delta
t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}$, we find

$\displaystyle\|\mathbf{u}_{h,\Delta t}-\mathbf{v}_{h,\Delta
t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}$
$\displaystyle\lesssim_{\sigma}\|\mathbf{u}-\mathbf{v}_{h,\Delta
t}\|_{\frac{1}{2},\frac{1}{2},\sigma,\ast}+\|\tilde{\lambda}_{H,\Delta
T}-\lambda\|_{\frac{1}{2},-\frac{1}{2},\sigma}+\|\tilde{\lambda}_{H,\Delta
T}-{\lambda}_{H,\Delta T}\|_{\frac{1}{2},-\frac{1}{2},\sigma}\ .$

The triangle inequality now shows

$\displaystyle\|\mathbf{u}-\mathbf{u}_{h,\Delta
t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}$
$\displaystyle\lesssim_{\sigma}\|\mathbf{u}-\mathbf{v}_{h,\Delta
t}\|_{\frac{1}{2},\frac{1}{2},\sigma,\ast}+\|\tilde{\lambda}_{H,\Delta
T}-\lambda\|_{\frac{1}{2},-\frac{1}{2},\sigma}+\|\tilde{\lambda}_{H,\Delta
T}-{\lambda}_{H,\Delta T}\|_{\frac{1}{2},-\frac{1}{2},\sigma}\,,$

and estimate (30) follows. ∎

######  Remark 8.

Note that the assumption $(\textbf{u},{\lambda})\in
H^{1/2}_{\sigma}(\mathbb{R}^{+},\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}\times{\widetilde{M}}^{+}$
in Theorem 7 allows jump discontinuities in time of the velocity
$\dot{\mathbf{u}}\in
H^{-1/2}_{\sigma}(\mathbb{R}^{+},\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}$, as
they occur in the case of an impact.
The existence of a unique solution $(\textbf{u},{\lambda})\in
H^{1/2}_{\sigma}(\mathbb{R}^{+},\tilde{H}^{1/2}(\Gamma_{\Sigma}))^{d}\times{\widetilde{M}}^{+}$
and the coercivity of $\mathcal{S}$ are known, for example, in the case of the
analogous Signorini contact problem for the wave equation in simple geometries
for sufficiently smooth right hand side $\textbf{f}\in
H^{3/2}_{\sigma}\left(\mathbb{R}^{+},H^{-1/2}(\Gamma_{\Sigma})\right)^{d}$ and
$g=0$, see p. 450 in [14] or (in a slightly larger function space for more
general data) in [36].

######  Remark 9.

If the solution $(\mathbf{u},\lambda)$ is more regular, we obtain convergence
rates which are optimal in $h$, but suboptimal in $\Delta t$. Specifically,
let us illustrate the results when the time step $\Delta t$ is sufficiently
small, $H=Ch$, $\Delta T=C\Delta t$ and in the spatial variables $\mathbf{u}$
is of regularity $H^{s}(\Gamma_{\Sigma})$, $\lambda$ of regularity
$H^{s-1}(\Gamma_{C})$, where $s\in(\frac{1}{2},\frac{3}{2}]$. Then there exist
positive constants $C_{\sigma}^{\prime}$ and $C_{\sigma}^{\prime\prime}$ such
that

$\|\mathbf{u}-\mathbf{u}_{h,\Delta
t}\|_{-\frac{1}{2},\frac{1}{2},\sigma,\ast}\leq
C_{\sigma}^{\prime}h^{s-\frac{1}{2}}$

and

$\|\lambda-\lambda_{H,\Delta T}\|_{0,-\frac{1}{2},\sigma}\leq
C_{\sigma}^{\prime\prime}h^{s-\frac{1}{2}}.$

Detailed, but suboptimal convergence rates for general $\Delta t,\Delta T,h,H$
can be derived from the approximation results in [24, 26].

##  5 Algorithmic details

In this section, to keep the notation as simple as possible and because the
numerical results reported in Section 6 are related to two-dimensional
elastodynamics, we will fix $d=2$.

###  5.1 Implementation of the Poincaré-Steklov operator

The numerical solution of the mixed formulation (25) of the contact problem
will require the solution of the weak equation (25a), involving the
discretized operator $\mathcal{S}_{h,\Delta t}$, at every iteration of the
Uzawa algorithm, as described in Algorithm 1. Here, in particular, we are
interested in the algorithmic details related to the numerical solution of
such an equation.
For the implementation of the Poincaré-Steklov operator, both the symmetric
and the non-symmetric formulations (10), (11), in terms of the layer
potentials introduced in Section 2.2, will be used. While they are equivalent
for the continuous problem, the resulting discretizations will lead to
slightly different numerical results, as we will comment in Section 6.
Moreover, for the stability of the time-step scheme, we consider an energetic
weak formulation which involves time derivatives of test functions.
In general, for a given right-hand side $\widetilde{\textbf{f}}\in
H^{1/2}\left([0,T],H^{-1/2}(\Gamma_{\Sigma})\right)^{2}$ we have to solve a
weak problem of the form:

find $\textbf{u}_{h,\Delta t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes
V^{0}_{\Delta t})^{2}$ such that

$\langle\mathcal{S}_{h,\Delta t}\textbf{u}_{h,\Delta
t},\dot{\textbf{v}}_{h,\Delta
t}\rangle_{0,\Gamma_{\Sigma},(0,T]}=\langle\widetilde{\textbf{f}},\dot{\textbf{v}}_{h,\Delta
t}\rangle_{0,\Gamma_{\Sigma},(0,T]},\quad\quad\forall\,\textbf{v}_{h,\Delta
t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t})^{2}.$ (32)

We first consider the definition of the Poincaré-Steklov operator in its
symmetric form (11): because of its dependence on the inverse operator ${\cal
V}^{-1}$, to handle the equation at hand we need to define the dummy variable
$\boldsymbol{\psi}_{h,\Delta
t}:=\mathcal{V}^{-1}\left(\mathcal{K}+\frac{1}{2}\right)\textbf{u}_{h,\Delta
t}$ , belonging to the space $(X^{-1}_{h,\Gamma}\otimes V^{-1}_{\Delta
t})^{2}$.

This translates in solving at every step of the Uzawa algorithm the following
system of equations:

find $\textbf{u}_{h,\Delta t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes
V^{0}_{\Delta t})^{2}$ and $\boldsymbol{\psi}_{h,\Delta
t}\in(X^{-1}_{h,\Gamma}\otimes V^{-1}_{\Delta t})^{2}$ such that

$\begin{array}[]{l}\left\langle\mathcal{V}\boldsymbol{\psi}_{h,\Delta
t},\dot{\boldsymbol{\eta}}_{h,\Delta
t}\right\rangle_{0,\Gamma,(0,T]}-\left\langle\left(\mathcal{K}+\frac{1}{2}\right)\textbf{u}_{h,\Delta
t},\dot{\boldsymbol{\eta}}_{h,\Delta
t}\right\rangle_{0,\Gamma,(0,T]}=0,\quad\quad\quad\quad\quad\quad\quad\forall\,\boldsymbol{\eta}_{h,\Delta
t}\in(X^{-1}_{h,\Gamma}\otimes V^{-1}_{\Delta t})^{2}\\\\[4.0pt]
\left\langle\left(\mathcal{K}^{*}+\frac{1}{2}\right)\boldsymbol{\psi}_{h,\Delta
t},\dot{\textbf{v}}_{h,\Delta
t}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}-\left\langle\mathcal{W}\textbf{u}_{h,\Delta
t},\dot{\textbf{v}}_{h,\Delta
t}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}=\left\langle\widetilde{\textbf{f}},\dot{\textbf{v}}_{h,\Delta
t}\right\rangle_{0,\Gamma_{\Sigma},(0,T]},\forall\,\textbf{v}_{h,\Delta
t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t})^{2}.\end{array}$ (33)

Taking into account the definition of the discrete polynomial spaces (21),
(22) and (23) , let us now consider the sets
$\left\\{w^{(\boldsymbol{\psi})}_{m}\right\\}_{m=1}^{M^{(\boldsymbol{\psi})}_{h}}$
and $\left\\{w^{(\textbf{u})}_{m}\right\\}_{m=1}^{M^{(\textbf{u})}_{h}}$: the
first contains the piece-wise linear basis functions of $X^{-1}_{h,\Gamma}$,
while the second is composed by the piece-wise linear basis functions of
$X^{0}_{h,\Gamma_{\Sigma}}$. The algebraic reformulation of (33) can be then
elaborated choosing the approximate components of the unknowns, in the spaces
$X^{-1}_{h,\Gamma}\otimes V^{-1}_{\Delta t}$ and
$X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t}$, of the following form:

$\psi_{i,h,\Delta t}(\textbf{x},t)=\sum_{\ell=0}^{N_{\Delta
t}-1}\sum_{m=1}^{M^{(\boldsymbol{\psi})}_{h}}\psi_{i,\ell,m}w^{(\boldsymbol{\psi})}_{m}(\textbf{x})v_{\ell}(t),\>\>u_{i,h,\Delta
t}(\textbf{x},t)=\sum_{\ell=0}^{N_{\Delta
t}-1}\sum_{m=1}^{M^{(\textbf{u})}_{h}}u_{i,\ell,m}w^{(\textbf{u})}_{m}(\textbf{x})r_{\ell}(t),\>\>i=1,2,$
(34)

where the time basis $v_{k}$ and $r_{k}$ are defined as

$v_{\ell}(t):=H[t-t_{\ell}]-H[t-t_{\ell+1}],\quad
r_{\ell}(t):=H[t-t_{\ell}]\frac{t-t_{\ell}}{\Delta
t}-H[t-t_{\ell+1}]\frac{t-t_{\ell+1}}{\Delta t}.$

These choices for the approximation of $\textbf{u}_{h,\Delta
t},\boldsymbol{\psi}_{h,\Delta t}$ leads to the algebraic reformulation of
(33) as a linear system $\mathbb{S}\textbf{X}=\widetilde{\textbf{F}}$, having
the form

$\left(\begin{array}[]{cccc}\mathbb{S}^{(0)}&\textbf{0}&\cdots&\textbf{0}\\\
\mathbb{S}^{(1)}&\mathbb{S}^{(0)}&\cdots&\textbf{0}\\\
\vdots&\vdots&\ddots&\vdots\\\ \mathbb{S}^{(N_{\Delta
t}-1)}&\mathbb{S}^{(N_{\Delta
t}-2)}&\cdots&\mathbb{S}^{(0)}\end{array}\right)\left(\begin{array}[]{c}\textbf{X}_{(0)}\\\
\textbf{X}_{(1)}\\\ \vdots\\\ \textbf{X}_{(N_{\Delta
t}-1)}\end{array}\right)=\left(\begin{array}[]{c}\widetilde{\textbf{F}}_{(0)}\\\
\widetilde{\textbf{F}}_{(1)}\\\ \vdots\\\ \widetilde{\textbf{F}}_{(N_{\Delta
t}-1)}\end{array}\right),$ (35)

where, for all time indices $\ell=0,...,N_{\Delta t}-1$, the unknown vectors
are structured as

$\textbf{X}_{(\ell)}=(\boldsymbol{\Psi}_{(\ell)},\textbf{U}_{(\ell)})=\left(\psi_{1,\ell,1},\ldots,\psi_{1,\ell,M_{h}^{(\boldsymbol{\psi})}},\psi_{2,\ell,1},\ldots,\psi_{2,\ell,M_{h}^{(\boldsymbol{\psi})}},u_{1,\ell,1},\ldots,u_{1,\ell,M_{h}^{(\textbf{u})}},u_{2,\ell,1},\ldots,u_{2,\ell,M_{h}^{(\textbf{u})}}\right)^{\top}\,$

and the right-hand side as

$\widetilde{\textbf{F}}_{(\ell)}=\left(0,\cdots,0,\,0,\cdots,0,\,\widetilde{F}_{1,\ell,1},\ldots,\widetilde{F}_{1,\ell,M_{h}^{(\textbf{u})}},\widetilde{F}_{2,\ell,1},\ldots,\widetilde{F}_{2,\ell,M_{h}^{(\textbf{u})}}\right)^{\top}\in\mathbb{R}^{2(M_{h}^{(\boldsymbol{\psi})}+M_{h}^{(\textbf{u})})}\,,$

with
$\widetilde{F}_{i,\ell,m}:=\left\langle\widetilde{f}_{i},w^{(\textbf{u})}_{m}\dot{r}_{\ell}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}$.
The matrix $\mathbb{S}$ shows a lower triangular Toeplitz structure, so that
the system (35) can be solved by backsubstitution. This leads to a marching-
on-in-time time stepping scheme, in which the solution of a system of
equations involving the first time block $\mathbb{S}^{(0)}$ is required in
every time step.

The structure of the linear system (35) and the expression of the matrix
entries as double spatial integrals are due to the analytical time integration
taken into account in formula (33). Let us note that the structure of a
generic block $\mathbb{S}^{(\ell)}$ has dimension
$2(M_{h}^{(\boldsymbol{\psi})}+M_{h}^{(\boldsymbol{\textbf{u}})})\times
2(M_{h}^{(\boldsymbol{\psi})}+M_{h}^{(\boldsymbol{\textbf{u}})})$ and assumes
the form:

$\mathbb{S}^{(0)}=\left(\begin{array}[]{cc}\mathbb{V}^{(0)}&-(\mathbb{K}^{(0)}+1/2\>\mathbb{M})\\\\[4.0pt]
({\mathbb{K}}^{(0)}+1/2\>\mathbb{M})^{\top}&-\mathbb{W}^{(0)}\end{array}\right),\quad\mathbb{S}^{(\ell)}=\left(\begin{array}[]{cc}\mathbb{V}^{(\ell)}&-\mathbb{K}^{(\ell)}\\\\[4.0pt]
{\mathbb{K}^{(\ell)}}^{\top}&\mathbb{W}^{(\ell)}\end{array}\right),\>\ell=1,...,N_{\Delta
t}-1.$

With obvious meaning of notation, the blocks $\mathbb{V}^{(\ell)}$,
$\mathbb{K}^{(\ell)}$ and $\mathbb{W}^{(\ell)}$ are obtained from the
discretization of the integral operators involved in (33) and $\mathbb{M}$ is
a mass matrix. For details about the numerical evaluation the entries of these
blocks the reader is referred to [15, 2].

We briefly also recall the problem we need to solve in case of non-symmetric
formulation of the Poincaré-Steklov operator (10). With similar arguments used
for the symmetric formulation, we obtain the following system of weak boundary
equations:

find $\textbf{u}_{h,\Delta t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes
V^{0}_{\Delta t})^{2}$ and $\boldsymbol{\psi}_{h,\Delta
t}\in(X^{-1}_{h,\Gamma}\otimes V^{-1}_{\Delta t})^{2}$ such that

$\begin{array}[]{l}\left\langle\mathcal{V}\boldsymbol{\psi}_{h,\Delta
t},\dot{\boldsymbol{\eta}}_{h,\Delta
t}\right\rangle_{0,\Gamma,(0,T]}-\left\langle\left(\mathcal{K}+\frac{1}{2}\right)\textbf{u}_{h,\Delta
t},\dot{\boldsymbol{\eta}}_{h,\Delta
t}\right\rangle_{0,\Gamma,(0,T]}=0,\quad\forall\,\boldsymbol{\eta}_{h,\Delta
t}\in(X^{-1}_{h,\Gamma}\otimes V^{-1}_{\Delta t})^{2}\\\\[4.0pt]
\left\langle\boldsymbol{\psi}_{h,\Delta t},\dot{\textbf{v}}_{h,\Delta
t}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}=\left\langle\widetilde{\textbf{f}},\dot{\textbf{v}}_{h,\Delta
t}\right\rangle_{0,\Gamma_{\Sigma},(0,T]},\quad\quad\quad\quad\quad\quad\quad\quad\forall\,\textbf{v}_{h,\Delta
t}\in(X^{0}_{h,\Gamma_{\Sigma}}\otimes V^{0}_{\Delta t})^{2}.\end{array}$ (36)

The discrete unknowns $\textbf{u}_{h,\Delta},\boldsymbol{\psi}_{h,\Delta}$ is
substituted again with (34), leading to an algebraic form
$\mathbb{S}\textbf{X}=\widetilde{\textbf{F}}$ of (36), analogous to the system
(35), with the only difference that the blocks $\mathbb{S}^{(\ell)}$, again
with dimension
$2(M_{h}^{(\boldsymbol{\psi})}+M_{h}^{(\boldsymbol{\textbf{u}})})\times
2(M_{h}^{(\boldsymbol{\psi})}+M_{h}^{(\boldsymbol{\textbf{u}})}),$ are
structured as

$\mathbb{S}^{(0)}=\left(\begin{array}[]{cc}\mathbb{V}^{(0)}&-(\mathbb{K}^{(0)}+1/2\>\mathbb{M})\\\\[4.0pt]
\mathbb{M}^{\top}&\textbf{0}\end{array}\right),\quad\mathbb{S}^{(\ell)}=\left(\begin{array}[]{cc}\mathbb{V}^{(\ell)}&-\mathbb{K}^{(\ell)}\\\\[4.0pt]
\textbf{0}&\textbf{0}\end{array}\right),\>\ell=1,...,N_{\Delta t}-1.$

###  5.2 Algebraic formulation of the Uzawa algorithm

For the numerical solution of contact problems, the right-hand side
$\widetilde{\textbf{f}}$ in (33) and (36) specifies depending on the Neumann
datum f and the discretized Lagrange multiplier $\lambda_{H,\Delta
T}(\textbf{x},t)$ as in (25a).
Hence, let us consider the vector ${\bf F}=\left({\bf
F}_{(\ell)}\right)_{\ell=0}^{N_{\Delta t}-1}$ with

$\textbf{F}_{(\ell)}=\left(0,\cdots,0,\,0,\cdots,0,\,F_{1,\ell,1},\ldots,F_{1,\ell,M_{h}^{(\textbf{u})}},F_{2,\ell,1},\ldots,F_{2,\ell,M_{h}^{(\textbf{u})}}\right)^{\top}\in\mathbb{R}^{2(M_{h}^{(\boldsymbol{\psi})}+M_{h}^{(\textbf{u})})}$

and $F_{i,\ell,m}:=\left\langle
f_{i},w^{(\textbf{u})}_{m}\dot{r}_{\ell}\right\rangle_{0,\Gamma_{\Sigma},(0,T]}$.
For the discretized Lagrange multiplier related to the simulations presented
in Section 6, in practice we have chosen $H=h$ and $\Delta T=\Delta t$ and we
have extended it trivially in vector form, in such a way that its Cartesian
components can be expressed as

$\lambda_{i,h,\Delta t}(\textbf{x},t)=\sum_{\ell=0}^{N_{\Delta
t}-1}\sum_{m=1}^{M^{(\boldsymbol{\lambda})}_{h}}\lambda_{i,\ell,m}w^{(\boldsymbol{\lambda})}_{m}(\textbf{x})v_{\ell}(t),\;i=1,2,$
(37)

with
$\left\\{w^{(\boldsymbol{\lambda})}_{m}\right\\}_{m=1}^{M^{(\boldsymbol{\lambda})}_{h}}$
the piece-wise constant basis functions of $X^{-1}_{h,\Gamma_{C}}$. For each
time step $\ell=0,...,N_{\Delta t}-1$, we introduce the vectors

$\boldsymbol{\Lambda}_{(\ell)}=\left(\lambda_{1,\ell,1},\ldots,\lambda_{1,\ell,M_{h}^{(\boldsymbol{\lambda})}},\lambda_{2,\ell,1},\ldots,\lambda_{2,\ell,M_{h}^{(\boldsymbol{\lambda})}}\right)^{\top}$

which are collected in the vector
$\boldsymbol{\Lambda}=\left(\boldsymbol{\Lambda}_{(\ell)}\right)_{\ell=0}^{N_{\Delta
t}-1}\in\mathbb{R}^{N_{\Delta t}2M_{h}^{(\boldsymbol{\lambda})}}$. We denote
by $\mathcal{J}_{\perp}$ the set of those indices $j$ of the vector
$\boldsymbol{\Lambda}$ which correspond to the (nontrivial) components of
$\boldsymbol{\lambda}_{h,\Delta t}$ normal to $\Gamma$, and by
$\mathcal{J}_{\parallel}$ the set of the remaining indices, corresponding to
the (trivial) components of $\boldsymbol{\lambda}_{h,\Delta t}$ tangential to
$\Gamma$. The Uzawa algorithm given in Algorithm 1 then translates into the
following algebraic procedure where the stopping criterion is specified:

Algorithm 2 (Algebraic formulation of Uzawa algorithm)

Fix $\rho>0$ and $\epsilon>0$.

$k=0$, $\boldsymbol{\Lambda}^{(0)}=\textbf{0}$ and
$\boldsymbol{\Lambda}^{(-1)}=\textbf{1}$

while
$\|\boldsymbol{\Lambda}^{(k)}-\boldsymbol{\Lambda}^{(k-1)}\|_{2}/\|\boldsymbol{\Lambda}^{(k)}\|_{2}>\epsilon$
do

solve
$\mathbb{S}\textbf{X}^{(k)}=\textbf{F}+\mathbb{M}^{*}\boldsymbol{\Lambda}^{(k)}$

extract  $\textbf{U}^{(k)}$ from $\textbf{X}^{(k)}$

compute
$\boldsymbol{\Lambda}^{(k+1)}=\text{pr}_{C}(\boldsymbol{\Lambda}^{(k)}-\rho\widetilde{\mathbb{M}}(\textbf{U}^{(k)}-\textbf{G}))$.

$k\leftarrow k+1$

end while

where $\text{pr}_{C}:\mathbb{R}^{N_{\Delta
t}2M_{h}^{(\boldsymbol{\lambda})}}\rightarrow\mathbb{R}^{N_{\Delta
t}2M_{h}^{(\boldsymbol{\lambda})}}$ is understood as the discretized version
of the projector considered in Algorithm 1, acting on a vector ${\bf W}$ as

$(\text{pr}_{C}\ {\bf
W})_{j}=\left\\{\begin{array}[]{lr}\max\,\\{W_{j},0\\},&j\in\mathcal{J}_{\perp}\\\
0,&j\in\mathcal{J}_{\parallel}\end{array}\,.\right.$ (38)

The matrix $\mathbb{M}^{*}\in\mathbb{R}^{N_{\Delta
t}\,2(M_{h}^{(\boldsymbol{\psi})}+M_{h}^{(\textbf{u})})\times N_{\Delta
t}2M_{h}^{(\boldsymbol{\lambda})}}$ represents the projection of discretized
Lagrange multipliers on the discretized displacement space, trivially extended
in such a way that the vector $\mathbb{M}^{*}\boldsymbol{\Lambda}^{(k)}$
matches the length of vector ${\bf F}$; the vector G contains the coefficients
of the interpolant of $g$ in $X^{0}_{h,\Gamma_{C}}\otimes V^{0}_{\Delta t}$,
suitably trivially extended to match the length of vector $\textbf{U}^{(k)}$,
while $\widetilde{\mathbb{M}}\in\mathbb{R}^{N_{\Delta
t}2M_{h}^{(\boldsymbol{\lambda})}\times N_{\Delta t}2M_{h}^{(\textbf{u})}}$ is
the mass matrix representing the interplay between the finite dimensional
spaces of discretized Lagrange multipliers and discretized displacements.

##  6 Numerical results

In the following we present numerical experiments which apply the proposed
methods to the approximation of two-dimensional ($d=2$) dynamic Signorini
contact problems. Different boundary geometries (slit, square boundary and
circumference) and contact conditions are considered, illustrating the
stability and efficiency of the approach.
Unless otherwise specified, the gap function in contact conditions (3) and the
mass density are respectively set to be $g=0$ and $\varrho=1$.

To study the convergence of the numerical solutions, we consider the error in
energy norm. Given a solution vector $\textbf{X}^{(k)}$ to the linear system
$\mathbb{S}\textbf{X}^{(k)}=\textbf{F}+\mathbb{M}^{*}\boldsymbol{\Lambda}^{(k)}$
in Algorithm 2, once the stopping criterion in this algorithm is satisfied, we
compute the square of the energy norm as
$\textbf{X}^{(k)\top}\mathbb{S}\textbf{X}^{(k)}$. In the examples where the
exact solution is unknown, we compare these energies with an extrapolated
benchmark value.

###  6.1 Example 1: non-symmetric and symmetric formulation of the Poincaré-
Steklov operator

In the time interval $[0,T]=[0,2]$, let us consider the square centred in the
origin of the fixed orthogonal reference system and with sides of unitary
length, i.e. $\Omega=[-0.5,0.5]^{2}$. We denote the square boundary by
$\Gamma:=\partial\Omega=\Gamma_{b}\cup\Gamma_{r}\cup\Gamma_{t}\cup\Gamma_{l}$,
with obvious meaning of notation for the bottom, right, top and left sides.
The elastodynamics velocities are fixed as $c_{S}=1/\sqrt{2}$ and $c_{P}=1$:
in this way the problem at hand results decoupled w.r.t. the horizontal and
vertical directions [18].

Here the partition of the boundary is trivial: $\Gamma=\Gamma_{N}$. Ipso
facto, we are dealing with a dynamics of interior elastic wave propagation
with Neumann boundary conditions, here determined by the vertical datum
$\textbf{f}=(0,f_{2})$, with

$f_{2}(\textbf{x},t)=\left\\{\begin{array}[]{c
l}0,&\textbf{x}\in\Gamma_{l}\cup\Gamma_{r}\\\ 1,&\textbf{x}\in\Gamma_{t}\\\
-2\,H(t-1),&\textbf{x}\in\Gamma_{b}\end{array}\right..$

Since $\Gamma_{C}=\emptyset$ the implementation of Algorithm 2 is trivial,
converging in fact in one iteration. For this assigned Neumann datum, the
displacement solution of (1) is analytically known [18] and reads:

$\begin{array}[]{l}u_{1}(\textbf{x},t)=0,\\\
u_{2}(\textbf{x},t)=H(c_{P}\,t-x_{2}+0.5)\,(c_{P}\,t-x_{2}+0.5)-H(c_{P}\,t-x_{2}-1.5)\,(c_{P}\,t-x_{2}-1.5)\,.\end{array}$
(39)

The related discrete energetic weak boundary integral problem (33) or (36) is
solved setting uniform meshes in space and time, generated by fixing space and
time steps of equal value: $h=\Delta t$.

Making use of the non-symmetric formulation of the Poincaré-Steklov operator
(10) and setting $h=\Delta t=0.05$, we obtain a good approximation
$u_{1,h,\Delta t}$ of the trivial horizontal displacement. The approximation
of the vertical displacement, in Figure 2(a), is instead more intriguing. The
comparison between the numerical solution $u_{2,h,\Delta t}$ and the
corresponding analytical one in the midpoint of the top, left (right) and
bottom sides of the square boundary is optimal, given that the respective
three lines are overlapped. In Figure 2(b) the surface of $u_{2,h,\Delta t}$
on the entire space-time domain is reported.

(a)  (b)

Figure 2: Comparison between $u_{2,h,\Delta t}$ and the corresponding
analytical solution (a) and $u_{2,h,\Delta t}$ surface on the entire space-
time domain (b).

Concerning Figure 3, results about the error analysis are collected,
considering a resolution of the problem by both non-symmetric and symmetric
couplings (see (10) and (11)). This analysis is partially summarized in Figure
3(a), which shows the space-time error in $L^{2}((0,T)\times\Gamma)$-norm
between the discrete components of the solution $u_{1,h,\Delta
t},u_{2,h,\Delta t}$ and the analytical ones reported in (39). For both types
of coupling, the errors behave as $\mathcal{O}(h^{1.5})$.
Considering the same levels of space-time discretization, Figure 3(b) collects
the squared energy error committed with both non-symmetric and symmetric
couplings of the problem: this error is characterized by a linear linear
slope, behaving in fact as $\mathcal{O}(h)$.

(a)  (b)

Figure 3: $L^{2}((0,T)\times\Gamma)$ space-time error committed w.r.t to the
analytical solution (39) (a) and squared error in energetic norm (b).

###  6.2 Example 2: dynamic contact problem on a segment

Here we consider the slit
$\Gamma=\left\\{(x_{1},0)\>|\>x_{1}\in[-0.5,0.5]\right\\}$ with a contact
region identified by the sub-interval
$\Gamma_{C}=\left\\{(x_{1},0)\>|\>x_{1}\in[-0.2,0.2]\right\\}$. For this
example we set $\Gamma_{\Sigma}=\Gamma_{C}$ and, since on
$\Gamma_{D}=\Gamma\setminus\Gamma_{\Sigma}$ the datum $\textbf{u}=0$ is
imposed, the slit results partially fixed on both sides and the dynamics
purely depends on an assigned contact force f. The considered wave speeds are
$c_{S}=1,c_{P}=2$ and the time interval of analysis is $[0,T]=[0,3]$.
Two different vertical contact forces f have been considered on $\Gamma_{C}$,
fixing for both a null horizontal impulse $f_{1}=0$:

$f_{2}(t)=\left\\{\begin{array}[]{l r}-\sin(8\pi
t)H(0.1875-t)+H(t-0.1875)&{\rm for\,Test\,1},\\\ -\sin(16\pi
t)H(0.21875-t)+H(t-0.21875)&{\rm for\,Test\,2}.\\\ \end{array}\right.$

Note that the slit $\Gamma$ considered in this Example does not quite satisfy
the assumptions described in Section 2.1. It approximates the contact between
an infinite upper half-space and a rigid lower half-space. By truncating the
infinitely extended interface to a finite segment $\Gamma$, we incur geometric
truncation errors in the numerical discretizations of $\mathcal{S}$.
The algebraic Uzawa Algorithm 2 is used with a tolerance $\epsilon=10^{-5}$,
while the parameter for the projector $\textrm{pr}_{C}$ has been chosen as
$\rho=10^{5}$. The discretization space-time steps can assume different values
but maintaining always the ratio $h/\Delta t=2$. The numerical results are
compared against an extrapolated benchmark for the symmetric, respectively
non-symmetric coupling.
Setting the space step $h=0.025$, in Figure 4 we report the corresponding
vertical component of displacement $u_{2,h,\Delta t}$ in the midpoint of the
contact region $\Gamma_{C}$, together with the vertical component of the
contact force f taken into account. This force, when positive, causes a not
trivial displacement solution on $\Gamma_{C}$, since the contact region is
pulled up by the contact force. Otherwise, when the vertical contact force is
negative, the resultant vertical displacement decreases to $0$ and the contact
region turns to its original flat form since the vertical displacement cannot
become negative. In this Figure, $\lambda_{2,h,\Delta t}$ is also shown.
With reference to Test 1, having fixed the same discretization parameters as
before, in Figure 5 the deformation of the contact region is represented
through the corresponding approximated horizontal and vertical components of
${\bf u}_{h,\Delta t}$ at three time instants when the prescribed vertical
force is positive. With reference to Test 2, Figure 6 shows the approximated
$u_{1,h,\Delta t}$ and $u_{2,h,\Delta t}$ at three time instants at the
beginning of the simulation. Looking in particular at the vertical component,
at first the displacement remains trivial, then it starts assuming positive
values, before decreasing, with more evidence towards the endpoints of the
contact region, when the vertical contact force becomes negative again.
In Figure 7, the approximation $\lambda_{2,h,\Delta t}$ is shown, both for
Test 1 and Test 2, on the global space-time domain $\Gamma_{C}\times[0,3]$,
having set $h=0.025$. As one can see, in Figure 7 (a), related to Test 1, the
non-trivial part of the surface corresponds to a nonzero contact force for
$t\in[0,0.125]$ to compensate the negative forcing term $f_{2}(t)$ applied
uniformly on $\Gamma_{C}$. The shape of $\lambda_{2,h,\Delta t}$ surface given
by Test 2, Figure 7 (b), presents similar features in the initial phase, but,
at the second negative pulsation of the contact force, $\lambda_{2,h,\Delta
t}$ is non-trivial just towards the endpoints of $\Gamma_{C}$, since in the
middle of the segment the vertical displacement has a decreasing, but still
positive value.
Let us remark that both ${\bf u}_{h,\Delta t}$ and $\lambda_{h,\Delta t}$ are
retrieved from the final Uzawa iteration of Algorithm 2, once the stopping
criterion is satisfied.

(a)  (b)

Figure 4: Time history of $u_{2,h,\Delta t}$ in the midpoint of $\Gamma_{C}$
for Test 1 (left), Test 2 (right) together with the corresponding vertical
contact force and $\lambda_{2,h,\Delta t}$. Figure 5: Vertical and horizontal
deformation of $\Gamma_{C}$ at three different time instants (Test 1). Figure
6: Vertical and horizontal deformation of $\Gamma_{C}$ at three different time
instants (Test 2).

(a)  (b)

Figure 7: Surface of $\lambda_{2,h,\Delta t}$ for Test 1 (a) and Test 2 (b) on
the entire space-time domain $\Gamma_{C}\times[0,3]$.

The corresponding error analysis for this experiment is reported in Figure 8:
the squared error in energy norm is indicated distinguishing between Test 1
and Test 2 (Figure 8(a)) and (Figure 8(b)) and for a subsequent refinement of
the $h$ step. For both experiments the straight line in the plots correspond
to convergence of order $\mathcal{O}(h)$. Following Remark 9, the convergence
rate is compatible with a solution $\mathbf{u}$ of regularity
$H^{1-\varepsilon}(\Gamma_{\Sigma})$ in the spatial variables.

(a)  (b)

Figure 8: Squared error in energy norm related to the Example 2 (Test 1 top,
Test 2 bottom).

###  6.3 Example 3: dynamic contact problem in non-flat geometry

For this example, we take into account again the square
$\Omega=[-0.5,0.5]^{2}$, with $\Gamma_{N}=\Gamma_{t}\cup\Gamma_{r}$ and the
same physical and discretization parameters considered in the Example 1. Let
us define the contact region as $\Gamma_{C}=\Gamma_{b}\cup\Gamma_{l}$. The
dynamic will be investigated using two contact forces, always fixing
$\textbf{f}=\textbf{0}$ on $\Gamma_{C}$, $f_{2}(\textbf{x},t)=-0.1H[t]$ on
$\Gamma_{t}$, $f_{2}(\textbf{x},t)=0$ on $\Gamma_{r}$ and

$f_{1}(\textbf{x},t)=\left\\{\begin{array}[]{c r}0\;{\rm
on}\,\Gamma_{t}\cup\Gamma_{r},&{\rm for\,Test\,1}\\\ -0.1H[t]\,{\rm
on}\,\Gamma_{r};\;0\,{\rm on}\,\Gamma_{t},&{\rm for\,Test\,2}\\\
\end{array}\right..$

Figure 9 perfectly shows how the two types of contact forces interact with the
square: in Test 1 this is pushed down from the top side, while for Test 2 the
square is pressed both from the top and from the left side.

The Uzawa algorithm has been set up with an exit test fixing
$\epsilon=10^{-5}$, while the parameter for the projector has been fixed as
$\rho=10^{2},10^{3},10^{4}$ for decreasing mesh parameters.

(a)

(b)

Figure 9: Schematic representation of the two types of force having role in
Test 1 and Test 2, respectively (a) and (b), in Example 3.

Using the non-symmetric formulation of $\mathcal{S}$ given in (10), having
fixed $h=\Delta t=0.05$ we obtain the approximate displacements shown in
Figure 10 related to Test 1 and in Figure 11 for Test 2.

Figure 10: Time history of $\textbf{u}_{h,\Delta t}$ components in the
midpoint of $\Gamma$ sides for Test 1, together with the corresponding
vertical Neumann datum.

Figure 11: Time history of $\textbf{u}_{h,\Delta t}$ components in the
midpoint of $\Gamma$ sides for Test 2, together with the corresponding
vertical Neumann datum.

Figure 12 shows the squared energy error decay for $h$-refinements. For both
Test 1 and Test 2, the considered representations (10) and (11) of the
Poincaré-Steklov operator lead to a linear slope for the squared error,
parallel to the line $\mathcal{O}(h)$. As in Example 2, the convergence rate
is compatible with a solution $\mathbf{u}$ of regularity
$H^{1-\varepsilon}(\Gamma_{\Sigma})$ in the spatial variables.

(a)  (b)

Figure 12: Squared error in energy norm related to the Example 3: figures (a)
and (b) include errors obtained by non-symmetric and symmetric coupling of the
problem.

With the aim of studying the conservation of the energy after multiple
bounces, let us consider the configuration depicted in Figure 9 (a). The
elastic properties take the same values as above, but we change the applied
Neumann vertical datum at the top side of the square boundary, namely

$f_{2}({\bf
x},t)=0.1\;\sum_{k=0}^{8}(-1)^{k}\left(H[t-0.1\,k]-H[t-0.1\,(k+1)]\right),\quad{\bf
x}\in\Gamma_{t}.$

For fixed $h=\Delta t=0.0125$, Figure 13 shows the time history of the
displacements components in the midpoints of the top (a) and bottom (b) sides
of the square boundary, together with $f_{2}({\bf x},t)$ in the same points.
In particular, in Figure 13 (b), after a delay due to the time for the elastic
wave to reach the bottom of the square from the top side, multiple bounces
against the contact region are highlighted. No spurious oscillations are
observed even after multiple contacts. For the same discretization parameters,
Figure 14 (a) shows the energy of the elastic system as a function of time.
The approximately linear increase for times $<1$ and constant behavior for
$t>1$ agrees with the depicted behavior of the solution and the applied
forces. Figure 14 (b) zooms in around the stationary behavior. It shows that
the energy is conserved up to a relative error of less than $10^{-4}$ for
$t>1$.

(a)  (b)

Figure 13: Displacement as a function of time in the midpoints of the top (a)
and bottom (b) sides, together with the assigned vertical Neumann datum
$f_{2}$.

(a)  (b)

Figure 14: Energy of the numerical solution as a function of time (a), zoomed-
in plot for $t>0.8$ (b).

###  6.4 Example 4: dynamic contact problem in a circular geometry

In the last example we take into account the disk
$\Omega=\left\\{(x,y)^{\top}\>:\>x^{2}+y^{2}\leq 0.2\right\\}$. The contact
region is given by the entire boundary of the disk, namely
$\Gamma=\Gamma_{C}$. The fundamental elastodynamic velocities are set as
$c_{S}=1$ and $c_{P}=2$. We have no Neumann data assigned and the contact
force is null, but the dynamic is not trivial since we consider in (3) the gap
function

$g(t)=(4\>\sqrt{1-1.5(t-0.5)^{2}}-4)\>H[1.3-t]-3.2\>H[t-1.3],$

modeling a contact that push up the disk from the bottom in the time interval
$[0.245,0.755]$. A schematic representation of the dynamic contact is reported
in Figure 15.

Figure 15: Schematic representation of disk and contact boundary: its movement
in time is governed by the given gap function $g$.

Using the symmetric formulation (11), we show in Figure 16 the squared energy
error, which decays with the halving of the parameters $h$ and $\Delta t$ (for
each level of discretization the ratio $h/\Delta t\simeq 2$ is maintained).
The exit test for stopping the Uzawa algorithm consists in fixing
$\epsilon=10^{-4}$, and we set the parameter $\rho=10^{4}$ in the Uzawa
update.

Figure 16: Squared error in energy norm for the test on the disk (symmetric
formulation).

Having fixed $h\simeq 0.02$, Figure 17(a) shows the time history of the
vertical coordinates of the bottom point of $\Omega$, $B$, and its barycenter,
$C$, that, before the contact, are identified respectively with the Cartesian
coordinates $(0,-0.2)$ and $(0,0)$. It is immediate to observe that during the
contact, the bottom point behaves exactly like the gap function $g(t)$, which
physically represents a lower bound for the vertical coordinates of the points
belonging to $\Omega$, that are dragged up and then released at the end of the
contact. Moreover, the approximated displacement at the barycenter, as
expected, is characterized by an initial exponential growth and then it
becomes linear in time, since after the contact there are no longer external
forces applied to the disk $\Omega$. Figure 17(b) represents instead the
global displacement to which the disk is subjected at some time instants:
after an initial phase of compression, the points globally move upwards and
the deformations are related to the elastic constants that characterize the
material of the disk.

(a)  (b)

Figure 17: Time history of the bottom and the barycenter points $B$ and $C$
(Figure (a)) and global displacement of the disk $\Omega$ at some time
instants (Figure (b)).

##  7 Conclusions

In this work we propose and analyze a Galerkin space-time boundary element
method to solve elastodynamic contact problems. Boundary elements provide a
natural and efficient formulation, as the contact takes place at the interface
between two bodies.
The article presents this approach in the case of a unilateral Signorini
contact problem without friction. In particular, the algebraic formulation and
implementation of an energetic space-time boundary element method are detailed
for a mixed formulation of the problem. Analytically we obtain an a priori
error analysis, based on ideas for a scalar variational inequality [22]. The
analysis of the method is crucially based on an inf–sup condition in
space–time. Moreover, the space-time Uzawa iterative algorithm, used to solve
the nonlinear problem, is shown to be provably convergent.
As a key point, the numerical results indicate stability and convergence for
different two-dimensional geometries, including both straight and curved
contact boundaries. Future work will extend the current approach to two-sided
and frictional contact problems, to three-dimensional geometries and space-
time adaptive mesh refinements [25]. While the formulation of the proposed
approach readily generalizes to Tresca or Coulomb friction [7, 43, 20, 29],
additional challenges are expected for the nonlinear solver in this case.

## Appendix A

This appendix introduces space–time anisotropic Sobolev spaces as a convenient
setting for the analysis of time dependent boundary integral operators. In the
case of the wave equation, a detailed exposition may be found in [31], and we
refer to [4, 8, 9] for elastodynamics. When $\Gamma$ is an open screen or line
segment, so that $\partial\Gamma\neq\emptyset$, we first extend $\Gamma$ to a
closed, orientable Lipschitz manifold $\widetilde{\Gamma}$ of dimension $d-1$.
We recall the definition of Sobolev spaces of supported distributions on
$\Gamma$:

$\widetilde{H}^{s}(\Gamma)=\\{u\in H^{s}(\widetilde{\Gamma}):\mathrm{supp}\
u\subset{\overline{\Gamma}}\\}\ ,\quad\ s\in\mathbb{R}\ .$

The space ${H}^{s}(\Gamma)$ is defined as the quotient space
$H^{s}(\widetilde{\Gamma})/\widetilde{H}^{s}({\widetilde{\Gamma}\setminus\overline{\Gamma}})$.
We now define a family of Sobolev norms. Let $\alpha_{i}$, $i=1,\dots,p$, be a
partition of unity subordinate to a covering of $\widetilde{\Gamma}$ by open
sets $B_{i}\subset\mathbb{R}^{d-1}$. If for each $i$, $\varphi_{i}$ is a
diffeomorphism from $B_{i}$ to the unit cube in $\mathbb{R}^{d-1}$, we define
a family of Sobolev norms involving a parameter
$\omega\in\mathbb{C}\setminus\\{0\\}$:

$||u||_{s,\omega,{\widetilde{\Gamma}}}=\left(\sum_{i=1}^{p}\int_{{\mathbb{R}^{d-1}}}(|\omega|^{2}+|\boldsymbol{\xi}|^{2})^{s}|\mathcal{F}\left\\{(\alpha_{i}u)\circ\varphi_{i}^{-1}\right\\}(\boldsymbol{\xi})|^{2}d\boldsymbol{\xi}\right)^{\frac{1}{2}}\
.$

Here, $\mathcal{F}=\mathcal{F}_{\mathbb{x}\mapsto\boldsymbol{\xi}}$ denotes
the Fourier transform $\mathcal{F}\varphi(\boldsymbol{\xi})=\int
e^{-i\mathbb{x}\cdot\boldsymbol{\xi}}\varphi(\mathbb{x})\ d\mathbb{x}$.
Different parameters $\omega\in\mathbb{C}\setminus\\{0\\}$ lead to equivalent
norms on $H^{s}(\Gamma)$,
$\|u\|_{s,\omega,\Gamma}=\inf_{v\in\widetilde{H}^{s}(\widetilde{\Gamma}\setminus\overline{\Gamma})}\
\|u+v\|_{s,\omega,\widetilde{\Gamma}}$ and on $\widetilde{H}^{s}(\Gamma)$,
$\|u\|_{s,\omega,\Gamma,\ast}=\|e_{+}u\|_{s,\omega,\widetilde{\Gamma}}$.
$e_{+}$ here denotes the extension by $0$, which extends a distribution on
$\Gamma$ to a distribution on $\widetilde{\Gamma}$. When a specific, fixed
$\omega$ is considered, we write $H^{s}_{\omega}(\Gamma)$ for $H^{s}(\Gamma)$,
and $\widetilde{H}^{s}_{\omega}(\Gamma)$ for $\widetilde{H}^{s}(\Gamma)$. One
observes that $\|u\|_{s,\omega,\Gamma,\ast}\geq\|u\|_{s,\omega,\Gamma}$.

The space-time anisotropic Sobolev spaces relevant to this article can now be
defined as follows:

######  Definition 10.

For $\sigma>0$ and $r,s\in\mathbb{R}$ we set

$\displaystyle H^{r}_{\sigma}(\mathbb{R}^{+},{H}^{s}(\Gamma))$
$\displaystyle=\\{u\in\mathcal{D}^{{}^{\prime}}_{+}(H^{s}(\Gamma)):e^{-\sigma
t}u\in\mathcal{S}^{{}^{\prime}}_{+}(H^{s}(\Gamma))\textrm{ and
}||u||_{r,s,\Gamma}<\infty\\}\ ,$ $\displaystyle
H^{r}_{\sigma}(\mathbb{R}^{+},\widetilde{H}^{s}({\Gamma}))$
$\displaystyle=\\{u\in\mathcal{D}^{{}^{\prime}}_{+}(\widetilde{H}^{s}({\Gamma})):e^{-\sigma
t}u\in\mathcal{S}^{{}^{\prime}}_{+}(\widetilde{H}^{s}({\Gamma}))\textrm{ and
}||u||_{r,s,\Gamma,\ast}<\infty\\}\ .$ (40)

Here, $\mathcal{D}^{{}^{\prime}}_{+}({H}^{s}({\Gamma}))$ denotes the space of
all distributions on $\mathbb{R}$ with support in $[0,\infty)$, with values in
the real-valued subspace of the Hilbert space ${H}^{s}({\Gamma})$, and
$\mathcal{D}^{{}^{\prime}}_{+}(\widetilde{H}^{s}({\Gamma}))$ is defined in an
analogous way.
$\mathcal{S}^{{}^{\prime}}_{+}({H}^{s}({\Gamma}))\subset\mathcal{D}^{{}^{\prime}}_{+}({H}^{s}({\Gamma}))$
and
$\mathcal{S}^{{}^{\prime}}_{+}(\widetilde{H}^{s}({\Gamma}))\subset\mathcal{D}^{{}^{\prime}}_{+}(\widetilde{H}^{s}({\Gamma}))$
denote the subspaces of tempered distributions. The Sobolev spaces are
equipped with the norms

$\displaystyle\|u\|_{r,s,\sigma}:=\|u\|_{r,s,\Gamma,\sigma}$
$\displaystyle=\left(\int_{-\infty+i\sigma}^{+\infty+i\sigma}|\omega|^{2r}\
\|\hat{u}(\omega)\|^{2}_{s,\omega,\Gamma}\ d\omega\right)^{\frac{1}{2}}\ ,$
$\displaystyle\|u\|_{r,s,\sigma,\ast}:=\|u\|_{r,s,\Gamma,\sigma,\ast}$
$\displaystyle=\left(\int_{-\infty+i\sigma}^{+\infty+i\sigma}|\omega|^{2r}\
\|\hat{u}(\omega)\|^{2}_{s,\omega,\Gamma,\ast}\
d\omega\right)^{\frac{1}{2}}\,.$ (41)

They are Hilbert spaces. For $r=s=0$ they correspond to the weighted
$L^{2}$-space with the scalar product $\langle
u,v\rangle_{{\sigma,\Gamma,\mathbb{R}^{+}}}=\int_{0}^{\infty}e^{-2\sigma
t}\int_{\Gamma}u\,v\,d\Gamma_{\mathbb{x}}\ dt$ (see (13) in Subsection 2.2).

The boundary integral operators for the elastodynamic problem obey the
following mapping properties between these spaces:

######  Theorem 11.

Let $\sigma>0$. Then following operators are continuous for
$r\in\operatorname{\mathbb{R}}$:

$\displaystyle\mathcal{V}:{H}^{r+1}_{\sigma}(\operatorname{\mathbb{R}}^{+},\tilde{H}^{-\frac{1}{2}}(\Gamma))^{d}\to{H}^{r}_{\sigma}(\operatorname{\mathbb{R}}^{+},{H}^{\frac{1}{2}}(\Gamma))^{d}\
,$
$\displaystyle\mathcal{K}^{\ast}:{H}^{r+1}_{\sigma}(\operatorname{\mathbb{R}}^{+},\tilde{H}^{-\frac{1}{2}}(\Gamma))^{d}\to{H}^{r}_{\sigma}(\operatorname{\mathbb{R}}^{+},{H}^{-\frac{1}{2}}(\Gamma))^{d}\
,$
$\displaystyle\mathcal{K}:{H}^{r+1}_{\sigma}(\operatorname{\mathbb{R}}^{+},\tilde{H}^{\frac{1}{2}}(\Gamma))^{d}\to{H}^{r}_{\sigma}(\operatorname{\mathbb{R}}^{+},{H}^{\frac{1}{2}}(\Gamma))^{d}\
,$
$\displaystyle\mathcal{W}:{H}^{r+1}_{\sigma}(\operatorname{\mathbb{R}}^{+},\tilde{H}^{\frac{1}{2}}(\Gamma))^{d}\to{H}^{r}_{\sigma}(\operatorname{\mathbb{R}}^{+},{H}^{-\frac{1}{2}}(\Gamma))^{d}\
.$

Acknowledgements: This work has been partially supported by the University of
Parma with the project Fil2020 - Action A1 “Time-domain Energetic BEM for
elastodynamic problems, with advanced applications”. Early stages of this
research were supported through the “Oberwolfach Research Fellows” program in
2020.

## References

  * [1] T. Abboud, P. Joly, J. Rodriguez, I. Terrasse, _Coupling discontinuous Galerkin methods and retarded potentials for transient wave propagation on unbounded domains_ , J. Comp. Phys. 230 (2011), 5877–5907.
  * [2] A. Aimi, G. Di Credico, M. Diligenti, C. Guardasoni, _Highly accurate quadrature schemes for singular integrals in energetic BEM applied to elastodynamics_ , Journal of Computational and Applied Mathematics 410 (2022), 114186.
  * [3] A. Aimi, M. Diligenti, A. Frangi, C. Guardasoni, _Energetic BEM-FEM coupling for wave propagation in 3D multidomains_ , Internat. J. Numer. Methods Engrg. 97 (2014), 377-394.
  * [4] A. Aimi, G. Di Credico, H. Gimperlein, E. P. Stephan, _Higher-order time domain boundary elements for elastodynamics – graded meshes and hp-versions_ , Numerische Mathematik (2023), https://doi.org/10.1007/s00211-023-01355-x
  * [5] A. Bamberger, T. Ha Duong, _Formulation variationnelle espace-temps pour le calcul par potentiel retarde d’une onde acoustique_ , Math. Meth. Appl. Sci. 8 (1986), 405-435 and 598-608.
  * [6] L. Banz, H. Gimperlein, Z. Nezhi, E. P. Stephan, _Time domain BEM for sound radiation of tires_ , Computational Mechanics 58 (2016), 45-57.
  * [7] L. Banz, H. Gimperlein, A. Issaoui, E. P. Stephan, _Stabilized mixed hp-BEM for frictional contact problems in linear elasticity_ , Numer. Math. 135 (2017), 217-263.
  * [8] E. Bécache, _A variational boundary integral equation method for an elastodynamic antiplane crack_ , International Journal for Numerical Methods in Engineering 36 (1993), 969-984.
  * [9] E. Bécache, T. Ha Duong, _A space-time variational formulation for the boundary integral equation in a 2d elastic crack problem_ , ESAIM: Mathematical Modelling and Numerical Analysis 28 (1994), 141-176.
  * [10] E. Burman, M. Fernandez, S. Frei, _A Nitsche-based formulation for fluid-structure interactions with contact_ , ESAIM Mathematical Modelling and Numerical Analysis 54 (2020), 531-564.
  * [11] F. Chouly, P. Hild, _A Nitsche-based method for unilateral contact problems: numerical analysis_ , SIAM J. Numer. Anal. 51 (2013), 1295-1307.
  * [12] F. Chouly, P. Hild, Y. Renard, _A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments_ , ESAIM: Math. Model. Numer. Anal. 49 (2015), 503-528.
  * [13] M. Cocou, _Existence of solutions of a dynamic Signorini’s problem with nonlocal friction in viscoelasticity_ , Z. Angew. Math. Phys. 53 (2002), 1099-1109.
  * [14] R. G. Cooper, Jr. III, _Two variational inequality problems for the wave equation in a half-space_ , J. Math. Analysis Appl. 231 (1999), 434-460.
  * [15] G. Di Credico, _Energetic Boundary Element Method for 2D Elastodynamics Problems in Time Domain_ , Ph.D. thesis, University of Modena and Reggio Emilia (2022), available at http://hdl.handle.net/11380/1265215.
  * [16] D. Doyen, A. Ern, _Analysis of the modified mass method for the dynamic Signorini problem with Coulomb friction_ , SIAM J. Numer. Anal. 49 (2011), 2039-2056.
  * [17] C. Eck, J. Jarusek, M. Krbec, _Unilateral contact problems. Variational methods and existence theorems_ , Pure and Applied Mathematics 270. Chapman & Hall/CRC, Boca Raton, 2005.
  * [18] A. C. Eringen, E. S. Suhubi, _Elastodynamics_ , Academic Press, New York, (1975).
  * [19] S. Falletta, G. Monegato, L. Scuderi _Two boundary integral equation methods for linear elastodynamics problems on unbounded domains_ , Comput. Math. Appl 78 (2019), 3841–3861.
  * [20] H. Gimperlein, M. Maischak, E. Schrohe, E. P. Stephan, _Adaptive FE–BE coupling for strongly nonlinear transmission problems with Coulomb friction_ , Numer. Math. 117 (2011), 307–332.
  * [21] H. Gimperlein, M. Maischak, E. P. Stephan, _Adaptive time domain boundary element methods and engineering applications_ , J. Integral Equations Appl. 29 (2017), 75-105.
  * [22] H. Gimperlein, F. Meyer, C. Özdemir, E. P. Stephan, _Time domain boundary elements for dynamic contact problems_ , Comput. Methods Appl. Mech. Engrg. 333 (2018), 147-175.
  * [23] H. Gimperlein, Z. Nezhi, E. P. Stephan, _A priori error estimates for a time-dependent boundary element method for the acoustic wave equation in a half-space_ , Mathematical Methods in the Applied Sciences 40 (2017), 448-462.
  * [24] H. Gimperlein, C. Özdemir, D. Stark, E. P. Stephan, _hp-version time domain boundary elements for the wave equation on quasi-uniform meshes_ , Comput. Methods Appl. Mech. Engrg. 356 (2019), 145-174.
  * [25] H. Gimperlein, C. Özdemir, D. Stark, E. P. Stephan, _A residual a posteriori estimate for the time-domain boundary element method_ , Numerische Mathematik 146 (2020), 239-280.
  * [26] H. Gimperlein, C. Özdemir, E. P. Stephan, _Time domain boundary element methods for the Neumann problem: error estimates and acoustic problems_ , J. Comput. Math. 36 (2018), 70-89.
  * [27] H. Gimperlein, J. Stocek, _Space-time adaptive finite elements for nonlocal parabolic variational inequalities_ , Comput. Methods Appl. Mech. Engrg. 352 (2019), 137-171.
  * [28] T. Gustafsson, R. Stenberg, J. Videman, _Nitsche’s method for unilateral contact problems_ , preprint.
  * [29] J. Gwinner, N. Ovcharova, _Coupling of finite element and boundary element methods with regularization for a nonlinear interface problem with nonmonotone set-valued transmission conditions_ , Comput. Math. Appl. 134 (2023), 45-54.
  * [30] J. Gwinner, E. P. Stephan, _Advanced Boundary Element Methods – Treatment of Boundary Value, Transmission and Contact Problems_ , Springer Series in Computational Mathematics, vol. 52, 2018.
  * [31] T. Ha Duong, _On retarded potential boundary integral equations and their discretizations_ , in: Topics in computational wave propagation, pp. 301-336, Lect. Notes Comput. Sci. Eng., 31, Springer, Berlin, 2003.
  * [32] P. Hauret, P. Le Tallec, _Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact_ , Comp. Methods Appl. Mech. Engrg. 195 (2006), 4890-4916.
  * [33] P. Joly, J. Rodriguez, _Mathematical aspects of variational boundary integral equations for time dependent wave propagation_ , J. Integral Equations Appl. 29 (2017), 137–187.
  * [34] B. Kager, M. Schanz, _Fast and data sparse time domain BEM for elastodynamics_ , Eng. Anal. Bound. Elem. 50 (2015), 212–223.
  * [35] H. B. Khenous, P. Laborde, Y. Renard, _On the discretization of contact problems in elastodynamics_ , Lect. Notes Appl. Comput. Mech. 27 (2006), 31-38.
  * [36] G. Lebeau, M. Schatzman, _A wave problem in a half-space with a unilateral constraint at the boundary_ , J. Differential Equations 53 (1984), 309-361.
  * [37] P. Le Tallec, P. Hauret, B. Wohlmuth, C. Hager, _Solving dynamic contact problems with local refinement in space and time_ , Comp. Methods Appl. Mech. Engrg. 201-204 (2012), 25-41.
  * [38] A. I. Luré, _Theory of elasticity. Foundations of engineering mechanics_ , Springer, Berlin, 2005.
  * [39] R. Mlika, Y. Renard, F. Chouly, _An unbiased Nitsche’s formulation of large deformation frictional contact and self-contact_ , Comput. Methods Appl. Mech. Engrg. 325 (2017), 265-288.
  * [40] F.-J. Sayas, _Retarded Potentials and Time Domain Boundary Integral Equations: A Road Map_ , Springer Series in Computational Mathematics 50 (2016).
  * [41] M. Schanz, W. Ye, J. Xiao , _Comparison of the convolution quadrature method and enhanced inverse FFT with application in elastodynamic boundary element method_ , Comp. Mech. 57 (2016), 523–536.
  * [42] O. Steinbach, C. Urzúa-Torres, _A New Approach to Space-Time Boundary Integral Equations for the Wave Equation_ , SIAM Journal on Mathematical Analysis 54 (2022), 1370-1392.
  * [43] E. P. Stephan, _Coupling of boundary element methods and finite element methods_ , In: _Encyclopedia of Computational Mechanics, Fundamentals_ , E. Stein, R. de Borst, T. J. R. Hughes (eds.), vol. I, Wiley, Amsterdam, 2004, 375-412.
  * [44] P. Wriggers, _Computational Contact Mechanics_ , Wiley, New York, 2002.