焦　健

(北京城建設(shè)計(jì)發(fā)展集團(tuán)股份有限公司，北京　100037)

?

數(shù)值流形程序開發(fā)設(shè)計(jì)及其應(yīng)用

焦健

(北京城建設(shè)計(jì)發(fā)展集團(tuán)股份有限公司，北京100037)

針對(duì)數(shù)值流形程序開發(fā)面臨覆蓋系統(tǒng)剖分算法、基于覆蓋的平衡方程的建立、接觸問題處理等幾個(gè)難題，提出了具體的解決對(duì)策，完成了數(shù)值流形程序的編制，并通過算例驗(yàn)證了該程序的可靠性，展現(xiàn)了該程序在處理不連續(xù)變形計(jì)算方面的獨(dú)特優(yōu)勢(shì)。

數(shù)值流形方法，不連續(xù)變形計(jì)算，程序開發(fā)

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

自數(shù)值流形方法提出以來一直頗受關(guān)注與好評(píng)，但石根華博士的數(shù)值流形程序并未在國內(nèi)作大力推廣，而科研工作者獨(dú)自將其編程實(shí)現(xiàn)又面臨諸多困難[3,4]，因此數(shù)值流形方法在國內(nèi)的研究和應(yīng)用還僅限于為數(shù)不多幾個(gè)科研機(jī)構(gòu)。在前人工作的基礎(chǔ)上[5]，完成了數(shù)值流形程序的開發(fā)設(shè)計(jì)，并成功地用于幾個(gè)實(shí)際問題的解決，取得了滿意的結(jié)果。本文重點(diǎn)介紹數(shù)值流形程序編制中幾個(gè)難點(diǎn)問題的解決方法。

1　覆蓋系統(tǒng)自動(dòng)剖分算法

數(shù)值流形方法中采用兩套獨(dú)立的網(wǎng)格，物理網(wǎng)格和數(shù)學(xué)網(wǎng)格。物理網(wǎng)格為材料固有邊界和內(nèi)部裂紋；數(shù)學(xué)網(wǎng)格可以由常規(guī)的網(wǎng)格，如有限元網(wǎng)格、規(guī)則的格子或級(jí)數(shù)收斂域轉(zhuǎn)化而得到。采用三角形或四邊形有限元網(wǎng)格，可以方便地取節(jié)點(diǎn)的形函數(shù)作為覆蓋函數(shù)的權(quán)函數(shù)，是目前主要的數(shù)學(xué)網(wǎng)格形式。二者相比，四邊形有限元網(wǎng)格具有拓?fù)溥\(yùn)算簡單、同階覆蓋函數(shù)條件下精度更高的優(yōu)勢(shì)。此外，數(shù)值流形方法并不要求數(shù)學(xué)網(wǎng)格與求解域完全吻合，而是能夠覆蓋整個(gè)求解域即可。所以，即使在邊界不規(guī)則處也無需調(diào)整網(wǎng)格形狀以適應(yīng)邊界條件。基于上述考慮，本文采用矩形數(shù)學(xué)網(wǎng)格，在此基礎(chǔ)上完成劃分流形單元和覆蓋系統(tǒng)自動(dòng)剖分的工作。

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

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

再對(duì)每一個(gè)數(shù)學(xué)網(wǎng)格節(jié)點(diǎn)進(jìn)行循環(huán)，獲得各個(gè)節(jié)點(diǎn)所對(duì)應(yīng)的物理覆蓋。算法如下：針對(duì)當(dāng)前數(shù)學(xué)網(wǎng)格節(jié)點(diǎn)，搜索所有包含當(dāng)前節(jié)點(diǎn)的流形單元并求其并集，構(gòu)成該節(jié)點(diǎn)的數(shù)學(xué)覆蓋。然后，令物理網(wǎng)格對(duì)該數(shù)學(xué)覆蓋進(jìn)行剖面，剖分處的每一個(gè)單獨(dú)區(qū)域即構(gòu)成一個(gè)物理覆蓋，將其編號(hào)記錄。

獲得物理覆蓋之后，便可以搜索每個(gè)流形單元所屬的物理覆蓋。如前述，每一個(gè)流形單元均記錄了其所在數(shù)學(xué)單元的四個(gè)節(jié)點(diǎn)，這四個(gè)節(jié)點(diǎn)均對(duì)應(yīng)若干個(gè)物理覆蓋，依次對(duì)該流形單元與四個(gè)節(jié)點(diǎn)的物理覆蓋進(jìn)行分析，判斷流形單元完全包含于哪個(gè)物理覆蓋中，共可以找到四個(gè)物理覆蓋，這四個(gè)物理覆蓋即為該流形單元所屬的物理覆蓋。

2　覆蓋函數(shù)和權(quán)函數(shù)

取覆蓋函數(shù)為常函數(shù)，矩形有限元的形函數(shù)作為權(quán)函數(shù)，即：

(1)

(2)

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

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

(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)個(gè)物理覆蓋上的常覆蓋函數(shù){ui,vi}。

3　平衡方程的建立

數(shù)值流形方法依據(jù)最小勢(shì)能原理建立平衡方程?？倓?shì)能包括：應(yīng)變勢(shì)能、初應(yīng)力勢(shì)能、荷載勢(shì)能、慣性力勢(shì)能、彈簧勢(shì)能和接觸摩擦力勢(shì)能。對(duì)它們分項(xiàng)求導(dǎo)使勢(shì)能極小化則可得總體平衡方程:

(7)

其中,總剛度矩陣K由單元?jiǎng)偠染仃嚒⒐潭c(diǎn)矩陣、慣性力矩陣和接觸彈簧矩陣迭加得到，各矩陣的具體形式及計(jì)算某些矩陣元素涉及到的單純形積分問題見文獻(xiàn)[2]。荷載矩陣由F初應(yīng)力矩陣、點(diǎn)荷載矩陣、分布荷載矩陣、體荷載矩陣、速度矩陣以及接觸摩擦力矩陣迭加組成。文獻(xiàn)[2]給出了除分布荷載矩陣外其他矩陣的具體形式；而在實(shí)際問題中，分布荷載也是很常見的荷載形式。本文推導(dǎo)建立了分布荷載矩陣。

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

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

(8)

(9)

繼而，分布荷載在該單元上的勢(shì)能有如下表達(dá)式:

(10)

4　接觸問題的處理

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

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

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

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

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

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

故本文除了在裂縫頂端及裂縫交點(diǎn)處設(shè)置初始接觸之外，還在裂縫的整個(gè)長度范圍內(nèi)的每個(gè)單元頂點(diǎn)處設(shè)置了“初始假接觸”。與裂縫頂端處的初始接觸不同，這些接觸的初始狀態(tài)不是閉合而是張開的，即并沒有設(shè)置法向和切向彈簧。只有在計(jì)算過程中當(dāng)這些接觸位置發(fā)生嵌入時(shí)，才添加接觸彈簧，令其從張開變?yōu)殒i定。之后，這些接觸也可以隨著變形的發(fā)展，在鎖定、滑動(dòng)和張開三種模式下進(jìn)行轉(zhuǎn)換。圖6為上述雙層簡支梁采用“初始假接觸”的辦法計(jì)算后的變形圖?？梢姡@種處理方式可以有效地避免在不連續(xù)面任意位置發(fā)生嵌入。

5　應(yīng)用算例

用有解析解的連續(xù)體變形問題驗(yàn)證了本文程序計(jì)算結(jié)果的可靠性，薄板計(jì)算模型見圖7，表1給出一帶圓孔薄板的應(yīng)力結(jié)果比較。

表1　計(jì)算結(jié)果對(duì)比

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

6　結(jié)論及展望

采用本文所述主體思路完成的數(shù)值流形程序能夠用于材料變形分析，經(jīng)驗(yàn)證，計(jì)算出的變形及應(yīng)力場(chǎng)符合理論解。特別是接

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

[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]曹文貴，速寶玉.流形元覆蓋系統(tǒng)自動(dòng)形成方法之研究[J].巖土工程學(xué)報(bào)，2001，23(20):187-190.

[5]蔡永昌，朱合華，夏才初.流形方法覆蓋系統(tǒng)自動(dòng)生成算法[J].同濟(jì)大學(xué)學(xué)報(bào)，2004，32(5):585-590.

[6]石根華.數(shù)值流形方法與非連續(xù)變形分析[M].裴覺民，譯.北京：清華大學(xué)出版社，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