\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[greek,romanian,ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{Modelling of Cavity Nucleation, Early-stage Growth and Sintering in
Polycrystal under Creep-fatigue Interaction}



\author[1]{J.-D. Hu}%
\author[1]{C.-J. Liu}%
\author[1]{Fuzhen Xuan}%
\author[2]{Bo Chen}%
\affil[1]{East China University of Science and Technology School of Mechanical Engineering and Power Engineering}%
\affil[2]{University of Leicester}%


\vspace{-1em}



  
  \date{June 17, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
A mechanistic based cavitation model that considers nucleation,
early-stage growth and sintering under creep-fatigue interaction is
proposed. The number density of cavities \selectlanguage{greek}\emph{ρ} \selectlanguage{english}and their evolution
during multi-cycle creep-fatigue loading are predicted. Both the cavity
nucleation and early-stage growth rates, controlled by grain boundary
(GB) sliding mechanism during the tension phase, are formulised as a
function of local normal stress \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}. The cavity
sintering that occurs during the compression phase is described as a
function of \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}, but the mechanism switches to
the unconstrained GB diffusion. By examining various load waveform
parameters, results provide important insights into experimental design
of studying the creep-dominated cavitation process under creep-fatigue
interaction. First, creep-fatigue test with initial compression will
promote higher \selectlanguage{greek}\emph{ρ} \selectlanguage{english}value compared to that with initial tension, if
the unbalanced stress hold time in favour of tension is satisfied.
Second, the \selectlanguage{greek}\emph{ρ} \selectlanguage{english}value does not have a monotonic dependence on
either the compressive hold time or stress level, because of their
competing effect on nucleation and sintering. Third, the optimum value
of stress variation rate exists in terms of obtaining the highest
\selectlanguage{greek}\emph{ρ} \selectlanguage{english}value due to sintering effect.%
\end{abstract}%



\sloppy


\textbf{Modelling of Cavity Nucleation, Early-stage Growth and Sintering
in Polycrystal under Creep-fatigue Interaction}

J.-D. Hu\textsuperscript{a,b,}11Present Address: Shanghai Nuclear
Engineering Research and Design Institute, Shanghai 200233,
China\textsuperscript{*}, C.-J. Liu\textsuperscript{a}, F.-Z.
Xuan\textsuperscript{a}, B. Chen\textsuperscript{b**}

\emph{\textsuperscript{a}Key Laboratory of Pressure Systems and Safety,
Ministry of Education, School of Mechanical and Power Engineering, East
China University of Science and Technology, Shanghai 200237, China}

\emph{\textsuperscript{b}School of Engineering, University of Leicester,
Leicester, LE1 7RH, UK}

*Corresponding authors. E-mail address:\emph{*jh788@leicester.ac.uk};

\emph{**bo.chen@leicester.ac.uk}

\textbf{Abstract}

A mechanistic based cavitation model that considers nucleation,
early-stage growth and sintering under creep-fatigue interaction is
proposed. The number density of cavities \selectlanguage{greek}\emph{ρ} \selectlanguage{english}and their evolution
during multi-cycle creep-fatigue loading are predicted. Both the cavity
nucleation and early-stage growth rates, controlled by grain boundary
(GB) sliding mechanism during the tension phase, are formulised as a
function of local normal stress \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}. The cavity
sintering that occurs during the compression phase is described as a
function of \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}, but the mechanism switches to
the unconstrained GB diffusion. By examining various load waveform
parameters, results provide important insights into experimental design
of studying the creep-dominated cavitation process under creep-fatigue
interaction. First, creep-fatigue test with initial compression will
promote higher \selectlanguage{greek}\emph{ρ} \selectlanguage{english}value compared to that with initial tension, if
the unbalanced stress hold time in favour of tension is satisfied.
Second, the \selectlanguage{greek}\emph{ρ} \selectlanguage{english}value does not have a monotonic dependence on
either the compressive hold time or stress level, because of their
competing effect on nucleation and sintering. Third, the optimum value
of stress variation rate exists in terms of obtaining the
highest\selectlanguage{greek}\emph{ρ} \selectlanguage{english}value due to sintering effect.

\textbf{Keywords:} Creep; Cavity; Viscoplastic; Analytic solutions;
Cavity sintering

\section*{1. Introduction}

{\label{introduction}}

Creep-dominated creep-fatigue interaction is a long-term material
failure mode that involves cyclic loading \textsuperscript{1}. Lack of
fundamental understanding on this topic prevents us to commercialise the
Generation IV high-temperature nuclear power plants with a design life
of 60 years \textsuperscript{2}. Although the mechanistic-based
descriptions of cavity nucleation, growth and coalescence under creep
have been established \textsuperscript{3}, little attempt has been made
to reveal the mechanism of cavity nucleation and its early-stage radius
change under creep-fatigue interaction \textsuperscript{4, 5}.

The empirical relationship between cavity nucleation rate under creep
and Monkman-Grant constant has been established by
Davanas\textsuperscript{6}, providing a good agreement with experimental
data. But it is unclear about its suitability for creep-fatigue
interaction. Cavitation under fatigue is known to be different from that
under creep, in a sense that fatigue induced cavities are smaller in
size but higher in their number density \textsuperscript{5}. There are
creep-fatigue lifetime prediction models that consider the creep
cavitation mechanism, e.g. Nam \textsuperscript{7}. However, the model
cannot explain the positive relationship between the cavity nucleation
and tensile hold time \textsuperscript{8}. Recent work by Wen,
Srivastava\textsuperscript{9} and Barbera, Chen \textsuperscript{10}
used creep damage models to predict crack growth rate under
creep-fatigue, by incorporating the late-stage cavity growth. However,
they did not consider the physical process of cavity nucleation or
early-stage growth.

Moreover, little knowledge has been gained regarding the cavity
sintering. Compressive loading can cause the nucleated creep cavities
shrink in radius and forces them to be removed completely under some
circumstances \textsuperscript{11, 12}. This means that their number
density after certain number of cycles is not the simple sum of
nucleated cavities from each cycle. To the best of the authors'
knowledge, creep cavitation model considering both the nucleation and
sintering events under creep-fatigue interaction does not exist so far.

The mechanism of cavity nucleation was initially proposed by
Greenwood\textsuperscript{13}, and developed further by Raj and
Ashby\textsuperscript{14}. It has been accepted that vacancies
agglomerate and form stable nuclei assisted by local normal stress. The
local normal stress can be significantly higher than the far-field
stress\textsuperscript{15}. One of the causes is grain boundary (GB)
sliding induced stress concentration. Extensive cavitation was found in
copper bi-crystals which had been subjected to pre-strain in favour of
GB sliding followed by creep loading \textsuperscript{16, 17}. Min and
Raj\textsuperscript{18} proposed a model to predict the local normal
stress under creep-fatigue loading based on the GB sliding mechanism. By
considering the creep effect on local normal stress, our model is
capable of predicting the cavity nucleation under one-cycle
creep-fatigue \textsuperscript{19}.

Modelling the early-stage radius change is key to understanding the
sintering process under creep-fatigue. Note that sintering shares a high
degree of commonality with the cavity growth. The classical growth
models were established under either or both of vacancy diffusion and
matrix deformation (Cocks \textsuperscript{20,} Cocks and
Ashby\textsuperscript{21,} Needleman and Rice \textsuperscript{22,}
Chuang, Kagawa \textsuperscript{23,} Chen and Argon
\textsuperscript{24}), and the growth rate is controlled by the
far-field stress. Nevertheless, these models might not be suitable for
the early-stage growth. First, the radius of nucleated cavities
(\textasciitilde{}5 nm) is much smaller than those commonly defined
`small' cavities (\textasciitilde{}1 \selectlanguage{greek}μ\selectlanguage{english}m)\textsuperscript{3}. This
indicates that the vacancy flow near the cavities is highly sensitive to
the local normal stress\textsuperscript{25}. Second, the cavity growth
mechanism map\textsuperscript{26} does not consider the role of GB
sliding on the early-stage growth, despite its importance
\textsuperscript{25}.

In this paper, the number density of cavities and their time evolution
during multi-cycle creep-fatigue loading are modelled. The
stress-controlled load waveform is considered because a higher creep
damage would be generated in comparison with the strain-controlled
one\textsuperscript{27, 28}. Type 316 stainless steel is selected for
twofold reasons: material parameters are available \textsuperscript{18,
29-31} and its wide application to the power generation
industry\textsuperscript{2, 29}.

\section*{2. Modelling approach}

{\label{modelling-approach}}

\subsection*{2.1 Theoretical modelling}

{\label{theoretical-modelling}}

Modelling cavitation under creep-fatigue contains two parts. The first
part focuses on describing the relationship between the applied stress
and local normal stress, and the cavity nucleation governed by vacancy
accumulation. The second part concerns about the radius change of
nucleated cavities during creep-fatigue loading. These two parts are
integrated through a numerical framework to give the final prediction to
the number density of cavities. The governing equations are described
below.

The rate of cavity nucleation () on the particles at GB can be described
as \textsuperscript{14}:

\par\null

where \selectlanguage{greek}Ω \selectlanguage{english}is the atomic volume and k is the Boltzmann constant. \selectlanguage{greek}δ\selectlanguage{english}\emph{D}
\textsubscript{B} is the GB thickness multiplied by its self-diffusion
coefficient. \selectlanguage{greek}γ \selectlanguage{english}is the free surface energy.\emph{F} \textsubscript{v} is
a shape factor related to the cavity volume, and it has a value of
0.1585 for cavity nucleation at the top of GB particles
\textsuperscript{14}. \selectlanguage{greek}\emph{ρ} \selectlanguage{english}\textsubscript{max} is the maximum number
density of potential nucleation sites, and its value can be worked out
through the relation of \textsuperscript{32}, by assuming that grain has
an idealised hexagonal shape. \emph{f} \textsubscript{b}is the area
fraction of GB particles, is their average radius (in \selectlanguage{greek}μ\selectlanguage{english}m), and \emph{d}
is the grain size (in \selectlanguage{greek}μ\selectlanguage{english}m). A coefficient of 10\textsuperscript{6}
appears due to the unit conversion when expressing the number density of
cavities in mm\textsuperscript{-2}.

The cavity nucleation rate in Eq. depends on three variables: number
density of nucleated cavities, local normal stress and temperature,
symbolised as \selectlanguage{greek}\emph{ρ} , \selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} and \emph{T} ,
respectively. Accordingly, the relationship between the local normal
stress and applied shear stress is formulised as follow:

\par\null

where , and are the applied shear stress, its change rate and the local
normal stress at time \emph{t} . The derivation of Eq. is described in
Appendix A. is the local normal stress at time \emph{t} +\selectlanguage{greek}Δ\selectlanguage{english}\emph{t} .
Note \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} is the time interval chosen for the numerical computation
that will be described in Section 2.3. \selectlanguage{greek}\emph{τ} \selectlanguage{english}is time constant defined
as:

\par\null

where \selectlanguage{romanian}\emph{ƞ} \textsubscript{b} is the damping coefficient of GB
sliding. \emph{k} \textsubscript{e} is the elastic modulus, defined
as\emph{G} /{[}0.57\emph{d} (1-\emph{v} ){]}. is the transient creep
damping coefficient that has a form of\emph{ƞ}
\textsubscript{p}/{[}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s}-\emph{f}
\textsubscript{b}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}{]}\textsuperscript{\emph{n}
-1}.\emph{G} and \emph{v} are shear modulus and Poisson's ratio.\selectlanguage{romanian}\emph{ƞ}
\textsubscript{p} is the reciprocal of pre-exponential factor and
\emph{n} is the exponent in the Norton power-law creep equation. More
details about the cavity nucleation model can be found in Hu,
Xuan\textsuperscript{19}.

To obtain the actual size of the nucleated cavities, the radius change
rate needs to be determined during the subsequent loadings. For cavity
growth, Chen \textsuperscript{25} pointed out that GB sliding played an
important role in the early-stage growth. The transient GB sliding
activity could immediately wedge open the supercritical nuclei along the
particles \textsuperscript{33}. According to this `crack sharpening'
mechanism, the cavity tip velocity is limited by surface
diffusion\textsuperscript{3}.

Fig. 1a shows a schematic diagram of the cavity nucleation and
early-stage growth assisted by GB sliding with a rate of . Because of
insufficient surface diffusion, the nucleated cavity tends to form an
irregular crack-like shape. This implies that the growth rate is
independent of the cavity radius, given that the process is limited by
the vacancy diffusion rate near the cavity tip. Based on the work of
Chuang, Kagawa \textsuperscript{23}, equations were formulised by Nix,
Yu \textsuperscript{34} to describe cavity growth rate controlled by the
vacancy surface diffusion under low and high stresses:

\par\null

where \emph{D} \textsubscript{S} is the surface diffusion
coefficient.\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S} is the width of surface
diffusion, which is related to atom volume \selectlanguage{greek}Ω, \selectlanguage{english}through\selectlanguage{greek}\emph{δ}
\selectlanguage{english}\textsubscript{S}=\selectlanguage{greek}Ω\selectlanguage{english}\textsuperscript{1/3}. \emph{b} is the limiting
cavity radius, equalling the half of average cavity spacing, \selectlanguage{greek}λ/\selectlanguage{english}2
\textsuperscript{23}. \selectlanguage{greek}λ \selectlanguage{english}is the average distance between two nucleation
sites, and its value can be estimated as below for an idealised
hexagonal grain shape:\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image17/image17}

\end{center}
\end{figure}

Fig. 1: Schematic diagrams of a cavity that (a) grows by the surface
diffusion limited GB sliding; (b) shrinks by the unconstrained GB
diffusion.

In the growth model \textsuperscript{34}, the cavity has been assumed to
have a crack-like shape, Fig. 1a. The growth process is driven by the
accumulation of vacancies at the cavity tip, leading to the phenomenon
similar to a crack-tip continuous extension. Meanwhile, the surface
diffusion is not very much faster than GB diffusion under low stress,
while the surface diffusivity becomes much higher than that of GB under
high stress. Refer to Appendix B for the derivation of stress level that
separates the two conditions.

The reason for ignoring the creep effect on the early-stage growth can
be justified as follows. The growth mechanism map proposed by Miller,
Hamm \textsuperscript{26} informs that the cavity radius change is
controlled by unconstrained GB diffusion when the ratio of cavity
spacing (\selectlanguage{greek}λ\selectlanguage{english}) to its diameter (2\emph{r} ) is greater than 10. In the case
of early-stage cavitation, the value of \selectlanguage{greek}λ/\selectlanguage{english}2\emph{r} falls into the range
of 200 to 3000, given that the value range of \emph{r} is 3 to 5
nm\textsuperscript{3} and that of \selectlanguage{greek}λ \selectlanguage{english}is 1 to 10 \selectlanguage{greek}μ\selectlanguage{english}m \textsuperscript{23}.
Therefore, the unconstrained GB diffusion is the rate-limiting factor
for the early-stage growth.

Because of the analogy of shrinkage to growth, shrinkage occurs by the
GB vacancy diffusion that pushes vacancies out of the cavity under the
compressive \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}. The reversed direction of GB
sliding, coupled with the reduced sharpness of the cavity tip, causes a
more uniform distribution of vacancies within the cavity, as shown in
Fig. 1b. Since the vacancies are no longer concentrated at the tip, more
surface areas become the diffusional interface. Ultimately, the cavity
shape changes back to sphere, that is the equilibrium shape when the
surface diffusion rate is fast enough \textsuperscript{23}.

To calculate the shrinkage rate, equation proposed by
Riedel\textsuperscript{35} for cavity growth under the unconstrained GB
diffusion is adapted:

\par\null

When the cavity radius \emph{r} becomes larger than \selectlanguage{greek}λ/\selectlanguage{english}4.24, is set to
zero to prevent the unrealistic positive value under compression. This
upper limit value has been used to define the largest cavity size. The
assumption behind the Riedel model is that cavity is not perfectly
spherical, which aligns very well with the considered crack-like
cavities in the present work. If the cavity radius \emph{r} becomes less
than the critical nuclei size defined by the value of 2\selectlanguage{greek}γ/\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} \textsuperscript{14}, cavities are considered as
sintered \textsuperscript{12, 36}.

\subsection*{2.2 Material parameters and temperature
effect}

{\label{material-parameters-and-temperature-effect}}

Except for the surface diffusion factor \emph{D} \textsubscript{S}, all
the material parameters have been sourced from the authoritative
literature listed in Table 1. \emph{D} \textsubscript{S} values in the
temperature range of 500-600 @\selectlanguage{english}C are not available; hence the diffusivity
ratio of\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S}\emph{D} \textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english}
\textsubscript{B}=0.001 has been used for calculations. This meets the
criterion of\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S}\emph{D}
\textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english} \textsubscript{B}\textless{}\textless{}1 for
a crack-like cavity \textsuperscript{34}. Justification of the\emph{D}
\textsubscript{S} value selection through a sensitivity study is
described in Appendix B. In this context, our model has been established
based on the state-of-the-art mechanistic understanding of cavitation
and the prediction results are reliable (no approximation or
extrapolation).

Creep database for Type 316H stainless steel \textsuperscript{29} has
been used to derive the power-law creep related parameters (1/\selectlanguage{romanian}\emph{ƞ}
\textsubscript{p} and n, Table 1) for two reasons. First, this material
was subjected to 65,015 h in service at temperatures between 490 and 530
\selectlanguage{english}@C, prior to creep testing \textsuperscript{37}. Thus, it is unlikely
that GB particle evolution occurs during the test (hence a fixed
\emph{f} \textsubscript{b} value). Second, the operating pressure of
less than 20 MPa would be too low to generate any noticeable creep
cavities in this ex-service material, and hence no need to consider
pre-existing cavities in our calculations. Note that the maximum
allowable carbon content for Type 316 stainless steel is \textless{}0.08
wt.\% according to \textsuperscript{38}. This means that the chosen
material with a carbon content of 0.06 wt.\%\textsuperscript{39} falls
into the category of Type 316 stainless steel.

Table 1. Material parameters in the proposed cavity nucleation and
sintering model for Type 316 stainless steel\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
Parameter & Symbol & Value & Unit & Temperature range/ @\selectlanguage{english}C &
References\tabularnewline
\midrule
\endhead
Average radius of the particle at GB & & 1.0 & \selectlanguage{greek}μ\selectlanguage{english}m & - & Min and Raj
\textsuperscript{18}\tabularnewline
Average size of grains & \emph{d} & 40 & \selectlanguage{greek}μ\selectlanguage{english}m & - & Min and Raj
\textsuperscript{18}\tabularnewline
Shear modulus & \emph{G} & 7.6\selectlanguage{ngerman}×10\textsuperscript{10} & Pa & - & Min and
Raj \textsuperscript{18}\tabularnewline
Poisson's ratio & \emph{v} & 0.31 & - & - & Min and Raj
\textsuperscript{18}\tabularnewline
Atom volume & \selectlanguage{greek}Ω & 1.21\selectlanguage{ngerman}×10\textsuperscript{-29} & m\textsuperscript{3} &
- & Min and Raj \textsuperscript{18}\tabularnewline
Area fraction of the particles & \emph{f}\textsubscript{b} & 0.2 & - & -
& Min and Raj \textsuperscript{18}\tabularnewline
Shape factor related to the cavity volume & \emph{F}\textsubscript{v} &
0.1585 & - & - & Raj and Ashby \textsuperscript{14}\tabularnewline
Free surface formation energy & \selectlanguage{greek}γ & 0.80 & J.m\textsuperscript{-2} & - &
Mortimer and Nicholas \textsuperscript{31}\tabularnewline
Boltzmann constant & k & 1.3806\selectlanguage{ngerman}×10\textsuperscript{-23} &
J.K\textsuperscript{-1} & - & Universal constant\tabularnewline
Gas constant & R & 8.314 &
J.mole\textsuperscript{-1}.K\textsuperscript{-1} & - & Universal
constant\tabularnewline
GB diffusion activation energy & Q\textsubscript{c} &
167\selectlanguage{ngerman}×10\textsuperscript{3} & J.mole\textsuperscript{-1} & - & Frost and
Ashby \textsuperscript{30}\tabularnewline
GB thickness multiplied by its self-diffusion coefficient &
\selectlanguage{greek}\emph{δD}\selectlanguage{english}\textsubscript{B} &
2.0\selectlanguage{ngerman}×10\textsuperscript{-13}exp(-Q\textsubscript{c}/R\emph{T}) &
m\textsuperscript{3}.s\textsuperscript{-1} & - & Frost and Ashby
\textsuperscript{30}\tabularnewline
Free surface diffusion coefficient\textsuperscript{*} &
\emph{D}\textsubscript{S} &
9.0\selectlanguage{ngerman}×10\textsuperscript{-7}exp(-167×10\textsuperscript{3}/R\emph{T}) &
m\textsuperscript{2}.s\textsuperscript{-1} & 500-600 & Not available in
literature\tabularnewline
Creep exponent & n & 11.83 & - & 500-600 & Chen, Flewitt
\textsuperscript{29}\tabularnewline
Pre-exponential factor & 1/\selectlanguage{romanian}\emph{ƞ}\textsubscript{p} &
2.199\selectlanguage{ngerman}×10\textsuperscript{-32} &
h\textsuperscript{-1}.MPa\textsuperscript{-n} & 500 & Chen, Flewitt
\textsuperscript{29}\tabularnewline
& & 8.754\selectlanguage{ngerman}×10\textsuperscript{-32} & & 525 &\tabularnewline
& & 6.489\selectlanguage{ngerman}×10\textsuperscript{-31} & & 550 &\tabularnewline
& & 3.035\selectlanguage{ngerman}×10\textsuperscript{-30} & & 575 &\tabularnewline
& & 2.005\selectlanguage{ngerman}×10\textsuperscript{-29} & & 600 &\tabularnewline
\bottomrule
\end{longtable}

*Note: Value has been calculated through the relation of\selectlanguage{greek}\emph{δ}
\selectlanguage{english}\textsubscript{S}\emph{D} \textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english}
\textsubscript{B}=0.001.

Fig. 2a shows the relationship between \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} and
cavity growth rate at three different temperatures of 500/550/600 @\selectlanguage{english}C, as
predicted by Eq.(4). The cavity growth rate increases by one order of
magnitude with temperature increasing from 500 to 600 @\selectlanguage{english}C. However, the
transition value of \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} (\textasciitilde{}400
MPa), which differentiates the cavity growth mechanism under low and
high stresses, is not affected by the temperature. This is as expected
because the GB diffusion has the same activation energy as the surface
diffusion (Table 1). Effects of the radius \emph{r} and\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} on cavity shrinkage rate are shown in a contour plot
of Fig. 2b, as predicted by Eq.(6). Different colours signify
temperatures of 500/550/600 @\selectlanguage{english}C, respectively. For each temperature, the
contour curves with the largest value are always towards the top left,
indicating that a higher compressive stress\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} or
smaller cavity radius \emph{r} results in a higher shrinkage rate. By
comparing different coloured curves with the same value, the shrinkage
rate is found to be positively correlated with temperatures.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image21/image21}

\end{center}
\end{figure}

Fig. 2: Effects of: (a) local normal stress \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}
on cavity growth rate; (b) cavity radius \emph{r} and\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} on shrinkage rate at 500 @\selectlanguage{english}C (blue), 550 @\selectlanguage{english}C (black) and
600 @\selectlanguage{english}C (red).

\subsection*{2.3 Numerical
implementation}

{\label{numerical-implementation}}

Numerical computation of the cavity number density during creep-fatigue
was performed in a sectionalised manner. Let's consider a time
section\emph{i} , the time period starts from (\emph{i} -1)\selectlanguage{greek}Δ\selectlanguage{english}\emph{t} and
ends at\emph{i} \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} , and a group of cavities can be nucleated
during the time interval \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} . This group of cavities are named as
`cavities of the \emph{i} -th group'. All the cavities belonging to this
group would have the same radius and hence they are assigned as the
nucleated cavities within the time section \emph{i} . The number density
of cavities of the \emph{i} -th group is marked
as\selectlanguage{greek}\emph{ρ\textsubscript{i}}\selectlanguage{english} , with its value equalling the product of
nucleation rate and \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} . The cavity is considered as nucleated
when its radius reaches the value of 2\selectlanguage{greek}\emph{γ} /\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n}, and we assign the
symbol\textsuperscript{0}\emph{r\textsubscript{i}} representing the
critical nucleation radius.

Since the growth rate is related to their radius, the cavities
of\emph{i} -th group would have their characteristic rate at later time
sections. For a time section \emph{k} (\emph{k} \textgreater{}\emph{i}
), the cavity radius of the \emph{i} -th group, symbolised
as\textsuperscript{\emph{k} -\emph{i}} \emph{r\textsubscript{i}} , can
be worked out:

where is the growth rate at time section \emph{j}(\emph{j}
\textgreater{}\emph{i} ). This implies that cavities of the\emph{i} -th
group are considered as sintered, when the value
of\textsuperscript{\emph{k} -\emph{i}} \emph{r\textsubscript{i}} is less
than that of \textsuperscript{0}\emph{r\textsubscript{i}} . Thus, the
sintered condition is defined as:

\par\null

A flag variable (e.g. Flag(\emph{i} ) for cavities of the \emph{i} -th
group) is assigned to each cavity group as soon as cavities are
nucleated, and this will be updated at following time sections. The flag
value of 1 means that this group of cavities still exists, whilst 0
means that it is as sintered.

The number density of cavities \selectlanguage{greek}\emph{ρ} \selectlanguage{english}(\emph{k} ) at the end of time
section \emph{k} (\emph{k} \textgreater{}\emph{i} ) after taking
sintering into account can be formulised as:

\par\null

where \selectlanguage{greek}\emph{ρ\textsubscript{k}}\selectlanguage{english} is the number density of cavities
nucleated within the time section \emph{k} . The cavitated GB area
fraction \emph{f} (\emph{k} ) at the end of time section \emph{k} can be
then derived with the assumption of idealised hexagonal grain shape:

\par\null

The coefficient of 10\textsuperscript{6} has been generated due to the
unit conversion that involves the parameters of \selectlanguage{greek}\emph{ρ} \selectlanguage{english}in
mm\textsuperscript{-2}, \emph{r} in \selectlanguage{greek}μ\selectlanguage{english}m, and \emph{d} in \selectlanguage{greek}μ\selectlanguage{english}m.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image28/image28}

\end{center}
\end{figure}

Fig. 3: Overall flow chart of the numerical computation

The numerical algorithm is depicted in Fig. 3 using a flow chart. First,
the material parameters in Table 1 are read by our self-written MATLAB
script. Afterwards, \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} is read from input, \selectlanguage{greek}\emph{ρ} \selectlanguage{english}and\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} are set to zero, and the iteration begins at the time
section \emph{k} =1. Within each iteration, \emph{T} and\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{s} at time \emph{k} \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} are read first,
then\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} can be calculated by Eq.. Next, the
nucleation rate is computed by Eq., from which the initial \selectlanguage{greek}\emph{ρ} \selectlanguage{english}that
does not consider the sintering effect can be worked out. After this
step, the script enters into the sintering module, which is ruled by
Eq.. Note that the flag variable will be updated with a value of 0 if
the considered cavity group is found to be as sintered using the
criterion given in Eq.. Afterwards, the cavitated GB area fraction
\emph{f} can be calculated by Eq.. At the end of each iteration, time
section \emph{k}is updated by \emph{k} +1, and a new iteration begins.
Refer to Appendix C for the limited numerical error generated by using
the forward Euler method together with \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} =0.5 s.

\section*{3. Model prediction}

{\label{model-prediction}}

\subsection*{3.1 Key outputs}

{\label{key-outputs}}

A typical load waveform in stress-controlled creep-fatigue with initial
compression under \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s}[?]{[}-150, 150{]} MPa
and\emph{T} =550 @\selectlanguage{english}C is shown in Fig. 4. The first cycle is presented for
illustration purpose, and the following cycles are the same. It consists
of a pre-compressive stress hold for a period \emph{t}
\textsubscript{c}and a tensile stress hold \emph{t} \textsubscript{t}.
The magnitude of stress applied during compression is designated as
\selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c}, and the stress applied during tension is
\selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{t}. The stress variation rate during load
reversal is constant, symbolised as . Table 2 summarises the default
value for load waveform parameters and test temperature. They are used
to elaborate the key model outputs. A total number of 50 cycles have
been calculated.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image31/image31}

\end{center}
\end{figure}

Fig. 4: A typical stress-controlled creep-fatigue load waveform

Table 2. Load waveform parameters and temperature in their default
values\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Parameters & Default value & Units\tabularnewline
\midrule
\endhead
\emph{t}\textsubscript{c} & 20 & s\tabularnewline
\emph{t}\textsubscript{t} & 100 & s\tabularnewline
\selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}\selectlanguage{english}\textsubscript{c} & 150 & MPa\tabularnewline
\selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}\selectlanguage{english}\textsubscript{t} & 150 & MPa\tabularnewline
& 50 & MPa.s\textsuperscript{-1}\tabularnewline
\emph{T} & 550 & @\selectlanguage{english}C\tabularnewline
\bottomrule
\end{longtable}

Fig. 5a shows the model predicted \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} as a
function of the elapsed time/cycle. The maximum\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} in each cycle is the sum of the relaxed\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} in pre-compression and the increment of\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} induced by the load reversal. The maximum\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} starts from a value of 761 MPa in the 1st cycle, and
then decreases cycle by cycle. After 40 cycles, the maximum\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} reaches its steady state with a value of 310 MPa, as
highlighted by the hatched region in Fig. 5a. The evolution of\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} during the 1st cycle is shown in Fig. 5b, where the
magnitude change due to the load reversal from compression to tension
appears to be the same as that from tension to compression. In addition,
the inset shows that the magnitude of relaxed\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}
in tensile hold (54 MPa) is larger than that in compressive hold (11
MPa). This stress difference is caused by the longer tensile hold time
\emph{t} \textsubscript{t} than the compressive time \emph{t}
\textsubscript{c}. Thus, the net \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}change during
the non-saturation cycle is negative (e.g. -42 MPa in the 1st cycle),
leading to gradually decreased maximum value of\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} (Fig. 5a).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image32/image32}

\end{center}
\end{figure}

Fig. 5: Predicted time evolution curves of (a) the local normal
stress\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}; (b) \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} in the
1st cycle; (c) \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} in the 40th cycle (d) the
number density of cavities \selectlanguage{greek}\emph{ρ} \selectlanguage{english}and cavitated GB fraction \emph{f} .

When cyclic loading continues, the magnitude of relaxed\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} in tensile hold decreases, whereas that in compressive
hold increases. This is as expected because Eq.(2) includes an
exponential decay function, indicating that an increased initial stress
level would cause a faster stress relaxation over a fixed period of
time. When the steady state is reached, the magnitude of relaxed\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} in the tensile and compressive holds becomes
identical. This can be clearly seen in Fig. 5c for the 40th cycle, where
the tensile and compressive holds introduce 23 MPa stress relaxation,
respectively. This means that the net stress change is 0 MPa in the
steady state, being consistent with the observation in Fig. 5a.

Since nucleation rate is highly sensitive to \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n},
the maximum \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} is key to cavity nucleation (i.e.
the number density of cavities \selectlanguage{greek}\emph{ρ} \selectlanguage{english}). The model predicted
\selectlanguage{greek}\emph{ρ}\selectlanguage{english}evolution as well as cavitated GB fraction \emph{f} during 50
cycles of creep-fatigue loading are shown in Fig. 5d. The \selectlanguage{greek}\emph{ρ} \selectlanguage{english}curve
in black shows that most of the cavities are nucleated in the tensile
hold of the 1st cycle, while the following compressive hold contributes
to the major part of cavity sintering. After the first two cycles, the
change in\selectlanguage{greek}\emph{ρ} \selectlanguage{english}becomes negligible. By comparison, the \emph{f} curve
in red with its value read from the right axis (Fig. 5d) shows an
overall increasing trend. This implies that the increased \emph{f} value
in later creep-fatigue cycles is most likely related to the cavity
growth rather than the nucleation of new cavities.

To show this more evidently, two enlarged views of the \emph{f} curve
exhibiting a zigzag shape are highlighted in Fig. 5d. The zigzag shape
can be explained by the fact that cavities nucleate and grow up
when\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}\textgreater{}0, and they shrink
when\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}[?]0. The zigzag characteristic is more
noticeable in the 2nd cycle than the 43th cycle, indicating that the
cavity radius change rate reduces with the cycling. This implies that
the cavity sintering is less likely to occur in the later cycles of
creep-fatigue loading under the present waveform condition (i.e.\emph{t}
\textsubscript{t}\textgreater{}\emph{t} \textsubscript{c} in Fig. 4).
Since the model predicts that the later creep-fatigue cycles neither
nucleate many cavities (black curve in Fig. 5d) nor cause the already
nucleated cavities to be sintered (red curve in Fig. 5d), it is
appropriate to reduce fatigue cycles from 50 to 10 to save the
computational cost.

In the work by Min and Raj \textsuperscript{18}, one purpose-designed
creep-fatigue cycle at 625 @\selectlanguage{english}C was employed to create the increased
fatigue crack growth rate from 0.05 to 0.17 mm/cycle on the Type 316
stainless steel. This specific load-waveform included a 15 mins stress
hold under compression (-382 MPa), followed by a rapid load reversal
from compression to tension (382 MPa) with a strain rate of
2\selectlanguage{ngerman}×10\textsuperscript{-4} s\textsuperscript{-1}, and then applying a 1 hr
stress hold under tension. The cavitation damage at grain boundaries was
responsible for the increased crack growth rate. They further showed
that imposing more than one fatigue cycle could not generate a further
increase in the crack growth rate. Therefore, our model prediction in
terms of the importance of the 1st fatigue cycle agrees with the
experimental observation.

According to the predicted \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} with a magnitude
of 761 MPa (for a tensile hold stress \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{t} of 160
MPa in Fig. 5a) in the present work, the cavity nucleation rate would be
2.1\selectlanguage{ngerman}×10\textsuperscript{-10}mm\textsuperscript{-2}.s\textsuperscript{-1}.
When the magnitude of\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} increases to 1016 MPa,
by applying a higher tensile stress hold of 200 MPa, the nucleation rate
would change to 4.0\selectlanguage{ngerman}×10\textsuperscript{-1}
mm\textsuperscript{-2}.s\textsuperscript{-1}. Recall the work by Raj
(1978), the predicted cavity nucleation rate was found to vary from
10\textsuperscript{14} to
10\textsuperscript{-16}mm\textsuperscript{-2}.s\textsuperscript{-1}, if
the local normal stress\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} is decreased by an
order of magnitude from 10\textsuperscript{3} MPa to
10\textsuperscript{2} MPa. Thus, our model prediction regarding the
nucleation rate is as expected. It is well recognised (e.g. Evans
\textsuperscript{15}) that cavity nucleation under the GB sliding
mechanism would require a high stress concentration (in the order of
10\textsuperscript{3} to 10\textsuperscript{4} MPa) particularly
associated with the classical nucleation theory. Also, the nucleation
rate (and cavity density) is strongly dependent on the material
parameters \textsuperscript{19, 40}.

In this context, it is the predicted trend, rather than the specific
value of cavity number density, that provides an important but missing
guideline in terms of optimising the load-waveform design to produce
more creep cavitation damage with reduced creep-fatigue experimental
cost.

\subsection*{3.2 Effect of loading
sequence}

{\label{effect-of-loading-sequence}}

Creep-fatigue test can start either with the initial compression or
tension. Fig. 6a compares these two loading sequences in terms of
the\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} evolution. The load waveform parameters
and temperature used here are the same as that listed in Table 2,
i.e.\emph{t} \textsubscript{t}\textgreater{}\emph{t} \textsubscript{c}.
It can be seen that the test with initial compression has a higher
maximum\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} than that with initial tension. This
is because the pre-compression period can relax the\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n}, contributing to the higher maximum\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} in the tension phase. The difference of maximum
\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} between the two loading conditions is 12 MPa,
which is in line with the relatively short\emph{t} \textsubscript{c} of
20 s. The number density of cavities\selectlanguage{greek}\emph{ρ} \selectlanguage{english}after 10 cycles was
calculated as 2.6\selectlanguage{ngerman}×10\textsuperscript{-9}mm\textsuperscript{-2} for the
test with initial compression, which is approximately one order of
magnitude higher than that with initial tension (\selectlanguage{greek}\emph{ρ}
\selectlanguage{english}=4.9\selectlanguage{ngerman}×10\textsuperscript{-10} mm\textsuperscript{-2}), Fig. 6b.

Now let's consider a completely different scenario in which the
compressive hold is prolonged, i.e.\emph{t}
\textsubscript{c}\textgreater{}\emph{t} \textsubscript{t}. To this end,
the \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} evolution for a creep-fatigue loading
with \emph{t} \textsubscript{t}=20 s and\emph{t} \textsubscript{c}=100 s
was calculated. The test with initial compression is compared with that
with initial tension. For both cases, the \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}
increases cycle by cycle, and the maximum\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}
appears at the last cycle, Fig. 6c. The test with initial compression
has a higher maximum \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}than that with initial
tension. At the 1st cycle, the stress difference is 54 MPa, reducing to
23 MPa at the 10th cycle (the last cycle of calculation). This implies
that the maximum \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}difference between the two
loading conditions will eventually diminish to 0 MPa if the cyclic
loading continues.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image33/image33}

\end{center}
\end{figure}

Fig. 6: Creep-fatigue test with initial compression is compared with
that with initial tension: (a) and (b) the time evolution of\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} and \selectlanguage{greek}\emph{ρ} \selectlanguage{english}for the test condition of\emph{t}
\textsubscript{t}=100 s, \emph{t} \textsubscript{c}=20 s (i.e.\emph{t}
\textsubscript{t}\textgreater{}\emph{t} \textsubscript{c}); (c) and (d)
the time evolution of \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} and \selectlanguage{greek}\emph{ρ} \selectlanguage{english}for the
test condition of \emph{t} \textsubscript{t}=20 s,\emph{t}
\textsubscript{c}=100 s (i.e.\emph{t}
\textsubscript{t}\textless{}\emph{t} \textsubscript{c})

Fig. 6d depicts the time evolution of \selectlanguage{greek}\emph{ρ} \selectlanguage{english}curves for these two test
conditions. After each cycle, the \selectlanguage{greek}\emph{ρ} \selectlanguage{english}curves drop to 0, indicating
that all of the cavities nucleated in the tension phase are sintered
during the subsequent long compressive hold (\emph{t}
\textsubscript{c}=100 s). By comparison of Fig. 6b with 6d, it can be
known that the nucleation rate in the tests with\emph{t}
\textsubscript{c}=100 s is almost ten orders of magnitude higher than
that in the tests with \emph{t} \textsubscript{c}=20 s. This is due to
the fact that the maximum \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} is higher in the
test with \emph{t} \textsubscript{c}=100 s (1043 MPa in Fig.6a) when
compared with that test with \emph{t} \textsubscript{c}=20 s (761 MPa in
Fig.6c). However, this does not necessarily lead to the overall increase
of cavitation damage as the prolonged compressive hold would also cause
sintering and close all cavities nucleated in the previous tension
phase.

\subsection*{3.3 Combined effect of tensile and compressive
holds}

{\label{combined-effect-of-tensile-and-compressive-holds}}

Combined effect of \emph{t} \textsubscript{t} and\emph{t}
\textsubscript{c} on the final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}after 10 cycles is presented in
Fig. 7a by the contour plot. Note that except for the\emph{t}
\textsubscript{t} and \emph{t} \textsubscript{c}, the default value as
given in Table 2 is used for all the other load-waveform parameters and
temperature. The contour curves of 0.1 and 0.0001 in Fig. 7a are very
close to each other, suggesting that the contour curve of 0.1 can be
used as the lower bound value. In other words, below this value, the
final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}is negligible. There is a threshold value of\emph{t}
\textsubscript{t} for cavitation when\emph{t} \textsubscript{c}=0 s,
indicating that the nucleated cavities can be sintered by the load
reversal from tension to compression. In addition, under the condition
of \emph{t} \textsubscript{c}=0 s, the contour curves with values of 0.5
and beyond do not intercept with the x-axis. This indicates that there
is an upper bound value for\emph{t} \textsubscript{t}; any prolonged
tensile hold\emph{t} \textsubscript{t} than this critical value would
not further contributing to the final \selectlanguage{greek}\emph{ρ} . \selectlanguage{english}This is probably
attributed to the rapid decrease of \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} during
tensile hold as informed by the exponential decay function in Eq.(2).

Now let's focus on the changing values of \emph{t} \textsubscript{c}.
When \emph{t} \textsubscript{t} is less than 60 s, increasing\emph{t}
\textsubscript{c} helps to reduce the final \selectlanguage{greek}\emph{ρ} . \selectlanguage{english}But when\emph{t}
\textsubscript{t} is greater than 60 s, an increased\emph{t}
\textsubscript{c} promotes the increased final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}until reaching
the allowable value as highlighted by the magenta dash curve in Fig. 7a.
It is also evident that the allowable \emph{t}
\textsubscript{c}increases with the increasing \emph{t}
\textsubscript{t}. The presence of allowable \emph{t} \textsubscript{c}
is the result of two competing effects: (i) the increased \emph{t}
\textsubscript{c} can promote the cavity nucleation rate by affecting
the maximum\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} in the next tension phase; (ii)
the prolonged\emph{t} \textsubscript{c} increases the likelihood of
cavity sintering process, given the compressive stress.

Furthermore, the combined effect of \emph{t} \textsubscript{t}
and\emph{t} \textsubscript{c} on the final \emph{f} is presented in Fig.
7b. The distribution of contour curves looks similar by the comparison
of Fig. 7a with 7b, suggesting that the dependence of final \selectlanguage{greek}\emph{ρ}
\selectlanguage{english}and\emph{f} on the load waveform parameters are highly consistent, i.e.
both are positively correlated with the maximum\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n}.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image34/image34}

\end{center}
\end{figure}

Fig. 7: Dependences of the cavity number density \selectlanguage{greek}\emph{ρ} \selectlanguage{english}after 10
cycles in (a) and the cavitated GB area fraction \emph{f} in (b) on the
tensile and compressive hold times of \emph{t} \textsubscript{t}
and\emph{t} \textsubscript{c}

\subsection*{3.4 Effect of stress magnitude in tensile and compressive
holds}

{\label{effect-of-stress-magnitude-in-tensile-and-compressive-holds}}

Fig. 8a presents the effect of \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{t} on
final\selectlanguage{greek}\emph{ρ} \selectlanguage{english}(i.e. calculated after 10 cycles); the values of \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{t} and \emph{t} \textsubscript{t} are varied for this
purpose. Other input parameters are kept the same as those listed in
Table 2. A positive correlation is found between \selectlanguage{greek}\emph{ρ} \selectlanguage{english}and \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{t} for both \emph{t} \textsubscript{t}=100 s and 60 s,
suggesting that a higher tensile stress promotes a higher cavity
nucleation rate. However, the breakdown of this correlation occurs when
\selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{t} is lower than \textasciitilde{}145 MPa under
the condition of \emph{t} \textsubscript{t}=60 s. This sharp decrease in
final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}is related to the effect of compression phase on cavity
sintering. When the tensile stress reduces to lower level, the
associated \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} under a fixed\emph{t}
\textsubscript{t} does not allow the nucleated cavities to grow large
enough to be survived from the following compression phase where the
sintering would occur. As shown in the enlarged view of Fig. 8a, such a
threshold \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{t} appears at the stress level of
\textasciitilde{}105 MPa under the condition of\emph{t}
\textsubscript{t}=100 s. This means that the threshold value for
\selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{t} to generate the cavitation is inversely
proportional to \emph{t} \textsubscript{t} (i.e. 145 MPa for\emph{t}
\textsubscript{t}=60 s as compared with 105 MPa for\emph{t}
\textsubscript{t}=100 s)\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image35/image35}

\end{center}
\end{figure}

Fig. 8: The dependence of final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}on (a) \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{t} and (b) \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c} with\emph{t}
\textsubscript{t}= 60/100 s; the evolution of \selectlanguage{greek}\emph{ρ} \selectlanguage{english}in (c) and
\emph{f} in (d) as a function of fatigue cycles calculated with
different stress magnitudes of \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c}=30/60/140/170
MPa but under the same tensile hold time of \emph{t}
\textsubscript{t}=60 s.

In terms of the effect of \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c} on the
final\selectlanguage{greek}\emph{ρ} , \selectlanguage{english}a positive correlation is found between the two until
\selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c} reaching the breakdown value, Fig. 8b. Such
breakdown of the positive correlation can be seen when \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{c}\textgreater{}120 MPa under the condition of\emph{t}
\textsubscript{t}=60 s. This is followed by a sharp drop of final
\selectlanguage{greek}\emph{ρ} \selectlanguage{english}when \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c}\textgreater{}160 MPa. With the
increase of tensile hold time to \emph{t} \textsubscript{t}=100 s, the
breakdown occurs at a higher compressive stress of when \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{c}\textgreater{}270 MPa.

To explain this interesting phenomenon, the change of \selectlanguage{greek}\emph{ρ}
\selectlanguage{english}and\emph{f} as a function of cycles are presented in Fig. 8c and 8d,
respectively. Four different stress levels of \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{c}=30/60/140/170 MPa are considered, and they represent
different characteristic regions as highlighted by using different
colours in Fig. 8b. Comparing the values of \selectlanguage{greek}\emph{ρ} \selectlanguage{english}at the end of the
1st cycle, the test with \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c}=170 MPa has the
largest \selectlanguage{greek}\emph{ρ} \selectlanguage{english}of 3.6\selectlanguage{ngerman}×10\textsuperscript{-9}mm\textsuperscript{-2},
whereas that with \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c}=30 MPa shows the smallest
\selectlanguage{greek}\emph{ρ} \selectlanguage{english}of 0.7\selectlanguage{ngerman}×10\textsuperscript{-9}mm\textsuperscript{-2}, Fig. 8c.
This indicates that a higher \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c} helps to
produce more creep cavities in the subsequent tension phase. However, in
the next compression phase of the 2nd cycle, these nucleated cavities
would be sintered to different extents depending on the stress level of
\selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c}. For example, the \selectlanguage{greek}\emph{ρ} \selectlanguage{english}value reduction
is found to be -0.2\selectlanguage{ngerman}×10\textsuperscript{-9}, -0.3×10\textsuperscript{-9},
-1.9×10\textsuperscript{-9} and
-3.6×10\textsuperscript{-9}mm\textsuperscript{-2} for \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{c}=30, 60, 140 and 170 MPa, respectively, Fig. 8c. Thus,
the competing effect of \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c} on \selectlanguage{greek}\emph{ρ} \selectlanguage{english}through
the cavity nucleation and sintering leads to the fact that an
intermediate level of \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c} is more likely to
create more cavities.

For the same reason, the breakdown of the positive correlation between
the final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}and \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c} as revealed in Fig.
8b, can be explained by the sintering effect that becomes predominant
when the compressive stress is sufficiently large. In terms of the
cavitated GB fraction \emph{f} , its value increases with the number of
cycles, Fig. 8d. The zigzag shape in the \emph{f} curves reflect the
alternating nucleation and sintering events. In summary, increasing the
stress level of \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{t} is more effective in
enhancing the final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}when compared to \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{c}. The underlying mechanism is the side effect
associated with the cavity sintering due to the increased \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{c}.

\subsection*{3.5 Effect of stress variation
rate}

{\label{effect-of-stress-variation-rate}}

After elaborating the effects of four input parameters (\emph{t}
\textsubscript{t}, \emph{t} \textsubscript{c}, \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{t} and \selectlanguage{greek}Δ\selectlanguage{english}\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{c}), we now consider the
effect of stress variation rate , Fig. 4. To limit the effect of hold
time on the prediction results, both\emph{t} \textsubscript{c} and
\emph{t} \textsubscript{t} are set to 0 s for calculations. Fig. 9a
shows the effect of on the final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}and\emph{f} after 10 cycles.
There is an optimum value of (\textasciitilde{}0.4 MPa/s) that leads to
the maximum final \selectlanguage{greek}\emph{ρ}\selectlanguage{english}and \emph{f} . To explain its presence, the
\selectlanguage{greek}\emph{ρ} \selectlanguage{english}values as calculated at the end of 1st cycle (i.e. the moment
of finishing the tension phase) are presented in Fig. 9b, while the
minimum \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} as calculated in the 2nd cycle (i.e.
the following compression phase) are shown in Fig. 9c. Note that the
magnitude of \selectlanguage{greek}\emph{ρ} \selectlanguage{english}at the end of 1st cycle can be regarded as the
evaluation of the nucleation ability in the considered load-waveform
shape. Equally, the absolute value of the minimum \selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} in the 2nd cycle can be regarded as the evaluation of
the sintering ability.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image37/image37}

\end{center}
\end{figure}

Fig. 9: (a) Dependence of final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}and \emph{f} on stress
variation rate under \emph{t} \textsubscript{c} and\emph{t}
\textsubscript{t}=0 s; (b) \selectlanguage{greek}\emph{ρ} \selectlanguage{english}value at the end of 1st cycle and
(c) the related minimum \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} in the compression
phase of 2nd cycle.

Fig. 9b shows that \selectlanguage{greek}\emph{ρ} \selectlanguage{english}decreases rapidly with the decrease of in
the range of \textless{}0.4 MPa/s. This implies that a too low stress
variation rate reduces the nucleation ability. This seems to agree with
the previous one-cycle creep-fatigue modelling work\textsuperscript{19},
where the local normal stress relaxation was found to occur during the
load reversal with low . On the other hand, the\selectlanguage{greek}\emph{ρ} \selectlanguage{english}curve decreases
with the increase of when \textgreater{}0.4 MPa/s. At first glance, this
might be attributed to the lack of nucleation time when becomes
sufficiently high. However, a closer examination of the change of final
\selectlanguage{greek}\emph{ρ} \selectlanguage{english}in Fig. 9a in comparison with the \selectlanguage{greek}\emph{ρ} \selectlanguage{english}curve in Fig. 9b
does not support the above explanation. It is evident that the final
\selectlanguage{greek}\emph{ρ} \selectlanguage{english}quickly decreases to a value of 10\textsuperscript{-20}
mm\textsuperscript{-2} with increased , Fig. 9a, whereas the \selectlanguage{greek}\emph{ρ} \selectlanguage{english}at
the end of 1st cycle decreases moderately to 10\textsuperscript{-12}
mm\textsuperscript{-2}, Fig. 9b.

We thereby hypothesised that the enhanced cavity sintering effect
contributed to the rapid decrease in the high range. Fig. 9c shows that
the absolute value of the minimum \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} has a
monotonically increasing trend with the increase of . This suggests that
the sintering effect becomes more pronounced in the high range, causing
the final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}to decrease, as revealed in Fig. 9a. To this end, our
hypothesis is supported by the observation.

Finally, it can be seen in Fig. 9b that the optimum does not coincide
with the value of producing the maximum \selectlanguage{greek}\emph{ρ} \selectlanguage{english}value. It is the
combined effect of sintering and nucleation over the course of 10 cycles
that determines the optimum value leading to the highest final
\selectlanguage{greek}\emph{ρ}\selectlanguage{english}in Fig. 9a (for brevity not shown here).

\section*{4. Model applications}

{\label{model-applications}}

\subsection*{4.1 Nucleation incubation
time}

{\label{nucleation-incubation-time}}

An incubation time is required to nucleate cavities \textsuperscript{5}.
The underlying mechanism is the time consumed for vacancies gathering
into the nuclei to form a stable cavity \textsuperscript{40}. In our
model predictions, there is a threshold value of\emph{t}
\textsubscript{t} for creep-fatigue tests without applying compressive
stress hold (\emph{t} \textsubscript{c}=0 s), Fig. 7a. This can be
considered as the equivalent nucleation incubation time. Fig. 10
presents the calculated values for Type 316 stainless steel under
compressive/tensile stresses of (in MPa): -150/+150, -100/+150, and
-150/+100 at different temperatures. The increase of temperature by 50
@\selectlanguage{english}C reduces the incubation time more than one-half. This agrees with the
experimental data trend on dispersion-strengthened copper alloy between
700 and 800 @\selectlanguage{english}C \textsuperscript{41}. If we take the prediction curve of
-150/+150 MPa as the comparison reference, decreasing the tensile stress
by 50 MPa results in doubling the incubation time, Fig. 10. On the other
hand, the incubation time is shortened by 25\%, when the compressive
stress is decreased by 50 MPa.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image38/image38}

\end{center}
\end{figure}

Fig. 10: Model predicted nucleation incubation time (with stress
amplitude and temperature indicated) in comparison with experimental
data. Experimental data points are collected from creep-fatigue tests on
Type 316 stainless steel.

Overall, our cavitation model informs that temperature and tensile
stress are two influencing factors for the nucleation incubation time.
For comparison purposes, experimental results are plotted as red data
points in Fig. 10. The `no cavity' data was obtained from the work by
Shi and Pluvinage \textsuperscript{42}. The creep-fatigue test was
conducted at 500 @\selectlanguage{english}C on Type 316L stainless steel under total strain
range of 1.60\% (stress amplitude of 350 MPa) together with tensile hold
of \emph{t} \textsubscript{t}=10 s. They observed no creep cavitation
damage, being consistent with the transgranular crack propagation mode.
In addition, the 600 @\selectlanguage{english}C `cavity' data point was obtained from the
creep-fatigue test on Type 316 stainless steel conducted by
Hales\textsuperscript{43}. They found cavitation damage under the test
condition of total strain range of 0.5\% and\emph{t}
\textsubscript{t}=60 s. No stress-strain hysteresis loop was given, and
hence the stress amplitude has been calculated as 140 MPa using relevant
database in \textsuperscript{38}. The 593 @\selectlanguage{english}C `cavity' data point was
obtained from the creep-fatigue crack growth test on Type 316 stainless
steel conducted by Michel and Smith \textsuperscript{44}. The test was
performed with the initial stress intensity factor range of 20
MPa.m\textsuperscript{-2}, and they reported the presence of
intergranular creep damage. In summary, our model prediction agrees with
the experimental observation.

\subsection*{4.2 Cycle frequency range}

{\label{cycle-frequency-range}}

Cavitation can also be found in high-temperature fatigue tests under
certain range of cycle frequencies \textsuperscript{45}. The cycle
frequency that is positively correlated with the stress variation rate
can be defined as . The optimum for final \selectlanguage{greek}\emph{ρ} \selectlanguage{english}has been calculated
under the condition of without stress hold (i.e.\emph{t}
\textsubscript{c} and \emph{t} \textsubscript{t}=0 s indicative of pure
fatigue) has been presented in Fig. 9a. Hence, it is possible to extend
our model prediction to examine the nucleation ability under different
cycle frequency range for high-temperature fatigue.

The frequency range for cavity nucleation depends on both the
temperature and stress amplitude. Fig. 11 presents the predicted
nucleation field in the stress amplitude and frequency space at 550 and
600 @\selectlanguage{english}C. Here, the nucleation criterion is defined as the final
\selectlanguage{greek}\emph{ρ}\selectlanguage{english}of [?]10\textsuperscript{-20} mm\textsuperscript{-2}. In other
words, if the condition of \selectlanguage{greek}\emph{ρ}
\selectlanguage{english}\textless{}10\textsuperscript{-20}mm\textsuperscript{-2} is met, the
combination of stress amplitude and frequency is judged as being unable
to create cavitation. The shape of the nucleation field looks like the
upside-down hills. This means that the load waveform with a higher
stress amplitude has a wider range of cycle frequencies to promote
cavitation. Also, an optimum frequency is found at each temperature,
allowing the lowest stress amplitude to meet the cavitation criterion.
The presence of optimum frequency shares the same mechanism as the
optimum in Fig. 9. In low frequency, the lack of nucleation ability is
the reason, whereas enhanced sintering is responsible for the reduced
cavities in the high frequency.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image41/image41}

\end{center}
\end{figure}

Fig. 11: Model predicted cavity nucleation field in stress amplitude and
frequency space for Type 316 stainless steel at 550 and 600 @\selectlanguage{english}C.

It is also interesting to note that the optimum frequency and the whole
frequency range are positively correlated with temperature, Fig. 11.
This can be explained by the time constant \selectlanguage{greek}τ, \selectlanguage{english}Eq.. \selectlanguage{greek}τ \selectlanguage{english}is inversely
proportional to temperature, given that \emph{k} \textsubscript{e} is
temperature independent but both the \selectlanguage{greek}\emph{η} \selectlanguage{english}\textsubscript{b} and are
inversely proportional to temperature. As commented
in\textsuperscript{18}, a lower \selectlanguage{greek}τ \selectlanguage{english}means that a higher stress variation
rate would be required to achieve the same level of\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} leading to cavitation. Therefore, the nucleation
frequency increases with the temperature (i.e. decreasing \selectlanguage{greek}τ\selectlanguage{english}).

Taplin, Tang \textsuperscript{46} and Tang, Taplin
\textsuperscript{47}developed a mechanism map for cavity nucleation in
the stress amplitude and frequency space. For GB sliding induced
cavitation, their predicted nucleation frequency for Type 304 stainless
steel was in the range from 3\selectlanguage{ngerman}×10\textsuperscript{-2} to 1 Hz under
stress amplitude of 180 MPa, and 3×10\textsuperscript{-2} to
5×10\textsuperscript{-1} Hz under stress amplitude of 140 MPa,
respectively. This agrees with our model prediction in a sense that a
higher stress amplitude leads to a wider range of nucleation
frequencies. Also, the frequency ranges obtained at 650 \selectlanguage{english}@C in their
predictions are generally higher than our model predictions at 600 @\selectlanguage{english}C
(3\selectlanguage{ngerman}×10\textsuperscript{-3} to 5×10\textsuperscript{-2} Hz), Fig. 11.
Thus, our model prediction successfully captures the positive
correlation between the nucleation cycle frequency and temperature.

\subsection*{4.3 Unequal ramp rate}

{\label{unequal-ramp-rate}}

In the past two decades, high-temperature fatigue testing with unequal
ramp rate has emerged as alternative method to study cavitation damage
during creep-fatigue interaction, as demonstrated in\textsuperscript{5}.
Yamaguchi and Kanazawa \textsuperscript{48} studied effects of unequal
ramp rate on the strain-controlled high-temperature fatigue in Type 316
stainless steel at 600 \selectlanguage{english}@C. It was found that the test with a slow-fast
cycling waveform exhibited an intergranular fracture mode, and the
fatigue life was reduced significantly due to creep cavitation. By
comparison, the test with a fast-slow waveform fractured in a
transgranular manner, and limited fatigue life reduction was found.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image43/image43}

\end{center}
\end{figure}

Fig. 12: (a) Applied stress as a function of time for the fast-slow and
slow-fast creep-fatigue load waveforms; (b) predicted time evolution
of\selectlanguage{greek}\emph{ρ} \selectlanguage{english}and the time integral of \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s} at 600
@\selectlanguage{english}C.

Now let's consider the unequal ramp rate scenarios from the modelling
perspective. Fig. 12a presents the time-evolution curves of\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{s} for both the fast-slow and slow-fast cycling, and they
were used as the model input. The corresponding \selectlanguage{greek}\emph{ρ} \selectlanguage{english}curves in black
are shown in Fig. 12b. It can be seen that the \selectlanguage{greek}\emph{ρ} \selectlanguage{english}in fast-slow
test gradually increases after a few cycles, whereas that of the
slow-fast test accumulates negligible number of cavities. This is in
conflict with the experimental results as the failure of slow-fast test
was dominated by creep cavitation. To reconcile this seemingly
contradictory phenomenon, it is important to recollect the concept of
our cavity nucleation and early-stage growth model.

Schematic diagram in Fig. 13 depicts the whole process of creep
cavitation. The respective nucleation and subsequent early-stage growth
in Fig. 13a and 13b are the main focus of the present work. Since the
nucleated cavities are extremely small (\textless{}0.1 \selectlanguage{greek}μ\selectlanguage{english}m) compared with
the particles at GB, they are surrounded by a highly localised and
time-dependent stress field \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}, induced by the
GB sliding with the rate of , Fig. 13b. However, the growth mechanism of
large sized cavities would be completely different. Fig. 13c shows that
the large cavity is surrounded by a uniform stress field under the
far-field stress \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s}. The cavity coalescence
shown in Fig. 13d is also affected by the stress concentration, but this
stress concentration is due to the reduction in effective load-bearing
area \textsuperscript{49}, which approximately equals\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{s}/(1-\emph{f} ) \textsuperscript{50}. Therefore, it can
be considered that the local normal stress\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} is
key to early-stage creep cavitation, whereas the far-field stress
\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s} controls its late-stage.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image45/image45}

\end{center}
\end{figure}

Fig. 13. A schematic diagram of cavitation with particular emphasis on
the difference between early- and late-stage cavity growth: (a)
nucleation; (b) early-stage; (c) late-stage and (d) coalescence.

In the late-stage cavitation, damage increase is manifested by the
radius increase of large sized cavities. According to the traditional
growth model, the growth rate is positively correlated with the cavity
radius \textsuperscript{24, 51, 52}. Also, the late-stage growth rate is
controlled by applied stress \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s}. Therefore, the
time integral of \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s} can be used to assess the
propensity of late-stage growth under a given load-waveform shape. In
this sense, a positive value of the integral of\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{s} indicates that the cavities would continue growing up,
whereas a negative value indicating shrinkage in size. Bearing the
difference between the early- and late-stage in mind, the\selectlanguage{greek}\emph{ρ} \selectlanguage{english}curves
in Fig. 12b for the fast-slow and slow-fast load-waveforms can be
reconciled as below.

The time integral of \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s} (red curves) for the two
scenarios are superimposed in Fig. 12b with their values indicated by
the secondary y-axis. The slow-fast one shows a positive value of
the\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s} integral, whereas the fast-slow one shows
a negative value. This implies that cavities in the case of fast-slow is
difficult to grow up under the negative value of\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{s} integral. By contrast, cavities formed in the
slow-fast cycling is more likely to continue growing up, although their
number density is much less compared to that in the fast-slow one. This
eventually leads to an increased cavitated GB fraction \emph{f} in the
case of slow-fast. The cavity coalescence occurs when \emph{f} becomes
sufficiently large, causing the formation of intergranular cracks in
macroscopic scale. Hence, cavitation associated intergranular fracture
was most likely to occur in the case of slow-fast compared to the
fast-slow, as experimentally observed in \textsuperscript{5, 48}.

\section*{5. Conclusions}

{\label{conclusions}}

A mechanistic based cavitation model that considers nucleation,
early-stage growth and sintering under multi-cycle creep-fatigue
interaction has been developed to provide predictions of both the number
density of cavities and the area fraction of cavitated grain boundaries
in Type 316 stainless steel. The model replies on the local normal
stress to connect the three concurrent or sequential cavitation events,
providing important insights about how to design the creep-fatigue load
waveform so that the creep cavitation can be enhanced. The following
conclusions can be reached:

\begin{enumerate}
\tightlist
\item
  When tensile hold time is longer than that of compression, the
  creep-fatigue test with initial compression is advantageous in terms
  of creating more cavities.
\item
  Unbalanced stress hold time in favour of compression most likely
  closes all of the cavities nucleated during the previous tension
  phase.
\item
  Effect of compressive hold time or stress level on the number density
  of cavities in creep-fatigue tests is not monotonic. The underlying
  mechanism is their competing effect on nucleation and sintering.
\item
  There is an optimum value for the stress variation rate to obtain the
  highest number density of cavities, and their presence can be
  attributed to sintering effect.
\item
  Our model can satisfactorily explain several interesting experimental
  observations, including nucleation incubation time, effective
  nucleation field in a stress-frequency space for high-temperature
  fatigue, and the effect of unequal ramp rate on cavitation.
\end{enumerate}

\section*{Acknowledges}

{\label{acknowledges}}

Bo Chen acknowledges financial supports by the UK's Engineering and
Physical Sciences Research Council, EPSRC Early Career Fellowship Scheme
{[}EP/R043973/1{]} as well as East China University of Science and
Technology through the Ministry of Education of the People's Republic of
China `the 111 Project'. Jingdong Hu is funded partly by China
Scholarship Council (CSC) {[}No.201906740075{]} through the PhD Exchange
Scheme. Bo Chen extends his sincere gratitude to Prof Shan-Tung Tu for
providing academic mentorship and endless support.

\section*{Appendix A. Equations for deriving the local normal stress in
Eq.}

{\label{appendix-a.-equations-for-deriving-the-local-normal-stress-in-eq.}}

To derive Eq.(2) for the local normal stress \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n},
let's consider the `spring-dashpot' framework adopted in the Chen-Hu
model \textsuperscript{19}, Fig. A1. The left arm describes the stress
applied on a grain boundary (GB), and the right arm describes the stress
applied on the neighbouring grains.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image46/image46}

\end{center}
\end{figure}

Fig. A1 Spring-dashpot model \textsuperscript{19}

The local normal stress \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n} is equivalent to the
stress concentration at GB particles, which gives:

where is the GB sliding rate. \selectlanguage{greek}\emph{η} \selectlanguage{english}\textsubscript{p} is damping
coefficient for the Norton power-law creep with the deformation rate as
formulised below:

\par\null

where \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{sa} is the stress shared in the
neighbouring grains. An intermediate parameter, transient damping
coefficient , has been introduced to reduce the complexity when
calculating the . Then, at the next time \emph{t} +\selectlanguage{greek}Δ\selectlanguage{english}\emph{t} can be
calculated by \selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{s}(\emph{t} +\selectlanguage{greek}Δ\selectlanguage{english}\emph{t} )
and\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}(\emph{t} ) from previous time \emph{t} :

For a sufficiently small time increment \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} , can be approximated
as a constant. Thus, the differential equation for\selectlanguage{greek}\emph{σ}
\selectlanguage{english}\textsubscript{n} can be derived as:

\par\null

Basically, Eq. is the solution of Eq., and that explicit form is more
convenient for numerical calculations.

\section*{\texorpdfstring{Appendix B. Surface diffusion factor\emph{D}
\textsubscript{S} and the derivation of
Eq.(4)}{Appendix B. Surface diffusion factorD S and the derivation of Eq.(4)}}

{\label{appendix-b.-surface-diffusion-factor-ds-and-the-derivation-of-eq.4}}

Free surface diffusion coefficient \emph{D} \textsubscript{S} is used
only in Eq.(4) to describe the crack-like cavity growth rate via the
diffusivity ratio of\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S}\emph{D}
\textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english} \textsubscript{B}=0.001. The condition that
determines the low or high stress, as reported in\textsuperscript{34},
is formulised as:

\par\null

A sensitivity study was conducted at 550 @\selectlanguage{english}C to examine the effect of
the\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S}\emph{D} \textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english}
\textsubscript{B}value selection on the cavity growth rate. Fig. B1
shows the calculated relationship between the local normal stress
\selectlanguage{greek}\emph{σ} \selectlanguage{english}\textsubscript{n}and cavity growth rate by using
different\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S}\emph{D} \textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english}
\textsubscript{B}values. Only the\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S}\emph{D}
\textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english} \textsubscript{B}values of less than 1 need
to be considered for the crack-like cavity growth \textsuperscript{34}.
Within the stress range of 100 to 1000 MPa, the calculated curve
with\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S}\emph{D} \textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english}
\textsubscript{B}=0.001 (in red) shows the highest cavity growth rate,
Fig. B1. The solid curves as predicted using Eq. are not continuous. A
modification was made to the function (see Eq.), to provide a continuous
curve prediction (the dashed curve). The stress range highlighted by red
box indicates the magnitude facilitating the cavity growth, according to
Figs. 5 and 6. For the maximum cavitation under creep-fatigue,\selectlanguage{greek}\emph{δ}
\selectlanguage{english}\textsubscript{S}\emph{D} \textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english}
\textsubscript{B}=0.001 has thus been selected to derive the value of
\emph{D} \textsubscript{S}.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image56/image56}

\end{center}
\end{figure}

Fig. B1 Effect of diffusivity ratio\selectlanguage{greek}\emph{δ} \selectlanguage{english}\textsubscript{S}\emph{D}
\textsubscript{S}/\selectlanguage{greek}\emph{δD}\selectlanguage{english} \textsubscript{B}on the cavity growth rate.
Dashed curves were predictions using Eq., while solid curves were based
on Eq..

\section*{Appendix C. Numerical error}

{\label{appendix-c.-numerical-error}}

The forward Euler method creates a limited numerical calculation error,
if the time interval \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} is small enough. This method has been
used in previous work on cavitation, e.g. Sanders,
Dadfarnia\textsuperscript{52}, Margolin, Gulenko \textsuperscript{53}.
Fig. C1 presents calculation results on the basis of different \selectlanguage{greek}Δ\selectlanguage{english}\emph{t}
values. It is clear that \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} =0.5 s generates a marginal
difference (less than 2.5\%) when compared with those using any smaller
values (e.g. \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} =0.1, 0.05 s). Thus the selection of \selectlanguage{greek}Δ\selectlanguage{english}\emph{t}
=0.5 is justified.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image57/image57}

\end{center}
\end{figure}

Fig. C1 Effect of \selectlanguage{greek}Δ\selectlanguage{english}\emph{t} value on the numerical calculation results
\selectlanguage{greek}Δ\selectlanguage{english}\emph{t} =0.5, under the condition of 550 @\selectlanguage{english}C, stress amplitude of 200
MPa with initial tension, zero mean stress, and tensile hold time of 60
s.

\section*{Reference}

{\label{reference}}

1. Hales R (1980) A quantitative metallographic assessment of structural
degradation of type 316 stainless steel during
creep-fatigue.\emph{Fatigue \& Fracture of Engineering Materials \&
Structures} .\textbf{3} : 339-356.

2. Murty KL, Charit I (2008) Structural materials for Gen-IV nuclear
reactors: Challenges and opportunities. \emph{Journal of Nuclear
Materials} . \textbf{383} : 189-195.

3. Kassner ME, Hayes TA (2003) Creep cavitation in
metals.\emph{International Journal of Plasticity} . \textbf{19} :
1715-1748.

4. Skelton RP, Gandy D (2008) Creep--fatigue damage accumulation and
interaction diagram based on metallographic interpretation of
mechanisms. \emph{Materials at High Temperatures} . \textbf{25} : 27-54.

5. Rodriguez P, Rao KBS (1993) Nucleation and growth of cracks and
cavities under creep-fatigue interaction. \emph{Progress in Materials
Science} . \textbf{37} : 403-480.

6. Davanas K (2020) Determination of creep cavity nucleation
rates.\emph{Materials at High Temperatures} . \textbf{37} : 75-80.

7. Nam SW (2002) Assessment of damage and life prediction of austenitic
stainless steel under high temperature creep-fatigue interaction
condition. \emph{Materials Science and Engineering A} . \textbf{322} :
64-72.

8. Cailletaud G, Nouailhas D, Grattier J, et al. (1984) A review of
creep-fatigue life prediction methods: identification and extrapolation
to long term and low strain cyclic loading. \emph{Nuclear engineering
and design} . \textbf{83} : 267-278.

9. Wen J-F, Srivastava A, Benzerga A, Tu S-T, Needleman A (2017) Creep
crack growth by grain boundary cavitation under monotonic and cyclic
loading. \emph{Journal of the Mechanics and Physics of Solids}
.\textbf{108} : 68-84.

10. Barbera D, Chen H, Liu Y (2017) Advances on creep--fatigue damage
assessment in notched components. \emph{Fatigue \& Fracture of
Engineering Materials \& Structures} . \textbf{40} : 1854--1867.

11. Wilkinson DS (1978) The Mechanisms of Pressure sintering. University
of Cambridge.

12. Stevens RA, Flewitt PEJ (1983) The sintering of creep-induced
cavities in a low alloy ferritic steel (1Cr1Mo0.75V).\emph{Metallurgical
Transactions A} . \textbf{14} : 679-686.

13. Greenwood JN (1952) Intercrystalline cracking of
metals.\emph{Journal of the Iron and Steel Institute} . \textbf{171} :
380-380.

14. Raj R, Ashby MF (1975) Integranular fracture at elevated
temperature. \emph{Acta Metallurgica} . \textbf{23} : 653-666.

15. Evans HE (1984) \emph{Mechanisms of creep fracture} . Elsevier
Applied Science, Essex, UK, 319.

16. Chen CW, Machlin ES (1957) On a mechanism of high temperature
intercrystalline cracking. \emph{JOM} . \textbf{9} : 829-835.

17. Dunand DC, Han BQ, Jansen AM (1999) Monkman-grant analysis of creep
fracture in dispersion-strengthened and particulate-reinforced
aluminum.\emph{Metallurgical and Materials Transactions A} . \textbf{30}
: 829-838.

18. Min BK, Raj R (1978) Hold-time effects in high-temperature
fatigue.\emph{Acta Metallurgica} . \textbf{26} : 1007-1022.

19. Hu J-D, Xuan F-Z, Liu C-J, Chen B (2021) Modelling of Cavity
Nucleation under Creep-fatigue Interaction. \emph{Mechanics of
Materials} . \textbf{156} : 103799.

20. Cocks ACF (1989) Inelastic deformation of porous
materials.\emph{Journal of the Mechanics and Physics of Solids} .
\textbf{37} : 693-715.

21. Cocks ACF, Ashby MF (1982) Creep fracture by coupled power-law creep
and diffusion under multiaxial stress. \emph{Metal Science} .\textbf{16}
: 465-474.

22. Needleman A, Rice JR (1980) Plastic creep flow effects in the
diffusive cavitation of grain boundaries. \emph{Acta Metallurgica}
.\textbf{28} : 1315-1332.

23. Chuang T-J, Kagawa KI, Rice JR, Sills LB (1979) Non-equlibrium
models for diffusive cavitation of grain interfaces. \emph{Acta
Metallurgica} . \textbf{27} : 265-284.

24. Chen I-W, Argon AS (1981) Diffusive growth of grain-boundary
cavities. \emph{Acta Metallurgica} . \textbf{29} : 1759-1768.

25. Chen I-W (1983) Mechanisms of cavity growth in creep. \emph{Scripta
Metallurgica} . \textbf{17} : 17-22.

26. Miller DA, Hamm CD, Phillips JL (1982) A mechanistic approach to the
prediction of creep-dominated failure during simultaneous
creep-fatigue.\emph{Materials Science and Engineering} . \textbf{53} :
233-244.

27. Tahir F, Liu Y (2017) A new experimental testing method for
investigation of creep-dominant creep-fatigue interaction in Alloy 617
at 950~\selectlanguage{ngerman}°C. \emph{International Journal of Pressure Vessels and Piping}
.\textbf{154} : 75-82.

28. Fournier B, Sauzay M, Caes C, et al. (2008)
Creep--fatigue--oxidation interactions in a 9Cr--1Mo martensitic steel.
Part I: Effect of tensile holding period on fatigue
lifetime.\emph{International Journal of Fatigue} . \textbf{30} :
649-662.

29. Chen B, Flewitt PEJ, Cocks ACF, Smith DJ (2015) A review of the
changes of internal state related to high temperature creep of
polycrystalline metals and alloys. \emph{International Materials
Reviews} . \textbf{60} : 1-29.

30. Frost HJ, Ashby MF (1982) \emph{Deformation Mechanism Maps: The
Plasticity and Creep of Metals and Ceramics} . Pergamon Press, Oxford,
UK.

31. Mortimer DA, Nicholas MG (1976) Surface and grain-boundary energies
of AISI 316 stainless steel in the presence of boron. \emph{Metal
Science} . \textbf{10} : 326-332.

32. Wang JS, Stephens JJ, Nix WD (1985) A statistical analysis of cavity
nucleation at particles in grain boundaries. \emph{Acta Metallurgica}
.\textbf{33} : 1009-1021.

33. Stroh AN (1957) A theory of the fracture of metals. \emph{Advances
in Physics} . \textbf{6} : 418-465.

34. Nix WD, Yu KS, Wang JS (1983) The effects of segregation on the
kinetics of intergranular cavity under creep
conditions.\emph{Metallurgical Transactions A} . \textbf{14} : 563-570.

35. Riedel H (1987) \emph{Fracture at high temperatures} . Springer,
Berlin, Heidelberg.

36. Fedelich B, Owen J (2009) Creep damage by multiple cavity growth
controlled by grain boundary diffusion. \emph{12th International
Conference on Fracture} , Ottawa, Canada.

37. Chen B, Hu JN, Flewitt PEJ, et al. (2015) Effect of thermal ageing
on creep and oxidation behaviour of Type 316H stainless
steel.\emph{Materials at High Temperatures} . \textbf{32} : 592-606.

38. ASME (2015) ASME BPVC III-Div.1 Subsection NH. \emph{Class 1
Components in Elevated Temperature Ser vice} . ASME, New York.

39. Chen B, Flewitt PEJ, Smith DJ, Jones CP (2011) An improved method to
identify grain boundary creep cavitation in 316H austenitic stainless
steel. \emph{Ultramicroscopy} . \textbf{111} : 309-313.

40. Raj R (1978) Nucleation of cavities at second phase particles in
grain boundaries. \emph{Acta Metallurgica} . \textbf{26} : 995-1006.

41. Fleck RG, Taplin DMR, Beevers CJ (1975) An investigation of the
nucleation of creep cavities by 1 MV electron microscopy. \emph{Acta
Metallurgica} . \textbf{23} : 415-424.

42. Shi HJ, Pluvinage G (1994) Cyclic stress-strain response during
isothermal and thermomechanical fatigue. \emph{International Journal of
Fatigue} . \textbf{16} : 549-557.

43. Hales R (1980) A quantitative metallographic assessment of
structural degradation of type 316 stainless steel during
creep\selectlanguage{english}-fatigue.\emph{Fatigue \& Fracture of Engineering Materials \&
Structures} .\textbf{3} : 339-356.

44. Michel DJ, Smith HH (1980) Accelerated creep-fatigue crack
propagation in thermally aged type 316 stainless steel. \emph{Acta
Metallurgica} . \textbf{28} : 999-1007.

45. Ghafouri SN, Faulkner RG, Chung TE (1986) Microstructural
developments in type 316L stainless steel during low-cycle fatigue at
350--550degC and their effects on cyclic strength and life.\emph{Materials
Science and Technology} . \textbf{2} : 1223-1232.

46. Taplin DMR, Tang NY, Leipholz HHE (1984) On
fatigue-creep-environment interaction and the feasibility of fatigue
maps. \emph{6th International Conference on Fracture} . Elsevier, New
Delhi, India, 127-142.

47. Tang NY, Taplin DMR, Plumtree A (1985) Schema for depicting cavity
nucleation during hightemperature fatigue. \emph{Materials Science and
Technology} . \textbf{1} : 145-151.

48. Yamaguchi K, Kanazawa K (1980) Effect of strain wave shape on high
temperature fatigue life of a type 316 steel and application of the
strain range partitioning method. \emph{Metallurgical Transactions A}
.\textbf{11} : 2019-2027.

49. Morin L, Leblond J-B, Benzerga AA, Kondo D (2016) A unified
criterion for the growth and coalescence of microvoids. \emph{Journal of
the Mechanics and Physics of Solids} . \textbf{97} : 19-36.

50. Cocks ACF, Ashby MF (1980) Intergranular fracture during power-law
creep under multiaxial stresses. \emph{Metal Science} . \textbf{14} :
395-402.

51. Hu J-D, Xuan F-Z, Liu C-J (2018) A void growth model of multiaxial
power-law creep rupture involving the void shape
changes.\emph{International Journal of Mechanical Sciences} .
\textbf{144} : 723-730.

52. Sanders JW, Dadfarnia M, Stubbins J, Sofronis P (2017) On the
fracture of high temperature alloys by creep cavitation under uniaxial
or biaxial stress states. \emph{Journal of the Mechanics and Physics of
Solids} . \textbf{98} : 49--62.

53. Margolin BZ, Gulenko AG, Kursevich IP, Buchatskii AA (2006) Modeling
for fracture in materials under long-term static creep loading and
neutron irradiation. Part 1. A physico-mechanical model. \emph{Strength
of Materials} . \textbf{38} : 221-233.

\selectlanguage{english}
\FloatBarrier
\end{document}

