DEIP：不连续单元插入程序-界面有限元建模的网格

SoftwareX 7（2018）162原始软件出版物DEIP，不连续单元插入程序-界面有限元建模的网格Timothy J. 委托人土木与环境工程系，田纳西大学，诺克斯维尔，318约翰D。Tickle Engineering Building，诺克斯维尔，TN 37996，美国ar t i cl e i nf o文章历史记录：2017年11月3日收到收到修订版，2018年5月10日接受，2018年保留字：界面要素间断Galerkin方法网格拓扑代码元数据a b st ra ct不连续单元插入程序是一个MATLAB/Octave工具箱，用于将零厚度界面单元插入到二维和三维有限元网格中。这些界面单元，在本文中称为底层算法是基于拓扑的，适用于复杂的，非结构化的混合型线性和二次元素的网格。插入是根据整个分析域中的区域或子域指定的，这是一种指定耦合器位置的几何直观方法©2018作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。当前代码版本DEIP v1.3此代码版本使用的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-17-00085伊利诺伊大学/NCSA开源许可证（NCSA）使用Mercurial的代码版本控制系统使用MATLAB的软件代码语言、工具和服务编译要求，操作环境依赖性MATLAB 2014 a如果可用，链接到开发人员文档/手册http://clmi.utk.edu/software/问题支持电子邮件ttruster@utk.eduCode Ocean compute capsulehttps://codeocean.com/2018/06/18/deip-discontinuous-element-insertion-program/code1. 动机和意义近年来，计算固体力学中的不连续公式越来越受欢迎.与断裂力学相关的物理应用包括脆性材料的碎裂[1，2]、复合材料中的分层[3类似地，不连续Galerkin（DG）方法[9数值实现这些方法的常用方法是使用所谓的零厚度界面有限元[13]，其中通过节点复制创建单元，使得界面的两侧最初重合，但随后可能分离。不幸的是，目前，标准的商业有限元代码电子邮件地址：ttruster@utk.edu。https://doi.org/10.1016/j.softx.2018.05.002不包含零厚度元素的网格生成功能因此，负担被放置在用户上以创建修改的网格连接性，这对于三维中的复杂网格而言可能是非常重要的各种研究人员已经提出了算法，用于生成零厚度元素的专门应用。一些最早的例子是为了在碎裂模拟过程中自适应地插入元素而开发的[1，14]。这些算法支持基于单元局部断裂准则的非本征内聚区（CZ）建模方法[15-这些算法还没有在开源软件中实现其他方法，如内聚区模型和不连续伽辽金方法，则要求所有界面元素从模拟开始时就存在于预定位置。这些方法可以受益于并发插入算法，其中基于直观的几何特征来确定界面小平面的列表物理领域。最近，已经开发了用于插入零厚度的算法和开放源代码[182352-7110/©2018作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxT.J. Truster / SoftwareX 7（2018）162163（）′结构域和HΩ是标准的希尔伯特空间，内聚区模型的一个重要条件是，Ω′使得对于所有的w∈H1（H′）：作为预处理步骤，该步骤已应用于使用混合DG-CZ公式[2]对晶界开裂、复合材料分层和基体开裂[19]进行但是，只有当所有区域间边界都需要界面元素而不是选择界面元素时，才能应用这种方法。因此，需要一种为用户提供选择界面元素的位置和类型的更大灵活性的方法在目前的工作中，从通用算法[20]开发了一个开源软件包，用于将零厚度界面元素（本文称为“耦合器”）插入二维和三维一致网格的指定区域。术语“耦合器”的引入是因为算法适用于具有不同连接模板的两种常见类型的界面元素：内聚区元素和不连续Galerkin元素。区别特征在于，用户通过（1）提供整个分析域的区域内的元素列表和（2）命名耦合器应插入的区域对来隐式地指定耦合器插入的位置。通过这种几何直观的手段，耦合器可以选择性地插入特定区域内或沿着特定界面插入。此外，不同类型的耦合器以及不同的材料特性可以直接分配给这些特定的集合，从而提供直观的手段来完成对界面修改的网格的描述以用于分析目的。这些界面网格的生成是纤维复合材料[5，21]中脱粘以及多晶金属[22，23]中晶界滑动和空穴的界面机械性能建模的关键因素2. 软件描述DEIP软件包包含一个包含耦合器插入和有限元分析脚本的源目录，以及一个演示模块的示例这些文件与MATLAB和Octave都兼容安装脚本将源目录加载到用户的搜索路径中，用户手册记录了软件包的用法。以下部分描述了用于指定耦合器位置的基本算法和区域概念，并列出了使用该包的模块和2.1. 拓扑定义驱动耦合器插入的区域概念DEIP是从由二维（2D）或三维（3D）空间中的有限元一致网格组成的域的拓扑结构导出的[24]。在下面的讨论中，拓扑实体用斜体字标识。网格中的每个元素由一组节点定义，这些节点是与其特定坐标相关联的域内的点;参见图中的2D示例。1.一、在2D中，包含三角形混合的网格包含四面体、棱锥体、楔形体或六面体元素的是允许的。 在图1中，节点由数字表示，在复合材料中，每个区域被认为具有不同的材料特性。然而，在本文中，区域是纯粹的几何构造，以实现完全通用的应用。域中所有元素的刻面可以分成三个不相交的集合。第一组是位于域边界上的那些小平面，它们恰好与一个元素相邻。第二组是位于两个不同区域的元素之间的那些元素，它们被称为属于界面。该集合根据区域对进一步划分，使得界面（A，B）是区域A和区域B之间的所有刻面的集合。第三组是位于同一区域的元素之间的那些。属于区域C的这种刻面将被表示为面内（C），等等。界面在图中以粗线段表示。 1，而面内显示为细线段。本文定义的最后一个相关拓扑概念是截面。如果两个元件之间的共享刻面不是指定由耦合器切割的刻面（界面或面内），则两个元件被定义为属于扇区扇区是满足这个定义的最大元素集合。例如，元素b、c和d在图5中形成围绕节点5的扇区。1.一、2.2. 耦合器插入算法本文提出了一种基于拓扑的耦合器插入算法。在有限元文献中，这样的计算实体通常被称为在这里，我们使用耦合器这个术语来区别于2.1节中的其他拓扑定义，并提供更广泛类型的计算实体。因此，耦合器被定义为拓扑单元，其由来自恰好两个元件的节点组成，所述两个元件在界面或面内相邻。耦合器通过复制位于由两个元件共享的小平面上的节点来生成，以沿着该小平面有效地分割网格。这些耦合器通常出现作为建模偏微分方程的不连续公式的数值实现。例如，为了模拟复合材料中的渐进脱粘，可以引入本征内聚区模型作为纤维元件和基体成分之间的界面耦合器耦合器从初始阶段就存在于分析中，以便通过牵引-分离关系捕获界面处断裂的开始和进展关于内在内聚区模型，读者可以参考[1，2，5，25]的数学方面和[18，26]的实施说明。简而言之，令f表示具有Lipschitz边界并且包含（可能多个和/或不相交的）界面ΓI的开有界域，界面Γ I将f划分为区域f’的集合，使得′的闭包等于此外，我们表示位移场为u，由位移导出的线弹性应力为：σ=C：ε（u），材料模量为C，对称梯度算子为：ε。弱f（或m）00元素由小写字母表示。例如，元素a由三个节点1、4和5组成。术语小平面是指元素边界上的2D边或3D面。除了上述与有限∫Ω′ε（w）：C：ε（u）d+[w]·t（[u]）dΓ=πw·bd（一）在元素离散化中，我们将区域定义为域内元素的连续集合，其通常可以形成非凸子域并且由空间上不相交的元素集合组成。域中的每个元素都是一个区域的成员。图在图1（a）中，区域用大写字母表示，属于每个区域的元素共享相同的颜色。在有限元建模的上下文中的区域的示例包括纤维和微纤维。其中b是体积力，t（[u]）是由穿过界面ΓI的位移跃变[u]定义的牵引场，10齐次Dirichlet边界条件在边界上，类似地，不连续Galerkin方法的公式[9，11]和实现[10，12]可以在所示考虑较大和四边形单元，在3D中，ΓIΩ′164T.J. Truster / SoftwareX 7（2018）162（）（）⊆0联系我们Fig. 1. 示例2D线性三角形有限元网格：（a）包含说明第2.1节中定义的界面和内面的区域的一致性网格;（b）根据第2.2节中的算法，沿切割界面和内面进行节点复制后得到的网格。参考文献，其弱公式可以表示为：u∈H1 使得对所有w∈H1 “你”：单元网格生成软件包，这意味着用户需要最少的数据准备接下来，数据（4）指定0∫′ε（w）：C：ε（u）d−Ω第一章0{C：ε（w）·n}在哪些拓扑位置插入耦合器。例如，纤维区域A和基质区域B之间的界面将通过标记界面（A，B）来指示;在图1中， 1（a），这表示为·（[u]−n）dr−n第一章ΓI[w]·{C：ε（u）·n}dΓ∫Ω（二）这些区域之间的虚线。在此，耦合器将被插入的这些刻面类似地，区域A内的面内（A）刻面在图1中以虚线表示。 1（a）说明[w]·τs·（[u]−ε）dΓ=′w·bdε其中Cε（u）n是界面牵引力的某个加权平均值，τs是稳定化参数，τs是通过使用牵引场和位移跃变的本构表达式导出的ΓI耦合器的作用是提供额外的拓扑数据，以评估在（1）和（2）中发现的ΓI上的界面积分以及其他不连续插值方法。注意，当在某些区域ω n′内需要不连续插值时，即单元内面，则这些区域中的函数空间H1（ω）需要由跨越有限单元的分段多项式的适当空间代替，并且将额外的面内积分添加到（2）;详见[9]。此外，一类广泛的不连续Galerkin方法通过允许边界值问题的连续离散化（例如具有悬挂节点的网格）提供更大的建模灵活性[11，27]。虽然本文所述的软件确实处理跨单元内面的不连续插值，但是该软件不适用于生成非协调耦合器;非协调离散化在专门的自适应网格程序中处理，如deal.II[28]。从最初一致的有限元网格开始，耦合器的插入通过参考各个区域之间的界面和内面来完成。该区域是用于分配耦合器的期望位置的自然几何实体。例如，当模拟纤维复合材料中的脱粘或沿多晶体中的晶界的空化时，不同材料区域之间的界面是期望的位置。类似地，当使用不连续Galerkin数值方法或当用内聚区模型模拟一般裂纹扩展时，面内是耦合器的自然位置该算法的输入数据是（1）节点的空间坐标，（2）连接到每个元件的节点的列表，（3）属于每个区域的元件的列表，以及（4）需要耦合器的接口和内表面的数据（1）(3) 是由几乎所有有限的他们也将被削减。随后，算法确定要插入的耦合器组和要复制的节点组基于网格拓扑的自动化过程耦合器插入算法由六个阶段组成，总结如下：1. 构造附加到每个节点的元素2. 将网格中的所有面分类为边界、界面或面内。3. 指定所有切割接口面，其中耦合器将被插入。4. 根据接口的扇区复制节点5. 内部面的重复节点6. 根据指定模板构造节点到耦合器这些阶段中的每一个都在[20]和用户手册中有更详细的描述。 在第六阶段结束时，得到的不连续网格显示在图。 1（b）适合分析。有意扩展网格以显示耦合器和复制节点的位置。请注意，节点5出现了四个单独的副本，并且区域A中的所有元素在第一阶段，确定连接到网格中每个节点的单元集在第二阶段，网格中的所有面被分类为边界、界面或面内。在第三阶段期间，如果相关联的接口被指定为具有插入的耦合器，则每个接口刻面被指示为切割。第四阶段是插入算法的显著特征，并规定了接口节点的复制，这依赖于围绕焦点节点的元素扇区的概念，如第2.1节所定义的。 一旦节点周围的扇区被识别，每个扇区被分配一个唯一的节点副本。 该步骤涉及三个拓扑变化：（i）添加全局+T.J. Truster / SoftwareX 7（2018）162165·=·=节点;（ii）节点坐标的复制;以及（iii）元件连接性的更新。参考图1中的网格，b、c和d属于一个扇区，因为耦合器沿着接口（A，B）和接口（B，D）添加。在第五阶段，对所有切割面内小平面。该操作导致在相关区域内的场的完全不连续的内插，例如图1中的区域A。 1（a）.在第六阶段期间，生成耦合器连接性。在此阶段之前已经执行了所有节点复制，使得网格中的所有邻接信息都准备好用于耦合器插入。与元素类似，每种耦合器类型都有一个特定的模板，用于在切割面的两侧布置局部节点。对于内聚带建模，例如（1），通常只需要位于界面上的节点用于计算，因为这些公式是基于位移的。然而，不连续Galerkin（DG）方法，例如（2），通常需要计算沿界面的应力和应变，这涉及使用相邻单元中的所有节点的位移梯度计算。2D线性元件和耦合器模板示例和3D在图2中示出;通过添加中间边缘、中间面和内部节点来自然地处理高阶对应物。体积元素由字母和每个元素的节点数表示，如耦合器上的不同元素类型也受支持，例如关于算法因此，操作计数和内存预计会随着元素和节点线性扩展，这一点通过第3.2节中的实验得到了证实。增加区域的数量对阶段一和阶段二没有影响，但是潜在地增加了切割小平面的数量，因此增加了阶段四中的计算量，这基本上与界面表面积与区域体积的比率有关此外，虽然若干局部计算（例如，集合的交集）可以使用内置MATLAB函数来进行向量化，但是当前算法的实现中的循环的数量对于中等规模的问题（见3.2节），网格修改算法的某些阶段可以并行化，加速。如果耦合器标识符是预先分配的，则第六阶段中的耦合器的构造可以容易地通过接口类型并行化。然而，阶段二到五中的新方面和节点的枚举目前需要是顺序的，以便在编号中不出现间隙，这限制了这些阶段的并行化放松这种编号限制将允许在阶段中并行处理元素和节点组，通过在阶段完成时重新编号结果来牺牲开销。这些权衡值得未来的发展，以评估其对整体执行时间的好处。2.3. 软件功能使 用 DEIP 插 入 耦 合 器 是 通 过 编 写 调 用 程 序 的 各 个 模 块 的MATLAB脚本来完成的唯一需要的输入是与有限元网格的连通性相关的信息，这些信息可以使用本机函数生成或导入，或者由用户提供。随后，可以执行其他模块以执行其他任务，诸如线性有限元分析、网格绘制、导出到可视化软件Paraview，或将包含插入耦合器的修改后的网格导出到支持的网格格式文件。模板脚本包含在程序中，每个模块的参数格式和用法都在用户手册中使用文本和屏幕截图进行记录。DEIP所需的最小数据结构（也在用户手册中记录）是以下三个数组：坐标（node，1：ndm）节点"node“的空间坐标NodesOnElement（elem，1：nen）连接到元素“elem”的节点列表• RegionOnElement（elem）=元素"elem"的区域ID其中Ndm是空间维度的数量（2或3），Nen是每个元素的节点的最大数量。随着内在内聚建模的日益普及，通过支持零厚度界面元素，它已被商业和开源有限元平台所采用。表1提供了具有固有零厚度单元建模能力的平台样本。然而，每个软件通常在ASCII文本或二进制输入文件的组合中具有其自己的细微差别的网格格式和分析命令结构为了覆盖商业和开源软件，DEIP目前支持ABAQUS的.inp文件格式和Gmsh的.msh文件格式[29]，用于导入和导出网格和耦合器数据，因为这两个程序的用户群很大DEIP包还包含示例模型文件，其中包含使用ABAQUS和Code_Aster模拟非线性牵引分离问题的网格文件。可以与其他程序建立连接，例如在表1中，通过开 发 文 本 扫 描 器 来 提 取 三 个 最 小 数 组 Coordinates ，NodesOnElement和RegionOnElement从其他网格文件格式。网格处理步骤的典型顺序概述如下如下所示(1) 导 入 网 格 几 何 形 状 （AbaqusInputReader或GmshInputReader）：如第2.2节所述，该算法需要节点坐标、节点到元素的连接性、元素区域和要切割的界面。这些数据阵列是非常简单的，因此用户可以编写MATLAB插件从标准 有 限 元 输 入 文 件 导 入 这 些 阵 列 此 外 ， 名 为AbaqusInputReader和GmshInputReader的实用程序脚本分别从ABAQUS .inp文件和Gmsh.msh文件导入此数据。(2) 节点复制（DEIPFunction）：通过调用函数DEIPFunction1-5实现沿接口和内部接口执行节点复制的算法阶段1-该函数的输出是具有重复节点的扩展节点坐标阵列、更新的元件连接性以及具有算法阶段6所需的面信息的数据结构。(3) 耦合器插入（InterFunction）：算法的第6阶段定义了沿内接口和内部接口的新耦合器连接性，通过调用函数InterFunction实现 该函数生成内聚区或不连续Galerkin类型的耦合器，如图所示。二、此外，还提供了一个列表，其中耦合器位于每个接口或内面，以便可以分配区域相关的材料属性。(4) 集合更新（UpdateNodeSetFunction）：可以通过调用函数UpdateNodeSetFunction来更新有限元网格中的其他节点数组，例如边界条件，以反映重复的节点。(5) 二次分析（FEA_Program）和可视化（PostParaview）：DEIP提供了一个配套的线性有限元代码FEA_Program，以简化166T.J. Truster / SoftwareX 7（2018）162图二、 线性单元和耦合器模板：（a）三角形二维单元;（b）四边形二维单元;（c）四面体三维单元;（d）六面体三维单元。表1具有零厚度界面元素的有限元软件平台，用于内部内聚区建模以及DEIP与其数据结构的当前兼容性。有限元软件网站许可证DEIP兼容性Abaqushttps://www.3ds.com/products-services/simulia/商业支持ANSYShttps://www.ansys.com/商业不支持代码_Asterhttps://www.code-aster.org/开源通过GmshCOMSOLhttps://www.comsol.com/商业不支持Gmshhttp://gmsh.info/开源支持LS Dynahttp://www.lstc.com/products/ls-dyna商业不支持驼鹿http://mooseframework.org/开源通过GmshWarp3dhttp://www.warp3d.net/开源不支持对生成的非连续单元网格进行验证。本程序具有图2中列出的所有单元类型的版本，用于求解线性各向同性弹性边值问题。此外，还提供了几个绘图功能，用于可视化网格和分析中计算的 位 移 / 应 力 ， 适 用 于 小 到 中 等 尺 寸 的 网 格 。 提 供 脚 本PostParaview，用于将DEIP和FEA_Program的结果导出到VTK文件与PARAVIEW兼容，PARAVIEW是一个流行的开源可视化程序。(6) 导 出 网 格 几 何 形 状 （ AbaqusInputWriter 或GmshInputWriter）：分析节点坐标的合适阵列、节点到元素和耦合器的连接性以及元素区域可以从MATLAB导出回标准有限元输入文件。此外，名为AbaqusInputWriter的实用程序脚本可以将修改后的网格导出到T.J. Truster / SoftwareX 7（2018）162167×=====图3.第三章。 包含多个区域的网格：（a）连续网格;（b）人工分离的不连续网格，指示耦合器位置。表2接口的区域标识符对接口1、3 1、5 1、8 2、6 3、5 4、6 5、6 5、8ABAQUS输入文件。如果为接口耦合器提供了材料属性，并且AbaqusInputReader读取的原始输入文件适合分析，则生成的文件适合或者，GmshInputWriter脚本将网格导出到新的.msh文件，然后需要编写Code_Aster的修改后的命令文件，其中包括接口耦合器材质属性和元素集。3. 说明性实例下面的数值例子验证了所提出的算法产生的分析合适的区域和界面的复杂模式的不连续网格。重点放在被复制的节点和被插入的零厚度界面元件（称为耦合器）上所有的有限元模拟使用的代码，只需要作为输入的传统阵列的节点坐标和元素的连接，这表明该算法3.1. 二维斑贴试验如图3（a）所示，考虑了具有线性三角形和线性四边形单元的矩形4 mm 3 mm域。这些元素被分组到一系列区域中，如图例中的数字和颜色所然后沿着域中的各种接口以及区域6的面内插入耦合器;识别接口的区域对列表在表1中给出。在图3（b）中，为了突出耦合器的位置，域的每个所得部分请注意网格的左下部分，其中五个区域相交于一个节点。该网格的一个激励物理问题是具有多个晶界的多晶材料，这些晶界是裂纹的主要来源，然后使用来自[11，12]的不连续Galerkin公式进行斑贴试验，以测试图3（b）中生成的拓扑的有效性。考虑平面应力条件，见图4。 DG法变形构形上的归一化应力场σxx/Σ b这些区域被赋予材料特性E 100，000 MPa和ν 0。25.沿该区域的左边界的水平位移被规定为零，并且向右边界施加1.2500 MPa的拉伸牵引力。这些边界条件导致右侧面有0.1 mm的水平位移计算的节点应力场如图4所示，为变形后的几何形状。观察到应力场在任何地方都基本恒定;而且，在网格中看不到间隙或重叠因此，我们得出结论，DG耦合器已被正确插入，使补丁测试是满意的变分一致DG方法。使用内聚区（CZ）耦合器代替DG耦合器进行第二次模拟。将内聚区的初始弹性刚度设定为50，000 MPa/mm，以扩大网格中耦合器的存在。因此，在2500 MPa的应力下，耦合器上的预期间隙为0.05 mm。在顶部和右侧域边界上规定了0.1 mm的法向位移，并将对称边界条件应用于底部和左侧表面。图1所示为变形构形上的归一化应力场σxx/Σ b，σ xx/σ b为2500 MPa。 5;变形被放大了两倍。零厚度单元在网格中也被可视化为实体单元。注意，在多个耦合器以多个角度相交的区域中出现了空隙。此外，由于界面的顺应性，区域6和8承载相对较低的应力，使得所施加的应力在整个域中重新分布。然而，没有间断168T.J. Truster / SoftwareX 7（2018）162图五. 用内禀CZ法求变形位形上的归一化应力场σxx/Σ b见图6。用二次四面体单元离散多晶材料域。在区域5和7之间存在位移或应力，因为这些界面没有解耦，并且采用了均匀的弹性特性;我们注意到，在两个图中，节点应力投影都应用在相邻区域内。4和5. 注意，区域5的单个三角形元素已经完全解耦，因为它的所有三个小平面都是切割界面（1，5）和界面（3，5）的一部分。因此，这些结果突出表明，所提出的算法产生的网格是适合CZ建模以及。3.2. 三维可扩展性研究一个数值可扩展性研究进行实际工程意义的几何形状。 该畴是多晶材料的代表性体积元（RVE），它包含100个具有一致晶界的晶粒，如图所示。 六、使用开源软件Neper [30]生成10节点二次四面体单元的网格考虑了三个级别的网格细化，每个网格中的节点和元素的数量列于表3中。因为每个网格当程序重新生成时，元素的数量并不正好以8的倍数增加。每个网格都包含非均匀拓扑;例如，在粗网格上，连接到节点的元素数量在1到40之间变化我们的目标是评估插入耦合器的算法的平均运行时间，类似于[15，17]中的研究。执行两种方案：耦合器仅沿着晶粒之间的界面插入，以及沿着网格中的所有界面和面内刻面在串行台式计算机上使用算法的MATLAB实现进行计时;记录的时间是第2.2节中六个阶段的每个阶段，不包括补片加载时间。节点和元素的随机表3中报告了脚本在每个网格上的总执行时间，这是三次分析的平均值请注意，还使用DG方法[11，12]对该三维几何结构进行了补丁测试，以验证修改后的网格是否适合分析近似线性缩放的计算时间观察相对于网格中的元素的数量或插入的耦合器的数量。在双对数标度上拟合表3耦合器和运行时间之间的线性缩放也通过[15，17]中的自适应耦合器插入来展示此外，所提出的算法的记录时间与其中报告我们希望运行时可以通过使用编译语言而不是脚本语言来稍微改进，如2.3节所讨论的。此外，Fig. 7 b报告模型的输入坐标、元素和区域标识符的内存使用情况，以及DEIP脚本生成的所有数组。我们注意到，从第一阶段到第五阶段的算法中的几个中间阵列是以稀疏格式生成的，并且这些稀疏矩阵的开销在中等网格和细网格之间呈现超线性增长请注意，最终输出网格数据小于初始输入数据库使用量的两倍在一般情况下，插入算法被认为具有最佳的缩放性能的执行时间。4. 影响该软件填补了当前开源和商业有限元代码中的空白，这些代码缺乏零厚度单元或耦合器的网格生成功能，从而使更多的研究人员能够模拟界面现象。表1中列出的几个具有内聚区建模功能的商业和开源有限元平台具有DEIP直接支持的网格格式;其他平台可以通过开发网格格式读取器来连接，以提取算法所需的最少三个阵列。虽然这种离散化技术可以应用于连续体的其他控制方程，但当前材料和结构中感兴趣的应用包括：纤维复合材料中的脱粘[5，21]，多晶金属中晶界处的空腔生长[22，23]，脆性材料中的裂纹扩展[2，19]以及结构中粘合剂粘结的分层[31，32]。代码所需的输入是最小的，并且元件区域的概念是用于描述耦合器位置的几何直观手段，使得代码可以容易地与现有的有限元分析程序结合。事实上，这个独立软件已经从http://clmi.utk.edu/software/下载了几百次。尽管DEIP提供了适用于固体或流体力学中许多类型的界面有限元建模的网格数据，但界面数据结构特别适合于高级微观结构建模。具有多相的材料微观结构因此，这些相之间的界面的机械和失效行为可能是复杂的。近年来的研究开始将材料的相相关性质纳入界面本构模型T.J. Truster / SoftwareX 7（2018）162169表3插入接口和接口耦合器的网格统计信息和运行时间网格元件节点接口耦合器接口时间（s）所有耦合器界面分数总时间（s）粗10,40215,7613,9217.5519,8960.197119.5介质68,71397,72012,11136.2134,4740.0901118罚款539,293742,16742,2202371,067,8850.0395797见图7。多晶畴耦合器插入的复杂性分析：（a）插入脚本执行时间;（b）输入数据和所有内部生成数组的内存使用。晶间裂纹、位错传递、摩擦滑动等。[33具有用于许多类型的接口耦合器的通用代码以及元件类型转换脚本还提供了用于比较不同配方和元件类型的平台，以更好地理解其准确性和稳定性属性。例如，先前的研究发现，具有特定数值积分规则的低阶单元为断裂力学分析提供了更好的稳定性[13，36，37]。其他研究已经证实了内聚区和混合不连续Galerkin公式用于模拟裂纹扩展的效率[2，19]。5. 结论不连续单元插入程序（DEIP）采用拓扑算法将零厚度界面单元（称为该工具箱非常适合于几类问题，包括：（i）内在的凝聚区建模，不连续Galerkin公式，或其他接口方法;（ii）线性和高阶有限元的各种形状;和（iii）二维和三维域。该方法的关键特征是使用域中元素的DEIP还提供了其他模块，以实现有限元分析和可视化中使用的数据输入/输出。这些模块包括一个独立的线性有限元分析代码，读写脚本的ABAQUS和Gmsh输入MATLAB的网格和轮廓绘图功能，以及PARAVIEW的提供了示例文件和图形这些模块以及软件中的错误检查使用户能够为复杂的三维模型生成适合分析的不连续网格。对二维和三维问题的数值试验表明，该算法与单元数成线性关系，并能产生正确的节点复制和耦合器插入模式，以保持区域内所需的连续特征。确认T. Truster通过橡树岭国家实验室的DE-AC 05 - 000 R22725项目得到了一个研究机构的支持感谢您的支持.引用[1] 潘多尔菲A，奥尔蒂斯M。三维碎片的实体建模方面。工程计算1998;14（4）：287-308.[2] [10]杨文，李文，李文.基于混合、不连续伽辽金、内聚元方法的可扩展三维断裂和破碎算法。 计算方法应用机 械 工程 2011;200（1-4）：326-44.[3] Alfano G，Crisfield MA.层压复合材料分层分析的有限元界面模型：力学和计算问题。国际数字方法工程杂志2001;50（7）：1701-36.[4] Raghavan P，Ghosh S.单向复合材料界面脱粘的连续损伤力学模型。机械材料2005;37：955-79.170T.J. Truster / SoftwareX 7（2018）162[5] 作 者： J.纤 维复 合材 料界 面 损伤 建模 的间 断/连 续Galerkin方 法。 ComputMech2013;52（3）：499-514.[6] Carrier B，Granet S.渗透性介质水力压裂问题的软聚带模型数值模拟。工程机械2012;79：312 - 28。[7] Sarris E，Papanastasiou P.粘性多孔弹塑性连续介质中流体驱动裂缝的数值模拟。Int J Numer Anal Methods Geomech2013;37（12）：1822-46.[8] 陈文辉，陈文辉，陈文辉.金属陶瓷界面结合区模型应用中的若干问题。国际固体结构杂志2002;39（10）：2827-55.[9] ArnoldDN，Brezzi F，Cockburn B，Marini LD.椭圆问题间断Galerkin方法的统一分析。SIAM J Numer Anal2002;39（5）：1749 - 79.[10] LiuR，Wheeler MF，Dawson CN.弹性力学问题一族间断Galerkin方法的三维数值实现。Comput Struct2009;87（3）：141-50.[11] 作 者 ： J. A Primal interface formulation for coupling multiplePDE ： Aconsistent derivation via the Variational Multiscale method. 计算方法应用机械工程2014;268：194-224.[12] MasudA，Truster TJ，Bergman LA. 界面耦合和摩擦接触建模的统一公式及嵌入误差估计。Internat JNumer Methods Engrg2012;92（2）：141-77.[13] DayRA，Potts DM. 零厚度界面元数值稳定性及其应用。国际岩土工程数值分析方法杂志1994;18（10）：689-708.[14] 潘多尔菲A，奥尔蒂斯M。三维碎裂模拟的有效自适应程序。工程计算2002;18（2）：148-59.[15] 放 大 图 片 作 者 ： A. 单纯形 有 限 元 网 格 的 断 裂 与 破 碎 。 InternatJ NumerMethods Engrg 2008;73（11）：1547-70.[16] Ortiz M，Pandolfi A.三维裂纹扩展分析的非可逆变形粘性单元。国际数字方法工程杂志1999;44（9）：1267-82.[17] PaulinoGH，Celes W，Espinha R，Zhang Z.基于拓扑的有限元网格自适应内聚单元插入框架。 EngComput2007;24（1）：59-78.[18] Nguyen VP. 一 个 开 源 程 序 ， 用 于 生 成 零 厚 度 内 聚 界 面 元 素 。 Adv EngSoftw2014;74：27-39.[19] NguyenVP.不连续Galerkin/非本征内聚区模型：计算断裂力学中的实施注意事项和应用。工程机械2014;128：37-68。[20] 信任TJ关于三维网格中的界面单元插入。工程机械2016;153：171-4.[21] 放大图片作者：Chen P，T.多组分材料中弱和强不连续性有限应变下的界面稳定性。计算机方法应用机械工程2017.[22] Ma R，Truster TJ，Puplampu SB，Penumadu D.用晶体塑性有限元法研究了5083铝合金在火灾下的力学性能退化. Int J Solids Struct2017.[23] 纳西夫 ·托拉斯 特·TJ用李代数求解梯 度晶体塑性的变分 投影法。 InternatJNumer Methods Engrg 2017;110（4）：303-32.[24] Celes W，Paulino GH，Espinha R.一种用于有限元网格表示的紧凑的基于邻接的拓扑数据结构。Internat J Numer MethodsEngrg2005;64（11）：1529-56.[25] 尼德曼河夹杂脱粘引起空洞形核的连续介质模型J. Appl. Mech. 1987;54（3）：525.[26] Park K，Paulino GH.基于PPR势的内聚模型在ABAQUS中的计算实现：教育视角. 工程机械2012;93：239-62。[27] 放大图片作者：J.二阶椭圆问题的间断Galerkin、混合Galerkin和连续Galerkin方法的统一杂交。SIAM J Numer Anal2009;47（2）：1319-65.[28] [10]杨文辉，陈文辉. 协议。II库，8.5版 JNumer Math2017;25（3）：137-45.[29] Geuzaine C，Remacle J-F. Gmsh：一个内置预处理和后处理设施的三维有限元网格生成器。Internat J Numer Methods Engrg2009;79（11）：1309-31.[30] 放大图片QueyR，Dawson PR，Barbe F.大规模三维随机多晶体有限元法：生成，网格化和重新网格化。 计算方法应用机械工程2011;200（17-20）：1729-45.[31] 迪尔·T基于惩罚的粘结区有限元法，第二部分：环氧树脂粘结铝带的非弹性剥离。Int J Adhes Adhes2008;28（4）：256-65.[32] Kowalewski，Gaznski M. 用粘聚力单元法确定砖墙的破坏模式。ProcediaEng2015;111：454-61.[33] Benedetti I，Aliabadi MH.多晶材料晶间退化和破坏的三维粘聚-摩擦晶界微观力学模型。计算方法应用机械工程2013;265：36 - 62.[34] Gottschalk D，McBride A，Reddy BD，Javili A，Wriggers P，HirschbergerCB. 考虑晶粒取向差和晶界取向的晶界模型的计算和理论方面。 ComputMaterSci2016;111：443-59.[35] AnnavarapuC，Settgast RR，Johnson SM，Fu P，Herbold EB.摩擦表面上小滑动接触的加权Nitsche稳定化方法。 计算方法应用机械工程2015;283：763-81。[36] SchellekensJCJ，De Borst R.界面元的数值积分国际数字方法工程杂志1993;36（1）：43-66.[37] 张文，张文，等.一种新的三维界面有限元法在 复合材料中的应用.北京：清华大学出版社，2001. 国际固体结构杂志2004;41（11-12）：2977-93.