扩展有限元法(XFEM)漫谈-我的XFEM学习经验分享

1/14扩展有限元法(XFEM)漫谈0.引言2012年笔者在图书馆偶然看到清华大学庄茁教授编写的《扩展有限单元法》一书[1]，意识到扩展有限元法（XFEM）可能是一种很有应用前景的数值方法，遂跟导师商量更换了研究方向，确定了“用扩展有限元方法模拟岩石破裂过程”这一方向，详见[2]。而在此之前，笔者的主要研究内容是“基于ANSYS的岩石材料本构开发”，见[3]。在博士研究课题进行过程中，笔者基于Matlab编写了扩展有限元程序并用于计算博士论文中的全部算例[2]。博士研究课题完成后，笔者花了近半年时间将原来基于Matlab编写的程序用Fortran进行了重写和优化，功能、计算速度和稳定性都得到了很大的提升。笔者目前的研究方向是基于扩展有限元法的水力压裂过程模拟。四年多的XFEM学习和科研过程中，不断的遇到并解决各种问题。尤其是在3万多行Fortran程序的编写过程中，笔者积累了不少关于XFEM的经验，而这些经验对于初学者而言或许有些益处。1.XFEM简介扩展有限元法由Belytschko教授提出于1999年[4]。有限元法在模拟断裂问题时的缺点：（1）裂缝每扩展一步，都要进行网格划分；（2）裂尖网格必须划分的很密，以便较为精确的模拟裂尖应力场以及计算裂尖应力强度因子；（3）某些情况下，新旧网格之间需要进行数据的转换，这会进一步增加模拟过程复杂度。而XFEM恰好能克服有限元法的上述缺点，对于XFEM：（1）裂纹扩展过程中无需重新划分网格；（2）裂尖不需要加密；（3）因为使用同样的网格，所以不存在数据传递的问题。对于扩展有限元法，单元内任意高斯点x的位移可表示为下式（该式是XFEM的核心）：xfem11411()()()()()()()()nmhjjhhhjhmtlkllkkklNNHHNFFuxxuxxxaxxxb(1-1)其中，n为单元常规有限元节点数，Nj为形函数，ju为常规有限元节点自由度向量，mh为裂纹面两边增强节点数，()Hx是高斯点x处的Heaviside函数值，()hHx是增强节点h处的Heaviside函数值，ha为裂纹面两边增强节点自由度向量；mt为裂尖增强节点数，()lFx是裂尖增强函数在高斯点x处的值，()lkFx是裂尖增强函数在增强节点k处的值，lkb为裂纹尖端增强节点自由度向量。Heaviside增强函数在裂纹面的一侧取值为1，另一侧取值为1。Fl是裂尖增强函数，是从线弹性断裂力学裂尖位移场的解析解表达式中提取出来的线性无关的一组基，其定义在裂尖极坐标系下，对于各向同性材料其表达式为：2/1441,sin,cos,sinsin,sincos2222llFrrrrr(1-2)其中，,r是高斯点x在裂尖极坐标系下的坐标。裂尖增强函数F1至F4的图像如图1所示，从其中可以看出，只有第一项，即F1与裂缝面两侧的位移跳跃有关（图1(a)是撕开的，其余三个都是连续的），因此，某些情况下为了减少计算量可以只取第一项[5,6]。(a)sin2r(b)2cosr(c)sinsin2r(d)sincos2r图1.裂尖增强函数F1至F4图像一个典型的裂缝增强例子见图2。其中包含6*9=54个FEM节点，对应54*2=108个自由度；22个节点被增强，对应92个自由度{（8个节点被裂尖增强*4+14个节点被Heaviside增强=46）*2=92}。所以，图2所示模型的总自由度数目是108+92=200。那么集成的系统整体刚度矩阵K的维度就是200*200。3/14图2.裂缝增强示意图如何确定增强节点是，是编写扩展有限元程序的关键之一。裂尖增强节点确定的原则是：支撑域（如图3所示）被裂缝完全切割成两部分的节点需要进行Heaviside增强，支撑域包含裂尖的节点则需要进行裂尖增强。图3.某节点k的支撑域2.关于XFEM的注意点2.1Heaviside增强节点的选择—裂缝两边的支撑域都应该含有Gauss点在选择Heaviside增强节点时，必须保证裂缝两侧的支撑域都含有Gauss点，否则不予增强。比如图4(b)中的节点i，该节点的裂纹上方的支撑域不包含任何高斯积分点，所以节点i不进行Heaviside增强。这一点尤其需要注意，否则编写的程序计算出来的裂缝面附近的某些节点的位移以及应力分布会出现莫名其妙的错误。详见文献[2]第20页。4/14(a)增强节点i(b)不增强节点i图4.Heaviside增强节点的选择2.2含折点的某些情况下必须要特殊处理（1-2）式中的裂尖增强函数仅适用于直线裂纹，对于带折点的裂纹，需要做特殊处理，否则将会因为裂纹面两侧Heaviside增强和裂尖增强的相互矛盾而导致计算出错。例如，图2-4为一折线裂纹，其中用“○”标记的是Heaviside增强节点，用“□”标记的是裂尖增强节点，对于阴影区域内的高斯积分点而言，位于裂纹片段BC的上侧，Heaviside增强函数值为1，但这些高斯点却位于裂纹片段AB延长线的下侧，二者是矛盾的。Belytschko[4]采用映射法，给出了这一问题的解决办法。图5.折线裂纹示意图以图5中折线裂纹AB和BC为例说明通过映射法修改裂尖增强函数的过程。（1-2）式裂尖增强函数仅适用于图5虚线（线段AB在B点的垂线）右侧的高斯积分点，虚线左侧的增强函数需要进行修正。如图6(a)所示，高斯点的坐标为(x,y)，α为BA和B之间的夹角，θR为BA和BC之间的夹角。定义夹角：/2,3/2/2,/2RRRRRR(2-1)A－A+i裂纹A－A+i裂纹iABC5/14然后通过下式将图6(a)中的高斯点(x,y)映射到图6(b)中的(,)xy：BAB*B*(,)(cos,sin)xylll(2-2)上式中，BAl是线段BA的长度，22BA2tip2tiplxxyy，Bl是线段B的长度22B22lxxyy。于是，22rxy，arctanyx，将r和代入（1-2）式即为折线裂纹的裂尖增强函数。(a)初始状态(b)映射点的确定图6折线裂纹的映射过程本注意点非常重要，在编写程序时必须要进行相应的处理才行，否则会经常遇到裂尖附近应力场计算错误的问题。详见文献[2]第19页。2.3裂缝开度的计算只跟增强自由度的值有关，而与FEM自由度无关在扩展有限元程序中，经常需要计算并保存裂缝的开度。如图7所示，裂纹上x点对应的裂纹面两侧的两个点分别记为x和x，则由（1-1）式可得：xfem11411()()()()()()()()nmhjjhhhjhmtlkllkkklNNHHNFFuxxuxxxaxxxb(2-3)xfem11411()()()()()()()()nmhjjhhhjhmtlkllkkklNNHHNFFuxxuxxxaxxxb(2-4)则x点的裂隙开度w为：xfemxfem(()())wnuxux(2-5)x1y1(x2,y2)(x1,y1)(xtip,ytip)(x,y)*αθRABCx1y1(x2,y2)(xtip,ytip)*AB(,)xyC6/14由于()jNx=()jNx=()jNx，x和x位移的区别仅仅是Heaviside增强函数值和裂尖增强函数的第一项（sin/2r）的值不同而已（见图1），于是x点的裂隙开度w可进一步写成：1112()2()mhmthhkkhkwNrNnxanxb(2-6)图7.裂隙开度计算示意图本小节内容详见文献[2]第61页。2.4对于裂缝面承受载荷的裂缝，这些载荷是作用于增强自由度上的对于水力压裂分析（裂缝面承受水压）、爆破分析（裂缝面承受高压气体作用）、裂缝面接触滑动（法向和切向接触力作用于裂缝面）、粘聚裂纹模型（粘聚内力作用于裂尖附近的裂缝面上）等裂缝表面承受载荷的分析，注意：这些载荷是作用于增强自由度上的，而不是FEM自由度。这一点对于初学者来说或许是个易犯的错误。2.5关于最大周向拉应力准则裂缝扩展路径的锯齿状震荡最大周向拉应力准则[7]是XFEM文献中使用最多的准则。但是利用该准则计算出来的裂缝扩展路径往往呈现出锯齿状（如图8(a)所示），与试验观察到的裂缝路径往往不符。出现锯齿状路径的原因是，最大周向拉应力准则预测的裂缝扩展方向与III/KK的绝对值直接相关，某些情况下III/KK的值较大，且III/KK的符号交替变化，造成锯齿状的裂缝扩展路径。虽然加密网格可以是锯齿看起来不那么明显，但是计算量会增大，得不偿失。所以，某些情况下，可以采用基于应力场的裂缝扩展准则，比如加权平均最大主拉应力准则来预测裂缝扩展方向。这类准则所得扩展路径光滑无锯齿，如图8(b)所示。(a)最大周向拉应力准则xxxw7/14(b)加权平均最大主拉应力准则图8.Mises应力云图及裂缝扩展路径2.6假如不进行裂尖增强，则裂尖必须位移单元边界上比如，在商业软件Abaqus中使用XFEM进行裂纹扩展模拟时，由于Abaqus不对裂尖进行增强，而只有Heaviside增强，因此裂缝每步扩展都使得裂尖位于单元边界上。在自己编写扩展有限元程序时，也应按照这一规则。2.7关于Gauss积分点的确定，有两种方案第一种方案，依据裂缝位置将增强单元划分成一系列的子三角形，然后依次对每个三角形进行积分。第二种方案将Gauss积分点均匀的分布于四边形单元中。第一种方案精度比后者高，但由于简单便捷，所以编程时建议选择后一方案，一般来说，积分点越多越精确，对于每个包含增强节点的单元选择8*8，即64个Gauss积分点即可。图9.两种积分方案示意图3.一个简单的例子本部分通过一个简单的二维平面应变的例子来说明扩展有限元的求解过程。3.1问题描述如图10所示，平面应变模型尺寸3m*3m，裂缝AB的长度为1.5m。共划分9个单元，每个单元的尺寸均为1m*1m，模型共包含16个FEM节点。单元编号和节点编号已在图中标出（节点编号和单元编号均由ANSYS得到）。增强单元的积分方案为8*8积分（见2.78/14节）。材料弹性模量为20GPa，泊松比为0.2。本例子的目标是计算A处的裂缝开度。图10.模型示意图及增强节点（红色圆点表示裂尖增强、黄色原点表示Heaviside增强）3.2组集刚度矩阵K由于模型包含16个节点，因此有16*2=32个FEM自由度。模型有四个节点被裂尖增强（每个节点对应4个额外的裂尖增强节点），对应4*4*2=32个增强自由度。两个Heaviside增强节点，对应2*2=4个增强自由度。因此模型的总自由度数目为32+32+4=68。由于模型底部约束了5个自由度（1个x方向，4个y方向），因此参与计算的系统的活动自由度数目是63个，也就是组集的总刚矩阵K的维数为63*63。按照下式逐个计算单元刚度矩阵：ijijijeijijijijijijijuuuaubauaaabbubabbKKKKKKKKKK(3-1)其中TΩdΩ,,,ersrsijijrsKBBuabD(3-2)上式中，D为应力-应变关系矩阵，B为形函数的导数矩阵：9/14,,,,,,,,,,,,00000014ixiiyiyixiixiiiyiiiiyxiixiiiyiiiiyxNNNNNHHNHHNHHNHHNFFNFFNFFNFFuabBBB(3-3)以4号单元为例进行说明。该单元自由度数目为（4+2+2*4）*2=28。