Numerical determination of Paris law constants for carbon steel using a two-scale model

For most engineering alloys, the long fatigue crack growth under a certain stress level can be described by the Paris law. The law provides a correlation between the fatigue crack growth rate (FCGR or da/dN), the range of stress intensity factor (ΔK), and the material constants C and m. A well-established test procedure is typically used to determine the Paris law constants C and m, considering standard specimens, notched and pre-cracked. Definition of all the details necessary to obtain feasible and comparable Paris law constants are covered by standards. However, these cost-expensive tests can be replaced by appropriate numerical calculations. In this respect, this paper deals with the numerical determination of Paris law constants for carbon steel using a two-scale model. A micro-model containing the microstructure of a material is generated using the Finite Element Method (FEM) to calculate the fatigue crack growth rate at a crack tip. The model is based on the Tanaka-Mura equation. On the other side, a macro-model serves for the calculation of the stress intensity factor. The analysis yields a relationship between the crack growth rates and the stress intensity factors for defined crack lengths which is then used to determine the Paris law constants.


Introduction
Structural health monitoring and detection of failures are highly important for fitness and service assessment of structures. Failure of structures under fatigue loading can occur at load levels below the yield stress of the used material. Therefore, it is of special importance to be able to make predictions of life cycles until catastrophic fracture occurs. A well-known example of such catastrophic failure is the huge train accident in Eschede where one wheel of the train broke due to cyclic loading (figure 1) that has not been considered during construction [1]. On the other hand not every single crack has to be critical immediately. Often structures can withstand cracks up to a certain crack length until unstable fracture occurs. In order to determine the point where crack growth becomes unstable, it is necessary to develop methods to describe crack growth mathematically. A very common and often used method for the characterization of the long crack growth is the Paris law which gives a relationship between the fatigue crack growth rate (FCGR or da/dN) and the stress intensity factor (ΔK) at the crack tip during the stable crack growth [2], equation (1). A typical fatigue growth rate curvealso known as da/dN versus ΔK curveis shown in figure 2. The curve is defined by regions I, II and III. The Paris law relationship can be visualized as a straight line for the region of stable crack growthregion II [3]. The accompanying mathematical equation contains two material parameters C and m, where m represents the slope of the line in figure 2 and C the y-axis-intercept [4]: The crack growth rates in the region II are typically in the order of 10 -9 to 10 -6 m/cycle and correspond to stable crack growth. The constants C and m are usually determined in experiments [5][6][7][8][9] and depend on the material and various influencing factors such as temperature, environmental medium and loading ratio [5,10]. The last is probably the most significant and usually results in closely spaced lines parallel to each other. For metals the exponent is typically of the order m = 2…4 [10].

Figure 2.
Geometry of the used specimen with a central crack [3].
Since the experimental determination of the Paris law constants is typically tedious and time consuming, the objective of this paper is to determine them numerically considering the influence of microstructure on the crack growth rate.

Modelling approach and details
The simulations for determination of the Paris law constants are done with the Finite Element (FE) software ABAQUS and by applying a two-scale model. The model contains the microstructural submodel of a material that is used to calculate the fatigue crack growth rate da/dN at a crack tip and the macro-model that serves for the calculation of the stress intensity factor ΔK. In order to get the typical da/dN versus ΔK curve, the crack growth rate and stress intensity factors for six different crack lengths are determined. The constants are then determined from the interpolated straight line, representing region II in figure 2, as the slope m and the y-axis-intercept C. The input data for the common diagram is calculated using the two aforementioned approaches, as shown in figure 3. On the one hand the stress intensity factor is determined using the full-scale macromodel with a center crack and on the basis of Linear Elastic Facture Mechanics (LEFM). The crack growth rate on the other hand is determined using the submodel based on the Tanaka-Mura equation [11,12]; the submodel is placed in the area behind the crack tip. The grain structure is assigned to the submodel to consider micro-crack initiation processes in the vicinity of the crack tip.
The Tanaka-Mura equation is typically used to estimate the number of cycles that are attributed to the crack nucleation in a single grain [11,12]. Jezernik et al. [13,14] introduced a modification of the original equation in the sense that the crack does not form instantaneously through the whole grain but it forms in segmental manner: where is the number of cycles required to form a crack segment in a single grain, i.e. crystal.
Furthermore, equation (2) contains two microstructure-related parameters, namely the length of slip line segment and the average shear stress range ∆ ̅ on it. Other material constants (such as shear modulus , crack initiation energy and Poisson's ratio ) can be either found in literature for most materials or can alternatively be obtained experimentally. The critical resolved shear stress (CRSS) is particularly interesting since it can be obtained by means of Molecular Dynamics (MD) simulations [15]. The CRSS represents a critical value of the shear stress along the slip direction that must be overcome for the dislocation to move, i.e. if the magnitude of the resolved shear stress is below the value of the CRSS, no dislocation movement is allowed and consequently no pile up at the grain boundary takes place [16].

Geometry of the macro-models
In order to accomplish the goal of the present paper it was necessary to create two slightly different full-scale macro-models. Both models represent the same geometry, however, the ways of modelling the crack for determination of stress intensity factor ΔK on the one side and fatigue crack growth rate da/dN on the other require certain adjustments at the regions of interest, i.e. on the crack-affected path. More specifically, the macro-model for the calculation of ΔK requires the usage of some special ABAQUS techniques to represent the crackin this case the seam crack is usedwhile the macromodel (further: global model) for determination of da/dN needs to be geometrically adjusted to the submodel that is placed at the tip of the structural crack. A seam defines an edge or a face in a model that is originally closed but can open as a crack during an analysis [17].
The models were built on the basis of a specimen with central pre-crack prepared for fatigue testing, shown in figure 4, considered in the paper of Božić et al. [3]. At the initial state the specimen has a notch and a pre-crack of 2a = 8 mm (figure 4, A -Crack detail). Only half of the specimen has to be modeled due to symmetry with respect to the vertical y-axis. The symmetry with respect to the x-axis is aligned with the crack and, therefore, cannot be used in either case; for ΔK determination the applied seam cannot extend along the boundaries of a part and must be embedded within a face of a two-dimensional (2D) part or within a cell of a solid part [17]. In the other case, the used submodel (microstructural model) at the crack tip area requires the transfer of boundary conditionsin this case displacementsfrom both parts of the global model, the upper and the lower one. The structural crack lengths which were used for this study were taken from [3], where ΔK values were determined using the FE software ANSYS. Nevertheless, ΔK values were calculated againthis time with ABAQUSfirst to compare the results with previous simulations and second to ensure that the models are built properly. The considered crack lengths for the simulations are listed in All models which are used here are built as 2D, in accordance to the geometry of the specimen and the fact that the main central part has a thickness of just 4 mm along the z-direction. The areas in the upper and lower part of the models with relatively higher thickness (54 mm) were also modelled as 2D. Those two parts are used to apply loading conditions and movement constraints, respectively, as shown in figure 5. The stress state in the middle part of the cracked plate is assumed to be plane stress. The thicknesses of the specimen in both regions, the central and outer, were considered by assigning plane stress thickness to their belonging sections.
As already mentioned only one half of the specimen has to be modelled. The symmetry is realized by using boundary conditions on the vertical y-axis which fix displacements in x-direction as well as rotations of any kind (figure 5). The material behaviour is assumed to be purely linear elastic; only a small plastic zone at the crack tip is expected ( figure 7), therefore, no plastic material data are necessary to be used. The isotropic material data for the specimen made of conventional mild carbon steel were adopted from the study of Božić et al. [3]. The material parameters applied to both models are: Young modulus E = 206 GPa, Poisson's ratio = 0.3, yield stress Re = 235 MPa and shear modulus G = 80 GPa. It is opportune to introduce the experimentally obtained Paris law constants for the selected material in order to provide validation data for the numerical results that follow later; m = 2.75 and C = 1.43x10 -11 [3].
As the material data, loading and boundary conditions were taken from [3]. The testing specimen was exposed to constant amplitude cyclic tension load in a hydraulic fatigue testing machine. The load was applied to the pin which was placed in the hole in the upper part of the model while the pin in the lower hole was fixed ( figure 5). The force range and the stress ratio applied to the half-model are denoted by F = Fmax -Fmin and R = Fmin ⁄ Fmax, and were F = 76 800 N and R = 0.0253 [3]. In contrast to the experiments, simulations on the macro-models were performed with static loading conditions, however, in combination with LEFM for the ΔKs determination and with the Tanaka-Mura equation for the determination of da/dN. For both models continuum plane stress 4-node bilinear elements with reduced integration and hourglass control (CPS4R) were used. In order to deal with the stress singularity at the crack tip in the case of ΔKs calculation, special elements have to be used that are able to indicate the infinite stresses properly. In ABAQUS this is done by collapsing one side of an 8-node isoparametric element connected to the crack tip [17], as can be seen in figure 7.

Submodel (microstructural model) details
The submodel or microstructural model is placed right at the tip of the global model structural cracks (table 1) where their extension is expected. As already mentioned, the Tanaka-Mura equation is typically used to estimate the duration of the short crack initiation stage. In this paper, however, the equation is applied within the microstructural model ( figure 6) with the aim to estimate the rate of the infinitesimal crack extension. The geometry of the submodel can be seen in figure 6, as well as the location where it is placed with respect to the global model. Its size is selected to 0.4 mm x 0.4 mm including the tip of the structural crack of 0.1 mm radius, which can be identified as a notch. As mentioned before, displacements of the global model are applied to the coinciding edges of the submodel, marked bold in figure 6. The submodel was very fine meshed with the same elements as the global model (CPS4R); what in the end gives smooth stress distribution in the undamaged as well as in the damaged submodel (figures [9][10][11][12]. To depict the mesh fineness, a single grain has 120 elements in average.

Results
The stress intensity factor ΔK as well as the J-Integral and other fracture mechanics characteristics, can be requested in ABAQUS as history output data. In order to get ΔK for each individual crack of six in total from table 1, the same number of variations of the macro-model were modelled. Figure 7 gives the von Mises stress distribution at the crack tip of a macro-model. As already mentioned in Chapter 2, the geometry itself with the structural pre-crack stays the same for all variations, however, the used seam length varies. By evaluating ΔK for each crack length and putting the results into a common diagram gives the expected linear relationship (see region II in figure 2) between ΔK and the crack length a (figure 8).  orientation and are distanced one from each other with an offset. Furthermore, each slip band is divided in four segments meaning that a crack does not form instantaneously through the whole grains but it forms in a segmental manner. The condition for the formation of crack segments is contained in the Tanaka-Mura equation (2) and it reads as follows: the average shear stress ∆ ̅ on a slip band segment needs to exceed two times the CRSS. Additionally, the equation (2) gives the number of cycles dN (i.e. Ns) that are spent for every formed crack segment. Furthermore, after the crack condition has been satisfied in one step, the model gets updated with the latest crack and remeshed in the following step where the condition is applied again and a new weakest crack segment is traced. The von Mises stress state before the first formed crack segment is shown in figure 9.     Generally what happens for all considered structural cracks is that the formation of cracks in the submodel ceases after a certain number of stepsthe number varies from one model setup to another. The reason for this is that grains that are favourable for cracking on the basis of the aforementioned condition fade out.
In order to determine the fatigue crack growth rate da/dN, both the cycles dN and the accompanying length da of each individual crack segment that formed in the microstructural model have to be determined. The crack length da can easily be quantified using the ABAQUS graphical interface or can be rather gathered from output data.
By taking the measured da and the correlated dN one can evaluate the growth rate da/dN for each step. It is necessary to indicate that the common growth rate da/dN is not a constant value, but it fluctuates as it can be seen in figure 13 for the case of the 20.1 mm long structural crack. It can be noticed that figure 12 contains 25 broken segments while only 11 are considered in figure 13  da/dN is plotted. An explanation for this is that starting from step 12 in figure 13 the number of required cycles for further formed cracks becomes significantly higher, meaning directly that the da/dN drops down to a negligible value. This leads to the assumption and consequence that the crack growth rate should be neglected in this relatively advanced stage. In order to get a typical da/dN versus ΔK plot as in figure 2, the fluctuating crack growth rate from, e.g. figure 13, needs to be averaged. For the structural crack length of 20.1 mm this results in da/dN = 8.49421x10 -9 , which seems to be appropriate according to literature [3,19]. This averaged growth rate can be considered as extensional growth rate for the structural crack. The same procedure was applied for remaining structural crack lengths. Finally, the results for da/dN and for ΔKs (figure 8) are put into a common log da/dN versus log ΔK diagram, figure 14. The resulting six single points were interpolated with a straight line from which the material constants of Paris equation (1) were determined, the slope of the line is m = 2.75 and the y-axis intercept C = 3.8x10 -12 , which agree quite well with the experimentally determined values of m = 2.75 and C = 1.43x10 -11 [3].

Discussion and conclusions
The aim of this paper was to numerically determine the Paris law material constants C and m and this was successfully accomplished by applying a two-scale model. Namely, stress intensity factors for six different structural cracks were calculated using the macro-model based on classical LEFM. In the second part, the fatigue crack growth rates at the tips of same structural cracks were determined by using the microstructural model enriched by the Tanaka-Mura equation.
Having determined the stress intensity factors and the crack growth rates separately, the values have been visualized afterwards in the common log-log diagram ( figure 14). Considering that the Paris law can be approximated by a straight line in the diagram for the case of stable crack growth, the material constants have been determined then as m = 2.75 and C = 3.8x10 -12 .
Since this study was the first, to the authors' knowledge, attempt to determine the Paris constants in a numerical way as well as taking into account the microstructure of the material, and considering that the literature values of the investigated material (m = 2.75 and C = 1.43x10 -11 ) match the numerically obtained ones quite good, the study and its working hypothesis seem to prove their worthiness. It is to be supposed that the results determined here may even become better during further research.