# Full 3D strategies for rotor-stator contact interaction in turbomachinery

Mathias Legrand, Christophe Pierre, Patrice Cartraud

► To cite this version:

Mathias Legrand, Christophe Pierre, Patrice Cartraud. Full 3D strategies for rotor-stator contact interaction in turbomachinery. 12th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Feb 2008, Honolulu, United States. hal-01008239

HAL Id: hal-01008239

https://hal.science/hal-01008239

Submitted on 16 Aug 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
FULL 3D STRATEGIES FOR ROTOR-STATOR CONTACT
INTERACTION IN TURBOMACHINERY

Mathias Legrand<sup>a</sup> Christophe Pierre<sup>a</sup> Patrice Cartraud<sup>b</sup>

ABSTRACT

In turbomachinery, direct contacts between bladed disks and surrounding casings occur through a variety of mechanisms that may lead to severe damages. These mechanisms still remain unclear and understanding their origins is a new challenge. On that subject, the present work deals with a new full 3D strategy in order to provide meaningful insights to designers. It involves reduced computational costs and a robust contact methodology able to account for high relative displacements at the contact interface. The first aspect of the problem is addressed through the reduction of the structures using the Craig-Bampton approach. The second one is treated with the use of spline concepts in order to smooth the contact surface and thus, avoid numerical issues. Accordingly, one set of interface nodes belonging to the area where contact is anticipated is defined for each structure. The contact constraints are enforced through a surface spline attached to the interface nodes of the casing. The respective contact forces are calculated using the Lagrange multiplier framework within an explicit time-marching procedure. Accuracy and computational costs are controlled by the size of the reduced-order models. The capabilities of this very versatile strategy encompass a wide range of finite element meshes and interaction problems that may occur in fan, turbine or compressor sections. First results show complex behaviors with coupling between the different modes of the structures and sensitivity to the friction coefficient.

INTRODUCTION

In rotordynamics, nonlinear coupling forces between
the rotating and surrounding stationary parts can give
rise to unexpected significant amplitudes of displace-

ment and subsequent high stresses. More specifically
in aircraft engines, several different mechanisms have
been identified and are usually classified in three main
categories:

* Interacting forces are due to variations of fluid pressure without structural contact.

* Interacting forces occur at a unique contacting point along the circumference of both structures.

* Interacting forces are induced by multiple simultaneous contacting points at different locations along the circumference.

From a structural point-of-view, the first two categories
are more or less well understood, and virtually all of the
prior literature focuses on these two specific problems
(see [Childs, 1993], [Muszynska and Goldman, 1995],
and [Lee, 1993] for instance). On the contrary,
the third category is an emerging field involv-
ing new research efforts (see [Schmiechen, 1997],
[Legrand et al., 2004] or [Lesaffre et al., 2007]). It
mainly comes from the need for better engine effi-
ciency and smaller aerodynamic losses, which are to
be achieved by reducing the tip-clearance between the
stationary and rotating components of the engine. One
can thus expect a significant increase of structural rub
with origins as different as: gyroscopic effects under
certain operating conditions; maneuvering loads dur-
ing take-off and landing of the aircraft; rotor imbal-
ance due to design uncertainties and manufacturing im-
perfections; bird strikes or blade failure; vibration due
to aerodynamic excitation; and outer casing distortion
caused by vibration or temperature gradient.

Depending on the nature of the induced contact,
these interactions between the bladed disks and sur-
rounding casing can give rise to long-lived phenomena

<footnote><sup>a</sup> Structural Dynamics & Vibration Laboratory, McGill University, Montreal QC H3A 2K6, Canada,
{mathias.legrand, christophe.pierre}@mcgill.ca</footnote>

<footnote><sup>b</sup> GeM, Équipe Calculs et Structures, École Centrale de Nantes, Nantes 44000, France,
patrice.cartraud@ec-nantes.fr</footnote>
characterized by initially intermittent soft contacts that
can lead to undesirably large amplitudes and high stress
levels, due to the excitation of the mode shapes of the
structures.

In the present work, a full 3D investigation of this
phenomenon is presented. Computational costs are
decreased by independently reducing each structural
model using the well-known Craig-Bampton procedure.
Accordingly, two sets of interface nodes belonging to
the respective surfaces where contact is anticipated are
defined. The precision of the numerical simulations is
governed by the number of internal degrees-of-freedom
(dof) retained in the reduced-order models. The con-
tact constraints are enforced by a bi-cubic B-spline at-
tached to the casing, in order to avoid the usual mathe-
matical inconsistencies in contact treatment due to the
finite element mesh approximation - i.e. exact distances
between the two structures and continuous definition
of the outward normal direction. Corresponding La-
grange multipliers are calculated using an explicit time-
marching procedure.

STRUCTURAL MODELS & STRATEGY

The considered compressor stage is depicted in Fig-
ure 1. The total diameter is 500 mm. 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,
thus reaching a total of 1,131,321 dof. The casing is
constituted of a total of 44,800 tetrahedric and hexa-
hedric finite elements connected through 68,768 nodes,
precisely 206,304 dof. The green circumferences in Fig-
ure 1 show zero-displacement boundary conditions.

This finite element model is quite large and compu-
tationally expensive. Consequently, a component mode
synthesis approach able to handle contact problems di-
rectly in the new reduced space is required. Due to this
central motivation, we adopt the two following ideas,
that may be visualized in Figure 2:

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 attached to the corresponding facing nodes.

2. each structure is reduced using the well-known Craig-Bampton 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.

<figure><img src="image_2.png" /><figcaption>Fig. 1: 3D view of the structures</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: 3D view of three sectors of the bladed disk and the corresponding part of the casing. Chosen interface nodes are denoted by ● on the bladed disk and • on the casing, to which a surface spline is attached.</figcaption></figure>

MODEL REDUCTION

Component mode synthesis

In most industrial applications, the finite element mod-
els contain a large of dof leading to cumbersome com-
putation times that can be drastically reduced by us-
ing component mode synthesis procedures. The latter
consist in subdividing a structure into non-overlapping
components called substructures for which the inter-
face dof<fnref n="1" /> is preserved. Each substructure is then represented by a set of specific shapes including vibration normal modes, rigid body modes, static modes, and interface modes used to defined a rectangular coordinate transformation in order to reduce the size of the original problem. The dynamics of the reduced-order full system is retrieved, up to a precision controlled by the size of the truncated sets, by imposing the compatibility conditions on the interface dof.

Only the reduction process is considered here: the bladed disk and the casing are already seen as separate substructures to be reduced and sharing a common nonlinear contact interface.

### Craig-Bampton procedure

The Craig-Bampton method [Craig and Bampton, 1968] is probably the most popular CMS method and as been [Hurty et al., 1971] recommend as the most trustworthy method for all vibration analyzes, well-suited to the determination of the vibration modes of large systems. Once the substructures are defined, the Craig-Bampton reduction basis of each substructure is composed of a truncated set of fixed interface modes and a full set constraint modes:

**Fixed interface normal modes:** they are normal modes of vibration of the component with all interface dof fixed<fnref n="2" />. They are stored in matrix $\Phi_L$ and called component modes in what follows.

**Constraint modes:** they are the sequential static responses of each component to a unit displacement of interface dof $i$ while all its interface counterparts are fixed. This set is denoted by $\Phi_R$.

For each substructure ($i$), the original equations of motion are projected onto the new reduced-order space such as the displacement vector $\mathbf{u}^{(i)}$ becomes:

$$ \begin{aligned} \mathbf{u}^{(i)} &= \Phi^{(i)} \mathbf{u}_{cb}^{(i)} \\ \left( \begin{matrix} \mathbf{u}_f \\ \mathbf{u}_i \end{matrix} \right)^{(i)} &= \left[ \begin{matrix} \mathbf{I} & \mathbf{0} \\ \mathbf{\Phi}_R & \mathbf{\Phi}_L \end{matrix} \right]^{(i)} \left( \begin{matrix} \mathbf{u}_f \\ \mathbf{q} \end{matrix} \right)^{(i)} \end{aligned} \quad (1) $$

The dynamics of both structures (i.e. casing and bladed disk) is now seen as a full set of constraint modes whose contributions $\mathbf{u}_f$ are the physical displacements of the interface dof, complemented with a truncated set of component modes of contribution $\mathbf{q}$ that controls the precision of the numerical simulations. The number of component modes (size of $\mathbf{q}$) is substantially smaller than the number of internal physical dof (size of $\mathbf{u}_i$). The initial dynamic equations of component ($i$):

$$ \mathbf{M}^{(i)} \ddot{\mathbf{u}}^{(i)} + \mathbf{D}^{(i)} \dot{\mathbf{u}}^{(i)} + \mathbf{K}^{(i)} \mathbf{u}^{(i)} = \mathbf{F}^{(i)} \quad (2) $$

is transformed to:

$$ \mathbf{M}_{cb}^{(i)} \ddot{\mathbf{u}}_{cb}^{(i)} + \mathbf{D}_{cb}^{(i)} \dot{\mathbf{u}}_{cb}^{(i)} + \mathbf{K}_{cb}^{(i)} \mathbf{u}_{cb}^{(i)} = \mathbf{F}_{cb}^{(i)} \quad (3) $$

where $\mathbf{X}_{cb}^{(i)} = \Phi^{(i)t} \mathbf{X} \Phi^{(i)}$ with $\mathbf{X} \equiv \mathbf{M}, \mathbf{D}, \mathbf{K}$ and $\mathbf{F}_{cb}^{(i)} = \Phi^{(i)t} \mathbf{F}^{(i)}$.

The contact constraints can directly be treated in the reduced space, thus avoiding permanent backward and forward mappings to the physical space.

### Convergence criterion

The Modal Assurance Criterion (MAC) is commonly used [Ting et al., 1993] for quantitatively comparing pairs of shapes obtained by different methods, numerical or experimental. In this study, this criterion is used knowing that the “difference” between the two models decreases with a higher number of component modes in transformation (1). From a mathematical point of view, given two real valued vectors $\mathbf{V}_i$ of frequency $f_i$ (from the original model) and $\mathbf{V}_j$ of frequency $f_j$ (from the reduced model), MAC$_{\mathbf{W}}(\mathbf{V}_i, \mathbf{V}_j)$ is defined by the following relation:

$$ \mathrm{MAC}_{\mathbf{W}}(\mathbf{V}_i, \mathbf{V}_j) = \frac{(\mathbf{V}_i^T \mathbf{W} \mathbf{V}_j)^2}{(\mathbf{V}_i^T \mathbf{W} \mathbf{V}_i)(\mathbf{V}_j^T \mathbf{W} \mathbf{V}_j)} \quad (4) $$

and often comes together with the frequency deviation factor:

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

where $\mathbf{W}$ is a weighting square matrix, generally the mass matrix. The MAC value, always varying between 0 and 1, is a quantitative correlation factor of the two vectors: if MAC is equal to unity, then the two vectors are parallel; if MAC is equal to zero, then the two vectors are orthogonal with respect to $\mathbf{W}$. Specializing the above mathematical definition to the case of two eigenvectors or mode shapes, such a number represents a meaningful quality factor of a mode shape with respect to the corresponding reference one.

For the sake of simplicity, the convergence study is here limited to the frequency deviation factors with respect to the size of the reduced order model. The latter is considered satisfactory when its $n$ first frequency deviation factors are smaller than a threshold value, chosen here as 0.3%.

### Casing

Due to the choice of the contact interface dof on the casing, its reduced model contains 1,200 constraint modes plus a set of component modes. The frequency deviation factors of the casing are depicted in Figure 3: they feature a pairwise pattern due to the axi-symmetry of the structure which has double frequency pairs. Concerning the first twenty eigenfrequencies of the casing, the reduced model including 200 component modes is

<footnote n="1">¹ boundary or junction dof common to adjoining substructures</footnote>

<footnote n="2">² fixed means zero-displacement in this approach</footnote>
satisfactory since the frequency error is less than 0.2%.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Frequency deviation factors (in %) of the casing with respect to the number of component modes in the reduced model</figcaption></figure>

## Bladed disk

The convergence study of the bladed disk must follow a specific approach because the structure features cyclic-symmetry. The convergence factors have to be calculated by modal family (or, equivalently harmonic family peculiar to system with cyclic-symmetry, see [Bladh, 2001] for details) so that no harmonic is omitted in the reduced order model. It is then necessary to consider the diagram by pairs of 28 spatial harmonics knowing that the structure is constituted of 56 blades. Respective diagrams are shown in Figure 4.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Bladed disk frequency deviation factors (in %) of the first modal family (first bending) with respect to the number of component modes in the reduced model</figcaption></figure>

It is proved that the Craig-Bampton procedure is converging with an increasing number of component modes. At least, three families of component modes are required to accurately approximate the first three families (just family 1 of first bending is shown here) of free normal modes of the bladed disk.

As for the casing, the deviation factors appear by pair because of the cyclic-symmetry of the bladed disk.

## Nonlinear terms

It has to be noted that the investigated system is highly nonlinear because of contact, and as a result the used criterion only provides an approximative assessment of the reduced model's accuracy. In the following phases of the study, numerical simulations will be conducted by increasing the number of component modes in order to verify this convergence still holds with the nonlinear terms.

## GENERAL BACKGROUND ON B-SPLINES

In classical contact algorithms, a master-slave approach is used such that contacting master nodes of the bladed disk are forced not to penetrate 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 [Chamoret et al., 2004], [Padmanabhan and Laursen, 2001]. They increase the convergence rate of the contact algorithms by removing the facetization issues while preserving the original meshes. In the present study, we focus on bicubic B-spline surfaces.

## B-spline curves

### Uniform B-spline curve

B-spline bases [de Boor, 1993] allow for the construction of complex parametric spatial curves $c(t)^3$ by multiplying the spline basis functions $B_{ni}$ by a set of control points $Q_i^4$:

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

They result by mapping a parametric space $t$ defined along a knot sequence $t_{i=0,N-1}$ into Cartesian space through control points and are completely specified by the curve's control points, the curve's order and the B-spline basis functions as seen in the Equation (6).

Two-dimensional splines will be later used to enrich the contact surface functional properties, thus avoiding numerical issues. This surface will be characterized by its degree $n$ and the knot sequence defined here by the following considerations:

* Choosing $n = 3$ provides a $C^1$ continuity to the normal of the contact surface, which is of primary importance to correctly calculate the distance between the structures and to avoid numerical jumps

<footnote n="3">³ $c(t)$ is a vector-valued function of a parametric value $t$</footnote>

<footnote n="4">⁴ in a three-dimensional cartesian space, $Q_i = [Q_{xi}, Q_{yi}, Q_{zi}]$</footnote>
in the estimation of the contact forces. Lower degree polynomials do not provide sufficient control of the surface's shape and higher degree polynomials are computationally more cumbersome and prone to numerical error.

* The knot sequence can be either uniform or nonuniform. A curve is uniform if the spacing between all the knots is the same, in other words $t_{i+1} = t_i + 1$. For many cases, this method is too simplistic and chord length parametrization may be preferred because it more precisely reflects the geometry of the data points by proportionally spacing the knot sequence to the distance between the data points. However, the external dof of the casing and supports of the surface spline are equally spaced along the circumferential direction and almost equally space along the axial direction (see Figure 1). Moreover, a uniform parametrization greatly simplifies the implementation.

Consequently, cubic uniform B-spline curves are used in this research. Each segment $p$ of the curve<fnref n="5" /> can be then written in a matrix form:

$$c_p(t) = \mathbf{T}\mathbf{M}\mathbf{Q}_p \quad (7)$$

where $\mathbf{T} = (t^3, t^2, t, 1)$, $\mathbf{Q}_p^t = (Q_{p-1}, Q_p, Q_{p+1}, Q_{p+2})$ and:

$$\mathbf{M} = \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)$$

By construction, B-splines only approximate and do not pass through, or interpolate, the control points since it lies within the convex hull of its control vertices. Nevertheless, it seems legitimate, in our study, to better control the position of our spline surface knowing that it has to be as close as possible to the original geometric surface of the casing.

## Interpolation

It is often desired that B-splines interpolate a set of data points. The inversion method addresses this issue by finding the $N+2$ control points $Q_i$, given a set of $N$ data points $P_i$ to be interpolated. $N$ linear equations are generated based on the fact that $c(t)$ has to pass through the data points.

In the case of a uniform parametrization, segment $p$ of the spline curve reaches its extremal points $P_p$ and $P_{p+1}$ for extremal values of parameter $t$, respectively $t=0$ and $t=1$ for $p=0, \dots, N-2$:

$$\begin{align} P_p &= c_p(0) = \frac{1}{6}(Q_{p-1} + 4Q_p + Q_{p+1}) \tag{9} \\ P_{p+1} &= c_p(1) = \frac{1}{6}(Q_p + 4Q_{p+1} + Q_{p+2}) \nonumber \end{align}$$

where $c_p(1) = c_{p+1}(0)$. The inversion problem (9) can be recast in a matrix form such as:

$$\begin{pmatrix} \text{end } 1 \\ P_0 \\ P_1 \\ \vdots \\ P_{N-2} \\ P_{N-1} \\ \text{end } 2 \end{pmatrix} = \frac{1}{6} \begin{bmatrix} \text{end conditions } 1 \\ 1 & 4 & 1 & 0 & \cdots & 0 \\ 0 & 1 & 4 & 1 & 0 & \cdots & 0 \\ \vdots & & \ddots & & \vdots & & \vdots \\ 0 & \cdots & 0 & 1 & 4 & 1 & 0 \\ 0 & \cdots & 0 & 1 & 4 & 1 & 0 \\ \text{end conditions } 2 \end{bmatrix} \begin{pmatrix} Q_{-1} \\ Q_0 \\ Q_1 \\ \vdots \\ Q_{N-2} \\ Q_{N-1} \\ Q_N \end{pmatrix}$$

or, in a shorter way:

$$\mathbf{P} = \mathbf{A}\mathbf{Q} \quad (10)$$

where the two end-conditions required to uniquely solve the problem have to be specified.

## End Conditions

Various methods are available for defining the end-conditions of a spline curve. Coming back to the casing problem, two different types of end-condition have to be investigated: free edges along the axial direction of the structure and continuity and closure conditions along the circumferential direction. The first one is dealt with the double vertice approach:

$$Q_{-1} = Q_0 \quad \text{and} \quad Q_N = Q_{N-1} \quad (11)$$

which assumes that the curvature of the spline curve at each end is zero. The second one requires the position of the spline, its first $c(t)'$ and second $c(t)''$ spatial derivatives to be equal at the end tips, mathematically leading to:

$$\begin{align} c_0(0) &= c_{N-2}(1) \Leftrightarrow P_0 = P_{N-1} \nonumber \\ c'_0(0) &= c'_{N-2}(1) \Leftrightarrow -Q_{-1} + Q_1 = \nonumber \\ & \qquad -Q_{N-2} + Q_N \tag{12} \\ c''_0(0) &= c''_{N-2}(1) \Leftrightarrow Q_{-1} - 2Q_0 + Q_1 = \nonumber \\ & \qquad Q_{N-2} - 2Q_{N-1} + Q_N \nonumber \end{align}$$

<figure><img src="image_6.png" /><figcaption>Fig. 5: two types of end condition</figcaption></figure>

<footnote n="5">⁵ generally vector-valued and written here in only one direction of the Cartesian space in what follows for the sake of</footnote>
## Bicubic uniform B-spline surfaces

### Construction

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} B_{ni}(s) B_{nj}(t) \quad (13) $$

involving two parametric values s and t defined on two knot vectors ($s_i$) and ($t_j$), and a set of control points $\mathbf{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, $\mathbf{S}(s,t)$ is a vector-valued function and each patch<fnref n="6" /> $pq$ of a bicubic uniform B-spline surface can be described in a matrix form by:

$$ S_{pq}(s,t) = \mathbf{S} \mathbf{M} \mathbf{Q}_{pq} \mathbf{M}^{\mathrm{T}} \mathbf{T}^{\mathrm{T}} \quad (14) $$

with $\mathbf{S} = (1, s^2, s^3)$, $\mathbf{T} = (1, 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} $$

Matrix $\mathbf{M}$ is defined in Equation (8). As shown in Equation (15), each bicubic B-spline patch is affected locally by sixteen control points.

### Interpolation

For the sake of simplicity $\mathbf{Q}$ (respectively $\mathbf{P}$) now denotes a reorganization of matrix $\mathbf{Q}_{pq}$ (respectively $\mathbf{P}_{pq}$) for $p = -1, \dots, M$ and $q = -1, \dots, N$ in a vector form such as, by choice<fnref n="7" />:

$$ \begin{aligned} \mathbf{Q}^{\mathrm{T}} &= (Q_{-1-1}, Q_{0-1}, \dots, Q_{M-1}, Q_{-10}, Q_{00}, \dots, Q_{MN}) \\ \mathbf{P}^{\mathrm{T}} &= (P_{-1-1}, P_{0-1}, \dots, P_{M-1}, P_{-10}, Q_{00}, \dots, P_{MN}) \end{aligned} \quad (15) $$

where vector $\mathbf{P}$ in Equation (15) contains non-physical points $P_{-1-1}$, $P_{0-1}$, ..., $P_{M-1}$, ... as end-conditions. The two dimensional counterpart of Equation (10) together with notations (15) and the appropriate end-conditions into Equation (13) leads to the following relationships between data and control points:

$$ \mathbf{P} = \mathbf{A} \mathbf{Q} \quad (16) $$

where matrix $\mathbf{A}$, of size $(N+2)^2 \times (M+2)^2$, is defined as the Kronecker product of matrices $\mathbf{A}_s$ and $\mathbf{A}_t$, counterparts of matrix $\mathbf{A}$ of Equation (10) with the appropriate end-conditions, respectively in u and v directions:

$$ \mathbf{A} = \mathbf{A}_s \otimes \mathbf{A}_t \quad (17) $$

Inverting matrix $\mathbf{A}$ becomes:

$$ \mathbf{A}^{-1} = \mathbf{A}_{\mathbf{s}}^{-1} \otimes \mathbf{A}_{\mathbf{t}}^{-1} \quad (18) $$

Accordingly, finding the control points of the surface

<figure><img src="image_7.png" /><figcaption>(a) free-edge end-conditions (first and second spatial derivatives) and compatibility in displacement (closure)
(b) continuity conditions (first and second spatial derivatives) and compatibility in displacement (closure)</figcaption></figure>

<caption>Fig. 6: 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.</caption>

spline reduces to inverting two constant smaller matrices $\mathbf{A}_s$ and $\mathbf{A}_t$ and is not computationally expensive. It can be carried out once at the beginning of the general algorithm.

## CONTACT DYNAMICS

### 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 [Laursen, 2002].

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 $\mathbf{x}$ belonging to the master surface $\Gamma_c^{(m)}$, its closest counterpart $\bar{\mathbf{y}}$ on the slave surface $\Gamma_c^{(s)}$:

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

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)}(\bar{\mathbf{y}}(\mathbf{x}))) \cdot \mathbf{n} \quad (20) $$

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

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

<footnote n="7">⁷ recasting matrix $\mathbf{Q}_{pq}$ by increasing first index q is also possible</footnote>
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 (21)$$

where $t_N$ stands for the contact pressure, assumed positive, acting on $\Gamma_c^{(s)}$. 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 (22)$$

for which $\mu$ is the coefficient of friction, $\mathbf{v}_T$, the tangential slip and $\mathbf{t}_T$, the tangential stress vector. The weak form of the contact problem can be then written in the following manner:

find the displacement field $\mathbf{u}$ such as for all admissible virtual displacement $\delta\mathbf{u}$:

$$\int_{\Omega} \rho \ddot{\mathbf{u}} \cdot \delta \mathbf{u} \, dV + \int_{\Omega} \bar{\bar{\sigma}} : \delta \bar{\bar{\varepsilon}} \, dV = \int_{\Gamma_{\sigma}} \mathbf{t}_d \cdot \delta \mathbf{u} \, dS \\ + \int_{\Gamma_c^{(s)}} (t_N \delta g + \mathbf{t}_T \cdot \delta \mathbf{u}_T) \, dS + \int_{\Omega} \mathbf{f}_d \cdot \delta \mathbf{u} \, dV \qquad (23)$$

where $t_N$ and $\mathbf{t}_T$ are constrained by conditions (21) and (22). Geometry and displacement of this formulation are discretized using the shape functions of the original prismatic, hexahedric, and tetrahedric finite elements of both structures. It is then reduced using the Craig-Bampton approach explained above. 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. 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 with Equations (16) and (18).

## Solution algorithm for interaction

Many different time-marching procedures dealing with contact have been developed in the past few years depending on the type of correction used (displacement or velocity) [Vola et al., 1998]. Because of its simplicity, the Forward Increment Lagrange Method [Carpenter et al., 1991] is used in our study. The matrix form of (23) is discretized in time using the explicit central differences scheme. The algorithm is then divided into three steps:

1. **prediction** at time step $n+1$ of the displacements $\mathbf{u}^8$ of the structures without considering any con-

tact. 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) \qquad (24)$$

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.

2. **determination** of the gap function vector $\mathbf{g}^p$ between the two structures using Equation (20). Each gap function where a penetration has been detected is kept in $\mathbf{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:

$$\mathbf{g}^{n+1} = \mathbf{C}_N^t \mathbf{u}^{n+1,c} + \mathbf{g}^p = \mathbf{0} \qquad (25)$$

where the superscript $^c$ means that the correction of the displacements is being calculated. $\mathbf{C}_N$ is the contact constraint matrix in the normal direction. 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 $\mathbf{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{aligned} \lambda &= \left( \mathbf{C}_N^t \left[ \frac{\mathbf{M}}{h^2} + \frac{\mathbf{D}}{2h} \right]^{-1} \mathbf{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} \mathbf{C}_{NT} \lambda \end{aligned} \right. \qquad (26)$$

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

Concerning the contact matrix $\mathbf{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 (see Figure 1) are neglected.

# RESULTS

## Conditions of interaction

In order to smoothly initiate and maintain contact between the casing and the blade tips, it is understood that a temperature gradient quasi-statically distorts the casing on a 2-nodal diameter mode. The rotational

<footnote n="8">⁸ 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>
velocity of the bladed disk is $\Omega = 500 \text{ rad/s}$, thus leading to 125 m/s relative velocities between the blade tips and the casing. In table 1 is given a comparison of the different structural matrices. It is worth noting that the casing is very stiff and heavy compared to the bladed disk. The friction coefficient as well as the number of

<table><caption>Tab. 1: Matrix maximum absolute row sum norm in order to compare the mechanical characteristics of each component.</caption><thead><tr><th>matrix</th><th>casing</th><th>bladed disk</th></tr></thead><tbody><tr><td>stiffness</td><td>2.25 · 10<sup>8</sup></td><td>1.58 · 10<sup>4</sup></td></tr><tr><td>mass</td><td>3.94 · 10<sup>-2</sup></td><td>5.4 · 10<sup>-6</sup></td></tr><tr><td>damping</td><td>1.39 · 10<sup>1</sup></td><td>7.8 · 10<sup>-4</sup></td></tr></tbody></table>

component modes are left to investigation.

## Gap distances and time-step size

It is well-known [Wriggers, 2002] that nonlinearities strongly constrain the time-step size of explicit algorithms that are conditionally stable due to the linear approximation of the displacements from one time-step to the next one. For linear simulations, the critical time-step is given by the Courant criterion $h_c = 2/\omega_{max}$ above which arises a numerical divergence and where $\omega_{max}$ is the largest eigenvalue of the system. Obviously, the latter increases with a higher number of component normal modes. Keeping all them for both structures (500 for the casing and 1,700 for the bladed disk), the linear critical time step size are independently found to be $h \le 1.15 \cdot 10^{-6}$ s for the casing and $h \le 7 \cdot 10^{-7}$ s for the bladed disk. Anyway, Courant criterion only provides a quick estimate of the critical time-step size in case of nonlinear problems. It is then required to conduct several simulations while decreasing the time-step size until convergence is met. In the present study, the time-step is chosen to be $h = 10^{-7}$ s, about seven times small than the critical linear time-step size.

<figure><img src="image_8.png" /><figcaption>Fig. 7: Blade-tips to casing distances for $h = 10^{-7}$ s.</figcaption></figure>

The distances between the blade-tips and the casing are illustrated in Figure 7 where the three families of interface nodes are clearly distinguishable. They remain positive for all times, confirming effectiveness of the surface spline as well as of the explicit time marching algorithm for the chosen time-step

## Reduced-order model convergence

In what follows, the friction coefficient is $\mu = 0.2$. Due to its own mechanical characteristics (mass and stiffness in Table 1), the casing is almost insensitive to the interaction. A parameter study concerning the number of kept component modes shows that a few of them are necessary for a correct approximation of its dynamics. Consequently, in what follows, only the results of the bladed disk are analyzed.

As already mentioned, the MAC provides valuable information for linear systems only. For nonlinear systems, it is necessary to check the reduced-order model convergence by increasing the number of component modes and comparing the subsequent results, namely the displacements at different locations of the structures. These results are depicted in Figures 8 and 9 with the color code described in Table 2. First of all, it is worth noting that a too low number (case 1) or a too high number of component modes (case 5) leads to numerical difficulties as observed in Figure 8 for instance where nonphysical vibrations appear. In Figure 8

<table><caption>Tab. 2: Computational normalized times for an increasing number of component modes. The color code is used in Figures 8 and 9</caption><thead><tr><th></th><th>q<sub>c</sub></th><th>q<sub>bd</sub></th><th>time</th><th>color</th></tr></thead><tbody><tr><td>1</td><td>0</td><td>0</td><td>0.68</td><td></td></tr><tr><td>2</td><td>20</td><td>112</td><td>0.69</td><td></td></tr><tr><td>3</td><td>100</td><td>168</td><td>0.79</td><td></td></tr><tr><td>4</td><td>200</td><td>336</td><td>0.87</td><td></td></tr><tr><td>5</td><td>300</td><td>504</td><td>1</td><td></td></tr><tr><td colspan="4">μ = 0, q<sub>c</sub> = 20, q<sub>bd</sub> = 112</td><td></td></tr></tbody></table>

where the considered interface node does not contact the casing, it is shown that convergence with respect to the number of component modes is very slow not to say nonexistent, even though the displacement in average are well predicted. On the contrary, the convergence is really good when focusing on a contacting interface node. More specifically, if the bladed disk reduced model includes at least two families of component modes (first bending and first torsion) in agreement with the preferential torsional motion of the blade-tips due to friction-only leading edge nodes comes into contact with the casing-, the convergence is correct since adding extra modes does not significantly improve it.

Additionally, the normalized times versus the size of
<figure><img src="image_9.png" /><figcaption>Fig. 8: Displacements of the second interface node of blade 1. This node does not impact the casing.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Displacements of the third interface node of blade 1. This node does impact the casing.</figcaption></figure>

the reduced model establish that the main part of the computation is dedicated to the contact treatment and the construction of the spline. Larger models would not be computationally prohibitive.

Friction coefficient

It is also shown that the sensitivity of the system to friction is not negligible. In fact, a very important be-

<figure><img src="image_11.png" /><figcaption>Fig. 10: Tangential plan tip displacements of blade 43 during interaction</figcaption></figure>

havior of the blades in general is observed. Figure 10 shows a strong coupling between the radial and tangential displacements of the contact node, and moreover when friction is omitted, the fact that each blade tends to bend in the positive direction -i.e. same as the rotational direction- when excited through direct normal contact. Consequently, when friction is accounted for, the tangential displacement of the contact node of the blade tends to increase in the direction opposite to the relative displacement due to rotation between the two structures. More importantly, it shows that the blade design under investigation tends to close the clearance between the structures because of friction: actually, the radial displacements depicted in Figure 9(a) are such that blades undergo more compression with friction and are then subjected to higher strains and subsequent stresses. Unfortunately, this design highly favors a divergent behavior. A design opening the clearance under tangential load is much more desirable in terms of stability. Any other friction coefficients larger than $\mu = 0.2$ have led to divergence.

As expected, a torsional behavior of the blade is also observed (displacement of node 3 is larger in magnitude than displacements of nodes 1 and 2) when friction is considered in Figure 10.

CONCLUSION

The emphasis of the study has been placed on the development of a numerical tool able to handle different kinds of direct contact interaction situations that can occur between bladed disks and casings in turbomachines, and aircraft engines more specifically. In order to capture the phenomenon of interest, the designed
tool combines the use of:

* an original full 3D finite element model.

* the Craig-Bampton procedure in order to reduce the size of the original problem together with a surface spline able to smoothly handle the contact constraints.

* the Lagrange multiplier approach embedded in an explicit time-marching procedure for general dynamic predictions.

First numerical simulations show the qualitative good behavior of our approach that represents an interesting compromise between low computational times and accuracy, with contact constraints perfectly captured by the surface spline. Results show complex behaviors involving coupling between the different modes of the bladed disk and a strong dependency to the friction coefficient.

The presented model holds promise for carrying out future investigations of the detected instability, including the examination of damping, blade tip geometric design, mistuning (large and small) or the use of abrable coating to name a few, all in very general industrial applications (front fan, turbine, compressor).

## ACKNOWLEDGEMENT

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

## REFERENCES

[Bladh, 2001] Bladh, R. (2001). *Efficient Predictions of the Vibratory Response of Mistuned Bladed Disks by Reduced Order Modeling*. PhD thesis, The University of Michigan, Ann Arbor, USA.

[Carpenter et al., 1991] Carpenter, N., Taylor, R., and Katonka, M. (1991). Lagrange constraints for transient finite element surface contact. *International Journal for Numerical Methods in Engineering*, 32:103–128.

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

[Childs, 1993] Childs, D. (1993). *Turbomachinery Rotor-dynamics Phenomena, Modeling and Analysis*. Prentice Hall.

[Craig and Bampton, 1968] Craig, R. R. and Bampton, C. C. (1968). Coupling of substructures for dynamics analyses. *AIAA Journal*, 6(7):1313–1319.

[de Boor, 1993] de Boor, C. (1993). *B(asic)-Spline Basics in Fundamental Developments of Computer-Aided Geometric Modeling*. (Piegl L., ed.) Academic Press, Washington D.C.

[Hurty et al., 1971] Hurty, W. C., Collins, J. D., and Hart, G. C. (1971). Dynamic analysis of large structures by modal synthesis techniques. *Computers & Structures*, 1(4):535–563.

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

[Lee, 1993] Lee, C.-W. (1993). *Vibration Analysis of Rotor-tors*. Kluwer Academic Publishers.

[Legrand et al., 2004] Legrand, M., Idoux, L., Pierre, C., and Cartraud, P. (2004). Rotor-stator nodal diameter mode interaction by means of contact in an aircraft engine. In *The 11-th of Symposium on Transport Phenomena and Dynamics of Rotating Machinery*, Honolulu, Hawaii.

[Lesaffre et al., 2007] Lesaffre, N., Sinou, J. J., and Thouzerez, F. (2007). Stability analysis of rotating beams rubbing on an elastic circular structure. *Journal of Sound and Vibration*, 299(4-5):1005–1032.

[Muszynska and Goldman, 1995] Muszynska, A. and Goldman, P. (1995). Chaotic responses of unbalanced rotor/bearing/stator systems with looseness or rubs. *Chaos, Solitons & Fractals*, 5:1683–1704.

[Padmanabhan and Laursen, 2001] Padmanabhan, V. and Laursen, T. A. (2001). A framework for development of surface smoothing procedures in large deformation frictional contact analysis. *Finite Elements in Analysis and Design*, 37:173–198.

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

[Ting et al., 1993] Ting, T., Chen, T. L. C., and Twomey, W. (1993). Correlating mode shapes based on the modal assurance criterion. *Finite Elements in Analysis and Design*, 14(4):353–360.

[Vola et al., 1998] Vola, D., Pratt, E., Raous, M., and Jean, M. (1998). Consistent time discretization for a dynamical contact problem and complementarity techniques. *Revue Européenne des Éléments Finis*, 7:149–162.

[Wriggers, 2002] Wriggers, P. (2002). *Computational contact mechanics*. John Wiley & Sons.