海面上下?lián)舯┝黠L場特性模擬及分析

尤韓煒，陳昌萍,3*

(1.廈門大學建筑與土木工程學院，福建 廈門 361005；2.廈門理工學院福建省風災害與風工程重點實驗室，福建 廈門 361024；3.廈門海洋職業(yè)技術學院，福建 廈門 361003)

下?lián)舯┝魇且环N強烈的下沉氣流，并在下沉氣流沖擊阻礙物時產生向周圍擴散的沖擊強風，下?lián)舯┝鞯恼麄€過程將威脅飛機、船舶與結構物的安全，其產生的災害影響越來越受到人們的關注.下?lián)舯┝魍ǔ０l(fā)生在雷暴期間，發(fā)生概率在60%至70%之間[1].已有學者[2-6]對下?lián)舯┝髯饔迷谄降孛?、山地、平屋面、高層建筑、輸電線塔等進行數(shù)值仿真模擬；另有學者[7-9]進行了大量的物理模擬試驗，得到與實際情況相似的結果.上述研究重點關注下?lián)舯┝鳟a生的直線風在陸地上形成的風場.2015年6月,“東方之星”客輪遭遇下?lián)舯┝髟谒蟹羀10]，針對該事件，段亞鵬等[11]通過數(shù)值模擬方法，使用中尺度區(qū)域氣象預報(ARPS)模式同化常規(guī)資料和事故發(fā)生地附近4部雷達資料，結合觀測結果從氣象學角度分析了下?lián)舯┝鬟^程中系統(tǒng)結構及發(fā)展的變化特點，但未涉及下?lián)舯┝髋c水域之間的相互作用.實際上，海洋約占地球表面的71%，海洋與大氣接觸范圍極廣，相互作用顯著，同樣存在雷暴[12-13].船舶在海上航行時，下沉氣流產生的瞬時強風可能會讓船上的工作人員判斷失誤，應對不當就可能造成船舶傾覆.但是由于海上觀測條件有限，目前雷暴與海洋相互作用研究并不多,針對下?lián)舯┝髯饔糜诤Ｑ笊蠈铀虻挠绊戇€未有針對性的探討.

本文結合氣-液兩相作用，利用CFD(computational fluid dynamics)技術，基于ANSYS Fluent，開展了下?lián)舯┝髯饔糜诤Ｑ笊蠈颖砻娴亩S模型數(shù)值仿真.根據數(shù)值仿真結果分析了下?lián)舯┝髯饔糜诤Ｑ笊蠈颖砻娴娘L場變化、風浪高度時程變化、特征時刻風速剖面和時均風速剖面，為實際應用提供一定的指導.

1 模型簡介

1.1 幾何模型

本文針對單個靜態(tài)下?lián)舯┝鏖_展研究，假定產生下?lián)舯┝髯饔脜^(qū)域為平靜海面.由于計算流域水平方向各向同性，取二維垂直面進行研究.模擬下沉氣流沖擊海洋表面過程的研究對象包含氣相和液相兩個部分，計算域在Kim等[14]的基礎上進行了修改，增加下墊面為液相的部分.下?lián)舯┝鞯南鲁翚饬鲾?shù)值模型可以使用沖擊射流模型模擬，已得到風洞試驗驗證[15].考慮重力作用，本文在計算域的近壓力出口處設置了一個坡度，令液相在重力作用下保持在計算域內，而氣相隨沖擊射流的發(fā)展，從兩側壓力出口溢出，以提高相對于氣相和液相共用一個壓力出口的求解收斂性.

本文假定下?lián)舯┝鞒跏忌淞髦睆紻為600 m，速度為29 m/s[14]，使用縮尺模型，設置幾何縮尺比例為1∶2 000，速度縮尺和時間縮尺采用Mason等[16]推薦的1∶3和1∶1 000，計算域尺寸與邊界類型如圖1所示，數(shù)值參數(shù)見表1.

圖1 計算域尺寸與邊界條件設置Fig.1 Computational domain size and boundary conditions settings

表1 參數(shù)設置

1.2 網格劃分

劃分結構化網格，在氣液兩相交界區(qū)域、沖擊射流入口下方這一重點關注區(qū)域，采用局部加密網格的方式，達到計算結果精細化的目的，合理分配計算資源.而在下墊面為陸地時要考慮的近壁面網格高度在本處不做過多考慮，劃分最小網格尺寸為0.016 m，網格劃分如圖2所示.

圖2 局部網格劃分方案示意圖Fig.2 Local meshing scheme

1.3 兩相流模型

本文研究對象為下沉氣流作用在海洋上層表面的風場整體，屬于氣-液兩相連續(xù)介質，可使用體積分數(shù)表達其分布關系.VOF(volume of fluid)模型[17]是一種多相流模型，可追蹤各相流體的體積分數(shù)，多應用于分層流、有自由表面的流動等.氣-液兩相問題正是屬于自由表面的流動，所以選用VOF模型.

VOF模型中引入了流體體積分數(shù)C，C表示單元內目標流體體積與單元體積的比值，可通過C的數(shù)值大小判斷流體的介質類型.根據單元內流量的變化可對C進行更新，同時通過兩相所占比例加權平均得到對應湍流模型封閉方程中的密度和黏度.由于質量守恒，單元內C的增長率等于流入單元的C的速度u的梯度，即：

(1)

此時，兩相流動問題可按單相流動問題方式進行求解，解得各相流體的速度、空間坐標等基本參數(shù).在ANSYS FLUENT中，即求解湍流模型.

1.4 湍流模型

雷諾應力模型(RSM)[18]是湍流模型中的一種，適用于自由剪切流、平均應變率陡然變換的流動、應變場復雜的流動.下沉氣流初始沖擊可視為自由剪切流；下沉氣流沖擊液相后向兩側繼續(xù)發(fā)展成為徑向氣流屬于平均應變率陡然變換的流動；下沉氣流引起不規(guī)則風浪，屬于應變場復雜的流動.因此本文的湍流模型選用雷諾應力模型.

雷諾應力模型的輸運方程針對其張量的各個分量的進行求解，根據雷諾應力的時間變化率與對流項的和等于湍流擴散項、分子擴散項、壓力產生項、浮力產生項、壓力應變項、耗散項、系統(tǒng)旋轉產生項與自定義源項的和，可表示為：

(2)

式中：t為時間；ρ為密度；ui為各個方向的平均速度；u′i為各個方向的脈動速度；xk(k=1,2,3)為空間坐標分量；p為壓力；δik為克羅內克符號；μ為動力粘度，β為熱膨脹系數(shù)；gi為各個方向的梯度向量；θ是輻射溫度；Ωk(k=1,2,3)為空間角速度分量；εjkm為單位完全反對稱張量.

2 結果與分析

2.1 速度云圖

下?lián)舯┝髯饔糜谄届o海洋上層的特征時刻速度云圖如圖3所示，由于該模型屬于對稱模型，這里取計算域的右半部分速度云圖進行分析.

圖3 下?lián)舯┝鞲魈卣鲿r刻速度云圖Fig.3 Velocity contour of the downburst at different instants of time

1) 由圖3(a)可知，當t=0.20 s時，此時刻為下沉氣流沖擊平靜海面之前，下沉氣流與原本存在的臨近氣體發(fā)生拖曳卷吸作用，在氣流前端兩側形成反向環(huán)形渦；

2) 由圖3(b)可知，當t=0.50 s時，下沉氣流到達海面，觀察到主渦也隨之下沉并貼近海面，下沉氣流向兩側發(fā)展形成徑向氣流，方向有所抬升，下沉氣流與兩側發(fā)展的徑向氣流交匯處速度值明顯小于原有下沉速度值，發(fā)展頂端速度值有明顯的放大；

3) 由圖3(c)可知，當t=0.65 s時，可以觀察到主渦開始脫離，側向發(fā)展的徑向氣流抬升較為明顯，氣流抬升發(fā)展的前端產生明顯的次生渦，徑向氣流速度值大小較前一時刻有所減小，但仍大于初始下沉速度值；

4) 由圖3(d)可知，當t=0.80 s時,主渦脫離并與次生渦一起向兩側移動，徑向氣流前端的最大速度值衰減至初始下沉氣流速度值以下；

5) 由圖3(e)可知，當t=1.50 s時，徑向氣流發(fā)展至邊緣附近，抬升趨勢減弱；

6) 由圖3(f)可知，當t=2.00 s時，徑向氣流發(fā)展至計算域外，整個風場趨于穩(wěn)定，徑向氣流速度值由沖擊中心至液相兩側逐漸降低，徑向氣流整體呈現(xiàn)小幅抬升趨勢.

2.2 風浪高度變化時程

下?lián)舯┝髯饔糜诤Ｑ笊蠈铀蛞痫L浪特征時刻的風浪高度分布如圖4所示，縱軸坐標Z/Djet表示液相所處高度的標準化處理結果，橫軸X/Djet表示液相所處水平位置的標準化處理結果，定義Z/Djet=0為液相表面所在高度.

圖4 下?lián)舯┝鞲魈卣鲿r刻風浪高度分布Fig.4 The height distribution of wind waves at different instants of time under the downburst

1) 當t=0.15 s時，下沉氣流開始發(fā)展，由于氣體壓縮與壓力傳遞，在下?lián)舯┝飨路街行某霈F(xiàn)微小變形；當t=0.30 s時，下沉氣流達到海洋上層表面，使得沖擊處海水表面下凹；

2) 當t=0.45 s時，此時下沉氣流往兩側發(fā)展為徑向氣流，沖擊處海洋表面進一步下凹，徑向氣流作用下凹處附近的海洋表面，形成風浪向兩側發(fā)展；

3) 當t=0.60 s時，沖擊處海洋表面下凹程度顯著增加，風浪高度加??；

4) 當t=1.50 s時，徑向氣流發(fā)展至液相兩側，在下?lián)舯┝鞑粩嘧饔孟?，沖擊處的海洋表面不斷下凹，兩側出現(xiàn)風浪特性，表現(xiàn)為許多高低長短不等的波；

5) 當t=2.00 s時，整體風浪趨于平穩(wěn)，由于沖擊處下沉氣流作用，致使沖擊處海洋表面下凹，下沉氣流向兩側發(fā)展為徑向氣流，隨著一定坡度的海洋平面向兩側流動，造成兩側海面出現(xiàn)高低長短不等的波，沿著兩側海平面總體為爬升趨勢，并在近海岸四分之一處海平面趨勢回歸平緩.

2.3 特征時刻徑向風速剖面

下?lián)舯┝髯饔糜诤Ｑ笊蠈拥挠绊憣ο蠖酁榇?，船舶受到波浪作用發(fā)生搖擺，加以強風影響，風浪聯(lián)合作用成為船舶穩(wěn)性在波浪運動中最重要的考慮因素，在此分析特征時刻的最不利徑向風速剖面.為探究高度、時間變化對不同位置徑向水平風速特性的影響，考慮到整個風場對稱分布，繪制了右半部分特征時刻的各個特征時刻、特征位置的徑向風速剖面圖，r表示距離下?lián)舯┝鳑_擊處中心的水平距離，V表示徑向氣流水平方向速度，數(shù)據作標準化處理，徑向風速剖面如圖5所示.

圖5 各特征時刻各徑向位置豎向風剖面Fig.5 Vertical wind profiles at each radial position at different instants of time

1) 由圖5(a)可知，當t=0.40 s時，下沉氣流產生的主環(huán)形渦接近海洋上層表面，下沉氣流沖擊海洋上層表面中心點處并開始向兩側發(fā)展，在1.0Djet處出現(xiàn)最大徑向風速值，接近下沉氣流速度，徑向風速值大小在水平方向1.0Djet兩邊分別遞減，上方存在主渦，所以在較高處出現(xiàn)與徑向氣流方向相反的徑向速度分量值；

2) 由圖5(b)可知，當t=0.50 s時，在水平方向2.0Djet處徑向最大速度值放大至下沉氣流速度值的1.2倍以上，出現(xiàn)徑向最大風速；

3) 由圖5(c)可知，當t=0.80 s時，在高度1.0～3.0Djet范圍內，同一高度處的徑向風速值隨著r的減小而減小，水平方向距離1.5Djet、2.0Djet、2.5Djet、3.0Djet各處最大徑向風速值接近初始下沉氣流風速值.

4) 由圖5(d)可知，當t=1.50 s時，整體徑向風速剖面趨于穩(wěn)定，各個水平位置處的徑向風速剖面表現(xiàn)出較強的相似性.

整體分析表明，徑向最大風速出現(xiàn)在下沉氣流沖擊液相表面的時刻附近,這一點與瞿偉廉等[3]研究下?lián)舯┝鳑_擊地面的數(shù)值模擬結果一致，并且最不利徑向風速值產生在下?lián)舯┝髯饔煤Ｑ笊蠈颖砻媲捌陔A段，從下?lián)舯┝鳑_擊中心1.0Djet處向遠端發(fā)展，隨時間推移最大徑向風速值逐漸降低，全局徑向風速剖面趨于穩(wěn)定，大小接近初始下沉氣流速度值.

2.4 時均風速剖面與液相流速剖面

為了探究水平位置徑向時均風速剖面的分布，觀測選取距離沖擊中心水平方向1.0Djet、1.5Djet、2.0Djet、2.5Djet、3.0Djet5處的豎向風剖面全時程平均數(shù)據，繪制時均風速剖面圖，如圖6所示，其中Vrad-mean表示時均徑向風速.進行數(shù)據標準化，可以觀察到各處風剖面的徑向風速最大值在較高位置處0.4附近，徑向水平距離1.0Djet處豎向風剖面較其他風剖面略有不同，原因是距離沖擊中心點較近，下沉氣流發(fā)展為兩側徑向氣流時還未發(fā)展完全，屬于兩者之間過渡階段，底部徑向風速值相對其他位置較小，而上部徑向風速值大于其他相對位置.徑向水平距離1.5Djet、2.0Djet、2.5Djet、3.0Djet處豎向風剖面形狀相近，區(qū)別在于各個位置處最大風速段呈現(xiàn)不同的區(qū)間大小，區(qū)間大小按照距離依次減少.整體而言，各個水平位置的徑向風剖面符合下?lián)舯┝髫Q向風剖面風速值沿高度由小到大再衰弱的特點.為了更好地說明本文數(shù)值模擬結果的可靠性，將水平徑向距離下?lián)舯┝鳑_擊中心1.5Djet處的豎向風剖面同文獻[8,19-20]對下?lián)舯┝鳑_擊壁面研究的結果進行比對，結果見圖7.其中Vmax為豎向風剖面中的速度最大值，Zmax為Vmax所在的高度.比對結果發(fā)現(xiàn)，在Zmax以下部分，本次數(shù)值模擬的結果與其他學者研究結果近似，同Vicroy[19]和Oseguera[20]結果基本吻合.但是在Zmax附近以及以上部分的結果存在偏差，具體表現(xiàn)為：Zmax大于其他學者研究結果；在Zmax以上部分速度衰減趨勢更快.對此的解釋為：由于本文研究的是下?lián)舯┝鳑_擊液面表面，在沖擊后液面產生了一定的下凹，這使得下沉氣流向兩側發(fā)展為徑向氣流時，徑向氣流會沿著液體表面發(fā)展而有一定的抬升，其表現(xiàn)出的Zmax相比于下?lián)舯┝鳑_擊壁面更高.標準化Z之后，Zmax增加使該豎向風剖面圖更加接近全局，從而表現(xiàn)出衰減的趨勢，也同樣驗證了該處的豎向風剖面符合下?lián)舯┝鞯撞匡L速快速發(fā)展，上部逐漸衰減的特點.

圖6 各徑向位置時均豎向風剖面Fig.6 Time-average vertical wind profiles at different radial positions

圖7 徑向位置1.5Djet處豎向風剖面與已有研究比對Fig.7 Comparisons of vertical wind profiles at 1.5Djet with different analytical models

提取了t=1.50 s時刻各特征位置處的氣、液相交界處徑向流速剖面，以確認液相在風場作用下是否發(fā)生相應的擾動，如圖8所示，可知，液相淺層受到下?lián)舯┝饔绊懏a生的徑向水流速度值遠小于徑向氣流的速度，隨著液相深度增加，液相徑向速度很快減小至0，說明液相深度足夠，不影響最終模擬結果.

圖8 各徑向位置處的氣、液相水平速度剖面Fig.8 Gas and liquid phase horizontal velocity profiles at different radial position

3 結 論

本文利用CFD數(shù)值模擬方法，完成了海洋上層受下?lián)舯┝髯饔玫姆抡?，具體分析了下?lián)舯┝髯饔糜诤Ｑ笊蠈訒r的風場變化、風浪高度時程變化、特征時刻的徑向風剖面和豎向時均風剖面的特性.仿真結果符合下?lián)舯┝鞯膶嶋H特征，主要結論如下：

1) 下?lián)舯┝髯饔糜诤Ｑ笊蠈铀驎r，下沉氣流與周圍空氣剪切形成環(huán)形渦，沖擊海洋上層表面后向兩側發(fā)展為徑向氣流，符合下?lián)舯┝靼l(fā)展特征，相對于下?lián)舯┝髯饔糜谄矫骊懙?，前期徑向氣流方向有所抬升，風場穩(wěn)定后徑向氣流中部仍表現(xiàn)為小幅抬升，至沖擊中心水平方向遠處逐漸平穩(wěn).

2) 沖擊中心附近的海洋上層表面隨下?lián)舯┝鞒掷m(xù)作用表面產生下凹，隨時間增加下凹加劇，產生的沖擊波向兩側傳遞發(fā)展，至1.50 s后趨于穩(wěn)定.早期以下沉氣流沖擊海面為主要影響，海洋上層表面波的傳遞規(guī)律明顯，隨后由于徑向氣流繼續(xù)發(fā)展作用于海洋上層表面產生風浪，穩(wěn)定后出現(xiàn)高低長短不一的波.沖擊中心處風浪所在水平高度沿兩側呈爬升趨勢，近岸處風浪所在水平高度趨于平緩.數(shù)值模擬直觀地展示了理想條件下下?lián)舯┝鲗ζ届o海面造成的影響，但是未考慮初始波浪的存在.根據波的疊加原理，下?lián)舯┝髯饔煤跏疾ɡ撕Ｃ嫘纬傻木植匡L浪波高將比作用平靜海面形成的風浪波高更大，后續(xù)研究需考慮初始波浪與風浪的疊加以及初始波浪與下?lián)舯┝黠L場的相互作用影響.

3) 早期徑向風速值將超過初始下沉氣流風速值，最大達到下沉氣流速度值的1.2倍，隨徑向氣流發(fā)展穩(wěn)定后，徑向風速值的最大值與下沉氣流風速值比值接近1.0，體現(xiàn)為下?lián)舯┝髯饔糜诤Ｑ笊蠈颖砻鏁r早期威脅性更大、更易聯(lián)合風浪使海上艦船發(fā)生傾覆，海上艦船應利用其特征防范與抵御下?lián)舯┝?

4) 各水平位置處的徑向時均風剖面符合下?lián)舯┝髫Q向風剖面風速值沿高度由小到大再衰弱的特點，由于下沉氣流沖擊海洋上層表面產生的下凹引導了徑向氣流發(fā)展向上方有所抬升，使得抬升處的徑向風剖面風速最大值在標準化高度0.4附近，相比下?lián)舯┝髯饔糜陉懙兀畲箫L速值所在高度有所增加.