胡海周，鄭金城，彭 剛，梁春華

（三峽大學土木與建筑學院，湖北 宜昌 443002）

0 引言

多重網(wǎng)格法從興起至今已有近50年了，此方法的研究與使用起初大多是數(shù)學界的專家[1-4]，楊強[5-7]等最先將多重網(wǎng)格法應用于工程結(jié)構(gòu)穩(wěn)定計算，采用剛體極限平衡方法計算得到的安全系數(shù)與現(xiàn)行規(guī)范[8]一致，計算結(jié)果易被設(shè)計部門理解與應用。因此，多重網(wǎng)格法開始被廣泛應用于工程結(jié)構(gòu)計算分析中。楊強[5]提出的多重網(wǎng)格法是:先將滑面重新剖分網(wǎng)格，搜索滑面中與某一節(jié)點距離最近的一個高斯點，取高斯點所在單元的所有高斯點對節(jié)點進行應力插值[9-10]，可得滑面上各節(jié)點應力和滑面上總應力，利用抗剪強度或抗剪斷公式計算安全系數(shù)。

楊強[5]提出的多重網(wǎng)格法只能考慮滑面P中與節(jié)點A距離最近的高斯點B所在的單元M內(nèi)的應力對A節(jié)點的貢獻，當與滑面節(jié)點A相對第二近的高斯點C不在單元M內(nèi)時，C點的應力將不能對A點的應力作貢獻，即它不能完全考慮與滑面節(jié)點第二近的高斯點應力對滑面應力所作的貢獻?；诖?，在楊強提出的理論基礎(chǔ)上對多重網(wǎng)格法進行了改進，搜索與滑面P中節(jié)點A距離相對較近的n個高斯點對A點進行應力映射，從而考慮了所有相對較近的高斯點應力對滑面節(jié)點應力的貢獻。

1 改進多重網(wǎng)格法的實現(xiàn)

1.1 基本思路

多重網(wǎng)格有有限元結(jié)構(gòu)模型網(wǎng)格和滑面網(wǎng)格。有限元分析中，在每個單元的高斯積分點處計算應力。三維有限元的應力成果與模型的結(jié)構(gòu)網(wǎng)格有關(guān)，本文采用多重網(wǎng)格法將三維有限元高斯積分點應力結(jié)果轉(zhuǎn)移到滑動面上。將滑面重新剖分網(wǎng)格 （見圖1），滑面上各個節(jié)點的應力由結(jié)構(gòu)網(wǎng)格上高斯積分點的應力進行插值得到。對于滑面P中的A節(jié)點，搜索與其距離較近的n個高斯點B1、B2、B3……Bn，然后將n個高斯點上的應力值進行插值得到A點應力狀態(tài)，公式為

圖1 滑面與結(jié)構(gòu)網(wǎng)格示意

1.2 計算程序的實現(xiàn)及步驟

運用ABAQUS軟件進行有限元分析，并結(jié)合FORTRAN語言共同實現(xiàn)多重網(wǎng)格法。多重網(wǎng)格法程序編制主要步驟如下:

（1）在ABAQUS中提取有限元分析后的所有單元高斯積分點坐標及對應的6個應力分量。

（2）給滑移面P重新劃分網(wǎng)格，提取滑面所有點坐標和組成單元的節(jié)點編號。

（3）假定滑移面其中一個節(jié)點A，搜尋與其較近的n個高斯積分點B1、B2、…、Bn。分別計算A點到n個高斯積分點的距離L1、L2、…、Ln，按照公式（2）計算權(quán)函數(shù)Sk。圖2為較近的高斯積分點個數(shù)為8時的權(quán)函數(shù)示意圖。由公式（1）計算A點的應力分量和滑面上其他節(jié)點所對應的應力分量。

（4）將滑面P上每個四邊形單元分成2個三角形 （見圖3），分別計算三角形①和②面積。以三角形①為例，其面積矢量為Aj=

（5）計算單元的力矢量。任意1個三角形上的力矢量Fi為

式中，Aj為某單元中1個三角形的面積矢量；為某單元中1個三角形所有節(jié)點應力矢量的平均值，其中，。單個四邊形單元的力矢量則為2個三角形力矢量之和。

圖3 面積矢量示意

（6）將滑面上各單元的法向矢量求和再平均，可得滑面的平均法向矢量，由各單元的所有力矢量求得合矢量，計算合矢量在滑面的法向合力和切向合力?；娴拿姘踩菿為

式中，f、c為滑面的剪摩系數(shù)；R、S分別為滑面的法向力與切向力。

2 算例分析

采用1個三維柱體進行有限元分析 （見圖4），柱體高l1=10 m，l2=5 m，l3=9 m，長h和寬b均為5 m。柱體的彈性模量E=20 GPa，泊松比μ=0.2，剪摩系數(shù)f=0.68、 c=0.6 MPa。 柱體重度γ=24 kN/m3。 其中，A為柱體中的一個面，假定為滑面。在z=10 m端部采用固支約束，考慮自重荷載。采用ABAQUS軟件進行模擬，單元為2000個，節(jié)點為2541個。有限元模型見圖5。

三維柱體僅受重力作用，塊體內(nèi)體力為0，滑面A上方塊體重力理論即為滑面A所受的合力

式中，F(xiàn)n為滑面A法向力；Fτ為滑面A切向力，方向為垂直于Y軸沿滑面向下；Vu為滑面A上方塊體體積。

圖4 三維柱體結(jié)構(gòu)尺寸

圖5 三維柱體有限元網(wǎng)格

由式（5）計算得到三維柱體內(nèi)應力分布并投影至滑面A上，可得A面上受力情況，可將此解近似為理論解，計算得滑面A的理論合力為4.2000 MN，法向力為3.2797 MN，垂直于y軸向下的切向力為2.6237 MN。

3 精度分析

3.1 高斯點個數(shù)

本次模擬采用C3D8單元類型進行分析，整個模型的高斯積分點的個數(shù)為16000個?；婢W(wǎng)格方案采用10×12。高斯點個數(shù)分別取為1到10時，滑面A的計算結(jié)果見表1。由表1可知，采用多重網(wǎng)格法進行插值的結(jié)果與理論解比較接近，驗證了多重網(wǎng)格法程序的正確性。隨著高斯點個數(shù)的增加，滑面A所受力的誤差率均是先減小后增大，當高斯點個數(shù)為5時，誤差率最小，誤差率均控制在1%以下。

3.2 滑面網(wǎng)格的疏密

對滑面網(wǎng)格進行重新劃分，選取6×8、8×10、10×12、 12×16、 15×20等 5套網(wǎng)格進行分析。 在計算分析時，高斯點個數(shù)取為5，單元類型為C3D8，計算結(jié)果見表2。

由表2可知，在5套網(wǎng)格方案中，以10×12網(wǎng)格方案最接近理論解，其誤差率最小。分析其原因是:10×12網(wǎng)格方案的滑面網(wǎng)格的密度和三維柱體有限元模型中的網(wǎng)格分布密度最接近。因此，當采用最近的高斯點進行插值時，誤差會相對較小。

表1 不同高斯點下滑面A受力計算結(jié)果

表2 不同的網(wǎng)格滑面A受力計算結(jié)果

4 結(jié)論

（1）利用改進的多重網(wǎng)格法所得的結(jié)果與理論解比較接近，誤差率均控制在1%以下，且當高斯點個數(shù)為5時結(jié)果最為接近理論值。

（2）滑面網(wǎng)格的疏密程度與有限元模型網(wǎng)格越接近，計算結(jié)果越理想。

［1］BRANDT A.Multi-Level Adaptive Solutions［J］.Mathematics of Computation,1977,31（138）:333-390.

［2］ZHANG J.Fast and High Accuracy Multigrid Solution of the Three Dimensional Poisson Equation［J］.Journal of Computational Physics，1998,143（2）:449-461.

［3］BRAESS D,VERFURTH R.Multigrid methods for nonconforming finite element methods［J］.SIAM Journal on Numerical Analysis， 1990,27（4）:979-986.

［4］BRENNER S C.An optimal-order nonconforming multigrid method for the biharmonic equation［J］.SIAM Journal on Numerical Analysis，1989,26（5）:1124-1138.

［5］楊強，朱玲，薛利軍.基于三維多重網(wǎng)格法的抗滑穩(wěn)定計算精度分析［J］.巖石力學， 2008， 29（1）:94-100.

［6］楊強，朱玲，薛利軍.基于三維多重網(wǎng)格法的極限平衡法在錦屏高邊坡穩(wěn)定性分析中的應用［J］.巖石力學與工程學報，2005，24（A02）:5313-5318.

［7］楊強，朱玲，翟明杰.基于三維非線性有限元的壩肩穩(wěn)定剛體極限平衡法機理研究［J］.巖石力學與工程學報，2005，24（19）:3403-3409.

［8］SL 282—2003 混凝土拱壩設(shè)計規(guī)范［S］.

［9］王潤英.計算三維溫度場及溫度應力的徑向點插值無網(wǎng)格法［J］.三峽大學學報， 2009， 31（2）:18-22.

［10］董曉華，薄會娟，鄧霞，等.降雨空間插值方法及在清江流域的應用［J］.三峽大學學報， 2009， 31（6）:6-10.