曹洪才，施圣賢，劉應征

（上海交通大學 動力機械及工程教育部重點實驗室，上海 200240）

0 引 言

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

1 任意運動變形邊界流場測量PIV算法

1.1 算法基本原理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 實驗應用

2.1 水槽及實驗模型

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

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

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

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

2.2 實驗結(jié)果及分析

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

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

圖9 柔性薄膜渦激振動不同時刻的速度云圖Fig.9 Contour plot of velocity field and streamline

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

3 結(jié) 論

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

參考文獻：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］.工程力學，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］.［碩士學位論文］.上海：上海交通大學，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.