- Research article
- Open access
- Published:
Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems
Advanced Modeling and Simulation in Engineering Sciences volume 5, Article number: 31 (2018)
Abstract
The aim of the present paper is to study theoretically and numerically the Verlet scheme for the explicit time-integration of elastodynamic problems with a contact condition approximated by Nitsche’s method. This is a continuation of papers (Chouly et al. ESAIM Math Model Numer Anal 49(2), 481–502, 2015; Chouly et al. ESAIM Math Model Numer Anal 49(2), 503–528, 2015) where some implicit schemes (theta-scheme, Newmark and a new hybrid scheme) were proposed and proved to be well-posed and stable under appropriate conditions. A theoretical study of stability is carried out and then illustrated with both numerical experiments and numerical comparison to other existing discretizations of contact problems.
Introduction and problem setting
Explicit time-marching schemes for the dynamics of deformable solids with impact has already been the subject of an abundant literature (see, e.g., [1,2,3] for some recent contributions). They are appealing since they can be easy to implement, fast and adapted to parallel architectures. Nevertheless, there still remains important difficulties to design robust explicit methods and to obtain reliable numerical simulations in this context (see, e.g., [4]). Among these difficulties, numerical stability and energy conservation remains one of the most important ones. Another one is to preserve the quality and the accuracy of the numerical solution, which can present spurious oscillations in the displacement, the velocity or the contact stress. A last one is to enforce properly the contact condition, particularly the non-penetration condition.
A precursory method is the one developed by Taylor and Flanagan [5] in the framework of PRONTO3D software (see also the description in [6]). Nevertheless, the method is not fully explicit, except in a node-to-node contact approximation, in the sense that the contact pressure is computed in an iterative process on the whole contact surface. To mention some other of the most important contributions, we can say that a widely resumed theoretical work in dynamic impact problems is due to Moreau [7, 8] for the impact of rigid body systems. The (implicit) schemes proposed by Moreau have been extended quite naturally to the elasticity case through finite element semi-discretization in space (for instance in [9]) which transforms the continuous impact problem into a discrete one very close to a rigid body system. These discrete impact problems, governed by a so-called measure differential inclusion are notoriously ill-posed and of very low regularity.
The ill-posedness can be fixed (for the most part) by the addition of an impact law with a restitution coefficient. As a matter of fact, standard schemes, such as the commonly used ones of Newmark’s family [10], have an erratic behavior when they are applied to dynamic contact problems. This is mainly because they select a solution corresponding to an arbitrary (and potentially very large) restitution coefficient (see [11]). Alternatively, a valuable scheme in this context is that of Paoli and Schatzman [12, 13] who implicitly takes into account this restitution coefficient. However, the addition of a restitution coefficient can be considered as artificial in the context of deformable solids. This does not diminish the interest for the Paoli–Schatzman scheme which will be a point of comparison with our proposed approach. The implicit inclusion of a restitution coefficient has also been considered in [14] to develop a wide range of schemes based on a time discontinuous Galerkin framework.
As noticed in [11], even in the case where the continuous problem is well-posed (see, e.g., [15, 16] for well-posedness results), the ill-posed measure differential inclusion that results from finite element semi-discretization in space has an infinite number of solutions, depending on the choice of a restitution coefficient on each node of the contact boundary. Moreover, it is not possible to decide which solution is more suitable than other. Indeed, the two most remarkable solutions are, on the one hand, the one for a unitary restitution coefficient which ensures conserving energy but which causes very important spurious oscillations of the contact nodes and unexploitable contact stress, and, on the other hand, the solution for a vanishing restitution coefficient which ensures stability and a better approximation of the contact stress but is energy dissipative, while the continuous problem is not. This resulted in [11] to design the mass redistribution method (generalized in [17, 18]) which allows a compromise in this context, i.e. well-posedness of the space semi-discretized problem, conservation of the energy and an improved quality of the contact stress. However, and this is also the case for the Paoli–Schatzman scheme, it introduces a global problem to be solved (at least on the contact nodes) when an explicit time-marching scheme is used. In the same spirit, a time-marching scheme has been designed in [19] for dynamic fracture problems, in which the cohesive forces are treated implicitly, while an explicit scheme is used for the dynamics of interior nodes.
For explicit time-integration, primal formulations of contact conditions are better suited. Indeed, since no additional unknown such as a Lagrange multiplier are introduced, they allow to enforce the contact conditions at the previous time-step, instead of the current one, so that the contact term appears at the right-hand side and does not require global (and non-linear) solving. A first possibility is to penalize / regularize the contact conditions (see, e.g., [20, 21]): the resulting penalty method is simple to implement and only an inversion of the mass-matrix is needed at each time-step to solve the resulting fully discretized problem (and the scheme becomes fully explicit when the mass matrix is lumped). Nevertheless, the penalty method is not consistent and the choice of the penalty parameter remains a difficulty (see, e.g., [22]). The alternative we explore in this paper is a Nitsche treatment of contact conditions, which is still a primal method, with the same advantages as penalty, but that remains consistent with the original problem, and more robust with respect to the Nitsche parameter. Nitsche’s method, originally designed to enforce weakly Dirichlet boundary conditions [23, 24], was adapted to unilateral contact in [25, 26] (see also [27] for an overview of recent results on this topic).
We studied previously in [28, 29] the behavior of Nitsche’s method for contact in elastodynamics, when combined to various implicit time-marching schemes. Particularly, when applied to contact-impact in elastodynamics, Nitsche’s method has the good property of leading to a well-posed semi-discrete problem in time (i.e., a system of Lipschitz differential equations) as it is shown in [28]. This feature is shared also by the penalty method and modified mass methods. Moreover the symmetric variant of Nitsche’s space semi-discretization conserves an augmented energy [28], as does the penalty method [30]. We studied as well theoretically the well-posedness, the stability and energy conservation properties of fully discrete schemes based on space semi-discretization with Nitsche’s method combined with the theta-scheme, the Newmark scheme and a new Hybrid scheme. This study was illustrated with some numerical experiments.
The aim of this paper is to study mathematically and numerically the approximation of contact problems in elastodynamics by Nitsche’s method combined with the explicit Verlet time-marching scheme. The choice of the Verlet scheme is motivated both by its simplicity and its attractive theoretical properties (symplecticity) [31]. We will also make comparisons with some of the existing methods mentioned above and with the approximation by penalized contact. The numerical comparison will be mainly performed on the one-dimensional problem introduced in [15] whose advantage is to present a known periodic solution and to make clear the occurrence of parasitic oscillations, the convergence and energy conservation properties. Comparisons for 2D and 3D problems will also be presented.
Let us introduce some useful notations. In what follows, bold letters like \(\mathbf{u },\mathbf{v }\), indicate vector or tensor valued quantities, while the capital ones (e.g., \(\mathbf{V },\mathbf{K }\ldots \)) represent functional sets involving vector fields. As usual, we denote by \((H^{s}(.))^d\), \(s\in \mathbb {R}, d=1,2,3\) the Sobolev spaces in one, two or three space dimensions (see [32]). The usual scalar product of \((H^{s}(D))^d\) is denoted by \((\cdot ,\cdot )_{s,D}\) and the corresponding norm is denoted by \(\Vert \cdot \Vert _{s,D}\)—we keep the same notation when \(d=1\) or \(d>1\). The letter C stands for a generic constant, independent of the discretization parameters.
We consider an elastic body \(\Omega \) in \(\mathbb {R}^d\) with \(d=1,2,3\). Small strain assumptions are made (as well as plane strain when \(d=2\)). The boundary \(\partial \Omega \) of \(\Omega \) is polygonal (\(d=2\)) or polyhedral (\(d=3\)). The normal unit outward vector on \(\partial \Omega \) is denoted \(\mathbf{n }\). We suppose that \(\partial \Omega \) consists in three nonoverlapping parts \(\Gamma _D\), \(\Gamma _N\) and the contact boundary \(\Gamma _C\), with meas\((\Gamma _D) > 0\) and meas\((\Gamma _C) > 0\). In its initial stage, the body is in contact on \(\Gamma _C\) with a rigid foundation and we suppose that the unknown contact zone during deformation is included into \(\Gamma _C\). The body is clamped on \(\Gamma _D\) for the sake of simplicity. It is subjected to volume forces \(\mathbf{f }\) in \(\Omega \) and to surface loads \(\mathbf{g }\) on \(\Gamma _N\).
We deal with the unilateral contact problem in linear elastodynamics during a period of time [0, T] where \(T>0\) is the final time. We denote by \(\Omega _{T}: = (0,T] \times \Omega \) the time-space domain, and similarly \(\Gamma _{DT}:= (0,T] \times \Gamma _D\), \(\Gamma _{NT}:= (0,T] \times \Gamma _N\) and \(\Gamma _{CT}:= (0,T] \times \Gamma _C\). The problem then consists in finding the displacement field \(\mathbf{u }: [0,T] \times \Omega \rightarrow \mathbb {R}^d\) verifying the equations and conditions (1, 2):
where the notation \({\dot{\mathbf{x }}}\) is used for the time-derivative of a vector field \(\mathbf{x }\) on \(\Omega _{T}\), so that \({\dot{\mathbf{u }}}\) is the velocity of the elastic body and \(\ddot{\mathbf{u }}\) its acceleration; \(\mathbf{u }_0\) and \({\dot{\mathbf{u }}}_0\) are initial displacement and velocity. The density of the elastic material is denoted by \(\rho \), and is supposed to be a constant to simplify the presentation (this is not restrictive and the results can be extended straightforwardly for a variable density). The notation \({\varvec{\sigma }}= (\sigma _{ij}), \;1\le i,j \le d,\) stands for the stress tensor field and \(\mathbf{div }\) denotes the divergence operator of tensor valued functions. The notation \({\varvec{\varepsilon }}(\mathbf{v }) = ({\varvec{\nabla }}\mathbf{v }+ {\varvec{\nabla }}\mathbf{v }^{^T})/2\) represents the linearized strain tensor field and \( \mathbf{A }\) is the fourth order symmetric elasticity tensor having the usual uniform ellipticity and boundedness property. For any displacement field \(\mathbf{v }\) and for any density of surface forces \({\varvec{\sigma }}(\mathbf{v }) \mathbf{n }\) defined on \(\partial \Omega \) we adopt the following notation
where \(\mathbf{v }_\mathbf{t }\) (resp. \({\varvec{\sigma }}_{\mathbf{t }}(\mathbf{v })\)) is the tangential component of \(\mathbf{v }\) (resp. \({\varvec{\sigma }}(\mathbf{v }) \mathbf{n }\)). The conditions describing unilateral contact without friction on \(\Gamma _{CT}\) are:
Note additionally that the initial displacement \(\mathbf{u }_0\) should satisfy the compatibility Condition \(u_{0n} \le 0\) on \(\Gamma _C\).
To our knowledge, the well-posedness of Problems (1), (2) is still an open issue. The few available existence results concern simplified model problems involving the (scalar) wave equation with Signorini’s conditions (see, e.g., [16, 33,34,35,36]) or thin structures like membranes, beams (see [37]) or plates (see [38]). Even in these simplified cases, obtaining uniqueness and energy conservation still involves difficulties in 2D or 3D. For a review on some of these results, one can refer to the book [39].
We introduce the Hilbert space
and the following forms:
for any \(\mathbf{u }\) and \(\mathbf{v }\) in \(\mathbf{V }\), for all \(t\in [0,T]\). The (total) mechanical energy associated with the solution \(\mathbf{u }\) of the dynamic contact problem (1, 2) is:
Let us take \(t\in [0,T]\). Formally, we get from (1), after multiplication by \({\dot{\mathbf{u }}}(t)\), integration by parts, with the boundary conditions on \(\Gamma _{DT}\), \(\Gamma _{NT}\) and the absence of friction:
Moreover, if we assume that the contact force does not dissipate any energy, i.e. satisfies the so called persistency Condition \(\sigma _n(\mathbf{u }(t)) \dot{u}_n (t) = 0\) (see, e.g., [30, 40, 41]) then we end up with:
Notably, when L vanishes, we get energy conservation: \(E(t) = E(0)\), for all \(t\in [0,T]\).
Note that, even if it is expected that solutions to Problems (1), (2) satisfy the persistency Condition \(\sigma _n(\mathbf{u }(t)) \dot{u}_n (t) = 0\) in order to respect the non-dissipative character of the frictionless contact condition, it has only been rigorously proved in a one dimensional framework (elastic bar) for instance in [36, Lemma 2.5].
The rest of our paper is outlined as follows. The first section is dedicated to the description of the fully discrete formulation for dynamic contact with Nitsche and Verlet explicit time-integration. Then, a stability analysis is carried out, and finally, some numerical comparisons with other classical methods are investigated and analysed.
Discrete setting: Nitsche’s method with Verlet scheme
We begin this section with preliminary notations and results. Then, we introduce our Nitsche-based finite element semi-discretization in space, and we recall its main properties of well-posedness and energy conservation. Finally we describe the fully discretized problem based on the Verlet explicit time-marching scheme.
Preliminary notations and results
We make use of the notation \(\left[ \cdot \right] _{_{\mathbb {R}^-}}\), that stands for the projection onto \(\mathbb {R}^-\) (\(\left[ x \right] _{_{\mathbb {R}^-}} = \frac{1}{2} ( x - |x| )\) for \(x\in \mathbb {R}\)). The notation \(H(\cdot )\) will stand for the Heaviside function \(H(x) = 1 \hbox { if } x > 0, \frac{1}{2} \hbox { if } x = 0, \hbox { and } 0 \hbox { if } x < 0\), which satisfies \(H(x) + H(-x)=1, \forall x \in \mathbb {R}\). Moreover we will make use of the equality \( H(-x) [ x ]_{_{\mathbb {R}^-}} = [ x ]_{_{\mathbb {R}^-}}, \quad \forall x\in \mathbb {R}\), and the following property of projection:
Let \(\mathbf{V }^h \subset \mathbf{V }\) be a family of finite dimensional vector spaces (see [42]) indexed by h coming from a family \(\mathcal{T}^h\) of triangulations of the domain \(\Omega \) (\(h = \max _{K \in \mathcal{T}^h} h_K\) where \(h_K\) is the diameter of the triangle K). The family of triangulations is supposed:
Regular, i.e., there exists \(\sigma >0\) such that \(\forall K \in \mathcal{T}^h, h_K / \rho _K \le \sigma \) where \(\rho _K\) denotes the radius of the inscribed ball in K,
Conformal to the subdivision of the boundary into \(\Gamma _D\), \(\Gamma _N\) and \(\Gamma _C\), which means that a face of an element \(K \in \mathcal{T}^h\) is not allowed to have simultaneous non-empty intersection with more than one part of the subdivision,
Quasi-uniform, i.e., there exists \(c>0\), such that, \(\forall h>0,~\forall K\in \mathcal{T}^h, ~h_K \ge c h\).
To fix ideas, we choose a standard Lagrange finite element method of degree k with \(k=1\) or \(k=2\), i.e.:
However, our results would be similar for any \(\mathscr {C}^0\)-conforming finite element method.
We consider in what follows \({\gamma _h}\), a positive piecewise constant function on the contact interface \(\Gamma _C\) which satisfies for every K that has a non-empty intersection of dimension \(d-1\) with \(\Gamma _C\)
where \(\gamma _{0}\) is a positive given constant (the Nitsche parameter). Note that the value of \({\gamma _h}\) on element intersections has no influence.
We next define convenient mesh-dependent norms, in fact weighted \(L^2(\Gamma _C)\)-norm (since \((\gamma _0 / {\gamma _h})|_K = h_K\)).
Definition 1
For any \(v \in L^2(\Gamma _C)\), we set
Additionally, it will be convenient to endow \(\mathbf{V }^h\) with the following mesh- and parameter-dependent scalar product:
Definition 2
For all \(\mathbf{v }^h , \mathbf{w }^h \in \mathbf{V }^h\) we set
and note \(\Vert \cdot \Vert _{\gamma _h}:= ( \cdot , \cdot )^{\frac{1}{2}}_{\gamma _h}\) the corresponding norm. Remark that the two norms \(\Vert \cdot \Vert _{\gamma _h}\) and \(\Vert \cdot \Vert _{1,\Omega }\) are equivalent on \(\mathbf{V }^h\), in the following sense (for a quasi-uniform mesh \(\mathcal{T}^h\)):
for any \(\mathbf{v }^h \in \mathbf{V }^h\). The positive constant C comes from the trace inequality and the constant of quasi-uniformity of the mesh \(\mathcal{T}^h\). For a mesh \(\mathcal{T}^h\) that is not quasi-uniform, the same relationship holds, replacing h by \((\min _{K \in \mathcal{T}^h} h_K)\).
We end this section with the following statement: a discrete trace inequality (see, e.g., [43]), that is a key ingredient for the whole mathematical analysis of Nitsche’s based methods.
Lemma 3
There exists \(C > 0\), independent of the parameter \(\gamma _0\) and of the mesh size h, such that, for all \(\mathbf{v }^h \in \mathbf{V }^h\)
Semi-discrete problem in space
Our Nitsche-based discretization of the contact condition comes from the following result (see [44] and as well [25] for a detailed formal proof).
Proposition 4
Let \(\gamma \) be a positive function defined on \(\Gamma _C\). The contact Condition (2) can be reformulated as follows:
As in [28, 29] we will consider a family of methods indexed by a parameter \(\Theta \in \mathbb {R}\) (with, in general, \(\Theta =-1,\,0,\,1\), see, e.g., [26]). Let us introduce the discrete linear operator
Define as well the bilinear form:
The space semi-discretized Nitsche-based method for unilateral contact problems in elastodynamics then reads (see, e.g, [27, 28]):
where \(\mathbf{u }^h_0\) (resp. \(\dot{\mathbf{u }}^h_0\)) is an approximation in \(\mathbf{V }^h\) of the initial displacement \(\mathbf{u }_0\) (resp. the initial velocity \(\dot{\mathbf{u }}_0\)), for instance the Lagrange interpolant or the \(L^2(\Omega )\) projection of \(\mathbf{u }_0\) (resp. \(\dot{\mathbf{u }}_0\)).
Remark 5
Note that, as in [27], we adopted in this presentation a different convention for notations compared to previous works [28, 29]. This is in order to get closer to the formulations provided in most of the papers on Nitsche’s method and on the augmented Lagrangian method.
We can reformulate (8) as a system of (non-linear) second-order differential equations. To this purpose, using Riesz’s representation theorem in \((\mathbf{V }^h , (\cdot ,\cdot )_{{\gamma _h}})\) we first introduce the mass operator \(\mathbf{M }^h : \mathbf{V }^h \rightarrow \mathbf{V }^h\), which is defined for all \(\mathbf{v }^h,\mathbf{w }^h \in \mathbf{V }^h\) by \( ( \mathbf{M }^h \mathbf{v }^h , \mathbf{w }^h )_{{\gamma _h}} = ( \rho \mathbf{v }^h , \mathbf{w }^h )_{0,\Omega }. \) Still using Riesz’s representation theorem, we define the (non-linear) operator \(\mathbf{B }^h : \mathbf{V }^h \rightarrow \mathbf{V }^h\), by means of the formula
for all \(\mathbf{v }^h,\mathbf{w }^h \in \mathbf{V }^h\). Finally, we denote by \(\mathbf{L }^h(t)\) the vector in \(\mathbf{V }^h\) such that, for all \(t\in [0,T]\) and for every \(\mathbf{w }^h\) in \(\mathbf{V }^h\): \( (\mathbf{L }^h(t),\mathbf{w }^h)_{{\gamma _h}} = L(t)(\mathbf{w }^h). \) With the above notation, Problem (8) reads:
Moreover, we recall the results of well-posedness and the energy estimate for the semi-discrete problem in space, that were established in [28]. First, the following theorem together with the boundedness of \( \Vert (\mathbf{M }^h)^{-1} \Vert _{{\gamma _h}}\) [see [28] show that Problem (8) or equivalently Problem (9)] is well-posed.
Theorem 6
The operator \(\mathbf{B }^h\) is Lipschitz-continuous in the following sense: there exists a constant \(C > 0\), independent of h, \(\Theta \) and \(\gamma _0\) such that, for all \(\mathbf{v }^h_1, \mathbf{v }^h_2 \in \mathbf{V }^h\):
As a consequence, for every value of \(\Theta \in \mathbb {R}\) and \(\gamma _0 > 0\), Problem (8) admits one unique solution \(\mathbf{u }^h~\in ~{{\mathscr {C}}}^2([0,T],\mathbf{V }^h)\).
Remark 7
Note that, conversely to the static case (see [25, 26, 45]) and the fully-discrete case there is no condition on \(\gamma _0\) for the space (semi-)discretization, which remains well-posed even if \(\gamma _0\) is arbitrarily small. However, this does not imply that the solution remains consistent when \(\gamma _0\) becomes small (see Remark 19 and Fig. 4 in the sequel).
We recall that the standard (mixed) finite element semi-discretization for elastodynamics with unilateral contact leads to ill-posed problems (see, e.g., [11, 22]), which is not the case of Nitsche’s formulation that leads to a well-posed (Lipschitz) system of differential equations. This feature is shared with the standard penalty method, the difference being that Nitsche’s method remains consistent in a strong sense (see [28]). Note that the standard (mixed) finite element semi-discretization is consistent as well as the singular dynamic method introduced in [18]. The mass redistribution method introduced in [11] is asymptotically consistent when h vanishes.
Now, we consider the energy estimates which are counterparts of the Eq. (3), in the semi-discretized case. Let us define the discrete energy as follows:
which is associated to the solution \(\mathbf{u }^h(t)\) to Problem (8). This is the direct transposition of the mechanical energy E(t) from the continuous system. As in [28], we define also a modified energy more suited to Nitsche’s method
in which a consistent term is added. This term denoted \(R^h(t)\) represents, roughly speaking, the non-fulfillment of the contact Condition (7) by \(\mathbf{u }^h\). The relationship between \(E^h(t)\) and \(E^h_{\Theta }(t)\) is provided in the Lemma below:
Proposition 8
For \(\Theta \ge 0\) and \(\gamma _0\) large enough, there exists \(C > 0\) independent of h, of \(\gamma _0\) and of the solution to Problem (8), such that, for all \(t \in [0,T]\):
Proof
This result is obtained using the coercivity of \(a(\cdot ,\cdot )\) and applying Lemma 3. \(\square \)
Remark 9
Proposition 8 states that the energy \(E^h_{\Theta }(t)\) remains always positive (if \(E^h(0)\) is) for \(\Theta \ge 0\) and \(\gamma _0\) large enough. For \(\gamma _0\) small, the existence of zero energy spurious modes cannot be excluded.
Remark 10
For \(\Theta < 0\), such a result with a constant independent of the mesh parameter h cannot be obtained. As a consequence, for \(\Theta < 0\), the quantity \(E^{h}_{\Theta } (t)\) cannot be used for optimal energy evolution estimates and might become even negative for h small.
Still in [28], the following evolution of \(E^h_{\Theta }\) is obtained:
Theorem 11
Suppose that the system associated to (1, 2) is conservative, i.e., that \(L(t) \equiv 0\) for all \(t\in [0,T]\). The solution \(\mathbf{u }^h\) to (8) then satisfies the following identity:
Notably, when \(\Theta = 1\), we get for any \(t\in [0,T]\): \( E^h_{1} (t) = E^h_{1} (0). \)
This result links the non-satisfaction of the energy conservation to the non-satisfaction of the so-called persistency condition. However, it appears in the present study that it would be preferable to use \(E^h_{1} (t)\) even for the variants \(\Theta \ne 1\) (see Remark 10), for which the following result can be established:
Theorem 12
Suppose that the system associated to (1, 2) is conservative, i.e., that \(L(t) \equiv 0\) for all \(t\in [0,T]\). The solution \(\mathbf{u }^h\) to (8) then satisfies the following identity:
Proof
Let us take \(\mathbf{v }^h = \dot{\mathbf{u }}^h(t)\) in (8):
which we reformulate as:
We split the second term, use \( \mathrm {P}^{\mathbf{n }}_{\Theta , {\gamma _h}}(\dot{\mathbf{u }}^h (t) ) = \Theta {\sigma _n}(\dot{\mathbf{u }}^h (t)) - {\gamma _h}\dot{u}^h_n = \mathrm {P}^{\mathbf{n }}_{1, {\gamma _h}}(\dot{\mathbf{u }}^h (t)) + (\Theta - 1) \, {\sigma _n}(\dot{\mathbf{u }}^h (t))\) and get:
The result is obtained by re-ordering the terms, using the property \(\frac{d}{dt} \frac{1}{2} [ x(t) ]_{_{\mathbb {R}^-}}^2 = [ x(t) ]_{_{\mathbb {R}^-}} \dot{x}(t)\) and the definition of \(E^h_1(t)\). \(\square \)
Remark 13
The above result still states that \(E_1^h(t)\) is conserved for the symmetric variant \(\Theta =1\), and, for \(\Theta \ne 1\) the variations of \(E_1^h(t)\) come from the non-fulfillment of the contact Condition (7) by \(\mathbf{u }^h\).
Verlet scheme
Let \(\tau > 0\) be the time-step, and consider a uniform discretization of the time interval [0, T]: \((t^0,\ldots ,t^N)\), with \(t^n = n\tau \), \(n=0,\ldots ,N\). Let \(\theta \in [0,1]\), we use the notation:
for arbitrary quantities \(\mathbf{x }^{h,n}, \mathbf{x }^{h,n+1} \in \mathbf{V }^h\). Hereafter we denote by \(\mathbf{u }^{h,n}\) (resp. \(\dot{\mathbf{u }}^{h,n}\) and \(\ddot{\mathbf{u }}^{h,n}\)) the resulting discretized displacement (resp. velocity and acceleration) at time-step \(t^{n}\).
The discretization of Problem (9) with the velocity–Verlet scheme reads:
with initial Conditions \(\mathbf{u }^{h,0} = \mathbf{u }^h_0\) , \(\dot{\mathbf{u }}^{h,0} = \dot{\mathbf{u }}^h_0\) and \(\ddot{\mathbf{u }}^{h,0} = \ddot{\mathbf{u }}^h_0\), and the notation \(\mathbf{L }^{h,n+1} = \mathbf{L }^h (t^{n+1})\), the initial value \(\ddot{\mathbf{u }}^h_0\) being obtained through \( \mathbf{M }^h\ddot{\mathbf{u }}^{h}_0 = \mathbf{L }^{h,0} - \mathbf{B }^h\mathbf{u }^{h,0}\).
This scheme corresponds to the variant of the Newmark scheme with \(\gamma = \frac{1}{2}\) and \(\beta =0\) (see, e.g., [22, 28]). As a result, this is a second order consistent scheme in \(\tau \). Note that, for practical implementation, the acceleration can be eliminated using the relationship \( \mathbf{M }^h\ddot{\mathbf{u }}^{h,n} = \mathbf{L }^{h,n} - \mathbf{B }^h\mathbf{u }^{h,n}\). This result into the following reformulation, where the only unknowns are the displacement and the velocity:
Since this scheme only requires the inversion of the mass matrix \(\mathbf{M }^h\) at each time-step, it becomes then fully explicit when the mass matrix \(\mathbf{M }^h\) is lumped.
We can transform the Scheme (13) into a two-steps scheme. Indeed, the first equation of (13) applied for time-steps n and \(n+1\) reads:
So, using the above relationships, the second equation of (13), at time-step n, can be written as
which can be simplified as:
We add the above equation to the second relationship in (15) and get
Using finally the third equation in (13), combined with the above relationship, we obtain:
This is a two-steps scheme, so called Störmer–Verlet scheme or central difference scheme, that involves only the displacement as an unknown (the first step \(\mathbf{u }^{h,1}\) is classically obtained via the first equation of (14) at \(n=0\)). Note finally that Leapfrog scheme is also equivalent to Verlet one (see e.g. [31, 46]): it suffices to define an intermediate velocity \(\dot{\mathbf{u }}^{h,\left[ n+\frac{1}{2}\right] }\) at half-time-steps \(t^{n+\frac{1}{2}}\) as follows:
where we used the notation \(\dot{\mathbf{u }}^{h,\left[ n+\frac{1}{2}\right] }\) so as to differentiate this new unknown from \(\dot{\mathbf{u }}^{h,n+\frac{1}{2}}\) defined earlier. Using this new intermediate velocity, the Scheme (13) is reformulated as:
with initial Conditions \(\mathbf{u }^{h,0} = \mathbf{u }^h_0\) , \(\dot{\mathbf{u }}^{h,\left[ 1/2\right] } = \dot{\mathbf{u }}^h_0 + \frac{\tau }{2}\ddot{\mathbf{u }}^h_0\).
Stability properties of Verlet scheme
First, we present different energies associated to the solution to Problem (13), and make explicit their relationships. Then, we derive energy estimates associated to the fully discrete Problem (13), and a (non-optimal) stability result is deduced. We conclude with some comments and arguments that a better result may be expected.
Discrete energies
We next define the following energy:
which is associated with the solution \(\mathbf{u }^{h,n}\) to Problem (13). Set also
Note that the energies \(E^{h,n}\) and \(E^{h,n}_{\Theta }\) are the fully discrete counterparts of the semi-discrete energies \(E^{h}(t)\) and \(E^{h}_{\Theta }(t)\). Note additionally that a result similar to Proposition 8 can be established, that allows to bound the physical energy by the augmented energy under appropriate conditions on the numerical parameters (and the statements of Remark 9 and Remark 10 still applies):
Proposition 14
For \(\Theta \ge 0\) and \(\gamma _0\) large enough, there exists \(C > 0\) independent of h, of \(\gamma _0\) and of the solution to Problem (13), such that, for all \(n=0,\ldots ,N\):
To simplify slightly the notations in the energy estimates below, we will make use of the convention: \(\mathrm {P}^n:= \mathrm {P}^{\mathbf{n }}_{1, {\gamma _h}}(\mathbf{u }^{h,n})\) for any \(n \in \mathbb {N}\). We will denote as well by \(\left[ \cdot \right] _{_{\mathbb {R}^+}}\) the projection onto \(\mathbb {R}^+\). Additionally, we define a modified energy adapted to the Verlet scheme:
Of course, this variant of the energy definition makes also sense and is usable for stability analysis only if it can be used to bound the physical energy \(E^{h,n}\). This is the aim of the following result:
Proposition 15
Suppose that \(L^n \equiv 0\) for all \(n=0,\ldots ,N\) and that the mesh \(\mathcal{T}^h\) is quasi-uniform. Then, for \(\gamma _0\) large enough and \(\frac{\tau }{h}\) small enough, there exists \(C > 0\), independent of h, of \(\gamma _0\) and of the solution to Problem (13), such that, for all \(n=0,\ldots ,N\):
Proof
We suppose \(L^n \equiv 0\) and take \(\mathbf{v }^h = \ddot{\mathbf{u }}^{h,n}\) in Problem (13):
We assume that \(\gamma _0\) is large enough, then use the Cauchy–Schwarz inequality, the Lemma 3 and the inverse inequality [47, Corollary 1.141, Remark 1.143] to bound the first right term:
with \(C > 0\) independent of \(\gamma _0\) and h. We use the Cauchy–Schwarz inequality for the second term:
Using once more Lemma 3 we bound:
still with \(C>0\) independent of h and \(\gamma _0\). We then use the trace inequality [48, Theorem 1.6.6] and the inverse inequality [47, Corollary 1.141, Remark 1.143] and get:
We insert the above bounds into (18), take into account that \(\gamma _0\) is large, apply once again the inverse inequality and get:
As a result, we obtain
which allows to conclude that \(E_1^{h,n} \le C E^{h,\tau ,n}\) for \(\frac{\tau }{h}\) small enough using the coercivity of \(a(\cdot ,\cdot )\). The estimate \(E^{h,n} \le C E^{h,\tau ,n}\) is then deduced from Proposition 14. \(\square \)
Energy evolution estimates
First, the straightforward adaptation of [29, Proposition 3.4], taking \(\gamma =\frac{1}{2}\) and \(\beta =0\) for Verlet scheme gives the following energy identity:
Proposition 16
Suppose that \(L^n \equiv 0\) for all \(n\ge 0\). The following energy identity holds for all \(n\ge 0\):
This result can be easily adapted as follows when the energy \(E^{h,n}_{1}\) is considered instead, even for \(\Theta \ne 1\):
Proposition 17
Suppose that \(L^n \equiv 0\) for all \(n\ge 0\). The following energy identity holds for all \(n\ge 0\):
Proof
This identity is deduced from
combined with the following rearrangement (where we use the identity \(\left[ \mathrm {P}^n \right] _{_{\mathbb {R}^-}} = {\sigma _n}(\mathbf{u }^{h,n}) - {\gamma _h}u_n^{h,n} - \left[ \mathrm {P}^n \right] _{_{\mathbb {R}^+}}\)):
and gathering the terms with the ones in the Expression (20) of Proposition 16.\(\square \)
As an interesting consequence, we obtain the following result for the discrete energy \(E^{h,\tau , n}\) by simplifying the previous one:
Proposition 18
Suppose that \(L^n \equiv 0\) for all \(n\ge 0\). The following energy identity holds for all \(n\ge 0\):
Proof
We use Eq. (13) to rewrite:
Then we just use the above identity in (21). \(\square \)
Remark 19
For \(\gamma _0\) small, the property of Proposition 15 can be lost and the energy \(E^{h,n}_{1}\) may become negative. In that case, some deformation corresponding to a negative energy may exist, which is of course a non-physical situation. This highlights that, even thought the semi-discrete problem (9) has a unique solution for \(\gamma _0\) small, the reliability of the discretization is guaranteed only for \(\gamma _0\) large enough.
Remark 20
Still referring to [29, Proposition 3.4], and instead of Verlet scheme, if we consider the explicit Newmark scheme \(\gamma =1\) and \(\beta =0\) and \(\Theta = 1\) as Nitsche’s variant, the pending energy evolution corresponding to Proposition 18 in that case involves the sole term \(\left[ \mathrm {P}^{n+1} \right] _{_{\mathbb {R}^-}} \left[ \mathrm {P}^n \right] _{_{\mathbb {R}^+}}\) (instead of \(\left[ \mathrm {P}^{n+1} \right] _{_{\mathbb {R}^-}} \left[ \mathrm {P}^n \right] _{_{\mathbb {R}^+}}- \left[ \mathrm {P}^n \right] _{_{\mathbb {R}^-}} \left[ \mathrm {P}^{n+1} \right] _{_{\mathbb {R}^+}}\) for Verlet scheme). This term being non-positive, the stability of this explicit scheme can be established thanks to Proposition 15 for \(\frac{\tau }{h}\) small enough.
Stability analysis in the case \(\Theta = 1\)
The main result of this section is the following (non-optimal) stability result for the Scheme (13) in the case \(\Theta =1\):
Proposition 21
Suppose that \(L^n \equiv 0\) for all \(n\ge 0\) and that the mesh \(\mathcal{T}^h\) is quasi-uniform. Then, for \(\Theta =1\), for \(\gamma _0\) large enough and for
with \(C > 0\) independent of \(\gamma _0\), h and \(\tau \), the energies \(E^{h,\tau ,n}\) and \(E^{h,n}\) remain bounded.
Proof
We already know from Proposition 18 that, for \(\Theta =1\):
We note that \(\left[ \mathrm {P}^{n+1} \right] _{_{\mathbb {R}^-}} \left[ \mathrm {P}^n \right] _{_{\mathbb {R}^+}} \le 0\) and use \(- \left[ \mathrm {P}^n \right] _{_{\mathbb {R}^-}} \left[ \mathrm {P}^{n+1} \right] _{_{\mathbb {R}^+}} \le \frac{1}{4} ( \left[ \mathrm {P}^{n+1} \right] _{_{\mathbb {R}^+}} - \left[ \mathrm {P}^n \right] _{_{\mathbb {R}^-}} )^2 \le \frac{1}{4} (\mathrm {P}^{n+1}- \mathrm {P}^n)^2\) to get:
Since the mesh is quasi-uniform, we then use Lemma 3, the trace inequality of [48, Theorem 1.6.6] and the inverse inequality [47, Corollary 1.141, Remark 1.143] to obtain, for \(\gamma _0\) large enough
Following exactly the same path as above, but using \(\Vert \cdot \Vert _{0,\Omega } \le \Vert \cdot \Vert _{1,\Omega }\) after the trace inequality, we bound also:
We combine the above result with the bound (19). This yields:
This results into:
from which we can deduce, using Proposition 15:
This means that, still with \(N = \frac{T}{\tau }\),
which remains bounded under the assumption (23). The boundedness of \(E^{h,N}\) is then deduced from Proposition 15. \(\square \)
Comments on the stability analysis
The stability result given by Proposition 21 is submitted to a CFL Condition \(\tau = \mathcal {O} (h^2)\). Of course, this is not the expected one which would be \(\tau = \mathcal {O}(h)\) in accordance with the result of Proposition 15 and the stability analysis of Verlet scheme for a linear evolution equation (see, e.g., [46]). The reason of the non-optimality of Proposition 21 is that we did not succeed to optimally bound the involved contact term \(\left( \left[ \mathrm {P}^{n+1} \right] _{_{\mathbb {R}^-}} \left[ \mathrm {P}^n \right] _{_{\mathbb {R}^+}} -\left[ \mathrm {P}^n \right] _{_{\mathbb {R}^-}} \left[ \mathrm {P}^{n+1} \right] _{_{\mathbb {R}^+}} \right) \). However, an important remark is that this term vanishes unless the terms \(\mathrm {P}^n\) and \(\mathrm {P}^{n+1}\) are of opposite signs, which occurs only when the contact status changes. Moreover it is positive only when the status changes from contact to non-contact. If we assume that the number of such transitions is finite during the simulation, a stability result with a Condition \(\tau = \mathcal {O}(h)\) may be recovered. However, an infinite number of changes of the contact status cannot be excluded. Another argument in favor of such a stability result is obtained via the definition of the following linear (but depending on \(\mathbf{u }^{h,n}\)) operator \(\mathbf{B }^{h,n} : \mathbf{V }^h \rightarrow \mathbf{V }^h\):
for all \(\mathbf{v }^h,\mathbf{w }^h \in \mathbf{V }^h\). Due to the relationship \(\left[ x \right] _{_{\mathbb {R}^-}} = H(-x)\, x\) for \(x\in \mathbb {R}\), there holds:
for all \(\mathbf{w }^h \in \mathbf{V }^h\), and therefore
The following bound on the norm of \(\mathbf{B }^{h,n}\) can be established, following an argument similar to [28, Theorem 2.8]:
Lemma 22
Let us suppose that \(\gamma _0\) is large enough. There exists a constant \(C > 0\), independent of \(\Theta \), \(\gamma _0\) and h, such that
where \(\Vert \cdot \Vert _{{\gamma _h}}\) is the operator norm induced by the norm \(\Vert \cdot \Vert _{{\gamma _h}}\) on \(\mathbf{V }^h\).
Proof
Let us take \(\mathbf{v }^h\) and \(\mathbf{w }^h\) in \(\mathbf{V }^h\). First, using Lemma 3 there holds
and the same bound holds for \(\Vert \mathrm {P}^{\mathbf{n }}_{1, {\gamma _h}}(\mathbf{v }^h) \Vert _{0,\Gamma _C}\), replacing \(| \Theta |\) by 1. Then, using Lemma 3 and the above result, we bound:
which proves the assertion (24). \(\square \)
Using \(\mathbf{B }^{h,n}\) we can rewrite the operator associated to the Scheme (14) as:
where \(\mathbf{I }\) is the identity operator. It links the successive values \(\mathbf{u }^{h,n+1}, \mathbf{u }^{h,n}\) and \(\mathbf{u }^{h,n-1}\) thanks to
when neglecting the source term \(\mathbf{L }^{h,n}\). For a linear problem this would correspond to the amplification matrix. The following result can be established for \(\mathbf{C }^{h,n}\):
Proposition 23
Let us suppose that the mesh \(\mathcal {T}^h\) is quasi-uniform, that \(\gamma _0\) is large and
where \(C > 0\) is a constant independent of \(\gamma _0\), h and \(\tau \). Then, \(\mathbf{C }^{h,n}\) is diagonalizable and its spectral radius \(\rho (\mathbf{C }^{h,n})\) is equal to 1. Furthermore, the same conclusion holds for \(\mathbf{C }^{h,n}\mathbf{C }^{h,n+1}\), which is diagonalizable with a spectral radius \(\rho (\mathbf{C }^{h,n}\mathbf{C }^{h,n+1})\) also equal to 1.
Proof
Let us consider \(\lambda \in \mathbb {C}\) an eigenvalue of the operator \(\mathbf{C }^{h,n}\). Denoting \(\mathbf{A }^{h,n} = 2 \mathbf{I }- \tau ^2(\mathbf{M }^h)^{-1}\mathbf{B }^{h,n}\), this means there exists a non-zero pair \((\mathbf{v }^h,\mathbf{w }^h) \in \mathbf{V }^h \times \mathbf{V }^h\) such that:
With help of the second equation, the first one can be written \(\lambda \mathbf{A }^{h,n} \mathbf{w }^h - \mathbf{w }^h = \lambda ^2 \mathbf{w }^h\), or equivalently
We then use [28, Lemma A.1] and Lemma 22 to bound
We now consider \(\gamma _0^{\frac{1}{2}} \tau / h\) small enough so that the eigenvalues of \(\mathbf{A }^{h,n} = (2 \mathbf{I } - \tau ^2(\mathbf{M }^h)^{-1}\mathbf{B }^{h,n} )\) are positive. We call \(\alpha _j\) these eigenvalues and \(\mathbf{w }^h_j\) the corresponding eigenvectors. Taking \(\mathbf{w }^h = \mathbf{w }^h_j\) in Eq. (25), we deduce that, for each index j, we can compute two eigenvalues \(\lambda ^{\pm }_j\) which are solution to
Since the eigenvalues of \((\mathbf{M }^h)^{-1}\mathbf{B }^{h,n}\) are all positive, we infer \(\alpha _j < 2\), and \(\Delta _j = \alpha _j^2 - 4 < 0\). Therefore the above algebraic equation has two imaginary solutions
Remark that these eigenvalues are such that
This allows to conclude that \(\mathbf{C }^{h,n}\) is diagonalizable and \(\rho (\mathbf{C }^{h,n}) = 1\). We can now make a similar computation for two successive iterations. Since
\(\lambda \in \mathbb {C}\) is an eigenvalue of \(\mathbf{C }^{h,n+1}\mathbf{C }^{h,n}\) if there exists a non-zero pair \((\mathbf{v }^h,\mathbf{w }^h) \in \mathbf{V }^h \times \mathbf{V }^h\) verifying:
which implies
If we denote by \(\beta _j\) the eigenvalues of \(\mathbf{A }^{h,n+1}\mathbf{A }^{h,n}\) which are close to 4 for \(\gamma _0^{\frac{1}{2}} \tau / h\) small enough, and \(\mathbf{v }^h_j\) the corresponding eigenvectors, we conclude that the eigenvalues of \(\mathbf{C }^{h,n+1}\mathbf{C }^{h,n}\) are
Since \(|\lambda _j^{\pm }| = 1\), the operator \(\mathbf{C }^{h,n+1}\mathbf{C }^{h,n}\) is diagonalizable with a unit spectral radius. \(\square \)
Remark 24
For a linear problem, we would conclude that the scheme is stable, under the Condition \(\frac{\tau }{h}\) small enough. However, in a nonlinear framework, the conclusion cannot be drawn since the matrix \(\mathbf{C }^{h,n}\) changes from an iteration to another. Moreover, it seems difficult to pursue the reasoning made on two iterations to an arbitrary number of iterations.
Numerical experiments
We first carry out numerical experiments in 1D, where we can compare our results with an exact solution. Then, numerical experiments in 2D/3D will be described. These numerical tests are performed with the help of our freely available finite element library GetFEM++ (see http://getfem.org).
1D numerical experiments: multiple impacts of an elastic bar
We first present the setting, and then the results obtained by combination of Nitsche’s contact discretization and Verlet scheme. These results are also compared with computations using other methods: the Paoli–Schatzman scheme, the Taylor–Flanagan scheme, the mass redistribution method and the penalty method. This section is ended with numerical convergence tests.
Setting
We first deal with the one-dimensional case \(d=1\) with a single contact point, namely an elastic bar \(\Omega = (0,L)\) with \(\Gamma _{C} = \{ 0\}\), \(\Gamma _{D} = \{ L\}\) and \(\Gamma _{N} = \emptyset \). The elastodynamic equation is then reduced to find \(u : \Omega _{T}= (0,T] \times (0,L) \rightarrow \mathbb {R}\) such that:
where E is the Young modulus and the Cauchy stress tensor is given by \( \sigma (u) = E ({\partial u}/{\partial x})\). Note that \(\sigma _n (u) = (\sigma (u) n) \cdot n = \sigma (u) \) on \(\Gamma _{C}\). In this case, Problem (1, 2) admits one unique solution (see e.g. [36]) for which the following energy conservation equation holds, for t a.e. in (0, T]:
We consider a finite element space using linear (\(k=1\)) or quadratic (\(k=2\)) finite elements and a uniform subdivision of [0, L] with M segments (so \(L=Mh\)). We denote the vector which contains all the nodal values of \(\mathbf{u }^{h,n}\) (resp. \(\dot{\mathbf{u }}^{h,n}\) and \(\ddot{\mathbf{u }}^{h,n}\)) by \(\mathbf{U }^n:= [ U^n_0 , \ldots , U^n_N]^T\) (resp. \(\dot{\mathbf{U }}^n\), \({\ddot{\mathbf{U }}}^n\)). The component of index 0 corresponds to the node at the contact point \(\Gamma _{C}\). We also note \(\mathbf{M }\), resp. \(\mathbf{K }\), the mass, resp. the stiffness, matrix that results from the finite element discretization. We introduce the Courant number which is defined as:
where \(c_0\) is the wave speed associated to the bar. For each simulation, we compute and plot the following time-dependent quantities:
- 1.
The displacement u at the contact point \(\Gamma _C\), given at time \(t^n\) by \(u^{h,n} (0) (= U^n_0)\).
- 2.
The contact pressure \(\sigma _C\), which, in the discrete case, is different from \(\sigma (u)\). If a standard (mixed) method is used for the treatment of contact, it is directly given by the Lagrange multiplier, i.e., \( \sigma _C^n := \lambda ^{h,n}\) at time \(t^n\). In the case of the Nitsche-based formulation, it can be computed as follows at time \(t^n\):
$$\begin{aligned} \sigma _C^n := \left[ \sigma _n (u^{h,n}) (0) - {\gamma _h}(-u^{h,n} (0)) \right] _{_{\mathbb {R}^-}} = \left[ \frac{E}{h} (U_1^n - U_0^n) + \frac{\gamma _0}{h} U^n_0 \right] _{_{\mathbb {R}^-}}, \end{aligned}$$which comes from the contact Condition (7).
- 3.
The discrete energy \(E^h\) which is at time \(t^n\)
$$\begin{aligned} E^{h,n} = \frac{1}{2} \left( (\dot{\mathbf{U }}^n)^{\mathrm {T}} \mathbf{M }{\dot{\mathbf{U }}^n} + (\mathbf{U }^n)^{\mathrm {T}} \mathbf{K }\mathbf{U }^n \right) , \end{aligned}$$and the discrete augmented energy \(E^h_{1}\):
$$\begin{aligned} E^{h,n}_{1} = E^{h,n} - R^{h,n}, ~~R^{h,n} = \frac{h}{2 \gamma _0} \left( (\sigma _n (u^{h,n}) (0))^2 - (\sigma ^n_C)^2 \right) . \end{aligned}$$
We propose a benchmark associated to multiple impacts. This allows to check both the presence of spurious oscillations and the long term energetic behavior of the method. In the absence of external volume forces, the bar is initially compressed. Then, it is released without initial velocity. It impacts first the rigid ground, located at \(x=0\), and then gets compressed again. We take the following values for the parameters: \(f=0\), \(E=1\), \(\rho =1\), \(L=1\), \(u_0 (x)=\frac{1}{2}-\frac{x}{2}\) and \(\dot{u}_0 (x)= 0\). This problem admits a closed-form solution u which derivation and expression is detailed in [15]. Notably, it has a periodic motion of period 3. At each period, the bar stays in contact with the rigid ground during one time unit (see Fig. 1). The chosen simulation time is \(T=12\), so that we can observe 4 successive impacts.
Multiple impacts of an elastic bar. The bar is clamped at \(x=L\) and the contact node is located at the bottom. The closed-form solution is periodic of period 3, with one impact during each period (here between \(t=1\) and \(t=2\))
Numerical results for Nitsche’s method
We discretize the bar with \(M=20\) linear finite elements (\(k=1\), \(h=0.05\)) and take \(\tau =0.01\). The resulting Courant number is \(\nu _C = 0.2\). Note that almost all the parameters have been taken identical as in [29] for comparison purposes. The number of element is smaller (\(M=20\) instead of 100 in [29]) and the time-step \(\tau \) is 0.01 for stability reasons. We first investigate the variant \(\Theta = 0\) with a parameter \(\gamma _0 = 1\). The mass matrix is computed in a standard fashion. The choice \(\gamma _0\) equal to 1 is guided by the concern to obtain a stiffness associated with the degree of freedom on the contact boundary comparable to the stiffnesses obtained by the finite element discretization inside the bar.
The results are depicted in Fig. 2 where the approximated solution corresponds to the solid blue line and the red dotted line is the exact solution. Note also that the dashed energy is \(E^h_1\), the modified energy given by expression (11).
Multiple impacts of an elastic bar. Nitsche’s method with Verlet scheme for \(h=1/20\), \(\tau = 0.01\), \(\Theta = 0\), \(\gamma _0=1\) and \(P_1\) Lagrange finite elements
For an element size which remains relatively rough, one notes the good approximation of the displacement. Significant oscillations can be deplored on the velocity of the contact point, but unfortunately they are very difficult to avoid. Indeed, since a velocity shock is propagating without attenuation or dissipation in the proposed test case, the Gibbs phenomena are inevitable in the finite element approximation. It would be the case even without a contact condition. The approximation of the contact stress is, although polluted by oscillations too, of good quality given the discontinuous character of the exact solution and the relatively coarse mesh size (one can compare with the results obtained in [29] for elements five times smaller). The evolution of the energy reveals some variations which are far from being negligible, but remain moderate, and without appearance of instabilities. Moreover these variations tend to decrease for a finer element size.
Multiple impacts of an elastic bar. Nitsche’s method with Verlet scheme for \(h=1/20\), \(\tau = 0.01\), \(\Theta = -1\), \(\gamma _0=1\) and \(P_1\) Lagrange finite elements
Multiple impacts of an elastic bar. Nitsche’s method with Verlet scheme for \(h=1/20\), \(\tau = 0.01\), \(\Theta = 1\), \(\gamma _0=1\) and \(P_1\) Lagrange finite elements
The calculation for the variant \(\Theta = -1\) and for Nitsche’s parameter \(\gamma _0\) still equal to 1 is presented on Fig. 3. It can be seen that the non-penetration condition is slightly better respected, which indicates that the additional terms compared to the variant \(\Theta = 0\) reinforce the consistency of the method. However, this is at the price of stronger oscillations on the velocity at the contact point. The approximation of the contact stress remains comparable to the \(\Theta = 0\) variant, as well as the energy evolution.
For the same test with the variant \(\Theta = 1\) and Nitsche’s parameter \(\gamma _0 = 1\), as shown in Fig. 4, one gets a non convergent approximation of the displacement. This is due to the loss of coercivity arising when the assumption of Proposition 14 is not satisfied and probably to the existence of zero-energy modes. Indeed, the variant \(\Theta = 1\) is the most restrictive from this point of view and needs a larger parameter \(\gamma _0\). Figure 5 represents the simulation for \(\gamma _0 = 2\) which allows to recover the coercivity. We observe the very good conservation of the discrete energy together with a good approximation of the non-penetration condition. The level of oscillations on the contact point velocity and on the contact stress is similar with the case \(\Theta = 0\).
Multiple impacts of an elastic bar. Nitsche’s method with Verlet scheme for \(h=1/20\), \(\tau = 0.01\), \(\Theta = 1\), \(\gamma _0=2\) and \(P_1\) Lagrange finite elements
Multiple impacts of an elastic bar. Nitsche’s method with Verlet scheme for \(h=1/20\), \(\tau = 0.01\), \(\Theta = 1\), \(\gamma _0=2\), \(P_1\) Lagrange finite elements and a lumped mass matrix
Figure 6 displays a simulation still for (\(\Theta = 1\), \(\gamma _0 = 2\)) but using a lumped mass matrix. Following a standard strategy for \(P_1\) Lagrange finite elements, on each row, the extra-diagonal components of the mass matrix are set to zero and added to the diagonal component. The comparison with Fig. 5 allows to notice more pronounced oscillations on the displacement, the velocity and the contact stress, but still with a very good energy conservation.
Multiple impacts of an elastic bar. Convergence of the displacement for Nitsche’s method for \(P_1\) Lagrange finite elements, \(\Theta = 1\), \(\gamma _0=2\), and for discretization parameters \((h=1/10, \tau =1/50)\), \((h = 1/40, \tau =1/200)\), \((h = 1/160, \tau =1/800)\) and \((h = 1/320, \tau =1/3200)\), respectively
Multiple impacts of an elastic bar. Convergence of the contact stress for Nitsche’s method for \(P_1\) Lagrange finite elements, \(\Theta = 1\), \(\gamma _0=2\), and for discretization parameters \((h=1/10, \tau =1/50)\), \((h = 1/40, \tau =1/200)\), \((h = 1/160, \tau =1/800)\) and \((h = 1/320, \tau =1/3200)\), respectively
Finally, Figs. 7 and 8 show the evolution of the solution for decreasing discretization parameters, and for the variant \(\Theta =1\), \(\gamma _0=2\), and the standard mass matrix. We note in Fig. 7 a rapid decrease of the oscillations on the displacement with the refinement of the discretization. Conversely, the convergence of the contact stress as depicted in Fig. 8 is rather slow, as it could be expected from the very low regularity of the solution. Indeed we observe a very gradual decrease in the amplitude of the oscillations.
Comparison with Paoli–Schatzman scheme
The Paoli–Schatzman scheme directly addresses the measure differential inclusion that results from the finite element semi-discretization of the dynamic contact problem. Following [12, 13, 49,50,51], it can be summarized as follows. Let \(\mathcal{U}_{adm}^h\) be the approximated set of admissible displacements such that \(\mathbf{U }\in \mathcal{U}_{adm}^h\) means that the vector of degrees of freedom \(\mathbf{U }\) satisfies a chosen approximated non-penetration condition. In our one-dimensional test case, this can be simply written \(U_0 \ge 0\), where \(U_0\) is the displacement of the contact point. Then, still denoting \(\mathbf{M }\), resp. \(\mathbf{K }\), the mass, resp. stiffness, matrices that result from the finite element discretization, the (generally ill-posed) measure differential inclusion resulting from the semi-discretization of the dynamic contact Problem (1), (2) can be written:
where \(\partial I_{\mathcal{U}_{adm}^h}(\mathbf{U })\) denotes the normal cone to \(\mathcal{U}_{adm}^h\) at \(\mathbf{U }\). Then, the Paoli–Schatzman scheme can be written for a given restitution coefficient \(e \in [0,1]\):
Of course, one easily recognize the Verlet scheme except for the contact constraint which is taken in an implicit manner and prescribed to the intermediate displacement \(\mathbf{U }^{n+1,e} = \displaystyle \frac{\displaystyle \mathbf{U }^{n+1}+e\mathbf{U }^{n-1}}{\displaystyle 1+e}\).
Multiple impacts of an elastic bar. Paoli–Schatzman scheme for \(h=1/20\), \(\tau = 0.01\), \(e = 0\) and \(P_1\) Lagrange finite elements
When implementation is considered, one usually introduces a multiplier to prescribe the constraint. Denoting by \(\mathbf{B }\) the matrix having \(N_c\) lines such that the non-penetration condition reads \( (\mathbf{B }\mathbf{U })_i \le 0, ~~~i = 1\ldots N_c\), the scheme may be re-written
or, as a one-step scheme
still with the addition of the complementarity conditions (29). The proof that the restitution coefficient e is asymptotically reached for a vanishing time-step is detailed in [12, 13].
Note that, even though Verlet scheme is an explicit scheme, the implicitation of the contact force in Paoli–Schatzman scheme results in a global problem to be solved on the contact surface at each time-step. Of course, this corresponds to a scalar algebraic equation which can be explicitly solved in the one-dimensional test.
Multiple impacts of an elastic bar. Paoli–Schatzman scheme for \(h=1/20\), \(\tau = 0.01\), \(e = 1/2\) and \(P_1\) Lagrange finite elements
Multiple impacts of an elastic bar. Paoli–Schatzman scheme for \(h=1/20\), \(\tau = 0.01\), \(e = 1\) and \(P_1\) Lagrange finite elements
The numerical tests for \(h=0.05\) and \(\tau =0.01\) are presented in Figs. 9, 10 and 11 for a restitution coefficient equal to 0, 1 / 2 and 1, respectively. The results for \(e=0\) and \(e=1/2\) are very similar to each other despite the difference between the restitution coefficients, and we observe a very similar loss of energy for each impact. The approximation of the displacement and of the non-penetration condition are quite good. The results for \(e=1\) show an excessive bounce of the contact point which leads to very noisy contact point velocity and contact stress.
Comparison with Taylor–Flanagan scheme
Taylor and Flanagan [5] developed an explicit scheme in the framework of PRONTO3D software which rapidly became a reference for explicit integration of contact and impact problems. To summarize the principle of the method, it is based on the Leapfrog form of Verlet scheme (17). When contact occurs, the persistency condition is prescribed at the half time-step by enforcing the relative velocity to vanish (see Eq. (21) in [6], for instance). To this aim, a Lagrange multiplier is introduced which is taken into account in an implicit way. The equation associated to the Lagrange multiplier can be solved locally only in a node-to-node contact approximation. For a more general contact condition, the Lagrange multiplier is obtained by solving a global problem on the contact surface. However, Taylor and Flanagan propose an iterative method to compute the Lagrange multiplier which needs only a few iterations.
Multiple impacts of an elastic bar. Taylor–Flanagan scheme for \(h=1/20\), \(\tau = 0.01\), and \(P_1\) Lagrange finite elements
Since the Taylor–Flanagan scheme prescribes the contact condition with an implicited Lagrange multiplier and enforces the persistency condition, it is very close to the Paoli–Schatzman scheme with a restitution coefficient \(e=0\) even if the contact condition is prescribed in a slightly different way. The consequence is that the results of the simulations shown on Fig. 12 for the Taylor–Flanagan scheme are almost identical to the results shown on Fig. 9 for the Paoli–Schatzman scheme with \(e=0\). Particularly, a loss of energy occurs at each impact.
Comparison with the mass redistribution method
The mass redistribution method, introduced in [11], considers a semi-discretization that comes from the finite element approximation of the dynamic contact problem combined with a Lagrange multiplier method to enforce the contact condition:
for a.e. \(t \in (0, T]\), still with \(\mathbf{K }\) the stiffness matrix, with \(\mathbf{B }\) the matrix representing the discrete normal trace operator on the contact boundary, and with \(\mathbf{M }_\mathrm {r}\) a modified mass matrix with a vanishing contribution on the contact boundary. The matrix \(\mathbf{M }_\mathrm {r}\) is simply built from the standard mass matrix \(\mathbf{M }\) by setting to zero the lines and columns corresponding to the degrees of freedom on the contact boundary and redistributing the removed mass on the internal degrees of freedom (see a possible redistribution algorithm in [11] and other strategies in [17, 52, 53]). In the one-dimensional test-case, this just means setting to zero the first column and first row and adding the removed mass on the first degree of freedom, which as been proved to be the most optimal strategy in [52]. It is proved in [11] that the mass redistribution method allows to recover the well-posedness of the discretization and that the solution to the approximated problem is energy conserving. Some convergence results can be found in [15, 17, 52, 54, 55]. Note also that such singular mass matrices can be obtained with the method introduced in [18] by considering different approximations for the displacement and the velocity instead of using a post-modification of the mass matrix.
Multiple impacts of an elastic bar. Mass redistribution method for \(h=1/20\), \(\tau = 0.01\) and \(P_1\) Lagrange finite elements
Since the mass matrix admits a kernel containing the vectors being only nonzero on the contact boundary, the system (32), (33) consists in an algebraic variational inequality when reduced on this kernel. Due to the Lipschitz continuity, with respect to the data, of the solution to this variational inequality, Problems (32), (33) reduces to a system of ordinary differential equations on the orthogonal of the kernel. This property, detailed in [11] allows to use quite arbitrary time-marching schemes to approximate (32), (33), among others the Verlet scheme. Of course, the method is not strictly an explicit one since a global solving has to be done on the kernel of the modified mass matrix. However, in the one-dimensional test case, this kernel is one-dimensional which allows an explicit solving.
The corresponding simulations can be seen on Fig. 13. One characteristic of the mass redistribution method is to produce low oscillating velocity and contact stress compared to other discretizations. One can see that the energy is conserved, but slightly modified compared to the standard energy.
Comparison with the penalty method
The penalty method is one of the simplest and most popular way to approximate the contact condition (see for instance [21]). With the notations used to write system (9) for Nitsche’s method and with the non-linear operator \(\mathbf{B }^h_\mathrm {p} : \mathbf{V }^h \rightarrow \mathbf{V }^h\), defined for all \(\mathbf{v }^h,\mathbf{w }^h \in \mathbf{V }^h\) by
the dynamic contact problem approximated by a penalty method reads:
This problem corresponds to a nonlinear system of ordinary differential equations on which a Verlet scheme can be applied. The parameter \(\gamma _h\) is now the penalty parameter and, following for instance the analysis in the static case of [56] and similarly to Nitsche’s method, we will still consider that 
for each finite element K.
Multiple impacts of an elastic bar. Penalty method for \(h=1/20\), \(\tau = 0.01\), \(\gamma _0 = 5\) and \(P_1\) Lagrange finite elements
Figures 14, 15 and 16 depict three simulations for the one-dimensional test case for \(\gamma _0 = 5\), \(\gamma _0 = 1\) and \(\gamma _0=0.25\), respectively. The augmented energy associated to penalty is the following one [30]:
For \(\gamma _0 = 0.25\) the contact interface is clearly too soft and a large penetration occurs, which makes the approximated solution being far from the exact one. Conversely, for \(\gamma _0 = 5\), the non-penetration condition is better approximated, but some important oscillations on the velocity of the contact point and on the contact stress occur. Similarly to Nitsche’s method, an acceptable compromise seems to set \(\gamma _0 = 1\), which corresponds to comparable stiffnesses on the contact point due to both the penalty term, on one side, and to the interior elasticity terms, on the other.
Multiple impacts of an elastic bar. Penalty method for \(h=1/20\), \(\tau = 0.01\), \(\gamma _0 = 1\) and \(P_1\) Lagrange finite elements
Multiple impacts of an elastic bar. Penalty method for \(h=1/20\), \(\tau = 0.01\), \(\gamma _0 = 0.25\) and \(P_1\) Lagrange finite elements
It is worth comparing Fig. 15 to Figs. 2, 3 and 4 for Nitsche’s method and the same value of the parameter \(\gamma _0\). The non-penetration condition is better satisfied with Nitsche’s method, which highlights its consistency. However, energy conservation is better preserved by the penalty method except when the variant \(\Theta =1\) of Nitsche’s method is used.
Numerical convergence
Simulations in the previous sections allow a qualitative comparison of the different studied methods on the one-dimensional test-case. The aim of this section is to complete this comparison with a convergence study, still for the same test-case. This is done for both linear finite elements (\(P_1\) Lagrange) on Fig. 17 and quadratic finite elements (\(P_2\) Lagrange) on Fig. 18. For M the number of elements, \(h=1/M\) is the element size and the time-step is chosen to be \(\tau = h/10\) for \(P_1\) elements and \(\tau = h/20\) for \(P_2\) elements.
Multiple impacts of an elastic bar. Convergence tests for \(P_1\) Lagrange finite elements
A comparison of Figs. 17 and 18 leads to the conclusion that despite the very low regularity of the exact solution (velocity and stress are discontinuous), there is a substantial gain in using quadratic elements. It even improves the convergence rate for the \(L^2(0,T,L^2(\Omega ))\)-norm of the error of the displacement. Globally, the mass redistribution with quadratic elements provides the best compromise. However, as the Paoli–Schatzman scheme, it necessitates to solve a global problem on the contact surface.
So, if we limit the comparison to primal discretizations, which do not require to solve such a global problem (except the inversion of the mass matrix if the mass matrix is not lumped), we can compare only Nitsche and penalty methods. What can be observed, especially on the \(L^2(0,T,L^2(\Omega ))\)-norm of the error in displacement, is that Nitsche’s method is less sensitive to the parameter \(\gamma _0\) due to its consistency. The difficulty to achieve a good compromise for penalty method is illustrated on Fig. 17: a low \(\gamma _0\) allows a good approximation of the contact stress, but the worst approximation of the \(L^2(0,T,H^1(\Omega ))\)-norm of the displacement for coarse meshes. Conversely, a large \(\gamma _0\) leads to a better approximation for the \(L^2(0,T,H^1(\Omega ))\)-norm but a too much oscillatory solution which prevents the \(L^2(0,T)\) convergence of the contact stress. There is, however, also a constraint for the choice of \(\gamma _0\) with Nitsche’s method because it has to be chosen sufficiently large to preserve the coercivity of the formulation (see Remark 19 and Fig. 4).
Multiple impacts of an elastic bar. Convergence tests for \(P_2\) Lagrange finite elements
2D/3D numerical experiments: multiple impacts of a disc / a sphere
Numerical experiments are then carried out in 2D and 3D, to assess the behavior of Nitsche’s method in a more realistic situation. We study the impact of a disc and a sphere on a rigid support. The physical parameters are the following: the diameter of the disc is \(D=40\), the Lamé coefficients are \(\lambda =30\) and \(\mu =30\), and the material density is \(\rho =1\). The total simulation time is \(T=120\).
The volume load in the vertical direction is set to \(\Vert \mathbf{f }\Vert = 0.1\) (gravity, oriented towards the support). On the upper part of the boundary is applied a homogeneous Neumann condition \(\mathbf{g }=\mathbf{0 }\) and the lower part of the boundary is the contact zone \(\Gamma _{C}\). There an initial vertical displacement (\(\mathbf{u }_0 = (0, 4)\)) and no initial velocity (\(\dot{\mathbf{u }}_0 = \mathbf{0 }\)). In such a situation, there is up to our knowledge no analytic solution to validate the numerical results.
\(P_2\) meshes used for the ball and the sphere
First bounce of the disc. Von Mises stress distribution
First bounce of the sphere. Von Mises stress distribution. One view over two is a sectional one
Comparison of the penalty, singular mass and Nitsche methods in the 2D case
Comparison of the penalty, singular mass and Nitsche methods in the 3D case
For space semi-discretization, Lagrange isoparametric finite elements of order \(k=2\) have been used. The mesh size is \(h=4\) for the ball and \(h=8\) for the sphere (see Fig. 19). Integrals of the non-linear term on \(\Gamma _C\) are computed with standard quadrature formulas of order 4. A snapshot of the evolution of the disc and the sphere during the first impact can be seen on Figs. 20 and 21.
The comparison of the simulations for the different methods is depicted Fig. 22 for the two-dimensional case and Fig. 23 for the three-dimensional case. For the sake of shortness, only the penalty, the singular mass and Nitsche methods are compared. First of all, a conclusion that can be drawn from these numerical experiments is that the tested methods are all capable of reliably approximating two and three-dimensional dynamic contact problems. An important difference between simulations in dimension 2 and 3 is a much smaller oscillation of the contact stress in dimension 3, except for the mass redistribution method which is not subjected to spurious oscillations. The energy is conserved more strictly with the penalty method, the mass redistribution method and the variant \(\Theta = 1\) of Nitsche’s method, the other two variants presenting significant disturbances in the energy evolution. The mass redistribution method appears to give the best compromise between energy conservation and the low level of oscillation on the contact boundary. Note however that it produces a weakening of the rebound, mainly in dimension 3, which we do not explain. The lake of consistency of the penalty method is illustrated on the normal displacement graph where we can note a larger interpenetration compared to the other methods. Finally, among the variants of Nitsche’s method, the symmetric variant \(\Theta = 1\) is the one that achieves the best compromise between energy conservation and the level of oscillations of the contact stress, which remains moderate.
Concluding remarks
In this paper, we studied the application of an explicit Verlet scheme for the approximation of elastodynamic contact problems with Nitsche’s method. The explicit method being commonly used in elastodynamic contact problems, it seemed important to complete the study that had been performed in [28, 29] for implicit schemes. We tried to characterize the stability properties of the different variants of Nitsche’s method for the Verlet scheme and we introduced a number of necessary tools for this analysis. Of course, we are aware that the stability result we establish is very partial (only for \(\Theta = 1\)) and certainly suboptimal: a stability condition such as \(\tau = \mathcal {O}(h)\) would be more satisfactory and would correspond to what we noted in numerical tests. This result remains to be refined. Moreover, it would certainly be possible to prove a convergence result in dimension one, as in [15], because in this context the existence and uniqueness of the solution is theoretically proven.
We numerically compared the Nitsche method to the main existing methods that can support an explicit scheme: the Paoli–Schatzman scheme, the Taylor–Flanagan scheme, the mass redistribution method and the penalty method.
We first performed this comparison on a one-dimensional test case whose exact solution consists of a shock wave indefinitely travelling between the two ends of a bar. We can see globally, by comparing the Figs. 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 and 16 that Nitsche’s method, especially the variant corresponding to \(\Theta = 1\), yields an approximation of a comparable quality as the one obtained with the schemes using an implicitation of the contact force (Paoli–Schatzman scheme, Taylor–Flanagan scheme, mass redistribution). Only the mass redistribution method results in lower oscillation levels. Compared to the penalty method, the oscillation levels are of the same magnitude, but contact penetration is more limited. This is reflected in Figs. 17 and 18 by quite similar convergence rates for all the methods, except for penalty. For this latter, a compromise remains difficult to find between a large penalization coefficient, which corresponds to a good approximation of the displacement but a poor approximation of the contact stress, and a small penalization coefficient, for which the interpenetration becomes too large. The decisive advantage of Nitsche’s method over the other methods, with the exception of the penalty method, is to be a primal method for which there is no need for an implicit resolution of the contact force. This allows a really explicit resolution in case of lumped mass matrix. The 2D and 3D test cases we performed also confirm the good behavior of Nitsche’s method. We can see in Figs. 22 and 23 the advantage in comparison to the penalty method in term of interpenetration, which is less. Still some better approximation results are obtained for the variant \(\Theta = 1\).
We can thus conclude that, among the variants of Nitsche’s method, the symmetric variant \(\Theta =1\) seems to be the most suitable for solving dynamic contact problems mainly because of its energy conservation properties. For the other variants a gain of energy can be observed, especially for low values of Nitsche’s parameter \(\gamma _0\). Some perspectives of this work could be to gain further insight into the properties of energy conservation, for instance using other definitions of the discrete energies, such as in [57, Theorem 4.1, Remark 4.1] in which an energy that remains positive irrespectively of the value of numerical parameters is introduced. Also some new explicit time-marching schemes endowed with appealing properties of energy conservation could be considered (see, e.g., [58]). Moreover, further study of the effect of the mass matrix lumping, particularly on the stability of the method, and of the proper choice of the Nitsche’s parameter \(\gamma _0\) are other perspectives of this work.
References
Cirak F, West M. Decomposition contact response (DCR) for explicit finite element dynamics. Int J Numer Methods Eng. 2005;64(8):1078–110.
Ryckman RA, Lew AJ. An explicit asynchronous contact algorithm for elastic body-rigid wall interaction. Int J Numer Methods Eng. 2012;89(7):869–96.
Wolff S, Bucher C. Asynchronous collision integrators: explicit treatment of unilateral contact with friction and nodal restraints. Int J Numer Methods Eng. 2013;95(7):562–86.
Stewart JR, Gullerud AS, Heinstein MW. Solution verification for explicit transient dynamics problems in the presence of hourglass and contact forces. Comput Methods Appl Mech Eng. 2006;195(13–16):1499–516.
Taylor LM, Flanagan DP. PRONTO3D: a three-dimensionnal transient solid dynamics program. Albuquerque: Sandia National Laboratories; 1989. p. 87185.
Heinstein MW, Mello FJ, Attaway SW, Laursen TA. Contact-impact modeling in explicit transient dynamics. Comput Methods Appl Mech Eng. 2000;187:621–40.
Moreau JJ. Unilateral contact and dry friction in finite freedom dynamics. Vienna: Springer; 1988. p. 1–82.
Moreau JJ. Numerical aspects of the sweeping process. Comput Methods Appl Mech Eng. 1999;177:329–49.
Vola D, Raous M, Martins JAC. Friction and instability of steady sliding: squeal of a rubber/glass contact. Int J Numer Methods Eng. 1999;46(10):1699–720.
Newmark NM. A method of computation for structural dynamics. J Eng Mech Div. 1959;85(3):67–94.
Khenous HB, Laborde P, Renard Y. Mass redistribution method for finite element contact problems in elastodynamics. Eur J Mech. 2008;27(5):918–32.
Paoli L, Schatzman M. A numerical scheme for impact problems. I. The one-dimensional case. SIAM J Numer Anal. 2002;40(2):702–33.
Paoli L, Schatzman M. A numerical scheme for impact problems. II. The multidimensional case. SIAM J Numer Anal. 2002;40(2):734–68.
Schindler T, Acary V. Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook. Math Comput Simulation. 2014;95:180–99.
Dabaghi F, Petrov A, Pousin J, Renard Y. A robust finite element redistribution approach for elastodynamic contact problems. Appl Numer Math. 2016;103:48–71.
Lebeau G, Schatzman M. A wave problem in a half-space with a unilateral constraint at the boundary. J Differ Equ. 1984;53(3):309–61.
Hager C, Hüeber S, Wohlmuth BI. A stable energy-conserving approach for frictional contact problems based on quadrature formulas. Int J Numer Methods Eng. 2008;5(27):918–32.
Renard Y. The singular dynamic method for constrained second order hyperbolic equations: application to dynamic contact problems. J Comput Appl Math. 2010;234(3):906–23.
Doyen D, Ern A, Piperno S. Quasi-explicit time-integration schemes for dynamic fracture with set-valued cohesive zone models. Comput Mech. 2013;52(2):401–16.
Belhytschko T, Neal MO. Contact-impact by the pinball algorithm with penaly and lagrangian methods. Internat J Numer Methods Eng. 1991;31:547–72.
Kikuchi N, Oden JT. Contact problems in elasticity: a study of variational inequalities and finite element methods. In: SIAM studies in applied mathematics. vol. 8. Society for Industrial and Applied Mathematics (SIAM), Philadelphia; 1988. p. 495.
Doyen D, Ern A, Piperno S. Time-integration schemes for the finite element dynamic Signorini problem. SIAM J Sci Comput. 2011;33(1):223–49.
Nitsche J. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg. 1971;36:9–15.
Stenberg R. On some techniques for approximating boundary conditions in the finite element method. J Comput Appl Math. 1995;63(1–3):139–48.
Chouly F, Hild P. A Nitsche-based method for unilateral contact problems: numerical analysis. SIAM J Numer Anal. 2013;51(2):1295–307.
Chouly F, Hild P, Renard Y. Symmetric and non-symmetric variants of Nitsche’s method for contact problems in elasticity: theory and numerical experiments. Math Comp. 2015;84(293):1089–112.
Chouly F, Fabre M, Hild P, Mlika R, Pousin J, Renard Y. An overview of recent results on Nitsche’s method for contact problems. Lect Notes Comput Sci Eng. 2018;121:93–141.
Chouly F, Hild P, Renard Y. A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes. ESAIM Math Model Numer Anal. 2015;49(2):481–502.
Chouly F, Hild P, Renard Y. A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments. ESAIM Math Model Numer Anal. 2015;49(2):503–28.
Hauret P, Le Tallec P. Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact. Comput Methods Appl Mech Eng. 2006;195(37–40):4890–916.
Hairer E, Lubich C, Wanner G. Geometric numerical integration, 2nd edn. Springer Series in computational mathematics. Berlin : Springer; 2006. p. 31
Adams R-A. Sobolev spaces. Pure and applied mathematics, vol. 65. New York: Academic Press; 1975.
Schatzman M. A hyperbolic problem of second order with unilateral constraints: the vibrating string with a concave obstacle. J Math Anal Appl. 1980;73(1):138–91.
Schatzman M. Un problème hyperbolique du 2ème ordre avec contrainte unilatérale: la corde vibrante avec obstacle ponctuel. J Differ Equ. 1980;36(2):295–334.
Kim JU. A boundary thin obstacle problem for a wave equation. Comm Partial Differ Equ. 1989;14(8–9):1011–26.
Dabaghi F, Petrov A, Pousin J, Renard Y. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary. M2AN Math Model Numer Anal. 2014;48:1147–69.
Ahn J, Stewart DE. Existence of solutions for a class of impact problems without viscosity. SIAM J Math Anal. 2006;38(1):37–63.
Pozzolini C, Renard Y, Salaün M. Vibro-impact of a plate on rigid obstacles: existence theorem, convergence of a scheme and numerical simulations. IMA J Numer Anal. 2013;33(1):261–94.
Eck C, Jaru\(\check{s}\)ek J, Krbec M. Unilateral contact problems. Pure and applied mathematics (Boca Raton). Boca Raton : Chapman & Hall/CRC; 2005. vol. 270, p. 398.
Laursen TA, Chawla V. Design of energy conserving algorithms for frictionless dynamic contact problems. Int J Numer Methods Eng. 1997;40(5):863–86.
Armero F, Petőcz E. Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems. Comput Methods Appl Mech Eng. 1998;158(3–4):269–300.
Ciarlet P-G. The finite element method for elliptic problems. Handbook of numerical analysis. Ciarlet PG, Lions JI. vol. II. Amsterdam : North-Holland Publishing Co.; 1991.
Thomée V. Galerkin finite element methods for parabolic problems. Springer Series in computational mathematics. vol. 25. Berlin: Springer; 1997. p. 302.
Alart P, Curnier A. A generalized Newton method for contact problems with friction. J Mech Theor Appl. 1988;7(1):67–82.
Chouly F. An adaptation of Nitsche’s method to the Tresca friction problem. J Math Anal Appl. 2014;411(1):329–39.
Cohen G, Pernet S. Finite element and discontinuous Galerkin methods for transient wave equations. Scientific Computation. Dordrecht: Springer; 2017. p. 381
Ern A, Guermond J-L. Theory and practice of finite elements. Applied mathematical sciences, vol. 159. New York: Springer; 2004.
Brenner S-C, Scott L-R. The mathematical theory of finite element methods. Texts in applied mathematics, vol. 15. New York: Springer; 2007.
Paoli L, Schatzman M. Schéma numérique pour un modèle de vibrations avec contraintes unilatérales et perte d’énergie aux impacts, en dimension finie. C R Acad Sci Paris Sér I Math. 1993;317(2):211–5.
Paoli L, Schatzman M. Approximation et existence en vibro-impact. C R Acad Sci Paris. 1999;317:1103–7.
Paoli L, Schatzman M. Numerical simulation of the dynamics of an impacting bar. Comput Methods Appl Mech Eng. 2007;196:2839–51.
Dabaghi F, Krejči P, Petrov A, Pousin J, Renard Y. A weighted finite element mass redistribution method for dynamic contact problems. J Comp Appl Math. 2019;345:338–56.
Wohlmuth BI. Variationally consistent discretization schemes and numerical algorithms for contact problems. Acta Numer. 2011;20:569–734.
Doyen D, Ern A. Convergence of a space semi-discrete modified mass method for the dynamic Signorini problem. Commun Math Sci. 2009;7(4):1063–72.
Doyen D, Ern A. Analysis of the modified mass method for the dynamic Signorini problem with Coulomb friction. SIAM J Numer Anal. 2011;49(5):2039–56.
Chouly F, Hild P. On convergence of the penalty method for unilateral contact problems. Appl Numer Math. 2013;65:27–40.
Burman E, Fernández MA, Frei S. A Nitsche-based formulation for fluid-structure interactions with contact. Research Report RR-9172, Inria (2018). hal-01784841. https://hal.inria.fr/hal-01784841
Marazzato F, Ern A, Mariotti C, Monasse L. An explicit energy-momentum conserving time-integration scheme for Hamiltonian dynamics. hal-01661608 (2017). https://hal-enpc.archives-ouvertes.fr/hal-01661608
Authors' contributions
The two authors collaborated at each stage of this study: obtaining mathematical results, developing test codes, analyzing numerical results, drafting the manuscript and finalizing it. Both authors read and approved the final manuscript.
Acknowledgements
We thank both referees for their constructive comments that helped to improve the presentation of the results. We thank Patrick Letallec for the rewarding discussions on existing explicit strategies in contact dynamics.
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
The datasets used and analysed during the current study are available from the corresponding author on reasonable request.
Funding
Not applicable.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Chouly, F., Renard, Y. Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems. Adv. Model. and Simul. in Eng. Sci. 5, 31 (2018). https://doi.org/10.1186/s40323-018-0124-5
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s40323-018-0124-5