陳金鵬 李貴山 李雙科 馬宏鋒

(蘭州工業(yè)學(xué)院電氣工程系，甘肅 蘭州 730050)

0 引言

不管是普通的流體視覺化，還是由計(jì)算機(jī)輔助的流體視覺化，它們?cè)诹黧w動(dòng)力學(xué)中對(duì)于觀察和監(jiān)控流體結(jié)構(gòu)都是非常有用的工具。流體視覺化通過在流體中引入染色物或煙氣形式的污染物(即被動(dòng)示蹤體)，使流體的結(jié)構(gòu)變得可視。

根據(jù)不同流體的密度，只有選擇合適的污染物，才能確保很好地跟蹤流體的流動(dòng)，視覺化才能集中于被動(dòng)示蹤體以及它在時(shí)間和空間上的濃度的變化，從而識(shí)別和確認(rèn)特定流體的特性。

對(duì)于用數(shù)碼相機(jī)采集的流體/被動(dòng)示蹤體流動(dòng)的圖像，可以進(jìn)一步對(duì)其進(jìn)行處理，以用于不同目的的研究，如渦流分離分析［1］、回流區(qū)位置確認(rèn)［2］、監(jiān)控氣蝕結(jié)構(gòu)［3］、流體流動(dòng)對(duì)沉積的影響［4］等。在這個(gè)階段出現(xiàn)的問題是能否借助計(jì)算機(jī)視覺獲取流體的定量特性，也就是說，能否采用被動(dòng)示蹤體將速度場(chǎng)這樣的流體運(yùn)動(dòng)量從一組流體圖像中抽取出來。

目前，流體速度場(chǎng)測(cè)量主要采用粒子圖像速度測(cè)量?jī)x。該測(cè)量?jī)x通過在流體中引入固體微粒，采用機(jī)器視覺的方法進(jìn)行測(cè)量。借助機(jī)器視覺對(duì)流體進(jìn)行定性和定量的研究，至少應(yīng)同時(shí)使用兩種方法或系統(tǒng)［5－6］，但這樣的組合系統(tǒng)往往很復(fù)雜。此外，大部分普通的測(cè)量方法都不能測(cè)量非平穩(wěn)或周期性流體結(jié)構(gòu)的流體速度。

本文僅通過采用視覺方法對(duì)加入示蹤體的流體的速度場(chǎng)進(jìn)行測(cè)量，并在實(shí)際使用中檢驗(yàn)這種方法的可行性。該方法是利用一種已知的物理關(guān)系，把污染物濃度和流體速度結(jié)合起來的方法。

1 研究方法

在所描述的研究方法中，最重要的步驟是把速度向量場(chǎng)與濃度的標(biāo)量場(chǎng)聯(lián)系起來，后者用灰度圖像表示。通過濃度N(即單位體積流體中污染體的分子數(shù))表征污染體引入流體的情況。如果采用如下速度供給污染體分子，則可以使x方向上N的梯度保持不變。

式中:Φ為X方向的擴(kuò)散流量;D為分子擴(kuò)散系數(shù)。

式(1)為菲克擴(kuò)散定律。如果考慮非穩(wěn)態(tài)擴(kuò)散過程，基本體積中濃度N的變化率可以用流出和流入該體積流體邊界的污染體的微分表示，則式(1)可以改寫為［7］:

當(dāng)有流體運(yùn)動(dòng)(即對(duì)流)時(shí)，用總的導(dǎo)數(shù)代替時(shí)間偏導(dǎo)數(shù)，則式(2)可進(jìn)一步展開為:

式(3)即為著名的對(duì)流擴(kuò)散方程，它表示示蹤體/污染體與流體動(dòng)力學(xué)的基本關(guān)系。其中，等式左邊第二項(xiàng)包括流體對(duì)流的影響。這樣，如果濃度N的值(即N的導(dǎo)數(shù)已知)和流體中污染體的擴(kuò)散系數(shù)D已知，則唯一的未知量就是速度v，即只有xi方向上的分量vi為未知量。

分子擴(kuò)散系數(shù)D是物質(zhì)的屬性。在一定壓力和溫度下，D與另一種物質(zhì)有關(guān)，污染體濃度N在不斷變化，因此，應(yīng)該即時(shí)測(cè)量。在這種情況下，應(yīng)該將它的值從流體灰度圖像中抽取出來。

引入跟蹤物的流體灰度圖像可形成具有不同灰度的像素的矩陣。一般來說，灰度的整數(shù)值從0(黑)到255(白)，它可以在0～1區(qū)間被歸一化。如果引入流體的污染物是用外光源照明的，那么它的灰度值一般比流體的灰度值高。流體中污染體的濃度越高，在灰度圖像中污染體的灰度值越高。因此，有:

式中:A為灰度圖像中選定像素窗口的平均灰度值。文獻(xiàn)［8］和文獻(xiàn)［9］對(duì)此已有論述，式(4)中灰度值和污染體濃度的比值是線性的。選定的像素窗口的平均灰度可以用如下關(guān)系［2］得到:

式中:E(i，j)為l×m像素窗口中實(shí)測(cè)的(i，j)點(diǎn)單個(gè)像素的灰度強(qiáng)度。

濃度的時(shí)間導(dǎo)數(shù)?N/?t可以利用如圖1所示的連續(xù)圖像中固定窗口的平均灰度計(jì)算得到。

圖1 流體連續(xù)圖像窗口Fig.1 Continous image windows of the fluid

利用已知的兩連續(xù)圖像的時(shí)間間隔Δt(通過相機(jī)的圖像采集頻率得到)，則?N/?t項(xiàng)可以用式(6)估算:

式中:ΔA為固定窗口中相距Δt的兩連續(xù)圖像之間平均灰度的差值。污染體濃度在空間上的偏導(dǎo)數(shù)?N/?xi可以通過計(jì)算流體的單個(gè)圖像獲得。為了得到平均灰度A在空間上的差值，將現(xiàn)在的窗口向選定的xi方向移動(dòng)，其示意圖如圖2所示。

圖2 在窗口中移動(dòng)單個(gè)圖像示意圖Fig.2 Diagram of moving single image in the window

空間導(dǎo)數(shù)可以用標(biāo)準(zhǔn)數(shù)字技術(shù)完成。式(7)為忽略截?cái)嗾`差的中心差分法［10］:

式(3)在兩維系統(tǒng)中表示兩個(gè)線性微分方程，它有兩個(gè)未知數(shù)，即兩個(gè)速度分量。求解這樣一組方程需要明確的初始條件，非常繁瑣和困難，因?yàn)榱黧w的運(yùn)動(dòng)特性一般都未知。但是，在速度分量和其空間導(dǎo)數(shù)仍未知的情況下，這組微分方程可以轉(zhuǎn)換為線性方程組。這組方程即使是超定的(即方程的個(gè)數(shù)多于未知數(shù)的個(gè)數(shù))，也可以通過用數(shù)字的方法很容易地解出。與微分方程組相比，線性方程組的未知數(shù)的個(gè)數(shù)增加了一倍，這就意味著對(duì)兩維系統(tǒng)至少需要四個(gè)一般方程。每一個(gè)方程需要兩個(gè)流體連續(xù)的圖像(因?yàn)橐?jì)算時(shí)間差)，即至少需要五個(gè)連續(xù)的流體圖像才能得到流體的兩個(gè)速度分量的信息。

上述研究方法是基于已知的物理關(guān)系，反映了借助機(jī)器視覺確定流體動(dòng)力學(xué)特性的基本原理。它的有效性可以用傳統(tǒng)的測(cè)量方法驗(yàn)證。為了達(dá)到這個(gè)目的，采用熱絲風(fēng)速計(jì)進(jìn)行試驗(yàn)驗(yàn)證。

2 試驗(yàn)結(jié)構(gòu)

測(cè)量來自垂直管空氣射流速度場(chǎng)的測(cè)量裝置如圖3所示。

圖3 測(cè)試部分的試驗(yàn)裝置Fig.3 Experimental device of the test section

流體的雷諾數(shù)在管道的出口處進(jìn)行測(cè)定，記為Re=1300。由于管道中的流體是層流，可以采用熱線風(fēng)速計(jì)法(hot-wire anemometer，HWA)和計(jì)算機(jī)輔助視覺法(computer-aided visualization，CAV)完成測(cè)量。在選定的測(cè)量點(diǎn)，用兩分支熱絲風(fēng)速計(jì)測(cè)定速度的分布。兩分支熱絲風(fēng)速計(jì)測(cè)量流向(y)和側(cè)向(x)瞬時(shí)速度波的步驟與 Bruun［11］和 Jφrgenson［12］的研究步驟一致。在9×17網(wǎng)孔中選定153個(gè)測(cè)量點(diǎn)，每?jī)蓚€(gè)點(diǎn)之間的距離用定位表設(shè)置為5 mm。

為了用熱絲風(fēng)速計(jì)測(cè)量流體速度，本文采用DANTEC 55P62恒溫?zé)峤z風(fēng)速計(jì)傳感器系統(tǒng)。該傳感器用鍍鉑鎢絲制成，長(zhǎng)度和直徑分別為1.25 mm、φ5 μm。風(fēng)速計(jì)放大器的截止頻率設(shè)置為10 kHz，風(fēng)速計(jì)的工作溫度為250℃。風(fēng)道中的兩分支熱絲風(fēng)速計(jì)傳感器的位置要求能從兩個(gè)方向測(cè)量瞬時(shí)速度。因此，風(fēng)速計(jì)傳感器采用PC控制的精確定位裝置進(jìn)行定位。熱絲風(fēng)速計(jì)的信號(hào)采用NI公司的16位數(shù)據(jù)采集板獲得。數(shù)據(jù)采集時(shí)間為6 s，采集頻率為50 kHz。

風(fēng)速計(jì)的輸出在傳送到具有SCXI模塊的A/D轉(zhuǎn)換器之前，先用頻率為10 kHz的低通四階貝塞爾濾波器予以濾波。系統(tǒng)采用LabVIEW軟件采集并存儲(chǔ)數(shù)據(jù)。熱絲風(fēng)速計(jì)的校準(zhǔn)是在專門校準(zhǔn)風(fēng)速計(jì)的測(cè)量臺(tái)上完成的。為了與國王定律［11］保持一致，對(duì)第一和第二個(gè)熱絲的常數(shù)設(shè)定如下:

此外，在校準(zhǔn)和測(cè)量過程中還對(duì)溫度進(jìn)行了補(bǔ)償。

溫度的測(cè)量采用A級(jí)四線電阻溫度計(jì)Pt100和Agilent 34970A數(shù)據(jù)采集器。根據(jù)國王定律方程式［11］和Jφrgenson方程［12］，從熱絲風(fēng)速計(jì)的輸出計(jì)算得到實(shí)際的流向流速vy和側(cè)向流速vx。測(cè)量也是用同樣的設(shè)備完成的，且在測(cè)量之前進(jìn)行了校準(zhǔn)［13］。測(cè)量的總誤差［14］是由工作點(diǎn)的選擇和波動(dòng)的不確定性、定位起始點(diǎn)的選擇、風(fēng)扇的轉(zhuǎn)速、定位表對(duì)熱絲溫度傳感器的定位、線性化、A/D的分辨率、數(shù)據(jù)采集、溫度補(bǔ)償、濕度以及熱絲風(fēng)速計(jì)的頻率響應(yīng)有限等因素引起的。風(fēng)扇由一臺(tái)分離變壓器供電，便于與測(cè)量設(shè)備的供電電源分開。瞬時(shí)速度的整體測(cè)量誤差估計(jì)為測(cè)量值的2.8%。

在CAV測(cè)量中，為了使流體可視化，在充分發(fā)展流中引入含有氣化石蠟油的被動(dòng)示蹤煙氣［15］。測(cè)量點(diǎn)與上述熱絲風(fēng)速計(jì)的測(cè)量點(diǎn)完全相同。通過流體入口將被動(dòng)示蹤體引入管道。相機(jī)放置于與風(fēng)道的測(cè)量段垂直并且距離測(cè)量段表面1 m的位置。相機(jī)的圖像采集頻率為400 Hz，拍攝圖像的快門速度為0.6 ms，所記錄的圖像具有8位灰度值、600×250像素的分辨率。本文采用具有直線導(dǎo)光的不間斷無閃爍光源(Vega Velum DC 150 W)進(jìn)行照明;利用LabVIEW軟件包實(shí)現(xiàn)對(duì)相機(jī)的設(shè)置、圖像的采集和保存。在每一個(gè)測(cè)量點(diǎn)采集1000個(gè)連續(xù)的圖像，圖像的個(gè)數(shù)受被動(dòng)示蹤體煙氣的發(fā)生時(shí)間限制。被動(dòng)示蹤體瞬時(shí)濃度的總體測(cè)量誤差估計(jì)為測(cè)量值的 3.8%［15］。

3 試驗(yàn)結(jié)果

試驗(yàn)中采集的圖像序列用Matlab程序包進(jìn)行處理。兩個(gè)速度分量的測(cè)定與試驗(yàn)測(cè)量(熱絲法)的情況一樣，都是在同樣的測(cè)量點(diǎn)、在x-y平面利用機(jī)器視覺測(cè)定的。根據(jù)文獻(xiàn)［16］，空氣(T=23℃、p=100 kPa)中的癸烷(C10H22)蒸汽 D=7.01 ×10－12m2/s。用 30 個(gè)連續(xù)圖像序列來確定x-y平面上的流速場(chǎng)，序列中圖像的數(shù)量保證了方程式(3)即使在被動(dòng)示蹤體稀少(測(cè)量平面x-y外沿)的情況下，對(duì)所有的測(cè)量點(diǎn)都有效。另一方面，當(dāng)圖像序列很大時(shí)，求取速度分量的平均值需要非常多的計(jì)算時(shí)間，有可能導(dǎo)致給定采集頻率上的信息丟失。為了減小導(dǎo)數(shù)數(shù)值計(jì)算的截?cái)嗾`差，式(7)中步長(zhǎng)Δxi應(yīng)選得盡可能小。因此，步長(zhǎng)接近于1個(gè)像素，它的大小取決于圖像分辨率。

本文采用HWA和CAV兩種方法測(cè)量流速場(chǎng)。測(cè)量結(jié)果顯示，在大部分測(cè)量點(diǎn)，HWA法和CAV法的測(cè)量結(jié)果一致。

兩種方法的測(cè)量結(jié)果表明，離開節(jié)流孔的射流流速減小，射流外邊沿的流體沿徑向朝外側(cè)運(yùn)動(dòng)(射流沿軸向和徑向散開)。盡管如此，兩種方法所測(cè)得的速度矢量的方向和大小都有差別，尤其在射流的外邊沿。為了正常工作，序列中的每個(gè)圖像都要有足夠的被動(dòng)示蹤體濃度值，而在這些位置，由于石蠟油蒸汽很少，影響了CAV法的效果。

對(duì)兩種方法測(cè)得的速度分量分別從流向方向(y)和正常方向(x)進(jìn)行定量分析。用x方向和y方向的歸一化值來表示結(jié)果，這里 x*、y*∈［0，1］。x*=0 和 x*=1分別對(duì)應(yīng)于圖3最左邊和最右邊的測(cè)量點(diǎn);類似地，y*=0和y*=1分別對(duì)應(yīng)于距離管道節(jié)流孔最遠(yuǎn)和最近的測(cè)量點(diǎn)。

經(jīng)歸一化流向，當(dāng)觀測(cè)點(diǎn)y*的值分別為0.125、0.375、0.625 和 0.875 時(shí)，兩種測(cè)量方法所測(cè)得的流速分量vy的值如圖4所示。

圖4 不同流向觀測(cè)點(diǎn)流速分量vyFig.4 Velocity component vyin different directional observing points

通過比較兩種測(cè)量結(jié)果可知，不管采用何種測(cè)量方法，y方向的流速具有相似的外形。但是，采用CAV法時(shí)vy的波動(dòng)比HWA法的測(cè)量結(jié)果要大得多，在遠(yuǎn)離管道節(jié)流閥處表現(xiàn)得更為突出。這種現(xiàn)象很可能是流體的漩渦和特殊測(cè)量方法的數(shù)據(jù)采集頻率共同作用的結(jié)果。當(dāng)選定雷諾數(shù)，且觀測(cè)點(diǎn)y*＞0.3時(shí)，管道中的流體形成了強(qiáng)烈的漩渦。

HWA法的數(shù)據(jù)采集頻率大約為50 kHz，在每個(gè)點(diǎn)的測(cè)量歷時(shí)約2 s。因此，HWA法獲得的結(jié)果是這段時(shí)間的平均值，這個(gè)時(shí)間足夠長(zhǎng)，足以保證大漩渦的旋轉(zhuǎn)對(duì)流體速度的平均值影響不大。對(duì)于同樣的時(shí)間段，一方面CAV法大約需要700個(gè)圖像序列，這將極大地影響處理/計(jì)算時(shí)間;另一方面，這種情況下CAV法的采集頻率僅為HWA法采集頻率的1/140。為了獲得更精確的結(jié)果，應(yīng)該增加圖像序列的數(shù)量，但該方法只適用于流體平均速度(波動(dòng))偏差較大、持續(xù)時(shí)間足夠長(zhǎng)且被CAV法檢測(cè)到的情況。即與HWA法相比，CAV法需要相機(jī)具有較高的采集頻率，兩者才有更好的可比性。

根據(jù)HWA法，圖4中每一個(gè)流向觀測(cè)點(diǎn)的vy(沿側(cè)向距離x的計(jì)算值)平均值的相對(duì)差如表1所示。

表1 vy平均值的相對(duì)差Tab.1 Relative difference of vyaverage values

表1表明四個(gè)選定的流向觀測(cè)點(diǎn)的平均值相對(duì)差大約為10%。因此，可以得出結(jié)論:與HWA法相比，即使在圖像序列很短的情況下，CAV法也能很好地預(yù)測(cè)vy的值。

根據(jù)CAV法和HWA法，vx平均值的相對(duì)差比較如表2所示。

表2 vx平均值的相對(duì)差Tab.2 Relative difference of vxaverage values

表2中，vx平均值相對(duì)差比表1中vy平均值相對(duì)差稍大一些。其主要原因是:與vy的值相比，vx值更小(接近于0，最多到0.04 m/s)。由于相對(duì)差再大也不會(huì)超過15%，所以，對(duì)基本CAV法而言，這可以認(rèn)為是較好的結(jié)果。

4 結(jié)束語

本文闡述了借助于機(jī)器視覺對(duì)流體動(dòng)力學(xué)性能量化分析的基本原則。通過流體中污染體(被動(dòng)示蹤體)的連續(xù)圖像序列，測(cè)量了流體的兩維速度場(chǎng)。這是一種通過對(duì)流擴(kuò)散方程把速度場(chǎng)和濃度場(chǎng)聯(lián)系起來的方法。所需參數(shù)的計(jì)算是用圖像處理和普通的數(shù)字技術(shù)來實(shí)現(xiàn)的。該方法本身是用傳統(tǒng)的熱絲風(fēng)速計(jì)法與機(jī)器視覺法同時(shí)進(jìn)行而予以驗(yàn)證的。

兩種測(cè)試方法的結(jié)果對(duì)比表明:視覺法(即使最基本的形式)能夠相當(dāng)精確地計(jì)算出流體的大小和方向;兩種方法在不同流向觀測(cè)點(diǎn)的速度分量外形是相似的，兩種結(jié)果的差別非常小。由此可以得出利用機(jī)器視覺測(cè)量流體速度場(chǎng)是可行的結(jié)論。

從現(xiàn)階段來看，視覺法的精度劣于傳統(tǒng)方法如HWA、PIV和LDA。因此，未來的研究不應(yīng)該僅僅專注于數(shù)字計(jì)算領(lǐng)域(如步長(zhǎng)大小、不同導(dǎo)數(shù)計(jì)算方法的應(yīng)用對(duì)最終結(jié)果的影響等)，而是更應(yīng)該專注于測(cè)量準(zhǔn)備工作領(lǐng)域(如照明、三維流體結(jié)構(gòu)、非定常性、序列中圖像的個(gè)數(shù)等對(duì)結(jié)果誤差的影響)。根據(jù)現(xiàn)有的算法和設(shè)備要求(合適的圖像分辨率和采集頻率)，視覺方法僅可應(yīng)用于有很高雷諾數(shù)的湍流，具有一定的有效性和局限性。

［1］Sirok B，Potocar E，Novak M.Analysis of the flow kinematics behind a pulsating adaptive airfoil using computer-aided visualization［J］.Strojniski Vestnik-Journal of Mechanical Engineering，2000，46(6):330－341.

［2］Sirok B，Bajcar T，Dular M.Reverse flow phenomenon in a rotating diffuser［J］.Journal of Flow Visualization and Image Processing，2002，9(2＆3):193 －210.

［3］Sirok B，Blagojevic B，Bajcar T，et al.Simultaneous study of pressure pulsation and structural fluctuations of a cavitated vortex core in the draft tube of a Francis turbine［J］.Journal of Hydraulic Research，2003，41(4):541 －548.

［4］Bajcar T，Blagojevic B，Sirok B，et al.Influence of flow properties on a structure of a mineral wool primary layer［J］.Experimental Thermal and Fluid Science，2007，32(2):440 －449.

［5］Hu H，Saga T，Kobayashi T，et al.Analysis of a turbulent jet mixing flow by using a PIV-PLIF combined system［J］.Journal of Visualization，2004，7(1):33 －42.

［6］Reungoat D，Riviere N，F(xiàn)aure J P.3C PIV and PLIF measurement in turbulent mixing［J］.Journal of Visualization，2007，10(1):99 －110.

［7］McComb W D.The physics of fluid turbulence［M］.Oxford:Clarendon Press，1996.

［8］Aider J L，Westfreid J E.Visualizations and PDF of the fluctuations of a passive scalar in a turbulent Goertler flow［J］.Experimental and Numerical Flow Visualization，1995，218:123 －130.

［9］Simoens S，Ayrault M.Concentration flux measurements of a scalar quantity in turbulent flows［J］.Experiments in Fluids，1994，16(3 －4):273－281.

［10］Mathews J H，Kurtis D F.Numerical methods using Matlab［M］.4th Edition.London:Prentice-Hall，2004.

［11］Bruun H H.Hot-wire anemometry-principles and signal analysis［M］.New York:Oxford University Press，1995.

［12］Dantec Dynamics A/S.How to measure turbulence with hot-wire anemometers-a practical guide［EB/OL］.［2002 －02 －01］［2012 －03 － 02］.Http://www.dantecdynamics.com/Admin/Public/Download.aspx?file=files%2fsupport_and_download%2fresearch_and_education%2fpracticalguide.pdf.

［13］Coleman H W，Steele W G.Experimentation and Uncertainty Analysis for Engineers［M］.New York:John Wiley ＆ Sons Inc.，1999.

［14］Eberlinc M，Sirok B，Hocevar M，et al.Numerical and experimental investigation of axial fan with trailing edge self-induced blowing［J］.Journal of Engineering Research，2009，73(3):129 －138.

［15］Eberlinc M，Dular M，Sirok B，et al.Influence of blade deformation on integral characteristic of axial flow fan［J］.Journal of Mechnical Engineering，2008，54(3):159 －169.

［16］Won D，Shaw C Y.Determining coefficients for mass-transfer models for volatile organic compound emissions from architectural coatings［C］//Proceedings of the 7th International Conference of Healthy Buildings，2003:1－6.