萬 發(fā),蔣中明,李海峰,梁賢浩

(1.長沙理工大學 水利與環(huán)境工程學院,長沙 410114; 2.珠江水利委員會 珠江水利科學研究院,廣州 510000)

0 引 言

壓縮空氣儲能(Compressed Air Energy Storage,CAES)技術(shù)是解決當前清潔能源間歇性和不穩(wěn)定性的關(guān)鍵技術(shù)[1-2],在風能、光伏發(fā)電持續(xù)迅速推廣下,CAES技術(shù)在國內(nèi)外一直備受關(guān)注[3]。目前CAES儲氣庫多采用地下空腔(鹽巖或硬巖洞穴),由于鹽巖選址局限,硬巖洞穴、廢棄礦洞或人工內(nèi)襯洞室成為解決該局限問題的較理想方案[4]。

CAES儲氣庫的安全穩(wěn)定性及密封性是研究CAES技術(shù)的關(guān)鍵問題,而解決關(guān)鍵問題的前提是了解運行過程中儲氣庫內(nèi)壓縮空氣的熱力學過程,在此方面Kushnir等[5],Raju等[6]分別基于圍巖傳熱和絕熱假設(shè)建立了壓縮空氣熱力學解析解數(shù)學模型,He等[7]在研究儲氣庫能量存儲時推導了等溫、絕熱和傳熱洞壁條件下的壓縮空氣熱力學解析模型,Xia等[8]在Kushnir等[5]研究的基礎(chǔ)上提出了洞室內(nèi)部溫度壓力的簡化解析解,劉澧源等[9]基于Kushnir等[5]成果提出了壓縮空氣熱力學過程的差分算法,萬發(fā)等[10]則提出了熱力學過程的數(shù)值計算方法。同時,在儲氣庫襯砌圍巖在熱力耦合作用下的穩(wěn)定性問題也有Rutqvist等[11]建立了儲氣庫在熱流固耦合作用下的儲氣庫穩(wěn)定性數(shù)值模型,周舒威等[12],蔣中明等[13],夏才初等[14]等基于研究了熱力耦合作用下的儲氣庫穩(wěn)定性。此外,在密封性研究上,Kim等[15]指出儲氣庫內(nèi)氣體會在襯砌密封初始缺陷或反復膨脹收縮產(chǎn)生的裂縫中發(fā)生泄漏,周瑜等[16]則對儲氣庫內(nèi)襯洞室的泄漏進行了研究。

上述研究過程極少考慮地下水對儲氣庫壓縮空氣滲流的影響,也未對襯砌密封層開裂后裂隙圍巖中兩相流的滲流傳熱規(guī)律進行研究。當襯砌密封層開裂時,壓縮氣體和地下水會在襯砌裂縫、薄弱連接處以及圍巖裂隙中發(fā)生兩相滲流和傳熱作用。

為研究CAES儲氣庫襯砌開裂對裂隙圍巖內(nèi)兩相滲流傳熱特性影響。本文基于COMSOL數(shù)值仿真軟件中的弱形式PDE(Partial Differential Equation)二次開發(fā)了雙重介質(zhì)兩相流模型,并聯(lián)合雙重介質(zhì)傳熱模型和壓縮空氣熱力學解析模型,建立了CAES地下儲氣庫裂隙圍巖雙重介質(zhì)兩相滲流傳熱模型。計算了包括初始充氣循環(huán)在內(nèi)的210個運行循環(huán)中,襯砌開裂和無開裂2種工況下裂隙圍巖水氣兩相滲流傳熱過程。

1 儲氣庫裂隙巖體兩相滲流傳熱理論

裂隙巖體中的空隙空間包括巖塊中的空隙空間和裂隙空間兩部分?？障逗土严犊臻g的流體可以相互交換,因此,裂隙巖體的真實滲流狀態(tài)可以采用裂隙孔隙雙重介質(zhì)滲流模型進行模擬。對于壓縮空氣地下儲氣庫來說,裂隙圍巖中的流體包括水和氣體兩種介質(zhì),其滲流為兩相流形態(tài),因此,應采用兩相滲流理論進行分析。兩相達西理論模型是較為便利的分析多孔介質(zhì)兩相滲流的理論模型,廣泛被運用于石油工程[17-19],電池化學[20]等領(lǐng)域?？紫督橘|(zhì)中兩相滲流過程可采用COMSOL現(xiàn)有模塊建模,但裂隙兩相滲流過程需進行二次開發(fā)。

1.1 孔隙介質(zhì)兩相流理論

由多孔介質(zhì)內(nèi)流體的質(zhì)量守恒和達西定律可以得出

(1)

其中:

式中:θ為多孔介質(zhì)孔隙率;s1、s2為流體1和流體2的飽和度,s1+s2=1;u為流速矢量矩陣(m/s);ρ為流體密度(kg/m3);ρ1和ρ2分別為不同相的流體密度(kg/m3);μ為流體動力黏度(Pa·s);μ1和μ2為不同相流體的動力黏度(Pa·s);k為多孔介質(zhì)滲透系數(shù)(mD);kr1和kr2分別為相對滲透率;λ為Corey模型孔隙大小分布指數(shù)。

對單元體積內(nèi)濕潤相流體流動,根據(jù)毛管力驅(qū)動可得

(2)

其中:

式中:c1為流體含量(kg/m3);pc為毛細管壓力(Pa),為飽和度s1的函數(shù);pec為毛細管入口壓力(Pa)。

由式(1)、式(2)和狀態(tài)變量方程化簡組合可得壓力和流體含量的方程組為:

?p(s1ρ1+s2ρ2)]=0 ;

(3)

(4)

1.2 裂隙介質(zhì)兩相滲流理論

根據(jù)兩相流體在裂隙介質(zhì)中質(zhì)量守恒和毛細管壓力理論可得

(5)

式中:θf為裂隙通道內(nèi)孔隙率;df為裂隙開度(m)。

裂隙通道內(nèi)兩相流體之間關(guān)系和輔助方程均同上,因此可得關(guān)于裂隙介質(zhì)流體滲流壓力以及飽和度方程

(7)

裂隙中的兩相滲流方程組求解p和c1(s1),采用弱形式PDE(Weak Form)方程組進行二次開發(fā)建模,流程如下:

(1)先定義邊界計算變量,該計算變量用于計算裂隙內(nèi)兩相滲流相關(guān)各種狀態(tài)變量,變量的定義需注意為孔隙介質(zhì)可識別變量,這樣可使得孔隙介質(zhì)與裂隙介質(zhì)之前進行流量交換,變量計算邏輯為對兩相滲流方程本構(gòu)的拆解;

(2)再通過邊界弱形式及試函數(shù)test(p)和test(c1)來構(gòu)建自定義變量之間計算關(guān)系。

1.3 雙重介質(zhì)滲流傳熱理論

滲流過程往往伴隨著熱量傳遞過程,此處采用局部熱平衡假設(shè),假設(shè)基質(zhì)和氣體溫度在極短時間內(nèi)達到熱平衡狀態(tài),滲流傳熱控制方程為

?[θsks+(1-θs)kef]?T=Q。

(8)

式中:θs為多孔介質(zhì)內(nèi)基質(zhì)體積分數(shù),無因次;ρs為基質(zhì)材料密度(kg/m3);Cps為基質(zhì)材料恒壓熱容(J/(kg·K));Cpef為流體有效比熱容(J/(kg·K)),Cpef=s1Cp1+s2Cp2,Cp1和Cp2分為流體相1和2的比熱容;T為溫度(K);u為孔隙內(nèi)氣體流動速度(m/s);?T為溫度梯度(K/m);ks為固體基質(zhì)導熱率(W/(m·K));kef為流體有效導熱率(W/(m·K)),kef=s1k1+s2k2,k1和k2分別為流體相1和2的導熱系數(shù);Q為熱源(J/(m3· s))。

裂隙內(nèi)兩相滲流傳熱理論方程為

?t[df(1-θp)kef+θpks]?tT=0 。

(9)

式中:df為裂隙寬度(m);Qs為流入熱量(J);θp為裂隙內(nèi)固體材料體積分數(shù);▽tT為裂隙內(nèi)溫度梯度(K/m)。

1.4 儲氣庫壓縮空氣熱力學理論

CAES儲氣庫充放氣循環(huán)會導致壓縮空氣熱力學狀態(tài)處于循環(huán)變化過程中,充氣時氣體壓縮,空氣壓力和溫度升高,放氣時壓力和溫度迅速降低。根據(jù)He等[7]研究結(jié)果,CAES儲氣庫壓縮空氣熱力學過程可采用以下理論模型進行計算:

(10)

(11)

(κ-1)hcA(Tw-T)] 。

(12)

式中:T為儲氣庫內(nèi)氣體溫度(K);min為充氣速率(kg/s);CP為氣體恒壓熱容(J/(kg·K));Tin為充氣溫度(K);V為儲氣庫容積(m3);p為儲氣庫內(nèi)氣體壓力(Pa);hc為氣體與洞壁對流換熱系數(shù)(W/(m2·K));Ac為洞壁對流換熱面積(m2);R為氣體常數(shù)(J/(kg·K));κ為絕熱系數(shù);Tw為洞壁溫度(K);m為儲氣庫內(nèi)氣體質(zhì)量(kg);mout、min為放氣速率和充氣速率(kg/s)。

1.5 雙重介質(zhì)兩相流計算程序二次開發(fā)驗證

為驗證雙重介質(zhì)兩相流數(shù)值模型正確性,以Fahad等[21]進行的可視化裂隙性多孔介質(zhì)兩相滲流試驗為驗證對象,采用COMSOL兩相達西和弱形式PDE模塊建立該物理試驗相對應的數(shù)值仿真模型,并將數(shù)值計算結(jié)果與實際物理試驗所得結(jié)果進行對照。試驗設(shè)置如圖1所示,多孔介質(zhì)中存在兩條交叉裂隙,其中A工況為一條裂隙豎直,一條裂隙與之成45°夾角,B工況為兩條裂隙成90°夾角。初始時多孔介質(zhì)和裂隙內(nèi)被水充滿,同時在模型頂部以等流量向內(nèi)注油,兩邊為不透水邊界,下邊界為常壓流出邊界。驗證模型采用參數(shù)如表1所示。

表1 雙重介質(zhì)油水兩相滲流驗證模型參數(shù)

圖1 雙重介質(zhì)油水兩相滲流驗證模型

某時刻裂隙孔隙雙重介質(zhì)油驅(qū)水驗證模型結(jié)果與文獻物理試驗對比如圖2所示,在經(jīng)過上層多孔介質(zhì)后,油在裂隙介質(zhì)中流速較快,形成侵入,因此裂隙及其周圍區(qū)域油相飽和度高,水相飽和度低。同時在兩條裂隙夾角處出現(xiàn)水相飽和度較高區(qū)域,數(shù)值模擬結(jié)果與物理試驗所得現(xiàn)象一致性較好。

圖2 數(shù)值模型兩相滲流結(jié)果與物理試驗對比

流體相對滲透率擬合結(jié)果如圖3所示,數(shù)值模型水相相對滲透率與試驗擬合良好,油相相對滲透率與試驗擬合有所差異,其原因為物理模型本身為較薄的三維結(jié)構(gòu),且裂隙為真實三維結(jié)構(gòu),與本文簡化后的二維模型存在差異。且文獻中油相Corey指數(shù)與水相不同,而本文中則使用了同一Corey指數(shù),因此造成了油相相對滲透率擬合差異。綜合可證明該雙重介質(zhì)兩相流數(shù)值模型的正確性。

圖3 數(shù)值模型流體相對滲透率與物理試驗對比

2 襯砌開裂對裂隙巖體兩相滲流和傳熱特性影響

2.1 數(shù)值模型的建立

為研究地下水影響下CAES儲氣庫襯砌開裂對裂隙巖體兩相滲流和傳熱特性影響,參考某壓縮空氣儲能電站工程及水文地質(zhì)條件,設(shè)計了如圖4所示的計算模型。儲氣庫洞室斷面為圓形,計算模型采用1/2斷面,尺寸為300 m×300 m,儲氣庫埋深150 m,地下水位線位于地表,圍巖為花崗巖,采用低滲混凝土作為襯砌密封層,厚度為50 cm。在儲氣庫圍巖中分別設(shè)置了3條長度均為27 m的裂隙,裂隙寬度為1 mm。由于襯砌密封層外側(cè)與圍巖在施工過程中存在接觸作用,并非整體,此處將其視作裂隙考慮,寬度為0.2 mm。此外,根據(jù)文獻[23]研究結(jié)果與本課題組研究成果,襯砌密封層在運行過程中會出現(xiàn)微小裂縫。為研究襯砌開裂對圍巖滲流和傳熱影響,模型分2種工況:① 襯砌密封層設(shè)置一條開裂裂縫,寬度為0.4 mm,位于襯砌密封層右上方;② 襯砌密封層完整無開裂。

圖4 CAES儲氣庫襯砌圍巖兩相滲流傳熱數(shù)值模型示意圖

計算過程分為兩步:① 儲氣庫開挖后,儲氣庫處滲流壓力為0 MPa,形成初始穩(wěn)定滲流壓力場;② 儲氣庫運行時,壓縮空氣的溫度和壓力隨充放氣循環(huán)變化,驅(qū)替襯砌圍巖中的水的同時進行傳熱過程。模型計算參數(shù)如表2所示。

表2 CAES儲氣庫襯砌圍巖兩相滲流傳熱模型計算參數(shù)

模型左邊界為對稱邊界,右邊界、上邊界和下邊界為靜水壓力邊界和水相飽和度邊界;對傳熱模型,右邊界、上邊界和下邊界為溫度邊界;儲氣庫內(nèi)壁為空氣溫度和壓力變化的非穩(wěn)態(tài)邊界,采用ODE(Ordinary Differential Equation)建模,假定運行過程中儲氣庫內(nèi)壁氣體飽和度始終為1。CAES電站運行周期一般為1 d,一個周期分為4個運行階段,分別為充氣段(8 h)-高壓儲氣段(4 h)-放氣段(4 h)-低壓儲氣段(8 h)。初始充氣段充氣時間24 h,充氣速率為52 kg/s,穩(wěn)定運行充氣速率47.5 kg/s,放氣速率為95 kg/s。本文模型計算過程為初始階段和隨后的正常運行210 d,其非穩(wěn)態(tài)邊界計算結(jié)果如圖5所示。

壓縮空氣壓力始終處于在7～10 MPa的運行區(qū)間內(nèi)循環(huán)變化,初始充氣時空氣被劇烈壓縮,溫度迅速升高,放氣時溫度迅速降低,運行過程中再次充氣時由于氣體壓力變化幅度低于首次充氣,氣體壓縮比例降低,溫度升高幅度較低,同時由于洞壁導熱作用,空氣溫度逐步降低,最后形成穩(wěn)定波動狀態(tài),穩(wěn)定運行時溫度變化范圍為272～303 K。

2.2 襯砌開裂對裂隙圍巖滲流壓力和飽和度影響

210 d時典型運行階段2種工況下裂隙圍巖的滲流壓力分布如圖6所示。

圖6 210 d典型運行階段裂隙圍巖內(nèi)滲流壓力分布

襯砌開裂的模型中圍巖裂隙和孔隙內(nèi)滲流壓力差異顯著,距洞壁相同距離的最大壓力差異為1.9 MPa;而不考慮開裂的模型圍巖裂隙與孔隙內(nèi)滲流壓力分布較為均勻,相同距離最大差異為0.4 MPa。開裂模型中,充氣末圍巖裂隙滲流壓力為9.6 MPa,滲流壓力>8 MPa的區(qū)域從裂隙盡頭向圍巖延伸10.2 m,充氣和放氣段末圍巖裂隙滲流壓力相差3 MPa;在無開裂模型中,充氣末圍巖裂隙滲流壓力為8.03 MPa,>8 MPa滲流壓力范圍僅延伸1.2 m,充氣和放氣段圍巖裂隙最大差異只有0.18 MPa。值得注意的是放氣階段儲氣庫氣體壓力迅速降低,但2種模型圍巖內(nèi)仍然保持較高滲流壓力,存在滲流遲滯效應。對比兩個模型的計算結(jié)果,表明襯砌開裂會顯著影響圍巖內(nèi)滲流壓力分布,影響主要體現(xiàn)在壓力分布均勻性、高滲壓傳遞范圍,以及圍巖裂隙和周圍孔隙內(nèi)滲流壓力的高低。

分別對圍巖裂隙方向和多孔介質(zhì)方向取截線,兩種模型在不同運行階段各截線滲流壓力分布如圖7所示。

圖7 圍巖裂隙和多孔介質(zhì)滲流壓力分布

在時間、距離均相同的條件下,開裂模型的裂隙和孔隙內(nèi)滲流壓力總體高于無開裂模型。截線1的充氣和放氣段滲流壓力分布線在開裂模型中于52 m之后重合,而在無開裂模型中則于37.5 m后重合,截線2的充放氣滲流壓力在開裂模型中距離約30 m后重合,在無開裂模型重合距離縮小為16 m。由此可知開裂后模型圍巖內(nèi)滲流壓力隨充放氣過程循環(huán)變化的范圍為52 m,無開裂時該變化的范圍為37.5 m,開裂模型中多孔介質(zhì)滲流壓力受充放氣影響范圍為30 m,無開裂模型中該影響范圍為16 m。因此,襯砌開裂后不僅會造成圍巖內(nèi)裂隙和孔隙內(nèi)滲流壓力升高,而且會使?jié)B流壓力受充放氣過程而波動變化的范圍大幅擴大,在本模型計算中,增長幅度可達46.7%。襯砌開裂會導致周圍孔隙內(nèi)滲流壓力隨運行過程起伏波動更大。

運行210 d后2種模型襯砌和裂隙圍巖內(nèi)水相飽和度分布變化如圖8所示,壓縮空氣主要經(jīng)裂隙通道驅(qū)替地下水,當裂隙通道完全被空氣占據(jù)后,壓縮空氣會在裂隙盡頭處向附近圍巖內(nèi)入侵。由于圍巖裂隙盡頭高程不同,滲流壓力不同,進而導致滲流速度不同,最終使得3條裂隙盡頭處氣相飽和度為1的范圍有所差異,最上方的裂隙盡頭處范圍最大,最下方的范圍最小。襯砌開裂對氣體飽和度影響主要集中在圍巖裂隙及其周圍。無開裂模型水平方向飽和度整體低于開裂模型,運行210 d后氣體擴散范圍為54 m,開裂模型運行210 d后水平方向氣體擴散影響范圍為66.7 m,較無開裂模型擴大了12.7 m。

圖8 開裂對裂隙圍巖水相飽和度的影響

2.3 襯砌開裂對儲氣庫整體泄漏速率的影響

對儲氣庫洞室內(nèi)邊界通過流體進行積分,可得儲氣庫每時刻泄漏速率,計算公式為

me=∮uρDds。

(13)

式中:me為漏氣速率(kg/s);D為洞室軸向長度(m);s為洞室橫截面周長(m)。

以第200—第210天計算結(jié)果為分析對象,兩種模型儲氣庫泄漏速率變化如圖9所示,儲氣庫充氣時內(nèi)部壓力高于雙重介質(zhì)內(nèi)滲流壓力,此時泄漏方向向外,放氣時泄漏方向向內(nèi)。襯砌開裂后的儲氣庫整體泄漏速率低于無開裂時,其原因為襯砌開裂后,由于襯砌裂隙與襯砌圍巖的間隙相連,直接導致襯砌外側(cè)滲流壓力高于無開裂模型,進而導致洞壁內(nèi)外壓差降低,洞壁處孔隙介質(zhì)滲流速度相比無開裂模型降低約10倍,盡管開裂點泄漏速率較高,但是其開裂寬度較小導致其流量與儲氣庫整體泄漏率相比仍然低了一個數(shù)量級。因此在該模型中,襯砌開裂會使儲氣庫整體泄漏速率降低。但值得注意的是,該結(jié)果規(guī)律僅在圍巖裂隙不與大氣相連時適用,當圍巖裂隙與大氣相連時,襯砌外側(cè)滲流壓力與大氣壓力形成穩(wěn)定梯度時,開裂泄漏不對襯砌外側(cè)滲流壓力造成影響,襯砌開裂將不會對洞壁孔隙介質(zhì)滲流速度造成影響,此時開裂會使泄漏速率增加,增加部分主要為襯砌裂隙處泄漏量。

圖9 開裂對儲氣庫泄漏速率的影響

2.4 襯砌開裂對襯砌及裂隙圍巖傳熱特性影響

2種模型在第210天充氣段末的襯砌及裂隙圍巖溫度分布如圖10所示,無論是否考慮開裂,襯砌和裂隙圍巖溫度分布均呈圓形圈層分布,水平截線和圍巖孔隙介質(zhì)截線溫度分布線相重合,因此可知,開裂和裂隙的分布對圍巖溫度的傳遞和分布影響非常微弱,其原因是滲流速度較低,同時空氣的比熱容較小,低流速的空氣所含有熱能變化較低,因此裂隙的存在不會影響圍巖溫度的分布。圍巖溫度與初始值相比最大差異為1 K,平衡運行時溫度變化影響范圍為20.8 m,影響范圍內(nèi)的圍巖溫度低于初始溫度,遠處溫度保持初始溫度。相比圍巖溫度分布,襯砌溫度梯度較大,充氣段末襯砌內(nèi)外溫差為16 K。

圖10 第210天充氣段末襯砌和裂隙圍巖溫度分布規(guī)律

溫度的分布變化主要集中襯砌密封層內(nèi),襯砌不同位置的測點溫度變化如圖11所示,由圖11可知,充放氣循環(huán)時溫度起伏波動變化規(guī)律與充放氣運行過程相一致,變化幅度則各有差異。襯砌內(nèi)側(cè)溫度變化幅度為31 K,襯砌中間25 cm處溫度變化幅度為15 K,襯砌外側(cè)溫度變化幅度為1.5 K,起伏變化稍滯后于內(nèi)側(cè)溫度的變化。整個運行過程中,溫度的變化都集中于襯砌密封層內(nèi),襯砌內(nèi)側(cè)會受到初始循環(huán)的高溫作用,此時襯砌內(nèi)外溫差較高,需注意由此產(chǎn)生的溫度裂縫。

圖11 襯砌不同位置溫度變化規(guī)律

3 結(jié) 論

以多孔介質(zhì)兩相滲流理論為基礎(chǔ),采用弱形式二次開發(fā)了雙重介質(zhì)兩相滲流模型,并進行了驗證。同時利用數(shù)學模塊全局ODE建立壓縮空氣熱力學模型,并耦合多孔介質(zhì)傳熱模型,建立了CAES硬巖儲氣庫的裂隙圍巖兩相滲流傳熱模型。以某壓縮空氣儲能電站儲氣庫為研究背景,研究了襯砌開裂對電站運行過程中裂隙圍巖兩相滲流和傳熱特性的影響。結(jié)果表明:

(1)基于弱形式PDE二次開發(fā)的雙重介質(zhì)兩相流模型可較好地模擬雙重介質(zhì)內(nèi)兩相滲流互驅(qū)過程。

(2)襯砌開裂會顯著影響圍巖內(nèi)滲流壓力分布,主要體現(xiàn)在圍巖裂隙和孔隙內(nèi)滲流壓力差增大,最大增加1.5 MPa;圍巖裂隙盡頭孔隙介質(zhì)高滲壓范圍擴大;以及使?jié)B流壓力受充放氣過程而波動變化的范圍放大。

(3)壓縮空氣從裂隙中向外擴散,并在裂隙盡頭處向圍巖內(nèi)驅(qū)替地下水,襯砌開裂會導致氣體擴散范圍顯著增加,本模型中210 d氣體水平擴散范圍由未開裂前的54 m擴大至66.7 m。

(4)襯砌開裂和圍巖裂隙的存在對襯砌和圍巖的溫度分布基本沒有影響,2種條件下相同距離裂隙和孔隙介質(zhì)溫度分布一致。圍巖的溫度變化平穩(wěn),穩(wěn)定運行后其溫度起伏變化在1 K內(nèi)。溫度變化主要集中襯砌密封層內(nèi),初始循環(huán)是襯砌內(nèi)外側(cè)溫差極大,本模型中>100 K,需注意此時不均勻溫度造成的溫度裂縫。