Determination of the Paris law constants in round bars from beach marks on fracture surfaces

Polytechnic Institute of Coimbra, ISEC, Department of Mechanical Engineering, Coimbra, Portugal bCEMU cDepart the crack shape evolution and the number of fatigue cycles, is employed. Finally, a compar- cular cracks tend to be flat [1–8]; besides, the production of standard specimens from mechanical components with small circular cross-sections is not easy. Another example of the inadequacy of the standard procedure can be seen in the study of high temperature fatigue in nickel-base superalloys in which corner crack specimens are widely used to replicate typical turbine component geometries and stress levels [9,10]. Therefore, alternative methods to determine both constants of the Paris law are required and desirable. 0013-7944/$ - see front matter � 2012 Elsevier Ltd. All rights reserved. ⇑ Corresponding author. Tel.: +351 239 790 200; fax: +351 239 790 201. E-mail address: [email protected] (R. Branco). Engineering Fracture Mechanics 96 (2012) 96–106 Contents lists available at SciVerse ScienceDirect Engineering Fracture Mechanics http://dx.doi.org/10.1016/j.engfracmech.2012.07.009 engineering. In this kind of components, fatigue, as a consequence of crack initiation and propagation, is recognised as one of the most common failure modes. As is well-known, fatigue crack propagation life can be estimated using the Paris law, da dN ¼ CðDKÞm ð1Þ where C and m are material constants; and DK is the stress intensity factor range. These two material constants are usually calculated using a well-established procedure, based on standard specimens, notched and pre-cracked (BS ISO 12108:2002 and ASTM 647-08e1). The specimens commonly employed (the middle-crack tension and the compact tension geometries) have mode-I through cracks with nearly straight stable shapes without significant shape changes during the propagation. Nevertheless, when considerable shape changes occur [1–3], the above mentioned procedure is not totally adequate. For example, in round bars subjected to mode-I loading, small cracks are practically part-elliptical or part-circular while longer Keywords: Paris law constants Round bar Crack shape evolution Reverse engineering 1. Introduction Mechanical components with cir ison between numerical and experimental results in terms of crack shapes and fatigue lives is carried out in order to find the Paris law constants that best fit the numerical predictions to the experimental results. This technique has been successfully applied to circular cross- section specimens made of S45 steel. Differences between experimental and predicted C and m constants less than 5% and 3%, respectively, were found. � 2012 Elsevier Ltd. All rights reserved. cross-sections, such as shafts, bolts, screws, wires, etc. play an important role in a r t i c l e i n f o Article history: Received 16 May 2011 Received in revised form 9 December 2011 Accepted 7 July 2012 a b s t r a c t In this article, a 3-step reverse engineering technique capable of determining the Paris law constants from the analysis of final fracture surfaces of fatigued round bars is proposed. Ini- tially, at least two crack shapes and the number of cycles between them are obtained experimentally. Then, a 3D-FE automatic fatigue crack growth technique, able to predict C, University of Coimbra, Department of Mechanical Engineering, Coimbra, Portugal ment of Engineering Mechanics, Shanghai Jiaotong University, Shanghai 200240, China Determination of the Paris law constants in round bars from beach marks on fracture surfaces R. Branco a,⇑, F.V. Antunes b, J.D. Costa b, Feng Peng Yang c, Zhen Bang Kuang c a journal homepage: www.elsevier .com/locate /engfracmech Nomenclature a/b crack aspect ratio C, m material constants D diameter of round bar E Young’s modulus E’ modified Young modulus FEM finite element method K stress intensity factor L1 radial size of crack front elements N number of fatigue cycles r radial distance from the crack tip r2 square correlation factor R stress ratio vp crack opening displacement WE work of external forces Dai (j) crack growth increment of ith node for jth iteration Da(j)max maximum crack growth increment for jth iteration DK range of stress intensity factor DKi (j) range of SIF value of ith node for jth iteration DK(j)max maximum range of SIF value R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 97 In this article, a 3-step reverse engineering technique, able to determine the Paris law constants in fatigued round bars from the observation of fatigue crack shapes on fracture surfaces, is presented. Firstly, an experimental test is performed in order to obtain, at least, two different crack shapes on fracture surface of the round bar as well as the number of cycles between them. Secondly, the crack shape evolution is predicted by employing a 3D-FE automatic crack growth technique [11–14]. Finally, a comparison between the numerical and the experimental crack shapes is carried out in order to determine the Paris law constants that fit best the numerical predictions to the experimental results i.e. different Paris law constants are used in the numerical model and the predictions of final crack shape and number of load cycles are compared with the experimental results. m Poisson’s ratio dx, dy, dz displacement components in x, y, z directions 2. Experimental work The experimental fatigue crack growth tests were conducted using 12 mm-diameter and 190 mm-long round bar spec- imens (Fig. 1a). The initial crack was a straight edge crack with 1 mm-depth which was created by using a linear cutting (a) (b) Fig. 1. (a) Specimen geometry; (b) initial crack shape (dimensions in millimetres). machine (Fig. 1b). Both ends of the specimen, with a diameter of 15 mm and length of 50 mm, were fixed to the grips of the machine using button-head connections. In order to minimise bending strains, the major axis of specimen was aligned with the load axis. The material used in this study was the carbon steel S45. Its mechanical properties are given in Table 1 [1]. A MTS809 servo-hydraulic test machine was used to apply the cyclic loading to the specimen. The tests were performed at room temperature, in load control, under a maximum cyclic tension of 25 kN, with a 15-Hz sinusoidal-waveform and stress ratio R = 0.1. The shape and depth growth of fatigue cracks were monitored using a zoom stereomicroscope and beach marking. The beach marks were produced by changing the stress ratio, i.e. the applied load was reduced to one-half for sev- eral cycles. Fig. 2a exhibits the typical beach marks obtained in the fatigue crack growth tests. From the experimental data, several visible crack shapes, at different places of cross-section, were selected. This selection was based on a careful analysis of the crack shape in terms of symmetry. The experimental crack shapes used in the determination of C and m constants (1, 2, A and B) are presented in Fig. 2b. Table 2 shows the polar coordinates (r,h) of each one of them measured according to the Table 1 Mechanical properties of S45 steel at room temperature [1]. Ultimate tensile strength, rm 775 � 106 Pa Monotonic tensile yield strength, r0 635 � 106 Pa Young’s modulus, E 206 � 109 Pa Poisson’s ratio, v 0.3 Fracture toughness, KIC 104 � 106 Pa m0.5 (b)(a) Fig. 2. (a) Fracture surface (initial crack depth of 1.0 mm and 25 kN cyclic tension loading); (b) sketch of several visible crack shapes used. 98 R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 Table 2 Polar coordinates of the experimental crack shapes used. Visible crack shapes 1 (a/D = 0.253) 2 (a/D = 0.418) A (a/D = 0.299) B (a/D = 0.383) Point, i h (�) r (mm) h (�) r (mm) h (�) r (mm) h (�) r (mm) 1 0 3.036 0 5.020 0.000 3.597 0.000 4.598 2 5.036 3.041 4.520 5.019 4.878 3.609 4.465 4.593 3 10.081 3.065 8.642 5.054 9.666 3.635 8.857 4.615 4 15.092 3.081 12.855 5.089 14.451 3.661 13.251 4.654 5 20.065 3.121 16.951 5.133 19.235 3.700 17.644 4.694 6 25.184 3.179 21.128 5.191 24.024 3.752 22.037 4.757 7 30.103 3.235 25.287 5.279 28.809 3.813 26.429 4.819 8 35.125 3.285 29.415 5.370 33.596 3.892 30.822 4.902 9 40.068 3.364 33.558 5.469 38.382 3.978 35.215 5.004 10 45.082 3.468 37.752 5.584 43.168 4.074 39.609 5.102 11 50.019 3.569 41.873 5.717 47.954 4.196 44.002 5.213 12 55.028 3.679 46.023 5.862 52.741 4.305 48.396 5.346 13 60.008 3.809 50.116 6.024 57.525 4.448 52.788 5.482 14 62.791 3.883 52.099 6.103 59.918 4.513 54.985 5.541 15 65.000 3.949 54.234 6.199 62.312 4.581 57.182 5.602 16 67.643 4.038 56.124 6.276 64.514 4.650 59.343 5.680 17 70.060 4.109 58.068 6.381 66.717 4.718 61.505 5.752 referential schematised in Fig. 2a. Only half of the crack front was analysed from which seventeen points were obtained. The numbers of load cycles applied experimentally between the various crack shapes analysed are presented in Table 3. Based on the experimental data, the fatigue crack growth rate calculated in the depth direction is given by Yang et al. [1]: da dN ¼ 1:9037� 10�9ðDKÞ3:256 ð2Þ i.e.,m = 3.256 and C = 1.9037 � 10�9 (da/dN [mm/cycle] andDK [MPa m0.5]). The crack length was measured by analysing the beach marks. For further details of the experimental determination of the Paris law constants, please see Yang et al. [1]. steps (Fig. 3 and fa � � where being Table 3 Number of cycles between the experimental crack shapes. Crack shape combinations 1–2 1–B A–2 Number of cycles 46523 40474 22953 R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 99 Fig. 3. calculat growth model is applied in order to define the crack front advances (Fig. 3d) and to estimate the corresponding number of fatigue cycles. The propagation at each crack front node, under remote mode-I loading, is assumed to occur along the direc- tion normal to the crack front. The crack increment, at an arbitrary node along the crack front, derived from the Paris law, can be expressed as DaðjÞi ¼ DKðjÞi DKðjÞmax " #m DaðjÞmax ð4Þ being DaðjÞi the crack growth increment of ith node for jth iteration, Da ðjÞ max the maximum crack growth increment for jth iter- ation, DKi(j) the SIF range of ith node for jth iteration, DK ðjÞ max the maximum SIF range for jth iteration, and m the Paris law exponent. According to Eq. (4), it becomes clear that the m exponent of the Paris law is the unique material parameter re- Vp is the crack opening displacement, r is the radial distance from the crack tip, E0 is the modified Young modulus, E’ = E/(1 � v2) in plane strain state and E’ = E in plane stress state, and v is Poisson’s ratio. Fourthly, an adequate crack processor GeoStar , included in the commercial FEM package Cosmos/M , is used to solve the FE model. Thirdly, the stress intensity factors along the crack front are computed (Fig. 3c) using the two-point extrapolation method [16]. The stress intensity factors are calculated at two points and extrapolated to the crack tip (r = 0). At a generic point P, located on the crack surface, the mode-I stress intensity factor is defined according to the following expression K ¼ ffiffiffiffiffi p 8r r � E0 � vP ð3Þ cyclically repeated, as schematised in Fig. 3. Firstly, a numerical model representative of the cracked body is created a). This step encompasses the definition of the geometry, boundary conditions, loading, crack shape, elastic constants tigue crack growth rate. Secondly, the displacement field of the crack front nodes is calculated (Fig. 3b). The pre- 3. Automatic fatigue crack growth technique The numerical simulation was performed using Lynx, a specific software capable of simulating fatigue crack growth of planar cracks under mode-I constant amplitude loading [15]. The automatic crack growth technique comprises five main (a) (b) (c) (d) (e) Automatic crack growth technique: (a) 3D FE model definition; (b) calculation of displacement field; (c) calculation of stress intensity factors; (d) ion of nodal advances; and (e) new crack front definition. quired to calculate the crack increment. The number of loading cycles can be calculated by employing the following expression. Nðjþ1Þ ¼ NðjÞ þ DNðjÞ () Nðjþ1Þ ¼ NðjÞ þ Da ðjÞ max C½DKðjÞmax�m ð5Þ Finally, the positions of corner and mid-side nodes of the crack front are moved to their final positions (Fig. 3e) by employing a third order cubic spline function that passes through the provisional corner nodes. This approach gives more realistic crack front shapes and therefore more accurate numerical results. The final crack front is then used directly as the input data of the next iteration and the procedure is repeated as long as no critical values of fracture toughness or crack length are reached. The physical model considered in this research is shown in Fig. 4. Due to material, loading and geometry symmetries, only a quarter of the specimen was simulated. The geometry of ends was simplified taking into consideration their remote posi- tion relatively to the crack front. Movements along x and z (dx = dz = 0) were restrained to simulate the constraint imposed by the high rigidity of the loading machine grips. The crack was assumed to be plane, normal to the axis of the specimen, and existing in its middle-section. Therefore, mode-I loading is expected along the whole crack front. The material was assumed to be homogeneous, isotropic and with linear elastic behaviour. Fig. 5 presents the finite element model developed. The crack front was divided into seventeen corner nodes (Fig. 5d) and sixteen mid-side nodes. A refined region nearby the free surface of round bar was considered. A spider web mesh consisting of three concentric rings with five elements arranged surrounding the crack front was created (Fig. 5e). The collapsed 20- node isoparametric element was used in the first concentric ring. The mid-side nodes were moved to quarter-point positions (Fig. 5c). The standard 20-node isoparametric element (Fig. 5a) was employed in the other two concentric rings of the spider web mesh as well as in the transition and regular meshes. The assembled model (Fig. 5g) had 71743 nodes and 7232 ele- ments. The stress intensity factors at corner nodes were computed using the points A and B identified in Fig. 5e. A plane strain condition was assumed for all positions except at the free boundary where plane stress state prevailed. The numerical procedure was carefully optimised. Several modelling decisions were made taking into account literature results and authors experience [11–21] such as: the spider web mesh topology (Fig. 5e) centred on the crack front; the use of singular elements (Fig. 5c) to simulate the stress singularity; the use of a cubic spline to define the crack front; and the use of a well-known K calculation method. Nevertheless, other important variables were optimised through a specific parametric study. The radial size of crack front elements (L1) plays an important role in the simulation of crack tip stress singularity and 100 R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 Fig. 4. Physical model. R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 101 (a) (d) its opt find th with t ofWE. ent te result proble do no were Lin an sion a calcul Th b) and with t lated d Fig. 5. nodes p imum value depends on the average size of the singular region at the crack tip [19–21]. The criterion adopted here to e optimum size of L1 was the maximisation of the work of external forces. In fact, the discretisation error associated he FEM produces over-rigidity, i.e. decreasesWE. Therefore, the minimisation of error corresponds to the maximisation The optimum value of L1 was found to be D/240. Accurate predictions of fatigue crack growth carried out by the pres- chnique depend greatly on the maximum crack growth increment (Damax) used in Eqs. (4) and (5). The inaccuracies from considering that C(DK)m is constant within each iteration which is false since DK increases continuously. This m can be solved by selecting a relatively small crack increment, small enough to guaranty that numerical results t suffer significant changes. Different maximum crack growth increments (D/50, D/100, D/150, D/200, D/250, D/300) tested. An excellent agreement among the predictions obtained for D/200, D/250, D/300 was observed. For example, d Smith [5,6] used a ratio D/250 in their studies about fatigue crack growth of surface cracks in round bars under ten- nd bending. Smaller ratios were considered here to guarantee the quality of the numerical predictions i.e. D/500 in m ations and D/2000 in C calculations. e crack shape and the crack length were analysed by using the well-known dependent parameters crack aspect ratio (a/ dimensionless crack length (a/D). The variables a and b represent the semi-axes of an ellipse whose centre is coincident he origin of the coordinate system schematised in Fig. 6. The prediction of the m constant was based on the accumu- ifference (ad) defined by the following expression ad ¼ Xn i¼1 jdij ð6Þ (b) (c) (e) (f) (g) Finite element mesh: (a) 20-node hexahedral isoparametric element; (b) 20-node collapsed-element; (c) 20-node collapsed-element with mid-side ositioned at quarter positions; (d) crack front definition; (e) spider web mesh; (f) spider web and transition meshes; and (g) assembled model. being di the difference between the numerical and the experimental coordinates at the ith node of crack front and n the num- ber of nodes. This global parameter analyses the entire crack front and is equal to zero only when both experimental and Fig. 6. Definition of the dependent parameters used. 102 R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 numerical crack shapes are superimposed. Fig. 7 presents typical fatigue crack shape developments obtained after the optimisation of the numerical procedure (m = 3). Three different surface crack configurations subjected to constant amplitude cyclic tension are examined, namely a part-circular crack (a0/b0=1), part-elliptical crack (a0/b0=0.4) and straight crack (a0/b0=0). A strong dependence on the initial crack shape is observed. It is clear that for the initial straight shape (a0/b0 = 0), the crack grows much more rapidly in depth direction than along the free surface whilst for the part-circular shape (a0/b0 = 1) the growing is more balanced along the whole crack front. However, this relevant effect of the initial crack configuration existing in the early stage is gradually weak- ened as the crack extends and consequently the crack fronts become similar. The amount of crack growth needed to achieve this part of propagation also depends on the initial crack configuration and therefore shapes closer to the preferred propa- gation path reaches it faster than the others. The numerical predictions were compared with the literature results for validation. Fig. 8a plots a/b against a/D for dif- ferent initial crack configurations (a0/b0=0, a0/b0=0.2, a0/b0=0.4, a0/b0=0.6, a0/b0=0.8, a0/b0=1) with the same length (a0/D=0.1). As can be seen, the trajectory drawn by the crack strongly depends on the initial crack aspect ratio. Nevertheless, this high dependence on initial crack shape gradually weakens, as the crack propagates, leading the crack shape to preferred propa- gation paths. Additionally, it is possible to observe a good agreement between the numerical results presented here and those found in the scientific literature [2–5,7] for the same propagation conditions (a0/D = 0.1, m = 3, m = 0.3). The effect of the Paris law exponent on crack aspect ratio, already discussed by several authors [3–7], is also exhibited in Fig. 8a. In the early propagation period, no relevant changes are distinguished. However, as the crack grows, different propagation paths are clearly observed. On the other hand, after a certain time, the propagation paths cross-over (a/D � 0.48) following from then on divergent trajectories. The cross-over point depends on several variables such as the value of m, loading type and initial crack shape [3–5,7,8]. Such a study can be used to select the most adequate experimental cracks shapes. Ideally, experimental defects in a region in which the crack shape is more sensitive tom are desirable. In this regard, crack shapes in the early period of propagation are not recommended because the shape modifications are not yet obvious. For example, in (b) (a) (c) Fig. 7. Fatigue crack shape development for a surface crack subjected to tension: (a) a0/b0 = 1, a0/D = 0.1; (b) a0/b0 = 0.4, a0/D = 0.1; (c) a0/b0 = 0, a0/D = 0.1. R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 103 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 a/ b a/D Carpinteri et al. [4] Lin et al. [5] Shin et al. [2, 3] Toribio et al. [7] a0/D=0.1 m=3 ν=0.3 m=4 m=2 0.9 1.0 a0/b0=1 (a) the case presented in Fig. 8a (a0/b0 = 0 and a0/D = 0.1), the three curves tend to follow an unequivocal path for values of a/D greater than 0.2. In that sense, all experimental crack shapes used in this study were selected, as indicated in Table 2, within that range. Fig. 8b shows the variation of Kmin/Kmax along the crack front during the crack growth for a surface crack of a round bar subjected to constant amplitude cyclic tension. Different initial crack configurations (a0/b0 = 0, a0/b0 = 1) with the same dimensionless crack length (a0/D = 0.1) were considered. This ratio (Kmin/Kmax) is interesting to characterise the SIF variations along the crack front with the crack growth. It increases suddenly at the early stage and then goes down slightly to values about 0.85–0.88. After that, the Kmin/Kmax ratio rises up slightly towards values close to unity (Kmin/Kmax � 0.98). Additionally, at the early stage, the gradient of ratio Kmin/Kmax is less intense for the initial part-circular crack shape (a0/b0 = 1) than for the initial straight crack shape (a0/b0 = 0). This explains the more significant shape changes observed in the latter case (Fig. 7c). Moreover, the results predicted in this study are in excellent agreement with those of Lin and Smith [5] for the same prop- agation conditions (a0/D = 0.1,m = 3, m = 0.3). Unlike other cases, the propagation of surface cracks in round bars subjected to mode-I loading does not tend to iso-K profiles. 4. Determination of the Paris law constants The first step encompassed the determination of them constant. A set of numerical simulations was performed using val- ues of m contained in the interval 2.6–3.6 (da/dN [mm/cycle], DK [MPa m0.5]). The crack shape 1 (a/D = 0.253), exhibited in Fig. 2b, was defined as the initial crack configuration. Each numerical simulation was interrupted when the first node of crack front (node 1 of Fig. 5d) reached the length of crack shape 2 (a/D = 0.418) of Fig. 2b. Then, for each case, the values of the accumulated difference (ad) between the numerical and experimental crack shapes were computed using Eq. (6). This dependent parameter has proven to be very sensitive to any shape change therefore its use is recommended for this purpose. 0.6 0.7 0.8 0.0 0.2 0.4 0.6 0.8 K m in /K m ax a/D Lin et al. [5] a0/b0=0 a0/D=0.1 m=3 ν=0.3 (b) Fig. 8. Evolution of: (a) a/b with a/D; (b) Kmin/Kmax with a/D. 104 R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 0.00 0.35 0.70 1.05 1.40 2.6 2.8 3.0 3.2 3.4 3.6 ad [m m] Exponent of Paris law, m ad(m) = 1.354m2 - 8.375m + 13.372 3.102 da/dN [mm/cycle], ΔK [MPa m0.5] 70 -1.006 (a) In corner crack specimens, a different global parameter able to characterise asymmetrical crack shapes has been used [22]. Fig. 9a plots the accumulated difference (ad) against the Paris law exponent. A well-defined tendency emerges from the re- sults. This is clear evidence of the suitability of the selected dependent parameter. The resulting data were fitted by the least square method to a second order polynomial function with a relatively high square correlation (r2 = 0.995). In theory, the correct value of m can be found by minimising the value of ad. So, it means that the derivative of ad (ad0) must be equal to zero: ad0ðmÞ ¼ 0() m ¼ 3:102 (da/dN [mm/cycle], DK [MPa m0.5]). The error involved in the prediction of m, in relation to the experimental values, is of 4.73% which is clearly acceptable in this context. Note that due to the relatively low stress ratio (R = 0.1) used, crack closure can exist which can affect the crack shape and therefore them prediction. In order to avoid such a phenomenon, higher stress ratios can be adopted. The second step comprised the determination of the C constant. The value of m previously predicted was fixed. Then, a new set of numerical simulations was performed using values of C within the range of 7 � 10�9 to 2.4 � 10�9 (da/dN [mm/ cycle], DK [MPa m0.5]). These simulations were computed using a maximum crack growth increment equal to D/2000. In a similar manner, the crack shape 1 (a/D = 0.253) of Fig. 2b was considered as the initial crack configuration. Each numerical simulation was interrupted when the first node of crack front (node 1 of Fig. 5d) reached the length of crack shape 2 (a/ D = 0.418) of Fig. 2b. The number of loading cycles between these two prescribed crack shapes, computed using Eq. (5), was the result of the numerical simulation. The determination of C could also be done by integrating the Paris law (Eq. (1)). However, this approach requires a SIF-based solution which, in general, is computed with a certain level of error. Due to this fact and in order to take advantage from the numerical procedure developed here, the implementation of a SIF-based solution was not followed. Fig. 9b compares the number of loading cycles predicted with the number of cycles obtained experimentally. An exponential function was fitted to the data by employing the least square method (r2 = 0.990). The C constant was found by equalising the fitted function to the experimental number of cycles, i.e. solving 30 40 50 60 1.5 1.7 1.9 2.1 2.3 2.5 N um be r o f c yc le s, N [ 10 3 ] Constant of Paris law, C [×10-9] N=46523 cycles (exp.) 1.8489 da/dN [mm/cycle], ΔK [MPa m0.5] N(C) = 85.895C (b) Fig. 9. Determination of: (a) m constant: (b) C constant. the following equation: NðCÞ ¼ 46523() C ¼ 1:8489� 10�9 (da/dN [mm/cycle], DK [MPa m0.5]). Numerical and experi- mental results are in good agreement. Note that the difference between both values is less than 2.87%. Naturally, this error magnitude can be considered acceptable. In order to evaluate the robustness of the proposed technique, the procedure was carefully repeated for other combina- tions of experimental crack shapes exhibited in Fig. 2b. Two different situations, corresponding to the combinations 1 to B (1–B) and A to 2 (A–2), were studied. For each of them, as described previously, the accumulated difference was computed through Eq. (6) for various values ofm. From the values of ad, a second order polynomial function was fitted. Then, the expo- nent of the Paris law was achieved by minimising this function. Next, with the predicted value of m, new numerical simu- lations were performed for different C constants aiming at obtaining the numerical number of cycles between both crack shapes. These data were fitted to an exponential function again. Equalising this function to the corresponding experimental number of cycles presented in Table 3, the C constant was calculated. Finally, the predicted constants were compared with the experimental values. Table 4 presents the achieved C andm constants as well as the errors relatively to the experimental results. Regardless of the combination analysed, the error magnitude is similar. For example, in the combination 1–B, no sig- nificant effects were introduced due to the use of a different first visible crack shape; in the combination B–2, the influence of a different second visible crack shape was not relevant. These results reinforce the robustness of the proposed approach. The means of the C and m values (C = 1.8629 � 10�9 and m = 3.101) were considered the final predictions of the Paris law con- stants. Notwithstanding, taking into account the consistency of the predictions presented, this assumption is not expected to have a strong influence on the following analysis. The final numerical values of C and m were used to carry out an entire simulation from the initial straight crack shape used in the experimental tests. Then, the numerical and experimental crack shapes were compared. This comparison was done through the difference di, schematised in Fig. 6, assuming that both crack shapes were superimposed at the symmetry line of the cross-section of specimen. For each crack front, the evolution of di/ri with hi was computed, being ri the experi- mental radius. The results calculated are presented in Fig. 10 for crack fronts X, 1, A, B, 2 and Y of Fig. 2b. As can be seen, these values have well-defined limits that vary between �1% and 5%. Besides, although some exceptions are observed, it is possible to distinguish a dominant tendency for the curves in which the differences increase progressively towards the surface. These results demonstrate that the reverse engineering technique proposed here is able to obtain the Paris law con- stants from materials in the form of round bars. Table 4 Predicted C and m constants (da/dN [mm/cycle], DK [MPa m0.5]) for different combinations of the experimental crack shapes used. Crack shape combination C m Error of C relatively to the experimental value (%) Cexp�CnumC Error of m relatively to the experimental value (%) mexp�mnummexp R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 105 exp 1–2 1.8489 � 10�9 3.102 2.87 4.73 1–B 1.8983 � 10�9 3.101 2.19 4.76 A–2 1.8615 � 10�9 3.099 1.35 4.82 -9 -6 -3 0 3 6 9 0 15 30 45 60 75 d i /r i [% ] θi [º] X 1 A B 2 Y 5% -1% Fig. 10. Evolution of di/ri with hi for several crack fronts. 106 R. Branco et al. / Engineering Fracture Mechanics 96 (2012) 96–106 5. Conclusions In the present paper, a reverse engineering technique capable of determining both Paris law constants from the analysis of fracture surfaces of fatigued round bars was proposed. The technique consists of three main tasks. Firstly, experimental work aiming at obtaining at least two visible crack shapes and the number of loading cycles between them is performed. Secondly, a three-dimensional fatigue crack growth technique able to predict the crack shape evolution and the fatigue life is employed. Thirdly, a set of numerical simulations is carried out by using different values of m. The m constant is found by minimising the shape differences between the experimental and numerical crack fronts. Then, using the calculated m value, new numerical simulations are done for different C values. The correct C value is found by equalising both the numerical and the experimental numbers of loading cycles. This technique was successfully applied to a 12 mm-diameter and 190 mm-long circular cross-section specimen made of S45 steel subjected to tension loading. Experimental work was developed at constant amplitude ratio (R = 0.1) to obtain visible crack shapes and the number of cycles between them. Marking was produced by changing the stress ratio. A three-dimensional finite element model was created to predict the crack shape evolution and the fatigue life by using the FCG software Lynx. The numerical procedure was optimised and the results of crack shape evolution and stress intensity fac- tors were successfully compared with values published by other authors. Additionally, as reported in the literature, it was demonstrated that the initial crack configuration affects significantly the early propagation period while the subsequent propagation is influenced by the Paris law exponent. Three predictions of C andm from three different combinations of experimental crack shapes were carried out. The values found were clearly convergent. The maximum errors in C and m, relatively to the experimental constants, were 4.82% and 2.87%, respectively. Although not perfect, both errors are entirely acceptable in this context. The proposed values of C and m were calculated by means of three predictions obtained from different crack shapes, being C = 1.8629 � 10�9 and m = 3.101, respectively. The differences between the numerical and the experimental results evidence accurate predictions in C and m and therefore the reverse engineering technique proposed here is very promising for the determination of the Paris law constants of materials in the form of round bars. Acknowledgements The authors are indebted to the Portuguese Foundation for the Science and Technology (FCT) through COMPETE program from QREN and to FEDER (European Regional Development Fund) for the financial support (Project PTDC/EME-PME/114892/ 2009). References [1] Yang Feng Peng, Kuang Zhen Bang, Shlyannikov VN. Fatigue crack growth for straight-fronted edge crack in a round bar. Int J Fatigue 2006;28:431–7. [2] Shin CS, Cai CQ. Experimental and finite element analyses on stress intensity factors of elliptical surface crack in a circular shaft under tension and bending. Int J Fract 2004;129:239–64. [3] Shin CS, Cai CQ. Evaluating fatigue crack propagation properties using a cylindrical rod specimen. Int J Fatigue 2007;29:397–405. [4] Carpinteri A. Shape change of surface cracks in round bars under cyclic axial loading. Int J Fatigue 1993;15:21–6. [5] Lin XB, Smith RA. Shape growth simulation of surface cracks in tension fatigued round bars. Int J Fatigue 1997;19:461–9. [6] Lin XB, Smith RA. Fatigue growth simulation for cracks in notched and unnotched round bars. Int J Mech Sci 1998;40:405–19. [7] Toribio J, Matos JC, González B, Escuadra J. Numerical modelling of crack shape evolution for surface flaws in round bars under tensile loading. Engng Fail Anal 2009;16:618–30. [8] Couroneau N, Royer J. Simplified model for the fatigue crack growth analysis of surface cracks in round bars under mode I. Int J Fatigue 1998;20:711–8. [9] Tong J, Byrne J, Hall R, Aliabadi MHA. Comparison of corner notched and compact tension specimens for high temperature fatigue testing. In: Proceedings of the conference engineering against fatigue, 17–21 March, University of Sheffield, UK; 1997. [10] Brown CW, Hicks HA. Fatigue growth of surface cracks in nickel-base superalloys. Int J Fatigue 1984;4:73–81. [11] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates. Part I: the numerical technique. Engng Fract Mech 1999;63:503–22. [12] Branco R, Antunes FV. Finite element modelling and analysis of crack shape evolution in mode-I fatigue middle-cracked tension specimens. Engng Fract Mech 2008;75:3020–337. [13] Lee WY, Lee JJ. Successive 3D analysis technique for characterization of fatigue crack growth behaviour in composite-repaired aluminum plate. Compos Struct 2004;66:513–20. [14] Gardin C, Courtin S, Bézine G, Bertheau D, Ben Hadj Hamouda H. Numerical simulation of fatigue crack propagation in compressive residual stress fields of notched round bars. Fatigue Fract Engng Mater Struct 2007;30:231–42. [15] Branco R, Antunes FV, Costa JD. A user-friendly computer application for simulating fatigue crack growth of planar cracks by using the FEM. Comput Appl Engng Educ. doi: http://dx.doi.org/10.1002/cae.20578, in press. [16] Guinea GV, Planas J, Elices M. KI evaluation by the displacement extrapolation technique. Engng Fract Mech 2000;66:243–55. [17] Antunes FV. Influence of frequency, stress ratio and stress state on fatigue crack growth in nickel base superalloys at elevated temperature. PhD thesis, Department of Mechanical and Manufacturing Engineering, University of Portsmouth, United Kingdom; 1999. [18] Nykänen TJ. Fatigue crack growth simulations based on free front shape development. Fatigue Fract Engng Mater Struct 1996;19:99–109. [19] Murti V, Valliappan S. A universal optimum quarter point element. Engng Fract Mech 1986;25:237–58. [20] Branco R, Antunes FV, Martins RF. Modelling fatigue crack propagation in CT specimen. Fatigue Fract Engng Mater Struct 2008;31:452–65. [21] Branco R, Antunes FV, Ricardo LCH, Costa JD. Extent of surface regions near corner points of notched cracked bodies subjected to mode-I loading. Finite Elem Anal Des. 2012;50:147–60. [22] Branco R, Antunes FV, Martins Ferreira JA, Silva JM. Determination of Paris law constants with a reverse engineering technique. Engng Fail Anal 2009;16:631–8. Determination of the Paris law constants in round bars from beach marks on fracture surfaces 1 Introduction 2 Experimental work 3 Automatic fatigue crack growth technique 4 Determination of the Paris law constants 5 Conclusions Acknowledgements References