技术领域
[0001] 本
发明涉及自适应非结构化三
角形网格生成及两相界面构造方法，尤其涉及一种高精度的两相流体界面捕获方法。
背景技术
[0002] 两相或
多相流动广泛存在于自然界及工程中，随着计算机技术的快速发展，数值模拟已经成为研究多相流动的有效方法。在研究气液/液液两相流动时，由于存在运动的相界面，同时在自由界面上流体物性和流动状态可能发生剧变，所以考虑相界面的
变形和位移以及界面的处理就显得至关重要，能否精准捕捉两相界面关系到整个流场计算的正确性和合理性。VOF(Volume ofFluid)方法作为一种界面捕捉类方法，由于其较好地满足
质量守恒特性该方法已经被广泛地应用于多相流体数值模拟。现有的关于相界面捕获的VOF方法研究大多是基于结构化网格上的，然而实际问题中经常遇到复杂的几何区域，由于结构化网格的适用范围有限，当计算区域不规则的时候，非结构化网格可以适应各种形状的计算区域，发展基于非结构化网格的VOF方法，可以有效模拟复杂计算区域中的多相流动，比结构化网格具有更加普遍的适用性，因此发展基于非结构化网格的界面捕获
算法非常重要。
[0003]
发明人前期研究工作中开展了基于非结构化网格的SLIC-VOF方法(可参考HUANG M.，CHEN B.，WU L.L.A SLIC-VOF method based on unstructuredgrid.Microgravity Science and Technology，2010，22(33)：305-314)，发现对于基于非结构化网格上的多相流动界面捕获方法，网格的质量以及网格对于界面的自适应性会直接影响计算结果的精度。尽管气泡堆积法生成的高质量非结构化网格使得基于非结构化网格上的VOF方法计算的精度得到了显著的改进，但由于相界面构造算法本身的精度问题，捕获的相界面不够光滑。为了提高所捕获界面的精度构建较为光滑的相界面，可以通过以下两方面的改进：1)增加计算区域内网格数量；2)发展高精度的界面构造方法。对于网格单元数量的增加，如果是增加整个计算区域内的网格数量，则需要额外的引入大量的计算资源，对于一般的计算机可能难以满足计算的要求，同时大量增加计算网格单元势必会降低程序的执行效率，然而在计算中我们最关心是两相流体的相界面
位置区域，因此考虑仅仅将相界面处的网格单元进行加密，有必要发展一种适用于两相流的自适应非结构化网格生成方法。当前国内外学者已经开展关于自适应网格生成方法的研究并提出了相关的实施方法，但大多数是基于结构化网格上的算法，对于非结构化网格上的适用于两相或多相流动的自适应网格生成技术刚刚开始初步研究，如Ito等人在2010年提出了一种适用于
气液两相流体的自适应非结构化网格技术，但在其提出的方法中需要反复地从细分前的双亲网格单元到细分后的子网格单元的对相函数等参数的插值分配，增加了算法的计算复杂性。
[0004] 对于VOF高精度界面构造方法，国内外也有很多学者作出了自己的尝试，但是目前基于非结构化网格的界面重构方法仍然不多，研究人员主要针对两个方面进行发展和改进：界面的构造和流体的输运。界面的构造是计算流体输运和
可视化的前提，现今较常见的是用一条具有任意斜率的线段逼近真实界面线，Barth在1995年提出了利用最小二乘法
迭代求出线段位置的方法，得到了较广泛的应用，但是这样每次构造界面都要在网格上进行数次迭代计算，而精度提高有限；对于流体输运量的计算，鉴于非结构化网格的不规则性，很多学者避开了流量的欧拉格式计算而采用Mosso在1997年提出的欧拉-拉格朗日结合的方法计算流量，该方法可以达到较高的精度，但是要用到欧拉和拉格朗日两种格式，此外还涉及到复杂的几何计算。
发明内容
[0005] 本发明的目的在于针对
现有技术的不足，提供了一种高精度的两相流体界面捕获方法。
[0006] 为达到上述目的，本发明采用的技术方案是：
[0007] 1)对计算区域进行建模，生成初始的高质量非结构化三角形网格，根据计算条件，对不同的流体相初始化VOF相函数F值，同时设置计算工况参数；
[0008] 2)采用“预估-校正”的两步计算实现网格的自适应加密：首先计算VOF相函数F在流场中的传输预估下一时刻流体相界面的位置，进而对F∈(0，1)的相界面网格进行细分加密将高质量非结构化三角形网格细分为四个子网格，然后再在细分加密后的网格上重新计算得到在细分后网格上下一时刻相函数F的分布；
[0009] 3)完成上述计算后，
对流体全部流出或者全部充满
指定相流体的四个子网格重新合并为一个网格单元，并计算该网格的相函数F值；
[0010] 4)采用PLIC方法重新构造流体相界面；
[0011] 5)反复迭代执行步骤(2)～(4)，直至完成计算终止。
[0012] 所述的步骤1)包括对计算区域建模生成初始的高质量非结构化三角形网格，然后根据计算条件对不同的流体相初始化VOF相函数F值，F表示流体某相在网格单元内的体积与网格体积的比，其取值介于0～1之间，F为0代表空网格，为1代表满网格，介于0和1之间则代表存在相界面，同时给定计算的工况参数，包括不同相流体速度的给定。
[0013] 所述的步骤2)采用“预估-校正”两步计算实现网格的自适应加密其具体步骤为：
[0014] a)预估相界面位置并加密：采用VOF方法计算流体相函数F的分布，预估下一时刻流体相界面的位置，进而对相函数F值介于0和1之间的相界面网格单元进行细分加密；
[0015] b)校正计算相函数F的传输：重新赋予网格单元预估计算前的F值，同时对细分加密后新添加的子网格初始化相函数F值，进而在细分后的子网格上重新计算流体相函数F在每个网格单元中的分布。
[0016] 所述的步骤3)完成在细分后网格单元上的计算后，将指定相流体全部流出(即F＝0)或者全部充满指定相流体(即F＝1)的四个子网格单元重新合并为细分前的网格单元，并得到该网格的相函数F值。
[0017] 所述的步骤4)采用如下方案构造相界面：
[0018] a)网格中分段线性界面的构造先根据周边相邻网格中F值分别求解出共用当前网格的三个
顶点的所有网格单元在该顶点处流体体积比函数的均值Fi，再用这三个顶点上F的均值Fi求其对于当前网格中心点相函数F0的梯度值，进而确定该网格界面线的法线在x轴和y轴上的分量nx和ny如下
[0019]
[0020] 其中(x0，y0)为三角形网格中心点的坐标，(xi，yi)为三角形网格顶点的坐标，接下来再由该网格单元的F值大小确定界面线的位置，从而构造出一条具有任意斜率和任意位置的界面线来逼近穿过当前网格的真实流体自由面的位置，完成界面的精细重构；
[0021] b)流量输运依靠几何关系进行计算，对于自由面上网格来说，根据每一边与界面线的位置关系以及该边沿外法向的速度分量，计算通过该边的流体体积输运量；对于流体内部的网格来说，则直接计算每条边上的流体输运，对所有网格进行遍历从而得到下一时刻流体体积比函数F分布，据此进行下一轮的界面构造。
[0022] 本发明针对界面构造和流体输运分别提出了解决方案，界面线的位置依靠网格定点平均的流体体积比函数在中心处的梯度矢量加权平均值求出的具有任意斜率线段来近似，流体输运依靠纯欧拉流量计算通过网格和界面线的几何关系来确定，最终实现相界面的重构。本发明拟在气泡堆积法生成的初始非结构化网格
基础上，根据VOF方法相函数的特征对两相流体的相界面网格进行自适应加密，提高界面捕获的精度得到光滑的相界面。
附图说明
[0024] 图2是本发明流体相界面处网格单元的细分示意图；
[0025] 图3是本发明相邻网格单元的边界处理；
[0026] 图4是本发明相邻网格均被细分的特殊情形；
[0027] 图5是本发明自适应网格VOF算法实施步骤；
[0028] 图6是本发明细分前后相界面重构的对比；
[0029] 图7是PLIC方法中界面法向的确定示意图；
[0030] 图8是PLIC方法中网格间流体输运的计算图；
[0031] 图9是PLIC-VOF方法实施流程图；
[0032] 图10是计算区域流体分布示意图；
[0033] 图11是PLIC-VOF计算结果图。
具体实施方式
[0034] 下面结合附图对本发明作进一步详细说明。
[0035] 如图1所示，其包括以下步骤：
[0036] 1)对计算区域进行建模，生成初始的高质量非结构化三角形网格；
[0037] 2)根据计算条件，对不同的流体相初始化VOF相函数F值，同时给定流体运动速度场；
[0038] 3)采用“预估-校正”的两步计算实现网格的自适应加密：首先计算VOF相函数F在流场中的传输预估下一时刻相界面的位置，进而对F∈(0，1)的相界面网格进行细分加密，然后再在细分后的网格上对预估计算前的流体状态重新计算，得到在细分后网格上下一时刻相函数F的分布；
[0039] 4)完成上述计算后，对指定相流体全部流出(即F＝0)或者全部充满指定相流体(即F＝1)的四个子网格单元重新合并为细分前的网格单元，并得到该网格的相函数F值；
[0040] 5)采用PLIC(Piecewise Linear Interface Capturing)方法重新构造界面；
[0041] 6)反复迭代执行步骤3)～5)，直至完成计算终止。
[0042] 所述的步骤1)对计算区域建模进而采用气泡堆积法生成初始的高质量非结构化网格(可参考WU L.L.，CHEN B.，ZHOU G.L.An improved bubblepacking method for unstructured grid generation with application tocomputational fluid dynamics.Numerical Heat Transfer，Part B，2010，58(5)：343-369.)。
[0043] 完成步骤1)后，如步骤2)所述根据计算条件对不同的流体相初始化VOF相函数F值，为了得到F值，首先在整个离散的计算区域内定义VOF函数f，对于某一相取为1，另一相取为0。f在单个网格单元内的积分与网格面积的比记为相函数F，如公式(1)所示，其取值介于0～1之间。
[0044]
[0045] F为0代表空网格，为1代表满网格，介于0和1之间则代表存在相界面。在计算前要尽可能地保证初始流体相界面的光滑，以减少由于初始相界面的不光滑性对后续计算的影响，于是要根据给定的工况参数计算初始的相函数分布；同时在该步骤中还需设置计算条件，包括不同相流体速度的给定。
[0046] 接下来，如步骤3)所述采用“预估-校正”的两步计算实现网格的自适应加密，首先计算VOF相函数F在流场中的传输预估下一时刻相界面的位置，对这些位置处的网格单元进行自适应细分加密，然后再在细分后的网格上重新计算下一时刻相函数F的分布。
[0047] 对界面位置处网格单元的细分准则是：根据VOF方法原理可知，相界面网格单元的相函数F值介于0和1之间，而对于两相流体我们最关心的区域也是相界面所处的位置，因此可根据相函数F值对三角形网格单元进行自适应细分。对于图2(a)所示的三角形网格单元，由于该单元存在两种不同流体，该网格单元的相函数F∈(0，1)，可以对该三角形网格单元进行细分。具体做法是联接该三角形网格单元各条边的中点，即可将该网格单元细分为四个小的网格单元，如图2(b)所示。细分操作可以一直进行下去，直至满足计算精度需求为止。
[0048] 对于执行了细分操作的网格单元ABC，原有的边界BC被一分为二，与其共用同一边界的网格单元BCD会出现一个边界对应细分后两个边界的问题，如图3(a)所示，计算时需要特殊处理，这样势必会破坏算法的一致性，本发明对细分单元的邻接网格单元进行了特殊处理。如图3(b)所示，假设网格单元ΔABC经过细分成为四个子单元，对于其共用BC边的相邻三角形单元ΔBCD，可以联接BC边的中点E和顶点D将ΔBCD分为两个子单元，这样可以保证相邻网格间界面一一对应的关系，进而保证了原始算法的一致性避免了在细分后单元界面BC上计算的特殊处理。如果某个网格单元有两个或者两个以上的邻接单元被一分为四，如图4(a)所示ΔBCD，在对该网格单元界面进行上述边界处理时会出现CF边和DE边的重叠。对于这种情形，本发明中采用了将该单元进一步细分为四个子网格单元，如图4(b)所示。同时，由于上述边界处理将一个网格单元细分为两个网格单元后势必会造成网格单元质量的降低，为了将该连接单元对界面的影响尽量减小，我们在计算中可以将相界面处网格单元的细分向周围邻接网格单元延续扩展几层。
[0049] 在自适应算法实施过程中，每一时刻流体会由一个网格单元进入其相邻网格单元中，该过程伴随着网格的进一步细分和合并。进一步细分是指网格单元当前时刻没有指定相流体(F＝0)或者充满指定相流体(F＝1)，随着时间的推移有流体的流进或流出，所以在下一时刻相界面将位于该网格单元内，于是将该网格单元细分为四个小的子网格单元；合并是指经过一个时间步长Δt的计算后，原来细分为四个子网格单元内不存在或者全部充满了指定相流体，此时可将这四个子网格单元再合并为细分前的网格单元。同时在网格细分、合并过程中需要分配相函数F，现有的文献中常用的做法是通过插值计算得到细分后小网格上的F或者是合并后网格上的F，而本发明中为了算法实施的简便，采用了“预估-校正”的思想，即先预估流体相界面在下一时刻的位置对其细分加密，接着再在细分后网格单元上重新计算相函数F进而重构界面，从t时刻到t+Δt的具体实施流程如下：
[0050] (1)由初始的网格和对应的F值，通过VOF方法计算下一时刻的F值分布，预估下一时刻流体相界面的位置，将新的相界面网格单元细分为四个子网格单元，如图5(a)和(b)所示；
[0051] (2)将当前时刻初始时的相函数F值在细分后的网格单元上重新分布，对于预估下一时刻流体相界面将要流到的位置处的网格单元在细分后，由于当前相界面不在这些网格单元内，所以它们的相函数F是0或者1，于是细分后的小网格单元的F值也同样是0或者1，如图5(c)所示；
[0052] (3)接下来再在细分后的子网格上重新计算流体相函数F在每个网格单元中的分布，完成上述计算后，将指定相流体全部流出(即F＝0)或者全部充满指定相流体(即F＝1)的四个子网格单元重新合并为细分前的网格单元，并得到该网格下一时刻的相函数F值
[0053] (4)采用PLIC方法重新构造界面，如图5(d)所示。
[0054] 由于本发明提出的先预估再计算的思想，所以每次被细分的网格单元在下一时刻将会是相界面网格单元，而在当前时刻不是相界面所处的位置，于是这些网格单元的相函数F值为0或者1，将其细分后新添加的子网格单元的相函数F值同样为0或者1，这样可以有效地避免从细分前网格分配F到细分后子网格单元时复杂的插值计算。对于流体相界面处网格单元，细分前后相界面的重构与实际相界面的对比如图6(a)和6(b)所示。
[0055] 对于界面构造和流体输运计算，本发明提出了基于非结构化网格的PLIC(Piecewise Linear Interface Capturing)方法，相比于SLIC方法更加精细，可以得到很高的界面捕捉精度。该方法满足VOF的基本方程
[0056]
[0057] 在三角形网格上离散得到如下形式
[0058]
[0059] 界面重构的核心是对式(3)的求解，实施起来主要分为以下步骤：
[0060] (1)确定流场中初始时刻的速度场以及流体体积比函数F值的分布；
[0061] (2)在位于自由面上的网格中构造一条直线段来近似实际的界面线；
[0062] (3)结合流场中流体体积比函数F值的分布和网格中构造的界面线位置计算当前时刻到下一时刻网格间的流体输运，从而确定下一时刻F值分布；
[0063] (4)在需要的时刻对界面进行可视化的输出。
[0064] 具体需要着重注意的是第(2)和第(3)步。在第(2)步中，界面线的构造依靠求解界面上网格与其周边网格中的流体体积比函数F在当前网格中心点处的梯度值来确定。如图7所示，假设网格ABC位于自由面上，对该网格单元每个顶点所属的全部网格中的F值在该点求均值，再由三个顶点的F均值在网格ABC的中心处求梯度，则该梯度矢量的方向就是界面线在该网格内的法线方向，其在x轴和y轴上的分量nx和ny如下：
[0065]
[0066] 其中，i为网格单元ABC的顶点编号，下标0为ΔABC中心点处的值，Fi由下式给出
[0067] N＝m，n，k (5)
[0068] 已知网格ABC内的F值大小，再根据界面线的法向，即可通过几何方式得到网格中界面线的具体位置，即界面线与网格的两个交点坐标。至此，自由面网格中的界面线构造就完成了，对流场中所有的相界面网格进行上述计算，就得到了整个流场中的相界面分布。
[0069] 在第(3)步中，网格间的流体输运由网格边界和网格中的界面线的几何关系求出。对于满网格来说，可以将F＝1代入式(3)直接得出通过三条边的流体输运量；对于自由面上的网格ABC来说，每条边和网格内界面线的位置关系都可以归纳为图8中所示的六种情形，其中，U是沿当前边外法向的速度分量(此处仅考虑正的分量，如果速度分量为负则直接跳过，因为在考虑其相邻网格上的同一条边时，此处速度分量必为正值，届时再予以计算)，我们在距该边U·dt的距离上做一条辅助线，则位于该辅助线与该网格边之间的流体体积即为dt时间内通过该边的流体输运，在网格ABC中剪掉这部分体积，并在其相邻网格中加上这部分体积，就完成了一条边的计算，对ΔABC的三条边均执行该步骤，继而对流场中所有网格进行遍历，最终得到每个网格中下一时刻的流体体积比函数分布。至此，一个时间步长内的流体输运计算完成，得到新的界面分布。
[0070] 上述PLIC-VOF方法的实施流程如图9所示。
[0071] 使用经典的运动界面数值模拟实验算例对本发明算法进行验证，如图10所示，边长为1的正方形计算区域内分布着按正方形同心环状分布的两种流体，初始t＝0时刻用实线标示，流体在均匀平直速度场的作用下运动到t＝T位置处(如图中虚线所示)，然后设置反向速度场，计算时间2T后正方形同心环将回到初始位置。计算所用的初始网格单元数为7600，图11比较了固定网格和自适应网格上使用PLIC界面重构方法得到的2T时刻的计算结果。从结果直观上可以观察到，使用固定网格时内环的正方形经过2T的输运已经变形成为圆形，产生了很大的失真；而采用自适应网格算法用PLIC方法捕获的相界面优于固定网格的结果，经过2T后仍可保持为正方形。
[0072] 为了进一步定量比较，分析了2T时刻计算结果相对于初始时刻计算结果的形变误差，其定义如下：
[0073]
[0074] 其中 是结束时刻网格单元i中的F值， 是初始时刻网格单元i中的F值，Si是网格单元i的面积。对自适应网格加密算法的结果与未做自适应的固定网格结果进行了对比，结果如表1所示。可以看出，本文算法在没有大量增加网格单元数量的情况下减小了相对形变误差，验证了本发明算法的有效性。
[0075] 表1计算结果对比
[0076]
