基于COMSOL的土石壩流-熱耦合影響研究

趙巖松

(大連河海水利水電勘測設計有限公司，遼寧 大連 116021)

1 概述

土石壩滲流引起的穩定性常常導致重大工程問題，造成巨大的損失。如1954年美國提堂水庫潰決造成重大人員財產損失[1]。甘肅某水庫潰壩共造成死亡288人，失蹤40人的重大事故[2]。目前關于土石壩滲流主要依靠達西定律進行計算。宋興亮[3]基于GeoStudio研究了新疆阿勒泰地區土石壩滲流響應問題，結果表明，壩體穩定性隨庫水位的升高而降低。朱趙輝等[4]基于現場監測數據，研究了百色某土石壩滲流特性，詳細分析了滲壓水位變化過程與上游水位相關性。邱德俊等[5]基于ANSYS開展了滲流場與應力變形場計算，結果表明，防滲結構對水位突降或突升均具有良好的防滲效果。劉悅[6]采用GeoStudio研究了深圳某雙面擋水土石壩滲流穩定性。李洪強等[7]基于正交分析原理，研究了土石壩滲流安全監測問題，證明該方法模擬精度較高，為土石壩的穩定性分析提供了新方法。

本文基于COMSOL數值分析軟件，建立了土石壩滲流-溫度耦合數學模型，分析了土石壩在不同滲透工況下滲流場和溫度場的變化過程，本文研究可為相關工程提供經驗。

2 滲流-溫度耦合模型

2.1 滲流模型

本文基于理查茲滲流方法考慮土石壩的滲流過程[8]：

(1)

其中土體函數率采用Van Genuchten模型描述，計算公式為：

θ=θr+Se(θs-θr)

(2)

(3)

式中，α、nv—無量綱參數；hp—壓力水頭，m。

土體容水度和相對滲透率計算表達式為：

(4)

kr=Se1/2[1-(1-Se1/m)m]2

(5)

2.2 溫度模型

溫度場采用熱對流傳熱方程描述[9]，表達式為：

(6)

式中，T—溫度，℃；Cw、Cs—流體和土體的比熱容，J/(kg·°C)；λ—導熱系數，W/(m·K)；DH—彌散系數，m2·s-1；V—水的平均速度，m/s；Qs—熱量源匯項。

3 實例應用

3.1 模型概況

本文研究的土石壩位于大連市。壩體心墻材料為黏土質，頂長134m，頂寬7m，壩高23m，坡率為1∶0.15。土石壩的排水材料主要采用高2m，寬22m的褥墊式水平排水礫石層。為考慮滲流場和溫度場的影響，假定在心墻的6.1、10.1、14.1m處存在寬度為0.02m的裂縫。具體模型如圖1所示。

圖1 土石壩結構圖

3.2 數值模型與參數

3.2.1模型參數

采用工程類比法得到本文模型所使用材料的物理力學參數見表1。其中模型劃分了3個區域，分別代表實際土石壩的壩體、心墻、排水褥墊層。

表1 流體-溫度耦合模型參數

3.2.2邊界條件

本文模型假設定水頭為104m，邊界范圍A-E-D；下游定水頭87m，邊界范圍B-C。其余邊界為零通量邊界。

溫度邊界采用簡化邊界，其中邊界A-E取當地平均溫度10℃，A-B邊界為絕熱邊界，B-C-D-E邊界為大氣邊界，取值為隨時間變化的傅里葉級數：

(7)

式中，T0—平均溫度，℃；k—傅里葉階數；Ak、Bk—系數；L—周期，s。

式7中的系數采用當地氣溫擬合數據，其中T0=16.4，A1=10.3，A2=-0.54，A3=0.27，A4=-0.29，A5=-0.088，B1=-7.1，B2=-0.22，B3=0.3，B4=-0.043，B5=0.03。

4 結果與分析

4.1 模型驗證

為例驗證本文模型的合理性，采用吳志偉等[10]的模型與本文模型進行對比驗證分析，得到不同模型心墻裂縫出口處的溫度和滲透流速變化曲線，如圖2所示。結果表明，采用2種模型得到的結果基本一致，在模型計算的2200d之前，溫度變化基本一致。當發生集中滲漏后，大約40d溫度同時出現最小值，且最小溫度值與庫水溫度基本一致。觀察滲流速度發現，2種模型在2200d以前速度基本為0m/s，之后速度發生突變，且維持在1.3×10-3m/s，采用本文模型速度大小為1.5×10-3m/s，隨后逐漸降低。滲流速度降低的主要原因是由于集中滲漏后心墻上游水頭逐漸降低所導致的。

圖2 集中滲漏時滲透流速及溫度曲線對比

通過與既有模型的比較發現，采用本文模型可以較好的模擬土石壩的滲流-溫度耦合作用下的壩體溫度場和滲流場的動態響應變化過程。

4.2 土石壩溫度場及滲流場變化特征分析

為研究不同位置處集中滲流下土石壩的溫度場和滲流場的變化規律，計算得到了8a之后的集中滲流結果，其中假定土石壩底6.1、10.1、14.1m處形成水平貫穿裂縫，集中滲漏計算時間為300d。

匯總得到不同工況下土石壩的壓力水頭和溫度場變化特征，如圖3—4所示。圖3(a)表明，在無集中滲漏時，浸潤線穿過心墻后快速降低，證明庫水基本可由排水層排至下游；但當心墻發生集中滲漏時(圖3(b)—(d))，心墻裂縫距壩底越高，壩體浸潤線隨著裂縫距壩底越高，到達下游時位置越高。

圖3 土石壩壓力水頭分布

圖4結果表明，在不發生集中滲漏時(圖4(a))，溫度場由心墻庫水低溫區逐漸過渡至下游高溫區，且等溫線變化均勻，發生集中滲漏后(圖4(b)—(d))，溫度等溫線在裂縫處非常集中，并且在低溫區由滲透通道使得壩體溫度下降。綜上分析可知，集中滲漏會嚴重影響土石壩的滲流場和溫度場，具體表現為，不同高度處的集中滲漏對壩體的滲流場和溫度場影響有較大差異。實際工程中，可采用測溫技術確定集中滲漏的位置。

圖4 土石壩溫度場分布

4.3 裂縫處流速及溫度變化特征

匯總得到不同高度裂縫位置的滲流速度對土石壩的滲流場和溫度場的影響強弱對比，如圖5所示。圖5(a)表明，在2920d之前，不同高度裂縫處的滲流速度基本保持0m/s不變，此后，3種工況速度均發生迅速增大，速度達到0.001m/s，隨著時間的增長速度逐漸降低并在300d后保持穩定。因此，心墻裂縫越靠近壩底，速度突變越劇烈且穩定速度也越大。

圖5 不同高度裂縫出口處溫度與流速對比

圖5(b)結果表明，裂縫位置不同處在集中滲漏下，溫度場均發生了突變，在溫度分別為14.5、11.5、11.3℃工況下。對于工況1而言，裂縫出口位置溫度迅速減小至10℃，此后隨氣溫變化而變化，對于工況2和3來說，裂縫出口溫度先增大而后減小，隨后保持10℃。這主要是由于氣溫在土石壩內傳播滯后所導致的。

以上分析表明，當心墻出現裂縫后，下游的滲流場和溫度場均會發生較大的突變，且裂縫高度越大，溫差差異性越大，滲流速度越小，反之，流速變化越大，溫差變化越小。

5 結語

基于COMSOL開展土石壩流體-溫度耦合影響研究，得到以下結論：

(1)與既有模型對比表明，本文提出的模型可以準確的模擬滲流-溫度耦合作用下土石壩滲流場和溫度場的動態變化規律。

(2)集中滲漏可以明顯改變土石壩的滲流場和溫度場，利用此原理可以通過測溫技術準確定位集中滲漏的位置，保證土石壩的安全運行。

(3)心墻裂縫出現位置距壩底越大，集中滲漏對滲流場和溫度場的影響越大。無集中滲漏時，不同高度裂縫處的滲流速度基本保持0m/s不變，發生集中滲漏后不同工況下的速度迅速增大，穩定速度可達0.001m/s。