计算多孔介质的孔隙率和渗透率

在模拟多孔介质中的流动时，使用均质化的宏观方法可以有效地简化真实多孔材料的几何复杂性。但是，如果我们不知道有效的宏观特性是什么，该怎么做？本文我们将研究如何从一个完全解析的微观子模型中提取孔隙率和渗透率宏观流动特性。

在之前的文章中，我们讨论了可用于模拟多孔介质中宏观流动的接口，包括达西定律接口。求解达西定律可以深入洞察许多不同的物理系统，在这些系统中，模拟完全解析的微尺度系统几乎是不可能的。在石油和天然气、土木工程以及生物和医学工程等应用领域中，由于微观和宏观系统的长度尺度之间存在巨大的差异，因此带来了求解的困难。

宏观方法假设孔隙空间的行为由两个平均量来量化：

1. 渗透率

2. 孔隙率

渗透率是指描述流体流过孔隙的阻力。孔隙度的定义为孔隙空间的体积分数，决定了表观平均流速。至于表观流速，是指通过均质域的等效速度，也可以理解为：通过孔隙空间的微观流动在宏观尺度上是均匀分布的。

如果孔隙率和渗透率未知，就必须通过实验结果将这些材料特性量化。通过模拟进行的数值实验也可以分析完全解析的包括空隙和固体颗粒的几何结构。通过在微尺度几何上求解纳维-斯托克斯方程（或其对低雷诺数的线性近似，称为斯托克斯流或蠕动流），可以提取多孔介质的孔隙率和渗透率。

在研究多孔介质的孔隙率和渗透率之前，我们必须讨论微观几何结构的创建。这不一定是一个简单的过程！创建这类几何结构通常需要使用专门的第三方软件（如 Simpleware 或 Mimics®）重建扫描图像数据，尤其是复杂的 3D 几何图形。本文，我们将重点研究 2D 横截面几何结构，例如由扫描电子显微镜(SEM)图像生成的横截面。

• 如何使用拓扑优化结果创建几何模型

• 在 COMSOL 中处理不规则形状的方法

我们来看一个求解了完全解析的孔隙空间几何的孔隙尺度流动示例模型。我们可以使用 COMSOL Multiphysics 中的后处理工具来提取孔隙率并使用达西定律计算渗透率。我们来快速浏览一下这个案例。

由整个几何形状的尺寸 640μm×320μm（孔隙空间中的通道宽度更窄）可知，特征长度尺度 L 很小。此外，考虑到流动是由 2Pa 压力梯度驱动的缓慢流动，典型的速度值U应该很低。因此，给定类水流体的密度 ρ=1000kg/m^3，动力黏度 μ=0.001Pa*s，我们可以忽略 Navier-Stokes 方程的惯性项并使用蠕动流接口求解，因为雷诺数 Re 明显小于 1。

如果想获取更多信息，请阅读文章：如何正确模拟不同类型的流体流动。

使用入口和出口条件从右到左施加压降。由于该模型是我们希望描述的多孔介质的代表性横截面，因此沿边界指定了一个由顶部和底部的几何形状截断引起的对称条件。左图显示了这些边界，右图为使用包含速度大小的颜色图和流动方向的箭头矢量图显示的结果。

使用蠕动流接口求解的完全解析的微尺度模型的几何形状和边界条件。

首先计算孔隙率 por。在这个 2D 示例中，我们需要计算总面积和由计算域表示的面积。

我们可以简单地把总面积作为一个变量：A_tot = L*H 来计算。为了确定孔隙空间的面积 A_por，我们可以在流域上对表达式“(1)”进行积分。这可以通过积分组件耦合的方法轻松实现，这是一个可以在域、边界等处对任何表达式进行积分的自定义算子。下图中显示了用于计算 A_tot、A_por 和 por 的设置。

上图还显示了渗透率 k0 是如何计算的。达西定律指出：

其中，u 是达西或表观流速，κ 是渗透率，μ 是动力黏度，∇p 是压力梯度。

我们可以在蠕变流接口借助预定义的变量计算达西速度。

变量 spf.out1.Mflow 定义了通过出口边界的质量流量，我们可以用它除以恒定密度 ρ0 来获得体积流量。然后除以多孔介质的高度，再除以 1m 来考虑 2D 近似，获得 x 方向的达西速度。

通过重新排列表达式(1)并替换 ，我们可以使用已知的压降 p0，通过多孔区域的长度 L 来计算变量 k0 的渗透率。

结果表明，这种微观尺度代表多孔介质的孔隙率为 0.553 和渗透率为 4.59×10^-12m^2。

在这篇文章中，我们讨论了如何使用仿真推导通过多孔介质的流动的宏观特性——在完全解析的微观尺度子模型中求解自由流动。有了这些信息，我们就可以将这些参数作为更具描述性的宏观模型的输入。更重要的是，这是了解在一个用户界面友好的仿真 App 中输入什么的理想方式，例如这个考虑射孔井生产率和安全性的仿真 App 示例。