# Consistent time discretization for a dynamical frictional contact problem and complementarity techniques

Didier Vola, Elaine Pratt, Michel Jean, Michel Raous

► To cite this version:

Didier Vola, Elaine Pratt, Michel Jean, Michel Raous. Consistent time discretization for a dynamical frictional contact problem and complementarity techniques. Revue Européenne des Éléments Finis, 1998, 10.1080/12506559.1998.11690471. hal-03178272

HAL Id: hal-03178272

https://hal.science/hal-03178272

Submitted on 23 Mar 2021

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.
# Consistent time discretization for a dynamical frictional contact problem and complementarity techniques

Didier Vola, Elaine Pratt, Michel Raous

## ► To cite this version:

Didier Vola, Elaine Pratt, Michel Raous. Consistent time discretization for a dynamical frictional contact problem and complementarity techniques. Revue Européenne des Éléments Finis, HERMÈS / LAVOISIER, 1998. hal-03178272

HAL Id: hal-03178272

https://hal.archives-ouvertes.fr/hal-03178272

Submitted on 23 Mar 2021

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

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.
Consistent time discretization for a
dynamical frictional contact problem and
complementarity techniques

Didier Vola, Elaine Pratt, Michel Jean and Michel Raous

Laboratoire de Mécanique et d'Acoustique
31, Chemin Joseph Aiguier, 13402 Marseille Cedex 20, France

RÉSUMÉ. Le problème dynamique de contact unilatéral avec frottement de Coulomb en viscoélasticité linéaire est formulé de manière à prendre en compte avec précision les éventuelles discontinuités des vitesses relatives lors d'un impact. Plusieurs schémas numériques implicites du premier ordre sont proposés pour l'approximation en temps et ils sont comparés à un schéma standard d'ordre plus élevé. Le problème discret, formulé en vitesse, est écrit en terme de complémentarité et résolu par la méthode de programmation mathématique de Lemke.

ABSTRACT. The dynamical problem of a viscoelastic body involving unilateral contact and Coulomb friction is set so as to take into account accurately eventual discontinuities of the relative velocities when impacts occur. Several first order implicit numerical schemes for the time discretized equations of the problem are proposed. The results are compared to those of a higher order standard numerical scheme. The discrete problem, with the velocity as unknown, is set in terms of a complementarity problem which is solved by Lemke's mathematical programming method.

MOTS-CLÉS : contact unilatéral, impact, frottement de Coulomb, dynamique des structures, θ-méthode implicite, éléments finis, problème de complémentarité linéaire, programmation mathématique, méthode de Lemke.

KEY WORDS : unilateral contact, impact, Coulomb friction, structure dynamics, implicit θ-method, finite elements, linear complementarity problem, mathematical programming, Lemke's method.
# 1. Introduction

Dynamical problems involving frictional contact are not regular because relative velocities are liable to become discontinuous when impulses are exerted between contacting bodies during impact. In this paper, only finite dimensional systems are considered. A finite element method, for instance, can be used to obtain a finite dimensional model of the continuous dynamical problem, however, such a discretization is not within the scope of this paper. The discrete model consists of a second-order differential equation with respect to time governing the displacement of the mesh nodes and of the frictional contact relations. Since impacts are expected, such a differential equation should be understood in the sense of distributions or of measure differential equations. Any velocity $v$ is chosen as a function with bounded variation, $v^{-}$ is the left limit representing the velocity before some possible impact time $t$ whereas $v^{+}$ is the right limit representing the velocity after the impact time. This implies that the acceleration is a differential measure. The reaction is represented by some positive measure involving, for instance, the Dirac measure at impact time. Several numerical schemes for the time discretization of this measure differential equation may be adopted, each one depending on different choices for the approximation of the displacement and of the terms representing internal and external forces. To the dynamical equation are added the frictional contact relations involving impulses and the relative velocities at the impact time. The unilateral condition appears as a complementarity condition, involving the gap and the reaction force, referred to as the Signorini condition. The Coulomb friction law is written using the right sliding velocity and the density of impulse. The use of the right sliding velocity, is motivated by numerical experiments on simple examples.

Most authors apply standard schemes, for instance, Newmark's scheme, to integrate the dynamical equation, as if relative velocities were sufficiently regular. Hughes et al. have proposed to correct the solution obtained with Newmark's scheme when an impact (or a release) occurs using a wave propagation analysis [?]. More recently, Taylor and Papadopoulos have preferred to add a rate and a second rate form of the displacement impenetrability condition [?] and add two Lagrange multipliers to enforce these new conditions (augmented lagrangian formulation).

As the time discretized problem for each time iteration is very similar to a stationary frictional contact problem, except that the primary unknown here is the velocity, the Lemke's mathematical programming method is used to solve at each time step the complete set of equations, the approximated dynamical equation, together with the complementarity frictional contact equations, [?].

Though the ideas developed in this paper may be applied straightforwardly to general large deformation dynamical problems with frictional contact (see section 6), for simplicity's sake, the case of linear viscoelasticity under small perturbations is developed. The simple example of the impact of two bars is presented [?], and used as a benchmark: the results obtained by the different numerical schemes proposed by the authors to solve the dynamical equation are compared to those obtained by using Newmark's method and with the exact solution. The results suggest that the proposed numerical schemes are quite
effective in coping with the discontinuities of the velocity. Other less academic examples are also presented showing the feasability of this approach on large deformation problems.

## 2. Space Discretization

The configuration of a continuous medium is represented by a n dimensional vector $q(t)$ whose components are the co-ordinates of the mesh nodes obtained by a finite element discretization (P1 in our case). The local gap $g(t)$, the relative velocity $V(t)$ and the local impulse reaction $R(t)$ exerted at some point candidate for contact, for instance a mesh node, are components of some $2n_c$ or $3n_c$ dimensional vectors where $n_c$ is the number of candidates. Contact and friction conditions are expressed through relations between $g(t)$, $V(t)$ and $R(t)$. Normal (respectively tangential) components are noted with a subscript N (respectively T). The discrete problem to be solved is the following:

**Problem P_n:** Find $q$ such that $\forall t \in [0, T]$, $q(0) = q_0$, $\dot{q}(0) = V_0$ and:

$$ M d\dot{q} + K q dt + C \dot{q} dt = F dt + r d\nu, \quad [1] $$

and for the contact nodes:

$$ \begin{gathered} g(t) \le 0, \quad R_N(t) \le 0 \quad \text{and} \quad g(t) R_N(t) = 0, \\ \|R_T(t)\| \le \mu |R_N(t)|, \quad [2] \\ \text{if } V_T^+(t) \ne 0 \text{ then } \|R_T(t)\| = -\mu |R_N(t)| V_T^+(t) / \|V_T^+(t)\|, \end{gathered} $$

where $dt$ is the Lebesgue measure, $d\dot{q}$ is a differential measure representing the discretized acceleration, $dv$ is a non-negative real measure relative to which $d\dot{q}$ happens to possess a density function, and $r$ is a representative of the local density of impulses exerted where contact occurs.

## 3. Time Discretization of the Measure Differential System

It may be useful here to remind the reader that a differential measure $dv$ can be associated to any function $v$ with bounded variation, and that $\int_{[a,b]} dv = v^+(b) - v^-(a)$ and $\int_{[a,b]} dv = v^+(b) - v^-(a)$ where $v^+$ and $v^-$ refer respectively to the right-limit and left-limit of $v$. For more details on the subject see J.J. Moreau [?]. We shall consider that the velocity is a right continuous function of time, so that $v = v^+$. A study of such measure differential systems as [1] can be found in [?]. The authors show that [1] may be written in the following equivalent form: for all $t$ in $]0, T]$

$$ M(\dot{q}(t) - \dot{q}(0)) = \int_{0}^{t} (F - Kq - C\dot{q})ds + \int_{]0,t]} r d\nu, \quad [3] $$

where $ds$ represents the Lebesgue measure, and also: for all $\tau$ in $]0, t[$
$$ egin{cases} M(\dot{q}(t) - \dot{q}(\tau)) = \int_{\tau}^{t} (F - Kq - C\dot{q})ds + \int_{[\tau,t]} r d\nu, \\ q(t) = q(\tau) + \int_{\tau}^{t} \dot{q}ds. \end{cases} \quad [4] $$

In order to obtain a numerical approximation of the solution of [3], the time interval $[0, T]$ is divided into sub-intervals. Considering some sub-interval $[t_i, t_{i+1}]$ of length $h$, [4] yields the following relations:

$$ egin{cases} M(\dot{q}(t_{i+1}) - \dot{q}(t_i)) = \int_{t_i}^{t_{i+1}} (F - Kq - C\dot{q})ds + \int_{[t_i, t_{i+1}]} r d\nu, \\ q(t_{i+1}) = q(t_i) + \int_{t_i}^{t_{i+1}} \dot{q}ds. \end{cases} \quad [5] $$

The next step consists in choosing a numerical scheme to approximate the two integrals $\int_{t_i}^{t_{i+1}} \dot{q}ds$ and $\int_{t_i}^{t_{i+1}} (F-Kq-C\dot{q})ds$. A different numerical method is obtained for each different choice of an integral approximation scheme. However, bearing in mind the fact that the velocity may be discontinuous, the use of high-order approximations should be avoided since they are accurate only for the integration of sufficiently regular functions. Therefore, the classical $\theta$-method is proposed in this paper (i.e., $\int_{t_i}^{t_{i+1}} f ds \approx h(\theta f(t_{i+1}) + (1-\theta)f(t_i))$, $\theta=1$ yields the implicit Euler method). This method is unconditionally stable for $\theta \ge \frac{1}{2}$. The mean value impulse $r^{i+1} = \frac{1}{h} \int_{[t_i, t_{i+1}]} r d\nu$ appears as a primary variable. Setting $\dot{q}^i, q^i, \dot{q}^{i+1}, q^{i+1}$, as approximations of respectively $\dot{q}(t_i)$, $q(t_i)$, $\dot{q}(t_{i+1})$, $q(t_{i+1})$, the following three methods to compute the unknowns $\dot{q}^{i+1}, q^{i+1}$ and $r^{i+1}$ at each time step may be obtained from expression [5].

**$\theta$-Method:** both integrals are approximated by the $\theta$-method, so that one has:

$$ egin{cases} M(\dot{q}^{i+1} - \dot{q}^i) = h[\theta(F^{i+1} - Kq^{i+1} - C\dot{q}^{i+1}) + (1-\theta)(F^i - Kq^i - C\dot{q}^i)] \\ q^{i+1} = q^i + h[\theta\dot{q}^{i+1} + (1-\theta)\dot{q}^i], \end{cases} $$

or

$$ egin{cases} \tilde{M}\dot{q}^{i+1} = \tilde{F}^{i+1} + hr^{i+1}, \quad q^{i+1} = q^i + h[\theta\dot{q}^{i+1} + (1-\theta)\dot{q}^i], \\ \tilde{M} = M + h\theta C + h^2\theta^2 K, \\ \tilde{F}^{i+1} = [M - h(1-\theta)C - h^2\theta(1-\theta)K]\dot{q}^i - hKq^i + h[\theta F^{i+1} \\ \qquad + (1-\theta)F^i]. \end{cases} \quad [6] $$

**$\theta$-Euler-Method:** the first integral is approximated by the $\theta$-method and the other by the implicit Euler method, so that one has:
$$ \left\{ \begin{array}{l} M(\dot{q}^{i+1} - \dot{q}^i) = h[\theta(F^{i+1} - Kq^{i+1} - C\dot{q}^{i+1}) + (1-\theta)(F^i - Kq^i - C\dot{q}^i)] \\ q^{i+1} = q^i + h\dot{q}^{i+1}, \end{array} \right. + hr^{i+1}, $$

or

$$ \left\{ \begin{array}{l} \tilde{M}\dot{q}^{i+1} = \tilde{F}^{i+1} + hr^{i+1}, \quad q^{i+1} = q^i + h\dot{q}^{i+1}, \\ \tilde{M} = M + h\theta C + h^2\theta K, \\ \tilde{F}^{i+1} = [M - h(1-\theta)C]\dot{q}^i - hKq^i + h[\theta F^{i+1} + (1-\theta)F^i]. \end{array} \right. \qquad [7] $$

**modified θ-Method:** both integrals are approximated by the θ-method, but in the contact relations, the displacement $q^{i+1}$ is replaced by $\dot{q}^{i+1} = q^{i+1} + h(1 - \theta)\dot{q}^{i+1}$. So that this method shall only differ from the θ-method when the contact relations are added.

In section 6.1, these three methods have been compared to the exact solution of a simple test example and to the results obtained by Newmark's method for which similar relations to [6] and [7] are obtained.

## 4. Discretization of the Frictional Contact Conditions

The relative velocities $V(t)$ are related to $\dot{q}(t)$ by kinematic relations,

$$ V(t) = H^*(q(t)) \dot{q}(t). $$

The densities of impulse $r(t)$ and $R(t)$ satisfy,

$$ r(t) = H(q(t)) R(t), $$

where $H^*(q(t))$ is the transpose of the linear mapping $H(q(t))$.

The following approximations may be adopted,

$$ V^{i+1} = H^*(q^i) \dot{q}^{i+1}, $$

$$ r^{i+1} = H(q^i) R^{i+1}, $$

where $V^{i+1}$ shall represent an approximation of $V(t_{i+1})$. The approximation of the gap is crucial for the behaviour of the algorithm. It is based on the following property: the time derivative of the gap function $t \mapsto g(t)$ is the normal component of the relative velocity $V_N(t)$,

$$ \frac{d}{dt} g(t) = V_N(t). $$

Using the above property and the different expressions obtained in Section 3 according to the discretization method adopted, the gap at $t_{i+1}$ is approximated by $\dot{q}^{i+1}$ and given by the following expressions:
for the θ-method:

$$
\begin{align*}
g^{i+1} &= g^i + h\theta V_N^{i+1} + h(1-\theta)V_N^i, \\
\frac{1}{h\theta}g^{i+1} &= V_N^{i+1} - \tilde{G}^i, \\
\tilde{G}^i &= \left(1-\frac{1}{\theta}\right)V_N^i - \frac{1}{h\theta}g^i, \tag*{[8]}
\end{align*}
$$

for the θ-Euler-method:

$$
\begin{align*}
g^{i+1} &= g^i + hV_N^{i+1}, \\
\frac{1}{h}g^{i+1} &= V_N^{i+1} - \tilde{G}^i, \\
\tilde{G}^i &= -\frac{1}{h}g^i, \tag*{[9]}
\end{align*}
$$

and for the **modified θ-method** (as the displacement is estimated from the relation $\hat{q}^{i+1} = q^{i+1} + h(1-\theta)\hat{q}^{i+1}$ whereas we still have $q^{i+1} = q^i + h\theta\hat{q}^{i+1} + h(1-\theta)\hat{q}^i$):

$$
\begin{align*}
\hat{g}^{i+1} &= \hat{g}^i + hV_N^{i+1}, \\
\frac{1}{h}\hat{g}^{i+1} &= V_N^{i+1} - \tilde{G}^i, \\
\tilde{G}^i &= -\frac{1}{h}\hat{g}^i. \tag*{[10]}
\end{align*}
$$

In the **modified θ-method**, the estimated gap is not the gap at the end of the time step but at some near future time $t_{i+1}+h(1-\theta)$. For the first two methods, the unknowns involved in the discretized Signorini condition are $g^{i+1}$, $R_N^{i+1}$, whereas, for the **modified θ-method** the unknowns are $\hat{g}^{i+1}$, $R_N^{i+1}$. In all three cases, the discretized Signorini condition is written as:

$$
V_N^{i+1} - \tilde{G}^i \le 0, \quad R_N^{i+1} \le 0 \quad \text{and} \quad (V_N^{i+1} - \tilde{G}^i) R_N^{i+1} = 0. \quad [11]
$$

The advantage of the **θ-Euler-method** and of the **modified θ-method** over the **θ-method**, is the compatibility between gap and velocity approxima-
tions. It can easily be observed from the above relations that if some contact
occurs at $t_i$, i.e., if $g^i = 0$ (or $\hat{g}^i = 0$), then $\tilde{G}^i = 0$ so that if $g^{i+1} = 0$ (or
$\hat{g}^{i+1} = 0$) then $V_N^{i+1} = 0$ also. This is not the case for the **θ-method** where, if
$g^i = 0$, $\tilde{G}^i$ is different from zero so that if the approximated gap at $t_{i+1}$ is equal
to zero, i.e. $g^{i+1} = 0$, one has $V_N^{i+1} = (1 - \frac{1}{\theta})V_N^i - \frac{1}{h\theta}g^i$, which is different
from zero and this introduces some sort of restitution energy.

Finally, the discretized Coulomb’s law involves the unknowns $V_T^{i+1}$, $R_T^{i+1}$,
and is written as:

$$
\begin{align}
\| R_T^{i+1} \| &\le \mu |R_N^{i+1}|, \tag{12} \\
\text{if } V_T^{i+1} \ne 0 \text{ then } R_T^{i+1} &= -\mu |R_N^{i+1}| \frac{V_T^{i+1}}{\|V_T^{i+1}\|}. \nonumber
\end{align}
$$

The advantage of expressing the contact conditions in this manner is that in
both conditions [11] and [12] the velocity appears as a primary unknown.
# 5. Complementarity Formulation and Resolution

## 5.1. The Complementarity Formulation

The unilateral Signorini condition appears naturally as a complementarity relation. Klarbring has shown that the Coulomb friction conditions can also be expressed in these terms (see [?]). In the two-dimensional case, it is natural to introduce four new variables, $\lambda_1$ and $\lambda_2$ (respectively the positive and negative parts of the tangential velocity for the contact nodes) and, $\Phi_1$ and $\Phi_2$ with $\Phi_1 = -R_T + \mu R_N$ and $\Phi_2 = R_T + \mu R_N$ (the three inequations $\Phi_1 \ge 0$, $\Phi_2 \ge 0$ and $-R_N \ge 0$ define the Coulomb cone). For each node in contact, the friction conditions [2] are then equivalent to [14]. When three-dimensional cases are considered the Coulomb cone has to be approximated by a polygonal cone and the number of new variables introduced depends on the number of facets of the polygon. The two complementarity conditions [14] are replaced by a number of similar conditions involving all the new variables introduced.

After having adopted one of the three integration methods of Section 3 and their corresponding approximations of the contact conditions introduced in Section 4, we obtain:

**Problem P_d** : Find $\dot{q}^{i+1}$ and $r^{i+1}$ such as:

$$ \tilde{M} \dot{q}^{i+1} = \tilde{F}^{i+1} + h r^{i+1}, \quad [13] $$

and for the contact nodes:

$$ \left\{ \begin{array}{l} -R_N^{i+1} \ge 0, \\ \tilde{G}^i - V_N^{i+1} \ge 0, \\ -R_N^{i+1} (\tilde{G}^i - V_N^{i+1}) = 0, \end{array} \right. \\ \left\{ \begin{array}{l} \Phi_1^{i+1} \ge 0, \lambda_1^{i+1} \ge 0 \text{ and } \Phi_1^{i+1} \cdot \lambda_1^{i+1} = 0, \\ \Phi_2^{i+1} \ge 0, \lambda_2^{i+1} \ge 0 \text{ and } \Phi_2^{i+1} \cdot \lambda_2^{i+1} = 0. \end{array} \right. \quad [14] $$

The matrix $\tilde{M}$ and the vectors $\tilde{F}^{i+1}$ and $\tilde{G}^i$ depend on the integration method one has adopted and have been defined in the preceding sections. The $\theta$-method corresponds to [6] and [8], the $\theta$-Euler-method to [7] and [9] and the modified $\theta$-method corresponds to [6] and [10].

Finally, in order to formulate the problem in terms of complementarity, the linear system [13] is condensed on all the degrees of freedom of the contact nodes (see [?]). The linear complementarity problem is solved using the pivoting method of Lemke (see [?]). The dimension $\tilde{n}$ of this problem is equal to $3n_c$ in the two-dimensional case. The solution of the linear complementarity problem is then introduced into the linear system [13], thus reducing the dimension of the system to the total number of degrees of freedom minus $\tilde{n}$. The resulting system is solved by a standard method for linear systems.
## 5.2. Mathematical Programming Solver

After condensation, problem $P_d$ is written as a linear complementarity problem of the following general form:

**Problem P_c:** Let $f$ be a given vector $\in \mathbb{R}^{\bar{n}}$, find $w$ and $z \in \mathbb{R}^{\bar{n}}$ such that:

$$ \begin{cases} w - Az = f \\ w_i \ge 0, z_i \ge 0, w_i z_i = 0 \quad \forall i=1, \dots, \bar{n} \end{cases} $$

There are several methods to solve this sort of problem. All are based on the following remark: if $f_i \ge 0, \forall i=1, \dots, \bar{n}$, then the solution is $w_i = f_i$ and $z_i = 0, \forall i=1, \dots, \bar{n}$. We have chosen the method introduced by Lemke (see [?]). A pivoting algorithm is used to construct, by linear combinations, a sequence of non-negative vectors $z$. Details of this algorithm can be found in [?] and [?]. Due to the change of variables (introduction of the new variables $\lambda_1, \lambda_2, \Phi_1$ and $\Phi_2$), the matrix $A$ is non-symmetrical and depends on the friction coefficient. From a mathematical point of view, the copositivity of the matrix $A$ ensures the convergence of the pivoting algorithm. In practice, no convergence problems have been observed. Chabrand et al. (see [?]) have shown the efficiency of the Lemke algorithm in the quasistatic case and compared it to the augmented lagrangian method and the Gauss-Seidel method accelerated by the Aitken procedure.

# 6. Numerical Examples

## 6.1. Impact of Elastic Bars

### 6.1.1. One Dimensional Test

This is a simple test example [?] modeling the impact of two identical elastic bars (see Figure 1). The material properties are given in Table 1.

<figure><img src="image_1.png" /><figcaption>Figure 1. Impact of two bars</figcaption></figure>

Figure 2 represents the evolution of the contact velocity for each of the four methods together with the exact solution which is known in this very simple case. The Newmark method with $\gamma = 0.5$ and $\beta = 0.25$ generates undesirable oscillations of the contact velocity during the contact phase. The same oscillations would have appeared if we had used the $\theta$-method with $\theta = 0.5$.
<table><caption>**Table 1.** *Data for the impact of two bars*</caption><tr><td>Young's modulus</td><td>0.20684.10<sup>12</sup> Pa</td><td>initial velocity</td><td>5.13588m/s</td></tr><tr><td>density</td><td>0.78957.10<sup>4</sup> Kg/m<sup>3</sup></td><td>evolution time</td><td>0.2.10<sup>-3</sup>s</td></tr><tr><td>length L</td><td>0.254m</td><td>time step</td><td>0.2226.10<sup>-5</sup>s</td></tr></table>

These oscillations are damped out by using coefficients that introduce dissipation (see the results obtained here by the $\theta$-Method with $\theta = 0.55$). For the $\theta$-Euler-method, the transition from the zero velocity of the contact phase to the negative velocity is much too slow. The modified $\theta$-method seems to be the most efficient in this simple case but this shall not be so in the two-dimensional case.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** *Contact velocity*</figcaption></figure>
### 6.1.2. Two-Dimensional Test

The data for this example is the same as for the previous one, except for the Poisson coefficient which was then equal to zero and which is now equal to 0.45. No exact solution is available, in this case, to compare the numerical results to. Figure 3 represents the normal velocity, computed by the three methods of Section 3, at the centre of the contact zone of both bars, the dotted line representing the velocity of the impacting bar. In this case, the $\theta$-Euler-method seems a good compromise between the $\theta$-method, where the velocity is too smooth and does not reach the value minus one, and the modified $\theta$-method which generates improbable oscillations.

<figure><img src="image_3.png" /><figcaption>Figure 3. Normal velocity at the centre of the contact zone for both bars</figcaption></figure>

## 6.2. Impact of a Cylinder

This is a two-dimensional frictionless impact test proposed by Taylor and Papadopoulos [?]. The plain strain hypothesis is adopted and the gravity effects are neglected. The geometrical and material properties are given on Figure 4. The mesh has 613 nodes (576 four-node finite elements). The results for
<figure><img src="image_4.png" /><figcaption>Figure 4. Geometry of the cylinder and data</figcaption></figure>

<table><tr><td>Young's modulus</td><td>5.10<sup>2</sup></td></tr><tr><td>Poisson's ratio</td><td>0.3</td></tr><tr><td>density</td><td>1</td></tr><tr><td>total evolution time</td><td>6</td></tr></table>

this example were computed with the $\theta$-Euler-Method. Figure 5 depicts the evolution of the deformed shapes until the release. The total contact force is plotted on Figure 6 for two time-steps ($\Delta t$ = 0.02 and $\Delta t$ = 0.15). For the smaller time-step, the results are similar to those obtained by Taylor and Papadopoulos. Moreover, when a larger time-step is used, the $\theta$-Euler-Method produces a smoother solution than theirs. The computation lasted about 1400s on an IBM RS/6000 590 system.

<figure><img src="image_5.png" /><figcaption>Figure 5. Deformed shapes at $t=0$, $t=2.1$, $t=3.3$ and $t=4.3$</figcaption></figure>

### 6.3. Frictional Oblique Impact

This is a two-dimensional example proposed by Kim and Kwak [?] dealing with the oblique impact of a plate with a round side against a rigid surface. The initial horizontal and vertical velocities are $V_h = 3m/s$ and $V_v = -5m/s$, the material properties and the geometric data are given in the following table. The
<figure><img src="image_6.png" /><figcaption>Figure 6. Evolution of the total contact force</figcaption></figure>

mesh (see Figure 7) is composed of 99 nodes (80 four-node finite elements). The plane stress hypothesis is adopted. Figures 8 and 9 present the nodal normal

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

<table><caption>Figure 7. Geometry of the plate and data</caption><tr><td>Young's modulus</td><td>10<sup>7</sup> Pa</td></tr><tr><td>Poisson's ratio</td><td>0.25</td></tr><tr><td>density</td><td>10<sup>3</sup> Kg/m<sup>3</sup></td></tr><tr><td>friction coefficient</td><td>0.1</td></tr><tr><td>length</td><td>0.08m</td></tr><tr><td>width</td><td>0.04m</td></tr><tr><td>depth</td><td>0.01m</td></tr><tr><td>radius of the round side</td><td>0.101m</td></tr><tr><td>time step</td><td>10<sup>-5</sup> s</td></tr><tr><td>total evolution time</td><td>3.10<sup>-3</sup> s</td></tr></table>

reactions and the evolution of the deformed shapes obtained with the modified θ-Method with θ = 0.55. Kim and Kwak have used the Newmark method associated to a correction procedure when impact occurs. The curves of Figure 8 have the advantage of not presenting the oscillations obtained by Kim and Kwak (p.4621 of [?]) during the contact time. Indeed, these oscillations cannot be explained by the propagation of an acoustic wave in the plate. For this test, the computing time on an IBM RS/6000 590 system was of about 300s.
<figure><img src="image_8.png" /><figcaption>Figure 8. Evolution of the nodal normal contact forces</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9. Deformed shapes at $t=10^{-3}$s, $t=2.10^{-3}$s and $t=2.2510^{-3}$s</figcaption></figure>

6.4. Impact of an Elastoplastic Cylinder

In order to show the feasibility of our approach in the context of large deformation problems, we present in Figure 10 the final shape obtained after the impact of an elastoplastic cylinder against a rigid surface. This is a classical benchmark [?] which is often treated by replacing the contact conditions by bilateral ones. A forthcoming article by Dubois et al. shall present the generalization of our approach to large deformation dynamical problems.

<figure><img src="image_10.png" /><figcaption>Figure 10. Initial mesh and final deformed shape of the elastoplastic cylinder</figcaption></figure>
<figure><img src="image_11.png" /></figure>
## 7. Bibliography

### Acknowledgements:

The authors thank Frédéric Dubois from the Laboratoire de Mécanique et d'Acoustique for the numerical treatment of the last example.
REFERENCE :

Vola, D., Pratt, E., Jean, M., and Raous, M. (1998). Consistent time discretization for a dynamical frictional contact problem and complementarity techniques. Revue Européenne des Eléments Finis, 7, pp.149-162.