基于離散-連續(xù)耦合方法的隧道壓力拱特性研究

高 峰，譚緒凱*，陳曉宇，丁其樂，李 星

(1.重慶交通大學 土木工程學院，重慶 400074; 2.中機中聯(lián)工程有限公司，重慶 400039;3.林同棪國際工程咨詢(中國)有限公司，重慶 401121)

1 引 言

隧道壓力拱理論是隧道深淺埋界限劃分的基礎,對隧道工程計算及設計具有重要作用。自Kovari[1]首次提出隧道開挖過程中的成拱效應以來,相關學者就隧道壓力拱問題展開了大量的研究。在壓力拱的定義上,目前普遍認為隧道壓力拱的形成為力的傳遞路徑的改變所致,認為壓力拱是指洞室開挖后的卸荷效應破壞了初始應力場的平衡,導致圍巖不均勻變形和應力重分布,圍巖體系的荷載傳遞路線發(fā)生偏移,主應力大小和方向隨著荷載釋放逐漸發(fā)生改變,最終在洞室周邊形成類似拱形的圍巖保護圈[2-6]。

目前,針對隧道壓力拱問題的研究主要集中在壓力拱成拱機理、壓力拱演化規(guī)律及壓力拱范圍等方面[7-9],其研究手段主要為模型試驗方法和數(shù)值模擬方法,而模型試驗手段易受研究成本、邊界效應、模型尺寸、試驗周期及人為因素等影響,數(shù)值模擬手段成為最為普遍的研究方法。目前,常用的數(shù)值方法主要分為連續(xù)介質數(shù)值方法[10-12](如有限元法和有限差分法等)和非連續(xù)介質數(shù)值方法(如離散元法等)。從研究成果發(fā)現(xiàn),連續(xù)介質方法無法從微觀角度直觀分析和觀察隧道壓力拱的形成機理及發(fā)展特性。

隨著離散元理論的發(fā)展,其廣泛應用于巖土工程的微觀力學行為研究,在壓力拱研究方面也有初步嘗試。如扈世民等[13]采用顆粒流程序PFC2D研究了大斷面黃土隧道的塌落拱和壓力拱的發(fā)展過程,直觀得出了黃土隧道破壞演化過程,并從應力場角度推出了壓力拱的外邊界。但是地下工程為半無限體,進行研究時需要選取較大的模型降低人工邊界對計算結果的影響,此時采用離散元程序計算時需要成千上萬的顆粒,需要大量的計算機存儲空間和計算時間,對于計算機水平要求太高。

耦合算法自20世紀80年代初期由Felippa等[14]提出以來,得到了長足的發(fā)展。有限元法能很好地反應連續(xù)介質力學性能,且對計算機存儲空間和計算時間要求較低；離散元法能很好地模擬非連續(xù)介質力學性能。因此,連續(xù)-離散耦合模型能充分發(fā)揮兩種方法的優(yōu)勢,且能直觀反映隧道開挖后的圍巖受力發(fā)展機理。目前,其已成為國內(nèi)外研究的一個重點方向。嚴瓊等[15]分別采用FLAC和PFC建立了連續(xù)模型和離散模型,從連續(xù)-離散模型的數(shù)據(jù)一致性實現(xiàn)耦合,研究了邊坡工程的穩(wěn)定性；徐國文等[16]基于離散元-有限差分耦合算法,研究了層狀軟巖隧道開挖后的圍巖破壞機理；金煒楓等[17]基于改進后的離散-連續(xù)耦合模型對地下結構地震過程中的破壞過程進行了模擬。目前,隧道及地下工程的離散-連續(xù)耦合模型研究已經(jīng)取得了一些進展,但其準確性及適用性依然值得進一步研究。

離散-連續(xù)耦合的關鍵在于交界處的模擬方法及模擬效果。周健等[18,19]通過保證外圍連續(xù)單元和內(nèi)部離散顆粒在耦合過渡區(qū)的力或速度的一致，實現(xiàn)了離散-連續(xù)模型的耦合。張銳等[20]通過編程,基于離散-連續(xù)模型耦合過渡區(qū)的動量傳遞的連續(xù)性及空間-時間的雙耦合,實現(xiàn)了離散-連續(xù)模型的耦合。張鐸等[21]從兩個方面實現(xiàn)了離散-連續(xù)模型的耦合,一方面通過速度實現(xiàn)連續(xù)單元向離散顆粒的過渡(耦合過渡區(qū)連續(xù)單元的節(jié)點速度通過形函數(shù)分配到離散顆粒),另一方面通過力的傳遞實現(xiàn)離散顆粒向連續(xù)單元的過渡(耦合過渡區(qū)離散顆粒的平均應力施加在連續(xù)單元上)。

本文采用二維顆粒流程序PFC2D,基于離散-連續(xù)耦合方法對隧道施工的壓力拱進行了研究。在距離隧道洞室較遠區(qū)域，通過連續(xù)模型模擬總體的受力及位移；在洞室周邊圍巖，采用離散元模型研究了隧道施工過程中周邊圍巖力的傳遞路徑的變化,進而直觀展現(xiàn)了隧道開挖后壓力拱的成拱過程及成拱范圍。

2 離散-連續(xù)耦合計算方法確定

2.1 離散-連續(xù)耦合程序的實現(xiàn)思路

基于現(xiàn)有離散-連續(xù)耦合思路,本文提出以下優(yōu)化思路。離散-連續(xù)耦合模型的關鍵是離散元模型與連續(xù)介質模型之間的邊界處理和過渡方式,離散顆粒和有限單元之間存在接觸,兩者間的接觸力可通過接觸點的位移進行計算。該接觸力可通過荷載形式施加在離散顆粒上,也可以外荷載形式施加在有限單元上,進而實現(xiàn)交界面上相互作用的連續(xù)性?？紤]到有限單元的受力特點,本文對于施加在有限單元上的荷載進行一定處理,取一定范圍離散顆粒的平均應力?；谏鲜鲴詈纤悸?通過力的傳遞,編寫連續(xù)單元算法和耦合算法嵌入到PFC2D程序中,實現(xiàn)耦合計算。

2.2 平均應力計算方法

在離散-連續(xù)模型交界面上，一側是有限單元，一側是離散單元,如圖1所示。

圖1 顆粒與單元接觸示意圖

Fig.1 Contact geometry between ball and continuum element

目前，顆粒與連續(xù)單元間的力及重疊厚度求解方法已有成熟的研究成果。顆粒與連續(xù)單元是否接觸通過顆粒圓心到臨近單元邊界距離與顆粒半徑之間的大小來確定。圖1(b)為顆粒與連續(xù)單元之間接觸的一種特例,顆粒接觸到兩個單元的交點處,該情況通過顆粒圓心與交點間的距離進行接觸與否的判定,接觸力的方向為結點2到顆粒圓心的方向。為了避免此種情況計算多次接觸力,規(guī)定此種情況下,顆粒只與邊界的起點相接觸,而與邊界的終點不算接觸,即圖1只與2-3段接觸,與1-2段不算接觸。接觸力只要2-3段與顆粒計算,方向為結點2指向顆粒圓心。

考慮到連續(xù)單元的應力是單元面積內(nèi)內(nèi)力的均值,因此,接觸單元的應力應是耦合過渡區(qū)域一定范圍內(nèi)離散顆粒的接觸力的平均值。一個體積為V的區(qū)域的平均應力可由式(1)求得[22]，

(1)

一定區(qū)域的孔隙比為

n=1-Aball/Aregion

(2)

(3)

2.3 平均應力計算區(qū)域面積的確定

耦合過渡區(qū)域的合理區(qū)域的選定對力的傳遞準確性至關重要。本文對比了如下3種方案。

(1) 按照過渡區(qū)域連續(xù)模型單元大小確定區(qū)域面積,遍歷該面積內(nèi)的顆粒,計算其平均應力,孔隙率近似采用初始孔隙率。該方案不能保證該區(qū)域內(nèi)顆粒與單元的質量相等。

(2) 在接觸單元的相鄰位置放置測量圓,其直徑按照該單元的最大邊長取值,由測量圓獲取該區(qū)域的平均應力。該方案可保證孔隙率的準確性。

(3) 方案3近似于方案2,但其主要通過測量圓獲取該區(qū)域的平均應變,進而通過本構方程得出平均應力。

為了選取最佳區(qū)域確定方法,對上述3種方案與全連續(xù)模型進行對比計算。平面模型計算范圍取10 m×10 m,在模型正中間區(qū)域的2 m×2 m范圍建立離散元模型,如圖2所示。

圖2 耦合計算模型

Fig.2 Model of coupling calculation

計算采用contact-bond模型,目前微觀參數(shù)主要通過試湊方法獲得,由于本文主要驗證程序的有效性,所以參數(shù)設定思路是給定確定的微觀參數(shù),然后利用FISHTank的模擬室內(nèi)壓縮試驗的程序計算彈性模量和泊松比。給定的微觀計算參數(shù)列入表1。

表1 耦合模型微觀參數(shù)

Tab.1 Microcosmic parameter of coupling calculation

參數(shù)參數(shù)值參數(shù)參數(shù)值最小半徑/m0.01切向粘結力/N1e6粒徑比1.66接觸面摩擦系數(shù)0.5法向剛度/N·m-11e8顆粒密度/kg·m-32200切向剛度/N·m-16.6e7初始孔隙率0.16法向粘結力/N1e6

由FISHTank程序模擬無側限單軸壓縮試驗,得到模擬壓縮曲線,如圖3所示?？梢钥闯?顆粒模型的宏觀參數(shù)為彈性模量Es=44 MPa,泊松比υ=0.14,平均密度ρ=1848 kg/m3。

圖3 應力-應變曲線

Fig.3 Stress -strain curve

模型計算荷載選取自重荷載,得到重力作用下體系的應力場和位移場,并與FLAC計算結果進行對比。各方案計算結果相對于FLAC計算結果的相對誤差列入表2。

表2 3種方案相對誤差最大值對比

Tab.2 Maximum relative error comparison of three kinds of solution

對比方案相對誤差/%豎向位移豎向應力方案113.4519.9方案212.6516.6方案331.5059.6

由表2可知,方案3誤差最大,這是由于在計算應力時,顆粒間隙的應力為0,而在顆粒間隙間的速度不為0。利用測量圓的平均應變計算得到的應力值較小,與對應的連續(xù)單元的應力相差很大,導致誤差較大。方案1和方案2的誤差幾乎相當,方案2的誤差略小。故最終通過方案2確定平均應力計算區(qū)域的面積。

由表2可知,耦合軟件計算結果與FLAC2D程序計算結果存在一定的誤差,將兩者沉降量和應力值的相對誤差整理得出誤差等值線，如圖4所示。可以看出,耦合程序的沉降結果普遍大于FLAC2D程序計算結果。其中,空單元附近誤差值最大,遠離空單元的位置誤差相對較小。由于保證了總體圍巖的受力特性,并且通過調(diào)整內(nèi)部離散元顆粒的微觀參數(shù)保證內(nèi)部離散顆粒計算精度,因此,耦合程序的交界面處的誤差對后續(xù)研究影響較小。

2.4 耦合程序遍歷顆粒方法的確定

在計算連續(xù)單元的交界面上接觸力時,耦合程序需要遍歷所有的顆粒,確定與連續(xù)單元存在接觸的顆粒,通過其接觸點的力求解邊界平均應力。程序如果每次都需遍歷所有的離散顆粒,將耗費大量的計算時間。本文考慮到隧道開挖過程中,遠離隧道開挖區(qū)域的圍巖穩(wěn)定性較好,其內(nèi)部顆粒一般不會竄跑到離散-連續(xù)模型的交界面。因此,提出遍歷的離散顆粒只包含模型交界面向內(nèi)的三層顆粒,如圖5所示。

圖4 相對誤差等值線

Fig.4 Contour map of relative error

3 程序執(zhí)行過程

離散-連續(xù)耦合模型計算流程如圖6所示,計算流程中的主要程序函數(shù)如下。

(1) FP_startup函數(shù)。 通過該函數(shù)程序實現(xiàn)連續(xù)單元程序擬FLAC2D和離散元程序PFC2D的初始化設置。

(2) FP_flacstatic函數(shù)。 在離散元程序PFC2D計算完成后運行該函數(shù),提取PFC2D中顆粒接觸應力,并將應力結果作為邊界條件施加到連續(xù)單元上,通過擬FLAC2D程序進行motion計算。進而進行連續(xù)單元程序的計算。

(3) FP_flacdyn函數(shù)。 該函數(shù)與FP_flacstatic函數(shù)屬于一個計算過程,主要負責擬FLAC2D中的動力計算模塊。

圖5 耦合算法中遍歷顆粒示意

Fig.5 Sketch of ball needed to ergodic in coupling algorithm

(4) FP_scanball函數(shù)。該函數(shù)主要是計算連續(xù)單元的交界面(本文稱之為段)上的接觸力，并將其施加到相應的顆粒上,其計算方法是通過連續(xù)單元交界面節(jié)點位移計算出相應位置的離散顆粒的等效力,將該等效力施加到相應離散顆粒上。進而進行離散單元程序的計算。

離散-連續(xù)耦合模型計算中,時間步設置對程序正常運行影響較大,主要根據(jù)以下幾點進行設置。① 初步應滿足同時小于FLAC和PFC的臨界時間步； ② 隧道開挖模型中假定顆粒與接觸面之間不承受拉力,如果時間步設置過大,會由于重力加速度的作用，單個顆粒的速度過大,導致隧道頂部顆粒與接觸面出現(xiàn)分離,因此,需要通過試算調(diào)整合理的時間步,直到顆粒與接觸面不出現(xiàn)分離現(xiàn)象為止。

4 壓力拱的形成過程及原理分析

4.1 模型及參數(shù)

通過離散-連續(xù)耦合程序模擬計算了隧道開挖過程。模型范圍取60 m×60 m,隧道取半徑為 2 m 的圓形隧道,隧道埋深取18 m。隧道區(qū)域及隧道附近的20 m×20 m區(qū)域建立離散模型,其他區(qū)域建立連續(xù)模型。計算模型如圖7所示。通過FISHTank程序模擬無側限壓縮試驗得到離散顆粒的宏觀壓縮模量為E=4.0e8 Pa,泊松比為 0.14。其宏觀密度為1848 kg/m3。離散模型計算參數(shù)列入表3。

圖6 離散-連續(xù)耦合計算流程

Fig.6 Flowchart of discontinuum-continuum coupling calculation

圖7 數(shù)值模型

Fig.7 Numerical model

表3 離散顆粒的微觀參數(shù)

Tab.3 Microcosmic parameter of discrete particle

參數(shù)參數(shù)值參數(shù)參數(shù)值最小半徑/m0.1切向粘結力/N1e6粒徑比1.66接觸面摩擦系數(shù)0.5法向剛度/N·m-11e9顆粒密度/kg·m-32200切向剛度/N·m-16.67e8初始孔隙率0.16法向粘結力/N1e6

4.2 壓力拱計算結果分析

隧道體系在自重作用下,最大為豎向應力,隧道開挖后壓力拱效應也主要表現(xiàn)在拱頂區(qū)域豎向應力的變化。因此,分析時主要調(diào)取豎向應力分析壓力拱效應,如圖8所示。可以看出,隧道開挖后,在拱頂和仰拱底附近應力很小。這是由于該處圍巖發(fā)生塑性變形或者開裂,荷載由深部圍巖承擔。

圖8呈現(xiàn)了隧道開挖后,圍巖最終的應力狀態(tài),在拱頂和仰拱底形成了拱形的小應力區(qū),體現(xiàn)出壓力拱效應。為了更形象地呈現(xiàn)壓力拱的形成過程,將圖8中顯示的應力最終狀態(tài)減去開挖前的初始應力狀態(tài),繪制出隧道開挖過程中的應力變化,如圖9所示。可以看出,隧道開挖后,拱頂和仰拱底附近的圍巖應力急劇減小,邊墻附近圍巖應力升高,承載拱圈逐漸形成。

隧道施工后的壓力拱形成的主要原因是圍巖力的傳遞路徑的變化。在PFC2D程序中,可通過顆粒間接觸力的大小和方向的變化分析隧道開挖的壓力拱效應,如圖10所示。圖中,深色為壓力,淺色為拉力,線寬表示力的大小。可以看出,隧道開挖后,隧道周邊圍巖的力的傳遞方向發(fā)生偏轉,以隧道上方圍巖為例,力傳遞向了拱腰區(qū)域。從總體來看,隧道開挖后,拱頂區(qū)域力鏈很稀疏且線寬很細,表明該區(qū)域圍巖對于圍巖承載能力的貢獻很小,拱頂區(qū)域出現(xiàn)壓力拱現(xiàn)象。

圖8 豎向等效應力云圖

Fig.8 Equivalent vertical stress

圖9 豎向等效應力變化

Fig.9 Variety of equivalent vertical stress

進一步通過隧道開挖前后周邊圍巖離散顆粒接觸力的變化呈現(xiàn)圍巖力的變化和成拱效應,如 圖11 所示。圖中,紅色表示顆粒的接觸力降低,黑色表示顆粒的接觸力升高。

可以看出,

(1) 隧道開挖后,邊墻附近的圍巖顆粒的接觸力全部升高,且升高程度較大；拱頂和仰拱底附近的圍巖顆粒的接觸力既有升高也有降低,但升高程度小,降低程度大；拱頂和仰拱底顆粒升高的接觸力力鏈在邊墻處交匯。

(2) 在拱頂和仰拱底區(qū)域的圍巖顆粒的接觸力變化幅度沿徑向逐漸減小；降低的接觸力的連線基本上呈現(xiàn)拋物線形。

(3) 在靠近開挖面附近的拱頂和仰拱底區(qū)域,圍巖接觸力只有徑向接觸力降低而無環(huán)向接觸力升高,表明該區(qū)域圍巖自重荷載不能通過壓力環(huán)自身承載,只能通過顆粒間的粘結力承擔。而巖土材料的粘結能力有限,故該區(qū)域已出現(xiàn)失穩(wěn)坍塌。進而通過已有文獻研究成果，定義該區(qū)域為準塌落區(qū)域,其邊界為壓力拱的內(nèi)邊界。

4.3 壓力拱隨埋深的變化特點研究

壓力拱是當前判定隧道深淺埋的重要依據(jù)。根據(jù)文獻[2]對壓力拱的定義，壓力拱是力的傳遞路徑發(fā)生偏轉,應力集中的區(qū)域?？梢妷毫笆撬淼篱_挖后的主要承載圍巖。結合本文離散-耦合程序計算結果,目前已確定壓力拱的內(nèi)邊界,為了研究壓力拱隨埋深的變化特點,特定義壓力拱應力集中等效區(qū)域,即以1.15倍為應力增大系數(shù)確定外部邊界,以壓力拱內(nèi)邊界為內(nèi)部邊界,中間區(qū)域定義為壓力拱應力集中等效區(qū)域。

圖10 開挖后離散顆粒的力鏈

Fig.10 Force chain of discrete particles after excavation

圖11 開挖前后接觸力變化

Fig.11 Variety of contact force before and after excavation

對比計算了埋深分別為10 m,14 m,18 m,22 m,26 m,30 m和34 m情況下的隧道開挖產(chǎn)生的壓力拱應力集中等效區(qū)域,見圖12和表4。

由圖12和表4可見,隨著隧道埋深的增加,壓力拱應力集中等效區(qū)域內(nèi)邊界幾乎不變,外邊界不斷變淺,等效區(qū)域的面積逐漸減小,可見,隨著隧道埋深的增加,隧道周邊圍巖的壓力拱承載能力越來越大；當隧道埋深超過30 m時,壓力拱應力集中等效區(qū)域隨埋深的變化程度減小,表明當埋深達到一定程度后,隧道壓力拱高度區(qū)域穩(wěn)定。

表4 不同埋深下的壓力拱應力集中等效區(qū)域

Tab.4 Stress-concentration equivalent region of pressure arch under different buried depths

隧道埋深/m內(nèi)邊界深度/m外邊界深度/m103.319.73143.329.36183.328.13223.338.10263.337.73303.346.90343.356.89

圖12 不同埋深下壓力拱范圍的變化

Fig.12 Variation of region of pressure arch under different buried depths

5 結 論

本文采用二維顆粒流程序PFC2D,基于離散-連續(xù)耦合方法對隧道開挖后的圍巖受力特性進行了模擬,從微觀角度分析了圍巖壓力拱形成機理及發(fā)展規(guī)律,得到結論如下。

(1) 基于現(xiàn)有離散-連續(xù)耦合程序,提出通過把一定范圍的離散顆粒的平均應力作為外荷載的形式施加到有限單元上的方法實現(xiàn)離散-連續(xù)之間的耦合；平均應力計算區(qū)域面積通過在接觸單元的相鄰位置放置測量圓來確定,其直徑按照該單元的最大邊長取值；耦合程序遍歷顆粒過程通過將離散顆粒半徑進行放大,并且只遍歷連續(xù)單元附近三層的顆粒的方法降低計算耗時。

(2) 通過優(yōu)化的離散-連續(xù)耦合程序對隧道開挖后的圍巖受力特性計算結果可見,開挖面周邊范圍出現(xiàn)層狀壓力環(huán),但在拱頂和仰拱底附近區(qū)域的圍巖壓力環(huán)的環(huán)間接觸力降低,圍巖承載能力急劇降低,力向深部圍巖傳遞,導致淺部圍巖容易出現(xiàn)失穩(wěn),甚至發(fā)生塌落,形成壓力拱。

(3) 隨著隧道埋深的增加,隧道周邊圍巖的壓力拱承載能力越來越大；但當隧道埋深達到一定界限后,隧道壓力拱高度趨于穩(wěn)定。

本文主要研究了離散-連續(xù)耦合模型在隧道壓力拱研究中的可行性,該耦合程序既提高了計算效率，又能從微觀角度研究圍巖失穩(wěn)機理。但耦合程序的參數(shù)直接影響計算結果的精度。因此,在下一步研究中應結合隧道工程特點及圍巖特性,加強對離散部分材料及耦合邊界參數(shù)的研究。