具体实施方式

以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。

在低渗透储层内的油、水两相流动过程中，通常含有大量分散的油滴和气泡。当这些个体运动到储层中的狭小孔喉(例如孔径为2.5μm的孔隙)时，由于这些个体的直径大于孔喉直径，其流动受到阻力，迫使油滴或气泡等发生变形，故将油滴或气泡等在直径分布不均匀的储层中运动时受到的附加的毛细管阻力的现象称为贾敏效应，如图1所示。与水锁效应相比，微小孔喉中受到阻力被迫变形的是油滴或者气泡，此时润湿相是水，非润湿相是油滴或气泡。

通常情况下，产生贾敏效应需同时满足以下条件：首先，原油应该处于分散状态，不能以连续的油流形式运动；其次，地质储层的孔道直径应该比较细小。当分散的油滴运动到油层狭小孔道时，由于其直径大于孔道直径，就会受到阻力而致使其无法继续运动，或者油滴在变直径孔道中流动时会产生毛细管效应，这样就会导致贾敏损害(因贾敏效应导致储层的渗透率下降的情形)。

具体地，如图1所示，在储层左右两侧的驱替压力P的作用下，液滴从左向右运动。当油滴运动到孔隙1处时，因孔隙1的阻挡，对于油滴前端(2R₁＝λ，其中R₁为前端液面的曲率半径，λ为孔隙1的孔径)，其受到泡前压力P₁的作用；对于油滴后端(其直径其中R₂为后端液面的曲率半径，为孔隙的平均孔径)，其受到P₃的作用；以及对于油滴内部，其泡内压力为P₂，并且，P₁、P₂、P₃满足以下两个等式：

由此，因贾敏效应油滴受到的附加阻力为：

考虑液滴前端液滴后端

其中σ为液滴的表面张力。当驱替压力(即孔喉外压力差)至少达到P_c时，气泡或液滴才能通过孔喉，否则就会被堵住。

贾敏效应受储层的孔隙结构、岩性、物性、侵入流体等多种因素的影响。贾敏损害同岩石的孔隙介质的几何特性关系密切，不同的孔喉结构分布模式和复杂性会导致水润湿相在岩石中的展布模式显著变化，从而影响储层的渗透率。因此，本发明各实施例的核心是要建立储层中的孔隙内的水相饱和度变化的动力学模型(即储层中的水相运动方程及渗透率的分布方程)。具体地，基于储层中的孔隙内的流体的对流扩散关系(即质量平衡方程)等建立储层中的水相运动方程；基于孔隙的孔径分布特征及预设渗透率模型，建立渗透率的分布方程；然后根据所述水相运动方程及所述渗透率的分布方程，确定贾敏效应影响待诊断井周围的储层中的渗透率分布的时空演化控制唯象模型，由此就可诊断渗透率等储层损害特征参数的时空场分布。

需要说明的是，为了简单描述起见，在本发明的各个实施例中的随时空演化的物理量、化学量可省略变量例如可简写为K。

图1是本发明一实施例提供的贾敏效应损害储层的建模方法的流程图。所述建模方法可包括步骤S101-S104。

步骤S101，确定待诊断井的预设区域内的储层中的流体的达西表观速度。

其中，所述储层的渗透率低于预设渗透率；以及所述待诊断井可例如为采油井。

对于步骤S101，所述确定储层中的流体的速度可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述压力传导方程及达西公式，确定所述流体的达西表观速度。

具体地，压力是驱动固-液混合液从注水井的井筒持续侵入待诊断井周围的储层的动力，由此可建立如公式(1)的所述流体进入储层的压力传导方程：

其中，所述压力传导方程的初始条件为以及边界条件为(也就是说，在注水井的井壁处的压力为P_w)。

根据公式(1)及达西公式(如下式(2))可确定所述流体的达西表观速度，

其中，为所述流体的压力；φ₀为所述储层的孔隙度的初始值；μ为流体粘度；c_t为流体-岩石综合压缩系数以及为所述储层的渗透率。

步骤S102，根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程。

在储层条件下，所述储层中的孔隙内的不同位置的含水量满足质量守恒方程。其中，外来水相在所述储层内的运动主要由对流与扩散两个过程决定。具体地，对于步骤S102，所述建立所述储层的水相运动方程可包括：根据所述流体的达西表观速度u与所述水分子的扩散系数D_w，建立下式表示的所述质量平衡方程：

其中，φ₀为所述储层的孔隙度的初始值；为所述储层内的孔隙被水相占据的绝对孔隙度；以及为所述储层内的任意点的空间位置(例如，以所述待诊断井的中心为原点)。

根据所述质量平衡方程及所述储层的水相饱和度的时空分布函数建立下式(3)表示的所述水相运动方程：

所述水相运动方程的初始条件为以及所述水相运动方程的边界条件为(也就是说，在注水井的井壁处的储层孔隙完全被水充满，即孔隙中的水相饱和度为1)。其中，φ₀为所述储层的孔隙度的初始值；r_w为所述待诊断井的井筒半径；以及S_wc为所述储层中的束缚水饱和度。

步骤S103，根据所述储层的孔隙的孔径分布特征及所述储层的预设渗透率模型，建立所述储层的渗透率的分布方程。

对于步骤S103，如图2所示，所述建立所述储层的渗透率的分布方程可包括步骤S201-S202。

步骤S201，根据所述储层的孔隙的孔径分布特征，确定孔径为λ的孔隙的体积密度函数与所述储层的水相饱和度的孔径分布方程。

储层岩石的孔隙结构具有强烈的不规则性，受到多种因素的控制，使得经典几何学为基础的孔隙结构模型与真实的孔隙相差太远。为了定量描述储层的孔隙结构，采用分形理论来研究孔隙结构中两相流的贾敏效应，并忽略水相在孔隙通道内表面的润湿薄膜影响。根据分形的几何原理，若储层的孔径分布具有分形特征，则储层中孔径大于λ的孔隙的数量N(>λ)与λ有如下幂函数关系：

其中，D为孔隙的分形维数(2<D<3)。

在所述储层的孔隙的孔径分布特征为所述储层中的孔径大于λ的孔隙的数量N(＞λ)满足上式的情况下，如图3所示，所述步骤S201中的确定孔径为λ的孔隙的体积密度函数可包括步骤S301-S302。

步骤S301，根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，确定所述储层中的孔隙的总体积为

具体地，根据上式可得到储层中的孔隙的总数目N(＞λ_min)及孔径小于λ的孔隙的数目：

根据公式(4)与(5)可得到下式(6)

上式(6)中的f(λ)为储层的孔径分布密度函数，则孔径小于λ_pc的孔隙的数量N(＜λ_pc)与λ之间的关系为下式(7)所表示的幂函数关系：

上述相关公式中的λ、λ_min、λ_max与λ_pc分别为所述孔隙的孔径、最小孔径、最大孔径(λ_min与λ_max可以由平均孔径与孔径分布标准差得出，一般而言有)与发生贾敏效应的最大孔径(即非水相受到所述最小毛细管阻力处的孔隙的特定孔径)；a为比例常数。

接着，由公式(7)可得到储层的孔径分布密度函数f(λ)满足下式(8)：

式中a'＝-Da为比例常数。

由上式(8)表示的所述孔径分布密度函数可以得到储层中的孔隙的总体积分形表达式：

其中α为与孔隙的形状有关的常数(孔隙的形状为立方体时α＝1，或者孔隙的形状为球体时α＝π/6)，通过积分可得：

类似地，所述储层中的孔径小于λ的孔隙的体积

步骤S302，根据所述储层中的孔隙的总体积Φ_max及所述储层中的孔径小于λ的孔隙的体积确定所述孔径为λ的孔隙的体积密度函数为

其中，D为所述孔隙的分形维数；以及λ、λ_min与λ_max分别为所述孔隙的孔径、最小孔径与最大孔径；以及A＝αa'/(3-D)(常数)。

步骤S202，根据所述预设渗透率模型、所述孔径为λ的孔隙的体积密度函数及所述水相饱和度的孔径分布方程，建立所述储层的渗透率的分布方程。

如图4所示，步骤S202中的确定所述储层的水相饱和度的孔径分布方程可包括步骤S401-S402。

步骤S401，根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，确定被非水相占据的孔隙的体积为

其中，λ_pc为非水相受到所述最小毛细管阻力处的孔隙的特定孔径(即发生贾敏效应时的最大孔喉直径)。

在储层为水湿(即亲水性)的情况下，向所述储层的两侧施加一个驱替压力P(从零逐渐增加到P)，根据贾敏效应产生条件，油滴或气泡可顺利通过孔径处于λ_max到λ_pc范围之间的孔隙(即连通孔径由最大λ_max到λ_pc逐次打开)。当压力P未到P_c时，λ大于λ_pc的孔喉完全被水占据，而λ小于λ_pc的孔喉还被油滴或气泡堵塞。

由于根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，可确定上式(8)确定的孔径分布密度函数f(λ)，则进而可根据f(λ)得到压力P未到P_c时的非水相占据的孔隙体积Φ_nw(λ)，

步骤S402，根据所述储层中的孔隙的总体积Φ_max与所述被非水相占据的孔隙的体积Φ_nw(λ)，确定下式表示的所述水相饱和度的孔径分布方程

其中，D为所述孔隙的分形维数；λ、λ_min与λ_max分别为所述孔隙的孔径、最小孔径与最大孔径；以及A＝αa'/(3-D)。

具体地，根据公式(10)与公式(12)可确定水相占据的孔隙体积为：然后再结合公式(10)可得到下式：

由于λ_min＜＜λ_max，故上式可写为下式(13)，即可确定所述水相饱和度的孔径分布方程(13)，

根据流体力学理论，在压差ΔP作用下，粘度为μ的流体流过毛细管束总流量可以由Hagen-poiseuille方程描述，

上式中，λ_i为孔隙通道i的直径(即孔隙i的孔径)，V_i为孔隙通道i的体积，V_i＝Φ₀ML_mξ_i，Φ₀为储层的孔隙度，ξ_i为孔隙通道i所占孔隙体积的分数，L_c为孔隙通道i的弯曲长度，L_m为孔隙通道i的直线长度，M为平均孔喉截面积为孔隙通道i的迂曲度(同时可由经验公式计算)，将上述各量代入上式得到，

结合达西定律，可得到毛管束模型的渗透率可表示为在本发明实施例中，可将储层的渗透通道看做是多个毛管束累积起来的。由于孔隙的尺寸分布具有连续性，因此可以将所述毛管束模型的渗透率的表达式用积分形式写出：

在所述储层的预设渗透率模型满足的情况下，所述建立所述储层的渗透率的分布方程可包括：根据所述储层的预设渗透率模型所述孔径为λ的孔隙的体积密度函数dξ及所述水相饱和度的孔径分布方程，建立下式表示的所述渗透率的分布方程，

其中，所述压力传导方程的初始条件为K(t＝0)＝1。K(S_w)为渗透率关于含水饱和度的函数(如图5所示，该函数可简称为无因次渗透率)；φ₀为储层的初始孔隙度，S_w为含水饱和度；D为所述孔隙的分形维数，其可由计算(D越大说明小孔喉所占比例越大，常规储层λ_max＞20μm，D<2.6，低渗储层λ_max＜0.2μm，D>2.8)，τ为迂曲度，λ_max为最大孔喉直径，λ_min为最小孔喉直径)。

具体地，首先，逐渐增加驱替压力直至P，此时直径大于λ的孔喉被水水占据，被水相占据的绝对孔隙度为φ_w，则将公式(11)代入公式(14)可得出所述储层的渗透率的孔径分布函数：然后，根据式(13)中的λ与含水饱和度S_w(λ)之间的关系式将上述渗透率的孔径分布函数K(λ)中的变量λ替换为S_w，并结合φ_w＝φ₀S_w，可得到公式(15)(即得到无因次渗透率的表达式，该无因次渗透率为含水饱和度的函数)。

步骤S104，根据所述渗透率的分布方程及所述水相运动方程，确定贾敏效应损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由贾敏效应引起的储层损害特征的四维时空演化过程。

具体地，根据公式(3)可求解到所述储层的水相饱和度的时空分布函数(如图6所示)，并将代入公式(15)可得到所述储层的渗透率在四维时空上的分布形式，即贾敏效应损害储层的时空演化模拟方程。

综上所述，本发明创造性地确定待诊断井的预设区域内的储层中的流体的达西表观速度，其中所述储层的渗透率低于预设渗透率；根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程；建立所述储层的渗透率的分布方程；以及根据所述渗透率的分布方程及所述水相运动方程，确定贾敏效应损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由贾敏效应引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。

相应地，本发明另一实施例还提供一种确定储层损害程度的方法，所述方法可包括：基于所述的贾敏效应损害储层的建模方法建立的时空演化模拟方程，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。

对于上述公式(15)所示的贾敏效应损害储层的时空演化模拟方程的求取，需要根据公式(3)计算得到对于公式(3)而言，在一维情形下，该类方程可以整理为以下的一般形式：

其中，a_a,b_b,c_c可为常数(如扩散系数)，也可为函数(如所述流体的速度)；f可为压力、水相饱和度、物质浓度(例如体积分数)、应力等。对时间采用向后差分，空间采用中心差分。则上述方程可以有如下差分格式：

其中，i＝1,2,3…N_i，n＝1,2,3…，t＝nΔt，N_i为离散空间点个数。

求解区间为x∈(0,x_max)，Δx、Δt为空间、时间步长。同时，考虑初始条件f_i ⁿ|_n＝0＝f_i ⁰,i＝1,2,3…,N_i和边界条件(f_i ⁿ|_i＝1＝f₀,n＝1,2,3...(井壁处)以及n＝1,2,3...)(构造了一个虚拟网格i+1，预设范围的边界处或距井壁数米处)。

首先，对于i＝2,3,...,N_i-1整理上述差分格式有：

其中，A1_i，A2_i，A3_i分别为，

同时，据公式(3)可确定a_i、b_i与c_i。并将所确定的a_i、b_i与c_i代入公式(19)可得到迭代关系式(18)的具体表现形式，由于该迭代关系式(18)的具体表现形式复杂，故在此不对其进行列出。然后，利用初始条件和边界条件进行迭代计算就可得到场f的值。

接着，对说明边界条件的差分求解过程进行说明。

上述迭代关系式(18)适用于非边界网格。而对于i＝1(井壁处)而言，因为采用的是点中心网格，且其为狄利克雷(Dirichlet)边界条件，故直接可得到以下关系式：

f₁ ⁿ＝f₀(常数)，i＝1 (20)

对于i＝N(预设范围的边界处距井壁数米处)而言，其为诺伊曼或第二类(Neumann)边界条件，增加一个虚拟网格i＝N_i+1，由n＝1,2,3...知将其代入式(18)可知：

根据上述过程可求解出场函数f的时空变化情况。由于上述数值模型是针对待诊断井(注水井)的井筒附近储层建立的，在求解某物理量f在井周的分布时，需要采用柱坐标系。由此，式需要变换为这种形式不利于等距差分，可以引入坐标变换：r＝r_we^x′，其中，r_w为井筒半径，x′为一个无量纲的空间坐标。将这个变换代入一般方程中，可以得到关于x′的方程：

如果将和作为新的方程系数，则上式和相比，本质上是一样的。因此，便可以在x′坐标进行等距差分并沿用前述的迭代格式。计算完f的值后，再将空间坐标从x′映射回r即可得到f(r,t)。

在通过上述方法计算得到所述储层的水相饱和度之后，再根据上述公式(15)可计算得到所述储层的渗透率由此通过上述贾敏效应损害储层的建模方法建立的时空演化模拟方程综合考虑了贾敏损害时多种物理化学因素对储层损害的影响，由此通过本实施例求解得到的储层的渗透率非常精确。

在得到所述储层的渗透率的基础上，表征待诊断井的预设区域内的储层的损害程度的特征参数。

在一实施例中，所述特征参数可为所述储层的渗透率损害率。

相应地，所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：

基于所述时空演化模拟方程，确定所述储层的渗透率以及基于所述储层的渗透率及公式(23)，确定所述储层的渗透率损害率

其中，为的最大值。

在另一实施例中，所述特征参数可为所述储层的表皮系数。

所述确定表征待诊断井的预设区域内的储层的损害程度的特征参数可包括：基于所述时空演化模拟方程，确定所述储层的渗透率以及基于所述储层的渗透率及公式(24)，确定所述储层的表皮系数

其中，为所述储层的渗透率的初始值，以及r_w为所述待诊断井的井筒半径，以及r_sw为所述储层的损害半径。

通过上述各个实施例得到的特征参数(例如所述储层的渗透率损害率I(S_w)是时空演化4D定量模拟的结果(分别如图7所示)。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。

图9是本发明一实施例提供的贾敏效应损害储层的建模系统的结构图。如图9所示，所述建模系统包括：速度确定装置10，用于确定待诊断井的预设区域内的储层中的流体的达西表观速度，其中所述储层的渗透率低于预设渗透率；第一建立装置20，用于根据所述流体的达西表观速度与所述流体中的水分子的扩散系数，建立所述储层的水相运动方程；第二建立装置30，用于根据所述储层的孔隙的孔径分布特征及所述储层的预设渗透率模型，建立所述储层的渗透率的分布方程；以及模拟方程确定装置40，用于根据所述渗透率的分布方程及所述水相运动方程，确定贾敏效应损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由贾敏效应引起的储层损害特征的四维时空演化过程。

可选的，所述速度确定装置10包括：压力传导方程建立模块(未示出)，用于所述流体进入所述储层的压力传导方程；以及速度确定模块(未示出)，用于根据所述压力传导方程及达西公式，确定所述流体的达西表观速度。

可选的，所述第一建立装置20包括：质量平衡方程建立模块(未示出)，用于根据所述流体的达西表观速度u与所述水分子的扩散系数D_w，建立下式表示的所述质量平衡方程：其中，φ₀为所述储层的孔隙度的初始值；为所述储层内的孔隙被水相占据的绝对孔隙度；以及为所述储层内的任意点的空间位置；以及水相运动方程建立模块(未示出)，用于根据所述质量平衡方程及所述储层的水相饱和度的时空分布函数建立下式表示的所述水相运动方程：

可选的，所述第二建立装置30包括：第一分布方程建立模块(未示出)，用于根据所述储层的孔隙的孔径分布特征，确定孔径为λ的孔隙的体积密度函数与所述储层的水相饱和度的孔径分布方程；以及第二分布方程建立模块(未示出)，用于根据所述预设渗透率模型、所述孔径为λ的孔隙的体积密度函数及所述水相饱和度的孔径分布方程，建立所述储层的渗透率的分布方程。

可选的，在所述储层的孔隙的孔径分布特征为所述储层中的孔径大于λ的孔隙的数量N(＞λ)满足的情况下，所述第一分布方程建立模块(未示出)包括：第一体积确定单元，用于根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，确定所述储层中的孔隙的总体积为以及体积密度函数确定单元，用于根据所述储层中的孔隙的总体积Φ_max及所述储层中的孔径小于λ的孔隙的体积确定所述孔径λ的孔隙的体积密度函数为其中，D为所述孔隙的分形维数；以及λ、λ_min与λ_max分别为所述孔隙的孔径、最小孔径与最大孔径；以及A＝αa'/(3-D)。

可选的，所述第一分布方程建立模块(未示出)包括：第二分布方程建立模块(未示出)，用于根据所述储层中的孔径大于λ的孔隙的数量N(＞λ)，确定被非水相占据的孔隙的体积为其中λ_pc为非水相受到所述最小毛细管阻力处的孔隙的特定孔径；以及水相饱和度分布方程确定单元，用于根据所述储层中的孔隙的总体积Φ_max与所述被非水相占据的孔隙的体积Φ_nw(λ)，确定下式表示的所述水相饱和度的孔径分布方程，其中，D为所述孔隙的分形维数；λ、λ_min与λ_max分别为所述孔隙的孔径、最小孔径与最大孔径；以及A＝αa'/(3-D)。

所述贾敏效应损害储层的建模系统与上述贾敏效应损害储层的建模方法相对于现有技术所具有的优势相同，在此不再赘述。

图8是本发明一实施例提供的确定储层损害程度的系统的结构图。如图8所示，所述系统可包括：接收装置50，用于接收基于所述的贾敏效应损害储层的建模系统建立的时空演化模拟方程；以及特征参数确定装置60，用于基于所述时空演化模拟方程，确定表征待诊断井的预设区域内的储层的损害程度的特征参数。

可选的，所述特征参数为所述储层的渗透率。相应地，所述特征参数确定装置60包括：渗透率计算模块(未示出)，用于基于所述时空演化模拟方程，确定所述储层的渗透率

可选的，所述特征参数为所述储层的表皮系数。相应地，所述特征参数确定装置60包括：渗透率计算模块(未示出)，用于基于所述时空演化模拟方程，确定所述储层的渗透率以及表皮系数计算模块(未示出)，用于基于所述储层的渗透率及公式确定所述储层的表皮系数

所述确定贾敏效应程度的系统与上述确定贾敏效应程度的方法相对于现有技术所具有的优势相同，在此不再赘述。

相应地，本发明一实施例还提供一种机器可读存储介质，该机器可读存储介质上存储有指令，该指令用于使得机器执行所述的贾敏效应损害储层的建模方法和/或所述的确定储层损害程度的方法。

所述机器可读存储介质包括但不限于相变内存(相变随机存取存储器的简称，Phase Change Random Access Memory，PRAM，亦称为RCM/PCRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体(Flash Memory)或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁盘存储或其他磁性存储设备等各种可以存储程序代码的介质。

上述步骤S101-S104、步骤S201-S202、步骤S301-S302及步骤S401-S402均可通过计算机来执行。并且，步骤S101-S104所涉及的各种物理化学量的处理过程实现了对贾敏效应损害储层的时空演化场的模拟，步骤S201-S202、步骤S301-S302及步骤S401-S402所涉及的各种物理化学量的处理过程实现了对贾敏效应损害储层的时空演化场的具体模拟。

以上结合附图详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。