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

吳 鵬，張 華

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

壩身泄洪孔或泄洪洞泄洪時(shí)的挑流水舌會(huì)發(fā)生霧化現(xiàn)象，其霧源有三個(gè)部分，分別是水舌擴(kuò)散摻氣、水舌入水噴濺、水舌空中相碰而產(chǎn)生的霧化[1- 5]。其中第一霧源主要來(lái)自水舌表面破碎，雖然其霧源量較小，但第一霧源的特性與第二、第三霧源的特性密切相關(guān)，并且是第二與第三霧源的基礎(chǔ)。因此，高速矩形水舌的表面破碎機(jī)理是一個(gè)重要的科學(xué)問(wèn)題。

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

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

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

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

1.1 基本假定

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

(2)不計(jì)初始風(fēng)場(chǎng)速度。

1.2 水舌運(yùn)動(dòng)基本方程

研究水舌在氣相環(huán)境下的運(yùn)動(dòng)，滿足如下連續(xù)性方程：

(1)

動(dòng)量方程：

(2)

式中，ui—速度張量，m/s；ρ—密度，kg/m3；p—壓強(qiáng)，Pa；ν—運(yùn)動(dòng)粘度，m2/s；fi—單位質(zhì)量力，m/s2。

1.3 體積分?jǐn)?shù)法與液滴轉(zhuǎn)換模型

流體體積的相方程為

(3)

式中,α—相體積分?jǐn)?shù)。

為統(tǒng)計(jì)與分析破碎后液滴的分布，采用流體體積法-離散顆粒法的液滴轉(zhuǎn)換模型[14]。該方法結(jié)合了歐拉-歐拉法的VOF(Volume of Fluid)模型和歐拉-拉格朗日的DPM(Discrete Phase Model)模型，如果液體塊滿足設(shè)置條件，即當(dāng)其體積小于體積標(biāo)準(zhǔn)或其球形度大于球形度標(biāo)準(zhǔn)時(shí)，則將其轉(zhuǎn)化為離散的液滴粒子，然后將這些粒子從連續(xù)相的歐拉體系中轉(zhuǎn)移到離散相的拉格朗體系中。該模型常用于模擬燃?xì)廨啓C(jī)、內(nèi)燃機(jī)和其他類似應(yīng)用中的液體霧化[15]。

(4)

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

1.4 大渦模擬模型

(5)

過(guò)濾操作后產(chǎn)生的亞格子模型主要捕捉小尺度渦，其亞格子尺度應(yīng)力τij為

(6)

因?yàn)樘袅魉嘣谏涑鰰r(shí)會(huì)產(chǎn)生渦結(jié)構(gòu)，所以采用LES中的WMLES(Algebraic Wall-Modeled LES)模型[16]，其允許在高雷諾數(shù)下計(jì)算壁面有界流動(dòng)，而無(wú)需像傳統(tǒng)LES那樣在高雷諾數(shù)下大幅提高網(wǎng)格分辨率。同時(shí)選擇WMLES中的S-Ω模型，該模型提供零渦流粘度，可以計(jì)算轉(zhuǎn)捩效應(yīng)，在分裂計(jì)算層計(jì)算中不會(huì)產(chǎn)生過(guò)大的渦流粘度[17]。其亞格子尺度的渦流粘度μt為

(7)

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

1.5 速度邊界條件

挑流水舌在通過(guò)泄洪孔或泄洪洞射出的這個(gè)過(guò)程中，由于墻體與水流的相互作用，會(huì)產(chǎn)生強(qiáng)烈的擾動(dòng)效果，導(dǎo)致在進(jìn)入自由氣相時(shí)，渦結(jié)構(gòu)其實(shí)已經(jīng)存在，在后續(xù)的水舌破碎中發(fā)揮重要作用。所以該模型的速度邊界采用渦流法速度邊界，如圖1所示，在入口截面上給予隨機(jī)分布的點(diǎn)渦[18]。

圖1 隨機(jī)點(diǎn)渦位置示意圖

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

(8)

(9)

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

將這個(gè)隨機(jī)點(diǎn)渦陣的擾動(dòng)速度場(chǎng)速度離散化后

(10)

(11)

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

故瞬時(shí)速度為

(12)

在水舌入口邊界上，將隨機(jī)點(diǎn)渦陣所誘導(dǎo)的擾動(dòng)速度，添加到指定的平均速度上，以實(shí)現(xiàn)渦結(jié)構(gòu)對(duì)水舌速度邊界條件的影響。

2 水舌數(shù)值模型的計(jì)算

2.1 網(wǎng)格劃分與計(jì)算條件

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

圖2 計(jì)算域網(wǎng)格

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

2.2 模型驗(yàn)證

為便于對(duì)比，采用圓射流的數(shù)值計(jì)算結(jié)果來(lái)驗(yàn)證本文的數(shù)值模型。應(yīng)用本文的數(shù)值模型對(duì)小尺度圓射流進(jìn)行模擬，入口直徑為2mm，平均速度為120m/s，湍流強(qiáng)度為5%，表面的液相體積分?jǐn)?shù)為0.2，時(shí)間為0.0003s，如圖3(a)所示。

表1 液體參數(shù)

表2 計(jì)算工況

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

通過(guò)計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果的對(duì)比，在噴口附近，現(xiàn)有一段表面相對(duì)光滑的階段，沿流向的后一段表面出現(xiàn)鱗狀突起。兩者射流表面的斑圖基本一致，中心速度相差不大，說(shuō)明本文的數(shù)值模型正確。

3 各工況計(jì)算結(jié)果與分析

使用本數(shù)值模型對(duì)上述計(jì)算工況1～5進(jìn)行數(shù)值模擬，并對(duì)各工況結(jié)果中的水舌形變、摻氣、渦結(jié)構(gòu)、液滴分布進(jìn)行分析。

3.1 水舌形變

觀察計(jì)算工況1的液相體積分?jǐn)?shù)為0.2的水舌，其整體形態(tài)如圖4(a)所示，根據(jù)形變可以將水舌分為如下三個(gè)階段。

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

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

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

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

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

圖3 本文計(jì)算結(jié)果與文獻(xiàn)21的對(duì)比圖

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

(13)

式中,ζ為矩形“波浪”邊界的形變參數(shù)；S—徑向截面液相體積分?jǐn)?shù)為0.98的等高線面積，m2；L—徑向截面液相體積分?jǐn)?shù)為0.98的等高線周長(zhǎng)，m。

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

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

圖5 t=0.03s，徑向截面液相體積分?jǐn)?shù)沿流向距離的變化

圖6 形變參數(shù)沿流向的變化

3.2 水舌摻氣

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

(15)

式中,C—摻氣面積比；S0.20—徑向截面中液相體積分?jǐn)?shù)大于0.20的面積，m2；S0.98—徑向截面中液相體積分?jǐn)?shù)大于0.98的面積，m2。

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

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

3.3 渦結(jié)構(gòu)與形變的關(guān)系

渦結(jié)構(gòu)是影響水舌表面破碎的重要因素之一，可以通過(guò)渦量的大小來(lái)追蹤渦結(jié)構(gòu)，如圖8所示。

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

隨著水舌流動(dòng)，渦結(jié)構(gòu)開(kāi)始發(fā)展擴(kuò)散，因?yàn)闅庀嗟拿芏缺纫合嗟拿芏刃。尤菀资艿綔u結(jié)構(gòu)的影響，所以渦結(jié)構(gòu)主要是從液相向氣相擴(kuò)散，集中于水舌表面。從圖8(c)和圖8(d)可以明顯發(fā)現(xiàn)，在水舌的過(guò)渡階段中，形變部位附近渦結(jié)構(gòu)集中，渦量較大，說(shuō)明渦結(jié)構(gòu)是導(dǎo)致水舌表面初次破碎的重要原因之一。

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

在矩形水舌的破碎階段，即圖8(e)至圖8(f)，渦結(jié)構(gòu)已經(jīng)完全擴(kuò)散至氣相中，部分的渦結(jié)構(gòu)開(kāi)始脫落，并且在重力的影響下，主要集中于水舌的兩側(cè)與上側(cè)。

3.4 液滴分布

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

圖9 液滴直徑分布

圖10 液滴的空間分布

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

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

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

4 結(jié)論

本文針對(duì)高速矩形水舌的表面破碎過(guò)程做了數(shù)值模擬，得出以下結(jié)論：

(1)建立了高速矩形水舌運(yùn)動(dòng)的數(shù)學(xué)模型。在水舌入口處設(shè)置隨機(jī)點(diǎn)渦陣來(lái)體現(xiàn)水舌入口處的隨機(jī)擾動(dòng)，模型能夠刻畫矩形水舌表面形變、表面破碎、液滴形成等過(guò)程。

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

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