趙 亮， 張 巍， 賀治國,2， 談利明， 蔣后碩

(1. 浙江大學(xué) 海洋學(xué)院, 浙江 舟山 316000； 2. 國家海洋局第二海洋研究所, 杭州 310012； 3. 伍茲霍爾海洋研究所, 美國 伍茲霍爾 02543)

浮力羽流指在不同流體介質(zhì)中受浮力驅(qū)動(dòng)的柱狀流動(dòng)形態(tài)[1].在自然環(huán)境及工程應(yīng)用中普遍存在著浮力羽流運(yùn)動(dòng)形態(tài)，如火山噴發(fā)形成的氣流、火箭噴射出的高溫氣流、汽車尾氣、煙氣羽流、深海熱液羽流等.浮力羽流在運(yùn)動(dòng)發(fā)展過程中不斷夾帶周圍環(huán)境流體并發(fā)生物質(zhì)交換，引起的物質(zhì)輸移對(duì)周圍環(huán)境具有重大影響：火源上方形成的煙氣羽流在火災(zāi)科學(xué)中為一重要概念，其輸移質(zhì)量決定了煙氣生成量，是防排煙工程中的重要參考量[2]；深海熱液羽流運(yùn)動(dòng)過程中伴隨的質(zhì)量輸移對(duì)海洋的物質(zhì)循環(huán)和生物活動(dòng)有重要影響[3]，噴射時(shí)攜帶大量礦物粒子與金屬粒子，這些粒子少部分直接在噴口附近沉積，大部分則隨著羽流運(yùn)動(dòng)到一定高度并向周圍擴(kuò)散，參與大洋物質(zhì)循環(huán)[4].

國內(nèi)外對(duì)于浮力羽流已形成完備的理論體系，且有充足的數(shù)值模擬和物理模型實(shí)驗(yàn)作為支持[5-10].Morton等[5]在1956年提出經(jīng)典MTT模型(以Morton、Taylor和Turner首字母命名)，給出了層結(jié)環(huán)境中點(diǎn)源浮力羽流的流場解析解，并詳細(xì)分析羽干卷吸過程，不僅適用于熱液羽流，而且適用于氣相羽流等大多數(shù)流態(tài).魏文禮等[6]采用歐拉雙流體模型和混合模型結(jié)合標(biāo)準(zhǔn)k-ε紊流模型研究氣泡浮力羽流的動(dòng)力特性，認(rèn)為歐拉雙流體模型的結(jié)果更優(yōu).何標(biāo)等[7]采用大渦模擬對(duì)氣體羽流分層特性進(jìn)行數(shù)值模擬，發(fā)現(xiàn)密度較小的組分主要分散在空間上部80%的范圍內(nèi).Pham等[8]應(yīng)用二維及三維粒子圖像測速(PIV)技術(shù)研究純羽流，觀測到羽流的渦結(jié)構(gòu)，揭示了其在卷吸過程中的作用.Han等[9]同樣采用PIV技術(shù)測得伴隨反應(yīng)發(fā)生的羽流流場，并指出反應(yīng)放熱會(huì)降低羽流的卷吸速率.Mirajkar等[10]研究層結(jié)環(huán)境中浮力羽流的擴(kuò)散過程，并且建立了羽流橫向擴(kuò)散速度隨時(shí)間的變化關(guān)系.目前，對(duì)于層結(jié)環(huán)境中浮力羽流所引起的質(zhì)量輸移還有待深入研究，而密度層結(jié)是常見的羽流發(fā)生環(huán)境，尤其在海洋環(huán)境中.

本文通過層結(jié)鹽水發(fā)生裝置，在靜止的線性密度層結(jié)鹽水中開展一系列浮力羽流實(shí)驗(yàn).采用時(shí)間解析PIV技術(shù)獲得羽流的時(shí)空高分辨率流場，通過計(jì)算羽流的卷吸系數(shù)及垂向質(zhì)量通量，定量分析羽流在演化發(fā)展過程中所產(chǎn)生的質(zhì)量輸移過程；提出了羽流最大輸移質(zhì)量計(jì)算公式并通過實(shí)驗(yàn)數(shù)據(jù)進(jìn)行驗(yàn)證，為進(jìn)一步深入研究自然界及工程活動(dòng)中的羽流過程所產(chǎn)生的物質(zhì)輸運(yùn)與混合奠定了理論基礎(chǔ).

1 實(shí)驗(yàn)方案與研究方法

1.1 實(shí)驗(yàn)裝置

基于“雙缸法”[11]，采用課題組前期搭建的層結(jié)鹽水發(fā)生裝置[12]，在有機(jī)玻璃實(shí)驗(yàn)水槽(50 cm×50 cm×50 cm)中生成穩(wěn)定靜止的線性密度層結(jié)鹽水.實(shí)驗(yàn)時(shí)將噴口垂直向下固定在水槽正上方中心位置，并通過以高密度鹽水從頂部向下注入線性層結(jié)鹽水的辦法，進(jìn)行浮力羽流實(shí)驗(yàn).實(shí)驗(yàn)水槽中線性層結(jié)鹽水深度為47 cm，噴口位于水面下3 cm(圖1).實(shí)驗(yàn)時(shí)浮力羽流的運(yùn)動(dòng)方向雖然與實(shí)際情況相反，但兩者等效驅(qū)動(dòng)力與運(yùn)動(dòng)方向一致，即兩者動(dòng)力學(xué)過程相同，浮力羽流的射流方向并不會(huì)影響其流體動(dòng)力學(xué)過程[13-14].

應(yīng)用PIV獲得浮力羽流發(fā)展與演變過程中的可視化二維圖像及流場，解析流場的演變過程.PIV系統(tǒng)分為圖像采集與數(shù)據(jù)后處理系統(tǒng)兩部分，其中圖像采集系統(tǒng)由激光器和高頻相機(jī)組成.實(shí)驗(yàn)采用連續(xù)波(波長532 nm，激光源功率10 W)垂直照射實(shí)驗(yàn)水槽的側(cè)表面，形成厚度約3 mm的薄光面，使之經(jīng)過噴口中心線，即通過羽流的中心線.在實(shí)驗(yàn)水槽正前方布置高頻CMOS相機(jī)(采樣頻率200 幀/s)，相機(jī)采用16 mm F2.0 鏡頭，在垂直方向z-y平面28 cm(z)×21 cm(y)范圍內(nèi)獲取羽流發(fā)展的二維可視化圖像，坐標(biāo)及示意圖如圖1所示.采集的PIV圖像分辨率為 2 320 像素×1 726 像素，8 bit灰度圖；添加到噴口射流鹽水及密度層結(jié)鹽水的PIV粒子直徑15 μm，為與鹽水密度相近的尼龍微粒，質(zhì)量濃度3%(圖2).需要指出的是，本文浮力羽流在達(dá)到充分演化發(fā)展時(shí)尺度均不超過實(shí)驗(yàn)水箱尺寸的3/4，即羽流到水槽底邊界及邊壁的距離均超過實(shí)驗(yàn)水箱尺寸的1/4，因此邊界影響較小，可以忽略邊界效應(yīng)[8,14-15].采用PIV通用后處理軟件對(duì)羽流二維圖像進(jìn)行處理，算法為多通迭代法，診斷窗口尺寸為32像素×32像素，初始與最終診斷窗口的重疊區(qū)為50%.獲得流場數(shù)據(jù)后，再通過MATLAB進(jìn)一步分析數(shù)據(jù).

圖1 實(shí)驗(yàn)裝置示意圖Fig.1 Experimental setup

圖2 PIV采集的原始圖像(工況10)Fig.2 The initial image collected by PIV (Case 10)

1.2 工況設(shè)計(jì)

浮力頻率(N)定量表征環(huán)境的密度層結(jié)性，初始浮力通量(B)為噴口入射條件的度量，兩者是浮力羽流的基本參數(shù).計(jì)算表達(dá)式分別為

(1)

(2)

式中：g表示重力加速度；ρ0表示噴口出流鹽水密度；ρa(bǔ)表示與噴口處于相同高度的周圍環(huán)境鹽水密度；ρb表示水槽底部的環(huán)境鹽水密度；Q表示噴口出流鹽水的體積通量(噴口流量)；?ρ1/?z表示層結(jié)環(huán)境的密度梯度.表1所示為根據(jù)N、Q、ρ0和噴口直徑d等參數(shù)，設(shè)計(jì)的18組不同的實(shí)驗(yàn)工況.

表1 實(shí)驗(yàn)工況相關(guān)參數(shù)Tab.1 Experimental parameters

在實(shí)驗(yàn)水槽的壁面等較深位置處布置8個(gè)針狀采水器(分別距離底部0、6、12、18、24、30、36、42 cm).在線性密度層結(jié)形成后， 羽流實(shí)驗(yàn)開始前， 抽取少量水樣，確定羽流發(fā)生環(huán)境的鹽度分布，并根據(jù)狀態(tài)方程[16]計(jì)算鹽水密度.應(yīng)用最小二乘法對(duì)各組工況的鹽水密度隨水深分布進(jìn)行線性回歸分析，保證浮力羽流的發(fā)生環(huán)境具有良好的線性層結(jié)性[12].

1.3 參數(shù)計(jì)算方法

所有實(shí)驗(yàn)工況中，CMOS相機(jī)記錄的羽流演變發(fā)展持續(xù)時(shí)間約為50 s.羽流在發(fā)展過程中的瞬時(shí)流速具有隨機(jī)波動(dòng)性，已有研究[17]指出：當(dāng)羽流演變發(fā)展至某一時(shí)刻t*后，羽流流場便會(huì)達(dá)到穩(wěn)定階段，t*為與N相關(guān)的浮力時(shí)間尺度，且t*=2π/N.對(duì)于所有實(shí)驗(yàn)工況，在3t*≤t≤5t*穩(wěn)定階段，由羽流瞬時(shí)流場數(shù)據(jù)計(jì)算得到羽流時(shí)均速度場：

(3)

式中：U和W分別表示時(shí)均速度場水平(x軸)和垂直(z軸)方向速度分量；u和w分別表示通過PIV后處理得到的瞬時(shí)速度場水平和垂直方向速度分量；符號(hào)〈〉表示對(duì)所有瞬時(shí)速度場進(jìn)行時(shí)間平均處理.通過得到的羽流速度場計(jì)算卷吸率和輸移質(zhì)量，定量研究浮力羽流的質(zhì)量輸移過程.

1.3.1卷吸系數(shù) 卷吸率指一種流體在運(yùn)動(dòng)過程中吸入另一種流體的速率，通常用卷吸系數(shù)α表示.卷吸過程是羽流在演化發(fā)展過程中產(chǎn)生物質(zhì)輸移的根本原因.對(duì)于單噴口浮力羽流，卷吸系數(shù)的計(jì)算表達(dá)式為：

αe=uh/ve

(4)

式中：uh表示羽流邊界處的水平速度分量；ve表示同一高度處的羽流最大垂向速度分量.

經(jīng)典MTT模型理論[5]認(rèn)為：對(duì)于單噴口浮力羽流，卷吸作用主要發(fā)生在羽干處；對(duì)于羽干同一高度zk處的羽流質(zhì)點(diǎn)，其垂向速度分量W(x,zk)呈現(xiàn)高斯分布.

W=Wme-(x-xm)2/b2

(5)

式中：Wm表示位于(xm,zk)處羽干中線的最大垂向速度分量；b表示位于高度zk處的羽干半徑.

羽流在發(fā)展階段，羽干半徑b與羽流高度呈線性相關(guān)，即

b=cz+b0

(6)

式中：b0為截距常數(shù)；斜率c表示羽干半徑的垂向變化梯度-db/dz，與αe關(guān)系如下[18]：

αe=5c/6

(7)

利用b隨羽流穿透距離z的垂向變化梯度計(jì)算羽流的卷吸系數(shù)αe.由式(6)和(7)可得：

(8)

1.3.2垂向質(zhì)量通量 浮力羽流在對(duì)于海底成礦、海洋勘探等具有重要影響.其輸移質(zhì)量可通過垂向質(zhì)量通量M來表示為[17]

(9)

式中：r表示圓柱坐標(biāo)系中的羽流水平半徑，且最大值rmax=b；ρ表示羽流密度.

2 實(shí)驗(yàn)結(jié)果與分析

2.1 卷吸過程

為了便于顯示，將羽干平均分成11層，采用式(5)對(duì)每層羽流垂向時(shí)均速度分量進(jìn)行曲線擬合(圖3)，從而得到b，再通過有限差分法計(jì)算-db/dz，由式(8)計(jì)算得到αe.羽干各層垂向速度分量的擬合結(jié)果如表2所示，羽干的垂向時(shí)均速度分量符合高斯分布，表明本文實(shí)驗(yàn)滿足MTT模型卷吸理論中關(guān)于羽干垂向速度分布的假設(shè).

采用文獻(xiàn)[12]中的方法確定本實(shí)驗(yàn)中羽流的最大穿透距離Zmax，并將其作為無量綱化參數(shù)，對(duì)羽流垂向穿透深度z進(jìn)行標(biāo)準(zhǔn)化處理.所有工況的卷吸系數(shù)隨標(biāo)準(zhǔn)化穿透深度的平均分布曲線如圖4所示.αe隨著穿透距離明顯變化：從噴口原點(diǎn)至z≈-0.3Zmax，αe逐漸增加；z≈(-0.3～-0.5)Zmax處，αe在 0.11 附近穩(wěn)定波動(dòng)，該段深度對(duì)應(yīng)羽干區(qū)；之后，αe從z≈-0.5Zmax處的 0.1 減小至在z≈-0.7Zmax處的0，此處對(duì)應(yīng)中性浮力層[12]，即羽流主要向周圍擴(kuò)散，而卷吸過程相對(duì)較弱；在z≈(-0.7～-1)Zmax處，αe為負(fù)值，這段深度對(duì)應(yīng)羽帽區(qū)[17].本實(shí)驗(yàn)得到的羽干區(qū)卷吸系數(shù)為 0.09～0.13，與文獻(xiàn)[17]范圍相一致，表明了可靠性.

圖3 垂向時(shí)均速度分量的高斯曲線擬合(工況10)Fig.3 Gaussian fitting of mean vertical velocity (Case 10)

表2 羽干各層垂向時(shí)均速度分量的高斯曲線擬合(工況10)

Tab.2 Results of Gaussian fitting for mean vertical velocity component at each level (Case 10)

z/cmWm/(cm·s-1)xm/mmb/mm擬合殘差-0.85-1.460.504.356.20×10-7-1.70-3.690.924.752.69×10-6-2.55-5.841.015.563.27×10-6-3.40-5.781.116.652.17×10-6-4.25-5.381.207.956.83×10-6-5.10-4.771.269.042.17×10-5-5.95-3.961.369.905.38×10-5-6.80-3.931.4210.24.98×10-5-7.65-3.621.2010.36.98×10-5-8.50-2.871.6010.26.89×10-5-9.35-1.401.707.561.28×10-5

圖4 所有工況的αe隨標(biāo)準(zhǔn)化穿透深度分布Fig.4 αe for all cases versus normalized penetration

2.2 輸移質(zhì)量

將式(5)代入式(9)，并將直角坐標(biāo)系轉(zhuǎn)化為相應(yīng)的圓柱坐標(biāo)系，可得：

(10)

變形可得：

(11)

本文實(shí)驗(yàn)中噴口出流鹽水密度與層結(jié)環(huán)境的鹽水密度差異均在3%以內(nèi)，可認(rèn)為滿足Boussinesq假設(shè)，即密度變化對(duì)羽流發(fā)展的影響可忽略[5].假定在同一深度的羽流內(nèi)部密度恒定且為常數(shù)，因此式(11)中的密度ρ(z)只是z的函數(shù)，與r無關(guān).對(duì)式(11)進(jìn)行積分求解，

(12)

代入并整理可得：

M(z)=πρWmb2

(13)

式(13)與文獻(xiàn)[2]中氣體羽流水平截面的質(zhì)量流量表達(dá)式相同，垂向質(zhì)量通量與最大垂向速度分量和羽流半徑這兩個(gè)參數(shù)密切相關(guān).采用式(13)得到標(biāo)準(zhǔn)化垂向質(zhì)量通量隨標(biāo)準(zhǔn)化穿透深度分布，如圖5所示.羽流垂向質(zhì)量通量從噴口處開始迅速增加，到噴口附近z≈-0.1Zmax處增速變緩，再到z≈-0.65Zmax處達(dá)到最大值，此后逐漸減小直到接近于0.當(dāng)z=-0.65Zmax，該處屬于中性浮力層，這說明在中性浮力層以上，羽流由于卷吸作用不斷從周圍環(huán)境吸收質(zhì)量；在中性浮力層及以下，攜帶大量物質(zhì)的羽流向四周擴(kuò)散，同時(shí)將自身攜帶和環(huán)境吸收的物質(zhì)輸運(yùn)到更遠(yuǎn)處.

圖5 所有工況M/Mmax隨標(biāo)準(zhǔn)化穿透距離分布Fig.5 M/Mmax versus normalized penetration

根據(jù)π定理，對(duì)最大質(zhì)量通量Mmax進(jìn)行量綱分析，得到其函數(shù)表達(dá)式：

Mmax=CMρa(bǔ)B3/4N-5/4

(14)

式中：CM為待定常數(shù)，可由實(shí)驗(yàn)確定；所有工況的最大垂向質(zhì)量通量相關(guān)參數(shù)見表3.對(duì)所得數(shù)據(jù)進(jìn)行線性相關(guān)分析，結(jié)果如圖6所示，矩形符號(hào)表示原始數(shù)據(jù)，實(shí)線表示擬合直線.

表3 Mmax線性相關(guān)分析的相關(guān)參數(shù)Tab.3 Parameters with linear correlation for Mmax

圖6 Mmax線性相關(guān)分析結(jié)果Fig.6 Linear correlation for Mmax

通過最小二乘法計(jì)算得到Pearson相關(guān)系數(shù)，評(píng)估線性相關(guān)程度，用r表示，r值越大表示相關(guān)性越強(qiáng).所有工況數(shù)據(jù)擬合得到r=0.708，待定常數(shù)CM=0.276.通過假設(shè)檢驗(yàn)，計(jì)算P值為 1.01×10-3(P<0.01，樣本量n=18)，這表明Mmax/ρa(bǔ)與B3/4N-5/4有顯著的線性相關(guān)性.線性層結(jié)環(huán)境中的羽流Mmax由B、N和ρa(bǔ)共同決定：

Mmax=0.276ρa(bǔ)B3/4N-5/4

(15)

結(jié)合羽流卷吸系數(shù)分布、質(zhì)量通量分布及羽流的發(fā)展演化形態(tài)[12]可知：從噴口噴出的浮力羽流在羽干處卷吸入周圍環(huán)境流體后沿射流方向被輸移到更遠(yuǎn)處；在羽帽區(qū)，羽流與吸入的環(huán)境流體發(fā)生混合后，最終在中性浮力層被輸移到周圍，這是浮力羽流的物質(zhì)輸移過程.

2.3 誤差分析

本文層結(jié)環(huán)境中的浮力羽流是典型的變密度流動(dòng)，折射率會(huì)隨著流體局部密度的變化而改變.折射率的非定常會(huì)引起圖像中示蹤粒子的模糊，導(dǎo)致PIV測量誤差.折射率變化引起的速度脈動(dòng)現(xiàn)象被稱為偽湍動(dòng)，有必要評(píng)估其大小.

沿著中心線方向，在噴口下方2 cm處垂直固定一根直徑為2 mm的黑色細(xì)桿，并在桿表面標(biāo)出7個(gè)白點(diǎn)，如圖7所示.在相同的實(shí)驗(yàn)設(shè)置基礎(chǔ)上，不在層結(jié)鹽水及射流鹽水中添加PIV示蹤粒子，并重復(fù)實(shí)驗(yàn)工況10，得到層結(jié)環(huán)境中浮力羽流中心線7個(gè)固定點(diǎn)的瞬時(shí)速度時(shí)間序列.其中，位于羽干中部點(diǎn)B的瞬時(shí)速度隨時(shí)間波動(dòng)最為劇烈，其時(shí)間序列如圖8所示， 時(shí)均速度分量為UB=1.1×10-3m/s,WB=-2.0×10-3m/s.工況10的相同位置的垂向時(shí)均速度分量為W=6.2×10-2m/s，即由式(13)計(jì)算羽流最大輸移質(zhì)量，最大誤差約3%.因此，偽湍動(dòng)對(duì)于羽流的真實(shí)時(shí)均速度場影響有限，對(duì)于浮力羽流最大輸移質(zhì)量的影響可忽略.

圖7 初始時(shí)刻通過固定點(diǎn)評(píng)估偽湍動(dòng)的PIV圖Fig.7 The image of fixed particles for the initial moment collected by PIV to assess pseudo turbulence

圖8 折射率變化引起的羽流流場中固定點(diǎn)B的瞬時(shí)流速時(shí)間序列Fig.8 Measured time series of velocity components of fixed Particle B due to variations of refractive index

3 結(jié)語

在線性層結(jié)鹽水中開展浮力羽流實(shí)驗(yàn)，通過PIV技術(shù)獲得流場數(shù)據(jù)，分析其卷吸及質(zhì)量輸移過程，提出浮力羽流最大輸移質(zhì)量計(jì)算公式.

(1) 羽流位于相同高度處的垂向速度分量符合高斯分布.從噴口原點(diǎn)至z≈-0.3Zmax范圍內(nèi)，αe逐漸增加；在z≈(-0.3～-0.5)Zmax時(shí)，αe在 0.11 附近波動(dòng)；在z≈(-0.5～-0.7)Zmax時(shí)，αe逐漸減小至0；在z≈(-0.7～-1)Zmax時(shí)，αe為負(fù)值.

(2) 從噴口噴出的浮力羽流在羽干處卷吸入周圍環(huán)境流體后沿射流方向被輸移到更遠(yuǎn)處，在 -0.65Zmax處羽流垂向質(zhì)量通量達(dá)到最大值；在羽帽區(qū)，羽流與吸入的環(huán)境流體發(fā)生混合后，在中性浮力層被輸移到周圍，垂向質(zhì)量通量減小至0.Mmax由浮力通量、浮力頻率和與噴口處于相同高度的周圍環(huán)境密度共同決定.

(3) 折射率變化引起的偽湍動(dòng)對(duì)于羽流真實(shí)速度場的影響可忽略，其導(dǎo)致的浮力羽流最大輸移質(zhì)量計(jì)算誤差約為3%.

本文實(shí)驗(yàn)在靜止恒溫的鹽水環(huán)境中進(jìn)行，有待進(jìn)一步開展考慮側(cè)向流、溫度變化等環(huán)境條件下的羽流實(shí)驗(yàn).本文得到的最大輸移質(zhì)量計(jì)算公式及相關(guān)結(jié)論為系統(tǒng)研究羽流演化過程中產(chǎn)生的物質(zhì)混合與輸移提供了理論依據(jù).