Irreversible phase field models for crack growth in industrial applications: thermal stress, viscoelasticity, hydrogen embrittlement

Three new industrial applications of irreversible phase field models for crack growth are presented in this study. The phase field model for irreversible crack growth in an elastic material is derived as a unidirectional gradient flow of the Francfort–Marigo energy with the Ambrosio–Tortorelli regularization, which is known to be consistent with classic Griffith fracture theory. The obtained compact parabolic-elliptic system of PDEs including two regularization parameters for space and time enables us to simulate various kinds of crack behaviors with standard finite element software, without any geometric restriction on the crack path. We extend the irreversible phase field model to thermal cracking in solder and to cracking in a viscoelastic material, keeping the compact forms of the PDEs and the energy consistency. On the other hand, for hydrogen-assisted cracking in metal, we propose a compact phase field model focusing on a kinematic jamming effect of the hydrogen by a weak coupling approach. Several numerical experiments for these three models show not only their reasonableness and usefulness but also flexible extendability of the phase field approach in fracture mechanics. An irreversible phase field model for crack growth realizes a non-healing property and an energy gradient structure simultaneously. New compact and energy-consistent phase field models for thermal cracking in solder and cracking in a viscoelastic material are proposed. As a novel approach to hydrogen-assisted cracking, a weak coupling model with a kinetic jamming effect is proposed and numerically studied. An irreversible phase field model for crack growth realizes a non-healing property and an energy gradient structure simultaneously. New compact and energy-consistent phase field models for thermal cracking in solder and cracking in a viscoelastic material are proposed. As a novel approach to hydrogen-assisted cracking, a weak coupling model with a kinetic jamming effect is proposed and numerically studied.


Introduction
In this paper, we present new applications of an irreversible phase field model for crack growth in three industrial problems; solder cracking under thermal stress, cracking of viscoelastic materials, and metal cracking with hydrogen embrittlement (HE). Preceding those applications, we give a detailed review on the irreversible phase field model for crack growth focusing on its mathematical background such as energy consistency, and on several technical issues in the numerical implementation.
The basic phase field model is mathematically derived as a gradient flow (or a time relaxation process) of the Francfort-Marigo energy [1], which is defined as a sum of the elastic energy and the surface energy of a crack with Ambrosio-Tortorelli space regularization [2]. The irreversible growth condition of the crack is represented as a unidirectional gradient flow [3] that has properties of both the energy gradient structure and the irreversibility at the same time. The energy variational structure and a compact PDE formulation are the key issues for the further development of the crack phase field model.
Major methods of crack growth simulation in engineering are based on the computational methods such as the extended finite element method (XFEM) [4,5], rigid body spring model (RBSM) [6], discrete element method (DEM) [7,8], and particle discretization scheme-FEM (PDS-FEM) [9,10]. In those simulations, the predicted crack behaviors are highly dependent on the computational mesh and other numerical parameters or algorithms. From a mathematical point of view, a mathematically closed continuous PDE model is also preferable.
On the other hand, some phase field approaches for cracks growth have been developed recently. Cahn and Hilliard [11], in a pioneering study on the phase field model, introduced an order parameter for the modeling of phase separation phenomena of binary alloys, and the so-called Cahn-Hilliard equation was derived as a gradient flow of a free energy. Since then, such diffuse interface approaches have often been called phase field models and have been widely applied to various kind of interfacial phenomena [12][13][14].
Even in fracture mechanics, the phase field approach has been applied and recently recognized as a powerful modeling and simulation tool. Here, we give a brief and partial review of some important studies on the phase field approach for the fracture phenomena. More details and a comprehensive list of related papers can be found in those references.
In [15][16][17][18][19][20], Bourdin and his collaborators formulated their phase field model by applying the Ambrosio-Tortorelli regularization [2] to the time discrete Francfort-Marigo variational fracture model [1]. The Bourdin model consists of a global minimization problem at each time step. They also proposed an effective numerical algorithm called the backtracking method and show its applicability in various examples. The Bourdin model has been widely extended to the nonpenetration condition [21,22], dynamic fracture for a mode III crack [23], complex morphology of cracks under thermal stress [24], crack nucleation [25], etc.
On the other hand, Karma et al. [26] proposed a phase field model for crack growth analogous to the phase separation dynamics with a double-well potential, and their formulation has also been widely adopted in many applications, such as the helical crack front in dynamic fractures [27]. Using a slightly different approach, Miehe et al. proposed a rate-independent model for crack growth [28,29]. This approach has also been studied through comparisons with experiments [30] and has been generalized as a cohesive zone and softening model [31].
Concerning the mathematical justification for the Ambrosio-Tortorelli regularization, we can find proof of the convergence of the Ambrosio-Tortorelli energy (in the sense of the Γ-convergence) to an energy with sharp discontinuity [2], which is essentially the same as the Francfort-Marigo energy in the mode III crack setting. We also refer to [32] for convergence of the minimizing movement and to [33] for convergence in a nonpenetrating setting.
In [34,35], the authors proposed an irreversible phase field model for mode III crack growth in a two-dimensional (2D) plate. In the present paper, we discuss some generalizations of our irreversible phase field model and consider further extension of the model to three kinds of industrial problems.
The phase field approach provides several advantages in the mathematical modeling of crack growth.
1. This approach is consistent with classic Griffith theory [36] of the energy release rate and Irwin's theory of the stress intensity factors [37]. 2. In contrast with XFEM, this model does not require a crack path search, as the crack path selection is implicitly included in the phase field model. 3. There are no geometric restrictions on the crack shape.
For example, complicated merging and tip splitting of the crack can be naturally treated even in a threedimensional (3D) setting. 4. It is easy to use finite element numerical code to solve the phase field model without any special modification. In particular, the adaptive mesh finite element method (FEM) works very effectively, as the crack profile is very localized [34].
This paper aims to demonstrate the flexible extendability of the irreversible phase field model for crack growth through several new examples related to industrial applications. The first example is a cracking problem under thermal stress. We propose a natural extension of the phase field model based on thermoelastic stress, strain, and energy. The proposed model (4.2) is numerically examined through a solder cracking problem and is compared with some experimental pictures. Similarly, we also propose a phase field model for crack growth in a viscoelastic material (5.8). In contrast to other mathematical models [38][39][40], as shown in (5.13), our model is regarded as a gradient flow of a total energy including a natural viscoelastic energy (5.7).
The third example is a hydrogen-assisted cracking problem. A phase field model for a hydrogen-assisted crack was first proposed by [41] and intensively studied in [42], etc. Their models are based on the gradient flow approach and are consistent with the total energy. On the other hand, we derive our new phase field model (6.4) based on weak coupling with a kinetic jamming effect of hydrogen. As a consequence, the system no longer has a complete gradient structure but maintains a relatively simple form as a system of PDEs, and numerical simulation can still reproduce the hydrogen accumulation phenomenon. These new extended models presented in this paper strongly suggest wide applicability of the phase field model for crack growth in additional fracture phenomena and related industrial applications.
The outline of this paper is as follows. In Sect. 2, we derive our basic phase field model in a 2D or 3D setting as the gradient flow of the Francfort-Marigo type energy with the Ambrosio-Tortorelli regularization. We also impose the irreversible growth condition for the crack without disrupting the gradient flow structure. Several miscellaneous remarks concerning the real implementation of our model are stated in Sect. 2.4.
In Sect. 3, we give some numerical examples of crack growth in 2D and 3D settings. For the simulation, we use the P2 adaptive mesh FEM with FreeFem++ software [43].
In Sects. 4, 5 and 6, we consider a solder cracking problem under thermal stress, a crack growth problem in a viscoelastic material, and a hydrogen-assisted cracking problem for metal, respectively. In each of these sections, we derive our model and give several numerical results with a detailed discussion and comparisons with previous studies.
Finally, we provide a short discussion on the phase field modeling for crack growth in Sect. 7, and we present concluding remarks in Sect. 8.

Theory of variational fracture
Let Ω be a bounded domain in ℝ d (d = 2, 3) ; we suppose that Ω represents a reference domain of an elastic material at a reference temperature. We suppose that the boundary Γ ∶= Ω is decomposed into two piecewise smooth portions Γ D and Γ N = Γ⧵Γ D , where Γ D is a nonempty (i.e. the (d − 1)-dimensional volume of Γ D is positive) piecewise smooth portion of Γ with the Dirichlet boundary condition, and a surface traction boundary condition is assumed on Γ N .
In this paper, we denote the gradient, divergence, and Laplacian operators with respect to x = (x 1 , … , x d ) T ∈ ℝ d by ∇ , div , and δ , respectively. We denote the Lebesgue space on Ω by L p (Ω) and the ℝ d -valued Lebesgue space by L p (Ω;ℝ d ) . We also often use the Sobolev space defined by H 1 (Ω) ∶= {u ∈ L 2 (Ω); ∇u ∈ L 2 (Ω;ℝ d )} , and its trace space on the boundary Γ D is denoted by H 1 2 (Γ D ) . Please see [44,45], etc., for precise definitions of the Sobolev spaces.
We assume that, for each time t ∈ [0, T ] , a body force f (x, t) ∈ ℝ d on Ω and a surface traction q(x, t) ∈ ℝ d on Γ N are given and that a boundary displacement g( We fix t ∈ [0, T ] and omit (t) temporarily We consider a crack Σ in Ω , where Σ is a closure of a (d − 1) dimensional smooth hypersurface. We suppose that Σ ∩ Γ 1 N = � . The outward unit normal vector on Γ ∪ Σ ± is denoted by ν ∈ ℝ d , where Σ ± represent both sides of the crack Σ (see Fig. 1).
According to the theory of linear elasticity, a small displacement of the elastic material Ω with a crack Σ can be described by the following boundary value problem sym denotes the space of real-valued symmetric d × d matrices. We use the Einstein summation convention for spatial indices. The elasticity tensor C(x) = (c ijkl (x)) is where σ ∶ e = σ ij e ij and At time t ∈ [0, T ] , we denote the weak solution to (2.1) by u(t). Then, the weak solution u(t) is characterized by the following variational principle, i.e., it is obtained as a unique minimizer of the following elastic energy including the body and surface forces under a suitable boundary condition: where argmin stands for the argument of the minimum.
The surface energy of the crack Σ is defined by According to the theory of the variational fracture [1,19], the crack growth is characterized by the following total energy: In [1], Francfort and Marigo started their discussion with the following time-discrete evolution model for sharp cracks: For a fixed time increment τ > 0 , starting from a given initial crack Σ 0 , Σ k ≈ Σ(kτ) for k = 1, 2, … is determined as The crack evolution Σ(t) is formally obtained as a limit evolution of {Σ k } k as τ → 0 . It is shown that this energy variational approach is compatible with classic Griffith theory [1].

Phase field model
We apply the phase field approach [14,26] for the crack growth to avoid the numerical difficulty in the original variational model with the sharp crack setting. We introduce the following smooth function z(x, t) for x ∈ Ω to represent the approximate profile of crack Σ(t) (Fig. 2). We assume that 0 ≤ z(x, t) ≤ 1 and z(x, t) ≈ 1 around crack Σ , and that z(x, t) ≈ 0 for the other region. The damaged stress tensor is defined as The function z can be considered a damage variable that represents the damage ratio of the material in the sense of (2.8) (see Remark 6 in Sect. 2.4). We call z(x, t) a damage variable or a phase field of the crack.
We define For v ∈ H 1 (Ω;ℝ d ) and z ∈ Z , we consider a modified elastic energy, .  We remark that these regularized energies approximate (2.6) as E 1 ≈ E el , E 2 ≈ E s , and E ≈ E tot . This is called the Ambrosio-Tortorelli regularization, and it is proved that the energy E 2 (z) approximates the surface energy ∫ Σ γ * (x)ds in the sense of the Γ-convergence for some special cases [2]. Similar to the case of E tot , for ζ ∈ Z , we define where For fixed t ∈ [0, T ] and z ∈ Z , as u = u(t, z) ∈ V (g(t)) is a minimizer of E(t, v, z) among all v ∈ V (g(t)) , the first variation of the energy E(t, u, z) with respect to u has to be zero. For all v ∈ V (0) , assuming that u ∈ H 2 (Ω;ℝ d ) , we obtain This implies that the force balance equations in Ω and on Γ N are as follows: We usually suppose that z ∈ [0, 1) (see Remark 2 in Sect. 2.4). As q(t) = 0 on Γ 0 N and z = 0 on Γ 1 N , follows from (2.12). Similarly, the first variation of E * (t, z) for z is formally derived as follows. Choosing any ζ ∈ H 1 (Ω) with ζ = 0 on Γ 1 N , we suppose that Then, using (2.11) and (2.12), we have Then, we derive our crack growth model as a gradient flow of the energy E * (t, z) with respect to the damage variable z(x, t) with a boundary condition z ν = 0 on Γ⧵Γ 1 N as follows: where α [J m −3 s = Pa s, for d = 3 ] is a positive constant that is related to the time relaxation (see Remark 9 in Sect. 2.4). Hence, we have reached the following initial-boundary value problem of the phase field model for crack growth.
The second equation expresses the crack evolution due to the magnitude of the elastic energy density σ∶e . A material property γ * (x) > 0 prescribes the critical value of the energy release rate in the Griffith criterion (see Remark 7 in Sect. 2.4). As a crack once generated can no longer be healed, we apply ( ) + to the right hand side, where (a) + = max(a, 0) . This guarantees the irreversible growth condition for the crack: z t ≥ 0 . We remark that the second equation is a fully nonlinear parabolic equation, and it is called an irreversible diffusion equation or a unidirectional gradient flow [3]. Please see remark 3 in Sect. 2.4 for more detail.

Energy estimate
Under some suitable regularity assumptions, formally we have the following energy decay property for the solution (u(x, t), z(x, t)) of (2.14). We define where the first, second and third terms in r.h.s represent the rates of energy injection through the body force, surface traction and boundary displacement, respectively. We remark that Ḟ = 0 if f, q, g are independent of t. Then, similarly to (2.13), we have where we used the equality aa + = (a + ) 2 with a + = α z t . This suggests that the gradient flow structure of our phase (2.14) in Ω.
field model (2.14) is valid even with the irreversible growth condition.

Some remarks
(1) Numerical implementation Using our mathematical model (2.14), we do not experience any difficulty numerically simulating crack growth. This model equation enables us to simulate the crack growth phenomena by using any standard numerical method such as FEM on the fixed domain Ω . The adaptive mesh FEM technique can also work very effectively, as the cracked region is highly localized. We show some numerical examples of crack growth using our phase field model in the next section. Similar mathematical models based on the Francfort-Marigo type energy and their FEM simulations can also be found in [15,16,19,26], etc.
(2) Scalar mode III crack growth model In [34], the authors proposed a similar phase field model to (2.14) in two dimensions with scalar anti-plane displace- where μ > 0 denotes the shear modulus, one of the Lamé constants. This is a mode III crack growth model, and several numerical examples of mode III crack phenomena, such as crack branching and merging, are shown in [34].
(3) Irreversible growth condition and irreversible diffusion equation The second equation of (2.14) belongs to a class of irreversible diffusion equations (on unidirectional gradient flow). A typical irreversible diffusion equation is u t = (δu + f (x)) + , for example. This is fully nonlinear, but the theory of viscosity solution of parabolic type [46] can be applied. The problem admits a unique viscosity in Ω, solution, which is a kind of weak solution, and satisfies the comparison principle. One of the authors also recently established global existence of a unique strong solution for the irreversible diffusion equation u t = (δu + f (x, t)) + and its gradient flow structure in [3].
(4) Solvability of u Let (u, z) be a sufficiently smooth solution to (2.14). When z 0 (x) ∈ [0, 1) , the solution z also satisfies z(x, t) ∈ [0, 1) for t > 0 . This property follows from the comparison principle. Actually, when z(x, t) = 1 occurs at some x ∈ Ω , the unique solvability of u(t) is not clear, as the coercivity of the linear elasticity equation is not guaranteed. To avoid this mathematical difficulty, in many studies, such as [18], the term (1 − z) 2 in (2.14) is replaced by (1 − z) 2 + η 2 with a small η > 0.
In this paper, however, we omit the term η 2 for simplicity, as we do not have any difficulty, at least in our numerical simulations.
(5) Profile of z 0 (x) for initial cracks The crack width in our phase field model is O(ϵ) . One of the possible choices for the initial condition z 0 (x) is as follows. For simplicity, we consider a one-dimensional ini- Fig. 3a. In this case, we set For the other configuration of the initial crack in 2D or 3D, we define z 0 similarly to a suitable shift and superposition of the above ζ 0 (x) (see Figs. 1, 2). In our numerical simulations, the initial condition of z is always chosen similarly. Therefore, we omit the precise description of z 0 (x) for simplicity in the following numerical examples.
(6) Damaged elastic modulus for isotropic material For a homogeneous isotropic material, the elasticity tensor and the stress tensor are given by where λ and μ are Lamé's constants. The positivity condition (2.3) is satisfied for c * = 2μ + d min(λ, 0) if μ > 0 and λ > −(2∕d)μ . As is well known, when d = 3 , Young's modulus E Y > 0 [Pa] and Poisson's ratio ν P ∈ (−1, 1∕2) are represented by The definition of the damaged stress tensor (2.8) means that the elasticity tensor c ijkl is modified as or, in the isotropic case, by However, (2.16) is written in the following form in terms of E Y and ν P : This means that the damage variable z affects the Young modulus but does not change the Poisson ratio. In other words, the following relations have been implicitly supposed in (2.8): -The damaged variable is given by -On the other hand, the Poisson ratio does not change while the crack growth, i.e., ν P = ν P .
(7) Choice of critical energy release rate γ * The energy release rate is the strain energy that is released during crack growth per unit surface area of the newly created crack. In classic Griffith fracture theory, the crack can grow if the energy release rate becomes greater than or equal to a critical energy release rate usually denoted by G c = 2γ , which is denoted by γ * [Pa m = J m −2 , for d = 3 ] in this paper. The critical energy release rate γ * is a material property, and represents the toughness of the material against the fracture. For the mode I crack in 2D, the relation is the critical stress intensity factor of mode I and where E � Y = E Y for plane stress and [47]. On the other hand, to apply our model to real-world applications, we have to take into account the difference between the ideal simplified situation and the real one. While our phase field model was derived from the theory of variational fracture for a brittle material, in a real-world situation, there are miscellaneous factors that weaken or strengthen the material, such as impurities, dislocations, defects and voids, microcracks, fatigue, vibration, shock, and impact, HE, plasticity, and so on. In those cases, we sometimes have to choose an effective critical energy release rate γ eff * . An example of γ eff * is shown in Sect. 4.

(8) Inhibition of crack nucleation
According to conventional (nonregularized) linear elasticity theory, the solution of the elasticity equations (2.1) for a sharply cracked domain, as in Fig. 1, has a stress singularity at the crack tip; that is, the amplitude of the stress |σ(x)| becomes infinity, as the position x approaches the crack tip. This means that fracture is most likely to occur at the crack tip. On the other hand, in our phase field model, the stress singularity at the crack tip is regularized on a length scale O(ϵ) , and the distinct priority for fracture at the crack tip is lost. Depending on the boundary conditions (domain shape and distribution of boundary loading), a position x 0 ( ∈ Ω⧵Σ ) far from the crack tip may have a stress amplitude |σ 0 | comparative to or larger than that around the crack tip, which leads to the initiation of a new crack from x 0 . To avoid initiation of a crack and its interference with the pre-existing crack, we have to carefully choose the small parameter ϵ > 0 to satisfy To interpret the above inequality, we suppose λ ≥ 0 and consider a situation whereby an undamaged material (satisfying z = 0 everywhere for t < 0 ) is exposed to a (spatially and temporally) constant stress σ 0 for t ≥ 0 . Let e 0 be the strain tensors at t = 0 . The time evolution of z(t) obeys the following ODE: The solution z(t) to (2.18) is given by The condition for the damage accumulation by |σ 0 | to be negligible is z ∞ ≪ 1 , that is, , where the factor 1 + ν P is omitted for simplicity. Conversely, we can deliberately cause crack initiation by choosing ϵ as where σ 0 (x 0 , t) is the stress of uncracked but strongly deformed positions. This is another interesting feature of the phase field model, but we do not discuss whether crack initiation in our model can adequately describe that in actual materials.
For a small parameter α > 0 , we further apply (2.20) as where the new regularization term is O(τ) as τ → 0 , since Although this regularization term vanishes as τ → 0 , it has a localization effect; i.e., z k in (2.20) is a global minimizer of E * (kτ, ⋅) , while z k in (2.21) is an approximation of a local minimizer of E * (kτ, ⋅) near z k−1 (see Remark 10 below). The constrained minimizing movement [48] defined by (2.21) is formally equivalent to the Euler-Lagrange equation: For more detail on the derivation of (2.22), we refer the reader to [3], in which a similar constrained minimizing movement problem is considered. The second equation of (2.14) is naturally derived as a limit of (2.22) as τ → 0. The coefficient α is mathematically introduced as a time constant that governs the time scale of the gradient flow dynamics of z. On the other hand, physically, α is the coefficient of resistance against the time change of z. When a crack propagates at a nonzero velocity V, the original critical energy release rate γ * is slightly modified as γ mod * (V ) , and contains an excess term related to α and V. A simple model calculation shows that γ mod * (V ) = γ * + (const.)αV . This relation is implied by dimensional analysis of . Details of the mechanical meaning of α will be reported elsewhere in a separate study.
(10) Advantage of the gradient flow model Although variational fracture theory has many strong points compared with other models, some difficulties were also pointed out in [1]. For example, (i) In the presence of a body force or boundary traction (f or q), the global minimizer might not exist or be unrealistic in general if the crack path is not prescribed. (ii) Even when the crack path is prescribed, the use of a global minimizer is not realistic, and a "local minimizer" often seems to be preferable. (iii) A numerical search for the global minimizer is not easy even for the regularized energy. In [18,19], an effective algorithm called the backtracking method was proposed. However, this algorithm is based on an assumption of a monotonically increasing load (MIL), which is a special type of loading.
On the other hand, our approach enables us to track local minimizers and to avoid the above difficulties. This is another major advantage of our phase field model approach.

Phase field model under unilateral contact conditions
In the case where the boundary condition is compressive, the linear crack problem (2.1) has an unrealistic overlapped deformation near the crack. To avoid this deformation, we have to impose the unilateral contact condition: where ν denotes the unit normal vector on Σ in the direction from Σ − to Σ + . In the weak form, we add the unilateral constraint (2.23) to the minimizing problem (2.5). Even in the phase field approach, we have the same problem. To avoid it, Amor et al. [21] proposed the following modification to the damaged elastic energy for isotropic elasticity.
We define the deviatoric part of the strain tensor by Then, we have where ( div u) + and ( div u) − = (− div u) + are the positive and negative parts of div u , respectively, and we set β ∶= λ + 2μ∕d . The three terms on the right hand side of (2.24) represent the components of stress related to expansion, compression, and shear. Instead of the damaged stress tensor σ[u] and the modified elastic energy E 1 defined in (2.8) and (2.9), we introduce the following modified damaged stress tensor and elastic energy: Then, the first and second equations in (2.14) are replaced by We remark that Chambolle et al. [33] proved Γ-convergence of the above type of energy to the energy of a sharp crack with the unilateral contact condition (2.23).
In this paper, we do not give any numerical example with the unilateral contact condition. However, the above discussion reveals that attributes of the phase field model with the unilateral contact condition proposed by [21,33] can be easily included in our gradient flow type phase field model.

Dynamic fracture model
We can also extend the phase field model (2.14) to model a dynamic fracture. For this purpose, we simply replace the first equation of (2.14) by with additional initial conditions for u: Research Article SN Applied Sciences (2021) 3:781 | https://doi.org/10.1007/s42452-021-04593-6 Even with this modification, the energy dissipation property (2.15) still holds in the following sense: where ρ > 0 is the density of the material and where the newly added term dx represents the kinetic energy. We refer the reader to Bourdin et al. [23] for numerical experiments of similar dynamic fracture models and to Abe-Kimura [49] for a study of discrete analogs of the above energy dissipation property of a vibration-fracture model based on a one-dimensional spring-mass system.

Summary of our phase field model
In this section, starting with a brief review on the variational fracture theory [1,19], we derive an irreversible phase field model as a unidirectional gradient flow [3], which was originally proposed in [34] for the mode III crack case. It should be stressed that, unlike the other previous phase field models, our model naturally and simultaneously realizes the irreversible crack growth condition z t ≥ 0 (nonhealing property of the crack) and the energy dissipation property (2.15). Following several remarks in Sect. 2.4 from the viewpoints of mathematical analysis and numerical implementation, we further extend our models to the unilateral contact condition (Sect. 2.5) and to dynamics fracture (Sect. 2.6). Since our derived model (2.14) becomes an initial and boundary value problem of an elliptic-parabolic PDE system, this extension allows us a relatively easy numerical implementation (Sect. 3) and flexible extension to more complicated industrial applications (Sects. 4-6).

Time discretization
We show some numerical results for the phase field model (2.14) of crack growth derived in the previous section. We use the software FreeFem++ [43] for our computation.
In the following examples, for simplicity, we suppose that f ≡ 0 , q ≡ 0 and Γ N = Γ N 0 (i.e. Γ N 1 = � ). We fix τ > 0 as a constant time step and denote the approximated solution of u and z at t = kτ (k = 0, 1, 2, …) by u k (x) and z k (x) , respectively. We also define g k ∶= g(⋅, kτ) on Γ D and consider the following time discretization for (2.14): dx +Ḟ(t, u(t), z(t)) ≤Ḟ(t, u(t), z(t)), At each time step, the displacement u k and the intermediate valuable z k are computed by solving linear elliptic equations for u k (3.1) and for z k (3.2), respectively. This scheme, proposed in [34], allows us to compute z k explicitly from z k and z k−1 by (3.3) without solving any nonlinear equation. We remark that 0 ≤ z k ≤ 1 is guaranteed when z 0 ∈ [0, 1] at t = 0 from the comparison principle for (3.2).
The weak forms of (3.1) and These weak forms are easily computed by the standard FEM. However, as the profiles of u and z have small O(ϵ)-scale spatial patterns, we have to choose the mesh size h to be as small as h < ϵ around the crack region.
In   We also remark that the deformed domain Ω deformed = {x + βu(x); x ∈ Ω} is drawn with a suitable magnification factor β > 0 in the figures. The other nondimensionalized values of parameters for each numerical example are listed in Table 1.

Numerical examples in 2D
Here, we show some numerical results of 2D crack growth. We set Ω ∶= (−1, 1) 2 ∈ ℝ 2 and The initial crack is set as shown in Fig. 3a for the examples of Figs. 4 and 6, and (b) for Fig. 5. The simulation uses the adaptive mesh P2 FEM on FreeFem++, where the mesh is adapted at each time step to the phase field variable z k . The minimum mesh size (hmin) is set to be more than 2 × 10 −3 , and the maximum number of vertices (nbvx) of the triangular mesh is set be less than 40000. Figure 4 shows a straight crack growth when the boundary displacement is given as g = ±t (0, 10) on Γ ± D . Figure 5 shows the final shape of the crack under three different boundary conditions. The boundary displacement is given as g = ±10t (sin θ, cos θ) on Γ ± D for three different angles θ = 0 , π∕3 , and 4π∕9 . The numerical results show that the kink shape appears when θ > 0 , while the crack propagation is straight when θ = 0. Figure 6 shows a straight crack growth when the boundary displacement is given as In the upper figures in Fig. 6, the crack first grows until approximatelyt t ≈ 0.05 but then stops and keeps its length owing to the irreversible growth condition α z∕ t = (⋯) + . However, the crack line disappears in case we adopt a model without it such as α z∕ t = (⋯) (Fig. 6  lower). These results reveal that the irreversible growth condition is necessary when the applied load does not always monotonically increase.

Numerical examples in 3D
The phase field model based on the variational fracture theory exhibits its advantage, especially in a 3D crack growth simulation. In this subsection, we give several numerical examples of crack growth in 3D brittle materials using our phase field model (2.14).
In the first two examples, we set Ω , Γ D , and an initial crack as shown in Fig. 7a, b: The boundary displacement is given as g = ±t (0, 2, 0) on Γ ± D . This configuration is similar to the geometry of the compact tension specimen experiment in the ASTM E647-00 standard test [50], chapter 7 of [47].
In Figs. 8 and 9, numerical results of mode I (opening mode) crack growth from a horizontal initial crack (Fig. 7a) are shown. We set the finite element tetrahedral mesh on Ω of (3.5) with division of the box boundary as 20 × 20 × 5 , 30 × 30 × 7 , 40 × 40 × 10 , and 50 × 50 × 12 . Table 2 lists the mesh settings, number of vertices N v , and number of tetrahedrons N t of each mesh. Figure 9 shows that the temporal evolution of the crack surface area is The figure suggests that the crack speed is affected by the mesh size, but the crack surface area approaches the limit, roughly speaking, when N v ≥ 4868 . Figure 8 shows the temporal evolution of the deformed shapes when N v = 4868 (upper) and N v = 3578 (lower), where the crack speed of N v = 3578 is a slightly slower than in the case of N v = 4868.
In Fig. 10, we use the mesh with N v = 4868 and set the initial crack to be slanted, as shown in Fig. 7b. The initially slanted crack rotates and becomes horizontal as the crack propagates, thus exhibiting a mixed mode of crack growth (mode I + mode II + mode III). Next, we set and assume a horizontal initial crack (Fig. 7c) and the boundary condition g = ∓t (0, 0, 6) on Γ ± D . The finite element tetrahedral mesh is generated from the boundary division of 40 × 40 × 10 . The results are shown in Fig. 11, where we can observe a crack growth of mode III (shear mode).
The last example shown in Fig. 12 gives the dynamic behavior of crack propagation. We set and assume a slanted initial crack inside of Ω and the boundary condition g = ±t (0, 0.5, 0) on Γ ± D . The finite element tetrahedral mesh is generated from the boundary division 50 × 50 × 25 . Around time t = 1.2 , the initial hidden crack appears on the box surface, and simultaneously, the right half of the bottom boundary Γ − D is broken. The crack dynamically rotates ( t = 1.6 ), ultimately leading to rapture ( t = 2).

Discussion
In this section, several numerical examples of the phase field model (2.14) for 2D and 3D crack propagation were shown. Several types of complex patterns for crack propagation in 3D were reproduced under various boundary conditions; in particular, a dynamically turning crack pattern in mixed-mode (modes I + III) was obtained in Fig. 12.
We confirmed the accuracy of the numerical computation by numerically checking that the crack area converges to some value when the mesh size tends to zero in Fig. 9. On this point, let us give a brief comparison with previous models for crack propagation from the viewpoints of modeling and numerics. As explained in the introduction, the phase field model has great advantages (items 1-4) over other discrete models such as XFEM and DEM, especially in 3D problems. This advantage has been shown through our numerical examples in this section.
On the other hand, even among other crack phase field models which have been rapidly developed in the past decade, our model has the following advantages. The Bourdin model [16] is also based on the Ambrosio-Tortorelli approximation of the total energy, which is essentially the same as our E = E 1 + E 2 , and it consists of energy minimization at each time step. Compared with our model, it is closer to the original idea of the variational fracture theory [1] since our model substitutes a gradient flow for the minimizing sequence. As a drawback, Bourdin's model requires a special numerical implementation, but our model does not since it is an elliptic-parabolic system of PDEs.
The model proposed by Karma et al. in [26] seems to be motivated by a bistable-type phase field model for phase transition phenomena, and it also allows us relatively easy numerical implementation. However, it slightly lacks consistency with fracture theory compared with Bourdin's model or ours. The irreversibility condition is also not explicitly included in their phase field model, unlike ours.
As seen in this section, our phase field model can be computed by a standard finite element code without any other special modifications. As the crack is very thin and localized in the domain Ω , the adaptive mesh FEM is often effective in these kinds of simulations with the phase field model, especially in the 3D case. For example, we can use adaptmesh() of FreeFem++ in 2D case. For a   Crack surface area Fig. 9 Relationship between the number of vertices and the crack surface area similar purpose, we used the ALBERTA toolbox [51] in [34,35]. For numerical studies on the efficiency of the adaptive mesh FEM for pattern formation with a reaction diffusion system in 2D and 3D, we refer the readers to the numerical studies in [52,53].

Crack propagation model under thermal stress
If the temperature of the material is not constant in space or time, the thermal expansion effect should be taken into account. We denote a reference temperature by Θ 0 ∈ ℝ and the coefficient of linear thermal expansion of the material by a L > 0 . For example, if the temperature field is given by Θ = Θ(x, t) , the original length l 0 changes to l = l 0 (1 + a L (Θ − Θ 0 )) at (x, t).
We define the strain and stress tensors with the thermal expansion effect as follows: where e[u] is defined by (2.2) and I = (δ ij ) ∈ ℝ d×d sym . For a fixed time, when a temperature field Θ(x) , a body force f(x), a surface traction q(x) on Γ N , and a boundary displacement g(x) on Γ D are given, we consider the following thermal stress problem: there exists a unique weak solution for (4.1), and it is characterized as a unique minimizer of the following energy: where V(g) is defined by (2.10).
Modifying the crack propagation model (2.14), we naturally extend it to the case of crack propagation under thermal stress. For simplicity, we assume that a temperature field Θ(t) ∈ L 2 (Ω) is given. The displacement u(t) ∈ H 1 (Ω, ℝ d ) and the damage variable z(t) ∈ Z solve the following system: We remark that Bourdin et al. studied the complex morphology of cracks under thermal stress using the extension of the Bourdin-type phase field model in [24].

Nondimensionalization
In (4.2), we assume that the space dimension is d = 3 , and we use [m] for the length unit and [Pa = J m −3 ] for the stress. Then, the material parameters of a typical solder are as shown in Table 3.
We consider the following transform of variables by scaling. For simplicity, we suppose a homogeneous situation: f = 0 , g = 0 , and q = 0 . Let x, u, Θ, C, α, ϵ, and γ * (4.1) in Ω. be the original variables and parameters of (4.2) as in the middle column of Table 3. We denote the scaled nondimensionalized variables by symbols with tildes: where c 1 , … , c 4 are positive scaling constants. Using the nondimensionalized variables, (4.2) becomes where In the following numerical example, we use c 1 = 10 −3 , c 2 = 1 × 10 −6 , c 3 = 1 × 10 9 , and c 4 = 1 . The resulting nondimensionalized values appear in the right column of Table 3.

Numerical simulation of solder cracking
One of the severe problems caused by cracking under thermal stress is the fatigue of a solder joint. This fatigue occurs mainly due to cyclic loading stemming from mechanical vibrations or periodic temperature changes.
In this section, we focus on the solder cracking problem under a periodic change in thermal stress corresponding to a standard thermal cycling test. In our numerical examples, we consider two types of solder joint shapes: (a) fillet shape and (b) volcano shape (see Fig. 13). Both of them are common soldering techniques for making electronic packages. In these setups, we suppose that the electronic components, such as the electric pin, integrated circuit (IC) component, and printed circuit board (PCB), are rigid bodies, and the solder material is a homogeneous and isotropic elastic body. We set Ω as the solder region and define the Dirichlet and Neumann boundaries as shown in Fig. 13. For the several physical properties, supposing a typical tin-lead (Sn-Pb) solder, we choose the values as in the middle column of Table 3.
We remark that the effective critical energy release rate γ eff * = 4.85 [Pa m] is adopted instead of the theoretical value γ * = 86.31 [Pa m] in the numerical examples, as explained at remark 7 in Sect. 2.4. Actually, in the case of solder cracking, thermal effects on the material properties [56] and inhomogeneity of the solder joint due to the poor soldering process [57] are the main reasons for this difference between γ * and γ eff * . As a thermal cycling test, we replace the temperature between −55 ∼ 125 [ • C] with the reference temperature Θ 0 = 35[ • C] . Here, we consider 100 cycles over 10 hours (i.e. 0 ≤ t ≤ 10 ) with a time step of τ = 10 −3 , which means the external temperature is Θ(t) = 90 sin(20πt) In addition, the important assumptions made in our simulation are as follows: first, the solder material does not contain an initial crack z 0 = 0 ; second, we omit the body force f = 0 ; and last, the given displacement and surface force are neglected or mathematically denoted as g = 0 and q = 0.
The crack phase field model under thermal stress (4.2) is discretized in time in a similar way to (3.1)-(3.3), and the space discretization is implemented by the adaptive mesh FEM with P2 element on FreeFem++ [43], where the mesh is adapted to the profile of z at each time step. The minimum mesh size (hmin) is 2 × 10 −3 , and the maximum number of vertices (nbvx) of the triangular mesh is 40,000.
The temporal evolution of the numerical solution for the fillet shape is shown in Fig. 14, where one crack appears near the IC component and the other crack appears near the PCB ( t = 5 ). Those cracks grow simultaneously; however, the crack near the PCB becomes more dominant at t = 5 and grows continuously until t = 10 toward an angle of approximately 45 • .
In solders with a fillet shape, the numerical result of solder cracking is qualitatively very similar to that obtained by experimental investigations [54,58]. Figure 16 shows an example of a real cracking experiment conducted with a fillet shape solder, where the crack path starts from the zone between the IC component and the PCB, and the crack direction changes at an angle of 45 • . The crack path is not straight in Fig. 16 owing to voids and inhomogeneity of the solder as a result of a poor soldering process [57].
In the case of solder with a volcano shape (Fig. 15), the initial cracks appear on both sides of the solder border at t = 5 and grow simultaneously ( t = 10 ) at an angle of approximately 45 • degrees. This result is similar to the experimental results obtained by Hillman et al. (Fig. 3 of [55]).

Discussion
As an application of our phase field model of crack growth (2.14), we consider a solder cracking problem under cyclic thermal stress. This is a typical example in which the irreversibility condition is essentially required.
The simulated crack propagation exhibits good qualitative agreement with an experimental observation by [58] (Fig. 16), and it convinces us of the usefulness of our approach to the solder cracking problem. Inheriting the merit of our phase field model (2.14), our approach construct a mathematically clear and easily implementable Fig. 13 Computational domain of the solder problem: a Fillet shape [54] and b volcano shape [55] mathematical model for the thermal stress problem in solder. This is a good advantage over previous approaches using the XFEM-based model.
On the other hand, the value of the effective critical energy release rate γ eff * that we introduced here should be more precisely investigated through quantitative comparisons between numerical simulation and experimental data. This is one of the important open problems in this approach.

Maxwell-type viscoelasticity model
Polymeric materials are indispensable in modern society. Various polymers with a wide range of physical (mechanical, thermal, and optical) properties have been developed to meet industrial requirements. Glassy and semicrystalline polymers used as structural materials possess all the following favorable properties: lightness, water resistance, a certain level of hardness, deformability and toughness. On the other hand, rubbery and sticky polymers are used for damping, shock absorbing and detachable gluing. It is well recognized that viscoelasticity substantially influences the fracture behavior of polymers (especially soft polymers). However, a comprehensive understanding of the effect of viscoelasticity on fractures has not been established. This is because it is not an easy task to suitably incorporate viscoelasticity into linear fracture mechanics. In addition, the viscoelastic effect in fractures is often coupled with other complicated effects such as geometrical nonlinearities (crack blunting) and the localized temperature change around the crack tip (viscoelastic dissipation can increase the local temperature, and an increase in temperature modifies the viscoelastic properties).
In this section, we describe how to incorporate the bulk viscoelastic effect in the framework of our phase field fracture model.
In phenomenological descriptions of viscoelasticity, it is convenient to refer to symbolic mechanical elements consisting of springs and dashpots, instead of considering a concrete microscopic structure of a particular material. Figure 17 shows one of the simplest cases; the Maxwell element. We denote the stress and strain tensors with respect to the elastic spring and the viscous dashpot by (σ 1 , e 1 ) and (σ 2 , e 2 ) , respectively. Figure 17 suggests that the elastic and viscous elements sustain a common stress and that the total strain is the sum of e 1 and e 2 : We also denote the elastic modulus of the spring and the viscosity of the dashpot by C and η > 0 , respectively. The constitutive equations are given as We denote e 2 by ê =ê(x, t) ∈ ℝ d×d sym and call ê a viscous component of the strain. This componet represents an internal strain related to the viscous dashpot. Summarizing (5.1) (5.1) σ = σ 1 = σ 2 , e = e 1 + e 2 .   where f is a given body force. We emphasize that (5.3) has a gradient flow structure with respect to the following total energy: We refer the readers to [59] for more detail on the formulation (5.3) and its variational structure. Owing to the existence of such system energy, we can extend the phase field crack growth model to viscoelastic materials while maintaining their variational structure. As explained in many texts, in the scalar setting, using the stress σ and the strain e, the constitutive equation (i.e., the relationship among stress, strain, and their time derivatives) of a Maxwell element is given as where μ and η are the elastic constant and viscosity of the spring and dashpot elements, respectively. The above model (5.3) is nothing but another expression of (5.6) in terms of the viscous component of strain ê.

Crack growth model
To obtain the phase field crack growth model for the linear viscoelastic material, we replace the elasticity tensor C with the damaged tensor, and derive the gradient flow of E . Similar to the purely linear elastic models, we obtain In the next part, we provide numerical results of a mode I fracture of a Maxwell-type material based on (5.8) under a condition of d = 2 and f (x, t) ≡ 0.

Numerical results
Choosing an adequate time step τ , we discretize (5.8) with respect to time as follows: in Ω.
We solve (5.9)-(5.12) by the adaptive mesh P2 FEM implemented on FreeFem++. Figure 18 compares the deformation and fracture behaviors for different η . For η = 4 , the initial crack gradually opens but never extends; instead, the viscous component of strain ê increases in the uncracked region ahead of the crack tip. For η = 6 , the initial crack extends in a quasibrittle manner. For η = 5 , the system seems to be near the marginal state of fracture; very slow crack extension and development of a large-ê region occur at the same time. Figure 19 represents the profiles of z and |ê| 2 on the line along the initial crack ( x 2 = 0 ). The "cliff" of the z-profile corresponds to the crack tip. For smaller η , the z-profile is almost stationary; moreover, |ê| 2 in the region ahead of the crack tip gradually increases with time. For larger η , the crack tip rapidly propagates toward the right end of the system at a time interval of 0.6 < t < 0.8 , leaving the |ê| 2 -profile ahead of the initial crack unchanged. 6 Hydrogen-assisted cracking in metal

Effect of HE
Actual fracture processes are frequently influenced by environmental factors, especially the permeation of chemically active species into the so-called process zone around crack tips and by the resultant change in the mechanism of cutting atomic bonds at these locations. Stress corrosion cracking and HE are typical examples of deterioration of mechanical performance by anodic and cathodic reactions, respectively. The HE is considered to be caused by accumulated hydrogen atoms around the crack tip. The detailed mechanism is, however, still a controversial issue because the entire HE process involves many mutually affecting processes, such as adsorption and penetration of hydrogen atoms and their interaction with dislocations. Several mathematical models have been proposed to describe HE, some of which assumes specific mechanisms for gathering and trapping hydrogen atoms around the crack tip, e.g., the decrease in the chemical potential of hydrogen by the negative pressure (i.e., reduction in the hydrostatic component of the stress tensor) produced by the strong elongation around the crack tip.
In this paper, we propose a phase field model for HE by modifying our previous model [61]. In this model, instead of assuming a gathering/trapping mechanism, hydrogen atoms are assumed to be able to accumulate around the crack tip by a sort of "jamming" effect. They are transported by fast diffusion along the crack surface and slowly diffuse to the crystal lattice of the metal via the near-crack tip region because the diffusion constant of hydrogen atoms in the vacancy is much larger than in the metal crystal. Hence, the transported hydrogen is concentrated around the "entrance", that is, the nearcrack tip region of the metal.
To construct a phase field fracture model of HE, we introduce a normalized density of hydrogen atoms w(x, t) ∈ [0, 1] and make two assumptions: (i) Here, γ * is a decreasing function of w and (ii) the diffusion coefficient of hydrogen is an increasing function z. According to (i), we assume a linear dependence of γ * on w where γ * 0 and γ * 1 are constants with γ * 0 > γ * 1 > 0 . With the aid of the Arrhenius law, ii) is formulated as where A is a constant (with the units of the diffusion coefficient) and D j = D w (j) > 0 for j = 0, 1 . Note that the factor E a (z) = E 0 a (1 − z) + E 1 a z over k B T in (6.3) denotes an activation energy, where E 0 a > E 1 a > 0 . As the activation energy E a (z) is decreasing in z, D w (z) becomes an increasing function of z.
We employ the following set of equations as a model of HE (for simplicity, we set f ≡ 0 and q ≡ 0).

Numerical results
With an adequate time step τ and denoting the approximated solutions of u, z, and w at t = kτ ( k = 0, 1, 2, … ) by u k (x) , z k (x) , and w k (x) , respectively, we discretize (6.4) as follows: As shown in Fig. 20, we test two setups: (a) hydrogen supply from the cracked left side and (b) hydrogen supply from a hole. In setup (a), the system is a 2D square domain Ω ∶= (−1, 1) 2 with an initial crack from the left boundary Γ 1 . We set w 1 = 1 on Γ 1 , which means that a large hydrogen reservoir makes contact with the left boundary to supply hydrogen to the inside of Ω . Note that D w (z) steeply increases as z approaches 1, which means that hydrogen can diffuse very rapidly along the initial crack where z ≈ 1.
In setup (b), the domain is given as Ω ∶= (−1, 1) 2 ⧵H with the same initial crack as (a), where H is a hole that is far from the initial crack defined by We set Γ 1 ∶= H and w 1 = 1 on Γ 1 , which means that the hydrogen is supplied from the hole H and permeates into the inside of Ω.
In both cases, we set the Poisson ratio ν P = 0.29 , the nondimensionalized Young modulus E Y = 50 , γ * 0 = 0.5 , γ * 1 = 0.01 , ϵ = 0.05 , α = 0.01 , and τ = 0.0001 . We perform numerical simulations with different D 0 , D 1 ; the employed values are described below. We also impose the Dirichlet boundary condition for u as g(x, t) ∶= ±t(0, 1) on We solve (6.5)-(6.7) by the adaptive mesh P2 FEM implemented on FreeFem++. Figures 21 and 22 show numerical results in setup (a). Under the supply of hydrogen, rapid crack growth occurs in an earlier period between t = 0.1 and 0.2 (Fig. 22). This rupture occurs at an earlier stage than in the case without (6.5) in Ω, HE (Fig. 21). This is due to the reduction in γ * by increasing w, that is, hydrogen-assisting cracking. To see how the hydrogen-assisting cracking occurs in our model, we compare, in Fig. 23, the profiles of z and w along the line of x 2 = 0 (the extended line of the initial crack) at different times for (D 0 , D 1 ) = (1, 100) (upper) and (D 0 , D 1 ) = (0.1, 100) (lower). Note that the cliff of the z profile corresponds with the crack tip. For a relatively large D 0 of D 0 = 1 , hydrogen permeates into the region ahead of the initial crack at an earlier time of t = 0.1 to increase w at the position of the tip of the initial crack; see the thick arrows in the upper-left and upper-middle graphs in Fig. 23. The realization of such profiles of w triggers the subsequent rapid crack extension. The crack tip has already arrived at the right end of the system at t = 0.2 . On the other hand, for the smaller D 0 of D 0 = 0.1 , the permeation of hydrogen hardly occurs by the time t = 0.1 . After the stretching deformation reaches a higher level, the crack tip suddenly propagates to break the material without the assistance of hydrogen. Figure 24 shows a numerical result for setup (b). Although at earlier times ( t ≈ 0.01 ), the hydrogen is axisymmetrically distributed around hole H, the distribution is elongated along the x 1 direction at approximately t = 0.02 . This corresponds to the initiation of new cracks from the left and right ends of the hole; i.e., z and D w (z) increase there. At t ≈ 0.03 , the left and right secondary cracks merge with the initial crack and the right boundary of the system, respectively, resulting in the complete system fracture.
For comparison, we investigate, by setting γ * = 0.5 (independent of w), the situation whereby the fracture behavior (the time evolution of z(x, t)) is decoupled from hydrogen distribution w(x, y) (Fig. 25). (w(x, t) is influenced by z(x, t) because of the z dependence of D w .) Although the final shape of the fracture surface (the upper right figure in Fig. 25) may seem similar to that in Fig. 24, its formation process is fairly different, as described below. In this decoupled condition, the crack tip propagates in a straight line in the earlier stage of the fracture and then turns toward the hole to disappear. When the system is further stretched, a new crack tip initiates at the right end of the hole and propagates toward the right boundary of the system. Note the difference in the kinking position of the fracture surfaces (crack tip trajectory) in Figs. 24 and 25.

Discussion
We have taken a different view on the mechanism of hydrogen condensation near the crack tip from that of previous studies [41,42,[62][63][64]. Yokobori et al. attributed the hydrogen condensation to a "thermodynamic gathering effect", that is, the decrease in the chemical potential of hydrogen due to the strong elongation and accompanying large negative pressure by stress concentration near the crack tip. This effect can be formulated by including the coupling terms of hydrogen concentration and displacement field (volume strain) in the total system energy. The recent phase field modeling of HE [41] is based on the essentially same idea (termed a "strong coupling model"). In contrust, when the system energy is included such that the coupling term can conserve the gradient flow structure, complicated evolution equations are produced with several mathematical drawbacks such as a risk of unsolvability.
Our model regards hydrogen condensation as a "kinetic jamming effect"; that is, hydrogen is transported rapidly along open cracks but accumulates near the crack tip because of the difference in the diffusion constants between the crack and the metal lattice. This physical picture can be formulated merely by making the diffusion coefficient of hydrogen an increasing function of the damage field and the fracture toughness a decreasing function of the hydrogen concentration ("weak coupling model").
Since the system energy function does not include the coupling term between the elastic field and hydrogen concentration, the resultant evolution equations are relatively simple and have mathematical clarity. This is the advantage of our method.
We do not deny the thermodynamic gathering mechanism. What we have shown here is that the simple kinetic effect also can explain the hydrogen accumulation and resultant embrittlement. To clarify the mechanism of hydrogen accumulation around the crack tip in the actual system, more detailed experimental research must be undertaken.

Discussion
Historically, computational fracture mechanics has been developed as a branch of the FEM, where physical modeling and construction of numerical schemes have been performed in complimentay manner. In this paper, we focus on modeling fracture dynamics by the gradient flowtype PDE. Numerical schemes are constructed via discretization of the PDE. The merits of this approach are that we can apply physical or intuitive considerations to the PDE (such as dimensional analysis and order estimation) and that we may utilize the mathematical theory of PDEs to evaluate the suitability of the models. We also emphasize that the system energies for the gradient flow have clear physical meanings. The exception is the modeling of hydrogen-assisted embrittlement; we assume the weak coupling between the mechanical fields (u and z) and the transport of the hydrogen molecules without deriving a "strict" gradient-flow system from total system energy consisting of mechanical terms (depending on z and u) and the coupling terms between hydrogen concentration and mechanical fields ("strong coupling model"). This is because the strong coupling model leads to complicated evolution equations. There is another strong reason for the choice of the weak coupling model. The hydrogen concentration obeys the conservation law, whereas the mechanical fields (u and z) are nonconservative. It is unreasonable to derive evolutionary equations of conservative and nonconservative field quantities from a single energy function.
Finally, we point out the possibility of a theoretical extension of the present model. The evolution equations in the phase field fracture model contain two model parameters ϵ and α , representing the thickness of the diffused crack surfaces. In actual numerical simulations, it is impossible to reduce ϵ down to a micron or smaller scale. Thus the shapes of cracks (that is, crack paths) are always fairly diffuse; our model is inadequate for simulating crack propagations that produce highly intricate crack paths such as fractal structures [65]. On the other hand, we may regard the length scale ϵ as corresponding to the size of the so-called process zone [47], the mesoscale region around the crack tips where dissipative nonelastic deformations take place. It is well known that the size of the process zone depends on the condition of fracture, such as the specimen size [66] and temperature, and affects the fracture behavior, e.g., the effective fracture toughness and crack propagation mode [67]. To test the possibility of such a theoretical extension, a suitable state equation or evolution equation for the variable ϵ(x, t) must be added.

Conclusion
We have considered several extensions of the phase field model (2.14), as introduced in Sect. 2 with an energy-consistent derivation. Owing to the adoption of the phase field approach, there are no geometric restrictions on the crack shape. The crack path search is implicitly included in the model, and complex crack patterns are naturally treated even in a 3D setting.
The characteristics of our model are summarized as follows.
1. The model is described as an initial-boundary value problem of an elliptic-parabolic PDE system with a unidirectional property that represents the irreversible growth condition for the crack. This allows us to simulate the model with standard finite element code without any special numerical technique. In particular, the adaptive mesh FEM works well in the 3D case owing to the spatially localized profile of the crack. 2. Based on the Griffith and Francfort-Marigo theories of variational fracture [1,36], the mathematical theories of the Ambrosio-Tortorelli regularization [2], and the unidirectional evolution equation of diffusion type [3], we derived our phase field model as a gradient flow of the total energy. Consequently, this model exhibits a natural energy dissipation property, as shown in Sect. 2.3. 3. As the model is an energy-consistent PDE model, we can naturally extend it to more complex and real situations, such as thermal stress cracking (Sect. 4), crack propagation in viscoelastic materials (Sect. 5), and hydrogen-assisted cracking in metals (Sect. 6).
Fracture mechanics were first explored in a pioneering study by Griffith [36], who succeeded in obtaining the fracture criterion for simple brittle materials such as silicon glass by describing fractures from an energy viewpoint. Francfort and Marigo [1] extended the idea proposed by Griffith and successfully established a basic concept of variational fracture. We utilized the ideas proposed by Griffith and Francfort-Marigo to derive a phase field model for practical application in Sect. 2. The flexible extendability of the phase field model is shown in the subsequent sections, and the highlighted examples prove the usefulness of the variational approach that was proposed in Griffith's work, even for fracture dynamics and crack propagation phenomena.

Data availability
The raw dataset can be shared upon reasonable request. Coda availability The numerical codes can be shared upon reasonable request for non-commercial purposes.

Conflict of interest
The authors declare that they have no conflict of interest.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.