The Wave Finite Element Method applied to a one-dimensional linear elastodynamic problem with unilateral constraints

The Wave Finite Element Method (WFEM) is implemented to accurately capture travelling waves propagating at a ﬁnite speed within a bouncing rod system and induced by unilateral contact collisions with a rigid foundation; friction is not accounted for. As opposed to the traditional Finite Element Method (FEM) within a time-stepping framework, potential discontinuous deformation, stress and velocity wave fronts are accurately described, which is critical for the problem of interest. A one-dimensional benchmark with an analytical solution is investigated. The WFEM is compared to two time-stepping solution methods formulated on a FEM semi-discretization in space: (1) an explicit technique involving Lagrange multipliers and (2) a non-smooth approach implemented in the Siconos package. Attention is paid to the Gibb’s phenomenon generated during and after contact occurrences together with the time evolution of the total energy of the system. It is numerically found that the WFEM outperforms the FEM and Siconos solution methods because it does not induce any spurious oscillations or dispersion and diffusion of the shock wave. Furthermore, energy is not dissipated over time. More importantly, the WFEM does not require any impact law to close the system of governing equations.


Introduction
Collisions between structural components can be observed in a variety of industrial processes, such as rock breaking equipment, vehicles crash testing, drilling tools, and many others.Through contact forces preventing non-physical inter-penetration of matter, a collision initiates a disturbance in the form of a stress wave which travels away from the region of contact.Corresponding governing equations are derived by the laws of mass conservation, momentum conservation, and energy conservation.Additional complementarity conditions are incorporated to properly reflect the impenetrability of bodies.Analytical solutions are only known for simple cases such as longitudinal collision of simple rods or transverse impacts on beams.For more general and challenging configurations, the solution should be numerically approximated.
As numerical procedures produce approximate solutions, they inherently exhibit limitations that might be unacceptable.For instance, the solution may feature non-physical spurious oscillations, commonly known as Gibb's phenomenon.The approximate wave propagation velocity may also be different from its physical counterpart due to period lengthening and amplitude decrease resulting in dispersion and dissipation errors [1].These potential numerical discrepancies should be avoided since they may lead to wrong design decisions.
A considerable amount of research has been devoted to the development of numerical techniques which attempt to overcome the mentioned limitation [2].These methods traditionally assume a solution that is separated in time and space, involving three ingredients [3]: (1) spatial discretization, (2) enforcement of the unilateral contact constraint, and (3) time integration.
In structural dynamics, the most popular technique for spatial discretization is the Finite Element Method (FEM).The displacement u.M; t / of an arbitrary point M at time t is approximated by a finite sum u.M; t / D P i i .M /u i .t/where each term of the sum is separated in space and time: u i .t/ is the participation of the shape function i .M / [4].The initial local governing equation commonly given as a partial differential equation is thus transformed into a set of coupled ordinary differential equations (ODE) which are solely time-dependent.For wave propagation problems, this form yields numerical difficulties because information propagates at an infinite velocity through the considered domain [5].
Unilateral contact conditions are mathematically expressed as a complementarity relation between the contact force and the clearance (gap) separating mechanical bodies [6] which can equivalently be seen as a set-valued function [7].This mathematical framework drastically complicates the development of adequate and robust solution methods.The most popular approaches enforcing unilateral contact conditions are the penalty method [8], the Lagrange multipliers method [9], and non-smooth techniques [10].
The semi-discretized problem in space in the form of ODEs with unilateral contact conditions turns out to be ill-posed [7].In order to recover uniqueness of the solution, an impact law is incorporated to close the system of equations.This approach is mainly used by the rigid-body dynamics community [11].Researchers usually consider the Newton impact law which relates the velocity after the impact v C with the velocity before impact v written as v C D ev , where e is a non-negative parameter called the coefficient of restitution [10].For a non-zero coefficient, non-physical oscillations might arise in the displacement of the contacting surface [12].Choosing e D 0 leads to a zero-velocity phase of the contacting surface but results in a non-acceptable loss of kinetic energy [13].
Time integration of the semi-discretized problem can be performed in various ways which vary according to the approach considered for enforcing the constraints.The idea of numerical time integration is to compute the state of the system advancing from a set of initial conditions in a step-wise fashion through time, going to a future state at time t i C1 from a present state at time t i and separated by a time-step t [13].Classical time integration schemes in elastodynamics such as Newmark and HHT [3] can be implicit or explicit procedures.In this framework, the formulations which make no explicit use of a restitution coefficient e essentially hide it within the procedure handling the contact constraints.It is known that the penalty approach possesses significant numerical limitations since the resulting ODE are known to be stiff and the time-step has to be very small [2,13].The integration in time of penalized models can instead be performed by solvers devoted to stiff ordinary differential equations [14].Newmark and HHT schemes can also be adapted to incorporate Lagrange multipliers [2].The introduction of unilateral contact constraints in an implicit scheme may produce divergent solutions due to the additional high frequency content [15].This can be overcome by adding high frequency numerical dissipation, as proposed in [16], however doing this causes conservative systems to dissipate energy [15].Explicit schemes require small time-steps to ensure stability and an associated high computational cost for large-scale simulations [4].
Time integration of non-smooth discretized formulation target systems where the contact constraint is treated by a nonsmooth approach, and are classified into: event driven methods and time-stepping techniques [17].The event-driven methods are formally based on an accurate contact event detection mechanism and the variable time-step is adapted such that the end of the step coincides with such event.They accurately predict when a contact event occurs; nevertheless, they become inefficient when frequent transitions arise in a short period of time [13].Most of the implementations of these methods have dealt with discrete systems with a few degrees-of-freedom [18].On the other hand, in time-stepping techniques, similar to the Newmark and HHT families mentioned above, detection of events is considered on the same step as the rest of the integration.They have been proven to be convergent and robust but a small time-step is often required [18].Two main approaches have been proposed: (1) the Schatzman-Paoli scheme is based on a position formulation of the unilateral contact constraint [19] and (2) the Moreau-Jean scheme relies on a velocity formulation of the unilateral contact constraint renamed the velocity Signorini condition [6].Recently, it was proposed to enforce contact conditions at the position and velocity levels simultaneously [20].This leads to minimal energy dissipation but non-physical numerical oscillations appear in the simulations.
In the present study, a technique called the Wave Finite Element Method (WFEM) [21] is implemented in order to numerically solve a one-dimensional linear elastodynamics problem with unilateral constraints 1 .Within the framework of continuum mechanics, the WFEM discretizes the system into nodes and elements of finite length x.Then, contrary to traditional finite element techniques, the corresponding impulse-momentum principle is enforced iteratively at each node and element in order to establish a set of time-dependent algebraic relations that provide the mechanical state (i.e.displacements, velocities and forces) of all the elements and nodes.An elastic wave propagating at a finite speed c is reproduced through iterations in time advancing by a known prescribed time interval t D x=c.
A one-dimensional benchmark problem offering a known analytical solution is investigated.The WFEM is compared to (1) an explicit time-stepping technique combining the FEM with forward Lagrange multipliers [9] and (2) a non-smooth approach implemented in the software Siconos [13].Attention is paid to the Gibb's phenomenon generated during and after contact occurrences together with the time evolution of the total energy of the system.
First, the benchmark problem is described.In a second section, the theoretical background of the WFEM is detailed but limited to the equations for the one-dimensional case only, and the appropriate enforcement of the contact conditions in the WFEM is described.Lastly, the results of the benchmark problem simulations are presented and discussed.

Benchmark Problem
In order to study the properties of the WFEM, a one-dimensional elastodynamics benchmark including unilateral constraints and exhibiting a closed-form solution is used in this work.It was first proposed in [3] to explore the numerical properties of various time integration schemes.
The problem consists of a homogeneous elastic rod of length L and constant cross-sectional area A, as shown in Fig. 1, bouncing against a rigid foundation.Its initial displacement at time t 0 D 0 is u 0 .x/and its velocity, v 0 .x/where x 2 OE0 I L is the coordinate of a point of the rod.It is subjected to a constant external body force q 0 which can be seen as "horizontal" gravity here.The contacting end x D L of the rod is initially separated from the rigid foundation by a distance g 0 .Its mass density is and E stands for Young's modulus.An elastic wave will thus propagate at velocity c D p E= .The unknown displacement at position x in the rod and at time t is denoted by u.x; t/; quantity r.t / is the contact force.The gap function is defined as g.t / D g 0 u.L; t /.The formulation is then: where .R / is a second derivative in time while ./ ;xx is a second derivative in space.
Boundary conditions at x D 0: Complementarity conditions at x D L: Initial conditions: This problem has a unique solution [3] and the variation of the energy is equal to the work of the gravity force q 0 , 1 2 The rod is initially at rest, that is u 0 .x/D v 0 .x/D 0, 8x 2 OE0 I L. By choosing a suitable set of parameters, the motion of the rod is periodic in time and known in a closed form.

The Wave Finite Element Method
This section introduces the background of the one-dimensional WFEM with emphasis on the appropriate enforcement of unilateral contact conditions.The solution algorithm is also provided.For a more detailed description of the method, the reader is referred to [21,23].

Formulation
The formulation targets longitudinal elastic wave propagation in a rod with the assumption of a known and constant wave velocity c.For t > 0 in the most general case, the rod is subjected to time dependent boundary forces F ḃ .t/ together with the distributed load q 0 .The superscripts signs (C or ) refer to the location of the nodal force on the boundary, "negative" on the left side and "positive" on the right side, as depicted in Fig. 1.The rod Homogeneous continuous rod of interest is discretized into n elements of equal length x D L=n.For a homogeneous rod, the time t D x=c of wave propagation along each element is space-independent.Quantities referring to the state of the elements will be labeled with a subscript j D 1; 2; : : : ; n.In addition, quantities referring to the two nodes associated to this element will be labeled with the same subscript j together with a superscript that defines the location of the node: "negative" sign on the left side of element j and "positive" sign on the right side.Lastly, external quantities acting on node j will be labeled with a subscript External equivalent nodal forces j D 1; 2; : : : ; n C 1 and an "asterisk" (*) as superscript.The distributed load q 0 is replaced by equivalent nodal forces F j D q 0 x acting on nodes j D 1; : : : ; n C 1, see Fig. 2. On the boundary nodes j D 1 and j D n C 1, nodal forces induced by the distributed load take the following form: and the total nodal forces acting on the boundary of the discretized rod can be written as: In Fig. 3, the discretized rod is depicted along with the boundary forces in Eq. (7).To perform numerical time integration, discrete The boundary nodal forces in Eq. ( 7) are approximated by step-wise values which are assumed to remain constant during a time-step t, and are denoted F 1 .ti / and F C n .ti /.Furthermore, it is assumed that every internal nodal quantity for j D 2; 3; : : : ; n (i.e.displacement, velocity and stress) remain constant during t, and are denoted v j .ti / and j .ti /.Only the states (i.e.velocity and stress) of the elements may vary during a time interval.Two adjacent elements j 1 and j are illustrated in Fig. 4.They are virtually separated for illustration purposes only.The law of momentum conservation at the connecting node j yields: where the pairs .j 1 .ti /; j .ti // and .vj 1 .ti /; v j .ti // are the stresses and velocities of the two elements at time t i .Moreover, 3 inner nodal equilibrium and continuity must be satisfied, that is: where j D F j =A.Substituting Eq. ( 8) into Eq.( 9) leads to a relationship for the stresses and velocities at the internal nodes, for j D 2; 3; : : : ; n: From Eq. ( 10), the state of each element at t iC1 is calculated as follows: The values computed by Eq. ( 11) are employed to determine the nodal quantities at time t iC1 denoted by C j 1 .ti C1 / and v C j 1 .tiC1 / using Eq.(10).Finally, Eq. ( 8) is applied to determine the boundary velocities by taking into account the states of elements j D 1 and j D n.For the given boundary stresses Equations ( 9) to ( 12) establish a set of simple algebraic relations whose solution for each element and node, marching through time with prescribed time-step t, approximates a longitudinal elastic wave travelling in the rod.

Floating boundary conditions
Unilateral contact conditions are such that the contacting node shall either stick to the rigid foundation or remain free.The switch between the two configurations might take place during any time-step t.This contradicts the assumption of the WFEM which says that internal and boundary nodal quantities must remain constant during such time intervals in order to ensure energy and momentum conservation.This can be seen as the compensation of ensuring the exact wave velocity.To enforce contact constraints in the WFEM equations, Shorr developed a concept named Floating Boundary Conditions (FBC) [21].This concept relies on the prediction of the position of the contacting node.If the unilateral constraint of impenetrability is violated by the node during t, a temporary change is performed to the position of the rigid foundation into a local admissible position in order to avoid penetrating condition during timestep t. / > 0, "Situation 1" occurs and the computation procedure continues with a free node n C 1, as depicted in Fig. 6a.On the other hand, if g.t i C1 / < 0, a penetration is detected during the ongoing time-step: this is the "Situation 2" illustrated in Fig. 6b and a temporary change to the position of the rigid foundation needs to be performed into a local admissible position.To this end, a temporary gap g 0 .ti / is calculated if contact is assumed to arise at t i or a temporary gap g 0 .tiC1 / if contact is assumed to occur at t i C1 .Both cases are expressed as: where g 0 .ti / denotes the distance that the rigid foundation has to undergo to be in a temporary admissible position.This distance is calculated by comparing jg.t i C1 /j with g.t i /: 1.If jg.t i C1 /j g.t i /, the rigid foundation virtually travels the distance jg.t iC1 /j such that g 0 .ti / D jg.t i C1 /j.Then, a virtual gap g 0 .tiC1 / is calculated with Eq. ( 15) and contact will occur at the subsequent time t iC1 , see Fig. 7a.

Solution Algorithm
The proposed solution algorithm 1 is a time marching procedure where Eq. ( 9) to (12) are solved iteratively at each node and element, starting from node and element j D 1.To simulate an elastic wave propagating with finite speed c, these equations are iterated in time advancing by a prescribed time-step t D x=c.Node n C 1 is the contacting node on which the FBC principle is implemented.The procedure at time t i can be divided in the following steps: 1. Calculation of the velocity v 1 .ti / through Eq. ( 12) and calculation of the velocity v 1 .ti C1 / and stress 1 .ti C1 / of element j D 1 for the subsequent time t iC1 using Eq. ( 11).
2. A recurrence procedure for j D 2; 3; : : : ; n is performed to calculate all nodal quantities and states of the elements, employing Eq. ( 9) to (11).end -end of floating boundary conditionsend Output: displacements, stresses, velocities, contact forces contact is activated, the position of the rigid foundation is temporarily changed and the velocity of the contacting node is set to v C n .ti / D 0. Otherwise, the calculation proceeds on node n C 1 with a free boundary condition.

Results
Dimensionless parameters identical to [3] are considered and listed in Tab. 1.The rod is expected to have a periodic motion in time with periodic bounces against the foundation.The WFEM is compared to (1) an explicit time-marching technique combining the FEM with forward Lagrange multipliers [9] and (2) a nonsmooth approach implemented in the package Siconos [13]  The time-step of the FEM and Siconos satisfies the CFL condition [5] t x=c.In Siconos, two coefficients of restitution are considered: e D 0 (inelastic) and e D 1 (perfectly elastic).Two computations were performed to evaluate the influence of the number of elements: (1) n D 100 and (2) n D 500 with corresponding time-steps listed in Tab. 2. Contact node displace- ments are shown in Figs. 8 and 12, respectively.The spurious oscillations occurring at the first impact of the rod are depicted in Figs. 9 and 13.Energy levels normalized with respect to the exact constant energy are depicted in Fig. 14.

Coarse spatial discretization
The contact node displacement and energy behavior for n D 100 is depicted, respectively, in Fig. 8 and Fig. 14.FEM and Siconos (e D 0) present an unsatisfactory energy behavior due to dissipation when the contact node impacts the foundation; such behavior is aggravated with subsequent impacts.This is caused by the energy dissipation properties of the time integration schemes [13,3].In addition, FEM does not exhibit Gibb's phenomenon in the displacement approximation, as depicted in Fig. 9, which is due to the exact enforcement of contact constraints in displacement [9].Further, for Siconos (e D 0) spurious oscillations are not present in the displacement approximation since the node remains stuck to the foundation because of v C D 0; however, the complementarity impenetrability condition is violated.Siconos (e D 1) represents a perfectly elastic collision and the total energy is conserved, as shown in Fig. 14.Nevertheless, a coefficient e D 1 physically assumes that the velocity of the contacting node after the impact v C is the exact opposite of the velocity v before the impact: accordingly, the contacting node will hardly be in a constrained position on the foundation.This is reflected in the approximation of the contact node displacement which features non-physical spurious oscillations after the first impact as depicted in Fig. 9.The approximation of the contact node displacement deteriorates after several bounces of the rod.Contrary to the other methods, the WFEM does not dissipate energy during unilateral contact occurrences.Furthermore, the approximation of the contact node displacement does not present numerical oscillations during the collision with the foundation, and the unilateral constraint is exactly enforced as displayed in Fig. 9 where the contact node velocity goes from a non-zero velocity v to exactly a zero velocity during the contact phase, before it retrieves a non-zero post-contact velocity v C .In other words, it is numerically shown that the implemented formulation of the WFEM does not "hide" neither necessitate any impact law to be well-posed.The velocity and contact force are depicted in where the solution has Gibb's phenomenon in the free flying phase without any oscillations during the contact phase.Instead, the WFEM is capable of capturing discontinuities in the velocity field with no undesirable oscillations: this is true for all times.

Fine spatial discretization
By increasing the number of elements to n D 500, an improved displacement approximation is obtained for the FEM and Siconos (e D 0) solutions, see Fig. 12. Their energy behavior still presents dissipation as multiple impacts occur in time, as observed in Fig 14 .Siconos (e D 1) provides a displacement approximation similar to its counterpart but energy is not dissipated.In Fig. 13, Siconos (e D 1) still introduces non-physical oscillations on the displacement of the contact node: although not as pronounced as for n D 100, this badly affects the dynamics of the whole system.For e D 0, the contact position is erro-  neously predicted which causes the contacting node to violate the complementarity conditions.For the WFEM, the time-step decreases as spatial discretization gets finer and again it is observed how the WFEM is capable of properly approximating displacements with higher accuracy than the other methods, and without energy dissipation as depicted in Fig. 14.This comes together with the fact that the unilateral constraint is exactly enforced without numerical oscillations.

Computation time
The computation time for each method is depicted in Fig. 15.All simulations were performed with parameters listed in Tab. 1 on a desktop computer with processor Intel Core i7-2600, CPU@3.4GHz and 16 GiB of RAM memory.It is clear that the WFEM outperforms its competitors.From the results depicted in Fig. 14, it is evident that the approaches based on FEM discretization present numerical energy dissipation that clearly affects the approximate solution and provides an erroneous prediction of the rod dynamics.Increasing the number of elements reduces the numerical dissipation, at the cost of a higher computational effort.

Conclusions
This study focused on the application of the WFEM to a onedimensional elastodynamics problem subjected to unilateral constraints.The WFEM was compared to (1) an explicit timestepping technique combining the FEM with forward Lagrange multipliers and (2) a non-smooth time-stepping approach implemented in the software "Siconos".Attention was paid to the Gibb's phenomenon generated during and after contact occurrences together with the time evolution of the total energy of the system.Additionally, the computation time of each method is included.
From the results, it is clear that the WFEM properly captures unilateral contact induced waves propagating at a finite speed with discontinuous deformation, stress and velocity wave fronts.This approach provides an accurate approximation of the solution without presenting any spurious oscillations.Energy is not dissipated over time as opposed to what is observed with the FEM and Siconos solvers.These advantageous numerical properties are not affected as time increases and get improved with a finer spatial discretization.This numerically shows that FEM-based numerical approximations assuming a solution in the form of a finite sum u.M; t / D P i i .M /u i .t/separated in space and time are inappropriate to accurately capture travelling waves, which is crucial in elastodynamics involving unilateral contact constraints.Additionally, the WFEM does not require any impact law to retrieve the exact solution as opposed to all FEM formulations.This is a spectacular outcome of the present study since these impact laws are questionable in the context of continuum mechanics.
These promising results should be complemented by further WFEM investigations in a multi-dimensional framework.This is possible by enforcing the impulse-momentum principle at the connecting nodes of adjacent elements, assuming small strain, and small element rotation.A discussion on two-dimensional problems can be found in [21], however it is limited to rectangular domains and elastic materials.Further studies should also focus on the implementation of the WFEM in complex geometries, three-dimensional domains, and large deformations problems.

Figure 3 .
Figure 3. Boundary forces acting on the rod

Figure 4 .
Figure 4. Interaction between elements j 1 j time t i

Figure 10 .
Figure 10.Contact node velocity for n D 100: Exact Sol.[ ], WFEM [ ], FEM [ ], Siconos (e D 0) [ ] and Siconos (e D 1) [ ] Figs. 10 and 11.Siconos (e D 1) presents unphysical oscillations during the contact phase and this behavior is aggravated with time.FEM and Siconos (e D 0) predict a similar response,where the solution has Gibb's phenomenon in the free flying phase without any oscillations during the contact phase.Instead, the WFEM is capable of capturing discontinuities in the velocity field with no undesirable oscillations: this is true for all times.

Figure 14 .Figure 15 .
Figure 14.Normalized energy for n D 100: FEM [ ], Siconos (e D 0) [ ]; n D 500: FEM [ ], Siconos (e D 0) [ ]; Siconos (e D 1) [ ]; WFEM [ ] Figure 5. approaching the rigid foundation at t i Consider element n of the rod impacting a rigid foundation as displayed in Fig.5.The gap between the rod and the foundation at instant t i is g.t i / and the gap at the subsequent instant t iC1 is g.t iC1 /.During the corresponding time-step t, four situations may arise: 1. Open gap remains, g.t i / > 0 and g.t i C1 / > 0 2. Open gap g.t i / > 0 changes to penetration g.t iC1 / < 0 3. Permanent contact g.t i / D g.t i C1 / D 0 4. Initial contact g.t i / D 0 changes to open gap g.t iC1 / > 0 An open gap occurs at t i , such that g.t i / > 0, as illustrated in Fig. 5. Accordingly, contact node n C 1 is free, that is C n .ti / D nC1 .Employing Eq. (12) to compute the contact node velocity v C n .ti / yields: v C n .ti / D v n .ti / C .nC1 n .ti //=.c/ (13) The gap g.t i C1 / is then calculated through: g.t iC1 / D g.t i / C g.t i / with g.t i / D v C n .ti /t (14) where the two states g.t iC1 / > 0 or g.t i C1 / < 0 might arise.If g.t iC1 2. If jg.t i C1 /j > g.t i /, the rigid foundation virtually travels the distance g.t i /, hence g 0 .ti / D g.t i /.Consequently, contact arises at time t i and the virtual gap g 0 .ti / is calculated with Eq. (15), see Fig. 7b.In this case, calculation at time t i must be performed considering that the contacting node is not moving due to the contact condition, such that v C n .ti / D 0. The position of the rigid foundation can be subsequently restored when an open gap arises again; in such a case, g 0 .ti / D 0. g 0 .ti / D jg.t iC1 /j Temporary change to the position of the rigid foundation Situations 3 or 4 A contact occurs at t i such that g.t i / D 0.
g 0 .ti / D g.t i / Change (b) Contact occurring at t i Figure 7. n .ti / < 0, the contact condition is true and the calculation can continue with v C n .ti / D 0. On the other hand, if C n .ti / > 0 a recalculation of the step considering a free boundary condition is performed.
Simulation Procedure Input: number of elements n, total number of steps N , initial gap g 0 , boundary conditions at node 1 and n C 1, initial conditions, external load q 0 , wave velocity c for i D 0 W N [Time Loop] do Discrete time instant, t i D it ; Stress and velocity at node 1; 3. Assumed contact g.t i / D 0. The contact stress is calculated to determine if such an assumption is satisfied.If Algorithm 1:

Table 2 .
Time-step for simulations