FM models of creep crack growth (CCG)

The FM approach uses the energy needed for crack development as the main source of information about the process. For a given crack length, it is possible to calculate the critical loading or alternatively the critical crack length for a given loading.

In the first stage of failure development, there is no observable crack, so it is impossible to describe this stage in terms of FM. However, it is possible to correlate this period with FM parameters and determine the time of crack initiation (cf. [21,22]).

In the second stage, FM is mainly used in the case of a single crack. If the stress field around the crack tip is known, the J-integral or its creep equivalents C* or C(t) can be calculated. The definition of the C* contour integral is similar to the J-integral but strains and displacements are replaced by their rates (cf. [22,23]):

When the stationary creep does not occur, the stress field around the crack tip can be described by the C(t) parameter [24]:

This integral is not path independent, it has to be evaluated near the crack tip. For a large t, when creep becomes stationary, C(t) approaches the C*-integral, for any t: C(t) > C*. It was shown that in stationary creep conditions, the C*-integral can be correlated with the crack growth rate (cf. e.g. [21,25,26]).

FM estimation of lifetime

The time for crack initiation t_i is not directly described by FM models. There are some estimations, and the proposition of Tan, Celard et al. [21] is applied in the current work. They defined a lower and upper approximation. The first formula uses the typical relation between the C* parameter and crack growth rate:

The upper approximation t_iU is achieved under assumption that the initial crack growth rate is (n+1) times smaller than the rate in a stationary state, where n is creep index in the Norton creep eq. (8):

The comparison of the results of experiments for different steels with the approximations defined above indicated that the crack initiation times fall within the expected range. The lower approximation was closer to the experimental results in most of the examined cases (cf. [21]).

The next step in the determination of the safety margin is the estimation of the crack growth time period t_f-t_i. This is usually determined though the integration of eq. (3), starting from the initial crack length a_i=da up to the final crack length a_f which is equal to the total width of the element (cf. [27]). To perform this integration, the relationship between the value of the C* parameter and the current crack length has to be established. This relationship can either be obtained analytically (cf. [22]) or numerically. This approach is limited by the assumption that the crack grows in a stationary environment, i.e. the stress distribution in the vicinity of the crack tip is described by stationary creep equations. This is possible only when the rate of stress redistribution is faster than the rate of crack growth. If this assumption is not satisfied, the C*-integral cannot be used and instead, other parameters, specific to the transition creep period, like C(t) or C_t, can be applied (cf. [28]).

The material parameters D, φ (eq. (3)) can be found experimentally or determined by some analytical formulas (see [22]). There is also an approximate solution which describes the creep crack growth quite well for most of the cases (cf. [5]):

Uniaxial tension creep damage model

All of the formulas of FM are valid only in situations where an initial crack exists. Parameters like C*, C(t) are undefined for uncracked specimens. There is no such limitation in the case of damage mechanics models. The damage develops in material that is initially undamaged, which leads to the initiation of a microcrack. Therefore, DM can be used for determination of the time t_i. It is also possible to model the flaws or microcrack systems in original material by setting an appropriate initial damage field.

When the crack starts to develop, the application of DM is more troublesome. However, attempts to describe the crack development by growth of the damage field up to its critical value and further development are often made (cf. e.g. [5,25,29,30,31]).

This approach, known as the local approach to fracture (LAF) encounters a series of problems [32]. Applications of continuous equations in the solving of DM problems (named continuous damage mechanics [33]) are limited by the size of the smallest element in which system variables can be considered continuous, called representative volume element (RVE). The requirements for RVE are contradictory. RVE should be large enough to contain the representative number of defects and sufficiently small that the variations of stress and strain values are relatively small. These requirements can be more or less satisfied for a problem with undamaged or randomly distributed damage to the material (cf. [34]). The damage distribution is usually corelated with the microstructure of the material and the size of RVE is then adopted to the characteristic dimensions of this microstructure.

In the original Kachanov model [2] and in its successors (e.g. [4,18,35]), the development of damage is dependent upon the stress state. There are also strain-based damage models (e.g. [15,36,37]) in which the critical value of the damage parameter is achieved when the ductility of material is exhausted. The current work is mainly based on the Kachanov proposition of creep damage development. Additionally, it is assumed that the creep strains are associated with damage (cf. [3]). The uniaxial constitutive equations are as follows:

Influence of stress triaxiality

The equations of the development of creep damage are fitted to experimental results in the uniaxial tension creep test. While the triaxial stress state is achieved in the vicinity of the crack tip (cf. e.g. [15,18,35,38]). Hayhurst, who analysed the growth of damage in the multiaxial state, assumed that effective stress σ_eff is responsible for damage development and it takes the form [39]:

In models where the damage parameter is based on strain, the ductility is described by the Manjoine equation (cf. e.g. [4,15,40]):

Here, the ductility is a function of the triaxiality coefficient defined as ratio σ_m/σ_e. In the vicinity of the crack tip, the value of this coefficient is large, so the strain at failure is small (cf. [18]). In Fig. 1 it can be seen that the strain at failure calculated from the Kachanov model (eqs. (7–8, 11)) is close to the Manjoine model in this range.

Fig. 1.

Comparison of the ductility obtained on the basis of the Kachanov and Manjoine models

Non-local damage model

The main problem in the application of continuous DM equations is strong localised damage and strain occurring in the numerical solution of crack problems. It results in the spurious mesh dependence of the solution (cf. [41]).

There are many possibilities to prevent excessive damage localisation, such as the mesh-dependent softening modulus, the limitation of mesh size, artificial viscosity, the Cosserat continuum, the gradient method and the non-local integral method (see [42] for review). The non-local damage theory seems to be one of the most popular theories. The idea is that the constitutive equation contains not only local values but also some non-local parameters obtained by averaging the chosen value in a particular volume. The motivation for the use of non-local theory is the observation that the growth of the crack does not depend on the state at a given point but on energy release in its vicinity (cf. [43]). This method is also used in the current work.

The definition of the non-local value of the state variable is:

The value of weighted function depends on the relative position of points x and r. The most popular version of it is Gaussian distribution:

Theoretically, the size of the neighbourhood V(x) can be unlimited, but for practical reason, it is limited to some extension of characteristic length. The smallest distance between points x and r, at which the weighted function vanishes or becomes negligible, is called the non-local interaction radius (cf. [43]). There are also methods for which the interaction radius is exactly equal to the characteristic length (e.g. grid method [13,48]).

Experiment description

The analysis was performed for a compact tension (CT) specimen made of 316 stainless steel creeping at a temperature of 550°C. Its results were compared with experimental data [18]. The dimensions of the specimen are presented in Fig. 2 and Table 1. The CT specimen was prepared according to the standard procedure given in ASTM E-399. The loading was 19,620 N which gives a time to rupture of about 40 200 hours at a temperature of 550°C. The reference stress level calculated according to [18] is 164 MPa. The fracture toughness of 316 stainless steel is about 240 MPam^0.5 at room temperature and it is reduced by about half at 550°C (cf. [49]). A similar reduction can also be found in [50].

Fig. 2.

Compact tension specimen

The specimen was modelled in the Abaqus Finite Element system in plane strain. Due to symmetry, only half of each specimen was modelled. The example of the mesh with symmetric boundary conditions is presented in Fig. 3.

Fig. 3.

Finite element mesh with symmetric boundary conditions

The material parameters were determined on the basis of uniaxial creep tests at a temperature of 550°C (cf. [14,18]). They are as follows: elasticity modulus E=169617 MPa, Poisson ratio ν=0.3, creep deformation parameters: B=6.561E-24 (MPa)^-n h^-1, n=7.778, damage development parameters: A=3.983E–23 (MPa)^-mh^-1, m=7.622. Additionally, the dependency of instantaneous strain ε₀ was modelled on the basis of plastic deformation and primary creep.

It can be observed in Fig. 4 that the Kachanov eqs. (7–8) assure a good approximation of strains in the secondary creep period but exceed the experimental values for tertiary creep. It affects the strain at failure but does not change the time to failure which is the main target of the analysis.

Fig. 4.

Comparison of results of the Kachanov model with chosen experimental points of uniaxial tension creep of 316 steel at a temperature of 550°C [14]

To model the influence of stress triaxiality, equation (11) is used. Parameter α was determined by Hayhurst et al. [14] in the notched bar tension experiments and the value of α=0.75 was used in the current paper.

FM crack growth simulation for stationary creep

In order to simulate the crack growth, the C*-integral was used. The dependency of C* on crack length was found using finite element analysis. The built-in procedure of Abaqus system (cf. [51]) was used for the determination of the C(t) parameter and then the limit value of it in the stationary state was established as the C* value. The Norton creep model – eq. (8) – without damage development (ω=0) was used as a constitutive equation. The results of this analysis are presented in Fig. 5 (blue line).

Fig. 5.

C*-integral for stationary creep (blue circles), and C(t) for a model with damage (orange squares) as a function of crack length

Next, the value of t_i (crack initiation) was established. Eqs. (4–5) require knowledge of material parameters D and φ. The approximate values of these parameters are defined in eq. (6). The index φ is equal to 0.85 and the value of parameter D is equal to 692.3, assuming a plane strain state and ductility of ε_f=0.13 (cf. [14]). Using these values, the lower estimation of initiation time is 12.6E3 h and the upper estimation is 110.8E3 h for da=0.5mm. The upper estimation is much larger than the experimental value of the time of final failure; therefore, the lower estimation was chosen for future analysis. Such an approach also corresponds to the results presented in [22].

The simulation of crack growth was obtained by the integration of eq. (6). The simplified flowchart of the applied procedure is presented in Fig. 6. The result is presented in Fig. 7 – solid blue line. The time to failure t_f obtained in this method (sum of the initiation time and crack growth time) is 84.7E3 h, which is two times greater than the experimental value. The safety margin according to this method is t_f/t_i=84.7/12.6=6.7.

Fig. 6.

Flowchart of numerical simulation of crack growth

Fig. 7.

Crack growth simulation for stationary creep (solid blue line) and for the creep damage model (dashed orange line) obtained on the basis of FM equations

FM crack growth simulation for the creep damage model

The method described above can be applied when in the vicinity of the crack tip, the stress distribution reflects the stationary state. This is the case when stress redistribution is faster than crack growth. The duration of stress redistribution can be calculated from the formula (cf. [1]):

The second model considers the damage field developing prior to the crack formation. As a consequence, the crack growth is faster and the time to failure shorter. The full set of eqs. (7–9) is used in this case. When the damage parameter reaches its critical value (0.99) at an integration point, this point is marked as damaged and excluded from the problem domain. The cluster of such points forms a crack or a void ahead of the crack tip.

In these conditions, it is impossible to reach the stationary distribution and determine the C* parameter. The C(t) defined by eq. (2) can be used instead (cf. [1,22]). There are not many works using this approach (cf. [52,53]); therefore, a custom procedure based on the scheme presented in Fig. 6 was used. The value of C(t) integral is determined at the moment just before the damage parameter reaches its critical value. The obtained results for different crack lengths are presented in Fig. 5 (orange line). These values are greater than the corresponding C* values, but for small crack lengths, the increment of C(t) is smaller, which generates a more stable rate of crack growth of order 1E-3 mm/h. The change of dependency function can be observed when the crack length achieves about 25 mm, which can be considered as the critical length.

The parameters of eq. (3) are the same as in the previous model. The index φ in the damaged environment is often equal to (cf. [54]):

However, in the examined case, this expression gives a value very close to the approximate value given by eq. (6), which was finally chosen (0.868 vs. 0.85).

The curve of crack growth obtained on the basis of the above assumptions is presented in Fig. 7 (dashed orange line). The time of crack growth up to failure is 10.2E3 h, and total time including the initiation time is 22.8E3 h. This value is about two times smaller than the experimental time and the safety margin obtained in this model is t_f/t_i=1.8. It is much smaller than predicted by the previous analysis and can be used as its lower approximation.

DM local model

In the next step, the local Kachanov model described by eqs. (7–9) was used. The material parameters were listed in Section 4.1. As the solution is strongly mesh-dependent, the simulations were performed for four meshes with a minimum mesh size from 0.5 to 0.05 mm. The times t_i and t_f obtained in these simulations are presented in Fig. 8. Despite the large spread of results, it can be noticed that the safety margin t_f/t_i is more repetitive and it falls between 2.0 and 2.6, except for one case (a very rare mesh). The lower bound of this range can be used as an appropriate approximation of the ratio t_f/t_i.

Fig. 8.

Crack initiation time ti and the time to failure tf as a function of minimum mesh size for local and non-local damage models

For the examined model, the value of time to failure parameter t_f is within the range of 160 h to 11E3 h, indicating a very large spread. All calculated values are smaller than the experimental value. The crack was developing in a plane of initial crack, as it was expected. Its growth rate is initially stationary at the level of 2E-2 mm/h (see Fig. 9). From Fig. 9, the approximate value of critical crack length can also be determined. Its value depends on the mesh density and varies from 20 to 32 mm.

Fig. 9.

Simulation of crack growth by local damage model for different mesh sizes (the time scale for mesh size 0.5 mm is ten times greater than for the other meshes)

DM non-local model

In the non-local model, the state of point depends not only on the previous state of the point but also on the state of neighbouring points. This causes the time to rupture of the cracked specimen to be greater than that obtained in the local method as the factors influencing the damage growth are volume averaged.

Different state variables can be averaged (see [43] for review). The most popular formulation assumes that the damage variable depends on averaged strains (cf. e.g. [45,55,56,57,58]). In the current work, damage development is a function of stress and the damage increment was chosen as the averaged parameter. Such an approach was proposed by Chaboche (cf. e.g. [13,44,59]). First, the local damage increment is calculated:

In the next stage, the equation of damage growth is formulated in terms of the average value of damage growth. It is calculated with use a formula derived from eq. (13):

The application of the non-local method enabled the obtaining of the time to crack initiation t_i which is much less dependent on mesh size than in the local model. Its dispersion was from 1.6E3 to 2.3E3 h. The value of the t_f parameter was much more scattered, from 5E3 to 17E3 h, but it was still more stable than in case of the local method. According to the performed simulations, the estimation of the safety margin is between 2 and 4, which is slightly more than that obtained from previous methods.

By comparison of the growth curves for different mesh sizes, one can notice that they are close to each other in the initial period. Thus, the goal of the application of the non-local method was achieved. Further development of the crack follows by step increments and the differences between individual solutions are observed. This indicates that an unambiguous solution cannot be achieved for time t_f using the presented method. Despite it can be noticed that all solutions give a common starting point for fast crack development, this point can be approximated as 25 mm of crack length. This result aligns very well with the previously obtained solution using FM in the damaged environment.

Fig. 10.

Simulation of crack growth by the non-local damage model for different mesh sizes