Alain Batailly
email: alain.batailly@polymtl.ca

Mathias Legrand

Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for flexible casings with friction and abradable coating

Keywords: rotor/stator interaction, aeronautical impeller, rubbing, nonlinear dynamics, unilateral contact dynamics, friction, abradable coating interaction rotor/stator, rouet aéronautique, frottement aube/carter, dynamique non-linéaire, dynamique du contact, frottement, revêtement abradable

Unilateral contact induced blade/casing vibratory interactions in impellers: Analysis for flexible casings with friction and abradable coating

Nomenclature

Introduction

In a previous article [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF], a numerical strategy for the analysis of unilateral contact induced blade/casing vibratory interactions in industrial impellers was derived assuming that the surrounding casing is perfectly rigid. This assumption is frequently considered in the modeling of blade-tip/casing rubbing phenomena. Existing publications on experimental rigs involving axial compressors commonly underscore that the dynamics of the casing is often negligible [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. However, this is in contradiction with a recent study [START_REF] Almeida | Experimental analysis of dynamic interaction between a centrifugal compressor and its casing[END_REF] which suggests that significant displacements may be observed in such systems, possibly induced by a critical vibratory coupling known as modal coincidence [START_REF] Schmiechen | Travelling wave speed coincidence[END_REF]. In this context, modeling flexible casings in aeronautical engines is crucial and this article extends the numerical strategy presented in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF]. The structural flexibility of the casing significantly increases the numerical complexity of the contact treatment and a thorough analysis must be undertaken in order to validate the presented results. In that sense, this paper is complementary to [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF]: all the qualitative developments regarding the aliasing effect, the harmonic tracking procedure or the two-dimensional Fourier transform explained in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] are still relevant and this work extends their scope to the general framework of unilateral contact occurrences between two flexible structures. A fully nonlinear contact detection procedure is incorporated in the solution method in order to handle any contact interface, irrespective of its geometrical curvature and shape. In addition, this article extends the numerical developments of [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] to 3D friction forces using a Coulomb law assuming only sliding occur. The removal of the sacrificial abradable coating deposited on the casing inner surface in order to avoid contact-induced damages is also accounted for. Accordingly, the proposed developments shall be adapted with no supplementary effort to a variety of systems such as automotive turbochargers, micro-rotors, or gas turbines, and are suitable to configurations ranging from the rubbing of a single blade to modal interaction coupling the dynamics of a full bladed disk and the surrounding casing.

Recent publications dealing with abradable coatings show the need for an improved understanding of the structural behaviour of such materials during impact. Among the open questions is the mechanism of material transfers between the impacting structure and the abradable coating during rubbing events [START_REF] Strock | Aluminum based abradable material with reduced metal transfer to blades[END_REF][START_REF] Langenbrunner | Dynamic response of a metal and a cmc turbine blade during a controlled rub event using a segmented shroud[END_REF], which is one of the motivations for blade repair technological enhancements [START_REF] Nagy | Automated brazing of abrasive tips onto superalloy turbine blades[END_REF]. Experimenters mostly explore the static and dynamic characterizations of the abradable coating mechanical behaviour while numericists are investigating workable ways to model it. From an experimental standpoint, several studies may be found on the static characterization of abradable coatings [START_REF] Bounazef | The wear, deterioration and transformation phenomena of abradable coating BN-SiAlbounding organic element, caused by the friction between the blades and the turbine casing[END_REF]9,[START_REF] Johnston | The sensitivity of abradable coating residual stresses to varying material properties[END_REF]. Similarly, several abradability investigations carried out statically-or at very low speed-are available [START_REF] Yi | Friction and wear behaviour and abradability of abradable seal coating[END_REF][START_REF] Yi | Abradability evaluation and tribological behaviour of abradable seal coating[END_REF] in addition to a few scratch tests [START_REF] Ma | Investigation of abradable seal coating performance using scratch testing[END_REF][START_REF] Ma | Evaluation of abradable seal coating mechanical properties[END_REF].

In the context of turbomachinery, one should mention the pioneering work [START_REF] Borel | The wear mechanisms occurring in abradable seals of gas turbines[END_REF] on wear mechanisms in gas turbines. Recently, experimental investigations on the dynamic characterization of abradable coatings within a simulated turbomachine environment were reported with respect to: a high relative interaction speed [START_REF] Cuny | Experimental measurement of dynamic forces generated during short-duration contacts: application to blade-casing interactions in aircraft engines[END_REF][START_REF] Baiz | Experimental investigation of the blade/seal interaction[END_REF][START_REF] Sutter | Dynamic analysis of the interaction between an abradable material and a titanium alloy[END_REF], elevated temperatures [START_REF] Dadouche | Experimental evaluation of abradable seal performance at high temperature[END_REF] as well as a combination of high interaction speed and elevated temperatures [START_REF] Mandard | Interacting force estimation during blade/seal rubs[END_REF][START_REF] Padova | Development of an experimental capability to produce controlled blade tip/shroud rubs at engine speed[END_REF]. While the influence of thermal conditions has long been investigated [START_REF] Novinski | Modified zirconia abradable seal coating for high temperature gas turbine applications[END_REF], experimental work carried out by industrials is rarely published. The only study which suggests a numerical model of the abradable coating is [9]: it aims at reproducing HR15Y test results. Based on the finite element method, a micro model is also proposed in [START_REF] Faraoun | Modelling route for abradable coatings[END_REF] in order to correlate observations made in the micro-structure to the mechanical properties of the macro-structure. However, in the context of an expensive time marching strategy, such as the one employed in the present article, this model is not appropriate. To the authors' knowledge, only two models are available in the literature for the modeling of abradable coating in the context of full 3D simulations [START_REF] Legrand | Numerical investigation of abradable coating removal through plastic constitutive law in aircraft engine[END_REF][START_REF] Williams | Simulation of blade casing interaction phenomena in gas turbines resulting from heavy tip rubs using an implicit time marching method[END_REF]. The model proposed in [START_REF] Legrand | Numerical investigation of abradable coating removal through plastic constitutive law in aircraft engine[END_REF] is used in this paper. It relies on mono-dimensional plastic finite elements-called abradable elements in the following.

The second section of the article centers on the modeling of the casing, and attention is paid to the mathematical smoothing of its contact interface: the corresponding formulation and its integration within the time marching algorithm are extensively detailed. The third section describes the contact algorithm and its validation. The nonlinear procedure dedicated to the computation of blades/casing distances is outlined as the central ingredient of the proposed time-marching strategy and attention is also paid to the contact constraints matrices. The fourth section introduces the modeling and associated numerical implementation of the abradable coating removal process. The possibility of simulating sophisticated interactions featuring a localized yet full removal of the abradable coating is discussed. Finally, the last two sections introduce illustrative case studies of the problem of interest. The first case study explores the response sensitivity to the friction coefficient and the abradable coating mechanical properties and the second case study targets a structural vibration analysis in the vicinity of a critical angular speed.

Note Amplitudes are normalized with respect to reference data for confidentiality purposes.

Contact treatment between two flexible structures

The model of the impeller extensively described in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] is used in this work and thus not recalled. Instead the focus is on the casing, the mathematical smoothing of its contact interface and the nonlinear contact treatment procedure devoted to the accurate detection of penetrations between the structures.

Casing modeling

The casing, assumed perfectly tuned, exhibits axial symmetry and is meshed with usual families of 3D finite elements. In this context, the advantageous mathematical properties brought by the theory of cyclic symmetry are considered for the efficient computation of its reduced-order model: accordingly, it is first divided in N = 100 identical sectors as pictured in Fig. 1. The boundary conditions are prescribed in agreement with typical assembly conditions within an aeronautical engine: nodes belonging to the red hatched area in Fig. 1 Surrounding casing of the impeller used in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] sector is presented in Fig. 3. This mesh contains 4,685 nodes and 784 elements so that the full 360-degree structure has 374,800 nodes and 78,400 elements. As illustrated in Fig. 3, the M = 5 boundary nodes of the sector define the contact interface of this sector. In other words, N × M nodes are used for the definition of the full casing contact interface. This finite element model is reduced in the exact same fashion as the impeller in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF]: this combines cyclic symmetry to the Craig-Bampton component mode synthesis method [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF]. The only difference between the reduction procedure of the casing and the impeller lies in the fact that centrifugal stiffening is omitted for the casing. The reader may refer to [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] for more details regarding the notion of nodal diameters.

Contact surface functional smoothing

As evidenced in Fig. 3, all together the boundary nodes of each sector draw the supporting frame of the casing contact surface. It is explained in previous works [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF][START_REF] Batailly | A comparative study between two smoothing strategies for the simulation of contact with large sliding[END_REF] that there is a need for contact surfaces smoothing when dealing with nonlinear contact simulations with industrial finite element models. Such models are typically made of quadratic finite elements that feature a discontinuity of the normal to the outer surface when sliding from an element to another. This discontinuity may trigger numerical instabilities-commonly known as the chatter effect-with 3D node-to-surface contact algorithms [START_REF] Puso | A mortar segment-to-segment contact method for large deformation solid mechanics[END_REF][START_REF] Muñoz | Modelling unilateral frictionless contact using the null-space method and cubic B-Spline interpolation[END_REF]. Solutions exist to bypass such discontinuity through a modification of the finite element type, with Hermite-like elements [START_REF] Sauer | Enriched contact finite elements for stable peeling computations[END_REF][START_REF] Batailly | A comparative study between two smoothing strategies for the simulation of contact with large sliding[END_REF] for instance. They are promising yet not fully satisfying from an industrial standpoint as they usually lead to an increased computational cost. Competitive techniques involve contact surface smoothing through a mathematical artefact, be it a B-spline [START_REF] Muñoz | Modelling unilateral frictionless contact using the null-space method and cubic B-Spline interpolation[END_REF], a non-ubiform rational B-spline (NURBS) or Nagata patches [START_REF] Neto | Applying nagata patches to smooth discretized surfaces used in 3d frictional contact problems[END_REF], as an additional virtual layer on top of the finite element contact surface. This mathematical layer is solely dedicated to the computation of gap functions and the direction of contact forces. In the present work, a bicubic B-spline is used. The formulation of these splines can be efficiently combined with the considered contact algorithm. Additionally, they are particularly well suited to the reduced-order models and ensure a perfect symmetry of the casing inner contact surface at rest. The numerical developments detailed in the remainder form an extension of a previous work [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF]. The two-dimensional spline is here intended to reflect any contact interface geometry irrespective of its curvature in the axial direction for instance. Contrary to [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF], more elaborate impeller casings may be accommodated and a fully nonlinear contact detection procedure is implemented. The improvements with respect to the developments presented in [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF] are detailed in Tab. 1.

Bicubic B-spline and associated notations

A bicubic B-spline is typically defined by a set of control points, denoted Q = Q(i, j) with i = 1, . . . , N + 3 and j = 1, . . . , M + 2 in this paper. One may compute the location of the data points, denoted P = P(i, j) with i = 1, . . . , N + 1 and j = 1, . . . , M which actually define the surface, based on:

P ▾ = AQ ▾ (1)
P(1, 1) = P(N + 1, 1) P(5, 1) P(4, 1) P(3, 1) P(2, 1)

P(N, 1) P(N -1, 1)
P(N -2, 1) where A is a rectangular matrix referred to as the interpolation matrix. It should be noted that Q(i, j) ∈ R 3 and P(i, j) ∈ R 3 . The data points are the boundary nodes, pictured in Fig. 3, of the finite element model of the casing. The ▾ superscript indicates that points are arranged in a single column vector following the convention:

P(N -3, 1) v = 0 v = M u = 0 u = N S
P(i, j + 1) S i,j P(3, M ) P(2, M ) P(1, M ) = P(N + 1, M ) P(N, M ) P(N -1, M ) sector 1 v ր u ր boundary nodes
X ▾ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ X(1, 1) X(1, 2) X(1, 3) ⋮ X(i, j) ⋮ X(n, m -2) X(n, m -1) X(n, m) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , with: X = X(i, j), i = 1, . . . , n, j = 1, . . . , m (2)
The control points may be divided into two complementary sets: namely the interior points Q int = Q(i, j) for i = 2, . . . , N + 2 and j = 2, . . . , M + 1 and the boundary points

Q bo = Q(i, j) for i = 1, N + 3 and j = 1, M + 2,
and by construction:

Q bo ∩ Q int = {∅} Q bo ∪ Q int = Q (3)
From the definition of the boundary conditions, the following relation can be established:

Q ▾ = CQ ▾ int (4)
Note that the number of interior points Q int is equal to the number of data points P so that matrices A and C have the same size. Equations ( 1) and ( 4) shall thus be combined as follows:

Q ▾ = C [AC] -1 P ▾ = BP ▾ (5)
Equation ( 5) reverses the usual definition of a bicubic B-spline given in Eq. ( 1) as control points Q are obtained from a set of known data points P, that is the casing boundary nodes. In this article, the computation of the control points Q based on Eq. ( 5) is referred to as the computation of the spline surface. Accordingly, at each time-step of the solution algorithm, the spline is supported by the predicted coordinates of the casing boundary nodes.

Computation of the interpolation matrix A

Equation ( 6) defines a single patch S i,j (u, v) depicted in Fig. 3. Notations u and v refer to the intrinsic coordinates of any point belonging to the spline. For a bicubic B-spline, it may be shown that for u, v ∈ [0 ; 1]:

S i,j (u, v) = 1 36 v 3 v 2 v 1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -1 3 -3 1 3 -6 3 0 -3 0 3 0 1 4 1 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Q i,j ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -1 3 -3 1 3 -6 0 4 -3 3 3 1 1 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ u 3 u 2 u 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (6)
where:

Q i,j = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Q(i, j) Q(i, j + 1) Q(i, j + 2) Q(i, j + 3) Q(i + 1, j) Q(i + 1, j + 1) Q(i + 1, j + 2) Q(i + 1, j + 3) Q(i + 2, j) Q(i + 2, j + 1) Q(i + 2, j + 2) Q(i + 2, j + 3) Q(i + 3, j) Q(i + 3, j + 1) Q(i + 3, j + 2) Q(i + 3, j + 3) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (7)
with i = 1, . . . , N and j = 1, . . . , M -1. Data points are located at each corner of a patch for (u, v) = (0, 0), (u, v) = (0, 1), (u, v) = (1, 1), and (u, v) = (1, 0) which yields the coordinates of the data points:

P(i, j) =S i,j (0, 0) = 1 36 Q(i, j) + 4Q(i, j + 1) + Q(i, j + 2) + 4Q(i, j + 3) + 16Q(i + 1, j + 1) + 4Q(i + 1, j + 2) + Q(i + 2, j) + 4Q(i + 2, j + 1) + Q(i + 2, j + 2) P(i, j + 1) =S i,j (0, 1) = 1 36 Q(i + 1, j) + 4Q(i + 1, j + 1) + Q(i + 1, j + 2) + 4Q(i + 2, j) + 16Q(i + 2, j + 1)4Q(i + 2, j + 2) + Q(i + 3, j) + 4Q(i + 3, j + 1)
+ Q(i + 3, j + 2)

P(i + 1, j) =S i,j (1, 0) = 1 36 Q(i, j + 1) + 4Q(i, j + 2) + Q(i, j + 3) + 4Q(i + 1, j + 1) + 16Q(i + 1, j + 2) + 4Q(i + 1, j + 3) + Q(i + 2, j + 1) + 4Q(i + 2, j + 2) + Q(i + 2, j + 3)
P(i + 1, j + 1) =S i,j (1, 1) = 1 36 Q(i + 1, j + 1) + 4Q(i + 1, j + 2) + Q(i + 1, j + 3) + 4Q(i + 2, j + 1) + 16Q(i + 2, j + 2) + 4Q(i + 2, j + 3) + Q(i + 3, j + 1) + 4Q(i + 3, j + 2) + Q(i + 3, j + 3) (8)
From these four equations arises the pattern of the general linear relationship between the data and control points in Eq. ( 1) with:

A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (9)
The matrix A is independent of the coordinates of the data points P and does not have to be updated during computation.

Spline boundary conditions: computation of matrix C

As mentioned above, matrix C is related to the boundary conditions imposed on the spline surface. Based on the notations introduced in Fig. 3, four edges to the spline surface are considered:

• v = 0 and v = M : these two edges-which refer respectively to the control points Q(i, 1) with i = 2, . . . , N and Q(i, M + 2) with i = 2, . . . , N -are defined as free edges:

Q(i, 1) = 2Q(i, 2) -Q(i, 3), i = 2, . . . , N Q(i, M + 2) = 2Q(i, M + 1) -Q(i, M ), i = 2, . . . , N (10)
• u = 0 and u = N : this is where the spline surface is closed and refers to the control points Q(1, j) with j = 1, . . . , M + 2 and Q(N + 3, j) with j = 1, . . . , M + 2. Along the closing line, it is ensured that first and second derivatives of the position with respect to u are equal one to another:

∂S i=1,j ∂u (u = 0, v) = ∂S i=N,j ∂u (u = 1, v), j = 1, . . . , M -1 ∂ 2 S i=1,j ∂u 2 (u = 0, v) = ∂ 2 S i=N,j ∂u 2 (u = 1, v), j = 1, . . . , M -1 (11)
These conditions lead to the following 2(M + 2) equations:

Q(1, 1) = 2Q(2, 2) -Q(2, 3) + 2Q(N + 1, 2) -Q(N + 1, 3) . . . ⋅ ⋅ ⋅ -2Q(N + 2, 2) + Q(N + 2, 3) Q(1, 2) = Q(2, 2) + Q(N + 1, 2) -Q(N + 2, 2) Q(1, 3) = Q(2, 3) + Q(N + 1, 3) -Q(N + 2, 3) ⋮ Q(1, M ) = Q(2, M ) + Q(N + 1, M ) -Q(N + 1, M ) Q(1, M + 1) = Q(2, M + 1) + Q(N + 1, M + 1) -Q(N + 2, M + 1) Q(1, M + 2) = -Q(2, M ) + 2Q(2, M + 1) -Q(N + 1, M ) . . . ⋅ ⋅ ⋅ + 2Q(N + 1, M + 1) + Q(N + 1, M ) -2Q(N + 2, M + 2) Q(N + 3, 1) = -2Q(2, 2) + Q(2, 3) + 2Q(3, 3) -Q(3, 3) + 2Q(N + 2, 2) -Q(N + 2, 3) Q(N + 3, 2) = -Q(2, 2) + Q(3, 2) + Q(N + 2, 2) Q(N + 3, 3) = -Q(2, 3) + Q(3, 3) + Q(N + 2, 3) ⋮ Q(N + 3, M ) = -Q(2, M ) + Q(3, M ) + Q(N + 2, M ) Q(N + 3, M + 1) = -Q(2, M + 1) + Q(3, M + 1) + Q(N + 2, M + 1) Q(N + 3, M + 2) = Q(2, M ) -2Q(2, M + 1) -Q(3, M ) . . . ⋅ ⋅ ⋅ + 2Q(3, M + 1) -Q(N + 2, M ) + 2Q(N + 2, M + 1) (12)
All these equations sum up to Eq. ( 4) with:

C = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2 -1 0 . . . 0 2 -1 0 . . . 0 -2 1 0 . . . 0 1 0 0 . . . 0 1 0 0 . . . 0 -1 0 0 . . . 0 0 1 0 . . . 0 0 1 0 . . . 0 0 -1 0 . . . 0 ⋮ ⋱ ⋮ 0 . . . 0 1 0 0 . . . 0 1 0 0 . . . 0 -1 0 0 . . . 0 0 1 0 . . . 0 0 1 0 . . . 0 0 -1 0 . . . 0 -1 2 0 . . . 0 -1 2 0 . . . 0 1 -2 2 -1 0 0 0 0 0 . . . 0 0 0 0 0 0 0 1 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 1 0 0 0 0 0 . . . 0 0 0 0 0 0 0 ⋮ ⋱ ⋮ 0 0 0 0 0 0 0 . . . 0 0 0 -1 2 0 0 0 0 0 0 0 0 0 . . . 0 0 0 0 -1 2 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 -1 2 -2 1 0 . . . 0 2 -1 0 . . . 0 2 -1 0 . . . 0 -1 0 0 . . . 0 1 0 0 . . . 0 1 0 0 . . . 0 0 -1 0 . . . 0 0 1 0 . . . 0 0 1 0 . . . 0 ⋮ ⋱ ⋮ 0 . . . 0 -1 0 0 . . . 1 0 0 . . . 0 0 1 0 0 . . . 0 0 -1 0 . . . 0 1 0 . . . 0 0 0 1 0 . . . 0 1 -2 0 . . . -1 2 0 . . . 0 0 -1 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (13)

Contact algorithm

Similarly to the numerical strategy considered in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF], the present work relies on an explicit time integration scheme combined with a Lagrange multiplier-based contact treatment procedure [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF]. Because the casing is now flexible, the complexity of the algorithm-detailed in Fig. 4-is significantly increased. In the following, the considered time step n + 1 appears as a subscript of the quantities of interest. The superscript * refers to the predicted quantities for which possible blade/casing contacts are ignored. As pictured in Fig. 4, predicted displacements of the casing v * c,n+1 and the impeller v * i,n+1 are first computed and potential blade/casing contact forces are omitted. From these displacements, the absolute coordinates of the casing boundary nodes P * n+1 are retrieved to compute the spline surface Q * n+1 through Eq. ( 5). Then, a nonlinear iterative procedure, detailed in section 3.1, finds the impeller/casing distances d * n+1 for each boundary node of the impeller and detects potential impeller/casing penetrations. From this point, two scenarii are to be considered: (1) one or more penetrations are detected, min(d * n+1 ) < 0, and non-admissible displacements displacements prediction for the impeller and the casing

v * i,n+1 = M i h 2 + D i 2h -1 2M i h 2 -K i v i,n + D i 2h -M i h 2 v i,n-1 + F ext i,n+1 v * c,n+1 = M c h 2 + D c 2h -1 2M c h 2 -K c v c,n + D c 2h -M c h 2 v c,n-1 + F ext c,n+1
coordinates of the casing boundary nodes (data points)

P * n+1 = P + v * c,n+1
construction of the spline surface along the casing contact surface

Q * n+1 = f (P * n+1 )
impeller/casing distance computation for each boundary node -contact and contact constraints matrices of the impeller : MC i,n+1 , VC i,n+1 and VCT i,n+1

d * n+1 = f (Q * n+1 , v * i,n+1 ) case
-contact forces :

λ n+1 = MC c,n+1 + MC i,n+1 -1 d * n+1
-impeller displacements correction :

v i,n+1 = v * i,n+1 -M i h 2 + D i 2h -1 VCT i,n+1 λ n+1
-casing displacements correction :

v c,n+1 = v * c,n+1 -M c h 2 + D c 2h -1 VCT c,n+1 λ n+1 - case 2 : no penetration detected min(d * n+1 ) ≥ 0 -predicted displacements are admissible : v i,n+1 = v * i,n+1 v c,n+1 = v * c,n+1
time update : ) > 0, and the predicted displacements are admissible; time is then updated.

t = t + h n = n + 1

Impeller/casing distances d * n+1

The computation of the impeller/casing distances d * n+1 is the corner stone of the contact treatment procedure. At a given time step n + 1, the predicted position of the k th impeller boundary node Γ * ,k n+1 is pictured in Fig. 5. The formulation of a bicubic B-spline is

starting point : S(u k n , v k n ) final iteration : S * (u k n+1 , v k n+1 )
iteration 1 P(i, j) P(i + 1, j) such that the determination of the closest counterpart S * (u k n+1 , v k n+1 ) of Γ * ,k n+1 , see Fig. 5, which belongs to the spline contact surface, is a nonlinear algebraic problem which is addressed by the iterative Newton-Raphson procedure detailed in Fig. 6. Since the explicit time integration scheme features by nature small time steps, the result S(u k n , v k n ) of this iterative procedure at the previous time step n provides an ideal starting point and convergence is reached within a handful number of iterations. For the very first time step, the initial guess is defined by the center of the closest patch that is easily determined through the distance of the impeller boundary node to the casing boundary nodes. The stop criterion mentioned in Fig. 6 is ε NR = 10 -5 ; the sensitivity of the results to this parameter is assessed in section 3.3.

P(i + 1, j + 1) P(i, j + 1) S i,j Γ * ,k n+1 v ր u ր

Contact matrices MC c,n+1 and MC i,n+1

When penetrations are detected, that is min(d * n+1 ) < 0, the correction of the displacements mentioned in Fig. 4 relies on the computation of contact matrices MC i,n+1 and MC c,n+1 respectively related to the impeller and the casing. These matrices are obtained following a four-step procedure:

1. computation of the outer normal vector n k out,n+1 to the spline contact surface at S * (u k n+1 , v k n+1 ), see Fig. 7:

n k out,n+1 = - ∂ → S * (u k n+1 , v k n+1 ) ∂u ∧ ∂ → S * (u k n+1 , v k n+1 ) ∂v (14)
2. computation of the tangential vector t k n+1 , see Fig. 7, to the spline contact surface at S * (u k n+1 , v k n+1 ):

t k n+1 = n k out,n+1 ∧ n k out,n+1 ∧ → Γ k n Γ * ,k n+1 (15)
3. computation of the contact constraints matrices: impeller: the contact constraints matrices of the impeller are obtained exactly as in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF][START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF] and is thus not recalled for the sake of brevity. casing: the spline surface governs the computation of the normal and tangential contact constraint matrices VC c,n+1 and VCT c,n+1 :

• first of all, the patch S * ,i,j n+1 to which belongs the point S * (u k n+1 , v k n+1 ) is determined, • the restriction of matrix B in Eq. ( 5) to the 16 control points Q i,j involved in the definition of patch S * ,i,j n+1 is denoted B, the k th boundary node of the impeller is denoted :

Γ * ,k n+1 = (x * ,k n+1 , y * ,k n+1 , z * ,k n+1 ) finding S * (u k n+1 , v k n+1
) the closest point from Γ * ,k n+1 on the spline -Newton-Raphson iterative procedure -solution of f n+1 (u k n+1 , v k n+1 ) = 0 where :

f n+1 (u k n+1 , v k n+1 ) =     -------------→ Γ * ,k n+1 S * (u k n+1 , v k n+1 ) • ∂ ----------→ S * (u k n+1 , v k n+1 ) ∂ u -------------→ Γ * ,k n+1 S * (u k n+1 , v k n+1 ) • ∂ ----------→ S * (u k n+1 , v k n+1 ) ∂ v     -leads to (u k n+1 , v k n+1 ), stop criterion: f n+1 (u k n+1 , v k n+1 ) ≤ ε NR computation of the bladed disk/casing distances vector dav * n+1 dav * n+1 (k) = -------------→ Γ * ,k n+1 S * (u k n+1 , v k n+1 )
determination of the sign of dav * n+1 (k) : negative if penetration, positive otherwise P(i + 1, j) • from (u k n+1 , v k n+1 ) and Eq. ( 6), the coefficient vector c i,j to be applied on the control points Q i,j so that:

d * n+1 (k) = min(0, dav * n+1 (k))
P(i + 1, j + 1) P(i, j + 1) S i,j Γ k n Γ * ,k n+1 (penetrated) S * (u k n+1 , v k n+1 ) n k out,n+1 t k n+1 ----→ Γ k n Γ * ,k n+1 ∂ -----------→ S * (u k n+1 , v k n+1 ) ∂u ∂ -----------→ S * (u k n+1 , v k n+1 ) ∂v
→ Γ * ,k n+1 S * (u k n+1 , v k n+1 ) = 0 (16)
is computed,

• contact constraint matrices are then constructed:

VC c,n+1 = VC 1 c,n+1 . . . VC k c,n+1 . . . VC N M c,n+1 VCT c,n+1 = VCT 1 c,n+1 . . . VCT k c,n+1 . . . VCT N M c,n+1 (17)
with:

VC k c,n+1 = B ⋅ c i,j ⋅ n k ext,n+1 VCT k c,n+1 = B ⋅ c i,j ⋅ (n k ext,n+1 + µt k ext,n+1 ) (18)
4. Finally, the contact matrices for each structure are [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF]:

MC c,n+1 = VC T c,n+1 M c h 2 + D c 2h -1 VCT c,n+1 MC i,n+1 = VC T i,n+1 M i h 2 + D i 2h -1 VCT i,n+1 (19)
and the displacement fields are corrected and updated as explained in Fig. 4.

Numerical validation

The exposed numerical procedure is highly nonlinear and thus requires a thorough validation. The main parameters include: (1) the stop criterion ε NR of the Newton-Raphson procedure given in section 3.1, (2) the time step h of the time integration scheme, (3) the friction coefficient µ in Eq. ( 18), ( 4) the basis used in the computation of the reduced-order models of the casing and the impeller and ( 5) the M boundary nodes retained on each sector of the casing describing the contact interface.

Stop criterion ε NR

Two values of the stopping criterion ε NR are considered in this section: default ε NR = 10 -5 and ε NR = 10 -10 to accurately check for convergence. Displacements of the casing and the impeller are drawn in Figs. 8(a) and 8(b), respectively. The results for both stop criteria are perfectly superimposed. In the context of 3D finite element models with linear elasticity considered in this study, the expected levels of vibration remain low. Typically, the radial displacement of a casing boundary node is two to three orders of magnitude smaller than the casing diameter which implies that the curvature of one spline patch has a negligible effect. Were larger radial displacements to be considered, the proposed strategy would still be efficient but a much lower stopping criterion would be required for converged results.

Friction forces and numerical sensitivity

Friction is usually sensitive from a numerical standpoint. Two numerical simulations are carried out, one with no friction µ = 0 and one with a very low friction coefficient µ = 0.01. Associated results are depicted in Figs. 9. The point is to ensure that the introduction of friction does not generate numerical instabilities. Due to friction, the results for each simulation are obviously slightly different yet almost identical which emphasizes that there is no numerical instability related to µ. As the friction coefficient increases, significant discrepancies shall arise.

Convergence analysis

The convergence of the results with respect to the other numerical parameters such as the time step h and the dimension of modal reduction bases were assessed but are not detailed. Convergence is observed with respect to each of these parameters and the numerical simulations carried out in the following consider h = 10 -7 s, and ten modes are retained for each nodal diameter of the casing and the impeller in their respective reduction basis, see [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] for more details. The number M of boundary nodes for each sector of the casing significantly influences the results but convergence is found for M ≥ 5. It should be mentioned here that increasing M has a fairly low impact on the computation times (multiplying M by 2 leads to an increase of computation times of about 20 %).

Other tests carried out in order to validate the numerical strategy feature the respect of symmetry in the displacements when a symmetric contact scenario is considered. Finally, the residual penetrations within the casing after correction of the displacements are systematically negligible (about 8 orders of magnitude lower than the level of vibration of the associated boundary node).

Abradable coating modeling

In addition to blade/casing unilateral contact and friction, the proposed developments also include the modeling of abradable coating, a sacrificial material deposited on the casing contact surface in order to mitigate blade/casing contacts. The removal of this material is accounted for through a residual plastic compression ∆ε p of the abradable elements as presented in [START_REF] Legrand | Numerical investigation of abradable coating removal through plastic constitutive law in aircraft engine[END_REF]. By choice, one element is used over the thickness of the abradable layer which is uniform over the casing contact surface. With this relatively simple model a real-time update of the contact interface is realized and recent numerical investigations have highlighted its relevance for the modeling of the interaction forces, be it abradable coating removal or unilateral contact, for both full-scale tests in turbomachines [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF] and smaller experimental set-up [START_REF] Batailly | Numerical-experimental confrontation in the simulation of tool/abradable material interaction[END_REF]. Numerically, the abradable elements are rigidly attached to the spline surface as depicted in Fig. 10.

General algorithm

Note: in the following, two contact scenarii are distinguished: the impeller may penetrate (1) the coating or (2) the casing. In the case the impeller penetrates the casing, it is implied that the abradable coating has been locally removed.

The algorithm embedding abradable coating removal is presented in Fig. 11 and generalizes the algorithm in Fig. 4. Once impeller/casing distances d * n+1 are computed, possible impeller/casing penetrations are dealt with following the procedure in Fig. 4. If a boundary node of the impeller penetrates the casing, the surrounding abradable coating is automatically removed: the plastic deformation of the related abradable elements is set to the initial thickness of the coating ∆ε p = e abr . Consequently, the node is quarantined for the current time step, in other words, it is excluded from the blade/abradable coating penetration detection procedure.

displacements prediction for the impeller and the casing

v * i,n+1 = M i h 2 + D i 2h -1 2M i h 2 -K i v i,n + D i 2h -M i h 2 v i,n-1 + F ext i,n+1 v * c,n+1 = M c h 2 + D c 2h -1 2M c h 2 -K c v c,n + D c 2h -M c h 2 v c,n-1 + F ext c,n+1
coordinates of the casing boundary nodes (data points)

P * n+1 = P + v * c,n+1
construction of the spline surface

Q * n+1 = f (P * n+1 )
impeller/casing distance computation for every boundary node -impeller displacements correction : For each boundary node of the impeller, the distance to the abradable coating d * ,abr n+1 is computed from d * n+1 and the local thickness of the abradable coating e abr -∆ε p . If min(d * ,abr n+1 ) < 0, the plastic constitutive law described in [START_REF] Legrand | Numerical investigation of abradable coating removal through plastic constitutive law in aircraft engine[END_REF] yields the contact forces λ abr n+1 . The computation of the contact constraints matrices and the contact matrices is identical to the procedure in section 3.2.

d * n+1 = f (Q * n+1 , v * i,n+1 )
v i,n+1 = v i,n+1 -M i h 2 + D i 2h -1 VCT abr i,n+1 λ abr n+1 -casing displacements correction : v c,n+1 = v c,n+1 -M c h 2 + D c 2h -1 VCT abr c,n+1 λ abr n+1 time update : t = t + h n = n + 1

Validation and convergence

Numerical investigations similar to the ones detailed in section 3.3 are not recalled here but identical conclusions are obtained.

Spatial convergence with respect to the number of abradable elements

Convergence is checked against the number η of abradable elements: η = 600, 3000, 6000, and 12000. This parameter has low to no impact on the computation cost of the procedure, no significant time differences were found. The radial displacement of the impeller sector 1 leading edge is pictured with the final abradable wear profiles in Fig. 12. Minor differences are visible between displacements and contact forces for various η in Figs. 12(a), 12(b), and 12(c). However, the number of abradable elements has a strong influence on the coating wear levels as depicted in Fig. 12(d). Ultimately, for η = 6000 and η = 12000, results are almost perfectly superimposed for all quantities of interest and η = 6000 abradable elements are used in the following.

0 t f = 0.9 2 3 4 -0.2 -0.1 0 0.1 0.2 time displacement (a) 0 t f = 0.9 2 3 4 -0.2 0 0.2 time displacement (b) 0 t f = 0.9 2 3 4 -1.5 -1 -0.

Localized removal of the abradable coating

As mentioned previously, the solution method is compatible with hybrid contact scenarii where the abradable coating may be locally completely removed. Such interaction scenario is encouraged by reducing the abradable coating thickness (e abr = 0.15 compared to typical values e abr = 2 or 3). The initial clearance between the impeller and the casing at the trailing edge is set to zero so that contacts and removal are activated from the very beginning of the simulation as the casing is progressively distorted. Results are pictured in Fig. 13. Fig. 13(a) indicates that the abradable coating is locally totally removed since a portion of the distances separating a few blades and the casing reaches zero. Related radial displacements at the trailing edge of each sector are plotted in Fig. 13(b). The corresponding contact forces depicted in Fig. 13(c) can be classified as follows: the ones stemming from the plastic constitutive law when the abradable coating is impacted (

) and the ones computed with Lagrange multipliers when direct unilateral contact is active after the abradable coating has been completely worn off (

). Evidence of total abradable removal is also displayed in Fig. 13(d) where the wear level reaches the initial coating thickness e abr .

Case study 1: two-lobe distortion of the casing

In this case study, a progressive loading is applied for t ∈ [0 ; t f ] along the casing circumference at the trailing edge. The deformed casing features two dominant symmetric contact areas with respect to the axis of rotation. For t > t f , the two structures are left free to interact while the impeller rotates at the angular speed Ω. The choice is made to initiate contact at the trailing edge where the casing exhibits high flexibility, far from the clamped nodes at the leading edge, shown in Fig. 1. Exemplary, the displacement and stress fields within the impeller and the casing are depicted in Figs. 14 and 15 respectively. ) and with direct contact on the casing ( ), (d) trailing edge wear level on the casing.

In agreement with the symmetry of the loading applied on the casing, the perfect symmetry of these fields is noticeable and underlines the accuracy of the employed numerical strategy. Also, while contact is initiated at the trailing edge, Fig. 14 undoubtedly indicates that the largest displacements arise in the neighbourhood of the leading edge, where the greater blade length yields higher amplitude displacements. Significant stresses are found within the impeller both in the blades and the disk sections.

Influence of friction

The influence of the friction coefficient µ is assessed considering three values: µ = 0, 0.15, and 0.3. Related results are pictured in Fig. 16. First of all, it is noticeable that the lower µ, the higher the wear level in the abradable coating, as pictured in Fig. 16(d). For the angular speed of interest Ω, µ has no effect on the casing displacements in Fig. 16(b). However, no clear conclusion can be drawn from Fig. 16(c) concerning the amplitude of the contact forces.

Sensitivity to the mechanical properties of the abradable coating

Previous work illustrated the key role of the mechanical properties of the abradable coating in the rotor/stator interactions [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF]. From the nominal values of the Young's and plastic modulii E n and K n obtained from empirical data, two extreme scenarii are considered:

(1) a lower bound E = E n 10, K = K n 10 and (2) an upper bound E = 10E n , K = 10K n . The changes in the mechanical properties have a visible effect on the displacements of both the impeller and the casing, as pictured in Figs. 17(a) and 17(b). For low Young's and plastic modulii(

), the amplitude of the boundary nodes displacements becomes lower for the impeller and higher for the casing. In that case, the removal of the coating is made easier and thus generates smaller contact forces on the impeller. On its side, the casing is less constrained in its vibration for t > t f after the external loading is not applied anymore.

The conclusions drawn in this section depend on Ω and results may be different in the vicinity of a critical speed as explained in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF].

Vibratory energy analysis

As mentioned earlier, the casing features two symmetric privileged contact areas induced by the external loading exerted on the casing. Analytical considerations exposed in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] predict that odd nodal diameters have a negligible role in the vibratory response of the system. However, the aliasing effect may generate a non negligible participation of all even nodal diameters. To confirm these predictions, the methodology described in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] could be implemented through harmonic tracking or 2D Fourier transforms, for instance. In order to present other post-processing options, our interest is directed on the impeller strain energy E i (t) and its participation within each of the six nodal diameters of the structure. If v i (t) stands for the displacement field of the impeller and K i for its stiffness matrix, then:

E i (t) = 1 2 v i (t) ⊺ K i v i (t) (20)
The displacement v i (t) can be expanded as a finite sum of free vibration modes φ p :

∃q p (t), p = 1, . . . , n, v i (t) = n p=1 q p (t)φ p (21)
After sorting the free vibration modes according to their number of nodal diameters n d , the strain energy associated to each nodal diameter takes the form:

E(n d = j, t) = 1 2 p,n d (φp)=j q p (t)φ p ⊺ K i p,n d (φp)=j q p (t)φ p (22)
The strain energy E i is depicted with all strain energies related to a given nodal diameter E(n d ) in Fig. 18. The strain energies associated with odd nodal diameters are found negligible with respect to those associated with even nodal diameters. Dominating nodal diameters 0 and 2 and a non-negligible participation of nodal diameter n d = 4 are found, which is consistent with the symmetry of the interaction scenario.

Case study 2: investigations around a critical speed

The linear interaction criterion introduced in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] can predict potential critical angular speeds. In particular, an angular critical speed is predicted for Ω = Ω 5 1 ≃ 0.49 as pictured in Fig. 19. This interaction point supposedly involves two 5-nodal diameter free vibration modes: accordingly, the casing is now distorted on a five-nodal diameter shape. Because of contact-induced stiffening, it is assumed that the interaction may actually arise for slightly higher angular speeds than Ω 5 1 . Numerical simulations are carried out from Ω = 0.48 to Ω = 0.58. The Fourier transform of the time responses obtained for each angular speed is depicted in Fig. 20. It is remarkable that no significant amplitudes of vibration is found in the vicinity of Ω 5 1 even though the interaction condition is satisfied. In order to clarify this unexpected absence of vibratory resonance, modeshapes are provided in Fig. 21. The comparison highlights a geometric incompatibility: significant deformations are localized around the leading edge for the impeller and around the trailing edge for the casing. It seems very improbable to see these two free-vibration modes dangerously interact through contact.

This very brief example shows a limitation of the linear interaction condition. As mentioned in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF], this criterion fails to predict super-or sub-harmonics related interactions. The criterion also predicts unimportant critical angular speeds as the modal shapes of the two structures are not geometrically compatible. As a result, the proposed numerical strategy advantageously complements the linear interaction condition by providing additional insight on the predicted critical angular speeds. It also brings the foundation for the development of optimized impellers robust to contact conditions.

Conclusion

The numerical strategy exposed in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] is here extended to rotor/stator interactions involving a flexible casing, full three-dimensional friction forces and the modeling of the abradable coating deposited on the casing. The related algorithms are detailed and validated. Their versatility is promising and the presented numerical strategy is applicable to a variety of rotor/stator systems potentially undergoing interaction phenomena from rubbing to modal interaction. The qualitative analysis exposed in [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] is complemented by the quantitative developments of the present work.

The proposed solution method allows for the numerical simulation of various interaction scenarii where abradable coating removal may be accounted for along the flexible casing. Accounting for the casing flexibility allows for broadening the scope of simulated interactions from rubbing phenomena to modal interaction and detecting critical angular speeds for each type of interaction. More specifically, hybrid scenarii featuring a localized yet complete removal of the coating are also explored. It is shown that the linear interaction condition, which stands as an important guideline for the design of impellers and casings, may be advantageously complemented by the numerical strategy exposed in both [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] and the present article: (1) super-and sub-harmonics related interactions that are not captured by this criterion but may be accurately predicted numerically instead and (2) the relevance of linearly predicted critical angular speed may be assessed. As shown in the last case study exposed in this article, a critical speed may be predicted for free vibration modes that are not geometrically compatible and thus, may not interact. ) and the casing ( ).

[9] F. Peyraut, J. L. Seichepine, C. Coddet, M. Hertter, Finite element modeling of abradable materials -identification of plastic

  α p internal hardening variable Γ a boundary node along a blade tip of the impeller λ contact forces vector ∆ε p plastic deformation of an abradable element η number of abradable elements along the casing circumference contact matrix of the time-stepping procedure M mass matrix n out outward normal vector to the contact surface P data points of the spline surface: boundary nodes of the casing model Q control points of the spline surface Q bo control points on the edges of the spline surface Q int control points inside the spline surface t tangential vector to the contact surface VCT contact constraints matrix in radial and tangential directions VC contact constraints matrix in the radial direction v displacement field E strain energy Q i,j control points related to the patch S i,j S i,j a patch of the spline surface µ friction coefficient φ free vibration mode of the impeller ε NR Newton-Raphson iterative procedure stop criterion

Figure 1 .

 1 Figure 1.Surrounding casing of the impeller used in[START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF]

  z

Figure 2 .

 2 Figure 2. Free vibration modes of the casing: (a) one nodal diameter, (b) two nodal diameters, (c) three nodal diameters, (d) four nodal diameters.The reader may refer to[START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: analysis for rigid casings[END_REF] for more details regarding the notion of nodal diameters.

Figure 3 .

 3 Figure 3. Partial representation of the spline attached to the casing contact surface. Boundary nodes of the first sector are denoted P(1, i), i = 1, . . . , M (corresponding control points Q are not plotted)

  . 1 4 1 0 ... 4 16 4 0 ... 1 4 1 0 ... 0 ... 0 1 4 1 0 ... 4 16 4 0 ... 1 4 1 0 0 ... 0 0 1 4 1 0 ... 4 16 4 0 ... 1 4 1

  1 : detection of penetrations min(d * n+1 ) < 0 -contact and contact constraints matrices of the casing : MC c,n+1 , VC c,n+1 and VCT c,n+1

Figure 4 .

 4 Figure 4. Time integration algorithm with unilateral contact conditions

Figure 5 .

 5 Figure 5. Newton-Raphson iterative procedure for the identification of the closest point on the spline from the bladed disk boundary node Γ * ,k n+1 .

Figure 6 .

 6 Figure 6. Bladed disk/casing distances algorithm

Figure 7 .

 7 Figure 7. Determination of the t k n+1 unit vector along which acts the tangential friction force

Figure 8 .

 8 Figure 8. Sensitivity to stopping criterion in Newton-Raphson algorithm: εNR = 10 -5 ( ) and εNR = 10 -10 ( ): (a) radial displacement of the leading edge of impeller sector 2, (b) radial displacement of a boundary node of the casing

Figure 9 .

 9 Figure 9. Sensitivity to friction: µ = 0 ( ) and µ = 0.01 ( ): (a) tangential displacement of the leading edge of impeller sector 2, (b) radial displacement of a boundary node of the casing

Figure 10 .

 10 Figure 10. Impeller sector with abradable coating (yellow) and supporting spline surface (green).

case 1 :

 1 detection of penetrations min(d * n+1 ) < 0 →impacted abradable elements are removed : ∆ε p = e abr →proceed with steps of Fig. 4, case 1 → impacting nodes of the impeller are quarantined case 2 : no detected penetration min(d * n+1 ) ≥ 0 →proceed with steps of Fig. 4, case 2 for the impeller boundary nodes that are not quarantined : if penetrations are detected in the abradable coating min(d *,abr n+1 ) < 0 -contact force computation λ abr n+1 for each element -parameters ∆ε p and α p are updated -contact constraints matrices computation VCT abr c,n+1 et VCT abr i,n+1

Figure 11 .

 11 Figure 11. Time integration algorithm contact treatment and abradable coating management.

Figure 12 .

 12 Figure 12. Convergence analysis with respect to the number of abradable elements η: η = 600 ( ), η = 3000 ( ), η = 6000 ( ) and η = 12000 ( ), (a) radial displacement of the impeller sector 1 leading edge, (b) radial displacement close to the casing trailing edge, (c) radial contact force on the impeller sector 1 trailing edge, (d) abradable removal profile.

Figure 13 .

 13 Figure 13. Locally full removal of the abradable coating: (a) distances between the trailing edge of each sector of the impeller and the casing, (b) radial displacements of the trailing edge of each sector, (c) radial contact forces on the trailing edge of each sector when impacting the abradable coating () and with direct contact on the casing ( ), (d) trailing edge wear level on the casing.

Figure 14 .

 14 Figure 14. Impeller during the interaction: (a) displacement field at t = 0.3 < t f , (b) Von Mises stress field at t = 0.3 < t f , (c) displacement field at t = 1.5 > t f , (d) Von Mises stress field at t = 1.5 > t f . Dsiplacements are magnified for better readability.

  z

Figure 15 .

 15 Figure 15. Casing during the interaction: (a) displacement field at t = 0.3 < t f , (b) Von Mises stress field at t = 0.3 < t f , (c) displacement field at t = 1.5 > t f , (d) Von Mises stress field at t = 1.5 > t f . Displacements are magnified for better readability.

Figure 16 .

 16 Figure 16. Influence of the friction coefficient: µ = 0 ( ), µ = 0.15 ( ) and µ = 0.3 ( ): (a) radial displacement of the impeller sector 1 trailing edge, (b) radial displacement at the casing trailing edge, (c) radial contact force on the impeller sector 1 trailing edge, (d) abradable removal profile.

Figure 17 .Figure 18 .Figure 19 .

 171819 Figure 17. Sensitivity to the abradable coating properties: (En 10, Kn 10) ( ), (En, Kn) ( ) and (10En, 10Kn) ( ): (a) radial displacement of the impeller sector 1 trailing edge, (b) radial displacement at the casing trailing edge, (c) radial contact force on the impeller sector 1 trailing edge, (d) abradable removal profile.

5 1Figure 20 . 1 z

 5201 Figure 20. Fourier transform of the impeller sector 1 radial displacement around Ω 5 1

Figure 21 .

 21 Figure 21. Modal shapes of the investigated 5-nodal diameter free vibration modes.

Table 1 .

 1 Differences between present numerical developments and reference[START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF]

	blade/casing distances computation hypotheses on the casing displacement friction	present work iterative nonlinear procedure none full 3D	reference [27] simplified procedure radial only tangential direction

-Laboratoire d'Analyse Vibratoire et d'Acoustique, École Polytechnique de Montr éal,

Boulevard Édouard Montpetit, H3T 1J4 2 -Laboratoire de Dynamique des Structures et Vibrations, Universit é McGill, 817 Sherbrooke Ouest, b âtiment McConnell, pi èce 122, H3A 0C3, Montr éal, Qu ébec, Canada

Acknowledgement

Thanks go to Turbomeca and Snecma for their technical and financial support.