Fatigue surface crack growth in cylindrical specimen under combined loading

Engineering Fracture Mechanics 131 (2014) 439–453

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Fatigue surface crack growth in cylindrical specimen under combined loading R. Citarella a,⇑, M. Lepore a, V. Shlyannikov b, R. Yarullin b a b

Dept. of Industrial Engineering, University of Salerno, via Giovanni Paolo II, 132, Fisciano, SA, Italy Research Center for Power Engineering Problems of Russian Academy of Sciences, Lobachevsky Street, 2/31, 420111 Kazan, Russia

a r t i c l e

i n f o

Article history: Received 8 February 2014 Received in revised form 25 August 2014 Accepted 30 August 2014 Available online 15 September 2014 Keywords: Surface flaw Tension and torsion Crack growth rate Crack path DBEM

a b s t r a c t The subject for studies is a steel bar of circular cross-section with straight-fronted edge notch undergoing fatigue loads. Both the optical microscope measurements and the crack opening displacement (COD) method are used to monitor and investigate both crack depth and crack length during the tests. The variation of crack growth behavior is studied under cyclic axial and combined tension + torsion fatigue loading. Results show that cyclic Mode III loading superimposed on the cyclic Mode I leads to a fatigue life reduction. In parallel to the experimental activity, numerical calculations are performed based on threedimensional DBEM analysis to determine the stress intensity factors along curvilinear surface crack front and fatigue life prediction. The experimental fatigue crack growth results obtained from round bar specimens have been compared with the numerical predictions. The computational DBEM results are found to be in satisfactory agreement with the experimental findings. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction In order to implement the operation with the constraint of a required fatigue life, it is necessary to carry out fracture mechanics assessment of a structural component under cyclic loading. The crack growth analysis of surface flaws is one of the most important parts for structural integrity prediction of the circular cylindrical metallic components (bars, wires, bolts, shafts, etc.) in the presence of initial and accumulated in-service damages. Fairly often, part-through flaws appear on the free surface of the cylinder and defects can be considered as semi-elliptical cracks. Multiaxial loading conditions including tension/compression, bending and torsion are typical for the circular cylindrical metallic components of engineering structures. The problem of residual fatigue life prediction of such type of structural elements is complex and analytical solutions are often not available because surface flaws are three-dimensional in nature. In the present study, firstly, experimental results of fatigue crack growth are presented with reference to a crack starting from a straight-fronted edge notch in a cylindrical specimen under axial tension loading with or without superimposed cyclic torsion. The influence of different loading conditions on fatigue life of tested specimens is discussed. The relations between crack opening displacement and crack length on the free surface of specimens are determined and it is found that the position and change of the curvilinear crack fronts is dependent on the initial straight-fronted edge notch depth. Using the aforementioned experimental relations, the crack front shape and crack growth rate in the depth direction can be predicted. Secondly, crack growth numerical calculations are performed based on the Dual Boundary Element Method ⇑ Corresponding author. E-mail addresses: [email protected] (R. Citarella), [email protected] (V. Shlyannikov). http://dx.doi.org/10.1016/j.engfracmech.2014.08.017 0013-7944/Ó 2014 Elsevier Ltd. All rights reserved.

440

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

Nomenclature a crack length on the free surface b crack depth C, m paris equation constants D specimen diameter E Young’s modulus G shear modulus In numerical constant J J-integral K1, K2, K3 Mode I , II and mode III elastic stress intensity factors K⁄ stress intensity factor at crack growth rate da/dN = 107 m/cycle MP Mode mixity parameter n strain hardening exponent N number of cycles R cyclic stress ratio r,h polar coordinates rc characteristic distance from the crack tip S strain energy density factor u, v , w displacement components V elementary volume W strain energy density B crack angle h⁄ angle of crack propagation g biaxial nominal stress ratio m Poisson ratio re Von Mises equivalent stress r~ e dimensionless angular function for Von Mises equivalent stress r0 yield stress FEM finite element method DBEM Dual Boundary Element Method SED Strain Energy Density MSED Minimum Strain Energy Density MPS Maximum Principal Stress SIF stress intensity factor COD crack opening displacement

(DBEM) [1–3]. The simulations for the crack path assessment are based on the DBEM code BEASY [4]. DBEM is an extremely powerful approach for three dimensional fracture mechanics problems also in comparison with FEM [5–7] or XFEM [8,9], due to its intrinsic accuracy and flexibility (the three dimensional crack propagation proceeds in a fully automatic way), but the best results are nowadays obtained by a coupled usage of the two procedures [10–12]. The computational 3D fracture analyses deliver variable mixed mode conditions along the crack front [13,14]. Different criteria for the crack path assessment like the Minimum Strain Energy Density (MSED) [15], Maximum Principal Stress (MPS) [16] and Q-Plain [17]) are adopted. In the latter approach, the crack will propagate in a direction where the multiaxiality factor has a minimum value (the multiaxiality factor is defined as the ratio of the Von Mises equivalent stress to the hydrostatic stress). The stress intensity factor (SIF) evaluation is based on the J-integral approach [18,19].

2. Test material and specimens The test material is carbon steel R2 M. Its main mechanical properties, including constants of the Paris equation (Eq. (1)), are listed in Table 1 where E is the Young’s modulus, rb is the nominal ultimate tensile strength, r0 is the monotonic tensile yield strength, rf is the true ultimate tensile strength, n is the strain hardening exponent, C and m are the Paris constants which are determined using standard compact specimen.

Table 1 Main mechanical properties of steel R2M. E (MPa)

rb (MPa)

r0 (MPa)

rf (MPa)

n

C ((mm/cycle)/(MPa mm0.5)m)

m

209000

810.3

540.8

890

6.134

4.1  1013

2.818

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

da ¼ C DK m dN

441

ð1Þ

The specimen geometry configuration is shown in Fig. 1. The diameter is equal to 10 mm in the test section and the length is equal to 100 mm. Using a linear cutting machine surface edge notches were cut with initial depth b0 = 1.0 mm. The geometric parameters of the specimen test section and of the growing crack are shown in Fig. 2; in this figure, b is the current crack depth, with the crack front approximated by an elliptical curve with major axis 2c and minor axis 2b. The crack length a is obtained by measuring the distance between the advancing crack break through point and the notch break through point, as shown in Fig. 2a. The depth of the initial straight edge notch is denoted by h and the initial notch length by L. The crack opening displacement is measured on the free specimen cylindrical surface, in the central plane of symmetry as shown in Fig. 2b. The Axial Torsion Test System Bi-00-701 is used for axial–torsional fatigue and fracture testing of the round bar specimens. This system is equipped with: fatigue rated axial–torsional dynamic load cell with axial capacity 100 kN and torsional capacity 2 kN-m; Bi-06-3XX series axial extensometers and torsional strain measurement fixture. The crack length on the specimen lateral surface are monitored using the optical instrumental zoom microscope whereas, to correlate the crack opening displacement with the gauge length, a pulley arrangement with an externally axial encoder is introduced (Fig. 3). All tests are carried out at room temperature, with sinusoidal loading shape, under imposed traction and torque. For the simple cyclic tension fatigue tests, the specimens are tested with an applied maximum nominal stress equal to 250 MPa and with a stress ratio R = 0.1. The combined tension/torsion tests are performed with the same stress ratio, applying synchronous and in-phase tensile and shear stresses whose maximum values are respectively equal to 250 and 100 MPa. Two different frequency values (10 and 7 Hz) are used during the fatigue test in order to highlight the crack front geometry during propagation: during each test, beach marks are produced on each specimen by reducing the applied frequency from 10 to 7 Hz, when the surface crack length is approximately equal to a ffi 1 mm. The typical beach marks on the post mortem cross section of different specimens are shown in Fig. 4a and b. From the crack front shape obtained in such manner, the relations between the relative crack depth b/D and the surface crack chord length a/D can be measured using an optical zoom microscope. Furthermore, based on periodically measured increments of surface crack length Da, the curve of surface crack propagation versus cycle numbers can be determined. Finally, using the relation of crack depth as a function of surface crack length, it is possible to calculate the crack growth rates db/dN in the depth direction. 3. Strain energy density distribution In order to study the influence of the different loading conditions on material cyclic mixed mode fracture resistance characteristics, it is necessary to use some equivalent functions for stress intensity factors (SIFs): to this aim the Strain Energy Density (SED) is selected. The surface flaw in cylindrical specimen undergoing combined loading exhibits mixed mode fracture conditions, with all the three components of stress intensity factors, mode I (KI), mode II (KII) and mode III (KIII). The combined tension and torsion remote loading causes such mode mixity. In a steel bar of circular cross-section with straight-fronted edge notch (Fig. 4) SIFs change from point to point along the curvilinear crack front. It is of interest to obtain the crack kinking angle h* at two main points – the break through point and the deepest point of crack front. To this end, the general SED equation proposed by Sih [15] can be applied. For linear elastic material behavior the strain energy density function dW/dV can be written as

dW 1 ¼ ðrxx exx þ ryy eyy þ rzz ezz þ rxy cxy þ ryz cyz þ rzx czx Þ dV 2        1 @u @v @w @ v @u @w @ v @u @w þ ryz þ rzx rxx þ ryy þ rzz þ rxy þ þ þ ¼ 2 @x @z @y @z @z @x @y @x @y

ð2Þ

where rij, eij are the stress and strain components, u, v, w are displacement components. Consider the first terms of singular solution, i.e. terms proportional to r1/2 and use the following stress and displacement expansion in form of equations [16]

Fig. 1. Details of specimen geometry.

442

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

Fig. 2. (a and b) Crack geometric parameters.

Fig. 3. Test equipment for measuring crack length and COD.

Fig. 4. Photograph of the cross section of specimens (a-pure tension, b-tension + torsion).

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

K

h



h

 3h

K

h

K

h



h

 3h

K

h

K

h

K

h



h

3h

1 2 ffi cos  pffiffiffiffiffiffiffiffiffi sin rxx ¼ pffiffiffiffiffiffiffiffi 1  sin sin 2 þ cos cos 2 2 2 2 2 2 2pr 2pr

h



3h

1 2 ffi cos þ pffiffiffiffiffiffiffiffiffi sin cos cos 1 þ sin sin ryy ¼ pffiffiffiffiffiffiffiffi 2 2 2 2 2 2 2pr 2pr

h

3h

K

h



h

443

ð3Þ

 3h

1 2 ffi sin cos cos þ pffiffiffiffiffiffiffiffi ffi cos 1  sin sin rxy ¼ pffiffiffiffiffiffiffiffi 2 2 2 2 2 2 2pr 2pr

K

h

3 3 ffi cos ; rzx ¼  pffiffiffiffiffiffiffiffiffi sin ryz ¼ pffiffiffiffiffiffiffiffi 2 2 2pr 2pr

K1 G

rffiffiffiffiffiffiffi rffiffiffiffiffiffiffi       r h 1 h K2 r h 1 h 2 þ ðj  1Þ þ sin ðj þ 1Þ þ cos2 cos sin 2p 2 2 2 2 2 2 G 2p

K1 ¼ G

rffiffiffiffiffiffiffi rffiffiffiffiffiffiffi       r h 1 h K2 r h 1 h 2 2 þ ðj þ 1Þ  cos ð1  jÞ þ sin sin cos 2p 2 2 2 2 2 2 G 2p

u¼

v

w¼

2K 3 G

ð4Þ

rffiffiffiffiffiffiffi r h E : sin ; G ¼ 2p 2 2ð1 þ mÞ

Applying Eqs. (3) and (4) to Eq. (2), and performing the necessary algebra, the expansion of the strain energy density field is resulting for inclined crack subjected to mixed mode loading

dW ða11 K 2I þ a12 K I K II þ a22 K 2II þ a33 K 2III Þ ¼ dV 4Gpr dW S ¼ ; dV 4Gpr

ð5Þ

where S is the SED factor with the specified angular functions adopted for a surface flaw by authors [20,21]; actually in [20] the strain energy density expression includes also higher order terms (the purpose of such article was to assess the effect of non-singular terms of the elastic stress expansion on the direction of crack propagation) that are neglected here

S ¼ a11 ðhÞK 2I þ a12 ðhÞK I K II þ a22 ðhÞK 2II þ a33 ðhÞK 2III

ð6Þ

a11 ¼ F x1 ex1 þ F y1 ey1 þ F xy1 ðg 11 þ g 21 Þ; a22 ¼ F x2 ex2 þ F y2 ey2 þ F xy2 ðg 12 þ g 22 Þ a12 ¼ F x1 ex2  F x2 ex1 þ F y1 ey2 þ F y2 ey1 þ F xy1 ðg 12 þ g 22 Þ þ F xy2 ðg 11 þ g 21 Þ; a33 ¼ 1: F x1 ¼ cos

    h h 3h h h 3h F x2 ¼ sin 1  sin sin 2 þ cos cos 2 2 2 2 2 2

F y1 ¼ cos

  h h 3h h h 3h F y2 ¼ cos sin cos 1 þ sin sin 2 2 2 2 2 2

F xy1 ¼ cos

h h 3h sin cos 2 2 2

F xy2 ¼ cos

  h h 3h 1  sin sin 2 2 2

ex1 ¼

     1 1 h h h 1 1 h h h 2 h 2 h cos  sin sin h cos2 ; ex2 ¼  sin þ cos sin h sin ðj  1Þ þ sin ðj þ 1Þ þ cos2 2 2 2 2 2 2 2 2 2 2 2 2

ey1 ¼

     1 1 h h h 1 1 h h h 2 h 2 h cos þ cos cos h sin  sin þ sin cos h cos2 ðj þ 1Þ  cos2 ; ey2 ¼ ð1  jÞ þ sin 2 2 2 2 2 2 2 2 2 2 2 2

g 11 ¼

     1 1 h h h 1 1 h h h 2 h 2 h  sin  cos sin h sin cos  sin sin h cos2 ðj þ 1Þ  cos2 ; g 12 ¼ ð1  jÞ þ sin 2 2 2 2 2 2 2 2 2 2 2 2

g 21 ¼

     1 h h h 1 1 h h h 2 h 2 h  sin þ sin cos h cos2 ; g 22 ¼ cos  cos cos h sin ðj  1Þ þ sin ðj þ 1Þ þ cos2 : 2 2 2 2 2 2 2 2 2 2 2

j = (3  m)/(1 + m) – plane stress; j = 3  4m – plane strain; m is the Poisson ratio; r, h are polar coordinates centered at the crack tip and G = E/(2(1 + m)) the shear modulus.

444

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

The third hypothesis of SED-criterion is used to determine the angle of crack propagation h* (Eq. (7)). Thus, a line drawn from each point on the crack front, lying in the plane normal to the crack front, with an angle h* with respect to the crack plane, indicates the directions in which the strain energy density has its minimal value

 @W   @h 

r¼rC

 @ 2 W  ¼ 0;  @h2 

> 0; W ¼ dW=dV:

ð7Þ

r¼r C

where rc is a characteristic distance from the crack tip. The general hypotheses of the strain energy density theory are associated with a concept of a critical distance, a fundamental parameter that provides a relation between the processes occurring on both micro- and macro-level scales. When applying any of fracture criteria to predict crack propagation, the stress–strain characteristics are not determined at the crack tip itself, but at some distance rc from it. To take advantage of Eq. (7), it is necessary to define the sense of the radial distance rc contained in this equation. The most modern of the fracture mechanics approaches are based on a critical distance from the crack tip. In the present study the critical distance rc ahead of the crack tip, according to the SED-concept, is assumed to pinpoint a location where the stress strain state in the specified element reaches a certain critical value that can be determined from uniaxial test data. A relative fracture damage zone size (FDZ) or critical distance  dc ¼ rc =a was introduced by Shlyannikov [22]

h   i92 8
Static loading

W c ¼



r0 ryn

2   1 2 an  nþ1 r f þ rf ; 2 nþ1

ð8Þ

cyclic loading



ðWÞc ¼ 4rf ef ð2Nf Þm

In these equations, ef is the fatigue ductility, rf is the fatigue strength coefficient, n is the strain hardening exponent, DK th is the threshold value of stress intensity range, Nf is the fatigue life. In Eq. (8) Si ði ¼ 1; 2; 3Þ and Sp are elastic and plastic coef~ eÞ ficients respectively. The work [22] provides more details about the determination of Si ¼ Si ðh; j; b; gÞ; Sp ¼ Sp ðn; m; In ; Mp ; r and about the employment of the strain energy density theory in relation to the critical distance or fracture damage zone size for the solution of mixed mode problems. Finally, substituting into Eq. (6) the numerical values of SIFs [23] related to the considered loading conditions of cylindrical specimens, the values of h* at all points of the curvilinear crack front can be calculated, using Eq. (7), for different combinations of mode mixity and surface flaw geometry. With reference to the crack front position defined by the point A on the free surface of the cylindrical specimen (Fig. 2), the variations of the strain energy density factor are displayed (Fig. 5) as a function of the polar angle h measured from the line extending from the crack plane (Fig. 6). According to the strain energy density criterion, the direction of surface crack growth initiation h = h* is assumed to correspond to the minimum of (dW/dV). Using the SED-criterion, i.e. Eqs. (5) and (6), the crack extension angles h* at the break through point A, under combined tension and torsion loading, are plotted for the surface flaws in Fig. 6 as functions of dimensionless crack length. It should be noted that at the same loading conditions the crack kinking angle in the deepest point B of the crack front is equal to zero (h* = 0).

Fig. 5. SED factor variations calculated at break through point A under combined loading.

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

445

Fig. 6. Crack deviation angle at break through point A.

4. Experimental results The evolution, during the tests, of the crack growth rate of the elliptical crack front is determined using COD and the zoom optical microscope. In order to study the crack growth in cylindrical specimens under fatigue tension loading with superimposed cyclic torsion, two specimens are tested for each loading condition, with an initial notch depth equal to 1 mm. Fig. 7a and b shows plots of the break through point advances a and of COD against the number of cycles N under pure tension and tension–torsion combined loads, respectively. As shown in Fig. 7a and b, in-phase cyclic torsion loading superimposed on cyclic tension leads to decreasing of fatigue life. Fig. 8(a and b) shows the superficial crack growth rate da/dN versus crack length a and versus COD, under pure cyclic tension and combined fatigue loading. It can be seen that there is not a significant increase of the crack growth rates along the external surface direction when the cyclic Mode III loading is superimposed on the cyclic tension. However, looking at Fig. 7 and considering changes in the general durability of the specimens under pure tension rather than under combined loading, significant increase in the crack growth rate in the depth direction b under the latter type of loading condition are expected. Fig. 9 shows the relationship between superficial crack length and COD on the cylindrical specimens undergoing pure tension and combined loading. It is found that the crack growth rates along the external surface direction are similar for different loading conditions. On the base of this experimental data, polynomial functions can be used to express the COD as a function of the superficial crack length. The evolution of surface crack shape of the straight-fronted edge cracks during the tests is obtained using beach marks and optical instrumental microscope. The results are observed for different loading conditions. An example of the shape development is shown in Fig. 10 which coincide with trends mentioned in [24].

Fig. 7. (a) Fatigue crack growth and (b) COD curves on free surface of specimens versus cycles.

446

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

Fig. 8. Crack growth rate against superficial crack length a (a) and COD (b).

Fig. 9. Relationship between a size and COD.

In the well-known Paris equation, that can also be written as

da ¼ CK m 1;max dN

ð9Þ

there exists a strong correlation between the parameters C and m (correlation coefficient not less than 0.96)

log C ¼ am þ b; ða; b < 0Þ:

ð10Þ

Therefore, the parameters C and m cannot be considered independent constants to describe the material resistance to crack growth under cyclic loading. Some authors [25] have proposed to use a different equation to describe the rate of crack propagation

   m da da K 1;max ¼ ;  dN dN K1

ð11Þ

where the exponent m is the same as in formula (6); K* is a parameter with the dimension of stress intensity factor; da/dN* is a coefficient having the dimension of the crack growth rate. It is assumed that da/dN* = 107 m/cycle because such crack growth rate is always localized in the linear part of the fatigue failure diagram for most of the structural materials. There

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

447

Fig. 10. Crack growth shape under cyclic loading.

is no correlation between the parameters m and K* and therefore the parameter K*, along with the constant m is an independent characteristic of cyclic crack resistance of the material. The parameter K* has a clear physical meaning: it is the greatest stress intensity factor at a value of the crack growth rate da/dN* = 107 m/cycle. From the compatibility solutions of Eqs. (9) and (11) we can find the relationship between the parameters C and K*

 C¼

da dN



ðK 1 Þ

m

:

ð12Þ

In case of mixed mode fracture instead of the mode I stress intensity factor (KI) it is necessary to use a particular equivalent parameter. In this paper, such parameter is the strain energy density function (S)

S ¼ a11 ðhÞK 21 þ a12 ðhÞK 1 K 2 þ a22 ðhÞK 22 þ a33 ðhÞK 23 :

0

ð5 Þ

Then the crack growth rate Eqs. (6) and (8) can be rewritten as

da ¼ CSnmax dN

ð13Þ

   n da da Smax ¼ ; dN dN S

ð14Þ

where n = m/2. Table 2 shows the experimental values of the constants of Eqs. (9), (11), (13) and (14) for the standard compact specimens and cylindrical specimens, the latter tested under pure tensile loading or combined tension and torsion. There are obvious advantages when using Eqs. (11) and (14) to characterize the material resistance to cyclic crack growth. The quantitative values of parameters K* and S* (Table 2) are also evident from Fig. 11; in particular, the minimum value of S* is related to the minimum specimen durability. Thus, from Table 2 it follows that the minimum value of S* corresponds to the crack growth rates calculated at the deepest point of the crack front under combined tension and torsion loading. This conclusion is confirmed by the relative position of crack growth curves in Fig. 11. The results of Table 2 indicate that the constant C, m and K* do not coincide with each other for compact and cylindrical specimens tested under a simple uniaxial tension. Moreover, the values of the constants C, m and S* are different from each other for cases of pure tensile loading and combined loading of cylindrical specimens. In our case it is obvious that the minimum value of S* = 153 MPa2 m corresponds to the crack growth rate in the deepest point of the crack front under combined loading. For this variant we have the lowest durability of the specimen. The values of S* = 469–479 MPa2 m correspond to test specimen under pure tension, with the greatest durability. This durability series are shown in Fig. 11b. 5. DBEM model The crack propagation in the cylindrical specimens N.1 and N.2, undergoing a combined traction-torsion fatigue load, with the geometry of the notch plane section detailed in Fig. 12a and b, is simulated by the DBEM model under linear elastic hypothesis (Fig. 13).

448

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

Table 2 Constants of crack growth rate equations. Compact specimen

Round bar

Constants of Eq. (6) in Constants of Eqs. (10) and SIF-terms (11) in SEDF-terms

Crack front position

C Pure tension

m

6.706  10

Tension + torsion –

12

C

2.818 2.438  10 –

–

n 11

S*

Constants of Eq. (6) in Constants of Eqs. (10) and SIF-terms (11) in SEDF-terms C

m

C

12

n

S*

12

1.409 366.5 Free surface (point A) 2.028  10 3.2046 7.193  10 1.553 465.1 Deepest point (point b) 1.489  1012 3.127 7.114  1012 1.547 479.3 –

–

Free surface (point A) – Deepest point (point B) –

– –

2.745  1011 1.423 317.7 1.189  1010 1.337 153.4

Fig. 11. Crack growth rate in: compact specimens and round bars (a) and only round bars (b).

Fig. 12. (a and b) Geometry of specimens N. 1 (a) and 2 (b) (length expressed in mm).

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

449

Fig. 13. DBEM model with highlight of initial crack.

Beasy STEP 0

Beasy STEP 5

Beasy STEP 10 Fig. 14. Contour plot of Maximum Principal Stresses (MPa) at different stages of crack propagation, with external and internal views of the growing crack.

The geometric model and the related mesh are realized directly in BEASY preprocessing environment. The specimen has a length equal to 60 mm, with initial and final mesh (after 10 steps of crack propagation) respectively based on 2258 and 2649 quadratic elements, with one end clamped and the other end constrained along the in plane radial directions and loaded along the axial and tangential directions. The initial crack, as indicated by experimental measurements (beach mark technique) has the following sizes: a = 0.855 mm, b = 0.852 mm, c = 3.705 mm. During the propagation, the average crack advance is set equal to 0.2 mm. As previously said, the SIF’s along the crack front are calculated by the J-integral approach (the J-paths along the crack front are visible from Fig. 13). The crack growth rate is calculated with the Paris formula (Eq. (1)), whose calibration parameters are shown in Table 1. The crack path is calculated using different criteria like MSED, MPS and Q-plane, in such a way to provide a comparison between them. The comparison between the numerical and experimental crack shape is still only qualitative but in the future will be made quantitative at least for the break through points.

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

Crack length a

450

5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 60000

Experimental Specimen 1 Beasy Specimen 1

80000

Experimental Specimen 2 Beasy Specimen 2

100000

120000

140000

Cycles Number Fig. 15. Crack length [mm] versus cycles for the two analyzed specimens.

0.7

Step 0 - MSED Step 0 - MPS Step 0 - Qplain

Growth angle

0.5 0.3 0.1 -0.1 0.0

0.2

0.4

0.6

0.8

1.0

-0.3 -0.5 -0.7

Crack front normalised abscissa Fig. 16. Growth angle [rad], from the three different criteria, as varying along the first crack front.

0.3

Growth angle

0.2 0.1 0.0

0.0

0.2

0.4

0.6

0.8

1.0

-0.1

Step 1 - MSED Step 1 - MPS Step 1 - Qplain

-0.2 -0.3

Crack front normalised abscissa Fig. 17. Growth angle [rad], from the three different criteria, as varying along the second crack front.

6. DBEM numerical results and crack growth prediction In Fig. 14 it is possible to see the crack propagation using the Q-Plain approach, starting from the initial configuration (step 0) and proceeding through the step 8 up to the final step 15: the crack kinking induced by the mode III (coming from the torsion load) superposition is evident. In Fig. 15 it is possible to appreciate an acceptable level of correlation between experimental and numerical crack growth rates for both the analyzed specimens. The numerical simulation starts from the first recorded crack front, so that there is no comparison in the initial stage of experimental monitored crack growth. In Fig. 16 the growth angle from the initial cracked configuration is shown against its variation along the crack front and against the analyzed crack growth criteria: it is possible to see that the all criteria provide similar predictions. In the next step (Fig. 17) the prediction of MSED criterion seems to significantly deviate from the predictions of the other two criteria: actually such deviation started in the previous crack propagation step even if with a much lower intensity.

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

Experimental Specimen 3 Experimental Specimen 2

Crack length

5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 110000

130000

150000

Numerical Specimen 3 Numerical Specimen 2

170000

190000

210000

Cycles Number Fig. 18. Crack length a [mm] versus number of cycles for the two specimens undergoing pure traction load.

70 60 50

KI

40 30 20 10

incr 0

0.0

0.2

incr 2

0.4

incr 4

0.6

0.8

1.0

0.8

1.0

0.8

1.0

Normalized abscissa 10 8 6 4

KII

2 0

-2 0.0

0.2

0.4

0.6

-4 -6 -8

incr 0

-10

incr 2

incr 4

Normalized abscissa 6 4

K III

2 0

0.0

0.2

0.4

0.6

-2 -4 -6

incr 0 incr 6

incr 2 incr 8

incr 4 incr 10

Normalized abscissa Fig. 19. SIF’s [MPa mm0.5] along the crack front related to the MSED usage, for mode I (KI), mode II (KII) and mode III (KIII).

451

452

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

70 60

K eff

50 40 30 20 10

incr 0

0.0

0.2

incr 2

incr 4

0.4

incr 6

0.6

incr 8

0.8

incr 10

1.0

Normalized abscissa Fig. 20. Equivalent SIF’s [MPa mm0.5] along the crack front related to the MSED usage.

In Fig. 18 the numerical-experimental correlation is presented with reference to the two specimens loaded in pure traction (originally there were three of them but the N. 1 failed): in this case the level of correspondence is poor. A possible reason, could be the following: the crack growth rate was measured on the surface of a cylindrical specimen (conditions close to plane stress), but the Paris constants (C and m) are defined on a compact specimen in conditions close to plane strain. Combined 3D tension and torsion loading under mixed mode I + III is probably less sensitive to the differences between 2D and 3D than under pure mode I. At the beginning of crack growth, the growth rate of the central point of the crack front is faster than that at the intersection with the surface (break through points). The reason is that the maximum stress intensity factor (Figs. 19 and 20) is attained at the deepest point on the crack front for the straight-fronted edge crack (b/c = 0). Consequently, a straight-fronted crack tends to become curved, and the flaw aspect ratio b/c increases. Furthermore, under combined loading conditions the crack growth rate is higher with respect to pure tension in the deepest point of the curvilinear crack front. 7. Conclusions The computational DBEM results are found to be in good qualitative (the crack path) and quantitative (the crack growth rates) agreement with experimental findings, which show a rather complex 3D crack growth behavior in case of traction and torsion. The fatigue life is decreased when an in-phase cyclic torsion is added to the cyclic tension. This is related to the increase of the mode mixity effect. Another interesting result pointed out in the present paper is the relative sensitivity, when considering the analyzed multi-axial loading condition, of growth angles and effective SIFs against a variation of the tested propagation criteria among the three considered. Moreover, the DBEM approach allow a reductions of the preprocessing time because the crack insertion and the whole crack propagation is fully automatic, with repeated remeshing realized at each crack step without user intervention. Acknowledgment The work was partly supported by the Russian Scientific Foundation (Project No 14-19-01716). References [1] Mi Y, Aliabadi MH. Dual boundary element method for three dimensional fracture mechanics analysis. Engng Anal Boundary Element 1992;10:161–71. [2] Citarella R, Cricri G, Armentani E. Multiple crack propagation with dual boundary element method in stiffened and reinforced full scale aeronautic panels. Key Engng Mater 2013;560:129–55. [3] Armentani E, Citarella R, Sepe R. FML full scale aeronautic panel under multiaxial fatigue: experimental test and DBEM simulation. Engng Fract Mech 2011;78(8):1717–28. Special Issue on ‘Multiaxial Fracture’. [4] BEASY V10r14, Documentation, C.M. BEASY Ltd; 2011. [5] Citarella R, Cricri G. Comparison of DBEM and FEM crack path predictions in a notched shaft under torsion. Engng Fract Mech 2010;77:1730–49. [6] Citarella R, Buchholz F-G. Comparison of crack growth simulation by DBEM and FEM for SEN-specimens undergoing torsion or bending loading. Engng Fract Mech 2008;75:489–509. [7] Citarella R, Cricrì G, Lepore M, Perrella M. In: Öchsner A, da Silva LFM, Altenbach H, editors. Materials with complex behaviour – advanced structured materials DBEM and FEM analysis of an extrusion press fatigue failure, vol. 3 part 2. Berlin (Germany): Springer-Verlag; 2010. p. 181–91. [8] Goangseup Z, Jeong-Hoon S, Budyn E, Sang-Ho L, Belytschko T. A method for growing multiple cracks without remeshing and its application to fatigue crack growth. Model Simul Mater Sci Eng 2004;12:901–15. [9] Chopp DL, Sukumar N. Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite element/fast marching method. Int J Engng Sci 2003;41:845–69.

R. Citarella et al. / Engineering Fracture Mechanics 131 (2014) 439–453

453

[10] Citarella R, Cricrì G, Lepore M, Perrella M. Thermo-mechanical crack propagation in aircraft engine vane by coupled FEM-DBEM approach. Adv Engng Software 2014;67:57–69. [11] Citarella R, Cricrì G. Three-dimensional BEM and FEM submodelling in a cracked FML full scale aeronautic panel. Appl Compos Mater 2014. http:// dx.doi.org/10.1007/s10443-014-9384-5. [12] Citarella R, Lepore M, Fellinger J, Bykov V, Schauer F. Coupled FEM-DBEM method to assess crack growth in magnet system of Wendelstein 7-X. FratturaedIntegritàStrutturale 2013;26:92–103. [13] Qian J, Fatemi A. Mixed mode fatigue crack growth. A literature survey. Engng Fract Mech 1996;55:969–90. [14] Rozumek D, Macha E. A survey of failure criteria and parameters in mixed-mode fatigue crack growth. Mater Sci 2009;45(2):190–210. [15] Sih GC. Strain energy density factor applied to mixed mode crack problems. Int J Fract 1974;10:305–21. [16] Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse shear. J Basic Engng 1963;86:519–27. [17] Labeas G, Kermanidis Th. Stress multiaxiality factor for crack growth prediction using the strain energy density theory. Theoret Appl Fract Mech 2006;45:100–7. [18] Rigby RH, Aliabadi MH. Mixed-mode J-integral method for analysis of 3D fracture problems using BEM. Engng Anal Boundary Elem 1993;11:239–56. [19] Rigby RH, Aliabadi MH. Decomposition of the mixed-mode J-integral – revisited. Int J Solids Struct 1998;35(17):2073–99. [20] Shlyannikov VN, Kislova SY, Tumanov AV. Inclined semi-elliptical crack for predicting crack growth direction based on apparent stress intensity factors’’. Theoret Appl Fract Mech 2010;53:185–93. [21] Shlyannikov VN, Tumanov AV. An inclined surface crack subject to biaxial loading. Int J Solids Struct 2011;48:1778–90. [22] Shlyannikov VN. Modelling of crack growth by fracture damage zone. Theoret Appl Fract Mech 1996;25:187–201. [23] Thompson KD, Sheppard SD. Stress intensity factors in shafts subjected to torsion and axial loading. Engng Fract Mech 1992;42:1019–34. [24] Feng Peng Yang, Zhen Bang Kuang, Shlyannikov VN. Fatigue crack growth for straight-fronted edge crack in a round bar. Int J Fatigue 2006;28:431–7. [25] Panasiuk VV, Ya Yarema S. General features of fatigue fracture diagrams of metals. Strength Mater 1996;1:30–5.