胡海周，郑金城，彭 刚，梁春华

（三峡大学土木与建筑学院，湖北 宜昌 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）将滑面上各单元的法向矢量求和再平均，可得滑面的平均法向矢量，由各单元的所有力矢量求得合矢量，计算合矢量在滑面的法向合力和切向合力。滑面的面安全度K为

式中，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个。滑面网格方案采用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.