一种基于四棱柱单元三维边坡网格的计算方法推导

常天龙

(中国铁路设计集团有限公司，天津 300251)

随着三维数值模拟逐渐得到推广与应用，三维地质建模技术快速发展，在岩土工程和地质工程领域中得到了广泛应用，但是现在的地质建模仍是主要依靠一些商用软件如：Midas、ANSYS、ABAQUS等软件中提供的建模模块实现三维建模。人们起初对三维地质建模的要求是实现建立的三维地质模型可以真实反映出研究对象包括的地层及不良地质体，而随着数值模拟的发展，为了提高数值计算结果的精确度，要求所建立的模型不仅要能真实反映研究对象的地质特征，而且模型还要具有较高质量的网格剖分度[1]。为了实现建模的效率与计算的准确性，往往为了突出研究的问题而采取了部分“有重点”的建模方法。

在三维数值模拟地质建模方面，许多学者提出了自己的方法，Jones提出了利用钻孔与剖面结合的实体建模方法为MODFLOW提供计算网格[4]，王春祥等提出L-W拓扑模型建立三维地质模型并实现地质信息与有限元计算的结合[6]，王明华等通过研究层状岩体的分布特征，提出基于松散模式的三维规则网格与FLAC3D之间的模型转化方法[5]。李明超等结合数值模拟计算的特点提出地质建模与数值模拟耦合的方法来进行前处理的思路[7]。以上研究成果为三维地质建模服务于数值模拟计算提供了很多新思路、新方法，为未来通用程序的开发提供了不同的思路。

不同学者开发的计算程序多将滑面近似圆弧形去处理，对于某些非圆弧滑面边坡往往很难得到准确的计算结果，本文结合前人的研究成果，提出了一种基于四棱柱网格单元的建模方法，结合Matlab矩阵计算的特点，实现数据的矩阵计算与存储。依据地质体地层呈层分布特点，本文采用矩阵计算插分得到网格坐标，通过单元节点三维坐标间的运算求取四棱柱的体积、单元表面积等参数，从而建立一种基于四棱柱单元可以求取任意形状滑面边坡稳定性的地质模型。

1 构建三维地质模型

1.1 网格节点坐标计算

通过现场测量得到边坡表面地面点及勘探钻孔坐标，将同一地层内的控制点三维坐标放入同一文本文件中作为插值点备用。Matlab中提供了四种插值函数1.nearest函数：最邻近插值方法(nearest neighbor interpolation)；2.linear函数：线性插值(linear interpolation)；3.spline函数：三次样条插值(cubic spline interpolation)；4. v5cubic函数：Matlab程序自带插值计算函数。根据不同情况，可以根据地形起伏变化的大小选用不同的插值函数进行网格插值划分，笔者在程序设计过程中对比不同的插值函数结果发现使用“v5cubic”函数进行插值可以得到比较平滑的曲面，减少表面凸点的出现，使得数据奇点更少，提高计算精度。使用自编的插值程序，各个地层按照自定义网格尺寸插值计算得出各个曲面的网格角点三维坐标。由于Matlab是矩阵运算程序，所有插值得到的网格角点坐标以矩阵的形势保存，方便读取修改(见图1)。

图1 插值得到的曲面网格

1.2 网格单元数据计算

在同一边坡不同位置上滑面所在地层不尽相同，计算某滑面以上不同地层的总重量既为该滑面单元的滑体重量，对于滑面所在单元，滑面切割该单元为两部分，滑面以上部分计入滑体重量，滑面以下需要剔除。针对不同的滑面与四棱柱单元相交情况可把单元分为三类12种单元。笔者按照滑面四个角点的相对位置分为：“四点同层”，“一点凸出”，“两点凸出斜交”三种单元类型。(图2中a为四点同层型，b为一点凸出型，c为两点凸出斜交型)从图中可以看出，假设滑面为H层，A、B、C层为地层面，H层平面矩形单元角点坐标沿顺时针方向分别为H(i,j)、H(i+1,j)、H(i+1,j+1)、H(i,j+1)其他层角点坐标角标相同不予赘述，每个层面网格坐标构成规则的坐标矩阵。为了满足不出现滑面跨层现象，网格尺寸需满足：maxΔH﹤2h或2h/a>tanθ且2h/b>tanθ，其中a、b为网格长宽，h为分层层厚；θ为滑面最大倾角；为滑面相邻两个网格角点高程之差；只需要将网格划分足够细，避免出现网格跨度过大从而导致跨层网格出现即可。下面就对上述三类网格体积计算公式进行推导：

未被滑面切割的一般单元体积为：

V=0.25ab[∑A-∑B]

(1)

1.3 滑面所在单元体积计算方法

1.3.1 四点同层类网格

判断条件为：B(i，j)

(2)

1.3.2 一点凸出类网格

判断条件为：H(i,j)>B(i,j),H(i,j+1)

六棱柱A(i,j)A(i+1,j)A(i+1,j+1)A(i,j+1)H(i,j)H(i+1,j)H(i+1,j+1)H(i,j+1)体积为：

(3)

滑面H(i,j)H(i+1,j)H(i+1,j+1)H(i,j+1)面积为：

(4)

1.3.3 两点凸出斜交型

判断条件为：B(i,j)>H(i,j),B(i,j+1)>H(i,j+1),B(i+1,j)

滑面单元H(i,j)H(i+1,j)H(i+1,j+1)H(i,j+1)与单元面B(i,j)B(i+1,j)B(i+1,j+1)B(i,j+1)交点为O1、O2。

V=VAB-Vu+Vd

(5)

其中：VAB为滑面以上A、B单元体积，Vu为上三棱柱H(i+1,j)H(i+1,j+1)B(i+1,j)B(i+1,j+1)O1O2体积，Vd位下三棱柱H(i,j)H(i,j+1)B(i,j)B(i,j+1)O1O2体积

V=0.25ab[∑A-∑B]

(6)

(7)

(8)

a为四点同层型，b为一点凸出型，c为两点凸出斜交型

以上给出了三类网格某种情况下的体积计算公式，其他同类型网格只需要调整各个角点的在公示中的顺序即可得到相应的计算公式，在此不予以一一赘述。通过Matlab中循环判断语句实现上文中的判断计算过程，计算出各个网格的体积后乘以该单元密度即可得到单元重度。

1.4 滑面单元面积及倾角

在求解安全系数时还需要输入每个滑面网格的表面积，滑面在XOY面上的投影即为a*b(a、b分别为设定的网格长宽)，通过求取每个滑面单元与水平面的夹角即可求得滑面的面积。面积和夹角公式为式9、10。

(9)

(10)

式中：α为滑面与水平面夹角

通过本节中各项计算公式，可以得到稳定分析中所需的单元重量，滑面倾角，滑面面积等所需数据，从而实现三维地质建模，数据以矩阵形式存储与调用。

2 结语

Matlab是基于矩阵计算的一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式算法开发程序，本文介绍了一种基于Matlab自编程序实现地层网格化的方法，从而得到网格节点坐标，经过文中推导的一系列坐标运算公式，可以快速的实现网格单元特性计算，从而实现数值分析前处理三维地质建模的目的，本方法建立的三维地质模型不受滑面形状所限，适用于任意形状滑面边坡计算，与边坡实际地层较吻合。该算法具有较好的实用性和可拓展性，通过自编不同功能的算法可

以实现不同的计算目的。

[1]徐能雄. 适于数值模拟的三维工程地质建模方法[J]. 岩土工程学报.2009.31(11): 1710-1715.

[2]胡斌,张倬元,黄润秋,等. FLAC3D 前处理程序的开发与仿真效果检验[J]. 岩石力学与工程学报.2002.21(9):1381-1397.

[3]张 煜, 温国强, 王笑海. 三维体绘制技术在工程地质可视化中的应用[J]. 岩石力学与工程学报.2002.21(4): 563-567.

[4]LEMONAM, JONES. Building solid models from boreholes and user-defined cross-sections[J]. Computers & Geosciences, 2003.29(5): 547-555.

[5]王明华, 白 云. 层状岩体三维可视化构模与数值模拟的集成研究[J]. 岩土力学.2005.26(7): 1123-1126.

[6]王纯祥, 白世伟. 三维地层信息系统与有限元方法集成研究[J]. 岩石力学与工程学报.2004.23(21): 3695-3699.

[7]李明超, 钟登华, 等. 基于三维地质模型的工程岩体结构精细数值建模[J]. 岩石力学与工程学报.2007.26(9): 1893-1898.

[8]Xing Z. Three-dimensional stability analysis of concave slopes in plan view[J]. Journal of the Geotechnical Engineering, 1988.114(6): 658-671.

[9]徐能雄, 武 雄, 汪小刚,等. 基于三维地质建模的复杂构造岩体六面体网格剖分方法[J]. 岩土工程学报.2006.28(8): 957-961 .

[10]谷天峰, 王家鼎, 王念秦. 吕梁机场黄土滑坡特征及其三维稳定性分析[J]. 岩土力学.2013.