传统脆性断裂相场模型的三维UEL理论及代码

1 引言

本部分介绍来自于《断裂相场法》书籍。

“1998年Francfort和Marigo根据Griffith脆性断裂理论，提出了一种断裂力学变分原理，他们以结构内可能的位移场和裂纹面作为自变量，将变形能与断裂面之和定义为结构总能量，并且认为真实的位移场与裂纹面使得该总能量最小。然而在数值模拟中将离散的裂纹面作为未知量来求解是非常困难的。因此2000年Bourdin等提出了一种相场模型，其中引入了一个连续的标量场，即相场，来近似地描述裂纹。相场值为1和0分别代表材料完全破坏和完好两种极限状态，而它们之间的值代表了一种损伤状态，并且裂纹的弥散程度由相场特征宽度来控制，其值越大弥散宽度越大，反之则越小。然后通过一个与相场相关的裂纹面密度泛函来重构结构内的断裂能，并将因损伤而退化的变形能与重构的断裂能代入Francfort-Marigo变分原理就得到了相场模型的基本列式。相场模型中的自变量为两个连续变化的场，即位移场和相场，因此它可以很方便地由不同数值方法实现。直观来看，相场模型将一个结构内裂纹萌生与演化问题，转化为了一个多场耦合情况下求最小能量的优化问题，因此它可以用于直接求解（例如分叉、交叉、融合、扭结等）复杂断裂问题，而不需要额外的裂纹路径追踪方法。”

2 理论

将系统的总势能表示为如下两项：

式中第一项能量为：

考虑损伤带来的退化，弹性能的表达式为：

式中

k为一个小值，用于防止数值不稳定现象。另一项断裂能为：

因此代入具体表达式可将系统总势能表达为：

对上述能量进行一阶变分可得：

即可得弱形式方程为：

具体外力虚功为：

式中本构方程为：

该弱形式方程是后续推导有限元方程的基础。同时，通过弱形式方程也可推导得到强形式的控制方程，即位移场和相场的控制方程。对上述弱形式进行分部积分可得：

因次位移场和相场的强形式控制方程为：

以及相应的边界条件为：

3 有限元离散

为推导有限元离散方程，对位移场和相场控制方程的弱形式进行处理：

对位移场和相场进行插值可得：

m指单元节点的个数。因此相应的梯度场可以插值为：

B矩阵的是由形函数对物理坐标的导数组成的。同理有：

代入到弱形式方程中可得残值方程；

使用牛顿迭代法求解上述非线性系统。更新格式为:

刚度矩阵为：

为了保证损伤不能愈合，即：

需要做出一些修改，即取历史上最大的弹性应变能，即：

4 代码

《断裂相场法》书中提供了传统脆性断裂相场模型的二维UEL代码。本文将其拓展为三维情况。

UEL需要更新单元的刚度矩阵和右端项，公式在理论部分已详细给出。

5 测试

5.1 一个单元拉伸破坏

对一个单元施加的边界条件为，x=0面约束所有位移自由度，即u=0,v=0,w=0;在x=1面进行位移加载，即u=0.1，v=0,w=0。单个单元拉伸破坏时具有解析解。因为相场在单元内是均匀的，因此相场梯度为0，因此可以直接求解相场控制方程得到：

因此可以得到真实应力的表达式为：

数值计算得到的应力应变曲线与解析解的对比结果如下：

取单元上一个积分点，绘制其相场值随加载时间的变化曲线如下：

5.2 单边裂纹拉伸

对于单边裂纹板的拉伸案例，我们选取两种相场特征裂纹宽度的情况，来展示特征宽度对于结果的影响。

第一个取l=0.015，相场分布图如下：

第二个取l=0.0075场分布图如下：

6 参考文献

[1]. 胡小飞, 张鹏, 姚伟岸. 断裂相场法. 北京: 科学出版社; 2022.

[2] Kristensen PK, Martínez-Pañeda E. Phase field fracture modelling using quasi-Newton methods and a new adaptive step scheme. Theoretical and applied fracture mechanics. 2020;107:102446.