王永虎， 林天龍

(1. 重慶交通大學(xué) 航空學(xué)院，重慶 400074; 2. 中國民航飛行學(xué)院 飛行技術(shù)學(xué)院，四川 廣漢618307)

0 引 言

目前結(jié)構(gòu)物在水中跌落問題的研究涉及到許多領(lǐng)域，如水下航行器入水、失事船舶或航空器碎片水中跌落、深彈水下彈道等。結(jié)構(gòu)物在水中跌落是連續(xù)的，水中的環(huán)境比較復(fù)雜，其軌跡很難準(zhǔn)確預(yù)測。在2014年3月發(fā)生的馬航MH370失蹤事件中，為了搜尋墜毀的殘骸，馬來西亞政府及世界各國耗費相當(dāng)大的財力和資源。然而海洋環(huán)境相當(dāng)復(fù)雜并且瞬息萬變，直至2017年1月，馬來西亞政府不得已宣布對失蹤的MH370客機進行搜查。此次失敗的搜查結(jié)果也印證了結(jié)構(gòu)物在跌入洋流之后，其真實的運動軌跡是難以人為預(yù)測的。物體跌落到水中會導(dǎo)致復(fù)雜的流體現(xiàn)象[1-2]，擾動也會對物體造成沖擊載荷，從而改變了物體的運動軌跡。針對物體在水中跌落現(xiàn)象的研究有很長的歷史[3]。K. T. HOLLAND[4]等在水中運動實驗中觀察到，當(dāng)參數(shù)不一的圓柱體在體心、重心、方向以及下降路徑等情況相同時，向重心的偏移致使物體發(fā)生豎直方向的變化等諸多復(fù)雜的實驗現(xiàn)象。許多學(xué)者致力于圓柱體自由跌落到水中的現(xiàn)象研究，未考慮到水流的影響[5-6]。A. V. ABELEV等[7]做了大量的關(guān)于全比圓柱體的實驗，同時對流體力學(xué)的特性進行了進一步的研究和爭論。J. MANN等[8]建立計算模型來預(yù)測圓柱狀的魚雷機體撞擊水表面和跌落到海底的水中運動情況。M. HOROBITZ等[9]對結(jié)構(gòu)物(主要為圓柱體)在只受重力狀況下，其入、出水方面動力學(xué)進行了更加深入的探究。陳宇翔[10]等用VOF方法對圓柱體入水進行了二維模擬仿真并與實驗對比，討論流體粘性的影響。屈秋林[11]等采用整體動網(wǎng)格法和VOF方法模擬分析飛機水上迫降情況。總之，針對物體在靜止水中跌落的研究較多，而對運動水流情況下的結(jié)構(gòu)物水中跌落問題研究較少。

主要通過ANSYS FLUENT軟件對圓柱體在有水流的水中跌落進行三維模擬仿真，采用了可實現(xiàn)的模型的湍流模型，運用了動網(wǎng)格和六自由度模型對圓柱體與水面相對運動進行有效處理；同時采用VOF對自由液面附近的變化進行了較好的捕捉，對水中的運動進行了較好的模擬；通過對整個跌落過程進行模擬，重點針對水流跌落軌跡進行了詳細分析。

1 問題描述

水流中跌落試驗裝置示意圖如圖1。

試驗對象為圓柱體，具體參數(shù)如表1。實驗水池的長寬高為2.0 m×0.60 m×0.60 m，水流速度為0.1 m/s，結(jié)構(gòu)物入水時的初始速度設(shè)置為0.5 m/s。采用電磁鐵控制夾具是否分離的方式，使圓柱體在12.739 mm的高度下以自由落體方式，落進實驗水池內(nèi)。利用高速攝像機，準(zhǔn)確清晰地將圓柱體落水時的運動軌跡和現(xiàn)象完整的記錄下來。

圖1 跌落實驗示意Fig. 1 Sketch of falling experiment

參數(shù)名稱數(shù)值參數(shù)名稱數(shù)值質(zhì)量/kg0.084 67Ixx/(kg·m2)2.843 5×10-4半徑/m0.01Iyy/(kg·m2)4.233 4×10-6高度/m 0.1Izz/(kg·m2)2.843 5×10-4

采用ICEM CFD16.0建立長寬高為0.2 m×0.16 m×0.25 m的三維數(shù)值實驗水池。非結(jié)構(gòu)四面體網(wǎng)格作為數(shù)值模型的整體網(wǎng)格的基礎(chǔ)網(wǎng)格，結(jié)構(gòu)物入水跌落軌跡作為實驗的側(cè)重點。為了提高實驗效率，對圓柱體周圍和氣液交界面用進一步加密的密度盒網(wǎng)格，而在結(jié)構(gòu)物撞擊區(qū)之外的區(qū)域(也就是流固耦合區(qū)域)，采用粗網(wǎng)格的方式進行實驗。為了模擬使得水流速度恒定，將流速為0.1 m/s的速度入口作為水流入口的邊界條件，將壓力出口作為出口水流的邊界條件，數(shù)值水池的上下面的入口邊界條件設(shè)為速度值為0，其余邊界條件均設(shè)為對稱邊界條件。

2 軟件介紹、理論基礎(chǔ)和計算方法

結(jié)構(gòu)物水中跌落涉及流固耦合和多相流等方面問題。結(jié)構(gòu)物從空氣中進入水中并在水中運動會改變流場、速度和壓強場等，并且結(jié)構(gòu)物運動會隨水流的運動而改變。圓柱體水中跌落現(xiàn)象復(fù)雜，不僅僅受到重力和浮力，還受到湍流阻力和水流推力等的影響。使用VOF方法[12]對自由表面進行捕捉，且修改邊界條件，打開明渠模型，用來制造水流。采用了可實現(xiàn)的k-ε模型的湍流模型，運用了有限體積法求解方程。有限體積法是CFD中比較成熟的一種運算方法。筆者將三維圓柱體視為剛體進行運算。三維圓柱體自由落體跌落到水中，并在水中下落是一個復(fù)雜的過程。筆者采用了六自由度動網(wǎng)格技術(shù)對圓柱體運行進行分析，并對其運動位置進行記錄。

2.1 FLUENT軟件介紹

ANSYS FLUENT是一款模擬流場、熱傳遞和化學(xué)反應(yīng)比較成熟的CFD軟件。在采用非結(jié)構(gòu)網(wǎng)格處理復(fù)雜模型時，F(xiàn)LUENT用其完整的網(wǎng)格適應(yīng)性來解決采用非結(jié)構(gòu)網(wǎng)格的流體問題。FLUENT支持的網(wǎng)格類型主要為二維網(wǎng)格(三角形與四邊形網(wǎng)格)、三維網(wǎng)格(四面體、六面體、棱柱等多面體)和混合網(wǎng)格，并且能夠根據(jù)具體情況細化或加粗網(wǎng)格。

2.2 控制方程

拉格朗日有限體積分析的質(zhì)量守恒方程：

(1)

式中：ρ為密度；t為時間；u、v和w為物體x、y和z方向的速度。

動量守恒方程為：

(2)

(3)

因為入水和水中跌落問題幾乎不考慮熱量交換等熱力學(xué)問題，所以只采用質(zhì)量守恒方程和動量守恒方程。

2.3 VOF模型

用VOF模型來模擬多相流的方式是先求解一組動量方程，然后在區(qū)域內(nèi)追蹤各種流體的體積分數(shù)。VOF方程主要適用于兩種或兩種以上不相容流體。另外一種流體在添加到模型之后，其體積分數(shù)也將進入到計算網(wǎng)格中。在任何一個網(wǎng)格中的參數(shù)和屬性不僅是單一相的代表，也是混合相的代表，其大小主要依賴體積分數(shù)。

假設(shè)第q種流體體積分數(shù)在網(wǎng)格下為αq，αq=0表示網(wǎng)格的第q種流體為空；αq=1表示第q種流體充滿網(wǎng)格的條件下，主要包括第q種流體和其他流體。體積分數(shù)的離散公式：

(4)

式中：n+1為新時間步長的指數(shù)；n為前一時間步長的指數(shù)；αq,f為第q個相流體的體積分數(shù)面數(shù)值；V為網(wǎng)格體積；Uf為以一定速度流經(jīng)面網(wǎng)格的體積通量。

2.4 明渠模型

水的流速設(shè)置通常情形下用明渠模型[13]。實驗中為了防止自由液面崩潰現(xiàn)象的發(fā)生，在壓力出口處調(diào)整了自由液面和底部的高度，用壓力規(guī)范法、密度插值法來處理自由液面。

自由液面ylocal在控制方程中為：

(5)

試驗仿真表明，自由液面的邊界變化在使用明渠模型后并未出現(xiàn)崩潰的情形，結(jié)果如圖2。

圖2 自由液面變化Fig. 2 Changing diagram of free surface level

2.5 動網(wǎng)格的使用

一般采用動網(wǎng)格來把握流場的運動和追蹤邊界隨時間的變化動態(tài)[14]。網(wǎng)格方法選用彈簧近似光順模型和局部重構(gòu)模型，其中彈性常數(shù)設(shè)置為0.5，網(wǎng)格重畫方法為局部網(wǎng)格和局部面，打開尺寸函數(shù)，其他系數(shù)使用默認值，并且打開六自由度模型以便獲取圓柱體的重心位移。使用UDF定義圓柱體的質(zhì)量和轉(zhuǎn)動慣量。

慣性坐標(biāo)系中，重心移動的運動控制方程為：

(6)

2.6 湍流模型

使用可實現(xiàn)的k-ε模型的湍流模型，可實現(xiàn)的k-ε模型包含湍流黏度的可選擇公式，在改良的運輸方程中，湍流耗散率ε來自于有關(guān)平方渦度波動運輸?shù)木_方程。

其運輸方程如式(7)、式(8)：

ρε-YM+Sk

(7)

(8)

式中：

(9)

(10)

(11)

式中：Gk為來自于平均速度的梯度湍流動能；Gb為浮力產(chǎn)生的湍流動能；YM為在可壓縮流波動的膨脹值；C2與C1ε為常數(shù)；σk和σε是k與ε的湍流普朗特數(shù)。

2.7 SIMPLEC算法

(12)

(13)

(14)

式中：

(15)

(16)

式中：p′為網(wǎng)格壓力修正值。

將通量修正方程(14)、(15)帶入方程(12)得到：

(17)

式中：

(18)

2.8 二階迎風(fēng)格式

采用二階迎風(fēng)格式，其具有較小的擴散性，方程為：

(19)

式中：

(20)

(21)

a+=max(a,0)

(22)

a-=min(a,0)

(23)

為了判斷式(19)是否穩(wěn)定，需滿足：

(24)

2.9 計算方法

選用SIMPLC算法求解壓力和速度耦合、最小二乘法求解梯度方程、PRESTO!求解壓強差值、體積重構(gòu)求解體積分數(shù)，其他均為二階迎風(fēng)求解。

3 計算結(jié)果及分析

對三維圓柱體在水流中跌落的相關(guān)問題，以利用動網(wǎng)格技術(shù)與VOF方法相結(jié)合方式來進行數(shù)值模擬試驗，先分析整體階段，進而深入分析圓柱體水中跌落階段，將實驗與模擬的成果相互對比，最后得出結(jié)論。

3.1 整體階段的速度分析

空中圓柱體只在重力的作用下跌落水中，然后在水中持續(xù)運動。在整個過程中圓柱體先受到重力的作用，入水之后逐漸轉(zhuǎn)變?yōu)橹亓Α⒏×σ约八鳟a(chǎn)生的水動力等復(fù)雜的合力，單一力到復(fù)雜力的演變過程，致使其各方向的運動速度也隨其情況的變動產(chǎn)生變動。

假設(shè)將試驗物體在剛與水接觸時設(shè)為0 s，此時速度為0.5 m/s。圖3顯示了結(jié)構(gòu)物從入水后運動過程的豎直方向速度隨跌落水中時間變化而發(fā)生的變化?？梢杂^察到物體由于受到流體的浮力以及湍流阻力的影響，在與水面接觸后，豎直方向速度變化由平緩逐漸遞增；在試驗物體完全入水中后，豎直方向的速度仍在增加，隨后豎直方向加速度突然增加，然后隨時間變化逐漸變得平緩。故物體在整個過程中受力逐漸趨于平穩(wěn)，且其在水中豎直方向速度的變化對物體水中跌落的軌跡產(chǎn)生重大的影響。

圖3 豎直方向速度的時間歷程Fig. 3 Vertical velocity time history diagram

圖4顯示的變化曲線為水平方向速度隨時間的變化狀況。由于水流的作用，圓柱體剛開始入水時與水的接觸面略小，水流作用小，其速度也就越平緩。隨著時間推移，在圓柱體全部浸入水中后，逐漸增加的水流作用力使得速度增加變快。當(dāng)圓柱體水平方向速度約70 mm/s時，水平速度增加變得緩慢，結(jié)構(gòu)物水平方向的位移也在不斷地變大。

圖4 水平方向速度的時間歷程Fig. 4 Horizontal velocity time history diagram

因為結(jié)構(gòu)物密度大于流體的密度，在水流中跌落運動中的速度一直增大，此時會受空穴的影響，使得空泡出現(xiàn)，減小其在水中運動過程受到的阻力。結(jié)構(gòu)物在持續(xù)跌落運動中因為氣泡發(fā)生潰散，圓柱體承受力也隨之不斷發(fā)生變化，甚至發(fā)生短暫的震蕩。圖5展示了結(jié)構(gòu)物水中跌落加速度隨時間的變化。在60 ms到70 ms間出現(xiàn)最大加速度時，結(jié)構(gòu)物表面的空穴與其脫離，水流作用力達到最大。隨著結(jié)構(gòu)物周圍氣泡持續(xù)破散，隨后物體的加速度持續(xù)降低，水流作用力逐漸減小，最后趨于平緩。當(dāng)結(jié)構(gòu)物的水平速度和水流速度相等時，水平方向的位移也呈現(xiàn)出直線的狀態(tài)。

圖5 水平方向加速度的時間歷程Fig. 5 Horizontal acceleration time history diagram

3.2 水中跌落階段分析

圓柱體入水相關(guān)的實驗研究相當(dāng)復(fù)雜，物體水中跌落運動是研究的重點。筆者主要研究完全入水后結(jié)構(gòu)物跌落軌跡的狀況。實驗截取了4個典型入水時刻的數(shù)值結(jié)果和試驗結(jié)果，具體實驗結(jié)果如圖6。

圖6 數(shù)值仿真與實驗記錄圖像對比Fig. 6 Comparison of numerical simulation and experimental results

因為流體會粘附在圓柱體的周遭，“空穴”效應(yīng)和氣泡現(xiàn)象就會在其入水的過程中顯現(xiàn)出來，如圖6。在結(jié)構(gòu)物跌落入水時，周圍的氣體也一并被壓進水里，在55 ms的時候，會有明顯的空穴效應(yīng)顯現(xiàn)。圓柱體運動到75 ms時，空穴與圓柱體分開并破散，此時水面出現(xiàn)波動現(xiàn)象。受從右向左的水流速度的影響，在95 ms時，圓柱體周身的氣泡在自由液面波浪與圓柱體最初位置關(guān)系在水平方向抵消時慢慢減小。對比發(fā)現(xiàn)，數(shù)值模擬很好演示出圓柱體在水流中跌落的運動過程，同時跌落空化的生成、瓦解等主要的水動力特征也很好的被再現(xiàn)出來。

通過準(zhǔn)確的判斷結(jié)構(gòu)物落點的區(qū)域，分析其水平方向的位移變化，確定其在水底分散情形，更加科學(xué)的探究結(jié)構(gòu)物在水流中跌落的運動軌跡。圖7展示的是結(jié)構(gòu)物水平方向的位移和時間變化的關(guān)系。采用工程偏差及實驗值表現(xiàn)方式展現(xiàn)出結(jié)果。工程偏差通過總結(jié)大量的物體落水運動試驗得出。

圖7 水平方向位移的時間歷程Fig. 7 Horizontal displacement time history diagram

通過采用FLUENT的數(shù)值模擬方法對圓柱體在水流中跌落過程進行仿真，仿真結(jié)果出現(xiàn)在試驗結(jié)果工程誤差的范圍內(nèi)，符合工程誤差的要求，因此仿真結(jié)果與試驗結(jié)果趨于相同。

結(jié)構(gòu)物在水流中下落時，隨著水平方向速度的增加，曲線變化率也由明顯到漸漸趨于一致，當(dāng)速度達到0.1 m/s，其速度不再增加，最后曲線漸漸趨于直線。從圖7的仿真數(shù)據(jù)能大致得出每個時間段圓柱體的落點區(qū)域范圍。

圓柱體水中跌落的運動軌跡見圖8。結(jié)構(gòu)物位置被清晰地展示出來，進而確定水流中結(jié)構(gòu)物具體跌落位置的范圍，對運動軌跡可以進行明確地分析和預(yù)測。

據(jù)圖8，仿真的運動軌跡和實驗結(jié)果非常切合，且處于合理誤差內(nèi)。

圖8 水平方向與豎直方向位移之間的關(guān)系Fig. 8 Relation between horizontal and vertical displacement

通過對圓柱體在實驗水流中軌跡進行擬合，可以得出結(jié)構(gòu)物水流中跌落的軌跡方程為：

y=0.091 6x3-1.648 7x2+20.221 3x+4.38

(24)

式(24)雖然只涉及具體情形下圓柱體水中跌落的軌跡方程，但針對跌落速度與來流速度等因素不同狀況的結(jié)構(gòu)物水中跌落過程試驗可提供借鑒，另一方面也有利于結(jié)合相似準(zhǔn)則，用于不同的水中跌落情況探究。

4 結(jié) 論

為了獲得結(jié)構(gòu)物水流中跌落軌跡和散布規(guī)律，采用CFD對圓柱體自由落體跌入水中情況進行了數(shù)值模擬。采用明渠模型真實捕捉自由液面，利用可實現(xiàn)的k-ε模型來模擬湍流，數(shù)值模擬結(jié)果與實驗結(jié)果吻合得較好，符合水動力特征，驗證了模型的可行性。研究結(jié)果表明，圓柱體在水流中脫離空穴時受水動力較大，影響跌落軌跡和散布落點。通過采集圓柱體重心位置的三維數(shù)據(jù)，結(jié)合誤差分析，擬合出水流中跌落軌跡方程，為后續(xù)研究不同入水速度和不同流速多工況的結(jié)構(gòu)物水中跌落提供參考依據(jù)。