考虑固液分子作用的多孔介质动态网络模型

第36卷第2期 北京科技大学学报 Vol.36 No.2 2014年2月 Journal of University of Science and Technology Beijing Feb.2014 考虑固液分子作用的多孔介质动态网络模型 王小锋)区，朱维耀”，邓庆军)，刘启鹏”，隋新光》，娄钰”，高英”， 张雪龄) 1)北京科技大学土木与环境工程学院，北京1000832)东北石油大学，大庆163000 3)大庆油田第一采油厂，大庆163000 ☒通信作者，Emai:bkwxiaofeng@163.com 摘要根据质量守恒定律，从固液分子力作用出发，结合岩心的孔喉分布特征，建立考虑固壁与流体分子作用的动态网络 模型.利用模型模拟得到不同时间步的孔隙压力变化，分析了固液分子间力、喉道半径、润湿性等因素对剩余油饱和度分布和 水驱相渗曲线的影响.研究结果表明，考虑固液分子间力作用时剩余油饱和度增大，岩石骨架与流体间的分子作用在多孔介 质的孔隙流动中不可忽略 关键词油田开发：多孔介质：渗流：分子间力：孔压：数学模型 分类号TE312:0359 Dynamic network model considering solid-iquid molecule interaction in porous media WANG Xiao feng",ZHU Wei-yao",DENG Qing jun?,LIU Qi-peng",SUI Xin-guang,LOU Yu,GAO Ying", ZHANG Xue-ing” 1)Civil &Environment Engineering School,University of Science and Technology Beijing,Beijing 100083,China 2)Northeast Petroleum University,Daqing 163000,China 3)The First Oil Production Plant,Daqing Oilfield,Daqing 163000,China Corresponding author,E-mail:bkwxiaofeng@163.com ABSTRACT According to the mass conservation law,a dynamic network model which considers interaction between solid walls and liquid molecules was established.This model is on the basis of solid-liquid intermolecular force and in combination with the distribution characteristics of pore throats in porous media.The change in pore pressure at different time steps was simulated through the model. The influences of solid-liquid intermolecular force,pore throat radius and wettability on the saturation distribution of remaining oil and the permeability curve of water flooding phase were analyzed.The results show that when considering solid-liquid intermolecular force the saturation of residual oil increases,thus interaction between rock and fluid molecules cannot be ignored in the study of porous flow in porous media. KEY WORDS oil field development;porous media:porous flow;molecular force:pore pressure:mathematical models 多孔介质储集层孔隙结构复杂极不规则，流体 和规律.在国外，1956年Fat0最早提出了网络模 在储层内的流动规律很难用几个参数来描述.为了 型的概念，并利用网络模型研究了多孔介质中流体 探知流体在多孔介质内的微观流动机理，研究人员 的动、静态性质：Dullien等p)使网络模型与逾渗理 把多孔介质的孔喉结构抽象为具有理想形状的几何 论相结合，研究了流体在孔隙空间内的窜通规则 空间，建立孔隙网络模型模拟多孔介质的流动过程 国内外有很多学者-国利用网络模型研究润湿性对 收稿日期：201308-22 基金项目：国家重大基金资助项目(50934003,11372033)：教育部专项资助项目(FRF-MP-B12006B) DOI:10.13374/j.issn1001-053x.2014.02.002:http://journals.ustb.edu.cn

第 36 卷 第 2 期 2014 年 2 月 北京科技大学学报 Journal of University of Science and Technology Beijing Vol． 36 No． 2 Feb． 2014 考虑固液分子作用的多孔介质动态网络模型 王小锋1) ，朱维耀1) ，邓庆军2) ，刘启鹏1) ，隋新光3) ，娄 钰1) ，高 英1) ， 张雪龄1) 1) 北京科技大学土木与环境工程学院，北京 100083 2) 东北石油大学，大庆 163000 3) 大庆油田第一采油厂，大庆 163000  通信作者，E-mail: bkwxiaofeng@ 163． com 摘 要 根据质量守恒定律，从固液分子力作用出发，结合岩心的孔喉分布特征，建立考虑固壁与流体分子作用的动态网络 模型． 利用模型模拟得到不同时间步的孔隙压力变化，分析了固液分子间力、喉道半径、润湿性等因素对剩余油饱和度分布和 水驱相渗曲线的影响． 研究结果表明，考虑固液分子间力作用时剩余油饱和度增大，岩石骨架与流体间的分子作用在多孔介 质的孔隙流动中不可忽略． 关键词 油田开发; 多孔介质; 渗流; 分子间力; 孔压; 数学模型 分类号 TE 312; O 359 Dynamic network model considering solid-liquid molecule interaction in porous media WANG Xiao-feng1)  ， ZHU Wei-yao1) ， DENG Qing-jun2) ， LIU Qi-peng1) ， SUI Xin-guang3) ， LOU Yu1) ， GAO Ying1) ， ZHANG Xue-ling1) 1) Civil ＆Environment Engineering School，University of Science and Technology Beijing，Beijing 100083，China 2) Northeast Petroleum University，Daqing 163000，China 3) The First Oil Production Plant，Daqing Oilfield，Daqing 163000，China  Corresponding author，E-mail: bkwxiaofeng@ 163． com ABSTＲACT According to the mass conservation law，a dynamic network model which considers interaction between solid walls and liquid molecules was established． This model is on the basis of solid-liquid intermolecular force and in combination with the distribution characteristics of pore throats in porous media． The change in pore pressure at different time steps was simulated through the model． The influences of solid-liquid intermolecular force，pore throat radius and wettability on the saturation distribution of remaining oil and the permeability curve of water flooding phase were analyzed． The results show that when considering solid-liquid intermolecular force the saturation of residual oil increases，thus interaction between rock and fluid molecules cannot be ignored in the study of porous flow in porous media． KEY WOＲDS oil field development; porous media; porous flow; molecular force; pore pressure; mathematical models 收稿日期: 2013--08--22 基金项目: 国家重大基金资助项目( 50934003，11372033) ; 教育部专项资助项目( FＲF-MP-B 12006B) DOI: 10． 13374 /j． issn1001--053x． 2014． 02． 002; http: / /journals． ustb． edu． cn 多孔介质储集层孔隙结构复杂极不规则，流体 在储层内的流动规律很难用几个参数来描述． 为了 探知流体在多孔介质内的微观流动机理，研究人员 把多孔介质的孔喉结构抽象为具有理想形状的几何 空间，建立孔隙网络模型模拟多孔介质的流动过程 和规律． 在国外，1956 年 Fatt［1］最早提出了网络模 型的概念，并利用网络模型研究了多孔介质中流体 的动、静态性质; Dullien 等［2--3］使网络模型与逾渗理 论相结合，研究了流体在孔隙空间内的窜通规则． 国内外有很多学者［4--13］利用网络模型研究润湿性对

·146 北京科技大学学报 第36卷 剩余油分布的影响并预测毛管压力曲线。在国内， 在岩石表面上的分布，对孔隙中毛管力的大小和方 王金勋等4利用网络模型研究了孔喉结构参数 向也有影响，表面张力越大，毛管力越大.因而，直 对两相相对渗透率的影响，胡雪涛和李允6利用随 接影响到流体在多孔介质中的渗流规律 机网络模型研究了孔隙结构对微观剩余油分布的影 响.目前，网络模型内流体流动都遵循泊肃叶定律， 2考虑分子力作用的微圆管运动方程 没有考虑在孔喉尺度下流体与岩石骨架的相互作 考虑流体与管壁间的分子作用，流体间的分子 用9，但是研究发现，流体在直径小于几百微米 引力增大，流体黏度增大，假设流体黏度由两部分组 的微管中流动时，实验与经典理论泊肃叶定律预测 成，一部分是不考虑固液分子作用时流体的黏度 的结果产生了偏离。这是因为孔喉的尺寸微小，孔 ,另一部分是由于固液长程范德华力增大产生的 喉壁面与流体之间分子作用己到了不可忽略的程 黏度，流体黏度可表示为0 度.为此，本文建立考虑流体与岩石骨架之间分子 力作用的水驱动态网络模型，模拟水驱后剩余油的形 u=4+6A4-A) (1) 成及分布，研究固液分子间力、喉道半径、润湿性等因 建立极坐标系下的NS方程，结合式(1)，推导 素对剩余油的影响，揭示多孔介质内的微观流动机理 可得流体在微圆管内的运动方程： 以及微观剩余油形成机制，获取宏观流动的重要参数， 为多孔介质细观流动理论的建立奠定理论基础 Q=-日，E(1-8) (2) 80 d 1多孔介质中的微观力 其中： 46(A,A-A)4 1.1分子力 In BoR+b(A.A.-A.) R b(JAA.-A.) 分子间力是指分子与分子之间或高分子化合物 分子内官能团之间的作用力.分子间力又可以分为 12 (AA-A.)I R+（A,A-A) 三种力，即取向力、诱导力和色散力.流体在微通道 MoR （JAA-A) 内流动时，由于固体表面与流体分子间力作用的影 463 响，微通道中的流体产生附加黏度，使流体的流动速 器(A) 度降低.这是实验测量数据与泊肃叶定律的理论值 偏差的根本原因 红-4)h么R+6(红-A b(A.A.-A.) 1.2静电力 12b3 静电力是带电分子或粒子间的作用力，它的大 (A不-A.)nR+6(AA-A) b(A.A.-A.) 小与距离的平方成反比.溶液在微通道中流动时， 22b （A.A-A). 固壁面与含有带电离子的溶液接触，溶液中会形成 3uoR 净电荷，当溶液受到压力差驱动时，这些净电荷随溶 若不考虑管壁与流体的固液分子作用时，即 液运动，在微通道两端形成流动电位势，产生流动诱 A,=A,则式(2)可退化为泊肃叶定律的形式，即 导电场，流动诱导电场作用在微通道的净电荷上，使 Q=-π.R虫 (3) 其向流动反方向运动，并通过黏性力带动溶液一起 8 Lo dx 反向流动，这样使溶液的正向流速减慢，就像是溶液 式中：为不考虑固液分子作用时水的黏度，mPa· 的黏度增大那样，形成“电黏性效应”，微通道尺寸 s;b是管壁与水分子作用起黏度增加的系数，它 越小，电黏性效应越明显 与固液间长程范德华力相关，本文中取值为5×103 1.3表面张力 Pasm;A.和A,分别代表水相和固体管壁的哈默 表面张力是由于液体表面层分子引力不均衡引 克常数，J;d为离固体表面的距离，m;E是分子作用 起的.表层分子比内部分子储存有多余的自由能， 系数：R为微圆管半径，m;dp/x是圆管的压力梯 如果要把内部的分子举升到表面上，需要付出能量 度，Pa·m-1;Q是通过圆管截面的流量，m. 做功，这种能量转化为表面自由能.表面自由能有 3多孔介质随机网络模型 趋于最小的趋势，液体表面就会有自动收缩，这种收 缩力称为表面张力.表面张力与温度、界面性质有 3.1网络模型的几何结构 关，油层中流体的表面张力将直接影响油层中流体 多孔介质的真实结构可以用三维的网络模型来

北 京 科 技 大 学 学 报 第 36 卷 剩余油分布的影响并预测毛管压力曲线． 在国内， 王金勋等［14--15］利用网络模型研究了孔喉结构参数 对两相相对渗透率的影响，胡雪涛和李允［16］利用随 机网络模型研究了孔隙结构对微观剩余油分布的影 响． 目前，网络模型内流体流动都遵循泊肃叶定律， 没有考虑在孔喉尺度下流体与岩石骨架的相互作 用［17--19］，但是研究发现，流体在直径小于几百微米 的微管中流动时，实验与经典理论泊肃叶定律预测 的结果产生了偏离． 这是因为孔喉的尺寸微小，孔 喉壁面与流体之间分子作用已到了不可忽略的程 度． 为此，本文建立考虑流体与岩石骨架之间分子 力作用的水驱动态网络模型，模拟水驱后剩余油的形 成及分布，研究固液分子间力、喉道半径、润湿性等因 素对剩余油的影响，揭示多孔介质内的微观流动机理 以及微观剩余油形成机制，获取宏观流动的重要参数， 为多孔介质细观流动理论的建立奠定理论基础． 1 多孔介质中的微观力 1. 1 分子力 分子间力是指分子与分子之间或高分子化合物 分子内官能团之间的作用力． 分子间力又可以分为 三种力，即取向力、诱导力和色散力． 流体在微通道 内流动时，由于固体表面与流体分子间力作用的影 响，微通道中的流体产生附加黏度，使流体的流动速 度降低． 这是实验测量数据与泊肃叶定律的理论值 偏差的根本原因． 1. 2 静电力 静电力是带电分子或粒子间的作用力，它的大 小与距离的平方成反比． 溶液在微通道中流动时， 固壁面与含有带电离子的溶液接触，溶液中会形成 净电荷，当溶液受到压力差驱动时，这些净电荷随溶 液运动，在微通道两端形成流动电位势，产生流动诱 导电场，流动诱导电场作用在微通道的净电荷上，使 其向流动反方向运动，并通过黏性力带动溶液一起 反向流动，这样使溶液的正向流速减慢，就像是溶液 的黏度增大那样，形成“电黏性效应”，微通道尺寸 越小，电黏性效应越明显． 1. 3 表面张力 表面张力是由于液体表面层分子引力不均衡引 起的． 表层分子比内部分子储存有多余的自由能， 如果要把内部的分子举升到表面上，需要付出能量 做功，这种能量转化为表面自由能． 表面自由能有 趋于最小的趋势，液体表面就会有自动收缩，这种收 缩力称为表面张力． 表面张力与温度、界面性质有 关，油层中流体的表面张力将直接影响油层中流体 在岩石表面上的分布，对孔隙中毛管力的大小和方 向也有影响，表面张力越大，毛管力越大． 因而，直 接影响到流体在多孔介质中的渗流规律． 2 考虑分子力作用的微圆管运动方程 考虑流体与管壁间的分子作用，流体间的分子 引力增大，流体黏度增大，假设流体黏度由两部分组 成，一部分是不考虑固液分子作用时流体的黏度 μ0，另一部分是由于固液长程范德华力增大产生的 黏度，流体黏度可表示为［20］ μ = μ0 + b ( 槡AsAw － Aw ) d ． ( 1) 建立极坐标系下的 N-S 方程，结合式( 1) ，推导 可得流体在微圆管内的运动方程: Q = － π 8 ·Ｒ4 μ0 ( 1 － ε) dp dx ． ( 2) 其中: ε = 4b 4 ( 槡AsAw － Aw) 4 μ4 0Ｒ4 ln μ0Ｒ + b( 槡AsAw － Aw) b( 槡AsAw － Aw) + 12b 2 ( 槡AsAw － Aw) 2 μ2 0Ｒ4 ln μ0Ｒ + ( 槡AsAw － Aw) ( 槡AsAw － Aw) － 4b3 μ3 0Ｒ3 ( 槡AsAw － Aw ) 3 － 10b2 μ2 0Ｒ2 ( 槡AsAw － Aw ) 2 + b μ0Ｒ( 槡AsAw － Aw ) ln μ0Ｒ + b( 槡AsAw － Aw ) b( 槡AsAw － Aw ) + 12b3 μ3 0Ｒ3 ( 槡AsAw － Aw ) 3 ln μ0Ｒ + b( 槡AsAw － Aw ) b( 槡AsAw － Aw ) － 22b 3μ0Ｒ( 槡AsAw － Aw ) ． 若不考虑管壁与流体的固液分子作用时，即 As = Aw，则式( 2) 可退化为泊肃叶定律的形式，即 Q = － π 8 ·Ｒ4 μ0 ·dp dx ． ( 3) 式中: μ0 为不考虑固液分子作用时水的黏度，mPa· s; b 是管壁与水分子作用引起黏度增加的系数，它 与固液间长程范德华力相关，本文中取值为 5 × 1013 Pa·s·m; Aw 和 As 分别代表水相和固体管壁的哈默 克常数，J; d 为离固体表面的距离，m; ε 是分子作用 系数; Ｒ 为微圆管半径，m; dp /dx 是圆管的压力梯 度，Pa·m － 1 ; Q 是通过圆管截面的流量，m3 ． 3 多孔介质随机网络模型 3. 1 网络模型的几何结构 多孔介质的真实结构可以用三维的网络模型来 · 641 ·

第2期 王小锋等：考虑固液分子作用的多孔介质动态网络模型 ·147· 描述，网络模型由喉道及其相连的孔隙体构成.如 2所示. 图1是一个二维的5×5网络模型示意图，图中黑色 圆形截面部分代表孔隙，与孔隙相连的部分为喉道， 相邻两个孔隙中心之间的距离为L.多孔介质孔隙 和喉道尺寸的大小一般满足数理统计上的几种分 布，比如常见的有瑞利分布、威布尔分布和对数正态 分布.本文中网络模型中孔喉尺寸满足的分布为威 布尔分布 P.(x)=B (x-y-e，x≥y (4) 图2湿相与非湿相占据的孔隙截面积 0, x<y. Fig.2 Wetting fluid and non-wetting fluid residing in crevices in the 式中，形状参数m>0,尺度参数B>0,y为位置 pore space 参数 非湿相所占据的孔隙截面积： 孔隙 AA-[cotarcos (a+- (a+0m)+cos(a+0m)sin(a+ (8) 式中，A为孔隙截面积，A为非湿相截面积，A.为湿 相截面积，「为湿相和非湿相两相交界面的曲率 半径，0a为润湿角，α为角隅半角.图3表示的是 喉道 矩形截面孔隙的角隅示意图. 图1二维孔隙与喉道网络模型 Fig.1 A network of 2-dimensional pores and throats 3.2毛管力模型 毛管力指的是两相界面上的压力差，数值上等 于界面两侧非湿相压力减去湿相压力，通常用 Young-aplace方程来表示： -n.=a(片+片) (5) 图3矩形截面孔隙的角隅 式中，P为非湿相压力，P.为湿相压力，σ为两相间 Fig.3 Comer geometry in pores with the rectangular cross-section 界面张力，”1和2分别为任意曲面的两个主曲率半 径.对于圆形截面的孔隙，毛管力可表示为 若考虑岩石与流体之间的分子力作用，那么网 2ocose P。s5 (6) 络模型中孔隙之间的流动遵循式(2)微管内的运动 r 方程，两个相邻孔隙i和j之间的导流系数G可表 对于矩形截面的孔隙，毛管力可表示为 示为 0+cos20-元-sin0cos0 Pe= 4 (9) (7) 6,=-g1-e)4 cos0-/4 -0+sin0cos0 一8 hroat,j 式中，0为润湿接触角 (10) 3.3有效导流系数 在真实的多孔介质当中孔隙和喉道有很多是由 那么相邻孔隙间的流量Q可表示为 非湿相与湿相流体共同占据的，湿相流体占据孔喉 Qg=G‘△p. (11) 的角隅位置，非湿相对于非圆形截面的孔喉，比如矩 式中，rH为有效半径，Lhmt.g为喉道长度，Thm.是喉 形截面的孔隙，湿相与非湿相所占据的截面积如图 道半径，ε为分子作用系数，为平均黏度

第 2 期 王小锋等: 考虑固液分子作用的多孔介质动态网络模型 描述，网络模型由喉道及其相连的孔隙体构成． 如 图 1 是一个二维的 5 × 5 网络模型示意图，图中黑色 圆形截面部分代表孔隙，与孔隙相连的部分为喉道， 相邻两个孔隙中心之间的距离为 L． 多孔介质孔隙 和喉道尺寸的大小一般满足数理统计上的几种分 布，比如常见的有瑞利分布、威布尔分布和对数正态 分布． 本文中网络模型中孔喉尺寸满足的分布为威 布尔分布． pw ( x) = m β ( x － γ) m － 1 e － ( x － γ) m β ， x≥γ; 0， x ＜ γ { ． ( 4) 式中，形状参数 m ＞ 0，尺度参数 β ＞ 0，γ 为位置 参数． 图 1 二维孔隙与喉道网络模型 Fig． 1 A network of 2-dimensional pores and throats 3. 2 毛管力模型 毛管力指的是两相界面上的压力差，数值上等 于界面两侧非湿相压力减去湿相压力，通 常 用 Young-Laplace 方程来表示: pc = pnw － pw = σ ( 1 r1 + 1 r ) 2 ． ( 5) 式中，pnw为非湿相压力，pw 为湿相压力，σ 为两相间 界面张力，r1 和 r2 分别为任意曲面的两个主曲率半 径． 对于圆形截面的孔隙，毛管力可表示为 pc = 2σcosθ r ; ( 6) 对于矩形截面的孔隙，毛管力可表示为 pc = σ ( r θ + cos2 θ － π 4 － sinθcosθ cosθ － 槡π/4 － θ + sinθcos ) θ ． ( 7) 式中，θ 为润湿接触角． 3. 3 有效导流系数 在真实的多孔介质当中孔隙和喉道有很多是由 非湿相与湿相流体共同占据的，湿相流体占据孔喉 的角隅位置，非湿相对于非圆形截面的孔喉，比如矩 形截面的孔隙，湿相与非湿相所占据的截面积如图 2 所示． 图 2 湿相与非湿相占据的孔隙截面积 Fig． 2 Wetting fluid and non-wetting fluid residing in crevices in the pore space 非湿相所占据的孔隙截面积: Anw = A － [ ncornersr 2 w/ nw [ cotαcos( α + θw/ nw ) － ( π 2 － ( α + θw/ nw ) ) + cos( α + θw/ nw ) sin( α + θw/ nw ] ] ) ． ( 8) 式中，A 为孔隙截面积，Anw为非湿相截面积，Aw 为湿 相截面积，rw/ nw为湿相和非湿相两相交界面的曲率 半径，θw/ nw为润湿角，α 为角隅半角． 图 3 表示的是 矩形截面孔隙的角隅示意图． 图 3 矩形截面孔隙的角隅 Fig． 3 Corner geometry in pores with the rectangular cross-section 若考虑岩石与流体之间的分子力作用，那么网 络模型中孔隙之间的流动遵循式( 2) 微管内的运动 方程，两个相邻孔隙 i 和 j 之间的导流系数 Gij可表 示为 Gij = － π 8 ( 1 － ε) r 4 eff μLthroat，ij ， ( 9) reff = ( 1 2 Anw 槡π + rthroat，ij ) ． ( 10) 那么相邻孔隙间的流量 Qij可表示为 Qij = Gij·Δp． ( 11) 式中，reff为有效半径，Lthroat，ij为喉道长度，rthroat，ij是喉 道半径，ε 为分子作用系数，μ 为平均黏度． · 741 ·

·148 北京科技大学学报 第36卷 3.4网络模型流动机理 假定流体是可不压缩的，网络模型结构如图1 K=- nuLo 所示，孔隙节点数为n×n,相邻两个孔隙中心的距 =An(P-P)=A(P.-P)· (18) 离为L,孔隙喉道半径的大小满足威布尔分布，网络 假设湿相为水相，非湿相为油相，那么由两相流 模型充满非湿相和湿相，非湿相占据中心位置，湿相 动时的模拟结果可以得到模型出口端的水相流量 占据角隅，忽略孔隙内的毛管压力，考虑喉道内的毛 Q。和油相流量Q。,对比单相流动时的流量Q可以 管力作用，那么根据质量守恒原理，注人孔隙的流量 得到模型的相对渗透率 之和应等于流出孔隙的流量，即 水相相对渗透率： 月9=0 (12) Ki= (19) Q 式中：9指的是孔隙i与相邻喉道连接的孔隙j之间 油相相对渗透率： 的流量；c是与孔隙i连接的喉道个数，即孔隙的配 (20) 位数. 飞-号 若相邻孔隙i和之间的喉道内只有单相流动 式中：S是孔隙i的含水饱和度；q是第n+1个 时，由导流系数的定义可知，通过该喉道的流量可表 时间步时水相进入孔隙i的流量，m3;V:是孔隙i的 示为 体积，m3;Am是网络模型的截面积，m2;q:是出口端 9i=G(p:-P）. (13) 孔隙i的流量，m3;P,和P,分别是模型入口端压力 湿相驱替非湿相时，若孔隙被驱替相充满，那 和出口端压力，Pa. 么连接孔隙i和j的喉道存在油水界面，此时通过该 3.5网络模型压力计算 喉道的流量可表示为 研究网络模型最重要的过程在于求解孔隙的压 Gi (pi-pj+Pe),Pi-pi+Pe>0; 力.图4是一个4×4二维网络型模.图中P表示对 9= 0 P:-P+P.g≤0. 应的孔隙压力，G和G分别表示空隙间的水平方 (14) 向和垂直方向的导流系数 式中，P,指的是孔隙i与孔隙j之间的毛管压力. 2 4 每经过一个时间步t,随着驱替相进入孔隙， 模型内的流动阻力发生变化，孔隙间的压力发生变 化，孔隙i与相邻喉道连接的孔隙广之间的平均黏度 以可表示为 2 u=aSg +u.(1-S） (15) G G 式中，u.是湿相黏度，u是非湿相黏度，S,是孔隙内 非湿相的饱和度. 9 G G 孔隙中湿相饱和度随着时间步t变化，任意一 个孔隙i在第n+1个时间步的湿相饱和度可表 示为 S=+9” 图4二维网络模型示意图 (16) V Fig.4 Schematic diagram of the 2-dimensional network model 模型的时间步并不是一个常量，每经过一个时 假设网络模型大小n×n,入口边界压力和出口 间步只有一个孔隙被驱替相填满，最先被填满的孔 边界压力已知，入口压力P1为第一列孔隙，即P:= 隙所需要的时间就是tm: P=…=Pm=P1,出口压力P2为第n列孔隙，即 (min =min (17) Pm=P2n=…=Pm=P2,对于模型中任意一个孔隙P (这里用压力表示)，与其相邻的孔隙压力及导流系 当网络模型只有单相流动时，根据网络模型 数如图5所示. 模拟得到的流量Q可以计算多孔介质的绝对渗 单相流动时，由质量守恒定律可知： 透率K: (Pj-P GH-+(P-P G+

北 京 科 技 大 学 学 报 第 36 卷 3. 4 网络模型流动机理 假定流体是可不压缩的，网络模型结构如图 1 所示，孔隙节点数为 n × n，相邻两个孔隙中心的距 离为 L，孔隙喉道半径的大小满足威布尔分布，网络 模型充满非湿相和湿相，非湿相占据中心位置，湿相 占据角隅，忽略孔隙内的毛管压力，考虑喉道内的毛 管力作用，那么根据质量守恒原理，注人孔隙的流量 之和应等于流出孔隙的流量，即 ∑ c j = 1 qij = 0． ( 12) 式中: qij指的是孔隙 i 与相邻喉道连接的孔隙 j 之间 的流量; c 是与孔隙 i 连接的喉道个数，即孔隙的配 位数． 若相邻孔隙 i 和 j 之间的喉道内只有单相流动 时，由导流系数的定义可知，通过该喉道的流量可表 示为 qij = Gij( pi － pj ) ． ( 13) 湿相驱替非湿相时，若孔隙 i 被驱替相充满，那 么连接孔隙 i 和 j 的喉道存在油水界面，此时通过该 喉道的流量可表示为 qij = Gij( pi － pj + pc，ij) ， pi － pj + pc，ij ＞ 0; 0， pi － pj + p { c，ij≤0． ( 14) 式中，pc，ij指的是孔隙 i 与孔隙 j 之间的毛管压力． 每经过一个时间步 tmin，随着驱替相进入孔隙， 模型内的流动阻力发生变化，孔隙间的压力发生变 化，孔隙 i 与相邻喉道连接的孔隙 j 之间的平均黏度 μ 可表示为 μ = μnw Sij + μw ( 1 － Sij) ． ( 15) 式中，μw 是湿相黏度，μnw是非湿相黏度，Sij是孔隙内 非湿相的饱和度． 孔隙中湿相饱和度随着时间步 tmin变化，任意一 个孔隙 i 在第 n + 1 个时间步的湿相饱和度可表 示为 Sn + 1 w，i = Sn w，i + qn + 1 w，i ·tmin Vi ． ( 16) 模型的时间步并不是一个常量，每经过一个时 间步只有一个孔隙被驱替相填满，最先被填满的孔 隙所需要的时间就是 tmin : tmin = min { Snw，iVi qw， } i ． ( 17) 当网络模型只有单相流动时，根据网络模型 模拟得到的流量 Q 可以计算多孔介质的绝对渗 透率 K: K = nμLQ Am ( P1 － P2 ) = nμL∑ n i = 1 qi Am ( P1 － P2 ) ． ( 18) 假设湿相为水相，非湿相为油相，那么由两相流 动时的模拟结果可以得到模型出口端的水相流量 Qw 和油相流量 Qo，对比单相流动时的流量 Q 可以 得到模型的相对渗透率． 水相相对渗透率: Krw = Qw Q ． ( 19) 油相相对渗透率: Kro = Qo Q ． ( 20) 式中: Sn + 1 w，i 是孔隙 i 的含水饱和度; qn + 1 w，i 是第 n + 1 个 时间步时水相进入孔隙 i 的流量，m3 ; Vi 是孔隙 i 的 体积，m3 ; Am 是网络模型的截面积，m2 ; qi 是出口端 孔隙 i 的流量，m3 ; P1 和 P2 分别是模型入口端压力 和出口端压力，Pa． 3. 5 网络模型压力计算 研究网络模型最重要的过程在于求解孔隙的压 力． 图4 是一个4 × 4 二维网络型模． 图中 pij表示对 应的孔隙压力，GH i，j和 GV i，j分别表示空隙间的水平方 向和垂直方向的导流系数． 图 4 二维网络模型示意图 Fig． 4 Schematic diagram of the 2-dimensional network model 假设网络模型大小 n × n，入口边界压力和出口 边界压力已知，入口压力 p1 为第一列孔隙，即 p11 = p21 = … = pn1 = p1，出口压力 p2 为第 n 列孔隙，即 p1n = p2n = … = pnn = p2，对于模型中任意一个孔隙 pij ( 这里用压力表示) ，与其相邻的孔隙压力及导流系 数如图 5 所示． 单相流动时，由质量守恒定律可知: ( Pi，j － Pi，j － 1 ) GH i，j － 1 + ( Pi，j － Pi，j + 1 ) GH i，j + · 841 ·

第2期 王小锋等：考虑固液分子作用的多孔介质动态网络模型 ·149· P-t 力为阻力，此时孔隙节点的压力局部增大；对于亲水 的孔隙，毛管压力为动力，孔隙节点的压力局部 下降 G Pux, 10 9 图5孔隙喉道流动平衡示意图 8 Fig.5 Schematic diagram of the flow equilibrium of pores and throats 7 6 (P-P-）G-+(P-P+1）G=0. (21) 联立需要求解压力的孔隙节点的质量守恒方 程，可以得到一个大型的线性方程组，这个方程组用 15 10 20 矩阵形式可以表示为 图6水驱时孔隙的压力变化 Ap =b. (22) Fig.6 Pressure changes of pores in water flooding 矩阵A与水平方向及垂直方向的导流系数有 关，常数b与入口和出口压力相关 4.3水驱时剩余油分布 (1)喉道半径对剩余油分布的影响.由图7可 4网络模型模拟实例 知，入口端和出口端如图中箭头所示，图中色阶代表 4.1模拟参数 孔隙的含油饱和度，当含水率达到98%时，平均喉 网络模型模拟过程所需要的参数如表1所示. 道半径为10μm的网络模型，含剩余油的孔隙比例 表1网络模型模拟参数取值 只有18%，平均喉道半径为5μm的网络模型，含剩 Table 1 Parameters of the network model 余油的孔隙比例多达44%.这是因为平均喉道越大 参数 取值 的网络模型，其绝对渗透率要大，孔隙间流通性越 网络节点 20×20 强，不易形成剩余油. 孔隙半径/μm 15.8-210.8 (2)润湿性对剩余油分布的影响.由图8可 喉道半径/μm 1.819.5 知，入口端和出口端如图中箭头所示，图中色阶代表 喉道长度/μm 孔隙的含油饱和度，当含水率达到98%时，油水润 80~160 孔喉比 湿比例为3：1的网络模型，含剩余油的孔隙比例多 1.5-20.8 达41%，油水润湿比例为1：1的网络模型，含剩余油 配位数 2-4 的孔隙比例只有27%.这是因为孔隙亲油时，毛管 分布函数参数 m,3~4:B,3~5 力为阻力，油湿孔隙比例越大，那么流动阻力越大， 油水界面张力/(mNm) 20 孔隙间流动速度越小，容易形成剩余油. 润湿接触角/(~) 30-160 (3)分子力对剩余油分布的影响.由图9可 油相黏度/(mPas) 20 知，入口端和出口端如图中箭头所示，图中色阶代表 水相黏度/(mPa·s) 1 孔隙的含油饱和度，当含水率达到98%时，考虑岩 入口压力/MPa 1 出口压力MPa 石壁面与流体分子作用影响的网络模型，含剩余油 0.1 的孔隙比例只有46%，不考虑岩石壁面与流体分子 相邻孔隙中心长度/μm 200 作用影响的网络模型，含剩余油的孔隙比例只有 4.2水驱时的动态压力分布 33%.这是因为考虑岩石壁面与流体分子作用影响 图6模拟的是水驱时网络模型孔隙的压力变 的网络模型，流体的黏滞力越大，孔隙间流动速度越 化，网络模型的入口端孔隙压力为1MPa,出口端孔 小，越容易形成剩余油 隙压力为0.1Pa.由图6可以看出，在网络模型的 4.4水驱时的相渗曲线 入口端，压力的波动比较平缓，而随着驱替的深入， (1)喉道半径对相渗曲线的影响.由图10可 孔隙节点的压力会发生较大的波动，这是由于两相 知，平均喉道半径越大，残余油饱和度越小，等渗点 流动时毛管压力的影响.对于油湿的孔隙，毛管压 右移，两相共渗区变大.这是因为喉道半径比较大

第 2 期 王小锋等: 考虑固液分子作用的多孔介质动态网络模型 图 5 孔隙喉道流动平衡示意图 Fig． 5 Schematic diagram of the flow equilibrium of pores and throats ( Pi，j － Pi － 1，j ) GV i － 1，j + ( Pi，j － Pi + 1，j ) GV i，j = 0． ( 21) 联立需要求解压力的孔隙节点的质量守恒方 程，可以得到一个大型的线性方程组，这个方程组用 矩阵形式可以表示为 Ap = b． ( 22) 矩阵 A 与水平方向及垂直方向的导流系数有 关，常数 b 与入口和出口压力相关． 4 网络模型模拟实例 4. 1 模拟参数 网络模型模拟过程所需要的参数如表 1 所示． 表 1 网络模型模拟参数取值 Table 1 Parameters of the network model 参数 取值 网络节点 20 × 20 孔隙半径/μm 15. 8 ～ 210. 8 喉道半径/μm 1. 8 ～ 19. 5 喉道长度/μm 80 ～ 160 孔喉比 1. 5 ～ 20. 8 配位数 2 ～ 4 分布函数参数 m，3 ～ 4; β，3 ～ 5 油水界面张力/( mN·m － 1 ) 20 润湿接触角/( °) 30 ～ 160 油相黏度/( mPa·s) 20 水相黏度/( mPa·s) 1 入口压力/MPa 1 出口压力/MPa 0. 1 相邻孔隙中心长度/μm 200 4. 2 水驱时的动态压力分布 图 6 模拟的是水驱时网络模型孔隙的压力变 化，网络模型的入口端孔隙压力为 1 MPa，出口端孔 隙压力为 0. 1 MPa． 由图 6 可以看出，在网络模型的 入口端，压力的波动比较平缓，而随着驱替的深入， 孔隙节点的压力会发生较大的波动，这是由于两相 流动时毛管压力的影响． 对于油湿的孔隙，毛管压 力为阻力，此时孔隙节点的压力局部增大; 对于亲水 的孔 隙，毛管压力为动力，孔隙节点的压力局部 下降． 图 6 水驱时孔隙的压力变化 Fig． 6 Pressure changes of pores in water flooding 4. 3 水驱时剩余油分布 ( 1) 喉道半径对剩余油分布的影响． 由图 7 可 知，入口端和出口端如图中箭头所示，图中色阶代表 孔隙的含油饱和度，当含水率达到 98% 时，平均喉 道半径为 10 μm 的网络模型，含剩余油的孔隙比例 只有 18% ，平均喉道半径为 5 μm 的网络模型，含剩 余油的孔隙比例多达 44% ． 这是因为平均喉道越大 的网络模型，其绝对渗透率要大，孔隙间流通性越 强，不易形成剩余油． ( 2) 润湿性对剩余油分布的影响． 由图 8 可 知，入口端和出口端如图中箭头所示，图中色阶代表 孔隙的含油饱和度，当含水率达到 98% 时，油水润 湿比例为 3∶ 1的网络模型，含剩余油的孔隙比例多 达 41% ，油水润湿比例为 1∶ 1的网络模型，含剩余油 的孔隙比例只有 27% ． 这是因为孔隙亲油时，毛管 力为阻力，油湿孔隙比例越大，那么流动阻力越大， 孔隙间流动速度越小，容易形成剩余油． ( 3) 分子力对剩余油分布的影响． 由图 9 可 知，入口端和出口端如图中箭头所示，图中色阶代表 孔隙的含油饱和度，当含水率达到 98% 时，考虑岩 石壁面与流体分子作用影响的网络模型，含剩余油 的孔隙比例只有 46% ，不考虑岩石壁面与流体分子 作用影响的网络模型，含剩余油的孔隙比例只有 33% ． 这是因为考虑岩石壁面与流体分子作用影响 的网络模型，流体的黏滞力越大，孔隙间流动速度越 小，越容易形成剩余油． 4. 4 水驱时的相渗曲线 ( 1) 喉道半径对相渗曲线的影响． 由图 10 可 知，平均喉道半径越大，残余油饱和度越小，等渗点 右移，两相共渗区变大． 这是因为喉道半径比较大 · 941 ·