Microstructurally Based Modeling of Creep Deformation and Damage in Martensitic Steels

Abstract

This chapter deals with modeling the microstructural evolution, creep deformation, and pore formation in creep-resistant martensitic 9–12% Cr steels. Apart from the stress and temperature exposure of the material, the input parameters for the models are as-received microstructure and one single-creep experiment of moderate duration. The models provide predictive results on deformation rates and microstructure degradation over a wide stress range. Due to their link to the underlying fundamental physical processes such as classical nucleation theory, Gibbs energy dissipation, climb, and glide of dislocations, etc., the models are applicable to any martensitic steel with similar microstructure to the presented case study. Note that we section the chapter into part 1: creep deformation and part 2: pore formation.

1. Introduction

When trying to simulate deformation rates during creep, one can basically follow two different approaches: (i) phenomenologically based—and or (ii) physically inspired models. In phenomenological models, the deformation (rate) is usually stated as analytical function with system parameters such as temperature and stress as input [1, 2, 3, 4]. These kinds of models are easily and quickly employed; however, the approaches carry some disadvantages as well: (1) they give little or no insight into the actual underlying physical processes governing the creep rate, and (2) the model parameters cannot be usually determined independently from the creep experiment the model is actually aiming to predict. For these reasons, we choose to focus on a physically based model instead [5], which is a reviewed, corrected, and extended version of the seminal work of Ghoniem [6]. In addition to avoiding the mentioned drawbacks of phenomenological approaches, our physical model has the advantage of including a variety of microstructural elements such as dislocations and subgrain boundaries (SGBs), and their interactions. This allows for getting a deeper understanding of the creep process and opens the opportunity of rating individual material badges by taking their as-received microstructure as starting condition for a creep simulation. Finally, we end up with an assembly of rate equations for the microstructural elements along with some side equations modeling the physical processes. In summary, the model gives us insights into the specific reasons why a material badge features good or bad creep behavior, as long as its microstructure can be considered homogeneous. The model also has the potential to rate the impact of individual microstructural phenomena.

Our model includes the microstructure by mean values of specific microstructural elements (e.g. dislocation density, grain boundary precipitates, etc.) instead of a spatially resolved features. This allows for a simpler construction of a “representative volume element,” which speeds up computation and is more easily compared to microstructural investigations.

In our work, we focus on martensitic 9–12% Cr steels. We select the material P91 to demonstrate the validity of the model due to its widespread use and industrial significance. Nevertheless, the concept can be adapted to other material groups by including their specific microstructural elements and their interactions.

Please note that all symbols used in the equations are explained in Table 6 at the end of the chapter.

2. Creep model: our microstructurally based approach

The aim of our creep model is to predict the creep rate and microstructural evolution, based on the initial microstructure and the system parameters stress and temperature. A very simple, yet useful approximation of the creep rate $\dot{\epsilon }$ has been introduced by Orowan [7] and later extended by Yadav [8] to include damage. While Yadav considers damage due to cavities D_cav and precipitates D_ppt, damage in our model is based on cavities alone, since the effect of precipitates is dealt with other equations in the framework in a physical manner. The resulting modified Orowan equation calculates the creep rate from the physical inputs of mobile dislocation density ${\rho }_m$ and an effective dislocation velocity ${\upsilon }_{\mathit{eff}}$; see Eq. (1). Latter term has been introduced by Riedlsperger [5] incorporating glide and climb processes. Since both quantities, ${\rho }_m$ and ${\upsilon }_{\mathit{eff}}$, result from interactions within the microstructure, we take a closer look at those, see Figure 1. Within Figure 1, we indicate following interactions:

1. Frank—read sources (generation of mobile dislocations ${\rho }_m$)

2. Emission from static dislocations ${\rho }_s$into mobile dislocations ${\rho }_m$

3. Immobilization of mobile dislocations ${\rho }_m$—generation of static dislocations ${\rho }_s$

4. Recovery by climb processes for mobile (d₁) and static dislocations (d₂)

5. Spontaneous annihilation of mobile (e₁) and static dislocations (e₂)

6. Subgrain boundaries produced from static dislocations

7. Subgrain growth minus Zener—pinning of boundary precipitates

Figure 1.

Microstructural interactions within a subgrain. See text for individual interactions. “X” signifies annihilation of dislocations and arrows mark transformations.

These interactions (a–g) are also integrated into the rate equations for the microstructural evolution of mobile dislocations ${\rho }_m$, static dislocations ${\rho }_s$, dislocations within subgrain boundaries ${\rho }_b$, and subgrain boundaries $R_{\mathit{sgb}}$; see Eqs. (2)–(5) in Table 1. Table 1 also indicates the original source for each equation.

Creep model equations	Eq.
Creep strain rate (modified Orowan equation, see text):	(1)
$\frac{\mathit{d\epsilon }}{\mathit{dt}}=\frac{b\bullet {\rho }_m\bullet {\upsilon }_{\mathit{eff}}}{M\bullet \left(1-D_{\mathit{cav}}\right)}$
Mobile dislocation density rate [5]:
$\frac{{\mathit{d\rho }}_m}{\mathit{dt}}=v_{\mathit{eff}}\bullet \left\{\underset{a}{\underbrace{{\rho }_m^{3/2}}}+\underset{b}{\underbrace{\frac{\beta {\rho }_s{R}_{\mathit{sgb}}}{h_b^2}}}-\underset{c}{\underbrace{\frac{{\rho }_m}{2R_{\mathit{sgb}}}}}-\underset{d1}{\underbrace{8{\rho }_m^{3/2}\frac{{\upsilon }_c}{{\upsilon }_{\mathit{eff}}}}}-\underset{e1}{\underbrace{{\delta }_{\mathit{anh}}\left({\rho }_m+{\rho }_s\right){\rho }_m}}\right\}$	(2)
Static dislocation density rate [6]:
$\frac{{\mathit{d\rho }}_s}{\mathit{dt}}=v_{\mathit{eff}}\bullet \left\{\underset{c}{\underbrace{\frac{{\rho }_m}{2\bullet R_{\mathit{sgb}}}}}-\underset{d2}{\underbrace{8\bullet \frac{{\upsilon }_c}{v_{\mathit{eff}}}\bullet \frac{{\rho }_s}{h_b}}}-\underset{e2}{\underbrace{{\delta }_{\mathit{anh}}\bullet {\rho }_s\bullet {\rho }_m}}\right\}$	(3)
Boundary dislocation density rate [6]:
$\frac{{\mathit{d\rho }}_b}{\mathit{dt}}=\underset{f}{\underbrace{8\bullet \left(1-2\zeta \right)\bullet {\rho }_s\bullet \frac{{\upsilon }_c}{h_b}}}-\underset{g}{\underbrace{M_{\mathit{sgb}}\bullet P_{\mathit{eff}}}}\bullet \frac{{\rho }_b}{R_{\mathit{sgb}}}$	(4)
Subgrain growth rate (h: nucleation of subgrains) [5]:
$\frac{{\mathit{dR}}_{\mathit{sgb}}}{\mathit{dt}}=\underset{g}{\underbrace{M_{\mathit{sgb}}\bullet P_{\mathit{eff}}}}-\underset{h}{\underbrace{G\bullet {\eta }_V\bullet K_c^2\bullet \left(\sqrt{{\rho }_m+{\rho }_s}-\frac{K_c}{2R_{\mathit{sgb}}}\right)\bullet \frac{D_s\mathrm{\Omega }}{2k_{B}T}}}$	(5)
Effective subgrain growth pressure (only precipitates at boundaries) [6]:
$P_{\mathit{eff}}=\frac{4}{3}\bullet G\bullet b^2\bullet {\rho }_b\bullet \left[1-\frac{\pi }{2}\bullet \left({\sum }_{i=1}^m{r}_{p,i}^2\bullet N_{V,i}\right)\bullet R_{\mathit{sgb}}\right]$	(6)
Subgrain boundary mobility – case 1, if $P_{\mathit{eff}}\ge 0$ [6]:
$M_{\mathit{sgb}}=\frac{2\bullet \pi \bullet {\eta }_v\bullet D_s\bullet \mathrm{\Omega }}{b\bullet k_B\bullet T\bullet {\theta }_i{f}_i}+\frac{2\bullet \pi \bullet b\bullet D_{\mathit{vp}}\bullet \mathrm{\Omega }\bullet {\theta }_i{f}_i}{h_b^2\bullet k_B\bullet T}$	(7)
Subgrain boundary mobility – case 2, if $P_{\mathit{eff}}<0$ [6]:
$M_{\mathit{sgb}}=\frac{D_s\bullet \mathrm{\Omega }}{2\bullet \pi \bullet k_B\bullet T\bullet {\sum }_{i=1}^m{r}_{p,i}^4\bullet N_{V,i}\bullet {\theta }_i{f}_i}$	(8)
Mean dislocation spacing within the subgrain boundary [6]:
$h_b=\frac{1}{\left({\rho }_b+{\rho }_s\right)\bullet R_{\mathit{sgb}}}$	(9)
Effective dislocation velocity (only precipitates in grain interior) [5]:
$\frac{1}{v_{\mathit{eff}}}=\frac{1}{v_g}+{\sum }_i\frac{\pi }{2}\bullet N_{V,i}\bullet r_{p,i}^3\bullet \frac{1}{v_c}$	(10)
Dislocation glide velocity [9]:
${\upsilon }_g=a_1\bullet \mathit{exp}\left\{-\frac{Q}{k_B\bullet T}\right\}\bullet \mathit{exp}\left\{-\frac{{\sigma }_i\bullet V_r}{k_B\bullet T}\right\}\bullet 2\bullet \mathit{sinh}\left\{\frac{{\sigma }_{\mathit{app}}\bullet V_r}{k_B\bullet T}\right\}$	(11)
Corrected applied stress [5]:
${\sigma }_{\mathit{app}}={\sigma }_{\mathit{app},0}\bullet \left(1+\nu \bullet \epsilon \right)$	(12)
Dislocation climb velocity [6]:
${\upsilon }_c={\upsilon }_{\mathit{cl}}+{\upsilon }_{\mathit{cp}}$	(13)
Dislocation climb velocity-lattice diffusion share [5]:
${\upsilon }_{\mathit{cl}}=\frac{2\bullet \pi \bullet {\eta }_v\bullet D_s}{\left[1-{\eta }_v\bullet ln\left(L_{\alpha }\bullet \sqrt{{\rho }_t}\right)\right]\bullet b}\bullet \mathit{exp}\left\{-\frac{{\sigma }_i\bullet \mathrm{\Omega }}{k_B\bullet T}\right\}\bullet 2\bullet \mathit{sinh}\left\{\frac{{\sigma }_{\mathit{app}}\bullet \mathrm{\Omega }}{k_B\bullet T}\right\}$	(14)
Parameter for elastic interactions between dislocations/defects [6]:
$L_{\alpha }=\frac{\left(1+\nu \right)\bullet G\bullet b\bullet 0.5\bullet \mathrm{\Omega }}{3\bullet \pi \bullet \left(1-\nu \right)\bullet k_B\bullet T}$	(15)
Dislocation climb velocity-pipe diffusion share [5]:
${\upsilon }_{\mathit{cp}}=\frac{2\bullet \pi \bullet b\bullet D_{\mathit{vp}}}{L_p^2}\bullet \mathit{exp}\left\{-\frac{{\sigma }_i\bullet \mathrm{\Omega }}{k_B\bullet T}\right\}\bullet 2\bullet \mathit{sinh}\left\{\frac{{\sigma }_{\mathit{app}}\bullet \mathrm{\Omega }}{k_B\bullet T}\right\}$	(16)
Distance of core-vacancy diffusion before evaporation into lattice [10]:
$L_p=\sqrt{2}\bullet a_g\bullet \mathit{exp}\left\{\frac{∆W}{2\bullet k_B\bullet T}\right\}$	(17)
Internal stress [11]
${\sigma }_i=\alpha \bullet M\bullet G\bullet b\bullet \sqrt{{\rho }_m+c_s\bullet {\rho }_s}$	(18)
Damage parameter for cavities, see text:
${\dot{D}}_{\mathit{cav}}=\mathrm{A}\bullet \epsilon \bullet {\dot{\epsilon }}^{1/2}$	(19)
Particle coarsening:
$r_i^3\left(t\right)-r_{i,t=0}^3={\mathrm{k}}_i\bullet \mathrm{t}$	(20)
Reduction of number density of particles due to coarsening:
$\frac{N_i\left(t\right)}{N_{i,t=0}}=\frac{r_{i,t=0}^3}{r_i^3\left(t\right)}$	(21)

Table 1.

Creep model equations and their primary literature sources.

In addition to the “rate equations” of the microstructural evolution of the material, Table 1 also assembles the framework of the underlying physical phenomena. Within this paragraph, we only give a brief overview. Detailed discussions can be found in the cited sources.

The effective subgrain growth pressure P_eff, Eq. (6), is rewritten, but it is equivalent to Ghoniem’s work [6] and includes dissipation of grain boundary energy as well as Zener pinning from precipitates. The mobility of the subgrains, M_sgb, Eqs (7) and (8), is governed mostly by diffusion coefficients and the boundary misorientation [6]. Dislocation spacing within a subgrain boundary, h_b, Eq. (9), has been deducted by geometrical means [6]. The effective velocity of mobile dislocations $v_{\mathit{eff}}$ includes glide processes as well as climbing of dislocations over precipitates within the subgrain interior [5]. The glide velocity of mobile dislocations ${\upsilon }_g$considers forward and potential backward movements according to their jump probabilities, which are linked to mechanical and thermal activation [9]. The corrected applied stress ${\sigma }_{\mathit{app}}$, as shown in Eq. (12), considers the reduced cross section of a creep sample due to its poisson ratio [5]. Since the climb of mobile dislocations depends on local diffusion, the climb velocity, ${\upsilon }_c$, as shown in Eq. (13) [6], is split into a lattice diffusion share, ${\upsilon }_{\mathit{cl}}$, as shown in Eq. (14) [5], and a pipe diffusion share, ${\upsilon }_{\mathit{cp}}$, as shown in Eq. (16) [5], with additional terms regarding the vacancy-dislocation interactions from Hirth and Lothe [10]; see Eq. (17). The internal stress ${\sigma }_i$, which is an important input for the glide and climb velocities and considers interactions of mobile-mobile and mobile-static dislocations, is taken from Basirat [11]; see Eq. (18). Finally, the phenomenological factor “cavitation damage” D_cav, as shown in Eq. (19), has been introduced by Basirat [11] as well but slightly adapted here for a better agreement with the observed master creep curve.

3. Creep model: experimental input and simulation setup

The input parameters needed to simulate creep over a range of stresses can be divided into three groups/types:

1. General material data from the material group;

2. Microstructural data in as-received condition of the material badge of interest;

3. Parameters which have to be adapted to one single master-creep-curve.

Group (i) is applicable to a wider range of materials and has already been collected in [5]. Table 2 gives an overview on the findings.

Input	Value	Details	Source
ag	2.866 × 10^−10 m	Fe-bcc	[12]
$\alpha$	0.02	9Cr-1Mo steel at 650°C	[11]
b	2.48 × 10^−10 m	Fe-bcc	[13]
cs	0.3	9Cr-1Mo steel at 650°C	[11]
${\delta }_{\mathit{anh}}$	5b	Edge dislocations	[11]
Ds	2 × 10^−19 m^2/s	Fe in Fe-1.5Mo at 650°C	[14]
Dvp	4.75 × 10^−19 m^2/s	P91	[5]
${\eta }_V$	2 × 10^−4	Ferritic steels	[6]
G	61.9 GPa	9–12% Cr steels at 650°C	[15]
kB	1.38065 × 10^−23 J/K	Fundamental constant	[16]
Kc	2.1	P91	[5]
M	3	Crystals	[11]
$\Omega$	1.177 × 10^−29 m	Volume/atom	[5]
Q	4.01 × 10^−19 J	HT-9 at 550°C	[6]
${\theta }_i{f}_i$	3°=0.0524 rad	9–12% Cr steels	[17]
Vr	35$\mathrm{\Omega }$	P91 at 650°C	[5]
$\nu$	0.317	9–12% Cr steels at 650°C	[15]
$∆W$	1.26 × 10^−19 J	P91 at 650°C	[5]
$\zeta$	0.034	9–12% Cr steels	[18]

Table 2.

Input parameters for the material group. Input parameters type (i).

Group (ii) is the parameter for a specific material badge. This group accounts for the different creep behavior stemming from specific processing routines, e.g. chemical composition and heat treatment. Result from the processing routine is the as-received microstructure, which also acts as an input for our simulation. In our case, we directly measured the subgrain size R_sgb and the boundary dislocation density ${\rho }_b$via electron backscatter diffraction (EBSD) [5]. The typical as-received mobile dislocation density ${\rho }_m$was taken from Panait [19]. Simulations reveal that ${\rho }_m$ quickly converges toward generic values during creep, so the results are not very sensitive to the initial value of this parameter. The same is true for${\rho }_s$, which has been adopted from [11]. The precipitate data have to be split into particles at subgrain boundaries (SGBs) and particles within subgrains (VN, NbC, and AlN). Both have been calculated by the thermodynamic software MatCalc [20] with the chemical composition of the heat, the actual heat treatment, and matrix data [5]. The results of the particle simulation, coarsening of particles from as-received condition onward, have been mimicked by a simple coarsening law, as shown in Eqs. (20) and (21), which is now used as input for the creep simulation. Table 3 summarizes the microstructure-specific input data for the as-received condition, which are the starting values of our creep simulation.

Input	Value	Note	Source
Rsgb	4.0 × 10^−7 m	EBSD	[5]
${\rho }_b$	5.9 × 10^14 m^−2	EBSD	[5]
${\rho }_m$	4.5 × 10^14 m^−2	TEM	[19]
${\rho }_s$	4.5 × 10^13 m^−2	Estimate 1/10 of ${\rho }_m$	[11]
ri (SGB)	1.49 × 10^−7 m	At boundaries	[5]
ri (VN)	4.17 × 10^−8 m	Grain interior	[5]
ri (NbC)	7.55 × 10^−8 m	Grain interior	[5]
ri (AlN)	4.43 × 10^−8 m	Grain interior	[5]
Ni (SGB)	9.52 × 10^18 m^−3	At boundaries	[5]
Ni (VN)	1.09 × 10^20 m^−3	Grain interior	[5]
Ni (NbC)	4.04 × 10^18 m^−3	Grain interior	[5]
Ni (AlN)	8.85 × 10^18 m^−3	Grain interior	[5]
ki (SGB)	6.90 × 10^−28 s^−1	At boundaries	[5]
ki (VN)	2.21·10^−30 s^−1	Grain interior	[5]
ki (NbC)	0 s^−1	Grain interior	[5]
ki (AlN)	1.19 × 10^−29 s^−1	Grain interior	[5]

Table 3.

Microstructure-specific input data for the as-received state. Input parameters type (ii).

After these parameters have been set, only the variables a₁, A, and $\beta$ remain (group (iii) of input parameters), whereas A is of phenomenological nature and thus cannot be measured directly as a matter of principle, and the parameters a₁ and $\beta$ have physical interpretations but are extremely inconvenient for direct assessment. We therefore have to set up one single (!) creep experiment at set temperature T and stress ${\sigma }_{\mathit{app},0}$ to fit these parameters.

We conducted/performed an instrumented creep test at 650°C and a nominal stress of 70 MPa which led to a rupture time of 8740 h [5]. The simulated creep curve has been validated against the experimental results and the missing parameters a₁, A, and $\beta$ were optimized using a least-squares fit of the time-dependent creep deformation. Table 4 summarizes the results which are also used for all other creep simulations presented in the following text.

Parameter	Value
a1	3.725 m/s
A	0.035
$\beta$	0.07

Table 4.

Parameters found by fitting against one creep curve. Input parameters type (iii).

Since the creep model already contains the impact of stress and temperature explicitly within its network of equations, the found input parameters are stress- and temperature-independent. We can thus use the same input for other stresses to produce time-to-rupture (TTR) diagrams. To do so, we have left all input parameters unchanged as indicated in Tables 2–4 and generated creep simulations in the stress range of 50–120 MPa. The next section comprises the result of the simulated master creep curve including the calculated microstructural evolution, as well as the changes of the creep behavior with altering stresses, leading to the construction of the TTR diagram.

4. Creep model: simulation results and discussion

First simulated result is the master-creep curve at 650°C and 70 MPa, indicating the creep deformation and deformation rate. Figures 2 and 3 demonstrate the agreement between the simulated and the experimental result of the creep deformation: the simulated creep curve is very close to the experiment including primary, secondary, and tertiary creep stage, and also the final fracture of the sample. The experimental minimum of the creep strain rate is about 2.5 × 10⁻⁶ h⁻¹ in the range between 1.000 and 3.000 h, whereas the simulated result is 3.0 × 10⁻⁶ h⁻¹ at 900 h, suggesting good agreement as well.

Figure 2.

Experimental and simulated master-creep curve at 650°C, 70 MPa.

Figure 3.

Experimental and simulated creep-strain rates at 650°C, 70 MPa.

Regarding the microstructural evolution, the simulation predicts a quick recovery of the mobile dislocation density ${\rho }_m$from 4.5 × 10¹⁴ m⁻² to 1.5 × 10¹³ m⁻² within the first 500 h of creep, which is exactly mirroring the continuous decrease in the creep rate during the primary creep regime. According to the (modified) Orowan equation (Eq. (1)), two potential reasons can be responsible for the creep rate: the mobile dislocation density and the effective velocity of the dislocations. Our model demonstrates the dominating role of the dislocation density in this regard. After reaching its minimum, ${\rho }_m$ increases slowly up to about 2 × 10¹³ m⁻² at the time of fracture. This level agrees well with several literature data on similar material under comparable conditions, indicating dislocation densities at time of fracture of 2 × 10¹³ m⁻² [21] or 2.7 × 10¹³ m⁻² to 3.5 × 10¹³ m⁻² [22].

${\rho }_b$ and R_sgb were verified using EBSD on the fractured sample of the master-creep experiment, showing very good agreement with the simulated data, whereas simulation indicates a subgrain radius of 0.95 μm at the time of fracture and experimental verifications reveal a mean value of 0.7 μm [5]. The experimental value of boundary dislocation density of 3.4 × 10¹⁴ m⁻² also compares very well to the simulated result of 3.0 × 10¹⁴ m⁻² (Figure 4) [5].

Figure 4.

Simulated microstructural evolution at 650°C, 70 MPa.

In summary, the simulation results of the test case of the master-creep experiment are good enough to apply the model to multiple stresses. We carried out creep simulation in the stress range of 50–120 MPa in steps of 10 MPa, with a resulting creep curve and accompanying microstructural evolution in each simulation. One side result, the rupture time t_R, could then be used for reconstructing a time-to-rupture (TTR) diagram. In the simulation, the sample ruptures when the damage parameter D_cav reaches a level of 1, which is basically a result of extensive creep strain and/or strain rate [5]. This typically occurs at a strain between 3% (low stresses) and 10% (high stresses). Please note that the fracture elongation is not an input, but a result of the simulation. Please also note that the simulation does not consider local necking but only deals with sample areas of uniform cross sections. Figure 5 shows the individual simulated creep curves within the investigated stress range.

Figure 5.

Simulated microstructural evolution at 650°C, 70 MPa.

Please note that the creep curves appear to look different from Figure 2 because the logarithmic time-scale is used in order to simultaneously show all results. Each of the creep curves feature primary, secondary, and tertiary creep regimes. Figure 6 now finally shows the constructed TTR stemming from the simulated creep data and compares them against the standard literature data from European Creep Collaborative Committee (ECCC) [23], ASME [24], and NIMS [25].

Figure 6.

Simulated TTR diagram compared against standard literature data.

The agreement is excellent, and the simulated curve lies right between the data from the three standard literature sources for creep rupture data of P91. Once again, it is important to mention that all model input data (except for the system stress) were identical for all creep simulations. This detail is very important, because the simulation allows for a predictive extrapolation from a single-creep experiment carried out for 8.740 h to up to six times longer creep times. As it appears, the simulation also allows for extrapolating to shorter running times by a factor of about 40 in our case. We thus motivate to use and test our model for even shorter reference experiments in order to extrapolate for long running times.

5. Pore formation model: introduction

It has been long established [26] and is now well accepted [27, 28] that failure during creep loading results mainly due to intergranular rupture. Cavities nucleate predominantly at grain boundaries, grow during creep exposure, and coalesce to form microcracks. In tertiary creep, these cracks are so numerous that they significantly weaken the microstructure, and the remaining available cross section is put under more stress which further promotes damage and accelerates the strain rate.

In some cases, the remaining creep life can be directly correlated with the degree of cavitation [29, 30].

6. Pore formation model: nucleation and growth

It is not well established by what mechanism cavities nucleate [31]; however, the linear relationship between cavity nucleation rate, first observed by Needham et al [32], still holds true to this day [33]. Grain boundary sliding, as necessitated to maintain contact between the grains when they elongate during diffusional creep, is one proposed mechanism [34]. This sliding generates cavities at ledges that are pulled apart at the grain boundary.

We propose a model based on the physics of diffusion and fluctuational theory, known as classical nucleation theory (CNT). CNT was formulated at the beginning of the twentieth-century by the works of Volmer and Weber [35], Becker and Döring [36], Frenkel [37], and Zeldovich [38]. It has been prominent and successful in modeling the nucleation of new phases, precipitates, and similar phenomena.

Balluffi [39] was the first to explain the nucleation of holes (cavities) by vacancy supersaturation. However, Raj and Ashby [40] were the first to consider the mechanical stress as the driving force for nucleation, a theory which was further developed by Hirth and Nix [41] and Riedel [28] and forms the basis for our nucleation model.

While CNT generally speaks of nuclei, which may form new phases, we specify these as clusters of vacancies which may form cavities.

The free energy change on formation of such a cluster in the bulk encompasses the pressure-volume work done by the external stress, σ, and the energy required by the newly formed surface between the cluster and the matrix, excluding the elastic energy [42]. Eq. (22) shows the relation between the free energy change, ΔF, and the cluster radius, r, with the specific surface energy, γ:

$$\mathit{\Delta F}=-\frac{4}{3}\pi r^3\sigma +4\pi r^2\gamma$$E22

Plotting this relation over the cluster radius, as shown in Figure 7, we see that the free energy reaches a maximum at a certain cluster size, r*, designated as the critical radius. This energy barrier needs to be overcome for a cluster to become a stable cavity which will then continue to grow as its total free energy decreases.

Figure 7.

Free energy change vs. radius for a spherical cluster of vacancies.

From Eq. (22), we derive the critical radius and critical (maximum) free energy as follows:

$$r^{\ast }=\frac{2\gamma }{\sigma }$$E23

$$∆F^{\ast }=\frac{16}{3}\pi \frac{{\gamma }^3}{{\sigma }^2}$$E24

Clusters below the critical size are naturally/always present in the microstructure due to thermal fluctuations [35, 43]. Their concentration is determined by the number of possible nucleation sites, N_S, and an Arrhenius term comprised of the free energy change over the atomic thermal energy, kT. These are the first two terms in Eq. (25) for the nucleation rate, and they quantify the number of critical cavities available in quasi-static equilibrium. The next term, β^*, is the vacancy attachment rate and Z, is the Zeldovich factor:

$$I=N_s\mathit{exp}\left(-\frac{∆F^{\ast }}{\mathit{kT}}\right){\beta }^{\ast }Z$$E25

Some vacancies, which exist throughout the microstructure and are more prevalent at higher temperatures, may find themselves on the surface of a critical cluster and only one atomic jump away from joining it. The number of these vacancies jumping toward the critical cluster per unit time is described by the vacancy attachment rate, β^*. In Eq. (26), this is shown to depend on the diffusion coefficient, D, the concentration of vacancies, X_V, and the surface area of a critical cluster, A^*, with a being the average interatomic distance:

$${\beta }^{\ast }=\frac{D}{a^2}{X}_v\frac{A^{\ast }}{a^2}$$E26

The Zeldovich factor is explained by its namesake [38] and other literature [28, 44] to reduce the nucleation rate, since steady-state nucleation artificially removes supercritical clusters and because slightly supercritical clusters are still more likely to dissolve rather than grow. It is defined in Eq. (27) with n^* signifying the number of vacancies in a critical cluster:

$$Z=\frac{1}{n^{\ast }}\sqrt{\frac{{\mathrm{\Delta }F}^{\ast }}{3\mathit{\pi kT}}}$$E27

While the nucleation rate is only directly proportional to most of the physical parameters in Eq. (25), the Arrhenius term dominates. Small changes in the height of the nucleation barrier lead to large variations in the equilibrium number of critical cavities available for nucleation and therefore the final nucleation rate. Smaller critical clusters are more likely to nucleate, such as in the case for clusters formed on grain boundaries as shown in Figure 8. The dihedral angle, δ, formed between the surface and the grain boundary as a result of the competing grain boundary surface energy, γ_gb, and the cluster surface energy, γ, reduce the volume and surface area [45] of a critical cluster, even though the curvature, r^*, is unchanged. Removing the previous grain boundary area (dashed gray line in Figure 8) also reduces the critical free energy:

Figure 8.

Cluster formed between two grains.

$$\delta =\mathit{\text{acos}}\left(\frac{{\gamma }_{\mathit{gb}}}{2\gamma }\right)$$E28

Nucleation is further boosted by the quicker diffusion of vacancies along the grain boundaries and the effect of converting the multiaxial stress state to an average stress on the grain boundary [46]. Also, real defects, such as dislocations, interacting with the grain boundary supply additional vacancies which can effectively increase the driving force by several gigapascals [47]. These effects predict cavity nucleation almost exclusively at grain boundaries and do not require extreme threshold stresses for nucleation, which is an enduring criticism of classical nucleation of cavities [48].

Finally, a theory based on generalized broken bonds (GNBBs) [49] is used to calculate the free surface energy from the energy of vacancy formation and a correction is applied [50] when dealing with nanosized critical clusters.

Diffusional cavity growth is less controversial and commonly assumed to follow the rate of radial growth in Eq. (29) by Hull and Rimmer [51]. Its resemblance to the Svoboda, Fischer, Fratzl, Kozeschnik (SFFK) model [20] used in precipitate growth simulations further strengthens its prestige:

$$\dot{r}=\frac{D{X}_v{a}^3}{\mathit{kT}r}\left(\sigma -\frac{2\gamma }{r}\right)$$E29

The sintering stress, 2γ/r, that opposes the driving stress, σ, also explains the shrinking of subcritical clusters predicted by CNT.

7. Pore formation model: model implementation

The equations for nucleation and growth are integrated into a Kampmann-Wagner framework [52] at a constant temperature and external stress state. At regularly spaced time intervals, a class of newly formed cavities with a population derived from Eq. (25) is formed, each with a radius slightly (20%) above the critical radius from Eq. (23). During the intervals, the respective radii of all classes grow according to Eq. (29). As the available nucleation sites are used, the nucleation rate diminishes. The simulation ends when all nucleation sites are consumed and there is no more uncavitated grain boundary area. The number of nucleation sites at grain boundaries is calculated [53] from average grain diameters, assuming all grains to be tetrakaidekahedral (Table 5).

Var.	Description	Unit	Source
γ	Free surface energy	J m^−2	[49]
σ	External stress	Pa	Experiment
Ns	Number of nucleation sites	—	Experiment and [53]
T	Temperature	K	Experiment
σD	Driving force of defects	Pa	[47]
γgb	Grain boundary energy	J m^−2	[54]
Qv	Vacancy formation energy	J mol^−1	[55]
D	Diffusion coefficient	m^2 s^−1	[56]
a	Interatomic spacing	m	[57]
dg	Grain diameter	m	Experiment

Table 5.

Model parameters and physical constants.

8. Pore formation model: results and comparison with experiments

Figures 9 and 10 compare simulated results of nucleated cavities with experimental results obtained from secondary electron microscopy and density measurements. Case studies comparing our model to experimental investigations have been published [42, 58].

Figure 9.

Comparison between simulated (red) and measured (blue) histograms of cavities in Nickel-based alloy 625 after 5500 h at 700°C and 183 MPa.

Figure 10.

Comparison between simulated (red) and measured (blue) histograms of cavities in P23 steel after 9000 h at 600°C and 90 MPa.

Var.	Description	Unit
a1	Parameter for glide velocity	[m/s]
ag	Size of unit cell	[m]
A	Material-specific constant	[–]
$\alpha$	Dislocation interaction factor	[–]
b	Burgers vector ½ <111> for bcc	[m]
$\beta$	Parameter for density of sources	[–]
cs	Weighting factor for statics in the internal stress	[–]
${\delta }_{\mathit{anh}}$	Dislocation annihilation length	[m]
Dcav	Cavitation damage parameter	[–]
Ds	Lattice diffusion coefficient	[m^2/s]
Dvp	Pipe diffusion coefficient	[m^2/s]
$\epsilon$	Creep strain	[–]
${\eta }_v$	Transfer coefficient for defects into jogs	[–]
G	Shear modulus	[Pa]
$\gamma$	Free surface energy	[J m^−2]
${\gamma }_{\mathit{GB}}$	Specific grain boundary energy	[J m^−2]
hb	Mean dislocation spacing within the subgrain wall	[m]
kB	Boltzmann constant	[J/K]
Kc	Holt constant	[–]
ki	Precipitate coarsening parameter of precipitate class i	[s^−1]
$L_{\alpha }$	Parameter for elastic interactions between dislocation/defects	[m]
$L_p$	Diffusion path of core-vacancy before evaporation into lattice	[m]
M	Taylor factor	[–]
Msgb	Subgrain boundary mobility	[m/Pa s]
Ni	Precipitate number density of class i	[m^−3]
Ns	Number of nucleation sites	–
$\nu$	Poisson’s ratio	[–]
Peff	Effective subgrain growth pressure	[Pa]
Q	Activation energy for dislocation glide	[J]
Qv	Vacancy formation energy	[J mol^−1]
ri	Precipitate radius of particle class i	[m]
Rsgb	Subgrain radius	[m]
${\rho }_b$	Boundary dislocation density	[m^−2]
${\rho }_m$	Mobile dislocation density	[m^−2]
${\rho }_s$	Static dislocation density	[m^−2]
${\rho }_t$	Total dislocation density	[m^−2]
${\sigma }_{\mathit{app},0}$	Applied initial stress	[Pa]
${\sigma }_{\mathit{app}}$	Corrected applied stress	[Pa]
${\sigma }_i$	Internal stress	[Pa]
σD	Driving force of defects	[Pa]
t	Time	[s]
tR	Rupture time	[s]
T	Temperature	[K]
${\theta }_i{f}_i$	Mean misorientation of subgrains	[rad]
vc	Climb velocity	[m/s]
vcl	Lattice diffusion share of vc	[m/s]
vcp	Pipe diffusion share of vc	[m/s]
veff	Effective velocity of mobile dislocations	[m/s]
vg	Glide velocity	[m/s]
Vr	Apparent activation volume for thermally activated glide	[m^3]
$\mathit{\Delta W}$	Activation energy difference between pipe and lattice diffusion	[J]
$\Omega$	Atomic volume (unit cell size divided by two atoms in bcc)	[m^3]
$\zeta$	Fraction of statics annihilating when transforming to boundaries	[–]

Table 6.

List of symbols and variables used within the creep model.

9. Conclusion

We have introduced a complex physically based creep model and demonstrated its capabilities in the case of the martensitic steel P91. The model is capable of simulating the creep deformation as well as the microstructural evolution during creep. As soon as some final parameters have been set, based on a single-creep experiment, those parameters can be used for simulating the creep behavior over a wide range of stresses allowing for extrapolating the creep behavior. Current results suggest an extrapolation of the creep lifetime by a factor of at least 6 over a reference experiment. Furthermore, we have introduced a physically based model for the formation of creep pores due to vacancy diffusion, which is also showing very good agreement with experiments.

Acknowledgments

The authors gratefully acknowledge funding from Austrian Science Fund (FWF) within project “Software development on dislocation creep in alloys” (P-31374).