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

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

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

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

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

1 實驗方案與研究方法

1.1 實驗裝置

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

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

圖1 實驗裝置示意圖Fig.1 Experimental setup

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

1.2 工況設計

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

(1)

(2)

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

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

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

1.3 參數(shù)計算方法

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

(3)

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

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

αe=uh/ve

(4)

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

經(jīng)典MTT模型理論[5]認為：對于單噴口浮力羽流，卷吸作用主要發(fā)生在羽干處；對于羽干同一高度zk處的羽流質(zhì)點，其垂向速度分量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的垂向變化梯度計算羽流的卷吸系數(shù)αe.由式(6)和(7)可得：

(8)

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

(9)

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

2 實驗結(jié)果與分析

2.1 卷吸過程

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

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

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

表2 羽干各層垂向時均速度分量的高斯曲線擬合(工況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隨標準化穿透深度分布Fig.4 αe for all cases versus normalized penetration

2.2 輸移質(zhì)量

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

(10)

變形可得：

(11)

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

(12)

代入并整理可得：

M(z)=πρWmb2

(13)

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

圖5 所有工況M/Mmax隨標準化穿透距離分布Fig.5 M/Mmax versus normalized penetration

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

Mmax=CMρaB3/4N-5/4

(14)

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

表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

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

Mmax=0.276ρaB3/4N-5/4

(15)

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

2.3 誤差分析

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

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

圖7 初始時刻通過固定點評估偽湍動的PIV圖Fig.7 The image of fixed particles for the initial moment collected by PIV to assess pseudo turbulence

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

3 結(jié)語

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

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

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

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

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