1. What exactly is the phase-field method?
The phase-field method converts a geometrically infinitely sharp crack surface into a continuous scalar field \(\phi(\boldsymbol{x})\).
This project uses
\[ \phi=0:\text{intact material}, \qquad \phi=1:\text{fully broken material}. \]
Therefore, the crack is no longer a geometric line that must be tracked. Instead, it is a narrow band in which \(\phi\) transitions smoothly from \(\phi=0\) to \(\phi=1\). The width of this band is controlled by \(l_0\).
The central idea of the phase-field method is:
Do not track the crack surface directly. Let the crack initiate, propagate, branch, and merge automatically through minimization of the total energy.
2. Starting from Griffith fracture theory
Classical Griffith theory can be written as
\[ \Pi(\boldsymbol{u},\Gamma) = \int_{\Omega\setminus\Gamma} \psi_e(\boldsymbol{\varepsilon})\,d\Omega + G_c|\Gamma| - W_{\mathrm{ext}}, \]
where:
- \(\boldsymbol{u}\) is displacement;
- \(\boldsymbol{\varepsilon}\) is strain;
- \(\Gamma\) is the actual crack surface;
- \(G_c|\Gamma|\) is the fracture energy required to create the crack.
The difficulty is that \(\Gamma\) is a discontinuous geometric surface that changes continuously during crack growth. The crack would need to be tracked and the mesh would need to be reconstructed.
The phase-field method replaces \(\Gamma\) with the continuous variable \(\phi\) and approximates the crack surface length or area using a crack-surface density function:
\[ \gamma_{l_0}(\phi,\nabla\phi) = \frac{1}{2l_0} \left( \phi^2+l_0^2|\nabla\phi|^2 \right). \]
Thus,
\[ |\Gamma| \approx \int_\Omega \gamma_{l_0}(\phi,\nabla\phi)\,d\Omega. \]
The fracture energy becomes
\[ \mathcal{E}_{\mathrm{frac}} = \int_\Omega \frac{G_c}{2l_0} \left( \phi^2+l_0^2|\nabla\phi|^2 \right)d\Omega. \]
This is the fundamental reason that \(G_c\), \(l_0^2\), and the gradient term appear in the phase-field equation.
3. Why does the stiffness degradation function \(g(\phi)\) appear?
The crack region cannot retain the tensile stiffness of intact material, so the following degradation function is defined:
\[ g(\phi)=(1-k)(1-\phi)^2+k. \]
Here,
\[ g(0)=1, \qquad g(1)=k\approx0. \]
Therefore:
- \(\phi=0\): intact material, with the full stiffness retained;
- \(\phi=1\): broken material, with tensile stiffness approaching zero;
- \(k\): a small residual stiffness used to prevent the stiffness matrix from becoming exactly singular.
The tensile and compressive energies are separated:
\[ \psi_e = g(\phi)\psi_e^+ + \psi_e^-. \]
Only the tensile energy \(\psi_e^+\) is degraded. The compressive energy is usually not degraded, in order to avoid generating tensile cracks incorrectly under a purely compressive state.
4. How is the phase-field governing equation derived?
To prevent crack recovery during unloading, the current tensile energy is replaced by the history field \(H\):
\[ H(\boldsymbol{x},t) = \max_{s\in[0,t]}Y(\boldsymbol{x},s). \]
First consider the energy related to the phase field:
\[ \mathcal{E}_\phi = \int_\Omega \left[ g(\phi)H + \frac{G_c}{2l_0} \left( \phi^2+l_0^2|\nabla\phi|^2 \right) \right]d\Omega. \]
Take the variation with respect to \(\phi\).
Because
\[ g'(\phi)=-2(1-k)(1-\phi), \]
we obtain
\[ -2(1-k)(1-\phi)H + \frac{G_c}{l_0}\phi - G_cl_0\nabla^2\phi =0. \]
Multiplying both sides by \(l_0/G_c\) gives
\[ -l_0^2\nabla^2\phi + \phi - \frac{2l_0(1-k)H}{G_c}(1-\phi) =0. \]
Define
\[ A_\phi= \frac{2l_0(1-k)H}{G_c}, \]
and rearrange:
\[ -l_0^2\nabla^2\phi + (1+A_\phi)\phi = A_\phi. \]
This is the Helmholtz equation used in COMSOL:
\[ \nabla\cdot(-l_0^2\nabla\phi) + (1+A_\phi)\phi = A_\phi. \]
Therefore, the COMSOL settings naturally become
c = l0^2
a = 1+A_phi
f = A_phiwhere
A_phi = 2*H*l0*(1-k)/GcThis coefficient is not an empirical fitting parameter. It is obtained by taking the variation of the total energy functional.
5. Where does fai_p come from?
The project definition of fai_p is
fai_p =
lanta*tra1_p^2/2
+ mu*(e1_p^2+e2_p^2+e3_p^2)
- (Sh*ex+Sv*ey+Sn*ez)It corresponds to the current crack-driving energy:
\[ Y_{\mathrm{current}} = \psi_e^+ + \boldsymbol{\sigma}^0: \boldsymbol{\varepsilon}. \]
The first part is
\[ \psi_e^+ = \frac{\lambda}{2} \langle\operatorname{tr}\boldsymbol{\varepsilon}\rangle_+^2 + \mu\sum_a\langle\varepsilon_a\rangle_+^2, \]
which is the tensile elastic energy.
The second part is
\[ \boldsymbol{\sigma}^0:\boldsymbol{\varepsilon} = -S_h\varepsilon_x -S_v\varepsilon_y -S_n\varepsilon_z, \]
which is the energy contribution from the initial in-situ stress and the incremental strain.
Therefore,
fai_p = current driving energy Y_currentand the history field is updated as
\[ H^{n+1} = \max\left(H^n,\texttt{fai\_p}^{n+1}\right). \]
6. What is H-nojac(if(fai_p>H,fai_p,H))?
It writes the maximum-value condition as a COMSOL algebraic residual:
H-nojac(if(fai_p>H,fai_p,H))This is equivalent to
\[ R_H = H-\max(\text{fai\_p},H)=0. \]
Its intended behavior is
fai_p > H: update H = fai_p
fai_p <= H: keep the original HThe function nojac() only prevents COMSOL from differentiating the inner if/max expression when forming the Jacobian. It does not automatically provide the value from the previous time step.
A stricter engineering form should explicitly use the old history:
H-nojac(if(fai_p>H_old,fai_p,H_old))or
H-nojac(max(H_old,fai_p))where H_old comes from the Previous Solution or another explicitly defined old-state variable.
7. What is the initial \(H_0\) for a pre-existing crack?
The expression
H0 = 0.5*B*Gc/l0*(1-2*abs(y-Dd)/l0)is not derived as the normal phase-field governing equation. It is an initial history field used to artificially seed a pre-existing crack.
It is a triangular distribution:
\[ H_0(y) = \frac{BG_c}{2l_0} \left( 1-\frac{2|y-D_d|}{l_0} \right). \]
It reaches its maximum at the crack center and becomes zero at
\[ |y-D_d|=\frac{l_0}{2}. \]
The symbol y here is a geometric coordinate, not the crack-driving energy \(Y\).
At the center,
\[ A_{\phi,\max} = \frac{2l_0(1-k)}{G_c} \frac{BG_c}{2l_0} = B(1-k). \]
With the current definition
B = 1000/(1-k)we obtain
\[ A_{\phi,\max}=1000, \qquad \phi_{\max}\approx\frac{1000}{1001}\approx0.999. \]
Therefore, this initial \(H_0\) makes the specified region almost fully broken from the start, which creates the pre-existing crack.
8. How does pressure couple into the hydraulic-fracture model?
After adding pore pressure, the mechanical stress is usually written as
\[ \boldsymbol{\sigma} = \boldsymbol{\sigma}^e + g(\phi)\boldsymbol{\sigma}^0 - \alpha p\boldsymbol{I}. \]
Pressure affects displacement and strain through \(-\alpha p\boldsymbol{I}\), while displacement affects fluid pressure through volumetric strain:
\[ \varepsilon_{\mathrm{vol}} = \nabla\cdot\boldsymbol{u}. \]
The fluid equation can be written as
\[ \rho S\frac{\partial p}{\partial t} - \nabla\cdot \left( \frac{\rho K_{\mathrm{eff}}}{\mu_{\mathrm{eff}}} \nabla p \right) = q_m - \rho\alpha\chi_r \frac{\partial\varepsilon_{\mathrm{vol}}}{\partial t}. \]
The phase field \(\phi\) changes
\[ K_{\mathrm{eff}}, \qquad \alpha, \qquad g(\phi). \]
Therefore, the complete coupling is
\[ p\rightarrow\boldsymbol{u} \rightarrow\boldsymbol{\varepsilon} \rightarrow Y \rightarrow H \rightarrow\phi \rightarrow K_{\mathrm{eff}},\alpha,g(\phi) \rightarrow p. \]
9. Final numerical solution procedure
At every time step, perform:
Fix \(H\) and \(\phi\), and solve for displacement \(u\) and pressure \(p\).
Compute the tensile energy and the initial-stress energy to obtain fai_p.
Update the history field:
\[ H^{n+1} = \max(H^n,\texttt{fai\_p}). \]
Fix \(H\) and solve the phase field:
\[ \mathbf{K}_\phi\boldsymbol{\Phi} = \mathbf{F}_\phi. \]
Use the new \(\phi\) to update stiffness, permeability, and the Biot coefficient.
Repeat the above process until \(u,p,H,\phi\) converge simultaneously.
Therefore, the essence of the phase-field method can be summarized as:
Use the continuous variable \(\phi\) to regularize the crack geometry, derive the phase-field equation through energy minimization, use the history field \(H\) to enforce crack irreversibility, and obtain crack propagation through staggered finite-element solutions of displacement, pressure, and phase field.