一种多孔介质中预交联凝胶悬浮液微观流动的模拟方法与流程

本发明涉及一种多孔介质中预交联凝胶悬浮液微观流动的数值模拟方法，属于油气田开发提高原油采收率、油藏数值模拟和计算流体动力学的交叉领域。

背景技术：

预交联凝胶颗粒驱是高含水油田进一步提高原油采收率的新兴深部调驱技术之一。然而，预交联凝胶悬浮液在多孔介质中的流动十分复杂，不仅具备刚性颗粒悬浮液的基本渗流特征，同时还具有弹性变形等特性。目前，考虑这些特性的油藏数值模拟技术多基于经典渗滤理论或尺寸排斥理论，属于宏观数值模拟的范畴，受限于连续介质假设，无法模拟预交联凝胶颗粒在多孔介质中“运移、封堵、变形、再运移”的真实过程。

为此，本发明提出一种多孔介质中预交联凝胶悬浮液微观流动的模拟方法，兼顾了颗粒变形的精细化处理、颗粒之间的接触作用模拟以及颗粒流体间的高效耦合，为深入认识预交联凝胶悬浮液的复杂渗流机理提供了一种研究手段。

技术实现要素：

为弥补现有技术的不足，本发明将燃烧算法和转折点算法结合，提出一种可以识别孔喉末端的孔隙喉道识别方法，其识别过程能准确识别末端并且显著降低识别后的喉道长度与实际长度的相对误差。

本发明技术方案具体步骤如下：

(1)对图像进行二值化处理，区分岩石颗粒与孔喉空间；

(2)以岩石颗粒为前景像素对孔喉空间进行欧式长度变换，长度变换计算公式为：

其中p为孔喉空间内一点，q为p点周围任意的岩石颗粒点。

(3)由基于长度变换的转折点算法提取图像中线，此时的中线未考虑末端；

(4)基于转折点中线，进行不考虑末端时的孔喉识别；

(5)利用燃烧算法提取中线，查找燃烧中线尾端端点，确定末端所在区域；

(6)借助燃烧中线尾端端点，利用转折点算法提取末端中线；

(7)将步骤(4)中的转折点中线与步骤(6)中的末端转折点中线合并，对末端区域进行孔喉识别，然后将末端识别结果与步骤(4)中未考虑末端区域的孔喉识别结果合并，得到考虑末端的孔喉识别结果。

所述步骤(4)中，孔喉识别的步骤如下：

①删除位于孔隙内部的转折点中线路径，将剩余中线路径所对应像素按顺序编号，并以每个像素的长度变换结果作为此像素对应的孔隙半径；

②以像素编号为横坐标，对应的孔隙半径为纵坐标，做出每条路径的沿程曲线；

③将沿程曲线看作信号，选择离散dmey波为小波基函数，对信号进行二尺度分解，选择合理阈值截取小波系数，然后重构信号完成小波消噪(参考文献：孙延奎.小波分析及其应用[m].dynomediainc.,2005.)，阈值选择公式为：

其中σ为噪声标准方差，n是信号长度；

④在去噪后的沿程曲线上查找最靠近两端的两个局部极小值，将其对应的像素作为孔隙和喉道的分界点即喉道的瓶颈，在瓶颈对应的位置画出孔喉识别线；

⑤孔喉识别线和岩石边界将非基质空间分成了若干封闭区域，将包含中线结点的区域作为孔隙，其余部分作为喉道。

所述步骤(6)中，提取末端中线步骤如下：

①从步骤(4)识别后的孔隙或喉道中，提取包含燃烧中线端点像素的喉道或孔隙区域即为末端区域；

②以燃烧中线的末端像素为原点，沿正交方向将提取的末端区域切分为四个子图，应用转折点算法分别提取四个子图的中线；

③将四个子图的中线合并在一起并清理分支，所得结果即为转折点中线在末端处的分支。

本发明具有以下有益效果及优点：

(1)以转折点算法为主要算法提取孔隙空间的中线，使中线的整体居中效果得到改进；

(2)利用小波消噪后的沿程曲线查找瓶颈，识别后的喉道长度与实际喉道长度相对误差显著降低；

(3)由燃烧算法获取的端点引导转折点算法提取末端部分的中线，可以将末端部分单独识别出来。

附图说明

图1为多孔介质中预交联凝胶悬浮液微观流动的模拟方法流程图。

图2为多孔介质模拟区域示意图。

图3为多孔介质中数值网格节点的类型标记示意图。

图4为多孔介质中预交联凝胶悬浮液微观流动的模拟结果图。

具体实施方式

图1给出了本发明的详细实施流程，下面将结合其他附图和具体实施例作进一步的说明，以便更好的理解本发明，但不限定本发明的范围。

(1)基于二维岩心ct切片图像，建立多孔介质模拟区域，如图2所示。该模拟区域多孔介质部分长度为260μm，为了模拟预交联凝胶颗粒进入和离开多孔介质的过程，分别在多孔介质入口端和出口端设置长度为40μm的管道，因此模拟区域总长度为340μm，宽度为100μm。设定网格步长为0.2μm，将该模拟区域进行数值网格离散，则共有1700×500＝850000个计算节点。

(2)设定预交联凝胶颗粒的初始位置在多孔介质入口端，原始形状为圆形，直径为16μm。将颗粒的真实曲线边界离散生成18个虚拟边界标记点，并由其连线组成的18边形近似表征预交联凝胶颗粒。

(3)标记数值网格节点类型，若某数值网格节点及其相邻8个节点均被预交联凝胶颗粒覆盖则为颗粒内部节点(图3中实心圆点)，若某数值网格节点被预交联凝胶颗粒覆盖但存在未被覆盖的相邻节点则为颗粒边界节点(图3中实心三角形)，若某数值网格节点及其相邻8个节点均未被预交联凝胶颗粒覆盖则为流体内部节点(图3中空心圆点)，若某数值网格节点未被预交联凝胶颗粒覆盖但存在被覆盖的相邻节点则为流体边界节点(图3中空心三角形)。图3中数值网格节点a为流体内部节点，b为颗粒内部节点，s为流体边界节点。

(4)输入网格步长δx＝0.2μm，时间步长δt＝0.0025s，松弛时间τ＝0.65，根据预交联凝胶颗粒边界上数值网格节点的动量变化计算流体对各虚拟边界标记点的作用力，其计算公式为：

其中，δx表示网格步长，δt表示时间步长，vt表示颗粒边界上以某数值网格节点为中心边长为δx的控制单元的总体积，vs表示上述控制单元中被颗粒覆盖部分的体积；τ表示松弛时间；cj为附加碰撞项；lj表示j方向上的离散速度；p表示各虚拟边界标记点控制范围内属于颗粒边界类型的数值网格节点个数。

(5)计算虚拟边界标记点之间的接触重叠量δ，输入接触刚度k＝10⁴kg/s²，阻尼系数λ＝0.3，计算各虚拟边界标记点受到的接触作用力，其计算公式为：

其中，k表示接触刚度，λ表示阻尼系数，δ表示各虚拟边界标记点之间的接触重叠量。

(6)输入预交联凝胶颗粒的弯曲恢复刚度γb＝10²kg/s²，拉伸恢复刚度γs＝10kg/s²，根据预交联凝胶颗粒的当前形状和原始形状，计算各虚拟边界标记点受到的弯曲恢复力和拉伸恢复力，其计算公式分别为：

fb＝γb(α-α⁰)

其中，α⁰和α分别表示预交联凝胶颗粒变形前后相邻三个虚拟边界标记点之间的向量夹角，和分别表示预交联凝胶颗粒变形前后相邻两个虚拟边界标记点之间的距离，γb表示弯曲恢复刚度，γs表示拉伸恢复刚度。

(7)基于牛顿第二定律，计算各虚拟边界标记点的速度u和位移x，将各虚拟边界标记点相连即得到预交联凝胶颗粒的整体形状和位置；

x＝xold+uδt

其中，uold和u分别表示迭代前后虚拟边界标记点的速度大小，m表示预交联凝胶颗粒的质量，n表示虚拟边界标记点的个数，xold和x分别表示迭代前后虚拟边界标记点的位移大小。

(8)重复步骤(4)至步骤(7)，当前后两次模拟结果相对误差小于0.01或者迭代次数大于100000时，结束模拟过程，输出模拟结果，如图4所示。

由上述实施例可以看出，本发明公开的一种多孔介质中预交联凝胶悬浮液微观流动的模拟方法能够直观反映预交联凝胶颗粒在多孔介质中的真实变形和运移过程，并能够定量表征预交联凝胶悬浮液在多孔介质中的微观渗流机制。因此，本发明为预交联凝胶悬浮液等含弹性颗粒悬浮液的微观渗流机制研究提供了一种有效的数值模拟手段。