高速矩形水舌表面破碎的大渦模擬

吳 鵬，張 華

(華北電力大學水利與水電工程學院，北京 102206)

壩身泄洪孔或泄洪洞泄洪時的挑流水舌會發生霧化現象，其霧源有三個部分，分別是水舌擴散摻氣、水舌入水噴濺、水舌空中相碰而產生的霧化[1- 5]。其中第一霧源主要來自水舌表面破碎，雖然其霧源量較小，但第一霧源的特性與第二、第三霧源的特性密切相關，并且是第二與第三霧源的基礎。因此，高速矩形水舌的表面破碎機理是一個重要的科學問題。

對于水舌表面初次破碎機理這個問題，可追溯到小尺度自由射流的研究，對其的研究已取得如下的研究進展。施紅輝等人推廣了Ryhming模型，確定了噴嘴內部結構對脈沖噴流特性的影響[6]。金晗輝等人對三維小寬高比矩形噴嘴射流流場進行了數值研究，著重分析了流場中各典型截面上流場速度的分布特點、脈動速度分量的特征[7]。Abbas等人重點研究了主動與被動擾動同時作用下噴流的時空動力學[8]。Jaberi等人通過理論和實驗研究了矩形和橢圓液體射流的軸變換現象的振蕩波長和頻率[9]。Zhang等人通過實驗發現了受外部擾動的矩形噴嘴噴射出的微噴流的破裂特性[10]。李子豐等人使用OpenFOAM中的VOF-LPT耦合算法模擬了射流霧化過程，為更大尺度的模擬射流提供了可行的高效模擬方案[11]。

但在這些研究中，大部分忽略了重力因素或其影響微乎其微，而重力對于挑流霧化的大尺度水舌表面破碎的影響不可忽略。除了重力因素，渦結構也是最主要的影響因素之一。劉宣烈認為，入口處由于流速脈動等因素，使水舌圍繞其運動跡線擺動，促使渦產生，加劇渦之間的相互作用，加速了水舌的破碎[12]。吳持恭綜合了表面波破碎理論、邊界層發展理論以及紊動強度理論，認為水流內部不同尺度的渦體是導致水舌表面初次破碎的主要原因之一[13]。

為研究高速矩形水舌的初次破碎，采用流體體積法與大渦模擬法對其進行數值模擬計算。速度邊界使用渦流法，在入口處加入隨機點渦陣，凸顯隨機點渦陣作用條件下，矩形水舌的破碎和霧化特性。為進一步研究挑流水舌的后續破碎、霧化奠定基礎。

1 數學模型

1.1 基本假定

(1)水舌為不可壓縮性流體；

(2)不計初始風場速度。

1.2 水舌運動基本方程

研究水舌在氣相環境下的運動，滿足如下連續性方程：

(1)

動量方程：

(2)

式中，ui—速度張量，m/s；ρ—密度，kg/m3；p—壓強，Pa；ν—運動粘度，m2/s；fi—單位質量力，m/s2。

1.3 體積分數法與液滴轉換模型

流體體積的相方程為

(3)

式中,α—相體積分數。

為統計與分析破碎后液滴的分布，采用流體體積法-離散顆粒法的液滴轉換模型[14]。該方法結合了歐拉-歐拉法的VOF(Volume of Fluid)模型和歐拉-拉格朗日的DPM(Discrete Phase Model)模型，如果液體塊滿足設置條件，即當其體積小于體積標準或其球形度大于球形度標準時，則將其轉化為離散的液滴粒子，然后將這些粒子從連續相的歐拉體系中轉移到離散相的拉格朗體系中。該模型常用于模擬燃氣輪機、內燃機和其他類似應用中的液體霧化[15]。

(4)

式中,Q—球形度；Vp—液體塊體積，m3；Sp—液體塊表面積，m2。

1.4 大渦模擬模型

(5)

過濾操作后產生的亞格子模型主要捕捉小尺度渦，其亞格子尺度應力τij為

(6)

因為挑流水舌在射出時會產生渦結構，所以采用LES中的WMLES(Algebraic Wall-Modeled LES)模型[16]，其允許在高雷諾數下計算壁面有界流動，而無需像傳統LES那樣在高雷諾數下大幅提高網格分辨率。同時選擇WMLES中的S-Ω模型，該模型提供零渦流粘度，可以計算轉捩效應，在分裂計算層計算中不會產生過大的渦流粘度[17]。其亞格子尺度的渦流粘度μt為

(7)

式中,κ—von Kármán常數，取值為0.4187；CSmag—Smagorinsky常數，取值為0.2；dw—離固壁的垂直距離，m；Δ—修改后的網格比例尺；y+—垂直于固壁的內部比例；S—應變率。

1.5 速度邊界條件

挑流水舌在通過泄洪孔或泄洪洞射出的這個過程中，由于墻體與水流的相互作用，會產生強烈的擾動效果，導致在進入自由氣相時，渦結構其實已經存在，在后續的水舌破碎中發揮重要作用。所以該模型的速度邊界采用渦流法速度邊界，如圖1所示，在入口截面上給予隨機分布的點渦[18]。

圖1 隨機點渦位置示意圖

通過擾動渦度場(即垂直于流向方向的平面上的二維平面)在指定的平均速度剖面上添加擾動，該渦旋法是基于渦旋量的二維演化方程的拉格朗日形式和畢奧-薩伐爾定律[19]。

(8)

(9)

式中,ω—渦量大小，s-1；x—位置向量；xi—點渦的位置向量；η—渦旋點的空間分布；Гi—速度環量，m2/s；A—入口截面面積，m2；k—湍流動能，m2/s2；N—渦旋點的數量。

將這個隨機點渦陣的擾動速度場速度離散化后

(10)

(11)

式中,z—流向方向的單位向量;σ—渦旋點的大小，m；ε—平均湍流耗散率，m2/s3。

故瞬時速度為

(12)

在水舌入口邊界上，將隨機點渦陣所誘導的擾動速度，添加到指定的平均速度上，以實現渦結構對水舌速度邊界條件的影響。

2 水舌數值模型的計算

2.1 網格劃分與計算條件

基于有限體積法，采用對矩形水舌與空氣初始交界面加密后的蝶形網格，如圖2(a)所示，矩形入口為0.01m×0.01m，這種網格可以更好地刻畫水舌表面的發展變化，易于獲得重點關注區域的數據。

圖2 計算域網格

現設計5種入口平均速度不同的計算方案，其分別為20、35、50、65、80m/s，速度變化在1.0～0.8u，為凸顯渦結構的作用，取I=10%[20]。它們液體參數與計算工況如表1與表2所示。

2.2 模型驗證

為便于對比，采用圓射流的數值計算結果來驗證本文的數值模型。應用本文的數值模型對小尺度圓射流進行模擬，入口直徑為2mm，平均速度為120m/s，湍流強度為5%，表面的液相體積分數為0.2，時間為0.0003s，如圖3(a)所示。

表1 液體參數

表2 計算工況

而文獻21中，通過實驗拍攝得到的圓形射流如圖3(b)所示，其入口直徑為2mm，時均速度約為98m/s[21]。

通過計算結果和實驗結果的對比，在噴口附近，現有一段表面相對光滑的階段，沿流向的后一段表面出現鱗狀突起。兩者射流表面的斑圖基本一致，中心速度相差不大，說明本文的數值模型正確。

3 各工況計算結果與分析

使用本數值模型對上述計算工況1～5進行數值模擬，并對各工況結果中的水舌形變、摻氣、渦結構、液滴分布進行分析。

3.1 水舌形變

觀察計算工況1的液相體積分數為0.2的水舌，其整體形態如圖4(a)所示，根據形變可以將水舌分為如下三個階段。

初始階段：入口附近的水舌呈現光滑表面，隨后四條直角邊開始出現起伏的“波浪”形狀，并且形變逐漸加劇，成為“鋸齒”形狀，如圖4(b)所示。

過渡階段：表面的“鋸齒”加深，逐漸形變成“絲狀”形狀，存在少量的液滴脫落，如圖4(c)所示。

破碎階段：直角邊的“絲狀”從中間斷裂，破碎分離出大量的小液滴，如圖4(d)所示。

高速矩形水舌的徑向截面形狀會沿流向距離發生變化。以計算工況1為例，在水舌的初始階段，如圖5(a)和圖5(b)所示，其整體還能夠保持著接近矩形的形狀。在過渡階段，即圖5(c)和圖5(d)，原來的矩形邊界由直線逐漸形變成為波浪狀，波峰波谷分布逐漸不規則，波峰越來與凸出，波谷也越來越凹陷，其四個頂角的形變較大。進入破碎階段，即圖5(e)和圖5(f)，液滴大量脫落后，水舌截面邊界線條逐漸圓潤。由于重力的影響，截面整體形狀接近水滴形狀，上窄下寬，摻氣程度也是上小下大。在整個水舌表面破碎的過程中，縱向的水舌截面未產生長軸與短軸的變換，即軸變換現象[22]，說明表面張力產生的毛細管效應不明顯。

為了描述這種矩形水舌截面邊界的形變，本文提出如下無量綱形變關系式：

圖3 本文計算結果與文獻21的對比圖

圖4 t=0.03s，水舌的形變(水舌方向：自左向右)

(13)

式中,ζ為矩形“波浪”邊界的形變參數；S—徑向截面液相體積分數為0.98的等高線面積，m2；L—徑向截面液相體積分數為0.98的等高線周長，m。

應用式(13)可以計算得出，圓，ζ=π/4；而正方形，ζ=1，可見形變參數ζ反映了正方形比圓的形變要大。當水舌截面產生“波浪”型邊界時，其面積一定比相同周長下的正方形更小，ζ也隨之變大，即“波浪”波高越大，“波浪”波長越小，ζ的值也越大。故無量綱形變參數ζ能夠評價水舌截面的不規則形變。

5個計算工況下的形變參數ζ沿流向方向的變化如6所示。在計算工況1、2中，形變參數ζ都是先震蕩增大后震蕩減小，存在一個最大形變峰值，意味著水舌在運動過程中存在一個紊動到穩定的過程。而隨著速度的增加，雷諾數也增大，最大形變發生的位置向后移動，從x/d=5后移到x/d=50，并且最大形變程度增大，從ζ=1.4346增大到ζ=2.9282，同時震蕩強度也增強。

圖5 t=0.03s，徑向截面液相體積分數沿流向距離的變化

圖6 形變參數沿流向的變化

3.2 水舌摻氣

摻氣程度也是水舌表面初次霧化過程中很重要的組成部分，本文提出無量綱的摻氣面積比公式：

(15)

式中,C—摻氣面積比；S0.20—徑向截面中液相體積分數大于0.20的面積，m2；S0.98—徑向截面中液相體積分數大于0.98的面積，m2。

各工況下，摻氣面積比沿流向距離的變化如圖7所示。由于液流速度與氣相速度的差值大，兩相交界面發生強烈摩擦，水舌摻氣程度快速上升。其中，隨著速度的增加，摻氣程度的波動也越來越劇烈，峰值也隨之增大。

圖7 摻氣面積比沿流向的變化

3.3 渦結構與形變的關系

渦結構是影響水舌表面破碎的重要因素之一，可以通過渦量的大小來追蹤渦結構，如圖8所示。

在矩形水舌處于初始階段時，即圖8(a)至圖8(c)，渦結構主要集中在入口內表面與水舌的表面附近，并且渦量極大，而水舌內部也存在部分離散的渦結構。

隨著水舌流動，渦結構開始發展擴散，因為氣相的密度比液相的密度小，更加容易受到渦結構的影響，所以渦結構主要是從液相向氣相擴散，集中于水舌表面。從圖8(c)和圖8(d)可以明顯發現，在水舌的過渡階段中，形變部位附近渦結構集中，渦量較大，說明渦結構是導致水舌表面初次破碎的重要原因之一。

圖8 t=0.03s時，渦量大小沿流向的變化

在矩形水舌的破碎階段，即圖8(e)至圖8(f)，渦結構已經完全擴散至氣相中，部分的渦結構開始脫落，并且在重力的影響下，主要集中于水舌的兩側與上側。

3.4 液滴分布

使用液滴轉換模型后，分離的液滴被轉換為離散的粒子，圖9與圖10是其直徑與空間的分布。

圖9 液滴直徑分布

圖10 液滴的空間分布

通過液滴的直徑分布圖，如圖9所示，可以發現水舌表面在初次破碎后形成的小液滴的直徑粒徑范圍在0.3～2.6mm，是一個“單峰”集中分布，其中主要集中在1.0mm左右。隨著速度的增加，粒徑分布更加集中，粒徑占比的峰值也越大。

液滴在截面上的方向占比如圖10(a)所示，破碎而出的液滴在水舌截面的上下部分的占比差異較大，在0°～180°內的液滴數量占比遠遠大于180°～360°的數量占比。說明在重力的影響下，水舌上表面比下表面更容易破碎，飛濺出更多的液滴。并且水舌截面的上部分中，45°與135°左右的液滴數量占比最高，而下部分的225°與315°左右的數量占比較大。主要原因在于頂點位置的液相與氣相的接觸面積更大，所受擾動更大，再加上渦結構的發展，導致頂點更容易形變凸出后分離液滴。說明矩形水舌在破碎過程中，最主要的噴濺點位于上表面的兩個頂點附近的區域。

而圖10(b)是液滴在流向方向上的數量分布。流向距離越遠，分離出的液滴越多；速度越快，破碎距離越小，液滴也越早分離。

4 結論

本文針對高速矩形水舌的表面破碎過程做了數值模擬，得出以下結論：

(1)建立了高速矩形水舌運動的數學模型。在水舌入口處設置隨機點渦陣來體現水舌入口處的隨機擾動，模型能夠刻畫矩形水舌表面形變、表面破碎、液滴形成等過程。

(2)矩形水舌表面的初次破碎分為三個主要過程：初始階段、過渡階段、破碎階段。渦結構最開始主要集中于水舌的內外表面附近，從水舌表面向氣相發展、擴散、脫落，并在重力的影響下，集中于水舌的左右兩側以及上側。速度越大，越早破碎。

(3)在重力與擾動的作用下，矩形水舌的初次破碎位置主要位于上表面的兩條指直角邊附近，液滴粒徑為0.3～2.6mm，集中于1.0mm。入口平均速度越大，粒徑分布越集中。