多種流形覆蓋方式耦合的數值流形法及其應用

摘要：為了更好地利用數值流形法（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.

（編輯 陳移峰）