This note records the paper-to-COMSOL correspondence for the crack-driving energy fai_p, the irreversible history field H, the prescribed initial history field for a pre-existing crack, and the Helmholtz coefficients used to solve the phase field.
The project convention is
\[ \phi=0\quad\text{intact}, \qquad \phi=1\quad\text{fully broken}. \]
Related notes: displacement equation and COMSOL mapping, pressure equation and COMSOL mapping, COMSOL workflow, and theory and equations.
1. Three primary fields and the history variable
The coupled hydrofracture model solves three primary physical fields:
| Field | Physical role |
|---|---|
| \(\boldsymbol{u}\) | displacement; produces strain and mechanical stress |
| \(p\) | pore-fluid pressure; produces Darcy flow and pore-pressure stress |
| \(\phi\) | phase field; describes the diffuse crack state |
The history field \(H\) is an irreversible internal variable rather than a fourth primary physical field. Stress, strain, and crack-driving energy are derived quantities.
The staggered dependency is
\[ \boldsymbol{u},p,\phi \longrightarrow \boldsymbol{\varepsilon},Y \longrightarrow H \longrightarrow \phi. \]
For one staggered iteration, the displacement/pressure state is solved first, the current crack-driving energy is evaluated, \(H\) is updated, and the phase field is then solved with \(H\) held fixed.
2. Paper definition of the crack-driving energy
2.1 Tensile elastic energy
The 2018 COMSOL paper introduces the tensile/compressive split and defines the tensile elastic energy as
\[ \psi_e^+ = \frac{\lambda}{2} \left\langle\operatorname{tr}\boldsymbol{\varepsilon}\right\rangle_+^2 + \mu\sum_{a=1}^{d} \left\langle\varepsilon_a\right\rangle_+^2, \]
where \(\varepsilon_a\) are principal strains and \(\langle x\rangle_+=\max(x,0)\). This is the tensile part that drives phase-field evolution.
The history field in the 2018 formulation is the maximum tensile energy reached in the loading history:
\[ H(\boldsymbol{x},t) = \max_{s\in[0,t]} \psi_e^+(\boldsymbol{x},s). \]
2.2 Initial-stress hydrofracture extension
The 2020 hydrofracture paper adds the contribution of the initial in-situ stress. The current driving quantity can be written as
\[ Y(\boldsymbol{x},t) = \psi_e^+(\boldsymbol{\varepsilon}) + \boldsymbol{\sigma}^{0}:\boldsymbol{\varepsilon}. \]
The corresponding history field is
\[ H(\boldsymbol{x},t) = \max_{s\in[0,t]} Y(\boldsymbol{x},s). \]
At discrete time steps,
\[ H_i^{n+1} = \max\left(H_i^n,Y_i^{n+1}\right). \]
This maximum operation is the irreversibility condition. During unloading, \(Y\) can decrease but \(H\) must not decrease, so the phase field cannot heal through a lower driving energy.
3. Project variable fai_p
The project defines fai_p in Phase_Field_study_Definitions.txt as
fai_p = lanta*tra1_p^2/2
+ mu*(e1_p^2+e2_p^2+e3_p^2)
- (Sh*ex+Sv*ey+Sn*ez)In the COMSOL model this variable belongs under Component 1 > Definitions > Variables 1. It is a user-defined variable, not a COMSOL built-in quantity.
The one-to-one mapping is
| Project variable | Mathematical meaning | Paper term |
|---|---|---|
| lanta*tra1_p^2/2 | \(\frac{\lambda}{2}\langle\operatorname{tr}\boldsymbol{\varepsilon}\rangle_+^2\) | tensile volumetric energy |
| mu*(e1_p2+e2_p2+e3_p^2) | \(\mu\sum_a\langle\varepsilon_a\rangle_+^2\) | tensile principal-strain energy |
| -(Shex+Svey+Sn*ez) | \(\boldsymbol{\sigma}^{0}:\boldsymbol{\varepsilon}\) | initial-stress contribution |
| fai_p | \(Y_{\mathrm{current}}\) | current crack-driving energy |
The project definitions used to assemble this expression are
e1_p = if(solid.ep1>0,solid.ep1,0)
e2_p = if(solid.ep2>0,solid.ep2,0)
e3_p = if(solid.ep3>0,solid.ep3,0)
tra1 = solid.ep1+solid.ep2+solid.ep3
tra1_p = if(tra1>0,tra1,0)
ex = solid.eel11
ey = solid.eel22
ez = solid.eel33If compressive initial stresses are stored as positive magnitudes, the initial stress tensor is
\[ \boldsymbol{\sigma}^{0} = \begin{bmatrix} -S_h&0&0\\ 0&-S_v&0\\ 0&0&-S_n \end{bmatrix}, \]
and therefore
\[ \boldsymbol{\sigma}^{0}:\boldsymbol{\varepsilon} = -S_h\varepsilon_x-S_v\varepsilon_y-S_n\varepsilon_z. \]
Thus the project expression is
\[ \text{fai\_p} = \psi_e^+ + \boldsymbol{\sigma}^{0}:\boldsymbol{\varepsilon}. \]
The units of fai_p are energy density, normally Pa or J/m^3. It is not the phase field phi.
4. Prescribed initial history field for a pre-existing crack
The expression
H0 = 0.5*B*Gc/l0*(1-2*abs(y-Dd)/l0)is an initial condition used to seed a pre-existing crack. It is not the normal time-history update. Its direct source is the 2018 COMSOL paper, Section 3.2, Eq. (28), where an initial history field is prescribed for an induced crack.
The paper’s coordinate form is adapted in the project to a horizontal crack centered at \(y=D_d\):
\[ H_0(y) = \frac{B G_c}{2l_0} \left( 1-\frac{2|y-D_d|}{l_0} \right). \]
This is a triangular, or tent-shaped, history field:
\[ H_0(D_d)=\frac{B G_c}{2l_0}, \]
and
\[ H_0=0 \quad\text{when}\quad |y-D_d|=\frac{l_0}{2}. \]
The factors have distinct roles:
| Factor | Role |
|---|---|
| \(G_c/l_0\) | energy-density scale |
| \(B\) | dimensionless amplitude control |
| \(1/2\) | gives the amplitude \(B G_c/(2l_0)\) |
| \(1-2|y-D_d|/l_0\) | linear taper across the crack band |
| \(D_d\) | center coordinate of the pre-existing crack |
| \(l_0\) | width of the diffuse crack band |
The reason for choosing a large amplitude is seen from the phase-field equation. Define
\[ A_{\phi} = \frac{2l_0(1-k)H}{G_c}. \]
If the gradient term is ignored locally,
\[ \phi\approx\frac{A_{\phi}}{1+A_{\phi}}. \]
At the center of the seeded crack,
\[ A_{\phi,\max} = \frac{2l_0(1-k)}{G_c} \frac{B G_c}{2l_0} = B(1-k). \]
With the current project definition
B = 1000/(1-k)the center of the seed gives
\[ A_{\phi,\max}=1000, \qquad \phi_{\max}\approx\frac{1000}{1001}\approx0.999. \]
The seed therefore makes the pre-existing crack almost fully broken under the project convention \(\phi=1\) fully broken.
The symbol in this spatial profile must be the geometric coordinate, normally COMSOL’s lowercase y. It is not the crack-driving energy \(Y\). Using the energy \(Y\) in abs(Y-Dd) would be dimensionally invalid because \(Y\) has units of Pa while \(D_d\) has units of length.
Use a cutoff so the initial history field does not become negative outside the band:
if(abs(y-Dd)<l0/2,
0.5*B*Gc/l0*(1-2*abs(y-Dd)/l0),
0)The horizontal length of the pre-existing crack must additionally be selected through the geometry, domain selection, or an x condition. The width expression alone controls the transverse band profile.
5. COMSOL residual for the history field
The paper requires the discrete update
\[ H^{n+1}=\max(H^n,Y^{n+1}). \]
The project expression
H-nojac(if(fai_p>H,fai_p,H))is an algebraic residual form of
\[ R_H = H-\max(\text{fai\_p},H)=0. \]
The if expression is the COMSOL equivalent of a maximum:
if(fai_p>H,fai_p,H) = max(fai_p,H)Its intended active-set behavior is
if fai_p > H: H = fai_p
if fai_p <= H: H is not reducednojac(…) only removes the derivative of the inner if/max expression from the Jacobian. It does not create a previous-time value and does not, by itself, guarantee irreversibility.
For an explicit time-step implementation, the old history should be represented separately:
H-nojac(if(fai_p>H_old,fai_p,H_old))or
H-nojac(max(H_old,fai_p))where H_old is the previous accepted history supplied by the Previous Solution mechanism or another explicitly defined old-state variable.
If the current unknown H is used on both sides, then in the branch fai_p<=H the residual becomes identically zero for any \(H\geq fai_p\). This is a fixed-point-like implementation and can be non-unique unless COMSOL’s time-history and initial-value handling supplies the intended old state. Verify pointwise that
\[ H^{n+1}-H^n\geq0. \]
6. Helmholtz phase-field mapping
The phase-field equation used by the papers and the project is
\[ \nabla\cdot\left(-l_0^2\nabla\phi\right) + \left[ 1+\frac{2l_0(1-k)H}{G_c} \right]\phi = \frac{2l_0(1-k)H}{G_c}. \]
Define
\[ A_{\phi}=\frac{2l_0(1-k)H}{G_c}. \]
COMSOL’s Helmholtz form is
div(-c*grad(phi)) + a*phi = fso the correct settings are
c = l0^2
a = 1+A_phi
f = A_phior directly,
c = l0^2
a = 1+2*H*l0*(1-k)/Gc
f = 2*H*l0*(1-k)/GcThe repeated expression is the same dimensionless coefficient in both places. The absorption coefficient has an additional +1 because the local phase-field term contains \((1+A_{\phi})\phi\), while the right-hand side contains only \(A_{\phi}\).
With the current parameters,
l0 = 0.1 m
Gc = 500 N/m
k = 1e-9the coefficient is approximately
\[ A_{\phi}\approx4\times10^{-4}H[\mathrm{Pa}]. \]
The diffusion coefficient must be \(l_0^2=10^{-2}\,\mathrm{m^2}\) for the current \(l_0\). A literal 102[m2] would represent a different and much larger regularization length.
7. Staggered implementation and checks
At time step \(n+1\):
Hold \(\phi^j\) and \(H^j\) fixed and solve displacement and pressure.
Compute principal strains and fai_p, including the initial-stress contribution.
Update the history field from the old accepted value:
\[ H^{j+1}=\max(H^n,\text{fai\_p}^{j+1}). \]
Hold \(H^{j+1}\) fixed and solve the Helmholtz equation for \(\phi^{j+1}\).
Update \(g(\phi)\), permeability, Biot coefficient, and the mechanical tangent.
Repeat until displacement, pressure, history, and phase-field changes satisfy tolerance.
The minimum engineering checks are:
1. fai_p has units Pa or J/m^3.
2. H has the same units as fai_p.
3. H_new-H_old >= 0 pointwise.
4. phi remains consistent with phi=0 intact and phi=1 broken.
5. A_phi is dimensionless.
6. c=l0^2 with unit m^2.
7. a-f=1 for the unscaled Helmholtz form.
8. The seeded H0 is nonnegative and restricted to the intended crack band.
9. Mesh size resolves l0; the current study uses hmax=l0/2.Sources
- Zhou, Rabczuk, and Zhuang (2018), Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies, Sections 2.2–2.3 and 3.1–3.3, especially Eqs. (7), (14), (15), (27), and (28).
- Zhou, Zhuang, and Rabczuk (2020), Phase field method for quasi-static hydro-fracture in porous media under stress boundary condition considering the effect of initial stress field, Sections 2.1–2.5, especially Eqs. (5), (10), (12), (20), and (21).
- Project definitions: C:/Users/cupk2026/Documents/Codex/Phase field/Reference/Phase_Field_study_Definitions.txt.
- Project parameters: C:/Users/cupk2026/Documents/Codex/Phase field/Reference/Phase_Field_study_Parameters.txt.