基于k-ε湍流模型的結構物水流中跌落過程數值分析

王永虎， 林天龍

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

0 引 言

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

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

1 問題描述

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

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

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

參數名稱數值參數名稱數值質量/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的三維數值實驗水池。非結構四面體網格作為數值模型的整體網格的基礎網格，結構物入水跌落軌跡作為實驗的側重點。為了提高實驗效率，對圓柱體周圍和氣液交界面用進一步加密的密度盒網格，而在結構物撞擊區之外的區域(也就是流固耦合區域)，采用粗網格的方式進行實驗。為了模擬使得水流速度恒定，將流速為0.1 m/s的速度入口作為水流入口的邊界條件，將壓力出口作為出口水流的邊界條件，數值水池的上下面的入口邊界條件設為速度值為0，其余邊界條件均設為對稱邊界條件。

2 軟件介紹、理論基礎和計算方法

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

2.1 FLUENT軟件介紹

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

2.2 控制方程

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

(1)

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

動量守恒方程為：

(2)

(3)

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

2.3 VOF模型

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

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

(4)

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

2.4 明渠模型

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

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

(5)

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

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

2.5 動網格的使用

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

慣性坐標系中，重心移動的運動控制方程為：

(6)

2.6 湍流模型

使用可實現的k-ε模型的湍流模型，可實現的k-ε模型包含湍流黏度的可選擇公式，在改良的運輸方程中，湍流耗散率ε來自于有關平方渦度波動運輸的精確方程。

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

ρε-YM+Sk

(7)

(8)

式中：

(9)

(10)

(11)

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

2.7 SIMPLEC算法

(12)

(13)

(14)

式中：

(15)

(16)

式中：p′為網格壓力修正值。

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

(17)

式中：

(18)

2.8 二階迎風格式

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

(19)

式中：

(20)

(21)

a+=max(a,0)

(22)

a-=min(a,0)

(23)

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

(24)

2.9 計算方法

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

3 計算結果及分析

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

3.1 整體階段的速度分析

空中圓柱體只在重力的作用下跌落水中，然后在水中持續運動。在整個過程中圓柱體先受到重力的作用，入水之后逐漸轉變為重力、浮力以及水流產生的水動力等復雜的合力，單一力到復雜力的演變過程，致使其各方向的運動速度也隨其情況的變動產生變動。

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

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

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

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

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

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

3.2 水中跌落階段分析

圓柱體入水相關的實驗研究相當復雜，物體水中跌落運動是研究的重點。筆者主要研究完全入水后結構物跌落軌跡的狀況。實驗截取了4個典型入水時刻的數值結果和試驗結果，具體實驗結果如圖6。

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

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

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

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

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

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

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

據圖8，仿真的運動軌跡和實驗結果非常切合，且處于合理誤差內。

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

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

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

(24)

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

4 結 論

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