基于連續(xù)介質(zhì)離散元的雙重介質(zhì)滲流應(yīng)力耦合模型

基于連續(xù)介質(zhì)離散元的雙重介質(zhì)滲流應(yīng)力耦合模型

1數(shù)值方法與數(shù)值模型由于水庫水和降雨的滑動(dòng)，不可避免地會(huì)涉及水的作用，反映在計(jì)算值時(shí)的滲流電壓耦合模型中。它分為2種:第一種是非雙重介質(zhì)滲流應(yīng)力耦合模型,它包括等效連續(xù)介質(zhì)滲流應(yīng)力耦合模型和裂隙網(wǎng)絡(luò)滲流應(yīng)力耦合模型,前者適用于孔隙介質(zhì)以及多組裂隙分割的破裂巖體,后者適用于水流只在裂隙網(wǎng)絡(luò)中的定向流動(dòng)。第二種是雙重介質(zhì)滲流應(yīng)力耦合模型,它適用于把巖土體看作由孔隙和裂隙組成的雙重介質(zhì)空隙結(jié)構(gòu)。大量工程實(shí)例發(fā)現(xiàn),庫區(qū)的古滑坡具有一條非常明顯的滑帶,把邊坡分為了上部松軟、透水的滑體和下部堅(jiān)硬、隔水的基巖。對(duì)于這類庫岸型古滑坡,人們進(jìn)行了很多的研究。理論方面,S.W.Qi等探討了由于庫水漲落所誘發(fā)滑坡的破壞機(jī)制。同普通的山地滑坡相比,除降雨外,由于水利工程而導(dǎo)致的外部水環(huán)境改變是這類滑坡所特有的。強(qiáng)度軟化作用,浮托作用,動(dòng)水壓力作用,固結(jié)作用,庫岸再造,降雨等都會(huì)不同程度影響庫岸型古滑坡的穩(wěn)定性。當(dāng)然,庫水漲落和降雨等都是外因,滑坡變形的最本質(zhì)原因是其特殊的地質(zhì)結(jié)構(gòu)和巖土體特性。數(shù)值模擬能夠較好地揭示滑坡破壞的動(dòng)態(tài)過程,同時(shí)數(shù)值模擬也能夠用來進(jìn)行敏感性分析,例如滑坡地形地貌,巖土體物理力學(xué)參數(shù),庫水漲落的速率等。庫水漲落和降雨所產(chǎn)生的水流通過坡體表面的裂縫迅速入滲,沿著滑帶迅速地流動(dòng),因此可以把裂縫和滑帶看成是裂隙介質(zhì)。但是滑體還是具有一定的滲透性,需要把其看成是孔隙介質(zhì)。所以模擬庫區(qū)邊坡在降雨和庫水漲落作用下的響應(yīng),選取雙重介質(zhì)滲流應(yīng)力耦合模型是比較符合實(shí)際情況的。目前國內(nèi)外對(duì)于雙重介質(zhì)滲流應(yīng)力耦合研究取得了一定的進(jìn)展。劉耀儒等采用基于elementby-element策略的有限元,開展了雙重介質(zhì)滲流應(yīng)力耦合模型并行數(shù)值分析;黎水泉和徐秉業(yè)考慮固結(jié)對(duì)滲流的影響,同時(shí)也考慮固體變形對(duì)滲流參數(shù)(孔隙度和滲透率)的影響,推導(dǎo)出了雙重孔隙介質(zhì)非線性流固耦臺(tái)滲流有限元計(jì)算格式;孔亮等利用雙重孔隙介質(zhì)模型,在多孔介質(zhì)彈性理論和廣義有效應(yīng)力的基礎(chǔ)上,給出了雙重孔隙介質(zhì)非飽和流固耦合方程的數(shù)學(xué)推導(dǎo)。G.I.Barenblatt等在1960年提出了均質(zhì)、各向同性的雙重孔隙介質(zhì)模型;V.Reichenberger等提出了一種基于點(diǎn)–中心有限體積法的全耦合、全隱式離散的雙重孔隙介質(zhì)模型。J.Rutqvist等結(jié)合Tough2和FLAC3D實(shí)現(xiàn)了雙重孔隙介質(zhì)的熱–水–力–化學(xué)耦合。連續(xù)介質(zhì)離散元(continuum-mediumdistinctelementmethod,CDEM)是近期發(fā)展起來的一種數(shù)值計(jì)算方法,它將任意塊體(平面)單元化作一系列四面體(三角形)單元,根據(jù)各個(gè)四面體(三角形)的單元?jiǎng)偠染仃?就可以寫出任意塊體(平面)的剛度矩陣,而塊體(平面)之間不組成總體剛度矩陣。采用動(dòng)態(tài)松弛技術(shù)來求解各個(gè)塊體(平面)單元的節(jié)點(diǎn)水頭、力和位移等變量。針對(duì)上述情況,本文發(fā)展了一種基于連續(xù)介質(zhì)離散元的雙重介質(zhì)滲流應(yīng)力耦合數(shù)值模型,它分為固體計(jì)算模型、孔隙滲流應(yīng)力耦合模型和裂隙滲流應(yīng)力耦合模型3個(gè)部分。根據(jù)定義,它包含2個(gè)互相覆蓋的連續(xù)場(chǎng)(見圖1):基質(zhì)系統(tǒng)和裂隙系統(tǒng)。每個(gè)連續(xù)場(chǎng)都擁有自己的滲流壓力場(chǎng)。2個(gè)系統(tǒng)之間有流體的交換,交換的程度取決于2個(gè)系統(tǒng)之間的水頭差。在該模型中,任意一點(diǎn)處都有2個(gè)水頭,即基質(zhì)系統(tǒng)水頭和裂隙系統(tǒng)水頭。固體計(jì)算模型分為材料單元和接觸單元。材料單元被看作是彈性體或者彈塑性體,在其內(nèi)部根據(jù)力邊界條件用有限元求各點(diǎn)的位移;接觸單元引入法向和切向彈簧根據(jù)位移和相對(duì)速度求力。固體計(jì)算的模型能夠反映地質(zhì)體的破壞規(guī)律,可以模擬從連續(xù)到非連續(xù)的破壞過程。在孔隙、裂隙滲流應(yīng)力耦合模型計(jì)算中,由于引入了CDEM,自由水位線(浸潤線)的計(jì)算也會(huì)變得更加方便,并且可以避免由于不連通裂隙(孤立裂隙)存在所導(dǎo)致的收斂性問題。同時(shí),由于采用了動(dòng)態(tài)松弛技術(shù),即迭代弱耦合方法,以上一時(shí)步孔隙滲流場(chǎng)(裂隙滲流場(chǎng))的結(jié)果作為下一時(shí)步裂隙滲流場(chǎng)(孔隙滲流場(chǎng))的輸入條件,從而實(shí)現(xiàn)孔隙滲流場(chǎng)和裂隙滲流場(chǎng)的更新。當(dāng)時(shí)步足夠小時(shí),計(jì)算能夠收斂到精確解。為驗(yàn)證模型的可靠性及適用性,本文將對(duì)一個(gè)典型算例進(jìn)行計(jì)算及討論。2單元?jiǎng)偠染仃嚨慕⒛Ｐ椭袎K體系統(tǒng)力學(xué)分析是指單個(gè)塊體的力學(xué)分析和相鄰塊體之間的相互作用。塊體小位移彈性動(dòng)態(tài)方程為應(yīng)變與位移關(guān)系:本構(gòu)關(guān)系:式(1)~(3)中:σij,j為應(yīng)力張量對(duì)長度的一階偏導(dǎo)數(shù);fi為單位體積的體力;ρ為密度;ui,tt為位移對(duì)時(shí)間的二階偏導(dǎo)數(shù);ui,t為位移對(duì)時(shí)間的一階偏導(dǎo)數(shù);α為阻尼系數(shù);ui,j,uj,i均為位移對(duì)長度的一階偏導(dǎo)數(shù);εij,εst均為應(yīng)變;σij為應(yīng)力;Cijst為Hooke張量。塊體被離散為有限數(shù)量的四面體單元,通過節(jié)點(diǎn)相連。在每個(gè)四面體單元中,包含一個(gè)塊體中心節(jié)點(diǎn),一個(gè)面中心節(jié)點(diǎn),棱邊端部的2個(gè)節(jié)點(diǎn)。塊體剛度矩陣由各個(gè)四面體單元?jiǎng)偠染仃嚱M合而成。在固體計(jì)算模型中,采用動(dòng)態(tài)松弛方法,由于引入慣性項(xiàng),力學(xué)過程事實(shí)上是動(dòng)態(tài)的,但可以通過引入人工阻尼項(xiàng)使得這一過程變?yōu)闇?zhǔn)靜態(tài)過程。塊體節(jié)點(diǎn)力的計(jì)算不像有限元方法那樣將節(jié)點(diǎn)周圍的塊體剛陣組合成總剛陣計(jì)算,而是按照各個(gè)塊體的剛度矩陣分別計(jì)算出屬于各自塊體的節(jié)點(diǎn)力。相鄰塊體之間依靠彈簧連接(見圖2),采用的破壞準(zhǔn)則為莫爾–庫侖準(zhǔn)則。由連續(xù)的單元邊界面轉(zhuǎn)化為非連續(xù)斷裂面的計(jì)算方法。計(jì)算單元界面兩側(cè)單元應(yīng)力(等參單元),并計(jì)算出在單元邊界上的法向和切向應(yīng)力。如果其中一個(gè)單元的應(yīng)力滿足破壞條件:式中:σn為法向應(yīng)力;σr為結(jié)構(gòu)面材料的容許拉應(yīng)力;σt為切向應(yīng)力;c,?分別為黏聚力和內(nèi)摩擦角。將單元界面轉(zhuǎn)化成為斷裂面,斷裂面上的接觸點(diǎn)的相對(duì)位移滿足:式中:Δun,Δut分別為法向位移和切向位移;nF,tF分別為法向力和切向力;σn1,σn2均為相鄰接觸點(diǎn)的法向應(yīng)力;Kn,Kt分別為彈簧的法向剛度和切向剛度;A為接觸點(diǎn)所代表的面積;σt1,σt2為相鄰接觸點(diǎn)的切向應(yīng)力。新定義的斷裂面改變了原來連續(xù)界面的特性,應(yīng)力場(chǎng)將會(huì)在繼續(xù)計(jì)算的過程中重新調(diào)整,此時(shí)的斷裂面可能還沒有完全破壞,直至滿足斷裂面拉伸破壞的條件:此時(shí),若斷裂面上的受力狀態(tài),滿足剪切破壞條件:此時(shí),由于動(dòng)態(tài)松弛方法不需要形成總體剛度矩陣,與傳統(tǒng)有限元法(FEM)相比,所需內(nèi)存較小。動(dòng)態(tài)松弛方法還可以求解裂隙巖體的變形和應(yīng)力,而傳統(tǒng)有限元在處理此類問題時(shí)將會(huì)遇到很大的困難。3孔絡(luò)流的開口3.1單個(gè)四面體滲流單元的滲流矩陣解析在高飽和度非飽和巖土體中,氣體以封閉的氣泡形式存在于孔隙水中,可認(rèn)為是一種可壓縮的流體,則式中:S為單位貯水量,表示下降單位水頭時(shí),為由于骨架壓縮和水的膨脹所釋放出的貯存水量,與巖土體和孔隙流體的壓縮系數(shù)有關(guān);kx,ky和kz分別為x,y和z方向的滲透系數(shù);h為總水頭。式(12)所示的微分方程在復(fù)雜的邊界條件下無法直接得到解析解,數(shù)值計(jì)算的方法是首先建立h(x,y,z,t)的泛函,一定邊值問題的解就是這個(gè)泛函的極小值,這個(gè)求解過程就是變分。在單個(gè)四面體滲流單元Ωe中,考慮液體的可壓縮性,對(duì)于非穩(wěn)定滲流,泛函Ee(h)可表示為式中:q為流量,s為表面積,Γ2e為表面區(qū)域。由于所要求的h(x,y,z,t)為泛函式(13)的極小函數(shù),因此必須滿足:因此可以導(dǎo)出單個(gè)四面體滲流單元的滲流矩陣方程為式中:[Se]為單元的貯存系數(shù)矩陣,[Ke]為單元滲透矩陣,{Qe}為根據(jù)流量邊界得到的節(jié)點(diǎn)流量矩陣。程序計(jì)算時(shí),同樣不組成總體的滲透矩陣,以單個(gè)孔隙單元為單位,采用動(dòng)態(tài)松馳技術(shù)來求解每個(gè)孔隙單元的節(jié)點(diǎn)水頭。3.2離散單元法原理由于滲流自由面的位置是待求的,故必須迭代求解,滲流自由面問題的求解方法主要分為水動(dòng)力學(xué)法、流網(wǎng)法、試驗(yàn)法和有限單元法。在這里采用離散單元法:(1)第一次孔隙滲流計(jì)算,全區(qū)域的滲透系數(shù)為給定的滲透系數(shù);(2)將節(jié)點(diǎn)水頭小于其位置勢(shì)的所有節(jié)點(diǎn)進(jìn)行標(biāo)識(shí),然后判斷處于自由面之上的所有單元;(3)將自由面以上單元的滲透系數(shù)給定一個(gè)小值:式中:kn為修正滲透系數(shù),kn-1為初始滲透系數(shù)。在計(jì)算過程中不再考慮自由面以上單元滲流的影響;(4)將本次求出的結(jié)點(diǎn)勢(shì)與上一次迭代求出的結(jié)點(diǎn)勢(shì)比較,判斷:式中:i為迭代計(jì)算次數(shù),j為節(jié)點(diǎn)編號(hào),δ為水頭誤差。若所有節(jié)點(diǎn)都滿足式(17),滲流場(chǎng)中所有滿足|H-z|≤δ的節(jié)點(diǎn)的連線即為自由面。這種方法借鑒了有限單元法中的變單元滲透系數(shù)固定網(wǎng)格法,離散單元法在處理浸潤線問題上有著比較突出的優(yōu)勢(shì),由于各個(gè)單元之間不需要形成總滲透系數(shù)矩陣,可以很方便地改變?cè)谧杂擅嬷系膯卧獫B透系數(shù)。假設(shè)滲流沿著正y方向(見圖3),從水平滲流的飽和土體中取出一個(gè)巖土體單元,長、寬、高分別為dx,dy,dz;設(shè)點(diǎn)1,2,3,4有相等的總水頭HA,點(diǎn)5,6,7,8有相等的總水頭HB,面3478為基準(zhǔn)面。則各點(diǎn)的孔隙水壓力分別為hw1=hw2=HA-dz,hw3=hw4=HA,hw5=hw6=HB-dz,hw7=hw8=HB。3.3fy、fz的計(jì)算沿著x,y,z方向的合力分別為:Fx=0,Fy=ρwg(HA-HB)dxdz,Fz=ρwgdxdydz。換算成力:,fz=ρwg,其中fy為滲透力,fz為浮力。3.4影響土滲透系數(shù)的因素孔隙滲流的滲透系數(shù)由以下公式給出:對(duì)于飽和土,滲透系數(shù)k取決于流體的動(dòng)力黏滯系數(shù)η,而η又與溫度、介質(zhì)濃度和水中氣泡等因素有關(guān),γw為流體的重度,與可溶性鹽的含量有關(guān)。此外,用毛細(xì)管代替孔隙通道,半徑為R,n為孔隙率。R和n受土粒的顆粒大小、級(jí)配等因素影響,這些也是影響土的滲透系數(shù)的因素。在程序中,根據(jù)下式計(jì)算不同時(shí)步巖土體的滲透系數(shù),實(shí)現(xiàn)應(yīng)力場(chǎng)對(duì)滲流場(chǎng)的影響:式中:k0為σ=0時(shí)的滲透系數(shù),σ為有效應(yīng)力,α為待定系數(shù)。4裂縫分辨率的力耦合4.1裂隙水單元單元假設(shè)研究區(qū)域初始為連續(xù)區(qū)域,劃分完塊體單元以后,當(dāng)單元界面內(nèi)所有彈簧都發(fā)生斷裂破壞時(shí),就把這個(gè)單元界面變?yōu)榱严端畣卧?。通過數(shù)值計(jì)算中單元界面斷裂計(jì)算模型,實(shí)現(xiàn)模擬地質(zhì)體實(shí)際的演化過程。4.2維滲流計(jì)算對(duì)于單裂隙來說,裂隙的隙寬t遠(yuǎn)遠(yuǎn)小于裂隙面的長度和寬度,因此可以將三維問題轉(zhuǎn)化為二維問題進(jìn)行計(jì)算。根據(jù)質(zhì)量守恒和達(dá)西定律推出非飽和滲流的基本微分方程:式中:D為二維滲流計(jì)算區(qū)域。類似于孔隙滲流場(chǎng)的計(jì)算過程,將裂隙水單元離散化,采用動(dòng)態(tài)松馳技術(shù)來求解每個(gè)裂隙水單元的節(jié)點(diǎn)水頭。4.3裂隙壁面滲透靜水壓力的表現(xiàn)由于在實(shí)際情況中,裂隙水沿著裂隙迅速地流動(dòng),這個(gè)過程相對(duì)于孔隙滲流來說相當(dāng)快。只考慮裂隙水流對(duì)裂隙壁法向的滲透靜水壓力作用,不考慮切向的拖曳力(滲透動(dòng)水壓力)作用,當(dāng)單裂隙中沒有被充填時(shí),裂隙水流對(duì)裂隙壁的滲透靜水壓力是面力,其方向垂直于裂隙壁面,對(duì)裂隙產(chǎn)生擴(kuò)張作用。裂隙壁面滲透靜水壓力p的表達(dá)式為式中:H,Hz分別為滲流的總水頭和位置水頭;ρ為水的密度;g為重力加速度常數(shù)。4.4應(yīng)力場(chǎng)滲流場(chǎng)的計(jì)算式中:q為裂隙的單寬流量;b為裂隙寬度;γ,μ分別為流體的容重和動(dòng)力黏滯系數(shù);Jf為水力梯度。因此可根據(jù)下式實(shí)現(xiàn)應(yīng)力場(chǎng)對(duì)裂隙滲流場(chǎng)的影響:式中:k0,b0分別為初始時(shí)刻的滲透系數(shù)和裂隙寬度;b為任意時(shí)刻的裂隙寬度。5孔隙滲流場(chǎng)和裂隙滲流場(chǎng)耦合的計(jì)算在實(shí)際計(jì)算中,庫水漲落和降雨引起水頭改變的同時(shí),也可作為孔隙滲流和裂隙的變邊界條件。相比較于孔隙滲流而言,裂隙滲流或者是庫水漲落的速度非常快。上一個(gè)時(shí)步中,在庫水漲落和降雨引起水頭改變的情況下只計(jì)算裂隙滲流,產(chǎn)生的裂隙滲流場(chǎng)水頭分布和由庫水漲落所引起的水頭改變作為孔隙滲流的變邊界條件,在后面的時(shí)步中,裂隙滲流不參與計(jì)算,而只有孔隙滲流參與計(jì)算,所以可以做如下假設(shè):只考慮降雨和庫水漲落的最終狀態(tài),裂隙滲流產(chǎn)生的水頭分布和由庫水漲落所引起的水頭改變作為孔隙滲流的變邊界條件,從而實(shí)現(xiàn)孔隙滲流場(chǎng)和裂隙滲流場(chǎng)的耦合。相當(dāng)于第一類邊界條件中變水頭邊界:式中:為孔隙滲流水頭邊界條件,H1為裂隙滲流產(chǎn)生的水頭分布,B1為孔隙滲流場(chǎng)和裂隙滲流場(chǎng)耦合的區(qū)域。以庫區(qū)邊坡(見圖4)為例,滲流邊界條件為(1)距離熱水的極限(2)流量限制(3)自由水域(4)滲面邊界面qn式(26)~(28)中:h為總水頭,為給定水頭,kij為二階滲透張量,ni為所在面的外法線方向的單位向量,qn為給定流量,Γh為水頭邊界區(qū)域,Γq為流量邊界區(qū)域,Γf為自由面邊界區(qū)域,Γs為潛在出滲面邊界區(qū)域。6降雨模擬結(jié)果為驗(yàn)證程序的正確性,引入一個(gè)庫水漲落和降雨影響下的古滑坡模型。不同于純土體坡體,古滑坡的滑帶通常就是滑動(dòng)面,而不需要通過計(jì)算得出圓弧滑動(dòng)面,當(dāng)然古滑坡滑體部分也能夠通過計(jì)算出現(xiàn)次級(jí)滑動(dòng)面。本文中弧形滑帶是事先網(wǎng)格劃分好的,代表古滑坡的滑動(dòng)面(見圖5)。棋盤式網(wǎng)格劃分,并不代表裂隙。裂隙的產(chǎn)生由單元界面斷裂計(jì)算模型實(shí)現(xiàn),在計(jì)算中,當(dāng)達(dá)到一定的破壞準(zhǔn)則時(shí),裂隙是沿著單元網(wǎng)格邊界產(chǎn)生,而不通過塊體內(nèi)部破壞來產(chǎn)生。當(dāng)網(wǎng)格劃分足夠密時(shí),即使是棋盤式的網(wǎng)格劃分,也能夠很好地反映坡體內(nèi)部的破壞規(guī)律,更為重要的是這種方法易于程序的編寫,并且可節(jié)省大量的計(jì)算時(shí)間。降雨的模擬,只考慮降雨的最終狀態(tài),不模擬降雨的過程,即整個(gè)過程是穩(wěn)態(tài)流過程。裂隙(孔隙)出露在坡體外部的節(jié)點(diǎn)為邊界節(jié)點(diǎn),邊界條件為水頭邊界條件。本文參與滲流計(jì)算的總水頭由2個(gè)部分組成:位置水頭和壓力水頭。降雨時(shí),裂隙(孔隙)邊界節(jié)點(diǎn)的水頭邊界條件為位置水頭(零水頭)。不管是降雨還是庫水漲落,只要存在水的作用,巖土塊體模型和彈簧模型的黏聚力和內(nèi)摩擦角都降低30%。計(jì)算模型的參數(shù)見表1,2所示,計(jì)算模型網(wǎng)格劃分如圖6所示。利用程序分別對(duì)古滑坡在自重作用下、庫水驟降情況下以及庫水驟降和暴雨聯(lián)合作用下進(jìn)行數(shù)值模擬,得出以下結(jié)論:(1)從圖7(a),8(a)和9(a)可以看出,庫水驟降時(shí)邊坡內(nèi)部彈簧破裂較多,而庫水上漲時(shí)邊坡內(nèi)部彈簧破裂較少。這是由于相比較于庫水上漲,庫水驟降更不利于庫區(qū)邊坡的穩(wěn)定性。(2)庫水上漲對(duì)邊坡會(huì)產(chǎn)生2個(gè)影響:一是付托力,二是降低黏聚力、內(nèi)摩擦角等物理參數(shù);后者肯定不利于邊坡穩(wěn)定性;而前者則不能一概而論。如果浮托力只是影響邊坡的抗滑段,則不利于邊坡穩(wěn)定,如果付托力影響下滑段的范圍遠(yuǎn)遠(yuǎn)大于阻滑段,那么就有利于邊坡穩(wěn)定。(3)庫水驟降除產(chǎn)生浮托力和降低巖土體物理參數(shù)外,對(duì)邊坡還會(huì)產(chǎn)生第三個(gè)影響,即滲透力;一般情況下,滲透力對(duì)邊坡的影響占主導(dǎo)作用。(4)比較圖7(b),8(b)和9(b),可以發(fā)現(xiàn)與自重作用下相比,邊坡在庫水上漲和驟降情況下,位移變化不大,邊坡只是在內(nèi)部產(chǎn)生漸進(jìn)破壞,沒有最終失穩(wěn)。(5)在漲水時(shí),邊坡體內(nèi)浸潤線位置大致與庫水水位線齊平(見圖8(c)),水對(duì)巖土體材料只有浮托力,無滲透力;而驟降時(shí),浸潤線位置呈一條曲線分布(見圖9(c)