Physics informed neural networks for phase field fracture modeling enhanced by length-scale decoupling degradation functions

The paper proposed a novel framework for efficient simulation of crack propagation in brittle materials. In the present work, the phase field represents the sharp crack surface with a diffuse fracture zone and captures the crack path implicitly. The partial differential equations of the phase field models are solved with physics informed neural networks (PINN) by minimizing the variational energy. We introduce to the PINN-based phase field model the degradation function that decouples the phase-field and physical length scales, whereby reducing the mesh density for resolving diffuse fracture zones. The numerical results demonstrate the accuracy and efficiency of the proposed algorithm.


Introduction
Predicting material and structural failure due to crack nucleation and extension is central to many engineering applications. Several numerical methods are developed to capture complex fracture phenomena, such as the finite element method [1], the boundary element method [2], the cohesion zone method [3,4], the extended finite element methods [5], the peridynamics Ha and Bobaru [6] and the phase field methods. In these methods, the phase field method has demonstrated advantages in description of complex fracture patterns. By introducing additional continuous field variables to track discrete discontinuities with diffuse fracture zones, the phase field method Bourdin et al. [7] unifies crack initiation, propagation, branching, and merging in structures.
In phase field method, the diffuse fracture region must be resolved with sufficient degrees of freedom to obtain an accurate solution. The length scale of the physical process zone is proportional to the ratio of the fracture energy to the square of the material strength. For most problems modeled by linear elastic fracture mechanics, the length scale of the physical process zone is very small compared to structures, but the commonly used formulas in fracture phase field approach does not distinguish between the phase field crack length scale the physical process zone length scale, which leads to prohibitive meshing requirement when extending the models with nanoscale phase field length to large structures. To alleviate the meshing burden, Wu et al. proposed a length-scale insensitive phase-field model [8], but the method is practically applicable to the scenario that the phase-field length scale and physical process zone length scale are in the same order of magnitude. Lo et al. [9] proposes a degradation function that separates the phase field length scale from physical length-scale, which enables one to simulate crack propagation in large scale structures.
As a branch of artificial intelligence, the deep learning based on Artificial Neural Networks has made tremendous progress with the explosive growth of data in the past decade [10][11][12] and applied in a wide range of areas. ANN were initially trained for Computer Visualization and Natural Language Processing tasks [13,14]. When it comes to the field of physical simulation, two major problems arise: firstly, it is cumbersome to obtain the data from complex engineering systems; secondly, the result predicted by the training data cannot ensure that the physical principles of the problem are satisfied. Hence, the effectiveness and reliability of ANN in characterizing physical phenomena are questionable. To overcome the difficulties, physics informed neural networks (PINN) [12,[15][16][17][18][19] were put forward, which comply with the distribution of the training data as much as possible and meanwhile obey the laws of physics which are commonly formulated by partial differential equations. Compared to purely data-driven neural network learning, PINN allows learning with fewer data samples and obtaining models with greater generalization capability. Goswami et al. [20], Goswami et al. [21] firstly applied PINN to solve the phase field model for simulating the growth and propagation of fractures in brittle materials. Their results offer improved accuracy and efficiency, demonstrating the advantages of PINN in simulating moving boundary problems.
In the present work, we introduce the degradation function proposed by Lo et al. [9] to the PINN-based phase field method [20,21], in order to analyse crack growth in large structures with a data driven and mechanism based hybrid approach. We use three numerical simulation examples to demonstrate that the physical neural network combined with the proposed degradation function is correct and advantageous. The remaining o the paper is organized as follows. Section 2 introduces fracture phase field method and shows how the phase field length scale and the physical process zone length are decoupled. Section 3 elaborates on how the fracture phase field model is discretized with the PINN. The numerical examples are given and analysed in Section 4, and Section 5 is the conclusion and outlook of future work.
2 Phase field model of cracks 2.1 Conventional phase field crack modeling The phase-field crack modeling of brittle fractures involves the integration of two fields: the elastic field u and the phase field d, and the crack propagation is determined according to the free energy minimization principle. Based on the energy decomposition proposed by Francfort and Marigo [22], the free energy of the fracture system is given as: where: G c is the critical energy release rate, ψ e is the elastic energy density, Γ is the fracture surface, and Ω is the problem domain. The right hand side of the equation is the sum of the elastic strain energy and the fracture surface energy, and the fracture phase field method in mechanics assumes that the crack should follow the direction of minimum free energy and be irreversible.
Since the boundary integral involved in the surface energy is not easy to handle, it is replaced by the crack density equation in the following: where: d is the order parameter; l is the diffuse crack width. Theoretically, the model draws on the elliptic regularization method of the Mumford-Shah generalization Lie et al. [23] in computer image segmentation (hence the model is also known as AT2).
Since the elastic strain energy cannot distinguish between positive and negative for stress and strain, Amor et al. [24] and Miehe et al. [25] proposed a tension-compression split of the elastic strain energy, which is defined as: where the elastic energy is determined by both the strain and the order parameter d. g(d) represents the degradation function, which is used to reduce the strength of material around the cracked region.
To model material failure and crack propagation, the material should be fully elastic when it is intact and disappears when it is completely cracked. This is achieved by multiplying the elastic energy by the degradation function in the phase field description of the damaged materials.
In addition, to prevent the release of elastic energy after unloading from causing crack healing, Miehe et al. [26] proposed a very concise and effective solution by introducing a historical strain function: Shape of the proposed degradation function in this work. As the parameter q increases the degradation function approaches the standard quadratic, and as q decreases towards 1 the peak stress predicted by the model increases.
Frontiers in Physics frontiersin.org Eq. 4 indicates that the crack driving force is taken as the historical maximum tensile elastic energy. Even after re-unloading, the tensile elastic energy maintains its maximum value and thus prevents crack healing.
Substituting the crack density in Eq. 2 and the elastic energy in Eq. 5 into the free energy of the fracture system in Eq. 1: Finally, we seek the solution to minimize the free energy E, which consists of the elastic strain energy ψ e and the fracture surface energy ψ c .

Degradation function
The degradation function needs to satisfy the following conditions: The most widely used degradation function in the literature is a quadratic function, With this degradation function, for the AT2 model (d 2 ) with uniform uniaxial tension, the peak stress is obtained as: Eq. 8 shows that the phase-field length, l 0 is not an independent parameter but dependent on the material strength, σ c , fracture energy, G c , and Young's modulus, E. This length scale, which represents the size of the physical process zone, l p , can be quite small for real materials. In the past, the typical approach has been to set l p equal to l 0 . However, this means that numerical solutions must have the ability to resolve the small l p , which can be a challenge for large engineering structures. This is because the mesh size around cracks must be a fraction of l 0 , at most l 0 /2, leading to computationally intensive simulations. To overcome the difficulties, Lo et al. [9] proposed a new degradation as: This degradation function meets all the requirements of Eq. 6. Its shape depends on the value of q. When q is 200, the shape of the proposed degradation function in Eq. 23 is close to that of the classical degradation function in Eq. 7. When the value of q is slightly larger than 1, the shape of the degradation function g(d) will dramatically alter following marginal change of the value of q. This can be seen in Figure 1. The peak stress for a bar under homogeneous uniaxial tension can be determined using the proposed degradation function, as follows: It's worth mentioning that we've introduced the notation σ c * to represent the peak stress calculated using the alternative degradation function, differentiating it from the conventional peak stress (σ c ). Additionally, the physical length scale l p in the case of the alternative degradation function is given as:

FIGURE 2
Schematic representation of the proposed physics informed neural network. X represents the input of the neural network, Y represents the output of the neural network, f(e) represents the elastic strain energy, and f(c) represents the fracture energy. For training, we have used the ADAM optimizer followed by L-BFGS.

FIGURE 3
Geometrical setup of a one-dimensional elastic bar with crack.
Frontiers in Physics frontiersin.org Frontiers in Physics frontiersin.org 04 This allows us to also express the peak stress with the proposed degradation function as follows: q can be used to maintain the same peak stress σ c * in numerical simulations with the increase of the phase-field length, which allows for a coarser mesh in the proximity of cracks. Supposing G c , E, and σ c * are given as constants, when q is large, the sizes of the physical process zone, l p , and the phase-field process zone, l 0 , are comparable. Moreover, the finite element mesh must be refined according to l 0 . If we would like to increase the l 0 by 4 times without changing the peak stress, σ c * or the l p . q needs to be 1.107 according to Eq. 10. If we aim to increase the length scale l 0 by 16 times while still maintaining the same peak stress, we find that q ≈ 1.0148. If we increase l 0 by 100 times, we find that q ≈ 1.00155. This demonstrates that the proposed degradation function allows for decoupling l 0 from l p . Note that these results assume homogeneous post-peak behavior and that localized deformations have been suppressed. Importantly, it's worth noting that the resulting peak stress remains the same for different phase-field length scales. Also note that the material is not damaged, i.e., d begins to increase from 0, until the strains exceed peak strains.

PINN-based fracture phase-field modeling
This section details the implementation of the phase field method based on PINN [20,21]. Firstly, the deep neural network will be described, and secondly, the synthesis between the PINN phase field method with the degradation function mentioned in the above section is explained.

Deep neural networks
Deep learning is a branch of machine learning based on deep neural networks. In comparison with shallow neural networks, deep learning has more hidden layers and is more capable of fitting non-linearities. Neural network training consists of two steps, i.e., forward propagation and backward propagation. The forward propagation is used to train the neural network parameters, while the backward propagation is used to update the neural network parameters and find the optimal solution. In this paper, a feed-forward deep neural network is used for the study. The hidden layer of the deep neural network contains its main parameter weights W and biases b, and the parameters are continuously optimized by optimization algorithms (e.g., adaptive moment estimation (ADAM), proposed Newton method (L-BFGS), etc.) to find the optimal solution. Supposing that the network consists of L hidden layers, with layer 0 denoting the input layer and layer (L + 1) denoting the output layer, the expression of the neural network can be written as: The calculation of the output Y in the feed-forward algorithm can be represented as follows: the activation function σ l−1 in layer l is used in conjunction with the number of neurons m l−1 in layer l − 1.
where x represents the input to the neural network in Eq. 14. Forward propagation is used to train the neural network parameters, (w, b), which represent all the W l and b l parameters that appear in Eq. 14. We evaluate a neural network prediction result by the loss function, and if the loss value reaches low enough, (which is generally not possible to be 0) we end the neural network training. The common loss functions for regression are MAE loss, MSE loss, and smooth L1 loss. Their expressions are given by the following equations.
Y i represents the target value and N(x i ; W, b) represents the predicted value of the neural network. Backward propagation is automatically performed based on the value of the loss function. Automated differentiation and weight update implemented in TensorFlow or PyTorch allow algorithm designers to perform their tasks without coding the back propagation from scratch. The reader is referred to LeCun et al. [27] for details of the adopted gradient calculation method and the layer-by-layer chain rule gradient derivation process. In the following section we focus on an algorithmic approach for solving problems using the variational energy-based physics informed neural network.

Variational energy based PINN
We can define the problem by a one-dimensional timeindependent differential equation: together with the Dirichlet boundary condition, which is defined as: The deep neural network can be utilized to compute the field variable, represented by u, by considering variables such as the Dirichlet boundary point, represented by x D , and the source term, represented by f(x). Furthermore, the first, second, and higher-order derivatives with respect to the independent variable x can be calculated using u x , u xx and u x. . .x respectively. The variational energy principle has the distinctive benefit of automatically satisfying homogeneous Neumann boundary conditions. At this point, the expression for the variational energy is: where y is the differentiable functional. Ω is the problem domain for which a solution is required, so the problem for which we require a solution can be expressed as follows: The steps for solving the differential equation using a neural network are detailed as follows: First, we construct a neural network, N(x; W, b) with initial parameters. Second, we revise the neural network according to the boundary conditions. In terms of addressing the boundary conditions in the context of neural network, weak form and the strong form. When dealing with boundary conditions in a weak form, a penalty term reflecting the boundary condition is added to the loss function. But this may cause interference among different loss terms, which slows or even cripples the convergence of the problem. So, we prefer the strong form here and we redesign the output of the neural network so as to comply with the Dirichlet boundary condition.To this end, we set: The functionũ D is selected such that it matches the value of u D at the Dirichlet boundary points. Meanwhile, B(x) is equal to zero at the same boundary. There is no requirement for a boundary component in the calculation of the loss function. Thirdly, the output of the neural network existing in variable energy expression is calculated using an automatic differential technique. Having computed the first and higher-order derivatives of u, we can now calculate the variational energy at each quadrature point. This is done by using u x from Eq. 22 and its derivatives obtained previously. The total variational energy of the system is then obtained by summing the energies at each point, as defined in Eq. 20. In this step we use the PDE (partial differential equation), i.e., the total variational energy, as the loss function of the neural network. The neural network constructed above is named physics informed neural network. The final step is to update the neural network parameters W and b by adjusting them so that the next neural network output value is closer to the target value. According to the above steps, Goswami et al. [20], Goswami et al. [21] proposed a PINN combined with phase field method to simulate the fracture problem. The displacement control is applied for the loading process. To train the PINN, a fixed displacement step, denoted by △u, is taken and the strain-history function is updated at each displacement increment. Before training the network, the weights are initialized randomly from a Gaussian distribution with the Xavier initialization technique [11]. In order to calculate the parameters of the neural network at the first Frontiers in Physics frontiersin.org displacement step i, the Gaussian quadrature points and their corresponding weights are generated using the above process, and in the subsequent step, the automatic differentiation method is employed to calculate the displacement increment, △u, and the eigenvalues of the strain, (λ 1 , . . ., λ d ), where d denotes the number of spatial dimensions. These eigenvalues are then utilized to determine ψ + e and ψ − e . Next we modify the degradation function in the initial PINN model as Lo et al. [9]: The initial crack is defined by the initial strain-history function, H(x, 0). This function is determined based on the closest distance between a point x in the domain and the initial crack, which represents the discrete crack [28]. The local strain-history functional approach allows for the specification of initial cracks in the system [26]. In particular, we set where G c is a material property known as the critical energy release rate, representing the energy needed to produce a fracture surface with a unit area. The parameter l 0 controls how the crack spreads. In Eqs. 24, B is a scalar parameter that controls the magnitude of the scalar history field and is calculated as: As is shown in Eq. 22, the output of the neural network is modified to meet the boundary conditions. As explained in Section 2, the solution to the crack problem is obtained by minimizing the free energy E. which consists of the elastic strain energy ψ e and the fracture surface energy ψ c . To use the variational energy based PINN approach to study the growth of fracture, the problem is formulated as: Minimize: E ψ e + ψ c subject to: Note that integration is needed to calculate ψ e and ψ c in Eq. 26 for obtaining the elastic strain energy and the surface energy over the whole domain. Figure 2 shows a schematic diagram of the proposed PINN framework. For simplicity, only one layer of hidden layers is explicitly shown.

Numerical examples 4.1 One dimensional cracked elastic bar
We are analyzing a bar with a crack located at the center (x = 0) and it is fixed at both ends (x = −1, x = 1). This bar is subjected to a sinusoidal load. The geometric setup is shown in Figure 3. For a given material, we take into account the elastic property E, material strength σ c , and fracture energy G c as the primary parameters. This implies that for the commonly used quadratic degradation function in Eq. 7, the length scale l 0 is not a separate variable and is valued as l 0 27G c E/256σ 2 c . We choose the following parameters for this model: σ c = 0.15MPa, G c = 2.7N/m, E = 1MPa, and l 0 = l p = 12.5 μm. The mesh size h in the vicinity of the crack is h = l 0 . The crack at center is imposed by the following initial strain-history function: where l 0 is the length scale parameter. The displacement field satisfies the Dirichlet boundary conditions, i.e., The analytical solution of the displacement field Schillinger et al. [29] is discontinuous and is given by the following equation: where u θ is the displacement field obtained as the output of the neural network. To optimize the proposed Physics-Informed Neural Network (PINN), we minimize the total variational energy of the system, which is defined in Eq. 20. To evaluate the precision of the results derived using this approach, we adopt two metrics: the relative error L rel 2 and the root mean square error (RMSE). Next, we study the effect of the proposed degradation function on the model. First, for the proposed degradation function, it can decouple between the length of the coupling field, l 0 , and the physical length, l p . By setting different q values, we can make l 0 multiple times larger than l p . For instance, l 0 = 100l p when q = 1.00155. Figure 4 shows the simulation results in cases of l 0 = l p = 0.0125, 1/10L 0 = l p = 0.00125 and 1/100L 0 = l p = 0.000125. To ensure an impartial evaluation, both methods were tested using the same neural network architecture and the number of integral points. The results from the proposed method were then compared to those obtained from the AT2-based PINN approach. Figures 4A,B show the displacement field u and phase field d obtained using l p = 0.0125. Visually, the obtained results overlap with the analytical solutions obtained using formulas Eqs. 29, 30. To quantify the accuracy of the   Figures 4E,F is actually due to a lack of resolution at the material length scale. Therefore increasing the integration points appropriately and adjusting the number of training iterations may narrow the difference. The results shown in Figures 4E,F are not perfect. However, as the linear elastic fracture mechanics must meet the constraint of the phase-field process on a small scale zones, the constraint of the phase-field process on a small scale zones also applies to the proposed modeling method. Nevertheless, the phase field length scale can be selected independent of the inherent material process zone length scale, provided that the material processing zone length scale is sufficiently refined.

Single-edge notched plate subjected to tension
In this case, we take into consideration a unit square plate with a horizontal crack running from the midpoint of the left outer edge to the center of the plate. The problem's geometry and boundary conditions are illustrated in Figure 5. The plate material has a tensile strength of σ c = 2.54MPa, a fracture energy of G c = 2.7N/ m, a Young's modulus of E = 282.69MPa, and a phase field length scale of l 0 = l p = 12.5 μm. The mesh size, h, near the crack is h = l 0 /4. The constant displacement increment applied during the computation is δv = 10 −3 mm. The crack path for the single-edge notched plate under tension was obtained using a fully connected neural network with three hidden layers, each of which has 50 neurons. The activation function for the first two layers adopt the tanh function, and for the last layer adopt the linear function. The initial crack was determined using the strain history functional described in Eq. 24. The boundary conditions imposed were Dirichlet boundary conditions.
where u and v are the solutions of the elastic field along the x and y-axes, respectively. To satisfy the Dirichlet boundary conditions, the neural network outputs for the elastic field are modified as: whereû andv are obtained from the neural network. Figure 6 illustrates the crack propagation process in case of l p = l 0 = 0.0125 mm. The result is very close to that of Goswami et al. [20]. As in the previous case, we first study the effect of the proposed degradation function on different physical length scale, l p for the same size model. The ratio between l 0 and l p is varied by setting different values of q. Figure 7 shows the simulation results for the cases l 0 = l p = 0.0125mm, l 0 = 4l p = 0.05 mm and l 0 = 10l p = 0.125 mm. In the above simulation, we observe that when the length of the phase field increases, the number of required Gaussian integration points decreases considerably, when l 0 = l p = 0.0125mm, the number of required integration points is: 20*20*16 (20 represents the mesh of elements and 16 represents the number of required integration points per element), when l 0 = 4l p = 0.05mm, the number of required integration points is: Frontiers in Physics frontiersin.org 15*15*16, When l 0 = 10l p = 0.125mm, the number of integration points required is: 12*12*16, and the same simulation speed is also increased. Therefore, assuming that the model is enlarged and the l p is unchanged, does it mean that we can simulate the harsh conditions in the large-scale model where the physical length scale, l p has to be small. We will therefore next investigate increasing the size of the model while keeping the length scale of the physical process zone constant. The previous example in Section 4.1 shows that decoupling the physical length scale from the phase field length is achievable albeit with limitations. The proposed degradation function allows to solve large-scale problem with much less computational cost. Suppose we multiply the length of the model edge by a factor of 4 and 10. In the meanwhile, the value of q is modified so that l 0 = 4l p and l 0 = 10l p . The value of σ c remains constant all the time. The proposed degradation function allows to increase the size of l 0 and decrease the mesh, 1 h (which is closely related to the size of l 0 , usually h = l 0 /4), while keeping the physical length scale l p constant. Simulation results for l 0 = 4l p and l 0 = 10l p are shown in figures Figure 8 and Figure 9. With the classical degradation function, if we want to increase the model by a factor of 10, we would need to add more integration points, whereas with the proposed degradation function, we can keep the number of integration points almost unchanged while increasing the model size by a factor of 10. We observe that the crack path shape is accurately captured by the larger length scale of the phase field. This adjustment to the phase field model formulation, therefore, allows for an increased sample size without the need for a significant increase in the computational volume of the length scale of the material processing zone.

Symmetrically double-edge notched plate sujected to tension
To verify the proposed method, we add another case, in which a double-edge notched plate subjected to a tensile load (see Figure 10) is considered. The material parameters are σ c = 2.54Mpa, G c = 2.7*10N/m, E = 282.69Mpa, and l 0 = l p = 12.5 μm.
To determine the crack path, we employed a fully connected neural network that includes four hidden layers, each of which has 50 neurons. The first three layers use the tanh activation function, while the final layer uses the linear activation. Both of cracks were initiated using the strain history functional, and the Dirichlet boundary conditions are specified: The solution for the elastic field along the x and y-axes are represented by u and v, respectively. In order to determine the crack path, a constant displacement increment of δu = 0.5p10 −4 mm has been applied. The output from the neural network for the elastic field is adjusted to conform with the Dirichlet boundary conditions: whereû andv are obtained from the neural network. The value of q is set to be 200. Figure 11 show the predicted crack propagation as a function of the tensile displacement. The predicted results are Frontiers in Physics frontiersin.org compared with the results in Goswami et al. [20], and our findings are highly consistent with the results reported in the literature. We multiply the lengths of the model edges by factors of 4 and 10, and modify the value of q so that l 0 = 4l p and l 0 = 10l p . The value of σ c is kept constant. The simulation results are shown in Figure 12 and Figure 13. Compared to the classical degradation function, there is no substantial increase in the running time when the model is expanded by a factor of 10 (mainly a small increase in the number of neural network iterations). We can find that with by introducing the length-scale separating degradation function to the PINN-based phase field model, maintaining a given mesh density sufficient to solve a larger scale model.

Conclusion and outlook
In this paper, we enhanced PINN-based phase field method (Goswami et al. [20], Goswami et al. [21]) with the degradation functions proposed by Lo et al. [9] in order to simulate fracture propagation in large structures. In conventional phase field method, the length scale of the phase field length scale is set to equate the physical length scale, which is very small compared to the structure size. Because a very refined mesh is needed to represent the diffuse fracture zone, it is impractical to apply PINN-based phase field method to large-scale models. To solve this problem, we introduce to PINN-based phase field simulation the degradation function that decouples the length scales of the phase field and the physical process zone. The merits and limitations of the modified phase field approach are discussed through several numerical examples. The proposed method combines the advantages of data driven and mechanism based computational approaches and improved its efficiency in fracture simulation. Although these examples are only for 2D problems, it is expected that the method could be applied to 3D situations without issue. In the future, we will investigate the uncertainty quantification of fractures with the present method. In addition, we will extend the present method to the hydraulic fracture problems.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.