蒙鐸,陳志文,黃俊,張菖(.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都 6756; .成都理工大學(xué)管理科學(xué)學(xué)院，成都 60059；.成都理工大學(xué)環(huán)境與土木工程學(xué)院，成都 60059)

一種基于CATIA軟件進(jìn)行地質(zhì)建模與有限元計算的前處理方法及其應(yīng)用

蒙鐸1,陳志文2,黃俊1,張菖3
(1.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都 611756; 2.成都理工大學(xué)管理科學(xué)學(xué)院，成都 610059；3.成都理工大學(xué)環(huán)境與土木工程學(xué)院，成都 610059)

使用CATIA軟件所建立的復(fù)雜地質(zhì)模型，很難直接用于有限元計算，需要對所建立的模型進(jìn)行前處理。提出了一種前處理方法：應(yīng)用巖土參數(shù)加權(quán)平均法，并結(jié)合工程地質(zhì)分析、有限元計算原理和拓?fù)鋵W(xué)分析方法，對CAITA軟件所建立的復(fù)雜地質(zhì)模型進(jìn)行前處理。該前處理方法既能近似地反映出地質(zhì)條件，又能達(dá)到有限元分析的要求，從而可保證分析結(jié)果的準(zhǔn)確可靠性。通過對四川某滑坡的地質(zhì)建模與有限元計算，對該前處理方法進(jìn)行了實例驗證，取得了良好效果。

地質(zhì)建模; 有限元計算; 拓?fù)浞治? 前處理

1 引言

CATIA雖然不是專業(yè)的地質(zhì)建模軟件，但是由于其強大的曲面造型功能和三維協(xié)同設(shè)計功能已經(jīng)在行業(yè)內(nèi)使用，特別是水電行業(yè)對CATIA三維的協(xié)同設(shè)計應(yīng)用已經(jīng)有了一個較高的水平[1-3]。CAITA建立的地質(zhì)模型相對于有限元分析來說過于復(fù)雜，但是直接使用CAE軟件建立地質(zhì)模型，又太過于復(fù)雜和困難，而且三維可視化程度低，因此需要對模型進(jìn)行必要的前處理。筆者將地質(zhì)模型分層進(jìn)行參數(shù)加權(quán)平均處理，并簡化了模型，使模型能導(dǎo)入ANSYS進(jìn)行計算，并且計算結(jié)果接近實際工況。

國內(nèi)已經(jīng)有很多關(guān)于CAITA地質(zhì)建模的論文和書籍[4]，筆者在本文中就不對CATIA如何建立地質(zhì)模型做過多的介紹，只是對四川某滑坡實例容做一些簡介。

2 滑坡地形地貌

該滑坡實例所處地域?qū)贅?gòu)造剝蝕中低山斜坡地貌，地勢西南高而東北低?；挛挥谏襟w中下部，整體地形坡度約為15°，局部坡度較大可達(dá)50°。整個滑坡在縱向上略顯階狀，呈條帶狀小平臺與斜坡相間狀地形。在滑坡后緣見明顯陡坎，坎高約2～5 m，坡度約40°，在勘查期間因施工影響，后緣階坎被破壞并見一人工開挖而成的溝壑。在滑坡中部西側(cè)因人工堆砌的形成一高約4 m的土質(zhì)平臺。在滑坡前緣滑體物質(zhì)大量堆積，現(xiàn)因前緣施工開挖多形成陡坡。

3 CATIA地質(zhì)模型

根據(jù)勘察資料，建立了如圖1的地質(zhì)三維地質(zhì)模型。

圖1 滑坡三維地質(zhì)模型

筆者將選出中風(fēng)化灰?guī)r模型(圖2)為讀者展示CATIA三維地質(zhì)模型直接導(dǎo)入ANSYS的困難。

圖2 中風(fēng)化灰?guī)r

從圖2可以看到建立該模型用了大量的網(wǎng)格，而這些大量的復(fù)雜的網(wǎng)格在導(dǎo)入ANSYS的時候會引起ANSYS的直接崩潰。這也就是為什么在用ANSYS分析CATIA三維地質(zhì)模型的時候需要對CATIA模型做一定的簡化的主要原因。下面筆者將先展示簡化后的模型(圖3)，再說明簡化的依據(jù)和簡化方法。

圖3 簡化后模型

4 分層參數(shù)加權(quán)平均法簡化三維地質(zhì)模型

勘察資料顯示滑坡覆蓋體是以前發(fā)生過滑坡的滑體，整個滑坡覆蓋體上支離破碎，地表大量發(fā)育因拉裂、錯動形成的裂縫。在滑坡覆蓋上還覆蓋有一部分的人工堆積體，這一部分的人工堆積體是滑坡體的一部分。加之人工堆積體和滑坡覆蓋體的參數(shù)接近，所以講這兩部分合并，是可行的。以下是該滑坡巖土體物理力學(xué)參數(shù)建議(表1)。

表1 巖土體物理力學(xué)參數(shù)建議值表

4.1 有限元分析原理對分層參數(shù)平均法模型的簡化優(yōu)點的說明

有限元法是把連續(xù)體結(jié)構(gòu)離散，近似地用有限個離散體(單元)代替連續(xù)結(jié)構(gòu)。然后建立約束方程等邊界條件簡化或等效，利用應(yīng)變位移關(guān)系和本構(gòu)關(guān)系，將平衡方程轉(zhuǎn)化為線性方程組來求解。有限元計算的步驟我們可以看出，與模型結(jié)構(gòu)緊密相關(guān)的是結(jié)構(gòu)離散化。

結(jié)構(gòu)區(qū)域離散的方法也就是通常所說的網(wǎng)格劃分來說，網(wǎng)格的劃分既是力學(xué)的標(biāo)準(zhǔn)，也是美學(xué)的標(biāo)準(zhǔn)。由于CATIA建模的方式和實際的地質(zhì)情況，會有一些角度較小的尖角，這樣的網(wǎng)格劃分出來是不能滿足有限元分析的要求。小角度的網(wǎng)格，會使某些點在使用線性求解位移值的時候帶來影響計算結(jié)果的誤差，甚至引起結(jié)果不收斂，無法完成有限元計算。分層參數(shù)平均法，分層后的各地質(zhì)體，這就避免了實際情況或者由于建模誤差產(chǎn)生的模型小角度[5、6]。

4.2 利用拓?fù)鋵W(xué)簡化模型

由于CAD和CAE軟件使用要求的不同，CAD軟件要將所設(shè)計或者是想要展現(xiàn)的東西以盡量詳細(xì)的方式表現(xiàn)出來，而CAE軟件的要求是簡單、明了，能夠以最簡單的方式來展現(xiàn)出所需要計算的模型。

一個三維多邊形表面模型由其表面上的頂點組成的,通過其拓?fù)浣Y(jié)構(gòu)和幾何位置兩方面的信息來描述,其中拓?fù)浣Y(jié)構(gòu)是指各頂點之間的相互連接關(guān)系,而幾何位置是指各頂點在三維空間中的坐標(biāo)信息[7，8]。在利用CAITA的曲面生成功能生成的曲面會由方格片構(gòu)成，這會使模型復(fù)雜程度增加，而是由于在建模時，為了查看滑坡主軸和副軸斷面，要提前對模型進(jìn)行切割，這就可能產(chǎn)生微小面或者小角度。在圖4、圖5中我們可以看到，由于分割模型的原因，產(chǎn)生了微小角度和微小面，這樣的模型形式會給ANSYS在進(jìn)行網(wǎng)格劃分的時候帶來極大的工作量，甚至網(wǎng)格劃分失敗。

圖4 未使用曲面簡化的模型

圖5 使用曲面簡化后的模型

使用CAITA創(chuàng)成式外形設(shè)計模塊中的曲面簡化功能，使曲面簡化，這樣就避免了上述建模方式對ANSYS網(wǎng)格劃分帶來的影響。進(jìn)一步來講，由于這樣的簡化，使得實體成為由4個面組成的單元，便可以使用映射網(wǎng)格劃分。

綜述：使用分層參數(shù)加權(quán)平均法，將地質(zhì)體的主要地質(zhì)體分層，本文算例是以強風(fēng)化灰?guī)r的2個分界面為分界，將整體的地質(zhì)模型分為3層。這樣的分層方法既使模型變的較為簡單，又避免了由于CATIA建模方式所帶來的工作量大或者模型無法識別的影響。

4.3 操作方法

(1) 分層

根據(jù)將巖土體所占的比例和地層分界面作為分層依據(jù)，使用CAITA的布爾操作將其合并。

(2) 切割

主要是根據(jù)滑坡主軸和副軸斷面進(jìn)行切割，以便于查看斷面的情況。

(3) 參數(shù)處理

式中，xi為i巖土體的某一項參數(shù)；Vi為i巖土體的體積。

筆者在結(jié)合該滑坡的勘查報告和《巖石力學(xué)參數(shù)手冊》，對本次模擬需要用到的參數(shù)做了一個統(tǒng)計(表2)。

表2 巖土體物理力學(xué)參數(shù)計算表

需要說明的是，由于《巖石力學(xué)參數(shù)手冊》中，提到的泥巖與灰混合的的粘聚力值是取為310 kPa，所以筆者在此將其取值為310 kPa?？辈橘Y料中的泥巖和灰?guī)r的內(nèi)摩擦角過高，筆者在《巖石力學(xué)參數(shù)手冊》中并未發(fā)現(xiàn)有如此大的內(nèi)摩擦角值。最下部地層由于鉆孔為涉及該區(qū)域，又不影響滑坡穩(wěn)定性，所以筆者將其簡化彈性體。

5 有限元計算

使用solid45單元模擬土體，參數(shù)選用分層加權(quán)平均處理的結(jié)果，對邊走向方向和兩側(cè)分別固定y和x方向位移，對模型最底部的面固定x,y,z三個方向的位移來進(jìn)行求解。使用系數(shù)折減法求得滑坡剪出口的位置。

ANSYS中的模型：見圖6。

圖6 導(dǎo)入ANSYS后的模型

劃分后的網(wǎng)格：見圖7。

圖7 劃分網(wǎng)格的結(jié)果

折減系數(shù)1時的結(jié)果：

依次為變形圖(圖8)，y方向位移圖(圖9)和塑性應(yīng)變圖(圖10)。

圖8 整體變形圖

圖9 y方向位移圖

圖10 塑性應(yīng)變圖

一直到折減系數(shù)為3.0時候，ANSYS提示計算結(jié)果不收斂。下面筆者將展示折減系數(shù)為2.8時候結(jié)果：

總體的變形(圖11)，y方向位移(圖12)，塑性應(yīng)變(圖13)。

圖11 折減系數(shù)2.8時總體變形圖

圖12 折減系數(shù)2.8時y方向位移圖

圖13 折減系數(shù)2.8時塑性變形圖

主軸斷面情況：見圖14、圖15。

圖14 折減系數(shù)2.8主軸斷面y方向位移圖

圖15 折減系數(shù)2.8主軸斷面塑性應(yīng)變圖

實際勘查結(jié)果：見圖16。

圖16 勘查資料顯示的滑坡影響范圍

紅色的虛線是勘查的滑坡范圍。

對比可以看出，模擬的結(jié)果和實際勘查結(jié)果在地表的剪出口位置很接近。這一點可以說明分層巖土參數(shù)加權(quán)平均法還是能較為準(zhǔn)確的模擬出滑坡剪出口的位置。y方向的位移的趨勢和通常對滑坡的認(rèn)識一致沒在剪出口位置最大[9-12]。塑性位移最大的位置，在勘查資料中表明，是強風(fēng)化的泥巖破壞的位置，但是由于勘查資料所給的泥巖、灰?guī)r的內(nèi)摩擦角太大的原因，筆者又將泥巖和灰?guī)r的內(nèi)摩擦角取小，在計算結(jié)果不收斂的前一個折減系數(shù)取值的情況下，得到了如下的塑性應(yīng)變圖(圖17)。

圖17 更改灰?guī)r參數(shù)塑性應(yīng)變圖

該圖依然有和勘查結(jié)果接近的剪出口，和最大塑性應(yīng)變位置。

從上面結(jié)果可以看出，分層巖土參數(shù)加權(quán)平均法，對于計算滑坡的總體趨勢還是比較有效的，能夠近似的模擬出剪出口的位置。

6 結(jié)論

在結(jié)合實際滑坡例子的基礎(chǔ)上，利用CATIA三維設(shè)計軟件建立了地質(zhì)模型。使用分層巖土參數(shù)加權(quán)平均法對巖土參數(shù)做了處理，并對模型進(jìn)行了可導(dǎo)入ANSYS進(jìn)行有限元計算的分析處理。在ANSYS分析中，巖土體采用Drucker-Prager模型，對簡化后的滑坡模型進(jìn)行了有限元計算。計算的結(jié)果顯示與實際情況比較接近。因此可以得出以下結(jié)論：采用分層巖土加權(quán)平均法處理后的模型和巖土參數(shù)進(jìn)行有限元分析計算可以近似的反映出滑坡情況，推測滑坡發(fā)展趨勢，因此該方法是一個較有效的對CATIA軟件建立的地質(zhì)模型進(jìn)行有限元分析的一個方法。

[1] 苑齋明.基于CATIA平臺重力壩三維設(shè)計及抗滑穩(wěn)定分析[D].大連理工大學(xué),2015.

[2] 蘇小寧，王野，韓旭，等.基于CATIA的三維地質(zhì)建模及可視化應(yīng)用研究[J].人民長江，2015，(19)：98-104.

[3] 錢驊,喬世范,許文龍,等.水利水電三維地質(zhì)模型覆蓋層建模技術(shù)研究[J].巖土力學(xué)，2014，(7)：2103-2108.

[4] 錢驊,喬世范,許文龍,等.水利水電三維地質(zhì)模型覆蓋層建模技術(shù)研究[J].巖土力學(xué),2014,(7):2103-2108.

[5] 盧波,葛修潤,孔祥禮.有限元法、無單元法及自然單元法之比較研究[J].巖石力學(xué)與工程學(xué)報,2005,(5):780-786.

[6] 鄭穎人,趙尚毅.巖土工程極限分析有限元法及其應(yīng)用[J].土木工程學(xué)報,2005,38(1):91-98+104.

[7] Paul S.Heckbert, Miehael Garland. Multiresolution Modeling for Fast Rendering[J].Graphics Intersacu,1970:43-50.

[8] Gabriel Taubin,JarekRossignac.3D Geometry ComPression[J].Computer Graphics,1995,(22):13-20.

[9] 鄭穎人,趙尚毅.有限元強度折減法在土坡與巖坡中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報,2004,23(19):3381-3388.

[10]梅婷婷.復(fù)雜地下巖體應(yīng)力分析三維建模技術(shù)及工程應(yīng)用[D].武漢工程大學(xué),2012.

[11]Dawson，E.M.，Roth W.H. & Drescher A. Slope stabilityanalysis by strength reduction [J].Geotechnique，1999，49(6)：835-840.

[12]D.V.Griffiths and P.A.lane. Slope stability analysis byfinite elements[J].Geotechnique, 1999, 49(3): 387-403.

作者簡介：蒙鐸(1991- )，男，四川南部縣人，碩士，從事地質(zhì)災(zāi)害防治工作。E-mail：342781655@qq.com

A PREPROCESSING METHOD AND IT’S IMPLEMENTATION FOR THE GEOLOGICAL MODELING AND THE FINITE ELEMENT CALCULATION BASED ON CATIA SOFTWARE

MENG Duo1, CHEN Zhi-wen2, HUANG Jun1, ZHANG Chang3
(1.Faculty of Geosciences and Environmental Engineering,Southwest Jiao Tong University，Chengdu 611756,China;2.Faculty of Management Sciences,Chengdu University of Technology,Chengdu 610059,China;3.College of Environment and Civil Engineering,Chengdu University of Technology,Chengdu 610059,Chin)

Complex geological models established by CATIA software, it is difficult to directly used in finite element calculation, needs pretreatments on the established model. Proposes a pretreatment method: application of geotechnical parameters of the weighted average method, and connecting with the engineering geological analysis, finite element calculation principle and topology analysis method to do pretreatments of the complex geological models established by CAITA software. The pretreatment method can not only approximately reflect the geological conditions, but also meet the requirements of finite element analysis, which can guarantee the accuracy of the analysis results. Through the geological modeling and finite element calculation for a landslide in Sichuan, verified this method by an example, it has obtained a good effect.

geological modeling; finite element calculation; topology; pretreatments

1006-4362(2016)04-0091-07

2016-07-20改回日期： 2016-08-25

2013年四川省戰(zhàn)略性新興產(chǎn)業(yè)發(fā)展專項資金項目(SC2013510109139)

X141;P642.22

A