基于VOF模型的結(jié)構(gòu)物出水過程數(shù)值模擬

鄒 星，李海濤，宗 智

(大連理工大學(xué)船舶工程學(xué)院，遼寧 大連 116024)

結(jié)構(gòu)物出水在許多工程領(lǐng)域都具有非常重要的意義，一直是流體力學(xué)工作者研究的重要內(nèi)容。當(dāng)物體從水中上升，穿過自由液面并進(jìn)入空氣時，周圍介質(zhì)的密度會發(fā)生一個很大的階躍，從水的密度998.2 kg/m3突然變?yōu)?.225 kg/m3，這種突變類似于氣體動力學(xué)中的強間斷問題，就像穿過了一層很強的激波一樣[1]，在水面附近會發(fā)生一系列復(fù)雜的物理現(xiàn)象。

出水問題的理論研究可追溯到 KARMAN[2]的工作，他在1929年提出了計算物體撞擊載荷的理論模型，這一思想至今不失其光彩。HAVELOCK[3]用浸沒于水中的運動的偶極子來計算水平放置的圓柱體在自由表面下的水平運動，給出了瞬時的受力結(jié)果。LIJU等[4]利用邊界元方法(BEM)模擬了一個對稱剛體在垂直穿過自由面時的液面變化，并利用實驗進(jìn)行了驗證。

國內(nèi)學(xué)者在結(jié)構(gòu)物出水方面也做了很多工作。魯傳敬[5]使用勢流模型得到細(xì)長體出入水問題的解析解。蔡軍輝等[6]通過解N-S方程，采用MAC法進(jìn)行平頭圓柱體垂直出水黏性流動的數(shù)值模擬。在實驗方面，張軍等[7]采用 TRPIV技術(shù)對鈍頭回轉(zhuǎn)體垂直出水進(jìn)行測量，揭示了出水過程中流動結(jié)構(gòu)及其演變。

VOF模型通過求解單獨的動量方程和處理穿過區(qū)域的每一流體的體積分?jǐn)?shù)來模擬兩種或3種不能混合的流體[8]。在VOF模型中，兩種流體沒有互相穿插，對增加到模型里的每一附加相，就引進(jìn)一個變量，即計算單元里的相的體積分?jǐn)?shù)[9]。在每個控制體積內(nèi)，所有相的體積分?jǐn)?shù)為1。所有變量及其屬性的區(qū)域被各相共享并且代表了體積平均值，只要每一相的體積分?jǐn)?shù)在每一位置是可知的，就能得出不同位置不同時刻下各相的分布情況。在結(jié)構(gòu)物穿過水和空氣的交界面時，涉及到空氣和水兩相流動的問題。在每個單元中，用a1代表空氣的體積分?jǐn)?shù)，a2代表水的體積分?jǐn)?shù)，通過求解任何一相體積分?jǐn)?shù)的連續(xù)方程即可實現(xiàn)對水和空氣之間界面的追蹤，知道不同時刻下液面的變形情況。

筆者基于VOF模型，在Fluent中運用動網(wǎng)格技術(shù)，將控制回轉(zhuǎn)體運動的自定義函數(shù)(UDF)與流動控制方程耦合，模擬得到回轉(zhuǎn)體在出水過程中的壓力場和速度場。然后與實驗及文獻(xiàn)比較，在物體頭部帶水和尾部拖水，渦量分布方面都吻合得較好，為出水過程的進(jìn)一步研究奠定了基礎(chǔ)。

1 數(shù)值計算方法

1.1 計算模型及邊界條件

為了便于與實驗進(jìn)行比較，筆者的計算模型類似張軍的TR-PIV實驗，只是運動速度改為5.1 m/s恒定不變，比其擴大了10倍，但這不影響壓力、速度等場量的分布情況。

計算模型是一鈍頭回轉(zhuǎn)體，頭部直徑D=30 mm，高 h=60 mm，計算域為 500 mm×800 mm，水深500 mm。網(wǎng)格劃分在 Fluent專業(yè)前處理器Gambit上完成。網(wǎng)格采用非結(jié)構(gòu)化三角形網(wǎng)格，在回轉(zhuǎn)體周圍進(jìn)行網(wǎng)格加密，網(wǎng)格總數(shù)為81 822個，計算時間為17.5 h，如圖1所示。

邊界條件設(shè)定:回轉(zhuǎn)體、計算域兩側(cè)和底部為固壁面無滑移條件，上端開口為壓力出口條件，通過外加UDF控制回轉(zhuǎn)體勻速直線運動。

1.2 控制方程

結(jié)構(gòu)物出水是一種非定常繞流問題，控制方程主要有連續(xù)性方程和N-S方程。

連續(xù)性方程為:

N-S 方程為[10]:

由于回轉(zhuǎn)體在出水過程中可能導(dǎo)致負(fù)的正壓力，為了更好地模擬湍流的物理規(guī)律，選擇Realizable κ-ε模型，關(guān)于κ和ε的輸運方程如下[11]:

其中，σκ=1.0，σε=1.2，C2=1.9，由 Fluent湍流模型在控制面板上提供。C1=max[0.43，η/

VOF模型通過求解水或空氣體積分?jǐn)?shù)的連續(xù)方程來實現(xiàn)相與相間的界面追蹤，對第q相，其方程如下:

式中，mqp為第q相對第p相的質(zhì)量傳輸;mpq為第p相對第q相的質(zhì)量傳輸;Saq為源項。主相容積比率的計算基于如下約束:

如果水和空氣分別用下標(biāo)1、2表示，則每一單元的密度由式(7)給出:

1.3 計算參數(shù)設(shè)置

控制動網(wǎng)格變化采用彈簧近似光順模型和局部重劃模型，并激活尺寸函數(shù)控制，防止計算過程中出現(xiàn)扭曲度過大的網(wǎng)格。壓力和速度耦合采用PISO算法。壁面處理采用增強的壁面函數(shù)法。

2 計算結(jié)果和分析

圖2為鈍頭回轉(zhuǎn)體出水過程的速度矢量場，按等時間間隔選取了6個不同時刻。從圖2中可以看出，在開始階段，速度的分布范圍很廣，囊括了以回轉(zhuǎn)體最大尺度為直徑的圓，頭部由于排擠作用，速度呈放射狀，最大速度出現(xiàn)在底部中心區(qū)域，左右兩側(cè)速度較小，并有明顯的回流區(qū)，左側(cè)為逆時針，右側(cè)為順時針，這是回轉(zhuǎn)體對水的拖帶作用結(jié)果。隨著物體的不斷上升，尾部速度分布范圍拉長，兩側(cè)的回流區(qū)也越來越不明顯，在物體穿越水面進(jìn)入空氣后(t=0.082 s以后)，沿兩側(cè)壁面出現(xiàn)較大速度，是因為回轉(zhuǎn)體將水帶出，水受黏性及慣性力作用繼續(xù)向上運動。

圖3為鈍頭回轉(zhuǎn)體出水過程的渦量分布圖，同速度矢量場一樣，選取了6個時刻，觀察其變化過程。從圖3中可以看出，在回轉(zhuǎn)體底部兩個拐角處分別分布一渦，并隨著運動繼續(xù)，渦變?yōu)閹?，并逐漸從壁面脫落，形成渦列。在頭部也可以看到左右分布著兩個渦。

回轉(zhuǎn)體在出水的瞬間，流體會因為黏性及慣性作用，附著在壁面上隨物體一起運動，形成薄薄的一層水膜，之后由于側(cè)面壓力減小，水受到的粘壓力不足以彌補重力，流體會慢慢滑落、凹陷，從而水膜破滅，壁面與空氣接觸。圖4為回轉(zhuǎn)體出水瞬間及以后的兩相圖。從圖4中可以看出回轉(zhuǎn)體拖帶水的情況及變化趨勢，開始階段，頭部帶出大量液體，呈山體狀，隨著物體繼續(xù)上升，水體逐漸變細(xì)，中部向內(nèi)凹，最后呈帶狀。水膜破滅是從兩側(cè)壁面開始，先是出現(xiàn)一些不規(guī)整的小突起，然后逐漸擴大到頭部，頭部的水分又由于重力的作用滑落到側(cè)壁，最后同底部帶狀水體融合。

圖2 不同時刻下回轉(zhuǎn)體的速度矢量圖

圖3 不同時刻下回轉(zhuǎn)體的渦量圖

圖4 回轉(zhuǎn)體出水瞬間及以后的兩相變化示意圖

回轉(zhuǎn)體在整個出水過程中，會同時受到壓力和流體黏性力的作用。計算過程中發(fā)現(xiàn)，兩側(cè)的黏性系數(shù)及底面的壓力系數(shù)一直很小，幾乎為零，而回轉(zhuǎn)體在出水過程中所受總的流體力變化及頂部的壓力系數(shù)如圖5、圖6所示。從圖中可以看出，總流體力和頂部壓力系數(shù)變化趨勢一致，在水中起伏很大，當(dāng)回轉(zhuǎn)體穿過水面進(jìn)入空氣后，兩者趨于穩(wěn)定，在很小范圍內(nèi)浮動。

同張軍等的TR-PIV實驗結(jié)果比較可以看出，鈍頭回轉(zhuǎn)體在出水過程中所引起的速度場、渦量場分布是吻合的，出水瞬間頭部、尾部拖帶水分情況也符合實際情況，這說明計算結(jié)果是準(zhǔn)確的，可以為出水過程的進(jìn)一步研究提供依據(jù)。

圖5 回轉(zhuǎn)體在出水過程中所受總流體力變化情況

圖6 回轉(zhuǎn)體在出水過程中頂部壓力系數(shù)變化情況

3 結(jié)論與展望

筆者利用VOF模型，通過編寫自定義函數(shù)(UDF)控制回轉(zhuǎn)體運動，計算得到了鈍頭回轉(zhuǎn)體在出水過程中的速度矢量場、渦量場的分布情況，以及出水后物體所帶出水分的變化情況，為結(jié)構(gòu)物出水的深入研究提供了依據(jù)。

結(jié)構(gòu)物出水是一個涉及到兩相流、流固耦合的非線性復(fù)雜問題，有很多情況有待深入研究。從渦量分布結(jié)果看，開始階段及往后較長時間內(nèi)，底部兩個拐角處一直分布著一對點渦，這是否與物體形狀、流體性質(zhì)等有關(guān)還有待研究。

[1] 楊曉光.潛射導(dǎo)彈水下發(fā)射及出水過程三維數(shù)值研究[D].哈爾濱:哈爾濱工業(yè)大學(xué)圖書館，2009.

[2] KARMAN V T，WATTENDORF F L.The impact on seaplane floats during landing[J].National Advisory Committee for Aeronautics，1929(TN321):353-358.

[3] HAVELOCK T.The forces on a circular cylinder submerged in a uniform stream[J].Proc R Soc London Ser A，1936(157):526-534.

[4] LIJU P Y，MACHANE R，CARTELLIER A.Surge effect during the water exit of an axis-symmetric body traveling normal to a plane interface:experiments and BEM simulation[J].Experiments in Fluids，2001(31):241-248.

[5] 魯傳敬.軸對稱細(xì)長體垂直出入水[J].水動力學(xué)研究與進(jìn)展，1990，5(4):387-393.

[6] 蔡軍輝，何友聲，葉取源.平頭圓柱體垂直出水粘性流動的數(shù)值模擬[J].水動力學(xué)研究與進(jìn)展，1990，5(4):115-121.

[7] 張軍，李英浩，金朋壽.垂直及斜出水流場的二維及三維 TR-PIV 試驗[J].船舶力學(xué)，2005，9(2):9-17.

[8] 童亮，于罡，彭政，等.基于VOF模型與動網(wǎng)格技術(shù)的兩相流耦合模擬[J].武漢理工大學(xué)學(xué)報:信息與管理工程版，2008，30(4):525-528.

[9] 江帆，黃鵬.Fluent高級應(yīng)用與實例分析[M].北京:清華大學(xué)出版社，2008:140-164.

[10] 夏國澤，馬乾初.船舶流體力學(xué)[M].武漢:華中科技大學(xué)出版社，2003:42-93.

[11] 王福軍.計算流體動力學(xué)分析:CFD軟件原理與應(yīng)用[M].北京:清華大學(xué)出版社，2004:105-126.