Heiko Gimperlein
email: h.gimperlein@hw.ac.uk.

Fabian Meyer

Özdemir Ceyhun

Ernst P Stephan

Ceyhun Özdemir

Time domain boundary elements for dynamic contact problems

Keywords: boundary element method, variational inequality, mixed method, a priori error estimates, wave equation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Contact problems play an important role in numerous applications in mechanics, from fracture dynamics and crash tests to rolling car tires [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. As the contact takes place at the interface of two materials, for time-independent problems boundary elements and coupled finite / boundary elements provide an efficient and much-studied tool for numerical simulations [START_REF] Gwinner | Advanced Boundary Element Methods -Treatment of Boundary Value, Transmission and Contact Problems[END_REF][START_REF] Stephan | Coupling of boundary element methods and finite element methods[END_REF]. The analysis of such problems is well-understood in the context of elliptic variational inequalities.

While contact for time-dependent problems is of clear practical relevance, neither its analysis nor rigorous boundary element methods have been much explored. There is an extensive computational literature, including [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF][START_REF] Hauret | Mixed interpretation and extensions of the equivalent mass matrix approach for elastodynamics with contact[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF][START_REF] Tallec | Solving dynamic contact problems with local refinement in space and time[END_REF], but analytically even the existence of solutions to these free boundary problem is only known for flat contact area [START_REF] Cooper | Two variational inequality problems for the wave equation in a halfspace[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]. Some rigorous results have recently been obtained for Nitsche stabilized finite elements [START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF].

In this work we propose a time domain boundary element method for a dynamic contact problem in the case of the scalar wave equation, as a model problem for elasticity. We provide a priori error estimates for our numerical scheme in the case of a flat contact area, and our numerical experiments indicate the convergence and efficiency also for curved contact geometries. Motivating references from the time-independent setting include [START_REF] Babuska | On the mixed finite element method with Lagrange multipliers[END_REF][START_REF] Eck | A symmetric boundary element method for contact problems with friction[END_REF].

For the precise statement of the problem, consider a Lipschitz domain Ω ⊂ R n with boundary Γ = ∂Ω. Let G be a bounded Lipschitz subset of Γ. We consider a unilateral contact problem for the wave equation for the displacement w : R × Ω → R. It corresponds to a simplified model for a crack in G between Ω and a non-penetrable material in R n \ Ω. The contact conditions for non-penetration are described in terms of the traction -µ ∂w ∂ν G and prescribed forces h:

{ w| G ≥ 0 , -µ ∂w ∂ν G ≥ h , w| G > 0 =⇒ -µ ∂w ∂ν G = h .
The full system of equations for the contact problem is given by:

                 ∂ 2 w ∂t 2 = c 2 s ∆w , for (t, x) ∈ R × Ω , w = 0 , on Γ \ G ,
w ≥ 0 , -µ ∂w ∂ν ≥ h , on G , ( -µ ∂w ∂ν -h) w = 0, on G , w = 0, for (t, x) ∈ (-∞, 0) × Ω .

(

Like for time-independent contact, a formulation as a nonlinear problem on the contact area G leads to efficient numerical approximations.

In this article we formulate (1) as a variational inequality on G in terms of the Dirichlet-to-Neumann operator for the wave equation. Similar to time-independent problems, the Dirichletto-Neumann operator may be computed in terms of boundary integral operators. In order to improve the conservation of energy in the numerical scheme, we replace the variational inequality by an equivalent mixed system, which we discretize with a time domain Galerkin boundary element method. The resulting discretized nonlinear inequality in space-time simultaneously approximates the displacement w and the contact forces -µ ∂w ∂ν on G. It is solved with a Uzawa algorithm, either as a time-stepping scheme or in space-time.

The resulting boundary element method is analyzed in the case of a flat contact area, a situation where the existence of solutions to the contact problem (1) is known. We obtain a priori estimates for the numerical error, both for the variational inequality and a mixed formulation. For the mixed problem, a key part of the proof is an inf-sup condition for the space-time discretization [START_REF] Babuska | On the mixed finite element method with Lagrange multipliers[END_REF]. We also demonstrate the convergence of the Uzawa algorithms. Numerical experiments for the mixed formulation confirm the theoretical results and indicate the efficiency of our approach, beyond flat geometries.

In addition to the contact problem [START_REF] Abboud | Coupling discontinuous Galerkin methods and retarded potentials for transient wave propagation on unbounded domains[END_REF], as a simpler test case we also consider a punch problem, which models a rigid body indenting an elastic half-space. The contact boundary conditions for the wave equation are given by:

{ -µ ∂w ∂ν G ≥ 0 , w G ≥ h , -µ ∂w ∂ν G > 0 =⇒ w G = h .
For the half-space they lead to a variational inequality for the single layer operator V , instead of the Dirichlet-to-Neumann operator, and we obtain similar theoretical and numerical results in this case.

Our approach also relates to the recent interest in coupled and nonlinear interface problems for wave propagation, solved by time domain boundary element methods. In particular, we refer to the fundamental articles [START_REF] Abboud | Coupling discontinuous Galerkin methods and retarded potentials for transient wave propagation on unbounded domains[END_REF][START_REF] Banjai | Stable numerical coupling of exterior and interior problems for the wave equation[END_REF] for the coupling of FEM and BEM, as well as [START_REF] Aimi | Energetic BEM-FEM coupling for wave propagation in 3D multidomains[END_REF] for an energetic Galerkin formulation of the coupling. Reference [START_REF] Banjai | Convolution quadrature for the wave equation with a nonlinear impedance boundary condition[END_REF] considers a nonlinear boundary problem. A first analysis of the time domain Dirichlet-to-Neumann operator goes back to [START_REF] Banjai | Time-domain Dirichlet-to-Neumann map and its discretization[END_REF].

The current work provides a first step towards efficient boundary elements for dynamic contact. For both stationary and dynamic contact, the relevance of adaptive methods to approximate the non-smooth solutions is well-known [START_REF] Gwinner | Advanced Boundary Element Methods -Treatment of Boundary Value, Transmission and Contact Problems[END_REF][START_REF] Tallec | Solving dynamic contact problems with local refinement in space and time[END_REF][START_REF] Stephan | Coupling of boundary element methods and finite element methods[END_REF]. Recent advances in the a posteriori error analysis and resulting adaptive mesh refinement procedures for time domain boundary elements [START_REF] Gimperlein | A residual a posteriori estimate for the time-domain boundary element method[END_REF][START_REF] Gimperlein | Adaptive time domain boundary element methods and engineering applications[END_REF] will therefore be of interest for the dynamic contact considered here, with a particular view towards tire dynamics [START_REF] Banz | Time domain BEM for sound radiation of tires[END_REF].

The article is organized as follows: Section 2 recalls the boundary integral operators associated to the wave equation as well as their mapping properties between suitable space-time anisotropic Sobolev spaces. Section 3 reduces the contact problem to a variational inequality on the contact boundary. Section 4 discusses the existence and uniqueness of solutions to the contact problem and a related, simpler Dirichlet-to-Neumann equation in a half space. Section 5 describes the discretization and proves a priori error estimates for the contact problem as well as for the Dirichlet-to-Neumann equation. The a priori error analysis for the mixed formulation is then presented in Section 6. A simpler contact problem, which involves the single layer operator V only, is analyzed in Section 7. Section 8 derives a time stepping scheme for the Dirichlet-to-Neumann equation, which forms the basis for both a time stepping and a space-time Uzawa algorithm for the nonlinear contact problem. Section 9 presents numerical experiments based on the mixed formulation.

Boundary integral operators and Sobolev spaces

We introduce the single layer potential in time domain as

Sφ(t, x) = ∫ R + ×Γ G(t -τ, x, y) φ(τ, y) dτ ds y ,
where (t, x) ∈ R + × Ω and G is a fundamental solution to the wave equation. Specifically in 3 dimensions, it is given by

Sφ(t, x) = 1 4π ∫ Γ φ(t -|x -y|, y) |x -y| ds y .
We similarly define the double-layer potential as

Dφ(t, x) = ∫ R + ×Γ ∂G ∂n y (t -τ, x, y) φ(τ, y) dτ ds y .
For the Dirichlet-to-Neumann operator, we require the single-layer operator V , its normal derivative K ′ , the double-layer operator K and hypersingular operator W for x ∈ Γ, t > 0:

V φ(t, x) = ∫ R + ×Γ G(t -τ, x, y) φ(τ, y) dτ ds y , Kφ(t, x) = ∫ R + ×Γ ∂G ∂n y (t -τ, x, y) φ(τ, y) dτ ds y , K ′ φ(t, x) = ∫ R + ×Γ ∂G ∂n x (t -τ, x, y) φ(τ, y) dτ ds y , W φ(t, x) = ∫ R + ×Γ ∂ 2 G ∂n x ∂n y (t -τ, x, y) φ(τ, y) dτ ds y . Remark 1. When Ω = R n + , the normal derivative of G vanishes on Γ = ∂R n + . Therefore, Kφ = K ′ φ = 0 in this case.
The boundary integral operators are considered between space-time anisotropic Sobolev spaces H s σ (R + , H r (Γ)), see [START_REF] Gimperlein | A priori error estimates for a time-dependent boundary element method for the acoustic wave equation in a half-space[END_REF] or [START_REF] Ha-Duong | On retarded potential boundary integral equations and their discretisations, Topics in computational wave propagation[END_REF]. To define them, if ∂Γ ̸ = ∅, first extend Γ to a closed, orientable Lipschitz manifold Γ.

On Γ one defines the usual Sobolev spaces of supported distributions:

H r (Γ) = {u ∈ H r ( Γ) : supp u ⊂ Γ} , r ∈ R . Furthermore, H r (Γ) is the quotient space H r ( Γ)/ H r ( Γ \ Γ).
To write down an explicit family of Sobolev norms, introduce a partition of unity α i subordinate to a covering of Γ by open sets B i . For diffeomorphisms φ i mapping each B i into the unit cube ⊂ R n , a family of Sobolev norms is induced from R n :

||u|| r,ω, Γ = ( p ∑ i=1 ∫ R n (|ω| 2 + |ξ| 2 ) r |F { (α i u) • φ -1 i } (ξ)| 2 dξ ) 1 2 .
The norms for different ω ∈ C \ {0} are equivalent and F denotes the Fourier transform. They induce norms on H r (Γ), ||u|| r,ω,Γ = inf v∈ H r ( Γ\Γ) ||u+v|| r,ω, Γ and on H r (Γ), ||u|| r,ω,Γ, * = ||e + u|| r,ω, Γ .

e + extends the distribution u by 0 from Γ to Γ. It is stronger than ||u|| r,ω,Γ whenever r ∈ 1 2 + Z. We now define a class of space-time anisotropic Sobolev spaces: Definition 2. For σ > 0 and s, r ∈ R define .

H s σ (R + , H r (Γ)) = {u ∈ D ′ + (H r (Γ)) : e -σt u ∈ S ′ + (H r (Γ)) and ||u|| s,r,σ,Γ < ∞} , H s σ (R + , H r (Γ)) = {u ∈ D ′ + ( H r (Γ)) : e -σt u ∈ S ′ + ( H r (Γ))
For |r| ≤ 1 the spaces are independent of the choice of α i and φ i . We further introduce the set of nonnegative distributions H r σ (R + , Hs (G)) + . The boundary integral operators obey the following mapping properties between these spaces: Theorem 3 ([21]). The following operators are continuous for r ∈ R:

V : H r+1 σ (R + , H-1 2 (Γ)) → H r σ (R + , H 1 2 (Γ)) , K ′ : H r+1 σ (R + , H-1 2 (Γ)) → H r σ (R + , H -1 2 (Γ)) , K : H r+1 σ (R + , H 1 2 (Γ)) → H r σ (R + , H 1 2 (Γ)) , W : H r+1 σ (R + , H 1 2 (Γ))) → H r σ (R + , H -1 2 (Γ)) .
Furthermore, as noted by Bamberger and Ha Duong [START_REF] Bamberger | Formulation variationnelle espace-temps pour le calcul par potentiel retarde d'une onde acoustique[END_REF], V and W are coercive in a weak norm when composed with a time derivative:

∥ϕ∥ 2 0,-1 2 ,σ, * ⟨V ϕ, ∂ t ϕ⟩, resp. ∥ψ∥ 2 0, 1 2 ,σ, * ⟨W ψ, ∂ t ψ⟩.
This is not sufficient for the analysis of variational inequalities, which no longer hold after differentiation in time.

When Ω = R n + , Fourier methods yield improved estimates for V and W , see also Section 4:

Theorem 4 ([25], pp. 503-506). The following operators are continuous for r, s ∈ R:

V : H r+ 1 2 σ (R + , Hs (Γ)) → H r σ (R + , H s+1 (Γ)) , W : H r σ (R + , Hs (Γ)) → H r σ (R + , H s-1 (Γ)) .
Theorems 3 and 4 imply the corresponding mapping properties for the composition with the restriction p Q to Q = R×G. For example, from Theorem 3 we obtain p Q V :

H r+1 σ (R + , H-1 2 (Γ)) → H r σ (R + , H 1 2 (G)) and p Q V : H r+1 σ (R + , H-1 2 (G)) → H r+1 σ (R + , H-1 2 (Γ)) → H r σ (R + , H 1 2 (G)
). Also, the crucial coercivity without composition with ∂ t can be shown in this flat geometry. In particular, ∥ϕ∥ 2

-1 2 ,-1 2 ,σ, * σ ⟨p Q V ϕ, ϕ⟩.

Boundary integral formulation

For some σ > 0, we write w σ = e -σt w, resp. h σ = e -σt h. Using appropriate units, we may also assume c s = 1. Multiplying ( 1) by e -σt , we then obtain

                 ( ∂ ∂t + σ ) 2 w σ = ∆w σ , for (t, x) ∈ R × Ω , w σ = 0 , on Γ \ G , w σ ≥ 0 , -µ ∂wσ ∂ν G ≥ h σ , on G , ( -µ ∂wσ ∂ν G -h σ ) w σ | G = 0, on G , w σ = 0, for (t, x) ∈ (-∞, 0) × Ω .
(

) 2
The Dirichlet-to-Neumann operator, defined as

S σ w σ | Γ := -µ ∂w σ ∂ν Γ , ( 3
)
allows to reduce the contact problem (1) to an equivalent inequality in terms of the trace u σ := w σ | Γ on the boundary: Find

u σ with supp u σ ⊂ Q 0 = R + × G such that u σ ≥ 0 , S σ u σ ≥ h σ , (S σ u σ -h σ ) u σ = 0 on Q = R × G . ( 4
)
More precisely, in terms of the restriction p Q to Q, we obtain the following weak formulation as a variational inequality for suitably smooth h σ : Find

u σ ∈ H 1 2 σ (R + , H 1 2 (
G)) such that:

u σ ≥ 0 and ∀v ∈ H 1 2 σ (R + , H 1 2 (G)) with v ≥ 0 : ⟨p Q S σ u σ , v -u σ ⟩ ≥ ⟨h σ , v -u σ ⟩ . ( 5
)

Existence and uniqueness for contact problems

We now restrict ourselves to the half space Ω = R n + , where sharp continuity and coercivity results allow to show the wellposedness of the contact problem.

We illustrate the use of Fourier methods in this geometry. As shown in [START_REF] Bennish | Mixed initial-boundary value problems for hyperbolic equations with constant coefficients[END_REF][START_REF] Ha-Duong | On the transient acoustic scattering by a flat object[END_REF] the Dirichlet-Neumann operator S σ is a pseudo-differential operator on

R t × ∂Ω x = R × R n-1 × {0}: S σ u σ (t, x) = 1 (2π) n ∫ R n e itξ 0 +ixξ ′ (-iµΓ(ξ 0 + iσ, ξ ′ ))û σ (ξ 0 , ξ ′ )dξ 0 dξ ′ , with symbol Γ(ξ 0 + iσ, ξ ′ ) = √ (ξ 0 + iσ) 2 -|ξ ′ | 2 .
From this explicit formula, one notices the estimates ( [START_REF] Ha-Duong | On the transient acoustic scattering by a flat object[END_REF], p. 499):

σ 1 2 |ξ 0 +iσ| 1 2 (|ξ 0 + iσ| 2 + |ξ ′ | 2 ) 1 2 ≤ |Γ(ξ 0 + iσ, ξ ′ )| ≤ (|ξ 0 + iσ) 2 + |ξ ′ | 2 ) 1 2 , Im Γ(ξ 0 + iσ, ξ ′ ) ≥ σ |ξ 0 +iσ| (|ξ 0 + iσ| 2 + |ξ ′ | 2 ) 1 2 .
They translate into the following coercivity and mapping properties:

Theorem 5. p Q S σ : H s σ (R + , H 1 2 (G)) → H s σ (R + , H -1 2 (G)) continuously and ∥ϕ∥ 2 -1 2 , 1 2 ,σ, * σ ⟨p Q S σ ϕ, ϕ⟩ ≤ C∥ϕ∥ 2 0, 1 2 ,σ, * .
These estimates are a basic ingredient to prove the well-posedness of the contact problem:

Theorem 6 ([13], p. 450). Let h ∈ H 3 2 σ (R + , H -1 2 (G)).
Then there exists a unique solution

u σ ∈ H 1 2 σ (R + , H 1 2 (G)) + of (5).
In terms of the original problem (2) we obtain:

Theorem 7 ([13], p. 451). Let h ∈ H 3 2 σ (R + , H -1 2 (G)). Then there exists a unique w(•, x n ) ∈ C(R + xn ; H 1 2 σ (R + , H 1 2 (R n-1 )) ∩ H 0 σ (R + , H 1 (R n )) satisfying (2).
We also note the (simpler) existence of solutions to the corresponding equality from [START_REF] Sako | A model for the Crack and Punch Problems in Elasticity[END_REF], p. 48:

Theorem 8. Let h ∈ H 3 2 σ (R + , H -1 2 (G)).
Then there exists a unique u σ ∈ H

1 2 σ (R + , H 1 2 (G)) which solves: ⟨p Q S σ u σ , v⟩ = ⟨h, v⟩ ∀v ∈ H -1 2 σ (R + , H 1 2 (G)).

Discretization and a priori error estimates

We assume that n = 2 or n = 3, and if Γ is not polygonal, we approximate it by a piecewise polygonal curve resp. surface. The approximation is again denoted by Γ. When n = 3, we may assume that the surface Γ is composed of N closed triangular faces

S i such that Γ = ∪ N i=1 S i . When n = 2, we assume Γ = ∪ N i=1 S i is composed of line segments S i .
The triangulation T S := {S 1 , . . . , S N } is assumed to be quasi-uniform and compatible with the area of contact G:

If S i ∩ G ̸ = ∅, int S i ⊂ G.
Associated to the triangulation, we obtain the space V q h (Γ) of piecewise polynomial functions of degree q. Due to the compatibility of the meshes we have V q h (G) ⊂ V q h (Γ). Moreover we define Ṽ q h (G) as the subspace of those functions in V q h (G), which vanish on ∂G for q ≥ 1. For the time discretization we consider a uniform decomposition of the time interval [0, ∞) into subintervals [t n-1 , t n ) with time step ∆t, such that t n = n∆t (n = 0, 1, . . . ). Associated to this mesh, the space V p t consists of piecewise polynomial functions of degree of p (continuous and vanishing at t = 0 if p ≥ 1).

Let T T = {[0, t 1 ), [t 1 , t 2 ), • • • , [t N -1 , T )} the
time mesh for a finite subinterval [0, T ). In spacetime we consider the tensor product of the approximation spaces, V q h and V p t , associated to the space-time mesh T S,T = T S × T T , and we write

V p,q t,h := V p t ⊗ V q h . ( 6
)
We analogously define

Ṽ p,q t,h := V p t ⊗ Ṽ q h (7)
We further define the subspace K + t,h ⊂ V p,q t,h as the subspace of nonnegative piecewise polynomials.

For u ∆t,h ∈ V p,q t,h we thus may write

u ∆t,h (t, x) = Nt ∑ i=0 Ns ∑ j=0 c i j β i ∆t (t)ξ j h (x) .
in terms of the basis functions β i ∆t in time and ξ j h in space.

We use the following notation for piecewise linear or constant functions:

• γ n ∆t (t) for the basis of piecewise constant functions in time,

• β n ∆t (t) for the basis of piecewise linear functions in time,

• ψ i h (x) for the basis of piecewise constant functions in space,

• ξ i h (x) for the basis of piecewise linear functions in space.

We recall the formulation as a continuous variational inequality:

Find u ∈ H 1 2 σ (R + , H 1 2 (G)) + such that ⟨p Q S σ u σ , v -u σ ⟩ ≥ ⟨h, v -u σ ⟩ (8) holds for all v ∈ H 1 2 σ (R + , H 1 2 (G)) + . The discretized variational inequality reads as follows. Find u ∆t,h ∈ K+ t,h such that ⟨p Q S σ u ∆t,h , v ∆t,h -u ∆t,h ⟩ ≥ ⟨h, v ∆t,h -u ∆t,h ⟩ (9)
holds for all v ∆t,h ∈ K+ t,h .

Note that we assume that S σ is computed exactly. We refer to [START_REF] Banjai | Time-domain Dirichlet-to-Neumann map and its discretization[END_REF] for the challenges of analyzing its discretization. While σ > 0 is required for the theoretical analysis, practical computations use σ = 0 [START_REF] Bamberger | Formulation variationnelle espace-temps pour le calcul par potentiel retarde d'une onde acoustique[END_REF][START_REF] Gimperlein | Adaptive time domain boundary element methods and engineering applications[END_REF].

Using a conforming ansatz space, we are able to derive an a priori estimate for the variational inequality. It is the hyperbolic counterpart of the elliptic estimate proved by Falk [START_REF] Falk | Error estimates for the approximation of a class of variational inequalities[END_REF].

Theorem 9. Let h ∈ H 3 2 σ (R + , H -1 2 (G)) and let u ∈ H 1 2 σ (R + , H 1 2 (G)) + , u ∆t,h ∈ K+ t,h ⊂ H 1 2 σ (R + , H 1 2 (
G)) + be the solutions of [START_REF] Banz | Stabilized mixed hp-BEM for frictional contact problems in linear elasticity[END_REF], resp. [START_REF] Banz | Time domain BEM for sound radiation of tires[END_REF]. Then the following estimate holds:

∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ σ inf ϕ ∆t,h ∈ K+ t,h (∥h -p Q S σ u∥ 1 2 ,-1 2 ,σ ∥u -ϕ ∆t,h ∥ -1 2 , 1 2 ,σ,⋆ + ∥u -ϕ ∆t,h ∥ 2 1 2 , 1 2 ,σ,⋆ ). ( 10
)
Proof. Rewriting ( 8) and ( 9), we note that

⟨p Q S σ u, u⟩ ≤ ⟨h, u -ϕ⟩ + ⟨p Q S σ u, ϕ⟩ (11)
and

⟨p Q S σ u ∆t,h , u ∆t,h ⟩ ≤ ⟨h, u ∆t,h -ϕ ∆t,h ⟩ + ⟨p Q S σ u ∆t,h , ϕ ∆t,h ⟩. ( 12
)
Using the coercivity in the

∥ • ∥ -1 2 , 1 2 ,σ,⋆ -norm (Theorem 5), we obtain ∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ σ ⟨p Q S σ (u -u ∆t,h ), u -u ∆t,h ⟩ ≤ ⟨h, u -ϕ⟩ + ⟨h, u ∆t,h -ϕ ∆t,h ⟩ + ⟨p Q S σ u ∆t,h , ϕ ∆t,h ⟩ + ⟨p Q S σ u, ϕ⟩ -⟨p Q S σ u, u ∆t,h ⟩ -⟨p Q S σ u ∆t,h , u⟩ = ⟨h, u -ϕ⟩ + ⟨h, u ∆t,h -ϕ ∆t,h ⟩ + ⟨p Q S σ u, ϕ -u ∆t,h ⟩ + ⟨p Q S σ u ∆t,h , ϕ ∆t,h -u⟩.
We rewrite

⟨p Q S σ u ∆t,h , ϕ ∆t,h -u⟩ = ⟨p Q S σ u -u ∆t,h , u -ϕ ∆t,h ⟩ -⟨p Q S σ u, u -ϕ ∆t,h ⟩ , so that ∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ σ ⟨h -p Q S σ u, u -ϕ ∆t,h ⟩ + ⟨h -p Q S σ u ∆t,h , u ∆t,h -ϕ⟩ + ⟨p Q S σ u -u ∆t,h , u -ϕ ∆t,h ⟩.
Because of the conforming discretization, we may choose ϕ = u ∆t,h and conclude

∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ σ ⟨h -p Q S σ u, u -ϕ ∆t,h ⟩ + ⟨p Q S σ u -u ∆t,h , u -ϕ ∆t,h ⟩. ( 13
)
We estimate both terms by duality:

∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ σ ∥h -p Q S σ u∥ 1 2 ,-1 2 ,σ ∥u -ϕ ∆t,h ∥ -1 2 , 1 2 ,σ,⋆ + ∥p Q S σ (u -u ∆t,h )∥ -1 2 ,-1 2 ,σ ∥u -ϕ ∆t,h ∥ 1 2 , 1
2 ,σ,⋆ . From the mapping properties of p Q S σ , Theorem 5, one then sees that

∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ σ ∥h -p Q S σ u∥ 1 2 ,-1 2 ,σ ∥u -ϕ ∆t,h ∥ -1 2 , 1 2 ,σ,⋆ + ∥u -u ∆t,h ∥ -1 2 , 1 2 ,σ,⋆ ∥u -ϕ ∆t,h ∥ 1 2 , 1
2 ,σ,⋆ . We conclude with the help of Young's inequality

∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ σ ∥h -p Q S σ u∥ 1 2 ,-1 2 ,σ ∥u -ϕ ∆t,h ∥ -1 2 , 1 2 ,σ,⋆ + ∥u -ϕ ∆t,h ∥ 2 1 2 , 1
2 ,σ,⋆ . Taking the infimum over all ϕ ∆t,h yields the assertion.

The theorem implies explicit convergence rates, using approximation properties in anisotropic Sobolev spaces as stated e.g. in [START_REF] Gläfke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF]:

Corollary 10. Let u ∈ H 1 2 +ϵ σ (R + , H 1 2 +ϵ (G))
for some ϵ > 0. There exists a constant C > 0 (independent of ∆t and h) such that

∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ ≤ C(h ϵ + (∆t) 1 2 +ϵ )∥u∥ 1 2 +ϵ, 1 2 +ϵ,σ + C(h 2ϵ + (∆t) 2ϵ )∥u∥ 2 1 2 +ϵ, 1 2 +ϵ,σ . Proof. We estimate the first term ∥u -ϕ ∆t,h ∥ -1 2 , 1
2 ,σ,⋆ on the right hand side in Theorem 9 by ∥u -ϕ ∆t,h ∥ 0, 1 2 ,σ,⋆ . Then we apply Proposition 3.56 in [START_REF] Gläfke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF]. We also consider the continuous and discrete variational equalities for the Dirichlet-to-Neumann operator.

Find u ∈ H 1 2 σ (R + , H 1 2 (G)) such that ⟨p Q S σ u σ , v⟩ = ⟨h, v⟩ (14)
holds for all v ∈ H

-1 2 σ (R + , H 1 2 (G)). Find u ∆t,h ∈ Ṽ p,q t,h such that ⟨p Q S σ u ∆t,h , ϕ ∆t,h ⟩ = ⟨h, ϕ ∆t,h ⟩ ( 15
)
holds for all ϕ ∆t,h ∈ Ṽ p,q t,h .

An a priori estimate for the variational equality is obtained from the previous arguments as a special case.

Theorem 11. Let u ∈ H 1 2 σ (R + , H 1 2 (G)) and u ∆t,h ∈ Ṽ p,q t,h ⊂ H 1 2 σ (R + , H 1 2 (G)
) be the solutions of (14)), resp. [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. We have the following a priori estimate:

∥u -u ∆t,h ∥ -1 2 , 1 2 ,σ,⋆ σ inf ϕ ∆t,h ∈ V p,q t,h ∥u -ϕ ∆t,h ∥ 1 2 , 1 2 ,σ,⋆ .
Proof. For the conforming test space Ṽ p,q t,h we have

⟨p Q S σ u, ϕ ∆t,h ⟩ = ⟨h, ϕ ∆t,h ⟩,
for all ϕ ∆t,h ∈ Ṽ p,q t,h . Therefore [START_REF] Cooper | Two variational inequality problems for the wave equation in a halfspace[END_REF] becomes

∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ σ ⟨h -p Q S σ u, u -ϕ ∆t,h ⟩ + ⟨p Q S σ (u -u ∆t,h ), u -ϕ ∆t,h ⟩ = ⟨p Q S σ (u -u ∆t,h ), u -ϕ ∆t,h ⟩ ≤ ∥p Q S σ (u -u ∆t,h )∥ -1 2 ,-1 2 ,σ ∥u -ϕ ∆t,h ∥ 1 2 , 1 2 ,σ,⋆ . The continuity of p Q S σ , ∥p Q S σ (u-u ∆t,h )∥ -1 2 ,-1 2 ,σ ∥u-u ∆t,h ∥ -1 2 , 1
2 ,σ, * , yields the assertion. Analogous to Corollary 10 we derive the following rate of convergence.

Corollary 12. Let u ∈ H 1 2 +ϵ σ (R + , H 1 2 +ϵ (G))
for some ϵ > 0. There exists a constant C > 0 (independent of ∆t and h) such that

∥u -u ∆t,h ∥ 2 -1 2 , 1 2 ,σ,⋆ ≤ C(h 2ϵ + (∆t) 2ϵ )∥u∥ 2 1 2 +ϵ, 1 2 +ϵ,σ .

Mixed formulation

To assure conservation of energy in the numerical approximation, it proves useful to impose the constraints on the displacement only indirectly. We therefore reformulate the variational inequality as an equivalent mixed system. The Lagrange multiplier λ = S σ u -h in this formulation provides a measure to which extent the variational inequality is not an equality; physically, it corresponds to the contact forces and indicates the contact area within the computational domain.

Theorem 13 (Mixed formulation). Let h ∈ H

3 2 σ (R + , H -1 2 (G)).
The variational inequality formulation ( 8) is equivalent to the following formulation:

Find (u, λ) ∈ H 1 2 σ (R + , H 1 2 (G)) × H 1 2 σ (R + , H -1 2 (G)) + such that { (a) ⟨S σ u, v⟩ -⟨λ, v⟩ = ⟨h, v⟩ (b) ⟨u, µ -λ⟩ ≥ 0, ( 16
)
for all (v, µ) ∈ H 1 2 σ (R + , H 1 2 (G)) × H 1 2 σ (R + , H -1 2 (G)) + .
Proof. We first note, that ( 8) is equivalent to the following problem. Find u ∈ H

1 2 σ (R + , H 1 2 (G)) + solving { (a) ⟨S σ u, u⟩ = ⟨h, u⟩ (b) ⟨S σ u, v⟩ ≥ ⟨h, v⟩, (17)
for all v ∈ H 1 2 σ (R + , H 1 2 (G)) + . Setting v = 2u, v = 0 in the variational inequality (8) we obtain ⟨S σ u, u⟩ ≥ ≤ ⟨h, u⟩ ⇒ ⟨S σ u, u⟩ = ⟨h, u⟩ (*)
If we add ( * ) to [START_REF] Banz | Stabilized mixed hp-BEM for frictional contact problems in linear elasticity[END_REF] we obtain the second line in [START_REF] Falk | Error estimates for the approximation of a class of variational inequalities[END_REF] To get ( 8) from [START_REF] Falk | Error estimates for the approximation of a class of variational inequalities[END_REF] we subtract (17) (a) from (17b).

We now show the equivalence of ( 17) and ( 16): (17) ⇒ (16): If we set λ = S σ u -h we have by (17b) :

⟨S σ u -h, v⟩ ≥ 0 for all v ∈ H 1 2 σ (R + , H 1 2 (G)) + and therefore λ ∈ H 1 2 σ (R + , H -1 2 (G)) + .
The first line in [START_REF] Eck | A symmetric boundary element method for contact problems with friction[END_REF] holds trivially. By (17a) we have that ⟨λ, u⟩ = 0. Therefore, ⟨u, µ -λ⟩ = ⟨u, µ⟩ ≥ 0, as u and µ are positive. ( 16) ⇒ (17

): Now let (u, λ) ∈ H 1 2 σ (R + , H 1 2 (G)) × H 1 2 σ (R + , H -1 2 (G)
) + be the solution to [START_REF] Eck | A symmetric boundary element method for contact problems with friction[END_REF]. Setting µ = 2λ and µ = 0 yields ⟨u, λ⟩ ≥ 0, ⟨u, λ⟩ ≤ 0. Therefore ⟨u, λ⟩ = 0. Assume that u is not ≥ 0. Then there exists µ ∈ H

1 2 σ (R + , H -1 2 (G))
+ such that ⟨u, µ⟩ < 0, and we obtain the contradiction 0 ≤ ⟨u, µ -λ⟩ = ⟨u, µ⟩ -⟨u, λ⟩ = ⟨u, µ⟩ < 0 .

Therefore u ∈ H 1 2 σ (R + , H 1 2 (G)) + . We now insert ṽ = v -u for u, v ∈ H 1 2 σ (R + , H 1 2 (G)) + into (16a). Note ⟨v -u, λ⟩ = ⟨v, λ⟩-⟨u, λ⟩ = ⟨v, λ⟩ ≥ 0. Inserting v -u in (16a), we have ⟨S σ u, v -u⟩ -⟨λ, v -u⟩ = ⟨h, v -u⟩ , or equivalently ⟨S σ u -h, v -u⟩ = ⟨λ, v -u⟩ ≥ 0.
The discrete formulation reads as follows:

Find (u ∆t 1 ,h 1 , λ ∆t 2 ,h 2 ) ∈ Ṽ 1,1 t 1 ,h 1 × (V 0,0 t 2 ,h 2 ) + such that { (a) ⟨S σ u ∆t 1 ,h 1 , v ∆t 1 ,h 1 ⟩ -⟨λ ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩ = ⟨h, v ∆t 1 ,h 1 ⟩ (b) ⟨u ∆t 1 ,h 1 , µ ∆t 2 ,h 2 -λ ∆t 2 ,h 2 ⟩ ≥ 0 (18) holds for all (v ∆t 1 ,h 1 , µ ∆t 2 ,h 2 ) ∈ Ṽ 1,1 t 1 ,h 1 × (V 0,0 t 2 ,h 2 ) + .
Like for the variational inequality (9), we assume that S σ is computed exactly. Note that, as in the elliptic case, we allow possibly different meshes for the displacement and the Lagrange multiplier. If the meshes for the Lagrange multiplier and the solution are sufficiently different, as in the elliptic case [START_REF] Babuska | On the mixed finite element method with Lagrange multipliers[END_REF], we obtain a discrete inf-sup condition in the space-time Sobolev spaces: Theorem 14. Let C > 0 sufficiently small, and max{h 1 ,∆t 1 } min{h 2 ,∆t 2 } < C. Then there exists α > 0 such that for all λ ∆t 2 ,h 2 :

sup µ ∆t 1 ,h 1 ⟨µ ∆t 1 ,h 1 , λ ∆t 2 ,h 2 ⟩ ∥µ ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * ≥ α∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ .
Proof. Let z be the solution to the equation z -∂ 2 t z -∆z = 0 in R + ×Ω with boundary conditions

∂z ∂t = 0 at t = 0, ∂z ∂ν = λ ∆t 2 ,h 2 on R + ×G, z = 0 on R + ×Γ \ G, z → 0 as t → ∞ .
We know from the coercivity of the Neumann-to-Dirichlet operator that

⟨z, λ ∆t,h 2 ⟩ = ⟨z, ∂ ν z⟩ ≥ C∥∂ ν z∥ 2 0,-1 2 ,σ = C∥λ ∆t 2 ,h 2 ∥ 2 0,-1 2 ,σ . (19)
Let δ such that ∥z∥ 0, 1 2 +δ,σ < ∞. From the approximation properties we note that there exists a function μ∆t 1 ,h 1 with

∥z -μ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * (max{h 1 , ∆t 1 }) δ ∥z∥ 0, 1 2 +δ,σ
. By the continuity of the Neumann-to-Dirichlet operator, the right hand side is estimated by

(h 1 + ∆t 1 ) δ ∥z∥ 0, 1 2 +δ,σ = (max{h 1 , ∆t 1 }) δ ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 +δ,σ . Using an inverse inequality for λ ∆t 2 ,h 2 , we conclude ∥z -μ∆t 1 ,h 1 ∥ 0, 1 2 ,σ,G, * (max{h 1 , ∆t 1 }) δ (min{h 2 , ∆t 2 }) δ ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ . ( 20
)
From the continuity of the Neumann-to-Dirichlet operator and (20), we obtain

∥μ ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * ≤ ∥z -μ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * + ∥z∥ 0, 1 2 ,σ, * (21)
(max{h 1 , ∆t 1 }) δ (min{h 2 , ∆t 2 }) δ ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ + ∥z∥ 0, 1 2 ,σ, * (max{h 1 , ∆t 1 }) δ (min{h 2 , ∆t 2 }) δ ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ + ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ . ( 22
)
Using μ∆t 1 ,h 1 from above and ( 21), we now estimate:

sup

µ ∆t 1 ,h 1 ⟨µ ∆t 1 ,h 1 , λ ∆t 2 ,h 2 ⟩ ∥µ ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * ≥ ⟨μ ∆t 1 ,h 1 , λ ∆t 2 ,h 2 ⟩ ∥μ ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * ⟨μ ∆t 1 ,h 1 , λ ∆t 2 ,h 2 ⟩ ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ = 1 ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ (⟨z, λ ∆t 2 ,h 2 ⟩ -⟨z -μ∆t 1 ,h 1 , λ ∆t 2 ,h 2 ⟩) .
The first term ⟨z, λ ∆t 2 ,h 2 ⟩ is estimated from below by ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ , with the help of ( 19), while for the second we have

⟨z -μ∆t 1 ,h 1 , λ ∆t 2 ,h 2 ⟩ ≤ ∥z -μ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * ∥λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ (max{h 1 , ∆t 1 }) δ (min{h 2 , ∆t 2 }) δ ∥λ ∆t 2 ,h 2 ∥ 2 0,-1 2 ,σ , (23)
where we used [START_REF] Gimperlein | Adaptive time domain boundary element methods and engineering applications[END_REF]. The assertion follows.

Theorem 15. The discrete mixed formulation (18) admits a unique solution. The following a priori estimates hold:

∥λ -λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ inf λ∆t 2 ,h 2 ∥λ -λ∆t 2 ,h 2 ∥ 0,-1 2 ,σ + (∆t 1 ) -1 2 ∥u -u ∆t 1 ,h 1 ∥ -1 2 , 1 2 ,σ, * , ( 24
) ∥u -u ∆t 1 ,h 1 ∥ -1 2 , 1 2 ,σ, * σ inf v ∆t 1 ,h 1 ∥u -v ∆t 1 ,h 1 ∥ 1 2 , 1 2 ,σ, * + inf λ∆t 2 ,h 2 { ∥ λ∆t 2 ,h 2 -λ∥ 1 2 ,-1 2 ,σ + ∥ λ∆t 2 ,h 2 -λ ∆t 2 ,h 2 ∥ 1 2 ,-1 2 ,σ } . ( 25
)
Proof. For a fixed mesh, the weak coercivity implies that the discretization of S σ is positive definite. The existence therefore follows from standard results for elliptic problems.

For the a priori estimate we first note that for arbitrary λ∆t 2 ,h 2 the following identity holds:

⟨λ ∆t 2 ,h 2 -λ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩ = ⟨S σ u ∆t 1 ,h 1 , v ∆t 1 ,h 1 ⟩ -⟨h, v ∆t 1 ,h 1 ⟩ -⟨ λ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩ = ⟨S σ u ∆t 1 ,h 1 , v ∆t 1 ,h 1 ⟩ -⟨S σ u, v ∆t 1 ,h 1 ⟩ + ⟨λ, v ∆t 1 ,h 1 ⟩ -⟨ λ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩ = ⟨S σ (u ∆t 1 ,h 1 -u), v ∆t 1 ,h 1 ⟩ + ⟨λ -λ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩, ( 26
)
where we made use of ( 16) and [START_REF] Gimperlein | Time domain boundary element methods for the Neumann problem: Error estimates and acoustic problems[END_REF]. By the inf-sup condition ( 14) and ( 26), we have:

α∥λ ∆t 2 h 2 -λ∆t 2 ,h 2 ∥ 0,-1 2 ,σ ≤ sup v ∆t 1 ,h 1 ⟨λ ∆t 2 ,h 2 -λ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩ ∥v ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * = sup v ∆t 1 ,h 1 ⟨S σ (u ∆t 1 ,h 1 -u), v ∆t 1 ,h 1 ⟩ + ⟨λ -λ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩ ∥v ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * .
We estimate both terms separately. From duality and an inverse inequality in time we obtain for the first term

|⟨S σ (u ∆t 1 ,h 1 -u), v ∆t 1 ,h 1 ⟩| ≤ ∥S σ (u ∆t 1 ,h 1 -u)∥ -1 2 ,-1 2 ,σ ∥v ∆t 1 ,h 1 ∥ 1 2 , 1 2 ,σ, * ∥u ∆t 1 ,h 1 -u∥ -1 2 , 1 2 ,σ, * (∆t 1 ) -1 2 ∥v ∆t 1 ,h 1 ∥ 0, 1 2
,σ, * . A similar argument applied to the second term yields:

|⟨λ -λ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩| ≤ ∥λ -λ∆t 2 ,h 2 ∥ 0,-1 2 ,σ ∥v ∆t 1 ,h 1 ∥ 0, 1 2 ,σ, * . We obtain the a priori estimate ∥λ -λ ∆t 2 ,h 2 ∥ 0,-1 2 ,σ σ inf λ∆t 2 ,h 2 ∥λ -λ∆t 2 ,h 2 ∥ 0,-1 2 ,σ + (∆t 1 ) -1 2 ∥u ∆t 1 ,h 1 -u∥ -1 2 , 1 2 ,σ, * . ( 27
)
Next we combine the Galerkin orthogonality

⟨S σ (u -u ∆t 1 ,h 1 ), V ∆t 1 ,h 1 ⟩ = ⟨λ -λ ∆t 2 ,h 2 , V ∆t 1 ,h 1 ⟩
with the coercivity of the Dirichlet-to-Neumann operator to obtain

∥u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ∥ 2 -1 2 , 1 2 ,σ, * σ ⟨S σ (u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ), u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ⟩ = ⟨S σ (u -v ∆t 1 ,h 1 ), u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ⟩ + ⟨S σ (u ∆t 1 ,h 1 -u), u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ⟩ = ⟨S σ (u -v ∆t 1 ,h 1 ), u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ⟩ + ⟨ λ∆t 2 ,h 2 -λ + λ ∆t 2 ,h 2 -λ∆t 2 ,h 2 , u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ⟩
for all v ∆t 1 ,h 1 and λ∆t 2 ,h 2 . From the mapping properties and the continuity of the dual pairing, we conclude

∥u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ∥ 2 -1 2 , 1 2 ,σ, * ∥u -v ∆t 1 ,h 1 ∥ 1 2 , 1 2 ,σ, * ∥u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ∥ -1 2 , 1 2 ,σ, * + ∥ λ∆t 2 ,h 2 -λ∥ 1 2 ,-1 2 ,σ ∥u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ∥ -1 2 , 1 2 ,σ, * + ∥ λ∆t 2 ,h 2 -λ ∆t 2 ,h 2 ∥ 1 2 ,-1 2 ,σ ∥u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ∥ -1 2 , 1 2 ,σ, * . Therefore, ∥u ∆t 1 ,h 1 -v ∆t 1 ,h 1 ∥ -1 2 , 1 2 ,σ, * σ ∥u -v ∆t 1 ,h 1 ∥ 1 2 , 1 2 ,σ, * + ∥ λ∆t 2 ,h 2 -λ∥ 1 2 ,-1 2 ,σ + ∥ λ∆t 2 ,h 2 -λ ∆t 2 ,h 2 ∥ 1 2 ,-1 2 ,σ . It follows that ∥u -u ∆t 1 ,h 1 ∥ -1 2 , 1 2 ,σ, * σ ∥u -v ∆t 1 ,h 1 ∥ 1 2 , 1 2 ,σ, * + ∥ λ∆t 2 ,h 2 -λ∥ 1 2 ,-1 2 ,σ + ∥ λ∆t 2 ,h 2 -λ ∆t 2 ,h 2 ∥ 1 2 ,-1 2 ,σ .

A variational inequality for the single layer operator

As a test case, we also consider the simpler punch problem, which models a rigid body indenting an elastic half-space. We consider the following formulation as a variational inequality for the single-layer operator V :

Find u ∈ H 1 2 σ (R + , H-1 2 (G)) + such that for all v ∈ H 1 2 σ (R + , H-1 2 (G)) + : ⟨p Q V u, v -u⟩ ≥ ⟨h, v -u⟩ . ( 28
)
As for the variational inequality for the Dirichlet-to-Neumann operator, a solution exists in the case where Ω = R n + is the half space:

Theorem 16. Let σ > 0 and h ∈ H 3 2 σ (R + , H 1 2 (G)).
Then there exists a unique classical solution u ∈ H [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF].

1 2 σ (R + , H-1 2 (G)) + of
Proof. See [START_REF] Cooper | Two variational inequality problems for the wave equation in a halfspace[END_REF], p. 456.

The corresponding discretized variational inequality reads as follows. Find u ∆t,h ∈ K+ t,h such that:

⟨p Q V u ∆t,h , v ∆t,h -u ∆t,h ⟩ ≥ ⟨h, v ∆t,h -u ∆t,h ⟩ (29)
holds for all v ∆t,h ∈ K+ t,h . The relevant a priori estimate reads:

Theorem 17. Let h ∈ H 3 2 σ (R + , H 1 2 (G)), and let u H 1 2 σ (R + , H-1 2 (G)) + , u ∆t,h ∈ K+ t,h
be the solutions of [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF], resp. [START_REF] Tallec | Solving dynamic contact problems with local refinement in space and time[END_REF]. Then the following estimate holds:

∥u -u ∆t,h ∥ 2 -1 2 ,-1 2 ,σ,⋆ σ inf ϕ ∆t,h ∈ K+ t,h (∥h -p Q V u∥ 1 2 , 1 2 ,σ ∥u -ϕ ∆t,h ∥ -1 2 ,-1 2 ,σ,⋆ + ∥u -ϕ ∆t,h ∥ 2 1 2 ,-1 2 ,σ,⋆ ) . ( 30
)
The proof proceeds analogous to the proof of Theorem 9. It replaces Theorem 5 for the Dirichlet-to-Neumann operator by the mapping properties of V in the half space, Theorem 4, and the coercivity ∥ϕ∥ 2

-1 2 ,-1 2 ,σ, *
σ ⟨p Q V ϕ, ϕ⟩ noted there. Similarly to the contact problem, for the numerical implementation a mixed formulation of the variational inequality (28) proves useful. Its analysis is analogous to the contact problem.

Theorem 18 (Mixed formulation)

. Let h ∈ H 3 2 σ (R + , H 1 2 (G)).
The variational inequality formulation ( 28) is equivalent to the following formulation:

Find (u, λ) ∈ H 1 2 σ (R + , H-1 2 (G)) × H 1 2 σ (R + , H 1 2 (G)) + such that { (a) ⟨V u, v⟩ -⟨λ, v⟩ = ⟨h, v⟩ (b) ⟨u, µ -λ⟩ ≥ 0, ( 31
)
for all (v, µ) ∈ H 1 2 σ (R + , H-1 2 (G)) × H 1 2 σ (R + , H 1 2 (G)) + .
The discrete formulation reads as follows:

Find (u ∆t 1 ,h 1 , λ ∆t 2 ,h 2 ) ∈ V 1,1 t 1 ,h 1 × (V 1,1 t 2 ,h 2 ) + such that { (a) ⟨V u ∆t 1 ,h 1 , v ∆t 1 ,h 1 ⟩ -⟨λ ∆t 2 ,h 2 , v ∆t 1 ,h 1 ⟩ = ⟨h, v ∆t 1 ,h 1 ⟩ (b) ⟨u ∆t 1 ,h 1 , µ ∆t 2 ,h 2 -λ ∆t 2 ,h 2 ⟩ ≥ 0 (32) holds for all (v ∆t 1 ,h 1 , µ ∆t 2 ,h 2 ) ∈ V 1,1 t 1 ,h 1 × (V 1,1 t 2 ,h 2 ) + .
8 Algorithmic considerations

Marching-on-in-time scheme for the variational equality

A first step towards the solution of the variational inequality is to find the solution of the Dirichletto-Neumann equation. Using the representation W -

(K ′ -1 2 I)V -1 (K -1 2 I
) for the Dirichlet-to-Neumann operator in terms of layer potentials, the Dirichlet-to-Neumann equation is equivalently reformulated as follows: For given h ∈ H

3 2 σ (R + , -1 2 (G)), find u ∈ H 1 2 σ (R + , H 1 2 (G)), v ∈ H 1 2 σ (R + , H-1 2 (G)) such that ∞ ∫ 0 ⟨W u -(K ′ -1 2 )v, ϕ⟩ G dt = ∞ ∫ 0 ⟨h, ϕ⟩ G dt , ( 33
) ∞ ∫ 0 [⟨V v, Ψ⟩ G -⟨(K -1 2 )u, Ψ⟩ G ] dt = 0, (34)
holds for all ϕ ∈ H

1 2 σ (R + , H 1 2 (G)),Ψ ∈ H 1 2 σ (R + , H-1 2 (G)). For the discretization, let u ∆t,h = ∑ m,i c m i β m ∆t (t)ξ i h (x, y) ∈ Ṽ 1,1 t,h , v ∆t,h = ∑ m,i d m i β m ∆t (t)ξ i h (x, y) ∈ V 1,1 t,h
be ansatz functions that are linear in space and time. To obtain a stable marching-on-in-time scheme we test the first equation against constant test functions in time and the second equation against the time derivative of constant test functions. We let

ϕ ij ∆t,h := γ [t i-1 ,t i ] (t)ξ j h (x, y) =: γ i (t)ξ j (x, y), Ψij
∆t,h = γi (t)ξ j (x, y) be test functions that are constant in time and linear in space. Thus, the discrete system reads as follows:

∞ ∫ 0 ⟨W u ∆t,h -(K ′ -1 2 )v ∆t,h , γ n ξ j ⟩ G dt = ∞ ∫ 0 ⟨h, γ n ξ j ⟩ G dt , ∞ ∫ 0 [⟨V v ∆t,h , γn ξ j ⟩ G -⟨(K -1 2 )u ∆t,h , γn ξ j ⟩ G ] dt = 0,
for all n = 1 . . . , N t , j = 1, . . . N s . Setting I j = I, if j = 0, 1, I j = 0 otherwise and Î0 = (-I), Î1 = I, Îj = 0 otherwise, where I is the mass matrix, we may rewrite the system as

M j := ( W j -((K j ) ′ -1 2 ∆t 2 I j ) -(K j -1 2 Îj ) V j
) , for all j = 2, . . . , N t (see [START_REF] Gimperlein | Time domain boundary element methods for the Neumann problem: Error estimates and acoustic problems[END_REF] for further details). From the convolution structure in time, we obtain: For arbitrary n ∈ {1, . . . , N t }:

∞ ∑ m=1 M n-m ( c m d m ) = ( ∆t 2 I(h n-1 + h n ) 0 ) . ( 35
)
By causality, M j = 0 when j < 0, so that the sum on the left hand side ends at m = n. This results in the time stepping scheme:

M 0 ( c n d n ) = ( ∆t 2 I(h n-1 + h n ) 0 ) - n-1 ∑ m=1 M n-m ( c m d m ) . ( 36
)
Hence, if we save the matrices from previous time steps, we only need to calculate one new matrix M n-1 in time step n to obtain the vector [c n , d n ] T .

Space-time Uzawa algorithm

The solution of the discrete mixed formulation (18) may be computed using a Uzawa algorithm in space-time.

Algorithm 1 Space-time Uzawa algorithm choose ρ > 0: k = 0 : y (0) = ⃗ 0 while stopping criterion not satisfied do solve:

Sx (k) = h + y (k) compute: y (k+1) = Pr K (y (k) -ρx (k)
), where (Pr K y) i = max{y i , 0} k ← k + 1 end while Lemma 19. The space-time Uzawa algorithm converges, provided that 0 < ρ < 2C σ . Here C σ is the coercivity constant in Theorem 5.

Proof. From the algorithm and the contraction property of the projection Pr K , ∥λ

(k+1) ∆t,h -λ ∆t,h ∥ 2 0,0,σ = ∥Pr K (λ (k) ∆t,h -ρu (k) ∆t,h ) -Pr K (λ ∆t,h -ρu ∆t,h )∥ 2 0,0,σ ≤ ∥λ (k) ∆t,h -λ ∆t,h -ρ(u (k) ∆t,h -u ∆t,h )∥ 2 0,0,σ = ∥λ (k) ∆t,h -λ ∆t,h ∥ 2 0,0,σ -2ρ⟨λ (k) ∆t,h -λ ∆t,h , u (k) ∆t,h -u ∆t,h ⟩ + ρ 2 ∥u (k) ∆t,h -u ∆t,h ∥ 2 0,0,σ .
We conclude that

∥λ (k) ∆t,h -λ ∆t,h ∥ 2 0,0,σ -∥λ (k+1) ∆t,h -λ ∆t,h ∥ 2 0,0,σ ≥ 2ρ⟨λ (k) ∆t,h -λ ∆t,h , u (k) ∆t,h -u ∆t,h ⟩ -ρ 2 ∥u (k) ∆t,h -u ∆t,h ∥ 2 0,0,σ . Further note that ⟨λ (k) ∆t,h -λ ∆t,h , u (k) ∆t,h -u ∆t,h ⟩ = ⟨S σ (u (k) ∆t,h -u ∆t,h ), u (k) ∆t,h -u ∆t,h ⟩ ≥ C σ ∥u (k) ∆t,h -u ∆t,h ∥ 2 -1 2 , 1
As ∥u

(k) ∆t,h -u ∆t,h ∥ -1 2 , 1 2 ,σ, * ≥ ∥u (k) ∆t,h -u ∆t,h ∥ 0,0,σ , we conclude that ∥λ (k) ∆t,h -λ ∆t,h ∥ 2 0,0,σ -∥λ (k+1) ∆t,h -λ ∆t,h ∥ 2 0,0,σ ≥ (2ρC σ -ρ 2 )∥u (k) ∆t,h -u ∆t,h ∥ 2 0,0,σ .
The right hand side is non-negative provided 0 < ρ < 2C σ . We conclude that ∥λ

(k) ∆t,h -λ ∆t,h ∥ 0,0,σ is a decreasing sequence. As ∥λ (k)
∆t,h -λ ∆t,h ∥ 0,0,σ ≥ 0, it therefore converges, and it follows that ∥u

(k) ∆t,h -u ∆t,h ∥ 0,0,σ → 0.
Practically, we choose the test functions of the mixed discretized formulation (18) to be piecewise constant in time to obtain a marching-on-in-time scheme. In this case, computing the coefficients x (k) = ({x m i } Ns,Nt i,m=1 ) (k) corresponds to solving the equality for the Dirichlet-to-Neumann operator in Section 8.1, with an augmented right hand side:

⟨S ( ∑ m,i (x m i ) k ) , γ n ϕ j ⟩ = ⟨ ∑ m,i (y m i ) k-1 β m ξ i , γ n ϕ j ⟩ + ⟨ ∑ m,i h m i β m ξ i , γ n ϕ j ⟩, (37)
for all j = 1, . . . N s , n = 1, . . . , N t . Comparing this system with the system (36) we observe that (37) may be written as:

n ∑ m=1 M n-m ( (x m ) k (d m ) k ) = ( ∆t 2 I(h n-1 + h n ) 0 ) + ( ∆t 2 I((y n-1 ) k-1 + (y n ) k-1 ) 0 ) , (38)
for n = 1, . . . , N t .

Remark 20. For a Lagrange multiplier that is constant in time, i.e.

λ ∆t,h = ∑ m,l y m l γ m (t)ξ l (x),
the second term on the right hand side of (38) becomes ∆t

( I(y n ) k-1 0
) . This will be relevant in the following section.

Time-step Uzawa algorithm

The space-time Uzawa algorithm solves the whole space-time system in every Uzawa iteration. While it is provably convergent, a time-stepping Uzawa algorithm turns out to be more efficient and will be useful for future adaptive computations.

A crucial observation to derive a time-stepping algorithm is that L2 (Q) continuously embeds in

H 1 2 σ (R + , H -1 2 (G)).
One may therefore use piecewise constant ansatz functions in time,

λ ∆t,h = ∑ m,i y m i γ m (t)ϕ i (x) ∈ L 2 (Q) ⊂ H and note that ∑ m,i y m i ∞ ∫ 0 ∫ Γ γ m ϕ i γ n ϕ j ds x dt = ∆t ∑ i y n i ∫ Γ ϕ i ϕ j ds x .
The resulting space-time variational inequality reads:

      M 0 M 1 M 0 . . . M Nt • • • M 0             [c 1 , d 1 ] t [c 2 , d 2 ] t . . . [c Nt , d Nt ] t       -∆t       [ Ĩy 1 , 0] t [ Ĩy 2 , 0] t . . . [ Ĩy Nt , 0] t       = ∆t 2       [I(h 0 + h 1 ), 0] t [I(h 1 + h 2 ), 0] t . . . [I(h Nt-1 + h Nt ), 0] t       , ( 39
) ⟨ ∑ m,i c m i β m ξ i , µ ∆t,h - ∑ m,i y m i γ m ϕ i ⟩ ≥ 0, ∀µ ∆t,h ∈ (V 0,0 t,h ) + . ( 40
) With µ ∆t,h = ∑ m,i µ m i γ m (t)ϕ i (x)
, we first consider the space-time variational inequality (40) for a fixed n = 1 . . . , N t :

⟨ ∑ m,i c m i β m ξ i , ∑ j (µ n j -y n j )γ n ϕ j ⟩ = ∆t 2 ⟨ ∑ i c n i ξ i , ∑ j (µ n j -y n j )ϕ j ⟩ + ∆t 2 ⟨ ∑ i c n-1 i ξ i , ∑ j (µ n j -y n j )ϕ j ⟩ , = (c n ) ⊤ Î(µ n -y n ) + (c n-1 ) ⊤ Î(µ n -y n ).
Here we used

⟨β m ξ i , γ n ϕ j ⟩ = ∫ Γ ξ i (x)ϕ j (x) ∞ ∫ 0 β m (t)γ n (t)dt ds x = ∫ Γ ξ i (x)ϕ j (x) ∆t 2 [δ n,m + δ n-1,m ]ds x ,
with δ n,m = 1 if n = m and = 0 otherwise. Also Îij = ∆t 2 I ij . Therefore (40) may be written as (c 1 ) ⊤ Î(µ 1 -y 1 ) + (c 0 ) ⊤ Î(µ 1 -y 1 ) + (c 2 ) ⊤ Î(µ 2 -y 2 ) + (c 1 ) ⊤ Î(µ 2 -y 2 ) + . . .

+ (c Nt ) ⊤ Î(µ Nt -y Nt ) + (c Nt-1 ) ⊤ Î(µ Nt -y Nt ) ≥ 0 ∀µ j ∀j ,
with c 0 = 0. Setting µ = (y 1 . . . , y k-1 , µ ′ , y k+1 , . . . , y Nt ) for k = 1, . . . , N t yields

(c k ) ⊤ Î(µ ′ -y k ) + (c k-1 ) ⊤ Î(µ ′ -y k ) ≥ 0, ∀µ ′ ,
and we see that the solution to the space-time variational inequality satisfies the following timestepping scheme:

For k = 1, . . . , N t , find (c k , d k , y k ) such that M 0 ( c k d k ) -∆t ( Ĩy k 0 ) = ( ∆t 2 I(h k-1 + h k ) 0 ) - k-1 ∑ m=1 M k-m ( c m d m ) , ( 41
) (c k ) ⊤ Î(µ k -y k ) ≥ -(c k-1 ) ⊤ Î(µ k -y k ) ∀µ k . ( 42
)
Conversely, if we have solutions to (41) and (42) for every k = 1, . . . , N t , summing (42) yields (40). We conclude:

Lemma 21. The solution to the space-time variational inequality is also a solution to the timestep variational inequality and vice versa. the error between the numerical approximations and the exact solution converges to 0 at a rate 0.53, as shown in Figure 3. Here, we compute the experimental convergence rate α as

α = log E(u 1 ) -log E(u 2 ) log DOF 1 -log DOF 2 .
In terms of h, the convergence rate is 1.6. The results exemplify the convergence of our proposed method to approximate the Dirichlet-to-Neumann operator.

Contact problem: Dirichlet-to-Neumann operator

We now consider the discretization of the nonlinear contact problem (5) for both flat and more general contact areas. In this case, no exact solutions are known, and we compare the numerical approximations to a reference solutions on an appropriately finer space-time mesh. We use the discretisation from Section 8.1, in particular with a Lagrange multiplier that is piecewise linear in space and time. The numerical solutions are compared to a reference solution on a mesh with 12800 triangles, and we use ∆t = 0.05.

The inequality is solved using the space-time Uzawa algorithm. We stop the solver when subsequent iterates have a relative difference of less than 10 -11 .

Figure 4 shows the solution u ∆t,h to the contact problem (left column) and the corresponding Lagrange multiplier λ ∆t,h (right column) for several time steps. The solution gradually develops a smooth bump in the center, which gives rise to a radially outgoing wave. A nonvanishing Lagrange multiplier λ = S σ u -h indicates the emergent contact forces at the depicted times t = 4.25 and t = 5.

Figure 5 considers the relative error to the reference solution in L 2 ([0, T ] × G). The numerical approximations converge at a rate of α = 0.7 with increasing degrees of freedom. This compares to the rate of convergence for the Dirichlet-to-Neumann equation in Example 23, where also the geometry needed to be approximated. Algorithmically, the computational cost of the nonlinear solver is dominated by the cost of computing the matrix entries. Next we compare the above space-time Uzawa algorithm with the time-step Uzawa variant. In a given time step, we stop the Uzawa iteration when subsequent iterates either have a relative difference of less than 10 -12 or if the ℓ ∞ -norm is less than 10 -10 . For the space-time Uzawa algorithm we use the stopping criterion from before.

Example 25. We consider the contact problem [START_REF] Banjai | Time-domain Dirichlet-to-Neumann map and its discretization[END_REF] with the geometry and right hand side from Example 24. On a fixed mesh of 3200 triangles, and with ∆t = 0.1 and T = 6, we investigate the difference of the approximate solutions obtained from the space-time and time-step Uzawa algorithms.

Figure 6 shows the temporal evolution of the relative difference in L 2 (Γ) between the two methods in a semi-logarithmic plot. The difference is smaller than 0.01% for all times, but increases sharply around the onset of contact shortly after time t = 4. To illustrate our method for non-flat contact geometries, we consider a cube with three contact faces. Physically, one may think of a rigid cube which is tightly fixed to an elastic surrounding material material on its three other faces. The error between the benchmark and approximate solutions on coarser meshes is shown in Figure 8. It shows a convergence with an approximate convergence rate of α = 0.6.

Single-layer potential

We finally address the punch problem described by the variational inequality [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF]. Unlike for the problems involving the Dirichlet-to-Neumann operator, in this case we do not need to approximate the integral operator. As in the previous numerical examples, we consider both flat and more general contact geometries.

Example 27.

For the punch problem [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF] we choose Γ = [-2, 2] 2 × {0} and the area of contact G = [-1.2, 1.2] 2 × {0}. We set T = 5 and keep the CFL-ratio fixed at ∆t h ≈ 0.7. As right hand side we again consider h(t, x) = e -2t t 4 cos(2πx) cos(2πy)χ [-0.25,0.25] (x)χ [-0.25,0.25] (y), and we look for a numerical solution and Lagrange multiplier which are piecewise constant in time, linear in space. The approximate solutions are compared to a benchmark obtained from a mesh with 12800 triangles and ∆t = 0.05.

We solve the variational inequality using the time-step Uzawa algorithm. In a given time step, we stop the Uzawa iteration when subsequent iterates either have a relative difference of less than 10 -12 or if the ℓ ∞ -norm is less than 10 -10 .

Figure 9 shows the solution u ∆t,h to the punch problem and its Lagrange multiplier λ ∆t,h . Contact is observed for most of the considered time interval. Because of the potentially low spatial regularity of the solution, which a priori only belongs to a Sobolev space with negative exponent, we do not consider the error of the numerical solutions in L 2 ([0, T ]×Γ). As a weaker measure, we consider an approximate energy norm. An approximation for the energy is computed via E(u) = 1 2 u ⊤ V u -u ⊤ h. Figure 10 shows the convergence of the numerical solutions in energy. The energy similarly shows a convergence rate of α = 1. We finally consider the punch problem with the entire surface of the cube as contact area.

Example 28. For the punch problem [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF] we choose G = Γ to be the surface of [-2, 2] 3 . We set T = 4 and keep the CFL-ratio fixed at ∆t h ≈ 0. The variational inequality is solved using the time-step Uzawa algorithm as in Example 27. Figure 11 shows the solution u ∆t,h to the punch problem and its Lagrange multiplier λ ∆t,h on the surface of the cube. Again contact is observed for most times. As in the previous example, in Figure 12 we show the convergence in energy. We here compute a convergence rate of α = 0.5.

Conclusions

In this work we propose and analyze a Galerkin boundary element method to solve dynamic contact problems for the wave equation. Boundary elements provide a natural and efficient formulation, as the contact takes place at the interface between two materials.

Analytically, we obtain a first a priori error analysis for a variational inequality involving the Dirichlet-to-Neumann operator, as well as a similar analysis for a mixed formulation. The analysis and the stability of the method are crucially based on the weak coercivity of the formulation, and for the mixed method an inf-sup condition in space-time. The proof requires a flat contact area, as only in this case the existence of solutions to the continuous problem is known. Also a variational inequality and a mixed formulation for the single layer operator are considered, which do not require the approximation of the operator.

Numerical experiments demonstrate the efficiency and convergence of the proposed mixed method. A time-stepping Uzawa method for the solution of the variational inequality proves more efficient, but also potentially less stable than a similar solver for the space-time system. The latter is shown to be provably convergent. As a key point, the numerical experiments indicate stability and convergence beyond flat geometries.

The current work provides a first, rigorous step towards efficient boundary elements for dynamic contact. Future work will focus on the a posteriori error analysis, which is essential for adaptive mesh refinements to resolve the singularities of the solution in space and time [START_REF] Gimperlein | A residual a posteriori estimate for the time-domain boundary element method[END_REF], as well as on stabilized mixed space-time formulations [START_REF] Banz | Stabilized mixed hp-BEM for frictional contact problems in linear elasticity[END_REF][START_REF] Burman | Galerkin finite element methods with symmetric pressure stabilization for the transient Stokes equations: stability and convergence analysis[END_REF]. For applications to traffic noise [START_REF] Banz | Time domain BEM for sound radiation of tires[END_REF], also the nonsmooth variational inequalities for frictional contact will be of interest.

Figure 1 :

 1 Figure 1: L 2 (Γ)-norm of the solution to Su = h for fixed CFL≈ 0.6.

Figure 2 :Figure 3 :

 23 Figure 2: Absolute error |∥u∥ L 2 (Γ) -∥u ∆t,h ∥ L 2 (Γ) | as a function of time for fixed ∆t h .

Example 24 .

 24 We choose Γ = [-2, 2] 2 × {0} with contact area G = [-1, 1]2 × {0} for times up to T = 5, with the CFL-ratio ∆t h = 0.7. The right hand side of the contact problem (5) is given by h(t, x) = e -2t t 4 cos(2πx) cos(2πy)χ [-0.25,0.25] (x)χ [-0.25,0.25] (y) .

Figure 4 :

 4 Figure 4: Evolution of u and λ in G = [-1, 1] 2 × {0} for the contact problem, Example 24.

Figure 5 :

 5 Figure 5: Relative L 2 ([0, T ] × Γ)-error vs. degrees of freedom of the solutions to the contact problem for fixed ∆t h , Example 24.

  vs. Time-Step Uzawa on 3200 triangles

Figure 6 :

 6 Figure 6: Relative L 2 (Γ)-error between the solutions of the space-time and time-step Uzawa algorithms, Example 25.

Example 26 .

 26 Let Γ be the surface of the cube [-2, 2] 3 , with contact area G consisting of the top, front and right faces. We set T = 6 and keep the CFL-ratio fixed at ∆t h ≈ 0.7. On each of the contact faces we prescribe a right hand sideh(t, x) = e -2t t 4 cos(2πx) cos(2πy)χ [-0.25,0.25] (x)χ [-0.25,0.25] (y) ,centered in the midpoint of each face. The reference solution is computed on a mesh with 19200 triangles, for ∆t = 0.1.The evolution of u ∆t,h and the Lagrange multiplier λ ∆t,h on the top face of the cube are depicted in Figure7. The nonzero displacement u ∆t,h spreads from the neighboring contact faces into the shown area by time t = 5 and eventually leads to strong contact near the upper left corner.

Figure 7 :

 7 Figure 7: Evolution of u and λ on [-2, 2] 2 × {2} for the contact problem on [-2, 2] 3 , Example 26.

Figure 8 :

 8 Figure 8: Relative L 2 ([0, T ] × Γ)-error vs. degrees of freedom for the contact problem for fixed ∆t h , Example 26.

Figure 9 :

 9 Figure 9: Evolution of u and λ for the punch problem with G = [-1.2, 1.2] 2 × {0}, Example 27.

Figure 10 :

 10 Figure 10: Relative energy error for the punch problem for fixed ∆t h , Example 27.

Figure 11 :

 11 Figure 11: Evolution of the density and the Lagrange multiplier for the punch problem, Example 28.

7 .

 7 As right hand side we consider the same function as in Example 26, h(t, x) = e -2t t 4 cos(2πx) cos(2πy)χ [-0.25,0.25] (x)χ [-0.25,0.25] (y) , centered in the midpoint of the top, front and right face. We approximate the reference solution on a mesh with 19200 triangles and ∆t = 0.1.

Figure 12 :

 12 Figure 12: Relative error of the energy for the solution to the punch problem for fixed CFL≈ 0.7.

,σ, * .

σ (R + , H -1 2 (G)), for the Lagrange multiplier. Set Ĩi,j = ∫ Γ ϕ i ϕ j ds x , I i,j = ∫ Γ ξ i ϕ j ds x ,

H. G. acknowledges support by ERC Advanced Grant HARG 268105 and the EPSRC Impact Acceleration Account.

The resulting time-step Uzawa algorithm is as follows:

Algorithm 2 Time-step Uzawa algorithm choose ρ > 0: for n=1,. . . , N t do k=0: (y n ) 0 = ⃗ 0 while stopping criterion not satisfied do solve:

The Uzawa algorithm converges in each time step, as long as M 0 is positive definite, and ρ is sufficiently small.

Remark 22.

As for the contact problem, the mixed formulation for the punch problem (31) may be discretized and solved with the above space-time or time-step Uzawa algorithms. The Dirichlet-to-Neumann operator is here replaced by the single-layer operator.

Numerical experiments

As stated above, we set σ = 0 in our computations and discretize the Dirichlet-to-Neumann operator as in Section 8.1.

Dirichlet-to-Neumann operator on unit sphere

Example 23. We solve the discretized variational equality [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF] for the Dirichlet-to-Neumann operator on Γ = S 2 , with a right hand side obtained from the Neumann data of a known, radially symmetric solution to the wave equation. Specifically,

where H is the Heaviside function. The solution u of the Dirichlet-to-Neumann equation Su = h corresponds to the Dirichlet data of the solution to the wave equation. Hence,

We fix the CFL ratio ∆t h ≈ 0.6 and set T = 5.

Figure 1 shows the L 2 (Γ)-norm of the exact, resp. numerical solution, while Figure 2 depicts the error in this norm. As the number of degrees of freedom increases, the L 2 ([0, T ] × Γ)-norm of


D M Zhang

H W Huang

P Y Hicher

Numerical Prediction of Long-term Settlements over Tunnels in Clay

The long-term settlements over tunnels can be rather significant, particularly when tunnels are embedded in soft and compressible soils. The long-term settlement can be predicted based on the stress-strain-time relationship of soft clay using multiple-step FEM. The influence of drainage condition of tunnel linings on the evolution of long-term surface settlement, as well as ground loss, are simulated with the multiple-step numerical procedure, and some conclusions are drawn from the numerical simulations. Finally, a case study is also performed and a satisfaction agreement between the predicted long-term settlement and the observed one is obtained.

INTRODUCTION

The long-term settlement can be significant, particularly when tunnels are embedded in soft and compressible soils. It is concluded by Shirlaw, based on comprehensive field data, that typically the increase in settlements over long-term is of the order of 30%-90% of the total settlement, and that in many cases a widened settlement trough develops [START_REF] Shirlaw | Observed and calculated pore pressures and deformations induced by an earth balance shield: Discussion1[END_REF]. This is confirmed, for example, by the observed settlements over a period of 11 years reported by [START_REF] O'reilly | Long-term settlements over tunnels: an elevenyear study at Grimsby[END_REF] for a 3m diameter tunnel constructed in normally consolidated silty clay in Grimsby, England (shown as Figure 1a). The final equilibrium of settlement was achieved after about 10 years, and the corresponding transverse surface profiles show a significant widening with time. Similar widening has been reported by a number of authors [START_REF] Glossop | Settlement caused by tunnelling through soft marine clay[END_REF][START_REF] Shirlaw | Observed and calculated pore pressures and deformations induced by an earth balance shield: Discussion1[END_REF]. The increase of long-term settlement reached 60-90% of the total settlement for No.2 metro tunnel in Shanghai, meanwhile the width of the settlement rose from 2.645a in the first month to 2.742a (a is the radius of tunnel) in the third month after construction (Lee et al. 1998, shown in Figure 1b). The long-term settlement can affect the function of the tunnel and the safety of the surroundings.

It is widely accepted that the effects of time on strength and deformation characteristics of soils, especially of saturated cohesive soft soils, plays an important role in a variety of geotechnical engineering activities where either rapid or long-term response is concerned [START_REF] Liang | A unified elasto-viscoplasticity model for clays, part I: theory[END_REF]. Accordingly, it is necessary to build a constitutive law that can be used to predict the time-dependent deformation around the tunnel. Consequently, an elastic-viscoplastic constitutive model coupled with Biot's consolidation theory was developed and used in the long-term settlement prediction over urban tunnels.

ELASTIC VISCO-PLASTIC CONSTITUTIVE MODEL

Visco-plastic flow rule

The visco-plastic constitutive model was established using the results of both triaxial and oedometer tests (the experimental results are given in [START_REF] Zhang | Time-dependent stress-strain behaviour of natural soft clay[END_REF]. The construction of the viscoplastic model was undertaken based on the framework of the Modified Cam-clay model [START_REF] Roscoe | On the generalized stress-strain behaviour of wet clay[END_REF]. The viscoplastic strain rate follows the framework proposed by Perzyna with associated flow rule [START_REF] Perzyna | The constitutive equation for rate sensitive plastic materials[END_REF][START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF].

The constitutive model is then implemented into the FEM code CESAR_LCPC coupled with Biot's consolidation theory. Consequently, the long-term response of urban tunnel can be numerically estimated. The numerical simulation on long-term settlement over urban tunnels is performed with the background of Shanghai metro line 2. The tunnel lining for Shanghai metro line 2 is 6.2m in external diameter and 5.5m in internal diameter. The average depth to the center-line of the tunnel is about 11m, which is about 1.8D (D is the tunnel diameter). The tunnel was excavated by EPB shield machine. The shield body is 6.24m long with 6.34m in diameter. Correspondingly, the clearance between the external diameters of shield body and the tunnel, which is usually named physical gap G p according to the definition of [START_REF] Lee | Subsidence due to tunnelling: Estimating the gap parameter[END_REF], will reach 140mm. The effective back grouting pressure is determined as 0.25MPa according to the in-situ measurements.

NUMERICAL SIMULATION PROCEDURE

The numerical simulation was performed with 8-node isoparametric finite elements under the assumption of plane strain conditions. Correspondingly, the FEM mesh of the tunnel is shown in Fig. 2.

Boundary conditions

The boundary conditions in the numerical simulation on long-term settlements contain the following two types, one is the displacement boundary condition, and the other is the drainage condition.

A free displacement boundary condition was adopted at the ground surface. It is assumed that no horizontal nor vertical displacement takes place at the lower boundary, for it is beyond the influence of tunnel construction. The lateral displacements at left-and right-hand boundary are both fixed as zero. The left-hand boundary is the line of symmetry.

The drainage condition is not a constant, but a variable during the different development stages of short-term and long-term settlements. The drainage condition at the ground surface is assumed to be free, hence the excess pore pressure will be kept as zero along the ground surface; meanwhile the lower boundary as well as the left-hand boundary condition are considered to be impermeable, in order to prevent the pore pressure dissipation across the boundary during the settlement development. However, the right-hand boundary is taken as impermeable, because of the low permeability of the soil around the tunnel during the development of short-term settlement, while the pore-pressure is assuned to maintain the hydrostatic pore-pressure at the right-hand boundary during the long-term settlement evolution. That means the right-hand boundary is taken as impermeable during short-term settlement and if the right-hand boundary is taken as a constant head boundary it is permeable in one direction only during the long-term settlement. The tunnel lining is also assumed to be an impermeable boundary during the short-term settlement, but permeable during the long-term settlement; however, this is questionable and was studied in the long-term settlement simulation.

The initial effective stresses and hydrostatic pore pressure were calculated based on the weight of the soil and the underground water condition.

Numerical simulation procedure

The evolution of long-term (or post-construction) settlement was simulated using 4 consecutive steps. The 4 consecutive steps correspond to different stages of the development of surface settlement over tunnels, and their correlations are as follows: STEP 0 means the initial state of the ground before taking account of the influence of tunnelling, STEP 1 represents the heave (or settlement) of the ground before the shield machine reaches, STEP 2 represents the shield machine just passes the studied position, STEP 3 is the final case in the simulation, where the long-term settlement occurs. The numerical results for each step are the initial state of the next step since the 4 steps are consecutive both in time and stress state. Among the 4 steps, STEP 0, 1 and 3 are simulated using the stress boundary conditions, while the displacement boundary condition is adopted in STEP 2 to simulate the closure of the tail void during the tunnel advances.

Determination of parameters

The parameters involved in the numerical simulation can be divided into two types. One type is the mechanical parameters of soils and tunnel segments. The other type is related to the tunnel construction process. The mechanical parameters of soils can be obtained through laboratory triaxial and oedometer tests, and are listed in table 1. Considering the influence of shield tunneling on the surrounding soil, the strength of the soils is discounted to some extent, based on the in-situ measurement, in step 2. The tunnel construction parameters include the advance rate of shield machine, the physical gap and the grouting pressure are presented in table 2.

4.DEVELOPMENT OF LONG-TERM SETTLEMENTS

The development of long-term settlement is dependent on many factors such as drainage condition of the tunnel linings, back-grouting, over-excavation, surcharge over the tunnel, time-dependent properties of soils around the tunnel and so on. Here, only the influence of drainage condition on the evolution of long-term settlements over tunnels are considered.

Influence of permeability of the tunnel on the evolution of long-term settlements

O' Reilly et al.(1991) and [START_REF] Shin | A numerical study of the effect of groundwater movement on long-term tunnel behavior[END_REF] concluded that the drainage condition of the tunnel had a significant influence on both the magnitude and distribution of long-term settlement. The permeability of the tunnel is described by the relative value to the permeability of surrounded soil in the multiplestep numerical simulation. Accordingly, the influence of permeability of tunnels is studied for the following four different cases: (1) impermeable,

∞ =

impermeable K t /K s =0.01 K t /K s =0.1 K t /K s =1
Fig. 3 Long-term surface settlements for Fig. 4 Long-term transverse surface settlement different lining permeability with time troughs for different lining permeabilities Fig. 3 reveals that the magnitude of maximum surface settlement is smaller at any time and reaches stability more quickly with an impermeable lining. The higher the permeability of tunnel lining, the more the surface settlement.

The long-term transverse surface settlement troughs, normalized to the corresponding maximum surface settlements, are shown in Fig. 4.

The normalized long-term transverse surface settlement troughs shown in Fig. 4 are independent of the permeability and exhibit uniform shape. Additionally, they can be described by Peck's empirical formu -lation (equation 1) [START_REF] Peck | Deep excavation and tunneling in soft ground[END_REF].

( )

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = 2 2 max 2 exp i x S x S (1)
However, the evolution of the width of the surface settlement troughs ( i ), estimated using equation 2, is closely associated with the permeability of tunnels linings, as shown in Fig. 5.

( ) max 606 . 0 S i S = (2)
It can be clearly found from Fig. 5 that i continuously increases with both time and the permeability of tunnel linings until it reaches the stability during the evolution of long-term settlements under drainage condition. However, i increases with time within a certain period, for example 300 days in Fig. 5, and thereafter it begins to decrease when time exceed this period prove the tunnel linings are impermeable.

The corresponding predicted ground loss with different permeability of tunnel linings is presented in Fig. 6 via time.

The evolution of ground loss with time corresponding to different permeability is similar with that of the width of long-term transverse surface settlement troughs i .

A BACK-ANALYSIS

A back-analysis was performed on the Shanghai metro No.2. The FEM meshes, the boundary and drainage conditions, and numerical procedures were as presented above. The calculated surface settlements and the observed ones are shown in Fig. 7.

In Fig. 7, the creep parameters used for the various analysis are:

, 10 1 Fig. 8 The ratio of max. immediate settlement Fig. 9 Long-term surface settlement troughs to max. long-term settlement corresponding to different tunnel excavation methods

The calculated and observed long-term settlements of Shanghai metro as well as those found in the literature are presented in Fig. 8 for a full understanding of the evolution of the long-term settlement of Shanghai metro.

Fig. 8 shows that the long-term settlements of Shanghai metro were significant, and couldn't be neglected. The transverse long-term surface settlement troughs corresponding to different tunnel excavation methods are presented in Fig. 9, to illustrate the development of surface settlement troughs in the long-term behavior of tunnels in clays.

Fig. 9 shows that the surface settlement troughs do not change much during the long-term if the EPB tunnelling is adopted, and can be well described by Peck's empirical formula (equation 2). However, the shape of long-term surface settlement troughs varies significantly between TBM and open shield methods. Hence, the evolution of long-term surface settlement troughs seems to be dependent on the tunnel excavation method, and it is consistent with the findings of [START_REF] Mair | Theme lecture: Bored tunneling in the urban environment[END_REF] and [START_REF] Shirlaw | Pore pressure around tunnels in clay: 1Discussion[END_REF].

CONCLUSIONS

Coupling the time-dependent stress-strain behavior of soils with a multiple-step numerical simulation method, the evolution of long-term settlement of tunnels in clays was studied and the following conclusions were obtained:

1. The Shanghai No.2 metro acts as a drain in the long-term. When the tunnel ling is permeable, surface settlement, ground loss and width of surface settlement troughs increase with time during the equilibrium period. In contrast, ground loss begins to decrease, after an initial increase, with time, under undrained condition.

Fig. 1

 1 Fig.1 Long-term transverse surface settlement troughs over (a) Grimsby Tunnel in silty clay in England (after O'Reilly et al., 1991), and (b) Shanghai No.2 Tunnel in soft clay (after[START_REF] Lee | Ground Response to the Construction of Shanghai Metro Tunnel-Line 2[END_REF]

  Fig.2 FEM mesh in numerical simulation

.

  The evolutions of long-term surface settlements over tunnels with the four different permeability cases are presented in Fig.3.

Fig. 5

 5 Fig.5 Width of long-term transverse surface settlement Fig. 6 Ground loss with time with varied troughs with time with different permeability permeability of tunnel

  Fig.7 demonstrates that: (1) the tunnel acts as a drain in the long-term. An encouraging agreement is found between the calculated and observed long-term surface settlement with 001 . 0 = v t k k . (2) The time-dependent properties of the soils around the tunnel enlarge the long-term settlements.

Fig. 7

 7 Fig.7 Comparison of calculated and in-situ measured surface settlements

Table 1

 1 Mechanical parameters of soils of Shanghai

	Depth (m)	0~2.40 2.40~4.0 4.0~8.0	8.0~17.7 17.7~28.7 28.7~31.7
	Young's modulus E (kPa) Poisson's ratio ν	2395 0.35	2032 0.35	1293 0.35	2710 0.35	3318 0.35	4216 0.35
	Water content W (%) Compression index λ	34.8 0.121	40.5 0.10	51.1 0.230	34.8 0.065	33.0 0.065	23.8 0.065
	Recompression index k	0.015	0.010	0.031	0.008	0.008	0.008
	Slope of C.S.L. M	1.32	1.14	1.01	1.20	1.47	1.20
	Permeability k (10 -9 m/s)	3.27	2.66	1.82	3.00	247	3.00
	Void ratio e	0.97	1.13	1.44	1.01	0.98	0.70
	Table 2 Tunnel construction parameters
	Advance rate of shield machine		Physical gap		Grouting pressure
	m/d		mm			MPa
	12		140			0.25

The long-term settlement of Shanghai No.2 was significant.

The evolution of long-term surface settlement troughs is dependent on the tunnel excavation methods. The shapes of surface settlement troughs do not change much in the long-term using the EPB tunnel excavation method.