适用二维切割体网格的THINC/QQ算法扩展研究

李超, 段文洋, 赵彬彬

(哈尔滨工程大学 船舶工程学院,黑龙江 哈尔滨 150001)

当采用计算流体力学(computational fluid dynamics,CFD)方法对多相流的自由界面进行模拟时,通常要对界面附近的网格进行局部加密,如波浪的二维传播、变形和破碎等问题。存在海底突起或斜坡等不规则边界时,采用非结构网格进行空间离散更加灵活。其中,切割体网格[1-2]既能支持悬挂点式的高效率加密方法,又能保证网格的整体正交性,因此具有很大优势。

在自由界面捕捉方面,流体体积法(volume of fluid,VOF)由于守恒性能好、能处理自由面大变形流动,故被广泛应用。而双曲正切函数界面捕捉(tangent of hyperbola for interface capturing,THINC)[3]类型VOF算法,通过引入双曲正切函数进行界面重构,有效地抑制了数值耗散。其中,THINC/QQ[4]作为一种混合几何/代数型VOF算法,在没有进行显式几何界面重构的前提下,实现了与几何型VOF算法相当的精度,具有很好的应用前景。对于二维网格,THINC/QQ算法只能处理常规四边形和三角形单元。然而,二维切割体网格还包含更为复杂的多边形单元,例如在粗密网格过渡处会存在悬挂点式四边形/三角形单元,导致THINC/QQ算法不能直接适用。

Xie等[5]提出将THINC/QQ与FVMS3(finite volume method based on merged stencil with 3rd-order reconstruction)[6]方法进行整合(本文将其记作THINC/QQ-FVMS3),能够在任意二维多边形单元中进行界面重构(也可用于三维问题);但是需要将目标单元分割为较多的子单元才能进行单元高斯积分,且三维的基本目标单元应是线性的[6]。Li等[7]提出了适用切割表面式非结构单元的扩展THINC/QQ算法,即THINC/QQ-SF,可用于常见形式的三维切割体网格,以及复杂三维切割体网格中的主要区域。其核心思想来自文献[8],根据复杂单元的主要特征并借助THINC/QQ进行界面重构,当计算VOF流量时再考虑真实单元的拓扑结构。这样就无需采用FVMS3方法,且单元高斯积分也更简便,尤其适合采用固定网格(或弹簧网格)的多相流数值模拟。

本文进一步将THINC/QQ-SF的原理应用到二维切割体网格中,同时设计了一系列测试网格用于典型二维VOF对流和溃坝两相流算例的验证。 为了更好地对比,本文还采用主流CFD软件常用的HRIC[9]算法同步进行模拟。

1 界面捕捉算法的基本原理

1.1 THINC/QQ算法基本原理

考虑物理上互不混溶的流体1和流体2,则流体1的体积分数φ遵循对流方程:

(1)

式中u是流体速度矢量。对于单元Ωi,φi为:

(2)

式中:|Ωi|表示单元Ωi的面积;x=x(x,y)是全局坐标系下的位置矢量;Hi(x)是流体1的指示函数,若x位于流体1中其值为1,否则为0。

(3)

式中:β是控制界面跳跃厚度的斜率因子(在本文中取为6);Pi(ξ)是一个二次多项式,其系数要根据分界面的法向和曲率进行计算,使得THINC/QQ隐式地具有几何型VOF算法的特点。

(4)

di根据式(4)构造非线性方程求解。对方程(1)应用有限体积离散。在Ωi内可以得到如下形式的半离散控制方程:

(5)

THINC/QQ算法默认对方程(5)采用三阶TVD龙格-库塔(RK3)法进行瞬态求解;在计算di时对常规四边形和三角形单元分别采用9点和6点高斯面积分;在计算φij时,在Eij上采用4点高斯线积分。

1.2 二维THINC/QQ-SF算法基本原理

图1 将单元节点分为悬挂节点和重构节点的流程Fig.1 The procedure of separating the cell nodes into hanging nodes and reconstructing nodes

THINC/QQ-FVMS3在计算di时,需要将目标多边形单元分解为J个三角形子单元,共需要3J个高斯积分点[6]。因此,THINC/QQ-SF在每一时间步内要更为简便;尤其是当网格拓扑关系不发生变化时(如固定网格或弹簧网格),单元虚拟简化步骤只需在网格前处理环节进行一次。对于一般的数据结构,线性四边形与含有一个悬挂点的三角形具有相同的存储特征[7];但THINC/QQ-FVMS3自身无法区分这2种单元,而THINC/QQ-SF则能区别对待。

1.3 HRIC及多相流求解算法

不同于THINC类型算法,HRIC作为代数型VOF算法,是基于归一化变量图法 (normalized variable diagram,NVD)保证对流格式的有界性[9],同时实现自由面的高分辨率求解。本文根据STAR-CCM+的默认参数来设置HRIC的相关公式[2]。并且与THINC/QQ一致,本文中HRIC也采用RK3方法对方程(5)进行瞬态求解。

值得说明的是,HRIC的公式中有低阶对流格式(一阶迎风)的贡献。当流动方向与网格的单元面不正交时,一阶迎风格式会使流场变量发生假扩散[10]。所以,在流体分界面发生大变形时,HRIC可能将界面由“清晰”变为“模糊”。

此外,本文还采用压力隐式算子分裂算法 (pressure implicit with splitting of operators,PISO)[11]求解不可压缩真实流动,具体按照文献[12]的方式将其应用于非结构同位网格两相流。对于动量方程,采用一阶隐式欧拉法进行瞬态求解。

在时间步长控制方面,本文的库朗数Co按照OpenFOAM的方式进行定义,对于Ωi有[7]:

(6)

并且同时采用2个最大库朗数,其中Coφmax适用于主体界面单元(0.01<φ<0.99),Comax适用于其他单元。

2 典型对流算例的模拟与分析

VOF对流算例采用定量形状误差[4]:

(7)

式中:Ncel是网格单元总数;φei代表 VOF值的初始解(或精确解);φni代表数值解。

2.1 开槽圆盘刚体旋转

首先以广泛使用的Zalesak开槽圆盘刚体旋转算例[13]进行测试。该圆盘的圆心位于(0.5,0.75),半径为0.15,开槽区域为:

{(x,y)|(|x-0.5|≤0.025)且

(0.6≤y≤0.85)}

(8)

旋转速度场为:

(u,v)=(ω(y-0.5),ω(0.5-x))

(9)

其中x、y∈[0,1]。在本文,ω=2π,且最大求解时间T=1, 相当于绕中心(0.5,0.5)旋转一周;而库朗数设置为:Coφmax=0.25且Comax=1.0×1010。

本例所用的测试网格包括A和B这2类二维切割体网格,每类网格又分为粗、细2个级别,分别记作A1、A2和B1、B2。每个网格都有2个内部加密区域,分别定义为:

{(x,y)|(|x-0.5|≤0.17)且

(|y-0.25|≤0.17)}

(10)

{(x,y)|(|x-0.5|≤0.17)且

(|y-0.75|≤0.17)}

(11)

图2对网格A1和B1及其初始流体分界面进行了展示。在粗、细网格的每个边上分别均匀分布50和100个单元。

图2 开槽圆盘刚体旋转的两类粗网格及其初始流体分界面Fig.2 Two types of coarse meshes for the solid-body rotation of a slotted disk with the initial interface

表1列出了每个网格的单元总数,以及在t=T时刻的E(L1)数值误差。在当前的计算条件下,HRIC的误差是THINC/QQ-SF的3倍以上。

表1 开槽圆盘刚体旋转算例在t=T时刻的E(L1)数值误差Table 1 The numerical errors of E(L1) for the solid-body rotation of a slotted disk at t=T

在定性方面,图3展示了在网格A1和B1上的代表性VOF等值线(0.05、0.5和0.95)。其中,THINC/QQ-SF算法能够在二维切割体网格上保持流体分界面的紧致性;而HRIC算法的数值耗散性则比较明显,尤其是对于较粗的网格。

注:在本文所有彩图中值为0.05、0.5和0.95的VOF等值线分别用蓝色、绿色和红色线表示图3 开槽圆盘刚体旋转算例在t=T时刻网格A1和B1上0.05、0.5和0.95的VOF等值线Fig.3 The VOF contour lines of 0.05, 0.5 and 0.95 on Mesh A1 and B1 for the solid-body rotation of a slotted disk at t=T

2.2 圆盘剪切变形

在二维刚体旋转算例的基础上,采用经典的single vortex圆盘剪切变形算例[14]做进一步测试。同样地,圆盘的圆心位于(0.5,0.75),半径为0.15;其速度场通过流函数进行定义:

(12)

其中:x、y∈[0,1];库朗数Coφmax为0.3且Comax为1.0×1010;T=8。

本例的网格类型也与开槽圆盘刚体旋转一致,只是2个内部加密区域的定义稍有差异,分别为:

{(x,y)|(|x-0.5|≤0.187 5)且(|y-0.25|≤

0.187 5)}、{(x,y)|(|x-0.5|≤0.187 5)且

(|y-0.75|≤0.187 5)}

(13)

在粗、细网格的每个边上分别均匀分布64和128个单元。表2列出了每个网格的单元总数,以及在t=T时刻的E(L1)数值误差。在本文计算条件下,HRIC的误差同样都大于THINC/QQ-SF的误差,且前者处于后者的1～2倍。

表2 圆盘剪切变形算例在t=T时刻的E(L1)数值误差Table 2 The numerical errors of E(L1) for the single-vortex shearing flow of a disk at t=T

图4定性地展示了在中间时刻t=T/2流体分界面(即VOF的0.5等值线)结果,可以看到,即使在分界面被严重拉伸变形的情况下,THINC/QQ-SF算法依然能在粗密网格之间实现分界面的光顺过渡,而且比HRIC的数值解耗散性更低。

图4 圆盘剪切变形算例在t=T/2时刻网格A2和B2上的流体分界面及其附近单元Fig.4 The interfaces along with the neighboring cells on mesh A2 and B2 for the single-vortex shearing flow of a disk at t=T/2

3 溃坝两相流算例模拟与分析

为了进一步验证二维THINC/QQ-SF算法的性能,采用典型二维溃坝算例[15]进行研究,包括有障碍物和无障碍物2种条件。计算域在水平x方向长0.584 m,在垂直y方向高0.34 m,初始水域为:

{(x,y)|(0≤x≤a)且(0≤y≤2a)},

a=0.146 m

(14)

在物理属性方面,水和空气的密度分别设置为1 000 kg/m3和1.0 kg/m3,而运动粘度分别为1.0×10-6m2/s和1.48×10-5m2/s;重力加速度值g=9.81 m/s2;且不考虑表面张力。库朗数设置为:Coφmax=0.25且Comax=0.5。

3.1 有障碍物溃坝流动

网格的单元总数为9 841,并进行局部加密。在边界条件方面,顶部为压力入口(压力等于0 Pa),其余均为不可滑移壁面。图5为典型时刻的VOF等值线对比结果,其中,THINC/QQ-SF和HRIC都能捕捉到流动的主要特征、所得的自由面在粗密网格之间可以光顺过渡;但THINC/QQ-SF的自由面更紧致,尤其是在界面发生大变形时。

图5 有障碍物溃坝的瞬时0.05、0.5和0.95 VOF等值线Fig.5 The transient VOF contour lines of 0.05, 0.5 and 0.95 for the dambreak flow with an obstacle

3.2 无障碍物溃坝流动

计算域外尺寸和边界条件与有障碍物溃坝一致,网格单元总数为7 424。图6展示了在t=1.0 s的VOF等值线,虽然HRIC同样能捕捉到典型位置的气泡,但是THINC/QQ-SF在控制界面厚度方面更有优势。

图6 无障碍物溃坝在1.0 s的0.05、0.5和0.95VOF等值线Fig.6 The VOF contour lines of 0.05, 0.5 and 0.95 for the dambreak flow without any obstacles at 1.0 s

图7 无障碍物溃坝流的自由面无量纲位移时历曲线Fig.7 Time histories of the non-dimensional interface displacements for the dambreak flow without any obstacles

在垂向位移方面,本文所用的2种VOF算法数值解与文献[12]的CICSAM数值解(采用8 400个网格单元)几乎重合,并都与文献[16]的实验值非常一致。

在水平位移方面,不同VOF算法的数值解与实验值[16-17]相比都存在一定的偏差,但整体趋势比较吻合。这与水密闸门的运动、实验测量方法等都有一定关系,有系统误差的影响。

整体来讲,在界面形变相对较小时,HRIC的界面质量与THINC/QQ-SF无显著差别。当界面出现较大翻卷(尤其是破碎)时,THINC/QQ-SF更有助于抑制VOF的数值耗散,并能很好地适用于二维切割体网格。

4 结论

1)本文将THINC/QQ-SF算法的原理应用到二维切割体网格中,使得原始THIINC/QQ算法能够用于悬挂点式四边形/三角形单元的界面重构。而且,与THINC/QQ-FVMS3相比,THINC/QQ-SF需要的面高斯积分点更少、更简便。

2)针对典型VOF对流算例(开槽圆盘刚体旋转和圆盘剪切变形)以及溃坝两相流算例,本文设计了一系列二维切割体网格,验证了所用扩展方法的可行性。同时,与代数型的HRIC算法相比,THIINC/QQ-SF的精度更高、得到的自由面更紧致,能够很好地继承原始THIINC/QQ的优点。

3)对于需要局部加密或计算域不规则的二维数值模拟,如波浪传播、经过海底凸起后的波浪变形、在斜坡上的波浪破碎等自由表面问题,THIINC/QQ-SF算法通过与二维切割体网格进行配合,具有很好的应用前景。