Khadija Masrour

Alain Batailly Florence

Jean-Philippe Ponthot

Marco Morandini

Alternative approach to model wear evolution of a blade undergoing contact interactions with the surrounding casing

Also, I would like to express my deep gratitude to my parents and my sisters Loubna, Botaina,

One of the most important requirements in modern engines design is the minimization of the clearance between rotating and static components. It has been shown in several studies that this clearance minimization increases the aerodynamic efficiency of the internal flow leading to significative improvement of the engine performances. During engine operation due to several phenomena such as thermal dilatations of the parts, misalignment of the shaft due to centrifugal effects, gravity effects during manouvers and unbalanced parts, contact interactions between the rotating blades and the abradable seal of the surrounding casing can occur. It has been observed that these contact interactions can lead to undesirable consequences that can affect significantly the performances of the engine. This thesis is concerned with the modeling of wear phenomena on an engine blade proposing an alternative approach for the geometry update. The basic idea of this approach is instead of re-meshing or eliminating at each iteration the worn out elements such as it is performed in the traditional remeshing and killing element methods for wear treatment problems, is to interpolate by a polynomial function the structural matrices in correspondence of a small subset of interpolating points that correspond to different amounts of wear and recover the structural matrices for all other configurations by the build polynomial function. The advantage of this technique is the simplicity of the operations involved at each iteration. Indeed at each iteration a simple algebraic sum is performed on the constant interpolating matrices to recover the solution for all the blade configurations due to wear evolution. Thus, there is no need to keep memory of the whole wear evolution history and the simplicity of the involved operations leads to significant decrease in the computational effort and memory storage. In the performed analysis this proposed approach is fully validated in case of a 1D structure. To perform such analysis a FEM model of a clamped beam has been exploited. In the test conditions the beam undergoes to contact interactions with a rigid moving wall. The wall is moving both in normal direction penetrating the blade tip and transversally with a constant upward velocity inducing blade vibrations in both the directions. The contact interactions are modelled with both Lagrange multipliers and Penalty methods. In a second moment the wear is added to the analysis. Before implementing the proposed approach, traditional re-meshing techniques have been applied to the problem in order to recover reliable results to use as term of comparison and validation for the proposed approach results. Accuracy and reliability of the proposed approach is demonstrated by different test case results involving different conditions, such as soft contact interactions and severe interactions. A focus have been reserved also to the discussion of the important savings in the computational cost registered. Finally, some suggestions and ideas for future developments and extension of such technique to 2D structures are given. All the above mentioned algorithms and simulations have been implemented exploiting Python language.

i Dedicato a tutti coloro che hanno contribuito direttamente o indirettamente alla realizzazione di questo lavoro. Alla mia famiglia, ai miei professori, ai miei colleghi e ai miei amici. È grazie al contributo di ognuno di loro che ho potuto raggiungere questo traguardo tanto significativo per me, nonchè un grande sogno.

Exactly one year ago, August 2017, I started searching for an internship to complete my Master's degree. Six months later, February 2018, I started this amazing journey towards an unknown destination but that soon became a new home, Montreal.

There I started my internship and along all the time at LAVA laboratory the guidance and advices of Dr. Alain Batailly and Dr. Florence Nyssen were fundamental for me in this journey, as well as their support. And from here I would like to express my sincere gratitude to both of them.

Next to this, my sincere thanks also go to my supervisors at both Université de Liège and Politecnico di Milano, prof. Jean-Philippe Ponthot and Prof. Marco Morandini, as well as all my professors in both the universities, my colleagues and my friends. It has been a unique path for me which principle gains have been double experiences, double emotions, double friendships and double memories that I will conserve forever. It is highlighted the advantage of having abradable linings where the induced clearance increase is locally concentrated, whereas, in case no abradable linings are used the induced gap increase interests the whole engine section [START_REF] Fois | Investigation and caharacterizaton of the wear mechanisms of abradable compressor linings[END_REF]. . . . . . . . . . . . . . . . . . . . . . 2.4 Shows the most used abradable materials in modern engines according to the compressor stage and thus to the environment temperature conditions [START_REF] Fois | Investigation and caharacterizaton of the wear mechanisms of abradable compressor linings[END_REF]. . . . . . . . . . 2.5 Flowchart of Newmark's method [START_REF]Lecture notes of Theory of vibrations course[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Stability Newmark plane [START_REF]Lecture notes of Theory of vibrations course[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Recovered time integration schemes by Newmark's method by different values of γ and β, with relative stability limit, amplitude error and periodicity error [START_REF]Lecture notes of Theory of vibrations course[END_REF]. . . . . . 2.8 Graphical comparison of the contact constraints fulfillment by Lagrange multipliers and Penalty method. Lagrange method does not allow any penetration between the two bodies, while Penalty does not fulfill perfectly the contact constraint allowing a certain amount of penetration [START_REF] Harmon | Robust, Efficient, and Accurate Contact Algorithms[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Are shown the contact conditions in normal direction. When the gap function is greater than zero (no contact), the normal contact pressure is zero. While when the gap is equal to zero (contact) a non-zero normal contact pressure occurs. No negative gap function is present because no bodies penetration is allowed [START_REF] Gitterle | A dual mortar formulation for nite deformation frictional contact problems including wear and thermal coupling[END_REF]. . . . . . . . . 2.10 Frictional contact conditions in tangential direction according to Coulomb's model [START_REF] Gitterle | A dual mortar formulation for nite deformation frictional contact problems including wear and thermal coupling[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Penalty regularization of normal contact conditions allowing a certain amount of penetration in fucntion of the normal contact pressure. The intensity of this penetration is expressed by the Penalty coefficient (C N ), higher is this value and lower is the allowed penetration and vice versa [START_REF] Gitterle | A dual mortar formulation for nite deformation frictional contact problems including wear and thermal coupling[END_REF]. . . . . . . . . . . . . . . . . . . . . . . 2.12 Abrasive wear of ductile material [START_REF] Kato | Wear Mechanisms[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Adhesion wear mechanism as result of rubbing of two surfaces [START_REF]Substances and Technologies, Knowledge source on Materials Engineering[END_REF]. . . . . . . . . . 2.14 Fatigue wear mechanism as result of rubbing of two surfaces [START_REF]Substances and Technologies, Knowledge source on Materials Engineering[END_REF]. . . . . . . . . . . viii LIST OF FIGURES ix 2.15 Single-layer re-meshing technique illustration [START_REF] Rezaei | Adaptive finite element simulation of wear evolution in radial sliding bearings[END_REF]. a) Initial configuration with selected the elements of the contact zone that will be interested by re-meshing. b) Remeshing steps during the wear simulation involving single layer remeshing at each iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Body configuration at the j th and (j +1) th configuration, showing its nodal displacements to accommodate the total amount of wear h tot,j+1 according to the equally displaced nodes technique [START_REF] Benabdallah | Finite element simulation of the wear of polyoxymethylene in pin-on-disc configuration[END_REF]. a is the critical element height H cr . . . . . . . . . . 2.17 Mesh distortions that appeared after few sliding cycles of a pin on a disc in which the wear has been accommodated following the equally displaced nodes technique [START_REF] Benabdallah | Finite element simulation of the wear of polyoxymethylene in pin-on-disc configuration[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 Body configuration in the unworn (a) and worn case (b). The nodes in the worn case are displaced according to the linear displaced nodes technique [START_REF] Tang | A multilayer nodes update method in FEM simulation of large depth fretting wear[END_REF]. . . . . . . . . . 2.20 Mapping of the physical coordinates of a quadrilateral element into natural coordinates. 2. [START_REF] Hirani | Wear Mechanisms-2[END_REF] Flowchart of elimination elements approach combined with re-meshing technique for wear simulation [START_REF] Békési | Modelling Friction and Abrasive Wear of Elastomers[END_REF]. At each iteration the contact calculations are performed allowing to evaluate the stress distribution. If the damage criterion for the analyzed element is fulfilled then the element is deactivated/killed otherwise a re-meshing is performed moving the element nodes. . . . . . . . . . . . . . . . . . . . . . . . . . 2.21 Shows the first four product rules for Gauss Integration with equal number of points along the two directions [START_REF] Nikishkov | Introduciton to the Finite ELement Method, Isoparametric transformations[END_REF] which starts with a clearance of 10µm and reaches a maximum penetration of 40µm into the beam. Once the wall leaves free the beam ( ∼ = 15.1ms), this one starts vibrating vertically with a frequency of approximately 840Hz, which corresponds to the its first frequency. These results are computed by Lagrange multipliers method applying a frictional coefficient of 0.3. . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Sensitivity of beam tip vertical displacement to the frictional coefficient. Higher is the frictional coefficient and higher is the frictional force, thus, higher is the energy with which the beam is pushed upward leading to higher vertical displacements. The results are obtained applying Lagrange multipliers method with contact conditions described in section3. 3.9 Comparison of the beam's tip axial displacement computed by Lagrange multipliers method and Penalty method. With the selected value of the Penalty parameter the Penalty method results overlap perfectly those given by Lagrange multipliers method. 48 3.10 Comparison of the beam's tip vertical displacement computed by Lagrange multipliers method and Penalty method. With the selected value of the Penalty parameter the Penalty method results overlap perfectly those given by Lagrange multipliers method. 49 3.11 Meshed jet-engine compressor blade in which only a subset of contact nodes is retained in CB reduced model while all other nodes are condensed [START_REF] Batailly | Redesign of a High-Pressure Compressor Blade Accounting for Nonlinear Structural Interactions[END_REF] In case of higher value of the wear coefficient (red line) the u penetration of the beam tip is lower because a higher value of material is removed from the beam. And also, since the beam length is decreased, the contact interaction starts later and ends earlier with respect to the other two simulations. . . . . . . . . . . . . . . . . . 76 4.25 Vertical displacement of the beam tip in function of different values of the wear coefficient. In case of higher value of the wear coefficient (red line), as shown in Fig. 4.24, the u penetration of the beam tip is lower because a higher value of material is removed from the beam, this results in lower vertical vibrations of the beam tip, that start later and end earlier in phase with the u displacement.

Chapter 1 Introduction

Modern engines are subjected to continuous and radical improvements in the design to match the basic requirement of operating in all atmosphere conditions, even the extreme ones, ensuring safety and reliability with an optimization of the performances.

The basic functionality concept of a jet engine is the acceleration of a sucked fluid mass from the external environment. This acceleration provide a thrust to the aircraft in the opposite direction propelling it forward. The air is sucked through an in front fan and its pressure is raised by a compressor, the compressed air is then guided to a combustion chamber in which fuel is sprayed in. An electrical spark initialize the mixture burning and the hot gases expands through a turbine and finally are ejected through a nozzle generating a forward thrust that propels the engine and thus the aircraft. A schematic illustration of a jet engine components is given in Fig. 1.1. The aim in an engine design is to keep as high as possible the ratio of such generated work with respect to the consumed fuel. Thus, the operating conditions of the mechanical parts, of the combustion chamber, of the injection and ejection parts and the conditions of the passage of the internal flow through each component reveals to be crucial for the designer to improve such operating performances.

One of the modern engine requirements and which will be a basic concept of such master thesis research is the reduction of the clearance between the static and rotor components of the engine, i.e. between the blade tips and the surrounding casing. From a mechanical point of view such clearance is desired to be as high as possible to avoid any contact interaction between rotor and static components which could lead to stress accumulation and thus to parts failure, while, from an aerodynamic point of view, it has been shown in different studies that such clearance has to be reduced as much as possible to lead to significant increase in the engine efficiency [START_REF] El-Abbasi | Stability and patch test performance of contact discretizations and a new solution algorithm[END_REF]. Under this basic requirement different experimental researches have been conducted to simulate the effects of such clearance decrease. But experimental investigations revealed to be time consuming and expensive, leading to the development of numerous numerical techniques. These studies can be classified in two main categories, models that predict the contact interactions and neglect the wear evolution in the structure, and models that take into account wear. The first one can be used when the wear is very limited on the parts otherwise it has to be taken into account.

But the main restriction for the developed techniques that take into account wear is how to model the abradable material and the high computational cost of the overall simulation limiting it to few blade revolutions. This is the main reason behind the proposition of the polynomial interpolation technique in this thesis.

Objectives of the research

The objective of this thesis research is to introduce and initialize the validation of an alternative technique for wear treatment in contact interactions of rotor blades with the surrounding casing as consequence of the clearance minimization. Due to these contact interactions material layers are continuously removed from the blade tips leading to its geometrical change which has to be updated at each iteration of the dynamic simulation.

Traditional wear treatment techniques that deals with this continuous update of the geometric configuration can be classified in two categories: re-meshing techniques and element killing techniques.

Both the approaches require perpetual intervention on the structural matrices of the part to update its geometry change leading to high computational effort and memory storage. Thus, the motivation behind this research is to propose an alternative approach of treating with wear problems that can be eventually extended to several applications and not limiting to aerospace field. The basic idea is instead of computing directly new structural matrices at each time iteration, compute them just in a small subset of geometrical configurations and perform an interpolation of such "exact" solutions to recover the structural matrices in all other configurations by a simple polynomial interpolation. Polynomial interpolation presents the important advantage of involving simple algebraic operations that will lead to important savings in the computational time.

Another important objective of such proposed approach is to apply a reduction technique to the model and interpolate the reduced structural matrices of the interpolating configurations leading to further decrease in the memory storage since, as known, reduction model techniques reduce significantly the model dimensions ensuring good accuracy of the results.

The proposed approach has to be first validated fully validated for a 1D structure. Once proved that such idea can be exploited to give accurate results in 1D case, then it can be thought to be extended to 2D structures and eventually to 3D structures.

In this master research the first point will be covered validating the approach for a 1D structure proving that such proposed technique can be reliable. Therefore, the second point will be initialized giving some inputs for the solution but further research is under development at LAVA laboratory at Polytechnique of Montreal. To reach such goal of validating the proposed approach for 1D structures some essential steps have to be covered to accurately model and analyze the vibratory response of the structure:

1. develop a simplified 1D model of the blade modeling it as a simple clamped beam exploiting FE method;

2. implement a contact algorithm to simulate the contact interactions of the beam with a second body. This second body has been selected to be a rigid moving wall.

3. select important parameters such as bodies penetration and frictional parameters to prepare a basis for the next step;

4. add the option of wear in the dynamic simulation;

5. implement traditional geometrical update algorithms such as re-meshing techniques;

6. implement the proposed polynomial interpolation technique and compare its results with those provided by the implemented re-meshing techniques algorithms that will be used as a term of validation of such technique;

7. apply reduction model technique to the studied system;

8. analyze the provided results and discuss the reliability of the polynomial interpolation technique;

9. provide some ideas for the extension of the proposed approach to 2D structures.

All the enumerated goals will be implemented exploiting Python language for the development of the algorithms.

Delimitations

As the aim of the thesis is to focus on the wear evolution and treatment of an engine blade, then the surrounding casing is modeled as a rigid component. A parallel research at LAVA laboratory is conducted dealing with abradable casing seals while modeling the blade as a rigid component. Once both the models are fully completed, their coupling will be performed. A second limitation is on the blade loads. Since the focus is on a mechanical vibratory response, a CFD analysis will not be performed. Hence, aerodynamical loads which include lift and drag force are not included. Furthermore, according to analytical calculations, the values of those forces are minimum with respect to mechanical loads, thus, they are neglected along the thesis. Finally, another delimitation is made restricting the analysis on 1D structures because only fully validating the proposed approach for 1D structures that can be thought to be extended to more complex models. However, a final chapter is dedicated to the introduction of such approach also to 2D structures with the proposition of some inputs for the initialization of the analysis and the further development.

Structure of the thesis

This master thesis research is divided into two main parts. The first part treats the wear modeling for 1D structures, in particular a clamped beam modeled by Finite Element Method (FEM). In this part the proposed approach of Polynomial interpolation technique is discussed and validated. In the second part the extension of such approach is initialized to 2D structures giving some inputs for future and further developments.

The layout of the thesis is the following:

• Chapter2 gives a global overview of theoretical concepts needed for the development of the project such as dynamic time integration techniques in particular those exploited for contact problems, available contact algorithms such as Lagrange multipliers, Penalty method and Augmented Lagrangian method, wear models and finally widely used geometry update techniques to adapt the FEM model to the wear evolution in the structure.

I Wear modeling in 1D structures:

• Chapter3 presents the model used and test conditions to simulate and analyze the desired frictional contact interactions. Both Lagrange multipliers and Penalty methods are implemented and relative results are presented.

• Chapter4 introduces wear in the problem using Archard's law to model it. Traditional wear treatment techniques as local and global re-meshing are applied and relative results are presented. Then the proposed approach of Polynomial interpolation of the structural matrices is introduced and its results are compared with those given by global and local re-meshing.

II Wear modeling in 2D structures:

• Chapter5 presents an attempt of extending the proposed approach to 2D structures giving some ideas for future developments.

• Chapter6 a summary of the founded results is presented with relative suggestions for future works.

Chapter 2

Literature review

In this chapter will be presented a generic theoretical overview of the concepts needed for the full development of this thesis project.

A first section is dedicated to the tip clearance minimization requirement in engines design which constitute the background for this performed research. One of the solutions adopted by designers to overcome or at least decrease the drastic effects that could have the contact events that occur between rotating and static components is the application of abradable linings on the internal casing surface and thus the interactions occur between the blades tip and the abradable linings. A synthesis of the properties and the abradable materials used in aerospace engines is given. Then a focus is deserved to time integration techniques and in particular to Newmark's method as it is one of the methods used to solve dynamic simulations. Some contact treatment methods and model reduction techniques are presented too.

The research is then extended to wear giving first a summary of the kind of wear phenomena that exist and how they are modeled, then a presentation of the available techniques exploited by analysts to adapt the FEM model to the evolving wear in the structure. Finally, a section is dedicated to 2D structures and how they are modeled in the FE method.

Engines design requirements

In an engine design different requirements have to be accomplished to ensure optimal efficiency and performances. Some of these requirements is a reduction in the consumed fuel, reduction in the emitted gases into the atmosphere, lower demand of working level of pressure and temperature and a maximization of the service life [START_REF] El-Abbasi | Stability and patch test performance of contact discretizations and a new solution algorithm[END_REF].

One of the solutions adopted to fulfill such requirements and increase the efficiency of a jet engine is the minimization of the distance between the rotor and stationary components. This parameter is normally called Tip Clearance in the literature. But this clearance minimization can lead to numerous contact interactions between blades and the surrounding casing. The criticity of these interactions is decreased applying abradable linings on the internal surface of the casing.

Tip clearance design

From a mechanical point of view an increase of the tip clearance is desired because it decreases the possibility of interaction between the engine components and thus decreases the stress accumulation that could lead to failure of the parts. However, such clearances are imposed from an aerodynamical point of view since it allows an increase of aerodynamic performance such has been proven in the experimental investigation conducted in Lewis Research Center [START_REF] Holeski | Effect of rotor tip clearance on the performance of a 5-inch single stage axial flow turbine[END_REF].

In this experiment the principal aim was to determine the effects of increasing rotor tip clearance on the performance of a 5inch (13cm) axial-flow turbine. The clearances studied range from 1.2% to 8.0% of passage height. The fluid used is Argon and the turbine was operated at constant design pressure ratio while the speed has been varied in a certain range. The results obtained by this experiment are presented in Fig. 2.1, in which are shown both the static and total efficiency in function of the blade-jet speed ratio and the tip clearance. The figure shows an approximately linear decrease in both the efficiencies with increasing rotor tip clearance. Indeed, for the analyzed turbine, a 1.0% increase in clearance to height ratio resulted in approximately a 3.0% reduction in efficiency. A comparison of the obtained results has been conducted with those of two reference reaction turbines which has indicated similar effects of clearance on efficiency, proving that for an optimal engine's design, a clearance minimisation is a must.

Abradable linings

The critical contact interactions that the tip clearance minimization could lead to between the rotating and stationary parts have to be taken into account as part of normal engine running conditions. One of the solutions adopted to decrease the criticity of these interactions is to apply abradable linings to the internal surface of the casing (Fig. 2.2). The reason behind this choice is the introduction of surface casing material that can wear more easily than the blade, leading to a local increase of the tip clearance and not global increase. As FOIS showed in his research [START_REF] Fois | Investigation and caharacterizaton of the wear mechanisms of abradable compressor linings[END_REF], if the blade undergoes wear, the induced variation of clearance will interest the whole circumference of the engine as the blade is a rotating component. While if the wear interests the lining, this increase in the clearance will be a local increase. This concept is shown in Fig. 2.2 where it is possible to see the advantage of having abradable linings in which the induced gap due to the contact interactions between the components is locally localized, while in case of no abradable linings the induced gap increase interests the whole engine section. The figure highlights three phases: prerubbing, interaction and clearance after rub interaction comparing the two cases in which an abradable lining is used and the second case in which no abradable coating is used. It is highlighted the advantage of having abradable linings where the induced clearance increase is locally concentrated, whereas, in case no abradable linings are used the induced gap increase interests the whole engine section [START_REF] Fois | Investigation and caharacterizaton of the wear mechanisms of abradable compressor linings[END_REF].

Abradable materials

The abradable materials used as casing linings have to present a good abradability property and good resistance to erosion. Among the materials that fulfill these requirements composites represents an optimum choice. They are exploited because they are soft enough to wear more easily than the blades and at the same time resistant enough to the damages that could occur as consequence of these intense contact interactions. Moreover, they are sufficiently resilient to withstand the high environment engine temperatures. Typically the composite materials used as abradable linings are made by the following structural components [START_REF] Fois | Investigation and caharacterizaton of the wear mechanisms of abradable compressor linings[END_REF]:

• metal matrix: responsible of the strength of the material and ensures a good resistance to oxidation;

• solid lubricant: shears easily absorbing the loads and thus avoiding their transfer to the blade. They ensure also that the debris are sufficiently small to not affect the engine operability conditions;

• porosity.

An overview of the materials used for the abradable linings in an aerospace engine, according to the compressor stage, is given in Fig.

Time integration techniques

The dynamic response of a structural body is normally described by a second ordinary differential equation. Different techniques are exploited to compute this response in function of the time, one of these is a direct integration of the equation that in case of a damped system assumes this form:

Mq(t) + C q(t) + Kq(t) = p(t) (2.2.1)
where q(t) is the vector of nodal displacements, M is a positive definite mass matrix, C is a nonnegative definite damping matrix, K is a non-negative definite stiffness matrix and p(t) is the vector of external dynamic loads.

There are two general classifications of direct time integration methods:

• Explicit schemes;

• Implicit schemes.

In explicit methods the solution at time step t i+1 can be determined in closed form from displacements, velocities and accelerations at t i , contrary to implicit schemes in which the solution at a certain time step depends on itself. Therefore, the equations involved in implicit schemes are more complex and require more computational effort in each solution step. The reason behind using implicit schemes is due to stability issues. Implicit schemes are indeed unconditionally stable, which means that for any arbitrarily large values of the time step the solution remains well behaved. This situation never occurs with explicit methods, which are mostly conditionally stable. Both these explicit and implicit schemes can be recovered by Newmark's method according to the selected values for the parameters γ and β.

An illustration of the method is given following the explanation given by GERADIN and RIXEN in [5]. Considering Taylor's expansion of a generic function f :

f (t n +h) = f (t n )+hf (t n )+ h 2 2 f (t n )+...+ h s s f s (t n )+ 1 s! tn+h tn f n+1 (τ )(t n +h-τ ) s dτ (2.2.2)
The equivalent Taylor expansion of the displacements and velocities at time t n is then:

q n+1 = q n + h qn + t n+1 tn q(τ )(t n+1 -τ )dτ (2.2.3) qn+1 = qn + tn+1 tn q(τ )dτ (2.2.4)
Taylor expansions of qn and qn+1 around τ ∈ [t n , t n+1 ]:

qn = q(τ ) + q (3) (τ )(t n -τ ) + q (4) (τ ) (t n -τ ) 2 2 + ... (2.2.5) qn+1 = q(τ ) + q (3) (τ )(t n+1 -τ ) + q (4) (τ ) (t n+1 -τ ) 2 2 + ... (2.2.6)
Multiplying eq.(2.2.5) by (1 -γ) and eq.( 2.2.6) by γ and combining them it is possible to extract q(τ ):

q(τ ) = (1 -γ)q n + γq n+1 + q (3) (τ )(τ -hγ -t n ) + o(h 2 q (4) ) (2.2.7)
Now, multiplying eq.(2.2.5) by (1 -2β) and eq.(2.2.6) by 2β and combining them the following expression is obtained:

q(τ ) = (1 -2β)q n + 2βq n+1 + q (3) (τ )(τ -2hβ -t n ) + o(h 2 q (4) ) (2.2.8)
Substituting the first expression of q(τ ) into t n+1 tn q(τ )dτ and the second expression of q(τ ) into t n+1 tn q(τ )(t n+1 -τ )dτ the following expressions are recovered:

t n+1 tn q(τ )dτ = (1 -γ)hq n + γhq n+1 + r n (2.2.9) t n+1 tn q(t n+1 -τ )(τ )dτ = ( 1 2 -β)h 2 qn + βh 2 qn+1 + r n (2.2.10)
where

r n = ( 1 2 -γ)h 2 q (3) (τ ) + o(h 2 q (4) ) (2.2.11) r' n = ( 1 6 -β)h 3 q (3) (τ ) + o(h 4 q (4) ) (2.2.12) and t n < τ < t n+1
So it is possible to write: qn+1 = qn + (1 -γ)hq n + γhq n+1 (2.2.13)

q n+1 = q n + h qn + h 2 ( 1 2 -β)q n + βh 2 qn+1 (2.2.14)
where γ and β are parameters associated with the quadrature scheme. A summarized scheme of Newmark's method is shown in Fig. 2.5.

Figure 2.5: Flowchart of Newmark's method [START_REF]Lecture notes of Theory of vibrations course[END_REF].

It is possible to notice from the scheme that by this method there are no more integration operators and they have been substituted by a discrete sum. However this scheme introduces some errors, so it is important to check if it is stable or not. To have a convergence of the method it is necessary to ensure both consistency and stability. An integration scheme is said to be consistent if:

lim h→0 q n+1 -q n h = q(t n ) (2.2.15)
It is possible to show that effectively Newmark's method is consistent [5].

Regarding the stability, a time-integration scheme can be considered stable if it exists a time step h 0 > 0 in such way that for any positive h lower than h 0 (h ∈ [0, h 0 ]), a finite variation of the state vector at time t n induces only a non-increasing variation of the state-vector q n+j calculated at a subsequent time t n+j [START_REF]Lecture notes of Theory of vibrations course[END_REF].

The stability conditions for Newmark's scheme, obtained analyzing the spectral properties of its amplification matrix, are shown in Fig. It is possible to notice from the graph that the algorithm is conditionally stable if γ ≥ 1 2 and it is unconditionally stable if furthermore β

≥ 1 4 (γ + 1 2 ) 2 .
But good stability of the solution does not ensure a good accuracy and this is what is summarized in Fig. 2.7.

Figure 2.7: Recovered time integration schemes by Newmark's method by different values of γ and β, with relative stability limit, amplitude error and periodicity error [START_REF]Lecture notes of Theory of vibrations course[END_REF].

It is possible to notice that for implicit schemes the stability limit is infinite, while for the other schemes such stability is bounded to a certain limit.

Reduction techniques

In many engineering applications, the model of the analyzed dynamic system can be rather complex. This complexity results in high order and complicated mathematical models, which accurately represent the problem, but could be unsuitable for the desired application. Under the limitations of computational, accuracy and storage capabilities, the necessity of developing simplified models that conserve the main features of the original dynamic systems have evolved during the years leading to important reduction methods available nowadays [START_REF] Saric | Dynamic Model Reduction: An Overview of Available Techniques with Application to Power Systems[END_REF].

Model reduction techniques are applied on the original complex high-order system, to find an equivalent simpler low-order system that approximates the first one conserving its features. Thus, the original system is replaced by a much smaller one, and if this approximated one describes the original one within a prescribed tolerance, then it will be the only one that is needed to be simulated, leading to a significant decrease in the computational time and memory storage [START_REF] Besselinka | A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control[END_REF].

The main characteristics of an order reduction model could be summarized as follows [START_REF] Saric | Dynamic Model Reduction: An Overview of Available Techniques with Application to Power Systems[END_REF]:

• "the approximation error is small, and there exists a global error bound;

• system properties (like stability or passivity) are preserved;

• the procedure is computationally stable and efficient".

The first reduction techniques development started in the structural dynamics field, due to the importance of structures dynamic analysis. The basic concept of reduction methods in structural dynamics field consists of building a subspace R of dimension n × m (m < n) so that the solution of the reduced model (q) can be approximated as:

q ≈ Rq (2.3.1)
Considering the variational approach [START_REF]Lecture notes of Theory of vibrations course[END_REF]:

δ(T max -V max ) = 0 (2.3.2)
in which T and V are respectively the kinetic and potential energy of the system, which expressions can be replaced leading to:

δ 1 2 q T Kq - 1 ω 2 q T Mq = 0 (2.3.3)
where K is the system stiffness matrix and M is the system mass matrix. Substituting 2.3.1 into 2.3.4, the following relation is recovered:

δ 1 2 qT R T KRq - 1 ω 2 qT R T MRq = 0 (2.3.4)
The reduced eigenvalue problem is then:

( K -ω 2 M)q = 0 (2.3.5) with: M = R T MR (2.3.6) K = R T KR (2.

3.7)

There exist several reduction techniques and they differ substantially by the choice of the matrix R.

In the following sections two model reduction methods will be described, which differ indeed by the choice of the matrix R and will be used in the studied case:

• Guyan-Irons (GI) static condensation method [Guyan 1968] [5];

• Craig-Bampton (CB) reduction method [Craig 1968] [5].

Guyon-Irons static condensation method

Since the reduction model techniques consist of decreasing the system dimensions, then it is opportune to partition the degrees of freedom of the system into remaining ones (kept), denoted by q R , and condensed ones (eliminated), denoted by q C . The undamped equations of motion are:

Mq + Kq = p (2.3.8)
The system with the reorganized degrees of freedom can be rewritten as:

M RR M RC M CR M CC qR qC + K RR K RC K CR K CC q R q C = p R p C (2.3.9)
The corresponding eigenvalue problem is:

(K -ω 2 M)q = 0 (2.3.10)
which can be rewritten as:

Kq = F (2.3.11)
with F the vector of inertia forces:

F = ω 2 Mq (2.3.12)
The system is then:

K RR K RC K CR K CC q R q C = p R p C (2.3.13)
If the masses in the condensed system are zero or negligible, then it is possible to neglect the condensed inertia forces:

F C ≈ 0 (2.3.14)
Under this assumption, it is possible to obtain this relation from the system:

q C = -K -1 CC K CR q R (2.3.15)
It is then possible to write:

q = q R q C = q R -K -1 CC K CR q R = I -K -1 CC K CR q R = Rq R (2.3.16)
The Guyan-Irons reduction matrix is then:

R = I -K -1 CC K CR = I R CR (2.3.17)
Practically, to avoid the inverse operation that appears in R, the second member is computed solving the following static problem:

K CC R CR = -K CR (2.3.18)

Craig-Bampton reduction method

Craig-Bampton method can be considered as an extension of Guyan-Irons method by giving a dynamic description of the condensed degrees of freedom. The transformation matrix indeed is an extension of Guyan-Irons transformation matrix where part of the condensed modes are added to it.

Thus, Craig-Bampton transformation matrix is:

R = I 0 -K -1 CC K CR XC (2.3.19)
where XC is a subset (m) of the eigenmodes associated to the following eigenvalue problem:

K CC q C -ωM CC q C = 0 (2.3.20)

Contact treatment

The development of methods for solving problems with moving interfaces is one of the major aims of continuum mechanics. Contact problems are inherently nonlinear since the contact area is a priori unknown. The main difficulty lies in the constitutive laws of contact and friction expressed by non-smooth multivalued force-displacement relationships [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF]. From a mechanical engineering viewpoint, three groups of methods are usually preferred for the numerical treatment of such laws since they directly provide a physical sense to the quantities of interest:

• Lagrange multipliers method [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF];

• Penalty method [START_REF] Mcguire | Matrix Structural Analysis:Second Edition[END_REF];

• Augmented Lagrangian method [START_REF] Harmon | Robust, Efficient, and Accurate Contact Algorithms[END_REF].

The main difference between the first two approaches consist in the fact that Lagrange multipliers method enforces exactly the non penetration conditions, while Penalty method allows a loosely controlled violation of the constraints between the structures, as shown in Fig. 2.8. Augmented Lagrangian method is instead a combination of these two approaches.

Figure 2.8: Graphical comparison of the contact constraints fulfillment by Lagrange multipliers and Penalty method. Lagrange method does not allow any penetration between the two bodies, while Penalty does not fulfill perfectly the contact constraint allowing a certain amount of penetration [START_REF] Harmon | Robust, Efficient, and Accurate Contact Algorithms[END_REF].

Following a contact event, the equations of motion of the system are modified in the following way:

Mq(t) + C q(t) + Kq(t) = p ext (t) + p c (t) (2.4.1)
where p c (t) is an added term and corresponds to the vector of contact loads. Before proceeding with the contact treatment, it is opportune first to distinguish between the master body and the slave body. Physically they cannot penetrate each other but numerically the master body can penetrate the slave body while the contrary is not allowed. Normally the rigid body is set to be the master.

Once made this distinction, the procedure consists in finding at each time iteration, for any material point x belonging to the master surface (Γ m ), its closest node y belonging to the slave surface (Γ s ). According to this condition, the gap function between the two bodies surface is:

g(x, t) = g 0 (x, t) + (q m (x, t) -q s (ȳ(x, t)) • n (2.4.2)

Normal contact condition

The normal contact conditions, widely discussed in the literature, very often are expressed according to Karush-Kuhn-Tucker conditions [START_REF] Gitterle | A dual mortar formulation for nite deformation frictional contact problems including wear and thermal coupling[END_REF]:

     g(x, t) ≥ 0 p c,N ≤ 0 p c,N g(x, t) = 0 (2.4.3)
The first condition of the system expresses the non-penetration constraint which does not allow the bodies to penetrate each other, the second condition implies that no adhesive stresses occur in the contact zone and the third condition ensures that the gap is closed when non-zero pressure occurs and the pressure to be zero when the gap is open. A graphic representing these conditions is given in Fig. 2.9.

Figure 2.9: Are shown the contact conditions in normal direction. When the gap function is greater than zero (no contact), the normal contact pressure is zero. While when the gap is equal to zero (contact) a non-zero normal contact pressure occurs. No negative gap function is present because no bodies penetration is allowed [START_REF] Gitterle | A dual mortar formulation for nite deformation frictional contact problems including wear and thermal coupling[END_REF].

Tangential contact condition

Models to describe frictional effects are several because it is a complicated phenomenon. In continuum mechanics, the most used model is the one introduced by Coulomb, called frictional Coulomb model. A graphical representation of such method is given in Fig. 2.10.

Figure 2.10: Frictional contact conditions in tangential direction according to Coulomb's model [START_REF] Gitterle | A dual mortar formulation for nite deformation frictional contact problems including wear and thermal coupling[END_REF].

Equivalently, Coulomb's model can be set as follows [START_REF] Gitterle | A dual mortar formulation for nite deformation frictional contact problems including wear and thermal coupling[END_REF]:

     p c,T -µ|p c,N | ≤ 0 p c,T -µ|p c,N | < 0 ⇒ v T,rel (x, t) = 0 p c,T -µ|p c,N | < 0 ⇒ ∃β ∈ R + 0 : v T,rel (x, t) + βp c,T = 0 (2.4.4)
where v T,rel is the relative velocity between the contact surfaces and µ is the frictional coefficient.

Lagrange multiplier method

The principal characterization of Lagrange multiplier method is a perfect compliance with the contact constraint limit, allowing no penetration of the slave into the master body. This contact constraint is introduced in the system adding a new unknown, the so-called Lagrange multipliers λ.

The new restrained equation of motion is [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF]:

Mq i+1 + C qi+1 + Kq i+1 + G T i+1 λ i+1 = p ext i+1 (2.4.5)
where G is the restrained displacements matrix on the surface contact which expresses the exact non penetration condition:

G i+1 = C T N q i+1,c + g i+1,p ≥ 0 (2.4.6)
where g is the gap function, C N the contact displacement constraint matrix in normal direction and q the nodal displacements vector in which the pedix c stands for "corrected". Setting to zero in the gap vector g i+1,p all the positive terms and keeping unchanged the negative ones, this leads to a displacement correction only for the penetrating nodes. The new expression for G i+1 becomes:

G i+1 = C T N q i+1,c + g - i+1,p = 0 (2.4.7)
Thus, the new system to solve is:

M 0 0 0 qi+1 λi+1 + C 0 0 0 qi+1 λi+1 + K G T G 0 q i+1 λ i+1 = p ext i+1 0 (2.4.8)
This is the system of equations corresponding to Lagrange multipliers approach for contact treatment, which can be combined with any time integration scheme in order to get the solution in terms of λ i+1 , which corresponds to the contact forces vector, and q i+1 , the vector of nodal displacements taken into account the correction term due to the contact constraints. Among the time integration schemes, the one more suitable for contact problems is the Central difference scheme [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF], which has been proven to be the one that allows a better tracking of the solution in time avoiding any spurious oscillations.

The advantage of the Lagrange multipliers method is the exact solution that gives since the nonpenetration condition is fully fulfilled. While the disadvantages are associated to the addition of new unknowns to the system (λ), the loose of the positive definiteness of the matrices , zero mass is associated to Lagrange multipliers which could lead to infinite eigenfrequencies and finally a higher complexity in the implementation with respect to Penalty method.

Penalty method

This approach is called Penalty method because of the penalization of the contact constraints by a parameter called Penalty coefficient. The normal contact condition for such method is shown in Fig. 2.11.

Figure 2.11: Penalty regularization of normal contact conditions allowing a certain amount of penetration in fucntion of the normal contact pressure. The intensity of this penetration is expressed by the Penalty coefficient (C N ), higher is this value and lower is the allowed penetration and vice versa [START_REF] Gitterle | A dual mortar formulation for nite deformation frictional contact problems including wear and thermal coupling[END_REF].

Thus, form Fig. 2.11 the following conditions can be stated:

p c,N = 0 for g > 0 (no contact) C N g for g ≤ 0 (contact) (2.4.9)
where p c,N is the contact pressure in normal direction, g the gap value and C N the penalty parameter. For a general case the contact force along the tangential direction, recalling the Coulomb frictional contact model, can be formulated as [START_REF] Gitterle | A dual mortar formulation for finite deformation frictional contact problems including wear and thermal coupling[END_REF]:

£ v p c,T = C T v T,rel -β p c,T p c,T (2.4.10)
where £ v p c,T is the LIE derivative of the frictional force, v T,rel the relative tangential velocity between the two body surfaces and C T the tangential penalty parameter. The approach suggested by Penalty method causes a not exact fulfillment of contact constraint, but leads to a direct connection between the kinematic quantities and the contact forces leading to a simpler formulation in which the introduction of additional unknowns, contrarily to Lagrange multipliers method, is not necessary and the problem can be formulated in displacement terms only when, at a later stage, the weak form is derived [START_REF] Gitterle | A dual mortar formulation for finite deformation frictional contact problems including wear and thermal coupling[END_REF].

The main disadvantage of this method is the violation of contact constraints. This drawback can be overcomed by increasing the penalty parameter value, thus enforcing better the contact constraint, but this leads to an ill conditioning of the system equations. Therefore, the penalty parameter has to be a compromise between these two limitations.

Augmented Lagrangian method

A third method to regularize the normal contact and friction conditions is the Augmented Lagrangian method. This third approach has been introduced to combine the advantages given by the two approaches of Lagrange multipliers and Penalty methods, overcoming their drawbacks. It consists in evaluating the normal contact traction in terms of both the Lagrange multiplier λ N in normal direction and the penalty parameter C N :

p c,N = 0 for λ N + C N g > 0 (no contact) λ N + C N g for λ N + C N g ≤ 0 (contact) (2.4.11)
while the tangential contact force is expressed as [START_REF] Gitterle | A dual mortar formulation for finite deformation frictional contact problems including wear and thermal coupling[END_REF]:

£ v (p c,T -λ T ) = C T v T,rel -β p c,T p c,T (2.4.12)
where λ T is the tangential Lagrange multiplier vector, v T,rel the relative tangential velocity between the two body surfaces and C T the tangential penalty parameter. By this formulation, the enforcement of the contact conditions is higher with respect to Penalty method because of the introduction of the Lagrange multiplier term, and at the same time the desired connection between the kinematic quantities and the contact forces is introduced, simplifying in such way the formulation with respect to Lagrange multipliers method.

Wear

In many engineering applications the contact and sliding between different parts is very common. In such conditions, the process of gradual removal of a material from the surfaces is initialized and it is what commonly is called wear.

The removal of material from the bodies surfaces as consequence of wear phenomena results in an increase of the gap between the contacting bodies, leading to changements in the application performances. Indeed almost all machines durability and reliability is reduced due to wear limiting in such way the possibility of developing new advanced machines [START_REF] Kato | Wear Mechanisms[END_REF]. Therefore, it is important to predict such wear depth evolution, as well as its pattern.

Kind of wear

In order to proceed with the wear control, it is first needed to classify the kind of wear that interests the analyzed application, in order to develop an adapt model to predict its wear evolution.

A first classification of wear patterns can be the following [START_REF] Zmitrowicz | Wear patterns and laws of wear, a review[END_REF]:

• abrasion;

• adhesion;

• fatigue;

• corrosion;

• erosion.

Abrasive wear

In the abrasion patterns, the material is removed from a soft surface after that a hard particle digs into it and plows out material. This ploughing takes place in sliding. An illustration of such phenomenological contact is given in Fig. 2.12.

Figure 2.12: Abrasive wear of ductile material [START_REF] Kato | Wear Mechanisms[END_REF].

As it is possible to see in the figure in this kind of wear, even if initially the contact occurs between flat surfaces parallel to the sliding direction, the interface rotates and becomes inclined and wavy [START_REF] Kato | Wear Mechanisms[END_REF]. This effect is a result of the combination of normal and tangential forces in sliding.

Adhesive wear

Adhesive wear is a result of the plastic deformation due to the local adhesion of the contact surfaces. In such mechanism the load applied to the contacting asperities that characterizes the surfaces is so high leading to their deformation and adhesion to each other forming micro-joints [START_REF] Kato | Wear Mechanisms[END_REF]. This material adhesion is shown in Fig. 2.13. Figure 2.13: Adhesion wear mechanism as result of rubbing of two surfaces [START_REF]Substances and Technologies, Knowledge source on Materials Engineering[END_REF].

To decrease adhesive wear between contact surfaces different solutions are adopted according to the studied application. For instance, one of these solutions is lowering the load, use of harder rubbing materials, use of contaminated rubbing surfaces, introduction of solid lubricants, use of a lubrication oil between the contact surfaces or introduction of anti-wear additives in oil [START_REF]Substances and Technologies, Knowledge source on Materials Engineering[END_REF].

Fatigue wear

In fatigue patterns, the surface is not damaged slowly and continuously but suddenly due to large failure of interacting surfaces [START_REF] Zmitrowicz | Wear patterns and laws of wear, a review[END_REF]. This sudden failure may occur as result of repeated loading and unloading cycles. Due to these cycles, it can be initialized the appearance and growth of surface and subsurface cracks leading to the formation of surface defects and eventually the sudden removal of material from the solids surface, as shown in Fig. 2.14.

Figure 2.14: Fatigue wear mechanism as result of rubbing of two surfaces [START_REF]Substances and Technologies, Knowledge source on Materials Engineering[END_REF].

When the number of such loading cycles is high, the expected wear phenomenon is the highcycle fatigue mechanism, while when it is low, the low-cycle fatigue mechanism is expected [START_REF] Kato | Wear Mechanisms[END_REF].

Corrosive wear

If the rubbing surfaces show a stadium of oxidation, then this will have effect on the wear mechanism leading to its acceleration [START_REF] Kato | Wear Mechanisms[END_REF]. The oxidation initialization is highly promoted by the friction phenomenon of the contact surfaces. The effects of the friction are indeed an increase of the temperature, as well as the removal of the protecting oxide films. Furthermore, the removed oxide particles from the surface can be trapped between the contact surfaces leading to the initialization of wear by abrasion mechanism of three bodies in contact [START_REF]Substances and Technologies, Knowledge source on Materials Engineering[END_REF].

Erosive wear

The Erosive wear mechanism is characterized by the impact of particles that could be solid, liquid or gaseous, against the surface of a body [START_REF] Hirani | Wear Mechanisms-2[END_REF]. These impacting particles remove material gradually from the impacted surface through repeated deformations and cutting actions. In engine bearings, the erosion is caused principally by cavitation in the lubrication oil. When the lubricant oil exits from the convergent gap between the bearing and journal surfaces its pressure decreases rapidly leading to the formation of voids. These formed bubbles can collapse producing a shock wave which removes particles from the bearing [START_REF]Substances and Technologies, Knowledge source on Materials Engineering[END_REF].

Wear modeling

When rubbing of surfaces is involved in the studied application, first it is important to define what is the kind of wear mechanism that occurs among the wear mechanisms discussed in section 2.5.1 in order to proceed with its modeling. For the wear modeling a large body of empirical data has been collected and some phenomenological models have been developed [START_REF] Zmitrowicz | Wear patterns and laws of wear, a review[END_REF]. Nevertheless, there is still a need of reliable and more robust procedures which can take into account complex phenomena of wear. One of the most used models to compute the wear in a structure is Archard's law due to its simplicity and at the same time applicability to several different cases.

Archard's law

The Archard's law is an extension of Holm's law, which was developed for electric contact cases [START_REF] Zmitrowicz | Wear patterns and laws of wear, a review[END_REF]. Holm formulation states a relationship between the depth of the material removed by wear (d w ) on both the sliding distance (s) and the true area of contact (A). Archard extended such relation observing that the removed material increases as the exchanged normal load (p n ) between the contact surfaces increases, and decreases as the hardness of the material increases (H) [START_REF] Luoa | Numerical analysis of die wear characteristics in hot forging of titanium alloy turbine blade[END_REF].

d w (T ) = K w (T ) p n s H(T ) (2.5.1)
K w is the Archard's wear coefficient and it is a function of the temperature, as well as it is the hardness. K w is an empirical coefficient and its expression varies from case to case.

In case in which the friction is taken into account, the relation becomes:

d w (T ) = µK w (T ) p n s H(T ) (2.5.2)
The coefficient of friction (µ) is used in this formulation of Archard's law in order to emphasize that the wear rate is related to the contact stress through the frictional work done. In a dynamical simulation, in case in which the relative sliding velocity between the two bodies, the contact pressure and temperature fields vary with the position and time, eq.(2.5.2) is re-expressed as:

∆d wij (T ) = µK w (T ) p nij v T ij H(T ) dt j (2.5.3)
in which v T ij is the relative sliding velocity, dt j the time step length and i and j repsectively the discretizations in space and in time.

At the end of the process, the wear depth of the i th node is given by:

d wi (T ) = n j=1 ∆d wij (T ) (2.5.4)
where n is the total number of steps of the simulation.

Adaptation of the FEM model to wear phenomena

In order to predict the pressure distribution at each time step in a simulation that involves wear phenomena on one or both the interacting bodies, a continuous adaptive meshing has to be performed to eliminate or at least reduce the deformation-induced element distortion. The methodology followed consists normally in computing first the pressure distribution by FEM in order to be able to apply Archard's law which allows to calculate the nodal wear increments. The next time step pressure distribution is obtained on the modified geometry and the cycle is repeated according to the simulated wear process. This methodology characterized by moving the contact nodes according to the nodal increment wear values is called re-meshing technique [START_REF] Békési | Wear simulation of a reciprocating seal by global remeshing[END_REF].

Another diffused technique consists in eliminating the worn out elements instead of moving their nodes and it is called killing elements technique.

In re-meshing method at each iteration the allowed material quantity to be removed from the interacting body is higher and thus requires a smaller number of contact iterations to simulate a certain sliding distance. While when re-meshing is not used this allowable amount of removable material is lower since the surface elements are continuously worn out and so the allowable sliding distance for each wear step is comparatively lower with respect to the remeshing technique case and is continuously decreasing, requiring more contact simulations for the same predetermined maximum sliding distance [START_REF] Hegadekatte | Finite element based simulation of dry sliding wear[END_REF].

To overcome this limitation in the elements elimination techniques it is possible to increase the mesh dimension to have enough height to accommodate more amount of wear. But such a strategy for meshing will affect the accuracy of the FE results because it is not suggested to have a coarse mesh in the contact region in order to obtain a good accuracy of the results. Thus, an important advantage with re-meshing techniques is the lower computational time needed for the wear simulation and at the same time a higher accuracy of the obtained results.

In the following two sections will be presented some developed remeshing and killing elements techniques fully used nowadays in different engineering applications.

Re-meshing

The developed re-meshing techniques along the years are several and differs according to the application field.

They can involve the re-meshing of the whole body part, global re-meshing, or more efficiently only the local area directly interested by the wear phenomena, local re-meshing.

Local re-meshing techniques differ in their turn in Single-layer re-meshing and Multi-layer remeshing.

Single-layer re-meshing techniques

Single-layer re-meshing techniques consist in accommodating the wear on a single FE layer in proximity of the contact surface at each iteration. Thus, they require the use of eventually a dynamic time step instead of a constant time step, based on the maximum amount of wear allowable for each iteration [START_REF] Podra | Simulating sliding wear with finite element method[END_REF]:

dt new j = dt j dh max dh j (2.6.1)
At each j th iteration, the time step is defined by the maximum allowable nodal wear increment dh max . If the computed wear increment for the j th iteration is higher than the maximum allowed wear increment dh max , then the time step is decreased and it is recomputed according to eq.(2.6.1) to respect such limitation. This maximum threshold wear depth is chosen based on the simulation conditions. It is important to note that a low dh max leads to a high results accuracy but also high computational time, while a high dh max leads to a low results accuracy but also low computational time. Therefore, it is very important to choose a proper dh max , according to the analyzed application.

Once the dynamical time step has been computed, the single-layer re-meshing followed is illustrated in Fig. 2.15.

Master Thesis -25-Khadija Masrour ). This movement of the contact node 1 is based on the Archard wear model. Once H i 13 reaches this critical height, the node 3 is moved towards node 5 (Fig. 2.15(b-3)). This movement is also incremental and it is about dh max at each increment. Once also element 2 reaches the critical height H cr 35 , the node 5 is moved towards node 7 and so on. The principal drawback of this technique is a fully dependence of i th step on the previous remeshing history. Therefore for the described procedure it has been thought important as a backup option the storage of a selected output data set into a special file after every solution step to save the previous steps data if, for any reason, the analysis is interrupted [START_REF] Rezaei | Adaptive finite element simulation of wear evolution in radial sliding bearings[END_REF].

Multi-layer re-meshing techniques

In such techniques the allowed incremental wear for each iteration is much higher than the Singlelayer re-meshing techniques, thus could be higher than the elements height and is accomodated simultaneously on multiple nodes layers. The registered wear is then partitioned among different elements or equally, or by a linear relation with nodes displacement that decreases as the nodes distance increases from the contact surface.

These two options of partitioning the wear equally or linearly among the above elements are described here below.

Equally displaced nodes: The described technique is based on a simple principle according to which each element can substain only a limited amount of wear that is defined as a fraction of its initial length, called critical height. Once this limit is reached, the wear is distributed to the above elements in an equal amount [START_REF] Benabdallah | Finite element simulation of the wear of polyoxymethylene in pin-on-disc configuration[END_REF]. Since the total amount of wear during a simulation goes increasing, it is possible to establish the maximum number of elements required to accommodate the reached total wear (d tot w ) at each time iteration by the following relation:

d tot w H cr H 0 = q + r (2.6.2)
where d tot w is the total amount of wear reached from the beginning of the simulation, H cr is the preset maximum critical height of the element and H 0 is its initial height. While q is the whole number of elements required to accommodate d tot w and r is the fraction of the remaining elements. Both q and r are unique to each contact node. An illustration of the described procedure is represented in Fig. 2.16, in which is shown the body configuration at two consequent iterations, j th and (j + 1) th . In the (j + 1) th configuration q element nodes are equally shifted above to accommodate the total amount of wear h tot w,j+1 and the first element is shifted taking into account also the remaining r fraction.

Figure 2.16: Body configuration at the j th and (j + 1) th configuration, showing its nodal displacements to accommodate the total amount of wear h tot,j+1 according to the equally displaced nodes technique [START_REF] Benabdallah | Finite element simulation of the wear of polyoxymethylene in pin-on-disc configuration[END_REF]. a is the critical element height H cr The principle advantage of such technique is that the solution at each iteration is not dependent on the body history, the cumulative wear data is the only information carried over from cycle to cycle. The drawback instead has been prooved to be an unstable procedure that not always lead to the convergence of the solution in case of a 2D or 3D body. BENABDALLAH and OLENEDER shew in [START_REF] Benabdallah | Finite element simulation of the wear of polyoxymethylene in pin-on-disc configuration[END_REF] that unrealistic discontinuities appeared on the wear profile after just a few sliding cycles (Fig. 2.17).

Figure 2.17: Mesh distortions that appeared after few sliding cycles of a pin on a disc in which the wear has been accommodated following the equally displaced nodes technique [START_REF] Benabdallah | Finite element simulation of the wear of polyoxymethylene in pin-on-disc configuration[END_REF].

Linearly displaced nodes: In this technique, since it is a Multi-layer re-mashing technique, the remesh is extended simultaneously to multiple rows of nodes. The updated nodes of the whole simulation are only the ones inside a pre-set region that extends until the height H. The depth H should be larger than the estimated value of maximum wear at the end of the simulation, but at the same time not too large to not increase much the computational time [START_REF] Tang | A multilayer nodes update method in FEM simulation of large depth fretting wear[END_REF]. Once selected the region of the nodes to update, the node displacement decreases linearly as its distance increases with respect to the contact surface as illustrated in Fig. 2.18. Following Fig. 2.18 and focusing on node Q which share the same x coordinate with the corresponding surface node C, according to the described technique, it will move of:

d Q = H -y Q H -y C d C (2.6.3)
in which H is the maximum selected depth to accommodate the whole wear, y Q and y C are respectively Q and C y coordinates from the contact surface, d Q and d C are their relative movement.

In a more general case in which the inner node is not equally aligned with a surface node with respect to y direction, so it does not share with it x coordinate, its wear displacement can be calculated as:

d P = H -y P H -y N d N (2.6.4)
in which the wear depth of the projected node N is approximated by a linear interpolation of the wear depths of nodes B and C:

d N = x C -x N x C -x B d B + x N -x B x C -x B d C (2.6.5)
The corresponding y coordinate can be obtained similarly as:

y N = x C -x N x C -x B y B + x N -x B x C -x B y C (2.6.6)

Elimination of elements

Another popular approach for modelling wear phenomena is the elimination of elements technique. This material deattachment is taken into account deactivating or deleting elements from the FE mesh.

The elimination process of such elements from the FE model is based on the characteristic damage criterion of the given material [START_REF] Békési | Modelling Friction and Abrasive Wear of Elastomers[END_REF]. This approach is more commonly used in manufacturing processes.But it can be used also in wear applications since, from the simulation point of view, it is similar to the manufacturing processes.

In wear problems the killing elements technique is commonly integrated with re-meshing techniques. Indeed at each iteration the nodes are displaced according to re-meshing process in order to accommodate the wear and if the damage criterion for the considered elements is fulfilled, then they are eliminated and the geometry is updated according to the new configuration before moving to the next iteration. Such methodology is illustrated in the flowchart represented of Fig. 2.19.

Elements for plate FEM modeling

In a 2D domain different kind of elements are available to mesh the structure. The most diffused ones are triangular and quadrilateral elements, thus involving respectively 3 nodes per element and 4 nodes per element. However the advantages of the second type of elements are fairly higher. Indeed these ones are associated to varying and not constant strain matrix, thus they give a distribution more realistic of strain and so of stress on the structure [START_REF]FEM for Two-Dimensional Solids (Finite Element Method) Part 2[END_REF]. Therefore an illustration and discussion on this kind of element will be performed. Dividing the whole structure into perfect rectangular elements is impossible, thus the exploitation of generic quadrilateral elements with straight but unparallel edges is almost always required and it is widely available and described in the literature.

In order to conserve the simplicity of quadrilateral elements with parallel edges, a mapping of the physical coordinates into new coordinates called Natural coordinates is performed as illustrated in Fig. 2.20. In FE method the displacement within an element is expressed through shape functions that interpolate its nodal displacements. Taking into account such mapping, the internal element displacement in fuction of the natural coordinates ξ and η is expressed as:

U(ξ, η) = N(ξ, η)d e (2.7.1)
Figure 2.19: Flowchart of elimination elements approach combined with re-meshing technique for wear simulation [START_REF] Békési | Modelling Friction and Abrasive Wear of Elastomers[END_REF]. At each iteration the contact calculations are performed allowing to evaluate the stress distribution. If the damage criterion for the analyzed element is fulfilled then the element is deactivated/killed otherwise a re-meshing is performed moving the element nodes.

where N is the matrix of shape functions, U the vector of element displacement in function of the natural coordinates and d e the vector of nodal displacements.

Similarly, the vector of internal element coordinates can be expressed as:

X(ξ, η) = N(ξ, η)x e (2.7.2)
where x e is the vector of nodal coordinates.

Denoting the natural coordinates of each quadrilateral node as (ξ j , η j ), the shape function have the following expression:

N j = 1 4 (1 + ξ j ξ)(1 + η j η) j = 1, 2, 3, 4 (2.7.3)
Hence, in case of quadrilateral elements in which each node has 2 DoFs the matrix of shape functions is the following:

N = N 1 0 N 2 0 N 3 0 N 4 0 0 N 1 0 N 2 0 N 3 0 N 4 (2.7.4)
For the FEM modeling the stiffness and mass element matrices are expressed as [START_REF] Griffiths | Stiffness matrix of the four-nodes quadrilateral element in closed form[END_REF]:

   K e = V B T DBdV M e = V ρN T NdV (2.7.5)
where D is the strain/stress matrix:

N =   E 1 E 2 0 E 2 E 1 0 0 0 G   (2.7.6)
with:

E 1 = E(1 -ν) (1 + ν)(1 -2ν) , E 2 = νE (1 -ν) , G = E 2(1 + ν) (2.7.7)
While matrix B is given by [START_REF] Nikishkov | Introduciton to the Finite ELement Method, Isoparametric transformations[END_REF]:

B = B 1 B 2 B 3 B 4 with : B i =      ∂N i ∂x 0 0 ∂N i ∂y ∂N i ∂y ∂N i ∂x      (2.7.8)
The partial derivative of the shape functions with respect to the physical coordinates x and y can be computed as follows:

   ∂N i ∂x = ∂N i ∂ξ ∂ξ ∂x + ∂N i ∂η ∂η ∂x ∂N i ∂y = ∂N i ∂ξ ∂ξ ∂y + ∂N i ∂η ∂η ∂y
(2.7.9)

The equivalent matricial form is:

   ∂N i ∂x ∂N i ∂y    =   ∂ξ ∂x ∂η ∂x ∂ξ ∂y ∂η ∂y      ∂N i ∂ξ ∂N i ∂η    = ∂(ξ, η) ∂(x, y)    ∂N i ∂ξ ∂N i ∂η    = J -1    ∂N i ∂ξ ∂N i ∂η    (2.7.10)
where J is the Jacobian matrix. The determinant of the Jacobian matrix is used to express the integral in dV from the physical coordinates to the natural coordinates:

dV = dxdy = |J|dξdη (2.7.11)
Thus the new expression of the stiffness ans mass element matrices is:

   K e = V B T DB|J|dξdη M e = V ρN T N|J|dξdη (2.7.12)
The integration is performed then following Gauss quadrature integration rule:

       K e = p i=1 n j=1 w ij B T DB|J|dξ i dη j M e = p i=1 n j=1 w ij ρN T N|J|dξ i dη j (2.7.13)
where p is the number of the integrating points along ξ and along η and i and j their indices. w ij are the relative weightning coefficients.

A graphical representation of the Gauss integration rules is given in Fig. 2.21, while in Tab.2.1 are shown the coordinates and the relative weights arriving until 3 × 3 integration points.

Chapter introduction

As first step in order to model the wear evolution on the blade tip it is required first to simulate the contact interactions with the surrounding casing. As already discussed, since the aim of this research project is the validation of an alternative technique to simulate the wear evolution on interacting 1D bodies, the system used corresponds to a one side clamped beam that goes in contact interactions with a wall that moves along both the directions, horizontal and vertical. The contact algorithms developed are the Lagrange multipliers method and the Penalty method. A focus will be addressed to discuss the differences between these two methods as well as their advantages and drawbacks. A model reduction according to Craig-Bampton method will be finally applied on both the algorithms with the aim of reducing the computational costs.

Beam model

As already mentioned, the research work is performed on the analysis of the contact interactions of a beam with a moving wall taking into account the wear phenomena. In this section will be presented the beam model used. The used beam is one side clamped while the other free extremity will be the one that will be subjected to contact inetractions with the wall. Its geometrical and material properties are resumed in Tab.3.1. The material used is a titanium alloy (Ti-6AI-4V) widely used for compressor blades in aerospace field. The beam has been modeled by the finite element method considering three degrees of freedom 3.2: The first six natural frequencies of the beam.

Properties

K el =         EA l 0 0 -EA l 0 0 0 12EIz l 3 6EIz l 2 0 -12EIz l 3 6EIz l 2 0 6EIz l 2 4EIz l 0 -6EIz l 2 2EIz l -EA l 0 0 EA l 0 0 0 -12EIz l 3 -6EIz l 2 0 12EIz l 3 -6EIz l 2 0 6EIz l 2 2EIz l 0 -6EIz l 2 4EIz l         (3.2.1) M el = ρAl 420         140 0 0 70 0 0 0 156 22l 0 54 -13l 0 22l 4l 2 0 13l -3l 2 70 0 0 140 0 0 0 54 13l 0 156 -22l 0 -13l -3l 2 0 -22l 4l 2         (3.2.2)

Damping model

Once modeled the conservative system by the mass and stiffness matrices, it is opportune now to introduce also a physical damping to model its non conservative behaviour. There exist different damping models, the one selected for this analysis is based on proportional assumption in order to get modal damping ratios of 0.2% for the two first modes. According to the proportional damping assumption the damping matrix can be written as:

C = aK + bM (3.3.1)
This gives:

β i = aγ i + bµ i (3.3.2)
where γ i is the modal stiffness and µ i the modal mass of the i th mode. So:

i = 1 2 (aω i + b ω i ) (3.3.3)
where i is the modal damping ratio.

From the following system knowing the value of the modal damping ratio for the first two modes it is possible to get a and b coefficients:

   1 = 1 2 (aω 1 + b ω 1 ) 2 = 1 2 (aω 2 + b ω 2 ) (3.3.4)
where:

• ω 1 = 5303.27 rad/s

• ω 2 = 33245.88 rad/s

• 1 = 0.002

• 2 = 0.002
Once computed a and b it is possible to build the model damping matrix by the relation 3.3.1.

Test case conditions

For the performed contact interaction analysis the clamped beam is considered fixed in space while the only moving body is the wall. The wall used is rigid and it is moving simultaneously in vertical direction with a constant upward velocity and horizontally following a harmonic law (Fig. 3.3). In a first performed simulation the wall penetrates the beam until a depth of 40µm in order to reproduce the contact conditions illustrated in [START_REF] Berthoul | Phenomenological modeling of abradable wear in turbomachines[END_REF] in which it has been estimated that normally the maximum incursion ranges from zero to about 50µm. Thanks to this harmonic motion the beam is left free after each interaction giving it some time to vibrate freely before the next contact interaction.

In a second simulation in order to simulate much more severe interacting conditions a penetration until 300µm will be simulated.

When the wall starts approaching the beam tip the positive gap function starts decreasing. It becomes zero when the wall is in perfect contact with the tip and negative when it starts penetrating the beam. In these conditions the contact algorithms correct such gap function setting it to zero (Lagrange multipliers method) or almost zero (Penalty method). The response of the beam to such penetrations along the axial direction and the upward motion of the wall are displacements in both the axial direction and the transversal direction which intensities depend on the contact interaction conditions. Such responses will be presented in the following sections but first it will be presented the time integration technique exploited to compute these solutions.

Time integration

To get the solution of the analyzed case in terms of displacements, velocities and accelerations in time it is opportune to integrate the equations of motion of the system at each time step. The damped system equations of motion are given by:

Mq(t) + C q(t) + Kq(t) = p(t) (3.5.1)
in which q(t) is the vector of nodal coordinates displacements, M the mass matrix get by the FEM discretization, C the damping matrix and K the stiffness matrix. On the right side of the equation appears the external dynamic loads p(t) vector, which in the studied case includes both the normal and tangential contact forces of the tip beam with the rigid wall. Such forces can be computed by Lagrange multiplier or Penalty contact methods, which will be fully described in section3.6.

Available time integration schemes are several and their selection depends on the studied application.

In different studies it has been prooven that using Implicit Newmark schemes to solve contact problems yields spurious oscillations and sizeable deviations from the exact solution [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. Therefore, explicit schemes are rather preferible than implicit schemes in solving contact problems, in particular the central time scheme. This scheme can be recovered by Newmark using γ = 1/2 and β = 0 or implemented directly. In this case the advantage is in the much simpler involved equations as it will be presented here below, which decreases significantly the computational effort.

In Central difference scheme the first and second derivatives of q(t) are approximated as:

q(t i ) = qi ≈ q i+1 -q i-1 2h (3.5.2) q(t i ) = qi ≈ q i+1 -2q i + q i-1 h 2 (3.5.3)
Substituting eq.(3.5.2) and eq.(3.5.3) into eq.(3.5.1), this leads to:

M h 2 + C 2h q i+1 ≈ 2M h 2 -K q i + C 2h - M h 2 q i-1 + p i (3.5.4)
which can be re-written in a more compact form as:

A 0 q i+1 ≈ A 1 q i + A 2 q i-1 + p i (3.5.5)
Since A 0 is a positive definite matrix, it may be factorized prior to the iterative solution by, for instance, a LDL T decomposition.

As it is possible to notice by this algorithm the iterative equation to solve at each time step is very simple. Although, since Central difference scheme is an explicit method, it is characterized by a conditional stability. Thus, the stability condition to respect for the time step selection is:

h max ≤ 1 πf max (3.5.6)
in which h max is the maximum allowed time step length and f max the maximum frequency of the discretized system.

Contact algorithms

In this section will be presented and described from a theoretical point of view the two contact algorithms used to solve the problem, the Lagrange Multipliers method and the Penalty method. Furthermore, a subsection will be dedicated to the presentation and validation of the results obtained by these two algorithms in which the contact conditions that occurs are the one described in section 3.4 with a maximum penetration of the wall into the beam of 40µm.

A last section finally will be dedicated to the comparison of the results recovered by these two methods.

Lagrange Multipliers method

In Lagrange Multipliers method the contact constraint between the two bodies is exactly enforced at each iteration, thus no penetration is allowed. According to such condition, in order to calculate the solution, the algorithm can be divided essentially into three fundamental steps [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF]:

(1) Prediction of the nodal displacements vector. Indeed at each time step the predicted solution of the equivalent non-contact case has to be computed applying the Central scheme difference expression:

q i+1,p ≈ M h 2 + C 2h -1 2M h 2 -K q i + C 2h - M h 2 q i-1 + p i (3.6.1)
(2) Determination of the gap function g i+1,p between the two surfaces. Applying the constraint contact condition leads to:

G i+1 = C T N q i+1,c + g i+1,p ≥ 0 (3.6.2)
where C N is the contact displacement constraint matrix in normal direction and q is the nodal vector displacements in which the subscript c stands for "corrected". The same constraint condition is fulfilled if all positive terms of g i+1,p are set to zero while all negative terms are kept unchanged, this leads to a displacement correction only for the penetrating nodes. The new constraint condition can be re-expressed as:

G i+1 = C T N q i+1,c + g - i+1,p = 0 (3.6.
3)

The new system equations according to Lagrange multipliers method taking into account the contact constraint becomes:

M 0 0 0 qi+1 λi+1 + C 0 0 0 qi+1 λi+1 + K G T G 0 q i+1 λ i+1 = p ext i+1 0 (3.6.4)
(3) Correction of the nodal displacements taking into account the contact constraint. Combining eq.(3.6.4) and the Central difference scheme velocity and acceleration approximations of eq.(3.5.2) and (3.5.3), the following solution is obtained:

λ i+1 = C T N M h 2 + C 2h -1 C N T -1 g i+1,p (3.6.5)
where C N T is the contact displacement constraint matrix that takes into account also the frictional effects. By this expression it is possible to compute the Lagrange multipliers vector which represent the contact forces.

Finally the vector of the nodal displacements taking into account the correction term is:

q i+1 = q i+1,p + q i+1,c = q i+1,p + M h 2 + C 2h -1 C N T λ i+1 (3.6.6)

Test case results

In the test case, according to the contact conditions described in section 3.4, the wall is moving vertically with a constant velocity of 176m/s and horizontally with a sinusoidal law of 30Hz. The maximum incursion of the wall is of 40µm and the frictional coefficient is of 0. In Fig. 3.4 it is shown the axial vibration of the external tip of the beam as response to these contact conditions. It is possible to notice how the beam tip under contact with the wall penetrates until a value of 40µm and once the beam is free from contact it registers axial free vibrations that are slowly damped out because of the low physical proportional damping introduced in the system.

It is important to notice also that the frequency of these free vibrations is approximately 12500Hz, which corresponds to the third frequency of the beam. The simulation duration is of 25ms and the wall penetration into the beam occurs from ∼ = 1ms to ∼ = 15.1ms. Since such penetration occurs relatively slowly, different vertical pushings of the beam have the time to occur with a frequency of approximately 870Hz. In Fig. 3.5 it is represented the vertical response of the beam tip to the above mentioned contact conditions. During the wall penetration phases the beam tip is pushed upward until a height of ∼ = 5mm, while once it is no more in contact with the wall, it registers free vibrations of approximately 840Hz, which corresponds to its first frequency.

Study of frictional coefficient

In this study the same test conditions described in section3.4 have been reproduced. The only varying parameter is the frictional coefficient, which it is expected to influence the vertical vibrations of the blade since the frictional forces depend directly on it. In Fig. 3.6 are shown different vertical displacements of the beam obtained for three different values of the frictional coefficient, 0.25, 0.3 and 0.35. It is possible to observe that higher is the frictional coefficient and higher is the vertical displacement since the beam is pushed away with higher energy.

In the performed study the reproduced free vibrations of the beam for the different frictional coefficients are very similar, although it is possible to have cases in which the free vibrations vary significantly and this could be explained by the conditions in which occurs the last contact with the wall. If the last contact occurs in a point in which the blade presents a high deviation from its neutral point then it is expected that it vibrates with higher energy and so higher amplitudes, and vice versa. Higher is the frictional coefficient and higher is the frictional force, thus, higher is the energy with which the beam is pushed upward leading to higher vertical displacements. The results are obtained applying Lagrange multipliers method with contact conditions described in section3.4.

Penalty method

Contrarily to Lagrange multipliers method the contact constraints in Penalty method are not exactly fulfilled but a certain amount of penetration between the bodies is allowed depending on the selected value of the Penalty coefficient. The contact conditions (shown graphically in Fig. 2.11) are the following:

p c,N = 0 for g > 0 (no contact) C N g for g ≤ 0 (contact) (3.6.7)
where p c,N is the contact force in normal direction, g the gap value and C N the penalty parameter.

In case of 1D problem, recalling the Coulomb frictional contact condition, the tangential component of the contact force can be formulated as:

p c,T = 0 for g > 0 (no contact) -sign(v T )µ|C N g| for g ≤ 0 (contact) (3.6.8)
where p c,T is the contact force in tangential direction which has to be in the opposite sense with respect to the relative tangential velocity v T of the two contact surfaces and µ is the frictional contact coefficient.

The contact condition between the two bodies is represented as a spring element which stiffness value is the Penalty parameter C N . Thus, higher is the penalty parameter value and higher will be the fulfillment of the contact constraint, and vice versa. But a high value of the Penalty parameter leads to an ill conditioning of the system structural matrices. Therefore the selected Penalty parameter value should be a compromise between these two limitations. Under such definitions, the equations of motion of the system can be re-written as follows:

Mq(t) + C q(t) + Kq(t) = p ext (t) + p c (t) (3.6.9)
Replacing eq.(3.6.7) and( 3.6.8) into eq.(3.6.9) leads to:

   Mq(t) + C q(t) + Kq(t) = p ext (t) Mq(t) + C q(t) + Kq(t) = p ext (t) + C N g(t) -sign(v T (t))µ|C N g(t)| (3.6.10)
This system can now be solved by the Central difference scheme expression for the displacement (eq.(3.6.12)).

In case of no contact (g < 0), the same formula is recalled:

q i+1 ≈ M h 2 + C 2h -1 2M h 2 -K q i + C 2h - M h 2 q i-1 + p ext i (3.6.11)
In case of detected contact (g > 0), the expression is updated as:

q i+1 ≈ M h 2 + C 2h -1 2M h 2 -K q i + C 2h - M h 2 q i-1 + +p ext i + C N g i -sign(v T,i )µC N g i ) (3.6.12)
It is possible to notice the simplicity of the formulated expressions to solve at each iteration which leads to a lower computational effort with respect to Lagrange multipliers method. Thus, lower computational time is expected applying this approach to solve a contact problem. Although, attention has to be paid in selecting a suitable value for the Penalty parameter in order to ensure both a good accuracy of the solution and avoid ill conditioning of the structural matrices.

Study of Penalty coefficient

As already described in section 3.6.1, also in this section will be discussed the beam tip response to the contact interactions of the test case described in section 3.4 which in this case is computed by Penalty method.

The difference with Lagrange multipliers method is that in this case there is the presence of a user chosen parameter, the Penalty coefficient, which has to be chosen properly in order to ensure both a good accuracy and stability of the solution. Penalty parameter indeed has to be selected high enough to ensure a good accuracy of the solution, but at the same time not very high to not lead to a bad conditioning of the system structural matrices and so to an explosion of the solution.

In Fig. 3.7 and 3.8 are shown the axial and vertical displacements of the tip beam as response to the interaction with the wall. It is possible to notice how the accuracy of the solution changes with respect to the Penalty parameter. In Fig. 3.7, when the wall penetrates the beam, i.e. from ∼ = 1ms to ∼ = 15ms, the beam tip axial displacement should reach a value of -40µm to ensure a fulfillment of the contact constraint since the wall penetrates the beam tip for a maximum value of -40µm. In the obtained results it is shown how increasing the value of the Penalty parameter a better compliance with the constraint of contact is ensured. The maximum Penalty parameter to select in such conditions is equal to C N = 10 11 N/m because higher values start showing some instabilities in the solution due to the bad conditioning of the structural matrices, while lower values does not ensure a good accuracy of the solution. Master Thesis -47-Khadija Masrour 3.6. Contact algorithms

Comparison of the two contact algorithms

Once checked the convergence of the results given by both the algorithms and selected an appropriate value for the Penalty parameter in case of the Penalty method, it is opportune to compare them and the advantages and drawbacks of the two approaches. In Fig. 3.9 and 3.10 is represented a comparison of the results given by the two techniques. It is possible to notice that the Penalty parameter selected ensures a good accuracy of the solution since the Penalty method results overlap completely the results given by the Lagrange multipliers method.

Figure 3.9: Comparison of the beam's tip axial displacement computed by Lagrange multipliers method and Penalty method. With the selected value of the Penalty parameter the Penalty method results overlap perfectly those given by Lagrange multipliers method. In Tab.3.4 it is shown a ratio between the computational costs of the two algorithms for the discussed simulation. It is possible to see, as it is expected, that the Penalty method is the more advantageous one from a computational effort point of view and at the same time ensures a good accuracy of the results. Therefore, the algorithm that will be selected for the wear treatment in section 4 is the Penalty method.

Computational cost ratio [-]

Lagrange multipliers method 1 Penalty method 0.72 Table 3.4: Computational cost comparison between Lagrange multipliers method and Penalty method for the simulation shown in Fig. 3.9 and 3.10.

Model reduction

In this section Craig-Bampton reduction technique is applied on the above described simulation exploiting Penalty method since it has been shown to be more advantageous with respect to Lagrange multipliers method.

It will be checked the convergence of the solution for a different number of retained modes in order to ensure a good accuracy of the solution keeping as low as possible the number of these retained modes, thus as low as possible the computational cost of the simulation.

The idea is to retain only the contact node of the beam and condense all the other nodes following the same choice made in different studies of retaining in the blade only a subset of contact nodes and condense all the other internal nodes as shown in Fig. 3.11.

Craig Bampton method

To apply Craig-Bampton reduction method it is needed to separate the boundary degrees of freedom from the internal ones. The aim is to consider as boundary node only the contact one, while all the other nodes are internal and can be added dynamically to the static condensation as represented in Fig. 3.12. The undamped equations of motion of the system are:

Mq + Kq = p (3.7.1)
in which M is mass system matrix, K the stiffness system matrix and p the vector of loads applied, e.g. in the analyzed case are the contact forces that result from the interaction with the wall.

Computational costs

As already discussed, the principal aim of reduction techniques is the gain in the computational cost. Therefore, in Tab.3.5 is shown a comparison of the computational cost ratios registered applying Craig-Bampton on the model and using Penalty method to simulate the contact interactions. The minimal number of modes required to reach a good accuracy of the solution, as discussed in section3.7.2, is equal to 10. In this case, the computational cost is reduced of almost 1/2 with respect to the case in which no reduction is applied, which represent an interesting save.

Computational cost ratio

Chapter conclusions

This section has been dedicated entirely to the development of contact interaction algorithms as basis for the incoming wear treatment which will be discussed in the following chapter (4). The contact algorithms developed and validated are the Lagrange multipliers method and Penalty method.

The test simulation performed corresponds to the contact of the beam tip with a rigid wall that moves in both axial and transverse directions with respect to the beam axis.

In transverse direction it moves with a constant upward velocity, while in axial direction it moves with a sinusoidal law, i.e. penetrating and leaving free the beam alternatively. It has been shown that Penalty method, selecting an adapt Penalty coefficient, allows to recover an accurate solution as well as important saving in the computational time with respect to Lagrange multipliers method. Therefore it has been selected as contact treatment algorithm for the incoming chapter in which wear will taken into account. Afterwards, a reduction technique according to Craig-Bampton method has been applied on both the algorithms in order to reduce the computational costs, concluding that a minimum number of retained modes equal to 10 is required in order to ensure a good accuracy of the solution.

Chapter 4

Wear modeling

Chapter introduction

At this step of the analysis, the algorithms to simulate the contact interactions between the beam tip and the moving wall have been developed and their results have been discussed. The aim of this section is now to take into account wear phenomena in these interactions.

During the contact, because of the friction of the beam with the wall a certain amount of material will be worn out from the tip at each interaction. This quantity can be modeled by Archad's law. Since the geometry of the beam changes step by step during the simulation, at each time iteration the finite element model has to be updated. The techniques most used nowadays for this finite element model update are the re-meshing ones and the killing elements ones.

The first one consists in keeping the same number of elements used to model the beam during the whole simulation. The update of geometry occurs changing their dimensions which will be decreased to accommodate the wear. Hence in remeshing techniques the total number of degrees of freedom of the system is preserved during the whole evolution of the simulation, contrarily to the second technique in which the wear is accommodated eliminating the worn out elements progressively. As pointed out, both the techniques require a modification of the system structural matrices at each bodies interaction which make them very costful methods. The alternative technique that will be introduced, discussed and validated in this section has the principal aim of avoiding this geometrical update at each time step. The structural matrices will be indeed recovered for each amount of wear directly by a polynomial interpolation of the exact structural matrices computed in correspondence of a small subset of desired interpolating points.

Once proved the validity of the results given by this alternative technique, another important step will consist in applying reduction model technique. The application of reduction would be advantageous because the polynomial interpolation technique involves constant interpolating matrices. Thus the reduction would be applied only once on these interpolating configurations and then the solution for each wear amount would be recovered involving only these reduced matrices. The gain in computational time that could lead to would be very interesting, therefore the validity of its application will be checked in the incoming analysis.

Archard's law

At each frictional contact interaction of the beam with the rigid wall a certain amount of material is removed from the beam, thus the gap between the beam and the wall is increased. Such gap increase has to be taken into account in the implemented contact algorithms.

In order to perform this gap correction at each iteration it is needed to predict the correct amount of material removed from the beam. For this purpose different wear models have been developed during the years but the most used one because of its simplicity and at the same time reliability is the Archard's law model. According to Archard's law, the material depth removed (d w ) is proportionally dependent on the sliding distance (s), the normal contact load (p n ) and eventually the frictional coefficient (µ) if the friction is taken into account, and inversely dependent on the material hardness (H) [START_REF] Luoa | Numerical analysis of die wear characteristics in hot forging of titanium alloy turbine blade[END_REF]:

d w (T ) = K w (T ) p n s H(T ) (4.2.1)
Such dependence is controlled by the empirical coefficient K w , called Archard's law coefficient. Higher is K w and higher will be the removed material from the beam and vice versa. Therefore, Archard's law coefficient will be varied in the performed simulations in order to simulate both soft and severe wear phenomena.

A discretization in time of eq.(4.2.1) leads to:

∆d wi (T ) = µK w (T ) p ni v T i H(T ) dt i (4.2.2)
where for the i th time iteration, ∆d wi (T ) is the incremental wear amount and v T i is the transversal relative velocity between the beam and the wall.

At the end of the process the wear depth of the beam tip is:

d tot w (T ) = n i=1 ∆d wi (T ) (4.2.3)
where n is the total number of steps of the simulation.

Geometry update

As presented in the literature review, different approaches exist in wear problems to update the interacting bodies geometry at each time iteration. The principal classification leads to the differentiation between re-meshing techniques and elements elimination techniques.

In case of elements elimination approach, called also elements killing approach, the number of elements within the model changes progressively during the simulation. Since the proposed alternative wear treatment technique, which constitute the core of such master thesis research and that will be presented and discussed in section 4.4, the preservation of a constant number of elements, i.e. constant matrices dimensions during the whole simulation evolution is a fundamental requirement, therefore, the selected method to implement and to use as term of comparison is the re-meshing approach.

Available and frequently used remeshing approaches are several. In the analyzed case it has been chosen to implement both the above mentioned approaches for wear treatment, i.e. global remeshing and local remeshing. The given results by such techniques will form a term of comparison and validation for the Polynomial interpolation wear treatment technique.

Global remeshing

A first approach implemented to simulate the wear evolution and update the geometry change during the interactions of the beam with the moving wall is a global remeshing of the elements after each contact.

A global remeshing approach consists in accommodating the incremental wear registered at each iteration on all the beam elements as shown in Fig. 4.1. As already discussed, the selected approach to update the system geometry depends on the studied application, thus the same procedure can be advantageous in some cases and totally disadvantageous in other cases.

A global remeshing technique in the analyzed case it can be predicted to be totally disadvantageous from a computational cost point of view. Indeed, the geometry deformation in such application is limited to a very small depth with respect to the beam length, and so a remesh of all the remaining beam elements it is not strictly needed, but leads only to further computations and so higher memory storage, which will be demonstrated from the incoming results.

Local remeshing

Since the wear in the analyzed case is limited to a very small depth on the free extremity of the beam, then a more clever choice is a remeshing limited to such area, thus a local remeshing approach. As pointed out, different local remeshing techniques have been validated and nowadays are used in several applications.

Among the local remeshing techniques presented and described in section 2.6.1, the equally displaced nodes multi-layer remeshing approach is exploited in order to decouple the solution at each iteration from the previous remeshing history and at the same time ensure a more uniform mesh near the contact area. In Fig. 4.2 it is shown an illustration of the selected approach. The total quantity of elements needed to accommodate the wear at each iteration is given by the relation:

d wear H cr H 0 = q + r (4.3.1)
in which d wear is the total amount of wear, H cr and H 0 are respectively the critical element height and initial configuration element height, while r is the fraction of the remaining wear which will be accommodated on the first surface contact element as shown in Fig. 4.3. The critical height H cr has been selected as 60% of the initial configuration element height H 0 .

Comparison of the two technique results

In this section it is presented the response of the beam tip to the interaction with the rigid wall under the same test conditions described in section 3, i.e. the wall moves in both horizontal and vertical directions. Along the first one it follows a sinusoidal law and penetrates the beam tip for a maximum value of 40µm and along the second direction by a constant upward velocity. The frequency of the harmonic motion is doubled with respect to the previous simulation of simple contact without wear, i.e. it is set equal to 60Hz. This choice is made in order to analyze the vibratory response of two consectutive interactions because wear is taken into account and thus it is desired to analyze its effect on more than one interaction. The parameters of this simulation are summarized in Tab.4.1.

Parameter Value Unit of measure

µ 0.3 [-] C N 10 11 [N/m] t max 25 [ms] D max N 40 [µm] f u wall 60 [Hz] V v wall 176 [m/s] K w 10 -6 [m 3 /N m]
Table 4.1: Parameters for simulation of contact interactions between the beam and the wall taking into account wear.

In this simulation the worn out of material from the beam tip due to the frictional interactions with the rigid wall is considered. In Fig. 4.4 and 4.5 it is represented the response of the beam tip to these interactions with the wall in terms respectively of axial displacement and vertical displacement of the beam tip, recovering in one case the solution by a global remeshing of the worn beam (yellow line) and in the other case by a local remeshing (red line). These two solutions, which are totally superimposed as it is possible to notice in the plots, are compared with the solution in case in which no worn out of material from the beam is considered (blue line). In Fig. 4.4, already at the first interaction with the wall (from ∼ = 1ms to ∼ = 7.5ms) it is possible to notice a worn out of material which increases significantly at the second interaction (from ∼ = 17.5ms to ∼ = 24ms). Indeed the beam tip does not penetrate until 40µm as in case of no worn out of material is considered, but it reaches a lower value. Simultaneously, in Fig4.5 it is possible to see how the beam, once it is pushed up after the wall penetration, it returns back earlier with respect to the not worn out beam due to the reduction of its length.

It is important to figure out also in this case that the beam, once it is no more in contact with the wall, vibrates freely along the axial direction according to its third eigenfrequency (approximately 12500Hz), and along its vertical direction according to its first eigenfrequency (approximately 840Hz). Other important quantity to analyze during the evolution of the simulation is the quantity of the removed material from the beam step by step and compare the match of the results given by the two techniques of global remeshing and local remeshing. In Fig. 4.6 it is represented the amount of material removed from the beam at each contact interaction. It is possible to notice that when no interaction with the wall is registered, no removal of material occurs. The removal of material occurs only during the interaction phases and the results given by the two techniques match perfectly. In Fig. 4.7 it is represented the evolution of the wear, i.e. of the total quantity of material removed from the beam during the simulation evolution.

In coherence with the results presented in Fig. 4.6, when no interaction with the wall occurs, no increase in the wear is registered. The increase occurs only during the interaction phases, i.e. from ∼ = 1ms to ∼ = 7.5ms and from ∼ = 17.5ms to ∼ = 24ms. Due to this removal of material from the beam, its length is reduced, and the consequence of this phenomenon is a stiffening effect of the beam which is shown in Fig. 4.8 for the first frequency. The same stiffening effect is registered for all the other beam frequencies.

Once again, it is important to make a focus point on the match of the results given by the two techniques of global and local remeshing, to state that the given results are validated and will form a basis of comparison for the alternative technique that will be introduced and discussed in section 4.4. About the computational costs, it is expected that the local remeshing technique is the more advantageous one with respect to the global remeshing one. In the local remeshing technique indeed, as already discussed, the only modification in the system model that is performed is inherent to the last finite elements, the ones more close to the contact node. Therefore only some terms of the structural matrices will be modified and updated according to the new beam geometric dimensions, while the other elements remain totally unchanged.

In the global remeshing instead the registered wear is equally distributed on each finite element of the system leading to an update of the whole structural matrices. As expected, in Tab.4.2 it is possible to notice that the computational costs registered by the local remeshing technique are more advantageous with respect to the one correspondent to the global remeshing technique. The aim with the alternative technique that will be introduced in the next section will be to reduce further the computational costs.

Computational cost ratio [-] Global remeshing 1

Local remeshing 0.75 Table 4.2: Comparison of the computational cost ratios for the above presented simulations involving in one case a global remeshing of the beam and in the other case a local remeshing of the beam.

Proposed alternative technique

The traditional wear treatment techniques require at each time step an update of the geometry. These operations require a high computational effort and memory storage, therefore, alternative techniques is possible to notice the simplicity of such operations to perform at each iteration which will lead to a significant decrease of the computational effort for wear treatment problems.

In most of the tested cases it has been verified that a number of 3 interpolating points is sufficient to have a good accuracy of the solution. In this case the interpolated equations are reduced to: . The interpolating points are in correspondence of δ 0 = 0, δ 1 = 1/2 and δ 2 = 1 for the proposed example. δ 2 = 1 corresponds to a maximum limit value of 5mm, so the final wear δ f inal on the beam has to be lower than this threshold (green circle).

     K(δ) = c 0 (δ)K 0 + c 1 (δ)K 1 + c 2 (δ)K 2 M(δ) = c 0 (δ)M 0 + c 1 (δ)M 1 + c 2 (δ)M 2

Validation of polynomial interpolation technique

In this section will be presented the results given by the Polynomial interpolation technique of the structural matrices for the test case described in section 3.4 in which the wall is moving in both the horizontal and vertical directions, along the first one following a sinusoidal law and penetrating the beam tip for a maximum value of 40µm, and along the second direction by a constant upward The obtained results will be compared with the already presented results given by the traditional techniques of global and local remeshing obtained in section4.3.3 to check the validity of the polynomial interpolation approach in a dynamic simulation.

Comparison of the Polynomial interpolation technique results with the traditional remeshing technique results

Once checked that the proposed approach presents a good match of the frequency evolution with the one given by global remeshing technique (Fig. 4.11), in this section will be applied such approach to simulate a contact case and compare the obtained results with those given by the traditional remeshing techniques. The aim of this comparison is a full validation of the proposed approach for a beam. Once validated the technique for a beam, then it will be possible to go further with this research on more complex systems checking if it can be trusted for giving reliable results. It would lead to significant advantages because of its simplicity in the formulation and the interesting gain in computational costs that could present.

In Fig. 4.12 and 4.13 are represented the axial and vertical displacements of the beam tip applying Penalty method for the contact interactions correction.

The wear of the beam has been taken into account by the three above discussed techniques, the global remeshing, the local remeshing and the polynomial interpolation technique.

It is important to focus on the found out behaviour. Indeed, as desired, the results of the polynomial interpolation technique match perfectly the results given by the traditional remeshing techniques. The same overlap of results is verified also for the amount of worn out material from the beam step by step, the total wear evolution in time and the frequency increase effect shown respectively in As already pointed out, one of the most important gains that the polynomial interpolation could lead to, and the main reason for which it has been introduced, is the interesting gain in time that could register. In Tab.4.3 are shown the computational times encountered for the simulation plotted on Fig. 4.12-4.16. The local remeshing, as it is expected, presents a computational time lower than the global remeshing.

The interesting result here and on which it is very worth to focus is the gain in time registered by the polynomial interpolation technique with respect to both the global and local remeshing techniques. The gain in time registered for the polynomial interpolation is approximately 1/2 with respect to the global remeshing and approximately 1/3 with respect to the local remeshing.

Computational cost ratio [-]

Validation of the method for more severe contact interactions

At this point, once checked the reliability of the polynomial interpolation technique for the above performed simulation, it is interesting to check the validity also for other contact interaction configurations, eventually more severe interactions.

Therefore, the simulation discussed in this section presents much more contact interaction phases between the beam tip and the wall, in particular 15 interactions. The wall indeed is still moving in both the directions, i.e. horizontally and vertically. The difference in this new simulation is that the frequency of the sinusoidal law with which is moving along the horizontal direction is no more 60Hz but 600Hz, and the maximum penetration is no more 40µm but 300µm. Another variation performed in this simulation is the Archard's law wear coefficient, which is no more 10 -6 mm 3 /N mm, but 10 -5 mm 3 /N mm , which can be translated in a more intense wear phenomenon, i.e. higher amount of removed material from the beam. These new simulation parameters are summarized in Tab. [START_REF]Lecture notes of Theory of vibrations course[END_REF] The aim at this point is to check the overlap of the results given by the polynomial interpolation technique for these more severe interaction conditions with the results given by the global and local remeshing. This comparison is shown in Fig. 4.17-4.21 in which are represented respectively the axial displacement of the beam tip as consequence of the contact with the wall, its vertical displacement, the amount of material removed from the beam step by step, the total amount of material removed from the beam during the whole evolution of the simulation and finally the stiffening effect of the beam in which it is possible to see that the frequency increases from approximately 844.47Hz to 847.36Hz for a total wear of approximately 185µm and the behaviour is still linear. It is possible to notice that also in these severe contact interactions the results given by the polynomial interpolation technique are fully reliable since they match perfectly the already validated results given by the global and local remeshing techniques. It is possible to conclude then that this alternative technique is validated for the discussed interaction conditions of a beam model and could present an interesting research development for other configurations and more complex systems for future applications.

Study of Archard's wear coefficient

For further analysis and study, Polynomial interpolation technique has been exploited also to observe the effect of Archard's wear coefficient on the wear evolution of the beam, and thus of its induced vibrations following each contact interaction.

According to Archard's law, higher is this coefficient and higher is the amount of material removed from the part at each interaction. The simulation performed in section4.5.1 has been implemented using a wear coefficient of 10 -6 m 3 /N m. In this section the same results will compared with those obtained exploiting a wear coefficient equal to 10 -5 m 3 /N m and 10 -7 m 3 /N m, thus in the first case more intense wear phenomena are expected while in the second case lower wear phenomena are expected. This is indeed what can be observed in Fig. wearNormalConditionsKwear, where the solution correspondent to k w = 10 -5 m 3 /N m (red line) reaches a much higher amount of removed material from the beam with respect to the other two cases (yellow and orange lines). This different intensity in the wear corresponds to a different amount of stiffening effect in the beam. Higher is the reached wear and lower will be the final dimension of the beam, thus higher will be the corresponding natural frequencies which can be accurately predicted by the polynomial interpolation technique as shown in Fig. 4.23. Such variations in the amount of material removed from the beam influence also the vibrations to which the beam is subjected to as well as the time in which such interactions occur. A more worn out beam is a beam that will be subject to contact with the wall later and such penetrations will become lower and lower in the deepness. In Fig. 4.24 indeed the first contact interaction occurs at the same instant for all the three simulations, but since with k w = 10 -5 m 3 /N m the amount of material removed from the beam is high already at this first interaction, thus the intensity of the wall penetration becomes lower and the second contact interaction starts later with respect to the other two simulations (yellow and orange lines) and finishes earlier. Simultaneously, as the u displacement penetration is lower in case of the red line and occurs for a lower range of time, also the vertical vibrations are lower in the intensity and occur in a lower range of time as it is possible to observe in Fig. 4.25. In case of higher value of the wear coefficient (red line) the u penetration of the beam tip is lower because a higher value of material is removed from the beam. And also, since the beam length is decreased, the contact interaction starts later and ends earlier with respect to the other two simulations. In case of higher value of the wear coefficient (red line), as shown in Fig. 4.24, the u penetration of the beam tip is lower because a higher value of material is removed from the beam, this results in lower vertical vibrations of the beam tip, that start later and end earlier in phase with the u displacement.

Model reduction on polynomial interpolation technique

At this stage, performed the validation of the polynomial interpolation technique and stated the reliability of its results in the wear treatment phenomena for a beam, a final step planned for further decrease of the computational cost and so to increase the profit of this alternative technique is the introduction of the reduction on the build finite element model. If the model reduction technique applicability is verified for the polynomial interpolation technique, then it could be considered the workhorse for this research work. Indeed, nowadays, the application of reduction techniques directly on the traditional remeshing techniques does not present an advantage from a computational time point of view, which is contrarily the basic aim for which they are applied. This is due to the fact that at each time step the geometrical dimensions of the system change due to the wear phenomena, and so the reduction procedure would have to be reapplied or modified step by step. Contrarly, in case of polynomial interpolation the model reduction could be easily applied because the structural matrices for each wear configuration are computed interpolating constant matrices, thus the reduction could be applied on these constant matrices and the interpolation performed directly on these reduced matrices. Therefore, the interpolated structural matrices of the system recovered step by step will be directly the reduced ones leading to a substantial decrease of the model dimensions and so of the computational costs.

The reduction technique applied, as in section 3.7, is the Craig-Bampton one. The statically retained modes are those of the contact node (beam tip), and dynamically are added progressively those correspondent to the other nodes starting from the external ones and going towards the clamped node.

More modes are added dynamically, and more accurate will be the computed solution but at the same time bigger will be the model, thus a balance between such limits has to be fulfilled in order to ensure both a good accuracy of the solution and low model dimensions.

The results relative to Craig-Bampton applied to the same simulation presented in section. The different curves plotted are computed applying the reduction on the system retaining statically the contact node modes and adding progressively the other modes in the above described way. In the yellow curve, which correspond to only 3 retained modes, i.e. only those of the contact node without other dynamically added modes (each beam node presents 3 DoFs), it is possible to notice that it presents a high error with respect to the solution given by the non reduced system (blue line), which can be considered the exact solution of reference for the convergence checking.

Increasing the number of the dynamically added modes, it is possible to observe that the solution moves towards the "exact" one, which is a good sign for the performed checking analysis. Furthermore, it is possible to note also that an overlap of the reduced model solutions with the non reduced one occurs starting from at least 10 retained modes, which is the minimum number of modes to retain to recover an accurate solution.

Comparison of computational costs

Once verified that Craig-Bampton reduction method is applicable on the polynomial interpolation technique, it is worth to compare and focus on the gain of computational times that are registered at this point of the research analysis. In Tab.4.5 is shown a comparison between the computational time ratios for the non reduced model results with the reduced ones in function of the total retained modes. As already pointed out from Tab.4.

Chapter conclusions

The principal stated aim of this section is to check the validity of the Polynomial interpolation technique. To check this validity, traditional algorithms of global and local remeshing techniques have been developed to constitute a term of comparison. From the found out results it has been possible to observe how the results given by the introduced technique match perfectly the ones given by the global and local remeshing, even in very severe contact interactions. A further analysis have been addressed to check the applicability of Craig-Bampton reduction method on the introduced technique and it has been shown its validity for the analyzed cases, leading to important savings in cost. The conclusion of the conducted analysis is the validation of the technique to simulate the contact between a moving rigid body and a flexible beam. The next step for future developments will be to check the validity of this technique in case of more complex systems with respect to a beam, starting from a plate model. In the next chapter will be presented some ideas of application of this proposed approach to a plate which will be continued and further developed at LAVA laboratory for future applications to simulate the wear evolution of a blade consequently to the contact interactions with the surrounding casing.

Part II

Wear modeling in 2D structures Chapter 5

Extension of polynomial interpolation technique to 2D systems

Chapter introduction

In this chapter it will initialized the extension of the proposed alternative approach of Polynomial interpolation technique to 2D structures such as a plate. Some ideas of such extension will be proposed but they represent just an initialization in such research. The presented approaches indeed require the assumption of the wear profile shapes and basing on this assumption, different interpolating geometric configurations will be used to build a model that from the obtained results will be shown to be acceptably accurate. Such results suggest that these configuration interpolations could represent a potential alternative technique to the traditional remeshing techniques, but further research is under development in order to come out with a more robust technique that can be applicable in more general cases.

Plate model

The plate used for the incoming simulations is made of titanium alloy Ti-6AI-4V and its geometric and material properties are shown in Tab.5.1. The plate, clamped on one side and free on the other side, is modeled by finite element method with 10 elements along the length and 4 elements along the width . The FEM elements used to model the plate are Quadrilater elements with two degrees of freedom per each node, u and v as illustrated in Fig. 5.1 and its relative first 6 natural frequencies are shown in Tab.5.2. In each element a 2 × 2 Gauss integrating points are used.

Properties

Proposed approaches

In Chapter4 it has been presented and validated the proposed approach of Polynomial interpolation technique in case of a beam. In this section the aim is to extend such approach to a plate model. It will be shown that there are some proposed ideas on how to apply it in a 2D case that leads to good results, but the research is still working on developing a robust approach able to solve more general cases.

The following proposed ideas are indeed based on assuming the final expected wear shape and basing on this assumption, build appropriate interpolating matrices. Such approach can be applied in applications in which the wear profile shape can be roughly predicted according to the direction of the relative mouvements that occur between the involved parts. An example can be, as the described test case for the beam, a clamped plate in contact with a wall that penetrates the plate along its u direction by a harmonic law and moves upward along its v direction. The expected final wear shape in this case is a more or less triangular shape with more worn out the lower nodes than the upper nodes. The upper nodes indeed will be subjected to less intense contact events with respect to the lower ones, thus less intense wear phenomena leading to this more or less triangular shape. Thus, in this case, basing on this final expected wear profile shape it is possible to build appropriate interpolating matrices. In the following three sections will be presented three cases of final expected wear profile and relative interpolating matrices. But this analysis can be extended to any other assumed wear profile.

Uniform wear profile

In many applications the wear profile that evolves is more or less uniform on the contact surface, such as the example of the moving pin on the disc presented in [START_REF] Benabdallah | Finite element simulation of the wear of polyoxymethylene in pin-on-disc configuration[END_REF]. In case of the plate an illustration of such more or less uniform wear shape evolution is given in Fig. 5.2. The interpolating plate geometric configurations in case of a quadratic polynomial interpolation (3 interpolating points) are shown in Fig. 5.3. Hence by these three geometric configurations (δ 0 , δ 1 and δ 2 ) it is possible to build the correspondent interpolating structural matrices: K 0 , M 0 , K 1 , M 1 and K 2 , M 2 . The wear at each iteration is more or less constant but not equal, thus in the polynomial interpolation equations an average of the wear at each contact node is performed. This average is a weighted average since the stiffening effect induced by the border nodes is less intense than the one induced by the internal nodes. Thus, the internal nodes will have a higher weight in such average:

                                     K(δ) = np i=1 np j =i δ-δ j δ i -δ j K(δ i ) = np i=1 c i (δ)K i M(δ) = np i=1 np j =i δ-δ j δ i -δ j M(δ i ) = np i=1 c i (δ)M i δ = dw dw,max d w = nc j=1 w j d wj nc j=1 w j (5.3.1)
where d w is the average of the amount of worn out material on each contact node, d w,max is the maximum expected wear depth in the plate at the end of the simula, n c is the number of contact nodes, n p the orde of the polynom and w j the weight given to each contact node in the average of the wear depth.

The solution results obtained by such approach interpolating the structural matrices at δ 0 = 0, δ 1 = 1/2 and δ 2 = 1 with a maximum limit value of 5mm wear depth are shown in Fig. 5.4. In the simulation the final average wear reached by the plate is equal to 4.4mm and the results for the first 6 frequencies are shown. It is possible to observe how the approximated solution matches properly the exact solution obtained by global remeshing.

Triangular wear profile

In this section will be presented the proposed approach in case in which the wear evolution profile is amenable to a more or less triangular shape as illustrated in Fig. 5.5. Also in this illustrated case it has been verified that a quadratic polynomial interpolation is sufficient to have a good accuracy of the results. The interpolating geometric configurations under this wear profile shape assumption are illustrated in Fig. 5.6. In the performed simulation, once built the interpolating model, the wear on each contact node is given as input and evolves from 0mm on each contact node (unworn plate) until the following final wear profile: 0mm, 0.9mm, 2.1mm, 2.9mm and 3.95mm starting from the right contact node in Fig. 5.5, which corresponds to approximately a triangular shape but not exactly. The results given in this case are shown in Fig. 5.7. Also in this case it is possible to conclude that the proposed approximation gives good results. Indeed even the more curved shapes of the second, fourth and fifth frequencies are almost perfectly overlapped. The developed simulation under this shape assumption has been implemented with a maximum wear limit of 6mm on the more worn out edge. The assumed final wear profile starting from the less worn out edge, for each contact node, is of: 1.8mm, 2.7mm, 3.9mm, 4.8mm and 5.7mm, which is not a perfect final trapezoidal shape but almost. The obtained results of this simulation are shown in Fig. 5.10. It is possible to see that the approximated solution also in this case matches the exact one with an acceptable accuracy. But the limitation imposed by the assumption on the wear profile is very restrictive, thus a more general approach is needed to model all possible wear profiles without limitations.

Chapter conclusions

In this chapter some ideas have been presented for the extension of the presented and validated Polynomial interpolation approach in the case of a beam to the case of a plate. The proposed approach in 2D case is based on the assumption of the expected wear profile evolution of the part depending on the interested application. The assumed wear profiles studied are uniform, triangular and trapezoidal shapes, but further wear profiles can be analyzed following the same methodology. In the developed simulations it has been shown from the obtained results that the approach guarantees a good accuracy since the approximated solutions matches perfectly in almost all the cases the exact solution. But it still presents this limit due to the priori assumption of the expected wear profile shape which represents a limit of the applicability of the approach to more general cases. Therefore, the research of alternative technique for 2D cases is planned because of the enormous advantages that could lead to, as shown for the 1D case. The new alternative technique for the 2D case should be characterized

Summary and discussion

In this master thesis research an alternative approach to deal with wear problems have been introduced and discussed. The principal aim is to replace the traditional techniques for wear treatment. These one indeed require a continuous and direct update of the structural matrices to take into account the changed geometric configuration because of the continuously removed material from the contact surfaces as result of contact and sliding interactions. This continuous direct update of the structural matrices leads to relatively high computational effort and memory storage making wear problems very costful simulations. In order to accomplish to such purpose a 1D structure model have been exploited to validate the possibility of application of the proposed approach ensuring reliable and accurate results. Only the validation of such technique for 1D structures that can be extended to 2D structures and eventually 3D structures. Instead of keeping memorized the whole history of wear evolution and thus updating the structural matrices of the model at each iteration according to the re-meshed or eventually removed elements, depending on if a re-meshing or killing elements technique is exploited, the basic concept of the proposed technique is to compute such matrices only in a small subset of geometric configurations with a different amount of wear, called interpolating configurations or interpolating points. Then recover the structural matrices in all other geometric configurations in function of the amount of applied wear by a simple polynomial interpolation of such interpolating configurations. The interpolating structural matrices will be computed only one time for the whole simulation while at each iteration only a simple algebraic sum of such matrices according to the polynomial interpolation formulation will be performed, leading to significant savings in computational time and memory storage since there is no need to memorize the whole wear evolution history of the part but just the total amount of wear for each contact node at each iteration. It is admissible to think that such approximations could lead to totally inaccurate solution, therefore, such master thesis research has been conducted to verify if this approach could lead to important results and the possibility of application of such technique to real wear problems or not. To conduct this study the following steps have been followed:

• Contact analysis of a flexible clamped beam with a rigid moving wall. From the widely available literature, the range of penetration to simulate the contact of a blade with the surrounding casing ranges from 0µm to a maximum of 50µm, thus, a maximum penetration of 40µm has been selected for a first simulation. The rigid wall moves both horizontally penetrating the beam, and vertically with a constant upward velocity selecting a value that respects the relative rotational velocity of a compressor engine blade with the surrounding casing.

• The algorithms implemented to simulate the contact interactions are both Lagrange multipliers method and Penalty method. An analysis of both the methods results have been made in order to select the more suitable one for the afterward analysis. Penalty method with a penalty coefficient of 10 11 N/m have been selected because of the good accuracy of the given results and the save in time that allows to reach with respect to Lagrange multipliers method. Furthermore, an analysis of the frictional coefficient have been performed to study the effect of such value on the transversal vibrations of the beam.

• Once analyzed the beam vibrations in both normal and transversal directions under the condition of coulomb frictional contact, wear has been added and modeled by Archard's law.

• Since, as already mentioned, wear consists in a continuous removal of material from the contact surfaces of the interacting parts, the geomtery changes continuously. To take into account this changement, global and local re-meshing techniques have been implemented and relative results have been compared to ensure the validity of the obtained results. These results have been compared also with the solution in case of wear is not considered showing how the vibrations change due to this continuous decrease in the beam length. Since re-meshing technique is a widely used approach and fully validated, its results are used as term of comparison and validation for the proposed technique.

• The proposed technique algorithm of polynomial interpolation have been implemented. A first check have been conducted to simulate the eigenfrequencies evolution in function of wear that consumes the beam until a maximum depth of 5% of its initial length. Because of wear, the length of the beam decreases and thus it is subject to a stiffening effect in which the natural frequencies increase progressively. These stiffening effect results obtained by global re-meshing technique have been compared with those of the proposed polynomial interpolation. It has been observed that a quadratic polynomial interpolation, thus only three interpolating configurations are sufficient to obtain an accurate solution. The results indeed match perfectly the "exact" one obtained by the remeshing technique, giving hopeness of the possibility of implementation and exploitation of such technique for future wear simulations and applications.

• The technique has been furthermore verified applying it to the case of contact interaction simulation of the clamped beam and the moving wall. In a first simulation the maximum penetration of the wall into the beam tip is fixed at a maximum of 40µm, with a wear coefficient of 10 -6 m 3 /N m and frictional coefficient of 0.3 to reflect the contact interaction conditions of a real blade with the surrounding casing. The frequency of penetration of the wall into the blade is of 60Hz, thus two contact interactions occur during a simulation of 25ms. A second test condition have been performed simulating much more severe contact interactions to verify if the proposed approach gives reliable and accurate results even in this case. In this second test the wall penetrates the blade until a depth of 300µm and the frequency of such wall penetrations is of 6000Hz. Thus, 15 contact interactions occur during the simulation of 25ms. In both the cases of normal contact conditions and severe contact interactions it has been observed that the proposed approach gives reliable and accurate results. The other important result is relative to the computational cost. As expected, the proposed approach allows to register significant savings in the computational time leading to conclude that, thanks to the accurate results and important advantages in terms of computational effort and memory storage, such technique can be an intresting focus for future research and developments not only for simulating the wear contact interactions between rotor and static components but also for other applications.

• A last step in the validation of the proposed approach for 1D structures is the application of Craig-Bampton reduction model. Another important advantage of the proposed approach is the possibility of exploiting the benefit of reduction model techniques even for wear problems.

It allows indeed to apply the reduction only on the interpolating structural matrices decreasing significantly their dimensions and performing the polynomial interpolation directly on this reduced model. The further decrease in computational cost have been analyzed and pointed out.

• Once conducted and completed the analysis for the possibility of application of the proposed technique for 1D structures, an initialization have been performed for 2D structures giving some inputs and suggestions for future works and research. The basic idea of the proposed approach for the extension of the polynomial interpolation technique to 2D structures is based on the studied application and thus on the expected wear shape evolution. Indeed, the principal complexity of applying such approach to 2D models is that more nodes are subjected to wear simultaneously. Thus the developed technique has to take into account the wear evolution on each node and build interpolating matrices that takes into account such information. A first attempt of solution is an assumption of the wear profile that could evolve on the contact parts according to the kind of application and direction of the contact interactions. Depending on such conditions it is possible to have a rough idea on the wear profile evolution that could be classified in a more or less uniform shape, triangular, trapezoidal, parabolic, ecc. According to the expected profile shape assumption it is possible to build appropriate interpolation matrices. It has been shown from the obtained results that such approach gives acceptably accurate results. But this shape assumption represents a limitation for this approach and thus more developments and researches are planned to give light to a more robust approach able to model simultaneously all wear profile shapes without limitations. Nevertheless, the accuracy of the obtained results under such limitation suggest that the possibility of application of this technique to 2D structures and eventually 3D structures is high and this encourages a further focus and research on the presented concept.

Recommendations for future work

As stated, the formulation of such approach for wear problem treatment has to be splitted into 3 main steps, first the validation for 1D structures, secondly the relative extension to 2D structures and eventually, once verified the reliability of the obtained results, a final extension to 3D structures. In this thesis it is covered and analyzed the application of this method to 1D strcutures and the dicussion in case of 2D structures is initialized giving open possibility of research and future developments. The proposed inputs in this thesis for the extension to 2D structures have the limitation on the wear profile shape evolution, reducing the generality of the approach. Thus, further improvements and developments have to be performed in order to build a robust model for 2D structures able to represent a good competitor for traditional wear treatment problems. An important aspect that should

1 . 1

 11 p c,N normal contact forces vector p c,T tangential contact forces vector q generalized DoFs vectorq C condensed modes in GI q R retained modes in GI v T,rel relative velocity vector v Schematic illustration of jet engine components [33]. . . . . . . . . . . . . . . . . . 2.1 Variation of static (a) and total (b) efficiency with blade-jet speed ratio at design pressure ratio for four tip clearances of 1.2, 3.1, 5.0 and 8.0, expressed in percentage of the passage height [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Cross-section of a compressor stage showing the rotating blade and the location of the abradable seal [3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Evolution of tip clearance following contact interaction events. Three phases are shown: a) Prerubbing; b) Rub interaction; c) Clearance after rub.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Beam element degrees of freedom. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Beam finite element model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Simplified beam model and rigid wall interactions. . . . . . . . . . . . . . . . . . . 3.4 Shows blade tip axial displacement (red line) due to the wall motion (black line), which starts with a clearance of 10µm and reaches a maximum penetration of 40µm into the beam. Once the wall leaves free the beam ( ∼ = 15.1ms), this one starts vibrating axially with a frequency of approximately 12500Hz which corresponds to the its third frequency. These results are computed by Lagrange multipliers method applying a frictional coefficient of 0.3. . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Shows blade tip vertical displacement (red line) due to the wall motion (black line),

  4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Comparison of the beam tip axial displacement as response to the frictional contact with the moving wall computed applying Penalty method and selecting different values of Penalty parameter C N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Comparison of the beam tip vertical displacement as response to the frictional contact with the moving wall computed applying Penalty method and selecting different values of Penalty parameter C N . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . 50 3.12 Boundary and internal degrees of freedom on the beam for Craig-Bampton reduction technique. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.13 Comparison of the beam tip axial displacement computed applying Craig-Bampton reduction technique retaining a different number of modes. . . . . . . . . . . . . . . 53 3.14 Comparison of the beam tip vertical displacement computed applying Craig-Bampton reduction technique retaining a different number of modes. . . . . . . . . . . . . . . 53 4.1 Global remeshing approach. a) Unworn beam. b) Globally remeshed worn beam to accommodate the total amount of wear d wear . . . . . . . . . . . . . . . . . . . . . . 57 4.2 Local remeshing approach. a) Unworn beam. b) Locally remeshed worn beam to accommodate the total amount of wear d wear following the equally displaced nodes multi-layer remeshing approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.3 Shows the elements dimensions after an equally displaced nodes multi-layer remesh. 58 4.4 Comparison of beam tip axial displacement in case of no worn out of material is considered (blue line) and in case in which wear phenomena are considered, applying in one case a global remeshing (yellow line) of the beam and in the other case a local remeshing (red line). Global and local remeshing techniques results are totally overlapped. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.5 Comparison of beam tip vertical displacement in case of no worn out of material is considered (blue line) and in case in which wear phenomena are considered, applying in one case a global remeshing (yellow line) of the beam and in the other case a local remeshing (red line). Global and local remeshing techniques results are totally overlapped. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.6 Amount of material worn out at each contact interaction, the results of the yellow line are computed applying a global remeshing of the beam, while the red line is obtained applying a local remeshing. . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.7 Total amount of material worn out during the evolution of the whole simulation of 25ms, the results of the yellow line are computed applying a global remeshing of the beam, while the red line is obtained applying a local remeshing. . . . . . . . . . 62 4.8 Stiffening effect of the beam due to the worn out material. The results of the yellow line are computed applying a global remeshing of the beam, while the red line is obtained applying a local remeshing . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.9 Illustration of Polynomial interpolation approach for wear treatment. The blue line represents the exact solution of the frequency evolution computed by global remeshing in function of δ, which represents the reached wear amount on the structure over the maximum expected wear value (d w /d w,max ). This maximum expected wear value can be arbitraly selected. The proposed approach consists in interpolating this exact solution in a small subset of points (red crosses) and recovering the solution in all other points by a polynomial interpolation. . . . . . . . . . . . . . . . . . . . . . 64 4.10 Illustration of proposed approach in case of Quadratic interpolation (3 interpolating points). The interpolating points are in correspondence of δ 0 = 0, δ 1 = 1/2 and δ 2 = 1 for the proposed example. δ 2 = 1 corresponds to a maximum limit value of 5mm, so the final wear δ f inal on the beam has to be lower than this threshold (green circle). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Shows the first 6 frequencies evolution of the beam for a maximum wear value of 3.9mm. The "exact" solution given by global remeshing technique is compared with the solution obtained by the proposed approach of Polynomial interpolation technique (Fig.4.10). It is possible to notice how this approximated solution matches perfectly the exact solution in case of all the 6 frequencies. The same match is found for all the other frequencies of the beam. . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Comparison of beam tip axial displacement computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Comparison of beam tip vertical displacement computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Amount of material worn out from the beam at each contact interaction computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Total amount of material worn out from the beam during the evolution of the whole simulation of 25ms computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . . . . 4.16 Comparison of the frequency evolution of the beam computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Comparison of beam tip axial displacement for more severe contact conditions computed by the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . . . . . . . . . . . . . . . . . . 4.18 Comparison of beam tip vertical displacement for more severe contact conditions computed by the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . . . . . . . . . . . . . . . . . . 4.19 Shows, for more severe contact interactions, the amount of material worn out from the beam at each interaction with the wall computed by the Penalty method taking into account the wear by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . 4.20 Shows, for more severe contact interactions, the total wear during the simulation evolution computed by the Penalty method and taking into account the wear by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . . . . . . . . . . . . . . . . . . . . . . 73 4.21 Shows, for more severe contact interactions, the frequency increase in function of the worn out material during the simulation evolution computed by the Penalty method and taking into account the wear by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line). . 74 4.22 Comparison of the wear evolution in simulations that involves a different value of the wear coefficient, 10 -7 m 3 /N m for the yellow line, 10 -6 m 3 /N m for the orange line and 10 -5 m 3 /N m for the red line. It is possible to notice the increase of the wear intensity as the value of the wear coefficient increases. Indeed the red line reaches an amount of wear of approximately 32µm while the yellow line is bounded into approximately only 2µm, while the orange line reaches an intermediate value. . . . . 75 4.23 Different intensity of the stiffening effect are encountered by the beam in case of the different values of the wear coefficient. Higher is this value and higher is the stiffening effect (red line) because it corresponds to a higher amount of material removed from the beam. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.24 Axial displacement of the beam tip in function of different values of the wear coefficient.

Fig. 3 .

 3 9 and 3.10. . . . . . . . . . . . . . . . . . 3.5 Computational cost ratios comparison for Penalty method in function of the retained modes according to Craig-Bampton reduction technique. The case of 10 retained modes is the selected one because, as discussed in section 3.7.2, it is the one that ensures a good convergence of the solution. . . . . . . . . . . . . . . . . . . . . . . 4.1 Parameters for simulation of contact interactions between the beam and the wall taking into account wear. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Comparison of the computational cost ratios for the above presented simulations involving in one case a global remeshing of the beam and in the other case a local remeshing of the beam. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Comparison of the computational cost ratios encountered with the global remeshing, local remeshing and polynomial interpolation techniques. . . . . . . . . . . . . . . . 4.4 Simulation parameters in case of severe contact interactions. . . . . . . . . . . . . . 4.5 Computational cost ratios comparison for Penalty method applying polynomial interpolation technique for wear treatment and retaining a different number of modes according to Craig-Bampton reduction method. . . . . . . . . . . . . . . . . . . . . 5.1 Plate geometrical and material properties. . . . . . . . . . . . . . . . . . . . . . . . 5.2 The first six natural frequencies of the plate. . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Schematic illustration of jet engine components [33].

Figure 2 . 1 :

 21 Figure 2.1: Variation of static (a) and total (b) efficiency with blade-jet speed ratio at design pressure ratio for four tip clearances of 1.2, 3.1, 5.0 and 8.0, expressed in percentage of the passage height [2].

Figure 2 . 2 :

 22 Figure 2.2: Cross-section of a compressor stage showing the rotating blade and the location of the abradable seal [3].

Figure 2 . 3 :

 23 Figure 2.3: Evolution of tip clearance following contact interaction events. Three phases are shown: a) Prerubbing; b) Rub interaction; c) Clearance after rub.It is highlighted the advantage of having abradable linings where the induced clearance increase is locally concentrated, whereas, in case no abradable linings are used the induced gap increase interests the whole engine section[START_REF] Fois | Investigation and caharacterizaton of the wear mechanisms of abradable compressor linings[END_REF].

2 . 4 .

 24 Starting from the fan and moving towards the high pressure compressor stage the engine temperature increases significantly, thus different composites with different properties are used.

Figure 2 . 4 :

 24 Figure 2.4: Shows the most used abradable materials in modern engines according to the compressor stage and thus to the environment temperature conditions [3].

  2.6.

Figure 2 . 6 :

 26 Figure 2.6: Stability Newmark plane [4].

Figure 2 . 15 : 13 (

 21513 Figure 2.15: Single-layer re-meshing technique illustration [24]. a) Initial configuration with selected the elements of the contact zone that will be interested by re-meshing. b) Remeshing steps during the wear simulation involving single layer remeshing at each iteration

Figure 2 . 18 :

 218 Figure 2.18: Body configuration in the unworn (a) and worn case (b). The nodes in the worn case are displaced according to the linear displaced nodes technique [26].

Figure 2 . 20 :

 220 Figure 2.20: Mapping of the physical coordinates of a quadrilateral element into natural coordinates.

Figure 2 . 21 :Table 2 . 1 :

 22121 Figure 2.21: Shows the first four product rules for Gauss Integration with equal number of points along the two directions [32].

3. 2 .

 2 Beam model for each node, u, v and θ, i.e. the node axial displacement, transversal displacement and the plane rotational angle (Fig.3.1).

Figure 3 . 1 :

 31 Figure 3.1: Beam element degrees of freedom.

Figure 3 . 2 :

 32 Figure 3.2: Beam finite element model.

Figure 3 . 3 :

 33 Figure 3.3: Simplified beam model and rigid wall interactions.

Figure 3 . 4 :

 34 Figure 3.4: Shows blade tip axial displacement (red line) due to the wall motion (black line), which starts with a clearance of 10µm and reaches a maximum penetration of 40µm into the beam. Once the wall leaves free the beam ( ∼ = 15.1ms), this one starts vibrating axially with a frequency of approximately 12500Hz which corresponds to the its third frequency. These results are computed by Lagrange multipliers method applying a frictional coefficient of 0.3.

Figure 3 . 5 :

 35 Figure 3.5: Shows blade tip vertical displacement (red line) due to the wall motion (black line), which starts with a clearance of 10µm and reaches a maximum penetration of 40µm into the beam.Once the wall leaves free the beam ( ∼ = 15.1ms), this one starts vibrating vertically with a frequency of approximately 840Hz, which corresponds to the its first frequency. These results are computed by Lagrange multipliers method applying a frictional coefficient of 0.3.

Figure 3 . 6 :

 36 Figure3.6: Sensitivity of beam tip vertical displacement to the frictional coefficient. Higher is the frictional coefficient and higher is the frictional force, thus, higher is the energy with which the beam is pushed upward leading to higher vertical displacements. The results are obtained applying Lagrange multipliers method with contact conditions described in section3.4.

Figure 3 . 7 :

 37 Figure 3.7: Comparison of the beam tip axial displacement as response to the frictional contact with the moving wall computed applying Penalty method and selecting different values of Penalty parameter C N .

Figure 3 . 8 :

 38 Figure 3.8: Comparison of the beam tip vertical displacement as response to the frictional contact with the moving wall computed applying Penalty method and selecting different values of Penalty parameter C N .

Figure 3 . 10 :

 310 Figure 3.10: Comparison of the beam's tip vertical displacement computed by Lagrange multipliers method and Penalty method. With the selected value of the Penalty parameter the Penalty method results overlap perfectly those given by Lagrange multipliers method.

Figure 3 . 11 :

 311 Figure 3.11: Meshed jet-engine compressor blade in which only a subset of contact nodes is retained in CB reduced model while all other nodes are condensed [15].

Figure 3 . 12 :

 312 Figure 3.12: Boundary and internal degrees of freedom on the beam for Craig-Bampton reduction technique.

Figure 3 . 13 :

 313 Figure 3.13: Comparison of the beam tip axial displacement computed applying Craig-Bampton reduction technique retaining a different number of modes.

Figure 3 . 14 :

 314 Figure 3.14: Comparison of the beam tip vertical displacement computed applying Craig-Bampton reduction technique retaining a different number of modes.

Figure 4 . 1 :

 41 Figure 4.1: Global remeshing approach. a) Unworn beam. b) Globally remeshed worn beam to accommodate the total amount of wear d wear .

Figure 4 . 2 :

 42 Figure 4.2: Local remeshing approach. a) Unworn beam. b) Locally remeshed worn beam to accommodate the total amount of wear d wear following the equally displaced nodes multi-layer remeshing approach.

Figure 4 . 3 :

 43 Figure 4.3: Shows the elements dimensions after an equally displaced nodes multi-layer remesh.

Figure 4 . 4 :

 44 Figure 4.4: Comparison of beam tip axial displacement in case of no worn out of material is considered (blue line) and in case in which wear phenomena are considered, applying in one case a global remeshing (yellow line) of the beam and in the other case a local remeshing (red line). Global and local remeshing techniques results are totally overlapped.

Figure 4 . 5 :

 45 Figure 4.5: Comparison of beam tip vertical displacement in case of no worn out of material is considered (blue line) and in case in which wear phenomena are considered, applying in one case a global remeshing (yellow line) of the beam and in the other case a local remeshing (red line). Global and local remeshing techniques results are totally overlapped.

Figure 4 . 6 :

 46 Figure 4.6: Amount of material worn out at each contact interaction, the results of the yellow line are computed applying a global remeshing of the beam, while the red line is obtained applying a local remeshing.

Figure 4 . 7 :

 47 Figure 4.7: Total amount of material worn out during the evolution of the whole simulation of 25ms, the results of the yellow line are computed applying a global remeshing of the beam, while the red line is obtained applying a local remeshing.

Figure 4 . 8 :

 48 Figure 4.8: Stiffening effect of the beam due to the worn out material. The results of the yellow line are computed applying a global remeshing of the beam, while the red line is obtained applying a local remeshing

( 4 . 4 . 2 )

 442 In Fig.4.11 are shown the solutions given by this proposed approach for the first 6 frequencies of the beam. The interpolating points are 3 and are the two extremity points and the middle point (δ = 0, δ = 1/2 and δ = 1) for a maximum wear of 5mm. Once build the interpolating matrices K 0 and M 0 , K 1 and M 1 , K 2 and M 2 by these three interpolating points, the wear simulated in the beam evolves until 3.9mm (Fig.4.10). Thus, it is given as input into eq.(4.4.2) a wear that evolves from 0mm until 3.9mm and at each iteration new matrices K(δ) and M(δ) are build, allowing to recover the relative eigenfrequencies evolution which are the ones shown in Fig.4.11 limiting to the first 6 frequencies. It is possible to notice that the approximated solution recovered by the polynomial interpolation technique matches perfectly the "exact" solution given by global remeshing technique.

Figure 4 . 10 :

 410 Figure 4.10: Illustration of proposed approach in case of Quadratic interpolation (3 interpolating points). The interpolating points are in correspondence of δ 0 = 0, δ 1 = 1/2 and δ 2 = 1 for the proposed example. δ 2 = 1 corresponds to a maximum limit value of 5mm, so the final wear δ f inal on the beam has to be lower than this threshold (green circle).

Figure 4 . 11 :

 411 Figure 4.11: Shows the first 6 frequencies evolution of the beam for a maximum wear value of 3.9mm. The "exact" solution given by global remeshing technique is compared with the solution obtained by the proposed approach of Polynomial interpolation technique (Fig.4.10). It is possible to notice how this approximated solution matches perfectly the exact solution in case of all the 6 frequencies. The same match is found for all the other frequencies of the beam.

Fig. 4 .

 4 14, 4.15 and 4.16.

Figure 4 . 12 :

 412 Figure 4.12: Comparison of beam tip axial displacement computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 . 13 :

 413 Figure 4.13: Comparison of beam tip vertical displacement computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 . 14 :

 414 Figure 4.14: Amount of material worn out from the beam at each contact interaction computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 . 15 :

 415 Figure 4.15: Total amount of material worn out from the beam during the evolution of the whole simulation of 25ms computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 . 16 :

 416 Figure 4.16: Comparison of the frequency evolution of the beam computed using the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 . 17 :

 417 Figure 4.17: Comparison of beam tip axial displacement for more severe contact conditions computed by the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 . 18 :

 418 Figure 4.18: Comparison of beam tip vertical displacement for more severe contact conditions computed by the Penalty method and taking into account the wear evolution by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 .

 4 Figure 4.19: Shows, for more severe contact interactions, the amount of material worn out from the beam at each interaction with the wall computed by the Penalty method taking into account the wear by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 . 20 :

 420 Figure 4.20: Shows, for more severe contact interactions, the total wear during the simulation evolution computed by the Penalty method and taking into account the wear by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 .

 4 Figure 4.21: Shows, for more severe contact interactions, the frequency increase in function of the worn out material during the simulation evolution computed by the Penalty method and taking into account the wear by the following three techniques: global remeshing (yellow line), local remeshing (orange line) and polynomial interpolation (red line).

Figure 4 . 22 :

 422 Figure 4.22: Comparison of the wear evolution in simulations that involves a different value of the wear coefficient, 10 -7 m 3 /N m for the yellow line, 10 -6 m 3 /N m for the orange line and 10 -5 m 3 /N m for the red line. It is possible to notice the increase of the wear intensity as the value of the wear coefficient increases. Indeed the red line reaches an amount of wear of approximately 32µm while the yellow line is bounded into approximately only 2µm, while the orange line reaches an intermediate value.

Figure 4 . 23 :

 423 Figure 4.23: Different intensity of the stiffening effect are encountered by the beam in case of the different values of the wear coefficient. Higher is this value and higher is the stiffening effect (red line) because it corresponds to a higher amount of material removed from the beam.

Figure 4 . 24 :

 424 Figure 4.24: Axial displacement of the beam tip in function of different values of the wear coefficient.In case of higher value of the wear coefficient (red line) the u penetration of the beam tip is lower because a higher value of material is removed from the beam. And also, since the beam length is decreased, the contact interaction starts later and ends earlier with respect to the other two simulations.

Figure 4 . 25 :

 425 Figure 4.25: Vertical displacement of the beam tip in function of different values of the wear coefficient. In case of higher value of the wear coefficient (red line), as shown in Fig.4.24, the u penetration of the beam tip is lower because a higher value of material is removed from the beam, this results in lower vertical vibrations of the beam tip, that start later and end earlier in phase with the u displacement.

  4.5.1 are shown in Fig.4.26 and 4.27 in terms of axial and vertical displacement of the beam tip.

Figure 4 . 26 :

 426 Figure 4.26: Beam's tip axial displacement computed applying Craig-Bampton reduction technique on Penalty method and using the polynomial interpolation technique for wear treatment; different number of retained modes have been compared with the non reduced model solution.

Figure 4 .

 4 Figure 4.27: Beam's tip vertical displacement computed applying Craig-Bampton reduction technique on Penalty method and using the polynomial interpolation technique for wear treatment; different number of retained modes have been compared with the non reduced model solution.

Figure 5 . 1 :

 51 Figure 5.1: Plate finite element model.

Figure 5 . 2 :

 52 Figure 5.2: Illustration of a more or less uniform wear shape evolution on a plate.

Figure 5 . 3 :

 53 Figure 5.3: Interpolating plate geometric configurations in case of a quadratic polynomial interpolation under the assumption of uniform wear profile evolution.

Figure 5 . 5 :

 55 Figure 5.5: Illustration of a more or less triangular wear shape evolution on a plate.

Figure 5 . 4 :

 54 Figure 5.4: Comparison of exact stiffening effect for the first 6 frequencies with respect to the solution obtained by the proposed polynomial interpolation technique in case of almost uniform wear evolution shape on the plate.

Figure 5 . 6 :

 56 Figure 5.6: Interpolating plate geometric configurations in case of a quadratic polynomial interpolation under the assumption of triangular wear profile evolution.

Figure 5 . 7 :

 57 Figure 5.7: Comparison of exact stiffening effect for the first 6 frequencies with respect to the solution obtained by the proposed polynomial interpolation technique in case of almost triangular wear evolution shape on the plate.

Figure 5 . 9 :

 59 Figure 5.9: Interpolating plate geometric configurations in case of a cubic polynomial interpolation under the assumption of almost a trapezoidal wear profile evolution.

Figure 5 . 10 :

 510 Figure 5.10: Comparison of exact stiffening effect for the first 6 frequencies with respect to the solution obtained by the proposed polynomial interpolation technique in case of almost trapezoidal wear evolution shape on the plate.

  Illustration of a more or less uniform wear shape evolution on a plate. . . . . . . . . 84 5.3 Interpolating plate geometric configurations in case of a quadratic polynomial interpolation under the assumption of uniform wear profile evolution. . . . . . . . . . . . 85 5.5 Illustration of a more or less triangular wear shape evolution on a plate. . . . . . . . 86 5.4 Comparison of exact stiffening effect for the first 6 frequencies with respect to the Coordinates and weights of Gauss integration rule for a number of integrating points p equal to 1, 2 and 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Beam geometrical and material properties. . . . . . . . . . . . . . . . . . . . . . . .

	List of Tables
	2.1

. . . . . . . . . . . . 77 4.26 Beam's tip axial displacement computed applying Craig-Bampton reduction technique on Penalty method and using the polynomial interpolation technique for wear treatment; different number of retained modes have been compared with the non reduced model solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.27 Beam's tip vertical displacement computed applying Craig-Bampton reduction technique on Penalty method and using the polynomial interpolation technique for wear treatment; different number of retained modes have been compared with the non reduced model solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.1 Plate finite element model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.2 solution obtained by the proposed polynomial interpolation technique in case of almost uniform wear evolution shape on the plate. . . . . . . . . . . . . . . . . . . . . 87 5.6 Interpolating plate geometric configurations in case of a quadratic polynomial interpolation under the assumption of triangular wear profile evolution. . . . . . . . . . . 88 3.2 The first six natural frequencies of the beam. . . . . . . . . . . . . . . . . . . . . . . 3.3 Parameters for simulation of simple contact interactions without wear. . . . . . . . . 3.4 Computational cost comparison between Lagrange multipliers method and Penalty method for the simulation shown in

Table 3 .

 3

		Symbol Value
	length	l	100 [mm]
	section radius	R	6 [mm]
	section area	A	113.04 [mm 2 ]
	density	ρ	4500 [kg/m 3 ]
	modulus of elasticity E	113.8 [GPa]
	Poisson's ratio	ν	0.3 [-]

1: Beam geometrical and material properties.

Table 3 . 3

 33

	Parameter Value Unit of measure
	µ	0.3	[-]
	t max	25	[ms]
	D max N	40	[µm]
	f u wall	30	[Hz]
	V v wall	176	[m/s]

[START_REF] Fois | Investigation and caharacterizaton of the wear mechanisms of abradable compressor linings[END_REF]

. Initially a clearance of 10µm is set between the wall and the blade tip (Tab.3.3). : Parameters for simulation of simple contact interactions without wear.

Table 3 .

 3 5: Computational cost ratios comparison for Penalty method in function of the retained modes according to Craig-Bampton reduction technique. The case of 10 retained modes is the selected one because, as discussed in section 3.7.2, it is the one that ensures a good convergence of the solution.

	[-]

Table 4 . 3 :

 43 Comparison of the computational cost ratios encountered with the global remeshing, local remeshing and polynomial interpolation techniques.

	Global remeshing	1
	Local remeshing	0.74
	Polynomial interpolation 0.53

Table 4 .

 4 .4. 4: Simulation parameters in case of severe contact interactions.

	Parameter Value Unit of measure
	µ	0.3	[-]
	C N	10 11	[N/m]
	t max	25	[ms]
	D max N	300	[µm]
	f u wall	600	[Hz]
	V v wall	176	[m/s]
	K w	10 -5	[m 3 /N m]

Table 4 .

 4 3, the polynomial interpolation technique without the application of the model reduction allows a gain of approximately 1/2 with respect to the global remeshing technique and of approximately 1/3 with respect to the local remeshing technique. Applying the reduction on the polynomial interpolation technique it is possible to increase this gain in time further and further. As discussed in section 4.6, the minimum number of modes to retain to ensure an accurate solution is equal to 10. With this reduced model dimensions the gain in time is approximately 1/2 with respect to the full model. This gain in time is without any doubt advantageous and interesting to point out to highlight the advantage of these results and to incite future researches and developments on it. 5: Computational cost ratios comparison for Penalty method applying polynomial interpolation technique for wear treatment and retaining a different number of modes according to Craig-Bampton reduction method.

		Quadratic interpolation technique
		computational cost ratios [-]
	No reduction applied 1
	23 retained modes	0.73
	18 retained modes	0.65
	10 retained modes	0.54
	8 retained modes	0.50
	5 retained modes	0.44
	4 retained modes	0.42
	3 retained modes	0.39

Table 5 .

 5 1: Plate geometrical and material properties.

		Symbol Value
	length	l	100 [mm]
	width	w	10 [mm]
	thickness	th	1 [mm]
	density	ρ	4500 [kg/m 3 ]
	modulus of elasticity E	113.8 [GPa]
	Poisson's ratio	ν	0.3 [-]

Table 5 .

 5 1 st natural frequency 848.84 2 nd natural frequency 5133.24 3 rd natural frequency 12383.15 4 th natural frequency 13722.56 5 th natural frequency 25462.08 6 th natural frequency 37274.01 2: The first six natural frequencies of the plate.

Khadija Masrour

Master Thesis

Acknowledgements Khadija Masrour

The partitioned system into boundary and internal degrees of freedom is:

The corresponding eigenvalue problem is:

Defining the vector of inertia forces as p = ω 2 Mq, then eq.(3.7.9) becomes:

Adding the assumption that the masses associated to the internal system degrees of freedom are zero or negligible, then it is possible to neglect the condensed inertia forces (p I ≈ 0) leading to:

The vector of nodal displacements can be re-written as:

Thus, the matrix of the static modes on the boundary is:

It is possible to notice that it is totally equivalent to Guyan-Irons reduction matrix(obtained in sec-tion2.3). Indeed Guyan-Irons method can be considered a particular case of Craig-Bampton reduction method in which only the boundary static modes are considered.

The representation of the solution can be completed adding the complementary n I = n -n B independent modes obtained by solving the vibration problem of the substructure fixed on its boundary:

The corresponding eigenmodes X I of this last equation allows to construct the second part of Craig-Bampton reduction matrix:

To operate a reduction only a subset m of such eigenmodes is retained with m << n:

Craig-Bampton reduction matrix is then:

The new stiffness and mass matrices are then given by:

Results and comparison of computational costs

In this section it will be shown the results obtained applying Craig-Bampton to the above discussed simulation. In Tab.3.4 it is shown that applying Penalty method leads to a save in computational effort with respect to Lagrange multipliers method. Further increase of this computational saving can be reached decreasing the model dimensions, thus applying a reduction technique such as Craig-Bampton. In Fig. 3.13 and 3.14 it is represented the solution of the contact interaction with the wall in terms of axial and vertical displacements of the beam tip in function of a different number of retained modes.

The results obtained applying the model reduction are compared with those obtained in case no reduction is applied to the system (blue line). In Fig. 3.14, which shows the beam tip vertical displacement, it is possible to notice even without zoom how retaining only three modes (yellow line), i.e. retaining only the contact node modes which are the boundary modes without any dynamically added internal modes, the solution is far from the exact one (blue line).

From the zoom in Fig. 3.13 and 3.14 it is possible to figure out that the number of retained modes required to have a good match of the reduced solution with the exact one is equal to 10 modes. Therefore it is possible to conclude that this is a minimum number of retained modes required to ensure a good accuracy of the solution and so to consider for the incoming analysis.

are potentially desired.

In this research work an alternative technique is presented and discussed. The basic idea of the approach is, instead of computing directly at each geometric configuration the new structural matrices, recover them by an approximation in some strategic points. An illustration of such idea is given in Fig. 4.9. ). This maximum expected wear value can be arbitraly selected. The proposed approach consists in interpolating this exact solution in a small subset of points (red crosses) and recovering the solution in all other points by a polynomial interpolation.

The structural matrices are then build as follows:

where n p is the order of the polynom, d w the amount in length of the worn out material, d w,max the maximum expected wear on the structure so that the real final wear is lower than this limit, K i and M i the interpolating stiffness and mass matrices in n p + 1 beam wear configurations.

In correspondence of each interpolating point the structural matrices K i and M i are computed exactly by any traditional re-meshing technique. Then the solution at any other point correspondent to a different amount of wear d w is recovered by these proposed polynomial interpolation equations. It

Trapezoidal wear profile

A final proposed shape is a more or less trapezoidal wear profile as illustrated in Fig. 5.8. The interpolating geometric configurations of the plate are shown in Fig. 5.9. In this case indeed, since a more complex shape is assumed, 4 interpolating points are needed to recover an accurate solution, thus a cubic polynomial interpolation.

Appendix A Algorithm flowchart