基于热−水−力耦合的电渗排水试验数值模拟

王柳江，刘斯宏，陈守开,，周斌

(1. 河海大学 水利水电学院，江苏 南京，210098；2. 水利部堤防民病害防治工程技术研究中心，河南 郑州，450003；3. 华北水利水电大学 水利学院，河南 郑州，450011)

基于热−水−力耦合的电渗排水试验数值模拟

王柳江1,2，刘斯宏1，陈守开2,3，周斌1

(1. 河海大学 水利水电学院，江苏 南京，210098；2. 水利部堤防民病害防治工程技术研究中心，河南 郑州，450003；3. 华北水利水电大学 水利学院，河南 郑州，450011)

基于多孔介质理论和热力学理论，考虑电渗产生的温度对土体变形及排水性质的影响，建立热−水−力耦合作用下的电渗排水多场耦合数学模型，并利用Fortran语言编制计算程序。在此基础上开展电渗排水现场试验的数值模拟，并与现场监测结果进行对比。研究结果表明：建立的模型可较好地模拟电渗过程中的温度场、渗流场及应力场变化和分布规律。土体温度随能耗增大而升高，尤其对于电极附近土体；阳极周围土体孔隙水压力减小，而阴极周围土体孔隙水压力增大，且考虑温度作用下的电渗流量随电渗的持续而增加；地基沉降由土体温度上升引起。

电渗；热−水−力耦合；排水；数值模拟；软土

电渗是指在土体中通直流电以带动阳极水分往阴极迁移的现象。由于电渗流速不受土体黏粒粒径的影响，因此，电渗法可显著提高水分在低渗透性黏土中的排水速度[1]。近年来，随着河道清淤及土壤修复工程的兴起，大量压缩性强、渗透性低及污染物含量高的软弱黏土亟需处理[2−3]。由于电渗法可有效克服传统排水固结法针对低渗透性黏土排水速率缓慢、处理周期长且效果不明显的缺点，在环境岩土领域得到了极大关注，并开始在一些高含水量软土及其地基的处理和污染土修复的环保工程中得到了应用[4−6]。电渗除带动水分迁移之外，通常还伴随着热能的产生和释放，从而导致土体温度上升[7]。相关研究表明，温度上升不仅对土壤的水力和力学性质具有一定影响，而且针对含有可挥发性有机化合物(VOCs)的污染土，还可有效去除土壤中污染物[8−9]；对电渗而言还会引起土体体积收缩使得电极与土体之间产生裂缝，进而增大两者之间的接触电阻并影响电渗效率；电渗过程中温度场的变化可对土体的渗透、力学特性产生一定的影响，同时渗流场与应力场的变化又反过来影响温度场的变化分布规律，三者之间是相互耦合的[10−11]。因此，考虑电渗过程中的土体温度变化对进一步了解土体渗流、变形及有机化学物分解变化规律，确定施工参数，优化电渗工法具有一定的理论指导意义。本文作者在文献[12−13]的基础上，建立热−水−力耦合的电渗数学模型，开发相应的有限元程序，并对电渗排水的现场试验进行数值模拟，通过与现场监测资料的对比验证了该电渗数学模型及程序的合理性，并进一步对电场、温度场、渗流场及应力场的变化分布规律进行分析。

1　数学模型

1.1基本假设

1) 忽略土体内化学、相变等因素的影响；

2) 忽略土体内细颗粒发生电泳而产生的电流；3) 设土体始终处于饱和状态；

4) 土体内电荷传递满足欧姆定律，渗流满足达西定律；

5) 热量传递满足傅里叶定律，并忽略对流影响；

6) 由电势梯度和水头梯度产生的流量可以叠加；

7) 由土体中离子移动而产生的电场相对外加电场忽略不计。

1.2静电场方程

根据电荷守恒原理，假设电流处于稳定状态，可建立计算电势分布的控制方程：

式中：V为电势；Re为电流汇源项；eσ为土体温度为T时的电导率(Ω−1·m−1)，与土体的含水率和温度有关，采用经验公式表示如下[14]：

式中：T为土体温度；25σ为25 ℃时的土体电导率；α为试验常数，反映土体的性质，随含水率变化而变化。

1.3热传导方程

根据能量守恒原理，可建立热传导方程表达式：

式中：cρ和wwcρ分别为土体和水的热容量；hk为土体渗透系数；λ为热传导系数；Q为热量汇源项，在该问题中主要为土体电阻产热量，其单位体积产热可表示为

1.4渗流连续性方程

根据质量守恒原理，假设电势和水头梯度产生的流量可以叠加，则电渗流连续性方程表示为

式中：

H为水头；u为位移；wγ为水的容重；n为孔隙率；sμ为贮水率或释水率；Ks与Kw分别为土颗粒与水的体积压缩模量；keo为电渗透系数；β为热膨胀系数；Rh为流量汇源项。

1.5静力平衡方程

土体的运动方程以总应力形式表示如下：

式中：σ为总应力；ρ为土体的密度。

假设土体始终处于饱和状态，土体中的孔隙气压等于大气压，则有效应力可表示为

式中：σ′为有效应力；m={1,1,1,0,0,0}；uw为孔隙水压力。

考虑温度及土体的非线性变形特性，土体的应力−应变关系可表示为

其中：DT为土体的切线刚度矩阵；dεT为热应变增量，dεT=m(β/3)dT 。

2　土体本构模型及有限元实现

2.1考虑温度影响的土体本构模型

由于电渗引起土体温度升高，而温度的变化又影响了土体的力学特性。因此，引用文献[15]中在修正剑桥模型基础上提出的考虑温度影响的正常固结饱和土热弹塑性模型的屈服函数：

式中：p′和q分别为平均有效应力和广义剪应力；M为临界状态应力比，假设M不随温度变化；为当前温度下饱和土的前期固结应力，

2.2应力−应变关系式推导

对于非等温弹塑性变形过程，土体应变增量可表示为弹性和塑性应变增量之和：

式中：De为弹性刚度矩阵；dλ为比例常数；Q为塑性势函数。

由非等温条件下的土体屈服函数可知应力、前期固结应力以及温度为状态变量，根据一致性条件，屈服函数的增量形式可表达如下：

式中：e0为初始孔隙比；λ和κ为试验参数，通过土体的等向压缩/回弹试验获得；根据塑性流动法则，塑性体应变增量可表示为

将式(14)代入式(13)后，再将式(11)和式(13)代入式(12)，即可得到比例常数dλ：

将比例常数式(15)代入式(11)后得到考虑温度影响的土体弹塑性应力应变关系式：

其中：

2.3有限元离散

将位移、水头、温度和电势作为基本未知量，在控制方程(1)，(3)，(5)和(6)中引入相应的边界条件，并采用Galerkin加权余量法对控制方程进行空间离散，得到以下有限元数值矩阵列式：

其中：

3　工程实例

3.1现场试验

电渗排水现场试验[16]布置如图1所示，试验区域主要位于两电极之间，阳极和阴极分别埋设在地表以下0.4 m和2.2 m处。阳极为长×宽为3.4 m×3.1 m的钛合金网，其表面镀有金属氧化物以防腐蚀，埋设时先开挖表层土体至地表以下0.4 m处，然后将阳极金属网放入并在其上方铺设厚度为2～3 cm的透水砂层，最后回填开挖土至地表；阴极采用石墨颗粒，通过水力劈裂的方法将其注入地表以下2.2 m，形成半径为2.5 m，厚度为0.2～0.6 cm的石墨导电层。试验时间为1996−08−02—1996−12−31，整个试验过程中两电极中的水头通过水泵控制保持不变，其中阳极水头固定在地表以上7 cm处，而阴极水头保持在地表以下73 cm处。通电开始时间为1996−08−06，在此之前主要为自重排水阶段，电渗初期的电源输出功率控制在2 kW左右(电势47 V，电流42 A)；当电能消耗达到1 500 kW·h后，输出功率增加至3.9 kW(电势73 V，电流53 A)。整个试验共耗能5 500 kW·h，试验过程中对排水量、孔隙水压力、电势、电流以及温度随能耗的变化规律进行了监测。

图1　电渗现场试验示意图Fig. 1 Schematic illustration of electro-osmosis field test

3.2计算模型及条件

为模拟电渗排水现场试验，本次计算所取范围如图2所示，其长×宽×高为25 m×20 m×20 m。设地表高程为0，且模型沿电极中心对称，阳极和阴极间距为1.8 m，阳极位于地下0.4 m处，其表层铺设有3 cm的砂层，由于阳极厚度仅0.9 mm，为简化模型，将阳极与砂层视为同种材料，其相应的材料特性可综合钛合金网与砂层的性质确定。阴极位于地表以下2.2 m处，为半径2.5 m、厚度0.5 m的薄层。

图2　计算范围Fig. 2 Diagram of simulated domain

电渗法是一个涉及电场、渗流场、温度场及应力场等多场耦合的复杂问题，计算时需设定相应的初始条件和边界条件。由于所取计算范围水平向距离为阴极直径的4倍，竖向深度约为两极间距1.8 m的14倍，因此可将模型四周及底部作水头和温度的固定边界，其中约束水头为地下水位(本次试验为−0.73 m)，边界温度随地基深度减小，其温度梯度为2～3 ℃/m。由于地表铺有土工膜防渗，对渗流和电场计算而言地表可视为不透水绝缘边界。此外，如图1所示电源通过导线与电极连接，其中阳极连接位置在钛合金网的边缘上，而阴极在石墨薄层的中心位置，因此接通电源时这2处的节点电势分别为47 V和0 V。同时，两电极处的水头通过水泵控制，则井点管与阳极和阴极的连接部位可设为水头边界，相应的节点水头分别为0.07 m和−0.73 m。对于位移边界，模型底部为三向位移固定边界，模型四周为法向位移约束边界，地表为自由边界。

根据土体性质沿深度分布情况，可将地层划分为如下3种材料：(Ⅰ) 0～1.8 m；(Ⅱ) 1.8～2.0 m；(Ⅲ) 2.0 m以下。本次计算所采用的模型参数主要包括该3种地层材料、阳极与砂层复合体以及阴极的电导率、电渗透系数、渗透系数、热容、导热系数以及热弹塑性本构模型相关的参数。表1所示为计算电场、渗流场及温度场所需的模型参数，取自文献[12]。由于计算中涉及地下水位的变动，表1中还给出了给水度μy，考虑到地基土体均为黏性土，其给水度设为0.05；而复合阳极与阴极材料性质类似砂性土，设为0.4。由于缺少各土层材料的力学特性资料，本次计算选取文献[15]中的热弹塑性本构模型参数：M=1.05，λ=0.14，κ=0.05，e0=0.8，γ=0.15，ν=0.3，β=0.9×10−6℃−1。由于电极相对较小，对地基变形影响可忽略，而各土层材料为性质相近的黏性土，可认为各材料的力学本构模型参数一致。

3.3计算结果分析

计算主要模拟了电渗排水现场试验的前半阶段：首先通过井点管控制电极中的水头，使水分在自重作用下由阳极流向阴极，使得两极间土体达到饱和状态；然后开始通电，控制电压为47 V，电源输出功率约为2 kW，持续通电750 h。

电势场等势线的分布情况如图3所示。其中r=0.6，1.2，2.1及3.0 m为试验时距阴极与电源导线接触点的距离，试验中在这些位置设置了竖向电压和温度测管。从图3可见：等势线在阳极与阴极之间较为密集且基本平行，从电极边缘往外开始弯曲并呈波状往外逐渐扩散，说明电极间的电渗作用最显著。此外还可以看出：阳极与阴极中的电势分布较为均匀，说明电极内部电势损失较小，而电极中的电势损失主要在电极与导线的接触部位。

图3　电场等势线分布Fig. 3 Distribution of electric potential

图4所示为不同位置的电势随深度分布的计算值与实测值的对比情况。由于电极与周围土体之间存在界面电阻，可以看到施加在土体中的有效电压有所减小，通过计算得到的两极间有效电压为39 V，较实测有效电压35.7 V略大，主要电势损失集中在阴极。此外，两极间电势分布由2段不同斜率的直线组成，其中地表以下0～1.8 m处的电势梯度较小，这是由于该土层材料的导电性比Ⅱ层和Ⅲ层的大。同时，计算得到的1.8 m以下土层电势梯度较实测值大，说明该土层材料的实际电导率大于室内试验值0.02 S/m，这主要是由于地基内土体受地应力的影响，土体所受的地应力越大，表明土体越密实，其对应的电导率也越大。r=0.3 m时所对应的电势梯度显著减小，这主要与该位置已超出了阳极尺寸范围有关。

表1　计算模型参数Table 1 Calculation model parameters

图4　电势随深度分布的计算值与实测值对比Fig. 4 Comparison of simulated and observed voltage distribution with depth

图5　温度随深度分布的计算值与实测值对比Fig. 5 Comparison of simulated and observed temperature distribution with depth

图5所示为土体温度沿地基深度的分布情况。从图5可以看到：温度的计算值与实测值基本吻合。地基内的温度变化分布规律如下：土体温度随着能耗的增加而升高；土体温度沿电极径向逐渐减小；阴极与电源导线连接处的接触电阻较大，相应的电势损失也较大，导致该位置附近的土体温度上升最为显著；温度沿地基深度基本呈减小趋势，而r=0.6 m时的温度沿深度方向却基本相同，说明受阴极与电源连接处持续产热的影响，其周围土体温度迅速升高，其相应的竖向温度梯度也迅速减小。通过上述现象可知，在进行电渗排水固结分析时，若电渗持续时间较长、能耗较大，则考虑土体的温度变化是十分必要的，尤其对于电极周围的土体，由于温度上升剧烈，其周围土体的力学、渗透特性的改变将十分明显。

图6所示为通电24 h后计算得到的孔隙水压力增量等值线分布规律。由图6可以看到：孔隙水压力的增量的分布规律与电势分布基本一致，这说明电渗排水过程中孔隙水压力的变化主要由电势梯度作用下的阳离子带动阳极水分往阴极流动引起的。在此作用下，阳极附近的孔隙水压力逐渐减小，而阴极的则增大。

图6　电渗后水头变化等值线分布Fig. 6 Distribution of hydraulic head increment after electro-osmosis

图7所示为电极周围土体孔隙水压力的变化情况。由图7可见：水头变化计算值与实测值十分接近，通电24 h后，阳极周围土体孔隙水头的减小了0～10 cm，阴极周围则增大了2～8 cm。为控制阴极水头保持地下0.73 m不变，需将通过电渗迁移至阴极的水分排出，实测的阴极电渗平均排水速率为0.6～0.8 L/h。图8所示为考虑和不考虑温度影响得到的阴极电渗排水速率，其平均值分别为0.8 L/h和0.68 L/h，均与实测值比较接近，说明本文采用的计算模型能够较好地预测电渗流速。由图7还可以看出：考虑温度影响时的电渗排水速率随能耗增大逐渐增大，而不考虑温度影响时电渗排水速率则基本不变，其主要原因是温度上升导致了土体塑性体积收缩，从而促使了土体孔隙水的排出，使得电渗排水速度逐渐增大。

图7　两极附近水头变化计算值与实测值对比Fig. 7 Comparison of simulated and observed hydraulic head increments around electrodes

图8　电渗排水速率随能耗的变化Fig. 8 Variation of the electro-osmotic flow rate with energy input

图9　电渗结束后沉降、温度分布等值线Fig, 9 Distribution of settlement and temperature after electro-osmosis

图9所示为试验结束时地基沉降和温度的分布规律。从图9可以看到：试验结束后地基主要表现为沉降变形，最大沉降发生在阳极中心位置，沿电极径向依次减小，这与温度在地基内的分布规律极其相似。考虑到本试验中孔隙水压力的变化较小，为0.2～ 1.0 kPa，对地基变形的影响可忽略，因此可以认为温度上升是导致地基产生沉降变形的主要因素。

4　结论

1) 电极间电势分布基本呈一维线性分布，电势损失主要发生在阴极与电源导线连接位置。

2) 由于阴极与导线接触部位产生大量热能，阴极周围土体温度上升最快，且温度沿电极径向往外扩散降低。

3) 土体受温度上升影响产生塑性体积收缩，地基产生沉降变形，其分布规律与温度分布规律基本相似。

4) 通电后孔隙水由阳极往阴极流动，电渗流量在不考虑温度影响时迅速增大并保持不变，而考虑土体温度上升的影响时，电渗流量呈逐渐增大趋势，这主要归因于土体温度上升引起的体积收缩变形。

[1] PROBSTEIN R F, HICKS R E. Removal of contaminants from soils by electric fields[J]. Science, 1993, 260: 498−503.

[2] 文海家, 严春风, 汪东云. 吹填软土的工程特性研究[J]. 重庆建筑大学学报, 1999, 21(2): 79−83. WEN Haijia, YAN Chunfeng, WANG Dongyun. Some engineering properties of the dredger fill[J]. Journal of Chongqing Jianzhu University, 1999, 21(2): 79−83.

[3] 朱伟, 张春雷, 刘汉龙, 等. 疏浚泥处理再生资源技术的现状[J]. 环境科学与技术, 2002, 25(4): 39−41. ZHU Wei, ZHANG Chunlei, LIU Hanlong, et al. The status quo of dredged spoils utilization[J]. Environmental Science and Technology, 2002, 25(4): 39−41.

[4] 房营光, 徐敏, 朱忠伟. 碱渣土的真空−电渗联合排水固结特性试验研究[J]. 华南理工大学学报(自然科学版), 2006, 34(11): 70−75. FANG Yingguang, XU Min, ZHU Zhongwei. Experimental research on the consolidation characteristic of the vacuumelectro osmosis combination drainage of soda residue soil[J]. Journal of South China University of Technology (Natural Science Edition), 2006, 34(11): 70−75.

[5] 龚晓南, 焦丹. 间歇通电下软黏土电渗固结性状试验分析[J].中南大学学报(自然科版), 2011, 42(6): 1725−1730. GONG Xiaonan, JIAO Dan. Experimental study on electroosmotic consolidation of soft clay under intermittent current condition[J]. Journal of Central South University (Science and Technology), 2011, 42(6): 1725−1730.

[6] LO K Y, MICIC S, SHANG J Q, et al. Electrokinetic strengthening of a soft marine sediment[J]. Int J Offshore Polar Eng, 2000, 10(2): 137−144.

[7] BURNOTTE F, LEFEBVRE G, GRONDIN G. A case record of electroosmotic consolidation of soft clay with improved soil electrode contact[J]. Can Geotech J, 2004, 41(6): 1038−1053.

[8] BUETTNER H M, DAILY W D. Cleaning contaminated soil using electrical heating and air stripping[J]. J Environ Eng, 1995, 121(8): 580−589.

[9] HERON G. Soil heating for enhanced remediation of chlo-rinated solvents: a laboratory study on resistive heating and vapor extraction in a silty, low-permeable soil contaminated with trichloro-ethylene[J]. Environmental Science Technology, 1998, 32: 1474−1481.

[10] TANG A M, CUI Y J. Modelling the thermomechanical volume change behaviour of compacted expansive clays[J]. Géotechnique, 2009, 59(3): 185−195.

[11] 蔡国庆, 赵成刚. 非饱和土渗流和变形强度特性的温度效应[J]. 力学进展, 2010, 40(2): 147−156. CAI Guoqing, ZHAO Chenggang. Temperature effect on seepage and strength deformation characteristics of unsaturated soils[J]. Advances in mechanics, 2010, 40(2): 147−156.

[12] CHEN Jianlong. Effects of electroosmosis on soil temperature and hydraulic head. Ⅱ: numerical simulation[J]. Journal of Environmental Engineering, 2002, 128(7): 596−603.

[13] 王柳江, 刘斯宏, 汪俊波, 等. 电场−渗流场−应力场耦合的电渗固结数值分析[J]. 岩土力学, 2012, 33(6): 1904−1911. WANG Liujiang, LIU Sihong, WANG Junbo, et al. Numerical analysis of electroosmostic consolidation based on the coupling of electrical field, seepage field and stress field[J]. Rock and Soil Mechanics, 2012, 33(6): 1904−1911.

[14] 李瑛, 龚晓南, 郭彪, 等. 电渗软黏土电导率特性及其导电机制研究[J]. 岩石力学与工程学报, 2010, 29(S2): 4027−4032. LI Ying, GONG Xiaonan, GUO Biao, et al. Research on conductivity characteristics of soft clay during electro-osmosis and its conductive mechanism[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(S2): 4027−4032.

[15] 姚仰平, 牛雷, 杨一帆, 等. 考虑温度影响的非饱和土本构模型[J]. 岩土力学, 2011, 32(10): 2881−2888. YAO Yangping, NIU Lei, YANG Yifan, et al. Constitutive model for unsaturated clays considering temperature effects[J]. Rock and Soil Mechanics, 2011, 32(10): 2881−2888.

[16] CHEN Jianlong. Effects of electroosmosis on soil temperature and hydraulic head. I: field observations[J]. Journal of Environmental Engineering, 2002, 128(7): 588−595.

(编辑 赵俊)

Numerical simulation of electro-osmostic drainage test based thermo-hydro-mechanical coupling model

WANG Liujiang1,2, LIU Sihong1, CHEN Shuokai2,3, ZHOU Bin1
(1. College of Hydropower and Water Conservancy, Hohai University, Nanjing 210098, China; 2. Dyke Safety and Ministry of Water Resources Engineering Research Center for Disease Control, Zhengzhou 450003, China; 3. School of Water Resources, North China University of Water Resources and Electric Power, Zhengzhou 450011, China)

The thermo-hydro-mechanical coupling mathematic model for electro-osmosis consolidation was established based on porous medium and thermodynamics theory, considering the effects of temperature induced by electro-osmosis on soil deformation and drainage properties. Then the computational program was compiled with Fortran language. The numerical simulation program was verified by comparison with field test of electro-osmosis drainage. The result shows that the proposed model for electro-osmosis consolidation can accurately simulate the variation and distribution of temperature, seepage, and stress field in the process of electro-osmosis. The temperature of the soil increases with the increase of energy input, especially in the vicinity of electrodes. The hydraulic head decreases near the anode, while it increases around the cathode. The simulated electro-osmotic flow rate increases with the duration of electro-osmosis considering the action of temperature. The settlement of surrounding soils is also induced due to the increase of temperature.

electro-osmosis; thermo-hydro-mechanical coupling; drainage; numerical simulation; soft clay

TU443

A

1672−7207(2016)03−0889−08

10.11817/j.issn.1672-7207.2016.03.023

2015−03−14；

2015−05−15

中国博士后科学基金资助项目(2014M561566)；江苏省博士后科研资助计划(1301015A)；水利部土石坝破坏机理与防控技术重点实验室开放基金资助项目(YK913007)；水利部堤防安全与病害防治工程技术研究中心开放课题基金资助项目(2014002) (Project(2014M561566) supported by China Postdoctoral Science Foundation; Project(1301015A) supported by Postdoctoral Research Funds of Jiangsu Province; Project(YK913007) supported by Key Laboratory of Earth-Rock Dam Failure Mechanism and Safety Control Techniques; Project(2014002) supported by Dyke Safety and Ministry of Water Resources Engineering Research Center for Disease Control)

王柳江，博士(后)，从事非饱和土本构理论及地基加固技术研究；E-mail: 15850514642@163.com