常天龍

(中國鐵路設(shè)計集團有限公司，天津 300251)

隨著三維數(shù)值模擬逐漸得到推廣與應(yīng)用，三維地質(zhì)建模技術(shù)快速發(fā)展，在巖土工程和地質(zhì)工程領(lǐng)域中得到了廣泛應(yīng)用，但是現(xiàn)在的地質(zhì)建模仍是主要依靠一些商用軟件如：Midas、ANSYS、ABAQUS等軟件中提供的建模模塊實現(xiàn)三維建模。人們起初對三維地質(zhì)建模的要求是實現(xiàn)建立的三維地質(zhì)模型可以真實反映出研究對象包括的地層及不良地質(zhì)體，而隨著數(shù)值模擬的發(fā)展，為了提高數(shù)值計算結(jié)果的精確度，要求所建立的模型不僅要能真實反映研究對象的地質(zhì)特征，而且模型還要具有較高質(zhì)量的網(wǎng)格剖分度[1]。為了實現(xiàn)建模的效率與計算的準確性，往往為了突出研究的問題而采取了部分“有重點”的建模方法。

在三維數(shù)值模擬地質(zhì)建模方面，許多學者提出了自己的方法，Jones提出了利用鉆孔與剖面結(jié)合的實體建模方法為MODFLOW提供計算網(wǎng)格[4]，王春祥等提出L-W拓撲模型建立三維地質(zhì)模型并實現(xiàn)地質(zhì)信息與有限元計算的結(jié)合[6]，王明華等通過研究層狀巖體的分布特征，提出基于松散模式的三維規(guī)則網(wǎng)格與FLAC3D之間的模型轉(zhuǎn)化方法[5]。李明超等結(jié)合數(shù)值模擬計算的特點提出地質(zhì)建模與數(shù)值模擬耦合的方法來進行前處理的思路[7]。以上研究成果為三維地質(zhì)建模服務(wù)于數(shù)值模擬計算提供了很多新思路、新方法，為未來通用程序的開發(fā)提供了不同的思路。

不同學者開發(fā)的計算程序多將滑面近似圓弧形去處理，對于某些非圓弧滑面邊坡往往很難得到準確的計算結(jié)果，本文結(jié)合前人的研究成果，提出了一種基于四棱柱網(wǎng)格單元的建模方法，結(jié)合Matlab矩陣計算的特點，實現(xiàn)數(shù)據(jù)的矩陣計算與存儲。依據(jù)地質(zhì)體地層呈層分布特點，本文采用矩陣計算插分得到網(wǎng)格坐標，通過單元節(jié)點三維坐標間的運算求取四棱柱的體積、單元表面積等參數(shù)，從而建立一種基于四棱柱單元可以求取任意形狀滑面邊坡穩(wěn)定性的地質(zhì)模型。

1 構(gòu)建三維地質(zhì)模型

1.1 網(wǎng)格節(jié)點坐標計算

通過現(xiàn)場測量得到邊坡表面地面點及勘探鉆孔坐標，將同一地層內(nèi)的控制點三維坐標放入同一文本文件中作為插值點備用。Matlab中提供了四種插值函數(shù)1.nearest函數(shù)：最鄰近插值方法(nearest neighbor interpolation)；2.linear函數(shù)：線性插值(linear interpolation)；3.spline函數(shù)：三次樣條插值(cubic spline interpolation)；4. v5cubic函數(shù)：Matlab程序自帶插值計算函數(shù)。根據(jù)不同情況，可以根據(jù)地形起伏變化的大小選用不同的插值函數(shù)進行網(wǎng)格插值劃分，筆者在程序設(shè)計過程中對比不同的插值函數(shù)結(jié)果發(fā)現(xiàn)使用“v5cubic”函數(shù)進行插值可以得到比較平滑的曲面，減少表面凸點的出現(xiàn)，使得數(shù)據(jù)奇點更少，提高計算精度。使用自編的插值程序，各個地層按照自定義網(wǎng)格尺寸插值計算得出各個曲面的網(wǎng)格角點三維坐標。由于Matlab是矩陣運算程序，所有插值得到的網(wǎng)格角點坐標以矩陣的形勢保存，方便讀取修改(見圖1)。

圖1 插值得到的曲面網(wǎng)格

1.2 網(wǎng)格單元數(shù)據(jù)計算

在同一邊坡不同位置上滑面所在地層不盡相同，計算某滑面以上不同地層的總重量既為該滑面單元的滑體重量，對于滑面所在單元，滑面切割該單元為兩部分，滑面以上部分計入滑體重量，滑面以下需要剔除。針對不同的滑面與四棱柱單元相交情況可把單元分為三類12種單元。筆者按照滑面四個角點的相對位置分為：“四點同層”，“一點凸出”，“兩點凸出斜交”三種單元類型。(圖2中a為四點同層型，b為一點凸出型，c為兩點凸出斜交型)從圖中可以看出，假設(shè)滑面為H層，A、B、C層為地層面，H層平面矩形單元角點坐標沿順時針方向分別為H(i,j)、H(i+1,j)、H(i+1,j+1)、H(i,j+1)其他層角點坐標角標相同不予贅述，每個層面網(wǎng)格坐標構(gòu)成規(guī)則的坐標矩陣。為了滿足不出現(xiàn)滑面跨層現(xiàn)象，網(wǎng)格尺寸需滿足：maxΔH﹤2h或2h/a>tanθ且2h/b>tanθ，其中a、b為網(wǎng)格長寬，h為分層層厚；θ為滑面最大傾角；為滑面相鄰兩個網(wǎng)格角點高程之差；只需要將網(wǎng)格劃分足夠細，避免出現(xiàn)網(wǎng)格跨度過大從而導致跨層網(wǎng)格出現(xiàn)即可。下面就對上述三類網(wǎng)格體積計算公式進行推導：

未被滑面切割的一般單元體積為：

V=0.25ab[∑A-∑B]

(1)

1.3 滑面所在單元體積計算方法

1.3.1 四點同層類網(wǎng)格

判斷條件為：B(i，j)

(2)

1.3.2 一點凸出類網(wǎng)格

判斷條件為：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為兩點凸出斜交型

以上給出了三類網(wǎng)格某種情況下的體積計算公式，其他同類型網(wǎng)格只需要調(diào)整各個角點的在公示中的順序即可得到相應(yīng)的計算公式，在此不予以一一贅述。通過Matlab中循環(huán)判斷語句實現(xiàn)上文中的判斷計算過程，計算出各個網(wǎng)格的體積后乘以該單元密度即可得到單元重度。

1.4 滑面單元面積及傾角

在求解安全系數(shù)時還需要輸入每個滑面網(wǎng)格的表面積，滑面在XOY面上的投影即為a*b(a、b分別為設(shè)定的網(wǎng)格長寬)，通過求取每個滑面單元與水平面的夾角即可求得滑面的面積。面積和夾角公式為式9、10。

(9)

(10)

式中：α為滑面與水平面夾角

通過本節(jié)中各項計算公式，可以得到穩(wěn)定分析中所需的單元重量，滑面傾角，滑面面積等所需數(shù)據(jù)，從而實現(xiàn)三維地質(zhì)建模，數(shù)據(jù)以矩陣形式存儲與調(diào)用。

2 結(jié)語

Matlab是基于矩陣計算的一種用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計算的高級技術(shù)計算語言和交互式算法開發(fā)程序，本文介紹了一種基于Matlab自編程序?qū)崿F(xiàn)地層網(wǎng)格化的方法，從而得到網(wǎng)格節(jié)點坐標，經(jīng)過文中推導的一系列坐標運算公式，可以快速的實現(xiàn)網(wǎng)格單元特性計算，從而實現(xiàn)數(shù)值分析前處理三維地質(zhì)建模的目的，本方法建立的三維地質(zhì)模型不受滑面形狀所限，適用于任意形狀滑面邊坡計算，與邊坡實際地層較吻合。該算法具有較好的實用性和可拓展性，通過自編不同功能的算法可

以實現(xiàn)不同的計算目的。

[1]徐能雄. 適于數(shù)值模擬的三維工程地質(zhì)建模方法[J]. 巖土工程學報.2009.31(11): 1710-1715.

[2]胡斌,張倬元,黃潤秋,等. FLAC3D 前處理程序的開發(fā)與仿真效果檢驗[J]. 巖石力學與工程學報.2002.21(9):1381-1397.

[3]張 煜, 溫國強, 王笑海. 三維體繪制技術(shù)在工程地質(zhì)可視化中的應(yīng)用[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]王明華, 白 云. 層狀巖體三維可視化構(gòu)模與數(shù)值模擬的集成研究[J]. 巖土力學.2005.26(7): 1123-1126.

[6]王純祥, 白世偉. 三維地層信息系統(tǒng)與有限元方法集成研究[J]. 巖石力學與工程學報.2004.23(21): 3695-3699.

[7]李明超, 鐘登華, 等. 基于三維地質(zhì)模型的工程巖體結(jié)構(gòu)精細數(shù)值建模[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]徐能雄, 武 雄, 汪小剛,等. 基于三維地質(zhì)建模的復雜構(gòu)造巖體六面體網(wǎng)格剖分方法[J]. 巖土工程學報.2006.28(8): 957-961 .

[10]谷天峰, 王家鼎, 王念秦. 呂梁機場黃土滑坡特征及其三維穩(wěn)定性分析[J]. 巖土力學.2013.