VOF方法中相函数的控制方程的求解方法

VOF方法中相函数的控制方程的求解方法一、VOF方法的基本原理VOF方法的基本原理是通过研究网格单元中流体和网格体积比函数F来构造和追踪自由界面的，而非追踪自由液面上质点的运动。VOF方法是在整个流场中定义一个函数F，在每个网格中，这个函数定义为一种流体（目标流体）的体积与网格体积的比值。只要求得每个网格上的F值，就可以构造出运动界面，然后求解物理方程时，可以在界面附近作特殊的精细处理，提高分辨率和精度。由于每个网格只需保存一个F值，从而可以大大减少储存量。设计算区域为，流体A所在区域记为A，流体B所在的区域记为B，首先定义一个函数(,)xt，1,(,)0,ABxxtx(1)对于两种不相容流体组成的流场，(,)xt满足拉格朗日流体体积的特性：0uutxy(2)其中(,)Vuv是流体的速度场，ijV为单个网格的体积，在每个网格ijI上定义ijF为(,)xt在网格上的积分，1(,)ijijijIFxtdVV(3)称之为VOF函数，同样它满足：0FFFuutxy(4)称之为VOF(VolumeofFluid)方程。由上可以看出每个流体上的体积函数实质上是F网格中的流体体积网格体积显然，F=1的网格充满流体A，我们称之为流体网格(fluid)；当F=0时为空网格(void)；而0F1的网格，则是含有流体界面的网格，称为界面网格(surface)。VOF方法是通过求解流体体积函数，实现对运动界面的追踪。由于流体体1积函数是单个网格中一种流体体积与网格体积之比，是具有物理意义的。但在求解VOF方程时，要保持它的这一物理意义是困难的。首先要求求解格式是守恒的，而且每一时间步得到的函数值在跨过界面时，即在含有界面的网格上，VOF函数要0F1，在其它网格上必须是1或0。二、流体体积方程的求解方法Hirt和Nichols针对界面追踪的一般方法以及通过装配格子内界面的斜率来确定运动界面的斜率和位置，提出了VOF方法。在这以后，AShgriz和Poo,Kim等人通过对格子内界面斜率更加精确的构造得到了其它不同的重构方法，尽管他们的精确性和形式是不同的。Rider和Kothe,Gueyffier和Lijie,Youngs等人利用了运动界面的法向矢量近似值来构造界面。Rudman,Vincent和Caltagirone,Ubbink和Issa等人采用了更精细的方法求解流体体积函数方程。对于流体体积方程的求解归纳起来主要有两种方法：一种是基于几何学的求解方法，这种方法又可以成为重构技术；另一种则是基于数学的求解方法。1.Hirt&Nichols法：i-1ii+1j-1jj+1图1Hirt-Nichols的重构模板这是一种基于几何学的求解方法，该方法用水平和竖直两种直线近似一个网格内的流体界面。采用施主--受主(Donor-Acceptor)的差分逼近格式计算运动界面的流体体积函数方程。界面的重构如下：将流体的自由表面看作局部的单值函Y(x)和X(x)采用9个网格的模板(如图1所示)，计算i-1,i,i+1网格的列Yl值和j-1,j,j+1网格列的Xl值，估算出每个网格上界面的斜率值dY/dx和dX/dy，然后根据流体体积函数和斜率的大小确定网格(i,j)上的自由面的位置和方向。11,1,,1jllkkkjYFyliii(5)211112()2iiiiiidYYYdxxxx(6)11,1,,1ilklkkiYFxljjj(7)11112()2jjjjjjXXdXdyyyy(8)如果：dYdXdXdY；则定义网格内的界面是水平的；否则，网格内的界面是竖直的。2.FLAIR法：这种求解流体体积方程的方法也是基于几何学理论。1991年，Ashgriz和Poo的FLAIR(FluxLine-segmentmodelforAdvectionandInterfaceReconstruction)方法通过构造跨过任意网格边界的两个相邻网格的带有倾角的直线段来近似界面。由于界面构造涉及到两个网格，它分成多种情况分别进行计算。首先利用网格边界的流体速度，确定施主单元和受主单元，分别用FD,FA表示其流体体积函数，然后按以下几种情况计算：（1）施主网格是满网格，即FD=1.0;（2）施主网格是空网格，即FD=0.0;（3）施主网格和受主网格都是半网格，即0.0FD1.0，0.0FA1.0;（4）施主网格是半网格，受主网格为空或者满网格，即0.0FD1.0，FA=0.0，或FA=1.0;对于前两种情况，单位时间内，流体体积流通量分别为ut，0.0。第三种情况又可分为十六中情形，经过对称和反转可以分为四类(如图2所示)。此时abff，且目标流体为图中阴影部分。近似界面直线的斜率可以通过两边的流体体积函数值来计算。从而可以计算出网格边界的流量运输，即ut的长度上阴影部分的面积(左边为施主网格)。3afbfut(I)(II)(III)(IV)afbfafbfafbf图2FLAIR的界面重构3.Youngs法：Youngs的PLIC(PiecewiseLinerInterfaceCalculation)方法，这种方法也是基于几何学原理的，是在单个网格上用直线段近似界面。直线的斜率由网格内界面的法向量确定，,,(,)xyijijnnn(9),1,11,1,11,11,1,1(22)/xijijijijijijijnFFFFFFx(10),1,1,11,11,1,11,1(22)/yijijijijijijijnFFFFFFx(11)与FLAIR方法类似地有十六种情形，通过对称和翻转可分为如图3所示的四种类型。运动界面与X轴的夹角，将其规格化为，则利用夹角和网格内的体积函数确定属于哪种类型，从而构造界面，再计算在一个时间步长内与四个相邻网格的体积流通量，来修改它们的流体体积函数值。arctan(),arctan(tan),2xynxny(12)(IV)(III)(II)(I)图3FLIC的四种类型44.积分平均型TVD格式这种求解流体体积方程的方法是基于数学理论，直接对流体体积方程积分求解。积分平均TVD格式是运动界面的捕捉的方法之一，它把求解函数构造成分片的线性函数，用特征线方法（式15）推进时间步，并利用Superbee限制器对斜率（式16）加以限制，使格式保持单调，具有TVD性质。1111221(,)11((),),02nnnnttiitinntFFxtdttFxxVtttdtVt(15),1()/LWnnniiiSFFx(16)对于二维和三维情形，我们可以采用分数步法进行计算，即在前半个时间步计算102FFutx，得到中间值1/2nF，再以此为初值在后半个时间步内计算102FFvty，得到n+1时刻的值1nF。对于求解函数的构造格式，还可以构造分片的二次或高次多项式近似（例如，PPM方法或ENO和WeightedENO方法），时间步的推进用Taylor展开方法，斜率的构造采用其它方法。5.FCT-VOF方法：FCT-VOF方法是基于FCT(flux-correctedtransportalgorithms)的思想求解流体体积函数方程的。该方法采用数值通量守恒形式（17），数值通量采用低阶通量LF（18）为upwind格式，和高阶通量HF（19）为downwind格式，两者相结合构造反扩散流量对中间值进行修正。11/21/2(),nniiiitFFrFFrx(17)1/21/21/21/211/200LiiiiiiiuFuFuFu(18)1/211/21/21/21/200HiiiiiiiuFuFuFu(19)6.Ubbink的界面捕捉方法--CICSAM方法：Ubbink和Issa提出的CICSAM(CompressiveInterfaceCapturingSchemeforArbitraryMeshes)突破了以往重构方法通过方向分裂发法向多维的推广，而是根据Crank-Nicolson格式离散（式13和14），采用半隐格式处理推广到多维，这5也为该方法应用到无结构网格中扫清障碍。11()(()())2nttttttpppIIIIfFFVFFFFt(13)(1)IIDIAFFF(14)式中，p,I,n,pV,IF分别代表控制体的中心，网格的边，网格的边界数目，控制体的体积，第I边上向外的通量。其中I通过NVD方法确定。它是在Hirt和Nichols方法根据自由面方向决定使用耗散格式还是色散格式的思想的基础上，提高耗散和色散的精度使得两者的转换更加的光滑。CICSAM方法通过修改网格边界上的体积流通量，既保持了体积函数有界性，又维持了体积函数的守恒。三、结束语本文对VOF方法中的几种主要的求解流体体积函数的方法分别从理论上做了简要的介绍。这些求解方法理论上的难易程度直接影响编程的难易程度，其中基于几何学理论的三种求解方法：Hirt&Nichols、Flair、Youngs，从在单个网格内用直线段近似界面到在相邻网格内用斜线段近似界面以及在单个网格内用斜线段近似界面，理论上越来越复杂，但求解精度越来越高。参考文献[1]C.W.Hirt,B.D.Nichols.VolumeofFluid(VOF)MethodfortheDynamicsofFreeBoundaries[J].ComputationalPhysics,1981,39:201-225.[2]孔珑,等.两相流体力学.北京：高等教育出版社,2004.[3]张健,方杰,范波芹.VOF方法理论与应用综述[J].水利水电科技进展,2005(4).[4]端木玉,朱仁庆.流体体积方程的求解方法[J].江苏科技大学学报(自然科学版),2007.[5]陈萍萍.VOF运动界面重构的流体体积分数保持法.河海大学硕士学位论文.2007.