辛建建 石伏龙 金秋

(武汉理工大学交通学院,武汉 430063)

一种径向基函数虚拟网格法数值模拟复杂边界流动

辛建建†石伏龙 金秋

(武汉理工大学交通学院,武汉 430063)

(2016年8月7日收到;2016年11月25日收到修改稿)

提出了一种径向基函数虚拟网格法浸入边界法以模拟复杂或多体浸入边界黏性绕流问题.在该方法中,采用有限差分法离散固定笛卡尔交错网格上的不可压缩Navier-Stokes方程,以分步法结合三阶Runge-Kutta格式进行时间积分,高阶TVDMUSCL(total variation diminishing monotonic upstream-centered scheme for conservation laws)格式离散对流项;一个边界连续的虚拟网格法施加物面边界条件以考虑尖锐边界对流场的影响;引入具多项式基的径向基函数描述和重构任意复杂浸入界面,并识别背景网格属性状态.采用Fortran 90语言开发相应的求解器,模拟了绕圆柱、机翼和交错布置圆柱群的黏性绕流问题,验证了本文方法的正确性、可靠性和对复杂边界流动问题的适用性.

浸入边界法,径向基函数,虚拟网格法,复杂边界

1 引 言

数值模拟绕复杂形状物体的黏性绕流问题一直是计算流体力学(computational fluid dynamics,CFD)领域的研究重点和难点.传统研究方法是应用贴体动网格技术,界面附近网格随物体一起运动以精确地描述运动界面.对于小幅运动或单物体时,该方法可得到较好的结果;对于复杂形状物体,生成单域贴体的动网格非常复杂,当物体发生大幅运动或拓扑结构大变形时,网格再生过程会显著降低计算效率,重构过程中也会引入插值误差.

近年来,针对复杂形状多体动边界的数值模拟问题,浸入边界法[1](immersed boundary method,IBM)受到极大的关注.其基于固定的背景笛卡尔网格,网格生成简单,无需进行网格重构,适于模拟涉及复杂几何形状、多体或大变形的动边界问题.根据流-固界面的表示不同,浸入边界法可分为扩散界面法和尖锐界面法[2].对于扩散界面法,由于引入光滑Delta函数或假设有限厚度界面,导致界面模糊化,计算精度受到一定的限制,难以模拟高雷诺数流动.近年来发展的一类尖锐界面法旨在弥补这个基本缺欠,维持界面尖锐是现代CFD技术的发展趋势,如虚拟网格法[3]、切割网格法[4]、嵌入边界法[5].

本文重点强调虚拟网格法,其最初由Majumdar等[6]提出,基本思想是构建物体内的虚拟网格以施加边界条件,考虑尖锐浸入边界对流场的影响.虚拟网格定义为至少有一个相邻流体网格的固体网格,对于每个虚拟网格,需要一个内插模板以隐式施加物面边界条件.本文采用双线性(三线性对于三维)插值[7]方法重构界面内插格式,该方法实施简单、稳定性好.然而对于高雷诺数流动,线性重构可能导致计算错误,一个更高精度的重构方法是在切向线性和在法向进行二次插值,或者采用基于泰勒展开的高阶插值方法[8].

相比需要引入人为刚度系数的反馈力法[9],虚拟网格法避免了时间步相关的稳定性问题;相比在动量方程中引入力源项的嵌入边界法[10],虚拟网格法维持界面尖锐,适用于高雷诺数流动问题;相比网格分类离散的切割网格法[11],虚拟网格法以统一格式离散背景网格,易扩展到三维和曲线网格.发展至今,虚拟网格法已经得到广泛的应用研究,如绕圆柱和机翼的可压缩流动[12]、绕圆柱和平板的不可压缩流动[13]、热传输[14]、昆虫飞行和游鱼运动[15].对于国内研究方面,胡偶等[16]模拟了绕圆柱和机翼的二维无黏可压缩流动.

本文提出一种简单灵活的径向基函数虚拟网格浸入边界法,应用于模拟任意复杂浸入界面的不可压缩流动问题.针对复杂或多体浸入界面流动的模拟,本方法有两个突出要点：1)采用一个边界连续的虚拟网格法施加物面边界条件以隐式考虑浸入边界对流场的影响,物体作为虚拟存在的边界,在整个计算域内采用统一的离散方法求解Navier-Stokes方程;2)以具多项式基的径向基函数(polynomial and radial basis function,PRBF)[17]描述和重构任意复杂浸入界面,并定位和追踪背景网格相状态.基于Fortran 90语言自主开发的二维浸入边界法求解器,模拟了圆柱绕流、机翼绕流和圆柱群的绕流问题.计算结果与已有的数值和试验数据符合较好,验证了本文方法的精度和可靠性,证明了本文方法模拟复杂或多体界面绕流问题的能力.

2 数学模型和数值方法

2.1 控制方程

非定常黏性不可压缩流体的控制方程如下[18]：

速度向量u=(u,v),u,v分别是直角网格x-,y-方向速度分量;ρ是密度,t是时间,p是压力,µ是流体的动力黏性系数.

2.2 数值离散方法

以上流体控制方程采用有限差分法在交错直角网格上离散,速度定义在网格面心,压力、密度、黏性、距离函数等标量定义在网格体心.在该求解器中,采用分步法解耦速度和压力[19].时间推进格式采用三阶Runge-Kutta(RK3)格式,对流项显式处理,黏性项以Crank-Nicolson格式半隐式处理.对于空间离散格式,对流项采用一种高阶TVD MUSCL(total variation diminishing monotonic upstream-centered scheme for conservation laws)[20]格式离散,黏性项以中心差分格式离散.

图1给出了本文流固耦合数值求解方法的流程图,对于该方法的两个关键数值技术,物面边界条件施加和复杂浸入界面描述分别在下文2.3节和2.4节详细介绍.

2.3 虚拟网格浸入边界法

在虚拟网格法中,在每个时间步,虚拟网格变量值通过浸入边界和临近的流体网格变量值插值计算,本文采用双线性插值方法.物面边界条件通过镜像方法直接满足,以确保数值稳定性.

虚拟网格法的关键一步是描述和重构物面边界,以定位识别背景网格相状态(流体网格、虚拟网格、固体网格),具体过程见2.4节.流体网格为网格中心位于流体区域的网格,固体网格为除去虚拟网格后网格中心位于固体区域内的网格.第二步是插值重构虚拟网格变量值,该过程由如下子步骤组成.

图1 浸入边界法数值算法流程图Fig.1.Flowchart of numerical algorithmfor immersed boundary method.

1)定义镜像点：对于每个虚拟网格,穿过浸入边界反射虚拟网格到流体区域惟一确定一个镜像点.本文镜像点定义方法类似于Pan和Shen[22]的方法,其实施简单,可惟一确保镜像点插值模板.定义镜像点到物体表面的间距为ΔL,ΔL是人为指定距离.对于二维问题,dx,dy分别为该流体网格x,y方向的网格尺寸.该ΔL的选择是确保镜像点在被四个流体网格节点形成的一个矩形区域内,而不被边界表面穿过.镜像点坐标计算为

(xIP,yIP)为镜像点网格坐标,(xGC,yGC)为虚拟网格点坐标,(xBI,yBI)为沿着物面法向从虚拟网格中心延伸与物面的投影点坐标,dGB=为虚拟网格点与投影点之间的距离,相关细节如图2.

2)插值计算镜像点流体变量：在上一步中围绕镜像点的四个流体网格节点即为插值模板点,分别用1,2,3,4表示.镜像点值是1,2,3,4节点值的权重叠加,镜像点变量值φIP以双线性插值方法计算：

[φ]T为1,2,3,4节点变量值;[B]为双线性内插格式的系数矩阵,表达式为

(x1,y1),(x2,y2),(x3,y3),(x4,y4)分别为节点1,2,3,4的直角网格坐标值.系数向量[A]从(8)式计算得到,通过(7)式就可计算出镜像点的变量值.

图2 虚拟网格法的二维插值格式Fig.2.Interpolation stencil of ghost cell method.

图3 邻近的网格节点表示Fig.3.Stencil for adjacent grid nodes.

3)计算虚拟网格值：镜像点值φIP计算后,沿着物面法向采用线性内插方法就可计算出虚拟网格节点值φGC：

对于速度,需满足无滑移边界条件,此时C1=-d1/d2和C2=φBI(d1+d2)/d2,φBI为物面节点速度值,d1为虚拟网格到投影点距离,d2为投影点到镜像点距离.对于压力,需要满足Neumann边界条件,此时C1=1和C2=N·(d1+d2),N为物面法向压力导数.

以虚拟网格法施加边界条件后,在全场(包括物面内)求解离散Navier-Stokes方程即可考虑浸入边界对流场的影响.区别于传统虚拟网格法只在物面外计算流场,本文虚拟网格法在整个计算区域求解.与固体区域重合的流体域称为虚拟流体区域,在该区域固体速度和流体网格速度满足无滑移边界条件.由于固体与流场的相互作用是通过边界条件实现,因此固体内虚拟流体的存在对外部流场没有影响.该方法虽然增加了些许计算量,但极大地简化了数据结构和浸入边界的处理,降低了编程难度.

2.4 PRBF界面描述方法

浸入边界法的网格生成过程得到极大简化,但由于浸入边界本身并未包含在网格中,精确表示在背景网格中的浸入边界极具挑战性.本文采用基于PRBF的隐式等值面方法重构任意界面,该方法以光滑基函数拟合离散的物面控制点构造任意复杂物体的等值曲面,并且根据等值面距离函数正负号识别场点属性.PRBF[23,24]界面描述方法简单、精度高,适于模拟复杂浸入界面.

针对任意浸入界面,基于有限数量的界面位置信息的离散控制点,引入PRBF构造如下场函数：

其中M为样品点总数,(xi,yi)为某点i的坐标值,rij表示两点间距离.φ(‖rij‖)在本文采用MQ核函数表示：

形状参数C为常数,这里取C=0.径向基函数权重系数αj满足正交条件：

多项式PF(xi)为

浸入界面的场函数d(xi)均具有等值距离λ(λ＞0),λ为大于零的任意常数,本文取λ=1.因此,等值面方程可转化为求解关于系数αj的(M+3)阶线性方程组.

MQ函数是条件正定的,增加多项式基后,方程(15)的系数矩阵即变为正定矩阵,确保了方程(15)存在惟一确定的解.另外,由于MQ核函数是在所有物面控制点进行(全域求解),其系数矩阵是满阵的.PRBF隐式界面描述方法以有限的物面控制点即可获得较高的精度,采用传统的高斯消元法即可快速地求解等值面方程(15).

进一步构造如下等值面函数ψ(x),以有效识别网格节点ψ(x)在场中属性：

等值面函数ψ(x)具有如下属性：对于任意坐标节点xP(见图3节点xP与邻近网格的符号表示),若ψ(xP)＞0,表示点(xP)在物体内;ψ(xP)＜0表示点(xP)在物体外.在求出全场的等值函数ψ(x)后,下一步是根据图3中网格节点相对位置的距离函数识别定位虚拟网格、流体网格、固体网格 (如图4).

1)如果ψ(xP)≥ 0和 (ψ(xE)≥ 0和ψ(xW)≥0和ψ(xN)≥0和ψ(xS)≥0),则为固体网格(S);

2)如果ψ(xP)≥ 0和 (ψ(xE)＜0或ψ(xW)＜0或ψ(xN)＜0或ψ(xS)＜0),则为虚拟网格(G);

3)如果ψ(xP)＜0和 (ψ(xE)＜0和ψ(xW)＜0和ψ(xN)＜0和ψ(xS)＜0),则为流体网格(F).

图4 背景网格相状态识别(F,流体网格;G,虚拟网格;S,固体网格)Fig.4.Identifying phase state of background grid(F,fl uid cell;G,ghost cell;S,solid cell).

3 算例与结果分析

浸入边界绕流是经典的流动问题,具有丰富的物理现象,广泛地存在于实际的工程应用中.本文将对雷诺数Re=40,100,200的单圆柱绕流,Re=500的NACA0012机翼绕流和Re=100,200的圆柱群绕流问题进行数值模拟,以验证其数值方法的正确性.

3.1 单圆柱的二维黏性流动

在该圆柱绕流算例中,计算区域设置为[-8D,16D]×[-8D,8D],D为圆柱直径,大的计算域使边界远离圆柱域以避免外部边界尾流影响中心流场,几何模型如图5.在入口处采用Dirichlet速度边界条件,顶部和底部均采用对称边界,出口处采用出流边界条件.对于压力边界条件,在出口指定压力为零,其余边界采用均匀Neumann边界条件.直径D设为1 m,入流速度U∞为1 m/s.Re=40,100,200时对应的运动黏性系数ν分别为0.025,0.01,0.005 m2/s(Re=U∞D/ν). 采用两套网格系统480×360和240×180(分别称为M1,M2)进行网格细化研究,并取固定时间步长Δt=0.001 s.

表1比较了三种雷诺数Re=40,100,200下圆柱绕流算例的时间平均阻力系数CD、平均升力系数CL和回流长度Lc,并与其他浸入边界法结果进行了比较.从中可看出,本文虚拟网格法计算的三个关键参数结果(升阻力系数和回流长度)与其他数值结果符合较好,验证了本文方法的精度.另外,两种网格系统的收敛性研究表明低雷诺数下粗细两种网格系统的数值结果趋近相同,说明了本文方法的可靠性.

图5 圆柱绕流的几何模型Fig.5.Geometry model of flow around cylinder.

图6 不同雷诺数下圆柱绕流的瞬时涡场图 (a)Re=40;(b)Re=100;(c)Re=200Fig.6.Instantaneous vorticity contours of flow around cylinder in three Reynolds number：(a)Re=40;(b)Re=100;(c)Re=200.

图6展示了采用网格系统M1计算在Re=40,100和200时圆柱尾流附近的瞬时涡场.在Re=40时,我们能看到一对稳定的对称涡附着在圆柱后侧(如图6(a)),此时流动为定常层流,升力系数为零.随着雷诺数增加,尾流变得不稳定,涡开始周期性交替地从圆柱后面脱落,形成Karman涡街,而且涡泄现象愈发明显,该现象可从Re=100,Re=200时的瞬时涡场图看出(图6(b)和图6(c)).由于Karman涡街的产生,导致圆柱体表面的升阻力发生周期性的变化.该现象可从图7中Re=200时M1网格系统计算的升阻力系数时间演化过程观察到.

表1 圆柱绕流算例中的升阻力系数比较Table 1.Comparison of lift and drag coefficient in the case of flow around cylinder.

图7 Re=200时升阻力系数的时间历程曲线Fig.7.Evolution history of drag and lift coefficients atRe=200.

3.2 翼型绕流

本算例研究雷诺数Re=500的NACA0012机翼绕流问题,以验证本文求解器模拟复杂浸入界面的能力.计算区域取为[-3b,9b]×[-3b,3b],采用600×400的非均匀网格系统,固定的时间步长取为dt=0.001 s,b为翼弦长度.机翼攻角取为0◦,并取机翼前缘点为坐标原点.入口边界为均匀流,在出口采用对流边界条件,自由流边界采用辐射边界.根据Imamura等[26]和李秋[27]的工作选择流体参数,雷诺数Re=U∞b/ν=500.NACA0012机翼的几何构型复杂,难以用解析函数表示,这里采用多项式径向基函数拟合141个样品点表示机翼构型.

当前方法计算的阻力系数 (CD=0.1653)比Imamura等[26]在最细网格计算的值(CD=0.1725)稍小,误差在接受范围内.表明本文CFD求解器和PRBF方法能较好地模拟复杂浸入边界流动问题.从图8的压力等值线图可看出,压力沿着翼弦方向对称布置,在机翼前缘点可观察到一个压力驻点.图9展示了y-轴方向的速度剖面u,v,当前结果与Imamura等[26]的结果总体符合较好.局部差别可能由于两者边界处理方式不同,本文采用在物面内布置虚拟网格的方式施加边界条件.

图8 机翼的瞬时压力场图Fig.8.Instantaneous pressure contour of airfoil.

3.3 交错布置的13个圆柱绕流

这里模拟了13个规则布置的圆柱黏性绕流问题,以验证当前方法对多体浸入边界的适用性.在计算区域[-10D,20D]×[-10D,10D]上布置均匀网格,在垂向和横向相邻的圆柱间距为D,交错圆柱间距为0.5D,时间步长取为0.001 s.图10和图11分别展示了Re=100和Re=200时圆柱群附近位置的压力和涡等值线流场图.从中可看出,在Re=100时,涡场和压力场比较稳定,但已不再是固定不动的对称涡.随着雷诺数增加,在Re=200时流场变成具有复杂涡泄的非定常流动.其与图12中圆心坐标(0,0)处的圆柱在Re=100,Re=200时阻力系数时间历程曲线反映的现象一致.在图10中,Re=100时的阻力系数趋向于稳定,周期性现象不明显,Re=200时阻力系数在经历较长一段时间不规则变化后,产生周期性的升阻力系数时间变化.而在单圆柱绕流中Re=100时周期性涡泄和阻力系数已经发生,表明圆柱之间的相互作用抑制了涡泄的产生.该结果与Gao和Tseng[28]在18圆柱群桩绕流中得出的结论是一致的.另外,在Re=100和Re=200时的阻力系数值明显小于单圆柱绕流计算的结果(对比表1),反映了多体之间的水动力干涉对圆柱表面压力载荷的影响.

图10 Re=100时的瞬时流场图 (a)压力场图;(b)涡场图Fig.10.Instantaneous fluid field of cylinder group atRe=100：(a)Pressure contours;(b)vorticity contours.

这里对多圆柱群桩绕流的水动力特性和多体间水动力干涉的影响进行了初步分析,对于其中的物理机理不在本文研究范围内.分析结果表明,本文方法能较好地考虑多个浸入物体对流场的影响,可应用于模拟多个物体绕流问题.

图11 Re=200时的瞬时流场图 (a)压力场图;(b)涡场图Fig.11.Instantaneous fluid field of cylinder group atRe=200：(a)Pressure contours;(b)vorticity contours.

图12 圆心坐标(0,0)的圆柱在Re=100,200时的阻力系数比较Fig.12.Comparison of drag coefficients of cylinder with center coordinates(0,0)atRe=100,200.

4 结 论

本文提出一种径向基函数虚拟网格浸入边界法模拟复杂边界流动.一个PRBF隐式界面描述方法以表面样品点描述任意复杂几何形状,高效地识别网格属性状态.在整个区域包括物体内部计算Navier-Stokes方程,以虚拟网格隐式考虑浸入边界对流场的影响.在圆柱绕流算例中验证了本文方法的可靠性和精度,在机翼绕流和圆柱群绕流中验证了本文方法模拟复杂几何形状或多体边界绕流的能力.本文方法能轻易地扩展到复杂动边界流动问题,这也是下一阶段我们的研究重点.

[1]Li Q,Li W M 2016Acta Phys.Sin.65 064601(in Chinese)[李强,李五明 2016物理学报 65 064601]

[2]Sotiropoulos F,Yang X 2014Prog.Aerosp.Sci.65 1

[3]Lee J,You D 2013J.Comput.Phys.233 295

[4]Dechristé D,Mieussens L 2015J.Comput.Phys.314 465

[5]Yang J,Stern F 2012J.Comput.Phys.231 5029

[6]Majumdar S,Iaccarino G,Durbin P 2001Annual Research Briefs30 353

[7]Tseng Y H,Ferziger J H 2003J.Comput.Phys.192 593

[8]Xia J,Luo K,Fan J 2015Int.J.Heat Mass Transfer89 856

[9]Yang Q,Cao S Y,Liu S Y 2014Acta Phys.Sin.63 214702(in Chinese)[杨青,曹曙阳,刘十一2014物理学报63 214702]

[10]Monasse L,Daru V,Mariotti C,Piperno S,Tenaud C 2012J.Comput.Phys.231 2977

[11]Schneiders L,Günther C,Meinke M,Schröder W 2016J.Comput.Phys.311 62

[12]Ghias R,Mittal R,Dong H 2007J.Comput.Phys.225 528

[13]Berthelsen P A,Faltinsen O M 2008J.Comput.Phys.227 4354

[14]Xia J,Luo K,Fan J 2014Int.J.Heat Mass Transfer75 302

[15]Mittal R,Dong H,Bozkurttas M,Najjar F M,Vargas A,von Loebbecke A 2008J.Comput.Phys.227 4825

[16]Hu O,Zhao N,Lin J M,Wang D H 2011Acta Aerodyn.Sin.29 491(in Chinese)[胡偶,赵宁,刘剑明,王东红2011空气动力学学报 29 491]

[17]Zhang X,Liu Y,Ma S 2009Adv.Mech.39 1(in Chinese)[张雄,刘岩,马上2009力学进展 39 1]

[18]Versteeg H K,Malalasekera W 2000An Introduction to Computational Fluiddynamics:the Finite Volume Method(Second edition)(Beijing：World Book Inc)p24

[19]Lee J,Kim J,Choi H,Yang K S 2011J.Comput.Phys.230 2677

[20]ShinBR2000EuropeanCongressonComputational Methods in Applied Sciences and Engineering,Barcelona,Spain,September 2-6,2000 p1

[21]Kershaw D S 1978J.Comput.Phys.26 1

[22]Pan D,Shen T T 2009Int.J.Numer.Method Fluid60 1378

[23]Rendall T C S,Allen C B 2007Int.J.Numer.Method Eng.74 10

[24]Ye J J,Li T Q 2013Proceedings of National Congress on HydrodynamicsZhoushan,China,September 20-25,2013 p343(in Chinese)[叶俊杰,李廷秋2013全国水动力学学术会议舟山,中国,2013年9月20-25日第343页]

[25]Frisani A,Hassan Y A 2015Comput.Fluids121 51

[26]Imamura T,Suzuki K,Nakamura T,Yoshida Masahiro 2005J.Comput.Phys.202 645

[27]Li Q 2010M.S.Thesis(Nanjing：Nanjing University of Aeronautics and Astronautics)(in Chinese)[李秋2010硕士学位论文(南京：南京航空航天大学)]

[28]Gao T,Tseng Y H,Lu X Y 2007Int.J.Numer.Method Fluid55 1189

PACS:47.85.Dh,47.11.Bc,45.50.Jf DOI:10.7498/aps.66.044704

Numerical simulation of complex immersed boundaryflow by a radial basis function ghost cell method

Xin Jian-Jian†Shi Fu-Long Jin Qiu

(School of Transportation,Wuhan University of Technology,Wuhan 430063,China)

7 August 2016;revised manuscript

25 November 2016)

A radial basis function ghost cell immersed boundary method of simulating flows around arbitrary complex or multiple immersed boundaries is proposed in this paper.In this method,incompressible Navier-Stokes equations are discretized on fixed Cartesian staggered gridby the finite difference method.A fractional step method is used for time integration,together with third order Runge-Kutta scheme.A high-order TVD MUSCL(total variation diminishing monotonic upstream-centered scheme for conservation law)scheme is used to discretize convective terms.Two salient features are emphasized in the present study.First,boundary conditions at the immersed interface are enforced by a continuous ghost cell method to consider the influence of immersed boundary on the flow field.The immersed bodies are treated as virtual boundaries immersed in the flow.And Navier-Stokes equations are solved in the entire computation domain,including solid domain.Therefore,programming complexity is greatly reduced and the treatment of immersed boundaries is simplified.Second,a polynomial and radial basis function is introduced to implicitly represent and reconstruct arbitrary complex immersed boundaries.Iso-surface distance functions about interface geometries are fitted with some sampling points of body surfaces.It is flexible and robust.Moreover,the information about interface positions on the background grid can be easily identified by the signed distance functions.Based on our in-house developed immersed boundary method solver,typical test cases are simulated to validate the proposed method.Theflows around a cylinder at Reynolds numbers of 40,100 and 200 are first simulated and a grid resolution study is carried out.Good agreement is achieved by comparing with previous numerical results,which shows that this method is accurate and reliable.In the second case of flow around airfoil,the good agreement with previous study shows that the present method has the ability to simulate complex immersed boundary flow.In the last case of flow around array of thirteen cylinders,the ability of present method for multiple immersed boundaries is well proved.And hydrodynamic interaction among multiple bodies is briefly analyzed.

immersed boundary method,radial basis function,ghost cell method,complex boundary

:47.85.Dh,47.11.Bc,45.50.Jf

10.7498/aps.66.044704

†通信作者.E-mail：xinjinjin1990@sina.com

†Corresponding author.E-mail：xinjinjin1990@sina.com