# Solving dynamic contact problems with local refinement in space and time

Corinna I Hager, Patrice Hauret, Patrick Le Tallec, Barbara I Wohlmuth

► To cite this version:

Corinna I Hager, Patrice Hauret, Patrick Le Tallec, Barbara I Wohlmuth. Solving dynamic contact problems with local refinement in space and time. Computer Methods in Applied Mechanics and Engineering, 2012, 201, pp.25 - 41. 10.1016/j.cma.2011.09.006. hal-01393141

HAL Id: hal-01393141

https://hal.science/hal-01393141

Submitted on 6 Nov 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Solving dynamic contact problems with local refinement in space and time

Corinna Hager<sup>a</sup>, Patrice Hauret<sup>b</sup>, Patrick Le Tallec<sup>c</sup>, Barbara I. Wohlmuth<sup>d</sup>

<sup>a</sup> IANS, Universität Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany

<sup>b</sup> Mécanique Numérique, Manufacture Française des Pneumatiques Michelin, 63040 Clermont-Ferrand Cedex 9, France

<sup>c</sup> Département de Mécanique, Ecole Polytechnique, 91128 Palaiseau Cedex, France

<sup>d</sup> M2 – Zentrum Mathematik, Technische Universität München, Boltzmannstraße 3, 85748 Garching, Germany

Frictional dynamic contact problems with complex geometries are a challenging task from the computational as well as from the analytical point of view since they generally involve space and time multi scale aspects.

To be able to reduce the complexity of this kind of contact problem, we employ a non conforming domain decomposition method in space, consisting of a coarse global mesh not resolving the local structure and an overlapping fine patch for the contact computation. This leads to several benefits: First, we resolve the details of the surface only where it is needed, i.e., in the vicinity of the actual contact zone. Second, the subproblems can be discretized independently of each other which enables us to choose a much finer time scale on the contact zone than on the coarse domain. Here, we propose a set of interface conditions that yield optimal *a priori* error estimates on the fine meshed subdomain without any artificial dissipation. Further, we develop an efficient iterative solution scheme for the coupled problem that is robust with respect to jumps in the material parameters. Several complex numerical examples illustrate the performance of the new scheme.

**Keywords:** Domain decomposition; Frictional contact; Semismooth Newton; Multi-scale discretization; Mortar coupling; Time subcycling

## 1. Introduction

The numerical simulation of dynamic contact problems plays an important role in many applications in mechanics, like the forming of sheet metal, crash tests and other examples of structures under impact, fracture dynamics, or tire rolling. Especially the latter application is a challenging task from the point of view of simulation, as the problem usually features a complex three dimensional geometry, nonlinear elastic materials as well as dynamic effects. In addition, the contact zone is usually quite small compared to the size of the tire but needs to be resolved very accurately to get a good picture of the evolution of the contact pressures during rolling contact. More importantly, many other multiscale contact simulations require local refinement in the vicinity of the contact zone which is a *priori* unknown.

In order to be able to perform an accurate simulation of a car tire, there is a huge demand for a sound numerical scheme that combines a suitable multi scale discretization in space and in time with an efficient solution algorithm for the frictional contact

conditions. The algorithm has to be robust with respect to jumps in the material parameters as well as to provide an energy consistent description of the dynamics. The aim of this work is to design such an algorithm by combining several well established mathematical methods with some new approaches. All of them are described in the following.

The basis of the algorithm is provided by a decomposition of the original structure into several overlapping subdomains which have different mesh sizes. The global *d* dimensional structure, *d* ∈ {2, 3}, is discretized with a relatively coarse mesh that does not resolve the details along the contact boundary, whereas at the contact area, an overlapping patch with a fine triangulation is introduced. A two dimensional example of such a geometry is sketched in Fig. 1. In order to avoid an expensive volume coupling, the transfer between the subdomains is only performed at the inner (*d* 1) dimensional interface. Here, we employ the variationally consistent tent mortar method (see, e.g., [1,2]) with dual Lagrange multipliers [3] to enforce the weak continuity of the traces.

The subject of domain decomposition methods is already well established in the literature; we refer to [4, 6] and the references therein for an overview of the topic. The construction of domain decomposition schemes which are robust with respect to the mesh size as well as to jumps in the material parameters has been the topic of several papers (e.g., [7, 9]). In this work, we make use of the overlapping decomposition in order to obtain an iterative
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sketch of overlapping domain decomposition; left: coarse domain; middle: fine patch resolving the details at the boundary and right: combined geometry.</figcaption></figure>

solution scheme whose convergence rate is bounded indepen-
dently of the mesh size or the Lamé parameters in the subdomains.

The next important item is the treatment of the contact inequality constraints (see [10 14] and the references therein for an overview of the topic). These conditions are enforced in a variationally consistent way using dual Lagrange multipliers, allowing for the application of a primal dual active set strategy [15 17]. This scheme can be interpreted as a semismooth Newton method [15] applied to a set of nonlinear equations describing the contact conditions [18,19]. In combination with the iterative subdomain coupling, an inexact Newton scheme is obtained which still shows superlinear local convergence provided that appropriate stopping criteria for the inner iteration are satisfied [20 22].

The incorporation of inertia effects into the formulation makes the simulation of the nonlinear contact problem even more chal lenging. Standard time stepping algorithms like the trapezoidal rule generally lose their energy conservation property if applied to nonlinear problems. Possible remedies are presented in, e.g., [23 25] for nonlinear material laws and in [26,27] for contact. But even with the energy consistent formulation of [27], the com puted results for the contact stresses show spurious oscillations in time. This is avoided employing a local modification of the mass matrix along the potential contact boundary [28 30].

The last feature of the algorithm is the possibility to use differ-
ent time step sizes in the subdomains. Here, the main challenge
lies in the construction of the interface constraints. A suboptimal
choice can lead to numerical instability or to artificial dissipation
at the interface, even if the time integrators in the subdomains
are stable and energy conserving [31 33]. In [34,35], a time sub-
stepping scheme with linear interpolation of the velocity and the
Lagrange multipliers has been proposed which is stable but dissi-
pative at the interface. An improved energy conserving scheme
with linear interpolation of the multipliers has been analyzed in
[36,37]. However, both methods rely on the expensive exact solu-
tion of the resulting coupled system. In contrast, we present a
time discrete coupled system that uses different time step sizes
in each subdomain, conserves the energy over a coarse time step
and does not require the exact solution of the coupled problem.

At this point, one might ask whether the algorithm presented in
this work possibly is just a combination of known numerical tech-
niques. As a return, we recall that we aim for formulating an algo-
rithm that is able to solve challenging dynamic contact problems
with complex local geometries. This cannot be achieved without
building upon the knowledge and experience gained in the above
mentioned topics of numerical simulation, i.e., we combine domain
decomposition with mortar coupling, contact modeling via semi
smooth Newton methods and energy consistent time integration.
However, the algorithm presented in this work also contains
important new aspects. The main novelty is the consistent use of
the overlapping domain decomposition approach for both the dis-
cretization and the solver, leading to a two scale formulation both
in space and in time which has been implemented with contact
constraints and dynamic terms. This is complemented by an exten-
sive theoretical and numerical investigation of the properties of the
resulting method.

We now turn to the structure of the rest of this work. In Sections
2 4, we present and investigate the proposed algorithm for a linear

dynamic problem. Section 2 introduces the notation and the gov-
erning equations as well as the algebraic formulation of the result-
ing fully coupled system. In Section 3, we describe an efficient
iterative solution scheme based on overlapping domain decompo-
sition and analyse its convergence rate. Section 4 numerically con
firms the theoretical results by means of several tests.

In Section 5, we extend the domain decomposition approach to nonlinear problems, leading to an inexact semismooth Newton method. A special focus is put on the approximation of the fric tional contact conditions. Section 6 contains several numerical examples illustrating the efficiency and the robustness of the resulting iterative scheme.

The case of different time step sizes is analyzed in Section 7,
where we present an energy conserving formulation that can effi-
ciently be solved by the iterative scheme used before. The numer-
ical results in Section 8 illustrate that the scheme can considerably
decrease the number of global systems to be solved and still pro-
vide a good local accuracy. Section 9 concludes the work with a
short summary.

## 2. Problem formulation for the linear setting

This section contains the problem formulation as well as the ba-
sic notation for the rest of this work. In Section 2.1, the governing
equations for the linear problem are stated in their strong and
weak forms, whereas Section 2.2 introduces the spatial and tempo-
ral discretization, including some properties of trace spaces and
mortar operators which will be used in the sequel. In Section 2.3,
the algebraic Schur complement formulation is presented.

Since Sections 2.2 and 2.3 only recall known results on mortar methods and Schur complement formulations, they can possibly be skipped by readers familiar with these topics.

2.1. Problem statement

In the following, we consider an elastic body Ω ⊂ R $^{d}$, where
d ∈ { 2 , 3 } denotes the number of spatial dimensions. The polyhe-
dral boundary ∂ Ω is partitioned into two nonoverlapping parts
Γ$_{D}$ , Γ$_{N}$ with meas(Γ$_{D}$) > 0. On Γ$_{D}$, we assume homogeneous Dirich
let boundary conditions for simplicity, whereas a surface load de
noted by g$_{N}$ acts on the boundary Γ$_{N}$ where the unit outer
normal n is defined almost everywhere. Then, the strong form of
the dynamic linear elasticity problem on Ω subject to the volume
load l reads (see, e.g., [38])

$$
\rho \ddot{u} - \text{div}\boldsymbol{\sigma}(\mathbf{u}) = \mathbf{l} \quad \text{in } \Omega, \qquad (1a)
$$

$$
\mathbf{u} = \mathbf{0} \quad \text{on } \Gamma_D, \tag{1b}
$$

$$
\boldsymbol{\sigma}(\mathbf{u})\mathbf{n} = \mathbf{g}_N \quad \text{on } \Gamma_N. \tag{1c}
$$

Here, the Cauchy stress **σ** = **σ**(**u**) is given by **σ** : **C**<sup>el</sup>**ε**(**u**) with the linearized strain tensor **ε**(**u**): 1/2(∇**u** + ∇**u**<sup>T</sup>). The Hooke tensor **C**<sup>el</sup> **C**<sup>el</sup>(**x**) can be defined either in terms of the Lamé parameters **λ**, **μ** or with respect to the Young modulus **E** and the Poisson coeffi- cient **v**. Further, **ρ** ≥ 0 denotes the density of **Ω**.

We consider the case for which there exists a given subregion
ω ⊂ Ω, with possibly different material parameters, where a more
accurate local resolution is desired. In the situation of Fig. 1, the
<figure><img src="image_3.png" /><figcaption>Fig. 2. Subdomains and interfaces.</figcaption></figure>

domain shown on the left picture is $\Omega$, and the union of the patches depicted on the right picture corresponds to $\omega$. The interface $\partial\omega \cap \Omega$ is denoted by $\Gamma$ and the domain $\Omega \setminus \bar{\omega}$ by $\Xi$, as sketched on the left side of Fig. 2. For simplicity, we assume that $\bar{\Gamma}_D \cap \bar{\Gamma} = \emptyset$, that the material is compressible, and that the material parameters $E, v, \rho$ are piecewise constant on each of the subdomain mains. The corresponding values are denoted by $\varrho\Theta, E\Theta, v\Theta, \Theta \in \{\Xi, \omega\}$. A more general formulation also valid for incompressible material is investigated in [39].

In order to obtain the variational form of (1), we introduce the vector valued function spaces

$$ \mathbf{H}^k(\Theta) : [H^k(\Theta)]^d, \quad \mathbf{V}_\Theta : \{\mathbf{v} \in \mathbf{H}^1(\Theta) : \mathbf{v}|_{\Gamma_D \cap \partial\Theta} = \mathbf{0}\} $$

for $\Theta \in \{\Omega, \Xi, \omega\}$ and the (bi) linear forms

$$ m_{\Theta}(\mathbf{u}, \mathbf{w}) : \int_{\Theta} \varrho \mathbf{u} \cdot \mathbf{w} \, d\mathbf{x}, $$

$$ a_{\Theta}(\mathbf{u}, \mathbf{w}) : \int_{\Theta} \mathbb{C}^{\mathrm{el}} \varepsilon(\mathbf{u}) : \varepsilon(\mathbf{w}) d\mathbf{x}, $$

$$ f_{\Theta}(\mathbf{w}) : \int_{\Theta} \mathbf{l} \cdot \mathbf{w} d\mathbf{x} + \int_{\Gamma_N \cap \partial\Theta} \mathbf{g}_N \cdot \mathbf{w} ds. $$

Further, let $\mathbf{W}_\Gamma : \mathbf{H}^{1/2}(\Gamma)$ be the space of traces on $\Gamma$. Then, the weak form of (1) can be formulated as two separate problems on the subdomains $\Xi, \omega$, where the continuity of the displacements across $\Gamma$ is enforced weakly by means of a Lagrange multiplier $\zeta_\Gamma \in \mathbf{M}_\Gamma : \mathbf{W}_\Gamma$ which is the dual space of $\mathbf{W}_\Gamma$. This leads to the fol- lowing coupled problem: find $(\mathbf{u}|_{\Xi}, \mathbf{u}|_{\omega}, \zeta_\Gamma) \in \mathbf{V}_\Xi \times \mathbf{V}_\omega \times \mathbf{M}_\Gamma$ such that for all $t \in (0, T]$

$$ m_{\Xi}(\ddot{\mathbf{u}}, \mathbf{w}) + a_{\Xi}(\mathbf{u}, \mathbf{w}) \langle \mathbf{w}, \zeta_{\Gamma} \rangle_{\Gamma} f_{\Xi}(\mathbf{w}), \quad \mathbf{w} \in \mathbf{V}_{\Xi}, \quad (2a) $$

$$ m_{\omega}(\ddot{\mathbf{u}}, \mathbf{w}) + a_{\omega}(\mathbf{u}, \mathbf{w}) + \langle \mathbf{w}, \zeta_{\Gamma} \rangle_{\Gamma} f_{\omega}(\mathbf{w}), \quad \mathbf{w} \in \mathbf{V}_{\omega}, \quad (2b) $$

$$ \langle \mathbf{u}|_{\omega}, \boldsymbol{\mu}_{\Gamma} \rangle_{\Gamma} \quad \langle \mathbf{u}|_{\Xi}, \boldsymbol{\mu}_{\Gamma} \rangle_{\Gamma} \quad \mathbf{0}, \quad \boldsymbol{\mu}_{\Gamma} \in \mathbf{M}_{\Gamma}, \quad (2c) $$

plus appropriate initial conditions $\mathbf{u}|_{t=0} \quad \mathbf{u}_0$, $\dot{\mathbf{u}}|_{t=0} \quad \mathbf{v}_0$. Above, we have used the $L_2$ duality pairing $\langle \cdot, \cdot \rangle_\Gamma$ on $\Gamma$. For the quasi static case, the well posedness of (2) follows from the condition meas $(\Gamma_D) > 0$ and Korn's inequality [40,41], whereas for the dynamic case, we need that both initial values $\mathbf{u}_0, \mathbf{v}_0$ satisfy condition (2c) [42].

## 2.2. Discretization

For the spatial discretization of (2), we first triangulate the global domain $\Omega$ in terms of a quasi uniform mesh $\mathcal{T}^H$ of simplicial or quadrilateral/hexahedral elements of size $H$. We assume that the Dirichlet boundary $\Gamma_D$ as well as the interface $\Gamma$ are re-solved by this triangulation. On $\mathcal{T}^H$, we consider lowest order con- forming finite element basis functions $\phi_p^H$, $p \in \mathcal{N}^H$, with $\mathcal{N}_\Gamma^H \subset \mathcal{N}^H$ denoting the subset of vertices of $\mathcal{T}^H$ on the interface $\Gamma$. On the patch $\omega$, we introduce a second regular and quasi uniform discretization $\mathcal{T}^h$ with elements of size $h < H$, the basis functions $\phi_p^h$, $p \in \mathcal{N}^h$, and the interface nodes $\mathcal{N}_\Gamma^h \subset \mathcal{N}^h$. We emphasize that the grids need not be nested on $\omega$.

With these triangulations, we define the following vector valued finite element spaces associated with the coarse and the fine discretization:

$$ \mathbf{V}_{\Xi}^{H} : \operatorname{span}\{\phi_{p}^{H}|_{\Xi}\}_{p \in \mathcal{N}_{\Xi}^{H}} \subset \mathbf{V}_{\Xi}, \quad \mathbf{V}_{\omega}^{H} : \operatorname{span}\{\phi_{p}^{H}|_{\omega}\}_{p \in \mathcal{N}_{\omega}^{H}} \subset \mathbf{V}_{\omega}, \quad (3a) $$

$$ \mathbf{V}^H : \{ (\mathbf{v}^H, \mathbf{w}^H) \in \mathbf{V}_{\Xi}^H \oplus \mathbf{V}_{\omega}^H : \mathbf{v}^H|_{\Gamma} \quad \mathbf{w}^H|_{\Gamma} \} \subset \mathbf{V}(\Omega), \quad (3b) $$

$$ \mathbf{V}^h \quad \mathbf{V}_\omega^h : \operatorname{span}\{\phi_p^h\}_{p \in \mathcal{N}_\omega^h} \subset \mathbf{V}_\omega. \quad (3c) $$

On the interface $\Gamma$, we define the discrete trace spaces $\mathbf{W}_\Gamma^m : \mathbf{W}^m|_{\Gamma} \subset \mathbf{W}_\Gamma$ as well as the Lagrange multiplier space $\mathbf{M}_\Gamma^m \subset \mathbf{M}_\Gamma$ spanned by the basis functions $\{\psi_p^m\}_{p \in \mathcal{N}_\Gamma^m}$, $m \in \{h, H\}$. There are several suitable possibilities to define these functions (see [43] and the references therein for more details); in the rest of this work, we use the so called dual basis functions $\{\psi_p^m\}_{p \in \mathcal{N}_\Gamma^m}$ which are piecewise (bi) linear, discontinuous and satisfy the locality requirement $\mathrm{supp}(\psi_p^m) \ \mathrm{supp}(\phi_p^m|_{\Gamma})$ as well as the biorthogonality property

$$ \int_{\Gamma} \psi_p^m \phi_q^m \, ds \quad \delta_{pq} \int_{\Gamma} \phi_q^m \, ds, \quad p, q \in \mathcal{N}_{\Gamma}^m, \quad m \in \{h, H\}. \quad (4) $$

As shown in [44], these basis functions can efficiently be constructed even for general grids in 3D. We remark that the above definition of the multiplier spaces yields $\dim(\mathbf{M}_\Gamma^m) = \dim(\mathbf{W}_\Gamma^m)$, i.e., we do not modify the basis functions of $\mathbf{M}_\Gamma^m$ near the locations where $\Gamma$ is nondifferentiable.

In the rest of this subsection, we formulate the discrete algebraic version of (2). For this, we introduce the mortar projection onto $\mathbf{W}_\Gamma^m$, $m \in \{h, H\}$, given by

$$ P^m : \mathbf{W}_\Gamma \to \mathbf{W}_\Gamma^m : \langle P^m \mathbf{w}_\Gamma, \boldsymbol{\mu}_\Gamma^m \rangle_\Gamma \quad \langle \mathbf{w}_\Gamma, \boldsymbol{\mu}_\Gamma^m \rangle_\Gamma, \quad \boldsymbol{\mu}_\Gamma^m \in \mathbf{M}_\Gamma^m. \quad (5) $$

As shown in [43], the mortar operator (5) is uniformly continuous in $h$ or $H$ with respect to the $L_2(\Gamma)$, the $\mathbf{H}^1(\Gamma)$ and the $\mathbf{H}^{1/2}(\Gamma)$ norm.

We define the matrices $D_\Gamma^{lm} \in \mathbb{R}^{d|\mathcal{N}_\Gamma^l| \times d|\mathcal{N}_\Gamma^m|}$ with $l, m \in \{h, H\}$ which are composed of the $d \times d$ submatrices

$$ (D_\Gamma^{lm})_{pq} \quad \mathrm{Id} \cdot \int_{\Gamma} \psi_p^l \phi_q^m ds, \quad p \in \mathcal{N}_\Gamma^l, \quad q \in \mathcal{N}_\Gamma^m \quad (6) $$

with $l$ denoting the $d \times d$ identity matrix. With this, the algebraic representation of the operator $P^l|_{\mathbf{w}_\Gamma^m} : \mathbf{W}_\Gamma^m \to \mathbf{W}_\Gamma^l$ is given by

$$ \Pi^{lm} \mathbf{w}_\Gamma^m \quad (D_\Gamma^l)^{1} D_\Gamma^{lm} \mathbf{w}_\Gamma^m. $$

**Remark 1.** Due to the biorthogonality property (4), the matrices $D_\Gamma^l$ in (6) are diagonal and can easily be inverted.

Above and in the following, the discrete functions and the cor- responding coefficient vectors are denoted with the same symbol for ease of notation. The extension of the matrices $D_\Gamma^{lm}$, $l \in \{h, H\}$, by zero to the vector space $\mathbf{V}^h$ is named $D_\Gamma^{lm} \in \mathbb{R}^{d|\mathcal{N}_\Gamma^l| \times d|\mathcal{N}_\Gamma^h|}$; similarly, the extension of $D_\Gamma^{lm}$ to the coarse upper space $\mathbf{V}_{\Xi}^H$ is denoted by $D^{lm}$.
Using this space discretization, standard matrix notation for the mass and stiffness matrices as well as the fine Lagrange multiplier space $M_I^h$ for the continuity constraints, the spatially discrete version of (2) reads

$$M_{\Xi}^{H} \dot{\mathbf{u}}^{H} + A_{\Xi}^{H} \mathbf{u}^{H} \quad (D^{hh})^{T} \xi_{I}^{h} \quad \mathbf{f}^{H}, \qquad (7a)$$

$$M_{\omega}^{h} \dot{\mathbf{u}}^{h} + A_{\omega}^{h} \mathbf{u}^{h} + (D^{hh})^{T} \xi_{I}^{h} \mathbf{f}^{h}, \qquad (7b)$$

$$D^{hh} \mathbf{u}^{h} \quad D^{hh} \mathbf{u}^{h} \quad \mathbf{0}. \qquad (7c)$$

together with suitable initial conditions

$$(\mathbf{u}^H, \mathbf{u}^H)|_t = \mathbf{0} \quad (\mathbf{u}_0^H, \mathbf{u}_0^H), \quad (\dot{\mathbf{u}}^H, \dot{\mathbf{u}}^H)|_t = \mathbf{0} \quad (\mathbf{v}_0^H, \mathbf{v}_0^H). \qquad (8)$$

**Remark 2.** For the case $W_I^h \subset W_I^h$, the definition of $M_I^h$ and (7c) imply that the finite element solution is continuous on the whole domain $\Omega$, and (7) results in a standard conforming scheme.

Next, we discretize (7) in time by partitioning the time interval $(0,T)$ into equidistant time steps $\Delta t$ with discrete times $t_j = j\Delta t$ and introducing the notation

$$\partial_{\Delta t} \mathbf{w}_j : \frac{1}{\Delta t} (\mathbf{w}_j - \mathbf{w}_j \mathbf{1}), \quad \mathbf{w}_j \mathbf{1} / 2 : \frac{1}{2} (\mathbf{w}_j + \mathbf{w}_j \mathbf{1}).$$

Using the implicit trapezoidal rule, i.e., the Newmark scheme with $\gamma = \frac{1}{2}$, $\beta = \frac{1}{4}$ [45], for the time stepping, we obtain the following problem: find sequences of vectors $(\mathbf{u}_j^h, \mathbf{u}_j^h \mathbf{1}_0, (\mathbf{v}_j^H, \mathbf{v}_j^H \mathbf{1}_0, (\xi_{Ij}^H \mathbf{1}_0 \mathbf{1}^T))$ satisfying the initial conditions (8) as well as

$$M_{\Xi}^{H} \partial_{\Delta t} \mathbf{v}_{j}^{H} + A_{\Xi}^{H} \mathbf{u}_{j}^{H} \mathbf{1}_{2} \quad (D^{hh})^{T} \xi_{Ij}^{h} \quad \mathbf{f}_{j}^{H} \mathbf{1}_{2} \quad \mathbf{0}, \qquad (9a)$$

$$\mathbf{v}_{j}^{H} \mathbf{1}_{2} \quad \partial_{\Delta t} \mathbf{u}_{j}^{H} \mathbf{0}, \qquad (9b)$$

$$M_{\omega}^{h} \partial_{\Delta t} \mathbf{v}_{j}^{h} + A_{\omega}^{h} \mathbf{u}_{j}^{h} \mathbf{1}_{2} + (D^{hh})^{T} \xi_{Ij}^{h} \quad \mathbf{f}_{j}^{h} \mathbf{1}_{2} \quad \mathbf{0}, \qquad (9c)$$

$$\mathbf{v}_{j}^{h} \mathbf{1}_{2} \quad \partial_{\Delta t} \mathbf{u}_{j}^{h} \mathbf{0}, \qquad (9d)$$

$$D^{hh} \mathbf{u}_{j}^{h} \quad D^{hh} \mathbf{u}_{j}^{h} \quad \mathbf{0}. \qquad (9e)$$

The discrete energy at time $t_j$ is the sum of the contributions from the subdomains $\Xi, \omega$ given by

$$\mathbb{E}_j \quad \mathbb{E}_{\Xi j}^{H} + \mathbb{E}_{\omega j}^{H} : \quad \left( \frac{1}{2} (\mathbf{v}_j^H)^T M_{\Xi}^H \mathbf{v}_j^H + \frac{1}{2} (\mathbf{u}_j^H)^T A_{\Xi}^H \mathbf{u}_j^H \quad (\mathbf{f}_j^H)^T \mathbf{u}_j^H \right) \\ + \left( \frac{1}{2} (\mathbf{v}_j^h)^T M_{\omega}^h \mathbf{v}_j^h + \frac{1}{2} (\mathbf{u}_j^h)^T A_{\omega}^h \mathbf{u}_j^h \quad (\mathbf{f}_j^h)^T \mathbf{u}_j^h \right). \qquad (10)$$

For time independent outer loads, the energy is conserved for the discrete system (9), i.e., $\mathbb{E}_j \mathbb{E}_j \mathbf{1}$ holds, provided that (9e) is valid for both time steps $t_j$ and $t_j \mathbf{1}$. This can easily be verified by multi plying (9ac) with the mean velocities $\mathbf{v}_j^m \mathbf{1}_{2}$ as test functions, adding the results and using (9bde).

Using (9bd) to eliminate the velocities from (9ac) and defining

$$K_{\omega}^{h} : \frac{2}{\Delta t^2} M_{\omega}^{h} + \frac{1}{2} A_{\omega}^{h}, \qquad (11a)$$

$$\mathcal{Q}_{j}^{h} : \mathbf{f}_{j}^{h} \mathbf{1}_{2} + \frac{2}{\Delta t} M_{\omega}^{h} \mathbf{v}_{j}^{h} \mathbf{1}_{2} \quad A_{\omega}^{h} \mathbf{u}_{j}^{h} \mathbf{1}_{2} + K_{\omega}^{h} \mathbf{u}_{j}^{h} \mathbf{1}_{2}, \qquad (11b)$$

as well as analogous expressions for the quantities on $\Xi$, the system to be solved for the displacement at time $t_j$ becomes

$$\left( \begin{array}{ccc} K_{\Xi}^{H} & \mathbf{0} & (D^{hh})^{T} \\ \mathbf{0} & K_{\omega}^{h} & (D^{hh})^{T} \\ D^{hh} & D^{hh} & \mathbf{0} \end{array} \right) \left( \begin{array}{ccc} \mathbf{u}_{j}^{H} \\ \mathbf{u}_{j}^{h} \\ \xi_{Ij}^{H} \end{array} \right) \left( \begin{array}{ccc} \mathcal{Q}_{j}^{H} \mathbf{1}_{1} \\ \mathcal{Q}_{j}^{h} \mathbf{1}_{1} \\ \mathbf{0} \end{array} \right). \qquad (12)$$

## 2.3. Schur complement formulation

The system (12) can equivalently be expressed in terms of the unknowns on the interface $\Gamma$ only. For this, we partition the vec tors $(\mathbf{u}^H, \mathbf{u}^H)$ and the corresponding matrices $K_{\Xi}^H, K_{\omega}^H$ into compo

nents associated with the interface nodes and the inner degrees of freedom:

$$\begin{align*} \mathbf{u}^H \quad \begin{pmatrix} \mathbf{u}_\Xi^H \\ \mathbf{u}_I^H \end{pmatrix}, \quad K_\Xi^H \quad \begin{pmatrix} K_{\Xi\Xi}^H & K_{\Xi I}^H \\ K_{I\Xi}^H & K_{II}^H \end{pmatrix}, \\ \mathbf{u}^H \quad \begin{pmatrix} \mathbf{u}_I^h \\ \mathbf{u}_\omega^h \end{pmatrix}, \quad K_\omega^h \quad \begin{pmatrix} K_{I\Gamma}^h & K_{\Gamma\omega}^h \\ K_{\omega\Gamma}^h & K_{\omega\omega}^h \end{pmatrix}. \end{align*}$$

Performing static condensation of the inner variables and introducing the symmetric positive semi definite Schur complement matrices $S_{\Xi}^H, S_{\omega}^H$ given by

$$S_{\Xi}^{H} : \quad K_{I\Gamma}^{H} \quad K_{\Gamma\Xi}^{H} (K_{\Xi\Xi}^{H}) \quad K_{\Xi\Gamma}^{H} \in \mathbb{R}^{d|\mathcal{N}_I^H| \times d|\mathcal{N}_I^H|}, \qquad (13a)$$

$$S_{\omega}^{H} : \quad K_{I\Gamma}^{h} \quad K_{\Gamma\omega}^{h} (K_{\omega\omega}^{h}) \quad K_{\omega\Gamma}^{h} \in \mathbb{R}^{d|\mathcal{N}_I^h| \times d|\mathcal{N}_I^h|}, \qquad (13b)$$

we can rewrite (9) as

$$\left( \begin{array}{ccc} S_{\Xi}^{H} & \mathbf{0} & (D_{I}^{hH})^{T} \\ \mathbf{0} & S_{\omega}^{h} & (D_{I}^{hh})^{T} \\ D_{I}^{hh} & D_{I}^{hh} & \mathbf{0} \end{array} \right) \left( \begin{array}{ccc} \mathbf{u}_{Ij}^{H} \\ \mathbf{u}_{Ij}^{h} \\ \xi_{Ij}^{H} \end{array} \right) \left( \begin{array}{ccc} \mathcal{Q}_{\Xi(j)1}^{H} \\ \mathcal{Q}_{\omega(j)1}^{h} \\ \mathbf{0} \end{array} \right), \qquad (14)$$

where the right hand side of (14) is defined as

$$\mathcal{Q}_{\Xi(j)1}^{H} : \quad \mathcal{Q}_{I(j)1}^{H} \quad K_{I\Xi}^{H} (K_{\Xi\Xi}^{H}) \quad 1 \mathcal{Q}_{\Xi(j)1}^{H}, \qquad (15a)$$

$$\mathcal{Q}_{\omega(j)1}^{h} : \quad \mathcal{Q}_{I(j)1}^{h} \quad K_{I\omega}^{h} (K_{\omega\omega}^{h}) \quad 1 \mathcal{Q}_{\omega(j)1}^{h}, \qquad (15b)$$

In many applications, the exact solution of the fully coupled problem (14) is very expensive. Instead, we consider an iterative solution scheme that incorporates a coarse mortar problem on $\mathbf{v}^H$. In Section 3, the corresponding algorithm is derived and theoretically analyzed. Numerical results are presented in Section 4.

## 3. Iterative coupling algorithm

In this section, we derive and investigate an iterative solution scheme for the coupled system (14), where each step consists of the solution of a coarse problem on the full domain $\Omega$ and a subse quent solution of a local fine grid problem on $\omega$. In Subsection 3.1, the algorithm is stated. The error propagation and the convergence rate of the iterative scheme are analyzed in Sections 3.2 and 3.3, respectively, followed by a short discussion on how its algebraic error can be measured.

### 3.1. Derivation

In order to formulate an iterative solution algorithm for the mortar system (14), we construct a coarse grid version of its second equation. For this, we introduce a Schur complement matrix

$$S_{\omega}^{H} : \quad K_{I\Gamma}^{H} \quad K_{\Gamma\omega}^{H} (K_{\omega\omega}^{H}) \quad K_{\omega\Gamma}^{H} \in \mathbb{R}^{d|\mathcal{N}_I^H| \times d|\mathcal{N}_I^H|} \qquad (16)$$

and a right hand side vector $\tilde{\mathcal{Q}}_{\omega j}^{H} \in \mathbb{R}^{d|\mathcal{N}_I^H|}$ which are defined similar as in (13b), (15b) but assembled with respect to the coarse grid $T_{I}^{H}|_{\omega}$ instead of $T^{-h}$. Further, we add an auxiliary coefficient vector $\mu_{\omega}^{H} \in \mathbb{R}^{d|\mathcal{N}_I^H|}$ which can be interpreted as a coarse Lagrange multiplier in $M_I^H$. These additional definitions allow to formulate an iterative scheme for the unknowns

$$\hat{\mathbf{z}} : \left( (\mathbf{u}_{Ij}^{H})^{T} \quad (\mu_{Ij}^{H})^{T} \quad (\mathbf{u}_{Ij}^{h})^{T} \quad (\xi_{Ij}^{h})^{T} \right)^{T}. \qquad (17)$$

The distinctive feature of the algorithm is that it contains a global coarse problem on $\Omega$ subject to an interface load on $\Gamma$, whose solution defines a coarse grid approximation of the fine grid quantities in $\omega$. The details of the scheme are summarized in Algorithm 1.

**Remark 3.** After one iteration of Algorithm 1, i.e., for $l \ge 1$, the residuals (20bcd) vanish.
**Remark 4.** Definition (16) ensures that the Dirichlet boundary conditions on $\Gamma_D \cap \partial\omega$ are respected in both $S_{\omega}^H$ and $S_{\omega}^h$ which is important for the performance of the iterative method (cf. [46] for a related algorithm). We come back to this point in Section 5, where we discuss contact problems.

**Algorithm 1.** Two way coupling scheme with augmented coarse grid problem

Starting from some initial guess $\hat{\mathbf{z}}^{(0)}$, compute sequentially for $l=0,1,\dots$

(i) Solve problem on global coarse space $\mathbf{V}^H$ with interface load on $\Gamma$ inherited from fine computation on $\omega$:

$$ (S_{\bar{\omega}}^H + S_{\omega}^H) \begin{pmatrix} \mathbf{0} \\ (D_{\Gamma}^{HH})^T \end{pmatrix} \begin{pmatrix} \delta \mathbf{u}_{\Gamma j}^{H,(l)} \\ \delta \boldsymbol{\mu}_{\Gamma j}^{H,(l)} \end{pmatrix} + \mathbf{r}_{\bar{\omega}}^{H,(l)} + \mathbf{r}_{\omega}^{H,(l)}. \quad (18) $$

(ii) Solve problem on local fine space $\mathbf{V}^h$ with weakly imposed trace on $\Gamma$ inherited from coarse computation on $\Omega$:

$$ \begin{pmatrix} S_{\omega}^{h} & (D_{\Gamma}^{hh})^{T} \\ D_{\Gamma}^{hh} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \delta \mathbf{u}_{\Gamma j}^{h,(l)} \\ \delta \boldsymbol{\zeta}_{\Gamma j}^{h,(l)} \end{pmatrix} + \begin{pmatrix} \mathbf{0} \\ D_{\Gamma}^{hh} \delta \mathbf{u}_{\Gamma j}^{H,(l)} \end{pmatrix}. \quad (19) $$

(iii) Update the solution vector:

$$ \hat{\mathbf{z}}^{(l+1)} : \hat{\mathbf{z}}^{(l)} + \delta\hat{\mathbf{z}}^{(l)}. $$

The residuals of (18), (19) are given by

$$ \mathbf{r}_{\bar{\omega}}^{H,(l)} \quad \mathbf{q}_{\bar{\omega}j}^{H} \quad 1 \quad S_{\bar{\omega}}^{H} \mathbf{u}_{\Gamma j}^{H,(l)} + (D_{\Gamma}^{hh})^{T} \boldsymbol{\zeta}_{\Gamma j}^{H,(l)}, \quad (20a) $$

$$ \mathbf{r}_{\omega}^{H,(l)} \quad \mathbf{q}_{\omega j}^{H} \quad 1 \quad S_{\omega}^{H} \mathbf{u}_{\Gamma j}^{H,(l)} (D_{\Gamma}^{HH})^{T} \boldsymbol{\mu}_{\Gamma j}^{H,(l)}, \quad (20b) $$

$$ \mathbf{r}_{\omega}^{h,(l)} \quad \mathbf{q}_{\omega j}^{h} \quad 1 \quad S_{\omega}^{h} \mathbf{u}_{\Gamma j}^{h,(l)} (D_{\Gamma}^{hh})^{T} \boldsymbol{\zeta}_{\Gamma j}^{h,(l)}, \quad (20c) $$

$$ \mathbf{v}_{\Gamma}^{h,(l)} \quad D_{\Gamma}^{hh} \mathbf{u}_{\Gamma j}^{h,(l)} \quad D_{\Gamma}^{hh} \mathbf{u}_{\Gamma j}^{h,(l)}. \quad (20d) $$

## 3.2. Error propagation

In order to analyse the error propagation of Algorithm 1, we reformulate the mortar system (14) as an equation with the only unknown $\mathbf{u}_{\Gamma j}^H \in \mathbb{R}^{d|\Lambda_{\Gamma}^H| \times d|\Lambda_{\Gamma}^H|}$ with

$$ S_{\omega}^{HH} : \begin{pmatrix} \mathbf{0} & (D_{\Gamma}^{hh})^{T} \\ D_{\Gamma}^{hh} & \mathbf{0} \end{pmatrix} \begin{pmatrix} S_{\omega}^{H} & (D_{\Gamma}^{hh})^{T} \\ D_{\Gamma}^{hh} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \mathbf{0} \\ D_{\Gamma}^{hh} \end{pmatrix} (\Pi^{hh})^{T} S_{\omega}^{H} \Pi^{hh}, \quad (21) $$

(14) can be rewritten as

$$ (S_{\bar{\omega}}^{H} + S_{\omega}^{HH}) \mathbf{u}_{\Gamma j}^{H} \quad \mathbf{g}_{\Gamma j}^{H} : (\Pi^{hh})^{T} \boldsymbol{\zeta}_{\omega j}^{h} + \boldsymbol{\zeta}_{\bar{\omega}j}^{H}, \quad (22) $$

The following lemma shows that Algorithm 1 can be reformulated as a fixed point iteration on $\mathbf{u}_{\Gamma j}^H$:

**Lemma 1.** Let $\mathbf{u}_{\Gamma j}^H$ denote part of the exact solution of the mortar system (14), and let $\mathbf{u}_{\Gamma j}^{H,(l)}$, $l=0,1,\dots$, be the sequence of trace vectors obtained from Algorithm 1. For $l \ge 1$, the error $\mathbf{e}_{\Gamma}^{H,(l)} : (\mathbf{u}_{\Gamma j}^{H,(l)} \quad \mathbf{u}_{\Gamma j}^{H})$ satisfies the relation

$$ \mathbf{e}_{\Gamma}^{H,(l+1)} \quad (\mathrm{Id} \quad (S_{\bar{\omega}}^{H} + S_{\omega}^{H})^{1} (S_{\bar{\omega}}^{H} + S_{\omega}^{HH})) \mathbf{e}_{\Gamma}^{H,(l)}. \quad (23) $$

**Proof.** The result follows from the definition of the iterative scheme by substitution. Details can be found in [39]. □

**Remark 5.** If the coarse and the fine grid on $\omega$ coincide, i.e., $\mathbf{V}^h \quad \mathbf{V}_{\omega}^H$ holds with the definitions from (3), then we have $S_{\omega}^{HH} \quad S_{\omega}^H$, and Lemma 1 implies that the exact coarse trace $\mathbf{u}_{\Gamma j}^H$ is obtained after at most two steps. If the starting vector $\hat{\mathbf{z}}^{(0)}$ satisfies

$$ (D_{\Gamma}^{HH})^{T} (\boldsymbol{\mu}_{\Gamma j}^{H,(0)} \quad \boldsymbol{\mu}_{\Gamma j}^{H}) \quad (D_{\Gamma}^{hh})^{T} (\boldsymbol{\lambda}_{\Gamma j}^{h,(0)} \quad \boldsymbol{\lambda}_{\Gamma j}^{h}) \quad \mathbf{0}, $$

then the exact solution is already obtained after the first iteration.

A damped version of (23) with a suitable damping parameter $\alpha_l > 0$ yields the error propagation

$$ \mathbf{e}_{\Gamma}^{H,(l+1)} \quad (\mathrm{Id} \quad \alpha_l (S_{\bar{\omega}}^{H} + S_{\omega}^{H})^{1} (S_{\bar{\omega}}^{H} + S_{\omega}^{HH})^{T}) \mathbf{e}_{\Gamma}^{H,(l)}. \quad (24) $$

If $\alpha_l = \alpha$ does not depend on $l$, (24) is equivalent to a preconditioned Richardson iteration for the Schur complement system (22). But in general, it is more effective to choose $\alpha_l$ in each step according to a conjugate gradient algorithm (see, e.g., [6,7]).

**Remark 6.** From (23), one can see that the difference between Algorithm 1 and the classical Dirichlet Neumann coupling with $\Xi$ as the Neumann subdomain is the additional term of $S_{\omega}^H$ in the factor $(S_{\bar{\omega}}^H + S_{\omega}^H)^{-1}$. The benefit of this term can be seen in the next subsection, where we show that Algorithm 1 is robust with respect to jumps in the material parameters in contrast to the Dirichlet Neumann iteration whose convergence rate degrades if $E_{\omega} \gg E_{\Xi}$ or $\Omega_{\omega} \gg \Xi_{\omega}$, i.e., if the local details are much more stiff than the rest of the domain (cf. [7,39]).

## 3.3. Condition number analysis

The convergence properties of the iteration (24) or its conjugate gradient version are determined by the condition number of the iteration matrix in (23). To this extend, we investigate the spectral equivalence of the matrices $S_{\omega}^H$ and $S_{\omega}^{HH}$ in the following theorem proved in Appendix A:

**Theorem 2.** Assume that the trace spaces $\mathbf{W}_{\Gamma}^H \subset \mathbf{W}_{\Gamma}^h$ are nested and that there exists a constant $C_{CFL}$ with

$$ \frac{\Omega_{\omega}}{\Delta t^2} \le C_{CFL} \frac{E_{\omega}}{H^2}. \quad (25) $$

Then, there exist constants $c^*$, $C^*$ independent of the diameter of $\omega$, $h$, $H$, $\Delta t$ and the values $E_{\omega}$, $\Omega_{\omega}$, such that the following estimates are satisfied for any $\mathbf{w}_{\Gamma}^H \in \mathbf{W}_{\Gamma}^H$:

$$ c^* (1 + C_{CFL})^{1} (\mathbf{w}_{\Gamma}^{H}, S_{\omega}^{H} \mathbf{w}_{\Gamma}^{H}) \le (\mathbf{w}_{\Gamma}^{H}, S_{\omega}^{HH} \mathbf{w}_{\Gamma}^{H}) \le C^{*} (\mathbf{w}_{\Gamma}^{H}, S_{\omega}^{H} \mathbf{w}_{\Gamma}^{H}). \quad (26) $$

If the finite element spaces (3) satisfy $\mathbf{V}_{\omega}^H \subset \mathbf{V}^h$, then the constant $C^*$ in (26) can be replaced by one.

**Remark 7.** The assumption $\mathbf{W}_{\Gamma}^H \subset \mathbf{W}_{\Gamma}^h$ is not necessary but the proof becomes more technical. We refer to [39] for details.

Thus, if the assumptions of Theorem 2 are satisfied, the condition number of the iteration matrix in (23) is bounded by

$$ \kappa \left( (S_{\bar{\omega}}^{H} + S_{\omega}^{H})^{1} (S_{\bar{\omega}}^{H} + S_{\omega}^{HH}) \right) \le \frac{\max(1, C^*)}{\min(1, c^* (1 + C_{CFL})^{1})} $$

with the constants $c^*$, $C^*$ from Theorem 2. Furthermore, the iterates of (24) converge to the solution $\mathbf{u}_{\Gamma j}^H$ of (22) if the damping parameter $\alpha_l$ is chosen within the set $0 < \alpha_{\min} \le \alpha_l < \frac{2}{C}$ for some fixed value of $\alpha_{\min}$ (see, e.g., [6]). Hence, if the finite element spaces $\mathbf{V}_{\omega}^H \subset \mathbf{V}^h$ are nested, Theorem 2 implies the convergence of (24) for $\alpha_l = 1$, i.e., the convergence of Algorithm 1.
**Remark 8.** The value of $C_{CFL}$ is related to the Courant Friedrichs Levy condition often used to analyse the convergence of explicit time integration schemes. In general, the CFL condition gives a mesh dependent upper bound on the time step size of the form $\Delta t^2 \le Ch^2$, where C is some constant depending on the maximum value of the ratio $E/\rho$. Implicit integration schemes like the one used in this work are generally employed for problems where larger time steps are used than permitted by the CFL condition. Hence, (25) is a reasonable assumption for the problems we are interested in.

However, if Assumption (25) is only satisfied with a large constant $C_{CFL}$, the lower bound in (26) degenerates. In this case, the proof in Appendix A can be modified, leading (neglecting logarithmic terms) to a lower bound of

$$c^* \left(1 + \frac{H}{h}\right)^{1} \left(\mathbf{w}_T^H, S_\omega^H \mathbf{w}_T^H\right) \le \left(\mathbf{w}_T^H, S_\omega^{HhH} \mathbf{w}_T^H\right)$$

instead of (26). The resulting dependence of the convergence rate of Algorithm 1 on the ratio $H/h$ can be observed from the numerical results in Section 4.3; furthermore, it is in agreement with the theoretical estimates of the spectrum of discrete Dirichlet to Neumann operators recently presented in [47].

### 3.4. Stopping criteria

In this subsection, we state a measure of the algebraic error introduced by solving (14) by means of Algorithm 1.

As the residuals (20bcd) vanish for $l \ge 1$, the only nonzero component of the residual vector is $\mathbf{r}_\Xi^{H,(l)}$ given in (20a). As the energy norm of this residual given by

$$\left( (S_{\Xi}^{H} + S_{\omega}^{HhH}) \mathbf{1}_{\Xi}^{H,(l)}, \mathbf{r}_{\Xi}^{H,(l)} \right) \qquad (27)$$

is too expensive to compute, it is approximated by the value

$$\left( (S_{\Xi}^{H} + S_{\omega}^{H}) \mathbf{1}_{\Xi}^{H,(l)}, \mathbf{r}_{\Xi}^{H,(l)}, \delta \mathbf{u}_{ij}^{H,(l)}, \mathbf{r}_{\Xi}^{H,(l)} \right),$$

where the last equality follows from (18) for $l \ge 1$. Thus, we propose to use the following relative algebraic error estimator for $l \ge 1$:

<figure><img src="image_4.png" /><figcaption>Fig. 3. Linear example; left: grid for $L=2$; middle and right: effective stress $\sigma_{eff}$ and pressure for the static case.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. True $e_{alg}^{(l)}$ (solid line) and estimated $n_{alg}^{(l)}$ (dashed line) relative algebraic error for the static case with $L \in \{1,2,3,4\}$ with respect to $l$; left: $E_\omega = E_\Xi$ and right: $E_\omega = 10^5 E_\Xi$.</figcaption></figure>

$$\left(\eta_{\text{alg}}^{(l)}\right)^2 = \frac{\left(\delta \mathbf{u}_{ij}^{H,(l)}, \mathbf{r}_{\Xi}^{H,(l)}\right)}{\left(\mathbf{u}_{ij}^{H,(l)}, (S_{\Xi}^{H} + S_{\omega}^{H}) \mathbf{u}_{ij}^{H,(l)}\right)}. \qquad (28)$$

As (28) is a norm of the error between the exact mortar solution and the approximate solution by means of Algorithm 1, it can be used to define a stopping criterion for the iterative process (see also Section 6).

**Remark 9.** In order to compute the denominator of (28) without solving a Schur complement problem in each step, one can use the relation

$$\left(S_{\Xi}^{H} + S_{\omega}^{H}\right) \mathbf{u}_{ij}^{H,(l)} \quad \left(S_{\Xi}^{H} + S_{\omega}^{H}\right) \mathbf{u}_{ij}^{H,(0)} + \mathbf{r}_{\omega}^{H,(0)} + \sum_{k=0}^{l-1} \mathbf{r}_{\Xi}^{H,(k)}. $$

### 4. Numerical results for the linear setting

#### 4.1. Geometry and parameters

For these first numerical tests, we consider the domain $\Omega = (0,2) \times (0,1)$ which is split into the patch $\omega = (0.5,1.5) \times (0,0.5)$ and the upper domain $\Xi = \Omega \setminus \bar{\omega}$. Both subdomains are initially discretized by quadrilaterals of size $H = h = 0.125$; afterwards, we perform $L$ additional refinements of the fine grid on $\omega$. The resulting grid for $L=2$ is depicted on the left side of Fig. 3. If not stated otherwise, we iteratively solve the equations of linear elasticity (1) with $l=0$ and do not apply damping, i.e., $\alpha = 1$. The coarse matrix $S_\omega^H$ and the right hand side $\bar{q}_j^H$ express as in (16) and (15b) with $h$ replaced by $H$, respectively, and the starting vector is taken to be $0$. Thus, all computations with $L=0$ yield the exact mortar solution after one iteration and are not presented in the following.

#### 4.2. Algebraic error for static case

For the first set of tests, we set $\rho = 0$, $\nu = 0.3$, and piecewise constant Young moduli $E_\Xi = 100$, $E_\omega = 10^{\text{par}} E_\Xi$, par $\in \mathbb{Z}$. Further, we enforce homogeneous Dirichlet conditions on the upper boundary, homogeneous Neumann boundary conditions on the left and right
<figure><img src="image_6.png" /><figcaption>Fig. 5. Algebraic error reduction factor $e_{alg}^{(l+1)}/e_{alg}^{(l)}$ for the static case with $E_{\omega} = 10^{\text{par}E_{\Xi}}$ with respect to $l$; left: $\text{par} = 0, L \in \{1, \dots, 6\}$; right: $L = 2, \text{par} \in \{0, \dots, 5\}$.</figcaption></figure>

side and a surface load of $g_N = 10^6 \max(0.25 |x_1 - 1|, 0)$ at the bottom. The effective stress $\sigma_{\text{eff}} \|\sigma\|_d \text{tr}(\sigma)\text{Id}\|$ and the pressure $\text{tr}(\sigma)$ of the corresponding mortar solution for $\text{par} = 0$ are depicted in Fig. 3.

In order to investigate the decrease of the algebraic error for Algorithm 1 with respect to the number of iterations, the difference between the values $(\mathbf{u}^{H,(l)}, \mathbf{u}^{H,(l)})$ obtained by Algorithm 1 and the solution $(\mathbf{u}^H, \mathbf{u}^h)$ of the discrete mortar system (14) is measured with respect to the relative energy norm. For this, we compute

$$ (e_{\text{alg}}^{(l)})^2: \frac{(\mathbf{u}^{H,(l)} \quad \mathbf{u}^{H})^T A_{\Xi}^{H} (\mathbf{u}^{H,(l)} \quad \mathbf{u}^{H}) + (\mathbf{u}^{h,(l)} \quad \mathbf{u}^{h})^T A_{\omega}^{h} (\mathbf{u}^{h,(l)} \quad \mathbf{u}^{h})}{(\mathbf{u}^{H})^T A_{\Xi}^{H} \mathbf{u}^{H} + (\mathbf{u}^{h})^T A_{\omega}^{h} \mathbf{u}^{h}} $$

the numerator being equal to the value of (27). The results for this relative algebraic error and the corresponding estimator (28) are shown in Fig. 4 for different values of $L \in \{1, 2, 3, 4\}$. The left picture shows the results for equal material parameters $E_{\omega} = E_{\Xi}$, whereas the right picture displays the convergence for discontinuous parameters $E_{\omega} = 10^5 E_{\Xi}$. One can see that the decay rate with respect to the number of iterations is the same for the true and the estimated algebraic error and that the difference between them is very small, indicating that $\eta_{\text{alg}}$ is well suited to measure the algebraic energy error due to the iterative solution of the coupled system.

The dependence of the algebraic error reduction on the ratio $H/h = 2^L$ as well as on the jump in the material parameters $E_{\omega} = 10^{\text{par}E_{\Xi}}$ is further investigated in Fig. 5, where the error reduction factor $e_{\text{alg}}^{(l+1)}/e_{\text{alg}}^{(l)}$ is plotted with respect to $l$. The results show that the algorithm converges best for small values of $L$ and $\text{par}$, but the reduction factor is limited from above by around 0.35 for $\text{par} = 0$ independently of $L > 0$ and by around 0.55 for $L = 2$ and $\text{par} \to \infty$. Hence, in the static case, Algorithm 1 is stable with respect to

<table><caption>**Table 1**
Asymptotic error reduction rates $e_{\text{alg}}^{(l+1)}/e_{\text{alg}}^{(l)}$ for the dynamic case with different values of $E_{\omega}, 10^{\text{par}E_{\Xi}}, \mathcal{Q}_{\Xi}, 10^{\text{par}E_{\omega}}, \mathcal{Q}_{\Xi}$ with $\text{par}_{\mathcal{Q}}, \text{par}_{E} \in \{-4, -2, 0, 2, 4, 6\}$.</caption><thead><tr><th>L</th><th>par<sub>E</sub>\par<sub>Q</sub></th><th>-4</th><th>-2</th><th>0</th><th>2</th><th>4</th><th>6</th></tr></thead><tbody><tr><td rowspan="6">2</td><td>-4</td><td>0.0001</td><td>0.0037</td><td>0.2572</td><td>0.7357</td><td>0.7498</td><td>0.7500</td></tr><tr><td>-2</td><td>0.0034</td><td>0.0052</td><td>0.2449</td><td>0.7352</td><td>0.7498</td><td>0.7500</td></tr><tr><td>0</td><td>0.2311</td><td>0.2315</td><td>0.2792</td><td>0.6933</td><td>0.7498</td><td>0.7500</td></tr><tr><td>2</td><td>0.5441</td><td>0.5442</td><td>0.5473</td><td>0.5568</td><td>0.7065</td><td>0.7495</td></tr><tr><td>4</td><td>0.5462</td><td>0.5463</td><td>0.5516</td><td>0.5568</td><td>0.5624</td><td>0.7066</td></tr><tr><td>6</td><td>0.5462</td><td>0.5464</td><td>0.5516</td><td>0.5560</td><td>0.5558</td><td>0.5624</td></tr><tr><td rowspan="6">4</td><td>-4</td><td>0.0001</td><td>0.0043</td><td>0.3210</td><td>0.9196</td><td>0.9373</td><td>0.9375</td></tr><tr><td>-2</td><td>0.0039</td><td>0.0059</td><td>0.2887</td><td>0.9177</td><td>0.9373</td><td>0.9375</td></tr><tr><td>0</td><td>0.2617</td><td>0.2621</td><td>0.3159</td><td>0.8144</td><td>0.9373</td><td>0.9375</td></tr><tr><td>2</td><td>0.6063</td><td>0.6063</td><td>0.6095</td><td>0.6240</td><td>0.8298</td><td>0.9355</td></tr><tr><td>4</td><td>0.6090</td><td>0.6091</td><td>0.6139</td><td>0.6193</td><td>0.6302</td><td>0.8300</td></tr><tr><td>6</td><td>0.6091</td><td>0.6091</td><td>0.6138</td><td>0.6195</td><td>0.6194</td><td>0.6303</td></tr></tbody></table>

jumps in the coefficients, in accordance with the theoretical results of Theorem 2 for the case $\mathcal{Q}_{\omega} = 0$.

### 4.3. Algebraic error decrease for dynamic case

Next, we take the inertia terms into account and consider one time step of a dynamic setting with $\Delta t = 10^{-3}$ and piecewise constant material parameters $E_{\Xi}, 100, \mathcal{Q}_{\Xi}, 0.01, E_{\omega}, 10^{\text{par}E_{\Xi}}, \mathcal{Q}_{\omega}, 10^{\text{par}E_{\omega}}, \mathcal{Q}_{\Xi}$. The parameters are chosen such that the minimal value of the constant $C_{\text{CFL}}$ in (25) is given by

$$ C_{\text{CFL}} \frac{100}{64} \cdot 10^{\text{par}_{\omega} \text{ par}_{E}}. \qquad (29) $$

We impose homogeneous Dirichlet boundary conditions on the top, the Neumann compression forces

$$ g_N \cdot 10^{f(\text{par}_E, \text{par}_Q)} \cdot \max(0.25 r), \quad \text{with } r^2 (0.5 x_1)^2 $$

on the bottom and homogeneous Neumann conditions elsewhere. The value of the function $f(\text{par}_E, \text{par}_Q) \in [0, 9]$ is chosen such that the deformation of the body is within a sensible range. Both the volume load and the initial velocity are set to zero.

In Table 1, the algebraic error reduction rates of Algorithm 1 are summarized for $L \in \{2, 4\}$, $\text{par}_E, \text{par}_Q \in \{4, 2, 0, 2, 4, 6\}$. The results illustrate the influence of the constant $C_{\text{CFL}}$ in (25) and the bound in Remark 8. As long as $\text{par}_Q \le \text{par}_E + 2$ and thus $C_{\text{CFL}} \le 200$, the error reduction factor is bounded from above independently of the ratio of the grid sizes or the material parameters $E_{\omega}$ and $\mathcal{Q}_{\omega}$, corresponding to Theorem 2. However, for $\text{par}_Q > \text{par}_E + 4$, we get a reduction factor of $(1 h/H)$, as predicted by Remark 8.

## 5. Extension to nonlinear problems

Motivated by the challenge of simulating the large deformation dynamics of structures with contact, we have presented and analyzed a flexible iterative scheme that allows for a locally improved spatial resolution. In this section, we extend the considerations of Section 3 to the nonlinear case by incorporating nonlinear effects like material nonlinearity and contact with friction.

In Section 5.1, a general nonlinear setting is sketched, as well as an iterative solution scheme using the overlapping domain decomposition approach. Section 5.2 treats the case of dynamic frictional contact as an example of a nonlinear effect.

### 5.1. Iterative solution scheme

In comparison with the mortar system (9), the matrices $K_{\Theta}^m, \Theta \in \{\Xi, \omega\}, m \in \{h, H\}$, are now replaced by nonlinear opera tors $K_{\Theta}^m$. As we want to state a general algorithm which can be
applied to different kinds of nonlinearities, we do not give a precise definition of these operators. Some characteristic examples are presented in the next subsection, but we refer to, e.g., [12,38,48,49] for further information on geometrical or material nonlinearities.

The general nonlinear version of the mortar coupled discrete problem (12) reads

$$ \begin{pmatrix} K_{\Xi}^{H}(\mathbf{u}_{j}^{H}) \\ K_{\omega}^{h}(\mathbf{u}_{j}^{h}) \\ D^{hh}\mathbf{u}_{j}^{h} \quad D^{hh}\mathbf{u}_{j}^{H} \end{pmatrix} + \begin{pmatrix} (D^{hh})^{T} \\ (D^{hh})^{T} \\ \mathbf{0} \end{pmatrix} \zeta_{Fj}^{h} \mathbf{0}. \quad (30) $$

Above, we have included the volume and surface forces as well as the terms from the last time step in the definition of the maps $K_{\Theta}^{m}$, i.e., $K_{\Theta}^{m}(\mathbf{u}_{j}^{m})$, $K_{\Theta}^{m}(\mathbf{u}_{j}^{m}, \mathbf{v}_{j}^{m}, \mathbf{u}_{j}^{m}, \mathbf{v}_{j}^{m})$, but we omit the dependence on the latter two arguments for ease of notation. Further, the exact form of $K_{\Theta}^{m}$ depends on the kind of nonlinearity considered. In Sub section 5.2, we state the conditions of frictional contact as an example for a nonlinear problem with inequality constraints and sketch how these conditions can equivalently be reformulated as a set of semismooth equations. Hence, we assume in the following that (30) can be solved by a generalized form of the Newton method for semismooth systems. We refer to [21] and the references therein for an overview of Newton methods and to [15,50] for the definition of semismoothness.

### Algorithm 2. Inexact nonlinear two way coupling scheme

Start from some initial guess $\hat{\mathbf{z}}^{(0, l_{\max})}$.

Newton loop: Compute sequentially for $k = 0, 1, \dots$

(1) Initialize the solution $\hat{\mathbf{z}}^{(k+1, 0)} = \hat{\mathbf{z}}^{(k, l_{\max})}$.

Compute the tangential stiffness matrices and the residuals on $\Xi$ and $\omega$ relating to $S_{\Xi}^{H,(k)}$, $S_{\omega}^{H,(k)}$ and $\mathbf{q}_{\Xi}^{H,(k)}$, $\bar{\mathbf{q}}_{\omega}^{H,(k)}$.

Define coarse grid approximations $S_{\omega}^{H,(k)}$, $\bar{\mathbf{q}}_{\omega}^{H,(k)}$ of $S_{\omega}^{h,(k)}$, $\bar{\mathbf{q}}_{\omega}^{h,(k)}$.

(2) Gauß Seidel loop: Compute sequentially for $l = 0, 1, \dots$.

(i) Solve problem on coarse space $\mathbf{V}^{H}$ with interface load on $\Gamma$

$$ \begin{pmatrix} S_{\Xi}^{H,(k)} + S_{\omega}^{H,(k)} \\ S_{\omega}^{H,(k)} \end{pmatrix} \begin{pmatrix} \mathbf{0} \\ (D_{F}^{hh})^{T} \end{pmatrix} \begin{pmatrix} \delta \mathbf{u}_{Fj}^{H,(k+1,l)} \\ \delta \boldsymbol{\mu}_{Fj}^{H,(k+1,l)} \end{pmatrix} = \begin{pmatrix} \mathbf{r}_{\Xi}^{H,(k+1,l)} + \mathbf{r}_{\omega}^{H,(k+1,l)} \\ \mathbf{r}_{\omega}^{H,(k+1,l)} \end{pmatrix}, \quad (31) $$

(ii) Solve problem on fine space $\mathbf{V}^{h}$ with weakly imposed trace on $\Gamma$:

$$ \begin{pmatrix} S_{\omega}^{h,(k)} \\ D_{T}^{hh} \end{pmatrix} \begin{pmatrix} (D_{F}^{hh})^{T} \\ \mathbf{0} \end{pmatrix} \begin{pmatrix} \delta \mathbf{u}_{Fj}^{h,(k+1,l)} \\ \delta \boldsymbol{\zeta}_{Fj}^{h,(k+1,l)} \end{pmatrix} = \begin{pmatrix} \mathbf{r}_{\omega}^{h,(k+1,l)} \\ \mathbf{v}_{F}^{h,(k+1,l)} \end{pmatrix} + D_{F}^{hh} \delta \mathbf{u}_{Fj}^{H,(k+1,l)}, \quad (32) $$

(iii) Update the solution vector:

$$ \hat{\mathbf{z}}^{(k+1, l+1)} : \hat{\mathbf{z}}^{(k+1, l)} + \delta \hat{\mathbf{z}}^{(k+1, l)}, $$

The residuals of (31), (32) are given by

$$ \mathbf{r}_{\Xi}^{H,(k+1,l)} \quad \mathbf{q}_{\Xi}^{H,(k)} \quad S_{\Xi}^{H,(k)} \mathbf{u}_{Fj}^{H,(k+1,l)} + (D_{F}^{hh})^{T} \boldsymbol{\zeta}_{Fj}^{h,(k+1,l)}, \quad (33a) $$

$$ \mathbf{r}_{\omega}^{H,(k+1,l)} \quad \mathbf{q}_{\omega}^{H,(k)} \quad S_{\omega}^{H,(k)} \mathbf{u}_{Fj}^{H,(k+1,l)} \quad (D_{F}^{hh})^{T} \boldsymbol{\mu}_{Fj}^{H,(k+1,l)}, \quad (33b) $$

$$ \mathbf{r}_{\omega}^{h,(k+1,l)} \quad \mathbf{q}_{\omega}^{h,(k)} \quad S_{\omega}^{h,(k)} \mathbf{u}_{Fj}^{h,(k+1,l)} \quad (D_{F}^{hh})^{T} \boldsymbol{\zeta}_{Fj}^{h,(k+1,l)}, \quad (33c) $$

$$ \mathbf{v}_{F}^{h,(k+1,l)} \quad D_{F}^{hh} \mathbf{u}_{Fj}^{H,(k+1,l)} \quad D_{F}^{hh} \mathbf{u}_{Fj}^{h,(k+1,l)}. \quad (33d) $$

(iv) If $l+1 = l_{\max}$ or (35) is satisfied, set $\hat{\mathbf{z}}^{(k+1, l_{\max})} = \hat{\mathbf{z}}^{(k+1, l)}$ and stop.

(3) Check convergence of Newton iteration.

In the following, we transfer the idea of Algorithm 1 to (30). Be cause of the nonlinearity, we have to combine two iterative pro cesses, namely the semismooth Newton loop for the nonlinear effects and the subdomain iteration presented in Section 3. The natural way of doing so is to use Algorithm 1 to solve the tangential systems of the semismooth Newton iteration applied to (30). The corresponding scheme is summarized in Algorithm 2, using the notation $\hat{\mathbf{z}}$ introduced in (17), with $S_{\Theta}^{m,(k)}$ denoting the Schur com plement of the tangential matrix

$$ K_{\Theta}^{m,(k)} : (\partial_{\mathbf{u}_j^m} K_{\Theta}^{m})(\mathbf{u}_j^{m,(k)}), \quad \Theta \in \{\Xi, \omega\}, \quad m \in \{h, H\}, \quad (34) $$

and the right hand side vectors

$$ \mathbf{q}^{H,(k)} : K_{\Xi}^{H,(k)} \mathbf{u}_{j}^{H,(k)}, \quad K_{\Xi}^{H}(\mathbf{u}_{j}^{H,(k)}), \quad \mathbf{q}^{h,(k)} : K_{\omega}^{h,(k)} \mathbf{u}_{j}^{h,(k)}, \quad K_{\omega}^{h}(\mathbf{u}_{j}^{h,(k)}). $$

The convergence rate of the inner Gauß Seidel iteration can be analyzed with the results from Section 3. Furthermore, the super linear local convergence of the outer Newton iteration can be pre served if the norm of the residual of the inner iteration is bounded in terms of the norm of the Newton residual

$$ \mathbf{R}^{(k)} : \begin{pmatrix} K_{\Xi}^{H}(\mathbf{u}_{j}^{H,(k)}) & (D^{hh})^{T} \boldsymbol{\zeta}_{Fj}^{h,(k)} \\ K_{\omega}^{h}(\mathbf{u}_{j}^{h,(k)}) + (D^{hh})^{T} \boldsymbol{\zeta}_{Fj}^{h,(k)} & \end{pmatrix} $$

(see, e.g., [20, 22]). For $l \ge 1$, the residual of the inner loop is given by the vector (33a), such that the Gauß Seidel iteration has to be solved until the condition

$$ \| \mathbf{r}_{\Xi}^{H,(k+1,l)} \| \mathcal{O}(\| \mathbf{R}^{(k)} \|) \quad (35) $$

is satisfied.

**Remark 10.** In principle, the two iterative processes can also be nested in a different way, leading to an outer fixed point loop and an inner Newton iteration. However, as a Newton step is, in general, more expensive than a linear fixed point step due to the reassembly of the stiffness matrix, this version is likely to be less efficient than the latter one, see also [39] for numerical results.

## 5.2. Frictional contact

In this subsection, we exemplify the formation of the nonlinear semismooth operators $K_{\Theta}^{m}$ used in (30) for the cases of a nonlinear displacement stress relationship and frictional contact. For the for mer effect, the derivation of the function $K_{\Theta}^{m}$ proceeds along the lines of Section 2.1, but with the Cauchy stress tensor $\boldsymbol{\sigma}$ in (1) being replaced by the first Piola Kirchhoff stress tensor $\mathbf{P}(\mathbf{u})$. The latter is given by $\mathbf{P}(\mathbf{u}) = (\mathrm{Id} + \nabla \mathbf{u})\mathbf{S}(\mathbf{u})$, with $\mathbf{S}(\mathbf{u})$ denoting the second Piola Kirchhoff stress tensor whose definition depends on the material. In some of the numerical results, we consider the compressible equivalent of the nonlinear Mooney Rivlin material law given in [39,51], where $\mathbf{S}$ is defined in terms of the Lamé parameters $\mu, \lambda$, the additional parameter $c_m$ and the right Cauchy Green tensor $\mathbf{C} = (\mathrm{Id} + \nabla \mathbf{u})^T (\mathrm{Id} + \nabla \mathbf{u})$ via

$$ \mathbf{S} = \frac{\lambda}{2} (\det(\mathbf{C})) \mathbf{1}\mathbf{C}^1 + \mu \left( (\mathbf{1} \quad c_m)(\mathrm{Id} \quad \mathbf{C}^1) + c_m ((\mathrm{tr}\mathbf{C})\mathrm{Id} \quad \mathbf{C} \quad (\mathrm{d} \quad \mathbf{1})\mathbf{C}^1) \right). \quad (36) $$

Due to the nonlinear relationship between $\mathbf{u}$ and $\mathbf{P}$, the stiffness terms in (7) are now nonlinear (i.e., $A_{\Xi}^{H}(\mathbf{u}^{H})$ and $A_{\omega}^{H}(\mathbf{u}^{H})$). Further, special care must be taken by applying the time discretization. We refer to, e.g., [12,23,25] for more details.

Next, we consider the constraints for frictional contact in more detail. Here, we designate a part $\gamma \subset \Gamma_N$ of the Neumann boundary where the body $\Omega$ can come into contact with a fixed smooth
obstacle $\Gamma_{\text{obs}}$; for the contact of several elastic bodies we refer to, e.g., [17,52]. On $\gamma$, we enforce the Neumann conditions $\mathbf{Pn} = \lambda$, where the contact stress $\lambda$ is a priori unknown and has to be deter mined by the contact conditions.

In general, a detailed resolution of the contact area is desired, and thus it is reasonable to assume that the patch $\omega$ is large en ough to cover the whole contact boundary, i.e., $\bar{\gamma} \subset \partial\omega$. Further, we assume $\bar{T} \cap \bar{\gamma} \neq \emptyset$, $\bar{T}_p \cap \bar{\gamma} \neq \emptyset$ for simplicity and denote the set of degrees of freedom of $\mathbf{V}_{\omega}^m$ on $\gamma$ by $\mathcal{N}_{\gamma}^m \subset \mathcal{N}^m$, $m \in \{h, H\}$.

The discrete conditions for normal contact with Coulomb or Tresca friction [11, 13] can be described in terms of the fine scale displacement $\mathbf{u}_{\gamma}^h$ ($\mathbf{u}_p^h$)$_{p \in \mathcal{N}_{\gamma}^h}$ and the Lagrange multipliers $\lambda_{pjn}^h$ ($\lambda_p^h$)$_{p \in \mathcal{N}_{\gamma}^h}$. For the discretization of $\lambda_{pjn}^h$, we employ dual basis functions $\{\tilde{\psi}_p^h\}_{p \in \mathcal{N}_{\gamma}^h}$ which are characterized by a biorthogonality condition similar to (4). This allows for a separation of the contact conditions into individual constraints for each contact node $p \in \mathcal{N}_{\gamma}^h$ (see [16,17] for more details). Since we consider unilateral contact with a fixed obstacle, these conditions can be formulated in terms of the local displacement and the corresponding multiplier. For this, we split the displacement vector $\mathbf{u}_p^h$ into its normal and tan gential part according to $\mathbf{u}_{pn}^h = \mathbf{u}_p^h \mathbf{n}_p$, $\mathbf{u}_{pt}^h = \mathbf{u}_p^h \mathbf{u}_{pn}^h \mathbf{n}_p$, where $\mathbf{n}_p$ denotes a suitable approximation of the unit outer normal of $\omega$ at $p \in \mathcal{N}_{\gamma}^h$. Further, we introduce the initial normal gap $g_{pn}$ be tween the node $p$ and the fixed obstacle $\Gamma_{\text{obs}}$ in the direction of $\mathbf{n}_p$ (see, e.g., [12,53] for a detailed definition), as well as the normal and tangential parts $\lambda_{pjn}^h = \lambda_p^h \mathbf{n}_p$, $\lambda_{pt}^h = \lambda_p^h \lambda_{pn}^h \mathbf{n}_p$ of the Lagrange multiplier. With this, the constraints in the normal direction at the time step $t_j$ read

$$ \mathbf{u}_{pjn}^h \quad g_{pn} \le 0, \quad \lambda_{pjn}^h \ge 0, \quad \lambda_{pjn}^h (\mathbf{u}_{pjn}^h \quad g_{pn}) \quad 0. \quad (37) $$

The first inequality of (37) comprises the non penetration con dition at $p$, whereas the second one ensures that the force in nor mal direction is compressive. The last complementarity equation states that a nonzero contact stress can only develop if the gap is closed.

The constraints for frictional contact with either Coulomb or Tresca friction are given by

$$ \begin{aligned} \|\lambda_{pjt}^h\| &\le g_{pt} + \tilde{\gamma} |\lambda_{pjn}^h|, \\ \|\lambda_{pjt}^h\| &< g_{pt} + \tilde{\gamma} |\lambda_{pjn}^h| \Rightarrow \partial_{\Delta t} \mathbf{u}_{pjt}^h \quad 0, \\ \|\lambda_{pjt}^h\| &\ge g_{pt} + \tilde{\gamma} |\lambda_{pjn}^h| \Rightarrow \beta\beta \ge 0 : \lambda_{pjt}^h \quad \beta\partial_{\Delta t} \mathbf{u}_{pjt}^h. \end{aligned} \quad (38) $$

If the friction coefficient $\tilde{\gamma} \ge 0$ is set to zero, the above inequalities describe the case of Tresca friction with the given friction bound $g_{pt} \ge 0$, whereas for $g_{pt} = 0$, the Coulomb friction law is obtained.

There are several possibilities to rewrite the above inequality constraints as a set of nonsmooth equations [19,50,54]; the form we propose employs the augmented trial

$$ \lambda_{pjn}^{\text{htr}} : \quad \lambda_{pjn}^{h} + c_n (\mathbf{u}_{pjn}^{h} \quad g_{pn}), \quad (39a) $$

$$ \lambda_{pjt}^{\text{htr}} : \quad \lambda_{pjt}^{h} + c_t (\mathbf{u}_{pjt}^{h} \quad \mathbf{u}_{p(j \ 1)t}^{h}) \quad (39b) $$

with some fixed constants $c_n, c_t > 0$, as well as the combined friction bound

$$ \tilde{g}_{pjt} : \quad g_{pt} + \tilde{\gamma} \max(0, \lambda_{pjn}^{\text{htr}}). $$

Then, we can rewrite the contact conditions (37) and (38) as

$$ \mathbf{C}_{\gamma}^{h}(\mathbf{u}_{jj}^{h}, \lambda_{\gamma j}^{h}) : \quad (\mathbf{C}_{p}^{h}(\mathbf{u}_{pj}^{h}, \lambda_{pj}^{h}))_{p \in \mathcal{N}_{\gamma}^{h}} \quad 0 $$

with the nonlinear complementarity (NCP) function

$$ \mathbf{C}_{p}^{h}(\mathbf{u}_{pj}^{h}, \lambda_{pj}^{h}) : \quad \left( \begin{array}{cc} \lambda_{pjn}^{h} & \max(0, \lambda_{pjn}^{\text{htr}}) \\ \max(\tilde{g}_{pjt}, \|\lambda_{pjt}^{\text{htr}}\|) & \lambda_{pjt}^{h} \\ \hline \end{array} \right) \quad (40) $$

The equivalence of (37), (38) and (40) follows from a case by case analysis according to the values of the max functions in (40). Further, this analysis directly leads to a partition of the set $\mathcal{N}_{\gamma}^h$ into subsets of active and inactive nodes denoted by $\mathcal{A}_n^h$, $\mathcal{A}_t^h$ and $\mathcal{I}_n^h$, $\mathcal{I}_t^h$, respectively, and defined via

$$ \mathcal{A}_n^h : \quad \{p \in \mathcal{N}_{\gamma}^h : \lambda_{pjn}^{\text{htr}} > 0\}, \quad \mathcal{I}_n^h : \quad \mathcal{N}_{\gamma}^h \setminus \mathcal{A}_n^h, \quad (41a) $$

$$ \mathcal{A}_t^h : \quad \{p \in \mathcal{N}_{\gamma}^h : \|\lambda_{pjt}^{\text{htr}}\| > \tilde{g}_{pjt} \text{ or } \tilde{g}_{pjt} \quad 0\}, \quad \mathcal{I}_t^h : \quad \mathcal{N}_{\gamma}^h \setminus \mathcal{A}_t^h. \quad (41b) $$

According to this partition, the NCP function (40) simplifies to

$$ C_p^h(\mathbf{u}_{pj}^h, \lambda_{pj}^h) = \left\{ \begin{array}{ll} \lambda_{pjn}^h, & p \in \mathcal{I}_n^h, \\ c_n(g_{pn}, \mathbf{u}_{pjn}^h), & p \in \mathcal{A}_n^h, \\ \tilde{g}_{pjt} c_t(\mathbf{u}_{pjt}^h, \mathbf{u}_{p(j \ 1)t}^h), & p \in \mathcal{I}_t^h, \\ \|\lambda_{pjt}^{\text{htr}}\| \lambda_{pjt}^h, & p \in \mathcal{A}_t^h, \end{array} \right. \quad (42) $$

Thus, the set $\mathcal{A}_n^h$ contains all nodes which are in contact with the obstacle, whereas $\mathcal{A}_n^h$ and $\mathcal{I}_t^h$ represent the slippery and sticky nodes, respectively. We remark that (42)$_{2,3}$ can be seen as partial Dirichlet conditions for the normal or tangential displacement, whereas (42)$_1$ is a partial Neumann and (42)$_4$ a Robin type condition.

**Remark 11.** If the normal $\mathbf{n}_p$ does not coincide with one of the co ordinate basis vectors, the partial boundary conditions are enforced by a local basis transformation with $\mathbf{n}_p$ being one of the new orthonormal basis vectors. We refer to [16] for details.

For the implementation of these boundary conditions within the $k$th iteration of the semismooth Newton loop, we compute the generalized derivative of the NCP function (42) which is given in Appendix B. This leads to intermediate active sets $\mathcal{A}_n^{h,(k)}$, $\mathcal{A}_t^{h,(k)}$ which correspond to the final active sets (41) after the Newton iteration has converged. The linearized version of (42) is then implemented by local static condensation of the dual variables. De tails can be found in [54,16].

**Remark 12.** As shown in [12], the normal contact conditions (37) in combination with the trapezoidal rule yield an unstable algorithm if inertia terms are present. In this case, it is advisable to change the time discretization of (37) like, e.g., described in [27,55]. But even then, the additional problem occurs that the computed contact stresses show spurious oscillations in time [30]. In order to avoid this, we employ a local modification of the mass matrix such that its entries associated with the potential contact nodes vanish. We refer to [28, 30] for more details.

In the rest of this subsection, we present a possible definition of the coarse grid quantities $S_{\omega}^{H,(k)}$, $\tilde{q}_{\omega}^{H,(k)}$ needed in Algorithm 2. Numerical tests like the one in Subsection 6.1 show that the con vergence rate of the inner iteration in Algorithm 2 degrades if the Dirichlet contact conditions, i.e., the active nodes $\mathcal{A}_n^{H,(k)}$ in nor mal direction and the sticky nodes $\mathcal{I}_t^{H,(k)}$ in tangential direction, are not respected within the coarse approximation. Hence, we propose to construct $S_{\omega}^{H,(k)}$, $\tilde{q}_{\omega}^{H,(k)}$ as their fine grid counterparts, but with re spect to sets $\mathcal{A}_n^{H,(k)}$, $\mathcal{I}_t^{H,(k)} \subset \mathcal{N}_{\gamma}^H$ approximating the fine contact zone. To construct the coarse active set $\mathcal{A}_n^{H,(k)}$, we define the finite ele ment function $\chi_{\gamma,n}^{H,(k)} \in \mathbf{W}_{\gamma}^h$ given by

$$ \chi_{\gamma n}^{H,(k)}(p) = \begin{cases} 1, & p \in \mathcal{A}_n^{H,(k)}, \\ 0, & \text{otherwise.} \end{cases} \quad (43) $$

Introducing the mortar operator $P_{\gamma}^H : \mathbf{W}_{\gamma} \to \mathbf{W}_{\gamma}^H$ defined according to (5), we choose an appropriate threshold value $\tau_n \in \mathbb{R}$ and set

$$ \mathcal{A}_n^{H,(k)}(\tau_n) : \quad \{p \in \mathcal{N}_{\gamma}^{H} : (P_{\gamma}^{H} \chi_{\gamma n}^{H,(k)})(p) > \tau_n\}. \quad (44) $$
The rest of the coarse nodes is set inactive, i.e., $\mathcal{I}_n^{H,(k)} : \mathcal{N}_\gamma^H \setminus \mathcal{A}_n^{H,(k)}$. We remark that smaller values of $\tau_n$ generally lead to a larger coarse contact set $\mathcal{A}_n^{H,(k)}$.

For the tangential part, we proceed similarly and define the set $\mathcal{I}_t^{H,(k)}$ of coarse sticky nodes by

$$ \mathcal{I}_t^{H,(k)}(\tau_t) : \left\{ p \in \mathcal{N}_\gamma^H : (P_\gamma^H \chi_{\gamma t}^{h,(k)})(p) > \tau_t \right\} \quad (45) $$

with $\tau_t \in \mathbb{R}$ and

$$ \chi_{\gamma t}^{h,(k)}(p) \begin{cases} \mathbf{1}, & p \in \mathcal{I}_t^{h,(k)}, \\ \mathbf{0}, & \text{otherwise.} \end{cases} \quad (46) $$

Due to the combination of boundary displacement and stress, the Robin conditions for the slippy nodes are not as easy to transfer to the coarse grid. Thus, we free the remaining nodes $\mathcal{A}_t^{H,(k)} : \mathcal{N}_\gamma^H \setminus \mathcal{I}_t^{H,(k)}$ in tangential direction by enforcing $\dot{\chi}_t^{H,(k+1)} = \mathbf{0}$.

**Remark 13.** The coarse representation of the Dirichlet constraints is similar to the construction of the coarse grid corrections in recent monotone multigrid methods; see, e.g., [56] and the references therein.

**Remark 14.** For the action of the operator $P_\gamma^H$ introduced above, we need to compute the integrals $\int_\gamma \psi_p^H \phi_q^H ds$ for all $p \in \mathcal{N}_\gamma^H$, $q \in \mathcal{N}_\gamma^H$. Depending on the local geometry, the evaluation of these integrals can be rather cumbersome (see, e.g., [44]); however, if the relative position of the coarse and the fine grid is not changed, they only have to be computed once.

## 6. Numerical tests for the nonlinear setting

In this section, we show some numerical results with two different three dimensional geometries.

### 6.1. Contact approximation with conforming geometry

First, we investigate the influence of the coarse grid approximations $\mathcal{A}_n^H(\tau_n)$, $\mathcal{I}_t^H(\tau_t)$ defined in (44) and (45) on the convergence rate. For this, we consider the simple three dimensional domain $\Omega = (0,1) \times (0,1) \times (0,2)$ with the patch $\omega = (0,1)^3$. We impose a fixed displacement of $(0,0, 0.03)^T$ on the top, unilateral contact with the obstacle

$$ \Gamma_{\text{obs}}(x_1, x_2) = (x_1, x_2, 0.25x_1 \cdot \sin(4\pi x_1) \cdot x_2 \cdot \sin(4\pi x_2)) $$

on the bottom (see the left picture of Fig. 6) and homogeneous Neumann conditions elsewhere. The material parameters for the linear elastic body are $E_\omega = E_\Xi = 100$, $\nu = 0.33$ and $\rho = 0$, and we use uni form hexahedral grids of mesh size $H = \frac{1}{4}$ and $h = \frac{1}{16}$.

We consider frictionless as well as frictional contact, the latter case with $g_t = 0$ and $\mathcal{F} \in \{0, 0.02, 0.2, 0.35, 0.5, 1, 5\}$. In Fig. 6, the fine active set $\mathcal{A}_n^h$ of the converged solution for $\mathcal{F} = 0$ is depicted, as well as the coarse grid approximations $\mathcal{A}_n^H(\tau_n)$ for $\tau_n \in \{0,0.4\}$. Table 2 summarizes the error reduction rates for Algorithm 2 with one inner Gauß Seidel step ($l_{\max} = 1$) and different values of $\mathcal{F}$, $\tau_n$ and $\tau_t$. The upper subtable contains the frictionless case with $\tau_n \in \{1, 0, 0.1, 0.2, 0.4, 2\}$, where the threshold values $\tau_n = 1$ and $\tau_n = 2$ correspond to the cases $\mathcal{I}_n^H = \emptyset$ and $\mathcal{A}_n^H = \emptyset$, respectively. One can see that the convergence rate is around 0.3 as long as $\tau_n \le 0.2$, but for larger values of $\tau_n$, the inner iterations of Algorithm 2 become less effective.

The remaining rows of Table 2 refer to the frictional case where we approximate the normal active set by (44) with $\tau_n = 0$ and compare the error reduction rates for different values of $\tau_t$. The notation div indicates that the method has diverged, whereas osc refers to an oscillation of the active sets. One can observe that the convergence rate is best for the threshold value $\tau_t = 0$ whose error reduction is rather independent of the value of the friction coefficient. In contrast, smaller coarse sets $\mathcal{I}_t^H$ obtained by larger values of $\tau_t$ yield a more unstable convergence behavior, and the extreme

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Frictionless contact example with conforming geometry; from left to right: isolines of obstacle, fine normal active set $\mathcal{A}_n^h$ (black nodes), and approximations $\mathcal{A}_n^H(\tau_n)$ for $\tau_n \in \{0,0.4\}$.</figcaption></figure>

<table><caption>**Table 2**
Asymptotic error reduction rates for Algorithm 2 and proportions of fine active nodes for $L=2$, $l_{\max} = 1$ and different values of $\mathcal{F}$; first row: $\mathcal{F} = 0$, $\tau_n \in \{-1, 0, 0.1, 0.2, 0.4, 2\}$; other rows: $\mathcal{F} > 0$, $\tau_n = 0$, $\tau_t \in \{-1, 0, 0.1, 0.2, 0.4, 2\}$.</caption><thead><tr><th>\tau_n</th><th>-1</th><th>0</th><th>0.1</th><th>0.2</th><th>0.4</th><th>2</th><th>\mathcal{A}_n^h/|\mathcal{A}_n^H|</th></tr></thead><tbody><tr><th>\mathcal{F} \tau_n</th><td>0</td><td>0.304</td><td>0.304</td><td>0.304</td><td>0.306</td><td>0.625</td><td>0.830</td><td>77/289 ≈ 0.27</td></tr><tr><th>\mathcal{F} \tau_n</th><td>-1</td><td>0</td><td>0.1</td><td>0.2</td><td>0.4</td><td>2</td><td>\mathcal{A}_n^h/|\mathcal{A}_n^H|</td><td>\mathcal{A}_n^h/|\mathcal{A}_n^H|</td></tr><tr><th>0.02</th><td>0.741</td><td>0.305</td><td>0.305</td><td>0.305</td><td>0.305</td><td>0.305</td><td>75/289 ≈ 0.26</td><td>0</td></tr><tr><th>0.2</th><td>0.371</td><td>0.364</td><td>Div</td><td>Div</td><td>Div</td><td>Div</td><td>51/289 ≈ 0.18</td><td>51/42 ≈ 0.12</td></tr><tr><th>0.35</th><td>0.305</td><td>0.302</td><td>Div</td><td>0.699</td><td>Div</td><td>Div</td><td>51/289 ≈ 0.18</td><td>51/51 ≈ 0.41</td></tr><tr><th>0.5</th><td>0.303</td><td>0.299</td><td>0.303</td><td>0.730</td><td>Div</td><td>Div</td><td>51/289 ≈ 0.18</td><td>51/51 ≈ 0.63</td></tr><tr><th>1</th><td>0.301</td><td>0.299</td><td>0.299</td><td>osc</td><td>osc</td><td>Div</td><td>51/289 ≈ 0.18</td><td>41/51 ≈ 0.80</td></tr><tr><th>5</th><td>0.302</td><td>0.300</td><td>0.301</td><td>0.343</td><td>0.943</td><td>Div</td><td>51/289 ≈ 0.17</td><td>48/48 ≈ 0.94</td></tr></tbody></table>
case $I_t^H \emptyset$ obtained for $\tau_t = 2$ hardly converges at all. This suggests that for a good convergence rate, one should choose either $\tau_t$ small enough (less than 0.1 in our example) or try to improve the coarse grid approximation for slippy nodes, possibly by using a coarse Ro bin condition or an approach similar to [56].

**Remark 15.** The results indicate that in the static case, the convergence rate of Algorithm 2 is quite sensitive to the threshold values $\tau_n$, $\tau_t$. Therefore, we suggest to choose these parameters small enough such that the supports of the fine grid indicator functions (43) and (46) are contained in their coarse counterparts.

However, if a mass term is included in the computation, the convergence rate is less sensitive to the approximation of the contact boundary conditions. For the frictionless contact problem with $\rho = 10^{-2}$, $\Delta t = 10^{-3}$, we have obtained a convergence rate of 0.29 independent of the threshold value $\tau_n$.

## 6.2. Tire application

Next, we apply our algorithm to a complex 3D geometry consisting of the lower half of a car tire centering on the $x_2$ axis and with an approximate radius of 318. In Fig. 7, the domain $\Omega$ as well as the multiply connected fine patch $\omega$ are sketched, with a difference in the mesh sizes of about $4 \le H/h \le 8$. In contrast to all previous numerical results, the geometry is nonconforming because the fine triangulation $T^h$ features some details at the potential contact boundary $\gamma$ which are not resolved by the coarse grid $T^H$ (see Fig. 8).

First, we restrict ourselves to the case of linear static elasticity with frictionless contact. The contact plane is located at $x_3 = 320$, and the coarse approximation for the normal active set computed with $c_n = 10^6$ is done using (44). On the cutting faces

of the tire (for $x_3 = 0$), a fixed displacement of (0,0, 1) is prescribed, whereas all other boundaries are free. Further, a volume force of (0,0, 1000) is applied to the tire, and the material parameters read $E_\Sigma = E_\omega = 2.5 \ 10^7$ and $\nu = 0.33$.

Local reinforcement, incompressibility and anisotropy have been voluntarily ignored, the aim of the present case mainly being to analyse the performance of the present method on a nonlinear contact problem with complex 3D geometry. The active set $A_n^h$ of the corresponding mortar solution is sketched on the left of Fig. 8, whereas its coarse grid approximation $A_n^H$ (0.15) can be seen on the right side.

Fig. 9 displays the error decay and the mean error reduction factor of the damped version of Algorithm 2 with the damping parameter $\alpha = 0.6$, $l_{max} = 4$ inner iterations, and different threshold values $\tau_n$. The results show that on the one hand, a too small threshold value can lead to oscillations in the coarse active set such that the algorithm does not converge. On the other hand, if $\tau_n$ is too large, the coarse contact set becomes smaller which degrades or even disables the convergence, as already observed in Section 6.1.

From now on, we fix $\tau_n = 0.15$ and investigate the performance of Algorithm 2 for different damping parameters $\alpha$. On the left side of Fig. 10, the error decay for $l_{max} = 4$ and $\alpha \in \{0.5,0.6,0.7,0.8\}$ is shown. For larger values of $\alpha$, generally more Newton steps are necessary in order to detect the correct active set. But as soon as the correct active set has been found, the mean error reduction factor is better for larger values of $\alpha$, varying between 0.49 for $\alpha = 0.7$ and 0.66 for $\alpha = 0.5$. This suggests an adaptive choice of $\alpha$ with respect to a change in the active sets.

In order to decrease the number of total iterations, we replace the condition $l \le l_{max}$ by the following stopping criterion in the spirit of (35):

<figure><img src="image_8.png" /><figcaption>Fig. 7. Geometry of the 3D tire example.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Fine contact set $A_n^h$ of exact solution with 555 nodes and approximation $A_n^H$ (0.15) with 36 nodes.</figcaption></figure>
$$ (\eta_{\text{abs}}^{(k+1,l)})^2 \le E_{\omega} \|\mathbf{R}^{(k)}\|^2. \quad (47) $$

The corresponding error decay is depicted on the right of Fig. 10. Comparing these results with those in the left picture, one can observe that the number of total iterations until convergence has been reduced by around 10, and that also the computation with $\alpha = 0.8$ is able to detect the correct active set within 15 Newton steps.

Finally, we show a dynamic example with both material and contact nonlinearities. For this, we employ the nonlinear material law previously stated in (36) with the parameters $E_{\Xi} = E_{\omega} = 2.5 \ 10^7$, $\nu = 0.33$, $\rho = 10^5$ and $c_m = 0.5$. We apply no volume forces and enforce frictional contact with the coefficient $\mathcal{F} = 0.5$. The linearized contact conditions are assembled in an updated Lagrangian manner, and we employ a redistribution of the mass near the contact boundary $\gamma$ as described in [29]. We compute 10 time steps with a step size of $\Delta t = 10^{-5}$ and prescribe

a time dependent displacement of $\frac{t}{10\Delta t}(0, 0, 0.75 + \text{sign}(x_1))$ disp, $\text{disp} \in \{0, 2\}$, on the top, whereas all other boundaries are free. The effective stress $\sigma_{\text{eff}}$ of the corresponding solutions at time $t_{10}$ is depicted in Fig. 11. The problem is solved using Algorithm 2 with $\alpha = 0.6$, $\tau_n = \tau_t = 0.15$ and $l_{\text{max}} = 1$. On the left of Fig. 12, the error reduction factor for the first 10 iterations of the time steps $t_2$ to $t_{10}$ is plotted, together with the geometric mean for each time step. The right figure shows the evolution of the fine active sets $|A_n^h|$ and $|T_t^h|$ with respect to time. One can observe that the second test setting with disp = 2 yields no sticky nodes and thus gives a slightly better convergence rate than the problem with disp = 0 where sticky and slippy nodes appear. However, the mean error reduction factor is between 0.4 and 0.5 for both settings, illustrating that Algorithm 2 is able to handle different types of nonlinearities within a single inexact Newton loop even on complex geometries.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Performance of Algorithm 2 with $\alpha = 0.6$, $l_{\text{max}} = 4$ w.r.t. number of total iterations for $\tau_n \in (0.05, 0.1, 0.15, 0.2, 0.3, 0.4)$; left: true and estimated rel. alg. error; right: Size of $A_n^H(\tau_n)$ and mean error red. factor after correct active set has been found.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. True $e_{\text{alg}}$ and estimated $\eta_{\text{alg}}$ relative alg. error of Algorithm 2 with $\tau_n = 0.15$ for $\alpha \in (0.5, 0.6, 0.7, 0.8)$; left: $l_{\text{max}} = 4$; right: inner iteration stopped according to (47).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Effective stress of dynamic tire problem at time $t_{10}$; left: disp = 0 and right: disp = 2.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Performance of Algorithm 2 with $\alpha = 0.6$ and $\tau_n = \tau_t = 0.15$ for the dynamic tire problem with frictional contact and nonlinear material. Left: error reduction factor for the first 10 iterations per time step and right: evolution of fine active sets.</figcaption></figure>

**Remark 16.** A crucial feature of the present approach is the fact that the global coarse grid solution can serve as a predictor of the effective contact zone. Using this prediction, a subrefinement of the coarse mesh can be generated automatically to provide a more accurate description of the deformation in the contact region. The contact detection and the corresponding refinement can, e.g., be reevaluated in each iteration after the solution of the coarse model and the corresponding update of the coarse active sets. Combined with the local time step refinement presented in the next section, this idea makes the efficiency of the present approach superior to most standard domain decomposition methods.

## 7. Local time subcycling

The domain decomposition approach introduced in the previous sections has been motivated by the fact that one is interested in a detailed resolution of the solution within the fine patch. Hence, in addition to the different spatial mesh sizes, it is desirable to use different time scales for the coarse and the fine problems, which is the topic of the current section. In Section 7.1, we present a mortar coupled solution with different time scales and analyse its asymptotic convergence order for $\Delta t \to 0$. Afterwards, we sketch in Section 7.2 how the discrete system can efficiently be solved using an iterative scheme based on the overlapping domain decomposition. Numerical tests are presented in Section 8.

For ease of presentation, we restrict ourselves to the linear case. The extension to nonlinear problems can be done as in Section 5. Further, we assume that the fine scale mass matrix $M_\omega^h$ is constricted in such a way that its entries associated with the interface $\Gamma$ vanish. Such a modification has already been successfully applied to dynamic contact problems (cf. Remark 12 and [30,28]), and the corresponding a priori error in space has been analyzed in [29]. The time subcycling algorithm with the standard mass matrix is presented in [57].

### 7.1. Interface coupling conditions

From now on, we denote the time step associated with the coarse space $V^H$ by $\Delta T$ and the time step on the patch $\omega$ by $\Delta t$. We assume that the time step sizes have an integer ratio, i.e., $\frac{\Delta T}{\Delta t}$ $J$ holds for some $J \in \mathbb{N}$. For ease of notation, the time indices refer to the fine time steps, such that the first macro time step lasts from $t_0$ to $t_j$. To derive continuity conditions on the interface for the case $J > 1$, we define intermediate values of the coarse solution at the discrete time steps $t_j = j\Delta t, j = 1, \dots, J$, by means of linear interpolation in time. The reasons for this choice will become clear in the following, where we show that the resulting coupled system

of equations provides conservation of the total energy as well as optimal a priori error estimates in time on the interior of the patch.

First, we look at the solution of one fine time step at the patch $\omega$, where we want to ensure the weak continuity of the displacement on $\Gamma$ with the linearly interpolated values from $\Xi$ (cf. (14)$_{2,3}$). If the coarse solution $u_{\Gamma J}^H$ at time $t_j$ is given, this yields

$$ S_{\omega}^{h} (D_{\Gamma}^{hh})^{T} \begin{pmatrix} \mathbf{u}_{\Gamma J}^{h} \\ \zeta_{\Gamma J}^{h} \end{pmatrix} = \begin{pmatrix} \mathbf{q}_{\omega j \ 1}^{h} \\ D_{\Gamma}^{hh} \begin{pmatrix} j & \mathbf{u}_{\Gamma J}^{H} \\ \mathbf{1} & j \end{pmatrix} + \begin{pmatrix} \mathbf{1} & \mathbf{u}_{\Gamma 0}^{H} \end{pmatrix} \end{pmatrix}. \quad (48) $$

The feedback to the coarse grid Eq. (14)$_1$ is carried out by the mean value of the Lagrange multipliers, leading to

$$ S_{\Xi}^{H} \mathbf{u}_{\Gamma J}^{H} = \mathbf{q}_{\Xi 0}^{H} + (D_{\Gamma}^{hh})^{T} \frac{1}{J} \sum_{j=1}^{J} \zeta_{\Gamma J}^{h}. \quad (49) $$

The Schur matrices $S_{\Xi}^H, S_\omega^h$ used above have been defined in (13), and the residuals $\hat{\mathbf{p}}_{\omega 0}^H, \hat{\mathbf{p}}_{\omega j}^H$ are given in (15).

The reason for implementing the backcoupling as stated in (49) is the following lemma which follows directly from the definitions:

**Lemma 3.** For time independent loads $f$, the coupled system (48), (49) conserves the total energy (10) after one coarse time step, i.e., $E_J E_0 = 0$.

Besides the fact that the solution of (48), (49) is energy conserving, it permits the a priori error estimates stated in the next theorem.

**Theorem 4.** Let the initial conditions $(\mathbf{u}_0^H, \mathbf{v}_0^H, \mathbf{u}_0^h, \mathbf{v}_0^h)$ be given. Let $(\mathbf{u}_J^H, \mathbf{v}_J^H, \mathbf{u}_1^h, \mathbf{v}_1^h, \zeta_1^h, \dots, \mathbf{u}_J^h, \mathbf{v}_J^h, \zeta_J^h)$

be the complete solution of (48), (49) with a mass matrix $M_\omega^h$ such that all entries associated with the nodes on $\Gamma$ vanish. Let

$$ (\hat{\mathbf{u}}^H(t), \hat{\mathbf{v}}^H(t), \hat{\mathbf{u}}^h(t), \hat{\mathbf{v}}^h(t), \hat{\zeta}_1^h(t)) \quad (50) $$

be the time continuous solution of the differential algebraic equation (7) with the same mass matrix. Then, we obtain the following local error estimates for $j \in \{1, \dots, J\}$:

$$ \| \mathbf{u}_j^H - \hat{\mathbf{u}}^H(t_j) \| + \| \mathbf{v}_j^H - \hat{\mathbf{v}}^H(t_j) \| \le C \Delta T^3, $$

$$ \| \mathbf{u}_{\omega j}^h - \hat{\mathbf{u}}_\omega^h(t_j) \| + \| \mathbf{v}_{\omega j}^h - \hat{\mathbf{v}}_\omega^h(t_j) \| \le C j \Delta t^3, $$

$$ \| \mathbf{u}_{\Gamma j}^h - \hat{\mathbf{u}}_\Gamma^h(t_j) \| + \| \mathbf{v}_{\Gamma j}^h - \hat{\mathbf{v}}_\Gamma^h(t_j) \| \le C \frac{j(J \cdot j)}{j^2} \Delta T^2 + \mathcal{O}(\Delta T^3), $$

$$ \| \zeta_{\Gamma j}^h - \hat{\zeta}_1^h(t_{j/2}) \| \le C \Delta t^2. $$

The constants $C$ of the leading order terms can depend on $h$, the material parameters and the time derivatives of (50) but are independent of $J, j$ and $\Delta T$.
**Proof.** The idea of the proof is to consider the time continuous solution (50) as a perturbed solution of the discrete system (48) and (49) and to estimate the perturbation terms via Taylor expansion. For the details, we refer to [57]. □

Theorem 4 implies that the local errors at time $t_j$, i.e., after one coarse time step, are of order $O(\Delta T^3)$ for both displacement and velocity, and thus the time stepping scheme is globally second or der convergent. Further, we obtain optimal local estimates for the fine grid quantities in the interior of the subdomain $\omega$.

However, a direct solution of the coupled system (48) and (49) is very expensive, as it involves the simultaneous solution of one coarse problem on $\Xi$ and $J$ fine subproblems on $\omega$. Here, we benefit from the iterative algorithm derived in Section 3 which provides an efficient way of solving this system, which is the topic of the next subsection.

## 7.2. Approximate solution scheme

The idea of Algorithm 1 can be transferred to the coupled system (48) and (49). One iteration of the resulting scheme consists of the solution of a coarse problem on $V^H$ with the macro time step size $\Delta T$, followed by $J$ intermediate fine scale problems with step size $\Delta t$ and Dirichlet boundary conditions on $\Gamma$. The backcoupling for the next iteration is performed via the difference in the coarse and fine boundary stress at $\Gamma$. The corresponding numerical scheme for the iterative computation of

$$ \hat{\mathbf{z}}_j : \left( \mathbf{u}_{\Gamma j}^H, \boldsymbol{\mu}_{\Gamma j}^H, \mathbf{u}_{\Gamma 1}^h, \zeta_{\Gamma 1}^h, \dots, \mathbf{u}_{\Gamma J}^h, \zeta_{\Gamma J}^h \right), $$
given all values at time $t_0$, is stated in Algorithm 3.

**Algorithm 3.** Two way coupling scheme with time subcycling
($t_0 \to t_j$)

Starting from $\hat{\mathbf{z}}_j^{(0)}$, compute sequentially for $l = 0, 1, \dots$,
($l_{\max} = 1$)

(i) Solve problem on coarse space $V^H$ with interface load on $\Gamma$:

$$ \begin{pmatrix} (S_{\Xi}^H + S_{\omega}^H) & \mathbf{0} \\ S_{\omega}^H & (D_{\Gamma}^{HH})^T \end{pmatrix} \begin{pmatrix} \delta \mathbf{u}_{\Gamma j}^{H,(l)} \\ \delta \boldsymbol{\mu}_{\Gamma j}^{H,(l)} \end{pmatrix} = \begin{pmatrix} \mathbf{r}_{\Xi}^{H,(l)} + \mathbf{r}_{\omega}^{H,(l)} \\ \mathbf{r}_{\omega}^{H,(l)} \end{pmatrix}. \quad (51) $$

(ii) Loop on time substeps: $j = 1, \dots, J$: Solve problem on fine space $V^h$ with weakly imposed trace on $\Gamma$:

$$ \begin{pmatrix} S_{\omega}^{h} & (D_{\Gamma}^{hh})^{T} \\ D_{\Gamma}^{hh} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \delta \mathbf{u}_{\Gamma j}^{h,(l)} \\ \delta \zeta_{\Gamma j}^{h,(l)} \end{pmatrix} = \begin{pmatrix} \mathbf{r}_{\omega j}^{h,(l)} \\ \mathbf{v}_{\Gamma j}^{h,(l)} \end{pmatrix} + \mathbf{0}, \quad \sum_{j} D_{\Gamma}^{hh} \delta \mathbf{u}_{\Gamma j}^{H,(l)}. \quad (52) $$

(iii) Update the solution vector:

$$ \hat{\mathbf{z}}_j^{(l+1)} : \hat{\mathbf{z}}_j^{(l)} + \delta \hat{\mathbf{z}}_j^{(l)}. $$

The residuals of (51), (52) are given by

$$ \mathbf{r}_{\Xi}^{H,(l)} \qquad \rho_{\Xi 0}^{H} \qquad S_{\Xi}^{H} \mathbf{u}_{\Gamma j}^{H,(l)} + \frac{1}{J} \sum_{j=1}^{J} (D_{\Gamma}^{hh})^{T} \zeta_{\Gamma j}^{h,(l)}, \quad (53a) $$

$$ \mathbf{r}_{\omega}^{H,(l)} \qquad \rho_{\omega 0}^{H} \qquad S_{\omega}^{H} \mathbf{u}_{\Gamma j}^{H,(l)} (D_{\Gamma}^{HH})^{T} \boldsymbol{\mu}_{\Gamma j}^{H,(l)}, \quad (53b) $$

$$ \mathbf{r}_{\omega j}^{h,(l)} \qquad \rho_{\omega j 1}^{h,(l)} \qquad S_{\omega}^{h} \mathbf{u}_{\Gamma j}^{h,(l)} (D_{\Gamma}^{hh})^{T} \zeta_{\Gamma j}^{h,(l)}, \quad (53c) $$

$$ \mathbf{v}_{\Gamma j}^{h,(l)} \left( 1 - \frac{j}{J} \right) D_{\Gamma}^{hh} \mathbf{u}_{\Gamma 0}^{H} + \frac{j}{J} D_{\Gamma}^{hh} \mathbf{u}_{\Gamma j}^{H,(l)} \quad D_{\Gamma}^{hh} \mathbf{u}_{\Gamma j}^{h,(l)}. \quad (53d) $$

**Remark 17.** For $l \ge 1$, the residuals (53bcd) vanish.

An appropriate stopping criterion for Algorithm 3 can be derived as in Section 3.4 based on the residual

$$ \mathbf{r}_{\Xi}^{H,(l)} = (S_{\Xi}^{H} + S_{\omega}^{H}) \delta \mathbf{u}_{\Gamma j}^{H,(l)}. $$

This leads to the following relative algebraic error estimator for $l \ge 1$:

$$ (\eta_{J,\text{alg}}^{(l)})^2 = \frac{(\delta \mathbf{u}_{\Gamma j}^{H,(l)}, \mathbf{r}_{\Xi}^{H,(l)})}{(\mathbf{u}_{\Gamma j}^{H,(l)}, (S_{\Xi}^{H} + S_{\omega}^{H})\mathbf{u}_{\Gamma j}^{H,(l)})}. \quad (54) $$

## 8. Numerical results for time subcycling

In this section, we illustrate the performance of Algorithm 3 applied to two different two dimensional geometries.

### 8.1. Conforming geometry

First, we investigate the energy conservation and the time discretization error of the coupled system (48), (49) by means of a simple 2D test setting with $\Omega = (0,1) \times (0,3)$ and $\omega = (0,1) \times (0,1.5)$. We impose homogeneous Dirichlet boundary conditions on the top, the time dependent Neumann forces

$$ g_N(t) = \begin{cases} \frac{7}{8} \cdot 10^9 \max(0, \frac{1}{2} x_1) t & t \le 8 \cdot 10^5, \\ \frac{7}{8} \cdot 10^9 \max(0, \frac{1}{2} x_1)(1.6 \cdot 10^4 t) & 8 \cdot 10^5 < t \le 1.6 \cdot 10^4, \\ 0, & \text{else} \end{cases} $$

on the bottom and homogeneous Neumann conditions elsewhere. The volume load as well as the initial velocity are set to 0. In order to measure only the error due to time discretization, we use a uni form quadrilateral grid of mesh size $H = h = 2^3$. The parameters for the linear elastic material are chosen according to $v = 0.33$, $E_{\Xi} = E_{\omega} = 10^5$ and a discontinuous density of $\rho_{\Xi} = 1$, $\rho_{\omega} = 0.01$.

The left picture in Fig. 13 illustrates the evolution of the total energy (10) after each coarse time step. One can observe that after the excitation on the bottom of the body, i.e., for $t > 1.6 \cdot 10^4$, the energy is constant with respect to time, in agreement with the theoretical result of Lemma 3.

Next, we fix $\Delta T = 4 \cdot 10^5$ and investigate the evolution of the relative time discretization error on the patch $\omega$ given by

$$ e_{\omega, \text{time}}^2 : = \frac{(\mathbf{u}^h - \hat{\mathbf{u}}^h)^T A_{\omega}^h (\mathbf{u}^h - \hat{\mathbf{u}}^h)}{(\hat{\mathbf{u}}^h)^T A_{\omega}^h \hat{\mathbf{u}}^h}, \quad (55) $$

where the “reference solution” $\hat{\mathbf{u}}$ has been computed with an over all fine time step of size $2.5 \cdot 10^7$. The corresponding results for different values of $J$ are shown on the right side of Fig. 13. One can observe that the time subcycling decreases the local energy error in all cases.

For a numerical test of the dependence of the convergence rate of Algorithm 3 on $J$, we refer to [57].

### 8.2. Nonconforming geometry

Finally, we apply Algorithm 3 to the two dimensional geometry sketched in Figs. 1 and 15. The domain $\Omega$ consists of a circular ring centered at the origin with the diameters $r_{\text{inner}} = 1.6$, $r_{\text{outer}} = 1.95$ and 60 additional salients of height 0.05. The fine domain $\omega$ is built by several separate patches which are associated with the salients but have an extended T like shape in order to include the corner singularities in the fine triangulation. We remark that these patches can be computed in parallel as they do not overlap. Though
<figure><img src="image_14.png" /><figcaption>Fig. 13. Conforming time subcycling example: left: evolution of total energy for $J = 5$ and right: Evolution of local rel. error $e_{time,\omega}$ for $\Delta T = 4 \cdot 10^{-5}$, $J \in \{1,2,4,8,16\}$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. 2D two-scale example with $L = 2$ and $\Delta t = 2 \cdot 10^{-5}$; upper row: $J = 1$; lower row: $J = 10$; from left to right: effective stress at times $t_{50}$, $t_{150}$ and $t_{250}$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. 2D two-scale problem with time subcycling; left: evolution of active sets and right: sketch of grids and convergence statistics.</figcaption></figure>

fairly simple, the dynamic of the corresponding structure is consid-
ered a good paradigm for space and time two scale applications.

On the potential contact boundary $\gamma$, each salient features four additional small sipes which are resolved by the fine grid $T^h$ but are not respected by the coarse grid $T^H$ (see Fig. 15). Hence, similar to Section 6.2, we have a geometrically nonconforming situation with $V_{\omega}^h \subsetneq V_{\omega}^h$. However, the triangulation with $H/h = 4$ is con-structed such that the trace spaces $W_r^h \subset W_r^H$ are nested.

We set $\Delta t = 2 \cdot 10^{-5}$ and use the nonlinear Neo Hooke material (which corresponds to the Mooney Rivlin material law in (36) with $c_m = 0$) with the parameters $E_\infty = E_\infty = 4.4 \cdot 10^6$ and $\nu = 0.33$. The density on the inner ring is given by $Q_\infty = 1$, whereas the sali-ents have a density of $10^{-3}$. The structure has an initial velocity of

$$v_0(\mathbf{x}) = 150 \begin{pmatrix} x_2 \\ x_1 \end{pmatrix} \begin{pmatrix} 0 \\ 20 \end{pmatrix}$$

and an initial displacement corresponding to stationary rolling [58]. We set the volume load to zero and apply homogeneous Neumann boundary conditions everywhere except for the potential contact

boundary $\gamma$ where frictionless contact with a flat obstacle at $x_2 = 2.04$ is enforced. The tangential matrices are updated in each inner iteration ($l_{max} = 1$), and the resulting quasi Newton method is solved in each time step subject to $\|R^{(k)}\| \le 10^{-5}$. Further, due to the inertia terms, we do not use active coarse nodes for the coarse approximate problem.

In Fig. 14, the results for the computations with $J = 1$ and $J = 10$ are compared at different time instances, whereas the evolution of the active set is shown on the left of Fig. 15. Although one can ob-serve small differences, the overall quality of the results with time subcycling is evident. We remark that the computation with an overall coarse time step of size $\Delta T = 2 \cdot 10^{-4}$ has not converged within 30 Newton iterations in the second and all subsequent time steps.

The total number of Newton iterations as well as the asymptotic error reduction factor are summarized on the right of Fig. 15. We emphasize that each Newton iteration consists of the solution of one global coarse and $J$ fine local systems. In comparison with $J = 1$, the number of fine solves has increased by 149% for the computation with $J = 10$, but the number of coarse global problems
has decreased by 75%. Hence, the subcycling scheme is extremely useful if the coarse global problems are much more expensive than the local solves.

## 9. Conclusion

In this paper, we have developed an efficient iterative algorithm for the solution of frictional dynamic contact problems with locally refined geometries. For this, we have employed an overlapping do main decomposition with an independent fine mesh around the contact zone. For linear problems, we have shown that under some reasonable assumptions, the convergence rate of the correspond ing algorithm is bounded independently of the material and dis cretization parameters. Further, we have investigated how the scheme can be extended to nonlinear contact problems as an inex act inner solver within the semismooth Newton iteration. Finally, we have shown that the overlapping discretization can be ex tended to a locally finer time scale as well as that the resulting cou pled problem is energy conserving in the linear case and can be solved efficiently by an iterative procedure.

## Appendix A. Proof of Theorem 2

First, we need to introduce some notation. We use the abbreviation $||_{k,\Theta} : ||_{k^h(\Theta)}$ and denote the space $\mathcal{R}B_{\omega}$ of rigid body modes on $\omega$ as well as the corresponding seminorm by

$$ \mathcal{R}B_{\omega} : \{ \mathbf{z} \in \mathbf{V}_{\omega} : a_{\omega}(\mathbf{z}, \mathbf{z}) = 0 \}, \quad |\mathbf{v}|_{\tilde{\mathbf{v}}_{\omega}} : \inf_{\mathbf{z} \in \mathcal{R}B_{\omega}} \|\mathbf{v} + \mathbf{z}\|_{1,\omega}, \mathbf{v} \in \mathbf{V}_{\omega}. $$

Further, according to (11a), we define the bilinear form

$$ k_{\omega}(\mathbf{v}, \mathbf{w}) : \frac{2}{\Delta t^2} m_{\omega}(\mathbf{v}, \mathbf{w}) + a_{\omega}(\mathbf{v}, \mathbf{w}), \quad \mathbf{v}, \mathbf{w} \in \mathbf{V}_{\omega}. $$

Due to the fact that the Schur complement matrices $S_{\omega}^H, S_{\omega}^{H\mathrm{H}}$ de fined in (16) and (21) correspond to discrete harmonic extensions onto $\omega$ with respect to $k_{\omega}(\cdot, \cdot)$, we obtain

$$ \begin{aligned} & (\mathbf{w}_I^H, S_{\omega}^{H\mathrm{H}} \mathbf{w}_I^H) \quad \inf_{\substack{\mathbf{v}^h \in \mathbf{V}^h \\ \mathbf{v}^h|_I \stackrel{p^h}{=} \mathbf{w}_I^H}} k_{\omega}(\mathbf{v}^h, \mathbf{v}^h), \\ & (\mathbf{w}_I^H, S_{\omega}^H \mathbf{w}_I^H) \quad \inf_{\substack{\mathbf{v}^H \in \mathbf{V}^H \\ \mathbf{v}^H|_I \stackrel{p^H}{=} \mathbf{w}_I^H}} k_{\omega}(\mathbf{v}^H, \mathbf{v}^H). \end{aligned} \quad \text{(A.1)} $$

As the spaces $\mathbf{V}^h, \mathbf{V}_{\omega}^H$ are finite dimensional, we can define the values

$$ \hat{\mathbf{v}}^h : \arg \inf_{\substack{\mathbf{v}^h \in \mathbf{V}^h \\ \mathbf{v}^h|_I \stackrel{p^h}{=} \mathbf{w}_I^H}} k_{\omega}(\mathbf{v}^h, \mathbf{v}^h), \quad \hat{\mathbf{v}}^H : \arg \inf_{\substack{\mathbf{v}^H \in \mathbf{V}^H \\ \mathbf{v}^H|_I \stackrel{p^H}{=} \mathbf{w}_I^H}} k_{\omega}(\mathbf{v}^H, \mathbf{v}^H). $$

The idea of the proof is to project the discrete function $\hat{\mathbf{v}}^h \in \mathbf{V}_{\omega}^h$ onto the space $\mathbf{V}_{\omega}^H$ and vice versa, and to construct an upper bound of the terms in (A.1) from these projections. A suitable projection operator $Z_{\omega}^m : \mathbf{V}_{\omega} \to \mathbf{V}_{\omega}^m$, $m \in \{h, H\}$, can be constructed similar to the Scott Zhang operator $Z_{\omega}^h$ defined in [59]. The only difference is that we change the original definition of [59] near the interface $\Gamma$ such that the relation $(Z_{\omega}^m \mathbf{v})|_{\Gamma}, P^m(\mathbf{v})|_{\Gamma}$ holds for $\mathbf{v} \in \mathbf{V}_{\omega}$, where the mortar projection $P^m$ has been defined in (5). With this, the operator $Z_{\omega}^m$ satisfies the approximation properties

$$ 0 \le l \le k \le 2, \quad k > \frac{1}{2} : \quad \|Z^h \mathbf{v} \|_{L,\omega} \le C h^k |\mathbf{v}|_{k,\omega}, \quad \mathbf{v} \in \mathbf{H}^k(\omega), \quad \text{(A.2)} $$

as well as the $\mathbf{H}^1(\omega)$ stability estimate

$$ |Z^h \mathbf{v}|_{1,\omega} \le C |\mathbf{v}|_{1,\omega}, \quad \mathbf{v} \in \mathbf{H}^1(\omega). \quad \text{(A.3)} $$

However, a uniform $\mathbf{L}_2(\omega)$ stability similar to (A.3) cannot be valid in general. As a counterexample, we look at the case that $\mathbf{v}^h$ is the extension by zero of $\mathbf{w}_I^h \in \mathbf{W}_I^h$ where we obtain

$$ \|Z_{\omega}^H \mathbf{v}^h\|_{0,\omega}^2 \le CH \|P^H \mathbf{w}_I^H\|_{0,\omega}^2 \le CH \|\mathbf{w}_I^H\|_{0,\omega}^2 \le C \frac{H}{h} \|\mathbf{v}^h\|_{0,\omega}^2. $$

Nevertheless, using (A.2) and an inverse estimate, the projection of $\mathbf{v}^h \in \mathbf{V}_{\omega}^H$ onto $\mathbf{V}^h$ gives the stability estimate

$$ \begin{aligned} \|Z_{\omega}^h \mathbf{v}^h\|_{0,\omega} &\le \|Z_{\omega}^h \mathbf{v}^H\|_{0,\omega} + \|\mathbf{v}^H\|_{0,\omega} \le Ch \|\mathbf{v}^H\|_{1,\omega} + \|\mathbf{v}^H\|_{0,\omega} \\ &\le (1 + C \frac{h}{H}) \|\mathbf{v}^H\|_{0,\omega} \le C \|\mathbf{v}^H\|_{0,\omega}. \end{aligned} \quad \text{(A.4)} $$

With these preliminaries, we can show the second inequality of (26), using the stability estimates (A.3), (A.4), the equation $(Z_{\omega}^h \mathbf{v}^H)|_{\Gamma}, P^h \mathbf{w}_I^H$ as well as the fact that $Z_{\omega}^h \mathbf{z}, \mathbf{z}$ holds for $\mathbf{z} \in \mathcal{R}B_{\omega} \subset (\mathbf{V}_{\omega}^H \cap \mathbf{V}_{\omega}^H):$

$$ \begin{aligned} k_{\omega}(\hat{\mathbf{v}}^h, \hat{\mathbf{v}}^h) &\le k_{\omega}(Z_{\omega}^h \hat{\mathbf{v}}^H, Z_{\omega}^h \hat{\mathbf{v}}^H) \le C \left( E_{\omega} \left|Z_{\omega}^h \hat{\mathbf{v}}^H\right|_{\tilde{\mathbf{v}}_{\omega}}^{\frac{2}{h}} + \frac{Q_{\omega}}{\Delta t^2} \|Z_{\omega}^h \hat{\mathbf{v}}^H\|_{0,\omega}^2 \right) \\ &\le C \left( E_{\omega} \left|\hat{\mathbf{v}}^H\right|_{\tilde{\mathbf{v}}_{\omega}}^{\frac{2}{h}} + \frac{Q_{\omega}}{\Delta t^2} \|\hat{\mathbf{v}}^H\|_{0,\omega}^2 \right) \le C^* k_{\omega}(\hat{\mathbf{v}}^H, \hat{\mathbf{v}}^H). \end{aligned} $$

For the other estimate, we make use of the nested trace spaces $\mathbf{W}_I^H \subset \mathbf{W}_I^h$, implying $P^h P^h \mathbf{w}_I^H, \mathbf{w}_I^H$. For any fixed rigid body motion $\mathbf{z} \in \mathcal{R}B_{\omega}$, we thus obtain

$$ \begin{aligned} k_{\omega}(\hat{\mathbf{v}}^H, \hat{\mathbf{v}}^H) &\le C \left( E_{\omega} \left|Z_{\omega}^H \hat{\mathbf{v}}^H\right|_{\tilde{\mathbf{v}}_{\omega}}^{\frac{2}{h}} + \frac{Q_{\omega}}{\Delta t^2} \|Z_{\omega}^H \hat{\mathbf{v}}^H\|_{0,\omega}^2 \right) \\ &\le C \left( E_{\omega} \left|\hat{\mathbf{v}}^H\right|_{\tilde{\mathbf{v}}_{\omega}}^{\frac{2}{h}} + \frac{Q_{\omega}}{\Delta t^2} \left( \|(Z_{\omega}^h \mathrm{Id})(\hat{\mathbf{v}}^H + \mathbf{z})\|_{0,\omega}^2 + \|\hat{\mathbf{v}}^H\|_{0,\omega}^2 \right) \right) \\ &\le C \left( E_{\omega} \left|\hat{\mathbf{v}}^H\right|_{\tilde{\mathbf{v}}_{\omega}}^{\frac{2}{h}} + \frac{Q_{\omega}}{\Delta t^2} \left( H^2 |\hat{\mathbf{v}}^H + \mathbf{z}|_{1,\omega}^2 + \|\hat{\mathbf{v}}^H\|_{0,\omega}^2 \right) \right), \end{aligned} $$

where we have used the estimates (A.2) and (A.3). Minimizing on $\mathbf{z} \in \mathcal{R}B_{\omega}$ and using Assumption (25), we get the upper bound

$$ \begin{aligned} k_{\omega}(\hat{\mathbf{v}}^H, \hat{\mathbf{v}}^H) &\le C(1 + C_{\mathrm{CFL}})E_{\omega} \|\hat{\mathbf{v}}^H\|_{\tilde{\mathbf{v}}_{\omega}}^{\frac{2}{h}} + C \frac{Q_{\omega}}{\Delta t^2} \|\hat{\mathbf{v}}^H\|_{0,\omega}^2 \\ &\le (c^*)^1 (1 + C_{\mathrm{CFL}}) k_{\omega}(\hat{\mathbf{v}}^H, \hat{\mathbf{v}}^H), \end{aligned} $$

implying (26).

Finally, we prove that the constant $C^*$ in (26) can be replaced by one if the finite element spaces $\mathbf{V}_{\omega}^H \subset \mathbf{V}^h$ are nested. For this, we need to show the inequality

$$ \inf_{\substack{\mathbf{w}^h \in \mathbf{W}^h \\ \mathbf{w}^h|_I \stackrel{p^h}{=} \mathbf{w}_I^H}} k_{\omega}(\mathbf{w}^H, \mathbf{w}^H) \quad \inf_{\substack{\mathbf{w}^h \in \mathbf{V}^h \\ \mathbf{w}^h|_I \stackrel{p^h}{=} \mathbf{w}_I^H}} k_{\omega}(\mathbf{w}^h, \mathbf{w}^h) \ge 0, $$

which follows from the nestedness of $\mathbf{W}_I^H \subset \mathbf{W}_I^h$ and

$$ \{\mathbf{w}^h \in \mathbf{V}_{\omega}^H : \mathbf{w}^h|_{\Gamma} \quad \mathbf{w}_I^H\} \subset \{\mathbf{w}^h \in \mathbf{V}^h : \mathbf{w}^h|_{\Gamma} \quad P^h \mathbf{w}_I^H\}. $$

## Appendix B. Linearization of contact conditions

The linearization of the NCP function (42) is given in terms of the (in) active sets $\mathcal{A}_n^{h,(k)}, \mathcal{I}_n^{h,(k)}$ and $\mathcal{A}_t^{h,(k)}, \mathcal{I}_t^{h,(k)}$ which are defined according to (41) but evaluated at the $k$th Newton iterate $(\mathbf{u}_t^{h}, \dot{\mathbf{u}}_t^{h})^{(k)}$. For the case of Tresca friction with $\tilde{\mathbf{g}} = 0$, it is implemented in the following way, omitting the indices $p, j$ and $h$ for ease of notation:

$$ \left\{ \begin{array}{ll} \lambda_n^{(k+1)} & 0, \quad p \in \mathcal{I}_n^{(k)}, \\ u_n^{(k+1)} & g_n, \quad p \in \mathcal{A}_n^{(k)}, \\ \mathbf{u}_t^{(k+1)} & \mathbf{u}_j \quad 1, t, \quad p \in \mathcal{I}_t^{(k)}, \\ \lambda_t^{(k+1)} & L_t^{(k)} \mathbf{u}_t^{(k+1)} \quad \mathbf{q}_t^{(k)}, \quad p \in \mathcal{A}_t^{(k)}. \end{array} \right. $$
Above, we have used the abbreviations

$$L_t: \frac{c_t g_t}{\|\lambda_t^{\text{tr}}\|} g_t \left( \text{Id} - \frac{a \|\lambda_t^{\text{tr}}\|}{g_t \xi_t} \lambda_t \otimes \lambda_t^{\text{tr}} \right),$$

$$\mathbf{q}_t: \frac{a \|\lambda_t^{\text{tr}}\|^2}{\xi_t} \lambda_t,$$

$$\xi_t: \|\lambda_t^{\text{tr}}\|^2 g_t \|\lambda_t^{\text{tr}}\| + a \lambda_t: \lambda_t^{\text{tr}}.$$

For the correct derivative of (42), the value of a above is one. But for stability reasons, we implement a quasi Newton scheme with $a \min(1, \frac{g_t}{\|\lambda_t^{\text{tr}}\|})$. We refer to [16,54] for details as well as for the case $\mathcal{R} > 0$.

## References

[1] F. Ben Belgacem, The mortar finite element method with Lagrange multipliers, Numer. Math. 84 (1999) 173-197.

[2] C. Bernardi, Y. Maday, A.T. Patera, A new nonconforming approach to domain decomposition: The mortar element method, in: H. Brezis, J.-L. Lions (Eds.), Nonlinear Partial Differential Equations and Their Applications, Collège de France Seminar, Vol. XI, Pitman, 1994, pp. 13-51.

[3] B.I. Wohlmuth, A mortar finite element method using dual spaces for the Lagrange multiplier, SIAM J. Numer. Anal. 38 (2000) 989-1012.

[4] A. Quarteroni, A. Valli, Domain Decomposition Methods for Partial Differential Equations, Oxford University Press, 1999.

[5] B. Smith, P.E. Bjerstad, W. Gropp, Domain decomposition, in: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge Univ. Press, 1996.

[6] A. Toselli, O.B. Widlund, Domain Decomposition Methods – Algorithms and Theory, Springer, Berlin Heidelberg, 2005.

[7] P.E. Bjerstad, O.B. Widlund, Iterative methods for the solution of elliptic problems on regions partitioned into substructures, SIAM J. Numer. Anal. 23 (1986) 1097-1120.

[8] P. Goldfeld, L.F. Pavarino, O.B. Widlund, Balancing Neumann-Neumann preconditioners for mixed approximations of heterogeneous problems in linear elasticity, Numer. Math. 95 (2003) 283-324.

[9] A. Klawonn, O.B. Widlund, FETI and Neumann-Neumann iterative substructuring methods: Connections and new results, Comm. Pure Appl. Math. 54 (2001) 57-90.

[10] C. Eck, J. Jarušek, M. Krbec, Unilateral Contact Problems: Variational Methods and Existence Theorems, CRC Press Taylor & Francis Group, Boca Raton, 2005.

[11] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM Studies in Applied Mathematics, vol. 8, Society for Industrial Mathematics, Philadelphia, 1988.

[12] T.A. Laursen, Computational Contact and Impact Mechanics, Springer, Berlin, Heidelberg, 2002.

[13] P. Wriggers, Computational Contact Mechanics, J. Wiley & Sons Ltd., 2002.

[14] B.I. Wohlmuth, Variationally consistent discretization schemes and numerical algorithms for contact problems, Acta Numer. 20 (2011) 569-734.

[15] M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semi-smooth Newton method, SIAM J. Optim. 13 (2003) 865-888.

[16] S. Hüber, G. Stadler, B.I. Wohlmuth, A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction, SIAM J. Sci. Comput. 30 (2008) 572-596.

[17] S. Hüber, B.I. Wohlmuth, A primal-dual active set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg. 194 (2005) 3147-3166.

[18] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Comput. Methods Appl. Mech. Engrg. 92 (1991) 353-375.

[19] P.W. Christensen, A. Klarbring, J.S. Pang, N. Strömberg, Formulation and comparison of algorithms for frictional contact problems, Int. J. Numer. Methods Engrg. 42 (1998) 145-173.

[20] R.S. Dembo, S.C. Eisenstat, T. Steinhaug, Inexact Newton methods, SIAM J. Numer. Anal. 19 (1982) 400-408.

[21] P. Deuflhard, Newton Methods for Nonlinear Problems, Springer, Berlin, Heidelberg, 2004.

[22] S.C. Eisenstat, H.F. Walker, Choosing the forcing terms in an inexact Newton method, SIAM J. Sci. Comput. 17 (1996) 16-32.

[23] J. Simo, N. Tarnow, The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics, Z. Angew. Math. Phys. 43 (1992) 757-792.

[24] P. Betsch, P. Steinmann, Conservation properties of a time FE method – Part II: Time-stepping schemes for non-linear elastodynamics, Internat. J. Numer. Methods Engrg. 50 (2001) 1931-1955.

[25] O. Gonzalez, Exact energy and momentum conserving algorithms for general models in nonlinear elasticity, Comp. Methods Appl. Mech. Engrg. 190 (2000) 1763-1783.

[26] P. Hauret, P. Le Tallec, Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact, Comput. Methods Appl. Mech. Engrg. 195 (2006) 4890-4916.

[27] V. Chawla, T.A. Laursen, Energy consistent algorithms for frictional contact problems, Internat. J. Numer. Methods Engrg. 42 (1998) 799-827.

[28] C. Hager, S. Hüber, B.I. Wohlmuth, A stable energy conserving approach for frictional contact problems based on quadrature formulas, Int. J. Numer. Methods Engrg. 73 (2008) 205-225.

[29] C. Hager, B.I. Wohlmuth, Analysis of a space-time discretization for dynamic elasticity problems based on mass-free surface elements, SIAM J. Numer. Anal. 47 (2009) 1863-1885.

[30] H.B. Khenous, P. Laborde, Y. Renard, Mass redistribution method for finite element contact problems in elastodynamics, Eur. J. Mech. A/Solids 27 (2008) 918-932.

[31] W.J.T. Daniel, Subcycling first- and second-order generalizations of the trapezoidal rule, Int. J. Numer. Methods Engrg. 42 (1998) 1091-1119.

[32] W.J.T. Daniel, A study of the stability of subcycling algorithms in structural dynamics, Comput. Methods Appl. Mech. Engrg. 156 (1998) 1-13.

[33] P. Smolinski, S. Sleigh, T. Belytschko, Stability of an explicit multi-time step integration algorithm for linear structural dynamics equations, Comput. Mech. 18 (1996) 236-244.

[34] A. Combescure, A. Gravouil, A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis, Comput. Methods Appl. Mech. Engrg. 191 (2002) 1129-1157.

[35] A. Gravouil, A. Combescure, Multi-time-step explicit-implicit method for non-linear structural dynamics, Int. J. Numer. Methods Engrg. 50 (2001) 199-225.

[36] A. Prakash, K.D. Hjelmstad, A FETI-based multi-time-step coupling method for Newmark schemes in structural dynamics, Int. J. Numer. Methods Engrg. 61 (2004) 2183-2204.

[37] N. Mahjoubi, S. Krenk, Multi-time-step domain coupling method with energy control, Int. J. Numer. Methods Engrg. 83 (2010) 1700-1718.

[38] P.G. Citart, Mathematical Elasticity, vol. I, North Holland, 1998.

[39] C. Hager, P. Hauret, P. Le Tallec, B.I. Wohlmuth, Overlapping domain decomposition for multiscale dynamic contact problems, 2010, IANS Preprint 2010/007, Universität Stuttgart.

[40] D. Braess, Finite Elements, Cambridge University Press, 1997.

[41] S.C. Brenner, Korn's inequalities for piecewise $H^1$ vector fields, Math. Comput. 73 (2003) 1067-1087.

[42] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Springer, Berlin, Heidelberg, 1991.

[43] B.I. Wohlmuth, Discretization Methods and Iterative Solvers based on Domain Decomposition, Springer, Berlin, 2001.

[44] B. Flemisch, B.I. Wohlmuth, Stable Lagrange multipliers for quadrilateral meshes of curved interfaces in 3D, Comput. Methods Appl. Mech. Engrg. 196 (2007) 1589-1602.

[45] T.J.R. Hughes, The Finite Element Method: Linear, Static and Dynamic Finite Element Analysis, Prentice-Hall, 1987.

[46] P. Hauret, P. Le Tallec, Two-scale Dirichlet-Neumann preconditioners for elastic problems with boundary refinements, Comput. Methods Appl. Mech. Engrg. 196 (2007) 1574-1588.

[47] X. Xu, L. Qin, Spectral analysis of Dirichlet-Neumann operators and optimized Schwarz methods with Robin transmission conditions, SIAM J. Numer. Anal. 47 (2010) 4540-4568.

[48] P. Le Tallec, Numerical methods for nonlinear three-dimensional elasticity, Handbook of Numerical Analysis, vol. 3, North Holland, 1994, pp. 465-622.

[49] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Springer, Berlin, Heidelberg, 1998.

[50] F. Facchinei, J.-S. Pang, Finite-Dimensional Variational Inequalities and Computational Problems, Springer, New York, 2003.

[51] H.-H. Dai, Model equations for nonlinear dispersive waves in a compressible Mooney-Rivlin rod, Acta Mech. 127 (1998) 193-207.

[52] S. Hüber, B.I. Wohlmuth, An optimal a priori error estimate for non-linear multibody contact problems, SIAM J. Numer. Anal. 43 (2005) 157-173.

[53] P. Wriggers, Nichtlineare Finite-Element-Methoden, Springer, Berlin, Heidelberg, 2001.

[54] C. Hager, B.I. Wohlmuth, Nonlinear complementarity functions for plasticity problems with frictional contact, Comput. Methods Appl. Mech. Engrg. 198 (2009) 3411-3427.

[55] T.A. Laursen, G.R. Love, Improved implicit integrators for transient impact problems – geometric admissibility within the conserving framework, Internat. J. Numer. Methods Engrg. 53 (2002) 245-274.

[56] R. Krause, A nonsmooth multiscale method for solving frictional contact problems in 2D and 3D with multigrid efficiency, SIAM J. Sci. Comput. 31 (2009) 1399-1423.

[57] C. Hager, Robust Numerical Algorithms for Dynamic Frictional Contact Problems with Different Time and Space Scales, Ph.D. Thesis, Universität Stuttgart, 2010.

[58] P. Hauret, J. Salomon, A. Weiß, B.I. Wohlmuth, Energy consistent co-rotational schemes for frictional contact problems, SIAM J. Sci. Comput. 30 (2008) 2488-2511.

[59] L.R. Scott, S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comput. 54 (1990) 483-493.