# Validation of a 3D contact algorithm for the study of blade-tip/casing contacts in turbomachines

Alain Batailly, Benoit Magnain, Mathias Legrand, Christophe Pierre

► To cite this version:

Alain Batailly, Benoit Magnain, Mathias Legrand, Christophe Pierre. Validation of a 3D contact algorithm for the study of blade-tip/casing contacts in turbomachines. 8th IFToMM International Conference on Rotordynamics, Sep 2010, Séoul, South Korea. hal-00533382

HAL Id: hal-00533382

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

Submitted on 5 Nov 2010

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
# Validation d'un algorithme de contact 3D utilisé pour l'étude de contacts aube/carter dans les turbomachines aéronautiques

**Alain Batailly**

Post-doctorant
Université McGill
Bâtiment d'ingénierie McDonald
Montréal H3A 2K6, Canada
alain.batailly@mcgill.ca

**Benoît Magnain**

Maître de conférences
ENSIB-Institut PRISME
88 boulevard lahitolle
18020 Bourges cedex, France
benoit.magnain@ensi-bourges.fr

**Mathias Legrand**

Post-doctorant
Université McGill
Bâtiment d'ingénierie McDonald
Montréal H3A 2K6, Canada
mathias.legrand@mcgill.ca

**Christophe Pierre**

Professeur
Université McGill
Bâtiment d'ingénierie McDonald
Montréal H3A 2K6, Canada
christophe.pierre@mcgill.ca

**Note 1** : cet article a été publié dans sa version anglaise parmi les actes de la conférence IFTOMM 2010.

**Note 2** : la version française de cet article est abrégée. Les figures et équations sont présentes dans la version anglaise uniquement.

## 1 Résumé

Afin d'optimiser le rendement des turbomachines aéronautiques, les concepteurs cherchent à réduire les pertes aérodynamiques entre les différents étages de la structure. Ces pertes sont directement liées aux jeux entre les parties tournantes et fixes qui sont donc réduits au minimum. La contre-partie de cette optimisation réside dans le fait que les contacts entre parties tournantes (rotor) et statiques (stator) deviennent plus fréquents. Parmi les phénomènes d'interaction rotor/stator pouvant être créés par ces contact, citons le problème d'interaction modale [1] dont les conséquences peuvent être critiques pour l'intégrité de la turbomachine.
Ce papier présente la validation d'un algorithme de contact 3D utilisé pour la simulation de telles interactions rotor/stator dans une configuration non-accidentelle. Cet algorithme de contact repose sur un schéma temporel explicite (différences finies centrées) et l'utilisation de la méthode des multiplicateurs de Lagrange pour le traitement du contact, qui permet notamment de respecter les conditions de non-pénétrabilité. Cet algorithme de contact est validé par comparaison à deux autres codes de contact :

* un algorithme avec schéma d'intégration temporelle implicite et traitant le contact à l'aide de la méthode du bi-potentiel;

* un code industriel, ANSYS.

Les résultats obtenus mettent en évidence la qualité de l'algorithme proposé dont il est égale- ment montré qu'il permet la conservation de l'énergie totale du système. À la fin de l'étude, un cas d'application est détaillé : l'algorithme de contact est utilisé pour simuler un cas de contact aube/carter avec des modèles éléments finis 3D.

## 2 Introduction

La simulation des problèmes de contact a donné lieu à de nombreuses études dans le but d'améliorer la précision et la stabilité des solutions obtenues lorsque des modèles éléments finis de très grande taille sont utilisés. D'un point de vue mathématique, un problème de contact est intrinsèquement non linéaire puisque la surface de contact est inconnue *a priori* et que et que les conditions limites sont déterminées à partir de la solution calculée. On distingue généralement trois types de méthodes pour le traitement du contact :

1. la méthode de pénalisation [2];

2. la méthode des multiplicateurs de Lagrange [3];

3. la méthode du Lagrangien augmenté [4, 5].

Une différence essentielle entre ces méthodes réside dans le fait que la méthode de pénalisa- tion repose sur le choix d'un paramètre numérique et autorise un certain niveau de pénétration alors que les multiplicateurs de Lagrange (qui représentent les efforts de contact) sont précisé- ment calculés pour annuler les pénétrations détectées. Parmi les stratégies admissibles, une formulation reposant sur la méthode du Lagrangien augmenté combinée à la méthode du bi- potentiel est proposée dans [6]. Avec cette formulation, le problème de contact avec frottement est résolu dans un espace réduit avec un algorithme de résolution de type prédiction/correction.
De plus, la méthode du bi-potentiel permet d'obtenir une seule inégalité couplant le contact unilatéral et le frottement à l'aide d'un bi-potentiel de contact. Dans notre étude, nous nous intéressons aux problèmes de contact aube/carter dans les turbomachines. Ces contacts sont la contre-partie inévitable de la réduction des jeux entre parties tournantes et fixes des turbomachines. En plus des dégâts qu'ils peuvent engendrer [7, 8], ces contacts peuvent provoquer l'apparition du phénomène d'interaction modale [1]. Des études précédentes réalisées sur des modèles simplifiés 2D [9] ont permis de mettre en évidence la grande sensibilité numérique des simulations d'interaction modale. Par conséquent, une très grande précision est requise pour le traitement du contact ce qui est incompatible avec l'utilisation d'une méthode de pénalisation. L'algorithme de contact présenté dans cette étude repose sur la méthode des multiplicateurs de Lagrange combinée à un schéma d'intégration temporelle explicite. Déjà utilisé avec des modèles simplifiés 2D [9], cet algorithme est ici détaillé et validé par comparaison avec deux autres codes. Un intérêt particulier est porté à la conservation de l'énergie au cours des simulations.

# 3 Éléments théoriques

## 3.1 Mécanique du contact

Les détails relatifs au traitement du contact, à l'utilisation des surfaces B-splines et à l'algorithme de résolution sont très brièvement décrits dans cette section. En particulier, les équations présentes dans la version anglaise ne sont pas réécrites mais il y est fait référence.

Le problème à résoudre est le système d'équations (1). Pour la résolution d'un problème de contact, il peut être utile de définir arbitrairement deux surfaces : l'une dite "maître" et l'autre "esclave". À chaque point matériel de la surface maître on associe le point le plus proche de la surface esclave en définissant la fonction $\bar{y}$ telle que décrite par l'équation (2). Il est alors possible d'expliciter la fonction écart (équation (3)) entre les deux structures.

Les conditions de contact utilisées sont les conditions dites de Signorini, détaillées dans l'équation (4) auxquelles est ajoutée une loi de frottement de Coulomb.

## 3.2 Surfaces B-splines

Dans notre étude la surface esclave est le carter. Afin de simuler avec une grande précision les contacts aube/carter, il est impératif d'éviter tout problème de discontinuité au passage d'un élément fini à un autre sur cette surface. Par conséquent, il est proposé de lisser la surface esclave à l'aide de surfaces B-splines.

Les surfaces B-splines sont usuellement définies comme le produit tensoriel de deux courbes
B-splines (équation (6)) et dépendent ainsi de deux paramètres *s* et *t* et d'un ensemble de points de contrôle $Q_{ij}$ qui forment un maillage sur la circonférence du carter.

## 3.3 Algorithme de résolution

La méthode utilisée dans notre étude [10] est l'algorithm de Carpenter qui repose sur la méthode des multiplicateurs de Lagrange. Un inconvénient mineur de cet algorithme, le décalage d'un pas de temps entre l'estimation de l'accélération et la prise en compte des efforts de contact, est contrebalancé par sa simplicité, le respect parfait des conditions de non-pénétration ainsi que la compatibilité des vitesses et des accélérations. Les trois étapes de cet algorithme sont :

1. **prédiction** : les déplacements au pas de temps *n* + 1 sont calculés sans considérer le contact (équation (9));

2. **calcul des pénétrations** : détermination de la fonction écart entre les deux structures;

3. **correction** : les déplacements sont corrigés grâce au calcul des multiplicateurs de Lagrange (11) et annulation de la fonction écart.

# 4 Validation

Les résultats fournis par le code de contact proposé dans notre étude sont confrontés à ceux donnés par un algorithme basé sur la méthode du bi-potentiel [6] avec un schéma d'intégration temporelle explicite et par le logiciel Ansys. Comme mentionné en introduction, le code de contact présenté est développé dans le but de simuler des cas d'interaction aube/carter. Du fait du nombre très élevé de degrés de liberté dans les modèles éléments finis des structures associées (une roue aubagée compte approximativement un million de degrés de liberté), il n'est pas possible d'utiliser un cas de contact aube/carter comme cas de validation. Par conséquent, le cas de validation considéré dans notre étude (représenté sur la figure 1) est un cas simplifié qui est comparable à un cas de contact aube/carter du fait des vitesses relatives et du frottement considérés.

Un cube est propulsé entre deux anneaux flexibles encastrés sur leurs extrémités. L'espace entre les anneaux étant inférieur à la largeur du cube, ce dernier entre en contact avec les deux anneaux qui sont repoussés. Le cas de contact étant symétrique par rapport au plan $(O; \vec{x}; \vec{y})$, une première validation du code proposé consiste à s'assurer de la symétrie des résultats obtenus, ce qui est mis en evidence sur les figures 3(a) et 3(b).
Les paramètres de simulation sont les suivants : le module d'Young du cube est cent fois plus faible que le module d'Young des anneaux ($E_c = 210 \cdot 10^7$ Pa et $E_a = 210 \cdot 10^9$ Pa), les masses volumiques et les coefficients de Poisson des deux structures sont identiques $\rho_c = \rho_a = 7800 \text{ kg} \cdot \text{m}^{-3}$ et $v_c = v_a = 0,3$, le coefficient de frottement vaut $\mu = 0,15$ et la vitesse initiale<fnref n="1" /> du cube vaut $v_c = 100 \text{ m} \cdot \text{s}^{-1}$.

Une comparaison des résultats obtenus entre les différents codes utilisés est présentée sur les figures 4(a), 4(b), 4(c), 4(d) et 5. Ces courbes montrent la bonne correspondance des résultats pour les déplacements et les efforts de contacts observés.

Enfin, il est montré sur un cas test similaire mais sans frottement que l'énergie totale du système est conservée au cours du temps (voir figures 6(a) et 6(b)).

# 5 Application

Le cas d'application proposé est un cas d'interaction aube/carter avec une roue aubagée de 56 aubes partiellement dessinée sur la figure 7 et un carter rigide sur lequel est construit la surface B-spline. Du fait de la taille du modèle éléments finis de la roue, une méthode de synthèse modale est utilisée pour réduire le nombre de degrés de liberté. Par ailleurs, certaines hypothèses simplificatrices sont faites :

1. pour chaque aube, trois nœuds sont utilisés pour la gestion du contact (sur le bord d'attaque, le milieu de corde et le bord de fuite);

2. il est supposé que la composante non-radiale des efforts de contact en sommet d'aube est uniquement tangentielle (composante axiale nulle);

3. la déformation du carter se fait selon un mode à deux diamètres (ce qui correspond à une ovalisation de la structure).

Après étude de convergence (figures 9 et 10) des modèles réduits calculés, les résultats temporels d'une simulation de contact sont tracés sur la figure 11.

Enfin, une analyse de l'influence des diamètres nodaux dans la déformation de la roue aubagée montre la prédominance du diamètre 2 par rapport aux autres diamètres, en accord avec le cas de contact considéré.

<footnote n="1">¹La vitesse initiale est supposée suffisamment élevée pour pouvoir négliger les effets de la gravité.</footnote>
# 6 Conclusion

Cette étude a permis de présenter et valider un algorithme de contact 3D conçu pour la simulation de cas d'interaction aube/carter. La validation a permis de montrer que cet algorithme est conservatif et le cas d'application a mis en évidence sa compatibilité avec les méthodes de synthèse modale. Il est ainsi adapté à la simulation de contact entre modèles comportant un très grand nombre de degrés de liberté. Les résultats présentés permettent d'envisager une très large gamme d'applications et de possibilités d'études pour mieux simuler les interactions aube/carter. Des développements sont en cours pour intégrer la modélisation d'un revêtement abradable sur la surface du carter.
# Validation of a 3d contact algorithm for the study of blade-tip/casing contacts in turbomachines

**Alain Batailly**

Post-doctoral fellow
McGill University
McDonald Engineering bldg
Montréal H3A 2K6, Canada
alain.batailly@mcgill.ca

**Benoît Magnain**

Associate Professor
ENSIB-Institut PRISME
88 boulevard lahitolle
18020 Bourges cedex, France
benoit.magnain@ensi-bourges.fr

**Mathias Legrand**

Post-doctoral fellow
McGill University
McDonald Engineering bldg
Montréal H3A 2K6, Canada
mathias.legrand@mcgill.ca

**Christophe Pierre**

Professor
McGill University
McDonald Engineering bldg
Montréal H3A 2K6, Canada
christophe.pierre@mcgill.ca

## 1 ABSTRACT

The optimization of modern turbomachines through the diminution of aerodynamic losses leads to more frequent contacts between the fixed and rotating components. Modal interaction is one of the phenomena that may arise from these contacts with potential dangerous consequences for the integrity of the turbomachine. This paper presents the validation of a 3D contact algorithm – within the context of the finite element method – allowing to efficiently study the phenomenon of interaction. It is based on the forward increment Lagrange multiplier method which allows the combination of explicit time integration with the respect of non penetration conditions. This contact algorithm is validated with respect to another algorithm based on an implicit time integration procedure and the commercial software Ansys. It is highlighted in this paper that results obtained with the proposed algorithms match reference results. Also, results emphasize that the presented algorithm allows for energy conservation. This contact algorithm is then used in combination with reduction techniques to simulate a blade-tip/casing contact case with 3D industrial finite element models.
## 2 INTRODUCTION

Contact problems have been extensively studied in the past in order to improve the accuracy and the stability of the computed solutions while dealing with larger finite element models. From the mathematical point of view, general contact problems are inherently nonlinear since the contact area is a priori unknown and the boundary conditions are determined as part of the solution. The main mathematical difficulty lies in the natural first order constitutive laws of contact and friction phenomena expressed by non-smooth multivalued force-displacement or force-velocity relations. Several numerical techniques in the finite element framework were developed for contact simulations. In general, three groups of methods are considered: (1) the penalty method [2], (2) the Lagrange multiplier method [3] and (3) the augmented Lagrangian method [4, 5]. A major difference between these methods is that while the penalty method relies on a choice of a parameter that allows penetrations between the structures, the Lagrange multipliers – which represent the contact reaction forces – are introduced to enforce exactly the non penetration conditions.

Among the possible strategies, an augmented Lagrangian formulation combined with the bi-potential method was proposed in [6]. In this formulation, the frictional contact problem is treated in a reduced system by means of a predictor-corrector solution algorithm. In addition, the bi-potential method leads to a single displacement variational principle and a unique inequality in which unilateral contact and friction are coupled via a contact bi-potential.

In our study, the focus is made on blade-tip/casing contact simulations in modern turbomachines. These contacts are the unavoidable counterpart of the optimization of the energy efficiency of modern aircraft engines leading to the reduction of structural clearances between the rotor and the stationary casing components. Beside the severe damages they may cause [7, 8], these contacts may be responsible for the modal interaction phenomenon [1] to arise. Previous studies conducted on 2D planar models [9] show that the simulation of blade-tip/casing interaction is very sensitive on a numerical point of view and requires a high precision on contact management which is incompatible with the use of the penalty method.

An algorithm based on the forward increment Lagrange multiplier method is presented in our study. This method combines the use of explicit integration operators while satisfying the non penetration constraints. Previously introduced for the study of 2D planar models [9], it is here presented for 3D contact configurations. In the first section of this paper, the contact algorithm is presented and validated considering other codes and focusing on the energy conservation aspect. Then, an application case with 3D industrial models is given in the last section.
## 3 NOMENCLATURE

**M** mass matrix

**K** stiffness matrix

**D** damping matrix

**F** external efforts

**u** solution vector

**h** time step

**E** Young modulus

**ρ** density

**ν** Poisson coefficient

**μ** friction coefficient

**λ** Lagrange multiplier

## 4 CONTACT ALGORITHM

### 4.1 General theory

The forces of particular interest in this study are the forces of contact acting between the blade tips and the casing. Equations of motion are derived using the Principle of Virtual Work within the kinematically linear framework following, in essence, the procedure described in [11]. After discretization of the weak formulation of the equations of motion, the general problem to be solved is stated as follows:

$$ egin{cases} 
\mathbf{M}\ddot{\mathbf{u}} + \mathbf{D}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{F}_{\text{ext}} + \mathbf{F}_{\text{c}} \\ 
\mathbf{u}(t = t_0) = \mathbf{u}_0 \\ 
\dot{\mathbf{u}}(t = t_0) = \dot{\mathbf{u}}_0 
\end{cases} 
\quad (1) $$

It is first convenient to arbitrary choose one surface subjected to contact as the master one so that the second one, commonly called the slave surface, can be parameterized. It is then possible to find for any material point **x** belonging to the master surface $\Gamma_c^{(m)}$, its closest counterpart $\tilde{\mathbf{y}}$ on the slave surface $\Gamma_c^{(s)}$:

$$ \tilde{\mathbf{y}} = \arg\min_{\mathbf{y} \in \Gamma_c^{(s)}} \|\mathbf{x} - \mathbf{y}\| \quad (2) $$

According to these notations, the gap function between the two structures can be stated:

$$ g(\mathbf{x}) = g_0(\mathbf{x}) + (\mathbf{u}^{(m)}(\mathbf{x}) - \mathbf{u}^{(s)}(\tilde{\mathbf{y}}(\mathbf{x}))) \cdot \mathbf{n} \quad (3) $$
where $g_0(\mathbf{x})$ represents the initial positive gap between the two structures, $\mathbf{n}$, the outward normal to $\Gamma_c^{(s)}$, and $\mathbf{u}$, the displacement of the master $^{(m)}$ and slave $^{(s)}$ structures. The contact conditions, referred to as the Kuhn-Tucker optimality conditions in the parlance commonly used in the literature, are such that for all $\mathbf{x} \in \Gamma_c^{(m)}$:

$$t_N \ge 0, g \ge 0, t_N g = 0 \qquad (4)$$

where $t_N$ stands for the contact pressure, assumed positive, acting on $\Gamma_c^{(m)}$. To these unilateral contact conditions, we add the well-known Coulomb friction law:

$$ \begin{aligned} \|\mathbf{t}_T\| &\le \mu t_N \\ \|\mathbf{t}_T\| &< \mu t_N \Rightarrow \mathbf{v}_T = 0 \\ \|\mathbf{t}_T\| &= \mu t_N \Rightarrow \exists \alpha > 0 \text{ such as } \mathbf{v}_T = \alpha \frac{\mathbf{t}_T}{\|\mathbf{t}_T\|} \end{aligned} \qquad (5) $$

for which $\mu$ is the coefficient of friction, $\mathbf{v}_T$, the tangential slip and $\mathbf{t}_T$, the tangential stress vector.

The virtual work of the contact forces is calculated on the slave surface $\Gamma_c^{(s)}$, in other words, the surface spline attached to the casing described in the next section. It means that the gap functions are evaluated using the surface spline -and not the original casing mesh- and that the virtual displacements of the interface nodes due to the contact forces acting on the surface spline are calculated using the shape functions of the surface spline together.

## 4.2 Bicubic uniform B-spline surfaces

In classical contact algorithms, a master-slave approach is used such that, in our study, the contacting master nodes of the bladed disk are not allowed to penetrate the bilinear surface facets formed by the slave elements of the casing. Problems occur when nodes must slide from one facet to an adjacent one due to the rotation of the engine: jumps in contact forces are typically encountered. These jumps are non-physical and can cause serious errors in the resulting simulated stresses and forces.

Consequently, the use of smoothing methods is required. Several of them have been developed in the past few years [12, 13]. They increase the convergence rate of the contact algorithms by removing the facetization issues while preserving the original meshes. In the present study, such a smoothing procedure is performed through the use of a bicubic B-spline surface.

B-spline surfaces are an extension of B-spline curves and are most commonly defined as the tensor product of B-spline curves:

$$ \mathbf{S}(s, t) = \sum_{i=0}^{M-1} \sum_{j=0}^{N-1} \mathbf{Q}_{ij} \mathbf{B}_{ni}(s) \mathbf{B}_{nj}(t) \qquad (6) $$
involving two parameters s and t defined on two knot vectors ($s_i$) and ($t_j$), and a set of control points $Q_{ij}$, organized in a grid and thus, not of arbitrary topology. Accordingly, the surface basis functions are products of two univariate bases. Similar to the uniform cubic B-spline curve, $S(s, t)$ is a vector-valued function and each patch<fnref n="2" /> $pq$ of a bicubic uniform B-spline surface can be described in a matrix form by:

$$S_{pq}(s, t) = \mathbf{S} \mathbf{H} Q_{pq} \mathbf{H}^{\mathrm{T}} \mathbf{T}^{\mathrm{T}} \quad (7)$$

with $\mathbf{S} = (1, s, s^2, s^3)$, $\mathbf{T} = (1, t, t^2, t^3)$ and:

$$\mathbf{Q}_{pq} = \begin{bmatrix} Q_{p-1,q-1} & Q_{p-1,q} & Q_{p-1,q+1} & Q_{p-1,q+2} \\ Q_{p,q-1} & Q_{p,q} & Q_{p,q+1} & Q_{p,q+2} \\ Q_{p+1,q-1} & Q_{p+1,q} & Q_{p+1,q+1} & Q_{p+1,q+2} \\ Q_{p+2,q-1} & Q_{p+2,q} & Q_{p+2,q+1} & Q_{p+2,q+2} \end{bmatrix}$$

and matrix $\mathbf{H}$ is defined as:

$$\mathbf{H} = \frac{1}{6} \begin{bmatrix} -1 & 3 & -3 & 1 \\ 3 & -6 & 3 & 0 \\ -3 & 0 & 3 & 0 \\ 1 & 4 & 1 & 0 \end{bmatrix} \quad (8)$$

As shown in equation (7), each bicubic B-spline patch is affected locally by sixteen control points.

## 4.3 Solution algorithm for interaction

Because of its simplicity, the Forward Increment Lagrange Method [10] is used in our study. The algorithm is divided into three steps:

1. **prediction** at time step $n+1$ of the displacements $\mathbf{u}^3$ of the structures without considering any contact. This predicted displacement, denoted with a superscript $p$, is analytically expressed as:

$$\mathbf{u}^{n+1,p} = \left[ \frac{\mathbf{M}}{h^2} + \frac{\mathbf{D}}{2h} \right]^{-1} \left( \left( \frac{2\mathbf{M}}{h^2} - \mathbf{K} \right) \mathbf{u}^n + \left( \frac{\mathbf{D}}{2h} - \frac{\mathbf{M}}{h^2} \right) \mathbf{u}^{n-1} + \mathbf{F}^n \right) \quad (9)$$

where $h$ is the time-step size of the explicit time marching procedure. Previous displacements $\mathbf{u}^n$ and $\mathbf{u}^{n-1}$ and external forcing $\mathbf{F}^n$ are known.

<footnote n="2">²two-dimensional counterpart of a segment, written here in one direction of the Cartesian space</footnote>

<footnote n="3">³in $\mathbf{u}$ are stored the displacements of both the casing and the bladed disk in their reduced-order Craig-Bampton space. Subsequent mass, stiffness, and damping matrices are written in a consistent manner</footnote>
2. **determination** of the gap function vector $g^p$ between the two structures using equation (3). Each gap function where a penetration has been detected is kept in $g^p$, all other coordinates of the latter being zero.

3. **correction** of the displacements through the calculation of the Lagrange multipliers. This step implies that the gap functions (linearized when necessary) vanish:

$$g^{n+1} = C_N^t u^{n+1,c} + g^p = 0 \qquad (10)$$

where the superscript c means that the correction of the displacements is being calculated. $C_N$ is the contact constraint matrix in the normal direction where superscript t stands for *transpose*. The new equations of motion taking the Lagrange multipliers (i.e. contact forces) into account and the contact constraints have to be solved simultaneously. To this end, a new contact matrix $C_{NT}$ containing the normal and tangential constraints is built by considering that, because of high relative velocities between the two structures, only sliding occurs:

$$\left\{ \begin{array}{l} \boldsymbol{\lambda} = \left( C_N^t \left[ \frac{\mathbf{M}}{h^2} + \frac{\mathbf{D}}{2h} \right]^{-1} C_{NT} \right)^{-1} \mathbf{d}^p \\ \mathbf{u}^{n+1} = \mathbf{u}^{n+1,p} + \left[ \frac{\mathbf{M}}{h^2} + \frac{\mathbf{D}}{2h} \right]^{-1} C_{NT} \boldsymbol{\lambda} \end{array} \right. \qquad (11)$$

Finally, time is incremented before going back to the first step until final time is reached.

Concerning the contact matrix $C_{NT}$, it is assumed that the direction of sliding between the two structures can be approximated by the predicted one where the axial displacements of the casing are neglected.

# 5 VALIDATION

## 5.1 Case studies

The validation of the contact algorithm is carried out by comparing the results obtained to the proposed algorithm and the ones provided by another code using the bi-potential method [6] combined with an implicit time integration procedure. As mentioned in the introduction, the proposed contact case aims at validating the contact algorithm used for the simulation of blade tip/casing contacts. Because of the very large size of the 3D industrial models – up to several million dofs – used in these simulations, these models are not compatible with the implicit time integration procedure proposed in [6]. Consequently, a simplified contact case is presented in Fig. 1, consistent with a typical blade tip/casing contact case in terms of friction and relative speed.
<figure><img src="image_2.png" /><figcaption>**Figure 1: Validation contact case**</figcaption></figure>

A cube is launched between two partial rings clamped on their extremities (---). The initial speed of the cube $\vec{v}_c$ is along $O\vec{y}$ axis. We define $C_l$ and $C_r$ the respective centers of the left partial ring ($x < 0$) and the right partial ring ($x > 0$). Both the cube and the rings are meshed with 8-

<figure><img src="image_3.png" /><figcaption>**Figure 2: Contact faces**</figcaption></figure>

node linear brick elements. The cube is discretized with 125 identical elements (5 elements on each length). Regarding to the rings, they are meshed regularly with 5 elements along their thickness, 8 elements along their width and 40 elements along their length, so that each ring is discretized with 1600 elements. For the sake of brevity, the convergence study – in terms of eigenfrequencies of each structure and displacements of each node for a contact simulation – leading to these mesh parameters is not presented in this paper. As pictured in Fig. 2 two set of nodes are considered for managing contacts:

1. the *left face* ($x < 0$) with 36 nodes (●)

2. the *right face* ($x > 0$) with 36 nodes (●)
In the sequel, subscripts r and c respectively refer to the rings and the cube. It must be noticed that for all the contact cases considered, the problem is symmetric with respect to the planes $(O, \vec{y}, \vec{z})$ and $(O, \vec{y}, \vec{x})$. However, the simulations are carried out without any symmetry assumption to show that the contact treatment is precise enough to obtain symmetric results.

The simulation parameters for each contact case are given in Tab. 1. The initial speed of the

<table><caption>Table 1: Simulation parameters</caption><thead><tr><th>E<sub>r</sub><br>(N·m<sup>-2</sup>)</th><th>E<sub>c</sub><br>(N·m<sup>-2</sup>)</th><th>ρ<sub>c</sub>=ρ<sub>r</sub><br>kg·m<sup>-3</sup></th><th>ν<sub>c</sub>=ν<sub>r</sub></th><th>μ</th><th>||<br>ν̇<sub>c</sub>||<br>(m·s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>210·10<sup>9</sup></td><td>210·10<sup>7</sup></td><td>7800</td><td>0.3</td><td><b>0.15</b></td><td><b>-100</b></td></tr></tbody></table>

cube is set to -100 m·s⁻¹ which is high enough to neglect the effect of gravity in our study.

## 5.2 Symmetry of the quantities of interest

It is of great importance to ensure the symmetry of both displacements and contact forces since the contact case is symmetric. This is achieved by comparing displacements $u_x$ and contact efforts $\lambda$ for nodes 1, 6, 85 and 90 on the left face and their respective symmetric counterparts on the right face: nodes 31, 36, 91 and 96 as pictured in Fig. 2. Results presented in Figs. 3(a) and 3(b) show the excellent symmetry of both displacements and contact force over time. The symmetry of the results is observed for all the contact cases proposed in this section.

## 5.3 Comparison

Results given by the proposed method and the bi-potential method match very well while structures are in contact. Displacements comparison for node 1 may be seen in Figs. 4(a), 4(b) and 4(c) and contact force comparison appears in Fig. 4(d). One may observe that the slopes of the curbs obtained with each method for the y displacement in Fig. 4(b) are identical. This means that the speed of the cube after going across the rings is equal for both methods.

Also, high frequency vibrations of significant amplitude appear for $u_{1x}$ and $u_{1z}$ in Figs. 4(a) and 4(c) due to the absence of damping and making it difficult to compare the results. However, it is noticeable that the amplitudes are very similar for both methods. A phase difference may be observed in Fig. 4(a) between both methods for $t > 0.75$ ms, when structures are not in contact anymore. The fact that implicit and explicit time integration methods are responsible for different periodicity errors [14] is well known. Also, in our study, there is a difference in the discretization of the rings between the explicit time integration method, for which bicubic B-splines are considered, and the implicit time integration method for which the rings are
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Displacements and contact force for nodes 1, 6, 85 and 90 (—) and their respective symmetric: nodes 31, 36, 91 and 96 (······)</figcaption></figure>

discretized following the finite element mesh. This directly impacts phase difference between results by affecting the time for which contact is lost between structures.

Comparison with Ansys<fnref n="4" /> results is given in Fig. 5 for the time structures are in contact. This figures clearly emphasizes that results match. In terms of computation times, the proposed algorithm and Ansys are very close (±5%) while the bi-potential method leads to computation times about 10 times longer due to the fact that this code takes into account large deformations.

Because of the large size of the finite element model used for 3D blade-tip/casing contact simulations, computation times may be drastically reduced with the proposed algorithm if combined with reduced order models [15]. That is the reason why the proposed algorithm is preferred to Ansys in our case.

## 5.4 Energy conservation

The notion of stability for a nonlinear algorithm is associated with the respect of conservative laws, i.e. momentum and energy conservation [11]. The design of nonlinear stable algorithms

<footnote n="4">⁴contact simulation with Ansys is carried out using augmented Lagrangian method and implicit time integration ($α = 0.25, β = 0.5$)</footnote>
<figure><img src="image_5.png" /><figcaption>(a) $u_{1x}$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $u_{1y}$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) $u_{1z}$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) $\lambda_1$</figcaption></figure>

**Figure 4:** Comparisons between the proposed algorithm (······) and the bi-potential method (—) for case 3
<figure><img src="image_9.png" /><figcaption>**Figure 5:** Superimposition of the results obtained with the proposed algorithm between the proposed algorithm (······), the bi-potential method (—) and Ansys (—)</figcaption></figure>

satisfying conservative laws has been widely studied [4] and debated [16] over the previous years. In the particular case of contact dynamics between flexible or rigid bodies, energy conserving algorithms have been developed for different applications using the penalty method [17] or the bi-potential method [18].

Due to the use of the explicit central difference method, total angular momentum conservation is ensured [10, 9] in our study. For this reason, the focus is made on energy conservation. Since no damping is considered in our study, the total energy of the system $E_T$ – the sum of kinetic and strain energies of each structure – must remain constant over time when the friction coefficient is $\mu = 0$. When $\mu > 0$, friction between the cube and the rings induces dissipation and the total energy of the system decreases with time. The energy balances depicted in Fig. 6(a) and 6(b) are obtained with a lower norm of the initial velocity of the cube which is now –50 m·s⁻¹. This allows us to observe significant variations of both kinetic and strain energies which do not occur for higher velocities. The total energy of the system is computed as the sum of the kinetic energy and the strain energy of the cube and the rings.

In the case of Fig. 6(a), no friction is considered between structures so that the cube goes across the rings, contrary to the case of Fig. 6(b) where friction between structures leads to a configuration in which the cube is blocked between the rings. This explains why the kinetic energy tends to zero in Fig. 6(b) while the strain energy equals the total energy.

Over time, the total energy of the system only decreases by 0.16% from 1,218.75 J to 1,216.83 J when no friction is considered. These results confirm that the contact algorithm presented respects energy conservation over time and show that only friction introduces energy decay.
<figure><img src="image_10.png" /><figcaption>Figure 6: Total (—), kinetic (······) and strain (—) energies evolution</figcaption></figure>

# 6 APPLICATION: BLADE-TIP/CASING INTERACTION

## 6.1 Structural models and strategy

The considered compressor stage is partially depicted in figure 7. The bladed disk is perfectly balanced and tuned, gyroscopic effects as well as centrifugal stiffening are neglected. It is constituted of 56 sectors, each of them discretized in 3,764 prismatic finite elements connected together through 6,734 nodes, reaching a total of 1,131,321 dof. The casing is constituted of 44,800 tetrahedric and hexahedric finite elements connected through 68,768 nodes, or equivalently 206,304 dof.

This finite element model is quite large and computationally expensive. Consequently, a component mode synthesis approach able to handle contact problems directly in the new reduced space is required. Due to this central motivation, we adopt the two following ideas:

1. for each blade, contact is anticipated on three nodes located in the leading edge, middle of the chord, and trailing edge. The slave contact surface on the casing will be built with the help of a spline connected to the corresponding facing nodes.

2. each structure is reduced using the well-known Craig-Bampton (CB) procedure. Accordingly, the three spatial displacements of the nodes where contact is anticipated are defined
as interface dof, all the remaining ones being internal dof<fnref n="5" />.

Structural damping is introduced in the model through a modal damping ratios, $\xi = 0.03$ for the casing and $\xi = 0.005$ for the bladed disk. The subsequent modal damping matrix is then projected back into the physical space. The contact surface of the casing is a closed surface to

<figure><img src="image_11.png" /><figcaption>**Figure 7: Representation of the spline used for contact management**</figcaption></figure>

the contrary of the rings considered for the validation case. Consequently, specific conditions – namely end-conditions – must be imposed on the B-splines in order to avoid numerical problems and to ensure the continuity of the normal to the surface as depicted in Fig. 8.

<figure><img src="image_12.png" /><figcaption>**Figure 8: Surface spline versus end-conditions along the circumferential direction. Free-edge end-conditions in the axial direction. Black dots • are the data points of the deformed casing. The control points and vertices are not shown.**</figcaption></figure>

<footnote n="5">⁵the computation of the reduced order models is not described in this paper for the sake of brevity. Details can be found in [19] regarding the theoretical background of the CB method and in [20] for its application to the same 3D industrial models.</footnote>
## 6.2 Convergence criterion

Once reduced order models are computed, their convergence toward the finite element model must be assessed before carrying out any simulation. The Residual Convergence Criterion (RCC) has been recently introduced [21] for quantitatively comparing eigenvectors obtained by the reduced models and the finite element model. From a mathematical point of view, given one real valued vectors $Z_i$ of frequency $f_i = \frac{\omega_i}{2\pi}$ from the reduced model, the RCC is defined as the ratio

$$ \epsilon_i = \frac{\mathbf{r}_{s,i}^T \mathbf{K} \mathbf{r}_{s,i}}{\mathbf{Z}_i^T \mathbf{\Phi}^T \mathbf{K} \mathbf{\Phi} \mathbf{Z}_i} \quad \text{with} \quad \mathbf{r}_{s,i} = \mathbf{K}^{-1} [\mathbf{K} - \omega_i^2 \mathbf{M}] \mathbf{\Phi} \mathbf{Z}_i \quad (12) $$

$\mathbf{r}_{s,i}$ is a non-zero vector due to the fact that $Z_i$ is computed from the reduced model. Along with this criterion, we consider the frequency deviation factor

$$ \Delta f = \frac{f_i - f_j}{f_i} \times 100(\%) \quad (13) $$

which is considered satisfactory when its $n$ first frequency deviation factors are smaller than 0.3%.

## 6.3 Bladed disk

For the sake of brevity – although both structures are reduced – convergence results are only presented for the bladed disk. Similar results are observed for the casing. Both the deviation

<figure><img src="image_13.png" /><figcaption>Figure 9: Frequency deviation factors (in %) of the bladed disk with respect to the number of component modes in the reduced model</figcaption></figure>

factors and the RCC appear by pair because of the cyclic-symmetry of the bladed disk. The RCC values depicted in Fig. 10 highlight the convergence of the reduced model since for any of the 56 first eigenvectors, RCC<1% when $\eta \ge 56$.
<figure><img src="image_14.png" /><figcaption>**Figure 10:** Residual convergence criterion (in %) of the bladed disk with respect to the number of component modes in the reduced model</figcaption></figure>

**6.4 Results**

It has to be noted that the investigated system is highly nonlinear because of contact, and as a result the criterion given in previous section only provides an approximative assessment of the reduced model's accuracy. The results presented in this section implicitly require that time marching simulations were conducted by increasing the number of component modes in order to examine the convergence of the reduced-order model involving nonlinear terms. All distances and time have been normalized in this section.

Contact simulations are carried out while a two nodal diameter effort is applied on the casing for several angular velocities. Physically, a two nodal diameter effort leads to an ovalization of the casing such as the one that may result – for instance – from thermal loads. Blade-tip/casing distances for the three nodes used to manage contact on the first blade – leading edge (1), middle of the chord (2) and trailing edge (3) – are depicted in Fig. 11 for one angular velocity. When

<figure><img src="image_15.png" /><figcaption>**Figure 11:** Blade-tip/casing distances</figcaption></figure>

the blade-tip/casing distance equals zero, then the blade is in contact with the casing. Thanks to the use of forward increment Lagrange multiplier method, residual penetrations are insignificant.

Cyclic symmetric structures usually present close natural frequencies [22] – particularly the
lowest ones – and a typical Fourier transformation of time results may be inefficient in revealing dominant modes. Consequently, it may be of great interest to analyze time results in terms of nodal diameter influence. As an example, for all the angular velocities considered, the influence

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

<figure><img src="image_17.png" /><figcaption>**Figure 12:** Influence of nodal diameters ($n_d$) 2 and 4 for different normalized angular velocities while a 2-nodal diameter effort is applied on the casing.</figcaption></figure>

of nodal diameters 2 and 4 is pictured in Figs. 12(a) and 12(b). It appears clearly that the participation of the fourth nodal diameter becomes negligible compared to the second nodal diameter when the angular velocity increases. Such result may be used for optimizing the computation of the reduced order models.

# 7 CONCLUSION

In this paper, a 3D contact algorithm based on the forward increment Lagrange multiplier method with explicit time integration has been proposed and validated with respect to a contact algorithm using a different contact management procedure and the commercial software Ansys. In
particular, the contact algorithm presented allows for energy conservation and is compatible
with the use of reduction techniques that are unavoidable when dealing with large industrial
finite element models. An application of this algorithm for a 3D blade-tip/casing contact sim-
ulation is given at the end of the paper. The results are pictured in terms of blade-tip/casing
distances and nodal diameter influence. They and allow to foresee a wide range of applications
for the study of blade-tip/casing contacts. Work is in progress to simulate abradable coating
wear on the casing.

8 ACKNOWLEDGEMENT

Thanks go to Snecma for its technical and financial support. This work takes place in the frame-
work of the MAIA mechanical research and technology program sponsored by CNRS, ONERA
and SAFRAN Group.

REFERENCES

[1] P. Schmiechen, *Travelling Wave Speed Coincidence*, Ph.D. thesis, Imperial College of Science, Technology and Medicine - University of London (1997).

[2] N. Kikuchi, *Penalty/finite element approximations of a class of unilateral contact problems*, Penalty Method and Finite Element MethodASME: New York.

[3] T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. Kanoknukulchai, *A finite element method for a class of contact-impact problems*, Computer Methods in Applied Mechanics and Engineering 8 (1976) 149-276.

[4] J. C. Simo, T. A. Laursen, *An augmented Lagrangian treatment of contact problems involving friction*, Computers & Structures 42 (1992) 97-116.

[5] P. Alart, A. Curnier, *A mixed formulation for frictional contact problems prone to Newton like solution methods*, Comp. Meth. Appl. Mech. Engng. 92 (1991) 353-375.

[6] G. De Saxcé, Z. Q. Feng, *The bi-potential method: a constructive approach to design the complete contact law with friction and improved numerical algorithms*, Mathematical and Computer Modeling 28 (1998) 225-245.

[7] C. Padova, J. Barton, M. Dunn, S. Manwaring, *Experimental Results From Controlled Blade Tip/Shroud Rubs at Engine Speed*, Journal of Turbomachinery 129 (2007) 713-723.
[8] A. Emery, J. Wolak, S. Etemad, S. Choi, *An Experimental Investigation of Temperatures due to Rubbing at the Blade-Seal Interface in an Aircraft Compressor*, Wear 91 (1983) 117-130.

[9] M. Legrand, B. Peseux, C. Pierre, *Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion*, in: Sixième Colloque National en Calcul des Structures, Giens, France, 2003.

[10] N. Carpenter, R. Taylor, M. Katona, *Lagrange constraints for transient finite element surface contact*, International Journal for Numerical Methods in Engineering 32 (1991) 103-128.

[11] T. Laursen, *Computational contact and impact mechanics*, Springer, 2002.

[12] E. Arnoult, *Modélisation numérique et approche expérimentale du contact en dynamique: application au contact aubes/carter de turboréacteur*, Ph.D. thesis, Université de Nantes (2000).

[13] D. Chamoret, P. Saillard, A. Rassineux, J.-M. Bergheau, *New smoothing procedures in contact mechanics*, Journal of Computational and Applied Mathematics 168 (1-2) (2004) 107-116.

[14] M. Géradin, D. Rixen, *Mechanical vibrations*, Wiley, 1997.

[15] A. Batailly, M. Legrand, P. Cartraud, C. Pierre, J. P. Lombard, *Evaluation of Component Mode Synthesis Methods for the Detection of Modal Interaction Through Rotor Stator Contacts*, Proceedings of the ASME 2009 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference: IDETC09, San Diego, California, 2009.

[16] M. Ortiz, *A note on energy conservation and stability of nonlinear time-stepping algorithms*, Computers & Structures 24 (1986) 167-168.

[17] M. L. Kaplan, J. H. Heegaard, *Energy-conserving impact algorithm for the heel-strike phase of gait*, Journal of Biomechanics 33 (2000) 771-775.

[18] B. Magnain, *Développement d'algorithmes et d'un code de calcul pour l'étude des problèmes de l'impact et du choc*, Thèse de doctorat (phd thesis), Université d'Évry-Val d'Essonne, Évry, France (2006).

[19] R. Craig, M. Bampton, *Coupling of substructures for dynamic analyses*, AIAA JOURNAL 6 (7).

[20] A. Batailly, *Simulation de l'interaction rotor/stator pour des turbomachines aéronautiques en configuration non-accidentelle*, Thèse de doctorat (phd thesis), École Centrale de Nantes, Nantes, France (2008).
[21] A. Bobillot, *Méthodes de réduction pour le recalage. Application au cas d'Ariane 5*, Thèse de doctorat (phd thesis), Ecole Centrale de Paris, France (2002).

[22] R. Bladh, M. P. Castanier, C. Pierre, *Component-mode-based reduced order modeling techniques for mistuned bladed diskd-part I: theoretical models*, Journal of Engineering for Gas Turbines and Power 123 (2003) 89-99.