基于Newmark隐式时间积分方案的裂纹动态扩展的数值计算方法

郭德平，李 铮，彭森林，曾志凯，吴岱峰

(1. 叙镇铁路有限责任公司， 云南 昭通 657900； 2. 武汉大学 土木建筑工程学院， 武汉 430072； 3. 重庆市城市建设投资(集团)有限公司， 重庆 400023； 4. 重庆大学 土木工程学院， 重庆 400045)

随着21世纪我国基础设施大规模建设的进行，西部大开发战略的实施以及世界经济危机以来国家对基础设施建设的投资，我国的铁路、公路、大中型水电站建设以及南水北调、西气东输等工程将有大量的长大隧道需要建设.因此，隧道掘进机(TBM)在我国的应用前景非常广阔，我国对掘进机及其技术的需求猛增.

在工程建设和运营过程中，结构或者裂隙岩体会受到多种形式的动力作用(如爆破、地震等)，在动荷载的作用下更容产生失稳破坏，故针对裂纹动态扩展的研究引起了广泛的关注[1-3].Sun等[4]基于有限元方法比较了显式时间积分和隐式时间积分在分析线性和非线性动力问题中的区别.2012年，Nilsson等[5]在黏结单元中嵌入裂纹，分析了动荷载作用于弹塑性厚板的动力响应.

扩展有限单元法(XFEM)[6-9]是近十五年发展起来的一种在常规有限元框架内求解不连续问题，特别是裂纹扩展问题的数值方法.其原理是在裂纹影响区域的单元结点上用裂尖渐近位移场函数及跳跃函数加强传统有限元的基，以考虑由于裂纹存在引起的位移不连续性，继承了标准有限元的所有优点，但其所使用的网格与结构内部的几何和物理界面无关，从而避免了常规有限元方法中的网格重构，不需要裂纹面与有限元网格一致，不需要在裂缝周围布置高密度网格，大大简化了裂纹扩展的分析过程.能够很容易刻画出裂纹面上位移的强不连续性质和裂纹尖端的应力奇异性，并且无需重新划分网格.2001年，Stolarska等[10]将水平集函数引入XFEM来描述裂纹的位置和裂纹扩展之后的更新位置.裂纹的几何位置能够很容易用两个零水平集函数来表达，这两个零水平集函数在裂纹尖端处彼此正交.同时，随着裂纹的扩展，裂纹面和裂纹尖端处需要富集的节点能够实时更新.

Zhou等[11]在推导扩展有限元算法的基础上，分析了应力强度因子的J积分计算方法及积分区域的选取.基于扩展有限元法对I型裂纹进行了计算，有限元网格独立于裂纹面，无需在裂纹尖端加密网格.Chen等[12]引入裂纹交叉汇合加强函数以分析多裂纹交叉汇合过程，并且在裂纹附近区域使用广义形函数，并引入线增函数消除混合单元，可有效提高裂纹附近的精度.黄宏伟等[13]在这些研究的基础上，采用扩展有限元研究了衬砌在主要影响因素作用下的裂缝分布规律、裂缝扩展过程、裂缝外观表现形式及发生机制.阮滨等[14]基于扩展有限元法对均质土坝坝顶的初始裂纹扩展路径进行了模拟，研究结果表明扩展有限元法对网格划分的要求比较高，网格须均匀，网格的疏密程度对计算的结果影响不大.Menouillard等[15]提出利用无网格近似方法的XFEM富集方案，该方法采用显式时间积分方案来分析动力过程.Wen等[16]基于改进的扩展有限单元法研究了裂纹动态扩展问题的计算精度和算法稳定性问题.

但是标准有限元在处理时间积分时，在裂纹不断扩展的过程中整体刚度矩阵的自由度也会不断的增大，从而导致迭代计算无法进行.本文基于扩展有限单元法模拟动态裂纹扩展的方法，提出了新的时间积分方案.将所有节点都富集Heaviside函数和裂纹尖端的渐近位移场函数，即每个节点都有12个自由度，从而使得总体刚度矩阵式保持一致，避免迭代计算式无法进行的情况.提出了一种稀疏矩阵技术来解决矩阵所占内存大和计算时间长的问题.

1 动力扩展有限元的控制方程

1.1 控制方程的强形式

图1 在动荷载作用下含有裂纹的二维区域示意图Fig.1 Diagram of two-dimensional domain containing a crack subjected to dynamic loads

动力平衡方程的强形式为

(1)

该二维区域的应力边界条件和位移边界条件为

σ·nt=T, 在Γt上

(2)

σ·nc=0, 在Γc上

(3)

U·nu=0, 在Γu上

(4)

式中：nt为边界Γt的单位外法向量；nc为边界Γc的单位外法向量；nu为边界Γu的单位外法向量.

在小应变和小位移的条件下，几何方程可以表达为

(5)

根据广义胡克定律，σ与应变张量的之间的关系可以表达为

σ=C∶ε

(6)

式中：C为Hooke张量.

1.2 控制方程的等效弱形式

根据虚位移原理，引入虚位移δU，式(1)可以改写为

(7)

再根据分部积分，式(7)可以化简为

(8)

由变分原理将式(8)的第一项进行变换为

(9)

因此，该二维区域的能量系统的泛函数可以用下式表达

(10)

对式(10)进行变分，并结合式(8)可以得到：

(11)

因此有，

(12)

1.3 区域的空间离散

水平集法是裂纹和非连续面追踪强有力的工具.Zhou等[17]首次引入水平集法来研究材料内部界面的追踪定位和形状描述.如图2所示，一根裂纹可以用2个在裂纹尖端彼此相互正交的零水平集函数φ(x)和ψ(x)表示.

图2 用于追踪裂纹的两个零水平集函数原理图Fig.2 Schematic diagram of two zero-level set functions used to track cracks

Stolarska等[10]给出了水平集法的更新算法，水平集函数φ(x)的更新算法为

(13)

式中：φi+1(x)为水平集函数φ(x)的更新值；Δt为时间增量；ui和vi分别为裂纹尖端沿x和y方向扩展的速度.

对于水平集函数ψ(x)，其更新算法为

ψi+1(x)=

(14)

式中：F=[FxFy]为裂纹尖端在裂纹面外法线方向的速度矢量；(xi,yi)为裂纹尖端的坐标.

如图3所示，区域Ω被离散化成ne个单元，I表示该区域所有节点，J表示裂纹面上被Heaviside函数富集的节点集合，由蓝色正方形标识，K表示裂纹尖端上被渐近位移场函数富集的节点集合，由红色圆圈形标识.

图3 扩展有限单元法的节点富集方案Fig.3 Enrichment scheme of XFEM

根据文献[6，18]的研究，单元位移场可以表达为

(15)

(16)

(17)

c4x1c4y1c4x2c4y2c4x3c4y3c4x4c4y4]T

(18)

式中：aix和aiy分别为第i个节点在x和y方向的自由度位移值；bix和biy分别为第i个Heaviside函数富集节点在x和y方向的附加自由度位移值；cixj和ciyj分别为第i个节点的第j个附加自由度沿x和y方向的附加自由度位移值.

根据Wu等[19]的研究，线弹性材料中的I型、II 型及 III 型裂纹尖端渐近位移场均可由几个特定的基函数组成的函数形式来表达，为了能体现出裂纹尖端渐近位移场的奇异性，将该组基函数引入裂纹尖端的位移场计算，如下：

(19)

(20)

(21)

式中：(r,θ)是以裂纹尖端为原点建立的极坐标系的坐标值；(xtip,ytip)是绝对坐标系中裂纹尖端的坐标值，坐标系的建立如图4所示，x′O′y′为以裂纹尖端建立的局部坐标系，其中α为裂纹中心线与绝对坐标系x轴的夹角.

图4 裂纹面上的两套坐标系Fig.4 Two coordinate systems of crack surface

裂纹尖端富集后的形函数Nφ(x)可以表达为

(22)

φj为尖端位移场.将式(15)代入式(5)，可以得到单元的应变张量εe(x)如下：

(23)

再将式(15)、(23)代入式(10)中可以得到：

(24)

将式(24)代入式(11)并化简，可以得到：

(25)

式中：

考虑到3种类型的自由度位移值Ua、Ub及Uc式彼此相互独立的，结合变分原理，将式(25)代入式(12)，可以得到：

(26)

根据式(26)，扩展有限单元法的运动学方程为

(27)

1.4 时间积分方案

在扩展有限单元法的计算中，时间积分会遇到困难.时间积分基于迭代的算法，而在裂纹不断扩展的过程中整体刚度矩阵的自由度也会不断的增大，从而导致迭代计算无法进行.在每次迭代中，时间积分会涉及当前步的位移场Ut和下一步的位移场Ut+Δt，具体表达式为

(28)

(29)

本文提出的时间积分方案是将所有节点都富集Heaviside函数和裂纹尖端的渐近位移场函数，即每个节点都有12个自由度，从而使得总体刚度矩阵式保持一致，而避免迭代计算式无法进行.但是，将所有节点都富集Heaviside函数和裂纹尖端的渐近位移场函数的代价是刚度矩阵的阶数增多，从而使得参与运算的矩阵所占内存和计算时间急剧增加.为此，本文提出了一种稀疏矩阵技术来解决矩阵所占内存大和计算时间长的问题.

对于两个常规有限元的自由度，若该节点的某个方向被约束(约束状态)，那么所对应的自由度值被定义为0，若该节点的某个方向没有被约束(激活状态)，那么所对应的自由度值作为未知数参与计算.对于两个Heaviside函数富集的自由度，若该节点不属于集合J时(约束状态)，这两个自由度位移值被定义为0，若该节点属于集合J时(激活状态)，这两个自由度位移值作为未知数参与计算.对于另外8个裂纹尖端渐近位移场函数富集的自由度，若该节点不属于集合K时(约束状态)，这8个自由度位移值被定义为0，若该节点属于集合K时(激活状态)，这8个自由度位移值作为未知数参与计算.

根据上述定义，可以得到一个将每个节点有两个自由度的体系(简称二维体系)映射到每个节点有12个自由度的体系(简称十二维体系)中的转换矩阵，该转换矩阵的构造方法如下：

(30)

式中：变量n和m根据模型的节点个数来定.二维体系中的第i个自由度映射于12维体系中的第j个自由度，当该自由度是激活状态时，该矩阵中所有元素为1，当该自由度是约束状态时，该矩阵中所有元素为0.

(31)

(32)

(33)

(34)

那么，扩展有限单元法在12维体系中的运动方程为

(35)

对于时间积分，本文采用Newmark隐式时间积分方案如下：

(36)

(37)

式中：ξ和η分别为Newmark隐式时间积分方案中的参数.

在Newmark隐式时间积分方案中，t+Δt时刻的位移场为

(38)

联合式(36)～(38)，可以得到12维体系的位移场的求解方程式：

(39)

该方程的初始条件是：

(40)

(41)

(42)

在12维体系中，t+Δt时刻的加速度场和速度场为

(43)

(44)

根据Newmark时间积分方法，式(39)求解稳定的条件是ξ和η应满足[20]：

ξ≥0.5，η≥0.25(0.5+ξ)2

(45)

2 实例与分析

2.1 实例1

如图5所示，平面板的长L=10 m，高H=10 m.预制裂纹在板的左侧中部，其长度l0=5 m，距离上边界高度h=2 m.板的下边界竖向被约束，左下角被水平方向约束.板的材料属性为：弹性模量E=210 GPa，泊松比ν=0.3，密度ρ=8 000 kg/m3，阻尼系数μ=0.05.上边界受到动荷载σ(t)的作用，作用力的函数表达式为

图5 实例1的几何布置和荷载条件(m)Fig.5 Geometric layout and load conditions of example 1 (m)

σ(t)=σ0f0(t)

(46)

根据Wang等[22-23]的研究，裂纹尖端的应力强度因子KI的解析解为

(47)

(48)

(49)

采用本文提出的方法计算不同时刻裂纹尖端的动态应力强度因子，再将其归一化.采用了3种不同的时间步Δt=10 μs，Δt=20 μs及Δt=50 μs，所得数值结果如图6所示.裂纹尖端动态应力强度因子在时间t=0～τc内几乎为0，因为这段时间应力波还没有到达裂纹的尖端.当t>τc时，裂纹尖端的动态应力强度因子开始逐渐增大.从如图6可以看出，本方法计算的结果与解析解的结果吻合.本方法计算的的动态应力强度因子具有一定的震荡性，但是震荡性较小，如图7所示(KI,S为震荡值)，并且震荡性随着时间的增长而逐渐衰弱.本方法计算的x方向应力σx的分布如图8所示，在裂纹尖端的应力集中比较明显.

图6 KI的历时曲线图Fig.6 Time history of KI

图7 不同时间步KI的震荡历时曲线图Fig.7 Time history of oscillation of KI at different time steps

图8 σx分布云图(Δt=10 μs，t=3τc)Fig.8 Contour of σx (Δt=10 μs，t=3τc)

2.2 实例2

图9 实例2的几何布置和荷载条件(m)Fig.9 Geometric layout and load conditions of example 2 (m)

σ(t)=σ0f0(t)

(50)

在本实例中，研究了25×100，50×200和70×280共3种不同的网格密度(行数×列数).不同网格密度下裂纹动态应力强度因子的历时曲线如图10所示，网格密度为25×100和50×200的裂纹动态应力强度因子的相关系数为 0.999 3，网格密度为50×200和70×280的裂纹动态应力强度因子的相关系数为 0.999 4，网格密度为25×100和70×280的裂纹动态应力强度因子的相关系数为 0.998 5.当t<0.5 ms时，裂纹尖端的应力强度因子几乎为0；当t≥0.5 ms时，裂纹尖端的应力强度因子随着时间的增加而逐渐增大，并呈现很小的震荡特性.变形后的简支梁如图11所示，裂纹基本沿着直线向上扩展.该简支梁变形后的应力分布如图12所示，裂纹尖端应力集中突出，水平方向应力值达2.09×103MPa.

图10 不同网格密度下KI的历时曲线图Fig.10 Time history of KI at different mesh densities

图11 变形后的简支梁(网格密度为50×200)Fig.11 Deformed mesh for mesh density (at a mesh density of 50×200)

图12 t=2 ms时刻的σx云图(网格密度为50×200)Fig.12 Contour of σx at t=2 ms (at a mesh density of 50×200)

3 结语

标准有限元在处理时间积分时，在裂纹不断扩展的过程中整体刚度矩阵的自由度也会不断增大，从而导致迭代计算无法进行.本文提出的基于扩展有限单元法模拟动态裂纹扩展的方法提出了新的时间积分方案，将所有节点都富集Heaviside函数和裂纹尖端的渐近位移场函数，即每个节点都有12个自由度，从而使得总体刚度矩阵式保持一致，避免迭代计算式无法进行.并且提出了一种稀疏矩阵技术来解决矩阵所占内存大和计算时间长的问题.数值计算的结果表明，利用空间变换理论计算动力问题时，施加的荷载以膨胀波的形式传递，裂纹尖端的应力强度因子比荷载施加的时间稍有延迟，数值解的结果能够与解析的结果很好地吻合.线性荷载的数值结果比瞬时脉冲荷载的数值结果的震荡性更小，应力分布更加平稳光滑.不同网格密度条件下的数值计算结果相差不大，说明该方法计算裂纹扩展对网格的依赖性较小.