摘要 巖體工程所處地質(zhì)條件復(fù)雜多變，在巖體內(nèi)部極易出現(xiàn)水力耦合現(xiàn)象。針對(duì)此現(xiàn)象進(jìn)行的數(shù)值模擬，可以及時(shí)反映該含水巖體狀態(tài)及變化，是保證工程正常施工與運(yùn)營(yíng)的重要手段。主要分析了不同巖體滲流數(shù)值模擬方法，闡述了各種數(shù)值模擬方法優(yōu)缺點(diǎn)，提出未來(lái)學(xué)者應(yīng)著眼于將多種數(shù)值模擬方法相耦合，以便達(dá)到快速高效解決實(shí)際工程問(wèn)題的目的。

關(guān)鍵詞 水力耦合；數(shù)值模擬；巖體滲流

中圖分類號(hào) TU457 文獻(xiàn)標(biāo)志碼 A

0引言

隨著國(guó)民經(jīng)濟(jì)的快速發(fā)展，大量的地質(zhì)問(wèn)題出現(xiàn)在露天礦開(kāi)采、鐵路公路建設(shè)和水利水電等諸多工程領(lǐng)域中。而地質(zhì)巖體在經(jīng)歷數(shù)百年演化作用后，將會(huì)產(chǎn)生孔隙、節(jié)理和裂紋等多種地質(zhì)情況。在此基礎(chǔ)上，地下水由于應(yīng)力作用會(huì)使巖體在裂縫中產(chǎn)生位移，從而以滲流的形式對(duì)巖體造成變形破壞，同時(shí)也影響了巖體內(nèi)部應(yīng)力場(chǎng)分布。該現(xiàn)象可能對(duì)裂隙結(jié)構(gòu)層形成沖擊，也可能造成巖塊本身滲透特性發(fā)生改變。這種滲流與應(yīng)力互相作用并改變彼此的現(xiàn)象被稱為水力耦合。如不及時(shí)處理此類問(wèn)題，將極易產(chǎn)生大壩失穩(wěn)垮塌、邊坡滑坡、隧道施工過(guò)程中的巖溶塌陷涌水等事故，傳統(tǒng)單一的仿真技術(shù)已經(jīng)不能滿足現(xiàn)代巖體工程施工的需求，發(fā)展更智能、高效的水力耦合數(shù)值模擬方法是必然趨勢(shì)。

近年來(lái)，計(jì)算機(jī)技術(shù)的不斷發(fā)展為國(guó)內(nèi)外學(xué)者拓寬了研究思路，在探索水力耦合時(shí)，主要著眼于數(shù)字模型理論分析與數(shù)值模擬相結(jié)合的方法，并取得了不錯(cuò)的進(jìn)展。因此，本文擬借助國(guó)內(nèi)外學(xué)者目前研究數(shù)據(jù)，對(duì)巖體在水力耦合作用下滲流特性的數(shù)值模擬方法進(jìn)行總結(jié)，并對(duì)未來(lái)發(fā)展做出預(yù)測(cè)，為今后相關(guān)工程項(xiàng)目提供一定指導(dǎo)。

1巖體的滲透定量分析方法

1.1DDA法（非連續(xù)變形法）

DDA方法，是一種分析不連續(xù)介質(zhì)力學(xué)性質(zhì)的數(shù)值仿真方法。該方案在考慮巖塊運(yùn)動(dòng)復(fù)雜性的基礎(chǔ)上，把由結(jié)構(gòu)面分割而成的塊體作為主要分析單位，聯(lián)合動(dòng)力學(xué)與靜力學(xué)問(wèn)題，并利用最小勢(shì)能原理，將塊體內(nèi)部的碰撞和塊體自身的運(yùn)動(dòng)變化聯(lián)合在矩陣中求解。該方法最早由美籍華人石根華［1］提出，其優(yōu)點(diǎn)是能夠模擬節(jié)理裂隙巖體產(chǎn)生大變形和大位移的效果。學(xué)者們?cè)谶吰鹿こ桃约傲鞴恬詈戏矫娲罅坎扇〈朔椒ǎ瑥垏?guó)新等［2］在利用DDA法的基礎(chǔ)上，考慮地下水的流動(dòng)，以及滲流壓力與巖石變形的耦合作用，對(duì)日本某隧道的坍塌事件進(jìn)行了仿真分析，再現(xiàn)了該隧道破裂的過(guò)程；李志環(huán)［3］以中國(guó)遼寧省白石水庫(kù)上游及右岸邊坡為例，在考察了水力耦合和邊坡巖體蠕動(dòng)變化的基礎(chǔ)上，在變形量、變化形態(tài)及其變化過(guò)程等方面開(kāi)展了DDA模擬計(jì)算，其結(jié)果與現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)相比較，有較大的重疊度。

1.2有限元分析法

目前，國(guó)內(nèi)外大多數(shù)學(xué)者采用有限元分析法分析巖體地下水問(wèn)題，其中包括有限單元法和有限差分法等，主要應(yīng)用的軟件有MidasGTSNX、FLAC3D、COMSOL、ABAQUS、ANSYS等。楊波等［4］依托杭州紫之隧道，利用ABAQUS有限元軟件進(jìn)行數(shù)值模擬，主要研究了滲透系數(shù)的取值對(duì)于隧道涌水量以及降深的影響。模擬結(jié)果表明，圍巖入滲系數(shù)、滲透系數(shù)、隧道的洞距大小都會(huì)對(duì)隧道內(nèi)涌水量產(chǎn)生影響；李廷春等［5］應(yīng)用FLAC3D軟件模擬了海底隧道開(kāi)挖后圍巖力學(xué)特性，為今后類似隧道工程的穩(wěn)定性提供參考。

其中，有限單元法的求解思路為，將計(jì)算區(qū)域分割為有限多個(gè)彼此互不交叉的單元，從每一單元中選取相應(yīng)的控制軸線作為求解函數(shù)的插值運(yùn)算節(jié)點(diǎn)，將求解區(qū)域分割成有限多個(gè)彼此互不交叉的小單元，從每一單元中選取相應(yīng)的控制軸線作為求解變量的內(nèi)插點(diǎn)，將變量改寫為線性表達(dá)式，并運(yùn)用變分基本原理，以實(shí)現(xiàn)常微分方程中離散型的求解。由于有限單元法中對(duì)地下水流動(dòng)、溶質(zhì)的數(shù)值模擬的求解方法基本相同，因此在處理各種溶質(zhì)之間的耦合問(wèn)題時(shí)也較為簡(jiǎn)便。

而有限差分法則不同，其基本思路是用網(wǎng)格節(jié)點(diǎn)的函數(shù)值的差商代替進(jìn)行離散，建立以網(wǎng)格節(jié)點(diǎn)上的值為未知數(shù)的代數(shù)方程組。有限差分法較適用于一維以及二維流等問(wèn)題，求解流程比較簡(jiǎn)潔清晰；同時(shí)，在地下水水流數(shù)值模擬過(guò)程中，精度較高。

1.3無(wú)單元法

有限元法在解決大變形、裂紋擴(kuò)展等問(wèn)題時(shí)，多采用網(wǎng)格重構(gòu)技術(shù)，而這將會(huì)導(dǎo)致計(jì)算精度嚴(yán)重受損，還易出現(xiàn)求解不收斂等問(wèn)題，算法效率也很低下。而無(wú)單元法可以巧妙避開(kāi)此做法，能保證仿真精度與效率。而無(wú)單元法在計(jì)算過(guò)程中隨意產(chǎn)生和加密結(jié)點(diǎn)，不受單元網(wǎng)格的限制，由Belytschko于1981年提出的，后經(jīng)寇曉東等［6］的完善和發(fā)展，已經(jīng)在裂隙巖體水力耦合得到了較廣的應(yīng)用。李曉春等［7］采用了有限元法和無(wú)網(wǎng)格法耦合求解裂隙巖體滲流問(wèn)題，但也只局限于少數(shù)裂隙的求解，對(duì)分布大量密集裂隙的巖體滲流特性還需要展開(kāi)深入研究［8］。

1.4BP神經(jīng)網(wǎng)絡(luò)法

在模型優(yōu)化上方面，大多數(shù)學(xué)者采用BP神經(jīng)網(wǎng)絡(luò)法，原因是其具有全局最優(yōu)求解、普適性等優(yōu)勢(shì)，因此，該算法也被廣泛應(yīng)用在巖體水力耦合方面。徐鎮(zhèn)凱等［9］通過(guò)模擬庫(kù)水位、降雨的變化情況以及對(duì)大壩滲流的作用，提出滯后影響函數(shù)來(lái)模擬該過(guò)程，結(jié)果表明，用自適應(yīng)BP神經(jīng)網(wǎng)絡(luò)法具有收斂速度快、執(zhí)行效率高等優(yōu)點(diǎn)。

但是，傳統(tǒng)的BP神經(jīng)網(wǎng)絡(luò)算法仍然存在訓(xùn)練效果低下、容易進(jìn)入局部極值狀態(tài)等問(wèn)題，許多研究人員對(duì)此方法進(jìn)行了改進(jìn)。吳云星等［10］采用LM算法，將同結(jié)構(gòu)的LMBP神經(jīng)網(wǎng)絡(luò)和標(biāo)準(zhǔn)BP神經(jīng)網(wǎng)絡(luò)應(yīng)用于某土石壩滲流壓力的預(yù)測(cè)中。應(yīng)用結(jié)果表明，前者具有更高的收斂速度與精度，可廣泛應(yīng)用于滲流壓力分析和預(yù)測(cè)問(wèn)題中。張海龍［11］引入遺傳蟻群混合算法（GACO），通過(guò)優(yōu)化后的方法調(diào)整了標(biāo)準(zhǔn)BP神經(jīng)網(wǎng)路的權(quán)值和閾值結(jié)構(gòu)，對(duì)土石壩的線性滲流壓力進(jìn)行了監(jiān)控，將預(yù)測(cè)結(jié)果與經(jīng)典BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)結(jié)果進(jìn)行對(duì)比后發(fā)現(xiàn)，改進(jìn)后的模型明顯更勝一籌，具有較好的實(shí)用性。

1.5反演分析方法

反演分析是根據(jù)巖土體在實(shí)際工程荷載作用下監(jiān)測(cè)到的性狀變化，采用數(shù)值模擬方法對(duì)巖土體的力學(xué)特性和初始應(yīng)力條件進(jìn)行分析的方法。而巖體的滲流特征成分復(fù)雜、歷史悠久和不確定性影響因素較多，迄今為止，還沒(méi)有一種可以快速、經(jīng)濟(jì)、高效確定其滲流系數(shù)張量的手段，所以需要運(yùn)用所掌握的大量觀測(cè)數(shù)據(jù)，來(lái)反演推斷巖體中的滲透系數(shù)與邊界條件。徐強(qiáng)［12］在三維非穩(wěn)定滲流反演中使用改進(jìn)后的粒子群算法，以抽水蓄能電站地下廠房為例，對(duì)非穩(wěn)定滲流進(jìn)行反演分析，并對(duì)其長(zhǎng)期運(yùn)行安全性做出評(píng)價(jià)。藺雅嫻［13］以“渭源至武都建設(shè)項(xiàng)目WWSY2標(biāo)段沙灣隧道”為實(shí)際工程背景，結(jié)合位移反分析方法，對(duì)圍巖的應(yīng)力-滲流耦合蠕變多參數(shù)進(jìn)行反演，并用反演所得參數(shù)進(jìn)行ABAQUS有限元分析，結(jié)果表明：隧道圍巖在蠕變期處于一個(gè)穩(wěn)定的狀態(tài)。未來(lái)，學(xué)者們的研究重點(diǎn)應(yīng)轉(zhuǎn)移到巖土工程反演與人工智能的結(jié)合方面，這具有很強(qiáng)的現(xiàn)實(shí)意義和研究?jī)r(jià)值。

1.6邊界元法（BEM）

邊界元法是近年來(lái)新興的一種仿真方法，其只在定義域的邊界上劃分單元，利用滿足條件的方程式去逼近邊界條件，相較于傳統(tǒng)有限元法，存在劃分單元個(gè)數(shù)少等明顯優(yōu)勢(shì)。而在邊坡巖體滲流問(wèn)題中，邊界元法可使自由面滲流問(wèn)題等大為簡(jiǎn)化，其優(yōu)越性相當(dāng)明顯。李清新［14］通過(guò)直接邊界元方法，對(duì)同質(zhì)設(shè)計(jì)各向同性滲流問(wèn)題、正交設(shè)計(jì)各向異性滲流問(wèn)題、非飽和土壤滲流問(wèn)題、分層土壤的線性滲流問(wèn)題等進(jìn)行了研究，除完成以上4種線性滲流問(wèn)題的基本求解流程之外，每一類問(wèn)題還都編制了matlab方法，并通過(guò)算例加以證明。雷衛(wèi)東等［15］采用邊界元法處理非飽和土穩(wěn)態(tài)滲流問(wèn)題，應(yīng)用邊界元法求解非飽和土滲流水頭，提出應(yīng)用Matlab自帶的Fsolve函數(shù)求解邊界積分方程中的非線性方程組，計(jì)算結(jié)果具有良好的穩(wěn)定性。

大多數(shù)情況下，邊界元法由于其僅考慮了線性滲流的邊界上部分單元的特點(diǎn)，被廣泛應(yīng)用與確定線性滲流自由面問(wèn)題方面。

1.7離散元法

離散元法是一種在動(dòng)力穩(wěn)定方面廣為通用的方法。它是由分析離散單元的塊間接觸入手，找出其接觸的本構(gòu)關(guān)系，建立接觸的物理力學(xué)模型，并基于牛頓第二定律進(jìn)行仿真計(jì)算，可分析巖塊內(nèi)部的變形與應(yīng)力分布，反映巖塊之間的滑移、分離與傾翻等大位移。目前，大多數(shù)學(xué)者在利用離散元數(shù)值模擬軟件研究滲流力學(xué)特性時(shí)，采用3DEC軟件。張玉福［16］探討了離散元法在涌水量預(yù)測(cè)中應(yīng)用的可行性，并對(duì)3DEC軟件實(shí)施的二級(jí)設(shè)計(jì)，為相關(guān)項(xiàng)目提供思路借鑒；唐忠林等［17］以立方體節(jié)理巖體為研究對(duì)象，在3DEC平臺(tái)上完成了立方體單裂隙節(jié)理及節(jié)理矩陣模型的構(gòu)建，并獲得了按水壓方向依次遞減的立方體單裂隙水壓分布云圖和節(jié)理矩陣水壓分布云圖。

2結(jié)束語(yǔ)

本文詳細(xì)歸納總結(jié)了國(guó)內(nèi)外學(xué)者在研究水力耦合巖體時(shí)所使用的數(shù)值分析方法，包括DDA法、有限元法、無(wú)單元法、神經(jīng)網(wǎng)絡(luò)法、反演分析方法、邊界元法、離散元法等。隨著科技發(fā)展，更多復(fù)雜裂隙水力耦合問(wèn)題亟需解決。然而，目前大部分單一數(shù)值方法還依然存在計(jì)算速度低效、精度不夠以及可操作性不夠強(qiáng)等問(wèn)題，因此，學(xué)者在未來(lái)進(jìn)行研究時(shí)應(yīng)著重考慮將多種數(shù)值方法結(jié)合起來(lái)，如：有限元-離散元耦合法，以便達(dá)到快速解決實(shí)際工程問(wèn)題的目的。同時(shí)，還需進(jìn)一步展開(kāi)室內(nèi)實(shí)驗(yàn)和現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)分析，驗(yàn)證該仿真方法的準(zhǔn)確性。未來(lái)，水力耦合數(shù)值模擬方法將日趨完善與優(yōu)化，并在巖體工程中發(fā)揮重要作用。

參考文獻(xiàn)

［1］石根華.數(shù)值流形方法與非連續(xù)變形分析［M］.裴覺(jué)民譯.北京：清華大學(xué)出版社，1997.97-114.

［2］張國(guó)新，武曉峰.裂隙滲流對(duì)巖石邊坡穩(wěn)定的影響——滲流、變形耦合作用的DDA法［J］.巖石力學(xué)與工程學(xué)報(bào)，2003（8）：1269-1275.

［3］李志環(huán).基于DDA分析的白石水庫(kù)庫(kù)岸硬巖邊坡傾倒變形研究［J］.水利科學(xué)與寒區(qū)工程，2020，3（1）：76-78.

［4］楊波，胥東.基于三維有限元的杭州紫之隧道運(yùn)營(yíng)期滲流場(chǎng)變化預(yù)測(cè)［J］.山西建筑，2019，45（19）：117-119.

［5］李廷春，李術(shù)才，陳衛(wèi)忠，等.廈門海底隧道的流固耦合分析［J］.巖土工程學(xué)報(bào)，2004（3）：397-401.

［6］寇曉東，周維垣.應(yīng)用無(wú)單元法追蹤裂紋擴(kuò)展［J］.巖石力學(xué)與工程學(xué)報(bào)，2000（1）：18-23.

［7］李曉春，陳劍平，石丙飛，等.交叉裂隙滲流無(wú)網(wǎng)格法初探［J］.巖土力學(xué)，2007，28（S1）：371-374.

［8］李九紅，程玉民.一種新的無(wú)網(wǎng)格方法與有限元耦合法［J］.工程數(shù)學(xué)學(xué)報(bào)，2008（6）：1035-1043.

［9］徐鎮(zhèn)凱，蔡磊，魏博文，等.考慮滲流時(shí)滯及參數(shù)不確定性的土壩壩坡失穩(wěn)概率分析［J］.水力發(fā)電，2018，44（6）：62-66+70.

［10］吳云星，周貴寶，谷艷昌，等.基于LMBP神經(jīng)網(wǎng)絡(luò)的土石壩滲流壓力預(yù)測(cè)［J］.人民黃河，2017，39（8）：90-94+148.

［11］張海龍.改進(jìn)的土石壩GACO-BP滲流監(jiān)控模型研究［J］.陜西水利，2020（5）：19-21+25.

［12］徐強(qiáng).地下洞室非穩(wěn)定滲流場(chǎng)模擬及反演分析［D］.武漢：武漢大學(xué)，2019.

［13］藺雅嫻.基于IA-BP智能算法的應(yīng)力-滲流耦合蠕變多參數(shù)反演方法及工程應(yīng)用［D］.沈陽(yáng)：沈陽(yáng)工業(yè)大學(xué)，2019.

［14］李清新.基于邊界元法對(duì)穩(wěn)態(tài)滲流問(wèn)題的研究［D］.哈爾濱：哈爾濱工業(yè)大學(xué)，2014.

［15］雷衛(wèi)東，李清新，陳銳.非飽和土穩(wěn)態(tài)二維滲流的邊界元法研究［J］.水利學(xué)報(bào)，2015，46（S1）：95-100.

［16］張玉福.三維離散元法在巷道涌水量預(yù)測(cè)中的應(yīng)用［J］.金屬礦山，2013（12）：123-126.

［17］唐忠林，劉志義，曾維成.純滲流條件下節(jié)理巖體的滲流數(shù)值模擬研究［J］.交通世界，2020（32）：83-85.