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

趙巖松

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

1 概述

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

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

2 滲流-溫度耦合模型

2.1 滲流模型

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

(1)

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

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

(2)

(3)

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

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

(4)

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

(5)

2.2 溫度模型

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

(6)

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

3 實例應(yīng)用

3.1 模型概況

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

圖1 土石壩結(jié)構(gòu)圖

3.2 數(shù)值模型與參數(shù)

3.2.1模型參數(shù)

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

表1 流體-溫度耦合模型參數(shù)

3.2.2邊界條件

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

溫度邊界采用簡化邊界，其中邊界A-E取當(dāng)?shù)仄骄鶞囟?0℃，A-B邊界為絕熱邊界，B-C-D-E邊界為大氣邊界，取值為隨時間變化的傅里葉級數(shù)：

(7)

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

式7中的系數(shù)采用當(dāng)?shù)貧鉁財M合數(shù)據(jù)，其中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 結(jié)果與分析

4.1 模型驗證

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

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

通過與既有模型的比較發(fā)現(xiàn)，采用本文模型可以較好的模擬土石壩的滲流-溫度耦合作用下的壩體溫度場和滲流場的動態(tài)響應(yīng)變化過程。

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

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

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

圖3 土石壩壓力水頭分布

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

圖4 土石壩溫度場分布

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

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

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

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

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

5 結(jié)語

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

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

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

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