VOF方法模擬水面振蕩流場

廖 斌，陳善群

VOF方法模擬水面振蕩流場

廖 斌，陳善群

（安徽工程大學(xué)機(jī)械與汽車工程學(xué)院，安徽蕪湖 241000）

應(yīng)用VOF（Volume of Fluid）方法結(jié)合不可壓縮流體的N-S方程和連續(xù)性方程對矩形水池中的水面振蕩流場進(jìn)行數(shù)值模擬。采用共軛剩余法結(jié)合差分法求解N-S方程和連續(xù)性方程，采用施主-受主法求解流體體積函數(shù)控制方程。計算得到一個振蕩周期內(nèi)典型時刻矩形水槽水面振蕩流場分布情況并與理論分析結(jié)果進(jìn)行比較。并將計算得到的前6個周期時刻的振蕩流場自由面分布曲線進(jìn)行比較。結(jié)果表明，數(shù)值模擬所得水面振蕩流場演化過程符合理論規(guī)律，VOF方法能夠較為精確地模擬水面振蕩流場。

水面振蕩；VOF方法；振蕩周期

常發(fā)生于水庫、湖泊及海灣的振蕩現(xiàn)象，是由于水波運(yùn)動受壩體、堤防等邊壁影響，使原本流動的流體受到阻擋后造成水面壅高。水流對壩體及堤防所造成的沖擊力，往往成為此種構(gòu)筑物遭到破壞的主要原因，因此對水面振蕩流場進(jìn)行數(shù)值模擬具有重要意義。

國內(nèi)外已有不少關(guān)于水面振蕩流場的研究，Onno Ubbink在其博士論文中采用CICSAM方法對水面振蕩自由面進(jìn)行了追蹤［1］，郭人豪采用等位函數(shù)法對水面振蕩流場進(jìn)行了數(shù)值模擬［2］，林東采用差分方法求解Boussinesq方程對封閉水池中的水面振蕩位移進(jìn)行了研究［3］。本文以矩形水槽中的水面振蕩流場為研究對象，用VOF方法［4］結(jié)合不可壓縮流體的N-S方程和連續(xù)性方程對其進(jìn)行數(shù)值模擬，試圖得到水面振蕩流場分布狀況，為模擬水面振蕩流場提供更加精確的數(shù)值方法。

1 VOF數(shù)值方法

1．1 控制方程［5］

VOF方法的控制方程為不可壓縮流體的N-S方程和連續(xù)性方程：

式中：t為時間；u，v分別為x，y方向上的流體速度分量；p為流體壓強(qiáng)；ρ為流體密度；υ為流體運(yùn)動粘滯系數(shù)；gx，gy為體積加速度分量；θ為部分單元體參數(shù)，其值在0～1之間。

為了表示自由面，VOF方法引入函數(shù)F（x，y，t），它的思想是跟蹤自由面的流體體積分?jǐn)?shù)，代表了網(wǎng)格中某種流體所占整個網(wǎng)格體積的百分比。如果F＝1，則這個網(wǎng)格單元充滿流體；如果F＝0，則該網(wǎng)格單元不含流體；如果0＜F＜1，則該網(wǎng)格單元一定含有自由面。因此通過求解F值，就可以確定自由面位于哪些單元內(nèi)。

F隨時間變化的控制方程為

方程（1）、（2）、（3）、（4）構(gòu)成了本文所用VOF方法的控制方程。

1．2 邊界條件

（1）壁面邊界條件：采用可滑移邊界條件。

（2）自由面邊界條件：滿足自由面動力邊界條件。

1．3 方程離散

采用的計算網(wǎng)格單元為交錯網(wǎng)格單元（Stag-gered Grid）。如圖1所示，速度的水平分量ui＋1／2，j定義在網(wǎng)格右界面，垂直分量vi，j＋1／2定義在網(wǎng)格上界面。壓強(qiáng)Pi，j和流體體積函數(shù)Fi，j定義在網(wǎng)格的中心。ARi＋1／2，j，ATi，j＋1／2為網(wǎng)格單元右側(cè)面與上側(cè)面可通過流體部分的邊長系數(shù)；ACi，j是單元可通過流體的面積系數(shù)，對于不含有邊界的單元上述幾何參數(shù)ARi＋1／2，j，ATi，j＋1／2和ACi，j都為1。對于部分單元體幾何參數(shù)來講，ARi＋1／2，j＝θi＋1／2，j，ATi，j＋1／2＝θi，j＋1／2，ACi，j＝θi，j。

圖1 交錯網(wǎng)格單元Fig．1 Staggered grid cell

x方向動量方程式（1）的差分形式為

其中對流項(xiàng)FUX，F(xiàn)UY和VISX分別表示如下：

式中sgn代表符號函數(shù)。

式（6）與式（7）中的系數(shù)α是控制迎風(fēng)差分量的參數(shù)。

y方向動量方程式（2）的差分形式可以仿x方向一樣寫出。

連續(xù)方程（4）的差分形式為

1．4 方程求解［6］

采用共軛剩余法［7］求解差分方程式（5）和式（9）。一個時間步的求解步驟如下：

（1）用前一時刻的流場結(jié)果代入N-S方程的顯式差分格式（5），求出新時刻的流場近似值。

（2）通過迭代調(diào)整每一單元的壓力，使得對內(nèi)部流體單元滿足連續(xù)方程的隱式差分方程式（9）。若式（9）的右端非零項(xiàng)用S表示，則每次迭代時壓力的調(diào)整量δp為

對表面單元滿足自由表面的動力邊界條件，該條件的滿足是通過表面單元的壓力pi，j等于自由表面處的壓力ps與插值單元的壓力pn的線性插值來實(shí)現(xiàn)的，這時的壓力調(diào)整量δp為

（11）式中的幾何參數(shù)dc，ds的定義如圖2所示。

圖2 表面單元壓強(qiáng)插值示意圖Fig．2 Schematic diagram of pressure interpolation in surface cell

（3）由于F函數(shù)是步進(jìn)函數(shù)，不能采用一般的差分方式。本文采用施主與受主單元模型計算每一單元新的F值，構(gòu)造新的自由面。

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

2．1 水面振蕩計算模型

綜合參照文獻(xiàn)［1］、［2］給定的水面振蕩計算模型，本文數(shù)值模擬計算在笛卡爾坐標(biāo)系下，給定水槽長度L＝0．10 m，靜水深H0＝0．050 m，初始時間步長δt＝1．0×10－4s。采用的計算網(wǎng)格為x方向160個，y方向104個，網(wǎng)格大小δx＝6．25×10－4m，δy＝6．25×10－4m，初始振幅為0．005 m。水的密度ρ＝1 g／cm3，重力加速度gy＝980 cm／s2，運(yùn)動粘滯系數(shù)υ＝1．002×10－2cm2／s，迭代控制精度ε＝10－4；迎風(fēng)差分量控制參數(shù)α＝1．0。水面振蕩的理論周期為T＝0．373 9［8］s。

初始水面的給定如圖3所示，其中左邊壁的水位給定HL＝0．055 m，右邊壁的水位給定HR＝0．045 m。整個計算區(qū)域水位H＝HL＋xi×；速度場u＝v＝0。

圖3 矩形水槽水面振蕩流場初始條件示意圖Fig．3 Initial condition for the liquid oscillation flow in rectangu lar tank

2．2 計算結(jié)果與分析

圖4給出第一個振蕩周期典型時刻的矩形水槽水面振蕩流場分布情況。t＝0時，振蕩水體處于靜止?fàn)顟B(tài)，此時水體只有勢能而沒有動能。水體起動之后，左邊水面逐漸下降而右邊水面逐漸上升，同時水體的勢能向動能轉(zhuǎn)化，并產(chǎn)生向右運(yùn)動的水流。t＝1／4周期時，水體的勢能完全轉(zhuǎn)化為動能，此時水體水面幾乎水平，向右的動能最大。由于慣性作用，水流將繼續(xù)向右擠壓，使得右邊水面逐漸升高而左邊水面逐漸下降，水體的動能又將逐漸轉(zhuǎn)化為勢能。t＝1／2周期時，水體的動能完全轉(zhuǎn)化為勢能，使得右邊水面達(dá)到最高而左邊水面達(dá)到最低。后半周期水面振蕩流場演化過程與前半周期形似但方向相反。

將前6個振蕩周期內(nèi)左邊壁水面高度隨時間的演化過程記錄，如圖5所示。與理論振蕩曲線作比較，結(jié)果顯示6個周期內(nèi)，計算曲線與理論曲線的振蕩周期基本相同，整個濕墻高度的變化趨勢與理論曲線相吻合，振幅方面與理論值存在些許誤差。究其原因，應(yīng)為考慮水體的運(yùn)動粘滯系數(shù)所致。圖6給出前6個周期時刻的振蕩流場自由面分布曲線，從圖中不難看出6個周期時刻的自由面曲線大部分相互重疊，從而說明本文采用的數(shù)值計算方法具有較高的精確度。

圖4 第一周期典型時刻矩形水槽水面振蕩流場分布圖Fig．4 Velocity distribution of the liquid oscillation flow in rectangular tank at the typical time of the first period

圖5 左邊壁濕墻高度隨時間演化圖Fig．5 Evolution of the wet wall height at the left boundary over time

圖6 前6個周期時刻的振蕩流場自由面分布曲線Fig．6 Distribution curves at free surface of the liquid oscillation flow for the first six periods

3 結(jié) 語

本文以矩形水槽水面振蕩流場為研究對象，用VOF方法結(jié)合不可壓縮流體的N-S方程和連續(xù)性方程對其進(jìn)行了數(shù)值模擬，得到一個振蕩周期內(nèi)典型時刻矩形水槽水面振蕩流場分布情況，整個流場演化過程符合理論水面振蕩規(guī)律。將前6個振蕩周期內(nèi)左邊壁水面高度隨時間的演化曲線與理論振蕩曲線相比較，結(jié)果顯示兩者振蕩周期基本相同，振幅存在些許誤差，應(yīng)為水體的運(yùn)動粘滯性所致。將前6個周期時刻的振蕩流場自由面分布曲線相比較，振蕩自由面曲線大部分相互重疊，從而說明本文采用的數(shù)值計算方法具有較高的精確度。

［1］ ONNO UBBINK．Numerical Prediction of Two Fluid Sys-temswith Sharp Interfaces［D］．Great Britain：The Uni-versity of London and Diploma of Imperial College，1997．

［2］ 郭人豪．以等位函數(shù)法求解含自由液面之流場［D］．中國臺灣：逢甲大學(xué)，2003．（GUO Ren-hao．A Level Set Technique Applied to Free Surface Flow［D］．China，Tai- wan：Feng Chia University，2003．（in Chinese））

［3］ 林 東．封閉水池中的水面振蕩［J］．海南師范大學(xué)學(xué)報，2004，17（1）：19－23．（LIN Dong．Water Wave Evo-lution in Closed Rectangular Basin［J］．Journal of Hainan Normal University，2004，17（1）：19－23．（in Chinese））

［4］ HIRT CW，NICHOLSB D．Volume of Fluid Method for the Dynamics of Free Boundary［J］．Journal of Computa-tional Physics，1981，93：201－225．

［5］ TORREY M D，CLOUTMAN L D，MJOLSNESSR C，et al．NASA-VOF2D：A Computer Program for Incompressi-ble Flowswith Free Surfaces［R］．United States：Los Al-mos Scientific Laboratory，LA-10612-MS，1985．

［6］ HOTCHKISSR S．Simulation of Tank Draining Phenome-na with the NASA SOLA-VOF Code［R］．United States：Los Almos Scientific Laboratory，LA-8163-MS，1979．

［7］ WILLIAM C S，PLOTR K S，JOSEPH B K．Precondi-tioned Conjugate-Residual Solvers for Helmholtz Equa-tions in Non-hydrostatic Models［J］．Monthly Weather Review，1997，（125）：587－599．

［8］ RAAD PE，CHEN S，JOHNSON D B．The Introduction of Micro Cells to Treat Pressure in Free Surface Fluid Flow Problems［J］．Journal of Fluids Engineering，1995，117（4）：683－690．

（編輯：羅玉蘭）

“水工巖石力學(xué)的理論與工程應(yīng)用實(shí)踐”項(xiàng)目成果通過鑒定

2011年4月12日，湖北省科學(xué)技術(shù)廳在武漢市召開了由長江科學(xué)院和長江勘測規(guī)劃設(shè)計研究有限責(zé)任公司完成的“水工巖石力學(xué)的理論與工程應(yīng)用實(shí)踐”成果鑒定會。

成果鑒定委員會的專家有總參科技委錢七虎院士、同濟(jì)大學(xué)孫鈞院士、解放軍后勤工程學(xué)院鄭穎人院士、長江水利委員會鄭守仁院士、中國地質(zhì)大學(xué)（武漢）唐輝明教授、湖北省水利廳陳斌教授、武漢大學(xué)談廣鳴教授、華中科技大學(xué)周建中教授、武漢科技大學(xué)鐘冬望教授。

鑒定委員會專家聽取了水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室副主任鄔愛清教授代表項(xiàng)目組所作的成果匯報，審閱了相關(guān)技術(shù)報告，經(jīng)質(zhì)詢與討論，最終形成了成果鑒定意見。與會專家一致認(rèn)為“水工巖石力學(xué)的理論與工程應(yīng)用實(shí)踐”的研究成果在理論與工程實(shí)踐結(jié)合方面總體上處于國際領(lǐng)先水平，建議加強(qiáng)成果凝煉，在國際上進(jìn)一步介紹與推廣。

（摘自《長江水利科技網(wǎng)》）

Simulation of W ater Oscillation Flow Field by VOF M ethod

LIAO Bin，CHEN Shan-qun
（College of Mechanical and Automotive Engineering，Anhui Polytechnic University，Wuhu 241000，China）

The VOF（Volume of Fluid）method associated with Navier-Stokes equations and continuity equations of incompressible fluid was used to simulate the oscillation flow field of liquid in a rectangular tank．Conjugate residu-al（CR）algorithm and finite differencemethod were employed to solve Navier-Stokes equations and continuity equa-tions，and donor-acceptor scheme was used to discretize the governing equation of fluid volume function．Thewave position and velocity vectors during one period obtained by the computationswere compared with theoretical results．Moreover，distribution curves of the wave position and velocity vectors for the first six oscillation periods of the liq-uid oscillation field were compared．The comparisons showed that the simulated evolutionary process ofwater wave oscillation field is consistentwith theoretical results．Therefore，VOFmethod can be applied to simulate the oscilla-tion flow field accurately．

water oscillation；VOFmethod；oscillation period

TV131．2

A

1001－5485（2011）05－0027－04

2010-06-28

安徽工程科技學(xué)院引進(jìn)人才科研啟動基金（2009YQQ009）

廖 斌（1985-），男，江西撫州人，碩士，助教，主要從事計算流體力學(xué)研究，（電話）13515534139（電子信箱）liaobinfluid＠126．com。