水下壓縮空氣儲能系統儲氣裝置的CFD數值模擬

王金舜，王 虎，熊 偉，王志文

(大連海事大學 船舶機電研究所，遼寧 大連 116026)

引言

隨著當前能源結構的變革，可再生能源正發揮著前所未有的重要作用。但可再生能源固有的“間歇性”、“隨機性”和“低能量密度”的特點極大地阻礙了其高效可靠的利用[1-2]。這些由可再生能源產生的波動的電能和波動的電力負荷之間存在著嚴重不平衡。壓縮儲能技術是解決這種不平衡問題最有效的方法之一[3]，水下壓縮空氣儲能技術是一種適合在沿海地區、海島、海洋平臺等區域建設的，規模化的可再生能源存儲技術。與陸上壓縮空氣儲能系統不同之處在于，水下壓縮空氣儲能系統中的儲氣裝置需要放置在水下環境中，可以保證即使發生失效事件也不會對整個系統或者周圍設施與居民造成嚴重的危害[4-5]，這雖然在實現等壓存儲的基礎上提高了系統的安全性，但復雜多變的海洋環境也對水下儲氣裝置的設計提出了不小的挑戰。因此，研究儲氣裝置周圍的流場分別在普通與極端惡劣水下環境時的變化特征是十分有必要的。

CFD數值模擬方法被廣泛應用于海洋工程結構物的設計過程當中，如水下立管、海上平臺樁腿等。在水下壓縮空氣儲能技術領域，數值模擬同樣是分析水下儲氣裝置周圍流場的有效研究方法之一。

VASEL-BE-HAGH等[6-9]使用了標準k-ω和LES Dyna-SM兩種湍流模型分析了在雷諾數為2.3×105的流動環境中氣球狀儲氣裝置周圍的流動結構及裝置的受力特性，闡明了裝置后旋渦可能的發生形式及脫落過程，并得到了該裝置的升阻力系數及渦脫頻率。通過對比兩種湍流模型的分析結果表明，使用標準k-ω模型對力系數的計算結果較大渦模擬結果偏大。WANG等[10-12]使用標準k-ω模型研究了在9種不同雷諾諾數下的氣球狀儲氣裝置的受力特性，并將數值模擬和實驗的力特性與流動結構結果進行了比較和討論。此外，他還使用k-ωSST湍流模型對3個不同自由端形狀的等效特征長度的鈍體進行了數值模擬，并對比分析以研究自由端的形狀效應，結果表明，k-ωSST湍流模型可以正確地預測尾流中氣球狀鈍體和流動結構的時均力特性，但不能有效捕捉到瞬態流場特性。繼續使用大渦模擬對氣球狀儲氣裝置流體動力學特性進行了數值模擬，結果表明大渦模擬既能夠準確預測流場時均特性，也能夠有效捕捉到流場瞬態特性。此外，由于自由端和錐形效應，還在鈍體后產生了獨特的渦旋結構。

本研究采用大渦模擬方法，模擬在不同海流速度條件下的剛性水下儲氣裝置的繞流動態過程。首先通過分析流場域內速度的分布特性，以了解儲氣裝置周圍不同位置處的流場信息以及該信息關于初始流速的敏感程度。然后以相同時間步長為間隔，分析一段時間內不同時刻的旋渦的發生及脫落的動態過程。最后與裝置的受力變化曲線相互驗證，分析旋渦脫落對裝置受力的影響并得到裝置所受流體力的主頻。

1 數學模型

1.1 控制方程

本研究采用大渦模擬(LES)的湍流模型對水下儲氣裝置的繞流流場進行數值模擬。大渦模擬通過濾波函數將旋渦在空間過濾，對尺度較大的旋渦進行直接模擬，再利用亞格尺度模型模擬小尺寸旋渦對大尺度旋渦的影響[13]。不可壓縮流體的N-S方程為：

(1)

式中，ρ—— 流體密度

p—— 壓力

ν—— 運動黏度

ui,uj—— 速度分量

xi,xj—— 位移分量

t—— 時間

經濾波函數處理N-S方程，可以得到大渦模擬的控制方程為：

(2)

(3)

經過濾所帶來的附加未知項τij為SGS應力，用來描述小尺度渦對大渦的影響。通常使用由SMAGORINSKY[14]提出的SM模型來計算SGS應力，如式(4)所示：

(4)

δij—— Kroneker符號

CS—— 常數，一般取0.1

Δ—— 網格濾波尺度

1.2 裝置的幾何模型與計算網格

本研究中所使用的儲氣裝置模型如圖1a、圖1b所示，模型經簡化處理為直徑D=10 m、總高度H=15 m的球頂圓柱體結構，可在內部儲存大約1000 m3的壓縮空氣以滿足設備儲存容量要求，如圖1c所示,裝置外徑10 m，同時也是計算不同工況下的雷諾數Re時所使用的特征長度。

圖1 儲氣裝置幾何模型與計算模型

在中國南海海域內，根據當地水文條件及洋流活動資料顯示，該海域內海底水流速度全年在0.5～2.0 kn 范圍內，冬季極端條件下海流速度可達2.8 kn以上。此外，水下溫度全年在8 ℃以下[15]。因此根據條件，本研究內設置4種工況以對應不同流速的海洋環境，各個工況由代號Un表示(代號中n對應當前工況下的海流速度)。不同工況的設置信息及雷諾數計算由表1給出。

表1 工況設置信息表

儲氣裝置計算域模型如圖2所示，采用笛卡爾坐標系為參考坐標，X，Y，Z分別代表順流方向、高度方向和橫流方向，坐標原點在速度入口面底部中心處，3個方向上的計算與尺寸分別為X∈(0，20D)，Y∈(0，3D)，Z∈(-5D，5D)。模型底面中心坐標為(5D，0，0)，到出口邊界的距離為15D以降低出口邊界對尾跡流場的影響，同時保證能夠監測到更大范圍的繞流尾跡。經計算，阻塞比為4.64%。根據不同工況設置對應流速的速度入口，流域底面及裝置外表面均為固定壁面，其余上、側表面為滑移壁面以降低壁面邊界對內部流場的影響，壁面移動速度與相應工況下的海流速度一致，流域出口邊界條件為自由流出口，流量權重為1。

圖2 計算域模型

圖3所示為裝置計算域網格劃分，在劃分網格時，將裝置周圍及后方流域做加密處理，并在裝置壁面以及流域底面劃分邊界層網格，用以保證裝置壁面受力的準確監測。圖3a～圖3c分別為網格的+X，+Y，+Z視圖，圖3d、圖3e為裝置邊界層網格。

圖3 計算域網格設置

2 結果與分析

2.1 流動結構

通過監測流場的時均速度沿各個方向上的分布情況，可以較為準確地劃分流場內不同流動特征的特殊流場區域，從而完整地描述整體流場，再根據所得信息，對特殊流場區域進行重點分析。

經過無量綱化處理后，圖4a、圖4b分別表示在裝置半高處，流向方向時均速度分量u與豎直方向時均速度分量v沿流向的分布情況。

圖4表明，裝置在不同流速的海洋環境中所形成的尾跡流場特征基本一致，但在形成位置以及速度量值上有所區別。圖4a所示的流向方向速度分量u在到達裝置迎風面之前迅速下降至0，在繞過裝置后，在距離背風面較近的尾流區內出現負速度，表示在這個流場范圍內，繞過裝置的流體的流動方向與來流速度相反，朝著裝置背風面流動，該區域稱之為回流區，隨后，向流場遠處發展的流體速度逐漸增加至正值并向來流速度接近直至達到穩定，并且通過4種工況的速度曲線對比，可以看出在0.5 m/s流速下，海流繞過裝置后形成的回流區較其他3種工況范圍更大，并且達到穩定的位置有所滯后。

圖4 不同流速下時均速度沿流速方向分布示意圖

不同工況下的u特征信息在表2中給出。從表中數據可以得出，4種工況下的最大回流速度umax均在0.15 Un附近，回流區的長度隨速度的增加而減小，不論是在何種工況下，在裝置后方均存在回流現象且逐漸穩定。0.5 m/s流速下達到穩定的位置在更遠的10D處。

表2 不同工況下的u特征信息表

圖4c表示豎直方向的時均速度分量v沿流向的分布，從圖中可以看出，與u的分布情況類似，裝置在不同流速的海洋環境中該方向速度分量分布情況基本一致，但在特征出現位置以及速度量值上有所區別。海流在到達裝置迎風面前方時有沿豎直方向向上運動的趨勢，流過裝置側面后，在背風面附近v迅速降低至負值，這表明在此區域內(影響區)，大量的海水向裝置底部運動(具體成因將在2.2節中說明)，隨著與裝置距離的增加，v逐漸從負值上升并趨于穩定。不同工況下的v特征信息如表3所示。從表3中可以看出，雖然4種工況下的v分布曲線形態相似，均在距離裝置12倍直徑的位置附近達到穩定，穩定值近似于0，這表明在此區域外的流場變化與裝置之間的關系很小，但在0.5 m/s工況下，海水繞流后向下運動的影響范圍較長且強度較弱。表4列出了在不同工況下、位于裝置的不同高度處，向下運動的海流的最大合速度信息，從數據中可以看出，位于裝置半高處附近，各個工況下的該速度達到最大值，表示在此區域內，下行海流強度最大。

表4 不同工況下、不同高度處的v特征信息表

根據圖4b與表3的數據，在v影響區內，距離裝置中心1D的位置，于流場中面沿裝置豎直方向再次監測時均速度分量v，雖然在流場區域內的v峰值不在此處，但也可以以此判斷向下運動的海流在豎直方向上的分布情況，如圖5a所示。在影響區外，距離裝置中心2D的位置，于裝置半高處沿橫流方向監測橫流方向時均速度分量w的分布情況，如圖5b所示。

表3 不同工況下的v特征信息表

由圖5a可知，v的分布沿裝置高度呈現拋物線的分布情況，在U1.0，U3.0，U5.0工況下，繞流后海水向下運動的速度峰值出現在裝置中高處偏下的位置，并顯示出與來流速度之間的良好的線性關系。U0.5工況下，該速度峰值出現的位置在中高處偏上，這是由于在低速流場下，繞過裝置頂端的下行流體在來流方向向上位置偏后，在此截面上(X=1D，YOZ截面)流動強度較弱、速度分量v較小所造成的。

圖5b顯示出流方向時均速度分量w在橫流方向上的分布關于流場中面是反向對稱的，表明在裝置后方有反向旋轉的旋渦產生，由旋渦產生的切向速度約為來流速度的0.3倍，在數值上，w的大小受來流速度影響較小。

圖5 不同流速下時均速度沿其他方向分布示意圖

2.2 旋渦結構

根據對各個工況下的時均速度場分析，可以判斷海水流過裝置時的整體過程及尾跡流場的時均形態，以此推測出在裝置后方，存在具有一定規則的旋渦結構，形成回流、海水下行、與橫流速度反向對稱現象等局部流場特征。

以0.5 m/s流速的工況為例，圖6是在裝置半高處、背風面附近的時均流線圖。如圖所示，裝置后方的回流區本質上是一對反向旋渦作用的結果，這是由于流體繞過曲面時速度先升后減，壓力先減后升，在曲面后方某一位置產生的邊界層分離現象所導致的，流體質點在流場中面附近與來流反向運動，因而速度為負值，形成了回流區。

圖6 Y=0.5 L，XOZ截面時均流線圖

圖7a和圖7b分別表示在流場中面及X=1D，YOZ截面的時均流線圖。如圖所7示，由于自由端的結構是球形，在流場中面內以及與此截面平行的各個截面內自由端的投影為圓形，所以流體繞過自由端時同樣會發生邊界層分離現象，形成一個順時針方向旋轉的旋渦結構，在背風面附近速度向下。與此同時，在與YOZ截面平行的截面內，如圖7b所示，同樣存在一對反向對稱的旋渦結構，在中面附近速度向下，因此，豎直方向的時均速度分量v的負值是圖7a、圖7b中旋渦結構共同作用的結果，形成了方向向下的海流。

圖7 其他截面時均流線圖

圖8a和圖8b分別表示在X=1D，X=2D的YOZ截面內，Q-criterion顯示的渦量與速度矢量合成圖的時間快照序列。如圖8a所示，隨著以1/4T為間隔的時間變化，在該截面雖然能夠清楚地觀察到旋渦結構，但由于旋渦結構復雜難以直觀地發現其發展規律。但是旋渦結構出現的位置相對固定，分別在裝置半高處流場中面兩側，以及裝置底部流場中面兩側。

在圖8b中，距離裝置稍遠處，在該截面能夠清楚地觀察到旋渦結構，如果以t=t0為特征流場，包括向左發展的下行海流以及在下行海流兩側存在一對反向旋渦，則當x=5/4T時刻，該特征流場再次出現。在2000～3600 s的時間內，特征流場出現的次數及出現的時間點如表5所示，將相鄰2個時刻做差得到的周期分布圖如圖9所示。經過計算，該周期分布的平均值為95.33 s，即為流場旋渦脫落的主頻。圖8a所示結果沒有規律的原因是距離裝置較近，截面處在回流區及下行海流影響區內，該區域是旋渦形成及發展的區域，流場雜亂無章。

表5 特征流場出現時間表

圖8 Q-criterion與速度矢量合成圖的時間快照序列

圖9 特征流場出現的周期分布圖

2.3 受力特征

圖10a～圖10d分別展示了4種工況下的裝置受力特征，阻力與升力通過式(5)、式(6)換算為相應的力系數。

圖10 不同工況下的裝置力系數變化曲線圖

(5)

(6)

隨著流動的發展，雖然在微觀上仍有波動，但裝置受到的阻力與升力逐漸趨于穩定。通過監測入口流量與出口流量的平衡與殘差收斂保證計算結果的收斂性。經過統計可以得出各個工況下阻力系數與升力系數的均值，見表6。當雷諾數超過3×106時，在自模區范圍內力系數不再隨Re變化，因此力系數應當不隨工況的變化而變化，對各個工況下的力系數進行對比的意義是以這樣一種方式進行“重復性實驗”，從而再次驗證計算結果的準確性，結果表明，儲氣裝置的阻力系數約為0.45，升力系數約為0.60。

表6 各個工況下阻力系數與升力系數的均值表

由于渦街的存在，裝置在橫流方向上所受的流體力隨著時間而周期性變化，同樣稱之為升力，用lift-z表示，以U0.5工況的計算結果為例，圖11是裝置所受橫流方向升力的升力系數時程圖，升力系數在0值附近正負波動證明了裝置側向存在一對交替產生的旋渦。經傅里葉變換后，換算成斯特羅哈爾數-幅值的關系曲線如圖12所示。

圖11 U0.5工況下橫流方向升力系數時程圖

圖12 U0.5工況下斯特羅哈爾數-幅值曲線圖

斯特羅哈爾數是表征流動周期性的相似準則，對于周期性的非定常流動，斯特羅哈爾數可以表示為速度與特征頻率之間的關系，從而確定旋渦脫落的頻率。圖12中圖像的最大幅值對應的頻率就是旋渦脫落的主頻，斯特羅哈爾數計算結果為0.17，表明在此工況下，渦脫的周期約為100 s，與特征流場出現的周期基本一致(見小節2.2)。

3 結論

使用LES模型對裝置處于不同流速的海洋環境中所形成的尾跡流場及受力特征進行了數值模擬，通過對流場內各個方向上的時均速度分量(u,v,w)、速度矢量、渦量，以及流體力等物理量的監測，對海流通過裝置后的流場及裝置的受力情況進行了研究，主要研究結論如下：

(1)不同速度的海流通過裝置后形成的流場特征相似，特征速度量值與特殊流域的范圍與海流速度有關，低流速海洋環境下的繞流影響區較廣，各個特殊流域的位置均有滯后；

(2)繞流形成的旋渦具有明顯的三維結構特征，經儲氣裝置側向流過的海流在裝置背風面附近逐漸向底部扭曲形成旋渦結構，背風面附近區域的兩側旋渦關于流場中面對稱分布；

(3)繞過裝置球頂后形成的分離渦導致海水向下運動，在裝置半高處附近與側向向下扭曲的旋渦合流，所形成的下行海流是尾跡流場中豎直方向負速度流體出現的主要原因之一；

(4)裝置所受阻力與升力隨流動的發展趨于穩定，但這種穩定是一系列無規則波動的力疊加而成。通過統計得出裝置的阻力與升力系數分別為0.45和0.6，以此可以計算在任意流速海洋環境下裝置所受到的流體力；

(5)側向形成的三維旋渦在距離中心2倍直徑的位置附近開始發生脫落，以橫流方向升力系數與流場速度矢量、渦量的變化規律為依據，對旋渦脫落周期定性定量分析，結果表明該旋渦脫落的周期在100 s左右，會對裝置造成低頻擾動。