Numerical analysis of the influence of liquid on propagation of a rolling contact fatigue crack

Numerical investigations of the propagation of rolling contact fatigue crack filled by the liquid have been conducted. Two models of fluid – crack interaction have been considered. In the first model called “hydrostatic” the assumption of incompressible, inviscid and weightless liquid was accepted. It was also assumed that due to the wheel load the trapped liquid could not get outside the crack and its volume remained constant until the rising pressure would open up the crack mouth again. On this assumption the analysis has a steady-state character. In the second model it has been assumed that the crack is filled by the viscous, incompressible fluid and the fluid motion as well as the resulting pressure distribution can be represented by one-dimensional form of the Reynolds equation. The method for solving the problem of the coupled motion of liquid and crack faces has been developed and series of calculation were made. The method has been employed for the predicting of crack deformation in the course of wheel rolling.


INTRODUCTION
he research presented in the paper is part of more comprehensive investigations aimed at understanding the propagation mechanisms of rolling contact fatigue cracks in railway rails.Liquid penetrating the crack interior may be considered as a very important factor affecting the crack propagation.The effect of liquid influence has been already investigated; e.g., M. Kaneta and Y. Murakami [1] considered a 3D model representing a semi-elliptical surface crack situated in the elastic half space.The load had a form of the Hertz stress of T ellipsoidal distribution moving along the half space.The effect of liquid presence in the crack was considered in terms of stress the magnitude of which was equal to the pressure acting on the crack mouth.Variations of non-dimensional stress intensity factors as functions of the load position were analysed.The contact between the crack faces was neglected, therefore the negative values of KI factor appeared.A.F. Bower [2] analysed a 2D model of a surface crack situated in the elastic half space using a dislocation distributions finding analytical method.The Hertz stress replaced the wheel action.The model considered included the contact effect and the tangential component of interaction between the crack faces as well as the tangential loads acting on the surface.The liquid interaction was considered using the following three different ways: (a) reducing the value of friction coefficient on the crack faces; (b) introducing a hydraulic pressure of the magnitude equal to the pressure acting on the crack mouth; (c) assuming a constant volume of the crack interior after the pressure acting on the crack mouth has closed it.Variations of stress intensity factors were analysed as functions of different parameters of the model and the direction of load motion.In this paper the Authors have once again attempted to define the role of liquid in the development of rolling contact fatigue (RCF) cracks in railway rails.The effect of liquid has been represented as a crack face pressure distribution.This distribution has been determined taking the following two ways (a) a hydrostatic approach in which a constant volume of liquid is assumed to be trapped in the crack interior and (b) a hydrodynamic approach based on the equations of viscous fluid motion.Additionally, both the wheel and rail have been modelled as finite, linear elastic solid bodies.The Authors used a unique algorithm that employs the flexibility matrices of the bodies in contact, which were determined using FEM [3].The algorithm employed allows for solving 2D and 3D problems of the wheel-rail contact taking into account the contact between the crack faces well as introducing many other phenomena, like pressure of the liquid filling the crack, wheel loads acting upon the contact zone, thermal loads, residual stresses in the rail, rail bending, variable friction coefficients over the contact surface and different types of the crack faces interaction [4,5].

THE MODEL DESCRIPTION
art of the analysis reported here was based on considerations of the 2D model, in which a rectangular prism with an oblique crack represented a rail head segment with a crack, and a cylinder represented the rail wheel.The value of cylinder radius was equal to that of the wheel.The magnitude of normal load R was selected in the way ensuring that maximal magnitudes of pressure in the cylinder-prism contact zone (remote from the crack) be equal to about 870 MPa, i.e. the magnitude revealed in the real wheel-rail contact zone.The values of Young modulus and Poisson's ratio assumed were equal to those revealed by steel.The boundary conditions for the prism are shown in Fig. 1.The conditions of plane strain-state were imposed on the model.In the course of investigations presented below the crack had a length of a = 5.81 mm and was inclined at α = 25.The existence of tangential forces was assumed in the contact zones, determined by the friction coefficients -μ p in the prism-cylinder contact zone and μ s at the crack faces, respectively.The current position of cylinder was represented by the coordinate x and the value x = 0 corresponded to the point at which the crack intercepted the rail raceway.On the diagrams the cylinder position is represented by the non-dimensional parameter x/b, where b = 6.7 mm is half of the contact length calculated using the Hertz formulae for a cylinder and a half space.The investigations conducted have proved that as the cylinder comes closer to the crack its mouth tends to open allowing the liquid to penetrate the crack interior (Fig. 2a) When the rolling wheel closes the crack mouth the space is formed in which the liquid is entrapped (Fig. 2b).

LIQUID MODELS
ne of the main problem is the use of a suitable fluid model.This model should provide high compatibility with reality and, on the other hand, be easy to solve.There are usually two main groups of models: hydrostatic and hydrodynamic.

Hydrostatic model
The first attempt the Authors made at investigations, in which the liquid presence in the crack was considered dates back to 1993 [6,7].In the analysis carried out the assumption of incompressible, inviscid and weightless liquid was accepted.It was also assumed that after the crack mouth had closed due to the wheel load the entrapped liquid could not get outside the crack and its volume remained constant until the rising pressure would open up the crack mouth again.This assumption was rather unrealistic since the real crack faces reveal rough surfaces making a tight closing impossible, the analysis however aimed at the determination of the limit of pressure and amplitudes of stress intensity factor variations that might be affected by the liquid presence in the crack.The aforementioned assumptions forced the quasi-static nature of the analysis, i.e. the time was not a parameter of the model.For that reason the Authors called their model the hydrostatic one.Additionally, zero pressure at the crack mouth was assumed, and in the "bubble", i.e. in the closed space with no contact, the magnitude of pressure was independent of a position, while along the contact zone in the crack interior the pressure magnitude was changed linearly from zero value to that in the "bubble".Similar models are still used in the calculations by different authors for their simplicity and quick calculation time.They assume a pressure distribution between a cylinder and a prism calculated according to the Hertzian theory [8 -11] or, like in the authors' solution, pressure distribution during rolling is calculated on an ongoing basis [12].Despite the fact that the assumption accepted for the hydrostatic model make it not very realistic and the results obtained cannot be fully reliable, one can get, however, an impression of how strongly the liquid could affect the propagation of crack.

Hydrodynamic model
The Authors have decided to apply the model in which the liquid filling the crack interior is represented in a more realistic way i.e. with its viscosity included.A similar approach was used at work [13].
Since the crack height is very small as compared to its length a 1D model was assumed for representation of a liquid flow in the crack.As the crack front is closed changes appearing in the crack height due to the wheel rolling generate a flow along the crack.The velocities of the crack faces directed along the crack were neglected.On the assumption that the crack is filled with incompressible liquid one can represent a pressure distribution using the Reynolds equations in the form: O where: p -liquid pressure h -local crack height η -coefficient of dynamic viscosity s -coordinate directed along the crack t -time The equation is completed by the following boundary conditions: It can be easily seen that the pressure distribution depends on the crack height and the speed at which the crack faces come closer to each other.The pressure magnitude at the crack mouth equals to the ambient pressure p amb while its value at the crack front results from the condition of zero pressure gradient that should be satisfied there.Since the applied algorithm for solving the contact problems deals with discrete systems it would be more convenient to transform the Reynolds equation also to a discrete form.Using the planes perpendicular to the crack plane one can divide it into the segments that correspond to those employed in the algorithm for solving the contact problems (Fig. 3).The sectional planes bear the numbers within (1, n).All quantities appearing in a given section have the indices equal to the plane number (Fig. 4).The section position is determined by the s-coordinate measured along the crack length.After replacing the derivatives in Eq (1) with the corresponding finite quantities the Reynolds equation has the form where: After transformation and substitution for being the velocity of crack faces motion we have   The discrete boundary conditions have the form After writing Eq (3) for all sections of the crack we obtain the following set of linear equations with a 3-diagonal matrix where: and from the boundary conditions we have For the sake of simplicity the ambient pressure p amb = 0 was assumed in calculations.The scheme of solution to the aforementioned problem is given below (the index in brackets represents the current iteration number while the index with no brackets stands for the section number).
Wheel path is divided into small fragments x small enough to fulfill assumption that the pressure can be treated as constant during motion along it.At the beginning of a new segment x pressure distribution along the crack and the crack height distribution are described by the vectors [p (0) ] and [h (0) ].A new pressure distribution vector [p (i) ] at i-th iteration is calculated in following way: a) Taking the pressure distribution vector [p(i-1)], the crack height distribution vector [h(i)] is calculated for the wheel position at the end of segment x b) First vector of crack faces transversal velocity [v(i)] is calculated from following formula: where: t is the time of wheel passing distance x c) Then by the help of Eqs. ( 4

), using vectors [v (i) ] and [h (i) ] vector describing the pressure distribution [p'
where: m -natural number from the range 10001000000 chosen experimentally (lower m-faster convergence, higher m -better stability) Modified pressure vector have been used to assure convergence of described iterative algorithm Iterative process is repeated until in all crack face nodes the pressure difference between actual pressure [p' (i) ] and modified pressure [p (i) ] is less than the arbitrary chosen value p eps = 0.002 MPa, i.e. .
If presented condition is fulfilled the pressure vector for the next wheel position is calculated.
In presented calculation segments x equal 0.025mm have been used.Two times reduction of the segment length x hasn't had any significant influence on the pressure distribution inside the crack.
The following assumptions were accepted in the model with liquid: -Liquid is weightless and incompressible -Liquid fills the whole crack interior even when subject to a negative pressure -Initial clearance between the crack faces changes linearly taking the value of 0.005 mm at the crack mouth and 0 at its front (existence of clearance have been find during investigation of real cracks) -Plane strain conditions -Tangential forces in the contact zones are neglected (μp = μs = 0) -Boundary conditions for the prism are shown in Fig. 2.
The following values of parameters were assumed: -Speed of wheel rolling 30 m/s (108 km/h) -Ambient pressure p amb) = 0 The dynamic viscosity of the pure water at a temperature of 20C has the value of η = 0.001 [Ns/m 2 ].But η varies increasing as the pressure increases and temperature decreases [14].Additionally, the water filling the crack can be oilcontaminated.Therefore, to ensure that the calculated stress intensity factors are not underestimated, a slightly higher value of viscosity has been accepted in calculations.

Comparison between the hydrostatic and the hydrodynamic model
o find out to how the form of the model used affects the obtained values of amplitudes of stress intensity factors the calculations were performed once again for the hydrostatic model, neglecting the tangential forces in contact zones and introducing a clearance between the crack faces and for hydrodynamic model.The results are presented on Fig. 5.In both cases the tangential forces in contact zones were neglected (friction coefficients μp = μs = 0.0) and a clearance between the crack faces changed linearly assuming the value of 0.005 mm at the crack mouth and 0 at its front, respectively  The value of amplitude KI obtained from the hydrodynamic model with viscous liquid was almost two times higher than the value obtained from the hydrostatic model with inviscid liquid while the value of amplitude K II remained unchanged.Substantial reduction in KI value observed when neglecting the tangential forces in the hydrostatic model and introducing at the same time a clearance in the crack resulted from both a smaller volume of liquid entrapped in the crack interior and reduction in the magnitude of wheel load acting on the rail.When the clearance appears in the crack its deflection arises due to a cylinder action and the load is transferred to the areas ahead of and behind the crack, respectively.

T
Results for the hydrodynamic model At the starting point of the analysis, i.e., when the cylinder is located ahead of the crack at a distance of x/b = -2.832,zero values of the pressure, crack faces speed and flows, respectively, are assumed.Despite the fact that the distance between the cylinder and the crack mouth equals almost 3b (18.975 mm) the crack height at its mouth due to the cylinder action reaches the value of about 10 μm.As the cylinder comes closer to the crack its faces open generating the negative pressure (down to -1.63 MPa) and the liquid is sucked into the crack interior.The process terminates when the cylinder comes into contact with the end of the wedge-shaped raceway part behind the crack.The process of crack mouth closing starts at that moment manifesting through high negative values of crack faces speed close to the crack mouth.As the cylinder moves on and the wedge-shaped raceway part behind the crack deflects the raceway points situated more far away from the crack come into contact with cylinder, while those close the crack mouth loos the contact, so the position at which the closing velocity of crack faces is highest (i.e.most negative) moves into interior of the crack.

CONCLUSIONS
n all the cases that have been investigated so far, with no liquid present in the crack, the amplitude of K II changes (KII ) dominates, while the value of KI remains low exerting only a weak influence on the fatigue crack propagation [2][3][4][5].The liquid presence in the crack causes a substantial growth of KI value, despite the model considered.The value of amplitude K I obtained from the hydrodynamic model with viscous liquid was almost two times higher than the value obtained from the hydrostatic model with inviscid liquid while the value of amplitude KII remained unchanged (Fig. 7).The tangential forces in the contact zones were neglected in both the models considered.It is known (see [5]) that the aspherities of the crack faces may come into contact creating a kind of shape joints that block the crack faces motion relative to each other, limiting therefore considerably the value of K II .If the liquid presence in the crack is considered KI may dominate and the crack may propagate by means of tearing of the crack front i.e. mode I of crack propagation may arise.It occurred in practice, however that the algorithm presented above reveals unsatisfactory convergence.When the cylinder is situated at a point far from the crack a few up to several dozen thousand of pressure iterations are required for the position change x = 0.025 mm.But when the cylinder comes above the crack the required number of pressure iterations grows up to several hundred thousand or even reaches a million.In view of the above it seems quite unreasonable to consider models of cracks filled with liquid with the tangential forces in the contact zones introduced into the model since a solution to the contact problem posed that way needs iterative procedures.To calculate the tangential forces and slips one should make several hundred up to several thousand iterations for each position change x.When using the model with the liquid and tangential forces introduced, after each iteration of tangential forces a smaller number of liquid pressure iterations is required than for the corresponding model with no tangential forces introduced, of course for the same cylinder position change x.One should, however, be conscious of the necessity for longer calculation time, even by over a dozen times.It should be noted that the analysis presented above still reveals some shortages.No limits imposed on negative pressure formation in liquid are introduced and as a result the negative pressure arising close to the crack front the magnitude of which reaches 449 MPa.The presence of negative pressure opens up the possibility that the phenomenon of cavitation may arise in the crack interior.Limiting the negative pressure value to -0.08MPa, i.e. the value at which water of a temperature of about 20 C starts to boil may restrain substantially the process of crack filling with liquid at the phase when the wheel comes closer to the crack.A smaller amount of liquid present in the crack results in a lower value of K I .Moreover, at the pressure of 800 MPa the liquid cannot be assumed as incompressible.Taking into account the liquid compressibility would probably result in additional reduction of the value of K I .
It should be noted that however the zero value of ambient pressure was assumed in the presented analysis, one should rather assume the hydrodynamic pressure that might be generated over the wheel-rail contact zone e.g.under wet weather conditions.In the future, the liquid model used here can be used to calculate the development of real cracks with geometry measured on real object.An example of such measurements can be found at work [15].I

Figure 1 :
Figure 1: Scheme of the 2D model of a rail with the crack of together with the rail wheel.

Figure 2 :
Figure 2: Part of a deformed FEM mesh representing the prism with the crack for two positions of the cylinder: a) x/b = -1.2corresponding to the maximum crack opening; b) x/b = -0.9-corresponding to the position at which the crack mouth is closed

Figure 3 :
Figure 3: Part of the FEM mesh showing the way the crack vicinity is divided into finite elements.The grid used in flow analysis is the same as boundary FEM nodes

Figure 4 :
Figure 4: The way of crack discretisation by means of its sectional division using n perpendicular planes If the vector of crack heights [h] = {h 1 , h 2 , ..., h n } is known ( i.e. the values h i are known for all sections) together with the crack faces velocity vector [v] = { v1, v2, ..., vn} one can determine the pressure vector [p] = { p1, p2, ..., pn} using Eq (4).In the case of crack in elastic material the pressure change [p] affects the height change [h] and the crack faces velocity change [v], respectively.Therefore, when dealing with the problem of a wheel rolling along a rail with a crack filled with liquid one should use iterative methods for calculation of the pressure distribution along the crack.

-
Crack length a = 5.81mm -Angle of the crack inclination α = 25 -For the material of prism and cylinder the Young modulus E = 210000 MPa and Poisson's ratio ν = 0.3 (the values for steel) -Cylinder radius r = 450 mm -Load acting upon the cylinder R = 9 kN/mm.-Coefficient of dynamic viscosity η = 0.01 [Ns/m 2 ]

Figure 5 :
Figure 5: Distributions of the stress intensity factors in the process of cylinder rolling along a prism with the crack filled with viscous liquid (hydrodynamic model) and inviscid liquid (hydrostatic model), respectively.

Fig. 6 Figure 6 :
Fig.6show the distributions of pressure, crack heights and liquid flow along the crack, respectively, for different cylinder positions obtain using hydrodynamic model.
Reduction of volume filled by liquid causes the pressure increase and induces the motion of liquid in direction of lowering pressure.Whereas, volume growth causes pressure decreases and makes possible of fluid inflow.Due to the use in calculation of simplified model of liquid omitting the cavitation (liquid evaporation), calculated pressure can drop down to values lower than zero.Such situation can be noticed in crack at wheel position equal x/b  -1.0 .Wheel, reducing the crack height near the crack mouth increases pressure in central crack segment and due to induced outward faces motion near the crack tip forces the dramatic pressure drop.Due to the viscosity of fluid the inflow of liquid to the low-pressure area is restrained and the pressure difference increases.When the cylinder is situated at the point x/b = -0.593 the pressure at the crack front reaches its minimum of -449 MPa.Starting from that point the pressure rises along the whole crack length and the growth is most rapid at the crack front, where it takes positive values as the cylinder comes as close as x/b = -0.556.The crack height reduces continuously at the crack mouth starting at the moment at which the cylinder comes into contact with the wedge-shaped part of raceway behind the crack (x/b = -1.146)and at the same time the flow rate of the liquid flowing out of the crack reduces, despite the pressure growth in the crack interior (near the crack front).The crack height reduces until the cylinder reaches the point x/b = 0.44, at which its value comes down to 0.32 μm at a distance of 0.33 mm from the crack mouth.At the point x/b = 0.701 the liquid pressure attains its maximal value equal to 907.8 MPa.As the cylinder rolls on the liquid pressure drops and the process is most rapid from the crack mouth side.It is caused by the process of crack mouth opening as the cylinder-prism contact zone rolls on behind the crack.When the cylinder crosses the point x/b  1.5 the opening rate of the crack mouth is high enough to enable generation of low negative pressure and the liquid flow into the crack (positive Q value at the crack mouth) while, at the crack front the faces come closer to each other and the pressure reaches the value of 700 MPa.The crack mouth opens until the cylinder-prism contact zone moves behind the crack (x/b  1.8) and then the crack height reduces.As the cylinder comes far away from the crack the pressure at the crack front drops, the velocities at which the crack faces come closer to each other reduce and the flows in the crack interior reduce as well.The analysis terminates at the cylinder position x/b = 5.224.During the rolling process the crack faces do not come into contact at any point.