Hydride corrosion kinetics on metallic surface: a multiphase-field modeling

The hydride growth on metallic surface can cause material failure, which is a significant form of pitting corrosion. A new multiphase-field model is constructed to study the hydride corrosion kinetics, in which three order parameters are introduced to represent the passive film, hydride and metal phase, respectively. Coupling with hydrogen concentration field and elastic strain field, this model not only presents the growth of hydride and the rupture of passive film, but also reveals the hydrogen diffusion mechanism and the effect of strain energy on pitting corrosion process. The simulation shows the semi-ellipsoidal cerium hydride forms and grows near the passive film/cerium interface. During this process, the passive film on the upper side of the hydride is also hydrogenated or peeled off, resulting in faster hydrogen transport, which in turn promotes the growth of hydride. The formation of cerium hydride causes volume expansion, and the strain energy is mainly distributed around the hydride, which inhibits its growth. The present study contributes to understanding the formation mechanism of hydride corrosion at mesoscale, especially the pitting corrosion kinetics of rare earth metals.


Introduction
The hydride growth on the metallic surface is a significant form of pitting corrosion, which can cause the destruction of materials. Hydride corrosion often occurs on the surface of rare earth metals, such as cerium [1][2][3], gadolinium [4][5][6], holmium [7,8], uranium [9][10][11], and Plutonium [12][13][14] etc, which can form stable metal-hydride phases [15][16][17][18] when they are exposed to low level hydrogen environment for a long time. A number of transition metals, for example, titanium [19,20], zirconium [21][22][23], are also known to form similar surface-breaking hydrides by the reaction of hydrogen and metals. As nuclear engineering materials, Plutonium and uranium are often troubled by hydride corrosion, so the research on the hydride growth mechanism is particularly important. Because pure cerium has similar chemical and physical properties to nuclear materials such as Plutonium and uranium, it is often used as a model metal to study the hydride corrosion kinetics of Plutonium and uranium. Understanding of the mesoscopic hydride corrosion kinetics can help to analyze and predict the hydriding process, which is beneficial for avoiding or reducing industrial hazards caused by the hydride corrosion.
The hydride corrosion kinetics on metallic surface generally involves the interaction between passive film, metal and hydride [9,11,12]. Composed of metal oxide, the passive film covers the metallic interface, which not only affects adsorption/chemisorption of the hydrogen but generally acts as a diffusion barrier to hydrogen [24]. Therefore, the damage of the passive film often makes the metallic surface partially exposed to the environment, resulting in a faster hydrogen adsorption and diffusion rate [13]. The nucleation and the precipitation of hydride usually occur near the passive layer/metal interface [25], and the initial stage of hydride formation depends on the adsorption, dissociation and transport of the hydrogen [26,27]. The hydriding reaction of metal M is suggested to occur via M+x[H]=MH x , where [H] is the hydrogen dissolved in metal or passive film, and x is the H/M ratio [18,28]. In addition, due to the hydride density is often less than that of the metal phase, the growth of the hydride can cause local expansion and generate internal stress around the hydride. As ionic crystals, the passive layer and hydride are brittle phases, which are prone to crack under the internal stress. Strained by hydride, the metal phase may undergo plastic deformation and bulge with the hydride [29].
Based on previous experimental observations, the hydriding process involves multiphase interactions of passive film, metal and hydride, the diffusion of hydrogen, and the effect of internal stress. However, previous experiments were difficult to reveal the dynamic effects of these various factors. For the study of corrosion kinetics, the mesoscopic simulation [30] was introduced to understand the mechanism of corrosion process and provide guidance for corrosion prevention. In recent years, Mai et al [31][32][33] designed a pitting corrosion model based on phase field and concentration field to simulate the activation and diffusion controlled pitting corrosion phenomena in metallic materials, which can reproduce different portions of the electrochemical corrosion process associated with the activation-controlled, diffusion-controlled, and mixed-controlled corrosion kinetics. Ansari et al [34,35] introduced a phase-field model to study the pitting corrosion kinetics in metallic materials, which couples the phase field, electric field and solute concentration field. Through this model, both the ion transport in the electrolyte and the electrochemical reactions at the electrolyte/metal interface are explicitly taken into consideration. Tsuyuki et al [36] proposed a new phase-field model to simulate the pHdependent corrosion of iron, which is based on Bockris's iron dissolution mechanism to describe the pH dependence of the corrosion rate. All these corrosion models are mainly suitable for electrochemical corrosion of metal-electrolyte systems.
However, the pitting corrosion caused by hydride growth is quite special, especially the kinetic characteristics of the passive film and the strain energy effect caused by hydride precipitation should be taken into consideration, which were not achieved by the previous corrosion models. The multiphase-field model [37] coupled with hydrogen concentration field and strain energy field [38] can be used to solve the problem of hydride pitting corrosion modeling, which can be implemented to simulate the morphology evolution of hydride, the distribution evolution of hydrogen and the strain energy density around hydride. Therefore, the present model can not only reveal the hydride pitting corrosion mechanism at mesoscale, but also predict the hydride growth on metallic surface. And the understanding of the hydride corrosion kinetics on metallic surface is beneficial to the corrosion prevention and control. As a typical case, the cerium hydride growth on the surface of cerium can be simulated as an application example of this model. This numerical simulation shows the evolution and interrelationship of the above three physical fields. In the early stage of hydride growth, the passive layer has a barrier effect on the hydrogen diffusion from the surface to the inside. And this effect is reflected in this simulation. The passive layer is destroyed in the later stage of hydride growth, which increases the rate of hydrogen diffusion and growth rate of hydride. In addition, the present work also explores the interaction of two adjacent cerium hydride grains, which reveals their mutual influence during growth process.

Method
The hydride corrosion kinetics on metal surfaces involves a variety of mesoscopic physical factors, including the interfacial reaction achieved by interfacial movement, the hydrogen transport driven by diffusion potential and the internal stress caused by hydride precipitation. Multi-physics coupling based on kinetic equations can be used to simulate these effects and present the kinetic characteristics. In the present study, three physical fields are coupled to construct the hydride corrosion model: multiphase field, hydrogen concentration field and elastic strain energy field. These fields are necessary for tracking the movement of the diffuse interfaces, calculating the diffusion/redistribution of hydrogen, and revealing the effect of strain energy field on the hydride growth.
For the multiphase field, three order parameters, which can be marked as f , p f h and f , m are introduced to represent the volume fractions of passive film, hydride and metal phase, respectively. As expressed by equation (1). the governing equation of the order parameters is implemented to calculate distribution evolution of three phases and their interfaces. As kinetics of diffuse interface, this equation is an extension of the Allen-Cahn equation [39] in the multiphase case, and the Gibbs-Thomson effect [40] is included implicitly, which resolves the interfacial movement into two factors: curvature effect and driving force [41][42][43]. The curvature effect is related to the first three terms of equation (1), which is dependent on the local morphology of the hydride interface, and the driving force is reflected in the last term of equation ( [44][45][46], while the latter, D b a  G el can be calculated from the local elastic strain energy density f el [47,48]. In addition, because the sum of volume fractions f , p f h and f m is 1, a constraint equation, equation (3) is introduced to ensure self-consistency.
where α, β, θ indicate the permutation of different phases and n is the quantity of total phases. Concretely, = n 3 and α, β, θ represent the permutation of passive film (p), hydride (h) and metal phase (m).
f is the order parameter of multiphase field.
M is the kinetic coefficient of multiphase field. e is the gradient energy coefficient.
H is the barrier energy of double well function.
is the driving force of phase transformation from β to α. D b a  G ch is the chemical driving force of phase transformation from β to α. D b a  G el is the strain-energy driving force. Based on the phase transformation from β to α. If the volume of α is larger than β, el el where f el is the strain energy density. For the hydrogen concentration field, the variable can be marked as C, which can be expressed in the mole fraction of hydrogen. As a conserved variable, the governing equation of hydrogen concentration field is the expressed by equation (4), which follow the basic form of diffusion equation. This equation not only governs the hydrogen transport in passive film, hydride and metal phase, but also achieves the hydriding reaction simulation by hydrogen redistribution in the interfacial region. The diffusion coefficients D in these three phases can be different. Especially, the passive film generally acts as a barrier to hydrogen transport, so the diffusion coefficient in passive film is often relatively small [28].
C is the concentration of hydrogen. α indicates the permutation of passive film (p), hydride (h) and metal phase (m). a c is the phase composition of hydrogen in a phase, which is derived from the quasi-equilibrium condition of KKS model. a D is the hydrogen diffusion coefficient in a phase. According to the approach of KKS model, the hydrogen concentration can be separated into fictitious phase compositions c , p c , h c m in passive film, hydride and metal, respectively. The introduction of phase compositions can not only be used to solve equation (4), but also be used to express the chemical driving force D b a  G ch in equation (2). The hydrogen diffusion potentials in passive film, hydride and metal can be marked as ( ) which can be expressed as the first derivative of free energy density functions of these three phases Based on the quasi-equilibrium assumption of KKS model, these fictitious phase compositions c , p c , h c m need to satisfy both diffusion equilibrium and a constraint condition [46], which can be expressed as As a system of equations, equation (5) can be solved to obtain phase compositions c , p c , h c , m and then the chemical driving force D b a  G ch in equation (2) can be calculated by these phase compositions [44], which can be expressed as, where a f and b f are the free energy density functions of a and b phases. The interaction between the hydrogen diffusion and the hydride interface movement is achieved by the chemical driving force and the hydrogen redistribution in the interfacial region: the hydrogen concentration around the hydride interface determines the chemical driving force D b a  G , ch and then the chemical driving force, the interface curvature and the strain-energy driving force together determine the velocity of hydride interface movement. On the other hand, this interfacial velocity determines the hydrogen consumption rate in the interfacial region, which affects the hydrogen concentration gradient near the interface, thereby affecting the hydrogen diffusion. Through this feedback mechanism, the hydrogen diffusion and the hydride interface movement can restrict each other to achieve a steady state.
For the strain energy field, the Khachaturyan's microelasticity theory [49,50] is used to reveal the physical effects of strain energy caused by hydride precipitation. Due to the local volume expansion, the hydride forms in which it does not quite fit. The hydride strains the matrix and is strained by it. At the hydride interface regions, the strain energy affects the hydriding reaction through strain-energy driving force D b a  G el in equation (2). Based on the Khachaturyan's method, the elastic strain energy fields can be solved from the Eshelby cycle [50], which separates the elastic strain energy density f el as two parts: f el is the elastic strain energy density.
f self and f relax are self-energy density and relaxation energy density defined by Eshelby cycle. e ij 0 and e ij are transformation strain and relaxation strain field defined by Eshelby cycle. l ijkl are the elastic moduli. s ij 0 is the transformation stress, which can be calculated by s l e = .
ij ijkl kl 0 0 According to the Eshelby cycle, the relaxation strain e ij is the sum of homogeneous strain field e ij and heterogeneous strain field e D , ij which can be expressed as, ( ) e e e = + D 10 As the uniform macroscopic strain, the homogeneous strain e ij determines the macroscopic shape deformation of whole system, which is produced by the internal stress due to the presence of hydride. Therefore, the e ij can be calculated by, where V h is the total volume of hydride, and V total is the total volume of whole system. The heterogeneous strain e D ij can be expressed as the first derivative of heterogeneous displacement u, which can be expressed as, Based on Eshelby cycle, the total internal stress field .. can be calculated by multiplying the elastic moduli and strain, which can be expressed as, For the calculation of heterogeneous strain field e D , ij the total internal stress s ij total expressed by equation (13) needs follow the conservation of momentum, which can be expressed as, Substituting equation (12) into equation (14), the heterogeneous displacement u can be solved, and then the heterogeneous strain e D ij can be calculated. According to equations (10), (11) and the heterogeneous strain e D , ij the relaxation strain e ij can be obtained. On this basis, the calculation of elastic strain energy density f el can be carried out by equations (7), (8) and (9), and this energy density field is not only dependent on the boundary condition, but more importantly, dependent on the shape and volume of hydride. There is an interaction between the strain energy density field and the interfacial movement of the hydride: the interfacial movement changes the shape and volume of the hydride, which in turn changes the strain energy distribution. As a part of the driving force, strain energy also affects the velocity of hydride interfacial movement, and this mechanism is reflected by equations (1) and (2).
The multi-physical simulation of the hydride growth on metallic surface is to solve a set of time-dependent equations consisting of equations (1) and (4). We apply a numerical scheme based on explicit finite difference algorithm to calculate these time-dependent equations. Accordingly, the iteration at each time step also needs the calculation of driving force D b a  G in equation (1) and phase compositions c , p c , h c m in equation (4), which can be obtained indirectly by solving equation (5) and equation (13). We apply the parabolic approximation scheme [46] to solve equation (5), and the fast Fourier transform method to solve equation (13). A Fortran program is designed for the calculation of these equations and a FFTW (Fastest Fourier Transform in the West) interface based on C programming is specially designed to solve equation (13). In addition, the distributed memory parallel computing based on MPI (Message Passing Interface) library is used to improve the efficiency of the calculation.
Due to the typical characteristics of hydride corrosion on cerium surface, we implement the simulation of cerium hydride growth as an example to test this kinetic model. Cerium exposed to hydrogen-containing atmosphere can produce cerium hydride pitting corrosion on its surface [51]. And the growth process of cerium hydride is similar to that of other rare earth metals and transition metals, so it can be used as an example to explore the hydride corrosion kinetics on metallic surface. Four allotropic forms of cerium are known to exist at standard pressure, which are δ-Ce (bcc), γ-Ce (fcc), β-Ce (hcp) and α-Ce (fcc). The γ-Ce that is stable at 300°C is set as the metal phase in this simulation. The surface of cerium is covered with a passive film composed of cerium oxide CeO 2 and Ce 2 O 3 [52,53], which acts as a diffusion barrier to hydrogen. The cerium hydride phase can be suggested as a heterogeneous solid solution of non-stoichiometric hydrides CeH x . When the hydride forms at the cerium/hydride interface, the H/Ce ratio of cerium hydride can fluctuate around 2:1 according to the local hydrogen supply. This hydriding reaction can be expressed as Ce+x[H]=CeH x . If the hydrogen concentration at the reaction interface is low, H/Ce ratio x may be less than 2, otherwise, it may be greater than 2, or even x=3 [54]. The ideal cerium hydride CeH 2 crystal [15] has the same structure as CaF 2 crystal, which is based on the fcc lattice and a basis with Ce atom at ( ) 0, 0, 0 and H atoms at ( ) , , The cerium hydride CeH x has a non-stoichiometric structure, which can be regarded as CeH 2 lattice with some H atoms missing (x<2), or CeH 2 lattice with more H atoms inserted in the gaps (x>2).
The initial distribution of order parameters f , p f h and f m is shown by figure 1. In order to simulate the interaction between the hydride, metal and the passive film during hydriding process, we place the crystal nucleus of the cerium hydride near the passive film/cerium metal interface. For the effectiveness of the multiphase-field model, the thickness of the passive film must be greater than the thickness of the diffuse interface. Therefore, the artificial thickness of the passive film is much greater than the actual thickness. The problem caused by this artificial thickness can be corrected by modifying the diffusion coefficient of the passive film phase. The boundary conditions of hydrogen concentration field are the key to determining the source of hydrogen. As shown in figure 1, periodic boundary conditions are used on the left and right boundaries, and Dirichlet boundary conditions are used on the upper and lower boundaries: the hydrogen mole fraction of both upper and lower boundaries is set to 0.05. The initial condition of the hydrogen concentration field is consistent with its boundary conditions. The initial hydrogen mole fraction in passive film and cerium metal is 0.05, and 0.667 in the cerium hydride. In addition, the physical parameters used in the simulation of Ce-H system are listed in table 1. Especially, the elastic moduli l ijkl used in the present simulation study are isotropic because the cerium surrounding the hydride is usually polycrystalline, and the complex simulation case of polycrystalline cerium and grain boundary diffusion will be studied in the future.
The free energy density functions of γ-Ce metal, cerium hydride and passive film can be marked as f , m f h and f , p respectively. As shown in figure 2, f m and f h are represented by the blue and red curves, which are calculated from the pressure-temperature-composition data at 300°C [16,54]. As the barrier of hydrogen diffusion, passive film phase only reacts with hydride. Therefore the free energy density of passive film f p can be set as same as that of cerium phase f . m Based on the quasi-equilibrium assumption of KKS model, the slopes of the parallel tangent lines in figure 2 reflect the hydrogen diffusion potentials in different phases, which are equal according to equation (5). The phase compositions c , m c h are the concentrations at the tangent points, which are dependent on the local order parameters and hydrogen concentration. The free energy density curve of the cerium hydride f h is steep around the extreme point, which makes the phase composition in cerium hydride always close to the extreme point / = c 2 3. Therefore, the hydrogen concentration in the cerium hydride obtained by the simulation is relatively stable and close to the hydrogen mole fraction of CeH 2 . The chemical driving force expressed by equation (6) can be regarded as the drop between two parallel tangent lines. In Figure 2. The free energy density of Ce-H system at 300°C: the blue solid curve represents cerium metal (γ-Ce), and the red solid curve represents cerium hydride. The blue and red dotted lines are the parallel tangent lines of free energy density curves of cerium metal and cerium hydride, respectively. D  G m h ch is the chemical driving force from cerium metal to cerium hydride.  In the multiphase-field maps, the brown, orange and red regions represent the passive film, cerium metal and cerium hydride phase, respectively. In the hydrogen concentration maps, continuous hue change from blue to red indicates varying of hydrogen mole fraction from 0 to 0.667, respectively. In order to show the gradient in the low concentration region, logarithm of mole fraction, ln(C) is used to correspond to the hue. In the strain energy density maps, continuous hue change from blue to red indicates varying of energy density from 0 to 1.7×10 6 kJ m −3 , respectively.).

Results and discussion
The simulation results reveal the evolution of multiphase field, hydrogen concentration field and strain energy field. As shown in figures 3(a)-(c), cerium hydride grows up near the passive film/cerium interface and gradually destroys the overlying passive film, which makes the hydride directly contact the upper boundary. After breaking through the passive film, the morphology of the hydride evolves into a semi-ellipsoid. Since the hydrogen diffusion coefficient of the passive film is smaller than that of hydride, the destruction of the passive film promotes the hydrogen transport from the upper boundary to the inside, thereby promoting the growth of cerium hydride. The hydrogen concentration evolution maps shown in figures 3(d)-(f) present the maximum gradient of hydrogen distribution appears at the regions of passive film/hydride interface and cerium/hydride interface, and hydrogen is redistributed in these interfacial regions to achieve the equilibrium of diffusion potential, which can be expressed by equation (5). The hydrogen mole fraction in the cerium hydride is relatively stable, which is slightly less than the hydrogen proportion in pure CeH 2 . The concentration of dissolved hydrogen in the cerium metal and passive film is much lower than that in cerium hydride, and the lowest concentration appears around the hydride, where the hydrogen diffuses to the hydride interface for hydriding reaction. Figures 3(g)-(i) present the distribution and evolution of strain energy density, and this strain energy comes from the local volume expansion caused by the formation of cerium hydride. The high-density region of strain energy is mainly distributed around the hydride interface, especially the highest strain energy density appears at the passive film/cerium/hydride three-phase junction. According to previous experimental results [2,3], the strain energy may not only cause cracking of the two brittle phases, cerium hydride and passive film, but also squeeze the cerium metal around the hydride to cause plastic bulge.  showing the position of the cerium hydride interface. In the early stages, the hydride is covered by passive film, and the growth of hydride is driven by chemical driving force but restricted by three factors: the passive-film barrier to hydrogen transport, the internal pressure caused by small radius of curvature, and the internal stress caused by the formation of hydride. After the hydride grows up and breaks through the passive film, several changes are conducive to the growth of hydride: not only hydrogen diffusion on the upper boundary becomes faster, but the radius of curvature of the cerium/hydride interface becomes larger, which makes the internal pressure caused by the curvature become smaller. Figure 4(b) shows the time-dependent changes in damage diameter of passive film and the maximum height of cerium hydride. When the evolution time is 200 s, the cerium hydride breaks through the passive film, after which the d(t) and h(t) curves become steeper which indicates that the growth of the hydride becomes faster due to the improvement of hydrogen diffusion. As the hydride grows up, the hydrogen transport distance from the upper boundary to the cerium/hydride boundary becomes longer, which makes hydrogen transport and interfacial movement gradually reach a steady state, so the curves d(t) and h(t) present a linear trend in the later stage. Figure 5(a) presents the contour lines of hydrogen concentration when the evolution time is 1800 s, and the diffusion direction of hydrogen is perpendicular to these contour lines and points to the direction of decreasing concentration. Therefore, both the hydrogen at the upper boundary and the hydrogen dissolved in cerium are driven to the hydride interface by the concentration gradient, where the process of interfacial movement and hydrogen redistribution takes place. The hydrogen concentration gradient in the cerium hydride is relatively small, which indicates the stability of the H/Ce ratio in cerium hydride. According to the free energy density curve of cerium hydride and its phase composition shown in figure 2, the hydrogen concentration in the cerium hydride is close to hydrogen mole fraction of CeH 2 , which is about 0.667. Although the ideal H/Ce ratio of CeH 2 is 2, the H/Ce ratio achieved by the simulation using the diffusion equation is slightly lower than the ideal ratio, which also makes the front edge of the interface deviate from the equilibrium state to ensure that the driving force is not zero. In the actual process, the hydride growth is controlled by hydrogen diffusion, and the H/Ce ratio of the hydride is usually slightly lower than 2. Figure 5(b) shows hydrogen concentration distribution along the center line of the hydride at different times. The region with high hydrogen concentration is cerium hydride, and the region with low concentration is cerium metal, and the largest gradient appears at the cerium/hydride interface, which is the result of hydrogen redistribution at the interfacial region. As the hydride grows up, the local radius of curvature of the cerium/hydride interface becomes larger and larger, a larger concentration gradient is required to generate a larger flow density to supply hydrogen. This trend is also revealed by figure 5(b): the hydrogen concentration gradient near the cerium/hydride interface becomes larger over time.
As shown in figure 6(a), the contour lines present the distribution of strain energy density when the evolution time is 1800 s. And the enlarged view shows the white box region where both strain energy and strain energy gradient are the largest in the field. As shown in figure 6(b), the comparison of two sets of simulations reveals the restrictive effect of strain energy on hydride growth: one set of simulation is coupled with strain energy, and the other set is uncoupled with strain energy. The coupling of strain energy and multiphase field is achieved through the strain-energy driving force D b a  G el in equation (2). Under the same growth time, the simulation uncoupled with strain energy results in a larger size of cerium hydride than the simulation coupled with strain energy. Therefore, the comparison in figure 6(b) indicates that the hydride grows faster without strain-energy driving force D b a  G .
el This strain energy is produced by the volume expansion caused by the growth of the cerium hydride, which in turn restricts the growth of the cerium hydride. Figures 6(c) and (d) show the changes in cross-sectional area and width/height ratio of cerium hydride over time, and the blue and red curves are based on the hydride growth simulations coupled with and uncoupled with strain energy, respectively. The comparison between these two simulation cases present that, not only the crosssectional area, but also the shape of the cerium hydride can be affected by the strain energy. As shown in figure 6(c), regardless of coupled with or uncoupled with strain energy, the cross-sectional area of the cerium hydride increases as the hydride grows, and the cross-sectional area obtained by hydride growth uncoupled with strain energy is always larger than the case coupled with strain energy. As shown in figure 6(d), the timedependent curves present the changes in width/height ratio of cerium hydride. In the initial stage, the spherical cerium hydride nucleus preferentially grows perpendicular to the passive film/cerium interface, and the width/ height ratio decreases from 1.0. Then the cerium hydride preferentially grows parallel to the passive film/cerium interface and evolves from a spherical shape to a semi-ellipsoid. During this process, width/height ratio of cerium hydride keeps increasing and finally stabilizes. It can be seen from the comparison of the simulation cases coupled/uncoupled with strain energy that the width/height ratios of hydride obtained from both two simulation cases are basically the same from the beginning to 400 s. And in the later stage of hydride growth, the cerium hydride coupled with strain energy has a larger width/height ratio than that uncoupled with strain energy, which indicates that the effect of strain-energy driving force D b a  G el makes the cerium hydride flatter. In order to reveal the interaction of cerium hydride grains growth on the surface of cerium, the growth simulation of two adjacent hydride grains is implemented to present the evolution of multiphase field, hydrogen concentration field and strain energy density field. As shown in figures 7(a)-(c), the morphology evolution of each phases is presented by the multiphase field, which show the growth and collision of two hydride grains, and the destruction of the passive film. The evolution of hydrogen concentration field in figures 7(d)-(f) reveals that the hydrogen dissolved in cerium forms an obvious diffusion gradient around the cerium hydride, and the its concentration is the lowest at the junction of the hydride grains. And the evolution of strain energy density field in figures 7(g)-(i) presents that the strain energy density near the hydride interface is larger. Before the two grains collide, there is the maximum strain energy density around the passive film/cerium/hydride three-phase junction, and after the collision, the maximum strain energy density is located at the junction of the two hydride In the multiphase-field maps, the brown, orange and red regions represent the passive film, cerium metal and cerium hydride phase, respectively. In the hydrogen concentration maps, continuous hue change from blue to red indicates varying of hydrogen mole fraction from 0 to 0.667, respectively. In order to show the gradient in the low concentration region, logarithm of mole fraction, ln(C) is used to correspond to the hue. In the strain energy density maps, continuous hue change from blue to red indicates varying of energy density from 0 to 4.0×10 6 kJ m −3 , respectively.).
grains. In particular, figure 7(h) shows that when two hydride grains are about to collide with each other, the cerium region between these two hydride grains has a very large strain energy density, which may cause local plastic deformation.
The morphology evolution of cerium hydride can also be visualized by figure 8. (a): two spherical hydride nuclei first grow independently into two semi-ellipsoids, and then become a half hourglass after collision. The distance between the initial hydride nuclei is defined as D, and figure 8(b) shows the time-dependent changes in the cross-sectional area of cerium hydride grains, in which different colors and line types represent the simulation results obtained by ten groups of different initial distance D. The cross-sectional area of hydride grains not only increases gradually with time, but also increases faster and faster. And the two black curves limit the range of all other colored curves: the black solid line on the lower side is obtained by the simulation case of a single hydride grain growth, and the black dotted line on the upper side is obtained by the growth simulation of two completely independent hydride grains. Therefore, these two black curves are the asymptotes of the two extreme cases: = D 0 and = ¥ D . Since two cerium hydride grains with a certain distance have little mutual influence in the initial stage of growth, the colored curves generally coincide with the black dotted line in the initial stage.

Conclusions
In short, the hydride corrosion kinetics on metallic surface is studied by a multiphase-field model, which takes account of the multiphase field, hydrogen concentration field and the strain energy density field. According to the governing equations of each field, this model can be used to simulate the multiphase interaction of passive film, hydride and metal, the diffusion and redistribution of hydrogen, and the effect of elastic strain energy. Using cerium hydride growth on the cerium surface as an application example of the model, the numerical simulation reveals that the morphology of the cerium hydride evolves from the initial (b) presents the changes in cross-sectional area of cerium hydride over time, in which different cases of center distance D are marked with different colors. The black solid line represents the growth of a single hydride grain, and the black dotted line represents the growth of two hydride grains that do not interact (center distance is infinite).