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

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

引言

隨著當(dāng)前能源結(jié)構(gòu)的變革，可再生能源正發(fā)揮著前所未有的重要作用。但可再生能源固有的“間歇性”、“隨機(jī)性”和“低能量密度”的特點(diǎn)極大地阻礙了其高效可靠的利用[1-2]。這些由可再生能源產(chǎn)生的波動(dòng)的電能和波動(dòng)的電力負(fù)荷之間存在著嚴(yán)重不平衡。壓縮儲(chǔ)能技術(shù)是解決這種不平衡問題最有效的方法之一[3]，水下壓縮空氣儲(chǔ)能技術(shù)是一種適合在沿海地區(qū)、海島、海洋平臺(tái)等區(qū)域建設(shè)的，規(guī)?；目稍偕茉创鎯?chǔ)技術(shù)。與陸上壓縮空氣儲(chǔ)能系統(tǒng)不同之處在于，水下壓縮空氣儲(chǔ)能系統(tǒng)中的儲(chǔ)氣裝置需要放置在水下環(huán)境中，可以保證即使發(fā)生失效事件也不會(huì)對(duì)整個(gè)系統(tǒng)或者周圍設(shè)施與居民造成嚴(yán)重的危害[4-5]，這雖然在實(shí)現(xiàn)等壓存儲(chǔ)的基礎(chǔ)上提高了系統(tǒng)的安全性，但復(fù)雜多變的海洋環(huán)境也對(duì)水下儲(chǔ)氣裝置的設(shè)計(jì)提出了不小的挑戰(zhàn)。因此，研究?jī)?chǔ)氣裝置周圍的流場(chǎng)分別在普通與極端惡劣水下環(huán)境時(shí)的變化特征是十分有必要的。

CFD數(shù)值模擬方法被廣泛應(yīng)用于海洋工程結(jié)構(gòu)物的設(shè)計(jì)過程當(dāng)中，如水下立管、海上平臺(tái)樁腿等。在水下壓縮空氣儲(chǔ)能技術(shù)領(lǐng)域，數(shù)值模擬同樣是分析水下儲(chǔ)氣裝置周圍流場(chǎng)的有效研究方法之一。

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

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

1 數(shù)學(xué)模型

1.1 控制方程

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

(1)

式中，ρ—— 流體密度

p—— 壓力

ν—— 運(yùn)動(dòng)黏度

ui,uj—— 速度分量

xi,xj—— 位移分量

t—— 時(shí)間

經(jīng)濾波函數(shù)處理N-S方程，可以得到大渦模擬的控制方程為：

(2)

(3)

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

(4)

δij—— Kroneker符號(hào)

CS—— 常數(shù)，一般取0.1

Δ—— 網(wǎng)格濾波尺度

1.2 裝置的幾何模型與計(jì)算網(wǎng)格

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

圖1 儲(chǔ)氣裝置幾何模型與計(jì)算模型

在中國南海海域內(nèi)，根據(jù)當(dāng)?shù)厮臈l件及洋流活動(dòng)資料顯示，該海域內(nèi)海底水流速度全年在0.5～2.0 kn 范圍內(nèi)，冬季極端條件下海流速度可達(dá)2.8 kn以上。此外，水下溫度全年在8 ℃以下[15]。因此根據(jù)條件，本研究?jī)?nèi)設(shè)置4種工況以對(duì)應(yīng)不同流速的海洋環(huán)境，各個(gè)工況由代號(hào)Un表示(代號(hào)中n對(duì)應(yīng)當(dāng)前工況下的海流速度)。不同工況的設(shè)置信息及雷諾數(shù)計(jì)算由表1給出。

表1 工況設(shè)置信息表

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

圖2 計(jì)算域模型

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

圖3 計(jì)算域網(wǎng)格設(shè)置

2 結(jié)果與分析

2.1 流動(dòng)結(jié)構(gòu)

通過監(jiān)測(cè)流場(chǎng)的時(shí)均速度沿各個(gè)方向上的分布情況，可以較為準(zhǔn)確地劃分流場(chǎng)內(nèi)不同流動(dòng)特征的特殊流場(chǎng)區(qū)域，從而完整地描述整體流場(chǎng)，再根據(jù)所得信息，對(duì)特殊流場(chǎng)區(qū)域進(jìn)行重點(diǎn)分析。

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

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

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

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

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

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

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

根據(jù)圖4b與表3的數(shù)據(jù)，在v影響區(qū)內(nèi)，距離裝置中心1D的位置，于流場(chǎng)中面沿裝置豎直方向再次監(jiān)測(cè)時(shí)均速度分量v，雖然在流場(chǎng)區(qū)域內(nèi)的v峰值不在此處，但也可以以此判斷向下運(yùn)動(dòng)的海流在豎直方向上的分布情況，如圖5a所示。在影響區(qū)外，距離裝置中心2D的位置，于裝置半高處沿橫流方向監(jiān)測(cè)橫流方向時(shí)均速度分量w的分布情況，如圖5b所示。

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

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

圖5b顯示出流方向時(shí)均速度分量w在橫流方向上的分布關(guān)于流場(chǎng)中面是反向?qū)ΨQ的，表明在裝置后方有反向旋轉(zhuǎn)的旋渦產(chǎn)生，由旋渦產(chǎn)生的切向速度約為來流速度的0.3倍，在數(shù)值上，w的大小受來流速度影響較小。

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

2.2 旋渦結(jié)構(gòu)

根據(jù)對(duì)各個(gè)工況下的時(shí)均速度場(chǎng)分析，可以判斷海水流過裝置時(shí)的整體過程及尾跡流場(chǎng)的時(shí)均形態(tài)，以此推測(cè)出在裝置后方，存在具有一定規(guī)則的旋渦結(jié)構(gòu)，形成回流、海水下行、與橫流速度反向?qū)ΨQ現(xiàn)象等局部流場(chǎng)特征。

以0.5 m/s流速的工況為例，圖6是在裝置半高處、背風(fēng)面附近的時(shí)均流線圖。如圖所示，裝置后方的回流區(qū)本質(zhì)上是一對(duì)反向旋渦作用的結(jié)果，這是由于流體繞過曲面時(shí)速度先升后減，壓力先減后升，在曲面后方某一位置產(chǎn)生的邊界層分離現(xiàn)象所導(dǎo)致的，流體質(zhì)點(diǎn)在流場(chǎng)中面附近與來流反向運(yùn)動(dòng)，因而速度為負(fù)值，形成了回流區(qū)。

圖6 Y=0.5 L，XOZ截面時(shí)均流線圖

圖7a和圖7b分別表示在流場(chǎng)中面及X=1D，YOZ截面的時(shí)均流線圖。如圖所7示，由于自由端的結(jié)構(gòu)是球形，在流場(chǎng)中面內(nèi)以及與此截面平行的各個(gè)截面內(nèi)自由端的投影為圓形，所以流體繞過自由端時(shí)同樣會(huì)發(fā)生邊界層分離現(xiàn)象，形成一個(gè)順時(shí)針方向旋轉(zhuǎn)的旋渦結(jié)構(gòu)，在背風(fēng)面附近速度向下。與此同時(shí)，在與YOZ截面平行的截面內(nèi)，如圖7b所示，同樣存在一對(duì)反向?qū)ΨQ的旋渦結(jié)構(gòu)，在中面附近速度向下，因此，豎直方向的時(shí)均速度分量v的負(fù)值是圖7a、圖7b中旋渦結(jié)構(gòu)共同作用的結(jié)果，形成了方向向下的海流。

圖7 其他截面時(shí)均流線圖

圖8a和圖8b分別表示在X=1D，X=2D的YOZ截面內(nèi)，Q-criterion顯示的渦量與速度矢量合成圖的時(shí)間快照序列。如圖8a所示，隨著以1/4T為間隔的時(shí)間變化，在該截面雖然能夠清楚地觀察到旋渦結(jié)構(gòu)，但由于旋渦結(jié)構(gòu)復(fù)雜難以直觀地發(fā)現(xiàn)其發(fā)展規(guī)律。但是旋渦結(jié)構(gòu)出現(xiàn)的位置相對(duì)固定，分別在裝置半高處流場(chǎng)中面兩側(cè)，以及裝置底部流場(chǎng)中面兩側(cè)。

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

表5 特征流場(chǎng)出現(xiàn)時(shí)間表

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

圖9 特征流場(chǎng)出現(xiàn)的周期分布圖

2.3 受力特征

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

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

(5)

(6)

隨著流動(dòng)的發(fā)展，雖然在微觀上仍有波動(dòng)，但裝置受到的阻力與升力逐漸趨于穩(wěn)定。通過監(jiān)測(cè)入口流量與出口流量的平衡與殘差收斂保證計(jì)算結(jié)果的收斂性。經(jīng)過統(tǒng)計(jì)可以得出各個(gè)工況下阻力系數(shù)與升力系數(shù)的均值，見表6。當(dāng)雷諾數(shù)超過3×106時(shí)，在自模區(qū)范圍內(nèi)力系數(shù)不再隨Re變化，因此力系數(shù)應(yīng)當(dāng)不隨工況的變化而變化，對(duì)各個(gè)工況下的力系數(shù)進(jìn)行對(duì)比的意義是以這樣一種方式進(jìn)行“重復(fù)性實(shí)驗(yàn)”，從而再次驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確性，結(jié)果表明，儲(chǔ)氣裝置的阻力系數(shù)約為0.45，升力系數(shù)約為0.60。

表6 各個(gè)工況下阻力系數(shù)與升力系數(shù)的均值表

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

圖11 U0.5工況下橫流方向升力系數(shù)時(shí)程圖

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

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

3 結(jié)論

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

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

(2)繞流形成的旋渦具有明顯的三維結(jié)構(gòu)特征，經(jīng)儲(chǔ)氣裝置側(cè)向流過的海流在裝置背風(fēng)面附近逐漸向底部扭曲形成旋渦結(jié)構(gòu)，背風(fēng)面附近區(qū)域的兩側(cè)旋渦關(guān)于流場(chǎng)中面對(duì)稱分布；

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

(4)裝置所受阻力與升力隨流動(dòng)的發(fā)展趨于穩(wěn)定，但這種穩(wěn)定是一系列無規(guī)則波動(dòng)的力疊加而成。通過統(tǒng)計(jì)得出裝置的阻力與升力系數(shù)分別為0.45和0.6，以此可以計(jì)算在任意流速海洋環(huán)境下裝置所受到的流體力；

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