基于FLUENT土壤养分渗流扩散的数值模拟研究

2021-07-30赵崤隆初燕芳李云川

江西农业学报 2021年7期

赵崤隆,初燕芳,李云川,王静

(1.云南农业大学资源与环境学院,云南昆明 650201;2.云南省高校城乡水安全与节水减排重点实验室,云南昆明 650201;3.云南农业大学水利学院,云南昆明 650201)

土壤作为以疏松状态存在于陆地表面并具有肥力的动态自然体,其养分的运动行为以质流和渗流为主,养分随下渗水流下渗融入土壤是植物根系拉伸生长和吸收营养物质的前提。在降雨量大、强降雨集中的湿润、半湿润地区,土壤养分传输、运移作用明显,这一方面使植物更容易吸收养分,另一方面没有被植物吸收的一部分养分会在运动过程中流失，导致土地质量和土壤肥力下降。降雨入渗是养分渗流的重要环节,从渗透机理出发研究土壤的入渗规律和特性,将所得到的渗透曲线应用于所建立的养分运移模型来模拟养分渗流特征,可以为养分有效利用、减少养分流失提供科学依据。目前对于土壤养分渗流规律的研究方法主要是试验和数值模拟相结合的方法,例如：Viotti等建立了一维非饱和土壤中溶质迁移方程,研究了输入不同参数情况下溶质浓度的分布特点[1];张家炳等建立了水分与溶质运移的数学模型，并用该模型模拟了入渗过程中离子的运移规律[2-3];Kumar等通过无网格法分析了二维非饱和多孔介质中溶质迁移问题[4];陶亚奇等采用Voronoi图分裂法来构造与真实土壤相仿的虚拟土壤,并输入到随机行走模型中对土壤溶质迁移进行数值模拟，从而分析了土壤质地对溶质迁移的影响[5-6];孙露露等利用Richard方程建立溶质迁移模型进行数值模拟研究，得出了入渗速度增强可以促进溶质对流和弥散的结论[7-8];魏新平等研究了土壤养分随水分运移过程中水流强度与入渗速度及溶质扩散的关系[9-11]。土壤颗粒在一般情况下处于离散状态(尤其是表层土壤),水分溶液在通过土壤颗粒时,渗流规律的数值分析可以通过有限元软件来进行。本研究模拟了降雨装置试验测定的土壤渗透系数随土壤深度的变化曲线,将Forchhcimer公式修改后得出了适用于描述入渗速率随土壤孔隙度变化的公式,并借助CFD-PDE耦合,采用离散元软件FLUENT的流体模拟功能进行二次开发，探究了磷酸二氢钾溶液在离散的多孔介质土壤中渗流速度和有效扩散系数随土层深度的变化规律。

1 土壤渗透系数的测定

1.1 渗透系数的物理模型

Richard方程能够模拟入渗过程中反应时间与含水率以及导水率的关系。含水量是整个测试过程中最为主要的测试变量,所以该数学模型能够反映整个渗流过程的渗流规律[12-13]。

Richard方程如下:

(1)

式(1)中:θ为土壤体积含水率(cm3/cm3);ψm为非饱和土壤总水势;D(θ)为非饱和土壤扩散率(cm3/cm3);K(θ)为非饱和土壤导水率(cm3/min);t为时间(min)。

1.2 渗透系数测试的机理

渗透系数测定及堵塞装置的工作原理:在试验装置上安装了2个电子水压力传感器和1个超声波流速器,所有传感器通过模数转换器连接电脑。通过水压力传感器可测得试件上、下表面的水头损失,同时通过流速传感器测出水管内水的流速,可将其导入公式(2)。在渗透试验过程中应注意避免土壤试件的侧壁渗漏,且必须在线连续记录渗透系数的变化过程,模拟土壤在一段时间内的渗流过程并且记录动态渗透率变化。渗透系数的测试装置见图1。

图1 渗透系数测试装置以及各组成部件

(2)

式(2)中:k为渗透系数(mm/s);v1为土柱内部的平均流速(mm/s);i为水头梯度;v2为出水管内出口的平均流速(mm/s);Aeff为土柱试件的有效截面积(mm2);Aou为出水管的有效截面积(mm2);L为试件长度(mm); △h为水头损失(mm)。

1.3 渗透率试验设计

为了研究土壤在不同水压下土壤水分入渗的规律，设置3个压力水头(110、160、210 mm)。供试土壤为云南省昆明市梁王山鲜花小镇的红壤土。首先制作符合试验设备的规格为50 mm×100 mm、100 mm×100 mm、150 mm×100 mm的有机玻璃磨具,并在磨具内部涂抹凡士林(起阻止侧壁渗流以及增加粘性的作用，使土壤颗粒与壁面更加紧密，防止取样时因有机玻璃内壁光滑而导致试件分离),将模具插入土壤形成试件,然后去掉周围的土壤取出试件(见图2和图3)。取样日期为2019年12月28日。

图2 供试土壤样品试件

图3 试验整体流程

1.4 渗透系数测定实验结果

如图4所示，3种颜色的曲线分别代表50 mm×100 mm、100 mm×100 mm、150 mm×100 mm试件的渗透系数,渗透系数分别由0.008、0.011、0.0075 mm/s逐渐降低；这些渗透系数的变化具有相同的变化特征,反映了土壤渗透系数随入渗时间增加而变慢的趋势,而50 mm高度的土柱试件因其土层厚度较薄,其土壤渗透系数的变化最为显著。实验还设置了不同高度的压力水头,而在不同高度压力水头下土壤渗透系数随时间的变化大致相同。

入渗速率随土层深度的不同而发生变化,对直径为100 mm,高度为150 mm的试件进行固相和液相的耦合仿真模拟,结果如图5所示，最大入渗速率为0.09 cm/s,最小入渗速率为0.35 cm/s,土壤的入渗速率随土层深度增加而降低。

1.5 分析与校正模型

将孔隙平均渗流速度乘以孔隙度可计算出土壤的渗流速度；通过Forchhcimer公式可以将渗流速度与渗透系数联系起来。由于在初始阶段土壤中含有一定量的空气,所以对Forchhcimer公式中的相关参数进行校正时无法做到准确定位,并且其渗透系数处于一个变化过程,因此很难准确地分析具体的渗透系数的转换,因此可以引进一个关于非饱和度的参数S,使得在整个渗流过程中存在一个非饱和度S，以此对渗透系数进行校正，从而得到了以下基于Darcy-Forchhcimer多孔渗流的渗透系数变化的公式，亦即改进的Forchhcimer公式:

(3)

(4)

式(3)中:J为压力梯度(Pa/m);S为非饱和度(%); μ为流体动力粘滞系数1.005;V为渗流速度(m/s);β为非达西系数(m-1)。

式(4)中：K为渗透率(10-10m2);k为水力传导系数(m/s)；μ为流体动力粘滞系数1.005。

修正后Forchhcimer公式成为一个与孔隙度有关的函数,从而使之更加符合土壤模型的适用性。将渗透曲线与改进后的Forchhcimer公式结合,能反映整体渗流过程中孔隙体积的变化规律。

图4 土壤渗透系数的变化曲线

图5 入渗速率随土层深度的变化曲线

2 养分运移的离散元土壤模型设置

土壤是典型的松散介质,利用离散元模型模拟土壤水化后聚团的形成及破裂是目前较为合理的方法,目前常用离散元法模拟土壤的模型有整合延迟弹性模型(hysteretic spring contact model, HSCM)和线性内聚力模型(liner cohesion model, LCM)。最近提出的EEPA建立云南红土农田土壤模型最为合适,其弹性塑性形变模拟能反映颗粒之间在受力作用下更为紧密的现象,能够表现颗粒之间孔隙压缩的过程和反映不同压实程度下的养分渗透机理。本实验选取的梁王山鲜花小镇的表层原状土在长时间形成过程中受到压实固化以及沉积作用，其颗粒粘结较为紧密,因此在进行实验设计时需要采用DEM中模拟抗压的设计进行一定的压实，进而确保模拟的试件与实验试件保持相同的密度[14-17]。

2.1 颗粒参数的设置

由于在模拟时不好掌握原状土层颗粒的分布情况,因此本实验先将取到的原状土层压碎，然后筛分土壤颗粒粒径,选取2.0、1.0、0.5 mm三种孔径的筛子进行筛分,其质量占比分别为20%、65%、15%。设置红黏土颗粒的相关离散元参数如下:恢复系数0.55,静摩擦系数0.20,动摩擦系数0.10,密度2600 kg/m3,剪切模量1107 Pa,泊松比0.25。

2.2 数学模型

组分运移方程如下:

(5)

式(5)中:U为流速;ρ为密度;Q为某种流体的属性,可以是标量(如浓度C)，也可以是向量(如流速U)。

2.3 物理模型的边界条件和网格划分

在渗流过程中,养分渗流扩散的运动分为水平运动和垂直运动。本文以多孔介质土壤中的养分渗流扩散为研究对象来建立模型,沿X轴正负方向长度为250 mm,沿Z轴负方向的深度为250 mm,在模型中距离Z轴0 mm和80 mm处各设置1条监测线。设计边界条件为:上界面为速度进口,在以坐标原点为中心处设40 mm×40 mm区域,该区域含有10%磷酸二氢钾养分释放点；下界面为速度出口,设为自由出入边界,与Z轴对称的左右侧边界设为压力远场,结构如图6所示。

根据土体试件模型内部的特性,将各计算域划分为四面体结构化和非结构化网格并进行分块。为了能够获得更加精确的数据和观察流场的特性,对质量分数网格进行细化,将试件以2.5 cm为标准尺度进行等距离剖分,初试的迭代控制参数设为默认值,计算区域的网格的质量百分数为95%以上，其中流体网格和介质网格划分的结果如图7所示。

图6 养分渗流扩散物理模型

图7 扩散模型的网格划分

3 数值模拟及结果分析

多孔介质扩散是指气体或液体进入固态物质孔隙的扩散。液态分子在多孔介质土壤中扩散是一种常见的现象,是由于浓度梯度的存在而发生的质量转移过程。数值模拟以长500 mm、宽500 mm、高250 mm的土体为研究对象,在中心点处和距离中心点80 mm处设置两个监测点。结合公式(3)～公式(5)通过计算流体力学CFD软件中的FLUENT进行液固两相的耦合仿真,以模拟磷酸二氢钾溶液渗流过程速度、质量分数和有效扩散系数的变化特征。

3.1 模拟流场内流体速度的变化特征

在溶液下渗扩散过程中两个监测点的速度变化如图8所示,两个监测点的扩散速度随土层深度增加总体上呈下降趋势,监测点1的扩散速度总体上处于0.00500～0.00506 m/s,在Z轴方向-0.02 m处达到峰值，这是因为从释放点自由落体到达上界面处存在始速度，其比扩散介质的阻力大;监测点2距离中心扩散源0.08 m,由于多孔介质的阻力项的存在,其质量分数扩散速度整体小于监测点1且呈缓慢下降趋势;模拟场中磷酸二氢钾溶液的渗流速度在中心点最大,在土壤0.02 m深处达到峰值,并随土层深度增加向下呈辐散式递减。

图8 两个监测点渗流速度随土层深度的变化

3.2 模拟流场内磷酸二氢钾溶液相对质量分数的变化

相对质量分数的变化可以反映溶液浓度在渗流过程中的变化状态。两个监测点磷酸二氢钾溶液相对质量分数的变化如图9、图10所示,其中监测点1的相对质量分数随土壤深度增加呈下降趋势,在Z轴方向0～0.03 m处下降幅度明显,之后下降幅度趋于平缓;监测点2的相对质量分数随土壤深度增加呈上升趋势,上升幅度在Z轴方向0.05～0.15 m处最明显。结合图11磷酸二氢钾溶液相对质量分数云图得其浓度在中部高,并沿土层的纵深方向和横向递减。

图9 监测点1相对质量分数随土层深度的变化

图10 监测点2相对质量分数随土层深度的变化

图11 模拟场内不同时间相对质量分数扩散云图

3.3 模拟流场内磷酸二氢钾溶液有效扩散系数的变化

有效扩散系数表示流体在多孔固体中的扩散状态。从图12可以看出，磷酸二氢钾溶液在监测点1的有效扩散系数先增加后减少,总体处于0.001～0.0012 m2/s，在土壤0.04 m深处达到峰值,之后有效扩散系数逐渐减少。由图13可见,监测点2的有效扩散系数整体上随土壤深度的增加呈上升趋势，在Z轴方向10～20 cm处上升幅度大。结合图14的云图得模拟场内中心点土层下0.02～0.03 m处扩散系数最大,向下和向两侧逐渐减小,最后呈弥散状态。

4 结论与讨论

同高度压力水头下土壤渗透系数随时间的变化大致相同,随时间增加呈下降趋势;在渗流过程中入渗速率随土层深度增加而降低,土壤渗透系数会出现动态下降的过程。拟合出与土壤渗透系数变化趋势一致的渗流曲线,并引入非饱和度参数S,得到修正后的能反映整体渗流过程中孔隙体积变化规律的Forchhcimer公式。

图12 监测点1有效扩散系数随土层深度的变化

修正的Forchhcimer公式与FLUENT结合模拟液相的养分溶液填充土壤孔隙的渗流过程,通过渗流速度、相对质量分数以及有效扩散系数的变化曲线和液相分布的数值模拟流场云图发现,多孔介质土壤中的养分在中心扩散源0.2～0.3 m处的渗流扩散能力最强,且随着土层深度和距中心点横向距离的增加呈辐散性减弱。多孔介质土壤作为阻力项的存在导致养分扩散系数在两个监测点的变化不同。