饱和土的变渗透系数固结理论及其合理性验证

党发宁，宋靖宇，周 玫，张 乐

（西安理工大学 西北旱区生态水利国家重点实验室，西安 710048）

1 研究背景

饱和土体的固结，是指土体受到外荷载后土骨架压缩变形、内部孔隙收缩、孔隙水排出的过程。描述饱和土体固结过程的Biot固结理论，将渗透系数视为常数。而在实际固结中，土体骨架受压变形，孔隙会随着土体内部水的排出而减小，渗流通道也会逐渐减小，致使表征土中水运动能力的渗透系数相应减小，特别是在淤地坝、尾矿库固结分析中，渗透系数随固结过程的变化不容忽视。

对于土体固结的非线性渗流研究，学者们做了很多贡献。朱泽明等[1]将渗透系数和孔隙比拟合为半对数关系，在Biot固结理论中考虑了渗透系数的非线性变化。王华敬等[2]用概率方法描述了土体的非均质特点，将固结沉降时土体的弹性模量和渗透系数视为随机变量，建立了Biot平面固结随机方程。Luo等[3]根据黏弹塑性本构关系和流变理论，考虑渗透系数和流变因素的动态响应，构建了基于Biot固结理论的三维耦合模型。Ai等[4]在Biot固结理论中考虑渗透系数的各向异性，利用所提出的传递矩阵推导出多层土的整体刚度矩阵方程，并通过拉普拉斯-傅里叶变换求出解析解。邓岳保等[5]在Biot固结理论的基础上引入了Hansbo非达西渗流定律，应用加权残数法进行三维固结方程的推导。刘忠玉等[6]用Hansbo渗流表达式代替达西定律来描述软土地基固结过程中的渗流状态，以此修正了Biot固结理论，分析了Hansbo渗流参数对孔压和沉降的影响。刘乐飞等[7]考虑了渗流过程中的起始比降，建立了基于Biot固结理论的非达西流固结方程。雷长春[8]在Slepicka的基础上认为孔隙水的渗流为指数渗流模型，用渗流速度与水力梯度之间的指数关系修正了Biot固结方程。梅海峰[9]根据渗透系数经验公式和孔隙比与压缩量的关系，结合太沙基一维固结公式，得到一维变渗透系数固结计算式，且在尾矿库渗流计算中得到了较好的应用。徐进等[10]基于Biot固结方程考虑了渗流各向异性和土体变形的横观各向同性，采用Galerkin法推导了其半解析解。李莎等[11]在Biot固结方程中引入了柯森-卡门渗透系数公式，建立了考虑固结过程中渗透系数变化的固结模型，依此研究了抽水井井壁附近渗流场的变化。沙莎等[12]在提出的应力-渗流-损伤耦合模型中，认为当单元不发生损伤时，渗透系数与有效应力呈指数关系；当单元产生损伤时，认为渗透系数与裂缝宽度有关。高俊等[13]通过考虑软土的初始固结状态及压缩曲线方程，将渗透系数表示为孔隙比、进一步表示为压缩应力的函数，实现了在固结过程中渗透系数随应力状态的变化。朱炳等[14]通过离散元（DEM）和计算流体力学（CFD）相耦合的方法，模拟了细颗粒在土体中迁移时发生的渗流淤堵以及出入口细颗粒流失现象，并且分析了土体内部局部水力坡降和局部细粒含量分布规律。可见，很多学者已经关注到固结过程中渗透系数是非线性变化的，但关于变渗透系数固结理论及其算法的研究仍鲜有报道。

本文利用土体固结时体应变和孔隙比之间的关系，结合达西渗透系数公式以及柯森-卡门公式，推导出能够反映固结过程中随孔隙及体应变变化的渗透系数公式，并将其应用于Biot固结公式中，得到可考虑固结过程中孔隙比与渗透系数变化的变渗透系数固结方程，在COMSOL中实现了求解计算，并通过工程实例进行了计算验证。

2 变渗透系数固结理论

2.1 Biot固结理论Biot固结理论利用弹性力学基本平衡方程，结合有效应力原理，通过位移法得出了平衡方程：

式中：G为剪切模量，G=E/2（1+μ），E、μ分别为弹性模量和泊松比；为拉普拉斯算子；ux、uy、uz分别为x、y、z方向的位移；p为孔压；γ为土体的重度。

此时三个方程无法求解四个变量，因此补充渗流连续方程：

式中：k为渗透系数；γw为水的重度；t为时间。

联立式（1）和式（2）则是Biot固结方程。尽管Biot固结理论公式是由弹性理论推导得来，但在数值分析的时步计算时，采用相应的切线模量，可以考虑土体非线性特性。

2.2 考虑固结过程中孔隙比变化的渗透系数公式已知土体受压时孔隙比e和体应变εv的关系：

式中：e0为土体的初始孔隙比；εv=-（εx+εy+εz）为土体的体应变。可以通过Biot固结方程有限元求解过程的每一求解步获得，依据式（3）就可解出Biot固结方程每一求解步的实时孔隙比。

下面仿照文献[13]用初始孔隙比和初始渗透系数修正达西渗透系数公式和柯森-卡门渗透系数公式，使其可反映固结过程中孔隙比的变化。

达西渗透系数经验公式：

式中：λ为临近颗粒影响系数；β为土颗粒体积系数；d为土颗粒直径，cm；μ为流体动力黏滞系数，g·s/cm2。

柯森-卡门渗透系数公式：

式中：c2为颗粒形状及水的实际流动方向有关的系数，约0.125；ρw为水的密度，g/cm3；s为土颗粒比表面积，cm-1。

假设土体在压缩过程中土颗粒粒径、比表面积以及水的重度等基本物理参数为常数，则由式（4）和式（5），可以将它们用初始渗透系数k0和初始孔隙比e0来表示为：

再将式（6）、式（7）反代入式（4）、式（5），进而得到土体压缩过程中，采用初始渗透系数 k0和初始孔隙比e0表示的渗透系数与固结过程孔隙比的关系式：

将式（3）代入式（8）和式（9）：

将εv与位移的关系式代入得：

式（12）和式（13）为可反映固结过程中孔隙比变化的渗透系数公式。

2.3 变渗透系数固结理论将式（12）或式（13）代入式（2）中，并联立式（1）和式（2）得到本文考虑固结过程中孔隙比与渗透系数变化的变渗透系数固结方程。可以看出，此时的渗透系数是用位移（即ux，uy，uz）表示的函数，而这三个位移参量又为固结方程的未知量，故将该渗透系数函数代入固结方程中，不引入新的未知量。因此，采用上述方法就可在Biot固结理论中考虑渗透系数随固结过程的变化。

3 分析和讨论

3.1 两种修正渗透系数公式分析由式（10）、式（11）可知，渗透系数作为体应变的函数，随着固结压缩的过程而改变。假定初始渗透系数k0=1.0×10-4m/s，初始孔隙比e0=1.0，考证渗透系数随体应变的变化规律，如图1所示。由图1可知，随着体应变（绝对值，下同）的增大，两种方法计算的渗透系数均以不同速率减小，相比之下，柯森-卡门渗透系数下降的速度比达西渗透系数下降的更快。当体应变增大至0.1（即10%）时，孔隙比减小至0.8，减小20%，达西渗透系数减小至7.11×10-5m/s，柯森-卡门渗透系数减小至5.69×10-5m/s，分别降低了28.9%和43.1%；当体应变增大至35%时，孔隙比减小至0.3，减小70%，两种渗透系数依次分别为1.38×10-5m/s和4.15×10-6m/s，分别下降86.2%和95.9%。可见，随着土体的压缩，渗透系数随体应变降低的幅度相当大，且这种渗透系数的非线性衰减是随土体固结过程逐渐增长的，符合多孔介质受压变形、孔隙逐渐减小的客观现象。

图1 渗透系数与体应变之间的关系

3.2 二维工况分析和讨论采用文献[13]计算模型，欠固结黏土层厚5 m，初始孔隙比为1.19，初始渗透系数为1.82 cm/a，压缩系数为0.00023 kPa-1。采用堆载压实固结法，如图2所示。利用COMSOL多场耦合有限元软件实现变渗透系数固结理论计算，考察渗透系数变化前后在不同荷载下的软土固结时沉降和超静孔压的变化，分析不同变形条件下变渗透系数对固结过程的影响。

图2 模型简化图

在COMSOL中建立长40 m、高5 m的二维模型。应力/位移边界条件设为：上边界为荷载边界条件或位移边界条件；左右两侧为法向约束；下边界为固定全约束。流体边界条件设为：上边界排水；左右及下边界不排水。材料属性中除上述取值外，泊松比取0.3，杨氏模量利用压缩系数计算得9.522 MPa，密度取1900 kg/m3。

渗透系数计算公式的实现依赖于COMSOL软件丰富的接口。固结过程中的渗透系数由式（12）、式（13）确定，为初始固结参数e0、k0以及体应变 εv的函数。在 COMSOL软件中，可以调用内部函数solid.evol来获取体应变值，在 “全局定义”里将其定义为变量的表达式，并代入渗透系数的表达式中，即可实现上述公式计算。

3.2.1 不同荷载下修正渗透系数对固结过程的影响 分析不同荷载作用下修正渗透系数对固结过程的影响，在上边界设置荷载边界条件，荷载分别为50 kPa、100 kPa、400 kPa、800 kPa、1000 kPa，分析固结过程中沉降量和超静孔压的变化，计算结果如图3、图4以及图5所示。图例中 “D、K分别表示考虑孔隙比变化的达西渗透系数和柯森-卡门渗透系数的计算结果；B、FDB、FKB分别表示渗透系数不变、修正达西渗透系数公式和修正柯森-卡门渗透系数公式的计算结果。后缀数字表示所受压力，单位kPa。为了方便表述，用修正的达西渗透系数公式计算的结果称之为方案D，用修正的柯森-卡门渗透系数公式计算的结果称之为方案K。

图3 两种修正系数在不同荷载作用下随时间的变化趋势

图4 不同荷载下的沉降计算结果

图5 不同荷载下超静孔压计算结果

图3是黏土体中心点处两种修正系数在不同荷载作用下随时间的变化曲线。由此看出，无论是达西渗透系数还是柯森-卡门渗透系数，均随上覆荷载的增大而减小，且荷载越大，渗透系数的减小幅度也越大。这是因为荷载越大，土体沉降变形也越大，致使土体体应变增大，根据式（3）可知孔隙比也随之减小，故渗透系数相应减小。这也符合固结压密过程中随着土体体积的压缩导致孔隙减小，从而使得渗流通道逐渐减小的现象。同时，对比达西渗透系数和柯森-卡门渗透系数，相同荷载状态下同一时刻达西渗透系数小于柯森-卡门渗透系数，且二者的差值会随上覆荷载的增大而变大，此系式（8）、式（9）形式所定。当荷载为50 kPa时，固结稳定时的达西渗透系数和柯森-卡门渗透系数分别为1.795 cm/a和1.764 cm/a，相比初始渗透系数分别下降1.36%和3.10%；当荷载为1000 kPa时，达西渗透系数和柯森-卡门渗透系数分别为1.452 cm/a和1.246 cm/a，相比初始渗透系数分别下降20.20%和31.52%。可见，随着固结沉降的发展，尤其是在荷载较大时，渗透系数的变化是相当大的。

图4为透水面上不同荷载作用下沉降量的计算结果，图5为土体底部相同工况所计算的超静孔压消散曲线。由图4、图5可以看出，固结过程中同一时刻相同荷载作用下的沉降量大小关系为：原方案＞方案D＞方案K，而超静孔压的规律相反。原因在于，当位移增大后，渗透系数降低，同时刻下的超静孔压消散变缓，从而在平衡方程中使得下一时刻的位移增量减小，进而渗流连续方程中位移增量随时间的变化率减小，并且会导致下一步渗透系数的继续降低，往复循环，彼此耦合。因此，在固结过程中考虑渗透系数变化计算的沉降量和沉降速率均小于未修正的（对于超静孔压变化规律与之相反），并且修正的渗透系数只影响沉降速率和超静孔压消散速率，并不影响最终沉降量，超静孔压最终也会完全消散，故沉降曲线和超静孔压线会在固结过程中与各自未修正的曲线出现差别，但最终将会趋于一致。对应于土体固结过程的解释为：饱和土的固结过程是土骨架变形、孔隙水排出的相互转化过程，土体压缩会导致孔隙比减小，渗流通道减小，土体排水能力降低。当渗流通道受阻后，根据土体有效应力原理，未消散的超静孔压会承担一部分总应力，土骨架的变形将会因超静孔压的消散滞缓而变慢，即渗透系数越小，超静孔压消散越慢，沉降速率越慢。

图6为土体竖向不同荷载下固结一年后孔压（含静水压力）沿深度方向的分布状态。可以看出，在固结一年后，荷载较小时的孔压沿垂向分布已经接近静水压力线，而荷载较大时的孔压正处于消散状态。在相同荷载工况下，原方案所计算的孔压垂向分布均小于方案D和方案K所计算的结果，并且该差距随着荷载的增大而越发明显。当荷载为400 kPa时，较原方案而言，方案D、方案K所计算的孔压差量（深度为2.5 m）分别为9.8%和16.9%；在1000 kPa时为38.2%和62.3%。这表明了在变形较大的工况中，考虑渗透系数变化所计算的孔压在固结过程中体现出的较大差异。

图6 固结一年后孔压沿土体深度方向变化图

图7为土体同一荷载作用下不同时刻点孔压（含静水压力）沿深度方向的分布状态。在800 kPa荷载工况下，孔压在固结初期沿深度方向呈现非线性分布，深度越深孔压变化越小，随着固结时间的增长，该孔压逐渐消散为静水压力，非线性特征也随之消失。以深度为2.5 m处进行分析，对比不同时间修正前后的孔压值，修正后与修正前的孔压差值随固结时间的发展呈现出先增大后减小的趋势，相比于未修正理论所计算的孔压，方案D的增量依次为：3.2%、13.0%、26.4%和0.1%；方案K的增量依次为5.2%、21.7%、43.3%和0.4%。这也进一步体现了修正的渗透系数对固结过程的影响。

图7 800 kPa下不同时间点孔压沿深度方向变化图

3.2.2 不同应变下修正渗透系数对固结过程的影响 分析不同应变下修正的渗透系数对超静孔压的影响，依次在上边界设置初始位移为-0.25 m、-0.5 m、-1.0 m、-1.25 m（分别对应5%、10%、20%、25%的应变）的边界条件，其他计算条件保持不变，透水面上超静孔压的计算结果如图8所示。

图8 不同应变下超静孔压计算结果

由图8可以看出，考虑渗透系数变化所计算的超静孔压大于同时刻未经修正的计算结果，且修正前后超静孔压之间的差值随着应变的增大而增大。例如当时间为0.5年时，原方案、方案D和方案K计算的超静孔压相对于初始孔压下降幅度如表1所示。可以看出，半年后原方案的孔压消散均超过90%以上，而方案D和方案K计算的孔压消散幅度均低于原方案计算结果。当应变为25%时，原方案计算的超静孔压消散率为90.26%，超静孔压值为8.76 kPa；方案 D计算的超静孔压消散率为71.49%，超静孔压值为24.33 kPa，为原方案计算结果的2.8倍；方案K计算的超静孔压消散率为58.78%，超静孔压值为35.21 kPa，为原方案计算结果的4.0倍。此外，当应变较小时（如5%），原方案、方案D和方案K所计算的超静孔压值分别为0.88 kPa、1.10 kPa和1.37 kPa，修正后的超静孔压虽然也有消散滞缓的现象，但差值并不明显，因此本文的修正方法建议在变形较大（如10%以上）的土体中使用，并且应变较大的土体中孔压相差非常明显。

表1 t=0.5 a时超静孔压下降幅度

图9是应变为25%时不同时间点上孔压沿深度方向（土体垂向轴线）的变化线。相比于荷载边界条件，位移边界条件作用下的孔压垂向分布规律更接近线性，且随着时间的增长线性特征越明显，直至固结完成所有曲线都趋近于静水压力线。同样也可以看出修正前后的渗透系数对固结过程中的孔压消散的影响，规律同图7类似。

图9 25%应变时不同时间孔压沿深度方向变化趋势

4 变渗透系数固结理论的合理性验证

采用参考文献[15]的试验数据来验证本文所提理论的合理性。文献中试验所用土样为淤泥质粉质黏土重塑样，颗粒比重为2.71，干密度1284 kg/m3。通过一维固结试验所得土体基本物理参数如表2所示。三维试验采用常规三轴标准试样（高80 mm、直径39.1 mm）。

表2 文献[15]试样固结试验结果

在COMSOL软件中进行数值计算，其中泊松比μ取0.47，弹性模量由压缩模量 Es换算而来，即E=βEs，其中 β=1-（2μ2）/（1-μ）。边界条件模拟常规三轴试验，固体力学边界为：圆柱底面采用固定约束，顶面、侧面受均压荷载；流体力学边界条件设置上下两面排水。固结计算结果如图10所示。其中 “T”表示试验数据，其它符号同上。

图10 不同固结应力下的三种方案计算结果与试验数据对比

图10为原方案、方案D和方案K的计算结果以及文献[15]的试验数据。可以看出，在固结过程中同一时刻下固结度从大到小的排序依次为：原方案＞方案D＞方案K＞试验结果。相比于方案D和方案K，原方案固结度计算偏离试验数据最大，而考虑了土体变形对渗透系数影响的计算结果更接近实测数据，并且修正的柯森-卡门渗透系数公式所计算的更接近试验结果。

图11为不同固结应力下原方案、方案D和方案K所计算的超静孔压消散曲线。在固结初期，将超静孔压略高于固结压力的现象称为曼德尔效应。根据汪江[16]对曼德尔效应与渗透系数的关系研究可知，渗透系数对曼德尔效应的影响体现在其发生的时间上，渗透系数越小，发生的时间越推迟。本文计算所得的渗透系数是与体应变相关的，固结初期的体应变变化量较小，故三种渗透系数差别不大，计算结果均发生了曼德尔效应，但在文献[15]中无论是试验还是其三维WG固结理论中均未出现，可见本文在Biot固结理论中考虑渗透系数变化的方法更为合理。

图11 不同固结应力下超静孔压计算结果

5 结论

本文建立了考虑孔隙比与渗透系数变化的变渗透系数固结方程，在COMSOL中实现了求解计算，通过算例分析和与试验结果的比较，得出以下研究结论：

（1）针对Biot固结理论的渗透系数在固结压缩过程中为常数这一不合理现象，利用固结过程中的体应变修改孔隙比，再利用孔隙比与渗透系数的关系修正达西渗透系数、柯森-卡门渗透系数公式，得到可反映固结过程中孔隙比变化的渗透系数公式，并将其应用于Biot固结理论，得到考虑固结过程中孔隙比与渗透系数变化的变渗透系数固结方程，改善了Biot固结理论中渗透系数不随压缩过程变化这一不合理现象。

（2）通过两个案例分析了运用两种修改的渗透系数公式计算的固结结果，均比未考虑渗透系数变化的方案更接近试验结果，用修改的柯森-卡门渗透系数公式所计算的结果更接近真实解；并且地基的压缩变形越大，原方案与修正方案的计算结果差值越大，修改的必要性越大。

（3）基于变渗透系数固结理论的计算能够更为合理的解释曼德尔效应，虽然在固结初期计算的超静孔压相差甚小，但修正方案能够考虑固结过程中孔隙比以及渗透系数的变化，其结果更接近真实解。