佘文軒,郭春雨,周廣利,吳鐵成,徐 鵬
(哈爾濱工程大學 船舶工程學院,哈爾濱 150001)
物體入水問題在船舶與海洋工程、航空航天等領域具有廣泛的工程應用背景和科學研究價值,由于其涉及氣、液、固三者之間的耦合作用,近自由液面處的強非線性及強非定常性流動給研究者帶來較大的技術挑戰[1].自Wagner[2]進行入水問題的開創性研究后,出現了大量具有針對性和專業性的入水問題研究手段與分析結果[3-5].
船舶航行過程中,由于較大的縱搖和垂蕩運動會導致危害性較大的底部砰擊作用[6],大多數學者將船體橫剖面簡化為楔形體模型進行深入分析[7-9].Wu等[4]應用邊界元法進行二維楔形體自由入水的數值模擬,并輔以試驗相驗證.Wang等[7]分析了楔形體入水過程中砰擊階段、過渡階段以及空穴閉合階段的流體力學特性.隨著計算機技術迅速發展,多種數值算法被提出并應用于入水問題研究中,如光滑粒子(SPH)法[10]、流體體積(VOF)法[11].但在實驗研究方面,大量學者仍采用加速度、位移、拉力、壓力傳感器等傳統實驗技術手段[4,9],該方法存在安裝繁瑣、測試物理量單一、單點接觸式測控等不足的問題.20世紀80年代發展起來的粒子圖像測試(PIV)技術能夠無干擾地測量全局流場的瞬態速度信息,不僅為數值方法提供了更為全面的數據驗證,還為流場壓力分布、渦量與湍流度等關鍵流動參數的間接測量提供了可能[12].張志榮等[13-14]對楔形體入水時的砰擊瞬態流場進行PIV測試.對于入水問題,研究者往往更關注入水過程中的砰擊壓力,Oudheusden等[15-16]總結了多種基于PIV的壓力重構方法,并指出壓力重構結果的準確性取決于特定的流動問題,目前尚無最優的解決方案.特別是對于具有自由液面遷移和物體邊界運動的入水問題,一些方法的適用性與精確性有待證實.Nila[17]、Panciroli[18]和Jalalisendi[19]等采用PIV技術對楔形體入水的流場進行測試和砰擊壓力的重構,表明了基于PIV技術進行楔形體入水砰擊壓力間接評估的可行性,但上述研究僅采用單一的楔形體模型進行砰擊壓力重構方法的驗證,其適用性有待進一步探究.
本文應用高頻響的時間解析PIV(TR-PIV)技術對不同底升角θ楔形體入水過程中的流場進行測試,并依據TR-PIV瞬態測試結果進行入水過程中瞬時流場壓力分布以及砰擊載荷重構,同時進行相應試驗工況下的數值模擬,對比分析基于TR-PIV數據進行砰擊壓力重構的準確性和適用性,探討楔形體入水過程中的細節流動特征.
試驗裝置如圖1所示,鋁型材框架、直線滑軌以及亞克力材質的透明水箱構成試驗裝置的主要部分,水箱長800 mm,寬500 mm,高500 mm,由鋁型材框架支撐,長800 mm的直線滑軌固定在水箱的正上方.滑軌中間設置有電磁裝置,楔形體上端設有連接滑塊,楔形體與連接滑塊組成1個整體,受電磁裝置的控制,能夠在滑軌間自由滑動,并以垂直的速度入水.
圖2所示為4個底升角分別為20°、25°、30° 及35° 的楔形體模型,由剛性聚乳酸(PLA)材料經3D打印制成,寬度為200 mm,用配重調整模型質量,使4個楔形體模型與連接滑塊的總重均為0.6 kg,從距離水面30 cm處自由落下.忽略空氣阻力的影響,根據理論得到其對應的入水初速度為2.425 m/s,激光片光照射楔形體的中間位置切面,可認為中間切面處的流動近似為二維流動[18-19].試驗時水溫為20 ℃,水的密度為998.16 kg/m3,重力加速度g為9.8 m/s2.
本文中所采用的TR-PIV系統主要由1臺連續激光器、NAC Memrecam HX-6高速CMOS相機以及計算機組成.10 W釔鋁石榴石晶體連續激光器(波長532 nm)作為試驗光源,測試區域的片光厚度約為1 mm.相機的內存為8 GB,試驗時相機的空間分辨率為 1 280 像素×1 000 像素,采集速率為 5 000 Hz,圖像深度為16 bit,楔形體底端初始接觸水面的時刻設置為0時刻.
楔形體入水過程具有流動對稱性,進行PIV測試時,只觀測入水過程中的一半流場區域.采用密度約為1.03 g/mm3、粒徑約20 μm的聚酰胺微珠PSP-20作為實驗示蹤粒子.網格大小為10 mm×10 mm 的標定板用于PIV系統標定.PIV系統測量區域的大小約為165 mm×130 mm,對應于7.63 像素/mm.記錄每次入水砰擊過程前20 ms的流場信息,即采集101張PIV粒子圖像.粒子圖像的互相關計算采用基于MATLAB平臺的開源軟件PIVlab,使用基于快速傅里葉變換(FFT)的多重網格迭代技術,設置有64像素×64像素、32像素×32像素、16像素×16像素多重判讀窗口,相鄰窗口重疊率為50%,采用3點高斯亞像素插值進行互相關峰值擬合,擬合結果精度約為0.1像素[20],流場矢量的空間分辨率為1.05 mm×1.05 mm.
進行PIV測試試驗的同時,使用CF0320-500型加速度傳感器和DH5922型數據采集器對楔形體入水過程中的加速度進行測試,加速度計的量程為0~1 000 m/s2,采集頻率為 4 000 Hz.
本文應用STAR-CCM+對相應試驗工況下的楔形體入水過程進行數值模擬.由于入水過程中的流速遠小于聲速,因此可將流動視為不可壓縮流動.設空氣和水互不相溶,其流動滿足質量和動量守恒定律:
(1)
(2)
式中:u為流體(水和空氣)的速度矢量;ρ為流體密度,大小為998 kg/m3;t為時間;p為壓力;μ為動力黏性系數,大小為0.889×10-5Pa·s;g為重力加速度.采用k-ε湍流模型對上述方程封閉并求解,采用速度和壓力場耦合求解,其中對流項計算采用二階迎風格式,耗散項計算采用二階中心差分格式.采用VOF法處理自由液面兩相流動問題,采用重疊網格技術來實現楔形體入水過程中運動的模擬.

對楔形體入水過程的流場進行PIV分析時,PIV原始圖像包含固體和氣體區域,會導致在流固和氣液交界面附近處的粒子圖像相關程度較低,可能會產生錯誤的流場矢量分析結果.因此,可應用 Radon 變換對采集的多張原始PIV圖像進行圖像掩膜,依次去除圖像中的固體部分和空氣部分,如圖4所示.
在得到楔形體入水過程的瞬態流場矢量信息后,依據N-S方程可對流場壓力進行重構,進一步計算入水過程中的砰擊載荷.對于平面二維流動,用u、v代替速度矢量u在x、y方向上速度分量,從歐拉角度將式(2)中物質導數展開,并忽略體積力項,僅考慮砰擊時的水動壓力,得到方程組:
(3)
(4)
式中:υ為運動黏性系數,對于楔形體入水問題一般可以忽略[4,5,19].采用中心差分方法對式(3)和(4)進行離散化差分,得到流場壓力梯度.隨后,依據壓力梯度進行空間積分獲得全場壓力.如圖5所示,圖中(i,j)為PIV數據在x,y方向的坐標位置.由于砰擊過程中存在楔形體的運動邊界和堆積區域,根據未擾動的自由液面將流域劃分為區域1和區域2,對于區域1內壓力采用多路徑壓力積分方法[21],令A點的水動壓力為零,由邊界AB和邊界BC向內積分獲得,對于區域2內壓力采用空間侵蝕方法[18]獲得.
對比不同楔形體在入水過程中加速度響應試驗結果與數值模擬結果,如圖6所示,圖中a為加速度.可以看出,發生砰擊的前20 ms時間內,加速度試驗測試結果與數值結果吻合良好,并且底升角為20°、25°、30°和35°楔形體入水過程中的實測加速度峰值依次約為205 m/s2、164 m/s2、126 m/s2和101 m/s2,出現的時間依次約在2.5 ms、3.0 ms、4 ms和5 ms,表明隨著楔形體底升角的不斷增大,砰擊的劇烈程度逐漸降低,砰擊加速度峰值的發生時間逐漸往后推移.
根據PIV原始圖像的Radon變換檢測出楔形體壁面邊界特征可得入水過程中的楔形體位移,如圖7所示,圖中h為位移,表示楔形體底端與未擾動自由液面間的距離.令自由液面位置的坐標為 (0,0),但是在入水的前2 ms內位移并未給出,這是由于在初始時刻,進入水中的楔形體部分較小,邊界檢測位置的誤差較大.由圖7可知,試驗測試進行邊界特征檢測出的楔形體位移值具有良好的穩定性,且與數值結果匹配度高,表明本文中PIV圖像邊界特征檢測算法具有較好的準確性.
圖8所示為底升角分別為20°、25°、30°及35°楔形體入水過程中t=10 ms時的流場速度云圖分布.可以看出,整體的速度分布以及局部堆積區域中細節的速度分布對于不同底升角楔形體的砰擊流場,PIV結果與數值結果均具有良好的一致性,表明本文PIV試驗具有良好的精確性.在所有結果中,速度峰值均出現在堆積區域頂端,即射流的根部,并且在堆積區域的速度等值線分布較密,表明該區域具有較大的速度梯度.t=10 ms時,20°、25°、30°和35°底升角楔形體的最大流速均約為2 m/s,,并且整體的速度分布大致相似,表明在該范圍內底升角對楔形體入水過程中的流動結構影響較小.隨著底升角由20°逐漸增加到35°,速度分布范圍逐漸收縮,0速度等值線0.2 m/s輪廓的橫坐標約由0.11 m收縮至0.09 m,這與入水過程中的加速度曲線的趨勢對應.對于相同質量的楔形體,底升角越小,產生的砰擊加速度越大,楔形體向流體傳遞的能量越多,從而導致流場流動分布范圍越廣.但在射流區域的速度分布中具有較大差別,主要是由于射流中的示蹤粒子無法捕捉,且堆積區域頂端粒子會不斷進入射流中,導致對該區域流速的低估.
在基于TR-PIV進行入水砰擊壓力重構前,需驗證所提出壓力重構方案的準確性,選擇底升角為25°楔形體入水第10 ms時刻的流場,利用CFD數據對砰擊壓力場進行重構,探究不同網格間距(Δl)和時間步長(Δt)對壓力重構結果的影響.圖9(a)為CFD計算的原始流場壓力云圖,將CFD數據插入大小為120×160單元網格內,對應的網格間距 Δl=1 mm,并構造與PIV系統時間解析能力相同的時間步長Δt=0.2 ms,依據上述方案重構出的砰擊壓力云圖如圖9(b)所示.該設置與TR-PIV測試結果(Δl=1.05 mm和Δt=0.2 ms)具有較好的相似性,隨后保持Δl=1 mm不變,分別將時間步長倍增至Δt=0.4 ms和Δt=0.8 ms,進行壓力重構,如圖9(c)和9(d)所示.最后保持Δt=0.2 ms不變,分別倍增網格間距至Δl=2 m和Δl=4 m,即網格大小分別為60×80和30×40,進行壓力重構,結果如圖9(e)和9(f)所示.對比CFD計算結果與重構出的流場壓力,可知不論是在壓力結構分布和堆積區域,壓力峰值的預測均具有非常良好的一致性,表明本文所提出壓力重構方案的準確性,即使在較大的網格間距和空間步長情況下,該方案依舊能較好地重構出瞬時砰擊壓力,但是由于截斷誤差的存在,會在一定程度上低估砰擊壓力峰值.
采用上述壓力重構方案對楔形體入水過程中的瞬時砰擊壓力進行重構,圖10所示為底升角分別為20°、25°、30°和35°楔形體入水過程中t=10 ms時的砰擊壓力云圖,圖中cp為無量綱化的壓力系數,cp=p/(0.5ρv′2) (v′為楔形體速度).首先,對于不同底升角楔形體的瞬時砰擊壓力場,其重構結果與CFD結果吻合十分良好,表明本文中的壓力重構方案針對楔形體入水問題具有良好的準確性和一定的普適性.其次,隨著底升角逐漸增大,砰擊壓力系數峰值隨之降低,峰值分別約為11.1、7.9、5.8及3.7,這與楔形體加速度曲線體現的砰擊劇烈程度相吻合.最后,不同楔形體的砰擊壓力分布與其速度云圖分布特點十分相似,堆積區域始終是壓力峰值與速度峰值的所在,這是由于該區域是射流形成的根部,速度梯度大,流動的物質導數較大等原因引起的.但在遠場流域中,PIV重構的壓力等值線輪廓與數值結果存在一定偏差,這主要是由于壓力積分過程中原始PIV測試誤差會沿積分路徑擴大,對低壓區域產生顯著影響,而對關鍵的堆積區域影響較小.
圖11所示為不同楔形體入水過程中t=10 ms時的壁面砰擊壓力曲線,圖中橫坐標為y/ξ,ξ為入水深度,y/ξ=0對應于楔形體底端位置.楔形體壁面處的壓力重構結果依然與CFD結果吻合良好.值得注意的是,隨著楔形體底升角由20°增至35°,基于TR-PIV重構所預報的壁面砰擊壓力峰值與CFD預測結果之間的差值逐漸降低,重構方案的壓力峰值預測更為準確.這主要是由于隨著底升角增加,堆積區域相應擴大,該區域內的粒子數目與PIV流場矢量增多,PIV測試誤差相對較小,同時流場矢量基數的增加也會使得重構誤差進一步降低[15,18],從而提高基于TR-PIV重構的精度.
除了砰擊壓力,入水過程中的砰擊載荷也可以通過基于TR-PIV重構的砰擊壓力沿楔形體壁面積分而獲得.圖12所示為不同楔形體入水過程中t<20 ms時的單位寬度砰擊載荷f的演變趨勢.對于不同楔形體,通過TR-PIV進行砰擊載荷評估,結果依然與CFD結果吻合良好,并且砰擊載荷曲線與加速度曲線反應的趨勢十分類似,這是動量定理的直觀體現.另外,隨著時間的推移,PIV重構結果存在明顯的振蕩現象,這可能是由于試驗過程中楔形體運動的不穩定性隨著砰擊深入發展逐漸增大所導致的.
應用高頻響TR-PIV技術對不同底升角的楔形體自由入水過程中的瞬時流場進行測試,同時應用CFD技術進行了相應工況下的數值模擬,分析了入水過程中流場特性以及基于TR-PIV的砰擊壓力重構方法,結論如下:
(1) 不同楔形體入水過程中加速度與位移的獨立多次試驗結果具有良好的一致性,且與數值結果吻合良好,隨著楔形體的底升角不斷增大,砰擊劇烈程度逐漸降低,砰擊加速度峰值的發生時間逐漸往后推移.
(2) 楔形體入水流場的TR-PIV測試結果與數值結果吻合良好,在不同楔形體測試結果中,速度峰值均出現在堆積區域頂端即射流的根部,并且底升角對楔形體入水過程中的流動結構分布影響較小.
(3) 采用本文中的壓力重構方案,利用CFD數據重構砰擊壓力場,原始數值結果與重構結果在壓力分布和堆積區域內壓力峰值預測上均具有良好的一致性,即使在較大的網格間距和空間步長情況下,依然有較好的重構結果,但是由于截斷誤差的存在,會在一定程度上低估砰擊壓力峰值.
(4) 基于TR-PIV重構不同底升角楔形體的砰擊壓力場結果與數值結果吻合良好,本文中的壓力重構方案針對楔形體入水問題具有良好的準確性和一定的普適性.隨著底升角增大,其砰擊壓力系數峰值降低.此外,通過TR-PIV間接評估砰擊載荷的結果依舊與數值結果吻合良好,其砰擊載荷曲線與加速度曲線趨勢一致.