[发明专利]一种用于求解VOF对流方程的切割体网格THINC方法

权利要求书

1.一种用于求解VOF对流方程的切割体网格THINC方法，其特征在于，包括以下步骤：

步骤1：获取分布两种互不混溶流体A和B的计算域的CFD网格；计算CFD网格中每个单元必要的拓扑信息和几何参数，包括单元与表面、单元与节点、表面与节点、表面与边、边与节点的关系，以及单元体积、表面面积和法向、单元和表面的质心；

步骤2：读入初始速度场和目标时间；读入计算域的CFD网格中两种互不混溶流体A和B的流体分界面的表达式H，获取计算域的CFD网格中各单元Ω_i的初始流体体积分数φ_i；

单元Ω_i的流体体积分数φ_i表示互不混溶流体A和B的流体体积分布比例；φ_i＝1表示单元Ω_i中仅分布一种流体A；φ_i＝0表示单元Ω_i中仅分布一种流体B；具体方法为：

步骤2.1：对计算域的CFD网格中每个单元Ω_i进行标记；

若单元Ω_i的节点全部位于流体分界面内，则将单元Ω_i标记为满单元；若单元Ω_i的节点全部位于流体分界面外，则将单元Ω_i标记为空单元；其余情况则将单元Ω_i标记为界面单元，获取计算域的CFD网格中界面单元的总数N_jm；

步骤2.2：将界面单元Ω_i分割为子单元Ω_ip；

若计算域的CFD网格为3D网格，则将3D界面单元Ω_i分割为四面体子单元Ω_ip；若3D界面单元Ω_i包含多边形表面，则先将多边形表面分割为若干三角形表面，再将3D界面单元Ω_i根据质心和各三角形表面分割为四面体子单元Ω_ip；

若计算域的CFD网格为2D网格，则根据2D界面单元Ω_i的质心和各边节点将2D界面单元Ω_i分割为三角形子单元Ω_ip；

步骤2.3：对子单元Ω_ip进行标记；

若子单元Ω_ip的节点全部位于流体分界面内，则将子单元Ω_ip标记为满子单元；若子单元Ω_ip的节点全部位于流体分界面外，则将子单元Ω_ip标记为空子单元；其余情况则将子单元Ω_ip标记为界面子单元；

步骤2.4：计算界面子单元Ω_ip的流体体积分数φ_ip；

步骤2.4.1：设定最大“分割-判断”级别R、3D单元的体积等分比或2D单元的面积等分比m；初始化r＝1，子单元Ω_ip为被分割的目标单元；

步骤2.4.2：将每个被分割的目标单元等分为m个子单元，对所有子单元进行标记，获取分割生成的满子单元数量N_ipr，找出所有界面子单元；

步骤2.4.3：若r＜R，则将“分割-判断”级别r中分割生成的所有界面子单元作为下一“分割-判断”级别r+1中被分割的目标单元，令r＝r+1，返回步骤2.4.2；

步骤2.4.4：计算界面子单元Ω_ip的流体体积分数φ_ip；

步骤2.5：计算界面单元Ω_i的流体体积分数φ_i；

其中，满子单元的流体体积分数为φ_ip＝1；所述的空子单元的流体体积分数为φ_ip＝0；N_sub为界面单元Ω_i分割出的子单元Ω_ip的数量；

步骤2.6：输出计算域的CFD网格中各单元Ω_i的初始流体体积分数φ_i；

满单元的流体体积分数为φ_i＝1；空单元的流体体积分数为φ_i＝0；

步骤3：判断计算域的CFD网格中各单元Ω_i的类型；

若计算域的CFD网格为3D网格，则3D单元Ω_i由J_i个表面S_ij和K_i个节点δ_ik组成；若计算域的CFD网格为2D网格，则2D单元Ω_i由J_i条边S_ij和K_i个节点δ_ik组成；j＝1，2，...，J_i；k＝1，2，...，K_i；

若3D单元Ω_i为四面体、六面体、三棱柱、四棱锥中的一种，则判定3D单元Ω_i为常规非结构单元；否则，判定3D单元Ω_i为多面体单元；若2D单元Ω_i为三角形或四边形，则判定2D单元Ω_i为常规非结构单元；否则，判定2D单元Ω_i为多边形单元；

步骤4：对于常规非结构单元Ω_i，直接将常规非结构单元Ω_i的各个节点δ_ik按照非结构THINC方法的惯例进行有序编号，计算常规非结构单元Ω_i中各个节点δ_ik在局部坐标系ξ(ξ，η，ζ)下的坐标ξ(δ_ik)，计算常规非结构单元Ω_i中每个表面或每条边S_ij的高斯积分点局部坐标ξ_ijq，q＝1，2，...，Q_ij；

步骤5：对于多面体单元Ω_i，根据各表面S_ij的法向和相对位置关系，将多面体单元Ω_i共面的表面合并为虚拟表面，得到由M个表面Γ_im和N个节点θ_in组成的虚拟单元且虚拟单元为常规非结构单元；

对于多边形单元Ω_i，将多边形单元Ω_i合并为由M条边Γ_im和N个节点θ_in组成的虚拟单元且虚拟单元为常规非结构单元；

其中，m＝1，2，...，M；n＝1，2，...，N；若虚拟单元为四面体，则M＝4，N＝4；若虚拟单元为六面体，则M＝6，N＝8；若虚拟单元为三棱柱，则M＝5，N＝6；若虚拟单元为四棱锥，则M＝5，N＝5；若虚拟单元为三角形，则M＝3，N＝3；若虚拟单元为四边形，则M＝4，N＝4；

步骤6：将虚拟单元的各个节点θ_in按照非结构THINC方法的惯例进行有序编号，再将虚拟单元投影变换到局部坐标系ξ(ξ，η，ζ)下并计算各个节点在局部坐标系ξ(ξ，η，ζ)下的坐标ξ(θ_in)；

步骤7：将多面体单元或多边形单元Ω_i中不属于虚拟单元的原始节点重新标记为悬挂点h_il，l＝1，2，...，L_i；通过对ξ(θ_in)进行插值来计算悬挂点h_il在局部坐标系下的坐标ξ(h_il)；

步骤8：根据ξ(θ_in)和ξ(h_il)计算多面体单元Ω_i的每个表面S_ij的高斯积分点局部坐标ξ_ijq；若S_ij不是三角形或四边形，则先将S_ij的共线边进行虚拟合并，得到虚拟的三角形或四边形表面再计算的高斯积分点坐标赋给S_ij；

对于多边形单元Ω_i，由于S_ij就是由两个节点组成的边，因此直接计算多边形单元Ω_i的每条边S_ij的高斯积分点局部坐标ξ_q；

步骤9：重复步骤4至步骤8，直到遍历计算域的CFD网格中所有单元；

步骤10：标记计算域的CFD网格中所有界面单元Ω_a和需要更新体积分数的目标单元Ω_b；a＝1，2，...，A；b＝1，2，...，B；所述的界面单元Ω_a满足充要条件：0＜φ_a＜1；

步骤11：在局部坐标系ξ(ξ，η，ζ)下重构每个界面单元Ω_a的双曲正切函数若Ω_a为多面体单元，则基于其虚拟单元对进行重构；

步骤12：计算每个目标单元Ω_b中各表面或边S_bj的流体体积分数φ_bj：

其中，φ_bup为S_bj上游单元内的流体体积分数；|S_bj|为表面S_bj的面积或边S_bj的长度；ξ_bjq为S_bj上第q个高斯积分点的局部坐标，ω_q为对应的权重；为S_bj上游单元内的双曲正切函数；

步骤13：计算每个目标单元Ω_b的VOF对流方程的有限体积半离散方程中的非瞬态项φ_b，完成一个时间步长内的计算；

其中，u_nbj为S_bj的外法向速度；

步骤14：判断是否到达目标时间；若未到达目标时间，则读入更新速度场，返回步骤10进行下一个时间步长的计算；若到达目标时间，则输出当前时刻计算域的CFD网格中各单元Ω_i的流体体积分数φ_i，得到计算域中两种互不混溶流体A和B的分布，完成VOF对流方程的求解。

2.根据权利要求1所述的一种用于求解VOF对流方程的切割体网格THINC方法，其特征在于：所述的步骤2.6中输出计算域的CFD网格中各单元Ω_i的流体体积分数φ_i之前需校验是否满足预设精度，具体方法为：

步骤2.6.1：计算所有界面单元包含的流体3D体积或2D面积V_jm在相邻两个“分割-判断”级别k和k-1之间的相对误差

其中，2≤k≤R；代表第u个界面单元；为界面单元在“分割-判断”级别k中的流体体积分数；为界面单元的子单元Ω_up在“分割-判断”级别k中的流体体积分数；

步骤2.6.2：计算CFD网格中所有单元包含的流体3D体积或2D面积V_all在相邻两个“分割-判断”级别k和k-1之间的相对误差

其中，N_cell为计算域的CFD网格中包含的所有单元的数量；

步骤2.6.3：判断与是否满足预设精度；若或不满足预设精度，则增加最大“分割-判断”级别R的值，返回步骤2.4；若与均满足预设精度，则输出计算域的CFD网格中各单元Ω_i的流体体积分数φ_i，得到计算域中两种互不混溶流体A和B的初始分布。