基于Ren模型的土石坝集中渗漏流热耦合研究

倪 枫 任 杰,2 张文兵 王 帆

（1.西安理工大学 省部共建西北旱区生态水利国家重点实验室,西安 710048；2.水利部堤防安全与病害防治工程技术研究中心,郑州 450003；3.河海大学 水利水电学院,南京 210098）

1 研究背景

土石坝渗漏一直以来都是土石坝隐患诊断的重点问题,对于土石坝的渗漏监测主要方法有同位素示踪法、电阻率成像法、渗流热监测技术以及弹性波探测等,这几种方法各自有优缺点[1].同位素示踪虽然可以准确地示踪出土石坝渗漏轨迹,但是布置造价过高,难以实现;电阻率成像法等方法缺乏试验支撑,并且图像难以解译,如果应用到实际工程中比较困难.相比之下温度作为天然的示踪剂可以更方便地应用到实际工程中,因此渗流热监测技术是目前研究土石坝渗漏的一个重要手段.渗流热监测技术的理论基础是岩土体中渗流场与温度场的耦合作用[2],该技术的原理是介质中的水分运移控制着土石坝温度场的分布,坝体温度场的时空分布反过来能够反映坝体渗流场的特征.

土壤导热系数是土壤热力学特性的一个重要指标,是影响土壤传热过程中温度分布和能量变化的主要因素,也是研究其物理过程,如流-热耦合、物质传递、气体扩散等的基础[3-4].Lu等[5]根据前人的研究,在2007年将土壤大致划分为粗土和细土,并引入了和砂粒土含量相关的经验常数.2015年,Lu和Dong等[6]提出了闭合土体的热导率经验方程,该方程将土壤水分特征曲线与导热系数结合起来,可以较准确地预测不同类型土体的导热系数.考虑到有机质对热导率的影响和土壤颗粒组分对形状因子的影响,Ren等[7]于2019年提出一个新的经验模型,与之前的热传导系数模型相比,Ren（2019）模型充分反映了土体有机质成分和颗粒成分对形状因子的影响.

本文基于COMSOL Multiphysics软件建立Ren（2019）模型的土石坝流-热耦合数学模型,以典型二维心墙土石坝的数值模拟作为算例,将模拟结果与吴志伟等[8]基于Lu（2007）模型的土石坝流热耦合模型进行对比,在此基础上模拟了土石坝心墙在不同集中渗漏情况下温度场和渗流场的变化,研究结果可为新的导热系数在渗流热监测技术的应用提供参考.

2 数学模型

2.1 渗流场模型

渗流场采用理查兹方程表示[9-10]:

式中:ρw为水的密度;g为重力加速度;Se为土体相对饱和度;Ss为弹性贮水率;p为压强;θ为含水量;ks为介质饱和渗透系数;kr（θ）为非饱和带相对渗透率[0≤kr（θ）≤1];μ（T）为水的动力黏度[10];z代表需要计算位置的高程;Cm为土体容水度,代表水分特征曲线斜率的负倒数;Qm代表水流源汇项;∇代表拉普拉斯算子.

用van Genuchten模型描述非饱和带土壤水分特征曲线[10]:

式中:θr为残余含水率;θs为饱和含水率.其中Se为:

式中:Se为土壤相对饱和度;α和nv为VG 模型参数;hp为压力水头;α是水分特征曲线进气值的倒数;nv是分特征曲线坡度的指示参数,m＝1－1/nv,在非饱和带压力水头hp等于负压水头hc[11].

土体容水度Cm和相对渗透率kr分别采用下述经验公式表示:

2.2 温度场模型

温度场采用热对流传热方程描述[12]:

式中:T为温度;Cw、Cs分别为流体和土体的比热容;λ为有效导热系数;DH为水动力弥散系数;u为平均速度;Qs为热量源汇项.

COMSOL中自带计算传热模块是用热量运移方程进行计算的,用PDE 模块来代替COMSOL 自带的传热模型进行计算,可灵活修改导热系数经验模型.

COMSOL中提供的系数型PDE方程如下[14]:

式中:ea、da、c、α、γ、β、a、q和h均为自定义系数;Ω为求解域,Ω的偏微分是其外边界;n是外法线方向;f、g、r分别为求解域和边界上的源项.

2.3 Ren（2019）土体导热系数模型

Ren等[7]考虑有机质含量对导热系数的影响及土壤粒组成分对形状因子的影响,提出一个新的导热系数经验模型:

式中:λsat和λdry与砂土含量和干土容重（γd）有关,当砂含量0＜Csand＜1 时,干土容重11＜γd＜20.γd（k N/m3）和ρb（g/cm3）之间的关系式为:

式中:g＝9.8 m2·s－1,是重力加速度.

Ke与θ的新关系如下:

式中:α和β是λ（θ）的形状因子,它们的取值会影响导热系数曲线的斜率和增长速度,表达式如下:

式中:Com为有机质的质量比（g·kg－1）;Cclay、Csilt、Csand分别表示黏粒、粉粒、砂粒的质量分数（%）.

3 实例应用

3.1 研究区域概况

选取湖北省黄冈市大山背水库进行土石坝的集中渗漏研究[15].大山背土石坝工程的基本概况如下:坝体为黏土心墙代料组合坝,坝顶长134 m、宽7 m;坝顶高程108.12 m,坝高23.2 m;心墙顶高程106.32 m,顶宽2.0 m,边坡为1∶0.15.坝体填筑材料假设为粉质黏土,心墙材料为黏土,坝体排水体采用褥垫式水平排水,排水褥垫层高2 m,宽22 m,褥垫层材料为砾石.为了说明不同渗漏方式对坝体温度场和渗流场的影响,本文假设在某一特定时刻,心墙形成0.02 m 宽的裂缝,并假设心墙裂缝分别存在3个不同高度（6.10、10.10、14.10 m）,如图1所示.

图1 二维心墙土石坝结构示意图

3.2 模型边界与具体参数

3.2.1 边界条件

渗流场方面,坝体闪游定水头为104.92 m,边界为A-E-D,坝体下游定水头86.92 m,边界为B-C,其余边界均为零通量边界.

温度场方面,由于湖北省属于华中地区且本研究对象坝高在23.2 m,根据文献[8,16]的数据分析和论证,库水温度随高程变化基本可忽略,因此简化边界A-E温度按当地水库的平均水温取10℃;A-B作为绝热边界;B-C-D-E作为大气边界,温度取日平均气温.当地温度随时间周期性变化的傅里叶级数函数[14]如下:

式中:T0为平均温度;k为傅里叶级数的阶数;Ak和Bk为傅里叶级数的系数;L为周期.

这里采用湖北省黄冈市的气温数据（数据来源于中国气象数据网站）并根据傅里叶函数来取5阶傅里叶级数就可以满足精度要求.拟合得气温回归模型系数为:T0＝16.37、A1＝－10.28、A2＝－0.541 7、A3＝0.266 7、A4＝－0.291 7、A5＝－0.088 22、B1＝－7.133、B2＝－0.216 5、B3＝0.3、B4＝－0.043 3、B5＝0.083 44.将系数代入式（15）即可求得温度随时间的周期性规律.

3.2.2 模型参数

结合当地坝体材料并参考文献[15],给出模型渗流场及温度场的计算参数,见表1.区域1、2、3分别代表图1中的坝体、心墙、排水褥垫层,3个区域的材料主体分别是粉质黏土、黏土、砾石,心墙裂缝区域的饱和水力传导率参照排水褥垫层取10－2m/s.

表1 心墙土石坝流-热耦合模型计算参数

4 结果与讨论

4.1 模型对比分析

为了开展基于Ren（2019）模型的土石坝流-热耦合研究,有必要验证该模型的可靠性和准确性,因此参照吴志伟等[8]进行的基于Lu模型土石坝集中渗漏研究,在与之模型边界条件和参数都相同的情况下模拟基于Ren（2019）模型的土石坝集中渗漏研究.基本工况为水库从0 d时刻开始蓄水,在0～2 190 d之间处于蓄水之后的正常工作期,随后在第2 190 d时心墙裂缝在距坝底10.08 m 处开始发生集中渗漏,并持续200 d.因此设置第一个计算区段为2 190 d,之后发生集中渗漏,计算时间为200 d,取时间步长1 d.通过模拟计算出心墙裂缝出口处的温度和渗透流速变化,将结果与吴志伟等[8]计算出的结果进行对比如图2所示.

图2 心墙裂缝出口处渗透流速及温度曲线对比图

从图2可以看出,基于两个不同导热系数模型计算出的心墙裂缝出口处温度和渗透流速变化趋势接近一致,两个模型在第2 190 d之前的心墙裂缝出口处温度变化趋势基本一致;发生集中渗漏以后,两个模型心墙裂缝出口处温度变化基本一致,在经过40 d左右温度同时达到最低,最低温度与库水温度基本一致,之后受大气温度影响逐渐上升;裂缝出口处渗透流速,两个模型在第2 190 d之前基本都接近于0 m/s,在第2 190 d时突变后的渗透流速同属于一个数量级10－3m/s.吴志伟等[8]计算出的心墙裂缝出口处渗透流速在第2 190 d突增后一直维持在1.33×10－3m/s,基于Ren模型计算出的心墙裂缝出口处渗透流速在第2190 d突增到1.48×10－3m/s,之后200 d逐渐降低稳定到1.31×10－3m/s,渗透流速在突增之后逐渐降低的原因是发生集中渗漏之后心墙上游面水头逐渐降低,后渗透流速达到稳定.

通过对上述两个模型计算结果的综合分析,得出基于Ren（2019）模型的土石坝流-热耦合模型可以准确地反映土石坝温度场和渗流场的变化过程.

4.2 不同集中渗漏条件下土石坝温度场及渗流场的变化特征

为研究不同集中渗漏条件下土石坝温度场及渗流场的变化规律及特征,建立基于Ren（2019）模型的土石坝饱和-非饱和流-热耦合模型.设置土石坝在正常运行8a（2920 d）之后发生集中渗漏,假设土石坝心墙分别在距坝底6.08、10.08、14.08 m 形成水平贯穿裂缝,对应工况1、2、3,集中渗漏计算时间为300 d,时间步长1 d.

各工况以及无集中渗漏的坝体渗流场如图3所示,图中红线代表坝体的浸润线.由图3（a）可知,浸润线在穿过心墙后急剧降低,之后基本与排水褥垫层持平,说明此时坝体的渗流库水基本都能从排水褥垫层排到下游,有效地保证了土石坝的安全运行.由图3（b）、（c）、（d）与图3（a）对比可知,当心墙发生集中渗漏时,浸润线在穿过心墙下游表面时会随着心墙裂缝的高度变化而发生变化,心墙裂缝距坝底越高,坝体浸润线在穿过心墙下游表面时的位置越高.各工况以及无集中渗漏的坝体温度场如图4所示.在没有集中渗漏的情况下,坝体温度场是从心墙前库水的低温区逐步向下游的高温区过渡,等温线分布比较均匀,而在集中渗漏之后,心墙前库水通过图4（b）、（c）和（d）所示的不同高度的渗漏通道分别进入心墙下游非饱和带,等温线在心墙裂缝进口处较为密集,且低温库水渗过集中渗漏通道后使心墙后坝体温度降低,相应的大气温度对坝体温度影响的区域减少.

图3 不同时刻坝体压力水头分布图（单位:m）

图4 不同集中渗漏条件下坝体温度场的变化（单位:℃）

比较上述4种情况的土石坝温度场和渗流场,得出坝体温度场分布在一定程度上受集中渗漏的影响,而且不同高度的集中渗漏通道对坝体温度场的影响不同,可以通过测温来反演渗流场去确定集中渗漏通道的具体位置.

4.3 心墙裂缝出口处渗透流速及温度变化规律分析

图5（a）为不同高度的心墙裂缝出口处渗透流速变化对比图.由图5可以看出,不同高度的心墙裂缝出口处渗透流速在2 920 d之前一直保持接近0 m/s,在2 920 d时3种工况出口处渗透流速发生突变,且都达到10－3m/s以上,之后随着时间推移逐渐降低并在300d后基本全部达到稳定.从图5可以看出对不同高度心墙裂缝突变后的渗透流速大小不同,且心墙裂缝越靠近坝底,渗透流速突变越大且达到稳定之后的渗透流速也越大.

图5（b）为不同高度的心墙裂缝出口处温度变化对比图.由图5（b）可以看出,不同高度的心墙裂缝出口处温度在发生集中渗漏的瞬间都发生了温度突降.通过分析图4中无渗漏情况的坝体温度场,可以得出工况1、2、3在集中渗漏发生前的温度分别为14.4℃、11.5℃、11.3℃.当集中渗漏在3种工况分别发生时,对于工况1来说,裂缝出口处温度突降以后随着时间推移一直降低到10℃左右,而后温度随大气温度变化而变化;对于工况2和工况3来说,裂缝出口处温度突降以后随着时间推移先升高后降低到10℃左右,而后温度也随大气温度变化成正弦变化.工况2和工况3在裂缝出口处温度在突降以后先升高后降低是由于外界大气温度在坝体内传递具有滞后性造成的.

图5 不同高度心墙裂缝出口处温度与渗透流速变化对比

综上,当心墙出现裂缝时,心墙下游表面的温度和渗透流速都会发生突变.裂缝距坝底高度越高,发生集中渗漏时温差变化越大,渗透流速变化越小,裂缝距坝底高度越低,发生集中渗漏时渗透流速变化越大,温差变化越小.

5 结论

1）基于COMSOL Multiphysics软件的Ren（2019）模型适用于土石坝饱和-非饱和流-热耦合研究,可以准确地描述土石坝温度场和渗流场的动态变化过程.

2）集中渗漏通过改变渗流场来改变土石坝温度场的分布,可以通过温度场去反演渗流场去确定集中渗漏通道的具体位置.

3）心墙裂缝距坝底高度越高,发生集中渗漏时温差变化越大,渗透流速变化越小.在实际工程中可以在坝体内部铺设光纤来监测坝体的温度场变化,从而精确地确定发生集中渗漏的位置,为实际工程安全提供保障.