{'fraction_non_alphanumeric': 0.04421331463519146, 'fraction_numerical': 0.016945191071637257, 'mean_word_length': 4.123129115944921, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A weighted finite element mass redistribution method for dynamic contact problems

This paper deals with a one-dimensional wave equation being subjected to a unilateral boundary condition. An approximation of this problem combining the finite element and mass redistribution methods is proposed. The mass redistribution method is based on a redistribution of the body mass such that there is no inertia at the contact node and the mass of the contact node is redistributed on the other nodes. The convergence as well as an error estimate in time are proved. The analytical solution associated with a benchmark problem is introduced and it is compared to approximate solutions for different choices of mass redistribution. However some oscillations for the energy associated with approximate solutions obtained for the second order schemes can be observed after the impact. To overcome this difficulty, an new unconditionally stable and a very lightly dissipative scheme is proposed.


Introduction
The present paper highlights some new numerical results obtained for a one-dimensional elastodynamic contact problem. Dynamical contact problems play a crucial role in structural mechanics as well as in biomechanics and a considerable amount of engineering and mathematical literature has been dedicated to this topic last decades. One of the main difficulties in the numerical treatment of such problems is the physically meaningful non-penetration condition that is usually modeled by using the so-called Signorini boundary condition. Basically, the lack of well-posedness results mainly originates from the hyperbolic structure of the problem which gives rise to shocks at the contact interfaces. Then the resulting nonsmooth and nonlinear variational inequalities lead to fundamental difficulties in mathematical analysis as well as in the development of numerical integration schemes. In view to avoid these difficulties, the non-penetration condition is quite often relaxed in the numerical integration schemes. We may also observe that most of unconditionally stable schemes for the linear elastodynamic problems lose their unconditional stability in the presence of contact conditions. Among them the classical Newmark method is the most popular one. However, its unsatisfactory handling of the non-penetration conditions may lead to artificial oscillations at the contact boundary and even give rise to an undesirable energy blow-up during the time integration, the reader is referred to [KLR08, DEP11, DP*13] as well as to the references therein for further details. To overcome these difficulties, some numerical methods based on the Newmark scheme for solving impact problems are proposed in [CTK91]. However, these methods lead to some important energy losses when the contact takes place even if the time step is taken sufficiently small. On the other hand, the energy conserving time integration schemes of Newmark type are introduced in [LaL02, LaC97,ChL98] as well as in the monograph [Lau03], but these schemes are unable to circumvent the undesirable oscillations at the contact boundary. These unphysical oscillations are avoided by the numerical methods developed in [DKE08] but these methods are still energy dissipative.
Another approach consists in removing the mass at the contact nodes and it was originally investigated in [KLR08] and later on used in [HHW08,Hau10,Ren10,CHR14]. This approach prevents the oscillations at the contact boundary and leads to well-posed and energy conserving semi-discretization of elastodynamic contact problems (see [LiR11,DP*12]). However, some numerical experiments, exhibited in [DP*13], highlight a phase shift in time between analytical and approximate solutions. Note that an analytical piecewise affine and periodic solution to our problem can be obtained by using the characteristics method while approximate solutions are exhibited for different time discretizations. This phase shift in time comes from the removed mass at the contact nodes for approximate problems which is unacceptable for many applications. Therefore a variant of the mass redistribution method is proposed in this work. More precisely, this new method consists in transferring the mass of the contact node on the other nodes meaning that the total mass of the considered material is preserved. Numerical experiments presented in this work show that the undesirable phase shift between the approximate and analytical solutions disappears and all the properties of the mass redistribution method mentioned above are preserved. They highlight that the weighted mass redistribution method is particularly well adapted to deal with contact problems.
The paper is organized as follows. In Section 2, the mathematical formulation of a one dimensional elastodynamic contact problem is presented. The contact is modeled by using the Signorini boundary conditions in displacement, which are based on a linearization of the physically meaningful non penetrability of the masses. Then a space semi-discretization based on a variant of the mass redistribution method is presented in Section 3. This variant of the mass redistribution method consists in transferring the mass of the contact node on the other nodes while the inertia vanishes at the contact node. The error estimate in time as well as the convergence result are established. A benchmark problem is introduced in Section 4 and its analytical solution is exhibited. Then numerical experiments for some space-time discretizations like the Crank-Nicolson or the backward Euler methods are reported. These numerical experiments highlight that the choice of the nodes where the mass is transferred plays a crucial role to get a better approximate solution. However some oscillations for the energy associated with approximate solutions for the second order schemes like Crank-Nicolson scheme can be observed after the impact. To overcome the difficulty, a hybrid scheme mixing the Crank-Nicolson as well as the midpoint methods and having the properties to be an unconditionally stable scheme is proposed in Section 5.

Mathematical formulation
The motion of an elastic bar of length L which is free to move as long as it does not hit a material obstacle is studied, see Figure 1. The assumptions of small deformations are assumed and the material of the bar is supposed to be homogeneous. Let u(x, t) be the displacement at time t ∈ [0, T ] , T > 0 of the material point of spatial coordinate x ∈ [0, L] . Let f (x, t) denotes a density of external forces, depending on time and space. The mathematical problem is formulated as follows: with Cauchy initial data and Signorini and Dirichlet boundary conditions at x = 0 and x = L , respectively, Here u t def = ∂u ∂t and u x def = ∂u ∂x . The orthogonality has a natural meaning: an appropriate duality product between two terms of relation vanishes. It can be alternatively stated as the inclusion where I [0,+∞) is the indicator function of the interval [0, +∞) , and ∂I [0,∞) is its subdifferential. Let us describe the weak formulation associated with (2.1)-(2.3). For that purpose, it is convenient to introduce the following notations: Thus the weak formulation associated with (2.1)-(2.3) obtained by multiplying (2.1) by v − u and by integrating formally this result over Q T def = (0, L) × (0, T ) reads: Theorem 2.1. Let u 0 ∈ H 3/2 (0, L) ∩ V , v 0 ∈ H 1/2 (0, L) , f ∈ H 3/2 (Q T ) be given. Then there exists a unique solution u ∈ L ∞ (0, T ; H 3/2 (0, L) ∩ V ) ∩ W 1,∞ (0, T ; H 1/2 (0, L)) of Problem (2.5) and the energy balance equation holds for all τ ∈ [0, T ] .
Existence and uniqueness results are obtained for a similar situation of a vibrating string with concave obstacle in one dimensional space in [Sch80] and also for a wave equation with unilateral constraint at the boundary in a half-space of R N in [LeS84]. An existence result for a wave equation in a C 2 -regular bounded domain constrained by an obstacle at the boundary in R 2 is proven in [Kim89]. The reader is also referred to [DP*12].

Finite element discretization and convergence of the mass redistribution method
This section is devoted to semi-discrete problems in space associated with (2.5) by using the mass redistribution method, see [KLR08, DP*12], assuming that the hypotheses of Theorem 2.1 are satisfied. More precisely, the weighted mass redistribution method consists in transferring the mass of the contact node on the other nodes implying that the node at the contact boundary evolves in a quasi-static way. To this aim, we choose an integer m > 1 , and put h def = L m (mesh size) with the goal to let m tend to +∞ . We introduce the spaces V h where P 1 is the space of polynomials of degree less than or equal to 1. We consider the following discretized problem: where u 0 h and v 0 h belong to V h and they are the approximations of the initial displacement and velocity u 0 and v 0 , respectively, and λ h is the Lagrange multiplier representing the contact force. The inclusion in (P u h ) (cf. also (2.4)) can be written as a variational inequality in the form The approximation u h are taken in the following form where the basis functions ϕ k are assumed piecewise linear, namely we have It is convenient for numerical computations to redistribute the mass and modify the problem (P u h ) as follows: where w h are weight functions which converge to 1 in suitable sense as h tends to 0 . We choose them to be piecewise constant for all i = 0, · · · , m − 1 . This is a problem of the type The symmetric matrices M = (M ik ) and S = (S ik ) can be computed directly from the formulas In matrix representation, we have Note that M and S are usually called mass and stiffness matrices, respectively. Consider first the problem (P mod* u h ) for i = 0 . We have (3.1) For w 0 > 0 , this produces oscillations of u 0 which are not observed in the limit. To eliminate these unphysical oscillations which are purely due to the numerical method, we assume w 0 = 0 , so that (3.1) becomes (note that where u + 1 denotes the positive part of u 1 . Then for i = 1 , we obtain from (P mod* u h ) that We have thus eliminated the singularities and problem (P mod* u h ) can be equivalently stated as with a Lipschitz continuous nonlinearity on the right hand side, and with matrices Note that (P mod** u h ) is related to a more general problem of Fučík spectrum (or 'jumping nonlinearity' in the old terminology); the reader is referred to [Fuc76] for further details. Furthermore, we may observe that (P mod** u h ) can be rewritten as follows: . . , f m−1 ) T , U 0 and V 0 approximate the initial position and velocity. Finally, the discrete energy associated with problem (P mod U h ) is given by We assume that the weights w i are chosen in such a way that M * is invertible and the matrix norm |||(M * ) −1 ||| of its inverse (M * ) −1 is bounded above by a constant independent of h . Below, we consider the following situations: In these cases, the condition on (M * ) −1 is satisfied.
Under this hypothesis, ( P mod U h ) can be rewritten as follows: Existence and uniqueness results for the problem (P mod U h ) follow from the Lipschitz continuity of G(U h ) , for further details the reader is referred to [CrM84]. In particular, we have Lemma 3.1. Let N ∈ N be given and let ∆t = T N be the time step. Then the time discretization error for the Crank-Nicolson method to solve the semi-discrete problem (P mod U h ) is of the order ∆t .
Proof. Keeping h fixed, we define discrete times t n = n∆t for n = 0, . . . , N and define the Crank-Nicolson discretization of Problem (P mod U h ) by the recurrent formula continuously extended to t N = T . We have by (3.6) that We cannot expect to obtain a higher order estimate, since G is not continuously differentiable because of the presence of the term u + 1 . On the other hand, U h is of class C 1 , and we may denote We thus have for a. e. t ∈ (0, T ) . By virtue of (3.5) we obtain for ∆t < h/C that and the assertion follows from the Gronwall argument.
The next goal is to prove the convergence of solutions to Problem (P mod u h ) (in the form (P mod U h ) ) as h → 0 . We first observe that (P mod u h ) is equivalent to We assume that the initial data u 0 h and v 0 h satisfy The convergence of the solution u h of (P mod var ) to the solution of (2.5) is proved below. To this aim, the same techniques detailed in the proof of Theorem 4.3 in [DP*12] are used. Here, we allow for general weight functions including the above cases (Mod 1)-(Mod 3), while in [DP*12], only the case (Mod 1) was considered. The reader is also referred to [ScB89].
Choosing in (P mod u h ) the test function v h (x) = u h,t (x, t) , we see that for all τ ∈ [0, T ] the following energy relation holds.
Theorem 3.2. Assume that the hypotheses of Theorem 2.1 and condition (3.11) hold. Let there exist two constants C w > c w > 0 such that c w ≤ w h (x) ≤ C w for all h > 0 and a. e. x ∈ (h, L) , and let Then the solutions u h of (P mod var ) converge in the strong topology of V to the unique solution u of (2.5) as h tends to 0 .
Proof. The energy relation (3.12) and the Gronwall lemma imply the existence of a constant hence, by virtue of (3.3), with a constant C independent of h . We thus have We conclude that there existsū ∈ W and a subsequence, still denoted by u h , such that Then, we may deduce from (3.16) that Notice that for all α < 1 2 , we have W → C 0, 1 2 (Q T ) → → C 0,α (Q T ) hold (see [ScB89]), where → and → → denote the continuous and compact embeddings, respectively. Finally we find for all α < 1 2 . Furthermore, u h andū belong to K . Our aim is to establish thatū = u , that is, the limitū coincides with the solution u of (2.5). However the elements of K are not smooth enough in time, then they should be approximated before being projected onto V h . Indeed this projection violates the constraint at x = 0 , and therefore, the elements of K need another approximation in order to satisfy the constraint strictly.
Assume that v ∈ K is an admissible test function for (2.5), that is, v =ū for t ≥ T − ζ . For η ≤ ζ/4 we define an auxiliary function where ψ is a smooth and positive function with the property ψ = 1 on [0, T − η/2] and ψ = 0 on [T − η/4, T ] . We precise now how the parameter k(η) to ensure that v η ∈ K ∩ L ∞ (0, T ; V) holds. Sinceū ∈ C 0, 1 2 (Q T ) , it follows that there exists a constant C * > 0 such that Then for all Let D h : V → V h be the piecewise linear interpolation mapping defined by the formula (3.23) By (3.21), we have for h → 0 the strong convergences We now use (3.16) and (3.24) to pass to the limit as h → 0 in (3.23) and find that (3.25) The passage to the limit in (3.25) as η → 0 is easy, and we conclude thatū = u is the desired solution of (2.5). It remains to prove that u h converge strongly in V . Passing to the limit as h → 0 in (3.12) we obtain for a. e. τ ∈ (0, T ) that (3.26) Since u h u weakly-* in W and w h → 1 strongly in every L p with p < ∞ , we conclude that u h,t √ w h converge weakly to u t in L 2 (0, T ; H) . Since the norms of (u h,t √ w h , u h,x ) in L 2 (0, T ; H) converge to the norm of (u t , u x ) in L 2 (0, T ; H) , we see that u h converge strongly to u in V . The limit solution u is unique, hence the whole system {u h : h > 0} converges to u as h → 0 , which completes the proof.

The wave equation with Signorini and Dirichlet boundary conditions
We consider a bar of length L = 1 clamped at one end and compressed at t = 0 . The bar elongates under the elasticity effect; as soon as it reaches a rigid obstacle at time t 1 then it stays in contact during the time t 2 − t 1 and it takes off at time t 2 , see Figure 2. This problem can be described mathematically by (2.1)-(2.3) with the density of external forces f (x, t) = 0 . We first describe how an analytical piecewise affine and periodic solution to our problem can be obtained by using the characteristics method, the reader is referred to [DP*13] for a detailed explanation. Then approximate solutions for some time-space discretizations and for several mass redistributions are exhibited and their efficiency are discussed.

Analytical solution
The domains considered here are defined by (0, L) × (t i , t i+1 ) , i = 0, 1, 2 , corresponding to the phases before, during and after the impact, respectively. Each of them are divided into four regions as it is represented on Figure 2. We choose below t i = i , i = 0, . . . , 3 .  III  III  III  I  I  I   contact  before  during  contact  after    The domain (2, 3) × (0, 1) corresponding to the phase after the impact is split into four regions. By using (4.2), we get u(·, 2) = 0 and u t (·, 2) = 1 2 which leads to u(x, t) = t−t 2 2 in the regions I, II, 1−x 2 in the regions III, IV.

Comparisons between different mass redistributions for some time-space discretizations
The time discretization is introduced in this section. To this aim, we divide the time interval [0, T ] by n + 1 discrete time-points such that 0 = t 0 < t 1 < . . . < t n = T . Let U n h ,U n h ,Ü n h and λ n be the approximations of the displacement U h (t n ) , the velocityU h (t n ) , the accelera-tionÜ h (t n ) and the Lagrange multiplier λ(t n ) , respectively. We deal with some approximate solutions to Problem (2.1)-(2.3) obtained by using several time-stepping methods like the Newmark, backward Euler and Paoli-Schatzman methods. For each of these time-stepping methods, approximate solutions (U n h , λ n ) are exhibited for several mass redistributions and they are compared to the analytical solution (u, λ) introduced in Section 4.1. In the numerical experiments presented below, we distinguish the cases that the mass of the contact node that is not redistributed (Mod 1), or uniformly redistributed on all the other nodes (Mod 2), or redistributed only on the nearest neighbor (Mod 3), according to the classification given in the previous section. We show that the efficiency of the mass redistribution method depends on the position of the nodes where the mass is redistributed. Indeed, the numerical experiments highlight that the closer from the contact node the mass is transferred better the approximate solutions are obtained. Then, it is not surprising that the best approximate solution can be expected and indeed it is obtained in the case (Mod 3), where all the mass of the contact node is transferred on the node preceding the contact node, see Figures 3, 5 and 7. Note that the numerical simulations presented below were performed by employing the finite element library Getfem ++ (see [ReP]).

The Newmark methods
The Taylor expansions of displacements and velocities neglecting terms of higher order are the underlying concept of the family of Newmark methods, see [New59]. These methods are unconditionally stable for linear elastodynamic problem for γ ≥ 1 2 and β ≥ 1 4 ( 1 2 +γ) 2 , see [Hug87,Kre06], but they are also the most popular time-stepping schemes used to solve contact problems. The discrete evolution for the contact problem (2.1)-(2.3) is described by the following finite difference equations: where ∆t is a given times step and (β, γ) are the algorithmic parameters, see [Hug87,Lau03]. Note that U 0 h ,U 0 h and λ 0 are given andÜ 0 h is evaluated by using the third identity in (4.3). We are particularly interested in the case where (β, γ) = ( 1 4 , 1 2 ) . This method is called the Crank-Nicolson method, it is second-order consistent and unconditionally stable in the unconstrained case. However the situation is quite different in the case of contact constraints, indeed the order of accuracy is degraded; for further details, the reader is referred to [Hug87,Kre06,GrH07]. The analytical solution (u, λ) exhibited in Section 4.1 and the approximate solutions (U n h , λ n ) obtained for different mass redistributions are represented on Figure 3. The approximate solution obtained for the nearest neighbor redistribution (Mod 3 on Figure 3) gives much better accuracy than the mass redistribution on all the nodes preceding the contact node (see Mod 2 on Figure 3) and no mass redistribution. This highlighted that the choice for the mass redistribution plays a crucial role.

The backward Euler method
We are concerned here with backward Euler's method which for the contact problem (2.1)-(2.3) is described by the following finite difference equations:

The Paoli-Schatzman methods
We focus on the so-called Paoli-Schatzman method that consists to fix the contact constraint at an intermediate time step. Indeed the method proposed below is a slight modification of Paoli-Schatzman method (see [Pao01,PaS02]) which takes into account the kernel of the modified mass matrix. A simple application of Paoli-Schatzman method based on Newmark scheme to our problem with γ = 1 2 leads to U 0 and U 1 given. (4.5) Here e belongs to [0, 1] and is aimed to be interpreted as a restitution coefficient. Note that U 0 h and U 1 h are given data and U 1 h can be evaluated by a one step scheme. We may observe that taking M = M mod in (4.5), we are not able to resolve the problem on the kernel of M mod . That is the reason why, SU n−1 h as well as SU n h are projected on the orthogonal of the kernel of M . We are interested here in Paoli-Schatzman's method with (β, e) = ( 1 4 , 1) . Note that the stability result immediately follows from [DuP06].

A hybrid time integration scheme
Generally, the second order schemes illustrate one of the difficulties when solving contact problems, namely some oscillations for energy associated with approximate solutions obtained for different choices of mass redistribution can be observed after each impact takes place, for instance see Figure 4. To overcome this problem, a hybrid time integration scheme is introduced in this section. More precisely, the scheme (P mod U h ) is modified to be an unconditionally stable and a second order in time scheme; the linear part of (P mod U h ) is discretized by using the midpoint method while the non-linear part is discretized by using the Crank-Nicolson as well as the midpoint methods. Observe that the midpoint method for the linear problem is energy conserving. The proposed hybrid time integration scheme inspired from [CHR14] reads as follows: and H is defined by We assume that the density of external forces F does not depend on time. Observe that 1 2h ((u n 1 ) + + (u n+1 1 ) + )e 1 and 1 2h (u n 1 + u n+1 1 ) + e 1 correspond to the contribution of Crank-Nicolson and midpoint methods, respectively.
The discrete evolution of the total energy is preserved in the purely elastic case when the density of external forces vanishes, see [Lau03]. However, the situation quite different in the case of contact constraints, the order of accuracy is degraded and , for further details see [Hug87,Kre06,GrH07]. Let us define now the energy evolution by ∆E n where E n h is assumed to be given by an algorithmic approximation of the energy E h (t n ) defined in (3.4). We evaluate now ∆E n h by using the midpoint scheme. More precisely, we get Since M * is symmetric matrix, it comes that , which implies that We establish below that the energy evolution by ∆E n h is nonpositive, namely the energy associated with the hybrid scheme decreases in time. This result is summarized in the following lemma: Lemma 5.1. Assume that the density of external forces F does not depend on time. Then the energy evolution ∆E n h is nonpositive for all n > 0 .
5. If u n 1 = 0 then ∆E n h = 1 2h (u n+1 1 )(u n+1 This proves the lemma. The numerical experiments presented on Figures 9 and 10 are obtained by using a new hybrid scheme where the mass of the contact node is redistributed on the node preceding the contact. This scheme circumvents the undesirable oscillations at the contact boundary and it prevents as well the small oscillations of the evolution of total energy occurring for Newmark methods (see Figures 4 and 8). Finally, the energy evolution ∆E n h is nonpositive for all n > 0 and it is much smaller than the energy evolution obtained by using implicit Euler method (compare with Figure 6).

Conclusion
This manuscript focuses on the weighted mass redistribution method which is particularly well adapted to approximate elastodynamic contact problems. This method leads to well-posed and energy conserving semi-discretization of elastodynamic contact problems. Furthermore, it prevents some undesirable oscillations at the contact boundary as well as some phase shift between approximate and analytical solutions. The efficiency of the weighted mass redistribution method depends on the position of the nodes where the mass is redistributed; the closer the mass of the contact node is transferred, the better are the approximate solutions. These results seem also valid in higher space dimensions, and in particular in 2D space (see Table 1). Indeed the weight mass redistribution on the nodes before the contact (Mod 3) gives much better absolute error than in the cases where any redistribution is done (Mod 1) or where the mass is just eliminated from the contact nodes (Mod 2). The total error rates are evaluated for the space steps ∆x 1 = ∆x 2 = 0.05 . However the energy associated with the Newmark and Paoli-Schatzman methods in time have small oscillations, (see for instance Figure 4), which is unacceptable from a mechanical view point. Then a new hybrid scheme having the properties to be an unconditionally stable has been developed giving some promising numerical results. It allow to have a far better approximation compared to existing unconditionally stable scheme on the implicit Euler scheme (compare Figures 6 and 10).