魏英杰,楊 柳,王 聰,夏維學(xué),李佳川
(1.哈爾濱工業(yè)大學(xué),哈爾濱 150001;2.天津航海儀器研究所,天津 300131)
運(yùn)動(dòng)體從空氣中穿越自由液面進(jìn)入水中所引起的一系列問(wèn)題統(tǒng)稱為入水流動(dòng)問(wèn)題。對(duì)于剛性運(yùn)動(dòng)體的入水問(wèn)題而言,包括入水噴濺[1-2]、入水空泡[3]、入水抨擊[4-6]等問(wèn)題,已經(jīng)有100 多年的研究歷史。但是對(duì)于超彈性運(yùn)動(dòng)體的入水問(wèn)題而言,由于涉及材料大變形的流固耦合問(wèn)題,具有相當(dāng)?shù)难芯侩y度,近年來(lái)才開(kāi)始有相關(guān)研究。
球體作為一種典型的運(yùn)動(dòng)體,由于其具有對(duì)稱性且運(yùn)動(dòng)穩(wěn)定性較好,已成為入水問(wèn)題的首選研究對(duì)象。Abraham 等通過(guò)數(shù)值方法分析了小球垂直入水空泡流動(dòng)特性[7],以及阻力系數(shù)隨表面張力、流動(dòng)狀態(tài)和入水速度的變化規(guī)律。Aristoff等研究了流體屬性、球體密度以及疏水性對(duì)球體入水空泡的影響[8-10],獲得了入水空泡形態(tài)與韋伯?dāng)?shù)、邦德數(shù)的關(guān)系,發(fā)現(xiàn)了空泡深閉合時(shí)間與密度無(wú)關(guān)。Johnson[11-12]、Belden[13]、Moxnes[14]開(kāi)展了球體入水運(yùn)動(dòng)參數(shù)和液體介質(zhì)屬性對(duì)入水跳彈的影響研究,獲得了運(yùn)動(dòng)參數(shù)對(duì)跳彈現(xiàn)象的影響規(guī)律,以及跳彈現(xiàn)象的臨界入射角度及其阻力系數(shù)。李佳川等開(kāi)展了加熱球體入水過(guò)程的實(shí)驗(yàn)和數(shù)值仿真研究[15],分析了球體溫度對(duì)空泡形態(tài)以及空泡的形成、發(fā)展閉合潰滅的影響規(guī)律。Truscott等采用實(shí)驗(yàn)方法研究了旋轉(zhuǎn)小球入水后的流場(chǎng)特性[16-17],發(fā)現(xiàn)了旋轉(zhuǎn)小球入水后呈現(xiàn)馬格努斯效應(yīng)的曲線軌跡。夏維學(xué)等基于VOF多相流模型分析了低弗勞德數(shù)下高速旋轉(zhuǎn)的小球入水空泡及流體動(dòng)力特性[18]。
大多數(shù)經(jīng)典的入水問(wèn)題研究都忽略了運(yùn)動(dòng)體入水過(guò)程的流固耦合作用,然而近年來(lái)由于彈性材料的發(fā)展,人們也開(kāi)始逐漸開(kāi)始重視入水過(guò)程中的流固耦合問(wèn)題。Yan等[19]采用實(shí)驗(yàn)方法研究了平板入水后的抨擊壓力和水動(dòng)力彈性振動(dòng),發(fā)現(xiàn)沖擊載荷作用下,水動(dòng)力彈性振動(dòng)使紊流氣-水場(chǎng)的壓力周期性下降。Xie等[20]結(jié)合實(shí)驗(yàn)和流固耦合數(shù)值方法對(duì)復(fù)合材料船體結(jié)構(gòu)入水開(kāi)展了研究,通過(guò)對(duì)結(jié)構(gòu)固有頻率的分析,討論了其對(duì)水動(dòng)力壓力和應(yīng)力響應(yīng)的影響。Michael等[21]采用數(shù)值模擬和實(shí)驗(yàn)方法,研究了彈性球體在水面的跳躍問(wèn)題,分析了彈性球體水面上跳躍的內(nèi)在機(jī)理。
對(duì)于超彈性球體的研究問(wèn)題,任九生等[22]基于有限變形動(dòng)力學(xué)理論,分析了受到均勻拉伸載荷作用下可壓縮超彈性材料球殼的有限振動(dòng)問(wèn)題,得到了球殼振動(dòng)的時(shí)程和相圖,并且證明了可壓縮超彈性球殼的振動(dòng)具有擬周期性和非線性。袁學(xué)剛等研究了橫向各向同性的不可壓縮超彈性球表面受周期性階梯拉伸載荷時(shí)孔洞的形成和運(yùn)動(dòng)問(wèn)題[23],結(jié)果表明拉伸載荷超過(guò)某一臨界值時(shí)在球體表面會(huì)形成孔洞,并且在周期性的階躍載荷下,給出了周期振蕩形成孔洞的存在條件。Peyraut等開(kāi)展了受靜水壓力下超彈性球體的數(shù)值仿真研究[24],發(fā)現(xiàn)了牛頓-拉弗森算法在無(wú)方向保持下容易發(fā)散,且提出最優(yōu)時(shí)間步載荷算法。Fanning和Hurd等采用高速攝像方法開(kāi)展了可變形彈性體球入水實(shí)驗(yàn)[25],研究發(fā)現(xiàn)彈性球體入水后形成了獨(dú)特的嵌套空泡現(xiàn)象,并通過(guò)材料性能和沖擊條件預(yù)測(cè)了球體的變形和振動(dòng)。
通過(guò)上述分析發(fā)現(xiàn),目前球體入水研究主要集中于剛性球體的入水問(wèn)題,同時(shí),超彈性球體的研究主要集中于均勻或者靜止載荷下超彈性球體的結(jié)構(gòu)動(dòng)態(tài)研究,而涉及超彈性球體入水流固耦合問(wèn)題的研究相對(duì)較少。同時(shí),雖然少數(shù)學(xué)者開(kāi)展了可變形彈性體球入水實(shí)驗(yàn),但是由于實(shí)驗(yàn)條件的限制,對(duì)于空泡輪廓的準(zhǔn)確捕捉、流場(chǎng)參數(shù)的分布和應(yīng)力變形的測(cè)量等問(wèn)題暫時(shí)都無(wú)法很好地解決,因此結(jié)合實(shí)驗(yàn)開(kāi)展數(shù)值模擬研究具有重要意義。本文采用實(shí)驗(yàn)與數(shù)值模擬相結(jié)合的方法,開(kāi)展超彈性球體垂直入水空泡流動(dòng)研究,分析超彈性球體入水過(guò)程的流動(dòng)特性。
本文的超彈性球體入水空泡實(shí)驗(yàn)在室溫下進(jìn)行,實(shí)驗(yàn)系統(tǒng)的整體示意圖如圖1 所示。水箱尺寸為1.5 m×0.8 m×0.9 m,與球體相比有足夠大的尺度,故可以忽略邊壁效應(yīng)對(duì)入水空泡形態(tài)的影響。水箱材料均為鋼化玻璃,其中側(cè)壁厚度為10 mm,底面厚度為15 mm,各面之間由酸性硅酮玻璃膠膠合。超彈性球體由釋放裝置在水槽中心位置的正上方無(wú)干擾地釋放,以實(shí)現(xiàn)對(duì)超彈性球體入水速度的控制。
實(shí)驗(yàn)采用Photron FASTCAM SA-X 型高速攝像機(jī)對(duì)超彈性球體垂直入水過(guò)程進(jìn)行拍攝,拍攝幀率為1000 fps。實(shí)驗(yàn)時(shí)在水箱后方設(shè)置平面LED 燈陣作為光源,以保證采光要求。同時(shí)在后方光源和水箱之間豎直布置柔光屏,以提供較柔和的光線,從而能夠較好地捕捉流場(chǎng)細(xì)微結(jié)構(gòu),得到優(yōu)質(zhì)的圖像。本文實(shí)驗(yàn)所捕捉的空泡形態(tài)、位移及變形量都是基于四次以上重復(fù)性實(shí)驗(yàn)結(jié)果獲得的。
圖1 實(shí)驗(yàn)系統(tǒng)整體示意圖Fig.1 Schematic diagram of the experimental system
實(shí)驗(yàn)采用的超彈性球體是由不可壓縮硅橡膠制作而成。超彈性球體的材料參數(shù)為:直徑D =60 mm,密度ρs=988 kg/m3,剪切模量G=12.7 kPa。
本文采用Navier-Stocks 方程描述入水空泡流動(dòng),通過(guò)VOF多相流模型捕捉運(yùn)動(dòng)體入水過(guò)程中的氣液界面。
運(yùn)動(dòng)方程:
連續(xù)性方程:
其中,ui為流體速度,ρ 為流體密度,μ 為流體動(dòng)力黏度,p 為流體壓強(qiáng),gi為重力加速度。
其中,ρl為水的密度;ρg為空氣密度;μl為水的動(dòng)力黏度;μg為空氣的動(dòng)力黏度;αl為液相體積分?jǐn)?shù);αg為氣相體積分?jǐn)?shù)。根據(jù)VOF多相流模型:αl=1代表流域網(wǎng)格單元完全為液相;αl=0代表流域網(wǎng)格單元完全為氣相;0<αl<1代表流域網(wǎng)格單元將出現(xiàn)氣液交界面。
基于虛功原理,采用非線性有限元法分析結(jié)構(gòu)的變形和水動(dòng)力荷載,必須使得物體每個(gè)單元的力和力矩得到平衡。
離散化的平衡方程可以寫成:
其中,PN為外力;IN為內(nèi)力;MNM為質(zhì)量矩陣;uM為單元加速度為由材料慣性引起的力。
內(nèi)力方程可以表示為:
其中,V0為球體的未變形體積;σ 為應(yīng)力張量;βN為彈性體本構(gòu)模型定義的應(yīng)變率。
外力方程可以表示為:
其中NN為單元形狀函數(shù);t 為曲面上的牽引力矢量;F 為球面上的力。
由于結(jié)構(gòu)計(jì)算運(yùn)行是隱式的,因此需要對(duì)矩陣求逆,并在每一增量中求解存在的非線性平衡方程。用牛頓法求解了非線性平衡方程的內(nèi)解。
本文基于徑向基函數(shù)(RBF)的變形插值法和虛功原理的力映射法,開(kāi)展超彈性球體入水流固耦合問(wèn)題的數(shù)值模擬。流場(chǎng)計(jì)算通過(guò)在StarCCM+軟件中采用重疊網(wǎng)格技術(shù)來(lái)實(shí)現(xiàn),流體計(jì)算域尺度為45D×20D×20D,重疊域尺度約為3D×3D×3D。重疊域跟隨超彈性球體在流場(chǎng)計(jì)算域中運(yùn)動(dòng)。計(jì)算域頂端為滯止入口、側(cè)面為壁面,底端為壓力入口,其壓力值及液相體積分?jǐn)?shù)通過(guò)自定義函數(shù)進(jìn)行指定。同時(shí),為保證計(jì)算的準(zhǔn)確性,計(jì)算域的網(wǎng)格劃分方式均采用六面體結(jié)構(gòu)網(wǎng)格,在球體3D 范圍內(nèi)進(jìn)行漸進(jìn)加密,加密區(qū)外圍的網(wǎng)格稍稀疏。流場(chǎng)計(jì)算總體網(wǎng)格數(shù)量約為280 萬(wàn),流場(chǎng)網(wǎng)格流場(chǎng)計(jì)算域如圖2 所示。其中,球體的初始入水速度為u0,自由液面下方球體最低點(diǎn)的深度為yb,初始直徑為D ,變形后直徑為(其定義為通過(guò)幾何中心平行于x 軸與球面交點(diǎn)之間的距離),定義球體變形系數(shù)為
圖2 流場(chǎng)計(jì)算域示意圖Fig.2 Schematic diagram of flow field
結(jié)構(gòu)計(jì)算在Abaqus模擬軟件中采用隱性算法來(lái)實(shí)現(xiàn)。在結(jié)構(gòu)場(chǎng)中定義了彈性體材料、運(yùn)動(dòng)方式。為了與流場(chǎng)網(wǎng)格節(jié)點(diǎn)相對(duì)應(yīng),結(jié)構(gòu)場(chǎng)計(jì)算也采用六面體結(jié)構(gòu)網(wǎng)格,結(jié)構(gòu)單元采用C3D8RH,固體結(jié)構(gòu)網(wǎng)格約130萬(wàn),固體結(jié)構(gòu)網(wǎng)格如圖3所示。
圖3 固體結(jié)構(gòu)網(wǎng)格Fig.3 Grid of solid structure
在流固耦合計(jì)算中,SIMULIA協(xié)同仿真引擎負(fù)責(zé)Abaqus和Star CCM+之間的雙向耦合代碼通信,利用求解器之間的這種通信,即可以獲得整個(gè)流體/固體交界面的完全求解。
在初始入水速度u0=4.7 m/s條件下,針對(duì)直徑D=60 mm,密度ρs=988 kg/m3,剪切模量G=12.7 kPa的不可壓縮硅膠超彈性球體開(kāi)展入水實(shí)驗(yàn)。入水空泡和球體變形的實(shí)驗(yàn)結(jié)果如圖4所示。在入水撞擊階段(t=0 ms),超彈性球體由于經(jīng)歷了從空氣到水的跨介質(zhì)突變而引起空泡與噴濺的形成。同時(shí),球體自身也受到強(qiáng)烈的沖擊,由于超彈性球體剛度較小,與水面發(fā)生沖擊時(shí)發(fā)生了明顯的變形(t=6 ms)。隨著球體入水深度增大,流動(dòng)分離點(diǎn)上方已獲得動(dòng)能的流體在慣性作用下繼續(xù)向外排開(kāi),導(dǎo)致空泡直徑變大,空泡演化形成敞開(kāi)趨勢(shì)。此時(shí)球體變形為一個(gè)與自由液面近似平行的橫向扁狀球體,將水大幅度地排開(kāi),進(jìn)而也形成了圖中所示又寬又短的空泡(t=8 ms)。但是由于超彈性球體內(nèi)部彈性力的存在,使得超彈性球體由原來(lái)近似平行于自由液面的橫向橢球體反彈變形為一個(gè)縱向扁平橢球體(t=18 ms)。然后該縱向橢球體徑向膨脹穿透第一個(gè)空泡壁面,同時(shí),在水動(dòng)力的作用下,縱向扁平橢球體又變形回原來(lái)的橫向扁平橢球體,排開(kāi)附近流體,進(jìn)而出現(xiàn)了第二個(gè)空泡,由于第二個(gè)空泡沒(méi)有大量空氣進(jìn)入,因此第二個(gè)空泡相對(duì)第一個(gè)空泡要小很多,第二個(gè)小空泡被稱作嵌套空泡[22](t=24~40 ms)。該嵌套空泡產(chǎn)生的主要原因是超彈性球體入水后發(fā)生材料震蕩導(dǎo)致球體形態(tài)發(fā)生大變形而產(chǎn)生。隨著超彈性球體不斷下降,球體將穿透嵌套空泡,由于此時(shí)超彈性球體下降過(guò)程被空泡包裹,其所受阻力小,變形不大,小幅度地排開(kāi)水,如水珠一般下降(t=66 ms),直至最后發(fā)生空泡閉合(t=110 ms)。
圖4 入水空泡形態(tài)和球體變形的實(shí)驗(yàn)結(jié)果Fig.4 Experimental results of cavity shape and sphere deformation
在初始入水速度u0=4.7 m/s條件下,針對(duì)直徑D=60 mm,密度ρs=988 kg/m3,剪切模量G=12.7 kPa的超彈性球體入水過(guò)程開(kāi)展數(shù)值模擬研究。入水空泡的數(shù)值模擬結(jié)果如圖5所示。從圖中可以看出,數(shù)值模擬較好地模擬出了實(shí)驗(yàn)中的入水空泡演化過(guò)程。在相同的入水時(shí)刻下,入水空泡形態(tài)、球體的變形行為與實(shí)驗(yàn)結(jié)果具有較好的一致性。
圖5 入水空泡形態(tài)和球體變形的數(shù)值結(jié)果Fig.5 Numerical results of cavity shape and sphere deformation
為了進(jìn)一步驗(yàn)證數(shù)值模擬方法的有效性,選取彈性球體最低點(diǎn)位移作為研究對(duì)象,將位移的數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比,如圖6所示。從圖中可以看出球體位移的數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)具有較好的一致性。
圖6 球體位移Fig.6 Comparison of displacement
為了驗(yàn)證數(shù)值結(jié)果中球體變形模擬的準(zhǔn)確性,將球體變形系數(shù)λ 的數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比,如圖7所示。從圖中可以看出二者在變形周期上稍有差異,但是在變形結(jié)構(gòu)方面符合得很好,最大誤差約為3.5%。
圖7 球體變形系數(shù)Fig.7 Comparison of deformation coefficient
通過(guò)上述關(guān)于入水空泡形態(tài)、球體位移和變形系數(shù)的分析對(duì)比,驗(yàn)證了本文數(shù)值計(jì)算方法的有效性。
為了深入研究超彈性球體入水特有的入水流動(dòng)特性,在入水速度u0=4.7 m/s的條件下,本文采用流固耦合數(shù)值模擬方法,針對(duì)直徑D=60 mm,剪切模量G=12.7 kPa的超彈性球體與具有相同密度和直徑、剪切模量G=337 MPa的剛性球體的入水流動(dòng)特性開(kāi)展對(duì)比研究。
3.3.1 超彈性球體入水空泡特性
為了定量地對(duì)比超彈性球體與剛性球體空泡的尺寸,分別定義空泡特征長(zhǎng)度和特征寬度。其中空泡長(zhǎng)度l取自由液面到空泡與球體接觸點(diǎn)的距離來(lái)表征,空泡寬度分別取空泡口橫向?qū)挾菵z和縱向?qū)挾菵y來(lái)表征,如圖8所示。
圖8 空泡特征長(zhǎng)度和寬度定義Fig.8 Definition of cavity
圖9 給出了超彈性球體與剛性球體無(wú)量綱空泡長(zhǎng)度的數(shù)值結(jié)果。從圖中可以看出,在同一時(shí)刻,超彈性球體的空泡長(zhǎng)度要小于剛性球體的空泡長(zhǎng)度。其原因主要是超彈性球體在入水過(guò)程中由于變形會(huì)導(dǎo)致能量損耗,超彈性球體的變形能量損耗明顯大于剛性球體(如圖10所示),因此入水過(guò)程中超彈性球體速度降低幅度大于剛性球體,進(jìn)而導(dǎo)致空泡長(zhǎng)度縮短。
圖9 超彈性球體與剛性球體的空泡長(zhǎng)度比較Fig.9 Length comparison of cavity for elastic and rigid spheres
圖10 超彈性球體與剛性球體的變形能量損失比較Fig.10 Deformation energy loss of elastic and rigid spheres
圖11 給出了超彈性球體和剛性球體無(wú)量綱空泡橫向?qū)挾扰c縱向?qū)挾鹊臄?shù)值結(jié)果。從圖中可以看出,超彈性球體入水空泡的橫向?qū)挾群涂v向?qū)挾染笥趧傂郧蝮w的空泡寬度。其原因主要是超彈性球體入水后變形嚴(yán)重,從而導(dǎo)致流體向外排開(kāi)液體的幅度越大,故入水后超彈性球體形成的空泡寬度變大。
圖11 超彈性球體與剛性球體的空泡寬度比較Fig.11 Width comparison of cavity for elastic and rigid spheres
其次,本文對(duì)比分析了某一時(shí)刻不同剪切模量條件下,多個(gè)截面處超彈性球體的入水空泡口輪廓,如圖12所示。從圖中可以看出在h/D=0~3的多個(gè)不同截面處,剪切模量為G=12.7 kPa球體空泡口輪廓最大,G=25.3 kPa球體空泡口輪廓其次,G=42.0 kPa球體空泡口輪廓最小。故球體的剪切模量越小,入水空泡輪廓越大,其主要原因是剪切模量越小,入水后變形量就越大,故入水后空泡口輪廓越大。
圖12 不同剪切模量超彈性球體入水空泡輪廓對(duì)比Fig.12 Comparison of cavity profiles for elastic spheres with different material shear moduli
3.3.2 超彈性球體表面壓力特性
為了對(duì)比超彈性球體與剛性球體的表面壓力特性,提取了xoy 平面與球體表面交線的下半部分圓弧線壓力,從而獲得典型時(shí)刻下球體表面的壓力變化規(guī)律,如圖13所示。
從圖中可以看出剛性球體入水后球體表面的壓力隨y 軸位置的變化趨勢(shì)都呈現(xiàn)開(kāi)口向下的拋物線形式,球體頂端受到的正向壓力最大,隨著位移向軸線兩側(cè)移動(dòng)壓力逐漸降低,在球體表面與空泡接觸的位置處壓力達(dá)到零。
從圖中還可以看出超彈性球體入水后球體表面壓力隨y 軸位置的變化在不同時(shí)刻呈現(xiàn)不同的變化趨勢(shì)。當(dāng)t=20 ms時(shí),球體內(nèi)部彈性力使得超彈性球體反彈形成縱向扁平橢球體,此時(shí)球體下表面呈現(xiàn)回彈收縮趨勢(shì),底部形成低壓區(qū),所以較剛性球面的壓力低,且球體赤道處向內(nèi)收縮變形劇烈,該變形會(huì)產(chǎn)生反向運(yùn)動(dòng)與空氣尾流發(fā)生碰撞,使得彈性球體附近的流體流動(dòng)受到較大的擾動(dòng),球體兩側(cè)形成局部擾流,故出現(xiàn)局部負(fù)壓區(qū)域。當(dāng)t=30 ms時(shí),超彈性球體受到較剛性球體大2倍以上的水動(dòng)壓力,使得超彈性球體由縱向扁平橢球體變形為橫向扁平橢球體,進(jìn)而也形成了第一個(gè)嵌套空泡,同時(shí),局部擾流導(dǎo)致球體赤道兩側(cè)出現(xiàn)兩處局部負(fù)壓區(qū)。當(dāng)t=50 ms時(shí),超彈性球體受到不均勻的表面壓力,此時(shí),球體處于第二個(gè)嵌套空泡的形成階段,逐漸由縱向扁平橢球體向橫向扁平橢球體變形,同時(shí),由于局部擾流球體赤道兩側(cè)出現(xiàn)兩處局部負(fù)壓區(qū)。同理,當(dāng)t=70 ms時(shí),超彈性球體受到局部高壓,但球體兩側(cè)因湍流出現(xiàn)負(fù)壓區(qū)。
圖13 超彈性球體和剛性球體的表面壓力分布Fig.13 Surface pressure distribution for elastic and rigid spheres
圖14給出了超彈性球和剛性球體入水流場(chǎng)的壓力分布云圖。從圖中可以看出,剛性球體入水后表面壓力的變化范圍約為0 ~7 k Pa,而超彈性球體表面壓力的變化范圍約為-8~10 kPa。超彈性球體的表面壓力變化范圍要比剛性球體大,這是由于局部擾流與球體變形作用導(dǎo)致的。
圖14 超彈性球體和剛性球體入水流場(chǎng)壓力分布云圖Fig.14 Pressure cloud diagram of the flow field for elastic sphere and rigid spheres
本文采用實(shí)驗(yàn)與數(shù)值模擬相結(jié)合的方法對(duì)超彈性球體低速垂直入水問(wèn)題進(jìn)行了研究,獲得的主要結(jié)論如下:
1)通過(guò)實(shí)驗(yàn)研究發(fā)現(xiàn)了超彈性球體入水后特有的“橫向扁平球體-縱向扁平球體-橫向扁平球體”的周期性變形行為,同時(shí)入水空泡的發(fā)展過(guò)程中會(huì)形成嵌入空泡。
2)由于變形產(chǎn)生的能量損耗,從而造成超彈性球體入水空泡長(zhǎng)度要小于剛性球體的入水空泡長(zhǎng)度。同時(shí),由于超彈性球體的變形行為,導(dǎo)致超彈性球體入水后形成的空泡寬度要大于剛性球體的空泡寬度。超彈性球體入水后,球體剪切模量越小,入水空泡口輪廓越大。
3)超彈性球體不同時(shí)刻的表面壓力分布與其變形情況密切相關(guān),由于局部擾流與球體變形導(dǎo)致超彈性球體表面壓力的變化范圍較大。