數值流形程序開發設計及其應用

焦　健

(北京城建設計發展集團股份有限公司，北京　100037)

?

數值流形程序開發設計及其應用

焦健

(北京城建設計發展集團股份有限公司，北京100037)

針對數值流形程序開發面臨覆蓋系統剖分算法、基于覆蓋的平衡方程的建立、接觸問題處理等幾個難題，提出了具體的解決對策，完成了數值流形程序的編制，并通過算例驗證了該程序的可靠性，展現了該程序在處理不連續變形計算方面的獨特優勢。

數值流形方法，不連續變形計算，程序開發

數值流形方法(NMM)誕生于巖石力學領域，是由石根華博士在創立非連續變形分析[1](DDA)之后提出的一種新型的數值計算方法。它通過采用連續和非連續覆蓋函數的辦法，將有限元與DDA整合到一種計算模式下[2]。

自數值流形方法提出以來一直頗受關注與好評，但石根華博士的數值流形程序并未在國內作大力推廣，而科研工作者獨自將其編程實現又面臨諸多困難[3,4]，因此數值流形方法在國內的研究和應用還僅限于為數不多幾個科研機構。在前人工作的基礎上[5]，完成了數值流形程序的開發設計，并成功地用于幾個實際問題的解決，取得了滿意的結果。本文重點介紹數值流形程序編制中幾個難點問題的解決方法。

1　覆蓋系統自動剖分算法

數值流形方法中采用兩套獨立的網格，物理網格和數學網格。物理網格為材料固有邊界和內部裂紋；數學網格可以由常規的網格，如有限元網格、規則的格子或級數收斂域轉化而得到。采用三角形或四邊形有限元網格，可以方便地取節點的形函數作為覆蓋函數的權函數，是目前主要的數學網格形式。二者相比，四邊形有限元網格具有拓撲運算簡單、同階覆蓋函數條件下精度更高的優勢。此外，數值流形方法并不要求數學網格與求解域完全吻合，而是能夠覆蓋整個求解域即可。所以，即使在邊界不規則處也無需調整網格形狀以適應邊界條件?；谏鲜隹紤]，本文采用矩形數學網格，在此基礎上完成劃分流形單元和覆蓋系統自動剖分的工作。

如圖1所示，粗線為材料物理網格，細線為數學網格(圖中一個基本數學網格定義為一個數學單元，如數學單元1-2-5-4)。基于數學網格節點定義了數學覆蓋：包含某節點的所有數學單元的并集稱為一個數學覆蓋，如節點5所對應的數學覆蓋為矩形區域1-3-9-7。一個數學覆蓋可以被物理網格剖分為若干物理覆蓋，如數學覆蓋1-3-9-7被物理網格剖分出六個獨立的子區域，其中只有13-15-16-26-22和22-26-19-20兩個子區域位于材料內部，成為兩個物理覆蓋。流形單元是數值流形方法的最小積分單位，它是若干個物理覆蓋的交集。在矩形數學網格的條件下，流形單元必定是四個物理覆蓋的交集。流形單元上的位移函數是這四個物理覆蓋上的覆蓋函數加權之和。編制數值流形程序的第一個困難就是在兩套網格的疊加下搜索記錄所有流形單元，以及每個流形單元所屬的物理覆蓋。

本文采用的覆蓋剖分算法如下：循環每一條物理網格與數學網格線做求交運算，交點將每條物理網格劃分為若干段，記錄每一段物理網格。同時，數學單元在物理網格的剖分下出現若干個子區域。若子區域位于材料內部，則構成一個流形單元，并為流形單元添加一個成員變量，記錄其所在的數學單元的四個節點；位于材料外部則不予考慮。如圖1中數學單元1-2-5-4被剖分成23-5-24-22，21-22-24和1-2-23-21-4三個子區域，只有前兩個子區域位于材料內部，屬于流形單元，最后一個子區域無需處理。整個系統完成全部流形單元的搜索后如圖2所示。

再對每一個數學網格節點進行循環，獲得各個節點所對應的物理覆蓋。算法如下：針對當前數學網格節點，搜索所有包含當前節點的流形單元并求其并集，構成該節點的數學覆蓋。然后，令物理網格對該數學覆蓋進行剖面，剖分處的每一個單獨區域即構成一個物理覆蓋，將其編號記錄。

獲得物理覆蓋之后，便可以搜索每個流形單元所屬的物理覆蓋。如前述，每一個流形單元均記錄了其所在數學單元的四個節點，這四個節點均對應若干個物理覆蓋，依次對該流形單元與四個節點的物理覆蓋進行分析，判斷流形單元完全包含于哪個物理覆蓋中，共可以找到四個物理覆蓋，這四個物理覆蓋即為該流形單元所屬的物理覆蓋。

2　覆蓋函數和權函數

取覆蓋函數為常函數，矩形有限元的形函數作為權函數，即：

(1)

(2)

其中，ui，vi分別為x向，y向的覆蓋函數；We(i)為流形單元第i(i=1～4)個物理覆蓋的權函數，采用的是流形單元所在矩形數學單元的形函數[9]。

則流形單元上的總體位移函數可表示為:

(3)

其中，

Te=[Te(1)(x,y),Te(2)(x,y),Te(3)(x,y),Te(4)(x,y)]

(4)

(5)

De={De(1),De(2),De(3),De(4)}T

(6)

其中，De(i)為流形單元第i(i=1～4)個物理覆蓋上的常覆蓋函數{ui,vi}。

3　平衡方程的建立

數值流形方法依據最小勢能原理建立平衡方程。總勢能包括：應變勢能、初應力勢能、荷載勢能、慣性力勢能、彈簧勢能和接觸摩擦力勢能。對它們分項求導使勢能極小化則可得總體平衡方程:

(7)

其中,總剛度矩陣K由單元剛度矩陣、固定點矩陣、慣性力矩陣和接觸彈簧矩陣迭加得到，各矩陣的具體形式及計算某些矩陣元素涉及到的單純形積分問題見文獻[2]。荷載矩陣由F初應力矩陣、點荷載矩陣、分布荷載矩陣、體荷載矩陣、速度矩陣以及接觸摩擦力矩陣迭加組成。文獻[2]給出了除分布荷載矩陣外其他矩陣的具體形式；而在實際問題中，分布荷載也是很常見的荷載形式。本文推導建立了分布荷載矩陣。

作用于某條邊上的分布荷載通常會跨越多個流形單元，首先搜索獲得荷載通過的所有流形單元，再逐個為這些單元添加分布荷載矩陣。以圖3中單元A-B-C為例推導分布荷載矩陣。

設A點(x1,y1)和B點(x2,y2)處的荷載集度分別為(fx1,fy1)，(fx2,fy2)。AB邊上任一點的坐標和荷載集度可表示為:

(8)

(9)

繼而，分布荷載在該單元上的勢能有如下表達式:

(10)

4　接觸問題的處理

為了保證不連續邊界在變形過程中滿足兩接觸邊之間無拉力、無嵌入，數值流形方法在接觸位置設置數值上的“法向接觸彈簧”；同時，為了模擬不連續面的抗剪能力，還在接觸位置處設置了“切向接觸彈簧”。

在本文的程序中，設定所有接觸的初始狀態都是鎖定的，即法向彈簧和切向彈簧同時對接觸兩側的相對位移起著限制作用。施加荷載后，隨著變形的繼續發展，接觸可以在鎖定、滑動和張開三種模式之間轉換(稱為開合迭代)，以模擬不連續面閉合、滑動和張開三種不同的狀態。接觸三種模式間的轉換對平衡方程的影響見文獻[2]。

接觸位置的確定，文獻[6]給出如下原則：

接觸面的邊緣是：1)塊體的邊界；2)材料場中裂縫兩側的一邊。

接觸的頂角可以是：1)物理邊界或裂縫的頂端；2)塊體邊界或裂縫的三角形的交點。

在各種接觸情形中，接觸面的邊緣一定如上所述是塊體的邊界或裂縫的一邊；而接觸的頂角卻不僅限于上述兩種情況。如圖4所示，雙層簡支梁中部受集中荷載作用。如果只在裂縫的頂端，即A，B兩點設置接觸彈簧，雖能確保裂縫頂端處不發生嵌入，卻不能避免在不連續面的其他部位發生嵌入。計算后的變形如圖5所示。

故本文除了在裂縫頂端及裂縫交點處設置初始接觸之外，還在裂縫的整個長度范圍內的每個單元頂點處設置了“初始假接觸”。與裂縫頂端處的初始接觸不同，這些接觸的初始狀態不是閉合而是張開的，即并沒有設置法向和切向彈簧。只有在計算過程中當這些接觸位置發生嵌入時，才添加接觸彈簧，令其從張開變為鎖定。之后，這些接觸也可以隨著變形的發展，在鎖定、滑動和張開三種模式下進行轉換。圖6為上述雙層簡支梁采用“初始假接觸”的辦法計算后的變形圖。可見，這種處理方式可以有效地避免在不連續面任意位置發生嵌入。

5　應用算例

用有解析解的連續體變形問題驗證了本文程序計算結果的可靠性，薄板計算模型見圖7，表1給出一帶圓孔薄板的應力結果比較。

表1　計算結果對比

程序在模擬不連續變形方面也取得了滿意的結果。圖8，圖9為雙層構架受集中力的變形情況。

6　結論及展望

采用本文所述主體思路完成的數值流形程序能夠用于材料變形分析，經驗證，計算出的變形及應力場符合理論解。特別是接

觸彈簧及開合迭代算法的應用，使程序在進行不連續變形分析方面具有獨特的優勢。目前，數值流形程序中材料的本構方程仍然是基于彈性理論的線性方程，下一步的工作將研究改進現有的彈性本構，以期描述巖土材料的非線性永久變形和應變軟化特征。

[1]Shi Genhua. Block system modeling by discontinuous deformation analysis [M].Southampton:Computational Mechanics Publications,1993.

[2]Shi Genhua.Numerical Manifold Method.In:Proc of IFDDA’96.Berkeley,California,USA:Tsi press,1996:52-204.

[3]Chen G，Miki S，Ohnishi Y.Automatic Creation of Mathematical Meshes in manifold method of Material Analysis [A].In:Working Forum on the Manifold Method of Material Analysis[C].California， USA:[s n.]，1995:105-126.

[4]曹文貴，速寶玉.流形元覆蓋系統自動形成方法之研究[J].巖土工程學報，2001，23(20):187-190.

[5]蔡永昌，朱合華，夏才初.流形方法覆蓋系統自動生成算法[J].同濟大學學報，2004，32(5):585-590.

[6]石根華.數值流形方法與非連續變形分析[M].裴覺民，譯.北京：清華大學出版社，1997.

Numerical manifold program development design and application

Jiao Jian

(BeijingUrbanConstructionDesign&DevelopmentGroupCo.,Limited,Beijing100037,China)

In light of cover system subdivision algorithm of numerical manifold program, based on over balancing equation establishment and contact problems processing and other difficulties, the paper puts forward specific solving countermeasures, finishes numerical manifold program compilation, testifies the program reliability through checking, and finally shows its unique advantages in processing discontinuous deformation computation problems.

numerical manifold method, discontinuous deformation computation, program development

1009-6825(2016)08-0256-03

2016-01-05

焦健(1981- )，男，工程師

TB115

A