土石壩心墻接觸面單元參數的標定與數值試驗方法

徐 暢，詹程遠,2，何 吉

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

數值模擬方法已經成為工程設計的重要手段而被廣泛采用，但是該方法在使用時經常由于材料參數的不合理選取而引起模擬結果的顯著偏差，甚至造成嚴重錯誤[1]，因此如何獲取合理的材料參數一直是數值模擬的重點和難點。

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

本文建立了一套標定土石壩心墻接觸面單元參數的數值試驗方法。該方法的主要步驟包括：

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

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

1 接觸面物理試驗

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

圖1 直剪物理試驗示意圖

圖2 接觸面剪切位移與剪切應力的關系曲線

2 數值試驗及參數標定

2.1 本構模型

2.1.1 砂礫石及瀝青混凝土本構模型

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

表1 心墻及過渡料鄧肯-張E-μ模型參數

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

(1)

(2)

(3)

式中：c、φ為抗剪強度參數；Rf為破壞比；Pa為大氣壓力；σ1、σ3為大、小主應力，以壓為正；K、n、F、G、D為無量綱參數。

2.1.2 接觸面本構模型

本文為了模擬瀝青混凝土心墻與砂礫石過渡層之間的復雜力學行為，在兩者間設置了Desai接觸面單元。該單元的應力應變關系服從[12]

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

(4)

(5)

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

(6)

τf=σntanφs+cs

(7)

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

2.2 數值試驗模型

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

圖3 數值試驗模型圖

2.3 標準厚度接觸面參數反演

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

表2 不同厚度的接觸面參數

2.4 增厚接觸面參數標定

為了考慮不同厚度對接觸面參數的影響，將接觸面厚度依次增加至20、50、150 mm。基于不同接觸面厚度的數值模型，分別模擬直剪試驗。計算試件頂面與底面之間(即圖3中cd點之間)的平均水平錯動量，并且繪制其與剪切應力(即單位面積的平均水平推力)的關系曲線(見圖4)。以接觸面厚度為1 mm的曲線為基準，反復調整其他厚度接觸面的各項參數，直至所有厚度接觸面對應的曲線與基準曲線相近為止(見圖4)，此時即為標定所得參數，如表2所示。

圖4 試件剪切位移與剪切應力的關系曲線

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

3 工程應用

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

圖5 壩體材料分布圖

圖6 三維有限元網格圖

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

表4 大壩線彈性模型材料參數

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

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

圖8 心墻上游接觸面錯動量隨高程分布圖

4 結 語

本文針對土石壩瀝青混凝土心墻接觸面單元的參數取值問題，建立了一套數值試驗方法，研究了接觸面參數隨單元厚度的變化規律，并且進行了工程應用。主要結論如下。

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

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

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

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

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