一種基于四棱柱單元三維邊坡網(wǎng)格的計(jì)算方法推導(dǎo)

常天龍

(中國(guó)鐵路設(shè)計(jì)集團(tuán)有限公司，天津 300251)

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

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

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

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

1.1 網(wǎng)格節(jié)點(diǎn)坐標(biāo)計(jì)算

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

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

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

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

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

V=0.25ab[∑A-∑B]

(1)

1.3 滑面所在單元體積計(jì)算方法

1.3.1 四點(diǎn)同層類(lèi)網(wǎng)格

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

(2)

1.3.2 一點(diǎn)凸出類(lèi)網(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 兩點(diǎn)凸出斜交型

判斷條件為：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)交點(diǎn)為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為四點(diǎn)同層型，b為一點(diǎn)凸出型，c為兩點(diǎn)凸出斜交型

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

1.4 滑面單元面積及傾角

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

(9)

(10)

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

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

2 結(jié)語(yǔ)

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

以實(shí)現(xiàn)不同的計(jì)算目的。

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

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

[3]張 煜, 溫國(guó)強(qiáng), 王笑海. 三維體繪制技術(shù)在工程地質(zhì)可視化中的應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報(bào).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]. 巖土力學(xué).2005.26(7): 1123-1126.

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

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

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