呂愛民，王立偉，龍 濤，王傳睿，董云振，趙健男，李緒延

(1.中國(guó)石油大學(xué)(華東) 石油工程學(xué)院，山東 青島 266580;2.中石化 青島安全工程研究院，山東 青島 266000)

引 言

碳酸鹽巖縫洞型油藏主要是以構(gòu)造變形產(chǎn)生的裂縫和巖溶作用形成的孔、縫、洞為主，大裂縫大溶洞發(fā)育，具有較強(qiáng)的非均質(zhì)性和多尺度特征。依據(jù)縫洞型碳酸鹽巖油藏的開發(fā)特點(diǎn)，“縫洞單元”是指周圍被相對(duì)致密的隔檔層體遮擋,由裂縫網(wǎng)絡(luò)相互串通、由溶孔、溶洞組合而成的多個(gè)孤立或孔隙連通的水動(dòng)力流動(dòng)單元。每個(gè)單元都具有相對(duì)獨(dú)立的壓力系統(tǒng)或相對(duì)一致的壓力變化規(guī)律、壓力連通關(guān)系以及相似的流體性質(zhì)，在生產(chǎn)中可作為一個(gè)相對(duì)獨(dú)立的流體運(yùn)動(dòng)單元和油氣開采的基本單位?？p洞型油藏開發(fā)是以縫洞單元作為基礎(chǔ)的，典型縫洞型油藏即典型縫洞單元，是指具有一定規(guī)律的縫洞組合、分布模式形成的縫洞單元，如裂縫-溶洞型、地下河-裂縫-地下河型、孤立溶洞-大裂縫型。

典型縫洞油藏的開發(fā)指標(biāo)對(duì)于該類儲(chǔ)集體的開發(fā)具有一定的指導(dǎo)意義[1]。建立一套較為簡(jiǎn)便的開發(fā)指標(biāo)計(jì)算方法對(duì)于輔助油田開發(fā)、提高采收率研究有重大幫助。

目前縫洞型油藏的地質(zhì)特征描述已經(jīng)取得了重大成果[2]。黃朝琴和Arbogast等[3-5]人利用離散縫洞網(wǎng)絡(luò)模型建立了典型縫洞模型滲透率張量公式，為縫洞型油藏?cái)?shù)值模擬提供了基礎(chǔ)。姚軍等[6]人提出了縫洞型油藏試井解釋方法，對(duì)于縫洞型油藏的地層壓力預(yù)測(cè)提供了一種較好的方法。康志江等[7]人對(duì)于縫洞型油藏滲流特征進(jìn)行了研究，認(rèn)為大縫大洞中的流動(dòng)屬于Stokes流動(dòng)。然而對(duì)于碳酸鹽巖縫洞型油藏的開發(fā)指標(biāo)，目前尚缺乏一種有效的計(jì)算方法，而且滲透率張量公式較為復(fù)雜，求解困難。本文對(duì)滲透率張量公式進(jìn)行了簡(jiǎn)化處理，建立了等效滲透率公式，在保證一定精度的基礎(chǔ)上降低了其復(fù)雜程度，減少了計(jì)算量；并利用等效滲透率公式和Buckley-Leverett理論推導(dǎo)出了典型縫洞油藏的關(guān)鍵開發(fā)指標(biāo)計(jì)算方法。

1 等效滲透率

1.1 理論公式建立

針對(duì)縫洞型介質(zhì)等效滲透率，Todd Arbogast等[8]人進(jìn)行了諸多研究，建立了縫洞型油藏中流體控制方程

(1)

研究表明[8]，對(duì)于不同形狀介質(zhì)，滲透率張量表達(dá)形式不同，且較復(fù)雜縫洞介質(zhì)其滲透率張量沒有理論公式，只能通過數(shù)值解法求解。因而對(duì)于發(fā)育一條層狀裂縫的矩形介質(zhì)，考慮基質(zhì)的作用，根據(jù)式(1)可得其等效滲透率為

(2)

式中：Kfx,Kfy為裂縫在x,y方向的滲透率，μm2；η為交界面處滑移系數(shù)；bf為裂縫開度，cm；ly為正方形介質(zhì)寬度，cm；Km為基質(zhì)滲透率，μm2。

由式(2)可以看出該理論考慮了基質(zhì)的作用，同時(shí)需要交界面滑移系數(shù)等參數(shù)，較為復(fù)雜。巖心測(cè)試表明縫洞型油藏裂縫、溶洞發(fā)育，其基質(zhì)區(qū)域的孔隙度和滲透率較小(一般滲透率小于1×10-3μm2，孔隙度小于0.1)。且實(shí)驗(yàn)表明基質(zhì)在束縛水下油相最小啟動(dòng)壓力高達(dá)20.0 MPa，最小啟動(dòng)壓力梯度在10.0 MPa/cm以上，因而流體在基質(zhì)中的流動(dòng)能力遠(yuǎn)遠(yuǎn)小于在裂縫和溶洞中，相對(duì)大溶洞大裂縫來說基質(zhì)的貢獻(xiàn)一般可以忽略不計(jì)。因而等效滲透率計(jì)算時(shí)，只考慮裂縫和溶洞的影響，忽略基質(zhì)的作用。因此對(duì)于發(fā)育多條層狀裂縫的矩形介質(zhì)，只考慮裂縫的影響，根據(jù)公式(2)可得到裂縫等效滲透率為

(3)

式中：nf為裂縫條數(shù)；Af=l·ly為裂縫截面積，cm2。

對(duì)于大型溶洞，可以簡(jiǎn)化為等體積的橢圓管。根據(jù)劉新純[9]的研究，利用泊肅葉公式和牛頓黏滯力定律可以得到傾斜橢圓管狀溶洞的泊肅葉定律為

(4)

式中：Qv為沿橢圓管軸線方向流量，cm3/s；a、b為溶洞截面長(zhǎng)短軸半長(zhǎng)度，cm；pin,pout為進(jìn)出口壓力，MPa；hin,hout為進(jìn)出口高度，cm；L為橢圓管長(zhǎng)度，cm。

結(jié)合流量等效原理，可得到溶洞軸向等效滲透率為

(5)

1.2 數(shù)值模擬驗(yàn)證

利用COMSOL軟件分別對(duì)縫長(zhǎng)10 cm，縫寬1 cm，不同開度的裂縫和長(zhǎng)10 cm，橢圓截面長(zhǎng)軸長(zhǎng)4 cm，不同短軸長(zhǎng)的溶洞中流體單相流進(jìn)行數(shù)值模擬，得到速度剖面分布圖(圖1、 圖2)， 進(jìn)而計(jì)算出等效滲透率，得到等效滲透率曲線(圖3、 圖4)， 對(duì)理論公式(3)、式(5)進(jìn)行驗(yàn)證。

圖1 裂縫開度為0.01 cm時(shí)裂縫開度方向裂縫中流體速度剖面Fig.1 Velocity profile in crack width direction when crack opening width is 0.01 cm

圖2 溶洞橢圓截面短軸長(zhǎng)2 cm時(shí)溶洞中流體速度分布Fig.2 Velocity distribution of fluid in vaggy when short axis of elliptic vuggy is 2 cm

圖3 裂縫理論解和數(shù)值解對(duì)比曲線Fig.3 Comparison of numerical solution with analytical solution of equivalent permeability of crack

圖4 溶洞理論解和數(shù)值解對(duì)比曲線Fig.4 Comparison of numerical solution with analytical solution of equivalent permeability of elliptic vuggy

從圖1裂縫開度方向速度剖面圖中可以看出，裂縫中的速度分布呈現(xiàn)拋物線形，符合層狀裂縫中的流動(dòng)特征。由圖3可以看出，理論解(式(3))和數(shù)值解吻合得較好，但隨著裂縫開度和縫長(zhǎng)比值得增加，二者的差別會(huì)有所增大。這是因?yàn)閿?shù)值解是基于N-S方程進(jìn)行數(shù)值模擬計(jì)算得到的，考慮了壁面的黏滯力影響，其結(jié)果更加符合地層的實(shí)際情況；裂縫開度越大，速度邊界層發(fā)展越充分，黏滯力影響越明顯，因而理論解會(huì)偏離數(shù)值解。但對(duì)于雷諾數(shù)遠(yuǎn)大于1的流體流動(dòng)來說，其黏滯項(xiàng)相對(duì)于速度項(xiàng)可以忽略不計(jì)，且野外露頭和巖心測(cè)試表明，對(duì)于大部分裂縫，其開度和長(zhǎng)度的比值一般小于0.01，因此用式(3)計(jì)算出的滲透率值精確度較高，誤差也在可接受范圍內(nèi)(相對(duì)誤差<5%)。

從圖2速度分布圖中可以看出，在橢圓管中的各剖面速度分布均符合泊肅葉定律，中間速度大，越靠近壁面速度越小，壁面處速度為零。由圖4可知，理論解(式(5))和數(shù)值解之間吻合得很好。隨著溶洞截面長(zhǎng)軸和短軸長(zhǎng)度比值的增加，二者之間差距越來越小。由于泊肅葉定律忽略了慣性力的影響，因而式(5)計(jì)算出的值偏小。但對(duì)于實(shí)際溶洞中的流體流動(dòng)來說，其加速度較小，流體一般處于穩(wěn)定流狀態(tài)，流體慣性力可以忽略不計(jì)，因此式(5)計(jì)算出的滲透率值具有較高的精確度，誤差很小。

綜上，等效滲透率理論公式和數(shù)值模擬結(jié)果吻合得較好，可以利用式(3)和式(5)根據(jù)流量等效原理組合來計(jì)算縫洞系統(tǒng)的等效滲透率K。若裂縫和溶洞傾斜于流動(dòng)方向時(shí)，可以將其分解為垂直于流動(dòng)方向和平行于流動(dòng)方向兩部分，分別計(jì)算；若裂縫、溶洞同時(shí)存在時(shí)，可利用等值滲流阻力法來計(jì)算。

(6)

(7)

式中：K為沿流動(dòng)方向等效滲透率，μm2；Av為溶洞截面積，cm2；A為等效截面積，cm2；L、Lf、Lv分別表示總長(zhǎng)度、裂縫長(zhǎng)度、溶洞長(zhǎng)度，cm。

2 開發(fā)指標(biāo)計(jì)算

2.1 基本假設(shè)

假設(shè)典型縫洞油藏等效滲透率為K；在油水兩相流區(qū)中，油水同時(shí)流動(dòng)，且都服從廣義達(dá)西定律；流體不可壓縮；考慮重力作用。

此外，大量實(shí)驗(yàn)研究證明，縫洞型油藏的相滲曲線屬于水相上凸，油相下凹型[10]。由圖5可以看出，在不同水相錐進(jìn)指數(shù)下，含水率曲線不再是典型的S型曲線，而表現(xiàn)出Γ型曲線的形狀；含水上升率曲線單調(diào)遞減，沒有經(jīng)典曲線中拐點(diǎn)的存在。這些特點(diǎn)可以表明，縫洞型油藏水驅(qū)過程屬于漸變式非活塞水驅(qū)，見水前沒有明顯的臺(tái)階式水驅(qū)前緣，水驅(qū)前緣飽和度為殘余油飽和度Sw=Swc=0.2。也就是說典型縫洞油藏具有Γ型含水率曲線、漸變式非活塞水驅(qū)等特點(diǎn)，其水驅(qū)過程含水飽和度分布曲線如圖6所示。

圖5 不同水相錐進(jìn)指數(shù)時(shí)含水率及含水上升率曲線Fig.5 Water-cut and water-cut rising rate curves under different water coning indexes

圖6 含水飽和度分布圖Fig.6 Distribution curves of water saturation

2.2 理論公式推導(dǎo)

根據(jù)上文可知，縫洞型油藏相滲曲線表現(xiàn)出如圖7的水相上凸，油相下凹特征。

利用MATLAB CFTOOL工具箱對(duì)圖7曲線進(jìn)行擬合得到相滲曲線方程為

圖7 典型縫洞模型相滲曲線Fig.7 Relative permeability curve of typical fracture-vuggy model

Krw=a1eb1Sw-c1e-d1Sw,Kro=a2e-b2Sw-c2e-d2Sw。

(8)

式中：Krw,Kro為水相、油相相對(duì)滲透率；Sw為含水飽和度；a1、b1、c1、d1、a2、b2、c2、d2為擬合系數(shù)，僅與模型縫洞結(jié)構(gòu)組成有關(guān)，可由實(shí)驗(yàn)獲得。

縫洞系統(tǒng)的縱向跨度一般較大[11]，重力作用不能忽略，忽略毛管力考慮重力的含水率[1]為

(9)

考慮重力的油水兩相流微分方程為[12]

(10)

上述方程特征線為dSw=0線，即等飽和度線。則式(10)特征方程為

(11)

根據(jù)B-L理論，解上述特征方程可得到等飽和度面推進(jìn)方程為

(12)

對(duì)式(12)積分可得水驅(qū)前緣(含水飽和度Sw=Swc處)位置為

(13)

2.2.1 見水前開發(fā)指標(biāo)計(jì)算

設(shè)注水速度為Qi，則有：

(1)瞬時(shí)產(chǎn)油量、瞬時(shí)產(chǎn)水量

由注采平衡可得

Qo=Qi，Qw=0。

(14)

(2)累積產(chǎn)油量

Np=Qit。

(15)

(3)見水時(shí)間

由前文可知，井排見水的判斷依據(jù)是Sw2>Swf=Swc，水驅(qū)前緣含水上升率fw′(Swf)=fw′(Swc)，則有

(16)

2.2.2 見水后開發(fā)指標(biāo)計(jì)算

由式(16)可得到見水后時(shí)間和含水飽和度關(guān)系式為

(17)

結(jié)合相滲曲線方程(8)可以得到見水后不同時(shí)刻的井排處含水飽和度Sw2，然后根據(jù)體積平均法可得到油藏平均含水飽和度為

(18)

(1)瞬時(shí)產(chǎn)水量

由含水率定義可得到瞬時(shí)產(chǎn)水量為

Qw=Qi·[1-fwg(Sw2)]。

(19)

(2)瞬時(shí)產(chǎn)油量

Qo=Qi·fwg(Sw2)。

(20)

(3)累計(jì)產(chǎn)油量

由油藏平均含水飽和度可得到

Np=ALφ×{Sw2-Swc-

(21)

(4)累積產(chǎn)水量

(22)

(5)采收率

(23)

2.3 數(shù)值模擬驗(yàn)證

利用數(shù)值模擬軟件COMSOL對(duì)水驅(qū)油過程進(jìn)行模擬，對(duì)開發(fā)指標(biāo)計(jì)算公式(18)、式(19)、式(23)進(jìn)行驗(yàn)證。

應(yīng)用COMSOL軟件中的CFD模塊對(duì)縫洞中的流體流動(dòng)進(jìn)行數(shù)值模擬計(jì)算。流體在較大裂縫和溶洞中的流動(dòng)滿足N-S方程，屬于Navier-Stokes流動(dòng)。COMSOL軟件利用相場(chǎng)法(Phase Field)對(duì)縫洞系統(tǒng)中的兩相流進(jìn)行求解。計(jì)算所用幾何模型由2組相交裂縫和2組橢圓柱狀溶洞交替排列而成，如圖8所示。幾何模型參數(shù)為：縫長(zhǎng)3.54 m，縫開度0.1 m，縫寬0.8 m，傾斜角度45°；溶洞長(zhǎng)2.5 m，橢圓截面長(zhǎng)軸長(zhǎng)2 m，短軸長(zhǎng)1 m。邊界條件為：注入端，定流量Qi=0.1 m3/s；采出端，定壓力，p=0 MPa。根據(jù)式(6)、式(7)可計(jì)算出其等效滲透率為K=5×106μm2。

圖8 典型縫洞幾何模型Fig.8 Schematic diagram of typical fracture-vuggy geometry model

流體基本參數(shù)為：水相黏度μw=1.2 mPa·s，油相黏度μo=5 mPa·s。

2.4 結(jié)果分析

利用COMSOL軟件對(duì)實(shí)際模型進(jìn)行水驅(qū)油模擬，得到產(chǎn)油量、產(chǎn)水量及采收率曲線，與理論公式得到的結(jié)果進(jìn)行對(duì)比，結(jié)果如圖9—圖11所示。

由圖9—圖11中曲線可以看出，理論曲線和兩種數(shù)值模擬曲線吻合得較好。瞬時(shí)產(chǎn)油量隨時(shí)間增加而下降，且下降趨勢(shì)越來越慢，這是由于隨著油藏開采，含油飽和度降低，油相相對(duì)滲透率下降，導(dǎo)致瞬時(shí)產(chǎn)油量減少；同時(shí)隨著含水飽和度的增加，油相相對(duì)滲透率下降速率減小，從而瞬時(shí)產(chǎn)油量下降趨勢(shì)減緩。

瞬時(shí)產(chǎn)水量趨勢(shì)與此相反，隨時(shí)間增加而增加，最終達(dá)到一個(gè)較為穩(wěn)定的狀態(tài)，這是由于縫洞型油藏水相相滲曲線表現(xiàn)為上凸型，隨著開采得進(jìn)行，含水飽和度增加，但水相相對(duì)滲透率增加速率越來越慢，導(dǎo)致產(chǎn)水速度隨時(shí)間增加而減緩。

采收率隨時(shí)間增加而增加，最終達(dá)到一個(gè)較高的平穩(wěn)水平。由于模型復(fù)雜程度不高，且沒有考慮溶洞充填情況，水驅(qū)過程進(jìn)行得較為徹底，因而采收率相對(duì)較高。

圖9 瞬時(shí)產(chǎn)油量曲線Fig.9 Instantaneous oil production curves of typical fracture-vuggy model

圖10 瞬時(shí)產(chǎn)水量曲線Fig.10 Instantaneous water production curves of typical fracture-vuggy model

圖11 采收率曲線Fig.11 Recovery factor curves of typical fracture-vuggy model

2.5 不同典型縫洞組合模式下理論和數(shù)值模擬結(jié)果對(duì)比

設(shè)計(jì)了兩種不同典型縫洞組合模式(圖12、圖13)，分別利用理論公式和COMSOL數(shù)值模擬軟件對(duì)水驅(qū)油過程進(jìn)行計(jì)算模擬，得到瞬時(shí)產(chǎn)油量曲線如圖14、圖15所示。

圖12所示模型為縫-洞結(jié)構(gòu)與孤立溶洞通過一條大裂縫相連接；圖13所示模型為縫-洞結(jié)構(gòu)與豎井-地下河結(jié)構(gòu)通過地下河相連接。

從圖14、圖15中可以看出對(duì)于不同縫洞組合模式，理論曲線和數(shù)值模擬曲線吻合得很好，這表明開發(fā)指標(biāo)理論計(jì)算公式具有較高準(zhǔn)確度。

圖12 孤立溶洞-大裂縫型和正交裂縫-溶洞型組合模式幾何模型Fig.12 Schematic diagram of isolated vug-large crack and orthogonal fractures-vug combination model

圖13 地下河-豎井型和正交裂縫-溶洞型組合模式幾何模型Fig.13 Schematic diagram of underground river-vertical vug and orthogonal fractures-vug combination model

圖14 孤立溶洞-大裂縫型和正交裂縫-溶洞型組合模式下瞬時(shí)產(chǎn)油量曲線Fig.14 Instantaneous oil production curves of isolated vug-large crack and orthogonal fractures-vug combination model

圖15 地下河-豎井型和正交裂縫-溶洞型組合模式下瞬時(shí)產(chǎn)油量曲線Fig.15 Instantaneous oil production curves of underground river-vertical vug and orthogonal fractures-vug combination model

3 結(jié) 論

(1)對(duì)滲透率張量公式進(jìn)行了簡(jiǎn)化處理，在計(jì)算時(shí)忽略了基質(zhì)的影響，只考慮了裂縫和溶洞的滲透率，建立了等效滲透率公式，既保證了一定的精度，又減少了計(jì)算量。

(2)綜合運(yùn)用Buckley-Leverett理論、等效滲透率理論對(duì)典型縫洞油藏模型的開發(fā)指標(biāo)進(jìn)行推導(dǎo)計(jì)算，得到了瞬時(shí)產(chǎn)油量、瞬時(shí)產(chǎn)水量、采收率等開發(fā)指標(biāo)理論公式，并與數(shù)值模擬結(jié)果做了對(duì)比，所得結(jié)果與數(shù)值模擬結(jié)果吻合較好。

(3)隨著開采進(jìn)行，含油飽和度降低，油相相對(duì)滲透率減小，瞬時(shí)產(chǎn)油量隨時(shí)間增加而減少，但下降速度逐漸減?。凰矔r(shí)產(chǎn)水量隨時(shí)間增加而增加，由于水相相滲曲線為上凸型，導(dǎo)致其增加速度逐漸減緩，最終達(dá)到較平穩(wěn)的水平；由于模型水驅(qū)過程進(jìn)行的較為徹底，采收率較高。

(4)對(duì)于實(shí)際油藏中多種縫洞組合模式，提出的開發(fā)指標(biāo)理論計(jì)算公式具有較高準(zhǔn)確度，可以應(yīng)用于不同縫洞組合模式下的開發(fā)指標(biāo)計(jì)算。

[1] 姜漢橋.油藏工程原理與方法[M].北京:石油大學(xué)出版社,2003.

JIANG Hanqiao.Principles and Methods ofReser-voir Engineering[M].Beijing:China University of Petroleum Press,2003.

[2] 龔洪林,袁剛,田彥燦,等.塔中碳酸鹽巖縫洞型儲(chǔ)集體精細(xì)成像[J].地球物理學(xué)進(jìn)展,2012,27 (2):555-561.

GONG Honglin,YUAN Gang,TIAN Yancan,et al.Accurate imaging way of carbonate fracture cave type’s reservoir in Tazhong[J].Progress in Geophysics,2012,27(2):555-561.

[3] 黃朝琴,姚軍,李亞軍,等.基于均勻化理論的縫洞型介質(zhì)滲透性分析[J].中國(guó)科學(xué)(技術(shù)科學(xué)),2010,40(9):1095-1103.

HUANG Zhaoqin,YAO Jun,LI Yajun,et al.Permeability analysis of fractured vuggy porous media based on homogenization theory[J].Science China (Technological Sciences),2010,40(9):1095-1103.

[4] ARBOGAST T,BRUNSON D S.A computational method for approximating a Darcy-Stokes system governing a vuggy porous medium[J].Computational Geosciences,2007,11(3):207-218.

[5] 楊堅(jiān).縫洞型碳酸鹽巖油藏縫洞單元數(shù)值模擬理論與方法研究[D].東營(yíng):中國(guó)石油大學(xué),2007.

YANG Jian.Study of numerical simulation theory and method of fractured-vug unit in fractured vuggy carbonate reservoirs[D].Dongyiing:China Univer-sity of Petroleum,2007.

[6] 姚軍,王子勝.縫洞型碳酸鹽巖油藏試井解釋理論與方法[M].東營(yíng):中國(guó)石油大學(xué)出版社,2007.

YAO Jun,WANG Zisheng.Theory and Method for Well Test Interpretation in Fractured-vuggy Carbonate Reservoirs[M].Dongying:China University of Petr-oleum Press,2007.

[7] 康志江,李江龍,張冬麗,等.塔河縫洞型碳酸鹽巖油藏滲流特征[J].石油與天然氣地質(zhì),2005,26(5):634-640.

KANG Zhijiang,LI Jianglong,ZHANG Dongli,et al.Percolation characteristics of fractured-vuggy carbonate reservoir in Tahe oilfield[J].Oil and Gas Geology,2005,26(5):634-640.

[8] ARBOGAST T,BRUNSON D S,BRYANT S L,et al.A preliminary computational investigation of a macro model for vuggy porous media[J].Develop-ments in Water Science,2015,55(1):267-278.

[9] 劉新純.傾斜橢圓柱管中流體的Poiseuille公式[J].數(shù)理醫(yī)藥學(xué)雜志,2007,20(5):735-736.

LIU Xinchun,The fluid Poiseuille’s formula in Ell-ipse cylindrical of inclination[J].Journal of Mathe-matical Medicine,2007,20(5):735-736.

[10] 李陽.碳酸鹽巖縫洞型油藏開發(fā)理論與方法[M].北京:中國(guó)石化出版社,2014.

LI Yang,Development Theory and Method of Carbo-nate Fractured vuggy Reservoir[M].Beijing:Sinopec press,2014.

[11] 呂愛民.碳酸鹽巖縫洞型油藏油藏工程方法研究[D].東營(yíng):中國(guó)石油大學(xué),2007.

LYU Aimin,Study of reservoir engineering meth-ods for fractured vuggy carbonate reservoirs[D].Dongying:China University of Petroleum,2007.

[12] 張建國(guó),雷光倫,杜殿發(fā),等.油氣層滲流力學(xué)[M].東營(yíng):中國(guó)石油大學(xué)出版社,2009.

ZHANG Jianguo,LEI Guanglun,DU Dianfa,et al.Oil and Gas Seepage Mechanics[M].Dongying:China University of Petroleum Press,2009.