图1物理模型示意图
Fig. 1 Schematic diagram of the physical model
我国西北塔里木盆地广泛发育碳酸盐岩缝洞型油藏，缝洞型油气藏的储集介质由溶洞、裂缝和基质组成，其中溶洞是主要的储集空间，而裂缝不但是储集空间，也是流动通道，二者的压力变化规律不同[1-2]．由于缝洞型油藏的渗透性好、产量高，以往以研究孔渗介质间窜流能力为主的研究方法不再适用，油田开发工作者更倾向于研究储层的缝-洞组构模式、溶洞发育规模等更加实用的储层参数，合理制定开发方案． 同时由于缝洞型油藏的非均质性极强，目前已经应用的静态方法如对地震、测井、录井资料的研究，都受精度或探测范围的影响而对井眼以外的储层性质知之甚少；因此，应用开发动态中录取的试井资料，研究井筒、溶洞、裂缝中的压力变化特征，探索出一种可以得到溶洞体积的试井方法具有重要的意义．
20世纪60年代以来，针对多重介质油藏渗流开展了大量的研究．最早由Barenblatt等[3]提出了裂缝中均质流体渗流的概念；在此基础上，Warren和Root[4]提出了基岩与裂缝的双重介质模型；Kazemi[5]分析了双重介质中裂缝均匀分布时的井底压力变化．以上模型都没有考虑溶洞的存在．随着碳酸盐岩油藏的开发，人们发现溶洞对于油藏的开发与生产也有巨大作用，溶洞的研究受到越来越多的重视，在此基础上形成了三重介质以及多重介质模型．1986年，Abdassah和Ershaghi[6]提出了缝洞型油藏中三重孔隙的概念，三重孔隙包括裂缝、溶洞和基质3种孔隙．Camacho-Velazquez等[7]在此基础上，通过求解三重介质渗流方程获得了井底压力及导数的双对数曲线，利用双对数曲线图详细分析了三重介质油藏井底压力变化特征．尽管研究有了很多的成果，但缝洞型油藏内部复杂的孔隙系统和流动规律仍然是研究中的难题．
自进入21世纪以来，由于在很多国家(如中国、加拿大)发现了大量的缝洞型油藏，缝洞型油藏试井分析受到重视．一些学者在三重介质模型的基础上，提出了一些新的试井模型，分析了井底压力变化的趋势．2003年，薛承谨等[8]提出了当井筒与裂缝相连通时的试井模型，分析了地层参数和边界条件对井底压力的动态影响．2004年，Wu等[9]针对三重介质油藏提出了一种描述流体流动和传输过程的模型．同年，常学军等[10]在文献[8]提出模型的基础上提出溶洞和裂缝同时和井筒相连的模型，并也对地层参数进行了敏感性分析；王子胜和姚军等[11-12]还考虑了非Newton(牛顿)液体和三重压敏介质的缝洞型油藏井底压力响应特征；李成勇等[13]提出了三重介质水平井的试井模型并进行了试井分析；Kang等[14]提出了三重介质中多相流流动的模型，并利用数值分析研究了压力的响应特征．2007年Wu等[15]在Camacho-Velazquez的基础上再次研究了三重介质模型的压力动态响应特征，求出了井底压力并运用数值反演技术得到了井底压力与时间的无量纲曲线，同时进行了敏感性分析．2007年，蔡明金等[16]提出了三重介质有限导流垂直裂缝的试井模型并进行了压力动态分析；张冬丽和李江龙等[17-19]提出了缝洞型油藏的数学模型及其研究进展，随后他们提出了三重介质油藏的数值试井模型，并进行了试井分析，研究了各个因素对试井结果的影响．2011年，Wu等[20]提出了多重介质数值模拟模型来研究缝洞型油藏中的单相流和多相流问题，取得了良好的效果．2012年，Guo等[21]研究了缝洞型油藏水平井中的流动问题．
目前对于缝洞型油藏井底压力的计算研究中，通常采用三重介质模型，将溶洞视作一种介质，采用了窜流方程来描述基质-裂缝-溶洞3种孔隙之间的流动，然而在实际测试的压力恢复曲线中，其双对数压力及导数曲线并未出现三重介质曲线特征，所以忽略洞中的流动而将其简化为窜流与生产实际并不吻合．本文通过联立流体力学中的质量守恒方程、动量守恒方程和能量守恒方程，提出了井筒中的压力是以波的形式向井底传播的，在此基础上建立了洞中波动和流动相耦合的模型，同时外部地层仍然考虑三重介质渗流方程，形成了缝洞型油藏新的控制方程组，并且得到了井底压力的双对数曲线．通过压力瞬时分析可以得到溶洞的体积，结合实例分析证明，采用这样的模型进行试井分析更符合生产实际，对于缝洞型油藏开发具有重要的指导意义．
为获得缝洞型油藏流动规律，本文采用以下假设：
1) 假设大溶洞为圆柱，并且与井筒同轴，仅考虑竖直方向的流动；
2) 地层为各向同性圆形油藏，圆心处有一定产量生产的油井；
3) 地层等厚且流体为单相，且岩石及流体物性等参数都是常数；
4) 地层是由基质、裂缝和溶洞组成，基质会向裂缝和溶洞窜流，同时裂缝和溶洞之间存在窜流；地层中的裂缝和溶洞都会向井筒和大溶洞渗流供给流体[3]；
5) 流体由储层向裂缝和地层中的溶洞渗流后再流入井筒，流体在储层中为等温渗流且符合Darcy(达西)渗流规律，忽略流体的毛管压力的影响[22]；
6) 油井测试前地层中各点的压力均匀．
基于上述假设，可将建立的物理模型简化为图1所示的理想模型．为将与井筒相连的溶洞和地层中的溶洞相区别，本文将与井筒相连的溶洞称为大溶洞，涉及该部分的公式中的参数下标取lv．
图1物理模型示意图
Fig. 1 Schematic diagram of the physical model
流体从溶洞流入井筒，再由井筒流出地面，这一过程的流体流动要满足的连续性方程、动量守恒和能量守恒方程为
(1)
(2)
(3)
式中，ρ为流体密度(kg/m3)；v为流体流动速度(m/s)；x轴为由井筒圆心向下建立的一维坐标轴；p为压力(Pa)；f为流体受到的摩擦阻力系数；D为井筒的直径(m)；pwf和plv分别为井筒和大溶洞中的压力(Pa)；vwf为井筒和大溶洞连接处流体的速度(m/s)．
按本文的假设，井筒与溶洞组成一个巨大的高压流体储集空间(如新疆缝洞型油藏，井底压力在80～100 MPa间),在高压流体作用下，储集空间的流体、管壁及岩石都是可压缩的．当开井生产时，井底压力变化由两部分组成：一部分是流体流动产生的压力；另一部分是高压流体泄压时，由于流体、管壁及岩石压缩性产生的波动压力(流体中水击效应也有类似的方程)[23],下面给出波传播压力的详细推导过程．
图2流体微元中流体流入和流出示意图
Fig. 2 Schematic diagram of the infinitesimal fluid inflow and outflow
在井筒中取流体微元如图2所示，由质量守恒可得到
(4)
式中，A为微元面积(m2)；δx为微元高度(m)．
在高压状况下，流体存在压缩性，油管也是一个弹性体，其变形由油管直径、壁厚及油管材料的弹性模量决定，方程(4)展开变为
(5)
根据流体力学中全导数和偏导数的关系：
(6)
方程(5)可修改为
(7)
由流体的可压缩性，可将方程(7)中的密度项表示成压力的函数：
(8)
式中，G为流体的体积模量(Pa)．假设油管为弹性变形，对薄壁的圆管，当压力增加dp时，其径向变形dD与dp的关系为[23]
(9)
式中，D为油管直径(m)；e为油管壁厚(m)；Ew为油管弹性模量(Pa)．由油管面积公式：
(10)
联合方程(8)和(9)，则式(7)可变为
(11)
定义
(12)
式中，C为管道及流体系统中的波速(m/s)．
根据方程(6)，方程(11)可变为
(13)
可见压力在x方向的传导是以波的形式进行的，波速为C．
将连续性方程和动量守恒方程联合，可得
(14)
考察流体在大溶洞中的流动，由于速度v较小，而摩擦力为v的二次方项，可略去二阶小量，则有
(15)
方程(15)的解即为溶洞中的流体流速：
(16)
式中，v0为初始时刻的速度(m/s)，可由地面产量确定．方程(16)中重力值出现在速度表达式中的常数项部分，考虑到大溶洞的半径rlv很大，而
更大，常数项应该是小量，为方便试井分析中无量纲量的定义，可将该常数项修正到附加压降中，在本文的后半部分定义了修正系数，因此方程(16)可表示为
(17)
由此得到溶洞提供的产量：
(18)
于是在井筒与大溶洞相连处流体的速度为
(19)
将方程(19)代入方程(3)，可得如下表达式：
(20)
在缝洞型油藏试井模型中，流动方程是由井筒、洞及地层三部分组成，对于外部地层的流动，本文仍然采用渗流方程．地层中的压力变化方程可表示为
(21)
(22)
(23)
考虑井筒和大溶洞中的表皮系数[24]：
(24)
(25)
考虑定产量井，由图1可知产量q一共由三部分组成：井筒存储和大溶洞存储；裂缝向井筒渗流和裂缝向大溶洞渗流部分；溶洞向井筒和大溶洞渗流部分．可得边界条件为
(26)
其中,μ为流体的黏度(Pa·s);Cw和Clv分别为井筒及大溶洞存储常数(m3/Pa);sw和slv分别为井筒和溶洞的表皮系数；q为总产量(m3/s);B为流体体积系数．
根据附录A的无量纲定义，可以将方程(21)～(26)化为无量纲方程组：
(27)
结合前面的无量纲定义，同时考虑ρ是压力的函数，并考虑重力、摩擦力及洞为非圆柱体等，无量纲井筒压力及溶洞压力之间的关系为
(28)
式中，α为前文所述的修正系数．
如假设地层为无限大(或圆形封闭或圆形定压)，结合方程(26)和(27)求解，可得到Laplace空间上的无量纲井底压力统一表达式为
(29)
附录B中详细给出推导过程及参数定义，对Laplace空间上的井底压力使用Stehfest算法[25]进行Laplace数值反演可得到实空间上给定无量纲时间下，无量纲压力及其导数数值．
本文的重点是研究缝洞型油藏井底压力计算，计算实例中，地层也只考虑无限大或圆形封闭或定压，并且kD取1．
无量纲参数β，γ和α分别是洞体积、波速、洞存储系数及洞中流动阻力等组合而形成的无量纲参数，这些参数取不同值会影响井底压力及导数曲线的形态．从β及γ的表达式可以看出，γ不仅与洞体积有关，还与波速有关，它对井底压力的影响是按时间指数式衰减，本文将γ命名为波动系数．β主要与洞体积及井筒半径有关，由于流体从洞流入井筒产生节流效应，对井底压力的贡献相当于阻力，本文将β命名为阻尼系数．随着时间增加，阻力效应会越来越小，所以将α命名为衰减指数．β，γ，α对典型曲线都有贡献．
采用无限大地层的井底压力表达式，可得到井底压力及其导数与时间的双对数曲线如图3所示．上面一条曲线为压力曲线，下面一条为压力导数曲线，下文中曲线均是按照这个规则排列．无限大地层中井与溶洞相连，井底压力及其导数双对数典型曲线可以分为以下6个特征段：
1) Ⅰ区为井筒存储段．此时井口产量主要由井筒存储供给，是典型的井筒存储段，所以压力及其导数双对数曲线均为斜率为1的直线段．
2) Ⅱ区为井筒存储向溶洞存储的过渡段．随着流体的采出，井筒存储逐渐向溶洞存储过渡，由于溶洞与井筒间流动通道突然变小，产生附加压降(类似于表皮效应)，所以压力导数开始下降，出现了第一个下凹段．
3) Ⅲ区为溶洞存储段，此时井筒存储完全结束，进入溶洞存储阶段，在压力导数图中再次出现一条斜率为1的直线段，随着流体不断采出，地层中的流体开始向溶洞补充流量，表皮效应开始显现．
4) Ⅳ、Ⅴ和Ⅵ区是地层渗流段，Ⅳ区为溶洞向裂缝窜流段，这一部分一开始为裂缝向井筒和大溶洞供给流体，此时裂缝系统中的压力下降，而基质和溶洞系统中的压力不变，这样在溶洞-裂缝系统和基质-裂缝系统之间形成了压差．当压差较小时，只有溶洞向裂缝窜流，导数曲线出现了第二个下凹．
5) Ⅴ区为基质向缝洞总系统窜流段，这是由于当压力差进一步增加时，基质向缝、洞的总系统窜流，导数曲线出现了第三个下凹．
6) Ⅵ区为径向流段，压力导数曲线变平缓，当流动到达径向流后，压力导数曲线为一条水平线．
图3井底压力与压力导数双对数曲线图
Fig. 3 The log-log graph of the wellbore pressure and the pressure derivative
为了更好地分析缝洞型典型曲线的特征，本文全面分析各个参数对曲线的影响．
2.2.1 阻尼系数β对曲线的影响分析
β选取不同值时压力及其导数双对数曲线如图4所示．
从图中可以发现：当时间较小时，压力及导数曲线只受井筒存储影响，β对井储(Ⅰ区)没有影响，只对曲线的其他区域产生影响．β越大曲线的第一个下凹越深；同时还会对流动进入Ⅲ区的时间产生影响，β越大越晚进入Ⅲ区溶洞存储段；β的取值还会影响Ⅲ区的最大值，β越大，Ⅲ区的最大值越小．β的取值会使导数曲线在Ⅳ、Ⅴ、Ⅵ区产生上下的偏移，但对曲线的形态没有影响，β越大曲线越往下偏移．此处将曲线上下平移后如图5所示．可以发现β取不同值时，导数曲线在Ⅳ、Ⅴ、Ⅵ区可以完全重合，这说明β对Ⅳ、Ⅴ、Ⅵ曲线的形态没有影响．因为β是一个与溶洞体积有关的参数，所以会影响到Ⅱ区(井筒存储向溶洞存储过渡段)、Ⅲ区(溶洞存储段)；因为大溶洞的存在，所以也会对Ⅳ、Ⅴ和Ⅵ区(地层渗流段)产生影响，使得出现径向流后其压力导数并非趋近0.5的水平线．
图4β取不同值时井底压力与压力导数双对数曲线图
Fig. 4 The log-log graph of the wellbore pressure and the pressure derivative for different values ofβ
图5β取不同值时曲线平移后井底压力与压力导数双对数曲线图
Fig. 5 The log-log graph of the wellbore pressure and the pressure derivative for different values ofβ after translation
2.2.2 波动系数γ对曲线的影响分析
γ取不同值时压力及其导数与时间的双对数曲线如图6所示．
从图6的压力及导数与时间双对数曲线图上可以发现：γ对井筒存储段(Ⅰ区)的导数曲线形态没有影响，但会影响Ⅰ区结束的时间，γ越大，Ⅰ区越晚结束；同时γ对其他区也会产生影响．γ越大，压力导数曲线第一个下凹段越浅；γ取值对Ⅲ区的影响体现在进入Ⅲ区的时间上，γ越大，越晚进入Ⅲ区流动阶段．γ的取值对Ⅳ、Ⅴ和Ⅵ区的影响同β一致，只会造成曲线的上下偏移，γ越大曲线越往上平移，但曲线的形态不变．γ是一个与溶洞的体积以及波速有关的参数，所以也会对Ⅱ区(井筒存储向溶洞存储过渡段)、Ⅲ区(溶洞存储段)的流动产生影响．因为大溶洞的存在，所以对地层渗流段也会产生影响，使得出现径向流后压力导数也不是一条为0.5的水平线．
图4和图6说明波动系数γ与阻尼系数β都共同影响压力及其导数的曲线形态，因为β和γ的无量纲定义上都有溶洞半径rlv，两者都反映溶洞的体积．利用本文的模型可以由试井方法得到溶洞的体积．当地层渗流达到径向流动以后由于洞的存在，压力导数曲线并非是一条0.5的水平线．
图6γ取不同值时井底压力与压力导数随时间变化曲线图
Fig. 6 The double logarithmic graph of the wellbore pressure and the pressure derivative for different values ofγ
2.2.3 衰减指数α对曲线的影响分析
由图7可以发现α也对曲线的Ⅰ区(井筒存储段)无影响．α越大，导数曲线的第一下凹段越浅；α的取值还会影响Ⅲ区的最大值，α越大，Ⅲ区的最大值越大．α的取值对Ⅳ、Ⅴ和Ⅵ区的影响同β和γ一样，这里不再赘述．
图7 α取不同值时井底压力与压力导数随时间变化曲线图
Fig. 7 The double logarithmic graph of the wellbore pressure and the pressure derivative for different values ofα
2.2.4 边界条件对曲线的影响分析
不同边界条件下的井底压力及压力导数曲线图如图8所示，可以发现，边界条件的不同仅仅影响曲线的Ⅳ、Ⅴ、Ⅵ区(地层渗流段)．不同边界条件下的地层渗流段曲线形态不同，当是无限大地层边界时，曲线末段为一条水平线；当是圆形封闭边界时，曲线末段上扬；当是定压边界时，曲线末段下落．
图8不同边界条件时井底压力与压力导数随时间变化曲线图
Fig. 8 The log-log graph of the wellbore pressure and the pressure derivative with different boundary conditions
2.2.5 裂缝储容比ωf对曲线的影响分析
由图9可以发现ωf对曲线的Ⅲ区(溶洞存储段)有影响，ωf越大，导数曲线的第二个下凹段越深，且持续时间越长．随着ωf的减小，第二个下凹段越来越不明显．
图9ωf取不同值时井底压力与压力导数随时间变化曲线图
Fig. 9 The double logarithmic graph of the wellbore pressure and the pressure derivative for different values ofωf
2.2.6 溶洞储容比ωv对曲线的影响分析
由图10可以发现ωv只对曲线的Ⅳ、Ⅴ区(地层渗流段)有影响，对于Ⅳ区(溶洞向裂缝窜流段)：ωv越大导数曲线在Ⅳ区越早到达下凹处，且第二个下凹段越浅；对于Ⅴ区(基质向缝洞总系统窜流段)的影响则与其对Ⅳ区的影响效果相反．
图10ωv取不同值时井底压力与压力导数随时间变化曲线图
Fig. 10 The double logarithmic graph of the wellbore pressure and the pressure derivative for different values ofωv
2.2.7 基质裂缝窜流系数λmf对曲线的影响分析
由图11可以发现：随着λmf的减小，Ⅳ区(溶洞向裂缝窜流段)持续时间越长，同时下凹越深；并且随着λmf的减小，压力导数曲线的第三个下凹段越深．当λmf越来越小的时候第二个下凹会越来越不明显，当λmf足够小的时候会和第三个下凹合并为一个下凹．
图11λmf取不同值时井底压力与压力导数随时间变化曲线图
Fig. 11 The double logarithmic graph of the wellbore pressure and the pressure derivative for different values ofλmf
2.2.8 基质溶洞窜流系数λmv对曲线的影响分析
由图12可以发现：随着λmv的减小，Ⅴ区(基质向缝洞总系统窜流段)持续时间越长，即第三个下凹部分出现越晚，下凹部分曲线向右平移但形状不变．
图12λmv取不同值时井底压力与压力导数随时间变化曲线图
Fig. 12 The double logarithmic graph of the wellbore pressure and the pressure derivative for different values ofλmv
我国在西部地区有大量的缝洞型油藏，井深都在7 000 m以上，自开发以来，已对大量井进行了关井压力恢复测试，采用现有的商业试井软件，仅能获得渗透率、窜流系数等参数，这些参数无法与现有的地质模型相吻合．本文提出的方法可有效解决缝洞型油藏的开发难题，可由试井结果反演出洞的体积等参数．为了进一步说明缝洞型油藏井底压力曲线在试井中的应用，本文对中国新疆某井进行了试井解释．表1为该井基本参数，表2为解释结果．根据该井双对数压力及其导数曲线图，采用井筒与洞相连，地层均质双渗模型，进行试井解释：图13给出了该井压力及其导数双对数曲线拟合图．从图中可以看出：理论计算与实例数据拟合质量都非常高．
表1 实例井基本参数
Table 1 Reservoir properties of the example well
parametervalueformation thickness h/m12wellbore radius rw/m0.067porosity ϕ0.15oil viscosity μ/(Pa·s)0.002oil formation volume factor B1.1oil density ρ/(kg/m3)860compressibility of rock Ct/MPa-10.001
表2 解释结果
Table 2 Results of interpretation
parametervaluestorage constant of the wellbore Cw/(m3/MPa)3.760 4skin factor of the wellbore sw-1.62volume of the large vug V/m332 280.99storage ratio ω0.267 6inter-porosity flow coefficient λ6.03×10-12fluctuation coefficient γ0.292 2friction coefficient β16.993 1
图13实例井井底压力和压力导数拟合图
Fig. 13 Fitting curves of the wellbore pressure and the pressure derivative for the example well
本文以力学三大守恒定律为基础，提出了压力在溶洞中是以波的形式传播的，建立了溶洞中波动与流动耦合方程，给出了井筒和溶洞中的压力表达式，以此作为内边界条件对渗流方程进行了无量纲化并求解,主要结论如下：
1) 利用力学三大定律，建立了溶洞中波动和流动的耦合方程，结合外部地层渗流方程形成了新的缝洞型油藏试井模型，再通过定义无量纲参数、Laplace变换及数值反演技术，得到试井分析所需的图版，分析了各参数对典型曲线的影响．
2) 井筒与溶洞相连的典型曲线可以分为6个特征段，早期为井筒存储段，第二阶段为井筒存储向溶洞存储过渡段，第三阶段为溶洞存储段，第四、五、六阶段为地层渗流段，其中第四阶段为溶洞向裂缝窜流段，第五阶段为基质向缝洞总系统窜流段，第六阶段为径向流段．
3) 本文为缝洞型油藏试井提供了一个新的建模方法，在原试井参数基础上增加了β，γ，α这3个参数并分析了不同边界条件下的曲线形态，其中β是与初始速度、溶洞体积及密度相关的参数；γ是与溶洞体积及波速相关的参数；α则是阻力效应衰减指数．文中详细给出了这3个以及外部地层参数和不同边界条件对典型曲线的影响．
4) 本文针对各个参数进行了敏感性分析，其中波动系数、阻尼系数和衰减系数主要影响曲线的前半段，而窜流系数和储容比主要影响曲线的后半段．本文还利用新疆油田某井的压力恢复数据进行了实例分析，可以发现压力及其导数的双对数曲线均拟合较好，解释结果与实际符合．
附录A
无量纲时间
无量纲压力
无量纲半径
无量纲高度
无量纲渗透率
无量纲窜流系数
无量纲储水系数
井筒存储系数
大溶洞存储系数
文中设定的无量纲系数
m 基质 f 裂缝 v 溶洞 w 井筒 lv 大溶洞 D 无量纲
附录B
对方程(28)和(29)进行Laplace变换可得
(B1)
(B2)
(B3)
(B4)
(B5)
(B6)
(B7)
由方程(B2)得
(B8)
将其代入方程(B1)和(B2)可得
(B9)
(B10)
其中
(B11)
(B12)
(B13)
方程(B9)和(B10)的通解可设为
(B14)
(B15)
将方程(B14)、(B15)代入方程(B9)、(B10)，通过比较Bessel(贝塞尔)函数I和K各自的系数可以得到一个关于α的方程：
α4kD(1-kD)-α2[kDm2(u)+(1-kD)m1(u)]-m2(u)+m1(u)m2(u)=0．
(B16)
这个方程可以解出关于α的4个解：α1，α2，α3，α4，易知
(B17)
(B18)
由于α只能取正，所以这里的α3取作α1，α4取作α2．而且可以通过赋值的方法得到下面的关系:
C1=1时，
(B19)
C2=1时，
(B20)
E1=1时，
(B21)
E2=1时，
(B22)
可知A1=B1，A2=B2．所以可设为
(B23)
(B24)
定义：
M1=A1I0(α1),
(B25)
M2=A2I0(α2),
(B26)
M3=B1K0(α1),
(B27)
M4=B2K0(α2),
(B28)
M5=A1α1I1(α1),
(B29)
M6=A2α2I1(α2),
(B30)
M7=-B1α1K1(α1),
(B31)
M8=-B2α2K1(α2),
(B32)
N1=I0(α1rlvD),
(B33)
N2=I0(α2rlvD),
(B34)
N3=K0(α1rlvD),
(B35)
N4=K0(α2rlvD),
(B36)
N5=α1rlvDI1(α1rlvD),
(B37)
N6=α2rv2DI1(α2rlvD),
(B38)
N7=-α1rv2DK1(α1rlvD),
(B39)
N8=-α2rv2DK1(α1rlvD),
(B40)
(B41)
则一些表达式可以使用上述定义简化为
(B42)
(B43)
(B44)
(B45)
(B46)
(B47)
以下分3种边界条件分别求解．
此时边界条件为
(B48)
易知此时λ1，λ2为0．
与方程(B4)～(B7)联立可得
λ3a1+λ4a2=Y,
(B49)
(B50)
其中
a1=M1-swM5+svN5-N1,
(B51)
a2=M2-swM6+svN6-N2,
(B52)
a3=M3-swM7+svN7-N3,
(B53)
a4=M4-swM8+svN8-N4．
(B54)
解得
(B55)
(B56)
所以井底压力可表示为
(B57)
将上文中定义的各个参数代入即可．
此时边界条件为
(B58)
结合方程(B4)～(B7)可得一个矩阵方程：
(B59)
其中
b1=uCwD(M1-swM5)+uCv2D(N1-svN5)-hDswM5-(1-hD)svN5,
(B60)
b2=uCwD(M2-swM6)+uCv2D(N2-svN6)-hDswM6-(1-hD)svN6,
(B61)
b3=uCwD(M3-swM7)+uCv2D(N3-svN7)-hDswM7-(1-hD)svN7,
(B62)
b4=uCwD(M4-swM8)+uCv2D(N4-svN8)-hDswM8-(1-hD)svN8,
(B63)
c1=α1A1I1(α1reD1),
(B64)
c2=α2A2I1(α2reD1),
(B65)
c3=-α1B1K1(α1reD1),
(B66)
c4=-α2B2K1(α2reD1),
(B67)
c5=α1I1(α1reD2),
(B68)
c6=α2I1(α2reD2),
(B69)
c7=-α1K1(α1reD2),
(B70)
c8=-α2K1(α2reD2)．
(B71)
运用Cramer法则即可解出这个矩阵方程，求得λ1，λ2，λ3，λ4，代入式(B46)即可得到井底压力．
边界条件为
(B72)
结合方程(B10)～(B13)可得一个矩阵方程：
(B73)
其中
d1=A1I0(α1reD1),
(B74)
d2=A2I0(α2reD1),
(B75)
d3=B1K0(α1reD1),
(B76)
d4=B2K0(α2reD1),
(B77)
d5=I0(α1reD2),
(B78)
d6=I0(α2reD2),
(B79)
d7=K0(α1reD2),
(B80)
d8=K0(α2reD2)．
(B81)
运用Cramer法则即可解出这个矩阵方程，求得λ1，λ2，λ3，λ4，代入式(B46)即可得到井底压力．
致谢 本文作者衷心感谢中石化重大专项(ZDP7004)对本文的资助．
[1] 李阳. 塔河油田碳酸盐岩缝洞型油藏开发理论及方法[J]. 石油学报, 2013,34(1): 115-121. (LI Yang. The theory and method for development of carbonate fractured-cavity reservoirs in Tahe oilfield[J].Acta Petrolei Sinica, 2013,34(1): 115-121.(in Chinese))
[2] 邢翠巧, 尹洪军, 李兴科, 等. 缝洞型碳酸盐岩油藏试井分析模型研究[J]. 河北工业科技, 2018,35(1): 31-36.(XING Cuiqiao, YIN Hongjun, LI Xingke, et al. Research of the fractured vuggy type carbonate reservoirs well test analysis model[J].Hebei Journal of Industrial Science and Technology, 2018,35(1): 31-36.(in Chinese))
[3] BARENBLATT G, ZHELTOV I, KOCHINA I. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks(strata)[J].Journal of Applied Mathematics &Mechanics, 1960,24(5): 1286-1303.
[4] WARREN J, ROOT P. The behavior of naturally fractured reservoirs[J].Society of Petroleum Engineers Journal, 1963,3(3): 245-255.
[5] KAZEMI H. Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution[J].Society of Petroleum Engineers Journal, 1969,9(4): 451-462.
[6] ABDASSAH D, ERSHAGHI I. Triple-porosity systems for representing naturally fractured reservoirs[J].SPE Formation Evaluation, 1986,1(2): 113-127.
[7] CAMACHO-VELAZQUEZ R, VASQUEZ-CRUZ M, CASTREJON-AIVARA R, et al. Pressure transient and decline curve behaviors in naturally fractured vuggy carbonate reservoirs[J].SPE Reservoir Evaluation &Engineering, 2005,8(2): 95-112.
[8] 薛承瑾, 戴卫华, 张皎生. 裂缝与井筒连通的三重介质油藏试井解释方法[J]. 油气地质与采收率, 2003,10(3): 39-41.(XUE Chengjin, DAI Weihua, ZHANG Jiaosheng. Well test interpretation method for triple media oil reservoir in which fracture communicated with wellbore[J].Petroleum Geology and Recovery Efficiency, 2003 ,10(3): 39-41.(in Chinese))
[9] WU Y S, LIU H H, BODVARSSON G S. A triple-continuum approach for modeling flow and transport processes in fractured rock[J].Journal of Contaminant Hydrology, 2004,73(1): 145-179.
[10] 常学军, 姚军, 戴卫华, 等. 裂缝和洞与井筒连通的三重介质油藏试井解释方法研究[J]. 水动力学研究与进展, 2004,19(3): 339-346.(CHANG Xuejun, YAO Jun, DAI Weihua, et al. The study of well test interpretation method for a triple medium reservoir[J].Journal of Hydrodynamics, 2004,19(3): 339-346.(in Chinese))
[11] 王子胜, 姚军, 张凯. 三重介质油藏非牛顿液不稳定渗流压力变化特征研究[J]. 中外能源, 2006,11(2): 13-17.(WANG Zisheng, YAO Jun, ZHANG Kai. The characteristic of transient pressure of non-Newtonian fluid’s transient flow in triple-media reservoirs[J].China Foreign Energy, 2006,11(2): 13-17.(in Chinese))
[12] 王子胜, 姚军. 缝洞向井筒供液时三重压敏介质油藏压力响应特征研究[J]. 水动力学研究与进展, 2006,21(1): 84-89.(WANG Zisheng, YAO Jun. Study of pressure-transient characteristic for stress-sensitive triple-medium reservoirs with fractures and vugs conveying fluids to wellbore[J].Journal of Hydrodynamics, 2006,21(1): 84-89.(in Chinese))
[13] 李成勇, 刘启国, 张燃, 等. 三重介质油藏水平井试井解释模型研究[J]. 西南石油大学学报(自然科学版), 2006,28(4): 32-35.(LI Chengyong, LIU Qiguo, ZHANG Ran, et al. The research of well test of horizontal well in triple medium reservoir[J].Journal of Southwest Petroleum University(Science &Technology Edition), 2006,28(4): 32-35.(in Chinese))
[14] KANG Z J, WU Y S, LI J L, et al. Modeling multiphase flow in naturally fractured vuggy petroleum reservoirs[C]//Proceedings of SPE Annual Technical Conference and Exhibition. San Antonio, 2006.
[15] WU Y S, EHLIG-ECONOMIDES C, QIN G, et al. A triple-continuum pressure-transient model for a naturally fractured vuggy reservoir[C]//Proceedings of SPE Annual Technical Conference and Exhibition. Anaheim, 2007.
[16] 蔡明金, 贾永禄. 三重介质有限导流垂直裂缝井压力动态分析[J]. 油气井测试, 2007,16(5): 12-15.(CAI Mingjin, JIA Yonglu. Dynamic analysis for pressure in limit conductivity vertical fracture wells of triple-porosity reservoir[J].Well Testing, 2007,16(5): 12-15.(in Chinese))
[17] 张冬丽, 李江龙. 缝洞型油藏流体流动数学模型及应用进展[J]. 西南石油大学学报(自然科学版), 2009,31(6): 66-70.(ZHANG Dongli, LI Jianglong. Application progress of fluid flow mathematic model for fracture-vug type reservoir[J].Journal of Southwest Petroleum University(Science &Technology Edition), 2009,31(6): 66-70.(in Chinese))
[18] 张冬丽, 李江龙, 吴玉树. 缝洞型油藏三重介质数值试井模型影响因素[J]. 西南石油大学学报(自然科学版), 2010,32(2): 113-120.(ZHANG Dongli, LI Jianglong, WU Yushu. Influencing factors of the numerical well test model of the triple-continuum in fractured vuggy reservoir[J].Journal of Southwest Petroleum University(Science &Technology Edition), 2010,32(2): 113-120.(in Chinese))
[19] 李江龙, 张冬丽, 吴玉树. 缝洞型油藏三重介质数值试井解释方法研究[J]. 水动力学研究与进展, 2012,27(6): 640-648.(LI Jianglong, ZHANG Dongli, WU Yushu. Triple-continuum numerical well test interpretation method for a naturally fractured vuggy reservoir[J].Journal of Hydrodynamics, 2012,27(6): 640-648.(in Chinese))
[20] WU Y S, YUAN D, KANG Z J, et al. A multiple-continuum model for simulating single-phase and multiphase flow in naturally fractured vuggy reservoirs[J].Journal of Petroleum Science &Engineering, 2011,78(1): 13-22.
[21] GUO J C, NIE R S, JIA Y L. Dual permeability flow behavior for modeling horizontal well production in fractured-vuggy carbonate reservoirs[J].Journal of Hydrology, 2012,464/465(13): 281-293.
[22] 刘启国, 徐有杰, 刘义成, 等. 夹角断层多段压裂水平井试井求解新方法[J]. 应用数学和力学，2018,39(5): 558-567.(LIU Qiguo, XU Youjie, LIU Yicheng, et al. A new well test analysis method for multi-stage fractured horizontal wells with angle faults[J].Applied Mathematics and Mechanics, 2018,39(5): 558-567.(in Chinese))
[23] DOUGLAS J F, GASIOREK J M, SWAFFILED J A. 流体力学[M]. 北京: 世界图书出版社, 1995.(DOUGLAS J F, GASIOREK J M, SWAFFIELD J A.Fluid Dynamics[M]. Beijing: World Publishing Corporation, 1995.(in Chinese))
[24] XING C Q, YIN H J, LIU K X, et al. Well test analysis for fractured and vuggy carbonate reservoirs of well drilling in large scale cave[J].Energies, 2018,11(1): 80-94.
[25] STEHFEST H. Numerical inversion of Laplace transform-algorithm[J].Communications of the ACM, 1970,13(1): 47-49.
应用数学和力学编委会，ISSN 1000-0887
杜鑫(1995—)，男，硕士生(E-mail: duxin67@mail.ustc.edu.cn)；
卢德唐(1965—)，男，教授，博士生导师(通讯作者. E-mail: dtlu@ustc.edu.cn).
杜鑫，卢志炜，李冬梅，徐燕东，李培超，卢德唐. 缝洞型油藏波动和流动耦合模型井底压力分析[J]. 应用数学和力学, 2019,40(4): 355-374.
DU Xin, LU Zhiwei, LI Dongmei, XU Yandong, LI Peichao, LU Detang. Pressure transient analysis of the fractured vuggy reservoir model coupling oil flow and wave propagation[J].Applied Mathematics and Mechanics, 2019,40(4): 355-374.