基于逆向工程建模方法的危巖體穩(wěn)定性數(shù)值模擬分析*

高相波 李麗慧 廖小輝 王學(xué)良 陳子干

(①中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院頁巖氣與地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室，北京 100029，中國) (②中國科學(xué)院地球科學(xué)研究院，北京 100029，中國) (③中國科學(xué)院大學(xué)，北京 100049，中國) (④衢州學(xué)院，衢州 324004，中國) (⑤浙江神仙居旅游集團(tuán)有限公司，臺州 317300，中國)

0 引 言

自然界地質(zhì)體的穩(wěn)定性分析是復(fù)雜的地質(zhì)力學(xué)問題，地質(zhì)工作者需要在野外地質(zhì)調(diào)查、工程測量等工作的基礎(chǔ)上，建立相應(yīng)的地質(zhì)力學(xué)模型。并采用一定的理論與方法，分析、模擬、再現(xiàn)地質(zhì)體的災(zāi)害發(fā)生發(fā)展過程，以解決實(shí)際工程問題(Huang et al.，1999)。

巖質(zhì)邊坡危巖體失穩(wěn)的研究一直是學(xué)者們關(guān)注的重點(diǎn)問題。傳統(tǒng)的穩(wěn)定性評價(jià)方法包括定性分析方法(李靖等，1998)，如自然歷史分析法、工程類比法和圖解法；定量分析方法(鄭穎人等，2001；盧坤林等，2013；唐樂人，2014)，如極限平衡法、數(shù)值模擬方法等；以及非確定性分析方法(夏元友等，2002；唐慧群，2005；楊天鴻等，2011)，如可靠度分析法、聚類分析法、灰色系統(tǒng)分析法等。目前工程實(shí)踐中最常用的危巖體穩(wěn)定性分析方法是極限平衡法，通過計(jì)算潛在滑移破壞面上的抗滑力(矩)與滑動力(矩)之比，得到穩(wěn)定性系數(shù)來判斷其穩(wěn)定性(陳祖煜，2003)。并在二維極限平衡法的基礎(chǔ)上發(fā)展了三維極限平衡法(Huang et al.，2000；Chen et al.，2001；張均鋒等，2005；鄭宏，2007)。也有學(xué)者實(shí)現(xiàn)了極限平衡法的有限元計(jì)算(Griffiths et al.，1999；邵龍?zhí)兜龋?001；曾亞武等，2005；于斯瀅等，2013)。但是極限平衡法適用于幾何形狀規(guī)則和荷載作用簡單的邊坡穩(wěn)定性評價(jià)，這與巖質(zhì)邊坡變形破壞受結(jié)構(gòu)面控制的實(shí)際情況存在一定的差距，僅憑極限平衡法的計(jì)算結(jié)果難以揭示危巖體失穩(wěn)破壞機(jī)理。無從得知巖體內(nèi)部應(yīng)力分布狀態(tài)，因此得到的穩(wěn)定性評價(jià)結(jié)果也是不準(zhǔn)確的。

近年來，數(shù)值模擬分析作為一種重要的邊坡穩(wěn)定性評價(jià)方法，可以對邊坡內(nèi)部應(yīng)力場、變形場等的分布特征進(jìn)行定量化分析，在工程地質(zhì)領(lǐng)域中得到越來越廣泛的關(guān)注和應(yīng)用，解決了許多理論和實(shí)際問題。目前較為流行的數(shù)值模擬方法主要包括有限單元法(鄭穎人等，2002，2004)、離散單元法(Lu et al.，2014；石秋俠，2014；王穎等，2018)、有限差分法(胡卸文等，2019)、非連續(xù)變形分析法(裴覺民等，1993；Maclaughlin et al.，2001；吳建宏等，2003；鄔愛清等，2008)等。由美國Itasca公司開發(fā)的三維有限差分軟件FLAC3D以其在工程地質(zhì)問題分析計(jì)算方面快速而穩(wěn)定的求解優(yōu)勢，成為研究人員理想的數(shù)值模擬計(jì)算工具。然而FLAC3D軟件在復(fù)雜地質(zhì)模型的建立與不規(guī)則模型網(wǎng)格的劃分等前處理模塊中仍然力不從心，存在地質(zhì)模型過于簡化、網(wǎng)格劃分過程經(jīng)常出錯以及復(fù)雜建模工作量太大、耗時(shí)過長等缺陷。進(jìn)而影響了數(shù)值模擬結(jié)果的準(zhǔn)確性和可靠性。因此有學(xué)者通過ANSYS有限元程序完成復(fù)雜地質(zhì)體建模、網(wǎng)格劃分，再利用接口程序?qū)隖LAC3D，實(shí)現(xiàn)建模的直觀化、快速化和自動化(廖秋林等，2005)；基于AutoCAD平臺構(gòu)建三維可視化模型(鄭文棠等，2007)；將三維地質(zhì)建模與數(shù)值模擬分析進(jìn)行耦合，開發(fā)地質(zhì)建模與分析系統(tǒng)VisualGeo，簡化數(shù)值模擬的前處理工作(李明超等，2007)；基于Surfer平臺提取、轉(zhuǎn)換地表及巖層分界面三維地質(zhì)信息，生成模型數(shù)據(jù)文件(崔芳鵬等，2008)等。

除此之外，在高陡邊坡的工程地質(zhì)勘察中，研究人員往往受地形因素的影響難以開展實(shí)地踏勘工作，制約了邊坡巖體結(jié)構(gòu)特征等信息的獲取。甚至在潛在危巖體區(qū)域的地質(zhì)災(zāi)害調(diào)查工作還存在人身安全的問題。因此，有研究人員利用無人機(jī)和三維激光掃描儀(董秀軍等，2006；Fekete et al.，2013)等新技術(shù)代替調(diào)查人員進(jìn)行地質(zhì)結(jié)構(gòu)的高精度測量和地質(zhì)結(jié)構(gòu)信息獲取工作。然后利用逆向工程建模軟件將獲取的地形等高線和地層界面、節(jié)理裂隙面等點(diǎn)云數(shù)據(jù)進(jìn)行CAD曲面模型重建(徐文杰等，2008；鄧小龍等，2017)，可以實(shí)現(xiàn)三維地質(zhì)模型的構(gòu)建?？梢钥闯觯延械臄?shù)值建模方法基本都是在野外調(diào)查結(jié)果的基礎(chǔ)上，通過建模軟件構(gòu)建地質(zhì)體模型，再導(dǎo)入數(shù)值模擬軟件中計(jì)算分析巖體穩(wěn)定性，操作較為繁瑣，且過程中可能丟失地質(zhì)體本身的巖體結(jié)構(gòu)信息。本文基于三維激光掃描技術(shù)，快速便捷地獲取邊坡巖體結(jié)構(gòu)點(diǎn)云數(shù)據(jù)，在逆向工程建模軟件中構(gòu)建精細(xì)化數(shù)值模型，盡可能還原巖體本身地質(zhì)結(jié)構(gòu)特征，最后通過數(shù)值模擬軟件分析巖體穩(wěn)定性，并計(jì)算其失穩(wěn)概率。

以浙江神仙居景區(qū)地質(zhì)災(zāi)害調(diào)查為依托，選取景區(qū)標(biāo)志性景點(diǎn)神龍瀑為研究對象，在野外地質(zhì)調(diào)查的基礎(chǔ)上，基于逆向工程建模方法利用FLAC3D有限差分?jǐn)?shù)值計(jì)算軟件進(jìn)行危巖體穩(wěn)定性分析。通過三維激光掃描技術(shù)獲取神龍瀑邊坡的巖體結(jié)構(gòu)幾何信息，利用逆向工程軟件Geomagic Studio對點(diǎn)云數(shù)據(jù)進(jìn)行處理，重構(gòu)坡體的曲面模型；再利用Hypermesh進(jìn)行有限單元網(wǎng)格劃分，將曲面模型實(shí)體化和網(wǎng)格化；最后基于FLAC3D軟件進(jìn)行數(shù)值模擬分析，并根據(jù)結(jié)構(gòu)面塑性破壞單元數(shù)目計(jì)算危巖體失穩(wěn)概率。為當(dāng)?shù)芈糜喂芾聿块T的災(zāi)害防治工作提供理論指導(dǎo)和依據(jù)，避免人員傷亡和財(cái)產(chǎn)損失。

圖 1 神龍瀑邊坡Fig. 1 Shenlong waterfall slope

1 研究區(qū)工程地質(zhì)條件

本文研究區(qū)域神龍瀑所在景區(qū)屬于中、低山區(qū)，相對高差200～600im，地形切割強(qiáng)烈，最大切割深度達(dá)800im，地勢陡峭，局部坡度達(dá)到80°以上。常見石峰、石嶺、石壁、懸崖等地貌。神龍瀑兩側(cè)為高陡邊坡，中間為溝谷。東側(cè)山體坡角50°，坡向230°，坡長20im；西側(cè)山體邊坡反傾坡角80°，坡向50°，坡長26.5im，坡寬25im。研究區(qū)主體巖性為流紋巖、流紋質(zhì)晶玻屑熔結(jié)凝灰?guī)r，局部夾沉凝灰?guī)r。

神龍瀑兩側(cè)邊坡節(jié)理裂隙十分發(fā)育，巖體較為破碎，呈碎裂狀結(jié)構(gòu)。受多組結(jié)構(gòu)面切割影響，已形成多處危巖體，并頻有塊石墜落。溝谷中堆積有數(shù)百塊歷史崩塌產(chǎn)生的塊石，其中數(shù)十塊直徑在1im以上，最大塊石直徑達(dá)10im。圖 2給出了研究區(qū)節(jié)理走向玫瑰花圖，從圖中可知北東東向節(jié)理最發(fā)育，與區(qū)域內(nèi)斷層的走向50°基本一致，說明節(jié)理發(fā)育主要受斷層構(gòu)造作用控制。

圖 2 研究區(qū)節(jié)理走向玫瑰花圖Fig. 2 Rose diagram of joint strikes in study area

神龍瀑危巖區(qū)有地表水匯集且水量較大，兩側(cè)山體表面潮濕，局部有細(xì)小水流滲出。地下水類型為基巖裂隙水，補(bǔ)給來源主要為大氣降水。坡體表面風(fēng)化作用強(qiáng)烈，為中風(fēng)化-強(qiáng)風(fēng)化，局部巖塊用手即可掰碎。坡面上部長有直徑5～10icm的樹木22棵，樹木沿裂縫和節(jié)理生長，根劈作用明顯。如圖 3所示神龍瀑坡面植被發(fā)育情況。

圖 3 神龍瀑邊坡植被發(fā)育情況Fig. 3 Vegetation development on Shenlong waterfall slope

2 危巖體形成機(jī)制

神龍瀑邊坡坡向55°，坡面反傾，傾角達(dá)80°，可見多處拉張裂縫，在坡面臨空方向形成多處潛在危巖體。且頻有掉塊和落石墜于坡底。

據(jù)野外調(diào)查發(fā)現(xiàn)，神龍瀑邊坡兩側(cè)出露寬約1im的斷層破碎帶，走向50°。受斷層構(gòu)造作用影響，控制該區(qū)巖體穩(wěn)定性的主控結(jié)構(gòu)面走向大致為北東東向(參見圖 2研究區(qū)節(jié)理走向玫瑰花圖)。經(jīng)過野外結(jié)構(gòu)面統(tǒng)計(jì)與測量，神龍瀑邊坡巖體主要受3組結(jié)構(gòu)面切割，包括一組層理和兩組節(jié)理，巖體結(jié)構(gòu)呈碎裂狀。其中層理面產(chǎn)狀225°∠9°，可見范圍內(nèi)跡長10im，間距0.2im；節(jié)理1產(chǎn)狀340°∠84°，可見范圍內(nèi)跡長10im，間距0.5im，節(jié)理面平直，張開度0.003im，局部有泥質(zhì)填充；節(jié)理2產(chǎn)狀210°∠78°，局部較為發(fā)育，可見范圍內(nèi)跡長5im，間距0.3im，節(jié)理面較粗糙，張開度0.005im，局部有泥質(zhì)和碎石塊填充。根據(jù)巖體穩(wěn)定性受巖體結(jié)構(gòu)控制理論，神龍瀑邊坡危巖體穩(wěn)定性受上述3組結(jié)構(gòu)面的控制。

在降雨、地震、風(fēng)化等自然營力作用下，節(jié)理裂隙不斷擴(kuò)展，巖體強(qiáng)度逐漸弱化，極有可能誘發(fā)危巖體失穩(wěn)。降雨、地震和長期風(fēng)化作用是危巖體發(fā)生崩塌的主要誘因，危巖體破壞運(yùn)動方式主要表現(xiàn)為滑移-墜落。

3 危巖體穩(wěn)定性分析

有限差分?jǐn)?shù)值計(jì)算軟件FLAC3D(Fast Lagrangian Analysis of Continua in 3 Dimensions)即連續(xù)介質(zhì)快速拉格朗日分析方法，基于顯式差分法對運(yùn)動方程和動力方程進(jìn)行求解，且程序內(nèi)置命令可以設(shè)置具有一定接觸剛度的接觸面以模擬斷層和節(jié)理等結(jié)構(gòu)面，能夠快速而穩(wěn)定地求解連續(xù)介質(zhì)大變形、大位移問題。本文借助數(shù)值計(jì)算軟件FLAC3D對神龍瀑邊坡危巖體進(jìn)行穩(wěn)定性分析，為當(dāng)?shù)鼐皡^(qū)地質(zhì)災(zāi)害防治提供參考。

3.1 逆向工程建模方法

為了得到更為精確的數(shù)值計(jì)算結(jié)果，軟件本身前處理模塊的建模工具顯然難以勝任地質(zhì)體模型的準(zhǔn)確還原，因此采用逆向工程建模方法對神龍瀑邊坡的巖體結(jié)構(gòu)特征進(jìn)行精細(xì)刻畫。

首先利用三維激光掃描技術(shù)對神龍瀑邊坡坡面進(jìn)行高精度掃描。采用Leica scanstation C10三維激光掃描儀，該儀器單次測量點(diǎn)位精度為6imm，單次測量距離精度為4imm，標(biāo)靶掃描精度為2imm，模型表面精度為2imm。將獲得的點(diǎn)云數(shù)據(jù)在逆向建模軟件Geomagic Studio中進(jìn)行三維重建，經(jīng)過網(wǎng)格清理、模型封裝等處理手段，將邊坡發(fā)育的植被部分點(diǎn)云進(jìn)行降噪、清理，構(gòu)造高質(zhì)量NURBS曲面片模型(圖 4)?？梢娭亟Ｐ突颈Ａ袅诉吰聨r體結(jié)構(gòu)特征信息。

圖 4 Geomagic Studio重建曲面片模型圖Fig. 4 Reconstructive surface model by Geomagic Studio

然后將重建的坡面曲面片模型進(jìn)行四周延伸并三維實(shí)體化，構(gòu)建神龍瀑邊坡三維模型(圖 5)。此時(shí)封裝完成的模型由于數(shù)據(jù)格式的限制尚無法直接導(dǎo)入FLAC3D軟件中進(jìn)行數(shù)值計(jì)算，且不規(guī)則坡面形態(tài)的網(wǎng)格劃分過程也十分困難。

圖 5 神龍瀑邊坡三維重建模型圖Fig. 5 3D reconstructive model of Shenlong waterfall slope

圖 6 Hypermesh有限單元網(wǎng)格模型圖Fig. 6 Finite element mesh model by Hypermesh

因此借助有限單元網(wǎng)格劃分軟件Hypermesh進(jìn)行處理。根據(jù)現(xiàn)場結(jié)構(gòu)面調(diào)查統(tǒng)計(jì)結(jié)果，將封裝模型分為神龍瀑邊坡、斷層破碎帶和節(jié)理面3類。為盡可能準(zhǔn)確刻畫模型特征，對重點(diǎn)考察區(qū)域的網(wǎng)格進(jìn)行細(xì)化，以提高數(shù)值計(jì)算精度。本次建模設(shè)置巖體有限單元網(wǎng)格劃分的尺寸為0.2im，斷層和節(jié)理單元網(wǎng)格劃分的尺寸為0.1im。該模型共包括節(jié)點(diǎn)數(shù)118i588個，單元數(shù)645i805個。

將Hypermesh所生成的模型節(jié)點(diǎn)和單元信息文件導(dǎo)入FLAC3D軟件中，從而建立數(shù)值計(jì)算模型(圖 7)。在FLAC3D中設(shè)置邊界條件，模型底部約束豎向位移，四周邊界約束水平向位移，模型上表面為自由邊界。

表 1 模型分組及其力學(xué)參數(shù)取值Table1 Numerical model groups and its parameters

模型部位容重/kN·m-3黏聚力/MPa內(nèi)摩擦角/(°)抗拉強(qiáng)度/MPa體積模量/GPa剪切模量/GPaknks神龍瀑邊坡25370420.822.7斷層破碎帶251.3702.51011節(jié)理面0.12506.9e96.9e10

圖 7 FLAC3D 數(shù)值計(jì)算模型圖Fig. 7 Numerical model by FLAC3D

3.2 力學(xué)參數(shù)的選取

本次計(jì)算中力學(xué)參數(shù)主要依據(jù)經(jīng)驗(yàn)值，并根據(jù)模型計(jì)算結(jié)果做一定修正。對結(jié)構(gòu)面的剛度系數(shù)采用反演分析方法確定。計(jì)算模型的本構(gòu)關(guān)系選擇Mohr-Coulomb彈塑性模型。數(shù)值計(jì)算中使用體積模量K和剪切模量G來代替楊氏模量E和泊松比ν。模型介質(zhì)分組及其力學(xué)參數(shù)取值如表 1所示。

3.3 數(shù)值模擬分析

圖 8給出了神龍瀑邊坡在主控結(jié)構(gòu)面切割作用下的危巖體分區(qū)。圖 9為對應(yīng)危巖體的塑性破壞區(qū)。從塑性破壞區(qū)可以明顯看出，神龍瀑邊坡潛在破壞區(qū)域主要分布在上部危巖體。包括一處明顯的受拉破壞危巖體和局部剪切破壞的小型掉塊和落石。

圖 8 神龍瀑邊坡危巖體分區(qū)圖Fig. 8 Unstable rock mass zones of Shenlong waterfall slope

圖 9 危巖體塑性破壞區(qū)Fig. 9 Plastic damage areas of unstable rock mass

圖 10 受拉破壞危巖體主應(yīng)力分布剖面圖Fig. 10 Profile of principal stress in tensile damaged unstable rock massa. 最小主應(yīng)力；b. 最大主應(yīng)力

為了進(jìn)一步對邊坡巖體內(nèi)部應(yīng)力狀態(tài)進(jìn)行分析，在FLAC3D中輸出受拉破壞危巖體的主應(yīng)力分布剖面圖(圖 10)。由圖 10a最小主應(yīng)力分布剖面圖中應(yīng)力分布狀態(tài)可知，在自然狀態(tài)下，神龍瀑邊坡巖體內(nèi)部的最小主應(yīng)力為拉應(yīng)力，最大值超過1.74iMPa。且坡體本身反傾，在節(jié)理裂隙的切割作用下形成臨空危巖體，后緣裂隙貫通，危巖體出現(xiàn)拉應(yīng)力集中區(qū)。神龍瀑坡體除上部巖體分布拉應(yīng)力集中區(qū)外，下部巖體最大主應(yīng)力主要為壓應(yīng)力(圖 10b)，坡腳位置和節(jié)理裂隙接觸處出現(xiàn)壓應(yīng)力集中區(qū)，最大值在25iMPa左右。

圖 11 受拉破壞危巖體垂向應(yīng)力分布剖面圖Fig. 11 Profile of vertical stress in tensile damaged unstable rock mass

圖 11給出了受拉破壞危巖體垂向應(yīng)力分布剖面圖。由圖分析可知，坡體上部危巖體在層理和后緣裂隙切割作用下產(chǎn)生拉應(yīng)力集中區(qū)，拉應(yīng)力值超過1.6iMPa。在降雨、風(fēng)化、樹木根劈作用等影響下，節(jié)理裂隙進(jìn)一步發(fā)育擴(kuò)展，下部巖體在壓應(yīng)力作用下破碎失穩(wěn)，進(jìn)而可能導(dǎo)致上部危巖體發(fā)生滑移-拉裂破壞。

圖 12給出了自然狀態(tài)下受拉破壞危巖體位移矢量分布圖，神龍瀑邊坡整體位移主要分布在受拉破壞危巖體，該危巖體位移矢量向坡體臨空方向發(fā)展，極有可能產(chǎn)生危巖體失穩(wěn)，威脅沿路游客的人身安全與景區(qū)設(shè)施。

圖 12 受拉破壞危巖體位移矢量分布圖Fig. 12 Displacement vector in tensile damaged unstable rock mass

3.4 基于FISH語言的危巖體失穩(wěn)概率計(jì)算

在野外現(xiàn)場調(diào)查結(jié)果的基礎(chǔ)上，結(jié)合上文對神龍瀑邊坡危巖體形成機(jī)制的工程地質(zhì)分析結(jié)果，危巖體主要受3組結(jié)構(gòu)面切割影響。其自身發(fā)育的節(jié)理裂隙等結(jié)構(gòu)面控制了危巖體穩(wěn)定性和失穩(wěn)破壞方式。因此，考慮計(jì)算危巖體結(jié)構(gòu)面的破壞概率作為危巖體失穩(wěn)概率，評價(jià)危巖體穩(wěn)定性。

在FLAC3D數(shù)值分析結(jié)果的基礎(chǔ)上，利用FISH語言進(jìn)行二次開發(fā)，利用循環(huán)語句輸出神龍瀑邊坡數(shù)值模型中各個危巖體接觸面單元的應(yīng)力狀態(tài)，拉應(yīng)力值超過抗拉強(qiáng)度的單元定義為拉張破壞單元，剪應(yīng)力值超過抗剪強(qiáng)度的單元定義為剪切破壞單元。然后統(tǒng)計(jì)各個接觸面的塑性破壞單元數(shù)目，進(jìn)而計(jì)算出塑性破壞單元數(shù)目在接觸面單元總數(shù)中的占比作為相應(yīng)危巖體的破壞概率。

經(jīng)計(jì)算，出現(xiàn)較多塑性破壞單元的部位在受拉破壞危巖體的后緣裂隙面和層理面，兩處接觸面破壞概率均達(dá)到20.4%。結(jié)合上文危巖體穩(wěn)定性評價(jià)結(jié)果和FLAC3D的數(shù)值計(jì)算結(jié)果，可以預(yù)測神龍瀑邊坡的潛在破壞區(qū)域極有可能發(fā)生在受拉破壞危巖體，應(yīng)及時(shí)采取防治措施，避免景區(qū)游客人員傷亡和財(cái)產(chǎn)損失。

4 結(jié) 論

(1)經(jīng)過現(xiàn)場地質(zhì)調(diào)查，神龍瀑邊坡巖體在節(jié)理裂隙的切割作用下呈碎裂狀結(jié)構(gòu)，在外力作用下極易脫離母巖，產(chǎn)生掉塊。且由于長期風(fēng)化作用，節(jié)理裂隙進(jìn)一步擴(kuò)展，巖體強(qiáng)度不斷弱化，可能誘發(fā)危巖體失穩(wěn)，其破壞運(yùn)動方式主要表現(xiàn)為滑移-墜落。

(2)FLAC3D數(shù)值分析結(jié)果表明，神龍瀑坡體受節(jié)理裂隙的切割作用形成多處臨空危巖體。其中受拉破壞危巖體后緣裂隙位置產(chǎn)生拉應(yīng)力集中，其最大拉應(yīng)力值(1.7iMPa)接近巖石抗拉強(qiáng)度(2.5iMPa)，且該危巖體存在明顯的位移，極有可能發(fā)生滑移-墜落，威脅景區(qū)內(nèi)游客人身安全，破壞景區(qū)公共設(shè)施。

(3)基于FLAC3D進(jìn)行FISH語言二次開發(fā)，根據(jù)結(jié)構(gòu)面塑性破壞狀態(tài)，計(jì)算出危巖體的失穩(wěn)概率。結(jié)果表明，出現(xiàn)較多塑性破壞單元的部位在受拉破壞危巖體的后緣裂隙面和層理面，其破壞概率均達(dá)到20%以上。結(jié)合野外地質(zhì)調(diào)查結(jié)果和FLAC3D的數(shù)值模擬分析，可以預(yù)測該危巖體為潛在破壞區(qū)域，建議景區(qū)管理部門立即對神龍瀑邊坡潛在危巖體進(jìn)行清理或加固。