This note establishes a one-to-one correspondence between the displacement equations in the phase-field hydrofracture papers and the equations displayed by COMSOL’s Solid Mechanics interface.
The key distinction is:
COMSOL S = total stress
COMSOL S_el = elastic stress sigma^e
COMSOL C(D) = constitutive tangent DThe pore-pressure stress and the degraded initial stress are added separately. They are not additional terms hidden inside the elastic tangent D.
The project convention is
\[ \phi=0\quad\text{intact}, \qquad \phi=1\quad\text{fully broken}. \]
Related notes: theory and equations, COMSOL workflow.
1. Three primary fields
The coupled hydrofracture model solves three primary 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. Stress, strain, and force are derived quantities, not additional primary fields.
For one staggered iteration, the computational dependency is
\[ \boldsymbol{u},p,\phi \longrightarrow \boldsymbol{\varepsilon},Y \longrightarrow H \longrightarrow \phi. \]
During the displacement subproblem, \(p\) and \(\phi\) are held fixed. Consequently, the displacement tangent is the derivative of the elastic stress with respect to strain; the pressure and initial-stress terms contribute to the residual but not to the elastic \(u\)-\(u\) tangent.
2. Displacement equation in the hydrofracture paper
2.1 Small strain
The 2020 hydrofracture paper defines the infinitesimal strain as
\[ \varepsilon_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right). \]
In tensor notation,
\[ \boldsymbol{\varepsilon} = \frac{1}{2} \left[ \nabla\boldsymbol{u} + (\nabla\boldsymbol{u})^{\mathsf T} \right]. \]
This is the same strain expression shown in the COMSOL Equation panel:
epsilon = 1/2*((grad(u))^T + grad(u))2.2 Elastic stress
The tensile/compressive spectral split is
\[ \boldsymbol{\varepsilon}^{+} = \sum_{a=1}^{d} \langle\varepsilon_a\rangle_{+} \boldsymbol{n}_a\otimes\boldsymbol{n}_a, \]
\[ \boldsymbol{\varepsilon}^{-} = \sum_{a=1}^{d} \langle\varepsilon_a\rangle_{-} \boldsymbol{n}_a\otimes\boldsymbol{n}_a, \]
where \(\varepsilon_a\) and \(\boldsymbol{n}_a\) are the principal strains and principal directions.
The elastic energy is
\[ \psi_e = g(\phi)\psi_e^{+}(\boldsymbol{\varepsilon}) + \psi_e^{-}(\boldsymbol{\varepsilon}), \]
with
\[ g(\phi)=(1-k)(1-\phi)^2+k. \]
Therefore the elastic stress used in the displacement equation is
\[ \boldsymbol{\sigma}^{e} = g(\phi)\frac{\partial\psi_e^{+}}{\partial\boldsymbol{\varepsilon}} + \frac{\partial\psi_e^{-}}{\partial\boldsymbol{\varepsilon}}. \]
For isotropic elasticity, the paper writes this in the form
\[ \boldsymbol{\sigma}^{e} = g(\phi) \left[ \lambda\langle\operatorname{tr}\boldsymbol{\varepsilon}\rangle_{+}\boldsymbol{I} +2\mu\boldsymbol{\varepsilon}^{+} \right] + \lambda\langle\operatorname{tr}\boldsymbol{\varepsilon}\rangle_{-}\boldsymbol{I} +2\mu\boldsymbol{\varepsilon}^{-}. \]
2.3 Total stress and equilibrium
The hydrofracture paper adds the degraded initial stress and the Biot pressure stress:
\[ \boldsymbol{\sigma} = \boldsymbol{\sigma}^{e} + g(\phi)\boldsymbol{\sigma}_0 - \alpha(p-p_{\mathrm{ref}})\boldsymbol{I}. \]
The displacement equilibrium equation is
\[ \nabla\cdot\boldsymbol{\sigma} + \boldsymbol{b} = \boldsymbol{0}. \]
In component notation,
\[ \frac{\partial\sigma_{ij}}{\partial x_j}+b_i=0. \]
The natural traction condition is
\[ \boldsymbol{\sigma}\boldsymbol{n}=\boldsymbol{t} \quad\text{on the Neumann boundary}. \]
3. The equation shown by COMSOL
The supplied COMSOL Equation panel displays the quasi-static form
0 = div(S) + f_vThis is the same equilibrium equation as
\[ \nabla\cdot\boldsymbol{\sigma}+\boldsymbol{b}=\boldsymbol{0}, \]
with the following notation change:
| Paper | COMSOL | Meaning |
|---|---|---|
| \(\boldsymbol{\sigma}\) | \(\boldsymbol{S}\) | total stress |
| \(\boldsymbol{b}\) | \(\boldsymbol{f}_v\) | body-force contribution |
| \(\boldsymbol{\sigma}^{e}\) | \(\boldsymbol{S}_{\mathrm{el}}\) | elastic stress |
The COMSOL material equation is displayed as
S = S_inel + S_el
epsilon_el = epsilon - epsilon_inel
S_el = C:epsilon_el
S_inel = S_0 + S_ext + S_q
epsilon = 1/2*((grad(u))^T + grad(u))
C = C(D)Written explicitly in LaTeX, the same COMSOL material equations are
\[ \begin{aligned} \boldsymbol{S} &= \boldsymbol{S}_{\mathrm{inel}} + \boldsymbol{S}_{\mathrm{el}}, \\[4pt] \boldsymbol{\varepsilon}_{\mathrm{el}} &= \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_{\mathrm{inel}}, \\[4pt] \boldsymbol{S}_{\mathrm{el}} &= \boldsymbol{C}:\boldsymbol{\varepsilon}_{\mathrm{el}}, \\[4pt] \boldsymbol{S}_{\mathrm{inel}} &= \boldsymbol{S}_{0} + \boldsymbol{S}_{\mathrm{ext}} + \boldsymbol{S}_{q}, \\[4pt] \boldsymbol{\varepsilon} &= \frac{1}{2} \left[ \left(\nabla\boldsymbol{u}\right)^{\mathsf T} + \nabla\boldsymbol{u} \right], \\[4pt] \boldsymbol{C} &= \boldsymbol{C}(\boldsymbol{D}). \end{aligned} \]
For this phase-field implementation, the relevant interpretation is
\[ \boldsymbol{S}_{\mathrm{el}} = \boldsymbol{\sigma}^{e}, \qquad \boldsymbol{C}=\boldsymbol{D}, \]
while the total stress is assembled as
\[ \boldsymbol{S} = \boldsymbol{S}_{\mathrm{el}} + \boldsymbol{S}_{\mathrm{ext}}^{(0)} + \boldsymbol{S}_{\mathrm{ext}}^{(p)} + \text{other enabled contributions}. \]
The exact names \(S_0\), \(S_q\), and \(S_{\mathrm{inel}}\) are generic COMSOL bookkeeping. They should not be interpreted as new terms in the paper unless a corresponding physical contribution has been intentionally enabled.
4. One-to-one mapping of the two External Stress nodes
4.1 External Stress 1: pore pressure
The COMSOL node displays
S_ext = -alphaB*(p_A - p_ref)*IThis is exactly the paper term
\[ \boldsymbol{S}_{\mathrm{ext}}^{(p)} = -\alpha_B(p_A-p_{\mathrm{ref}})\boldsymbol{I}. \]
The settings in the screenshot are:
| COMSOL setting | Meaning |
|---|---|
| Stress input | Pore pressure |
| \(p_A\) | user-defined pressure variable \(p\) |
| \(p_{\mathrm{ref}}\) | reference pressure |
| \(\alpha_B\) | phase-dependent Biot coefficient |
The phase-dependent coefficient in the project is
\[ \alpha_B(\phi) = \alpha_r\,\mathrm{COE}_r(\phi) + \alpha_f\,\mathrm{COE}_f(\phi), \]
with \(\alpha_r=0.05\) and \(\alpha_f=1\) in the current definitions.
If p is a gauge pressure, using \(p_{\mathrm{ref}}=0\) is consistent. If p is an absolute pressure, \(p_{\mathrm{ref}}\) must be the absolute reference level.
This pressure stress is part of the total stress and the displacement residual. It is not part of the elastic tangent \(\boldsymbol{D}\) when \(p\) is held fixed during the displacement substep.
4.2 External Stress 2: initial in-situ stress
The screenshot gives the diagonal tensor
[-Sh*cc, 0, 0]
[0, -Sv*cc, 0]
[0, 0, -Sn*cc]This corresponds to
\[ \boldsymbol{S}_{\mathrm{ext}}^{(0)} = g(\phi)\boldsymbol{\sigma}_0, \]
with
\[ \boldsymbol{\sigma}_0 = \begin{bmatrix} -S_h&0&0\\ 0&-S_v&0\\ 0&0&-S_n \end{bmatrix}, \qquad cc=g(\phi). \]
The negative signs mean compression when tensile stress is positive. This is the same \(g(\phi)\boldsymbol{\sigma}_0\) term in the paper’s total stress equation.
The initial stress must not also be entered as an equivalent inelastic strain unless the two representations have been derived to be exactly equivalent. Otherwise the initial stress is counted twice.
5. Origin and use of the large D matrix
The constitutive tangent is
\[ D_{ijkl} = \frac{\partial\sigma^{e}_{ij}} {\partial\varepsilon_{kl}}. \]
The tensile/compressive split gives
\[ D_{ijkl} = \overline{D}_{ijkl} + \widetilde{D}_{ijkl}, \]
where
\[ \overline{D}_{ijkl} = \lambda \left[ g(\phi)H_e(\operatorname{tr}\boldsymbol{\varepsilon}) + H_e(-\operatorname{tr}\boldsymbol{\varepsilon}) \right] \delta_{ij}\delta_{kl}, \]
and
\[ \widetilde{D}_{ijkl} = 2\mu \left[ g(\phi)P^{+}_{ijkl} + P^{-}_{ijkl} \right]. \]
Here \(H_e\) is the Heaviside function. It is not the history field \(H\).
The terms in Phase_Field_study_Definitions.txt have the following meaning:
| Text-file variable | Mathematical role |
|---|---|
e1_p, e2_p, e3_p |
\(\langle\varepsilon_1\rangle_+\), \(\langle\varepsilon_2\rangle_+\), \(\langle\varepsilon_3\rangle_+\) |
tra1 |
\(\operatorname{tr}\boldsymbol{\varepsilon}=\varepsilon_1+\varepsilon_2+\varepsilon_3\) |
tra1_p |
\(\langle\operatorname{tr}\boldsymbol{\varepsilon}\rangle_+\) |
n11–n33 |
components of the three principal directions \(\boldsymbol{n}_a\) |
d11_p–d33_p |
tensile Heaviside indicators \(H_e(\varepsilon_a)\) |
d11_n–d33_n |
compressive indicators \(H_e(-\varepsilon_a)\) |
l1, l2, l3 |
perturbed principal strains for repeated-eigenvalue regularization |
g1111_p etc. |
components of \(P^+_{ijkl}\) |
g1111_n etc. |
components of \(P^-_{ijkl}\) |
cc |
degradation factor \(g(\phi)\) |
dd |
volumetric degradation switch |
g1111 etc. |
final components \(D_{ijkl}\) |
For example,
g1111 = 2*mu*(cc*g1111_p + g1111_n) + lanta*ddmeans
\[ D_{1111} = 2\mu\left[ g(\phi)P^{+}_{1111}+P^{-}_{1111} \right] + \lambda\,dd. \]
The \(\lambda\) term appears in D1111, D1122, D1133, D2222, D2233, and D3333, because \(\delta_{ij}\delta_{kl}=1\) for these components. It does not appear in D1112 or D1212 because a shear index makes the Kronecker-delta product zero.
The 6-by-6 matrix is then only a Voigt arrangement of this fourth-order tensor:
\[ \boldsymbol{D}_e= \begin{bmatrix} D_{1111}&D_{1122}&D_{1133}&D_{1112}&D_{1123}&D_{1113}\\ D_{2211}&D_{2222}&D_{2233}&D_{2212}&D_{2223}&D_{2213}\\ D_{3311}&D_{3322}&D_{3333}&D_{3312}&D_{3323}&D_{3313}\\ D_{1211}&D_{1222}&D_{1233}&D_{1212}&D_{1223}&D_{1213}\\ D_{2311}&D_{2322}&D_{2333}&D_{2312}&D_{2323}&D_{2313}\\ D_{1311}&D_{1322}&D_{1333}&D_{1312}&D_{1323}&D_{1313} \end{bmatrix}. \]
The matrix is used in the element stiffness
\[ \boldsymbol{K}^{e}_u = \int_{\Omega_e} \boldsymbol{B}_u^{\mathsf T} \boldsymbol{D}_e \boldsymbol{B}_u\,\mathrm{d}\Omega. \]
COMSOL performs this assembly through the Solid Mechanics material node after the constitutive matrix is supplied through C=C(D).
5.1 A simple principal-axis check
Suppose the principal directions coincide with the Cartesian axes and
\[ \varepsilon_1>0, \qquad \varepsilon_2<0, \qquad \varepsilon_3<0, \qquad \operatorname{tr}\boldsymbol{\varepsilon}>0. \]
Then the model gives approximately
\[ D_{1111}=(\lambda+2\mu)g(\phi), \]
\[ D_{1122}=\lambda g(\phi), \]
\[ D_{2222}=\lambda g(\phi)+2\mu. \]
The first two terms are degraded because the volumetric state is tensile. The compressive deviatoric contribution \(2\mu\) in \(D_{2222}\) is not degraded. This is the physical meaning of the different cc and dd factors in the text file.
6. Finite-element residual and matrix correspondence
Let \(\boldsymbol{d}\) be the nodal displacement vector. Then
\[ \boldsymbol{\varepsilon} = \boldsymbol{B}_u\boldsymbol{d}. \]
The element internal mechanical force is assembled from the total stress:
\[ \boldsymbol{F}^{e}_{u,\mathrm{int}} = \int_{\Omega_e} \boldsymbol{B}_u^{\mathsf T} \boldsymbol{S}\,\mathrm{d}\Omega. \]
The local tangent of the elastic part is
\[ \boldsymbol{K}^{e}_{uu} = \frac{\partial\boldsymbol{F}^{e}_{u,\mathrm{int}}} {\partial\boldsymbol{d}} = \int_{\Omega_e} \boldsymbol{B}_u^{\mathsf T} \boldsymbol{D}_e \boldsymbol{B}_u\,\mathrm{d}\Omega. \]
In a staggered displacement solve, \(p\) and \(\phi\) are fixed. Thus the pressure stress and \(g(\phi)\sigma_0\) affect the residual through \(\boldsymbol{S}\), while \(\boldsymbol{D}_e\) is the elastic \(u\)-\(u\) tangent.
The Newton form is
\[ \boldsymbol{K}^{(j)}_{uu} \Delta\boldsymbol{U} = \boldsymbol{R}^{(j)}_u. \]
For a fully monolithic solution, additional off-diagonal derivatives such as \(\partial R_u/\partial p\) and \(\partial R_u/\partial\phi\) may be included. They do not change the meaning of the local elastic tangent \(D\).
7. Complete staggered implementation
At time step \(n+1\), use the previous converged state as the initial guess:
Hold \(p^j\) and \(\phi^j\) fixed and solve Solid Mechanics for \(\boldsymbol{u}^{j+1}\) using the current \(D\) and the two external-stress contributions.
Compute principal strains, principal directions, \(Y\), and update
\[ H^{j+1}=\max(H^n,Y^{j+1}). \]
Hold \(H^{j+1}\) fixed and solve the Helmholtz phase-field equation for \(\phi^{j+1}\).
Update \(g(\phi)\), \(\alpha_B(\phi)\), permeability indicators, and the next mechanical tangent.
Check displacement, pressure, history, and phase-field changes. If they do not satisfy tolerance, repeat the staggered iteration.
After convergence, advance to the next time step.
The paper’s phase-field equation is solved as
\[ \mathbf{K}_{\phi}(H^{j+1})\boldsymbol{\Phi}^{j+1} = \mathbf{F}_{\phi}(H^{j+1}). \]
The phase field is therefore obtained from a matrix solve. \(Y\) and \(H\) enter the phase-field residual and matrix; they do not directly prescribe \(\phi\).
8. Implementation checks
- Confirm that
C=C(D)uses the same Voigt ordering as the matrix assembled fromg1111–g1313. - Confirm whether COMSOL uses tensorial shear strain or engineering shear strain in the material matrix. A factor-of-two mismatch produces incorrect shear stiffness.
- Keep
S_el=C:epsilon_elseparate from the pressure and initial-stress contributions. - Check that
cc=(1-k)*(1-phi)^2+kuses the solved phase field, not displacementu. - Check that
p_ref=0matches the pressure variable’s gauge/absolute convention. - Do not enter the initial stress both as an inelastic strain and as External Stress 2.
- Do not infer dynamic inertia merely from a Time Dependent study. The displayed equation is quasi-static unless the mass term is enabled.
- Verify that the
n1andm1definitions in the current text file are not being relied upon; they useuand are not referenced by the final tangent components. - Verify boundary tractions and displacement constraints separately from the domain stress equation.
Sources
- Zhou, Rabczuk, and Zhuang (2018), Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies, especially Sections 2.2–2.3 and 3.1–3.4, Eqs. (10)–(16), (20)–(26), and (30)–(41).
- 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, especially Sections 2.1–2.5 and the finite-element discretization section, Eqs. (13)–(18) and (32)–(39).
- Local COMSOL definitions:
Phase_Field_study_Definitions.txt.