應(yīng)用VOF方法的水平圓柱入水?dāng)?shù)值模擬

陳宇翔，郜冶，劉乾坤

(哈爾濱工程大學(xué)航天與建筑工程學(xué)院，黑龍江哈爾濱150001)

物體入水現(xiàn)象如船舶航行砰擊、水上飛機降落、空投魚雷入水及垮塌橋梁落水廣泛存在于自然界.物體入水問題是典型的非線性問題，其入水過程特別是初始瞬間物體會受到巨大的沖擊載荷，以風(fēng)浪中航行的船舶為例，船體砰擊會造成船體的震動、減緩船速，嚴(yán)重的話還會導(dǎo)致結(jié)構(gòu)損壞;空投入水的魚雷則可能因為入水沖擊造成魚雷結(jié)構(gòu)損壞、彈道失控等嚴(yán)重問題.物體入水問題在現(xiàn)代許多軍事及民用領(lǐng)域的廣泛存在及其重要性激發(fā)了眾多科學(xué)家及工程人員的研究熱情［1］.因為圓柱、楔等形狀結(jié)構(gòu)簡單，?？勺鳛樗柩芯课矬w的簡化外形，所以研究者們常常以它們作為研究對象.20世紀(jì)90年代早期以前，大部分的研究者［2-4］均基于勢流理論假設(shè)研究過楔等形狀物體的入水問題.然而，邊界元方法在處理復(fù)雜自由表面問題方面不夠強?。?］.自1981年Hirt與 Nichols［6］等提出 VOF方法以來，帶自由表面的兩相粘性不可壓縮流體運動的數(shù)值計算技術(shù)得到了迅速的發(fā)展.近30年來，VOF方法在數(shù)值計算物體入水領(lǐng)域得到了廣泛的應(yīng)用.Arai［7］假設(shè)流動為無粘不可壓流動，應(yīng)用VOF方法對二維楔、圓柱、船艏的砰擊入水問題進行了模擬;Schumann［8］用類似的方法模擬了艏的砰擊入水.Xing-Kaeding［9］基于VOF方法，應(yīng)用Comet軟件對二維水平圓柱的粘性不可壓出入水運動進行了模擬，取得了與實驗結(jié)果符合較好的結(jié)果.本文應(yīng)用Fluent軟件中的VOF方法結(jié)合動網(wǎng)格技術(shù)對二維圓柱的入水問題進行了數(shù)值模擬計算，實現(xiàn)了圓柱運動的剛體動力學(xué)方程和氣液兩相流動的N-S方程的耦合求解.本文數(shù)值模擬計算結(jié)果與Greenhow和Lin［10］的實驗對比表明，圓柱入水的重要現(xiàn)象如飛濺現(xiàn)象、空氣墊效應(yīng)［11］等自由表面形狀的變化得到了很好的呈現(xiàn)，圓柱運動軌跡得到了很好的追蹤.同時，還對粘性對圓柱入水運動問題的影響進行了研究.

1 問題的描述及數(shù)值水池的建立

如同文獻［10］中Greenhow和Lin設(shè)計的實驗，一個直徑11 cm的零浮力(neutral buoyancy)水平圓柱由靜水上方某一確定位置自由落體進入水池中.而且當(dāng)圓柱在接觸水面的瞬間，其速度為3.13 m/s，隨后圓柱繼續(xù)向下運動直至到達水深為30 cm處.

由以上描述的問題，本文建立了一個寬60 cm，高60 cm(水深30 cm)的二維數(shù)值水池.為準(zhǔn)確地捕捉圓柱入水沖擊引起的自由表面變化等重要現(xiàn)象，整個計算區(qū)域采用正交四邊形網(wǎng)格，為此在圓柱邊界采用階梯近似［12］處理，并對圓柱附近的網(wǎng)格進行了加密處理，見圖1.圓柱邊界應(yīng)用無滑移壁面邊界條件，水池上邊界應(yīng)用壓力出口邊界.為了節(jié)省計算網(wǎng)格并盡可能模擬真實流動區(qū)域，采用文獻［9］中的方法，水池的左右邊界(包括靜水和空氣區(qū)域)應(yīng)用UDF定義邊界條件使壓力保持為自由表面未變化前的靜壓;水池下邊界應(yīng)用速度為零的速度入口邊界條件.

圖1 圓柱入水問題計算網(wǎng)格Fig.1 Numerical grid for the circular cylinder waterentry case

2 理論基礎(chǔ)和計算方法

圓柱入水問題是復(fù)雜的兩相流體流動和固體運動的耦合過程.圓柱的入水沖擊會突然改變兩相流場中水的壓力場和速度場，而流場的改變同時又會反作用來影響圓柱的運動.本文中，有限體積法被用于計算存在自由表面的粘性不可壓兩相流體流動;二維圓柱視為剛體，其在落水過程中為一自由度(1DOF)運動，流場形狀隨時間變化，需要應(yīng)用動網(wǎng)格.

2.1 流體流動與剛體運動方程

式中:V是流體運動速度矢量，Vb是有限體積即網(wǎng)格單元表面的運動速度.動量方程為

對于每一個有限體積單元，連續(xù)性方程為

式中:T為粘性應(yīng)力張量，I為單位張量，b為體積力矢量.當(dāng)網(wǎng)格單元移動或改變形狀時，有限體積單元滿足空間守恒律:

由于自由表面的形狀和位置是未知的，需要由求解過程給出，因此數(shù)值求解帶自由表面的兩相流動問題變得十分困難.為了模擬有自由表面的兩相流流動，本文應(yīng)用VOF(volume of fluid)方法.VOF是一種捕捉自由表面的方法，它同時對水和空氣建立控制域，通過求解附加的方程來捕捉自由面，這種方法在處理復(fù)雜的自由表面流動問題具有很大優(yōu)勢.引入的附加體積分?jǐn)?shù)輸運方程:

式中:c=Vl/V為液相體積分?jǐn)?shù).VOF方法將0＜c＜1的兩相混合流體處理為單相流，認(rèn)為其有同樣的壓力和速度值.其物性參數(shù)如密度處理為

由牛頓第二定律，流體中運動的圓柱剛體動力學(xué)方程為

式中:VC為運動時的圓柱重心的速度矢量.

2.2 網(wǎng)格運動方案

為了適應(yīng)由于圓柱的剛體運動而不斷改變的流場形狀，需要應(yīng)用動網(wǎng)格模型.為了獲得很好的數(shù)值解，特別是在對于圓柱入水沖擊引起的的強非線性自由表面變化模擬的情況下，需要很好的網(wǎng)格質(zhì)量，因此本文采用整體網(wǎng)格移動而不變形的方案［13］.本文應(yīng)用Fluent的動態(tài)層模型，由定義整個流動區(qū)域和圓柱為剛體區(qū)域，而將數(shù)值水池上下邊界定義為靜止區(qū)域來實現(xiàn)網(wǎng)格整體運動而不變形.每一當(dāng)前時間步的網(wǎng)格運動速度為上一時間步網(wǎng)格速度與當(dāng)前時間步圓柱所受的兩相流體壓力、粘性力以及自身重力所引起的速度增量之和.

2.3 計算方法

本文中，PISO算法用來耦合速度場和壓力場;VOF兩相流模型的自由面重構(gòu)方法用來捕捉自由表面的變化情況;對于流場中的湍流的描述則采用標(biāo)準(zhǔn)k-ε兩方程模型和壁面函數(shù)法;壓強差值采用PRESTO!差分格式.

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

本文應(yīng)用之前描述的Fluent動網(wǎng)格技術(shù)結(jié)合VOF方法對二維零浮力水平圓柱入水問題進行了計算，對比分析了數(shù)值計算結(jié)果與實驗結(jié)果，并且對于粘性對圓柱入水運動的影響進行了分析.

3.1 自由表面的變化

圓柱入水的瞬間，改變了其周圍壓力等流場參數(shù)，圓柱的速度突然下降，而同時與其接觸的水相應(yīng)地被突然加速，在圓柱兩側(cè)激起兩股射流，射流向斜上方上升直至其變得不穩(wěn)定.由于重力的影響，隨著時間的發(fā)展，射流形狀由直線逐漸彎曲向下.

圖2 零浮力圓柱入水自由表面變化:實驗(左)和數(shù)值(右)結(jié)果對比Fig.2 Free surface deformation of a neutrally buoyant cylinder water-entry:simulation(right)and experiment(left)

在圓柱入水的過程中，一些空氣也隨同圓柱一起被壓入水中，會形成所謂的＂空氣墊＂.由于圓柱底部傾角較大，導(dǎo)致底部的空氣逃逸較多，可壓縮空氣墊主要存在于圓柱兩側(cè)，將圓柱兩側(cè)表面與水隔開，隨著圓柱入水深度的不斷增加，空氣不斷排出，圓柱的沾濕面積不斷增加.

為了方便比較觀察實驗照片和數(shù)值計算結(jié)果，圖2給出了5個不同時刻實驗和數(shù)值模擬的水池自由表面變化情況對比圖.由圖2可知，本文計算的零浮力圓柱入水過程自由表面變化情況和實驗對比符合很好，成功地模擬顯示了圓柱入水沖擊引起的水射流形成、行進、彎曲以及入水圓柱周圍的空氣墊等自由表面的變化現(xiàn)象.

3.2 圓柱的剛體運動

進入水池之前，圓柱在空氣中自由下落，在氣流場中受到的阻力較小.圓柱接觸水面的瞬間，兩相流場中的水由于圓柱的沖擊被大大加速，壓力場、速度場特別是圓柱附近的壓力分布突然改變.相應(yīng)地，圓柱受到流場很大的阻力(包括壓差阻力、摩擦阻力)，速度大大降低.圖3給出了圓柱入水過程的無量綱阻力(流場數(shù)據(jù)積分得到的入水阻力與重力的比值)與時間的關(guān)系，零時刻圓柱受到很大的一個沖擊力，而隨著時間的進行，進入水中的圓柱周圍的壓力分布逐漸均勻，圓柱速度越來越慢，阻力的變化逐漸平緩.

圖3 圓柱入水過程的無量綱阻力時間曲線Fig.3 Time histories of non-dimensional drag force during cylinder water-entry

圖4 圓柱豎直運動時間曲線Fig.4 Time histories of vertical motion of the cylinder

圖4給出了圓柱入水過程沒入水中的圓柱位置(y軸坐標(biāo)為入水后圓柱底部最低點的位置)和入水時間關(guān)系曲線，隨著圓柱在水中的下降，圓柱受到的阻力不斷減小，速度也逐漸降低，因此圓柱位置的變化曲線也逐漸平緩.由圖4中的曲線與實驗記錄的數(shù)據(jù)點對比可知，數(shù)值解同樣和Greenhow和Lin的實驗相當(dāng)吻合(圖中偏差較大、標(biāo)記問號的一個點同時也是實驗中認(rèn)為存在問題的數(shù)據(jù)點).

3.3 粘性對入水圓柱運動的影響

由于水的粘性，入水過程中圓柱除了受到由于兩相流場的壓力場分布不均帶來的壓差阻力，還受到粘性引起的摩擦阻力.

圖5 圓柱入水豎直運動過程時間曲線(湍流和無粘流)Fig.5 Time histories of vertical motion of the cylinder during cylinder water-entry(turbulence and without turbulence model)

圖6 圓柱入水過程的無量綱阻力(湍流和無粘流對比)Fig.6 Time histories of non-dimensional drag force during cylinder water-entry(turbulence and without turbulence model)

為了考察粘性對圓柱入水問題的影響，本文計算了無粘流模型下的圓柱入水過程，并與之前基于湍流模型的數(shù)值計算結(jié)果進行比較.圖5、6分別給出了基于湍流模型和無粘流模型計算的圓柱入水過程圓柱運動軌跡和圓柱所受阻力時間曲線.可以看出，2種模型的數(shù)值模擬結(jié)果基本一致，所以入水過程中圓柱主要受氣液兩相流場壓力分布不均引起的壓差阻力，而其所受的摩擦阻力對圓柱入水過程受力和圓柱運動軌跡影響較小.

4 結(jié)束語

本文以Fluent軟件為計算工具，應(yīng)用VOF結(jié)合動網(wǎng)格技術(shù)的數(shù)值方法對水平圓柱入水過程氣液兩相流動和圓柱剛體運動問題進行了模擬，數(shù)值模擬的自由表面變化過程以及圓柱的剛體運動軌跡與經(jīng)典文獻的實驗符合得很好.從中可以看出，本文應(yīng)用的數(shù)值方法可以用來模擬物體入水特別是伴隨著強非線性自由表面變化的問題，而之前常應(yīng)用的邊界元方法不能精確模擬此類問題;本文的數(shù)值方法可以應(yīng)用來計算物體入水所受的沖擊力，進而對物體如船體結(jié)構(gòu)的設(shè)計提出指導(dǎo)意見和提供參考數(shù)據(jù).

［1］王永虎，石秀華.入水沖擊問題研究的現(xiàn)狀與進展［J］.爆炸與沖擊，2008，28(3):276-282.

WANG Yonghu，SHI Xiuhua.Review on research and development of water-entry impact problem［J］.Explosion and Shock Waves，2008，28(3):276-282.

［2］Von KARMAN T.The impact on seaplane floats during landing.NACA TN 321［R］.Washington:National Advisory Committee for Aeronautics，1929.

［3］WAGNER H.Uber stoss-und gleitvorgange an der oberflache von flussigkeiten［J］.Zeitschrift Für Angewandte Mathematik und Mechanik，1932，12(4):193-215.

［4］ZHAO R，F(xiàn)ALTINSEN O M.Water-entry of two-dimensional bodies［J］.Journal of Fluid Mechanics，1993，246:593-612.

［5］ZHU Xinying，F(xiàn)ALTINSEN O M，HU Changhong.Water entry and exit of a horizontal circular cylinder［J］.Journal of Offshore Mechanics and Arctic Engineering，2007，129:253-264.

［6］HIRT C W，NICHOLS B D.Volume of fluid(VOF)method for dynamics of free boundaries［J］.Journal of Computational Physics，1981，39:201-225.

［7］ARAI M，CHENG L Y，INOUE Y.A computing method for the analysis of water impact of arbitrary shaped bodies［J］.Journal of the Society of Naval Architects of Japan，1995，176:233-240.

［8］SCHUMANN C.Volume-of-fluid computations of water entry of bow sections［C］//Proceedings of Euromech 374.Poitiers，1998:209-217.

［9］XING-KAEDING Y，JENSEN G，PERIC M.Numerical simulation of water-entry and water-exit of a horizontal circular cylinder［C］//Proceedings of the 6th International Conference on Hydrodynamics.Perth，2004:663-669.

［10］GREENHOW M，LIN W M.Nonlinear free-surface effects:Experiments and theory.MIT report 83-19［R］.Cambridge:Massachusetts Institute of Technology，1983.

［11］陳震，肖熙.空氣墊在平底結(jié)構(gòu)入水砰擊中作用的仿真分析［J］.上海交通大學(xué)學(xué)報，2005，39(5):670-673.

CHEN Zhen，XIAO Xi.Simulation analysis on the role of air cushion in the slamming of a flat-bottom structure［J］.Journal of Shanghai Jiao Tong University，2005，39(5):670-673.

［12］XING T，STERN F.Grid generation and post-processing for computational fluid dynamics(CFD)［Z］.［s.l.］:IIHR-Hydroscience and Engineering，2006.

［13］HAD?ICA I，HENNIGB J，PERICM，et al.Computation of flow-induced motion of floating bodies［J］.Applied Mathematical Modelling，2005，29(12):1196-1210.