1、第九讲：VOF ( Volume-Of-Fluid )方法及其应用nLatticeGasLatticeBoltzmanEulerianEulerianEulerianLagrangianFront分子气体动力学模型格子气体（）法离散的介观层次格子类方法格子波尔兹曼（）方法均相模型单颗粒动力学模型气体颗粒两相流分相模型 颗粒群轨道模型拟流体模型均相模型连续介质模型方法/方法界面追踪（气体液体两相流分相模型(双流体模型)() PICTrackingFLICMACFrontCapturingVOFLevelSetEulerianLagrangian方法）类方法方法方法高度函数法线段法界面捕捉类方法方
2、法方法混合方法GhostFluidMethodPhaseField影子流体（）法其它（辅助）方法相场（）方法两相流数值模拟方法的分类分类方法很多气(液)液两相流的相界面描述方法可分为两大体系：即（）界面捕捉法(Front-Capturing Methods)（）界面跟踪法(Front-Tracking Methods)。汽液两相流界面描述方法的分类 界面追踪方法，直观的理解，属于Lagarange类方法。如果有足够多双眼睛盯住“锁定”界面上的各个点，通过对各个点的观察、记录就可以确定整个界面的运动过程和轨迹。这就是界面追踪方法产生的最初的出发点和本质。这类方法中最具代表性的是：PIC (Par
3、ticle-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
4、bore）及Rayleigh Taylor不稳定性现象进行了成功的数值模拟模拟。 一、VOF方法的基本思想VOF方法的基本思想 VOF方法在整个流场定义了一个相函数f，在每个网格中，这个函数的定义为一种流体（目标流体）的体积与网格体积的比值。相函数满足一定的输运方程，在得出相函数输运方程之前先定义一个染色函数从而在控制体上，目标流体的体积函数表达式为 由拉格朗日流体体积的特性可知 1, ,0 x ya x y t在点上有目标流体否则 , ,a x y t d 0ddt 将上述全导数展开可得输运方程为相函数f的定义式为因此可得相函数满足的输运方程为当f=0时，网格内没有目标流体，称为空网格；当f
5、=1时，网格被目标流体充满，称为满网格；当0f1时，网格含有目标流体，但未被充满，称之为界面网格。0Vtf0fVft （）VOF方法用相函数(Phase Function)F取代了MAC方法中的虚拟无质量彩色粒子。从这个意义上说，VOF方法可看作是MAC方法的一个变种。（）在一种流体相（比如说，液体）中，相函数F取值为1，而在另一种流体（比如，气体或另一种液体）中取值为0；在相函数取0到1之间的数值的地方即为相界面位置。（）相界面的取向可由界面附近各点上的F值来确定。一、VOF方法的基本思想 （4） VOF方法避开了采用工程浩大的Marker点的方案，比MAC方法对计算机内存的要求低，更容易实
6、施。（5） Hirt和Nichols在最初的VOF模型中，设计了类似于MAC方法和PIC方法的Staggered型差分格式即将压力和相函数（体积分数）定义在格子的中心处，而将X方向的速度定义在格子的左、右格边中点，将Y方向的速度定义在格子的上、下格边中点。一、VOF方法的基本思想关于相函数的概念：关于相函数的概念：（）相函数类似于气液两相流中的截面含气率（容积含气率），表示某一相介质占据网格面积（二维）或体积（三维）的分数。（）相函数是以一个网格为单元来定义的，与含气率的概念不同。一、VOF方法的基本思想关于相函数的概念：关于相函数的概念：（）相函数就是一个介质指针，指示着占据某一望个的介质种
7、类。但这种指示不仅是定性的，而且是定量的。（）对应于同一个相函数值，气液界面在网格内的形状和方位是多值的。如图所示。（）F必须而且只能在0和1之间。一、VOF方法的基本思想VOFVOF方法中的基本问题：方法中的基本问题：1.气体的动量运动控制方程2.液体的动量运动控制方程3.相函数F的控制方程及其求解4.如何由相函数F的分布获得气液相界面相界面的构造问题。二、VOF方法中的基本问题：三、气液两相流的动量方程两相流运动控制方程三、气液两相流的动量方程两相流运动控制方程三、气液两相流的动量方程两相流运动控制方程三、气液两相流的动量方程两相流运动控制方程四气液相界面的控制方程四气液相界面的控制方程五
8、气液相界面控制方程的求解方法五气液相界面控制方程的求解方法五气液相界面控制方程的求解方法与与Level SetLevel Set方法不同的是：方法不同的是：在Level Set方法中，求得了Level Set函数，由Level Set函数的等值面可立即得到相界面的位置和形状；但是，在VOF方法中，得到了相函数F的分布之后，还有一个必须解决的关键关键问题？问题？ 即如何根据相函数即如何根据相函数F的分布准确地确定相界面在每一的分布准确地确定相界面在每一个时间层上的空间位置，也即如何实现所谓的个时间层上的空间位置，也即如何实现所谓的“相界面重新构相界面重新构造造”。六气液相界面的构造方法六气液相界
9、面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法方程（）中的最小函数方程（）中的最小函数“Mins”主要是为了防止从主要是为了防止从“施主施主”网格中过多网格中过多的流出流体的流出流体;在一个时间单元内，流过（流出流入）一个网格的距离不可能在一个时间单元内，流过（流出流入）一个网格的距离不可能超过单元格的宽度！超过单元格的宽度！方程（）中的最大函数方程（）中的最大函数“Maxs”主要是为了考
10、虑在计算中可能出现的、由主要是为了考虑在计算中可能出现的、由于非目标流体的流动量超过了可能获得的非目标流体量而导致的额外附加的流于非目标流体的流动量超过了可能获得的非目标流体量而导致的额外附加的流体流动。体流动。六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法 施主-受主的思想： Hirt和Nichols在提出VOF方法的同时，给出了相函数的求解方法，两位学者设计了类似于MAC方法和PIC方法的Staggered型差分格式，即相函数定义在网格中心处，而速度定义在网格边界上。在求解的过程中采
11、用了施主-受主的思想，具体实施过程如下： 采用施主-受主型的差分对对流项进行逼近得：其中1,1,1/2/xijiji ji jiuuFuuuuuxx,1/2,1/2,/ 2i jijijuuu1/2,1/2,00iji ji jiji juuuuu1/21/2iiixxx六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法CFL条件，以Courant,Friedrichs,Lewy三个人的名字命名 界面重构的实施范围六气液相界面的构造方法六气液相界面
12、的构造方法aynxn21六气液相界面的构造方法22,/xxyxi ji ji jnnnn,xyi jjnnn则单位法向为22,/yxyyi ji ji jnnnn六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法6.3 FLAIR6.3 FLAIR型重构技术型重构技术FLAIR型重构技术由Ashgriz和J. Y. Poo两位教授首先于1990年提出。FLAIR一词的英文原文是a Flux-Line-SegmentModel for Advection and Interface Recon
13、struction，即对流和界面重构的通量线段模型。该方法运用直线段来近似两个相邻网格内的界面。基本的思路：由其名称本身已可以看到其基本的思路：即对任意与界面相关的两个相邻网格，通过构造一条带有倾角的直线段作为跨过该网格边界的近似界面。然后计算单位时间内流过该网格边界的流体体积通量通量(用用flux表示表示)，并作为修改流体体积函数的数值流通量。六气液相界面的构造方法由于界面的构造要涉及到相邻的两个网格，因此需要分成多种情况进行讨论计算。根据相邻网格的相含率的不同及施主与受主网格的判断，即按通过网格边界的流体速度确定施主单元和受主单元，分别用FD和FA表示，根据FD和FA的值可以分为以下五种情
14、况。 值得首先说明的是，进行这种分类的目的，主要是为了计算通过网格边界的流体体积通量。这种通量主要是由施主网格的状这种通量主要是由施主网格的状态决定的态决定的，因此这里的分类也主要以施主网格的状态为依据。这种分类主要为了厘清计算思路、方便计算过程的实施。六气液相界面的构造方法六气液相界面的构造方法（1）施主网格是满网格，受主网格为满网格或空网格，即fa=1.0，fb=0.0/1.0（2）施主网格是空网格，受主网格为满网格或空网格，即fa=0.0，fb=0.0/1.0六气液相界面的构造方法（1）施主网格是满网格，受主网格为满网格或空网格，即fa=1.0，fb=0.0/1.0（2）施主网格是空网格
15、，受主网格为满网格或空网格，即fa=0.0，fb=0.0/1.0六气液相界面的构造方法六气液相界面的构造方法 针对施主和受主网格均为半网格时的上述四种情况，其判断依据为： （1）若为第一种情况应满足： （2）若为第二种情况应满足： （3）若为第三种情况应满足： （4）若为第四种情况应满足：3.and. 23baabffff23.0.5baababbffandfff ff121.512/34/3aababaababffffffffff10.511/0.5abbbabffffff六气液相界面的构造方法 以下将简单介绍第一种情况的输运： （1）确定界面的斜率 界面方程为： 通过比较相邻网格内的相含率
16、fa、fb可得界面表达式的系数为： （2）确定界面的流量 根据相邻网格相邻界面的速度以及界面方程可得界面的输运为: 式中：yaxb132baabaffbhff*/ 20fs abasu*/ 20fs abasu*/bb h/su t hf, Fs, f+,f-六气液相界面的构造方法六气液相界面的构造方法 上述四种情况的判断依据为： （1）若为第一种情况应满足： （2）若为第二种情况应满足： （3）若为第三种情况应满足： （4）若为第四种情况应满足：11.122aafandf 11.22aafandf 1.122aafandf 1.22aafandf直线斜率 ? 上述4种界面的构造最为复杂，其主
17、要体现在界面斜率的计算。 网格（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,
18、1*0,1，初始界面为圆，圆心为(0.5,0.3)、半径为0.2。计算网格为200*200时间步长为0.0001s。 0000,cossin,sincosux yxxyyvx yxxyy 图3.3 剪切速度场数值结果（上为旋转1s，下为旋转2s）六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法PLIC界面重构技术在二维直角坐标中的实现 PLIC界面重构技术也是应用直线段对界面进行重构，但其与FLAIR不
19、同，PLIC方法是在单个网格内对界面进行近似。PLIC界面存在20种情况，其中4种情况较为简单，界面与坐标轴平行，其界面流量的计算非常简单。其它16种情况的界面与两个坐标轴均有交点，此界面的法向为：22,/xxyxi ji ji jnnnn,xyi jjnnn,1,11,1,11,11,1,122/xi jijijijijijijnffffffx ,1,1,11,11,1,11,122/yi jiji jijiji jijnffffffy 则单位法向为22,/yxyyi ji ji jnnnn分母上系数，二维时为8；三维时为32？0,0 xynn 通过坐标轴的旋转以及坐标转换可将各种界面转换为
20、 的情况， 的界面存在4种情况，分别为：0,0 xynn图3.4 PLIC最终简化类型 上述4种情况是以图3.5中的A、B两点为分界的。图3.5中界面的方程为：xy0dydxnnxnyx1y1BA图3.5 PLIC界面示意图1xxyynnyxxnn 当 时，只存在图3.4中的1、2、4情况，此时界面沿法线方向移动会先后经过A、B点。界面经过A点时，阴影面积为：界面经过B点时，阴影面积为：xyndxdyn2112xynSdxn22122yxxynnSdydydxnn方程？若 ，则界面为情况1；若 ，则界面为情况2；若 ，则界面为情况4。 1CdxdyS12SCdxdyS2SCdxdy 当 时，只
21、存在图3.4中的1、3、4情况，此时界面沿法线方向移动会先后经过B、A点。界面经过B点时，阴影面积为：界面经过A点时，阴影面积为：xyndxdyn2112yxnSdyn22122yxyxnnSdxdxdynn若 ，则界面为情况1；若 ，则界面为情况2；若 ，则界面为情况4。 1CdxdyS12SCdxdyS2SCdxdyC=?相函数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。 xy0dydxnnxnyx1y1BA图3.6 PLIC界面示意图22
22、211112yxyxnnFsxM xdxN xdyCdxdynn 通过上式可以求得各情况时，界面与x轴的交点为： 第一种情况： 第二种情况： 第三种情况： 第四种情况： 12.0yxndxxCdyn12.0yxndyxCdxn12.0 1.0yyxxnnxdxdyCdxdynn12.0yxnxCdxdyn,10.0i jf 现以情况4为例介绍目标网格向周围4个网格的流量输运。xy0vi,j+1vi,jui+1,jui,jx1y1图3.6 PLIC界面的第4种情况(1)网格i,j向i,j+1网格的流量输运：如果Vi,j+10，则存在输运： IF Vi,j+1*dt0，不存在输运； 如果Vi,jy
23、1 else,1i jfC2,1,0.5*1*/*yi ji jxnfCyvdtdx dyn （3）网格i,j向i+1,j网格的流量输运： 如果Ui+1,j0，则存在输运： IF Ui+1,j*dt0，不存在输运； 如果Ui,jx1 else1.ijfC21,0.5*1*/*xiji jynfCxudtdx dyn 其它几种情况的界面输运与情况4相近，求得各界面的输运以后就可以得到新时刻的相函数分布。 PLIC算法的验证算例如下： （1）剪切流场0000,cossin,sincosux yxxyyvx yxxyy 取(x0,y0)=(0.5,0.5)，计算区域取0,1*0,1，初始界面为圆，圆
24、心为(0.5,0.3)、半径为0.2。计算网格为200*200时间步长为0.0001s。 图3.7 剪切速度场数值结果（上为旋转1s，下为旋转2s） （2）旋转流场：00,ux yyxvx yxy 取(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，Y
25、L=0.2m；计算所采用的网格为102*102；计算时间步长为0.0001s；初始气泡的半径为0.02m；圆心为（0.1,0.1）。气水两相的物性分别为：水的密度为1000.0kg/m3、水的粘度为0.001Pa.s；空气的密度为1.266kg/m3、空气的粘度为17.8E-6Pa.s；气水表面张力为0.0728N/m。界面重构方法图3.9 气泡上升的数值结果（依次0ms、50ms、100ms、130ms、170ms、200ms）图3.10 气泡上升的流场（依次100ms200ms） （2）PLIC方法在油水两相流动中的应用r界面锥 油水混合区不渗透固壁油水相界面计算区域(a) 油洞抽象示意图
26、thLd(b) 物理模型及尺寸示意图图3.11 溶洞油藏的物理模型与计算模型六气液相界面的构造方法 3 3、 PLICPLIC界面重构技术在三维直角坐标中的实现界面重构技术在三维直角坐标中的实现 3 3、 PLICPLIC界面重构技术在三维直角坐标中的实现界面重构技术在三维直角坐标中的实现 在三维网格中同样存在满网格、空网格、相网格，空网格不需要计算界面的流量输出，而满网格的界面输运计算非常简单，工作的重点为相网格。 根据相网格中界面的法向可将界面分为以下三种： （1） ，即界面与网格的一个面一个面平行，此界面的输运比较好计算； （2） ，即界面与网格的一个坐标轴一个坐标轴平行，此界面的输运与
27、二维PLIC的界面输运相近，因此求解并不是很困难； （3） ，界面与所在的平面与三个坐标轴三个坐标轴均有交点，此界面存在情况较多，求解也较为困难。,0,0nC nn1,2,0nC nCn1,2,3nC nCnC 第一种类型界面的输运求解非常简单； 第二种类型界面与二维PLIC相近；因此下面将不对这两种类型的界面进行分析。yzxnx0,ny=0,nz=0yzxnx0,ny0,nz=0yzxnx0,ny0,nz0图3.16 三种界面类型的示意图 的界面情况非常多，因此需要简化。首先可通过坐标轴旋转、坐标转换及目标流体与非目标流体的转化将各种情况转为 ，此时可将界面分为5种情况，如下：1,2,3nC
28、 nCnC0zyxnnnyzxyzxyzx情况1情况2情况3yzxyzx情况4情况5图3.17 3D-PLIC的5种界面情况 界面方程为 式中， 上述5种情况的判断准则为： （1）当界面为第一种情况时，应满足下式： （2）当界面为第二种情况时，应满足下式： （3）当界面为第三种情况时，应满足下式：312316n n n Cn33311232216nn n n Cnnn333331232211233313233333123221123121266if nnnnnnn n n Cnnnnnif nnnnnnn n n Cnnnn1 12233nnndiiixx （4）当界面为第四种情况时，应满足下
29、式： （5）当界面为第五种情况时，应满足下式：1233333333333132123121231123.6nnnandnnnnnn n n Cnnnnnnnnn1233331231212.6nnnandn n n Cnnnn 以情况1为例来分析目标网格向周围6个的流量输运 目标流体所占的体积为： 因此： 所以可求得界面在x轴上的截距为：yzxABCEFGKPMQN图3.18 界面的第一种情况 2316xABCEyznVAEn n2316xyznAECdxdydzn n326yzxn nAECdxdydzn 相网格向周围6个网格的输运如下：(1)网格i,j,k通过PQGF面向i+1,j,k网格的
30、流量输运：如果Ui+1,j,k0，则存在输运： IF Ui+1,j,k*dt0，不存在输运；如果Ui,j,k=AE else 1, ,ij kfC231, , ,1/6xij ki j kyznfCAEudtdxdydzn n(3)网格i,j,k通过BMKG面向i,j+1,k网格的流量输运：如果Vi,j+1,k0，则存在输运： IF Vi,j+1,k*dt0，不存在输运；如果Vi,j,k=BE else ,1,i jkfC23,1, ,1/6yi jki j kxznfCBEvdtdxdydzn n(5)网格i,j,k通过PNMQ面向i,j,k+1网格的流量输运：如果Wi,j,k+10，则存在
31、输运： IF Wi,j,k+1*dt0，不存在输运；如果Wi,j,k=CE else , ,1i j kfC23, ,1, ,1/6zi j ki j kxynfCCEwdtdxdydzn n 其它几种情况的界面输运与情况1相近，求得各界面的输运以后就可以得到新时刻的相函数分布。 PLIC界面重构技术的验证算例： （1）常数速度场 算例为立方体在上述流场中的变化，三维计算区域为 ，所计算立方体的下角点为 、边长为0.4，计算网格分别为 、 ，时间步长为0.001s。, ,0.5,0.5,0.5u v w 0,10,10,10.06,0.06,0.0652 52 52102 102 102图3.
32、19 常数速度场的验证（网格为52*52*52、时间分别为0.0s、0.5s、1.0s） 图3.20 常数速度场的验证（网格为102*102*102、时间分别为0.0s、0.5s、1.0s） （2）旋转速度场 考察了Zalesaks球在流场中的变化，此计算区域为 ，Zalesaks球的球心为 ，半径为0.15， 计算网格分别为 、 ，计算时间步长为0.001s。 0.50.50.0uyvxw 0,10,10,10.5,0.75,0.5102 102 102152 152 152 图3.21 旋转速度场的验证（网格为102*102*102、时间分别为0.0s、1.0s、2.0s） 图3.22 旋
33、转速度场的验证（网格为152*152*152、时间分别为0.0s、1.0s、2.0s） （3）剪切速度场 算例为圆球在上述流场中的变化，三维计算区域为 ，所计算球体的求心为 、半径为0.15，计算周期为T=3s，计算网格为 、 ，计算时间步长为0.001s。 2222sinsin 2sin 2cossin 2sinsin 2cossin 2sin 2sincostuxyzTtuxyzTtuxyzT 0,10,10,10.35,0.35,0.35127 127 127152 152 152图3.23 剪切速度场的验证（网格为127*127*127、时间分别为0s、0.5s、1.0s、1.5s、3
34、.0s、2.5s、2.0s ） 图3.24 剪切速度场的验证（网格为152*152*152、时间分别为0s、0.5s、1.0s、1.5s、3.0s、2.5s、2.0s ） 3D-PLIC技术在气泡动力学中的应用： （1）单个气泡的上升过程t=0mst=10mst=20mst=30mst=35mst=40ms图3.25 单个气泡的上升过程 （2）两个气泡的融合过程t=0mst=20mst=40mst=60mst=80ms图3.26 两个气泡的融合过程六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构
35、造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法六气液相界面的构造方法七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较七数值模拟算例与比较八关于相界面的构造问题的几点说明八关于相界面的构造问题的几点说明九VOF方法中气、液两相流场中的物性表示方法十VOF方法的优点十一VOF方法的缺点十一VOF方法的缺点十二VOF方
36、法的应用VOF方法在气泡动力学中的应用：图1.1 VOF方法在气泡动力学中的应用 （M.van Sint Annaland 2004）VOF方法模拟液滴撞击薄液膜：图1.2 VOF方法模拟液滴撞击液膜 （N.Nikolopoulos 2007）VOF方法对管内弹状流的模拟：图1.3 VOF方法模拟垂直管与倾斜管内弹状流 （Taha Taha 2005）十二VOF方法的应用更多的应用在实践中体会！十三结束语气液两相流的数值模拟方法有很多，本文根据所使用的数学物理原理的不同对相关的方法进行了分类和介绍。总而言之，实现气液两相流数值模拟的关键是相界面位置实现气液两相流数值模拟的关键是相界面位置及其运
37、动过程的确定。及其运动过程的确定。关于相界面描述方法的研究也经历了由简单到复杂、循序改进的发展过程；新方法总是在原有方法的完善或启迪下形成的。十三结束语值得一提的是，数值模拟毕竟不是物理过程本身。一个完整的数值模拟过程实际包含着从物理问题到数学模型，从控制方程（数学模型）到离散方程，再到离散方程的求解等环节；这些环节会涉及数学模型的准确性、离散方程的截断误差、离散方程解的离散误差及数值计算过程的舍入误差等误差；因此，要获得能反映所描述现象的物理本质的数值解是一件非常复杂而艰巨的工作，对气液两相流的数值模拟更是如此。十三结束语为了实现界面（包括气液两相流中的气液界面）的捕捉或描述，数学家和计算物
38、理学者设计和开发了许多中各具特色的格式（或方法）。应该看到的是，数学家数学家追求的效果是被描述的界面尽可能接近物理上追求的效果是被描述的界面尽可能接近物理上的间断本质（阶跃式变化，或陡峭变化）的间断本质（阶跃式变化，或陡峭变化）；但是在实际气液两相流的数值模拟中，这样的高阶格式不仅计算量大，而且很容易引发数值不稳定性。正如早些时候为了实现对激波的求解而有意识地引入“人工粘性”一样，只要能达到计算要求，不一定要追求高阶的数值格式。只要能达到计算要求，不一定要追求高阶的数值格式。十三结束语基于上述认识，我们不难理解，每一种方法都有自己的优点和缺点，就像每一种实验方法也有各自的优缺点一样。就气液两相
39、流这一复杂就气液两相流这一复杂的物理现象而言，大多数的数值模拟方法仍在进一步的完善和检验之中，的物理现象而言，大多数的数值模拟方法仍在进一步的完善和检验之中，要说那一种方法是完善的或最好的尚为时太早。要说那一种方法是完善的或最好的尚为时太早。 Level Set方法和VOF方法是目前较为流行的两种方法，具有很好的应用前景。而由于其自身的特点，Level set方法具有更大的优越性。新的尝试是将VOF方法和Level Set方法相结合，形成所谓的VOSET方法。 十三结束语各类方法共同的共同的不足不足之处之处表现在下述几个方面第一、两相界面两侧的物性差别不能太大，不能满足实际的气液两相流问题的要求（例如，气液密度比为1:1000）；第二，不能精确的求解气液固三相交接处的移动界面问题；第三，对于气液两相流这样的问题，当气液界面上存在能量与质量传输时，如何准确地确定界面的物理与几何特性变化仍值得深入研究。这些都是现有的数值模拟方法值得改进的地方，也是新的数值模拟方法应着力解决的问题。第九讲结束！
