孫致學(xué)，姚 軍，樊冬艷，王月英，張礦生

(1.中國(guó)石油大學(xué)石油工程學(xué)院，山東青島266580;2.長(zhǎng)慶油田分公司油氣工藝研究院，陜西西安710021)

水平井結(jié)合水力壓裂技術(shù)已成為油田增儲(chǔ)上產(chǎn)的必要措施，特別是對(duì)于低滲透以及非常規(guī)油氣藏的開(kāi)采［1-3］。國(guó)內(nèi)外大量學(xué)者對(duì)壓裂水平井的流動(dòng)形態(tài)及非穩(wěn)態(tài)產(chǎn)能進(jìn)行了研究，其中解析解主要以實(shí)空間及Laplace空間點(diǎn)源解為出發(fā)點(diǎn)［4-9］，假定裂縫面流量均勻分布，再用附加壓降的方式建立有限導(dǎo)流壓裂水平井模型，該方法并非把基巖和裂縫真正耦合求解，僅對(duì)均質(zhì)地層的橫向或縱向人工裂縫解析方法較為實(shí)用;由于地層構(gòu)造及壓裂工藝的復(fù)雜多樣，人工裂縫并非規(guī)則延伸，而是在周?chē)嬖诜种?、傾斜和不對(duì)稱(chēng)等復(fù)雜分布的情況［10-12］，此時(shí)解析方法往往較難，求解過(guò)程繁瑣。其次，對(duì)于常用的有限差分?jǐn)?shù)值方法［13-15］，需要對(duì)人工裂縫進(jìn)行局部網(wǎng)格加密處理，導(dǎo)致計(jì)算量巨大，且結(jié)構(gòu)化網(wǎng)格無(wú)法準(zhǔn)確刻畫(huà)人工裂縫的復(fù)雜分支結(jié)構(gòu)，因此筆者針對(duì)裂縫網(wǎng)絡(luò)采用有限元非結(jié)構(gòu)化網(wǎng)格進(jìn)行剖分，且利用離散裂縫模型的思想［16-20］對(duì)人工裂縫進(jìn)行降維處理，減小數(shù)值計(jì)算量;考慮人工裂縫內(nèi)流動(dòng)，對(duì)基巖和裂縫的壓力值進(jìn)行同時(shí)耦合求解，驗(yàn)證算法的正確性并分析復(fù)雜裂縫系統(tǒng)下壓裂水平井的非穩(wěn)態(tài)流動(dòng)。

1 復(fù)雜裂縫系統(tǒng)水平井模型的建立

1.1 物理模型

由于實(shí)際地質(zhì)構(gòu)造及壓裂工藝的多樣性，水平井進(jìn)行多段壓裂改造，周?chē)霈F(xiàn)復(fù)雜的裂縫系統(tǒng)。為了更準(zhǔn)確地模擬水平井壓裂改造后的流動(dòng)形態(tài)及產(chǎn)量變化，模擬三維盒狀封閉油藏內(nèi)一口壓裂水平井，其中裂縫存在分支、交錯(cuò)、傾斜等復(fù)雜形態(tài);油藏在x、y和z方向的長(zhǎng)度分別為Xe、Ye和h，水平井長(zhǎng)度為L(zhǎng)。

假設(shè):對(duì)水平井段進(jìn)行人工壓裂，裂縫條數(shù)為Nf，在裂縫擴(kuò)展過(guò)程中可能產(chǎn)生復(fù)雜的分支小裂縫;巖石和流體均微可壓縮，且壓縮系數(shù)為常數(shù)Ct;裂縫為有限導(dǎo)流，假定地層和裂縫內(nèi)流體均服從達(dá)西滲流規(guī)律;水平井裸眼或射孔完井，若定井底流量生產(chǎn)時(shí)可分析水平井的不穩(wěn)定壓力動(dòng)態(tài)，或者定井底流壓生產(chǎn)可得到水平井的動(dòng)態(tài)產(chǎn)能;油藏等溫滲流，油層較薄不考慮重力的影響。

1.2 數(shù)學(xué)模型

基于離散裂縫模型的思想，分別建立流體在基巖和裂縫系統(tǒng)的數(shù)學(xué)模型，由于均服從達(dá)西滲流規(guī)律，故基巖和裂縫系統(tǒng)的控制方程可統(tǒng)一寫(xiě)為

其中

式中，l=m，f分別代表基巖系統(tǒng)和裂縫系統(tǒng);φl(shuí)為基巖或裂縫系統(tǒng)的孔隙度，當(dāng)裂縫系統(tǒng)無(wú)填充時(shí)φf(shuō)=1，否則0＜φf(shuō)＜1，因此該模型可以考慮人工裂縫被支撐劑充填的影響;Ct為綜合壓縮系數(shù)，10 MPa－1;pl為基巖或裂縫系統(tǒng)的壓力值，10-1MPa;Kl為基巖系統(tǒng)或人工裂縫系統(tǒng)滲透率張量，μm2;ql為點(diǎn)源單位體積流量，s－1;δ(M-M′)為 Delta函數(shù)，M=M′時(shí)等于1，其他情況等于0。

定義無(wú)因次參數(shù)為

式中，kc為參考滲透率值，本文取 kc=kxx，m，即基巖在x方向的主滲透率值。

定井底流量生產(chǎn)時(shí)，無(wú)因次壓力和流量的定義為

定水平井井底流壓生產(chǎn)時(shí)，無(wú)因次流量和壓力的定義為

化簡(jiǎn)得到基巖系統(tǒng)的無(wú)因次數(shù)學(xué)模型為

人工裂縫系統(tǒng)的無(wú)因次數(shù)學(xué)模型為

式中，KfD為人工裂縫內(nèi)的無(wú)因次滲透率值;φf(shuō)D=φf(shuō)/φm為無(wú)因次裂縫孔隙度，即裂縫和基巖孔隙度比值;qfD為人工裂縫流入水平井筒的流量值。

2 復(fù)雜裂縫系統(tǒng)水平井模型的求解

在巖石力學(xué)中通常把裂縫描述為光滑平行板模型，裂縫中流體符合Navier-Stokes方程，當(dāng)流動(dòng)為層流時(shí)，可把其流量寫(xiě)成等效的達(dá)西定律形式，即著名的立方定律［21］，由此得到裂縫的等效滲透率以及沿裂縫開(kāi)度的滲流速度，此時(shí)沿壓裂縫開(kāi)度方向流動(dòng)參數(shù)及相關(guān)的物理量為常數(shù)，因此可對(duì)裂縫進(jìn)行降維處理，建立離散裂縫模型。在三維問(wèn)題中，裂縫簡(jiǎn)化為二維的面單元，如圖1所示。為了保證積分值相等，須在面積分前乘以裂縫開(kāi)度，該方法大大減少了網(wǎng)格剖分?jǐn)?shù)目，提高了數(shù)值計(jì)算效率。

圖1 離散裂縫模型示意圖Fig.1 Schematic of discrete-fracture model

此時(shí)，油藏Ω分為基巖Ωm和裂縫Ωf兩部分，用F代表油藏的流動(dòng)方程，則其積分形式可表示為

對(duì)于離散裂縫模型，其相應(yīng)的積分形式簡(jiǎn)化為

式中，af為無(wú)因次裂縫開(kāi)度值。

首先，對(duì)基巖系統(tǒng)進(jìn)行單元特性分析，采用三維四面體單元進(jìn)行剖分，并選取一次單元形函數(shù)，則基巖系統(tǒng)任一點(diǎn)的壓力值可近似為

其中

式中，Ne，m為單元的形函數(shù);Pe，m為該單元節(jié)點(diǎn)處的壓力值。

對(duì)方程(2)中的控制方程單元積分，并利用高斯分部積分公式可得到四面體單元的特性矩陣為

其中

同理，對(duì)復(fù)雜裂縫系統(tǒng)采用二維三角形單元進(jìn)行剖分，可得裂縫的單元特性矩陣為

其中

把基巖和人工裂縫單元矩陣和列陣進(jìn)行組合得到油藏的整體矩陣和列陣，如圖2所示。油藏結(jié)點(diǎn)總數(shù)為 Np，P= ［p1，p2，……，pNp］T，為了方便起見(jiàn)，油藏的整體方程組記為

圖2 基巖和裂縫單元組合示意圖Fig.2 Schematic of implementation for matrix and fracture element

系數(shù)矩陣包括基巖和裂縫兩部分，即

在時(shí)間上對(duì)方程(9)采用隱式向后差分格式，則得到相應(yīng)的有限元法支配方程為

若水平井定產(chǎn)量生產(chǎn)，則方程組(10)右邊項(xiàng)已知，直接求解方程組得到水平井井底的壓力值。若定井底流壓生產(chǎn)，則在井點(diǎn)處采用Peaceman方程［22］對(duì)產(chǎn)量進(jìn)行處理，首先得到各點(diǎn)的壓力值，再代入方程求解得到井底對(duì)應(yīng)各節(jié)點(diǎn)處的流量值。

3 模型的驗(yàn)證及動(dòng)態(tài)分析

3.1 模型的正確性驗(yàn)證

為了驗(yàn)證新模型有限元求解的正確性，對(duì)常規(guī)的壓裂水平井模型進(jìn)行求解，如圖3所示。在此分別計(jì)算了水平井壓裂前和壓裂后的流動(dòng)形態(tài)，其中盒狀封閉油藏水平井和人工裂縫的無(wú)因次參數(shù)為:油藏XeD=12，YeD=12，hD=0.1，水平井長(zhǎng)度LD=1，三條裂縫在x方向無(wú)因次坐標(biāo)(-0.4，0，0.4)，無(wú)因次裂縫半長(zhǎng)xfD=0.05，裂縫滲透率kf/km=1000，無(wú)因次裂縫開(kāi)度af=10-5。

圖3 常規(guī)壓裂水平井示意圖Fig.3 Schematic of conventional fractured horizontal well

首先，水平井壓裂前無(wú)需考慮人工裂縫，此時(shí)與經(jīng)典的Odeh水平井模型［23］解析解進(jìn)行對(duì)比，設(shè)置新模型中人工裂縫滲透率為0，計(jì)算結(jié)果如圖4所示，可見(jiàn)新模型的有限元數(shù)值解與解析解結(jié)果一致，一方面說(shuō)明了編制的有限元數(shù)值算法正確，同時(shí)新模型也適用于未壓裂水平井的模擬計(jì)算。

其次，對(duì)水平井進(jìn)行人工壓裂，假設(shè)只有裂縫向水平井井筒供液，與Zerzar［5］利用Laplace空間點(diǎn)源得到裂縫流量均勻分布時(shí)壓裂水平井解析解進(jìn)行對(duì)比，其中新模型假定裂縫內(nèi)流體服從達(dá)西滲流規(guī)律，此時(shí)得到的壓裂水平井流動(dòng)形態(tài)對(duì)比如圖5所示，由圖5可見(jiàn)兩模型的流動(dòng)形態(tài)基本一致，都存在裂縫線(xiàn)性流、裂縫徑向流、地層線(xiàn)性流及地層徑向流4個(gè)流動(dòng)階段，但由于新模型中裂縫內(nèi)部為達(dá)西滲流，所以存在裂縫內(nèi)部的徑向流動(dòng)階段，因此在新模型下壓裂水平井的流動(dòng)過(guò)程包括邊界影響為6個(gè)流動(dòng)階段:A——裂縫內(nèi)徑向流、B——裂縫線(xiàn)性流、C——裂縫徑向流、D——地層線(xiàn)性流、E——地層徑向流、F——到達(dá)邊界的擬穩(wěn)態(tài)。

圖4 本文模型與Odeh模型對(duì)比Fig.4 Comparison of presented model and Odeh model

圖5 本文模型與Zerzar模型的壓裂水平井流動(dòng)形態(tài)對(duì)比Fig.5 Comparison of presented model and Zerzar model for fractued horizontal well flowing performance

3.2 復(fù)雜裂縫系統(tǒng)下水平井動(dòng)態(tài)分析

依據(jù)實(shí)際地質(zhì)構(gòu)造及分支裂縫延伸情況，本文中構(gòu)造了一口水平井周?chē)嬖趶?fù)雜的分支裂縫，該模型更符合實(shí)際油藏中水平井壓裂后的裂縫展布，其剖面如圖6所示，各分支裂縫由2階貝塞爾曲線(xiàn)在z方向上延伸得到，由于在z方向裂縫全部穿透地層，故只給出各節(jié)點(diǎn)平面上的無(wú)因次坐標(biāo)值(表1)。若能通過(guò)微地震等現(xiàn)代檢測(cè)技術(shù)得到水平井周?chē)鷮?shí)際的裂縫分布，也可以構(gòu)造相應(yīng)的水平井裂縫系統(tǒng)，并對(duì)其進(jìn)行壓力及產(chǎn)能動(dòng)態(tài)分析。

圖6 復(fù)雜裂縫系統(tǒng)下水平井剖面圖Fig.6 Profile map of horizontal well with complex fracture system

表1 各裂縫分支的無(wú)因次平面坐標(biāo)Table 1 Dimensionless plane coordinates of fractures

3.2.1 壓力動(dòng)態(tài)分析

若假定水平井定產(chǎn)量生產(chǎn)，水平井射孔完井，流體只能通過(guò)裂縫流入水平井井筒，且裂縫和地層均服從達(dá)西滲流規(guī)律。盒狀封閉油藏水平井和人工裂縫的無(wú)因次參數(shù)為:油藏 XeD=12，YeD=12，hD=0.1，水平井長(zhǎng)度LD=1，3條裂縫在x方向無(wú)因次坐標(biāo)(－0.4，0，0.4)，無(wú)因次裂縫半長(zhǎng) xfD=0.05，裂縫滲透率kf/km=1000，無(wú)因次裂縫開(kāi)度af=10－5?；陔x散裂縫模型把裂縫簡(jiǎn)化為面單元，水平井為線(xiàn)單元，利用四面體單元和三角形單元對(duì)模型進(jìn)行剖分，且為了在較小的時(shí)間步長(zhǎng)下結(jié)果更準(zhǔn)確，對(duì)裂縫周?chē)M(jìn)行加密處理，得到壓裂水平井周?chē)木W(wǎng)格剖分圖，如圖7所示。

為了方便對(duì)比，把復(fù)雜裂縫系統(tǒng)下水平井和常規(guī)壓裂水平井的無(wú)因次流量及壓力動(dòng)態(tài)進(jìn)行對(duì)比得到圖8、9。由圖8可見(jiàn)，隨裂縫分支的增加，外部裂縫所占流量的比例增大，內(nèi)部裂縫的流量比例相對(duì)減小，總體來(lái)說(shuō)相對(duì)流量的變化并不是很大;但由圖9可見(jiàn)復(fù)雜裂縫壓力值一直小于常規(guī)壓裂水平井，由于裂縫分支增加了地層的改造體積，減少了滲流阻力，因此復(fù)雜裂縫條件下壓裂水平井總產(chǎn)量增大;其次相對(duì)于常規(guī)壓裂水平井的流動(dòng)形態(tài)，復(fù)雜裂縫情況下裂縫內(nèi)的徑向流動(dòng)和裂縫線(xiàn)性流不再出現(xiàn)，而是很快過(guò)渡到裂縫周?chē)膹较蛄鲃?dòng)，到達(dá)地層線(xiàn)性流后，與常規(guī)壓裂水平井壓力及壓力導(dǎo)數(shù)曲線(xiàn)逐漸靠近，且同時(shí)達(dá)到地層徑向流動(dòng)階段，說(shuō)明壓裂水平井裂縫擴(kuò)展的復(fù)雜性，對(duì)壓裂水平井晚期的流動(dòng)影響不大，早期影響比較明顯。

圖7 復(fù)雜裂縫系統(tǒng)下水平井網(wǎng)格剖分Fig.7 Meshing diagram of horizontal well with complex fracture system

圖8 水平井各裂縫無(wú)因次流量對(duì)比Fig.8 Dimensionless flux distribution of horizontal well fractures

圖9 復(fù)雜與常規(guī)壓裂水平井流動(dòng)形態(tài)對(duì)比Fig.9 Inflow performance comparison of complex and conventional fractured horizontal well

3.2.2 不穩(wěn)定產(chǎn)能

為了說(shuō)明分支裂縫對(duì)水平井壓裂產(chǎn)量的影響，假定復(fù)雜裂縫系統(tǒng)下水平井定井底流壓生產(chǎn)，水平井射孔完井，流體只能通過(guò)裂縫流入水平井井筒，且裂縫和地層均服從達(dá)西滲流規(guī)律，得到復(fù)雜和常規(guī)壓裂水平井的產(chǎn)量(圖10);此時(shí)油藏及裂縫的基本參數(shù)為:油藏尺寸600 m×600 m×5 m，水平井長(zhǎng)300 m，利用無(wú)因次定義可得到分支裂縫參數(shù)，人工裂縫開(kāi)度3 mm，地層初始?jí)毫?4.6 MPa，井底流壓10 MPa，地層平均滲透率 0.1 ×10－3μm2，人工裂縫內(nèi)的滲透率為10 μm2。

圖10 復(fù)雜與常規(guī)壓裂水平井日產(chǎn)量和累積產(chǎn)量對(duì)比Fig.10 Productivity and cumulative yield comparison of complex and conventional fractured horizontal well

由圖10可見(jiàn)，壓裂水平井日產(chǎn)量隨生產(chǎn)時(shí)間逐漸減小，同時(shí)復(fù)雜裂縫系統(tǒng)下產(chǎn)量明顯高于常規(guī)壓裂水平井的產(chǎn)量，主要由于裂縫分支擴(kuò)大了儲(chǔ)層改造體積(SRV)，減小了流體流入井筒的滲流阻力，大大增加了壓裂水平井的產(chǎn)量值;常規(guī)壓裂水平井生產(chǎn)三個(gè)月后的產(chǎn)量大約有4 m3/d，而復(fù)雜裂縫系統(tǒng)下水平井產(chǎn)量為8 m3/d，可見(jiàn)分支裂縫不但增加了油井產(chǎn)量也延緩了水平井產(chǎn)量遞減的速度，水平井壓裂過(guò)程中對(duì)裂縫分支的研究非常必要。

4 結(jié)論

(1)與常規(guī)壓裂水平井Zerzar模型解析解對(duì)比結(jié)果說(shuō)明了數(shù)值算法的正確性，同時(shí)指出壓裂水平井到達(dá)邊界前，除了裂縫周?chē)€(xiàn)性流、裂縫周?chē)鷱较蛄?、地層線(xiàn)性流和地層徑向流動(dòng)之外，還可能在最早期出現(xiàn)裂縫內(nèi)的徑向流動(dòng)。

(2)定井底流量生產(chǎn)時(shí)，裂縫分支對(duì)壓裂水平井晚期的流動(dòng)影響不大，早期影響比較明顯，外部裂縫產(chǎn)能所占比重相對(duì)增加;定井底流壓生產(chǎn)時(shí)，復(fù)雜裂縫大大增加了儲(chǔ)層改造體積，減小了滲流阻力，提高了水平井產(chǎn)量，同時(shí)延緩了產(chǎn)量遞減速度。

［1］姜瑞忠，蔣廷學(xué)，汪永利.水力壓裂技術(shù)的近期發(fā)展及展望［J］.石油鉆采工藝，2004，26(4):52-56.

JIANG Rui-zhong， JIANG Ting-xue， WANG Yong-li.Present development and prospecting of hydraulic fracturing technology［J］.Oil Drilling ＆ Production Technology，2004，26(4):52-56.

［2］張衛(wèi)東，郭敏，楊延輝.頁(yè)巖氣鉆采技術(shù)綜述［J］.中外能源，2010，15(6):35-40.

ZHANG Wei-dong，GUO Min，YANG Yan-hui.Shale gas drilling technologies at a glance［J］.Sino-Global Energy，2010，15(6):35-40.

［3］唐穎，張金川，張琴，等.頁(yè)巖氣井水力壓裂技術(shù)及其應(yīng)用分析［J］.天然氣工業(yè)，2010，30(10):33-38.

TANG Ying，ZHANG Jin-chuan，ZHANG Qin，et al.An analysis of hydraulic fracturing technology in shale gas wells and its application［J］.Natural Gas Industry，2010，30(10):33-38.

［4］ZERZAR A，BETTAM Y.Interpretation of multiple hydraulically fractured horizontal wells in closed systems［R］.SPE 84888，2003.

［5］ZERZAR A，TIAB D，BETTAM Y.Interpretation of multiple hydraulically fractured horizontal wells［R］.SPE 88707，2004.

［6］孫海，姚軍，廉培慶，等.考慮基巖向井筒供液的壓裂水平井非穩(wěn)態(tài)模型［J］.石油學(xué)報(bào)，2012，33(1):117-122.

SUN Hai，YAO Jun，LIAN Pei-qing，et al.A transient reservoir/wellbore coupling model for fractured horizontal wells with consideration of fluid from base rocks into wellbores［J］.Acta Petrolei Sinica，2012，33(1):117-122.

［7］ LIN Jia-jing，ZHU Ding.Predicting well performance in complex fracture systems by slab source method［R］.SPE 151960，2012.

［8］RBEAWI S A，TIAB D.Transient pressure analysis of a horizontal well with multiple inclined hydraulic fractures using type-curve matching［R］.SPE 149902，2012.

［9］樊冬艷，姚軍，王子勝，等.基于不同傾角的壓裂水平井試井解釋?zhuān)跩］.水動(dòng)力學(xué)研究與進(jìn)展，2009，24(6):705-712.

FAN Dong-yan，YAO Jun，WANG Zi-sheng，et al.Well testing on fractured horizontal well with different dip angles［J］.Chinese Journal of Hydrodynamics，2009，24(6):705-712.

［10］張廣清，陳勉.水平井水力裂縫非平面擴(kuò)展研究［J］.石油學(xué)報(bào)，2005，26(3):95-97.

ZHANG Guang-qing，CHEN Mian.Non-planar propagation of hydraulica fracture near horizontal wellbole［J］.Acta Petrolei Sinica，2005，26(3):95-97.

［11］翁定為，雷群，胥云，等.縫網(wǎng)壓裂技術(shù)及其現(xiàn)場(chǎng)應(yīng)用［J］.石油學(xué)報(bào)，2011，32(2):280-284.

WENG Ding-wei，LEI Qun，XU Yun，et al.Network fracturing techniques and its application in the field［J］.Acta Petrolei Sinica，2011，32(2):280-284.

［12］陳作，薛承謹(jǐn)，蔣廷學(xué)，等.頁(yè)巖氣井體積壓裂技術(shù)在我國(guó)的應(yīng)用建議［J］.天然氣工業(yè)，2010，30(10):30-32.

CHEN Zuo，XUE Cheng-jin，JIANG Ting-xue，et al.Proposal for the application of fracturing by stimualted reservoir volume(SRV)in shale gas wells in China［J］.Natural Gas Industry，2010，30(10):30-32.

［13］HEGRE T M.Hydraulically fractured horizontal well simulation［R］.SPE 35506，1996.

［14］MAVADDAT M，SOLEIMANI A，RASAEI M R，et al.Well test analysis of multiple hydraulically fractured horizontal wells(MHFHW)in gas condensate reservoirs［R］.SPE 144345，2011.

［15］RASHEED O B.Rate transient analysis and investigation of unit slopes in multistage hydraulicalically fractured tight carbonate reservoirs［R］.SPE 147798，2011.

［16］MOINFAR A，WAYNE N，HUI M H，et al.Comparison of discrete-fracture and dual-permeability models for multiphase flow in naturally fractured reservoirs［R］.SPE 142295，2011.

［17］姚軍，黃朝琴，王子勝，等.縫洞型油藏的離散縫洞網(wǎng)絡(luò)流動(dòng)數(shù)學(xué)模型［J］.石油學(xué)報(bào)，2010，31(5):815-824.

YAO Jun，HUANG Zhao-qin，WANG Zi-sheng，et al.Mathematical model of fluid flow in fractured vuggy reservoirs based on discrete fracture-vug network［J］.Acta Petrolei Sinica，2010，31(5):815-824.

［18］ZENG F，ZHAO G，LIU H.A new model for reservoirs with a discrete-fracture system［J］.Journal of Canadian Petroleum Technology，2012，51(2):127-136.

［19］姚軍，王子勝，張?jiān)?，?天然裂縫性油藏的離散裂縫網(wǎng)絡(luò)數(shù)值模擬方法［J］.石油學(xué)報(bào)，2010，31(2):284-288.

YAO Jun，WANG Zi-sheng，ZHANG Yun，et al.Numerical simulation method of discrete fracture net work for naturally fractured reservoirs［J］.Acta Petrolei Sinica，2010，31(2):284-288.

［20］LI L Y，SEONG H L.Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media［J］.SPE Reservoir Evaluation ＆ Engineering，2008，11(4):750-758.

［21］WITHERSPOON P A，WANG J S，IWAI K，et al.Validity of cubic law for fluid flow in a deformable rock fracture［J］.Water Resources Research，1980，16(6):1016-1024.

［22］CHEN Zhangxin，HUAN Guanren，MA Yuanle.Computational methods for multiphase flows in porous media［M］.USA:Society for Industrial and Applied Mathematics，2006:445-475.

［23］ODEH A，BABU D K.Transient flow behavior of horizontal well:pressure drawdown and buildup ananlysis［R］.SPE 18802，1992.