一种改进的相场模型裂纹模拟方法

2022-05-13 05:04张晨露程勇刚李通盛

土木工程与管理学报 2022年2期

高宇，张晨露，岳强，王桥，周伟，程勇刚，李通盛

(1. 大唐宣威水电开发有限公司，云南宣威 655400； 2. 武汉大学水资源与水电工程科学国家重点实验室，湖北武汉 430072)

在实际工程应用中，对断裂的数值计算有较高的要求，目前研究比较多的方法包括内聚力模型(Cohesive Zone Model, CZM)[2]、连续断裂力学法(Continuum Damage Mechanics, CDM)[3]、扩展有限元法(Extended Finite Element Method, XFEM)[4]。内聚力模型能够模拟材料中裂纹萌生、演化直至完全失效的连续过程[5]，但是由于需要事先通过内聚力单元来定义断裂路径，因此内聚力模型很难模拟复杂受力情况下和含有不同方向的多裂纹的裂纹扩展规律。连续断裂力学法是另一种常用于有限元法中求解断裂力学问题的方法，为了得到较高的计算精度并有效预测裂纹的扩展规律，该类型的方法需要在裂纹尖端处加密网格。扩展有限元法是传统有限元法的改进，扩展有限元法虽然模拟裂纹扩展问题较为有效，但同样需要将裂纹面处理为不连续面，在模拟多裂纹交叉时存在一定难度。除了以上研究较多的方法，还有许多其他方法被应用于裂纹扩展问题中。例如Tsang等[6]采用离散裂隙网络模型(Discrete Fracture Network Model, DFN)能够综合考虑裂纹中的断裂力学和流体流动，但该方法计算精度较差，并且对裂纹的角度有较大的限制。Zhang等[7]采用了离散元法(Discrete Element Method, DEM)对散体结构中的水力劈裂过程进行了模拟，该方法模拟真实问题时计算量过大，无法实时得到计算结果。

近年来，基于相场模型(Phase-Field Model, PFM)[8]的裂纹模拟方法得到了广泛研究。传统的断裂力学中，存在着裂纹面和弹性体之间的不连续性问题，甚至还会有裂纹尖端的奇异性问题[9]。相场模型通过引入一个相场参数或函数来描述被近似为一个裂纹场(Crack Field)的裂纹面，将完全破坏的裂纹面与没有破坏的弹性体之间光滑化，消除了不连续性[10]。该方法的优势在于可以非常方便地模拟复杂裂纹的扩展、贯通、交叉和分叉等行为，能够预测裂纹开裂和扩展而不需要额外的判据[11]。一些不同的相场模型分别由不同的学者在物理界和力学界提出，Miehe等[12]提出了一种热力学一致性相场模型来模拟断裂问题，该模型基于连续力学和热力学。本文应用的相场模型主要是基于力学界中的脆性断裂变分方程[13]和正则化方程[14]，该方程是Griffith断裂理论[15]的延伸。对于相场模型的求解，目前大部分的方法均是基于有限元法。

从上文介绍中可以发现，相场模型是一种求解断裂问题较为理想的方法，该方法可以非常方便地模拟多裂纹的扩展、贯通、开叉等规律，并相对容易由二维问题扩展到三维问题，是一种非常有前景的求解断裂问题的方法。

1 准静态问题的相场模型

1.1 断裂变分准则

相场模型是基于Griffith断裂理论的变分准则，认为弹性体的总势能包括弹性体能和裂纹表面能[16]。

(1)

(2)

式中：Ω为弹性体的求解域；Γc⊂Ω为弹性体内的裂隙面(图1)；u为位移向量；ε为应变张量；Gc为材料的断裂韧度；ψ0为弹性体的能量密度；λ，μ为Lame常数。应变张量和位移向量的关系为：

(3)

式中：∇为梯度符号。

图1 弹性体中的裂纹面

1.2 正则化方程

式(1)中包含对裂纹面的积分，在时间数值计算中存在两大困难，一是位移场在裂纹处存在不连续性，二是从能量方程中无法直接得到最佳的裂纹扩展路径。为了解决这一问题，一些学者提出了相应的正则化方程[14]：

(4)

式中：s为用来表征裂纹场的场变量，其值在0～1之间变化，如果是0，表示完全破坏，如果是1，表示没有破坏；参数l0>0用来控制s的过渡区域的宽度(图2)；η为大于0且远远小于1的无量纲参数，用来模拟断裂时，即s=0，产生的残余刚度。

图2 s表示的裂纹用相场

于是应力张量可以定义为：

=(s2+η)(λtr(ε)I+2με)

(5)

1.3 控制方程

为了得到相场模型的最终控制方程，可以通过最小势能原理建立的泛函推导出相场模型的弱积分形式：

(6)

(7)

式(6)(7)分别对应弹性力学部分和相场部分，分别对式(6)左端域积分和式(7)的第一项域积分进行分部积分，可得：

(8)

(9)

式中：n是边界∂Ω上的外法相向量。由于式(8)对任意δu都成立，于是式(8)(9)对应的控制方程分别为：

divσ=0 onΩ

(10)

(11)

相应的边界条件分别为：

(12)

(13)

式中：∂Ωt为面力已知的边界(图3)；∂Ωs′为相场模型的Neumann边界(图4)。

图3 弹性部分的边界条件

图4 相场部分的边界条件

1.4 历史场变量

(14)

εi和ni分别为主应变和主应变方向，则

(15)

其中尖括号的含义为：

〈x〉±=(x± |x|)/2

(16)

式(5)便可以写为：

(17)

相场方程可改写为：

(18)

式中：ћ为历史常变量，其定义为：

(19)

最终的控制方程和边界条件可以归纳为：

(20)

弹性部分的边界条件和普通无裂纹的弹性力学问题的边界条件一样(图3)，∂Ωu为位移已知边界，相场部分的边界条件为在整个求解域的边界∂Ω上为Neumann边界条件，即在∂Ω=∂Ωs′上有sini=0，且在裂纹面Γc上s=0(图4)。

2 分步解耦算法

在采用有限元法数值计算过程中，对于每一加载时间步，可以直接对式(20)进行联立求解，为避免直接求解非线性方程组，本文采用解耦的方法分别求解弹性方程和相场方程。对于每一时间步，将式(20)解耦为弹性方程和相场方程：

(21)

(22)

式中：σn，un，sn，ћn分别为第n时间步的应力张量值、位移场、相场和历史场变量。分步解耦算法的具体流程图见图5。

图5 分步解耦算法的具体流程图

3 应变-相场法预设裂纹的改进方法

本文拟采用有限元法求解相场模型，位移场通过有限元法形函数近似可得

(23)

以二维情况为例，有

u=[u1(x)u2(x)]T

(24)

(25)

(26)

能量密度函数ψ0可以写为：

(27)

当模型有初始裂纹时，一般的方法是采用如下的应变-相场法来预设初始裂纹[1]：

(28)

式中：L为裂纹曲线或裂纹面；d(x,L)为点x到L的最短距离；H0为初始历史场变量的值。

但参照已有的数值算例表明，采用该方法的计算结果与几何不连续法在某些情况下差别较大，主要原因为几何不连续法可以认为是将裂纹面上下单元弱化，裂纹面上下单元完全分离。本文在此基础上，提出了改进的应变-相场法来预设裂纹，期望得到与几何不连续法预设裂纹相近似的效果，其基本思想是假设在裂纹面附近的单元一直弱化，无论其应变状态如何。

其中，对于平面应变问题有

(29)

在改进的应变 - 相场法中，首先通过式(28)来计算裂纹面附近积分点的初始历史场变量，且在所有计算过程中，采用下式来计算弱化后的D矩阵

(30)

于是式(21)的弱积分形式可以写为

(31)

最终可得的弹性方程的有限元离散格式为：

(32)

(33)

(34)

4 数值算例

4.1 含单边裂纹的受拉方板

考虑如图6所示包含初始裂纹的受拉方板，边长为1 mm，划分为如图7所示的54040个三角形单元，在预计裂纹扩展路径处局部加密，该处网格大小约为0.001 mm。设该问题为平面应变问题，弹性模量为210 kN/mm2，泊松比为0.3，断裂韧度为Gc=2.7×10-3kN/mm2。采用位移加载法，前500步位移增量为Δu=1×10-5mm，之后位移增量调整为Δu=1×10-6mm直到完全破坏。在计算过程中，取l0=0.004 mm，η=0。

首先采用几何不连续法来预设初始裂纹，即在裂纹处将上下单元在几何上分离，计算得到的裂纹扩展过程如图8所示。

图6 受拉方板/mm

图7 受拉方板网格划分

图8 受拉方板裂纹扩展路径(裂纹几何建模)

若采用应变 - 相场法来预设初始裂纹，得到的裂纹扩展路径如图9所示。可以发现，采用该方法来预设裂纹时，模型会早一些开裂，主要原因在于对于该模型，采用应变 - 相场法预设裂纹时，把裂纹附近在l0/2范围内的单元均弱化了，而采用几何不连续来预设裂纹时，并未弱化任何单元，或者说只将包含裂纹面的单元弱化了，因此对于该模型而言，采用应变 - 相场法来预设裂纹时，更多的单元被弱化，模型会更早开裂。进一步采用改进的应变 - 相场法来预设裂纹，其得到的裂纹扩展路径如图10所示。通过比照受拉方板的力 - 位移曲线(图11)，可以发现采用改进的方法会进一步提前裂纹的开裂，但该计算结果是合理的，与采用几何不连续来建立的模型结果差别较大，根本原因如前所述，数值模型有一些差别。

图9 受拉方板裂纹扩展路径(裂纹应变-相场建模)

图10 受拉方板裂纹扩展路径(改进的裂纹应变-相场建模)

图11 受拉方板的力-位移曲线

4.2 含单边裂纹的受剪切方板

采用与算例4.1同样的材料参数研究方板受剪切作用下的裂纹扩展过程，其计算模型如图12所示，所有的边界上竖向位移均约束为零。采用位移加载法，每一加载步位移增量为Δu=1×10-4mm，采用非均匀网格，在预计裂纹扩展的位置加密网格，最小网格大小在0.001 mm左右，网格模型如图13所示。

图12 受剪切方板/mm

图13 受剪切方板不均匀网格

首先采用几何不连续法来预设裂纹，其计算得到的裂纹扩展路径如图14所示，若采用应变-相场法来预设裂纹，计算得到的裂纹扩展路径如图15所示，可以发现与图14差别较大，材料很久才开裂，且开裂路径也不一致，主要原因是二者计算模型在该算例下差别较大，采用应变 - 相场法预设裂纹时，初始裂纹附近并未完全弱化，而采用几何不连续法预设裂纹时，裂纹面上下节点是分离的，可以有错位。