Development of a numerical model for simulating stress corrosion cracking in spent nuclear fuel canisters

Prediction and detection of the chloride-induced stress corrosion cracking (CISCC) in Type 304 stainless steel spent nuclear fuel canisters are vital for the lifetime extension of dry storage canisters. This paper conducts a critical review that focuses on the numerical modeling and simulation on the research progress of the CISCC. The numerical models emphasizing the residual stress, susceptible microstructure, and corrosive environment are summarized individually. Meanwhile, the simulation studies on the role of hydrogen-assisted cracking are reviewed. Finally, a multi-physical numerical model, which combines the different fields is proposed based on our recent investigation.


INTRODUCTION
Generally, after being cooled in the water pool for a few years, spent nuclear fuel (SNF) would be removed and transferred to other sites for longer-term storage. The helium-filled stainless steel canisters, which reside in passively ventilated dry storage systems, are the current solution for resolving this issue. These dry storage systems were originally licensed for up to 20 years, and renewals also up to 20 years. In 2011, 10 CFR 72.42(a) was modified to allow for extending the initial license periods and renewals for up to 40 years. However, as the Department of Energy (DOE) in United States noted, before the final disposal of the SNF, the fuel may need to remain in storage for many decades beyond the original design intent. Before the extension of the interim dry storage facilities for long-term storage, the safety concerns must be identified and the potential degradation of these SNF canisters should be assessed. As highlighted in several recent reports 1-4 , chloride-induced stress corrosion cracking (CISCC) is considered as the major knowledge gap for the long-term performance of these dry storage systems.
CISCC is a type of degradation that leads to cracks in certain stainless steel materials 5 . Generally, three necessary conditions, i.e., susceptible metal material, aggressive environment and sufficient amount of tensile stress, must be presented before the CISCC can occur. As for those near-marine SNF dry storage canisters, analyses have shown that chloride-rich salts could be present on the canister surface [6][7][8] . Thus, when the canister surface temperature cools down, the salts would deliquesce to form a chloride-rich aqueous brine layer, which could offer the corrosive environment for the CISCC to happen. Meanwhile, the welded interim storage canisters are made of austenitic stainless steels, in which the Types 304/304L stainless steel (SS) and 316/ 316L SS are usually container material candidates, and Types 308 SS and 316 SS are weld metals for Types 304 SS and 316 SS, respectively. Previous studies have shown that austenitic stainless steels are susceptible to SCC when exposed to chlorides in marine atmospheres [9][10][11] . The overall susceptibility of the material is a function of many factors, including element composition and gradient, microstructure evolution during cold work and welding process, processing parameters and surface conditions. Tensile stress is also essential for the CISCC to happen in the SNF canisters. Though for dry storage canisters, external loads usually contribute small amount of stress to the metal, which is unlikely to provide sufficient driving force for the initiation and propagation of crack, the welding process during the construction of canister could result in high through-wall tensile stress to support CISCC 12,13 . Therefore, it is reasonable to expect that during the period of interim storage in near-marine environments, the three conditions would be fulfilled, at least at some Independent Spent Fuel Storage Installations sites, especially if we have to delay the development of a repository for final disposal.
CISCC of a dry storage canister would affect its safety functions. To investigate the potential of CISCC in stainless steel canister, many experimental studies have been conducted. Parts of these studies focus on the actual environmental conditions of the in-service canisters, like the surface temperature variables [14][15][16] , the relative humidity 17,18 , the sea salts deliquesce and concentration 16,[19][20][21] , and the residual stress for the canister structure 12 . Meanwhile, some CISCC related experimental tests on the canister materials based on the environmental conditions were conducted in the laboratory 10,[22][23][24][25][26][27] . Even though these experimental studies could offer important information directly or indirectly relevant to the CISCC behavior of in-service canisters, the results can only partially cover the actual in-service conditions. Also, because the corrosion in most cases happens at a very slow rate, the experimental tests usually need several weeks, or months to get some distinguishable data even the accelerated conditions were used. The high labor and financing cost make it hard to carry out any systematic experimental study. Besides, some of the underlying mechanisms about the initiation and growth of the potential pits and cracks related to the CISCC of actual canisters have still not been revealed.
In recent years, in addition to the experimental investigations, numerical methods that are capable of simulating the factors influencing the corrosion processes of SNF canisters and thus reliably predicting the service life have gained increasing 1 attention 28 . Modeling and simulation methods play significant roles in predicting the degradation behavior of stainless steel containers in corrosive environment. For example, the atmospheric corrosion of stainless steel in stores (ACSIS) model developed by NDA/RWM 29 , which aims at evaluating the ignition possibility and extent of localized corrosion and SCC, could greatly help assess the extent of corrosion damage to intermediate level waste (ILW) containers for the required in-service period. The ACSIS model has been successfully implemented through a software tool and applied in atmospheric corrosion in ILW stores. We also expect a similar numerical model could assist us better understanding the corrosion behavior of the SNF canisters efficiently and economically.
In this paper, a critical review about the application of numerical methods in the study of CISCC for the in-service canister was conducted. The factors of mechanical stress field, electrochemical reactions and the material microstructure, hydrogen-assisted cracking were considered. Then the coupling effects between different fields in several recent studies were discussed. In the last part of this paper, a multi-physical deterministic model was proposed for the future study of CISCC behavior of the SNF canisters.

INFLUENCE OF STRESS, ENVIRONMENT, AND MICROSTRUCTURE Mechanical stress and strain field
Tensile stress is essential for the CISCC occurrence. As discussed in ref. 30 , the degradation caused by pitting corrosion is slow, while residual tensile stresses will lead the degradation rate largely increased. Residual tensile stresses are those tensile stresses that remain in the material after the original causes of the stresses have been removed. If the remained stresses are above the yield point, then plastic deformation would be induced, which has great impact on the SCC behavior of the canister. For the in-service canister, the stresses can be induced by the manufacturing stages (e.g., machining and rolling), fabrication stages (e.g., welding), and in-service conditions (e.g., temperature and pressure). It was revealed that welding process during the canister construction would contribute most of the residual tensile stress 12 .
To understand the magnitude and distribution of weld residual stress, nuclear regulatory commission (NRC) used finite element method (FEM) to simulate the weld residual stresses in a representative dry storage austenitic stainless steel canister 31 . This research adopted a 2-dimensional (2D) sequentially coupled thermal-structural FE model and considered the influence of isotropic and kinematic hardening laws. The results show that the applied constraints had little effect on the stress distribution when compared to the large differences exhibited by the different hardening rules. In this investigation, the longitudinal weld would result in high through-thickness axial direction tensile stress in the weld and heat affected zone (HAZ) of the canister. While the circumferential weld (which is not shown here) would result in high through-thickness hoop direction tensile stress. Both the tensile stresses present a value near or above the yield point of the canister material. They predicted that the through-thickness tensile residual stress exists in at least one component direction in the circumferential, longitudinal, base plate, and lid closure welds of the canister, which would offer the information that there are sufficient stresses to allow for SCC initiation and potential through-wall growth if exposed to a corrosive environment. However, a 2D model was used in their study, by which the boundary conditions have to be modified from the threedimensional (3D) condition to adapt to the 2D model. Besides, the 2D model cannot consider the interaction between different welds, as well as the overall welding sequence for different parts, all of which have a significant influence on the measured results. Meanwhile, some assumptions such as the ignorance of the fabrication stresses and the simplification of the weld geometry, interpass cooling time, are made, which means that corresponding measurement data on the canister is necessary to make their conclusions more reliable. Regarding to these issues, our recent studies proposed to use a 3D sequentially coupled thermal-structural FE model to simulate the weld residual stresses in an in-service canister, and comparison to the actual measurement data was conducted for validation 13 .
The geometry model and the related welding parameters in our studies were used according to the DOE report 12 . For the research in ref. 12 . a full-diameter mockup canister was first constructed according to the major manufacturing procedures used for building the fielded SNF interim storage canisters, followed by the measurement of residual stresses in the mockup canister using different experiment methods. The geometry feature of the weld beads for the longitudinal and circumferential welds was experimentally characterized and applied in the FE model. The Goldak double-ellipsoidal heat source model was used for the simulation of the arc heat distribution, and the parameters were adjusted to make the total fusion zone in the numerical model be close to the experimental measurement. Detailed boundary conditions and the simulation parameters could be obtained from ref. 13 . Though there are still some assumptions, such as the ignorance of the fabrication stress, which was demonstrated to be with small magnitude by ref. 12 and simplification of the weld bead morphology due to small transition angle, similar stress distributions could be predicted for longitudinal weld and circumferential weld when compared to the measurement data in ref. 12 . Fig. 1(b), 1(c) show the comparison between the finite element analysis (FEA) results and measurement data for longitudinal weld. A balancing compressive stress, which is expected for hoop residual stress, is existed in the radial direction. Both of the figures show an overall good match between the simulation data and measurement data. The much lower stress at the sites close to the surface from the deep-hole drilling (DHD) test is induced by the grinding process after the finishing of welding. It demonstrates a through-wall high tensile stress for the longitudinal weld, which is the same as that indicated in ref. 31 . The maximum value is close to 420 MPa along the axial direction, much larger than the yield strength (~250 MPa) of the material, indicating a large plastic deformation has happened for the weld material during the solidification. The tensile stress region extends~25-30 mm from the weld centerline. For the circumferential weld (not shown here), through-wall high tensile hoop stress was found, and the maximum value is close to 310 MPa. The tensile stress region, which is a little wider than the longitudinal weld, extends around 40-50 mm from the weld centerline.
Under marine environment, there is potential pitting corrosion happening on the surface of the canister. The pits formed on the surface of the structure tend to intensify the local stress field and hence induce the transition from pit to crack if the critical condition is satisfied. It is realized that crack readily nucleates at the pits having the maximum stress concentration factor (SCF) value. Moreover, the severity of pits varies with the depth, shape and size. We could not find any studies related to the influence of pitting on the SCF for the storage canister structure. While this type of investigation has been done on other structures.
In 2009, Cerit investigated the stress concentration at corrosion pit by using a series of 3D semi-elliptical pitted models 32 . He claimed that the magnitude of SCF is mainly affected by the value of pit aspect ratio (a/2c). The occurrence of maximum stress could X. Wu et al. happen either at the pitting bottom or just below the pitting mouth. The nucleation of a secondary pit around the bottom of the primary one would lead to a large increase of the SCF value. Meanwhile, the interacting effect between two neighboring pits was also observed 33 , and it was found that the increasing angle of two transverse pits would result in nonlinear decrease of the SCF value. The interacting effect can be ignored under an ultra-small distance condition. A model with random arrangement of multipits is close to the actual pitting conditions of the canister, which was conceptually introduced in ref. 34 .
For the canister, the same concept was used in our studies to build a series of FE models to calculate the stress concentration phenomenon induced by pitting in the canister. The geometry and material database were taken from the actual canister data. Moreover, a uniaxial tensile load equal to 250 MPa (maximum tensile stress in the HAZ) was applied at the two ends of the canister. The pits with different shapes and depths, open widths were considered. A similar stress concentration results as that found in refs. 32-34 could be expected. Through this model, the distribution of principle stress and principle strain for different types of pits were also derived. Figure 2 gives the analysis results for a pit with semi-ellipsoidal shape as an example. It is seen that the existence of a tensile load would induce a maximum stress distribution along the pit base. The maximum value of the stress is 365.2 MPa, which is above the yield point of the 304 SS, and a SCF value as 1.46 was predicted. For the plastic strain distribution ( Fig. 2(b)), maximum value locates in the pit shoulder region, and is close to 3.72e-3. This type of stress and strain distribution of a specific pit was also revealed in ref. 35,36 and similar distribution trend was predicted. It was explained that the plastic deformation causes the stress to redistribute towards the pit base. From linear elastic assumption, it is known that the pit-to-crack transition would happen in the stress concentration region. While if the plasticity is considered, an energy-based criterion should be applied. The localized large plastic deformation could offer the energetically most favorable site for the nucleation of crack. In the experiment, some cracks were found to initiate at the pit shoulder rather than pit base, as shown in Fig. 2(c). This research could provide a mechanic-based explanation. The initiation of stress corrosion cracks at the shoulders of the surface pits/defects were widely reported [37][38][39][40][41] , which could be complementary to the classical Kondo theory 42 . Undoubtedly, the mechanic-based explanation is meaningful to its further applications.

Electrochemical influences
Generally speaking, for corrosion modeling the coupled ionic species transport equations have to be solved first under designated electrolyte solution and potential distribution. Then the corresponding current density distribution is determined from the potential, based on which the anodic dissolution rate was calculated via Faraday's law 43 . Transport of species i in the dilute solution is usually represented by Nernst-Planck equation in the numerical models, for which the basic form is formulated as 44 : where c i is the concentration, D is the diffusion coefficient, z is the charge, F is the Faraday constant, Φ is the electric potential, ν is the relative velocity, R is the ideal gas constant, T is the temperature, and N is rate of homogeneous production of species i. In most cases, the electrolyte solution is assumed as no concentration gradient and electro-neutral, which could result in the use of the Laplace equation to model the current and potential distributions 45 : For the canister material in the marine corrosive environment, the electrochemical anodic and cathodic reactions of the stainless steel are mainly caused by the iron oxidation (the oxidation of other elements such as Cr, Ni 46 are ignored to simplify the simulation process) and hydrogen evolution or oxygen reduction 47,48 : At equilibrium, rate of the forward reaction for Eqs. (3)(4)(5) is equal to the reverse case. The potential and current density at the equilibrium state are called equilibrium potential and exchange current density. When the equilibrium state is disturbed, there is a net current accompanied by new potential at the electrodes, which is termed polarization. Corrosion of the steel is activationcontrolled. The activation polarization of the electrodes is described by: Cathode where η is the over potential, β is the Tafel slope (V/dec), i is the current density (A/mm 2 ) and i 0 is the exchange current density (A/mm 2 ), a and c represent anodic and cathodic reactions, respectively.
Then the polarized potential for the cathodic and anodic reaction can be written as: where ϕ 0 is the equilibrium potential (V). Based on the above governing equations, several numerical models were proposed to calculate the corrosion in steel structure. In 2010, Luc et al. 49 developed a simple and significant improved inverse relation for the calculation of current density from measured potential for cathodic reaction. The proposed inverse relation was determined as: where i L is the limiting current density of the cathodic reaction (A/ mm 2 ). i L can be calculated as follows: where D O2 is the effective oxygen diffusion coefficient (m 2 /s), δ( = 0.005 mm) is the thickness of the stagnant layer of electrolyte. C O2 is the concentration of oxygen around the steel (mol/l). By using this inverse relation, the current density was demonstrated to be determined accurately from the potential. Then this inverse relation was applied in a unified and straightforward algorithm to perform different types of steel corrosion modeling 50 . The proposed algorithm can be computationally more efficient and capable of modeling complex geometries and incorporating parameters that vary over the domain. From the comparison of the proposed model with those based a macroand-micro-cell modeling, it is obvious that they predicted very Then, in 2016, Lu et al. 52 used a numerical model to predict the corrosion rate of weld joint, which is interesting to the canister corrosion case study because in most possibility, the SCC would happen in the vicinity of canister weld region. While in their model, the mechanical stress was not taken into consideration. The corrosion rates obtained from the FE model were calculated from current density based on Faraday's law, and were compared with those estimated from mixed potential theory and experimental techniques. As shown in Fig. 3(a), the numerical model can predict a 4.8 μm corrosion depth at the base metal side of weld joint, and for the immersion test, the corrosion depth is 5.9 μm, which means the numerical prediction is within 20 pct. of that obtained from the immersion experiment. For Fig. 3(b), the numerical simulation can predict a corrosion depth as 1.72 mm at the anode zone of the weld joint, and the constant potential polarization test could obtain a value as 1.62 mm. Therefore, the maximum corrosion depth predicted by FE model is within 10 pct. with respect to the experimental test.
For the numerical simulation of the canister corrosion, a similar methodology can be adopted, while the influence of the residual stress must be considered. Besides the general FE models based on electrochemical reaction, the corrosion problem was also investigated by extended finite element methods 53 , finite volume method 54 , phase-field method 55 , atomic simulation method 56 , and cellular automata 57 , etc., which could offer as instructive references for the research of SNF canister corrosion.

Coupling between different fields
For the simulation process of the CISCC in canister, besides the electrochemical factors, stress is an important factor that can influence the corrosion behavior significantly. Thus, for the actual CISCC modeling, the synergistic effects of local stress concentration, corrosion electrochemical reactions, and the pitting geometry should be considered. It was described previously that the pitting geometry could influence the local stress, and the governing functions for the corrosion reactions were clearly indicated. Besides, the existence of localized stress could interact with the corrosion reactions. According to the mechanochemical interaction theory developed by Gutman 58,59 , the elastic and plastic deformations could affect the electrochemical activity, resulting in changes of electrode potential. The change in electrode potential induced by deformation is described by: Plastic deformation : where ΔP is the applied load (Pa), V m is the molar volume of the material (m 3 /mol), z is the valence of the metal ions, F is Faraday's constant (96,485 C/mol), T is the temperature, R is ideal gas constant (8.314 J/mol K), ν is the orientation-dependent factor, α is a constant of 10 9 cm -2 to 10 11 cm -2 , ε p is the plastic strain, N 0 is the initial density of dislocations prior to plastic deformation. For the canister case, both elastic and plastic deformation could contribute the change in electrode potential by: where ΔP m is the external pressure within elastic deformation limit, which is equal to 1/3 of the yield strength of the material. The anodic dissolution current density, i, of the metal can be expressed by: where i a and i c are the anodic and cathodic current densities, respectively. Based on these functions, there are some numerical studies conducted to investigate the corrosion of the material under coupling effect of mechanical stress field and electrochemical field. In 2013, Wang et al. 60 coupled a cellular automata model with FE method to describe the growth and the stability transition of pitting corrosion under mechanical stress effect in stainless steel. Under the stress load condition, the current was found to be far higher than that without stress loading. As the corrosion growth is determined by the anodic dissolution current by Eq. 14, the pitting growth rate for stress-accompanied condition is then determined to be much larger than that of stress-free condition. They also demonstrated that exist of stress would cause the pitting corrosion to reach a stable growth more easily. By using the same methodology, they further investigated the interactions between two pits under mechanoelectrochemical (M-E) effects 61 . The interactions within double pits would lead to a higher plastic strain, generating a quicker increase of transient current when compared to a single pit. The pit interactions could also result in easier steps for the metastable pits entering the stable growth regime.
In 2013, Xu et al. 62 developed a multi-physics FE model to study the corrosion of pipeline. In their model, the mechanical effect was coupled into the electrochemical reactions by using the above governing functions. The corrosion potential and corrosion current density were well compared with the experimental measurements. They figured out that the electrochemical corrosion behavior was little influenced by the mechanical stress under an elastic deformation, while it was remarkably increased under plastic deformation condition. Figure 4 shows the effect of M-E on anodic dissolution current density, which can be used to determine the corrosion rate according to Eq. 15. As shown in Fig. 4(a), in the absence of tensile strain, the metal dissolutes uniformly inside the defect with a dissolution current density as 3.49 μA/cm 2 . When a tensile stain is applied, the current density at the center is increased, while the current densities at the sides are still low. The increase will be enhanced by larger tensile strain. Also, by increasing the depth of the defect, the corrosion rate at the defect center would be further accelerated ( Fig. 4(b)). In Xu's thesis 63 , based on the same M-E coupling concept, a novel FE model was further proposed to enable the assessment and prediction of the time-dependent growth of corrosion defect on pipelines, which could provide a promising alternative for assessing the life of the corroded structure.
In 2016, Wang et al. 64 also developed a similar numerical methodology to evaluate the M-E performance of the corroded steel structure under the conditions of external loading and weld residual stresses. They predicted a dynamical response of the corrosion behavior to the loading conditions as well, and the M-E corrosion is more affected by plastic deformation. Though these numerical models are idealized with the exclusion of coatings and passivation and the simplified corrosion processes under specific loading and boundary conditions, they are able to offer important corrosion information to the structures once these models are calibrated by experiment. Meanwhile, it is very important to get instructive inspiration from these models to study the CISCC behavior of the SNF canisters. However, it should be noted that the assumption of low-concentration gradient in these studies may be not applicable once a pit/defect generates because it would result in locally concentrations that can be substantially different from bulk conditions. The future studies should pay special attention to this issue.

Influence of susceptible microstructure
Corrosion behavior in metals was experimentally demonstrated to be strongly depended upon microstructure features at the corroding surface [65][66][67][68][69][70][71][72][73][74] . For the canister material in different regions, the microstructures (inclusions, phases, grain boundaries, grain sizes, crystallographic orientations, flaws, dislocation density) would change the CISCC possibility significantly. To expand the knowledge base of the microstructural effects on pitting corrosion, several models under computational framework have been recently developed.
From 2012, Qidwai's group began to focus on the numerical simulation of pitting corrosion with the consideration of microstructure [75][76][77][78][79] . In the work in 2017, they systematically used a numerical model to investigate the growing process of stable pits under the consideration of the influence of crystallographic plane orientation. They took use of the orientation data measured from experiment, and then used a relation based on the experimental evidence to describe how the corrosion potential and polarization behavior depend on the crystallographic plane orientation. As shown in Fig. 5(a, b), serial sectioning and electron backscatter diffraction (EBSD) mapping methods were first adopted to uncover the microstructure information, then ten random locations (one is shown in Fig. 5(a) as an example) were selected from the microstructure obtained to set the positions of original pits. The boundary conditions used in the model are also depicted ( Fig. 5(b)). After determining the local potential distribution with respect to the crystallographic orientations, the final pit shapes after some time corrosion are predicted in Fig. 5(c). It is seen that the crystallographic orientation significantly affects local pit growth rates. The higher potential variation (10%) generally causes more irregular shapes to form after 100 seconds growth time. Additionally, the local changes in pit shapes would change the stress concentrations under mechanical loading, which could further affect the pit-to-crack transition behavior. Though the computational technique developed in their work needs to be further improved (Some important microstructural features were ignored), it can pave the way for studying complicated corrosion process involving electrochemical reactions and microstructural characteristics, which is of great significance for the research of CISCC in SNF canisters. Also, considering the above described M-E coupling concept, the future model should involve this microstructure effect into the M-E coupling.
Besides the above crystallographic orientation, the incorporation of some other realistic microstructures of the material into a numerical simulation model for pitting and cracking corrosion studies have also been recently explored.
In 2007, Jivkov et al. 80 applied a 3D finite element model to analyze the effects of resistant grain boundary on the intergranular stress corrosion cracking (IGSCC) growth rate. They demonstrated that the propagation rates of short crack could be significantly retarded by boundary crack bridging ligaments.
Then in 2010, Arafin and Szpunar 81 coupled a Markov Chain-Monte Carlo model with the features of grain boundary and info of orientation-dependent vulnerability statistics to model the IGSCC propagation in pipeline steel. Their model could predict data in good agreement with experiment observations. Chen and Bobaru 82 investigated the influence of microstructural heterogeneities on pitting corrosion in metals based on   5 Microstructure-sensitive modeling of pitting corrosion. a shows a random location in the microstructure (from EBSD mapping) for pit growth and b shows its associated boundary value problem within the electrolyte domain. c Comparison of grown pits at 100s at four locations in the microstructure for 5 pct. and 10 pct. maximum variations in potential with respect to the crystallographic orientation. c is metal ion concentration 79 . Reproduced with permission from ref. 79 , copyright (Elsevier, 2017). peridynamic models. The pit shapes were found to strongly depend on the particular microstructures.
Meanwhile, Mei et al. 83 employed a phase-field model that incorporated microstructural information to uncover the effect of crystallographic orientation in the activation-controlled pitting corrosion phenomena.
Yu and Pan 84 developed a time-dependent FE model that takes the complicated microstructure phases in iron into consideration to investigate the mechanisms of pitting corrosion in steels. These models could significantly inspire the construction of the potential CISCC numerical model for the SNF canister.
Hydrogen-assisted stress corrosion cracking For the SCC mechanisms of the metal under corrosion environment, there are usually two groups can be divided into: anodic dissolution (AD) and hydrogen-assisted cracking (HAC) [85][86][87][88] . Hydrogen uptake in the canister material can initiate or speedup the degradation. There are a lot of demonstrations showing the hydrogen-induced degradation, such as external crack, internal crack and blisters [89][90][91] . During the anodic dissolution of the canister material in CISCC, the cathodic reactions could produce atomic hydrogen, which can combine to form the dihydrogen or get into the iron lattice, leading to the increase of hydrogen pressure, because of which microcracks can be easily nucleated and eventually catastrophic failures would happen. To simulate the process of HAC, two issues including diffusion of hydrogen and hydrogen-induced cracking are important. For the simulation of hydrogen diffusion, it was demonstrated that the mechanical stress field would significantly influence the hydrogen distribution [92][93][94] . Therefore, it is essential to consider the hydrogen diffusion process in the canister material under the residual stress field. The stress-assisted diffusion of hydrogen was expressed by a modified Fick's second law 93,95,96 : where c is the hydrogen concentration, D is the diffusion coefficient, V H is the partial molar volume of hydrogen in steel. R is the mole gas constant (8.314 J/mol K), T is the absolute temperature (K), and σ h is the hydrostatic stress.
With the existence of hydrogen, the intrinsic fracture resistance would be reduced. Based on the fracture mechanics, the effect of hydrogen localization on the threshold of K IC for cracking propagation prediction can be expressed as 97 : where α, α", and β' are constants, c H is the diffusible hydrogen concentration, σ ys is the yield stress, and K IG is the intrinsic Griffith toughness and can be expressed as: In fracture mechanics, the criterion for the crack to propagate can be described as: (19) Therefore, the existence of hydrogen could reduce the threshold value of the stress intensity factor, which makes it easier to meet the criterion described in Eq. 19.
Based on the above equations, Lee et al. 88 derived a function to describe the SCC growth rate for the HAC part: where a is the initial crack depth, ξ and η are the displacement and the time variables, respectively. The SCC growth rate for AD part can be described by Faraday's law as da dt where i a is the anodic dissolution current density, M is the molar mass, n is the valence and ρ is the density. Then the total corrosion rate is: There are some numerical models developed to investigate the HAC process. In 2011, Raykar et al. 98 presented a strategy that combined analytical and FE solution to simulate the HAC growth. They built a relation to describe the influence of plastic deformation rate on the concentration dependent decrease in cohesive strength, and their results compared well with the experimental data. Traidia et al. 89 applied a comprehensive FE model to investigate the HAC in steel pipelines under the corrosive media of sulfurous compounds. In this model, they considered the coupling effects between the diffusion of non-Fickian hydrogen, building-up of pressure and extension of cracks. The values of HAC growth rate under the effects of pipe wall thickness, initial crack radius, position, and fracture toughness were studied. The most referable model for our research is the work by Lee et al. 88 in 2015. They proposed a comprehensive numerical model to analyze the growth rate of SCC for aluminum alloy 7050-T6 by a combined AD and HAC mechanism. In their model, they coupled the crack chemistry and different types of hydrogen diffusion. Their results demonstrate that AD provides the critical conditions for HAC to happen, which could lead to the stepwise propagation of the crack as observed in experiment. The total crack growth rate due to the combined effects of AD and HAC has also been determined. The results indicated that for pH > 3, the calculated (da/dt) AD well presents the crack growth in this part, while for pH <3, the calculated (da/dt) HAC presents a better fit. It is because at lower pH, the hydrogen concentration at the crack front will be higher, which makes the trigger of the HAC part easier to happen. After the HAC part is triggered, the corrosion rate of this part is much faster than the AD part. The crack growth rate shows a linear relationship with the cathodic current density on a log-log plot. Thus through their model, it was proved that both the AD and HAC parts contribute to the crack growth, but they dominate individually at different pH levels. However, their model ignored the influence of mechanical stress field on the electrochemical corrosion reactions, which is very important for the CISCC modeling process of the SNF canister. As a result, our proposed numerical model will involve the factor of mechanical stress based on their concept.
In some cases, the concentration of hydrogen may be too small to cause HAC due to the requirement of a critical hydrogen content. Nonetheless, hydrogen accumulating in the sites of cracktip could interact with the stresses, generating great charges for free-energy at the vicinity of crack-tip. Based on thermodynamics analysis, Mao et al. [99][100][101][102][103][104] conducted a series of studies and proposed that the increase in free-energy of the metals enhances the anodic dissolution cracking rate. They concluded that the SCC of steels in low and near-neutral pH solution was controlled by hydrogen-facilitated anodic dissolution, in which the synergistic effect of hydrogen and stress on anodic dissolution rate was determined by: where i A is the anodic dissolution current density without hydrogen and stress, and where ΔU is the internal energy change, ΔS is the entropy change, M is the molar weight, σ 1 , σ 2 , σ 3 are the principal stresses, E is the Young's modulus, ρ is the density, σ H is a volume stress given by (σ 1 + σ 2 + σ 3 )/3, V H is the partial molar volume of hydrogen in the metal. κ H , κ σ , and κ Hσ represent the effect factors induced by the exist of hydrogen, stress and their synergistic function, respectively. κ H is the effect factor for the influence of hydrogen on anodic dissolution under stress-free condition; κ σ is the effect factor for the stress under hydrogenfree condition; κ Hσ represents the synergistic effect of the presence of both hydrogen and stress. Based on the above governing functions, the interaction between the mechanical stress and hydrogen accumulated around the crack-tip will bring large changes to free-energy, which leads to an accelerated SCC growth rate for the anodic dissolution.
In 2007, Cheng et al. 105 confirmed that the effect factor of hydrogen is only 1.53 for enhanced anodic dissolution of steel when the anodic potential is −0.4 V (SCE) at a solution with near-neutral condition (pH is close to 7), which means hydrogen-enhanced anodic dissolution by itself may contribute little to the high crack growth rate in pipe steel for nearneutral solution condition. To further investigate the mechanism of near-neutral pH SCC of pipelines, Cheng updated Mao's model by adding an expression of hydrogen concentration factor. The synergistic effect of hydrogen and stress on anodic dissolution was given as 106 : where x, y represent the amount of hydrogen atoms permeating into the unstressed and stressed steel, respectively. G H2 is the formation free-energy of the hydrogen molecule, G H abs is the formation free-energy of the absorbed hydrogen atoms in the steel. κ ΔCH represents the effect of the different hydrogen concentration for the stressed and unstressed material on the reaction of anodic dissolution.
Based on Eq. (27), Cheng pointed out that to quantitatively predict the SCC growth rate in pipeline steel, the key relies on how to determine the relationship between anodic dissolution rate at crack-tip and the concentration of hydrogen.
However, in 2009, Lu et al. 107 proposed that the SCC of pipeline steel in the anaerobic groundwater of near-neutral pH is unlikely induced by the hydrogen-facilitated anodic dissolution mechanism. Their thermodynamics analysis indicated that the synergistic effect of mechanical stress and dissolved hydrogen on active dissolution can be neglected. Meanwhile, the effect of dissolved hydrogen, mechanical deformation (either elastic or plastic) on anodic dissolution of pipeline steel in groundwater with near-neutral pH condition is quite limited. Then in another study 108 , they revealed that the hydrogen-enhanced SCC should be caused by the hydrogen-induced passivity degradation. Nonetheless, they also pointed out the anodic dissolution is nearly insusceptible to the tensile stress the material experienced, which may be doubtable because there are some researches demonstrating that the M-E effect developed on the steel is significant under plastic strain, and this phenomenon was well explained by the Gutman's mechanochemistry theory 58,59,[62][63][64][109][110][111] .
Undoubtedly, there is still uncertainty for the role of hydrogen in SCC of materials, and the suitability of using thermodynamics for analyzing this question needs more validations. More experiments and theories are expected, and our numerical model may contribute to it. We plan to focus our model on the influence of hydrogen on the corrosion potential, and add hydrogen role into the Gutman's mechanochemistry theory. Together with the experiment validation, we hope this model could uncover the effect of hydrogen, mechanical stress, and their synergistic effect on the SCC of stainless steel canister.
The newly proposed multi-physical model for SNF canisters Based on the research progress described above, we are trying to develop a multi-physics numerical model, which could further involve more influence factors and corrosion mechanisms to simulate and predict the CISCC in the Type 304 SS SNF canisters. In our proposed model, we plan to consider the influence of weld residual stress, corrosion environment (temperature, pH, O 2 ), electrochemical reactions, microstructures, and hydrogen embrittlement, so as to predict the corrosion behavior more accurately. In this section, we describe the theory of the newly proposed model, and the primary results based on the assumptions we already have right now. Figure 6 describes the simulation model for the CISCC prediction. a shallow defect is assumed to exist on the surface of the canister (Fig. 6(a)), and the canister is in the corrosive environment. A 2D model was used, and the thickness and length parameters were taken from the realistic canister structure ( Fig. 6(b)). The solution interacts with the outer surface of the canister with a thickness equal to half of the outer diameter. An elliptical defect exists in the center site of the canister with dimensions noted in Fig. 6(c). The future model will update to further study the influence of the shape and geometry parameters of the defects. The corrosion reactions happen in the interface between metal and electrolyte. The anodic and cathodic reactions were described as the oxidation of iron and reduction of hydrogen, respectively, as described in Eqs. (3) and (5). A uniaxial tensile load was applied at the edge of the canister to simulate the influence of weld residual stress of the canister, and the other edge of the metal part was fixed. The electrically isolated boundary condition was used for the solution borders, except that the solution/metal interface was set as a free boundary. The canister bottom was set as electric grounding. The model was meshed with triangular elements, with a much finer mesh in the vicinity of the defect.
By using this FE model, the AD corrosion behavior of the canister under the effect of electrochemical reactions and mechanical stress field can be considered, which is similar to that reported in ref. 62 . More importantly, the diffusion of hydrogen in the metal, which could induce the HAC, was also considered, and the total final cracking growth rate is determined by the combination of AD part and HAZ part. The hydrogen used in the HAC part is generated due to the cathodic reaction during the AD part, by which these two parts can be combined together. The role of dissolved hydrogen in the anodic dissolution part, which may be important for the corrosion process as noted above, has not been considered in the current model.
For the solution procedures of this multi-physics model, the first step is to determine the stress distribution under the tensile load, so as to obtain the stress/strain concentration induced by the pitting corrosion under weld residual stress. The elastoplastic solid stress simulation was conducted by using the Type 304 SS material database. The isotropic hardening model embedded in COMSOL was used for the stress calculation, the hardening function σ h is formulated as (http://cn.comsol.com/model/ elastoplastic-analysis-of-holed-plate-244): where σ e is the effective (von Mises) stress, E is Young's modulus, σ exp is stress function derived from measured stress-strain curve, ε pe is the effective plastic strain, σ yso is the initial yield stress. All the parameters used here are the same as that used in the weld residual stress simulation model. The electrochemical reactions are described in Eqs. (3) and Eq. (5), and the electrode kinetics is represented by Eqs. (6)(7)(8)(9). The equilibrium potentials can be derived by the standard equilibrium potentials and ion concentration based on Nernst equations 112 : where ϕ a0 (SHE) and ϕ c0 (SHE) represent the equilibrium potentials under standard hydrogen electrode. The solution of Eqs. (29) and (30) need the values of the concentration of ferrous ions and pH, which were taken as 10 −6 M and 8.07, respectively. The other electrochemical parameters, like Tafel slope and exchange current density, were obtained through the potentiodynamics polarization curve as shown in Fig. 7 (only part of the polarization curve is shown here). The corrosion test was conducted using the Gamry Instruments@Reference 600. The test material is Type 304L SS from mockup canister, and the electrolyte is ASTM seawater. Based on the curves obtained from the corrosion test, we can derive the Tafel slope information for anodic and cathodic reactions. The exchange current densities are the current densities in the curves corresponding to the equilibrium potentials. The final derived values are summarized in Table 1.
The influence of mechanical stress and strain on the electrochemical reactions is described by Eqs. (12)(13)(14)(15). Then the corrosion rate induced by the AD part can be described by Faraday's law as shown in Eq. (21).
For the HAC part, several assumptions were made in our present model: (1) the diffusion of hydrogen into the electrolyte was ignored; (2) the diffusion of hydrogen from electrolyte into metal was ignored; (3) the original hydrogen in the metal was ignored. The hydrogen is generated by the cathodic reaction in AD part as represented by Eq. (5). The flux of hydrogen N H is derived by the hydrogen reduction current at the electrode surface as: where υ is the stoichiometric coefficient, i c is the local current density, n is the number of participating electrons, F is the Faraday constant.
The diffusion of hydrogen is considered under the influence of mechanical stress field by Eq. (16). Then based on the fracture mechanics equations (Eqs. (17)(18)(19)), the crack growth rate by the HAC part is represented by Eq. (20). The final total crack growth rate is the combination of the AD and HAC part (Eq. (22)). The overall main parameters used in the calculation of the functions are provided in Table 2.
Based on the assumptions and parameters we used, some initial results can be derived from the current model.  density in the electrolyte and the corresponding plastic strain in the metal. It can be seen that the pitting defect grows under the corrosive environment. The tensile load would induce plastic strain into the metal, and the plastic strain is concentrated in the defect center. The generated potential and current density concentrate in the high plastic strain region. Figure 8(b) gives the moving boundary results of the defect under different corrosion time. It shows that after 5 years, the corrosion is only induced by the AD part, while the HAC part is not triggered. After 30 years, the HAC part has been triggered and the grow rate induced by HAC part is much larger than the AD part. This is because the role of HE part can only be effective when the concentration of hydrogen reaches a critical value. We are working to determine this critical hydrogen concentration.

CONCLUSIONS
This paper summarized the recent studies about the numerical modeling and simulation work to build a model relevant to CISCC for the stainless steel dry storage canister. The published literatures demonstrate the shortage of data on this area, especially for the coupling between different factors, like the mechanical stress, electrochemical reactions, hydrogen, microstructures, temperature, pH, ion concentrations. Based on this situation, we proposed a multi-physics model trying to involve these factors into one single model to study the CISCC behavior of the SNF canisters. In our present model, the mechanical stress field has already been successfully coupled with the electrochemical reactions, and the AD part and HAC part was combined to consider the SCC growth rate. Although there are some simplifications and assumptions in the proposed model, e.g., the simplified anodic reactions, the simplified diffusion of cathodic hydrogen, the ignorance of passivity & pit stability, the neglect of oxygen reduction. The proposed model can reflect the actual conditions of the in-service canister and provide useful information. Our group is still working to involve more aspects to make it more mature, like the consideration of the microstructures in the defect front (phases, grain orientations, grain size and so on), the consideration of passivity & pit stability, the cathodic oxygen reduction, and the consideration of the surface tension term during the interaction between the metal and solution. Besides, validation of the modeling data is urgently needed from different aspects to demonstrate the reliability and capacity of the proposed model. Currently, we have designed four-point-bend specimens that consider the factors from microstructure, stress, and corrosion environment and placed the specimens into the humidity chamber to investigate the pitting and pit-to-crack transition behavior, though months, even years are needed before some comparable results could be obtained.