多种流形覆盖方式耦合的数值流形法及其应用

武鑫 郭璇 康歌 甯尤军

摘要：为了更好地利用数值流形法（NMM）模拟连续和非连续问题，根据NMM数学覆盖选择的灵活性，将有限元网格作为一种数学覆盖方式生成流形覆盖系统以避免细小流形单元的出现，并提出了将单个材料体描述为一个流形单元的独立覆盖方法。运用多种覆盖方式相耦合的方法模拟了混凝土的细观拉压模型，以及一个复杂结构岩质边坡模型，模拟结果较好地再现了相应的力学变形破坏过程。多种覆盖方式的耦合使得构建的NMM模拟模型更加合理；独立覆盖方式的使用降低了构建离散体系NMM模型的复杂程度，并减少了模型中的单元数量以及不连续面上的接触数量，提高了模拟计算效率。

关键词：数值流形法（NMM）；独立覆盖；有限覆盖；耦合使用；离散体系

中图分类号：O347.7 文献标志码：A文章编号：1000-582X（2018）05-060-08

Abstract： To better use the numerical manifold method （NMM） to simulate complex continuous and discontinuous deformation problems， for the flexibility of the choice of mathematical covers in the NMM， we adopt the finite element mesh to generate NMM mathematical covers so as to avoid the appearance of too-small manifold elements， and put forward an independent cover method， in which a single material domain is described as one manifold element. Different cover types are coupled in the simulations of the tensile and compressive concrete meso-scale models， as well as a rock slope model with complex structures. The corresponding mechanical deformation and failure processes are well reproduced by NMM simulations. The coupling use of different cover types makes the generated NMM simulation models more reasonable. The independent cover method decreases the complexity to construct NMM simulation models for discrete systems， and reduces the numbers of manifold elements and the contact numbers along discontinuity surfaces， which helps to improve the NMM computational efficiency.

Keywords： numerical manifold method （NMM）；independent cover；finite cover；coupling use；discrete system

混凝土、岩石是典型的非均匀、非连续工程及地质材料，其变形、破坏是一个复杂的连续、非连续固体介质力学行为问题。数值模拟是研究这类复杂力学问题的重要手段。基于连续理论的有限元法（FEM），在模拟材料不连续变形和运动力学行为时能力存在不足。在基于非连续理论的方法中，离散元法（DEM）运用显式有限差分法模拟多体不连续系统的大位移和旋转过程[1]；非连续变形分析（DDA）方法则基于最小势能原理，较好地解决了多材料体的大变形和大位移问题。对于岩石力学行为的模拟，DEM和DDA方法适用于对离散岩体结构的稳定状态和变形模式及过程的模拟，但对于连续体内裂纹的产生和扩展模拟则存在一定不足[2]。显然，2种数值方法对于混凝土这种颗粒基体材料的细观变形破坏模拟也同样存在难度。连续非连续数值方法的耦合，如FEM-DEM[3]，以及连续非连续理论在数值方法中的统一，如数值流形法（NMM）[4]，则是更好的解决途径。

数值流形法（NMM， numerical manifold method）通过2套独立的覆盖系统很好地将连續与非连续问题统一在一个理论框架内进行描述，将FEM、DDA和解析法统一起来。覆盖系统是NMM的重要组成部分，对此已有了较多的研究。蔡永昌等[5]论述了二维弹性问题里使用矩形覆盖流形方法的理论和实施技术，并给出了其覆盖系统的全自动生成算法；陈刚等[6]将有限元网格作为生成流形覆盖系统的数学网格，利用有向图的关联矩阵结合块体搜索算法，形成最终的流形法覆盖系统；李海枫等[7]研究了三维流形单元的生成方法；蔡永昌等[8]提出了一种完全独立覆盖的高阶NMM，在独立覆盖之间假设几何厚度可以忽略的连接弹簧，以解决高阶流形法的线性相关性问题和物理覆盖系统的生成复杂性问题；苏海东等[9]研究了矩形独立覆盖及其覆盖加密，在独立覆盖中使用解析级数，并提出了任意形状的独立覆盖形式。NMM在类混凝土材料的细观变形破坏，以及岩质边坡失稳等岩石力学与工程问题的研究中都已得到了一定的应用。黄涛[10]和赵妍[11]分别采用NMM模拟研究了混凝土以及高聚物粘结炸药这两类颗粒基体材料的细观变形破坏问题；张国新等[2]验证了采用NMM可以正确计算边坡倾倒安全系数和模拟边坡倾倒破坏过程；黄涛等[12]应用NMM对岩石双孔爆破过程进行了模拟；Ning和An等[13-14]采用基于莫尔库仑开裂准则的NMM模拟了岩质边坡中的裂纹扩展和边坡破坏全过程。

根据NMM中数学覆盖和物理覆盖的定义，提出了将单个材料体描述为单个流形单元的独立流形覆盖方法，并与传统的规则网格有限覆盖以及由有限元网格生成的有限覆盖相结合，以便对包含离散体系的复杂材料结构的力学行为进行模拟。通过混凝土细观变形破坏行为模拟和复杂结构岩质边坡失稳破坏模拟算例，验证和讨论了这一数值模拟途径的可行性及优势。

1 数值流形法的有限覆盖与独立覆盖

1.1 数值流形法的基本思想

采用2套独立的覆盖系统：数学覆盖系统和物理覆盖系统[3]。数学覆盖（MC， mathematical cover）系统由许多部分重叠的覆盖片组成，它覆盖住整个物理区域。在二维NMM中，一般多采用规则三角形网格生成数学覆盖系统。如图1（a）所示，当采用三角形网格生成数学覆盖系统时，对于网格中的任一节点，含有这一节点的所有三角形组成1个数学覆盖。在图1（a）中，1-6号三角形组成的正六边形就是该节点处的数学覆盖。

物理覆盖（PC， physical cover）是数学覆盖和物理区域的交集，当将数学覆盖置于物理区域上时，数学覆盖与物理区域重叠的部分就形成物理覆盖系统。如果材料边界或内部不连续面将1个数学覆盖切割成多个完全独立的区域，则位于材料域内的每一独立区域就是1个物理覆盖。在图1（a）中，以节点a、b为中心的六边形数学覆盖Ma、Mb完全位于物理区域内，并且均没有被不连续面切割成完全独立的区域，从而使得Ma、Mb分别成为1个物理覆盖，记为P1a、P1b；以节点c为中心的数学覆盖Mc只有一部分与物理区域相重叠，只有位于物理区域内的这一部分才形成1个物理覆盖，记为P1c，同理，以节点e为中心的数学覆盖Me处形成1个物理覆盖，记为P1e；而对于以节点d为中心的数学覆盖Md而言，其位于物理区域内的部分被物理区域中的不连续面切割成2个完全独立的区域，所以，在节点d处形成2个物理覆盖，记为P1d、P2d。

对于由以上方法产生的物理覆盖，他们也是相互部分重叠的，这些相互重叠的物理覆盖所形成的交集则构成流形单元（ME， manifold element）。如图1所示，图1（a）中的物理覆盖P1a、P1b和P1c相交形成流形单元e1，物理覆盖P1b、P1d、P1e相交形成流形单元e2，物理覆盖P1b、P2d、P1e相交形成流形单元e3。显然，每个流形单元是3个物理覆盖的交集。

1.2 有限覆盖与独立覆盖

如图1所示的由规则三角形网格生成数学覆盖，是NMM广泛采用的一种覆盖形式。由于不需要网格边与物理边界相重合，特别是对于复杂形态的物理区域，前处理变得易于实现，但同时也会在物理区域的边界产生一些面积极小或细长的流形单元，如图2（a）所示，从而对计算性能造成一定影响。生成NMM覆盖系统的另外一种方式是采用传统有限元网格作为数学覆盖，如图2（b）所示，当采用三角形有限元网格时，对于其中任一节点，含有这一节点的所有三角形单元形成1个数学覆盖，进而生成相应的物理覆盖。对于任意一个有限单元，其3个顶点处的物理覆盖相互重叠形成的流形单元就是该有限单元本身。这时NMM方法实质上退化为了有限单元法，只是在NMM理论框架下，材料开裂破坏的描述较有限元更加易于实现。无论是采用规则三角形网格还是采用有限元非规则网格生成的流形覆盖，都可以称为NMM的有限覆盖。采用有限覆盖可以很方便地描述连续材料体的变形和开裂破坏行为。

对于离散体系而言，主要关心离散体的位移和相互之间的作用，而离散体自身的变形不重要时，采用如上所述的有限覆盖则不具必要性和计算经济性。这里提出独立覆盖的概念，即每一个离散体物理区域采用独立的覆盖进行描述，每个块体形成一个流形单元，并拥有属于自身的独立的位移函数，实现离散体间位移不连续的描述。这时，NMM实质上退化为了DDA方法，只是采用了不同的具体未知量和位移函数。值得注意的是，这里提出的独立覆盖的概念和文献[8-9]中独立覆盖的概念是不一样的。

当采用三角形网格时，独立覆盖的构造方法是生成1个包络单个离散体的三角形，并以该三角形的3个顶点为数学覆盖的中心，产生3个数学覆盖，进而生成相应3个物理覆盖。如图3（a）所示，1个四边形离散体物理区域被覆盖其上的Δijk完全包含。根据数学覆盖的概念，节点i、j、k处的数学覆盖分别为以节点i、j、k为中心的3个六边形， 记为Mi、Mj和Mk。它们与物理区域的交集分别构成3个四边形物理覆盖Pi、Pj、Pk（区域均为物理区域本身），这3个四边形物理覆盖相互重叠形成1个四边形流形单元，即该物理区域本身。这样以来，单个离散体物理区域自身形成1个流形单元。对于如图3（b）所示的离散体系统，其中每一个四边形离散体均拥有属于自身的独立覆盖，各自形成独立的流形单元。

以上独立覆盖生成方法略去了数学覆盖被物理区域再分割的过程，简化了物理覆盖与流形单元的生成过程。在保证相互对应的前提下，覆盖系统和流形单元可以分开建立。具体可以先建立离散体系统模型，每个离散體形成1个流形单元，进而根据离散体的几何和坐标信息，对每个离散体引入1个三角形并形成相应的数学与物理覆盖。如图3（b）所示，由此生成的覆盖系统，保证了属于任何2个离散体的物理覆盖间的独立性。

对于具体问题的模拟，在覆盖系统的生成过程中，可以根据问题的几何及物理特征，将由规则三角形或有限元网格生成的有限覆盖，与独立覆盖进行耦合使用，得到更加合理的NMM模拟模型。在文中，有限覆盖和独立覆盖采用相同的覆盖位移函数，自然地将各种覆盖方式耦合到了统一的NMM数学描述框架下，具有简洁性和实用性。

2 不同流形覆盖方式的耦合应用

2.1 混凝土细观模型的单轴拉压模拟

混凝土在细观尺度上可看成是由骨料、水泥砂浆及其交界面组成的颗粒基体复合材料。其中，骨料与水泥砂浆的交界面是混凝土材料内部最薄弱的相，在外荷载作用下，微裂纹往往率先在交界面上开始萌生，并逐渐演化、发展和贯穿，最终形成宏观裂纹导致混凝土材料的破坏。混凝土拉伸破坏后的断面形态表明，绝大部分破坏面是粗骨料与水泥砂浆的界面粘结破坏（拉脱），而小部分是骨料间的水泥砂浆被拉断。即使是在压缩载荷作用下的破坏，一般也是以水泥砂浆的开裂和界面开裂为主，较少发生骨料破坏的情况。因此，对于混凝土这种颗粒基体复合材料的细观模拟，骨料的变形和破坏一般可以认为并不重要。

在混凝土细观模型的NMM模拟中，可以采用独立覆盖描述骨料，采用有限覆盖描述砂浆，界面则通过粘接强度加以考虑。如图4（a）所示，其中的骨料、加载板以及测量板采用了独立覆盖，砂浆则采用了有限元网格来生成覆盖系统，以避免采用规则三角形网格生成覆盖时在砂浆骨料界面处产生细长或极小尺寸流形单元。测量板位于加载板和细观模型之间，因采用独立覆盖而具有常应力，它在竖直方向上的应力反映了混凝土细观模型在竖向方向上的平均应力。为加以对比，在图4（b）的模型中，对骨料也采用有限元网格进行了覆盖。

在方件上下两端均通过加载板进行准静态线性拉伸或压缩位移加载，加载速率为2 mm/s。假定骨料与水泥砂浆为线弹性材料，各类材料参数如表1所示。NMM模拟中的最大步位移比、时间步长和接触弹簧刚度分别取为5×10-4，1×10-5 s和600 GPa。

图5（a）与图5（b）所示为单轴拉伸的细观模拟结果。骨料采用2种不同覆盖方式，模型中主贯通裂纹均产生于模型中下部，并伴有多处分散裂纹，其中，砂浆骨料界面开裂占整个开裂面的较大部分。2种覆盖条件下模型的细观破坏形态极为相似。图6（a）给出了由测量板得到的模型在竖直方向上的平均应力变化曲线。曲线在初期表现为线性，随着加载位移的增大，界面的逐步脱粘和砂浆的逐步开裂导致软化现象的出现。在主裂纹逐渐贯通的过程中，曲线开始出现明显下降。由于材料内部结构的不均匀性，曲线在下降过程中出现了一定的起伏现象。可以看出，骨料分别采用有限覆盖和独立覆盖时，模型竖直方向的应力位移曲线在前期极为近似，只是在后期的起伏阶段才产生了较明显差异。

图5（c）与图5（d）所示为单轴压缩的模拟结果。在压缩载荷作用下，细观结构的非均匀性导致大量分散裂纹的产生，没有形成明显优势裂纹带，但单条分散裂纹基本都为倾斜方向的剪切裂纹。可以看出，在骨料采用2种不同覆盖的条件下，模型中细观破坏形态仍极为相似。图6（b）所示的模型竖直方向平均应力变化曲线表明，初期模型整体响应呈线性，随着分散裂纹的逐步产生，软化现象逐步显现。对比骨料采用2种不同覆盖方式得到的应力曲线，可以发现它们仍非常接近。

综上所述，对于混凝土这样的颗粒基体材料，当骨料颗粒分别采用有限覆盖和独立覆盖时，得到了基本一致的模拟结果。但当骨料采用独立覆盖时，覆盖系统的生成更加简单，并且由于界面处接触数量的减少，以及整个模型中单元数目的减少，模拟计算效率也得到了一定的提高。

2.2 复杂岩质边坡的失稳破坏模拟

岩体是由连续岩石介质和不连续面构成的天然地质材料，岩质边坡复杂的组成结构对其破坏形式有着直接重要影响。根据Goodman[15-16]关于岩质结构边坡失稳的研究，建立如图7（a）所示的由连续岩石介质与离散岩块构成的边坡模型。在由连续岩石介质和少量不连续面组成的节理欠发育岩体中，连续岩石介质的开裂破坏是非常重要的破坏形式。而在包含大量不连续面的节理较充分的发育岩体中，岩体可看做是由独立岩块构成的离散岩块系统，岩体的变形破坏主要沿着已有的不连续面进行。下面通过NMM方法模拟该边坡在重力作用下的失稳破坏过程。基于前文所提及的有限覆盖和独立覆盖2种不同覆盖方式，边坡中节理欠发育岩体部分采用规则三角形网格进行有限覆盖以描述其变形破坏，边坡中节理充分发育岩体部分采用独立覆盖，以便于描述离散岩块系统的变形破坏及运动。同时，用作边界约束条件的模型边框也采用独立覆盖。NMM网格模型如图7（b）所示。

模拟中，节理充分发育部分岩体（记为岩层）的材料参数为：密度2 500 kg/m3，弹性模量30 GPa，泊松比0.25。节理欠发育部分岩体（记为岩块）的材料参数为：密度2 700 kg/m3，弹性模量30 GPa，泊松比0.25。岩层的强度参数取为：内摩擦角27°，内聚力0.5 MPa，抗拉强度0.5 MPa，岩块部分不考虑开裂破坏。模型中不连续面的强度参数均取为：摩擦角10°，内聚力0.05 MPa，抗拉强度0.05 MPa。这里取相对一般实际更小的强度参数，是为了使得边坡在重力作用下即发生破坏失稳。

图7（c）和图7（d）给出了边坡破坏过程的NMM模拟图像。在自重和上部岩体的压载作用下，节理发育良好的边坡下部岩体沿既有不连续面发生破坏，岩块系统发生松动。下部离散岩块系统承载能力的降低，导致上部岩层的不稳定并沿竖向贯通节理面发生下错，岩层中既有的非贯通节理面进而发生扩展，将岩层逐步折断。随后，下部离散岩块被向右挤出并向前方倾倒，上部破坏岩层下端与下部离散岩块一起向右侧运动，整个上部岩层则发生向后倾倒并下滑。这种破坏称为滑动倾倒的“二次倾倒模式”[15]，当上部岩层滑动趋势被下端岩块阻止时发生。以上模拟结果表明，采用有限覆盖和独立覆盖耦合的方式能很好地模拟复杂岩质边坡的失稳破坏问题。

3 結 论

基于NMM有限覆盖理论，采用有限元网格作为一种数学覆盖方式生成NMM流形覆盖系统，避免采用规则网格生成流形覆盖系统时可能产生的细长或极小尺寸流形单元，并提出了将单个材料体描述为单个流形单元的独立覆盖方法。通过多种流形覆盖方式相耦合，模拟了混凝土方件的细观变形破坏，以及复杂结构岩石边坡的失稳破坏过程。结果表明，独立覆盖方式的采用，降低了模型的生成过程以及最终生成模型的复杂程度，减少了计算过程中的接触数量以及流形单元数量，提高了计算效率。独立覆盖与有限覆盖的耦合使用，合理地描述了混凝土的细观结构和岩质边坡的复杂构造，NMM模拟较好地再现了相应的物理力学过程。文中提出的独立覆盖方法，适用于单体变形不重要的离散体的模拟，也可用于模型中加载板、约束边框等单体的描述。

参考文献：

[1] 刘凯欣，高凌天.离散元法研究的评述[J].力学进展，2003，33（4）：483-490.

LIU Kaixin， GAO Lingtian. A review on the discrete element method[J]. Advances in Mechanics， 2003， 33（4）： 483-490.（in Chinese）

[2] 张国新，赵妍，石根华，等.模拟岩石边坡倾倒破坏的数值流形法[J].岩土工程学报，2007，29（6）：800-805.

ZHANG Guoxin， ZHAO Yan， SHI Genhua， et al. Toppling failure simulation of rock slopes by numerical manifold method[J]. Chinese Journal of Geotechnical Engineering， 2007， 29（6）： 800-805. （in Chinese）

[3] Mahabadi O K， Lisjak A， Munjiza A， et al. Y-Geo： New combined finite-discrete element numerical code for geomechanical applications[J]. International Journal of Geomechanics， 2012， 12（6）： 676-688.

[4] Ma G W， An X M， He L. The numerical manifold method： a review[J]. International Journal of Computational Methods， 2010， 7（1）： 1-32.

[5] 蔡永昌，張湘伟.流形方法的矩形覆盖系统及其全自动生成算法[J].重庆大学学报（自然科学版），2001，24（1）：42-46.

CAI Yongchang， ZHANG Xiangwei. Rectangular cover system of manifold method and its auto mesh algorithm[J]. Journal of Chongqing University （Natural Science Edition）， 2001，24（1）： 42-46. （in Chinese）

[6] 陈刚，刘佑荣.流形元覆盖系统的有向图遍历生成算法研究[J].岩石力学与工程学报，2003，22（5）：711-716.

CHEN Gang， LIU Yourong. Generation of cover system for numerical manifold in travel theory of the oriented graph[J]. Chinese Journal of Rock Mechanics and Engineering， 2003， 22（5）： 711-716. （in Chinese）

[7] 李海枫，张国新，石根华，等.流形切割及有限元网格覆盖下的三维流形单元生成[J].岩石力学与工程学报，2010，29（4）：731-742.

LI Haifeng， ZHANG Guoxin， SHI Genhua， et al. Manifold cut and generation of three-dimensional manifold element under FE mesh cover[J]. Chinese Journal of Rock Mechanics and Engineering， 2010， 29（4）： 731-742.（in Chinese）

[8] 蔡永昌，刘高扬.基于独立覆盖的高阶流形方法[J].同济大学学报（自然科学版），2015，43（12）：1794-1800.

CAI Yongchang， LIU Gaoyang. High-order manifold method with independent covers[J]. Journal of Tongji University （Natural Science）， 2015， 43（12）： 1794-1800.（in Chinese）

[9] 苏海东，颉志强，龚亚琦，等.基于独立覆盖的流形法的收敛性及覆盖网格特性[J].长江科学院院报，2016，33（2）：131-136.

SU Haidong， XIE Zhiqiang， GONG Yaqi， et al. Characteristics of convergence and cover mesh in numerical manifold method based on independent covers[J]. Journal of Yangtze River Scientific Research Institute， 2016， 33（2）： 131-136.（in Chinese）

[10] 黄涛.聚合物粘结炸药损伤与断裂的流形元数值模拟[D].北京：北京理工大学，2006.

HUANG Tao. Numerical simulation of damage and fracture of polymer explosive by manifold method[D]. Beijing： Beijing Institute of Technology， 2006. （in Chinese）

[11] 赵妍.基于流形元法的裂缝追踪技术在结构破坏模拟中的应用研究[D].北京：中国水利水电科学研究院，2008.

ZHAO Yan. Crack tracking techniques based on the manifold method and its applications in structure failure simula-tions[D]. Beijing： China Institute of Water Resources and Hydropower Research， 2008. （in Chinese）

[12] 黄涛，陈鹏万，张国新，等.岩石双孔爆破过程的流形元法模拟[J].爆炸与冲击，2006，26（5）：434-440.

HUANG Tao， CHEN Pengwan， ZHANG Guoxin， et al. Numerical simulation of two hole blasting using numerical manifold method[J]. Explosion and Shock Waves， 2006， 26（5）： 434-440. （in Chinese）

[13] Ning Y J， An X M， Ma G W. Footwall slope stability analysis with the numerical manifold method[J]. International Journal of Rock Mechanics and Mining Sciences， 2011， 48（6）： 964-975.

[14] An X M， Ning Y J， Ma G W， et al. Modeling progressive failures in rock slopes with non-persistent joints using the numerical manifold method[J]. International Journal for Numerical and Analytical Methods in Geomechanics， 2014， 38（7）： 679-701.

[15] Goodman R E， Kieffer D S. Behavior of rock in slopes[J]. Journal of Geotechnical and Geoenvironmental Engineering， 2000， 126（8）： 675-684.

[16] Goodman R E . Toppling-A fundamental failure mode in discontinuous materials—description and analysis[C]∥[S.l.]Geo-Congress 2013： Stability and Performance of Slopes and Embankments III. 2013： 2338-2368.

（編辑 陈移峰）