自适应笛卡尔网格Ghost Cell方法研究

刘剑明,赵 宁,胡 偶,王东红

(1.南京航空航天大学航空宇航学院,江苏南京 210016；2.徐州师范大学数学系,江苏徐州221116)

0 引 言

众所周知,计算流体力学中一个持续的障碍是复杂几何外形的网格生成。现今存在的流场离散方法主要包括非结构网格方法、结构网格方法和笛卡尔网格方法[1]。虽然如今网格生成技术得到了进一步的发展,但网格生成过程仍然是一个费时费事的工作。非结构网格主要在二维流场使用三角网格,三维使用四面体或棱柱,其主要优点是易于生成复杂外形的网格,但非结构网格生成费时,计算时间存储量一般都比结构网格高。结构网格单独使用单一的贴体网格很难处理复杂外形,而且会导致网格的高度倾斜,因而一般不会使用单一的结构网格处理高度复杂的几何外形,通常使用嵌套网格和多块对接网格等,虽然这些方法已经得到了成功的使用,但它们需要交换不同网格间的数据,以及处理其他的一些复杂问题。使用笛卡尔网格方法,物体与背景网格切割,或者使用物体内部的虚拟控制体利用嵌入边界来处理与物体的相交。切割网格需要考虑不同情形的复杂切割单元形状,而且切割单元可能会任意小,从而会产生刚性以及时间步长稳定性问题。使用浸入边界Ghost Cell方法可以克服小网格稳定性限制问题,而且更易于实现。笛卡尔网格相对于结构网格、非结构网格来说,具有更多的优势,其网格生成简单,易于实现自动化、具有更强的自适应能力,更适合于处理复杂几何外形的绕流和由于物体运动或变形等产生的非定常问题,故本文研究笛卡尔网格方法。

最近Dadone和Grossman在他们的一系列文章[2-4]中,系统研究了结构网格中的Ghost Cell方法,并把他们命名为ST(Symmetry Technique)方法以及CCST(Curvature Corrected Symmetry Technique)方法,而且这些有效的边界处理方法也被他们推广用于笛卡尔网格[5-7],并取名为GBCM方法。后来,Dadone等[6]利用远场网格加粗,以及Iblanking方法[8]进行自适应。在Wang的文章[9]中,此类方法也被用于二维非结构网格,得到了很好的结果。此外,Forrer[10,11]等提出了一种有趣的浸入Ghost Cell边界方法求解二维静止与运动物体。本文把ST方法以及GBCM方法用于自适应叉树结构的笛卡尔网格,比较了Forrer的浸入边界方法,Dadone的GBCM方法,以及其他的一些可用的边界处理方法。数值结果显示GBCM方法熵误差最小,但总压误差和其他的边界方法基本差不多,而且ST方法与GBCM方法的阻力系数相对于其他方法都较小。此外,本文利用基于叉树结构的自适应笛卡尔网格Ghost Cell方法,求解一个强激波问题以及翼型,数值结果显示本文所使用的自适应方法具有很好的分辨率,且易于实现。

1 控制方程和数值方法

本文仅研究无粘可压流体,考虑二维欧拉方程

其中计算使用MUSCL方法加MINMOD限制器获得高阶精度,采用HLLC离散数值通量[12]

其中

这里K=L,R,q≡unx+vny中间波速估计以及最大最小波速估计采用如下方法[13]

其中λ1(URoe),λ4(URoe)分别是 Roe矩阵[14]的最小最大特征值。在本文的数值模拟中,利用空间分裂法求解多维问题,使用最优二阶 TVD Runge-Kutta方法[15]进行时间推进。

2 Ghost Cell方法

利用笛卡尔网格,物体与网格相交,如图1所示。本文使用Ghost Cell方法需要给出浸入边界(Immersed Boundary)内部网格点处的物理量,比如点A处的值。最直接最简单的方法就是使用一般的对称反射边界条件(ST)

这里点B是A的镜像点。这样给出的边界条件能够满足无穿透边界条件vwall◦n=0。Forrer等人[11]对于浸入固体内部的Ghost Cell值建议速度仍然使用一般的反射,而压强密度使用如下公式

近来,Dadone与Grossman[5-7]系统的给出了一个新颖的Ghost Cell方法(GBCM)处理笛卡尔网格下的静止物体。在他们的文章中,Ghost Cell物理量在墙附近通过一个法向假想的具有局部对称分布的熵(S)与总焓(H)涡流场来得到。这种流体模型满足如下法向动量方程

这里R是带符号的局部曲率半径,如果曲率中心在物体内部为正,反之为负。是切向速度,满足无穿透边界条件vwall◦n=0,此外沿着物体的表面法向强加反对称的熵与总焓的法向导数。这种熵与总焓分布使得当流动是无旋时,产生零法向导数,甚至在墙壁出现涡时也能满足Crocco关系。比如对于二维问题,如图1所示,采用如下边界条件

图1 笛卡尔网格[5-7]Fig.1 Cartesian grid[5-7]

我们称此方法为ECFGCM(Entropy Correction Forrer's Ghost Cell Method)方法,给出这个边界条件的目的在于研究加了熵修正后是否真的能达到熵误差的改善。

在数值计算中Ghost Cell中心A所对应的法向对称点B变量的值可以通过一般插值得到,比如可以使用双线性插值[5,6],或者利用权系数为距离倒数的线性插值公式,对于三维使用三线性插值[7]等,在我们的代码中始终使用流体内部网格点插值。

3 各种Ghost Cell方法的数值比较

为了能够比较各种不同的边界方法,考虑亚声速无粘流体吹向单位圆柱,来流马赫数取为0.38,初始背景网格取为1/4个圆半径,在圆柱附近作四次局部加密,如图2(a)所示。图2(b)给出了利用GBCM 方法得到的等马赫图,其他边界处理方法的等马赫图基本相似,故省略。表1中给出了各种不同边界方法所产生的熵与总压相对误差,从表中可以发现使用GBCM方法熵相对误差要比其他方法约减少40%-70%,加了熵修正后的ECFGCM方法比FGCM稍好约减少20%。总压误差这几种方法都差不多,ST情形稍微大一些。此外表1中还给出了不同边界条件的阻力系数,从表1可以推知ST与GBCM方法阻力系数最小,分别是 3.01e-3与 4.58e-3。此外ST,FGCM,ECFGCM方法都很容易推广到三维,而GBCM就不那么直接了。从计算复杂度考虑,ST方法计算最简单,GBCM方法最复杂,尤其是在三维时,使用GBCM方法[7]需要计算物面与由流线在物面投影方向与物面法向所组成的平面交线的曲率,此外每一步需要进行坐标变换,Dadone的文章[7]给出了详细的叙述。虽然Dadone的文章中指出GBCM方法具有稍好的结果,但对于一些复杂的几何体,为了减少处理曲率的复杂性,尤其是三维问题,以及一些繁琐的坐标变换,可以直接使用ST方法,一般就能得到可以接受的结果。为了和贴体网格的结果比较,我们直接引用文献[9]中的结果。在文献[9]中,Wang等人在非结构网格下使用Roe格式对同样问题进行了计算,在他们的文章中利用128×32个四边形贴体网格得出ST方法相对熵误差与阻力系数分别为3.6e-3与8.0e-3。从我们的计算可以知道,笛卡尔网格的数目要比贴体网格多,计算时间要长一些。使用笛卡尔网格方法虽然网格不贴体,但得到的结果和传统贴体网格方法的结果基本相当,能达到贴体网格的数值精度。

图2 (a)局部加密圆柱笛卡尔网格,(b)GBCM方法(6)等马赫线Fig.2 (a)Local refined Cartesian gird,(b)Mach contours of GBCM(6)

表1 相对误差及阻力系数Table1 Therelativeerrors and drag coefficients

4 叉树结构的自适应笛卡尔网格

为了更有效地处理一些复杂问题,本文把上述一些边界处理方法用于基于叉树结构的自适应笛卡尔网格,为了自动进行解自适应,使用如下自适应判据[17]

应用如下条件判定网格是否需要自适应:

(1)如果 τci＞σc或者 τdi＞σd,控制体 i被标记需要加细。

(2)如果 τci＜1/10σc且τdi＜1/10σd,控制体 i需要加粗。

为了验证本文自适应加密代码的有效性,考虑来流马赫数为3的超声速圆柱绕流问题,使用GBCM方法。图3(a)显示了几何局部加密三次,解自适应加密三次后的网格图,图中精确捕捉了激波的位置,以及尾部需要加密的区域。图3(b)显示了压强等值线图。从图中可以发现GBCM方法以及叉树自适应方法能够有效结合提高解的分辨率,尤其是激波的分辨率。此外为了显示处理复杂几何的能力,用本文方法计算了NACA0012翼型跨声速绕流问题,这里取来流马赫数 Ma∞=0.8,攻角 α=1.25°,自适应判据仅仅使用了速度散度。图4(a)是解自适应加密三次后的网格图,图形清晰地显示了强激波以及弱激波需要加密的位置,图4(b)给出了等压线图,图4(c)是压力系数图,图形具有很好的激波分辨率。

图3 (a)三次解自适应加细后的网格,(b)压强等值线Fig.3 (a)Adaptive Cartesian grid,(b)Pressure contours

图4 翼型跨声速绕流Fig.4 Transonic flow through airfoil

5 结 论

本文主要研究了Ghost Cell方法,比较了各种不同的Ghost Cell边界条件的数值结果,并将其应用于叉树结构的自适应笛卡尔网格,得到了理想的结果。数值实验显示本文方法是十分有效的,并且这些方法可以直接推广到三维问题。

[1]KOH E PC,TSAI H M.Euler solution using Cartesian grid with a gridless least-squares boundary treatment[J].AIAA Journal,2005,43(2):246-255.

[2]DADONE A,GROSSMAN B.Surface boundary conditions for the numerical solution of the Euler equations[R].AIAA 93-3334.1993.

[3]DADONE A.Symmetry techniquefor the numerical solution of the2d Euler equations at impermeable boundaries[J].Int.J.Numer.Meth.Fluids.,1998,28(7):1093-1108.

[4]DADONE A,GROSSMAN B,Surface boundary conditions for thenumerical solution of the Euler equations in three dimensions[J].Lect.Notes.Phys.,1995,453:188-94.

[5]DADONE A,GROSSMAN B.Ghost-Cell method for inviscid two-dimensional flows on cartesian grids[J].AIAA Journal,2004,42(12):2499-2507.

[6]DADONE A,GROSSMAN B.Ghost-cell method with far field coarsening and mesh adaptation for Cartesian grids[J].Computers&Fluids,2006,35(7):676-687.

[7]DADONE A,GROSSMAN B.Ghost-cell method for inviscid three-dimensional flows on cartesian grids[J].Computers&Fluids,2007,36(10):1513-1528.

[8]DURBIN PA,IACCARINO G.An approach to local refinement of structured grids[J].J.Comput.Phys.,2002,181(2):639-53.

[9]WANG Z J,Yuzhi Sun,Curvature-based wall boundary condition for the Euler conditions on unstructured grids[J].AIAA Journal,2003,41(1):27-33.

[10]FORRER H,JELTSCH R.A higher-order boundary treatment for Cartesian-grid methods[J].J.Comput.Phys.,1998,140(2):259-277.

[11]FORRER H,BERGER M.Flow simulations on Cartesian grids involving complex moving geometries[A].Proc.Int.Conf.Hyp.Problems[C].Zurich,Switz,Feb.1998.

[12]BATTEN P,LESCHZINER M A,GOLDBERG U C.Average-state Jacobians and implicit methods for compressible viscous and turbulent flows[J].J.Comput.Phys.,1997,137(1):38-78.

[13]BATTEN P,CLARKE N,LAMBERT C,CAUSON D M.On the choice of wavespeeds for the HLLC Riemann solver[J].SIAM J.Sci.Comput.,1997,18(6):1553-1570.

[14]ROE P L.Approximate riemann solvers,parameter vectors and difference schemes[J].J.Comput.Phys.,1981,43(2):357-372.

[15]SHU C W.Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[R].NASA ICASE Report 97-65,1997.

[16]LÖHNER R,BAUM J,MESTREAU E,SHAROV D.Adaptiveembedded unstructured grid methods[J].Int.J.Numer.Meth.Engng.,2004,60(3):641-660.

[17]DARREN L.DE ZEEUW.A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations[D].[PhD thesis].The University of Michigan,1993.