裂隙白云巖流固耦合力學(xué)特性研究

劉德鵬， 饒軍應(yīng)， 王石光

(1.貴州省交通規(guī)劃勘察設(shè)計研究院股份有限公司， 貴陽 550081； 2.貴州大學(xué)土木工程學(xué)院， 貴陽 550081； 3.中鐵十八局第二工程有限公司， 唐山 064000)

近年來，隨著交通建筑行業(yè)的快速發(fā)展，導(dǎo)致了巖土、邊坡及隧道等工程方面的事故頻發(fā)。巖體本身是一種天然地質(zhì)材料，其內(nèi)部常蘊(yùn)藏裂隙、裂紋及節(jié)理等原始缺陷。地下巖體周圍通常分布有應(yīng)力場和滲流場，它們之間的相互作用使得裂隙、孔隙等萌生、擴(kuò)展、貫通，從而威脅工程的安全。然而，如何通過應(yīng)力-滲流耦合來描述巖體的受力狀態(tài)或力學(xué)特性的異常變化是保障工程設(shè)計、施工、管養(yǎng)安全與質(zhì)量的重要途徑之一。

流固耦合模型是一種將應(yīng)力場與滲流場結(jié)合的有效分析方法，表示巖體裂隙內(nèi)部流體與應(yīng)力的相互作用，學(xué)者們在這方面做了諸多多研究。劉佳佳[1]研究發(fā)現(xiàn)流固耦合模型可以描述瓦斯應(yīng)力變化規(guī)律；周曉敏等[2]研究發(fā)現(xiàn)降低滲透系數(shù)與泊松比對圍巖有效應(yīng)力控制具有重要作用；陳勇等[3]研究發(fā)現(xiàn)邊坡穩(wěn)定性與流固耦合作用有關(guān)；馬勇新[4]研究表明，流固耦合模型比傳統(tǒng)的非耦合模型更能表征高氣壓藏滲流與巖體變形特生；王偉[5]研究表明，流固耦合模擬結(jié)果與實際值接近，且該模型可以描述滑坡現(xiàn)象；高虎軍[6]研究表明，隧道襯砌應(yīng)力的變化與流固耦合現(xiàn)象有關(guān)；張建營[7]研究發(fā)現(xiàn)邊界壓應(yīng)力趨于穩(wěn)定時，耦合模型產(chǎn)生的動力效應(yīng)與實際相符；黃啟翔[8]研究發(fā)現(xiàn)耦合模型可以描述煤與瓦斯突出的力學(xué)作用機(jī)理。綜上所述，目前學(xué)者們多基于彈塑性原理、介質(zhì)的有效應(yīng)力原理、張量定理等開展對應(yīng)力-滲流耦合模型的研究，側(cè)重考慮水在裂隙中的滲流及其與裂隙巖體間應(yīng)力的耦合關(guān)系，而針對注漿液在巖體裂隙中的擴(kuò)散、凝結(jié)狀態(tài)對工程的影響及注漿液與巖體裂隙間的相互作用等方面的研究較少。

為此，基于多孔連續(xù)介質(zhì)的雙重有效應(yīng)力原理[9-12]，建立應(yīng)力-滲流的計算模型，采用控制變量法對比分析滲流時間t、流體密度ρ及應(yīng)力σ間的耦合作用規(guī)律，從而給出裂隙白云巖流固耦合作用下的力學(xué)特性，可為工程安全施工提供指導(dǎo)意見。

1 雙重有效應(yīng)力原理

劉忠玉等[13]、胡丹等[14]在一維固結(jié)模型基礎(chǔ)上并繼續(xù)研究基質(zhì)的力學(xué)和物理特性，得到了更為合理的三維固結(jié)模型理論；李傳亮[15-16]在研究土力學(xué)時發(fā)現(xiàn)了雙重有效應(yīng)力原理，即多孔介質(zhì)巖體變形具有結(jié)構(gòu)變形和本體變形，在進(jìn)行多孔介質(zhì)的應(yīng)力分析時，規(guī)定如下：一是多孔介質(zhì)發(fā)生結(jié)構(gòu)變形時，骨架體積視為不變，考慮其排列方式的變化，變形量用；二是多孔介質(zhì)發(fā)生本體變形時，骨架的體積發(fā)生變化，但其介質(zhì)分子的排列方式不變。以雙重有效應(yīng)力原理為依據(jù)開展對裂隙巖體應(yīng)力-滲流-時間耦合模型的研究。

基于多孔介質(zhì)的本體有效應(yīng)力，可表示為

(1)

基于多孔介質(zhì)的結(jié)構(gòu)有效應(yīng)力，可表示為

(2)

有效應(yīng)力σa為本體有效應(yīng)力與結(jié)構(gòu)有效應(yīng)力的疊加，可表示為

(3)

(1)當(dāng)φ→0時，介質(zhì)結(jié)構(gòu)較密，其物理化學(xué)性質(zhì)與相同介質(zhì)的固體材料性質(zhì)相似，此時有

(4)

(2)當(dāng)φ→1時，此時介質(zhì)材料比較疏松，它的物理化學(xué)性質(zhì)就不能采用固體的化學(xué)性質(zhì)來代替，必須通過其他手段來確定其力學(xué)性質(zhì)。則有

σa=δ-P

(5)

(3)根據(jù)雙重有效應(yīng)力原理可知，φ=0的介質(zhì)是不存在的。

(4)雙重有效應(yīng)力原理主要用于多孔介質(zhì)的巖體類結(jié)構(gòu)，將其應(yīng)用到流固耦模型之中，它能充分的考慮巖體的結(jié)構(gòu)變形和本體變形，比單純采用Terzaghi及Biot理論有效應(yīng)力更能反映巖體的力學(xué)性質(zhì)及變形機(jī)理。

2 應(yīng)力場方程

假設(shè)多孔介質(zhì)材料為各向同性、且有相等的彈性模量，其本構(gòu)關(guān)系方程[17]為

εa=Mδa，a=x、y、z

(6)

設(shè)ux、uy、uz分別為固體骨架x、y、z方向的位移，則得到幾何方程為

(7)

將式(7)代入式(6)可得

(8)

(1)本體有效應(yīng)力。利用矩陣疊加思想，視本體有效應(yīng)力為三維坐標(biāo)系下的應(yīng)力組合，即式(1)，將式(8)代入式(1)得

(9)

(2)結(jié)構(gòu)有效應(yīng)力。利用矩陣論思想，將結(jié)構(gòu)有效應(yīng)力看成是三維坐標(biāo)系下的應(yīng)力組合，即式(2)，將式(8)代入式(2)可得

(10)

由多孔連續(xù)介質(zhì)的雙重有效應(yīng)力原理知，總的有效應(yīng)力等于本體有效應(yīng)力與結(jié)構(gòu)有效應(yīng)力的之和，即將式(9)、式(10)代入式(3)得總的應(yīng)力分量[式(11)]。

(11)

式(11)中：σa為本體有效應(yīng)力與結(jié)構(gòu)應(yīng)力疊加應(yīng)力，即雙重有效應(yīng)力；σax、σay、σaz分別為x、y、z方向上的雙重有效應(yīng)力分量。

式(11)是以ux、uy、uz、P為變量的平衡方程。從而可知連續(xù)多孔介質(zhì)材料的雙重有效應(yīng)力是ux、uy、uz、P的函數(shù)，即σa=σa(ux、uy、uz、P)。

3 滲流場方程

設(shè)h為滲流場中各點(diǎn)的測管水頭，根據(jù)Darcy定律，則有一點(diǎn)水力坡度降在3個坐標(biāo)軸上的分量ix、iy和iz可表示為[18]

(12)

式(12)中：負(fù)號表示水力坡度降的正值對應(yīng)測管水頭降低的方向；ix、iy和iz為x、y、z方向水力坡度。

對于三維空間的滲流場，且?guī)r土體為各項同性，故可得到流體的滲流速度為

(13)

式(13)中：k為土體的滲透性系數(shù)；υx、υy、υz為x、y、z滲流速度。

在某一小微元內(nèi)，流體在滲流場中流速是穩(wěn)定的，其體積為dxdydz，在x、y、z方向的流速分別為υx、υy、υz。根據(jù)能量守恒原理，在其微元內(nèi)部必然存在一個內(nèi)源壓力q，使得整個滲流系統(tǒng)能量守恒。

根據(jù)雙重有效應(yīng)力原理，可得

(14)

流體連續(xù)性方程為

(15)

式(15)中：Θ為流體流量；t為時間。

流體狀態(tài)方程為

ρ=βρ0e-βP

(16)

測管漿液水頭等于壓力水頭與位置水頭之和，可表示為

(17)

流體體積含水率與孔隙水壓力關(guān)系為

(18)

將式(12)～式(18)代入Darcy定律可得

(19)

式中：ρ為流體密度，g/cm3；υ為滲流速度，m/s；ψ為水的動力黏度；ρ0為水的標(biāo)準(zhǔn)密度，ρ0=1 g/cm3；k為裂隙滲透系數(shù)，cm/s；β為水的壓縮系數(shù)，Pa-1，取β=4×10-10Pa-1；Z為標(biāo)高，m；mw為阻滯系數(shù)，一般情況取mw=1。

4 應(yīng)力-滲流-時間耦合模型

為便于研究應(yīng)力-滲流耦合模型與流體凝結(jié)或擴(kuò)散時間之間的關(guān)系，視流體具有流動性、可壓縮性、連續(xù)性等特點(diǎn)，且由式(12)可知，流體密度與流體壓強(qiáng)之間存在一定理論關(guān)系，從而根據(jù)雙重有效應(yīng)力原理知，應(yīng)力場與滲流場之間是相互影響的，即應(yīng)力場的存在會影響流體凝結(jié)效果及在裂隙里的擴(kuò)散效應(yīng)，而滲流場的存在也會反過來影響應(yīng)力場的性質(zhì)，具體如圖1所示。

圖1 應(yīng)力-滲流之間關(guān)系Fig.1 The relationship between stress and fluid seepage

結(jié)合式(17)、式(19)可推出式(20)，即

(20)

由式(14)得，依據(jù)矩陣間運(yùn)算，解矩陣方程得各應(yīng)力分量為

(21)

由式(16)流體狀態(tài)方程得

(22)

式(22)中：ρ0、β均為常數(shù)，在式(22)等式兩邊分別對x、y、z求二階導(dǎo)可得

(23)

結(jié)合式(20)～式(23)，記：

(24)

式(24)中：?、Δ為換算符號。

化簡得

(25)

式(25)即為流固-耦合方程，裂隙巖體在應(yīng)力-滲流耦合時是關(guān)于時間t、坐標(biāo)x、y、z的函數(shù)，即σ=σ(x、y、z、t)，可通過邊界條件及積分進(jìn)行求解。

5 模型求解及典例計算

5.1 模型求解

于實際工程項目而言，β、k、φ、φC、mw、ρ0、q、E及μ都是常數(shù)，由滲流知識可知，流體是位置坐標(biāo)的一次函數(shù)[19]，即

(26)

設(shè)A、B、C、D4個參數(shù)，其取值分別為

(27)

從而將式(25)化簡得

(28)

由E=σ/ε，將式(27)兩邊同時除以ε，即

(29)

取H=E2/D+BE/D，再次化簡得

H=σΔρ

(30)

(31)

流體微分表示在極短時間內(nèi)，可以將流體看作不變量，在某一瞬間，當(dāng)流體流動時間Δt→0時有

(32)

將式(32)代入式(17)并積分得

化簡得

(33)

從式(33)可知，σ與ρ、t有關(guān)，且流體密度越大，時間越大，流固耦合產(chǎn)生的應(yīng)力就越??；流體密度小，時間越小，此時耦合產(chǎn)生的應(yīng)力就越大。這與實際工程相符合，流體密度小、流體與圍巖作用時間短，此時產(chǎn)生的基本是圍巖初始應(yīng)力，被流體分擔(dān)的部分就變小，而流體具有稀釋作用，流體密度及作用時間的變少，就意味著，圍巖處于平衡狀態(tài)，此時初始應(yīng)力是很大的。

5.2 典例計算

以貴州境內(nèi)某隧道為依托案例，該隧道為雙洞單線式，全長925.411 m。資料顯示：ZDK26+143.2～+300段覆蓋層為塊石層與紅黏土層，下伏基巖為松子坎二段泥質(zhì)白云巖，巖體較破碎，且節(jié)理裂隙較發(fā)育，呈淺灰色、灰白色，中厚層狀，存在大量裂隙，其巖體物理力學(xué)性質(zhì)如表1所示。

根據(jù)依托工程可知，裂隙白云巖的孔隙率φ為2%，有效孔隙率φC為0.5%，阻滯系數(shù)mw為1，將其已知參數(shù)代入式(27)，可得相關(guān)的物理參數(shù)計算如表2所示。

表1 巖體物理力學(xué)參數(shù)Table 1 Physical and mechanical parameters of rock mass

根據(jù)已知數(shù)據(jù)及σ、ρ、t三者之間的關(guān)系，以依托工程為基礎(chǔ)進(jìn)行數(shù)值模擬分析，對比分析在不同條件下的應(yīng)力隨時間、流體密度的變化規(guī)律及找出注漿液的最佳凝結(jié)時間，故采用相關(guān)編程軟件，繪制出應(yīng)力云圖(圖2)。

M1～M9、E1～E10、F1～F10為應(yīng)力轉(zhuǎn)點(diǎn)圖2 應(yīng)力-滲流耦合云圖Fig.2 Stress-seepage coupling cloud diagram

由應(yīng)力云圖(圖2)可知，耦合模型的應(yīng)力大小與流體密度及時間均有關(guān)。應(yīng)力值隨著時間t、流體密度ρ的增加而增大，應(yīng)力增大到一定范圍后，應(yīng)力值波動范圍很小，甚至不變。

為更直觀的判斷應(yīng)力-滲流耦合時間t及流體密度ρ對耦合應(yīng)力的影響情況，這里采用控制變量法進(jìn)行研究，即控制流體密度不變，改變耦合時間t，繪制σ-t曲線，如圖3所示；控制耦合時間不變，改變流體密度，繪制σ-ρ曲線，如圖4所示。

由圖3可知：流體密度不變時，耦合模型的應(yīng)力值隨耦合時間的增大而增大，且作用時間越長，產(chǎn)生的應(yīng)力峰值就越小。當(dāng)t=0.02 d時，每種流體密度在σ-t曲線上均出現(xiàn)應(yīng)力最大值，在此之后，應(yīng)力值基本不隨時間變化而變化。

σ1max、σ2max、σ3max、σ4max分別為在t=0.02 d時轉(zhuǎn)折點(diǎn) 1、2、3、4對應(yīng)的峰值應(yīng)力圖3 σ-t曲線示意圖Fig.3 Schematic diagram of σ-t curve

σA′、σB′、σC′、σD′分別為在分界點(diǎn)A′、B′、C′、D′對應(yīng)的耦合應(yīng)力； ρA′、ρB′、ρC′、ρD′分別為在分界點(diǎn)A′、B′、C′、D′對應(yīng)的流體密度圖4 σ-ρ曲線示意圖Fig.4 Schematic diagram of σ-ρ curve

由圖4可知，在控制時間不變情況下，耦合模型產(chǎn)生的應(yīng)力值隨流體密度的增加而減少，當(dāng)增加ρ為0.8 g/cm3時，應(yīng)力值基本不變。當(dāng)ρ為0.2 g/cm3時，每種時間對應(yīng)下的σ-ρ曲線均產(chǎn)生應(yīng)力最大值。

表2 裂隙白云巖相關(guān)參數(shù)Table 2 Related parameters offractured dolomite

為確定應(yīng)力開始變緩時的時間t及流體密度ρ，進(jìn)行數(shù)值分析時分別取了3個不同應(yīng)力增幅進(jìn)行研究，且將數(shù)值解和理論解進(jìn)行對比，如圖5所示。

ρF1、 ρF2分別為在標(biāo)記點(diǎn)F1、F2對應(yīng)的流體密度； σF1、σF2、σH3分別為在標(biāo)記點(diǎn)F1、F2、H3對應(yīng)的耦合應(yīng)力； tH3為標(biāo)記點(diǎn)H3對應(yīng)的時間圖5 各區(qū)域應(yīng)力變化Fig.5 Stress changes in each area

圖5(a)表明：在負(fù)值增值幅區(qū)域，應(yīng)力的大小隨流體密度的增加而增加，當(dāng)流體密度增加到0.8 g/cm3時，此后應(yīng)力增幅變緩，應(yīng)力基本不變。說明應(yīng)力值σ隨流體密度ρ的增加不是無限制的，結(jié)合圖3可知，無論耦合時間時長多少，流體密度ρ為0.8 g/cm3為該耦合模型應(yīng)力增幅基本不變的分界點(diǎn)，即H1點(diǎn)。

由圖5(b)可知，在正值增值幅區(qū)域，應(yīng)力值隨流體密度的增加而增加，當(dāng)流體密度ρ為0.8 g/cm3后，應(yīng)力增幅變緩，應(yīng)力趨勢線接近水平線。說明應(yīng)力值σ隨流體密度ρ的增加不是無限制的，結(jié)合圖3可知，無論耦合時間為多少，流體密度ρ為0.8 g/cm3為該耦合模型應(yīng)力增幅基本不變的分界點(diǎn)，即H2點(diǎn)。

由圖5(c)可知，在正負(fù)應(yīng)力增值區(qū)域，應(yīng)力的大小隨耦合時間的增加而增加，當(dāng)t為0.02 d后，此時應(yīng)力達(dá)到峰值，應(yīng)力增幅基本不變，應(yīng)力趨勢線也接近水平線。說明應(yīng)力值σ隨耦合作用時間t的增加是有限的，且結(jié)合圖2可知，無論流體密度ρ為多少，t為0.02 d時為該耦合模型應(yīng)力增幅基本保持不變的分界點(diǎn)，即H3點(diǎn)。

圖3～圖5共同說明，應(yīng)力值隨流體密度ρ、時間t的增加而增大，且在t為0.02 d、ρ為0.8 g/cm3之后，應(yīng)力值基本保持不變，應(yīng)力趨勢線也近似看作是水平線。

6 結(jié)論

(1)應(yīng)力場與滲流場方程理論推導(dǎo)過程，考慮巖體變形特征、流體密度、位置坐標(biāo)、流體壓力及滲流時間等因素，采用流體密度變分方法，建立了應(yīng)力-滲流-時間等多場耦合理論分析模型。

(2)經(jīng)理論分析發(fā)現(xiàn)，應(yīng)力-滲流-時間多場耦合時，應(yīng)力、流體密度及滲流時間是相互制約、相互影響的，且存在數(shù)值對應(yīng)關(guān)系。

(3)通過實例分析，流體密度及時間使裂隙白云巖力學(xué)特性發(fā)生較大變化，且當(dāng)0.01 d

(4)該模型在前人研究基礎(chǔ)上，增加了時間變量、巖體結(jié)構(gòu)變形及本體變形本構(gòu)關(guān)系和流體密度對時間、位置坐標(biāo)的變分等對模型建立的影響，但該模型視滲透系數(shù)、裂隙寬度、孔隙度為常數(shù)，后續(xù)研究可將滲透系數(shù)、裂隙寬度、孔隙度等考慮在內(nèi)，其圍巖應(yīng)力變化情況有待研究。

(5)雖然對貴陽某地區(qū)裂隙白云巖力學(xué)特性進(jìn)行數(shù)值分析，但本文模型的理論推導(dǎo)是基于一般裂隙巖體，所以該耦合模型也可為其他裂隙巖體流固耦合模型理論推導(dǎo)提供參考。