A mesoscale phase-field model of intergranular liquid lithium corrosion of ferritic/martensitic steels

Abstract

A phase-field model is developed to simulate intergranular corrosion of ferritic/martensitic steels exposed to liquid lithium. The chromium concentration of the material is used to track the mass transport within the metal and liquid (corrosive) phase. The framework naturally captures intergranular corrosion by enhancing the diffusion of chromium along grain boundaries relative to the grain bulk with no special treatment for the corrosion front evolution. The formulation applies to arbitrary 2D and 3D polycrystalline geometries. The framework reproduces experimental measurements of weight loss and corrosion depth for a 9 wt% Cr ferritic/martensitic steel exposed to static lithium at 600 $\degree$C. A sensitivity analysis, varying near-surface grain density, grain size, and chromium depletion thickness, highlights the microstructural influence in the corrosion process. Moreover, the significance of saturation is considered and evaluated. Simulation results show that near-surface grain density is a deciding factor, whereas grain size dictates the susceptibility to intergranular corrosion.

1 Introduction

The commercialization of fusion energy remains hindered by many overwhelming engineering problems that need addressing prior to its implementation. A key one is the feasibility of breeding tritium in situ using lithium-bearing materials. Liquid lithium (Li) provides unique benefits [1, 2]. The concern of radiation damage is eliminated as the breeder material is in a liquid state, which offers greater flexibility and efficacy in the extraction of tritium as the breeder material can circulate away from the extreme conditions of the reactor. Furthermore, the inherent efficacious thermal properties of liquid metals indicate that one can operate at higher temperatures, leading to greater thermal efficiencies [3, 4]. Several potential candidates for the optimum structural material for future fusion power plants are currently being investigated, including ferritic/martensitic (F/M) steels and vanadium-based alloys [5, 6, 7]. Each candidate material possesses its respective advantages, yet owed to the worldwide expertise and generations of use in the fission industry, many of the proposed designs have opted for F/M steels [8], thereby exploiting the high strength and low creep rate at elevated temperatures plus desirable irradiation performance [9].

A limiting concern to the implementation of liquid Li as a breeder material is the harsh corrosive environment it produces when in contact with structural alloys, ultimately leading to weight loss and surface recession (i.e., wall thinning) influencing the longevity of the structural material and heightening the risk of accumulated activated material [10]. The corrosion of alloys with liquid Li is solely based on physio-chemical processes, where the principal mechanisms are solubility-driven dissolution, intergranular corrosion (IGC), and mass transfer [11, 12]. As such, the composition of the structural alloy in contact with the liquid metal is a substantial factor in the corrosion resistance, owed to the varying solubility limits of alloying elements differing by several orders of magnitude in some cases [13]. Furthermore, Li preferentially attacks grain boundaries (GBs) due to their heterogeneity, tendency for chemical segregation, and precipitation of secondary phases. The compatibility of F/M steels with liquid Li is heavily correlated with the amount of Cr in the composition of the steel [14]. Cr is preferentially leached from the structural alloy attributed to its comparatively heightened solubility limit [14, 15]. The amount of Cr regulates the chemical activity of carbon in the metal by forming metallic carbides (i.e., Cr₂₃C₆). Liquid Li preferentially corrodes GBs as these regions are decorated with carbides, relative to a tempered martensitic microstructure [16]. Once exposed to liquid Li, these carbide precipitates become unstable and are subsequently leached [17] via the following reaction

$$\text{Cr}_{23}\text{C}_{6\,(\mathrm{s})}+6\text{Li}_{(\mathrm{l})}\rightarrow 3\text{Li}_{2}\text{C}_{2\,(\mathrm{sol})}+23\text{Cr}_{(\mathrm{sol})}$$ (1)

resulting in a heavily chemically altered region [18], commonly observed through a phase transformation of the matrix from martensite to ferrite [19]. Once the alloying elements in the immediate vicinity are depleted, resulting in a chemical gradient, bulk alloying elements diffuse towards the corroding interface to replenish the lost material. Although the compatibility of F/M steels with liquid Li at elevated temperatures has been experimentally explored in depth [20, 21], showcasing its applicability, there still remains little consensus on the dominant corrosion mechanism and the principal microstructural features which govern the corrosion process.

Evaluating IGC resistance and the long-term effects of materials exposed to liquid Li presents significant challenges due to the unique conditions of this corrosive environment. Numerical modelling can expand the reach of experimental studies by providing insights and extending the duration of experiments involving liquid Li under various conditions. Ideally, these models can guide the design of materials optimized for liquid Li environments and suggest strategies to minimize susceptibility to IGC. Modelling the corrosion resistance of structural materials has been detailed in abundance through various computational techniques, including Lagrangian-Eulerian approaches [22, 23], peridynamics [24, 25], cellular automata models [26, 27] and phase-field methods [28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. However, these are mainly used in the context of uniform and localized corrosion and stress corrosion cracking. As such, these models are limited in capturing IGC, assessing the compatibility of liquid Li with F/M steels, and incorporating the microstructural influence on IGC. The synergistic effect of aggressive environments and underlying microstructures has been recently resolved in the context of IGC [39, 40, 41, 42, 43]. Albeit these IGC models provide mechanistic predictions and insight, they do not properly resolve the differing corrosion kinetics in the grain bulk and GBs [39, 40, 41], a key feature of liquid Li IGC. Other models utilize multiphase-field formulations [42, 43] to distinguish between GBs, grain interiors, and the corrosive environment. However, these multiphase-field approaches inherently possess high computational costs. Moreover, neither model surveyed the correlation between microstructural properties and their influence on the IGC process. The present work aims to develop a phase-field model to assess the IGC of F/M steels in contact with liquid Li. As the underlying physical mechanics in the corrosion process, the diffusion process is differentiated between GBs and grain bulk, effectively and naturally addressing the varying corrosion kinetics while keeping the phase-field equation unaltered. An additional stationary parameter is introduced to distinguish mass transport between GBs and grain interiors, which reduces the computational cost intrinsically associated with multiphase-field formulations.

The remainder of the paper is organized as follows. The underlying corrosion mechanism of F/M steels exposed to liquid Li and the associated modelling assumptions are presented in the following section. The chromium (Cr) concentration of the material is used to track mass transport within the metal and liquid (corrosive) phase. The phase-field model is subsequently derived from a generalized thermodynamic free energy functional. The governing equations for the phase-field parameter and Cr concentration in the system, as the primary variable of the model, are derived in a thermodynamically consistent way. GBs are differentiated from grain interiors using an additional independent stationary parameter that allows for incorporating representative diffusivities of Cr along GBs and grain bulk, exploiting physical material parameters to capture IGC. The constructed framework is calibrated and validated against experimental measurements on a 9 wt% Cr F/M steel in contact with static liquid Li at 600 $\degree$C in Section 4.2, demonstrating its ability to capture IGC phenomena. After validation, a sensitivity analysis for different microstructural features, including grain density at the exposed surface, grain size, and Cr depletion thickness, is performed to investigate their role in the corrosion resistance of the simulated material in Section 2. The main conclusions are subsequently discussed in Section 3 as well as the investigation findings and recommendations for future work.

2 Results

2.1 Influence of grain density at the exposed surface

The microstructures employed to calibrate the model had 6 GBs at the exposed surface. To understand the influence of grain density at the near-surface, the grain structure distribution is varied while keeping an average grain size of 20 $\upmu$m. In addition to the reference geometry in Section 4.2, two further microstructures, with 5 and 7 GBs at the exposed surface, were analysed. As previously, the average behaviour was taken from ten equiaxed microstructures for each case. Fig. 1 displays the representative microstructure for each scenario while the ten microstructures used for the 5 GB and 7 GB microstructures can be viewed in Figs. S.2 and S.3 (Supplementary Information). The resultant phase-field predictions can be seen in Fig. 2 where the standard deviation (SD) data at 100, 250, and 500 hours exposure time are depicted in Fig. 3.

The weight loss data highlights a positive correlation with the number of grains at the exposed surface. The microstructures that possess 7 GBs exhibited the most significant weight loss over the 500-hour simulation. Contrastingly, the microstructures with 5 GBs give rise to the lowest weight loss. Overall, all three microstructures with 5, 6, and 7 GBs at the exposed surface produce largely identical profiles, whereby the difference in weight loss between the three types of microstructures increases with increasing exposure time. As such, it is apparent during the first 500 hours of exposure that the rate of weight loss is most severe for the 7 GB microstructure, resulting in a comparatively greater weight loss of 1.34 g/m². The 5 GB microstructure, on the other hand, exhibited 1.01 g/m² of material weight loss. On average, the total weight loss after 500 hours increases by $\sim$15% each time, increasing the grain density at the interface by a single grain. Consequently, the 7 GB microstructures overestimated the weight loss at 100 and 250 hours compared to the experimental measurements (Section 4.2). Yet, it remains within the relative error for the 250-hour weight loss experimental measurement. Alternatively, the 5 GB microstructures predicted with a degree of high accuracy the weight loss at 100 hours. However, it underestimates the weight loss at 250 hours, while remaining within the relative error. Furthermore, after 500 hours, the 5 and 7 GB microstructures give rise to almost identical SDs, which are more than twice as large compared to the 6 GB microstructures, Fig. 3(a). The corrosion depth data displays zero influence on the near-surface grain density as all three types of grain structures give rise to predominantly identical corrosion depth profiles. Overall, the three microstructures simulated generate an average corrosion depth of 12.3 $\upmu$m. Contrary to the weight loss data, the SDs regarding the corrosion depth exhibit a clear order. The 5 GB microstructures produce the greatest SD, followed by the 6 GB microstructures, where the 7 GB microstructures return the smallest SD after 500 hours of exposure time, Fig. 3(b).

Fig. 4 shows the progression of corrosion of the representative microstructures (Fig. 1) following a 30,000-hour exposure time. The onset of corrosion initiates at the GBs exposed to the liquid metal at the upper most boundary. As time increases, the evolution of $\phi$ continues along the GBs deeper into the material. Congruently, Fig.5 shows the diminished concentration of Cr intensifies as Li advances into the steel. The spread of concentration of Cr along the GBs indicates that the corrosion process operates in the diffusion-controlled regime. The near-surface GBs are seen to be completely leached of Cr, with this process increasing in depth as time progresses. In addition, the degree of corrosion regarding the near-surface GBs visibly heightens, as seen through the increased thickness of the $\phi=0$ region. In some instances, two corrosion fronts are seen advancing at either end of the same GB, where they ultimately meet, engulfing the entire grain. The expansive nature of Li to diffuse across all GBs is evident, as it proceeds deeper into the material, clearly leaving behind a deteriorated material. After 30,000 hours of exposure time, the liquid Li corrodes most GBs in the simulated solid phase. However, it only causes marginal corrosion to the grains on the exposed surface.

2.2 Influence of grain size

The grain size plays a governing role in IGC due to the synergistic relationship with GB density within the bulk of the material. As such, the average grain size is altered to observe this dependency. Two additional average grain sizes of 10 and 40 $\upmu$m are modelled. To isolate the effect of grain size from near-surface grain density, all three average grain size microstructures should possess the same number of GBs at the exposed surface. Nonetheless, the feasibility of a 10 and 40 $\upmu$m equiaxed microstructures possessing the same number of grains at the exposed surface is unattainable. Therefore, although the near-surface grain density may differ between the three grain sizes (10, 20, and 40 $\upmu$m), they remain constant within the ten microstructures for each case. As such, the 10 and 40 $\upmu$m microstructures possess 10 and 2 GBs at the exposed surface. Representative microstructures for all three grain sizes are shown in Fig. 6, while the entire set of microstructures for the 10 and 40 $\upmu$m microstructures can be seen in Figs. S.4 and S.5 (Supplementary Information). The resultant weight loss and corrosion depth following a 500-hour simulation are given in Fig. 7. The SD data at 100, 250, and 500 hours exposure time can be seen in Fig. 8.

The weight loss shown in Fig. 7(a) showcases that the 10 $\upmu$m average grain size microstructures exhibits the most severe weight loss after 500 hours of exposure time. Alternatively, the largest grain size simulated of 40 $\upmu$m displays the highest corrosion resistance. Unsurprisingly, neither the 10 $\upmu$m nor the 40 $\upmu$m microstructures coincide with the experimental data extracted from a 20 $\upmu$m average grain size F/M specimen. Fig. 7(b) portrays that the average corrosion depth reached by liquid Li is largely independent of the grain size. Although there is little difference between the corrosion depths, there is nonetheless a noticeable pattern in the SDs for each grain size simulated, whereby the 10 $\upmu$m grain size possesses the smallest SDs and the 40 $\upmu$m microstructures give rise to the largest SD, more than twice as large as the SDs generated from the 20 $\upmu$m microstructures, after 500 hours (Fig. 8(b)). The SDs regarding the weight loss show an exceeding variation for the 10 $\upmu$m microstructures, whereas the 20 $\upmu$m and 40 $\upmu$m produce similar minute variances.

Fig. 9 highlight the IGC evolution across the representative microstructures (Fig. 6) for 30,000 hours of exposure time. Analogous to varying the near-surface grain density, the detrimental impact of Li corrosion is apparent, given by the extensive reach during the exposure time. The correlation between the average grain size and severity of corrosion is similarly made clear. The 40 $\upmu$m microstructure experiences considerably less degradation via the reduced presence of Li within the material. The low density of GBs, resulting in comparatively long GBs with little GB branching, limited the deterioration following exposure to liquid Li. Alternatively, the 10 $\upmu$m microstructure suffered extreme IGC over the exposure time, established through the vast area of the simulated domain whereby $\phi=0$. Furthermore, due to the relatively finer grain size, grains are more easily engulfed by Li containing GBs, leaving behind an intrigue network of corroded GBs with vastly diminished concentrations of Cr, Fig. 10.

2.3 Influence of the thickness of the smeared GB region

As detailed in Section 4.1, the computational GB thickness $l_{p}$ dictates the thickness of the smeared GB region that has increased diffusivity relative to the grains. In addition, implementation of the physical thickness of Cr depletion at the GB, $\delta_{gb}$, circumvents the issue with modelling nanometer-sized features far from their physical dimensions and thus avoids computationally expensive simulations. As such, $\delta_{gb}$ is kept constant while $l_{p}$ is altered to analyse its impact on the corrosion process of the simulated specimen. For this case, the representative microstructures for the 5, 6, and 7 GB microstructures (Fig. 1) were selected. Only a single run is completed for each change in $l_{p}$, as running all ten microstructures for each case is not needed to establish a correlation. Fig. 11 shows the weight loss and corrosion depth for the three microstructures and the values of $l_{p}$ chosen: 50, 100, and 200 nm.

For all three types of grain structures with 5, 6, and 7 GBs at the exposed surface, the 50 nm thickness gives rise to the lowest weight loss and 200 nm thickness results in the most extensive weight loss. That said, the number of GB entry points remains a deciding factor in weight loss. This is most clearly seen for the 7 GB microstructure with 100 nm thickness results in more significant weight loss than for the 5 GB microstructure with 200 nm thickness. This is additionally valid for the 7 GB microstructure with 50 nm thickness when compared against the 5 GB microstructure with 100 nm thickness. On average between all three microstructures simulated, decreasing the computational GB thickness from 100 nm to 50 nm, results in a final weight loss drop of 19.7$\%$ whereas when increased from 100 nm to 200 nm the final weight loss increases by 24.7$\%$. As presented in Fig. 2, there is no obvious relationship between the corrosion depth and the near-surface grain density. Consequentially, varying $l_{p}$ produces three distinctive groups as opposed to the weight loss data, which is clearly spread. It is clear that the microstructures with 50 nm thickness suffer the greatest corrosion depth, and those with 200 nm thickness have the highest resistance to Li penetration.

2.4 Influence of saturation

As stated in Section 4.2, a concentration sink is implemented to encapsulate the leaching effect of liquid Li in static conditions. This, however, omits the impact of saturation, which hinders and ultimately halts the corrosion process, leading the corrosion data to plateau [44]. To prove the capability of the current model to capture the effect of saturation and to observe whether a concentration sink might influence the corrosion data, a small liquid phase in contact with the solid phase is introduced whereby the concentration sink boundary condition is removed. The incorporation of a liquid phase necessitates a liquid phase diffusivity. That said, as the IGC process depends on the solid-state diffusivity, the value of the liquid phase diffusivity is not of significant importance as long as it is greater than the GB diffusivity. The liquid phase diffusivity of Cr in static liquid Li has been reported to be $2.2\times 10^{-8}$ m²/s at 570 $\degree$C [45]. However, to keep the diffusion coefficient gradient between phases at a manageable difference to avoid unnecessary computational cost, the liquid phase diffusivity is assigned to be $D_{\mathrm{gb}}$ as this will always be greater than $D^{\prime}_{\mathrm{gb}}$.

As the normalized concentration of Cr in the liquid phase ($\overline{c}$) approaches the normalized equilibrium concentration of the liquid phase ($\overline{c}_{\mathrm{l,eq}}$), the chemical free energy density decreases thereby reducing the chemical driving force of the corrosion process, Eqs. (4) and (5). The diminished chemical driving force decreases the rate of dissolution, thereby producing an asymptotic relationship between $\overline{c}$ and $\overline{c}_{\mathrm{l,eq}}$. As such, herein saturation of the system is considered met when $\overline{c}\geq 99\%\cdot\overline{c}_{\mathrm{l,eq}}$. It is consequentially found when simulating a 1 $\upmu$m thick liquid phase in contact with the reference microstructure from Section 4.2 (i.e., 20 $\upmu$m average grain size with 6 GBs at the exposed surface) that saturation is achieved after 6000 hours. This is shown in Fig. 12, where the normalized concentration of the liquid phase quickly approaches and eventually reaches the condition of saturation. The thickness of 1 $\upmu$m is chosen to reach saturation relatively quickly. Consequentially, all subsequent simulations involving the liquid phase are set with a 6000-hour exposure time with a 1 $\upmu$m thick liquid phase. To survey the influence of the liquid phase, and thus, the effect of saturation on the corrosion process, the liquid phase is introduced to the models described in Section 2.1 and Section 2.2.

The corrosion data for varying the near-surface grain density can be seen in Fig. 13. As all the microstructures with 5, 6, and 7 GBs are in contact with an identical quantity of liquid phase, the amount of dissolved Cr needed to reach saturation is constant, thereby producing, if not identical weight loss profiles resulting in largely consistent final weight loss of 2.22 g/m² following 6000 hours exposure time, Fig. 13(a). Nonetheless, the weight loss experienced by the 7 GB microstructures plateaus before the other two grain structures, indicating the system reached saturation quicker. Therefore, the corrosion process of leaching halted earlier. Alternatively, the 5 GB microstructures take noticeably longer for the weight loss to eventually plateau and terminate. Fig. 14(a) reinforces this notion, with the 7 GB microstructures clearly reaching saturation within a shorter exposure time than that of the 6 GB, followed by the 5 GB microstructures. Interestingly, contrary to the results displayed in Section 2.1 where a concentration sink is implemented, the SDs for the weight loss when saturation is considered are firstly too narrow to have any influential significance and secondly are all consistent across the three types of microstructural conditions simulated. The effect of saturation on the corrosion process is also seen in Fig. 13(b), where corrosion depth eventually plateaus. In this case, there is, although marginal, a clear trend emerging. The 5 GB microstructures produce the deepest depth, followed by the 6 GB, and finally, the 7 GB microstructures give rise to the shallowest corrosion depth. The SDs of the corrosion depths remain in the same order as that found with the concentration sink, where the 5 GB microstructures possess the greatest followed by 6 GB and finally 7 GB microstructures.

Additional analysis is carried out in which the 1 $\upmu$m thick liquid phase is implemented while varying the average grain size from 10, 20, and 40 $\upmu$m. Initially, the weight loss data, displayed in Fig. 15(a), shows that the grain size and rate of dissolution are positively correlated. Moreover, the 10 $\upmu$m microstructures reaches saturation first, indicative of the comparatively early plateau of the weight loss profile. With respect to the 40 $\upmu$m microstructures, the weight loss profile has yet to plateau, indicating the system did not reach saturation. This is confirmed by Fig. 14(b) whereby the finer grain-sized microstructures bring the system to saturation at the fastest rate. In addition, it is further evident from Fig. 14(b) that the 40 $\upmu$m microstructures did not reach saturation within the 6000 hours exposure time. Interestingly, however, is the progression of the weight loss profile for the 10 $\upmu$m microstructures, as it peaks around 3000 hours, where thereafter begins to dip and drop below the weight loss profile for the 20 $\upmu$m microstructures. As the same quantity of liquid phase is implemented throughout, the total weight loss experienced strongly resembles the final weight loss when varying the near-surface grain density. Another similarity between Fig. 13(a) and Fig. 15(a) is the negligible SDs. Fig. 15(b) displays a clear pattern in the corrosion depth when altering the average grain size of the microstructure. As depicted, the 40 $\upmu$m microstructures show the deepest corrosion depth, whereas the 10 $\upmu$m microstructures experiences the shallowest corrosion depth. Furthermore, similar to the weight loss data regarding the 10 $\upmu$m microstructures, the corrosion depth peaks and subsequently recedes. Analogous to when the concentration sink is implemented, the SDs for the corrosion depth remain in the same order starting at 10, 20, and 40 $\upmu$m in increasing order.

3 Discussion

It is clear that the most susceptible microstructures to weight loss possess the greatest number of near-surface grains giving rise to more GB entry points to infiltrate into the material, Fig. 2(a). To emphasize this trend, the inverse projected grain size ($G$), where $G$ is the ratio of the exposed surface length against the number of near-surface grains, was plotted against the total weight loss (not shown here for brevity). This yielded an $R^{2}$ of 0.9973 closely aligning with Bhave et al. [43] reported $R^{2}$ of 0.9966, consolidating the significance of GB entry points in IGC. The initial weight loss is dictated by the number of GB entry points, where thereafter, the total length of GBs available to corrode, governed by the grain size, offers a more comprehensive indicator to the degree of dissolution. Reducing the average grain size dramatically increases the density of GBs and thus amplifies the weight loss. Furthermore, the greater variation in the distribution of GBs among the 10 $\upmu$m microstructures, impacting the amount of accessible GBs, resulted in high variances in the weight loss, Fig. 8(a), thus behaving more inconsistently relative to the larger grain sizes. The inverse relationship between grain size and IGC susceptibility has been similarly observed experimentally [46, 47], highlighting the capabilities of the model.

As the corrosion process can continue indefinitely, via the concentration sink, the relationship between corrosion depth and microstructural features is not evident, Figs. 2(b) and 7(b). Nonetheless, it can offer insight into governing traits that dictate the susceptibility of structural materials to corrosion depth. The large variance in corrosion depths for the 40 $\upmu$m microstructures suggests that penetration depth is highly sensitive to the microstructure when the average grain size is large, Fig. 8(b). The relative orientation of GBs greatly effects the ability of Li to penetrate into the material, whereby the significance of their orientation increases with fewer GB entry points. Furthermore, the degree of GB branching governs the length Li must diffuse before changing its course, placing additional significance on their orientation as it strictly controls the depth Li can reach. Therefore, larger grain sizes produce a greater spread in corrosion depths as opposed to finer grain sizes. This suggest, in light of predictability, that finer grained structural materials are more desirable to optimize consistency in relation to IGC. That said, a holistic approach need consider the increased susceptibility of IGC as a consequence of the reduced grain size in order to determine the prime microstructure [48]. Likewise when varying the near-surface grain density; the 7 GB microstructures initially possess, albeit marginally, more entry points through which Li penetrates and, therefore, the relative orientation of the GBs has less influence on the depth of corrosion resulting in a lower variance, Fig. 3(b).

Replacement of the concentration sink with the effect of saturation in the corrosion process shifts microstructural influence on the corrosion behaviour. Due to the limitless corrosion process brought by the concentration sink, the relationship between the microstructural features (i.e., near-surface grain density and grain size) and weight loss is evident whereas the corrosion depth is independent of any microstructural changes. On the other hand, incorporating a minor quantity of liquid phase, thereby implementing saturation, clearly highlights the synergy between corrosion depth and microstructural features, yet consequentially dampens the correlation with weight loss. It should be noted, the conditions brought by the concentration sink strongly mimic the conditions of a dynamic breeder loop. The thermal gradients of dynamic loops produce detrimental corrosive conditions created by the perpetual dissolution and deposition of material; mass transfer, thus inhibiting saturation [12, 10]. Consequentially, evaluating the corrosion behaviour between static (i.e. saturation) and dynamic (i.e. concentration sink) conditions offers insight into the varied compatibility of F/M steels with liquid Li when transitioning from a controlled environment to higher fidelity conditions.

It is apparent that the near-surface grain density dictates the rate of saturation, as it offers more entry points to initiate penetration. Therefore, a greater number of GB entry points yields a faster rate of dissolution thus reaching saturation quicker, Fig. 14(a). This, in turn, controls the corrosion duration before the process is halted, ultimately dictating the depth Li is able to penetrate. The average grain size, via the bulk GB density, has greater influence on the rate of saturation, impacting more heavily the corrosion depth. Moreover, the 40 $\upmu$m microstructures had yet to reach saturation as seen through the sustained severity of the corrosion depth, Figs. 14(b) and 15(b). Regarding the variance in the corrosion behaviour, the introduction of saturation does not affect the order of corrosion depth SDs remaining as 7 GBs $<$ 6 GB $<$ 5 GB and 10 $\upmu$m $<$ 20 $\upmu$m $<$ 40 $\upmu$m. Nonetheless, the SDs on weight loss are significantly less prominent when saturation is accounted for. Overall, high density of near-surface grains coupled with small grain size maximizes the rate of saturation, thereby reducing the corrosion duration, and in turn, producing shallower corrosion depths, aligning with experimental reports [49]. It is additionally observed that the corrosion depth and the weight loss for the 10 $\upmu$m microstructures peak where thereafter it decreases, Fig. 15. As mentioned, near saturation the chemical free energy density decreases considerably, and for the 10 $\upmu$m microstructures, it becomes smaller than the interfacial free energy term, thus making the latter the dominant term in Eq. (8). Consequentially, the phase-field interface becomes more rigid and thus becomes inapt in simulating complex morphologies causing the IGC pathways to become narrower and shorter during the later stages of exposure, resulting in an apparent diminished weight loss and corrosion depth.

Figs. 4 and 9 illustrates the IGC evolution with respect to $\phi$, showcasing the diffusive behaviour of Li through the material. As the exposure time progresses, near-surface grains are engulfed by corroded GBs, yet the grains remain largely uncorroded. This becomes more apparent as the average grain size decreases. In reality, as corroded GBs surround a grain, they become susceptible to detachment from the bulk matrix, as observed experimentally [50, 51], producing a surface pebbled morphology. Although the corrosion depth was independent of microstructural features, it is overwhelming apparent that the degree of Li infiltration is strongly correlated to the near-surface and bulk GB density. The greater density of GBs additionally exposes more material to leaching, exacerbating the extent of Cr dissolution in the material, Figs. 5 and 10.

The value of $l_{p}$ governs the thickness of the Cr depletion region, thus dictating the carbide-rich region whereby Li can penetrate and leach out alloying elements influencing the weight loss. Naturally, a positive correlation arises between $l_{p}$ and weight loss, Fig. 11(a). Additionally, attributed to the constant product approach, the value of $l_{p}$ directly affects the GB diffusivity ($D^{\prime}_{\mathrm{gb}}$) ultimately altering the kinetics of the reaction. Consequentially, in increasing $l_{p}$, Li is unable to penetrate as deep into the material leading to a decreased average corrosion depth, Fig. 11(b). It should be stressed that the width of GB Cr depletion should not impact the kinetics of Li penetration; the observed correlation is purely artificial and attributed to the constant product approach. Nonetheless, the influence on weight loss is valid, which emphasizes the importance of regulating the degree of alloy segregation at the GBs to optimize the resistance of high Cr F/M steels to liquid Li. Moreover, slight kinks, resulting in a sudden jump in corrosion depth, are observed in Fig. 11(b).

Albeit displaying key microstructural features to improve the IGC susceptibility of steel alloys, there remain limitations that dampen the fidelity of the model. To align with the experimental work of Xu et al. [44], a computationally intensive liquid phase model would be required, necessitating a concentration sink instead. Although a valid assumption prior to saturation, it permits a limitless corrosion process ultimately resulting in unrealistic behaviours over longer exposure times. Furthermore, the addition of a 1 $\upmu$m liquid phase required 6000 h for the reference microstructure to saturate, Fig. 12, conflicting with experimental data. The chemical free energy density curvature parameter $A$ is selected based on the accuracy of the phase-field predictions to the experimental data. This proportionality constant should be determined and applied to liquid Li systems to improve the capability of the current framework. The microstructures employed are not fully representative of a tempered martensite microstructure due to the absence of GB angles, thus omitting low-angle GBs (i.e., lath and block boundaries) which should be included to fully capture the microstructural influence of IGC. In addition, the model precluded the nucleation of corrosion products, namely the ternary nitride complex, a prevalent mechanism in liquid Li corrosion with steels [52]. Its formation goes beyond perturbing the chemical potential of the system, thus exacerbating the corrosion process, yet additionally is thought to amplify the kinetics of penetration by fostering a favourable diffusive environment for Li [53].

Regarding the long-term performance of the material following exposure with liquid Li, confidence can be placed on the concentration sink model as this best resembles the continual dissolution conditions of a dynamic breeder blanket. Plus, the long-term behaviour of steels with respect to static systems is evidently dependent on the volume of liquid metal used. As detailed, the continual and indefinite corrosion process brought by the concentration sink results in largely identical corrosion depths independent of the grain structure of the specimen. As such, an average corrosion depth from the results displayed in Section 2.1 and Section 2.2, following 30,000 hours, equating to over 3 years of exposure, results in 65.8 $\upmu$m. This extensive penetration into the material would inevitably cause severe degradation of the mechanical integrity of the material thereby necessitating routine replacement of the material within this time period. Additionally, the current model assumes the primary degradation mechanism to be intergranular penetration of Li via solubility-driven dissolution. Although this is a governing process in the compatibility between liquid Li and F/M steels, there remains other phenomena not included in the presented model which would likely exacerbate the penetration rate and thus need be considered to assess the longevity of these steels. A few notable examples include the impact of non-metallic impurities thus forming corrosion products, the influence of fluid flow and the effect of dissimilar metallic systems. Compounding these processes in a comprehensive model to accurately evaluate the performance of structural materials with liquid metals is paramount in identifying the ideal candidate for future fusion reactors.

The clear relationship between the microstructure of the specimen and its susceptibility to intergranular penetration detailed through the model predictions grants perspective on the key traits to alter, thereby maximizing the performance of structural alloys in contact with liquid Li. Under an indefinite dissolution process, the grain size should be enlarged, considering desired mechanical properties, to limit the amount of available GBs to diffuse through and thus corrode. Greater attention should be placed on expanding the size of PAGs as the corresponding GBs, due to their high angle nature, foster greater density of Cr₂₃C₆ type carbides. Reports have detailed that higher austenitization temperatures over a longer duration facilitate greater growth of austenite grains [54, 55]. Another potential mitigation strategy is to limit the degree of Cr segregation to the GBs either by reducing the carbon content in the steel or by employing high-affinity carbide-forming elements such as vanadium, niobium, or tantalum [56, 17]. This would additionally diminish the presence of free carbon within the matrix of the steels, reducing the tendency for non-metallic impurity interaction with liquid Li.

To demonstrate the capability of the current model as well as compare the implication in modelling varying dimensions, the results for a 3D simulation following 500 hours exposure time are shown in Fig. 16. The geometry and computational conditions of the 3D model mimicked that of the reference 2D model detailed in Section 4.2, i.e., 100 $\upmu$m in length and width with a depth of 15 $\upmu$m. The average grain size of the microstructure was 20 $\upmu$m with a concentration sink implemented on the upper-most boundary positioned at $z$ = 15 $\upmu$m. The depth of the material was selected to minimize the computational expense and was guided through the final corrosion depth reached following 500 hours exposure time with a concentration sink from the 2D models equating to $\sim$12 $\upmu$m, Fig. 2 and Fig. 7. The evident evolution of the order parameter $\phi$ along the GBs can be seen in Fig. 16. In accordance to the 2D model, the enhanced corrosion is isolated solely at the boundaries that possesses the heightened diffusivity, enabling the rapid intergranular extraction of alloying elements. The final weight loss of the 3D microstructure following 500 hours exposure time equated in 80.93 g/m², far exceeding the 1.16 g/m² experienced for the 2D microstructures detailed in Section 4.2. It is important to emphasize that in the current formulation the interface kinetics coefficient $L$ is tailored to 2D microstructures. Therefore, $L$ ought to be derived relative to material properties to generate a more direct linkage between model predictions and experimental data. Bhave et al. [43] similarly reported greater mass loss from the 3D models relative to the 2D models in predicting the IGC behaviour of Ni-Cr alloys exposed to molten Li salt. Other similar studies [40, 41, 42, 57, 58, 59] have been primarily focused on proof-of-concept of their respective formulations rather than analysing the significance regarding the corrosion data between 2D and 3D simulations. Nonetheless, 2D models have shown to yield representative data, that supports experimental theories, thus formulating our understanding of the microstructural influence on IGC. Future work should consider the role of the flow of the liquid phase, the nucleation of corrosion products from non-metallic impurities, the crystallographic orientation of grains, and GB angle (high and low-angle GBs), in mass transport and IGC of polycrystalline steels. It would additionally be warranted to conduct a parametric study consisting of a larger spectrum of microstructural dimensions (i.e., near-surface grain density and grain size) thus extracting more comprehensive data allowing deeper analysis of their impact on the IGC process.

In conclusion, a numerical framework based on the phase-field method is developed for assessing the IGC of polycrystalline steels with corrosive media. Herein, the model was applied to capture the leaching effects of liquid Li when in contact with a 9 wt% Cr F/M steel. The key findings can be summarized as follows:

1. The governing microstructural features that dominate and ultimately dictate the severity of corrosion are the number of GB entry points, i.e., near-surface grain density. The grain boundary density in the bulk, i.e., grain size, governs the susceptibility to Li corrosion.

2. The thickness of Cr depletion along the GBs $l_{p}$ plays a deciding role in the weight loss of the simulated material. This emphasizes the importance of limiting alloying GB segregation and consequential depletion in these regions, ideally hindering the ability for Li to diffuse intergranularly.

3. The effect of saturation is considered and compared against the use of a concentration sink. The former highlights the dependence of the corrosion depth on the microstructural features, whereas the latter amplifies the relationship with weight loss. Although the 10 $\upmu$m microstructures exhibit the fastest weight loss, this in turn results in the shallowest corrosion depth. With the idea of using grain engineering to maximize the compatibility with Li, it is worth dictating which corrosion behaviour is more detrimental to the structural integrity of the alloy.

4. The SDs regarding the corrosion depths shed light on the importance of the orientation of the GBs as well as their total length (i.e., GB density) in facilitating the penetration of Li into the material. It is shown that the 40 $\upmu$m microstructures yield varying corrosion depths as a result.

4 Methods

4.1 Phase-field model of intergranular corrosion

It is assumed in the current work the primary degradation mechanism is driven by the bulk diffusion of Cr in the metal phases. This assumptions aligns with [52, 60, 61]. The impact of other phases present in the metal (i.e., precipitates and inter-metallic species) and the effect of the Li₂C₂ by-product (Eq. 1) on the corrosion process is not considered in this study. The following kinematic variables are introduced to characterize the two phases. A continuous phase-field variable $\phi(\mathbf{x},t)$ is implemented to track each phase and the evolution of the corroding interface. $\phi=1$ represents the solid phase (i.e., F/M steel), $\phi=0$ defines the liquid corrosive agent (i.e., liquid Li), and $0<\phi<1$ describes the thin diffuse interface separating the opposing phases, Fig. 17. As the composition of Cr in the F/M steels plays arguably the most central role in the corrosion process, modelling the entire composition of the steel is unwarranted. As such, the current model adopts a simplified F/M steel in the form of a simple binary Fe-Cr alloy. The metal phase in the present model is assumed to be Fe-9 wt% Cr steel. It is further assumed in the present investigation that the material has a uniform composition and is composed of equiaxed grains separated by GBs. The composition of the material and its evolution during corrosion is characterized by the normalized concentration of Cr $\overline{c}(\mathbf{x},t)$ = $c/c_{\mathrm{solid}}$, where $c$ is the absolute concentration of Cr and $c_{\mathrm{solid}}$ the concentration of Cr initially present in the material. Two independent diffusion coefficients between the GB ($D_{\mathrm{gb}}$) and metal grains ($D_{\mathrm{mg}}$) are introduced to enhance the corrosion process along GBs, Fig. 17 whereby the employed values are displayed Table 1. The GB interpolation is achieved via an additional stationary parameter $\eta(\mathbf{x})$, that takes value $\eta=1$ at the GBs and $\eta=0$ elsewhere, as discussed below in Section 4.1. In the present work, the grains and GBs maintain an isotropic diffusion behaviour. The contribution of crystallographic orientation and GB angle (high and low-angle GBs) will be addressed in future work.

For the underlying corrosion mechanism considered, the total free energy functional of the heterogeneous system in Fig. 17 can be expressed as

$$\mathcal{F}=\int_{\Omega}[f^{\mathrm{chem}}(\overline{c},\phi)+f^{\mathrm{grad}}(\nabla\phi)]d\Omega,$$ (2)

where $f^{\mathrm{chem}}$ and $f^{\mathrm{grad}}$ are the chemical and interfacial free energy densities detailed further below. $\Omega$ in the previous expression represents the whole system domain that includes both the corrosive liquid Li environment and the polycrystalline material.

The chemical free energy density is expressed as the weighted sum of free energy density from each contributing phase [62]

$$f^{\mathrm{chem}}(\overline{c},\phi)=h(\phi)f^{\mathrm{chem}}_{\mathrm{s}}(\overline{c}_{\mathrm{s}})+(1-h(\phi))f^{\mathrm{chem}}_{\mathrm{l}}(\overline{c}_{\mathrm{l}})+\omega g(\phi),$$ (3)

where $f^{\mathrm{chem}}_{\mathrm{s}}(\overline{c}_{\mathrm{s}})$ and $f^{\mathrm{chem}}_{\mathrm{l}}(\overline{c}_{\mathrm{l}})$ are the chemical free energy densities with respect to the normalized concentrations in the liquid ($\overline{c}_{\mathrm{l}}$) and solid ($\overline{c}_{\mathrm{s}})$ phases. $g(\phi)=16\phi^{2}(1-\phi)^{2}$ is the double-well free energy function employed to describe the two equilibrium states for the solid ($\phi=1$) and the liquid ($\phi=0$) phases. $\omega$ is the constant that determines the height of the energy barrier at $\phi=1/2$ between the two minima at $\phi=0$ and $\phi=1$. $h(\phi)=\phi^{3}(6\phi^{2}-15\phi+10)$ is a monotonously increasing interpolation function that interpolates the chemical free energy density between the two phases. The chemical free energy density of each phase can be reasonably approximated by a simple parabolic function [63] around equilibrium concentrations with the same free energy density curvature parameter $A$ as [30]

$$f^{\mathrm{chem}}_{\mathrm{s}}(\overline{c}_{\mathrm{s}})=\frac{1}{2}A(\overline{c}_{\mathrm{s}}-\overline{c}_{\mathrm{s,eq}})^{2}\quad\quad f^{\mathrm{chem}}_{\mathrm{l}}(\overline{c}_{\mathrm{l}})=\frac{1}{2}A(\overline{c}_{\mathrm{l}}-\overline{c}_{\mathrm{l,eq}})^{2},$$ (4)

where $\overline{c}_{\mathrm{s,eq}}=c_{\mathrm{solid}}/c_{\mathrm{solid}}=1$ and $\overline{c}_{\mathrm{l,eq}}=c_{\mathrm{sat}}/c_{\mathrm{solid}}$ are the normalized equilibrium concentrations for the solid and liquid phase. Here, $c_{\mathrm{sat}}$ represents the saturation limit of the metal species in the liquid phase. Each material point in the present model is characterized as a mixture of both solid and liquid phases with different compositions yet the same diffusion chemical potentials [64]. This assumption renders the following expression for the chemical free energy density $f^{\mathrm{chem}}(\overline{c},\phi)$ [65]

$$f^{\mathrm{chem}}(\overline{c},\phi)=\frac{1}{2}A[\overline{c}-h(\phi)(\overline{c}_{\mathrm{s,eq}}-\overline{c}_{\mathrm{l,eq}})-\overline{c}_{\mathrm{l,eq}}]^{2}+\omega g(\phi).$$ (5)

The interfacial free energy density is commonly expressed as

$$f^{\mathrm{grad}}(\nabla\phi)=\frac{1}{2}\kappa|\nabla\phi|^{2}$$ (6)

where $\kappa$ is the isotropic gradient energy coefficient. The phase-field parameters $\omega$ in Eq. (3) and $\kappa$ in Eq. (6) are connected to the interfacial energy $\Gamma$ and the chosen nominal interface thickness $\ell$ as [66]

$$\kappa=\frac{3}{2}\Gamma\ell\quad\quad\omega=\frac{3\Gamma}{4\ell}.$$ (7)

The governing equations for the independent kinematic fields $\phi(\mathbf{x},t)$ and $\overline{c}(\mathbf{x},t)$ are derived by minimizing the total energy of the system and conserving the total concentration of Cr composition within the system [67]. The evolution of the non-conserved phase-field parameter $\phi$ follows the Allen-Cahn equation [68]

$$\frac{\partial\phi}{\partial t}=-L\frac{\delta\mathcal{F}}{\delta\phi}=-L\Big{(}\frac{\partial f^{\mathrm{chem}}}{\partial\phi}-\kappa\nabla^{2}\phi\Big{)}\quad\text{in}\quad\Omega\qquad\kappa\mathbf{n}\cdot\nabla\phi=0\quad\text{on}\quad\partial\Omega,$$ (8)

where $L$ is the kinetic coefficient that characterizes the interfacial mobility and controls the motion of the solid$-$liquid interface. The magnitude of this parameter determines the underlying corrosion mechanisms, such as activation-controlled and diffusion-controlled, and regulates the corrosion rate [29, 30, 31, 37, 69]. The condition for $L$ to capture both mechanisms is given in the following section.

The transport of Cr composition in the system is subjected to the conservation law

$$\frac{\partial\overline{c}}{\partial t}=-\nabla\cdot\mathbf{J}\quad\mathrm{in}\;\Omega\qquad\mathbf{J}=-M\nabla\Big{(}\frac{\delta\mathcal{F}}{\delta\overline{c}}\Big{)}\qquad\mathbf{n}\cdot\mathbf{J}=0\quad\text{on}\quad\partial\Omega,$$ (9)

where $\mathbf{J}$ stands for the diffusional flux and $M$ the mobility parameter that characterizes the motion of Cr composition. Here, the mobility parameter is expressed as: $M=D/(\partial^{2}f^{\mathrm{chem}}/\partial\overline{c}^{2})=D/A$ [30], where $D$ is the effective diffusion coefficient of Cr composition. After substituting the expression for the mobility parameter into Eq. (9), the diffusional flux is written as

$$\mathbf{J}=-D\nabla\overline{c}-Dh^{\prime}(\phi)(\overline{c}_{\mathrm{l,eq}}-\overline{c}_{\mathrm{s,eq}})\nabla\phi.$$ (10)

The effective diffusion coefficient of Cr composition is interpolated between the grain bulk and GBs as

$$D=D_{\mathrm{gb}}^{\prime}\eta+(1-\eta)D_{\mathrm{mg}}\quad\quad D^{\prime}_{\mathrm{gb}}=\frac{\delta_{gb}}{l_{p}}D_{\mathrm{gb}},$$ (11)

where $D_{\mathrm{mg}}$ and $D_{\mathrm{gb}}$ are the diffusion coefficient of Cr composition in the grain bulk and along GBs, respectively. $D_{\mathrm{gb}}\gg D_{\mathrm{mg}}$ is enforced to enhance the diffusion of Cr composition along GBs. $D^{\prime}_{\mathrm{gb}}$ in the previous equation is the effective diffusion coefficient of Cr composition along GBs defined using the constant product approach [70]. This approach relates the actual diffusion coefficient $D_{\mathrm{gb}}$ and the ratio between the physical $\delta_{gb}$ and computational thickness $l_{p}$ of the Cr depletion region along GBs. The expression proportionally alters the GB diffusivity with respect to an experimentally determined GB thickness, which in turn, rectifies the relative contribution of the otherwise unaltered GB diffusivity and additionally decreases the computational expense of simulating nanometre-sized features. The governing equation for the interpolating parameter $\eta$ is given as [71]

$$\nabla\cdot(-l_{p}^{2}\nabla\eta)+\eta=0\quad\text{in}\quad\Omega\qquad l_{p}^{2}\mathbf{n}\cdot\nabla\eta=0\quad\text{on}\quad\partial\Omega\quad\text{and}\quad\eta=1\quad\text{on GBs}.$$ (12)

GBs are given a value of $\eta=1$ and smoothly transitions to $\eta\rightarrow 0$ further away from the GBs, Fig. 17. The sensitivity of the evolution of IGC to the computational thickness of the smeared GB region is investigated in Section 2.3. It is assumed in the present formulation that the Cr GB depletion region is given a uniform thickness, neglecting any confinement effects.

The framework developed is implemented into the multi-purpose finite element software package COMSOL Multiphysics [72]. The computational domain is discretized using triangular finite elements with second-order Lagrangian interpolation functions. All regions expected to corrode in the metal grains are given a characteristic maximum element size at least five times smaller than the interfacial thickness $\ell$ to ensure a smooth transition between the metal and corrosive agent. This setting proved sufficient based on previous literature [30, 65, 67]. Moreover, as the evolution of the interface is expected to be most prominent at the GBs compared to the metal grain, a maximum element size of $\ell$/20 is applied to all GBs. Finally, to limit the computational cost, the remaining domain of the solid phase is given a maximum element size of $\ell$. These conditions are fulfilled in all the simulations. The depictions of the finite element mesh for a representative case study are shown in Fig. 18(c) and Fig. 18(d). Each simulation consists of a two-step study. The governing equation (12) for the interpolating parameter $\eta$ that defines the smeared GB thickness is solved in the first step using a steady-state (time-independent) solver. The governing equations (8) and (9) for the evolution of the phase-field parameter $\phi$ and Cr composition $\overline{c}$ are then solved in the second step using a time-dependent study. An implicit time-stepping method is used for temporal discretization in the time-dependent step. A fully coupled solution algorithm is selected to solve the governing equations. The maximum time step is 2 hours. The solver accuracy in each time step is controlled by a relative tolerance of 10^-4. The code developed together with example case studies and documentation are available at https://mechmat.web.ox.ac.uk/codes (after article acceptance).

4.2 Model calibration and validation

The model developed is calibrated and validated against experimental data given in [44] where greater detail regarding the experimental set-up can be found. The steel specimen underwent a normalized and tempered heat treatment resulting in a tempered martensite microstructure with a prior austenitic grain (PAG) size of 20 $\upmu$m [20]. The experimental apparatus used by Xu et al. is displayed in Fig. 18. 100 mL of high purity Li was used, whereby the concentration of nitrogen was estimated to never surpass 100 ppm. Nitrogen has been experimentally reported to exacerbate the corrosion process through chemical reactions producing the stable ternary nitride corrosion product, thereby intensifying the leaching and degradation of structural materials [52, 73]. Nonetheless, following the equilibrium concentration of nitrogen required to form the ternary nitride corrosion complex [74], the propensity for the corrosion complex to have nucleated and thus contribute to the penetration process was unlikely, consolidating the bulk diffusion mechanism adopted herein. The geometric ratio between the volume of Li to the total surface area of the specimen was $\sim$ 4 cm. The specimen was exposed to static liquid Li at 600 $\degree$C for 750 hours. Experimental measurements in terms of weight loss and corrosion depth are used to calibrate the model.

The numerical simulation is conducted by considering a 2D computational domain with 100 $\upmu$m in length and height. Only 2D simulations are performed for simplicity whereby the implications of this assumption are discussed later. The dimensions are selected based on minimizing computational cost while maintaining an adequate region of microstructural features to observe its influence. The computational domain is characterized by a microstructure with an average grain size equal to the PAG size, 20 $\upmu$m from the experiment in [44], Fig. 18. Solely simulating prior austenitic GBs is a reasonable approximation based on literature that details the extensive corrosion and significant quantity of liquid Li at high-angle GBs following compatibility experiments [73, 75]. Nonetheless, incorporating the entire tempered martensitic constituents (i.e., lath and block boundaries) is important to fully capture the microstructural effects of liquid Li IGC and, as such, this will be investigated in the future. Ten different microstructures are constructed to have a statistically significant sample and generate an average. The ten microstructures used can be viewed in Fig. S.1 (Supplementary Information). The simulation starts by solving for the smeared GB thickness using the governing equation and corresponding boundary conditions in Eq. (12). $\eta=1$ is enforced on selected GBs that are expected to corrode. Afterward, the resulting set of the governing equations (8) and (9) is solved with accompanying initial and boundary conditions, Fig. 18. The computational domain initially consists solely of the metal phase. Thus, the initial values for the phase-field variable and Cr concentration are set to $\phi=1$ and $\bar{c}=1$. The liquid Li phase is represented through a concentration sink implemented by prescribing the Cr concentration to zero. The sink is positioned on the uppermost boundary of the simulated domain. As the simulated solid phase is only exposed from the upper boundary, the other three edges possess zero flux boundary conditions ($\mathbf{n}\cdot\nabla\phi=0$ and $\mathbf{n}\cdot\mathbf{J}=0$) for both phase-field variable and Cr concentration. Substituting the liquid phase with a concentration sink does discount the effect of saturation, which ultimately halts the corrosion process in static conditions [44]. Consequentially, using a concentration sink to approximate the liquid phase facilitates endless IGC. Yet, it is a justifiable approximation before saturation is reached.

The material properties and model parameters employed in the simulation are summarized in Table 1. The diffusion coefficients of Cr in the metal grain bulk $D_{\mathrm{mg}}$ and along GBs $D_{\mathrm{gb}}$ are taken from Čermák et al. [76]. Due to a lack of experimental work, the interfacial energy $\Gamma$ is obtained based on GB interfacial energies for BCC metals and alloys [77, 78]. A value of $\Gamma$ = 4 J/m² is utilized in the present work. The interfacial thickness $\ell$ is selected to be significantly smaller than the computational domain and characteristic grain size ($\ell$ = 4 $\upmu$m). The phase-field parameters $\omega$ (Eq. (5)) and $\kappa$ (Eq. (6)) are then expressed using $\ell$ and $\Gamma$, as indicated in Eq. (7). The thickness of the Cr-depleted region $\delta_{gb}$ is adopted based on experimental work [79], which found that the thickness varied between 10$-$15 nm. As such, a conservative approach is taken in the present study, setting $\delta_{gb}$ = 15 nm. The computational GB thickness $l_{p}$ is set to $l_{p}$ = 100 nm. The two aforementioned thicknesses; $\delta_{gb}$ represents the experimentally determined parameter and $l_{p}$ signifies the thickness of the simulated GB region, are employed via the constant product approach, to link the artificial microstructure to physical material properties. The chemical free energy density curvature parameter $A$ in Eq. (5) is chosen from similar studies [30, 65, 67, 29]. Employing the first-principles calculations and the CALPHAD approach can yield a more accurate description of the free energy functional with thermodynamic properties over various compositions and temperatures [80], which can be fed into the current phase-field framework. The chemical driving force in the present model originates from the difference in Cr composition and the equilibrium values $c_{\mathrm{solid}}$ and $c_{\mathrm{sat}}$ in each phase, Eq. (5). The former represents the initial known concentration of Cr in a steel specimen before exposure to liquid Li, $c_{\mathrm{solid}}$ = 13.4 mol/L [44]. The latter signifies the equilibrium concentration of Cr in the steel specimen following exposure to liquid Li under static conditions. Due to the extended exposure time tested in the compatibility experiments conducted by Xu et al. [44], along with the knowledge that chemical equilibrium in static conditions with 100 mL of liquid Li is reached within 250 h [19, 44], it is warranted to take the post-corrosion test surface concentration of Cr as the equilibrium value, $c_{\mathrm{sat}}$ = 10.3 mol/L [20]. The interface kinetics coefficient $L$ governs the motion of the metal$-$liquid interface and, as such, dictates the rate-limiting step of the corrosion process. Comparatively large values give rise to diffusion-controlled corrosion and transitions to activation-controlled as the value decreases [29]. Dimensional analysis [67] indicates that $L$ needs to be $L\gg D_{\mathrm{gb}}/(\ell^{2}\omega)$ for diffusion-controlled and $L\ll D_{\mathrm{gb}}/(\ell^{2}\omega)$ for activation-controlled regimes. To ensure the corrosion process possesses adequate kinetics to operate under diffusion control, it is given a value of $L=$ 1 m²/(N$\cdot$s), which satisfies the previous condition.

The weight loss of the simulated specimen is determined in the simulation as

$$\Delta m=\int_{0}^{\overline{c}_{\mathrm{l,eq}}}\overline{c}\cdot c_{\mathrm{solid}}\cdot A_{Cr}\,d\Omega,$$ (13)

where $\Delta m$ represents the mass loss in grams and $A_{Cr}$ is the atomic weight of Cr in grams/mole. In accordance with the experimental procedure in evaluating the susceptibility of specimens with liquid Li, the obtained weight loss is normalized with respect to the exposed surface area of the specimen. As such, Eq. (13) is divided by the length of the concentration sink region, 100 $\upmu$m. The corrosion depth is determined by taking the mean corrosion depth from all advancing corrosion fronts ($\phi<0.5$) along GBs that are diffusing deeper towards the bulk. Fig. 19 illustrates the weight loss and corrosion depth averaged from ten varying microstructures. Fig. 19(a) showcases the strong resemblance between the phase-field predictions and experimental data regarding weight loss. The weight loss predicted by the model passes through the upper relative error at 100 hours and matches the experimental data at 250 hours with high accuracy. Moreover, the corrosion predictions displayed in Fig. 19(b) similarly shows a high correlation to the experimental data, resulting in a corrosion depth of 9.4 $\upmu$m after 250 hours. After 500 hours, the model predicted an average weight loss of 1.16 g/m² where the liquid Li reached an average depth of 12.5 $\upmu$m. The logarithmically varying profile for the weight loss and corrosion depth signifies that the corrosion process is in the diffusion-controlled region, further aligning with the literature and proving that the phase-field mobility parameter $L$ is adequately selected. The close comparison between the phase-field predictions and experimental measurements [44] showcases the capabilities of the current model to simulate IGC. Following model calibration and validation, microstructural and GB properties are altered to observe their dependence on the corrosion resistance of the simulated F/M specimen to static liquid Li in Section 2. The data is compared against the reference geometry described in this section.

Data Availability

The relevant data are available from the corresponding authors upon reasonable request.

Code Availability

The code developed together with example case studies and documentation is available at https://mechmat.web.ox.ac.uk/codes.

Author Contributions

A.L.: Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. S.K.: Writing – review & editing, Supervision, Software, Methodology, Investigation, Formal analysis, Conceptualization. D.N.-M.: Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization. J.L.: Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization. E.M.-P.: Writing – review & editing, Supervision, Software, Resources, Project administration, Methodology, Funding acquisition, Conceptualization. M.W.: Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization. All authors have read and approved the manuscript.

Competing Interests

The authors declare no competing interests.

Figure 1: Representative microstructures for the 5, 6 and 7 surface GB microstructure simulated. Each microstructure has an average grain size of 20 $\upmu$m, while varying the number of GBs at the exposed surface located at the upper most boundary.

Figure 2: Phase-field predictions following a 500-hour simulation for 20 $\upmu$m average grain size microstructure with 5, 6, and 7 GBs present at the exposed surface compared to experimental measurements [44]. Phase-field predictions for a weight loss and b corrosion depth. The data presented for each case is taken from an average of ten microstructures where the shaded coloured regions represents the SD of the ten simulations.

Figure 3: Standard deviations at 100, 250, and 500 hours intervals for 20 $\upmu$m average grain size microstructure with 5, 6, and 7 GBs present at the exposed surface. Standard deviations for a weight loss and b corrosion depth. The data presented for each case is taken from an average of ten microstructures.

Figure 4: Long-term IGC evolution while varying the near-surface grain density. The phase-field variable $\phi$ in the simulated metallic specimen during a 30,000-hour simulation for the three representative microstructures possessing 5, 6, and 7 GBs at the exposed surface.

Figure 5: Long-term intergranular leaching of Cr while varying the near-surface grain density. The normalized concentration of Cr $\overline{c}$ in the simulated metallic specimen during a 30,000-hour simulation for the three representative microstructures possessing 5, 6, and 7 GBs at the exposed surface.

Figure 6: Representative microstructures for the three average grain sizes simulated. The GBs present at the exposed surface are 10, 6, and 2 for the 10, 20, and 40 $\upmu$m microstructures.

Figure 7: Phase-field predictions following a 500-hour simulation for 10, 20, and 40 $\upmu$m average grain size microstructures compared to experimental measurements [44]. Phase-field predictions for a weight loss and b corrosion depth. The data presented for each case is taken from an average of ten microstructures where the shaded coloured region represents the SD of the ten simulations.

Figure 8: Standard deviations at 100, 250, and 500 hours intervals for 10, 20, and 40 $\upmu$m average grain size microstructures. Standard deviations for a weight loss and b corrosion depth. The data presented for each case is taken from an average of ten microstructures.

Figure 9: Long-term IGC evolution while varying the grain size. The phase-field variable $\phi$ in the simulated metallic specimen during a 30,000-hour simulation for the three representative microstructures with average grain sizes of 10, 20, and 40 $\upmu$m.

Figure 10: Long-term intergranular leaching of Cr while varying the grain size. The normalized concentration of Cr $\overline{c}$ in the simulated metallic specimen during a 30,000-hour simulation for the three representative microstructures with average grain sizes of 10, 20, and 40 $\upmu$m.

Figure 11: Phase-field predictions following a 500-hour simulation for 20 $\upmu$m average grain size geometry with 5, 6, and 7 GBs present at the exposed surface while varying $l_{p}$ by 50, 100 and 200 nm compared to experimental measurements [44]. Phase-field predictions for a weight loss and b corrosion depth. The coloured lines in the legend denote the microstructure (i.e., number of near-surface GBs) whereas the black lines of varying line-type illustrate the computational GB thickness $l_{p}$. In conjunction they signify the conditions of each dataset.

Figure 12: Normalized concentration of Cr in the liquid phase as a function of exposure time. With a 1 $\upmu$m thick liquid phase above a 20 $\upmu$m average grain size microstructure with 6 GBs, it takes 6000 hours for the system to reach saturation. The data presented here is an average from ten microstructures where the shaded coloured region represents the SD of the ten simulations. Note, the consistency across the ten microstructures resulted in a negligible SD. Sub-graph displays a magnified version to clearly highlight the liquid concentration nearing saturation. All subsequent concentration plots possess a range and domain akin to the sub-graph.

Figure 13: Phase-field predictions following a 6000-hour simulation for 20 $\upmu$m average grain size microstructure with 5, 6, and 7 GBs present at the exposed surface with a 1 $\upmu$m liquid phase in contact with the solid phase. Phase-field predictions for a weight loss and b corrosion depth. The data presented for each case is taken from an average of ten microstructures each where the shaded coloured regions represents the SD of the ten simulations. Sub-graph displays a magnified version to clearly distinguish between the weight loss profiles.

Figure 14: Normalized concentration of Cr in the liquid phase as a function of exposure time while varying microstructural properties. Normalized concentration of Cr in the liquid while varying a the near-surface grain density for a 20 $\upmu$m average grain size microstructure and b the average grain size. The data presented for each case is an average of ten microstructures where the shaded coloured regions represents the SD of the ten simulations.

Figure 15: Phase-field predictions following a 6000-hour simulation for 10, 20, and 40 $\upmu$m average grain size with a 1 $\upmu$m liquid phase in contact with the solid phase. Phase-field predictions for a weight loss and b corrosion depth. The data presented for each case is taken from an average of ten microstructures each where the shaded coloured regions represents the SD of the ten simulations.

Figure 16: 3D simulation of IGC. Depiction of the a geometrical dimensions of the 3D model and b the intergranular corrosion evolution at 25, 100, 250 and 500 hours exposure time. Only the solid phase $\phi\geq 1/2$ is shown for clarity.

Figure 17: Polycrystalline material in contact with a corrosive environment highlighting the diffuse interface between the liquid ($\phi=0$) and solid ($\phi=1$) phases. The GBs possess a heightened diffusivity $D_{\mathrm{gb}}$ compared to the metal grain $D_{\mathrm{mg}}$ via an additional parameter that differentiates between grain boundary ($\eta=1$) and metal grain ($\eta=0$). The corrosion mechanism is based on the bulk diffusion of Cr towards the exposed surface.

Figure 18: Computational conditions, domain and meshing. Schematic of the a experimental apparatus used by Xu et al. [44] to expose F/M specimen to static liquid Li, b corresponding computational domain consisting of a polycrystalline material with the initial values and boundary conditions, c finite element mesh of the whole computational domain, and d enlarged area corresponding to the red square in c highlighting the finite element size within the grains and along the GBs in the expected area of interface propagation. The exposed surface is the upper most boundary.

Figure 19: Phase-field predictions following a 500-hour simulation with a 20 $\upmu$m average grain size microstructure compared to experimental measurements from Xu et al. [44]. Phase-field predictions for a weight loss and b corrosion depth. The data presented is taken from an average of ten microstructures where the shaded coloured region represents the SD of the ten simulations.

Table 1: Material and computational parameters employed in the model.

Parameter	Value	Unit	Ref
Chemical free energy density curvature parameter $A$	$5\times 10^{9}$	N/m2
Interface kinetics coefficient $L$	1	m2/(N$\cdot$s)
Concentration in the solid phase $c_{\mathrm{solid}}$	13.4	mol/L	[44]
Saturated concentration in the liquid phase $c_{\mathrm{sat}}$	10.3	mol/L	[20]
Interfacial energy $\Gamma$	4	N/m	[78, 77]
Computational GB thickness $l_{p}$	100	nm
Physical Cr depletion thickness $\delta_{gb}$	15	nm	[79]
Interfacial thickness $\ell$	4	$\upmu$m
Grain boundary diffusivity $D_{\mathrm{gb}}$	$1.70\times 10^{-15}$	m2/s	[76]
Metal grain diffusivity $D_{\mathrm{mg}}$	$5.11\times 10^{-21}$	m2/s	[76]