# Numerical Solution of Dynamic Contact Problems

Patrick Le Tallec, Franck Dambakazi, Patrice Hauret

► To cite this version:

Patrick Le Tallec, Franck Dambakazi, Patrice Hauret. Numerical Solution of Dynamic Contact Problems. Computational Methods in Structural Dynamics and Earthquake Engineering, Jun 2007, Rethymno, Greece. hal-00175626

HAL Id: hal-00175626

https://hal.science/hal-00175626v1

Submitted on 16 Jan 2020

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 CC BY 4.0 - Attribution - International License
# NUMERICAL SOLUTION OF DYNAMIC CONTACT PROBLEMS

Patrick Le Tallec¹, Franck Dambakizi², and Patrice Hauret²

¹Ecole Polytechnique
91 128 Palaiseau Cedex, France
e-mail: patrick.letallec@polytechnique.fr

² Ecole Polytechnique and CEA DAM
91 128 Palaiseau Cedex, France
e-mail: franck.dambakizi@lpolytechnique.edu

³MFP Michelin
63040 Clermont Ferrand Cedex, France
e-mail: patrice.hauret@fr.michelin.com

**Keywords:** contact dynamics, subscale modelling, conservative schemes, lumped modified mass matrix.

**Abstract.** The present paper describes modelization and discretization issues which seem to play an important role in the numerical modeling of contact problems in solid dynamics.

The first aspect concerns subscale modeling, which is crucial when dealing with frictional contact. We will see how the use of one dimensional implicit thermomechanic models resolved in a subgrid inside each boundary cell leads to practical and physically accurate solutions. This model is controlled by the average velocity of the cell, but can reproduce the detailed aspect of the temperature field and elastoplastic strains in the interface layer.

The next issue concerns time discretisation scheme. The proper approximation in time of elastic and contact efforts is a major issue in order to enforce energy conservation or persistency during contact in an exact or controlled way. We review in the paper a technique proposed in [1] which treats the contact by penalty and enforces energy correction techniques to all penalty and energy terms present in the problem formulation.

A last key point concerns the treatment of inertia terms. The numerical solution of a dynamic contact problem often predicts contact pressures which have spurious oscillations both in space and in time at a scale related to the discretisation grid. As described in [2], the stability of standard algorithms and the regularity of the contact pressures is improved by introducing a modified mass matrix in which the nodes in potential contact will have no mass, the corresponding mass being affected to the neighboring internal elements. One can prove that this technique guarantees the regularity in time of the space discrete solution [2], and does not affect the convergence rate in the linear case. The numerical results to be presented indicate that these schemes are surprisingly good in practice. The authors would like to thank J.P. Perlat for many helpful discussions.
# 1 INTRODUCTION

Contact problems in dynamics lead to a large variety of open problems and numerical challenges. Many new solutions have been proposed in the recent literature to handle different aspects of this field. In particular, a lot of progress has been made in the understanding of the physical aspect of contact at a nanoscopic scale. The situation is more disappointing when considering the numerical aspect of the problem. Most of the numerical work has been focused on the development of numerical algorithms able to handle the algebraic structure of contact problems. It appears that such algorithms are of limited efficiency when they overlook modelisation and discretisation issues. The present paper describes three such issues which play an important role in the numerical modeling of contact problems in solid dynamics.

The first aspect concerns subscale modeling. This appears to be crucial when dealing with any type of frictional contact. For example, as indicated in [3], when dealing with high speed frictional contact between metallic objects, a very thin layer of hot metal is produced at a sub-micrometric length due to the heat generated by friction. This brings locally the metal close to melting and therefore strongly affects its resistance to friction. Because of this small scale in space, this layer cannot be described by a standard finite element technique. As illustrated in the second section of this paper, the use of one dimensional implicit thermomechanic models resolved in a subgrid inside each boundary cell leads to practical and physically accurate solutions. This model is controlled by the average velocity of the cell, but can reproduce the detailed aspect of the temperature field and elastoplastic strains in the interface layer. The same philosophy could be applied if one would like to resolve microrugosities when modelling rolling tyres on roads, provided there is no strong two-way coupling between the two scales in space.

The next issue concerns space and time discretisation scheme. Monotonicity and stability can be greatly improved by using adequate mortar elements formulations imposing in a weak way the non penetration condition on positive nodal functions [4], [5]. But the proper approximation in time of elastic and contact efforts is also a major issue in order to enforce energy conservation or persistency during contact in an exact or controlled way. This is important for the robustness and long term accuracy of the scheme. In section 3, different possible solutions are discussed, including one proposed in [1] which treats the contact by penalty and enforces energy correction techniques to all penalty and energy terms present in the problem formulation.

A last key point concerns the treatment of inertia terms. It was known for a rather long time [6], [1] that a slight delay in the expression of inertia forces could lead to numerical schemes with controllable dissipation terms. What has been recently identified in [2] deals with the space discretisation of the mass operator. Indeed, the numerical solution of a dynamic contact problem often predicts contact pressures which have spurious oscillations both in space and in time at a scale related to the discretisation grid. But, as described in [2], the stability of standard algorithms and the regularity of the contact pressures is improved by introducing a modified mass matrix in which the nodes in potential contact will have no mass, the corresponding mass being affected to the neighboring internal elements (Section 4). One can prove that for regular problems not involving contact, this will not affect the convergence rate of the discretisation scheme [7]. On the other hand in presence of contacts, this technique guarantees the regularity in time of the space discrete solution [2]. As shown in section 4, for a good combination of time and space step (which is still governed by intuition), these schemes are surprisingly good in practice, although the theory is far to be fully understood.
# 2 SUBGRID MODELLING OF HIGH SPEED FRICTION

## 2.1 The basic problem

High speed dynamic friction refers to the physics that governs the tangential force acting across a material interface say after the passage of a shock wave. It plays a key role in explosively-driven systems where metal interfaces are submitted to large sliding velocities, large shear forces and large contact pressures. From a theoretical point of view, high speed frictional properties remain largely unknown although various experimental techniques have recently been developed [8]. Severe loading conditions lead to extremely localized thermomechanical processes. These local phenomena may profoundly change the state of interface during the dynamic slip process. Since the sheared layer is very thin (0.1µm) and the sliding speeds are relatively high, important strain rates are obtained near the interface with major irreversible plastic deformation. Metal interfaces are also submitted to a warm-up phase due to frictional heat production, and plastic dissipation. It leads to an important and localized increase in surface temperature, that may lead in extremely short time to the fully melt temperature regime and to the formation of a thin molten metal film.

Although they are widely used, the current generation of hydrocodes designed to handle problems involving high energy, explosives and shock propagation either neglect friction or rely on simple empirical models in which the frictional stress is related to the normal stress and/or the sliding velocity through a user's defined Coulomb's law. A solution would be to apply the regularisation law of [9]. Nevertheless, such laws do not really handle situations where the shear stress is limited by yielding conditions strongly depending on the thermomechanical history at the interface. Because of strong localisation effects, this history is not represented inside standard finite element techniques. A simple solution is to integrate 1D implicit finite difference sub-grid models accounting for frictional contact, elastoplastic behavior, yielding and work hardening, heating by friction and plastic work, heat diffusion, thermal softening and melting. The subgrid model introduces a thin shear layer of micrometric thickness *h* at the interface (Figure 1)

<figure><img src="image_2.png" /><figcaption>Figure 1: Interface shear layer (out of [10])</figcaption></figure>

Inside the layer, we want to predict the time and space evolution of the slip velocity $u_{\tau}(t)$, the stress tensor $\underline{\sigma}(\underline{x}, t)$ whose deviatoric part is supposed to be in pure shear $\underline{s} = \tau(\underline{e}_n \otimes \underline{e}_{\tau} + \underline{e}_{\tau} \otimes \underline{e}_n)$, the strain rate tensor $\underline{d}(\underline{x}, t)$ with its elastic and plastic components $\underline{d}_e$ and $\underline{d}_p$ and the temperature field $\Theta(\underline{x}, t)$ using the equations of motion, the heat equation, the elastoplastic
constitutive laws including pressure and kinematic hardening and thermal softening

$$ \dot{\underline{s}}(\underline{x}, t) = 2G \left( \underline{d}(\underline{x}, t) - \underline{d}_p(\underline{x}, t) \right) \quad (\text{elasticity}) \quad (1) $$

$$ f(\underline{x}, t) = \sigma_{eq}(\underline{x}, t) - Y(\underline{x}, t) \le 0, \quad \text{with } \sigma_{eq} = \sqrt{\frac{3}{2}} \|\underline{s}(\underline{x}, t)\| \quad (\text{yield condition}) \quad (2) $$

$$ \underline{d}_p(\underline{x}, t) = \dot{\lambda} \frac{\partial f}{\partial \underline{\sigma}(\underline{x}, t)} = \frac{3}{2} \left( \frac{\underline{s}(\underline{x}, t)}{\sigma_{eq}(\underline{x}, t)} \right)_+ \cdot \dot{\varepsilon}_p(\underline{x}, t) \quad (\text{plastic flow rule}) \quad (3) $$

$$ Y(\underline{x}, t) = Y_0 (1 + \beta \varepsilon_p(\underline{x}, t))^\eta (1 + gP - h(\Theta(\underline{x}, t) - \Theta_0)) \exp \left( -0.001 \frac{\Theta(\underline{x}, t)}{\Theta_f - \Theta(\underline{x}, t)} \right) \quad (\text{temperature dependent yield condition tending to zero at melting temperature } \Theta_f), \quad (4) $$

and adequate initial and boundary conditions. Temperature has a strong influence on the yield limit, and therefore strongly affects the mechanical problem. On the other hand, the friction heat governs the temperature evolution through the boundary conditions

$$ \dot{q}(0, t) = u_{\tau}(t) \cdot \tau(0, t), \quad \dot{q}(-h, t) = 0. \quad (5) $$

To simplify the local model, we neglect in a first step inertial effects and tangential derivatives which implies from the local equations of motion that the stresses are constant through the layer. We also assume that the velocities at both ends of the cell (slip regime) are given by the macroscopic velocity of the cell. By integrating the elastic constitutive law through the thickness, we then obtain a relation between the time derivative of the shear stress and the integral of the plastic deformation (when assuming equality of velocities at both ends of the cell)

$$ h \dot{\underline{s}}(t) = -2G \underline{d}_p = -3G \int_{-h}^{0} \frac{\underline{s}(t)}{\sigma_{eq}} \dot{\varepsilon}_p(z, t), \quad (6) $$

the local evolution of the plastic deformation being governed by the plastic flow rule when the shear reaches its yield limit.

## 2.2 Numerical scheme

Using a time implicit discretisation scheme, inheriting the contact pressure $P^{n+1}$ and the slip velocity $v_{slip}^{n+1}$ from the global finite element model, and knowing the values of the internal variables $\underline{s}^n, \varepsilon_p^n(z), \Theta^n(z)$ at the previous time step, the local problem reduces to the following scalar problem in $|\underline{s}^{n+1}|$

$$ h (\underline{s}^{n+1} - \underline{s}^n) + 3G \cdot \frac{\underline{s}^{n+1}}{\sigma_{eq}^{n+1}} \int_{-h}^{0} (\Delta \varepsilon_p(z))^{n+1} dz = 0. \quad (7) $$

In this equation, the equivalent stress is given by $\sigma_{eq}^{n+1} = \sqrt{\frac{3}{2}} \|\underline{s}^{n+1}\|$, the local increase $\Delta \varepsilon_p(z)^{n+1} = \varepsilon_p^{n+1}(z) - \varepsilon_p^n(z)$ is a function of the stress intensity through the yield criteria

$$ \begin{aligned} \sigma_{eq}^{n+1} &\le Y^{n+1}(z), \quad \Delta \varepsilon_p(z)^{n+1} \ge 0, \quad (\sigma_{eq}^{n+1} - Y^{n+1}(z)) \Delta \varepsilon_p(z)^{n+1} = 0, \\ Y^{n+1}(z) &= \frac{Y_0}{C_T^{n+1}} (1 + \beta \varepsilon_p^{n+1}(z))^\eta (1 + gP^{n+1} - h(\Theta^{n+1}(z) - \Theta_0)), \\ C_T^{n+1} &= \exp \left( 0.001 \frac{\Theta^{n+1}(z)}{\Theta_f - \Theta^{n+1}(z)} \right). \end{aligned} $$
On the other hand, the temperature is obtained by solving the following one dimensional partial differential equation in space with homogeneous Neumann boundary condition at the lower boundary

$$c \frac{\Theta^{n+1} - \Theta^n}{\Delta t}(z) - k \Delta \Theta^{n+1}(z) = s^n(z), \forall z \in (0, h), \quad -k \frac{\partial \Theta^{n+1}}{\partial z}(-h) = 0, \quad (8)$$

and heterogeneous Neumann boundary condition at the interface

$$-k \frac{\partial \Theta^{n+1}}{\partial z}(0) = \dot{q}_{int}^{n+1} = \beta v_{slip}^{n+1} \cdot \left( \frac{\sigma_{eq}^{n+1}}{\sqrt{3}} \right).$$

The coefficient $\beta$ corresponds to the part of the heat flux which is transmitted to the lower layer.

This scalar equation in equivalent stress could be solved in theory by a secant method, each evaluation of the function requiring the solution of a one dimensional partial differential equation in space to compute the temperature and an integration of plastic strain through the thickness. In practice, the differential equation in temperature is solved by a finite difference technique, and the integration uses a trapezoidal quadrature rule using the same nodes as for the temperature.

This numerical strategy turns out to be unstable next to the melting temperature (taken here to be at 1900 degrees Kelvin), because there is no guarantee that the temperature would not exceed the melting temperature during the iterative process. A much more stable technique is to solve the same scalar equation with respect to the inverse of the thermal softening factor

$$C_T^{n+1} = \exp \left( \frac{0.001}{\frac{\Theta_f}{\Theta^{n+1}(0)} - 1} \right) \in (0, \infty).$$

The interface temperature is an explicit function of this new unknown and the internal temperature can be obtained by solving the differential equation (8) with Dirichlet boundary condition

$$\Theta^{n+1}(0) = \frac{\Theta_f \cdot \log(C_T^{n+1})}{0.001 + \log(C_T^{n+1})}.$$

From this, we deduce the heat flux $\dot{q}_{int}^{n+1} = -k \frac{\partial \Theta^{n+1}}{\partial z}(0)$, the shear stress $\tau^{n+1} = \frac{\dot{q}_{int}^{n+1}}{\beta v_{slip}^{n+1}}$, and the plastic strain

$$\sigma_{eq}^{n+1} \le Y^{n+1}(z), \quad \Delta\varepsilon_p(z)^{n+1} \ge 0, \quad (\sigma_{eq}^{n+1} - Y^{n+1}(z)) \Delta\varepsilon_p(z)^{n+1} = 0, \\ Y^{n+1}(z) = \frac{Y_0}{C_T^{n+1}} (1 + \beta\varepsilon_p^{n+1}(z))^{\eta} (1 + gP^{n+1} - h(\Theta^{n+1}(z) - \Theta_0)).$$

Hence the scalar function (7) is indeed a scalar function of $C_T$ whose evaluation also requires the solution of a one dimensional partial differential equation in space and an integration of plastic strains through the thickness, and whose root is easily obtained by a secant method.

## 2.3 Numerical results

The above technique has been applied to a case with constant normal pressure and slip velocity, taking as initial conditions a state of stress in shear reaching the elastic limit at temperature $\Theta_0$. The simulation was carried with an imposed slip velocity of 34m/s, on a time interval of 0, 6µs with a time step of $\Delta t = 10^{-9}$s, on a shear layer of $h = 10\mu m$ with two hundred discretization points inside the layer. The normal pressure was of 1.3 GPa.
<figure><img src="image_3.png" /><figcaption>Figure 2: Cross section in space of the yield limit at different times. The top part corresponds to region where the material is in plastic flow. Results are displayed in international units (SI).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Time evolution of the interface temperature and shear between 0 and 0.6µs.</figcaption></figure>
# 3 ENERGY CONSERVING SCHEMES

## 3.1 Basic time integration schemes

The next issue in the numerical simulation of contact problem in dynamics deals with the design of an appropriate time discretisation scheme. To address this issue, we consider the simpler case of an elastic body in frictionless contact with a rigid obstacle, whose governing equations define the position $\underline{x}(M, t)$ of the different material points $M$ at time $t$ by:

$$m(\ddot{\underline{x}}, \underline{\hat{U}}) + a(\nabla \underline{x}, \nabla \underline{\hat{U}}) = \int_{\Omega} \underline{f} \cdot \underline{\hat{U}} + \int_{\partial\Omega} \underline{g} \cdot \underline{\hat{U}} + \int_{\partial\Omega_c} \lambda \underline{\nu} \cdot \underline{\hat{U}}, \forall \underline{\hat{U}} \in \mathcal{U}, \quad (9)$$

$$\underline{x} \cdot \underline{\nu} \geq g_0, \lambda \geq 0 \text{ and } \lambda(\underline{x} \cdot \underline{\nu} - g_0) = 0 \text{ on } \partial\Omega_c. \quad (10)$$

Above, the structural mass operator $m$ has the usual linear expression encountered in Lagrangian dynamics

$$m(\ddot{\underline{x}}, \underline{\hat{U}}) = \int_{\Omega} \rho \ddot{\underline{x}} \cdot \underline{\hat{U}} dx.$$

The stiffness term $a(\nabla \underline{x}, \nabla \underline{\hat{U}})$ is defined as $a(\nabla \underline{x}, \nabla \underline{\hat{U}}) = \int_{\Omega} \underline{F} \cdot \underline{\Sigma}(\nabla \underline{x}) \cdot \nabla \underline{\hat{U}}$ with $\underline{\Sigma}$ the second Piola-Kirchhoff stress tensor given by:

$$\underline{\Sigma} = 2W_{,C}, \quad (11)$$

where $W$ denotes the stored elastic energy, which is a given function of the right Cauchy-Green strain tensor $\underline{C} = \underline{F}^t \cdot \underline{F}$, and $\underline{F} = \nabla \underline{x}$ denotes the deformation gradient. The frictionless contact constraint $\underline{x} \cdot \underline{\nu} \geq g_0$ is imposed on a part $\partial\Omega_c$ of the domain boundary where the distance $\underline{x} \cdot \underline{\nu}$ of the material point to a given obstacle cannot get below a given threshold $g_0$. In practice, this frictionless contact constraint is often handled by a penalty approach giving the normal reaction $\lambda$ as a function of the interpenetration distance $|\underline{x} \cdot \underline{\nu} - g_0|_-$ = $\max(0, g_0 - \underline{x} \cdot \underline{\nu})$ by $\lambda = \frac{1}{\epsilon_c}|\underline{x} \cdot \underline{\nu} - g_0|_-$, where $\epsilon_c$ is a small penalty coefficient.

A standard implicit scheme in elastodynamics uses a trapezoidal rule for time integration of the acceleration $\frac{\underline{u}^{n+1} - \underline{u}^n}{\Delta t_n} = \frac{1}{2}(\ddot{\underline{x}}^{n+1} + \ddot{s}\underline{x}^n)$, of the velocity $\frac{\underline{x}^{n+1} - \underline{x}^n}{\Delta t_n} = \frac{1}{2}(\underline{u}^{n+1} + \underline{u}^n)$ and of stresses. For nonlinear problems, Simo or Crisfield [12, 13] have proposed to use in addition a transport averaging, which reduces the time integral of the stiffness term in (9) to $\frac{1}{2}(\nabla \underline{x}^{n+1} + \nabla \underline{x}^n) \cdot \underline{\Sigma}^{n+1/2}$. In theory, the above time integration schemes have good properties with respect to energy conservation, achieving second order accurate conservation, with an error vanishing at the linear limit. In practice, for both the trapezoidal and the mid point schemes, such a second order conservation is not good enough for nonlinear structures, and numerical instabilities are often observed in real life simulations [1].

Nonlinear corrections are then needed, with different choices proposed in the literature. We have tested and adopted a nonlinear and non symmetric correction term proposed by Gonzalez [14], where the elastic stress appearing in the time integration of the stiffness term is given by:

$$\underline{\Sigma}^{n+1/2} = 2 \frac{\partial W}{\partial \underline{C}} (\underline{C}^{n+1/2}) + 2 \left( W(\underline{C}^{n+1}) - W(\underline{C}^n) - \frac{\partial W}{\partial \underline{C}} (\underline{C}^{n+1/2}) : \delta \underline{C}^{n+1/2} \right) \frac{\delta \underline{C}^{n+1/2}}{\delta \underline{C}^{n+1/2} : \delta \underline{C}^{n+1/2}}, \quad (12)$$
with $\delta \underline{C}^{n+1/2} = \underline{C}^{n+1} - \underline{C}^n$. By construction, we then directly have

$$ \frac{1}{2}(\nabla \underline{x}^{n+1} + \nabla \underline{x}^n) \cdot \underline{\underline{\Sigma}}^{n+1/2} : (\nabla \underline{x}^{n+1} - \nabla \underline{x}^n) = \frac{1}{2} \underline{\underline{\Sigma}}^{n+1/2} : \delta \underline{\underline{C}}^{n+1/2} = W(\underline{\underline{C}}^{n+1}) - W(\underline{\underline{C}}^n), $$

implying exact energy conservation, even at the incompressible limit. The resulting numerical tests observed on a simple incompressible beam are then quite convincing and in sharp contrast with the diverging results of the original trapezoidal rule.

## 3.2 Extension to contact problems

But, even after these first two corrections, the proposed scheme does not handle well contact conditions. In the framework of frictionless contact, both Laursen and Chawla [15] and Armero and Petocz [16] observing such difficulites, have shown the interest of the persistency condition $\lambda(t, x) \frac{d}{dt}(\underline{x} \cdot \nu - g_0) = 0$ to obtain energy conservation in the discrete framework. This difficulty can be overcome as in [17] by introducing a discrete jump in velocities during impact. Another solution, proposed in [1], adapts the energy correction (12) in order to enforce energy conservation by averaging separately the geometric update (transport) of the normal $\nu$ and the kinetic force $\lambda$, and by replacing local derivatives by divided differences. In more details, this technique defines the normal vector during the time step by

$$ \underline{\nu}^{n+1/2} = \nu(\underline{x}^{n+1/2}) + [\underline{x}^{n+1} \cdot \underline{\nu}(\underline{x}^{n+1}) - \underline{x}^n \cdot \underline{\nu}(\underline{x}^n) - \underline{\nu}(\underline{x}^{n+1/2}) \cdot \delta \underline{x}] \frac{\delta \underline{x}}{\delta \underline{x} \cdot \delta \underline{x}}, $$

where $\underline{\nu}(\underline{x}_{n+1/2})$ is the normal outward unit vector to the obstacle at mid point $\underline{x}^{n+1/2}$ and $\delta \underline{x} = \underline{x}^{n+1} - \underline{x}^n$ is the displacement update between two successive time steps. Observe that we always have by construction

$$ \underline{\nu}^{n+1/2} \cdot \delta \underline{x} = \underline{x}^{n+1} \cdot \underline{\nu}(\underline{x}^{n+1}) - \underline{x}^n \cdot \underline{\nu}(\underline{x}^n) := \delta g, $$

and that for a plane obstacle for which $\underline{\nu}(\underline{x}^{n+1/2}) = \underline{\nu}(\underline{x}^{n+1}) = \underline{\nu}(\underline{x}^n) = \underline{\nu}$, the above construction simply reduces to $\underline{\nu}_{n+1/2} = \underline{\nu}$. Similarly, the reaction force is updated by

$$ \lambda^{n+1/2} = \frac{1}{\epsilon_c \delta g} \left( \frac{1}{2} |\underline{x}^{n+1} \cdot \underline{\nu}^{n+1} - g_0|^2 - \frac{1}{2} |\underline{x}^n \cdot \underline{\nu}^n - g_0|^2 \right), $$

which implies by construction

$$ \int \lambda^{n+1/2} \underline{\nu}^{n+1/2} \cdot \delta \underline{x} = \int \frac{1}{2\epsilon_c} |\underline{x}^{n+1} \cdot \underline{\nu}^{n+1} - g_0|^2 - \frac{1}{2\epsilon_c} |\underline{x}^n \cdot \underline{\nu}^n - g_0|^2, $$

that is perfect conservation of the penalty energy. The numerical efficiency of this technique is assessed in [1].

# 4 MODIFIED MASS LUMPING METHODS

## 4.1 Motivation

Even when using stable monotone finite element discretisation of the frictionless contact problem (that we may write in a mixed form introducing mortar elements to represent contact stresses as in [4]) and energy conserving schemes, the numerical solutions may exhibit large spurious oscillations of the contact pressures, both in time and in space, directly related to the
discretisation steps in time and space [2, 7]. Such a behavior indicate a weak convergence of these contact pressures. In fact, this can be related to a deep mathematical problem : there is no general theory on the solvability of variational inequalities involving dynamic terms [18]. In specific cases, existence of solutions can be obtained but with very little regularity [19]. After space discretisation, using a standard discretization of the inertia terms, and imposing the contact persistency condition [15], solutions of the time continuous problem can be proved to exist, but they are not uniquely defined and have little regularity [20].

We therefore have a theoretical problem even in the absence of friction and of time discretiza-
tion errors. A practical solution has been recently proposed by H.B. Khenous, P. Laborde et Y.
Renard [21, 22] by redistributing nodal masses next to the contact region in order to remove
all mass from nodes in potential contact with the obstacle. This decomposes the initial time
evolution problem into a standard elastodynamic problem set on all internal nodes coupled to
static variational inequalities to be satisfied at nodes in potential contact.

4.2 Modified lumped mass matrix

In a variational framework, a consistent finite element mass redistribution can be obtained
by defining the new mass operator by

$$
m_h(\underline{x}_h, \underline{\hat{U}}_h) = \sum_T m_T(\mathcal{I}_T \underline{x}_h, \mathcal{I}_T \underline{\hat{U}}_h), \quad (13)
$$

the discrete interpolation operator $\mathcal{I}_T$ being defined on each finite element $T$ from the local finite element space $(\mathcal{U}_h)_T$ into itself with the following properties :

1. $\mathcal{I}_T$ is equal to the identity except for triangles having a node in potential contact,

2. $\mathcal{I}_T \underline{\psi}_k = 0$ for all nodal functions $\underline{\psi}_k$ associated to nodes $k \in \mathcal{K}$ in potential contact.

3. $\mathcal{I}_T$ is $L^2$ stable and preserves constant functions.

In matricial form, by introducing the matrices $I_T$ and $M_T$ of the interpolation operator and
local mass matrix in the finite element nodal basis, the local modified mass matrix is given by

$$
\bar{M}_T = I_T M_T I_T^t.
$$

The first property guarantees that the sparsity of the modified mass matrix $\bar{M}$ is unchanged compared to $M$. It imposes that for all nodes $i$ and $j$ belonging to the set $\mathcal{J}$ of nodes which are not in potential contact, the coefficient $(I_T)_{ij}$ of $I$ must be zero if $i \neq j$ and 1 for $i = j$. The second property guarantees that when ordering the nodes in order to place the nodes in potential contact at the end of the list, the modified mass matrix $\bar{M}$ takes the block form

$$
\bar{M} = \begin{pmatrix} \bar{M}_{\mathcal{J}\mathcal{J}} & 0 \\ 0 & 0 \end{pmatrix}
$$

The $L^2$ stability is satisfied by taking all coefficients of $I_T$ to be non negative. Preserving constant functions requires that the coefficients of a given column be of unit sum

$$
\sum_i (I_T)_{ij} = 1.
$$

It guarantees mass conservation at the finite element level.
Our choice uses as interpolating coefficients the weighted inverses of the distances to the contact surface. If $k \in \mathcal{K}$ (respectively $i \in \mathcal{J}$) is the number of a node M in (respectively not in) potential contact, we define

$$ (I_T)_{ik} = dist(M_i, M_k)^{-1} / \left( \sum_{j \in \mathcal{J} \cap T} dist(M_j, M_k)^{-1} \right). $$

More complex choices are proposed in [2, 7].

## 4.3 Theoretical results

After space discretisation and modification of the mass matrix, the problem takes the matrix form

$$ \bar{M}_{\mathcal{J}\mathcal{J}} \ddot{X}_{\mathcal{J}}(t) = -K_{\mathcal{J}\mathcal{J}} X_{\mathcal{J}}(t) - K_{\mathcal{J}\mathcal{K}} X_{\mathcal{K}}(t) + F_{\mathcal{J}}(t) \quad (14) $$

$$ K_{\mathcal{J}\mathcal{K}} X_{\mathcal{K}}(t) - \sum_{k \in \mathcal{K}} \lambda(k) \underline{\psi}_k \cdot \underline{\nu}(M_k) = -K_{\mathcal{K}\mathcal{J}} X_{\mathcal{J}}(t) \quad (15) $$

$$ \underline{x}(M_k, t) \cdot \underline{\nu}(M_k) \leq g_0(M_k), \lambda(k) \geq 0, \quad (16) $$

$$ (\underline{x}(M_k, t) \cdot \underline{\nu}(M_k) - g_0(M_k))\lambda(k) = 0, \forall k \in \mathcal{K}. \quad (17) $$

Here, we have separated all vectors and matrices in components $\mathcal{J}$ (respectively $\mathcal{K}$) associated to nodes $j$ (respectively $k$) which are not in (respectively in) potential contact. Because of the absence of any inertia term, the last three lines define a well posed elliptic variational inequation with unknown $X_{\mathcal{K}}(t)$ and data $X_{\mathcal{J}}(t)$. Hence $X_{\mathcal{K}}(t)$ is a Lipschitz continuous function of $X_{\mathcal{J}}(t)$. The first equation is thus a second order Lipschitz ordinary differential equation with respect to the unknown $X_{\mathcal{J}}(t)$. Its solution exists, is unique, is differentiable in time, with a Lipschitz continuous time derivative [2].

Moreover, for linear problems in dynamics, hence for problems not involving contacts, one can prove that the local modification of mass matrix on the boundary $\partial\Omega_c$ does not affect the optimal convergence rate of the underlying finite element approximation [7], if the local interpolation operator satisfies the properties enumerated in the previous section. For first order elements in space with shape regular elements of maximal diameter $h$ and second order time discretization schemes with time step $\Delta t$, the error between the fully discrete solution and the continuous solution stays of order $h + \Delta t^2$.

## 4.4 Numerical simulations

As a simple test case, we consider a cube of 10cm width impacting a plane wall at 10 km/h. The tests deal with nonlinear materials, using a trapezoidal rule for integrating the inertia terms and a backward Euler scheme for the integration of stiffness and contact terms.

The figures show the time evolution of the penetration energy during time with and without mass redistribution. Without redistribution, the results display chaotic oscillations in time which stay alive even after refining the time steps. Animations to be presented at the conference confirm the presence (respectively the absence) of time and space instabilities on the contact zone in the case without (respectively with) mass redistribution.

# 5 CONCLUSIONS

We can see from the above that the important points in contact problems more or less identified and deal with subscale modeling, approximation strategies and treatment of inertia terms.
<figure><img src="image_5.png" /><figcaption>Figure 4: Penetration energies for a cube impacting a plane wall without mass redistribution (left) and with mass redistribution (right).</figcaption></figure>

Partial remedies are available in order to handle each of this point. But a global picture of the problem is lacking. In particular, we would need to know in our numerical simulations which local regularisation should be used in time and which is the right limiting mathematical model.

## REFERENCES

[1] P. Hauret, P. Le Tallec. Energy controlling time integration methods for nonlinear elastodynamics and low velocity impact, *Comput. Methods Appl. Mech. Engrg, 27*, 31-38, 2006.

[2] H.B. Khenous, P. Laborde, Y. Renard, A energy conserving approximation for elastodynamic contact problems. C.A Mora Soares and al. eds., *Third European Conference on Computational Mechanics, Solids, Structures and Coupled Problems in Engineering, Lisbon, Portugal, June 5-8 2006*.

[3] M.A. Irfan, V. Prakash, Time resolved friction during dry sliding of metal on metal, *Int. Journal of Solids and Structures, 37*, 2859-2882, 2000.

[4] S. Hueber and B.I. Wohlmuth, A primal dual active set strategy for nonlinear multibody contact problems. *Comp. Meth. in Appl. Mech. and Eng., 194*, 3147-3166, 2005.

[5] T.A. Laursen, *Computational contact and impact mechanics*, Springer Verlag, 2002.

[6] H.M. Hilber and T.J.R. Hughes and R.L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics. *Earthquake Engin. and Struct. Dynamics*, **5**, 283-292, 1977.

[7] C. Hager, B.I. Wohlmuth, Analysis of a modified mass lumping method for the stabilisation of frictional contact problems. submitted to CMAME

[8] G.J. Ball, R.E. Winter and P.T. Keightey, Mechanics of shock induced dynamic friction, *J. of Physics D.*, 2006.

[9] J.T. Oden and J.A.C. Martins. Models and computational methods for dynamic friction phenomena. *Comput. Methods Appl. Mech. Engrg., 52*, 527-634, 1985.
[10] A. Juanicotena, Experimental investigation of dynamic friction at high contact pressure applied to an aluminium/stainless steel tribo pair, *J. of Physics IV*, **134**, 559-564 2006.

[11] K-J. Bathe. Finite Element Procedures in Engineering Analysis, Prentice-Hall, 1982.

[12] J.C. Simo and N. Tarnow. The discrete energy-momentum method. conserving algorithms for non linear elastodynamics. *Z angew Math Phys*, **43**, 757–792, 1992.

[13] M.A. Crisfield. Nonlinear Finite Element Analysis of Solids and Structures, Wiley, 1997.

[14] O. Gonzalez, Exact energy and momentum conserving algorithms for general models in nonlinear elasticity. *Comput. Methods Appl. Mech. Engrg.*, **190**(13-14), 1763–1783, 2000.

[15] T.A. Laursen and V. Chawla. Design of energy conserving algorithms for frictionless dynamic contact problems. *International Journal for Numerical Methods in Engineering*, **40**, 863–886, 1997.

[16] F. Armero and E. Petcz. Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems. *Comp. Meth. in Appl. Mech. and Eng.*, **158**, 269–300, 1998.

[17] T.A. Laursen and G.R. Love. Improved implicit integrators for transient impact problems; geometric admissibility within the conserving framework. *Int. J. Num. Meth. Engr.*, **53**(2), 245–274, 2002.

[18] G. Duvaut and J.L. Lions. Les inéquations en Mécanique et en Physique, Dunod, 1972.

[19] P. Ballard and E. Basseville. Existence and uniqueness for dynamical unilateral contact with Coulomb friction : a model problem. *M2AN*, **39** (1), 59–77, 2005.

[20] J.J. Moreau, Numerical aspects of the sweeping process. *Comput. Methods Appl. Mech. Engrg.*, **177**, 329–349, 1999.

[21] H.B. Khenous, P. Laborde, and Y. Renard., Comparison of two approaches for the discretization of elastodynamic contact problems. *CRAS Mathématiques*, **342/10**, 791–796, 2006.

[22] H.B. Khenous, P. Laborde, and Y. Renard., On the discretization of contact problems in elastodynamics. *Lecture Notes in Applied Comp. Mech.*, **27**, 31–38, 2006.