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

Carlos Yoong, Anders Thorin, Mathias Legrand

► To cite this version:

Carlos Yoong, Anders Thorin, Mathias Legrand. The Wave Finite Element Method applied to a one-dimensional linear elastodynamic problem with unilateral constraints. The ASME 2015 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, Aug 2015, Boston, United States. 10.1115/DETC2015-46919. hal-01194922

HAL Id: hal-01194922

https://hal.science/hal-01194922

Submitted on 26 Nov 2016

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

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

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

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# The Wave Finite Element Method applied to a one-dimensional linear elastodynamic problem with unilateral constraints

Carlos Yoong*, Anders Thorin, Mathias Legrand

**Abstract**

The Wave Finite Element Method (WFEM) is implemented to accurately capture travelling waves propagating at a finite 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.

**Keywords**

wave finite element method — unilateral contact constraints — elastic wave propagation — non-smooth dynamics

Department of Mechanical Engineering, McGill University, 817 Sherbrooke Street West, Montréal, Québec, Canada H3A 0C3

* Corresponding author: carlos.yoong@mail.mcgill.ca

## 1. 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) = \sum_i \phi_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 $\phi_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^+$ with the velocity before impact $v^-$ written as $v^+ = -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 = 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+1}$ from a present state at time $t_i$ and separated by a time-step $\Delta t$ [13]. Classical time
Nomenclature

<table><tr><td>ρ</td><td>mass density of the rod</td></tr><tr><td>A</td><td>cross-sectional area of the rod</td></tr><tr><td>c</td><td>wave velocity</td></tr><tr><td>E</td><td>Young's Modulus of the rod</td></tr><tr><td>L</td><td>length of the rod</td></tr><tr><td>g</td><td>time-dependent clearance function</td></tr><tr><td>g<sub>0</sub></td><td>initial clearance between bar tip and rigid foundation</td></tr><tr><td>q<sub>0</sub></td><td>external body force</td></tr><tr><td>u, ˙u, ˙ü</td><td>space-time displacement, velocity, and acceleration fields</td></tr><tr><td>u<sub>0</sub>, v<sub>0</sub></td><td>space initial displacement and velocity fields</td></tr><tr><td>Δt</td><td>time-step</td></tr><tr><td>Δx</td><td>element length</td></tr><tr><td>N</td><td>total number of steps (in time)</td></tr></table>

<table><tr><td>n</td><td>number of elements (in space)</td></tr><tr><td>t<sub>i</sub></td><td>time instant</td></tr><tr><td>F<sub>h</sub><sup>+</sup>, F<sub>b</sub><sup>-</sup></td><td>boundary force acting on right/left side of the rod</td></tr><tr><td>v<sub>j</sub><sup>+</sup>, v<sub>j</sub><sup>-</sup></td><td>nodal velocity of the right/left side of element j</td></tr><tr><td>σ<sub>j</sub><sup>+</sup>, σ<sub>j</sub><sup>-</sup></td><td>nodal stress acting on right/left side of element j</td></tr><tr><td>F<sub>j</sub><sup>*</sup>, σ<sub>j</sub><sup>*</sup></td><td>external nodal force/stress acting on node j</td></tr><tr><td>σ<sub>j</sub></td><td>stress acting on element j</td></tr><tr><td>v<sub>j</sub></td><td>velocity of element j</td></tr><tr><td>F<sub>1</sub><sup>-</sup></td><td>boundary force acting on left side of element 1</td></tr><tr><td>v<sub>1</sub><sup>-</sup></td><td>boundary velocity of the left side of element 1</td></tr><tr><td>F<sub>n</sub><sup>+</sup></td><td>boundary force acting on right side of element n</td></tr><tr><td>v<sub>n</sub><sup>+</sup></td><td>boundary velocity of the right side of element n</td></tr><tr><td>g'</td><td>temporary gap</td></tr><tr><td>Δg'</td><td>temporary distance to rigid foundation</td></tr></table>

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 non-smooth 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<fnref n="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 = Δ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.

## 2. 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 = 0$ is $u_0(x)$ and its velocity, $v_0(x)$ where $x \in [0; 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 = 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 = \sqrt{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

<footnote n="1">¹It should be noted that a different method under the same name Wave Finite Element Method was proposed by Mace and collaborators [22]. The WFEM considered in the present work has a completely different ground.</footnote>
contact force. The gap function is defined as $g(t) = g_0 - u(L, t)$. The formulation is then:

* Local equation:

$$ \rho\ddot{u} - Eu_{,xx} = \rho q_0, \quad \forall x \in ]0; L[, \quad \forall t > 0 \quad (1) $$

where (•) is a second derivative in time while (•)$_{xx}$ is a second derivative in space.

* Boundary conditions at $x = 0$:

$$ u_{,x}(0, t) = 0, \quad \forall t > 0 \quad (2) $$

* Complementarity conditions at $x = L$:

$$ g(t) \ge 0, \quad r(t) \ge 0, \quad r(t)g(t) = 0, \quad \forall t > 0 \quad (3) $$

* Initial conditions:

$$ u(x, 0) = u_0(x), \quad \dot{u}(x, 0) = v_0(x), \quad \forall x \in ]0; L[ \quad (4) $$

This problem has a unique solution [3] and the variation of the energy is equal to the work of the gravity force $q_0$.

$$ \frac{1}{2} \frac{d}{dt} \left( \int_{0}^{L} (\rho \dot{u}^2 + Eu_{,x}^2) dx \right) = \int_{0}^{L} -\rho q_0 \dot{u} dx, \quad \forall t > 0 \quad (5) $$

The rod is initially at rest, that is $u_0(x) = v_0(x) = 0, \forall x \in [0; L]$. By choosing a suitable set of parameters, the motion of the rod is periodic in time and known in a closed form.

### 3. 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].

#### 3.1 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_b^\pm(t)$ together with the distributed load $q_0$. The superscripts $+$ 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

<figure><img src="image_2.png" /><figcaption>Figure 1. Homogeneous continuous rod of interest</figcaption></figure>

is discretized into $n$ elements of equal length $\Delta x = L/n$. For a homogeneous rod, the time $\Delta t = \Delta 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 = 1, 2, \dots, 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

<figure><img src="image_3.png" /><figcaption>Figure 2. External equivalent nodal forces</figcaption></figure>

$j = 1, 2, \dots, n + 1$ and an “asterisk” (*) as superscript. The distributed load $q_0$ is replaced by equivalent nodal forces $F_j^* = q_0 \Delta x$ acting on nodes $j = 1, \dots, n + 1$, see Fig. 2. On the boundary nodes $j = 1$ and $j = n + 1$, nodal forces induced by the distributed load take the following form:

$$ F_1^* = F_{n+1}^* = \frac{1}{2} q_0 \Delta x \quad (6) $$

and the total nodal forces acting on the boundary of the discretized rod can be written as:

$$ F_1^-(t) = F_b^-(t) + F_1^* \quad \text{and} \quad F_n^+(t) = F_b^+(t) + F_{n+1}^* \quad (7) $$

In Fig. 3, the discretized rod is depicted along with the boundary forces in Eq. (7). To perform numerical time integration, discrete

<figure><img src="image_4.png" /><figcaption>Figure 3. Boundary forces acting on the rod</figcaption></figure>

time instants $t_i = i \Delta t, i = 0, 1, \dots, N$ are defined, $N$ being the total number of time-steps. The time-step is assumed constant such that $\Delta t = \Delta x/c$ holds true where $c$ is the wave velocity.

The boundary nodal forces in Eq. (7) are approximated by step-wise values which are assumed to remain constant during a time-step $\Delta t$, and are denoted $F_1^-(t_i)$ and $F_n^+(t_i)$. Furthermore, it is assumed that every internal nodal quantity for $j = 2, 3, \dots, n$ (i.e. displacement, velocity and stress) remain constant during $\Delta t$, and are denoted $v_j^\pm(t_i)$ and $\sigma_j^\pm(t_i)$. Only the states (i.e. velocity and stress) of the elements may vary during a time interval.

<figure><img src="image_5.png" /><figcaption>Figure 4. Interaction between elements $j-1$ and $j$ at time $t_i$</figcaption></figure>

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:

$$ \rho c(v_{j-1}^{+}(t_i) - v_{j-1}(t_i)) = \sigma_{j-1}^{+}(t_i) - \sigma_{j-1}(t_i) \\ \rho c(v_j^{-}(t_i) - v_j(t_i)) = -(\sigma_j^{-}(t_i) - \sigma_j(t_i)) \quad (8) $$

where the pairs $(\sigma_{j-1}(t_i), \sigma_j(t_i))$ and $(v_{j-1}(t_i), v_j(t_i))$ are the stresses and velocities of the two elements at time $t_i$. Moreover,
inner nodal equilibrium and continuity must be satisfied, that is:

$$ 
\begin{aligned}
\sigma_j^* + \sigma_j^-(t_i) - \sigma_{j-1}^+(t_i) &= 0 \quad &\text{Equilibrium} \\
v_j^-(t_i) &= v_{j-1}^+(t_i) \quad &\text{Continuity}
\end{aligned}
\qquad (9) $$

where $\sigma_j^* = F_j^*/A$. Substituting Eq. (8) into Eq. (9) leads to a relationship for the stresses and velocities at the internal nodes, for $j = 2, 3, \dots, n$:

$$ 
\begin{aligned}
\sigma_{j-1}^{+}(t_i) &= \frac{1}{2}(\sigma_j(t_i) + \sigma_{j-1}(t_i) + \sigma_j^* + \rho c(v_j(t_i) - v_{j-1}(t_i))) \\
v_{j-1}^{+}(t_i) &= v_j(t_i) + (\sigma_{j-1}^{+}(t_i) - \sigma_j(t_i))/(\rho c)
\end{aligned}
\qquad (10) $$

From Eq. (10), the state of each element at $t_{i+1}$ is calculated as follows:

$$ 
\begin{aligned}
\sigma_{j-1}(t_{i+1}) &= \sigma_{j-1}^{-}(t_i) + \sigma_{j-1}^{+}(t_i) - \sigma_{j-1}(t_i) \\
v_{j-1}(t_{i+1}) &= v_{j-1}^{-}(t_i) + v_{j-1}^{+}(t_i) - v_{j-1}(t_i)
\end{aligned}
\qquad (11) $$

The values computed by Eq. (11) are employed to determine the nodal quantities at time $t_{i+1}$ denoted by $\sigma_{j-1}^+(t_{i+1})$ and $v_{j-1}^+(t_{i+1})$ using Eq. (10). Finally, Eq. (8) is applied to determine the boundary velocities by taking into account the states of elements $j = 1$ and $j = n$. For the given boundary stresses $\sigma_1^-(t_i) = F_1^-(t_i)/A$ and $\sigma_n^+(t_i) = F_n^+(t_i)/A$, the velocities become:

$$ 
\begin{aligned}
v_1^-(t_i) &= v_1(t_i) - (\sigma_1^-(t_i) - \sigma_1(t_i))/(\rho c), \\
v_n^+(t_i) &= v_n(t_i) + (\sigma_n^+(t_i) - \sigma_n(t_i))/(\rho c)
\end{aligned}
\qquad (12) $$

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 $\Delta t$, approximates a longitudinal elastic wave travelling in the rod.

## 3.2 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 $\Delta 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 $\Delta 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 time-step $\Delta t$.

<figure><img src="image_6.png" /><figcaption>Figure 5. Rod approaching the rigid foundation at $t_i$</figcaption></figure>

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_{i+1}$ is $g(t_{i+1})$. During the corresponding time-step $\Delta t$, four situations may arise:

1. Open gap remains, $g(t_i) > 0$ and $g(t_{i+1}) > 0$

2. Open gap $g(t_i) > 0$ changes to penetration $g(t_{i+1}) < 0$

3. Permanent contact $g(t_i) = g(t_{i+1}) = 0$

4. Initial contact $g(t_i) = 0$ changes to open gap $g(t_{i+1}) > 0$

<figure><img src="image_7.png" /><figcaption>Figure 6. Possible contact configurations during time step $\Delta t$</figcaption></figure>

**Situations 1 or 2** An open gap occurs at $t_i$, such that $g(t_i) > 0$, as illustrated in Fig. 5. Accordingly, contact node $n+1$ is free, that is $\sigma_n^+(t_i) = \sigma_{n+1}^*$. Employing Eq. (12) to compute the contact node velocity $v_n^+(t_i)$ yields:

$$ v_n^+(t_i) = v_n(t_i) + (\sigma_{n+1}^* - \sigma_n(t_i))/(\rho c) \qquad (13) $$

The gap $g(t_{i+1})$ is then calculated through:

$$ g(t_{i+1}) = g(t_i) + \Delta g(t_i) \quad \text{with} \quad \Delta g(t_i) = -v_n^+(t_i)\Delta t \qquad (14) $$

where the two states $g(t_{i+1}) > 0$ or $g(t_{i+1}) < 0$ might arise. If $g(t_{i+1}) > 0$, “Situation 1” occurs and the computation procedure continues with a free node $n+1$, as depicted in Fig. 6a. On the other hand, if $g(t_{i+1}) < 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'(t_i)$ is calculated if contact is assumed to arise at $t_i$ or a temporary gap $g'(t_{i+1})$ if contact is assumed to occur at $t_{i+1}$. Both cases are expressed as:

$$ 
\begin{aligned}
g'(t_i) &= g(t_i) + \Delta g'(t_i) \quad \text{at time } t_i \\
g'(t_{i+1}) &= g(t_i) + \Delta g'(t_i) \quad \text{at time } t_{i+1}
\end{aligned}
\qquad (15) $$

where $\Delta g'(t_i)$ denotes the distance that the rigid foundation has to undergo to be in a temporary admissible position. This distance is calculated by comparing $|g(t_{i+1})|$ with $g(t_i)$:

1. If $|g(t_{i+1})| \le g(t_i)$, the rigid foundation virtually travels the distance $|g(t_{i+1})|$ such that $\Delta g'(t_i) = |g(t_{i+1})|$. Then, a virtual gap $g'(t_{i+1})$ is calculated with Eq. (15) and contact will occur at the subsequent time $t_{i+1}$, see Fig. 7a.
2. If $|g(t_{i+1})| > g(t_i)$, the rigid foundation virtually travels the distance $-g(t_i)$, hence $\Delta g'(t_i) = -g(t_i)$. Consequently, contact arises at time $t_i$ and the virtual gap $g'(t_i)$ 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_n^+(t_i) = 0$.

The position of the rigid foundation can be subsequently restored when an open gap arises again; in such a case, $\Delta g'(t_i) = 0$.

<figure><img src="image_8.png" /><figcaption>(a) Contact occurring at $t_{i+1}$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Contact occurring at $t_i$</figcaption></figure>

<caption>**Figure 7.** Temporary change to the position of the rigid foundation</caption>

**Situations 3 or 4** A contact occurs at $t_i$ such that $g(t_i) = 0$. Then, assuming that contact remains during the time-step, the contact stress $\sigma_n^+(t_i)$ is calculated with $v_n^+(t_i) = 0$ employing Eq. (12):

$$ \sigma_n^+(t_i) = \sigma_n(t_i) - \rho c v_n(t_i) + \sigma_{n+1}^* \quad (16) $$

where both $\sigma_n^+(t_i) < 0$ (compression) and $\sigma_n^+(t_i) > 0$ (tension) can occur. If $\sigma_n^+(t_i) < 0$, the contact condition is true and the calculation can continue with $v_n^+(t_i) = 0$. On the other hand, if $\sigma_n^+(t_i) > 0$ a recalculation of the step considering a free boundary condition is performed.

### 3.3 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 = 1$. To simulate an elastic wave propagating with finite speed $c$, these equations are iterated in time advancing by a prescribed time-step $\Delta t = \Delta x/c$. Node $n + 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^-(t_i)$ through Eq. (12) and calculation of the velocity $v_1(t_{i+1})$ and stress $\sigma_1(t_{i+1})$ of element $j = 1$ for the subsequent time $t_{i+1}$ using Eq. (11).

2. A recurrence procedure for $j = 2, 3, \dots, n$ is performed to calculate all nodal quantities and states of the elements, employing Eq. (9) to (11).

3. Assumed contact $g(t_i) = 0$. The contact stress is calculated to determine if such an assumption is satisfied. If

**Algorithm 1: Simulation Procedure**

**Input:** number of elements $n$, total number of steps $N$, initial gap $g_0$, boundary conditions at node 1 and $n + 1$, initial conditions, external load $q_0$, wave velocity $c$

**for** $i = 0 : N$ **[Time Loop]** **do**

Discrete time instant, $t_i = i \Delta t$;

Stress and velocity at node 1;

Stress and velocity at element 1;

**for** $j = 2 : n$ **[Element Loop]** **do**

Stress and velocity at node $j - 1$;

Stress and velocity at element $j - 1$;

Apply continuity and equilibrium;

**end**

— floating boundary conditions —

**if** $g(t_i) = 0$ **then**

Contact occurs at node $n + 1$;

Calculation of contact stress: $\sigma_n^+(t_i)$;

**if** $\sigma_n^+(t_i) > 0$ **then**

Free boundary condition at node $n + 1$;

Calculate gap: $g(t_{i+1})$;

Change rigid foundation position;

**else**

Contact occurs at node $n + 1$;

**end**

**else**

Free boundary condition at node $n + 1$;

Calculate gap: $g(t_{i+1})$;

Change rigid foundation position;

**end**

— end of floating boundary conditions —

**end**

**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_n^+(t_i) = 0$. Otherwise, the calculation proceeds on node $n + 1$ with a free boundary condition.

## 4. 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 non-smooth approach implemented in the package Siconos [13]<fnref n="2" />.

<table><caption>**Table 1.** Dimensionless simulation parameters</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Simulation Time</td><td>20</td></tr><tr><td>Young's Modulus, <i>E</i></td><td>900</td></tr><tr><td>Density, <i>ρ</i></td><td>1</td></tr><tr><td>Rod Length, <i>L</i></td><td>10</td></tr><tr><td>Initial Gap, <i>g</i><sub>0</sub></td><td>5</td></tr><tr><td>External Body Force, <i>q</i><sub>0</sub></td><td>10</td></tr><tr><td>Wave Velocity, <i>c</i> = &radic;<i>E</i>/<i>ρ</i></td><td>30</td></tr></tbody></table>

The time-step of the FEM and Siconos satisfies the CFL condition [5] $\Delta t \le \Delta x/c$. In Siconos, two coefficients of restitution

<footnote n="2">²Moreau's Time Stepping Scheme based on the $\theta$-Method is employed. For the proposed simulations, an implicit version is considered with $\theta = 1/2$.</footnote>
are considered: $e = 0$ (inelastic) and $e = 1$ (perfectly elastic). Two computations were performed to evaluate the influence of the number of elements: (1) $n = 100$ and (2) $n = 500$ with corresponding time-steps listed in Tab. 2. Contact node displacement 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.

<table><caption>**Table 2.** Time-step for simulations</caption><thead><tr><th>Method</th><th>Δt [n = 100]</th><th>Δt [n = 500]</th></tr></thead><tbody><tr><td>FEM</td><td>10<sup>-3</sup></td><td>10<sup>-4</sup></td></tr><tr><td>Siconos</td><td>10<sup>-3</sup></td><td>10<sup>-4</sup></td></tr><tr><td>WFEM</td><td>Δx/c = 3 × 10<sup>-3</sup></td><td>Δx/c = 6 × 10<sup>-4</sup></td></tr></tbody></table>

## 4.1 Coarse spatial discretization

The contact node displacement and energy behavior for $n = 100$ is depicted, respectively, in Fig. 8 and Fig. 14. FEM and Siconos ($e = 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

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Contact node displacement for $n = 100$: Exact Sol. [---], WFEM [—], FEM [—], Siconos ($e = 0$) [—] and Siconos ($e = 1$) [—]</figcaption></figure>

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 = 0$) spurious oscillations are not present in the displacement approximation since the node remains stuck to the foundation because of $v^+ = 0$; however, the complementarity impenetrability condition is violated. Siconos ($e = 1$) represents a perfectly elastic collision and the total energy is conserved, as shown in Fig. 14. Nevertheless, a coefficient $e = 1$ physically assumes that the velocity of the contacting node after the impact $v^+$ 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,

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Gibb's phenomenon for $n = 100$: WFEM [—], FEM [—], Siconos ($e = 0$) [—] and Siconos ($e = 1$) [—]</figcaption></figure>

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^+$. 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

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Contact node velocity for $n = 100$: Exact Sol. [---], WFEM [—], FEM [—], Siconos ($e = 0$) [—] and Siconos ($e = 1$) [—]</figcaption></figure>

Figs. 10 and 11. Siconos ($e = 1$) presents unphysical oscillations during the contact phase and this behavior is aggravated with time. FEM and Siconos ($e = 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.

## 4.2 Fine spatial discretization

By increasing the number of elements to $n = 500$, an improved displacement approximation is obtained for the FEM and Siconos ($e = 0$) solutions, see Fig. 12. Their energy behavior still presents dissipation as multiple impacts occur in time, as observed in Fig. 14. Siconos ($e = 1$) provides a displacement approximation similar to its counterpart but energy is not dissipated.

In Fig. 13, Siconos ($e = 1$) still introduces non-physical oscillations on the displacement of the contact node: although not as pronounced as for $n = 100$, this badly affects the dynamics of the whole system. For $e = 0$, the contact position is erro-
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Contact force for $n = 100$: Exact sol. [---], WFEM [—], FEM [—], Siconos ($e = 0$) [—] and Siconos ($e = 1$) [—]</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Contact node displacement for $n = 500$: Exact sol. [---], WFEM [—], FEM [—], Siconos ($e = 0$) [—] and Siconos ($e = 1$) [—]</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Gibb's phenomenon for $n = 500$: WFEM [—], FEM [—], Siconos ($e = 0$) [—] and Siconos ($e = 1$) [—]</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Normalized energy for $n = 100$: FEM [---], Siconos ($e = 0$) [---]; $n = 500$: FEM [—], Siconos ($e = 0$) [—]; Siconos ($e = 1$) [—]; WFEM [---]</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Computation time: WFEM [—], FEM [—] and Siconos ($e = 0, e = 1$) [—]</figcaption></figure>

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.

## 4.3 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.4 GHz 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.

## 5. Conclusions

This study focused on the application of the WFEM to a one-dimensional elastodynamics problem subjected to unilateral constraints. The WFEM was compared to (1) an explicit time-stepping 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) = \sum_i \phi_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.

## References

[1] S. HAM, K.J. BATHE. "A finite element method enriched for wave propagation problems". *Computers & Structures* 94-95 2012, pp 1-12.
DOI: 10.1016/j.compstruc.2012.01.001.

[2] P. WRIGGERS. *Computational Contact Mechanics*. Springer, 2006.
DOI: 10.1007/978-3-540-32609-0.

[3] D. DOYEN, A. ERN, S. PIPERNO. "Time-integration schemes for the finite element dynamic Signorini problem". *SIAM Journal of Scientific Computing* 33(1) 2011, pp 223-249.
DOI: 10.1137/100791440.
OAI: hal-00440128.

[4] T.J.R. HUGHES. *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*. Prentice Hall, Inc., 1987.

[5] K.J. BATHE. *Finite Element Procedures*. Prentice Hall, 1996.

[6] M. JEAN. "The non-smooth contact dynamics method". *Computer Methods in Applied Mechanics and Engineering* 177(3-4) 1999, pp 235-257.
DOI: 10.1016/S0045-7825(98)00383-1.

[7] C. ECK, J. JARUSEK, M. KRBEC. *Unilateral Contact Problems: Variational Methods and Existence Theorems*. CRC, 2005.

[8] F. ARMERO, E. PETOCZ. "Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems". *Computer Methods in Applied Mechanics and Engineering* 158(3-4) 1998, pp 269-300.
DOI: 10.1016/S0045-7825(97)00256-9.

[9] N.J. CARPENTER, R.L. TAYLOR, M.G. KATONA. "Lagrange constraints for transient finite element surface contact". *International Journal for Numerical Methods in Engineering* 32(1) 1991, pp 103-128.
DOI: 10.1002/nme.1620320107.

[10] B. BROGLIATO. *Nonsmooth Mechanics: Models, Dynamics and Control*. Springer, 1999.
DOI: 10.1007/978-1-4471-0557-2.

[11] J.J. MOREAU. "An introduction to unilateral dynamics". *Novel Approaches in Civil Engineering*. Ed. by M. Frémond, F. Maceri. Vol. 14. Lecture Notes in Applied and Computational Mechanics. Springer, 2004, pp 1-46.
DOI: 10.1007/978-3-540-45287-4_1.

[12] J.J. MOREAU. "Unilateral contact and dry friction in finite freedom dynamics". *Nonsmooth Mechanics and Applications*. Ed. by J.J. Moreau, P.D. Panagiotopoulos. Vol. 302. International Centre for Mechanical Sciences. Springer, 1988, pp 1-82.
DOI: 10.1007/978-3-7091-2624-0_1.

[13] V. ACARY, B. BROGLIATO. *Numerical Methods for Nonsmooth Dynamical Systems*. Springer, 2008.
DOI: 10.1007/978-3-540-75392-6.

[14] M.T. DARVISHI, F. KHANI, A.A. SOLIMAN. "The numerical simulation for stiff systems of ordinary differential equations". *Computers & Mathematics with Applications* 54(7-8) 2007, pp 1055-1063.
DOI: 10.1016/j.camwa.2006.12.072.

[15] R. KRAUSE, M. WALLOTH. "Presentation and comparison of selected algorithms for dynamic contact based on the Newmark scheme". *Applied Numerical Mathematics* 62(10) 2012, pp 1393-1410.
DOI: 10.1016/j.apnum.2012.06.014.

[16] P. DEUFLHARD, R. KRAUSE, S. ERTEL. "A contact-stabilized Newmark method for dynamical contact problems". *International Journal for Numerical Methods in Engineering* 73(9) 2008, pp 1274-1290.
DOI: 10.1002/nme.2119.

[17] L. PAOLI. "Time discretization of vibro-impact". English. *Philosophical Transactions: Mathematical, Physical and Engineering Sciences* 359(1789) 2001, pp 2405-2428.
DOI: 10.1098/rsta.2001.0858.

[18] C. STUDER. *Numerics of Unilateral Contacts and Friction: Modeling and Numerical Time Integration in Non-Smooth Dynamics*. Springer, 2009.
DOI: 10.1007/978-3-642-01100-9.

[19] L. PAOLI, M. SCHATZMAN. "A numerical scheme for impact problems I: The one-dimensional case". *SIAM Journal on Numerical Analysis* 40(2) 2002, pp 702-733.
DOI: 10.1137/S0036142900378728.

[20] O. BRÜLS, V. ACARY, A. CARDONA. "Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized $\alpha$-scheme". *Computer Methods in Applied Mechanics and Engineering* 281 2014, pp 131-161.
DOI: 10.1016/j.cma.2014.07.025.
OAI: hal-01059823.

[21] B.F. SHORR. *The Wave Finite Element Method*. Springer, 2004.
DOI: 10.1007/978-3-540-44579-1.

[22] B.R. MACE, D. DUHAMEL, M.J. BRENNAN, L. HINKE. "Finite element prediction of wave motion in structural waveguides". *The Journal of the Acoustical Society of America* 117(5) 2005, pp 2835-2843.
DOI: 10.1121/1.1887126.

[23] B.F. SHORR. "Analysis of wave propagation in elastic-plastic rods of a variable cross section using direct mathematical modelling". *Archive of Applied Mechanics* 65(8) 1995, pp 537-547.
DOI: 10.1007/BF00789095.