半疏水-半親水球體垂直入水空泡數(shù)值仿真研究

孫釗， 曹偉， 王聰， 路中磊

(哈爾濱工業(yè)大學(xué) 航天學(xué)院, 黑龍江 哈爾濱 150001)

?

半疏水-半親水球體垂直入水空泡數(shù)值仿真研究

孫釗， 曹偉， 王聰， 路中磊

(哈爾濱工業(yè)大學(xué) 航天學(xué)院, 黑龍江 哈爾濱 150001)

采用流體體積多相流模型耦合連續(xù)表面力模型，對(duì)具有非對(duì)稱表面潤(rùn)濕性的半疏水- 半親水球體垂直入水空泡形態(tài)發(fā)展過(guò)程進(jìn)行數(shù)值仿真，分析了半疏水- 半親水球體垂直入水過(guò)程中受到的流體動(dòng)力。研究結(jié)果表明：半疏水-半親水球體垂直入水后，產(chǎn)生非對(duì)稱入水空泡及“心”型噴濺，同時(shí)球體運(yùn)動(dòng)軌跡將偏離原豎直運(yùn)動(dòng)軌道，由疏水半球一側(cè)向親水半球一側(cè)偏斜；入水初期，在球體表面形成液體薄層運(yùn)動(dòng)，在疏水半球一側(cè)，液體薄層與球體表面分離，導(dǎo)致空氣進(jìn)入形成敞開(kāi)空泡；在親水半球體一側(cè)，液體薄層沿球體表面向上運(yùn)動(dòng)最終在球體頂點(diǎn)匯聚；液體薄層匯聚后形成楔形流，楔形流在球體頂點(diǎn)與球體表面分離，繼續(xù)向疏水半球一側(cè)產(chǎn)生的入水空泡壁面運(yùn)動(dòng)并撞擊，形成“心”型噴濺。

流體力學(xué)； 多相流； 半疏水- 半親水球體； 垂直入水； 表面潤(rùn)濕性； 數(shù)值仿真； 空泡

0 引言

入水問(wèn)題是一個(gè)非常復(fù)雜的涉及跨介質(zhì)的固體及液體相互作用問(wèn)題[1-2]，在運(yùn)動(dòng)體入水過(guò)程中伴隨許多現(xiàn)象，例如撞擊瞬間、噴濺現(xiàn)象、入水空泡產(chǎn)生、空泡掐斷、空泡潰滅等。許多工程問(wèn)題及自然現(xiàn)象與入水問(wèn)題存在重要關(guān)系，例如跳彈[3]、空投魚雷[4]、水上飛機(jī)水上著陸[5]、水上行走生物[6]等．

近年來(lái)，隨著高速攝像以及高速粒子成像技術(shù)的發(fā)展，研究者們開(kāi)展了大量有關(guān)入水問(wèn)題的實(shí)驗(yàn)研究，例如Dulax等[7]、Grumstrup等[8]、Bergmann等[9]、Aristoff等[10]、Truscott等[11-12]、Techet等[13]、Sudo等[14]。Yilmaz等[15-16]通過(guò)粒子圖像測(cè)速(PIV)技術(shù)對(duì)球體垂直入水過(guò)程展開(kāi)了實(shí)驗(yàn)研究，獲得了球體入水過(guò)程空泡形態(tài)發(fā)展規(guī)律，并發(fā)現(xiàn)入水空泡閉合位置、閉合時(shí)間與弗勞德數(shù)Fr呈線性關(guān)系。Marston等[17]通過(guò)高速攝像技術(shù)拍攝了球體垂直入水過(guò)程在自由液面上方產(chǎn)生的噴濺，觀察到了噴濺輪廓的褶皺條紋，并提出了噴濺穹頂?shù)那什环€(wěn)定性。Benedict等[18]通過(guò)實(shí)驗(yàn)研究了剛性球體垂直進(jìn)入雙層液體的入水空泡現(xiàn)象，該雙層液體由不互溶的兩種液體組成。這種入水空泡，其空泡壁面不再保持光順，而是出現(xiàn)波紋狀漣漪，且空泡閉合時(shí)間、閉合深度與球體垂直落入均一純水中的入水空泡顯著不同。

入水物體的表面潤(rùn)濕性同樣對(duì)入水空泡有顯著影響。路中磊等[19]開(kāi)展了圓柱空心殼體的低速垂直入水過(guò)程，實(shí)驗(yàn)結(jié)果呈現(xiàn)了波動(dòng)狀態(tài)的入水空泡輪廓，入水速度以及空腔殼體幾何結(jié)構(gòu)對(duì)該波動(dòng)現(xiàn)象有較大影響。Duez等[20]研究了不同表面潤(rùn)濕性球體的垂直入水空泡，發(fā)現(xiàn)只有在入水速度超過(guò)某一臨界值的情況下才會(huì)產(chǎn)生入水空泡，該臨界值與球體的表面潤(rùn)濕性有重要關(guān)系。Do-Quang等[21]通過(guò)耦合Navier-Stokes方程以及Cahn-Hilliard方程方法，利用數(shù)值仿真方法研究了不同表面潤(rùn)濕性球體的垂直入水運(yùn)動(dòng)，獲得了大量的入水空泡矢量場(chǎng)圖片。Kintea等[22]通過(guò)數(shù)值仿真方法研究了旋轉(zhuǎn)剛性球體以及非旋轉(zhuǎn)剛性球體的入水過(guò)程，數(shù)值仿真方法建立在流體體積(VOF)多相流模型基礎(chǔ)上，引入表面張力以及接觸角，且數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好。

本文通過(guò)數(shù)值方法研究半疏水- 半親水球體的垂直入水空泡發(fā)展過(guò)程。所謂半疏水- 半親水球體，即球體兩側(cè)表面潤(rùn)濕性不對(duì)稱，一側(cè)半球表面為疏水性表面，另一側(cè)為親水性表面。一般情況下，球體以一定速度垂直撞擊自由液面后，將形成典型的“沙漏狀”入水空泡及“皇冠狀”噴濺。然而，當(dāng)球體的表面潤(rùn)濕性不對(duì)稱時(shí)，將產(chǎn)生較復(fù)雜的入水現(xiàn)象，例如非對(duì)稱入水空泡以及心型噴濺，這正是本文的研究重點(diǎn)。

1 數(shù)值方法

1.1 控制方程

本文所研究的垂直入水問(wèn)題，可簡(jiǎn)化為三維不可壓縮二相流問(wèn)題。數(shù)值計(jì)算采用VOF多相流模型描述計(jì)算流域中氣相及液相體積分?jǐn)?shù)，令液相體積分?jǐn)?shù)為α1，則α1=0代表該網(wǎng)格單元沒(méi)有液相分布；α1=1代表該網(wǎng)格單元全部為液相；0<α1<1代表相界面。對(duì)于二相流問(wèn)題，氣相的體積分?jǐn)?shù)即為1-α1.

流動(dòng)控制方程為

動(dòng)量守恒方程

(1)

質(zhì)量守恒方程

(2)式中：ρ為混合物密度；u為速度矢量；μ為混合物動(dòng)力黏度；p為壓強(qiáng)；g為重力加速度；F為表面張力項(xiàng)。

動(dòng)量方程中的混合物密度ρ以及混合物動(dòng)力黏度μ由計(jì)算流域中各相的體積分?jǐn)?shù)決定：

ρ=ρlαl+ρg(1-αl)，

(3)

μ=μlαl+μg(1-αl)，

(4)

式中：ρl、ρg、μl、μg分別為液相及氣相的密度和動(dòng)力黏度，具體數(shù)值如表1所示。

表1 數(shù)值過(guò)程液相及氣相屬性Tab.1 Properties of liquid and gas phases in simulation

本文采用Brackbill等[23]提出的連續(xù)表面力(CSF)模型，將表面張力項(xiàng)作為源項(xiàng)加入到動(dòng)量方程中:

(5)

式中：αi、αj分別為i相和j相的體積分?jǐn)?shù)；ρi、ρj分別為i相和j相的密度。

(6)

(7)

1.2 計(jì)算域及網(wǎng)格劃分

由于計(jì)算域及入水物體在幾何結(jié)構(gòu)上具有對(duì)稱性，為降低計(jì)算成本，本文采用1/2模型進(jìn)行數(shù)值計(jì)算。

計(jì)算域、邊界條件以及網(wǎng)格劃分如圖1所示。D為球體直徑，計(jì)算域軸向長(zhǎng)度40D，自由液面被設(shè)置在軸向長(zhǎng)度的正中間，即20D處。球體位置設(shè)定為剛好接觸自由液面前一瞬間的位置，并給予其初始速度u0，計(jì)算域整體尺寸為40D×20D×10D. 重力g方向與z軸正向一致，z軸通過(guò)球體球心，坐標(biāo)零點(diǎn)z=0定義為自由液面處，時(shí)間零點(diǎn)t=0定義為球體接觸自由液面瞬間時(shí)刻，數(shù)值計(jì)算過(guò)程以及后續(xù)處理選擇豎直向下為z軸正向。

圖1 計(jì)算域、邊界條件以及網(wǎng)格劃分Fig.1 Computational domain, boundary condition, and mesh generation

全計(jì)算域采用結(jié)構(gòu)化網(wǎng)格。球體周圍采用O-block網(wǎng)格對(duì)球體壁面進(jìn)行包裹，O-block尺度為4D×4D×2D，并對(duì)球體附近網(wǎng)格進(jìn)行加密，網(wǎng)格劃分策略采用BiGeometric，控制比率因子為1.1，如圖1(c)所示。

1.3 計(jì)算設(shè)置

本文應(yīng)用的數(shù)值計(jì)算軟件為ANSYS Fluent，邊界條件設(shè)置、計(jì)算流域初始化、自編程序的編譯等均在軟件界面中設(shè)置。數(shù)值計(jì)算過(guò)程中，采用VOF法對(duì)流動(dòng)控制方程離散；壓力場(chǎng)與速度場(chǎng)的耦合選用壓力隱式分裂算子(PISO)算法；壓力場(chǎng)的空間離散采用PRESTO！策略；各相體積率離散采用CICSAM策略；各求解變量的離散采用2階迎風(fēng)策略；對(duì)流項(xiàng)采用QUICK離散策略；相界面的幾何重構(gòu)采用Geo-Reconstruct策略。

2 結(jié)果與分析

2.1 數(shù)值與實(shí)驗(yàn)結(jié)果對(duì)比

本文通過(guò)與文獻(xiàn)[12]實(shí)驗(yàn)結(jié)果對(duì)比，驗(yàn)證數(shù)值方法的有效性。數(shù)值計(jì)算過(guò)程中入水運(yùn)動(dòng)參數(shù)設(shè)置與實(shí)驗(yàn)一致，球體直徑D=5.72 cm，球體質(zhì)量ms=0.17 kg，入水速度u0=1.72 m/s，親水性球體表面接觸角θ=60°，疏水性球體表面接觸角θ=120°.

首先，針對(duì)計(jì)算域幾何模型建立了不同網(wǎng)格密度的結(jié)構(gòu)化網(wǎng)格。通過(guò)逐漸降低網(wǎng)格單元尺寸，進(jìn)行網(wǎng)格細(xì)化研究：Finer mesh(包含1 538 340網(wǎng)格單元)；Fine mesh(包含987 600網(wǎng)格單元)；Coarse mesh(包含455 670 網(wǎng)格單元)。

圖2給出了球體垂直入水過(guò)程運(yùn)動(dòng)位移及運(yùn)動(dòng)速度隨時(shí)間變化。如圖2所示，隨網(wǎng)格數(shù)量增加，運(yùn)動(dòng)位移及速度差別逐漸減小。擬選定Finer mesh作為本文數(shù)值計(jì)算網(wǎng)格。在此情況下，數(shù)值計(jì)算在一個(gè)20核心處理器的服務(wù)器上進(jìn)行，時(shí)間步長(zhǎng)選取經(jīng)驗(yàn)性的設(shè)置1×10-6，平均每一個(gè)入水過(guò)程算例需要約48 h. 最終認(rèn)為Finer mesh可以保證計(jì)算精度的同時(shí)，所需的計(jì)算成本是可以接受的，因此作為本文數(shù)值計(jì)算的最終網(wǎng)格。

圖2 疏水性球體不同網(wǎng)格密度計(jì)算結(jié)果Fig.2 Numerical results of different mesh densities for hydrophobic sphere

兩種不同表面潤(rùn)濕性球體的垂直入水空泡數(shù)值仿真結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比分別示于圖3及圖4中。如圖3所示，在這種親水性球體的垂直入水過(guò)程中，球體撞擊自由液面后僅有一個(gè)垂直向上的濺射流產(chǎn)生。隨球體下落，在自由液面以下沒(méi)有形成入水空泡，球體完全被液體包裹，數(shù)值仿真結(jié)果與實(shí)驗(yàn)結(jié)果有較高的吻合度。圖4為另一組數(shù)值仿真結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比。在這種疏水性球體的入水現(xiàn)象中，入水初期液體在球體表面某一位置分離并形成入水空泡，隨后入水空泡在自由液面以下某一深度掐斷(深閉合)形成經(jīng)典的“沙漏狀”空泡形態(tài)，數(shù)值仿真結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)應(yīng)良好。

圖3 親水性球體入水空泡形態(tài)數(shù)值 仿真結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比Fig.3 Comparison of the numerical and experimental results of cavities created by hydrophilic sphere

圖4 疏水性球體入水空泡形態(tài)數(shù)值仿真 結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比Fig.4 Comparison of the numerical and experimental results of the cavities created by hydrophobic sphere

圖3及圖4所示的兩球體入水過(guò)程，其中球體直徑、入水速度以及液體屬性等外界參數(shù)完全一致，僅有球體表面潤(rùn)濕性的差別，卻產(chǎn)生了完全不同的入水空泡現(xiàn)象，而本文數(shù)值方法實(shí)現(xiàn)了這種因表面潤(rùn)濕性不同而產(chǎn)生的入水現(xiàn)象差異，且與實(shí)驗(yàn)結(jié)果對(duì)比良好。

入水空泡閉合現(xiàn)象是入水空泡發(fā)展過(guò)程的重要過(guò)程之一，本文分別提取了數(shù)值仿真結(jié)果與實(shí)驗(yàn)結(jié)果的入水空泡參數(shù)(見(jiàn)圖5)：閉合時(shí)間tp，閉合深度hp，閉合位置zp，閉合時(shí)刻的空泡開(kāi)口直徑Dp，具體量值如表2所示。

2.2 非對(duì)稱入水空泡形態(tài)

親水性球體與疏水性球體垂直入水后將產(chǎn)生差異較大的入水空泡形態(tài)(見(jiàn)圖3和圖4)，當(dāng)球體表面潤(rùn)濕性非對(duì)稱時(shí)，例如一側(cè)半球?yàn)橛H水性表面，另一側(cè)為疏水性表面，即半疏水- 半親水球體的垂直入水空泡形態(tài)如何？這正是本文的研究重點(diǎn)。

圖5 入水空泡閉合時(shí)刻坐標(biāo)示意Fig.5 Coordinate for cavity pinching-off表2 入水空泡閉合參數(shù)Tab.2 Pinch-off data of water entry cavity

參數(shù)數(shù)值仿真值實(shí)驗(yàn)值誤差/%tp/ms94597531hp/mm85587523zp/mm1716175120Dp/mm1298136348

圖6 不同表面潤(rùn)濕性球體垂直入水空泡形態(tài)Fig.6 Water entry cavities created by spheres with different wettabilities

圖6(a)為親水性球體的入水過(guò)程，球體表面接觸角θ=60°，俯視圖及側(cè)視圖均說(shuō)明該球體在入水過(guò)程中沒(méi)有入水空泡產(chǎn)生。

圖6(b)為疏水性球體垂直入水空泡發(fā)展過(guò)程，球體表面接觸角θ=120°，該球體在入水過(guò)程中產(chǎn)生了入水空泡，空泡在自由液面下方某一位置掐斷，形成典型“沙漏狀”入水空泡。

圖6(c)為半疏水- 半親水球體的垂直入水空泡發(fā)展過(guò)程，該球體左側(cè)半球表面為疏水性表面，θL=120°，右側(cè)半球表面為親水性表面，θR=60°. 可以預(yù)見(jiàn)，在該球體親水性半球一側(cè)，將呈現(xiàn)類似圖6(a)中的無(wú)入水空泡現(xiàn)象；而在該球體疏水性半球一側(cè)，將呈現(xiàn)類似圖6(b)的入水空泡；最終，將會(huì)形成非對(duì)稱入水空泡。圖6(c)示出了這種非對(duì)稱入水空泡的發(fā)展過(guò)程，并且從圖中可以觀察到該球體入水過(guò)程中不僅產(chǎn)生非對(duì)稱空泡，同時(shí)球體運(yùn)動(dòng)軌跡不再保持豎直向下，而是發(fā)生了橫向位移。該橫向位移方向由疏水側(cè)指向親水性一側(cè)，即從有空泡產(chǎn)生一側(cè)指向無(wú)空泡產(chǎn)生一側(cè)。

值得注意的是，這種半疏水- 半親水球體垂直入水后不僅產(chǎn)生非對(duì)稱入水空泡，同時(shí)形成“心”型噴濺狀態(tài)，如圖6(c)俯視圖所示。疏水性半球一側(cè)排開(kāi)液體形成敞開(kāi)入水空泡，親水性一側(cè)并不產(chǎn)生入水空泡，而是形成液體薄層。液體薄層在球體頂點(diǎn)匯聚后形成尖銳的楔形流，楔形流由親水半球一側(cè)向疏水半球一側(cè)運(yùn)動(dòng)，即圖6(c)中從右向左，最終撞擊到疏水半球一側(cè)產(chǎn)生的入水空泡壁面。隨球體下落，入水空泡及楔形流同時(shí)發(fā)展形成較明顯的“心”型噴濺。

圖7為楔形流產(chǎn)生初期入水空泡形態(tài)及相應(yīng)矢量場(chǎng)，空泡壁面在數(shù)值后處理過(guò)程中被設(shè)置為具有一定透明度，因此可以觀察到空泡壁面內(nèi)部的楔形流。圖7(a)和圖7(c)分別為空泡形態(tài)俯視圖及側(cè)視圖；圖7(b)和圖7(d)分別為距離自由液面0D、1/2D深度處橫斷面空泡輪廓及矢量場(chǎng)。如圖7所示，楔形流由親水性一側(cè)產(chǎn)生，在球體頂點(diǎn)位置附近與球體表面分離，并從親水半球一側(cè)向疏水半球一側(cè)空泡壁面運(yùn)動(dòng)，此時(shí)楔形流還沒(méi)有撞擊到左側(cè)空泡壁面；從俯視圖中可以觀察到，由于楔形流的產(chǎn)生，液面上方的噴濺不再呈圓周對(duì)稱，而是形成“心型”噴濺。

圖7 球體入水后t=22.75 ms入水空泡形態(tài)及矢量場(chǎng)Fig.7 Cavity formation and velocity field taken for t=22.75 ms after impact

這種楔形流的產(chǎn)生與親水性一側(cè)球體表面形成的液體薄層有重要關(guān)系，如圖8所示。對(duì)于親水性球體，產(chǎn)生的液體薄層沿球體表面向上運(yùn)動(dòng)，并最終在球體頂部匯合，液體薄層匯合后形成向上的濺射流，不產(chǎn)生入水空泡；對(duì)于疏水性球體，液體薄層在運(yùn)動(dòng)到球體赤道附近處與球體表面分離，從而導(dǎo)致空氣進(jìn)入并形成入水空泡。半疏水- 半親水球體由于球體表面的非對(duì)稱潤(rùn)濕性，液體薄層在兩側(cè)球體表面具有不同的運(yùn)動(dòng)狀態(tài)：在親水性一側(cè)，液體薄層沿表面向上運(yùn)動(dòng)，至球體頂點(diǎn)匯聚，形成楔形流尖端，并向左側(cè)繼續(xù)運(yùn)動(dòng)；在疏水性表面一側(cè)，液體薄層極早和球體表面分離，形成開(kāi)口入水空泡。

圖8 不同表面潤(rùn)濕性球體入水初期空泡形態(tài)Fig.8 Initial stage of water impact by spheres with different wettabilities

2.3 軌跡、速度、加速度

相對(duì)于表面潤(rùn)濕性均一的普通球體，具有非對(duì)稱表面潤(rùn)濕性的半疏水- 半親水球體垂直入水后，其最明顯以及最讓人感興趣的特征不僅是產(chǎn)生了非對(duì)稱入水空泡及“心”型噴濺，其運(yùn)動(dòng)軌跡也將發(fā)生偏斜，如圖9所示。

圖9 不同表面潤(rùn)濕性球體入水運(yùn)動(dòng)軌跡Fig.9 Trajectories of spheres during water entry

非對(duì)稱入水空泡的產(chǎn)生說(shuō)明球體兩側(cè)與周圍水域相互作用的動(dòng)量傳遞不對(duì)稱，導(dǎo)致球體產(chǎn)生水平方向位移(x方向)。球體垂直入水過(guò)程是一個(gè)球體運(yùn)動(dòng)動(dòng)能及動(dòng)量與周圍流體持續(xù)轉(zhuǎn)移變化的過(guò)程。隨球體下落，球體給每一層水以一個(gè)基本上是橫向的動(dòng)量。相應(yīng)地，流體將反作用力施加于球體本身。

在球體疏水性半球一側(cè)產(chǎn)生入水空泡，說(shuō)明疏水性半球一側(cè)固體- 液體之間的動(dòng)量轉(zhuǎn)移比較顯著，球體給予周圍流體橫向動(dòng)量，相應(yīng)地，流體傳遞給球體等值反向的沖量。而在親水性半球一側(cè)沒(méi)有入水空泡產(chǎn)生，球體和周圍流體幾乎不發(fā)生動(dòng)量傳遞。最終，導(dǎo)致球體兩側(cè)受到的流體動(dòng)力在水平方向上不平衡。因此，球體運(yùn)動(dòng)軌跡由疏水性一側(cè)向親水性一側(cè)偏移，并且入水空泡相應(yīng)發(fā)生一定彎曲，同時(shí)產(chǎn)生心型入水噴濺。

圖9示出了5種不同表面潤(rùn)濕性球體運(yùn)動(dòng)軌跡，球體相對(duì)密度ρ*=2.7，入水速度u0=2.37 m/s. 5種不同表面潤(rùn)濕性球體分別為：θ=30°；θL=60°，θR=30°；θL=120°，θR=30°；θL=150°，θR=30°；θ=150°. 球體球心位置的x和z坐標(biāo)均用球體直徑D進(jìn)行無(wú)量綱化，z/D=0代表自由液面。圖9中空心三角符號(hào)標(biāo)記了入水空泡閉合時(shí)刻球體球心位置。θ=30°和θL=60°，θR=30°兩種球體表面兩側(cè)均不產(chǎn)生入水空泡，因此不涉及空泡閉合。θ=30°及θ=150°兩種球體示出了一條筆直的運(yùn)動(dòng)軌跡，其運(yùn)動(dòng)方向?yàn)樨Q直向下，入水過(guò)程中球體兩側(cè)流體流動(dòng)狀態(tài)沒(méi)有發(fā)生不對(duì)稱現(xiàn)象，因此球體運(yùn)動(dòng)軌跡沒(méi)有偏離豎直方向。隨球體兩側(cè)表面潤(rùn)濕性差異增加，球體運(yùn)動(dòng)軌跡偏斜程度越加明顯。θL=60°,θR=30°球體運(yùn)動(dòng)軌跡僅輕微偏離了豎直方向。θL=120°，θR=30°；θL=150°，θR=30°兩種球體，均產(chǎn)生非對(duì)稱入水空泡，隨球體運(yùn)動(dòng)軌跡的偏斜，入水空泡也相應(yīng)發(fā)生一定彎曲，并且空泡閉合位置偏離原豎直運(yùn)動(dòng)軌道。球體運(yùn)動(dòng)軌跡由疏水性一側(cè)向親水性一側(cè)偏斜，即圖9中由左向右，空泡閉合位置也向親水性一側(cè)偏離原豎直軌道。

為更好地理解不同非對(duì)稱表面潤(rùn)濕性差異的半疏水- 半親水球體垂直入水空泡，圖10給出了3種不同半疏水- 半親水球體入水空泡發(fā)展過(guò)程，運(yùn)動(dòng)參數(shù)與圖9所述一致。3種球體分別為：θL=60°,θR=30°;θL=120°,θR=30°;θL=150°，θR=30°. 如圖10(a)所示，對(duì)于θL=60°,θR=30°球體，雖然球體兩側(cè)表面潤(rùn)濕性不對(duì)稱，但是垂直入水后球體兩側(cè)均不產(chǎn)生入水空泡，因此球體運(yùn)動(dòng)軌跡的偏斜并不明顯。隨半疏水-半親水球體兩側(cè)的非對(duì)稱表面潤(rùn)濕性差異增加，產(chǎn)生的非對(duì)稱空泡欲加明顯，從而導(dǎo)致球體運(yùn)動(dòng)軌跡的偏斜增大，如圖10(b)和圖10(c)所示。

圖10 不同半疏水- 半親水球體垂直入水空泡形態(tài)Fig.10 Water entry cavities created by three asymmetric spheres

2.4 流體動(dòng)力分析

半疏水- 半親水球體垂直入水后產(chǎn)生與普通球體完全迥異的入水空泡形態(tài)，尤其非對(duì)稱入水空泡的產(chǎn)生導(dǎo)致球體兩側(cè)受到的流體動(dòng)力不平衡，最終使其運(yùn)動(dòng)軌跡偏離豎直運(yùn)動(dòng)軌道。本小節(jié)討論非對(duì)稱潤(rùn)濕表面球體垂直入水過(guò)程中的流體動(dòng)力。首先，給出半疏水- 半親水球體入水過(guò)程位移、速度、加速度隨時(shí)間的變化，如圖11所示。

圖11 半疏水- 半親水球體垂直入水位移、速度以及加速度隨時(shí)間變化Fig.11 Position, velocity and acceleration in z and x directions as a function of non-dimensional time

圖11示出了5種不同表面潤(rùn)濕性球體的位移、速度、加速度隨時(shí)間的變化，球體密度、入水速度均與圖9相同。z/D、x/D分別代表豎直方向位移及水平方向位移；uz/u0、ux/u0分別為豎直方向速度及水平方向速度；az/g、ax/g分別代表豎直方向加速度及水平方向加速度。

如圖11(a)所示，θ=30°球體豎直方向位移衰減最小，表明其下落速度最快。隨球體整體疏水性增加，z/D衰減越加明顯。如圖11(b)所示，θ=30°以及θ=150°兩種球體由于表面潤(rùn)濕性均一，不產(chǎn)生非對(duì)稱入水空泡，沒(méi)有橫向位移。θL=60°，θR=30°球體橫向位移極小，表明其微弱的橫向偏移。θL=150°，θR=30°球體具有最明顯的橫向位移，說(shuō)明球體兩側(cè)表面潤(rùn)濕性非對(duì)稱差異越大，橫向位移越明顯。

圖11(c)和圖11(d)給出了與圖11(a)、圖11(b)對(duì)應(yīng)的速度變化。如圖11(c)所示，θ=30°球體豎直方向速度衰減相對(duì)平緩，隨球體表面疏水性增加，速度衰減逐漸增加，θ=150°球體速度uz/u0衰減最大。值得注意的是，擁有疏水性表面的球體，其豎直方向速度在入水空泡閉合期間u0t/D=4～6出現(xiàn)了拐點(diǎn)。通過(guò)前面對(duì)入水空泡形態(tài)發(fā)展過(guò)程的分析可知，速度出現(xiàn)拐點(diǎn)主要是由于入水空泡深閉合導(dǎo)致球體尾部空泡壓強(qiáng)突增，對(duì)球體產(chǎn)生一定沖擊作用。對(duì)于θ=30°及θL=60°，θR=30°兩種親水性表面球體，無(wú)論表面潤(rùn)濕性是否對(duì)稱，球體在下落過(guò)程中完全被周圍流體包裹，不產(chǎn)生入水空泡，因此不會(huì)出現(xiàn)拐點(diǎn)。這種突然加速的過(guò)程在圖11(e)和圖11(f)加速度隨時(shí)間變化曲線中可以進(jìn)一步觀察到。

軌跡、速度、加速度都包含了重力加速度的作用，為了分析球體在入水過(guò)程中受到的流體動(dòng)力，對(duì)球體進(jìn)行受力分析：球體在入水過(guò)程中受到的力為重力mg以及與周圍全部的流體(氣體及液體)相互作用產(chǎn)生的總流體動(dòng)力Fh，

(8)

所謂總流體動(dòng)力Fh，來(lái)源于與周圍流體相互作用所產(chǎn)生的壓力、黏性力、表面張力的綜合作用：

(9)

(10)

(11)

式中：u=u(t)為不同時(shí)刻的瞬時(shí)速度；ρ=ρw=998.2 kg/m3為液體密度；Fhz和Fhx分別為Fh在z和x方向上的分量。

圖12給出了不同表面潤(rùn)濕性球體垂直入水過(guò)程的總流體動(dòng)力系數(shù)隨時(shí)間變化。如圖12(a)所示，球體觸水瞬間，u0t/D=0～0.5時(shí)刻，經(jīng)歷幾乎相同的入水沖擊作用，說(shuō)明入水沖擊作用主要與球體的入水速度有關(guān)，和表面潤(rùn)濕性沒(méi)有顯著關(guān)系。隨后，因不同球體產(chǎn)生的入水空泡發(fā)展過(guò)程不同，形成一定差異的流體動(dòng)力。θL=30°，θR=30°和θL=60°，θR=30°兩種球體由于不產(chǎn)生入水空泡，在入水初期u0t/D=0.5～2，由于液體薄層沿表面的運(yùn)動(dòng)，其流體動(dòng)力形成一定波動(dòng)；隨后，球體一直處于被周圍液體完全包裹的狀態(tài)，其流體動(dòng)力趨于恒值。θL=120°，θR=30°和θL=150°，θR=30°兩種球體由于液體薄層與球體表面分離，形成入水空泡，在入水初期沒(méi)有發(fā)生類似的波動(dòng)，而是在空泡深閉合時(shí)刻(u0t/D=4～6)出現(xiàn)明顯峰值，說(shuō)明空泡閉合對(duì)球體有明顯沖擊作用。如圖12(b)所示，非對(duì)稱表面潤(rùn)濕性產(chǎn)生的水平方向流體動(dòng)力先增加后減小，并且在空泡深閉合發(fā)生時(shí)刻出現(xiàn)峰值，說(shuō)明非對(duì)稱空泡的深閉合對(duì)球體水平方向的運(yùn)動(dòng)同樣具有沖擊作用。對(duì)于非對(duì)稱差異較小的θL=60°，θR=30°球體，由于球體兩側(cè)均不產(chǎn)生空泡，這種水平方向的流體動(dòng)力微乎其微。

圖12 半疏水- 半親水球體垂直入水流體 動(dòng)力系數(shù)隨時(shí)間變化Fig.12 Hydrodynamic force coefficients as a function of non-dimensional time

3 結(jié)論

本文通過(guò)數(shù)值方法研究了半疏水- 半親水球體的垂直入水過(guò)程，對(duì)非對(duì)稱入水空泡形態(tài)以及入水過(guò)程中球體軌跡、速度、加速度、流體動(dòng)力進(jìn)行了分析。研究結(jié)果表明，通過(guò)改變球體表面潤(rùn)濕性，可以有效的改變球體垂直入水現(xiàn)象，尤其可以改變?nèi)胨张菪螒B(tài)。

具有非對(duì)稱潤(rùn)濕表面的半疏水- 半親水球體，垂直入水后產(chǎn)生非對(duì)稱入水空泡，說(shuō)明球體與周圍水域相互作用的動(dòng)量傳遞不對(duì)稱，最終導(dǎo)致由疏水半球一側(cè)向親水半球一側(cè)的橫向位移。通過(guò)不同半疏水- 半親水球體的垂直入水過(guò)程數(shù)值模擬，分析了球體位移、速度、加速度。研究結(jié)果表明，這種水平方向的位移和球體兩側(cè)的表面潤(rùn)濕性差異有重要關(guān)系，對(duì)于差異較小的θL=60°，θR=30°球體，橫向位移微乎其微；對(duì)于差異較大的θL=150°，θR=30°球體，橫向位移非常明顯，說(shuō)明非對(duì)稱潤(rùn)濕性差異越大，產(chǎn)生的非對(duì)稱入水空泡及軌跡偏轉(zhuǎn)越顯著。

半疏水- 半親水球體的垂直入水運(yùn)動(dòng)，展現(xiàn)了具有非對(duì)稱潤(rùn)濕性表面的運(yùn)動(dòng)體入水后更加復(fù)雜的入水空泡形態(tài)發(fā)展過(guò)程，進(jìn)一步說(shuō)明表面潤(rùn)濕性在入水過(guò)程中的重要作用。該方面研究可應(yīng)用于許多流體動(dòng)力問(wèn)題，尤其涉及需要考慮彈道、流體動(dòng)力、多相流場(chǎng)結(jié)構(gòu)等軍事領(lǐng)域及工程應(yīng)用。在入水運(yùn)動(dòng)過(guò)程中，考慮表面潤(rùn)濕性的研究目前較少，對(duì)于更寬廣入水條件的研究，例如更高的入水速度以及更復(fù)雜的幾何運(yùn)動(dòng)體，有待于進(jìn)一步探索。本文數(shù)值仿真方法與相應(yīng)實(shí)驗(yàn)結(jié)果吻合良好，在一定程度上奠定了用數(shù)值仿真方法研究考慮表面潤(rùn)濕性的入水多相流動(dòng)過(guò)程的可行性。

References)

[1] Worthington A M, Cole R S. Impact with a liquid surface, studied by the aid of instantaneous photography[J]. Philosophical Transactions of the Royal Society of London, 1900, 189:137-148.

[2] Worthington A M. A study of splashes[M]. NY, US: Longmans Green and Company, 1908.

[3] Rosellini L, Hersen F, Clanet C, et al. Skipping stones[J]. Journal of Fluid Mechanics , 2005, 543:137-146.

[4] May A. Water entry and the cavity-running behavior of missiles, ADA020429[R]. Silver Spring, MD, US: NAVSEA Hydroballistics Advisory Committee, 1975.

[5] Von Karman T. The impact on seaplane floats during landing, TN321[R]. Washington DC, US: National Advisory Committee on Aeronautics , 1929.

[6] Bush J W M, Hu D L. Walking on water: biolocomotion at the interface[J]. Annual Review of Fluid Mechanics, 2006, 38: 339-369.

[7] Duclaux V, Caillé F, Duez C,et al. Dynamics of transient cavities[J]. Journal of Fluid Mechanics, 2007, 591: 1-19.

[8] Grumstrup T, Keller J B, Belmonte A. Cavity ripples observed during the impact of solid objects into liquids[J]. Physical Review Letters, 2007, 99(11):114502-1-114502-4.

[9] Bergmann R, van der Meer D, Gekle S,et al. Controlled impact of a disk on a water surface: cavity dynamics[J]. Journal of Fluid Mechanics, 2009, 633:381-409.

[10] Aristoff J M, Truscott T T T, Techet A H, et al. The water entry of decelerating spheres[J]. Physics of Fluids, 2010, 22:032102-1-032102-8.

[11] Truscott T T, Techet A H. A spin on cavity formation during water entry of hydrophobic and hydrophilic spheres[J]. Physics of Fluids, 2009, 21: 121703-1-121703-4.

[12] Truscott T T, Techet A H. Water entry of spinning spheres[J]. Journal of Fluid Mechanics, 2009, 625:135-165.

[13] Techet A H, Truscott T T. Water entry of spinning hydrophobic and hydrophilic spheres[J]. Journal of Fluids and Structures, 2011, 27(5/6):716-726.

[14] Sudo S, Takayanagi H, Kamiyama S. Water entry of a magnetic fluid coated sphere[J]. Journal of Magnetism and Magnetic Materials, 2011, 323(10):1348-1353.

[15] Yilmaz N, Nelson R. Cavity dynamics of smooth sphere and golf ball at low Froude numbers, part 1: high-speed imaging and quantitative measurements[J]. Journal of Visualization, 2015, 18(2): 335-342.

[16] Yilmaz N, Nelson R. Cavity dynamics of smooth sphere and golf ball at low Froude numbers, part 2: PIV analysis[J]. Journal of Visualization, 2015, 18(4): 679-686.

[17] Marston J O, Truscott T T, Speirs N B,et al. Crown sealing and buckling instability during water entry of spheres[J]. Journal of Fluid Mechanics, 2016, 794: 506-529.

[18] Benedict C W T, Vlaskkamp J H A, Denissenko P, et al. Cavity formation in the wake of falling spheres submerging into a stratified two-layer system of immiscible liquids[J]. Journal of Fluid Mechanics, 2016, 790: 33-56.

[19] 路中磊，魏英杰，王聰，等. 基于高速攝像實(shí)驗(yàn)的開(kāi)放腔體圓柱殼入水空泡流動(dòng)研究[J]. 物理學(xué)報(bào)，2016,65(1):301-315. LU Zhong-lei, WEI Ying-jie, WANG Cong, et al. An experimental study of water-entry cavitating flows of an end-closed cylindrical shell based on the high-speed[J]. Acta Physica Sinica, 2016,65(1):301-315. (in Chinese)

[20] Duez C, Ybert C, Clanet C, et al. Making a splash with water repellency[J]. Nature Physics, 2007, 3:180-183.

[21] DoQuang M, Amberg G. The splash of a solid sphere impacting on a liquid surface: numerical simulation of the influence of wetting[J]. Physics of Fluids, 2009, 21:022102.

[22] Kintea D M, Breitenbach J, Gurumurthy V T, et al. On the influence of surface tension during the impact of particles on a liquid-gaseous interface[J]. Physics of Fluids, 2016, 28: 012108.

[23] Brackbill J U, Kothe D B, Zemach C. A continuum method for modeling surface tension[J].Journal of Computational Physics, 1992, 100:335-354.

Numerical Investigations on Water-entry Cavity of HalfHydrophobic-half Hydrophilic Sphere

SUN Zhao， CAO Wei， WANG Cong， LU Zhong-lei

(School of Astronautics, Harbin Institute of Technology, Harbin 150001, Heilongjiang, China)

The water entry of a solid sphere impacting on a liquid surface has challenged researchers for centuries and remains of interest to the researchers today. A simulation study of the water entry cavity of half hydrophobic-half hydrophilic sphere is performed. Particular attention is given to the simulation method based on solving the Navier-Stokes equations coupled with VOF model and CSF model. The numerical results are in agreement with the experimental results, thus validating the suitability of the numerical approach to simulate the water entry of sphere under different wetting conditions. Based on this method, the development of cavity created by the half hydrophilic-half hydrophobic sphere is investigated. Results show that the water entry of half hydrophobic-half hydrophilic sphere creates an asymmetric cavity and “cardioid” splash, which causes the sphere to travel laterally from the hydrophobic side to the hydrophilic side. Further investigations show that the fluid film presents during initial stage of impact, and on the half hydrophobic sphere, the fluid film detaches from the sphere to lead to cavity formation; on the half hydrophilic sphere, the fluid film moves up on the sphere surface and gathers at the vertex of the sphere, forming a wedge flow. The wedge flow moves and finally impacts on the opposite side of the cavity so as to cause “cardioid” splash. In addition, the total hydrodynamic force coefficient is investigated as a result of the forces acting on the sphere during water entry dictated by the cavity formation.

fluid mechanics; multiphase flow; half hydrophobic-half hydrophilic sphere; vertical water-entry; surface wettability; numerical simulation; cavity

2016-12-02

國(guó)家自然科學(xué)基金項(xiàng)目(11672094)； 哈爾濱市科技創(chuàng)新人才專項(xiàng)基金項(xiàng)目(2013RFLXJ007)

孫釗(1985—)， 男， 博士研究生。 E-mail: flame_1985@163.com

王聰(1966—)， 男， 教授， 博士生導(dǎo)師。 E-mail: alanwang@hit.edu.cn

TV131.2+2

A

1000-1093(2017)05-0968-10

10.3969/j.issn.1000-1093.2017.05.017