Phase-field modelling: effect of an interface crack on precipitation kinetics in a multi-phase microstructure

Premature failures in metals can arise from the local reduction of the fracture toughness when brittle phases precipitate. Precipitation can be enhanced at the grain and phase boundaries and be promoted by stress concentration causing a shift of the terminal solid solubility. This paper provides the description of a model to predict stress-induced precipitation along phase interfaces in one-phase and two-phase metals. A phase-field approach is employed to describe the microstructural evolution. The combination between the system expansion caused by phase transformation, the stress field and the energy of the phase boundary is included in the model as the driving force for precipitate growth. In this study, the stress induced by an opening interface crack is modelled through the use of linear elastic fracture mechanics and the phase boundary energy by a single parameter in the Landau potential. The results of the simulations for a hydrogenated (α+β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha +\beta $$\end{document}) titanium alloy display the formation of a precipitate, which overall decelerates in time. Outside the phase boundary, the precipitate mainly grows by following the isostress contours. In the phase boundary, the hydride grows faster and is elongated. Between the phase boundary and its surrounding, the matrix/hydride interface is smoothened. The present approach allows capturing crack-induced precipitation at phase interfaces with numerical efficiency by solving one equation only. The present model can be applied to other multi-phase metals and precipitates through the use of their physical properties and can also contribute to the efficiency of multi-scale crack propagation schemes.


Introduction
The loss of function and usefulness of a component can be related to the deterioration of mechanical and physical properties of structures, depending on the environment and the stress conditions in which it operates. This degradation of materials can be the result of damage mechanisms, such as stress corrosion cracking and hydrogen embrittlement, for which the combination of a corrosive environment and stress can enhance the propagation of cracks (Jones 1992). Hydrogen embrittlement typically occurs in metals operating in hydrogen-rich environment, e.g. in nuclear power reactors or rocket engines. Diverse damages induced by the interaction of hydrogen with the materials can be observed such as hydrogen attack, blistering and hydride formation (Cramer 2003). The latter phenomenon is characterized by the precipitation of non-metallic phases, hydrides, usually more brittle than the metallic matrix. The relatively low fracture toughness of the hydrides can contribute to the reduction of the load bearing capacities of metals such as titanium-and zirconium-based alloys (Coleman and Hardie 1966;Chen et al. 2004;Coleman et al. 2009;Luo et al. 2006). The delayed hydride cracking (DHC) mechanism, which is observed in these materials as they operate in hydrogen-rich environment under stress, is characterized by a stepwise sequence of complex phenomena involving diffusion of hydrogen, hydride precipitation, material deformation and crack propagation (Northwood and Kosasih 1983;Singh et al. 2004;Coleman 2003;Puls 2012). Hydrogen diffusion in the material is generally driven by positive stress gradients, which can be induced by the presence of stress concentrators such as dislocations, cracks, notches or residual stresses (Birnbaum 1976;Takano and Suzuki 1974;Grossbeck and Birnbaum 1977;Shih et al. 1988;Cann and Sexton 1980). When the solubility limit is exceeded hydrides can precipitate and expand-a material swelling takes place in the reacting area resulting in a volume increase of the material (Coleman 2003;Varias and Massih 2002). In regions of the phase diagram, where only single-phase solid solutions are stable in stress-free conditions, hydride precipitation can be triggered by the presence of stresses (Varias and Massih 2002;Birnbaum 1984;Allen and Vander Sande 1978), inducing a shift of the terminal solid solubility. In fact, the solubility limit of hydrogen in the metal depends not only on temperature but also on the level of pressure/stress. Some titanium alloys such as Ti-6Al-4V and Ti-0.3Mo-0.8Ni can present two phases, α (hexagonal close-packed-HCP) and β (body-centered cubic-BCC) and be subjected to hydride forming (Liu et al. 2018;Sun et al. 2015). The crystallography of the δ hydride phase, considered in this paper, is face-centered cubic-FCC. Hydride formation commonly occurs preferentially in phase α because of its low solubility and diffusivity in hydrogen compared to that of phase β (Banerjee and Arunachalam 1981;Manchester 2000). Another reason for this is the difficulty for BCC structures to transform into FCC (Liu et al. 2018). In these materials, grain and phase boundaries, more energetic than the rest of the material, have been observed to be preferential sites for hydride formation (Banerjee and Mukhopadhyay 2007;Coleman 2003;Tal-Gutelmacher and Eliezer 2004). For instance, hydride regions are observed to precip-itate as large lamellar at α/α interfaces and to form as slender plates at the α/β interfaces in Ti-0.3Mo-0.8Ni (Liu et al. 2018). In grain and phase boundaries, the energy required for hydride precipitation is usually lower. The presence of an externally induced local stress in proximity of a grain or phase boundary can further facilitate hydride growth.
A number of models related to second-phase formation at a flaw tip in various crystalline materials have been developed in the past years (Varias and Massih 2002;Deschamps and Bréchet 1998;Gómez-Ramírez and Pound 1973;Boulbitch and Korzhenevskii 2016b;Léonard and Desai 1998;Hin et al. 2008;Massih 2011;Bjerkén and Massih 2014;Jernkvist and Massih 2014;Jernkvist 2014;Nigro et al. 2018). The phasefield method (Provatas and Elder 2010), based on the Ginzburg-Landau theory and employed in some of the cited works, has been not only extensively used within the magnetic field (Cyrot 1973;Berger 2005; Barba-Ortega et al. 2009;Cao et al. 2013;Gonçalves et al. 2014) but also and particularly to predict microstructure evolutions in material. Some studies can be found in Chen (2002), Moelans et al. (2008), Steinbach (2009), Bair et al. (2017), Hektor et al. (2016), and Tourret et al. (2017). In the phase field theory, conserved and non-conserved phase-field parameters, also referred to as order parameters, are utilized to represent the microstructure and their time variation corresponds to the microstructure evolution (Desai and Kapral 2009). The conserved and the non-conserved phase-field variables are employed to describe diffusional and diffusionless phase transformation respectively. In phasefield approaches, solid solution and precipitates have different degrees of order, which can be represented by non-conserved phase-field variables. The concentration of solute is a typical conserved order parameter. Diffusion is commonly slower than the reordering of the microstructure and is, therefore, the limiting aspect for precipitation. Thus, the phase transformation kinetics is usually mainly driven by the diffusional process. However, in a material subjected to an external load, high stress concentration can arise in the vicinity of stress concentrators and induce a local shift of the solubility limit for a given concentration of solute as mentioned earlier. Consequently, a diffusionless phase transformation can occur in proximity of flaws. This type of situation can justify the choice of some authors to solely employ the nonconserved parameters to represent the microstructural changes. For instance, to study a quasi-static phase transformation in the process zone of a propagating crack, Boulbitch and Korzhenevskii (2016a) employ a non-conserved order parameter, which interacts with the crack-induced stress field. A parametric study of second-phase formation in presence of a crack based on a similar set-up was carried out in Nigro et al. (2018).
The present paper provides a phase-field approach in order to model second-or third-phase precipitation kinetics induced by a crack lying along a phase interface. The model makes use of a two-component and non-conserved phase field variable. The spatiotemporal evolution of the microstructure is driven by an energy minimization scheme using the time-dependent Ginzburg-Landau, TDGL, (or Allen-Cahn (Allen and Cahn 1979) equation only. The formulation includes a coupling between the phase-transformation swelling, the terminal solid solubility and the applied stress, as a driving force for phase transformation. The presence of a interface crack is modelled by using linear elastic fracture mechanics (LEFM) (Rice et al. 1990). In addition, the bulk free energy density is formulated such that the energy of the grain or phase interfaces can be controlled to promote phase transformation therein. The aforementioned model features contribute to low computational costs and are discussed later in the paper.
Since a number of titanium alloys are used in industry and interact with hydrogen through welding or direct exposition to hydrogen, e.g. turbines pumping the hydrogen from the tank to the combustion chamber of a rocket, it is relevant to apply the model on such metals. In this paper, a titanium alloy with α and α + β regions and containing cracks along α/β and α/α interfaces is considered. The concentration of hydrogen is set below the nominal terminal solid solubility of phase α. The effects of the load and the energy of the α/α and α/β interfaces on the phase transformation kinetics are studied. The impact of isotropic and anisotropic phase transformation-induced swelling is also regarded. This approach is suitable to investigate hydride formation in Ti-alloys but is also applicable to other types of precipitation, multiphase microstructure, morphologies of grain/phase boundaries, defects and loading modes.
The organization of the paper is as follows; first, the approach for modelling the growth kinetics of a precipitate in the presence of an interface crack is described in Sect. 2. Thereafter, the numerical methodology employed in the simulations is presented in Sect. 4. In Sect. 5, the parameters used in the simulations are provided and explained before presenting and interpreting the results. In Sect. 6, some features of the present model are discussed.

Model description
In this work, a study of the evolution of a multiphase system of volume V is carried out. The system matrix is considered to be initially composed of two prevailing phases: a phase α and a phase β. In specific conditions described later, a third phase, δ, can exist and be stable. The present model is formulated such that the phase evolution of the multiphase system is represented by the change of a phase field η. It is defined as a two-component vector η = η i with i ∈ {1, 2} and depends on time and space. The components of η i are real scalar functions that are related to the order of a crystal structure and, consequently, its morphology. In this work, phases α, β and δ are characterized by (η 1 , η 2 ) = (−1, −1), (η 1 , η 2 ) = (−1, 1) and (η 1 , η 2 ) = (1, η 2 ) respectively. The α/δ and β/δ interfaces are defined by (|η 1 | = 1, η 2 = −1) and (|η 1 | = 1, η 2 = 1) respectively. The α/β interface is represented by (η 1 = −1, |η 2 | = 1). Quantities, relative to phase α and phase β are respectively denoted with the superscript α and β. In order to keep these quantities continuous and differentiable through the α/β interface, an interpolation function h α β varying in η 2 is employed. In this paper, this smoothing functions is chosen as, where X α , X β are quantities related to phase α and phase β respectively, and p = 3 4 (X β − X α ) and q = 1 2 (X α + X β ). This function is formulated such that In order to represent the effect due to the presence of a third phase, an additional interpolation function, h m δ , which varies in η 1 and is defined such that h m δ (η 1 = −1) = 0 and h m δ (η 1 = 1) = 1, is utilized and is expressed as, Based on the evolution of the free energy of the system, the model is formulated such that the system evolves toward an energetic minimum, which corresponds to a stable or meta-stable state, through the change of the order parameters in time. This is achieved by using the TDGL equation expressed as where Γ i j is a diagonal 2×2-matrix, which can also be referred to as the mobility matrix. The diagonal terms Γ 11 and Γ 22 are kinetic coefficients and are respectively related to the phase transformation occurring between the phases α and β, and between matrix and hydride. The total free energy of system F can be expressed as where The term f grad stands for the gradient free energy density and accounts for the spatial variation of the order parameters in the phase interfaces of the system (Desai and Kapral 2009). It can be expressed as where g mδ = h α β (η 2 , g αδ , g βδ ). The quantities g αβ , g αδ and g βδ are material parameters related to the interfacial energy, which resides at the α/β, α/δ and β/δ interfaces respectively. These quantities affect the width of the interfaces, where the gradients in η 1 and η 2 vary significantly (Provatas and Elder 2010). For simplicity, the interfacial energies associated with the α/β, α/δ and the β/δ interfaces are assumed isotropic by choosing g αβ , g αδ and g βδ constant. The term f land accounts for the bulk free energy density and is also known as the Landau potential, which was historically formulated by Landau and Lifshitz (1980) as a phenomenological contribution to the free energy in the form of a polynomial. In this paper, the historical expression is modified to agree with the definitions of the order parameters given above and is given by where Landau potential for a αδ = 0.4, a βδ = 0.64, a αβ = 1 and s = 0.5. The color scale goes from dark blue, representing the minimum, to yellow, representing the maximum An appearance of the Landau potential is presented in Fig. 1 to have a better understanding of the different parameters used in this function and its variation with respect to the phase field variables. The parameters a αδ , a βδ and a αβ are the respective energy barrier coefficients of the α/δ, the β/δ, and the α/β transitions. The relative nucleation energy barriers are obtained by multiplying these coefficients by the positive energy constant P 0 . The term P 0 f a f b stands for the Landau potential associated to the transition between the matrix and the hydride phases. The height of the energy barrier relative to the α/δ and β/δ transformations is continuously taken into account through the α/β interface by using the interpolation function f b . The second term of the bulk free energy density accounts for the presence of the interface between the matrix phases with the height of the energy barrier P 0 (a αβ − s h m δ ). The positive scalar s is used to control the energy of the α/β interface and the associated nucleation energy barrier.
Through the use of f c , the impact of the parameter s is maximum for η 2 = 0 and is progressively attenuated as |η 2 | −→ 1. When s = 0, for an α/α or a β/β interface configuration, the energy of a grain boundary is the same as that of the rest of the system. This can be understood as a situation where there is no grain boundary. Thus, the function f c f d allows to model the α/β interfaces with higher energy compared to the bulk of the material. The bulk free energy density presented in Eq. (7) possesses minima for (η 1 , η 2 ) = (−1, −1), (−1, +1), (+1, −1) and (+1, +1) for f d > 0. When f d ≤ 0 and for (η 1 , η 2 ) = (+1, 0), one more minimum appears and becomes a global minimum as f d decreases. The minima of the system's free energy correspond to the stability or metastability of specific phases. The stability of the different phases considered in the system is examined in Sect. 3 as s and the external load are varied. The sum f st + f int designates the elastic-strain free energy density f el , which includes an external work, as described in Chen (2002). Its terms can be expressed as and where σ i j , σ A i j , ε el i j and ε i j denote the internal stress, the applied stress, the elastic strain and the total strain tensors respectively. In this paper, all material phases are assumed homogeneous, isotropic, linear elastic, and are supposed to undergo small deformations in presence of stresses. This suggests that Hooke's law is applicable in all regions of the system. For simplicity, phase δ is considered to have the same elastic constants as the parent material region where it forms from. The energy functionals F st and F int are named strain free energy and interaction free energy respectively in this paper. The former corresponds to the elasticstrain energy for a linear elastic material, and the latter is comparable to the interaction free energy in Bulbich (1992), Léonard and Desai (1998), and Massih (2011) and coupling potential energy in Li and Chen (1998). The term F int includes a coupling between the applied stress and the phase transformation-induced dilatation effect. The variation of the interaction free energy through a change of the applied stress corresponds to a modification of the terminal solid solubility (Nigro et al. 2018). The total strains ε i j are expressed as i j ) denote the eigenstrains, stress-free strains, or phasetransformation strains, which are induced by the lattice mismatch between matrix and hydride phases, and directly connected to the volume change of the transforming material. The tensors ε s,α i j and ε s,β i j represent the phase-transformation strain tensors in phase α and phase β respectively.
In this model, the elastic strains ε el i j are assumed independent of η 1 , which implies that the total strains are modified solely by the transformation strains as the Fig. 2 Geometry of the crack third phase forms and mechanical equilibrium is considered satisfied at all times. The parameter Q = C/C s reflects the effect of hydrogen concentration C and solubility limit C s = h α β (C α s , C β s ) for a solid solution in stress-free conditions. The concentration of solute is assumed constant in the whole system as diffusion is disregarded.
The present study regards a planar problem, where an opening crack lies in and along the interface α/β as displayed in Fig. 2. The system can be described through the use of a polar coordinate system (r, θ) or a Cartesian one (x, y) = (r cos θ, r sin θ). The coordinate related to the out-of-plane direction is denoted z and the position of the origin coincides with that of the crack tip. The crack direction is set along the x-axis, where y = 0, and points towards the positive values of x. In the studied configuration, the phase boundary splits the systems in two regions, i.e. two semi-infinite planes, which are dominated by phase α for y < 0 and by phase β for y > 0. The applied stress field σ A i j is chosen to be an analytical description of the stress induced by an interface crack lying between two dissimilar materials provided by LEFM. The near cracktip stress field can be expressed as Rice et al. (1990), represent the angular functions in phase α and phase β respectively and their expressions can be found in polar coordinates in Rice et al. (1990) and in Cartesian coordinates in Deng (1993). The functions Σ k, (α) i j are obtained from those of Σ k,(β) i j by replacing π by −π . The angular functions are interpolated through the phase boundary to ensure continuity of the stress field through the phase boundary. For a bi-material with an α phase and a β phase, the oscillatory parameter ω is defined as where β 0 is a Dundurs parameter (Dundurs 1969) expressed as where μ ξ = E ξ /[2(1+ν ξ )] for ξ = {α, β}. The parameters ν α and ν β are Poisson's coefficients relative to phase α and phase β respectively. The parameter K is referred to as the complex interface stress intensity factor and is formulated as where 2 a 0 is the length of the crack and σ ∞ yy and σ ∞ xy are in-plane components of the remote stress tensor. The earlier is a tensile/compressive stress normal to the crack direction and the latter is a shear stress in the direction of the crack. In this paper, for practicality, the structure is chosen to undergo a pure tensile stress at infinity, i.e. σ ∞ xy is set to 0. In the literature, it is shown that the term K r J ω oscillates as r −→ 0, predicting zones of contact or interpenetration between the crack lips (Rice 1988;Rice et al. 1990;Hutchinson and Suo 1991). The problem of the formulation validity was studied in Comninou (1977) and Comninou and Schmueser (1979). It has been argued that the use of Eq. (10) can be valid for an interpenetration or contact zone size r con sufficiently small, e.g. less than an atomic size. An estimate of r con is given in Rice (1988) assuming a small |ω| and that ψ = arg (σ ∞ yy + i σ ∞ xy ) is taken in [−π/2, π/2]. Elaborations of the estimation of r con can be found in Hutchinson and Suo (1991) and Wang and Suo (1990). When the crack lies along a grain boundary, i.e. along an α/α interface or a β/β interface, ω = 0 and the stress field results in being naturally continuous through the interface assuming the same material properties inside and outside the grain boundary. Thus, for a grain boundary, the expression of the in-plane stress field σ A i j boils down to which is that of a classical opening crack in mode I lying in an infinite plane made of an homogeneous and isotropic material. In this situation, the stress intensity factor K = σ ∞ yy √ π a 0 is real and is generally denoted or found in the literature, e.g. in Tada et al. (2000).
Additionally, the migration of the α/β phase boundary is assumed to be much slower than the transformation of the matrix into hydride, i.e. Γ 22 /Γ 11 ≈ 0. Thus, it is enough to solve Eq. (3) for (i, j) = (1, 1) only since the distribution of η 2 is considered to remain unchanged in time with respect to the evolution of η 1 . In these terms, the equation describing the phase transformation between the matrix phases into hydride can be expressed as where τ ≡ Γ 11 t. The derivative of the different considered energy densities are determined analytically and can be written as with ∂h mδ /∂η 1 = 3 (1 − η 2 1 )/4. By assuming the eigenstrain tensor diagonal in the global coordinate system and the transformation dilations isotropic, i.e.

Analysis of the model
In this section, an analytical study of model is made in order to predict the trend of the numerical results and give them a meaning. To this end, the form of the total free energy of the system modified by the variation of the grain/phase boundary energy and the stress field is examined. Simple calculations are also made through the use of classical nucleation theory to approximate the critical size of a nucleus lying ahead of the interface crack tip. This allows to predict if a simulation will result in a growth or disappearance of the precipitate. In order to estimate the phase situation the system is evolving towards, it is necessary to identify the extrema of F provided that the total free energy of the system evolves towards a minimum. With this in mind, the functional derivative of the system's total free energy with respect to η 1 is set equal to 0, as where the functional derivative of the gradient energy is neglected. This simplification implies that the value of the phase field components in a material point remains unaffected by the surrounding points. By omitting this term, singularities, discontinuities and sharp transitions/interfaces can appear as noted in Nigro et al. (2018). Thus, this analysis boils down to the examination of the Landau potential functional with respect to η 1 while being modified by the variation of s and the applied stress. The minimization of the total free energy of the system with respect to η 2 is disregarded in this study as explained in Sect. 2. Three solutions, are found by solving Eq. (19). By considering the roots of F , three phase diagrams, where the applied load and the energy of the phase/grain boundary is varied, are drawn and presented in Fig. 3. The x-axis represents the variation of the grain/phase energy and the y-axis represents the applied load. These quantities are normalized and scaled such that the phase diagrams are valid for both matrix phases, whose material properties can be different. The regions of stability of the considered phases are characterized by a minimum of the system's total free energy. A global minimum and a local minimum indicate the stability and the metastability of a given phase respectively. In region I, the energy of the system presents a global minimum for η 1 = 1, i.e. in this situation phase δ is stable and the phases of the solid solution are unstable. In region II, one local minimum is seen for η 1 = −1 and a global one is found for η 1 = 1. This means that the phase α and phase β are metastable, and phase δ is stable. By reasoning in the same manner, the matrix phases are found to be stable and phase δ to be metastable in region III. In region IV, the matrix phases are expected to be stable and phase δ unstable. Far from the phase/grain boundary, i.e. η 2 = ±1, there is no effect of the parameter s unlike for η 2 = ±1 as formulated. Along the stability line lying between region II and III, the total free energy of the system presents two equal global minima with respect to η 1 ,

Fig. 3
Phase diagrams related to the present model. The continuous black line displays the stability line, corresponding to a material state where all present phases are stable. The blue dotted/dashed and the red dashed lines are transition lines, representing the threshold between the metastability and the instability of a phase. The appearance of the system's free energy is sketched with respect to η 1 in each region of the phase diagram i.e. the matrix phases and phase δ are equally stable. This occurs, far away from the α/β interface or in the interface for s = 0, when no load is applied. In the α/β interface, this transition is linear and is characterized by It is noticed that even under compressive stresses precipitation can occur for s = 0 in the α/β interface. More generally, third-phase formation can take place spontaneously in the grain/phase boundary whose energy density is sufficiently large with respect to the interaction free energy. The two other lines indicate the limits beyond which a metastable phase becomes unstable. The equations for these lines are . The variation of the minimum of the modified Landau potential for η 1 = 1 corresponds to a shift in the terminal solid solubility of the system, promoting or hindering phase transformation accordingly to the phase diagrams presented in Fig. 3.
Initially, if the matrix phases contain a third-phase nucleus, which is sufficiently large such that the energy barrier can be overcome, and is subjected to positive stresses, then the growth of the precipitate will occur infinitely. The same reasoning can be made for a nucleus lying in the grain/phase boundary with a higher energy than the rest of the material. Consequently, the combination of a positive stress and an energetic α/β interface makes it easier for a nucleus to grow. Simple calculations based on the classical nucleation theory (Porter et al. 2009b) are performed in Appendix A to estimate the size of nucleus needed for the phase transformation to occur.
Equation (3) is such that the time change of the phase field variable is proportional to the functional derivative of the total free energy of the system. In addition, the larger the phase/grain boundary energy and/or the elastic-strain energy the lower the minimum of F for η 1 = 1. Around the minimum being lowered, the curve representing F becomes steeper, i.e. the functional derivative of the total free energy of the system results higher. Thus, an increase of the tensile applied stress and energy of the grain/phase boundary is expected to enhance the overall growth rate of the precipitate.

Numerical method
In this work, Eq. (15) is solved through the use of the finite volume method within the open-source partial differential equation solver package FiPy (Guyer et al. 2009), well-suited to study phase transformation kinetics. The solution of this equation is the spatiotemporal description of η 1 , i.e. the description of the matrix/hydride transition, within a defined computing domain. A 1.2 µm × 0.6 µm mesh composed of equally-sized square elements with a side length l is employed. The element size is chosen such that several elements are lying along the phase interface thickness. This choice is found in the literature as a requirement to capture the profile and motion of the interfaces (Qin and Bhadeshia 2010; Moelans et al. 2008). Moreover, in order to ensure numerical stability of the solutions during the simulations, the element size and time step Δτ are taken such that Δτ ≤ l 2 /(4 max{g αδ , g βδ }) (Provatas and Elder 2010). A convergence study is also performed in order to select an optimum value for l and Δτ as a compromise between precision of the results and the numerically performable character of the simulations. The width (length in x) and height (length in y) of the forming hydride are the main features, which were monitored to verify the convergence of the results. Different time step values and element sizes were used to determine that convergence had been reached by calculating the relative error over a time period character-  Table 1, over a portion of the considered domain ized by large gradients and time variation of η 1 . The selected time step is Δτ = 9.58 × 10 −10 J/m 3 giving an average relative error minor to 1%. The selected element size is l = 2.5 nm, which induces an average relative error less than 4%. In this work, all simulations are performed on a duration of 8000 Δτ .
The value of η 1 is set to be symmetric across the domain boundary by applying the condition n · ∇η 1 = 0, where n is a unit vector perpendicular to the domain limits. In order to prevent boundary effects on hydride formation in the simulations, the domain is taken to be large enough. Thus, the minimization of the energy within the domain and the mathematical validity of the model are ensured by Eq. (15) associated with the boundary conditions. The initial value of η 1 is randomly distributed all over the mesh within [−1, −0.9]. If no extra information is provided the same seeding is used for the simulations. This ensures consistency within the comparisons of the results when a given parameter is varied. In order to model the initial system configuration, i.e. two phase regions, α (η 2 = −1) and β (η 2 = −1), separated by a smooth interface, the value of η 2 is distributed throughout the mesh satisfying the relation η 2 (y) = tanh y 2 a αβ P 0 /g αβ . The latter is the steady-state solution of Eq. (3) with j = 2 for an elastically unconstrained system. The distribution of η 2 is illustrated in Fig. 4.
Behind the interface crack tip, i.e. for x < 0, it is irrelevant to model the presence of the phase/grain boundary through a variation of its energy because of the material discontinuity. Therefore, in this situation, the following conditions are imposed: s = 0, for x < 0, and s = s 0 , for x ≥ 0. In order to keep continuity of x 0 is the abscissa of the crack tip and l sub is set as a sub-atomic length.

Input parameters
The values for the parameters used in the analysis and simulations are given in Table 1 and are commented in this section. This paper aims at presenting a model for interface crack-induced precipitation within metals. In this paper, the example of hydride forming in Ti-6Al-4V (Ti64) is treated. The material properties employed in the simulations are the elastic moduli and Poisson's ratios relative to the phases of the Ti64 microstructure, i.e. phase α and phase β, which are taken in Fan (1993). In this reference, Poisson's ratio is considered the same for phase α and phase β, i.e. ν α = ν β = ν. The selected values and boundary contitions lead to |ω| ≈ 0.01766 and ψ = 0 o . According to Rice et al. (1990), it leads to a r con much smaller than an atomic size and, therefore, interpenetration or contact zone can be ignored.
Hydride growth at an α/α grain boundary and an α/β interface is investigated with the present model. A typical interface between an HCP and a BCC crystal structure is characterized by the plane relationship {110}||{0001} (Ojha and Sehitoglu 2016). This is the  Pederson et al. (2003) and, for phase β and δ, they are chosen from Shih and Birnbaum (1986) crystallographic configuration chosen in this paper for an α/β interface. The expansion components due to lattice mismatch at α/δ and β/δ interfaces are calculated by using the material data provided in In phase-field modelling, the interfacial energy between the different existing phases and the interface thickness is represented by the gradient energy coefficients and the considered energy barriers (Provatas and Elder 2010;Chen 2002;Moelans et al. 2008;Desai and Kapral 2009;Kim et al. 1999). Between a phase ζ and a phase ξ , the interfacial energy and thickness are denoted γ ζ ξ and w ζ ξ . The expression of these quantities can be obtained from the steady-state solution of Eq. (15) in 1d without any elastic-strain energy contributions as in Kim et al. (1999). For the present model, they are given in Eqs. (21)-(22), as where ζ ∈ {α, β} and ξ ∈ {δ, β} with ζ = ξ . The quantities associated with the indices ζ ξ are related to the ζ /ξ interface. The parameter α 0 is set to 2.944 such that the interface is defined for values of η 1 comprised in [−0.90, 0.90]. This study deals with a 2d problem where the parameters s 0 and F el have non-zero values in most studied cases and are expected to have an impact on γ ζ δ and w ζ δ with ζ ∈ {α, β}, which can slightly vary within the considered domain. The phase interface is usually 1 − 10 nm thick (Provatas and Elder 2010) whereas the hydrides can commonly be found between 0.1 and 1000 µm (Shih et al. 1988;Daum et al. 2009). The quantities a ζ /ξ P 0 and g ζ /ξ are calculated through Eqs. (21) and (22) and w ζ ξ is set equal to 10 nm as in Bair et al. (2016). Because of a lack of data, the values of the different interfacial energies are chosen as follows. In the literature, the interfacial energy of a semi-incoherent interface like α/β interface is found to approximately lie between 0.2 and 0.5 J/m 2 and it increases with the reduction of dislocation spacing (Smallman and Ngan 2014). The chosen value is taken within this range. The value of γ αδ is taken to be the same as that for α-Zr, another hydride forming metal, which possesses a similar crystallography as α-Ti (Bair et al. 2017). The interfacial energy γ βδ is set such that it reflects the difference in volume mismatch between hydride and the matrix phases by satisfying the relation γ βδ /γ αδ = ε . The terminal solid solubilities in hydrogen for phase α and β are taken in the Ti-H phase diagram at the eutectoid temperature, 298 • C, (Manchester 2000). The concentration of hydrogen, is set to be half the solid solubility limit for phase α.
In order to focus the study on the growth of the precipitate, a nucleus whose diameter is 2 R = 25 nm is placed at the crack tip. The nucleus size is of the order of magnitude of those in Bair et al. (2016) and, with the help of the relations given in Appendix A, R is chosen to be sufficiently large so that it does not disappear in phase α regardless of the external load and energy of the grain boundary considered in the next sections. Furthermore, the crack size is of the same order of magnitude as in Shih et al. (1988) and the external load, i.e. σ ∞ yy is arbitrary but is supposed to be realistic. For simplification notations in the rest of the paper, σ ∞ is written in lieu of σ ∞ yy .

Isotropic expansion of hydride
The results for an isotropic expansion of the system during hydride precipitation are presented in this section. The transformation-strain tensor are presumed diagonal in both matrix phases. The non-zero components for the isotropic transformation-strain tensor, denoted

Precipitation within α/α interface
The sum σ A ii , also referred to as the applied hydrostatic stress, present in Eq.(18) and residing around the crack tip in the α/α interface, is illustrated in Fig. 5. Provided the parameters given in Sect. 5.1, the hydrostatic stress increases as r decreases and is even found greater than 1 GPa ahead and in the close vicinity of the crack tip. Since the stress varies in cos(θ/2)/ √ r , the hydride formation takes place in an non-uniform stress field.
First, the formation of hydride is investigated in an homogeneous α phase, i.e. s 0 = 0.0. In order to analyse the evolution of the phase transformation, the change in the distribution of η 1 , as a marker of the presence of phase δ, is regarded and, in the present case, is depicted  Fig. 6. Initially, the distribution of the phase parameter is random in the range [−1, −0.9] and a circular nucleus with the radius R is centred at the crack tip, as explained in Sect. 5.1 and displayed in Fig. 6a. As time goes, the nucleus grows with respect to the crack direction and its changes shape. At τ = 500Δτ , it can be seen that the shape of the precipitate, depicted in Fig. 6b, resembles that of an isostress contour as in Fig. 5. In the next steps and until the end of the simulation, the δ-phase growth appears self-similar ahead of the crack tip as illustrated in Fig. 6c, d. Thus, the geometry of the hydride region is similar to that described by the isostress contours except behind the crack tip.
Regarding the case for which s 0 = 2.0, the nucleus also develops symmetrically with respect of the crack direction while changing its shape. In fact, at τ = 500 Δτ , the δ-phase region looks elongated towards the positive x as displayed in Fig. 7b. This elongation forms a tail, which starts from the lowest and highest points of the second-phase region with respect of the y-direction and converges into one point in the grain boundary. However, a portion of the α/δ interface, especially at the upper and the lower side of the precipitate, still follows the isostress contours. As phase δ continues broadening symmetrically with respect of the crack direction until the end of the simulation, the precipitate grows faster along the grain boundary for x > 0 than the rest of the hydride region as displayed in Fig. 7c, d. The development of the left hand side of the second-phase region looks the same as for the case with s 0 = 0.0. A quantitative study is now carried out by analysing the time evolution of the height and the width of the precipitate as the applied load and the energy of the α/α phase boundary are varied. The hydride width W αα and height H αα are defined as the distance along the x and the y directions respectively between the extreme points of the precipitate as illustrated in Fig. 7d. For practicality in this analysis, the material is considered to be hydride from the middle of the α/δ interface, i.e. η 1 = 0. The time evolutions of W αα and H αα are presented in Figs. 8 and 9 with respect the parameter s 0 , and the applied load respectively. The information brought from these figures is completed by the Figs. 10, 11, 12, and 13 where W αα , H αα , and their time derivativeẆ αα andḢ αα are presented with respect to s 0 and σ ∞ for τ = {1000Δτ, 8000Δτ }. The superscripts "inter" and "end" are employed to define the instants τ = 1000Δτ and τ = 8000Δτ respectively.
The parameters W αα and H αα are seen to increase with time regardless of a change in grain boundary energy or applied load as shown in Figs. 8,9,10,11,12,and 13. Nevertheless, the slope of the curves in Figs. 8 and 9 decreases with time similarly to a logarithmic or square root function. For relatively low grain boundary energy, i.e. s 0 ≤ 1.0, the hydride is larger and grows faster vertically than horizontally. The opposite trend is observed, for s 0 > 1.0 as displayed in Figs. 8 and 10. The vertical growth of the hydride is observed to be somewhat constant with respect to s 0 in the same figures. At τ = 8000 Δτ , the growth rate of the hydride height is constant with respect to s 0 . A slight increase oḟ H αα with respect to s 0 is found only at early times of the simulations, for instance at τ = 1000 Δτ as illustrated in Fig. 12. This explains why the curves representing  Fig. 8 are not exactly overlapping but have the same slope for large τ . The hydride width is seen to nonlinearly increase with s 0 in Fig. 10. The same observations are made regardingẆ αα in Fig. 12. For s 0 = 3.0, and at τ ≈ 130 Δτ and τ ≈ 480 Δτ , jumps in the hydride width value are noticed in the grain boundary as depicted in Fig. 8. These events are associated with coalescence of two hydride regions. For larger values of s 0 , the results, not shown in this paper, exhibit mul-  MPa. The width and height of the hydride region are observed to linearly increase with respect to σ ∞ as displayed in Fig. 11. These two geometric quantities are also seen to vary faster as the applied load increases as illustrated in Fig. 13. In the same figure,Ḣ αα is slightly larger thanẆ αα with respect to σ ∞ . However, it seems that this small difference tends to vanish in time sincė H inter αα −Ẇ inter αα <Ḣ end αα −Ẇ end αα . For s 0 = 0.0, the ratio H αα /W αα > 1 is maintained, i.e. the hydride region results larger in the y direction than in the x direction. This is in agreement with the observations made earlier regarding the results related to Fig. 6, i.e. it has been observed that the hydride growth follows the isostress contours, which are larger along the y-direction than along the x-direction for s 0 = 0.0. In Fig. 11, it appears that the curves representing H αα and W αα do not have the same variation at τ = 1000 Δτ and τ = 8000 Δτ . The same thing can be said for Fig. 10. Thus, the results for H αα and W αα are not scalable in time from one load case to another and/or from one value s 0 to another with a single linear coefficient. Some results not shown in this paper display the evolution of η 1 towards -1 in the whole considered domain, i.e. the nucleus disappears, for s 0 = 0.0 and σ ∞ < 30 MPa approximately. This is in agreement with the nucleation study made in Appendix A, which predicts, for σ ∞ < 30 MPa, approximate critical radii R c > 15.4 nm, which are larger than the radius of the nucleus used in the simulations. The expression given in Appendix A not only predicts that the critical radius is smaller for larger σ ∞ but also for larger s 0 . This is verified by the simulations. For example, the nucleus is numerically seen to grow for σ ∞ = 20 MPa and s 0 = 0.6 with the radius used for all the simulations presented in this paper, i.e. R = 12.5 nm. With the same conditions, the expressions in Appendix A gives R c ≈ 11.8 nm whereas, for s 0 = 0, R c ≈ 35 nm.
The analysis of the results for the cases studied until this point reflects well the fact that, in the present model, the energy of the grain boundary and the derivative of the elastic-strain energy with respect to η 1 , here through the crack-induced stress field, are the driving forces of the precipitation. In fact, the results display a hydride region growth, which follows the crackinduced isostress contours, and is faster as the energy in the grain boundary and the applied load increase. In abscence of grain boundary, i.e. s 0 = 0.0, the geometry of the second-phase region, growing self-similarly and symmetrically with respect to the x-axis, is, by formulation, directly related to the geometry of the isostress contours, described by the applied hydrostatic stress. This has also been shown in Nigro et al. (2018) for a specific set of parameters. The difference between the isostress contour geometry and the left-hand side of the precipitate can be imputed to the presence of the interfacial energy, which tends to reduce the interfaces, through the gradient free energy term. Without the latter, singularities and sharp transition can arise as mentioned in Sect. 3. In presence of a grain boundary, i.e. s 0 > 0.0, the same reasoning can be made to explain the fact that the formation of the precipitate is not confined solely in the phase boundary away from the crack tip. Instead, the α/δ interface forms a smooth transition between its separation point with the isostress contour and its position ahead of the crack tip in the grain boundary. Besides, the interfacial energy has a major role regarding the stability of a nucleus as it can be responsible for the disappearance of the second phase for too small nucleii. The deceleration of the precipitate development can mainly be explained by the form of the applied hydrostatic stress, which decreases in 1/ √ r increases as seen in Eq. (14) and illustrated in Fig. 5. As the transition front moves away from the crack tip, the minimum of system's the total free energy in η 1 = 1 increases as seen in Fig. 3, resulting in a lower driving force for phase transformation and, therefore, slower precipitation. In the analysis made in Sect. 3, which predicts an infinite growth of the second phase as long as the applied load is positive and the grain boundary energy is higher than that of the rest of the system, the Laplacian term in Eq. (15) has been neglected, and, therefore, the effect of the interfacial energy has not been fully taken into account and is believed to represent an additional contribution to the precipitation slowing down. In addition, it is possible that, after a sufficiently long time, the broadening of the precipitate can momentarily stop locally, i.e. in the area of the hydride tail tip, or reach a steady state when the interfacial energy dominates. The variation of the grain boundary energy is comparable to a variation of interaction energy except that s 0 affects only the grain boundary and is constant with respect to x. This explains why the height of the precipitate is affected for small τ but is similar for larger τ as s 0 varies. It also gives an explanation to the fact that multiple nucleii can spontaneously appear and grow along the whole grain boundary independently of the applied load for large s 0 as observed in Liu et al. (2018). For low values of s 0 , i.e. less energetic grain boundaries, the interfacial energy can be responsible for the disappearance of existing nucleii. Yet, for any value of s 0 > 0.0, hydride growth is enhanced in the grain boundary. In the results, the case for which s 0 = 3.0 is an intermediate case where the nucleii forms close to the crack tip due to the decrease of the nucleation energy barrier in the grain boundary and the intense crack-induced stress field, and coalesce. This type of coalescence may correspond to that observed in (Shih et al. 1988). For an α/α interface, the variation of σ ∞ corresponds to a variation of K I and can be translated in terms of variation in crack length 2 a 0 .
An attempt to estimate the results depicted in Figs. 8 and 9 has been made by collecting the phase-transformation front speed obtained by simulation for different constant applied loads on small domains. The height and the the width of the isostress contours are H iso = 3 r 0 √ 3/4 and W iso = 9 r 0 /8 respectively, where r 0 is the distance from the crack tip to the isostress contours along the x-axis for x > 0. These geometric lengths are considered to correspond to the hydride geometry, and by incrementing r 0 and, therefore, the hydrostatic stress, it has been possible to reconstruct the curves for the hydride width and height with respect to the time. The results of the simulation and those of the reconstruction displayed overlapping for some cases but discrepancies for others. The differences are imputed to the role of the interfacial energy in the phase transformation. In addition, the observations made above have demonstrated that no simple scaling factor was found between the different studied cases. The complex and non-linear nature of Eq. (15) and the non-trivial manner the interfacial energy is taken into account, through the gradient free energy and the bulk free energy, may explain the difficulty in estimating long simulation time results with simple functions and the non-scalability of the results.

Precipitation within α/β interface
In this section, the interface crack-induced hydride formation is investigated for a structure with dissimilar materials, i.e. with different material properties taken in Table 1, remotely subjected to pure tensile stress. For such study, the β phase is set to reside in the upper half part of the domain and α phase in the lower one. The same cases as in the previous section are examined here.
Some isostress contours representing the applied hydrostatic stress residing around the interface crack tip are presented in Fig. 14. Overall the stress varies in 1/ √ r as described by Eq. (10) and as in the previous section. Here, there is no symmetry with respect to the x-axis since the stress is more intense in the α phase than in the β phase as the isostress contours surrounds a larger region in phase α than in phase β. In the α/β interface the interpolation described by Eq. (1) in term of stress is visible.
First, the results obtained with isotropic transformation strains are regarded. Thereafter, the effect of anisotropy on the eigenstrains is examined.  = 20, 50, 100, 200, 500, 1000, 2000, 5000, 8000 for an α/β interface configuration with isotropic phase transformation strains. The dotted line refers to position of the α/δ and the β/δ interfaces at τ/Δτ = 8000 for for an α/β interface configuration with anisotropic eigenstrains. The dashed line represents the position of the α/δ interface at τ/Δτ = 8000 for an α/α configuration with isotropic phase transformation strains. The crack lies along the x-axis and points towards the positive values of x

Isotropic expansion of hydride at an α/β interface
The spatio-evolution of the hydride formation at an α/β interface is illustrated in Figs. 15 and 16 for s 0 = 0.0 and s 0 = 2.0 respectively. In this section, the precipitate is presumed to grow considering isotropic eigenstrains.
The formation of the hydride phase region from a pre-existing nucleus follows the same pattern in phase  = 20, 50, 100, 200, 500, 1000, 2000, 5000, 8000 for an α/β interface configuration with isotropic phase transformation strains. The dotted line refers to position of the α/δ and the β/δ interfaces at τ/Δτ = 8000 for for an α/β interface configuration with anisotropic eigenstrains. The dashed line represents the position of the α/δ interface at τ/Δτ = 8000 for a α/α configuration with isotropic phase transformation strains. The crack lies along the x-axis and points towards the positive values of x α as in the previous section. Overall, the hydride growth appears to take place in both phases but is much slower in phase β than in phase α. In the very first steps, the area covered by the nucleus is found to recede in the β phase until τ ≈ 350 Δτ . After this time point, the hydride phase starts growing in the α/β interface for η 1 > 0, but much slower than in phase α and for η 1 ≤ 0. At the end of the simulation, on the β phase side, the portion of the hydride height measured for y ≥ 0 is d end βδ ≈ 5 nm. On the α phase side, the portion of the hydride height, measured for y ≤ 0, d end αδ ≈ 160 nm. For non-zero values of s 0 and in the simulation time span, phase δ is observed to grow in phase β while the hydride growth in phase α along the y-direction is not found to be significantly affected by the change in s 0 . For example, for s 0 = 2.0, d end βδ ≈ 30 nm and d end αδ ≈ 170 nm. The precipitates gets elongated along the x-direction in phase α and in the α/β interface for s 0 > 0. The results for an α/α and an α/β interface configurations exhibit overall differences in hydride shape in phase α depending on the type of interface at τ = 8000 Δτ . For y ≤ 0, the δ-phase region appears broader in phase α for an α/β interface than for an α/α one regardless of the value of s 0 . The lower part of the hydride phase approximately follows the same isostress contour for both type of matrix interface. In the interface between matrix phases, the α/δ interface is found to reach a further point from the crack tip for an α/α interface than for an α/β one. The smoothening of the α/δ interface from its separation point with the isostress contour to its position in the interface between the matrix phases results larger for an α/α interface configuration than for an α/β one.
The trend of precipitate width W αβ and height H αβ and their time derivatives are observed to be similar to those for the previous studied case. However, the values of these parameters with respect to time, applied remote stress and energy of the phase boundary differ from those related to an α/α interface. The parameters W αβ , H αβ and their respective time derivatives are presented in Figs. 8,9,10,11,12, and 13 with respect to s 0 and applied load together with the results for the α/α interface cases. In those figures, it appears that W αβ ≥ H αβ andẆ αβ ≥Ḣ αβ for all studied cases. The precipitate height H αβ is constant with respect to s 0 but appears to be sensitive to the change of the applied load in a similar manner as H αα previously. When comparing the results for the α/α interface and the α/β one, it can be noticed that H αα > H αβ and W αα > W αβ regardless of s 0 and σ ∞ . Nevertheless, the difference between is the W αα and W αβ is relatively small with respect to σ ∞ and seems to decrease with increasing σ ∞ . Regarding the time derivatives, it appears that, overall,Ẇ αα >Ẇ αβ at the beginning of the simulation but, after some time, W αα ≈Ẇ αβ for s 0 < 1 and regardless of the remotely applied stress. By noticing that the curves forẆ αα anḋ W αβ in Fig. 12 have the same appearance for s 0 > 0 but are shifted from one another, it is possible to find the values for s 0 such that the hydride width rate is the same in both types of interface at a given time of simulation. In the present case, the shift between the curves is Δs 0 ≈ 0.5.
Overall, the development kinetics of the hydride phase for a crack lying in an α/β interface follows the same pattern as for an α/α interface. However, the hydride growth is observed to be non-symmetric and occurs much faster in phase α than in phase β. These observations are expected when considering the ratio of interfacial energy γ β δ /γ α δ ≈ 1.60 and the parameters involved in the driving force for phase-transformation related to the stress field, i.e. here, the right-hand side term in Eq.(18). Because of a difference in solubility in hydrogen, Q is one order of magnitude larger in phase α than in phase β. The stresses are also larger in phase α because of the difference in value of the Young's Modulus in the considered phases, i.e. E β /E α ≈ 0.70. However, the phase-transformation strain is greater in phase β than in α, i.e. ε s,(β) 0 /ε s,(α) 0 ≈ 1.60. Thus, the minimum of the system's total free energy for η 1 = 1 results higher in phase β than in phase α at a same distance from the crack tip. This also corresponds to a lower position in the phase diagrams in Fig. 3 for phase β than for phase α along the y-axis.
Moreover, the interfacial energy, larger in phase β and interpolated through the phase boundary, induces a lower energy of the α/β interface than that of the α/α for a same value of s 0 . The parameters involved in ∂ f el ∂η 1 , different in one phase from the other, are also interpolated through the interface. Thus, the driving force for phase transformation results lower in the α/β interface than in the α/α one. Therefore, the kinetics of the hydride development is attenuated in the α/β interface compared with that in the α/α interface regardless of the value of s 0 while it is somewhat similar in the rest of the α phase. This explains why, for y ≤ 0, the precipitate results larger and more elongated for an α/α interface configuration than for an α/β one at a given time and for the same value of s 0 .
In addition, it has been observed that the nucleus initially retracts in phase β. This is believed to come from the fact that the interfacial energy is larger in phase β than in phase α. The critical size of the nucleus results, therefore, to be larger in phase β than in phase α. For example, for s = 0.0 and σ ∞ = 50 MPa, and by using the method shown in Appendix A, R c ≈ 6 nm in phase α and R c ≈ 460 nm in phase β. At τ ≈ 350 Δτ , since the parameters composing the driving forces and the interfacial energy are interpolated throught the interface, the precipitate can grow in the α/β interface even though the critical radius calculated for phase β is not reached yet. In the cases for positive s 0 , the growth of the δ phase is facilitated not only along the phase boundary for x > 0 but also in phase β by the energy of the α/β interface, higher than that of the rest of the material, as illustrated in Fig. 16. Further, no coalescence has been observed for this type of interface and for the different studied cases.
In the literature, it is uncommon to find hydride formation forming out of phase β in an (α + β)-Ti alloy at the hydrogen concentration considered in this paper and is not expected for pressure P ≤ 1 MPa as seen in the phase diagrams given in Manchester (2000) or Sun et al. (2015). In this study, the effect of the stress on hydride formation is studied. This also means that pres-sure is varied. In this manner, the binary phase diagrams for the levels of pressure considered in the studied cases might be different and display the possibility for phase β to transform directly or indirectly, e.g. β → α → δ, into phase δ at the selected hydrogen concentration. This can be also be referred to as a stress-induced shift of terminal solid solubility.

Anisotropic expansion of hydride at an α/β interface
In this section, the results considering anisotropic expansion due to the α/δ transition are regarded. The growth of the hydride is similar to that for which isotropic eigenstrains are assumed as illustrated, for s 0 = 0.0 and s 0 = 2.0 in Figs. 15 and 16 respectively. The only difference lies in the precipitate size, which results slightly larger with anisotropic expansion at a given time as depicted in Figs. 15 and 16 at τ = 8000 Δτ . The height and the width of the precipitate increase slightly faster than for isotropic expansion. However, the hydride height is observed to develop more quickly than the hydride width for the cases with anisotropic expansion compared to those with isotropic expansion. In phase β, the hydride broadening is identical as in Sect. 5.4.1. The larger growth of the thirdphase region along the y direction can be explain by the fact that ε s,(α) x x > ε s,(α) 0 . However, in the middle of the interface, where the width of the δ-phase region is maximum, ε x x = 10.87 in the cases with anisotropic expansion and ε x x = 11.19 in that of isotropic expansion. Here, the hydride width growth is expected to be faster in the latter than in the former. Yet, the opposite result is observed. This comes from the fact that the sum σ A x x ε s x x + σ A yy ε s yy , present in Eq. (17), is slightly larger in phase α for the cases with anisotropic expansion than for the cases with isotropic expansion. Overall, the differences between the results for isotropic and anisotropic swelling are not significant within the time frame used for in the simulations. However, they might be relevant for larger time spans.

Summary of the results
Hydride growth takes place in proximity of the interface crack tip by following the isostress contours of the hydrostatic stress in case of isotropic expansion, but is enhanced in the direction of the crack when the energy of the grain/phase boundary is larger than that of the rest of the material. Hydride formation is observed to occur faster in phase α than in phase β. This is understood to be induced by the variation of the different considered quantities involved in the driving force for phase transformation, ∂ f el ∂η 1 , from one phase to the other and is symmetric in case of an α/α interface. The transition front appears to non-linearly decelerate in all directions as it moves away from the crack tip in all studied cases reflecting the decrease of the stress field with r and the action of the interfacial energy. Overall, the increase of the grain/phase boundary energy and the applied load induces faster hydride growth. The variation of the grain/phase boundary energy effects only the interface resulting in the aforementioned elongation of the precipitate. The α/δ and β/δ interfaces undergo a smoothening between their separation point with the isostress contour and their position in the interface between the matrix phases, understood to be due to the presence of interfacial energy through the gradient free energy. In phase α, the precipitate is found to grow faster in the matrix interface and results more elongated at a given time and for the same value of s 0 for an α/α interface configuration than for an α/β one. This is explained by the fact that for a same s 0 , the interfacial energy is larger in an α/β interface than in the α/α interface, and that the stress, the terminal solid solubility and the elastic constants are interpolated through the matrix interface. The height and width of the hydride region are found to be slightly larger by assuming isotropic eigenstrains than anisotropic ones for an α/β interface. However, in the time period of the simulation, these differences are considered insignificant. Initially, the part of the nucleus, which lies in phase β, disappears before growing again at a certain time. Coalescence of several hydride regions are observed for a relatively large energy of an α/α interface but has not been noticed for an α/β interface.

Further discussions
The approach described in this paper presents the capacity to model stress-induced precipitation in proximity of an opening interface crack lying between two phases, e.g. a monomaterial with a polycrystalline microstructure and bimaterials. The model is formulated such as the presence of the grain/phase interface, the anisotropic expansion due to phase transformation, and the effect of the solute concentration and termi-nal solid solubility are included. Thus, precipitation at grain/phase boundary, as observed in Liu et al. 2018;Banerjee and Mukhopadhyay 2007) and in the vicinity of the crack tip as seen in Shih et al. (1988), is captured. The implementation of the aforementionned features in the present approach represents major improvements to the previous model described in Nigro et al. (2018). The present model is also capable to represent precipitation for other kind of flaws, such as dislocations, by using analytical expressions as in Bjerkén and Massih (2014). Crack-induced precipitation and second/thirdphase formation at grain/phase boundary can also be captured separately.
One benefit of this approach is that most parameters of the present model are physical and measurable. They can, therefore, be calibrated through experiments. The artificial parameter s 0 can also be fitted through the use of experimental data. This parameter is a constant in the present study but could be taken as a function of various features of the interface such as the thickness and degree of anisotropy (Provatas and Elder 2010). In addition, in some cases, e.g. for δ-hydride formation in a α-Zr matrix, the interfacial energy is found to be strongly anisotropic (Han et al. 2019). For more realistic results when modeling such phase transformation cases, it is necessary to make the gradient free energy of the present model anisotropic by, for example, making the parameters g αβ , g αδ and g βδ function of the angle between the direction normal to the interface and a reference axis (Provatas and Elder 2010).
Moreover, for known stress field distributions through an analytical expression or data, only one equation, the TDGL one, is necessary to account for the microstructure evolution. Here, the crack-induced stress field is implicitly represented through the use of LEFM. Thus, the computational resources are reduced compared to a model where a crack is explicitly represented and, consequently, where the mechanical equations need to be solved numerically. In other situations more complex than the one treated in the paper, the mechanical equilibrium might have to be solved partially or wholly at all time, reducing drastically the numerical efficiency of the present methodology.
The stress-based driving force employed in the model includes the effet of the ratio between the concentration of solute and the terminal solid solubility. The latter is believed to have an impact on nucleation and the precipitation rate as highlighted in (Liu et al. 2018;Banerjee and Mukhopadhyay 2007) for Ti-alloys and is shifted in presence of stress allowing phase transformation for lower concentration of solute than in stress-free conditions. In this pilot project, the concentration of solute is assumed constant and uniform in the two-phase system. The results presented in this paper provide an idea of the potentialities of the model. In order to obtain more realistic results it might be useful to map the distribution of the solute concentration over the system and utilize it as an input to the simulations. As seen in Sect. 3, the model is formulated such as the precipitate may grow infinitely in a material subjected to an applied load. This can happen on unreasonable orders of magnitude compared to the lifetime of a component, its size or even its embrittlement process, which can occur on smaller time scales such as DHC. In the context of DHC, the diffusionless phase transformation studied in this paper is believed to take place as a first step of the damaging process, i.e. from the time when the structure is subjected to a load, phase transformation can take place in proximity of a crack for a low hydrogen concentration. Thereafter, other processes, acting on larger time scales, occur and are combined to crystal reordering such as diffusion and propagation of cracks through the hydride-rich paths. Thus, for simulations on larger time scales, diffusion can be taken into account by solving the Cahn-Hilliard equation while coupled to the TDGL one (Provatas and Elder 2010;Moelans et al. 2008).
The information provided from this paper aims at studying precipitation kinetics induced by the coupling between an external stress field, transformation strains and terminal solid solubility regardless of the different phases of the system. As a first approximation, the material properties of the precipitate are assumed the same as the matrix phase, which it develops from. Nevertheless, the difference between the phases are taken into account through their respective eigenstrains, solid solubility and interfacial energy. Local deformations can arise from heterogeneities (Chen 2002;Moelans et al. 2008;Li and Chen 1998) and can be taken into account to make the results more realistic. However, this implies that the mechanical equilibrium is solved numerically although the crack-induced stress field can still be present implicitly. In the example presented in this paper, as heterogeneities are not represented, the hydride region arising from the crack tip can be seen as a cluster of hydrides instead of a single second/thirdphase region. The distribution of hydride phase around a notch in a Zr-2.5Nb cantilever found in Ma et al. (2006) displays a geometry, which resembles that found with the present approach for an α/α interface with no increase of the energy in the grain boundary. This suggests that the present model can be used as a fairly good approximation in terms of stress-induced hydride precipitation. The inhomogeneities induced by the appearance of a precipitate with different material properties can be addressed in the next steps of the model improvement. Additionally, the study was carried out by making a simplification consisting in considering all phases of the system isotropic in terms of elastic constants. In order to have a more realistic model, the anisotropic stiffness tensor can be employed and treated as in Nigro et al. (2017). Thus, the anisotropy due to the crystallography of an (α +β)-Ti alloys, i.e. hexagonal closedpacked for phase α (Tromans 2011) and body-centred cubic for the phase β, can be captured.
Plastic deformation, as a consequence of the presence of dislocations, is involved in the total free energy of the system as it can influence precipitation (Porter et al. 2009a). In this study, by using the expression by Irwin, the crack tip plastic zone size is estimated to be of the order of magnitude of the initial hydride nucleus size considering the yield stress of Ti and the chosen loads. Thus, as a first approximation, the impact of the non-linear zone is considered negligible for the subsequent second-phase broadening. Nevertheless, the model is not limited to LEFM as other stress fields can be superposed. For example, plasticity can be included in the present approach by simply using the stress fields induced by dislocations as in Bjerkén and Massih (2014).
Interpolation functions are used for the stress field, material properties and the transformation strains through the interface in order to keep continuity of these parameters in the system. These functions are such that the value of a parameter in the middle of the interface is the mean of its values in each matrix phases. It can be useful to use ab-initio calculation data related to the variation of these parameters through the phase interface so that more suitable interpolation functions can be used.
The present approach is fairly efficient on the computing point of view. Nevertheless, some existing techniques can be used to optimize the computing resources further. For example, adaptive meshing alogrithms as in Bair et al. (2017) or the increase of the artificial interface width together with the decrease of the element size, e.g. through the use of thin-interface asymptotics (Karma and Rappel 1998), can be employed. Nevertheless, when using adaptive meshing the algorithms usually reduce the elements size in the region of the interfaces where the phase variable gradients are the largest. This suggests that the interfaces are tracked down as in the Stefan problem (Provatas and Elder 2010) and, therefore, the numerical efficiency of the phase field method can be cut down.
Over the years, phase field models have been applied to crack propagation, (Kiendl et al. 2016;Schneider et al. 2016;Shanthraj et al. 2016;Spatschek et al. 2006Spatschek et al. , 2007Spatschek et al. , 2011. A few studies have recently been made in regard of phase-field models coupling fracture and phase transformation aspects without using LEFM. For instance, a PF model coupling crack propagation and martensitic transformation has been developed by Zhao et al. (2016). In addition, Wu and Lorenzis (2016) have presented a phase-field model fracture and its coupling with diffusion. These works could be considered in the further development of the present model in terms of crack propagation.
In hydride forming metals, where DHC can occur, crack propagation preferentially takes place through hydride-rich paths (Puls 2012; Shih et al. 1988). The present approach, which takes into account microstructural features such as the presence of grain and phase boundaries, which are known to promote hydride formation (Liu et al. 2018;Banerjee and Mukhopadhyay 2007), can be utilized to contribute to the efficiency of multi-scale crack propagation prediction schemes.

Conclusions
In this paper, we present a phase-field approach for modelling of stress-induced phase transformation kinetics in multi-phase microstructures including defects. The aim of this study is to model the phenomenon occurring in proximity of stress concentrators by taking into account the external stress, the phase transformation-induced expansion of the system, the terminal solid solubility in stress-free conditions, the effect of the interfacial energy and the energy of the grain/phase boundaries, with a numerically efficient approach. In particular, the flexibility of the formulation allows to represent isotropic and anisotropic eigenstrains, the effect of an interface crack in a bi-material through the use of linear elastic fracture mechanics and the energy of the grain/phase boundary through the use of a single parameter. With this configuration, only one evolution equation, the TDGL equation is solved, which contributes to the numerical efficiency of the present approach.
The results displaying δ hydride growth induced by an opening interface crack in a grain/phase boundary are obtained for an (α + β)-Ti alloy which contains a concentration of hydrogen lying below the stress-free terminal solid solubility. Precipitation appears to occur in proximity of the interface crack tip, following the isostress contours away from the grain/phase boundary. In the latter, the total free energy is larger and, therefore, hydride growth is enhanced inducing the elongation of the precipitate along this interface. The difference in material parameters in either side of the phase boundary is reflected through the lower hydride growth rate in phase β than in phase α. Coalescence of several hydride regions are observed for relatively large energy of the interface. Thus, these results reflect the possibility of the present approach to capture the dependency of the material properties on the precipitation kinetics at a microstructural level.
The present model can include different kinds of defects, multi-phase microstructures, morphologies of grain and phase boundaries and loading modes. We believe that the approach presented in this paper can contribute to the efficiency of multi-scale crack propagation prediction schemes.

Funding Open Access funding provided by Malmö University
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix A: Approximation of the critical size of a nucleus
In this section, the competition between the interfacial energies γ αδ and γ βδ , and the volumetric free energies is investigated in order to estimate the critical size of a δ-phase nucleus. The critical size of a nucleus defines the size above which the precipitate grows and below which the nucleus disappears. In this analysis, homogeneous nucleation is assumed considering a circular nucleus and isotropic material properties. The radius and critical radius of the nucleus are designated by R and R c . To simplify, two grains of the same phase separated by an interface with a thickness w αα or w ββ are regarded. By examining Eqs. (15), (16) and (18), it is noticed that the driving force for the phase transformation resides in the term of the Landau potential, −P 0 f c s h mδ , which controls the energy of the phase/grain boundary, and the coupling term −σ A i j ε s i j h mδ Q of the interaction free energy density. For h mδ = 1, the former term is active on the nucleus portion, which lies within the phase/grain boundary, and the latter is associated to the whole nucleus volume. The integral of these two terms forms the volumetric energies, which compete with the interfacial energy, integrated on the surface of the nucleus. To simplify the calculations, the parameter s = s 0 is taken to be uniform in the thickness of the interface and to be zero for x < 0, i.e. behind the crack tip. By putting the considered energy terms in equation as in the classical nucleation theory for an homogeneous nucleation (Porter et al. 2009b), an approximate critical radius can be found such that Eq. (23) is satisfied for R ≤ w ζ ζ /2 where ζ = α or ζ = β depending on the phase considered for the solid solution.
2 πγ ζ δ − P 0 π R c s 0 For R > w ζ ζ /2, an approximation of R c is found numerically by solving Eq. (24), 2 πγ ζ δ − P 0 s 0 From these expression it results that the critical radii are lower for larger σ ∞ and for larger s 0 . Physically, the increase of the parameters σ ∞ and s 0 helps decreasing the minimum of the system's total free energy and reducing the nucleation energy barrier.