Orientation of propagating crack paths emanating from fretting-fatigue contact problems

In this work, the orientation and propagation of cracks in fretting fatigue problems is analyzed numerically using the finite element method (FEM) and the extended finite element method (X-FEM). The analysis is performed by means of a 2D model of a complete-contact fretting problem, consisting of two square indenters pressed onto a specimen subjected to cyclic fatigue. For the simulation, we allow for crack face contact in the implementation during the corresponding parts of the fatigue cycle. The problem is highly nonlinear and non-proportional and we make use of the so-called minimum shear stress range orientation criterion, min(), proposed by the authors in previous works. This criterion is introduced to predict the crack path in each step of the crack growth simulation. The objective of the work is to detect which is the relevant parameter affecting the crack path orientation. A parametric study of some a priori relevant magnitudes is carried out, such as normal load on the indenters, bulk load on the specimen, stress ratio and relative stiffness of the indenter and specimen materials. Contrary to previous expectations, it is shown that the relative magnitude of the applied loads has no significant effect. However, it is found that the stiffness of the indenter material with respect to the specimen material has the greatest effect. A simple explanation of this behavior is also provided.


INTRODUCTION
retting fatigue problems involve two or more solids in contact that experience relative displacements of small amplitude.The contact region acts as a stress raiser causing crack initiation and subsequent crack propagation until the eventual failure.In these problems, loads are usually nonproportional, since the normal load acting on the F contacting body is constant with time, whereas the substrate is subjected to cyclic fatigue load.As shown in experimental tests carried out by the authors [1] (see Fig. 1) and in many works in the literature, cracks emanating from the edge of a contact pressing onto a surface tend to grow with a slight deviation inwards beneath the contact and not fully perpendicular to the applied bulk stress.This slight deviation from the normal direction cannot be predicted using a conventional orientation criterion, such as the maximum tangential stress criterion (MTS) and this is the main motivation of this research.
Figure 1: Micrographs for the propagation of non-failure cracks of four tests [1], emanating from the edge of contact.
As a result of previous investigations [1], the nonproportional character of the load evolution precludes the application of conventional crack orientation criteria in LEFM for propagating cracks.In a previous work [1], we proposed an orientation criterion based on the minimum range of the shear stress ahead the crack tip, min().This criterion has been successfully applied to the prediction of the orientation of long cracks in nonproportional configurations found in fretting fatigue.
Continuing with this study, we analyze in this work the influence of several a priori relevant parameters on the crack orientation, such as the normal contacting load, fatigue bulk load, stress ratio R, stiffness of the materials involved, etc.We have found that the amount of contacting normal load is not a relevant magnitude affecting the crack orientation.Contrary to previous expectations, the relative difference in stiffness between indenter and substrate is the parameter that has the greatest effect on the direction of crack propagation.The numerical analysis is performed using the standard finite element method (FEM) and extended finite element method (X-FEM), taking into account the crack face contact along the loading cycle, as implemented by the authors [2].To predict the crack direction in each step of the crack growth simulation, we use the criterion of min() [1], as reviewed in the next section.

THE CRITERION OF THE MINIMUM SHEAR STRESS RANGE
rom the numerical analyses and for the geometric and loading configuration considered in this work, it is found that the crack remains closed during a large part of the loading cycle.The application of some of the criteria previously reported and reviewed in [1] did not lead to good predictions of the actual crack path as observed in the experimental tests performed.Consequently, it can be argued that the stress state existing under a crack face contact condition has an important influence to be considered.Assuming an elastic behaviour, the stress state under crack face contact conditions must be essentially controlled by K II , the only stress intensity factor that can exist for a totally closed crack in 2D.The criterion proposed here is a generalization for non-proportional loading of the so-called "criterion of local symmetry" well established for proportional loading, see Cotterell and Rice [3].The criterion of local symmetry states that the crack F will propagate in the direction that KII=0.For non-proportional loading, the condition of KII=0 cannot be reached in general, and therefore, the proposed criterion seeks the angle for which the range KII is minimized.This obviously reduces to the condition KII=0 when applied to proportional loading problems.We note in passing that, for proportional loading, the criteria of K II =0, Nuismer and MTS lead to the same result [3,4].In practice, computing KII values under crack face contact must include the effect of friction tractions on crack faces, as in [2,5], which can be cumbersome and prone to inaccuracies when using domain and contour integrals.Instead, and equivalently, we will seek the angle for which the shear stress range  at the crack tip is minimized.Shear stresses develop always in two orthogonal planes and there are two orthogonal planes on which  is minimum, min().From these two potential crack growth directions, we choose the plane with the maximum  n , because it will be the plane where less frictional energy is lost and there is more energy available for propagating the crack.This approach is in line with the principle that a crack will grow in the direction which maximizes the strain energy release rate G [3,4].As verified in the results section, the min() direction coincides with the direction of the maximum range of normal stress, max(n).This is due to the in-plane stress tensor transformation that yields both extremes in the same direction, although this may not be the general case.However, the direction predicted by the maximum range of the effective normal stress, max( n,eff ), does not lead to good results, at least in the problems studied here, despite the intuitive idea that only the positive normal stresses (effective) will govern the crack behaviour under an elastic material behaviour.Fig. 2 sketches the convention used in the procedure.When applying the proposed criterion to the propagation stage,  is evaluated ahead the current crack tip.For each crack growth increment, stresses are evaluated ahead the crack tip and the prospective local direction  is searched for which  is minimum (see example of the estimation of the third increment direction in Fig. 2, left).In the results section, the predicted angle is reported with respect to a fixed reference:  is the predicted angle measured from the specimen surface.This way, a crack segment growing inwards (with respect to the indenter contact zone) has an angle 0º< <90º and -90º< <0º indicates a crack segment growing outwards.

Crack in a simple bar subjected to tension-compression cycle
As an example of a simple application of the minimum shear stress range criterion, a cracked bar of uniform section loaded in tension is analyzed.Fig. 3 shows the geometry and loads of the model and a contour plot of the von Mises stress field.This preliminary analysis is performed using standard FEM with ABAQUS (no X-FEM is considered at this stage).The tensile load is cyclic with R = -1 and is applied in the same fashion as the blue curve of Fig. 6 (no indenter load exists in this simple example).Fig. 4, left, shows the variation of the normal stress  (in blue) and shear stress  (in red) on a plane forming an angle  with respect to the horizontal surface.Stresses are evaluated at the finite element located ahead the crack tip and transformed according to the angle  of the prospective plane.The successive curves plot the variation along the time of the last step (step 4).Note that the normal stress  is maximum at the end of the step (curves located at the top of the figure) and at a plane orientation of ±90º, as expected.This tensile stress is much higher than the corresponding compressive stress due to the effect of the crack opening (mode I of fracture), whereas the closing stage does not concentrate such high stresses (in magnitude).The evolution of the shear stresses (in red) is analogous.Note that the shear stresses  are zero when the normal stress is maximum or minimum (for any given increment of time).The fact that all the maxima and minima are attained at the same angles (±90º for  and ±45º for ) is indicative that the stresses evolve in a proportional way.Should the maxima and minima be shifted along the  axis, then it means that the loads are nonproportional.Fig. 4 right, shows the application of the min() criterion.The same  curves of Fig. 4, left, are replotted and the maximum and minimum with time along the step 4 are marked in black.Then, the range of variation of  is computed simply as  max - min .The minimum shear stress range criterion predicts that the prospective propagation angles are either 0º or 90º.Due to the nature of the shear stresses, there are always two prospective angles with a difference of 90º.The discrimination between both angles is done choosing the angle that also leads to the maximum normal stress to that plane.The method sharply indicates that the predicted angle of propagation for this case is 90º, as expected in such a simple problem.

NUMERICAL MODEL
ue to symmetry conditions, a quarter 2D finite element model has been considered to represent the fretting fatigue tests, as shown in Fig. 5.The rectangle LxB corresponds to the portion of the analyzed specimen and has a length of L = 4B = 20 mm, the half length of the indenter c is 5 mm, and the distance between the contact plane and the point of the indenter at which loads are applied is h = 10 mm.Four node, plane strain quadrilateral elements were used with a thickness t = 5 mm.The smallest element size considered is 5m at the right end of the contact zone.The friction model assumed for the contact zone is a Coulomb model and the ABAQUS contact formulation based on Lagrange multipliers is used to model the contact between the indenter and the specimen.The friction coefficient is taken as  = 0.8 [6].The material behaviour is assumed linear elastic, despite the high stress concentration at the contact edge.

D
The specimen material is an aluminium alloy 7075-T6, with E = 72GPa,  = 0.3.For some of the cases analyzed, the indenter material is changed according to Tab. 2 included in the results section.The application of the linear regime is deemed valid, due to the very small edge radius of the indenter and the relative high yield stress of the aluminium alloy.As a consequence, the existing plasticity is very localized and a small scale yielding assumption can be applied, analogous to the small scale yielding assumption admitted in LEFM around the crack tip.This is confirmed by the observation of the tested specimens, which showed no macroscopic evidence of plasticity (see micrographs in Fig. 1 in which crack faces match very well each other and a view of the specimen contact surface in Fig. 7 of our previous work [7]).The loading sequence is represented in Fig. 6 for one of the examples analyzed (case 3, see results section), where four load steps have been considered in the analysis.Due to the non-linearity of the contact problem, loads were applied in sufficiently small time increments.At time t = 2.00 the maximum  Bulk is being applied, which produces a clear opening of the crack.When the bulk load is decreased in the first half of step 3 (2.00 < t < 2.50), mode I is reduced and a clear mixed mode condition appears, which has been observed through FE analyses.Note that the vertical load due to the indenter is kept constant during the cycle and mode II increases its dominance over mode I as  Bulk is reduced.At time t = 2.50 crack face contact is produced and a mode II condition is present at the crack tip.At time t = 3.00 the bulk load is completely reversed (since the stress ratio in case 3 is R=-1) and the load is transmitted through the crack faces.The end of the contact zone acts now as a strong stress raiser, as the specimen is compressed against the contact corner.Results in the following section are presented for the load step 4 (3.00 < t < 4.00), when ''shakedown'' of the numerical model response is produced).It has been verified that the stress states at t = 3.50 and t = 4.00 are very similar to those at t = 2.50 and t = 2.00, respectively.

RESULTS
he geometrical model with an initial crack of length a 0 =0.3 mm and initial orientation of  =90º has been analyzed under 24 different cases.The aim is to detect those parameters that have the largest influence in the crack orientation, since the growth is always directed slightly inwards beneath the contact zone, as shown in Fig. 1.

Variation of normal load, bulk load and stress ratio
Tab. 1 shows 13 load cases performed changing either the constant (with time) normal load P applied on the indenter and the corresponding P, the variable bulk load on the specimen Bulk and the stress ratio R. The material stiffness is 72 GPa in all cases, both for the indenter and specimen.The last column indicates the predicted angle using the minimum shear stress range criterion, min().Contrary to was initially expected, the first thing that draws attention is that there is no practical variation on the predicted angle, since all cases lead to an orientation angle of 78º-79º.Even for the cases with negligible contacting normal load,  P =10 -6 , the prediction leads to angles pointing inwards.The influence of the wide ranges tested for  Bulk and R is also negligible.This is in full agreement with the experimental evidence collected by the authors [1,7], summarized in Fig. 1   Fig. 7 shows the variation of  versus the prospective crack orientation angle  for the last step of the loading cycle.This enables the application of the minimum shear stress range criterion.Fig. 7, left, shows the results for case 1.The high proportionality of the loads is demonstrated by the same location of the maxima and minima (no shifting of the curves).The load proportionality is caused by the extremely low value of  P considered in this case 1.However, even under this T situation, the effect of the indenter causes the deflection of the crack to 79º.The results for the case 7 are presented in Fig. 7, middle.Here the nonproportionality is evident due to the high value of P, which is equal to Bulk, being the curves shifted ones with respect to the other.Due to the dominance of the constant normal load, the effect of the cyclic load is less evident and the range between  min and  max is not so important.This range is even less for case 8 (Fig. 7, right), due to the change in R from -1 to 0. In all cases, the predicted angle is around 79º.We remark that by application of a conventional orientation criterion, such as the MTS, at the instant of maximum bulk load, an incorrect prediction of the crack direction is obtained (pointing outwards, see [1]).Other examples of inaccurate growth orientations using the MTS criterion under non-proportional fretting loading can be found in Figs. 6 and 7 of one of our former works [8].In that work, the MTS criterion was applied at the instant of maximum bulk load.All of our analyses predicted growth paths outwards the contact zone, which is in contradiction with experimental evidence.In [9], the application of the MTS criterion to a fretting fatigue analysis at the instant of maximum load is also reported, yielding growth directions outwards the contact zone as well.

Variation of stiffness of the indenter and specimen materials
Tab. 2 presents the results obtained when changing the material stiffness of the indenter, i.e. considering the indenter and specimen as made of dissimilar materials.It can be seen that the effect is much more evident that in the above cases.Despite the range of variation is not large (between 76º and 90º), this proves that the stiffness of the indenter (relative to the stiffness of the specimen) is the main cause of the crack deflection inwards the contact zone and allows to gain insight into the mechanisms that cause this effect.

Case
Eindenter   Cases 14 and 15 illustrate well this behaviour and the application of the min () to them is shown in Fig. 8.In case 14 a negligible stiffness is given to the indenter.The behaviour is in practice proportional, as if there is no indenter.The predicted angle is 90º.In case 15 the same loads are considered, but the Young's modulus of the indenter is changed to steel.The nonproportionality of the behaviour is more noticeable and the predicted angle is 78º.Fig. 9 shows cases 16 and 22 (E indenter =20000 MPa and E indenter =1000000 MPa, respectively) with a non-negligible normal load on the indenter.It is verified that the larger the material stiffness of the indenter, the larger the deflection (84ºand 77º, respectively).The reason governing crack deflection due to the material stiffness can be explained as follows.The indenter acts as a contacting solid through which the force lines escape and deviate, just due to its stiffness and geometry, since a stiff solid tends to carry a higher load than a compliant solid (assuming a parallel configuration).This can be visualized by the directions followed by the maximum principal stresses shown in Fig. 10 for case 15.It can be seen that the principal directions (that can be assimilated to local force lines) tend to escape to the indenter just behind the crack.Therefore, it is expected a deflection angle approximately normal to the directions of the force lines in this region.This angle has been correctly estimated by minimum shear stress range criterion.Note also that the amount of deflection becomes saturated, despite a high increase of E indenter , due to the geometric configuration of the model that does not allow further deviation of the force lines.

Crack propagation using X-FEM
In this section, the extended finite element method [10] is used to model crack propagation.The great advantage of the X-FEM method is that the crack faces do not need to conform to the element sides of a mesh.Therefore, a single mesh can be used for virtually any arbitrary crack intersecting the mesh.This avoids remeshing and it becomes especially useful when modeling crack propagation in fatigue problems.This is accomplished through a special mathematical formulation of the FE method that includes the enrichment of the standard finite elements with additional degrees of freedom (DOFs) at the nodes.These additional DOFs are associated with the nodes of the elements that are geometrically intersected by the crack location (called enriched nodes and elements, respectively).Thus, the discontinuity is included in the numerical model without modifying the discretization.The X-FEM formulation allows for a further type of enrichment for those nodes that surround the crack-tip.These nodes are enriched with additional DOFs to represent the first term of the classical Williams series expansion in linear elastic fracture mechanics in terms of the displacement field.We have applied the X-FEM method to two problems with different indenter materials.The aim is to predict numerically the crack paths and the analysis is carried out using the X-FEM implementation developed by the authors [11,12].The implementation is performed as a user's subroutine in the commercial code ABAQUS and can take into account crack face contacts along the loading cycle, which have been proved to be essential for the correct crack prediction.Fig. 11 shows the result of the crack propagation in the vicinity of the contact corner after 10 increments using the min() criterion after each increment to predict the orientation of the next crack segment.In Fig. 11, left, the Young's modulus for the indenter is Eindenter= 10 3 MPa, whereas in the simulation of Fig. 11, right, Eindenter=10 8 MPa.In both cases Especimen= 72•10 3 MPa.As expected, it can be observed that the crack growth path is clearly deviated inwards when the material stiffness of the indenter is large with respect to the specimen.The procedure presented in this work enables the prediction of more accurate crack paths in components that can be used for the integration of crack growth laws in order to estimate fatigue lives.

CONCLUSIONS
rack propagation directions and paths have been predicted for fretting fatigue tests under complete contact conditions.This type of problem is subjected to non-proportional loading, which invalidates the application of conventional orientation criteria usual in LEFM and that are only useful for proportional loading.The criterion of the minimum shear stress range has been proposed, and it is based on the minimum value of  evaluated ahead the crack tip and along the entire cycle.The prediction has been performed numerically using FEM and X-FEM including a formulation that allows for crack face contact, which is essential to take into account the effects during the compressive part of the cycle.A parametric study has been accomplished by variation of the normal load on the indenter, the bulk load on the specimen, the stress ratio and the elasticity modulus of the indenter.It has been shown that the parameters related to the loading have no effect on the crack deflection, whereas the change in the material stiffness of the indenter has a significant effect on the predicted crack path direction.The numerical results are in good agreement with the experimental observations, validating the procedure.This can lead to more accurate fatigue life estimations once the crack path is predicted using the proposed procedure.

Figure 2 :
Figure2: Application of the min() criterion to predict the third crack-growth increment direction.Sign convention for direction angles of a crack growth increment.

Figure 3 :
Figure 3: Left, geometry and loads of the model of a cracked bar in tension.Right, detailed view of a von Mises contour plot.

Figure 4 :
Figure 4: Left, variation of the normal stress  (in blue) and shear stress  (in red) on a plane forming an angle  with respect to the horizontal surface.Stresses are evaluated at the finite element located ahead the crack tip.Right, application of the min() criterion.

Figure 5 :
Figure 5: Left, model geometry.Right, complete contact testing rig, showing the contact elements.

Figure 6 :
Figure 6: Loads applied to the numerical model for one of the cases studied (case 3).Evolution with time.

Figure 10 :
Figure 10: Maximum principal stress directions for case 15.

Table 1 :
, with growing directions around 79º. Predicted orientation angles for different load cases, generated by variation of P, Bulk and R.

Table 2 :
Predicted orientation angles for different material stiffness of the indenter and some variations of  P and  Bulk .