巖體工程數(shù)值流形方法的固定邊界約束處理方法

1、第24卷 第19期 巖石力學(xué)與工程學(xué)報(bào) Vol.24 No.19 2005年10月 Chinese Journal of Rock Mechanics and Engineering Oct.，2005巖體工程數(shù)值流形方法的固定邊界約束處理方法朱愛軍1，鄧安福1，曾祥勇1，唐樹名2(1. 重慶大學(xué) 土木工程學(xué)院，重慶 400045；2. 重慶交通科研設(shè)計(jì)院，重慶 400067)摘要：在巖體工程數(shù)值流形方法分析中，往往存在較多的固定約束邊界。目前，數(shù)值流形方法一般采用罰函數(shù)方法處理固定邊界問題，但罰彈簧的布置與大小對(duì)數(shù)值模擬的效果有一定的影響，且處理方式也比較復(fù)雜?；诹餍螁卧衔灰坪瘮?shù)的組成提

2、出固定約束處理的新方法，將固定邊界的約束處理轉(zhuǎn)化為對(duì)廣義節(jié)點(diǎn)的約束處理，改變廣義節(jié)點(diǎn)上的覆蓋函數(shù)使固定邊界的約束條件得到嚴(yán)格滿足，并推導(dǎo)相應(yīng)的流形單元?jiǎng)偠染仃?。該方法在物理意義上嚴(yán)格滿足固定邊界的約束條件，同時(shí)簡(jiǎn)化了處理工作，有利于數(shù)值流形方法程序的實(shí)現(xiàn)和工程應(yīng)用。關(guān)鍵詞：數(shù)值分析；數(shù)值流形方法；廣義節(jié)點(diǎn)；覆蓋函數(shù)；固定邊界；約束處理中圖分類號(hào)：O 241 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：10006915(2005)19358206CONSTRAINTS TREATMENT OF FIXED BOUNDARY IN NUMERICALMANIFOLD METHODZHU Ai-jun1，DENG A

3、n-fu1，ZENG Xiang-yong1，TANG Shu-ming2(1. College of Civil Engineering，Chongqing University，Chongqing 400045，China；2. Chongqing Communications Research and Design Institute，Chongqing 400067，China)Abstract：In numerical manifold method，fixed boundary is treated by penalty function method. The arrangeme

4、nt and value of penalty spring will affect the effect of numerical modeling；and work is complex. Based on the composition of displacement function on manifold element，a new method of constraints treatment of fixed boundary is proposed；and it changes the constraints treatment of fixed boundary into t

5、he constraints treatment of generalized node；and the constraints condition will be met strictly by using special cover function on physical cover；and then corresponding stiffness matrix is deduced. New method had been proved theoretically that it is correct；and it not only meets strictly the constra

6、ints condition in physical meaning，but it has simplified work of treatment. It is helpful for programming and the application of numerical manifold method to engineerings.Key words：numerical analysis；numerical manifold method；generalized node；cover function；fixed boundary；constraints treatment1 引 言數(shù)

7、值流形方法是文1，2在研究DDA和有限元的數(shù)學(xué)基礎(chǔ)后提出的一種新的數(shù)值方法，與傳統(tǒng)數(shù)值方法相比，該方法具有連續(xù)與非連續(xù)問題統(tǒng)一分收稿日期：20040523；修回日期：20040702基金項(xiàng)目：西部交通建設(shè)科技項(xiàng)目(200231800026) 析及高精度求解的特點(diǎn)。國(guó)內(nèi)外對(duì)流形方法的研究中，有很多將其應(yīng)用于巖體工程不連續(xù)變形問題的模擬分析及材料裂隙擴(kuò)展的跟蹤模擬3類研究中以梁等小規(guī)模的算例居多。 在數(shù)值流形方法分析中，無(wú)論是巖體工程這類9。同時(shí)，12對(duì)流形方法的高精度特性也有較多的研究10，該作者簡(jiǎn)介：朱愛軍(1975)，男，2005年于重慶大學(xué)土木工程學(xué)院獲博士學(xué)位，主要從事巖土工程方面的研究

8、工作。E-mail：zhuai juncqu。第24卷 第19期 朱愛軍等. 巖體工程數(shù)值流形方法的固定邊界約束處理方法 3583 規(guī)模大、結(jié)構(gòu)復(fù)雜的問題還是梁等小規(guī)模問題，都涉及到固定約束邊界的處理，固定約束邊界的處理方法將對(duì)計(jì)算的效率與精度產(chǎn)生影響。數(shù)值流形方法中，固定約束邊界一直沿用罰函數(shù)方法進(jìn)行處理，即在固定邊界上設(shè)置剛性彈簧，在平衡方程中加入相應(yīng)的固定點(diǎn)矩陣。從流形方法發(fā)展至今的研究證明，罰函數(shù)方法是有效和可行的。但從理論上講，罰函數(shù)方法不是一種嚴(yán)格的固定約束處理方法，罰彈簧的布置與大小對(duì)數(shù)值模擬的效果有一定的影響，同時(shí)也存在彈簧取值較小而引起較大誤差或彈簧取值過大導(dǎo)致解奇異的可能性

9、。甚至當(dāng)分析的重點(diǎn)為固定邊界上的應(yīng)力時(shí)，罰函數(shù)方法可能得不到準(zhǔn)確的結(jié)果。目前對(duì)固定邊界約束處理方面所開展的研究較少，文13對(duì)原有罰函數(shù)的方法進(jìn)行了一定的改進(jìn)。本文則通過改變廣義節(jié)點(diǎn)上的覆蓋函數(shù)，并采用相應(yīng)的剛度矩陣計(jì)算方法來對(duì)固定邊界的約束進(jìn)行處理，該方法從理論上嚴(yán)格符合固定約束的物理意義，同時(shí)簡(jiǎn)化了處理工作，并經(jīng)筆者大量算例證實(shí)是有效和準(zhǔn)確的。2 物理覆蓋的覆蓋函數(shù)和流形單元的位移函數(shù)2.1 物理覆蓋的覆蓋函數(shù)物理覆蓋上的覆蓋函數(shù)通常取為坐標(biāo)(x，y)的級(jí)數(shù)形式，目前多采用(x，y)的0，1和2階完全多項(xiàng)式。對(duì)平面問題，如果物理覆蓋為Ui，則其覆蓋函數(shù)為ui(x，y)mfij(x，y)0di

10、，2j1vi(x，y)=j=10fij(x，y)d (1) i，2j式中：fij(x，y)為完全多項(xiàng)式的基本級(jí)數(shù)；di，2j1，di，2j均為物理覆蓋的覆蓋自由度，對(duì)應(yīng)完全多項(xiàng)式的各系數(shù)。流形單元的覆蓋函數(shù)可取0階到高階的完全多項(xiàng)式，對(duì)應(yīng)0，1和2階情況，式(1)中m分別為1，3和6。2.2 流形單元的位移函數(shù)各物理覆蓋的覆蓋函數(shù)在其交集流形單元上通過權(quán)函數(shù)對(duì)其進(jìn)行加權(quán)平均構(gòu)成流形單元的位移函數(shù)。流形方法里的權(quán)函數(shù)與有限元法里的插值形函數(shù)的物理意義和定義原則一致。如果流形單元e由q個(gè)物理覆蓋Ue1，Ue2，L，Ueq的交集形成，令物理覆蓋Uei上的權(quán)函數(shù)為wei(x，y)。則流行單元e的位移函

11、數(shù)為ue(x，y)qv=e(x，y)wuei(x，y)ei(x，y)=i=1vei(x，y)qmwei(x，y)feij(x，y)00w)f i=1j=1ei(x，yeij(x，y)dei，2j1 d(2) ei，2j令物理覆蓋Uei的變形矩陣為Tei=Tei1Tei2LTeim (3)其中，Tte1，2j1(x，y)te1，2j(x，y)eij=te2，2j1(x，y)te2,2j(x，y)= wei(x，y)feij(x，y)00wei(x，y)feij(x，y) (4) 令Dei=dei，1dTei，2Ldei，2m1dei，2m (5)式中：Dei為物理覆蓋自由度，為2m階向量。由此，

12、式(2)可表示為ue(x，y)qv=e(x，y)TeiDei (i = 1，2，L，q) (6)i=13 流形單元的剛度矩陣與有限元法一樣，應(yīng)力與應(yīng)變之間的關(guān)系可表示為=E (7)式中：為應(yīng)力向量，為應(yīng)變向量，E為勁度矩陣。對(duì)平面問題，由流形單元的總體位移函數(shù)和應(yīng)變的關(guān)系可以得到應(yīng)變表達(dá)式為De1xDe2y=Be1Be2LBeq=BMeDe (8)xyDeq式中：Bei為物理覆蓋Uei的應(yīng)變矩陣，且Ke11Ke12Ke13Ke21Ke22Ke23 (11) Ke31Ke32Ke33其中，Kers= ABerTEBesdxdy(r = 1，2，3；s = 1，2，3) (12)如果在組成流形單元

13、e的各物理覆蓋上使用相同階次的覆蓋函數(shù)，則式(11)中各子矩陣的階數(shù)為2m×2m階。4 固定邊界約束問題的處理本文將固定邊界的約束處理轉(zhuǎn)化為對(duì)廣義節(jié)點(diǎn)的約束處理，通過對(duì)廣義節(jié)點(diǎn)的處理來實(shí)現(xiàn)對(duì)固定邊界約束條件的滿足。該方法具有一般性，即對(duì)各階次的覆蓋函數(shù)及各種形狀的數(shù)學(xué)網(wǎng)格都是適用的。本文采用三角形數(shù)學(xué)網(wǎng)格的情況進(jìn)行論述。圖1為開挖隧道計(jì)算模型的約束邊界。設(shè)模型沿豎線AD對(duì)稱，因此分析時(shí)可以只取模型的一半，則計(jì)算模型的約束條件為：AB邊受雙向固定約束；BC，DE和FA三邊受x向固定約束，y向自由。對(duì)其進(jìn)行流形方法分析時(shí)，如果采用三角形的數(shù)學(xué)網(wǎng)格，則對(duì)該模型構(gòu)造的流形元覆蓋系統(tǒng)如圖1所示

14、。圖1中，網(wǎng)格節(jié)點(diǎn)上的數(shù)字為廣義節(jié)點(diǎn)的編號(hào)，即物理覆蓋的編號(hào)，網(wǎng)格中央的數(shù)字為流形單元的4.1 只包含單向位移的覆蓋函數(shù)如果物理覆蓋為Ui，采用(x，y)的級(jí)數(shù)構(gòu)造其覆蓋函數(shù)，則只包含x向位移的覆蓋函數(shù)為 ui(x，y)ui(x，y)mfijvi(x，y)=0=(x，y)dij (13)j=10物理覆蓋Ui上只包含y向位移的覆蓋函數(shù)為 ui(x，y)0m0vi(x，y)=vi(x，y)=j=1fij(x，y)dij (14)4.2 含雙向固定約束條件的流形單元?jiǎng)偠染仃噲D2中流形單元e與圖1中一條邊受雙向固定約束的單元相對(duì)應(yīng)，例如單元15，17等。為便于敘述，設(shè)圖2中單元的廣義節(jié)點(diǎn)號(hào)為1，2，3

15、，對(duì)應(yīng)的物理覆蓋為Ue1，Ue2和Ue3，其中廣義節(jié)點(diǎn)1，2位于固定邊上，受雙向固定約束，節(jié)點(diǎn)3為自由節(jié)點(diǎn)。3 e 1圖2 一條邊受雙向固定約束的流形單元 Fig.2 Manifold element that a side is in two-directionconstraint第24卷 第19期 朱愛軍等. 巖體工程數(shù)值流形方法的固定邊界約束處理方法 3585 流行單元e的位移函數(shù)為ue(x，y)3uei(x，y)ve(x，y)=wei(x，y) (15) i=1vei(x，y)令ue1(x，y)ue2(x，y)v=0e1(x，y)=ve2(x，y) 0則式(15)可寫為ue(x，y)v

16、=wue3(x，y)e (16) e(x，y)3(x，y)ve3(x，y)根據(jù)權(quán)函數(shù)定義，在流行單元e的邊12上有 we3(x，y)=0 (x，y)12因此，由式(16)表示的單元位移函數(shù)在單元邊12上嚴(yán)格滿足位移為0的條件，即ue(x，y)v=0 (x，y)1e(x，y)2 0此時(shí)，與物理覆蓋Ue3所對(duì)應(yīng)的應(yīng)變矩陣為Be3，所以流行單元e的剛度矩陣為KTe= ABe3EBe3dxdy (17)與式(11)表達(dá)的單元?jiǎng)偠染仃囅啾?，?17)只由一個(gè)子矩陣構(gòu)成，其大小為式 (11)的1/9。4.3含單向固定約束條件的流形單元?jiǎng)偠染仃噲D3中的流形單元與圖1中一條邊受x向固定約束的單元相對(duì)應(yīng)，如單元

17、40，209等。圖3 一條邊受x向固定約束的流形單元Fig.3 Manifold element that a side is in x-direction constraint圖3中，廣義節(jié)點(diǎn)1，2位于x向固定約束邊上，受x向固定約束。流行單元e的位移函數(shù)可用式(15)表達(dá)。在物理覆蓋Ue1，Ue2上采用只包含y方向位移的覆蓋函數(shù)，即式(13)，則單元位移函數(shù)為ue(x，y)0ve(x，y)=we1(x，y)ve1(x，y)+ w0ue3(x，y)e2vwx，e2(x，y)+e3(y)v (18) e3(x，y)在單元的邊12上，有we3(x，y)=0 (x，y)12所以，由式(18)表達(dá)的

18、單元位移函數(shù)在邊12上嚴(yán)格滿足x方向位移為0的條件，即ue(x，y)=0(x，y)12對(duì)物理覆蓋Ue1，其加權(quán)覆蓋函數(shù)可表達(dá)為wue1(x，y)e1v=Te1De1 e1(x，y)(19)根據(jù)式(14)，式(19)有T00L0e1=t= e2，1(x，y)te2，2(x，y)Lte2，m(x，y)T0we1(x，y)fe1，1(x，y)0w(x，y)f(x，y)e1e1，2 M (20) 0we1(x，y)fe1，m(x，y)De1=de1，1de1，2Lde1，mT (21)式(20)為只包含y向位移的物理覆蓋變形矩陣。式(21)為相應(yīng)的覆蓋自由度向量，其階數(shù)為m，為式(5)所表達(dá)的覆蓋自由

19、度階數(shù)(2m階)的一半。由此，物理覆蓋Ue1的應(yīng)變矩陣為00L0Bte2，1te2，2te2，mLe1= yyy(22)te2，1tte2，2e2，mxxLx同理，Ue2的應(yīng)變矩陣Be2可通過以上過程得到。式(22)為只包含y向位移的物理覆蓋應(yīng)變矩陣，對(duì)于只包含x向位移的物理覆蓋應(yīng)變矩陣可用同 樣的步驟得到，其與雙向自由的物理覆蓋應(yīng)變矩陣(式(9)是不同的。流行單元e中，廣義節(jié)點(diǎn)3為雙向自由節(jié)點(diǎn)，其對(duì)應(yīng)的應(yīng)變矩陣Be3同式(9)。得到Ue1，Ue2和Ue3的應(yīng)變矩陣Be1，Be2和Be3后即可按式(10)和(11)的步驟得到單元的剛度矩陣Ke。如果在單元的3個(gè)廣義節(jié)點(diǎn)上使用相同階次的覆蓋函數(shù)，

20、則Ke的9個(gè)子矩陣中，Ke33為2m×2m階矩陣；Ke13，Ke23為m×2m階矩陣；Ke31，Ke32為2m×m階矩陣，其他為m×m階矩 3586 巖石力學(xué)與工程學(xué)報(bào) 2005年陣。以上討論了含固定約束條件單元的剛度矩陣的計(jì)算，對(duì)于相應(yīng)的荷載矩陣或其他矩陣(如慣性力矩陣等)的計(jì)算方法與計(jì)算剛度矩陣的原理是一樣的。以上以圖2，3所示的2種單元情況為例論述了對(duì)固定邊界的處理，說明了對(duì)固定邊界的約束處理最終可以轉(zhuǎn)化為對(duì)廣義節(jié)點(diǎn)的約束處理。就圖1的模型而言，還有如圖4所示的幾種含約束條件的單元類型。圖4與圖1中單元的對(duì)應(yīng)關(guān)系為：?jiǎn)卧?圖4(a)單元1；單元(圖

21、4(b)單元4等；單元(圖4(c)單元22，39等；單元(圖4(d)單元2；單元(圖4(e)單元20。(a) (b)(c) (d)(e)圖4 受不同約束條件的流形單元Fig.4 Manifold elements in different constraint conditions由此可見，單元受約束的類型很多，但再?gòu)?fù)雜的單元約束條件轉(zhuǎn)化為廣義節(jié)點(diǎn)約束條件，只有3種情況：(1) 廣義節(jié)點(diǎn)受雙向固定約束；(2) 廣義節(jié)點(diǎn)受x向的固定約束；(3) 廣義節(jié)點(diǎn)受y向的固定約束。由此，將固定邊界約束條件轉(zhuǎn)化為廣義節(jié)點(diǎn)約束條件后，約束處理則是簡(jiǎn)單的。對(duì)約束處理一般化的敘述為：當(dāng)廣義節(jié)點(diǎn)受雙向固定約束時(shí)，取

22、消該廣義節(jié)點(diǎn)覆蓋函數(shù)在單元位移函數(shù)中的組成部分；當(dāng)廣義節(jié)點(diǎn)只受x或y向固定約束時(shí)，在該節(jié)點(diǎn)上使用只包含自由方向位移的覆蓋函數(shù)。由此得到流形單元的位移函數(shù)來進(jìn)行其他公式(剛度矩陣、荷矩陣、慣性力矩陣等)的推導(dǎo)，則可準(zhǔn)確地對(duì)固定約束問題進(jìn)行處理。如果所取的分析對(duì)象的約束條件是正確的，則采用上述方法計(jì)算得到的系統(tǒng)總剛度矩陣為對(duì)稱正定陣，不需要進(jìn)行其他處理。需要說明的是，對(duì)于受固定約束的邊界，在生成覆蓋系統(tǒng)的前處理中，必須設(shè)置數(shù)學(xué)網(wǎng)格線通過，以使在固定邊界上生成廣義節(jié)點(diǎn)，通過對(duì)廣義節(jié)點(diǎn)的處理實(shí)現(xiàn)邊界固定約束的處理。5 算 例為了便于觀察，本文采用一簡(jiǎn)單的計(jì)算模型作為算例，如圖5所示。圖5 計(jì)算模型及流

23、形單元?jiǎng)澐?Fig.5 Model used in calculation梁的參數(shù)為：彈性模量E = 1.5 MPa，泊松比v =0.24；荷載f1= 150 N，f2=1 000 N。梁的約束條件為：點(diǎn)A受雙向固定約束；AB邊受x向固定約束。用數(shù)值流形方法對(duì)該梁進(jìn)行計(jì)算，將梁劃分為4個(gè)流形單元。在物理覆蓋上采用2階覆蓋函數(shù)。取消節(jié)點(diǎn)A對(duì)應(yīng)的物理覆蓋上的覆蓋函數(shù)；節(jié)點(diǎn)B對(duì)應(yīng)的物理覆蓋上，采用只有y向位移的覆蓋函數(shù)，由此可嚴(yán)格滿足上述約束條件。表1為梁在荷載f1作用下，BC邊撓度的解析解與流形方法計(jì)算結(jié)果的對(duì)比。圖6中繪出了梁在荷載f2作用下，流形方法計(jì)算得到的梁AB截面的x向應(yīng)力分布圖。表1

24、梁BC邊撓度Table 1 Deflections of the beam BC mm距B點(diǎn)的距離/m方法2 4 6 8 10解析解 5.60 20.80 43.20 70.40 100.00 流形方法5.51 20.47 42.77 69.66 99.20表1中，流形方法僅對(duì)梁劃分了4個(gè)單元，得到的撓度計(jì)算結(jié)果與解析解非常接近，最大誤差在2%以下。圖6中，AB截面為受x向固定約束的邊界，其應(yīng)力的大小為流形方法計(jì)算得到的直接結(jié)第24卷 第19期 朱愛軍等. 巖體工程數(shù)值流形方法的固定邊界約束處理方法 3587 m/高梁水平應(yīng)力/kN圖6 梁AB截面的x向應(yīng)力分布圖Fig.6 Stress di

25、stribution of AB cross-section in x-direction果，沒有經(jīng)過任何處理及位數(shù)的舍入，其數(shù)值與解析解完全一樣。表1和圖6顯示，高階流形方法有較高的計(jì)算精度，同時(shí)說明了本文提出的對(duì)固定邊界約束處理的方法是有效的、精確的。6 結(jié) 論(1) 本文方法在理論上嚴(yán)格符合固定約束的物理意義，相對(duì)于罰函數(shù)方法，它是一種更精確的處理方法，且不存在彈簧取值較小而引起較大誤差及彈簧取值過大導(dǎo)致解奇異的可能性，而且不用引進(jìn)固定點(diǎn)矩陣，簡(jiǎn)化了處理工作。(2) 該方法由于不考慮固定邊界或固定點(diǎn)與固定方向?qū)?yīng)的自由度，減小了平衡方程的規(guī)模，可以提高計(jì)算效率，對(duì)如邊坡工程、路堤工程等含

26、有較多固定邊界的巖土工程問題的求解具有一定意義。參考文獻(xiàn)(References)：1石根華. 數(shù)值流形方法與非連續(xù)變形分析M. 裴覺民譯. 北京：清華大學(xué)出版社，1997. 125253.(Shi Genhua. Numerical Manifold Method and Discontinuous Deformation AnalysisM. Translated by Pei Juemin. Beijing：Tsinghua University Press，1997. 125253.(in Chinese) 2Shi G H. Modeling rock joints and block

27、s by manifold methodA.In：Proceedings of 33rd U S Symposium on Rock MechanicsC. Santa Fe. New Mexico：s. n.，1992. 639648. 3Techin K. Artificial joint-based manifold methodA. In：Working Forum on the Manifold Method of Material Analysis(Vol.1)C.s. l.：s. n.，1995. 2137.4 Tadeshi S，Morikawa S J，Ishii D K.

28、Elastoplastic analysis of jointedrock models by manifold methodA. In：Proceedings of the Second International Conferece on Analysis of Discontinuous DeformationC. New York：JSCE，1997. 309316.5 Tsay R J，Chiou Y J，Chuang W L. Crack growth prediction bymanifold methodJ. Journal of Engineering Mechanics，A

29、SCE，1999，125(8)：884890. 6王芝銀，王思敬，楊志法. 巖石大變形分析的流形方法J. 巖石力學(xué)與工程學(xué)報(bào)，1997，16(5)：399404.(Wang Zhiyin，Wang Sijing，Yang Zhifa. Manifold method in analysis of large deformation for rockJ. Chinese Journal of Rock Mechanics and Engineering，1997，16(5)：399404.(in Chinese) 7張國(guó)新，王光綸，裴覺民. 基于流形方法的結(jié)構(gòu)體破壞分析J. 巖石力學(xué)與工程學(xué)報(bào)，

30、2001，20(3)：281287.(Zhang Guoxin，Wang Guanglun，Pei Juemin. Structure failure analysis based on manifold methodJ. Chinese Journal of Rock Mechanics and Engineering，2001，20(3)：281287.(in Chinese) 8莫海鴻，陳尤雯. 流形法在巖石力學(xué)研究中的應(yīng)用J. 華南理工大學(xué)學(xué)報(bào)(自然科學(xué)版)，1998，26(9)：4853.(Mo Haihong，Chen Youwen. Application of the manifold method to rock mechanicsJ.