Modelling the growth and filling of creep-induced grain-boundary cavities in self-healing alloys

A set of numerical and analytical models is presented to predict the growth and contraction of grain-boundary creep cavities in binary self-healing alloys. In such alloys, the healing is realised by preferential precipitation of supersaturated solutes at the free surface of the cavity. The cavity grows due to the diffusional flux of vacancies towards the cavity, which is driven by the stress gradient along the grain boundary. Upon deposition of healing solute atoms on the cavity wall, effectively vacancies are removed from the cavity due to the inverse Kirkendall effect. The competition between the inward and outward vacancy fluxes results in a time-dependent filling ratio (i.e. the fraction of the vacancies removed from the original cavity) of the creep cavity. It is found that for stress levels lower than a critical stress σcr, the filling ratio can proceed to unity, i.e. to complete filling and annihilation of the pore. For applied stresses higher than σcr, complete filling is not achieved and the open volume of the creep cavity will continue to grow once a maximum filling ratio is reached at the critical time tcr. The critical stress σcr, critical time tcr, and time for complete filling th (if fully filling is achievable) are derived from the models for different combinations of parameters. The results from the analytical model and from previous nanotomography experiments are compared and are found to be in good agreement.


GRAPHICAL ABSTRACT Introduction
Creep takes place when metals are exposed to a combination of a constant mechanical stress below the yield stress and a relatively high temperature (higher than 0.4T m , where T m is the melting temperature). During creep loading of polycrystalline metals, quasi-spherical micron-sized cavities form preferentially at those grain boundaries oriented perpendicular to the load direction. Upon prolonged loading, such pores grow and coalesce into microand subsequent macrocracks, which leads to complete fracture of the sample or even the installation [1]. Traditionally, the strategy to improve the creep resistance and the high-temperature properties has been to create a microstructure which retards the cavity formation and minimises the cavity growth as long as possible [2][3][4][5]. However, irrespective of the metallurgical strategy used, cavities once formed will continue to grow, coalesce and the largest ones will finally induce catastrophic failure. Some years ago, the concept of self-healing has been proposed as an alternative strategy [6][7][8]. In this approach, the occurrence of creep cavities triggers dedicated solute atoms (so-called self-healing solutes) to move towards these cavities, to fill them and make them harmless before they grow to catastrophic dimensions. Provided that the diffusivity of the selected solute atoms is faster than the diffusivity of the host atoms, then due to the Kirkendall effect, the diffusion of the healing agent towards the cavity surface results in a net diffusion of vacancies in the opposite direction, i.e. the vacancies are removed from the creep-induced cavity [9] and the empty volume of the cavity becomes smaller. Furthermore, once the cavity is completely filled, the driving force for directed vacancy flux is removed and the damage can be considered to have been healed. As a result, the coalescence of neighbouring creep cavities is prevented. In case the growth of the cavities is faster than the rate of filling, then only partial cavity filling will take place and cavities will continue to grow but at a lower rate than in non-self-healing system. This is expected to have a much smaller effect on the desired life time extension.
The self-healing concept based on selective precipitation has been tested and verified for multiple binary and more complex metallic systems. Laha and co-workers [10][11][12][13] reported the unintended experimental finding that the addition of boron and nitrogen is effective in suppressing the cavity growth rate of austenitic stainless steels (with an fcc lattice structure). During creep, BN preferentially precipitates on the free cavity surfaces and thereby partially heals the cavities, leading to an increased creep rupture strength. Lumley and co-workers [14,15] demonstrated that underaged Al alloys showed a lower strain rate compared to a fully aged counterparts, which is contributed by the dynamic precipitation of free solute atoms and the subsequent retardation of dislocation motion during the creep test. In order to demonstrate intentional filling of creep damage by tailored precipitation on the cavity wall, Van Dijk and Van der Zwaag initiated a number of creep damage studies using selected binary Fe-X J Mater Sci (2022) 57:12034-12054 and ternary Fe-X-Y alloys, including Fe-Cu [16][17][18], Fe-Au [19][20][21][22][23], Fe-Mo [24], Fe-W [25], and Fe-Au-W [26]. These systems were selected on the basis of the following considerations: (1) they can be brought into a supersaturated state (typically 1 at.%) at a creep relevant temperature at which the alloy is in a ferritic state (here 550°C) by prior homogenisation and rapid quenching, (2) the energy barrier for precipitation on a free surface is considerably lower than that in the grain interior, (3) the diffusion coefficient of the solute is higher than that of the iron atoms, and finally (4) the solute atoms are not consumed in other reactions while migrating to the cavity [8,22]. These alloys serve as model alloys for future self-healing multicomponent ferritic creep steels as they have the advantage that no other metallurgical reactions than the intended healing reaction can take place.
Multiple characterisation techniques have been applied to monitor the damage formation and the subsequent precipitation healing in these alloys, including electron microscopy techniques (SEM, TEM, EBSD), atom probe tomography (APT), and advanced synchrotron X-ray nanotomography. In an early research by Zhang and co-workers [22], it is demonstrated that a healing efficiency of 80% can be achieved for low stress levels, while for high stress levels, the healing precipitation rate cannot catch up with the damage rate and that leads to a lower healing ratio. Recent studies by Fang and co-workers [19,25] revealed detailed 3D distribution of creep cavities and the healing precipitates in Fe-Au and Fe-W systems with a voxel size down to 30 nm. It is found that the creep cavities show a strong preference for the stress-affected grain boundaries and that the precipitates nucleate at one or multiple sites at the surface of the cavities. The early-stage cavities are usually found to be isolated and small in size (tens of nanometres in diameter) with spherical or equiaxed shapes, while at a late stage of creep, the grainboundary cavities are likely to link with their neighbours, resulting in more complexed shapes and larger dimensions (equivalent diameter larger than 10 lm). It is also demonstrated that isolated cavities can be filled completely, while the filling ratio of the linked cavities first increases, but then decreases after the linkage takes place. It is therefore crucial to heal the cavities before the linkage occurs. In order to predict the pore filling behaviour of the creep cavities in such binary alloys, a semi-quantitative model has been proposed by Versteylen and co-workers [9]. In this model, the growth or shrinkage of the cavity is determined by the competition between the inward and outward vacancy fluxes. The inward flux, which is controlled by the stress-induced gradient in the chemical potential on the grain boundary surrounding the creep cavity, diffuses towards the cavity and therefore leads to further growth of the cavity. The opposite flux, resulting from the transport of the supersaturated solute into the cavity, removes vacancies from the cavity by the Kirkendall effect and results in a shrinkage of the cavity. If the outward vacancy flux is higher than or equal to the stressdriven inward flux, the cavity can be fully healed, and a critical stress can be calculated accordingly. A full healing can be achieved when the applied stress is lower that the critical value. In the case where the stress-driven inward exceeds the Kirkendall flux, the cavity cannot be fully healed, but a desirable reduced strain rate can be expected. The Versteylen model [9] was used to evaluate the healing efficiency of the iron based binary systems, and it was concluded that Au is the most healing-efficient solute. However, in the model by Versteylen and co-workers the normal stress at the grain boundary was not calculated explicitly and the solute diffusion was simplified as a stationary time-independent flux. In reality, however, with the evolution of the solute diffusion profile in the matrix, the solute flux decreases with time [22]. A critical time, where the stress-induced inward vacancy flux exceeds the healing flux and the healing cannot keep up with the cavity growth, is expected. Therefore, the time dependence should be taken into consideration for a more comprehensive understanding of the self-healing process.
The present work aims to investigate the time dependence of the self-healing behaviour and predict the time evolution of the cavity under different conditions. We propose a model to predict the evolution of the open volume of a grain-boundary cavity as a function of time. The inward vacancy flux is driven by the stress distribution along the grain boundary, while the time-dependent outward vacancy flux results from the solute segregation on the cavity surface. Both the stress distribution and the solute flux are calculated with a multi-physics finite element package (COMSOLÒ [27]). Based on the accumulated inward and outward vacancy fluxes, the model can be used to estimate the evolution of the filling ratio, the critical stress below which a fully filling can be achieved, the critical time after which the cavity growth exceeds the healing, and the time needed for a complete filling. The influences of the grainboundary-to-bulk diffusivity ratio, the level of solute supersaturation, and the spacing between neighbouring cavities on the healing efficiency are addressed explicitly and solutions are proposed as a function of the imposed stress. The model does not address the initial nucleation of the precipitate on the cavity wall nor does it take into account the kinetics of the internal atomic transport within the precipitate formed. The paper ends with a simple analytical model for the pore filling process taking the results from the numerical model to justify the assumptions made in the analytical model. The analytical model provides a quick estimation of the critical stress, the critical time, and the healing time.

Model geometry
As shown in Fig. 1a, the rotationally symmetric simulation box is a cylinder with radius k and height H. The radius k reflects half the distance between the centres of neighbouring creep cavities on the same creep affected grain boundary and the height H reflects half of the grain size. The vertical axis (z) and the radial axis (r) are indicated by the red dashed lines. The z axis is set as the symmetry axis, while the horizontal plane with z = 0 is set as a symmetry plane. It is assumed that a pre-existing cavity is located at the centre. According to Raj and Ashby [28], creep cavities are lens-shaped with a cavity radius a, an equilibrium opening angle of w & 75°, and a cavity height h that scales with the cavity radius as h/a = (1 -cos(w))/sin(w) & 0.77 [29]. In the present model, the cavity is assumed to have an ellipsoidal shape to simplify the calculation of the stress distribution (see 2.2 for details). The cavity radius is a, with a corresponding cavity height of h & 0.77a and an ellipsoidal cavity volume of Note that for a lens-shaped cavity the cavity volume is slightly smaller: In the experimental tomography studies on the Fe-Au system [19], the evolution of the average volume of the cavities is fitted to a power law with the form ofV ¼ kt n . Taking k = 0.33 lm 3 h -0.8 and n = 0.8 (corresponding to a stress level of 100 MPa), the average radius of the creep cavities with an ellipsoidal shape yields 0.4 and 0.5 lm after 10 and 20 h creep, respectively. The cavity radius a is taken as 0.5 lm in the present research. In all simulations, the grain-boundary width is assumed to be d = 2 nm and the height of the simulation box is fixed at H = 10 lm. The radius of the simulation box k varies from 2.5 to 25 lm. During the healing process, the solute in the matrix will be depleted by diffusion towards the cavity. The vertical diffusion length (parallel to z direction) can be calculated by , where D s m is the solute diffusivity in the matrix. Given D s m = 7.39 9 10 -19 m 2 /s [30], this combination of pore and matrix dimensions ensures that the reservoir of solute atoms (i.e. the number of supersaturated solute atoms in the total matrix volume considered) is not depleted during the healing process.

Cavity growth by stress-driven vacancy diffusion
Multiple theoretical models have been proposed to describe the growth of creep cavities. Depending on the applied stress, temperature, and the creep stage, the dominant cavity growth mechanism can be diffusion, plasticity, grain-boundary sliding, or a combination thereof [31]. A coupled model was proposed [32][33][34][35], which proposed that the creep cavities, their growth being controlled by grain-boundary diffusion of vacancies, are embedded in a matrix controlled by power-law deformation. Generally, diffusive growth dominates cavity behaviour in case of a small cavity, a low stress level r 0, and a low temperature T, while plasticity effects dominate otherwise. A diffusion kT _ e 1=3 [36] has been introduced to estimate the conditions in which diffusive growth dominates, where _ e is the strain rate, k is the Boltzmann's constant, D gb is the grain-boundary diffusivity, and X is the atom volume. As (a/K) increases, the creep flow becomes more important, and the power-law limit is reached when (a/K) ? ?. Taking the creep data of the Fe-Au, Fe-W, and Fe-Au-W model alloys [19,22,25,26], at a stress level of 300 MPa, the diffusion lengths K in these systems correspond to 17,19, and 26 lm, respectively. These diffusion lengths are comparable to or larger than the size of the simulation box. Therefore, it is reasonable J Mater Sci (2022) 57:12034-12054 to assume that the growth of the cavity is properly captured by the vacancy diffusion through the grain boundary.
In the model, it is assumed that the cavity growth at elevated temperatures results from vacancy diffusion and that the vacancy diffusion is driven by the stress gradient. The vacancy flux through the grain boundary in contact with the cavity can be written as [37][38][39][40] in which D v gb is the vacancy diffusivity at the grain boundary, x v is the vacancy concentration, and l is the chemical potential. Assuming that the vacancy concentration is at equilibrium, the stress-dependent contribution to the chemical potential is Dl r n ð Þ ¼ Àr n X, in which r n is the local normal stress on the grain boundary [39]. The stress distribution along the grain boundary with an ellipsoidal creep cavity is calculated by the finite element method using COMSOLÒ [27]. The Linear Elastic Material module with stationary node and nodal serendipity (quadratic) elements is used, and the elastic equilibrium equation can be expressed as r Á r þ F ¼ 0, in which F is the body force per unit volume, while r is the Cauchy stress tensor and can be calculated with the Hooke's law r ¼ C : e, where e is the elastic strain tensor and C is the fourth order of stiffness tensor. A constant stress with the magnitude of r 0 is applied uniformly to the top edge of the simulation volume (Fig. 1a). The material is regarded as isotropic with a density q of 7800 kg/m 3 . At a temperature of T = 550°C (823 K), the Young's modulus E is taken to be 165 GPa [40] and the Poisson's ratio m is set to 0.33 [40]. The governing equations and boundary conditions are also summarised in the supplementary material. Figure 2 shows the normal stress r n (i.e. the z component of the stress tensor) distribution as a function of (a/r) for different (k/a) ratios, in which r is the distance from the centre of the (empty) cavity with radius a. Note that this stress is normalised by the applied stress r 0 . A stress concentration occurs at the edge of the cavity (a/r) ? 1 while r n approaches the applied stress r 0 at the far end of the grain boundary with (a/r) ? (a/k). The stress distribution shows similar profiles along the grain boundary for the cases where (k/a) C 5. Therefore, the stress distribution is fitted with a simple phenomenological expression independent of (k/a): where A = 1.65(2) 9 10 -3 and B = 6.82(3) are dimensionless constants. The numbers in brackets are the uncertainty in the last digit. The inward vacancy flux in Eq. (1) can now be estimated using Eq. (2) and is proportional to the applied stress r 0 .
The equilibrium vacancy concentration at the grain boundary corresponds to [41]: where Q is the vacancy formation energy at the grain boundary, a = 1.28 is a dimensionless proportionality constant, and s ¼ ð1ÀsÞ X boundary Q = Q gb & Q m /2 is taken as half of that in the bulk with Q m = 0.58 eV [43]. The vacancy diffusivity at the grain boundary is taken as the Fe grainboundary diffusivity of the Fe host atoms [43], i.e. D v gb = D h gb = D Fe gb = 1.74 9 10 -12 m 2 s -1 at a temperature of 550°C. The equilibrium vacancy concentration at the grain boundary of bcc Fe then corresponds to x eq v & 7.4 9 10 -4 at 550°C.

Solute diffusion and cavity closing
When a grain-boundary cavity, and hence an open volume, is present, then healing solute atoms start to diffuse towards the cavity in order to precipitate at the cavity surface leading to a reduction in the level of supersaturation. The diffusion of solute atoms is driven by the solute concentration gradient, which can be described by where D s and x s are the diffusivity and the concentration of the solute (atom fraction), respectively. In a binary alloy system containing also vacancies, the mass balance gives x s ? x h ? x v = 1, where x h is the concentration of host atoms. Assuming that the vacancy concentration is low everywhere in the system, and at equilibrium, the flux of the host atoms can be written as The difference in diffusivity between the host and the solute atoms results in a net vacancy flux, i.e. the Kirkendall flux: The vacancy flux is oriented opposite to the direction of the fastest diffusing component. In a selfhealing system, the solute atoms need to show a higher diffusivity than the host atoms in order to generate a vacancy flux oriented outwards of the cavity. 'Vacancies' making up the open volume of the pre-existing cavity can be removed from the cavity when the Kirkendall vacancy flux J K v (oriented outwards) is higher than the stress-driven vacancy flux J r v (oriented inwards).
The Transport of Diluted Species module with linear Lagrange elements in COMSOLÒ [27] is used to simulate the solute transport to the cavity. The governing diffusion equation without source is ox s ot þ XrJ s ¼ 0, where x s is the concentration of the solute s (atom fraction) and J s is the flux (atom•m -2 •s -1 ) and can be calculated by The solute concentration at the cavity edge is maintained at 0.07 at.%, while the initial solute concentration in the bulk and the grain boundary is set at 0.25, 0.5, 1, 2, 3, and 4 at.%, respectively. The Au diffusivity in bcc Fe at 550°C (7.39 9 10 -19 m 2 /s) [30] is used as the solute diffusivity in the bulk, while the self-diffusivity for bcc Fe at 550°C (1.50 9 10 -21 m 2 /s) [30] is set as the host diffusivity in the bulk. In order to estimate the effect of (1) the grain-boundary diffusivity and (2) the interspacing between two neighbouring cavities on the solute transport efficiency, the grain-boundary diffusivity is set to be 10 n times the diffusivity in the bulk, in which n = 1-9, while the length of the simulation box is set to 2.5, 5, 10, and 25 lm, corresponding to the relative inter-cavity spacing (k/a) of 5, 10, 20, and 50. The simulation time is 1000 h for all cases. In total 8319, 14,527, 21,601, and 26,113 elements with 7395, 13,456, 20,398, and 24,800 vertices are used for the simulation box with the length of 2.5, 5, 10, and 25 lm, respectively. The governing equations and boundary conditions are also summarised in the supplementary material. After calculating the solute transport to the cavity, the net vacancy flux (i.e. the Kirkendall flux) is calculated by Eq. (5).

Results and discussion
Solute transfer profile Figure 3a-c shows the solute concentration profile for different D gb /D m ratios after 1000 h, for a condition where the supersaturation is 1 at% and the relative inter-cavity spacing is k/a = 20. The white arrows indicate the transport direction of the solute. In Fig. 3a, it can be seen that for a relatively low grainboundary diffusivity (D gb /D m = 10) the solute transport has a 3D nature: the solute concentration contour reflects the geometry of the cavity, indicating a more or less uniform solute transport in both the matrix and in the grain boundary. As shown in Fig. 3c, for a relatively high value for the grainboundary diffusivity (D gb /D m = 10 8 ) the concentration profile reflects primarily the grain-boundary geometry (instead of the cavity geometry), indicating the 1D nature of the solute transport towards the grain boundary. In this case, the grain boundary provides a fast diffusion path. As a result, the solute in the matrix tends to diffuse vertically towards the grain boundary with a vertical diffusion length L ? ¼ . The 1D nature of the diffusion pattern can also be induced from the direction of the diffusion indicated by the white arrows: in Fig. 3c, the solute diffusion is approximately perpendicular to the grain boundary, while in (a), the diffusion direction is perpendicular to the cavity edge. As shown in Fig. 3b, an intermediate value for the relative grain-boundary diffusivity (D gb /D m = 10 4 ) provides a diffusion pattern that reflects a crossover between 1 and 3D solute diffusion. The present findings are in line with previous results by Versteylen and co-workers [29].

Vacancy fluxes
As a result of the Kirkendall effect, the solute diffusion generates a vacancy flux in the opposite direction. Thereby the vacancies present in the cavity can be removed through either the cavity/matrix interface, or the cavity/grain-boundary interface. Figure 3d shows the percentage of the vacancy flux through the cavity/grain-boundary interface. It is clear that the higher the grain-boundary diffusivity, the larger fraction of vacancies that is removed via the grain boundary. In previous work [22], the ratio for Au diffusivity in the grain boundary D Au gb and in the bcc Fe matrix D Au m À Á was estimated to be D Au gb =D Au m ¼ 10 5 À10 6 , meaning that the grain boundary carries approximately 99% of the transport capacity for (k/a) C 20 (Fig. 3d), 91% for (k/a) = 5, and 97% for (k/a) = 10. Figure 4 shows the evolution of the inward stressdriven vacancy flux (dashed lines) and the outward Kirkendall vacancy flux (solid lines) as a function of time. The inward vacancy flux is constant in time and proportional to the applied stress, while the outward Kirkendall flux decreases with time and is affected by the grain-boundary diffusivity D gb , the supersaturation Dx, and the half inter-cavity spacing k. The Kirkendall flux, which is proportional to the solute flux transported to the cavity, shows a power-law relationship with time J K v / t n , with an time exponent n ranging from n = -1/2 to n = 0, depending on (D gb /D m ) and (k/a). From the value of the time exponent n, one can estimate whether the solutedriven vacancy diffusion has a 1D, 2D, or 3D character [29]. For instance, in Fig. 4a, for a low diffusivity ratio value (D gb /D m = 10), the comparable diffusivity in the bulk and in the grain boundary leads to a solute diffusion field that is almost uniform in the bulk and in the grain boundary (as shown in Fig. 3a), indicating the 3D character of the diffusion. This is consistent with the time exponent n, which shows an increase from an initial value of about -1/2 and approaches 0 for longer times. It is worth to note that for an ideal 3D diffusion, the solute flux is constant for longer time scales, but even for a low diffusivity ratio, the 3D diffusion character will eventually breakdown once the depletion zone reaches the edge of the simulation box. After that, the depletion of the solute is no longer uniform throughout the matrix and the grain boundary, and the solute flux (as well as the Kirkendall vacancy flux) can no longer remain constant. For a high grain-boundary diffusivity (D gb / D m = 10 6 -10 9 ), the solute in the grain boundary depletes first, and during this stage, the outward vacancy flux is determined by the grain-boundary diffusivity D gb and the width of the grain boundary d. This process can be simplified as a 2D diffusion, where the solute atoms are mainly transported from the grain boundary (instead of the matrix) to the cavity surface and the time exponent for J K v corresponds to n = 0. After that, the grain boundary acts as a fast diffusion path, through which the solute in the bulk can be transported towards the cavity, and the diffusion pattern experiences a transition from 2 to 1D. Once the 1D diffusion pattern is developed, the bulk solute diffusivity becomes the limiting parameter for the solute transport, while the grain-boundary diffusivity no longer limits the diffusion. Therefore, for high D gb /D m ratios (at fixed values of k/a and supersaturation Dx) on longer time scales, the flux into the cavity approaches to the same level and the time exponent stabilises at n = -1/2.
The 2D to 1D transition can be seen clearly in Fig. 4b: for a relatively large value of k/a (see e.g. k/ a = 50), the time exponent n decreases gradually from 0 to -1/2. For a fixed grain-boundary diffusivity D gb , the transition time for the crossover from 2 to 1D behaviour is determined by the time required for the depletion of the solute in the grain boundary. This crossover time from 2 to 1D behaviour t 2D-1D can be estimated by comparing the diffusion length 2 ffiffiffiffiffiffiffiffiffi D gb t p with the distance between the cavity and the edge of the simulation box k-a, resulting in t 2DÀ1D % 4 k À a ð Þ 2 =D gb % 4k 2 =D gb . For relatively small values of k/a (see k/a = 5 and 10), the 1D diffusion pattern is developed within relatively short time scales, so that the time exponent starts from -1/2 in our simulations. During the 2D diffusion stage, the outward vacancy flux, which is limited by the grainboundary diffusivity D gb , is the same for different k/ a ratios as the diffusion profile has not reached the edge of the diffusion box. While during the 1D diffusion stage, where the matrix diffusivity D m is the limiting parameter, the outward vacancy flux is proportional to the (k/a) 2 . It is observed from Fig. 4c that the supersaturation Dx does not change the nature of the diffusion, but a higher supersaturation leads to a higher concentration gradient, and therefore a larger solute flux into the cavity.
As explained in Sect. 2, the open volume of the cavity shrinks when the Kirkendall vacancy flux is larger than the stress-driven vacancy flux and grows otherwise. The time integrated difference between the two vacancy fluxes equals to the amount of vacancies removed from the pre-existed cavity at the corresponding time (in mol). Since the stress-driven vacancy flux is assumed to be constant, while the outward vacancy flux decreases with time, a crossover takes place at a certain time t cr , at which the possible maximum filling is achieved and the cavity, if not fully healed, starts to grow although solutes remain to be transported to the cavity, as illustrated in Fig. 4d. This is different from the previous research [9]. In the previous study, both the Kirkendall flux and the stress-induced vacancy flux were assumed to be time independent. Therefore, fully filling can be achieved as long as the (initial) Kirkendall flux is higher than the stress-induced flux. In the present research, however, it is possible that the cavity cannot be fully filled, even if the Kirkendall flux initially exceeds the stress-induced flux before t cr . As illustrated in Fig. 4d, the cavity can only be fully filled if the net amount of vacancies (indicated by the green shadow) is higher than the volume of the original cavity.

Filling ratio
Starting from a pre-existing cavity with an initial volume V 0 , hypothetical 'vacancies' making up the cavity are being removed by the Kirkendall flux J K v [ 0 À Á and real vacancies coming from the solid matrix are added by the stress-driven vacancy flux J r v \0 À Á . In the current model, we do not actually calculate the displacement of the original pore-matrix boundary, but keep it constant and only calculate the total volume of solute atoms entering the cavity and forming the precipitate. We make no assumptions on the shape of the precipitate, nor do we describe the precipitate/remaining empty pore interface. In this sense, the model is a degenerate moving boundary or Stephan problem [44].
The filling ratio can now simply be defined as the ratio between the net amount of the removed vacancies and the amount vacancies in the original cavity (V/X): where S is the cavity surface and t the time (note that the positive direction for the flux through the closed surface S is pointing outwards). It is assumed that the stress-driven vacancy flux only enters the cavity J Mater Sci (2022) 57:12034-12054 through the grain-boundary/cavity interface, while the Kirkendall flux can enter the cavity from both the grain boundary and from the matrix. Figure 5 compares the time evolution of the filling ratio for different values of the supersaturation Dx, the grain-boundary/matrix diffusivity ratio D gb /D m , the relative inter-cavity spacing k/a, and the applied stress r 0 . We can clearly distinguish two sorts of overall behaviour: for the first one, the cavity can become fully filled (i.e. proceeds to a FR of 1.0) at the healing time t h and the reaction stops. For the second type, the filling ratio initially increases and then peaks at a critical time t cr at a level below FR = 1, which corresponds to the maximum filling ratio that can be achieved. For t [ t cr , the filling ratio again decreases, indicating that the open volume of the cavity starts to grow faster than can be made up by the solute transport. When the filling ratio eventually reaches a negative value (FR \ 0), the open volume of the partially filled cavity is larger than its original size. For the sake of completeness, the dashed line segments in the Figs. 5 for conditions at which the condition FR = 1 is reached show the hypothetical behaviour assuming that achieving the completely filled state does not lead to a change in transport rate. The filling ratio demonstrates the competition between the outward Kirkendall vacancy flux and the inward stress-driven vacancy flux. Generally, at a given time, a higher filling ratio can be achieved with a higher supersaturation, a higher grain-boundary diffusivity, and a larger relative inter-cavity spacing, while a higher stress level results in a lower filling ratio. Figure 6 shows the maximum filling ratio for various combinations of parameters. In Fig. 7, the time of healing t h (solid line, for the case where the cavity can be fully filled) and the critical time t cr (dotted line, for the case where the cavity cannot be fully filled) are presented. Since the stress-driven inward vacancy flux eventually will exceed the outward Kirkendall vacancy flux, for proper healing it is crucial that the cavity becomes free of any remaining 'empty' volume before t cr is reached.

Critical stress
From Figs. 6 and 7, it is clear that there is an upper threshold for the applied stress, above which the cavity cannot be fully filled. In Fig. 8, the calculated critical stress as a function of the grain-boundary/ bulk diffusivity ratio D gb /D m is shown for different values of the relative inter-cavity spacing k/a and the supersaturation Dx. It can be seen that both the grainboundary diffusivity and the inter-cavity spacing play a significant role in the critical stress. The critical stress increases with the grain-boundary diffusivity D gb , until it reaches a level where a quasi-1D diffusion pattern exists and the bulk diffusivity becomes the limiting diffusion rate. When D gb is high enough so that D m becomes the limiting diffusion rate, a larger inter-cavity spacing k results in a higher critical stress, since the solute atoms are collected and transported towards the cavity over a larger grainboundary surface area p(k 2a 2 ). With a small relative inter-cavity spacing (e.g. k/ a = 5), the critical stress level is expected to be low even if the grain-boundary diffusivity is sufficiently high. This is the situation for the late creep stages: more cavities are located at the grain boundary, resulting in a decrease in the cavity interspacing. This further leads to a fast growth of the open volume of the cavity and the subsequent coalesce of neighbouring cavities, which causes the catastrophic failure of the material. In the initial creep stage, where the grain-boundary cavities have a small size and the inter-cavity interspacing is large, the cavities can still be filled relatively easy. In the fully filled cavities, the open volume disappears and the open volume will need to nucleate somewhere else at the grain Maximum filling ratio Figure 6 Maximum filling ratio as a function of applied stress r 0 for different values of the supersaturation Dx, the grain-boundary diffusivity to matrix diffusivity ratio D gb /D m , and the relative inert-cavity spacing k/a. For a maximum filling ratio equal to or higher than 1, the cavity can be fully filled under the corresponding combination of parameters. Note that the curves overlap for D gb / D m [ 10 6 .
boundary, in which case the vacancy and solute diffusion start over again. From Fig. 8, it can be deduced that the critical stress can be raised by increasing (a) the supersaturation level, and (b) the diffusivity difference between the healing solute and the host. In practice, the intended supersaturation level should not be higher than a certain limit dictated by the criteria that the healing atoms are fully soluble in the matrix at the homogenising temperature and that the strain penalty between the matrix and the precipitate is large enough to suppress solute consuming precipitation in the bulk. Furthermore, the diffusivity for the solute should always be higher than the diffusivity of the host, while a higher solute ratio is beneficial.
For the situation of 3D diffusion, the grainboundary diffusivity and the bulk diffusivity should be comparable, and the solute atoms can be transported to the cavity if their distance to the cavity surface is less than the diffusion length L ¼ 2 ffiffiffiffiffiffiffi D s t p (D s being the diffusivity of the solute atoms), resulting in a depletion zone of the healing atoms resembling the shape of the cavity, and the critical stress is independent of the interspacing distance k. For the situation of 1D diffusion, as explained in Sects. 3.1 and 3.2, all the solute atoms within a volume p k 2 À a 2 À Á L ? on both sides of the grain boundary can be transported to the cavity (where L ? is the vertical diffusion length determined by the solute diffusivity in the bulk D s m ), and therefore a large k significantly increases the healing capacity. In practice, a quasi-1D diffusion pattern is observed [22] since the grainboundary diffusivity is usually orders of magnitude higher than the bulk diffusivity.
The effect of the supersaturation on the critical stress is shown in Fig. 9. A scaling relation r Critical / Dx ð Þ m applies, with m & 1 for a lower D gb /D m ratio (indicating a 3D diffusion), and m & 2 for a higher D gb /D m ratio (indicating a 1D diffusion), respectively. The 3D to 1D transition is especially obvious with a lower interspacing/pore radius ratio k/a = 10, since less time is required to develop the 1D diffusion pattern for a lower interspacing-radius ratio. In fact, when k/a = 10, the bulk diffusivity becomes the limiting parameter when D gb /D m C 10 6 , while for a higher interspacing-radius ratio k/a = 50, the grainboundary diffusivity D gb still plays a role on the critical stress when D gb /D m = 10 8 .
The interspacing-radius ratio k/a has a crucial influence on the critical stress, as demonstrated in Fig. 10. A scaling relationship r Critical / k=a ð Þ k is obtained with k = 0-4. For a lower D gb /D m ratio, which corresponds to a 3D diffusion regime, the critical stress is almost independent of k/a. In this situation, the healing atoms are transported through the bulk, and whether a solute atom can be transported towards the cavity is determined by its distance to the cavity surface. Therefore, the interspacing distance k does not significantly affect the healing kinetics. However, a higher D gb /D m ratio results in a 1D diffusion, where a solute atom can be transported towards the cavity through the grain boundary (instead of through the bulk) as long as its vertical distance to the grain boundary is within the diffusion length L ? ¼ 2 ffiffiffiffiffiffiffiffi ffi D m x t p . As a result, a larger k means that more solute atoms can be used to heal the cavity, leading to a higher critical stress. It is worth to note that for stress levels up to 300 MPa, the growth of the cavity can be regarded as being controlled by vacancy diffusion, as discussed in Sect. 2.2. This  stress range is larger than the stress range we are interested in, and also covers the applied stresses in previous experimental research [19,22,25,26]. With an increase in applied stress, the diffusion length decreases and the cavity growth mechanism shifts towards creep flow. The results for high stresses above 300 MPa are also presented in Figs. 6,7,8,9 and 10 for completeness.

Analytical model
From the finite element model simulations, we can identify the dominant mechanism for the self-healing of creep cavities by supersaturated solute corresponds to 1D matrix diffusion of the solute towards the grain boundary and subsequently transport of this solute though the grain boundary towards the free cavity surface connected to the grain boundary. This 1D solute diffusion is observed for the following conditions: The first condition indicates that the grain boundary acts as a much faster transport route for solute than the matrix. In order to reverse the vacancy flux by the Kirkendall effect, the diffusivity of the solute atoms needs to be higher than that of the host atoms in the matrix. The second condition requires that the grain boundary is largely depleted from supersaturated solutes as the diffusion length exceeds the distance between the creep cavity surface and the edge of the simulation box (half of the inter-cavity spacing). This condition can only be met for a finite inter-cavity spacing, which is generally the case for practical situations. The second condition can also be translated in a time requirement: For k = 20 lm and D gb-= 2.4 9 10 -14 m 2 s -1 (at 550°C) [22] this characteristic time for the onset of 1D diffusion amounts to t 2D-1D & 400 s. Under these conditions, an approximate analytical model will be constructed, for which the results can than subsequently be compared to those of the full finite element model.
The time evolution of the cavity depends on the competition between the stress-driven vacancy flux and the Kirkendall vacancy flux. Combining Eqs. (1) and (2), the stress-driven vacancy flux at the interface between the cavity surface and the grain boundary ðr ¼ aÞ can be written as: where A and B are the fitting parameters for the stress distribution along the grain boundary.
This inward stress-driven vacancy flux is constant in time given a stable stress distribution. The integrated stress-driven flux provides the time evolution of the open volume of the cavity: where V 0 / a 3 is the initial cavity volume and a ¼ For 1D solute diffusion in the matrix towards the grain boundary, the difference in diffusivity between the solute and the host atoms results in a 1D Kirkendall vacancy flux at the interface between the matrix and the grain boundary. When the solute transport is much faster in the grain boundary than in the matrix, the solute concentration in the grain boundary is nearly constant and close to its equilibrium value. In that case, the 1D Kirkendall vacancy flux at the interface between the matrix and the grain boundary can be estimated by [29]: where D s m and D h m are the matrix diffusivities for the solute and the host atoms, respectively. This Kirkendall vacancy flux is time dependent and for D s m [ D h m oriented from the grain boundary towards the matrix. The integrated vacancy flux provides the time evolution of the open volume of the cavity: where the factor arises from the vacancy transport at the top and bottom interfaces of the grain boundary Combining all contributions to the time evolution of the open volume of the cavity gives: Healing of a creep cavity is achieved when the open volume of the cavity vanishes as a result of the solute transport at V = 0. The lowest time solution is The critical maximum stress for which a solution can be obtained is: Equation 12 agrees well with the results in Figs. 8, 9 and 10. It is worth to note that the numerical results are all based on a fixed initial cavity radius, i.e. a = 0.5 lm, which corresponds to the average size of the creep cavities in a Fe-Au system with a stress level of 100 MPa after 10-20 h creep [19]. For a cavity with smaller size (which is usually the case at early stage of creep), the healing time for the cavity will be shorter, and the critical stress increases accordingly. From Eq. 12, the critical stress scales as Dx ð Þ 2 / a À3 . When the radius of the cavity a is half of the value, the critical stress will be increased by a factor of 8 and the healing time to fully fill the cavity will be decreased by a factor of 8.
The corresponding critical time is related to the critical stress by at cr r cr ¼ V 0 . As shown in Fig. 11, the value of at cr r cr =V 0 is approximately constant for a low inter-spacing ratio k/a. For an inter-spacing ratio k/a = 5 and a diffusivity ratio D gb /D m higher than 10 5 , a value of at cr r cr =V 0 = 2.3 is found, which indicates that the simplified analytical model is a reasonable approximation. Assuming that the region where the processes are dominated by 1D diffusion is marked by a value of at cr r cr =V 0 \ 10, then it can be seen that for diffusivity ratio D gb /D m higher than 10 7 , the 1D diffusion behaviour is found for k/a \ 30. For lower values of the diffusivity ratio (D gb /D m \ 10 7 ), the critical k/a value for 1D diffusion behaviour continuously decreases and reaches a value of k/a & 10 for D gb /D m = 10 4 . This behaviour is qualitatively in line with the 1D diffusion results in the absence of stress [29]. The critical time is the maximum time for which a solution can be obtained and this scales as: For a stress below the critical stress r 0 \r cr ð Þ , the time to complete filling t h can be expressed in terms of the critical time t cr and the critical stress r cr : The healing time grows continuously to t h ¼ t cr for r 0 ¼ r cr . These results of the simplified analytical model for the case of 1D solute diffusion are qualitatively in line with the results of the finite element model calculations. It is worth to recall that these time scales need to be large compared to the characteristic time scale t 2DÀ1D % k 2 =ð4D gb Þ, which is required to deplete the solute in the grain boundary and establish the 1D diffusion pattern.
The predictions of the analytical model are compared with the experimental tomography data in ref [19]. The results for the volume of precipitates, volume of open cavities, and the ratio between the volume of the precipitate and the sum of that volume and the remaining empty cavity V prec with time is shown in Fig. 12. The cavity radius at t = 0 is taken as 0.5 lm. The inter-cavity spacing is estimated based on the number density of cavities and the grain size in the Fe-Au alloy. Assuming that the cavities are uniformly spaced on the grain boundaries, the half inter-cavity spacing k, the grain size d g , and the number density of cavities N 0 (lm -3 ) fulfil the following relationship 1 4k 2 ¼ N 0 d g . Taking the grain size d g = 57 lm and number densities from ref [19], the half inter-cavity spacing yields 7.9, 7.6, 4.8, and 6.9 lm, respectively. In Fig. 12, the half intercavity spacing is taken as 5.7 lm, i.e. k = d g /10. The stress level is taken as 100 MPa. The calculated volumes of both the precipitates and the open cavities Figure 11 The value of at cr r cr =V 0 for as a function of the relative inter-cavity spacing k/a and the grain-boundary/bulk diffusivity ratio D gb /D m for a supersaturation Dx of 4 at.%. A value close to 1 represents a 1D diffusion pattern. The contour where at cr r cr =V 0 ¼ 10 is plotted to indicate the boundary for the 1D diffusion regime.
have the same magnitude as the experimental data, while the decrease in the V prec V prec þ V open À Á ratio with time agrees well with the model. It is worth to note that the tomography data were taken from the creep-failed experiments. At an earlier stage of creep, when the cavities are smaller and the inter-cavity spacing is larger, it is likely that the cavities can be fully filled.
In Fig. 13, the results for V prec , V open and V prec V prec þ V open À Á are calculated using the experimental data (i.e. half inter-cavity spacing k = 7.9, 7.6, 4.8 and 6.9 lm; the stress level and lifetime are 60 MPa (642 h), 80 MPa (376 h), 100 MPa (210 h), and 117 MPa (57 h), respectively) and are compared with the experimental results. Again, the calculation is in good agreement with the experimental data.

Conclusions
We have developed a set of models which predicts the growth and filling of grain-boundary creep cavities in self-healing binary alloys as a function of time, stress, level of supersaturation, and inter-cavity spacing and diffusion ratios. The competition between the inward and outward vacancy fluxes results in the growth or shrinkage of the open volume of the creep cavity. It is found that the filling ratio shows a maximum value at a critical time t cr , which corresponds to the time when the inward vacancy flux exceeds the outward vacancy flux (integrated over the surface area). Two conditions can be distinguished: one where the cavity becomes fully filled before this critical time is reached (and the process stops) and one where the only partial filling will take place and the growth of the cavity will continue. For each combination of parameters, the critical applied stress r cr is calculated, below which the cavity can be fully filled. The analytical model, fed by insights from the numerical model, shows that for conditions leading to quasi-1D solute transport, the critical stress approximately scales as r cr / k=a ð Þ 4 D s m =D gb À Á Dx ð Þ 2 , where k=a is the relative inter-cavity spacing, D s m =D gb the ratio of the solute diffusivity in the matrix and the vacancy diffusion in the grain boundary and Dx is the supersaturation of the solute. The critical healing time obtained at the critical stress scales as t cr / 1=r cr , while for lower stresses the pore filling time t h can be expressed in terms of t cr and r cr . The results from the analytical model are in good agreement with the results from the nanotomography experiments.