具体实施方式

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

储层内部的孔隙结构和固体基质(即岩石基质或骨架颗粒)的分布及其复杂。微观的一个固体基质(即团簇)形成一条渗流通道，多个固体基质(即团簇)之间夹杂着多个孔隙(相当于水力管束)，而所述多个固体基质与所述多个孔隙形成如图1A所示的毛管束(也称为毛细管束)。这些孔隙具有特定的迂曲度，使固体基质与孔隙的实际长度L_c0比由毛管束形成的多孔介质的长度L₀(即毛管束的表观长度)更长，如图1A所示。

并且，一方面因储层的上覆岩层的作用而使储层受到向下的压力，另一方面因储层内部的流体的支撑作用而使储层受到向上的压力，由此，储层在这两个压力(即有效应力)作用下，毛管束收缩(如图1B所示)、毛管束延展(如图1C所示)，进而会影响液体在储层中的渗透率(即简称为所述储层的渗透率)。应力敏感性受储层的孔隙结构、岩性、物性、侵入流体等多种因素的影响。因此，本发明各实施例的核心是要建立有效应力作用下的储层中的流体的流量与渗透率的分布方程。具体地，根据储层的孔隙的孔径分布特征、有效应力下的储层的各个毛管束的直径与长度及流体流量公式，确定所述储层中的流体的流量，再结合渗透率模型，就可诊断渗透率等储层损害特征参数的时空场分布。

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

图2是本发明一实施例提供的应力敏感性损害储层的建模方法的流程图。如图2所示，所述建模方法包括步骤S201-S203。

步骤S201，确定待诊断井的预设区域内的储层受到的有效应力。

对于步骤S201，所述确定待诊断井的预设区域内的储层受到的有效应力可包括：建立所述流体进入所述储层的压力传导方程；以及根据所述储层的压力传导方程与所述储层的上覆岩层的压力，确定所述储层受到的有效应力。

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

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

然后根据公式(1)可得到所述流体对所述储层产生的压力(即孔隙压力)，再根据与上覆岩层的压力确定所述储层受到的有效应力

步骤S202，根据所述储层的孔隙的孔径分布特征、所述有效应力下的所述储层的各个毛管束的直径与长度、及流体流量公式，确定所述储层中的流体的流量。

其中，所述毛管束由多个固体基质及所述多个固体基质之间的孔隙组成。

对于步骤S202，如图3所示，所述确定所述储层中的流体的流量可包括步骤S301-步骤S302。

在执行步骤S301之前，所述确定所述储层中的流体的流量还可包括：根据所述储层的弹性模量，确定所述有效应力下的所述储层的各个毛管束的直径；以及根据所述储层的弹性模量与泊松比，确定所述有效应力下的所述储层的各个毛管束的长度。

下面，分别介绍关于所述有效应力下的所述储层的各个毛管束的直径与长度的确定过程。

待诊断井的横波速度是确定储层的泊松比、弹性模量等重要数据的关键特征。而在现场测井时，一般不会针对待诊断井进行横波测井，而是通过含泥量和关于含泥量的线性公式确定泊松比。其中，所述线性公式通过特定岩石样本拟合而成，故针对不同的储层，通过所述线性公式确定的泊松比的精确性非常差。在本实施例中，根据已获取的待诊断井周围的多个邻井的横波速度、纵波速度及对待诊断井进行现场测井所获取的纵波速度，可间接计算得到所述诊断井的横波速度。

在执行所述确定所述有效应力下的所述储层的各个毛管束的长度的步骤之前，所述确定所述储层中的流体的流量还可包括：根据位于所述预设区域内的多个特定相邻井中的每一者的横波速度与纵波速度，确定所述储层的横波速度；以及根据所述储层的横波速度，确定所述储层的所述弹性模量与所述泊松比。

其中，所述特定相邻井可为已通过测井方式获取横、纵波速度的井。

具体地，所述待诊断井的横波速度υ_t可由纵波速度υ_l及关系式υ_t＝aυ_l-b 确定。其中，所述关系式中的系数a、b可基于统计的邻井的横波速度与纵波速度计算得到。例如，通过最小二乘法来线性回归出所述系数a、b：

其中，n为通过测井得到其横、波速度的所述预设区域内的邻井的数目；υ_ti为第i个邻井的横波速度；υ_li为第i个邻井的纵波速度；为n个邻井的横波速度的平均值；为n个邻井的纵波速度的平均值。

与上述现有的方式获取的泊松比相比，由于测井数据非常精确，故本实施例通过邻井的横、纵波速度等测井数据获取的待诊断井的横波速度的结果非常精确，从而后续获取的泊松比等数据也非常准确。接着下面将介绍如何通过待诊断井的横波速度获取储层的泊松比与弹性模量。

具体地，可通过下式(3)与所述待诊断井的横波速度，得到储层的泊松比ν与弹性模量E：

其中，ρ为所述储层的体积密度；Δt_t为横波时差(可由单发射双接收声速测量装置测量得到，即垂直方向上相距L的两个接收器接收到声波所形成的滑行波横波的时间差)；以及Δt为声波时差(可由单发射双接收声速测量装置测量得到，即垂直方向上相距L的两个接收器接收到声波所形成的滑行波的时间差)。

下面基于所获取的弹性模量与泊松比，确定一定应力作用下的储层的任一毛管束的直径与长度。

第一，所述确定所述有效应力下的所述储层的各个毛管束的直径可包括：根据所述储层的弹性模量E及下式，确定所述有效应力下的所述储层的特定毛管束的直径λ_p，

其中，F为所述特定毛管束的形状因子；以及λ_c0为无有效应力下的所述储层的特定毛管束的初始直径。

第二，所述确定所述有效应力下的所述储层的各个毛管束的长度可包括：根据无有效应力下的所述储层的特定毛管束的表观长度、无有效应力下的所述储层的特定毛管束的初始直径及分形尺度法则，确定无有效应力下的所述储层的特定毛管束的初始长度L_c0；以及根据所述储层的弹性模量E、所述储层的泊松比ν、无有效应力下的所述储层的特定毛管束的初始长度L_c0及下式，确定所述有效应力下的所述储层的特定毛管束的长度L_p，

其中，所述确定无有效应力下的所述储层的特定毛管束的初始长度L_c0可包括：根据无有效应力下的所述储层的特定毛管束的表观长度L₀、无有效应力下的所述储层的特定毛管束的初始直径λ_c0及分形尺度法则确定所述特定毛管束的初始长度L_c0，其中，D_cT为所述特定毛管束的迂回分形维度。

步骤S301，根据所述有效应力下的所述储层的各个毛管束的直径与长度及流体流量比率公式，确定所述流体的流量比率。

根据上述得到的所述有效应力下的所述储层的特定毛管束的直径λ_p与长度L_p、及流体流量比率公式确定所述流体的流量比率

其中，G为所述特定毛管束的横截面的几何因子；Δp为所述储层的横截面两端的驱替压力；以及F为所述特定毛管束的形状因子。

步骤S302，根据所述储层的孔隙的孔径分布特征、所述储层中的流体的流量比率及所述流体流量公式，确定所述流体的流量。

对于步骤S302，所述确定所述流体的流量包括：根据所述储层的孔隙的孔径分布特征所述储层中的流体的流量比率及所述流体流量公式确定所述流体的流量

其中，λ_c0为所述储层的特定毛管束的初始直径；λ_c0max为所述储层的最大毛管束与最小毛管束的初始直径；D_cf为所述孔隙的分形维度；以及为所述有效应力。

具体地，基于岩石分形几何理论，直径为λ_c0的固体基质的孔隙的数密度N_c(即单位体积内某一直径固体基质的数量累积分布，或孔隙的孔径分布特征)满足关系：对于二维(储层)模型，0＜D_cf＜2；对于三维(储层)模型，0＜D_cf＜3。

根据上述储层的孔隙的孔径分布特征，可得到然后将dN_c及公式(4)代入所述流体流量公式，可得到下式(5)表示的所述流体的流量

步骤S203，根据所述储层的渗透率模型及所述储层中的流体的流量，确定应力敏感性损害储层的时空演化模拟方程。

其中，所述时空演化模拟方程用于模拟由应力敏感性引起的储层损害特征的四维时空演化过程。

对于步骤S203，所述确定应力敏感性损害储层的时空演化模拟方程可包括：根据所述储层的渗透率模型及所述储层中的流体的流量确定应力敏感性损害储层的时空演化模拟方程。

其中，μ为所述流体的粘度；L₀为无有效应力下的所述储层的特定毛管束的表观长度；A为所述储层的横截面的面积；以及Δp为所述储层的横截面两端的驱替压力。具体地，然后代入上述dN_c可得到A。

具体地，根据所述储层的渗透率模型及上式(5) 表示的流量可得到下式(6)表示的应力敏感性损害储层的时空演化模拟方程，

综上所述，本发明创造性地确定待诊断井的预设区域内的储层受到的有效应力；根据所述储层的孔隙的孔径分布特征、所述有效应力下的所述储层的各个毛管束的直径与长度、及流体流量公式，确定所述储层中的流体的流量；以及根据所述储层的渗透率模型及所述储层中的流体的流量，确定应力敏感性损害储层的时空演化模拟方程。由此，通过所确定的时空演化模拟方程可定量模拟由应力敏感性引起的储层损害特征的四维时空演化过程，从而对未发生储层损害的井进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义，以及对已损害井优化设计解堵措施、提高或恢复油井产量和水井注水能力，以及提高油藏数值模拟精度都具有十分重大意义。

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

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

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

其中，i＝1,2,3...N_i，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...(井壁处)以及)(构造了一个虚拟网格i+1，预设范围的边界处或距井壁数米处)。

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

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

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

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

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

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

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

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

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

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

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

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

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

其中，为的最大值。

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

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

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

通过上述各个实施例得到的特征参数(例如所述储层的渗透率表皮系数)与渗透率损害率是时空演化4D定量模拟的结果(如图4 所示)。更具体地，图5示出了由储层渗透率损害率表征的在第40天储层应力敏感性损害储层的半径的示意图，相关工作人员可通过该图5直观地确认储层被损害的程度。因此，可根据渗透率或表皮系数的演化特点进行储层损害定量预测和损害规律时空推演，为预防或避免储层损害、制定油藏的开发方案以及之后增产措施具有科学指导意义。

图6是本发明一实施例提供的应力敏感性损害储层的建模系统的结构图。如图6所示，所述建模系统包括：应力确定装置10，用于确定待诊断井的预设区域内的储层受到的有效应力；流量确定装置20，用于根据所述储层的孔隙的孔径分布特征、所述有效应力下的所述储层的各个毛管束的直径与长度、及流体流量公式，确定所述储层中的流体的流量，其中，所述毛管束由多个固体基质及所述多个固体基质之间的孔隙组成；以及模拟方程确定装置30，用于根据所述储层的渗透率模型及所述储层中的流体的流量，确定应力敏感性损害储层的时空演化模拟方程，其中所述时空演化模拟方程用于模拟由应力敏感性引起的储层损害特征的四维时空演化过程。

优选地，所述应力确定装置10可包括：压力方程建立模块，用于建立所述流体进入所述储层的压力传导方程；以及应力确定模块，用于根据所述储层的压力传导方程与所述储层的上覆岩层的压力，确定所述储层受到的有效应力。

优选地，所述流量确定装置20可包括：流量比率确定模块，用于根据所述有效应力下的所述储层的各个毛管束的直径与长度及流体流量比率公式，确定所述流体的流量比率；以及流量确定模块，用于根据所述储层的孔隙的孔径分布特征、所述储层中的流体的流量比率及所述流体流量公式，确定所述流体的流量。

优选地，所述流量确定装置20还包括：直径确定模块，用于根据所述储层的弹性模量，确定所述有效应力下的所述储层的各个毛管束的直径；以及长度确定模块，用于根据所述储层的弹性模量与泊松比，确定所述有效应力下的所述储层的各个毛管束的长度。

优选地，所述长度确定模块包括：初始长度确定单元，用于根据无有效应力下的所述储层的特定毛管束的表观长度、无有效应力下的所述储层的特定毛管束的初始直径及分形尺度法则，确定无有效应力下的所述储层的特定毛管束的初始长度L_c0；以及长度确定单元，用于根据所述储层的弹性模量E、所述储层的泊松比ν、无有效应力下的所述储层的特定毛管束的初始长度L_c0及下式，确定所述有效应力下的所述储层的特定毛管束的长度L_p，

优选地，所述流量确定装置还包括：横波速度确定模块，用于根据位于所述预设区域内的多个特定相邻井中的每一者的横波速度与纵波速度，确定所述储层的横波速度；以及泊松比确定模块，用于根据所述储层的横波速度，确定所述储层的所述弹性模量与所述泊松比。

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

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

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

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

所述确定储层损害程度的系统与上述确定储层损害程度的方法相对于现有技术所具有的优势相同，在此不再赘述。

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

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

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

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

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

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