壓縮空氣儲(chǔ)能填充床儲(chǔ)熱特性研究

姜 萍，趙振家

(1.河北大學(xué)電子信息工程學(xué)院，河北保定 071000;2.河北大學(xué)羅克韋爾自動(dòng)化實(shí)驗(yàn)室，河北保定 071000)

0 引言

在我國當(dāng)前階段，能源需求不斷增長難以滿足工業(yè)生產(chǎn)，大力發(fā)展可再生能源是轉(zhuǎn)變能源消費(fèi)觀念，發(fā)展低碳經(jīng)濟(jì)的必由之路。太陽能和風(fēng)能以分布廣泛，取材豐富因而作為可再生能源開發(fā)的重點(diǎn)。但是風(fēng)能和太陽能分布不均勻，受到時(shí)間和地區(qū)的限制，并且風(fēng)能和太陽能具有間歇性和隨機(jī)性，因而不能夠?yàn)樯a(chǎn)提供穩(wěn)定的能量輸出。為此人們提出如果能夠把風(fēng)能和太陽能存儲(chǔ)起來，便能夠改善風(fēng)光天然的間歇性和隨機(jī)性，從而能夠提高電能品質(zhì)，保證用電的可靠性，為生產(chǎn)生活提供穩(wěn)定的能量輸出。儲(chǔ)能技術(shù)便成為了當(dāng)今研究的焦點(diǎn)，發(fā)展儲(chǔ)能技術(shù)能夠有效地解決能源需求和能源供應(yīng)之間的矛盾。

壓縮空氣儲(chǔ)能系統(tǒng)是一種能夠?qū)崿F(xiàn)大容量和長時(shí)間電能存儲(chǔ)的電力儲(chǔ)能系統(tǒng)，在峰谷電能回收調(diào)節(jié)、平衡負(fù)荷、頻率調(diào)制、分布式儲(chǔ)能等方面具有很大應(yīng)用潛力。利用電網(wǎng)低谷或可再生能源供給過剩電能帶動(dòng)壓縮機(jī)，通過壓縮空氣儲(chǔ)存多余的電能。在用電高峰期或可再生能源供給不足時(shí)，高壓氣體經(jīng)過填充床溫度升高，高溫高壓氣體通過膨脹機(jī)做功發(fā)電。壓縮空氣儲(chǔ)能系統(tǒng)的原理圖如圖1所示。儲(chǔ)熱技術(shù)是有效提高壓縮空氣儲(chǔ)能效率的重要方式，從而既擁有遠(yuǎn)高于壓縮空氣儲(chǔ)能的能量密度，又大大提高空氣儲(chǔ)能的效率。蓄熱的主要方式分為3 種:顯熱蓄熱、潛熱蓄熱、和化學(xué)反應(yīng)蓄熱，填充床是最常見的顯熱蓄熱。顯熱蓄熱的運(yùn)行方式簡(jiǎn)單，成本低，熱效率高，因此適合大規(guī)模應(yīng)用。

文獻(xiàn)［1］分析了以巖石，鋼球和相變硝酸鹽作為儲(chǔ)熱材料時(shí)的出熱量、熱分層、相間溫差等特性以及儲(chǔ)熱—釋熱周期對(duì)熱效率的影響，表明填充床是解決可再生能源儲(chǔ)熱問題的儲(chǔ)熱方式。文獻(xiàn)［2～3］開展了不同壓力下填充床蓄熱周期實(shí)驗(yàn)，測(cè)量了填充床內(nèi)溫度分布和進(jìn)、出口空氣溫度隨時(shí)間的變化，實(shí)驗(yàn)表明提高罐壁保溫性能，以較大空氣質(zhì)量流量和較低空氣壓力條件下可以得到較高的蓄熱效率。文獻(xiàn)［4～5］為了提高壓縮空氣儲(chǔ)能系統(tǒng)效率建立了基于渦旋壓縮機(jī)的壓縮空氣儲(chǔ)能系統(tǒng)的動(dòng)態(tài)模型，為適應(yīng)控制策略的研究，建立了兼顧過/欠膨脹機(jī)損耗的渦旋膨脹機(jī)平均模型。文獻(xiàn)［6］提出了一種模型自激式發(fā)電機(jī)，仿真實(shí)驗(yàn)通過連續(xù)不斷的壓縮空氣放電來驅(qū)動(dòng)感應(yīng)發(fā)電機(jī)，表明此設(shè)備可以為可再生能源直接連接電網(wǎng)供電。文獻(xiàn)［7～8］將壓縮空氣儲(chǔ)能系統(tǒng)與風(fēng)力發(fā)電廠相結(jié)合進(jìn)行了討論研究，表明壓縮空氣儲(chǔ)能系統(tǒng)能夠有效改善風(fēng)力發(fā)電廠輸出功率的穩(wěn)定性，起到了削峰填谷的作用。文獻(xiàn)［9～10］在MATLAB/SIMULINK 中討論了壓縮空氣儲(chǔ)能系統(tǒng)的動(dòng)態(tài)建模問題，給出了壓縮機(jī)，儲(chǔ)氣室和膨脹機(jī)的數(shù)學(xué)模型，設(shè)計(jì)了電網(wǎng)側(cè)的控制系統(tǒng)，建立了電力電子裝置調(diào)節(jié)氣動(dòng)馬達(dá)轉(zhuǎn)子的轉(zhuǎn)速進(jìn)行充放電。文獻(xiàn)［11］在流體動(dòng)力學(xué)基礎(chǔ)上建立的壓縮空氣儲(chǔ)能系統(tǒng)仿真模型進(jìn)行了熱力分析。但是對(duì)于帶有蓄熱裝置的壓縮空氣儲(chǔ)能系統(tǒng)的建模和熱力學(xué)分析工作鮮有報(bào)道，與傳統(tǒng)壓縮空氣儲(chǔ)能系統(tǒng)相比，帶有蓄熱裝置的壓縮空氣儲(chǔ)能系統(tǒng)不依賴于化石燃料且無污染，可以大規(guī)模推廣。

本文采用MATLAB 數(shù)值計(jì)算方法求解了傳熱模型的偏微分方程，得到了模型的數(shù)值解和模擬圖形，并基于ANSYS FLUENT 14.0 對(duì)石子填充床內(nèi)的流體動(dòng)力學(xué)和傳熱特性進(jìn)行了仿真，對(duì)空氣通過填充床的儲(chǔ)熱過程進(jìn)行了數(shù)值模擬，通過計(jì)算得到了填充床出口溫度與入口空氣質(zhì)量流量的傳遞函數(shù)，實(shí)驗(yàn)結(jié)果表明計(jì)算所得的傳遞函數(shù)能夠完整表述填充床儲(chǔ)熱動(dòng)態(tài)特性。

圖1 壓縮空氣儲(chǔ)能原理圖Fig.1 Schematic diagram of compressed air energy storage

1 基于一維非穩(wěn)態(tài)導(dǎo)熱－擴(kuò)散傳熱模型的驗(yàn)證

Schumann 最早通過數(shù)學(xué)推導(dǎo)得到了填充床內(nèi)流體溫度與實(shí)踐距離的關(guān)系，Littman 和Gunn 發(fā)現(xiàn)在低雷諾數(shù)下流體和固體之間的擴(kuò)散現(xiàn)象明顯，因而采用了一維非穩(wěn)態(tài)導(dǎo)熱擴(kuò)散兩相模型來描述填充床內(nèi)的傳熱現(xiàn)象。

經(jīng)過壓縮機(jī)和膨脹機(jī)后的空氣為高溫高壓氣體，換熱介質(zhì)為多孔介質(zhì)。考慮空氣與顆粒之間的溫差，在假設(shè)流體速度分布均勻，填充床為絕熱系統(tǒng)，孔隙率均勻，氣體與固體的物性參數(shù)不隨溫度變化，固相與流體相具有連續(xù)性，忽略熱輻射和自然對(duì)流，并且固體顆粒和流體無徑向溫度梯度的情況下，Kunii 和Smith 通過數(shù)學(xué)推導(dǎo)得到理想填充床內(nèi)流體及固體溫度與時(shí)間和距離的函數(shù)關(guān)系式，建立了流動(dòng)方向上固體與流體的一維、非穩(wěn)態(tài)和常物性能量微分方程，式(1)為流體相方程，式(2)為固體相方程。

式中:Ts和Tf分別為填充床內(nèi)石子和空氣溫度;x為沿流體流動(dòng)方向距離;t 為時(shí)間;cs和cf分別為石子和空氣比熱;ε 為孔隙率;mf為單位截面空氣質(zhì)量流量;a 為單位體積填充床內(nèi)石子總面積;h為空氣與顆粒之間的傳熱系數(shù)。和分別為石子和空氣的導(dǎo)熱系數(shù)。取ε =0.4，ρs=2 650 kg/m3，cs=0.92 kJ/(kg℃)，h =2w/(m·k)，a =0.001 m3，L =1.5 m，mf=0.05kg/s，ρf=17 kg/=4 w/(m·k)，帶入偏微分方程整理得

初始條件和邊界條件為

Tf(0，t)= 430(K)Ts(x，0)= 298(K)

對(duì)于偏微分方程的某些定解問題不能通過解析解嚴(yán)格求出，只能用近似方法來求出滿足實(shí)際需要的近似解，即它的數(shù)值解。偏微分方程的數(shù)值解法主要有有限差分法、有限元法以及邊界元法等。本文采用MATLAB 中的pdepe 函數(shù)來求解偏微分方程。

MATLAB 中的pdepe 函數(shù)可以求解一般的PDEs，具有較大的通用性，pdepe 函數(shù)的調(diào)用格式為sol=pde(m，@ pdefun，@ pdeic，@ pdebc，x，t)，@pdefun 是問題的描述函數(shù)，@pdeic 是PDE 的初值條件，@ pdebc 是PDE 的邊界條件描述函數(shù)。輸出參數(shù)sol:是一個(gè)三維數(shù)組，sol(:，:，i)表示ui的解，通過sol，就可以使用pdeval()計(jì)算某個(gè)點(diǎn)的函數(shù)值。帶入數(shù)值后得到固體和流體溫度與時(shí)間和流體流動(dòng)方向距離關(guān)系如圖2、3所示。

2 基于ANSYS 填充床儲(chǔ)熱特性仿真

計(jì)算流體力學(xué)是近代流體動(dòng)力學(xué)的一個(gè)重要分支，F(xiàn)LUENT 是國際知名的CFD 軟件，集成了非常豐富的參數(shù)化方案和物理模型，可用于多種流體問題的數(shù)值模擬。由于FLUENT 穩(wěn)定性能好，運(yùn)算精度高，適用范圍廣，高效省時(shí)，有先進(jìn)的動(dòng)/變形網(wǎng)格技術(shù)和強(qiáng)大的網(wǎng)格支持能力而備受人們親睞。

圖2 固體溫度與時(shí)間和流體流動(dòng)方向距離關(guān)系圖Fig.2 The diagram of the solid temperature for time and the fluid flow direction distance

圖3 流體溫度與時(shí)間和流體流動(dòng)方向距離關(guān)系圖Fig.3 The diagram of the fluid temperature for time and the fluid flow direction distance

填充床是內(nèi)部還有眾多空隙的多孔介質(zhì)，蓄熱的主要方式分為3 種:顯熱蓄熱、潛熱蓄熱、和化學(xué)反應(yīng)蓄熱。本文選擇具有顯熱儲(chǔ)熱特性的球形石子顆粒為內(nèi)部填充滿材料，填充床進(jìn)口管道直徑為159 mm，高為1 500 mm，球形石子顆粒粒徑為10 mm，石子填充區(qū)域直徑為500 mm，高度為1 200 mm。本文利用ANSYS FLUENT14.0 對(duì)填充床的儲(chǔ)熱特性進(jìn)行了仿真研究，為了降低數(shù)值模擬的復(fù)雜度并縮短計(jì)算時(shí)間，假設(shè)填充床與環(huán)境絕熱，不考慮援助罐體熱容的影響，填充床內(nèi)孔隙率分布均勻，忽略熱輻射的影響。

本文利用FLUENT 前處理軟件Gambit 進(jìn)行網(wǎng)格的劃分，讀入并檢查網(wǎng)絡(luò)之后，設(shè)置求解器參數(shù)，湍流模型中選擇k－ epsilon(2eqn)，定義材料屬性，設(shè)置材料為巖石，導(dǎo)熱系數(shù)為2 W/(m·K)，密度為2 650 kg/m3;設(shè)置區(qū)域條件，孔隙率為0.4;設(shè)置邊界條件，對(duì)空氣進(jìn)口邊界條件進(jìn)行設(shè)置，本實(shí)驗(yàn)系統(tǒng)所采用填充床入口直徑為6 cm，所用壓縮機(jī)的壓縮比約為1.44，空氣入口質(zhì)量流量為0.05 kg/s，入口空氣溫度為430 K，填充床為絕熱邊界;最后進(jìn)行求解計(jì)算，設(shè)置求解松弛因子和收斂臨界值，設(shè)置流場(chǎng)初始化，用Patch 按鈕對(duì)填充床的初始溫度進(jìn)行設(shè)置，填充床內(nèi)石子和壁面初始溫度為298 K，監(jiān)視出口溫度的變化曲線，得到出口溫度曲線如圖4所示。

圖4 空氣出口溫度與時(shí)間的關(guān)系Fig.4 The relationship between outlet temperature and time

由出口溫度變化曲線可以看出，出口熱風(fēng)達(dá)到最高溫度所需時(shí)間很短，由于采用的是氣固換熱熱平衡方程，所以空氣最高溫度與入口空氣溫度相等，均為430 K，300 s 之前的時(shí)間，出口空氣溫度均保持為298 K，隨著填充床的不斷升溫，出口空氣溫度逐漸升高，最終達(dá)到600 多秒時(shí)，換熱結(jié)束，達(dá)到熱平衡。表(1)顯示了空氣入口質(zhì)量流量為0.05 kg/s，入口溫度為430 K 時(shí)，出口空氣溫度與時(shí)間的關(guān)系。

圖5～8 分別顯示了當(dāng)時(shí)間為300 s，350 s，400 s，500 s 時(shí)填充床的溫度分布情況。

由溫度場(chǎng)云圖可以看出，填充床是自下而上逐層升溫的，流動(dòng)方向上溫度梯度大，橫向溫度梯度幾乎為0。

表1 填充床出口溫度與時(shí)間關(guān)系表Tab.1 The relation table for outlet temperature and time of the packed bed

圖5 300 s 時(shí)填充床的溫度分布圖Fig.5 The temperature distribution of the packed bed on 300 second

3 填充床儲(chǔ)熱特性數(shù)值分析

建立被控過程數(shù)學(xué)模型的方法包括機(jī)理法和測(cè)試法，壓縮空氣儲(chǔ)能系統(tǒng)填充床的熱力特性為復(fù)雜的偏微分方程組描述形式，不便于采用機(jī)理法建模。為了獲得填充床的動(dòng)態(tài)特性，通入熱空氣之后，得到填充床的空氣出口溫度呈現(xiàn)S 形單調(diào)曲線，本文采用階躍響應(yīng)曲線兩點(diǎn)法來求得填充床動(dòng)態(tài)特性的特征參數(shù)，得到填充床出口空氣溫度與入口空氣質(zhì)量流量之間的傳遞函數(shù)。

圖6 350 s 時(shí)填充床的溫度分布圖Fig.6 The temperature distribution of the packed bed on 350 second

先將階躍響應(yīng)曲線化為無量綱的階躍響應(yīng)y*(t)，階躍響應(yīng)的穩(wěn)態(tài)值y(∞)與階躍輸入的幅值x0之比為被控過程的靜態(tài)放大系數(shù)，K =y(∞)/x0，K = 13.2，取y*(t1)=0.4，y*(t2)=0.8 兩點(diǎn)，從曲線上確定對(duì)應(yīng)的t1=367(s)，t2=395(s)，由t1/t2= 0.785，確定系統(tǒng)為16 階，T =。

圖7 400 s 時(shí)填充床的溫度分布圖Fig.7 The temperature distribution of the packed bed on 400 second

圖8 500 s 時(shí)填充床的溫度分布圖Fig.8 The temperature distribution of the packed bed on 500 second

最終得到填充床的傳遞函數(shù):

填充床在壓縮過程中吸收壓縮空氣的熱量，而在膨脹過程中填充床要對(duì)高壓空氣進(jìn)行加熱，本文得到的填充床蓄熱過程的數(shù)學(xué)模型并不適用于膨脹發(fā)電填充床蓄冷過程的數(shù)學(xué)模型，膨脹過程空氣出口溫度的變化曲線如圖9所示。

圖9 膨脹發(fā)電時(shí)空氣出口溫度與時(shí)間的關(guān)系Fig.9 The relationship between outlet temperature and time in expansion system

經(jīng)過計(jì)算得到膨脹發(fā)電過程填充床出口空氣溫度與入口空氣質(zhì)量流量的傳遞函數(shù)為

4 實(shí)驗(yàn)結(jié)果分析

本實(shí)驗(yàn)系統(tǒng)所采用填充床入口直徑為6 cm，所用壓縮機(jī)的壓縮比約為1.44，空氣入口質(zhì)量流量為0.05 kg/s，入口空氣溫度為430 K，填充床內(nèi)石子和壁面初始溫度為298 K，巖石孔隙率為0.4，通過ANSYS FLUENT 仿真實(shí)驗(yàn)得到了填充床出口溫度的變化曲線，并在MATLAB 中繪制了填充床邊界部分固體溫度的實(shí)驗(yàn)曲線與理論推導(dǎo)曲線的對(duì)比圖，分析了實(shí)驗(yàn)所得的數(shù)據(jù)曲線與傳遞函數(shù)響應(yīng)曲線存在差異的原因。

圖10所示為填充床邊界部分固體溫度的實(shí)驗(yàn)曲線與理論推導(dǎo)曲線的對(duì)比圖，由于計(jì)算模型較簡(jiǎn)單，而仿真所用模型相對(duì)精確，導(dǎo)致由理論推導(dǎo)所得曲線與實(shí)驗(yàn)的數(shù)據(jù)曲線存在一定誤差，但所反映填充床邊界部分固體溫度變化大體一致。圖11所示為實(shí)驗(yàn)所得的數(shù)據(jù)曲線與傳遞函數(shù)響應(yīng)曲線的對(duì)比圖，由于傳遞函數(shù)為有限階函數(shù)，不能完全描述實(shí)驗(yàn)所得數(shù)據(jù)曲線，導(dǎo)致實(shí)測(cè)值偏小于傳遞函數(shù)值，從圖中可以看出實(shí)驗(yàn)所得傳遞函數(shù)能夠近似描述熱空氣流過填充床后空氣溫度的變化。

圖10 填充床邊界部分固體溫度的實(shí)驗(yàn)曲線與理論推導(dǎo)曲線的對(duì)比圖Fig.10 The comparison between the experimental curve and theoretical derivation curve for the solid temperature of boundary part

壓縮空氣儲(chǔ)能系統(tǒng)熱力特性多為復(fù)雜的偏微分方程組描述形式，難以用于控制系統(tǒng)的設(shè)計(jì)，本文研究了填充床的蓄冷蓄熱過程中熱力特性參數(shù)的動(dòng)態(tài)變化，利用仿真軟件模擬及數(shù)值擬合，采用測(cè)試法進(jìn)行數(shù)學(xué)建模，得到了適用于控制系統(tǒng)設(shè)計(jì)的數(shù)學(xué)模型。

圖11 實(shí)驗(yàn)所得的數(shù)據(jù)曲線與傳遞函數(shù)響應(yīng)曲線的對(duì)比圖Fig.11 The comparison between the experimental curve and the transfer function response curve

5 結(jié)論

與電池、超級(jí)電容、飛輪等儲(chǔ)能方式相比，壓縮空氣儲(chǔ)能是唯一可在綜合效益方面和抽水蓄能相媲美的儲(chǔ)能方式［12］。但傳統(tǒng)壓縮空氣氣體膨脹需要與化石燃料混合燃燒，對(duì)化石燃料和地理環(huán)境的依靠和排放污染等方面限制了大規(guī)模推廣。帶蓄熱裝置的壓縮空氣儲(chǔ)能系統(tǒng)可以解決傳統(tǒng)壓縮空氣儲(chǔ)能系統(tǒng)需要補(bǔ)燃這一缺點(diǎn)，填充床蓄熱的運(yùn)行方式簡(jiǎn)單，成本低，熱效率高，適合大規(guī)模應(yīng)用。

研究非補(bǔ)燃式壓縮空氣的建模對(duì)于控制系統(tǒng)的設(shè)計(jì)研究具有重要意義。通過對(duì)被控過程及相關(guān)設(shè)備的數(shù)學(xué)模型進(jìn)行仿真和模擬試驗(yàn)，在計(jì)算機(jī)上進(jìn)行計(jì)算、分析，以獲取代表接近真實(shí)過程的定量關(guān)系，可以為控制系統(tǒng)的設(shè)計(jì)和調(diào)試提供數(shù)據(jù)，從而大大降低設(shè)計(jì)實(shí)驗(yàn)成本并節(jié)約時(shí)間、加快設(shè)計(jì)進(jìn)程。本文以石子填充床作為研究對(duì)象首先利用MATLAB 數(shù)值計(jì)算方法，求解得到了填充床一維非穩(wěn)態(tài)導(dǎo)熱－擴(kuò)散兩相模型相應(yīng)的偏微分方程，得到了該模型的數(shù)值解和模擬圖形，并基于ANSYS FLUENT 14.0 對(duì)石子填充床內(nèi)的流體動(dòng)力學(xué)和傳熱特性進(jìn)行了仿真，對(duì)熱空氣通過填充床儲(chǔ)熱器的儲(chǔ)熱過程進(jìn)行了數(shù)值模擬，通過計(jì)算得到了填充床出口溫度與入口空氣質(zhì)量流量的傳遞函數(shù)。

實(shí)驗(yàn)結(jié)果表明，計(jì)算所得到的傳遞函數(shù)曲線可以完整描述實(shí)驗(yàn)所得曲線，并且與理論公式計(jì)算所得曲線基本吻合，具有較強(qiáng)的實(shí)用性，可以運(yùn)用于非補(bǔ)燃式壓縮空氣儲(chǔ)能系統(tǒng)的建模中。

［1］趙巖，王亮，陳海生，等.填充床顯熱及相變儲(chǔ)熱特性分析［J］.工程熱物理學(xué)報(bào)，2012，33 (12):2052－2057.

［2］楊亮，劉佳，岳雷，等.巖石填充床蓄熱性能試驗(yàn)研究［J］.可再生能源，2013，31 (8):88－92.

［3］劉佳，楊亮，盛勇，等.壓縮空氣儲(chǔ)能的多級(jí)填充床蓄熱實(shí)驗(yàn)研究［J］.科學(xué)技術(shù)與工程，2013，13(6):1462－1468.

［4］褚曉廣，張承慧，李珂，等.基于渦旋壓縮機(jī)的新型壓縮空氣儲(chǔ)能系統(tǒng)動(dòng)態(tài)建模與效率分析［J］.電工技術(shù)學(xué)報(bào)，2011，26 (7):126－132.

［5］褚曉廣，張承慧，李珂，等.渦旋膨脹機(jī)發(fā)電系統(tǒng)效率優(yōu)化控制策略［J］.電工技術(shù)學(xué)報(bào)，2012，27(6):25－31.

［6］Vongmanee V，Monyakul V.A Modeling of Self Excited Induction Generators Driven by Compressed Air Energy Based on Field Oriented Control Principle［C］.IEEE 2nd International Conference on Power and Energy，Johor Baharu，Malaysia，1－3 December，2008:834－838.

［7］Saadat M，Li P Y.Modeling and Control of a Novel Compressed Air Energy Storage System for Offshore Wind Turbine［C］.2012 American Control Conference，Montreal，QC，27－ 29 June，2012:3032－3037

［8］Dahraie M V，Najafi H R，Azizkandi R N.Study on Compressed Air Energy Storage Coupled With a Wind Farm［C］.2012 Second Iranian Conference on Renewable Energy and Distributed Generation，Tehran，Iran，6－8 March，2012:147－152.

［9］Martínez M，Molina M G，Mercado P E.Dynamic Performance of Compressed Air Energy Storage (CAES)Plant for Applications in Power Systems［C］.2010 IEEE/PES Transmission and Distribution Conference and Exposition，Latin America，Sao Paulo，8－ 10 Nov.，2010:496－503.

［10］Hasan N S，Hassan M Y，Majid M S，et al.Mathematical Model of Compressed Air Energy Storage in Smoothing 2MW Wind Turbine［C］.2012 IEEE International Power Engineering and Optimization Conference，Melaka，Malaysia，6－ 7 June，2012:339－343.

［11］Grazzini G，Milazzo A.A Thermodynamic Analysis of Multistage Adiabatic CAES［C］.Proceedings of the IEEE，2012，100 (2):461－472.

［12］盧強(qiáng).淺談壓縮空氣儲(chǔ)能—專訪中國科學(xué)家盧強(qiáng)院士［N］.中國能源報(bào)，2012－06－02.