楔形體入水的時(shí)間解析-粒子圖像測(cè)速測(cè)試及數(shù)值研究

佘文軒，郭春雨，周廣利，吳鐵成，徐 鵬

(哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001)

物體入水問(wèn)題在船舶與海洋工程、航空航天等領(lǐng)域具有廣泛的工程應(yīng)用背景和科學(xué)研究?jī)r(jià)值，由于其涉及氣、液、固三者之間的耦合作用，近自由液面處的強(qiáng)非線性及強(qiáng)非定常性流動(dòng)給研究者帶來(lái)較大的技術(shù)挑戰(zhàn)[1].自Wagner[2]進(jìn)行入水問(wèn)題的開(kāi)創(chuàng)性研究后，出現(xiàn)了大量具有針對(duì)性和專(zhuān)業(yè)性的入水問(wèn)題研究手段與分析結(jié)果[3-5].

船舶航行過(guò)程中，由于較大的縱搖和垂蕩運(yùn)動(dòng)會(huì)導(dǎo)致危害性較大的底部砰擊作用[6]，大多數(shù)學(xué)者將船體橫剖面簡(jiǎn)化為楔形體模型進(jìn)行深入分析[7-9].Wu等[4]應(yīng)用邊界元法進(jìn)行二維楔形體自由入水的數(shù)值模擬，并輔以試驗(yàn)相驗(yàn)證.Wang等[7]分析了楔形體入水過(guò)程中砰擊階段、過(guò)渡階段以及空穴閉合階段的流體力學(xué)特性.隨著計(jì)算機(jī)技術(shù)迅速發(fā)展，多種數(shù)值算法被提出并應(yīng)用于入水問(wèn)題研究中，如光滑粒子(SPH)法[10]、流體體積(VOF)法[11].但在實(shí)驗(yàn)研究方面，大量學(xué)者仍采用加速度、位移、拉力、壓力傳感器等傳統(tǒng)實(shí)驗(yàn)技術(shù)手段[4,9]，該方法存在安裝繁瑣、測(cè)試物理量單一、單點(diǎn)接觸式測(cè)控等不足的問(wèn)題.20世紀(jì)80年代發(fā)展起來(lái)的粒子圖像測(cè)試(PIV)技術(shù)能夠無(wú)干擾地測(cè)量全局流場(chǎng)的瞬態(tài)速度信息，不僅為數(shù)值方法提供了更為全面的數(shù)據(jù)驗(yàn)證，還為流場(chǎng)壓力分布、渦量與湍流度等關(guān)鍵流動(dòng)參數(shù)的間接測(cè)量提供了可能[12].張志榮等[13-14]對(duì)楔形體入水時(shí)的砰擊瞬態(tài)流場(chǎng)進(jìn)行PIV測(cè)試.對(duì)于入水問(wèn)題，研究者往往更關(guān)注入水過(guò)程中的砰擊壓力，Oudheusden等[15-16]總結(jié)了多種基于PIV的壓力重構(gòu)方法，并指出壓力重構(gòu)結(jié)果的準(zhǔn)確性取決于特定的流動(dòng)問(wèn)題，目前尚無(wú)最優(yōu)的解決方案.特別是對(duì)于具有自由液面遷移和物體邊界運(yùn)動(dòng)的入水問(wèn)題，一些方法的適用性與精確性有待證實(shí).Nila[17]、Panciroli[18]和Jalalisendi[19]等采用PIV技術(shù)對(duì)楔形體入水的流場(chǎng)進(jìn)行測(cè)試和砰擊壓力的重構(gòu)，表明了基于PIV技術(shù)進(jìn)行楔形體入水砰擊壓力間接評(píng)估的可行性，但上述研究?jī)H采用單一的楔形體模型進(jìn)行砰擊壓力重構(gòu)方法的驗(yàn)證，其適用性有待進(jìn)一步探究.

本文應(yīng)用高頻響的時(shí)間解析PIV(TR-PIV)技術(shù)對(duì)不同底升角θ楔形體入水過(guò)程中的流場(chǎng)進(jìn)行測(cè)試，并依據(jù)TR-PIV瞬態(tài)測(cè)試結(jié)果進(jìn)行入水過(guò)程中瞬時(shí)流場(chǎng)壓力分布以及砰擊載荷重構(gòu)，同時(shí)進(jìn)行相應(yīng)試驗(yàn)工況下的數(shù)值模擬，對(duì)比分析基于TR-PIV數(shù)據(jù)進(jìn)行砰擊壓力重構(gòu)的準(zhǔn)確性和適用性，探討楔形體入水過(guò)程中的細(xì)節(jié)流動(dòng)特征.

1 試驗(yàn)裝置和TR-PIV系統(tǒng)

1.1 試驗(yàn)裝置

試驗(yàn)裝置如圖1所示，鋁型材框架、直線滑軌以及亞克力材質(zhì)的透明水箱構(gòu)成試驗(yàn)裝置的主要部分，水箱長(zhǎng)800 mm，寬500 mm，高500 mm，由鋁型材框架支撐，長(zhǎng)800 mm的直線滑軌固定在水箱的正上方.滑軌中間設(shè)置有電磁裝置，楔形體上端設(shè)有連接滑塊，楔形體與連接滑塊組成1個(gè)整體，受電磁裝置的控制，能夠在滑軌間自由滑動(dòng)，并以垂直的速度入水.

圖2所示為4個(gè)底升角分別為20°、25°、30° 及35° 的楔形體模型，由剛性聚乳酸(PLA)材料經(jīng)3D打印制成，寬度為200 mm，用配重調(diào)整模型質(zhì)量，使4個(gè)楔形體模型與連接滑塊的總重均為0.6 kg，從距離水面30 cm處自由落下.忽略空氣阻力的影響，根據(jù)理論得到其對(duì)應(yīng)的入水初速度為2.425 m/s，激光片光照射楔形體的中間位置切面，可認(rèn)為中間切面處的流動(dòng)近似為二維流動(dòng)[18-19].試驗(yàn)時(shí)水溫為20 ℃，水的密度為998.16 kg/m3，重力加速度g為9.8 m/s2.

1.2 TR-PIV系統(tǒng)

本文中所采用的TR-PIV系統(tǒng)主要由1臺(tái)連續(xù)激光器、NAC Memrecam HX-6高速CMOS相機(jī)以及計(jì)算機(jī)組成.10 W釔鋁石榴石晶體連續(xù)激光器(波長(zhǎng)532 nm)作為試驗(yàn)光源，測(cè)試區(qū)域的片光厚度約為1 mm.相機(jī)的內(nèi)存為8 GB，試驗(yàn)時(shí)相機(jī)的空間分辨率為 1 280 像素×1 000 像素，采集速率為 5 000 Hz，圖像深度為16 bit，楔形體底端初始接觸水面的時(shí)刻設(shè)置為0時(shí)刻.

楔形體入水過(guò)程具有流動(dòng)對(duì)稱(chēng)性，進(jìn)行PIV測(cè)試時(shí)，只觀測(cè)入水過(guò)程中的一半流場(chǎng)區(qū)域.采用密度約為1.03 g/mm3、粒徑約20 μm的聚酰胺微珠PSP-20作為實(shí)驗(yàn)示蹤粒子.網(wǎng)格大小為10 mm×10 mm 的標(biāo)定板用于PIV系統(tǒng)標(biāo)定.PIV系統(tǒng)測(cè)量區(qū)域的大小約為165 mm×130 mm，對(duì)應(yīng)于7.63 像素/mm.記錄每次入水砰擊過(guò)程前20 ms的流場(chǎng)信息，即采集101張PIV粒子圖像.粒子圖像的互相關(guān)計(jì)算采用基于MATLAB平臺(tái)的開(kāi)源軟件PIVlab，使用基于快速傅里葉變換(FFT)的多重網(wǎng)格迭代技術(shù)，設(shè)置有64像素×64像素、32像素×32像素、16像素×16像素多重判讀窗口，相鄰窗口重疊率為50%，采用3點(diǎn)高斯亞像素插值進(jìn)行互相關(guān)峰值擬合，擬合結(jié)果精度約為0.1像素[20]，流場(chǎng)矢量的空間分辨率為1.05 mm×1.05 mm.

進(jìn)行PIV測(cè)試試驗(yàn)的同時(shí)，使用CF0320-500型加速度傳感器和DH5922型數(shù)據(jù)采集器對(duì)楔形體入水過(guò)程中的加速度進(jìn)行測(cè)試，加速度計(jì)的量程為0～1 000 m/s2，采集頻率為 4 000 Hz.

2 數(shù)值模擬方案

本文應(yīng)用STAR-CCM+對(duì)相應(yīng)試驗(yàn)工況下的楔形體入水過(guò)程進(jìn)行數(shù)值模擬.由于入水過(guò)程中的流速遠(yuǎn)小于聲速，因此可將流動(dòng)視為不可壓縮流動(dòng).設(shè)空氣和水互不相溶，其流動(dòng)滿(mǎn)足質(zhì)量和動(dòng)量守恒定律：

(1)

(2)

式中：u為流體(水和空氣)的速度矢量；ρ為流體密度，大小為998 kg/m3；t為時(shí)間；p為壓力；μ為動(dòng)力黏性系數(shù)，大小為0.889×10-5Pa·s；g為重力加速度.采用k-ε湍流模型對(duì)上述方程封閉并求解，采用速度和壓力場(chǎng)耦合求解，其中對(duì)流項(xiàng)計(jì)算采用二階迎風(fēng)格式，耗散項(xiàng)計(jì)算采用二階中心差分格式.采用VOF法處理自由液面兩相流動(dòng)問(wèn)題，采用重疊網(wǎng)格技術(shù)來(lái)實(shí)現(xiàn)楔形體入水過(guò)程中運(yùn)動(dòng)的模擬.

3 圖像處理與基于PIV的壓力重構(gòu)方法

對(duì)楔形體入水過(guò)程的流場(chǎng)進(jìn)行PIV分析時(shí)，PIV原始圖像包含固體和氣體區(qū)域，會(huì)導(dǎo)致在流固和氣液交界面附近處的粒子圖像相關(guān)程度較低，可能會(huì)產(chǎn)生錯(cuò)誤的流場(chǎng)矢量分析結(jié)果.因此,可應(yīng)用 Radon 變換對(duì)采集的多張?jiān)糚IV圖像進(jìn)行圖像掩膜，依次去除圖像中的固體部分和空氣部分，如圖4所示.

在得到楔形體入水過(guò)程的瞬態(tài)流場(chǎng)矢量信息后，依據(jù)N-S方程可對(duì)流場(chǎng)壓力進(jìn)行重構(gòu)，進(jìn)一步計(jì)算入水過(guò)程中的砰擊載荷.對(duì)于平面二維流動(dòng)，用u、v代替速度矢量u在x、y方向上速度分量，從歐拉角度將式(2)中物質(zhì)導(dǎo)數(shù)展開(kāi)，并忽略體積力項(xiàng)，僅考慮砰擊時(shí)的水動(dòng)壓力，得到方程組：

(3)

(4)

式中：υ為運(yùn)動(dòng)黏性系數(shù)，對(duì)于楔形體入水問(wèn)題一般可以忽略[4,5,19].采用中心差分方法對(duì)式(3)和(4)進(jìn)行離散化差分，得到流場(chǎng)壓力梯度.隨后，依據(jù)壓力梯度進(jìn)行空間積分獲得全場(chǎng)壓力.如圖5所示，圖中(i,j)為PIV數(shù)據(jù)在x,y方向的坐標(biāo)位置.由于砰擊過(guò)程中存在楔形體的運(yùn)動(dòng)邊界和堆積區(qū)域，根據(jù)未擾動(dòng)的自由液面將流域劃分為區(qū)域1和區(qū)域2，對(duì)于區(qū)域1內(nèi)壓力采用多路徑壓力積分方法[21]，令A(yù)點(diǎn)的水動(dòng)壓力為零，由邊界AB和邊界BC向內(nèi)積分獲得，對(duì)于區(qū)域2內(nèi)壓力采用空間侵蝕方法[18]獲得.

4 結(jié)果與分析

4.1 楔形體運(yùn)動(dòng)響應(yīng)分析

對(duì)比不同楔形體在入水過(guò)程中加速度響應(yīng)試驗(yàn)結(jié)果與數(shù)值模擬結(jié)果，如圖6所示，圖中a為加速度.可以看出，發(fā)生砰擊的前20 ms時(shí)間內(nèi)，加速度試驗(yàn)測(cè)試結(jié)果與數(shù)值結(jié)果吻合良好，并且底升角為20°、25°、30°和35°楔形體入水過(guò)程中的實(shí)測(cè)加速度峰值依次約為205 m/s2、164 m/s2、126 m/s2和101 m/s2，出現(xiàn)的時(shí)間依次約在2.5 ms、3.0 ms、4 ms和5 ms，表明隨著楔形體底升角的不斷增大，砰擊的劇烈程度逐漸降低，砰擊加速度峰值的發(fā)生時(shí)間逐漸往后推移.

根據(jù)PIV原始圖像的Radon變換檢測(cè)出楔形體壁面邊界特征可得入水過(guò)程中的楔形體位移，如圖7所示，圖中h為位移，表示楔形體底端與未擾動(dòng)自由液面間的距離.令自由液面位置的坐標(biāo)為 (0,0)，但是在入水的前2 ms內(nèi)位移并未給出，這是由于在初始時(shí)刻，進(jìn)入水中的楔形體部分較小，邊界檢測(cè)位置的誤差較大.由圖7可知，試驗(yàn)測(cè)試進(jìn)行邊界特征檢測(cè)出的楔形體位移值具有良好的穩(wěn)定性，且與數(shù)值結(jié)果匹配度高，表明本文中PIV圖像邊界特征檢測(cè)算法具有較好的準(zhǔn)確性.

4.2 楔形體入水過(guò)程流場(chǎng)分析

圖8所示為底升角分別為20°、25°、30°及35°楔形體入水過(guò)程中t=10 ms時(shí)的流場(chǎng)速度云圖分布.可以看出，整體的速度分布以及局部堆積區(qū)域中細(xì)節(jié)的速度分布對(duì)于不同底升角楔形體的砰擊流場(chǎng)，PIV結(jié)果與數(shù)值結(jié)果均具有良好的一致性，表明本文PIV試驗(yàn)具有良好的精確性.在所有結(jié)果中，速度峰值均出現(xiàn)在堆積區(qū)域頂端，即射流的根部，并且在堆積區(qū)域的速度等值線分布較密，表明該區(qū)域具有較大的速度梯度.t=10 ms時(shí)，20°、25°、30°和35°底升角楔形體的最大流速均約為2 m/s,，并且整體的速度分布大致相似，表明在該范圍內(nèi)底升角對(duì)楔形體入水過(guò)程中的流動(dòng)結(jié)構(gòu)影響較小.隨著底升角由20°逐漸增加到35°，速度分布范圍逐漸收縮，0速度等值線0.2 m/s輪廓的橫坐標(biāo)約由0.11 m收縮至0.09 m，這與入水過(guò)程中的加速度曲線的趨勢(shì)對(duì)應(yīng).對(duì)于相同質(zhì)量的楔形體，底升角越小，產(chǎn)生的砰擊加速度越大，楔形體向流體傳遞的能量越多，從而導(dǎo)致流場(chǎng)流動(dòng)分布范圍越廣.但在射流區(qū)域的速度分布中具有較大差別，主要是由于射流中的示蹤粒子無(wú)法捕捉，且堆積區(qū)域頂端粒子會(huì)不斷進(jìn)入射流中，導(dǎo)致對(duì)該區(qū)域流速的低估.

4.3 壓力重構(gòu)算法驗(yàn)證

在基于TR-PIV進(jìn)行入水砰擊壓力重構(gòu)前，需驗(yàn)證所提出壓力重構(gòu)方案的準(zhǔn)確性，選擇底升角為25°楔形體入水第10 ms時(shí)刻的流場(chǎng)，利用CFD數(shù)據(jù)對(duì)砰擊壓力場(chǎng)進(jìn)行重構(gòu)，探究不同網(wǎng)格間距(Δl)和時(shí)間步長(zhǎng)(Δt)對(duì)壓力重構(gòu)結(jié)果的影響.圖9(a)為CFD計(jì)算的原始流場(chǎng)壓力云圖,將CFD數(shù)據(jù)插入大小為120×160單元網(wǎng)格內(nèi)，對(duì)應(yīng)的網(wǎng)格間距 Δl=1 mm，并構(gòu)造與PIV系統(tǒng)時(shí)間解析能力相同的時(shí)間步長(zhǎng)Δt=0.2 ms，依據(jù)上述方案重構(gòu)出的砰擊壓力云圖如圖9(b)所示.該設(shè)置與TR-PIV測(cè)試結(jié)果(Δl=1.05 mm和Δt=0.2 ms)具有較好的相似性，隨后保持Δl=1 mm不變，分別將時(shí)間步長(zhǎng)倍增至Δt=0.4 ms和Δt=0.8 ms，進(jìn)行壓力重構(gòu)，如圖9(c)和9(d)所示.最后保持Δt=0.2 ms不變，分別倍增網(wǎng)格間距至Δl=2 m和Δl=4 m，即網(wǎng)格大小分別為60×80和30×40，進(jìn)行壓力重構(gòu)，結(jié)果如圖9(e)和9(f)所示.對(duì)比CFD計(jì)算結(jié)果與重構(gòu)出的流場(chǎng)壓力，可知不論是在壓力結(jié)構(gòu)分布和堆積區(qū)域，壓力峰值的預(yù)測(cè)均具有非常良好的一致性，表明本文所提出壓力重構(gòu)方案的準(zhǔn)確性，即使在較大的網(wǎng)格間距和空間步長(zhǎng)情況下，該方案依舊能較好地重構(gòu)出瞬時(shí)砰擊壓力，但是由于截?cái)嗾`差的存在，會(huì)在一定程度上低估砰擊壓力峰值.

4.4 楔形體入水過(guò)程砰擊壓力與載荷

采用上述壓力重構(gòu)方案對(duì)楔形體入水過(guò)程中的瞬時(shí)砰擊壓力進(jìn)行重構(gòu)，圖10所示為底升角分別為20°、25°、30°和35°楔形體入水過(guò)程中t=10 ms時(shí)的砰擊壓力云圖，圖中cp為無(wú)量綱化的壓力系數(shù)，cp=p/(0.5ρv′2) (v′為楔形體速度).首先，對(duì)于不同底升角楔形體的瞬時(shí)砰擊壓力場(chǎng)，其重構(gòu)結(jié)果與CFD結(jié)果吻合十分良好，表明本文中的壓力重構(gòu)方案針對(duì)楔形體入水問(wèn)題具有良好的準(zhǔn)確性和一定的普適性.其次，隨著底升角逐漸增大，砰擊壓力系數(shù)峰值隨之降低，峰值分別約為11.1、7.9、5.8及3.7，這與楔形體加速度曲線體現(xiàn)的砰擊劇烈程度相吻合.最后，不同楔形體的砰擊壓力分布與其速度云圖分布特點(diǎn)十分相似，堆積區(qū)域始終是壓力峰值與速度峰值的所在，這是由于該區(qū)域是射流形成的根部，速度梯度大，流動(dòng)的物質(zhì)導(dǎo)數(shù)較大等原因引起的.但在遠(yuǎn)場(chǎng)流域中，PIV重構(gòu)的壓力等值線輪廓與數(shù)值結(jié)果存在一定偏差，這主要是由于壓力積分過(guò)程中原始PIV測(cè)試誤差會(huì)沿積分路徑擴(kuò)大，對(duì)低壓區(qū)域產(chǎn)生顯著影響，而對(duì)關(guān)鍵的堆積區(qū)域影響較小.

圖11所示為不同楔形體入水過(guò)程中t=10 ms時(shí)的壁面砰擊壓力曲線，圖中橫坐標(biāo)為y/ξ，ξ為入水深度，y/ξ=0對(duì)應(yīng)于楔形體底端位置.楔形體壁面處的壓力重構(gòu)結(jié)果依然與CFD結(jié)果吻合良好.值得注意的是，隨著楔形體底升角由20°增至35°，基于TR-PIV重構(gòu)所預(yù)報(bào)的壁面砰擊壓力峰值與CFD預(yù)測(cè)結(jié)果之間的差值逐漸降低，重構(gòu)方案的壓力峰值預(yù)測(cè)更為準(zhǔn)確.這主要是由于隨著底升角增加，堆積區(qū)域相應(yīng)擴(kuò)大，該區(qū)域內(nèi)的粒子數(shù)目與PIV流場(chǎng)矢量增多，PIV測(cè)試誤差相對(duì)較小，同時(shí)流場(chǎng)矢量基數(shù)的增加也會(huì)使得重構(gòu)誤差進(jìn)一步降低[15,18]，從而提高基于TR-PIV重構(gòu)的精度.

除了砰擊壓力，入水過(guò)程中的砰擊載荷也可以通過(guò)基于TR-PIV重構(gòu)的砰擊壓力沿楔形體壁面積分而獲得.圖12所示為不同楔形體入水過(guò)程中t<20 ms時(shí)的單位寬度砰擊載荷f的演變趨勢(shì).對(duì)于不同楔形體，通過(guò)TR-PIV進(jìn)行砰擊載荷評(píng)估，結(jié)果依然與CFD結(jié)果吻合良好，并且砰擊載荷曲線與加速度曲線反應(yīng)的趨勢(shì)十分類(lèi)似，這是動(dòng)量定理的直觀體現(xiàn).另外，隨著時(shí)間的推移，PIV重構(gòu)結(jié)果存在明顯的振蕩現(xiàn)象，這可能是由于試驗(yàn)過(guò)程中楔形體運(yùn)動(dòng)的不穩(wěn)定性隨著砰擊深入發(fā)展逐漸增大所導(dǎo)致的.

5 結(jié)論

應(yīng)用高頻響TR-PIV技術(shù)對(duì)不同底升角的楔形體自由入水過(guò)程中的瞬時(shí)流場(chǎng)進(jìn)行測(cè)試，同時(shí)應(yīng)用CFD技術(shù)進(jìn)行了相應(yīng)工況下的數(shù)值模擬，分析了入水過(guò)程中流場(chǎng)特性以及基于TR-PIV的砰擊壓力重構(gòu)方法，結(jié)論如下：

(1) 不同楔形體入水過(guò)程中加速度與位移的獨(dú)立多次試驗(yàn)結(jié)果具有良好的一致性，且與數(shù)值結(jié)果吻合良好，隨著楔形體的底升角不斷增大，砰擊劇烈程度逐漸降低，砰擊加速度峰值的發(fā)生時(shí)間逐漸往后推移.

(2) 楔形體入水流場(chǎng)的TR-PIV測(cè)試結(jié)果與數(shù)值結(jié)果吻合良好，在不同楔形體測(cè)試結(jié)果中，速度峰值均出現(xiàn)在堆積區(qū)域頂端即射流的根部，并且底升角對(duì)楔形體入水過(guò)程中的流動(dòng)結(jié)構(gòu)分布影響較小.

(3) 采用本文中的壓力重構(gòu)方案，利用CFD數(shù)據(jù)重構(gòu)砰擊壓力場(chǎng)，原始數(shù)值結(jié)果與重構(gòu)結(jié)果在壓力分布和堆積區(qū)域內(nèi)壓力峰值預(yù)測(cè)上均具有良好的一致性，即使在較大的網(wǎng)格間距和空間步長(zhǎng)情況下，依然有較好的重構(gòu)結(jié)果，但是由于截?cái)嗾`差的存在，會(huì)在一定程度上低估砰擊壓力峰值.

(4) 基于TR-PIV重構(gòu)不同底升角楔形體的砰擊壓力場(chǎng)結(jié)果與數(shù)值結(jié)果吻合良好，本文中的壓力重構(gòu)方案針對(duì)楔形體入水問(wèn)題具有良好的準(zhǔn)確性和一定的普適性.隨著底升角增大，其砰擊壓力系數(shù)峰值降低.此外，通過(guò)TR-PIV間接評(píng)估砰擊載荷的結(jié)果依舊與數(shù)值結(jié)果吻合良好，其砰擊載荷曲線與加速度曲線趨勢(shì)一致.