Primary implementation checklist.
| Step | Module | Main objects |
|---|---|---|
| 1 | Definitions | parameters, functions, variables |
| 2 | Solid Mechanics | \(\boldsymbol{u}\), material, stress |
| 3 | Fluid Transport | \(p\), storage, flow |
| 4 | Phase-Field | \(\phi\), \(H\), \(g(\phi)\) |
| 5 | Coupling | \(\chi_r,\chi_f,\alpha,\boldsymbol{K}_{\mathrm{eff}}\) |
| 6 | Study / Solver | segregated or coupled iterations |
| 7 | Mesh / Diagnostics | \(h\), \(\ell_0\), \(\Delta t\), residuals |
1. Definitions
1.1 Parameters
| Group | Parameters |
|---|---|
| Phase field | \(\phi_0=0\), \(k\), \(G_c\), \(\ell_0\) |
| Transition | \(\mathrm{fai\_cr1}\), \(\mathrm{fai\_cr2}\), \(\mathrm{dfai\_cr}\) |
| Pressure | \(p_{\mathrm{ref}}\), \(\alpha_B\), \(S\) |
| Material | \(E\), \(\nu\), \(\lambda\), \(\mu\), \(\boldsymbol{k}_r\), \(\boldsymbol{k}_f\) |
| Stress | \(S_v\), \(S_h\), \(c_c\) |
\[ \mathrm{dfai\_cr} = \mathrm{fai\_cr2}-\mathrm{fai\_cr1}, \qquad \mathrm{fai\_cr1}<\mathrm{fai\_cr2}. \]
1.2 Analytic functions
Function argument: \(x\). Evaluate with \(x=\phi\).
\[ \mathrm{COE}_r(\phi) = \begin{cases} 1,&\phi\le\mathrm{fai\_cr1},\\ \dfrac{\mathrm{fai\_cr2}-\phi}{\mathrm{dfai\_cr}}, &\mathrm{fai\_cr1}<\phi<\mathrm{fai\_cr2},\\ 0,&\phi\ge\mathrm{fai\_cr2}, \end{cases} \]
\[ \mathrm{COE}_f(\phi) = \begin{cases} 0,&\phi\le\mathrm{fai\_cr1},\\ \dfrac{\phi-\mathrm{fai\_cr1}}{\mathrm{dfai\_cr}}, &\mathrm{fai\_cr1}<\phi<\mathrm{fai\_cr2},\\ 1,&\phi\ge\mathrm{fai\_cr2}. \end{cases} \]
\[ \mathrm{COE}_r(\phi)+\mathrm{COE}_f(\phi)=1. \]
1.3 Variables
\[ \chi_r=\mathrm{COE}_r(\phi), \qquad \chi_f=\mathrm{COE}_f(\phi). \]
\[ \alpha=\alpha_r\chi_r+\chi_f, \qquad \boldsymbol{K}_{\mathrm{eff}} = \boldsymbol{k}_r\chi_r+\boldsymbol{k}_f\chi_f. \]
\[ g(\phi)=(1-k)(1-\phi)^2+k, \qquad H_{\mathrm{new}}=\max(H_{\mathrm{old}},Y). \]
2. Solid Mechanics
2.1 Material matrix
Voigt order:
\[ (11,22,33,12,23,13). \]
Full matrix shown in the material node:
\[ \mathbf{D}=\begin{bmatrix} g_{1111} & g_{1122} & g_{1133} & g_{1112} & g_{1123} & g_{1113}\\ g_{1122} & g_{2222} & g_{2233} & g_{2212} & g_{2223} & g_{2213}\\ g_{1133} & g_{2233} & g_{3333} & g_{3312} & g_{3323} & g_{3313}\\ g_{1112} & g_{2212} & g_{3312} & g_{1212} & g_{1223} & g_{1213}\\ g_{1123} & g_{2223} & g_{3323} & g_{1223} & g_{2323} & g_{2313}\\ g_{1113} & g_{2213} & g_{3313} & g_{1213} & g_{2313} & g_{1313} \end{bmatrix}. \]
\[ \boldsymbol{\sigma}=\mathbf{D}\boldsymbol{\varepsilon}, \qquad \mathbf{D}\in\mathbb{R}^{6\times6}. \]
| Symbol | Meaning |
|---|---|
| \(g_{1111}\) | \(D_{1111}\) |
| \(g_{1122}\) | \(D_{1122}\) |
| \(g_{ijkl}\) | elasticity-tensor component |
| \(g(\phi)\) | phase-field degradation |
Use isotropic data for a homogeneous isotropic material. Anisotropic selection alone does not create fracture anisotropy.
2.2 Pore pressure
\[ \boldsymbol{\sigma}_p = -\alpha_B \left(p-p_{\mathrm{ref}}\right) \boldsymbol{I}. \]
\(p_{\mathrm{ref}}=0\) is a gauge-pressure convention; choose consistently for absolute pressure.
2.3 External stress
\[ \boldsymbol{S}_{\mathrm{ext}} = \begin{bmatrix} -S_hc_c & 0 & 0\\ 0 & -S_hc_c & 0\\ 0 & 0 & -S_nc_c \end{bmatrix}. \]
The screenshot shows units of \(\mathrm{N/m^2}\). The negative sign denotes compression under the tensile-positive convention. Here \(S_n\) is kept exactly as shown in the screenshot; verify whether it is an intentional normal-stress parameter or a naming typo for \(S_h\).
2.4 Boundaries
| Group | Pressure | Constraint |
|---|---|---|
| \(S_v\) | \(S_v\) | \(u_x=0\), \(u_y\) free |
| \(S_h\) | \(S_h\) | \(u_x\) free, \(u_y=0\) |
Check each boundary ID:
- external boundary vs internal interface;
- normal constraint vs tangential restraint;
- physical wall/symmetry/roller meaning.
3. Fluid Transport
\[ \rho S\frac{\partial p}{\partial t} - \nabla\cdot \left( \frac{\rho\boldsymbol{K}_{\mathrm{eff}}}{\mu_{\mathrm{eff}}} \nabla p \right) = q_m - \rho\alpha\chi_r \frac{\partial\varepsilon_{\mathrm{vol}}}{\partial t}. \]
\[ \varepsilon_{\mathrm{vol}} = \nabla\cdot\boldsymbol{u}. \]
Verify pressure diffusion and mass balance before coupling.
3.1 Storage and material checks
\[ S = \varepsilon_p c_r+ \frac{(\alpha_B-\varepsilon_p)(1-\alpha_B)}{K_{Vr}}. \]
| Fluid setting | Value |
|---|---|
| Type | incompressible |
| \(\rho_{\mathrm{ref}}\) | \(1000\,\mathrm{kg/m^3}\) |
| \(\mu\) | \(10^{-3}\,\mathrm{Pa\,s}\) |
| Check | \(S>0\) in solved domains |
The incompressible-fluid option does not remove porous-medium storage. A zero or negative \(S\) can make the initial pressure problem singular or ill-conditioned.
Check:
\[ \min(S)>0, \qquad \mathrm{COE}_r(\phi)+\mathrm{COE}_f(\phi)=1. \]
4. Phase-Field
4.1 Helmholtz Equation interface
Dependent field:
\[ \mathrm{comp1.phi}=\phi. \]
Interface form:
\[ \nabla\cdot(-c\nabla\phi)+a\phi=f. \]
Define
\[ A_{\mathrm{pf}}=\frac{2\ell_0(1-k)H}{G_c}. \]
| Coefficient | Unscaled form | Unit |
|---|---|---|
| \(c\) | \(\ell_0^2\) | \(\mathrm{m^2}\) |
| \(a\) | \(1+A_{\mathrm{pf}}\) | \(1\) |
| \(f\) | \(A_{\mathrm{pf}}\) | \(1\) |
\[ \nabla\cdot(-\ell_0^2\nabla\phi) +(1+A_{\mathrm{pf}})\phi =A_{\mathrm{pf}}. \]
Equivalent scaled form:
| Coefficient | Scaled form | Unit |
|---|---|---|
| \(c\) | \(1\) | \(1\) |
| \(a\) | \((1+A_{\mathrm{pf}})/\ell_0^2\) | \(\mathrm{m^{-2}}\) |
| \(f\) | \(A_{\mathrm{pf}}/\ell_0^2\) | \(\mathrm{m^{-2}}\) |
Do not mix the two scalings.
4.2 Discretization and boundary
| Setting | Choice |
|---|---|
| Element order | Linear Lagrange |
| Boundary | Zero Flux |
| Initial \(\phi\) | \(\phi_0=0\) |
| Initial \(\partial_t\phi\) | \(0\); no governing time term |
\[ \boldsymbol{n}\cdot(-\ell_0^2\nabla\phi)=0. \]
Check whether an all-boundaries selection includes internal partitions.
4.3 Driving energy variable
The user-defined variable \(\mathrm{fai\_p}\) represents \(Y\), not \(\phi\):
\[ \mathrm{fai\_p} = \frac{\lambda}{2}\,\mathrm{tra1}_p^2 + \mu\left(e_{1p}^2+e_{2p}^2+e_{3p}^2\right) - \left(S_h\varepsilon_x+S_v\varepsilon_y+S_n\varepsilon_z\right). \]
Tensile parts:
\[ e_{ip}=\max(\varepsilon_i,0), \qquad \mathrm{tra1}_p = \max(\varepsilon_1+\varepsilon_2+\varepsilon_3,0). \]
| Symbol | Meaning |
|---|---|
| \(\mathrm{fai\_p}\) | current crack-driving energy \(Y\) |
| \(\phi\) | phase field |
| \(H\) | irreversible history field |
The current variable definitions correspond to
$ e_{1p}=(,0), e_{2p}=(,0), e_{3p}=(,0). $
$ =++, _p=(,0). $
4.4 History field and initial crack
History update:
\[ H_{\mathrm{new}}=\max(H_{\mathrm{old}},\mathrm{fai\_p}). \]
For a distributed ODE/DAE history container:
- \(d_a=0\), \(e_a=0\) gives an algebraic residual;
- dode is a numerical container, not a new physical law;
- nojac(expr) removes expr from the symbolic Jacobian;
- nojac() alone does not enforce irreversibility.
Initial seed:
\[ H_0 = \frac{B\,G_c}{2\ell_0} \left( 1-\frac{2|y-D_d|}{\ell_0} \right). \]
Use a cutoff:
\[ H_{0,\mathrm{cut}} = \begin{cases} \dfrac{B\,G_c}{2\ell_0} \left(1-\dfrac{2|y-D_d|}{\ell_0}\right), & |y-D_d|<\ell_0/2,\\[8pt] 0,&\text{otherwise}. \end{cases} \]
Check the coordinate name \(y\); do not confuse it with the energy variable \(Y\).
For the supplied values:
\[ \ell_0=0.1\,\mathrm{m}, \qquad G_c=500\,\mathrm{J/m^2}, \qquad D_d=5\,\mathrm{m}, \qquad B\approx1000, \]
\[ H_{0,\max} = \frac{B\,G_c}{2\ell_0} \approx 2.5\times10^6\,\mathrm{J/m^3}. \]
Locally,
\[ \phi=\frac{A_{\mathrm{pf}}}{1+A_{\mathrm{pf}}}, \qquad A_{\mathrm{pf}}=\frac{2\ell_0(1-k)H}{G_c}. \]
At the seed center, \(A_{\mathrm{pf}}\approx1000\) and \(\phi\approx0.999\).
A seeded \(H\) is more persistent than an initial \(\phi\) alone. Exact \(\phi=1\) requires an additional constraint or prescribed initial field.
| Seed item | Current check |
|---|---|
| \(l_{\mathrm{ini}}\) | \(0.8\); geometry/selection control, not amplitude |
| domains | verify domains 11 and 12 are the crack band |
| coordinate | verify lowercase \(y\), not energy variable \(Y\) |
| cutoff | prevent negative \(H_0\) outside the intended band |
5. Coupling
\[ \phi \rightarrow \{\chi_r,\chi_f\} \rightarrow \{\alpha,\boldsymbol{K}_{\mathrm{eff}}\} \rightarrow \{p,\boldsymbol{u}\} \rightarrow Y \rightarrow H \rightarrow \phi. \]
| Function | Role |
|---|---|
| \(g(\phi)\) | tensile stiffness degradation |
| \(\chi_r(\phi)\) | reservoir properties |
| \(\chi_f(\phi)\) | fracture properties |
6. Study / Solver
Per time step:
- update \(\chi_r,\chi_f,\alpha,\boldsymbol{K}_{\mathrm{eff}}\);
- solve \(\boldsymbol{u}\);
- solve \(p\) or a coupled \((\boldsymbol{u},p)\) block;
- compute \(\mathrm{fai\_p}=Y\);
- update \(H\);
- solve \(\phi\);
- repeat to convergence.
| Setting | Recommendation |
|---|---|
| Solid Mechanics | Study controlled |
| Darcy’s Law | Study controlled |
| hzeq | time-dependent study, zero phase time coefficient |
| dode | distributed history container |
| Solver | Fully Coupled, Segregated, or outer staggered |
Study controlled follows the active study type. Time dependent forces the transient equation form. These are related but distinct settings.
The phase equation is elliptic at each time step; it is re-solved because \(H\) changes. Do not introduce an unintended \(\partial_t\phi\) term.
Use continuation/ramping if needed.
6.1 Verified time-dependent settings
| Setting | Current choice |
|---|---|
| Study | Time Dependent |
| Output | \(\operatorname{range}(0,0.05,10)\) |
| Relative tolerance | \(0.01\) |
| Time integrator | BDF |
| Absolute tolerance | manual scaled: \(0.001\) |
| Jacobian update | every iteration |
| Nonlinear strategy | Segregated |
\(\operatorname{range}(0,0.05,10)\) controls stored output times, not internal solver steps.
6.2 Segregated order
\[ (\boldsymbol{u},p) \rightarrow H \rightarrow \phi. \]
| Step | Variables | Notes |
|---|---|---|
| 1 | \(\boldsymbol{u},p\) | Direct solver; one inner iteration acceptable only if outer convergence is reached |
| 2 | \(H\) | previous accepted value + irreversible update |
| 3 | \(\phi\) | hzeq with \(H\) fixed |
Use BDF for the stiff quasi-static DAE system. Update the Jacobian every iteration because \(\phi\), \(H\), permeability, and degraded stiffness change strongly.
6.3 Scaling and direct solver
| Variable | Approximate scale |
|---|---|
| \(p\) | \(8.6\times10^8\) |
| \(\boldsymbol{u}\) | \(0.14\) |
| \(\phi\) | \(1\) |
| \(H\) | \(1.7\times10^6\) |
Manual scaled absolute tolerance is numerical, not physical. Compare \(0.001\) with \(0.0001\) for final verification.
The active direct solver is MUMPS. Reuse preordering and sparsity pattern when appropriate. Keep out-of-core off while memory is sufficient; switch to Automatic if allocation or swapping becomes problematic.
6.4 History node
| Node | Role |
|---|---|
| Previous Solution \(H\) | supplies \(H^n\) |
| Lower Limit \(H\ge0\) | non-negativity only |
| dode | distributed history container |
Previous Solution does not implement the maximum update by itself, and \(H\ge0\) does not enforce \(H^{n+1}\ge H^n\).
7. Mesh
| Item | Check |
|---|---|
| \(\ell_0\) | diffuse-crack width |
| \(h\) | resolves the crack band |
| \(\Delta t\) | time convergence |
| crack-tip region | local refinement |
Studies:
\[ h\text{-refinement}, \qquad \ell_0\text{-study}, \qquad \Delta t\text{-study}. \]
8. Diagnostics
| Quantity | Check |
|---|---|
| residual | convergence |
| mass balance | conservation |
| \(\max\phi\) | bounded, irreversible |
| fracture energy | no abnormal jump |
| \(\mathrm{COE}_r+\mathrm{COE}_f\) | equals \(1\) |
| boundary IDs | intended entities |
| \(\mathrm{comp1.phi}\) warning | disappears or is explained |
History checks:
- verify that \(H_{\mathrm{old}}\) is a previous-time or previous-iteration value;
- confirm that \(\max(H_{\mathrm{old}},\mathrm{fai\_p})\) is not using the current unknown \(H\) on both sides;
- confirm that nojac() changes only the Jacobian and does not provide the previous value automatically.
- compare pressure, displacement, \(\max(\phi)\), and crack path for scaled absolute tolerances \(0.001\) and \(0.0001\);
- check \(\min(S)\) and pointwise \(H_{\mathrm{new}}-H_{\mathrm{old}}\ge0\);
- distinguish small linear errors from convergence of the outer segregated iteration.
Record: version, model file, physics, variables, parameters, units, boundaries, mesh, solver, tolerances, and continuation settings.