CN113178011B - 一种用于求解vof对流方程的切割体网格thinc方法 - Google Patents

一种用于求解vof对流方程的切割体网格thinc方法 Download PDF

Info

Publication number
CN113178011B
CN113178011B CN202110467955.8A CN202110467955A CN113178011B CN 113178011 B CN113178011 B CN 113178011B CN 202110467955 A CN202110467955 A CN 202110467955A CN 113178011 B CN113178011 B CN 113178011B
Authority
CN
China
Prior art keywords
unit
omega
interface
grid
phi
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110467955.8A
Other languages
English (en)
Other versions
CN113178011A (zh
Inventor
李超
段文洋
赵彬彬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Engineering University
Original Assignee
Harbin Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN202110467955.8A priority Critical patent/CN113178011B/zh
Publication of CN113178011A publication Critical patent/CN113178011A/zh
Application granted granted Critical
Publication of CN113178011B publication Critical patent/CN113178011B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • Computer Graphics (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明属于CFD技术领域，具体涉及一种用于求解VOF对流方程的切割体网格THINC方法。本发明针对切割体网格单元的特点，对于其中包含多边形表面的多面体单元，在保留单元实际拓扑的条件下，通过合并单元的共面表面形成属于常规非结构类型的虚拟单元，可以直接应用现有THINC方法进行双曲正切函数的重构；在基于重构的双曲正切函数计算原始实际单元各表面的VOF值时，对于多边形表面，先将其共线边进行合并后形成虚拟三角形或四边形表面，再进行表面高斯数值积分，因此更简单高效。本发明能够在3D切割体网格上基于非结构THINC方法的基本原理实现VOF对流方程的高精度求解。

Description

一种用于求解VOF对流方程的切割体网格THINC方法
技术领域
本发明属于CFD技术领域，具体涉及一种用于求解VOF对流方程的切割体网格THINC方法。
背景技术
在CFD(Computational Fluid Dynamics)计算中，VOF(Volume of Fluid)方法被广泛应用于互不混溶多相流体的分界面模拟，其原理是通过求解流体体积分数VOF的对流方程来实现对流体分界面的动态捕捉；对于一个给定的物理计算域，CFD技术需要先将其划分为网格，而网格由一系列特定类型的单元组成；所谓流体体积分数VOF就是某相流体在网格单元中所占的体积比(2D时为面积比)；若将流体A和B所在的区域分别记作ΩA和ΩB，并定义一个流体指示函数H：
Figure BDA0003044044900000011
其中x＝(x,y,z)代表位置矢量，则在网格单元Ωi内，流体A的体积分数VOFi(简记为φi)定义为：
Figure BDA0003044044900000012
其中，|Ωi|代表一维单元的长度或二维单元的面积或三维单元的体积，图2展示了流体A和B的真实分布情况和流体A在网格单元中的VOF值；在物理上，φ遵循下述对流方程：
Figure BDA0003044044900000013
其中，t表示时间，u为速度矢量；因此，借助VOF方法捕捉流体运动分界面的关键就是如何高精度地求解上述VOF(φ)对流方程，即如何根据当前时刻tn的流体分布φn和速度场vn来准确得到下一时刻tn+1的流体分布φn+1。
求解VOF对流方程的方法主要分为几何重构类和代数类两种；通常，几何重构类方法的精度高，但算法非常复杂，如PLIC(Piecewise Linear Interface Calculation)方法；代数类VOF方法一般相对简单、易于拓展到非结构网格中，但精度有限，如HRIC(High-Resolution Interface Capturing)方法；而THINC(Tangent of Hyperbola InterfaceCapturing)方法，通过引入双曲正切函数和充分考虑分界面的几何信息，在不进行显式几何重构的条件下，实现了与几何重构类方法相当的精度，并保持了代数类方法的简洁性，因此具有很好的应用前景；THINC方法的基本原理首先以一维重构进行说明，对于如图3所示的一维网格，真实流体指示函数H(x)为阶梯型，而THINC方法通过构造一维双曲正切函数
Figure BDA0003044044900000014
来对H(x)进行近似：
Figure BDA0003044044900000021
其中，β为斜率因子，用于控制
Figure BDA0003044044900000022
曲线的陡峭程度；γ控制着
Figure BDA0003044044900000023
跳跃的方向；参数di代表跳跃中心的位置，通过下述积分等式求得：
Figure BDA0003044044900000024
而THINC方法在多维重构时，采用的多维双曲正切函数形式如下：
Figure BDA0003044044900000025
其中，Pi(x)为一个多项式，其系数依赖于单元Ωi及其周边单元的φ值分布，Pi(x)+di＝0代表单元Ωi内流体分界面的表面方程，类似地，参数di由下述积分等式求解：
Figure BDA0003044044900000026
在所有界面单元(即满足0＜φ＜1的单元)的
Figure BDA0003044044900000027
重构完成后，可以用
Figure BDA0003044044900000028
求解VOF对流方程(8)：
Figure BDA0003044044900000029
首先，对公式(3)采用有限体积离散，则在单元Ωi内有如下形式的半离散控制方程：
Figure BDA00030440449000000210
其中，Sij为单元Ωi的第j个表面(或2D单元时的边)，|Sij|为Sij的面积(或2D单元时的长度)，J为总表面数，unij为Sij的外法向速度，φij为Sij上的体积分数值，通过下述公式进行计算：
Figure BDA00030440449000000211
其中，φiup为Sij上游单元内的流体体积分数，xq＝(xq,yq,zq)为Sij上第q个高斯积分点的坐标，Q为Sij上高斯积分点的总数，ωq为对应的权重，
Figure BDA00030440449000000212
为Sij上游单元内的双曲正切函数，当unij＞0时，iup＝i；通常，对方程(8)采用三阶TVD Runge-Kutta法进行时间积分，即在每一时间步长内对计算域内的φij和φj重复更新三次。
由于实际的物理计算域形状一般比较复杂，故通常采用非结构网格进行求解；非结网格THINC方法都采用多维双曲正切函数重构，但在重构
Figure BDA0003044044900000031
的方式上有所不同，因此导致不同THINC方法能适用的网格单元类型也存在差异；其中，UMTHINC(UnstructuredMTHINC)和THINC/QQ(THINC method with Quadratic surface representation andgaussian Quadrature)方法，需要将目标单元从全局坐标系x(x,y,z)投影变换到局部坐标系ξ(ξ,η,ζ)下，通常ξ,η,ζ的取值范围为[0,1]或[-1,1]，图4展示了三角形单元的投影变换示意图，此时
Figure BDA0003044044900000032
和Pi(x)应分别记作
Figure BDA0003044044900000033
和Pi(ξ)，公式(9)中的高斯积分点xq也应记作ξq＝(ξq,ηq,ζq)；由于受到投影变换公式的限制，UMTHINC和THINC/QQ只能处理如图5所示的常规2D三角形、四边形单元和常规3D四面体、六面体、三棱柱、四棱锥单元类型，而且在计算多项式Pi(ξ)的系数和ξq的坐标时需要将单元的无序节点θ按图5所示进行有序排列，详见文献“Bin Xie,Feng Xiao,Toward efficient and accurate interface capturing onarbitrary hybrid unstructured grids:the THINC method with quadratic surfacerepresentation and gaussian quadrature,2017”。
文献“Bin Xie,Peng Jin,Yanping Du,ShiJun Liao,A consistent andbalanced-force model for incompressible multiphase flows on polyhedralunstructured grids,2020”中提出的THINC方法，可以直接在全局坐标系x(x,y,z)下对目标单元的双曲正切函数
Figure BDA0003044044900000034
进行重构，不仅能处理图5所示的单元类型，还能够处理如图6所示的2D任意多边形单元和图7所示的表面为三角形或四边形的3D任意多面体单元类型；详见该方法所在文献的参考文献“Bin Xie,Xi Deng,Shijun Liao,High-fidelity solveron polyhedral unstructured grids for low-Mach number compressible viscousflow,2019”中的“The linear elements are used where the boundary surfaceΓijcomprises two vertices in 2D and three or four vertices in 3D”相关内容。
在不同类型的非结构网格中，切割体网格(Trimmed Mesh)由于其主体采用正交直角切割网格，靠近物体边界时可以采用贴体网格，既能够保证求解质量，又便于控制网格单元数量，因此被广泛应用；在切割体网格中，3D单元可以是表面为任意多边形的任意多面体，2D单元可以是任意多边形，图8(a)和图8(b)分别对典型的3D和2D切割体网格进行了展示，图9为3D切割体网格中包含的典型单元类型示意图；因此，现有的非结构网格THINC方法已经能够适用于2D切割体网格，但不能用于3D切割体网格中。
发明内容
本发明的目的在于克服现有非结构THINC方法不能在3D切割体网格上求解VOF对流方程的不足，提供一种用于求解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；
步骤3：判断计算域的CFD网格中各单元Ωi的类型；
若计算域的CFD网格为3D网格，则3D单元Ωi由Ji个表面Sij和Ki个节点δik组成；若计算域的CFD网格为2D网格，则2D单元Ωi由Ji条边Sij和Ki个节点δik组成；j＝1,2,…,Ji；
k＝1,2,…,Ki；
若3D单元Ωi为四面体、六面体、三棱柱、四棱锥中的一种，则判定3D单元Ωi为常规非结构单元；否则，判定3D单元Ωi为多面体单元；若2D单元Ωi为三角形或四边形，则判定2D单元Ωi为常规非结构单元；否则，判定2D单元Ωi为多边形单元；
步骤4：对于常规非结构单元Ωi，直接将常规非结构单元Ωi的各个节点δik按照非结构THINC方法的惯例进行有序编号，计算常规非结构单元Ωi中各个节点δik在局部坐标系ξ(ξ,η,ξ)下的坐标ξ(δik)，计算常规非结构单元Ωi中每个表面或每条边Sij的高斯积分点局部坐标ξijq，q＝1,2,…,Qij；
步骤5：对于多面体单元Ωi，根据各表面Sij的法向和相对位置关系，将多面体单元Ωi共面的表面合并为虚拟表面，得到由M个表面Γim和N个节点θin组成的虚拟单元
Figure BDA0003044044900000041
且虚拟单元
Figure BDA0003044044900000042
为常规非结构单元；
对于多边形单元Ωi，将多边形单元Ωi合并为由M条边Γim和N个节点θin组成的虚拟单元
Figure BDA0003044044900000043
且虚拟单元
Figure BDA0003044044900000044
为常规非结构单元；
其中，m＝1,2,…,M；n＝1,2,…,N；若虚拟单元
Figure BDA0003044044900000045
为四面体，则M＝4,N＝4；若虚拟单元
Figure BDA0003044044900000046
为六面体，则M＝6,N＝8；若虚拟单元
Figure BDA0003044044900000047
为三棱柱，则M＝5,N＝6；若虚拟单元
Figure BDA0003044044900000048
为四棱锥，则M＝5,N＝5；若虚拟单元
Figure BDA0003044044900000049
为三角形，则M＝3,N＝3；若虚拟单元
Figure BDA00030440449000000410
为四边形，则M＝4,N＝4；
步骤6：将虚拟单元
Figure BDA0003044044900000051
的各个节点θin按照非结构THINC方法的惯例进行有序编号，再将虚拟单元
Figure BDA0003044044900000052
投影变换到局部坐标系ξ(ξ,η,ζ)下并计算各个节点在局部坐标系ξ(ξ,η,ζ)下的坐标ξ(θin)；
步骤7：将多面体单元或多边形单元Ωi中不属于虚拟单元
Figure BDA0003044044900000053
的原始节点重新标记为悬挂点hil，l＝1,2,…,Li；通过对ξ(θin)进行插值来计算悬挂点hil在局部坐标系下的坐标ξ(hil)；
步骤8：根据ξ(θin)和ξ(hil)计算多面体单元Ωi的每个表面Sij的高斯积分点局部坐标ξijq；若Sij不是三角形或四边形，则先将Sij的共线边进行虚拟合并，得到虚拟的三角形或四边形表面
Figure BDA0003044044900000054
再计算
Figure BDA0003044044900000055
的高斯积分点坐标赋给Sij；
对于多边形单元Ωi，由于Sij就是由两个节点组成的边，因此直接计算多边形单元Ωi的每条边Sij的高斯积分点局部坐标ξq；
步骤9：重复步骤4至步骤8，直到遍历计算域的CFD网格中所有单元；
步骤10：标记计算域的CFD网格中所有界面单元Ωa和需要更新体积分数的目标单元Ωb；
a＝1,2,…,A；b＝1,2,…,B；所述的界面单元Ωa满足充要条件：0＜φa＜1；
步骤11：在局部坐标系ξ(ξ,η,ζ)下重构每个界面单元Ωa的双曲正切函数
Figure BDA0003044044900000056
若Ωa为多面体单元，则基于其虚拟单元
Figure BDA0003044044900000057
对
Figure BDA0003044044900000058
进行重构；
步骤12：计算每个目标单元Ωb中各表面或边Sbj的流体体积分数φbj：
Figure BDA0003044044900000059
其中，φbup为Sbj上游单元内的流体体积分数；|Sbj|为表面Sbj的面积或边Sbj的长度；ξbjq为Sbj上第q个高斯积分点的局部坐标，ωq为对应的权重；
Figure BDA00030440449000000510
为Sbj上游单元内的双曲正切函数；
步骤13：计算每个目标单元Ωb的VOF对流方程的有限体积半离散方程中的非瞬态项φb，完成一个时间步长内的计算；
Figure BDA00030440449000000511
其中，unbj为Sbj的外法向速度；
步骤14：判断是否到达目标时间；若未到达目标时间，则读入更新速度场，返回步骤10进行下一个时间步长的计算；若到达目标时间，则输出当前时刻计算域的CFD网格中各单元Ωi的流体体积分数φi，得到计算域中两种互不混溶流体A和B的分布，完成VOF对流方程的求解。
本发明还可以包括：
所述的步骤2中读入计算域的CFD网格中两种互不混溶流体A和B的流体分界面的表达式H，获取计算域的CFD网格中各单元Ωi的初始流体体积分数φi的方法具体为：
步骤2.1：对计算域的CFD网格中每个单元Ωi进行标记；
若单元Ωi的节点全部位于流体分界面内，则将单元Ωi标记为满单元；若单元Ωi的节点全部位于流体分界面外，则将单元Ωi标记为空单元；其余情况则将单元Ωi标记为界面单元，获取计算域的CFD网格中界面单元的总数Njm；
步骤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个子单元，对所有子单元进行标记，获取分割生成的满子单元数量Nipr，找出所有界面子单元；
步骤2.4.3：若r＜R，则将“分割-判断”级别r中分割生成的所有界面子单元作为下一“分割-判断”级别r+1中被分割的目标单元，令r＝r+1，返回步骤2.4.2；
步骤2.4.4：计算界面子单元Ωip的流体体积分数φip；
Figure BDA0003044044900000071
步骤2.5：计算界面单元Ωi的流体体积分数φi；
Figure BDA0003044044900000072
其中，满子单元的流体体积分数为φip＝1；所述的空子单元的流体体积分数为φip＝0；Nsub为界面单元Ωi分割出的子单元Ωip的数量；
步骤2.6：输出计算域的CFD网格中各单元Ωi的初始流体体积分数φi；
满单元的流体体积分数为φi＝1；空单元的流体体积分数为φi＝0。
所述的步骤2.6中输出计算域的CFD网格中各单元Ωi的流体体积分数φi之前需校验是否满足预设精度，具体方法为：
步骤2.6.1：计算所有界面单元包含的流体3D体积或2D面积Vjm在相邻两个“分割-判断”级别k和k-1之间的相对误差
Figure BDA0003044044900000073
Figure BDA0003044044900000074
Figure BDA0003044044900000075
Figure BDA0003044044900000076
其中，2≤k≤R；
Figure BDA0003044044900000077
代表第u个界面单元；
Figure BDA0003044044900000078
为界面单元
Figure BDA0003044044900000079
在“分割-判断”级别k中的流体体积分数；
Figure BDA00030440449000000710
为界面单元
Figure BDA00030440449000000711
的子单元Ωup在“分割-判断”级别k中的流体体积分数；
Figure BDA00030440449000000712
步骤2.6.2：计算CFD网格中所有单元包含的流体3D体积或2D面积Vall在相邻两个“分割-判断”级别k和k-1之间的相对误差
Figure BDA00030440449000000713
Figure BDA0003044044900000081
Figure BDA0003044044900000082
其中，Ncell为计算域的CFD网格中包含的所有单元的数量；
步骤2.6.3：判断
Figure BDA0003044044900000083
与
Figure BDA0003044044900000084
是否满足预设精度；若
Figure BDA0003044044900000085
或
Figure BDA0003044044900000086
不满足预设精度，则增加最大“分割-判断”级别R的值，返回步骤2.4；若
Figure BDA0003044044900000087
与
Figure BDA0003044044900000088
均满足预设精度，则输出计算域的CFD网格中各单元Ωi的流体体积分数φi，得到计算域中两种互不混溶流体A和B的初始分布。
本发明的有益效果在于：
本发明提供了一种用于求解VOF对流方程的切割体网格THINC方法，能够在3D切割体网格上基于非结构THINC方法的基本原理实现VOF对流方程的高精度求解。本发明能够针对切割体网格单元的特点，对于其中包含多边形表面的多面体单元，在保留单元实际拓扑的条件下，通过合并单元的共面表面形成属于常规非结构类型的虚拟单元，可以直接应用现有THINC方法进行双曲正切函数的重构。在基于重构的双曲正切函数计算原始实际单元各表面的VOF值时，对于多边形表面，先将其共线边进行合并后形成虚拟三角形或四边形表面，再进行表面高斯数值积分，因此更简单高效。本发明克服了现有非结构THINC方法不能在3D切割体网格上求解VOF对流方程的不足，且本发明的双曲正切函数重构方式同样也适用于二维切割体网格中的多边形单元。
附图说明
图1是本发明的总体流程图。
图2是两种流体A和B的真实分布情况和流体A在网格单元中的VOF值示意图。
图3是THINC方法的构造一维双曲正切函数的示意图。
图4是现有非结构THINC方法中将2D三角形单元投影变换为局部坐标系下的标准单元示意图。
图5是常规2D三角形、四边形单元和常规3D四面体、六面体、三棱柱、四棱锥单元示意图。
图6是2D任意多边形单元示意图。
图7是表面为三角形或四边形的3D任意多面体单元示意图。
图8(a)是典型3D切割体网格示意图。
图8(b)是典型2D切割体网格示意图。
图9是3D切割体网格中包含的典型单元类型示意图。
图10是本发明将3D多面体单元的共面表面虚拟合并的示意图。
图11是本发明将3D单元的多边形表面虚拟转化为三角形或四边形表面的示意图。
图12是本发明将2D多边形单元的共线边虚拟合并的示意图。
图13(a)是本发明测试算例的物理计算域、边界条件和流体分类。
图13(b)是本发明测试算例的CFD切割体网格和初始VOF云图(即初始流体分布)。
图13(c)是初始VOF值的0.05、0.5、0.95等值线(分别用蓝、绿、红线表示)。
图13(d)是本发明THINC方法计算得到的t＝0.5s时刻VOF上述等值线。
图13(e)是本发明THINC方法计算得到的t＝1.0s时刻VOF上述等值线。
图13(f)是本发明THINC方法计算得到的t＝1.5s时刻VOF上述等值线。
图13(g)是本发明THINC方法计算得到的t＝2.0s时刻VOF上述等值线。
图13(h)是本发明THINC方法计算得到的t＝4.0s时刻VOF上述等值线。
图13(i)是MULES方法计算得到的最终时刻VOF上述等值线。
图13(j)是本发明THINC方法计算得到的最终时刻VOF＝0.5的等值线与精确解的对比图。
图13(k)是MULES方法计算得到的最终时刻VOF＝0.5的等值线与精确解的对比图。
具体实施方式
下面结合附图对本发明做进一步描述。
本发明提供了一种用于求解VOF对流方程的切割体网格THINC方法，能够在3D切割体网格上基于非结构THINC方法的基本原理实现VOF对流方程的高精度求解。本发明能够针对切割体网格单元的特点，对于其中包含多边形表面的多面体单元，在保留单元实际拓扑的条件下，通过合并单元的共面表面形成属于常规非结构类型的虚拟单元，可以直接应用现有THINC方法进行双曲正切函数的重构。在基于重构的双曲正切函数计算原始实际单元各表面的VOF值时，对于多边形表面，先将其共线边进行合并后形成虚拟三角形或四边形表面，再进行表面高斯数值积分，因此更简单高效。本发明克服了现有非结构THINC方法不能在3D切割体网格上求解VOF对流方程的不足，且本发明的双曲正切函数重构方式同样也适用于二维切割体网格中的多边形单元。
一种用于求解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；
步骤3：判断计算域的CFD网格中各单元Ωi的类型；
若计算域的CFD网格为3D网格，则3D单元Ωi由Ji个表面Sij和Ki个节点δik组成；若计算域的CFD网格为2D网格，则2D单元Ωi由Ji条边Sij和Ki个节点δik组成；j＝1,2,…,Ji；k＝1,2,…,Ki；
若3D单元Ωi为四面体、六面体、三棱柱、四棱锥中的一种，则判定3D单元Ωi为常规非结构单元；否则，判定3D单元Ωi为多面体单元；若2D单元Ωi为三角形或四边形，则判定2D单元Ωi为常规非结构单元；否则，判定2D单元Ωi为多边形单元；
步骤4：对于常规非结构单元Ωi，直接将常规非结构单元Ωi的各个节点δik按照非结构THINC方法的惯例进行有序编号，计算常规非结构单元Ωi中各个节点δik在局部坐标系ξ(ξ,η,ζ)下的坐标ξ(δik)，计算常规非结构单元Ωi中每个表面或每条边Sij的高斯积分点局部坐标ξijq，q＝1,2,…,Qij；
步骤5：对于多面体单元Ωi，根据各表面Sij的法向和相对位置关系，将多面体单元Ωi共面的表面合并为虚拟表面，得到由M个表面Γim和N个节点θin组成的虚拟单元
Figure BDA0003044044900000101
且虚拟单元
Figure BDA0003044044900000102
为常规非结构单元；
对于多边形单元Ωi，将多边形单元Ωi合并为由M条边Γim和N个节点θin组成的虚拟单元
Figure BDA0003044044900000103
且虚拟单元
Figure BDA0003044044900000104
为常规非结构单元；
其中，m＝1,2,…,M；n＝1,2,…,N；若虚拟单元
Figure BDA0003044044900000105
为四面体，则M＝4,N＝4；若虚拟单元
Figure BDA0003044044900000106
为六面体，则M＝6,N＝8；若虚拟单元
Figure BDA0003044044900000107
为三棱柱，则M＝5,N＝6；若虚拟单元
Figure BDA0003044044900000108
为四棱锥，则M＝5,N＝5；若虚拟单元
Figure BDA0003044044900000109
为三角形，则M＝3,N＝3；若虚拟单元
Figure BDA00030440449000001010
为四边形，则M＝4,N＝4；
步骤6：将虚拟单元
Figure BDA00030440449000001011
的各个节点θin按照非结构THINC方法的惯例进行有序编号，再将虚拟单元
Figure BDA00030440449000001012
投影变换到局部坐标系ξ(ξ,η,ζ)下并计算各个节点在局部坐标系ξ(ξ,η,ζ)下的坐标ξ(θin)；
步骤7：将多面体单元或多边形单元Ωi中不属于虚拟单元
Figure BDA00030440449000001013
的原始节点重新标记为悬挂点hil，l＝1,2,…,Li；通过对ξ(θin)进行插值来计算悬挂点hil在局部坐标系下的坐标ξ(hil)；
步骤8：根据ξ(θin)和ξ(hil)计算多面体单元Ωi的每个表面Sij的高斯积分点局部坐标ξijq；若Sij不是三角形或四边形，则先将Sij的共线边进行虚拟合并，得到虚拟的三角形或四边形表面
Figure BDA0003044044900000111
再计算
Figure BDA0003044044900000112
的高斯积分点坐标赋给Sij；
对于多边形单元Ωi，由于Sij就是由两个节点组成的边，因此直接计算多边形单元Ωi的每条边Sij的高斯积分点局部坐标ξq；
步骤9：重复步骤4至步骤8，直到遍历计算域的CFD网格中所有单元；
步骤10：标记计算域的CFD网格中所有界面单元Ωa和需要更新体积分数的目标单元Ωb；
a＝1,2,…,A；b＝1,2,…,B；所述的界面单元Ωa满足充要条件：0＜φa＜1；
步骤11：在局部坐标系ξ(ξ,η,ζ)下重构每个界面单元Ωa的双曲正切函数
Figure BDA0003044044900000113
若Ωa为多面体单元，则基于其虚拟单元
Figure BDA0003044044900000114
对
Figure BDA0003044044900000115
进行重构；
步骤12：计算每个目标单元Ωb中各表面或边Sbj的流体体积分数φbj：
Figure BDA0003044044900000116
其中，φbup为Sbj上游单元内的流体体积分数；|Sbj|为表面Sbj的面积或边Sbj的长度；ξbjq为Sbj上第q个高斯积分点的局部坐标，ωq为对应的权重；
Figure BDA0003044044900000117
为Sbj上游单元内的双曲正切函数；
步骤13：计算每个目标单元Ωb的VOF对流方程的有限体积半离散方程中的非瞬态项φb，完成一个时间步长内的计算；
Figure BDA0003044044900000118
其中，unbj为Sbj的外法向速度；
步骤14：判断是否到达目标时间；若未到达目标时间，则读入更新速度场，返回步骤10进行下一个时间步长的计算；若到达目标时间，则输出当前时刻计算域的CFD网格中各单元Ωi的流体体积分数φi，得到计算域中两种互不混溶流体A和B的分布，完成VOF对流方程的求解。
所述的步骤2中读入计算域的CFD网格中两种互不混溶流体A和B的流体分界面的表达式H，获取计算域的CFD网格中各单元Ωi的初始流体体积分数φi的方法具体为：
步骤2.1：对计算域的CFD网格中每个单元Ωi进行标记并计算其流体体积分数；
若单元Ωi的节点全部位于流体分界面内，则将单元Ωi标记为满单元；若单元Ωi的节点全部位于流体分界面外，则将单元Ωi标记为空单元；其余情况则将单元Ωi标记为界面单元，获取计算域的CFD网格中界面单元的总数Njm；
所述的满单元的流体体积分数为φi＝1；所述的空单元的流体体积分数为φi＝0；对于界面单元Ωi，其流体体积分数φi的计算方法为：
步骤2.1.1：将界面单元Ωi分割为子单元Ωip；
若计算域的CFD网格为3D网格，则将3D界面单元Ωi分割为四面体子单元Ωip；若3D界面单元Ωi包含多边形表面，则先将多边形表面分割为若干三角形表面，再将3D界面单元Ωi根据质心和各三角形表面分割为四面体子单元Ωip；
若计算域的CFD网格为2D网格，则根据2D界面单元Ωi的质心和各边节点将2D界面单元Ωi分割为三角形子单元Ωip；
步骤2.1.2：对子单元Ωip进行标记并计算其流体体积分数；
若子单元Ωip的节点全部位于流体分界面内，则将子单元Ωip标记为满子单元；若子单元Ωip的节点全部位于流体分界面外，则将子单元Ωip标记为空子单元；其余情况则将子单元Ωip标记为界面子单元；
所述的满子单元的流体体积分数为φip＝1；所述的空子单元的流体体积分数为φip＝0；对于界面子单元Ωip，其流体体积分数φip的计算方法为：
步骤2.1.2.1：设定最大“分割-判断”级别R、3D单元的体积等分比或2D单元的面积等分比m；初始化r＝1，子单元Ωip为被分割的目标单元；
步骤2.1.2.2：将每个被分割的目标单元等分为m个子单元，对所有子单元进行标记，获取分割生成的满子单元数量Nipr，找出所有界面子单元；
步骤2.1.2.3：若r＜R，则将“分割-判断”级别r中分割生成的所有界面子单元作为下一“分割-判断”级别r+1中被分割的目标单元，令r＝r+1，返回步骤2.1.2.2；
步骤2.1.2.4：计算界面子单元Ωip的流体体积分数φip；
Figure BDA0003044044900000131
步骤2.1.3：计算界面单元Ωi的流体体积分数φi；
Figure BDA0003044044900000132
其中，Nsub为界面单元Ωi分割出的子单元Ωip的数量；
步骤2.2：计算所有界面单元包含的流体3D体积或2D面积Vjm在相邻两个“分割-判断”级别k和k-1之间的相对误差
Figure BDA0003044044900000133
Figure BDA0003044044900000134
Figure BDA0003044044900000135
Figure BDA0003044044900000136
其中，2≤k≤R；
Figure BDA0003044044900000137
代表第u个界面单元；
Figure BDA0003044044900000138
为界面单元
Figure BDA0003044044900000139
在“分割-判断”级别k中的流体体积分数；
Figure BDA00030440449000001310
为界面单元
Figure BDA00030440449000001311
的子单元Ωup在“分割-判断”级别k中的流体体积分数；
Figure BDA00030440449000001312
步骤2.3：计算CFD网格中所有单元包含的流体3D体积或2D面积Vall在相邻两个“分割-判断”级别k和k-1之间的相对误差
Figure BDA00030440449000001313
Figure BDA00030440449000001314
Figure BDA00030440449000001315
其中，Ncell为计算域的CFD网格中包含的所有单元的数量；
步骤2.4：判断
Figure BDA00030440449000001316
与
Figure BDA00030440449000001317
是否满足预设精度；若
Figure BDA00030440449000001318
或
Figure BDA00030440449000001319
不满足预设精度，则增加最大“分割-判断”级别R的值，返回步骤2.1；若
Figure BDA00030440449000001320
与
Figure BDA00030440449000001321
均满足预设精度，则输出计算域的CFD网格中各单元Ωi的流体体积分数φi，得到计算域中两种互不混溶流体A和B的初始分布。
本发明不限制于求解两种互不混溶流体A和B在计算域的分布，也可用于求解计算域中存在多种互不混溶时，各流体在计算域的分布。只需先取一种流体，将该种流体视为流体A，其余所有类型的流体均视为流体B，先利用本发明求取流体A在计算域中分布，然后再从B中取一种流体视为A，再次计算流体A在计算域中分布。重复上述步骤，直至得到全部流体在计算域中的分布。
与现有技术相比，本发明的有益效果是：
(1)本发明能够在3D切割体网格上基于非结构THINC方法的基本原理实现VOF对流方程的高精度求解；
(2)本发明能够针对切割体网格单元的特点，在保留单元实际拓扑的条件下，通过合并单元的共面表面形成属于常规非结构类型的虚拟单元，可以直接应用现有THINC方法进行双曲正切函数的重构，在基于重构的双曲正切函数计算原始实际单元各表面的VOF值时，对于多边形表面，先将其共线边进行合并后形成虚拟三角形或四边形表面，再进行表面高斯数值积分，因此更简单高效。
(3)本发明的双曲正切函数重构方式同样也适用于二维切割体网格中的多边形单元。
实施例1：
本发明涉及一种求解VOF对流方程的方法，尤其涉及一种在切割体网格上求解VOF对流方程的THINC方法。
步骤1，读入计算域CFD网格，如Ansys Fluent的.msh格式的网格文件，并计算每个单元必要的拓扑信息和几何参数，包括单元与表面、单元与节点、表面与节点、表面与边、边与节点的关系，以及单元体积、表面面积和法向、单元和表面的质心；
步骤2，读入初始的流体分布(即VOF值)和速度场；
步骤3，判断网格中各单元的类型，并按常规非结构单元(四面体、六面体、三棱柱、四棱锥)和多面体单元进行分类；
步骤4，对于多面体单元Ωi，其由表面Sij(j＝1,2,…,J)组成，其节点为δik(k＝1,2,…,K)，根据其各表面Sij的法向和相对位置关系，将Ωi共面的表面合并为虚拟表面，得到新的虚拟单元
Figure BDA0003044044900000141
该虚拟单元由表面Γim(m＝1,2,…,M)和节点θun(n＝1,2,…,N)组成，且
Figure BDA0003044044900000142
属于常规非结构四面体(M＝4,N＝4)、六面体(M＝6,N＝8)、三棱柱(M＝5,N＝6)、四棱锥(M＝5,N＝5)单元中的一种，图10中步骤①和步骤②展示了
Figure BDA0003044044900000143
为六面体的情形，即在数据结构上，先将待合并表面的公共边去掉，再把共线边合并；
步骤5，将虚拟单元
Figure BDA0003044044900000144
的各个节点θin，按照非结构THINC方法的惯例进行有序编号，如图10中步骤③所示，再将
Figure BDA0003044044900000145
投影变换到局部坐标系ξ(ξ,η,ζ)下并计算各个节点在局部坐标系下的坐标ξ(θin)，而每个ξ(θin)的具体值可以参考THINC/QQ或UMTHINC相关文献很方便地得到；
步骤6，将Ωi中不属于
Figure BDA0003044044900000151
的原始节点重新标记为悬挂点hil(l＝1,2,…,L)，如图10中步骤④所示，并通过对ξ(θin)进行插值来计算悬挂点hil在局部坐标系下的坐标ξ(hil)，具体采用线性插值的方法，例如对于图10中的hi1，有：
ξ(hi1)＝λξ(θi3)+(1-λ)ξ(θi4)
其中，插值因子λ根据全局坐标系下各点的实际距离计算如下：
Figure BDA0003044044900000152
同理可以求得其他位于
Figure BDA0003044044900000153
各边上的悬挂点坐标，而对于处在
Figure BDA0003044044900000154
的表面上的悬挂点，以图10中的hi2节点为例，可以借助其他辅助点来计算，例如：
ξ(hi2)＝αξ(hi1)+(1-α)ξ(hi3)
类似地，插值因子α为：
Figure BDA0003044044900000155
也可以构造新辅助点来计算；
步骤7，根据ξ(θin)和ξ(hil)计算Ωi的每个表面Sij的高斯积分点局部坐标ξq(q＝1,2,…,Q)，若Sij不是三角形或四边形，则先将Sij的共线边进行虚拟合并，得到虚拟的三角形或四边形表面
Figure BDA0003044044900000156
如图11所示，再计算
Figure BDA0003044044900000157
的高斯积分点坐标赋给Sij，而Q和ξq的具体值及其对应的权重ωq可以参考现有THINC/QQ或UMTHINC或其他数学方面的文献很方便的得到。
步骤8，重复步骤3至步骤7，直到遍历所有单元，其中，对于常规非结构单元，直接进行单元节点的有序编号、节点局部坐标和每个单元表面的高斯积分点局部坐标的计算；
步骤9，标记所有界面单元Ωa(a＝1,2,…,A)和需要更新体积分数的目标单元Ωb(b＝1,2,…,B)，其中A为界面单元的总数，B为目标单元的总数，某一单元成为界面单元的充要条件为0＜φ界面单元＜1，而成为目标单元的充分条件包括以下三条：
(1)若0＜φc＜1，则Ωc为目标单元
(2)若φc＝1且Ωc至少存在一个相邻单元Ωcn满足0＜φcn＜1，则Ωc为目标单元
(3)若Ωd至少存在一个相邻单元Ωc满足条件(1)或条件(2)，则Ωd为目标单元；
步骤10，重构每个界面单元Ωa的双曲正切函数
Figure BDA0003044044900000158
若Ωa为多面体单元，则基于其虚拟单元
Figure BDA0003044044900000159
对
Figure BDA00030440449000001510
进行重构，默认在局部坐标系ξ(ξ,η,ξ)下进行重构，得到双曲正切函数
Figure BDA00030440449000001511
如下所示：
Figure BDA0003044044900000161
也可以在全局坐标系x(x,y,z)下进行重构双曲正切函数，即：
Figure BDA0003044044900000162
其中，Pa(ξ)或Pa(x)以及da的计算方法和公式可以参考现有非结构THINC算法的相关文献而直接得到；
步骤11，计算每个目标单元Ωb的所有表面Sbj(j＝1,2,…,Jb)的流体体积分数φbj：
Figure BDA0003044044900000163
其中，Jb为Ωb的总表面数，φbup为Sbj上游单元内的流体体积分数，|Sbj|为Sbj的面积，ξr为Sbj上第r个高斯积分点的坐标，若Sbj不是三角形或四边形则ξr为其虚拟表面
Figure BDA0003044044900000164
的高斯积分点坐标，R为Sbj高斯积分点的总数，ωr为对应的权重，
Figure BDA0003044044900000165
为Sbj上游单元内的双曲正切函数；
步骤12，计算每个目标单元Ωb的VOF对流方程
Figure BDA0003044044900000166
的有限体积半离散方程
Figure BDA0003044044900000167
中的非瞬态项，其中，|Ωb|为Ωb的体积，unbj为Sbj的外法向速度；
步骤13，重复步骤9至步骤12，直到完成一个时间步内的时间积分；按照非结构THINC算法的惯例，在每一时间步内采用三阶TVD龙格库塔法进行时间积分；
步骤14，更新速度场u和时间步长，并判断是否输出当前时刻的VOF值，若是则输出；可以根据CFL条件数调整时间步长，也可以采用固定时间步长；计算结果可以按照等时间步长输出，或按照等时间间隔输出，结果文件可以采用Tecplot格式进行输出；
步骤15，重复步骤9至步骤14，直到完成VOF对流方程的求解；
本发明还包括这样一些技术特征：
所述步骤10，也可以直接在全局坐标系x(x,y,z)下进行重构，得到双曲正切函数
Figure BDA0003044044900000168
当采用全局坐标系重构时，只需将ξ(θin)、ξ(hil)、ξq和ξr分别替换为相同节点在全局坐标系下的坐标值x(θin)、x(hil)、xq和xr；
本发明的双曲正切函数重构方式同样也适用于二维切割体网格中的多边形单元，即：在步骤3中，将网格中各单元并按常规非结构单元(三角形、四边形)和多边形单元进行分类，在步骤4中，Sij和Γim分别代表多边形单元Ωi和虚拟单元
Figure BDA0003044044900000171
的边，且
Figure BDA0003044044900000172
属于三角形(M＝3,N＝3)或四边形(M＝4,N＝4)，由于Sij就是由两个节点组成的边，因此在步骤7中可以直接对Sij计算其高斯积分点坐标而不必将其转化为
Figure BDA0003044044900000173
其他方面与三维网格的原理相同，图12展示了
Figure BDA0003044044900000174
为三角形的情形；
为了进一步展示本发明的技术效果，考虑如图13(a)所示的算例进行测试，该算例所用CFD数值水槽长5m、高1m，槽内充满两种互不混溶流体A和B，其中流体A所在的区域形状为圆形，以水槽左下角点为原点，则圆心坐标为(0.5,0.5),半径为0.25m，入口和内部流体的速度为：
u＝(ux,uy)＝(1,0)
定常流动，且只有流体B不断从入口边界流入水槽中，流动总时间为4.0s；图13(b)中展示了该算例所用的切割体网格以及流体A的初始VOF云图(即初始分布)，其中网格为三维单层网格(宽0.1m)、z方向边界设置为对称条件；通过求解VOF对流方程来对水槽内流体A和B的运动界面进行捕捉；
对于本发明的切割体网格THINC方法，在步骤9中采用THINC/QQ方法的二次表面重构(且在目标单元的局部坐标系下进行)，并设置CFL条件数为0.25；为了与本发明进行对比，在主流开源CFD软件OpenFOAM的环境中采用广泛使用的MULES算法进行相同条件的计算；由于z向速度为0，所以为了方便观察，只对xy平面内的VOF等值线进行展示；
通常在VOF方法中，以VOF＝0.5的等值线(或等值面)来代表真实的流体分界面，同时也关注0.05和0.95等值线的分布情况；图13(c)是初始VOF值的0.05、0.5、0.95等值线(分别用蓝、绿、红线表示)，图13(d)至图13(h)是本发明计算得到的不同时刻VOF上述等值线，图13(i)是MULES方法计算得到的最终时刻VOF上述等值线，图13(j)和图13(k)分别是本发明和MULES方法计算得到的最终时刻VOF＝0.5的等值线与精确解的对比图；可以明显看到，本发明的切割体网格THINC方法比MULES方法精度更高、对流体运动分界面的捕捉效果更好。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

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网格中界面单元的总数Njm；
步骤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个子单元，对所有子单元进行标记，获取分割生成的满子单元数量Nipr，找出所有界面子单元；
步骤2.4.3：若r＜R，则将“分割-判断”级别r中分割生成的所有界面子单元作为下一“分割-判断”级别r+1中被分割的目标单元，令r＝r+1，返回步骤2.4.2；
步骤2.4.4：计算界面子单元Ωip的流体体积分数φip；
Figure FDA0003626288090000021
步骤2.5：计算界面单元Ωi的流体体积分数φi；
Figure FDA0003626288090000022
其中，满子单元的流体体积分数为φip＝1；所述的空子单元的流体体积分数为φip＝0；Nsub为界面单元Ωi分割出的子单元Ωip的数量；
步骤2.6：输出计算域的CFD网格中各单元Ωi的初始流体体积分数φi；
满单元的流体体积分数为φi＝1；空单元的流体体积分数为φi＝0；
步骤3：判断计算域的CFD网格中各单元Ωi的类型；
若计算域的CFD网格为3D网格，则3D单元Ωi由Ji个表面Sij和Ki个节点δik组成；若计算域的CFD网格为2D网格，则2D单元Ωi由Ji条边Sij和Ki个节点δik组成；j＝1，2，...，Ji；k＝1，2，...，Ki；
若3D单元Ωi为四面体、六面体、三棱柱、四棱锥中的一种，则判定3D单元Ωi为常规非结构单元；否则，判定3D单元Ωi为多面体单元；若2D单元Ωi为三角形或四边形，则判定2D单元Ωi为常规非结构单元；否则，判定2D单元Ωi为多边形单元；
步骤4：对于常规非结构单元Ωi，直接将常规非结构单元Ωi的各个节点δik按照非结构THINC方法的惯例进行有序编号，计算常规非结构单元Ωi中各个节点δik在局部坐标系ξ(ξ，η，ζ)下的坐标ξ(δik)，计算常规非结构单元Ωi中每个表面或每条边Sij的高斯积分点局部坐标ξijq，q＝1，2，...，Qij；
步骤5：对于多面体单元Ωi，根据各表面Sij的法向和相对位置关系，将多面体单元Ωi共面的表面合并为虚拟表面，得到由M个表面Γim和N个节点θin组成的虚拟单元
Figure FDA0003626288090000023
且虚拟单元
Figure FDA0003626288090000024
为常规非结构单元；
对于多边形单元Ωi，将多边形单元Ωi合并为由M条边Γim和N个节点θin组成的虚拟单元
Figure FDA0003626288090000025
且虚拟单元
Figure FDA0003626288090000026
为常规非结构单元；
其中，m＝1，2，...，M；n＝1，2，...，N；若虚拟单元
Figure FDA0003626288090000027
为四面体，则M＝4，N＝4；若虚拟单元
Figure FDA0003626288090000031
为六面体，则M＝6，N＝8；若虚拟单元
Figure FDA0003626288090000032
为三棱柱，则M＝5，N＝6；若虚拟单元
Figure FDA0003626288090000033
为四棱锥，则M＝5，N＝5；若虚拟单元
Figure FDA0003626288090000034
为三角形，则M＝3，N＝3；若虚拟单元
Figure FDA0003626288090000035
为四边形，则M＝4，N＝4；
步骤6：将虚拟单元
Figure FDA0003626288090000036
的各个节点θin按照非结构THINC方法的惯例进行有序编号，再将虚拟单元
Figure FDA0003626288090000037
投影变换到局部坐标系ξ(ξ，η，ζ)下并计算各个节点在局部坐标系ξ(ξ，η，ζ)下的坐标ξ(θin)；
步骤7：将多面体单元或多边形单元Ωi中不属于虚拟单元
Figure FDA00036262880900000315
的原始节点重新标记为悬挂点hil，l＝1，2，...，Li；通过对ξ(θin)进行插值来计算悬挂点hil在局部坐标系下的坐标ξ(hil)；
步骤8：根据ξ(θin)和ξ(hil)计算多面体单元Ωi的每个表面Sij的高斯积分点局部坐标ξijq；若Sij不是三角形或四边形，则先将Sij的共线边进行虚拟合并，得到虚拟的三角形或四边形表面
Figure FDA0003626288090000038
再计算
Figure FDA0003626288090000039
的高斯积分点坐标赋给Sij；
对于多边形单元Ωi，由于Sij就是由两个节点组成的边，因此直接计算多边形单元Ωi的每条边Sij的高斯积分点局部坐标ξq；
步骤9：重复步骤4至步骤8，直到遍历计算域的CFD网格中所有单元；
步骤10：标记计算域的CFD网格中所有界面单元Ωa和需要更新体积分数的目标单元Ωb；a＝1，2，...，A；b＝1，2，...，B；所述的界面单元Ωa满足充要条件：0＜φa＜1；
步骤11：在局部坐标系ξ(ξ，η，ζ)下重构每个界面单元Ωa的双曲正切函数
Figure FDA00036262880900000310
若Ωa为多面体单元，则基于其虚拟单元
Figure FDA00036262880900000311
对
Figure FDA00036262880900000312
进行重构；
步骤12：计算每个目标单元Ωb中各表面或边Sbj的流体体积分数φbj：
Figure FDA00036262880900000313
其中，φbup为Sbj上游单元内的流体体积分数；|Sbj|为表面Sbj的面积或边Sbj的长度；ξbjq为Sbj上第q个高斯积分点的局部坐标，ωq为对应的权重；
Figure FDA00036262880900000314
为Sbj上游单元内的双曲正切函数；
步骤13：计算每个目标单元Ωb的VOF对流方程的有限体积半离散方程中的非瞬态项φb，完成一个时间步长内的计算；
Figure FDA0003626288090000041
其中，unbj为Sbj的外法向速度；
步骤14：判断是否到达目标时间；若未到达目标时间，则读入更新速度场，返回步骤10进行下一个时间步长的计算；若到达目标时间，则输出当前时刻计算域的CFD网格中各单元Ωi的流体体积分数φi，得到计算域中两种互不混溶流体A和B的分布，完成VOF对流方程的求解。
2.根据权利要求1所述的一种用于求解VOF对流方程的切割体网格THINC方法，其特征在于：所述的步骤2.6中输出计算域的CFD网格中各单元Ωi的流体体积分数φi之前需校验是否满足预设精度，具体方法为：
步骤2.6.1：计算所有界面单元包含的流体3D体积或2D面积Vjm在相邻两个“分割-判断”级别k和k-1之间的相对误差
Figure FDA0003626288090000042
Figure FDA0003626288090000043
Figure FDA0003626288090000044
Figure FDA0003626288090000045
其中，2≤k≤R；
Figure FDA0003626288090000046
代表第u个界面单元；
Figure FDA0003626288090000047
为界面单元
Figure FDA0003626288090000048
在“分割-判断”级别k中的流体体积分数；
Figure FDA0003626288090000049
为界面单元
Figure FDA00036262880900000410
的子单元Ωup在“分割-判断”级别k中的流体体积分数；
Figure FDA00036262880900000411
步骤2.6.2：计算CFD网格中所有单元包含的流体3D体积或2D面积Vall在相邻两个“分割-判断”级别k和k-1之间的相对误差
Figure FDA00036262880900000412
Figure FDA00036262880900000413
Figure FDA00036262880900000414
其中，Ncell为计算域的CFD网格中包含的所有单元的数量；
步骤2.6.3：判断
Figure FDA0003626288090000051
与
Figure FDA0003626288090000052
是否满足预设精度；若
Figure FDA0003626288090000053
或
Figure FDA0003626288090000054
不满足预设精度，则增加最大“分割-判断”级别R的值，返回步骤2.4；若
Figure FDA0003626288090000055
与
Figure FDA0003626288090000056
均满足预设精度，则输出计算域的CFD网格中各单元Ωi的流体体积分数φi，得到计算域中两种互不混溶流体A和B的初始分布。
CN202110467955.8A 2021-04-28 2021-04-28 一种用于求解vof对流方程的切割体网格thinc方法 Active CN113178011B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110467955.8A CN113178011B (zh) 2021-04-28 2021-04-28 一种用于求解vof对流方程的切割体网格thinc方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110467955.8A CN113178011B (zh) 2021-04-28 2021-04-28 一种用于求解vof对流方程的切割体网格thinc方法

Publications (2)

Publication Number Publication Date
CN113178011A CN113178011A (zh) 2021-07-27
CN113178011B true CN113178011B (zh) 2022-08-02

Family

ID=76926893

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110467955.8A Active CN113178011B (zh) 2021-04-28 2021-04-28 一种用于求解vof对流方程的切割体网格thinc方法

Country Status (1)

Country Link
CN (1) CN113178011B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115344902B (zh) * 2022-10-18 2023-03-10 中国空气动力研究与发展中心计算空气动力研究所 自由面重构方法、装置、设备及介质
CN116522828B (zh) * 2023-07-04 2023-10-20 中国空气动力研究与发展中心计算空气动力研究所 非结构线性三棱柱网格单元重构方法、系统、设备及介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102129517A (zh) * 2011-03-10 2011-07-20 西安交通大学 一种高精度的两相流体界面捕获方法
CN109726465A (zh) * 2018-12-26 2019-05-07 电子科技大学 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2317348B1 (en) * 2009-10-30 2014-05-21 Services Pétroliers Schlumberger Method for building a depositional space corresponding to a geological domain
CN106445882A (zh) * 2016-07-12 2017-02-22 南京航空航天大学 一种由vof函数快速构造符号距离函数的改进clsvof方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102129517A (zh) * 2011-03-10 2011-07-20 西安交通大学 一种高精度的两相流体界面捕获方法
CN109726465A (zh) * 2018-12-26 2019-05-07 电子科技大学 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法

Also Published As

Publication number Publication date
CN113178011A (zh) 2021-07-27

Similar Documents

Publication Publication Date Title
CN113178011B (zh) 一种用于求解vof对流方程的切割体网格thinc方法
CN109377561B (zh) 一种基于共形几何的数模表面网格生成方法
RU2519331C2 (ru) Способ вычисления физического значения, способ численного анализа, программа вычисления физического значения, программа численного анализа, устройство вычисления физического значения и устройство численного анализа
CN102306396B (zh) 一种三维实体模型表面有限元网格自动生成方法
CN111489447B (zh) 一种适用于格子Boltzmann方法的直角网格自适应建模方法
Avila et al. A parallel CFD model for wind farms
CN108447120A (zh) 一种构建三维属性体模型的方法及装置
Lahnert et al. Towards lattice-Boltzmann on dynamically adaptive grids–minimally-invasive grid exchange in ESPResSo
CN115758938A (zh) 面向粘性边界流场数值模拟的附面层网格生成方法
Jude et al. An Octree-based, Cartesian CFD Solver for Helios on CPU and GPU Architectures.
Reiter et al. Mesh generation for thin layered domains and its application to parallel multigrid simulation of groundwater flow
Gerritsma et al. The geometric basis of numerical methods
CN113868931B (zh) 复合材料有限元建模方法、系统及存储介质
Fries Higher-order accurate integration for cut elements with Chen-Babuška nodes
CN111047687B (zh) 一种基于三维t样条的异质材料实体建模方法
Loch The level set method for capturing interfaces with applications in two-phase flow problems
Koren Euler flow solutions for a transonic windtunnel section
CN113177373B (zh) 一种基于vof原理的流体分布计算方法
Scholz et al. High-order quadrature on planar domains based on transport theorems for implicitly defined moving curves
Dreyer The local discontinuous galerkin method for the advection-diffusion equation on adaptive meshes
Chang et al. Automatic multigrid generation for an unstructured parallel overset-grid solver
CN115374567B (zh) 叶盘轮毂加工路径生成方法、装置、介质和电子设备
WO2021024295A1 (ja) データ構造、数値計算方法、および数値計算プログラム
Ferrer et al. A robust arbitrarily high order transport method of the characteristic type for unstructured tetrahedral grids
US11282274B1 (en) Systems and methods for constructing conformal connections between meshes

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant