圓柱體落水過(guò)程多相流體動(dòng)力特性數(shù)值分析

王 聰，夏維學(xué)，陳超倩

(1. 哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱 150001；2. 海軍工程大學(xué)兵器工程系，武漢 430033；3. 西北機(jī)電工程研究所，咸陽(yáng) 712099)

0 引言

運(yùn)動(dòng)體落水廣泛存在于海洋工程、航空航天和日常生活，如船舶撞擊、魚(yú)雷入水、飛機(jī)水上迫降、滑跳等。本文采用數(shù)值方法開(kāi)展具有初始傾角的圓柱體落水過(guò)程多相流動(dòng)特性和多自由度運(yùn)動(dòng)特性研究。

針對(duì)運(yùn)動(dòng)體落水問(wèn)題的研究最早見(jiàn)于19世紀(jì)末，Worthington[1]針對(duì)不同參數(shù)的球體、液滴垂直下落撞擊自由液面落水的實(shí)驗(yàn)研究，分析了落水參數(shù)對(duì)空泡、水面波紋、回射流水柱、液滴反彈等現(xiàn)象的影響。運(yùn)動(dòng)體抨擊自由液面瞬間受到瞬時(shí)抨擊載荷作用[2-3]，該抨擊載荷可能導(dǎo)致運(yùn)動(dòng)體結(jié)構(gòu)破壞，內(nèi)部?jī)x器設(shè)備失效等問(wèn)題。為了解決落水撞擊問(wèn)題，May等[4-9]較為系統(tǒng)地研究了小球落水撞擊、空泡生成、空泡閉合，空泡潰滅等流動(dòng)現(xiàn)象，并分析了落水流動(dòng)機(jī)理及其水動(dòng)力特性[8-9]。Mcgehee等[10]開(kāi)展了返回艙原型和模型落水實(shí)驗(yàn)，定量分析了不同接觸角和攻角落水的加速度峰值。何春濤等[11]、路中磊等[12]針對(duì)圓柱體低速入水空泡形態(tài)特性進(jìn)行了實(shí)驗(yàn)研究，獲得了不同速度下空泡閉合與演化規(guī)律。Grumstrup等[13]、Yang等[14]使用不同材料球體，開(kāi)展了低弗勞德數(shù)條件下落水實(shí)驗(yàn)，解釋了運(yùn)動(dòng)體落水過(guò)程的空泡波紋現(xiàn)象。Bodily等[15]通過(guò)在模型內(nèi)嵌入內(nèi)測(cè)單元，獲得不對(duì)稱親疏水表面細(xì)長(zhǎng)體落水加速度和角加速度的時(shí)變特性。Truscott等[16-17]、夏維學(xué)等[18]進(jìn)行了大量的旋轉(zhuǎn)小球落水實(shí)驗(yàn)，得到了旋轉(zhuǎn)小球落水非對(duì)稱的空泡演化現(xiàn)象，并分析了低弗勞德數(shù)旋轉(zhuǎn)小球產(chǎn)生的水楔、馬格努斯效應(yīng)曲線軌跡等。Bocquet[19]、Rosellini等[20]采用實(shí)驗(yàn)和理論方法研究了剛性圓片落水彈跳現(xiàn)象，獲得了不同彈跳周期內(nèi)空泡演化、再落水抨擊時(shí)間的變化規(guī)律，并建立了預(yù)測(cè)圓片運(yùn)動(dòng)和能量耗散的理論模型。Belden等[21]、Farouk等[22]分別通過(guò)實(shí)驗(yàn)和數(shù)值模型開(kāi)展了小球落水彈跳研究，獲得了小球落水多相流場(chǎng)演化，給出了基于高弗勞德數(shù)到低弗勞德數(shù)的漸進(jìn)臨界角的修正經(jīng)驗(yàn)?zāi)Ｐ汀?/p>

本文通過(guò)具有初始傾角的圓柱體在低弗勞德數(shù)條件下落水多自由度運(yùn)動(dòng)的數(shù)值分析，開(kāi)展不同直徑圓柱體大角度落水過(guò)程空泡演化規(guī)律研究，對(duì)比分析尺寸效應(yīng)對(duì)落水空泡的影響，并獲得圓柱體落水過(guò)程流體動(dòng)力演化特性。

1 數(shù)值分析方法

本文針對(duì)圓柱體在低弗勞德數(shù)條件下落水過(guò)程，忽略流體介質(zhì)的可壓縮特性、空化和熱效應(yīng)，假設(shè)圓柱體為剛體?？紤]到具有初始傾角的圓柱體落水過(guò)程的復(fù)雜性和非定常多相流動(dòng)特性，采用基于VOF多相流模型的大渦模擬方法，以及重疊網(wǎng)格方法開(kāi)展圓柱體落水過(guò)程復(fù)雜多自由度運(yùn)動(dòng)及流場(chǎng)演化特性研究。

1.1 基本控制方程

VOF多相流模型是適用于多種互不相溶的流體介質(zhì)界面追蹤技術(shù)，并忽略相間的相對(duì)滑移。對(duì)于僅涉及水和氣兩相的流動(dòng)，定義αw為控制單元內(nèi)水相的體積分?jǐn)?shù)，其表達(dá)式為

αw=Vw/Vc

(1)

式中，Vw為控制單元內(nèi)水介質(zhì)的體積,Vc為控制單元的體積。根據(jù)以上定義， 1-αw為空氣相體積分?jǐn)?shù)，則給定控制單元內(nèi)流體體積分?jǐn)?shù)為：αw=1表示控制單元內(nèi)全為水介質(zhì),αw=0表示控制單元內(nèi)全為空氣,0<αw<1表示控制單元由空氣和水的混合介質(zhì)組成。

控制單元內(nèi)混合介質(zhì)密度ρm和混合動(dòng)力黏度μm表達(dá)式為

ρm=(1-αw)ρa(bǔ)+αwρw

(2)

μm=(1-αw)μa+αwμw

(3)

式中，ρa(bǔ)為空氣相的密度，μa為空氣相的運(yùn)動(dòng)黏度，ρw為水相的密度，μw為水相的運(yùn)動(dòng)黏度。

(4)

(5)

式中，xi為笛卡爾坐標(biāo)，下標(biāo)i=1,2,3分別代表笛卡爾坐標(biāo)系x,y,z這3個(gè)坐標(biāo)分量，下同；ui為笛卡爾坐標(biāo)系下的速度分量；t為落水物理時(shí)間；p為控制體正壓力；νm為混合運(yùn)動(dòng)黏度，定義為νm=μm/ρm；τij為亞格子尺度應(yīng)力，針對(duì)不可壓縮湍流流動(dòng)的亞格子尺度應(yīng)力張量模型為[24]

(6)

(7)

為了獲得精確的近壁面湍流，本文數(shù)值分析采用WALE SGS模型[25]，則νt的定義為

(8)

1.2 計(jì)算模型

本文基于Star-CCM+開(kāi)展圓柱體落水過(guò)程多相流動(dòng)特性研究。圖1為圓柱體大角度落水的計(jì)算域及邊界條件，該計(jì)算域與圓柱體落水試驗(yàn)的水箱一致[26]。計(jì)算域頂面的邊界條件為壓力入口，其余表面均為壁面邊界，如圖1所示。

圖1 值模擬計(jì)算域和邊界條件Fig.1 Computational domain and boundary condition

重疊網(wǎng)格技術(shù)被用于捕捉圓柱體落水后的復(fù)雜多自由度運(yùn)動(dòng)，圖2給出了采用重疊網(wǎng)格技術(shù)的網(wǎng)格分布，其中圓柱體直徑D=29 mm，長(zhǎng)度L=180 mm。為了保障計(jì)算精度的同時(shí)提高計(jì)算效率，采用切割體網(wǎng)格法生成如圖2所示的非結(jié)構(gòu)背景網(wǎng)格。為了準(zhǔn)確捕捉圓柱體落水后的水動(dòng)力變化和空泡瞬態(tài)演化，在背景域中建立了10D×3D×12D和25D×3D×24D兩個(gè)網(wǎng)格預(yù)加密區(qū)Ⅰ和Ⅱ。同時(shí)自由液面區(qū)域也進(jìn)行了加密處理，以提高自由液面變化捕捉精度。為了更好地控制圓柱體壁面邊界層網(wǎng)格的增長(zhǎng)以及圓柱體域的網(wǎng)格分布，采用ICEM軟件生成圓柱體域網(wǎng)格，如圖2(c)所示,保證了在本文研究的圓柱體落水速度條件下滿足y+<1的約束條件。采用直徑為3D、長(zhǎng)度為8D的圓柱計(jì)算域作為子計(jì)算域，數(shù)值分析過(guò)程將圓柱體表面邊界條件設(shè)置為無(wú)滑移壁面，并給定固定接觸角為67°。

為了研究尺寸效應(yīng)對(duì)圓柱體落水多相流動(dòng)的影響，以圓柱體直徑D=29 mm為參考，在保證圓柱體落水速度v0= 2.5 m/s、傾角α0=120°、長(zhǎng)細(xì)比、當(dāng)量密度以及質(zhì)量分布等初始參數(shù)相同的條件下，開(kāi)展不同直徑的圓柱體落水過(guò)程數(shù)值分析，表1中給出了圓柱體的詳細(xì)參數(shù)。

(a)xoy截面網(wǎng)格

(b)yoz截面網(wǎng)格

(c)圓柱體附近網(wǎng)格圖2 計(jì)算域網(wǎng)格分布Fig.2 Distribution of the computational grid

表1 圓柱體參數(shù)

1.3 數(shù)值分析方法驗(yàn)證

圖3 圓柱體落水過(guò)程空泡演化實(shí)驗(yàn)和數(shù)值結(jié)果對(duì)比Fig.3 Comparison of the experimental and numerical cavity for cylinder entering water

圖4 圓柱體落水加速度和角加速度變化時(shí)程數(shù)值和實(shí)驗(yàn)結(jié)果對(duì)比Fig.4 Comparison of cylinder acceleration and angular acceleration in water entry

2 結(jié)果分析

2.1 落水空泡演化過(guò)程

表2對(duì)比了不同直徑圓柱體落水不同時(shí)刻的空泡形態(tài)，為了凸顯直徑變化對(duì)落水空泡的影響，通過(guò)定長(zhǎng)高比縮放調(diào)節(jié)，保證圖像中所有圓柱體直徑相同。Truscott等[17,27]研究表明，運(yùn)動(dòng)體落水空泡閉合時(shí)間隨其特征尺度的變化而變化。為了對(duì)比相對(duì)閉合時(shí)刻不同空泡形態(tài)，表2中圖像提取時(shí)間間隔為0.2tp，同時(shí)給出了尾空泡閉合時(shí)刻tpt的圖像。

對(duì)于直徑D*≥0.75D的圓柱體，落水過(guò)程頭空泡在空泡敞開(kāi)、膨脹、分離、深閉合等演化過(guò)程高度一致，但是自由液面附近空泡口寬度(dcw)與圓柱體直徑之比(dcw/D*)隨直徑增加而逐漸減小，而尾空泡的演化形式差異明顯。首先，當(dāng)水楔侵入尾空泡后，空氣填充水楔尾部形成明顯的尾空泡凸起現(xiàn)象，且圓柱體直徑越小，凸起越明顯(如表2中tp時(shí)刻)。其次，當(dāng)D*=2D的圓柱體落水尾空泡直徑較大，因此在當(dāng)前落水速度下，沿圓柱體背流面的相對(duì)流動(dòng)水射流沒(méi)有切斷尾空泡。圖5(a)為圓柱體落水過(guò)程流場(chǎng)演化。圖5(b)給出了表2中A-A截面空泡輪廓在tpt時(shí)刻的矢量分布。從圖5中可以看到，尾空泡在水楔沖擊作用下分裂為兩個(gè)渦旋方向相反的空泡，最后對(duì)流旋轉(zhuǎn)導(dǎo)致尾空泡閉合，而直徑D*=1.5D,1D和0.75D的尾空泡均被沿圓柱體的水射流切斷而形成閉合空泡。

不同直徑圓柱體在落水相同時(shí)刻的速度高度一致，即直徑越大的圓柱體壁面曲率越小，越不容易發(fā)生流動(dòng)分離，反之亦然。因此，D*≥0.75D的圓柱體落水后均發(fā)生空泡從圓柱體表面分離的現(xiàn)象。然而對(duì)于直徑D*=0.5D和0.25D的圓柱體落水空泡演化呈現(xiàn)明顯的尺寸效應(yīng)。圓柱體側(cè)空泡形成后與頭空泡連在一起形成連接空泡層。當(dāng)圓柱體尾部浸沒(méi)在自由液面下時(shí)，空泡層發(fā)生拉脫現(xiàn)象。隨后，對(duì)于直徑D*=0.5D的圓柱體，空泡層從拉脫位置自上而下破裂，直至空泡閉合潰滅。側(cè)空泡在空泡層破裂后沿圓柱體壁面向尾部運(yùn)動(dòng)，空氣不斷逃逸至尾空泡中，導(dǎo)致尾空泡在閉合后仍形成較大的駐留空泡。

針對(duì)D*=0.25D的圓柱體，空泡層從頭空泡閉合點(diǎn)附近自下而上破裂?？张萜屏押罂焖傧蛭膊亢拖蛏狭鲃?dòng)，破裂射流攜眷空泡沖擊尾空泡使其閉合?？张輰釉陬^空泡和尾空泡閉合后并未完全分裂(如表2中tpt時(shí)刻所示)?？諝庠诳张輰涌焖偈湛s過(guò)程需通過(guò)連接頭空泡的空泡層排出(如圖5(c)所示為1.12tp時(shí)刻空泡內(nèi)氣體流動(dòng))，所以空泡層沒(méi)有完全分裂。

(a)圓柱體落水過(guò)程流場(chǎng)演化

(b)A-A截面空泡輪廓速度矢量

(c)xoy平面空泡內(nèi)速度矢量圖5 局部速度矢量分布Fig.5 Local distribution of velocity vector

表2 不同直徑圓柱體落水空泡對(duì)比

圖6為無(wú)量綱頭空泡閉合時(shí)間(v0tp/D)隨弗勞德數(shù)Fr的變化。由表2空泡閉合形態(tài)可以看出，不同直徑圓柱體落水頭空泡均為拉斷深閉合。即使是直徑D*=0.25D的圓柱體，落水頭空泡也是在環(huán)境壓力作用下自中部附近向內(nèi)收縮，然后出現(xiàn)空泡層破裂，并最終發(fā)生拉斷閉合，即v0tp/D隨Fr增加而線性增大。尾空泡的閉合形式比較復(fù)雜，很難僅根據(jù)尾空泡的演化結(jié)構(gòu)去判斷空泡是否閉合。

圖6 空泡閉合時(shí)間隨Fr的變化Fig.6 Change of pinch-off time as a function of froude number Fr

2.2 閉合特征尺寸

圖7為不同直徑圓柱體落水空泡閉合瞬間特征尺寸隨Fr的變化。其中,圖7(a)為頭空泡閉合瞬間自由液面附近無(wú)量綱空泡口寬度dcw/D*隨Fr的變化，圖中橫坐標(biāo)以10為底的對(duì)數(shù)變化。從圖7(a)中可以看出，dcw/D*隨Fr的對(duì)數(shù)增加線性增長(zhǎng)。頭空泡閉合瞬間無(wú)量綱閉合點(diǎn)深度(hp/D*)以及圓柱體底面中心下降深(hb/D*)度隨Fr0.5增加而呈線性增長(zhǎng)，如圖7(b)所示。同時(shí)，圖7(b)中給出了hp/hb隨Fr的變化，在本文研究范圍內(nèi)，hp/hb在Fr≤5.62時(shí)隨Fr增加呈線性增長(zhǎng)，而在Fr≥5.6范圍內(nèi)保持定值hp/hb≈0.48。

(a)空泡口寬度變化

(b)空泡閉合深度和圓柱體底面深度變化圖7 不同直徑圓柱體落水空泡閉合瞬間特征尺寸Fig.7 Cavity critical sizes at pinch-off moment for the cylinders with different diameters

(a)浸沒(méi)深度

(b)水平偏移距離

(c)軌跡長(zhǎng)度圖8 空泡閉合瞬間圓柱體質(zhì)心位置隨Fr的變化Fig.8 Changes in cylinder center as a function of Froude number at pinch-off moment

2.3 水動(dòng)力特性

本文采用瞬時(shí)動(dòng)壓力對(duì)圓柱體阻力、升力和力矩進(jìn)行無(wú)量綱處理，得到圓柱體落水后的阻力系數(shù)Cd、升力系數(shù)Cl和力矩系數(shù)CM，其表達(dá)式為

(9)

式中，A= πD2/4為圓柱體橫截面積。

2.3.1 升力系數(shù)

圓柱體落水升力系數(shù)Cl變化特性如圖9所示。對(duì)于具有初始傾角的圓柱體撞擊自由液面瞬間，底面下緣壓強(qiáng)取得極大值[28]，如圖9(a)所示。隨后圓柱體側(cè)壁面撞擊水介質(zhì)，Cl逐漸增大，直至頭空泡發(fā)生閉合。空泡閉合瞬間在閉合點(diǎn)附近形成局部高壓，由表2中tp時(shí)刻空泡形態(tài)可知，閉合壓力導(dǎo)致圓柱體背流側(cè)壓力增加，使得壓差力的方向與升力方向相反，所以Cl在空泡閉合時(shí)出現(xiàn)負(fù)向脈沖?？张蓍]合后，圓柱體旋轉(zhuǎn)誘導(dǎo)升力起支配作用且隨圓柱體下降逐漸減小，但是圓柱體瞬時(shí)速度減小更快，因此Cl呈現(xiàn)緩慢增加的趨勢(shì)。通過(guò)Cl隨落水時(shí)間的變化發(fā)現(xiàn)，除去抨擊和空泡閉合階段外，Cl大致呈現(xiàn)線性變化，對(duì)空泡閉合前后的Cl進(jìn)行線性擬合，獲得如圖9(b)所示的變化特性。

(a)升力系數(shù)隨落水時(shí)間的變化

(b)空泡閉合前后升力系數(shù)漸進(jìn)特性圖9 不同直徑圓柱體落水升力系數(shù)變化特性Fig.9 Properties of lift coefficients for the cylinders with different diameters

2.3.2 力矩系數(shù)

(a)力矩系數(shù)隨落水時(shí)間的變化

(b)空泡閉合前后力矩系數(shù)漸進(jìn)特性圖10 不同直徑圓柱體落水力矩系數(shù)變化特性Fig.10 Properties of lift coefficients for the cylinders with different diameters

2.3.3 阻力系數(shù)

圖11(a)為不同直徑圓柱體落水阻力系數(shù)Cd隨落水無(wú)量綱時(shí)間t/tp的變化。圓柱體落水后Cd迅速增加，并在圓柱體底面中心沒(méi)入自由液面后達(dá)到極大值。當(dāng)落水速度一定時(shí)，到達(dá)極大值所需時(shí)間隨D*增加而逐漸增加，且極大值隨圓柱體直徑增加也緩慢增加，如圖11(a)中局部圖ip所示。落水撞擊階段結(jié)束受后，Cd隨圓柱體下降而逐漸增加。由圖11(a)可見(jiàn)，不同直徑圓柱體的Cd從撞擊至頭空泡閉合期間變化一致，但是在空泡閉合后出現(xiàn)明顯分化，直徑越小的圓柱體Cd在空泡閉合后增加越快。

空泡閉合前后Cd的平均增長(zhǎng)率隨Fr的變化如圖11(b)所示。不同直徑圓柱體落水的對(duì)數(shù)在空泡閉合前隨Fr的對(duì)數(shù)增加而線性減小，分析發(fā)現(xiàn)隨以自然數(shù)e為底Fr的指數(shù)變化的倒數(shù)1/eFr成線性關(guān)系。

通過(guò)上述分析可知，帶傾角圓柱體垂直落水尾空泡閉合形式非常復(fù)雜，且很難僅根據(jù)尾空泡的演化形態(tài)定量測(cè)量尾空泡的閉合時(shí)間。但是，根據(jù)上述分析可知，尾空泡閉合形成的局部高壓將導(dǎo)致圓柱體阻力減小，因此再結(jié)合Cd的變化，可以比較準(zhǔn)確地判斷尾空泡閉合時(shí)間tpt。

(a)阻力系數(shù)隨落水時(shí)間的變化

(b)空泡閉合前阻力系數(shù)漸進(jìn)特性圖11 不同直徑圓柱體落水力系數(shù)變化特性Fig.11 Properties of force coefficients for the cylinders with different diameters

從圖11(a)中可以看到，不同直徑圓柱體落水Cd在空泡閉合后存在兩次較為明顯的衰減，在圖中分別用pt和ps表示。沿圓柱體表面的水射流切斷尾空泡形成局部壓力增加是導(dǎo)致D*=2D,1.5D,1D和0.75D的圓柱體阻力在pt時(shí)刻衰減的原因；而D*=0.25D和0.5D的圓柱體則是由頭空泡閉合射流水楔沖擊圓柱體壁面導(dǎo)致阻力衰減，因此從圖中可以觀察到Cd衰減明顯(第一次衰減)。直徑D*=2D,1.5D,1D和0.75D的圓柱體尾空泡被水射流切斷后的上半尾空泡中間部分收縮、旋轉(zhuǎn)最終發(fā)生拉斷閉合，閉合形成的局部高壓導(dǎo)致Cd出現(xiàn)明顯的衰減現(xiàn)象；而D*=0.25D和0.5D的圓柱體落水，閉合射流不斷抨擊圓柱體壁面而釋放了能量，所以Cd變化不明顯，但是尾部空泡在閉合瞬間具有非常明顯特征。因此，空泡形態(tài)和閉合水動(dòng)力特征可以比較準(zhǔn)確地判斷尾空泡閉合。

圖12為不同直徑圓柱體落水空泡閉合后阻力系數(shù)第一次衰減及尾空泡閉合的無(wú)量綱時(shí)間隨弗勞德數(shù)的平方Fr2的變化規(guī)律。第一次衰減的無(wú)量綱時(shí)間tptptv0/D的對(duì)數(shù)隨Fr2的增加線性增長(zhǎng)，如圖12(a)所示。尾空泡無(wú)量綱閉合時(shí)間tptpsv0/D的對(duì)數(shù)隨Fr2的對(duì)數(shù)呈線性增加，如圖12(b)所示。

(a)第一次阻力系數(shù)衰減閉合時(shí)間

(b)尾空形閉合時(shí)間圖12 不同直徑圓柱體落水泡閉合后階段特征時(shí)間隨Fr2的變化Fig.12 The characteristic times after pinch off for water entry cylinders as a function of Fr2

3 結(jié)論

本文針對(duì)圓柱體大角度落水過(guò)程多自由度運(yùn)動(dòng)特性進(jìn)行數(shù)值計(jì)算，對(duì)比分析了不同直徑圓柱體落水空泡演化及水動(dòng)力特性，主要結(jié)論如下：

1)數(shù)值分析獲得的圓柱體落水空泡形態(tài)及受到的水動(dòng)力與實(shí)驗(yàn)結(jié)果均高度一致，驗(yàn)證了本文所建立的數(shù)值分析方法的有效性。

2)直徑大于臨界值的圓柱體的落水頭空泡敞開(kāi)、膨脹、分離、深閉合等演化過(guò)程相同，但自由液面附近無(wú)量綱空泡口寬度隨圓柱體直徑增加而逐漸減小。直徑小于臨界值的圓柱體落水空泡演化呈現(xiàn)明顯側(cè)空泡現(xiàn)象。不同直徑圓柱體落水頭空泡均呈現(xiàn)深閉合，且空泡閉合無(wú)量綱時(shí)間隨弗勞德數(shù)增加而線性增大。

3)除落水抨擊和空泡閉合階段外，空泡閉合前升力系數(shù)的平均增長(zhǎng)率總體比空泡閉合后快。不同直徑圓柱體的阻力系數(shù)從落水撞擊至空泡閉合階段變化一致，但是在空泡閉合后出現(xiàn)分化現(xiàn)象，直徑越小的圓柱體阻力系數(shù)在空泡閉合后增加越快。