萬(wàn) 發(fā),蔣中明,2,唐 棟,3
(1長(zhǎng)沙理工大學(xué)水利工程學(xué)院;2長(zhǎng)沙理工大學(xué)水沙科學(xué)與水災(zāi)害防治湖南省重點(diǎn)實(shí)驗(yàn)室;3長(zhǎng)沙理工大學(xué)洞庭湖水環(huán)境治理與生態(tài)修復(fù)湖南省重點(diǎn)實(shí)驗(yàn)室,湖南 長(zhǎng)沙410114)
間歇性和不穩(wěn)定性是太陽(yáng)能、風(fēng)能、潮汐能等清潔可再生能源并入電網(wǎng)使用的主要障礙,壓氣儲(chǔ)能(compressed air energy storage,以下簡(jiǎn)稱CAES)技術(shù)是解決該問(wèn)題的大規(guī)模儲(chǔ)能技術(shù)[1-5]。CAES 電站儲(chǔ)氣庫(kù)洞室作為能源存儲(chǔ)裝置,其安全性和穩(wěn)定性是電站穩(wěn)定運(yùn)行的基礎(chǔ),揭示儲(chǔ)氣庫(kù)洞室內(nèi)部熱力學(xué)特性[6],對(duì)于評(píng)價(jià)儲(chǔ)氣庫(kù)洞室的安全性和穩(wěn)定性意義重大。
在研究CAES 儲(chǔ)氣庫(kù)熱力學(xué)特性方面國(guó)內(nèi)外學(xué)者已經(jīng)進(jìn)行了很多研究,并取得較多成果,如Raju[7]、Kushnir[8]、Zhou[9]、Xia 等[10]學(xué)者用數(shù)學(xué)解析解的方式建立了溫度、壓力解析方程,研究了洞室內(nèi)氣體熱力學(xué)變化過(guò)程,但解析解的計(jì)算結(jié)果為洞室內(nèi)部溫度壓力的均值,并不能描述溫度、壓力、能量在洞室空間內(nèi)部分布規(guī)律。德國(guó)Huntorf 電站給出儲(chǔ)氣庫(kù)洞室溫度壓力實(shí)測(cè)值時(shí),指出洞室內(nèi)溫度和出口處溫度在運(yùn)行過(guò)程中存在差異[11],根據(jù)Geissbühler 等[12]在瑞士南部的實(shí)際試驗(yàn)和蔣中明等[13]在湖南平江儲(chǔ)氣室試驗(yàn)結(jié)果以及劉佳的研究[14]都表明地下儲(chǔ)氣庫(kù)洞室內(nèi)部溫度、換熱器TES 內(nèi)溫度在空間分布上有著較大差異性。儲(chǔ)氣庫(kù)洞室內(nèi)部不均勻溫度場(chǎng)會(huì)導(dǎo)致圍巖受到不均勻溫度應(yīng)力,極端情況下甚至?xí)诙词揖植繀^(qū)域形成極高溫區(qū)域,引起洞室圍巖和襯砌產(chǎn)生熱損傷破壞。因此對(duì)儲(chǔ)氣庫(kù)洞室熱力學(xué)特性分布規(guī)律的研究可有效指導(dǎo)儲(chǔ)氣庫(kù)洞室設(shè)計(jì),避免出現(xiàn)極端不均勻溫度場(chǎng)而引起的熱力學(xué)破壞。
He 等[15]通過(guò)存儲(chǔ)?值的計(jì)算了不同傳熱邊界下的洞室?guī)烊?。確定庫(kù)容后需要選擇適合的洞室形式,夏才初等[16]研究了隧道式和大罐式洞室的內(nèi)壓下圍巖力學(xué)特性。確定庫(kù)容和洞室形式后,仍有不同的長(zhǎng)寬設(shè)計(jì)方案,不同洞室設(shè)計(jì)方案對(duì)于洞室熱力學(xué)特性的影響卻未見(jiàn)報(bào)道。
為揭示CAES 儲(chǔ)氣庫(kù)洞室運(yùn)行過(guò)程熱力學(xué)特性分布規(guī)律,研究不同洞室設(shè)計(jì)方案對(duì)熱力學(xué)特性的影響,本文擬采用COMSOL 非等溫共軛傳熱模型,模擬Huntorf 電站單儲(chǔ)氣庫(kù)日循環(huán)運(yùn)行過(guò)程,并采用實(shí)測(cè)數(shù)據(jù)對(duì)模型方法進(jìn)行驗(yàn)證,同時(shí)進(jìn)一步對(duì)圓柱形洞室不同體型參數(shù)和入氣口布置位置工況進(jìn)行建模,模擬其充氣過(guò)程的洞室熱力學(xué)變化過(guò)程,對(duì)比分析不同設(shè)計(jì)參數(shù)對(duì)洞室熱力學(xué)時(shí)空分布特性的影響,為洞室的設(shè)計(jì)提供依據(jù)。
本文研究氣體在洞室內(nèi)部流動(dòng)和傳熱狀態(tài)時(shí),基于以下假設(shè):①由于充氣口進(jìn)入的氣體流速較高,因此假定儲(chǔ)氣庫(kù)內(nèi)氣體流動(dòng)狀態(tài)服從標(biāo)準(zhǔn)湍流k-ε 雙方程模型;②氣體具備可壓縮性,壓縮氣體的溫度壓力滿足氣體狀態(tài)方程;③熱對(duì)流過(guò)程除強(qiáng)制對(duì)流外考慮自然對(duì)流,儲(chǔ)氣過(guò)程中,溫差和重力導(dǎo)致儲(chǔ)氣庫(kù)內(nèi)氣體流動(dòng);④計(jì)算流動(dòng)加熱時(shí)忽略氣體的黏性耗散;⑤儲(chǔ)氣庫(kù)洞壁圍巖熱傳導(dǎo)過(guò)程滿足傅里葉定律。
流體連續(xù)性方程(1)和湍流k-ε 二方程(2)~(3)
式中,ρ 為氣體密度,kg/m3;▽為漢密頓算子;u 為速度矢量,m/s。
式中,pK為湍流速度變化所引起的壓力,Pa;μT=ρCμk2/ε 湍流黏度系數(shù);ε=ep 為湍流脈動(dòng)動(dòng)能耗散率;常數(shù)Cμ=0.09;σε=1.3;Cε1=1.44;Cε2=1.92。
洞室內(nèi)氣體流動(dòng)過(guò)程中的能量主要由流入過(guò)程中的動(dòng)量、氣體壓縮做功產(chǎn)生的能量變化、熱量向洞壁傳遞過(guò)程三個(gè)主要部分組成,其中氣體流動(dòng)的動(dòng)量控制方程為
氣體壓縮做功的能量控制方程
式中,Qp為氣體壓力做功,J;T 為溫度,K。
熱傳導(dǎo)方程
式中,Qw為氣體固體傳熱流量,J/S;q 為界面處熱流密度,J/s;uτ為壁面處綜合摩擦系數(shù);Tw為固體壁面溫度,K;T 為流體溫度,K。
為研究洞室內(nèi)能量存儲(chǔ)狀況,采用文獻(xiàn)[15]中?值計(jì)算方法來(lái)計(jì)算能量存儲(chǔ)?
式中,m 為充氣速率,kg/s;cp為氣體恒壓比熱容;T 為洞室內(nèi)氣體溫度,K;T0為氣體初始溫度,K;p 為洞室氣體壓力,Pa;p0為氣體初始?jí)毫?,Pa。
為驗(yàn)證理論模型的正確性,采用德國(guó)Huntorf電站給出的實(shí)測(cè)數(shù)據(jù)做對(duì)比。Huntorf 電站采用地下鹽巖溶腔作為儲(chǔ)氣庫(kù)洞室,由于其洞室邊界實(shí)際的不規(guī)則性,此處將洞室簡(jiǎn)化為圓柱形洞室。采用COMSOL 軟件湍流和流固傳熱模塊并耦合非等溫流進(jìn)行建模計(jì)算,洞室模型網(wǎng)格示意圖如圖1 所示,共有域單元77794 個(gè),模型半徑為67.24 m,其中圍巖厚50 m,模型分為兩個(gè)域,內(nèi)部洞室為流體域,外部圍巖為固體域。流體域頂部布置入氣口,針對(duì)湍流模型:充氣入口處為流量邊界,流體域外邊界為0 通量邊界;針對(duì)流固傳熱模型:充氣入口為溫度邊界(充氣段)、絕熱邊界(儲(chǔ)氣段)和流出邊界(放氣段),圍巖外邊界為絕熱邊界,采用“層殼”導(dǎo)熱的方式處理傳熱場(chǎng)的流體固體交界面。采用分離求解器PARDISO 求解算法進(jìn)行求解,收斂條件為容差控制,相對(duì)容差為0.005。洞室和圍巖物理性質(zhì)參數(shù)如表1 所示。Huntorf 電站的日循環(huán)運(yùn)行過(guò)程包含了三種運(yùn)行狀態(tài)的交替往復(fù),分別為放氣狀態(tài)、儲(chǔ)氣狀態(tài)、充氣狀態(tài),儲(chǔ)氣室各個(gè)運(yùn)行狀態(tài)的時(shí)間和充放氣體流率如圖2 所示。為簡(jiǎn)化計(jì)算,第一個(gè)小時(shí)的放氣速率簡(jiǎn)化為150 kg/s 均勻放氣,其中3 個(gè)充氣階段的充氣溫度分別為50.96 ℃、45.95 ℃、49.08 ℃。
圖1 洞室有限元模型網(wǎng)格示意圖Fig.1 Sketch of FEM meshes of Huntorf cavern
表1 Huntorf壓氣儲(chǔ)能電站儲(chǔ)氣庫(kù)參數(shù)[8,10]Table 1 Parameters of Huntorf CAES plant[8,10]
圖2 洞室日循環(huán)運(yùn)行充放氣速率Fig.2 Mass flow rate of air to the cavern in daily circulation
Huntorf 日循環(huán)過(guò)程洞室內(nèi)氣體溫度隨時(shí)間變化規(guī)律以及與文獻(xiàn)[16]給出的實(shí)測(cè)值對(duì)比結(jié)果如圖3 所示,由于三維洞室內(nèi)溫度場(chǎng)不均勻性的存在[16],其中溫度分別考慮了洞室內(nèi)部氣體溫度均值,以及洞室與圍巖界面處溫度均值;洞室內(nèi)氣體壓力隨時(shí)間變化規(guī)律以及與文獻(xiàn)[16]給出的實(shí)測(cè)值對(duì)比結(jié)果如圖4 所示,氣體壓力也考慮了洞室內(nèi)部壓力均值和界面處壓力均值。由圖可知,洞室溫度模擬結(jié)果中,界面溫度均值與實(shí)測(cè)值擬合結(jié)果比整個(gè)洞室內(nèi)溫度均值的擬合結(jié)果好,但在第二段充氣結(jié)束時(shí)刻界面溫度均值與實(shí)測(cè)值有著3 K 的差異,隨后的儲(chǔ)氣段和充氣段前期都有著2~3 K 差異。氣體壓力在洞室內(nèi)部均值和界面處均值與實(shí)測(cè)值之間的擬合結(jié)果較好,該模型可以較好地模擬洞室內(nèi)氣體壓力的變化。
圖3 洞室氣體溫度模擬結(jié)果與實(shí)測(cè)對(duì)比Fig.3 The comparison between simulated air temperature and measured value of the Huntorf cavern
圖4 洞室氣體壓力模擬結(jié)果與實(shí)測(cè)對(duì)比Fig.4 The comparison between simulated air pressure and measured value of the Huntorf cavern
由以上模擬結(jié)果與實(shí)測(cè)對(duì)比可知,該模型可以較好地模擬CAES 電站充放氣過(guò)程中洞室氣體壓力變化過(guò)程,對(duì)于溫度而言,由于三維模型溫度空間不均勻性,模擬計(jì)算結(jié)果與實(shí)測(cè)值之間的誤差是在可接受范圍之內(nèi)的??偠灾?,該模型可以較好地模擬計(jì)算CAES 洞室內(nèi)熱力學(xué)特性變化規(guī)律。
圖5 最佳擬合點(diǎn)溫度變化規(guī)律Fig.5 Validation of the best fitting point of air temperature
CAES 儲(chǔ)氣庫(kù)洞室涉及熱力學(xué)過(guò)程復(fù)雜,其中充放氣速率、溫度、運(yùn)行壓力,庫(kù)容估算等都有據(jù)可循,但是確定庫(kù)容之后的設(shè)計(jì)方案對(duì)洞室內(nèi)熱力學(xué)特性的影響結(jié)果卻仍不清楚。為研究不同設(shè)計(jì)方案對(duì)洞室內(nèi)熱力學(xué)特性影響,本文研究豎直圓柱形洞室不同體型參數(shù)K(洞室長(zhǎng)度與洞室半徑之間比值)和不同入氣口位置對(duì)洞室內(nèi)熱力學(xué)特性的影響。
為研究圓柱形洞室體型參數(shù)K和不同入氣口位置對(duì)洞室內(nèi)熱力學(xué)特性的影響,擬設(shè)置4 種不同K值和兩種入氣口布置形式共五種的工況,五種洞室設(shè)計(jì)方案變量參數(shù)如表2所示。五種方案下洞室體積為2500 m3,洞室圍巖研究半徑30 m,模型分為兩個(gè)域:內(nèi)部流體域,外部圍巖固體域。圍巖固體域外邊界采用開(kāi)放式邊界,等于大氣溫度20 ℃,初始溫度為20 ℃,初始?jí)毫? atm(1 atm=101.325 kPa),充氣速率為10 kg/s,充氣溫度為10 ℃,其他模型參數(shù)如表1 所示。K=10 時(shí)模型如圖6所示,模擬計(jì)算過(guò)程僅考慮充氣過(guò)程,充氣時(shí)長(zhǎng)為8 h。
4種K值設(shè)計(jì)方案洞室充氣8 h的內(nèi)部氣體溫度均值T和壓力p變化規(guī)律如圖7、圖8所示。由圖可知,四種洞室內(nèi)溫度均值隨充氣時(shí)間增加而升高,其中上升速率最快的洞室為K=20 的狹長(zhǎng)洞室,上升速率最小的洞室為K=5的洞室。氣體溫度均值在四種K值的洞室內(nèi)有著較大的差異,在充氣末時(shí)刻K=20 的狹長(zhǎng)洞室氣體溫度均值比K=5 的洞室溫度均值高35.8 K,溫度均值隨著K值的降低呈降低趨勢(shì),但不是K 值越小,溫度均值越小,K=2.5 的洞室體溫度均值比K=5 的洞室高1.5 K。洞室氣體壓力p 在充氣過(guò)程中區(qū)別不大,充氣8 h 時(shí),K=20 的洞室內(nèi)部氣體壓力比其余3個(gè)洞室高0.15 MPa。
表2 四種設(shè)計(jì)方案下體型參數(shù)及初始條件Table 2 Shape parameters and initial condition under four design schemes caverns
圖6 體型參數(shù)K=10時(shí)模型網(wǎng)格示意圖Fig.6 Grid diagram of the model when the shape parameter K=10
不同K值洞室氣體溫度最大值的變化規(guī)律如圖9 所示,K=20 和K=10 的洞室在充氣過(guò)程中均會(huì)形成較高的最大溫度,其中K=20 充氣末時(shí)刻氣體溫度最大值為649.62 K;K=2.5 和K=5 的洞室氣體溫度最大值上升較緩,末時(shí)刻分別為378.19 K 和365.66 K,合理的選擇K 值可以有效避免在充氣過(guò)程中洞室內(nèi)出現(xiàn)較高氣體溫度。充氣過(guò)程中洞室存儲(chǔ)的?值B隨時(shí)間變化如圖10所示,由圖可知,不同K 值對(duì)于?值的存儲(chǔ)影響不大,其中K=2.5、5、10 的洞室內(nèi)存儲(chǔ)的?值相同,只有K=20 時(shí),?值較其余3個(gè)洞室高2.5×10-9J。
圖7 不同K值洞室氣體溫度均值T變化規(guī)律Fig.7 Variation of average temperature T in the caverns with different K value
圖8 不同K值洞室氣體壓力p變化規(guī)律Fig.8 Variation of average pressure p in the caverns with different K value
圖9 不同K值洞室氣體溫度最大值Tmax變化規(guī)律Fig.9 Variation of maximum temperature Tmax in the caverns with different K value
圖10 不同K值洞室氣體存儲(chǔ)?變化規(guī)律Fig.10 Variation of exergy storage B in the caverns with different K value
不同K值的洞室在充氣過(guò)程中,熱力學(xué)特性影響較大的為溫度,K值較大的洞室極易形成較高的氣體溫度,對(duì)于壓力和存儲(chǔ)的?值影響較小幾乎可以忽略。較大的K值的洞室有較高的氣體溫度均值和最大值,但是不是K值越小,氣體溫度均值和最大值越小,當(dāng)K值小到一定時(shí),氣體溫度均值和最大值重新開(kāi)始上升。
充氣末時(shí)刻480 min 時(shí)刻4 種K 值洞室中截面的溫度分布等值線如圖11所示,對(duì)于K值較小的短粗型洞室,氣體溫度容易在入氣口兩側(cè)的頂部聚積,形成高溫區(qū)域;而K值較大的狹長(zhǎng)型洞室,由于洞室較長(zhǎng),氣體在底部進(jìn)行壓縮升溫,所以洞室底部容易進(jìn)行高溫區(qū)域,且K值越大,洞室越窄越長(zhǎng),洞室底部的氣體溫度越高。
對(duì)K=10 設(shè)計(jì)方案的洞室分別采用頂部入氣口布置和洞壁中部入氣口布置,同樣的邊界條件和充氣條件下模擬計(jì)算充氣8 h 洞室內(nèi)的熱力學(xué)過(guò)程。兩種入氣方案的模型洞室氣體溫度均值和壓力隨時(shí)間變化規(guī)律分別如圖12、圖13 所示。頂部入氣和中間入氣其內(nèi)部壓力差異極小,可以忽略不計(jì),最終時(shí)刻內(nèi)部壓力分別為3.79和3.81 MPa;而儲(chǔ)氣庫(kù)內(nèi)部體平均溫度則在45 min后開(kāi)始產(chǎn)生差距,且差距越來(lái)越大,頂部入氣體平均溫度最終為368 K,而中部入氣體平均溫度最終為352 K,表明中間入氣可有效降低儲(chǔ)氣室內(nèi)部溫度均值,該參數(shù)下可降低溫度均值16 K。
兩種不同入氣方式下的儲(chǔ)氣庫(kù)內(nèi)能量?B和焓H隨時(shí)間的變化如圖14、圖15所示,焓H是溫度的函數(shù),由內(nèi)部溫度均值可知,頂部入氣方式下的儲(chǔ)氣庫(kù)內(nèi)氣體焓值較高,8 h末時(shí)為2.36×1010J,而中間入氣的末時(shí)刻焓值為1.87×1010J。儲(chǔ)氣庫(kù)內(nèi)?值隨時(shí)間呈曲線上升,兩種充氣方式下的焓值雖然不同,但是?值卻相差不大,頂部入氣下儲(chǔ)氣庫(kù)內(nèi)最終?值為6.84×1010J,而中間入氣下儲(chǔ)氣庫(kù)內(nèi)最終?值為6.73×1010J。
圖11 不同K值洞室中界面溫度分布Fig.11 Contour of middle section temperature distribution in cavern with different K values
圖12 不同入氣口位置洞室氣體溫度均值變化規(guī)律Fig.12 Variation of average temperature in the cavern with different air inlet locations
圖13 不同入氣口位置洞室壓力變化規(guī)律Fig.13 Variation of pressure in the cavern with different air inlet locations
圖14 不同入氣口位置洞室存儲(chǔ)?值隨時(shí)間變化規(guī)律Fig.14 Variation of exergy storage in the cavern with different air inlet locations
圖15 不同入氣口位置洞室氣體焓值隨時(shí)間變化規(guī)律Fig.15 Variation of enthalpy in the cavern with different air inlet locations
充氣過(guò)程中的洞室內(nèi)氣體溫度最大值變化規(guī)律如圖16 所示,由圖可知,頂部入氣的布置方式會(huì)使洞室內(nèi)氣體溫度迅速上升,到220 min 時(shí)開(kāi)始平緩,最終時(shí)刻溫度最大值為528.87 K,中間入氣的布置方式洞室內(nèi)氣體溫度上升較緩,最終時(shí)刻溫度最大值為369.26 K。儲(chǔ)氣庫(kù)不同的入氣口位置對(duì)充氣過(guò)程中的壓力影響不大,但是對(duì)于儲(chǔ)氣庫(kù)內(nèi)部的溫度影響較大,內(nèi)部溫度最大值可降低159.61 K,中部入氣可以避免最內(nèi)部出現(xiàn)極高溫度,這對(duì)于鹽巖這類(lèi)軟巖儲(chǔ)氣庫(kù)圍巖穩(wěn)定性有極大幫助。
對(duì)中間入氣的洞室充氣末時(shí)刻480 min 時(shí)中截面的溫度分布等值線圖進(jìn)行繪制,結(jié)果如圖16 所示,由圖10 和圖16 的對(duì)比可知,中部入氣的洞室溫度分布更為均勻,氣體在從入氣口注入,向上下流動(dòng)壓縮,溫度在入氣口一側(cè)上下方累計(jì),頂部入氣會(huì)在底部形成高溫區(qū),中間入氣可避免洞室內(nèi)出現(xiàn)較高溫度區(qū),可有效避免洞壁出現(xiàn)較大的溫度應(yīng)力。
圖16 不同入氣口布置洞室溫度最大值變化規(guī)律Fig.16 Variationofmaximumtemperatureinthecavernwith differentairinletlocations
圖17 不同入氣口布置洞室中截面溫度分布等值線Fig.17 Contour of middle section temperature distribution in cavern with different air inlet
兩種入氣口布置方式對(duì)于洞室內(nèi)充氣過(guò)程的熱力學(xué)特性影響主要體現(xiàn)在溫度上,對(duì)于壓力和?值存儲(chǔ)幾乎沒(méi)有影響。頂部入氣會(huì)產(chǎn)生較高的溫度和較大的不均勻溫度場(chǎng),高溫區(qū)域主要在洞室底部;中間入氣方式溫度均值較低,極大值也較頂部入氣低199.5 K,洞室內(nèi)溫度場(chǎng)分布較均勻,對(duì)于鹽巖這類(lèi)性質(zhì)較軟、蠕變流變性較好的,較低的均勻溫度場(chǎng)對(duì)于其穩(wěn)定性和氣密性是有利的,因此,在條件允許下,應(yīng)當(dāng)盡可能采用中間入氣的充氣方式。
(1)該數(shù)值方法可以較好地模擬CAES 電站運(yùn)行過(guò)程中的洞室內(nèi)熱力學(xué)過(guò)程。在Huntorf 電站日循環(huán)運(yùn)行過(guò)程中,洞室氣體溫度場(chǎng)存在不均勻性,變溫層主要存在于洞壁附近,洞室內(nèi)壓力分布均勻。
(2)對(duì)于圓柱形洞室而言,設(shè)計(jì)方案中體型參數(shù)K 會(huì)顯著影響洞室熱力學(xué)過(guò)程。相同充氣條件下,K值對(duì)洞室壓力、存儲(chǔ)?影響極小,但是會(huì)顯著影響氣體溫度。K值越大,洞室內(nèi)氣體溫度最大值和均值越大,但不是K 值越小,溫度均值就越小,對(duì)一定的充氣速率和充氣溫度下,存在最佳K值使得洞室氣體溫度場(chǎng)分布均勻,溫度均值最小。K 值較大時(shí),洞室內(nèi)氣體高溫區(qū)域集中在洞室底部,而K值較小時(shí)的洞室高溫區(qū)域集中在入氣口一段的兩側(cè)。
(3)圓柱形洞室的入氣口布置位置會(huì)顯著影響洞室內(nèi)氣體溫度分布,但不會(huì)影響氣體壓力和存儲(chǔ)的?值。中間入氣的布置方式會(huì)明顯降低洞室內(nèi)溫度均值和最大值,分別可降低16 K和159.61 K。在洞室容積有限而充放氣速率較快的情況下,應(yīng)當(dāng)盡可能采用K值為5、中部入氣的設(shè)計(jì)方案。