土石壩滲流有限元分析與壩坡穩(wěn)定性分析

土石壩滲流有限元分析與壩坡穩(wěn)定性分析

1抗剪強(qiáng)度折減系數(shù)及壩體穩(wěn)定性分析近年來，隨著計(jì)算機(jī)軟件和硬件的發(fā)展，使用有限差的穩(wěn)定性分析越來越受到重視。許多工程研究人員表示，使用有限差分穩(wěn)定性分析容易受到使用。與基于塑性極限平衡原理的傳統(tǒng)方法相比,有限元法能夠考慮復(fù)雜的邊界條件,不受材料不均勻的限制,也不需預(yù)先假定滑裂面,能模擬滑坡的可能運(yùn)動(dòng)方向,確定破壞區(qū)的位置和破壞發(fā)展情況,同時(shí)在理想的彈塑性理論框架下全面而又嚴(yán)密地考慮了土的本構(gòu)關(guān)系。強(qiáng)度折減法是邊坡穩(wěn)定有限元分析中比較常用的方法,它將極限平衡法的安全系數(shù)定義與有限元法相結(jié)合,按照某一系數(shù)sF逐步將土的實(shí)際抗剪強(qiáng)度參數(shù)進(jìn)行折減,基于折減后的抗剪強(qiáng)度進(jìn)行一系列的有限元分析,直至所選定的抗剪強(qiáng)度折減系數(shù)sF使得邊坡達(dá)到臨界失穩(wěn)狀態(tài),此時(shí)的sF即為安全系數(shù)[2,3,4,5,6,7,8,9,10,11]。通過這種方法還可同時(shí)得到臨界滑裂面的位置。目前在通常的有限元程序中實(shí)現(xiàn)強(qiáng)度折減有限元法大致有兩種方法:一種常規(guī)方法是先折減后加載,這種方法概念明確,計(jì)算可靠,但需要反復(fù)修改輸入文件中的材料強(qiáng)度參數(shù)c值和f值進(jìn)行試算,計(jì)算量大;另一種方法是先加載后折減,這種方法以現(xiàn)實(shí)的應(yīng)力狀態(tài)為起點(diǎn),逐步降低土體的強(qiáng)度參數(shù)值進(jìn)行迭代計(jì)算,直至邊坡進(jìn)入極限狀態(tài),由此得出相應(yīng)的安全系數(shù)。第2種方法更加符合實(shí)際情況,整個(gè)分析計(jì)算一次完成,不需要調(diào)整輸入文件反復(fù)試算,但需要重新編制有限元程序才能實(shí)現(xiàn)。徐千軍等結(jié)合兩種方法的長處,在邊坡穩(wěn)定分析中提出采用溫控參數(shù)折減有限元法進(jìn)行邊坡穩(wěn)定分析。這種分析方法,參數(shù)一次輸入,不需要重新編寫程序,大大提高了計(jì)算效率。對于土石壩的實(shí)際邊坡穩(wěn)定性分析一般需要合理地考慮土與水的相互效應(yīng)和滲流作用。李湛等將滲流計(jì)算的有限元分析與邊坡穩(wěn)定性分析的強(qiáng)度折減彈塑性有限元法相結(jié)合,提出了滲流作用下邊坡穩(wěn)定性的強(qiáng)度折減有限元法。這種方法雖然能夠通過有限元滲流計(jì)算較為合理地考慮滲流作用力得到安全系數(shù),但仍屬于先折減后加載的常規(guī)有限元強(qiáng)度折減的做法,計(jì)算過程中需要反復(fù)修改輸入文件中的材料強(qiáng)度參數(shù)c值和f值進(jìn)行試算,計(jì)算量大,不易得到壩體穩(wěn)定安全系數(shù)。本文利用ABAQUS有限元軟件進(jìn)行滲流有限元分析,得到了壩體在穩(wěn)定滲流期的各單元節(jié)點(diǎn)水頭,迭代計(jì)算出下游逸出點(diǎn)和浸潤線,并根據(jù)水力梯度矢量計(jì)算出滲透力;然后采用溫控參數(shù)折減有限元法,計(jì)算了新疆克孜爾水庫土石壩壩坡在重力、浮力、滲透力以及地震力作用下的穩(wěn)定性。通過分析,認(rèn)為以壩頂水平位移突變作為結(jié)構(gòu)失穩(wěn)判定準(zhǔn)則較為合理。此外,計(jì)算結(jié)果表明,采用此方法進(jìn)行實(shí)際的工程穩(wěn)定分析是合理的。2流量分析與等載負(fù)荷2.1土骨架產(chǎn)生滲透力土石壩蓄水后,由于上下游存在較大的水頭差,此時(shí),水會在土骨架之間的孔隙中發(fā)生流動(dòng),產(chǎn)生滲流。滲流會對土骨架產(chǎn)生滲透力,這種滲透力與土骨架對水的滲流所產(chǎn)生的阻力構(gòu)成一對作用力與反作用力,對土體能產(chǎn)生一定的破壞作用。在極限平衡法中,地下水的作用是考慮成沿滑裂面上施加的孔隙水壓力;但在限元計(jì)算中,因?yàn)闆]有確定的滑裂面,地下水的作用應(yīng)該通過另一種方式來體現(xiàn),即對浸潤線以下的土體施加浮力和滲透力,這兩種方式是等價(jià)的。2.2單元傳導(dǎo)矩陣調(diào)整法數(shù)值分析在土骨架不可壓縮的穩(wěn)定滲流條件下,由水流的連續(xù)條件和達(dá)西定律可得到土石壩的二維滲流方程式中:kx、kz分別為水平、豎向的滲透系數(shù);h為總水頭。采用有限元法求解無壓滲流時(shí),滲流自由面的確定一直都是一個(gè)難題。由于浸潤線的位置事先是未知的,滲流自由面屬于混合邊界問題,必須同時(shí)滿足水頭邊界條件和流量邊界條件。合理確定滲流自由面,對于分析土石壩邊坡穩(wěn)定具有重要的意義。目前有自由面的滲流計(jì)算,就其計(jì)算方法主要分為變動(dòng)網(wǎng)格法和固定網(wǎng)格法兩種。變動(dòng)網(wǎng)格法為首先假定一個(gè)浸潤線,然后對浸潤線以下區(qū)域劃分網(wǎng)格,給定邊界條件計(jì)算此假定滲流域,接著根據(jù)計(jì)算結(jié)果調(diào)整浸潤線,重新劃分網(wǎng)格,逐步迭代求解,直到浸潤線上滿足h=z的條件為止。固定網(wǎng)格法又分為剩余流量法、單元傳導(dǎo)矩陣調(diào)整法、子單元法和初流量法等。而在目前的滲流分析中比較常用的是先對整個(gè)壩體劃分網(wǎng)格,假定自由面以上區(qū)域內(nèi)的滲透系數(shù)k遠(yuǎn)小于自由面以下區(qū)域內(nèi)的滲透系數(shù)k,然后進(jìn)行試算求解節(jié)點(diǎn)水頭的單元傳導(dǎo)矩陣調(diào)整法。變動(dòng)網(wǎng)格法在進(jìn)行目的為給出滲透力的滲流分析中,要不斷地改變計(jì)算網(wǎng)格,工作復(fù)雜,計(jì)算量大。而在單元傳導(dǎo)矩陣調(diào)整法中物理概念上有不妥之處:滲流不同于明渠流動(dòng),沒有嚴(yán)格意義上的自由面;對浸潤線同時(shí)給定水頭邊界和流量邊界不符合力學(xué)原理(在彈塑性力學(xué)中也不能對一段邊界同時(shí)給定位移和表面力)。本文結(jié)合以上兩種方法的優(yōu)點(diǎn),在土石壩的滲流有限元計(jì)算時(shí),先將整個(gè)壩體劃分網(wǎng)格,假設(shè)下游逸出點(diǎn),計(jì)算出零孔壓線。再根據(jù)零孔壓線與下游壩坡交點(diǎn)位置對逸出點(diǎn)進(jìn)行修正,直到修正的逸出點(diǎn)和計(jì)算出的逸出點(diǎn)位置相同。此時(shí),零孔壓線就是浸潤線。3溫度控制參數(shù)的降壓具有有限元法的基本原則3.1強(qiáng)度儲備安全系數(shù)強(qiáng)度折減法就是在理想彈塑性有限元計(jì)算中,將邊坡土體抗剪強(qiáng)度參數(shù)逐漸降低,直到其達(dá)到極限狀態(tài)為止,同時(shí)得到邊坡的強(qiáng)度儲備安全系數(shù)Fs。于是有式中:c′、?′為土體有效黏聚力和有效內(nèi)摩擦角;c′f、?′f分別為與強(qiáng)度儲備安全系數(shù)sF所對應(yīng)的有效黏聚力和有效內(nèi)摩擦角。有限元強(qiáng)度折減法不需要對滑動(dòng)面形狀和位置做假定,通過強(qiáng)度折減使邊坡達(dá)到極限平衡狀態(tài),此時(shí)所對應(yīng)的強(qiáng)度參數(shù)折減系數(shù)就是邊坡的穩(wěn)定安全系數(shù)sF。3.2土體應(yīng)力準(zhǔn)則同一模型采用不同的屈服準(zhǔn)則,會使計(jì)算結(jié)果有很大的差異。在本文的有限元計(jì)算中采用MohrCoulomb準(zhǔn)則作為屈服準(zhǔn)則。根據(jù)主應(yīng)力和以壓應(yīng)力為負(fù)的規(guī)定,該準(zhǔn)則也可以寫成式中:σ1′、σ3′分別為最大和最小主應(yīng)力。有限元計(jì)算中,采用關(guān)聯(lián)還是非關(guān)聯(lián)流動(dòng)法則,取決于膨脹角ψ值:ψ=?時(shí)(?為土體內(nèi)摩擦角)為關(guān)聯(lián)流動(dòng)法則;ψ=0時(shí)為非關(guān)聯(lián)流動(dòng)法則。本文計(jì)算中采用膨脹角ψ=0的非關(guān)聯(lián)流動(dòng)法則,以便于與圓弧滑裂面極限平衡分析的經(jīng)典算法進(jìn)行對比。3.3溫度場時(shí)變化時(shí)在通常的有限元分析程序里,采用傳統(tǒng)的數(shù)值離散算法時(shí),強(qiáng)度參數(shù)需要在輸入文件中給定,并不隨時(shí)間步而變化,每一組參數(shù)都要進(jìn)行一個(gè)完整的加載過程。顯然,要找到邊坡極限平衡狀態(tài)下的Fs并非易事。然而利用ABAQUS程序中的場變量功能,定義溫度場使材料強(qiáng)度參數(shù)隨著溫度場的變化而變化?？蓪?shí)現(xiàn)直接、方便地修改材料參數(shù),得到安全系數(shù)。此時(shí)溫度場為變量場,并不代表真實(shí)溫度,只起到改變材料參數(shù)的作用。如果給定其熱膨脹系數(shù)為0,那么溫度變化不會造成結(jié)構(gòu)應(yīng)力和應(yīng)變的變化。溫度參數(shù)θ隨時(shí)間t從0線性增長到1,用θ來定義安全系數(shù)為1~10的情況。此時(shí)定義材料參數(shù)隨著θ線性折減,即1/Fs由1折減到0.1。在分析過程中將有效黏聚力和有效內(nèi)摩擦角定義為與溫度參數(shù)相關(guān)的函數(shù)c′f(θ)和?′f(θ)。其數(shù)學(xué)表述關(guān)系式為溫度參數(shù)θ隨時(shí)間增量步t線性變化,有由式(2)、(4)、(5)便可得到安全系數(shù)sF與加載時(shí)間步t之間的關(guān)系這樣,可利用材料強(qiáng)度參數(shù)與加載時(shí)間步t一一對應(yīng),并隨時(shí)間t的增加而線性折減,計(jì)算得到安全系數(shù)sF。3.4壩體邊坡失穩(wěn)狀態(tài)目前,彈塑性有限元方法計(jì)算中對于邊坡的臨界失穩(wěn)或極限平衡狀態(tài)的評判,尚未得到統(tǒng)一的認(rèn)識。主要包括以下幾類:(1)以有限元數(shù)值計(jì)算不收斂作為邊坡失穩(wěn)的標(biāo)準(zhǔn);(2)以邊坡內(nèi)某特征點(diǎn)處的位移突變作為臨界失穩(wěn)的標(biāo)準(zhǔn);(3)以廣義塑性應(yīng)變或者等效塑性應(yīng)變從坡腳到坡頂貫通作為邊坡破壞的標(biāo)準(zhǔn)。本文采用壩坡控制點(diǎn)的水平位移發(fā)生突變時(shí)所對應(yīng)的狀態(tài)作為臨界失穩(wěn)狀態(tài)。此時(shí)與之相對應(yīng)的抗剪強(qiáng)度折減系數(shù)作為土石壩的整體安全系數(shù)。具體來說,根據(jù)計(jì)算可得到壩頂控制點(diǎn)的水平位移和時(shí)間步t的依賴關(guān)系,進(jìn)而根據(jù)這種依賴關(guān)系曲線確定控制點(diǎn)水平位移突變發(fā)生的時(shí)間步t,并將時(shí)間步t代入式(6)計(jì)算折減系數(shù),以此作為壩體安全系數(shù),此時(shí)所對應(yīng)的邊坡狀態(tài)為失穩(wěn)狀態(tài)。在滲流條件下,其采用前述方法計(jì)算滲流對壩體所產(chǎn)生的作用力,并將自重、地震等其他外荷載一起施加于壩體上,采用溫控強(qiáng)度參數(shù)折減有限元法計(jì)算邊坡的安全系數(shù)。考慮滲流作用時(shí),浸潤線以上部分的土體采用天然重度,浸潤線以下部分的土體采用浮重度進(jìn)行計(jì)算。4計(jì)算與分析4.1壩體及壩體材料特性以新疆克孜爾水庫大壩為例進(jìn)行數(shù)值計(jì)算?？俗螤査畮煳挥谛陆S吾爾自治區(qū)阿克蘇地區(qū)拜城縣境內(nèi)渭干河干流木扎提河與支流克孜爾河交匯處,是建造在活動(dòng)斷層帶上的一座大型I等水利工程。大壩為砂礫石壩殼黏土寬心墻壩,壩頂高程為1154.6m,壩頂寬為10.0m,壩頂長為920.6m,最大壩高為44.0m;上游壩坡坡比為1:2.25,在高程1136.0m處設(shè)有寬為2.0m的馬道;下游壩坡坡比為1:2.75,在高程1135.0m處設(shè)有寬為2.0m的馬道。黏土心墻為不對稱寬心墻,其軸線位于壩軸線上游1m處;心墻頂高程為1153.4m,頂寬為6.66m,最大底寬為38.0m,最大高度為42.4m,心墻上、下游側(cè)坡比分別為1:0.6和1:0.2心墻兩側(cè)設(shè)反濾層及反濾過渡層。水庫正常蓄水位1149.6m,下游為1115.0m,滲流計(jì)算考慮壩體為正常運(yùn)行期穩(wěn)定滲流情況。取河床部位最大橫斷面為計(jì)算的典型剖面,進(jìn)行有限元分析。滲流計(jì)算時(shí)基底為不透水邊界,上游水位以下的上游壩坡和下游水位以下的下游壩坡為給定水頭邊界條件,逸出點(diǎn)(其最終位置靠迭代計(jì)算確定)到下游水位之間的下游壩坡為零孔壓邊界條件。穩(wěn)定計(jì)算時(shí)基底采用剛性邊界,兩側(cè)為水平約束,上部自由。大壩典型斷面、材料分區(qū)如圖1所示,材料特性如表1所示。對邊坡穩(wěn)定分析問題,單元類型的選擇不盡相同。Griffiths和Lame采用減縮積分的八節(jié)點(diǎn)四邊形單元,趙尚毅和鄭穎人等采用四節(jié)點(diǎn)四邊形單元和六節(jié)點(diǎn)三角形單元。徐千軍等對各種單元進(jìn)行比較分析,認(rèn)為采用四節(jié)點(diǎn)四邊形單元能夠較好地保證數(shù)值解的準(zhǔn)確性且具有良好的收斂性。本文在邊坡穩(wěn)定計(jì)算中為了與滲流有限元計(jì)算相匹配,采用同一套的網(wǎng)格單元,均為四節(jié)點(diǎn)四邊形單元。4.2結(jié)果與分析4.2.1面逸出點(diǎn)計(jì)算在考慮滲流的邊坡穩(wěn)定分析時(shí),首先要采用迭代法進(jìn)行浸潤線的計(jì)算?？俗螤査畮齑髩蔚湫蛿嗝嬉莩鳇c(diǎn)計(jì)算的迭代過程如表2所示,孔隙水壓力的最終計(jì)算結(jié)果如圖2所示。圖中零孔壓線即為浸潤線,零孔壓線與下游壩坡的交點(diǎn)即為逸出點(diǎn)。從分析過程和計(jì)算結(jié)果不難看出采用此方法計(jì)算浸潤線具有簡化計(jì)算工作、原理可靠合理、可直接利用現(xiàn)成通用有限元軟件分析等優(yōu)點(diǎn)。4.2.2壩體壩體位移與加載時(shí)間步t的關(guān)系首先在不考慮滲流作用的情況下,計(jì)算壩體邊坡的穩(wěn)定安全系數(shù)。此時(shí)各料區(qū)的重度采用天然重度。在此基礎(chǔ)上,分別對考慮滲流作用、聯(lián)合考慮滲流與地震作用兩種情況進(jìn)行計(jì)算?？紤]滲流作用時(shí),對浸潤線以下土體施加浮力和滲透力;考慮地震荷載時(shí),取地震設(shè)計(jì)烈度為Ⅶ度(地震系數(shù)0.1g),在壩體上施加水平向擬靜力。在上述各種條件下進(jìn)行溫控強(qiáng)度折減有限元計(jì)算所得到的安全系數(shù)列于表3。表3中還給出了Bishop法的計(jì)算結(jié)果。由表3可見,上述3種方案溫控強(qiáng)度折減有限元法所確定的安全系數(shù)均高于Bishop法計(jì)算所得到得安全系數(shù),兩者差別在3%~6%之間。這一結(jié)果與張魯渝等認(rèn)為的在一般穩(wěn)定分析中有限元強(qiáng)度折減系數(shù)法所得穩(wěn)定安全系數(shù)比簡化Bishop法平均高出約5.7%的結(jié)果非常相近。單獨(dú)考慮滲流作用時(shí),圖3給出了采用本文方法計(jì)算的下游壩坡坡頂節(jié)點(diǎn)水平位移與加載時(shí)間步t的關(guān)系曲線。圖中水平位移在折減前(分析步1、2)因?yàn)槭艿綁误w自重、自重與滲流同時(shí)作用的影響,壩頂水平位移出現(xiàn)向上游發(fā)展的趨勢;折減開始后(分析步3),由于土體強(qiáng)度參數(shù)的降低,壩頂水平位移明顯地傾向下游。取壩頂水平位移突變時(shí)對應(yīng)的強(qiáng)度折減系數(shù)作為土石壩的整體安全系數(shù),扣除強(qiáng)度折減前的計(jì)算步2后,式(6)中的臨界加載時(shí)間步t=0.392,對應(yīng)的安全系數(shù)sF=1.544。在這種情況下,壩體的位移矢量、變形網(wǎng)格如圖4、5所示。由位移場和變形場的分布情況可知,用這種判定準(zhǔn)則判定邊坡是否達(dá)到極限平衡狀態(tài)是合理的。這種確定安全系數(shù)的方法避免了目前采用數(shù)值計(jì)算不收斂作為判斷標(biāo)準(zhǔn)的人為不確定性。事實(shí)上,本算例在t=0.429時(shí)才開始不收斂。如果按照計(jì)算不收斂作為判斷標(biāo)準(zhǔn),安全系數(shù)sF=1.629,對應(yīng)的水平位移為0.312m,可以認(rèn)為壩坡早己破壞。從圖4、5中可看出,臨界破壞時(shí)壩體的位移與變形發(fā)展情況,同時(shí)也能比較直觀、清晰地評判出臨界破壞面的位置和形狀。對于新疆克孜爾黏土心墻壩,所得到的破壞面接近于圓弧形狀,且滑動(dòng)面的位置比較深。對于不考慮滲流作用和聯(lián)合考慮滲流與地震效應(yīng)時(shí)的壩頂最大水平位移和時(shí)間步t的關(guān)系見圖6、7。從圖6可以看出,在不考慮滲流和地震效應(yīng)時(shí),壩頂?shù)乃轿灰圃趖=0.498時(shí)有明顯的突變,此時(shí)安全系數(shù)Fs=1.811;而計(jì)算在t=0.500時(shí)才不收斂,此時(shí)壩頂?shù)乃轿灰七_(dá)到5.078m,顯然此時(shí)壩坡早已破壞。而在聯(lián)合考慮滲流和地震效應(yīng)時(shí),與單獨(dú)考慮滲流作用時(shí)一樣,在前兩個(gè)分析步中因?yàn)槭艿綁误w自身重力、滲流作用力和地震擬慣性力的作用,在強(qiáng)度折減前壩頂就產(chǎn)生了部分水平位移,且水平位移突變也不如不考慮滲流和地震效應(yīng)時(shí)那樣明顯。但在t=2.152時(shí)(扣除強(qiáng)度折減前的計(jì)算步2后,式(6)中的臨界加載時(shí)間步t=0.152時(shí)),水平位