The Stress-Sensitivity, Mesh-Dependence, and Convergence of Continuum Damage Mechanics Models for Creep

The long-term exposure of industrial gas turbine (IGT) components such as steam pipe sections, pressure vessels, boilers heat exchangers, and disks at elevated temperature and pressure makes creep damage critically important to consider [1,2]. The growing demand to obtain higher thermal efficiency and reduced CO₂ emissions have led to higher operating temperatures and pressures. There is an increasing need to develop accurate creep deformation, damage evolution, crack initiation, crack propagation, and rupture life prediction models [1].

The effective stress concept of continuum damage mechanics (CDM) can be employed to describe the damage process in materials from crack initiation to rupture [3]. In the CDM approach, a continuous damage variable, $ω$⁠, is coupled with the viscous function of a constitutive model to incorporate the effects of microstructural damage into the constitutive response. The damage variable is assumed to evolve from zero (no damage) to unity (rupture). Kachanov originated the phenomenological CDM approach that was later extended by Rabotnov [4,5].

One of the potentially useful applications of the classical KR model is in crack growth analysis. Finite element analysis (FEA) using KR exhibits mesh shape, distribution, and size dependency such that upon mesh refinement local CDM does not converge to a single solution [6–11]. For KR, damage tends to localize around the crack tip [7,10,12]. Needleman and Tvergaard observed increased crack growth rates upon mesh refinement and proposed a microstructural model that specifies inclusion size, spacing, and initial crack tip radii to mitigate mesh sensitivity [11]. Penny introduced a modified equivalent stress to define critical damage as a function of the yield and ultimate tensile strengths (UTS) and thus match the high-stress region of the stress–rupture curve [13]. Liu and Murakami proposed a new model to mitigate damage localization using a micromechanics approach and a mathematically valid damage evolution equation using the calculus of variations [7], based on the work of Hutchinson [8] and Riedel [9]. The Liu–Murakami model preserves the original KR rupture prediction equation. Rouse et al. illustrated that the rupture prediction of the Liu–Murakami model is linear (on a log–log scale) and suggested that the Dyson Sinh model is superior [1]. An alternative Sinh model similar to the Dyson model has been validated for a wide range of operating stress [14].

While the local CDM approach represents an “old school” approach to the issues of creep crack growth and rupture, and nonlocal CDM approaches have been introduced to mitigate damage localization and mesh dependence [15,16], local CDM continues to be applied in the design against creep of structural components in the power generation industry [17–19]. In some organizations, local CDM has only recently been implemented to replace analytical master curve models, such as Larson–Miller, Monkman–Grant, etc. [20,21]. The MPC omega method, Theta projection, and local CDM are the in-practice constitutive models for simulating long-term (>100,000) creep deformation in complex fossil energy structural components [22]. The international codes, such as ASME B&PV III, French RCC-MR, and British R5, recommend a phenomenological approach to creep and creep-fatigue where damage/life fraction rules, regression analysis, and confidence bands are used to manage reliability and preserve conservatism [22–25]. The advantages and disadvantages of the in-practice models should be investigated and improved models suggested for the practice.

The trigonometric function, Sinh, is commonly employed in constitutive and life prediction models. Yu et al. used the Sinh function to define a stress index to accommodate multiaxial stress [2]. Damage is the void fraction, $fh$⁠, and the proposed void growth rate contains a $(1−fh)$ term in the denominator that may lead to infinite void growth rate as the void fraction, $fh$⁠, approaches unity. Rouse et al. further observed that having a $(1−fh)$ term in the denominator may increase the computing time in FEA and suggested that the Dyson model containing a Sinh function may mitigate this problem [1]. The Dyson model includes a state variable $ϕ$ from $(0<ϕ<1)$ representing aging and another state variable $ω2$ representing cavitation damage. It is reported that $ϕ$ is as low as 0.1 and $ω2=1/3$ for uniaxial cases. These critical damage values are not unity, are a function of stress and temperature, and are subject to scatter such that it can be difficult to produce an accurate prediction of rupture [1,26]. Haque and Stewart introduced an alternative Sinh constitutive model for Waspaloy under triaxial creep conditions [27]. Analytical solutions have been developed to determine the material constants [28]. A comparative and functional relationship study of the MPC Omega, Theta projection, and alternative Sinh model shows that Sinh offers more flexibility and better damage prediction [29,30].

In the current study, the creep deformation, damage, and crack growth of the KR and Sinh models are compared. A detailed explanation and analytical derivation of each model are performed followed by calibration to three isotherms of 304 stainless steel creep data. Using the calculus of variations, the damage variation as a function of stress variation is mathematically exercised. The critical damage calculated at each test condition is compared. Creep crack growth simulations of a thin finite-width plate under uniform tension with a center hole with radial cracks are performed at three mesh sizes. Contour plots of the damage field, bar graphs of the damage distribution on the X- and Y-axis, crack growth rate versus crack length, and time-step versus central processing unit (CPU) time versus simulated time plots are employed to compare the models.

The alloy 304 stainless steel (304SS) is an Fe–Cr–Ni alloy that has a high temperature, oxidation, and corrosion resistance. Experimental data are obtained from a study on the statistical creep properties of 304SS at elevated temperature [31]. Creep tests were repeated five times and performed according to ASTM E139 at 160, 180, 240, 260, 300, and 320 MPa stress from 600 °C to 700 °C.

When the creep strain rate is equal to the minimum creep strain, $ε˙cr=ε˙min=A⋅σn$⁠, the term $(ε˙cr/A)1/n$ of the analytical damage (Eq. (3)) becomes equal to the equivalent stress, $σ$⁠, such that damage (Eq. (3)) becomes zero. As soon as the creep strain rate is greater than the minimum creep strain rate, $ε˙cr>ε˙min$⁠, the term $(ε˙cr/A)1/n$ becomes larger than the equivalent stress, $σ$⁠, and irreversible damage begins.

In this form, the KR damage only depends on the $ϕ$ material constant.

When damage approaches unity, $(ωr→1)$⁠, a small variation of stress, $∂σ(t)$⁠, will create a near infinite damage variation, $∂ω(t)$⁠, when $ϕ$ is greater than zero. This will be graphically determined in Sec. 3.3.

where $ε˙final$ is the final creep strain rate and $ε˙min$ is the minimum creep strain rate calculated directly from creep data. The variable $ω$ represents damage evolution from zero to unity. In the damage rate (Eq. (10)), the material constants M, $ϕ$⁠, $χ$⁠, and $σt$ must be greater than zero.

The Sinh damage is always unity at rupture. Damage begins to accumulate when $ε˙cr>ε˙min$⁠. During steady-state creep, the strain rate, $ε˙cr$⁠, in the numerator of the analytical damage (Eq. (12)) becomes equal to the minimum creep strain rate, $ε˙min$⁠, in the denominator and analytical damage (Eq. (12)) becomes $ω(ε˙)={ln[1]/λ}2/3=0$⁠. When the final creep strain rate is reached, $ε˙cr=ε˙final$⁠, analytical damage (Eq. (12)) reduces to $ω(ε˙final)=( λ/λ)2/3=1$⁠. Thus, for any experimental dataset, the critical damage of Sinh is always unity.

In this form, the Sinh damage only depends on the $ϕ$ material constant.

When damage is critical, $(ωr→1)$⁠, the variation of damage, $∂ω(t)$⁠, does not become infinite and remains finite throughout the damage evolution.

The form of the damage equation and its time derivative are important for the performance of finite element (FE) simulations. This can be realized by considering a fracture simulation where the crack tip is in front of an element as shown in Fig. 1. Around the crack tip, a creep damage zone exists where the damage is distributed from zero to unity $(0≤ω≤1)$⁠. The crack propagates when damage becomes unity. The sharp edge at the crack tip creates a stress concentration. A stress gradient or “variation” exists across the element at the crack tip. A stress variation of $∂σ$ across an element will cause a damage variation of $∂ω(∂σ)$ as depicted in Fig. 1. The damage variation, $∂ω(∂σ)$⁠, for each time step, $∂t$ is determined using the time derivative of the damage equation and added to the previous damage. Stress concentration at the crack tip initiates creep deformation and damage accumulation assessed by the model equation. When damage approaches unity, the KR model (Eq. (8)) develops an infinite damage variation for an infinitesimal change of stress that can lead to localized and time-step sensitive damage growth. The Sinh model (Eq. (14)) does not have a damage term in the denominator, and stress variation does not lead to infinite damage variation.

An analytical exercise of the damage variation equation of KR and Sinh as a function of damage (on a log-normal scale) is depicted in Fig. 2. Using the KR model, an infinitesimal stress variation (for example, $δσ=$ 10 × 10⁻⁶, 10 × 10⁻⁴, or 10 × 10⁻³MPa) will produce a damage variation approaching infinity as damage approaches unity as shown in Fig. 2(a). The horizontal line in Fig. 2(a) indicates where damage variation is equal to unity. The KR model reports a damage variation that is the order of magnitude greater than unity at damage less than unity. Given an infinitesimal stress variation and/or element size, the KR model will always localize. The formed damage variation of KR will cause numerical instabilities. Other investigators have also identified this problem. Liu and Murakami found damage localization using the KR model in near homogeneous stress fields with a very small (10 × 10⁻⁸) stress gradient [7]. Peerling et al. found that as the mesh is refined, the fracture resistance is exceeded and crack growth becomes infinite [15]. Using the Sinh model, an infinitesimal stress variation (for example, $δσ=$ 10 × 10⁻⁶, 10 × 10⁻⁴, or 10 × 10⁻³MPa) will produce a damage variation that exhibits a power law trend and does not explode to infinity as shown in Fig. 2(b). When damage is unity, for an infinitesimal stress variation, $δσ$ produces a small damage variation, $δω$⁠. Given an appropriate element size, the Sinh model will always remain numerically stable. It is evident from the typical fracture simulation illustrated in Fig. 1 and the supported evidence depicted in Fig. 2 that when damage is unity, the KR model becomes numerically unstable, whereas the Sinh model remains stable.

Traditionally, critical damage is set to unity; however, recent investigations and suggestions by literature have demonstrated that critical damage can be assumed to be less than unity and exhibits a dependence on stress and temperature, $ωr=f(σ,T)$⁠, similar to the rupture data, $tr=f(σ,T)$ [7,32]. The obviousness of this relationship extends from Eq. (4). The uncertainty in stress-rupture data would suggest that obtaining accurate critical damage values is a major obstacle. Some reported practices to determine the nonunity critical damage are to

1. (1)

   use pure numerical optimization to determine an optimal, single critical damage value for a given dataset. This approach has the downside of providing no physical justification for the value obtained as well as restricting modeling to interpolation only [14].

2. (2)

   find a functional relationship between critical damage and monotonic tensile properties. In some cases, critical damage is fixed to a constant irrespective to stress and temperature (typically between 0.3 and 0.6) equivalent to the average creep ductility [33]. Unfortunately, creep ductility exhibits the same uncertainty as stress-rupture data. In other cases, critical damage is a function of temperature only. Penny defined critical damage as a function of the yield strength, $σy(T)$⁠, ultimate tensile strength, $σUTS(T)$⁠, and a scaling factor $0≤α(T)≤1$ through the assumption of a critical net stress above which rupture occurs [13]. Unfortunately, all the constants exhibit temperature dependence necessitating a large experimental dataset and careful selection of the regression analysis functions to avoid inflections and physically unrealistic predictions. In the final case, critical damage is a function of stress and temperature and often takes the form of $ωr=f(σ,T)⋅tr$ or $ωr=f(σ,T)⋅εr$⁠. The downfall of this approach is the additional complexity added to the modeling approach.

It would be advantageous to maintain critical damage as unity irrespective to stress and temperature.

where $E0$⁠, $E$⁠, and $Ecr$ are the initial, current, and critical Young's modulus, respectively. When critical damage is less than unity, the critical stiffness can remain large. This creates a numerical problem where the stiffness behind the crack tip needs to be relieved. Some methods to relieve the stiffness are as follows

1. (1)

   allow stiffness to remain behind the crack tip which leads to inaccurate simulations [34],

2. (2)

   step change the stiffness from the critical value $Ecr$ to a very small number which leads to numerical instability where periods of intense time-step bisection are encountered when critical damage is reached in each element [35],

3. (3)

   apply the element deletion approach which violates the conservation of mass [36], or

4. (4)

   apply the cohesive zone approach where the crack surface is separated using weak springs [37]; however, the crack path must be proscribed a priori.

These stiffness reduction methods can be avoided by introducing constitutive models where critical damage is always unity. When critical damage is unity, the stiffness behind the crack tip becomes a moot point.

The KR (Eqs. (1) and (2)) and Sinh (Eqs. (9) and (10)) models are implemented as implicit creep equations with von Mises potential using a simplified Hill's tensor. For both models, the equivalent stress in the damage evolution equation is replaced by the first principal stress, $σ1$⁠. The constitutive equations are implemented into ansys FEA software. The user material routine (USERMAT), an ansys user-programmable feature (UPF), is used to define the constitutive behavior. For every Newton–Rapson iteration and every material integration point, the USERMAT UPF is called. At the beginning of a time increment, the current stresses, strains, and state variables are inputs. The USERMAT must then provide updated stresses, inelastic strains, state variables, and the material Jacobian matrix as outputs [38]. The stress increment is determined using the radial return technique that consider creep as an incompressible process where volumetric stress has no effect on the creep strain and the stress has no effect on volumetric creep strain [39]. The critical damage is restricted to 0.99 instead of unity to prevent a divide by zero error when using the KR model (Eqs. (1) and (2)).

The exact solution to this derivative is available from the authors [28].

At each time step, damage accumulates and the elastic modulus degrades using Eq. (15).

A smooth creep specimen is simulated by simplifying the geometry to a single three-dimensional eight-node element. A uniform load is applied on the top surface of the element. Appropriate displacement constraints were applied to replicate the uniaxial stress observed in a smooth creep specimen. The simulation is run up to the experimental rupture time to facilitate comparison with experimental data. The material constants of 304SS for the KR and Sinh model are listed in Tables 1 and 2, respectively.

Histories of creep deformation, $εcr(t)$⁠, and damage evolution, $ω(t)$⁠, are compared to creep data collected at 700 °C, 650 °C, and 600 °C in Fig. 3. Both the KR and Sinh model are able to predict creep deformation within the boundaries of the repeated creep tests. When compared to the Sinh model, the KR model underpredicts the secondary creep regime and overpredicts the tertiary regime (depicted in the figures as bold drop lines). This under/over prediction is due to the KR equations. The KR creep strain rate (Eq. (1)) term $(σ/(1−ω))n$ increase rapidly at the onset of tertiary creep $(ω→1,n>1,σ>1)$⁠. When the KR model constants are optimized to fit the tertiary creep regime, the secondary creep regime is under predicted.

Both the KR and Sinh models fit the experimentally derived analytical damage (using Eqs. (3) and (12)) within the boundaries of the repeated creep tests. The KR model damage evolves from 0 to 0.4 and then undergoes a dramatic increase in damage rate to unity (solid black line in Figs. 3(b), 3(d), and 3(f)). The Sinh model damage exhibits a continuously increasing damage rate. Thus, the Sinh damage model (short dash black line in Figs. 3(b), 3(d), and 3(f)) follows a more realistic damage path.

Fracture simulations of a thin finite-width plate under uniform tension with a center hole with radial cracks are conducted. The dimensioned geometry of the plate is provided in Fig. 4(a). A uniform load of 33 MPa is applied to the edge of the plate at 700 °C. The Young's modulus at 700 °C is $E0=190$ GPa. A quarter section of the specimen is simulated with displacement constraints shown in Fig. 4(b). Three mesh sizes (0.125 mm, 0.05 mm, and 0.01 mm) are simulated. The mesh size, $Δe$⁠, is defined as the fixed element edge length sets along the crack path. Two-dimensional eight-noded PLANE183 elements with the plane stress option are performed. These simulations are conducted with the first principal stress acting as the driving force behind damage evolution and subsequent crack propagation. The mesh statistics including the number of elements and nodes and the stress concentration factor are provided in Table 3.

Contour plots of the damage distribution along the crack (2 mm from the center of the hole) of the 0.01 mm meshed KR and Sinh simulations are provided in Fig. 5. The KR and Sinh contours are taken at identical crack lengths. The Sinh model exhibits a continuous damage distribution where an equation can be written that describes the contours of damage. The KR model has a discontinuous damage field where inflections appear in the damage contours (as illustrated by the white arrows). As the mesh is refined, the inflection points intensify. From a mathematical perspective, the KR damage field is discontinuous.

An observation concerning the mesh sensitivity of the models is obtained when the length of damage distribution (with ω ranging from 0.11 to 0.99) along the X- and Y-axis is profiled with respect to mesh size as depicted in Figs. 6(a) and 6(b), respectively. Overall the Sinh damage distribution is larger (indicated by the bar height) and reduces at a higher rate as the mesh is refined (indicated by the solid line), while the KR damage distribution is smaller (indicated by the bar height) and reduces at a lower rate (indicated by the dotted line). Both models appear to localize upon mesh refinement; however, Sinh reaches an asymptote (marked in Fig. 6) beyond which mesh-size dependence vanishes. When the damaged area is calculated (the area where ω ranges from 0.11 to 0.99), the damaged area of KR is always smaller than Sinh. For meshes 0.125, 0.05, and 0.01 mm, the KR damage area is 75%, 60%, and 62.5% smaller than Sinh. These findings demonstrate that the Sinh model is less localized than KR.

An observation concerning the numerical instability of the models is obtained when the damage with 0.05 mm mesh are compared for the KR and Sinh in Fig. 7. In the model, the mesh size increases from 0.05 mm in the P region to 0.125 mm in the Q region. The mesh is shown in Fig. 4(b). The simulations are run to rupture where numerical instability due to the loss of stiffness prevented further convergence. The KR model became numerically unstable earlier than Sinh as indicated by a shorter crack length. The KR model exhibits a sudden change in damage distribution at the transition point from P to Q as circled in Fig. 7(a). Once, in the Q region, the damage field becomes tortuous and unstable. The Sinh model exhibits a smooth damage field with no sudden change in the damage distribution upon mesh size increase as shown in Fig. 7(b).

The crack growth rate using the KR model is provided in Fig. 8(a). For KR, as the mesh size decreases, the amplitude and frequency of oscillations in the crack growth rate increase. When examining the crack growth rate at initiation, the initial crack growth rate is inconsistent between the three meshes having no discernable trend. In addition, the initial crack growth rate is not the lowest observed with a drop observed between 0 and 1 mm of crack length. The cause is found analyzing the KR equations. When $ω→1$⁠, the KR creep strain and damage rate (Eqs. (1) and (2)) and damage variation (Eq. (8)) become infinite. The high stress and damage variation at the crack tip contribute to KR producing numerically unstable crack growth. As the mesh is refined, the stress concentration increases. This contributes to the inconsistencies and increasing oscillations upon mesh refinement.

The crack growth rate using the Sinh model is provided in Fig. 8(b). For Sinh, as the mesh size decreases, the amplitude and frequency of oscillations in the crack growth rate decrease. When examining the crack growth rate at initiation, the initial crack growth rate decreases upon mesh refinement maintaining a discernable trend. If a polynomial was fit to each mesh size, the deviations between the three polynomials will be minimum. As such, a single polynomial could be used to describe the crack growth rate irrespective of mesh size. The accuracy of the crack growth rate using Sinh is not mesh-sensitive for the given mesh sizes. When looking closely at the decreasing oscillations with mesh size, only the precision of the crack growth rate is mesh-sensitive.

In summary, the Sinh crack growth rate is smoother and more consistent than KR. The Sinh model produces self-accurate predictions of the crack growth rate that become more precise with mesh refinement. The maximum representative volume element (RVE) size of Sinh is larger than the RVE of KR. The increasing oscillations in the crack growth rate of KR upon mesh refinement make it difficult to determine an appropriate RVE for KR. The KR crack growth rate becomes numerically unstable while the Sinh model improves in stability and precision upon mesh refinement.

A good metric to analyze mesh sensitivity and convergence is the evolution of time-step size and CPU time with respect to simulated time at various mesh sizes as illustrated in Fig. 9. The time-step is the increment of simulated time taken during each iteration in a simulation. The initial time-step was set to 1 × 10⁻⁶ h. The time-step size increases until the maximum time-step size of 5 h is achieved. ansys bisects the time-step when numerical instabilities prevent convergence. The CPU time indicates the real computer time required to complete a simulation. If bisection occurs, more iterations are required to reach simulated rupture life and the CPU time will increase.

The Sinh model exhibits less time-step bisection and CPU time does not increase dramatically upon mesh refinement when compared to KR as presented in Fig. 9. For both models, time-step bisection is required to accommodate the high-stress near-instantaneous rupture at the end of life. The KR model bisects earlier than the Sinh model with bisections occurring at 62% and 13% of simulated rupture life in 0.05 and 0.01 mm meshes, respectively. The Sinh model experiences bisection at 89% and 65% of simulated rupture life in 0.05 and 0.01 mm meshes, respectively. The early bisection of the KR model contributes to a larger CPU time as the mesh is refined. For a 0.05 mm mesh, the CPU time required for the KR and the Sinh model is almost identical up to 5000 h of simulation time before deflection as shown in Fig. 9(a). Upon mesh refinement to 0.01 mm, the KR model becomes unstable earlier such that the CPU time required for the KR and the Sinh model is only identical up to 1000 h of simulation time as shown in Fig. 9(b). This is indicative of multiple bisections occurring at each iteration until an appropriate time-step is obtained for continued convergence. The failed-bisected-iteration attempts to increase the overall CPU time. ansys always attempts to return to the maximum time-step size after bisection. This can contribute to additional failed iteration attempts as the time-step size continuously oscillates from the maximum to the optimal. Overall, it is determined that the Sinh model requires less bisection and thus less CPU time.

Further observations concerning mesh sensitivity can be elucidated by examining the simulated rupture time as depicted in Fig. 10(a). Upon refinement from 0.05 to 0.01 mm, the KR rupture time increased by 20% (1079 h) while Sinh only increased by 0.6% (32 h). The root cause of this phenomenon involves the damage variation as a function of stress variation when damage approaches unity illustrated in Fig. 2(a) and depicted in Fig. 2(b).

The convergence of KR and Sinh can be examined by looking at the total CPU time required to reach rupture depicted in Fig. 10(b). Upon mesh refinement from 0.05 to 0.01 mm, the required CPU time increased by a factor of 10.86 and 5.50 for KR and Sinh, respectively. Logically, as mesh size decreases, the computational costs increase, resulting in more CPU time; however, bisections and failed-bisection attempts can dramatically increase CPU time. It is determined that the Sinh model exhibits better convergence than the KR model.

It is concluded that the Sinh model offers better creep damage and crack growth analysis. When the Sinh model is compared to the KR model, it is observed that

1. (1)

   The Sinh model is normalized by the experimental final creep strain rate (Eq. (12)) such that critical damage is unity irrespective of stress and temperature conditions while the KR critical damage (Eq. (3)) cannot be unity as the experimental final creep strain rate is not infinite

2. (2)

   The Sinh model exhibits a more distributed damage field when compared to KR model. The Sinh crack growth rate becomes smoother upon mesh refinement, while the KR crack growth rate becomes numerically unstable. Mathematically, it is proven that when damage approaches unity $(ω→1)$⁠, for a small variation in stress, the Sinh damage variation remains finite (Eq. (14)), while KR develops a near infinite damage variation (Eq. (8)) leading to numerical instability.

3. (3)

   The Sinh model is less stress-sensitive, less mesh-dependent, and exhibits better convergence when compared to the KR model. The discontinuous damage of the KR model creates time stepping issues leading to longer CPU time upon mesh refinement. The KR model simulations have difficulty reaching convergence as the crack grows and the stress concentration intensifies.

In future work, an analytical solution will be derived to transform the KR material constants into Sinh constants.

This material is based upon work supported by the Department of Energy, National Energy Technology Laboratory under Award No. DE-FE0027581.

$A,n$ =

Norton power law constants (%MPa⁻ⁿ h⁻¹, dimensionless)

$A,σs$ =

Mcvetty creep law constants (%1/h, MPa)

$Ci$ =

stiffness tensor, $i=$ type

$E$ =

Young's modulus (MPa)

M =

Hill compliance tensor

$M,σt,χ,ϕ,λ$ =

Sinh model constants (h⁻¹, MPa, dimensionless)

$M,χ,ϕ$ =

Kachanov–Rabotnov model constants (MPa^−χ h⁻¹, dimensionless)

s =

Cauchy stress vector

$Si$ =

compliance tensor, $i=$ type

$tr$ =

rupture life (h)

$α$ =

scaling factor

$ε˙cr$ =

creep strain rate (%1/h)

$ω$ =

analytical damage (dimensionless)