六面體覆蓋的高階數(shù)值流形方法的探討

劉登學(xué)，張友良，譚 飛，張禮仁

（中國(guó)科學(xué)院武漢巖土力學(xué)研究所巖土力學(xué)與工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430071）

1 引言

數(shù)值流形方法是石根華[1－2]在不連續(xù)分析方法（DDA）基礎(chǔ)上提出的，是利用“流形”數(shù)學(xué)概念的有限覆蓋技術(shù)建立起來(lái)的統(tǒng)一解決連續(xù)與非連續(xù)變形問(wèn)題的數(shù)值方法，連續(xù)體的有限元方法（FEM）和塊體系統(tǒng)的非連續(xù)變形分析方法（DDA），只是這一數(shù)值方法的特例。這種方法采用2個(gè)網(wǎng)格體系（物理網(wǎng)格和數(shù)學(xué)網(wǎng)格）分析問(wèn)題，在分析域內(nèi)的物理覆蓋上建立一般覆蓋位移函數(shù)加權(quán)求和形成總體位移函數(shù)。當(dāng)流形方法的覆蓋函數(shù)推廣到高階情況時(shí)，其求解精度將會(huì)有很大的提高。

目前二維流形方法的理論已經(jīng)比較完善，而且國(guó)內(nèi)許多學(xué)者也將該理論應(yīng)用到了隧道的開(kāi)挖模擬中[3－5]，并取得了很好的效果，如王水林等[16]、鄭宏等[7]利用數(shù)值流形方法進(jìn)行了二維裂紋擴(kuò)展問(wèn)題的研究。然而，實(shí)際的隧道和地下工程問(wèn)題都是三維問(wèn)題，任何計(jì)算方法最終都應(yīng)該發(fā)展成為三維并服務(wù)于工程實(shí)際，近些年來(lái)，國(guó)內(nèi)外部分學(xué)者圍繞三維數(shù)值流形方法理論及應(yīng)用做了大量的工作。姜清輝[8－11]等基于四面體覆蓋網(wǎng)格建立了0 階及高階的三維數(shù)值流形分析格式，還成功模擬了三維無(wú)壓滲流問(wèn)題，也對(duì)三維數(shù)值流形方法的點(diǎn)面接觸模型進(jìn)行了研究。鄭榕明等[12]對(duì)基于六面體覆蓋的0 階三維數(shù)值流形方法進(jìn)行了探討。何磊等[13]將三維數(shù)值流形方法應(yīng)用到裂隙巖質(zhì)邊坡的穩(wěn)定性分析中。郭旭陽(yáng)等[14]、林毅峰等[15]研究了高階數(shù)值流形方法中存在的線(xiàn)性相關(guān)問(wèn)題。

本文在六面體有限元覆蓋的三維流形方法基礎(chǔ)上，分別采用完全一階近似位移覆蓋函數(shù)和二階近似位移函數(shù)，建立相應(yīng)的分析格式，分析了高階數(shù)值流形法中線(xiàn)性相關(guān)問(wèn)題，編制了相應(yīng)的程序，并對(duì)該方法在隧道與地下工程應(yīng)用做了簡(jiǎn)單的分析和展望。

2 基本理論

2.1 權(quán)函數(shù)與位移函數(shù)

采用8 節(jié)點(diǎn)六面體有限元網(wǎng)格作為數(shù)學(xué)覆蓋網(wǎng)格，取六面體有限單元的形函數(shù)Ni(i=1～8)作為權(quán)函數(shù)，C8 型等參單元的形函數(shù)可以表示為

式中： ξi、ηi、 ζi為六面體第i(i=1:8)個(gè)頂點(diǎn)在局部坐標(biāo)系（ξ，η，ζ）中對(duì)應(yīng)的坐標(biāo)值。

采用兩種不同的局部近似函數(shù)：（1）{1，x，y，z}；（2）{1，x2，y2，z2}，完全一階近似位移函數(shù)帶來(lái)的線(xiàn)性相關(guān)問(wèn)題將在后面部分討論。兩種不同近似函數(shù)對(duì)應(yīng)的物理覆蓋i 的位移函數(shù)可表示為

本文采用8 節(jié)點(diǎn)六面體有限元網(wǎng)格作為數(shù)學(xué)覆蓋，一個(gè)流形單元由8個(gè)物理覆蓋（節(jié)點(diǎn)）組成，流形單元的位移函數(shù)有8個(gè)物理覆蓋的位移函數(shù)加權(quán)平均，可表示為

2.2 三維流形單元的剛度矩陣

對(duì)于流形方法，剛度矩陣的積分區(qū)域?yàn)榱餍螁卧?，與有限元一樣單元應(yīng)變可根據(jù)幾何方程求出：

式中：

類(lèi)似的可以寫(xiě)出局部近似函數(shù)為{1，x2，y2，z2}的應(yīng)變矩陣。

由單元e 彈性應(yīng)力產(chǎn)生的應(yīng)變能為

式中：J為雅克比矩陣；{E}為彈性矩陣。則單元?jiǎng)偠染仃嚍?/p>

2.3 荷載及固定約束子矩

2.3.1 點(diǎn)荷載矩陣

點(diǎn)荷載(Fx，F(xiàn)y，F(xiàn)z)T作用在單元e 中得點(diǎn)(x0，y0，z0)上，則對(duì)應(yīng)的點(diǎn)荷載矩陣為

式中：(ξ0,η0,ζ0)為(x0,y0,z0)在局部坐標(biāo)系中對(duì)應(yīng)的坐標(biāo)值。

2.3.2 體積荷載矩陣

設(shè)(fx,fy,fz)T是不變的體積力，它作用在單元e 的材料體積內(nèi)，則體積力引起的荷載矩陣為

2.3.3 固定約束矩陣

假定固定點(diǎn)是在單元e 的(x0,y0,z0)上，沿xyz方向分別設(shè)置彈簧，彈簧的剛度為p，則對(duì)應(yīng)的位移約束矩陣為

式中：(ξ0,η0,ζ0)為(x0,y0,z0)在局部坐標(biāo)系中對(duì)應(yīng)的坐標(biāo)值。

2.4 等參坐標(biāo)與整體坐標(biāo)的相互映射

因本文采用C8 型等參單元中的形函數(shù)作為三維流形單元的權(quán)函數(shù)，勢(shì)必要涉及坐標(biāo)的變換。等參坐標(biāo)到整體坐標(biāo)的一一對(duì)應(yīng)關(guān)系為

式中： xi、 yi、 zi為六面體結(jié)點(diǎn)在總體（笛卡爾）坐標(biāo)內(nèi)的坐標(biāo)值。

Ni對(duì)于x、y、z 的偏導(dǎo)數(shù)可利用自然坐標(biāo)顯示，表示為

本文數(shù)學(xué)覆蓋網(wǎng)格采用的是規(guī)則的六面體，因而整體坐標(biāo)中任意一點(diǎn)(x，y，z)在等參坐標(biāo)中對(duì)應(yīng)的(ξ,η,ζ)可以通過(guò)式（17）表示：

式中：x0、y0、z0為流形單元所對(duì)應(yīng)的8個(gè)物理覆蓋（節(jié)點(diǎn)）形成的六面體的中心坐標(biāo)；a、b、c為數(shù)學(xué)覆蓋網(wǎng)格的三個(gè)方向的尺寸大小。

3 高階數(shù)值流形法的線(xiàn)性相關(guān)問(wèn)題

數(shù)值流形法中，當(dāng)局部函數(shù)采用不低于一階的多項(xiàng)式時(shí)會(huì)導(dǎo)致總體剛度矩陣K 的過(guò)分虧秩，即使在引入位移邊界消除了模型的剛體位移后，K 仍然是虧秩的，即所謂的線(xiàn)性相關(guān)問(wèn)題。線(xiàn)性相關(guān)會(huì)導(dǎo)致平衡方程的解不惟一，但在物體施加完整的位移約束后，每個(gè)解所對(duì)應(yīng)的位移是惟一的[16]。此外，林毅峰、朱合華等還發(fā)現(xiàn)，采用完全一次多項(xiàng)式局部近似函數(shù)雖然線(xiàn)性相關(guān)，但求解仍然收斂，且精度高于線(xiàn)性無(wú)關(guān)的單元；近似函數(shù)為{1，x2，y2，z2}的流形方法形成的整體剛度矩陣在施加位移邊界條件后滿(mǎn)秩[15]。

針對(duì)完全一階近似位移函數(shù)導(dǎo)致的線(xiàn)性相關(guān)問(wèn)題，本文在求解平衡方程時(shí)采用郭朝旭等提出的改進(jìn)的LDLT 算法[14]，穩(wěn)定的求出方程的一組特解。

4 算 例

本文采用MATLAB 軟件基于所建立的分析格式，編制了相應(yīng)的計(jì)算程序，并對(duì)2個(gè)算例進(jìn)行了驗(yàn)證。

算例1.計(jì)算模型左端固定，右端受到大小為1 kPa 的均布荷載，模型尺寸為4 m×2 m×2 m，物體材料常數(shù)：E=1 GPa，v=0.25，不考慮物體的重力。算例1 的計(jì)算模型如圖1 所示。采用如圖2 所示的數(shù)學(xué)覆蓋網(wǎng)格，因此模型為連續(xù)體流形單元的物理覆蓋數(shù)等于有限元的節(jié)點(diǎn)數(shù)，共生成16個(gè)流形單元，45個(gè)物理覆蓋（節(jié)點(diǎn)）。圖3為物體沿x 方向的位移。

圖1 物理計(jì)算模型Fig.1 Physical computational model

圖2 數(shù)學(xué)覆蓋網(wǎng)格Fig.2 Mathematical covers for numerical problem

圖3 物體x 方向的位移Fig.3 Displacement in x direction

從圖3 可以看出，采用0 階近似位移函數(shù)時(shí)，計(jì)算結(jié)果與理論解有一定的誤差；采用完全一階近似函數(shù)時(shí)，剛度矩陣的虧秩數(shù)量為107，利用修正的LDLT 算法穩(wěn)定的得出了一組特解，計(jì)算結(jié)果能很好的逼近理論解；采用{1，x2，y2，z2}二階近似函數(shù)時(shí)，求解沒(méi)有出現(xiàn)線(xiàn)性相關(guān)問(wèn)題，結(jié)果也能很好地逼近理論解，但計(jì)算量及計(jì)算時(shí)間顯著增加，該結(jié)果也驗(yàn)證了本程序的正確性。

算例2.大小為7 m×6 m×6 m 矩形物體，中間有一個(gè)1 m×1 m×6 m 的矩形孔洞，物體下端，左右端固定，計(jì)算模型見(jiàn)圖4。材料常數(shù)：E=2 GPa，v=0.4，ρ=2 000 kg/m3。采用如圖5 所示的數(shù)學(xué)覆蓋網(wǎng)格，共生成252個(gè)流形單元，392個(gè)物理覆蓋（節(jié)點(diǎn)）。

圖4 物理計(jì)算模型Fig.4 Physical computational model

圖5 數(shù)學(xué)覆蓋網(wǎng)格Fig.5 Mathematical covers for numerical problem

圖6為在自重應(yīng)力作用下矩形孔洞兩側(cè)（x=3 m或4 m）的豎向（z 方向）的位移曲線(xiàn)。當(dāng)局部近似函數(shù)采用完全一階時(shí)，整體剛度矩陣的秩為4 149，虧秩數(shù)量為555；當(dāng)采用{1，x2，y2，z2}的二階局部近似函數(shù)時(shí)，整體剛度矩陣的秩為4 704，滿(mǎn)秩。從圖可以看出，采用高階的三維數(shù)值流形方法對(duì)該問(wèn)題的模擬計(jì)算結(jié)果與采用ABAQUS 等有限元的模擬計(jì)算結(jié)果基本一致。算例2 通過(guò)對(duì)簡(jiǎn)單地下洞室模型的計(jì)算，表明三維高階數(shù)值流形方法在計(jì)算該方面的問(wèn)題時(shí)是可行的。

圖6 孔洞兩側(cè)的豎向位移Fig.6 The vertical displacement along the side of the hole

5 結(jié)語(yǔ)

本文基于六面體的數(shù)學(xué)覆蓋網(wǎng)格，建立了高階的三維數(shù)值流形方法分析格式，利用MATLAB 編制了相應(yīng)的計(jì)算程序。通過(guò)對(duì)算例1 的計(jì)算，驗(yàn)證了所編制的計(jì)算程序的正確性，其計(jì)算結(jié)果還表明高階的數(shù)值流形方法相對(duì)于0 階的流形方法有更高的精度。算例2 通過(guò)對(duì)簡(jiǎn)單地下洞室模型的計(jì)算，驗(yàn)證了三維高階數(shù)值流形方法在計(jì)算地下洞室變形方面的有效性。

在隧道及地下工程分析中必需控制求解精度，在數(shù)值流形法中則通過(guò)控制數(shù)學(xué)覆蓋的網(wǎng)格的稀疏和覆蓋位移的階數(shù)來(lái)達(dá)到精度的要求。隧道和地下工程不可避免的會(huì)存在大量的節(jié)理裂隙，數(shù)值流形法中由于采用2 套獨(dú)立覆蓋系統(tǒng)，在不連續(xù)面處采用非連續(xù)的覆蓋函數(shù)，可實(shí)現(xiàn)對(duì)不連續(xù)、大變形問(wèn)題的分析。數(shù)值流形方法在隧道及地下工程的應(yīng)用前景是可預(yù)見(jiàn)的。

[1]SHI G H.Manifold method[C]//Proc.of the 1st Int.Forum on DDA Simulation of Discontinuous Media.Bekerley:[s.n.],1996:52－204.

[2]石根華.數(shù)值流形方法與非連續(xù)變形分析[M].北京:清華大學(xué)出版社,1997.

[3]焦健,喬春生,徐干成.開(kāi)挖模擬在數(shù)值流形方法中的實(shí)現(xiàn)[J].巖土力學(xué),2010,31(9):2951－2957.JIAO Jian,QIAO Chun-sheng,XU Gan-cheng.Simulation of excavation in numerical manifold method[J].Rock and Soil Mechanics,2010,31(9):2951－2957.

[4]曹文貴,蘇永華,程曄.開(kāi)挖模擬的數(shù)值流形分析方法研究[C]//第八次全國(guó)巖石力學(xué)與工程學(xué)術(shù)大會(huì)論文集.北京:科學(xué)出版社,2004.

[5]朱愛(ài)軍,鄧安富,曾祥勇.數(shù)值流形方法對(duì)巖土工程開(kāi)挖卸荷問(wèn)題的模擬[J].巖土力學(xué),2006,27(2):179－183.ZHU Ai-jun,DENG An-fu,ZENG Xiang-yong.Numerical manifold method for simulation of excavation unloading in geotechnical engineering[J].Rock and Soil Mechanics,2006,27(2):179－183.

[6]王水林,葛修潤(rùn).流形方法在模擬裂紋擴(kuò)展中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào),1997,16(5):405－410 WANG Shui-lin,GE Xiu-yun.Application of manifold method in simulating crack propagation[J].Chinese Journal of Rock Mechanics and Engineering,1997,16(5):405－410

[7]ZHENG H,XU D D.New strategies for some issues of numerical manifold method in simulation of crack propagation[J].International Journal for Numerical Methods in Engineering,2014,97(13):986－1010

[8]姜清輝,周創(chuàng)兵.四面體有限單元覆蓋的三維數(shù)值流形方法[J].巖石力學(xué)與工程學(xué)報(bào),2005,24(24):4455－4460 JIANG Qing-hui,ZHOU Chuang-bing.Three-dimensional numerical manifold method with tetrahedron finite element covers[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(24):4455－4460.

[9]姜清輝,鄧書(shū)申,周創(chuàng)兵.三維高階數(shù)值流形方法研究[J].巖土力學(xué),2006,27(9):1471－1474.JIANG Qing-hui,DENG Shu-shen,ZHOU Chuang-bing.Study of three-dimensional high-order numerical manifold method[J].Rock and Soil Mechanics,2006,27(9):1471－1474.

[10]JIANG Q H,DENG S S,ZHOU C B,et al.Modeling unconfined seepage flow using three-dimensional numerical manifold method[J].Journal of Hydrodynamics,2010,22(4):554－561.

[11]姜清輝,周創(chuàng)兵,張煜.三維數(shù)值流形方法的點(diǎn)-面接觸模型[J].計(jì)算力學(xué)學(xué)報(bào),2006,23(5):569－572.JIANG Qing-hui,ZHOU Chuang-bing,ZHANG Yu.A model of point-to-face contact for three-dimensional numerical manifold method[J].Chinese Journal of Computational Mechanics,2006,23(5):569－572.

[12]CHENG Y M,ZHANG Y H.Formulation of a three-dimensional manifold method with tetrahedron elements[J].Rock Mechanics and Rock Engineering,2008,41(4):601－628.

[13]HE L,AN X M,MA G W.Development of three-dimensional numerical manifold method for jointed rock slope stability analysis[J].International Journal of Rock Mechanics&Mining Sciences,2013,64:22－35.

[14]郭朝旭,鄭宏.高階數(shù)值流形方法中線(xiàn)性相關(guān)問(wèn)題的研究[J].工程力學(xué),2012,29(12):228－232.GUO Chao-xu,ZHENG Hong.Study of linear dependence problem in high-order numerical manifold method[J].Engineering Mechanics,2012,29(12):228－232.

[15]林毅峰,朱合華,蔡永昌.數(shù)值流形方法中線(xiàn)性相關(guān)問(wèn)題的研究[J].計(jì)算力學(xué)學(xué)報(bào),2012,29(5):753－758.LIN Yi-feng,ZHU He-hua,CAI Yong-chang.Research on linear dependence problem of numerical manifold method[J].Chinese Journal of Computational Mechanics,2012,29(5):753－758.

[16]BABUSKA I,OSBORN B U,OSBORN J E.Generalized finite element methods-main ideas,results and perspective[J].International Journal of Computational Methods,2004,1(1):67－103.

[17]朱愛(ài)軍.數(shù)值流形方法研究及其在巖土工程中的應(yīng)用[博士學(xué)位論文D].重慶:重慶大學(xué),2005.