有限差分法的坐标变换诱导误差

2021-07-07刘君魏雁昕韩芳

航空学报 2021年6期

刘君，魏雁昕，韩芳

大连理工大学航空航天学院，大连 116024

有限差分法应用于具有复杂外形的网格时需要进行坐标变换，即采用数学方法将流体力学控制方程由笛卡尔坐标系变换到与物体形状相吻合的贴体坐标系(Body-Fitted Coordinates，BFC)以后再进行格式离散。在实际工程问题中，很少能直接给出坐标变换的解析表达式，因此，网格生成技术成为计算流体力学(Computational Fluid Dynamics，CFD)的一个重要分支。离散是对连续的近似表达，在生成的网格点基础上计算坐标变换系数存在误差。例如，一维函数u=u(x)从物理空间变换到计算空间ξ=ξ(x)，导数ux变成ux=ξxuξ，可以看出，如果坐标变换系数ξx采用数值算法得到，则必然也会引入误差。算例表明，即使坐标变换系数ξx采用解析值，误差依然存在，本文把这种由坐标变换引起的误差称为坐标变换诱导误差。

早在1961年Trulio和Trigger[1]就注意到了这个问题，提出了一维问题中的“体积守恒(conservation of volume)”概念。1978年Steger等[2-4]研究了从守恒型控制方程出发计算具有复杂几何外形的流场时坐标变换系数的影响，即使是均匀流场，采用传统坐标变换系数求解方法也会产生非零源项，该源项由坐标变换引起，在流场中造成非物理解，严重时会引起数值振荡，影响计算稳定性[5-7]。1979年Thomas和Lombard[8-9]在前人研究基础上提出了守恒型的坐标变换系数计算方法，并正式提出了被后来广泛引用的几何守恒律(Geometric Conservation Law，GCL)概念。根据计算过程中网格是否随时间变形，GCL被进一步区分为“体积守恒律”(Volume Conservation Law，VCL)和“面积守恒律”(Surface Conservation Law，SCL)，前者解决“任意拉格朗日-欧拉”(Arbitrary Lagrange-Euler，ALE)方程动网格问题，应用范围远没有后者普遍，目前在高精度格式领域研究GCL常指后者。本文后续提到的GCL也指“面积守恒律”。

按照文献[9]中网格导数变换理论，如果使用二阶精度格式离散uξ，那么采用中心差分格式离散ξx不影响导数ux的精度。基于上述认识，对GCL研究工作基本停顿下来，随着高阶精度格式的发展与应用这个问题又引起重视。从1999年Gaitonde和Visbal[10]提出了紧致格式的GCL以来，近20年来对基于高精度有限差分的GCL的研究得到了充分发展，其中2011年邓小刚等证明了满足GCL的充分条件[11]并提出了坐标变换系数的守恒计算方法——C MM(Conservative Metric Method)算法；随后，在2013年进一步发展出了坐标变换系数计算结果唯一的对称守恒计算方法——SC MM算法(Symmetrical Conservative Metric Method)[12]，该研究成果对于紧致类格式具有普适性指导意义。近年来科研工作者们对构建WENO格式的GCL时常借鉴其思想，如将WENO格式分解成中心差分部分及数值耗散部分[13-15]，中心部分直接采用SC MM计算。

文献[15]的观点“面积守恒律的简化过程中，用到两个数学前提，即网格物理坐标对计算坐标连续可导和混合导数都连续。但是在实际计算过程中，不能保证计算网格坐标点等间距或二阶可导，简化过程中所做假设不成立，并且网格导数由数值差分格式计算得到，因此引入数值误差”在GCL研究领域具有普遍性。为了验证上述结论，针对坐标变换函数的各阶导数和混合导数均存在且连续的若干算例进行考查。

1 数值现象

目前常见的高精度格式大多包含随着流场变化的模板或限制器，对于参数变化的流场难以区分误差来源是差分格式还是几何诱导。因此，选择超声速均匀流动，来流马赫数M∞=2，无量纲流动参数(ρ,u,v,p)=(1.4,2,0,1)，其中ρ、u、v、p分别为无量纲密度、X方向速度、Y方向速度和压力。从二维Euler方程出发进行数值模拟，首先给出柱坐标算例。

1.1 柱坐标均匀网格的坐标变换诱导误差

由柱坐标和直角坐标之间的变换函数及其逆函数可以写出解析表达式：

(1)

式中：r为圆环区域半径；x为X轴方向坐标，x=rcosθ；y为Y轴方向坐标，y=rcosθ；θ为径向网格线与X轴正向夹角。

计算区域为r∈[1,2]的圆环，沿径向r和周向θ网格均匀分布，径向网格点总数Nξ=51，周向网格点总数Nη=(360/7.2)+Nη0，其中Nη0为θ=0°网格线沿圆周方向编号。为便于处理周向边界条件，进行计算域的拓展，例如，一阶迎风格式沿周向拓展2个点，编码η=Nη0=1和η=51 重合于x轴(即θ=0°网格线)，编码η=0和η=50 对应θ=-7.2°的网格线，编码η=2和η=52均是θ=7.2°的网格线，对于七点WENO格式需拓展更多的点，此时Nη0=3，但是网格夹角保持不变。

基于以上设定条件，计算空间坐标和柱坐标之间的线性关系为

(2)

计算中用到的5个坐标变换系数解析求出：

(3)

对流项计算采用一阶迎风离散格式及Van Leer通量分裂格式，计算收敛以后数值解和理论值的误差就是坐标变换诱导误差，其流动参数误差的分布规律类似，此处以密度误差为例进行分析。当流场径向内外圆边界根据Riemann不变量处理时，其密度误差δρ分布如图1(a)所示，由于超声速流动扰动不会向上游传播，流场中x<0区域(左侧半圆)的计算结果定性合理，但是在x>0 区域存在明显的边界影响。采用固定边界值的方法处理径向内外圆边界也可以得到收敛解，其密度误差分布云图如图1(b)所示，可以看出，流场参数误差分布左右两侧不对称特性依然存在。当径向外圆边界采用Riemann不变量处理时，内圆边界把x<0出口值作为在y相等高度上x>0对应点的入口值，其计算结果如图1(c)所示，可以看出流场参数误差在x=0左右两侧呈现对称分布。由此可以推断：图1(a)和图1(b)中流场密度误差分布出现左右两侧不对称的原因是在x>0内环进口边界处理中没有考虑坐标变换引起的误差。

图1 3种边界条件下密度误差云图

定量比较图1密度误差特性，可以看出圆柱上游误差分布相等且不受内圆边界处理影响，因此，仅对x<0区域的密度误差分布进行分析。为了后续研究，坐标变换系数除了采用式(3)求导的理论值，还采用二阶中心格式的数值解进行计算，图2为不同的系数计算方法和不同的差分格式下密度误差δρ比较，为便于对比两种格式，误差云图选取相同幅值。从结果看，WENO格式的误差比一阶迎风格式更大。

图2 柱坐标系密度误差云图

同心圆网格算例表明进出口边界变换成曲线以后，坐标变换也有影响，在超声速流动稳定性、气动声学、湍流转捩等研究中，上游微小的扰动会导致结果出现差异，数值模拟过程中应该分辨和评估这种进口边界误差的不均匀特性对结果造成的干扰。为了排除进口边界的影响，进行1.2节的数值实验。

1.2 正弦函数网格的坐标变换诱导误差

设计了如图3(a)所示的拼接网格，左侧第1块网格(x≤0.6)完全是直角坐标系下的均匀网格：

(4)

式中：0≤ξ≤30, 0≤η≤100。

第2块网格(x>0.6)按照如式(5)所示的变换函数生成：

(5)

式中：31≤ξ≤130, 0≤η≤100。

计算中用到的坐标变换系数可以解析求出，5个系数中有2个为常数：

(6)

另外3个系数J、ξx、ξy的分布如图3所示，JACB代表雅可比行列式，可以看出，在拼接线x=0.6处整体来说存在明显不连续，但是在y=-1,0,1坐标处系数基本上光滑过渡。

图3 正弦网格及其变换系数

采用Steger通量分裂格式，不同离散格式和坐标系数计算方法得到的密度误差特性如图4所示，可以看出边界处理对计算结果没有影响。同样WENO格式的计算误差比一阶迎风格式更大，其原因将在后文中讨论，这里主要关注坐标系数精度的影响。

图2和图4的计算结果均表明坐标系数计算方法影响坐标变换诱导误差，采用二阶中心格式数值解的误差比根据式(3)解析求得的误差小两个数量级。

图4 不同坐标变换系数下正弦网格密度误差云图

GCL研究的理论基础如前所述，从Steger等研究GCL开始[2-4]，就是针对差分格式构建更高精度的坐标系数计算方法，从而实现不影响方程整体离散精度的目标，例如，按照满足SC MM算法计算ξx，对于五阶精度WCNS-E-5格式，需要六阶中心插值[16]，同样对于WENO格式ξx也是形式复杂的高精度插值[13-15]。原则上坐标系数计算方法的精度极限是准确的理论值，图2和图4 密度误差云图中数值解的误差比解析解的还要大，这种现象与GCL的研究目标似乎相矛盾。

从同心圆算例和正弦网格算例计算结果可以看出，在“网格物理坐标对计算坐标连续可导和混合导数都连续”这两个数学前提成立条件下，即使坐标系数采用理论值误差依然存在，由此看来，认为坐标变换诱导误差来自变换导数不连续的观点不准确，不是坐标系数计算方法及其精度可以解决的问题。

目前建立的GCL算法针对不同的对流项离散格式，还没有考虑通量分裂格式这个因素。图5 是一阶迎风格式条件下坐标变换系数解析求解采用AUSM、HLLC、ROE、Van Leer 4种通量分裂格式的压力误差δp云图。采用WENO格式得到的误差更大，分布类似，这里不再显示。

从图5中可以看出，不同通量分裂格式的坐标变换诱导误差存在差异，研究GCL时应该考虑这个因素。GCL把坐标变换诱导误差归因于网格，但是调研未看到网格质量因素如何影响的定量研究文献。网格生成技术研究领域经常采用正交性、光滑性和合理分布性评价网格质量[17]，前述柱坐标算例不但变换函数的导数连续，而且均匀网格完全正交且光滑，依然存在坐标变换诱导误差，说明这些网格质量评价因素未能反映问题的本质，设计1.3节的算例对这个问题进行数值实验。

图5 不同通量分裂格式的压力误差

1.3 网格正交性和光滑性的影响

在柱坐标结果(图1和图2)分析中，注意到误差分布具有轴对称特征，考察计算中用到的5个坐标变换系数分布，只有雅可比行列式J是轴对称的。在正弦函数网格中沿着y=-1,0,1这3条线点J变化较小，对应的误差也较小。为考察J和坐标变换诱导误差之间的关系，设计和计算了正交网格和非正交网格。对于网格线与坐标轴平行的完全正交网格，坐标系数ξy=ηx=0，所有算例中的J、ξx和ηy同时不为常数，图6(a)和图6(b)为指数加密网格及J云图。可以看出存在明显的不均匀特性。除了常用的指数加密，还计算了Δxi+1=αΔxi和Δyi+1=βΔyi等比例加密，其中Δxi和Δyi为i点处网格尺度，α和β为参数，计算网格如图6(c)所示，计算中坐标变换系数采用二阶中心差分的计算值，该网格下J云图如图6(d)所示。图6(a)和图6(c)两种网格下的误差云图如图7所示，说明单纯J变化不一定产生坐标变换诱导误差。由于相邻网格之间的网格尺度Δx和Δy均不相等，从网格质量评价因素看光滑性不够好，因此选取α=β=2.00得到的网格光滑性较差，且网格节点数减少为图6(a)和图6(c)中网格的1/10，该网格实际应用很少采用，计算结果也没有发现误差，至少表明对均匀流动来说，光滑性不会引起坐标变换诱导误差。

将图6(a)中的直角网格变成如图8(a)所示的倾斜网格，在不光滑网格基础上考查正交性的影响。网格线倾斜以后除ηx以外其余4个坐标系数均不为常数，定义网格y轴与x轴负向的夹角为θ′，例如θ′=45°工况下J、ξx和ξy分布如图8(b)～图8(d)所示。在0°<θ′<180°范围内，选择夹角参数θ′=m×15°(m=1,2,…,11)进行数值实验，图9给出了夹角参数θ′=30°、θ′=45°和θ′=120° 的计算结果，所有算例的误差云图和图7相同，均没有坐标变换诱导误差，说明在光滑性较差的情况下，夹角保持常数时正交性较差的网格也不会产生误差。

图6 正交网格参数

图7 3种网格误差云图

图8 45°夹角倾斜网格及其坐标系数分布

图9 不同夹角倾斜网格密度误差云图

通过以上数值现象，可以得出：

1) 坐标变换诱导误差来自“网格物理坐标对计算坐标的变换函数的导数不连续”的观点不准确。

2) 不能仅采用网格的正交性和光滑性评估坐标变换诱导误差的影响特性。

2 坐标变换诱导误差的机理

将笛卡尔坐标系下(t,x,y,z)的三维守恒型Euler方程变换到贴体坐标系(τ,ξ,η,ζ)，其中t为时间，τ为贴体坐标系下的时间，ζ为三维空间贴体坐标下与ξ、η轴成右手系的坐标。其中得到原始形式：

UIt+FIx+GIy+HIz=S

(7)

如果不考虑“体积守恒律”，那么It=0。“面积守恒律”方程为

(8)

式中：ξz、ηz、ζx、ζy和ζz均为坐标变换度量系数。以Ix=0为例，根据逆变换关系式得

(9)

式中：yζ、zξ、zη和zζ均为逆变换度量系数。继续求导，以式(9)中第1项为例，得

y″η ξzζ+yηz″ζ ξ-z″η ξyζ-zηy″ζ ξ

(10)

式中：上标″表示求二阶偏导数。

原则上在给定网格点(x,y,z)ijk位置以后，可以按照式(3)计算ξx，但是由于存在计算误差，导致Ix=0不成立，即使对均匀来流也存在引起质量不守恒的现象，这就是产生坐标变换诱导误差的机理。

根据式(10)可以看出，即使yη和ξx采用准确的理论值，由于还有∂/∂ξ的离散误差，恒等式也是不成立的，依然存在坐标变换诱导误差，可以合理解释1.1节和1.2节算例采用理论值也出现坐标变换诱导误差的数值现象。坐标变换系数采用数值解抹平了yη和ξx的变化斜率，对应的∂/∂ξ也较小，形成的源项S也较小，由此可以解释图2和图4中出现的二阶中心格式数值解的诱导误差比准确的理论值还要小的数值现象。同样，考查1.3节网格正交性和光滑性的算例，如果在离散以后恒等式还成立，那么也没有误差。

认识到坐标变换诱导误差的机理和坐标变换系数的理论值无法消除误差，那么就通过构建yη、ξx和∂/∂ξ相互匹配的算法使在离散空间上恒等式Ix=Iy=Iz=0重新成立，这就是GCL的出发点。由于yη和ξx只能计算得到，在GCL研究领域有个无法摆脱的逻辑困境：坐标变换系数采用解析值是不满足GCL算法的。

3 几何守恒律算法

从量纲角度分析，yη等参数对应矢量，ξx等参数是矢量面积，J-1对应体积，如果采用二阶中心差商计算矢量，涉及离散点(x,y,z)ijk及其相邻网格点正好形成六面体，在网格分布非等距情况下，采用不同算法得到的yη结果不同，网格点不共面情况下矢量面积ξx也与算法相关，得到的体积也不唯一，因此，早期CFD应用中Steger等解释Ix≠0的原因就是矢量、面积和体积的算法不相容，采用一种加权平均算法计算坐标变换系数，解决了二阶中心格式离散流动方程时遇到的问题[2-4]。后来研究者发现，如果流动方程的离散格式非二阶中心格式，即使采用Steger等几何相容满足Ix=0等面积守恒律，依然存在均匀来流质量不守恒的现象，认识到计算式(10)中y″ηξ等二阶混合导数与流动方程的离散格式相关，不同格式需要相适应GCL算法，其中邓小刚等[12]提出的SC MM实际上是一个普适性的准则，应用该准则能够满足面积守恒律要求，即Ix=Iy=Iz=0。

按照SC MM准则，如果计算过程中流动方程算子发生变化，坐标变换系数的算子也要做相应调整。例如对于一阶迎风格式，在时间迭代过程中流动算子随着参数动态变化，满足SC MM要求的坐标变换系数算子也要具有迎风特性，且每步需要重新计算，这给具体应用带来困难。目前SC MM主要应用于格式模板系数确定的格式，从给出的均匀流动计算看，对于可以写成若干中心格式组合的高阶精度紧致类格式效果明显，对于WENO格式不能完全消除坐标变换诱导误差[14]，推测与其中的迎风型模板有关。

目前的GCL算法仅与格式相关，如果再将通量分裂格式、限制器等因素增加进来，沿着这条技术路线消除坐标变换诱导误差构建的算法会非常复杂，预测将给应用带来困难。

4 离散等价方程算法

除了第3节中所述GCL外，另一种思路是2018年刘君等提出离散等价方程算法[18]，放弃对恒等式Ix=Iy=Iz=0的约束条件，在离散对流项的同时也离散S≠0源项，以此消除坐标变换诱导误差。下面对该算法进行简单介绍。

坐标变换系数在以上三维守恒型Euler方程通量中是线性的，引入

(11)

对流项通量和源项表示为

(12)

(13)

(14)

根据Steger和Van Leer通量分裂格式的通用表达式，可以直接证明坐标变换系数在通量分裂以后也是线性的，算子分别作用后

(15)

同样使用算子δ1离散源项：

(16)

得到半离散形式：

(17)

(18)

对于Steger和Van Leer通量分裂格式可以推导出如式(18)所示的离散方程的算子形式，由于AUSM、HLLC和Roe通量分裂格式是在半点处重构出来的，直接推导难度较大，采用数值实验验证，计算结果也没有误差。图10和图11是采用离散等价方程算法后第1节部分算例的计算结果，可以看出离散等价方程算法完全消除了坐标变换诱导误差。

图10 柱坐标系密度误差云图(一阶迎风)

图11 采用离散等价方程算法后的正弦网格密度误差云图

采用算子形式证明离散等价方程算法适用于任意格式，近期比较一阶迎风格式和五阶ENO类格式[20]，计算结果和1.1节、1.2节定性一致，高精度格式的坐标变换诱导误差更大，通过分析误差来源，得出的结论是ENO模板涉及的周围相邻网格点较多，在源项较大的局部区域引入的误差较大，结论与本文一致。

式(11)～式(18)从带有S≠0源项的离散等价方程算法推导中，算子作用在对流项守恒通量上，有些高精度格式采用原始变量重构界面的形式，讨论这种基于原始变量的离散等价方程算法实现。

5 非均匀网格界面通量重构

以密度为例进行讨论，笛卡尔坐标系下均匀网格的MUSCL格式为

(19)

(20)

对应二阶精度的3点中心格式：

(21)

根据Taylor级数：

(22)

(23)

式中：Δx1和Δx2分别为i点与i+1点网格间距和i点与i-1点网格间距；O(·)代表截断误差。

可以看出在非均匀网格Δx1≠Δx2条件下，应用MUSCL的3点中心格式式(21)达不到二阶计算精度，在该种网格下二阶精度格式表示为

(24)

分析计算空间内3点中心格式的精度：

(25)

不论坐标变换函数如何变化，其实质是物理空间到网格编码的映射关系。以一维空间进行简要说明，通过变换函数ξ=ξ(x)把物理空间x∈[a,b]映射到计算空间ξ∈[ξ(xa),ξ(xb)]，其中xa和xb分别为物理空间中a、b两点坐标，把ξ空间离散成N个网格点，格点间距为常数：

Δξ=[ξ(xa)-ξ(xb)]/(N-1)

(26)

计算中网格编码i和ξ之间存在线性关系i=1+[ξ(xi)-ξ(xb)]/Δξ，常数Δξ合并至坐标变换系数中，因此在离散后本质还是物理空间x∈[a,b]到i∈[1,N]的映射关系。实际应用中很少引入上述中间变换，而是直接从x∈[a,b]映射到ξ∈[1,N]，在计算空间离散点上满足ξi=i。在离散后的计算空间上Δξ=1，决定计算空间大小的唯一参数是网格节点数N，如果网格加密过程中节点总数N保持不变，实际上只是调整离散点在物理空间[a,b]分布，与计算空间的映射关系没有任何变化。如果通过增加网格总数N加密网格，即使物理区域[a,b]保持不变，其在计算空间的映射关系也发生相应的变化。

根据以上分析可知通过式(25)难以直接得到截断误差量级，严格的精度分析还需要变换至物理空间，由Taylor级数可得

(27)

由于Ji+1同样需要Taylor展开，具体推导过程非常复杂，因此针对一维空间变换函数ξ=ξ(x)进行定性的量级分析。

根据微分关系式dξ=ξxdx可以得到坐标变换系数的量级ξx～O(Δx-1)，进一步推导出高阶导数量级：

(28)

式中：f代表任意变量，用于任意物理量的高阶导数量级分析。

在空间多维情况下，J至少与两个变量有关，证明过程较为复杂，按照有限体积法的量级定性分析。MUSCL类型的格式应用于有限体积法考虑了度量系数，常用形式为

(29)