本帖最后由 shifang157300 于 2018-6-15 09:09 编辑
扩展有限元法(XFEM)漫谈 发布日期：2016-08-21（2017-08-15更新）
本贴第3部分更新了我编写的扩展有限元程序PhiPsi程序（phipsi.top）相关介绍和应用实例。更新日期：2016-10-22
2016-10-27：新增三维裂缝XFEM实例及输入文件
2017-08-15：本贴第4部分新增PhiPsi GUI程序的介绍及下载链接
0. 引言 2012年笔者在图书馆偶然看到清华大学庄茁教授编写的《扩展有限单元法》一书[1]，意识到扩展有限元法（XFEM）可能是一种很有应用前景的数值方法，遂跟导师商量更换了研究方向，确定了“用扩展有限元方法模拟岩石破裂过程”这一方向，详见 [2]。而在此之前，笔者的主要研究内容是“基于ANSYS的岩石材料本构开发”，见[3]。在博士研究课题进行过程中，笔者基于Matlab编写了扩展有限元程序并用于计算博士论文中的全部算例[2]。博士研究课题完成后，笔者花了近半年时间将原来基于Matlab编写的程序用Fortran进行了重写和优化，功能、计算速度和稳定性都得到了很大的提升。笔者目前的研究方向是基于扩展有限元法的水力压裂过程模拟。四年多的XFEM学习和科研过程中，不断的遇到各种问题然后耐心的解决，尤其是近两万行Fortran程序的编写过程中，笔者积累了不少关于XFEM的经验，而这些经验对于初学者而言可能是十分有用的，看过这篇帖子，你或许可以少走一些我走过的弯路。把经验总结出来呈献给大家，正是本贴的初衷。由于帖子中有些公式，而在帖子编辑的时候没办法输入，只能截图，所以帖子看起来不是特别整洁。 欢迎大家通过QQ与我交流：1549221758。
1. XFEM简介
扩展有限元法由Belytschko教授提出于1999年[4]。有限元法在模拟断裂问题时的缺点：（1）裂缝每扩展一步，都要进行网格划分；（2）裂尖网格必须划分的很密，以便较为精确的模拟裂尖应力场以及计算裂尖应力强度因子；（3）某些情况下，新旧网格之间需要进行数据的转换，这会进一步增加模拟过程复杂度。而XFEM恰好能克服有限元法的上述缺点，对于XFEM：（1）裂纹扩展过程中无需重新划分网格；（2）裂尖不需要加密；（3）因为使用同样的网格，所以不存在数据传递的问题。对于扩展有限元法，单元内任意高斯点 x 的位移可表示为下式（该式是XFEM的核心）: (1-1)
其中， n为单元常规有限元节点数， Nj为形函数， uj 常规有限元节点自由度向量， mh为裂纹面两边增强节点数, H( x)是高斯点 x处的Heaviside函数值， H( xh)是增强节点 h处的Heaviside函数值, ah为裂纹面两边增强节点自由度向量； mt为裂尖增强节点数, Fl( x)是裂尖增强函数在高斯点 x处的值， Fl( xk)是裂尖增强函数在增强节点 k处的值， bk为裂纹尖端增强节点自由度向量。Heaviside增强函数在裂纹面的一侧取值为1，另一侧取值为-1 。 Fl 是裂尖增强函数，是从线弹性断裂力学裂尖位移场的解析解表达式中提取出来的线性无关的一组基，其定义在裂尖极坐标系下，对于各向同性材料其表达式为： (1-2)
其中， 是高斯点x在裂尖极坐标系下的坐标。裂尖增强函数F1至F4的图像如图1所示，从其中可以看出，只有第一项，即F1与裂缝面两侧的位移跳跃有关（图1(a)是撕开的，其余三个都是连续的），因此，某些情况下为了减少计算量可以只取第一项[5, 6]。
图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。
图2.裂缝增强示意图
如何确定增强节点是，是编写扩展有限元程序的关键之一。裂尖增强节点确定的原则是：支撑域（如图3所示）被裂缝完全切割成两部分的节点需要进行Heaviside增强，支撑域包含裂尖的节点则需要进行裂尖增强。
图3.某节点k的支撑域
2. 关于XFEM的注意点
2.1 Heaviside增强节点的选择:裂缝两边的支撑域都应该含有Gauss点
在选择Heaviside增强节点时，必须保证裂缝两侧的支撑域都含有Gauss点，否则不予增强。比如图4(b)中的节点 i，该节点的裂纹上方的支撑域不包含任何高斯积分点，所以节点 i不进行Heaviside增强。这一点尤其需要注意，否则编写的程序计算出来的裂缝面附近的某些节点的位移以及应力分布会出现莫名其妙的错误。详见文献[2]第20页。
(a)增强节点i (b)不增强节点i
图4.Heaviside增强节点的选择
2.2 含折点的某些情况下必须要特殊处理（1-2）式中的裂尖增强函数仅适用于直线裂纹，对于带折点的裂纹，需要做特殊处理，否则将会因为裂纹面两侧Heaviside增强和裂尖增强的相互矛盾而导致计算出错。例如，图2-4为一折线裂纹，其中用“○”标记的是Heaviside增强节点，用“□”标记的是裂尖增强节点，对于阴影区域内的高斯积分点而言，位于裂纹片段BC的上侧，Heaviside增强函数值为1，但这些高斯点却位于裂纹片段AB延长线的下侧，二者是矛盾的。Belytschko[4]采用映射法，给出了这一问题的解决办法。
图5.折线裂纹示意图 本注意点非常重要，在编写程序时必须要进行相应的处理才行，否则会经常遇到裂尖附近应力场计算错误的问题。详见我的博士论文[2]第19页。
2.3 裂缝开度的计算只跟增强自由度的值有关，而与FEM自由度无关在扩展有限元程序中，经常需要计算并保存裂缝的开度。如图7所示，裂纹上 x点对应的裂纹面两侧的两个点分别记为 和 ，则由（1-1）式可得： (2-3)
(2-4)
则 x点的裂隙开度 w为：
(2-5)
由 N( x+)= N( x-)= N( x) ， x+和 x- 位移的区别仅仅是Heaviside增强函数值和裂尖增强函数的第一项的值不同（见图1），于是 x点的裂隙开度 w可进一步写成：
(2-6)
图7.裂隙开度计算示意图
本小节内容详见文献[2]第61页。
2.4 对于裂缝面承受载荷的裂缝，这些载荷是作用于增强自由度上的
对于水力压裂分析（裂缝面承受水压）、爆破分析（裂缝面承受高压气体作用）、裂缝面接触滑动（法向和切向接触力作用于裂缝面）、粘聚裂纹模型（粘聚内力作用于裂尖附近的裂缝面上）等裂缝表面承受载荷的分析，注意：这些载荷是作用于增强自由度上的，而不是FEM自由度。这一点对于初学者来说或许是个易犯的错误。
2.5 关于最大周向拉应力准则裂缝扩展路径的锯齿状震荡
最大周向拉应力准则[7]是XFEM文献中使用最多的准则。但是利用该准则计算出来的裂缝扩展路径往往呈现出锯齿状（如图8(a)所示），与试验观察到的裂缝路径往往不符。出现锯齿状路径的原因是，最大周向拉应力准则预测的裂缝扩展方向与 的绝对值直接相关，某些情况下 的值较大，且 的符号交替变化，造成锯齿状的裂缝扩展路径。虽然加密网格可以是锯齿看起来不那么明显，但是计算量会增大，得不偿失。所以，某些情况下，可以采用基于应力场的裂缝扩展准则，比如加权平均最大主拉应力准则来预测裂缝扩展方向。这类准则所得扩展路径光滑无锯齿，如图8(b)所示。
(a)最大周向拉应力准则
(b) 加权平均最大主拉应力准则
图8.Mises应力云图及裂缝扩展路径
比如，在商业软件Abaqus中使用XFEM进行裂纹扩展模拟时，由于Abaqus不对裂尖进行增强，而只有Heaviside增强，因此裂缝每步扩展都使得裂尖位于单元边界上。在自己编写扩展有限元程序时，也应按照这一规则。
2.7 关于Gauss积分点的确定，有两种方案
第一种方案，依据裂缝位置将增强单元划分成一系列的子三角形，然后依次对每个三角形进行积分。第二种方将单元剖分成若干个小的四边形，每个四边形取4个积分点。第一种方案精度比后者高，但由于简单便捷，所以编程时建议选择后一方案，一般来说，积分点越多越精。
图9.两种积分方案示意图
PhiPsi是我在博士和博士后期间基于Fortran语言编写的通用扩展有限元程序，独立代码行数约三万行，可实现水力压裂、地震作用、爆破冲击、复合材料损伤等过程中裂缝扩展过程的模拟，此外还有传统有限元计算功能、温度场分析、热应力计算、复合材料物性参数预测等功能。该软件我至今仍在持续编写完善，最新版我会定期上传到对应网站phipsi.top上，该网站上还有早期版本的PhiPsi Fortran源代码供下载（需要最新完整功能源代码的可以联系我QQ:1549221758）。之所以取名PhiPsi是因为两方面原因：（1）XFEM框架下水平集函数用Φ和Ψ描述；（2）塑性力学中Φ和Ψ用来描述内摩擦角和剪胀角。
以下例子网站上均有相关介绍。
3.1 程序下载
PhiPsi的运行需要关键字文件（PhiPsi关键字手册：）的支持，类似于ANSYS的命令流，在其中定义好输入文件的路径，输入文件的文件名，分析类型，初始裂缝坐标等。对应2D问题，PhiPsi的输入文件包括6个，分别定义节点坐标、单元节点编号、x方向载荷、y方向载荷、x方向约束、y方向约束。笔者编写了一个ANSYS宏文件，可以在ANSYS中建模后快速导出PhiPsi所需的输入文件。笔者还基于Matlab编写了功能完善的后处理程序，可以绘制变形图、应力云图、绘制各种曲线（如应力强度因子）、生成各种动画，3.2小节中的应用实例均由该程序生成。该Matlab后处理程序源代码完全公开： 程序下载，适用于Windows64位系统：
一个演示算例（Toturial）及对应文件。
3.2 应用实例
（1）“T”型交叉裂缝的形成过程模拟
相关文件下载：
裂缝扩展过程
（2）含夹杂复合材料的拉伸模拟
相关文件下载：
模型示意图
网格变形图
（3）水力压裂模拟
相关文件下载：
注水点水压变化曲线
水力裂缝扩展过程
（4）地震作用下大坝裂缝的扩展过程模拟
相关文件下载：
裂缝尖端应力强度因子的变化情况（未考虑接触）
裂缝扩展路径
（5）3维裂缝的扩展有限元模拟
相关文件下载：
3D网格变形、增强节点（蓝色圆点）及裂缝面（红色）
(6)高内压作用下空心圆筒筒壁三维裂缝相关文件下载：
网格图 及变形图 (7)多裂缝扩展（9个初始裂缝），详见PhiPsi.top网站
(8)粘聚裂缝三点弯曲梁，详见PhiPsi.top网站 4 关于PhiPsi图形界面程序PhiPsi (更新时间2017-08-15) 最近作者采用vb.net语言开发了PhiPsi GUI (图形用户界面)，可直接在PhiPsi GUI中建模、计算并进行后处理。下载链接：http://phipsi.top/downloads_cn.html
5 一些关于扩展有限元的资料 (1)我最早看的一篇硕士论文（英文）：
(2)我的博士论文： (3)我的关于扩展有限元水力压裂的论文：
参考文献
[1] 庄茁，柳占立，成斌斌，等.扩展有限单元法[M]. 北京： 清华大学出版社,2012.
[2] 师访.岩石破裂过程的扩展有限元法研究[D]. 徐州：中国矿业大学,2015.
[3] 师访.ANSYS二次开发及应用实例详解[M]. 北京： 中国水利水电出版社,2012.
[4] Belytschko T, Black T. Elastic crack growth infinite elements with minimal remeshing[J]. International Journal for NumericalMethod in Engineering, 1999, 45(5)：601–620.
[5] Elguedj T, Gravouil A, Maigre H. An explicitdynamics extended finite element method. Part 1：Mass lumping for arbitrary enrichment functions[J].Computer Methods in Applied Mechanics and Engineering, 2009, 198(30-32)：2297-2317.
[6] Menouillard T, Song J H, Duan Q, et al. Timedependent crack tip enrichment for dynamic crack propagation[J]. InternationalJournal of Fracture, 2010, 162(1-2)：33-49. [7] Sih G C. Strain Energy density factor appliedto mixed mode crack problems[J]. International Journal of Fracture Mechanics,1974, 10(3)：305–321.
| 
IP卡
狗仔卡
发表于 2016-8-21 20:00:45
提升卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
发表于 2016-8-22 21:50:20
