土石壩心墻接觸面單元參數(shù)的標(biāo)定與數(shù)值試驗(yàn)方法

徐 暢，詹程遠(yuǎn),2，何 吉

(1. 武漢大學(xué)水利水電學(xué)院，湖北 武漢 430072； 2. 長(zhǎng)江科學(xué)院，湖北 武漢 430010)

數(shù)值模擬方法已經(jīng)成為工程設(shè)計(jì)的重要手段而被廣泛采用，但是該方法在使用時(shí)經(jīng)常由于材料參數(shù)的不合理選取而引起模擬結(jié)果的顯著偏差，甚至造成嚴(yán)重錯(cuò)誤[1]，因此如何獲取合理的材料參數(shù)一直是數(shù)值模擬的重點(diǎn)和難點(diǎn)。

對(duì)于瀝青混凝土心墻壩而言，壩體主要由堆石料壩殼、瀝青混凝土心墻、過(guò)渡層(位于心墻兩側(cè))等部分組成[2]。這些結(jié)構(gòu)通?；谟邢迒卧?，采用非線性彈性模型(鄧肯-張模型等)[3]或者彈塑性模型(K-G模型等)[4]等本構(gòu)模型進(jìn)行模擬，相應(yīng)的材料參數(shù)可以通過(guò)成熟的試驗(yàn)方法獲得。但是由于心墻與過(guò)渡層的力學(xué)性能存在顯著差異，因此現(xiàn)實(shí)中兩者間存在剪切、錯(cuò)動(dòng)等相對(duì)變形，并且伴隨局部侵入等復(fù)雜現(xiàn)象[5]。這些復(fù)雜力學(xué)行為的合理模擬對(duì)于正確判斷心墻的應(yīng)力狀態(tài)，避免防滲體的破壞至關(guān)重要。為了模擬心墻和過(guò)渡層之間復(fù)雜的力學(xué)行為，兩者間通常需要設(shè)置接觸面單元[6]。常用的接觸面單元分為無(wú)厚度單元和有厚度單元兩種，前者以Goodman單元[7]為代表，后者以Desai單元[8]為代表。Goodman單元在受壓時(shí)接觸面兩側(cè)存在相互嵌入的問題。為了減少嵌入的程度，通常需要顯著增大接觸面的法向剛度，但是這樣極易造成迭代計(jì)算時(shí)誤差的明顯積累[9-10]。土石壩的數(shù)值模擬無(wú)論在施工期、蓄水期，還是運(yùn)行期均涉及大量迭代計(jì)算，因此更適合采用Desai單元。

本文建立了一套標(biāo)定土石壩心墻接觸面單元參數(shù)的數(shù)值試驗(yàn)方法。該方法的主要步驟包括：

1)采用鄧肯-張非線性彈性模型和Desai接觸面單元，建立土石壩心墻接觸面直剪試驗(yàn)的數(shù)值模型；

3)基于數(shù)值試驗(yàn)，標(biāo)定增厚的接觸面單元參數(shù)，本文依托茅坪溪土石壩的物理試驗(yàn)成果，采用上述方法，標(biāo)定了心墻接觸面單元在不同厚度時(shí)的參數(shù)，研究了各參數(shù)隨單元厚度的變化規(guī)律，并將標(biāo)定所得參數(shù)成功應(yīng)用于某土石壩工程的三維數(shù)值模擬中。

1 接觸面物理試驗(yàn)

茅坪溪土石壩是三峽水利樞紐的重要組成部分。為了研究該壩瀝青混凝土心墻與砂礫石過(guò)渡層之間接觸面的力學(xué)性能，張治軍[11]等人采用大型疊環(huán)單剪儀進(jìn)行了試驗(yàn)。試驗(yàn)設(shè)備如圖1所示，圖1中下部為方形剪切盒，尺寸600 mm×600 mm×300 mm(長(zhǎng)×寬×高)，盒內(nèi)裝有砂礫石材料；上部為堆疊放置的10層方環(huán)，每環(huán)尺寸600 mm×600 mm×30 mm，底部?jī)蓪臃江h(huán)內(nèi)裝有瀝青混凝土材料，其余8層方環(huán)內(nèi)裝有砂礫石材料。試驗(yàn)中保持法向應(yīng)力1.5 MPa不變，逐步增大水平推力直至試件破壞。試驗(yàn)結(jié)果顯示，破壞面沿上接觸面或下接觸面隨機(jī)分布，接觸面剪切位移與剪切應(yīng)力的關(guān)系曲線如圖2所示。

圖1 直剪物理試驗(yàn)示意圖

圖2 接觸面剪切位移與剪切應(yīng)力的關(guān)系曲線

2 數(shù)值試驗(yàn)及參數(shù)標(biāo)定

2.1 本構(gòu)模型

2.1.1 砂礫石及瀝青混凝土本構(gòu)模型

根據(jù)專著《長(zhǎng)江三峽水利樞紐建筑物設(shè)計(jì)及施工技術(shù)》[12]，茅坪溪土石壩在大部分研究中瀝青混凝土心墻、砂礫石過(guò)渡料以及壩體各種堆石料均推薦采用鄧肯-張E-μ模型，并且提供了具體參數(shù)，如表1所示。本文在數(shù)值試驗(yàn)中繼承了上述研究成果。

表1 心墻及過(guò)渡料鄧肯-張E-μ模型參數(shù)

鄧肯-張E-μ模型為非線性彈性模型，相應(yīng)的加載切線變形模量Et、卸載切線變形模量Eur、切線泊松比μt分別為[12]：

(1)

(2)

(3)

式中：c、φ為抗剪強(qiáng)度參數(shù)；Rf為破壞比；Pa為大氣壓力；σ1、σ3為大、小主應(yīng)力，以壓為正；K、n、F、G、D為無(wú)量綱參數(shù)。

2.1.2 接觸面本構(gòu)模型

本文為了模擬瀝青混凝土心墻與砂礫石過(guò)渡層之間的復(fù)雜力學(xué)行為，在兩者間設(shè)置了Desai接觸面單元。該單元的應(yīng)力應(yīng)變關(guān)系服從[12]

{dε}=[C]{dσ}

(4)

(5)

式中：[C]為接觸面的柔度矩陣；Et、μt為接觸面的切線彈性模量和切線泊松比，假定與砂礫石過(guò)渡料相同(式(1)～式(3))；Gt為接觸面的切線剪切模量，其表達(dá)式為

(6)

τf=σntanφs+cs

(7)

式中：τf為接觸面的抗剪強(qiáng)度；σn為接觸面的法向正應(yīng)力；φs、cs為接觸面的摩擦角和黏聚力；Rfs為接觸面的破壞比；Gs、ns為無(wú)量綱參數(shù)。

2.2 數(shù)值試驗(yàn)?zāi)Ｐ?/h3>

根據(jù)圖1中物理試驗(yàn)的試件尺寸和材料分布，建立三維數(shù)值模型(見圖3)，同時(shí)施加約束和荷載如下：物理試驗(yàn)中方形剪切盒的約束作用，以法向約束的方式施加于數(shù)值模型下部砂礫石材料的所有外表面；物理試驗(yàn)的法向應(yīng)力和千斤頂反力，分別以面力的形式沿鉛直方向和水平方向施加于數(shù)值模型的頂面。

圖3 數(shù)值試驗(yàn)?zāi)Ｐ蛨D

2.3 標(biāo)準(zhǔn)厚度接觸面參數(shù)反演

假定一組接觸面參數(shù)(包括cs、φs、Rfs、Gs、ns)，基于數(shù)值模型(見圖3)，采用有限元法模擬直剪試驗(yàn)。計(jì)算接觸面邊界之間(即圖3中ab點(diǎn)之間)的平均水平錯(cuò)動(dòng)量，并且繪制其與剪切應(yīng)力(即單位面積的平均水平推力)的關(guān)系曲線(見圖2)。采用參數(shù)反演方法，反復(fù)調(diào)整接觸面的各項(xiàng)參數(shù)，直至該曲線與物理試驗(yàn)的結(jié)果相近為止(見圖2)，此時(shí)接觸面的參數(shù)即為反演所得參數(shù)，如表2中第一行所示。

表2 不同厚度的接觸面參數(shù)

2.4 增厚接觸面參數(shù)標(biāo)定

為了考慮不同厚度對(duì)接觸面參數(shù)的影響，將接觸面厚度依次增加至20、50、150 mm?；诓煌佑|面厚度的數(shù)值模型，分別模擬直剪試驗(yàn)。計(jì)算試件頂面與底面之間(即圖3中cd點(diǎn)之間)的平均水平錯(cuò)動(dòng)量，并且繪制其與剪切應(yīng)力(即單位面積的平均水平推力)的關(guān)系曲線(見圖4)。以接觸面厚度為1 mm的曲線為基準(zhǔn)，反復(fù)調(diào)整其他厚度接觸面的各項(xiàng)參數(shù)，直至所有厚度接觸面對(duì)應(yīng)的曲線與基準(zhǔn)曲線相近為止(見圖4)，此時(shí)即為標(biāo)定所得參數(shù)，如表2所示。

圖4 試件剪切位移與剪切應(yīng)力的關(guān)系曲線

由表2可知，各參數(shù)與接觸面厚度的相關(guān)性中，Gs最大，cs、φs及Rfs較小，ns基本無(wú)關(guān)；cs、φs及Gs與厚度正相關(guān)，Rfs為負(fù)相關(guān)。由于增厚以后的接觸面單元不僅反應(yīng)了接觸面的力學(xué)性能，也反應(yīng)了部分過(guò)渡層材料的性能，因此隨著厚度增加，過(guò)渡層的影響逐漸增大，考慮到過(guò)渡層的力學(xué)性能強(qiáng)于接觸面，造成接觸面單元的力學(xué)性能逐漸增強(qiáng)。表2中各參數(shù)隨厚度的變化規(guī)律符合上述原理。

3 工程應(yīng)用

由于無(wú)法獲取茅坪溪土石壩詳盡的設(shè)計(jì)資料，不能為其建立三維數(shù)值模型，本文選取某瀝青混凝土心墻壩作為工程應(yīng)用。該壩壩高64.5 m，壩軸線長(zhǎng)586 m，心墻厚0.4 m，過(guò)渡層厚2 m。茅坪溪與該壩的材料分區(qū)并非完全相同，為了避免不同分區(qū)的干擾，專注于接觸面的研究，因此假定壩體僅采用一種堆石料，同時(shí)假定壩基僅含有一種基巖材料，材料分區(qū)如圖5所示。

圖5 壩體材料分布圖

圖6 三維有限元網(wǎng)格圖

表3 大壩鄧肯-張E-μ模型材料參數(shù)

表4 大壩線彈性模型材料參數(shù)

基于三維有限元網(wǎng)格(見圖6)，采用有限單元法模擬了壩體的分層施工過(guò)程，總計(jì)16個(gè)計(jì)算步。施工結(jié)束時(shí)壩體的沉降分布如圖7所示，由圖7可知該壩施工期最大沉降量為591.55 mm，位于壩體中部。由于最大沉降量小于壩高1%，所以符合土石壩設(shè)計(jì)規(guī)范[13]的要求。心墻上游接觸面與過(guò)渡層之間的錯(cuò)動(dòng)量如圖8所示，由圖8可知心墻相對(duì)過(guò)渡層整體向下錯(cuò)動(dòng)，即心墻的沉降量大于過(guò)渡層，壩體上部錯(cuò)動(dòng)量較大、下部較小，最大錯(cuò)動(dòng)量在1.8 mm以內(nèi)。為了比較錯(cuò)誤的接觸面參數(shù)對(duì)錯(cuò)動(dòng)量的影響，此處選取表2中1 mm厚的參數(shù)進(jìn)行計(jì)算，得到心墻與過(guò)渡層間的錯(cuò)動(dòng)量如圖8所示。由圖8可知，采用1 mm厚的參數(shù)時(shí)會(huì)明顯夸大層間錯(cuò)動(dòng)量。這是因?yàn)镚s與厚度正相關(guān)，厚度越小時(shí)接觸面越容易發(fā)生變形。因此增厚以后的接觸面單元不能沿用增厚之前的參數(shù)，工程中需要根據(jù)接觸面單元的實(shí)際厚度標(biāo)定參數(shù)，否則會(huì)對(duì)結(jié)果產(chǎn)生顯著影響。

圖7 施工期最大壩高斷面沉降分布圖

圖8 心墻上游接觸面錯(cuò)動(dòng)量隨高程分布圖

4 結(jié) 語(yǔ)

本文針對(duì)土石壩瀝青混凝土心墻接觸面單元的參數(shù)取值問題，建立了一套數(shù)值試驗(yàn)方法，研究了接觸面參數(shù)隨單元厚度的變化規(guī)律，并且進(jìn)行了工程應(yīng)用。主要結(jié)論如下。

1)為了模擬瀝青混凝土心墻與砂礫石過(guò)渡層之間復(fù)雜的力學(xué)行為(剪切、錯(cuò)動(dòng)、局部嵌入等)，推薦在兩者間設(shè)置Desai接觸面單元。該單元是一種常用的有厚度接觸面單元，使用方法簡(jiǎn)便。

2)在土石壩的數(shù)值建模中，特別是三維建模，為了降低建模難度、減少單元總數(shù)，通常需要選擇較厚的接觸面單元設(shè)置于心墻與過(guò)渡層之間。接觸面單元在選取參數(shù)時(shí)必須考慮厚度的影響。

3)本文基于有限元方法和直剪物理試驗(yàn)，建立了一套標(biāo)定不同厚度接觸面單元參數(shù)的數(shù)值試驗(yàn)方法。依托茅坪溪土石壩的物理試驗(yàn)成果，驗(yàn)證了該方法的可行性。

4)本文采用上述數(shù)值試驗(yàn)方法，依托茅坪溪土石壩的物理試驗(yàn)成果，研究了接觸面單元各項(xiàng)參數(shù)與厚度的相關(guān)關(guān)系。研究結(jié)果顯示，相關(guān)性中Gs最大，cs、φs及Rfs較小，ns基本無(wú)關(guān)；cs、φs及Gs與厚度正相關(guān)，Rfs為負(fù)相關(guān)。究其原因，接觸面單元不僅反應(yīng)了接觸面的力學(xué)性能，也反應(yīng)了部分過(guò)渡層材料的性能，因此隨著接觸面單元厚度的增加，過(guò)渡層的影響也逐漸增大，考慮到過(guò)渡層的力學(xué)性能強(qiáng)于接觸面，造成接觸面單元的力學(xué)性能逐漸增強(qiáng)。

5)本文將數(shù)值試驗(yàn)標(biāo)定的接觸面參數(shù)，應(yīng)用于某瀝青混凝土心墻壩的三維數(shù)值模型中，成功模擬了心墻和過(guò)渡層之間的相互錯(cuò)動(dòng)。研究結(jié)果表明，接觸面參數(shù)對(duì)心墻與過(guò)渡層之間的錯(cuò)動(dòng)量具有顯著影響，進(jìn)一步論證了工程中需要根據(jù)接觸面單元的實(shí)際厚度標(biāo)定參數(shù)的重要意義。