Jikaï Guillaume

Benoît, Didier, Jing, Lucas Les Futurs Docteurs Alexandre

Nicolas

Samy Romains

Keywords: Modélisation numérique, méthode semi analytique, anisotropie, revêtement, contact

Merci à Daniel pour sa confiance et son encadrement au cours de cette thèse ainsi que pour m'avoir laissé autant de liberté.

Merci à Marie pour ses recommandations et son soutien dans les moments les plus délicats.

Thank you to Tim for his kindness, his hospitality and his being my jury's president. Merci aux rapporteurs d'avoir consacré du temps à la relecture et à la correction de cette thèse.

Résumé

Le fretting et l'usure sont des problèmes récurrents dans le domaine de l'aéronautique. Les contacts aube/disque au niveau du compresseur ou de la turbine haute pression des moteurs d'avion, par exemple, sont soumis à d'importantes sollicitations de fretting à de fortes températures. L'enjeu des industriels est d'optimiser la durée de vie de ces composants et d'être capable de prévoir l'amorçage de fissures.

Afin d'améliorer la tenue des pièces, des revêtements sont utilisés pour les protéger. Leurs propriétés mécaniques et matériaux ont un impact direct sur le contact et la durée de vie. Les choix de matériaux, du nombre de couches, de l'épaisseur, de l'ordre sont donc primordiaux.

De par leur composition (fibres, mono-cristaux), leur élaboration (extrusion) ou leur mode de déposition, l'hypothèse de considérer des matériaux homogènes isotropes s'avère trop réductrice. L'anisotropie est un paramètre important à prendre en compte au niveau du dimensionnement. Les matériaux composites sont de plus en plus utilisés dans l'aéronautique.

Dans cette optique, cette thèse a pour objectif l'étude du comportement des matériaux homogènes anisotropes, en s'intéressant à l'influence des principaux paramètres mécaniques caractéristiques d'un matériau afin de mieux appréhender leurs effets. On s'attardera sur le module de Young (ou module d'élasticité), le module de Coulomb (ou module de cisaillement) et le coefficient de Poisson, et leurs valeurs selon les différentes directions.

Comme attendu, le module de Young dans la direction normale au contact joue un rôle prépondérant dans la détermination du profil de pression. Néanmoins, l'influence du module de Young dans le plan tangent au contact n'est pas à négliger, il modifie aussi la forme de l'aire de contact. L'orientation du matériau par rapport au contact est par conséquent un paramètre à prendre en considération, il peut directement atténuer ou accentuer l'effet du module de Young dans une direction privilégiée. Les module de Coulomb et coefficient de Poisson ont aussi été analysés. Il en résulte qu'ils influent significativement sur le contact.

Ces résultats se confirment dans le cas d'un massif revêtu, à la différence que les effets du revêtement et du substrat peuvent se compenser. L'impact des propriétés du revêtement sera d'autant plus important que celui-ci sera épais.

L'échelle du contact par rapport aux matériaux utilisés importe aussi sur les profils de pression. Une comparaison entre le modèle anisotrope homogène et un modèle isotrope hétérogène a été réalisée. A l'échelle mesoscopique, le composite est composé d'une matrice avec des fibres qui induisent des pics de pression alors qu'à l'échelle macroscopique, le matériau composite est perçu comme un matériau homogène, les profils de pression sont lissés. ji and s (N) Orientation angle of the material Notation

List of figures

M Matrix M M Conjugate of M M -1 Inverse of M M T Transpose of M u Vector u u k,i
Partial derivative of u k by x i u k,i j Partial derivative of u k by x i and by x j ix

Introduction

Modeling of engineering problems is becoming more complex when trying to reduce the gap between the model and the actual application. It means that less restrictive assumptions should be formulated, or in other words more physics should be implemented in the model. Among the challenges to succeed in such a task, the material properties should be considered accurately. Supposing the material is isotropic is not enough. For most composite and mono-crystal materials, their composition or the elaboration and manufacturing processes imply that it exists one or two main directions or even a general anisotropy. Moreover, coatings are often used to prevent or control wear. Coatings do not have, generally, the same properties than the substrate and may have various thicknesses. The influence of the anisotropy orientations (in the coating and in the substrate) has to be taken into account to better predict the distribution of the contact pressure and the subsurface stress-field in order to optimize the service life of industrial components.

The purpose of this work is the extension of the semi analytical model developed at the laboratory Lamcos to analyse mechanical contacts between homogeneous anisotropic materials.

The first chapter presents the context and motivations. This thesis is part of an industrial project which aims to develop innovative coatings with a good temperature resistance. Several methods have proven their efficiency in mechanical contact. Semi analytical methods have the privilege to be fast and accurate methods.

The second chapter sets out the theoretical and numerical part of this work. The elementary solutions for an anisotropic material subjected to a unit point load are demonstrated. The influence coefficients, which link the loading to the displacements or to the strains, for an anisotropic material, are obtained using Green's functions. These coefficients are validated through a comparison performed under Hertzian loading between results obtained for an isotropic material using Hertz solutions and those corresponding to an orthotropic material based on the Finite Element Method (FEM).

The third chapter deals with parametric studies. A rigid sphere is loaded normally against an anisotropic half space, coated or not. The main material characteristics are considered, as the Young's modulus along the principal directions, the Poisson's ratio or the Coulomb's modulus. The orientation of the material is also studied. The stress field is analysed for different cases of anisotropy.

Finally, the fourth chapter presents a comparison with an isotropic heterogeneous model, which takes explicitly into account the fibers in composite materials. The homogenization method, to obtain the homogenized properties of an equivalent material Introduction from an isotropic heterogeneous material, is described. The pressure profile for different cases and the stresses are analysed.

2

Chapter 1

Contact mechanics

This chapter is an introduction to contact problems encountered in aeronautic with anisotropic materials.

Modeling in contact mechanics is briefly exposed. The friction and wear are recurring problems in the field of mechanics and especially in aeronautics. The project provides innovative solutions to address the problems of fretting, fatigue and fretting wear at high temperature for aeronautical structures such as fasteners blade / disk in compressors and turbines at high temperature, fasteners mast / motor, bearings and knuckles bleed air systems. These issues are very diverse in terms of thermomechanical stresses and contact geometry.

It is now necessary to develop coatings resistant to high temperatures and with good behavior in sliding friction under load while having a wear rate controlled to provide good tribological properties. In order to remain competitive, it is important to reduce the number of unscheduled maintenance operations, to reliable tribological scenarios, to pool specific and innovative resources and to take into account an eco-design into the project. The purpose of this partnership is manifold:

-development and optimization of innovative coatings and processes, -delay the boot up in fretting / fatigue, minimize wear, -taking into account the extreme thermomechanical conditions, -evolution of tribological tests at high temperature, -advanced modeling of life in fretting / fatigue. The project consists of five major industrial groups (Airbus, Liebherr, SKF, Snecma, Turboméca), four medium-size enterprises (APS, CITRA, Mecaprotec, Orapi) and six laboratories and research organizations (ARMINES, CIRIMAT, ENISE, LAMCOS, LMI, LTDS).

The LaMCoS, laboratory of contact and structure mechanics, is a joint research unit INSA Lyon / CNRS, whose main research interests are tribology, structural mechanics, system dynamics, damage and fatigue. The team involved in this project is the team MSE (solid and damage mechanics) whose activity is intended to explain and predict behavior, damage and failure in the mechanical parts subjected to extreme stress, under development or in service. This team has developed over the last twenty years some contact models based on semi-analytical methods that rely on basic analytical solutions. These models are based on numerical techniques adapted to specific tribological applications.

The project is organized into seven tasks. The LaMCoS is involved in three of these tasks: Task 2 (development of tools for characterization and modeling), task 5 (elementary tribological characterization of new coatings) and task 6 (specific tribological characterization: technology demonstrators and modeling). The contribution of the MSE team will consist of:

-carry out theoretical and numerical developments necessary to apply these models to the problem of high temperature coated contacts, -validation of these models compared to the experience, -compare and analyze the performance of these methods versus the finite element method.

Applications

The mechanical connections are omnipresent elements in aerospace structures. Their role is to introduce the degrees of freedom in rotation in the mechanism in which they are integrated. They can eliminate the degrees of hyperstatism in mechanical assemblies and intervene also at the kinematic assembly of mechanical subgroups. The functionality of the bond is conditioned by the good behaviour of its interface. Generally, the two main problems on connections in service are:

-the jamming of the pieces caused by unfavourable slip conditions leading to the transmission of unwanted effort in the mechanical environment of the connection, -a significant wear of the pieces, causing an abnormality in the functioning of mechanism.

The service life of the connection is related to the behaviour of the contact. Also, it should be selected for this interface a tribological material couple suitable for the application in order to minimize friction, jamming and wear and thus to optimize the lifetime of the connection. Contacts between blade and disk at the high pressure compressor (CoHP) of the turbojet engine (see Fig. 1.1) are subjected to stresses of fretting at temperatures between 200 ˚C and 450 ˚C. The materials used in this part of the engine are composed with titanium. The blades at the contact area are coated with a specific material layer. Figure 1.2 shows the prints due to fretting that occurs on the contact surfaces of the blade in the CoHP of a high power engine. Note that the four pictures of this part come from the industrial groups of this project. Contacts between blade and disk of the high pressure turbine (TuHP) are composed by nickel on nickel. Operating temperatures vary between 500 ˚C and 750 ˚C. The specificity of TuHP joints lies in their design of type "foot fir" with a double contact on both sides of the blade. Currently, the contact area of blades in TuHP are not coated. Figure 1.3 shows the prints due to fretting that occurs in the contact zone of the disk in the TuHP of helicopter's engine. The connection between the shaft and the engine has the distinction of being subjected to high temperatures. Figure 1.4 shows the damage related to fretting in the ball joints. This environmental context can lead to a loss of functionality resulting in a jamming at the spherical interface of the ball and socket joint. This product can undergo two types of loading: a radial loading and an axial loading, while allowing movements of rotation. The most important movements in terms of amplitude occur during engine starting and shutting off, with a rotation angle amplitude up to 1.5 degres. In flight, the joint will suffer from engine vibrations, these movements causing the fretting at the joint. Currently, the temperature reached in the environment of these joints during the flight is 160 ˚C. Changes made to the new generations of engines should raise this temperature up to 570 ˚C. These joints must withstand and be efficient at these extremely hot temperatures but also they must be used on the ground in very cold region (-70 ˚C). The joints of the valves located in a hot zone around the aircraft engine in the platform are also subjected to these problems. The thermal and vibration environment around aircraft engines becoming more severe. The air-conditioning systems, and in particular the valves which permit the extraction of the air on the engines, are among the first victims of the increasing severity in service. These joints are used to support heavy loads while allowing rotational movements. Parts are subject to severe constraints in term of pressure, amplitude of movement and also temperature (160 ˚C to 570 ˚C). The current trends in engine design aircraft (increased engine power) suggest a significant increase of temperature in the immediate vicinity of the joints. It is therefore essential to anticipate this problem.

The development of new coatings aimed at protecting the part more efficiently, with an improved resistance to high temperatures and good tribological properties (wear and friction), using lighter and more environmentally friendly new technologies.

The main aim in using these new coatings is an extension of the lifespan that would lead to a reduction in maintenance costs and benefits of designing with reductions in volume and mass. Indeed, they delay the initiation of fretting, fatigue, wear and thus improve the reliability and lifetime. We expect a multiplication of their current life by two or even by four for some components, with, as a consequence, the same lifetime extension of the coated parts. It is also possible to estimate a minimum of 25% earned over a cycle disassembly -reassembly and divide by two the number of inspections of certain areas in the life of an engine. Greater load capacity in the ball joints allow a reduction in the size of the bearing, thereby reducing the weight of the whole equipment of the order of 2%.

Contact mechanics

Methods used in contact mechanics

Analytical solutions

Hertzian contact Contact mechanics dates back to 1882 with the original work of Heinrich Hertz On the contact of elastic solids [HER 82] . Hertz solved the problem involving contact between two elastic bodies with curved surfaces, called the Hertzian contact. Surfaces are supposed continuous and non-conforming. The following assumptions are retained:

-the contact zone is elliptical, -the problem does not account for friction, -half-spaces are considered.

-the area of contact must be much smaller than the characteristic dimensions of the bodies so that stresses are concentrated in the contact region and are not altered by boundary conditions.

Because half-spaces are considered, theory of elasticity in elastic half-spaces can be used. Hertzian theory is quite restrictive because of those assumptions but is still relevant and provides a description of the contact pressure, contact dimension, displacements and stresses within the mating bodies.

Non-Hertzian contacts

Additional complications arise when some or all these assumptions are violated and such contact problems are usually called non-Hertzian. Those non-Hertzian contacts have been extensively studied and are presented in [ JOH 87]. Many solutions exist for conforming geometries, when contacts are not elliptical, even if solutions are still based on half-space theory. Some solutions exist for contact against a sheet or a shell. Finally, many solutions are given when geometric discrepancies are considered, such as sharp edges, etc. However, many of those solutions are only two-dimensional. Most of three-dimensional solutions are limited to axisymetric geometries or other particular geometries. When Hertzian assumptions are removed, solutions often have to be found using cutting edge mathematics. It is the case when considering « singular integral equations » used by Muskhelishvili 66] proposed a theory of elastic contact mechanics of rough surfaces which is today the foundation of many theories in tribology (friction, adhesion, thermal and electrical conductance, wear, etc.).

Non-elastic behavior or non-homogeneous bodies

Analytical solutions for a uniform thickness coating do exist, if both the coating and the substrate are elastic homogeneous and isotropic [MEI 68]. Integral methods are also used [GLA 80]. Non-elastic behaviors, such as plasticity, are also studied. However, most analytical studies are limited and consider a perfectly plastic behavior. Indentation process is particularly studied and solutions exist for a conical, spherical and pyramidal tip indenter [ JOH 87]. But spatial displacements in the media are assumed to be radial. It is clear that fully analytical methods are limited when considering inelastic behaviors and non-uniform coatings. and Mindlin [MIN 49] considered a sphere normally and tangentially loaded. Tangential force is lower or equal to the limit fixed by the Coulomb's law in fully-sliding conditions. Because the Coulomb's law must be observed at each point within the contact area, sliding will appear at the edge of the contact creating a slip annulus. Mindlin & Deresiewicz [MIN 53] also studied this spherical contact when tangential force is a linear function of the normal force. The very popular Cattaneo-Mindlin Concept has been extented to any two dimensional geometry by Ciavarella [CIA 98a, CIA 98b], but analytical solutions are still limited when dealing with any three-dimensional geometry and loading.

Frictional contacts

Numerical methods

Finite element method

The finite element method (FEM) is a numerical technique for determining approximate solutions of partial differential equations (PDE) as well as of integral equations. It is largely used in mechanics for complicated domains changing over time, and several modern FEM packages also include specific components (dynamic, thermal, electromagnetic, plasticity [PEI 05], viscosity [CHA 96], friction [RAO 99], etc). Anciaux [ANC 09] presented a multi scale coupling method to address contact problem.

But computation costs can be important, depending on the level of accuracy required. When contacts and inclusions are considered, using a finer mesh at the two body interface, imposed by the gradients that are expecting, drastically increases the computing time. Moreover, for multi scale problems, it is difficult to model at pertinent scales both the structure (macro scale) and the contact area (micro scale) while dealing with wear or cracking surface. For these reasons, FEM was not selected for solving complex contact problems.

Semi-analytical methods When analytical solutions are ways too complicated, it is possible to discretize the full-problem in a sum of elementary problems. Then, the solution is the numerical summation of analytical solutions for each elementary problem. Semi analytical methods, or SAM, have been intensively used in contact mechanics but numerical techniques can be different from an author [KAL 90] (Newton-Raphson algorithm) to an other [JAE 04] (Gauss-Seidel Algorithm). The multigrid techniques which enable fine mesh sizes coupled with local refinement techniques can deal with macroscopic and microscopic aspects, for heterogeneous material [

The contact solver

The model developed at the LaMCoS is based on the numerical summation of elementary analytical solutions (also called semi analytical method or SAM). The Finite Element (FE) method is used here only to validate SAM solutions.

Several methods can be used for contact simulation of anisotropic materials. SAM main advantage is the small computing time compared to the one required using FE method that is however widely used for many contact problems.

The contact problem solution in contact mechanics consists in finding the actual contact area, the contact pressure, the shear and slip distributions. The origin of all the theory is the famous paper of Heinrich Hertz, which gives the solution of the elastic contact between two ellipsoidal bodies without friction.

The semi analytical method consists in the summation of elementary solutions known analytically. One of the difficulties is the derivation or the identification of these elementary analytical solutions, such as the well known Boussinesq and Cerruti solutions in isotropic elasticity, available in [ JOH 87]. The framework of the three dimensional problem is simplified here by assuming the contact between one anisotropic elastic half space and a rigid body. The contact area is small in comparison to the dimensions of bodies justifying the assumption of half spaces. The curvatures of the surfaces in contact are also small. Each point of the surface is assigned a value of the pressure corresponding to the total load divided by the surface area (Fig. 1.5). Analytical solutions giving the contributions of normal and tangential loading assumed uniform over a single rectangular element will be used. By summation, the elastic deflection at each point within and near the contact area will be derived.

The elastic displacements are expressed by a double discrete convolution product between influence coefficients and the pressure or shear at the contact surface. The normal problem and the tangential problem in partial slip are therefore solved. The solution is performed by minimizing the complementary energy, so the contact pressure is constrained to be positive everywhere and there is no interpenetration. An algorithm is developed based on the conjugated gradient method (CGM). To accelerate the calculation, the Fast Fourier Transforms (FFT) are used to perform the double convolution product between the pressure and the influence coefficient matrix, at each iteration of the CGM (Eq. 1.1). The Consequently, no mesh reduction in areas of interest is possible.

Once the contact problem is solved, the strains in the half space are calculated.

u z (x i , y j ) = Nx ∑ k=1 Ny ∑ l=1 p(x k , y l )K(x l -x i , y k -y j ) (1.1)
General algorithm Two bodies are considered. A concentrated load or a rigid body displacement is imposed in one or several increments (in case of non linear constitutive law). The elastic contact is solved (CGM, FFT) which gives pressure distributions and elastic stresses. Small strain assumption is made. Then by a return mapping, plastic strains and residual stresses are calculated, if they exist. Residual stresses lead to residual displacements, so the surface geometry must be updated. The contact surface is modified to take into account residual displacements until convergence occurs. Finally the load or the interference is updated until end of loading. The general algorithm is presented in Fig. 1.6. Anisotropic materials are defined by the elastic stiffness tensor C i jkl , which satisfies the full symmetry C i jkl = C jikl = C kli j . This tensor links stresses to strains according to the relationship

Historical development of the code

σ i j = C i jkl ε kl (1.2)
Under symmetry constraints, the stresses and strains are reduced to σ = (σ 11 σ 22 σ 33 σ 12 σ 13 σ 23 ) T and ε = (ε 11 ε 22 ε 33 2ε 12 2ε 13 2ε 23 ) T . The materials used are often isotropic, cubic or orthotropic. Depending on these families of materials, two, three or nine parameters are necessary for defining completely the elastic tensor. A general anisotropic material has twenty one independent parameters

C =         C 1111 C 1122 C 1133 C 1112 C 1123 C 1131 C 2222 C 2233 C 2212 C 2223 C 2231 C 3333 C 3312 C 3323 C 3331 C 1212 C 1223 C 1231 Sym C 2323 C 2331 C 3131        
.

For an isotropic material, with G = E/(2(1 + ν)), or a cubic material the elastic stiffness tensor can be written as

C =           E(1-ν) (1+ν)(1-2ν) Eν (1+ν)(1-2ν) Eν (1+ν)(1-2ν) E(1-ν) (1+ν)(1-2ν) Eν (1+ν)(1-2ν) 0 E(1-ν) (1+ν)(1-2ν) G Sym G G          
.

Materials in contact

For an orthotropic material, the tensor C is defined by

C =         E 1 (1 -ν 23 ν 32 )δ E 1 (ν 21 + ν 31 ν 23 )δ E 1 (ν 31 + ν 21 ν 32 )δ E 2 (1 -ν 13 ν 31 )δ E 2 (ν 32 + ν 12 ν 31 )δ 0 E 3 (1 -ν 12 ν 21 )δ G 12 Sym G 23 G 31        
with δ = 1/(1 -ν 12 ν 21 -ν 23 ν 32 -ν 31 ν 13 -2ν 21 ν 32 ν 13 ) and ν i j /E i = ν ji /E j . E 1 , E 2 and E 3 are the Young's moduli along the three main directions, G 12 , G 23 and G 31 the shear moduli and ν 12 , ν 23 and ν 31 the Poisson's ratio.

Composite materials are designed to provide greater resistance under certain conditions. There are three main classes; composites with polymer matrix, metal matrix or ceramic matrix. For single crystals, the anisotropy comes from the organization of atoms in a unit cell. Wood and bones are anisotropic natural materials, because of their origin. Determining the mechanical constants can be done experimentally. Tensile tests can be carried out different main directions, allowing to obtain Young's moduli and Poisson's ratios. Shear tests can be performed to identify Coulomb's modulus. These parameters can change versus the temperature or even humidity in the room.

Orientation of the material Due to the anisotropy it is possible to modify materials according to local loading conditions, which allows a significant reduction of the quantitative material and therefore weight, a vital criterion in aeronautics. In this work, the depth corresponds to direction 3, which means that the surface is defined by directions 1 and 2. Note also that the main directions or axes of symmetry are not necessarily the same than those of the contact ones. In such a situation a basis change has to be carried out. Supposing C the elastic stiffness tensor in the material basis e i , the tensor C in the contact basis a m is obtained by

C mnop = Ci jkl ( e i a m )( e j a n )( e k a o )( e l a p ).
By using the Euler angles, the transformation matrix T m can be written as

T m =   cos ψ cos φ -cos θ sin φ sin ψ cos ψ sin φ + cos θ cos φ sin ψ sin ψ sin θ -sin ψ cos φ -cos θ sin φ cos ψ -sin ψ sin φ + cos θ cos φ cos ψ cos ψ sin θ sin θ sin φ -sin θ cos φ cos θ   ,
with φ the first angle of rotation around e 3 , θ the second angle of rotation around 1 -axis and ψ the third angle of rotation around a 3 , as shown in Fig. 1.7. Thereafter, we consider only one rotation around e 1 = a 1 , or φ = 0, ψ = 0 and θ any angle. The transformation matrix is simplified.

T m =   1 0 0 0 cos θ sin θ 0 -sin θ cos θ

Coated materials

Motivations Coatings can protect mechanical parts from oxidation, corrosion, temperature extremes, abrasion... They are often used to prevent or control wear which slowly modifies geometries of parts, affects the aerodynamic efficiency of engines and involves important servicing costs because of premature part changes. In many engineering applications, protective coatings are increasingly used to extend the fatigue life of mechanical components in contact and to provide low friction coefficients and wear resistance of tribological surface. In consequence, the study of protection against wear and corrosion induced by the use of coatings is an important research field in aeronautical engineering.

Characterization In many cases, coatings are used to improve surface properties of the substrate and extend the life span. They may be applied as liquids, gases or solids. The application processes for applying coatings are multiple; application of coatings in liquid form on a flat substrate or on a substrate of any shape, the electroplating (coating deposition method in ionic form on conductive coatings), the methods of vacuum deposition of thin films (physical or chimical vapor deposition, sol-gel process, epitaxy), extrusion coating, powder coating (deposition of coatings in the form of polymer melts)... Coatings do not have, generally, the same properties as the substrate. They are more or less compliant, rigid and resistant to high pressure and high temperature, with various friction coefficient. In this document, we will study the influence of the Young's modulus in the different main directions (in surface or along the depth), the Poisson's ratio and the Coulomb's modulus in the coating with an isotropic substrate or in the substrate when the coating is isotropic (Fig. 1.8). The anisotropy of material can also be oriented. The coating and the substrate do not have necessarily the same orientation.

Coatings may have various thicknesses. Depending on this parameter, the role played by the coating will change. Larger is the layer, more its properties will be predominant on the contact mechanics but the initial geometry of the part will be modified.

In order to protect components from fretting and improve performance, coatings (its properties, its thickness) have to be optimized.

In the literature

Chapter 2 A semi analytical method for contact between anisotropic materials

A contact model using semi analytical methods, relying on elementary analytical solutions, has been developed. The aim of this work is to extend this model to the consideration of the anisotropy, by quantify displacements and stresses of a coated

anisotropic material contacting both an isotropic or anisotropic material. This chapter presents the influence coefficients and their validations.

Introduction

According the semi analytical model developed by Nelias' team, the contact problem is solved for isotropic materials. The influence coefficients, which link the elastic displacement to the load, are those developed by Boussinesq and Cerruti. To study anisotropic materials, these coefficients need to be replace. It is the aim of this work. This chapter is devoted to the formulation of anisotropic influence coefficients and to their numerical expressions. First, the explicit expressions of the Green's functions for anisotropic infinite space by Ting and Lee [TIN 97] are recalled. They are used as starting point to express the half space anisotropic contribution. Functions for anisotropic half space are a sum of infinite space Green's functions and a complementary part which accounts for the free surface of the half-space, similarly to the Mindlin's superposition method. The influence coefficients for anisotropic half space are obtained with the three dimensional Green's functions for anisotropic bimaterials [PAN 00]. And starting from the three dimensional Green's functions in anisotropic trimaterials [YAN 02], these functions are derived for layered anisotropic half space. The Green's functions for anisotropic infinite space have an explicit expression whereas the complementary part needs to be integrated numerically. Finally, the validation of both cases, contact for an anisotropic half-space with or without coating, is led.

2. A semi analytical method for contact between anisotropic materials

Influence coefficients for anisotropic space

The influence coefficients, which link the loading to the displacements or to the strains, for an anisotropic material, are obtained with the Green's functions. The infinite part [TIN 97] determines these elements for the anisotropic infinite space.

Let u be the displacement vector in a three dimensional anisotropic elastic material. A concentrated force f = ( f 1 , f 2 , f 3 ) is applied at the origin of the reference frame, x = 0. The equation of equilibrium in terms of displacements u k is written as

C i jkl u k,l j = -δ(x 1 )δ(x 2 )δ(x 3 ) f i , (2.1)
where δ(x) is the Dirac function. Applying Fourier transforms to this equation with respect to x 1 , x 2 , x 3 , solving for the transformed u k , and then back transforming leads to

u = 1 (2π) 3 Q -1 (n * )fe -in * x dn * 1 dn * 2 dn * 3 (2.2)
where

Q ik (n * ) = C i jkl n * j n * l . (2.3) Equation 2.2 is reduced to u = 1 8π 2 r s Q -1 (n * )fds (2.4)
where r =| x | and the integral is taken around a unit circle | n * |= 1 on the plane normal to x. If we write u = G(x)f, (2.5)

G(x) = 1 8π 2 r s Q -1 (n * )ds.
(2.6)

where G(x) is the Green's function for the infinite space.

Let n and m be two orthogonal unit vectors in the oblique plane whose normal is the position vector x. The vectors [n, m, x/r] form a right handed triad. The unit vector n * in the oblique plane can then be represented as

n * = n cos ψ + m sin ψ, (2.7)
where ψ is an arbitrary parameter. Equation 2.3 can be written as

Q ik (ψ) = C i jkl (n j cos ψ + m j sin ψ)(n l cos ψ + m l sin ψ), (2.8)
and as Q(ψ) is periodic in ψ with periodicity ψ, Eq. 2.9 becomes

G(x) = 1 8π 2 r 2π 0 Q -1 (ψ)dψ = 1 4π 2 r π/2 -π/2 Q -1 (ψ)dψ.
(2.9) Three matrices 3 × 3 are defined with the elastic stiffness tensor C i jkl . Q, R and T are the double projections of the elastic stiffness tensor, independent of ψ.

Q ik = C i jkl n j n l , R ik = C i jkl n j m l , T ik = C i jkl m j m l .
(2.10)

Q ik (ψ) in Eq. 2.8 reduces to Q ik in Eq. 2.10 when ψ = 0. In view of Eq. 2.10, Eq. 2.8 can be written as

Q(ψ) = Q cos 2 ψ + (R + R T ) cos ψ sin ψ + T sin 2 ψ = cos 2 ψΓ(p), (2.11) with Γ(p) = Q + p(R + R T ) + p 2 T, p = tan ψ.
(2.12)

Six distinct eigenvalues p are obtained by calculating the roots of the sextic equation in p det(Γ(p)) = 0.

(2.13)

The roots are three pairs of conjugate complex numbers

p j = α j + iβ j , β j > 0 ( j = 1, 2, 3), (2.14)
with α j and β j real coefficients, p j the conjugate of p j . The determination of the roots' sextic equation is made numerically.

If we define

H[x] = 1 π π/2 -π/2 Q -1 (ψ)dψ, (2.15)
Eq. 2.9 is

G(x) = 1 4πr H[x].
(2.16)

H[x]
is independent of the direction of x, not on its magnitude. With the use of Eqs. 2.11 and 2.12, Eq. 2.15 can be written as

H[x] = 1 π ∞ -∞ Γ -1 (p)d p = 1 π ∞ -∞ Γ(p) det(Γ(p)) d p.
(2.17)

Γ is the adjoint of Γ such that Γ(p) Γ(p) = det(Γ(p))I.
(2.18) I represents the identity matrix. Let p ν be the roots with a positive imaginary part. The integral can be replaced by a contour integral over the upper half of the complex plane p. Using Cauchy's theory of residues we obtain

H[x] = 2i 3 ∑ ν=1 Γ(p ν ) det(Γ(p ν )) , ( 2
f (p) = (p -p 1 )(p -p 1 )(p -p 2 )(p -p 2 )(p -p 3 )(p -p 3 ).
(2.21)

Hence

H[x] = 2i det(T) 3 ∑ ν=1 Γ(p ν ) f (p ν ) , (2.22) where f (p) = d f (p)/d p. Since f (p 1 ) = (p 1 -p 1 )(p 1 -p 2 )(p 1 -p 2 )(p 1 -p 3 )(p 1 -p 3 ) (2.23)
and similar expressions exist for f (p 2 ) and f (p 3 ), Eq. 2.22 has the expression

H[x] = 1 det(T) Γ(p 1 ) β 1 (p 1 -p 2 )(p 1 -p 2 )(p 1 -p 3 )(p 1 -p 3 ) + ... . (2.24)
The three dots represent two more terms obtained from the first term by a cyclic permutation of the subscripts. The tensor H[x] in Eq. 2.24 is not valid for the degenerate cases p 1 = p 2 and p 1 = p 2 = p 3 . An alternative expression is presented below, which is clearly real and which remains valid for the degenerate cases.

The adjoint matrix Γ(p) is a polynomial of degree four in p

Γ(p) = 4 ∑ n=0 p n Γ (n) (2.25)
with Γ (n) the real matrices which are independent of p. H[x] can be written as

H[x] = 1 det(T) 4 ∑ n=0 q n Γ (n) .
(2.26)

The coefficients q n are defined by

q n = p n 1 β 1 (p 1 -p 2 )(p 1 -p 2 )(p 1 -p 3 )(p 1 -p 3 ) + ... . (2.27)
It remains to compute q n which should be real. The first term in Eq. 2.27 is

-p n 1 (p 2 -p 2 )(p 3 -p 3 ) 4β 1 β 2 β 3 (p 1 -p 2 )(p 1 -p 2 )(p 1 -p 3 )(p 1 -p 3 ) .
(2.28)

By writing the numerator as

-p n 1 [(p 2 -p 1 ) + (p 1 -p 2 )][(p 3 -p 1 ) + (p 1 -p 3 )], (2.29)
Influence coefficients for anisotropic space the first term becomes

-p n 1 4β 1 β 2 β 3 { 1 (p 1 -p 2 )(p 1 -p 3 ) - 1 (p 1 -p 2 )(p 1 -p 3 ) + 1 (p 1 -p 2 )(p 3 -p 1 ) - 1 (p 1 -p 2 )(p 3 -p 1 )
}.

(2.30) Similar expressions are obtained for the second and third terms in Eq. 2.27. It is then not difficult to compute the coefficients q n . The result can be summarize as

q n =      -1 2β 1 β 2 β 3 Re p n 1 (p 1 -p 2 )(p 1 -p 3 ) + ... -δ n2 -1 2β 1 β 2 β 3 Re p n-2 1 p 2 p 3 (p 1 -p 2 )(p 1 -p 3 ) + ... -δ n2 , for n = 0, 1, 2 2, 3, 4
where δ n2 = 1 if n = 2 and δ n2 = 0 if n = 2. q 2 has two alternative expressions. It can be noticed that for an isotropic material q 1 = q 3 = 0. The infinite part of the stresses is deduced from

σ i j = C i jks ε ks = 1 2 C i jks (u k,s + u s,k ).
(2.31)

Note about programming The general expression of the coefficients of det(Γ(p)),

which is a polynomial of degree six in p, and of Γ(p), a polynomial of degree four in p, is obtained with the commercial computer algebra system Maple. The difficulty consists in computing accurately the roots of these influence functions. For this, routines from Numerical Recipes, a book on algorithms and numerical analysis, are used. Results were compared to the ones obtained with the software library LAPACK (Linear Algebra PACKage) and they are similar.

Influence coefficients for anisotropic half space

Green's functions

The influence coefficients, which link the loading to the displacements or to the strains, for an anisotropic material, are obtained using the Green's functions. The calculation of the displacements u and the stresses σ anywhere in the material is carried out in two stages. An infinite part [TIN 97] determines these elements for the anisotropic infinite space and a complementary part [PAN 00] is used to correct this solution for the presence of the free surface.

u = u ∞ + u comp σ = σ ∞ + σ comp .
(2.32)

These complementary functions can be described explicitly in the Fourier domain, whereas it is more complicated to obtain their formulation in the physical domain due to the general anisotropy of the material. They are therefore obtained in the physical domain by applying the inverse Fourier transform to the explicit formulas of the Fourier domain.

A concentrated force f = ( f 1 , f 2 , f 3 ) is applied at the origin of the reference frame, which is located on the surface of the anisotropic elastic half space (x 3 ≥ 0). In the absence of body forces, the equations of equilibrium in terms of displacements u k are written as C i jkl u k,l j = 0.

(2.33)

The two dimensional Fourier transforms are applied

ũk (y 1 , y 2 , x 3 ) = u k (x 1 , x 2 , x 3 )e iy * x * dx 1 dx 2 , (2.34)
where y * = (y 1 , y 2 ) is the transformed vector, x * = (x 1 , x 2 ) and y * x * = y 1 x 1 + y 2 x 2 . In the transformed domain, Eq. 2.33 becomes for the subscript i equals 1 to 3 (the character i which is not a subscript is the imaginary unit)

3 ∑ k=1 2 ∑ α=1 2 ∑ β=1 C iαkβ y α y β ũk + i(C iαk3 +C i3kα )y α ũk,3 -C i3k3 ũk,33 = 0. (2.35)
Three matrices 3 × 3 are defined with the elastic stiffness tensor C i jkl and the vectors n and m, which form a right handed triad with the position vector x. Q, R and T are the double projections of the elastic stiffness tensor.

x =   r cos θ sin φ r sin θ sin φ r cos φ   , n =   cos θ sin θ 0   , m =   0 0 1   , y * = η cos θ sin θ , (2.36) Q ik = C i jks n j n s , R ik = C i jks n j m s , T ik = C i jks m j m s .
(2.37)

A general solution of Eq. 2.35 can then be expressed as ũk (y 1 , y 2 , x 3 ) = a k e -ipηx 3 .

(2.38)

The six distinct eigenvalues p are obtained by calculating the roots of

det(Γ(p)) = det(Q + p(R + R T ) + p 2 T) = 0.
(2.39)

They are either complex or purely imaginary (two by two conjugate) by the positive strain energy density requirement (shown by [ESH 53] and [TIN 96]). The complex eigenvectors a are obtained with Eq. 2.40. They are not trivial solutions.

[Q + p(R + R T ) + p 2 T]a = 0.
(2.40)

The traction vector t on the plane x 3 =constant and the in-plane stress vector s are deduced from t ≡ (σ 13 , σ 23 ,

σ 33 ) = (C 13ks u k,s ,C 23ks u k,s ,C 33ks u k,s ) (2.41) s ≡ (σ 11 , σ 12 , σ 22 ) = (C 11ks u k,s ,C 12ks u k,s ,C 22ks u k,s ).
(2.42)

Using the displacement solution (Eq. 2.38), the transformed traction and in-plane stress vectors can be found as tk = -iηb k e -ipηx 3 , (2.43) sk = -iηc k e -ipηx 3 .

(2.44)

The eigenvectors b are derived by Eq. 2.67

b = (R T + pT)a = - 1 p (Q + pR)a, ( 2.45)
and the eigenvectors c by c = Da.

(2.46)

The matrix D is defined by

D 1 j = C 11 jα n α + pC 11 j3 , D 2 j = C 12 jα n α + pC 12 j3 , D 3 j = C 22 jα n α + pC 22 j3 , with j, α = 1, 2, 3. (2.47)
If p j , a j and b j ( j = 1, 2, ... , 6) are the eigenvalues and the associated eigenvectors, it is defined that

Im(p j ) > 0 ( j = 1, 2, 3), p j+3 = p j a j+3 = a j , b j+3 = b j , c j+3 = c j , A = [a 1 , a 2 , a 3 ], B = [b 1 , b 2 , b 3 ], C = [c 1 , c 2 , c 3 ].
(2.48)

Supposing that the eigenvalues are distinct, it is helpful to normalize the eigenvectors so that a and b satisfy:

b i T a j + a i T b j = δ i j (2.49)
2. A semi analytical method for contact between anisotropic materials with δ i j the Kronecker delta, then the Green's functions are obtained by superposing the six eigensolutions of Eq. 2.38.

ũ(y 1 , y 2 , x 3 ) = ũ∞ -iη -1 A < e -ip * ηx 3 > GA T f t(y 1 , y 2 , x 3 ) = t∞ -B < e -ip * ηx 3 > GA T f s(y 1 , y 2 , x 3 ) = s∞ -C < e -ip * ηx 3 > GA T f (2.50) with < e -ip * ηx 3 >= diag[e -ip 1 ηx 3 , e -ip 2 ηx 3 , e -ip 3 ηx 3 ]. The matrix G is determined with G = B -1 B.
(2.51)

In order to obtain the Green's functions in the physical domain, the inverse Fourier transform is applied to Eq. 2.50.

u(x 1 , x 2 , x 3 ) = u ∞ (x 1 , x 2 , x 3 )- i 4π 2 η -1 A < e -ip * ηx 3 > GA T e -i(x 1 y 1 +x 2 y 2 ) dy 1 dy 2 f
(2.52) To handle the double infinite integrals, the polar coordinate transform (y 1 = η cos θ, y 2 = η sin θ and dy 1 dy 2 = ηdθdη) is introduced so that the infinite integral with respect to the radial variable can be carried out exactly. Thus, the final Green's functions in the physical domain can be expressed in terms of a regular line-integral over [0, 2π]. As the method is the same for the displacements and the stresses, only the displacement solution will be detailed. So, Eq. 2.52 becomes

u(x 1 , x 2 , x 3 ) = u ∞ (x 1 , x 2 , x 3 ) - i 4π 2   2π 0 ∞ 0 A < e -ip * ηx 3 > GA T e -iη(x 1 cos θ+x 2 sin θ dη   f
(2.53) Since the matrices A and G are independent of the radial variable η, the integral with respect to η can be performed analytically. Assuming x 3 = 0, the displacements are reduced to

u(x 1 , x 2 , x 3 ) = u ∞ (x 1 , x 2 , x 3 ) + 1 4π 2   2π 0 AG u A T dθ   f (2.54)
where the matrix G u is defined by

(G u ) i j = (G) i j -p i x 3 -(x 1 cos θ + x 2 sin θ) . (2.55)
Similarly, the Green's stresses can be written as

t(x 1 , x 2 , x 3 ) = t ∞ (x 1 , x 2 , x 3 ) + 1 4π 2   2π 0 BG t A T dθ   f s(x 1 , x 2 , x 3 ) = s ∞ (x 1 , x 2 , x 3 ) + 1 4π 2   2π 0 CG t A T dθ   f (2.56)
with the matrix G t as

(G t ) i j = (G) i j [-p i x 3 -(x 1 cos θ + x 2 sin θ)] 2 .
(2.57)

If the studied point is at the surface of the anisotropic half-space (x 3 = 0), the formulation is slightly different because of a singularity in the complementary part. Equation 2.53 becomes:

u(x 1 , x 2 , 0) = u ∞ (x 1 , x 2 , 0) - i 4π 2   2π 0 ∞ 0 AGA T e -iη(x 1 cos θ+x 2 sin θ dη   f (2.58)
By introducing the polar coordinate transform for the field point x 1 = r cos θ 0 and x 2 = r sin θ 0 , we obtain x 1 cos θ + x 2 sin θ = r cos(θ -θ 0 ). And by carrying out the integral with respect to η using relation ([BAR 75])

∞ 0 e -ikx dx = - i k + πδ(k) (2.59)
Equation 2.58 can be written as

u(x 1 , x 2 , 0) = u ∞ (x 1 , x 2 , 0)- i 4π 2    2π 0 AGA T - i r cos(θ -θ 0 ) + πδ(r cos(θ -θ 0 )) dθ    f (2.60) By virtue of ([BAR 75]) δ(r cos(θ -θ 0 )) = δ(θ -θ 0 ± π/2) r| sin(θ -θ 0 )| (2.61)
the Green's displacement on the surface of the half space can be finally expressed as

u(x 1 , x 2 , 0) = u ∞ (x 1 , x 2 , 0) - 1 4πr    1 π 2π 0 AGA T cos(θ -θ 0 ) dθ + i[AGA T ] θ=θ 0 ±π/2    f (2.62)
with the integral being understood as the Cauchy principal value. If x = 0, the function δ(x) is equal to 1, if not it is a null function. Green's stresses, where the integral is defined by the Hadamard finite part, can be derived as

t(x 1 , x 2 , 0) = t ∞ (x 1 , x 2 , 0) + 1 4πr 2    1 π 2π 0 BGA T cos 2 (θ -θ 0 ) dθ ∓ i d[BGA T ] dθ | θ=θ 0 ±π/2    f (2.63) s(x 1 , x 2 , 0) = s ∞ (x 1 , x 2 , 0) + 1 4πr 2    1 π 2π 0 CGA T cos 2 (θ -θ 0 ) dθ ∓ i d[CGA T ] dθ | θ=θ 0 ±π/2    f. (2.64)
There are some materials for which the general procedure outlined here breaks down, for instance isotropic materials. In an attempt to substitute materials constants corresponding to an isotropic material into the formulas for A and B, the term in the matrices are infinite. Some fully anisotropic materials may exhibit the same degeneracy. [CHO 03] have found a way to rewrite the complex variable formulation for isotropic materials into a form that is identical in structure to the Stroh formulation. This approach is very useful, because it enables to solve problems involving interfaces between isotropic and anisotropic materials. [BOW 09] deals with this problem too.

Note about programming

The general expression of eigenvectors is obtained with Maple and the eigenvalues are calculated with the routines of the Numerical Recipes. The calculation of integrals is made numerically. θ is taken from π/100 to 2π(100-1) 100 + π/100 with a step equals to 2π/100.

Cauchy principal value

In mathematics, the Cauchy principal value is a method for assigning values to certain improper integrals which would otherwise be undefined.

Considering c as a singularity of an one real variable function f and supposing that a < c < b, the following limit

lim ε→0 c-ε a f (x)dx + lim η→0 b c+η f (x)dx = L < +∞
exists and is finite. So, we say that the improper integral of f (x) on the interval exists and its value is defined by L.

If the above limit does not exist, it is possible that it exists when ε and η tend to zero simultaneously, that is if the limit

lim ε→0 ( c-ε a f (x)dx + b c+ε f (x)dx) = L

Hadamard finite part

The integral in the sense of Hadamard shows how to calculate an integral that has multiple singularities, as shown below.

If the Cauchy principal value of the improper integral

b a f (t) t -x dt exists, then the finite part of the divergent integral is b a f (t) (t -x) 2 dt = d dx b a f (t) t -x dt.
It can also be calculated using the definition

b a f (t) (t -x) 2 dt = lim ε→0    x-ε a f (t) (t -x) 2 dt + b x+ε f (t) (t -x) 2 dt - 2 f (x) ε    .

Influence coefficients for coated anisotropic half space 2.4.1 General formulation

The method is similar to this for influence coefficients for anisotropic half space. We work with the Green's functions in the Fourier domain to obtain analytical expressions.

The main steps are recalled.

In the absence of body forces, the equations of equilibrium in terms of displacements u k are written as C i jkl u k,l j = 0. Three matrices 3 × 3 are defined with the elastic stiffness tensor C i jkl and the vectors n and m, which form a right handed triad with the position vector x = (x 1 , x 2 , x 3 ) T . Q, R and T are the double projections of the elastic stiffness tensor.

Q ik = C i jks n j n s , R ik = C i jks n j m s , T ik = C i jks m j m s .
(2.65)

[Q + p i (R + R T ) + p 2 i T]a i = 0.
(2.66)

The subscript T denotes transpose of a vector or matrix. Six pairs of eigenvalues p i and eigenvectors a i are obtained by solving Eq. 2.66. Only three pairs of them are independent (A = (a 1 , a 2 , a 3 ) and p i with Im(p i ) > 0 and i = 1, 2, 3), the three others are their complex conjugates (A and

p i ). Matrices B = (b 1 , b 2 , b 3 ) and C = (c 1 , c 2 , c 3 ) are related to the matrix A by b i = - 1 p i (Q + p i R)a i , (2.67)
with the normalization relation b i T a j + a i T b j = δ i j , and

c i = D i a i , with D k ji = C 1k jα n α + p i C 1k j3 for k=1, 2, D 3 ji = C 22 jα n α + p i C 22 j3 .
(2.68) δ i j is the Kronecker delta. The matrix C is different from the fourth-rank elastic stiffness tensor C i jkl , which is always written in its component form in the text. The two matrices M and N are defined by the following expressions.

M = BA -1 , N = CA -1 .
(2.69)

The fundamental solution, also called Green's functions, are denoted by u * ji (x) for the displacement and σ * jki (x) for the stress. The concentrated force f, a vector, is applied at the origin of the reference frame, which is located on the surface of the anisotropic coating (x 3 ≥ 0). The direction of the force is indicated by the last subscript i. Solutions at field point x due to a point force can be written as

u j (x) = u * ji (x) f i , σ jk (x) = σ * jki (x) f i . (2.70)
While u is a vector and σ a tensor of second rank , u * is a tensor of the second rank and σ * is a tensor of the third rank. The boundary conditions, at the interface x 3 = Z C , require that the displacement and the traction vector, t = σ(0, 0, 1) T , are continuous. Therefore the six stress components are divided in two parts : t * the traction (out-of-plane stress) vector and s * the in-plane stress vector. The fundamental stress solutions can be expressed with the derivatives of displacement, taken with respect to the field point x, as

t * ≡ (σ * 13i , σ * 23i , σ * 33i ) = (C 13kl u * ki,l ,C 23kl u * ki,l ,C 33kl u * ki,l ), s * ≡ (σ * 11i , σ * 12i , σ * 22i ) = (C 11kl u * ki,l ,C 12kl u * ki,l ,C 22kl u * ki,l ).
(2.71)

The Green's functions in Fourier transformed domain (denoted by the tilde) can be written in the following series forms, with y 1 and y 2 the transform coordinates,

ũ * m (y 1 , y 2 , x 3 ) = ũ * (∞) m (y 1 , y 2 , x 3 ) + ∞ ∑ n=1 ũ * (n) m (y 1 , y 2 , x 3 ), t * m (y 1 , y 2 , x 3 ) = t * (∞) m (y 1 , y 2 , x 3 ) + ∞ ∑ n=1 t * (n) m (y 1 , y 2 , x 3 ), s * m (y 1 , y 2 , x 3 ) = s * (∞) m (y 1 , y 2 , x 3 ) + ∞ ∑ n=1 s * (n) m (y 1 , y 2 , x 3 ).
(2.72)

The subscript m defined the material in which the studied point is, 0 for the coating or 2 for the substrate. The infinite part of the displacement and the stress need to be carry out only for the coating, where the concentrated force is applied. In the substrate, this part is equal to zero, as no force is applied. u 0 * (∞) , t 0 * (∞) , s 0 * (∞) are obtained by [TIN 97] and the complementary parts by [YAN 02].

The displacement in the coating, when 0 < x 3 < Z C , is decomposed in two parts.

ũ * (N) 0 (x 3 ) = ũ * (N) 01 (x 3 ) + ũ * (N) 02 (x 3 ) with ũ * (N) 01 (x 3 ) = A 0 < e -ip 0 ηx 3 > A 0 -1 ũ * (N) 01 (0) and ũ * (N) 02 (x 3 ) = A 0 < e -ip 0 η(x 3 -Z C ) > A 0 -1 ũ * (N) 02 (Z C ).
(2.73)

According to the value of the order N, first order N = 1 or superior order N = 2, 3, . . ., ∞, the displacements in x 3 = 0 and

x 3 = Z C have different expressions. ũ * (1) 01 (0) = -M 0 -1 M 0 ũ * (∞) 0 (0), ũ * (N) 01 (0) = -M 0 -1 M 0 ũ * (N-1) 02 (0).
(2.74)

ũ * (1) 02 (Z C ) = (M 2 -M 0 ) -1 (M 0 -M 2 ) ũ * (∞) 0 (Z C ), ũ * (N) 02 (Z C ) = (M 2 -M 0 ) -1 (M 0 -M 2 ) ũ * (N-1) 01 (Z C ).
(2.75) 2. A semi analytical method for contact between anisotropic materials

The infinite part of displacement for x 3 = 0 and x 3 = Z C can be written explicitly as

ũ * (∞) 0 (0) = A 0 < e -ip 0 η0 > A 0 -1 (M 0 -M 0 ) -1 , ũ * (∞) 0 (Z C ) = A 0 < e -ip 0 ηZ C > A 0 -1 (M 0 -M 0 ) -1 .
(2.76)

The stresses are also expressed in two parts and derived expressly from the displacement.

t * (N) 0 (x 3 ) = t * (N) 01 (x 3 ) + t * (N) 02 (x 3 ) = -iηM 0 ũ01 * (N) (x 3 ) -iηM 0 ũ * (N) 02 (x 3 ), s * (N) 0 (x 3 ) = s * (N) 01 (x 3 ) + s * (N) 02 (x 3 ) = -iηN 0 ũ01 * (N) (x 3 ) -iηN 0 ũ * (N) 02 (x 3 ).
(2.77)

When x 3 > Z C
, the studied point is in the substrate. The displacement can be expressed as ũ * (N)

2 (x 3 ) = A 2 < e -ip 2 η(x 3 -Z C ) > A 2 -1 ũ * (N) 2 (Z C ), (2.78)
with, depending on the order

N ũ * (1) 2 (Z C ) = (M 2 -M 0 ) -1 (M 0 -M 0 ) ũ * (∞) 0 (Z C ), ũ * (N) 2 (Z C ) = (M 2 -M 0 ) -1 (M 0 -M 0 ) ũ * (N-1) 01 (Z C ).
(2.79)

The formulation for the stress is

t * (N) 2 (x 3 ) = -iηM 2 ũ2 * (N) (x 3 ), s * (N) 2 (x 3 ) = -iηN 2 ũ2 * (N) (x 3 ).
(2.80)

The transformed displacement, traction and in-plane stress tensors can be written as a sum of terms with each term having the following form

ũ * = ũ * (∞) + ∑ N iη -1 J N+1 < e -ir N η > J N • • • < e -ir n η > J n • • • < e -ir 0 η > J 0 , t * = t * (∞) + ∑ N J N+1 < e -ir N η > J N • • • < e -ir n η > J n • • • < e -ir 0 η > J 0 , s * = s * (∞) + ∑ N J N+1 < e -ir N η > J N • • • < e -ir n η > J n • • • < e -ir 0 η > J 0 .
(2.81) ũ * , t * and s * do not share necessarily the same vectors r n and tensors J n , see Tabs. 2.2 to 2.4.

By inserting expressions in Eq. 2.81 into the inverse-transform operator, the displacement u * , the stresses t * and s * are obtained. The in-plane stress s * has a similar expression Influence coefficients for coated anisotropic half space than the traction stress t * .

u * = u * (∞) + ∑ N 1 (2π) 2 ∞ 0 2π 0 η ũ * (N) e iy α (X α -x α ) dηdθ = u * (∞) + ∑ N 1 (2π) 2 ∞ 0 2π 0 iJ N+1 < e -ir N η > J N • • • < e -ir n η > J n • • • < e -ir 0 η > J 0 e iy α (X α -x α ) dηdθ (2.82) and t * = t * (∞) + ∑ N 1 (2π) 2 ∞ 0 2π 0 η t * (N) e iy α (X α -x α ) dηdθ = t * (∞) + ∑ N 1 (2π) 2 ∞ 0 2π 0 ηJ N+1 < e -ir N η > J N • • • < e -ir n η > J n • • • < e -ir 0 η > J 0 e iy α (X α -x α ) dηdθ.
(2.83) J n and r n are functions of θ but independent of η. As shown in Eq. 2.81, the transformed displacement and stress are written as the multiplication of a series of exponential functions of η and a factor of η. The double integrals are reducible to a 1D integral by carrying out the integral in η. The reduced integrals are given by

u * ji = u * (∞) ji + ∑ N 1 (2π) 2 2π 0 G (N) ji ( 1 s (N) + iδ(s (N) ))dθ,
(2.84)

t * ji = t * (∞) ji + ∑ N 1 (2π) 2 2π 0 G (N) ji (- 1 s (N)2 + iδ (s (N) ))dθ, (2.85)
where δ(k) is the Dirac delta, and the prime indicates the first derivative with respect to θ, with

G (N) ji = (J jk N+1 ) N+1 (J k N+1 k N ) N • • • (J k n+1 k n ) n • • • (J k 1 i ) 0 , s (N) = (r k N+1 ) N + • • • (r k n+1 ) n + • • • (r k 1 ) 0 + x 1 cosθ + x 2 sinθ.
(2.86)

Here, for the contact between a rigid sphere and an anisotropic substrate with a coating, the order N will be limited to N = 3, which gives a good approximation. Yang [YAN 02] studied the error for different orders and for this case, the order N = 4 does not provide any further significant contribution. This was also numerically verified using the present model. The calculation details are addressed in the following part.

As discussed just above, for the case of the contact between a rigid sphere and an anisotropic substrate with a coating, the order N will be limited to N = 3.

For 0 < x 3 < Z C , Eqs. 2.84 and 2.85 are reduced to

ũ * (1) 01 (x 3 ) = iη -1 A 0 < e -ip 0 ηx 3 > A 0 -1 (-M 0 -1 M 0 ) ũ * (∞) 0 (0) = iη -1 A 0 < e -ip 0 ηx 3 > A 0 -1 (-M 0 -1 M 0 )A 0 < e -ip 0 η0 > A 0 -1 (M 0 -M 0 ) -1 , ũ * (1) 02 (x 3 ) = iη -1 A 0 < e -ip 0 η(x 3 -Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 ) ũ * (∞) 0 (Z C ) = iη -1 A 0 < e -ip 0 η(x 3 -Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 )A 0 < e -ip 0 ηZ C > A 0 -1 (M 0 -M 0 ) -1 .
(2.87)

ũ * (2) 01 (x 3 ) = iη -1 A 0 < e -ip 0 ηx 3 > A 0 -1 (-M 0 -1 M 0 ) ũ * (1) 02 (0) = iη -1 A 0 < e -ip 0 ηx 3 > A 0 -1 (-M 0 -1 M 0 )A 0 < e -ip 0 η(-Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 )A 0 < e -ip 0 ηZ C > A 0 -1 (M 0 -M 0 ) -1 , ũ * (2) 02 (x 3 ) = iη -1 A 0 < e -ip 0 η(x 3 -Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 ) ũ * (1) 01 (Z C ) = iη -1 A 0 < e -ip 0 η(x 3 -Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 )A 0 < e -ip 0 ηZ C > A 0 -1 (-M 0 -1 M 0 )A 0 < e -ip 0 η0 > A 0 -1 (M 0 -M 0 ) -1 .
(2.88)

ũ * (3) 01 (x 3 ) = iη -1 A 0 < e -ip 0 ηx 3 > A 0 -1 (-M 0 -1 M 0 ) ũ * (2) 02 (0) = iη -1 A 0 < e -ip 0 ηx 3 > A 0 -1 (-M 0 -1 M 0 )A 0 < e -ip 0 η(-Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 )A 0 < e -ip 0 ηZ C > A 0 -1 (-M 0 -1 M 0 )A 0 < e -ip 0 η0 > A 0 -1 (M 0 -M 0 ) -1 , ũ * (3) 02 (x 3 ) = iη -1 A 0 < e -ip 0 η(x 3 -Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 ) ũ * (2) 01 (Z C ) = iη -1 A 0 < e -ip 0 η(x 3 -Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 )A 0 < e -ip 0 ηZ C > A 0 -1 (-M 0 -1 M 0 )A 0 < e -ip 0 η(-Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 )A 0 < e -ip 0 ηZ C > A 0 -1 (M 0 -M 0 ) -1 .
(2.89)

t * (N) 01 (x 3 ) = -iηM 0 ũ01 * (N) (x 3 ), s * (N) 01 (x 3 ) = -iηN 0 ũ01 * (N) (x 3 ). (2.90) t * (N) 02 (x 3 ) = -iηM 0 ũ02 * (N) (x 3 ), s * (N) 02 (x 3 ) = -iηN 0 ũ02 * (N) (x 3 ).
(2.91)

For x 3 > Z C , Eqs. 2.84 and 2.85 give

ũ * (1) 2 (x 3 ) = iη -1 A 2 < e -ip 2 η(x 3 -Z C ) > A 2 -1 ũ * (1) 2 (Z C ) = iη -1 A 2 < e -ip 2 η(x 3 -Z C ) > A 2 -1 (M 2 -M 0 ) -1 (M 0 -M 0 )A 0 < e -ip 0 ηZ C > A 0 -1 (M 0 -M 0 ) -1 , (2.92) ũ * (2) 2 (x 3 ) = iη -1 A 2 < e -ip 2 η(x 3 -Z C ) > A 2 -1 ũ * (2) 2 (Z C ) = iη -1 A 2 < e -ip 2 η(x 3 -Z C ) > A 2 -1 (M 2 -M 0 ) -1 (M 0 -M 0 )A 0 < e -ip 0 ηZ C > A 0 -1 (-M 0 -1 M 0 )A 0 < e -ip 0 η0 > A 0 -1 (M 0 -M 0 ) -1 , (2.93) ũ * (3) 2 (x 3 ) = iη -1 A 2 < e -ip 2 η(x 3 -Z C ) > A 2 -1 ũ * (3) 2 (Z C ) = iη -1 A 2 < e -ip 2 η(x 3 -Z C ) > A 2 -1 (M 2 -M 0 ) -1 (M 0 -M 0 )A 0 < e -ip 0 ηZ C > A 0 -1 (-M 0 -1 M 0 )A 0 < e -ip 0 η(-Z C ) > A 0 -1 (M 2 -M 0 ) -1 (M 0 -M 2 )A 0 < e -ip 0 ηZ C > A 0 -1 (M 0 -M 0 ) -1 .
(2.94)

t * (N) 2 (x 3 ) = -iηM 2 ũ2 * (N) (x 3 ), s * (N) 2 (x 3 ) = -iηN 2 ũ2 * (N) (x 3 ).
(2.95)

The following tables summarize the relevant coefficients to obtain the displacements.

Validations of the contact solution

The half space is defined by the coordinate system (x, y, z) where z represents the depth and so (x, y) the surface. As the material can be oriented to the contact, an other system (1, 2, 3) is defined for the material properties. If the angle θ m is equal to 0, then the direction 3 corresponds to the depth, as shown in Fig. 2.1.

Coefficients for anisotropic half space

This section deals with the validation of the influence coefficients derived for anisotropic half space. The material is assumed to be isotropic or orthotropic.

Validation for an isotropic case

The first validation is a comparison between the results obtained using the anisotropic model when dealing with an isotropic half space submitted to an Hertzian loading and the analytical solution given in Eqs. 2.96 and 2.97. Assuming a rigid indenter, with a spherical tip of radius R indenter , in contact with an isotropic half space, the Hertz pressure P Hertz and the Hertz contact radius a Hertz can be written as

P Hertz = 3P 2πa 2

Hertz

(2.96)

a Hertz = ( 3PR indenter (1 -ν 2 solid ) 4E solid ) (1/3) .
(2.97) P represents the load. E and ν are the Young's modulus and the Poisson's ratio of the isotropic half space, respectively. The material is assumed cubic, almost isotropic. It is not possible to adress isotropic materials with the anisotropic formulation based on the Stroh's formalism. It leads to a degenerate case where repeated roots occur from the characteristic equation of the material. Solution of Eq. 2.39 leads to double roots responsible for singularities. The Coulomb's modulus is equal to the isotropic Coulomb's modulus plus 0.7%. The following properties are E = 120GPa, ν = 0.3 and G = 46500MPa.

Two meshes are studied (the two first numbers refer to the directions parallel to the surface and the third one to the depth): 2. A semi analytical method for contact between anisotropic materials -a so-called standard mesh, 47 × 47 × 25 elements leading to 23 × 23 elements within the contact area, -a so-called refined mesh, 95 × 95 × 50 elements leading to 71 × 71 elements within the contact area. Note that for contact problems it is usually admitted that 10 points are usually sufficient to accurately describes a specific wavelength.

The differences between the anisotropic model results and the reference Hertz solution are:

-0.6% for the refined mesh, -2.1% for the standard mesh. The (anisotropic) numerical solution tends towards the (isotropic) analytical solution for isotropic properties. Table 2.5: The influence of the mesh and the model on the CPU time

The difference between computing time with isotropic and anisotropic models stems from the calculation of the influence coefficients, whose calculation is more complex and lengthy for the latter model. Note that the influence coefficients are computed only once for a given anisotropic material, and stored to be used for further simulations. According to Tab. 2.5 and Fig. 2.2, the standard mesh is a good compromise, it gives accurate results in a reasonable computing time with a personal computer (Intel R Core T M 2 Duo CPU T9600 @ 2.80GHz 1.59GHz, 3.48Go RAM). This standard will be used in the following parts.

Validation for an anisotropic case A finite element model has been developed with the commercial FE package Abaqus (version 6.9), in order to validate the anisotropic semi analytical method; a ball on plane contact is considered.

An isotropic sphere, with Young's modulus 10 6 times larger than the half space Young's modulus and identical Poisson's ratio, is in contact with a semi infinite homogeneous orthotropic half space. The mesh is shown in Fig. 2.3 and described in Tab. 2.6. It includes semi infinite elements CIN3D8 on the edges of the body (domain Ω 4 ) to mimic nil displacements like those obtained within the half space assumption, and quadratic elements with 20 nodes and 27 integration points (C3D20) in the contact area (domain Ω 1 ) for more precision. To decrease the computation cost, the model is discretized with other elements between this two areas. Elements in contact with semi infinite elements have to be parallelepiped, so C3D8R are used (domain Ω 3 ) and tetraedral elements (domain Ω 2 ) enable to reduce the number of elements and so the number of nodes and integration points. In order to respect the Hertz conditions, the sphere radius is 30 times larger than the contact radius. To optimize the computation time, properties of symmetry along direction 1 are used. Nil displacements are imposed at the bottom of the body. The first step consists in validating the FE mesh though a comparison of the contact solution obtained for an hertzian loading. The difference between the SAM pressure distribution (i.e. on the maximum pressure and the contact radius) and the FEM solution is smaller than 0.4%.

Then, two anisotropic cases are studied. The half space is considered orthotropic, with the following elastic properties: E 2 = 120GPa, the same Poisson's ratio ν = 0.3 and the same Coulomb's modulus G = 46500MPa.

Note that a Hertz corresponds to the half width of the contact obtained for an isotropic 2. A semi analytical method for contact between anisotropic materials

= E 2 = E 3 = 120GPa and ν = 0.3.
-First case: the Young's modulus along direction 3 (the depth) is twice larger than the Young's moduli in directions 1 and 2, E 3 = 240GPa. The FE and SAM pressure distributions are plotted in Fig. 2.4(a), in the plane x = 0. The difference on the maximum pressure is lower than 0.2%.

-Second case: the Young's modulus along direction 1 (on the surface) is twice larger than those in directions 2 and 3, E 1 = 240GPa. A difference of 2.2% for the maximum pressure is observed in Fig. 2.4(b), where the pressure is represented in the plane x = 0. As the results obtained for an anisotropy along the surface are not so close, a comparison for the finite element model between the two main directions of the surface has been carried out (see Fig. 2.5). The pressure profile is plotted in the plane x = 0 in Fig. 2.5(a) and in the plane y = 0 in Fig. 2.5(b). Two main points can be observed; results are different depending on the direction of the anisotropy along the surface (1 and 2) and results are the same for both cutting plane whereas the anisotropy along the surface has an effect on the contact area shape. After several verification, this difference is still misunderstood but we assume that it comes from the properties of symmetry used to simplify the model. By using this finite element model presented hereabove, results for an anisotropic half space are obtained with a CPU time of 144 hours on the same class of personal computer used for the SAM.

Coefficients for coated anisotropic half space

This second section aims to validate the influence coefficients for coated anisotropic half space. Two cases are studied, the first one when the coating and the substrate are isotropic and the second one when the half space is orthotropic, without coating.

The discretization of the half space using SAM is the so-called standard mesh, as described previously.

Validation for an isotropic substrate with an isotropic coating This validation consists in comparing the results obtained using the layered anisotropic model when dealing with a layered isotropic half space submitted to an Hertzian loading and the results given by [O'S 88].

E S , E C and ν S , ν C are the Young's moduli and the Poisson's ratio of the isotropic substrate and the isotropic coating, respectively. P represents the load. As in [O'S 88], the indenter radius, R indenter , was taken to be 10 times the coating thickness. Materials are assumed cubic, almost isotropic, because isotropic materials lead to double roots which create singularities. The Coulomb's modulus is equal to the isotropic Coulomb's modulus more or less 1%.

The gap between the O'Sullivan solutions and the numerical solutions for the pressure are inferior to 4%, except for the case E C = 0.25E S where the difference reaches 6%, see Fig. 2.6 in the plane y = 0. It can be concluded that the (anisotropic) numerical solution tends towards the (isotropic) analytical solution for isotropic properties, as expected.

Validation for an anisotropic half space

The validation is done by comparison with the finite element model used in the precedent section (2.5.1). An anisotropic case is hence studied. The half space is considered orthotropic, with the same Poisson's ratio and the same Coulomb's modulus. The Young's modulus along direction 3 (the depth) is twice larger than the Young's moduli in directions 1 and 2. In the numerical solution, the coating and the substrate are composed with the same material. Note that a Hertz corresponds to the contact half width obtained for an isotropic half space with E 1 = E 2 = E 3 . The FE and SAM pressure distributions are plotted in the plane x = 0 in Fig. 2.7. The error on the maximum value of the pressure is lower than 0.2%. A very good agreement is obtained for this case between SAM and FE results.

Conclusion

The influence coefficients for an anisotropic half space with or without coating are obtained with anisotropic Green's functions and by using the superposition method. They are expressed as a sum of an infinite part and a complementary part which extends the formulation to the half space case. These coefficients are deduced from the functions for anisotropic bimaterials or trimaterials. The three main points of the methods are the double projection of the elastic stiffness tensor according to the principal axes, the determination of eigenvalues and eigenvectors from an equation of degree six (solutions are in the complex domain) and the use of the Fourier transforms.

The model is validated by comparison with the Hertz solution for isotropic materials, with a FE model for anisotropic materials and with the results of O'Sullivan for layered isotropic materials. Isotropic materials can be studied by assuming the material cubic and modified lightly (of the order of 1%) the Coulomb's modulus, to avoid singularities.

The method is very efficient. The computation time is mainly required for the quantification of the influence coefficients. Note that an alternative is to compute them once for a given configuration and to store them. In Tab. 2.7, the CPU times, with or without the computation of the influence coefficients, are given. The mesh is composed by 77 × 77 × 40 elements and the contact area by 55 × 55 elements. Chapter 3

Normal contact for anisotropic materials

This chapter highlights the effects of anisotropy on the normal contact solution between anisotropic materials, the pressure distribution and the contact area (size and shape) and the corresponding stresses. 3. Normal contact for anisotropic materials

Parametric studies for anisotropic half space

The influence of the properties of an elastic anisotropic half space which is in contact with a rigid indenter, with a spherical tip, is studied here ([BAG 12]). The contact parameters, i.e. the contact area and the pressure distribution, will be more specifically investigated.

For the material properties, the depth corresponds to direction 3, which means that the surface is defined by directions 1 and 2. Nevertheless, the coordinate axes may be different from the material directions, as defined in Fig. 3.1. Except for the study of orientation, θ m is equal to 0. The parallelepiped is discretized using 97 × 97 × 50 elements and 53 × 53 elements in the contact area.

In these examples, materials are orthotropic, with the identical Poisson's ratio and the identical Coulomb's modulus along the three directions. The material data used in the numerical simulations are synthesized in Tab. 3.1. Note that all the result are normalized with respect to Hertzian parameters (a Hertz , P Hertz ) obtained for an isotropic material with ˚)

E = E 1 = E 2 = E 3 and ν = ν 12 = ν 13 = ν 23 .
E 1 /E 2 = 1 E 1 /E 2 = 1 E E 3 /E 2 = 1 E 3 /E 2 = 1 ν 12 = ν 13 = ν 23 , ν i j = ν ji .E i /E j ν G 12 = G 13 = G 23 = E 2 /(2(1 + ν 12 )) G
Table 3.1: Elastic properties of studied materials

The Young's modulus

Contact pressure and contact area The influence of E 1 and E 3 on the contact pressure distribution is shown in Figs. 3.2 to 3.4. It is observed that a change of the Young's modulus along a direction parallel to the surface (E 1 here) has a moderate effect (Fig. 3.2). The maximum contact pressure is increased by 3% when E 1 is doubled. Moreover it can be observed that the contact area is no more circular (Fig. 3.3) but becomes elliptical.

The effect of the Young's modulus along the depth (E 3 ) is more pronounced: increase of 32% of the maximum contact pressure and decrease of the contact radius by 11% when the Young's modulus is increased by a factor 2 (Fig. 3.4).

-2 -1.5 -1 -0.5 0 0.5 5 and3.6 summarize the influence of the Young's modulus on the maximum pressure. A slight modification of E 1 around E isotropic does not have a significant effect, conversely to E 3 that does have a pronounced effect. Nevertheless, a monotonic increase of the maximum contact pressure is observed when E 1 increases (see Fig. 3.5). The maximum pressure drops quickly when E 3 decreases (see Fig. 3.6).

1 1.5 2 0 0.2 0.4 0.6 0.8 1 1.2 y/a Hertz p/p Hertz E 1 =0.25 E isotropic E 1 =0.75 E isotropic E 1 = E isotropic E 1 =2.0 E isotropic E 1 =4.0 E isotropic E 1 =8.0 E isotropic

Normal contact for anisotropic materials

The material's orientation

The effect of the material's orientation in the plane normal to the contact is shown in Fig. 3.9. Here the Young's modulus of the material in direction 3 is equal to E 3 = 3 E isotropic , but the material main direction is different from that of the contact. The angle progressively to the solution corresponding to an anisotropic material along the surface, as seen previously when E 1 = E isotropic . This is confirmed in Fig. 3.10, where the maximum pressure is plotted versus the angle, for different ratios between E 3 and E isotropic . Hence, the method is valid whatever the orientation angle compared to the surface.

δ z /R indenter F n /(E isotropic R indenter 2 ) E 1 =0.5 E isotropic E 1 = E isotropic E 1 =2.0 E isotropic

The Poisson's ratio and the Coulomb's modulus

In Figs. 3.11 and 3.12, the influence of the Poisson's ratio ν and the Coulomb's modulus G are observed more precisely. The studied material is cubic. Dividing (multiplying) by 2 the Poisson's ratio increases by 6% (decreases by 14%) the contact radius. In terms of contact area it corresponds to an increase of 12% or a decrease of 30% for the contact area, respectively. A modification of the Coulomb's modulus by a factor 2 changes the

δ z /R indenter F n /(E isotropic R indenter 2 ) E 3 =0.5 E isotropic E 3 = E isotropic E 3 =2.0 E isotropic
E 3 =0.25 E isotropic E 3 =0.5 E isotropic E 3 =0.75 E isotropic E 3 = E isotropic E 3 =4/3 E isotropic E 3 =1.5 E isotropic E 3 =2.0 E isotropic E 3 =3.0 E isotropic E 3 =5.0 E isotropic

Parametric studies for layered anisotropic half space

An elastic anisotropic half space with an anisotropic coating is in contact with a rigid indenter, with a spherical tip ([BAG 13]). The influence of the material's properties are studied here. The contact parameters, i.e. the contact area and the pressure distribution, will be more specifically investigated. The depth corresponds to direction 3, which means that the surface is defined by directions 1 and 2. But the material main directions can be different from that of the contact by the angle θ m (see Fig. 3.13). We consider θ m = 0, except when the material orientation is studied. (x, y, z) corresponds to the contact reference axes and (1, 2, 3) to the material main direction system. The subscript C represents the coating and S the substrate.

E C E S y x, 1 z 3 2 θ m θ m E C E S y x, 1 z 3 2 θ m θ m

Figure 3.13: Coordinate systems

In the first part, materials for the coating are orthotropic and the substrate is cubic, with the same Poisson's ratio and the same Coulomb's modulus. The Coulomb's modulus is slightly modified relative to an isotropic material (G = G isotropic ± 1%), in order to avoid singularities. In the second part, roles are reversed, the coating is cubic, almost isotropic, and materials which defined the substrate are orthotropic.

The model is discretized using 77 × 77 × 50 elements and 53 × 53 elements in the contact area.

An anisotropic coating on an isotropic substrate

The Young's modulus The influence of E C1 and E C3 on the contact pressure distribution is analysed for E C1 /E S and E C3 /E S ratios ranging from 0.25 to 6 (or 8). Results are plotted in Figs. 3.14 to 3.16. The coating thickness is equal to the half of the Hertz contact radius, Z C = 0.5a Hertz . It is observed that a change of the coating Young's modulus along a direction parallel to the surface (E C1 here) has a limited influence (Fig. 3.14) on the maximum contact pressure with an increase by 2% when E C1 is doubled and by 6.8% when E C1 is multiplied by 6, but has a more pronounced effect on the shape of the contact area. From circular it becomes elliptical (Fig. 3.15).

The effect of changing the coating Young's modulus (E C3 ) in the depth direction is more pronounced: increase of 14% of the maximum contact pressure and decrease of the 3. Normal contact for anisotropic materials contact radius by 8% when the Young's modulus is increased by a factor 2 (Fig. 3.16).

The material's directions

The effect of the material's directions within the coating relative to the contact can be observed in Fig. 3.17 E C3 = E S and E C2 = 3E S . So the influence of small angles is important because of the effect of the anisotropy along the depth, whereas when the angle increases, the influence decreases since the effect of the Young's modulus along a direction parallel to the surface is moderate. Hence, the method is valid whatever the orientation angle compared to the surface. The coating thickness Figures 3.18 and 3.19 show the effect of the coating thickness on the contact pressure distribution. In the first case (Fig. 3.18), the Young's Modulus of the coating in direction 3 is equal to E C3 = 2E S . The case Z C = 0a Hertz corresponds to a half space without any coating (an isotropic half space) and the case Z C = ∞a Hertz to the half space whose properties are those of the coating (an orthotropic half space). When 3. Normal contact for anisotropic materials the thickness of the coating increases, the influence of the coating becomes predominant relative to the influence of the substrate as already observed for isotropic layered bodies ([KUL 10]). As E C3 > E S , the maximum contact pressure increases with the coating thickness.

In the second figure (Fig. 3.19), the normalized maximum pressure versus the coating thickness is plotted for different E C3 /E S ratios, when the coating is orthotropic (Fig. 3.19(a)), or different E C /E S ratios, when the coating is cubic almost isotropic (Fig. 3.19(b)). The maximum pressure increases quickly for small values of Z C then from Z C = a Hertz the curve tends slowly to an horizontal asymptote. When the coating is isotropic, the asymptote corresponds to (E C /E S ) 2/3 . For the case of an orthotropic coating, the asymptote is different but its determination will require further investigations.

An isotropic coating on an anisotropic substrate

The Young's modulus The influence of E S1 and E S3 , the Young's modulus of the substrate, on the contact pressure distribution is investigated here and results are summarized in Figs. 3.20 and 3.21. The coating thickness is equal to the half of the Hertz contact radius, Z C = 0.5a Hertz . It can be first observed that the effect of E S1 is limited. When E S1 is doubled, the maximum pressure increases by 0.7% only. If it is multiplied by 8, the maximum pressure is raised by 3%. Conversely the effect of the Young's modulus of the substrate along the depth, E S3 , is more pronounced. The maximum pressure is 13% higher when E S3 is twice than E S1 = E S2 . In addition the profile pressure is more rounded when the coating is anisotropic, whereas it is rather domed for the case of an isotropic coating.

The coating thickness Figure 3.22 shows the effect of the thickness of an isotropic coating on an anisotropic substrate on the pressure distribution. The Young's modulus of the substrate along direction 3 is twice the one of the coating, E S3 = 2E C . When the thickness of the coating increases, the layered half space tends to behave as an isotropic half space. Since here E C < E S3 , one may observe a decrease in the maximum contact pressure with the coating thickness. It can be also observed that from Z C > 2a Hertz , the difference between the maximum pressure of the isotropic half space and the layered half space is less than 1%, the coating erases the anisotropy effect of the substrate. It can be concluded that, for a normal and frictionless contact, the coating and/or substrate anisotropy along the direction normal to the interface has a strong influence on the contact solution. The influence of the coating increases with increasing thickness, both on the pressure distribution and on the contact area (shape and size). Conversely a change of the Young's modulus along a direction parallel to the surface has a limited effect on the maximum contact pressure, however it affects the shape of the contact area that is no longer circular but becomes elliptical. It reaches its maximum from a thickness greater than the contact radius. Therefore the elastic properties, the direction of anisotropy and the thickness of the coating have to be carefully chosen -and chosen together in a complementary manner -to efficiently protect the substrate.

Normal contact for anisotropic materials

Parametric studies for layered anisotropic half space

Stress field

This section analyses the stress field in an anisotropic material without coating in a first part and with a coating in a second part. After a validation, a parametric study is performed. Results are normalized with the Hertz parameters.

Coefficient for anisotropic half space

Validation of an isotropic case First, the comparison between the isotropic model and the anisotropic model (with a cubic, almost isotropic, material) is carried out. The half space is submitted to an Hertzian loading. The profiles of the dimensionless stress components are plotted along versus depth (z/a Hertz ) either at the contact center (x = y = 0), Fig. 3.23, or at an offset point (x = 0.5a Hertz , y = 0.25a Hertz ) in Fig. 3.24.

The numerical values are identical except for σ zz very close to the surface (z < 0.2a Hertz ) at the contact center only, and for σ yz when the dimensionless depth ranges between 0 and 0.8 whatever the surface point considered is. The reason why the numerical values diverge from the analytical solution near the surface is not clear. It looks like one bug remains in the calculation of the σ yz stress tensor component. For the σ zz component, it seems that the integral formulation becomes singular when z tends to zero for x = y = 0.

Validation of an orthotropic case Now a comparison between the FE model presented in the previous section (2.5.1) and the anisotropic SAM is done for the orthotropic case E 3 = 2E 2 (Fig. 3.25).

As for the comparison with the analytical solution in the isotropic case the error on the σ yz component remains whereas the one on σ zz is amplified. In addition a difference is also noticed for the σ xx and σ yy components, with also a slight difference between them while they should be identical since here E 1 = E 2 . Again the reason(s) why there is some differences is not clear. It will require further investigations.

Coefficient for coated anisotropic half space

Validation of an isotropic half space Coating and the substrate are here assumed anisotropic. They are made with the same material, a cubic one with G = G isotropic ± 1%. The results are presented here for two values of the coating thickness, a first one equals to zero and a second one chosen arbitrary equals to Z C = 0.5a Hertz . The numerical results obtained with SAM with the isotropic influence coefficients are considered as the reference solution. They are compared to the ones calculated with SAM with the anisotropic influence coefficients. The normalized stresses are presented for x = y = 0 in Fig. 3.26. The stresses are validated excepted for σ zz which, again here, differs near the surface (Z < 0.2a Hertz ). Note that the problem with the σ yz component vanishes when using the formulation with a coating.

Validation of an orthotropic half space

The coating is still taken in the same material than the substrate but their properties are orthotropic, with the same Poisson's ratio and the same Coulomb's modulus. For the first case (Fig. 3.27), the Young's modulus in direction 1 is twice larger than the Young's moduli in direction 2 and 3. For the second case (Fig. 3.28), E 3 is twice larger than the two other Young's moduli. Solutions obtained with the SAM are compared with solutions of the FEM presented in chapter 2. Two coating thicknesses are chosen, Z C = 0 and Z C = 0.5a Hertz .

Results are in a good agreement, except for the stress σ zz near the surface (Z C < 0.2a Hertz ). Moreover, a difference can be observed according to the value of the thickness. A small difference can also be observed for σ xx and σ yy components between the SAM and FEM numerical results.

E 1 =0.25 E 3 E 1 =0.5 E 3 E 1 =1 E 3 E 1 =1.5 E 3 E 1 =2 E 3 E 1 =4 E 3 0 0.5 1 1.5 2 2.5 3 3.5 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 z/a Hertz σ YY /P Hertz E 1 =0.25 E E 1 =0.5 E E 1 =1 E 3 E 1 =1.5 E E 1 =2 E 3 E 1 =4 E 3 0 0.5 1 1.5 2 2.5 3 3.5 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 z/a Hertz σ ZZ /P Hertz E 1 =0.25 E 3 E 1 =0.5 E 3 E 1 =1 E 3 E 1 =1.5 E 3 E 1 =2 E 3 E 1 =4 E 3 0 0.5 1 1.5 2
E 1 =0.25 E E 1 =0.5 E E 1 =1 E 3 E 1 =1.5 E E 1 =2 E 3 E 1 =4 E 3 0 0.5 1 1.5 2 2.5 3 3.5 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 z/a Hertz σ XZ /P Hertz E 1 =0.25 E 3 E 1 =0.5 E 3 E 1 =1 E 3 E 1 =1.5 E 3 E 1 =2 E 3 E 1 =4 E 3 0 0.5 1 1.5 2 2.5 3 3.5 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 z/a Hertz σ YZ /P Hertz E 1 =0.25 E E 1 =0.5 E E 1 =1 E 3 E 1 =1.5 E E 1 =2 E 3 E 1 =4 E 3
E 3 =0.25 E 1 E 3 =0.5 E 1 E 3 =1 E 1 E 3 =1.5 E 1 E 3 =2 E 1 E 3 =4 E 1 0 0.5 1 1.5 2 2.5 3 3.5 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 z/a Hertz σ YY /P Hertz E 3 =0.25 E E 3 =0.5 E 1 E 3 =1 E 1 E 3 =1.5 E 1 E 3 =2 E 1 E 3 =4 E 1 0 0.5 1 1.5 2 2.5 3 3.5 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 z/a Hertz σ ZZ /P Hertz E 3 =0.25 E 1 E 3 =0.5 E 1 E 3 =1 E 1 E 3 =1.5 E 1 E 3 =2 E 1 E 3 =4 E 1 0 0.5 1 1.5 2
E 3 =0.25 E E 3 =0.5 E 1 E 3 =1 E 1 E 3 =1.5 E 1 E 3 =2 E 1 E 3 =4 E 1 0 0.5 1 1.5 2 2.5 3 3.5 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 z/a Hertz σ XZ /P Hertz E 3 =0.25 E 1 E 3 =0.5 E 1 E 3 =1 E 1 E 3 =1.5 E 1 E 3 =2 E 1 E 3 =4 E 1 0 0.5 1 1.5 2 2.5 3 3.5 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 z/a Hertz σ YZ /P Hertz E 3 =0.25 E E 3 =0.5 E 1 E 3 =1 E 1 E 3 =1.5 E 1 E 3 =2 E 1 E 3 =4 E 1

Conclusion

The semi analytical model developed for the contact problem of elastic materials is used to study the effect of anisotropy. It is found that, for a normal and frictionless contact, the stiffness along the normal to the contact, in the coating or in the substrate, has a strong influence on the contact solution in terms of pressure distribution and contact size; an increase in E 3 leads to a higher maximum contact pressure and a smaller contact radius. Conversely a change of the Young's modulus along a direction parallel to the surface (plane (1,2)) does not significantly affect the contact pressure distribution, however the contact area is no more circular but becomes elliptical. The pressure profile is different, more or less domed, if the anisotropy affects the coating or the substrate.

The thickness of the coating changes also the pressure distribution. It should be emphasized that, when the coating thickness exceeds the contact radius, the effect of the substrate anisotropy vanishes quickly.

The performance of the method is highlighted by analysing the effect of the orientation of the material main directions compared to the surface normal.

Therefore the elastic properties, the direction of anisotropy and the thickness of the coating have to be carefully chosen -and chosen together in a complementary mannerto efficiently protect the substrate.

The stress components are also studied. They can not be all validated, some issues remain. Nevertheless, it is observed that the Young's modulus along the surface affects particularly the stress σ xx and the one along the depth has an important influence on σ zz . Results are confirmed by the ones on the influence of the material orientation.

Introduction

The aim of this part is to compare the anisotropic homogeneous model, presented in this thesis, to an isotropic heterogeneous model, developed in parallel by J. Leroux in the framework of his PhD thesis [LER 13]. The latter considers an isotropic half space, also called the matrix, with isotropic heterogeneities which form fibers. The study of composite materials is then possible at the mesoscopic scale. In order to compare both approaches, macroscopic and mesoscopic, properties of the heterogeneous material need to be homogenized to obtain an equivalent homogeneous material.

Brief introduction to the model with inclusions

During his PhD, Leroux [LER 13] worked on the semi analytical model and improved it in order to take into account inhomogeneities and simulate contact problem with composite materials which are composed by an isotropic matrix and isotropic fibers. The Eshelby's equivalent inclusion method is considered to quantify the surrounding material stresses. An iterative process is implemented to determine the displacements and stress fields caused by the eigenstrains of all inclusions, accounting for interactions between close inclusions. Each heterogeneity can be made of a single inclusion of simple shape (cuboid, sphere, ellipsoid) or discretized in many inclusions of simple shape. Curved cylindrical fibers are here reproduced by a series of oblate ellipsoidal inclusions. For more details on the method and algorithm the reader may refer to [LER 13]. The proposed method can be seen as an enrichment technique for which the effect of heterogeneous inclusions is superimposed on the homogeneous solution. 3D and 2D Fast Fourier Transforms are utilized to improve the computational efficiency. Displacements and pressure relationships account for this contribution, as well as the coupling between the normal and tangential contact problems. It is found that the presence of heterogeneities in the vicinity of the surface contact affects significantly the contact pressure distribution and subsequently the distribution of shear and slip at the interface.

In the article [LER 10], configurations such as stringer and cluster of spherical inclusions are analysed. The effects of Young's modulus, Poisson's ratio, size and location of the inhomogeneities are also investigated. From a numerical point of view the role of Poisson's ratio is found very important. One of the findings is that a relatively "soft" and nearly incompressible inclusion -for example a cavity filled with a liquid -can be more detrimental for the stress state within the matrix than a very hard inclusion with a classical Poisson's ratio of 0.3.

The stick-slip contact problem was investigated, [LER 11], when at least one of the contacting bodies behaves as an ideal composite material with long fibers perpendicular to the direction of movement. Cylindrical inhomogeneous inclusions within a homogeneous media and with axes parallel to the contact surface are considered.

Fiber volume fraction

The parameter V f determines the ratio of fibers in the matrix, in term of volume. It is an important parameter, it affects directly the effective properties of the composite.

The material is often composed by a cylindrical fiber piled very disorderly. Local properties of a composite are strongly influenced by the distance between fibers. Models based on a perfect arrangement of fibers (square or hexagonal arrangement) are not able to take into account this variability and remain inadequate. Couegnat [COU 08] has implemented a program named GENCELL, to characterize the "disorder" in the arrangement of fibers from the analysis of covariance, in order to determine the minimum size of a representative elementary volume (REV).

Fiber orientation

While the fiber volume fraction is primarily responsible for global mechanical behavior of the composite, the anisotropy directions are due to the orientation of fibers. Experimental measurements are based on the geometrical properties of fibers, which have circular section. Thus, by measuring the length of the axis' elliptical sections of inclined fibers, it is possible to obtain the fiber orientation [GOM 98, MAY 92]. This laborious method suggests that the fibers stay parallel to the midline.

Homogenization method for composite materials

To determine the effective elastic properties of the composite, two analytical approaches are often used for particulate reinforced composites and short fibers: the Mori-Tanaka method and the self-consistent method. Whatever the method used, the effective behaviour of the composite is given by:

C e f f = C m + n ∑ α=1 V α f (C α -C m ) A α (4.1)
where C e f f is the effective stiffness tensor of the composite, A α the strain concentration tensor for an ellipsoidal inclusion α, V α f the volume fraction of α, C α the stiffness tensor of the inclusions and C m the stiffness tensor of the matrix.

The strain concentration tensor takes a different form depending on the homogenization method used. These methods are based on the Eshelby transformation principle. Some equations useful for the developments of different homogenization methods are recalled. The strain field induced by one ellipsoidal heterogeneity Ω α with elastic tensor C α in an infinite elastic matrix defined by C m , subjected to strain at infinity ∞ , is given by:

ε α i j (x) = S α i jkl ε α * i jkl , ∀x ∈ Ω α , (4.2) ε α i j (x) = D α i jkl (x) ε α * i jkl , ∀x / ∈ Ω α . (4.3) α (x)
is the strain due to the inclusion α at the point x and α * is the eigenstrain of the inclusion α. This latter traduces the differences in the material properties between the inclusion and the matrix, it is as a strain of "incompatibility". S α is the Eshelby tensor of the inclusion α. It depends only on the Poisson's coefficient of the matrix and on the geometry of the inclusion. It is expressed in the inclusion's base. This tensor is constant if:

-the matrix is infinite, elastic and isotropic, -the inclusion has an ellipsoidal shape (a sphere, a cylinder, a flat disc, a prolate or an oblate ellipsoid), -the strain field at infinity is uniform. Thus, α and α * are uniform. As the Eshelby tensor exists only for points inside the inclusion, an other tensor, the influence coefficients D, is defined for outside points, depending on the coordinate of this point.

These two equations enable to establish the following equivalence equation:

C α i jkl (ε ∞ kl + S α klmn ε α * mn ) = C m i jkl (ε ∞ kl + S α klmn ε α * mn -ε α * kl ) (4.4) or also α * = [C α (S α -I) -C α S α ] -1 (C α -C m ) ∞ . (4.5)

Mori-Tanaka model

The approach proposed by Mori-Tanaka considers a volume fraction of inclusions V f of same shape and same orientation inside an infinite matrix [MOR 73]. The elastic stiffness tensor C is not necessarily the same for each inclusion. Inclusions can be subdivided into subfamilies such as ∑V α f = V f . The method to estimate the effective moduli of the material consists, for each family of inclusions, in considering an equivalent single ellipsoidal inclusion in the matrix, assumed to be infinite, presupposing that this inclusion is subjected at the infinite to an homogeneous strain field equal to the average strain field in the matrix m . S αα is the Eshelby tensor which represents the effect of the inclusion on itself. S αβ is the Eshelby tensor which represents the effect of the inclusion β on the inclusion α. The total strain field in the area composed by the matrix and by n inclusions is expressed as:

V m f m + n ∑ α=1 V α f α = 0. (4.8)
m represents the average strain field of the matrix, V α f and V m f are the volume fraction of the inclusions and the matrix, respectively.

V m f + n ∑ α=1 V α f = 1 (4.9)
The last term of Eq. 4.6 represents additional disturbances induced by neighbouring inclusions on the inclusion α. It is also called image strain im . In the Mori-Tanaka model, it is considered equal everywhere, in all the inclusions and in the matrix. In other words, all the inclusions are subjected to the same average strain field, although their elasticity tensor can be different. The first step in this process of homogenization is to define a concentration tensor and a dilute concentration tensor (see Fig. 4.2). The concentration tensor connects all strains induced in the inclusion relative to the strain applied at infinity. The dilute concentration tensor connects the total strain in the inclusion relative to the strain of the matrix. From Eq. 4.8, an explicit expression of the concentration tensor is established

A α = A α m V m f I + n ∑ β=1 V β f A β m . (4.15)
The dilute concentration tensor is obtained from the combination of Eqs 4.4 and 4.12 such that

A α m = I + S αα (C m ) -1 (C α -C m ) -1 . (4.16)
Interactions between inclusions are taken into account indirectly through an average of strain and are independent of the considered inclusions. Each inclusion is affected indirectly by the presence of adjacent inclusions through the total strain in the matrix. So, the maximum value of the fiber volume fraction should be less than V f = 0.5, to avoid too many interactions between fibers.

Self-consistent model

This method [HIL 65], originally developed for aggregates (polycrystals), establishes a concentration law by considering each phase of the material as an ellipsoidal inclusion interacting with a fictive matrix, infinite and homogeneous, which has the behaviour of intended effective homogeneous material, characterized by C e f f . The resolution of the equation is performed by an iterative procedure on the effective material tensor C e f f :

-the effective material tensor is initialized by the value of the matrix material tensor C e f f = C m , -the Eshelby tensor S αα is calculated for each individual inclusion. Its value depends on the effective material properties of the matrix, -the strain concentration tensor A α is determined by:

A α = A α m V m f I + n ∑ β=1 V β f A β m -1 (4.17)
where

A α m = I + S αα C e f f -1 C α -C e f f -1 , ( 4.18)
-the tensor of the effective composite material updated is calculated:

C e f f = C e f f + n ∑ α=1 V α f C α -C e f f A α , ( 4.19)
-an iteration over the effective material tensor is done until convergence of the latter.

Convergence is relatively fast, except for situations involving inclusions of material properties very different. The self-consistent approximation takes into account with a global way the interactions between inclusions of an aggregate with the hypothesis of "perfect disorder". In comparison, the Mori-Tanaka method considers them partially.

Results and discussions

To model the contact of a composite material, two methods are possible:

-taking into account explicitly heterogeneities of the material in contact, -use of anisotropic material properties obtained by homogenization technique. The first criterion to choose between both methods is the characteristic dimension of the heterogeneities with respect to the dimensions of the contact area.

To compare these two models in the case where heterogeneities have a size which can not be neglected compared to the contact scale, a frictionless normal indentation is performed on an isotropic elastic half space, composed with identical isotropic unidirectional cylindrical fibers (see Fig. 4.3). A normal load F n = 10kN is applied to a rigid sphere of radius R = 40mm. The Young's modulus and Poisson's ratio of the half space which is the matrix of the composite material are taken equal to E m = 120GPa and ν m = 0.3. For an isotropic homogeneous material without inclusions or fibers (E = E m and ν = ν m ), the Hertz theory gives a contact radius a Hertz = 1.315mm and a maximal contact pressure P Hertz = 2760MPa. In what follows each center of cylindrical fiber will be at the same position. The diameter of fibers changes with the fiber volume fraction. The first layer of fibers is located at the depth z = 0.120mm (position of the center). Note that the dimension and location of the fibers does not affect the global homogenized behaviour in an infinite space. However, they do affect the contact pressure distribution and subsurface stress field in the heterogeneous model.

The results for anisotropic homogeneous materials are obtained with the influence coefficients for coated anisotropic half space. Coating and substrate have the same properties and the thickness of the coating is nil. Two parametric studies are carried out: one on the ratio γ of the Young's modulus of the fiber and the matrix, and another one on the volume fraction of fibers V f . The Poisson's ratio for the fibers and the matrix are equal ν f = 0.3. Concerning the geometry features of transverse isotropic behaviour tends to minimize the local contact pressure. The contact area obtained using the homogeneous model is greater than that obtained with the heterogeneous model. It is further no more circular as obtained for homogeneous isotropic half space and becomes slightly shrinked in the direction where the Young's modulus is higher (always along x or L here). In the case γ = 0.5, the homogeneous model gives an average value of the fluctuation pressure produced by heterogeneities. Contact areas obtained by the two approaches are very similar. This is also explained by the close values of E L and E T , and the Poisson's ratios ν LT and ν T T . For both approaches, the contact area is enlarged compared to that obtained in the Hertz case, which represents a softer material than the Hertz reference case.

Conversely, in the case where fibers are more rigid than the matrix (γ > 1), the homogeneous transverse isotropic behaviour gives the upper limit of the contact pressure distribution which surrounded pressure peaks obtained by the heterogeneous approach. The contact area obtained with the homogeneous model is markedly reduced in all directions, again with the contact ellipse slightly elongated in the transverse direction (direction y or T ). These differences can be explained by the numerical value of the longitudinal modulus E L much higher than the transverse module E T , itself higher than the one of the matrix, and a Poisson's ratio ν T T greater than ν LT for ratios γ = 10 and γ = 20. 4.7 shows the influence of the fiber volume fraction on the distribution of contact pressure in the plane (Oyz) as well as on the contact area. For both models, the contact area along y is reduced when the fiber volume fraction increases. However, only the homogeneous model tends to have a reduced length along the fiber direction x, which is more remarkable for high volume fraction. Moreover, increasing the fiber volume fraction raises the amplitude of pressure peaks. For small fiber volume fractions, the contact pressure is overestimated by the approach of homogenization.

n ˚γ V f E L [MPa] E T [MPa] G LT [MPa] ν LT ν T T 2 -1 10 0.
A comparison is now carried out for the stresses. A normal load F n = 20kN, more important than previously, is applied to the spherical rigid indenter. The material properties of the matrix are identical to the one considers in the previous case studies (E m = 120GPa and ν m = 0.3). In the case of homogeneous material, the Hertz theory prescribes a contact radius a Hertz = 1.657mm and a maximum contact pressure P Hertz = 3478MPa. Fibers are twice harder than the matrix (γ = 2) with a Poisson's ratio ν f = 0.3 and the fiber volume fraction is V f = 0.4. Figure 4.8 shows the distribution of normal stresses (σ xx , σ yy , σ zz ) along the depth for both models. Section planes are: the plane x = 0, the plane y = 0 and the plane z = 0.22a Hertz which cuts across the center of the fibers constituting the second layer parallel to the contact surface. The homogenized anisotropic material tends 4. Comparison with an isotropic heterogeneous model to be more constrained than the heterogeneous isotropic material. As discussed earlier we encountered some difficulties in the calculation of some stress components for the layer close to the surface. That corresponds to the thick red lines in the left plots of Fig. 4.8.

With the observations above, it may be interesting to analyse these differences with a macroscopic point of view by looking at the indentation curve for configuration n ˚1-5 (γ = 10). Figure 4.9 shows the variation of the normal force normalized by E m R 2 indenter relative to the displacement response normalized by R indenter . When the effort increases, an important difference can be observed between both models. This is coherent with those previously done for the micro scale as the contact area is smaller for the homogeneous approach than for the heterogeneous approach. Conversely, and surprisingly, for the heterogeneous model the contact area becomes very elongated along the longitudinal direction (x or L) along which fibers are lying and for which the Young's modulus is the higher. It means that the geometrical features of the composite reinforcement dominates the behaviour, compared to the anisotropic elastic properties. However, this large difference between the two curves may come from the homogenisation method. Indeed, the Mori-Tanaka method assumes an infinite space composed by a matrix and heterogeneities which is reduced to a REV with a uniform or random heterogeneities distribution. To simulate the contact between two bodies, the condition of semi-infinite space imposes a restrictive framework, particularly in homogenization approaches where boundary conditions are obsolete. The Mori-Tanaka method has therefore limitations to its application in this kind of study. Further reflection and future developments will be required to introduce the concept of half-space in the homogenization method.

Conclusion

The anisotropic homogeneous model presented in this thesis is compared to an isotropic heterogeneous model which takes into account the inhomogeneities into the half space. Two methods are described to obtain the homogenized material properties from the properties of the heterogeneous materials. For the comparison, the Mori-Tanaka model is used. The volume fraction of fibers is taken inferior to 0.5. The influence of the Young's modulus ratio between the matrix and fibers has been investigated on the contact pressure and the contact area. When fibers are softer than the matrix, the homogeneous model tends to minimize the contact pressure and maximize the contact area. Inversely, when the fibers are harder, the homogeneous model tends to maximize the contact pressure. The homogeneous model gives smooth profile pressure whereas with the inhomogeneous model, the fibers can be visualized, the profile pressure presents some peaks and grooves. Interestingly it has been noticed that, while the contact area becomes elliptical with the semi-minor axis lying along the direction where the Young's modulus is the larger (direction x or L here) with the anisotropic homogenized model, the shape of the contact area with the isotropic heterogeneous model is mostly driven by the orientation of unidirectional fibers. The fiber volume fraction affects also the homogenized material properties. An increase on the fiber volume fraction raises up the amplitude of the pressure peaks.

It has been observed that the results obtained using these two models can be notably different, both at the mesoscopic scale (i.e. at the contact scale) and at the macroscale (i.e. for the indentation curve). This can be explained by the limitation of the Mori-Tanaka model that is valid for a representative element volume with symmetrical boundary conditions, which does not hold anymore near the contact surface (or close to any free surface).

Regarding the subsurface stress pattern, it is obvious that the isotropic heterogeneous model is the only one than could accurately describes the stress distribution within the material.

Conclusion and Perspectives

The semi analytical contact solver developed at the LaMCoS laboratory has been improved by introducing the anisotropy of the material, for uncoated and coated surfaces.

The influence coefficients linking the elastic displacements to the surface load have been formulated using the anisotropic Green's functions. Starting from three dimensional Green's functions established for anisotropic infinite bimaterials and trimaterials, these functions were rewritten for anisotropic half space and coated half space, respectively. The Fourier transforms are used and solutions are expressed in the complex domain. According to these new improvements, it is possible to address the normal contact of anisotropic half space, with or without an anisotropic coating. Contact solutions obtained for indentation conditions have been compared to those obtained with a finite element code (Abaqus version 6.9) and a very good agreement was found, which validates the influence coefficient formalism. Regarding the corresponding internal stresses, some components need more investigations since they present some divergence in the near surface. The source of this error remains to be identified.

Parametric studies have been performed on the effect of the Young's modulus along different directions, on the Poisson's ratio, on the Coulomb's modulus, on the material orientation relative to the contact normal, both for the substrate and the coating, and on the coating thickness. The pressure distribution, the contact area and the stresses have been analysed. It is observed that the Young's modulus in the direction normal to the contact has a strong influence on the contact parameters (pressure distribution and subsequently subsurface stresses, size of the contact area). In addition a difference in the Young's modulus along the directions parallel to the surface results in a change of the contact shape, which becomes elongated in the direction where the Young's modulus is higher and shrunk in the perpendicular direction (i.e. for a spherical indentation the contact area is no more circular but becomes elliptical).

The manuscript ends with a comparison between two models, the anisotropic homogeneous model detailed here and an isotropic heterogeneous model developed in parallel by J. Leroux in his PhD thesis. An application to an idealized composite material is carried out, i.e. a matrix reinforced by unidirectional and cylindrical fibers parallel to the surface. The homogenized anisotropic properties have been derived through the Mori-Tanaka model. According to the scale of the contact, composite materials can be seen as heterogeneous materials (which create fluctuations in the pressure and in stresses) or as anisotropic homogeneous materials (the pressure profile and stresses are smooth). It has been observed that the results obtained using these two models can be notably different, Conclusion and Perspectives both at the mesoscopic scale (i.e. at the contact scale) and at the macroscale (i.e. for the indentation curve). This can be explained by the limitation of the Mori-Tanaka model that is valid for a representative element volume with symmetrical boundary conditions, which does not hold anymore near the contact surface (or close to any free surface). This work will be extended in order to deal with the effect of the temperature. In several engineering applications, parts in contact are subjected to thermal loading or exposed at high temperature, which may affect material properties, and in case of coating/substrate a mismatch of thermal properties may induce localized strains and stresses. Future developments would be to account for plasticity for anisotropic material. Simulating the wear due to fretting would be also an interesting point to look at.

1. 1

 1 Blade of high pressure compressor . . . . . . . . . . . . . . . . . . . . . 1.2 Blade root of CoHP subjected to fretting . . . . . . . . . . . . . . . . . . 1.3 Contact zone of the disk in the TuHP subjected to fretting . . . . . . . . . 1.4 Outer ring ball damaged by fretting . . . . . . . . . . . . . . . . . . . . . 1.5 Pressure distribution discretized in cells of uniform pressure . . . . . . . 1.6 General algorithm of the contact solver by the semi analytical method . . 1.7 Euler angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Parameters used for an layered half space . . . . . . . . . . . . . . . . . 2.1 Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Validation of the pressure profile for an isotropic case . . . . . . . . . . . 2.3 Finite element model with a detailed view of the contact area . . . . . . . 2.4 Validation of the pressure profile for an orthotropic material . . . . . . . . 2.5 FEM and SAM pressure profiles, for E 1 = E 2 = 2E 3 . . . . . . . . . . . . 2.6 Comparison of the pressure profile for layered isotropic medium . . . . . 2.7 Validation of the pressure profile for layered orthotropic material . . . . . 3.1 Coordinate and material axes . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Influence of E 1 on the contact pressure, in the plane x = 0 . . . . . . . . . 3.3 Influence of E 1 on the contact area . . . . . . . . . . . . . . . . . . . . . 3.4 Influence of E 3 on the contact pressure, in the plane x = 0 . . . . . . . . . 3.5 Influence of E 1 on the maximum pressure . . . . . . . . . . . . . . . . . 3.6 Influence of E 3 on the maximum pressure . . . . . . . . . . . . . . . . . 3.7 Influence of E 1 on the indentation curve . . . . . . . . . . . . . . . . . . 3.8 Influence of E 3 on the indentation curve . . . . . . . . . . . . . . . . . . 3.9 Influence of the material's orientation on the contact pressure . . . . . . . 3.10 Influence of the material's orientation on the maximum pressure . . . . . 3.11 Influence of the Poisson's ratio on the contact area . . . . . . . . . . . . . 3.12 Influence of the Coulomb's modulus on the contact area . . . . . . . . . . 3.13 Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Influence of E C1 on the contact pressure . . . . . . . . . . . . . . . . . . 3.15 Influence of E C1 on the contact area shape . . . . . . . . . . . . . . . . . 3.16 Influence of E C3 on the contact pressure . . . . . . . . . . . . . . . . . . iii List of figures 3.17 Influence of the material's orientation on the maximum pressure . . . . . 3.18 Influence of Z C on the contact pressure . . . . . . . . . . . . . . . . . . . 3.19 Influence of Z C on the maximum pressure . . . . . . . . . . . . . . . . . 3.20 Influence of E S1 on the contact pressure . . . . . . . . . . . . . . . . . . 3.21 Influence of E S3 on the contact pressure . . . . . . . . . . . . . . . . . . 3.22 Influence of Z C on the contact pressure . . . . . . . . . . . . . . . . . . . 3.23 Stress tensor components versus depth for an isotropic material . . . . . . 3.24 Stress tensor components for an isotropic material, in the plane x = y = 0 3.25 Stress tensor components for an orthotropic material, E 3 = 2E 2 . . . . . . 3.26 Stress tensor components versus depth for an isotropic coated surface . . 3.27 Stress tensor components for an orthotropic material, E 1 = 2E 2 . . . . . . 3.28 Stress tensor components for an orthotropic material, E 3 = 2E 2 . . . . . . 3.29 Influence of E 1 on the stress profile along depth . . . . . . . . . . . . . . 3.30 Influence of E 3 on the stress profile along depth . . . . . . . . . . . . . . 3.31 Influence of θ m on the stress components σ xx and σ yz . . . . . . . . . . . 3.32 Influence of θ m on the stress components σ yy and σ zz . . . . . . . . . . . 4.1 Equivalence of material properties . . . . . . . . . . . . . . . . . . . . . 4.2 The strain concentration tensors . . . . . . . . . . . . . . . . . . . . . . 4.3 Normal indentation on materials composed with cylindrical fibers . . . . 4.4 Coordinate systems for a transverse isotropic material . . . . . . . . . . . 4.5 Influence of the Young's modulus ratio γ < 1 . . . . . . . . . . . . . . . 4.6 Influence of the Young's modulus ratio γ > 1 . . . . . . . . . . . . . . . 4.7 Influence of the fiber volume fraction V f . . . . . . . . . . . . . . . . . . 4.8 Normal stresses in the section planes x = 0, y = 0 and z = 0.22a Hertz . . . 4.9 Indentation curve of both models . . . . . . . . . . . . . . . . . . . . . . iv

Figure 1 . 1 :

 11 Figure 1.1: Blade of high pressure compressor

Figure 1 . 2 :

 12 Figure 1.2: Blade root of CoHP from an engine high power subjected to high temperature fretting

Figure 1 . 3 :

 13 Figure 1.3: Contact zone of the disk in the TuHP subjected to fretting in temperature

Figure 1 . 4 :

 14 Figure 1.4: Outer ring ball damaged by fretting in temperature

  Coulomb's law of friction is used in most existing analytical solutions. When contact is considered fully sliding, shears are directly obtained from the hertzian pressure. Stresses for a cylindrical contact are given by McEwen [MCE 49] while stresses for a spherical contact are given by Hamilton [HAM 63]. Those solutions are extended to elliptical contacts by Sackfields & Hills [SAC 83]. Cattaneo [CAT 38]

Figure 1 . 5 :

 15 Figure 1.5: Pressure distribution discretized in cells of uniform pressure

Figure 1 . 7 :

 17 Figure 1.7: Euler angles, e 1 e 2 e 3 the material system and a 1 a 2 a 3 the contact system

Figure 1 . 8 :

 18 Figure 1.8: Parameters used for a spheric indenter in contact with an layered half space

  Influence coefficients for anisotropic half space . . . . . . . . . . . . . . 26 2.3.1 Green's functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3.2 Cauchy principal value . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.3 Hadamard finite part . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4 Influence coefficients for coated anisotropic half space . . . . . . . . . . 32 2.4.1 General formulation . . . . . . . . . . . . . . . . . . . . . . . . . 32

(

  exists and is finite. In this case, the limit L is called the Cauchy principal value of the improper integral, what is writtenv.p. b a f (x)dx = L.The definition extends easily to the case with n singularities a < x 1 , ..., x n < b. If for ε > 0 the integrals x 1 -ε a f (x)dx and b x n +ε f (x)dx exist and are finite, and the limit lim ε→0 -θ 0 ) dθ can be calculated by avoiding the singularities equal to θ 0 + π/2 modulo n. The difficulty is to determine the value of ε.

Figure 2 . 1 :

 21 Figure 2.1: Coordinate systems

Figure 2 . 2 :

 22 Figure 2.2: Validation of the pressure profile for an isotropic case by comparison with the Hertz solution

Figure 2 . 3 :

 23 Figure 2.3: Finite element model with a detailed view of the contact area

Figure 2 . 4 :

 24 Figure 2.4: Validation of the pressure profile for an orthotropic material. Comparison between FE and SAM results.

Figure 2 . 5 :

 25 Figure 2.5: FEM and SAM pressure profiles. Comparison for E 1 = E 2 = 2E 3 .

Figure 2 . 6 :

 26 Figure 2.6: Comparison of the pressure profile for layered isotropic medium. Comparison between O'Sullivan 's solutions and SAM ones.

Figure 2 . 7 :

 27 Figure 2.7: Validation of the pressure profile for orthotropic material when the coating is in the same material than the substrate by comparison with FE model, E 3 = 2E 2

Contents 3 . 1

 31 Parametric studies for anisotropic half space . . . . . . . . . . . . . . . 48 3.1.1 The Young's modulus . . . . . . . . . . . . . . . . . . . . . . . . 48 3.1.2 The material's orientation . . . . . . . . . . . . . . . . . . . . . . 50 3.1.3 The Poisson's ratio and the Coulomb's modulus . . . . . . . . . . . 51 3.2 Parametric studies for layered anisotropic half space . . . . . . . . . . . 55 3.2.1 An anisotropic coating on an isotropic substrate . . . . . . . . . . . 55 3.2.2 An isotropic coating on an anisotropic substrate . . . . . . . . . . . 58 3.3 Stress field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.3.1 Coefficient for anisotropic half space . . . . . . . . . . . . . . . . 62 3.3.2 Coefficient for coated anisotropic half space . . . . . . . . . . . . . 66 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

Figure 3 . 1 :

 31 Figure 3.1: Coordinate and material axes

Figure 3 . 2 :-

 32 Figure 3.2: Influence of E 1 on the contact pressure, in the plane x = 0

Figure 3 . 3 :

 33 Figure 3.3: Influence of E 1 on the contact area

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Influence of E 3 on the contact pressure, in the plane x = 0

Figure 3 . 6 :

 36 Figure 3.6: Influence of E 3 on the maximum pressure

Figure 3 . 7 :

 37 Figure 3.7: Influence of the Young's modulus along axis 1 (parallel to the surface) on the indentation curve

Figure 3 . 8 :

 38 Figure 3.8: Influence of the Young's modulus along axis 3 (normal to the surface) on the indentation curve

Figure 3 . 10 :Figure 3 . 11 :Figure 3 . 12 :

 310311312 Figure 3.10: Influence of the material's orientation on the maximum pressure

Figure 3 . 14 :Figure 3 . 15 :

 314315 Figure 3.14: Influence of E C1 on the contact pressure (anisotropic coating over isotropic substrate, Z C = 0.5a Hertz )

Figure 3 . 16 :

 316 Figure 3.16: Influence of E C3 on the contact pressure (anisotropic coating over isotropic substrate, Z C = 0.5a Hertz )

Figure 3 . 17 :

 317 Figure3.17: Influence of the material's direction on the maximum pressure (anisotropic coating over isotropic substrate, Z C = 0.5a Hertz )

Figure 3 . 18 :

 318 Figure 3.18: Influence of Z C on the contact pressure (anisotropic coating over isotropic substrate, E C3 = 2E S )

Figure 3 . 19 :

 319 Figure 3.19: Influence of Z C on the maximum pressure (isotropic substrate)

Figure 3 . 20 :Figure 3 . 21 :

 320321 Figure 3.20: Influence of E S1 on the contact pressure (isotropic coating over an anisotropic substrate, Z C = 0.5a Hertz )

Figure 3 . 22 :

 322 Figure 3.22: Influence of Z C on the contact pressure (isotropic coating over an anisotropic substrate, E S3 = 2E C )

Figure 3 . 23 :Figure 3 . 24 :Figure 3 . 25 :

 323324325 Figure 3.23: Stress tensor components versus depth (z) at the center of the contact (x = y = 0) for an isotropic material

Figure 3 . 26 :Figure 3 . 27 :Figure 3 . 28 :

 326327328 Figure 3.26: Stress tensor components versus depth (z) at the center of the contact (x = y = 0) for a coated surface. The coating and the substrate behave isotropically and have the same elastic properties. Results for two coating thicknesses, Z C = 0 and Z C = 0.5a Hertz , are compared to the Hertz solution.

Figure 3 . 29 :

 329 Figure 3.29: Influence of E 1 (= /E 3 , E 2 = E 3 ) on the stress profile along depth at the contact center (x = y = 0)

Figure 3 . 30 :Figure 3 . 31 :Figure 3 . 32 :

 330331332 Figure 3.30: Influence of E 3 (E 1 = E 2 ) on the stress profile along depth at the contact center (x = y = 0)

Figure 4 . 1 :

 41 Figure 4.1: Equivalence between a material composed with heterogeneities and a material composed with inclusions by adding eigenstrains We seek to determine the average strain field in an inclusion α, surrounded by n -1 inclusions (see Fig. 4.1).

mΣFigure 4 . 2 :

 42 Figure 4.2: The strain concentration tensors

  with an isotropic heterogeneous model

Figure 4 . 3 :

 43 Figure 4.3: Normal indentation on an elastic half space composed with unidirectional cylindrical fibers, evenly distributed

0 - 0 -

 00 n ˚1-1 (γ = 1/20): Contact pressure in x = n ˚1-2 (γ = 1/10): Contact pressure in x = n ˚1-3 (γ = 1/2): Contact area

Figure 4 . 5 : 0 -

 450 Figure 4.5: Comparison between anisotropic homogeneous and isotropic heterogeneous materials: influence of the Young's modulus ratio γ < 1, with V f = 0.4.

Figure 4 . 6 : 0 - 0 - 0 -

 46000 Figure 4.6: Comparison between anisotropic homogeneous and isotropic heterogeneous materials: influence of the Young's modulus ratio γ > 1, with V f = 0.4.

Figure 4 . 7 :

 47 Figure 4.7: Comparison between anisotropic homogeneous and isotropic heterogeneous materials: influence of the fiber volume fraction V f , with γ = 10.

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Comparison between anisotropic homogeneous and isotropic heterogeneous materials: normal stresses in the section planes x = 0, y = 0 and z = 0.22a Hertz (the center of the second layer of unidirectional fibers).

  . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Coefficients for ũ * 01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Coefficients for ũ * 02 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Coefficients for ũ * 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 The influence of the mesh and the model on the CPU time . . . . . . . . 2.6 Number and type of elements in the FE model . . . . . . . . . . . . . . . 2.7 Contact problem CPU time . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Elastic properties of studied materials . . . . . . . . . . . . . . . . . . . 3.2 Summary of stresses' validation from anisotropic model . . . . . . . . . 3.3 Summary of stresses' validation from coated anisotropic model . . . . . . 4.1 Material parameters for the influence study of the Young's modulus ratio . 4.2 Material parameters for the influence study of fiber volume fraction . . .

	Nomenclature
	(x, y, z) Contact reference system
	(1, 2, 3) Material reference system
	δ(x)	Dirac delta function
	γ ν φ σ	Nomenclature Young's modulus ratio between fibers and matrix Strain tensor Poisson's ratio Orientation of the vector position according to the depth axis direction Stress tensor
	A B C C i jkl σ comp Stress tensor relative to the complementary part Matrix of the eigenvectors Matrix derived from the eigenvectors Matrix derived from the eigenvectors Elastic stiffness tensor σ ∞ Stress tensor relative to the infinite part θ Orientation of the vector position according to the surface axis direction θ m
	E	Young's modulus
	E i	Young's modulus in direction i
	E C	Young's modulus of the coating
	E Ci	Young's modulus of the coating in direction i
	E S	Young's modulus of the substrate
	E Si	Young's modulus of the substrate in direction i
	G	Coulomb's modulus
	I	Identity matrix
	P	Contact pressure distribution
	P Hertz Maximum Hertzian pressure
	Q	Double projections of the elastic stiffness tensor
	R	Double projections of the elastic stiffness tensor
	T	Double projections of the elastic stiffness tensor
	V f	Fiber volume fraction
	Z C	Coating thickness
	a Hertz Hertzian contact radius
	a i	Eigenvector associated to p i
	b i	Vector derived from a i
	c i	Vector derived from a i
	f	Concentrated force
	i	Imaginary unit
	p	Vector of eigenvalues
	p i	Eigenvalue
	r	Magnitude of the vector position
	s	In-plane stress vector
	t	Traction vector
	u	Displacement vector
	u comp Displacement vector relative to the complementary part
	u ∞	Displacement vector relative to the infinite part
	x i	Cartesian coordinates of the studied point
	viii	v

  [MUS 53] , then by Mikhlin [MIK 57], Galin [GAL 53] and Aleksandrov [ALE 86]. « Integral transforms » are also used, such as the Fourier transform [SNE 51]. Westergaard [WES 39] used it in rough contacts if roughness is a sine function. Greenwood & Willamson [GRE

.3 Materials in contact 1.3.1 Anisotropic materials

	1. Contact mechanics 1. Contact mechanics
	tangential effects, which is required for frictional contact problems between elastically dissimilar materials. The same group of researchers also improved the contact algorithm in several ways so that it becomes affordable to simulate wear for 3D contact problems, cycle after cycle [NÉL 06, GAL 06, GAL 07b]. Note that SAM could be also linked with 1Motivations Most substrates are assumed isotropic as a first approach. However for the FE method and used as a zoom on the contact to account for the effect of rough-most composite and monocrystal materials their compositions or the elaboration and man-ness or simulate wear [GAL 10a]. More recently Chaise and Nelias [CHA 11a] improved ufacturing processes imply that it exists one or two main directions or even a general the numerical model to account for kinematic and isotropic hardening and analyzed the anisotropy. In aerospace applications hot parts such as turbine blades have special crys-problem of a rolling load versus indentation. They also proposed a method to predict the coefficient of restitution when an elastic or rigid sphere is impacting an elastic-plastic half-space [CHA 11b]. This PhD is based on Jacq, Gallego, Boucly and Fulleringer semi analytical codes Initial state • Loading path Normal Contact Problem (CGM) Current displacement δ z (n) Or Load F c (n) Loading increment n Tangential Contact Problem (CGM) • Loading path Problem (CGM) Or Load F c (n) Problem (CGM) Initial state Normal Contact Current displacement δ z (n) tallographic directions and so behave as anisotropic materials. The material properties Loading increment n Tangential Contact should be considered accurately to optimize the service life of industrial components.
	developed at the LaMCoS, in the programming language Fortran. Structures of contact algorithms and numerical methods remain unchanged and are not detailed in this docu-ment (see [GAL 07a] and [FUL 11]). However, the new influence coefficients used for anisotropic materials, with or without coating, are detailed. Displacement δ z / Load F c history • Geometry • Hardening law • Initial plastic strains • Initial Stresses Elastic Stresses computation Pressure components p, q x , q y σ e Next iteration k=k+1 Displacement δ z / Load F c history • Geometry Characterization In this thesis, materials are assumed homogeneous. The validity of Pressure components p, q x , q y this assumption depends on the scale at which the model is used. Although a composite • Hardening law • Initial plastic strains • Initial Stresses Elastic Stresses computation σ e is by definition heterogeneous, it is generally possible to use a homogeneous model if the Next iteration k=k+1 scale of interest is much larger than that of heterogeneities [RÖS 07].
				σ tot = σ e + σ res + σ init σ tot = σ e + σ res + σ init
			Return-mapping algorithm Return-mapping algorithm
		σ res σ res	• Plastic strain increment δε p computation • Plastic strain increment δε p computation
			• ε p(k) = ε p(k) + λ.δε p(k-1) • ε p(k) = ε p(k) + λ.δε p(k-1)
				• Residual displacement • Residual displacement
	Residual Stresses computation Residual Stresses computation	increment δu res increment δu res
	FFT2D -FFT3D -Combination FFT2D -FFT3D -Combination	• Yield function f(x,y,z) • Yield function f(x,y,z)
			Convergence Convergence
	Update geometry • u res(k) = u res (k) + δu res (k-1) • u res(k) = u res (k) + δu res (k-1) Update geometry	No No	• Contact problem |δu res(k) | < tol. • Plastic problem • Contact problem • Plastic problem |δu res(k) | < tol.
			f(x,y,z)<tol. f(x,y,z)<tol.
				Yes Yes
	The semi analytical contact solver is based on the No Next loading increment n=n+1 End of loading history Next loading increment n=n+1 No End of loading history
	Yes Yes pioneering work of Jacq et al. [JAC 02] for elastic-plastic contacts. The solver has since Common stages Common stages Common stages Optional stages Optional stages Optional stages
	been developed and improved in several ways. Boucly et al. [BOU 05] and Chen et al. End End
	[CHE 08] introduced thermal aspects. Fulleringer and Nelias [FUL 10] focused on the
	influence of a cuboid of uniform plastic strain in a half space, on the normal and tangen-tial displacements of surface points, and derived the corresponding Green's functions in Figure 1.6: General algorithm of the contact solver by the semi analytical method
	an analytical form. Then Leroux et al. [LER 10] and Zhou et al. [ZHO 11] studied the ef-
	fects of the presence of inhomogeneous inclusions within a half-space, both on the contact
	pressure distribution and on the corresponding subsurface stress field. Discrete Convolu-
	tion (DC) and 3D FFT have been since then widely used in the contact solvers. Leroux
	and Nelias [LER 11] worked also on the stick-slip problem for a sphere in contact with a
	flat half-space containing unidirectional cylindrical fibers. Gallego et al. [GAL 10b] pro-
	posed an algorithm based on the CGM to account for the coupling between normal and

Table 2 .

 2

7: Contact problem CPU time, with and without the influence coefficient computation

Table 3 .

 3 2 summarises the components of the stress which need further investigations. Anisotropic model is not in a good agreement for Isotropic material σ yz for z < 0.8a Hertz σ zz for z < 0.2a Hertz Orthotropic material σ xx , σ yy , σ zz , σ yz

	Table 3.2: Summary of stresses' validation from anisotropic model

Table 3 .

 3 3 summarises the components of the stress which are not validated and which need further investigations.

	Coated anisotropic model is not in a good agreement for
	Isotropic material Orthotropic material	σ zz for z < 0.2a Hertz σ zz for z < 0.2a Hertz
		(Better results with Z C = 0)
	Coated orthotropic material	invalidated
	Table 3.3: Summary of stresses' validation from coated anisotropic model

Table 4 .

 4

	1 227111	142047	108767	0.300 0.344
	2 -2 10 0.2 334417	165320	128447	0.300 0.357
	2 -3 10 0.4 549622	226825	182164	0.300 0.361

2: Material parameters for the influence study of fiber volume fraction Figure

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2013ISAL0014/these.pdf © [C. Bagault], [2013], INSA de Lyon, tous droits réservés

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2013ISAL0014/these.pdf © [C. Bagault], [2013], INSA de Lyon, tous droits réservés 2. A semi analytical method for contact between anisotropic materials

Acknowledgment-Remerciements Acknowledgement

This work was supported by the French FUI Project Innolub. The author would like also to gratefully acknowledge the Rhone Alpes Region for providing a doctoral mobility grant.

N = 1 G

(1) ji = (J jk 2 )

(1)

(1)

1 + (r k 1 )

(1)

(2)

1 + (r k 1 )

(2)

(3)

1 + (r k 1 )

(3)

ji and s (N)

The thickness of the coating is taken equal to zero, so it corresponds to a homogeneous media.

In Figs. 3.29 and 3.30, the influence of E 1 and E 3 is studied on each stress components: σ xx , σ yy , σ zz , σ xy , σ xz , σ yz . If i is the direction studied, the ratio between E i and E 2 is 0.25, 0.5, 1, 1.5, 2, 4. The stress components are normalized with the Hertz pressure, corresponding to an isotropic half space with Young's modulus E = E 2 . Results are plotted versus the depth normalized by the Hertz's radius at x = y = 0. The scale is the same for all figures except for the stress σ xz , in order to see the difference between the configurations.

The Young's modulus E 1 affects particularly the stress σ xx (the one in direction 1). The maximal value is multiplied by 2.8 between E 1 = 0.25E 3 and E 1 = 4E 3 , whereas it is multiplied by 2.1 between E 3 = 4E 1 and E 3 = 0.25E 1 . The influence of E 1 is more important than E 3 influence on σ xx . Moreover, increasing E 1 causes an increase in σ xx too whereas increasing E 3 causes a drop in σ xx .

The effect of the Young's modulus E 3 is shown in Fig. 3.30. The principal change occurs for the σ zz component, which increases by 29% when E 3 is multiplied by two at the depth z = 0.3a Hertz . Note that as discussed earlier the numerical values for σ zz are incorrect between the surface and the depth z = 0.3a Hertz .

In Figs. 3.31 and 3.32, we analyse the influence of the material orientation relative to the contact, with θ m = 0 ˚, 30 ˚, 60 ˚, 90 ˚. The Young's modulus in direction 3 (the depth when θ m = 0 ˚and the surface when θ m = 90 ˚) is twice the Young's modulus in direction 1 and 2. Stresses are presented in the plane x = 0 and they are normalized by the Hertz's pressure.

The study of the influence of the material orientation is similar to that of the influence of Young's modulus. Indeed, σ xx is higher when the Young's modulus along the surface is more important. It is the same for σ yy . Similar results are also obtained for σ zz . Over the angle increases (Young's modulus along the depth decreases), the maximum value is lower and decreasing stress in depth is faster.

Chapter 4 Comparison with an isotropic heterogeneous model

In this last part, the anisotropic homogeneous semi analytical model is compared with an isotropic heterogeneous model, which describes the fibers in composite materials. 4. Comparison with an isotropic heterogeneous model the fibers, ordered arrangement, circular geometry in the plane (Oyz) and assumed infinite length in the x direction (the contact dimension) allow us to use a square REV. The Mori-Tanaka algorithm presented above provides effective homogenized material properties at the macro scale. The maximum value of the fiber volume fraction is less than V f = 0.5. This value defines the validity domain of the Mori-Tanaka method. Beyond this limit, the interactions between fibers are too high and make this method inadequate. The solution is to use the self-consistent algorithm but it can cause convergence problems for low ratios of Young's moduli.

x, 1, L z, 3, T y, 2, T The homogenization of this configuration generates a transverse isotropic material due to the symmetry properties. The elastic mechanical properties of the material are determined by five independent elastic constants: Recall that it has been shown previously that E 3 = E T has a great influence on the contact parameters (pressure and contact area). γ is a combination of E L and E T influences.

For each configuration, contact pressure and contact area are analysed and compared for the two models, the first one considering an equivalent anisotropic homogeneous material and the second one considering explicitly the heterogeneities in the material. Figures 4.5 and 4.6 show the influence of the ratio between Young's modulus on the contact pressure in the plane (Oyz) and the evolution of the contact area for either fibers more compliant or rigid than the matrix.

In the case where fibers are more compliant than the matrix (γ < 1), the homogeneous