相场断裂¶
把脆性断裂看作能量最小化问题。损伤场 \(\alpha \in [0, 1]\) 把裂纹弥散在长度 \(\ell\) 的范围内，于是裂纹在哪里萌生、沿哪条路径扩展、何时失稳扩展，都由同一个能量的最小化决定 [Francfort1998] [Bourdin2000]。脚本 examples/solid/phase_field_fracture/phase_field_fracture.py 在平面应变下运行单边缺口拉伸（SENT）基准算例 [Miehe2010]。
与 塑性（带各向同性硬化的 J2） 相比，新增的要素是 不等式约束：损伤不可愈合，\(\alpha_n \ge \alpha_{n-1}\)。因此每个载荷步都是一个带界约束的最小化问题。脚本用交替最小化求解它，两个子问题都使用 torch-sla 中的 vi_solve——一种带上下界的牛顿求解。
模型¶
位移 \(\mathbf{u}\) 与损伤场 \(\alpha\) 的能量为
其中 \(\boldsymbol{\varepsilon} = \tfrac12(\nabla\mathbf{u} + \nabla\mathbf{u}^T)\)，\(g_0\) 是一个小的残余刚度。
裂纹模型。 AT1 取 \(w(\alpha) = \alpha\)、\(c_w = 8/3\)，它有一个弹性阶段：应力低于临界值时不产生损伤 [Pham2011]。AT2 取 \(w(\alpha) = \alpha^2\)、\(c_w = 2\)，在任意载荷下都会产生损伤。用
--model选择。能量分解 [Amor2009]。只有受拉的体积部分和偏量部分会被退化，因此压缩不会产生裂纹：
\[\psi_D = \tfrac{\kappa}{2} \langle \operatorname{tr}\boldsymbol{\varepsilon} \rangle_+^2 + \mu\, |\operatorname{dev}\boldsymbol{\varepsilon}|^2, \qquad \psi_R = \tfrac{\kappa}{2} \langle \operatorname{tr}\boldsymbol{\varepsilon} \rangle_-^2.\]
脚本只在 PhaseFieldModel 中写一次这两个能量密度。energy() 对它们积分，autograd 再把积分转换为残差 \(\partial\mathcal{E}/\partial\mathbf{u}\) 和 \(\partial\mathcal{E}/\partial\alpha\)。
以界约束表达不可逆性¶
在载荷步 \(n\)，损伤在 \([\alpha_{n-1}, 1]\) 中求解。在这些界约束下的极小点逐节点满足
要么残差为零，要么该节点停在一个界上，且残差把它推向界外。这是一个变分不等式，正是 torch_sla.vi_solve(residual, x0, lower=..., upper=..., jacobian_fn=...) 所求解的问题。所用方法是降维空间的有效集牛顿法，即 PETSc 的 SNESVINEWTONRSLS 所用的算法 [Benson2006]。每个牛顿步先对不在界上的节点求解，再更新哪些节点停在界上。
给定值同样是界约束：lower == upper 即固定该分量。位移的狄利克雷边界条件就是这样施加的，因此位移求解不需要 Condenser。预制裂纹也可以用同样的方式给定，即在一条节点线上保持 \(\alpha = 1\)（alpha_fixed 参数）。下面的基准算例改为在网格中开一条狭缝，使缺口面从一开始就是自由的。
交替最小化¶
每个载荷步在两个场之间交替求解，直到二者相容：
位移，固定 \(\alpha\)。\(\mathcal{E}\) 关于 \(\mathbf{u}\) 是凸的，但只是分片二次的，因为能量分解随 \(\operatorname{tr}\boldsymbol{\varepsilon}\) 的符号切换。带线搜索的牛顿法只需几步。
损伤，固定 \(\mathbf{u}\)。对 AT1 和 AT2，\(\mathcal{E}\) 关于 \(\alpha\) 是二次的，因此这是一个带界约束的二次规划。它的 Hessian 不依赖于 \(\alpha\)，只需装配一次。
停止：用新的损伤计算出的、未约束自由度上的位移残差小于反力的 \(10^{-4}\) 倍时停止。
def solve_step(self, u_bc):
"""One load step with the Dirichlet values ``u_bc`` (full dof vector);
returns the number of staggered iterations."""
inf = torch.full_like(self.u, math.inf)
u_lower = torch.where(self.u_fixed, u_bc, -inf)
u_upper = torch.where(self.u_fixed, u_bc, inf)
a_lower = torch.where(self.alpha_fixed, self.alpha_value, self.alpha) # irreversibility
a_upper = torch.where(self.alpha_fixed, self.alpha_value, torch.ones_like(self.alpha))
scale = 1e-12
for it in range(1, self.max_staggered + 1):
alpha = self.alpha
self.u = vi_solve(lambda u: self.residual_u(u, alpha), self.u, lower=u_lower, upper=u_upper,
jacobian_fn=lambda u: self.tangent_u(u, alpha), line_search=True,
tol=1e-9, atol=1e-11, linear_solver=self.linear_solver)
u = self.u
hessian = self.hessian_alpha(u) # constant during the alpha solve
alpha_new = vi_solve(lambda a: self.residual_alpha(a, u), alpha, lower=a_lower, upper=a_upper,
jacobian_fn=lambda a: hessian, tol=1e-9, atol=1e-11,
linear_solver=self.linear_solver)
r = self.residual_u(u, alpha_new)
free = float(torch.linalg.vector_norm(r[~self.u_fixed]))
scale = max(scale, float(torch.linalg.vector_norm(r[self.u_fixed])))
if free <= self.tol * scale:
self.alpha = alpha_new
return it
if self.relaxation != 1.0:
alpha_new = (alpha + self.relaxation * (alpha_new - alpha)).clamp(a_lower, a_upper)
self.alpha = alpha_new
raise RuntimeError(f"staggered scheme did not converge in {self.max_staggered} iterations")
两个雅可比矩阵都是双线性形式。与 TensorMesh 的所有形式一样，forward 描述一个求积点和一对基函数。点数据 disp 和 alpha 以 graddisp 和 alpha 的形式传入：
class DisplacementTangent(ElementAssembler):
"""d2E/du2 with the split's tangent C = g C_0 + (1 - g) kappa [tr eps < 0] I x I."""
def __post_init__(self, model: PhaseFieldModel):
self.model = model
def forward(self, gradu, gradv, graddisp, alpha):
m = self.model
g = m.g(alpha)
compressed = (graddisp.diagonal().sum() < 0).to(gradu.dtype)
lam = g * m.lam + (1 - g) * m.kappa * compressed
I = torch.eye(gradu.shape[0], dtype=gradu.dtype, device=gradu.device)
return lam * torch.outer(gradu, gradv) + g * m.mu * ((gradu * gradv).sum() * I + torch.outer(gradv, gradu))
class PhaseFieldHessian(ElementAssembler):
"""d2E/dalpha2, which does not depend on alpha for AT1 and AT2 (g'' = 2)."""
def __post_init__(self, model: PhaseFieldModel):
self.model = model
def forward(self, u, v, gradu, gradv, graddisp):
m = self.model
reaction = 2 * m.psi_D(m.strain(graddisp)) + m.Gc / (m.cw * m.ell) * m.ddw
return reaction * u * v + 2 * m.Gc * m.ell / m.cw * (gradu * gradv).sum()
单边缺口拉伸¶
几何：单位正方形 \([0, 1]^2\) mm，沿 \(y = 0.5\) 从左边缘到中心开有一条狭缝。狭缝上方的单元拥有其节点的独立副本，因此两侧缝面可以自由分开。
材料：\(E = 210\) kN/mm²，\(\nu = 0.3\)，\(G_c = 2.7 \times 10^{-3}\) kN/mm，\(\ell = 0.015\) mm，AT1，\(g_0 = 10^{-6}\)。
加载：底边固支；顶边 \(u_x = 0\)、\(u_y = \delta\)，增量为 \(10^{-4}\) mm。
网格：用 NumPy 生成的对角线交错（union-jack）三角形网格。相邻单元的对角线像棋盘一样交替，因此网格关于裂纹线镜像对称，不偏向任何裂纹方向。在裂纹路径附近的条带内 \(h = \ell / 2\)，其上下更粗（约 6,400 个节点）。
图 63 phase_field_fracture.py 的输出。左：峰值载荷时狭缝尖端（白线）附近的损伤。中：断裂后的试样，位移经过放大。右：顶边上的反力。¶
在 \(\delta \approx 3.9 \times 10^{-3}\) mm（\(F \approx 0.55\) kN）之前，响应是线弹性的，此时狭缝尖端才首次出现损伤：AT1 有弹性阶段。随后尖端的过程区不断长大，而载荷继续上升，在 \(\delta = 6.2 \times 10^{-3}\) mm 时达到 0.87 kN 的峰值。在下一个载荷步中，裂纹笔直地贯穿韧带，反力降为零。在位移控制下这种扩展是不稳定的，因此在单个载荷步内完成。
弹性阶段每步只需一到两次交错迭代，过程区增长的步最多约三十次。裂纹扩展的那一步需要 138 次，整个计算共约 320 次。使用 --model at2 时，损伤从第一步就开始增长，响应在峰值之前就已软化（\(\delta = 5.9 \times 10^{-3}\) mm 时为 0.73 kN），但裂纹同样在一个载荷步内贯穿。
验证¶
tests/assemble/test_phase_field_fracture_example.py 检查以下内容：
Hessian：将装配得到的 Hessian 与 autograd 残差的中心差分进行比较。
均匀解：AT1 下单轴拉伸的条带在临界应变 \(\varepsilon_c = \sqrt{3 G_c (1 - \nu^2) / (8 E \ell)}\) 以下保持无损伤。超过临界应变后，损伤为 \(\alpha = 1 - (\varepsilon_c / \varepsilon)^2\)，反力为 \(g(\alpha)\, E \varepsilon / (1 - \nu^2)\) 乘以宽度。
不可逆性：损伤从不减小，卸载步之后也是如此。
给定损伤：通过
alpha_fixed保持在 \(\alpha = 1\) 的节点严格保持该值。粗网格 SENT 算例：裂纹从狭缝尖端出发，沿 \(y = 0.5\) 笔直扩展到右边缘。它耗散的能量为 \(G_c\) 乘以裂纹长度，相差一个与网格有关的有效韧性系数（见下文）。
运行它¶
cd examples/solid/phase_field_fracture
python phase_field_fracture.py # AT1, ell = 0.015 mm
python phase_field_fracture.py --model at2
python phase_field_fracture.py --ell 0.04 --h 0.02 # coarse and quick
默认算例在一个 CPU 节点上约需一分钟（AT2 约一分半），粗网格算例只需几秒。--device cuda 在 GPU 上装配，此时牛顿方程组在 torch-sla 支持 cuDSS 时用 cuDSS 求解，否则在主机上用 SciPy 求解。
说明¶
分辨率。 \(h\) 必须能分辨 \(\ell\)（这里 \(h = \ell/2\)）。此时耗散能会比 \(G_c\) 乘以裂纹长度多出一个与网格有关的系数，因为完全断裂的带至少有一个单元宽。这里 0.5 mm 的裂纹耗散 \(0.61\, G_c\)，系数为 1.22。
长度尺度是材料参数。 对 AT1，\(\ell\) 还决定了强度 \(\sigma_c = \sqrt{3 G_c E' / (8 \ell)}\)，平面应变下 \(E' = E / (1 - \nu^2)\) [Tanne2018]。改变 \(\ell\) 不仅改变裂纹的光滑程度，也改变裂纹何时萌生。
求解器。 交替最小化稳健，但在裂纹扩展时收敛较慢。对损伤更新做超松弛 [Farrell2017] 有帮助：这里
--relaxation 1.6得到相同的结果，迭代次数少 15–20 %；当裂纹在许多载荷步中逐步扩展时节省更多。复合牛顿法和整体（monolithic）牛顿法所需的迭代还可以更少，它们同样用vi_solve施加界约束。
下一步¶
塑性（带各向同性硬化的 J2）：每步更新一次的历史变量，这是让固体模型具有路径相关性的另一种方式。
G. A. Francfort and J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, Journal of the Mechanics and Physics of Solids, 46(8):1319–1342, 1998.
B. Bourdin, G. A. Francfort and J.-J. Marigo, Numerical experiments in revisited brittle fracture, Journal of the Mechanics and Physics of Solids, 48(4):797–826, 2000.
C. Miehe, M. Hofacker and F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Computer Methods in Applied Mechanics and Engineering, 199(45–48):2765–2778, 2010.
K. Pham, H. Amor, J.-J. Marigo and C. Maurini, Gradient damage models and their use to approximate brittle fracture, International Journal of Damage Mechanics, 20(4):618–652, 2011.
H. Amor, J.-J. Marigo and C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments, Journal of the Mechanics and Physics of Solids, 57(8):1209–1229, 2009.
S. J. Benson and T. S. Munson, Flexible complementarity solvers for large-scale applications, Optimization Methods and Software, 21(1):155–168, 2006.
E. Tanné, T. Li, B. Bourdin, J.-J. Marigo and C. Maurini, Crack nucleation in variational phase-field models of brittle fracture, Journal of the Mechanics and Physics of Solids, 110:80–99, 2018.
P. E. Farrell and C. Maurini, Linear and nonlinear solvers for variational phase-field models of brittle fracture, International Journal for Numerical Methods in Engineering, 109(5):648–667, 2017.