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。