CN104331579A - 一种低渗透储层原油边界层的模拟方法 - Google Patents

Abstract

一种低渗透储层原油边界层模拟方法，包括步骤如下：1)CT扫描法构建数字岩心；2)对上述已构建的数字岩心运用Delaunay三角剖分算法进行剖分，得到三维几何模型；3)考虑边界层影响的条件下运用有限元方法对流体在上述剖分后的三维描述几何模型进行数值模拟，即对流体在孔道中微观流动的不可压缩方程N-S方程进行求解，得到流体流动的速度场和压力分布场。本发明利用comsol模拟，可以节省人力物力及时间。本发明利用规则孔网模型，可以准确模拟真实岩心的孔喉情况，更加符合实际情况。本发明对数据结果进行处理，可以准确计算边界层厚度，更接近真实岩心的原油边界层厚度。

Description

一种低渗透储层原油边界层的模拟方法

技术领域

本发明涉及一种低渗透储层原油边界层的模拟方法，属于石油天然气开发的模拟实验的技术领域。

背景技术

目前的研究表明：现有的采收率机理研究多数都基于宏观渗流理论，宏观渗流理论都是以达西渗流方程为基础，所以在石油开发领域里应用以达西方程为基础所建立的渗流方程只能宏观描述油藏流体的渗流特征，无法给出更精确的孔隙级别上的描述。特别是在低渗透储层的研究中，许多学者已经通过低渗透岩心驱替实验和开发实践表明：低渗透储层中，油气渗流不符合达西定律。为表征低渗透储层的渗流特征，许多学者提出了各种数学模型，其中主要有拟启动压力梯度模型，分段模型，连续模型。上述数学模型都是基于实验数据的函数拟合，对渗流特征的解释都有所欠缺。储层的多数宏观性质，如毛管压力、渗透率、相对渗透率等实际上取决于岩石的微观结构。因此，仅仅通过对宏观性质的分析展开的研究只是停留在事物表面，无法对其内部的本质有所认识，由此提出的提高采收率的技术手段也有其本质上的局限性。

若使提高原油采收率技术大幅度提高，仅在宏观层次上进行理论研究和技术开发是远远不够的，必须从微观层面上突破，建立完善合理的微观渗流理论。通过研究多孔介质内部孔隙空间的空间分布、流体在孔道空间中的分布特点、相互作用机理等能决定宏观渗流现象的本质问题，才能从根本上认识宏观与微观的联系，为提高采收率找到更好的技术手段。数字岩心技术是开展微观渗流研究所必需的平台，是解决上述问题的关键所在。借助数字岩心技术能够充分认识多孔介质的孔隙结构，流体在多孔介质中流动的规律和相互作用机理等本质问题。同时CT技术在数字岩心构建中的运用是对数字岩心技术的完善和强化，通过CT技术建立的数字岩心具有最接近真实等特点。

对于中低渗储层，控制岩样渗流能力的主要因素是喉道特征而不是孔隙特征，因此，深入认识低渗透油藏的微观孔隙结构，是提高低渗储层动用率的关键。

致密砂岩介质的特点：岩性致密，孔隙尺寸细小，渗透率低；孔隙喉道几何形状复杂，毛细管现象严重；流体与油藏介质作用机理复杂，非牛顿现象相对显著。

目前，致密砂岩油藏已成为油气开采的重要对象，但是由于其中流体的渗流特征(启动压力梯度、变化的视渗透率等)不符合达西规律，很多微观机理问题尚未解决。因此，通过数字岩心技术对储层微观孔隙结构进行评价，并模拟边界层效应，进而分析微观孔隙结构和边界层效应对非达西渗流曲线和油水两相渗流曲线的影响。

目前，国内大多采用实验方法测量边界层厚度，主要包括测量玻璃毛细管的有效半径求得平均边界层厚度和测量岩石毛细管模型求取平均边界层厚度。这两种模型的不足之处分别在于，玻璃毛细管模型方法基于达西定律，对于低渗透、特低渗透油藏来讲，流体在岩心孔喉内的流动不符合达西定律；同时用玻璃毛细管代表岩心孔喉，不能模拟流体在孔喉内流动时与岩心表面的相互作用。岩石毛细管模型根据岩心驱替前后的原油质量差，推算边界层的平均厚度，由于实验条件的影响，误差较大，并且模型过于理想化。

发明内容

针对现有技术的不足，本发明提供一种低渗透储层原油边界层的模拟方法。该模拟方法可以模拟不同粘度的原油，在不同的压力梯度下，流过不同的孔喉半径的情况，借助数据处理可以计算原油边界层的有效厚度。

本发明的技术方案如下：

一种低渗透储层原油边界层模拟方法，包括步骤如下：

1)CT扫描法构建数字岩心包括步骤(a)-(c)：

(a)通过CT扫描获取岩样的多方位CT投影图像数据；其中岩样的形状对CT投影图像没有影响，并且，本发明使用的CT是MicroXCT-400高分辨3DX射线显微镜，可对岩样进行360度成像，所以拍摄岩样的角度对于岩样投影图像没有影响。理论上，本发明使用的MicroXCT-400CT扫描设备支持样品尺寸达300mm，重量达15Kg；利用CT进行扫描，扫描后的单张CT投影图像数据被图像获取软件自动捕获并存储；将所述岩样沿其轴向旋转固定角度后再次对其CT扫描，扫描后的单张CT投影图像数据被图像获取软件自动捕获并存储，直至所述的岩样沿轴向累计旋转360°后结束采集CT投影图像数据；

CT的成像原理：当X射线穿过任何物质时，它会与物质的原子相互作用而引起能量衰减，并且物质不同的组成成分对X射线具有不同的吸收系数(即衰减系数)，反而言之，通过测量物质对X射线的吸收系数可以判定物质的组成成分。当一束X射线穿过物体时，它所穿透路径上的所有物质对X射线吸收系数的总和都将反映在对X射线强度的测量结果中，如(1.1)所示：

$I=I_o{e}^{-\underset{i}{\mathrm{\Σ}}{\mathrm{\μ}}_i{x}_i}---\left(1.1\right)$

式中，I_o为X射线的初始强度，I为X射线穿过物体后的强度，即X射线被衰减后的强度，i表示在该X射线穿过的路径上物质中的某一组成成分，μ_i、x_i分别为第i组分对X射线的衰减系数和该组分在X射线当前途经路径上的长度；

CT的成像原理正是建立在这一基础上，通过对穿过物体截面的X射线进行测量，之后，采用一定的重建方法计算出与物体断层空间位置一一对应的吸收系数，从而恢复物体截面的结构信息；

本发明模拟计算时，岩样的具体尺寸要根据岩样的属性确定，以达到最好的技术效果。为保证岩心CT图像能够清晰反映岩心内部微观结构特征，实验岩心尺寸不可过大，原因主要有两个方面：一方面，由于岩心密度较大，射线穿透岩心过程中被大幅衰减导致其到达探测器***时信号已十分微弱，若岩心尺寸过大则可能导致射线完全被屏蔽；另一方面，由于软件***在求解时所采用的最大格点数量是一定的，岩心尺寸增大意味着每个格点对应的实际尺寸增大，即图像的分辨率降低了。对现有的工业用CT设备而言，要获得较清晰的CT图像，岩心柱体直径最大为厘米级；

(b)利用图像重建方法将所述CT投影图像数据转换为所述岩样的灰度图像；

CT成像的核心是如何由投影数据重建图像。图像重建的实质是由采集后的CT投影数据求解图像矩阵中各像素的吸收系数，然后重新构造图像的过程；图像重建问题的求解方法根据其特点可分为两大类：第一类是变换法(也叫解析法)，其特点是先在连续域解析处理，最后离散化利用计算机计算；第二类是迭代法(也叫代数重建法、级数展开法、优化技术等)，其特点是直接进行离散化分析得到数值解。相比而言，变换法优点突出：实现简单，速度快，对足够精确的投影数据能获得很好的重建质量。因此，目前实用CT***中广泛采用变换法。

所述的图像重建方法为变换法，包括不含反投影步骤的Fourier重建方法FRM和基于反投影的卷积或滤波方法：前者以直接Fourier变换方法DFM为代表，后者以卷积或滤波反投影C/FBP方法为代表；

DFM的基本思想是根据辐射源与成像目标之间的作用机理，寻找到目标函数f(x,y)的二维Fourier变换(FT)与投影函数的一维FT之间的关系，利用这种关系，尽管目标函数f(x,y)是未知的，但它的二维FT却能通过不同方位上投影函数的一维FT得到，一旦获得了f(x,y)的所有投影，也就获得了它的二维FT，f(x,y)则可由二维FT的反变换来求出；

基于反投影的卷积或滤波方法C/FBP方法源于Radon求反公式的直接实现；C/FBP公式利用FT和投影定理导出；设p(x',θ)的最高空间频率是B；

C/FBP的数学公式写成

$f(x,y)={\∫}_0^{\mathrm{\π}}q(x^{\′},\mathrm{\θ})\mathrm{d\θ}---\left(1.2\right)$

其中，

$q(x^{\′},\mathrm{\θ})={\∫}_{-\∞}^{\∞}\left|R\right|P(R,\mathrm{\θ})\exp \left(j2\mathrm{\π}{x}^{\′}R\right)\mathrm{dR}---\left(1.3\right)$

或

q(x',θ)＝c(x')*p(x',θ)   (1.4)

式(1.2)、(1.3)、(1.4)中，f(x,y)为图像重建的目标函数，即岩样对射线吸收系数的分布函数，x，y分别为岩样扫描断面上某点的直角坐标；P(R,θ)、p(x',θ)为f(x,y)沿θ方向上的投影函数，即CT探测器检测到的射线强度信号，x'为岩样扫描断面上某点(x，y)在旋转坐标系中的横坐标，R为该点到旋转中心，即坐标原点，的距离；c(x')为卷积函数，它由滤波器|R|(|R|≤B)的Fourier反变换求得；

式(1.3)、(1.4)中的q(x',θ)分别称为投影P(R,θ)、p(x',θ)的滤波投影和卷积投影；由式(1.2)、(1.3)实现的图像重建方法称为滤波反投影方法FBP，由式(1.2)、(1.4)实现的图像重建方法称为卷积反投影方法CBP；因此，如获得了图像的投影，则通过卷积或滤波反投影方法反向求解得到岩样断面对射线吸收系数的分布函数f(x,y)，将f(x,y)采用灰度图像的方式表达后就得到岩样断面的扫描图像；

(c)采用图像二值分割方法分离所述岩样的灰度图像中孔隙空间和岩样骨架，建立数字岩心；

在灰度图像中，孔隙空间和岩石骨架通过图像分割方法对其进行划分：图像分割在于分割阈值的选取，如果阈值选取不合理，很容易将目标对象理解成为背景图像，或将背景理解为目标对象，这样在对目标图像进行特征分析时就无法有效分析目标特征；

阈值定义如下：假设F(i,j)表示灰度图像中图像像素灰度值，经分离后灰度值大于阈值T的像素表示目标对象，即孔隙，用1表示；其它的像素点表示图像背景，即骨架，用0表示，则由下式即可获得目标对象：

$g(i,j)=\left\{\begin{array}{cc}1;& F(i,j)\≥T\\ 0;& \mathrm{else}\end{array}\right.---\left(1.5\right)$

式中，T即为阈值，又称门限值；

基于岩心孔隙度的图像二值分割方法

岩心孔隙度是开展油层物理实验时必须测定的基本物理量；在现有的实验条件下，借助气体膨胀法、液体饱和法等都可以方便准确的测量得到；由于孔隙度给出了岩心中孔隙空间所占的份额，因此，在已知的情况下，把孔隙度作为约束条件结合到岩心CT图像的分割过程中无疑会提高图像的分割质量；

设岩样的孔隙度为图像的最大、最小灰度值分别为I_Max、I_MIN，灰度值为i的像素数为p(i)，则满足式(1.6)的灰度值k^*即所求分割阈值：

基于岩样孔隙度的图像二值分割方法实际上是选择合适的分割阈值，使得通过该阈值分割得到的岩心图像中孔隙像素数占总像素数的比例等于岩心孔隙度。由于像素灰度值由0到255区间中的整数表示，因此，分割得到的图像孔隙度可能与岩心实验孔隙度存在细微差别，但这不影响图像分割质量。

通过对实验图像进行分割测试发现，基于岩心孔隙度的图像二值分割方法能够得到高质量的分割结果，因此，在岩心实验孔隙度已知的情况下，可以借助该方法对岩心CT图像进行二值分割。然而，并非所有用于CT实验研究的岩心都有实验孔隙度，当没有孔隙度资料时，需借助其它图像二值分割方法来分割岩心CT图像。在图像二值分割方法中，最大类间距法因原理简单、程序设计方便且能够获得较好的分割效果而被广泛应用。

当没有孔隙度时，利用最大类间距法分割灰度图像；将待分割灰度图像中所有像素点的灰度值所构成的总体视为一个集合，假定把灰度值的集合用阈值分成两组，则当两组的平均值的方差(组间方差)和各组的方差(组内方差)之比达到最大时，其所对应的分割阈值即为合理值；所述岩样内部的孔隙空间和骨架各成一组，每一组内部的差异很小而两组间的性质差异很大；当采用合理的阈值将两组分割开后，两组间的性质差异与组内差异之比将达到极大值，可见，岩样CT扫描投影图像本身的特点与最大类间距法对图像的分割思路十分吻合，这也是最大类间距法分割岩心CT图像能够获得良好分割结果而被广泛采用的重要原因；

最大类间距法选取阈值的原理如下。设图像的最大、最小灰度值分别为I_Max、I_MIN，灰度值为i的像素数为p(i)；灰度阈值为k，它将像素点分为两组：组1(假设为孔隙空间)和组2(假设为岩石骨架)。分别统计并计算两组的像素数N₁(k)、N₂(k)，平均灰度值μ₁(k)、μ₂(k)和方差及全部像素的平均灰度值μ_T，则组内方差及组间方差分别为：

${\mathrm{\σ}}_W^2=N_1{\mathrm{\σ}}_1^2+N_2{\mathrm{\σ}}_2^2---\left(1.7\right)$

${\mathrm{\σ}}_B^2=N_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2+N_2{({\mathrm{\μ}}_2-{\mathrm{\μ}}_T)}^2---\left(1.8\right)$

式中：

${\mathrm{\μ}}_1\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{k}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_1\left(k\right)}$ ${\mathrm{\μ}}_2\left(k\right)=\frac{\underset{i=k+1}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_2\left(k\right)}$

${\mathrm{\μ}}_T\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_1\left(k\right)+N_2\left(k\right)}$

${\mathrm{\σ}}_1^2\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{k}{\mathrm{\Σ}}}{[i-{\mathrm{\μ}}_1\left(k\right)]}^2\·p\left(i\right)}{N_1\left(k\right)}$ ${\mathrm{\σ}}_2^2\left(k\right)=\frac{\underset{i=k+1}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}{[i-{\mathrm{\μ}}_2\left(k\right)]}^2\·p\left(i\right)}{N_2\left(k\right)}$

选择最佳灰度值k^*，使其满足(1.9)式，则k^*即所求阈值。

$f\left(k^{*}\right)=\max \{f\left(k\right)={\mathrm{\σ}}_B^2\left(k\right)/{\mathrm{\σ}}_W^2\left(k\right)\}---\left(1.9\right);$

采用最大类间距法分割图像的算法流程，如图1所示；

2)对上述已构建的数字岩心运用Delaunay三角剖分算法进行剖分，得到三维几何模型；

其中，所述Delaunay三角剖分算法是也是一种面绘制方法，它的优点是在构建几何模型的过程中能通过控制几何单元的参数，使得形成的几何模型更利于之后的分析和计算。由于Delaunay三角剖分算法的诸多优点，这种方法已经被广泛应用于有限元的分析、医学等多个领域^[1]。[1]丁永祥，夏巨谌.任意多边形的Delaunay三角剖分[J].计算机学报，1994，17(4):270-275。

Delaunay三角剖分是优化的三角剖分，具有很多优良特性^[2]。[2]张花艳.基于Delaunay三角剖分与场表示的曲面重建[D]:郑州大学，2011.

(i)利用Delaunay三角剖分建立的几何模型是唯一的；

(ii)Delaunay三角剖分之后的几何模型最外层边界形成的是一个凸多边形的边缘；

(iii)对形成三角网进行处理某一个顶点时只会对临近三角形区域产生影响，而对其它的三角形没有影响；

(iv)Delaunay三角形的内部没有任何顶点存在，即三角形的内部是中空的；

(v)Delaunay三角剖分得到的三角形的最小角与其它三角剖分相比是最大的；

(vi)如果可以相互交换任意两个相邻凸四边形的对角线，最小的内角角度并不会因此而变大；

(vii)Delaunay三角剖分以最近的三点构成三角形。

与二维空间的Delaunay三角剖分相比，三维空间的Delaunay三角剖分无论在难度还是复杂性上远超过二维空间。调用CGAL算法库进行三维空间中的Delaunay剖分是很有效率的一种方法。Computational Geometry AlgorithmsLibrary，简称CGAL，是由以色列和欧盟的几个国家和地区的8个研究机构联合开发的一个C++算法库^[3]。[3]Fabri A，Pion S.CGAL:the computational geometryalgorithms library:proceedings of the Proceedings of the 17th ACM SIGSPATIAL internationalconference on advances in geographic information systems，2009[C].ACM.

CGAL提供了可直接用于三维图像的几何模型重构的3D面网格剖分包，该剖分包基于带限制的Dlaunay三角剖分算法。

3D面网格剖分包剖分算法进行三维图像几何重构的主要步骤如下:

1)在物体的表面选择一部分做样品点计算；

2)利用第1步计算的空间体中提取出一个曲面；

3)***点到提取的曲面网格中，并对形成的新的空间体进行优化处理；

4)判断新的空间体网络是不是满足构建准则，如不满足继续重复优化；

5)如满足条件，则继续之后的处理，一次完成所有的网格重构。

由此可见，构建准则是决定剖分质量和结果的关键，并影响网格最终的形状和尺寸。

3)考虑边界层影响的条件下，运用有限元方法对流体在上述剖分后的三维描述几何模型中进行数值模拟，即对流体在孔道中微观流动的不可压缩方程N-S方程进行求解，得到流体流动的速度场和压力场分布；

a.有限元法的基本思想：先将研究对象的连续求解区域离散为一组有限个且按一定方式相互联结在一起的单元组合体。由于单元能按不同的联结方式进行组合，且单元本身又可以有不同形状，因此可以模拟成不同几何形状的求解小区域；然后对单元(小区域)进行力学分析，最后再整体分析。这种化整为零,集零为整的方法就是有限元的基本思路^[4]。[4]陈锡栋，杨婕.有限元法的发展现状及应用[J].中国制造业信息化，2010，39(11):6-8

有限元法是常用的计算流体力学模拟方法之一，有着严格的数学理论，并且不需要对孔隙空间进行简化，但有限元的计算是在几何模型的基础上进行的，因此需要通过某种方法将我们得到的数字岩心进行几何重构。同时为了减小三维几何重构过程中的计算量，我们通常还要采取一定的岩心孔隙处理措施，比如特征单元体分析、孤立孔隙的去除等。

其中所述N-S方程：

宏观连续模型基于连续假设，流体遵循质量守恒、动量守恒和能量守恒，其描述方程为Navier-Stokes方程(N-S方程)：

$\frac{\∂\mathrm{\ρ}}{\∂t}+\▿\·\left(\mathrm{\ρu}\right)=0---\left(1.10\right)$

$\frac{\∂\left(\mathrm{\ρu}\right)}{\∂t}+\▿\·\left(\mathrm{\ρuu}\right)=\▿\·\mathrm{\σ}---\left(1.11\right)$

$\frac{\∂\left(\mathrm{\ρe}\right)}{\∂t}+\▿\·\left(\mathrm{\ρue}\right)=\mathrm{\σ}:\▿u-\▿\·q---\left(1.12\right)$

式中，ρ为流体宏观密度；u为流体速度；T为温度；e为内能；σ为应力张量；q为热通量；上述方程组并不封闭，要求解压力、速度等参数必须补充本构关系、状态方程等关系。宏观连续模型是目前发展比较成熟，应用最为广泛的方法。

求解的N-S方程的具体形式如下：

$\begin{array}{c}\mathrm{\ρ}\frac{\∂u}{\∂t}+\mathrm{\ρ}(u\·\▿)u=\▿\·[-\mathrm{\ρI}+\mathrm{\μ}(\▿u+{(\▿u)}^T)]\\ \mathrm{\ρ}\▿\·u=0\end{array}---\left(1.13\right)$

其中，ρ为流体的密度；u为流体的速度；μ为流体的粘度；Ι为单位矩阵；

边界条件为：

$\begin{array}{c}p=p_0\\ [\mathrm{\μ}(\▿u)+{(\▿u)}^T]\·n=0\end{array}---\left(1.14\right)$

其中，p为压力，p₀为初始压力；

所述模拟步骤包括步骤①-⑤：

①导入剖分好的三维几何模型；

②整体参数设置，包括液体密度，原油粘度及压力；

③选择对应的岩样出口、入口及壁面，进行计算，得到岩样出口和入口的速度；

④对岩样出口和入口的速度进行积分，得到岩样出口和入口的流量；

⑤对模拟的结果进行绘图，得到岩样速度场和压力场分布。

根据上面得到的数据，可以绘制速度场、压力场分布及流线等，通过处理可以看出流体在数字岩心内的流动情况。

通过上述方法，可以对比不同类型岩石中边界层的影响状况。

边界层的厚度和粘度对微观孔隙中流体的渗流有很大影响。边界层的阻力作用在特征尺寸比较小的岩石中更为明显，主要表现为速度下降更快，压力减小更快；边界层对低渗介质的影响更为明显，在低渗透介质中研究边界层更为重要。

相对于现有技术，本发明具有如下有益效果：

1、本发明利用Comsol模拟，可以节省人力物力及时间。

2、本发明利用数字岩心技术，可以准确模拟真实岩心的孔喉情况和压力、流体分布，更加符合实际情况。

3、本发明对数据结果进行处理，可以准确计算边界层厚度，更接近真实岩心的原油边界层厚度。

附图说明

图1为本发明采用最大类间距法分割图像的算法流程图。

图2为本发明模拟计算方法步骤的流程图。

具体实施方式

下面结合实施例对本发明做详细的说明，但不限于此。

实施例、

一种低渗透储层原油边界层模拟方法，包括步骤如下：

1)CT扫描法构建数字岩心包括步骤(a)-(c)：

(a)通过CT扫描获取岩样的多方位CT投影图像数据；其中岩样的形状对CT投影图像没有影响，并且，本发明使用的CT是MicroXCT-400高分辨3DX射线显微镜，可对岩样进行360度成像，所以拍摄岩样的角度对于岩样投影图像没有影响。理论上，本发明使用的MicroXCT-400CT扫描设备支持样品尺寸达300mm，重量达15Kg；利用CT进行扫描，扫描后的单张CT投影图像数据被图像获取软件自动捕获并存储；将所述岩样沿其轴向旋转固定角度后再次对其CT扫描，扫描后的单张CT投影图像数据被图像获取软件自动捕获并存储，直至所述的岩样沿轴向累计旋转360°后结束采集CT投影图像数据；

CT的成像原理：当X射线穿过任何物质时，它会与物质的原子相互作用而引起能量衰减，并且物质不同的组成成分对X射线具有不同的吸收系数(即衰减系数)，反而言之，通过测量物质对X射线的吸收系数可以判定物质的组成成分。当一束X射线穿过物体时，它所穿透路径上的所有物质对X射线吸收系数的总和都将反映在对X射线强度的测量结果中，如(1.1)所示：

$I=I_o{e}^{-\underset{i}{\mathrm{\Σ}}{\mathrm{\μ}}_i{x}_i}---\left(1.1\right)$

式中，I_o为X射线的初始强度，I为X射线穿过物体后的强度，即X射线被衰减后的强度，i表示在该X射线穿过的路径上物质中的某一组成成分，μ_i、x_i分别为第i组分对X射线的衰减系数和该组分在X射线当前途经路径上的长度；

CT的成像原理正是建立在这一基础上，通过对穿过物体截面的X射线进行测量，之后，采用一定的重建方法计算出与物体断层空间位置一一对应的吸收系数，从而恢复物体截面的结构信息；

本发明模拟计算时，岩样的具体尺寸要根据岩样的属性确定，以达到最好的技术效果。为保证岩心CT图像能够清晰反映岩心内部微观结构特征，实验岩心尺寸不可过大，原因主要有两个方面：一方面，由于岩心密度较大，射线穿透岩心过程中被大幅衰减导致其到达探测器***时信号已十分微弱，若岩心尺寸过大则可能导致射线完全被屏蔽；另一方面，由于软件***在求解时所采用的最大格点数量是一定的，岩心尺寸增大意味着每个格点对应的实际尺寸增大，即图像的分辨率降低了。对现有的工业用CT设备而言，要获得较清晰的CT图像，岩心柱体直径最大为厘米级；

(b)利用图像重建方法将所述CT投影图像数据转换为所述岩样的灰度图像；

CT成像的核心是如何由投影数据重建图像。图像重建的实质是由采集后的CT投影数据求解图像矩阵中各像素的吸收系数，然后重新构造图像的过程；图像重建问题的求解方法根据其特点可分为两大类：第一类是变换法(也叫解析法)，其特点是先在连续域解析处理，最后离散化利用计算机计算；第二类是迭代法(也叫代数重建法、级数展开法、优化技术等)，其特点是直接进行离散化分析得到数值解。相比而言，变换法优点突出：实现简单，速度快，对足够精确的投影数据能获得很好的重建质量。因此，目前实用CT***中广泛采用变换法。

所述的图像重建方法为变换法，包括不含反投影步骤的Fourier重建方法FRM和基于反投影的卷积或滤波方法：前者以直接Fourier变换方法DFM为代表，后者以卷积或滤波反投影C/FBP方法为代表；

DFM的基本思想是根据辐射源与成像目标之间的作用机理，寻找到目标函数f(x,y)的二维Fourier变换(FT)与投影函数的一维FT之间的关系，利用这种关系，尽管目标函数f(x,y)是未知的，但它的二维FT却能通过不同方位上投影函数的一维FT得到，一旦获得了f(x,y)的所有投影，也就获得了它的二维FT，f(x,y)则可由二维FT的反变换来求出；

基于反投影的卷积或滤波方法C/FBP方法源于Radon求反公式的直接实现；C/FBP公式利用FT和投影定理导出；设p(x',θ)的最高空间频率是B；

C/FBP的数学公式写成

$f(x,y)={\∫}_0^{\mathrm{\π}}q(x^{\′},\mathrm{\θ})\mathrm{d\θ}---\left(1.2\right)$

其中，

$q(x^{\′},\mathrm{\θ})={\∫}_{-\∞}^{\∞}\left|R\right|P(R,\mathrm{\θ})\exp \left(j2\mathrm{\π}{x}^{\′}R\right)\mathrm{dR}---\left(1.3\right)$

或

q(x',θ)＝c(x')*p(x',θ)   (1.4)

式(1.2)、(1.3)、(1.4)中，f(x,y)为图像重建的目标函数，即岩样对射线吸收系数的分布函数，x，y分别为岩样扫描断面上某点的直角坐标；P(R,θ)、p(x',θ)为f(x,y)沿θ方向上的投影函数，即CT探测器检测到的射线强度信号，x'为岩样扫描断面上某点(x，y)在旋转坐标系中的横坐标，R为该点到旋转中心，即坐标原点，的距离；c(x')为卷积函数，它由滤波器|R|(|R|≤B)的Fourier反变换求得；

式(1.3)、(1.4)中的q(x',θ)分别称为投影P(R,θ)、p(x',θ)的滤波投影和卷积投影；由式(1.2)、(1.3)实现的图像重建方法称为滤波反投影方法FBP，由式(1.2)、(1.4)实现的图像重建方法称为卷积反投影方法CBP；因此，如获得了图像的投影，则通过卷积或滤波反投影方法反向求解得到岩样断面对射线吸收系数的分布函数f(x,y)，将f(x,y)采用灰度图像的方式表达后就得到岩样断面的扫描图像；

(c)采用图像二值分割方法分离所述岩样的灰度图像中孔隙空间和岩样骨架，建立数字岩心；

在灰度图像中，孔隙空间和岩石骨架通过图像分割方法对其进行划分：图像分割在于分割阈值的选取，如果阈值选取不合理，很容易将目标对象理解成为背景图像，或将背景理解为目标对象，这样在对目标图像进行特征分析时就无法有效分析目标特征；

阈值定义如下：假设F(i,j)表示灰度图像中图像像素灰度值，经分离后灰度值大于阈值T的像素表示目标对象，即孔隙，用1表示；其它的像素点表示图像背景，即骨架，用0表示，则由下式即可获得目标对象：

$g(i,j)=\left\{\begin{array}{cc}1;& F(i,j)\≥T\\ 0;& \mathrm{else}\end{array}\right.---\left(1.5\right)$

式中，T即为阈值，又称门限值；

基于岩心孔隙度的图像二值分割方法

岩心孔隙度是开展油层物理实验时必须测定的基本物理量；在现有的实验条件下，借助气体膨胀法、液体饱和法等都可以方便准确的测量得到；由于孔隙度给出了岩心中孔隙空间所占的份额，因此，在已知的情况下，把孔隙度作为约束条件结合到岩心CT图像的分割过程中无疑会提高图像的分割质量；

设岩样的孔隙度为图像的最大、最小灰度值分别为I_Max、I_MIN，灰度值为i的像素数为p(i)，则满足式(1.6)的灰度值k^*即所求分割阈值：

基于岩样孔隙度的图像二值分割方法实际上是选择合适的分割阈值，使得通过该阈值分割得到的岩心图像中孔隙像素数占总像素数的比例等于岩心孔隙度。由于像素灰度值由0到255区间中的整数表示，因此，分割得到的图像孔隙度可能与岩心实验孔隙度存在细微差别，但这不影响图像分割质量。

通过对实验图像进行分割测试发现，基于岩心孔隙度的图像二值分割方法能够得到高质量的分割结果，因此，在岩心实验孔隙度已知的情况下，可以借助该方法对岩心CT图像进行二值分割。然而，并非所有用于CT实验研究的岩心都有实验孔隙度，当没有孔隙度资料时，需借助其它图像二值分割方法来分割岩心CT图像。在图像二值分割方法中，最大类间距法因原理简单、程序设计方便且能够获得较好的分割效果而被广泛应用。

当没有孔隙度时，利用最大类间距法分割灰度图像；将待分割灰度图像中所有像素点的灰度值所构成的总体视为一个集合，假定把灰度值的集合用阈值分成两组，则当两组的平均值的方差(组间方差)和各组的方差(组内方差)之比达到最大时，其所对应的分割阈值即为合理值；所述岩样内部的孔隙空间和骨架各成一组，每一组内部的差异很小而两组间的性质差异很大；当采用合理的阈值将两组分割开后，两组间的性质差异与组内差异之比将达到极大值，可见，岩样CT扫描投影图像本身的特点与最大类间距法对图像的分割思路十分吻合，这也是最大类间距法分割岩心CT图像能够获得良好分割结果而被广泛采用的重要原因；

最大类间距法选取阈值的原理如下。设图像的最大、最小灰度值分别为I_Max、I_MIN，灰度值为i的像素数为p(i)；灰度阈值为k，它将像素点分为两组：组1(假设为孔隙空间)和组2(假设为岩石骨架)。分别统计并计算两组的像素数N₁(k)、N₂(k)，平均灰度值μ₁(k)、μ₂(k)和方差及全部像素的平均灰度值μ_T，则组内方差及组间方差分别为：

${\mathrm{\σ}}_W^2=N_1{\mathrm{\σ}}_1^2+N_2{\mathrm{\σ}}_2^2---\left(1.7\right)$

${\mathrm{\σ}}_B^2=N_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2+N_2{({\mathrm{\μ}}_2-{\mathrm{\μ}}_T)}^2---\left(1.8\right)$

式中：

${\mathrm{\μ}}_1\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{k}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_1\left(k\right)}$ ${\mathrm{\μ}}_2\left(k\right)=\frac{\underset{i=k+1}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_2\left(k\right)}$

${\mathrm{\μ}}_T\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}i\·p\left(i\right)}{N_1\left(k\right)+N_2\left(k\right)}$

${\mathrm{\σ}}_1^2\left(k\right)=\frac{\underset{i=I_{\mathrm{MIN}}}{\overset{k}{\mathrm{\Σ}}}{[i-{\mathrm{\μ}}_1\left(k\right)]}^2\·p\left(i\right)}{N_1\left(k\right)}$ ${\mathrm{\σ}}_2^2\left(k\right)=\frac{\underset{i=k+1}{\overset{I_{\mathrm{MAX}}}{\mathrm{\Σ}}}{[i-{\mathrm{\μ}}_2\left(k\right)]}^2\·p\left(i\right)}{N_2\left(k\right)}$

选择最佳灰度值k^*，使其满足(1.9)式，则k^*即所求阈值。

$f\left(k^{*}\right)=\max \{f\left(k\right)={\mathrm{\σ}}_B^2\left(k\right)/{\mathrm{\σ}}_W^2\left(k\right)\}---\left(1.9\right);$

采用最大类间距法分割图像的算法流程，如图1所示；

2)对上述已构建的数字岩心运用Delaunay三角剖分算法进行剖分，得到三维几何模型；

其中，所述Delaunay三角剖分算法是也是一种面绘制方法，它的优点是在构建几何模型的过程中能通过控制几何单元的参数，使得形成的几何模型更利于之后的分析和计算。由于Delaunay三角剖分算法的诸多优点，这种方法已经被广泛应用于有限元的分析、医学等多个领域^[1]。[1]丁永祥，夏巨谌.任意多边形的Delaunay三角剖分[J].计算机学报，1994，17(4):270-275。

Delaunay三角剖分是优化的三角剖分，具有很多优良特性^[2]。[2]张花艳.基于Delaunay三角剖分与场表示的曲面重建[D]:郑州大学，2011.

(i)利用Delaunay三角剖分建立的几何模型是唯一的；

(ii)Delaunay三角剖分之后的几何模型最外层边界形成的是一个凸多边形的边缘；

(iii)对形成三角网进行处理某一个顶点时只会对临近三角形区域产生影响，而对其它的三角形没有影响；

(iv)Delaunay三角形的内部没有任何顶点存在，即三角形的内部是中空的；

(v)Delaunay三角剖分得到的三角形的最小角与其它三角剖分相比是最大的；

(vi)如果可以相互交换任意两个相邻凸四边形的对角线，最小的内角角度并不会因此而变大；

(vii)Delaunay三角剖分以最近的三点构成三角形。

与二维空间的Delaunay三角剖分相比，三维空间的Delaunay三角剖分无论在难度还是复杂性上远超过二维空间。调用CGAL算法库进行三维空间中的Delaunay剖分是很有效率的一种方法。Computational Geometry AlgorithmsLibrary，简称CGAL，是由以色列和欧盟的几个国家和地区的8个研究机构联合开发的一个C++算法库^[3]。[3]Fabri A，Pion S.CGAL:the computational geometryalgorithms library:proceedings of the Proceedings of the 17th ACM SIGSPATIAL internationalconference on advances in geographic information systems，2009[C].ACM.

CGAL提供了可直接用于三维图像的几何模型重构的3D面网格剖分包，该剖分包基于带限制的Dlaunay三角剖分算法。

3D面网格剖分包剖分算法进行三维图像几何重构的主要步骤如下:

1)在物体的表面选择一部分做样品点计算；

2)利用第1步计算的空间体中提取出一个曲面；

3)***点到提取的曲面网格中，并对形成的新的空间体进行优化处理；

4)判断新的空间体网络是不是满足构建准则，如不满足继续重复优化；

5)如满足条件，则继续之后的处理，一次完成所有的网格重构。

由此可见，构建准则是决定剖分质量和结果的关键，并影响网格最终的形状和尺寸。

3)考虑边界层影响的条件下，运用有限元方法对流体在上述剖分后的三维描述几何模型中进行数值模拟，即对流体在孔道中微观流动的不可压缩方程N-S方程进行求解，得到流体流动的速度场和压力场分布；

a.有限元法的基本思想：先将研究对象的连续求解区域离散为一组有限个且按一定方式相互联结在一起的单元组合体。由于单元能按不同的联结方式进行组合，且单元本身又可以有不同形状，因此可以模拟成不同几何形状的求解小区域；然后对单元(小区域)进行力学分析，最后再整体分析。这种化整为零,集零为整的方法就是有限元的基本思路^[4]。[4]陈锡栋，杨婕.有限元法的发展现状及应用[J].中国制造业信息化，2010，39(11):6-8

有限元法是常用的计算流体力学模拟方法之一，有着严格的数学理论，并且不需要对孔隙空间进行简化，但有限元的计算是在几何模型的基础上进行的，因此需要通过某种方法将我们得到的数字岩心进行几何重构。同时为了减小三维几何重构过程中的计算量，我们通常还要采取一定的岩心孔隙处理措施，比如特征单元体分析、孤立孔隙的去除等。

其中所述N-S方程：

宏观连续模型基于连续假设，流体遵循质量守恒、动量守恒和能量守恒，其描述方程为Navier-Stokes方程(N-S方程)：

$\frac{\∂\mathrm{\ρ}}{\∂t}+\▿\·\left(\mathrm{\ρu}\right)=0---\left(1.10\right)$

$\frac{\∂\left(\mathrm{\ρu}\right)}{\∂t}+\▿\·\left(\mathrm{\ρuu}\right)=\▿\·\mathrm{\σ}---\left(1.11\right)$

$\frac{\∂\left(\mathrm{\ρe}\right)}{\∂t}+\▿\·\left(\mathrm{\ρue}\right)=\mathrm{\σ}:\▿u-\▿\·q---\left(1.12\right)$

式中，ρ为流体宏观密度；u为流体速度；T为温度；e为内能；σ为应力张量；q为热通量；上述方程组并不封闭，要求解压力、速度等参数必须补充本构关系、状态方程等关系。宏观连续模型是目前发展比较成熟，应用最为广泛的方法。

求解的N-S方程的具体形式如下：

$\begin{array}{c}\mathrm{\ρ}\frac{\∂u}{\∂t}+\mathrm{\ρ}(u\·\▿)u=\▿\·[-\mathrm{\ρI}+\mathrm{\μ}(\▿u+{(\▿u)}^T)]\\ \mathrm{\ρ}\▿\·u=0\end{array}---\left(1.13\right)$

其中，ρ为流体的密度；u为流体的速度；μ为流体的粘度；Ι为单位矩阵；

边界条件为：

$\begin{array}{c}p=p_0\\ [\mathrm{\μ}(\▿u)+{(\▿u)}^T]\·n=0\end{array}---\left(1.14\right)$

其中，p为压力，p₀为初始压力；

所述模拟步骤包括步骤①-⑤：

①导入剖分好的三维几何模型；

②整体参数设置，包括液体密度，原油粘度及压力；

③选择对应的岩样出口、入口及壁面，进行计算，得到岩样出口和入口的速度；

④对岩样出口和入口的速度进行积分，得到岩样出口和入口的流量；

⑤对模拟的结果进行绘图，得到岩样速度场和压力场分布。

根据上面得到的数据，可以绘制速度场、压力场分布及流线等，通过处理可以看出流体在数字岩心内的流动情况。

通过上述方法，可以对比不同类型岩石中边界层的影响状况。

边界层的厚度和粘度对微观孔隙中流体的渗流有很大影响。边界层的阻力作用在特征尺寸比较小的岩石中更为明显，主要表现为速度下降更快，压力减小更快；边界层对低渗介质的影响更为明显，在低渗透介质中研究边界层更为重要。

Claims (1)

1.一种低渗透储层原油边界层模拟方法，其特征在于，所述方法包括步骤如下：

1)CT扫描法构建数字岩心包括步骤(a)-(c)：

(a)通过CT扫描获取岩样的多方位CT投影图像数据；利用CT进行扫描，扫描后的单张CT投影图像数据被图像获取软件自动捕获并存储；将所述岩样沿其轴向旋转固定角度后再次对其CT扫描，扫描后的单张CT投影图像数据被图像获取软件自动捕获并存储，直至所述的岩样沿轴向累计旋转360°后结束采集CT投影图像数据；

(b)利用图像重建方法将所述CT投影图像数据转换为所述岩样的灰度图像；

所述的图像重建方法为变换法，包括不含反投影步骤的Fourier重建方法FRM和基于反投影的卷积或滤波方法：前者以直接Fourier变换方法DFM为代表，后者以卷积或滤波反投影C/FBP方法为代表；

DFM的基本思想是根据辐射源与成像目标之间的作用机理，寻找到目标函数f(x,y)的二维Fourier变换(FT)与投影函数的一维FT之间的关系，利用这种关系，尽管目标函数f(x,y)是未知的，但它的二维FT却能通过不同方位上投影函数的一维FT得到，一旦获得了f(x,y)的所有投影，也就获得了它的二维FT，f(x,y)则可由二维FT的反变换来求出；

基于反投影的卷积或滤波方法C/FBP方法源于Radon求反公式的直接实现；C/FBP公式利用FT和投影定理导出；设p(x',θ)的最高空间频率是B；

C/FBP的数学公式写成

$f(x,y)={\∫}_0^{\mathrm{\π}}q(x^{\′},\mathrm{\θ})\mathrm{d\θ}---\left(1.2\right)$

其中，

$q(x^{\′},\mathrm{\θ})={\∫}_{-\∞}^{\∞}\left|R\right|P(R,\mathrm{\θ})\exp \left(j2\mathrm{\π}{x}^{\′}R\right)\mathrm{dR}---(1.3)$

或

q(x',θ)＝c(x')*p(x',θ)       (1.4)

式(1.2)、(1.3)、(1.4)中，f(x,y)为图像重建的目标函数，即岩样对射线吸收系数的分布函数，x，y分别为岩样扫描断面上某点的直角坐标；P(R,θ)、p(x',θ)为f(x,y)沿θ方向上的投影函数，即CT探测器检测到的射线强度信号，x'为岩样扫描断面上某点(x，y)在旋转坐标系中的横坐标，R为该点到旋转中心，即坐标原点，的距离；c(x')为卷积函数，它由滤波器|R|(|R|≤B)的Fourier反变换求得；

式(1.3)、(1.4)中的q(x',θ)分别称为投影P(R,θ)、p(x',θ)的滤波投影和卷积投影；由式(1.2)、(1.3)实现的图像重建方法称为滤波反投影方法FBP，由式(1.2)、(1.4)实现的图像重建方法称为卷积反投影方法CBP；因此，如获得了图像的投影，则通过卷积或滤波反投影方法反向求解得到岩样断面对射线吸收系数的分布函数f(x,y)，将f(x,y)采用灰度图像的方式表达后就得到岩样断面的扫描图像；

(c)采用图像二值分割方法分离所述岩样的灰度图像中孔隙空间和岩样骨架，建立数字岩心；

在灰度图像中，孔隙空间和岩石骨架通过图像分割方法对其进行划分：图像分割在于分割阈值的选取；

阈值定义如下：假设F(i,j)表示灰度图像中图像像素灰度值，经分离后灰度值大于阈值T的像素表示目标对象，即孔隙，用1表示；其它的像素点表示图像背景，即骨架，用0表示，则由下式即可获得目标对象：

$g(i,j)=\left\{\begin{array}{cc}1;& F(i,j)\≥T\\ 0;& \mathrm{else}\end{array}\right.---\left(1.5\right)$

式中，T即为阈值，又称门限值；

基于岩心孔隙度的图像二值分割方法

设岩样的孔隙度为图像的最大、最小灰度值分别为I_Max、I_MIN，灰度值为i的像素数为p(i)，则满足式(1.6)的灰度值k^*即所求分割阈值：

当没有孔隙度时，利用最大类间距法分割灰度图像；

2)对上述已构建的数字岩心运用Delaunay三角剖分算法进行剖分，得到三维几何模型；

3)考虑边界层影响的条件下，运用有限元方法对流体在上述剖分后的三维描述几何模型中进行数值模拟，即对流体在孔道中微观流动的不可压缩方程N-S方程进行求解，得到流体流动的速度场和压力场分布；

其中所述N-S方程：

宏观连续模型基于连续假设，流体遵循质量守恒、动量守恒和能量守恒，其描述方程为Navier-Stokes方程(N-S方程)：

$\frac{\∂\mathrm{\ρ}}{\∂t}+\▿\·\left(\mathrm{\ρu}\right)=0---\left(1.10\right)$

$\frac{\∂\left(\mathrm{\ρu}\right)}{\∂t}+\▿\·\left(\mathrm{\ρuu}\right)=\▿\·\mathrm{\σ}---\left(1.11\right)$

$\frac{\∂\left(\mathrm{\ρe}\right)}{\∂t}+\▿\·\left(\mathrm{\ρue}\right)=\mathrm{\σ}:\▿u-\▿\·q---\left(1.12\right)$

式中，ρ为流体宏观密度；u为流体速度；T为温度；e为内能；σ为应力张量；q为热通量；

求解的N-S方程的具体形式如下：

$\begin{array}{c}\mathrm{\ρ}\frac{\∂u}{\∂t}+\mathrm{\ρ}(u\·\▿)u=\▿\·[-\mathrm{\ρI}+\mathrm{\μ}(\▿u+{(\▿u)}^T)]\\ \mathrm{\ρ}\▿\·u=0\end{array}---\left(1.13\right)$

其中，ρ为流体的密度；u为流体的速度；μ为流体的粘度；Ι为单位矩阵；

边界条件为：

$\begin{array}{c}p=p_0\\ [\mathrm{\μ}(\▿u)+{(\▿u)}^T]\·n=0\end{array}---\left(1.14\right)$

其中，p为压力，p₀为初始压力；

所述模拟步骤包括步骤①-⑤：

①导入剖分好的三维几何模型；

②整体参数设置，包括液体密度，原油粘度及压力；

③选择对应的岩样出口、入口及壁面，进行计算，得到岩样出口和入口的速度；

④对岩样出口和入口的速度进行积分，得到岩样出口和入口的流量；

⑤对模拟的结果进行绘图，得到岩样速度场和压力场分布。