土質心墻壩防滲墻頂部土體剪切帶模擬研究

黃 寧,何 蕃 民,吳 夢 喜

(1.四川交通職業技術學院 道路與橋梁工程系,四川 成都 611130; 2.中冶成都勘察研究總院有限公司,四川 成都 610023; 3.中國科學院力學研究所,北京 100190; 4.中國科學院大學,北京 100049)

0 引 言

建造在深厚覆蓋層上的土質心墻壩,壩基中常設置混凝土墻防滲。防滲墻與其上部壩體的土質心墻多采用插入式或廊道式連接[1]。剛性、薄而深的防滲墻夾在柔性覆蓋層土體之中,合理地確定混凝土防滲墻應力是大壩防滲墻設計的重要內容,也是一個難題[2]。

防滲墻在垂直方向上的受力主要包括:作用在防滲墻或廊道頂部的土壓力,墻身因兩側土體與墻體差異沉降而產生的側壁摩擦力以及墻身的自重。從垂直方向看,防滲墻是偏心受壓結構[2]。防滲墻在水平方向的作用力主要是墻兩個側面的水壓力和土壓力。由于防滲墻一般嵌入巖槽一定深度,因此其周邊嵌入基巖的墻身水平方向位移受到較強約束。從水平方向看,防滲墻屬于周邊受約束的薄板,墻身的彎矩對防滲墻的應力影響很大[2-3]。土對防滲墻的作用力(包括防滲墻上的土壓力與摩擦力)引起防滲墻變形,防滲墻變形又影響作用力的大小。墻和土是相互作用的一個體系,合理的計算方法必須將地基、防滲墻和壩體視為整體,用變形協調將三者聯系起來[3]。覆蓋層和壩體填筑土體的變形模擬,以及土與結構相互作用的模擬,兩者均對防滲墻應力的合理計算產生重要影響。

土與結構相互作用的模擬,包括土與結構接觸面上的不連續變形(剪切滑移和接觸面張開)和結構物上部土體中剪切帶的模擬。對于土與結構物的接觸面,通過在土與結構物界面上設置接觸面單元,如Goodman接觸面單元[5]和Desai接觸面單元[6],可以較好地模擬出土體與結構接觸面上的非連續變形特性。對于土體內部剪切帶的模擬,則起源于邊坡穩定的有限元法研究[7]。土體內部剪切帶的萌生與擴展是導致邊坡滑動失穩的原因。20世紀60年代,Skepmton[8]和Bjerrum[9]對超固結黏土邊坡漸進破壞問題,最先開啟了剪切帶擴展理論的研究。在此基礎上Palmer[10],Rice[11-12],Puzrin等[13]基于斷裂力學方法推導了各種情況下臺階狀邊坡中預設滑裂面的擴展條件,而土體中剪切帶的萌生條件可由Mohr-Coulomb強度理論描述已經成為共識[14]。

早在1980年,Johnson等[15]就意識到,使用連續插值函數的傳統有限元方法(位移法)并不是特別適合于塑性問題的求解,因此建議修改形函數使非連續位移可以被更好描述。在四邊形單元中嵌入非連續線來模擬剪切帶的方法[16-17]被提出。越來越多的研究采用在單元內嵌入不連續位移場的處理方法,允許單元產生較大畸變,以便模擬整體的強不連續位移場[18-21]。大量的研究比較了在單元中嵌入非連續位移的各種方法[22-25],這些方法用于模擬因剪切帶形成的不連續位移場在網格依賴性和應力鎖死方面較彌散裂紋模型有改進。然而,這些方法對于不連續位移場的描述均基于網格單元,剪切帶尖端難以確定,也就無法描述和揭示由土體應變軟化特性導致的尖端應力集中和重分布對剪切帶擴展過程和方向的影響[25]。

擴展有限元等方法[26-29]可將物體的材料非線性和土體結構的邊界非線性區別對待,既可較好地描述內部運動的不連續位移場,又可相對獨立地選擇合適的本構模型,為剪切帶擴展過程的模擬提供了強有力的技術支持。與單元內嵌入非連續方法相比,擴展有限元法網格依賴性較小,對剪切帶形狀和界面接觸狀態的描述更準確。這類方法有可能考察剪切帶尖端應力集中對擴展過程和方向的影響,促進了已有剪切帶擴展分析方法,并給剪切帶擴展過程的模擬研究提供了新的思路。

土中剪切帶擴展機理和模擬方法等方面仍有很多問題需要深入研究。由于涉及大變形、應變局部化[30]、非連續[31]等力學研究的前沿理論,相關的數值模擬方法并不成熟,也難以應用于土體非線性彈性數值計算體系,同時還存在如網格的依賴性[28]等問題,將這些數值模擬方法直接應用于實際工程還比較困難。

對于深厚覆蓋層上的土質心墻壩中結構物頂部在壩體填筑過程中因差異沉降產生的土體內部剪切帶,本文提出一種在可能出現剪切帶的部位預設剪切面單元的方法,來探討模擬土中剪切帶及其擴展。此方法首先用于溝埋式涵洞頂部土壓力離心機試驗算例[30],驗證方法的合理性后,再給出了用于模擬土質心墻壩中結構物頂部土體內非連續變形的算例。根據大渡河、瀑布溝等深厚覆蓋層上的高心墻堆石壩應力變形計算與實際監測[32-33],常規數值分析計算防滲墻應力明顯被高估,因此有必要發展防滲墻頂部土體剪切帶的模擬技術,提高其計算分析精度。

1 剪切帶變形模擬方法

1.1 防滲墻頂部土體剪切帶產生機理

防滲墻的變形模量比覆蓋層和壩體填筑土體大幾個數量級,因而防滲墻身的垂直壓縮量遠遠小于兩側的覆蓋層。如圖1所示,在大壩填筑過程中,在上部填土的自重作用下覆蓋層和防滲墻之間的沉降存在較大的差異,墻土之間會產生剪切滑移;當填筑體達到防滲墻的頂部后,其上部土體的填筑將防滲墻頂部覆蓋起來,隨著壩體填筑高程的繼續上升,防滲墻頂部兩側土體的沉降量將顯著大于墻頂部正上方土體的沉降量(防滲墻或廊道的頂托作用),即防滲墻上部土體中會產生沉降量急劇變化的剪切帶。底部設置剛性涵洞結構物填土的離心機試驗表明[25,30],結構頂部土體內的剪切帶(滑移帶)是真實存在的,滑移帶兩側的差異沉降在結構頂部最大,向上逐步減小,并在上部出現0差異沉降的等沉面[4]。通過上述分析可見,土質心墻壩防滲墻頂部土體剪切帶相對確定,發生的起始地方一般在防滲墻頂部,但剪切帶不一定貫穿結構物上部的全部土體,即滑移帶的結束位置是未知的。

圖1 防滲墻頂部土體的變形示意Fig.1 Soil deformation on top of cut-off wall

1.2 剪切面單元

如圖2所示,在等厚度薄層四邊形單元(本文以二維為例)中,坐標系XOY為整體坐標,坐標系X′O′Y′為局部坐標,坐標軸X′方向為剪切帶土體滑移方向(簡稱“切向”),坐標軸Y′方向為法向。

圖2 剪切面單元(局部坐標與整體坐標)Fig.2 Shear band element(local coordinate and global coordinate)

在局部坐標系下,當前計算步i的剪切面單元應力應變關系與普通的實體單元在形式上基本一致,具體表達式如下(省略了計算步下標i):

(1)

(2)

(3)

根據有限單元法基本理論[34],局部坐標系下剪切面單元節點力{F′}e與節點位移{δ′}e之間的平衡關系可以寫成:

{F′}e=[K′]e{δ′}e

(4)

式中,[K′]e為剪切面單元的剛度矩陣,采用式(5)計算。

[K′]e=?Ω[B′]T[D′][B′]dΩ

(5)

式中:彈性矩陣[D′]=[C′]-1,[B′]為平面應變矩陣。

1.3 整體坐標系下剪切面單元剛度矩陣

剪切面單元節點力、節點位移在整體坐標與局部坐標系下滿足下列關系[35]:

{F′}e=[R]{F}e, {δ′}e=[R]{δ}e

(6)

式中:{F}e,{δ}e分別為整體坐標系下剪切面單元節點力和節點位移;{F′}e,{δ′}e分別為局部坐標系下剪切面單元節點力和節點位移。對于圖2所示的剪切面單元(四節點矩形單元),單元旋轉矩陣[R]為

(7)

式中,旋轉矩陣[α]寫成

(8)

將式(6)代入式(4),整理后為

{F}e=[K]e{δ}e=[R]-1[K′]e[R]{δ}e

(9)

根據式(9),可得出整體坐標系下剪切面單元剛度矩陣如下:

[K]e=[R]-1[K′]e[R]

(10)

整個剪切面單元剛度矩陣計算流程如圖3所示。

圖3 剪切面單元剛度矩陣計算流程Fig.3 Calculation flow of stiffness matrix of shear band element

2 涵洞算例

2.1 計算模型及參數

以溝埋式涵洞土壓力室內試驗[30]進行模擬分析,并將模擬結果與試驗結果進行對比。如圖4(a)所示,矩形溝槽寬B為0.5 m,高為1.1 m;涵洞高h、寬D均為0.1 m。模型填料為中細砂,容重14.0 kN/m3,內摩擦角31°,黏聚力為0,壓縮模量0.5 MPa,泊松比0.21。模型有限元網格劃分如圖4(b)所示。為探討剪切面單元的設置高度,初始時在涵洞頂部土體內沿涵洞兩側豎直向上設置兩個貫穿這個砂土層的剪切面單元,單元厚度為1 cm。剪切面單元彈性模量E、泊松比v與填料取值一樣,而剪切模量G按式(2)、(3)計算。

圖4 矩形溝埋式涵洞Fig.4 Rectangular trench buried culvert

在細砂與涵洞、木板的交界面上設置無厚度的Goodman接觸面單元[5],接觸面單元的本構關系采用Duncan-Clough雙曲線模型[31],計算參數見表1。

表1 涵洞算例接觸面計算參數Tab.1 Calculation parameters of contact element of culvert calculation example

2.2 涵洞周圍土體位移

涵洞頂部土體內不設和設置剪切面單元計算得到土體的豎向位移如圖5所示。由圖5可見,兩種方法計算得到的最大豎向位移都發生在土層的中部,最大值分別為7.33 mm和7.40 mm。設置剪切面單元計算的沉降最大值更大,且得到的涵洞頂部土體等值線明顯更加密集,剪切面單元(圖中紅色虛線)兩側錯動位移明顯,錯動位移達到2 mm。而不設置剪切面單元的計算模型則無法反映出土體內部剪切帶的錯動變形。

圖5 土體的豎直方向位移(單位:mm)Fig.5 Vertical displacement of soil

繪制涵洞頂部剪切面單元兩側土體沉降差與土體高度關系曲線,如圖6所示。由圖6可知,越靠近涵洞頂部沉降差最大(沉降差為2.16 mm),隨著土體高度增加,沉降差逐漸減小,當土體高度為0.275 m時,沉降差降低至0.01 mm以內,真實再現了試驗中觀測到的等沉面[4]。由此可見,預設剪切面單元具備模擬剪切破壞前土體連續的變形特性,以及土體剪切破壞發生后剪切帶土體的錯動變形特性,且兩種變形狀態模擬轉換在計算過程中自動判別。

圖6 涵洞頂部土體沉降差Fig.6 The subsidence difference of soil on top of culvert

2.3 涵洞頂部土壓力

填土高度為1.0 m時涵洞頂部豎向土壓力分布情況如圖7所示?？梢钥闯?設置與不設置剪切面單元計算得到的涵洞頂部土壓力值平均分別為15.48,16.36 kPa,分別較試驗值(15.35 kPa)大0.8%,6.2%。

圖7 涵洞頂部豎向土壓力Fig.7 Vertical earth pressure on top of culvert

從土壓力的分布來看,在涵洞頂部中間兩側各3.5 cm范圍兩種計算方案得出的土壓力值基本重合,而越靠近涵洞左右兩端,土壓力數值差值越大,說明應力集中效應更加明顯。將涵洞頂部左右端土壓力值列于表2。由表2可以看出,不設置剪切面單元計算結果要比試驗值更大,最大高出35.0%,由此可以看出僅僅設置土與結構接觸面單元是不夠的;設置剪切面單元計算結果更接近試驗值,計算得出土壓力值比試驗值,最大高出僅為2.4%。由此可見,設置剪切面單元計算結果更接近實際值。

表2 涵洞頂部土壓力值比較Tab.2 Comparison of earth pressure on top of culvert kPa

3 典型土質心墻壩

3.1 計算模型及參數

以一個典型土質心墻壩作為工程算例(如圖8所示)。大壩壩高100 m,覆蓋層厚度為50 m。壩基覆蓋層防滲采用厚1.4 m的混凝土防滲墻,墻底部嵌入基巖1 m,防滲墻上部與心墻的連接方式采用插入式連接,插入深度15 m;其中墻頂部5 m為變截面段,截面由寬1.4 m逐漸變為寬0.5 m;以下10 m及覆蓋層中防滲墻厚度均為1.4 m。

圖8 典型土質心墻壩結構(尺寸單位:m)Fig.8 Structure drawing of typical core rockfill dam

計算模型采用平面應變模型,混凝土防滲墻和基巖均采用線彈性模型,混凝土防滲墻彈性模量取30 GPa,泊松比取0.17,基巖彈性模量取20 GPa,泊松比取0.28;筑壩材料采用鄧肯E-v模型[36],切線變形模量Et和切線泊松比vt分別按式(11)、(12)計算,具體計算參數見表 3。

(11)

(12)

(13)

式中:K,n,Rf,G,F,D為試驗參數;Rf為剪應力破壞比;c為黏聚力;φ為內摩擦角;pa為大氣壓力。

表3 心墻壩計算參數Tab.3 Calculation parameters of core rockfill dam

表4 心墻壩接觸面計算參數Tab.4 Calculation parameters of contact element of core rockfill dam

在土體與混凝土防滲墻之間設置無厚度的Goodman接觸面單元,單元的本構關系采用Duncan-Clough雙曲線模型,計算參數見表 4。實際工程中,覆蓋層中防滲墻周邊裹挾著一層泥皮,而插入心墻的防滲墻周圍是高塑性黏土,因此覆蓋層中防滲墻接觸面剪切參數中黏聚力值較心墻中的要更小些。

計算中不考慮水荷載。下文對比不設剪切面單元和設置剪切面單元兩種方案的計算結果。

3.2 壩體位移

兩種方案計算的壩體整體位移規律相同,壩體上下游水平位移基本對稱,上下游壩體水平位移最大值分別出現在靠近堆石區中部建基面區域,最大值分別約為-0.20 m(向上游)和0.20 m(向下游);壩體最大沉降出現在心墻下部2/3部位,為1.22 m。防滲墻頂部設置剪切面單元時計算得到的壩體填筑完成后位移等值線如圖9所示。

圖9 壩體位移(單位:m)Fig.9 The displacement of dam

防滲墻頂部周圍土體豎向位移計算結果、墻頂部土體非連續位移與高程關系分別如圖10、11所示?？梢钥闯?防滲墻頂部設置剪切面單元后,防滲墻兩側土體的豎向位移明顯增大,使得墻頂兩側土體內出現高約10 m的剪切錯動帶;墻頂最大錯動位移約為0.4 m。由此可見,防滲墻頂部非連續變形的模擬與否直接影響著兩側土體的變形計算結果。此外,圖10顯示了插入心墻土體部分周邊土體沉降等值線較下部砂卵石層有更明顯的集中,其因有二:① 防滲墻插入心墻周圍設置有高塑性黏土,而高塑性黏土的黏聚力比覆蓋層中的砂卵石層更大;② 防滲墻上部的變截面也導致局部土體沉降變形受限。

圖10 防滲墻頂部周圍土體豎向位移(單位:m)Fig.10 Vertical displacement of soil around the top of cut-off wall

圖11 防滲墻頂部土體中非連續位移Fig.11 Discontinuous displacement of soil on top of cut-off wall

3.3 防滲墻應力

防滲墻等截面段應力沿高程分布情況如圖12所示。由圖12可見,兩種方案計算得到防滲墻的水平向正應力基本相同。對于豎直向正應力,兩種計算方案均呈現出應力值隨覆蓋層深度增加而增大的趨勢,應力值在8.0～22.0 MPa之間(見表5),與文獻給出的部分100 m級心墻壩防滲墻應力監測值基本吻合[32]。從設置剪切面單元計算結果來看,設置剪切面單元的計算結果明顯比不設時的計算結果小,其中墻體中部減小最大,最大減小了10.4%。

表5 100 m級心墻壩防滲墻應力監測值[32]Tab.5 Stress monitoring value of cutoff wall of 100 m core dams

圖12 防滲墻應力沿高程分布Fig.12 Stress distribution along elevation of cutoff wall

圖13給出了防滲墻變截面段頂部豎向應力的結果。設置剪切面單元后,豎向應力平均值8.30 MPa,與不設時豎向應力平均值(11.93 MPa)相比,設置后計算結果小30.43%左右?？梢?在土與結構相互作用時,若不考慮對結構物頂部土體的非連續變形進行模擬,結構的應力計算結果可能會與實際相差較大。

圖13 防滲墻頂部豎向應力Fig.13 Vertical stress at top of cut-off wall

4 結 論

針對結構物頂部土體內部剪切帶,如土質心墻壩混凝土防滲墻頂部土體在壩體填筑過程中因差異沉降產生的剪切帶,本文提出了一種在可能出現剪切帶的部位預設剪切面單元的方法來模擬,得出以下結論:

(1) 在填土1 m的矩形涵洞算例中,設置了剪切面單元后,計算得到涵洞頂部土體荷載結果分布更接近試驗值。由此說明,本文提出的預設剪切面單元,既可以模擬普通土體單元連續變形特性,也可以模擬剪切破壞時剪切帶土體的錯動變形特性。

(2) 典型土質心墻壩計算結果表明,設置剪切面單元后,防滲墻兩側土體內出現高約10 m的錯動變形區域,其中墻頂最大錯動位移約為0.4 m。從防滲墻頂部豎向應力計算結果來看,設置剪切面單元計算得到豎向應力平均值約8.30 MPa,較不設時降低約30%?？梢?在土與結構相互作用時,若不對結構物頂部土體的非連續變形進行模擬,結構的應力計算結果可能會與實際相差較大。

(3) 矩形涵洞算例和典型土質心墻壩工程算例的計算分析表明,在進行土與結構相互作用分析時,除應考慮土與結構接觸界面上的非連續變形外,還應充分考慮結構物頂部土體的非連續變形對結構物應力變形特性的影響。本文提出的在可能出現剪切帶位置預設剪切面單元的分析方法,克服了當前剪切帶數值模擬中繁瑣的計算流程以及對計算網格依賴性強等問題,且實現過程簡單、高效,具有廣闊的應用前景。