# A new contact finite element coupled with an analytical search of contact

Etienne Arnoult, Isabelle Guilloteau, Bernard Peseux, Jérôme Bonini

► To cite this version:

Etienne Arnoult, Isabelle Guilloteau, Bernard Peseux, Jérôme Bonini. A new contact finite element coupled with an analytical search of contact. ECCOMAS 2000, Sep 2000, Barcelone, Spain. hal-01007776

HAL Id: hal-01007776

https://hal.science/hal-01007776

Submitted on 9 Oct 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 4.0 International License
A NEW CONTACT FINITE ELEMENT COUPLED WITH
AN ANALYTICAL SEARCH OF CONTACT

E. Arnoult*, I. Guilloteau*, B. Peseux*, and J. Bonini†

*Laboratoire Mécanique et Matériaux
Ecole Centrale Nantes
BP 92101 F-44321 Nantes cedex 3, France
Etienne.Arnoult@ec-nantes.fr, http://www.ec-nantes.fr/MecaMater/Structures

†Département Méthodes
SNECMA Villaroche F-77550 Moissy-Cramayel, France
Jerome.Bonini@snecma.fr, http://www.snecma.fr

**Key words:** Contact element, Contact search, Newmark implicit scheme, Penalty function method, New contact stiffness, B-splines, analytical formulation

**Abstract.** This paper presents a new way to deal with contact within an implicit finite element code. The contact detection is realised with an analytical method using B-spline functions, and contact analysis is done with a contact element based on penalty function method using a new contact stiffness.
# 1 INTRODUCTION

The study of aircraft engines requires numerical modelisation to try to avoid costly and destructive experimental tests performed on the whole engine. The modelisation of contact phenomena is of particular importance because the results of simulations of blade-off events for instance greatly depend on the way contact is handled during calculations.

The two points we decided to look at precisely concern first the way contact is located during the movement of the structure (geometrical part of our study: see section 2), and then the way we evaluate the contact force (structural part of our study: see section 3).

The methods we propose are implemented in the finite element code Samcef, used by our industrial partner Snecma. The results we obtained are presented in section 4.

# 2 ANALYTICAL SEARCH OF CONTACT

## 2.1 Context

Dynamical contacts in turbojets are dealt with in many studies and numerous research and processing algorithms have been specially developed and implemented in finite element codes so as to improve the simulations of complex movements. However, within a code, the facetisation (which is a direct consequence of the discretisation of structures) stands for a good example of a kind of problems that decrease the performance of the numerical tools. During dynamic studies, especially when contact phenomenon arises, two situations particularly underline the generated difficulties: discontinuity of the normal vector between elements (see Figure 1(a)) and bad estimation of curvatures of the real geometry (see Figure 1(b)).

<figure><img src="image_2.png" /><figcaption>Figure 1: Problems due to spatial discretisation.</figcaption></figure>

Usually, these difficulties are got round using the following artifices: the problem of the indetermination of the normal vector is solved by attributing to node B a normal vector reached by the sum of the normal vectors of adjacent elements; as for effects of the second difficulty, they can be lessened by adding supplementary nodes.

However, these solutions are only mid-terms: in the first case, the normal vector, although defined in every point of the mesh, remains discontinuous when crossing element
boundaries, which is not a sufficient condition in some configurations (for example: ro-
tation of a transmission shaft in a bearing support). In the second case, the addition of
nodes decreases the relative size of the elements and increases the size of the mesh. So it
globally increases the time of calculation, particularly in rapid dynamics where the value
of the time step depends directly on the element size.

The method we suggest in this paper has been motivated by the study of the search
of intersection points of two geometrical entities in the particular framework of studies
of blades/casing contacts in turbojet engines. In this configuration indeed, problems of
facetisation involve a very bad simulation of the rotation of the bladed disk in the casing,
and can lead to a “numerical” break of the two parts due to the fact that the number of
nodes generally retained for the model is too small.

## 2.2 General overview

The suggested method takes place in the stage of research of contact points. More
exactly, it is fully related to the determination of the nodes of the mesh that have not
filled a given criterion. The next stage, that corresponds to the calculation of efforts to
apply to the structure to take into account the contact, will be described in paragraph 3.

Therefore, our process includes first of all a stage of identification of potential contact
areas on the two structures in movement. Geometrical entities are then constructed with
nodes contained in each area. As part of the studies of contact between blades and casing,
the constructed entities will be therefore on the one hand a surface (the casing) and on
the other hand a skew curve (profile of the blade tip).

A research of intersection is then realised on these geometrical entities by following a
recursive process. Finally, once contact points are identified, it remains to project the
results on the initial mesh, that is to say to attribute to each impacted shell and to each
impacting node some values that will allow the calculation of the contact force.

Later on, we will suppose that the first stage (identification of potential contact areas)
has been undertaken. We therefore dispose of a group of nodes upon which we are going
to construct some geometrical entities presenting a sufficient character of continuity to
get rid off problems of facetisation. The flexibility given by the theory of splines is a great
argument that justify its use. This theory was mainly developed by de Boor in the late
70's ([1], and also [2]).

### 2.2.1 Modelisation of skew curves

Skew curves are represented using parametric coordinates: each point $P$ of the curve is
therefore determined by its three coordinates $x, y$ and $z$ that depend on a same parameter
$t$ whose interval of variation is to be defined with care [3]. The data set of the problem is
a group of $npt$ points $P_i$, and a degree $n$ of modelisation. The principle of construction
of a spline function is the following: within the interval of definition of the parameter $t$,
the vectorial expression of the skew curve approaching the group of points $P_i$ is:

$$c(t) = \sum_{i=0}^{N} Q_i B_{ni}(t) \qquad (1)$$

where $B_{ni}$ are basis functions of the space of the splines of degree $n$ called B-splines (these basis functions are in fact polynomials of degree $n$ defined within the parametric interval), and the $N+1$ points $Q_i$ are called control points of the curve. Three methods are frequently used to determine these control points:

* take for the $Q_i$ the points of the data set ($P_i$): it is called the direct method;

* calculate $Q_i$ in order that the resulting curve passes exactly through points $P_i$ for some particular values of the parameter $t$: it is therefore an interpolation;

* calculate $Q_i$ so that the sum of squares of distances from each point $P_i$ to the resulting curve is minimal: it is called the least square smoothing method.

### 2.2.2 Modelisation of surfaces

As for curves, surfaces will be represented using parametric coordinates, and control points can be determined with one of the three methods presented above. The data set of the problem is a group of $nptu \times nptv$ points $P_{ij}$. The parametrisation is therefore undertaken in two preferential directions we will note $u$ and $v$. The group of data points will be therefore approached by a surface whose vectorial equation is:

$$s(u, v) = \sum_{i=0}^{N_u} \sum_{j=0}^{N_v} Q_{ij} B_{nui}(u) B_{nvj}(v) \qquad (2)$$

where $Q_{ij}$ are the $(N_u + 1) \times (N_v + 1)$ control points of the surface. As shown in this formulation, the modelisation of a surface is equivalent to the crossed modelisation of several skew curves: those defined by points $Q_{ij}$, $j$ being constant, and those defined by points $Q_{ij}$, $i$ being constant. It is thus possible to use B-splines of different degrees in the two directions (these degrees are noted $n_u$ and $n_v$ in the expression above).

To take into account the thickness of the surface we translate the nodes of the mesh by a quantity equal to the half-thickness of the shells in the normal direction to elements. This normal direction is classically determined by calculating the average of normal vectors of all the elements that surround each node. Note that this operation requires the definition of an orientation of the surface (giving a point external to the surface for example).

## 2.3 Principle of the method

Now we have an analytic formulation of two geometrical elements, the problem to solve is then to determine their possible common parts. Methods commonly used in finite
element codes (LS-DYNA [4], Plexus [5]...) generally begin by dividing nodes into two
groups called “masters” and “slaves”. Then a first stage consists in undertaking a loop
on the totality or part of slave elements: for each of them, a test is realised on the totality
or on a part of the master elements so as to determine whether some of them can match
a predefined criterion. Afterwards, if some group of master/slave nodes have satisfied the
first criterion, a second test is realised to precisely confirm or invalidate the existence of
contact. The stage of search of master/slave nodes couples can be made in a global way
(hierarchical approach) or in a local way (vicinity approach).

In the suggested method, the first stage of classification is retained: the curve is chosen
to be the slave entity; the other stage is realised by using the spline formulation. Due
to the non - linearity of this formulation, an analytic solution can not be found in the
general case. A numerical iterative process has therefore to be set.

The result of the search of intersection is defined for a given precision $\varepsilon$. The result
given by the algorithm should be understood as follows: inside a sphere of radius $\varepsilon$, it
exists at least two points, one belonging to the curve and the other to the surface. Each
of the involved parameters, that is $t$ for the curve, and $(u, v)$ for the surface, is then
discretised. The obtained increments $dt$, $du$ and $dv$ depend on the precision $\varepsilon$ [6].

The boundaries of the parameterised areas can be determined after a first rough sorting
using a hierarchical approach for example. Subsequently, $P_t$ will designate the point of
the curve associated to the parameter $t$, and $P_{uv}$ the point of the surface associated to
the parameters $(u, v)$. The curve is used as a support for the search (see Figure 2).

<figure><img src="image_3.png" /><figcaption>Figure 2: Initial situation.</figcaption></figure>

For each discretised value of the parameter $t$, the surface is swept, and to each couple
$(u, v)$, the distance $P_tP_{uv}$ is calculated (see Figure 3). If this distance becomes smaller
than the precision $\varepsilon$, parameters $t$, $u$ and $v$ are stored in memory (these are told “initial
parameters”).

Values of parameters are then increased. As soon as the distance becomes bigger than $\varepsilon$, parameters are stored again in memory (these are told “terminal parameters”). Thus, some potential intersection areas between the two geometrical entities can be exhibited (see Figure 4). It is then possible to proceed to a new scanning of these areas, with a new precision smaller than $\varepsilon$. The operation is repeated until a required precision is reached.
<figure><img src="image_4.png" /><figcaption>Figure 3: Principle of the method.</figcaption></figure>

The first simulations realised with this method [7] show that the search algorithm is not so heavy in term of time calculation and underline the fact that the possibility to use meshes with fewer elements provides with a suitable payment for this increase of time calculation.

<figure><img src="image_5.png" /><figcaption>Figure 4: Final situation.</figcaption></figure>

# 3 EVALUATION OF THE CONTACT FORCE

## 3.1 Discrete contact

Unilateral contact is usually characterised by two parameters : gap between bodies candidates for contact $g$ and contact pressure $\lambda_n$. Three constraints derive from the unilateral contact law :

* the non-penetrability constraint $g \ge 0$,

* the compressive constraint $\lambda_n \le 0$, meaning that there is no tensile force between the bodies in contact,

* a complementary constraint $g\lambda_n = 0$, meaning that the bodies are either in contact ($g=0$) or separated ($\lambda_n = 0$).
As using the finite element method, we focus on discrete mechanical equations. In case of linear elastodynamics, without structural damping, the discrete matrix equation is :

$$M\ddot{U} + KU = F^{\text{ext}} - F^{\text{contact}} \quad (3)$$

where $M$ = mass matrix, $K$ = structural stiffness matrix, $U$ = nodal displacement vector, $\ddot{U}$ = nodal acceleration vector, $F^{\text{ext}}$ = generalised external force vector, $F^{\text{contact}}$ = generalised contact force vector.

The contact impenetrability constraint becomes $(G)_i \ge 0$ with the linearised gap vector :

$$G = G_0 + QU \quad (4)$$

where $(G)_i$ = gap vector component, $Q$ = constraint contact matrix, depending on the discrete contact surfaces.

The main mathematical models to compute contact forces are the Lagrange multiplier method and the penalty function method.

The Lagrange multiplier method enforces the exact linear impenetrability constraint but adds new unknowns. These unknowns, called Lagrange multipliers $\Lambda$, represent contact forces :

$$F^{\text{contact}} = Q^T \Lambda \quad (5)$$

The penalty function method allows penetration between contact surfaces by approximating the linear impenetrability constraint. The contact forces, also noted $\Lambda$, are proportional to penetration distances $G$, and represent reactions of fictitious springs located between contact entities (nodes or segments) :

$$F^{\text{contact}} = Q^T \Lambda \quad (6)$$

$$\Lambda = K_c G \quad (7)$$

where $K_c$ is the contact stiffness matrix. The main problem of the penalty function method is the contact stiffness choice : a small value induces large penetrations, whereas a large value induces numerical vibrations and some convergence problems. The most widely used contact stiffness expression is the one proposed by Hallquist et al. [8] :

$$k_c = f_s \frac{KA^2}{V} \quad (8)$$

in case of solid element interaction, with $K$, $A$ and $V$ the Bulk modulus, area and volume of the contact master element, respectively. The parameter $f_s$ is an arbitrary parameter, that user has to change according to contact problem.
## 3.2 Kinematic contact stiffness

The aim of our project was to develop a contact element in the finite element code MECANO (part of SAMCEF software) that uses the classical Newmark implicit time integration scheme :

$$U_{n+1} = U_n + \Delta t \dot{U}_n + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{U}_n + \beta \Delta t^2 \ddot{U}_{n+1} \quad (9)$$

$$\dot{U}_{n+1} = \dot{U}_n + (1 - \gamma) \Delta t \ddot{U}_n + \gamma \Delta t \ddot{U}_{n+1} \quad (10)$$

The Lagrange multiplier method is difficult to use when implemented in an implicit finite element code. It leads to numerical problems (i.e. undesirable oscillations and non convergence of the Newton-Raphson iterative process) due to the non-linearity of the contact phenomenon. Moreover, the results are depending on the Newmark parameter values [9], [10].

The penalty method has been chosen because it is more adaptable. But the tricky problem of the contact stiffness choice is still to be solved. A dynamic analysis of the Lagrange multiplier method, associated with a Newmark implicit scheme [9] and an explicit central difference scheme [11], has led to express a new contact stiffness. Indeed, in case of linear elastodynamics, using Newmark time integration scheme and considering a predictor (without contact) - corrector process, the predicted equilibrium state is :

$$\left(\frac{M}{\beta \Delta t^2} + K\right) \Delta U^* = F_{n+1}^{ext} - F_n^{int} + \frac{M}{\beta \Delta t^2} \left(\Delta t \dot{U}_n + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{U}_n\right) \quad (11)$$

where $\Delta U^* = \text{predicted displacement increment at } t_{n+1}$, $F_n^{int} = \text{vector of internal forces at time } t_n$.

Gap value is $G_{n+1}^* = G_0 + Q_{n+1}^* U_{n+1}^*$, where the matrix $Q_{n+1}^*$ represents the predicted constraints at $t_{n+1}$. We assume that all components of $G_{n+1}^*$ are negative, i.e. there are penetrations, and that $Q_{n+1} = Q_{n+1}^*$, i.e. the contact entities at $t_{n+1}$ are the ones determined during the predictor step. Next, during the corrector step, the velocity and displacement vectors are calculated in order to satisfy :

$$G_{n+1} = G_{n+1}^* + Q_{n+1}^* \Delta U^c = 0 \quad (12)$$

with $U_{n+1} = U_n + \Delta U^* + \Delta U^c$, the superscript c meaning corrected values.

New equilibrium state is expressed as :

$$M \ddot{U}_{n+1} + F_{n+1}^{int} = F_{n+1}^{ext} - Q_{n+1}^{T} \Lambda_{n+1} \quad (13)$$

and the contact force is :

$$\Lambda_{n+1} = \left[ Q_{n+1} \left( \frac{M}{\beta \Delta t^2} + K \right)^{-1} Q_{n+1}^{T} \right]^{-1} G_{n+1}^* \quad (14)$$
with a displacement correction $\Delta U^c = - \left( \frac{M}{\beta \Delta t^2} + K \right)^{-1} Q_{n+1}^T \Lambda_{n+1}$.

The “Forward Increment Lagrange Multiplier” proposed by Carpenter et al. [11] leads to a similar expression of the contact force vector with an explicit central difference scheme :

$$ \Lambda_n = [\Delta t^2 Q_{n+1} M^{-1} Q_{n+1}^T]^{-1} G_{n+1}^* \quad (15) $$

For a contact correction at the same time $t_{n+1}$, implicit and explicit contact force vectors (Relations 14 and 15) can be regarded as the product of an equivalent contact stiffness and the predicted gap $G_{n+1}^*$. Assuming a Lagrange multiplier method behaviour similar to the penalty function method one, it is possible to conclude that the Lagrange multiplier convergence problems are due to a too large value of its equivalent contact stiffness. Therefore a new contact stiffness is proposed by truncating equivalent contact stiffness derived from the contact force expression (Relation 14) :

$$ K_c = \frac{1}{\beta \Delta t^2} [Q M^{-1} Q^T]^{-1} \quad (16) $$

The “influence contact matrix” $[QM^{-1}Q^T]$ is a condensation result of all contact candidate masses. This new contact stiffness is then called “kinematic contact stiffness”.

Kinematic contact stiffness, unlike Hallquist one [8], does not depend on material features. Furthermore, the Newmark parameter $\beta$, that takes part in time integration accuracy, has a new function, i.e. a contact hardness rate.

A simple analysis of two punctual mass impact, initially proposed by Chaudhary and Bathe [12], enables to deduce a shock hardness rate associated to kinematic contact stiffness [13] :

$$ 4\beta - \gamma - 1/2 \quad (17) $$

This rate depends on the two Newmark parameters. If it is positive, then contact releases too late, and if it is negative, contact releases earlier. A particular point of interest is that the trapezoidal rule, $\beta = 1/4; \gamma = 1/2$, leading to the most accurate time integration, gives a rate equal to zero. A similar study has been done with the Lagrange multiplier method. It has led to another shock hardness rate expression [13] : $2\beta - \gamma - 1/2$. Values recommended by Chaudhary and Bathe [12] ($\beta = \gamma = 1/2$) give a rate equal to zero, i.e. contact is well released.

# 4 IMPLEMENTATION AND RESULTS

The new contact element is developed in the user Fortran routine of the finite element software MECANO. It consists of a contact analytical detection phase, described in Section 2, and a contact handling with penalty function method associated with kinematic contact stiffness, described in Section 4.
## 4.1 Kinematic contact stiffness assessment

<figure><img src="image_6.png" /><figcaption>Figure 5: Non dimensional impact of two elastic rods.</figcaption></figure>

The one-dimensional impact of two different elastic rods [15] is used to analyse the kinematic contact stiffness behaviour. The impact features (without dimension) are listed in Table 1.

<table><caption>Table 1: Different rod impact features.</caption><thead><tr><th></th><th>rod 1</th><th>rod 2</th></tr></thead><tbody><tr><td>Young's modulus</td><td><i>E</i><sub>1</sub> = 0.49</td><td><i>E</i><sub>2</sub> = 1.0</td></tr><tr><td>density</td><td><i>ρ</i> = 1.0</td><td><i>ρ</i> = 1.0</td></tr><tr><td>wave speed</td><td><i>c</i><sub>1</sub> = 0.7</td><td><i>c</i><sub>2</sub> = 1.0</td></tr><tr><td>cross-section area</td><td><i>S</i> = 1.0</td><td><i>S</i> = 1.0</td></tr><tr><td>length</td><td><i>L</i> = 10.</td><td><i>L</i> = 10.</td></tr><tr><td>finite element number</td><td><i>n</i><sub>1</sub> = 100</td><td><i>n</i><sub>2</sub> = 70</td></tr><tr><td>initial velocity</td><td><i>V</i><sub>0</sub><sup>1</sup> = 0.1</td><td><i>V</i><sub>0</sub><sup>2</sup> = 0.0</td></tr></tbody></table>

Comparisons are made with the exact solution and the Lagrange multiplier method. Furthermore different Newmark parameter values are used :

* the trapezoidal rule, $\beta = 0.25$; $\gamma = 0.50$, leading to a constant acceleration approximation and satisfying the kinematic stiffness shock rate $4\beta - \gamma - 1/2 = 0$;

* a special couple, $\beta = 0.505$; $\gamma = 0.51$, satisfying the Lagrange multiplier shock rate $2\beta - \gamma - 1/2 = 0$.

Figures 6 and 7 present displacement evolution of the contact rod surfaces, Figures 8 and 9 present velocity evolution of the contact rod surfaces and Figures 10 and 11 present contact force evolution.

Using the new contact stiffness with a penalty function method in an implicit finite element code, has several advantages over the classical Lagrange multiplier method :

- a decrease of Newton-Raphson iteration number (more than 12%, see Table 2);

- fewer convergence problems ;
<figure><img src="image_7.png" /><figcaption>Figure 6: Contact surface displacement evolution, with trapezoidal rule ($\beta = 1/4$; $\gamma = 1/2$).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Contact surface displacement evolution, with Newmark parameters satisfying $(2\beta - \gamma - 1/2) = 0$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Contact surface velocity evolution, with trapezoidal rule ($\beta = 1/4$; $\gamma = 1/2$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Contact surface velocity evolution, with Newmark parameters satisfying $(2\beta - \gamma - 1/2) = 0$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Contact force between elastic rods, with trapezoidal rule ($\beta = 1/4$; $\gamma = 1/2$).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Contact force between elastic rods, with Newmark parameters satisfying $(2\beta - \gamma - 1/2) = 0$.</figcaption></figure>
<table><caption>Table 2: Performance of two different rod impact.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>(β = 0.25; γ = 0.50)</td><td>Kinematic stiffness</td><td>Lagrange multiplier</td></tr><tr><td>average iteration number</td><td>1.05</td><td>1.20</td></tr><tr><td>maxi. iter. number for one step</td><td>2</td><td>2</td></tr><tr><td>total iteration number</td><td>316</td><td>358</td></tr></tbody><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>(β = 0.505; γ = 0.51)</td><td>Kinematic stiffness</td><td>Lagrange multiplier</td></tr><tr><td>average iteration number</td><td>1.01</td><td>1.17</td></tr><tr><td>maxi. iter. number for one step</td><td>2</td><td>2</td></tr><tr><td>total iteration number</td><td>302</td><td>352</td></tr></tbody></table>

- a compatibility with all classical Newmark parameter values. Especially, the trapezoidal rule leads to good results, with fewer oscillations. Values satisfying $2β - γ - 1/2 = 0$, leading to better results with the Lagrange multiplier method than the trapezoidal rule [13], give very good results with the new contact element. Indeed, these values produce an useful damping effect on the numerous oscillations in contact forces and contact entity velocities.

## 4.2 Contact detection

The method of contact search with spline functions has been implemented in the new contact element.

We intend to compare the results given by a test involving contact when different methods are used. The test we use simulates the rotation of an unbalanced blade inside a casing. The blade is represented by a generalised spring, and the casing (part of a cylinder clamped at its basis) is modelised by shells of Mindlin. The new contact element is used to take into account contact phenomena between the blade tip and the casing. The mesh is represented in Figure 12.

The first method which is implemented will be designated under the term “classic method”: the search of contact uses geometrical entities of degree 1 built on finite elements. These entities are parts of the shells approximating plans for the casing, and segments for the blades. The other method we implement is the spline method in two formulations: interpolation and smoothing.

The test is used to observe the value of the normal force on the blade tip during the contact, and the CPU time to simulate 5 rotations of the blade. The varying parameters are on the one hand the choice of the method of search of the contact (classic, interpolation splines and smoothing splines) and on the other hand the number of elements used to discretise the circumference of the casing (three configurations have been studied: 72, 36 and 18 elements).
<figure><img src="image_13.png" /><figcaption>Figure 12: General view of the mesh.</figcaption></figure>

### 4.2.1 Contact force

The first presented results concern the contact force : these results are normalised with respect to the maximal amplitude of the first contact force obtained with a method of interpolation in the 72 element configuration. It can be observed in Figure 13 (72 elements in the circumference of the casing) that spline methods give very similar results with regard to the maximal amplitude of the force of the first contact (the difference between the methods is approximately 3%).

<figure><img src="image_14.png" /><figcaption>Figure 13: Configuration with 72 elements.</figcaption></figure>
With a classic method, the radius of the casing is systematically underestimated: as a consequence the contact force whose calculation depends on the radial penetration is always overestimated. In the present case, the contact force is about 30% bigger than the force given by a spline method. Moreover, the finite element approximation involves an early detection of the instant of the first contact (t=4 ms for the classic method and t=6 ms with a spline interpolation).

When the number of elements used to define the circumference of the casing decreases, these two problems greatly emphasize. When just 18 elements are used for example (see Figure 14), the classic method detects a strong initial penetration that induces an amplitude of the first contact force twice bigger than the one obtained with a spline interpolation.

<figure><img src="image_15.png" /><figcaption>Figure 14: Configuration with 18 elements.</figcaption></figure>

By representing on a same figure the results given by the three configurations (72, 36 and then 18 elements in the circumference of the casing), it is interesting to note that the interpolation spline method is very little sensitive to the fineness of the mesh concerning the amplitude of the first contact force peaks (see Figure 15).

### 4.2.2 CPU time

The second series of results focuses on the required CPU time to simulate 5 revolutions of the blade. All CPU times presented on Figure 16 are normalised with respect to the required CPU time with 72 elements and the classic method. It appears that spline methods are more expensive than the classic method (approximately 25% of supplementary time in the configuration with 72 elements). However, it is important to note that a spline
<figure><img src="image_16.png" /><figcaption>Figure 15: Interpolation method.</figcaption></figure>

interpolation method with 36 elements requires less CPU time (approximately 30%) than
the classic method with 72 elements, while providing a comparable result.

<figure><img src="image_17.png" /><figcaption>Figure 16: CPU time comparison.</figcaption></figure>
# 5 CONCLUSION

Two methods have been developed to deal with contact phenomena within an implicit finite element code. The first one uses a geometrical approach to precisely determine the intersection points using a continuous representation of the structure which allows to get rid off facetisation problems. The second one consists in a simplified contact element using a penalty function method, and based on the definition of a new kinematic contact stiffness that provides better results than the classical Lagrange multiplier method.

These methods have been successfully implemented in the finite element code Samcef, and used on a simplified test. Future works will be dedicated to the simulation of more complex configurations, including some validating tests given by the literature.
## REFERENCES

[1] C. de Boor, *A practical guide to Splines*, vol. 27 of Applied Mathematical Sciences, Springer Verlag, (1978).

[2] G. Nurnberger, *Approximation by spline functions*, Springer Verlag, (1989).

[3] M. Daniel, *Modélisation de courbes et surfaces par des B-splines - Application à la conception et à la visualisation de formes*, PhD Thesis (in French), Université de Nantes (1989).

[4] J.O. Hallquist, *LS-Dyna3D, Theoritical Manual*, Dynalis, (1995).

[5] CEA, *Plexus, Theoritical Manual*, CEA, (1997).

[6] E. Arnoult, *Modélisation numérique et approche expérimentale du contact en dynamique - Application au contact aubes/carter de turboréacteur*, PhD Thesis (in French), Université de Nantes (2000).

[7] E. Arnoult, B. Peseux and J. Bonini, "Recherche analytique du contact appliquée à un code éléments finis en dynamique", *Proceedings of the 14th French Congress of Mechanics AUM'99, Toulouse Aug. 30th - Sept. 3rd 1999, CD-ROM 218*, 1-6 (1999).

[8] J.O. Hallquist, G.L. Goudreau and D.J. Benson, "Sliding interfaces with contact-impact in large-scale lagrangian computations", *Comp. Meth. in Appl. Mech. and Engng.*, **51**, 107-137 (1985).

[9] I. Guilloteau, B. Peseux and J. Bonini, "Simplified dynamic contact model – Application to the rotor/stator interaction", *Proceedings of the European Conference on Computational Mechanics ECCM'99, CD-ROM 592*, 1-13 (1999).

[10] I. Guilloteau, B. Peseux and J. Bonini, "Modélisation du contact en implicite - Application à l'interaction rotor/stator" (in French), *Revue Européenne des Eléments Finis*, to be published (2000).

[11] N.J. Carpenter and R.L. Taylor and M.G. Katona, "Lagrange constraints for transient finite element surface contact", *Int. J. Num. Meth. Engng.*, **32**, 103-128 (1991).

[12] A.B. Chaudhary and K.J. Bathe, "A solution method for static and dynamic analysis of three-dimensional contact problems with friction", *Comp. & Struct.*, **24**, 855-873 (1986).

[13] I. Guilloteau, *Modélisation du contact en implicite - Application à l'interaction rotor/stator*, PhD Thesis (in French), Université de Nantes (1999).
[14] M. Géradin and D. Rixen, *Théorie des vibrations - Application à la dynamique des structures*, 2e édition, Masson, Paris (1996).

[15] T.J.R. Hughes, R.L. Taylor, J.L. Sackman, A. Curnier and W. Kanoknukulchai, A finite element method for a class of contact-impact problems, *Comp. Meth. in Appl. Mech. and Engng.*, **8**, 249-276 (1976).

[16] SAMCEF, *MECANO/ROTORT Manual, V7.1.*, SAMTECH S.A. (1997).