土工膜界面動力剪切特性數值模擬及應用

岑威鈞 ，盛希璇，溫朗昇

(1. 水能資源利用關鍵技術湖南省重點實驗室，湖南 長沙 410014；2. 河海大學 水利水電學院，江蘇 南京 210098)

土工膜是一類柔性薄膜材料，其自身物理力學特性與下墊層土石料差別較大，用于（較）高土石壩壩面防滲時，需要考慮靜動力條件下兩者之間的相互作用效應，即土工膜界面的剪切特性[1-2]。目前國內外已有一些土工膜與土體之間靜力接觸特性的相關研究。例如，Wu 等[3]通過直剪試驗分別對復合土工膜與砂礫石及與聚合物混合碎石接觸面剪切特性進行研究，試驗結果表明前者接觸面的剪應力與剪切位移呈非線性關系，而后者近似呈非線性-理想彈塑性關系。童軍等[4]對復合土工膜-砂礫料墊層界面進行直剪試驗研究，發現界面的摩擦因子與粗顆粒含量、相對密度呈正比。Cen 等[5-8]對土工膜與細砂、砂礫石、土工布、無砂混凝土等工程中常用的墊層材料開展了一系列界面直剪試驗，得到了各類界面的抗剪強度指標，并建立了土工膜-細砂界面靜力接觸面數學模型；同時，對土工膜與砂礫石循環剪切特性進行較深入研究，建立了動力接觸面本構模型。本文簡要介紹土工膜-砂礫石界面循環剪切試驗成果及相應的本構模型，通過在ADINA軟件中的二次開發，對土工膜界面循環剪切試驗進行數值模擬驗證，在此基礎上進一步開展土工膜防滲土石壩的抗震分析應用，得到地震作用下考慮界面動力剪切作用的壩面防滲土工膜的動力變形和應力結果。

1 土工膜界面循環剪切試驗

對室內土石料大型疊環式單剪試驗儀進行改裝，使之適用于土工膜界面循環剪切試驗。試驗采用的材料、試驗方案、試驗控制參數，以及得到的試驗數據、數學模型及參數確定詳見文獻[6-7]。循環剪切試驗設置50、100 和200 kPa 等3 組豎向壓力。剪切位移幅值分別取1.0、2.0、3.0、4.0 和5.0 mm。試驗結果表明：各級豎向壓力下，隨著剪切位移的增加，剪應力逐漸增大，但增大的幅值逐漸減小，呈非線性特性。不同剪切位移幅值下，一次循環剪切形成的滯回圈總體上關于原點呈中心對稱分布。同一豎向壓力下，隨著剪切位移幅值的不斷增加，滯回圈逐漸由狹長變得飽滿，同時剪應力峰值的連線（即主干線）也逐漸向剪切位移軸傾斜，表明接觸面的割線剪切剛度隨著剪切位移幅值增加逐漸減小。相同位移幅值時，豎向壓力越大，剪應力也越大，相應的滯回圈亦越大。

嚴格來說，動力本構建模時應考慮施加荷載的循環性及加載速率等因素。由于較大的加載速率對試驗儀器要求比較高，因此室內也常用循環剪切試驗模擬接觸面的動力特性（常見的地震荷載屬于低頻往復振動），根據低速循環剪切試驗建立的本構模型亦可認為屬于動力本構模型范疇。本文借鑒混凝土與土接觸面模型[9]及土動力學中黏彈性模型建模思路[10-11]，采用等效剪切勁度和阻尼比來刻畫土工膜界面動力剪切特性。根據試驗所得滯回圈及主干線，建立如下土工膜界面循環剪切勁度表達式[7]：

式中：Kmax為最大剪切模量，Kmax=K1γw(σn/Pa)n1，其中 K1和 n1為無因次系數， γw為水的重度， Pa為大氣壓強， σn為法向應力； Rf為破壞比； u為剪切位移； δc為循環剪切摩擦角。

循環剪切試驗得到的滯回圈代表著一個循環加載過程中的能量損失，可用等效阻尼比 λ刻畫這一能量耗損。圖1 給出了不同豎向壓力下土工膜界面等效阻尼比試驗點[7]。由圖1 可見，同一豎向壓力作用下，阻尼比隨剪切位移不斷增長，但增幅變小，最終趨于穩定。

圖 1 土工膜-砂礫料界面阻尼比與剪切位移關系Fig. 1 Relationship between damping ratio and shear displacement of geomembrane-gravel interface

為了考慮上述試驗現象和主要影響因素，建立了如下反映壓力相關的阻尼比數學模型[7]。

式中： λ0和λult分別為初始阻尼比和極限阻尼比； k、 α1和 α2為模型參數；a=1/K1γw(σn/Pa)n1，b=Rf/σntanδc。

根據圖1 中試驗數據，可以求得上述模型的相關試驗參數， λ0、λult、 k、 α1和 α2分別為0.32、0.72、1.25、1.84 和0.012 4[7]。圖1 中給出了由式（2）得到的擬合曲線，與試驗點吻合較好，驗證了所建阻尼比模型的模擬能力。

2 界面動力剪切特性的有限元模擬

對上述土工膜-砂礫料界面循環剪切試驗進行有限元數值模擬，計算模型尺寸同實際試驗值，即上盒中為直徑300.0 mm、高度180.0 mm 的砂礫石圓柱形試樣；下剪切盒為長方體，其邊長360.0 mm×360.0 mm，高度120.0 mm。在土工膜和砂礫料之間設置薄層接觸單元，厚度為5.0 mm。切向剪切勁度采用式（1），阻尼比采用式（2）計算。上剪切盒內土石料試樣周邊設置三向固定約束，下剪切盒為剛性材料，采用彈性模量設置大值進行模擬。對下剪切盒施加往復位移荷載，具體工況設置同模型試驗，即往復位移幅值分別取1.0、2.0、3.0、4.0 和5.0 mm。

圖2 給出了計算模型接觸面中心節點處豎向壓力50 kPa 和200 kPa 下的剪應力和剪切位移關系曲線。由圖2 可見，不同壓力和剪切位移幅值下各滯回圈總體上與相應的試驗點吻合較好，主干線自然也吻合較好，初步驗證了所建土工膜-砂礫料界面動力接觸本構模型的合理性。

圖 2 土工膜-砂礫石界面的剪應力和剪切位移關系Fig. 2 Relationship between shear stress and shear displacement of geomembrane-gravel interface

3 土工膜界面動力本構模型應用分析

3.1 計算模型與計算條件

某土工膜防滲堆石壩，壩高80 m，壩頂寬度9.3 m，上游壩坡1∶1.5，下游壩坡1∶1.8，壩面采用0.75 mm厚HDPE 土工膜防滲，土工膜上鋪塊石保護層，下設砂礫料墊層，計算水位78 m。大壩三維有限元計算模型節點13 006 個，單元11 704 個。在土工膜與墊層之間設置了382 個薄層接觸面單元。開展大壩地震反應計算之前先要進行靜力分析，其中前19 步模擬地基、壩體分級填筑、鋪設壩面土工膜和保護層，以及壩面施加水荷載等靜力加載過程，第20 步施加地震荷載。動力計算時，壩基采用無質量地基，地震輸入采用Taft 波，將順河向和壩軸向峰值加速度調為0.1g，豎向峰值加速度調為0.067g。靜力計算時壩體材料采用鄧肯E-B 模型，動力計算時采用Hardin-Drnevich 模型，相應的計算參數見表1。靜動力計算時土工膜與墊層之間采用相同接觸面本構模型。

表 1 壩料靜動力計算參數Tab. 1 Static and dynamic parameters of dam materials

3.2 結果分析

計算得到的壩體靜、動力應力和變形結果符合土石壩一般規律，限于篇幅不再詳細展示。由于接觸面模型的設置主要影響土工膜的受力變形特性，因此對土工膜的相關計算結果進行重點分析。圖3 為蓄水期壩面土工膜變形等值線分布圖。由圖3 可見，對稱河谷條件下壩面土工膜變形呈左右對稱分布，順河向、豎向和總變形均沿壩高呈先增大后減小的趨勢，最大值出現在河床壩段的40%壩高附近。壩軸向變形由兩岸指向河床。進一步對比墊層上表面位移，發現土工膜與墊層之間發生了相對位移，說明兩者之間出現了剪切滑移。由此可見，土工膜界面是否設置接觸面模型對計算結果有一定影響，主要減小了土工膜的剪切位移。

圖 3 蓄水期土工膜位移分布（單位：cm）Fig. 3 Displacement distribution of geomembrane during water storage period (unit: cm)

圖4 為蓄水期壩面土工膜主拉應變和主拉應力等值線分布。由圖4 可見，土工膜主拉應變分布與主拉應力分布規律基本相似，最大值均位于兩側岸坡偏底部，表明受邊界約束土工膜的可能受拉破壞主要集中于此。計算表明，蓄水期土工膜主拉應變和主拉應力極值分別為0.71%和0.80 MPa，遠小于0.75 mm 厚HDPE 土工膜的屈服伸長率14.7%和屈服強度12.78 MPa，說明土工膜不會拉壞。與不設置土工膜接觸面的結果相比，土工膜主拉應變和主拉應力極值均有所減小，這與前述的土工膜與墊層之間發生了相對剪切位移是互相吻合的。

圖 4 蓄水期土工膜靜力主拉應變和主拉應力Fig. 4 Static principal tensile strain and stress of geomembrane during water storage period

圖5 為地震作用時壩面土工膜的動位移等值線包絡圖（即絕對值的最大值）。由圖5 可見，土工膜三向動位移沿壩高逐漸增大，在河床中央附近壩頂處達到最大值，與壩體動位移極值處相吻合。圖6 為地震作用時壩面土工膜動主拉應變和動主拉應力等值線包絡圖。地震作用下土工膜主拉應變和主拉應力均呈現隨壩高的增加先增大后減小的分布規律，最大值出現在80%壩高附近。由于在土工膜與墊層之間設置了接觸面單元，地震作用時發生了相對位移，因此土工膜動主拉應變和動主拉應力值均較小，其值分別為0.05%和0.05 MPa。這與前述土工膜與墊層之間發生了剪切位移有關。設置接觸面模型后，計算得到的土工膜動拉應力和拉應變均有所減小。進一步將蓄水期和地震作用下壩面土工膜的靜動主拉應力進行疊加，可得到總的土工膜主拉應力包絡圖（圖7）。由圖7 可知，靜動疊加后土工膜主拉應力分布規律與靜力結果相似，最大值出現在岸坡偏底部，其值為0.81 MPa，小于土工膜的屈服拉伸強度。由此可見，對于土工膜防滲土石壩來說，無論靜力狀態還是地震作用時，應重點關注岸坡和河床底部錨固處土工膜的受力特性和工作狀態，防止該區域土工膜發生局部受拉破壞。

圖 5 地震作用下壩面土工膜動位移包絡圖（單位：cm）Fig. 5 Dynamic displacement of geomembrane under earthquake (unit: cm)

圖 6 地震作用下壩面土工膜動主拉應力和動主拉應變包絡圖Fig. 6 Dynamic principal tensile stress and strain of geomembrane under earthquake

圖 7 靜動力疊加后土工膜主拉應力包絡圖（單位：MPa）Fig. 7 Principal tensile stress of geomembrane after static and dynamic superposition (unit: MPa)

4 結 語

（1）土工膜-砂礫料界面循環剪切試驗表明，剪切勁度合理反映界面動力剪切的非線性特性，等效阻尼比反映界面動力剪切的滯后性。循環剪切試驗的有限元數值模擬結果與試驗值總體吻合較好，驗證了土工膜界面動力剪切模型的合理性。

（2）土工膜防滲土石壩三維靜動力有限元計算表明，蓄水期土工膜變形分布規律和極值均在合理范圍。地震作用下，土工膜動位移、加速度、動拉應變和動拉應力分布規律性好。靜動疊加后的壩面土工膜主拉應力分布大值區主要集中在岸坡偏底部和河床底部相錨固附近。由于土工膜與墊層之間設置了接觸面單元，土工與墊層之間發生了剪切位移，因此土工膜的主拉應力總體較小，極值為0.81 MPa，遠小于其抗拉屈服強度，土工膜安全可以保障。