Multi-scale modelling of creep cavity nucleation and growth in polycrystalline Type 316 stainless steel

ABSTRACT A common creep damage mode in Type 316 stainless steel under high-temperature power plant conditions is intergranular cavitation. A review of the literature has confirmed that cavitation in Type 316 is controlled by nucleation, which is not fully understood. In order to provide further insights into the physics of this process, existing strain-based empirical and stress-based (classical nucleation theory) nucleation models are modified in this study by considering experimentally-observed features of cavity nucleation in Type 316. The models are employed locally within a newly-developed crystal plasticity finite element (CPFE)-interface element framework. Modelling results suggest that the strain-based model as a function of local inelastic strain rate does not explain the physical nature of the nucleation process as observed experimentally. By contrast, the modified classical nucleation theory is able to capture features of the observed macroscopic failure response and the distribution of cavities in the microstructure. A number of missing features are identified in the mechanistic model, which need to be incorporated in future unified cavity nucleation theories. These findings highlight key aspects of the nucleation process, which need to be examined experimentally.

Volumetric factor of nucleated cavity/particle geometry j s , j Fitting parametershardening k Boltzmann constant L d , L p , L s Spacing of dislocation junctions, precipitates and solute dispersion atoms DL r Annihilated dislocation segment length N d , N 0 , N p , N s , N li Current/initial number density of forest junctions, precipitates, solute atoms, prismatic loops Number density of GB particles with f v (c ′ ); Number density of cavities and cavitation rate r m p , r p Intragranular and intergranular particle radius t, t f , t red Time and failure time; Redistribution time for unloading particle-free regions along boundaries T, T tk , T n , T e , T nuc Temperature; Tangential, normal, effective traction; Vacancy condensation stress v Poisson's ratio W c Fitting parameterstatic recovery a 0 , a d , a p , a s Obstacle bypassing strength; dislocation junction, precipitate, solute strengths a n , a t , a ′ Coupling parameters for normal and shear tractions; Nucleation constant g s ,ġ max Lattice surface energy per unit area; Maximum local resolved shear strain rate u gb Orientation of the grain boundary with respect to loading axis l p , L Spacing of intergranular particles; Diffusive length m Mean of log-normal distribution of f v (c ′ ) s ij , S,s Stress tensor; macroscopic stress tensor; Standard deviation of log-normal distribution of f v (c ′ ) t, t cr , t r m , t * , t 1. Introduction

Motivation
Under typical high-temperature power plant operating conditions (e.g. AGR nuclear plants), intergranular creep cavitation, leading to microcracking, is the common failure mode in polycrystalline Type 316 stainless steels [1,2]. Brittle creep fracture processes of this type can lead to unexpected plant closures. The main stages of the cavitation process involve nucleation, growth and coalescence of cavities. The process of cavity growth is well understood [3][4][5]. However, the detailed mechanisms of cavity nucleation are still not fully understood. Experimental studies demonstrate that the creep ductility is sensitive to the nucleation of cavities, as well as the evolution of secondphase particles along the grain boundaries [6][7][8][9]. Furthermore, intergranular creep fracture is a process of sequential failure along different grain boundaries, with the local stress state controlled by the deformation of the surrounding grains [1,10]. This suggests that creep failure is intimately linked to the deformation process. A number of damage assessment approaches employed in industry are based on measures of inelastic deformation, e.g. ductility exhaustion models [11]. These practical approaches do not explicitly account for any microstructural features, which influence and affect the process of cavity formation and onset of failure. Not accounting for such features often leads to a significant degree of conservatism in component life assessment. The development of robust physically-informed models for creep deformation and damage is required to provide further understanding of the process of intergranular creep failure, as well as its sensitivity to microstructural features. In this study, the multi-scale mechanistic modelling framework for polycrystalline materials under creep conditions described in [12] is extended to model cavity nucleation and growth along grain boundaries. The model framework is then applied to study the creep deformation and failure of polycrystalline Type 316 stainless steel which is a structural material used in existing and future high-temperature power plant designs.

Experimental observations of cavitation in Type 316 stainless steel
This section gives a short overview of intergranular creep cavitation in Type 316 under the temperature range of interest here: 500-700°C. This serves to inform the modelling activities in this study, where focus is on the short-term creep response of ex-service 316H stainless steel, for which extensive creep and cavitation data is available. Experimental studies on ex-service material after operation in plant for >65,000 h at temperatures 500-625°C and stresses in the range 10-300 MPa reveal creep cavitation to be associated with second-phase particles [1,2,13,14]. Recent findings from [1,9,15,16] of ex-service 316H at 550°C suggest that creep failure is dominated by intergranular microcracking, the precursor for which is the nucleation and growth of creep cavities. Chen et al. [6] suggest that creep cavity formation in Type 316H stainless steel is a result of stress concentrations at obstacles on grain boundaries due to localised inelastic deformation. Type 316H stainless steel, which is extensively used in the UK's AGR plant fleet, contains higher amounts of carbon and is more prone to carbide precipitation than lower carbon grades of Type 316. Ex-service Type 316H contains pre-existing populations of second-phase particles, mostly inter-and intragranular M 23 C 6 carbides. Creep cavities in this material were found mainly at intergranular M 23 C 6 carbides with an average particle diameter of ∼400 nm (min. 100 nm to max. 1500 nm). Lenticular cavities with average diameters between 1 and 3 μm were also noted at ferrite particles in Type 316H [17], but the role of ferrite is negligible since the spacing of this phase (∼10 μm) is too large to lead to interlinkage of cavities. Under certain conditions cavities are found at the intermetallic σ-phase [18]. The σ-phase has not been observed in Type 316/316H stainless steel at temperatures <600°C and creep exposures of <100,000 h [1,2,6,13,14,19]. Evidence from reference [18] shows that σ-phase does not form in Type 316H until after 100,000 h at 550°C. Based on the above observations, σ-phase will be ignored as potential nucleation sites in this study. Observations of cavity morphology across exservice 316H suggest that cavities are of approx. spherical shape with diameters between 200 and 1500 nm [1,2,7,9,13,14,17,20]. The nucleation of creep cavities is a continuous process with time [9], which affects growth and complicates the physical explanation of damage development, as discussed by Dyson [21] and Cho et al. [22]. It is also difficult to determine where the nucleation process commences, since a population of cavities had grown substantially in the above studies where samples were extracted from ex-service components. The studies described here suggest that cavities form due to development of high stress at intergranular particles as a result of grain-boundary sliding or localised slip near the particles. Similar findings were made in [19,23,24]. Since cavities are usually found at second-phase carbides in ex-service Type 316H, these locations will be treated as the dominant nucleation sites. However, experimental studies of cavities usually examine voids that have already undergone substantial growth, which does not aid the understanding of how they form. None of these studies successfully identified the mechanism of cavity nucleation, posing challenges to the development of physical models. This is discussed further in Section 2.

Modelling of intergranular damage development
As noted earlier, damage development on grain boundaries can be divided into a number of stages: cavity nucleation, growth, and coalescence of cavities to form microcracks on discrete grain boundaries, which then interact and link to cause failure. Since the seminal work on diffusive cavity growth by Hull and Rimmer [25] there have been extensive studies on cavity growth and coalescence. These growth and coalescence processes are now reasonably well understood [26] and there is general agreement on the governing processes and constitutive relationships describing them. There is much less agreement in the literature on the factors that control the nucleation of cavities. Models that have been developed generally fall within one of two categories: (a) empirical models based on the correlation of experimental observations; and (b) theoretical models based on classical nucleation theory.
The earliest empirical nucleation model is that by Greenwood [27] who examined the creep response of copper and determined what the rate of nucleation of cavities needed to be to predict the experimentally-determined failure time when combined with the cavity growth model from [25]. This provided a nucleation rate that depended on the (i) applied stress and (ii) creep strain rate. Later, Dyson [21] observed that for a wide range of materials the rate of nucleation of cavities is directly proportional to the creep strain rate. More recently, He and Sandstrom [28,29] have developed a creep strain-based nucleation model which also includes the role of grain-boundary sliding. Their model captures the general trend in the failure response and supports the strong link between lattice deformation, grain-boundary sliding and grain-boundary damage. Raj and Ashby [30] employed classical nucleation theory to determine the rate of nucleation of cavities on boundaries that lie normal to an applied stress σ (and on particles on these boundaries). The resulting nucleation rate is very sensitive to the applied stress, such that the process is well described in terms of a critical stress for nucleation, below which the rate is so slow that it can be taken as zero for all practical applications, and above which nucleation occurs almost instantaneously as soon as this stress is exceeded. The magnitude of the critical stress depends on the interfacial energies (surface, grain-boundary and particle/matrix interface) at the nucleation site. In practice, unrealistic interfacial energies need to be assumed for the critical nucleation stress to lie within typical ranges of macroscopic applied stress, and/or the stress experienced at the nucleation site is much greater than the local nominal stress. A number of studies have explored how these high stress concentrations can occur in practice, either as a result of grain-boundary sliding or localised deformation (slip) within the adjacent grains [34][35][36]. Despite these efforts there is no generally accepted model for nucleation of cavities that simultaneously explains the experimental observations described above.
A number of studies have been undertaken which combine models of cavity nucleation and growth with models of grain deformation within a polycrystalline framework. This allows tracking cavitation on individual grain boundaries, leading to the formation of facet cracks and how they interact to cause failure. Early examples of such studies include those of Tvergaard [32,33] and van der Giessen [34,35], who examined the effect of grain-boundary sliding on creep failure in hexagonal arrays of grains that deform by power-law creep. These multiscale frameworks have been extended to account for grain anisotropy and irregular shape effects through crystal plasticity finite element (CPFE) models. The CPFE method has been employed to study the high-temperature response of polycrystals [36][37][38]. An example is the work of Wei and Anand [39] who coupled the separation and sliding response of the boundary to the elastic-plastic anisotropic deformation response of the grains, which was described by a crystal plasticity finite element (CPFE) model. In crystal plasticity models, the grain crystallographic orientations (i.e. plastically harder/ softer) can result in significant variations of internal stress and strain within a material, which in turn leads to pronounced variations of cavitation and microcrack development and the pattern of damage development in the material, e.g. [40]. The framework from [39] was extended in [41,42] to couple the deformation of the grains to the grain-boundary response under coupled dislocation-climb and diffusion controlled creep. An extensive review of the development and application multi-scale modelling approaches which consider damage accumulation along grain interfaces and the inelastic deformation of anisotropic grains can be found in Roters et al. [43]. Ma et al. [44] and Simonovski and Cizelj [45] used CPFE modelling to study the effect of slip localisation at grain boundaries and local stress fields as precursors for intergranular damage during monotonic plasticity in polycrystals. Pham et al. [46] computed the local stress and strain fields during creep at intergranular regions in a Ni-based superalloy, and used these as inputs to phenomenological cavitation models. Of interest are the recent studies by Phan et al. [47] and Nassif et al. [48] who incorporated a CPFE deformation model for the grains with models for the grain boundaries to predict the intergranular creep failure of some common structural alloys. In [47], realistic tertiary creep curves of ferritic-martensitic steels were obtained by coupling the cavity nucleation and growth model from [49] with a phenomenological CPFE model.
Key issue identified in all these studies is the sensitivity of the predicted failure response to the nucleation model assumed. None of the latter studies, however, considers mechanistic models for cavity nucleation, which has the potential to provide further understanding of the nucleation process and its sensitivity to the observed morphology of grain boundary defect features linked to the cavitation process. Furthermore, the need is identified that such mechanistic models for cavity nucleation should be employed together with physically-informed micromechanical models for the creep response of the grains. Recent research studies [50][51][52] identify the need to develop robust cavity nucleation theories. In this study, various assumptions and modifications will be probed by application of strain-or stress-based nucleation models at the microscale within a multi-scale modelling framework. Of particular interest is the mechanistic stress-based nucleation model and its modification to account for the distribution of defects. The modelling framework employs physicallyinformed micromechanical models for the creep deformation response of the grains and cavitation process along the grain boundaries. These aspects are discussed in the following section.

Purpose of the present study
The main objectives of this study are: . Modify existing cavitation models to capture microstructural features and local processes affecting the nucleation response, and include them in an extended interface element framework, coupled with a CPFE model for 316H SS. . Evaluate the sensitivity of intergranular cavitation to the detailed form of the cavitation model and how it is influenced by the local stress-strain fields, grain boundary orientation with respect to the load axis and the crystallographic orientations of the neighbouring grains. . Explore the ability of the modified cavitation models to capture general features of the intergranular creep failure process and identify aspects to be considered in the future development of nucleation theories.
This study explores the relationship between assumptions made on the cavitation process at the microscale and the macroscopic failure response, and how this compares with available experimental data. In turn, we seek to identify areas where additional studies (both theoretical and experimental) are required to improve and further enhance the models. This is done by including the intergranular cavitation response within an existing multi-scale CPFE-interface element framework, described in [12], which has been shown to accurately predict the global creep deformation response, the evolution of microstructural state, and sliding behaviour of the grain boundaries of 316 stainless steel. In this study, only intergranular creep failure is modelled for the stress-range (280-320 MPa) based on experimental studies discussed in Section 1.2. Section 2 revisits the creep cavity nucleation theories briefly described in Sections 1.2-1.3. Important experimental observations are evaluated in order to enhance these models in describing the intergranular cavitation processes in Type 316H under the thermo-mechanical conditions of interest. Modified local strain-based (Dyson-type) and stress-based (classical vacancy condensation theory) nucleation models are introduced, which are coupled to diffusive growth models, and implemented within an interface element formalism in Section 2. The constraint at triple grain junctions in a polycrystal introduced in [12] is also modified to account for the potential volumetric opening due to damage along adjacent interfaces. The applicability of the local strain-and stress-based nucleation models in predicting (i) macroscopic failure, (ii) continuous nucleation, and (iii) cavity distributions and damage morphology within the microstructure as observed experimentally is quantified. Model predictions of grain-boundary cavitation, microcrack formation and their interaction leading to creep failure in simple bicrystals and polycrystalline aggregates are presented and analysed in Sections 3 and 4. A discussion of future modelling and experimental activities to support the development of more robust cavity nucleation models can be found in Section 5.

Proposed modifications to existing cavitation models and modelling framework
Despite a considerable research effort over the past 50 years, it is still not understood whether cavity nucleation is strain-or stress-controlled [26,53]. Here, modifications of existing general theories for cavity nucleation are proposed in order to assess this aspect further. The modifications are guided by experimental observations from the literature. An approach is described for extending Dyson's empirical macroscopic strain-based nucleation model to the microscale. Two modifications to the mechanistic, stress-based classical nucleation theory are then explored, by considering the dependence of the cavity nucleation process on particle morphology and local stress. A diffusive cavity growth model is presented and coupled to the nucleation models.

Proposed local Dyson-type nucleation model
It is demonstrated in [21,54,55] that prior plasticity increases the probability of intergranular cavitation during creep and cavity spacing is often found to be similar to the slip band spacing [56,57]. Dyson et al. [58] suggest that nucleation occurs via a Zener-Stroh decohesion mechanism at either grain-boundary particles or slip ledges, with an empirical model describing the process given in [21]. Early in creep life, cavity distributions in [58] showed peaks at boundaries, approx. parallel to the loading axis, similar to observations from [59,60] on Type 316 and Astroloy (FCC). Later in creep life, cavity growth was observed along boundaries transverse to the loading axis, while cavities along parallel boundaries did not readily grow, agreeing with [61]. The authors in [58] reason that higher cavity number densities are limited to boundaries enclosed by plastically-softer grains, where local slip activity in the vicinity of grainboundary obstacles is greater. Only preferentially-orientated boundaries for diffusive growth (i.e. approx. normal to principal load or enclosed by plastically-harder grains, which carry higher stresses) are prone to cavity coalescence. A recent 1-D micromechanical nucleation model, introduced by He and Sandstrom [28], also suggests a relationship between the creep deformation process of the grains, the grain-boundary morphology and the nucleation rate. While the model successfully predicts evolution of cavity densities in austenitic stainless steels, due to its continuum nature it fails to explain the above-mentioned distribution of cavities with respect to grain-boundary orientation and its effect on the subsequent events leading to complete failure. Based on the above, a simplified approach is employed here to relate the nucleation process to the localised slip activity in a grain in the vicinity of a grain boundary. In the macroscopic strain-based nucleation model it is assumed that the nucleation rateṄ h = a ′Ė c e (whereĖ c e is the macroscopic effective creep strain rate and a ′ is a material constant) [21]. In the local model employed here we express the nucleation rate as a function of the maximum resolved shear strain rate, g max , in the grains adjacent to a boundary aṡ We do not assign a physical mechanism to the nucleation processit could be local delamination, i.e. slip activity leading to a Zener-Stroh mechanism at grain-boundary obstacles, as described by Dyson et al. [58] or it could be by vacancy coalescence [30] or some other mechanism. We simply treat this as an empirical model here, that captures some features of the observed experimental behaviour. Dyson [21] suggested that the nucleation constant a ′ can be estimated in terms of local microstructural features, such as particle spacing. In [21], a ′ is suggested to decrease throughout the creep life as (a ′ − a n ) as nucleation sites a n are depleted. Nevertheless, the assumptions that this value is a constant is not unreasonable for Type 316 stainless steel, where intergranular creep failure is considered to be nucleation-controlled [9,62]. Employing Equation (1) within a polycrystalline crystal plasticity framework would assist in interpreting cavity nucleation as a function of crystallographic orientation of the grains, enclosing the boundaries, and the relative local inelastic strain accumulation. Section 4 will assess the suitability of Equation (1) in describing the intergranular failure process.

Proposed modification of classical nucleation theorydistribution of particle strengths
The mechanistic theory for cavity nucleation, originally proposed, by Raj and Ashby [30] is based on the assumption that vacancy condensation leads to the formation of creep voids. The stress driving the vacancy condensation process, T nuc , is the local microscopic stress normal to the nucleation site [26]. In [30], the nucleation rate of voids of critical size (i.e. critical pore radius r c = 2g s /T nuc ) is expressed aṡ where Ω is the atomic volume, d b D b is the grain-boundary thickness times the diffusivity, g s is the lattice surface energy per unit area, c max and c are the concentrations of potential nucleation sites and critical nuclei, f p is a perimeter shape factor and f v is a volumetric factor of the cavity/particle geometry, both expressed in terms of a dihedral angle c ′ , which depends on the relative magnitudes of the surface and interfacial energies per unit area.
Of particular relevance is f v , which takes the functional form f v = 4p/3(2 − 3 cos c ′ + cos 3 c ′ ). In engineering alloys with a population of intergranular particles, c max and c are determined by multiplying the number density of total (N p ) and already-cavitated particles (N h ) by the number of nuclei which can be formed on the particle surface, (r p T nuc /g s ) 2 , respectively. When the classical nucleation theory is employed within a continuum plasticity framework it essentially predicts a threshold stress above which the cavity nucleation rate is rapid and below which it is so slow that it can be ignored [26]. This disagrees with the continuous evolution of cavity density with inelastic strain observed experimentally, e.g. [21]. Detailed observations in [1,13,17,60] reveal that intergranular particles in exservice Type 316 are of highly-irregular shapes with multiple faces. Kim et al. [63,64] further found that planar carbides rather than faceted carbides offered higher cavitation resistance under creep. These findings suggest that the cavity nucleation process is sensitive to the detailed shape of intergranular particles. Furthermore, cavity growth will increase the stress on the remaining non-cavitated particles as the void fraction increases (Figure 1) [65]. This suggests that the boundary contains populations of particles with different nucleation strengths, agreeing with observations from Slater et al. [60]. Therefore, we assume that only when the stress acting on the remaining uncavitated particles has increased substantially, nucleation at the next set of particles with a higher nucleation threshold is possible. The concept of employing such distributions of particles with different nucleation thresholds within the mechanistic nucleation theory appears important in describing the nucleation process [66]. Here, the dihedral angle c ′ which describes the cavity/particle geometry following nucleation is identified as an important microstructural feature, associated with the nucleation energy barrier, agreeing with [26,51,67]. This suggests that a distribution of c ′ , or equivalently f v (c ′ ), for the existing population of particles along a boundary could be incorporated within the mechanistic theory of Equation (2). A log-normal distribution of f v (c ′ ) is assumed to exist along the boundary. The number density of GB particles on which cavities form with volumetric shape factors between f v and f v +δf v is where l p is the average spacing of intergranular particles, and where μ ands are the mean and standard deviation of the log-normal distribution. The total nucleation rate of critical nuclei is then obtained by integration over the rates for each population of particles aṡ Note that c max is a function of N p ( f v ) and will vary for different values of f v . For simplicity, it is assumed here that only one stable cavity can form on a particle, which is often observed in experimental studies, e.g. [13]. Hence, the number density of cavitated particles, N h (or simply number density of cavities in Figure 1) is N h = c h (r p T nuc /g s ) −2 . This modified version of the classical nucleation theory will be employed and evaluated in detail in Sections 3 and 4. Figure 1. Schematic of continuous nucleation at different GB cavity/particles of mean spacing, 2l h . In (a), cavity is nucleated at a lower-nucleation barrier configuration. (b) The initial cavity grows and load redistributes from cavitated particles (lower nucleation barriers) to those with higher nucleation barriers.

Proposed modification of classical nucleation theorylocal stress concentrations
The other feature we need to consider is the degree to which the stress is elevated at second-phase particles compared to that calculated by the continuum crystal plasticity theorywhich provides the material response averaged over a length scale much larger than the particle spacing. There are two major ways in which a higher level of stress can be generated at the particles: (a) diffusional relaxation of the stress across a boundary between particles; and/or (b) the development of dislocation structures in the matrix adjacent to the particles. We consider each of these in turn.

Stress concentration at intergranular particles due to diffusional relaxation
Closed form solutions for maximum local stress at particles due to relaxation of shear stresses and normal stress are given by Argon et al. [31] and Riedel [26], respectively, for cases when the accommodation around the particle is by plasticity and diffusional flow. According to Needleman and Rice [5], coupling between plastic flow and diffusion exists. This coupling is represented by a characteristic diffusive length, are the stress and associated strain rate in the grains either side of the boundary. The significance of this characteristic length is that when 2Λ is greater than the intergranular cavity spacing (2l h ) or particle diameter (d p ), diffusive transport dominates over plastic flow (i.e. moderate s 1 and T). At higher temperatures and for smaller particles, diffusion modifies the stress distribution [31]. Under such conditions, when 2L .. 2r p , the maximum local traction at the particle interface is approximated in [26,31] by t where T i is the average traction carried by the boundary, and subscript i denotes either normal (n) or shear (t) traction component (see Figure 2). Recall that inclined boundaries in FCC materials tend to have a higher number densities of cavities i ) are in equilibrium with the average traction T i (i = n,t). Normal or shear displacement incompatibilities arise along the particle interface, since the interface diffusivity is lower than the grain-boundary diffusivity, e.g. [70]. Slower accommodation of incompatibilities at particles, in highlighted regionsin (c), causes an increase in stress in that region. [59,60], with the highest cavity densities in 316H stainless steel on boundaries orientated at approx. 60°to the principal load direction [60]. This suggests that nucleation is most potent along boundaries, on which the normal and shear stresses acting on particles have the largest combined contribution. This suggests that the stress driving the vacancy condensation process (T nuc ) is in fact a resultant traction, T e = f (T n , T t ). Noting that the solution for t p i is similar for both the shear and normal traction limiting cases in [26,31], the maximum possible local stress t * in the vicinity of grain-boundary particles under the condition 2L .. 2r p is approximated by which is a scalar, since T e is also a scalar. Section 3 will demonstrate that 2L .. 2r p in Type 316 stainless steel over the temperature and stress range of interest. Hence, Equation (6) could be employed within the multi-scale modelling framework, which computes the tractions along grain boundaries in the presence and absence of sliding, to estimate the maximum possible stresses at particles. Note that a timescale, t red , associated with the increase of particle stress due to diffusive accommodation along particle-free regions is discussed in [31,68], and is given by where t b i is the traction along the particle-free region and E is the elastic modulus of the matrix. Details of the derivation of Equation (7) are given in [69]. For t red 0, the maximum stress at particles will be given by Equation (6), whereas for t red 1, the particles will carry the average boundary traction. In Section 3, the effect of this redistribution time and the coupling between the normal (T n ) and shear (T t ) tractions in relation to the local stress state in the vicinity of the boundary are assessed.

Stress enhancement at intergranular particles due to development of dislocation structures
It is believed that lattice dislocation structures, such as pile-ups at intergranular particles or ledges are responsible for cavity formation by increasing the local stress at these locations [21,71,72]. Mechanistic cavity nucleation models need to account for the resulting residual stress at these locations and how they influence the nucleation process. The simple analysis here aims to highlight a feature of the dislocation structure near the boundary, which could result in elevation of local stresses at intergranular particles.
Dislocations in the grains bypass precipitates via an Orowan type of process in Type 316 initially and then through a process that involves the punching out of prismatic loops [73]. This causes the formation and storage of geometrically-necessary dislocations around particles to accommodate local inhomogeneity [74]. These dislocation structures result in an increase in stress on the particles (and a concomitant reduction of stress in the matrix). We assume that a similar process occurs at particles either at or in the vicinity of a grain boundary. Following the approach from [73], the residual stresses at intergranular particles due to the presence of lattice prismatic loops in the vicinity of boundaries t GB p can be evaluated. Measurements of the dislocation density in the vicinity of grain boundaries and the grain interior in Type 316 at >400°C [75] show a difference of less than ∼15%. The simple assumption is made whereby the density of loops, N A li and N B li , and microscopic residual stresses in the lattice, t r m,A and t r m,B , are similar in the grain interior and local to the grain boundary, i.e. regions A and B in Figure 3(a). Considering force equilibrium in Region B in the unloaded state ( t r m,B (1 − f GB a ) + t GB p f GB a = 0) allows the residual stress at GB particles due to the presence of lattice prismatic loops to be associated with inelastic deformation near the boundary, according to˙ is the area fraction of GB particles and t GB p is a scalar that accounts for the on-average density of dislocation loops; other terms are described in Section 2.5. The prismatic dislocation loops can be assumed to contribute to T nuc . The CPFE framework, described in Section 2.5 predicts the evolution ofṄ li during creep within the lattice and the effect of such idealised dislocation features near the boundary can be tracked using Equation (8). This is discussed in Section 2.6. The presence of loops gives rise to microscopic residual stresses at lattice particles ( t r p ) and matrix ( t r m ). The stress at GB particles, t GB p , due to prismatic loops in the vicinity of the boundary (Region B) is estimated by assuming t r m,A = t r m,B . In (b), schematic of the contributions to T nuc is shown.

Diffusive cavity growth model and evolution of damage
Coupling between diffusion-and plasticity-controlled growth of cavities is possible [5,76]. In Type 316 under the thermo-mechanical conditions of interest, the local diffusive distances in the material are greater than both particle size and cavity spacing (see Section 3). Under these conditions, grain-boundary diffusion controls the growth process [26]. The rate of diffusion plating on a boundary, interpreted as an inelastic opening rateḋ cr n can be expressed [3] aṡ where f h is the cavity area fraction on the boundary (f h = r 2 h /l 2 h ; r h is the cavity radius and l h = N −1/2 h /2 is the half-spacing) and T n is the average normal traction at a material point along the boundary. Following the approach from [3], the evolution of cavity area fraction as a function ofḋ cr n is given bẏ where h(c ′ ) is a function of the cavity shape [30] (∼0.6 for spherical cavities). Coalescence of cavities occurs when the area fraction reaches a limiting value, f c h , often taken as 0.25 [3,50]. The models from Equations (9)-(10) are employed in Sections 3 and 4 to describe the cavity growth process. Note that in the model of Equation (9), the cavity growth rate is determined by the local traction across the grain boundary. The magnitude of this traction is determined within the finite element framework and is determined by the creep response and constraint imposed by any surrounding grains. In situations where a cavitated facet is surrounded by undamaged material the model is fully consistent with Dyson's classical constrained cavity growth model [77]. In the following sub-section, the theoretical modelling framework and aspects of the implementation of the cavitation models are presented.

Theoretical framework for modelling grains and grain boundaries
The computational framework employed here extends that described by Petkov et al. [12]. The solution of a boundary value problem of a body, discretised into finite continuum and interface elements with nodal displacements U (vector containing all the degrees of freedom), is given by [78] where B and B c are matrices that relate the local strain and separation (δ) across an interface to the nodal displacements, respectively; σ is the local stress, T c are the interface tractions and F is a vector of nodal forces. Solutions are obtained by using the commercial software ABAQUS. The constitutive response of the grains, and determination of stress and the material Jacobian matrix are represented by the physically-based CPFE model described in [12,79]. The deformation of the anisotropic grains is modelled using crystal plasticity (CP) theory [80], where the plastic strain rate in each crystal is determined by summing the slip rates,ġ, on each individual slip system, α. The assumed flow rule describes the thermally-activated dislocation glide process at a temperature T and resolved shear stress t (a) , and is given by [81] where DF 0 (which is equal to a 0 G 0 b 3 , G 0 : shear modulus at 0 K) is the required free activation energy, k is Boltzmann constant andġ 0 is the reference shear strain rate. The micromechanical model captures the contributions to internal resistance t (a) cr , assumed the same for every three slip systems on a given slip plane i, of the dominant obstacles to dislocation motionforest dislocation junctions, precipitates and solute atomsas where G is the shear modulus, b the Burgers vector, and a p , a d and a s are dimensionless obstacle strength parameters [82]. The number densities for the dislocation junctions, precipitates and discrete solute atoms are N di , N p and N s , respectively. The evolution law for the number density of dislocation junctions, including the hardening, static and dynamic recovery terms, and prismatic loop formation, is: where D c is the dislocation core diffusivity andṄ li is the production rate of prismatic loops. The free parameters j s , W c and DL r for the hardening, static and dynamic recovery models are calibrated according to the procedure described in [12]. (N.B. N p , N s = const since ex-service material is considered). The term t r m is the microscopic residual stress in the matrix due to the presence of prismatic loops [73]. The generation and annihilation of loops on the i-th slip plane during creep are determined by the resolved shear strain rates (ġ i ) and dislocation obstacle density (N di ) on that plane. The production rate of loopsṄ li , which adds a contribution to the overall obstacle density on the i-th plane, isṄ with f p v and r m p being the volume fraction and radius of intragranular precipitates, respectively. A fuller description of the model and its implementation can be found in [12,79].
The surface-type formulation of the interface element from [78], implemented in ABAQUS as a User Element routine, is described in Appendix B. Traction-separation-rate laws of Equations (16)-(17) represent the elasticcreep response of grain-boundary regions. The tangential creep response of the interface is described by the linear viscous law forḋ c t from [12]. Cavity nucleation and growth models have been implemented within the traction-separation-rate laws by combining Equations (9) and (10) to giveḋ c n . The normal and tangential separation rates are then given bẏ where, at a given integration point along the interface, f h is the cavity area fraction, N h is the number density of cavities, andḋ e i andḋ cr i represent elastic and viscous components, respectively. All other terms remain as described in [12]. The tangential and normal stiffnesses are sufficiently high (1 × 10 6 N/mm) to limit the elastic separation along the boundary. Within Equations (16)-(17), f h increases as a function of (i) the inelastic openingḋ cr n (Equation (9)), and (ii) cavity nucleation due to an increase in N h . The latter is described either by the Dyson-type nucleation model (Equation (1)) or the modified stressbased classical nucleation theory (Equation (5)). Cavity coalescence and failure at an integration point occur when f h reaches a critical value of f c h = 0.25 [3], when the interface tractions reduce to zero. The integration procedure for obtaining solutions and updating f h , T n and T tk at the end of each increment is briefly described in Appendix B.

Implementation of local cavity nucleation models
The local strain-based model of Equation (1) is employed as a function of the maximum resolved shear strain-rate in the lattice,ġ max , in the vicinity of a material point along the interface. The modelling framework extracts the maximum value ofġ (a) , computed by the CPFE scheme at increment i, at the two continuum-element integration points either side of an interface integration point (see Figure 4(a)). The higher value of the two is used to compute the current nucleation rate at this interface integration point. Explicit integration is used to determine the cavity number density at the end of each increment, according to N i+1 h is the nucleation rate at the i-th increment and Δt is the time step.
The modified vacancy condensation theory of Equation (5) assumes a lognormal distribution of volumetric shape factors f v (c ′ ). Here, the cavity number density N h at an interface integration point at the end of the i-th increment is determined by explicit integration of Equation (5) as where F n represents all material properties in Equation (5) andċ f v h is the individual nucleation rate of critical nuclei for a given population of cavity/particle geometries with volumetric shape factors between f v (c ′ ) and f v (c ′ + Dc ′ ) (Figure 4(b)). The integral in Equation (18) is solved numerically by integrating individual nucleation rates across the dihedral angles (5), the limiting cases of particles carrying the average resultant boundary traction (T e ; t red 1) and maximum particle stress (T e /(d p /l p ) 2 ; t red 0) are examined in Section 3. This provides insights into the influence of timescale for stress transfer on the nucleation response. Stress at particles can also be enhanced due to the presence of dislocation structures in the vicinity of particles. This process competes with that due to diffusional relaxation. Here, when considering enhancement of stress at particles we assume that diffusional relaxation dominates and employ We use Equation (8) to monitor the potential contribution from the dislocation loops to the stress experienced by the particle but we do not include it in T nuc in this study. The value of N max li , used to determine t GB p in Equation (8), is extracted from the CPFE model as illustrated in Figure 4(c), following the same approach to extractġ max .
Furthermore, T e in Equation (19) is the resultant boundary traction, as discussed in Section 2.3. T e is determined here by coupling of the normal (T n ) and shear (T tk ) tractions, according to where a n and a t determine the contribution from the tractions to the resultant traction. In reality, the values of a n and a t are likely to be a function of the faceted particle geometry, but for the purpose of this study these parameters are selected to explore the relationship between local nucleation condition and observed macroscopic trends in behaviour (see later). In addition, only tensile normal tractions contribute to the resultant traction in Equations (20)sintering of cavities is not considered. Section 4 demonstrates that a unique set of (a n , a t ), associated with the local stress state, captures the expected distribution of cavity number density in FCC materials. The inclusion of grain-boundary cavitation in the framework described in [12] requires modification of the sliding constraint at triple junctions, examined in the following section.

Triple grain junctions and inelastic opening of adjacent boundaries
In [12], artificial opening of the triple point as a result of sliding along adjacent boundaries was restricted in thin-sliced, quasi-3D grain aggregates. This was achieved by the introduction of a triple line element, formulated to restrict the change in area of a fictitious element perpendicular to the line of the element (L s Line = AU = 0). In the present article the damage response of full-3D bicrystals ( Figure 5(a)) and quasi-3D, thin-sliced polycrystalline aggregates ( Figure 5(b)) is examined. Triple lines do not form in the bicrystal, but they do in the polycrystalline geometry. The approach from [12] is extended here by considering the volume change of the triple line element due to opening at the boundaries connected to the triple line. Consider the geometry of the front face of the triple line element in Figure 5(c), where the boundary has a thickness h. The instantaneous rate of area change of the region N in Figure 5(c) due to opening of an individual boundary, taking h 12 as an example, is expressed aṡ where c 12 is the height of the triangle with respect to h 12 andu n i is the nodal displacement rate in a direction, perpendicular to the mid-surface of the adjacent boundary (i.e. direction of opening in the co-ordinate system of the interface element).
As in [12], the approach is applied to triple lines of arbitrary orientation, consisting of triangles with nodes (1-3) and (4-6) (for more detail, see [69]). Dividing by h 13 gives the rate quantity, independent of boundary thickness (not to be confused with the out-of-plane thickness of the polycrystal in Figure 4(b):  [12]. The instantaneous area change due to normal separation along a single interface element, enclosing the triple line, is illustrated in (c) with respect to a reference configuration. The case for opening along all three sides is given in (d).
where primed quantities are associated with the rear face, formed by nodes (4)(5)(6). Interface elements of zero thickness are employed and change in length of the triple line element is allowed. The rate quantityL n Line is a function of the triple line nodal displacements rates in the global co-ordinate system,U, aṡ L n Line = B nU . The matrix B n is given in Appendix C. The constraint in the presence of opening along boundaries, adopted here, requires that the area change due to opening is accommodated by area change due to sliding or opening along adjacent boundaries. The constraint can be expressed aṡ

Approach to modelling and probing effectiveness of cavitation models
The modelling framework from Section 2 is used in Section 3 to study intergranular creep failure of bicrystals of Type 316 stainless steel. Once assessed and examined along the isolated grain-boundary in the bicrystal, the damage models of Equations (1) and (5) are applied in Section 4 to model a polycrystal of ex-service Type 316H stainless steel at 550°C consisting of 39 grains, as shown in Figure 5(b). Effectiveness of cavitation models in predicting the failure process is judged by: (a) agreement with macroscopic creep data in the range 250-320 MPa; (b) ability to predict the macroscopic failure response of both bi-and polycrystalline geometries; and (c) ability to predict typical cavity distributions in FCC materials described in Sections 2.1-2.3, i.e. higher cavity densities along inclined boundaries (∼45-60°to load axis) and dominant failure along transverse boundaries (∼90°to load axis), needs to be captured. Note that the polycrystalline aggregate here is not treated as a RVE. Nevertheless, the study provides insights into the intergranular damage process in polycrystals of anisotropic grains with irregular shapes in relation to stress state, particle morphology and idealised dislocation features near grain boundaries.

Modelling failure of bicrystals of Type 316/316H
The modelling framework is employed to simulate intergranular creep failure of the bicrystal of Figure 5(a) for solution-treated (ST = annealed) Type 316 at 700°C and ex-service Type 316H stainless steels at 550°C. The compositions are given in [69]. Initially, we focus on the response at 700°C since creep strength (t f vs. Σ) data at this temperature for Type 316 in the range 60-160 MPa are available in [62]. A cell of the bicrystal sample from [62] with dimensions 15 × 20 × 10 μm is examined; the actual bicrystal grain size in [62] is ∼ 900-1000 μm. Constant stress is applied on the top face of the cell. Due to ease of implementation, a simple class of boundary conditions (BC) are used throughout this section, whereby the cell faces remain straight. Calibration of the damage models here is based on the measured time to failure t f . It is found that the imposed BC to the bicrystal do not significantly alter the result. Continuum C3D8 elements for the grains, which deform according to the CPFE model from [12], and 8-noded linear interface elements along the grain boundaries (Appendix B) are employed. Five different bicrystal meshes with 1, 20, 80, 216 and 480 interface elements along the boundary were used for a mesh convergence study. The mesh with 80 interface elements was selected, since predictions for macroscopic failure time, and distribution of N h and f h remain unchanged with further refinement (see example in Figure  6(a)). The crystallographic orientations of the grains either side of the boundary are arbitrarily selected (see Figure 6). Note that short-term plasticity and creep deformation data for the bicrystals are not available. For simplicity, their deformation response is calibrated to polycrystalline material deformation data, not presented for the sake of brevity. Based on observations from Section 1.2, only M 23 C 6 intergranular carbides are treated as cavity nucleation sites in this study. Information about the inter-and intragranular carbide morphology is obtained from experimental studies on solution-treated (ST) Type 316 [24,62] and exservice Type 316H [1,2,6,13]. The kinetics of precipitate formation and growth are rapid in ST Type 316 at 700°C [19,24], hence the values of L p , r m p , d p , l p for M 23 C 6 carbides in both ST and EX material are assumed constant and identical for each grain boundary for a given temperature. The parameters L p , r m p , d p , l p are extracted from experimental observations to reflect the specific cast/thermal treatments of the modelled material. Note that recent work by Xiong et al. [83] demonstrated how such parameters can be extracted by simulating the thermal ageing processes in 316H. Note that that such simulations may not capture the evolution of both precipitates/ solutes and the dislocation content accumulated in the ex-service material which was aged under load. Hence, we chose to base sizes of microstructural features on experimental observations here. The model parameters for the CPFE and interface elements are given in Table 1.
In the following sections, the failure response of the bicrystal is simulated for cases where intergranular cavitation is a function of diffusive cavity growth coupled to (i) the local Dyson-type model, and stress-based classical nucleation theory with (ii) a distribution of particle strengths and (iii) a distribution of particle strengths and maximum stress at particles, described in Section 2.

Local Dyson-type nucleation model and diffusive growth
Initially, we focus on Type 316 stainless steel at 700°C. No pre-existing cavities are assumed along the boundary, so that the initial density of cavities is N 0 h = 0. This also implies that f 0 h = 0, but to avoid numerical issues that arise with the use of Equation (16) in this limit, we assign an initial small value to f 0 h . The local Dyson-type nucleation model of Equation (1) is employed in this subsection to describe the evolution of cavity density and is coupled to the diffusive growth model of Equations (9) and (16). The sensitivity of the time to failure on the nucleation constant a ′ was assessed for an initial cavity area fraction of f 0 h = 10 −15 ( f 0 h 0) and failure fraction f c h = 0.25 in Figure 6(a). Note that the failure response was insensitive to the selected initial cavity area fraction  [51]. Values in brackets for bicrystals. Values of C ij , v,ġ 0 , G 0 , a 0 , b, Ω and α d are those from [12,79].
for f 0 h , 10 −2 , agreeing with findings from [86]. Hence, values of N 0 h = 0, f 0 h = 10 −15 and f c h = 0.25 are used in all subsequent sections. For the simplest scenario of N h = const. (i.e. pre-existing population of cavities along the boundaries with no subsequent nucleation), failure is driven by diffusive growth and the creep strength varies linearly with stress. This disagrees with the experimental response in Figure 6. Thus cavity nucleation needs to be considered in order to capture the non-linear trend in experimental data. Predicted failure times for different values of a ′ in the Dyson-type model are presented in Figure 6(a) for N 0 h = 0, and agreement with experimental data is found for a value of a ′ = 8.9 × 10 6 mm −2 . Note that the interface element framework for the Dyson-type model, as well as for diffusive growth only, was validated against analytical results from [3,26]. The validation procedure is not presented here. Cavities could also grow in a non-linear manner, controlled by creep of the surrounding matrix [4]. When Λ ≫ d p and L ≫ 2l h (cavity spacing), the contribution to cavity growth due to plasticity/creep is minimal. Following from the discussion in Section 2, the mean diffusive lengths Λ in Type 316 are 40.0-4.2 μm (60-160 MPa, 700°C) and 4.2-2.2 μm (250-320 MPa, 550°C). These are larger than values for d p from Table 1 and 2l h predicted by the model here with typical values of 0.7-1.8 μm (see later). Throughout this article, only diffusive cavity growth is considered.
The nucleation model from Equation (1) is dependent on the local strain field in adjacent grains. In Figure 6(b), the predicted failure response of 5 bicrystals with grains of different crystallographic orientations is simulated for a ′ = 8.9 × 10 6 mm −2 . Failure shifts to longer times for samples with plastically-harder grains with orientations away from the [001]-direction (e.g. Ori set 5/6) for this strain-based model, which creep at lower strain rates for the same Σ [12]. This is also supported by the lower cavity densities predicted between harder, as compared with softer grains, at a given time. These results can be found elsewhere [69].
Failure of ex-service Type 316H at 550°C is of interest to this study, but bicrystal data at this temperature are not available at present. According to Huang [51] and Arnaud [87], the nucleation constant a ′ for Type 316/316H stainless steels is constant in the range 525-700°C. Hence, the intergranular failure response of a hypothetical bicrystal of ex-service Type 316H at 550°C, deforming by the CPFE model with parameters given in Table 1, is simulated by the Dyson-type model with a ′ = 8.9 × 10 6 mm −2 . The predictions of t f over the stress range 230-320 MPa at 550°C are shown in Figure 7. The predicted trend by the Dyson-type model for a bicrystal of ex-service Type 316H at 550°C in Figure 7 lies below the polycrystalline experimental creep strength data for this material. This difference between bi-and polycrystalline data appears identical to that at 700°C, which suggests that the predicted trend of the Dyson-type model at 550°C is reasonable. Since bicrystal creep strength data for Type 316H at 550°C are not available, we assume that the predicted response with a ′ = 8.9 × 10 6 mm −2 can be used in the following sections as reference data to calibrate the modified stress-based nucleation model at this temperature. In Section 4.1, the calibrated Dyson-type model is again used to simulate the polycrystalline failure response.

Classical nucleation theory with distribution of particle strengths
In this section, the modified classical nucleation theory of Equation (5), with a log-normal distribution of nucleated cavity/particle geometries with different volumetric shape factors f v (c ′ ) (or dihedral angles c ′ ) is coupled to the diffusive growth model. As before, we first consider failure of a bicrystal of Type 316 at 700°C. In this section, the limiting case of t red 1 from Section 2.3 is examined, i.e. no relaxation of stress along the particle-free boundary regions, i.e. T nuc = T n . Before evaluating the failure response for different distributions of f v (c ′ ), we demonstrate the predicted trend in creep strength when the original form of the vacancy condensation theory from Equation (2) is employed (N.B. all cavity/particle geometries are assumed the same at a given material point). In this case, the predicted trend in creep strength for Type 316 at 700°C with c ′ = 5°lies close to experimental data in Figure 8, but does not capture the variation of t f with macroscopic applied stress Σ. For 90 < Σ < 140 MPa, failure times are approx. an order of magnitude shorter compared to the data, whereas at Σ < 75 MPa there is virtually no failure. This illustrates the deficiency of Equation (2) discussed in section 1.3. Next, the response of the bicrystal is simulated by the modified vacancy condensation theory of Equation (5)   found that a distribution of log (m(f v )) = −9 and log (s(f v )) = 0.6 captures the trend in experimental creep strength data at 700°C, as shown in Figure 8(a), with the parameters from Table 1. Examining the variation of the mean cavity number density along the boundary, predicted by Equations (2) and (5), in Figure 8(b) helps explain the predictions in Figure 8(a).
The sharp drop in failure time for Σ > 90 MPa from Figure 8(a), when Equation (2) is employed, is not present when a distribution of cavity/particle geometries, or equivalently cavity nucleation thresholds, is considered. In Figure 8(b), the original nucleation model predicts virtually no nucleation for Σ < 92 MPa, whereas approx. 80-100% of all possible nucleation sites are depleted for Σ > 92 MPa in less than 1 h. This explains the sharp increase in t f from Figure 8(a) for c ′ = 5°. On the contrary, the variation of N h with time is smoother when a distribution of f v (c ′ ) with log (m(f v ),s(f v )) = (−9, 0.6) is assumed. For the assumed distribution of f v (c ′ ), approx. 80% of all possible nucleation sites are depleted at t < 1 h at 160 MPa, but a continuous increase in cavity number density is found between 10 −3 and 300 h for Σ < 160 MPa. The findings from Figure 8(b) demonstrate that the stress-based nucleation theory could predict a continuous nucleation process when cavity formation at different populations of nucleation sites at a given material point is considered. Now, assume that the distribution of f v (c ′ ) is a physical characteristic of the grain boundary in Type 316/316H bicrystals. Then, the calibrated modified nucleation model from Figure 8(a) with log (m(f v ),s(f v )) = (−9, 0.6) can be employed to simulate the failure response of ex-service Type 316H bicrystals at 550°C and compared to that predicted by the Dyson-model of Section 3.1, assumed as the reference bicrystal response in the range Σ = 230-320 MPa. Predictions are given in Figure 9(a). The results show that failure times at 550°C are shorter by two orders of magnitude for log (m(f v ),s(f v )) = (−9, 0.6). In Figure 9 (a), for the set of parameters for ex-service Type 316H from Table 1, a distribution of log (m(f v ),s(f v )) = (−5, 1.4) predicts the same trend in creep strength at 550°C as the reference response from Section 3.1. To further support this, the evolution of N h along the boundary predicted by the modified vacancy condensation model for log (m(f v ),s(f v )) = (−5, 1.4) at 550°C agrees with that given by the Dyson-model from Section 3.1 (see [69]).
The determined distribution of f v with log (m(f v )) = −5 has a mean c ′ of 15°, which gives a lenticular profile of the nucleated cavity/particle geometry, compared to c ′ = 6°for log (m(f v )) = −9. The higher dihedral angle is more realistic [26,53]. It is, therefore, assumed that log (m(f v ),s(f v )) = (−5,1.4) characterises the physical distribution of nucleated cavity/particle geometries along this bicrystal boundary. This assumption would require that the distribution withlog (m(f v )) = −5 also captures the trend in creep strength in Type 316 at 700°C. Employing log (m(f v ),s(f v )) = (−5, 1.4) at 700°C with the parameters from Table 1 overpredicts failure times by ∼2 orders of magnitude and a change in the adopted parameters is required to capture the trend in response at 700°C. From the physical parameters in Table 1, the surface energy (g s ) and grain-boundary diffusivity could be affected by the microstructural state of the boundary [88]. Note that g s has a greater effect on the nucleation rate in Equation (5). We modify the value of g s and find that the experimental trend in creep strength of Type 316 at 700°C in Figure 9(a) for log (m(f v ),s(f v )) = (−5, 1.4) is captured when g s is set to 0.8 J/m 2 . It is likely that g s is reduced by the presence of impurities along particle/matrix interfaces [89,90]. However, results from [26] suggest that such drastic reductions in g s (∼60%) due to temperature change or presence of impurities is unlikely. In addition, the dihedral angle c ′ depends on the energy balance between the boundary and particle/matrix interface [26]. This balance could potentially change if there are differences in particle morphology or dislocation features at the boundary at 550°C and 700°C. Hence, it could be possible that different distributions of f v (c ′ ) characterise the grain boundary across the temperature range. Detailed information about such differences in Type 316/316H is not available at present. These aspects reveal limitations of the mechanistic nucleation model, which are discussed in Section 5.
In the following sub-section, the distribution of log (m(f v ),s(f v )) = (−5, 1.4) is probed where the vacancy condensation process is driven by the local particle stress (Section 2.3.2).
3.3. Classical nucleation theory using maximum particle stress with a distribution of particle strengths The failure of bicrystals of ex-service Type 316H at 550°C is simulated for Σ = 250-320 MPa. The stress driving the nucleation process in Section 3.2 was assumed to be the average normal traction T n (i.e. t red 1). In this section, the limiting case of rapid relaxation of stress in the matrix region and increase in particle stress is examined (t red 0) as represented by Equation (19). The possible increase in stress at the particle due to the generation of prismatic dislocation loops as the material creeps is also quantified here. Three limiting cases for the coupling of normal and shear tractions (Equation (20)) along the grain boundary are examined: (i) a n = 1, a t = 1 (normal only), (ii) a n = 1, a t = 1 (shear only) and (iii) a n = 1, a t = 1 (simple coupling). The analysis aims to provide insights into the effect of this coupling on damage development and failure. This will be examined further in Section 4.2.
Firstly, the modified model of Equation (5) is coupled to growth and T nuc is given by Equation (19) for a n = 1, a t = 1. When the distribution of f v with log (m(f v ),s(f v )) = (−5, 1.4) is employed the resulting failure times are shorter by an order of magnitude, compared to the reference response of the Dysontype model predictions (Section 3.1). This reference trend in the bicrystal creep strength at 550°C is only obtained with (a n = 1, a t = 1) and log (m(f v ),s(f v )) = (−2.8, 1.3), see Figure 9(b). Increasing the mean of the distribution to a higher dihedral angle offsets the increase in particle stress by a factor (d p /l p ) −2 (see [69] for more detail). In practice, t red will be finite and the maximum stress at particles will steadily increase, giving rise to a gradual acceleration of the nucleation rate along the boundary. Employing Equation (7), t red for γ-Fe at 550°C is estimated here to be ∼100 s, which is short compared to the creep life of the samples. Hence, the limiting case considered in this section of t red 0 is explored further. Next, we consider the other limiting case of traction coupling -(a n = 1, a t = 1). When the crystallographic orientations either side of the boundary are different a shear stress is induced along the grain boundary [12]. This shear stress would always be much smaller than the applied Σ in the bicrystal and, therefore, use of a shear based criterion, i.e. (a n = 1, a t = 1), employing the same distribution of f v , gives a negligible nucleation rate. Also, we find that the failure response for the simple coupling (a n = 1, a t = 1) is identical to that for (a n = 1, a t = 1), since T n dominates the coupling. Of interest is also the evolution of the residual stress at particles due to the development of local dislocation structures, t GB p . We have tracked the value of t GB p in these simulations and values less than 10 MPa were obtained, suggesting that the local increase of stress is dominated by diffusional relaxation of stress local to the particles. This is discussed further in Section 5.
In this section, we have demonstrated that both the Dyson-type model and the modified classical nucleation theory can be calibrated to fit experimental data. In order to probe more fully the suitability of the resulting models it is important to examine the failure response of polycrystalline materials. In Section 4, the models are used to simulate the failure of polycrystalline EX Type 316H at 550°C.

Modelling failure of polycrystalline Type 316H
Here, intergranular creep failure in polycrystalline ex-service Type 316H stainless steel at 550°C is simulated. The array of 39 irregularly-shaped grains from Figure 5(b) is used. The CPFE model is calibrated to macroscopic short-term plasticity and creep data in the range 250-320 MPa in the absence of damage. The employed model parameters are given in Table 1. The response of 10 sets of 39 crystallographic grain orientations is examined. Only the results for one orientation set are presented for the sake of brevity since the trends in response are similar in all 10 sets. The aggregate is loaded in uniaxial tension at its top face and the boundary conditions from Figure 5(b) are adopted. Four meshes of different degrees of refinement with interface element lengths (i) l c = 2.6 μm, (ii) 1.8 μm, (iii) 1.0 μm and (iv) 0.75 μm were examined for the mesh convergence study. Convergence was assumed when differences between the macroscopic creep strain and t f at the end of the simulations were <5% with further refinement. The average normal stress along a given boundary was compared between meshes (i)-(iv) as a separate criterion for mesh convergence to ensure a difference of <15% with further refinement. We find that the mesh with l c = 1.8 μm was selected, with an example of the mesh convergence study shown in Figure 10. Based on the limited amount of grain-boundary sliding in Type 316 [12], sliding here is not considered (ḋ 0t 0) with the aim of isolating the effect of intergranular cavitation on creep deformation. The average grain size is ∼40 μm. The Dyson-type model and modified vacancy condensation model, as calibrated in Sections 3.1 and 3.3, are coupled to diffusive growth in analysing the 39-grain aggregate response. The ability of the nucleation models to capture (i) macroscopic creep deformation, (ii) variation of cavity density with orientation of grain boundaries to the load axis, and (iii) the time to failure is evaluated. Damage development with respect to crystallographic grain orientations, and the relation between nucleation and local stress state in the vicinity of boundaries are also examined.

Local Dyson-type nucleation model and diffusive growth
The polycrystalline response of ex-service Type 316H at 550°C is simulated in this section where cavitation is described by the Dyson-type model as calibrated against biaxial data in Section 3.1, with a ′ = 8.9 × 10 6 mm −2 , coupled to diffusive growth. The macroscopic creep curves in the range Σ = 280-320 MPa in the presence and absence of damage are compared to experimental data in Figure 10. The predicted failure times of t f = 42, 35 and 25 h at 280, 300 and 320 MPa for the polycrystal, where t f is taken as the end of the tertiary regime, are shorter by ∼22 times compared to the experimental response.
The shorter failure times predicted are the result of rapid nucleation along boundaries in the early stages of creep. Secondary creep regimes are relatively short and the tertiary creep regime occurs at an early stage. This is illustrated by the time variation of cavity number density (N h ) and damage fraction (f h /f c h ), averaged individually along all 115 grain boundaries and plotted as a function of boundary orientation to the load axis, u GB , in Figure 11(a,b). The polycrystalline average value of N h across all boundaries at 1 and 30 h for Σ = 300 MPa are 1.8 × 10 5 mm −2 and 7 × 10 5 mm −2 , respectively. These values are larger by 1-2 orders of magnitude as compared to the average N h values in the bicrystal (Section 3.1) at 1 h (6700 mm −2 ) and 30 h (2 × 10 5 mm −2 ) for Σ = 300 MPa.
The rapid nucleation, predicted by the calibrated Dyson-type model in the early stages of creep, allows cavity coalescence and failure along boundaries to occur earlier in the polycrystal, compared to experimental data. The reason is that the local inelastic strain rate in the vicinity of boundaries is faster in the polycrystalline aggregate, compared to the bicrystal for a given macroscopic applied stress. This phenomenon is due to localisation of strain to accommodate incompatibilities between anisotropic grains in the polycrystal. Another important finding from Figure 11 is that N h is not a function of u GB when the Dyson-type model is employed. This does not resemble the general trend often observed in FCC materials, e.g. [58][59][60], where higher cavity densities are found on boundaries with u GB ≈ 45-60°as compared to boundaries, transverse to the load axis. The model of Equation (1) is assumed to be dependent on the local shear strain rate in the vicinity of the boundary (ġ max ), which in turn depends on the crystallographic orientation of neighbouring grains. This explains the limited variation of N h with u GB in Figure 11(a). Here, higher cavity densities are found along boundaries enclosed by plastically-softer grains with crystallographic orientations closer to the [001]-direction on the IPF (see [69] for more details). Boundaries, transverse to the load axis, tend to fail (i.e.f h /f c h = 1) as illustrated in Figure 11(b). These boundaries carry higher normal average tractions, as illustrated in Figure 11(c), and diffusive growth on these boundaries is dominant, leading to earlier cavity coalescence. This prediction is in agreement with the observed distribution of failed facets from [24,[91][92][93].
The polycrystalline failure times predicted by the Dyson-type model are shorter than those in bicrystals from Section 3.1 for the same value of a ′ . This disagrees with the trend in creep strength data for Type 316 bi-and  (1)) has a number of deficiencies and it does not provide further understanding of the microscopic nucleation process in polycrystalline materials. The focus of the remaining sections is on the modified classical nucleation theory.

4.2.
Classical nucleation theorya distribution of particle strengths, maximum stress at particles and diffusive growth 4.2.1. Case 1maximum particle stress for limiting cases a n = 1, a t = 1 or a n = 1, a t = 1 The failure response of the 39-grain aggregate is examined here by coupling diffusive growth to the modified classical nucleation theory with T nuc given by Equations (19) and (20). The cases (a n = 1, a t = 1) and (a n = 1, a t = 1) are evaluated first with the distribution log (m(f v ),s(f v )) = ( − 2.8, 1.3), determined in Section 3.3. The macroscopic creep response for Σ = 280-320 MPa in the presence and absence of damage is compared.
We find that for (a n = 1, a t = 1) the predicted failure times of 30, 38 and 52 h at Σ = 320, 300 and 280 MPa are ∼20 times shorter than the experimental response; the creep curves are shown elsewhere [69]. The reason is that local concentrations of stress arise along grain boundaries in the polycrystal. These stress concentrations increase the tractions carried by individual boundaries, and respectively the nucleation rate for the same parameters, compared to the bicrystal's isolated boundary for the same Σ. At 300 MPa and a n , a t = (1, 1), the average value of t * = T n /(d p /l p ) 2 ∼900 MPa along the bicrystal boundary. Along boundaries almost normal to the loading direction in the polycrystal, i.e. u GB = +75 90 • , maxima of up to ∼1600 MPa (i.e. T n 520MPa) are found. The nucleation model of Equation (5) is a function of local stress at particles. As a result, the average cavity number density along individual boundaries in Figure 12(b) follows closely the distribution of t * with u GB . The profile of N h , therefore, shows higher cavity densities along boundaries normal to the applied load (full set of results presented in [69]). Diffusive growth is again dominant along transverse boundaries enclosed by harder grains. The observations reveal that the distribution of N h with respect to u GB , is not consistent with the experimental trend of higher cavity densities on inclined boundaries in FCC materials [58][59][60]. The polycrystalline response is also simulated by adopting Equations (19) and (20) for the case (a n = 1, a t = 1) when the maximum particle stress is a function of shear tractions only. For log (m(f v ),s(f v )) = ( − 2.8, 1.3), we found that the boundary shear tractions are not sufficiently high to generate cavity densities which could lead to cavity coalescence and affect the global failure response.
To illustrate an important point, we could assume that the vacancy condensation process is indeed driven by the macroscopic shear tractions only (which can give rise to tensile stresses normal to the particle/matrix interface). A suitable distribution of f v (c ′ ) with a lower nucleation barrier (i.e. lower m( f v )) for the assumed lower T nuc (a n = 1, a t = 1) can be determined to match experimental creep strength data. We find that a distribution of log (m(f v ),s(f v )) = (−5.0, 1.4) approximates the macroscopic creep experimental data; N.B. creep curves are shown elsewhere [69]. For this scenario, examination of predicted distributions of averaged N h and f h /f c h along individual boundaries revealed that higher cavity densities develop along boundaries inclined at ∼45-60°to the applied load. However, virtually no cavities nucleate or coalesce on transverse boundaries (N h < 10,000 mm −2 at u gb = +75 90 • ). The profile of N h for (a n = 1, a t = 1) is qualitatively similar to that in FCC materials [58][59][60], but boundaries normal to the applied load remain intact. This disagrees with observations from experiments on different grades of Type 316, e.g. [24,91,93].
These results illustrate that the limiting cases of the coupling between (a n , a t ) are unable to capture the damage patterns in the material. The following section determines a suitable form of the coupling. 4.2.2. Case 2maximum particle stress with suitable couplinga n = 1, a t = 1 The failure response of the polycrystal is now examined using the modified vacancy condensation model for the case (a n , a t ) = 1. To calibrate the coupling between the tractions, the following assumptions are made. Firstly, the coupling should predict the highest cavity number densities along inclined boundaries (u GB = +45 60 • ) with lower cavity number densities along boundaries transverse to the applied load (u GB = +75 90 • ). Note that Figure 12. Predicted creep curves in the absence and presence of damage, EX Type 316H at 550°C, modified classical nucleation theory with a n , a t = (1.5, 0.3) and log (m,s) = −2.0, 1.3. Directions of resultant traction at particles (κ) and of maximum principal stress (β I ) to the boundary plane of selected facets for different combinations of a n , a t . Boundaries with N h >10 5 mm −2 at t/t f = 0.75 were examined. N.B. k = tan −1 ((T n /a n )/(T t /a t )) whereas β I is averaged as extracted from a row of continuum elements, either side of a boundary. detailed experimental data for the actual distribution of N h along boundaries with u GB are not available; the assumed distribution is informed qualitatively by results from 316H and other FCC materials [58][59][60]. As identified in the preceding section, modification of (a n , a t ) also requires that the assumed distribution of f v (c ′ ) is changed in order to capture the trend in macroscopic creep failure. Furthermore, observations reported in [1,93,94] suggest that cavities nucleate at GB particle interfaces, coinciding with the direction of the local maximum principal stress. Hence, we assume that the local maximum principal stress and the effective traction driving nucleation as defined in Equation (20) must coincide.
We find that a n , a t = (1.5, 0.3) and log (m(f v ),s(f v )) = (−2.0, 1.3) satisfy the assumptions above and capture the general trend in macroscopic creep response of Type 316H at 550°C in the range 280-320 MPa, as shown in Figure 12(a). A difference in failure time within a factor of two is acceptable for the purpose of this study, due to the scatter in experimental failure data for Type 316H [95]. Figure 12(b) reveals the direction of the resultant traction at particles (κ) nearly coincides with the orientation of the local maximum principal stress in the lattice with respect to the boundary (β I ).
The predicted distribution of averaged N h and f h /f c h along individual boundaries is examined in Figure 13(a,b) at 300 MPa. At 460 h, average cavity number densities of ∼4×10 5 mm −2 are found along boundaries with u GB = +45 60 • . These values are almost double the cavity number densities of ∼2.1×10 5 mm −2 along boundaries with u GB = +75 90 • . The results in Figure 13(b) show that f h /f c h is higher along facets with u GB = +75 90 • . The predicted cavity density profile follows closely the variation of maximum particle stress t * with u GB for the set a n , a t = (1.5, 0.3) from Figure 13(c). Figure 13(d) shows the evolution of the averaged residual stress due to the generation of prismatic loops, t GB p , which was tracked during the simulations. The magnitude of t GB p is small compared to t * and justifies ignoring this potential contribution to the local stress here. This finding is in agreement with observations from Section 3.3. The contribution of terms of this type is discussed further in Section 5. Model predictions demonstrate that higher cavity and intergranular facet crack densities are found between grains with harder crystallographic orientations, away from the [001]-direction on the IPF in Figure 14. This is expected since these grains will carry higher relative stresses, causing more rapid nucleation and growth rates according to the stress-based models considered in this section. Cavity coalescence by diffusive growth is dominant on facets normal to the load direction, e.g. boundaries 18, 83 in Figure 14.
The evolution of average damage fraction f h /f c h , cavity number density N h and normal boundary traction T n with time along selected grain boundaries is examined in Figure 15 at 300 MPa for a n , a t = (1.5, 0.3). In Figure 15(b), N h along boundaries nearly normal to the load direction (boundaries 7,18,83 with u GB = +75 90 • ) does not increase rapidly at t > 100 h and cavity coalescence is driven by diffusive growth due to the high normal tractions on these boundaries (Figure 15(c)). Along boundaries 23 and 43 (u GB = +45 60 • ), the normal traction is lower, but the nucleation rate is higher. The development of damage along boundaries 23 and 43 in Figure 15 (a) is correlated with the evolution and higher values of N h along these boundaries, as shown in Figure 15(b). This allows coalescence at the lower values of T n on these inclined boundaries. Hence, some inclined facets, e.g. boundaries 23 and 43, where a sufficiently high number densities of cavities have accumulated, also fail and lead to interlinkage of inclined and normal facets. This is an important observation since the facet interlinkage process is a precursor for microcrack formation in Type 316H [96].
Note that at failure, the stress carried by the failed boundaries reduces to zero. Figure 15(c) shows how the tractions from failed boundaries 7, 18, 23 and 49 are redistributed to boundaries in the polycrystal, which are intact, e.g. boundaries 14, 17, 41 and 44. In turn, nucleation and diffusive growth accelerates on the latter boundaries, as seen on Figure 15(a,b). The progressive failure and stress redistribution process causes an increase in stress in non- damaged regions, further illustrated in Figure 16(a-c), where contours representing the axial stress in the grains is plotted at different times. This increase in stress causes faster creep rates in non-damaged regions of the aggregate. We find that the accelerated creep rate in the grains due to this stress redistribution process constitutes ∼96% of the increase in macroscopic creep rate during tertiary creep.
Note that the distribution of f v (c ′ ) and set of (a n , a t ) have been arbitrarily selected in this section. It is important to examine the creep strength of the bicrystal from Section 3.3 with the set a n = 1.5, a t = 0.3; log (m(f v ),s(f v )) = ( − 2.0, 1.3). With these parameters, the predicted failure times for the ex-service Type 316H bicrystal are 1070, 683 and 423 h at 280, 300 and 320 MPa, respectively. This is illustrated in Figure 17(a). These failure times are longer by a factor of 1.5, compared to failure times for the polycrystal, and by a factor of 5, compared to the reference bicrystal t f − S response at 550°C (Section 3.1). А distribution of f v (c ′ ) with  log (m(f v ),s(f v )) = ( − 3.6, 1.3) is needed to capture the reference creep strength response from Section 3.1 of the bicrystal for the assumed coupling of a n , a t = (1.5, 0.3) [69]. This indicates that physical phenomena are possibly missing from the assumed model of Equation (5), since the failure response of bi-and polycrystalline samples is not captured simultaneously with the assumed sets of f v (c ′ ) and (a n , a t ). Additional comparisons to experimental data of the creep deformation and failure response of ex-service 316H are provided in Figures 17-19 to explore the ability of the model framework and associated assumptions to predict the multi-scale response. The predicted minimum creep rate and uniaxial creep ductility (failure strain) of the material agrees with the magnitude and trends in experimental data at 550°C in Figure  17(b-c). Furthermore, the predicted intergranular creep failure morphology from Figure 14(c) at Σ = 300 MPa is compared in Figure 18 to various experimentally observed failure morphologies in 316H at 550°C and comparable stress states. Findings reveal a qualitative agreement with the development of  creep microcracks predominantly along boundaries transverse to the applied stress which begin to progressively interlink with those on inclined boundaries.
In Figure 19, the average cavity spacing from the model predictions is compared to several SANS (small-angle neutron scattering) measurements of onaverage cavity spacings in as-received and aged pre-strained 316H specimen, crept at 550°C and 300-320 MPa. Note that the experimental results in Figure 19 do not resemble fully the microstructural condition of the exservice conditions examined here and are extracted from different heats for 316H Refs. [9,100,102]. Nevertheless, the comparison to data from as-received, non-pre-strained material suggests that the model and corresponding assumptions from Sections 2.2 to 2.3 overestimate the on-average spacing of nucleated Figure 18. (a) Simulated damage morphology at 550°C, 300 MPa. (b) SEM images of microcracking in as-received 316H at 550°C∼300 MPa near failure region in crept specimen (after [100], used with permission), (c) backscatter image of localised intergranular failure of exservice pre-strained 316H at 550°C, 320 MPa (after [9], used with permission), and (d) creep microcracking in ex-service 316H at ∼550C crept at 320 MPa up to 6.8% creep strain (after [101], used with permission). Applied maximum principal stress direction is indicated where given. Figure 19. Comparison of average cavity spacing with timemodel predictions for the polycrystal at 300 MPa at compared to data from Ref. [9,100]. The assumed limit for nucleation sites based on the l p is illustrated. cavities by a factor of 3-10. This further illustrates potential missing features in the nucleation model which would account for nucleation activity at additional sites. The cavity spacing data for 8% pre-strained and crept material shows even larger number density of cavities compared to non-pre-strained data which strongly suggests that these additional nucleation sites are associated with inelastic deformation at or near the grain boundaries. This observation is also aligned with recent findings in FCC copper in Ref. [103]. Note that incorporation of additional nucleation sites in the models of Sections 2.2-2.3 would accelerate the predicted failure times in Figure 17(a). This would suggest that the cavity diffusive growth model assumed in Section 2.4 must also consider additional phenomena in order to capture the macroscopic failure response correctly. This is discussed further in Section 5.

Suitability of examined cavity nucleation models
The Dyson-type model of Equation (1) was calibrated to bicrystal Type 316 stainless steel creep strength data at 700°C, giving a ′ = 8.9 × 10 6 mm −2 for the nucleation constant. The predicted macroscopic times to failure for the polycrystalline material using this value were shorter by a factor of 22, compared to experimental data. Fitting the model to polycrystalline material data gives values of a ′ closer to experimentally-recorded values of 6 × 10 5 mm −2 [87], and 2 × 10 5 mm −2 estimated from the mechanistic 1-D nucleation model of He and Sandstrom [28]. This type of model, however, predicts a distribution of cavity number density (N h ) that is independent of the orientation of a grain boundary with respect to the load axis. This is not consistent with the general experimental trend observed in FCC materials, which develop the highest cavity densities on boundaries inclined between about 45-60°with respect to the loading direction [58]. We, therefore, focus our attention on the mechanistic vacancy condensation theory.
The modified classical nucleation theory provides insights into the sensitivity of cavity nucleation to local stress fields and type of nucleation sites. A unique set of a n , a t = (1.5, 0.3) for the coupling between normal and shear tractions at intergranular particles, contributing to T nuc , was found to predict a distribution of N h along grain boundaries in the polycrystal that resembles the experimentally observed trend in FCC materials, with the highest values of N h along facets inclined at 45-60°to the applied load, and cavity coalescence dominating along transverse facets. An important finding is that the orientation of the resultant traction with respect to the boundary for a n , a t = (1.5, 0.3) is coincident with the direction of maximum principal stress in the vicinity of the boundary. This suggests that the maximum principal stress near the boundaries indeed determines the coupling of local tractions at nucleation sites and controls the nucleation, which is in qualitative agreement with observations from [94]. Predicted boundary tractions and local particle stresses are larger than the global applied stress. Higher stresses along boundaries and at particles on the boundaries allow distributions of nucleated cavities with mean dihedral angles of 15-30°to be employed. The predicted nucleated cavity profiles with c ′ . 15 • are closer to the lenticular profiles of newly-nucleated cavities hypothesised in [104]. Experimental observations of cavity geometries in the early stages of nucleation in Type 316 are needed to confirm this hypothesis. Different distributions of log (m(f v ),s(f v )) are, however, required to fit the polycrystal and bicrystal data. Furthermore, the model prediction for average cavity spacing in Figure 19 appeared to overestimate the spacing as compared to experimental data for as-received 316H crept under similar stresses. Qualitatively, faster nucleation rates and more nucleation sites would be needed to approximate these data. This would also require slower rates of diffusive cavity growth if the model predictions were to agree with macroscopic failure data since higher cavity densities would result in larger values of f h at a given time. Another key finding is that extensive experimental data for the evolution of cavity densities with time under creep are sensitive to the cast and heat treatment of the material Ref. [9,100]. Such data are also scarce which makes calibration or validation of mechanistic nucleation models challenging. These observations suggest that although the model captures some important features of the material response, important physical features are missing in the model. This is discussed further in Section 5.2.

Missing features in mechanistic nucleation model and the need for further targeted experiments
Individual grain boundaries have different misorientations, which is a function of the crystallographic orientation of the grains on either side of the boundary. The boundary energy, boundary diffusivity [88] and particle shapes [64] will be different for boundaries of different misorientations. The resulting variation of properties and particle morphology (which would give rise to different distributions of f v (c ′ ), where c ′ is a function of the local grain-boundary and particle/matrix interface energies [30]) will vary with boundary misorientation. This will affect the nucleation response on individual boundaries. This is confirmed by Slater et al. [60] and Pommier et al. [105] who found that boundaries with misorientations of 25-55°are most prone to cavity formation in Type 316H/-L. The presence of impurities (e.g. phosphorus) or other phases (e.g. σ) will further modify the boundary and surface energies, and grain-boundary diffusivity [26,51]. These local features are not considered in the modified classical nucleation model. Note that the effect of boundary misorientation could be more pronounced in bicrystals, since macroscopic failure in the bicrystal is dominated by the properties of a single boundary, while failure of the polycrystal is determined by the development of damage on a large number of boundaries of different misorientation and orientation with respect to the load axis. Interpreting the distribution of f v (c ′ ) as a function of boundary misorientation could explain the difference in log (m(f v ),s(f v )) for bi-and polycrystalline samples. Information about the misorientation of the boundary in bicrystal experiments modelled here is not available and apart from the studies in [60,96,105] we are not aware of any systematic studies of the effect of local crystallography on damage development. The results from these studies are insufficient to inform the extension of the model in Equation (5) to account for GB misorientations. Advanced characterisation techniques, such as EBSD and correlative tomography, can provide detailed information at the sub-grain scale with the opportunity to conduct a range of targeted interrupted experiments on bicrystal samples with different degrees of misorientation for the grain boundary and different orientations with respect to the loading direction, and on randomly orientated polycrystals. Coupling these experiments with CPFE models of the type described here will provide additional insights, aiding the calibration of log (m(f v ),s(f v )) and how this changes with misorientation of the boundary.
Within the modified vacancy condensation model we have made a number of simplifying assumptions, in particular that the traction acting on a boundary is completely relaxed between particles early in life. There are two possible mechanisms for concentrating stress on the particles in this way: diffusional relaxation and the generation of dislocation structures around the particles, due to inelastic lattice deformation. We have assumed that the first of these is rapid and we have simply monitored the potential contribution from the generation of a representative class of dislocation features using a model that is consistent with that for inelastic deformation of the grains. The stress generated by these loops is found to be small in comparison to the applied stress. It is possible that diffusional relaxation occurs over a longer length scale than assumed here, and leads to a slower relaxation of stress between the particles (and build-up of stress on the particles). Also, other dislocation mechanisms specific to the grain-boundary could generate high stresses rapidly, such as dislocation pile-ups, or grain-rotations arising from the constraint of the boundary. Each of these effects require further investigation experimentally. High resolution EBSD in particular can provide information about inelastic deformation within the grains and how this varies from grain to grain. As well as validating CPFE model predictions, observations of this type, when combined with crystallographic information, can provide important information about the factors that mainly control the nucleation of cavities. The assumed generalisation of the Dyson model employed here does not capture the microscopic information concerning the distribution of cavities within a material observed experimentally. However, there is compelling evidence that inelastic deformation plays an important role in the cavitation process, e.g. recent work by Das et al. [103] on measuring the evolution of cavity densities in pure Cu supports this. The authors in [103] discuss that the cavity spacing was well-aligned with grain-boundary ledges/dislocations estimated from the 1D model of He and Sandstrom [28]. The findings from Figure 19 here suggest that assuming only second-phase grain boundary particles as nucleation sites may not give representative cavity densities in 316H. The missing nucleation sites, not considered here, are likely associated with the dislocation features at or near the grain boundary discussed by Ref. [28,103]. Higher cavity densities are also noted for pre-strained 316H by Jazaeri et al. [9] which further suggests that accumulation and impediment of dislocations at grain interfaces would give rise to more nucleation sites. As well as enhancing the stress at grain-boundary precipitates and other boundary features, dislocations can be the additional source of changes in energy that drive the nucleation process as discussed by He and Sandstrom [28] and Lim [106]. Another aspect of the cavitation process, which is not considered here, is sintering [107]. Compressive normal tractions along boundaries could result in cavity sintering, which will alter some of the predicted cavity distributions in the aggregate in Section 4. Grain-boundary sliding is also not considered in this study. In [12], it was found that the stress along transverse boundaries increases by <3% in Type 316 for the limited amount of sliding in this material at stress levels comparable to those considered here. This suggests that cavitation along transverse facets, predicted in Section 4, will be accelerated and the a n − a t coupling will be altered when grain boundaries slide relative to each other. The additional features, discussed above, which could be incorporated to the mechanistic nucleation model in Sections 2.2-2.3 suggest that the nucleation rates predicted in Section 4 will increase. This in turn would result in faster failure times in Section 4, with everything else kept constant. It is worth mentioning that two other aspects of the framework would require more careful considerations to account for this: . The diffusive growth model from Section 2.4 does not consider any reactions which could impede the diffusive transport along an imperfect boundary during cavity growth. In particular, interface reactions at particles-boundary interfaces [65,108] would affect the local stress at particles but also alter the rate of transport of matter away from cavities. . The simulated polycrystalline geometry in this study is a simplified quasi-3D aggregate. This was assumed as an idealisation of realistic 3D geometries in order to illustrate essential features in the response. Full 3D polycrystalline aggregates of realistic grain shapes, such as those studies in e.g. [48,79], should employ the mechanistic models presented here. The evolution of damage morphology in such cases will be more representative which could provide higher-fidelity estimates of failure times, compared to those reported here.
This discussion has identified a number of features of the cavity nucleation process which need to be considered further in the mechanistic nucleation models considered in this study. Further development of the models must be guided by experiments that provide detailed information at the sub-grain scale. Capturing these features could provide insights into the dependence of nucleation rate on inelastic strain rate, local stress state and the physical processes of cavity formation. The modelling framework developed here can serve as a tool to incorporate these features and explore their effect on the general characteristics at the meso-and macro-scales, discussed in this study.

Conclusions
. Cavity growth models have been combined with modified strain-and stressbased nucleation models within a CPFE-interface element modelling framework to simulate intergranular creep failure of bi-and polycrystalline samples of Type 316 stainless steel over the range 550-700°C. . The local strain-based cavity nucleation model, which is a function of the resolved shear strain rate, does not capture the distribution of higher cavity densities on boundaries inclined to the applied load observed experimentally, nor does it offer insight into the mechanistic nucleation process. . The classical nucleation theory, modified to consider a log-normal distribution of nucleated cavity/particle geometries (distribution of f v (c ′ )), captures the trend in bicrystal creep strength of the material and predicted continuous nucleation along an isolated boundary. . When it is assumed that the stress on a boundary between intergranular particles is completely relaxed, the stress on a particle achieves the maximum possible value. If this is assumed to drive the nucleation process, the newly-nucleated cavities assume lenticular shapes, closer to the hypothesised cavity geometry during the early-stage growth. . The coupling between local tractions at particles obtained by fitting experimental data, predicts the resultant traction direction to be coincident with the direction of local maximum principal stress with respect to the boundary plane. With this model the highest cavity number densities are predicted along inclined boundaries, while growth and coalescence dominate on boundaries perpendicular to the applied load in the polycrystal, agreeing with experimental observations in Type 316 stainless steel [24,60,93]. . Missing features have been identified in the mechanistic modified classical nucleation theory employed here. Further work is required to incorporate the effects of boundary misorientations, impurities, and faceted cavity/particle shapes on the nucleation process. The development of stress at the intergranular particles needs to be accurately modelled, and combined with the effect of any additional energy provided by the local dislocation structure.
. Development of more robust nucleation models requires experimental studies of bi-and polycrystalline samples, focused on early-stage cavity/particle morphologies, dislocation features near boundaries, cavity number density as a function of boundary misorientations and orientation to the load axis.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
where the primed quantities are associated with the rear face of the triple line element.