Phase-Field Hydrofracture: Pressure Equation and COMSOL Mapping

This note establishes a one-to-one correspondence between the pressure-flow equations in the phase-field hydrofracture papers and the equations displayed by COMSOL’s Porous Matrix, Darcy-flow, and Mass Source nodes.

The central distinction is:

COMSOL S_p = local scalar pressure-storage coefficient
paper S   = local scalar pressure-storage coefficient
paper S_p = assembled finite-element pressure storage matrix

The symbol $S_p$ is therefore overloaded. In the COMSOL material panel it has units of $1/\mathrm{Pa}$ and is a scalar coefficient. In the 2019 paper Appendix, $\boldsymbol{S}_p$ is a matrix assembled from this scalar coefficient.

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}_s$	solid 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 subscript on $\boldsymbol{u}_s$ is used here to distinguish solid displacement from the Darcy velocity displayed by COMSOL. The pressure-flow velocity is denoted by $\boldsymbol{v}$ below.

The history field $H$ is an irreversible internal variable. It is not a fourth primary field, and it is not the Darcy velocity. Stress, strain, velocity, and force are derived quantities.

For one staggered iteration, the pressure-related dependency is

\[ \boldsymbol{u}_s,\phi \longrightarrow \varepsilon_{\mathrm{vol}},\chi_r,\chi_f,\alpha \longrightarrow S,\boldsymbol{K}_{p},\boldsymbol{q}_{p} \longrightarrow p. \]

The pressure solve is performed together with, or sequentially after, the displacement solve according to the chosen segregated implementation. The phase field is then updated and the pressure-dependent material coefficients are refreshed.

2. Pressure equation in the hydrofracture paper

2.1 Mass conservation and physical meaning

The 2020 hydrofracture paper writes the fluid mass balance as

\[ \begin{aligned} \rho S\frac{\partial p}{\partial t} +\nabla\cdot(\rho\boldsymbol{v}) &=q_m -\rho\alpha\chi_r \frac{\partial\varepsilon_{\mathrm{vol}}}{\partial t}. \end{aligned} \]

The terms have the following physical roles:

Term	Physical meaning
$\rho S\,\partial p/\partial t$	fluid storage caused by pressure change
$\nabla\cdot(\rho\boldsymbol{v})$	mass transport by Darcy flow
$q_m$	externally prescribed mass source or injection
$-\rho\alpha\chi_r\,\partial\varepsilon_{\mathrm{vol}}/\partial t$	fluid mass change caused by volumetric deformation of the porous skeleton

The volumetric strain is

\[ \varepsilon_{\mathrm{vol}} =\nabla\cdot\boldsymbol{u}_s. \]

The factor $\chi_r$ is the reservoir indicator. In the present implementation it is represented by $\mathrm{COE}_r(\phi)$.

The same governing equation is Eq. (23) of the 2019 multizone paper. In the 2020 paper it is Eq. (24).

2.2 Storage coefficient

The scalar storage coefficient is

\[ S = \varepsilon_p c + \frac{(\alpha-\varepsilon_p)(1-\alpha)}{K_{Vr}}. \]

This is Eq. (25) of the 2020 paper and Eq. (24) of the 2019 paper.

Here:

Symbol	Meaning
$\varepsilon_p$	porosity
$c$	fluid compressibility
$\alpha$	effective Biot coefficient
$K_{Vr}$	drained bulk modulus of the reservoir/calculation domain

The first term is fluid storage in the pore volume. The second term is the skeleton/grain-pore contribution associated with the drained bulk modulus.

The phase-dependent interpolation is

\[ \chi_r=\mathrm{COE}_r(\phi), \qquad \chi_f=\mathrm{COE}_f(\phi), \qquad \chi_r+\chi_f=1, \]

\[ c=c_r\chi_r+c_f\chi_f, \]

\[ \alpha=\alpha_r\chi_r+\alpha_f\chi_f =\alpha_r\chi_r+\chi_f, \qquad \alpha_f=1. \]

The 2019 paper also states that the fracture-domain porosity is $\varepsilon_p=1$. A phase-consistent interpolation is therefore

\[ \varepsilon_p =\varepsilon_{p,r}\chi_r+\chi_f. \]

The exact expression used by COMSOL’s dl.epsilon must be checked in the Porous Matrix node.

2.3 Darcy’s law

When gravity is neglected, the Darcy velocity is

\[ \boldsymbol{v} = -\frac{K_{\mathrm{eff}}}{\mu_{\mathrm{eff}}}\nabla p. \]

This is Eq. (26) of the 2020 paper and Eq. (25) of the 2019 paper when the gravity term is omitted.

The effective permeability and viscosity are interpolated as

\[ K_{\mathrm{eff}} =k_r\chi_r+k_f\chi_f, \]

\[ \mu_{\mathrm{eff}} =\mu_r\chi_r+\mu_f\chi_f. \]

The Darcy velocity $\boldsymbol{v}$ is not the solid displacement $\boldsymbol{u}_s$.

2.4 Combined pressure equation

Substituting Darcy’s law into the mass balance gives the pressure equation

\[ \begin{aligned} \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}. \end{aligned} \]

This is Eq. (27) of the 2020 paper and Eq. (26) of the 2019 paper for the no-gravity case.

The pressure Dirichlet and mass-flux boundary conditions are

\[ p=p_D\quad\text{on }\partial\Omega_D, \]

\[ -\boldsymbol{n}\cdot\boldsymbol{v}=M_N \quad\text{on }\partial\Omega_N. \]

3. The Porous Matrix equation shown by COMSOL

The supplied COMSOL Porous Matrix panel displays

rho*S_p*dp/dt + div(rho*v) = Q_m
v = -(kappa/mu)*grad(p)

In mathematical notation,

\[ \begin{aligned} \rho S_p\frac{\partial p}{\partial t} +\nabla\cdot(\rho\boldsymbol{v}) &=Q_m. \end{aligned} \]

\[ \boldsymbol{v} = -\frac{\kappa}{\mu}\nabla p. \]

After substitution,

\[ \begin{aligned} \rho S_p\frac{\partial p}{\partial t} -\nabla\cdot\left( \frac{\rho\kappa}{\mu}\nabla p \right) &=Q_m. \end{aligned} \]

This is the same structure as the paper equation. The notation mapping is

COMSOL	Paper	Meaning
S_p	$S$	local scalar storage coefficient
kappa	$K_{\mathrm{eff}}$	effective permeability
mu	$\mu_{\mathrm{eff}}$	effective dynamic viscosity
v	$\boldsymbol{v}$	Darcy velocity
Q_m	$q_m-\rho\alpha\chi_r\dot{\varepsilon}_{\mathrm{vol}}$	total mass source
p	$p$	pore pressure
rho	$\rho$	fluid density

The COMSOL label S_p is a scalar coefficient with units $1/\mathrm{Pa}$. It is not the assembled pressure matrix in the 2019 Appendix.

The finite-element matrix also denoted by $\boldsymbol{S}_p$ is assembled later as

\[ \boldsymbol{S}_p = \int_{\Omega} \boldsymbol{N}_p^{\mathsf T} \rho S \boldsymbol{N}_p\,\mathrm{d}\Omega. \]

Thus the scalar coefficient belongs in the Porous Matrix material settings, while the matrix is produced by finite-element assembly.

4. Mass Source 1: deformation coupling

The supplied Mass Source 1 node displays the generic mass-balance form

\[ \begin{aligned} \frac{\partial(\varepsilon_p\rho)}{\partial t} +\nabla\cdot(\rho\boldsymbol{v}) &=Q_m. \end{aligned} \]

with the user-defined source

Q_m = -dl.rho*alphaB*d(evol,TIME)*COE_r(phi)

The corresponding continuum expression is

\[ Q_m = -\rho\alpha\chi_r \frac{\partial\varepsilon_{\mathrm{vol}}}{\partial t}. \]

The variable mapping is

COMSOL variable	Model quantity
dl.rho	$\rho$
alphaB	$\alpha$
evol	$\varepsilon_{\mathrm{vol}}=\nabla\cdot\boldsymbol{u}_s$
COE_r(phi)	$\chi_r$
d(evol,TIME)	$\partial\varepsilon_{\mathrm{vol}}/\partial t$

This term represents fluid mass change caused by volumetric deformation of the porous solid. It is not an external injection. It acts as a source or sink depending on the sign of the volumetric strain rate and the sign convention used by Solid Mechanics.

In the present model,

evol = solid.ep1 + solid.ep2 + solid.ep3

when solid.ep1, solid.ep2, and solid.ep3 are the three principal strains. Therefore evol is the principal-strain trace

\[ \varepsilon_{\mathrm{vol}} =\varepsilon_1+\varepsilon_2+\varepsilon_3. \]

This deformation-coupling term is the second term on the right-hand side of the 2020 paper Eq. (27), and it is included in the 2019 paper Eq. (23) and Eq. (26).

Because $\mathrm{COE}_r(\phi)=\chi_r$, the term is active in the reservoir and transition regions and tends to zero in the fully fractured region when $\chi_r=0$.

5. Mass Source 2: deformation coupling plus injection

The supplied Mass Source 2 node uses

Q_m = -dl.rho*alphaB*d(evol,TIME)*COE_r(phi)
      + 10[kg/m^3/s]

In paper notation,

\[ Q_m = -\rho\alpha\chi_r \frac{\partial\varepsilon_{\mathrm{vol}}}{\partial t} + q_F, \]

with

\[ q_F =10\ \mathrm{kg}/(\mathrm{m}^3\mathrm{s}). \]

This source contains two physically different contributions:

$-\rho\alpha\chi_r\dot{\varepsilon}_{\mathrm{vol}}$ is the poroelastic deformation source/sink.
$q_F=10\ \mathrm{kg}/(\mathrm{m}^3\mathrm{s})$ is an externally prescribed fluid-injection mass source.

The Mass Source 2 node is selected only on domains 11 and 12. Therefore the constant injection term is restricted to these domains. It must not be applied to the entire porous medium unless the physical problem is intended to inject fluid everywhere.

The 2020 paper states that its numerical examples use a source term

\[ q_F=10\ \mathrm{kg}/(\mathrm{m}^3\mathrm{s}) \]

in the pre-existing notches. The general source $q_m$ and the deformation-coupling term are part of the flow equations in both the 2020 and 2019 papers.

The total source used by the pressure equation can therefore be written as

\[ Q_m =Q_{\mathrm{deformation}} + Q_{\mathrm{injection}}, \]

where

\[ Q_{\mathrm{deformation}} = -\rho\alpha\chi_r \dot{\varepsilon}_{\mathrm{vol}}, \]

and, only in the selected injection domains,

\[ Q_{\mathrm{injection}} =q_F. \]

6. Phase-dependent COMSOL properties

The indicator functions are evaluated from the solved phase field:

chi_r = COE_r(phi)
chi_f = COE_f(phi)

The project definitions use

alphaB = 0.05*COE_r(phi) + 1*COE_f(phi)
c_rf   = c_r*COE_r(phi) + c_f*COE_f(phi)
K      = lanta + 2/3*mu

The storage coefficient entered in the Porous Matrix node is

S_p = dl.epsilon*c_rf
      + (alphaB-dl.epsilon)*(1-alphaB)/K

The direct mapping is

\[ \mathrm{dl.epsilon}\longleftrightarrow\varepsilon_p, \qquad \mathrm{c\_rf}\longleftrightarrow c, \qquad \mathrm{alphaB}\longleftrightarrow\alpha, \qquad \mathrm{K}\longleftrightarrow K_{Vr}. \]

For the current parameters,

\[ E_0=60\ \mathrm{GPa}, \qquad \nu=0.3, \qquad K=\lambda+\frac{2}{3}\mu=50\ \mathrm{GPa}. \]

If the intact region uses $\varepsilon_p=0.05$ and $\alpha=0.05$, then

\[ S_{\mathrm{reservoir}} =0.05(10^{-8}) =5\times10^{-10}\ \mathrm{Pa}^{-1}. \]

If the fully fractured region uses $\varepsilon_p=1$ and $\alpha=1$, then

\[ S_{\mathrm{fracture}} =1(10^{-8}) =10^{-8}\ \mathrm{Pa}^{-1}. \]

These values require dl.epsilon to be assigned consistently in each region. A constant value 0.05 in the fracture would not reproduce the paper’s open-fracture porosity assumption.

7. Finite-element pressure residual and matrix correspondence

The 2019 paper’s Appendix gives the weak form of the pressure equation as

$$ _{} S p,

-_{} (p),

= _{} M_Np,S

_{} ( q_m-_r ) p,. $$

With pressure shape functions $p=\boldsymbol{N}_p\widehat{\boldsymbol{p}}$, the discrete residual is

\[ \boldsymbol{R}^p = \boldsymbol{S}_p \frac{\mathrm{d}\widehat{\boldsymbol{p}}}{\mathrm{d}t} + \boldsymbol{H}_p\widehat{\boldsymbol{p}} - \boldsymbol{q}^{\mathrm{ext}}_p =\boldsymbol{0}. \]

The pressure storage and permeability matrices are

\[ \boldsymbol{S}_p = \int_{\Omega} \boldsymbol{N}_p^{\mathsf T} \rho S \boldsymbol{N}_p\,\mathrm{d}\Omega, \]

\[ \boldsymbol{H}_p = \int_{\Omega} \boldsymbol{B}_p^{\mathsf T} \frac{\rho K_{\mathrm{eff}}}{\mu_{\mathrm{eff}}} \boldsymbol{B}_p\,\mathrm{d}\Omega. \]

After backward-Euler or BDF-type time discretization, the pressure block has the form

\[ \left( \frac{\boldsymbol{S}_p}{\Delta t} + \boldsymbol{H}_p \right) \widehat{\boldsymbol{p}}^{\,n+1} = \boldsymbol{q}^{\mathrm{ext}}_p + \frac{\boldsymbol{S}_p}{\Delta t} \widehat{\boldsymbol{p}}^{\,n} + \boldsymbol{q}^{\mathrm{coupling}}_p. \]

The source terms from Mass Source 1 and Mass Source 2 enter $\boldsymbol{q}^{\mathrm{ext}}_p$ or the equivalent COMSOL source residual. The matrices $\boldsymbol{S}_p$ and $\boldsymbol{H}_p$ are generated by finite-element assembly; they are not entered manually as 6-by-6 material matrices.

8. Pressure coupling in the staggered implementation

At time step $n+1$, the pressure-related staggered procedure is

Hold $\phi^j$ and $H^j$ fixed, update the phase-dependent hydraulic properties, and solve displacement/pressure for $(\boldsymbol{u}_s,p)^{j+1}$.
Compute $\varepsilon_{\mathrm{vol}}^{j+1}$ and its time derivative for the deformation source.
Use the intended source expression: deformation-only in porous domains and deformation-plus-injection only in domains 11 and 12.
Compute the crack-driving energy and update

\[ H^{j+1}=\max(H^n,Y^{j+1}). \]
Hold $H^{j+1}$ fixed and solve the phase-field equation for $\phi^{j+1}$.
Update $\chi_r$, $\chi_f$, $\alpha_B$, $K_{\mathrm{eff}}$, $S_p$, and the mechanical tangent.
Check displacement, pressure, history, and phase-field changes. If the tolerance is not satisfied, repeat the staggered iteration.

The pressure equation is therefore one block of the coupled three-field problem. The pressure source depends on the solid deformation, while the storage, Biot, and permeability coefficients depend on the phase field.

9. Implementation checks

Keep $\phi=0$ intact and $\phi=1$ fully broken.
Distinguish the scalar COMSOL S_p from the assembled finite-element matrix $\boldsymbol{S}_p$.
Distinguish Darcy velocity $\boldsymbol{v}$ from solid displacement $\boldsymbol{u}_s$.
Use evol=solid.ep1+solid.ep2+solid.ep3 only if these variables are the principal strains of the Solid Mechanics interface.
Check the sign of d(evol,TIME) against the Solid Mechanics tensile-positive/compressive-negative convention.
Use COE_r(phi) as $\chi_r$ and do not pass a geometric coordinate to the indicator function.
Keep -dl.rhoalphaBd(evol,TIME)*COE_r(phi) as the deformation coupling term, not as the injection term.
Restrict +10[kg/m^3/s] to the intended injection domains 11 and 12.
Confirm that the source is a mass source. Do not mix $\mathrm{kg}/(\mathrm{m}^3\mathrm{s})$ with a volumetric source without multiplying or dividing by $\rho$ consistently.
Confirm that fluid density is not multiplied twice by the COMSOL interface and the user-defined expression.
Confirm that dl.epsilon uses the intended reservoir-to-fracture interpolation.
Do not remove the pressure-storage term merely because the fluid model is labelled incompressible.
Use $K=\lambda+2\mu/3$ as the drained bulk modulus in the storage coefficient, not as the permeability.
Keep the pressure reference and pressure boundary conditions consistent with gauge or absolute pressure.
In the segregated solve, update the source and hydraulic coefficients at the same stage used by the active solver sequence.

Sources

Zhou, S., Zhuang, X., and Rabczuk, T. (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.5 and 3, Eqs. (23)–(28). The paper provides the storage coefficient, Darcy law, combined pressure equation, pressure boundary conditions, and the source value $q_F=10\ \mathrm{kg}/(\mathrm{m}^3\mathrm{s})$ used in the numerical examples.
Li, K. and Zhou, S. (2019), Numerical investigation of multizone hydraulic fracture propagation in porous media: New insights from a phase field method, Section 2.4, Eqs. (23)–(35), and Appendix, Eqs. (38)–(57). The paper provides the flow equation, storage coefficient, indicator functions, Darcy law, weak form, pressure residual, and matrices $\boldsymbol{S}_p$ and $\boldsymbol{H}_p$.
Local COMSOL definitions: Phase_Field_study_Definitions.txt and Phase_Field_study_Parameters.txt.