江晨琦，劉 磊，2，李世舟

（1.上海交通大學a.海洋工程國家重點實驗室；b.船舶海洋與建筑工程學院，上海 200240；2.上海交通大學三亞崖州灣深?？萍佳芯吭海Ｄ先齺?572024）

0 引言

含顆粒的兩相羽流和射流常見于各種自然及工程流動中，如深海采礦作業(yè)中的尾礦排放、工廠含顆粒廢水排放、海洋疏浚作業(yè)、火山噴發(fā)等。兩相流動時顆粒會在環(huán)境中擴散，對生態(tài)環(huán)境產(chǎn)生不利影響［1-2］。了解羽流中顆粒的擴散速度、濃度等特性，對于預測和控制羽流排放具有重要意義。

目前，對于單相羽流及射流的理論及試驗研究已比較成熟［3］，但對于多相羽流尤其是復雜環(huán)境下（如橫流、分層流等）的多相羽流特性研究尚不深入。國內(nèi)外學者對兩相羽流開展了試驗研究。Jiang 等［4］使用三維粒子圖像測試（PIV）技術(shù)研究了垂直含沙羽流的顆粒速度特性，但對于含較細顆粒（小于100 μm）的羽流，圖像處理很難區(qū)分PIV粒子和羽流顆粒。Lee等［5］研究了水平含沙羽流特性，采用圖像局部顆粒計數(shù)的方法確定羽流中顆粒濃度。Hall 等［6］采用接觸式的光學探針測量垂直羽流顆粒的速度和濃度，初步研究顆粒速度和濃度特性、羽流擴散率受初始參數(shù)的影響以及顆粒動量通量的變化，但對于試驗尺度羽流，此類接觸式測量方法精度有限并會對流場造成影響。Mcconnochie 等［7］研究了染色流體中浮力方向?qū)τ鹆鲾U散率的影響。James等［8］通過染色對羽流中的顆粒和流體分別進行可視化，基于圖像結(jié)果研究了顆粒羽流在橫流中的軌跡及顆粒與羽流流體的分離現(xiàn)象，并與單相羽流結(jié)果進行對比，但這類方法成本較高，并且難以獲得羽流內(nèi)部的顆粒速度特性。

設(shè)計了模擬固液兩相羽流排放及擴散的試驗裝置，分析了基于高速圖像采集和計算機圖像處理的羽流擴散特性，并提出了顆粒濃度及速度的測試分析方法，用于靜止及流動條件下羽流特性的準確提取。

1 試驗裝置和試驗流程

1.1 試驗裝置

羽流排放及擴散模擬試驗裝置如圖1 所示，主要由羽流發(fā)生裝置、輸送管路系統(tǒng)、標定池及試驗水槽組成。羽流發(fā)生裝置包括直徑0.6 m、高0.4 m 的圓柱形玻璃桶和雙槳葉攪拌裝置，頂部設(shè)有注水口及顆粒投放漏斗。羽流通過管路系統(tǒng)被輸送至試驗水槽，管路設(shè)有流量計，用于實時測量排出流量。標定池用于確定羽流的初始濃度。

圖1 試驗裝置

試驗水槽長2.5 m、寬0.6 m、高1.0 m，為全透明玻璃材質(zhì)。水槽上部設(shè)有拖曳系統(tǒng)，可拖曳羽流排放口做勻速直線運動，以模擬均勻橫流中的羽流排放［8-9］。試驗中，可根據(jù)試驗目的改變排放深度、角度、拖曳速度以及排放口的直徑、形狀，研究不同排放條件下的羽流擴散特性。

試驗中采用光學無接觸式測量分析方法，圖像數(shù)據(jù)使用SONYa7M3 4K 全畫幅微單數(shù)碼相機采集，分辨率為2 160 ×3 840；進一步地，為捕捉局部顆粒的運動特征，采用Photron FASTCAM 系列高速攝像機進行圖像采集，配備105 mm 微距鏡頭，分辨率為1 024 ×1 024。

試驗顆粒可根據(jù)試驗需求進行選取。本試驗采用密度約為2 500 kg／m3的白色玻璃微珠，粒徑范圍為100～150 目（目指每英寸篩網(wǎng)上的孔眼數(shù)目），如圖2所示。

圖2 模擬羽流顆粒

1.2 試驗流程

試驗流程如圖3 所示。向試驗水槽中注水至水面沒過排放口5 cm。在發(fā)生裝置中加入適量顆粒和水，開啟攪拌裝置至穩(wěn)定轉(zhuǎn)速，持續(xù)攪拌至混合物均勻、穩(wěn)定。打開管路閥門，啟動輸送泵，待管路中羽流流量和濃度穩(wěn)定后，將混合物引入標定池，標定排放羽流濃度和流量，并與流量計測量結(jié)果對比校驗。標定完成后，關(guān)閉標定池閥門，打開羽流出口閥門，將混合物排入試驗水槽，開展試驗。排放達到穩(wěn)定狀態(tài)后，采用數(shù)碼相機或高速相機采集羽流圖像數(shù)據(jù)，流量計記錄流量數(shù)據(jù)。對于每個工況，進行數(shù)次重復試驗以驗證重復性。待試驗結(jié)束后，基于圖像處理方法對采集的圖像數(shù)據(jù)進行處理，得到羽流中顆粒的擴散及運動特性。

圖3 試驗流程

2 測試分析方法

基于非接觸式的光學圖像處理方法，以靜止環(huán)境中向下排放的垂直顆粒羽流為例，闡述羽流擴散特性、顆粒速度和顆粒濃度的圖像分析方法。圖像處理算法均通過Matlab軟件實現(xiàn)。

2.1 羽流輪廓

羽流輪廓是基于時均圖像灰度閾值的方法進行的。以典型試驗工況（羽流出口管徑D＝8 mm、平均初始流速u0＝1.92 m／s、初始質(zhì)量濃度c0＝15 g／L）為例，首先對攝像機拍攝得到的連續(xù)圖像進行灰度值時間平均，然后將計算獲得的圖像去除背景灰度，最后通過二值化和邊緣檢測算法得到羽流的視覺邊緣及由視覺邊緣包絡(luò)的羽流視覺半寬b（即羽流在某一高度處的半徑），如圖4 所示。

圖4 顆粒羽流形態(tài)示意圖

在羽流分析中，較常使用的是速度高斯半寬bp和濃度高斯半寬bc。在某一羽流高度z上，bp定義為流場速度uf＝（1／e）ufm上的點與中心軸線的距離，bc定義為顆粒濃度c＝（1／e）cm上的點與中心軸線的距離，其中ufm和cm分別表示羽流中心軸線上的流速和顆粒濃度。速度高斯半寬bp與視覺半寬b的關(guān)系可以表示為：2 種高斯半寬的比值λ ＝bc／bp通常取1.2［10］。理論上，在靜止環(huán)境下羽流視覺半寬隨高度線性增加［10］。因此，根據(jù)圖像各高度處羽流邊緣的像素點坐標進行線性擬合，得到羽流的視覺半寬擴散率αb＝db／dz，進一步按照比例關(guān)系得到速度高斯半寬擴散率αp＝dbp／dz。

圖像處理過程及結(jié)果如圖5 所示，其中x為羽流視覺邊界與中心軸線的距離。

圖5 羽流輪廓提取圖像處理

2.2 羽流濃度

為確定羽流初始濃度，將均勻、穩(wěn)定的混合物注入標定池內(nèi)進行標定。標定池出口下方設(shè)置長方體水箱以及孔徑遠小于試驗用顆粒的薄濾網(wǎng)。將混合物注入標定池數(shù)秒后關(guān)閉，記錄水箱中混合物總體積；取出濾網(wǎng)中的顆粒進行清洗、烘干、稱重，獲得顆?？傎|(zhì)量，進而獲得混合物的初始質(zhì)量濃度。

利用我院自制的調(diào)查問卷對兩組孕婦進行調(diào)查，主要包括：（1）孕婦的基本信息，如：年齡、民族、身高和體重，同時還包括文化背景、職業(yè)、孕產(chǎn)次和收入水平等。（2）孕前半年至孕前3個月的日常生活、工作、居住環(huán)境，同時還包括是否發(fā)生負性生活事件和感冒[3]。孕前期是否發(fā)生以下感染性疾病，如：慢性疾病和生殖系統(tǒng)感染，還需對孕婦的藥物服用史進行和飲食狀況進行掌握。營養(yǎng)攝入為：乳制品、蔬菜、肉、水果，在此基礎(chǔ)上還包括葉酸和維生素。（3）對圍產(chǎn)兒的具體狀況進行了解，如：胎兒狀況、分娩情況和診斷結(jié)果等[4]。

羽流濃度基于時均圖像進行測試分析。Mingotti等［11］在對含顆粒的固液兩相流濃度進行圖像處理時指出，在相同視場寬度下，顆?？倲?shù)應與光強（圖像灰度）線性相關(guān)。對于三維圓形羽流，在不同高度和不同水平位置羽流顆粒所占的厚度均不相同。由于采集的圖像僅呈現(xiàn)二維信息，因此可將圖像中二維信息依據(jù)高斯分布假設(shè)進行轉(zhuǎn)換，獲得該視場中的羽流中心濃度或平均濃度，如圖6 所示。以羽流邊緣延長線與中心軸線的交點為原點O，z軸與羽流出口方向平行，向下為正；x、y為水平坐標，x、y軸分別與相機視場方向垂直和平行。cm0表示x＝x1、y＝0 處的顆粒濃度。

圖6 濃度圖像處理信息轉(zhuǎn)換

根據(jù)三維多相羽流的高斯?jié)舛确植技僭O(shè)［12］：

式中，r為該點距離羽流中心的位置。由于本試驗采用的是120～250 μm的細顆粒，在質(zhì)量濃度小于20 g／L條件下重疊不明顯，因此可以認為時均圖像上某點的灰度值與該處視場中的顆?？倲?shù)正相關(guān)。根據(jù)高斯分布函數(shù)的性質(zhì)，在任意高度z上，有

根據(jù)式（2），對于圖像中顆粒相范圍內(nèi)的任意一點，該點所在高度的羽流視覺寬度2b范圍內(nèi)的視場平均濃度與該位置中心最大濃度cm之比為一定值A(chǔ)，故將圖像各點的灰度值除以該高度處的羽流視覺寬度作為該處的相對粒子濃度強度，從而得到全局相對濃度分布，如圖7 所示。根據(jù)標定得到的羽流出口濃度c0，可以得到全局的顆粒濃度分布。

圖7 羽流濃度圖像處理結(jié)果

濃度處理結(jié)果通過Lai等提出的顆粒羽流微分方程模型［13］并結(jié)合抽吸法［14］進行驗證。抽吸法是一種精度較高的接觸式濃度測量方法，如圖8 所示。將內(nèi)徑為5 mm、外徑為7 mm的透明抽吸管伸入待取樣位置，在羽流試驗進行的過程中同步抽吸到燒杯中，通過后續(xù)質(zhì)量分析得到濃度。Lee［15］的研究表明，當抽吸管角度與流速平行、速度比（即抽吸速度us與流速uc的比值）q＝us／uc＝1 時效果最佳，且當0.8 ＜q＜2.0時，測量誤差均在10%以下。垂直顆粒羽流驗證試驗中，抽吸管與羽流排放方向平行布置，蠕動泵抽吸速度的選擇根據(jù)羽流在該處的流場時均理論速度值確定，使得q接近于1，以保證抽吸效率。

圖8 抽吸法測顆粒濃度

以表1 中的工況3 為例，羽流充分發(fā)展段圖像處理得到的相對濃度分布結(jié)果如圖9 所示?？梢园l(fā)現(xiàn)，在羽流中心軸線上，圖像處理得到的濃度分布與微分方程模型結(jié)果、抽吸法單點測量濃度（z／D＝10，20）均吻合良好；在各水平高度位置，圖像處理結(jié)果均與高斯分布擬合結(jié)果吻合良好，同時擬合得到的高斯半寬值與基于視覺半寬確定的高斯半寬值具有良好的一致性，如圖10 所示。

表1 試驗工況參數(shù)

圖9 羽流濃度圖像處理結(jié)果及驗證

圖10 2種方法所得速度高斯半寬值對比

2.3 顆粒速度

羽流中顆粒的運動測量是基于高速相機采集和改進的粒子追蹤測速（PTV）圖像處理算法。由于羽流中顆粒速度的三維高斯分布及紊動特性，通過視場中特定位置處的顆粒速度并不穩(wěn)定，而且整體向量場較為稀疏，因此應收集一定區(qū)域內(nèi)的所有速度值來確定顆粒速度。將高速相機采集的結(jié)果逐幀提取，按8 幀為一組進行后續(xù)處理，處理流程和結(jié)果如圖11 所示。具體過程為：①將灰度圖像使用灰度閾值去除景深外的虛化顆粒，并通過濾波減小光照不均及顆粒反射對顆粒識別的影響；②進行自適應二值化；③較低濃度下對二值圖像直接進行連通圖提取，較高濃度下使用基于灰度圖像的分水嶺算法對重疊顆粒進行分離；④根據(jù)連通圖所占像素點數(shù)，按一定閾值篩選符合條件的連通圖作為顆粒，并根據(jù)連通圖像素點平均坐標計算顆粒中心坐標；⑤對相鄰幀使用相關(guān)系數(shù)法進行顆粒匹配；⑥若某一顆粒在8 幀內(nèi)能夠連續(xù)捕捉，則利用標定像素長度與位移、時間關(guān)系，計算顆粒通過該區(qū)域的速度值；⑦記錄一定時間內(nèi)通過該區(qū)域的所有顆粒的速度值。

圖11 圖像速度處理過程

對羽流勢核段外（充分發(fā)展段）的局部區(qū)域進行拍攝，監(jiān)測范圍為3.35 cm×3.35 cm區(qū)域。部分結(jié)果如圖12 所示，其中f為顆粒速度的柱狀頻數(shù)，uz為顆粒垂向速度，顆粒的平均沉降速度為0.41 m／s。

圖12 羽流充分發(fā)展段顆粒速度圖像處理結(jié)果

進一步地，在羽流的邊緣處，通過對單個顆粒的持續(xù)追蹤，可以確定顆粒在一定時間內(nèi)的運動軌跡和速度，如圖13 所示。據(jù)此，可對羽流的卷吸現(xiàn)象［3］開展進一步研究。

圖13 羽流邊緣顆粒運動檢測（相鄰幀時間Δt ＝0.8 ms）

基于上述顆粒速度的測試分析方法，在試驗水槽內(nèi)對200～250 mm 顆粒低濃度自由沉降過程進行分析，穩(wěn)定段顆粒速度分布如圖14 所示。可以發(fā)現(xiàn)，97%顆粒的速度范圍為32～42 mm／s，均值為36 mm／s。作為對比，使用傳統(tǒng)沉降法測得顆粒速度范圍在30～45 mm／s，均值為39 mm／s；使用Stokes公式估算相應顆粒范圍內(nèi)的沉降速度區(qū)間為31～48 mm／s?？梢哉J為，所提出的圖像處理算法可以較好地分析顆粒的運動特征。

圖14 顆粒自由沉降速度圖像處理結(jié)果

需要說明的是，該速度圖像處理方法不適用于顆粒濃度較大的羽流（質(zhì)量濃度＞20 g／L），此類羽流顆粒重疊情況過多，難以使用算法對顆粒進行識別、分離、匹配等操作。

2.4 測量誤差分析

試驗中產(chǎn)生的誤差主要來源于試驗裝置本身以及圖像處理算法。

（1）羽流發(fā)生裝置輸送的混合物由于管道淤積等原因造成羽流初始濃度和流速存在小幅波動，影響羽流初始參數(shù)的精確確定。按照流程進行試驗，先將混合物排放入標定池一段時間，再排入試驗水槽中。采用質(zhì)量分析方法，每隔5 s 記錄一次標定池管口或試驗水槽管口的質(zhì)量濃度值，如圖15 所示。結(jié)果表明，混合物經(jīng)過管道輸送需要一段時間（約10 s）才能達到穩(wěn)定狀態(tài)，同時標定池與試驗水槽出口質(zhì)量濃度存在一定偏差，重復測試結(jié)果表明該偏差值均在6%以內(nèi)。因此，在試驗開始前應盡量清除管內(nèi)的殘余顆粒，并且需要羽流排放10 s以上達到穩(wěn)定狀態(tài)后再進行數(shù)據(jù)采集工作。該發(fā)生裝置通過調(diào)整輸送泵的功率（0～150 W）控制流量，隨著試驗的進行混合物體積流量均呈現(xiàn)先升高再緩慢降低的趨勢，如圖16 所示。一般輸送泵運行120 s 后，體積流量和顆粒質(zhì)量濃度比初始穩(wěn)定值降低7%以上，因此單次試驗持續(xù)時間不宜過長。

圖15 出口顆粒質(zhì)量濃度變化

圖16 不同泵功率下混合物體積流量曲線

（2）對于排放速度較大的工況，羽流觸底后會揚起底部的部分顆粒，影響顆粒速度和濃度的圖像分析結(jié)果。因此，每次試驗前盡可能去除底部沉積的顆粒，并保證試驗水槽中無其他懸浮雜質(zhì)。

（3）濃度處理過程中，光照不均勻會顯著影響濃度結(jié)果，因此試驗過程中應盡量減少外界光源的影響，并采用均一無頻閃強光源照射試驗測試區(qū)域。

（4）速度處理過程中，算法對于較高濃度下重疊顆粒的處理存在誤差，因此對于每一個工況應保證持續(xù)時間足夠長，從而獲得足夠多的正確顆粒速度樣本。

（5）在算法的二值化以及連通圖提取步驟中，部分閾值需要根據(jù)顆粒大小以及光照強度確定。Muste等［16］在PTV 算法中也是根據(jù)經(jīng)驗選擇二值化閾值。閾值選擇不當會導致顆粒識別出現(xiàn)遺漏或?qū)⒃朦c識別為顆粒，因此試驗參數(shù)及環(huán)境發(fā)生變化后，需要調(diào)整閾值以減少在該條件下試驗產(chǎn)生的誤差。

3 試驗結(jié)果分析

3.1 試驗工況

一般地，羽流數(shù)Np和弗勞德數(shù)Fr是影響顆粒羽流中顆粒相擴散率的重要參數(shù)［12］，分別表示為：

式中：Ws為顆粒的特征沉降速度；φ0為顆粒的體積分數(shù)；g′＝g（ρs－ρf）／ρf，其中ρs和ρf分別為顆粒和水的密度。由式（3）、（4）可知，u0同時控制Np和Fr這2 個無因次數(shù)。因此，u0對顆粒羽流特性影響更為顯著?？刂艱不變，通過改變u0及c0來調(diào)整羽流初始狀態(tài)，主要試驗工況如表1 所示。

3.2 試驗結(jié)果

Lai等［13］的試驗結(jié)果表明，Np越大，羽流擴散率越大。針對工況1～6 開展試驗，在其他條件不變的情況下通過改變初始流速改變Np，得到各工況下速度高斯半寬擴散率αp，如圖17 所示。

圖17 速度高斯半寬擴散率αp與羽流數(shù)Np關(guān)系

通過時均結(jié)果和瞬時結(jié)果對比可以發(fā)現(xiàn)，羽流存在局部湍流擴散。對于時均圖像，顆粒相邊界近似擬合為一條直線。瞬時圖像與該擬合邊界相比存在差異，如圖18 所示。這是由羽流的卷吸現(xiàn)象導致的。當羽流初始流速較低時，瞬時圖像與時均邊界基本吻合；當羽流初始流速較大時，瞬時圖像與時均邊界存在明顯的突出和凹陷，而且瞬時圖像顆粒分布不均勻且不對稱，羽流中心處的高濃度區(qū)域在邊緣環(huán)境流體卷吸的作用下向一側(cè)偏移?？梢哉J為，初始流速高的羽流，卷吸或湍流擴散效應更加顯著。

圖18 瞬時圖像與時均邊界對比

針對工況1、3、5，將羽流中心z＝6D位置處的濃度作為標準質(zhì)量濃度c0，則羽流中心相對濃度c／c0隨高度z／D的變化規(guī)律如圖19 所示。隨著與羽流出口距離的增加，羽流中心相對濃度不斷降低，不同工況下變化趨勢基本相似。對于初始流速較大的羽流，卷吸效應更明顯、顆粒相擴散率更大，因此顆粒的擴散濃度下降速度相對較快。

圖19 不同初始流速下羽流中心相對濃度變化

羽流中顆粒在流體阻力的作用下動量逐漸減小，最終速度接近穩(wěn)定自由沉降速度［12］。選擇工況7、8、9，對低初始質(zhì)量濃度（c0＝5 g／L）、不同初始流速（u0＝0.99～2.45 m／s）的羽流進行試驗，分析距離羽流出口位置z／D＝17.5～22.5 范圍內(nèi)中心軸線顆粒平均速度um，測試分析區(qū)域為0.78 cm×0.20 cm的矩形區(qū)域，結(jié)果如圖20 所示。在相同高度處，初始流速較高的羽流顆粒速度衰減率相對更高。

圖20 不同初始流速下羽流中心顆粒相對速度變化

4 結(jié)語

提出了模擬多相羽流排放和擴散的試驗裝置與測試分析方法?；跁r均灰度圖像，采用邊緣識別算法捕捉羽流邊界和輪廓，獲取顆粒相的擴散率；針對顆粒濃度，基于時均灰度圖像并應用高斯分布假設(shè)，將二維平面圖像中灰度值所表征的視場平均濃度轉(zhuǎn)化為該位置羽流中心點濃度；針對顆粒速度，使用改進的PTV算法，提取非高濃度區(qū)域中顆粒速度，實現(xiàn)羽流邊緣等較低濃度區(qū)域內(nèi)顆粒速度和運動軌跡的連續(xù)捕捉。根據(jù)典型工況下試驗，驗證了該方法的準確性和可靠性。進一步研究了不同初始流速下羽流的擴散特征。結(jié)果表明，羽流初始流速越大，邊緣卷吸效應越顯著，顆粒擴散率越大，顆粒濃度與速度隨高度的衰減越快。