任意運(yùn)動(dòng)變形邊界流場(chǎng)測(cè)量PIV算法的研究：薄膜渦激振動(dòng)

曹洪才，施圣賢，劉應(yīng)征

（上海交通大學(xué) 動(dòng)力機(jī)械及工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200240）

0 引 言

隨著人們對(duì)仿生機(jī)械、微型飛行器等研究的日益深入，以前只關(guān)注機(jī)械結(jié)構(gòu)設(shè)計(jì)的單一研究思路被擯棄，研究人員開始結(jié)合流體力學(xué)、結(jié)構(gòu)力學(xué)等多學(xué)科進(jìn)行優(yōu)化設(shè)計(jì)，深入研究邊界變形與流體流固耦合的物理機(jī)制及其對(duì)氣動(dòng)性能的影響［1］。此外，流固耦合現(xiàn)象在自然界和工程實(shí)踐中很常見，如鳥類振翅、旗幟飄動(dòng)、心臟瓣膜開合過程、機(jī)翼顫振等［2-4］。在針對(duì)這類問題的數(shù)值建模和計(jì)算分析中，由于涉及流體和固體控制方程，兩種非線性控制方程存在強(qiáng)烈耦合，使得這類問題的理論分析和數(shù)值計(jì)算面臨很大的挑戰(zhàn)，所建立的數(shù)值模型及其計(jì)算結(jié)果也需要高精度的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行驗(yàn)證［5］。因此，精確測(cè)得強(qiáng)烈非線性流固耦合過程中的流場(chǎng)信息，尤其是運(yùn)動(dòng)變形邊界處的局部流動(dòng)特征，對(duì)于流固耦合機(jī)制的理論分析和數(shù)值計(jì)算都尤為重要。激光粒子圖像測(cè)速技術(shù)（PIV）由于具有全場(chǎng)、非接觸測(cè)量等特點(diǎn)，被廣泛應(yīng)用于流場(chǎng)實(shí)驗(yàn)研究中。然而，將傳統(tǒng)的PIV 測(cè)量技術(shù)用于存在任意運(yùn)動(dòng)變形邊界問題的流固耦合流場(chǎng)測(cè)試時(shí)，互相關(guān)計(jì)算窗口中可能包含流固邊界，會(huì)導(dǎo)致不正確甚至錯(cuò)誤的測(cè)量結(jié)果［6］。此外，傳統(tǒng)PIV 算法在用于流固耦合流場(chǎng)測(cè)試時(shí)，無法獲得流固耦合邊界的運(yùn)動(dòng)特征。因此，將PIV 應(yīng)用于流固耦合這種復(fù)雜流場(chǎng)測(cè)試時(shí)，需要充分考慮這些特點(diǎn)，發(fā)展能智能處理任意運(yùn)動(dòng)變形邊界條件下的PIV 流場(chǎng)測(cè)試技術(shù)［7］。對(duì)于強(qiáng)流固耦合問題，由于存在邊界時(shí)變、變形大、邊界運(yùn)動(dòng)復(fù)雜等特點(diǎn)，主要難點(diǎn)在于大變形運(yùn)動(dòng)邊界的智能識(shí)別、精確數(shù)字化，以及邊界附近局部流場(chǎng)信息的精確定量分析。針對(duì)這一問題，采用基于Radon變換的滑動(dòng)窗口邊界識(shí)別方法以精確獲取任意運(yùn)動(dòng)變形邊界信息，并應(yīng)用貼體自適應(yīng)互相關(guān)計(jì)算窗口測(cè)得瞬態(tài)流場(chǎng)信息，大大提高了流固耦合問題的PIV 測(cè)試能力。最后，將提出的新型PIV 算法應(yīng)用于柔性薄膜渦激共振流場(chǎng)的PIV 實(shí)驗(yàn)測(cè)試，獲得了薄膜在各種大變形狀態(tài)下的高分辨率瞬態(tài)流場(chǎng)數(shù)據(jù)。

1 任意運(yùn)動(dòng)變形邊界流場(chǎng)測(cè)量PIV算法

1.1 算法基本原理

傳統(tǒng)PIV 算法一般只適合于模型固定的流場(chǎng)測(cè)量，而流固耦合相關(guān)問題具有模型邊界時(shí)變、變形量大等特點(diǎn)，對(duì)傳統(tǒng)PIV 流場(chǎng)測(cè)量提出了挑戰(zhàn)。采用基于邊界信息的貼體自適應(yīng)PIV 互相關(guān)計(jì)算窗口可以獲得高精度近壁流場(chǎng)信息，但這一方法的前提是需要從“高噪聲”粒子圖像中提取運(yùn)動(dòng)變形邊界的精確位置信息，即任意運(yùn)動(dòng)變形邊界的識(shí)別問題。對(duì)于運(yùn)動(dòng)邊界的識(shí)別而言，這里所說的“噪聲”就是示蹤粒子本身。Radon變換是一種檢測(cè)直線的圖像處理算法，對(duì)灰度圖像I（x，y）進(jìn)行Radon變換，即在灰度圖像平面內(nèi)沿不同的直線做線積分，可以理解為灰度圖像在ρ-θ空間的投影。

直接使用Radon變換難以處理信噪比低的邊界信息，因而需要對(duì)算法進(jìn)行必要的修正和改良。提出了基于Radon變換的滑動(dòng)窗口邊界識(shí)別方法（如圖1所示）：（1）根據(jù)初態(tài)條件確定粒子圖像的初始窗口，進(jìn)行Radon變換，根據(jù)灰度值積分偵測(cè)初始窗口中的直線特征；（2）基于前一個(gè)窗口識(shí)別的直線特征，預(yù)測(cè)下一個(gè)檢測(cè)窗口的中心位置，并基于模板匹配的互相關(guān)方法對(duì)預(yù)測(cè)的窗口中心位置進(jìn)行修正；（3）基于先期位置信息，建立窗口中邊界位置的概率模型，進(jìn)行圖像灰度值的加權(quán)，以降低圖像中示蹤粒子的干擾；（4）經(jīng)過預(yù)測(cè)—修正—加權(quán)的窗口圖像，圖像信噪比得到極大提升。對(duì)灰度值進(jìn)行Radon變換，檢測(cè)低信噪比圖像中的直線特征。如此循環(huán)，可以識(shí)別出邊界分割開的若干小直線線段，進(jìn)行曲線擬合，可以得到粒子圖像中運(yùn)動(dòng)變形邊界的精確位置。

圖1 基于Radon變換邊界識(shí)別方法示意圖Fig.1 Interface tracking based on Radon transformation

為驗(yàn)證邊界識(shí)別算法的精確性，根據(jù)已知流場(chǎng)信息和邊界運(yùn)動(dòng)信息數(shù)字合成粒子圖像對(duì)。在數(shù)字合成粒子圖像中，邊界位置已預(yù)先知道，可和采用Radon變換算法識(shí)別出來的邊界信息進(jìn)行誤差對(duì)比，驗(yàn)證邊界識(shí)別的正確性及精度。圖2所示為數(shù)字合成粒子圖像，紅色圓圈為邊界識(shí)別程序識(shí)別出來的邊界信息。

圖2 數(shù)字合成粒子圖像及邊界識(shí)別示意圖Fig.2 Sketch of synthetic PIV image and interface tracking

為考察粒子濃度大小也即“噪聲”對(duì)邊界識(shí)別精度的影響，生成了不同粒子濃度的數(shù)字合成圖像來模擬不同信噪比環(huán)境下邊界識(shí)別的適用性。如圖3所示為32pixel×32pixel窗口中粒子數(shù)量分別為5 個(gè)和20個(gè)的數(shù)字合成圖像示意圖。

圖3 不同信噪比的數(shù)字合成粒子圖像Fig.3 Synthetic PIV images with different SNR

圖4 粒子濃度與信噪比、識(shí)別精確性的關(guān)系曲線Fig.4 SNR and tracking accuracy under different particle density

如圖4所示為粒子濃度與數(shù)字合成圖像信噪比的關(guān)系曲線，以及邊界信息識(shí)別誤差大小。隨著32pixel×32pixel像素窗口中粒子濃度從5個(gè)逐漸增加到20個(gè)，數(shù)字合成圖像的信噪比隨之降低。在圖像中加入白噪聲后，圖像的信噪比進(jìn)一步降低，如紅色曲線所示。對(duì)邊界信息識(shí)別的誤差絕對(duì)值進(jìn)行了統(tǒng)計(jì)，可見識(shí)別誤差在0.45個(gè)像素以內(nèi)。

圖5 傳統(tǒng)矩形和貼體自適應(yīng)PIV 互相關(guān)計(jì)算窗口示意圖Fig.5 Sketch of traditional and body-fitted window for cross-correlation

傳統(tǒng)PIV 算法采用規(guī)則正交的矩形互相關(guān)計(jì)算窗口，計(jì)算具有變形邊界的流固耦合流場(chǎng)時(shí)，會(huì)出現(xiàn)對(duì)邊界附近不同側(cè)粒子進(jìn)行互相關(guān)計(jì)算的問題，引起當(dāng)?shù)亓魉俚恼`判。如圖5所示，為提高具有運(yùn)動(dòng)邊界流場(chǎng)的測(cè)量精度，在識(shí)別出流固耦合運(yùn)動(dòng)邊界信息之后，基于運(yùn)動(dòng)變形邊界信息生成貼體自適應(yīng)互相關(guān)計(jì)算窗口，然后采用標(biāo)準(zhǔn)互相關(guān)算法計(jì)算當(dāng)?shù)厮俣仁噶?。采用此算法，可精確獲取流固耦合邊界附近流場(chǎng)信息。

1.2 基于數(shù)字合成粒子圖像的算法驗(yàn)證

為驗(yàn)證基于邊界信息的貼體自適應(yīng)窗口改進(jìn)PIV 算法的正確性，首先用計(jì)算流體手段計(jì)算了一個(gè)具有彎曲邊界的流場(chǎng)，其速度分布是已知的，稱為標(biāo)準(zhǔn)速度場(chǎng)。然后，采用數(shù)字合成PIV 粒子圖像對(duì)的方法，數(shù)字合成了兩幀粒子圖像。通過互相關(guān)運(yùn)算計(jì)算出兩幀粒子圖像的速度場(chǎng)，與標(biāo)準(zhǔn)速度場(chǎng)進(jìn)行對(duì)比，可驗(yàn)證改進(jìn)PIV 算法的正確性。

圖6 改進(jìn)PIV 算法的數(shù)字合成圖像驗(yàn)證Fig.6 Velocity accuracy for improved PIV algorithm by synthetic PIV images

如圖6（a）所示為利用改進(jìn)PIV 算法計(jì)算出的速度云圖和矢量圖，而圖6（b）是通過數(shù)值計(jì)算獲得的標(biāo)準(zhǔn)流場(chǎng)速度云圖和矢量圖。圖中的小方格是生成的貼體自適應(yīng)互相關(guān)計(jì)算窗口，通過互相關(guān)運(yùn)算在每個(gè)窗口生成一個(gè)速度矢量，如圖中黑色小箭頭所示。對(duì)比兩圖對(duì)應(yīng)窗口的速度矢量，可以發(fā)現(xiàn)改進(jìn)PIV算法計(jì)算的流場(chǎng)流態(tài)與標(biāo)準(zhǔn)流場(chǎng)基本一致。圖中顏色代表該互相關(guān)計(jì)算窗口計(jì)算出的速度值，在運(yùn)動(dòng)邊界附近，貼體自適應(yīng)窗口經(jīng)過互相關(guān)運(yùn)算的速度值相對(duì)誤差在5%以內(nèi)。總體來看，改進(jìn)PIV 算法計(jì)算出的速度云圖與標(biāo)準(zhǔn)的速度場(chǎng)非常接近，計(jì)算結(jié)果基本正確，尤其在大變形邊界附近沒有出現(xiàn)傳統(tǒng)PIV 算法的明顯問題。

2 實(shí)驗(yàn)應(yīng)用

2.1 水槽及實(shí)驗(yàn)?zāi)Ｐ?/h3>

為展示所發(fā)展的任意運(yùn)動(dòng)變形邊界條件下流場(chǎng)的PIV 實(shí)驗(yàn)測(cè)量效果，選擇布置在鈍體尾跡區(qū)域中的柔性薄膜渦激共振現(xiàn)象作為實(shí)驗(yàn)對(duì)象。流場(chǎng)測(cè)量實(shí)驗(yàn)在低速閉式循環(huán)水洞［8］中進(jìn)行，由進(jìn)口水箱、整流段、收縮段、實(shí)驗(yàn)段、循環(huán)水箱、磁力泵、變頻儀以及管路系統(tǒng)等組成。進(jìn)口水箱處設(shè)置有孔板，以削弱水泵出水的振蕩，達(dá)到初步整流的效果。水槽收縮段長(zhǎng)450mm，收縮型線采用三次方曲線。實(shí)驗(yàn)段橫截面為正方形，截面尺寸為100mm×100mm。為了減小水泵振動(dòng)對(duì)流場(chǎng)的影響，采用變頻儀控制的磁力泵驅(qū)動(dòng)，可以通過改變變頻儀的頻率來調(diào)節(jié)流速。水槽的具體設(shè)計(jì)及結(jié)構(gòu)詳見文獻(xiàn)［9］。

實(shí)驗(yàn)時(shí)，在封閉水槽實(shí)驗(yàn)段的中央放置一個(gè)剛性鈍體，鈍體采用直徑D=8mm 的方柱。鈍體下游一倍D處放置一端固定的柔性薄膜，柔性薄膜的上下端和另一端自由。柔性壓電薄膜的尺寸為：長(zhǎng)度L=41mm，寬度b=10mm，厚度h=0.1mm。自由來流速度u0=1.2m／s，雷諾數(shù)Re=u0D／μ=9930。柔性薄膜與鈍體尾流之間的流固耦合作用示意圖如圖7所示。流場(chǎng)上游經(jīng)過方柱繞流后，在下游所自然形成的旋渦脫落造成柔性薄膜變形運(yùn)動(dòng)，驅(qū)使柔性薄膜振動(dòng)。由于柔性薄膜的結(jié)構(gòu)振動(dòng)使流場(chǎng)發(fā)生改變，流場(chǎng)變化和薄膜變形運(yùn)動(dòng)在此過程中形成強(qiáng)耦合作用。

圖7 鈍體尾流柔性薄膜渦激振動(dòng)示意圖Fig.7 Sketch of vortex induced vibration of flexible film

PIV 實(shí)驗(yàn)時(shí)，示蹤粒子采用密度為1.03g／mm3的空心玻璃微珠，直徑為20～30μm。照明采用半導(dǎo)體連續(xù)式激光器，測(cè)量區(qū)域片光厚度小于1mm。為提高空間分辨率，實(shí)驗(yàn)中采用了1600萬像素的CCD相機(jī)。通常，相機(jī)的時(shí)間分辨率和空間分辨率是矛盾的?？臻g分辨率高的相機(jī)，其時(shí)間分辨率一般較低。采用了高空間分辨率的CCD相機(jī)，其采樣頻率較低為1f／s，采集窗口大小為4904pixel×3280pixel。測(cè)得薄膜各種運(yùn)動(dòng)變形狀態(tài)時(shí)的粒子圖像后，采用Radon變換對(duì)粒子圖像中的薄膜形態(tài)進(jìn)行圖像處理精確獲取薄膜形態(tài)曲線的數(shù)字化信息，然后生成貼體自適應(yīng)圖像互相關(guān)分析窗口，采用北京立方天地公司的Micro Vec粒子圖像分析軟件計(jì)算得到速度矢量。計(jì)算中采用圖像偏置及迭代算法［10］?；ハ嚓P(guān)計(jì)算的判讀窗口大小為64pixel×64pixel，步長(zhǎng)為32pixel?；ハ嚓P(guān)運(yùn)算的峰值采用高斯擬合方法確定［11］，計(jì)算精度可達(dá)到±0.1pixel。

2.2 實(shí)驗(yàn)結(jié)果及分析

圖8所示為Re=9930下鈍體尾渦柔性薄膜渦激振動(dòng)不同時(shí)刻的PIV 粒子圖像，圖中藍(lán)色所示是程序所識(shí)別的邊界位置信息。這些PIV 粒子圖可以使我們對(duì)柔性薄膜渦激振動(dòng)這種強(qiáng)烈的流固耦合現(xiàn)象有直觀的認(rèn)識(shí)。從圖中可以看出，方柱脫落的旋渦，由于柔性薄膜的阻隔作用，上下剪切層無法交換能量，旋渦貼附薄膜向下游輸運(yùn)。這樣，上下剪切層能量存在差異，柔性薄膜兩側(cè)受到壓力差而發(fā)生彎曲變形，在不同時(shí)刻呈現(xiàn)不同的振動(dòng)狀態(tài)。

圖8 柔性薄膜連續(xù)振動(dòng)的不同時(shí)刻粒子圖Fig.8 Particle images of voetex induced vibration in different frames

圖9 柔性薄膜渦激振動(dòng)不同時(shí)刻的速度云圖Fig.9 Contour plot of velocity field and streamline

圖9為柔性薄膜渦激振動(dòng)過程中4個(gè)不同時(shí)刻的流場(chǎng)速度云圖（Re=9930）。圖9（a）鈍體尾流中出現(xiàn)兩個(gè)大尺度旋渦，薄膜上方為順時(shí)針旋渦，下方為逆時(shí)針旋渦。薄膜上方的旋渦處于發(fā)展階段，距離方柱的位置相對(duì)較遠(yuǎn)。圖9（b）中，上方的旋渦向下游輸運(yùn)過程，在其影響下，柔性薄膜發(fā)生向下的彎曲變形。同時(shí)，在上游已經(jīng)生成新的旋渦。圖9（c）中，鈍體上、下緣脫落的旋渦貼附柔性薄膜輸運(yùn)，在柔性薄膜前半段出現(xiàn)兩個(gè)尺度相同的反向旋渦。在下游，柔性薄膜后半段被下方旋渦迅速向上揚(yáng)起。圖9（d）中，在柔性薄膜上方存在兩個(gè)旋渦，旋渦向下游輸運(yùn)過程中迫使柔性薄膜發(fā)生彎曲變形。而柔性薄膜下方只有一個(gè)逆時(shí)針旋渦，其能量剛好可以平衡上方的順時(shí)針旋渦。這與文獻(xiàn)［12］的計(jì)算流體動(dòng)力學(xué)數(shù)值模擬結(jié)論是一致的。

3 結(jié) 論

提出了一種用于測(cè)量任意運(yùn)動(dòng)變形邊界流場(chǎng)的PIV 計(jì)算方法，可以智能識(shí)別強(qiáng)烈流固耦合問題中的運(yùn)動(dòng)變形邊界，繼而生成基于邊界信息的貼體自適應(yīng)互相關(guān)計(jì)算窗口，獲得任意運(yùn)動(dòng)變形邊界附近的流場(chǎng)精確信息。采用這一算法，對(duì)已知運(yùn)動(dòng)變形邊界流場(chǎng)信息的數(shù)字合成圖像進(jìn)行系統(tǒng)的PIV 分析，發(fā)現(xiàn)所獲得的變形邊界附近的流場(chǎng)信息與已知流場(chǎng)信息高度一致，證明了所提出的PIV 算法的合理性。最后，利用PIV 測(cè)量了低速閉式循環(huán)水洞中鈍體尾渦柔性薄膜渦激振動(dòng)現(xiàn)象，采用改進(jìn)PIV 算法對(duì)4 種典型薄膜形態(tài)下的流場(chǎng)進(jìn)行了分析。實(shí)驗(yàn)測(cè)試結(jié)果表明，鈍體尾流周期脫落的旋渦與尾流中的柔性結(jié)構(gòu)形成強(qiáng)烈的非線性耦合作用，流態(tài)特征與薄膜形態(tài)變化的關(guān)聯(lián)特征與相關(guān)文獻(xiàn)結(jié)果一致。這說明所提出的任意運(yùn)動(dòng)邊界流場(chǎng)測(cè)量PIV 算法可用于強(qiáng)流固耦合問題的流場(chǎng)實(shí)驗(yàn)測(cè)量，為流固耦合物理機(jī)制理論分析和數(shù)值計(jì)算方法的驗(yàn)證提供了有力的實(shí)驗(yàn)方法。

參考文獻(xiàn)：interaction［J］.Annual Review of Fluid Mechanics，2001，33（1）：445-490.

［3］ SARPKAYA T.A critical review of the intrinsic nature of vortex-induced vibrations［J］.Journal of Fluids and Structures，2004，19（4）：389-447.

［4］ GABBAI R D，BENAROYA H.An overview of modelingand experiments of vortex-induced vibration of circular cylinders［J］.Journal of Sound and Vibration，2005，282：575-616.

［5］ 王文全，張立翔，閆妍，等.方柱繞流誘發(fā)的彈性薄板流固耦合特性研究［J］.工程力學(xué)，2011，28（3）：17-22.

［6］ ADRIAN R J.Particle image techniques for experimental fluid mechanics［J］.Annual Review of Fluid Mechanics，1991，23：261-304.

［7］ JERRY Westerweel.Particle image velocimetry for complex and turbulent flows［J］.Annual Review of Fluid Mechanics，2013，45：409-436.

［8］ 陳建民.鈍體尾流壓電薄膜渦激共振流動(dòng)特性的實(shí)驗(yàn)研究［D］.［碩士學(xué)位論文］.上海：上海交通大學(xué)，2013.

［9］ SHI Sheng-xian，LIU Ying-zheng.Flapping dynamics of a low aspect-ratio energy-harvesting membrane immersed in a square cylinder wake［J］.Experimental Thermal and Fluid Science，2013，46：151-161.

［10］WESTERWEEL J，DABIRI D，GHARIB M.The effect of a discrete window offset on the accuracy of crosscorrelation analysis of digital PIV recordings［J］.Experiments in Fluids，1997，23（1）：20-28.

［11］YASUHIKO S，NISHIOS，OKUNO T.A highly accurate iterative PIV technique using a gradient method［J］.Measurement Science and Technology，2000，11（12）：1666-1673.

［12］LIEW K M.A computational approach for predicting the hydroelasticity of flexible structures based on the pressure Poisson equation［J］.International Journal for Numerical Methods in Engineering，2007，72：1560-1583.