孔隙网络中砂岩和流体模拟的三维孔隙尺度建模

1、孔隙网络中砂岩和流体模拟的三维孔隙尺度建模Stig Bakke and Pal-Ericren王少依 翻译；杨晓宁 校对摘要：本文提出了一种新方法，用来建立真实均质和非均质三维孔隙尺度的砂岩模型。此方法的本质是通过模拟砂岩形成的主要地质学过程砂粒的沉积、压实及成岩作用建立实际砂岩体的类比模型。模型所需的数据来源于实际砂岩薄片的图像分析。砂岩模型在X, Y和Z方向上的空间连续性是通过基于独立比例侵入渗透算法确定的。最终得到的空间连续性椭圆函数，可以用来对砂岩模型进行非均质性描述。非均质性分析表明，与垂向相比，压实作用在水平方向上能够更迅速的减少空间连续性。孔隙空间所呈现的网络结构和几何形态，是直

2、接在经过充分表征的砂岩模型中应用各种3-D图像分析算法得到的。从一个强水润湿的Bentheimer 砂岩薄片数据中得到的3-D孔隙网络被用作一个双相网络流动模拟系统的输入数据。模拟得到的砂岩模型的传导特性和实验得到的结果具有良好的一致性。前言高渗透率空隙或孔隙的存在是储集岩区别于其他岩石的特征。储集岩的孔隙高度不规则，其包含的孔隙空间网络由较大孔隙（孔隙本身）和连通这些孔隙的较窄孔隙（孔隙喉道）组成。孔隙空间及其相应颗粒基质的结构和几何形态决定了岩石的一些宏观特性，如绝对渗透率、相对渗透率、毛管压力、地层因素以及电阻率系数等。由微观结构和孔隙尺度的物理现象预测岩石的宏观特性问题，一直是被广泛研

3、究的主题。此项研究中得到普遍应用的一个工具是网络流体模拟系统（网络模型）。该方法要求对发生在孔隙尺度的物理过程的细致理解，以及对孔隙空间形态的完整描述。通过在微模型中识别出的物理现象，此步骤已经在简单或理想多孔介质的两相流体中得到成功应用。最近，许多三相流体孔隙尺度的物理现象的公开，导致了三相流体网路模型的发展。由于对真实多孔介质三维孔隙结构的精确描述仍然比较困难，所以针对真实多孔介质的技术扩展也十分复杂。关于储集岩孔隙结构的信息经常借助于压汞数据和薄片图像分析。Yanuka 和 Wardlaw 对表征孔隙结构的方法进行了很好的概述。压汞数据提供了孔隙喉道大小分布的静态信息，或者精确点儿说，是

4、被标定的孔隙喉道大小中可能被侵入的体积的分布。对薄片图像进行的传统的孔隙几何形态分析，一般基于Delesse 定律，其认为薄片中的2-D观察可以代表3-D样本。但是，这只适用于标量，例如面积和体积方面的考虑。因此，这个定律适用于孔隙度和矿物成分的测定，而不适用于孔隙和喉道的几何形态测定，因为这些是向量。为了在薄片图像中表征孔隙或喉道的几何形态，有必要定义2-D孔隙元素中哪些是孔隙，哪些是喉道。进一步，我们需要知道每一个孔隙元中流体的流向，才能测量垂直于流向的平面上孔隙喉道的几何形态。当前，没有方法可以直接的从薄片中提取这些信息。因此，传统孔隙结构分析得到的数据，不足以直接提供对真实岩石复杂孔隙

5、网络的完全3-D描述。高级技术例如微CT, 连续切片结合BSE切片成像，为岩石的3-D孔隙结构提供了详细的描述。然而这些技术是昂贵而又低效的。利用Finney 构建的完全定义的等球粒密封体，Bryant 经过模拟砂岩形成过程的结果，建造了3-D孔隙网络。重建的多孔介质的方法也已经被用来表征3-D孔隙结构。多孔介质建模也成为了许多化学工程领域的研究课题。本文提出了一种新的建立真实而又充分表征3-D孔隙网络模型的方法。此方法基于砂岩形成的主要的地质学过程结果的数字模型。模型砂岩的孔隙空间网络的表示，通过图像分析技术，从相应的矿物基质网络中分离出来。流体流模拟直接在不规则的孔隙空间进行，然后计算传导

6、特性和现有的实验数据进行对比。砂岩的3-D过程/随机建模砂岩是一系列复杂的地质学和流体动力学过程的最终产物，其过程可以简述如下：开始是含石英岩石被侵蚀，然后是砂粒通过空气、水和冰川的搬运，最后在沉积盆地沉积。沉积的颗粒可能会经过再次或多次侵蚀，搬运和沉积的循环而被重新改造。经过压实作用和多种成岩作用的过程，最终形成了砂岩。砂岩样本及其岩相学参数是其经历的所有地质学和流体动力学过程导致的结果，这些过程影响了沉积盆地及其环境。我们并不想建立这些极度复杂过程的详细动力学模型，而是试图通过使用薄片分析得到的岩相学数据，模拟砂岩形成的主要过程的结果，从而重建3-D砂岩。砂岩重建过程包含三个主要的建模阶段

7、：（）沉积阶段（砂粒沉积），（）压实阶段，（）成岩阶段。（）沉积阶段沉积过程的建模开始于粒度分布曲线的产生，而粒度分布曲线来源于真实砂岩薄片的BSE图像。我们利用图像分析技术测定图像中所有砂粒的直径。目前，所有的碎屑砂粒被当作球形的石英颗粒。复杂的地质学和流体动力学过程影响了砂床上每一个单独砂粒的沉积，为了模拟其随机性，从所给粒度分布中选取单独颗粒也是随机的。随机颗粒的选择方式取决于粒度分布曲线的形状。如果曲线接近于一个已知的数学函数，粒度则按照已有的分布方式进行随机分配。否则，所有从薄片分析中记录的粒度都要按照大小编号和排序。然后，在0和已记录的颗粒总数之间挑选一个随机实数，对应的颗粒大小通

8、过在离散曲线点之间插值得到（图1）。粒度分布没有大小限制，从粉砂岩到砾岩都可以进行建模。图1 举例说明如何从一个已经编号并排序的粒度列表中随机选取粒度值动力学和几何学控制每一个砂粒在已有的底床上沉积，用数字来描述它们是相当复杂而又枯燥的。为了简化这个问题，每一个在沉积之前随机选择的新颗粒的半径折合为0。底床上所有颗粒的半径增加一个相当于新颗粒半径之和的值。这样每一个新颗粒缩小为一个点，底床则变为一个固体面（图2），从而在几何学和动力学上相对容易进行描述。新的“颗粒点”沉积到一个稳定的位置后，在对下一个颗粒重复以上步骤之前，所有颗粒的半径都要恢复其初始值。这个步骤持续进行直到一个预先定义好的次级

9、盆地单元（极限盒子）被颗粒充填。所有颗粒的中心坐标和半径都要被记录并且用球形展现（图3图7）。图2 沉积模型（低能沉积）2-D示意图。一个新的半径为R的砂粒沉积到底床表面。新颗粒简化为一个点，底床所有颗粒的半径增加R。新颗粒沉积到一个局部稳定的位置后，所有颗粒的半径恢复其初始值。图3 低能（局部极小值）沉积形成的开放底床，孔隙度=37.8%每一个新颗粒在底床上沉积的精确位置取决于其沉积过程中环境能量的强弱。对于低能环境，假设新颗粒的水平速度为0，且没有横向力影响底床。这种情况下，每一个随机选择的新砂粒（简化为一点）以一个随机的X,Y位置掉落在（扩张的）底床上。撞击到底床上以后，颗粒从一个坡度局

10、部极大值处滚动直到一个局部极小值的相对稳定位置（图2和3）。由此形成了具有局部峰谷的波状起伏的底床表面。然而，在极度低能的盆地中沉积的砂粒是罕见的。砂岩一般是高能沉积事件的产物（低能沉积形成泥岩和粉砂岩）。当水和空气等搬运介质的水平速度降低到某个临界值以下时，砂粒从中沉淀出来。这一临界值与粒度成比例。砂床经常受河流和波浪的水平作用力以及斜坡的重力作用影响。这些水平力搬运每一个颗粒，直到它们沉降到一个受搬运力影响最小的稳定位置。高能沉积建模即把每一个随机选择的砂粒放到底床上（图4）现有的最低位置（全局极小值）。与低能沉积相反，高能沉积下底床表面几乎是平坦的生长。图3和4中的底床孔隙度值表明高能（

11、全局极小值）沉积形成的底床比低能（局部极小值）沉积密度更大。通常，细粉砂沉积用低能沉积建模，而砂砾岩沉积用高能沉积建模。图4 与图3完全相同的输入参数的高能沉积（全局极小值）形成的密度更大的底床，孔隙度=35.3%图5 层理底床模型上面的模型展现的是均一的或块状砂岩模型。自然界的砂岩通常表现出不同规模一定程度的沉积非均质性。在m和mm级，孔隙空间本身都是非均质的。在mm到cm级，由于层理，递变层理和成岩作用过程，砂岩也可能是非均质的。层理砂岩建模即三种具有不同而又相互独立的粒度分布的纹层相互叠置而成（图5）。递变层理建模即向上变细或变粗的砂床，反应的是粒度随深度的渐变（图6）。图6 向上变粗的

12、底床模型图7 压实因子分别是0.0，0.1，0.2和0.3的线性压实模型产生的效应。（）压实阶段压实，或总体积减小，对上覆沉积垂直压力的响应是经历埋藏作用的石英砂岩减孔的重要因素。浅部压实是诸如颗粒重排、破碎等机械作用形成的。对于深部压实，富含K2O矿物的存在促进了压实作用，导致颗粒穿透是其主要的机制。即使在极度均一的砂岩中，孔隙空间的局部变化或非均质性也有可能是压实造成的。在连续切片模型，甚至在薄片中，有时可以观察到稳定的石英颗粒细胞结构是如何阻止压实力，并在其内部保存开放的颗粒结构；而在细胞外的颗粒结构有可能遭受深度的压实。压实作用的结果是一个类似于Bryant 描述的线性模型。每一个砂粒

13、中心的Z坐标与其初始Z坐标呈一个总比例的方式垂直向下移动。Z=Z0(1-)(1)其中Z是新的Z坐标，Z0是原始坐标，是一个范围从0到1之间的压实因子。图7和表1表示的是不同的压实因子对孔隙度的压实作用影响。（）成岩阶段成岩作用过程有时引入了非均质性概念，这些叠加在了砂岩原生构造的非均质性模式上。由于地质条件的改变，成岩作用过程可能十分复杂，涉及到不同时期的矿物溶解作用和次生加大胶结等，并且和压实作用相互影响。目前，仅已知的成岩作用中一些简化的模式得到模拟：石英的次生加大胶结和后来在自由表面生长的粘土膜。通过均匀地增加砂床上所有砂粒的半径，模拟石英的次生加大胶结（图8-10）。石英次生胶结总量可

14、能等于压实作用中颗粒互结生长损失的体积（遵守物质守恒定律），或者独立地通过胶结厚度输入参数来限定。后者模拟了SiO2运移进入或离开砂体模型。粘土膜在石英颗粒的自由表面、石英胶结物或已存在的粘土的表面随机沉淀下来（图8-11）。粘土沉淀的总量从岩相分析中得到。薄片中有时可以观察到，一些孔隙单元中粘土矿物相对富集，而另一些孔隙单元中则几乎没有粘土矿物。这种粘土膜现象可以通过聚合生长规律来解释，即新的粘土质点（3-D阵元或网格单元）沉淀到已有的性质相似的粘土质点上的可能性较大。现在，我们不区分不同类型的粘土矿物。模型假设粘土主要是由高岭石和一些伊利石组成的混合物。图8 由图4的底床模型建立的砂岩模型

15、的代表性等值面。石英颗粒为灰色，石英胶结物为浅灰色，粘土矿物为深灰色。边界条件；沉积和压实过程的模拟是基于目标的建模，随后是模型的3-D网格化。成岩作用模拟是在这个基于网格化的模型上进行。例如，一个100100110三维象元的3-D阵列定义了砂床中颗粒中心的空间限定或边界盒子。压实作用模拟时，边界盒子Z轴方向上的扩大取决于额外的空间要求。在自然界，这样一个接近于立方体的颗粒堆积是不稳定的，因为垂向角度超出了颗粒堆积的（最大）静止角。但这不是问题，因为我们只是模拟连续砂岩体的一小部分。表征砂岩模型的时候，为了减弱可能存在的边界效应，只取砂岩模型中心的部分做进一步的处理和分析。例如，如果工作阵列包

16、含100100110的三维象元，中间的707070象元才是我们感兴趣的区域。模型验证；沉积建模是模型处理中最耗费CPU的部分。对于一个具有285285315网格、平均粒度大小为30网格单位的模型，运用基于目标的平行计算软件在一个具4 个CPU的SGI Challenge R10000型计算机上一般需要运行3个小时。模型处理中的压实和成岩部分运用单循环计算，对于以上模型一般约需要5分钟。通过对比模型和薄片或（岩芯）柱塞的孔隙度，以及对比通过模型的可视化2-D任意截面和砂岩样品的BSE图像，可以验证砂岩模型（图12）。调整压实和成岩作用参数，直到2-D截面和薄片之间的孔隙度值和视觉误差均达到最小。

17、图9 图8中所示模型的代表性“椅子模式”等值面。石英为灰色，石英胶结物为浅灰色，粘土矿物为深灰色。模型可视化科学的可视化对于检查复杂的3-D数据和模型是一个非常有用的工具。目前，应用红外探测器对建立的砂岩模型及其孔隙空间进行可视化操作。图8展示的是由图4所示的模拟底床所建立的砂岩模型的代表性等值面。不同的灰度阴影代表了不同的矿物和胶结物。利用现在的可视化工具，可以把模型的一部分移除，研究其内部特征（椅子-模式）。例如图9和图10所示。运用二元砂岩模型及相互间的转换，孔隙空间和砂岩模型可以同时实现可视化（图10）。图11展示的是砂岩模型完整的孔隙空间。无论是3-D模型的正交或者随机的2-D截面（

18、模拟薄片）都可以提取并可视化（图12）。从而将砂岩模型中的模拟薄片和假设模拟的真实砂岩的薄片进行对比。这样才能确保模拟薄片和实际的岩石薄片具有相同的统计属性（包括岩相学、矿物学和几何学属性）。图10 相应的孔隙空间。孔隙空间为中灰色，石英灰色，石英胶结物浅灰色，粘土矿物深灰色。图11 图810中所示的模型中孔隙空间的代表性等值面图12 通过砂岩模型随机的任意方向的2-D截面表征空间连续性；孔隙空间的连续性在X、Y和Z方向上强烈地影响岩石中各种尺度的流体流动。目前，沿着指定流体流动方向的主要路径的最窄通道用来确定空间连续性。为了定量确定空间连续性，有必要对孔隙空间进行几何描述。侵蚀几何学分析；侵

19、蚀是一种众所周知的2-D图像分析功能，通过在双重图形中运用不同的像素遮罩和邻区连通性缩小目标的尺寸。通过3-D邻区连通性的概念，侵蚀功能已经延伸到3-D分析。图13是3-D侵蚀技术展示的孔隙本身和喉道的大小以及孔隙网络空间连续性的可视化信息。通过计算那些溶蚀特定象元所必需的溶蚀循环，可以获得定量的孔隙大小信息。图13 Bentheimer砂岩模型中5个侵蚀循环的最终结果。如右下角的图形所示，最大的孔隙半径大概为5.5 单位（箭头所示）。经过2次侵蚀，模型孔隙空间的连续性为0。最终的3-D量化侵蚀靠侵蚀几何学分析得以实现。图14表示的是侵蚀几何学分析的代表性2-D结果。这一侵蚀技术提供了孔隙本身

20、和喉道大小的量化信息，但是没有关于（孔隙）网络结构和通过孔隙空间的主要流体通道的直接信息。侵入渗透模拟通过在3-D孔隙空间模型的侵蚀几何学分析成果上输出的模拟侵入渗透类型过程，得到指定方向上通过孔隙空间的主要流体路径。侵入“流体”从一个侧面注入，按照侵入渗透规律渗透通过模型，也即在每一步骤，拥有最多与前缘接近的三维象元的孔隙空间被侵入。图15 展示了一个模拟侵入的过程。图14 经过理想孔隙结构2-D截面的二重图像的2-D量化侵蚀。0代表基质像素，1代表孔隙像素。右边的图像表示的是侵蚀每一个像素所需的侵蚀环数目。最大的侵蚀环数目是9。也就是说左边图像最大孔隙体积的半径是9单位。孔隙喉道的半径是2

21、单位。当侵入“流体”到达对面的时候（也即突破通过），侵入体积表明了指定方向上最适合的流体路径。 “流体”侵入的最窄通道（最低象元值，最小孔）可识别出。此通道决定了排替压力的大小，而排替压力可以确定通过指定方向孔隙空间的连续流体。利用最窄通道的大小可以测定指定流体方向上的空间连续性。其他两个主要方向上的孔隙空间的连续性可以通过同样的方法测定。最终的空间连续性函数可能是一个椭圆函数，并且可以用来描述孔隙空间的非均质性（表1）。表1 压实作用对孔隙度及其各向异性的影响表1 颗粒堆积的孔隙度随压实因子的增大呈近似线性下降，且X，Y方向的空间连续性比Z方向减小的更快。除去那些受有限的分离点值影响的小波动

22、，非均质性（各向异性）随压实的增强而增强。这和Bryant的结果相一致的。表1总结了图7中不断增强的压实作用对底床模型空间连续性的影响。底床模型中颗粒的半径在指数分布中位于10到20单位之间。胶结物厚度为1个单位。图15 从一个最终侵蚀孔隙空间模型（透明）的60步侵入渗透模拟（固体）中选取的6个图像。第4张图所示的是突破时的状态，也即实际方向中最适合的流体路径。孔隙结构网络结构，骨架目前，孔隙空间的体系结构定义为孔隙网络结构的完整3-D表征。Serra 和Adler把距离固体墙上面两个或更多点距离相等的点（三维象元）的集合定义为孔隙空间的骨架。因此，这个骨架可以认为是收缩变形核的空间展示，或者

23、孔隙网络的中心线。这条空间中心线将包含2条或者更多条线的交汇点。这些点即是网络节点（近似于孔隙体的中心），它们通过链接和其它网络节点连接或与喉道连接。在2-D图像分析中，通过薄化算法从双值图像中提取（孔隙）骨架。2-D薄化算法主要基于删除或保留经过分析的像素，其过程是分析它周围的8个像素的1和0的双值形式规律。在2-D（空间）中，一个像素周围的0和1的分布组合数目是：28=256在3-D（空间）中，一个像素周围的26个0和1的分布组合数目为：226=这样使得定义一个完全的3-D薄化算法规律的集合非常枯燥。Thovert等成功的建立了一个适用于他们模型的3-D薄化算法。本文试图建立一个相似的薄化

24、算法。但是，可视化检测表明，在一些复杂的象元交汇处，薄化算法可能会给孔隙网络骨架带来人为的空隙。对流体流来说，这可能是灾难性的，因为它会导致人为的液压回路。然而砂岩建模系统的逻辑解决了3-D（孔隙）骨架或薄化问题。与薄化函数相反的是膨胀。因此，孔隙网络薄化的最终结果和相应的颗粒网络膨胀的最终结果是精确的一致。如图16左所示，首先数字化每一个颗粒，然后数字化的底床阵列最终膨胀为图16右所示的颗粒形态。图16 数字化颗粒的底床最终膨胀为右图所示的数字底床。不同灰度表示不同数目的颗粒。在最初的矿物基质空间中，两个最终膨胀的颗粒间的边界面是粒间（点-点）接触的。但是，在最初的孔隙空间中，这个面包含了孔

25、隙网络骨架。因此，在两个或更多膨胀颗粒相遇的位置可以找到孔隙网络的骨架。那些来自于3个或更多最终膨胀颗粒的相邻象元的空间点定义了这条线（图16）。这条空间线上来自于4个或以上不同的最终膨胀颗粒的相邻象元的点即为孔隙网络的节点。图17 图11所示的孔隙网络骨架孔隙几何学如图11所示，孔隙空间的几何形态细节是相当不规则的。为了模拟孔隙网络中流体的流动，有必要把不规则几何形状的孔隙和喉道替换为等效的简化几何形状，这种简化几何形状保留与流体流动相关的本质特征。这可以通过记录下面的孔隙空间描述信息得以实现：（）孔隙体的最大内切圆半径（由体积计算得到），（）连接两个孔隙体的孔隙空间的最小内切半径（排驱压力

26、和液压传导系数），（）孔隙形状（为了描述在相同孔隙中同时流动的两种或更多流体）。骨架（图17）是描绘节点间孔隙喉道的基础。对于每一个孔隙喉道的骨架象元，孔隙喉道壁被映射在一个与喉道局部方向垂直的平面上。通过一个旋转矢量测量喉道半径，即在平面内以10度为间隔测量36次骨架象元到喉道壁之间的距离。由此可以记录下喉道最窄处的面积、周长以及内切半径（图18）。每一个孔隙体的体积通过平均半径计算得到，而平均半径则来源于骨架节点确定的几何中心。分析孔隙体几何形态的方式与分析喉道的几何形态类似，即在一个平面内以10度为间隔旋转半径矢量，共测量36次来研究孔隙壁。以10度为间隔旋转平面18次，即涵盖了整个空间

27、角度。研究向量所抵达的孔隙壁上的点的半径和坐标都会被记录。过长的半径被认为是测量到了喉道里面，这个半径值需要删除，而用一个平均值来代替。对于每一个孔隙体，记录其体积和内切球半径。图18 对孔隙喉道进行几何学分析的2-D示意图。上面的纵切面表示的是几何学分析里的骨架线和面。下面的横切面为了说明基于骨架的几何学分析在每一个横切面中是如何实现的。真实孔隙介质的孔隙结构的一个关键特征就是其孔隙具有棱角。孔隙棱角保留了润湿流体，且允许两种或多种流体同时通过孔隙和喉道。目前，每一个孔隙体和喉道的横截面形状用一个无量纲的形状因子G来描述。定义如下：G=A/p2 （2）式中A是孔隙体或喉道横截面的面积，p是相

28、应的周长。形状因子的范围可以从裂缝孔隙的0到等边三角形孔隙的0.0481。上述分析技术完整且定量的表征了孔隙网络结构及孔隙体和喉道的几何形态。图19展示的是一个简化后的球-柱数字化孔隙网络的画像，该画像是从2160毫达西的Bentheimer砂岩的薄片数据建立的。数字化孔隙网络的流动边界条件定义为：（流体）流动方向控制哪些面是入口和出口。在流入面（和流出面）上的所有节点和链接，均连接到这个面外部的一个虚拟节点上。现在，其余的四个面没有流动边界（正在进行建立周期性的流动边界）。流动方向可能平行于X,Y或Z轴。图19 一个具有1262个节点的Bentheimer砂岩孔隙网络的球柱形示意结构。球形的

29、半径与孔隙体的平均半径成正比，圆柱体的半径与孔隙喉道的最小半径成正比。立方体边长1.7mm。流体流的网络模型我们前面已经描述了重新构建的3-D砂岩网络结构中流体驱替的模拟，但为了本文的完整性现在简要复述一下。流动假设运用如下假设模拟孔隙网络中的两相（流体）流动：1. 流动形式都是层流（Poiseuille流动）；全部流体都是牛顿流，不可压缩，不混溶。注入流体是非润湿的。2. 孔隙体和喉道可能同时被一种或两种流体占据。假如是后者，则非润湿相流体是主体流，而润湿相流体以润湿膜的形式保存在孔隙和喉道的边角处。3. 孔隙体或节点被划分为两组完整节点和界面节点。完整节点（未被侵入的节点）被一种主体流占据

30、，而界面节点（被侵入节点）被两种流体占据。界面节点包含的流体随时间而改变。4. 流体压力仅在孔隙体中测定。孔隙体被认为足够大，可以使通过侵入界面的毛管压力可以忽略不计。目前的讨论，都是假设孔隙网络是强水润湿的。混合润湿系统不在我们当前的讨论范围。孔隙体和孔隙喉道的侵入当侵入的非润湿相流体位于孔隙I中，而与孔隙I相连的吼道IJ中充满了润湿相流体，毛管压力阻止非润湿相流体自发的进入喉道。只有非润湿相流体压力（Pnw）超过润湿相流体压力（Pw）的值达到初始毛管压力Pc时，非润湿相流体才可以通过喉道。Pnw，I Pw，J Pc，IJ = （3）式中是界面张力，rIJ是最窄喉道的内切圆半径。当满足上述条

31、件时，喉道才对侵入流体是开放的，流体可以立即侵入喉道。否则，喉道就是封闭的，喉道入口处的界面仍保持固定。当喉道IJ开放后，侵入界面立即进入到邻接的节点J，节点J变为界面节点。如果界面节点被非润湿相流体占据，就会发生聚合，则驱替是完全的。如果界面节点被润湿相流体占据，非润湿相流体从节点驱替润湿相流体。界面节点充填或倒空的速率取决于两种流体的粘滞压力梯度。流体传导率为了计算网络中的粘滞压力下降，我们首先必须赋值孔隙体和喉道的液压传导率。孔隙体和喉道的液压传导率取决于对流体开放的横截面的几何形态和优势的毛管压力。如图20所示，当孔隙体和喉道中同时存在两种流体时，受毛管压力支配，非润湿相流体占据孔隙主

32、体空间，而润湿相流体只在边角处得以保留。强润湿条件下，对两种流体流开放的横截面积可以用初等几何简单计算出来。= （4） = （5）式中r是孔隙体或喉道的内切半径；rw是膜界面的曲率半径，由杨拉普拉斯公式得出rw = （6）图20 三角形孔隙体或者喉道中润湿和非润湿相流体的示意性分布。当前处理中，我们可以简化占据孔隙体和喉道大部分空间的流体i的（压力）传导系数，即通过平均液压半径Ri的概念来表征，其定义如下：Ri = 0.5（r + rv，i） （7）式中rv，i是等效体积半径，定义为体积与占据孔隙体和喉道的流体i相等的圆柱体的半径，即rv，i= （8）如果孔隙体大部分被非润湿相流体占据，通过等

33、式（5）计算Ai；如果完全被润湿相流体占据，则Ai=r2/4G。流体i的（液压）传导系数gi 通过Poiseuille定理得到：gi = （9）式中i是流体粘度，l是孔隙体和喉道的长度。当早前非润湿相流体占据的孔隙体和喉道的边角主要部分存在润湿相流体时，我们采用gw = （10）式中是一个无量纲的流动阻力因子，它可解释导致接近孔隙壁的润湿相流体液压传导系数的减小。（孔隙体和喉道的）边角流体问题的有限元分析表明，取决于边角的几何形态以及润湿和非润湿相流体界面处的边界条件。例如，对于一个等边三角形的孔隙体或喉道，当流体界面没有流体流动时=5.3，而对于边界条件为没有压力的界面，=2.5。相位压力对

34、于层流来说，相互连通的孔隙体I和J中的流体i的体积流动速率Qi,IJ，可以通过以下公式得到Qi,IJ =gi,IJ (Pi,I - PiJ) （11）式中gi,IJ是孔隙体I和J中心间流体i的液压传导率。为了简化，我们假设gi,IJ是连接孔隙的喉道以及两个孔隙体本身液压传导率的调和平均值，即： （12）由于流体是不可压缩的，连续性要求每一个节点满足： （13）式中J涉及所有到节点I 的孔隙喉道。以公式（11）和（13）为基础，加上合适的初始值和边界条件，可以建立节点或孔隙中心压力的一系列线性代数公式，用矩阵形式表示如下GP=b （14）式中G是一个稀疏阵，其元素包含g的元素，P是（未知的）节点

35、压力，b是一个源向量（除流体进入、离开和界面节点以外，其值为0）。节点压力通过迭代共轭梯度法得到。驱替的前提运用准静态近似法处理。在t的时间内，与每种流体的计算的压力场相对应，两种流体在网络中以恒速流动。时间步长的选择必须满足在每一个步长内只能一个界面节点完全填满。因为可能有许多界面节点同时被侵入，我们需识别出所有接受填充的节点，并计算填充每个节点所需要的时间。确定最小侵入时间t，每一个界面节点的非润湿相流体的饱和度由下面公式得到：SS （15）式中VI是节点I的体积，并且总和只包括与界面节点I相连的非润湿相填充的孔隙喉道。在一些界面节点中，Snw可能会减小。如果Snw减小到0，非润湿相流体从

36、连接（孔隙体）的孔隙喉道自发缩回，该节点也从界面节点中除去。初始值和边界条件最初，网络中只有润湿相流体。非润湿相流体通过一个非润湿流体储层以恒速注入，此储层通过孔隙喉道与沿着网络入口侧的每一个孔隙体相接。在网络中润湿相压力任意设定为0 的地方，流体从网络相反方向的流出面逃逸。沿着与流动主方向平行的侧面，非渗透边界条件被加强。因为贯穿整个网络的润湿流体是连续的，公式（14）必须应用于所有孔隙体，而不是仅位于出口处压力任意设为0的孔隙体。非润湿流体也是连续的，但在每一瞬间它只占据了网络的一部分。因此，公式（14）只适用于被非润湿流体占据的孔隙体。两相流体压力的解决是反复的，因为在大量流体压力平衡的

37、界面节点，单相的流体压力与界面节点的压力相互耦合。有必要指定一个先验性估计：对于侵入的非润湿流体，哪些孔隙喉道是开启的，哪些是封闭的。对于每一个时间步长t，步骤如下：1. 每一个包含流体的节点的每相流体的压力运用公式（14）进行计算。2. 计算得到的压力再用来计算每一个封闭喉道的毛管压力。如果毛管压力超过门限压力，则喉道是开启的，且压力解决方案需要重新计算以确保喉道保持开启。此步骤不断重复，直到计算得到的毛管压力与门限毛管压力标准完全一致。3. 每一个节点的流动速率运用公式（11）计算。流速用来计算时间步长t，以确保在这个时间步长内仅一个界面节点完全填满。4. 所有界面节点中的饱和度通过公式（

38、15）更新。5. 对下一个时间步长重复以上步骤。宏观特性为了使边界效应最小化，宏观迁移特性的计算在孔隙网络的中心部位进行。当网络中饱和润湿相流体时，我们指定一个穿过网络的固定压力降，然后计算穿过研究域的（流体）流速Q和压力降P。绝对渗透率K通过达西定律得到：K = （16）式中A是研究区域的横截面积，L是其长度。相对渗透率kr运用相似方式计算得到。在驱替的各个阶段，我们计算经过研究域的非润湿和润湿流体的平均压力降。同时计算经过研究域的非润湿和润湿流体的饱和度和总流速。流体i的相对渗透率由达西定律计算得到Kr，i（S）= （17）研究域中的有效毛管压力定义为非润湿和润湿流体平均压力的差值，并且在模拟的各个阶段分别计算。当网络中饱和导电性卤水时，其电导率r可以通过公式（11-14）得到，但其中电流I取代Q