一种基于求解椭圆型方程的结构动网格生成方法

卢凤翎， 陈小前， 禹彩辉， 苗萌

1.国防科技大学 航天科学与工程学院， 长沙 410073 2.空间物理重点实验室， 北京 100076

一种基于求解椭圆型方程的结构动网格生成方法

卢凤翎1,2,*， 陈小前1， 禹彩辉2， 苗萌2

1.国防科技大学 航天科学与工程学院， 长沙 410073 2.空间物理重点实验室， 北京 100076

基于椭圆型网格生成法，实现了一种简单高效的贴体结构动网格生成方法，可用于具有移动边界问题的非定常流动数值模拟。该方法提出，在网格变形过程中，Poisson方程需要的控制网格间距和正交性的源项可以通过提取已知的静态网格源项直接得到，并在整个动网格生成过程中保持不变。因此，在椭圆型网格生成中需要通过外迭代确定源项的过程可以得到省略，而且该方法不需要人工指定参数。这使得方法具有高效和易于嵌入到已有程序中的特点。数值模拟结果证明，采用这种方法获得的网格能够较好地保持静态网格原有的正交性和光滑性，在相同迭代步数约束下，网格求解效率低于传统弹簧模拟法,但鲁棒性优于弹簧模拟法。

动网格生成； 椭圆型网格生成； 计算流体力学； 数值网格生成； 网格变形方法

工程应用中经常碰到具有移动边界的非定常流体动力学问题，例如火箭级间分离、飞机机弹分离、飞机俯仰运动和涡轮机内流动等。要处理这些问题，一般需要涉及到动网格生成技术。也就是说，网格需要跟随物体运动或产生变形以适应物体的运动。

对某些问题，例如进行飞机俯仰运动的非定常模拟，可以考虑采用网格不变形但随物体同时移动的所谓刚性网格方法去解决；但在处理具有相对运动部件的问题时，例如多段翼中的襟翼运动模拟以及火箭级间分离等问题的模拟，就不能采用刚性网格方法，而需要考虑其他的方法，例如重叠网格 (Overlapping Grid) ，也称为Chimera 网格或者嵌套 (Overset) 网格的方法[1-4]。这种方法允许网格间重叠、嵌套。当物体运动时，贴体的部件网格随之运动，数值计算在各个分块子网格上分别进行，并在重叠或嵌套的区域间通过插值来实现流场信息的传递。重叠网格虽然降低了网格生成难度，但是挖洞、找点使得网格构造变得极其耗时，这种方法的计算效率成为一个非常突出的问题。

此外,在求解上述问题时,还可以考虑采用动网格技术,即在每一时间步对网格进行部分或整体重新生成。网格重新生成时，如果网格规模和拓扑结构发生了变化，则称之为网格重构法；如果各个时间步上的网格与原始网格在规模和拓扑结构上均保持一致，则称之为网格变形法。按此分类，网格变形法只是将物面的刚性运动或者局部变形传递到空间网格，对空间网格的坐标进行适当更新，而不改变网格的规模和拓扑结构，效率更高，也更容易与流场求解相结合，因而更受欢迎，也是动网格技术研究的重要内容。本文以结构网格变形技术为研究对象，但只涉及物面作刚性运动情况下的网格变形方法。

实际上，有多种方法可以实施网格变形[5]，例如超限插值(Trans-Finite Interpolation，TFI)法、弹簧模拟法、弹性体法、Delaunay图映射法和径向基函数(Radial Basis Functions, RBF)法等。超限插值[6-7]法属于代数插值类方法，一般只用于结构网格，该方法计算量小，可以实现相对复杂的网格变形，但难以保证变形后的网格质量，特别是物面网格的正交性。而其他几种方法既可用于结构网格，也可用于非结构网格。其中，弹簧模拟法最初是由Batina[8]在1990年提出的，他把这种方法应用于求解一个受到强迫振动的翼型。在其原始形式中，单元的每一边被假设为一个具有刚度的弹簧，并且刚度与其边长成反比。物体运动引起的网格变形反映在网格节点的位移上，此位移可以通过求解基于静平衡状态给出的一个大型非线性系统方程获得。把节点位移与原节点位置相加即可得到变形后的网格。原始的弹簧模拟法鲁棒性较差，为此，文献[9-12]提出增加抗扭弹簧、张力弹簧等改进措施，取得了一些好的结果。弹簧模拟法的不足之处在于，这种方法不能对网格大小、形状进行控制，在位移较大的地方容易出现网格交错，产生负体积。对于结构网格，这种方法不能提供可靠的、自动保证网格正交性和光滑性的手段。弹性体法[13]是将计算域比拟为弹性体，通过求解弹性力学平衡方程得到新的网格坐标，这种方法具有较强的网格变形能力，但求解大型方程组带来较大计算量，限制了它的使用。Delaunay图映射法[14]是一种简单、高效、无迭代过程的网格变形技术，它是通过在物面与外边界之间生成Delaunay图，也称为背景网格的方式，将变形量映射到Delaunay背景网格上，然后再通过背景网格插值得到新网格坐标。该方法计算量小，在凸边界小变形问题中具有较大优势，但对于凹边界网格变形处理比较复杂，而且当背景网格出现交叉时，会导致此后的变形网格质量越来越差。RBF法[15-16]是近期变形网格方法研究的一个热点, 其原理是运用径向基函数对物面网格点的位移进行拟合，并以此来建立径向基函数插值模型，然后利用该模型计算空间网格节点变形量并更新网格坐标，此方法的优点是能够适应复杂外形的大变形问题，网格质量高。但在应用于大规模网格变形时，如果直接采用所有物面节点来构造径向基插值函数的话，也会存在计算量过大的问题。为此，一些作者提出采用数据精简的办法，例如，利用贪婪算法减少径向基基函数的方法[16-17]、峰值选点法[18]、RBF与Delaunay图映射法相结合的方法[19]等，以提高网格变形效率，取得了好的效果。

可以看到, 上述各种网格变形法几乎都是以代数或几何的方法去解决动网格生成问题，而不像在静态结构网格生成时那样，多以求解偏微分方程的方法来解决网格生成问题。这是否意味着偏微分方程法不适用于动网格生成呢，当然不是。事实上，早有研究人员尝试使用这种方法来解决网格变形问题，只是还没有获得实用的结果。例如，Johnson和Tezduyar[20]曾尝试通过求解由网格速度导出的线弹性泊松方程，以产生随流体边界或界面运动的网格。Ushijima[21]则尝试采用椭圆型网格生成法在每一时间步生成贴体网格，这种方法即便在较复杂的物体运动下也能生成高质量的网格。但在每一时间步求解泊松方程的计算代价实在太高，远远高于弹簧模拟法。为了减少求解泊松方程的计算开销，Grüber 和Cartens[22]采用了一种简化策略，不是在每一时间步，而是只在选定的有限步求解泊松方程。比如只在物体运动的最大和最小时间步两个位置给出由偏微分方程法生成的网格，其他中间位置的网格则用这两个网格进行插值。严格的讲，这种方法只能算是偏微分方程法与代数插值法的组合。

本文提出了一种基于求解椭圆型微分方程的网格变形法，可用于结构动网格的生成。数值算例表明，用该方法生成的动网格可以保持原始网格的光滑性、正交性，具有良好的鲁棒性和较高的计算效率。

1 一般椭圆型网格生成法

椭圆型网格生成法是静态结构网格生成中最受欢迎、使用最普遍的一种方法。这种方法可以生成高质量的网格，具有网格正交性可控、光滑的C2连续性、鲁棒性较好的优点。

用笛卡儿坐标r=[xyz]T代表物理空间Ω内的点，曲线坐标ρ=[ξηζ]T代表计算空间D内的点。计算空间内的点ρ与物理空间内的点r通过可逆映射T相联系：

T∶r→ρ⟹T-1∶ρ→r

因此，网格生成其实就是一种由式(1)向量值函数确定的映射：

(1)

映射函数被要求满足泊松方程：

(2)

式中：P、Q、R为强迫函数，也称为源项。

通过交换自变量(x，y，z)与因变量(ξ，η，ζ)的位置，在变换后的计算域D内，式(1)可写为

α1rξ ξ+α2rη η+α3rζζ+2(β12rξη+β23rη ζ+

β31rζ ξ)+J2(rξP+rηQ+rζR)=0

(3)

式中：

(4)

另有，物理空间边界∂D上的边界条件为

(5)

式(3)经整理后可写为

α1(rξ ξ+φPrξ)+α2(rηη+φQrη)+α3(rζζ+φRrζ)+

2(β12rξη+β23rη ζ+β31rζ ξ)=0

(6)

数值求解式(6)，方程的解即为笛卡儿坐标下的网格点。内点单元的大小和形状通过调节泊松方程中的源项来实现。显然，源项在生成适宜精确模拟流动问题的网格时起到了关键性的作用。

源项控制方法有多种形式，一般可分成2种类型：第1种类型是由式(6)导出源项表达式，然后在给定的坐标边界上以正交条件和指定的第一层网格与壁面的距离作为约束条件求出初始边界源项，初始内点源项则通过插值获得，然后迭代求解方程；第2种类型是在迭代过程中根据源项的变化情况，采用“人工”控制源项值，来得到所期望的方程。Thomas和Middlecoff[23]以及Sorenson 和 Steger[24-25]方法属于第一类，Hilgenstock[26]方法属于第二类。这些方法都要求将物面和外边界处的源项值内插至内部网格点处。

求解式(6)最实用的方法是高斯-赛德尔连续超松弛(Gauss-Seidel Successive Overrelaxation，GSSOR)迭代法。迭代循环一般由外迭代和内迭代组成，在外迭代循环中改变源项，在内迭代循环中源项保持不变。

另外，在结构网格生成中，还需要用到以下几个概念。譬如，一对一的映射也被称为规则映射。而采用规则映射得到的网格往往被称为规则网格。

从更实际的角度出发，规则网格也可以采用以下的定义[27]：如果一个网格的所有网格线均处于一个物理域的规定域内，并且同族网格线不相交，不同族的两条网格线如果相交，且只相交一次，则称这样的网格是规则的。

对于由式(6)生成的网格，可以做出以下论断[27]：假设已由式(6)生成了一个规则网格，那么以此网格为初始网格, 改变源项值，但保持新源项的符号分布与初始网格相同，则所生成的新网格也是规则的。

2 本文动网格生成算法

假设已有一个网格，准备用于带有移动边界的某非定常流动问题的模拟。该网格可以是由任何网格生成方法获得的网格，例如，该网格可以是由代数方法或者双曲方法生成，也可以直接由椭圆型网格生成法生成。以此网格为出发点构造后续每一时间步的动网格。

2.1 源项的获取

显然，无论采用何种方法生成的网格，它首先应该是规则的。否则，它将不适宜被用来模拟一个实际的流动问题。剩下的问题是，是否可以在泊松类方程的解空间中找到与此相同的一个网格。对这个问题先不做正面回答。观察式(6)的网格生成系统，可以看到其中所有导数项以及度量系数都可以由已知网格通过差分方法计算得到，系统中未知的是源项φP、φQ和φR。式(6)可以写为

α1rξφP+α2rηφQ+α3rζφR=-α1rξ ξ-

α2rηη-α3rζζ-2(β12rξη+β23rη ζ+β31rζξ)

(7)

如果令

b=-α1rξξ-α2rηη-α3rζζ-2(β12rξη+β23rη ζ+β31rζξ)

式(7)又可写成标量形式：

(8)

如果用i、j、k分别表示相应于坐标ξ、η、ζ方向的整型指标，式(7)和式(8)中的导数项可以表示成中心差分的形式：

(9)

其他2个坐标方向η和ζ的导数也有类似的表达。

式(8)的边界条件可以十分方便地采用有限体积法中的哑元(Dummy Cell)技术，或者称为虚拟单元技术来确定。简单说，就是通过在每个网格块的边界处附加另外的虚拟网格来实现,以i=Imax边界为例,第一层虚拟单元为

rImax+1,j,k=rImax,j,k+rImax,j,k-rImax-1,j,k

第二层虚拟单元可在第一层虚拟单元的基础上生成。当然，如果边界导数项采用单边差分算法求解的话，甚至可以不用虚拟单元，但是会造成编程上麻烦。

通过求解式(8)可以在初始网格的每一点上都获得一个源项φP、φQ和φR。同时，也可以看到，对于任意一个规则网格，有且只有唯一一个与式(8)的解相对应的源项。因此，可以作出以下论断：每一个规则网格都可视为是由具有某种源项分布的泊松类方程生成的结果。

显然，源项分布是椭圆型网格生成法中的一个关键因素。它包含了描述一个网格的重要信息，但又不是网格本身。观察物理域内的泊松方程式(2)，可以看到源项满足的是曲线坐标对笛卡儿坐标二阶导数项组成的方程，因此，源项的变化是通过影响这些导数项来改变网格位置的，不是直接对网格位置进行控制。网格本身也就是式(2)的离散解只有在每个坐标变换方向的离散点数、正交性要求、网格间距要求(可转化为对源项分布的要求)确定后，才能被确定。

通过上面的讨论，知道一个规则网格必有一个属于这个网格本身的源项分布。或者说，对任意一个规则网格总存在满足式(8)的一个源项分布。本文的重点在于如何在动网格生成中利用初始网格中已知的源项分布。

对某些移动边界问题，如飞机外挂物分离，虽然存在相对运动，但物面边界的外形是不变的。因此，如果采用一个物理空间内的解析函数f(x,y,z)来定义这个物面边界，并且假定它对应计算空间内的坐标面ξ=ξmin，如图1所示，尽管物面边界经过了位移和/或转动，但向量值函数f(x,y,z)仍然保持不变。这意味着此物面边界上的点在经历了从t0到t1的运动后，对曲线坐标的一阶和二阶以及交叉导数也没有变化。因此，有

(10)

(11)

式中：i，j=1，2，3；ξ1=ξ;ξ2=η;ξ3=ζ。

如前所述，在时间t=t0时刻的初始网格源项(或称为静态网格)可由式(8)确定。再考虑到式(10)和式(11)时，可知在ξ=ξmin边界面上的源项是时不变的。因此可以得到

(12)

式中：L=P，Q，R。事实上，式(12)对于其他外形不变的静止或运动的边界也同样成立。

除了物面边界和远场边界外，网格中一般还有内部边界，即块与块之间的对接面，或者是O形以及C形拓扑中具有相同的物理空间坐标，但分属不同的计算坐标面而具有边界含义的网格面(有的资料称为坐标切口, Coordinate Cut)。这种边界存在于网格内部，在新的网格生成过程中，不能证明式(12)的存在。但从源项控制的角度看，源项直接控制的是边界附近网格的角度和网格间距，如果只要求内部边界附近的网格角度和网格间距得到保持，而不关心内部边界的实际位置的话，完全可以在内部边界上人为指定式(12)成立。

这意味着，所有用于动网格生成的边界源项都可以通过继承静态网格的边界源项而被决定下来。至于内部源项，一般可以通过在同一时刻的两两边界之间插值获得。但在本文中，考虑到由式(8)获得的解不但包括了源项本身，而且也隐含地包括了它们的插值关系。因此，不需要重新插值，而是直接使用由式(8)获得的解作为动网格生成的源项。

2.2 动网格生成求解过程

为了更容易理解本文提出的动网格生成方法，有必要先回顾一下通常用于静态椭圆型网格生成的数值求解过程。这个迭代求解过程可以归纳为

1) 假定边界源项初值。

2) 改变边界源项值，即实施“源项控制”，使其更接近于对边界附近网格角度和间距控制所需要的源项值。

3) 将边界上新的源项值插值到内部网格点处。

4) 迭代求解出离散网格生成系统。

5) 重复步骤3)和4)，直到获得收敛解。

称步骤4)中的迭代过程为内迭代，而步骤2)～4)的迭代为外迭代。源项只在外迭代中改变，在内迭代中不变。

显然，在静态网格生成中，外迭代的主要功能是决定出源项值。在本文中，由于用于动网格生成的源项假定为从静态网格继承而来，因而是已知的。所以，新的动网格生成过程为

1) 读入已知静网格。

2) 通过求解式(8)，从已知静网格提取源项。

3) 在移动边界和其他边界上施加与静网格相同的网格点分布。

4) 采用上一时间步的网格坐标值作为当前时间步的初始解。

5) 迭代求解离散方程式(6)直至收敛。

在步骤5)的迭代收敛过程中，源项始终保持不变。因此不像静态网格生成中那样需要外迭代过程，这使得动态网格生成更加高效。

以上动网格生成过程，新的源项不仅继承了静态网格源项的符号和数值，也继承了对网格特性的控制。因此边界附近的网格正交性以及网格间距也得到了保持。

3 算法验证

3.1 带平移和转动的NACA0012翼型动网格生成

本算例中，采用NACA0012翼型模拟边界运动。所用静态网格是一个准二维单块结构网格，如图2所示,图中横纵坐标均为无量纲长度。该网格由大小为106×69的二维O型网格沿展向拉伸一个弦长得到，二维O型网格采用双曲型网格生成法生成，第一层网格与翼表面的距离为弦长的10-3。

本算例中，移动边界的位置随时间的变化关系为

(13)

(14)

式中：α(t)为绕1/4弦线的转动；α0为转动运动幅值；Tm为转动运动周期；xc(t)、yc(t)为翼型形心平动位移；hm为水平方向的最大位移。显然，这种平动与转动位移的组合运动可用于模拟飞机外挂物分离问题。

本算例的目的在于验证本文方法的正确性和收敛性，并对生成网格的质量进行评估。因此，只生成网格而不求解流场。

设α0=10.0°,Tm=0.15 s,hm=2.0c，c为翼型的无量纲弦长。在一个周期Tm内取300个点，因此时间步长Δt=0.000 5 s。

图3给出翼型运动过程中4个不同时刻翼面附近网格的局部视图，可以看到，翼面附近网格的质量得到了很好的保留。从定量角度，给出了最大长宽比和平均长宽比的迭代收敛历程，如图4和图5所示；并给出了网格单元角对正交角偏离的统计结果，如表1所示。由图4和图5可以看出，网格最大长宽比和平均长宽比的迭代过程是稳定的，在迭代600步时已收敛到与静态网格最大长宽比和平均长宽比相当的程度。从网格角度偏离的统计结果表1来看，内点以及外边界上的角度偏离比静态网格的角度偏离要大，但仍然处于可以接受的程度。情况最好的是翼型表面，几乎达到了和静态网格一样的程度。

GridInteriorOutboundarySurfaceofairfoilMaximumAverageMaximumAverageMaximumAverageStaticgrid(t=0s)12.82000.62180.78950.00870.00710.0002Dynamicgrid(t=0.1125s)15.16002.414011.45000.29490.00730.0002Dynamicgrid(t=0.15s)12.79000.89696.35800.11730.00710.0002

3.2 绕俯仰振荡翼型NACA0012的非定常无黏跨声速流动

本算例所用静态网格与3.1节所用静态网格完全一致。翼型绕1/4弦线的俯仰振荡由式(15)决定：

α(t)=αm+α0sin(ωt)

(15)

式中：α(t)为瞬时迎角；αm、α0和ω分别为平均迎角、迎角振荡幅度和角运动圆频率。角频率ω与减缩频率k的关系定义为

k=ωc/(2U∞)

(16)

式中：U∞为自由来流速度。

同时采用本文提出的动网格生成方法和线段弹簧(Segment Spring)模拟法对以下工况进行模拟。

工况AMa∞=0.76,αm=0.016°,α0=2.51°,k=0.081 4。

工况BMa∞=0.60,αm=3.16°,α0=4.59°,k=0.081 4。

图6～图9给出了瞬时升力系数CL和俯仰力矩系数Cm与试验结果[28-29]的比较。除了工况B中的Cm曲线外，2种动网格方法的预测结果与试验结果都具有很好的一致性。根据文献[28]的分析，Cm曲线的差异是由于实验模型的实际俯仰力矩中心与理论名义中心之间存在微小差异而引起的。对于实际力矩中心，文献[28]认为应该是27.3%而不是25%。图9也给出了采用实际力矩中心对数值结果进行修正后的情况，修正后的结果与实验结果具有更好的一致性。

为比较本文方法和线段弹簧模拟法在网格生成上的效率，采用工况B的运动参数让翼型做同样的俯仰振荡运动，但屏蔽流场求解代码，在每一时间步长上只生成网格，不求解流场。并指定2种方法的外迭代次数都为500次,内迭代次数都为200次。网格求解都采用Gauss-Seidel点松弛法。对2种算法花费的总CPU时间进行了统计,同时给出了每个网格点每次外迭代所用时间，结果如表2 所示。实测数据说明在同样的迭代次数下本文方法在计算效率上落后于弹簧模拟法。但需要注意,这是在指定相同迭代步数的前提下作出的比较，如果换为以某种网格质量指标作为约束条件的话，结果可能会有所不同。也就是说，如果是在达到同样网格质量的情况下来比较，椭圆型方法就可能因为需要的迭代步数更少，而表现出更高的计算效率。

表2 2种动网格生成方法执行时间比较(迭代步数=500,网格点数=34 845)

Table 2 Comparison of execution time for two dynamic grid generation methods (iterations=500, number of grid points=34 845)

MethodTotalcomputationtime/sGridpointiterationtime/sThepresent454.92.61×10-5Springanalogy317.51.82×10-5

为了研究方法的鲁棒性，选择工况B作为研究实例。让角运动幅度α0按照每次4.59° 的间隔递增，然后采用2种方法分别求解动网格，直到网格线出现交叉或负体积网格。当α0增加到7×4.59°=32.13° 时，弹簧模拟法首先出现了网格线交叉的情况，如图10所示。该图也同时展示了翼型前缘网格的变形情况。在同样条件下采用本文方法得到的网格如图11所示，可见，无论在前缘还是后缘附近，网格质量明显高于弹簧模拟法得到的结果。2种方法的迭代次数都是500次。这说明本文方法确实具有较高的鲁棒性，能够适应较大的边界运动。

3.3 再入返回舱的俯仰振荡

通过对再入返回舱(OREX)俯仰振荡的模拟，演示本文动网格生成方法处理三维多块结构网格的能力。

OREX的几何尺寸来自文献[30]。舱体附近的三维静态网格如图12所示，整个网格被划分为4块，各网格块的规模分别为51×21×21、51×21×21、67×41×51、67×41×51，其中第1指标指示法向，第2指标指示流向，第3指标指示周向。计算区域为直径的30倍，第1层网格距离壁面为直径的10-3。

算例采用无黏流模型，对返回舱施加式(15)形式的强迫运动，旋转轴为通过质心的z轴。根据Yoshinaga等[31]的实验结果，选择以下两组参数进行非定常流模拟：

1)Ma∞=0.9,αm=7.0°,α0=1.0°,7.0°，k=0.047 7。

2)Ma∞=1.0,αm=7.0°,α0=1.0°,7.0°，k=0.047 7。

图13是返回舱旋转7.0° 时在xOy平面内的网格。表3给出的是静导数Cmθ实验结果与计算结果的比较。考虑到实验结果没有明确给出平衡迎角αm，而且采用的是自由振动法，因此也没有对应的强迫振动迎角振幅α0以及减缩频率k，这几个参数是本文指定的, 可能与实验参数之间并没有准确的对应关系，表中的结果可用于定性比较。

此处，利用三维算例再次对网格更新时间进行了测算。同样，屏蔽流场求解，在每一时间步长上只生成网格，不求解流场。并指定本文算法和弹簧模拟算法的外迭代次数都为200次,内迭代次数仍然为200次。结果如表4所示，其中，每个网格点单次迭代的时间开销与表2的结果十分接近，验证了表2的结果。

图14是返回舱旋转45.0° 时在xOy平面内的网格。说明这种方法在三维大位移情况也能得到较好的结果。

表3 静导数Cmθ实验结果与计算结果的比较

Table 3 Comparison of computational and experimental results of static derivativeCmθ

NominalMachnumberα0/(°)StaticderivativeCmθ/rad-1ComputationExperiment[31]0.91.0-0.1437.0-0.127-0.1541.01.0-0.1137.0-0.098-0.148

表4 2种动网格生成方法执行时间比较(迭代步数=200,网格点数=325 176)

Table 4 Comparison of excution time for two dynamic grid generation methods (iterations=200, number of grid points=325 716)

MethodTotalcomputationtime/sGridpointiterationstime/sThepresent1783.62.73×10-5Springanalogy1345.52.06×10-6

3.4 黏性网格算例

原理上，使用式(8)提取静态网格源项的方法并不仅限于无黏网格，黏性网格同样适用。不同之处在于,黏性动网格的生成需要特别关注初始迭代向量,即每一时间步长下网格迭代初场的选取以及物面锐角边界的处理问题[32]。这是容易导致黏性动网格收敛失败的2个主要因素。第1个因素主要是由于黏性网格中高梯度区域的存在使求解网格的偏微分方程非线性程度加剧，收敛域缩小[27]而引起的。解决的办法一般是通过减小时间步长来控制物面位移量，采用第n步已收敛的网格作为第n+1步的迭代初场，尽量使迭代初场位于第n+1步网格的收敛域内。因此, 黏性动网格的生成往往需要考虑时间步长的约束,尤其是在物面存在凸锐角边界时这个问题会更突出。

此处，以一个带凸锐角边界的RAE2822翼型为例,验证本文方法对黏性网格的适应性。翼型绕1/4弦线的俯仰振荡关系仍然由式(15)决定。取振荡幅度α0=5°,振荡频率f=10 Hz(ω=2πf)。C形静态网格如图15和图16所示，网格数据来自文献[33]。

图17～图19显示了翼型绕1/4弦线旋转5°后前缘及后缘部分的局部网格与原始静态网格(图20和图21)相比，除后缘部分的网格质量有所下降外，其他物面网格质量仍得到了较好的保持。后缘部分网格质量下降的原因，主要是由于凸锐角边界(斜率间断)引起。而对间断较为敏感其实是微分方程法本身固有的弱点,这也是导致在间断附近数值解收敛十分困难的原因。

对这个问题，文献[32]建议采用人工固定凸锐角边界附近网格点的办法来解决。但由于动网格生成过程需要与流场求解器相融合，难以提供有效的人机交互接口来指定需要“固定”的网格点，因而难以实施。本文建议在原始静态网格生成时，对带有凸锐角的物面边界，首先应考虑生成具有O型拓扑的网格，并对凸锐角采取“打钝”的办法，预先用1～2 mm左右的圆弧对锐角进行圆角化，避免凸锐角的出现。对只能采用C型拓扑的网格，在凸锐角边界附近无法采用钝化处理。可考虑采用其他策略，例如，考虑在迭代初期采用最简单的线段弹簧法生成动网格“初场”，然后由椭圆型方程方法继续进行迭代，获得更加光顺的网格。弹簧法与椭圆型方程法的迭代步数可按3∶7 进行分配。图22是采用以上策略后得到的后缘部分的网格，与图19相比，网格质量得到了明显改善，甚至部分恢复到了与初始静态网格(图21)相当的程度。

4 结 论

1) 基于求解椭圆型偏微分方程方法，实现了一种简单高效的结构动网格重构方法，数值模拟结果证明了方法的正确性。

2) 在同样迭代步数约束下，本文方法花费的网格求解时间略长于传统弹簧模拟法，但网格质量较好，方法的鲁棒性优于弹簧模拟法。

3) 文中采用继承静态网格源项，并在动网格生成过程中保持不变的策略，既解决了动网格生成过程中如何确定网格源项的问题，又节约了搜索源项需要的外迭代过程，是方法成功的关键。

[1] STEGER J L, DOUGHERTY F C, BENEK J A. A chimera grid scheme[C]//Mini Symposium on Advances in Grid Generation， 1982.

[2] CHAN W M, PANDYAY S A, ROGERS S E. Effcient creation of overset grid hole boundaries and effects of their locations on aerodynamic loads[C]//21st AIAA Computational Fluid Dynamics Conference. Reston: AIAA, 2013.

[3] KIM N, CHAN W M. Automation of hole-cutting for overset grids using the x-rays approach: AIAA-2011-3052[R]. Reston: AIAA, 2011.

[4] 王文, 阎超. 新型重叠网格洞面优化方法及其应用[J]. 航空学报, 2016, 37(3): 826-835. WANG W, YAN C. Novel overlapping optimization algorithm of overlapping grid and its applications[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(3):826-835 (in Chinese).

[5] 张伟伟, 高传强, 叶正寅. 气动弹性计算中网格变形方法研究进展[J]. 航空学报, 2014, 35(2): 303-319. ZHANG W W, GAO C Q, YE Z Y. Reseach progress on mesh deformation in computational aeroelasticity[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(3): 303-319 (in Chinese).

[6] FOSTER N F. Accuracy of high-order cfd and overset interpolation in finite volumee/difference codes[C]//22nd AIAA Computational Fluid Dynamics Conference. Reston: AIAA, 2015.

[7] UCHIDA Y, KANDA H. Spacing control of 3-D transfinite interpolation grid generation: AIAA-2005-5243[R]. Reston: AIAA, 2005.

[8] BATINA J T. Unsteady Euler airfoil solutions using unstructured dynamic meshes[J]. AIAA Journal, 1990, 28(8): 1381-1388.

[9] BLOM F J. Considerations on the spring analogy[J]. International Journal of Numerical Methods in Fluid, 2000, 32(6): 647-668.

[10] CANTARITI D L, GRIBBEN W M, BADCOCK K J, et al. A grid deformation technique for unsteady flow computations[J]. International Journal of Numerical Methods in Fluids, 2000, 32(3): 285-311.

[11] 刘枫. 动网格技术研究及其在高超声速流动中的应用[D]. 长沙：国防科学技术大学， 2009: 54-55. LIU F. Investigations of dynamic grid generation and its applications for supersonic flow[D]. Changsha: National University of Defense Technology, 2009:54-55 (in Chinese).

[12] 郭正, 刘君, 翟章华. 用非结构动网格方法模拟有相对运动的多体绕流[J]. 空气动力学学报, 2001, 19(3): 310-316. GUO Z, LIU J, ZHAI Z H. Simulation of flows past multi-body in relative motion with dynamic unstructured method[J]. Acta Aerodynamic Sinica, 2001, 19(3): 310-316 (in Chinese).

[13] TEZDUYAR T E. Stability finite element formulations for incompressible flow computations[J]. Advances in Applied Mechanics, 1992, 28(1): 1-44.

[14] LIU X, QIN N, XIA H. Fast dynamic grid deformation based on delaunay graph mapping[J]. Journal of Computational Physics, 2006, 211(2): 405-423.

[15] BORE A, SCHOOT M S, FACULTY S H. Mesh deformation based on radial basis function interpolation[J]. Computers and Structures, 2007, 85(11): 784-795.

[16] RENDALL T C S, ALLEN C B. Efficient mesh motion using radial basis functions with data reduction algorithms[J]. Journal of Computational Physics, 2010, 229(8): 2810-2820.

[17] WANG G, HARIS H M, YE Z Y. An improved point selection method for hybrid unstructured mesh deformation using radial basis functions: AIAA-2013-3076[R]. Reston: AIAA, 2013.

[18] 魏其, 李春娜, 谷良贤, 等. 一种基于径向函数和峰值选择法的高效网格变形技术[J]. 航空学报, 2016, 37(7): 1401-1414. WEI Q, LI C N, GU L X, et al. An efficient mesh deformation method based on radial basis functions and peak-selection method[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(7): 1401-1414 (in Chinese).

[19] 孙岩, 邓小刚, 王光学, 等. 基于径向基函数改进的Delaunay图映射动网格方法[J]. 航空学报, 2014, 35(3): 727-735. SUN Y, DENG X G, WANG G X, et al. An improvement on Delaunay graph mapping dynamic grid method based on radial basis functions[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(3): 727-735 (in Chinese).

[20] JOHNSON A A, TEZDUYAR T E. Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 119(1-2): 73-94.

[21] USHIJIMA S. Three-dimensional arbitrary Lagrangian-Eulerian numerical prediction method for non-linear free surface oscillation[J]. International Journal for Numerical Methods in Fluids, 1998, 26(5): 605-623.

[22] GRÜBER B, CARTENS V. Computation of the unsteady transonic flow in harmonically oscillating turbine cascades taking into account viscous effects[J]. Journal of Turbomachinery, 1998, 120(1): 104-111.

[23] THOMAS P D, MIDDLECOFF J F. Direct control of the grid point distribution in meshes generated by elliptic equations[J]. AIAA Journal, 1979, 18(6): 652-656.

[24] SORENSON R L, STEGER J L. Numerical generation of 2D grids by the use of poisson equations with grids control: N81-14722[R]. 1981.

[25] SORENSON R L. A computer program to generate 2D grids about airfoil and other shapes by the use of poisson’s equation: N80-26266[R]. 1980.

[26] HILGENSTOCK A. A fast method for the elliptic generation of three-dimensional grids with full boundary control[C]//Numerical Grid Generation in Computational Fluid Mechanics, 1988: 137-146.

[27] SONAR T. Grid generation using elliptic partial differential equations: DFVLR-FB89-15[R]. 1989.

[28] GAITONDE A L, FIDDES S P. A moving mesh system for the calculation of unsteady flows: AIAA-1993-0641[R]. Reston: AIAA, 1993.

[29] LI J, HUANG S. Unsteady viscous flow simulations by a fully implicit method with deforming grid: AIAA-2005-1221[R]. Reston: AIAA, 2005.

[30] UNMEEL B M. Synthesis of contributed simulations for OREX test cases: NASA/TM-1998-112238[R]. Washington, D.C.: NASA, 1998.

[31] YOSHINAGA T, TATE A, WATANABE M. Dynamic test of the orbital reentry vehicle (OREX) in a transonic wind tunnel with comparison to flight data: AIAA-1995-1901-CP[R]. Reston: AIAA, 1995.

[32] THOMPSON J F. Numerical grid generation-foundations and applications[M]. North Holland: Amsterdam, 1985.

[33] SLATER J W. RAE2822 transonic airfoil: Study #1[EB/OL]. (2015-01-22)[2016-07-17]. http://www.grc.nasa.gov/WWW/wind/valid/raetaf/raetaf.html.

(责任编辑：李明敏)

URL：www.cnki.net/kcms/detail/11.1929.V.20161124.1034.002.html

＊Corresponding author. E-mail: lufengling126@126.com

A dynamic structured grid generation method based onsolving elliptic equations

LU Fengling1,2,*, CHEN Xiaoqian1, YU Caihui2, MIAO Meng2

1.CollegeofAerospaceScienceandEngineering,NationalUniversityofDefenceTechnology,Changsha410073,China2.ScienceandTechnologyonSpacePhysicsLaboratory,Beijing100076,China

A simple robust structured dynamic grid generation method based on the solution of elliptic partial differential equations is developed for computing the unsteady flows with moving boundaries. In the method, the source terms of the Poisson equation which can control the spacing and the orthogonality of the grid are inherited from the known static grid, and held fixed throughout the process of dynamic grid generation. With the process, the outer iterations for determining the source terms usually needed in the elliptic grid generation can thus be saved, and no adjustable parameters are required to be prescribed. This makes the method more efficient and easy to implement in an existing CFD code. The numerical results demonstrate that the proposed method can provide an efficient way of deforming the grid based on solving elliptic partial differential equations. The orthogonality and smoothness of original static grid can be maintained well by the proposed method. When the same number of iterations are given as the constraint condition, the grid generation efficiency of the method is lower than that of the spring analogy, but the robustness of the method is superior to the spring analogy.

dynamic grid generation; elliptic grid generation; computational fluid dynamics; numerical grid generation; grid deformation approach

2016-07-21； Revised：2016-08-11； Accepted：2016-11-20； Published online：2016-11-24 10:34

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0303

2016-07-21； 退修日期：2016-08-11； 录用日期：2016-11-20； 网络出版时间：2016-11-24 10:34

www.cnki.net/kcms/detail/11.1929.V.20161124.1034.002.html

＊通讯作者.E-mail: lufengling126@126.com

卢凤翎， 陈小前， 禹彩辉， 等． 一种基于求解椭圆型方程的结构动网格生成方法[J]． 航空学报, 2017, 38(3): 120632. LU F L, CHEN X Q, YU C H, et al. A dynamic structured grid generation method based on solving elliptic equations[J]. Acta Aeronau-tica et Astronautica Sinica, 2017, 38(3): 120632.

V211.3

A

1000-6893(2017)03-120632-14