姜 萍，趙振家

(1.河北大學電子信息工程學院，河北保定 071000;2.河北大學羅克韋爾自動化實驗室，河北保定 071000)

0 引言

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

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

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

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

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

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

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

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

式中:Ts和Tf分別為填充床內(nèi)石子和空氣溫度;x為沿流體流動方向距離;t 為時間;cs和cf分別為石子和空氣比熱;ε 為孔隙率;mf為單位截面空氣質(zhì)量流量;a 為單位體積填充床內(nèi)石子總面積;h為空氣與顆粒之間的傳熱系數(shù)。和分別為石子和空氣的導熱系數(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)

對于偏微分方程的某些定解問題不能通過解析解嚴格求出，只能用近似方法來求出滿足實際需要的近似解，即它的數(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:是一個三維數(shù)組，sol(:，:，i)表示ui的解，通過sol，就可以使用pdeval()計算某個點的函數(shù)值。帶入數(shù)值后得到固體和流體溫度與時間和流體流動方向距離關系如圖2、3所示。

2 基于ANSYS 填充床儲熱特性仿真

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

圖2 固體溫度與時間和流體流動方向距離關系圖Fig.2 The diagram of the solid temperature for time and the fluid flow direction distance

圖3 流體溫度與時間和流體流動方向距離關系圖Fig.3 The diagram of the fluid temperature for time and the fluid flow direction distance

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

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

圖4 空氣出口溫度與時間的關系Fig.4 The relationship between outlet temperature and time

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 實驗結果分析

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

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

圖10 填充床邊界部分固體溫度的實驗曲線與理論推導曲線的對比圖Fig.10 The comparison between the experimental curve and theoretical derivation curve for the solid temperature of boundary part

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

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

5 結論

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

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

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

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

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

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

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

［5］褚曉廣，張承慧，李珂，等.渦旋膨脹機發(fā)電系統(tǒng)效率優(yōu)化控制策略［J］.電工技術學報，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］盧強.淺談壓縮空氣儲能—專訪中國科學家盧強院士［N］.中國能源報，2012－06－02.