VOF方法模擬水面振蕩流場

廖 斌，陳善群

VOF方法模擬水面振蕩流場

廖 斌，陳善群

（安徽工程大學機械與汽車工程學院，安徽蕪湖 241000）

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

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

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

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

1 VOF數值方法

1．1 控制方程［5］

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

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

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

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

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

1．2 邊界條件

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

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

1．3 方程離散

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

圖1 交錯網格單元Fig．1 Staggered grid cell

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

其中對流項FUX，FUY和VISX分別表示如下：

式中sgn代表符號函數。

式（6）與式（7）中的系數α是控制迎風差分量的參數。

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

連續方程（4）的差分形式為

1．4 方程求解［6］

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

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

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

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

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

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

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

2 計算結果與分析

2．1 水面振蕩計算模型

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

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

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

2．2 計算結果與分析

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

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

圖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 結 語

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

［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］ 郭人豪．以等位函數法求解含自由液面之流場［D］．中國臺灣：逢甲大學，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］．海南師范大學學報，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．

（編輯：羅玉蘭）

“水工巖石力學的理論與工程應用實踐”項目成果通過鑒定

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

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

鑒定委員會專家聽取了水利部巖土力學與工程重點實驗室副主任鄔愛清教授代表項目組所作的成果匯報，審閱了相關技術報告，經質詢與討論，最終形成了成果鑒定意見。與會專家一致認為“水工巖石力學的理論與工程應用實踐”的研究成果在理論與工程實踐結合方面總體上處于國際領先水平，建議加強成果凝煉，在國際上進一步介紹與推廣。

（摘自《長江水利科技網》）

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

安徽工程科技學院引進人才科研啟動基金（2009YQQ009）

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