利用多组分格子玻尔兹曼方法预测部分饱和多孔介质中的毛细压力

《Advances in Water Resources》：Capillary pressure prediction in partially saturated porous media with a multicomponent lattice Boltzmann method

　　S. Joseph|H. Cheng|J. Harting|V. Magnanimo
荷兰恩斯赫德特文特大学工程技术学院土壤微观力学系主任

摘要
准确量化平均毛细压力对于理解部分饱和多孔介质中的流体流动和保持非常重要。在多组分伪势Shan–Chen格子玻尔兹曼方法中，扩散界面和虚假速度使得毛细压力对流体相的整体定义和采样方式非常敏感。在这项工作中，我们比较了传统的体积平均毛细压力测量方法与基于界面长度加权的局部测量方法，后者可以通过杨–拉普拉斯定律从界面曲率获得，或者通过对界面的局部采样获得。首先，我们讨论了一种确定临界孔径的程序，这被认为是进行准确毛细压力计算所需的关键几何特征，因为这样做可以捕捉到控制体积平均毛细压力的大多数流体–流体界面。然后，我们通过网格收敛性研究以及对局部尺度上杨–拉普拉斯关系的验证来评估毛细压力估计值。在有序多孔介质中，体积平均毛细压力与基于长度加权的局部毛细压力平均值非常接近。然而，在非均匀多孔介质中，两者之间存在差异，因为局部毛细压力分布不均，且不同簇对相体积和界面长度的贡献也不同。使用颜色梯度公式作为基准，我们展示了当仔细识别出整体区域时，Shan–Chen公式能够产生准确的毛细压力预测。

符号列表
符号 描述
声速 物理单位中的声速
最小粒子直径 最小粒子直径
与界面的间隙距离 与界面的临界间隙距离
时间步长 物理单位中的时间步长
格点间距 物理单位中的格点间距
离散速度 第i组分的离散速度
流体–流体相互作用力 第i组分的流体–流体相互作用力
固体–流体相互作用力 第i组分的固体–流体相互作用力
第i组分在第k方向上的粒子分布函数 第i组分在第k方向上的粒子分布函数
表面张力 物理单位中的表面张力
流体组分间的相互作用强度 流体组分间的相互作用强度
固体与非润湿流体间的相互作用强度 固体与非润湿流体间的相互作用强度
网格细化因子 三维域的长度
非润湿流体占据的格点数 非润湿流体占据的格点数
润湿流体占据的格点数 润湿流体占据的格点数
运动粘度 物理单位中的运动粘度
第i组分的伪势函数 第i组分的伪势函数
毛细压力 毛细压力
基于界面长度加权的平均毛细压力 来自界面周围局部采样的平均毛细压力
基于杨–拉普拉斯关系的平均毛细压力 来自杨–拉普拉斯关系的平均毛细压力
非润湿流体的平均压力 非润湿流体的平均压力
润湿流体的平均压力 润湿流体的平均压力
流体密度 物理单位中的流体密度
第i组分的流体密度 第i组分的流体密度
虚拟固体密度 虚拟固体密度
划分整体区域与临界区域的临界占位分数阈值 划分整体区域与临界区域的临界占位分数阈值
第i组分的松弛时间 第i组分的松弛时间
非润湿流体的接触角 非润湿流体的接触角
润湿流体的接触角 润湿流体的接触角
平衡速度 平衡速度
质心速度 质心速度
簇对其相界面长度的贡献 簇对其相体积的贡献

1. 引言
多孔介质指的是材料内部均匀或非均匀分布有孔隙的固体。这些孔隙可以包含一种或多种流体，并且大多数孔隙是相互连通的。如果这些流体不互溶，它们会在界面处形成明显的界面，该界面受界面张力的控制。流体对固体表面的润湿或非润湿行为通过接触角来量化。理解毛细压力在涉及多孔介质中流体流动的多个学科中都非常重要，包括岩土工程、油气勘探、水资源、建筑材料科学、碳储存以及土壤和植物科学（Rahardjo等，2019；Jennings，1987；Vavra等，1992；Schmidt和Slowik，2013；Ioannou等，2003；Tokunaga和Wan，2013；Clothier和Green，1994）。

格子玻尔兹曼方法（LBM）能够在介观层面上模拟流体流动（Harting等，2004；Aidun和Clausen，2010；Chen等，2014）。该方法适用于模拟多孔介质中的流体流动，因为它可以处理复杂的几何形状和固体–流体相互作用，并且适合并行计算（Schmieschek等，2017；Liu等，2016；Martys和Chen，1996；Yang和Boek，2013；Zudrop等，2017）。它可以模拟各种流体系统：单组分、单组分多相、多组分或多组分多相系统。单组分系统只包含一种流体（Cheng等，2019）。单组分多相系统（例如水和蒸汽）同时存在流体及其蒸汽相（Sukop和Or，2004；Delenne等，2015）。多组分系统包含两种不同的流体（例如水和油，Huang等，2007；Sukop等，2008）。多组分多相系统包含不同的流体及其蒸汽相（例如水、空气和水蒸气，Bao和Schaefer，2013）。本研究关注的是多组分类别中的部分饱和系统，其中两种不相溶的流体（润湿和非润湿）共存，并可能占据相同的格点位置。这些流体使用伪势Shan–Chen（SC）LBM进行数值模拟（Zudrop等，2017；Shan和Doolen，1995；Porter等，2012）。

在部分饱和的多孔介质中，毛细压力来源于流体–流体界面和固体–流体界面处的作用力。它取决于孔隙几何形状、饱和度、界面面积和界面曲率（Morrow，1970；McClure等，2016a；Miller等，2019）。除了组分的饱和度外，域内界面的分布对于定义部分饱和多孔介质的状态也是至关重要的（Hassanizadeh和Gray，1993）。因此，使用如LBM这样的基于扩散界面的数值方法获得精确的毛细压力估计值，取决于在多孔结构中精确定位扩散界面和流体整体边界。

多孔介质中的宏观毛细压力可以通过几种平均技术得出，如体积平均、界面平均或边界段平均。体积平均是在流体的整个体积上计算的；界面平均是在界面上计算的；边界段平均是在边界段上计算的，通常在域边界上的压力被明确控制时使用（McClure等，2016）。研究表明，在压力分布不均匀的多孔介质中，这些方法会产生不同的毛细压力（McClure等，2016；Starnoni和Pokrajac，2020）。体积平均简单易行，但表面平均需要仔细处理界面以避免扩散界面伪影。由于多孔介质的内部异质性，基于边界的方法（例如在SC LBM研究中通过入口–出口压力差来确定毛细压力（Porter等，2009；Warda等，2017；Pan等，2004；Vogel等，2005；Ahrenholz等，2008）可能与体积平均结果不同（Nordbotten等，2008）。Delenne等（2015）、Nekoeian等（2018）和Hosseini等（2024）的研究中使用了体积平均方法，尽管他们没有明确描述用于排除扩散界面区域的过滤过程。Li等（2018）和Li等（2019）的工作中描述了过滤过程，其中通过应用固定的密度阈值将格点分类为非润湿、润湿和界面类型。由于界面密度剖面取决于诸如流体–流体相互作用强度和分辨率等参数，固定范围缺乏通用性。这突显了需要一种自适应过滤方法，可以系统地调整阈值，直到毛细压力值稳定。

在具有狭窄通道和大孔隙的非均匀介质中，流体可能会分裂成多个不连续的簇，导致域内的压力场不均匀。此外，大孔隙和小通道的共存使得选择合适的模拟分辨率变得具有挑战性。一些毛细压力研究报道的网格分辨率为每个平均孔隙通道半径约2.7个格点单位（Pan等，2004）或每个最小孔隙半径4个格点单位（Li等，2019）。需要根据扩散界面宽度来确定特征长度并调整网格尺寸，以确保毛细压力的收敛性。

基于上述知识空白，我们提出了一种自适应过滤方法，帮助定位扩散界面和流体整体边界，从而确保只有无虚假速度的整体部分被采样用于毛细压力平均。在毛细条件空间均匀的多孔介质中，单一的体积平均毛细压力可以提供有用的流体配置描述。然而，在具有不连续簇的非均匀多孔介质中，毛细压力可能在不同流体–流体界面间变化，从而产生局部毛细压力的分布。这引发了体积平均毛细压力与非均匀多孔介质中基于界面的毛细压力测量方法之间的关系问题。通过在有序和非均匀多孔介质中进行网格收敛性研究，验证了所提出的毛细压力测量的准确性。最后，我们建立了一种确定临界孔径作为控制非均匀介质中毛细压力收敛性的特征长度的程序，并用扩散界面宽度来表达分辨率标准。

本文的结构如下：第2节介绍多组分SC LBM的理论方面、状态方程、毛细压力计算和分区方法；第3节讨论模拟设置；第4节研究了在简化几何形状（如拉普拉斯测试、静止液滴测试、板间液桥和部分饱和有序多孔介质）中的毛细压力计算，重点介绍了整体区域的识别；第5节将分析扩展到非均匀多孔介质，展示了各种毛细压力测量方法；第6节讨论了研究的范围、应用和局限性；第7节总结了研究。

2. 理论
2.1. 多组分伪势Shan–Chen格子玻尔兹曼方法
LBM的基本量是粒子分布函数（φ）。它描述了流体在空间和时间中的运动。与传统计算流体动力学（CFD）方法（Zawawi等，2018）不同，LBM不直接使用纳维–斯托克斯方程来模拟流体运动。相反，流体被描述为存在于格子上的虚拟粒子集合（图1）。这些粒子移动并碰撞，它们的集体行为决定了流体的宏观性质。

带有Bhatnagar–Gross–Krook算子的格子玻尔兹曼方程描述了第i组分的粒子分布函数φ的演化，代表了在时间t和离散速度vi处找到粒子在位置x的概率：

（公式1未在原文中给出，但根据上下文可以推断：φ(x, t, vi) = Γ(n(x, v(i), φ(y, t, vi-1))，其中Γ是递归函数。）

方程的左侧表示流动步骤，其中粒子分布沿着它们的速度方向在时间步长内传播到相邻的格点。右侧包含两个项：碰撞项和源项。碰撞项模拟分布函数φ向其平衡分布φeq的松弛，松弛速率由松弛时间τ决定。平衡分布φeq取决于局部组分密度ρi、平衡速度ui和离散格点速度vi：

（公式2未在原文中给出，但根据上下文可以推断：φeq = ρi * ui * exp(-τ * φ(x-xi)/N)）

该函数使用权重因子ωi和格点声速c来确保质量和动量守恒，以及各向同性。平衡速度ui等于质心速度v，表示所有组分的平均速度：

（公式3未在原文中给出，但根据上下文可以推断：v = (1/3) * (ρi * ui + ρj * vi) / (N * ρi + ρj)）

其中总密度N定义为所有格点上的粒子数。

表1介绍了离散格点速度和权重因子ωi，它们表示粒子在格子上可能的移动方向（图1）。源项表示相互作用力和外力对粒子分布的影响，如下所示：

（公式4未在原文中给出，但根据上下文可以推断：F = ∑i=1^N [ρi * ωi * φ(xi, t-μ) * f(xi, yj, ti)]，其中f是作用在粒子间的力。）

多组分多孔介质中时间步长t的总力F是流体–流体相互作用力Fff、固体–流体相互作用力Fsf和外力Fa的总和：

（公式5未在原文中给出，但根据上下文可以推断：F = Fff + Fsf + Fa）

表1中的ωi和f的值分别表示第i组分在图1中可能移动的方向及其相对重要性。这种力捕捉了多种流体相之间的非理想相互作用，使得能够对相分离和界面动力学等现象进行建模（Yang和Boek，2013年；Zudrop等人，2017年；Guo等人，2002年；Chen等人，1996年；Sukop和Thorne，2006年）。该力与固体-流体相互作用强度以及流体组分的赝势成正比，计算公式如（9）所示。它是通过对晶格方向进行加权求和得到的，其中固体指示函数用于识别相邻的固体晶格点。这一项允许在复杂几何形状中模拟润湿性和固体-流体相互作用（Shan和Doolen，1995年；Chen等人，1996年；Pelusi等人，2022年）。总力（公式（6）-（9）通过公式（5）中定义的源项进入公式（1），该源项考虑了扩散的流体-流体和固体-流体界面。通常使用标准反弹方法和中间位置反弹方法来管理固体-流体界面（Ginzbourg和Adler，1994年；Ladd，1994年）。当分布函数到达边界时，它会简单地弹回之前的节点，从而在边界处实现无滑移条件。在标准反弹方法中，固体边界与固体节点重合，反弹发生在时间步的末尾。在中间位置反弹方法中，边界位于流体节点和固体节点之间的一半位置，反弹发生在时间步的一半处。与标准反弹方法相比，这种方法提高了捕捉固体-流体界动态的准确性。周期性边界条件给人一种无限域的假象，即从域的一边出来的粒子群体会不受宏观速度或密度限制地重新进入另一边（Chen等人，1996年；Sukop和Thorne，2006年）。在本研究中，使用了具有周期性边界条件的中间位置反弹方案。2.2. 状态方程和压力计算多组分SC模型中的压力可以表示为理想贡献与相互作用贡献之和：（10）状态方程（EoS）包括一个与流体密度成正比的理想压力项，以及一个与相互作用强度成正比且依赖于赝势的相互作用压力项：（11）该公式同时捕捉了理想气体行为和多相系统中的非理想相互作用（Shan和Doolen，1995年；Krüger等人，2017年）。如公式（6）-（9）所定义的，组分之间的相互作用力负责压力方程中的非理想贡献。非润湿流体（）和润湿流体（）的固有平均相压力是通过对其各自体积内所有晶格点的压力进行平均计算得到的：（12）这里，和分别表示非润湿流体和润湿流体体积内的晶格点总数，而是从公式（11）得到的局部压力。在强相分离情况下，这与其主导组分的压力相等，因此平均值和可以解释为平均相压力。毛细压力（）定义为非润湿流体和润湿流体的平均相压力之差（Bottero等人，2011年；Whitaker，1977年）：（13）在本研究中，这个定义对应于体积平均量：（14）Pcvol=pΩnwbulk?pΩwbulk。此外，还考虑了两种局部度量方法。一种基于Young-Laplace关系的局部度量方法，通过为每个界面分配一个毛细压力Pc,iYL（根据其曲率），然后使用界面长度作为权重对所有界面进行平均：（15）Pc,iYL=γκi,PcYL=∑i=1NintLiPc,iYL∑i=1NintLi。第二种基于状态方程的局部度量方法，通过从界面附近的局部体积采样为每个流体-流体界面分配一个毛细压力Pc,ilocal，然后使用界面长度作为权重对所有界面进行平均：（16）Pc,ilocal=pΩnw,ibulk?pΩw,ibulk,Pclocal=∑i=1NintLiPc,ilocal∑i=1NintLi。这里，〈p〉Ω表示一组晶格点Ω上的平均局部压力，Ωnwbulk和Ωwbulk表示完整域中非润湿相和润湿相的体积区域，Ωnw,ibulk和Ωw,ibulk表示与界面i相关的相应局部体积区域，Li是界面长度，Nint是界面总数，γ是界面张力，κi是界面i的符号曲率。下载：下载高分辨率图像（254KB）下载：下载全尺寸图像图2. 二元流体系统的区域划分。晶格点被分类为临界区和体积区。在基于距离的标准中，与界面有间隙距离DI≤DcI的点被分配到临界区，而DI>DcI的点则构成体积区。在基于占据分数的标准中，χ≤χc的点被分配到临界区，而χ>χc的点构成体积区。NW：非润湿流体，W：润湿流体。2.3. 用于体积压力采样的自适应过滤过滤的目的是将晶格分类为临界区和体积区，如图2所示。临界区是扩散界面，包含具有虚假速度的晶格点。所有毛细压力估计都排除了临界区。体积区包括扩散界面之外的晶格，没有任何污染。比较了两种过滤技术：一种基于距离，另一种基于占据分数。在基于距离的过滤中，晶格点根据它们到相团边界的欧几里得距离进行分类。对于每个连接的润湿或非润湿团簇，如果一个点与其最近相边的距离超过规定的界面间隙距离DI，则将其保留为体积核心点，其中相边是指团簇的任何流体-流体或固体-流体边界。临界间隙距离阈值DcI决定了晶格点的分类：对于DI ≤ DcI，该点被分类为位于临界区内；而对于DI>DcI，则被分类为位于体积区内。在第二种过滤方法中，区域分类基于每个晶格点上主导组分的占据分数。在多组分LBM模拟中，每个晶格点包含两种共存的密度，属于密度较高的组分。晶格点的总密度ρtotal是两种密度之和。组分的占据分数定义为χσ=ρσ/ρtotal，代表其对总密度的贡献。主要占据分数χ=max(χw,χnw)表示两种组分分数中较大的一个。通过为主要占据分数χ设置特定阈值来划分晶格点到不同的区域。特定阈值χc决定了晶格点的分类。对于χ≤χc，该点被分类为位于临界区内；而对于χ>χc，则被分类为位于体积区内。2.4. 转换为物理单位LBM中的任何量（Q）都可以使用公式Qp=CQ从其晶格单位转换为其对应的物理量（Qp）。关键量的转换因子（C）值总结在表2中。物理单位下的晶格间距和时间步长分别为：（17）Δxp=ClΔx，Δtp=CtΔt。物理单位下的运动粘度νp与松弛时间τ和晶格声速cs,p的关系为（18）νp=cs,p2τ?12Δxp2Δtp，而物理和晶格表面张力之间的关系为（19）γp=γ×ρp×Δxp3Δtp2。表2. 转换因子。属性符号值长度ClΔxp时间CtΔtp密度Cρρp压力CpCρ(Cl)2/(Ct)2速度CuCl/Ct粘度Cν(Cl)2/Ct表面张力CγCρ(Cl)3/(Ct)23. 模拟设置在我们的研究中，我们探索了多种2-D系统，包括拉普拉斯测试、静止滴测试、板间的液桥以及部分饱和的有序和异质多孔介质。这些系统是使用3-D LBM晶格模拟的，在该晶格中，相同的2-D截面在垂直于平面的方向上均匀重复在几个晶格节点上。为了解决这种复杂系统的计算挑战，我们使用了LB3D软件，该软件以其在并行计算方面的效率而闻名（Schmieschek等人，2017年）。这些设置突出了在计算平均毛细压力时过滤临界区域的重要性。SC模拟与颜色梯度（CG）模型进行了基准测试，后者在这里被用作参考公式，因为它通过扰动和重新着色步骤保持了更清晰的流体-流体界面（Leclaire等人，2017年）。进一步使用Young-Laplace关系作为分析参考，以验证毛细压力与界面曲率的一致性。在继续进行多孔介质模拟之前，使用拉普拉斯和静止滴测试评估了SC公式的毛细压力网格收敛性。有序多孔介质的特点是两种不同孔径的颗粒以2的比率排列在重复的六角形晶格配置中。在孔内初始化液滴，导致相邻颗粒之间形成液桥，如图3(a)所示。下载：下载高分辨率图像（545KB）下载：下载全尺寸图像图3. 多孔介质系统。(a) 具有液桥的有序多孔介质，采用双分散六角形晶格配置。(b) 在异质多孔介质中的双层结构初始化。(c) 平衡相分离状态。(d) 固体-流体边界处中间位置反弹方案的示意图。箭头指示分布函数向边界流动并在中间点（半个晶格间距和一个时间步长）反射。PBC表示所有边的周期性边界条件。异质多孔介质由75个随机排列的圆形颗粒定义，并允许它们接触。系统通过交替的润湿和非润湿层进行初始化，如图3(b)所示。这种初始化确保了模拟的可重复性，这对于网格收敛性研究尤为重要。图3(c)展示了经过足够多的时间步后从双层结构初始化得到的相分离。图3(d)展示了固体-流体边界处的中间位置反弹方案。在每种情况下，分辨率由域维度Lx和Lz给出；二维平面上的相应网格点数为Lx×Lz。流体的润湿和非润湿性质分别由固体-流体相互作用强度gww、gwnw控制，这些强度代表接触角（公式（9））。流体的互溶性由参数Gσσ?控制，该参数与表面张力相关，此后称为G。我们使用了一个密度和粘度相等的双液模型系统，界面张力σ=70dyn/cm，代表水-空气界面。异质多孔介质的基线网格分辨率、模拟参数、物理缩放值、流体和颗粒属性在表3中进行了总结。图3(c)所示系统的孔隙几何结构以基线分辨率的1.2倍分辨率进行了处理，并用于第5.2节的Pc比较。从孔隙空间的中部骨架中提取边缘和顶点。在基于骨架的方法中，孔隙喉部位于距离场沿分支的局部最小值处，而孔隙体通常位于连接节点（Dong和Blunt，2009年；Lindquist等人，1996年；Lindquist等人，2000年；Prodanovi?等人，2006年；Jiang等人，2007年；Sheppard等人，2005年）。与这些研究一样，喉部定义在距离场沿骨架分支的局部最小值处，面积来自内切圆盘。然而，孔隙体是通过分割孔隙空间获得的（Vincent和Soille，1991年；Sheppard等人，2004年）。孔隙体的大小是从截面面积中计算得出的，排除了喉部的面积。图4显示了由此产生的喉部和孔隙体大小分布。在这种分辨率下，平均喉部宽度为13.5个晶格节点。表3. 基线网格分辨率、模拟参数、物理缩放因子、以晶格单位（lu）和相应物理单位报告的流体和颗粒属性。数量符号Lu值物理单位物理值基线网格分辨率计算域Lx×Lz750 × 750––平均喉部宽度w?t11.25––模拟参数流体-流体相互作用强度G5.6––松弛时间τ2.65––非润湿壁相互作用gwnw?0.06––润湿壁相互作用gww0.06––物理缩放因子长度尺度Δxp1m1.25×10?7时间尺度Δtp1s1.373×10?9质量尺度Δmp=ρpΔxp31kg1.95×10?18流体和颗粒属性参考密度ρ1.0kg/m31000运动粘度ν0.7167mm2/s4.1界面张力γ0.0375dyn/cm70非润湿接触角θnw–°118润湿接触角θw–°62最小颗粒直径dmin56μm基于状态方程精确计算毛细压力需要精确过滤临界区域和适当的分辨率。我们使用第2.3节讨论的两种自适应过滤方法，并在第3节引入的简单配置进行了比较。然后通过网格收敛性测试证明了这种方法的有效性。下载：下载高分辨率图像（222KB）下载：下载全尺寸图像图4. 异质介质中的孔隙喉部和孔隙体宽度分布。宽度通过最小颗粒直径dmin进行归一化。下载：下载高分辨率图像（237KB）下载：下载全尺寸图像图5. Young-Laplace曲率提取。SC和CG公式的流体-流体界面叠加图，SC的界面用圆形拟合（蓝色）表示。固体-流体界面用橙色表示。毛细压力测量值与过滤标准的比较
第2.2节中定义的毛细压力测量值Pcvol和PcYL针对简单配置进行了比较。对于每种设置，Pc均根据数值基准CG进行评估。图5显示，两种公式下的流体-流体配置结果非常接近，同时展示了SC的圆形拟合。提取的曲率用于估计两种公式的解析参考值PcYL。如图6(a)所示，两种公式的表面张力相似。图6(b)-6(d)展示了Pcvol与用于选择体核的界面间隙距离DI的依赖性。具体来说，对于每个规定的DI，体核由所有与界面间隙超过DI的位置组成（第2.3节）。在有序介质设置中，PcYL根据方程(15)计算为所有界面的界面长度加权平均值。在所有简单配置中，Pcvol(DI)达到平台值并趋近于PcYL。平台的开始点用于选择体采样的临界间隙距离阈值DcI和相应的Pcvol估计值。

下载：下载高分辨率图像（565KB）
下载：下载全尺寸图像

图6. 基于距离的过滤。从与所有界面间距大于DI的位置采样的相体压力获得Pcvol。DcI表示Pcvol达到平台值的起始点，此后Pcvol趋近于PcYL。SC公式作为数值基准与CG进行评估，而PcYL为相应公式提供解析参考。

图7显示了当使用占据分数准则选择体核时Pcvol估计值的变化。提高阈值χ可以仅保留体核中更纯由单一流体相占据的位置。得到的Pcvol值与基于距离的过滤得到的参考估计值进行比较。对于CG公式和没有固体的SC情况，占据分数准则得到的Pcvol估计值与PcYL一致。然而，对于含有固体的SC模拟，根据润湿实现的不同，基于χ的过滤可能无法恢复接近PcYL的Pcvol值，因为它可能在润湿相和非润湿相之间引入不对称过滤。因此，基于距离的准则在复杂多孔几何形状中更为稳健，并且在SC和CG公式中保持一致。

下载：下载高分辨率图像（592KB）
下载：下载全尺寸图像

图7. 基于占据分数的过滤。根据占据分数χ对体进行采样得到的毛细压力Pcvol。蓝色实线是基于距离的过滤得到的Pcvol估计值。

4.1. 网格分辨率的影响
网格收敛性研究对于确保LBM模拟中毛细压力的准确评估至关重要。提高模拟的分辨率会增加域中的晶格点数，从而更详细地表示物理系统，前提是晶格间距Δxp也相应减小。在SC模型中，流体属性由相互作用强度和松弛时间等参数决定，而不是明确定义的，因此流体属性会随分辨率变化。为了保持恒定的声速cs，需要根据方程(18)调整松弛参数。图8(a)显示了不同相互作用强度参数G下的Laplace测试中毛细压力与逆滴半径的线性拟合得到的表面张力值。图8(b)显示了两个代表性G值下的界面密度剖面，表明改变G会改变体平衡密度和界面宽度。图8(b)-8(d)显示了两个G值下平衡状态下非润湿相占据分数χnw的快照。为了保持不同分辨率下的表面张力和毛细效应，必须根据缩放关系（方程(19)相应调整G。

G与γ之间的关系由(20)给出：
γ=0.0310G?0.1352

下载：下载高分辨率图像（523KB）
下载：下载全尺寸图像

图8. 流体-流体相互作用强度G对表面张力和界面宽度的影响。(a) 不同G值下，Laplace测试中毛细压力作为逆滴半径1/r的函数；图例中的表面张力γ是从每个G的线性拟合斜率得到的。(b) 两个代表性相互作用强度G=4.8和G=6.4下的界面密度剖面。实线和虚线分别表示润湿相和非润湿相的密度。(c,d) 图(b)中两个G值下非润湿相占据分数χnw的快照。

对于网格细化因子λ=Δxcoarse/Δxfine，保持表面张力所需的缩放为：
(21) Gfine=λGcoarse?0.1352(λ?1)0.0310

因此，随着分辨率的提高，必须重新校准G和松弛参数以保持流体的物理属性。为了在Laplace测试中保持恒定的70达因/厘米的表面张力，根据方程(21)调整G。图9(a)显示了不同分辨率下毛细压力与逆半径1/r的关系。图9(b)显示了不同初始滴体积分数(Vi)下模拟中滴半径与分辨率的关系，表明随着分辨率的提高，滴尺寸趋于稳定。

对于静止滴测试，模拟了在非润湿流体中的润湿流体滴在板上静止的情况，接触角由固-流相互作用强度差gwnw?gww确定。图10显示了网格收敛性研究的结果。图10(a)展示了接触角随gwnw?gww的变化，而图10(b)展示了不同分辨率下的相应毛细压力。当G进行缩放后，超过某个分辨率时，接触角和毛细压力都不再随网格尺寸变化。

下载：下载高分辨率图像（424KB）
下载：下载全尺寸图像

图9. Laplace测试中的网格收敛性。(a) 通过随分辨率缩放G来观察不同网格分辨率(Lx×Lz)下的表面张力收敛性。(b) 随着网格分辨率的提高，不同初始体积分数(Vi)下的滴半径收敛性。

下载：下载高分辨率图像（449KB）
下载：下载全尺寸图像

图10. 静止滴测试中的网格收敛性。(a) 不同gwnw?gww值下不同网格分辨率(Lx×Lz)下的接触角收敛性。(b) 当G随分辨率缩放时毛细压力的收敛性。这是在不同分辨率下保持表面张力和接触角一致的结果。

图11显示了部分饱和有序多孔介质中的网格收敛性，展示了使用Laplace和静止滴测试中应用的G缩放方法，不同分辨率下毛细压力随gwnw?gww的变化。由于 pore 结构的均匀性（如孔径和形状），我们观察到毛细压力与接触角之间存在单调关系。相比之下，在异质多孔介质中实现网格收敛性更具挑战性，因为孔结构不规则。

下载：下载高分辨率图像（196KB）
下载：下载全尺寸图像

图11. 有序多孔介质中的网格收敛性。在不同分辨率(Lx×Lz)下，使用适当的G缩放方法得到的毛细压力与gwnw?gww的关系。

5. 异质多孔介质中毛细压力的计算
除了前一节讨论的因素外，异质多孔介质中毛细压力的准确性还取决于识别和解析影响它的特征。因此，详细检查了SC公式中网格分辨率的作用。通过将其与数值基准CG进行比较来评估体积平均测量值Pcvol（方程(14)的准确性。使用Young-Laplace关系验证了在单个界面处测量的局部毛细压力Pc,ilocal。然后，比较了两种公式下的局部毛细压力测量值Pclocal（方程(16)的准确性。最后，比较不同的毛细压力测量值，以评估单个毛细压力值在多大程度上能代表异质多孔介质。

下载：下载高分辨率图像（241KB）
下载：下载全尺寸图像

图12. 异质多孔介质中的网格收敛性。不同网格尺寸(Lx,Lz)下的毛细压力作为阈值函数的收敛性，通过物理缩放域来实现收敛。

5.1. 网格分辨率
确定异质多孔介质仿真的最佳分辨率可能具有挑战性。我们使用粗分辨率计算初始晶格间距估计值。在保持所有其他模拟参数不变的情况下，随着网格尺寸(Lx, Lz)的增加而缩放物理系统。图12显示了异质多孔介质中毛细压力随阈值的变化。观察到当Lx和Lz = 750时，毛细压力收敛，这被选为基准分辨率。

为了识别和解析对异质多孔介质中毛细压力至关重要的特征，我们对孔隙空间进行距离变换，测量每个点到最近固体边界的欧几里得距离。该距离场的局部最大值表示孔隙中心，每个孔径宽度估计为相应最大距离的两倍。孔径宽度有助于生成完整的孔隙结构特征，因为它们包括喉部、孔体和过渡宽度，后者是指喉部和孔体之间的中间宽度。图13显示了不同分辨率下归一化孔径宽度的分布，突出了临界孔径宽度与界面宽度之间的关系。对于我们的系统，临界孔径宽度被确定为0.19dmin，对应于毛细压力收敛时的孔径宽度超过界面宽度的分辨率。

下载：下载高分辨率图像（803KB）
下载：下载全尺寸图像

图13. 归一化孔径宽度的分布。直方图显示了随着分辨率的提高，归一化孔径宽度(pw/dmin)的分布，其中pw通过最小粒子直径进行归一化。y轴表示每个区间的孔隙百分比。第一个条形（灰色）表示亚临界孔隙。虚线黑色垂直线表示临界孔径宽度，实线红色表示全流体-流体界面宽度，两者均通过dmin进行归一化。在基准分辨率(d)下，毛细压力收敛，临界孔径宽度超过界面宽度。

在远低于临界宽度的尺度上解析孔隙狭窄导致毛细压力仅有2%的相对变化。分析表明，当临界孔径宽度相对于界面厚度得到充分解析时，控制全局体积毛细压力的流体-流体界面就得到了有效解析。对于特定的流体分布，控制狭窄的部分可以比平均孔径宽度更窄，并不需要是域中的全局最小值。

在确定基准分辨率后，重要的是在保持第4.1节中解释的相同流体属性和物理域的情况下进行网格收敛性测试。将静止滴测试中使用的G缩放应用于超出基准的分辨率。这种缩放保持了不同分辨率下一致的毛细压力平均值。然而，保持不同分辨率下一致的流体分布具有挑战性：研究表明，在异质样品中，每个分辨率下孔的形状和大小解析方式不同，即使在相同的指定接触角下，不同的孔几何结构也会表现出不同的行为（Mohanty等人，1981年；Rabbani等人，2018年；Scanziani等人，2017年）。分辨率的变化隐含地改变了孔的几何形状，导致有效接触角和流体分布的空间变化。附录说明了理想化孔结构和复杂孔结构之间的差异。

因此，对于异质介质，平均毛细压力等宏观参数是一个实用的收敛性标准，前提是在不同分辨率下尽可能一致地表示流体-流体配置，以最小化平均曲率半径的差异。在本研究中，随着分辨率从基准值提高，平均毛细压力的相对差异为2.4%。

下载：下载高分辨率图像（205KB）
下载：下载全尺寸图像

图14. 通过应用G缩放，在异质多孔介质中观察到的网格收敛性。

5.2. 比较毛细压力测量值
在本节中，我们分析和比较了异质多孔介质的体积平均毛细压力(Pcvol)和局部毛细压力(Pclocal)。图15(a)显示了使用基于距离的过滤（第2.3节）选择的体核评估得到的SC和CG公式的Pcvol。它展示了随着采样核心与固体-流体和流体-流体界面的间隙距离DI的增加，Pcvol的变化。图15(b)展示了用于平均的SC公式的润湿和非润湿体核，以及由平台开始时的间隙距离DcI定义的临界区域。图16直观地展示了异质多孔介质中同相簇之间的压力空间非均匀性。对于每个相，计算每个不连通簇的平均压力p?i。在图中，簇根据簇的压力差p?i?〈p?i〉进行着色，其中〈p?i〉表示该相所有簇的p?i的平均值。这些变化表明，同相的不连通簇可以具有不同的压力水平，这激发了接下来对基于体积平均值和基于界面的毛细压力测量方法的比较。下载：下载高分辨率图像（547KB）下载：下载全尺寸图像

图15. 异质多孔介质中不同公式之间的比较。(a) 通过基于距离的过滤选出的体积核心得到的Pcvol，与距离界面DI的距离绘制，以SC公式作为数值基准进行评估。(b) 由于基于距离的过滤过程而产生的域的分区。

局部毛细压力Pclocal是根据第2.2节中介绍的流体-流体界面平均值获得的。首先从每个界面周围构建的局部区域中的非润湿侧和润湿侧之间的压力差估算Pc,ilocal。然后从相场?=ρw?ρnw中提取平滑的?=0等高线，并以参数形式表示为x=x(s)和z=z(s)，其中s是弧长。接着评估局部曲率为(22)κ(s)=x′(s)z′′(s)?z′(s)x′′(s)x′(s)2+z′(s)2/3，这里κ(s)是沿界面的点态曲率，而κrep作为κ(s)的平均值，是用来定义该界面的单一有效1/R的代表性曲率。在获得各个局部压力后，可以使用公式(16)计算Pclocal。下载：下载高分辨率图像（779KB）下载：下载全尺寸图像

图16. 异质多孔介质中同一相内簇压力的空间非均匀性。(a) 润湿相和(b) 非润湿相簇根据它们与相应相的平均压力的差值进行着色。固体显示为橙色。

这种方法用于计算单个界面的Pc,ilocal在图17中进行了检验，使用了局部Young-Laplace关系。图17分别将Pc,ilocal与SC和CG公式的曲率1/R进行绘制，在这两种情况下，拟合的斜率都很好地符合Young-Laplace定律规定的0.0375。

接下来，图18比较了整个域和各个簇对的基于体积和基于界面的毛细压力测量结果。图18(a)使用箱形图总结了Pc,ilocal的异质分布，并将其与整个域的Pcvol进行比较。在两种公式中，Pcvol都高于局部毛细压力分布的中位数。图18(b)将成对的体积平均毛细压力Pcvol,pair与长度加权的局部毛细压力Pclocal,pair进行比较，适用于所有唯一的润湿和非润湿（W–NW）簇对。因此，每个点代表一个独特的W–NW簇对。Pcvol,pair是从相应对内的非润湿相和润湿相之间的压力差中获得的，而Pclocal,pair仅使用该对共享的界面计算得出，每个界面都按其界面长度进行加权。这两种测量结果对所有点都是一致的，表明对于个别簇对而言，总体压力差异与相应的局部界面毛细压力是一致的。下载：下载高分辨率图像（355KB）下载：下载全尺寸图像

图17. 异质多孔介质中局部毛细压力的Young-Laplace验证。(a) 和 (b) 将每个界面的局部毛细压力Pc,ilocal与提取的1/R以及数值基准CG进行比较。图19(a)显示了随着相邻簇逐渐加入主导的W–NW对（定义为具有最大共享界面长度的对）时，总体毛细压力测量的演变。由此产生的流体配置形成了越来越大的异质簇集，因此分析逐步向总体平均值过渡。Pcvol始终高于Pclocal，并且从早期阶段到生长最终阶段的偏差越来越大。图19(b)对于选定的几何形状，比较了在此生长过程中各个簇对总体相体积和总体W–NW界面长度的相对贡献。对于每个相，wV表示给定簇贡献的总体相体积的比例，wL表示该簇集内总体W–NW界面长度中属于该簇的比例。在每个生长阶段，箱形图是根据润湿和非润湿簇的wL/wV值组合构建的。从早期阶段到最终阶段，wL/wV分布的扩宽表明，随着簇的逐渐增加，簇对界面长度和相体积的贡献比例越来越不同。下载：下载高分辨率图像（350KB）下载：下载全尺寸图像

图18. 基于界面和基于体积的毛细压力比较。(a) 显示了整个域的Pc,ilocal分布。标记表示Pcvol，箱形线表示Pc,ilocal的中位数。(b) 对于SC，比较了成对的体积平均值Pcvol,pair和仅使用属于该簇对的界面计算的基于长度加权的局部毛细压力Pclocal,pair。下载：下载高分辨率图像（272KB）下载：下载全尺寸图像

图19. SC生长过程中总体毛细压力和簇分布的演变。(a) 在簇聚集生长的三个不同阶段，体积平均毛细压力与局部总体毛细压力之间的差异。(b) 在簇聚集生长的三个不同阶段，以簇界面长度分数与簇体积分数比的形式显示簇的异质性。

从图20(a)中对SC公式所示的Pc,ilocal的长度加权分布可以理解Pcvol和Pclocal的相对位置。观察到大部分总体界面长度与Pc,ilocal≤Pcvol的界面相关联。图20(b)进一步显示，Pc,ilocal>Pcvol的界面通常比Pc,ilocal≤Pcvol的界面更短。这些结果解释了为什么在本例中，Pcvol位于局部毛细压力的长度加权平均值之上。在这种情况下，差异大约为10%。

我们回顾了简单的有序多孔介质，在这种介质中，Pcvol收敛到界面长度加权的局部测量值PcYL（第4节）。这是均匀局部毛细压力和簇贡献的极限情况。然而，当同时满足两个条件时，体积平均值可能与典型界面毛细压力分布的值不同：(i) 局部毛细压力是异质的；(ii) 簇的分布是异质的，簇对选定集合的总体体积和总体界面长度的贡献不成比例。下载：下载高分辨率图像（402KB）下载：下载全尺寸图像

图20. SC公式的Pc,ilocal和界面长度分布。(a) 以Pcvol为参考的Pc,ilocal的长度加权分布。(b) 对于Pc,ilocal≤Pcvol和Pc,ilocal>Pcvol的界面，其界面长度分布。

6. 讨论

需要注意的是，本研究仅限于准静态情况，这常见于部分饱和多孔介质中的自发吸收或缓慢排水过程。在这些条件下，毛细力占主导地位，因此分离它们的效应以理解毛细压力-饱和关系非常重要。所提出的过滤方法结合网格分辨率要求和局部毛细压力分布，有助于改进使用扩散界面模拟获得的毛细压力-饱和关系的估计和解释。在土壤科学中，毛细压力-饱和关系用于灌溉的水分利用效率和调度（Clothier和Green，1994年，Obreza等人，1997年），而在岩土工程中，它用于通过将土壤毛细压力与剪切强度和安全系数联系起来进行边坡稳定性分析（Fredlund等人，1978年，Ng和Pang，2000年），用于非饱和土壤的承载能力评估（Vanapalli和Fmo，2007年），以及用于沉降预测（Jeong等人，2018年，Onyelowe等人，2022年）。在水库工程中，这种关系有助于水库特性分析、评估密封能力、识别含水层和估算碳氢化合物柱高度（Vavra等人，1992年）。最后，在材料设计中，它可以用来表征土工合成材料作为堤坝中的毛细屏障或吸水层时的非饱和流体响应（Zornberg等人，2010年，Jarjour等人，2024年）。

使用基于距离的标准进行体积核心识别的方法在准静态条件下是稳健的，但在动态场景中，对于宽度与扩散界面宽度相当的喉部，体积核心的检测变得敏感，需要比这里考虑的准静态情况更严格的网格分辨率。此外，相位平均压力差包括了粘性和惯性效应以及毛细压力，因此不能纯粹解释为毛细效应（Hassanizadeh和Gray，1993年，Sánchez-Vargas等人，2023年，Lasseux和Valdés-Parada，2023年）。流体的大密度和粘度对比会进一步放大虚假速度，而且在大粘度对比的情况下，达到真正的准静态平衡可能非常慢。如果模拟在较小的但仍有限的毛细数下停止，测量的毛细压力仍然会包含粘性效应。

当前的过滤方法依赖于在逐步过滤受污染的界面单元后检测平均毛细压力的平台，这只有在临界孔径得到良好解析时才有意义。确定临界孔径的程序是通用的，可以应用于任何多孔介质。然而，所得到的分辨率标准必须针对每种几何形状重新校准。对于自然多孔介质中看到的高度异质性几何形状，分辨率要求会过度解析大孔，浪费计算资源，使得通过将分辨率集中在狭窄的喉部来采用自适应网格细化成为可行的解决方案。

对于动态场景或大的流体对比，建议结合使用改进的方案，如MRT或累积（中心矩）碰撞算子（d’Humières，2002年，Geier等人，2015年），更精细的状态方程（如Peng–Robinson或Carnahan–Starling（Yuan和Schaefer，2006年），更高阶的各向同性或多范围相互作用（Sbragaglia等人，2007年，Qin等人，2023年），以及颜色梯度LBM（Leclaire等人，2017年，Joseph等人，2026年），这些都是有助于减少虚假电流和提高稳定性的方法。将这种方法扩展到代表自然岩石和土壤的完全三维几何形状，具有实际的密度和粘度对比（例如，水-空气），将是未来工作的重点。在这样的介质中，颗粒角度可能导致更不规则的孔几何形状，而更紧密的堆积会促进相捕获。这些效应可能会增加簇和压力的异质性，使得体积平均值和局部界面毛细压力之间的区别更加重要。

7. 结论

在这项研究中，使用Shan–Chen LBM来研究部分饱和多孔介质。主要目标是在扩散界面方法的框架内改进毛细压力（Pc）的估计。讨论了一种在复杂多孔材料模拟中改进这些估计的实际方法，以及网格收敛性和基线分辨率的指导原则。使用越来越复杂的二维设置，将体积平均Pcvol与基于界面长度加权的局部测量值进行比较，后者可以通过直接提取的Young–Laplace界面曲率获得（简单测试和有序多孔介质），或者在更复杂的情况下通过界面Pclocal的局部采样获得（异质系统）。

扩散界面和体积之间的边界位置对于计算Pc至关重要，该边界定义了临界区域。实际上，使用Shan–Chen公式获得的毛细压力与Young–Laplace分析参考值和高度准确的LBM颜色梯度方法的吻合度很好，只有在精确识别体积区域时才能实现。在简单测试（Laplace、静止液滴和板间的液桥）和有序多孔介质中，通过将体积平均压力Pcvol与长度加权的Young–Laplace测量值PcYL进行比较来评估其准确性。同样，对于异质多孔介质，Pcvol随时间、距离界面的距离或网格大小的收敛性本身并不能确保物理准确性；它还应该得到局部尺度上与Young–Laplace关系一致性的支持。

我们的网格收敛性研究表明，在异质多孔介质中，Pcvol受到临界孔径的控制，该临界孔径被定义为必须超过扩散界面宽度的孔径，以便总体体积平均毛细压力能够收敛。解决这个宽度可以确保捕获到最具有影响力的限制界面并控制Pcvol的入口，其大小可以显著小于域内的平均孔径。在这里考虑的异质系统中，临界孔径被确定为最小粒子直径的0.19倍。因此，为了获得可靠的毛细压力估计值，需要足够高的基线分辨率，以确保关键的孔径宽度大于界面宽度。相比之下，在有序多孔介质中，确定控制孔径宽度相对容易，因为只有少数几种孔喉宽度存在，并且相邻颗粒之间分离得很好。此外，在细化网格以保持表面张力时，对流体-流体相互作用参数G进行适当的缩放也是必不可少的。在比较有序多孔系统中的毛细压力测量结果时，Pcvol与局部毛细压力的长度加权平均值非常接近。由于在这种情况下界面和簇在几何形状上相同，这两种方法得出的有效毛细压力几乎相同。然而，在具有复杂簇分布的非均质多孔介质中，毛细压力通常不能仅用一个标量来完全表征，而需要用局部毛细压力的分布来表示。在本研究的配置中，Pcvol与局部毛细压力分布的长度加权平均值存在差异。这种差异是由于局部毛细压力的非均匀性以及簇对总相体积和内部界面长度的不同贡献所共同作用的结果。因此，体积平均毛细压力应被视为描述整体流体配置的宏观尺度指标，而不是界面毛细状态的直接表示。

贡献声明：
S. Joseph：撰写——原始草稿、可视化、软件开发、方法学研究、数据分析、概念化。
H. Cheng：撰写——审阅与编辑、监督、方法学研究、概念化。
J. Harting：撰写——审阅与编辑、方法学研究。
V. Magnanimo：撰写——审阅与编辑、项目管理工作、方法学研究、概念化。