\section{Introduction}

Over the years, brittle fracture simulation has received more attention in many fields of the sciences. The overarching goal of brittle fracture modeling is to predict the onset of brittle fracture and its characteristics, such as crack trajectory, crack velocity, crack branching, and coalescence~\cite{MIEHE_2010, BLEYER_2018, TANNE_2018, ULLOA_2019}.
The classical Griffith theory of brittle fracture provides a criterion for crack propagation, yet it is alone incapable of determining crack trajectories or crack branching~\cite{MIEHE_2010, CHARLOTTE_2006}, this motivates the development of numerous new methods and methodologies. An isotropic second-order phase-field fracture model was introduced in \cite{FRANCFORT_1998}, an anisotropic model derivated from the model in \cite{BOURDIN_2000} is studied in \cite{AMOR_2009}, the model considers the elastic energy density decomposed into volumetric and deviatoric components to prevent cracking in the compression domain. A fourth-order phase field model to improve convergence rates for numerical solutions using isogeometric finite elements is analyzed in \cite{BORDEN_2014}. In order to reduce computational time in simulations for the hybrid phase-field model, an optimization method of matrix assembly for the FEM method with refined mesh is proposed in \cite{Yosra_2020}.  


In recent research, methods based on variational approaches have been studied \cite{FRANCFORT_1998,BOURDIN_2008,MASO_2002,BULIGA_1998}, within this variational framework, the process of brittle fracture is fully characterized by directly minimizing an energy functional \cite{BOURDIN_2008,MESGARNEJAD_2015,AREIAS_2016}, and the numerical methods can be systematically derived through $\Gamma$-convergence analysis~\cite{MUNFORD_1989}.


\begin{table*}[ht!]
\indexspace
\begin{center}
\caption{List of principal parameters.
}\label{table_vars}
\begin{tabular}{|ll|}
\hline 
$\boldsymbol{u}$ & Displacement field\\
$\phi$ & Scalar damage variable\\
$\lambda$, $\mu$ & Lamm\'e constants (kN/mm$^2$)\\
$G_c$ & Critical energy release rate (kN/mm)\\
$\eta_n$ & Constant of order $O(\ell)$ \\
$\ell$ & Regularization length (mm)\\
$\psi$ & Elastic energy density\\
$\boldsymbol{\varepsilon}$ & Small strain tensor\\
$H$ & History variables \\
$H_n$ & Strain energy computed at load step n\\
$\Delta u$ & Incremental displacement (mm)\\
$\overline{\boldsymbol{t}}$ & Traction (kN/mm$^2$) \\
$C(\cdot)$ & Multidimensional covariance function \\
$S(\cdot)$ & Spectral desity function \\
$\mathcal{F}$ & Fourier transform \\
\hline
\end{tabular}
\end{center}
\end{table*}

In the last decade, many phase-field brittle fracture models have been introduced, yielding remarkably accurate predictions to simulate brittle fractures for various material classes and loading scenarios~\cite{MIEHE_2010,BORDEN_2012, HESCH_2014, WEINBERG_2017}, these models are derivated on the variational formulation~\cite{MIEHE_2010}. The popularity of the phase-field method stems from the fact that conventional numerical methods, such as finite elements, can be easily employed to seek solutions of phase-field equations, due that the cracks are represented only by a continuous field as is elucidated in~\cite{HIRSHIKESH_2019}. 

Models of brittle fracture often assume spatial homogeneity of the underlying material properties. This assumption is, however, seldom valid as the microstructure underlying material behavior tends to be endowed with randomness due to the presence of defects, such as grains, inclusions, voids, etc., hence, for modeling brittle fractures, non-homogeneous properties must be included. A few brittle fracture models directly incorporating randomness have been proposed, for instance, a heterogeneous cohesive (HC) crack model is developed to predict macroscopic strength of materials based on meso-scale random fields of fracture properties in \cite{YANG_2008}. Moreover, a study of potential cracks represented by pre-inserted cohesive elements with tension and shear softening constitutive laws are modeled by spatially-varying Weibull random fields in \cite{YANG_2009, SU_2010}.
 

In this article, we propose to employ a Gaussian random field with a Mat\'ern covariance matrix to model spatial non-homogeneity of the energy release rate $G_c$. 
We also propose to use a weighted-variational approach as a surrogate model for the hybrid phase-field model to expedite evaluations of the phase-field brittle fracture when we consider the energy release rate $G_C$ as a Gaussian random field.  
Our proposal of the weighted-variational model, for the selected examples, provide reasonable surrogate simulations to the hybrid phase-field model, reducing execution times by 90\%.  To solve the equations of the variational and phase-field models, we implement a code in Python with the open-source finite element software FEniCS,  following previous researches in ~\cite{HIRSHIKESH_2019} and~\cite{TANNE_2018}.



The paper is organized as follows: Section \ref{sec:2} presents an overview of the variational and phase-field mathematical models. In section \ref{sec:3}, the formulation of the Gaussian random field of the energy release rate with Mat\'ern covariance matrix is presented. Details of the implementation in FEniCS of the
mathematical approaches are presented in Section \ref{sec:4}. In Section \ref{sec:5}, we present the solution of the standard boundary value problems to validate the implementation in FEniCS, we show simulations using the Gaussian random fields for $G_c$, and we present an application of our proposal of surrogate model in a laboratory experiment. Finally, a conclusion of this work is presented in Section \ref{sec:6}.




\section{Mathematical models}\label{sec:2}

In this section, we present two mathematical models to describe brittle fracture propagations: the variational approach and the phase field approach.

\subsection{Variational approach for brittle fracture}

Variational phase-field models of brittle fracture have their origins in the research of Bourdin {\it et al.}~\cite{BOURDIN_2000}, these approaches are based on a regularization of the Mumford-Shah problem in image processing~\cite{AMBROSIO_1990, AMBROSIO_1992,BRAIDES_1998}. The variational model considers a linear-elastic body with a crack in the body $\Omega \subset \mathbb{R}^n$ , where $n=2,3$.  Let be $\boldsymbol{u}: \Omega \rightarrow \mathbb{R}^n$ the displacement field, and $\phi: \Omega \rightarrow [0, 1]$ the scalar damage variable, describes the trajectory of the crack on $\Omega$. The body is subjected to boundary conditions: tractions, $\bar{\tau}$, applied along the neumann boundary $\Gamma_{N}$; displacements, $\overline{\boldsymbol{u}}$, applied along the dirichlet boundary $\Gamma_{D}$, see Figure \ref{body_1}. The boundary value
problem for brittle fracture is: find $(\boldsymbol{u}, \phi): \Omega
\rightarrow \mathbb{R}^d\times[0,1]$ satisfying the minimization
problem,
\begin{equation}
  \label{eq:minimization}
 \min E_{\ell}(\boldsymbol{u}, \phi)
\end{equation}
where,
\begin{equation}
 E_{\ell}(\boldsymbol{u}, \phi)=\int_{\Omega} \left[a(\phi)+\eta_{\ell}\right]\psi(\boldsymbol{\varepsilon}, \phi ) d \Omega+\frac{G_{c}}{4 c_{w}} \int_{\Omega}\left(\frac{w(\phi)}{\ell}+\ell|\nabla \phi|^{2}\right) d\Omega,
\end{equation}
the symmetrized gradient of $\boldsymbol{u}$ defined as $\boldsymbol{\varepsilon}(\boldsymbol{u}):=\frac{1}{2}\left(\nabla \boldsymbol{u}+\nabla^{T}
\boldsymbol{u}\right)$  denotes the strain tensor; the elastic energy density function described in terms of the
elastic Lamm\'e constants is $\psi(\boldsymbol{\varepsilon}, \phi ) = \frac{1}{2}A\varepsilon(\boldsymbol{u}):\varepsilon(\boldsymbol{u}) = \frac{1}{2} \lambda tr^2(\varepsilon(\boldsymbol{u})) + \mu tr(\varepsilon^2(\boldsymbol{u}))$, with $\lambda>0$ and $\mu>0$; $a(\phi)$ and $w(\phi)$ are
continuous monotonic functions such that $a(0)=1$, $a(1)=0$, $w(0)=0,$ and
$w(1)=1$,  $\eta_{\ell}=o(\ell)$ is a small parameter, and
$c_{w}:=\int_{0}^{1} \sqrt{w(s)} d s$ is a normalization term. { The regularization length $\ell>0$ is the length scale associated with the phase-field regularization of the fracture surface, and its value is important and necessary for numerical simulations, that is, to capture the propagation of the crack, the element size for finite element method should be smaller than $\ell$. }

 \begin{figure}[H]
    \begin{center} 
  \includegraphics[width=7 cm]{crack_01.png}
\end{center} 
\caption{Specimen $\Omega$ with an inner crack $\Gamma_C$: a) sharp crack, b) approximated diffuse crack as function of $\ell$.}\label{body_1}
\end{figure}


The crack is described by the smooth transition function $\phi$, and must satisfy the following conditions ~\cite{MSEKH_2015, AMBATI_2015, MIEHE_2010a}:

\begin{itemize}
 \item is a symetric function respect to the crack,
 \item is a monotonic increasing function on time,
 \item the maximum value of $\phi$ on the crack is one, and decays to zero as it moves away from the fracture.
\end{itemize}

The most common variational mathematical model in the literature
\cite{TANNE_2018} employs $a(\phi) = (1-\phi)^2$ and $w(\phi) =
\phi^2$, and the functional becomes to,

\begin{equation}\label{eq21}
 E_{\ell}(\boldsymbol{u}, \phi)=\int_{\Omega} \left[(1-\phi)^{2}+\eta_{\ell} \right ]\psi(\boldsymbol{\varepsilon}, \phi ) d \Omega +\frac{G_{c}}{2} \int_{\Omega} \left( \frac{\phi^{2}}{\ell}+\ell|\nabla \phi|^{2} \right) d \Omega .
\end{equation}


\subsection{Phase field method}

In recent years, the phase-field mathematical model for brittle fracture has been used to simulate brittle cracks \cite{AMBATI_2015, HIRSHIKESH_2019}.
This approach has its foundations in the variational approach, due the total energy \eqref{eq:minimization} is convex with respect to $\mathbf{u}$ and $\phi$ separately, but not with respect to both of them \cite{MUIXI_2020}, this motivates the solution of the system by means of a staggered scheme of their principal functionals. The mathematical formulation is obtained from the first variation of \eqref{eq21} (see details in \cite{GULTEKIN_2016}). The resulting system consists of two coupled non-linear differential equations, the
first equation describes the phase-field to emulate the fracture, and the second equation describes the displacements on the domain $\Omega$. The phase field system is described by
\cite{HIRSHIKESH_2019},

\begin{equation}\label{sis_1}
\left\lbrace\begin{array}{ll} 
 \left[ (1 - \phi)^2 + \eta_{\ell} \right] \nabla \cdot \boldsymbol{\sigma} = 0, & \text{in } \quad \Omega, \\
 -G_c \ell \nabla ^2 \phi + \left[ \frac{G_c}{\ell} + 2 H \right] \phi = 2H, & \text{in } \quad \Omega,
\end{array}\right.
\end{equation}

with boundary conditions:

\begin{equation}
\left\lbrace\begin{array}{ll} 
 \left[ (1 - \phi)^2 + \eta_{\ell} \right] \nabla \cdot \boldsymbol{\sigma} = \overline{\boldsymbol{t}}, & \text{on } \quad \Gamma_N, \\
 \boldsymbol{u} = \overline{\boldsymbol{u}}, & \text{on } \quad \Gamma_D, \\
 \nabla \phi \cdot \boldsymbol{n} = 0, & \text{on } \quad \Gamma_N,
\end{array}\right.
\end{equation}

where $\boldsymbol{\sigma} = \frac{\partial \psi (\boldsymbol{\varepsilon})}{\partial \boldsymbol{\varepsilon}}$ is the Cauchy stress tensor. 
The isotropic model in \cite{AMBATI_2015} depends on the history variable $H$ defined as:
\begin{equation}
H = 
\begin{cases}
\psi(\varepsilon) & \psi(\varepsilon) < H_n \\
H_n & \text{otherwise},
\end{cases}
\end{equation}
where $H_n$ is the strain energy computed at load step n. 

A modification of the phase-field method known as hybrid phase field method (Hybrid P-F model) is analized in \cite{AMBATI_2015}, this model
is given by,
\begin{equation}
  \label{ec3}
 -G_c l_0 \nabla ^2 \phi + \left[ \frac{G_c}{l_o} + 2 H^+ \right] \phi = 2H^+, \quad \text{in } \quad \Omega,
\end{equation}
where $H^+ = max_{\tau \in [0,t]}\psi^+ \left( \varepsilon (x, \tau) \right)$.

To prevent crack face interpenetration, the model should satisfy the
constraint,

\begin{equation}
 \forall x: \psi^+ < \psi^- \Rightarrow \phi:=0,
\end{equation}
and 
\begin{equation}
 \psi^{\pm}(\varepsilon) = \frac12 \lambda \langle tr(\varepsilon) \rangle^2_{\pm} + \mu tr(\varepsilon^2_{\pm})
\end{equation}
with 
\begin{equation}
\varepsilon_{\pm} = \sum^3_{I = 1} \langle \varepsilon_I \rangle_{\pm}{\bf{n_I}}\otimes {\bf{n_I}}, \quad \text{and} \quad \varepsilon =  \sum^3_{I = 1} \langle \varepsilon_I \rangle \bf{n_I}\otimes \bf{n_I} ,
\end{equation}
where $\langle x \rangle_{+}:=(|x|+x)/2$ and $\langle x \rangle_{-}:=(|x|-x) / 2$, and $\varepsilon_I$ and $\bf{n_I}$ are the principle strains and principle strain directions, respectively.  

This model has the property that the linear momentum equation is retained as in the isotropic model and the evolution of the phase-field variable is controlled by the tensile elastic energy $\psi^+$, ensuring that the crack surfaces do not grow in the regions of compression, avoiding imaginary cracks. The hybrid formulation reduce the computational complexity reducing execution times for
simulations \cite{HIRSHIKESH_2019}.


\section{Random energy release rate $G_c$}\label{sec:3}

The behavior of materials under stress is governed primarily by their microstructure. For many materials, the microstructure can rarely be characterized purely from a deterministic viewpoint as the size and spatial distribution of defects constituting the microstructure tends to be random, such as glass, ceramic, concrete, and some alloys. Hence, many material properties of brittle materials, such as fracture strength or fracture energy, have to be modeled as random. However, material properties are very often spatially dependent, and thus, can be effectively handled employing random fields~\cite{YANG_2008, YANG_2009, Jeulin_1994}. In the context of brittle fracture, it is clear that a fracture path will be ultimately determined by the laws of physics, yet it will be strongly influenced by the inherent material randomness. 

We propose to model the randomness associated with a brittle fracture through Gaussian random fields, specifically, we will assume that the energy release rate $G_c$ is a stationary Gaussian random field. 

A Gaussian random field is defined as a random function $f:\Omega \rightarrow R$ such that any finite collection of random variables $f(\boldsymbol{x}_1),\ldots,f(\boldsymbol{x}_n)$ at points $\boldsymbol{x}_1,\ldots,\boldsymbol{x}_n \in \Omega$ has a
multidimensional Gaussian distribution \cite{YANG_2008, Sarkka_2013}. A Gaussian field is fully characterized in terms of its mean $m(\boldsymbol{x})$ and covariance function $k\left(\boldsymbol{x}, \boldsymbol{x}^{\prime}\right),$
\[
\begin{aligned}
m(\boldsymbol{x}) &=\mathrm{E}[f(\boldsymbol{x})] \\
k\left(\boldsymbol{x}, \boldsymbol{x}^{\prime}\right) &=\mathrm{E}\left[(f(\boldsymbol{x})-m(\boldsymbol{x}))\left(f\left(\boldsymbol{x}^{\prime}\right)-m\left(\boldsymbol{x}^{\prime}\right)\right)\right] .
\end{aligned}
\]
The joint distribution of an arbitrary finite collection of random variables $f\left(\boldsymbol{x}_{1}\right), \ldots, f\left(\boldsymbol{x}_{n}\right)$ is then multidimensional Gaussian
\[
\left(\begin{array}{c}
f\left(\boldsymbol{x}_{1}\right) \\
\vdots \\
f\left(\boldsymbol{x}_{n}\right)
\end{array}\right) \sim \mathrm{N}\left(\left(\begin{array}{c}
m\left(\boldsymbol{x}_{1}\right) \\
\vdots \\
m\left(\boldsymbol{x}_{n}\right)
\end{array}\right),\left(\begin{array}{ccc}
k\left(\boldsymbol{x}_{1}, \boldsymbol{x}_{1}\right) & \cdots & k\left(\boldsymbol{x}_{1}, \boldsymbol{x}_{n}\right) \\
\vdots & \ddots & \\
k\left(\boldsymbol{x}_{n}, \boldsymbol{x}_{1}\right) & & k\left(\boldsymbol{x}_{n}, \boldsymbol{x}_{n}\right)
\end{array}\right)\right)
\]
A Gaussian process is stationary if its mean is constant and the two-argument covariance function is of the form,
\[
k\left(\boldsymbol{x}, \boldsymbol{x}^{\prime}\right)=C\left(\boldsymbol{x}^{\prime}-\boldsymbol{x}\right)
\]
where $C(\boldsymbol{x})$ is another function, the stationary covariance function of the process.

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=9 cm]{example_GRF}
\end{center} 
   \begin{center} 
  \includegraphics[width= 7cm]{example_GRF_1}
\end{center}
\caption{A 2D Gaussian random field with Mat\'ern covariance (Top);  Mat\'ern covariance as function of distance $r$ (Botton).}\label{im1}
\end{figure}

To obtain a heterogeneity on $G_c$, we  use a Gaussian random field generated by the multidimensional Mat\'ern covariance function,
\[
C(r)=\sigma^{2} \frac{2^{1-v}}{\Gamma(\nu)}\left(\sqrt{2 \nu} \frac{r}{l}\right)^{v} K_{v}\left(\sqrt{2 v} \frac{r}{l}\right)
\]
 
where $r=\left\|\boldsymbol{\xi}-\boldsymbol{\xi}^{\prime}\right\|$,  $ \boldsymbol{\xi}=\left(x_{1}, x_{2}, \ldots, x_{d-1}, t\right) \in \mathbb{R}^{d}$,  $\Gamma$ is the gamma function, $K_{\nu}$ is the modified Bessel function of the second kind, and $l$ and $\nu$ are positive parameters of the covariance. 

A non-complex way to obtain a random field follows the classical Wiener-Khinchin theorem that states that the stationary
covariance function of a process can be expressed as the inverse Fourier transform of the spectral density~\cite{Sarkka_2013},
\[
C(t)=\mathcal{F}^{-1}[S(\omega)]=\frac{1}{2 \pi} \int S(\omega) \exp (i \omega t) \mathrm{d} \omega.
\]

Applying this idea, we can recover a Gaussian random field with Mat\'ern covariance using the corresponding spectral density for the Mat\'ern covariance matrix,
\[
S\left(\omega_{r}\right)=S\left(\omega_{x}, \omega_{t}\right) \propto \frac{1}{\left(\lambda^{2}+\left\|\omega_{x}\right\|^{2}+\omega_{t}^{2}\right)^{v+d / 2}},
\]
where $\lambda=\sqrt{2 v} / l$.


 
\section{FEniCS implementation}\label{sec:4}

In this section, we describe the numerical methodology to seek
solutions of deterministic problems~\eqref{eq:minimization}
and~\eqref{ec3}. In both cases, we employ FEniCS open-source
finite-element software to obtain numerical simulations. 

\subsection{Numerical methodology for variational approach}

In order to solve the variational problem~\eqref{eq:minimization}, a iterative algorithm known as alternate minimization \cite{BOURDIN_2000, MARIGO_2016} is presented in algorithm \ref{algoritmo1}, we assume a uniform time step $\Delta t$, this numerical technique uses the displacement and the damage field $\left(\mathbf{u}_{i-1}, \phi_{i-1}\right)$ at time step $t_{i-1},$ the solution at time step $t_{i}$ is obtained by solving the following bound-constrained minimization problem
\begin{equation}\label{ec:min}
 \inf \left\{ E_{\ell,t_{i}}(\mathbf{u}, \phi)|  \mathbf{u} \in \mathcal{C}_{t_{i}},  \phi \in \mathcal{D}_{t_{i}}\right\}
 \end{equation}
 
 
where $\mathcal{C}_{t_{i}}(\Omega)= \left\{u \in H^{1}\left(\Omega \right)| u=\mathbf{u}_{i-1}\right.$ on $\left.\partial_{\mathbf{u}_{i-1}} \Omega\right\}$, and $\mathcal{D}_{t_{i}}=\left\{\phi \in H^{1}(\Omega): \phi(x) \geq \phi_{i-1}\right. \text{a.e. in }\Omega, \phi=\phi_{i-1}$ on $\partial_{\phi_{i-1}} \Omega\}$. The constraint $\phi(x) \geq \phi_{i-1}$ is the time-discrete version of the irreversibility of damage. Following previous analysis in \cite{AMBROSIO_1992} and \cite{BOURDIN_2000}, for  $a(\alpha)=(1-\alpha)^{2}+\eta_{\ell}$ and $w(\alpha)=\alpha^{2},$ with $\eta_{\ell}=o(\ell)$ it is possible to show through asymptotic methods ( $\Gamma$ -convergence) that the solutions of the global minimization problem \eqref{ec:min} tend to the solutions of the global minimization problem \eqref{eq21} as $\ell \rightarrow 0$ \cite{GIACOMINI_2005}. 

\begin{algorithm}
\caption{Alternate minimization}\label{algoritmo1}
\begin{algorithmic}
\STATE { { Given } $\left(\mathbf{u}_{i-1}, \phi_{i-1}\right)$ , the state at the previous loading step. }\\
\STATE {{ Set } $\left(\mathbf{u}^{(0)}, \phi^{(0)}\right):=\left(\mathbf{u}_{i-1}, \phi_{i-1}\right)$}\\
\WHILE {not converged}
\STATE { { Find } $\mathbf{u}^{(p)}:=\underset{\mathbf{u} \in \mathcal{C}_{t_{i}}}\arg \min  {E_{\ell}}_{t_{i}}\left(\mathbf{u}, \phi^{(p-1)}\right)$ }\\
\STATE{{ Find } $\phi^{(p)}:=\underset{\phi \in \mathcal{D}_{t_{i}}}\arg \min  {E_{\ell}}_{t_{i}}\left(\mathbf{u}^{(p)}, \phi\right)$} \\
\ENDWHILE
\STATE {{ Set } $\left(\mathbf{u}_{i}, \phi_{i}\right)=\left(\mathbf{u}^{(p)}, \phi^{(p)}\right) $}\\
\STATE { End}\\
\end{algorithmic}
\end{algorithm}


\subsection{Numerical methodology for hybrid phase-field approach}\label{sec:hfpm}

 The strategy to solve the hybrid phase-field method is the staggered scheme, this algorithm considers at each load step $\Delta u$, we compute the displacement field $\mathbf{u}$ and the damage field $\phi$ alternately until convergence. 

The hybrid phase field method \eqref{ec3}, is solved using the finite element method. We consider the infinite dimensional trial $(U, P)$ and test spaces $(V^0, Q)$. Let $W(\Omega)$ include the linear displacement field and phase field variable
\begin{align*}
  (U, V^0) &:= \left \{ ({\bf{u,v}})\in [C^0(\Omega)] \right. : \\
 ({\bf{u,v}}) \in [W(\Omega)]^d \subseteq & \left. [H^1(\Omega)]^d, {\bf{u}} = \overline{{\bf{u}}}, {\bf{v}}= 0 \text{  on   } \Gamma_u \right\}\\
 (P, Q^0) & :=\left \{ (\phi, q) \in [C^0 (\Omega)]^d \right.: \\
 (\phi, q) \in & \left.  [W(\Omega)]^d , (\phi, q) = 0, \text{   on   } \Gamma_{\phi} \right \}
\end{align*}



The standard Bubnov-Galerkin procedure leads to the weak formulation:
Find ${\bf{u}}\in U$ and $\phi \in P$ such that, for all ${\bf{v}}\in V^0$ and $\theta \in Q,$
\begin{equation}
 A({\bf{u,v}}) = f({\bf{v}}), \quad B(\phi, \theta) = g(\theta),
\end{equation}
where, 
\begin{equation}
 A({\bf{u,v}}) = \int_{\Omega} \left[ (1 - \phi)^2 + \eta_{\ell} \right] {\bf{\sigma(u):\varepsilon(v)}}d\Omega 
\end{equation}
\begin{equation}
 f({\bf{v}}) = \int_{\Omega} {\bf{b \cdot v}} d\Omega + \int_{\Gamma_t} {\bf{\overline{t} \cdot v}} d\Gamma,
\end{equation}
\begin{equation}
 B({\bf{u,v}}) = \int_{\Omega} \left \{ \nabla \theta G_c l_0 \nabla \phi + \theta \left[ \frac{G_c}{l_o} + 2 H^+ \right] \phi  \right \}d\Omega 
\end{equation}
\begin{equation}
 g(\theta) = \int_{\Omega} {2 H^+} \theta d\Omega + \int_{\Gamma} \nabla \phi \cdot {\bf{n}}\theta d\Gamma,
\end{equation}

\begin{algorithm}
\caption{Staggered Scheme} \label{algoritmo2}
\begin{algorithmic}
\STATE { { Given } $\left(\mathbf{u}^n, \phi^{n}\right)$ , the state at the previous loading step. }\\
\STATE {Set $\left(\mathbf{u}^{0}, \phi^{0}\right)(\mathbf{x})=(\mathbf{0}, 0)$ for all $\mathbf{x}$ in $\Omega$ and $\left[\phi^{n+1}\right]^{0}=\phi^{n}$ for $n>0$}\\
\WHILE {not converged}
\STATE { Compute $\left[\mathbf{u}^{n+1}\right]^{i+1}$ by solving the equation
$$
\nabla \cdot \boldsymbol{\sigma}\left(\left[\mathbf{u}^{n+1}\right]^{i+1},\left[\phi^{n+1}\right]^{i}\right)=0 \text { in } \Omega
$$ with $\text { boundary conditions } \boldsymbol{\sigma} \cdot \mathbf{n}=\mathbf{t}^{n+1} \text { on } \Gamma_{N},\left[\mathbf{u}^{n+1}\right]^{i+1}=\mathbf{u}_{D}^{n+1} \text { on } \Gamma_{D} \text { . }$}\\
\STATE{Update the history field $\left[\mathcal{H}^{n+1}\right]^{i+1}=\max \left(\mathcal{H}^{n},\left[\Psi_{0}^{n+1}\right]^{i+1}\right)$.} \\
\STATE{Compute $\left[\phi^{n+1}\right]^{i+1}$ by solving
$$
-\ell^{2} \Delta\left[\phi^{n+1}\right]^{i+1}+\left[\phi^{n+1}\right]^{i+1}=\frac{2 \ell}{G_{c}}\left(1-\left[\phi^{n+1}\right]^{i+1}\right)\left[\mathcal{H}^{n+1}\right]^{i+1} \text { in } \Omega
$$
with boundary condition $\left(\nabla\left[\phi^{n+1}\right]^{i+1}\right) \cdot \mathbf{n}=0$ on $\partial \Omega$}
\ENDWHILE
\STATE { End}\\
\end{algorithmic}
\end{algorithm}

Generally, the linear matrix system is solved with the staggered scheme \cite{HIRSHIKESH_2019, HUYNH_2019}, see Algorithm \ref{algoritmo2}.
This strategy requires iterating over $i$ until some stopping criterion indicating convergence is satisfied, or alternatively we should to take sufficiently small increments $\Delta u$ \cite{MUIXI_2020}, this is, we require small steps to obtain a correct approximation of the crack, implying high computational costs. An analysis about  iterations in the staggered scheme is presented in \cite{AMBATI_2015}. 





\section{Results and simulations}\label{sec:5}

\subsection{Calibration of the phase field model}

To calibrate our algorithms, let us consider a square plate with a straight edge crack (standard for calibrating algorithms, \cite{HIRSHIKESH_2019}). We use two kinds of test problems: first, the pure mode-I fracture, the plate is subjected to displacement at the top in the y-direction to simulate a brittle fracture, the geometry and boundary conditions are in Figure \ref{figcf}-a). The expected fracture propagation is a crack in the direction of the initial fracture. The second test, the shear mode-II fracture, considers the same plate loaded in shear mode, see Figure \ref{figcf}-b) to observe boundary conditions, in this case, the expected propagation of the initial crack is a curve in the bottom-right side of the specimen.
  
  \begin{table}[h]
\indexspace
\begin{center}
\begin{tabular}{cc}\hline\hline\ 
Parameter & Value \\
\hline
$\lambda$ & 121.15 kN/mm$^2$ \\
$\mu$ & 80.0 kN/mm$^2$ \\
$G_c$ &  2.7 \\
$\ell$ & 0.01104 \\
\hline
\end{tabular}
\caption{Fixed parameters for simulations.}\label{t_mpt}
\end{center}
\end{table}
  
Simulations with the hybrid phase-field model was obtained with the set of parameters presented in Table \ref{t_mpt}. For the mode-I fracture, we use 131,072 triangular elements that correspond to $\ell = 0.011$, and we apply a displacement increments of $\Delta u = 1\times 10^{-5}$ mm up to $u = 5\times 10^{-3}$, and  $\Delta u = 1\times 10^{-6}$ up to $u= 6\times 10^{-3}$. The numerical propagation of the initial crack is similar to the expected propagation, see Figure \ref{frac1}.The displacement-force plot is presented in Figure \ref{frac_comp_}, we observe the similar behaviour of the simulated profile and the reported profile in \cite{AMBATI_2015}.

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{cuadro.png}
\end{center} 
\caption{Plate with edge crack: Geometry and boundary conditions for (a) Tension and (b) Shear test. Dimensions in mm.}\label{figcf}
\end{figure}

The second example, the shear mode-II fracture, used the same discretization of the pure Mode-I fracture. The resulting propagation of the initial crack is a special case in which we study the parameter dependencies over simulations (mesh size, $\ell$, $\Delta u$). It is well known that the solutions to brittle fracture problems are notoriously mesh-dependent, several simulations are shown in different researches \cite{MIEHE_2010, HIRSHIKESH_2019, AMBATI_2015,  HUYNH_2019}, each author fix values of mesh size, $\ell$, and $\Delta u$, in each research, the authors obtain different profiles that try to explain the behavior of the uncertain crack. Different reported parameter configurations for numerical simulations are shown in table \ref{tabla_sim}. For example, force-displacement profiles with a beak are found in~\cite{AMBATI_2015, HIRSHIKESH_2019, HUYNH_2019}, while mountain-like profiles are presented in~\cite{MIEHE_2010,  WEINBERG_2017}. 



\begin{figure}[H]
\begin{multicols}{2}
    \begin{center} 
  \includegraphics[width=5 cm]{phi_sr_1_}
\end{center} 
   \begin{center} 
  \includegraphics[width=5 cm]{phi_sr_2_}
\end{center}
\end{multicols}
\caption{Status of the mode-I fracture of the hybrid P-F model at displacements of $5.3\times10^{-3}$ mm (Left) and $5.4\times10^{-3}$ mm (Right).}\label{frac1}
\end{figure}

\begin{figure}[H]
\begin{center} 
  \includegraphics[width=8 cm]{comparative_Mod1}
\end{center}
\caption{Displacement-force of the hybrid P-F model and the reported in \cite{AMBATI_2015}. }
\label{frac_comp_}
\end{figure}

  \begin{table}[h]
\indexspace
\begin{center}
\caption{Fixed parameters for simulations. $^a$ = triangular elements, $^b$ = quiadrilateral elements, $^* =$  refinement in zones where crack is expected to grow, $^{**}$=Authors uses the anisotropic model from Amor {\sl et al,} 2009 \cite{AMOR_2009}.}\label{tabla_sim}
\resizebox{9cm}{!}{ 
\begin{tabular}{p{108pt}p{55pt}p{38pt}p{33pt}p{30pt}p{38pt}
\hline\hline\ 
Researches & Model & FEM & $h$ (mm) & $\ell$  (mm) & $\Delta u$ (mm)\\
\hline
Miehe, 2010 \cite{MIEHE_2010} & Hybrid P-F & 30000$^a$ & $\approx 0.002$ & $0.015$ & $1 \times 10^{-5}$\\

Ambati, 2015 \cite{AMBATI_2015} & Hybrid P-F & 20592$^b$ & $0.001$ & $0.004$ & $1 \times 10^{-5}$\\ 

Hirshikesh, 2017 \cite{HIRSHIKESH_2019} & Isotropic  & 131072$^a$ & $0.0055$ &  $0.011$ & $1 \times 10^{-5}$\\

Hirshikesh, 2017 \cite{HIRSHIKESH_2019} & Hybrid P-F & 131072$^a$ & $0.0055$ &  $0.011$ & $1 \times 10^{-5}$\\

Hirshikesh${**}$, 2017 \cite{HIRSHIKESH_2019, AMOR_2009} & Anisotropic & 131072$^a$ & $0.0055$ &  $0.011$ & $1 \times 10^{-5}$\\

Huynh, 2019 \cite{HUYNH_2019} & Hybrid P-F& 25000$^{a,*}$  &  $\approx 0.004$  & $0.008$ & $1 \times 10^{-5}$ \\

This work & Hybrid P-F 1 & 131072$^a$ & $0.0055$ & $0.0087$ & $1 \times 10^{-5}$ \\

This work & Hybrid P-F 2 & 131072$^a$ & $0.0055$ & $0.01104 $ & $1 \times 10^{-4}$\\
\hline
\end{tabular}
}
\end{center}
\end{table}

We present our results of two different configurations of parameters for the shear Mode-II fracture, see Table \ref{tabla_sim}. In our first simulation, we assume $\ell = 0.0087$ with a displacement increment of $\Delta u=1 \times 10^{-4} $ mm up to $u = 1.5\times ^{-2}$ mm. Here, the force-displacement profile agrees well with the results in ~\cite{MIEHE_2010, WEINBERG_2017}. For our second simulation we use $\ell = 0.0116$ with a displacement increment of $\Delta u=1 \times 10^{-5} \mathrm{mm}$, the resulting force-displacement profile contains a beak around the displacement of $0.011$~mm. Our simulation methodology is capable of capturing the essential features of brittle fracture (c.f. Figure~\ref{frac_shear}). Figure \ref{frac_comp} shows the comparative profiles reported by other authors \cite{HIRSHIKESH_2019, AMBATI_2015}, we observe different profiles after $u = 0.01$ mm. This behavior describes the manner, intensity, and velocity of the generated crack in the specimen. The displacement-force curve present a different profile if the mesh is refined over the region where the crack is expected \cite{HUYNH_2019}, and a maximum value around to $u = 0.08$ mm with a mountain-like profile is reported in \cite{MIEHE_2010a}.

\begin{figure}[H]
\begin{multicols}{2}
    \begin{center} 
  \includegraphics[width=5 cm]{phi_cr_1_}
\end{center} 
   \begin{center} 
  \includegraphics[width=5 cm]{phi_cr_2_}
\end{center}
\end{multicols}
\caption{Status of the shear mode-II fracture with the hybrid P-F model at displacements of $5.3\times10^{-3}$ mm (Left), and  $16\times10^{-3}$ mm (Right).}\label{frac_shear}
\end{figure}


\begin{figure}[H]
\begin{center} 
  \includegraphics[width=8 cm]{comparative_all
\end{center}
\caption{Displacement-force curves for the shear mode-II fracture in this work compared with the reported in the literature, each configuration of mesh size $h$, $\ell$, and $\Delta u$ provides a characteristic profile (Table \ref{tabla_sim}).}\label{frac_comp}
\end{figure}

\subsection{Variational approach as a surrogate mathematical model}

In this section, we propose to employ the variational approach to construct a surrogate model of the hybrid phase-field method to reduce the computational cost and complexity, as is described in section \ref{sec:hfpm}. 
Usually, the computational cost of the phase-field method depends on displacement increment $\Delta u$, mesh size $h$ and $\ell = k*h$.  The value of $k=2$ has been suggested~\cite{MIEHE_2010,HIRSHIKESH_2019}. The computational cost increases when $\Delta u$ is small, for instance, for a maximum displacement of $0.6$ and an incremental displacement of $\Delta u = 1\times10^{-4}$, we need to solve system \eqref{sis_1} 6,000 times.  If we change $\Delta u = 1\times10^{-5}$, 60,000 simulations should be realized.  In addition, the staggered scheme depends on the number of iterations to stabilize each solution according to the incremental load $\Delta u $, \cite{AMBATI_2015}; on the other hand, if the mesh size of the discretization of $\Omega$ decreases, time of execution increases.

In principle, the variational approach should be able to accurately capture fracture as compared to the phase-field model at a
considerably lower computational cost.  Therefore, the variational model presents a viable alternative as a surrogate mathematical model for the phase field model.

We propose the weighted-variational model (W-V model) given by,

\begin{equation}
 E^*_{\ell}(\boldsymbol{u}, \phi) = \xi E_{\ell}(\boldsymbol{u}, \phi),
 \end{equation}
 
where $\xi$ is a real constant, as a surrogate model of the hybrid phase field method.

We compare the results of the phase-field model and the variational model using $\xi = 1$ for shear mode-II fractures. We descretizate the body in 32768 triangular elements, using $\ell=0.016$. The applied displacements were $\Delta_u =t_M/N_{nodes}$, where $t_M =$ maximun load, and $N_{nodes} =$ number of nodes.  In figure \ref{sigma_111}, we present two simulations of shear cracks with $N_{nodes} = \{10, 20\}$, we observe that results are very similar between them, both has the same trajectories for the crack. We observe the imaginary crack discussed by other authors in previous works \cite{MIEHE_2010,  HIRSHIKESH_2019}. Nonetheless, the force-displacement curves are  useful to describe the behavior of the fracture as the displacement is aplied, in this case, the curves are similar, suggesting that the crack is generated around the displacement of $0.011$ mm, see Figure \ref{sigma_com1}.

\begin{figure}[H]
\begin{multicols}{2}
    \begin{center} 
  \includegraphics[width=5 cm]{Nuc_10}
\end{center} 
   \begin{center} 
  \includegraphics[width=5 cm]{Nuc_50}
\end{center}
\end{multicols}
\caption{Comparative results of the W-V model with discretization of the displacement axis of:  10 points - Left, 50 points - Right}\label{sigma_111}
\end{figure}

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{comparative_VM_02}
\end{center} 
\caption{Comparative results of displacement-force obtained from the W-V model using 20 points (blue line, $\approx$ 20 min.)
  and phase field method using 1024 points (red line, $\approx$ 8
  hrs.) .}\label{sigma_com1}
\end{figure}

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{comparative_VM_03}
\end{center} 
\caption{Results of the weighted-variational model with 30 nodes (black line) and 100 nodes (red line), compared with the resported in the literature.}\label{sigma_com2}
\end{figure}


We present a second set of simulations, W-V model 2 and W-V model 3, in this case, we approximate solutions reported in \cite{MIEHE_2010a} and \cite{HUYNH_2019}, see Figure \ref{tabla_sim}. W-V model 2, used the same triangulation than V-W model 1, with number of nodes $N_{nodes} = 30$, we fixed $\xi = 1.28$. W-V model 3 used $N_{nodes} = 100$ and $\xi = 1.26$, this model present the particular firs beak around $u = 0.008$ mm.

In our three configurations of parameters (W-V model a,2 and 3), we approximate the reported profiles in the literature, see Figures \ref{sigma_com1} and \ref{sigma_com2}.



\subsection{Gaussian random field $G_c$}


In this section, we consider a stationary Gaussian random field $G_c$ on $\Omega$ to understand the spatial effect of the randomness on the solutions of $\boldsymbol{u}$  and $\phi$. A
sample of the Gaussian random field with Mat\'ern covariance for
$G_c$ is shown in  Figure~\ref{sigma_1_}. The parameters for simulations of a pure mode I fracture and shear mode II fracture are in Table \ref{t_mpt}, both simulations used $l_0 =$ 0.011040. 

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=7 cm]{random
\end{center} 
\caption{Example of Gaussian random field with Mat\'ern covariance.}\label{sigma_1_}
\end{figure}

We carry out five simulations for each kind of fracture using different Mat\'ern covariances. In Figure \ref{RM_SR1}, we present a propagation path of the initial crack for the pure Mode-I fracture. We observe the perturbations caused by the random field $G_c$ compared with the base trajectory using a fixed $G_c$, see Figure \ref{frac1}. We plotted the baseline displacement-force curve (simulation without random field) compared with the displacement-forces caused by the Gaussian random field.

\begin{figure}[H]
\begin{multicols}{2}
    \begin{center} 
  \includegraphics[width=5 cm]{phi_sr_martern_1}
\end{center} 
    \begin{center} 
  \includegraphics[width=5 cm]{phi_sr_martern_2}
\end{center} 
\end{multicols}
\caption{Trayectory of a mode-I fracture with random field $G_c$ at $5.3\times10^{-3}$ mm (Left); $5.4\times10^{-3}$ mm (Right).}\label{RM_SR1}
\end{figure}

The second set of experiments consider the shear mode-II fracture with displacement increment of $\Delta u$ $=1 \times 10^{-4} \mathrm{mm}$ up to $u = 1.6\times 10^{-2}$ mm. Figure~\ref{RM_CR1} shows a perturbed crack path with the applied Gaussian random field for $G_c$ at the displacement of $u=5.3\times 10^{-3}$ mm (the crack begin its propagation), and when the complete perturbed crack is formed, $u = 1.6\times 10^{-2}$ mm, see Figure \ref{frac_shear} to compare the crack propagation without Gaussian random field. The perturbed displacement-force curves related to the five perturbed cracks from the five samples of the random fields are shown in Figure ~\ref{sigma_CR_M1}, we observe different intensities and velocities when the fracture is created and propagated.  In Figure ~\ref{over_02}, we present the overlap of the five possible paths for the fractures in shear mode. 

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{comparative_M1}
\end{center} 
\caption{Displacement-force cruves of the mode-I fracture of 5 simulations with Mat\'ern covariance (blue lines) compared with the  baseline curve (red line).}\label{sigma_SR_M1} 
\end{figure}

In Table \ref{table_times}, we present the execution times for simulations using the random Gaussian fields $G_c$ with different standard deviation for the covariance matrix, for this, we use the ratio $\frac{G^*_c}{\sigma}$, where $G^*_c$ is the fixed value without perturbation. We observe that the average of the execution times is approximately between 6.5 to 7.7 hours.


\begin{figure}[H]
\begin{multicols}{2}
   \begin{center} 
  \includegraphics[width=5 cm]{phi_cr_martern_1}
\end{center}
   \begin{center} 
  \includegraphics[width=5 cm]{phi_cr_martern_2}
\end{center}
\end{multicols}
\caption{Trayectory of a mode-II fracture with the hybrid P-F model at displacement of $5.3\times10^{-3}$ mm (Left);  $16\times10^{-3}$ mm (Right).}\label{RM_CR1}
\end{figure}

\begin{table}[h]
\indexspace
\begin{center}
\caption{Average of execution time in hours for simulation with different Gaussian random field for $G_c$.}\label{table_times}
\begin{tabular}{c|ccccc}
\hline \hline 
${G^*_c}/{\sigma}$& 1 & 5 & 10 & 50 & 100 \\
\hline
mean & 7.65 &	6.58 &	6.51 &	6.56 &	6.52 \\
standard deviation & 0.32 &	0.03 &	0.01 &	0.02 &	0.03 \\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{comparative_shear}
\end{center} 
\caption{Displacement-force curves of 5 perturbed shear mode-II fractures (blue lines), baseline curve (red line).}\label{sigma_CR_M1}
\end{figure}

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=7 cm]{overlap_SR}
\end{center} 
\caption{Perturbed paths with the hybrid P-F model for the mode-I fracture.}\label{over_01}
\end{figure}

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=7 cm]{overlap_CR}
\end{center} 
\caption{Perturbed paths with the hybrid P-F model for the shear mode-II fracture.}\label{over_02}
\end{figure}

We apply the Gaussian random field in the weighted-variational model. In Figure \ref{sigma_com_00}, we present the overlapping of the 5 paths of crack- propagations of five samples for $G_c$. The curve of displacement-force is plotted in Figure \ref{sigma_com_11}. The mean of the execution time for the hybrid phase-field model is approximately 7 hours, meanwhile, with the variational model we have a mean of the execution time of around 20 minutes, reducing the time complexity around 90\%, this implies a reduction of computational time, storage, and reducing the number of operations in memory.

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=7 cm]{overlap_V5}
\end{center} 
\caption{Overlap of 5 perturbed crack paths with the weighted-variational model.}\label{sigma_com_00}
\end{figure}

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{comparative_VM_04}
\end{center} 
\caption{Comparative paths of 20 simulations of the W-V model 1 for the shear mode-II fracture with different Gaussian random field with Mat\'ern covariance.}\label{sigma_com_11}
\end{figure}

Simulations have shown that our model reproduces the behavior of the crack propagations, similar to the phase-field model, even after adding the Gaussian perturbations in the parameter $G_c$. Figure \ref{sigma_com_11} shows the overlap to the displacement-force curves of 20 samples of perturbed W-V model 1 and compares it with our simulation of our Hybrid P-F 2 and W-V model 1. In Figure \ref{sigma_com_111} we plotted 20 perturbed samples of the W-V model 2, and we compare the results with the profiles presented in \cite{MIEHE_2010a} and \cite{HUYNH_2019}, and with our results of W-V model 2 and W-V model 3.


\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{comparative_VM_01}
\end{center} 
\caption{Comparative paths of 20 simulations of the W-V model 2 for the shear mode-II fracture with different Gaussian random field with Mat\'ern covariance.}\label{sigma_com_111}
\end{figure}





\subsection{Comparisons with a laboratory experiment}

In this subsection, we apply the variational and phase-field models with a Gaussian random field to simulate experimental example of Ambati {\sl et al.} \cite{AMBATI_2015} (Figures \ref{sigma_exp_2} and \ref{sigma_exp_3}). The geometry and boundary conditions are shown in Figure \ref{sigma_exp_2}-a). The specimen was made of cement mortar, composed of $22$\% cement (cement I 32.5: high alumina cement 4:1), $66$\% sand (grain size $<$ $1$~mm), and $12$\% water, leading to a water ratio of $0.55$, more details in \cite{AMBATI_2015}.  The specimen is a notched plate, subject to a load applied by a top pin and a fixed lower pin.  The sample has a hole offset from the center to induce mixed-mode fracture the displacement controlled loading was $0.1$~mm/min (Figure \ref{sigma_exp_2}-b).  

The material parameters $\lambda = 1.94$~kN/mm$^2$, $\mu = 2.45$~kN/mm$^2$, $G_c = 2.8 \times 10^{-3}$~kN/mm are used. Figure \ref{sigma_exp_2} shows the physical experimental with the resulting crack and the expected fracture presented in \cite{AMBATI_2015}.

In our simulations, we discretize the domain $\Omega$ into $4,947$ triangular elements, and we used linear test functions for FEM formulation, we set $\ell = 1.5*l = 2.8008$. One simulation of the fracture propagation using the phase field model with 80 displacement increments ($\Delta u = 0.0275$) takes approximately 45 minutes to complete, while the variational model with $\xi=1.8$, and 50 displacement increments ($\Delta u = 0.044$) has an execution time of about 5 minutes.  The surrogate model is comparable and describes the expected behaviour of the crack, as evident in Figure~\ref{sigma_3exp}.  The simulated displacement-force curves for the phase field and the variational surrogate model present similar profiles, see Figure \ref{sigma_PF_VAR}. These results show clearly that the variational model, calibrated with the adequate set of parameters ($\ell$, $\Delta_u$ and mesh size), can accurately reproduce results of the phase-field model. In our simulation, in both models, the first crack is generated around the applied load of $0.35$~mm.

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=9 cm]{ambati_example.png}
\end{center} 
\caption{Experiment setting of the specimen taken from \cite{AMBATI_2015}.}\label{sigma_exp_2}
\end{figure}


\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=7 cm]{ambati_example2.png}
\end{center} 
\caption{a) Experimental, and b) expected paths of the fracture presented in \cite{AMBATI_2015}.}\label{sigma_exp_3}
\end{figure} 

To simulate non-homogenieties in the material and obtain more realistic simulations, we carried out 20 simulations with different Gaussian random fields for $G_c$. Figure~\ref{GRF_con} shows the overlap of the perturbed fracture trajectories, the results are similar to the expected in Figure \ref{sigma_exp_3}. Comparative displacement-force curves are shown in Figure \ref{GRF_con2}. We observe that the profiles shows that the first crack is reach around 0.35 mm.
Numerical experiments using random fields provide more realistic information about the cracks. In real structures (concrete, ceramic, etc.) the physical properties are essential to obtain resistivity, and durability, for this reason, the study of the random properties in materials that can induce a fracture propagation is important for practitioners.

\begin{figure}[H]
\begin{multicols}{2}
    \begin{center} 
  \includegraphics[width=4 cm]{comparative_concrete_PF}
\end{center} 
    \begin{center} 
  \includegraphics[width=4 cm]{comparative_concrete_VM}
\end{center} 
\end{multicols}
\caption{Comparative results of displacement-force: phase field model - 80 points ($\Delta u = 0.0275$), $\approx$ 45 minutes (Left), weighted-variational model - 50 points ($\Delta u = 0.044$), $\approx$ 5 minutes (Right).}\label{sigma_3exp}
\end{figure}

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{comparative_concrete_0}
\end{center} 
\caption{Comparative results of displacement-force in this research.}\label{sigma_PF_VAR}
\end{figure}

The weighted-variational model as a surrogate model of the phase-field presents an alternative to obtaining simulations in a short time, reducing complexity, and computational storage. We have observed that we recover the expected trajectories of the fractures in our numerical examples.


\section{Conclusion}\label{sec:6}

 

In this work, we present the results of the spatial perturbed $G_c$ with Gaussian random fields using Mat\'ern covariances. We simulate the expected fracture paths using two mathematical models: the hybrid phase-field model and a variational model. We observed the parameter dependence of $\ell$, $\Delta_u$, and mesh size $h$ on the simulations. The calibration of these numerical parameters is necessary to obtain realistic results and comparable with real experiments. We propose the weighted-variational model as a surrogate model of the hybrid phase-field model to reduce the execution time and storage for simulations. With this approach, we can find the value of $\xi$ that can emulate the hybrid phase-field model.

\begin{figure}[H]
\begin{multicols}{2}
    \begin{center} 
  \includegraphics[width=4.1 cm]{comparative_concrete_M}
\end{center} 
   \begin{center} 
  \includegraphics[width=4 cm]{comparative_concrete_2}
\end{center}
\end{multicols}
\caption{Left: Gaussian random field with Mat\'ern covariance; Right:
  Overlap of the 20 simulations of the perturbed fractures.}\label{GRF_con}
\end{figure}

\begin{figure}[H]
    \begin{center} 
  \includegraphics[width=8 cm]{comparative_concrete}
\end{center} 
\caption{Comparative perturbed displacement-force curves 
  (cian lines) with the hybrid (red line) and the W-V model (blue line).}\label{GRF_con2}
\end{figure}

 

In our numerical examples, the weighted-variational model as a surrogate model of the hybrid phase-field approach reduces the execution times in 90\%. The use of random fields to analyze the uncertainty in fractures is a topic that should be analyzed to create new materials and study their strength and resistance. Each possible trajectory of the crack will be affected by the random nature of the specimen, caused by water, air bubbles, and fibers, these anomalies can be described by random fields to simulate more realistic properties on materials. Numerical experiments can provide information to describe crack propagation. This work provides an overview of the use of random fields to analyze the heterogeneous properties of a material. In future work, we will analyze the uncertainty quantification of the elastic parameters applying a Gaussian random field with Mat\'ern covariance as a priori distribution, due these parameters are also affected by the random nature of the specimen.


\section*{Acknowledgments}
 This research is partially founded by CONACYT CB-2016-01-284451 and RDECOM grants.
 
