两相流数值模拟（第9讲）VOF方法及其应用课件.ppt

第九讲：VOF ( Volume-Of-Fluid )方法及其应用n两相流数值模拟方法的分类分类方法很多气(液)液两相流的相界面描述方法可分为两大体系：即（）界面捕捉法(Front-Capturing Methods)（）界面跟踪法(Front-Tracking Methods)汽液两相流界面描述方法的分类 界面追踪方法，直观的理解，属于Lagarange类方法如果有足够多双眼睛盯住“锁定”界面上的各个点，通过对各个点的观察、记录就可以确定整个界面的运动过程和轨迹这就是界面追踪方法产生的最初的出发点和本质这类方法中最具代表性的是：PIC (Particle-in-Cell)MAC (Mark-and-Cell)界面追踪法（Front Tracking Methods）界面捕捉法中最具有代表性的是： Level SetVOF高度函数法界面捕捉法（Front Capturing Methods）1VOF方法的基本思想： 1981年，C. W. Hirt 和B. D. Nichols在 International Journal of Computational Physics 杂志上首先正式发表了著名的VOF论文，开创性地提出了用VOF方法进行运动相界面追踪的思想，并首先用VOF方法对溃坝和浪涌（Broken Dam, Breaking bore）及Rayleigh Taylor不稳定性现象进行了成功的数值模拟模拟。

一、VOF方法的基本思想（）VOF方法用相函数(Phase Function)F取代了MAC方法中的虚拟无质量彩色粒子从这个意义上说，VOF方法可看作是MAC方法的一个变种在一种流体相（比如说，液体）中，相函数F取值为1，而在另一种流体（比如，气体或另一种液体）中取值为0；在相函数取0到1之间的数值的地方即为相界面位置相界面的取向可由界面附近各点上的F值来确定一、VOF方法的基本思想 （4） VOF方法避开了采用工程浩大的Marker点的方案，比MAC方法对计算机内存的要求低，更容易实施5） Hirt和Nichols在最初的VOF模型中，设计了类似于MAC方法和PIC方法的Staggered型差分格式即将压力和相函数（体积分数）定义在格子的中心处，而将X方向的速度定义在格子的左、右格边中点，将Y方向的速度定义在格子的上、下格边中点一、VOF方法的基本思想关于相函数的概念：关于相函数的概念：（）相函数类似于气液两相流中的截面含气率（容积含气率），表示某一相介质占据网格面积（二维）或体积（三维）的分数相函数是以一个网格为单元来定义的，与含气率的概念不同一、VOF方法的基本思想关于相函数的概念：关于相函数的概念：（）相函数就是一个介质指针，指示着占据某一望个的介质种类。

但这种指示不仅是定性的，而且是定量的对应于同一个相函数值，气液界面在网格内的形状和方位是多值的如图所示F必须而且只能在0和1之间一、VOF方法的基本思想VOFVOF方法中的基本问题：方法中的基本问题：1.气体的动量运动控制方程2.液体的动量运动控制方程3.相函数F的控制方程及其求解4.如何由相函数F的分布获得气液相界面相界面的构造问题二、VOF方法中的基本问题：三、气液两相流的动量方程两相流运动控制方程三、气液两相流的动量方程两相流运动控制方程三、气液两相流的动量方程两相流运动控制方程三、气液两相流的动量方程两相流运动控制方程四气液相界面的控制方程五气液相界面控制方程的求解方法五气液相界面控制方程的求解方法五气液相界面控制方程的求解方法与与Level SetLevel Set方法不同的是：方法不同的是：在Level Set方法中，求得了Level Set函数，由Level Set函数的等值面可立即得到相界面的位置和形状；但是，在VOF方法中，得到了相函数F的分布之后，还有一个必须解决的关键关键问题？问题？ 即如何根据相函数即如何根据相函数F的分布准确地确定相界面在每一个的分布准确地确定相界面在每一个时间层上的空间位置，也即如何实现所谓的时间层上的空间位置，也即如何实现所谓的“相界面重新构造相界面重新构造”。

六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法CFL条件，以Courant,Friedrichs,Lewy三个人的名字命名 界面重构的实施范围六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法则单位法向为六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法6.3 FLAIR6.3 FLAIR型重构技术型重构技术FLAIR型重构技术由Ashgriz和J. Y. Poo两位教授首先于1990年提出FLAIR一词的英文原文是a Flux-Line-SegmentModel for Advection and Interface Reconstruction，即对流和界面重构的通量线段模型。

该方法运用直线段来近似两个相邻网格内的界面基本的思路：由其名称本身已可以看到其基本的思路：即对任意与界面相关的两个相邻网格，通过构造一条带有倾角的直线段作为跨过该网格边界的近似界面然后计算单位时间内流过该网格边界的流体体积通量通量(用用flux表示表示)，并作为修改流体体积函数的数值流通量六气液相界面的构造方法由于界面的构造要涉及到相邻的两个网格，因此需要分成多种情况进行讨论计算根据相邻网格的相含率的不同及施主与受主网格的判断，即按通过网格边界的流体速度确定施主单元和受主单元，分别用FD和FA表示，根据FD和FA的值可以分为以下五种情况 值得首先说明的是，进行这种分类的目的，主要是为了计算通过网格边界的流体体积通量这这种通量主要是由施主网格的状种通量主要是由施主网格的状态态决定的决定的，因此这里的分类也主要以施主网格的状态为依据这种分类主要为了厘清计算思路、方便计算过程的实施六气液相界面的构造方法六气液相界面的构造方法（1）施主网格是满网格，受主网格为满网格或空网格，即fa=1.0，fb=0.0/1.0（2）施主网格是空网格，受主网格为满网格或空网格，即fa=0.0，fb=0.0/1.0六气液相界面的构造方法 针对施主和受主网格均为半网格时的上述四种情况，其判断依据为： （1）若为第一种情况应满足： （2）若为第二种情况应满足： （3）若为第三种情况应满足： （4）若为第四种情况应满足：六气液相界面的构造方法 以下将简单介绍第一种情况的输运： （1）确定界面的斜率 界面方程为： 通过比较相邻网格内的相含率fa、fb可得界面表达式的系数为： （2）确定界面的流量 根据相邻网格相邻界面的速度以及界面方程可得界面的输运为: 式中：f, Fs, f+,f-六气液相界面的构造方法 上述四种情况的判断依据为： （1）若为第一种情况应满足： （2）若为第二种情况应满足： （3）若为第三种情况应满足： （4）若为第四种情况应满足：直线斜率 ? 上述4种界面的构造最为复杂，其主要体现在界面斜率的计算。

网格（i,j）界面的计算要在垂直方向上运用两次类型3的界面构造技术，首先对网格（i,j）、（i,j-1）进行界面构造，得到斜率1，然后再对网格（i,j）、（i,j+1）进行计算，得到斜率2，最后取两斜率的均值作为网格（i,j）界面的斜率 然后，再根据界面流速计算界面流量 六气液相界面的构造方法六气液相界面的构造方法再进一步，可按类似于情形（）中的方法构造直线，作为此网格内的界面近似，然后求出经过网格边界流入受主网格的流体体积流量，同样可以归结为图所示的种情形六气液相界面的构造方法 FLAIR算法的算例剪切速度场取(x0,y0)=(0.5,0.5)，计算区域取0,1*0,1，初始界面为圆，圆心为(0.5,0.3)、半径为0.2计算网格为200*200时间步长为0.0001s 图3.3 剪切速度场数值结果（上为旋转1s，下为旋转2s）六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法PLIC界面重构技术在二维直角坐标中的实现 PLIC界面重构技术也是应用直线段对界面进行重构，但其与FLAIR不同，PLIC方法是在单个网格内对界面进行近似。

PLIC界面存在20种情况，其中4种情况较为简单，界面与坐标轴平行，其界面流量的计算非常简单其它16种情况的界面与两个坐标轴均有交点，此界面的法向为：则单位法向为分母上系数，二维时为8；三维时为32？ 通过坐标轴的旋转以及坐标转换可将各种界面转换为 的情况， 的界面存在4种情况，分别为：图3.4 PLIC最终简化类型 上述4种情况是以图3.5中的A、B两点为分界的图3.5中界面的方程为：图3.5 PLIC界面示意图 当 时，只存在图3.4中的1、2、4情况，此时界面沿法线方向移动会先后经过A、B点界面经过A点时，阴影面积为：界面经过B点时，阴影面积为：方程？若 ，则界面为情况1；若 ，则界面为情况2；若 ，则界面为情况4 当 时，只存在图3.4中的1、3、4情况，此时界面沿法线方向移动会先后经过B、A点界面经过B点时，阴影面积为：界面经过A点时，阴影面积为：若 ，则界面为情况1；若 ，则界面为情况2；若 ，则界面为情况4 C=?相函数F，fC=?相函数F，f 图3.6中阴影面积为：式中，M、N为系数，若界面为情况1，M=1、N=1；若界面为情况2，M=1、N=0；若界面为情况3，M=0、N=1；若界面为情况4，M=0、N=0。

图3.6 PLIC界面示意图 通过上式可以求得各情况时，界面与x轴的交点为： 第一种情况： 第二种情况： 第三种情况： 第四种情况： 现以情况4为例介绍目标网格向周围4个网格的流量输运图3.6 PLIC界面的第4种情况(1)网格i,j向i,j+1网格的流量输运：如果Vi,j+10，则存在输运： IF Vi,j+1*dt0，不存在输运； 如果Vi,jy1 else （3）网格i,j向i+1,j网格的流量输运： 如果Ui+1,j0，则存在输运： IF Ui+1,j*dt0，不存在输运； 如果Ui,jx1 else 其它几种情况的界面输运与情况4相近，求得各界面的输运以后就可以得到新时刻的相函数分布 PLIC算法的验证算例如下： （1）剪切流场取(x0,y0)=(0.5,0.5)，计算区域取0,1*0,1，初始界面为圆，圆心为(0.5,0.3)、半径为0.2计算网格为200*200时间步长为0.0001s 图3.7 剪切速度场数值结果（上为旋转1s，下为旋转2s） （2）旋转流场：取(x0,y0)=(0.5,0.5)，计算区域取0,1*0,1，初始界面为开口圆，圆心为(0.5,0.5)、半径为0.25、开口宽度为0.08、开口底部距圆心的距离为0.04、开口向下。

计算网格为200*200时间步长为0.0001s图3.8 旋转速度场数值结果（依次旋转0.5s、1.0s、1.5s、2s） PLIC方法在多相流中的应用： （1）气泡在水中的上升过程 计算区域为：XL=0.2m，YL=0.2m；计算所采用的网格为102*102；计算时间步长为0.0001s；初始气泡的半径为0.02m；圆心为（0.1,0.1）气水两相的物性分别为：水的密度为1000.0kg/m3、水的粘度为0.001Pa.s；空气的密度为1.266kg/m3、空气。