Research on soil pore structure extraction model for microfluidic chips

LiuWu FU; YiDing WU; XiaoYu LIANG; HongXiang ZHOU

doi:10.6038/pg2024II0039

Progress in Geophysics >

2024 , Vol. 39 >Issue 6: 2483 - 2492

DOI: https://doi.org/10.6038/pg2024II0039

Research on soil pore structure extraction model for microfluidic chips

LiuWu FU ,
YiDing WU ,
XiaoYu LIANG ,
HongXiang ZHOU ^,*

Expand

College of Metrology Measurement and Instrument, China Jiliang University, Hangzhou 310018, China

Received date: 2024-01-26

Online published: 2025-01-14

Copyright

Fold

Abstract

Due to the opacity of soil, it is difficult to observe the fluid movement and ecological processes inside soil pores in detail, which has always been a difficulty and hot spot in research. Microfluidic chips developed in recent years provide a visualization platform for pore-scale observation. In order to depict a more realistic soil pore structure on the microfluidic chip, this study constructed two two-dimensional soil pore structure models based on the Monte Carlo algorithm and the watershed algorithm. By performing seepage simulation within the structure using the lattice Boltzmann method, the permeability and two-phase flow saturation of the structure were analyzed. Research results show that the pore structure generated by the Monte Carlo algorithm is closer to the real soil sample in terms of pore parameters and permeability, while the pore structure generated by the watershed algorithm is closer to the soil sample in terms of saturation. In addition, the goodness of fit of the multiple linear regression analysis between the simulation results of permeability and saturation and the pore parameters is 0.96 and 0.99 respectively. The hydraulic characteristics of the model are highly correlated with the pore parameters. Therefore, the soil microfluidic structure extraction scheme proposed in this article is better than the traditional cylindrical scheme and is more in line with the actual soil structural characteristics.

Key words： Pore structure; Microfluidic chip; Lattice Boltzmann method; Permeability; Two-phase flow saturation

Cite this article

LiuWu FU , YiDing WU , XiaoYu LIANG , HongXiang ZHOU . Research on soil pore structure extraction model for microfluidic chips[J]. Progress in Geophysics, 2024 , 39(6) : 2483 -2492 . DOI: 10.6038/pg2024II0039

0 引言

土壤作为一种农业生产资源，其内部的水分布影响着土壤中的物质运移和能量存储，对于作物的生长有着重要的作用(Kravchenko et al., 2019)，所以在孔隙尺度下的土壤水分运动研究对于农业资源的合理调配和根际生态的机理阐释具有重要意义.为了直观观察土壤孔隙内的活动，目前可采用透明材料制作的土壤微流控芯片来模拟土壤内部的多孔结构，基于芯片的精密加工和显微观察技术，直接观测水分在孔隙结构内部的流动情况.对于这类模拟土壤的微流控芯片，现在被学者们称为土壤微流控芯片(Soil-on-a-Chip，SOC).目前对于土壤微流控芯片的刻制主要是在二维层面上(Aleklett et al., 2018)，所以将三维土壤孔隙结构压缩到二维平台上的同时，又能尽可能多地保留土壤孔隙特征，在SOC内部构建合理的土壤多孔结构是一个研究热点.

近年来，有研究提出了用均匀圆柱阵列分布的多孔结构(Scheidweiler et al., 2019)来设计SOC，用于深入探索土壤微生物群落中相互作用的机制，与之相似的还有均匀圆柱交错分布的排列方式来模拟土壤孔隙网络.这两种设计方法关注了土壤孔喉半径这一参数，通过调整微柱直径大小、微柱间距、排布规律等特征来对孔隙率以及孔喉半径等参数进行拟合，在几项研究中都较好地发挥了土壤孔隙结构的作用，是一种较为理想化的构建方法.在此基础上进行改进的是随机半径的圆柱随机分布的构建方式，这种方式的优点是在圆柱构建方法的基础上较为简便地生成各向异性的孔隙结构(朱瑜明等，2023).除了用圆柱的方法构建孔隙网络，还有基于泰森多边形原理的构建方法，Wolf等(2022)利用泰森多边形建立了双孔隙度微观模型用于研究研究碳酸盐中的流动.为了使泰森多边形的几何形状更加接近土壤颗粒的不规则形状，有研究利用改进的算法，选取了伸长度、圆形度、圆整度、规则度来表征颗粒的形状，进行了改进，用该模型来模拟了土壤物理结构(Aufrecht et al., 2022)，来对片上的植物根际进行研究.

目前主流的SOC构建方法是圆柱均匀排列法，对于孔隙度较小的土壤只利用了孔喉参数来拟合，是缺少基于真实土壤数据的SOC构建方法；而泰森多边形的构建结果受到控制点的影响，如果控制点的分布不均匀，会导致泰森多边形的形状不准确.本文基于土壤的CT扫描图像和孔隙结构参数，分别通过分水岭算法和蒙特卡洛算法，生成能够刻制在SOC上的二维孔隙网络，并利用格子Boltzmann方法(LBM)对提取结构进行渗流模拟.本文从孔隙参数、渗透率和饱和度方三个方面，将提取结构与圆柱结构和真实土壤结构进行对比验证，来测试不同结构与真实土壤的拟合程度.此外，本文还在孔隙尺度下，分析了不同的结构参数对水力特性的影响.

1 材料与方法

1.1 土壤孔隙结构参数提取

试验土壤样品是采自浙江省温州市(28°11′23″N, 121°1′10″E)的水稻土，采用X射线断层扫描法，在CT扫描机中采集土壤图像，像素分辨率为26 μm.本文采用灰度分割的方法，将土壤孔隙部分与土壤颗粒部分分割开来.利用处理好后的图像切片进行三维重构，来得到真实土壤孔隙，如图 1所示.其中，红色部分代表土壤固体，白色部分代表土壤孔隙.

显示原图|下载原图ZIP|生成PPT

图1 土壤孔隙示意图

Fig 1 Soil pore diagram

对于土壤孔隙的孔隙参数，本文采用最大球算法进行提取(Dong and Blunt, 2009)，该算法认为孔隙空间是各个不同直径的球体和圆柱相结合，每个孔的直径由该孔可以占据的最大球体的直径得到.最终，土壤孔隙结构通过上述的最大球算法简化为几何模型并在ParaView平台上进行成像，如图 2所示，其中孔隙结构的孔隙度和平均孔喉直径分别为25.13%和160.35 μm.

显示原图|下载原图ZIP|生成PPT

图2 土壤孔隙模型

Fig 2 Soil pore model

1.2 孔隙结构提取模型

(1)均匀圆柱阵列的构建

由于土壤样品孔隙度较小且平均喉道半径与微通道渗流的关联度要高于孔隙率和平均孔隙半径(佘冬立等，2023). 所以本文基于土壤样品的平均喉道半径构建了均匀圆柱整列，平均孔径为160 μm，孔隙率为66.88%, 如图 3所示.其中，蓝色部分为孔隙区域，红色部分为固体区域.

显示原图|下载原图ZIP|生成PPT

图3 圆柱结构的设计

Fig 3 Design of cylindrical structures

(2)基于分水岭算法的结构设计

分水岭算法是图像处理中的一种图像分割方法.其计算流程是将灰度图像中的所有像素按照灰度值进行分区，并在各区内设定一个阈值，让阈值从最小灰度值开始增长，增长过程中会碰到领域的像素，读取其灰度值与设定的阈值进行比较，若小于阈值则被同化，若高于阈值则设置分水岭，这样就对这些邻域像素进行了分类，进而完成分割(Gaol et al., 2020).由于分水岭算法可以处理任意形状的孔隙，自适应地对图像进行分割，适用于复杂的孔隙结构.所以本文将利用该算法进行孔隙设计.

结构提取的具体步骤可总结如下：①将土壤CT切片图像每层叠加，如图 4a所示.其中，黑色区域为孔隙部分，白色区域为固体部分；②将每层对应像素点位的值进行叠加，再除以层数，即可得到颗粒密度矩阵，矩阵项的范围从零到一，可视化图像为图 4b.其中，深红色部分表示颗粒密度高的区域，深蓝色部分代表颗粒密度低的区域；③将颗粒密度矩阵以灰度图像的形式生成，如图 4c所示，用分水岭算法来进行分割；④最后生成的分水岭孔隙网络，如图 4d所示.

显示原图|下载原图ZIP|生成PPT

图4 基于分水岭算法的孔隙结构提取步骤示意图

(a)二值图像的堆叠示例；(b)二值图像堆叠生成的颗粒密度图；(c)灰度图展示的颗粒密度图；(d)利用颗粒密度灰度图分割的孔隙结构.

Fig 4 Schematic diagram of pore structure generation based on watershed algorithm.

(a) Example of stacking binary images; (b) Particle density map generated by stacking binary images; (c) Particle density map displayed in grayscale image; (d) Pore structure segmented using particle density grayscale map.

基于上述的结构提取过程，利用土壤样品的平均孔径参数来调整微通道的宽度为160 μm，同时利用孔隙率参数，生成与土壤样品相近的孔隙率，如图 5a所示.而仅考虑平均孔喉半径的结构如图 5b所示.其中，黑色部分为孔隙区域，白色部分为固体区域.

显示原图|下载原图ZIP|生成PPT

图5 基于分水岭算法生成的孔隙网络

(a)考虑土壤孔隙率和平均孔喉半径的结构图；(b)仅考虑平均孔喉半径的结构图.

Fig 5 Pore network generated based on watershed structure

(a) Structural diagram considering soil porosity and average pore throat radius; (b) Structural diagram considering only average pore throat radius.

(3)基于蒙特卡洛算法的结构设计

蒙特卡洛算法的步骤主要是构造随机的概率以及从已知的概率中分布抽样.本文所用的已知概率为配位数(孔隙空间中与单个独立孔隙相连的喉道数量)的概率分布(Katuwal et al., 2015).配位数由上述的孔隙网络模型中提取得到.

基于蒙特卡洛结构生成的具体步骤可总结如下：①初步指定吼道的长度，在选定的区域内生成等间距的点阵；②每个点的配位数基于样品配位数的概率分布随机生成，通过配位数的赋值来确定该点连通其他点的个数；③对比生成的孔隙网络模型和目标孔隙率来验证生成的模型是否合理.若生成的孔隙率和目标孔隙率相近，则成功生成；反之则返回步骤1修改孔喉长度.该算法流程如图 6所示.

显示原图|下载原图ZIP|生成PPT

图6 蒙特卡洛算法流程图

Fig 6 Monte Carlo algorithm flow chart

同样的，基于上述的结构生成过程，将微通道的宽度设为160 μm，同时利用孔隙率参数，生成与土壤样品相近的孔隙率，如图 7a所示.而仅考虑平均孔喉半径的结构如图 7b所示.其中，黑色部分为孔隙区域，白色部分为固体区域.使用该方法构建出的孔隙模型除了能参考土壤样品的孔喉半径和孔隙率这两种孔隙参数，还能将配位数引入到模型的构建中，在结构参数上更加靠近真实土壤.

显示原图|下载原图ZIP|生成PPT

图7 基于蒙特卡洛算法生成的孔隙网络

(a)考虑土壤孔隙率和平均孔喉半径的结构图；(b)仅考虑平均孔喉半径的结构图.

Fig 7 Pore network generated based on Monte Carlo algorithm

(a) Structural diagram considering soil porosity and average pore throat radius; (b) Structural diagram considering only average pore throat radius.

1.3 孔隙结构参数统计

本文生成的二维孔隙结构是用于微流控芯片内部孔隙微通道的刻画，并最终制作成土壤微流控芯片，如图 8所示.所以本文将以真实土壤结构作为标准，来研究生成的各种二维孔隙结构.

显示原图|下载原图ZIP|生成PPT

图8 土壤微流控芯片模型结构示意图

(a)模型主视图；(b)模型俯视图.

Fig 8 Schematic diagram of soil microfluidic chip model structure

(a) Front view of the model; (b) Top view of the model.

对于上述生成的各种结构，本文通过统计连通性和迂曲度这两个孔隙形态参数的计算来比较各个结构之间的区别.

本文使用了欧拉数来表征孔隙网络的连通性，通过计数孔隙间未互连的数目I减去冗余连接C的数目来计算得到欧拉数E(Ohser and Schladitz, 2009).由欧拉数的计算方法可知，当未互连的孔隙数目小于相互连通的孔隙数目时，整体孔隙的连通性较好，而此时的欧拉数计算为负值.所以，若欧拉数为正值，则表示孔隙的连通性较差；若欧拉数为负数，则表明连通性较好(Zhou et al., 2018).

而迂曲度的定义为水在施加压差方向上的实际行进距离与介质长度之间的比值.目前，迂曲度还不能直接测量.取而代之的是，它可以通过图像分析或流动模拟来计算.在本研究中，通过流体速度场来计算迂曲度值(Muljadi et al., 2016)，计算方法如式(1)所示：

$T_j=\frac{\langle u\rangle}{\left\langle u_j\right\rangle}$

其中〈u〉是局部流速的平均大小，〈u_j〉是速度在平均流动方向上的j分量，可以是x，y，或z方向上.由于水稻土壤的特点，本文主要在z方向上进行计算.

1.4 流动仿真模拟

结合研究机理和计算硬件的需求，本研究选用了LBM进行微通道内的单相饱和流模拟和两相流模拟，区别于传统的计算流体方法，它能够更有效的处理多孔介质的复杂边界.数值模拟过程是在高性能格子Boltzmann代码MF-LBM中实现的.格子排列的方法采用D3Q19模型，即采用三维空间和19个速度方向. 宏观流体密度ρ和速度u可以定义为:

$\rho=\sum\limits_{i=0}^{18} f_i$

$u=\frac{1}{\rho} \sum\limits_{i=0}^{18} f_i e_i$

其中，f是分布函数，e_i是沿着第i个方向的微观速度.文中对上述的孔隙结构分别进行了单相流和两相流动模拟，其中，单项流动的运动方程为：

$f_i\left(x+e_i \delta t, t+\delta t\right)-f_i(x, t)=-\frac{1}{\tau}\left[f_i(x, t)-f_i^{\text {eq }}(\rho, u)\right], i = 0, 1, 2, \cdots, 18.$

在上述的等式中，等号左边为流体粒子的流动过程，等号的右边为流体粒子的碰撞过程，其中平衡分布函数表示为：

$f_i^{\mathrm{eq}}(\rho, u)=w_i \rho(x)\left[1+\frac{3 e_i u}{c^2}+\frac{9\left(e_i u\right)^2}{2 c^4}-\frac{3 u^2}{2 c^2}\right]$

其中f ^eq_i是沿着第i方向在位置x和时间t处的平衡分布函数，w_i是D3Q9模型的权重系数.参数c的数值由式(6)给出，其中δx是晶格间距，δt是时间步长.τ是与运动黏度相关的无量纲弛豫时间，通过式(7)给出：

$c=\frac{\delta x}{\delta t}$

$v=\frac{1}{3}\left(\tau-\frac{1}{2}\right) \frac{\delta x^2}{\delta t}$

而在两项流动的模拟中，则采用双分布函数来表示，使用颜色(红色r和蓝色b)来标记不同的流体，分别由它们自己的粒子分布函数f ^r_i和f ^b_i表示，在位置x和时间t处的混合则表示为f_i(x, t)=f ^r_i(x, t)+f ^b_i(x, t)，流体r和流体b的格子Boltzmann方程均为：

$f_i^{\mathrm{s}}\left(x+e_i \delta t, t+\delta t\right)=f_i^{\mathrm{s}}(x, t)+\Omega_i^{\mathrm{s}(3)}\left\{\Omega_i^{\mathrm{s}(1)}+\Omega_i^{\mathrm{s}(2)}\right\}, {\text{s}} = {\text{r}}, {\text{b}}, i = 0 \cdots 18，$

其中上标s表示流体r或流体b，并且Ω_i^s(1)、Ω_i^s(2)和Ω_i^s(3)分别是负责黏性效应、表面张力效应、以及不同流体分离的碰撞算子.对于密度相同的流体，不需要分别计算每个组分的碰撞算子Ω_i^s(1)和Ω_i^s(2).

本文模拟了空气驱替水的流动.所采用的空气动态黏度为0.018 mPa · s，水的动态黏度为0.97 mPa · s，水和空气的表面张力为72 mN · m^-1.本文使用Latva-Kokko和Rothman(2005)给出的R-K重着色方案，将重着色算子应用于碰撞后的格子.本文将R-K重着色方案中的重着色参数设置为0.9，以形成相对清晰的界面.

2 结果与分析

2.1 仿真模型的验证

本研究中使用了GPU加速的LBM进行两相流模拟，并将模拟的渗透率值与之前的渗透率实验数据(聂绍凯等，2023)进行对比，来验证仿真数据的可靠性.在实验中，当注入微流控芯片中的水达到饱和后，以100 μL/min的速率向模型中连续注入无气水.通过压力传感器读取稳定的压力值，并据此计算微流控芯片的渗透率.实验采用的微流控芯片共设置四组，均为圆柱方阵排列的多孔结构，其直径分别为500 μm、1000 μm、1500 μm和2000 μm.本文导入了实验结构并进行了流动模拟，得出了各组结构的渗透率计算结果，并将仿真计算结果与实验结果进行比较，如表 1所示.

表1 渗透率仿真验证

Table 1 Permeability simulation verification

芯片模型	实验值/(10^-11 m²)	仿真值/(10^-11 m²)
500 μm	5.07	5.45
1000 μm	6.11	6.32
1500 μm	7.51	6.99
2000 μm	8.66	8.61

在得到模拟值后，本文分别利用了平均百分比误差(MAPE)、均方根误差(RMSE)以及最大相对误差(MRE)进行误差分析.其中，MAPE为4.6%, RMSE为0.34，MRE为0.076.上述的误差统计值都相对较小，表明渗透率的仿真结果能较好的匹配实验数据，误差分析如图 9所示.

显示原图|下载原图ZIP|生成PPT

图9 渗透率误差分析图

Fig 9 Permeability error analysis chart

模型的两相流动通过仿真值和实测值的结果对比，在收敛性和可靠性上得到了验证(Chen et al., 2019).此外，本文还进行了空气驱水的两相流实验.实验中，当微流控芯片中的染色水达到饱和时，以100 μL/h的速率向模型中注入空气.当微流控芯片中的两相分布达到准稳态时，使用显微镜拍照并使用图像计算气体的饱和度.实验中设置了六组不同的芯片模型，每组微流控模型均采用随机圆柱排列的多孔结构，本文分别对这六组模型进行了模拟.实验的气体饱和度值与仿真值对比如表 2所示.其中，MAPE为2.53%, RMSE为2.63，MRE为0.052，各误差统计值也相对较小，说明了饱和度的仿真结果也能较好地匹配实验数据，误差分析如图 10所示.

表2 饱和度仿真验证

Table 2 Saturation simulation verification

芯片模型	实验值/%	仿真值/%
组1	93.93	95.84
组2	90.24	94.96
组3	93.92	95.99
组4	92.27	94.79
组5	93.12	92.49
组6	92.25	94.25

显示原图|下载原图ZIP|生成PPT

图10 饱和度误差分析图

Fig 10 Saturation error analysis chart

2.2 孔隙结构参数的对比

土壤的微观结构是影响微通道内流动的重要因素，根据上述的结构参数计算方法，本文对各个孔隙结构参数进行统计.为了方便表示，后面用S表示土壤样品结构；C表示圆柱结构；F1表示考虑土壤孔隙率的分水岭算法结构；F2表示大孔隙率的分水岭算法结构；M1表示考虑土壤孔隙率的蒙特卡洛算法结构；M2表示大孔隙率的蒙特卡洛算法结构.最后，对各种结构的参数统计见表 3.

表3 孔隙结构参数统计

Table 3 Pore structure parameter data

结构名	S	C	F1	F2	M1	M2
孔隙率	25.13	66.88	25.51	67.18	25.64	66.81
连通性	-645	-1023	-137	-403	-56	-627
迂曲度	1.23	1.07	1.48	1.17	1.19	1.13

本文在连通性上使用了Euler数进行表征，值越小表示连通性越高.从表中可以看出，M2最接近土壤样品的连通性，可知在结构生成中采用了配位数的设计，在连通性上能够更加接近所模拟的土壤样品，但M1的连通性却与土壤样品的差距则较大，这是因为蒙特卡洛算法的生成需要较多的数据样品，而由于孔喉宽度、孔隙率以及生成范围的限定，导致在M1样品中的数据点较少，所以产生了较大的偏差.此外，M2结构也较为接近土壤样品，而C的连通性则过大，说明在同等孔隙率以及孔喉宽度的限制下，基于蒙特卡洛算法生成的结构在连通性上拟合效果最好.

迂曲度的数据表明，M1结构更加接近土壤样品.C结构的迂曲度为1.07，说明圆柱结构生成的孔隙通道与直线非常接近.而在同等孔隙率以及孔喉宽度的限制下，基于分水岭算法生成的结构迂曲度最大.

综上所述，在连通性和迂曲度上与土壤样品最接近的分别是M2和M1，说明蒙特卡洛算法生成的孔隙网络在结构拟合上相较于另外两种方法更有优势.

2.3 单相流渗透率的结果比较

基于提取的孔隙网络结构以及构建的孔隙网络模型，本文利用LBM来模拟饱和单相水在微通道的渗流情况，得到了各个结构的渗透率，如表 4所示.

表4 孔隙结构的渗透率模拟

Table 4 Permeability simulation of pore structure

结构名称	渗透率/m²
S	8.43×10^-11
C	2.93×10^-10
F1	4.99×10^-11
F2	3.99×10^-10
M1	1.08×10^-10
M2	4.03×10^-10

从表中可以看出M1的渗透率结果与土壤样品的渗透率最为接近，结构的拟合效果较好，表明蒙特卡洛算法结构能较好的反应原始土壤结构中的渗透率情况.而仅使用孔喉参数的C、F2以及M2结构的渗透率模拟结果与土壤样品相比均偏大，而加上了孔隙率参数的F1和M1是最为接近的两组，表明在构造的结构中孔隙率是拟合渗透率的重要参数.此外，将构造参数相同的三组C、F2以及M2进行比较，发现在相同的孔隙率以及孔喉宽度的条件下，圆柱结构的渗透率最小，而分水岭算法生成的结构以及蒙特卡洛算法生成的结构模拟出的渗透率非常接近.

而多孔介质内的微通道流动会受到孔隙结构的影响.前人的研究也表明，在孔隙尺度上孔隙分布以及孔隙拓扑关系对于土壤的水力特性有着主要的影响(Katuwal et al., 2015).为了探究各项孔隙参数对孔隙网络的作用，本文结合了上述的孔隙率，连通性以及迂曲度来对渗透率进行多元线性回归，来观察各个参数对渗透率的影响.由于各组孔隙参数的数据量级不一致，先对各组数据进行归一化处理后再进行线性回归.所以多项式中的系数大小能反应对结果的影响程度，多项式拟合结果如图 11所示，其拟合优度为0.96，拟合结果较为显著.式(9)为拟合多项式：

$\text { per }=(3.6+30.3 \mathrm{por}+10.9 E-7.3 \mathrm{tor}) \times 10^{-11}$

显示原图|下载原图ZIP|生成PPT

图11 渗透率的拟合结果图

Fig 11 Fitting result diagram of permeability

其中per表示渗透率；por表示孔隙度，E为欧拉数，即连通性；tor为迂曲度.多项式系数表明在相同的平均孔喉半径下，结构的渗透率主要受到孔隙率大小的影响.

上述的拟合结果表明，孔隙参数对渗透率的影响中，孔隙率占主导地位.并且本文发现，孔隙率与渗透率呈正相关，而迂曲度与渗透率呈负相关，这也与前人的微流控实验结果一致(聂绍凯等，2023).Song等(2018)也认为孔隙通道的结构是影响土壤结构渗透能力的主要决定因素，较大的孔隙度可以为结构内的流动提供通道.

2.4 两相流饱和度的结果对比

两相流是在同一空间内两种流体相的流动情况，本文借助了两相流驱替测试来反应多孔介质内的流动特性.在仿真中，本文模拟了空气驱替水的情况，利用了驱替过程中水的饱和度来量化流动情况.仿真过程是对各组孔隙结构注入空气，当注入量达到一个孔隙体积后停止注入，获取的饱和度结果如图 12所示.

显示原图|下载原图ZIP|生成PPT

图12 各组孔隙结构的饱和度结果图

Fig 12 Saturation result diagram of each group of pore structures

从结果中可以看出F1结构的饱和度结果与土壤样品最接近，而其余几组结构的饱和度均小于土壤样品的饱和度，说明在其余几组结构中，空气更容易将水从结构中排出，而基于分水岭算法的结构能较好的反应原始土壤结构中的饱和度情况.同时也发现，在相同的构造方法下，孔隙度较小的其饱和度较高.

同样的，将各个土壤结构参数与饱和度结果相结合，来观察探究各项参数对饱和度的作用.在进行相同的数据预处理后进行多元线性归回，回归结果如图 13所示.

显示原图|下载原图ZIP|生成PPT

图13 饱和度的拟合结果图

Fig 13 Fitting result graph of saturation

从拟合结果图中可以看出，基于孔隙参数的饱和度拟合效果较好，其拟合优度为0.99.数据归一化后生成的拟合多项式为：

$S a=0.3-0.13 \text { por }-0.11 E+0.14 \text { tor }$

多项式结果表明，迂曲度参数更能影响孔隙结构内水的饱和度.而相较于渗透率的回归结果中孔隙率占主导地位，这三种结构参数对饱和度的影响程度差距较小，表明在相同的平均孔喉半径下，两相流的饱和度结果受到迂曲度、孔隙率以及连通性的综合作用.而饱和度拟合结果与渗透率拟合结果相反，孔隙率与饱和度呈负相关相关，而迂曲度与渗透率呈正相关.表明孔隙结构的孔隙率越大，渗透率效果越好，且结构内的水相更容易被空气相排除.表现为孔隙率越大，结构的输送能力越强，这一结果与Wei等(2019)的观点一致，水分移动的难易程度与单位土壤内的孔隙体积密切相关.

3 结论

(1) 在相同的孔喉尺寸条件下，孔隙度是影响结构特性的关键因素.蒙特卡洛算法生成的孔隙结构在孔隙参数和渗透率结果上更加接近土壤样品，分水岭算法生成的孔隙结构在饱和度结果上更接近土壤样品.

(2) 利用孔隙率、连通性以及迂曲度拟合的结果与仿真值吻合程度较高，表明孔隙结构的水力特性与孔隙参数的关联程度较高.多孔介质内微观结构对于控制水分流动起到重要作用.

(3) 拟合多项式的系数表现出的规律与多孔介质实验表现出的规律相符，表明本文的结构模型具有一定的合理性，生成的孔隙网络也能反应多孔介质的结构特性.

感谢审稿专家提出的修改意见和编辑部的大力支持!

References

Publishing order | Descend order by publishing year | Descend order by cited within

Aleklett

, Kiers

E T

, Ohlsson

, et al. Build your own soil: exploring microfluidics to create microbial habitat structures. The ISME Journal. 2018, 12(2 312 319

模态框（Modal）标题

Abstract

Cite this article

0 引言

1 材料与方法

1.1 土壤孔隙结构参数提取

图1 土壤孔隙示意图

图2 土壤孔隙模型

1.2 孔隙结构提取模型

图3 圆柱结构的设计

图4 基于分水岭算法的孔隙结构提取步骤示意图

图5 基于分水岭算法生成的孔隙网络

图6 蒙特卡洛算法流程图

图7 基于蒙特卡洛算法生成的孔隙网络

1.3 孔隙结构参数统计

图8 土壤微流控芯片模型结构示意图

1.4 流动仿真模拟

2 结果与分析

2.1 仿真模型的验证

表1 渗透率仿真验证

图9 渗透率误差分析图

表2 饱和度仿真验证

图10 饱和度误差分析图

2.2 孔隙结构参数的对比

表3 孔隙结构参数统计

2.3 单相流渗透率的结果比较

表4 孔隙结构的渗透率模拟

图11 渗透率的拟合结果图

2.4 两相流饱和度的结果对比

图12 各组孔隙结构的饱和度结果图

图13 饱和度的拟合结果图

3 结论

References