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

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

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