Solving dynamic contact problems with local refinement in space and time

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.


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 simula tion, 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 roll ing contact. More importantly, many other multiscale contact sim ulations 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 consis tent description of the dynamics. The aim of this work is to design such an algorithm by combining several well established mathe matical methods with some new approaches. All of them are de scribed 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 {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 consis 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 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 vari ationally 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.

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.

Problem statement
In the following, we consider an elastic body X & R d , where d 2 {2, 3} denotes the number of spatial dimensions. The polyhe dral boundary oX is partitioned into two nonoverlapping parts C D , C N with meas(C D ) > 0. On C D , we assume homogeneous Dirich let boundary conditions for simplicity, whereas a surface load de noted by g N acts on the boundary C N where the unit outer normal n is defined almost everywhere. Then, the strong form of the dynamic linear elasticity problem on X subject to the volume load l reads (see, e.g., [38]) . € u div rðuÞ l in X; Here, the Cauchy stress r = r(u) is given by r : C el eðuÞ with the linearized strain tensor eðuÞ : 1 2 ru þ ru T À Á . The Hooke tensor C el C el ðxÞ can be defined either in terms of the Lamé parameters k, l or with respect to the Young modulus E and the Poisson coeffi cient m. Further, . P 0 denotes the density of X.
We consider the case for which there exists a given subregion x & X, with possibly different material parameters, where a more accurate local resolution is desired. In the situation of Fig. 1 [39].
In order to obtain the variational form of (1) Further, let W C : H 1/2 (C) be the space of traces on C. Then, the weak form of (1) can be formulated as two separate problems on the subdomains N, x, where the continuity of the displacements across C is enforced weakly by means of a Lagrange multiplier f C 2 M C : W 0 C which is the dual space of W C . This leads to the fol lowing coupled problem: plus appropriate initial conditions uj t 0 u 0 ; _ uj t 0 v 0 . Above, we have used the L 2 duality pairing h , i C on C. For the quasi static case, the well posedness of (2) follows from the condition meas (C D ) > 0 and Korn's inequality [40,41], whereas for the dynamic case, we need that both initial values u 0 , v 0 satisfy condition (2c) [42].

Discretization
For the spatial discretization of (2), we first triangulate the glo bal domain X in terms of a quasi uniform regular mesh T H of sim plicial or quadrilateral/ hexahedral elements of size H. We assume that the Dirichlet boundary C D as well as the interface C are re solved by this triangulation. On T H , we consider lowest order con With these triangulations, we define the following vector val ued finite element spaces associated with the coarse and the fine discretization: On the interface C, we define the discrete trace spaces C & M C spanned by the basis functions fw m p g p2N m C ; m 2 fh; Hg. 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 fw m p g p2N m C which are piecewise (bi) linear, discontinuous and satisfy the locality requirement suppðw m p Þ suppð/ m p j C Þ as well as the biorthogonality property As shown in [44], these basis functions can efficiently be con structed even for general grids in 3D. We remark that the above def inition of the multiplier spaces yields dimðM m C Þ dimðW m C Þ, i.e., we do not modify the basis functions of M m C near the locations where C is nondifferentiable.
In the rest of this subsection, we formulate the discrete alge braic version of (2). For this, we introduce the mortar projection onto W m C ; m 2 fh; Hg, given by As shown in [43], the mortar operator (5) is uniformly continuous in h or H with respect to the L 2 (C) , the H 1 (C) and the H 1/2 (C) norm. We define the matrices D lm with Id denoting the d Â d identity matrix. With this, the algebraic representation of the operator P l j W m C : W m C ! W l C is given by Remark 1. Due to the biorthogonality property (4), the matrices D ll C 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 lh C ; l 2 fh; Hg, by zero to the vector space V h is named D lh 2 R djN l C jÂdjN h j ; similarly, the extension of D lH C to the coarse upper space V H N is denoted by D lH . Using this space discretization, standard matrix notation for the mass and stiffness matrices as well as the fine Lagrange multiplier space M h C for the continuity constraints, the spatially discrete ver sion of (2) reads together with suitable initial conditions  (7) results in a standard conforming scheme.
Next, we discretize (7) in time by partitioning the time interval (0, T) into equidistant time steps Dt with discrete times t j = jDt and introducing the notation Using the implicit trapezoidal rule, i.e., the Newmark scheme with The discrete energy at time t j is the sum of the contributions from the subdomains N, x given by For time independent outer loads, the energy is conserved for the discrete system (9), i.e., E j E j 1 holds, provided that (9e) is valid for both time steps t j and t j 1 . This can easily be verified by multi plying (9ac) with the mean velocities v m j 1=2 as test functions, adding the results and using (9bde).
Using (9bd) to eliminate the velocities from (9ac) and defining . h as well as analogous expressions for the quantities on N, the system to be solved for the displacement at time t j becomes

Schur complement formulation
The system (12) can equivalently be expressed in terms of the unknowns on the interface C only. For this, we partition the vec tors (u H , u h ) and the corresponding matrices K H N ; K h x into compo nents associated with the interface nodes and the inner degrees of freedom: Performing static condensation of the inner variables and introduc ing the symmetric positive semi definite Schur complement matri ces S H N ; S h x given by we can rewrite (9) as where the right hand side of (14) 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 V H . In Section 3, the corresponding algorithm is derived and theo retically analyzed. Numerical results are presented in Section 4.

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 X and a subse quent solution of a local fine grid problem on x. 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.

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 The distinctive feature of the algorithm is that is contains a global coarse problem on X subject to an interface load on C, whose solu tion defines a coarse grid approximation of the fine grid quantities in x. The details of the scheme are summarized in Algorithm 1.
Remark 3. After one iteration of Algorithm 1, i.e., for l P 1, the residuals (20bcd) vanish. x and S h x 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ẑ ð0Þ , compute sequentially for (ii) Solve problem on local fine space V h with weakly imposed trace on C inherited from coarse computation on X:

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 u H Cj 2 R djN H C j . Introducing the matrix S HhH (14) can be rewritten as The following lemma shows that Algorithm 1 can be reformu lated as a fixed point iteration on u H Cj : Proof. The result follows from the definition of the iterative scheme by substitution. Details can be found in [39]. h Remark 5. If the coarse and the fine grid on x coincide, i.e., then the exact solution is already obtained after the first iteration. A damped version of (23) with a suitable damping parameter a l > 0 yields the error propagation e H;ðlþ1Þ If a l = a 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 a l in each step according to a conjugate gradient algorithm (see, e.g., [6,7]).
. 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 x ) E N or . x ) . N , i.e., if the local details are much more stiff than the rest of the domain (cf. [7,39]).

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 H x and S HhH x in the following theorem proved in Appendix A: Then, there exist constants c ⁄ , C ⁄ independent of the diameter of x, h, H, Dt and the values E x , . x , such that the following estimates are satis fied for any w H If the finite element spaces (26) can be replaced by one.

Remark 7.
The assumption W H C & W h C 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 condi tion number of the iteration matrix in (23) is bounded by with the constants c ⁄ , C ⁄ from Theorem 2. Furthermore, the iterates of (24) converge to the solution u H Cj of (22) if the damping parameter a l is chosen within the set 0 < a min 6 a l < 2 C Ã for some fixed value of a min (see, e.g., [6]). Hence, if the finite element spaces V H x & V h are nested, Theorem 2 implies the convergence of (24) for a 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 D t 2 6 Ch 2 , where C is some constant depending on the maximum value of the ratio E/.. 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 logarith mic terms) to a lower bound of (26). The resulting dependence of the convergence rate of Algorithm 1 on the ratio H/h can be observed from the numerical re sults in Section 4.3; furthermore, it is in agreement with the theo retical estimates of the spectrum of discrete Dirichlet to Neumann operators recently presented in [47].

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 P 1, the only nonzero com ponent of the residual vector is r H;ðlÞ N given in (20a where the last equality follows from (18) for l P 1. Thus, we pro pose to use the following relative algebraic error estimator for l P 1: 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  (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.

Algebraic error for static case
For the first set of tests, we set . = 0, m = 0.3, and piecewise con stant Young moduli E N 100; E x 10 par E N ; par 2 Z. Further, we enforce homogeneous Dirichlet conditions on the upper boundary, homogeneous Neumann boundary conditions on the left and right  In order to investigate the decrease of the algebraic error for Algorithm 1 with respect to the number of iterations, the differ ence between the values (u H,(l) , u h,(l) ) obtained by Algorithm 1 and the solution (u H , u h ) of the discrete mortar system (14) is mea sured with respect to the relative energy norm. For this, we compute e ðlÞ alg 2 : 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 2 {1, 2, 3, 4}. The left picture shows the results for equal material parameters E x = E N , whereas the right picture displays the convergence for discontinuous param eters E x = 10 5 E N . 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 g 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 x = 10 par E N is further investigated in Fig. 5, where the error reduction fac tor e ðlþ1Þ alg =e ðlÞ alg is plotted with respect to l. The results show that the algorithm converges best for small values of L and par, but the reduction factor is limited from above by around 0.35 for par = 0 independently of L > 0 and by around 0.55 for L = 2 and par ? 1. Hence, in the static case, Algorithm 1 is stable with respect to jumps in the coefficients, in accordance with the theoretical results of Theorem 2 for the case . x = 0.

Algebraic error decrease for dynamic case
Next, we take the inertia terms into account and consider one time step of a dynamic setting with Dt = 10 3 and piecewise con stant material parameters E N 100; . N 0:01; E x 10 par E E N ; . x 10 par . . N . The parameters are chosen such that the minimal value of the constant C CFL in (25) 9] is chosen such that the deformation of the body is within a sensible range. Both the vol ume load and the initial velocity are set to zero. In Table 1, the algebraic error reduction rates of Algorithm 1 are summarized for L 2 {2, 4}, par E , par . 2 { 4, 2, 0, 2, 4, 6}. The results illustrate the influence of the constant C CFL in (25) and the bound in Remark 8. As long as par . 6 par E + 2 and thus C CFL 6 200, the error reduction factor is bounded from above independently of the ratio of the grid sizes or the material parameters E x and . x , correspond ing to Theorem 2. However, for par . > par E + 4, we get a reduction factor of (1 h/H), as predicted by Remark 8.

Extension to nonlinear problems
Motivated by the challenge of simulating the large deformation dynamics of structures with contact, we have presented and ana lyzed 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 decom position approach. Section 5.2 treats the case of dynamic frictional contact as an example of a nonlinear effect.

Iterative solution scheme
In comparison with the mortar system (9)   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 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 Þ, but we omit the dependence on the latter two arguments for ease of notation. Further, the exact form of K m H depends on the kind of nonlinearity considered. In Sub section 5.2, we state the conditions of frictional contact as an exam ple 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]  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ẑ introduced in (17) 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.

Frictional contact
In this subsection, we exemplify the formation of the nonlinear semismooth operators K m H 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 m H proceeds along the lines of Section 2.1, but with the Cauchy stress tensor r in (1) being replaced by the first Piola Kirchhoff stress tensor P(u). The latter is given by P(u) = (Id + ru)S(u), with S(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 S is defined in terms of the Lamé parameters l, k, the additional parameter c m and the right Cauchy Green tensor Due to the nonlinear relationship between u and P, the stiffness terms in (7) are now nonlinear (i.e., A H N ðu H Þ and A h x ð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 c & C N of the Neumann boundary where the body X can come into contact with a fixed smooth obstacle C obs ; for the contact of several elastic bodies we refer to, e.g., [17,52]. On c, we enforce the Neumann conditions Pn = k, where the contact stress k 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 x is large en ough to cover the whole contact boundary, i.e., 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 2 N h c (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 u h p into its normal and tan gential part according to u h pn : u h p n p ; u h pt : u h p u h pn n p , where n p denotes a suitable approximation of the unit outer normal of x at p 2 N h c . Further, we introduce the initial normal gap g pn be tween the node p and the fixed obstacle C obs in the direction of n p (see, e.g., [12,53] for a detailed definition), as well as the normal and tangential parts k h pn : k h p n p ; k h pt : k h p k h pn n p of the Lagrange multiplier. With this, the constraints in the normal direction at the time step t j read u h pjn g pn 6 0; k h pjn P 0; k h pjn ðu h pjn g pn Þ 0: ð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 kk h pjt k 6 g pt þ Fjk h pjn j; kk h pjt k < g pt þ Fjk h pjn j ) @ Dt u h pjt 0; If the friction coefficient F P 0 is set to zero, the above inequalities describe the case of Tresca friction with the given friction bound g pt P 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 k htr pjn : k h pjn þ c n u h pjn g pn ; ð39aÞ with some fixed constants c n , c t > 0, as well as the combined friction bound g pjt : g pt þ F max 0; k htr pjn : Then, we can rewrite the contact conditions (37) and (38) The equivalence of (37), (38) and (40) Thus, the set A h n contains all nodes which are in contact with the obstacle, whereas A h t and I h t represent the slippy 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 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 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 kth iteration of the semismooth Newton loop, we compute the generalized derivative of the NCP function (42)  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]  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 H c introduced above, we need to compute the integrals 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.

Numerical tests for the nonlinear setting
In this section, we show some numerical results with two dif ferent three dimensional geometries.

Contact approximation with conforming geometry
First, we investigate the influence of the coarse grid approxima tions A H n ðs n Þ; I H t ðs t Þ defined in (44) and (45) on the convergence rate. For this, we consider the simple three dimensional domain X = (0, 1) Â (0, 1) Â (0, 2) with the patch x = (0, 1) 3 . We impose a fixed displacement of (0, 0, 0.03) T on the top, unilateral contact with the obstacle . We consider frictionless as well as frictional contact, the latter case with g t = 0 and F 2 f0; 0:02; 0:2; 0:35; 0:5; 1; 5g. In Fig. 6, the fine active set A h n of the converged solution for F 0 is depicted, as well as the coarse grid approximations A H n ðs n Þ for s n 2 {0, 0.4}. Table 2  The remaining rows of Table 2 refer to the frictional case where we approximate the normal active set by (44) with s n = 0 and com pare the error reduction rates for different values of s t . The nota tion 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 s t = 0 whose er ror reduction is rather independent of the value of the friction coef ficient. In contrast, smaller coarse sets I H t obtained by larger values of s t yield a more unstable convergence behavior, and the extreme  ; obtained for s t = 2 hardly converges at all. This suggests that for a good convergence rate, one should choose either s 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 s n , s 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 . = 10 2 , Dt = 10 3 , we have obtained a convergence rate of 0.29 independent of the threshold value s n .

Tire application
Next, we apply our algorithm to a complex 3D geometry con sisting 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 X as well as the multiply connected fine patch x are sketched, with a differ ence in the mesh sizes of about 4 6 H/h 6 8. In contrast to all pre vious numerical results, the geometry is nonconforming because the fine triangulation T h features some details at the potential con tact boundary c 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 pre scribed, whereas all other boundaries are free. Further, a volume force of (0, 0, 1000) is applied to the tire, and the material param eters read E N = E x = 2.5 10 7 and m = 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 h n of the corresponding mortar solution is sketched on the left of Fig. 8, whereas its coarse grid approximation A H n ð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 a = 0.6, l max = 4 inner iterations, and different threshold values s 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 s n is too large, the coarse contact set becomes smaller which de grades or even disables the convergence, as already observed in Section 6.1.
From now on, we fix s n = 0.15 and investigate the performance of Algorithm 2 for different damping parameters a. On the left side of Fig. 10, the error decay for l max = 4 and a 2 {0.5, 0.6, 0.7, 0.8} is shown. For larger values of a, 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 fac tor is better for larger values of a, varying between 0.49 for a = 0.7 and 0.66 for a = 0.5. This suggests an adaptive choice of a with re spect to a change in the active sets.
In order to decrease the number of total iterations, we replace the condition l 6 l max by the following stopping criterion in the spirit of (35): The corresponding error decay is depicted on the right of Fig. 10.
Comparing these results with those in the left picture, one can ob serve that the number of total iterations until convergence has been reduced by around 10, and that also the computation with a = 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 N = E x = 2.5 10 7 , m = 0.33, . = 10 5 and c m = 0.5. We apply no volume forces and enforce frictional contact with the coefficient 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 c as described in [29]. We compute 10 time steps with a step size of Dt = 10 5 and prescribe a time dependent displacement of t 10Dt ð0; 0; 0:75 þ signðx 1 Þ dispÞ; disp 2 f0; 2g, on the top, whereas all other boundaries are free. The effective stress r eff of the corresponding solutions at time t 10 is depicted in Fig. 11. The problem is solved using Algorithm 2 with a = 0.6, s n = s t = 0.15 and l 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 jA h n j and jI h t j 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.    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.

Local time subcycling
The domain decomposition approach introduced in the previ ous 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 asymp totical convergence order for Dt ? 0. Afterwards, we sketch in Sec tion 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 h x is con structed in such a way that its entries associated with the interface C vanish. Such a modification has already been successfully ap plied 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 ma trix is presented in [57].

Interface coupling conditions
From now on, we denote the time step associated with the coarse space V H by DT and the time step on the patch x by Dt.
We assume that the time step sizes have an integer ratio, i.e., The feedback to the coarse grid Eq. (14) 1 is carried out by the mean value of the Lagrange multipliers, leading to The Schur matrices S H N ; S h x used above have been defined in (13), and the residuals q H N0 ; q h xðj 1Þ 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 conserv ing, it permits the a priori error estimates stated in the next theorem.
The constants C of the leading order terms can depend on h, the mate rial parameters and the time derivatives of (50) but are independent of J, j and DT. 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 expan sion. For the details, we refer to [57]. h Theorem 4 implies that the local errors at time t J , i.e., after one coarse time step, are of order OðDT 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 x.
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 N and J fine subproblems on x. 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.

Approximate solution scheme
The idea of Algorithm 1 can be transferred to the coupled sys tem (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 DT, followed by J intermediate fine scale problems with step size Dt and Dirichlet boundary conditions on C. The backcoupling for the next iteration is performed via the difference in the coarse and fine boundary stress at C. The corresponding numerical scheme for the iterative computation of given all values at time t 0 , is stated in Algorithm 3. The residuals of (51), (52)  In this section, we illustrate the performance of Algorithm 3 ap plied to two different two dimensional geometries.

Conforming geometry
First, we investigate the energy conservation and the time discretization error of the coupled system (48) 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 10 4 , the energy is constant with respect to time, in agreement with the the oretical result of Lemma 3.
Next, we fix DT = 4 10 5 and investigate the evolution of the relative time discretization error on the patch x given by e 2 x;time : where the ''reference solution''û has been computed with an over all fine time step of size 2.5 10 7 . The corresponding results for dif ferent 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].

Nonconforming geometry
Finally, we apply Algorithm 3 to the two dimensional geometry sketched in Figs. 1 and 15. The domain X consists of a circular ring centered at the origin with the diameters r inner = 1.6, r outer = 1.95 and 60 additional salients of height 0.05. The fine domain x 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 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 c, 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 H x å V h x . However, the triangulation with H/h = 4 is con structed such that the trace spaces W H C & W h C are nested. We set Dt = 2 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 N = E x = 4.4 10 6 and m = 0.33. The density on the inner ring is given by . N = 1, whereas the sali ents have a density of 10 3 . The structure has an initial velocity of v 0 ðxÞ 150 x 2 x 1 0 20 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 c 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 kR (k) k 6 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 DT = 2 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.

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 abbrevi ation k k k;H : k k H k ðHÞ and denote the space RB x of rigid body modes on x as well as the corresponding seminorm by Further, according to (11a), we define the bilinear form k x ðv; wÞ : 2 Dt 2 m x ðv; wÞ þ a x ðv; wÞ; v; w 2 V x : Due to the fact that the Schur complement matrices S H x ; S HhH x de fined in (16) and (21) correspond to discrete harmonic extensions onto x with respect to k x ( , ), we obtain w H C ; S HhH The idea of the proof is to project the discrete functionv h 2 V h x onto the space V H x and vice versa, and to construct an upper bound of the terms in (A.1) from these projections. A suitable projection operator Z m x : V x ! V m x ; m 2 fh; Hg, can be constructed similar to the Scott Zhang operator Z h defined in [59]. The only difference is that we change the original definition of [59] near the interface C such that the relation ðZ m x vÞj C P m ðvj C Þ holds for v 2 V x , where the mortar projection P m has been defined in (5). With this, the operator Z m x sat isfies the approximation properties 0 6 l 6 k 6 2; k >