电渗透系数与土颗粒大小无关，因此利用电渗法对软土进行处理时能够得到不错的效果。土体的电渗现象最早由俄国科学家Reuss发现，但在相当长的时间内一直未得到工程应用，直到20世纪30年代，Casagrande^[1]首次将电渗排水法应用于铁路路基的加固并取得了良好的效果，这一新的方法立刻引起了工程界广泛的关注；Bjerrum等^[2]给出了电渗排水加固法中电学参数的设计方法，极大地指导了后续电渗法的应用研究。

基础性的试验研究^[3–4]是促进人们对电渗法认知的重要手段。电势差是电渗流发生的主要驱动力，不同的电压会导致不同强度的电势分布^[5–9]，使得土体中孔隙水的定向流动产生差异，产生不同的电渗加固效果。平行排布、梅花、矩形是最为常见的排布形式。王柳江等^[10]试验研究指出，电渗对象中阳极的数量越多，加固后土体强度越高。李一雯等^[11]的研究结果表明，平行排布产生的裂缝较为规律，平行错位形式下排水效果最好，而梅花形排布情况下土体裂缝开展情况较为严重。陶燕丽^[12]、张雷^[13]等试验研究指出阳极对于材料的敏感程度要高于阴极。郑凌逶等^[14]对电渗排水加固法的研究与应用进行了全面的总结，为后续的研究提供了方向性指导。众多电渗试验研究的成果，为后续电渗固结理论研究提供了必要支持。

Esrig^[15]最先提出了1维模型的电渗固结理论，水力渗流与电渗流满足线性叠加原理^[16]；后续学者在此基础上进一步拓展了电渗固结理论的研究内容 ^[17–19]。李瑛等^[20]建立了轴对称形式下的电渗固结理论；王军^[21]、王柳江^[22]等在1维电渗固结模型中引入了堆载作用的影响；吴辉^[23]、胡黎明^[24]等在电渗固结计算中考虑了电导率及土体参数的非线性，并进行了数值求解。也有学者开展了多物理场耦合作用下的电渗固结研究并进行了数值计算^[25–27]。苏金强等^[28]针对平行排布形式情况，基于“分块思路”建立了2维电渗固结理论。电渗固结理论研究为电渗法工程应用提供了必要的理论与技术支持。

然而，电渗试验研究^[5–12]表明电渗过程存在有效电势衰减现象，一方面，有效电势的衰减使得电渗后期效率低下，电渗法能耗较高在一定程度上限制了电渗法的应用；另一方面，电极间距相当时，电势场分布的2维效应显著，采用1维固结理论计算会高估电渗加固效果，而已有的电渗固结理论没有考虑有效电势衰减及电势场的2维特性带来的影响。针对上述问题，本文在前人研究的基础上，依据平行、矩形排布形式的组成规律，提出了点状电极单元的概念，并将其作为2维电渗固结理论的研究对象；针对研究对象中多电极的有效电势衰减情况，依据黎曼求和方式的思路，将有效电势衰减转化为离散、多态的阶梯曲线，并结合电渗与有效电势衰减在时间上的统一性，采取了初始条件“继承”的策略，将2维电渗固结全时间域上的求解转化为对单一黎曼状态段中间变量方程的反复迭代；基于Galerkin法给出了方程的有限元形式，并采用Python语言进行了数值实现；为了验证该算法及确定有效电势变化规律，进行了电渗试验研究，利用实测电势变化进行了多态下的计算，该方法能够体现出有效电势变化对电渗孔压的影响，且方法简洁、思路巧妙。

1.   点状电极单元

电渗排水固结法中，平行排布是最为常用的一种电极排布形式，当异性电极间距远大于同性电极间距时，通常会将其简化为1维模型，Esrig理论便是在这一基础上建立的。电渗机理指出土中孔隙水会沿着电势梯度下降的方向流动，即从阳极流向阴极，使得阳极区的土体加固效果较好。王柳江等^[10]对这一问题进行了试验研究，结果表明阳极数量越多则土体加固效果越好。基于这种考虑，将平行排布形式中相应的阴极替换为阳极，可以得到矩形排布形式，即1根阴极周围环绕8根阳极，如图1所示。

图  1  矩形电极排布形式

Fig.  1  Rectangular electrode arrangement

下载: 全尺寸图片

矩形排布的分布特性决定了该形式无法被简化为1维模型。无论是对平行形式还是对矩形形式的研究，依据排布形式的组成情况与重复规律，都可以归结为对最小电极单元的研究。每一个最小电极单元由4根电极组成，矩形排布形式下最小电极单元是非对称的，含有3根阳极与1根阴极，如图2所示，其中，b为单元宽度，h为单元高度。

图  2  非对称点状电极单元

Fig.  2  Asymmetric punctiform electrode unit

下载: 全尺寸图片

最小电极单元是一个完整的2维电渗单元，能够体现出电势场与渗流场的2维空间描述。最小电极组成单元中，电极自身的大小远远小于单元尺寸，所以可以将电极当作点状电极，这样单元的电势分布可以按照静电场的边值问题进行求解，电势场分布计算与电偶极子类似。

点状电极单元依据电极的极性可以分为对称和非对称两种形式，单元中边界分为同性边界与异性边界。电极极性的非对称特性使得单元电势场分布与均匀电场存在较大差异，此时不应简化为1维模型。

2.   有效电势衰减及多态处理

电渗过程中化学变化产生的气体及阳极区土体孔隙水的定向流动，使得阳极与土体之间接触电阻增大，造成施加在土体上的有效电势发生衰减，这一现象已在诸多电渗试验研究^[5–12]中得到了证实。

2.1   有效电势衰减下电渗固结计算的关键问题

事实上，有效电势衰减的现象早已被试验验证，但是电渗固结理论中一直缺乏对该现象的考量，根本原因是在电渗固结控制方程中引入中间变量时，如果电势分布与时间有关，那么中间变量中电势项对于时间的导数不为0，使得变换后的方程仍然具有两个场函数，无法进行求解。Esrig理论中电渗固结的控制方程如下^[15]：

式中： $u(x,t)$ 为孔隙水压力，Pa； $V(x)$ 为电势分布函数，V； ${k_{{\rm{e}}x}}$ 为 $x$ 方向的电渗透系数， ${{\text{m}}^{\text{2}}}{\text{/(V}} \cdot {\text{s)}}$ ； ${k_{{\rm{h}}x}}$ 为 $x$ 方向的水力渗透系数， ${\text{m/s}}$ ； ${\gamma _{\rm{w}}}$ 为水的重度， ${\text{N/}}{{\text{m}}^{\text{3}}}$ ； ${m_{\rm{v}}}$ 为土体体积压缩系数， ${\text{P}}{{\text{a}}^{{{ - 1}}}}$ 。

式（1）中，电势分布函数仅是关于空间坐标的函数，即表明电势分布与时间无关，为了进行方程的求解，引入中间变量 $\zeta (x,t)$ ，表示如下：

引入中间变量后，控制方程（1）等价于式（3）:

式中， ${C_{\rm{v}}}$ 为土体固结系数， ${C_{\rm{v}}} = {k_{{\rm{h}}x}}/{m_{\rm{v}}}{\gamma _{\rm{w}}}$ ， ${{\text{m}}^{\text{2}}}{\text{/s}}$ 。

变换后得到的中间变量方程式（3）是最简形式的抛物线型偏微分方程，此时该方程是可以进行解析求解的，在得到中间变量场的分布函数 $\zeta (x,t)$ 之后，利用式（2）可得到电渗的孔压分布。

上述是Esrig理论对电渗控制方程的求解思路与方法，可以看到控制方程能够进行求解的关键在于中间变量 $\zeta (x,t)$ 中电势分布项与时间无关，这样 $\zeta (x,t)$ 对时间 $ t $ 求导时，电势项可以当作常量，其对时间的导数为0，从而可以得到最简形式的抛物线型偏微分方程。然而，有效电势的衰减使得电势分布函数与时间相关，这一客观特性与中间变量函数引入时的要求相矛盾，这一矛盾关系制约着有效电势衰减下电渗固结方程的求解，是求解计算的难点也是关键点。此外，如果直接在中间变量表达中引入电势的时间效应，会导致数学逻辑的严密性不足。上述分析表明，有效电势衰减影响下的电渗固结计算的关键性问题在于如何处理有效电势衰减的数学描述与中间变量函数引入时的限制条件之间的矛盾关系。

2.2   有效电势衰减的离散多态方法

在2维模型中，有效电势是电势场边值问题的边界条件，有效电势 ${V_{\rm{e}}}(t)$ 随时间变化，引起电势场分布的变化，已有的研究成果中对于该问题的解答是完全迎合电势衰减的描述，而弱化中间变量引入过程中的数学要求。实际上，完全可以从一个相反的角度来考虑这个问题，即通过某种有效电势衰减的描述，未完全满足中间变量引入时的要求。

在微积分理论中，黎曼积分是一个十分重要、普遍的概念，它给出了定积分的精准定义，而左侧黎曼求和形式是定积分计算的基本方法和思路。对于任何一个有效电势变化函数 $ {V_{\rm{e}}}(t) $ ，该函数是时间的函数；当采用左侧黎曼求和形式时，如图3所示，在积分区段内进行子状态段的划分，每一个状态段内函数值均采用区间段左端点处的函数值，即图3中的 ${V_{{\rm{e}}1}}$ 、 ${V_{{\rm{e}}2}}$ 、 ${V_{{\rm{e}}3}}$ ，经过黎曼求和形式的划分，能够保证在子区间内函数值保持为常数，不会随自变量发生变化，使得在引入中间变量时对时间的导数为0，从而可以获得标准的抛物型偏微分方程。

图  3  左侧黎曼求和形式

Fig.  3  Left Riemann summation form

下载: 全尺寸图片

黎曼求和的处理方式使得任意一个状态段内的函数值是一个常数，同时在整个积分域上又能逼近原函数的变化。这一特性十分适合用于描述有效电势的衰减，因为在子状态段内有效电势是一个常数，那么中间变量引入时的条件自然可以得到满足，仅从这一点来看，在单一的子状态段内，该处理方式能够解决这一核心矛盾关系。

2.3   电势场分布

电渗法中施加的是直流电，因此电势场满足静电场理论，将 $ x - y $ 坐标系原点取在单元左下角的阴极处，单元模型如图2所示。对有效电势进行离散、多态处理后，记第 $ k $ 黎曼状态段内场状态变量电势分布函数为 ${V_k}(x,y)$ ，其满足Laplace方程，如式（4）所示：

有效电势衰减函数的离散、多态处理，是在时间域上的分割，而有效电势是电势场的边界条件，电势分布是一个边值问题，不同的黎曼段内的有效电势不同，导致各个状态段的电势场分布存在一定的差异。边界条件满足Dirichlet条件，见式（5）：

对式（4）采用Galerkin法进行离散，电势场的空间域记作 $\varOmega $ ，边界为 $\partial \varOmega $ ，单元类型采用三节点三角形单元，则单元内任一点处的场函数的近似值 ${\overline V_k}(x,y)$ 满足式（6）：

式中：N为形函数矩阵， ${\boldsymbol{N}} = \left[ {{N_1}}\quad {{N_2}}\quad {{N_3}} \right]$ ； ${{\boldsymbol{v}}_{{{\rm{k}}}}} = [ {{v_1}}\quad {{v_2}} {{v_3}} ]^{\text{T}}$ ，其中， ${v_1}、 v_2、 v_3$ 为单元节点处的电势值。

方程（4）对应的等效积分强的形式表示如下：

利用分部积分进行处理，离散时考虑边界条件为完全齐次，则分部积分的降阶项在边界 $\partial \varOmega $ 处为0，边界条件的影响在方程求解时利用“乘大数法”进行考虑。控制方程的等效积分弱形式表示如下：

将式（6）代入式（8），得到电势场控制方程的有限元形式表示如下：

通过式（9）可以获得每一个单元内的刚度矩阵；然后，通过集成得到总体的刚度矩阵，边界条件采用“乘大数法”引入。有效电势的衰减对于电势场的影响体现在边界条件的变化，而离散多态的处理方法能够为黎曼状态段内电势分布计算提供修正的边界条件。随着黎曼状态段数的增加，离散形式的描述会不断逼近原始的有效电势衰减函数，与此相应地需要付出更多次的稳态求解。

取单元大小为 $b = h = 0.4{\text{ m}}$ ，阳极处电势值为48 V，依据前述有限元形式（8）进行求解，可以得到非对称电极单元的电势场分布，如图4所示。

图  4  非对称电极单元的电势场分布

Fig.  4  Potential gradient direction field of asymmetric electrode element

下载: 全尺寸图片

图4中箭头为各点处的电势梯度的方向（场强方向），非对称情况下电势的分布与均匀电场差异很大，同时，与轴对称模型的电势分布也存在明显区别。

对称情况下，单元中存在2个阳极、2个阴极，电势场分布如图5所示，该情况下电势梯度在单元中部区域，梯度方向与均匀电场差异较大。无论是对称还是非对称形式，电极单元的电势场分布具有明显的2维特性，电渗中电势梯度决定了电渗流的大小及流动方向，因而采用2维电渗固结模型是较为合理的。

图  5  对称电极单元的电势梯度方向场

Fig.  5  Potential gradient direction field of symmetric electrode element

下载: 全尺寸图片

3.   2维电渗固结的多态继承方法

利用第2.2节黎曼求和的思路，对有效电势衰减函数进行离散多态处理后，可以得到多个黎曼状态段；在任意一个黎曼状态段内，有效电势 $ {V_{\rm{e}}} $ 是一个常量，每一个状态段都对应一个电势分布状态。此时，在状态段内电渗固结方程引入中间变量进行变换是完全成立的，但对于整个时间域上的求解，还需要考虑各个黎曼状态段内初始条件在时间维度上的先后顺序。

3.1   多态继承的基本思路

有效电势衰减是一个相对于时间而言的概念，因此离散多态的处理也只是在整个时间域上进行有限状态段的划分；电渗固结问题是一个初边值问题，从时间维度上来看，电渗孔压的非稳态与有效电势的衰减是同一个过程，在时间维度上具有统一性，对于有效电势离散多态的分割实际上也是对电渗固结非稳态的一个分割。

基于上述认识，可以认为多态处理后得到多个黎曼状态段，第 $k$ 状态段必然是发生在第 $k + 1$ 个状态段之前的。如果令每一个状态段的总时间长度均为 ${\text{d}}T$ ，那么在状态段 $k$ 内，场状态量在经历了 ${\text{d}}T$ 时间的变化之后，会得到该状态段末时刻的场函数值 ${\zeta _k}(x,y,{\text{d}}T)$ ；而对于下一状态段 $（k + 1）$ 来说，它的发生是在状态段 $k$ 的基础之上的，自然 ${\zeta _k}(x,y,{\text{d}}T)$ 就是 $k + 1$ 状态段内非稳态问题的初始条件 ${\zeta _{k + 1}}(x,y,0)$ 。若从第1个状态段开始进行各个状态段的遍历迭代计算，并记录每一个黎曼状态段的非稳态结果，最终可以获得整个电渗过程在全时间域上的解答。为了简化计算，可以取各个黎曼状态段的时间长度均为 ${\text{d}}T$ ，又因为中间变量场的边界条件全程保持不变，这样将对整个时间域上的求解转化为了对单一黎曼状态段上的反复迭代求解。

3.2   方法的有限元形式

由于该方法需要进行反复的迭代计算，因此，采用有限元法进行数值实现。点状电极单元中2维电渗固结控制方程的建立依然承认Esrig理论中的基本假设^[15]，此外加入以下假设：

1）不考虑单元外部的孔隙水补给；

2）单元中渗流场以电渗流为主，非对称单元中，异性边界上孔隙水的流动主要为电渗流，方向为电势梯度方向；

3）电极尺寸远小于单元几何尺寸，将其理想化为“质点”；

4）土体性质为各向同性介质，即 ${k_{{\rm{e}}x}} = {k_{{\rm{e}}y}} = {k_{\rm{e}}}$ ， ${k_{{\rm{h}}x}} = {k_{{\rm{h}}y}} = {k_{\rm{h}}}$ 。

水力渗流满足达西定律，电渗流与水力渗流满足线性叠加原理^[16]。经过多态处理后的，使得方程的求解域被限制在了子状态段内，子状态段内有效电势为常量，在状态段的时间范围内电势场是与时间无关的，本状态段内的电势分布记作 ${V_k}(x,y)$ 。依据连续性原理，得到第k黎曼状态段内单元的2维电渗固结控制方程式（10）：

式中，下标k表示物理量属于第k状态段。

由于在第 $k$ 状态段内有效电势值是一个常量，则该状态段内的电势分布函数 ${V_k}(x,y)$ 与时间无关，因此引入中间变量 ${\zeta _k}(x,y,t)$ 是完全成立的，见式（11）：

由此得到控制方程（10）等效变换后的形式，见式（12）：

网格划分及单元类型均与前述的电势场分布计算保持一致，单元的位移模式见式（13）：

式中： ${{\boldsymbol{\zeta }}_{\boldsymbol{k}}}(t) = {[ {{{\bar \zeta }_{k1}}}\quad {{{\bar \zeta }_{k2}}}\quad {{{\bar \zeta }_{k3}}} ]^{\text{T}}}$ ， ${\bar \zeta _{ki}}（i = 1,2,3）$ 为状态段k内t时刻单元节点处的场状态函数值； ${\bar \zeta _k}(x,y,t)$ 为状态段 $k$ 内，单元内一点处场函数近似值。

该问题为非稳态问题，采用Galerkin法进行空间项离散，将近似解 ${\bar \zeta _k}(x,y,t)$ 代入式（12）后，得到等效积分强形式，见式（14）：

采用分部积分进行处理，在齐次边界条件下，分部积分产生的降阶积分项在 $\partial \varOmega$ 处为0，得到中间变量方程对应的等效积分弱形式，见式（15）：

将单元位移模式式（13）代入式（15）后，得到式（16）：

式中， ${\boldsymbol{B}} = {C_{\rm{v}}}\displaystyle\iint {\left(\dfrac{{\partial {{\boldsymbol{N}}^{\rm{T}}}}}{{\partial x}}\cdot\dfrac{{\partial {\boldsymbol{N}}}}{{\partial x}} + \dfrac{{\partial {{\boldsymbol{N}}^{\rm{T}}}}}{{\partial y}}\cdot\dfrac{{\partial {\boldsymbol{N}}}}{{\partial y}}\right)}{\text{d}}x{\text{d}}y$ ， ${\boldsymbol{S}} = \displaystyle\iint ({{\boldsymbol{N}}^{\text{T}}} \cdot {\boldsymbol{N}}){\text{d}}x{\text{d}}y$ 。

这里采用的是三节点三角形单元，所以单元形函数为关于 $x、 y$ 的一次项，可以直接积分得出依据形函数的具体形式得到 ${\boldsymbol{B}}$ ；当积分点放在单元的重心处时，可以得到 ${\boldsymbol{S}}$ 的具体表述式如式（17）所示：

式中，A为有限单元的面积。

上述计算已经完成了空间坐标的半离散。为了对时间项进行离散。采用有限差分形式进行时间项的离散，考虑到迭代的稳定性，这里采用后向差分形式，对式（16）进行后向差分处理，黎曼状态段内时间步数为 $n$ ，因为每一个状态段的时间均为 ${\text{d}}T$ ，得到差分的时间步长为 $\Delta t = {\text{d}}T/n$ ，最终中间变量控制方程的有限元形式见式（18）：

式中， ${{\boldsymbol{\zeta }}_{\boldsymbol{k}}}(i - 1)$ 表示 $i - 1$ 时刻的场函数值向量， ${{\boldsymbol{\zeta }}_{\boldsymbol{k}}}(i)$ 为 $i$ 时刻的场函数值向量， $\Delta t$ 为差分的时间步长。

式（18）给出了单一状态段内中间变量控制方程的有限元形式，在边界条件及初始条件下，可以完成求解。第 $k$ 状态段内的初始条件需要继承第 $k - 1$ 黎曼状态段末时刻的值，即为式（19）：

式中： ${{\boldsymbol{\zeta }}_{\boldsymbol{k}}}(i = 0)$ 为第 $k$ 状态段内初始时刻时各节点的中间变量向量，从0开始计数； ${{\boldsymbol{\zeta }}_{{\boldsymbol{k - 1}}}}(i = n)$ 表示第 $k - 1$ 状态段末时刻的各节点处中间变量向量，从0开始计数，该状态段内时间步数为 $n$ ，末时刻 $i = n$ 。

依据第3.2节的假设1）～3），保持阴极开孔进行排水，而其他边界不排水，得到非对称单元中中间变量方程的边界条件，见式（20）：

通过上述有限元形式及（初始）边界条件，遍历每一个状态段可以获得全时间域上的有效电势衰减下的中间变量函数分布，然后通过式（13）求解出电渗孔压值。

4.   试验验证与算例分析

4.1   点状电极单元模型试验

为探究点状非对称单元的电渗固结特性，同时验证提出算法的可靠性，进行了非对称点状电极单元的电渗固结试验研究，模型尺寸为 $ 40{\text{ cm}} \times 40{\text{ cm}} $ ，电渗电源电压为48 V，阳极电极采用铁电极，阴极开孔为排水边界。设置了16个电势测点（包含4个电极），位置及编号如图6所示。图6中：a0、a3、a12、a15分别为电极界面电势探针，布置在距离电极表面约1 cm处，监测有效电势；孔压计采用动态微型负孔隙水压力计，型号为SCYG318，探头直径8 mm，量程–100～100 kPa，其距离上边界、右侧边界均为6.65 cm，利用RS485模块可读取孔压传感器信号，孔压采样间隔60 s。试验土体为重塑的粉质黏土，初始土样平均含水率约为29.4%，土样粒径分布如图7所示，试验装置如图8所示。

图  6  电渗模型及测点分布

Fig.  6  Electroosmotic model and distribution of measuring points

下载: 全尺寸图片

图  7  土样粒径分布

Fig.  7  Soil particle size distribution

下载: 全尺寸图片

图  8  非对称点状电极单元电渗试验装置

Fig.  8  Asymmetric punctiform electrode unit electroosmosis test

下载: 全尺寸图片

4.2   试验实测结果与分析

模型通电约51 h后，孔压计位置处土体裂缝的开展使得孔压读数为0，随即停止了试验。电渗过程中电流初期呈上升趋势，这是由土体的含盐量较高所导致的^[8]。电渗结束时电流下降了约45%，试验过程中总电源电压保持48 V不变，则总模型的电阻在增大，如图9所示。由图9可知：初始时刻模型总电阻约为35.87 $\Omega $ ，51 h电渗后模型总电阻增长到64.86 $\Omega $ ，并且模型的电阻在后期增大较快，导致电渗后期电渗效率较低；排水速度随着电渗的进行逐步降低，在48 h后，1 h内的排水量小于2.0 g。

图  9  电渗排水量、电流、电阻随时间的变化曲线

Fig.  9  Curves of electroosmotic discharge, current and resistance with time

下载: 全尺寸图片

由于电渗排水加固法具有加固不均匀的特性，试验后在模型阳极–阴极对角线上进行了含水率测试，结果见表1。电渗后，阳极附近土体含水率下降幅度较大，加固效果较好，测点15附近含水率下降9.85%；阴极由于是排水边界，在电势梯度作用下，模型中的孔隙水从阳极向阴极发生定向流动，使得阴极附近土体的含水率下降幅度较小，阴极处土体含水率下降约2.98%，排水加固效果不明显。

电渗过程中通过测量金属19个探针位置处的电势以确定点状电极单元的实际电势分布；其中，a0、a3、a12、a15处的电势值用于计算土体与电极之间因接触所损失的电压，从而得到真正施加在土体上用于电渗固结的有效电势。图10为界面电压随时间变化曲线。由图10可知：阴极处的界面电压变化幅度较小，总体呈下降趋势，中后期基本保持不变，这主要是由于阴极作为排水边界，土体中的水在阴极汇聚排出，使得阴极附近的土体含水率较高，土体电导率下降较小，从初始时刻的8.3 V下降到3.2 V。界面损失电压的降低，对于电渗效率具有积极作用。距离阴极较近的2个阳极a3与a12，变化趋势较为一致；界面电压随着电渗的进行，增长幅度较大，a3从3.7 V增加到21.7 V，增加幅度为486.5%；a12从6.6 V变化到23.5 V，增长幅度为256.1%。距离阴极较远阳极的界面电压a15总体呈下降趋势，但下降幅度较小，从初始时刻的13.2 V下降到6.4 V，下降幅度约为51.5%。

图  10  界面电压–时间变化曲线

Fig.  10  Curves of interface voltage with time change

下载: 全尺寸图片

将阴极与土体的接触界面确定为零势面，即可得到施加在土体上用于电渗作用的有效电势（电压），3个阳极处的有效电势变化曲线如图11所示。

图  11  有效电势–时间变化曲线

Fig.  11  Curves of effective potential with time change

下载: 全尺寸图片

由图11可知：电渗初始时刻，真正作用在土体上的电压分别为阳极a3为36.0 V，阳极a12为33.1 V，阳极a15为26.5 V；电渗结束后，阳极a3为23.1 V，阳极a12为21.3 V，阳极a15为38.4 V。有效电压的变化受到电极与土体接触界面的影响（图10）。土体在电场作用下失水，而受到渗透距离的影响，使得距离阴极（排水边界）较近的阳极区域土体含水率下降幅度较大，导致该区域的电阻增长较大，界面损失电压较高，有效电压随着电渗的进行发生衰减；阳极a15附近排水路径较远，土体电阻增加幅度较小，低于其他两个阳极处，而电源总电压保持不变，所以分担的有效电压是在增大的。点状电极单元具有显著的电势2维效应，并且电渗过程中土体受到物理、化学反应的影响，使得有效电势不断发生变化，从而导致电势分布时刻在变化，电渗的机理指出电渗固结过程中的负孔压幅值受到电势值的影响，如果忽略有效电势的变化，那么电渗固结计算后的负孔压值将被高估。总体上，有效电势近似呈线性变化，通过离散的数据电可以拟合出试验时间内有效电势的连续函数，进而采用多态方式进行处理。拟合得到的有效电势变化连续函数为：

图12为初始时刻实测电极单元电势分布。图13为电渗51 h后电极单元电势分布。电渗的进行伴随着土体电学特性的变化，有效电势发生变化，导致了点状电极单元的电势场分布发生了重分布现象。观察初始与电渗结束的电势分布，可以发现：阴极区附近土体电势随着电渗进行逐步降低，而阳极附近电势是在增加，两者的残差矩阵如图14所示。2维单元的电势变化主要沿着最远端阳极向阴极发生变化。

图  12  初始时刻实测电极单元电势分布

Fig.  12  Initial measured electrode element potential distribution

下载: 全尺寸图片

图  13  电渗51 h后单元电势分布

Fig.  13  Potential distribution after electroosmosis for 51 hours

下载: 全尺寸图片

图  14  电势残差矩阵分布

Fig.  14  Distribution of residual matrix of electric potential

下载: 全尺寸图片

4.3   算法验证

为了验证本文提出算法的可靠性及实用性，依据试验中提取到的有效电势变化曲线，依据实测的电势变化曲线进行电势场及电渗固结孔压的计算。土体电渗参数依据文献[29]的研究成果，土体电渗参数见表2。电渗中有效电势按照图11中曲线确定；多态分段数量m确定为1、10和20段，当状态段数为m=1时，就是不考虑电势变化下的计算结果，为了保证具有对比性，3种计算方案中均保证相同的时间步长，总时间取51 h，总时间步数为n=10 000，计算方案详细参数见表3。单元尺寸为 $b = h = 0.4{\text{ m}}$ ，网格单元采用三节点的Delaunay三角形单元，模型初始孔压值依据孔压计实测数据进行确定，初始孔压值为1.465 kPa。

采用Python语言自行编程开发了PyEcFem库，并进行电势多态的分段，以及多态下电势稳态分布和2维电渗固结的多态继承的计算，本系统具有完整的后处理功能。状态段数为10时的电极有效电势的多态分割曲线如图15所示。

图  15  有效电势多态分割

Fig.  15  Effective potential polymorphism segmentation

下载: 全尺寸图片

利用实测的有效电势数据，通过拟合构造得到近似的连续函数，如图11所示，具体表达式见式（21）。非对称单元的3电极处的有效电势变化接近一次函数，在获得拟合函数后，利用多态方法进行处理，具体如图15所示。多态方法的优势在于，对函数形式没有具体要求，在获得拟合函数的表达式基础上，均可以对其进行多态的分割处理。

在电极有效电势进行分割的基础上，进行每一个状态段的电势分布计算，由于篇幅有限，这里仅给出状态段数为10初始状态与末时刻（51 h）的电势残差分布，如图16所示。与电势残差矩阵相比，两者变化具有相似的规律：较远处的阳极附近电势值在增加，其余两个阳极处电势值下降，电势分布下降区域主要分布在单元中部。这一结果反映出多态继承算法下的电势场能够较好地符合实际情况。

图  16  有效电势残差矩阵分布（状态段数为10）

Fig.  16  Distribution of residual matrix of electric potential (status segment number is 10)

下载: 全尺寸图片

采用微型孔压传感器进行了电渗孔压的监测，受室内模型尺寸的影响，随着排水的进行，土体裂缝在持续的开展，在49.6 h时孔压计附近处土体裂缝贯通，使得孔压计读数回弹至0；此外，需要注意的是，电渗过程中的电化学反应会对传感器进行腐蚀，阴极附近对于金属的腐蚀强烈，腐蚀会使阴极处的孔压计读数在约10 h左右出现异常。

孔隙水压力（简称为孔压）实测值如图17所示，在49.6 h达到–21.53 kPa，阳极15处有效电势随着电渗的进行增加，如图11所示。因此，对于理论计算来说，如果不考虑这一种变化，那么得到的孔压幅值是偏小的，考虑电势变化得到的孔压更符合实际情况。电渗后多态继承算法在分段数为10时的孔压值达到–29.01 kPa，而不考虑衰减的结果为–14.11 kPa。

图  17  实测模型孔隙水压力曲线

Fig.  17  Measured model pore water pressure curves

下载: 全尺寸图片

4.4   多态方法的特性

由于第4.3节的计算是依据试验而进行的，时间较短，土体电渗固结过程不够充分，为了探究多态继承算法数值结果的特性，设置电渗时间为1×10⁸ s，有效电势变化仍然按照图11来确定，其他参数保持不变，进行状态数段数为1、10、20下的计算并对比。图18为状态段数为20时的电渗过程中第1 000步、5 000步和8 000步的负孔压分布云图。

图  18  有效电势衰减下孔压分布云图（状态段数为20）

Fig.  18  Pore pressure distribution under effective potential attenuation (status segment number is 20)

下载: 全尺寸图片

对于点状电极单元而言，由于电极的非对称性，使得电势分布及电渗孔压分布与1维理论存在较大的区别；距离排水边界越近，渗流路径越短，土体孔隙水排出越容易，从而失水导致土体电导率发生变化，进而导致界面接触电压改变，较远处阳极受到渗流路径的影响，有效电势在增大，使得该位置处的土体负孔压值最大，加固效果最好，其他两个阳极处次之。图19为不考虑电势衰减下的电渗固结孔压分布。

图  19  不考虑有效电势衰减下孔压分布云图

Fig.  19  Pore pressure distribution without considering effective potential attenuation

下载: 全尺寸图片

对比图19与图18的结果可知，不考虑有效电势的影响得到的负孔压幅值要高于衰减条件下的结果，说明出多态继承算法能够体现出电渗过程中有效电势变化对电渗固结的影响。

图20为不同位置处孔隙水压力变化曲线。由图20可知：距离阳极较近的点，电渗产生的负孔压值较大，符合经典电渗固结理论中电势值决定负孔压大小的原理；位于单元中部的点距离阳极较远，电势值较低，因此负孔压值要小于阳极处的。考虑有效电势变化下的计算结果负孔压值要低于有效电势不变的结果，即经典理论中如果不考虑电势变化，会高估电渗产生的负孔压值。随着分段数目的增多，由于继承所带来的突变现象逐渐减弱，曲线趋于平滑，高电势区域继承带来的突变越明显。因此，本文提出的算法能够反映出有效电势对于电渗的影响，并且计算逼近具有稳定性。

图  20  不同位置处孔压变化曲线

Fig.  20  Pore water pressure curves at different locations

下载: 全尺寸图片

电渗固结与水力排水固结不同，电渗排水过程会产生负孔压，因此，孔压的消散往往指的是从正的超孔压发展为负孔压，描述这个过程的固结度定义与普通水力排水固结中的固结度定义不同。固结度本质反映的是实时孔压的变化与孔压变化最大幅度之间的比值，由于有效电势的衰减变化，负孔压并非在最后时刻达到极值，使得电渗中土体固结度的定义需要被修正，见式（22）：

式中， $\overline U $ 为平均固结度， ${u_0}$ 初始时刻的孔压值， $u{(x,y,t)_{\max }}$ 为电渗孔压幅值的最大值。

式（22）分母表示总的孔压变化范围，即初始孔压到负孔压幅值最大值之间的距离；分子项反映实时孔压变化的程度。通过这样的定义可以反映出电势变化对固结度的影响。点状电极单元整体平均固结度及平均孔压曲线如图21所示。由图21可知：在电渗后期固结度曲线会出现回弹现象，说明在之前的电势分布下已经达到负孔压极值，由于有效电势的衰减使得后续难以维持这一种较高的幅值状态，会逐步降低回弹，土中的单元平均孔压曲线反映出这一现象，最终固结度必定小于1。m=10时最后时刻的固结度为0.980，m=20时最终固结度衰减到0.979；而不考虑衰减时，最终时刻固结度为1，曲线单调增长。

图  21  单元平均固结度及平均孔压曲线

Fig.  21  Unit average consolidation degree and average pore water pressure curve

下载: 全尺寸图片

图22为单元平均孔压与经典结果（分段数为1时的计算结果）的偏差曲线。m=10与m=20这两种情况下差别较小，与经典结果相比，主要在电渗后期残差较大，分别达到9.160 kPa与9.516 kPa，意味不考虑有效电势变化的结果比多态结果孔压高出9.160 kPa与9.516 kPa，如图22（a）所示。多态算法可以有效地体现出有效电势对电渗中孔压的影响，针对有效电势变化下的电渗固结求解是具有合理性的。状态段数的增加导致逼近误差的降低，为了体现状态段数目对计算的影响，进行了m=10与m=20下孔压偏差研究，如图22（b）所示。偏差表现出初始条件继承带来的突变现象，随着时间呈锯齿态增加，从10段增大到20段，计算结果的偏差最大值为0.353 kPa；曲线锋与谷之间的差值随着时间不断降低，从初期的0.12 kPa降低到0.07 kPa。因此，状态段数没有必要设定太多，在曲线已经较为平滑的情况下，再增加状态数目对提升计算精度贡献并不高。

图  22  单元平均孔压与经典结果的偏差曲线

Fig.  22  Deviation curve between unit average pore pressure and classical results

下载: 全尺寸图片

5.   结　论

针对同性电极间距较大排布形式下的电渗固结计算问题，将实际模型抽象为点状电极单元，有利于该问题的理论研究，并针对有效电势衰减条件下的电渗固结计算，提出了多态继承的算法，建立了点状电极单元的2维电渗固结的控制方程，利用Galerkin法进行了方程的离散，给出了算法完整的有限元格式，基于Python语言进行了数值实现；同时，进行了非对称形式下点状电极单元的电渗试验，并与所提出算法的数值结果进行了对比分析。主要结果如下：

1）异性电极间距与同性电极间距相差不大时，单元的电势分布具有明显的2维效应，与1维电渗固结理论有较大区别；点状电极单元能够描述不同极性匹配下的电渗模型，有利于电渗固结理论的研究。

2）利用黎曼求和思路进行有效电势变化曲线的处理，使得在某一状态段内满足中间变量引入的条件限制，并利用初始条件继承的方式，将全时间域内的求解转化为子状态段内的迭代计算。

3）非对称点状电极单元的电渗试验表明，单元整体上有效电势是在降低的，而渗流路径的长短会对有效电势的衰减变化产生影响。距离阴极（排水边界）较远的阳极处，有效电势随时间增长；而其余两个阳极距离阴极距离相同且较近，有效电势处于衰减状态。电渗过程中单元内电势分布随有效电势的变化而变化，主要衰减区域处于单元两个阳极对角线附近，阴极电势衰减幅度较小。

4）多态继承算法下的电渗固结计算结果能够反映出有效电势变化对电渗中负孔压分布的影响。初始条件继承时造成的孔压突变幅度受到状态段数的影响，越接近阳极，幅度越大。当孔压曲线较为平滑后，继续增大状态段数，数值解精度提升有限，不利于计算效率的提升。