曹仁义，程林松，杜旭林，时俊杰，杨晨旭

中国石油大学（北京）石油工程学院，北京 昌平 102249

引言

随着中国国民经济的快速发展，能源需求量越来越大，中国常规油气藏的产量已不能满足国家发展的需要。2020 年，中国石油对外依存度已超过70%，石油资源的保障压力越来越大，势必要发展新的石油资源战略接替区[1-3]。致密油藏主要是指与生油岩互层、紧邻的致密砂岩、致密碳酸盐岩等储集岩中，未经过大规模长距离运移的石油聚集[4]。根据现行的储层分类标准和国内外的勘探开发实践，在一般情况下，致密油储层孔隙度小于10%，基质覆压渗透率小于0.1 mD，单井无自然工业产能。致密油资源在中国的主要盆地广泛分布，大致上有陆相致密砂岩油藏、湖相碳酸盐岩油藏及泥灰岩裂缝油藏3 种类型[5]。相比低渗、特低渗储层，致密油藏的开发难度更大，由于其渗透率更低、物性更差，微纳米孔喉和微裂缝发育，加之目前采用水力压裂技术开发，其渗流环境更加复杂，传统渗流理论已不再适应，亟需发展与致密油藏高效开发相匹配的现代油气渗流理论，这同时也是中国致密油气开发领域中的重点攻关方向。

致密油气储集层的性质在很大程度上取决于其微观孔隙结构[6-7]。目前已用于储层微观孔隙结构分析的主要技术有：真实砂岩微观模型驱替实验技术、恒速压汞分析技术、核磁共振分析技术及CT 扫描技术等。杨正明等[8]利用核磁共振方法，研究了低渗透含水气藏岩芯的孔隙结构特征，结果表明渗透率不同，其孔道半径基本相同，而喉道半径不同。王瑞飞等[9]发现有效喉道半径、有效孔隙半径等特征参数与孔隙度、渗透率之间具有较好的相关性，在致密储层中影响流动的主要因素为：孔喉尺寸、孔喉分布、微观孔喉结构及矿物组成等。朱如凯等[10]采用扫描电子显微镜获得了致密砂岩的孔隙结构图片，发现致密储层微米级孔隙类型为粒间溶孔、颗粒溶蚀孔及微裂缝，孔隙直径主体在10～200µm，纳米级孔隙类型为原生粒间孔与自生矿物晶间孔，孔隙直径主体在70～400 nm。上述研究表明，致密储层中影响流动的主要因素包括孔喉尺寸、孔喉分布、微观孔喉结构及矿物组成等。

从高效开发方式来看，致密油藏的高效开发上依赖于“甜点区”或“甜点段”的识别与选取[11]。致密油藏“甜点区”是指在平面上成熟优质烃源岩分布范围内，具有工业价值的石油高产富集区。致密油藏“甜点段”是指在剖面上源储共生的黑色页岩层系内，经人工改造可形成工业价值的石油高产层段。一般具有较大分布范围、一定厚度规模、优质烃源岩、较好储层物性、较高含油气饱和度、较轻油质、较高地层能量、较高脆性指数及天然裂缝与局部构造发育等特征，在目前经济技术条件下，可优先进行勘探开发[12]。

致密油藏的高效开发倒逼致密储层复杂渗流理论的研究与创新。本文在分析致密油藏储层特征的基础上，阐述致密油藏的特殊渗流机理，并结合中国石油大学（北京）油藏数值模拟课题组的研究特色，综述致密油藏渗流规律与数学模型方面的最新研究进展及展望，对致密油藏的科学高效开发具有重要的理论指导意义。

1 微纳米孔喉流动机理及数学模型

由于致密油藏储层的渗流环境复杂，微纳米尺度下的流体已不再是连续态，而是以单个粒子的形式存在，并且在该尺度下粒子主要的运动方式为吸附、解吸及扩散，因此，不满足于泊肃叶方程和达西定律，需采用新方法和新模型来表征流体在微纳米尺度下的流动机理。

1.1 致密储层中的边界层效应

目前，研究普遍认为致密储层中出现非达西现象的根本原因在于固-液界面间的分子间作用力及液体分子间力的作用。国外方面，Tpe6ин[13]提出了石油渗流偏离达西定律的主要原因是表面活性组分（沥青质、胶质等）在岩石孔道表面上的吸附；杰弗里卡莫夫[14]用实验证实了流体与岩石表面产生吸附作用，使得渗透率急剧下降；Pertsin等[15]从理论上证明了圆柱体孔隙中的流体存在密度剖面；Wensink 等[16]采用分子模拟的方法论证了吸附层流体在轴向上存在密度和迁移率剖面。国内方面，黄延章[17]最早指出渗流流体的定义，强调渗流流体是渗流环境中的流体，包括体相流体和边界流体两个部分。体相流体是指其性质不受界面影响的流体，分布在多孔介质孔喉的中轴部位，边界流体是指其性质受界面影响的流体，它紧靠在孔喉壁上并表现为边界层形态。上述研究逐渐揭示了“边界层效应”的存在。对于致密油藏，其孔喉比较大，纳米喉道成为了制约流体流动的主要因素，由于原油与孔隙表面长期接触，会产生一种薄层流体吸附在孔喉表面难以流动的、性质不同于体相液体的吸附层，这种由原油中的极性物质组成的吸附层即为边界层[18]。致密油藏微纳米孔喉中的边界层形态如图1 所示。

图1 微纳米孔喉中的边界层形态Fig.1 Morphology of boundary layer in micro-nano pore throat

对于中、高渗透储层，由于原油边界层厚度小，所占比例极小，所以对渗流规律影响有限，在一定范围可将其视为牛顿流体。对于致密油藏，由于孔隙半径和边界层厚度几乎在同一数量级上或者更小，再加上多孔介质以及黏土的影响，这种固-液边界层的影响会大大增强，边界流体占孔隙中流体的比例增大，边界层的影响已经不容忽略，边界层的存在压缩了流体的有效流动空间，使得储层的有效流动半径减小，且边界流体难以流动，限制了体相流体的流动能力。虽然众多的学者对边界流体的属性和滑移长度做了大量研究，取得了丰富的成果，但是针对致密油藏储层微纳米孔喉与流体介质相互作用关系研究较少。由于微纳尺度下分子运动的影响对流体与岩石相互作用产生的壁面滑移即边界层效应不可忽略，如何定量表征边界层的影响并分析其中的主控因素，是研究致密油藏基质的关键渗流参数的重要基础。

1.2 微尺度流动物理模拟实验

早期的微尺度流动实验研究主要在半径100 µm 以上的圆管及狭缝中开展，材质为各类金属、塑料及玻璃等，流体为各类纯流体、溶液或气液两相流体，实验操作简单，可施加较高的压力梯度，因此，流动的雷诺数也较高，实验中会出现层流到紊流的流态变化。近几十年来，随着工业的发展，技术的进步和实验手段不断提高，生产与科研对微机电系统技术的要求也在不断提高，接连有各领域的学者开展了更微小尺度下的流体流动规律实验研究，采用的流体通道主要有长方形和微圆管两种，流体介质有去离子水、模拟油、石蜡及烷烃等，采用的实验材料、方法以及实验环境不同，实验结果的差异性也较大。李中锋等[19]、Cui 等[20]、岳湘安等[21]、徐绍良等[22]、李洋等[23]、宋付权等[24]、Wu 等[25-26]开展了水、烃类及聚合物等流体在半径1～30µm、内壁光滑的石英微管中的微尺度流动实验研究，论证了微尺度下的流量与压力梯度成非线性关系，并列出了多种非线性判定标准。典型微尺度流动实验研究参数见表1。

表1 典型微尺度流动实验研究Tab.1 Typical experimental study on microscale flow

物理模拟实验有其局限性，目前，微尺度流体流动实验研究面临以下挑战：（1）液体流动阻力大，在直接注入过程中难以维持稳定的压力和流量；（2）实验所需试件的加工及表征的难度大，此外还需要更高精度的测量设备；（3）受制于工业生产能力，微圆管的半径难以降低到1µm 以下，寻求能表征纳米级通道的最优材料十分重要。

基于微尺度流动物理模拟实验[25]，分析实验数据的统计学规律，认为边界层的比例（h/r）与喉道半径是指数关系，与压力梯度是乘幂关系，与流体黏度是线性关系，由此得出式（1）所示的边界层的厚度表征模型[36]，边界层效应的影响因素分析如图2所示。

图2 边界层效应的影响因素分析[36]Fig.2 Analysis of influence factors of boundary layer effect

1.3 耗散粒子动力学模拟

利用恒速压汞和高压压汞对致密油藏的孔喉分布进行研究发现，致密油藏的主流喉道半径一般小于2µm[9，37]。物理模拟实验存在一定局限性，由于实验尺度难以降低至微米级以下，达不到致密油藏基质的喉道分布范围，所以通过计算机模拟的方法来研究微尺度流动的重要性与紧迫性便凸显出来。纳米尺度下，比表面积的相对大小、流体与固壁粒子间的相互作用有着决定性作用，用传统岩芯尺度的连续介质力学无法揭示产生的滑移机理。因此，一些微观和介观尺度的研究方法得到了发展，目前用于模拟微尺度流动的常规方法，主要包括分子动力学方法（MD）、直接蒙特卡洛模拟方法（DSMC）、玻尔兹曼方程直接离散求解方法（BTE）及格子玻尔兹曼模拟方法（LBM）等[38]。不同的模拟方法基本原理及其优缺点对比如表2 所示。

表2 微尺度流动模拟方法对比Tab.2 Comparison of simulation methods for micro-scale flow

针对微尺度模拟计算量大的问题，本文介绍了一种新近发展起来的介观尺度模拟方法——耗散粒子动力学方法（Dissipative Particle Dynamics，DPD）。该方法以离散介质力学和牛顿运动定律为基础，与分子动力学同属于无网格方法，保留了粒子质量、速度、位置、不同粒子间的相互作用力等重要信息，忽略分子内部原子间相互作用等次要信息，大大减小了模拟的计算量，比格子玻尔兹曼方法更加灵活，且耗散粒子动力学方法在处理微尺度流动、复杂流体等粒子特性明显及壁面有滑移的情况比格子玻尔兹曼方法更有优势，因此，在模拟多相流动和复杂流体流动等多个方面表现突出。DPD 方法的应用主要集中于材料科学和医学领域[46-47]，在石油工程领域的应用主要集中于高聚物在微孔道中的运移和变形封堵机理研究[48-50]，在纳米喉道中边界负滑移流动方面研究较少，主要原因在于现有的DPD 模型大多没有考虑壁面分子对流体分子的特殊作用力，在模拟中未能充分考虑流固相互作用力的影响，因此，模型与真实物理情景存在一定差异。曹仁义等对原有DPD 模型进行改进，引入新的势函数表征壁面处的流固相互作用力，完善了纳米尺度下的油、气、水及岩石间的微观力场；优化了分子结构的粗化方法，提高了纳米离散介质流动模拟的尺度，从而进一步减小了模拟的计算量，目前已实现了不同流体在1～100 nm 尺度纳米喉道中的流动模拟[51-53]。

图3 展示了DPD 模拟得到的不同纳米喉道的密度剖面特征，可以看出，随着喉道半径的减小，密度分布的非均质性越强。x方向上，粒子密度越来越小，这是压力梯度在微观上的体现。y 方向上，固体壁面处有一层粒子密度明显超出其他位置，即为密度剖面角度观测到的边界层。图4 展示了不同喉道半径的速度剖面，可以看出，随孔喉尺度越小，速度分布逐渐不符合规则的抛物线形状，这是由于当喉道半径极小时，粒子运动速度很慢，分子热运动的影响逐渐显现。在靠近固体边界处存在厚度大约为几纳米的薄层，该范围内的粒子运动很不规则，为流体力学中的“双电层”，即为速度剖面角度观测到的边界层。边界层中的粒子其实也在一直运动，但受到了固体吸引力而始终在固体壁面周围做运动，无法顺畅地从喉道中流出。

图3 DPD 模拟的密度剖面[36]Fig.3 Density profile simulated by DPD

图4 DPD 模拟的速度剖面[36]Fig.4 Velocity profile simulated by DPD

根据Secchi 等[54]与Granick[55]的研究，边界层流体并非不流动，而是由于黏度升高，流体内摩擦力增加，流动速度较低，因此，学界一直存在另一种边界层假设，即动边界层假设。如果边界层流体参与流动，则其应有相应的流动剖面。假设体相流体的流动仍遵循纳维-斯托克斯方程，本文考虑边界层流体的黏度远高于体相流体且随着与固壁距离缩小而急剧升高的情况，基于DPD 模拟得到的幂函数规律来描述边界层流体的黏度剖面，根据流体内摩擦力与剪切力的平衡方程，推导了微尺度喉道（圆管）中考虑动边界层条件下的去离子水速度剖面模型[56]及流量公式。

2 致密油藏应力敏感及数学模型

在油藏开采前，储层岩石受上覆地层压力、孔隙流体压力以及岩石骨架应力的综合作用，一般处于平衡状态，但随着注入流体的进入或地层流体的采出，岩石骨架的受力平衡被打破，喉道受到有效应力时发生形变，从而会造成渗透率损失，这种性质称为储层的应力敏感性[57]。致密油藏储层普遍存在应力敏感现象，会造成油井生产一段时间后产量大幅递减，因此，定量分析应力敏感对油田生产的影响，有利于指导致密油藏的高效开发。

现阶段国内外学者所研究的应力敏感表征模型通常基于相关理论公式与实验结果的拟合，从渗透率与有效应力的关系表达形式来看，可以分为乘幂模型、指数模型、二项式模型及对数模型[58]，表征模型见表3。不同储层、不同岩性的多孔介质在相同的围压和孔隙流体压力作用下，表现出的有效应力亦有不同。针对不同渗流环境和不同假设条件，目前已发展了多种应力敏感数学模型，见表4。

表3 渗透率与有效应力关系表征模型Tab.3 Characterization model of the relationship between permeability and effective stress

表4 应力敏感数学模型对比Tab.4 Comparison of stress sensitive mathematical models

目前关于储层应力敏感的研究已有很多，但基于多孔介质微观孔喉结构形变及微尺度流动效应表征致密储层应力敏感的模型却极少。

本文结合鄂尔多斯盆地某致密油藏多组不同渗透率级别岩芯的恒速压汞实验结果，建立了式（4）所示考虑填隙物含量、填隙物类型及原始喉道变形的致密储层有效喉道压缩程度的表征模型[66-67]；考虑微尺度流动的边界层效应，建立了式（5）所示的致密储层有效渗透率的应力敏感数学模型[68-69]，该式将实验结果与孔喉变形机理相结合，实现了致密储层应力敏感的有效表征，反映了致密油藏有效应力、压力梯度、喉道分布及边界层性质对应力敏感特征的影响。

基于式（5）和行业标准SY/T 5358—2017 中的应力敏感损害率的定义，可得到式（6）所示的应力敏感的损害程度表征模型[68]。可以看出，致密油藏基质渗透率的损失程度由两部分组成，第一部分是由于喉道半径被压缩至临界喉道半径以下，导致喉道丧失全部渗流能力，见图5a；第二部分是由于喉道半径减小，但并未完全丧失渗流能力的部分渗透率损失，见图5b。

图5 应力敏感导致的渗透率损失示意图[68]Fig.5 Schematic diagram of permeability loss caused by stress sensitive

3 致密油藏非线性运动方程

致密储层孔喉细微，渗透率极低且渗流环境复杂，流体的流动规律某种程度上偏离达西定律，导致了低速非线性渗流特征的出现[70-72]。目前的非线性运动方程主要可分为分段式和连续式两类。拟启动压力梯度模型为典型的分段式[73]，是通过用直线代替非线性段来描述非线性渗流问题，本质上为启动压力梯度问题，只适用于非线性段不明显的常规低渗油藏。

Zeng 等[74]通过岩芯驱替实验证实了非线性的存在并确定了储层渗透率和拟启动压力梯度大小的关系；Bear 等[75]、冯文光等[76]通过将非线性段考虑成不同类型的分段函数来描述这一渗流规律，包括直线逼近、幂指函数逼近、二次函数逼近等[77-78]，但由于非线性段的划分难以提供明确的界限，并且由于其在表达式上具有不连续性和复杂性，因此，分段式不利于实际矿场应用。

针对分段式的局限性，邓英尔等[79]首次提出了考虑这种复杂的渗流规律的三参数连续式非线性运动方程；杨清立等[80]提出了精度更高的连续式运动方程，较好地反映了启动压力梯度和弯曲段；刘辉等[81]通过对长庆油田超低渗透油藏岩芯的实验结果进行拟合提出了新的运动方程，但是这种方法是通过实验方法回归参数得到，缺少对非线性渗流这种物理现象本质的探讨；时宇等[82]从毛细管模型出发，以边界层理论为基础，建立了能够定量反映造成非线性现象内在因素的非线性运动方程，并通过恒速压汞实验确定了非线性段的适应范围；Song等[83]同样基于边界层理论，提出了一种新的低渗透油藏指数型运动方程并利用实验数据进行了验证；姜瑞忠、杨仁锋等[84-86]论证了非线性渗流特征的必要性，并基于毛管束模型建立了低速非线性运动方程，同时结合相对渗透率的概念，将模型扩展到两相非线性渗流中，且分析了两相体系的局限性，在此基础上，建立了考虑非线性段连续变化的数值模拟方法并建立了考虑动边界的试井模型。除此之外，还有部分学者对非线性问题的数值求解进行了相关研究[87-89]。上述研究中的典型非线性运动方程，见表5。

表5 典型的非线性运动方程Tab.5 Typical nonlinear motion equations

在前人的研究中启动压力梯度的定义过分夸大了致密油藏中的渗流阻力，本文基于微纳米喉道边界层效应，提出了致密油藏中表征真实流动能力的有效喉道半径，见式（7）；基于毛管束模型，推导了考虑边界层效应与喉道迂曲度的圆管层流有效渗透率表达式，见式（8）；在此基础上，建立了可适用于描述致密油藏的非线性运动方程[54]，见式（9）。

表6 为不同喉道分布参数，图6 为致密油藏不同喉道分布下的非线性渗流特征，可以看出，随着喉道分布范围变小，致密油藏渗流的非线性程度增强，具体表现在：（1）开始有明显流量的压力梯度增大；（2）弯曲段曲率增大；（3）直线段出现的更晚。这是由于喉道分布的范围越小，边界层对渗流的影响越大，微尺度效应越强烈，宏观表现出来随着喉道分布范围缩小，非线性程度增强。

表6 不同喉道分布参数Tab.6 Parameters of different throat distribution

图6 致密储层非线性渗流特征[36]Fig.6 Nonlinear flow characteristics of tight formation

4 孔隙网络模型及非线性渗流规律

致密储层孔喉结构复杂，微纳米级孔喉发育。传统的毛管束模型在表征致密储层复杂孔喉结构方面存在缺陷，数字岩芯虽然能够准确表征致密储层的空间拓扑结构，但基于数字岩芯的流动模拟方法计算量大，计算时间长并且需要进行大量的并行计算，这对于非均质性很强的致密储层而言，流动模拟结果不一定可靠[90-91]。孔隙网络模型是将多孔介质用理想几何形状来表示，并用形状因子来表征孔喉结构的复杂性，其在模拟多孔介质孔隙级流动方面具有较大的优越性。

孔隙网络模型实质上是基于管流的渗透率模型，它是将多孔介质内流体的流动设想为在封闭管路网络内的流动。根据其不同的分类方式可分为多种类型，如据其发展方向分为三维随机网络模型[92-93]和真实拓扑孔隙网络模型[94-95]，据主导受力作用分为静态网络模型[96-97]和动态网络模型[98-99]。孔隙网络模型常用于分析微观剩余油分布规律及黏性指进现象，并获取相关的宏观物性参数，包括流速压差关系曲线、渗透率变化规律曲线、毛管力曲线及相渗曲线等[100-101]。近年来，不少学者将孔隙网络模型结合其他方法表征储层结构和渗流参数。Wang 等[102]将分子动力学与孔隙网络模型结合校正泊肃叶方程，提出了考虑页岩组成及双峰孔径分布（PSD）影响的页岩基质表观渗透率模型；杨永飞等[103]基于页岩储层富含有机孔隙与无机孔隙的特点，分别构建数字岩芯，将其整合为既包含无机孔隙信息又包含有机孔隙信息的双孔隙网络模型，并用于研究油水两相驱替流动过程[104-105]；孔隙网络模型还应用于其他跨尺度分析，王沫然等[106]基于孔隙网络模拟提出了一种跨尺度混合模拟算法，可用于研究宏观尺度的气井产量衰减曲线。目前，关于如何构建准确且高效的自适应孔隙网络模型已取得了较大进展，但由于致密地层中存在微纳米孔隙和高围压，在微观边界层与孔喉介质变形的双重影响下，亟需明确致密储层的非线性渗流规律及其主控因素。

基于此，本文采用最大球体法提取了规则的孔隙网络模型，实现了考虑边界层及应力敏感协同作用下微尺度流动效应的孔隙网络模拟[107-108]。模型中耦合边界层效应和应力敏感时，孔隙网络结构不再是固定的，孔隙网络的喉道半径随边界层厚度和有效应力的变化而变化。图7 展示了考虑边界层和应力敏感（介质变形）协同效应对孔隙网络压力分布的影响，能看出由于两者协同效应的存在，部分尺寸较小的喉道失效，会导致一定比例的流动空间损失。从图8a 可以看出，考虑边界层与介质变形情况下储层非线性渗流规律明显，同驱替压力梯度下的流体流速是最为缓慢的，流速只有理想情况下的1/4甚至更低，与单独考虑边界层或介质变形相对比，可以发现，两者效果均存在情况下会对储层储层流体流动造成更大的影响；当考虑边界层存在后，从图中可以清晰地看出存在明显的启动压力梯度，当驱替压力梯度小于此压力梯度时，流体处于不流动的状态。从图8b 能看出，在考虑边界层条件下，随着有效应力增加，流体流动空间压缩越多，流体流速也会减小。

图7 协同效应对孔隙网络压力分布的影响Fig.7 Influence of combined effect on pore-network pressure distribution

图8 协同效应对非线性渗流的影响[106]Fig.8 Influence of combined effect on nonlinear flow

致密储层非线性流速曲线的偏移程度往往取决于储层物性参数的好坏，当孔隙网络模型中考虑边界层效应及介质变形时，由于这两种机制对有效流动空间的压缩，导致储层流动阻力在低渗透率下显得更为明显。孔喉比指的是喉道与孔隙的半径比，孔隙大小决定了储层储量的多少，而喉道则决定了储层流动能力。

图9a 展示了不同孔喉比对非线性渗流的影响程度，当孔喉比为133 时，增加压力梯度，流速增加缓慢，说明在此孔喉比下，致密储层中流体流动困难。随着孔喉比增加，非线性渗流曲线斜率逐渐增大。孔喉比从133 提高至222 时，在同一驱替压力梯度下，流体流速增加达到数倍或甚至更大。配位数是反映多孔介质连通性的一个确定参数，从图9b所示的不同配位数对流速曲线的影响可以看出，模型的连通性随着配位数的增加而增大，导致绝对渗透率比增大，并且配位数增加也会增加储层微观波及系数，从而增大产量。

图9 储层参数对非线性渗流的影响[106]Fig.9 Influence of reservoir parameters on nonlinear flow

5 致密基质-裂缝耦合模型流动规律

致密油藏大规模压裂后产生毫米尺度的人工裂缝，同时储层可能天然发育有微米尺度的微裂缝，形成了致密基质-微裂缝-人工裂缝的多尺度复杂流动空间，因此，如何建立致密基质与裂缝耦合关系的数学模型，并准确地表征致密基质与裂缝之间的流动规律，这是模拟致密油藏开发过程的关键问题。致密基质与裂缝之间的耦合流动主要体现在两个方面，一是由于压差作用而导致的流动，二是由于渗吸作用而导致的流动。

5.1 致密基质-裂缝压差流动规律

由于致密基质与裂缝的流动能力差异极大，油藏开采过程中流体会在压差作用下在基质与裂缝间发生窜流。双重介质模型是刻画裂缝性油藏的传统方法，最早由Barenblatt[109]提出，基质与裂缝之间的压差窜流量由式（10）表示。

形状因子σ 是决定基质与裂缝之间压差流动的关键核心参数。国内外学者对形状因子的取值进行了不同的假设和研究。最初的研究基于拟稳态，给出的形状因子为恒定的、与时间无关的值[110-112]。而实际中基质-裂缝间的窜流过程并非完全满足拟稳态假设[113]（如较小速率的注水过程），因此，常量形状因子不能对此进行准确表征。随后，学者们提出形状因子是随时间变化的，周德华[114]等针对复杂的裂缝系统，采用解析方法建立了非均质裂缝与基岩间的窜流函数的表达式。He[115]在非稳态假设下推导了不同裂缝几何形状的形状因子并运用到产能预测和试井分析中。

Saboorian-Jooybari 等[116]考虑了毛管作用和重力的影响，研究了渗吸过程中随时间变化的基质-裂缝形状因子。Sun 等[117]通过分形压力扩散方程，导出了考虑分形维数和最大孔径的形状因子解析表达式。黄山等[118]假设基质的压力存在应力敏感，裂缝的压力随时间的增加指数递减，推导了考虑裂缝压力衰竭的时变形状因子。然而，与时间相关的形状因子的发展仍是当前的挑战，对于致密基质-裂缝系统需要考虑的情况应更加复杂。

致密储层由于孔喉尺度更小，流体流动更加困难，基质中流体向裂缝窜流过程中基质内部平均压力的位置会发生变化。由于致密油藏基质渗透率极低，压力在致密基质中传播速度较慢，导致从非稳态到稳态需要更长的过程，故致密油藏中稳态或拟稳态假设往往是不适用的。此外，致密油藏基质中流体的流动存在边界层、应力敏感等多种非线性因素。因此，本文认为研究致密油藏基质-裂缝之间的压差流动规律应该集中在3 个方面[119]：（1）由于致密基质与裂缝的渗透率差异极大，从而导致压力传播速度不同，因此，要对基质压力pm和裂缝压力pf的取值重新定义。（2）致密基质存在大量的纳米或微米尺度的喉道，需考虑边界层效应、应力敏感等非线性因素的影响，在窜流公式中不能用单一的渗透率值表征基质渗透率Km。（3）基于基质与裂缝之间不同的接触关系，需考虑瞬态压力的变化，推导致密油藏基质-复杂裂缝系统的时变形状因子。

对此，本文将基质与裂缝之间的不同接触关系分为3 种类型[119-120]：（1）基质与面缝（一维）；（2）基质与面缝-天然缝（二维）；（3）基质-体积压裂缝（三维）（图10）。

图10 致密基质与裂缝接触关系[119]Fig.10 Contact relationship between tight matrix and fracture

假设裂缝压力为裂缝面处的压力，基质压力取平均压力，距离选取基质宽度的一半来推导基质与裂缝的窜流方程，并引入窜流修正系数来弥补假设和真实情况的差距，模型考虑了致密储层渗流的非线性因素及基质中压力的非均匀分布，见式（11）。

不同接触情况下窜流修正系数的表达式如表7所示。图11 为不同渗透率下窜流修正系数随时间变化的曲线，可以看出，窜流修正系数随时间的增大而减小，最终趋于稳定，且基质渗透率越低，窜流修正系数达到稳定所需的时间越久，因此，对于致密油藏考虑基质-裂缝压差流动的瞬态变化影响是十分必要的。

图11 不同渗透率下窜流修正系数随时间变化的曲线[119]Fig.11 Curves of correction flow coefficient with time under different permeability

表7 致密基质-裂缝不同接触窜流修正系数[119]Tab.7 Correction flow coefficient in different contact between tight matrix and fracture

5.2 致密基质-裂缝渗吸流动规律

随着非常规油藏的兴起和压裂技术的广泛应用，渗吸作用受到越来越多学者的重视。对于致密基质-裂缝的渗吸流动规律，学者们已进行了大量探索，目前，关于渗吸流动的数学模型主要有归一化模型、解析与半解析模型以及数值模型。归一化模型是基于实验结果提出的经验模型，应用最早且最为广泛，Aronofsky 等[121]于1958 年最早将渗吸量表示为无因次时间的指数函数，但该模型过于简单，不能体现出渗吸过程的所有影响因素；在此基础上，Morrow 等[122]考虑岩芯尺寸、流体性质等影响因素，提出了新的渗吸流动模型；之后的众多学者[123-125]对无因次时间参数做了修正和改进，周林波等[126]将多孔介质中不同级别的孔隙考虑到了渗吸公式中，提出了三参数模型。由于归一化模型是源于实验结果的经验公式，缺少有效的理论依据，因此，解析及半解析模型方面，Kashchiev 等[127]根据假设的毛管力、相渗曲线函数推导得到了包含有渗吸项的控制方程；Li 等[128]提出了一个新的包含渗吸前期和后期的近似一维解析解；Cai 等[129-130]基于分形理论建立了包含重力因素的渗吸质量随时间变化的解析模型，并给出了裂缝性双重介质渗吸机理的分形判据模型和图版；Shi 等[131]建立了毛细管中存在非活塞式渗吸的分形模型，得到了不考虑流体静压的解析方程和考虑流体静压的半解析方程。解析及半解析模型的缺陷在于应用范围有限，对于多维及更复杂的问题可以采用数值模型解决。Pooladi-darvish 等[132]提出了渗吸数值模型并通过计算明确了正向和逆向渗吸的区别；Delijani 等[133]采用了格林元函数对一维条件下的油水渗吸采收率进行了数值计算；杨柳等[134]建立了压后焖井期间压裂液在毛细管力作用下自发渗吸进入致密油储层的数学模型并采用有限差分方法进行了求解；郑欢[135]提出了一种适用于毛管力作用下自发渗吸水油流动过程的双重介质模型窜流函数并对窜流量进行了精确计算。

事实上，分析致密油藏基质与裂缝间的渗吸流动过程，还要充分考虑致密孔喉的微尺度效应及致密基质-裂缝的复杂接触关系。对此，本文利用边界元和径向积分理论的径向积边界元方法，建立了考虑任意复杂形状基质块和边界层效应的渗吸数学模型[136-139]，并在岩芯和油藏尺度进行了验证，该方法无需将基质块进行网格划分，在保证数值模拟精度的同时提高了计算效率。

基质系统内饱和度的变化由毛管力控制，渗吸过程用非线性扩散方程描述

裂缝系统中，渗吸传质量作为油相和水相控制方程的源汇相

其中，使用逐次稳态替换法对真实毛管压力扩散系数进行处理，将强非线性项D(Sw)线性化，即在每一时间步内使用有效毛管力扩散系数De代替这段时间内真实的毛管力扩散系数D(Sw)。引入径向积边界元法得到离散格式的边界积分方程

图12 为致密基质-裂缝耦合流动含水饱和度场图，采用嵌入式离散裂缝模型结合PEBI 网格贴合人工水力压裂缝的方法，对鄂尔多斯盆地某致密油藏某多级压裂水平井进行周期50 d 的水吞吐开发模拟，分析油藏尺度下渗吸对致密油藏开发效果的影响。对比图12a 和图12b 可以看出，致密油藏水吞吐开采过程中，逆向渗吸的存在具有一定程度提高采油速度的作用；对比图12b 和图12c 可以看出，边界层的存在增大了毛管力，但也减小了储层渗透率且缩小了有效渗吸的范围，如果忽视了致密储层的边界层效应会严重高估油藏产能。

图12 致密基质-裂缝耦合流动含水饱和度场[139]Fig.12 Water saturation field of coupled flow in tight matrix and fracture

6 展望

6.1 致密油藏微纳米孔喉流动机理

现阶段关于微纳米级孔喉流动的分子或粒子模拟技术和方法已经取得了很大突破，但仍存在一定程度的局限性，一方面，从计算力学角度，表现为计算量大，模拟的分子或粒子数量有限，难以反映真实储层的空间维度和时间尺度，孔隙壁面岩石成分和流体组成复杂，多相多组分流体与复杂多孔介质构型的问题亟待解决。随着计算机水平的持续进步，对更为复杂的微纳米级孔喉流动机理的探索与认知是下一步攻关的重要方向。另一方面，关于物理模拟实验，目前的研究在实验模拟尺度至1µm 以下面临巨大挑战，主要原因在于扫描电镜等手段无法准确测量半径在1µm 以下的微管半径，且由于流量极低，现行的测量方法精度也不够。因此，未来在实验所用的微圆管材料及测量方法上需要进一步攻关，以便于实现更小空间尺度的实验探索，揭示致密油藏真实的微纳米孔喉流动机理。

6.2 纳米限域下的多种类分子孔隙网络模拟

孔隙网络模型相较于数字岩芯在致密油藏孔隙级流动规律研究方面具有较大的优势。致密储层孔喉尺寸极小，流体与岩石的相互作用对流体的相态、流体性质和流体的运移规律有很大的影响。目前，国内外学者已经针对致密油藏纳米孔隙中流体分子分布和运移情况进行了一定的研究，但是研究手段有限且缺乏实验验证，需要继续开展不同矿物组成的富有机质纳米孔隙中混合烷烃分子的运移规律模拟及实验研究，揭示致密储层基质中流体分子的运移规律，并利用孔隙网络模拟揭示纳米限域下的岩芯尺度流体分布情况及运移规律。致密油藏润湿性展布复杂，不同的润湿壁面流动机理存在较大差异，如何实现致密储层润湿性表征及微尺度流体在复杂润湿孔隙网络中的流动模拟也是研究致密油藏孔隙级流动规律的关键所在。

6.3 微观、非线性渗流机理与宏观数值模拟的结合

目前国内外学者已经开展了大量微观、非线性渗流机理研究，可以准确表征介质变形、边界层等复杂的渗流特征，但如何将此类微观、非线性渗流理论与宏观数值模拟准确结合，仍面临着诸多挑战，如当岩石内流动的空间缩小到微纳米尺度，传统的状态方程来表征流体的相态及PVT 属性的计算方法已经缺乏普遍适用性，这就意味着需进一步对微纳米尺度的流体相态特征的表征方法进行研究，并形成相应的宏观流体参数计算和数值模拟方法，以便于更准确地模拟致密油藏的生产过程。

符号说明

h—边界层厚度，µm；

r—喉道半径，µm；

a1、a2、a3—与流体性质相关的常量，当流体为去离子水时，a1=0.258，a2=-0.261，a3=-0.419；

∇p—压力梯度，mPa/m；

µ—流体黏度，mPa·s；

vmicrotube—微圆管内流体流速，mm/s；

x—圆管横截面上的某一坐标位置，µm；

r0—圆管半径某一个值，µm；

b1、b2—黏度剖面幂函数的两个待定系数，无因次；

qmicrotube—微圆管内流体流量，×10−6µL/s；

K—储层渗透率，D；

K0—储层参考渗透率，D；

σeff—有效应力，MPa；

σ0—参考有效应力，MPa；

α—乘幂模型应力敏感系数，无因次；

β—指数模型应力敏感系数，MPa−1；

z1—二项式模型二次项应力敏感系数，D/MPa2；

z2—二项式模型一次项应力敏感系数，D/MPa；

z3—二项式模型零次项应力敏感系数，D；

S—对数模型应力敏感系数，无因次；

pk—上覆地层压力，MPa；

ph—裂缝闭合应力，MPa；

p0—参考油藏应力，MPa；

H—裂缝表面粗糙度高度分布的标准差，µm；

a0—裂缝参考半开度，µm；

p—油藏压力，MPa；

C0—和岩石性质相关的常量，无因次；

pi—岩石持久变形的等价压力，MPa；

Cp—孔隙压缩系数，MPa−1；

Δσ—地层有效应力差，MPa；

φ0—参考孔隙度，无因次；

ν—泊松比，无因次；

r∗—圆形管束模型截面圆半径，µm；

E—储层弹性模量，MPa；

e—椭圆形管束模型长轴，µm；

ε—椭圆形管束模型纵横比，无因次；

Kg—气相渗透率，D；

pp—孔隙流体压力，MPa；

Rg—气体常数，8.314 J/mol·K；

T—热力学温度，K，K=°C+273.15；

M—气体分子质量，kg/kmol；

ci—半径为ri的喉道受到有效应力作用后的压缩程度，无因次；

ri—第i组实验所测得的喉道半径，≤ri≤rmax，µm；

K′—受到有效应力作用后的致密砂岩岩芯有效渗透率，D；

τ—迂曲度，无因次；

rc—原始喉道半径，µm；

rmax—最大喉道半径，µm；

φ′—受到有效应力作用后的致密砂岩岩芯有效孔隙度，，无因次；

ni—毛管束的个数；

Ac—岩芯的横截面积，µm2；

φ—致密砂岩岩芯原始孔隙度，无因次；

D—致密砂岩岩芯应力敏感系数，无因次；

hi—喉道半径ri处的原始边界层厚度，µm；

v—渗流速度，m/s；

η1、η2、η3、η4、η5、η6、η7、η8—单位换算系数，无因次；

λb—启动压力梯度，MPa/m；

Δp—渗流压差，MPa；

L—渗流长度，m；

f—指数系数，无因次；

λc—结束非线性段时的压力梯度，Pa/m；

A1、A2—实验参数，×1012（Pa·s）/m2；

A3—实验参数，×106s/m；

B1—非线性渗流影响因子，无因次；

B2—拟启动压力梯度的倒数，m/MPa；

d1、d2、d3、d4—与渗透率相关的影响参数，d1=−7×10−8K3.5382，d2=3×10−5K1.8958，d3=0.0013K2.3537，d4=−0.0018K+0.0035；

c—拟合系数，无因次；

rmax—最大喉道半径，µm；

rmin—最小喉道半径，µm；

∇pmin—最小压力梯度，MPa/m；

ξ1、ξ1—拟合参数，MPa/m；

reff—有效喉道半径，µm；

Keff—致密砂岩岩芯有效渗透率，D；

σHg—汞气表面张力，N/m；

θHg—汞气接触角，（°）；

Ss—不可动流体饱和度，Ss=+(1−SHgmax)，无因次；

SHg—进汞饱和度，无因次；

SHgmax—进汞最大饱和度，无因次；

qmf—基质与裂缝之间的窜流量，m3；

σ—形状因子，m−2；

Km—基质渗透率，mD；

V—基质块的体积，m3；

B—体积系数，m3/m3；

pm—基质压力，MPa；

pf—裂缝压力，MPa；

Cf—窜流修正系数，无因次；

A—基质的横截面面积，m2；

lm—基质的长度，m；

tD—无因次时间，无因次；

t—真实时间，s；

φm—基质孔隙度，无因次；

ct—总压缩系数，MPa−1；

αn—n阶贝塞尔函数的根，无因次；

n—贝塞尔函数的阶数，无因次；

R—等效球形半径，m；

D(Sw)—毛管力扩散系数，m2/s；

Sw—基质中的含水饱和度，无因次；

λo—油的流度，m2/(Pa·s)；

λw—水的流度，m2/(Pa·s)；

So—基质中的含油饱和度，无因次；

Cav—平均流体压缩系数，Pa−1；

λav—平均流体流度，m2/(Pa·s)；

pc—毛管压力，Pa；

Vf—裂缝单元的体积，m3；

Kf—裂缝渗透率，D；

Kro—油的相对渗透率，无因次；

µo—油的黏度，Pa·s；

pfo—裂缝中的油相压力，Pa；

qoim—因渗吸作用产生的“源”项，m3/s；

qfo，qfw—因采油等（非渗吸因素）产生的源汇项，m3/s；

φf—裂缝孔隙度，无因次；

Sfo—裂缝中的含油饱和度，无因次；

Krw—水的相对渗透率，无因次；

µw—水的黏度，Pa·s；

pfw—裂缝中的水相压力，Pa；

qwim—因渗吸作用产生的“汇”项，m3/s；

Sfw—裂缝中的含水饱和度，无因次；

qwbj—第j个边界单元的渗吸速度，m3/s；

Ci—第i个网格的函数值，无因次；

Si—第i个网格的含水饱和度，无因次；

N—基质块边界被划分的片段单元数，无因次；

j—第j个片段单元，无因次；

Swbj—第j个边界单元的含水饱和度，无因次；

G(x,y)—基本解；

x—场点；

y—源点；

De—有效毛管力扩散系数，m2/s；

Γ—计算区域的边界；

ni—第i个边界的外法向量，无因次；

m—第m个时间步，无因次；

G′(x,y)—基本解的一阶导数；

C(y)—点y 的函数值，无因次；

Ω—计算域；

ω—场点与源点间的距离，m。