Influence of contact parameters in fretting-fatigue contact 3D problems

A numerical study on the influence of contact parameters on fretting fatigue behavior is carried out by using a new method SFEM (Stretching Finite Element Method) is presented in this article. Several parameters are made to vary: geometric shapes of the mesh, materials and contact parameters. The three-dimensional parametric model is composed by specimen and a pad in full contact. A Fortran Code is used to generate the parametric mesh. The stress intensity factors are calculated by varying the above contact parameters and the stress intensity factor under modes I, II and III are computed..


INTRODUCTION
retting occurs when two bodies are in contact, where a load presses these two bodies together [1]. A damage process is created between them. Fatigue is a phenomenon that occurs when a metal component is subjected to variable loads over time and after a given number of loading, a crack can initiate and grow until the component reaches breaking point [2]. Fatigue can occur when there are parts in contact, called fretting fatigue [2]. This phenomenon may occur in several mechanical connections such as the turbines, bearing; steel cables etc. Fretting fatigue can reduce the lifetime of a component compared with plane fatigue [4]. Many researchers use the fretting fatigue models to understand the details of fretting fatigue and to classify the importance of the design parameters involved in this phenomenon [5]. In recent years, a major effort has been made to study the phenomenon of fretting fatigue. Nowell et al. [6] describes a series of fretting fatigue experiments carried out under closely controlled conditions of partial slip. Ochi et al. [7] investigated the effect of contact pressure and slip amplitude on fretting fatigue behaviour of rail steel with fishplate: the results shows that the fretting fatigue strength decrease with increasing the contact pressure and the pad length. Tesh et al. [8] investigated the titanium alloy Ti-6Al-4V behaviour in fretting fatigue using a sphere on flat contact, to study the influence of microstructure and F contact condition. Kond et al. [9] proposed a local stress concept to evaluate the fretting fatigue limit for contact edge cracks. Ciavarella et al. [10] solved the contact problem compared with plain of fatigue for a flat punch having rounded corners under condition of normal loading and a shearing force. Heung et al. [11] proposed hybrid layer method to examine the effect of finite width of contact finite on fretting fatigue according to this method traction distribution obtained from the coarse finite element three-dimensional analyses was used in two-dimensional plane strain finite element model to increase computational efficiency. Kataoka et al. [12] study the effect of contact conditions on growth of small crack in fretting fatigue. Juoksukangas et al. [13] compared the effect of both round contact edge and the sharp contact edge on fretting fatigue behavior of a complete contact using the finely multi-axial criterion and the theory of criterion distance to calculate the cracking. In another work [14] they compared between the measured displacement from digital image correlation and the computed displacements of corresponding finite elements method, where the displacement field and its derivate were compared using the friction coefficient as a variable in the numerical model. Eugenio et al. [23] use the finite element method (FEM ), the extended finite element method (XFEM) And a simple experimental to predicted crack propagation direction and path on fretting fatigue in complete contact condition. Noraphaiphipaksa et al. [15] used the maximum shear stress range criterion, the maximum relative slip amplitude and the maximum tangential stress range criterion to predict the location of fretting fatigue crack nucleation and the fretting fatigue crack path with cylindrical-onflat contact. In addition, using the effective stress intensity factor range estimated the fretting fatigue lives. In the same year, Tongyan et al. [16] explained the effects of roughness on fretting fatigue. They indicate that the higher roughness of surface reduces the life of fretting fatigue. Rodrigo et al. [24] studied the behavior of two cables (1350-H19 aluminum) in contact; subjected to the phenomenon of fretting fatigue using a numerical method such as a 3D finite element model, they also used an experimental method on the MTS3222.21 machine. Kyvia et al. [25] published a very interesting article on aspects of fretting fatigue, based on the experiences, analyzes and research of researchers in this field. He was able to present an article to fully understand the parameters that relate to fretting fatigue such as cylindrical contact, crack nucleation, crack propagation, and fretting fatigue modeling. In this work, stress Intensity factors are calculated by the Abaqus code using the EDI calculation [22]. Fortran program allows to create a parametric mesh using the Stretching Finite Element Method SFEM [18]. This mesh permits to change various model parameters, such as the applied load, mesh sizes, the shape and the dimensions of the crack as well as the contact parameters. The problem is solved by the Abaqus computer code to determine the main crack parameter, the stress intensity factor K, and its influence on fretting fatigue behavior.

FRETTING FATIGUE CONTACT TYPE
s Fig. 1 shows, the contact types can be classified depending on the geometry of the contacting bodies [17]. The first type is plan/plan contact (Fig. 1a) or the complete contact witch a flat contact surface is used for fretting pad [17]. This type produces contact area that is independent of load. Complete contact condition necessarily has singularities at its outside edges [18], where the stress is theoretically infinite, hence, plastic deformation is expected to occur at the pad edge. Fig. 2 presents the model utilized in this investigation; it clarifies the surface contact formed between the specimen and pad, the crack position and loads applied in two bodies. The second type is incomplete contact ( Fig. 1b and 1c) where at least of the two bodies is cylinder or spherical, which produces a contact area whose size depends on the magnitude the normal force magnitude [18]. The stress intensity factor is determined by considering the displacements and stress field near the crack front. It describes the stress state at crack tip. Historically, stress intensity factors K have been obtained from elastic solutions by considering 3D problems, under plain stress/strain condition. In the case of finite dimensional structure, i K (i = 1, 2 and 3 for mode 1, 2 and3) can be expressed in the following general form: where i F is the stress correction factor or form factor, taking into account the geometry of the crack as well as the type of stress applied on the structure, σ is the applied stress, and a is the crack length or the crack depth [20].

EQUIVALENT DOMAIN INTEGRAL (EDI)
he integral J is a measure of variation of potential energy, which in linear elasticity, can be directly related to the stress intensity factor SIF [21]: W is the strain energy density, Г is the contour around the end of the crack,  i u displacement vector at a point of the contour, i t stress vector given by  i t = ij j n . For an elastic material, this integral is identical to the rate of energy release G [19] where: is the shear modulus The Equivalent Domain Integral method use all the area A* showed in Fig.3, it is well applied to solve problems of solid in 3D with crack [18].
q is a regularized virtual displacement of the points of A * in direction x1 between a maximum value in Γ0 and a zero minimum value in Γ1 [2].

STRETCHING FINITE ELEMENT METHOD (SFEM)
he finite element method by stretching of the mesh SFEM is a new method; Bentahar et al. [3] have used it for a two dimensional crack propagation. Using a Fortran program to generate a parametric mesh, this non-remeshing method allows us to stretch the mesh elements with giving new node coordinates at each new mesh. [3]. A parametric mesh is created by using a fortran Code that keeps the same number of the elements and nodes but their coordinates are changed at each new mesh parameters (mesh size, crack orientation, crack size and singularity zone size). As Fig. 4 shows, this mesh is composed of the specimen and the pad, and they are exposed to two external forces σ and P, which means, the mesh is deformed during each new mesh parameters. In order to simulate the singularity that characterized the crack point, we use singular elements called 'quarter point' around the front of the crack, that is quadratic elements collapsed to obtain a triangular element. T The modeling of contact with Abaqus is based on the interaction of surfaces. The contact between the two bodies, specimen and pad is called an 'even contact' where a surface can interact on the other. A contact pair is defined by a slave surface (pad) which is the deformable surface, and the master surface (specimen) which is the most rigid surface, and the direction of contact is always perpendicular to this area. For each node of the slave surface, Abaqus looks for the nearest point on the master surface. The interaction is then discretized between the points obtained on the master surface and the nodes of the slave surface [2]. A contact pair is defined by the commands: *CONTACT PAIR, INTERACTION=interaction-name surface-slave, surface-master [2].

NUMERICAL MODEL
parametric mesh is composed of two bodies, which are in complete contact: the specimen that contains the crack is subjected to a tensile load σ = 100 MPa, and the pad is subjected to both a small displacement oscillation and a pressure P=10 MPa.The mesh for different sizes W/t = 1.0 with t = 25 mm. Tab. 1 shows different materials for the two bodies that are considered in the present study.  Due to the symmetry, only half of the structure, specimen and pad, is modeled. A fortran code is used to create a parametric mesh, which allows changing the values of the shape and dimension parameters. Abaqus 6.14 software is used to resolve the problem. The results of the stress intensity factor, for the three modes 1, 2 and 3, are obtained numerically from Abaqus.

RESULTS AND DISCUSSION
baqus calculates the stress intensity factor for 5 contours which are shown in Fig. 5. Each contour has 17 values from the integral J=1 to J=17 according to the angular division number n. We took the stress intensity factors SIF values of the third contour avoiding the first contour because of the singularity where the stresses are infinite and the more we move away from the point of singularity, the more we lose in precision. Only the second values and third values of the contour can give good results. For each contact parameter friction parameters f, mesh size H, plane inclination α, singularity zone size L, stress σ, materials steel, aluminum and brass), the results of the stress intensity factor for modes 1, 2 and 3, are represented as a function of angular division n = 17, see Fig. 7, Fig. 8, Fig. 10, Fig. 11, Fig. 12.

COEFFICIENT OF FRICTION
he effect of the coefficient of friction, ranging between 0.1 and 1.0 is investigated. Fig. 7 presents three graphs for three different materials, graph A for steel, graph B for aluminum and graph C for brass. For each graph, the evolution of the intensity factors of stress Ki / K 0 (Mode K1, K2 and K3) is presented as a function of the angular division n for different values of the friction coefficient. For the three different materials, the K1 mode has the greatest value; it is the most important mode in this study because it is considered as the most damaged mode because of the tensile force applied to the mesh. The more the coefficient of friction increase from 0.1 to 1.0, the more the value of K1 decreases. Aluminum has the greatest value of K1 / K0 = 20.36 at n = 7, followed by brass with K1/ K0 = 20.28 at n = 7, and steel K1/ K0 = 19.99 at n = 6. Regarding the second K II mode and for each material, its evolution as a function of the number of division n increases from n=1 to n=17, but the friction effect is weaker the same observation is made for third K3 mode. It is also observed that the value of the second mode K2 between the two points n = 1 and n = 2 is greater for aluminum and brass compared to steel. It is also noted that the first node of the zone of singularity, the value is not correct and that for the reason that this zone the constraint is infinite.      In the K2 mode curve, the values increase as a function of the number of angular division n, therefore we observe two parts, the first is before n = 4 or the values of K2 decrease each time the angle increases in contrast. The second part the evolution becomes proportional. The K3 mode has an inverse relationship as a function of the angle of inclination of the crack plane. Note that between the three materials aluminum has the highest value of (K1/ K 0 = 29.30; n=3).     Fig. 4 presents the singular zone dimension in the mesh; this zone is located at the front of the crack. It is represented by the dimension L considered as an important parameter to study its effect on the evolution of the stress intensity factor SIF. The values of L (mm) vary as {0.4, 0.5 and 0.6} where the other parameters are fixed at H = 8 mm, f = 1.0, α = 15 ° and σ = 100 MPa. The results obtained are presented in Fig. 11 for steel, Aluminum and Brass. The results of K i / K 0 concerning the effect of the singularity zone size L are presented in graphs A for steel, B for Aluminum and C for Brass. For each graph, the results are almost identical depending on the variation of L from 0.4 to 0.6. The evolution of K1 mode for the three materials increases for the first points as a function of the number of angular divisions in succession, it remains constant. It is observed that the mode curve K1 for the value of L = 0.6 mm is greater than L = 0.5 mm and L = 0.4 mm. That is to say the values of SIF in the curve increase with the increase in the dimension of the singularity zone L. Note that the greatest value is recorded for aluminum (K1/ K0 = 19.77 at n = 7 ). The values of stress intensity factor SIF in the K2 mode curve for the three material tests increase with the angular division angle n. We also observe that the effect of the dimension of the zone of singularity L is not very visible, up to point n = 9

EFFECT OF THE LOAD
n this case, the effect of external forces on the evolution of the stress intensity factor SIF was studied. We always note that σ is the tensile force, and p is the pressure force of the pad on the specimen and this is the force that is created by the fretting phenomenon. This force is very small compared to the tensile force, that is why we have fixed the relation between the two forces by the following relation σ = 10 p, 10 times greater. The other parameters remain constant during the three materials tests (f = 1.0, H = 8 mm, L = 0 .6 mm and α = 45 °). The results obtained are shown in Fig. 12.   Fig. 12 the effect of the tensile force σ on the evolution of the stress intensity factor SIF in the three graphs appears clearly. These graphs show a proportional relationship between the change in mode K1 and the tensile force values σ for the three material tests the highest K1 mode value appears in the aluminum graph B (K1=138.47 at n=7). The evolution of the second mode of K2 as a function of the angular division number n is increasing. We observe two parts of the evolution on the three graphs. The first is before n = 6, where the evolution of mode K2 decreases each time the tensile force σ increase from 100 to 700 Mpa. The second part is after n = 6 , where the relation becomes inverse. For Mode K3 between n = 2 and n = 16, the evolution of the three curves increases as a function of the number of divisions n while the evolution as a function of the tensile force σ decreases each time the force σ increases by 100 at 700 MPa.

CONCLUSIONS
numerical simulation of parametric mesh witch composed of two bodies in full contact is developed in this study. for that, a parametric mesh is analyzed with the Abaqus 6.14 software to study the effect of the different contact parameters of fretting fatigue phenomenon. This model is subjected to two cyclic loads, the first is a tensile stress, it is applied to the specimen and the second is a perpendicular pressure on the pad. In this study, we calculated the stress intensity factor in three modes K1, K2 and K3. The results are presented on graphs, which show the evolution of the factors as a function of the angle of division n for the different parameters. From the results obtained, the parameters are divided into two groups. The first group (the size mesh H, the angle of inclination of the plane of crack α and the tensile force σ) all have a significant effect on the evolution of the stress intensity factor SIF, and the second group (coefficient of friction and the dimension of the singularity zone L) have a weaker effect. Through the results obtained, K1 mode shows the higher value for the three material tests and for all different parameters. It is noted that the following parameters: the size mesh H, the inclination of plane of the crack α and the coefficient of friction have an inverse effect on the evolution of mode K1. Regarding the material effect, although the results obtained for the three material tests are close, Aluminum is more influential on the evolution of K1 mode. We conclude that the parameters which have been applied externally (σ, H and, α) have an efficiency compared to the parameters which have been applied internally (f, L and the different materials). Finally, with this SFEM method we can study a large field (range) of parameters and compare several results at the same time which allows us to make a rather rich study.