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

黃 寧,何 蕃 民,吳 夢 喜

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

0 引 言

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

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

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

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

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

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

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

1 剪切帶變形模擬方法

1.1 防滲墻頂部土體剪切帶產(chǎn)生機理

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

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

1.2 剪切面單元

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

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

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

(1)

(2)

(3)

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

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

(4)

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

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

(5)

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

1.3 整體坐標(biāo)系下剪切面單元剛度矩陣

剪切面單元節(jié)點力、節(jié)點位移在整體坐標(biāo)與局部坐標(biāo)系下滿足下列關(guān)系[35]:

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

(6)

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

(7)

式中,旋轉(zhuǎn)矩陣[α]寫成

(8)

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

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

(9)

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

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

(10)

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

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

2 涵洞算例

2.1 計算模型及參數(shù)

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

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

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

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

2.2 涵洞周圍土體位移

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

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

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

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

2.3 涵洞頂部土壓力

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

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

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

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

3 典型土質(zhì)心墻壩

3.1 計算模型及參數(shù)

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

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

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

(11)

(12)

(13)

式中:K,n,Rf,G,F,D為試驗參數(shù);Rf為剪應(yīng)力破壞比;c為黏聚力;φ為內(nèi)摩擦角;pa為大氣壓力。

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

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

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

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

3.2 壩體位移

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

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

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

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

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

3.3 防滲墻應(yīng)力

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

表5 100 m級心墻壩防滲墻應(yīng)力監(jiān)測值[32]Tab.5 Stress monitoring value of cutoff wall of 100 m core dams

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

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

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

4 結(jié) 論

針對結(jié)構(gòu)物頂部土體內(nèi)部剪切帶,如土質(zhì)心墻壩混凝土防滲墻頂部土體在壩體填筑過程中因差異沉降產(chǎn)生的剪切帶,本文提出了一種在可能出現(xiàn)剪切帶的部位預(yù)設(shè)剪切面單元的方法來模擬,得出以下結(jié)論:

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

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

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