王 聰，夏維學，陳超倩

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

0 引言

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

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

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

1 數(shù)值分析方法

本文針對圓柱體在低弗勞德數(shù)條件下落水過程，忽略流體介質的可壓縮特性、空化和熱效應，假設圓柱體為剛體。考慮到具有初始傾角的圓柱體落水過程的復雜性和非定常多相流動特性，采用基于VOF多相流模型的大渦模擬方法，以及重疊網格方法開展圓柱體落水過程復雜多自由度運動及流場演化特性研究。

1.1 基本控制方程

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

αw=Vw/Vc

(1)

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

控制單元內混合介質密度ρm和混合動力黏度μm表達式為

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

(2)

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

(3)

式中，ρa為空氣相的密度，μa為空氣相的運動黏度，ρw為水相的密度，μw為水相的運動黏度。

(4)

(5)

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

(6)

(7)

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

(8)

1.2 計算模型

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

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

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

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

(a)xoy截面網格

(b)yoz截面網格

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

表1 圓柱體參數(shù)

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

圖3 圓柱體落水過程空泡演化實驗和數(shù)值結果對比Fig.3 Comparison of the experimental and numerical cavity for cylinder entering water

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

2 結果分析

2.1 落水空泡演化過程

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

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

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

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

(a)圓柱體落水過程流場演化

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

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

表2 不同直徑圓柱體落水空泡對比

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

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

2.2 閉合特征尺寸

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

(a)空泡口寬度變化

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

(a)浸沒深度

(b)水平偏移距離

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

2.3 水動力特性

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

(9)

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

2.3.1 升力系數(shù)

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

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

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

2.3.2 力矩系數(shù)

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

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

2.3.3 阻力系數(shù)

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

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

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

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

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

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

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

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

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

3 結論

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

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

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

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