# Phase-field modeling of fracture with physics-informed deep learning

M. Manav<EMAIL_ADDRESS>R. Molinaro<EMAIL_ADDRESS>S.
Mishra<EMAIL_ADDRESS>[ L. De Lorenzis
<EMAIL_ADDRESS>[

###### Abstract

We explore the potential of the deep Ritz method to learn complex fracture
processes such as quasistatic crack nucleation, propagation, kinking,
branching, and coalescence within the unified variational framework of phase-
field modeling of brittle fracture. We elucidate the challenges related to the
neural-network-based approximation of the energy landscape, and the ability of
an optimization approach to reach the correct energy minimum, and we discuss
the choices in the construction and training of the neural network which prove
to be critical to accurately and efficiently capture all the relevant fracture
phenomena. The developed method is applied to several benchmark problems and
the results are shown to be in qualitative and quantitative agreement with the
finite element solution. The robustness of the approach is tested by using
neural networks with different initializations.

###### keywords:

Phase-field fracture , Physics-informed machine learning , Deep Ritz method ,
Non-convex optimization , Crack nucleation , Crack propagation

url]https://camlab.ethz.ch/

url]https://compmech.ethz.ch/

[1]organization=Department of Mechanical and Process Engineering, ETH Zurich,
addressline=Tannenstrasse 3, city=Zurich, postcode=8092, country=Switzerland

[2]organization=Seminar for Applied Mathematics, Department of Mathematics,
ETH Zurich, addressline=Ramistrasse 101, city=Zurich, postcode=8092,
country=Switzerland

##  1 Introduction

Predictive modeling and simulation of fracture processes is critically
important in design for many engineering applications and still poses
significant challenges. In the past two decades, phase-field fracture modeling
[1, 2] has emerged as a game changer; it offers a unified variational
framework and, upon discretization, a flexible computational method to
encompass various aspects of fracture including crack nucleation, propagation,
merging, and branching [2, 3, 4], while obviating the need for crack tracking
procedures as well as for ad-hoc parameters and criteria. For these reasons,
it has found applications and extensions to a variety of fracture problems
including thermal cracks [5], ductile fracture [6], drying-induced cracks [7],
and hydraulic fracture [8], among many others. Furthermore, it has been
demonstrated to accurately capture topologically complex crack paths, even in
three dimensions [5], while delivering quantitatively accurate predictions
[9].

The first variational phase-field model for brittle fracture was derived by
regularizing the free-discontinuity problem obtained from recasting Griffith’s
fracture criterion in a variational form [10, 1]. Regularization is performed
by introducing an additional field, known as the phase field, ranging between
0 (representing the intact material) and 1 (representing the fully cracked
material). Localization of this field in a narrow band, whose width is
controlled by a length-scale parameter inherent to the model, gives rise to a
smeared representation of a crack. Later, it was shown that a more flexible
but substantially analogous variational modeling framework can be constructed
from a special family of gradient damage models [11], whereby the phase field
is interpreted as a damage variable and is therefore postulated to be
irreversible. In the computational setting, the regularization of a sharp
crack to a smeared representation overcomes the computational challenges
encountered in handling an a priori unknown crack path. The crack is now
represented implicitly by the continuous phase field; along with the
displacement field, this field is obtained by solving the non-linear coupled
partial differential equations (PDEs) stemming from the stationarity
conditions of the governing energy functional. However, two challenging
aspects remain: the energy functional is non-convex, which needs special
attention in the numerics [12], and the resolution of the (small) length scale
in the discretized setting is computationally expensive. This limitation
restricts the applicability of the approach to e.g. optimization, uncertainty
quantification and inverse problems.

In recent years, considerable effort has been directed towards leveraging
advances in machine learning to accurately model physical phenomena [13];
various physics-informed deep learning approaches [14] have been employed to
solve PDEs and variational problems [15, 16, 17, 18]. One prominent approach,
known as physics-informed neural networks (PINNs), learns the solution of a
PDE by the unsupervised (or semi-supervised) training of a neural network (NN)
[14, 16, 19, 20]. PINNs build upon the universal approximation property of NNs
[21, 22, 23], i.e. the ability of NNs to approximate any continuous function,
which allows for the use of NNs as the ansatz space for PDE solutions.
Notably, PINNs can seamlessly incorporate any available data and compute
solutions to forward and inverse problems with the same NN architecture. They
have demonstrated accurate approximation of solutions to a variety of PDEs and
inverse problems, also in solid mechanics [24, 25, 26, 27, 28, 29, 30, 31].
Rigorous estimates of the generalization error of PINNs approximations have
also been established [32, 20, 33]. Furthermore, multiple variations of PINNs
have been proposed, such as, hp-VPINNs [17], cPINNs [34], xPINNS [35], FBPINNs
[36], among others, to tailor them for different problems and to improve their
performance.

In PINNs, training of the NN representing the solution field involves the
minimization of a loss function comprising the PDE residual and additional
terms, which enforce initial and boundary conditions and possibly other
constraints. An alternative approach known as the deep Ritz method (DRM) [18,
37, 38], directly minimizes the energy functional instead of the PDE residual.
The DRM is particularly useful in obtaining stable solutions when the
governing energy functional is non-convex, as in phase-field fracture
modeling. Additionally, since the computation of the energy requires a lower
order of derivatives compared to the corresponding strong-form PDE, the DRM is
expected to be computationally less expensive [37]. On the other hand, energy
calculation requires integration over the entire domain, and only a few
studies have investigated quadrature rules for NNs [39].

PINN approaches have been developed to solve phase-field problems (Cahn-
Hilliard and Allen-Cahn equations) involving the phase field as the only
scalar unknown [40, 41]. In contrast, phase-field fracture modeling involves
two coupled fields: the displacement (vector) field and the phase field.
Furthermore, while the steep variation of the phase field at localization is
governed by the regularization length scale, the displacement field undergoes
a much sharper variation wherever the phase field localizes. These differences
lead to additional and unique challenges in approximating the solution of
phase-field fracture problems using NNs. While the DRM has shown some promise,
its applicability has been demonstrated primarily in problems involving crack
propagation with simple crack paths [42, 43]. Crack initiation in the absence
of a notch has only been demonstrated by assuming a very large value for the
regularization length scale [42]. Additionally, these approaches have required
a large number of quadrature points, exceeding the number needed to resolve
the regularization length scale. Notably, a quantitative assessment of the
solution accuracy is also lacking. The applicability of approaches other than
the DRM has been partially investigated in [42, 44]. An operator learning
approach [45, 46], namely variational DeepONets, has also been applied to
predict the crack path in a problem involving crack propagation in quasi-
brittle materials [47]. However, the applicability of physics-informed deep
learning to phase-field fracture modeling hinges on the ability of NNs to
learn all the fracture processes accurately, especially considering that this
flexibility lies at the core of the phase-field approach and is also one of
the key reasons of its success.

In this work, we elucidate the challenges involved in learning the solution of
phase-field fracture problems, and investigate the design of NNs and of their
training strategies aimed at learning crack initiation, propagation, kinking,
branching, and coalescence. We also address the challenges associated with
learning the solution fields with the same level of domain discretization as
required in finite element analysis (FEA). We systematically demonstrate the
accuracy of the learned solution by comparing it with the FEA solution.
Additionally, we test the robustness of our approach by training NNs with
different random initializations.

The paper is organized as follows. A brief review of the phase-field brittle
fracture formulation is presented in Section 2. In Section 3, the DRM is
summarized, and the challenges in learning the solution of a phase-field
fracture problem are discussed along with the design of the NNs and their
training strategy. This is followed by numerical examples of crack nucleation,
propagation, kinking, branching, and coalescence in Section 4. Finally, we
summarize the main conclusions in Section 5.

##  2 Phase-field model of brittle fracture

In this section, we summarize the governing equations of the phase-field model
of brittle fracture which we aim at solving with the deep learning approach
proposed in the following sections.

###  2.1 Basics of the formulation

Let $\Omega\subset\mathbb{R}^{d},~{}d\in\\{1,2,3\\}$ be a bounded open domain
occupied by a $d$-dimensional body. $\Gamma_{D,0}$, $\Gamma_{D,1}$ and
$\Gamma_{N}$ are disjoint sections of the boundary $\partial\Omega$ of
$\Omega$ with prescribed homogeneous Dirichlet, non-homogeneous Dirichlet and
Neumann boundary conditions, respectively. The body is assumed to be linear
elastic with the strain energy density given by
$\Psi(\bm{\varepsilon})=\frac{1}{2}\lambda{\rm
tr}^{2}(\bm{\varepsilon})+\mu{\rm tr}(\bm{\varepsilon}\cdot\bm{\varepsilon})$,
where $\lambda$ and $\mu$ are the Lamé constants, and
$\bm{\varepsilon}=\frac{1}{2}(\bm{\nabla u}+\bm{\nabla u}^{T})$ is the
infinitesimal strain tensor, with $\bm{u}:\Omega\to\mathbb{R}^{d}$ as the
displacement field. The total energy functional for the body acted upon by
body force $\bar{\bm{b}}$ and surface traction $\bar{\bm{t}}$ is given by:

$\displaystyle\mathcal{E}(\bm{u},\alpha)=$
$\displaystyle\underbrace{\int_{\Omega}\left(g(\alpha)\Psi^{+}(\bm{\varepsilon}(\bm{u}))+\Psi^{-}(\bm{\varepsilon}(\bm{u}))\right){\rm
d}\Omega}_{\mathcal{E}^{el}}+\underbrace{\frac{G_{c}}{c_{w}}\int_{\Omega}\left(\frac{w(\alpha)}{l}+l|\bm{\nabla}\alpha|^{2}\right){\rm
d}\Omega}_{\mathcal{E}^{d}}$
$\displaystyle-\int_{\Omega}\bar{\bm{b}}\cdot\bm{u}\,{\rm
d}\Omega-\int_{\Gamma_{N}}\bar{\bm{t}}\cdot\bm{u}\,{\rm d}s.$ (1)

Here $\mathcal{E}^{el}$ and $\mathcal{E}^{d}$ represent the elastic and damage
parts of the energy, respectively. $\alpha:\Omega\to[0,~{}1]$ is the phase
field (or damage field), with $\alpha=0$ representing the undamaged state and
$\alpha=1$ the fully damaged state. The function $g(\alpha)$, known as
degradation function, modulates the degradation of the elastic energy with
increasing damage. $w(\alpha)$ is also known as local dissipation function, as
it dictates the energy dissipated through damage in the unit volume of the
material. The parameter $l$ with $0<l\ll{\rm diam}(\Omega)$ is the
regularization length controlling the thickness of the transition zone between
undamaged and fully damaged material regions; $G_{c}$ represents the fracture
toughness of the material, and $c_{w}$ is a normalization constant given by
$c_{w}=4\int_{0}^{1}\sqrt{w(t)}dt$. In this work, we adopt the following
choices [11]:

$\displaystyle{\rm AT1\,\,model}:\quad g(\alpha):=(1-\alpha)^{2}+\eta,\quad
w(\alpha):=\alpha,\quad c_{w}=8/3,$ (2) $\displaystyle{\rm AT2\,\,model}:\quad
g(\alpha):=(1-\alpha)^{2}+\eta,\quad w(\alpha):=\alpha^{2},\quad c_{w}=2.$ (3)

with $\eta=o(l)$. The elastic strain energy density $\Psi$ is decomposed into
active (i.e. damage-driving) and inactive (i.e. damage-resisting) parts,
$\Psi^{+}$ and $\Psi^{-}$, respectively, and the degradation is applied only
to $\Psi^{+}$ [48, 49, 50, 51]. While several energy decompositions are
available [48, 49, 50, 52], in this work we focus on the volumetric-deviatoric
decomposition [48] whereby

$\displaystyle\Psi^{+}=\frac{1}{2}K\langle\rm
tr(\bm{\varepsilon})\rangle_{+}^{2}+\mu{\rm tr}(\mathbf{e}\cdot\mathbf{e}),$
$\displaystyle\Psi^{-}=\frac{1}{2}K\langle\rm
tr(\bm{\varepsilon})\rangle_{-}^{2},$ (4)

with $\langle a\rangle_{+}={\rm max}\\{0,a\\}$, $\langle a\rangle_{-}={\rm
min}\\{0,a\\}$, the bulk modulus $K=\lambda+\frac{2}{3}\mu$, and the
deviatoric strain tensor
$\mathbf{e}=\bm{\varepsilon}-\frac{tr(\bm{\varepsilon})}{3}\mathbf{I}$, where
$\mathbf{I}$ is the third-order identity tensor.

###  2.2 Time-discrete evolution problem

In the time-discrete setting of the quasi-static evolution problem, the state
of the system at the (pseudo-)time or loading step $n\geq 1$ is obtained as
the solution of the minimization problem

$\operatorname*{argmin}_{\bm{u},\alpha}\\{\mathcal{E}_{n}(\bm{u},\alpha):~{}\bm{u}\in\bm{V}_{\bar{\bm{u}}_{n}},~{}\alpha\in\mathcal{D}_{\alpha_{n-1}}\\},$
(5)

where

$\mathcal{E}_{n}(\bm{u},\alpha)=\mathcal{E}^{el}(\bm{u},\alpha)+\mathcal{E}^{d}(\alpha)-\int_{\Omega}\bar{\bm{b}}_{n}\cdot\bm{u}\,{\rm
d}\Omega-\int_{\Gamma_{N}}\bar{\bm{t}}_{n}\cdot\bm{u}\,{\rm d}s,$ (6)
$\bm{V}_{\bar{\bm{u}}_{n}}:=\\{\bm{u}\in(H^{1}(\Omega))^{d}:~{}\bm{u}=\bm{0}~{}{\rm
on}~{}\Gamma_{D,0},~{}\bm{u}=\bar{\bm{u}}_{n}~{}{\rm on}~{}\Gamma_{D,1}\\}$
(7)

denotes the space of kinematically admissible displacement fields, and

$\mathcal{D}_{\alpha_{n-1}}:=\\{\alpha\in
H^{1}(\Omega):~{}\alpha\geq\alpha_{n-1}~{}{\rm in}~{}\Omega\\}$ (8)

denotes the space of admissible phase fields. The condition
$\alpha\geq\alpha_{n-1}$ represents the irreversibility of the phase field
evolution. The first-order necessary conditions of the minimization problem
deliver the strong form of the governing equations

$\displaystyle\bm{\nabla}\cdot\bm{\sigma}+\bar{\bm{b}}_{n}=0,$ (9)
$\displaystyle\alpha-\alpha_{n-1}\geq 0,$ (10) $\displaystyle
g^{\prime}(\alpha)\Psi^{+}+\frac{G_{c}}{c_{w}}\left(\frac{w^{\prime}(\alpha)}{l}-2l\Delta\alpha\right)\geq
0,$ (11)
$\displaystyle(\alpha-\alpha_{n-1})\left[g^{\prime}(\alpha)\Psi^{+}+\frac{G_{c}}{c_{w}}\left(\frac{w^{\prime}(\alpha)}{l}-2l\Delta\alpha\right)\right]=0,$
(12)

all valid in $\Omega$, and represent the local equilibrium of forces, damage
irreversibility, damage criterion and loading-unloading conditions,
respectively. The Cauchy stress is defined as
$\bm{\sigma}=g(\alpha)\frac{\partial\Psi^{+}}{\partial\bm{\varepsilon}}+\frac{\partial\Psi^{-}}{\partial\bm{\varepsilon}}$,
and $()^{\prime}=\frac{\partial}{\partial\alpha}()$. The whole set of boundary
conditions reads as follows:

$\displaystyle\bm{u}=\bm{0}\quad$ $\displaystyle{\rm on~{}}\Gamma_{D,0},$
$\displaystyle\bm{u}=\bar{\bm{u}}_{n}\quad$ $\displaystyle{\rm
on~{}}\Gamma_{D,1},$ $\displaystyle\bm{\sigma}\bm{n}=\bar{\bm{t}}_{n}\quad$
$\displaystyle{\rm on~{}}\Gamma_{N},$ $\displaystyle\alpha-\alpha_{n-1}\geq
0,\,\,\nabla\alpha\cdot\bm{n}\geq
0,\,\,(\alpha-\alpha_{n-1})(\nabla\alpha\cdot\bm{n})=0\quad$
$\displaystyle{\rm on~{}}\Gamma.$ (13)

###  2.3 Irreversibility of the phase field

One way to enforce the irreversibility of the phase field in the computations
is to add to the energy functional in (1) the following energetic penalty
[53]:

$\mathcal{E}^{ir}(\alpha)=\int_{\Omega}\frac{1}{2}\gamma_{ir}\langle\alpha-\alpha_{n-1}\rangle_{-}^{2}{\rm
d}\Omega,$ (14)

where $\gamma_{ir}$ is the penalty parameter, whose value can be chosen as
suggested in [53]:

$\displaystyle\gamma_{ir}=\frac{G_{c}}{l}\frac{27}{64\mathtt{TOL_{ir}^{2}}}\quad\quad\quad\,\,({\rm
AT1~{}model}),$ (15)
$\displaystyle\gamma_{ir}=\frac{G_{c}}{l}\left(\frac{1}{\mathtt{TOL_{ir}^{2}}}-1\right)\quad({\rm
AT2~{}model}).$ (16)

Here $0<\mathtt{TOL_{ir}}\leq 1$ is the prescribed irreversibility tolerance
threshold; in this work, we set it to $\mathtt{TOL_{ir}}=5\times 10^{-3}$.

###  2.4 Non-dimensionalization scheme

Non-dimensionalization not only reduces the number of relevant physical
parameters, but also aids in learning by ensuring that the inputs are $\sim
1$. So, we non-dimensionalize the energy functional using the following
scheme:

$\displaystyle\tilde{\bm{x}}=\frac{\bm{x}}{L},\quad\tilde{l}=\frac{l}{L},\quad\tilde{\bm{u}}=\frac{\bm{u}}{L}\left(\frac{G_{c}}{El}\right)^{-1/2},$
$\displaystyle\tilde{\lambda}=\frac{\lambda}{E},\quad\tilde{\mu}=\frac{\mu}{E},\quad{\rm
and}\quad\tilde{\mathcal{E}}=\frac{l\mathcal{E}}{L^{3}G_{c}},$ (17)

where $\bm{x}\in\Omega$ ($\tilde{\bm{x}}\in\tilde{\Omega}$) is a material
point coordinate, $L$ is the characteristic length of the body, $E$ is the
Young’s modulus, and $\sqrt{\frac{G_{c}}{El}}$ is of the order of the critical
strain for crack nucleation in a 1D bar [11]. Then the dimensionless form of
the energy functional in (1), with no surface tractions and body forces, and
including the energetic penalty to enforce the irreversibility of the phase
field, reads

$\displaystyle\tilde{\mathcal{E}}(\tilde{\bm{u}},\alpha)=$
$\displaystyle\underbrace{\int_{\tilde{\Omega}}\left(g(\alpha)\tilde{\Psi}^{+}(\tilde{\bm{\varepsilon}}(\tilde{\bm{u}}))+\tilde{\Psi}^{-}(\tilde{\bm{\varepsilon}}(\tilde{\bm{u}}))\right){\rm
d}\tilde{\Omega}}_{\tilde{\mathcal{E}}^{el}}$
$\displaystyle+\underbrace{\frac{1}{c_{w}}\int_{\tilde{\Omega}}\left(w(\alpha)+\tilde{l}^{2}|\bm{\tilde{\nabla}}\alpha|^{2}\right){\rm
d}\tilde{\Omega}}_{\tilde{\mathcal{E}}^{d}}$
$\displaystyle+\underbrace{\int_{\Omega}\frac{1}{2}\tilde{\gamma}_{ir}\langle\alpha-\alpha_{n-1}\rangle_{-}^{2}{\rm
d}\tilde{\Omega}}_{\tilde{\mathcal{E}}^{ir}},$ (18)

where $\bm{\tilde{\nabla}}=L\bm{\nabla}$,
$\tilde{\bm{\varepsilon}}=\bm{\varepsilon}\left(\frac{G_{c}}{El}\right)^{-1/2}$,
$\tilde{\Psi}^{\pm}=\frac{l{\Psi}^{\pm}}{G_{c}}$,
$\tilde{\gamma}_{ir}=\frac{l\gamma_{ir}}{G_{c}}$,
$\tilde{\mathcal{E}}^{el}=\frac{l\mathcal{E}^{el}}{L^{3}G_{c}}$,
$\tilde{\mathcal{E}}^{d}=\frac{l\mathcal{E}^{d}}{L^{3}G_{c}}$, and
$\tilde{\mathcal{E}}^{ir}=\frac{l\mathcal{E}^{ir}}{L^{3}G_{c}}$. Also, the
non-dimensional stress follows as
$\tilde{\bm{\sigma}}=g(\alpha)\frac{\partial\tilde{\Psi}^{+}}{\partial\tilde{\varepsilon}}+\frac{\partial\tilde{\Psi}^{-}}{\partial\tilde{\varepsilon}}$.
Furthermore, leveraging the relation among Lamé constants, Young’s modulus and
Poisson’s ratio $\nu$, we obtain $\tilde{\lambda}=\frac{\nu}{(1+\nu)(1-2\nu)}$
and $\tilde{\mu}=\frac{1}{2(1+\nu)}$. Notably, as a consequence, the only
parameters in the non-dimensional energy are $\nu$ and $\tilde{l}$. For
simplicity of notation, from here onwards, we only use non-dimensional
quantities but omit the $\tilde{(\cdot)}$ symbol.

##  3 A Deep Ritz method for phase-field fracture modeling

###  3.1 The Deep Ritz method

In the DRM, as in any other physics-informed deep learning approach for
solving PDEs and variational problems, the solution field is represented by an
NN, frequently a feedforward NN (also known as a multilayer perceptron). A
feedforward NN is a composition of affine transformations and scalar non-
linear activation functions. Given an input $\bm{y}\in\mathbb{R}^{n}$, the NN
maps it to the output $\bm{z}_{\theta}\in\mathbb{R}^{m}$ as follows

$\bm{z}_{\theta}(\bm{y})=C_{K}\circ\sigma\circ
C_{K-1}\cdot\cdot\cdot\circ\,\,\sigma\circ C_{1}(\bm{y}),$ (19)

where $C_{k}$ (with $1\leq k\leq K$) represents an affine transformation, also
known as the $k^{th}$ layer of the NN. Specifically,
$C_{k}\bm{z}_{k}=W_{k}\bm{z}_{k}+b_{k}$, where $W_{k}$ and $b_{k}$ are
trainable weights and biases and $\bm{z}_{k}$ is the input to the layer. The
symbol $\circ$ represents the composition of functions, and $\sigma$ is a
scalar activation function. Commonly used activation functions are the ReLU,
sigmoid, and hyperbolic tangent ($\tanh$) functions [54]. $\theta$ denotes the
set of all trainable parameters of the NN, i.e.
$\theta=\\{W_{k},~{}b_{k}\\},\forall~{}1\leq k\leq K$.

Training of the NN aims to find $\theta$ such that $\bm{z}_{\theta}$
approximates the target solution. It involves constructing an appropriate loss
function and then minimizing the loss with respect to the parameters in
$\theta$. The minimization is typically performed using stochastic gradient
decent algorithms such as Adam [55] or higher-order optimization algorithms
such as L-BFGS [56].

For our problem, we define the displacement and phase fields obtained from the
NN as $\bm{u}_{\theta}(\bm{x})$ and $\alpha_{\theta}(\bm{x})$. Then the loss
function is expressed as

$\mathcal{L}=\log(\mathcal{E}_{\theta}),$ (20)

where $\mathcal{E}_{\theta}=\mathcal{E}(\bm{u}_{\theta},\alpha_{\theta})$ is
obtained from (18). Taking the $\log$ of the energy in (20) makes the loss
function $\sim 1$, allowing us to set the same weight for weight
regularization (to be discussed later) and the same stopping criterion for the
optimizer for all the problems investigated in this work.

Note that, in our experience, the enforcement of the Dirichlet boundary
conditions using a soft constraint in DRM is challenging. Hence, we enforce
them by constructing an ansatz, following the approach employed in [14], see
Section 3.2.1. For a general domain shape for which an ansatz is not feasible,
general approaches for the enforcement of boundary conditions as in [57, 58]
can be employed.

###  3.2 Construction and training of the NN

Our aim is to obtain an NN-based approximation of the solution fields
employing the DRM in learning. This approach entails two key challenges. The
first challenge lies in ensuring that the NN-based approximation of the energy
surface, $\mathcal{E}_{\theta}$, is a faithful projection of the energy
landscape in the space of NN-parameters. Specifically, it is crucial that the
energy barriers and minima are appropriately represented. In turn, the
representation of the energy surface relies on three ingredients: the
construction of the NN(s) representing the fields, the computation of the
gradients of the fields, and the integration over the domain. The second
challenge lies in the choice of the optimization algorithm. This algorithm
needs to respect the energy barriers and reach the correct local energy
minimum regardless of the complexity of the solution associated with that
energy minimum. As follows, we discuss all these aspects and refer to some
Appendices for additional details.

Figure 1: Scheme of our deep Ritz method for the two-dimensional case (left);
function constraining the phase field between 0 and 1 (right)

####  3.2.1 Choice of the NN architecture

Physically distinct fields such as the displacement field and the phase field
are expected to be learned efficiently when represented by independent NNs.
However, we find that an intimate coupling between the fields, facilitated by
representing the displacement and phase fields as the outputs of the same NN,
helps in learning the nucleation of a crack as shown in A. Hence, we construct
an NN which takes as input the position vector of a point in the domain (2
components in 2D, 1 in 1D) and delivers outputs which correspond to the
components of the displacement vector and the value of the phase field at that
point (3 outputs in 2D, 2 in 1D), see Figure 1. The non-dimensionalization
scheme outlined in Section 2.4 ensures that the inputs are $\sim 1$. In
constructing an ansatz to apply Dirichlet boundary conditions, we ensure that
the outputs of the NN corresponding to the displacement field are also
approximately $\sim 1$ as follows:

$\displaystyle u_{\theta}=f_{u}(\bm{x},\hat{u}_{\theta})U_{p},$ (21)
$\displaystyle v_{\theta}=f_{v}(\bm{x},\hat{v}_{\theta})U_{p},$ (22)

where $\hat{u}_{\theta}$ and $\hat{v}_{\theta}$ are the direct NN outputs
corresponding to the displacement components, and $U_{p}$ is the prescribed
displacement. $f_{u}$ and $f_{v}$ are designed for each problem. To ensure
that $\alpha_{\theta}\in[0,1]$, the application of the sigmoid function to the
direct NN output corresponding to the phase field, $\hat{\alpha}_{\theta}$, is
an attractive choice; however, it turns out that its vanishing slope hinders
crack nucleation. Therefore, we design the following function, plotted in
Figure 1, to overcome this problem:

$\alpha_{\theta}=f_{\alpha}(\hat{\alpha}_{\theta})=\begin{cases}\frac{\hat{\alpha}_{\theta}}{4}+\frac{1}{2},&{\rm
for~{}}|\hat{\alpha}_{\theta}|\leq 2\\\ \beta(\hat{\alpha}_{\theta}+2),&{\rm
for~{}}\hat{\alpha}_{\theta}<-2\\\ \beta(\hat{\alpha}_{\theta}-2)+1,&{\rm
for~{}}\hat{\alpha}_{\theta}>2\end{cases}$ (23)

By setting $\alpha_{\theta}=f_{\alpha}(\hat{\alpha}_{\theta})$ in the 1D bar
problem, we find that the order of the smallest $\beta$ for which an NN can
learn the correct critical load for crack nucleation is $10^{-3}$. So, we set
$\beta=10^{-3}$ in all examples of this paper.

####  3.2.2 Choice of the activation function

Another important aspect of the NN design is the choice of the non-linear
activation function. In PINNs, smooth activations like $\tanh$ are frequently
employed. However, we find that an NN with this activation fails to learn the
correct crack path (see B). The evolution of the loss during training suggests
that an NN with $\tanh$ activation fails to represent the energy barrier in
the loss landscape which would prevent crack propagation in an incorrect
direction. Therefore, we employ ReLU activation with the following
modification:

$z_{k+1}=\max\\{0,~{}m_{k}(W_{k}z_{k}+b_{k})\\},$ (24)

where $m_{k}$ are learnable coefficients. While initializing the NN, we set
$m_{k}=m_{0}$, where $m_{0}$ is a constant. We observe that while the value of
$m_{k}$ does not change significantly during training, some choices of
$m_{0}$, as detailed in Section 4.1, facilitate learning the critical applied
displacement for crack initiation or propagation with higher accuracy.

####  3.2.3 Gradient computation and quadrature

Computation of the gradients poses another challenge. When employing
autodifferentiation to compute the gradients of the fields represented by an
NN with a smooth activation function like $\tanh$, the learned fields display
high-frequency oscillations (similar to those resulting from Gibbs phenomena)
leading to the localization of the stress in random directions near the crack
tip (see C). This can be viewed as an undesired modification of the energy
barrier, allowing for crack propagation in unexpected directions and leading
to an inaccurate crack path prediction.

While this is an additional argument against the use of tanh activation, the
problem is not solved by the use of the ReLU activation function. To address
it, we discretize the domain as we would do in FEA. The field values at the
nodes are obtained from the NN, whereas the field values and their gradients
within each element, in particular at the Gauss points, are computed using
shape functions like in FEA, as follows

$\displaystyle{\bm{u}}^{h}_{\theta}({\bm{x}})=\sum_{a=1}^{n}N^{a}\left(\bm{x}\right)\bm{u}_{\theta}^{a},\qquad$
$\displaystyle\alpha^{h}_{\theta}({\bm{x}})=\sum_{a=1}^{n}N^{a}\left(\bm{x}\right)\alpha_{\theta}^{a},$
$\displaystyle\nabla{\bm{u}}_{\theta}^{h}({\bm{x}})=\sum_{a=1}^{n}\bm{u}_{\theta}^{a}\otimes\nabla
N^{a}\left(\bm{x}\right),\qquad$
$\displaystyle\nabla\alpha^{h}_{\theta}({\bm{x}})=\sum_{a=1}^{n}\alpha_{\theta}^{a}\,\nabla
N^{a}\left(\bm{x}\right),$ (25)

where the superscript $h$ denotes the finite element approximation of the
fields, ${N}^{a}({\bm{x}})$ is the shape function corresponding to node $a$ of
the finite element discretization, and $\bm{u}_{\theta}^{a}$ and
$\alpha_{\theta}^{a}$ are the values of $\bm{u}_{\theta}$ and
$\alpha_{\theta}$ at the same node, respectively, with $a=1...n$ and $n$ as
the number of nodes. This approach to gradient computation makes the training
considerably faster.

We use the finite element discretization also for the purpose of approximating
the integral in (18) in order to compute the energy using the NN-based
representation of the displacement and phase fields. We evaluate the integral
employing Gauss quadrature at each element of the finite element
discretization.

####  3.2.4 Weight regularization

Weight regularization plays an important role in obtaining the correct
solution of a phase-field fracture problem using the same level of domain
discretization as needed in FEA. In the problem of crack nucleation in a 1D
bar in Section 4.2, when employing autodifferentiation for strain computation,
the NN learns incorrect solutions with approximately zero energy after crack
nucleation. An analysis of the governing PDEs (see D) reveals that the phase-
field model does not constrain the maximum slope of the displacement field
near the localized phase field, thus this slope depends on the expressivity of
the NN. Therefore, if the training points are not fine enough to resolve the
sharply changing displacement field expressible by the NN, and not only the
length scale $l$, the DRM yields incorrect solutions as shown in Figure 26.
Correct results can be obtained by a finer discretization of the domain than
needed to resolve the regularization length scale; this is however
computationally inefficient, as it requires an even finer resolution than
needed for FEA computations. Note that, while we first recognized the issue
when adopting autodifferentiation, the employment of gradient computation by
finite element discretization (as illustrated in Section 3.2.3) does not solve
this issue. Instead, a solution is to limit $\left(\frac{{\rm
d}u_{\theta}}{{\rm d}x}\right)_{max}$ by introducing weight regularization.
This leads to a computationally efficient approach, in which the
discretization of the domain needed for an accurate solution is comparable to
the one used in FEA. Therefore, we employ weight regularization throughout the
examples in this paper, see Section 4.1 for more details.

####  3.2.5 Choice of the optimization algorithm

Optimization algorithms play a key role in the training of an NN and selecting
one that respects energy barriers and reaches the correct energy minimum is
essential. In this study, we consider the L-BFGS algorithm [56] and the
resilient backpropagation (RPROP) algorithm [59]. While L-BFGS is the fastest
of the two, it is unable to reach a good energy minimum with the localized
phase field and the sharply varying displacement field (see E). In contrast,
the RPROP algorithm is observed to respect energy barriers and reach the
desired energy minimum. It should be noted that the performance of first-order
optimization algorithms like Adam for these problems was considerably poorer
compared to L-BFGS and RPROP, hence we do not consider them further in this
work.

####  3.2.6 Transfer learning

As loading is applied in incremental steps, we utilize transfer learning to
learn the solution quickly. The trained network from the previous loading step
is used as the initial network for training in the current loading step. We
observe that transfer learning enables learning the solution field in
typically fewer than 500 steps of the RPROP optimizer provided a load step
does not lead to a sudden crack propagation and associated jump in energy.

####  3.2.7 Nucleation and damage evolution

Our experience shows that crack nucleation is more challenging to learn than
crack propagation. In the AT1 model, the phase field does not evolve slowly
with increasing load but jumps abruptly from 0 to 1 at the location of crack
nucleation when the threshold loading for crack nucleation is reached. This
poses a formidable challenge for NNs to learn crack nucleation in 2D.
Therefore, for 2D problems involving crack nucleation we employ the AT2 model,
which allows for a slow evolution of the phase field to nonzero values before
it jumps to 1 at the threshold loading condition.

####  3.2.8 Summary

To summarize, our approach involves one NN with adaptive ReLU activation
expressing both the displacement and phase fields. Dirichlet boundary
conditions are enforced employing an ansatz which also ensures that the NN
outputs corresponding to the displacement field are approximately $\sim 1$. To
incentivize the phase field to lie between 0 and 1 without hindering crack
nucleation, we design a new function $f_{\alpha}(\hat{\alpha}_{\theta})$.
Field gradients are computed numerically like in FEA. To estimate the integral
in the computation of the energy, Gauss quadrature is employed for each
element in a discretized domain. In the training of the network, weight
regularization is employed to constrain the maximum strain expressible by the
NN. The RPROP algorithm is found to perform best in training the NN to learn
the solution fields. We also employ transfer learning to learn the solution
efficiently. Additionally, while the NN learns crack propagation, branching,
and coalescence using the AT1 model, it is able to learn crack nucleation in
2D only with the AT2 model.

##  4 Numerical examples

In this section, we first demonstrate crack nucleation in a 1D homogeneous bar
and an L-shaped panel. We then study crack propagation in a single-edge
notched (SEN) specimen under both tensile and shear loading conditions. To
illustrate crack branching, we again model a SEN specimen under shear loading
without strain energy decomposition. Crack coalescence is demonstrated by
subjecting a specimen with three preexisting cracks to tensile loading.

Since these examples are taken from the phase-field fracture literature, for
each example we use the same material properties as in the corresponding
paper. This not only enables an easier comparison with the results of the
referenced papers; it also demonstrates that the developed approach works
robustly for different material properties (and geometries). Throughout the 2D
examples, we assume plane-strain conditions.

###  4.1 Numerical setup

For the 1D problem, we use an NN of 4 hidden layers and 50 neurons in each
layer. The coefficient in the activation function is set to $m_{0}=1$, and we
do not include $m_{k}$ among the learnable parameters. The L-BFGS optimizer is
found to be sufficient to learn the solution in this case. Eight NNs with
different initializations obtained by using Xavier initialization [60] with
random seeds are trained for this problem.

For the 2D problems, we use an NN of 8 hidden layers and 400 neurons in each
layer. The size of the network is chosen based on the numerical experiment
conducted for the problem of crack nucleation in an L-shaped panel (see F). We
set $m_{0}=2$ for crack initiation in the L-shaped panel, and $m_{0}=3$ for
the remaining 2D problems. We train 8 NNs with different initializations for
each problem employing Xavier initialization with random seeds. In the first
loading step, an NN is first trained with L-BFGS before training it with
RPROP; for all the subsequent loading steps only RPROP is used. For RPROP, the
learning rate is set to $10^{-5}$ with the smallest and largest step sizes of
$10^{-10}$ and 50, respectively. The remaining optimizer parameters are kept
as in the standard implementation in Pytorch. The stopping criterion is set to
correspond to a relative change in the loss function of less than $5\times
10^{-6}$ for 10 consecutive optimization steps, with a maximum of 10000
optimization steps. Additionally, we apply $L_{2}$ weight regularization with
a weight of $10^{-5}$. The Pytorch library [61] is used for implementation and
training.

In the 1D problem, we discretize the domain into elements of size $l/5$ and
assume linear shape functions. In the 2D problems, we utilize triangular
elements with linear shape functions. Only the regions of the domain where the
crack is expected to propagate are finely discretized with elements of size
$l/5$, to limit the computational cost. Away from these regions, the element
size smoothly increases up to $4l$. Discretization is performed using Gmsh
[62]. We use one Gauss point per element for integration.

The numerical setup for the reference FEA computations is detailed in G.

###  4.2 Crack nucleation in a 1D homogeneous bar

Figure 2: 1D bar problem: geometry and boundary conditions.

We first study crack nucleation in a linear elastic bar shown in Figure 2.
This problem lends itself to an analysis that sheds light on the reason for
the occurrence of incorrect solutions when the domain is not discretized
sufficiently finely.

In this case, the expression of the energy in (18) simplifies to the following
form:

$\displaystyle\mathcal{E}(u,\alpha)=\int_{-0.5}^{0.5}\frac{1}{2}g(\alpha)\left(\frac{{\rm
d}u}{{\rm d}x}\right)^{2}{\rm
d}x+\frac{1}{c_{w}}\int_{-0.5}^{0.5}\left(w(\alpha)+l^{2}\left|\frac{{\rm
d}\alpha}{{\rm d}x}\right|^{2}\right){\rm d}x.$ (26)

The only parameter in the above expression is $l$, and we set $l=0.05$. We
employ the AT1 model (see (2)). The following ansatz is employed to apply the
boundary conditions:

$\displaystyle u_{\theta}=[(x+0.5)(x-0.5)\hat{u}_{\theta}+(x+0.5)]U_{p},$
$\displaystyle\alpha_{\theta}=(x+0.5)(x-0.5)\hat{\alpha}_{\theta}.$ (27)

Note that, since fracture at the boundaries requires less energy compared to
fracture inside the domain, the NN (as well as the FEA) has a tendency to
choose the solution with fracture at the boundary. To suppress this solution,
we set $\alpha_{\theta}=0$ at the boundaries using the above ansatz.

Figure 3: 1D bar: elastic and fracture energies vs prescribed displacement
from NN and FEA.

Figure 4: 1D bar: displacement and phase fields before and after crack
nucleation from NN and FEA.

For the 1D bar, the energies obtained from the NN are compared with the FEA
solution in Figure 3, showing close agreement. Additionally, the displacement
and phase fields before and after crack nucleation are compared with the FEA
solution in Figure 4, again showing that the DRM achieves a quite accurate
solution.

###  4.3 Crack nucleation in an L-shaped panel

Figure 5: L-shaped panel: geometry and boundary conditions.

This problem illustrates crack nucleation in 2D. The geometry of the panel and
the boundary conditions are as shown in Figure 5, whereas the material
properties are set to $\nu=0.18$ and $l=0.01$ [53]. To prescribe homogeneous
and nonhomogeneous Dirichlet boundary conditions, we construct distance
functions $d_{1}(x,y)$ and $d_{2}(x,y)$, respectively, such that $d_{1}(x,y)$
equals 1 at $y=-0.5$ and smoothly decreases to 0 away from it, and
$d_{2}(x,y)$ equals 1 at the section of the boundary where $U_{p}$ is applied
and smoothly decreases to 0 away from it (see H for more details):

$\displaystyle u_{\theta}=[(1-d_{1}(x,y))\hat{u}_{\theta}]U_{p},$
$\displaystyle
v_{\theta}=[(1-d_{1}(x,y))(1-d_{2}(x,y))\hat{v}_{\theta}+d_{2}(x,y)]U_{p},$
$\displaystyle\alpha_{\theta}=f_{\alpha}(\hat{\alpha}_{\theta}).$ (28)

Figure 6: L-shaped panel: displacement and phase fields at $U_{p}=1.2$ from NN
and FEA. Figure 7: L-shaped panel: elastic and fracture energies vs prescribed
displacement from NN and FEA.

As discussed in Section 3.2, the NN finds it difficult to learn crack
nucleation when the AT1 model is employed; therefore, we use the AT2 model in
this problem. Figure 6 compares the solution obtained from the NN with the FEA
solution at $U_{p}=1.2$. The crack path in the NN solution closely resembles
that in the FEA solution, although the FEA solution shows comparatively
sharper turning of the crack. Note that, due to the non-convexity of the
governing energy functional, even the FEA solution is in general not unique,
and changes in the obtained crack patterns may be induced by numerical
perturbations as discussed in [63]. Figure 7 illustrates the mean and standard
deviation of the energies obtained from 8 trained NNs with different
initializations, and compares them with the energies obtained from FEA. The
energies are in close agreement before crack nucleation, and the NN accurately
captures the critical $U_{p}$ at which nucleation takes place. However, in the
subsequent stage the damage energy obtained from the NN is higher than in FEA.
This discrepancy results from the difficulty in learning a curving crack and
can be improved by reducing the loading step sizes and imposing a stricter
stopping criterion for the optimizer, which however would increase the
computational cost.

###  4.4 Crack evolution in a notched specimen

Figure 8: SEN specimen: geometry and boundary conditions (left); phase field
from NN at the smallest applied tensile load (right).

The geometry of the SEN specimen and boundary conditions are as shown in
Figure 8, while the material properties are $\nu=0.3$ and $l=0.01$ [53]. We
define the following ansatz to apply the boundary conditions and to compute
$\alpha$:

$\displaystyle
u_{\theta}=[(y+0.5)(0.5-y)\hat{u}_{\theta}+(y+0.5)\cos(\omega)]U_{p},$
$\displaystyle
v_{\theta}=[(y+0.5)(0.5-y)\hat{v}_{\theta}+(y+0.5)\sin(\omega)]U_{p},$
$\displaystyle\alpha_{\theta}=f_{\alpha}(\hat{\alpha}_{\theta}).$ (29)

We define the notch in the sample by setting an initial phase field as
described in I; the NN learns the notch as a localized phase field as shown in
Figure 8.

####  4.4.1 Tensile loading

Figure 9: SEN specimen under tensile loading ($\omega=\pi/2$): displacement
and phase fields at $U_{p}=0.2$ from NN and FEA.

With this problem, we demonstrate the ability of the NN to learn crack
propagation. We subject the SEN specimen to tensile loading, i.e. the loading
angle is $\omega=\pi/2$, and adopt the AT1 model. Figure 9 compares the fields
obtained from NN and FEA at $U_{p}=0.2$, showing close agreement. In this
example, the simplicity of the crack path allows the NN to easily and
accurately learn the cracked solution. Figure 10 compares the energies
associated with the NN solution with the energies from FEA, exhibiting once
again a close agreement. Note that the comparatively higher standard deviation
near the critical load for crack propagation results from differences in the
prediction of the critical load for different initializations of the NN.

Figure 10: SEN specimen under tensile loading ($\omega=\pi/2$): elastic and
fracture energies vs prescribed displacement from NN and FEA.

####  4.4.2 Shear loading

This problem demonstrates the ability of the NN to learn crack kinking, i.e.
the change in direction of a crack. To this aim, we subject the SEN specimen
to shear loading, i.e. the loading angle is $\omega=0$, and again use the AT1
model. The displacement and phase fields obtained from the NN at $U_{p}=0.4$
are compared with the FEA solution in Figure 11, showing close agreement. Note
that FEA encounters convergence issues for $U_{p}>0.4$, while the NN does not
exhibit any problem even in that range, as demonstrated by Figure 12, which
shows the fields at $U_{p}=0.47$.

Figure 11: SEN specimen under shear loading ($\omega=0$): displacement and
phase fields at $U_{p}=0.4$ from NN and FEA. Figure 12: SEN specimen under
shear loading ($\omega=0$): displacement and phase fields at $U_{p}=0.47$ from
NN.

The energies associated with the NN solution are compared with the energies
from FEA in Figure 13. They are also very close, although the critical load
for crack propagation in the NN solution is slightly higher than in FEA.
Moreover, higher standard deviations in the energies result from slight
variations in the critical loads for different NN initializations.

Figure 13: SEN specimen under shear loading ($\omega=0$): elastic and fracture
energies vs prescribed displacement from NN and FEA.

###  4.5 Crack branching in a notched plate

Figure 14: Crack branching test: displacement and phase fields at $U_{p}=0.37$
from NN and FEA.

To demonstrate crack branching, we again prescribe the shear loading condition
on the SEN specimen, i.e. $\omega=0$. In this problem, in the absence of the
strain energy decomposition in (4), which differentiates between tensile and
compressive stresses, unphysical branching of the crack occurs [3]. Since in
this case we do not seek a physically realistic prediction of the behavior but
rather an accurate solution of the governing PDEs, we purposefully omit the
strain energy decomposition, i.e. we take $\Psi^{+}=\Psi$ and $\Psi^{-}=0$.
Figure 14 illustrates that the NN is able to learn crack branching. The fields
obtained from the NN at $U_{p}=0.37$ are compared with the FEA solution in
Figure 14, displaying again a good agreement, with some differences in the
prediction of the curvature of the crack path.

The evolution of the energies with $U_{p}$ obtained from the NN are compared
with the FEA energies in Figure 15. The critical load at which branching
initiates shows variability for different NN initializations, leading to a
higher standard deviation; however, the average value is close to the critical
load predicted by FEA. Also, the difficulty encountered by the NN in learning
the curvature of the crack path is reflected in the deviation of the energies
from those computed from FEA after the critical load. Note that one of the 8
NNs had difficulty in learning branching; the corresponding result is not
included in the computation of the mean and standard deviation of the
energies.

Figure 15: Crack branching test: elastic and fracture energies vs prescribed
displacement from NN and FEA.

###  4.6 Crack coalescence in a plate with cracks

To demonstrate crack coalescence, tensile loading is applied on a specimen
with preexisting cracks (Figure 16). In this case, we again adopt (29) with
loading angle $\omega=\pi/2$, and the material properties are $\nu=1/3$ and
$l=0.01$ [64]. We use a larger value for $l$ in comparison to [64] to reduce
the computational cost.

Figure 16: Specimen with preexisting cracks: geometry and boundary conditions
(left) and phase field at the location of the cracks from NN at the smallest
applied load in training (right).

Figure 17 demonstrates that the NN is able to learn crack coalescence. Unlike
in the FEA solution, cracks starting from the initial cracks and reaching the
domain boundaries are at an angle from the horizontal direction. On the other
hand, as discussed earlier, the FEA solution itself is not guaranteed to be
unique [63].

Figure 17: Specimen with preexisting cracks: displacement and phase fields at
$U_{p}=0.2$ from NN and FEA.

The energies obtained from the NN with increasing $U_{p}$ are compared with
the FEA energies in Figure 18. The critical load at which propagation and
merging occurs, while close to the FEA prediction, shows slight variability
for different NN initializations, leading to a higher standard deviation near
the critical load. The energies obtained from the NN after the crack
propagates through the sample are in close agreement with the FEA energies
even though the crack paths show differences. This suggests that different
crack paths may be competing at similar energy levels, which is a known
phenomenon in phase-field fracture modeling [63]. An alternative explanation
may be that the NN is not reaching a sufficiently deep energy minimum, in
which case results can be improved by prescribing a stricter stopping
criterion for the optimizer. Note that here two NNs (out of 8) exhibit
difficulty in learning the correct crack evolution, thus the corresponding
results are not included in the computation of the mean and standard deviation
of the energies.

Figure 18: Specimen with preexisting cracks: elastic and fracture energies vs
prescribed displacement from NN and FEA.

###  4.7 Computational cost

Training of the NNs is performed on NVIDIA GeForce RTX 2080 Ti GPUs, whereas
FEA computations are performed on Intel(R) Xeon(R) W-2223 CPUs. While our DRM
approach is robust and we are able to learn solutions involving various
fracture phenomena, we observe that the DRM is computationally more expensive
compared to FEA, with the training for NNs taking up to one order of magnitude
longer than the corresponding FEA computations. This observation is not
suprising; our aim in this work was not to deliver an alternative solution
framework for a single boundary value problem in phase-field fracture, but
rather to develop a robust approach which can lay the foundation to solve
parametric phase-field fracture problems with the DRM.

##  5 Conclusions

To harness physics-informed deep learning for phase-field fracture modeling,
we propose the design of an NN and of a learning approach based on the DRM,
aimed at learning complex fracture processes. The main conclusions of this
work are as follows:

  * 1.

Phase-field fracture modeling requires solving a coupled problem which
involves fields with sharp spatial variations. Training NNs to learn such
fields is challenging; successful training relies, on one hand, on the ability
of the NN to accurately approximate the energy surface and, on the other hand,
on the ability of an optimizer to reach the correct local minimum respecting
the energy barriers.

  * 2.

We empirically establish that fully connected feed-forward NNs, with a
modified ReLU activation and an appropriate ansatz in the construction of the
solution fields, along with weight regularization and numerical gradient
computation, is able to learn complex fracture processes. Furthermore, RPROP
is found to be the most effective optimization algorithm for training.

  * 3.

Our approach can solve examples of crack initiation, propagation, kinking,
branching, and coalescence taken from the phase-field fracture literature,
within one single numerical setup. Results in terms of solution fields and
energy vs. applied displacement curves are in quantitatively excellent
agreement with the respective FEA results. Additionally, our approach is
robust to different network initializations.

In this work, we focused on the design of a robust NN and learning approach to
learn the solution of boundary value problems involving diverse fracture
phenomena in phase-field fracture modeling. While our approach is
computationally expensive compared to FEA, our aim is not to deliver an
alternative solution framework for a single boundary value problem in phase-
field fracture; rather, the present work is intended as the first step in the
direction of learning solutions to parametric phase-field fracture models with
the DRM. In this setting, NNs can be trained on a few realizations of the
parameter space and results can then be inferred online for all other
realizations, leveraging the true potential of the DRM. Learning solutions to
parametric phase-field fracture models will be the goal of our future
research.

## Acknowledgements

We acknowledge funding from the Swiss National Science Foundation through
Grant No. 200021-219407 ‘Phase-field modeling of fracture and fatigue: from
rigorous theory to fast predictive simulations’.

## Code availability

The code used to produce the numerical results in the paper will be made
available at the time of publication.

## References

  * [1] B. Bourdin, G. A. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture, Journal of the Mechanics and Physics of Solids 48 (4) (2000) 797–826.
  * [2] B. Bourdin, G. A. Francfort, J.-J. Marigo, The variational approach to fracture, Journal of elasticity 91 (1) (2008) 5–148.
  * [3] M. Ambati, T. Gerasimov, L. De Lorenzis, A review on phase-field models of brittle fracture and a new fast hybrid formulation, Computational Mechanics 55 (2) (2015) 383–405.
  * [4] J.-Y. Wu, V. P. Nguyen, C. T. Nguyen, D. Sutula, S. Sinaie, S. P. Bordas, Phase-field modeling of fracture, Advances in applied mechanics 53 (2020) 1–183.
  * [5] B. Bourdin, J.-J. Marigo, C. Maurini, P. Sicsic, Morphogenesis and propagation of complex cracks induced by thermal shocks, Physical review letters 112 (1) (2014) 014301.
  * [6] R. Alessi, M. Ambati, T. Gerasimov, S. Vidoli, L. De Lorenzis, Comparison of phase-field models of fracture coupled with plasticity, Advances in computational plasticity: A book in honour of D. Roger J. Owen (2018) 1–21.
  * [7] C. Luo, L. Sanavia, L. De Lorenzis, Phase-field modeling of drying-induced cracks: Choice of coupling and study of homogeneous and localized damage, Computer Methods in Applied Mechanics and Engineering 410 (2023) 115962.
  * [8] Y. Heider, A review on phase-field modeling of hydraulic fracturing, Engineering Fracture Mechanics 253 (2021) 107881.
  * [9] T. Wu, A. Carpiuc, M. Poncelet, L. De Lorenzis, Phase-field simulation of interactive mixed-mode fracture tests on cement mortar with full-field displacement boundary conditions, Engineering Fracture Mechanics 182 (2017) 658–688.
  * [10] G. A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, Journal of the Mechanics and Physics of Solids 46 (8) (1998) 1319–1342.
  * [11] K. Pham, H. Amor, J.-J. Marigo, C. Maurini, Gradient damage models and their use to approximate brittle fracture, International Journal of Damage Mechanics 20 (4) (2011) 618–652.
  * [12] L. D. Lorenzis, T. Gerasimov, Numerical implementation of phase-field models of brittle fracture, in: Modeling in Engineering Using Innovative Numerical Methods for Solids and Fluids, Springer, 2020, pp. 75–101.
  * [13] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nature Reviews Physics 3 (6) (2021) 422–440.
  * [14] I. E. Lagaris, A. Likas, D. I. Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE transactions on neural networks 9 (5) (1998) 987–1000.
  * [15] J. Sirignano, K. Spiliopoulos, Dgm: A deep learning algorithm for solving partial differential equations, Journal of computational physics 375 (2018) 1339–1364.
  * [16] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics 378 (2019) 686–707.
  * [17] E. Kharazmi, Z. Zhang, G. E. Karniadakis, hp-vpinns: Variational physics-informed neural networks with domain decomposition, Computer Methods in Applied Mechanics and Engineering 374 (2021) 113547.
  * [18] B. Yu, et al., The deep ritz method: a deep learning-based numerical algorithm for solving variational problems, Communications in Mathematics and Statistics 6 (1) (2018) 1–12.
  * [19] S. Mishra, A machine learning framework for data driven acceleration of computations of differential equations, Mathematics in Engineering 1 (1) (2019) 118–146.
  * [20] S. Mishra, R. Molinaro, Estimates on the generalization error of physics-informed neural networks for approximating pdes, IMA Journal of Numerical Analysis (2022).
  * [21] G. Gybenko, et al., Approximation by superposition of sigmoidal functions, Mathematics of Control, Signals and Systems 2 (4) (1989) 303–314.
  * [22] K. Hornik, M. Stinchcombe, H. White, Multilayer feedforward networks are universal approximators, Neural networks 2 (5) (1989) 359–366.
  * [23] D. Yarotsky, Error bounds for approximations with deep relu networks, Neural Networks 94 (2017) 103–114.
  * [24] E. Haghighat, M. Raissi, A. Moure, H. Gomez, R. Juanes, A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics, Computer Methods in Applied Mechanics and Engineering 379 (2021) 113741.
  * [25] C. Rao, H. Sun, Y. Liu, Physics-informed deep learning for computational elastodynamics without labeled data, Journal of Engineering Mechanics 147 (8) (2021) 04021043.
  * [26] E. Zhang, M. Dao, G. E. Karniadakis, S. Suresh, Analyses of internal structures and defects in materials using physics-informed neural networks, Science advances 8 (7) (2022) eabk0644.
  * [27] S. Rezaei, A. Harandi, A. Moeineddin, B.-X. Xu, S. Reese, A mixed formulation for physics-informed neural networks as a potential solver for engineering problems in heterogeneous domains: comparison with finite element method, Computer Methods in Applied Mechanics and Engineering 401 (2022) 115616.
  * [28] A. Henkes, H. Wessels, R. Mahnken, Physics informed neural networks for continuum micromechanics, Computer Methods in Applied Mechanics and Engineering 393 (2022) 114790.
  * [29] J.-H. Bastek, D. M. Kochmann, Physics-informed neural networks for shell structures, European Journal of Mechanics-A/Solids 97 (2023) 104849.
  * [30] S. Niu, E. Zhang, Y. Bazilevs, V. Srivastava, Modeling finite-strain plasticity using physics-informed neural network and assessment of the network performance, Journal of the Mechanics and Physics of Solids 172 (2023) 105177\.
  * [31] S. A. Faroughi, N. M. Pawar, C. Fernandes, M. Raissi, S. Das, N. K. Kalantari, S. Kourosh Mahjour, Physics-guided, physics-informed, and physics-encoded neural networks and operators in scientific computing: Fluid and solid mechanics, Journal of Computing and Information Science in Engineering 24 (4) (2024).
  * [32] Y. Shin, J. Darbon, G. E. Karniadakis, On the convergence and generalization of physics informed neural networks, arXiv e-prints (2020) arXiv–2004.
  * [33] S. Mishra, R. Molinaro, Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for pdes, IMA Journal of Numerical Analysis 42 (2) (2022) 981–1022.
  * [34] A. D. Jagtap, E. Kharazmi, G. E. Karniadakis, Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems, Computer Methods in Applied Mechanics and Engineering 365 (2020) 113028.
  * [35] A. D. Jagtap, G. E. Karniadakis, Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations., in: AAAI Spring Symposium: MLPS, 2021, pp. 2002–2041.
  * [36] B. Moseley, A. Markham, T. Nissen-Meyer, Finite basis physics-informed neural networks (fbpinns): a scalable domain decomposition approach for solving differential equations, Advances in Computational Mathematics 49 (4) (2023) 62\.
  * [37] V. M. Nguyen-Thanh, X. Zhuang, T. Rabczuk, A deep energy method for finite deformation hyperelasticity, European Journal of Mechanics-A/Solids 80 (2020) 103874\.
  * [38] E. Samaniego, C. Anitescu, S. Goswami, V. M. Nguyen-Thanh, H. Guo, K. Hamdia, X. Zhuang, T. Rabczuk, An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications, Computer Methods in Applied Mechanics and Engineering 362 (2020) 112790.
  * [39] S. Berrone, C. Canuto, M. Pintore, Variational physics informed neural networks: the role of quadratures and test functions, Journal of Scientific Computing 92 (3) (2022) 100.
  * [40] C. L. Wight, J. Zhao, Solving allen-cahn and cahn-hilliard equations using the adaptive physics informed neural networks, arXiv preprint arXiv:2007.04542 (2020).
  * [41] R. Mattey, S. Ghosh, A novel sequential method to train physics informed neural networks for allen cahn and cahn hilliard equations, Computer Methods in Applied Mechanics and Engineering 390 (2022) 114474.
  * [42] S. Goswami, C. Anitescu, S. Chakraborty, T. Rabczuk, Transfer learning enhanced physics informed neural network for phase-field modeling of fracture, Theoretical and Applied Fracture Mechanics 106 (2020) 102447.
  * [43] S. Goswami, C. Anitescu, T. Rabczuk, Adaptive fourth-order phase field analysis using deep energy minimization, Theoretical and Applied Fracture Mechanics 107 (2020) 102527.
  * [44] Y. Ghaffari Motlagh, P. K. Jimack, R. de Borst, Deep learning phase-field model for brittle fractures, International Journal for Numerical Methods in Engineering (2022).
  * [45] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis, Learning nonlinear operators via deeponet based on the universal approximation theorem of operators, Nature machine intelligence 3 (3) (2021) 218–229.
  * [46] Z. Li, H. Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, A. Anandkumar, Physics-informed neural operator for learning partial differential equations, arXiv preprint arXiv:2111.03794 (2021).
  * [47] S. Goswami, M. Yin, Y. Yu, G. E. Karniadakis, A physics-informed variational deeponet for predicting crack path in quasi-brittle materials, Computer Methods in Applied Mechanics and Engineering 391 (2022) 114587.
  * [48] H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments, Journal of the Mechanics and Physics of Solids 57 (8) (2009) 1209–1229.
  * [49] F. Freddi, G. Royer-Carfagni, Regularized variational theories of fracture: a unified approach, Journal of the Mechanics and Physics of Solids 58 (8) (2010) 1154–1174.
  * [50] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations, International journal for numerical methods in engineering 83 (10) (2010) 1273–1311.
  * [51] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Computer Methods in Applied Mechanics and Engineering 199 (45-48) (2010) 2765–2778.
  * [52] F. Vicentini, C. Zolesi, P. Carrara, C. Maurini, L. De Lorenzis, On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states, Engineering Fracture Mechanics 182 (2017) 658–688.
  * [53] T. Gerasimov, L. De Lorenzis, On penalization in variational phase-field models of brittle fracture, Computer Methods in Applied Mechanics and Engineering 354 (2019) 990–1026.
  * [54] S. R. Dubey, S. K. Singh, B. B. Chaudhuri, Activation functions in deep learning: A comprehensive survey and benchmark, Neurocomputing 503 (2022) 92–108.
  * [55] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * [56] D. C. Liu, J. Nocedal, On the limited memory bfgs method for large scale optimization, Mathematical programming 45 (1) (1989) 503–528.
  * [57] J. Berg, K. Nyström, A unified deep artificial neural network approach to partial differential equations in complex geometries, Neurocomputing 317 (2018) 28–41.
  * [58] C. Leake, D. Mortari, Deep theory of functional connections: A new method for estimating the solutions of partial differential equations, Machine learning and knowledge extraction 2 (1) (2020) 37–55.
  * [59] M. Riedmiller, H. Braun, A direct adaptive method for faster backpropagation learning: The rprop algorithm, in: IEEE international conference on neural networks, IEEE, 1993, pp. 586–591.
  * [60] X. Glorot, Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in: Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, 2010, pp. 249–256.
  * [61] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala, Pytorch: An imperative style, high-performance deep learning library, in: H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, R. Garnett (Eds.), Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 8024–8035.
  * [62] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities, International journal for numerical methods in engineering 79 (11) (2009) 1309–1331.
  * [63] T. Gerasimov, U. Römer, J. Vondrejc, H. Matthies, L. De Lorenzis, Stochastic phase-field modeling of brittle fracture: computing multiple crack patterns and their probabilities, Computer Methods in Applied Mechanics and Engineering 372 (2020) 113353.
  * [64] S. Zhou, X. Zhuang, H. Zhu, T. Rabczuk, Phase field modelling of crack propagation, branching and coalescence in rocks, Theoretical and Applied Fracture Mechanics 96 (2018) 174–192.

##  Appendix A Independent NNs for the displacement and phase fields

Using independent NNs for the physically different fields to learn their
unique features separately can be expected to make training faster. We
construct two NNs of 8 hidden layers and 200 neurons, representing the
displacement and phase fields separately, with all the other settings as in
Section 4. Note that this NN configuration lowers the number of the learnable
parameters. We first study the SEN specimen under shear loading described in
Section 4.4.2, and conclude that the NNs are able to learn crack propagation,
as shown in Figure 19. Subsequently, we use the two NNs to learn crack
nucleation in an L-shaped panel as described in Section 4.3; here, the NNs are
unsuccessful, as shown in Figure 20. So, we use one NN to represent both the
fields for all the problems.

Figure 19: Phase field at $U_{p}=0.4$ after crack propagation (left) and
energies (right) for the SEN specimen under shear loading ($\omega=0$) when
the displacement and phase fields are represented by two separate NNs.

Figure 20: Phase field at $U_{p}=1.2$ (left) and energies (right) for the
L-shaped panel when the displacement and phase fields are represented by two
separate NNs. In this case, the crack fails to nucleate.

##  Appendix B Activation functions

We assess the ability of the NNs with different activation functions in
learning crack propagation in the SEN specimen under shear loading. We use the
same NN size and optimizer settings as in Section 4.4.2, but set the
activation function to be $\tanh$. As shown in Figure 21, the NN is able to
learn the displacement and phase fields which yield sharp stress fields near
the crack tip.

Figure 21: Stress field at $U_{p}=0.31$ in the SEN specimen under shear
loading ($\omega=0$).

It also learns the critical $U_{p}$ for crack propagation accurately. However,
it fails to learn the correct crack path as shown in Figure 22. Furthermore,
the smoothly decreasing loss and the comparison of the phase fields at two
different optimization steps shown in Figure 22 suggest that the NN-based
representation of the loss lacks the energy barrier which could prevent the NN
from learning the incorrect solution.

Figure 22: Phase field at $U_{p}=0.32$ in the SEN specimen under shear loading
($\omega=0$) (left) and evolution of the loss function during training
(right).

##  Appendix C Autodifferentiation for gradient computation

When employing a smooth activation and autodifferentiation for gradient
computation, the solution field exhibits spurious oscillation. The
displacement field in the 1D bar problem in Section 4.2, learned using an NN
with 2 hidden layers and 20 neurons with the $\tanh$ activation function,
exhibits oscillations as shown in Figure 23.

Figure 23: Spurious oscillations observed in the displacement field after the
localization of the phase field in the 1D bar problem.

We also train an NN with 6 hidden layers and 50 neurons with $\tanh$
activation to learn the solution of the SEN specimen under shear loading.
While spurious oscillations in the displacement field are not apparent, they
can be observed in the stress field as shown in Figure 24.

Figure 24: Displacement field, phase field, and stress field at $U_{p}=0.31$
in the SEN specimen under shear loading ($\omega=0$).

Furthermore, at the critical $U_{p}$, the crack propagates in an incorrect
direction as shown in Figure 25.

Figure 25: Phase field at $U_{p}=0.32$ in the SEN specimen under shear loading
($\omega=0$).

##  Appendix D The maximum strain in the 1D bar

For the bar described in Section 4.2, the equilibrium equation is given by

$\frac{{\rm d}\sigma}{{\rm d}x}=0,$ (30)

hence, the stress in the bar is uniform. Since $\sigma=g(\alpha)\frac{{\rm
d}u}{{\rm d}x}$, the strain can be expressed as $\frac{{\rm d}u}{{\rm
d}x}=\frac{\sigma}{g(\alpha)}=S(\alpha)\sigma$, where $S(\alpha)=1/g(\alpha)$
is the compliance.

During evolution of the phase field, eq. (11) (written as equality since
$\alpha>\alpha_{n-1}$), reads

$-\frac{1}{2}S^{\prime}(\alpha)\sigma^{2}+w_{1}\left[w^{\prime}(\alpha)-2l^{2}\frac{{\rm
d}^{2}\alpha}{{\rm d}x^{2}}\right]=0,$ (31)

where $w_{1}=\frac{1}{c_{w}}$ and the first term has been modified to express
it in terms of compliance and stress. Multiplying the above equation with
$\frac{{\rm d}\alpha}{{\rm d}x}$ and integrating with respect to $x$ yields

$-\frac{1}{2}S(\alpha)\sigma^{2}+w_{1}\left[w(\alpha)-l^{2}\left(\frac{{\rm
d}\alpha}{{\rm d}x}\right)^{2}\right]=C,$ (32)

with $C$ as the integration constant. Since, far away from the point of
localization of the phase field, $\alpha$ and its gradient vanish (thus $w=0$
and $S=1$), we obtain $C=-\frac{1}{2}\sigma^{2}$; further substitution of
$\frac{{\rm d}u}{{\rm d}x}=S(\alpha)\sigma$ yields

$\frac{{\rm d}u}{{\rm
d}x}=\sigma+\frac{2w_{1}}{\sigma}\left[w(\alpha)-l^{2}\left(\frac{{\rm
d}\alpha}{{\rm d}x}\right)^{2}\right].$ (33)

At the location in the bar where the maximum strain occurs, the following
condition holds

$\frac{{\rm d}}{{\rm d}x}\left(\frac{{\rm d}u}{{\rm
d}x}\right)=\left(w^{\prime}(\alpha)-2l^{2}\frac{{\rm d^{2}}\alpha}{{\rm
d}x^{2}}\right)\frac{{\rm d}\alpha}{{\rm d}x}=0,$ (34)

hence, $\frac{{\rm d}\alpha}{{\rm d}x}=0$ and

$\left(\frac{{\rm d}u}{{\rm
d}x}\right)_{max}=\sigma+\frac{2w_{1}w(\alpha)}{\sigma}.$ (35)

When the phase field localizes, the stress in the bar goes from a positive
value to 0 (however, for $\sigma=0$ the so-called optimal profile forms, where
$\frac{{\rm d}\alpha}{{\rm d}x}\neq 0$ at the center). When the stress
approaches 0, the maximum of strain, given by (35), is unbounded. In FEA, the
maximum strain is set by the element size which is chosen small enough to
resolve the length scale $l$. In contrast, the maximum strain in the DRM is
dictated by the expressivity of the NN. Consequently, if the training points
are not fine enough to resolve the sharp slope of the displacement field
expressible by an NN and not only the length scale $l$, the DRM yields an
incorrect solution as shown in Figure 26.

Figure 26: Incorrect displacement field, phase field and energies when the
maximum strain expressed by the NN solution is not controlled. The plots are
obtained using an NN with 2 hidden layers and 20 neurons with the $\tanh$
activation. The correct solution fields and energies are shown in Figure 4 and
Figure 3, respectively.

##  Appendix E The limitations of the L-BFGS optimizer

To probe the ability of the L-BFGS optimizer in learning the sharp fields in
the DRM, we apply it to learn the solution of the SEN sample under shear
loading as described in Section 4.4.2. We train the same NN used in Section
4.4.2. As shown in Figure 27, the NN fails to learn the sharp displacement
field even before crack propagation begins. It also does not learn crack
propagation as $U_{p}$ is increased.

Figure 27: Displacement and phase fields at $U_{p}=0.31$ for the SEN specimen
under shear loading ($\omega=0$) when the fields are learned using only the
L-BFGS optimizer.

##  Appendix F Network size estimation

We conduct numerical experiments with different network sizes for the problem
of crack nucleation in an L-shaped panel in Section 4.3. The energies
associated with the solutions obtained from networks with different sizes are
compared in Figure 28. Based on these experiments, we choose a network with 8
hidden layers and 400 neurons for all the 2D problems solved in this work.

Figure 28: Energy associated with the solutions learned with different network
sizes vs prescribed displacement. In the left plot, number of neurons is 400;
in the right plot, number of hidden layers is 8. (hl: hidden layers)

##  Appendix G Numerical setup for FEA

In the FEA computation for the 1D problem in Section 4.2, we discretize the
domain into elements of size $l/5$ and assume linear shape functions. To
obtain the solution of the 2D problems, we employ quadrilateral elements with
bilinear shape functions and two Gauss points per parametric direction. The
regions of the domain where the crack is expected to propagate are discretized
with elements of size $l/5$. The element size smoothly increases up to $4l$
away from these regions.

Initial cracks are modeled by prescribing $\alpha=1$ at the location of the
crack and solving a so-called recovery problem. We set an irreversibility
tolerance of $\texttt{TOL}_{ir}=10^{-3}$. The staggered solution algorithm
together with the Newton-Raphson procedure is employed to iteratively converge
to the solution [53]. We set an error tolerance on the residual norm of
$10^{-4}$ for the staggered scheme and $10^{-6}$ for the Newton-Raphson
procedure. Moreover, we set a maximum of 500 iterations for the Newton-Raphson
procedure and 1000 iterations for the staggered scheme.

##  Appendix H Distance functions for the L-shaped panel

To ensure the strict enforcement of the Dirichlet boundary conditions in the
problem in Section 4.3, we construct distance functions that satisfy $C^{1}$
continuity. For applying the fixed boundary condition at $y=-0.5$, we
construct the following distance function:

$d_{1}(x,y)=\begin{cases}\left(\frac{y+0.5}{a}-1\right)^{2},&{\rm
for~{}}y+0.5\leq a\\\ 0,&{\rm otherwise}\end{cases}$ (36)

To apply the prescribed displacement on a section of the boundary shown in
Figure 5, we construct the following distance function:

$d_{2}(x,y)=\begin{cases}\left(\frac{|y|}{a}-1\right)^{2},&{\rm for~{}}x\geq
0.44{\rm~{}and~{}}|y|\leq a\\\
\left(\frac{\sqrt{(0.44-x)^{2}+y^{2}}}{a}-1\right)^{2},&{\rm for~{}}x\leq
0.44{\rm~{}and~{}}\sqrt{(0.44-x)^{2}+y^{2}}\leq a\\\ 0,&{\rm
otherwise}\end{cases}$ (37)

In this work, we set $a=0.1$.

##  Appendix I Initial $\alpha$ in a notched sample

To define the initial $\alpha$ in a notched sample, we make use of the
analytical expression for $\alpha$ in a 1D bar upon crack nucleation [4].
Below we provide the expression for the initial $\alpha$ in a notched sample
shown in Figure 29 in the AT1 model. Note that $x^{\prime}-y^{\prime}$ is the
transformed coordinate system aligned with the crack and with its origin at
one end of the crack.

Figure 29: A sample with an initial crack

$\alpha_{0}=\begin{cases}\left(1-\frac{|y^{\prime}|}{2l}\right)^{2},&{\rm
for~{}}0<x^{\prime}<L_{c}{\rm~{}and~{}}|y^{\prime}|<2l\\\
\left(1-\frac{\sqrt{x^{\prime 2}+y^{\prime 2}}}{2l}\right)^{2},&{\rm
for~{}}x^{\prime}<0{\rm~{}and~{}}\sqrt{x^{\prime 2}+y^{\prime 2}}<2l\\\
\left(1-\frac{\sqrt{(x^{\prime}-L_{c})^{2}+y^{\prime 2}}}{2l}\right)^{2},&{\rm
for~{}}x^{\prime}>L_{c}{\rm~{}and~{}}\sqrt{(x^{\prime}-L_{c})^{2}+y^{\prime
2}}<2l\\\ 0,&{\rm otherwise}\end{cases}$ (38)