李云川 赵崤隆 石平平 李靖 南轩 王静

摘要：基于计算流体动力学与离散元（CFD-DEM）耦合的数值模拟方法，结合土壤入渗试验，进行土壤渗透系数随土层深度变化的研究。试验结果表明，随着入渗过程的持续，3、6、9 cm土柱渗透系数分别由0.008 0、0.011 0、0.007 5 mm/s逐渐降低。在液固耦合模型中，通过计算流体动力学（CFD）设置流体相参数，离散元（DEM）模拟土壤颗粒固体相，得出渗流速率与土层深度的定量变化关系式为y=ax2+bx+c（y为渗流速率;x为土层深度;a为常数，等于 -0.078 86;b为常数，等于-0.008 67;c为常数，等于0.092 96）。通过方差分析对模型拟合的回归性方程进行验证，得到相关系数R2=0.950 3，校正决定系数R2adj=0.999 77，都接近于1，表明拟合性好。

关键词：土壤渗透系数;数值分析;液固耦合;CFD;DEM;数值模拟

中图分类号：S152.7+2

文献标志码：A

文章編号：1002-1302（2020）16-0255-05

土壤和水是人类社会生产、生活必不可少的自然资源。土壤水的运动行为主要分为质流和入渗，入渗是一个复杂的过程，无论是植物根系拉力导致的质流作用，还是土壤中的原始水分含量，均来源于入渗过程。降水入渗是土壤水循环过程的重要环节，开展土壤水分入渗规律和特性的研究，可以合理判定地表径流，进行暴雨洪水预警，为制定因土壤侵蚀导致的水土流失治理措施提供依据[1-3]。有关土壤水分入渗的研究可以追溯到1856年，法国水力学家达西（Darcy，1856）提出了达西定律，可反映水在岩土空隙中的渗流规律[4-7]，该试验定律的提出为研究土壤水分入渗理论和入渗过程模型奠定了基础。Green-Ampt入渗物理模型仅适用于在薄积水或土质均匀条件下的土壤水分入渗研究[8]，将土壤在垂直方向上划分为入渗湿润层、入渗饱和层和湿润锋，该模型反映了入渗速率与土壤特性之间的关系。1931年，Richard在达西定律的研究基础上提出了可用于非饱和流土壤水分入渗分析的Richard方程[8-9]，并得出导水率与基质势之间的关系，用于描述土壤水分运动的规律，Richard入渗模型（Richard，1931）的建立是达西定律的继承与发扬。Philip入渗模型（Philip，1957）以Richard入渗模型为基础，描述入渗速率随渗透时间的变化趋势，但该模型存在模拟值与实测值吻合度不高的问题。在众多入渗模型中，Kostiakov入渗模型（Kostiakov，1932）的使用最为简单、便捷，因此在实际研究中被广泛应用，但在利用该模型进行参数率定时，各参数的物理意义表述仍有争议。以上数学模型的创建极大地扩展了土壤水分入渗的研究，根据长期的科学研究可知，土壤入渗率随土层深度的增加而呈现“高- 低-稳定”的变化趋势，但以往的研究很少涉及土层深度对入渗速率变化影响的定量描述[10]，本研究基于计算流体动力学与离散元（CFD-DEM）耦合的数值模拟方法，分析土壤入渗速率随土层深度增加的变化情况，以期获得不同深度土层土壤的渗透系数定量变化规律。

1 渗透率变化的试验

1.1 试件制作

供试土壤为云南省昆明市盘龙区松华坝水源区棕壤土，土壤样品为表层土壤，采样时剔除土壤表面植物残渣和固体石块，捏碎表面大颗粒黏结土块，使得清理后的土壤表面均匀平整，将相应尺寸的有机玻璃圆柱体打入土体取样。供试土样基本理化性质如下：容重、全氮含量、全磷含量、全钾含量、pH值、EC（可溶性盐浓度）值、氧化还原电位（Eh）分别为1.69 t/m3、0.14%、0.21%、0.75%、7.3、0.8 mS/cm、136 mV。

试验所用试件的直径为10 cm，高度为3、6、9 cm，材质为有机玻璃，为圆柱形筒体，渗透率测量仪器采用透水性渗透设备。分别设置直径为 10 cm，高度为3、6、9 cm的圆柱形土柱各3个，共计9个土柱试件，土壤试件如图1所示。

1.2 渗透系数测定和试验装置

渗透系数计算公式如下：

式中：K为渗透系数，mm/s;υ1为土柱试件内部的平均流速，mm/s;i为水梯度，表示土柱试件高度与长度的比值;υ2为出水管内出口的平均流速，mm/s;S1为出水管的有效截面积，m2;S2为土柱内的有效截面积，m2;L为土柱试件长度，mm;ΔH为水头损失，mm。

试验装置（图2）工作原理：透水性渗透测定设备上安装了2个电子水压力传感器和1个超声波流速器，压力水头高度为11、16、21 cm，出水口为自由出水口，传感器、超声波流速器与数据采集仪均连接测量设备，水压力传感器可测得土柱试件上下两端水头损失（ΔH），流速传感器可测得土柱试件内部流速（υ1）和出水管内水流速（υ2），将数据采集仪收集的数据导入公式（1），计算得出渗透系数。

1.3 渗透系数测定试验结果

图3显示了在压力水头高度为11 cm条件下，土柱试件渗透系数随时间的变化规律，3条曲线分别对应不同高度的土柱试件，它们具有相同的线型，均反映了土壤渗透系数随入渗时间增加而呈现由高到低的变化趋势，3、6、9 cm土柱渗透系数分别由0.008 0、0.011 0、0.007 5 mm/s逐渐降低。3 cm高度的土柱试件因其较薄的土层厚度，土壤渗透系数的变化最为明显。试验还设置了高度为16、21 cm 的压力水头，而在不同高度压力水头下土壤渗透系数随时间的变化规律大致相同。

2 入渗研究的数学模型

2.1 Green-Ampt入渗物理模型

改进前的Green-Ampt入渗模型适用条件范围单一，仅适用于土质均匀、入渗初期土壤水分含量低的土壤入渗过程研究，导致该模型的使用有很大的局限性。国内外学者出于拓宽Green-Ampt入渗模型适用范围的目的，对其进行了改进，修正后的Green-Ampt入渗模型可表示为

式中：K为渗透系数，mm/s;Ks为有效饱和导水率，mm/s;θ1为初始含水率，mm3/mm3;θ2为饱和含水率，mm3/mm3;I为累计入渗量，mm;SK为入渗土层湿润锋处的固定吸力。

2.2 土壤水分运移方程

水分在土壤中的运动轨迹可分为水平侧移和纵向垂直运移，Richard模型可以描述土壤水在一维垂直方向上的运动规律，数学方程式可表示为

式中：θ为土壤体积含水率，%;H为压力水头高度，mm;t为渗透时间，s;α为流向与垂直方向夹角，rad，该模型为一维垂直方向的一阶偏导，因此α=0;Kf为非饱和渗透系数，mm/s;S为源汇项。

3 数值模拟

3.1 设置网格和时间步长

根据土柱透水性渗透试验，选取9 cm高度土柱为模拟对象，将土柱均匀划分，以1 cm为尺度标准将土柱等距剖分，每隔1 cm设置1个观测点（图4）。初始时间步长为0.001 s，最小时间步长为0.001 s，最大时间步长为10 s，迭代控制参数采用默认值。

3.2 初始条件和边界条件

入渗过程中土壤水的运移轨迹分为水平和垂直方向，本研究建立模型时只考虑水分在垂直面内的运动。土壤渗透试验从土柱上端连续进水，压力水头设定为11 cm高度定水头，模拟模型上边界选择恒定压力水头边界（constant pressure head），下边界为自由出水边界（free drainage）。

3.3 参数设置

CFD设置：液相为试验用水，选用RNG k-ε湍流模型进行模拟，采用压力进口，自由出口，壁面条件选择增强壁面函数和无滑移边界条件，动量和湍流动能采用二阶迎风格式，采用压力耦合方程组的半隐式（SIMPLE）算法进行求解。

DEM设置：固相为土壤，进出口条件与液相相同，模拟过程为瞬时模拟，时间步长为0.001 s，模拟时间总时长为5 s。

4 模拟结果与分析

4.1 模拟场内流体受到的动态压强和颗粒总能量

根据实际试验条件，设置高度为11 cm的压力水头，选取直径为10 cm，高度为9 cm的土柱试件作为仿真模拟对象，进行CFD-DEM耦合数值模拟分析，通过流体力学软件FLUENT得到CFD-DEM耦合场下的流体动态压强和固相土壤颗粒总能量。

从图5流体受到的动态压强云图可以看出，当水流与土壤接触时，流体动态压强基本保持不变，在入渗过程持续一段时间后，土壤颗粒中的空隙压力会骤然上升，这是由于流体的持续入渗填充了原本土壤颗粒间的空隙，当水分入渗量达到饱和时，流体所受压强不再发生明显的变化，基本保持不变。

从图6可以看出，在渗流初期土壤颗粒总能量最大，这是因为流体从高处落下时携带一定的势能，在冲刷到土壤时产生曳力作用，导致土壤颗粒受力，从而增加土壤颗粒的动能，土壤颗粒总能量也随之增大。随着入渗过程的持续，土壤水分含量达到饱和，颗粒总能量减少，少部分区域的土壤颗粒因土壤吸水产生的黏结作用而呈现总能量增大的变化，但总的来说，颗粒总能量基本保持不变。

4.2 數值模拟耦合场中土壤颗粒相与流体相的分布

如图7、图8所示，结合颗粒模型可以发现，流体在多孔颗粒域整体所占体积极小，且大部分位置的体积分数都小于0.187，同时可以发现，上部的压实程度较小，流体相体积分数较大，而下部压实较严重，体积分数相对较小，这与渗透系数随深度增加而减少一致。

4.3 数值模拟耦合场内入渗速率的变化

进行数值模拟时选择直径为10 cm，高度为 9 cm 的土柱试件为研究对象。每隔1 cm土层深度设定1个监测面，通过FLUENT软件进行液固2相耦合仿真，得出每个监测面流体的平均流速。结合图9、图10可以看出，土壤模型上层孔隙度较大，监测面所监测的平均入渗速率大于下层。

5 模型验证

为了确保模型的准确性和适应性，需要对模型的预测能力进行评估，根据需要对响应面方程进行显著性检验，一般采用相关系数R2和校正决定系数R2adj来了解其逼近程度。

式中：ST=∑ni=1（y1-y），为总平方和;SA为回归平方和;SE为残差平方和，相关系数为完全拟合的度量值，反映了响应面与试验数据的符合程度，其值要在0.9以上才能保证拟合效果。

式中：校正决定系数R2adj为自变量与因变量的相关程度，其值越接近1拟合效果越好。对拟合回归性方程进行方差分析可以得到，渗透速率回归方程的相关系数R2=0.950 3，校正决定系数R2adj=0.999 77，两者都接近于1，表示方程拟合性良好，可认为拟合得到的回归方程能够反映参数之间的关系，可以用来拟合渗流速度与土壤深度的变化关系。得到如下拟合方程：

式中：y为渗流速率，m/s;x为土层深度，m;a为常数项，等于-0.078 86;b为常数项，等于 -0.008 67，c为常数项，等于0.092 96。

6 结论

通过设置11、16、21 cm高度的压力水头，对直径为10 cm，高度为3、6、9 cm土柱分别进行渗流试验，分析土壤渗透系数随土壤深度的变化规律，并基于CFD-DEM耦合进行土壤渗透性数值模拟研究，得出结论如下：

（1）通过相分布云图和渗透系数的曲线可知，入渗过程中土壤深度越低的区域，土壤含水量越高，土壤深度越高的区域，土壤含水量越低，这是由渗透系数随土壤深度的变化决定的。

（2）渗透系数受压力水头高度的影响较小，在11、16、21 cm高度的压力水头下，渗透系数无明显变化。

（3）通过数值模拟分析得出入渗速率随土层深度增加而降低的定量变化关系，可表示为y=ax2+bx+c，其中y为入渗速率，x为土层深度，a为常数项，等于-0.078 86;b为常数项，等于-0.008 67，c为常数项，等于0.092 96。

（4）對拟合回归性方程进行方差分析，得到相关系数R2=0.950 3，校正决定系数R2adj=0.999 77，二者均接近于1，表明拟合度好。

参考文献：

[1]高 岩，胡晓蕾，刘 然. 雨水渗流过程中土壤水分动态变化规律研究[J]. 建筑技术，2019，50（1）：4-7.

[2]杨 帆，武桂芝，冯增帅，等. 大沽河河床非饱和入渗模型研究[J]. 水土保持通报，2017，37（6）：152-156.

[3]王 俊，黄岁樑. 土壤水分特征曲线模型对数值模拟非饱和渗流的影响[J]. 水动力学研究与进展，2010，25（1）：16-22.

[4]李俊烨，苏宁宁，胡敬磊. 基于CFD-DEM耦合的磨粒流微小孔加工数值分析与试验[J]. 农业工程学报，2018，34（16）：80-88，299.

[5]杨庆璐，李子涵，李洪文. 基于CFD-DEM耦合的集排式分肥装置颗粒运动数值分析[J]. 农业机械学报，2019，50（8）：81-89.

[6]祁文燕. 纸坊沟流域坡面尺度土壤水分动态变化及其运移数值模拟[D]. 兰州：兰州大学，2018.

[7]丁海晶，姜 姜，张金池，等. 土壤渗透性的区域变化规律及因子分析[J]. 水土保持学报，2019，33（1）：51-56.

[8]孙丹焱，郑 涛，徐竟成，等. 城市绿地土壤渗透性改良对雨水径流污染的削减效果及去除规律[J]. 环境工程学报，2019，13（2）：372-380.

[9]Moret-Fernández D，Latorre B，Giner M L，et al. Estimation of the soil hydraulic properties from the transient infiltration curve measured on soils affected by water repellency[J]. Catena，2019，178：298-306.

[10]兰简琪，谢世友. 有机复合肥对土壤水分入渗特性的影响[J]. 江苏农业科学，2020，48（5）：280-286.

[11]Solekhudin I. Steady infiltration in heterogeneous soil[J]. Journal of Physics（Conference Series），2018，1108（1）：012030.