王 羿，王正中，劉銓鴻，劉 月

基于彈性薄層接觸模型研究襯砌渠道雙膜防凍脹布設

王 羿，王正中※，劉銓鴻，劉 月

（西北農(nóng)林科技大學旱區(qū)寒區(qū)水工程安全研究中心；旱區(qū)農(nóng)業(yè)水土工程教育部重點實驗室，楊凌，712100）

渠道襯砌下雙膜防滲層可解除渠基土凍結力，釋放襯砌結構位移以適應不均勻凍脹，但會影響襯砌結構的剛度與穩(wěn)定性。為研究最佳雙膜防滲層布設形式及其適應性，考慮雙膜防滲層接觸本構及襯砌結構與渠基土的相互作用，依據(jù)基土水熱力三場耦合理論建立了寒區(qū)襯砌渠道的凍脹模型，通過數(shù)值模擬仿真分析了不同斷面及不同雙膜防滲層布設形式的弧底梯形渠道襯砌的凍脹受力與變位。結果表明：寬淺斷面渠道采用雙膜襯砌結構時，膜間摩擦力過小造成剛度和穩(wěn)定性嚴重降低以及襯砌拉應力增加。適當?shù)哪らg摩擦力可以解除部分凍結約束，調(diào)整局部不均勻凍脹，而且發(fā)揮弧底反拱作用不產(chǎn)生拉應力。采用PE膜與無紡布作為層間接觸的雙膜襯砌結構，凍脹位移方差減小25%，同時整體位移增加不超過0.2 cm，削減凍脹應力50%以上。在不顯著降低結構整體穩(wěn)定性前提下，其整體剛度和結構強度提高程度最大，故同時適用于窄深渠道和寬淺渠道。研究為寒區(qū)襯砌渠道中雙膜防滲層的合理布設提供科學依據(jù)。

渠道；膜；防凍；凍脹破壞；雙膜防滲層；水熱力耦合；凍脹模型；接觸本構

0 引 言

中國寒旱區(qū)存在大范圍季節(jié)凍土，其周期性凍融造成土體表面不均勻變形，對上部建筑物產(chǎn)生破壞，導致新建渠道運行2 a后必須進行翻修，造成了嚴重的經(jīng)濟損失和輸水功能的喪失[1]。為此，開展寒區(qū)渠道防凍融破壞措施研究，對于保障工程正常輸水功能、降低造價具有重要意義。

寒區(qū)輸水渠道凍融破壞是土體凍融變形[2]以及土體襯砌之間相互作用力共同作用結果，這一作用力通過接觸面?zhèn)鬟f。早期研究在分析襯砌破壞時，將凍土與襯砌結構相互作用看作凍土對襯砌結構的法向與切向的凍脹力，而法向凍脹力產(chǎn)生的根源在于襯砌板與渠基凍土凍結為一體的切向凍結力的存在[3-4]，其值與凍土凍結溫度、含水率和襯砌材料結構相關，通?？捎稍囼灉y定[2,5-7]，基于此建立了襯砌渠道凍脹破壞的簡化結構力學模型。隨后的研究，將凍土與襯砌視為2種剛度不同的材料組成的一體化整體復合結構，按襯砌與凍土協(xié)調(diào)變形建立了有限元模型對襯砌變形與應力進行計算，與結構力學簡化模型相比有限元法可以考慮到凍土與結構協(xié)調(diào)變形過程中的相互作用隨時間與空間的變化[8-10]。以上研究均表明，引起渠道襯砌板破壞的主要外力是襯砌板下的切向凍結力及由此產(chǎn)生的法向凍脹力，而切向凍結力產(chǎn)生的必要條件是渠道邊坡襯砌板和其下凍土層之間凍結約束的存在[11]?；诖?，劉旭東等[12]提出在傳統(tǒng)單膜防滲的基礎上采用雙膜結構，以解除切向凍結、釋放切向位移（雙膜即雙層塑料薄膜），進而達到減小凍脹破壞的目的，并用ADINA有限元軟件對雙膜襯砌渠道的凍脹變形及凍脹力進行了數(shù)值模擬，提出了“采用雙膜防滲復合襯砌渠道解除凍結約束來削減襯砌板凍脹受力、均勻化凍脹分布來減輕凍脹破壞”的建議。隨后，李爽等[13-15]研究表明凍土與襯砌間存在接觸層，并通過設計試驗裝置對接觸層的力學性質(zhì)進行了研究。然而，過度解除基土凍結約束對于大型渠道及預制板襯砌都存在結構剛度與穩(wěn)定性不足的嚴重問題，特別是開春融化期[16-17]；而且雙膜襯砌渠道的凍脹研究尚存在以下3 個問題：1）對渠基土凍脹的計算沒有考慮水熱力三場耦合凍脹過程及原位水及遷移水凍結相變的熱效應和變形效應[18-19]；2）沒有考慮雙膜間摩擦力取值的影響，沒有準確探明不同雙膜布設形式摩擦約束不同對襯砌凍脹力與變形的影響規(guī)律[12]；3）沒有研究過度解除基土凍結約束對于大型渠道及預制板襯砌結構剛度與穩(wěn)定性的影響，也未探明斷面形狀對雙膜襯砌適用性的問題。針對以上問題，本文采用多物理場耦合軟件COMSOL Multiphysics建立考慮溫度-水分-變形3個物理場耦合的橫觀各向同性土體凍脹本構模型；建立考慮了單膜、無摩擦雙膜及實際工程中的3種雙膜布設形式等5種情況下的凍土-膜-襯砌結構的多界面接觸模型，研究弧底梯形斷面與雙膜結構相應的凍脹規(guī)律與防凍脹效果。

1 雙膜防滲層襯砌渠道凍脹模型

1.1 襯砌渠道凍脹破壞特點及雙膜防滲層防凍脹機理

渠道基土凍脹時，渠槽上部內(nèi)縮而底部上抬，混凝土襯砌板與基土發(fā)生相對位移從而產(chǎn)生凍脹力和凍結力[1]。其襯砌破壞過程可以概況為：1）由于沿渠道斷面基土的凍脹量分布不均，凍脹量大的弧底基土對襯砌產(chǎn)生法向凍脹力，導致弧底襯砌上抬同時對臨近的渠坡襯砌產(chǎn)生頂托作用。2）渠坡襯砌在渠底襯砌頂托下發(fā)生相對基土向上的切向位移趨勢，當這種趨勢受到凍結約束時產(chǎn)生了切向凍脹力，同時反作用于渠底襯砌約束其上抬位移。3）渠頂基土自由凍脹，與其凍結在一起的襯砌隨之受到拉力。4）綜合凍脹力與凍結力作用結果，渠底弧形段和臨近渠坡段為偏心受壓狀態(tài)，渠頂襯砌與臨近渠坡襯砌為偏心受拉狀態(tài)，呈現(xiàn)渠底抬升和渠坡鼓脹開裂的破壞特點[2,7,20-21]。

根據(jù)以上分析，基土的凍結約束限制了襯砌位移，降低襯砌對基土凍脹變形的適應性，導致結構下部呈現(xiàn)壓彎組合變形和頂部呈現(xiàn)拉彎組合變形，不利于襯砌結構的安全。通過在襯砌下鋪設2層防滲土工膜，利用膜間相對滑動消除凍結約束[12]，從而消減或調(diào)整不均勻法向凍脹力從而達到防止凍脹破壞的目的。工程上雙層土工膜的鋪設形式有3種（圖1）：1）雙層一布一膜土工膜采用膜-膜接觸；2）雙層一布一膜土工膜采用膜-布接觸；3）雙層兩布一膜土工膜采用布-布接觸，如圖1。3 種鋪設方式關鍵區(qū)別在于2層土工膜間摩擦系數(shù)不同，對凍結約束的解除程度也不同。文獻[12]探討了膜間完全無摩擦力情況下理想雙膜措施的凍脹變形適應能力。但是完全無摩擦力是理想狀態(tài)既不符合工程實際，而且完全無凍結約束也會對單薄的襯砌結構的剛度及穩(wěn)定性產(chǎn)生不利，影響結構安全。

1.2 雙膜防滲層襯砌渠道凍脹數(shù)值模型

從渠基凍土凍脹的基本理論出發(fā)[22-24]，首先建立考慮溫度、土質(zhì)、含水率和凍結過程等因素的土體凍脹模型，其次建立合理反映凍土-雙膜-襯砌之間相互作用的接觸本構模型，最后結合2個模型計算渠基土體凍脹量及襯砌結構對凍脹的力學響應。

1.2.1 渠基凍土凍脹基本理論與水熱力三場耦合模型

考慮到渠道為線狀工程，順水流方向（縱向）長度遠大于垂直水流方向（橫向），且縱向土質(zhì)、含水率變化較小，可以近似看作平面應變問題處理。二維平面渠基土熱傳導及熱量平衡方程[25]如下。

式中T為土體溫度，K；ρ為土體密度，kg/m3；Cp為土體等效定壓熱容，J/(kg·K)；λeq為土體等效導熱系數(shù)，W/(m2·K)；Q為熱源項，W；x、y為二維平面坐標軸，m；t為時間，s。土體密度、等效定壓熱容與等效導熱系數(shù)的定義見文獻[26]。

式（1）熱源發(fā)生于正凍土區(qū)，由水分遷移相變所產(chǎn)生，按文獻[18,27]其定義如下

在凍土溫度場-水分場耦合控制方程基礎上，為描述土體在溫度梯度方向上的凍脹[28]，在文獻[10,18]的基礎上，對凍脹本構進行修正，得到渠基凍土平面應變狀態(tài)下力學控制方程

式中σ和σ分別為凍土內(nèi)和方向上的正應力分量，MPa；τ為切應力分量，MPa；f和f為外力，N；ε0和ε0分別為凍土內(nèi)由于凍脹引起的、方向線應變分量，m/m；ε0為凍脹引起的土體剪應變分量，m/m；grandT為溫度梯度方向的單位厚度土體凍脹量，m/m；和為溫度梯度向量的方向余弦；為土體彈性模量，MPa；為土體泊松比。由式（1）～式（5）共同組成描述凍土凍脹的溫度-水分-應力（變形）三場耦合模型的控制方程。采用多物理場分析軟件Comsol Multiphysics求解耦合的多物理場控制方程組，從而可以數(shù)值模擬已知環(huán)境溫度、含水率和土質(zhì)條件下的混凝土襯砌渠道的凍脹過程。模型準確性驗證已在文獻[18]中詳細闡述，本文不再贅述，而是在原模型基礎上對凍脹量方向進行了改進，以表現(xiàn)土體在溫度梯度方向的凍脹行為。

1.2.2 渠道襯砌-雙膜-凍土相互作用接觸模型

土工膜厚約為0.25～0.5 mm，在計算時處理為凍土與襯砌結構間的接觸界面，采用界面法向反力和切向反力來表示兩者接觸作用。普通單膜情況下凍脹時土工膜上下表面分別與襯砌板底和渠基土體凍結在一起，將基土凍脹力和凍結力全部傳遞給襯砌[29]。雙膜情況下，上、下土工膜分別與襯砌和土體凍結，但雙膜之間按接觸界面考慮，可以分離也可以傳遞法向凍脹力，不能傳遞法向凍結力，而通過摩擦力傳遞切向凍結力約束。據(jù)此，本文通過對彈性薄層模型進行改進，用于描述雙膜情況下襯砌與凍土的接觸力學作用。

彈性薄層模型即將結構間的接觸層處理為具有一定剛度的法向和切向彈簧單元，根據(jù)相互接觸結構間的相對位移來計算兩者間的接觸反力，并根據(jù)反力調(diào)整結構間相對位移。與接觸界面算法相比，將動態(tài)的搜索過程簡化為靜態(tài)的彈塑性平衡問題，尤其將動摩擦力處理為切向彈簧剛度，使程序計算容易收斂且保持一定的真實性。模型表達式如下：