多重網格法抗滑穩定計算精度分析

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

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

0 引言

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

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

1 改進多重網格法的實現

1.1 基本思路

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

圖1 滑面與結構網格示意

1.2 計算程序的實現及步驟

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

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

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

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

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

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

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

圖3 面積矢量示意

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

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

2 算例分析

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

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

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

圖4 三維柱體結構尺寸

圖5 三維柱體有限元網格

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

3 精度分析

3.1 高斯點個數

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

3.2 滑面網格的疏密

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

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

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

表2 不同的網格滑面A受力計算結果

4 結論

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

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

［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］楊強，朱玲，薛利軍.基于三維多重網格法的抗滑穩定計算精度分析［J］.巖石力學， 2008， 29（1）:94-100.

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

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

［8］SL 282—2003 混凝土拱壩設計規范［S］.

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

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