- Email: [email protected]
rail contact systems
Wear 253 (2002) 35–41
Boundary element analysis of surface initiated rolling contact fatigue cracks in wheel/rail contact systems Makoto Akama a,∗ , Tadao Mori b,1 a
Railway Technical Research Institute, 2-8-38, Hikari-cho, Kokubunji-shi, Tokyo 185-8540, Japan b JR Souken Information System, 2-8-38, Hikari-cho, Kokubunji-shi, Tokyo 185-0034, Japan
Abstract This study performed the boundary element analysis of the cycle of stress intensity factors at the surface initiated rolling contact fatigue crack tip under Hertzian contact stress including an accurate model of friction between the faces of the crack and the effect of fluid inside the crack. A two-dimensional model of a rolling contact fatigue crack has been developed. The model includes the effect of Coulomb friction between the faces of the crack. The fluid in the crack was assumed not only to lubricate the crack faces and reduce the crack face friction coefficient but also to generate a pressure. The obtained results are compared with those of other researchers showing good consistency. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Boundary element method; Rolling contact fatigue; Stress intensity factor
1. Introduction Rolling contact fatigue cracks are one of the most important problems for the railway industry. Since Way [1] observed that lubrication was a key factor for the crack to propagate, many researchers proposed theoretical stress intensity cycles for the loading history [2–4]. Two possible effects of fluid have been proposed; hydraulic pressure mechanism [2] and fluid entrapment mechanism [3]. However, the studies by using these models did not seem to have confirmed yet. Boundary element method (BEM) is a well-established numerical technique in engineering fields and particularly suits to the problem of severe stress or strain gradients as there is no discretization of the interior data. BEM is also suitable for contact problems that have to be solved by iteration with incremental technique. When sliding occurs in the contact zone, the normal and tangential tractions have to be coupled with a friction parameter. Because the tractions are directly treated as nodal values in BEM, the friction conditions are easily treated. An extension to three-dimensional problems, including residual stresses or thermal stresses, can be easily performed. ∗ Corresponding author. Tel.: +81-42-573-7293; fax: +81-42-573-7289. E-mail addresses: [email protected] (M. Akama), [email protected] (T. Mori). 1 Tel.: +81-42-573-7682; fax: +81-42-580-7037.
This paper describes the development of a two-dimensional model of a surface initiated rolling contact fatigue crack under Hertzian contact by means of the BEM analysis and studies the change of stress intensity factors (K-values) at the crack tip. To be able to take into account residual and thermal stresses around the crack tip, boundary integral formulation including these stresses is also presented.
2. Boundary element analysis 2.1. Basic equations The boundary integral equation which permits to solve the general boundary value problem of elastostatics is as follows. i i clk uk + tlk∗ uk dΓ = u∗lk tk dΓ (1) Γ
Γ
where Γ is the boundary surface; uk and tk the k-components i a constant of displacement and traction, respectively; clk i = δ /2 that is decided by the shape of boundary and clk lk (δ lk is the Kronecker delta) when Γ is smooth at i; and uik represents the k-component of the displacement at the load application point i. The u∗lk and tlk∗ are k-components of displacement and traction due to a unit point load in the l direction, that are called fundamental solutions introduced by Kelvin.
0043-1648/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 3 - 1 6 4 8 ( 0 2 ) 0 0 0 8 0 - 7
36
M. Akama, T. Mori / Wear 253 (2002) 35–41
If a residual stress state exists, the elastic stresses are given by σij = σijt + σijr
(2)
where σijt represents the total stress and σijr the residual one. Taking into consideration the fundamental solution unit forces acting in the l direction, the following expression can be obtained from Eq. (1). i i ∗ ∗ ∗ clk uk + tlk uk dΓ = ulk tk dΓ − σijr εljk dΩ (3) Γ
Γ
Ω
∗ = ∂u∗ /∂x ; Ω is the domain on which the equawhere εljk j lk tion applies and is assumed to be bounded by. Thermoelastic problems may also be studied considering the temperature changes Θ in an elastic body as a body force equal to −γ Θ ,k at each point and increasing the tractions by γ Θ nk where
γ =
2µα(1 + ν) (1 − 2ν)
(4)
where α is the coefficient of thermal expansion; µ is the shear modulus; v is the Poisson’s ration; and nk is the direction cosine of the normal with respect to k-direction. Taking into account that for steady-state conduction, the boundary integral is i i clk uk + tlk∗ uk dΓ Γ (1 − 2ν)γ = u∗lk tk dΓ + u∗lk (Glk,kj Θnj 2(1 − ν) Γ Γ −Glk,k Θ,j nj ) dΓ (5) where Glk is the Galerkin vector for
u∗lk .
of the half-space is travelled by a Hertzian pressure distribution, with a contact half-width a and maximum pressure p0 . (6) P (x) = p0 1 − x 2 /a 2 The load moves over the surface from left to right. The contact area is also subjected to the distribution of tangential traction, with the maximum value q0 . The traction acts in the opposite direction of motion of the load. q(x) = f p(x)
(7)
where f is a frictional coefficient. The ratio q0 /p0 is set at −0.05. And the ratio of the length of the crack to the contact half-width was chosen as c/a = 0.5. Those were referred to in the work done by Bower [3]. 2.3. Modeling of elastic half-space The elastic half-space is represented by a finite sized mesh by appropriately displacing its boundaries. The boundary element mesh using multi-region modeling is shown in Fig. 2. The displaced portion of the mesh has a shape of half circle to avoid special treatment of corner nodes because tractions along the boundaries are unknown. There are 230 elements and 230 nodes. For the purpose of calculating displacements at the boundaries, the Hertzian pressure distribution is represented by the following equivalent concentrated forces. p0 π a P = (8) 2 Q=fP
(9)
The displacements on the boundary points are calculated for the concentrated forces using analytical solutions called
2.2. Coordinate system and representative input values The model used here is shown in Fig. 1. A two-dimensional inclined surface crack of length c is in an elastic half-space and the inclined angle β is 25◦ to the surface. The surface
Fig. 1. Analytical model, coordinate system and notation.
Fig. 2. Boundary element mesh.
M. Akama, T. Mori / Wear 253 (2002) 35–41
37
the ‘simple radial distribution’ by Flamant [5] taking into account the restrained following conditions. (a) Points on the z-axis have no lateral displacements when its symmetry is considered. (b) Let us assume that a point on the z-axis at a distance D from the origin does not move vertically. (c) There is no rigid body rotation of the solid at infinity. As the value of D, a distance between the loaded point and a point at the bottom of the mesh can be taken. The displacements due to P are ur =
2(1 − ν 2 ) P cos θ (ln D − ln r) πE (1 − 2ν)(1 + ν) − P θ sin θ πE
ur =
Fig. 3. Displacements due to concentrated forces P and Q.
2(1 − ν 2 )
P sin θ (ln r − ln D) πE (1 − 2ν)(1 + ν) + P (sin θ − θ cos θ ) πE 2ν(1 + ν) + P sin θ πE and those due to Q are uθ =
(10)
(11)
(1 − ν 2 ) (1 − 2ν)(1 + ν) 2Q cos θ ln r − Qθ sin θ πE πE (1 − 2ν)(1 + ν) + Q sin θ 2E (1 + ν) (1 − ν 2 ) 2 ln D + Q cos θ (12) + πE πE
(1 − ν 2 ) ν(1 + ν) 2Qsin θ ln r + Qθ sin θ uθ = πE πE (1 − 2ν)(1 + ν) (1 − 2ν)(1+ν) − Qθ cos θ + Q sin θ πE πE (13)
where r and θ are the polar coordinates centered at the point of application of the concentrated forces shown in Fig. 3; and E is the modulus of elasticity. These are given to the boundaries as forced displacements shown in Fig. 3. Fig. 4 shows a comparison of stress distributions under a concentrated normal force using the previous modeling without a crack with those predicted by the analytical solution under a Hertzian pressure distribution [6]. This figure indicates that the BEM mesh with displaced boundaries is an adequate representation of an elastic half-space. 2.4. Modeling of movement of Hertzian pressure distribution Consider a case where the Hertzian distribution moves one element to the right as shown in Fig. 5. To simulate the movement, an incremental distribution shown in (b) is obtained by subtracting the distribution before movement (a) from that after movement (c). Then a calculation is performed for the incremental distributed load to obtain all
Fig. 4. Stress distributions in an elastic half-space at depth z/a = 0.5 (P is applied to the center of the mesh).
38
M. Akama, T. Mori / Wear 253 (2002) 35–41
Fig. 6. Definition of tractions, displacements and normal gap between the crack faces.
2.6. Analysis of crack face contact The faces of the crack are divided into three parts: separated, locked and slipped regions. The tractions acting on the crack and displacements are defined as shown in Fig. 6. If the values are described in increments, the following relationships hold for each region. In the separated region, when tζ A0 + )tζoA > 0, tζ A0 + )tζ A = tζ B0 + )tζ B = pw0 + )pw , tξ A0 + )tξ A = tξ B0 + )tξ B = 0 Fig. 5. Modeling of movement of Hertzian pressure distribution.
displacements and tractions at nodes along the boundaries of the mesh. When those are superimposed on the accumulated values before movement, the displacements and tractions after movement are obtained. If the process is repeated for the number of elements that correspond to the distance of total movement, the simulation of translation completes. According to this modeling, the forced displacements that are given to the boundaries of mesh are also incremental. 2.5. The effect of fluid penetration Fluid may be trapped inside the crack when the mouth of the crack is sealed by the contact area. The fluid is assumed to be incompressible so that the pressure is determined by the condition that the volume of the crack remains constant. As the load moves over the surface, the fluid is forced toward the tip of the crack, and, subsequently, the crack begins to re-open until the mouth of the crack re-opens. The mouth is assumed to re-open at the instant when the condition that the volume of the crack can not remain constant by maintaining the element closest to the crack mouth contact. After that, the volume is not necessarily be constant and the pressure of the fluid inside the crack is assumed to be the same as that of the component of the contact pressure in the direction of the crack at the mouth.
(14)
where superscript o indicates the previous result of iterative calculation and subscript 0 is the convergent solution before the current iterative calculation is performed. Subscripts A and B denote corresponding elements on the two crack composing surfaces. The pw is the fluid pressure when the mouth of the crack is sealed. The ) denotes the rate of change with time. In the locked region, when tζ A0 + )tζoA ≤ 0, (15) tξ A0 + )tξoA |≤ µc | tζ A0 + )tζoA tζ A0 + )tζ A = tζ B0 + )tζ B , tξ A0 + )tξ A = tξ B0 + )tξ B
(16)
uζ A − uζ B − δp = 0,
(17)
uξ A − uξ B = 0
where µc is the friction coefficient of the crack face. The δ p is the normal gap between the faces of the crack. In the slipped region, when tζ A0 + )tζoA ≤ 0, (18) tξ A0 + )tξoA |= µc | tζ A0 + )tζoA tζ A0 + )tζ A = tζ B0 + )tζ B , tξ A0 + )tξ A = tξ B0 + )tξ B
(19)
uζ A − uζ B − δp = 0
(20)
The possibility of frictional locking between the crack faces makes the shape of the crack history dependent, so the load must be applied in stages. Therefore, the distributed load that
M. Akama, T. Mori / Wear 253 (2002) 35–41
is defined in Section 2.4 is applied in steps, while finding the locked and slipping parts of the crack using an iterative procedure. The location of the locked and slipping regions are guessed as follows. (a) With the crack faces assumed to be free, locations of node along the crack faces are analyzed when the incremental load is applied. Then the nodes that go into the body are searched based on the accumulated displacements. (b) If such regions exist, they are assumed to be locking regions and the BEM analysis is performed again. (c) If the assumed region is larger than the correct region, tensile stresses may be generated at some portions in the region. If there are such regions, they are excluded from the locking region and the same analysis is performed again. (d) If some portions of the locking region do not satisfy the condition expressed by the Eq. (15), these portions are assumed to be the slipping ones and tangential tractions that are calculated by the Eq. (18) are given as the boundary conditions. Normal tractions that are used in the calculation are the values of previous iteration. (e) If the mouth of the crack is sealed by the contact area, the fluid pressure is applied on the separated face of the crack so that the volume of the crack remains constant. When this condition is not met, the steps from (a) to (d) are performed again by changing the pressure. This process is repeated until convergence occurs. If convergence does not occur, the incremental load is halved and steps are repeated again. The calculational flow of the program is shown in Fig. 7.
3. Results The BEM model outlined in Section 2 has been used to calculate K-values for the crack tip. In this study, K-values are evaluated from the opening displacement at the seven nodes closest to the crack tip. By a least-square regression analysis of these data and extrapolation to the crack tip, K-values, the tensile opening mode KI (mode I) and the in-plane shear mode KII (mode II) are found by comparison of the values with known solutions in terms of K-values. G 2π KI = (21) lim )uζ 2(1 − v) r→0 r 2π G KII = (22) lim )uξ 2(1 − v) r→0 r where G is the modulus of shearing elasticity; r is the distance from the crack tip; )uζ and )uξ are relative displacements in each direction between the corresponding nodes on the crack faces. First, it is assumed that no fluid pressure acts on the faces of the crack but the faces of the crack are assumed to be well lubricated so that the crack face friction coefficient is 0.1.
Fig. 7. Calculational flow in an incremental step.
39
40
M. Akama, T. Mori / Wear 253 (2002) 35–41
Fig. 8. Variation of mode II stress intensities vs. the position of the load for a lubricated unpressurised crack.
Fig. 9. Variation of modes I and II stress intensities vs. the position of the load.
M. Akama, T. Mori / Wear 253 (2002) 35–41
The BEM model outlined in Section 2 has been used to calculate K-values for the surface breaking inclined crack. The variation of KII at the crack tip caused by the movement of the load over the half-space is shown in Fig. 8. The KII -value has been ploffed in a non-dimensional form as KII /p0√a . The results are shown as a function of the distance of the load from the mouth of the crack, x/a. The cycle may be regarded as a variation of the KII -value at the tip of crack as a function of time. The load rolls over the surface of the half-space in the positive x direction, so that x/a initially has a negative value and increases during the cycle. These are also compared with those by Bower [3] to show a good agreement. Next, the stress intensities generated by the process of fluid entrapment is considered. Fig. 9 shows the cycle of stress intensities when it is assumed that the fluid is trapped inside the crack when the mouth of the crack is sealed by the contact surface. The variation of KI and KII at the crack tip caused by the movement of the load over the half-space is shown in Fig. 9 in comparison with those by Kaneta et al. [2] and Bower [3]. Note that the results of Kaneta et al. are for a three-dimensional semi-circle crack and c/a = 0.5, q0 /p0 = −0.1, β = 45◦ without the effect of contact between the faces of crack. Again, the results are very similar to each other; the modes I and II stress intensity cycles overlapped and the mode II stress intensity is fully reversed. The relatively large differences between the results of Kaneta et al. and those of others are thought to have resulted from the fact that the former is for a three-dimensional semi-circle crack and without the effect of contact between the faces of crack.
41
4. Conclusions BEM analysis of the cycle of K-values was performed at the inclined surface crack tip under Hertzian contact stress by including an accurate model of friction between the faces and the effect of fluid trapped inside the crack. The resulting cycles of modes I and II stress intensities were calculated as the load rolls over the crack. The results were compared with those by other researchers, showing good agreement. Future work should include further theoretical fracture mechanics studies on the cracks including the effect of residual and thermal stresses. In the BEM analysis, these will be performed easily by using Eq. (3) or (5) instead of (1). References [1] S. Way, Pitting due to rolling contact, ASME J. Appl. Mech. 2 (1935) A49–A58. [2] M. Kaneta, H. Yatsuzuka, Y. Murakami, Mechanism of crack growth in lubricated rolling/sliding contact, ASLE Trans. 28 (3) (1985) 407– 414. [3] A.F. Bower, The influence of crack face friction and trapped fluid on surface initiated rolling contact fatigue cracks, transactions of the ASME, J. Tribol. 110 (1988) 704–711. [4] S. Bogdanski, M. Olzak, J. Stupnicki, Numerical modelling of 3D rail RCF squat type crack under operating load, Int. J. Fatigue Fracture Eng. Mat. Struct. 21 (1998) 923–935. [5] S.P.Timoshenko, J.N.Goodier, Theory of Elasticity, 3rd Edition, McGraw-Hill, 1985, pp. 97–104. [6] E. McEwen, Stresses in elastic cylinders in contact along a generatix, Phil. Mag. 40 (1949) 454–459.
Boundary element analysis of surface initiated rolling contact fatigue cracks in wheel/rail contact systems Makoto Akama a,∗ , Tadao Mori b,1 a
Railway Technical Research Institute, 2-8-38, Hikari-cho, Kokubunji-shi, Tokyo 185-8540, Japan b JR Souken Information System, 2-8-38, Hikari-cho, Kokubunji-shi, Tokyo 185-0034, Japan
Abstract This study performed the boundary element analysis of the cycle of stress intensity factors at the surface initiated rolling contact fatigue crack tip under Hertzian contact stress including an accurate model of friction between the faces of the crack and the effect of fluid inside the crack. A two-dimensional model of a rolling contact fatigue crack has been developed. The model includes the effect of Coulomb friction between the faces of the crack. The fluid in the crack was assumed not only to lubricate the crack faces and reduce the crack face friction coefficient but also to generate a pressure. The obtained results are compared with those of other researchers showing good consistency. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Boundary element method; Rolling contact fatigue; Stress intensity factor
1. Introduction Rolling contact fatigue cracks are one of the most important problems for the railway industry. Since Way [1] observed that lubrication was a key factor for the crack to propagate, many researchers proposed theoretical stress intensity cycles for the loading history [2–4]. Two possible effects of fluid have been proposed; hydraulic pressure mechanism [2] and fluid entrapment mechanism [3]. However, the studies by using these models did not seem to have confirmed yet. Boundary element method (BEM) is a well-established numerical technique in engineering fields and particularly suits to the problem of severe stress or strain gradients as there is no discretization of the interior data. BEM is also suitable for contact problems that have to be solved by iteration with incremental technique. When sliding occurs in the contact zone, the normal and tangential tractions have to be coupled with a friction parameter. Because the tractions are directly treated as nodal values in BEM, the friction conditions are easily treated. An extension to three-dimensional problems, including residual stresses or thermal stresses, can be easily performed. ∗ Corresponding author. Tel.: +81-42-573-7293; fax: +81-42-573-7289. E-mail addresses: [email protected] (M. Akama), [email protected] (T. Mori). 1 Tel.: +81-42-573-7682; fax: +81-42-580-7037.
This paper describes the development of a two-dimensional model of a surface initiated rolling contact fatigue crack under Hertzian contact by means of the BEM analysis and studies the change of stress intensity factors (K-values) at the crack tip. To be able to take into account residual and thermal stresses around the crack tip, boundary integral formulation including these stresses is also presented.
2. Boundary element analysis 2.1. Basic equations The boundary integral equation which permits to solve the general boundary value problem of elastostatics is as follows. i i clk uk + tlk∗ uk dΓ = u∗lk tk dΓ (1) Γ
Γ
where Γ is the boundary surface; uk and tk the k-components i a constant of displacement and traction, respectively; clk i = δ /2 that is decided by the shape of boundary and clk lk (δ lk is the Kronecker delta) when Γ is smooth at i; and uik represents the k-component of the displacement at the load application point i. The u∗lk and tlk∗ are k-components of displacement and traction due to a unit point load in the l direction, that are called fundamental solutions introduced by Kelvin.
0043-1648/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 3 - 1 6 4 8 ( 0 2 ) 0 0 0 8 0 - 7
36
M. Akama, T. Mori / Wear 253 (2002) 35–41
If a residual stress state exists, the elastic stresses are given by σij = σijt + σijr
(2)
where σijt represents the total stress and σijr the residual one. Taking into consideration the fundamental solution unit forces acting in the l direction, the following expression can be obtained from Eq. (1). i i ∗ ∗ ∗ clk uk + tlk uk dΓ = ulk tk dΓ − σijr εljk dΩ (3) Γ
Γ
Ω
∗ = ∂u∗ /∂x ; Ω is the domain on which the equawhere εljk j lk tion applies and is assumed to be bounded by. Thermoelastic problems may also be studied considering the temperature changes Θ in an elastic body as a body force equal to −γ Θ ,k at each point and increasing the tractions by γ Θ nk where
γ =
2µα(1 + ν) (1 − 2ν)
(4)
where α is the coefficient of thermal expansion; µ is the shear modulus; v is the Poisson’s ration; and nk is the direction cosine of the normal with respect to k-direction. Taking into account that for steady-state conduction, the boundary integral is i i clk uk + tlk∗ uk dΓ Γ (1 − 2ν)γ = u∗lk tk dΓ + u∗lk (Glk,kj Θnj 2(1 − ν) Γ Γ −Glk,k Θ,j nj ) dΓ (5) where Glk is the Galerkin vector for
u∗lk .
of the half-space is travelled by a Hertzian pressure distribution, with a contact half-width a and maximum pressure p0 . (6) P (x) = p0 1 − x 2 /a 2 The load moves over the surface from left to right. The contact area is also subjected to the distribution of tangential traction, with the maximum value q0 . The traction acts in the opposite direction of motion of the load. q(x) = f p(x)
(7)
where f is a frictional coefficient. The ratio q0 /p0 is set at −0.05. And the ratio of the length of the crack to the contact half-width was chosen as c/a = 0.5. Those were referred to in the work done by Bower [3]. 2.3. Modeling of elastic half-space The elastic half-space is represented by a finite sized mesh by appropriately displacing its boundaries. The boundary element mesh using multi-region modeling is shown in Fig. 2. The displaced portion of the mesh has a shape of half circle to avoid special treatment of corner nodes because tractions along the boundaries are unknown. There are 230 elements and 230 nodes. For the purpose of calculating displacements at the boundaries, the Hertzian pressure distribution is represented by the following equivalent concentrated forces. p0 π a P = (8) 2 Q=fP
(9)
The displacements on the boundary points are calculated for the concentrated forces using analytical solutions called
2.2. Coordinate system and representative input values The model used here is shown in Fig. 1. A two-dimensional inclined surface crack of length c is in an elastic half-space and the inclined angle β is 25◦ to the surface. The surface
Fig. 1. Analytical model, coordinate system and notation.
Fig. 2. Boundary element mesh.
M. Akama, T. Mori / Wear 253 (2002) 35–41
37
the ‘simple radial distribution’ by Flamant [5] taking into account the restrained following conditions. (a) Points on the z-axis have no lateral displacements when its symmetry is considered. (b) Let us assume that a point on the z-axis at a distance D from the origin does not move vertically. (c) There is no rigid body rotation of the solid at infinity. As the value of D, a distance between the loaded point and a point at the bottom of the mesh can be taken. The displacements due to P are ur =
2(1 − ν 2 ) P cos θ (ln D − ln r) πE (1 − 2ν)(1 + ν) − P θ sin θ πE
ur =
Fig. 3. Displacements due to concentrated forces P and Q.
2(1 − ν 2 )
P sin θ (ln r − ln D) πE (1 − 2ν)(1 + ν) + P (sin θ − θ cos θ ) πE 2ν(1 + ν) + P sin θ πE and those due to Q are uθ =
(10)
(11)
(1 − ν 2 ) (1 − 2ν)(1 + ν) 2Q cos θ ln r − Qθ sin θ πE πE (1 − 2ν)(1 + ν) + Q sin θ 2E (1 + ν) (1 − ν 2 ) 2 ln D + Q cos θ (12) + πE πE
(1 − ν 2 ) ν(1 + ν) 2Qsin θ ln r + Qθ sin θ uθ = πE πE (1 − 2ν)(1 + ν) (1 − 2ν)(1+ν) − Qθ cos θ + Q sin θ πE πE (13)
where r and θ are the polar coordinates centered at the point of application of the concentrated forces shown in Fig. 3; and E is the modulus of elasticity. These are given to the boundaries as forced displacements shown in Fig. 3. Fig. 4 shows a comparison of stress distributions under a concentrated normal force using the previous modeling without a crack with those predicted by the analytical solution under a Hertzian pressure distribution [6]. This figure indicates that the BEM mesh with displaced boundaries is an adequate representation of an elastic half-space. 2.4. Modeling of movement of Hertzian pressure distribution Consider a case where the Hertzian distribution moves one element to the right as shown in Fig. 5. To simulate the movement, an incremental distribution shown in (b) is obtained by subtracting the distribution before movement (a) from that after movement (c). Then a calculation is performed for the incremental distributed load to obtain all
Fig. 4. Stress distributions in an elastic half-space at depth z/a = 0.5 (P is applied to the center of the mesh).
38
M. Akama, T. Mori / Wear 253 (2002) 35–41
Fig. 6. Definition of tractions, displacements and normal gap between the crack faces.
2.6. Analysis of crack face contact The faces of the crack are divided into three parts: separated, locked and slipped regions. The tractions acting on the crack and displacements are defined as shown in Fig. 6. If the values are described in increments, the following relationships hold for each region. In the separated region, when tζ A0 + )tζoA > 0, tζ A0 + )tζ A = tζ B0 + )tζ B = pw0 + )pw , tξ A0 + )tξ A = tξ B0 + )tξ B = 0 Fig. 5. Modeling of movement of Hertzian pressure distribution.
displacements and tractions at nodes along the boundaries of the mesh. When those are superimposed on the accumulated values before movement, the displacements and tractions after movement are obtained. If the process is repeated for the number of elements that correspond to the distance of total movement, the simulation of translation completes. According to this modeling, the forced displacements that are given to the boundaries of mesh are also incremental. 2.5. The effect of fluid penetration Fluid may be trapped inside the crack when the mouth of the crack is sealed by the contact area. The fluid is assumed to be incompressible so that the pressure is determined by the condition that the volume of the crack remains constant. As the load moves over the surface, the fluid is forced toward the tip of the crack, and, subsequently, the crack begins to re-open until the mouth of the crack re-opens. The mouth is assumed to re-open at the instant when the condition that the volume of the crack can not remain constant by maintaining the element closest to the crack mouth contact. After that, the volume is not necessarily be constant and the pressure of the fluid inside the crack is assumed to be the same as that of the component of the contact pressure in the direction of the crack at the mouth.
(14)
where superscript o indicates the previous result of iterative calculation and subscript 0 is the convergent solution before the current iterative calculation is performed. Subscripts A and B denote corresponding elements on the two crack composing surfaces. The pw is the fluid pressure when the mouth of the crack is sealed. The ) denotes the rate of change with time. In the locked region, when tζ A0 + )tζoA ≤ 0, (15) tξ A0 + )tξoA |≤ µc | tζ A0 + )tζoA tζ A0 + )tζ A = tζ B0 + )tζ B , tξ A0 + )tξ A = tξ B0 + )tξ B
(16)
uζ A − uζ B − δp = 0,
(17)
uξ A − uξ B = 0
where µc is the friction coefficient of the crack face. The δ p is the normal gap between the faces of the crack. In the slipped region, when tζ A0 + )tζoA ≤ 0, (18) tξ A0 + )tξoA |= µc | tζ A0 + )tζoA tζ A0 + )tζ A = tζ B0 + )tζ B , tξ A0 + )tξ A = tξ B0 + )tξ B
(19)
uζ A − uζ B − δp = 0
(20)
The possibility of frictional locking between the crack faces makes the shape of the crack history dependent, so the load must be applied in stages. Therefore, the distributed load that
M. Akama, T. Mori / Wear 253 (2002) 35–41
is defined in Section 2.4 is applied in steps, while finding the locked and slipping parts of the crack using an iterative procedure. The location of the locked and slipping regions are guessed as follows. (a) With the crack faces assumed to be free, locations of node along the crack faces are analyzed when the incremental load is applied. Then the nodes that go into the body are searched based on the accumulated displacements. (b) If such regions exist, they are assumed to be locking regions and the BEM analysis is performed again. (c) If the assumed region is larger than the correct region, tensile stresses may be generated at some portions in the region. If there are such regions, they are excluded from the locking region and the same analysis is performed again. (d) If some portions of the locking region do not satisfy the condition expressed by the Eq. (15), these portions are assumed to be the slipping ones and tangential tractions that are calculated by the Eq. (18) are given as the boundary conditions. Normal tractions that are used in the calculation are the values of previous iteration. (e) If the mouth of the crack is sealed by the contact area, the fluid pressure is applied on the separated face of the crack so that the volume of the crack remains constant. When this condition is not met, the steps from (a) to (d) are performed again by changing the pressure. This process is repeated until convergence occurs. If convergence does not occur, the incremental load is halved and steps are repeated again. The calculational flow of the program is shown in Fig. 7.
3. Results The BEM model outlined in Section 2 has been used to calculate K-values for the crack tip. In this study, K-values are evaluated from the opening displacement at the seven nodes closest to the crack tip. By a least-square regression analysis of these data and extrapolation to the crack tip, K-values, the tensile opening mode KI (mode I) and the in-plane shear mode KII (mode II) are found by comparison of the values with known solutions in terms of K-values. G 2π KI = (21) lim )uζ 2(1 − v) r→0 r 2π G KII = (22) lim )uξ 2(1 − v) r→0 r where G is the modulus of shearing elasticity; r is the distance from the crack tip; )uζ and )uξ are relative displacements in each direction between the corresponding nodes on the crack faces. First, it is assumed that no fluid pressure acts on the faces of the crack but the faces of the crack are assumed to be well lubricated so that the crack face friction coefficient is 0.1.
Fig. 7. Calculational flow in an incremental step.
39
40
M. Akama, T. Mori / Wear 253 (2002) 35–41
Fig. 8. Variation of mode II stress intensities vs. the position of the load for a lubricated unpressurised crack.
Fig. 9. Variation of modes I and II stress intensities vs. the position of the load.
M. Akama, T. Mori / Wear 253 (2002) 35–41
The BEM model outlined in Section 2 has been used to calculate K-values for the surface breaking inclined crack. The variation of KII at the crack tip caused by the movement of the load over the half-space is shown in Fig. 8. The KII -value has been ploffed in a non-dimensional form as KII /p0√a . The results are shown as a function of the distance of the load from the mouth of the crack, x/a. The cycle may be regarded as a variation of the KII -value at the tip of crack as a function of time. The load rolls over the surface of the half-space in the positive x direction, so that x/a initially has a negative value and increases during the cycle. These are also compared with those by Bower [3] to show a good agreement. Next, the stress intensities generated by the process of fluid entrapment is considered. Fig. 9 shows the cycle of stress intensities when it is assumed that the fluid is trapped inside the crack when the mouth of the crack is sealed by the contact surface. The variation of KI and KII at the crack tip caused by the movement of the load over the half-space is shown in Fig. 9 in comparison with those by Kaneta et al. [2] and Bower [3]. Note that the results of Kaneta et al. are for a three-dimensional semi-circle crack and c/a = 0.5, q0 /p0 = −0.1, β = 45◦ without the effect of contact between the faces of crack. Again, the results are very similar to each other; the modes I and II stress intensity cycles overlapped and the mode II stress intensity is fully reversed. The relatively large differences between the results of Kaneta et al. and those of others are thought to have resulted from the fact that the former is for a three-dimensional semi-circle crack and without the effect of contact between the faces of crack.
41
4. Conclusions BEM analysis of the cycle of K-values was performed at the inclined surface crack tip under Hertzian contact stress by including an accurate model of friction between the faces and the effect of fluid trapped inside the crack. The resulting cycles of modes I and II stress intensities were calculated as the load rolls over the crack. The results were compared with those by other researchers, showing good agreement. Future work should include further theoretical fracture mechanics studies on the cracks including the effect of residual and thermal stresses. In the BEM analysis, these will be performed easily by using Eq. (3) or (5) instead of (1). References [1] S. Way, Pitting due to rolling contact, ASME J. Appl. Mech. 2 (1935) A49–A58. [2] M. Kaneta, H. Yatsuzuka, Y. Murakami, Mechanism of crack growth in lubricated rolling/sliding contact, ASLE Trans. 28 (3) (1985) 407– 414. [3] A.F. Bower, The influence of crack face friction and trapped fluid on surface initiated rolling contact fatigue cracks, transactions of the ASME, J. Tribol. 110 (1988) 704–711. [4] S. Bogdanski, M. Olzak, J. Stupnicki, Numerical modelling of 3D rail RCF squat type crack under operating load, Int. J. Fatigue Fracture Eng. Mat. Struct. 21 (1998) 923–935. [5] S.P.Timoshenko, J.N.Goodier, Theory of Elasticity, 3rd Edition, McGraw-Hill, 1985, pp. 97–104. [6] E. McEwen, Stresses in elastic cylinders in contact along a generatix, Phil. Mag. 40 (1949) 454–459.