CFD中统计误差的数值精度分析

闵耀兵，马燕凯，李松

中国空气动力研究与发展中心，绵阳 621000

计算流体力学(CFD)是一门利用数值方法模拟流动问题的科学，其所依赖的数值方法在理论设计时需满足一定的精度要求，而理论设计的精度在CFD应用中不一定能够达到，因此经常需要对算法的理论设计精度进行数值验证[1-6]。算法精度的数值验证过程多以加密网格的方式进行，即通过逐步加密网格，使得数值误差逐渐减小，由网格加密前后数值误差的比值和网格尺度的比值可以计算出算法的数值精度[1-6]。一般情况下，相比于选取某固定点的数值误差而言，对整个计算域内的数值误差进行统计赋范更具有代表性和普适性。通常情况下，数值误差的统计方式一般可取其L1范数、L2范数和L∞范数，并认为统计误差的各范数在描述算法的数值精度上具有等价性。因此，在CFD中对算法进行数值精度验证时，有些研究人员仅考察统计误差的L1范数[6-7]或者L2范数[8-10]，也有同时考察L1范数和L∞范数的[4-5]以及同时考察L2范数和L∞范数的[11]，当然也有更为严谨的研究者同时基于L1范数、L2范数和L∞范数等3种常用范数进行数值精度分析的[1-2,12]，却鲜有文献仅参考L∞范数。

然而，在对某具体的CFD算法的理论设计精度进行数值验证时，由于CFD算法的理论设计多基于光滑流场假设，在遭遇流场间断(如激波和接触间断等)[8]、网格不连续(如拐折和突然拉伸等)[12-13]等问题时，算法的数值精度可能会存在不同程度的降阶问题[8,11-13]。即便是在光滑流场中，在极值点附近采用非线性加权插值也可能会产生降阶问题[14]。诸如此类的降阶问题通常只会出现于计算区域中的某个局部，不会导致整个计算区域的数值精度全部降阶[13]。此时统计误差的各范数在数值精度上一般不再具有等价性，各范数的数值精度与其赋范方式有关，不同赋范方式可能会给出不同的数值精度[13]，且一般情况下L1范数的数值精度最高，而L∞范数的数值精度最低，其中的原因和关系本文将给出详细分析。

针对统计误差各范数的数值精度不一致的问题[13]，本文通过对统计误差的具体赋范方式进行理论分析，指出：当且仅当整个计算区域内各点的数值精度完全一致时，统计误差各范数的数值精度才相等，此时统计误差的各范数在数值精度上具有等价性；而当全场的数值精度不完全一致时，统计误差的各范数的数值精度并不一致，其具体精度与赋范方式有关，且一般情况下表现为L1范数的数值精度最高，而L∞范数的数值精度最低。本文的理论分析结果很好地解释了统计误差各范数的数值精度之间的关系，同时也为CFD算法精度验证时的误差统计方式提供了理论参考。

在理论分析的基础上，本文分别针对全场各点精度不完全一致和完全一致的情形进行了验证，验证的结果与理论分析结论相符。然后对极值点附近的非线性加权插值引起的降阶问题进行了具体分析，利用本文的理论分析结论，从数值精度上较好地解释了非线性加权插值在极值点附近的降阶问题。

1 统计误差的数值精度分析

当计算区域中离散的网格尺度为h时，记整个离散区域中误差统计的总点数为Nh，则随着计算网格的加密，网格尺度h减小而总点数Nh增大，即Nh与1/h成正变关系。将每个点的数值误差ei简单表示为

ei=cihrii=1,2,…,Nh

(1)

式中：ri为第i点的数值精度，一般情况下ri均为整数(∀i)；ci为与网格尺度h无关的误差常数，一般可以表示为求解变量及其导数的函数。

(2)

(3)

(4)

1.1 各点精度不完全一致的情形

在实际CFD计算中，由于离散边界、流场间断以及网格不连续等问题的存在，经常会出现各点数值精度ri并不完全相同的情况，即有

min(ri)

(5)

为便于分析，不妨设整个计算区域中总共存在n个数值精度(当各点数值精度不完全相同时，n≥2)，其数值由小到大分别记为R1，R2，…，Rn，即

R1=min(ri)

(6)

一般情况下，随着计算网格的加密，精度值Rj以及其数量n均保持不变。不妨记数值精度为Ri的数值误差的点数为Ni，则有

(7)

(8)

不妨记

(9)

由于ci为与网格尺度h无关的常系数，故随着计算网格的加密，Cj也基本保持不变。由式(9)有

(10)

在网格加密过程中，不妨设网格尺度h的加密倍数rfc为

(11)

式中：下标c和f分别对应粗网格和细网格，网格的加密过程使得rfc>1。当网格尺度加密倍数为rfc时，整个计算区域中误差统计的总点数Nh的变化关系为

(12)

式中：dN为计算网格的加密维数，dN=1,2,3分别对应一维、二维和三维网格。同时记数值精度为Rj阶的误差点数Nj随网格的加密满足：

(13)

由于式(7)及Nj

dj≤dNj=1,2,…,n

(14)

随着计算网格的加密，误差统计的总点数Nh要比各精度的点数Nj增加得更快。

根据定义(12)及(13)，容易得到

(15)

(16)

为便于分析，不妨记

(17)

进一步容易得到

(18)

(19)

(20)

将式(18)代入(20)中，并整理得到

(21)

由上述分析以及网格尺度关系式(11)，容易得到

(22)

由于dN≤3以及式(14)，容易知道dj≤3(∀j)，特别地，对于本文数值算例中考虑的二维问题(dN=2)有dj≤2，又m≥1，则式(22)中的数值精度Om在绝大多数情况下可以表示为

(23)

当m趋于正无穷大时，有

(24)

考虑到一般情况下式(14)及式(23)成立，故有

O1≥O2≥…≥O∞

(25)

即当各点计算精度不完全相同时，其统计误差L1范数的数值精度O1最高，L2范数的数值精度O2次之，以此类推，L∞范数的数值精度O∞最低。

1.2 各点精度一致的情形

当全场精度均一致时，有

min(ri)=max(ri)=R

(26)

类似于式(8)中的分析，容易得到

(27)

式中：

(28)

类似于Cj、Ch也为与网格尺度h无关的常数。此时统计误差Lm范数的数值精度Om可以表示为

(29)

考虑到网格间距关系式(11)容易得到

Om=R

(30)

与各点计算精度不完全相同时的情形不一样，当各点计算精度均相同时，其统计误差的各范数均为R阶精度，此时统计误差的各范数之间在数值精度上具有等价性。

2 数值算例

针对上述关于统计误差数值精度的理论分析，设计相应的数值试验进行验证。类似于Mao等[13]的做法，设计不同连续程度的网格，计算各点的几何守恒律误差并进行误差统计，通过逐步加密网格以得到统计误差的各范数的数值精度表现。

2.1 网格生成

本文中采用的周期性二维计算网格基于如下算法解析生成：

doj=1,Nj

doi=1,Ni

if(|x0|<0.4) then

y2=sinκ(x2)

else

y2=0

end if

xi,j=x0

yi,j=y0+A·y1·y2

end do

end do

(31)

式中：Ni和Nj分别为计算坐标下ξ方向和η方向的网格分布点数，实际网格生成过程中取Ni=Nj=N+1；A为波动幅值(A=0时为均匀直角网格)，在本文中均取A=0.2；函数y2的指数κ取不同值(在本文中指数κ取1,2,3)时可以得到在x=±0.4处不同连续程度的网格，函数y2描述的网格线如图1所示。

图1 函数y2分布

2.2 数值离散

二维网格中几何守恒律误差可以表示为(以Ix为例)

(32)

式中：

(33)

(34)

将式(34)的差分算子δ均取为如下的二阶精度中心差分格式：

(35)

(36)

2.3 结果分析

(37)

(38)

同理有

(39)

而当κ=3时，由于R1=Rn，则由式(26)对应的分析易知

(40)

按照算法(31)生成的网格均满足d1=1，下面本文考虑另外一种网格中统计误差的数值精度问题。取κ=3，在算法(31)的基础上，按照如下方法修正x=±0.4、y=0处的网格分布：

表1 几何守恒误差的统计精度(κ=1)

表2 几何守恒误差的统计精度(κ=2)

表3 几何守恒误差的统计精度(κ=3)

if(|x0|=0.4&y0=0) then

end if

(41)

式中：γ为网格间距指数；B为与网格间距无关的常数，其取值方式为

(42)

式中：参数B的选取使得在x=±0.4、y=0处的网格波动在加密过程中不至于被机器字长所湮没(双精度浮点运算)。按照算法(41)生成的计算网格在间断附近的分布如图2所示。

图2 局部网格示意图(γ=1)

表4 几何守恒误差的统计精度(γ=1)

(43)

(44)

(45)

(46)

基于算法(41)生成的网格的数值试验结果见表4～表6。当γ=1和γ=2时，表4和表5中所示的统计误差的各范数对应的数值精度分别与表达式(44)和式(45)中的完全相同。当γ=3时，表6中给出的统计误差的各范数的数值精度与理论分析(46)相同，值得注意的是表6中L1范数的数值精度为二阶，与理论分析(22)相同，而式(23)此时给出的数值精度应为三阶，此时式(22)和式(23)给出的数值精度并不相同。

2.4 算例应用

针对WENO格式的非线性加权插值[4]在极值点附近可能会降阶的问题，Henrick等[14]设计了相应的数值算例进行验证，并指出WENO格式[4]中为避免分母为零而引入的小量ε的取值对数值精度的影响：ε取值相对较大时(如ε=10-6)，非线性权近似于其线性最优权，WENO格式在极值点附近一般不会降阶；而ε取值较小时(如ε=10-40)，非线性权远远偏离其线性最优权，WENO格式在极值点附近会存在降阶问题。针对目前已得到广泛应用且同样采用非线性加权思想构造的WCNS-E-5格式[3,18-19]，有必要测试其在极值点附近的数值精度情况[18]，为了准确反映WCNS-E-5格式[3]的非线性权在极值点附近的特性，取ε=10-300(四倍精度浮点运算)。

依据算法(31)取A=0生成均匀网格，参照Henrick等[14]的测试方法给定周期性初值：

(47)

表5 几何守恒误差的统计精度(γ=2)

表6 几何守恒误差的统计精度(γ=3)

采用WCNS-E-5格式[3]离散ux并将其相对于解析值的数值误差统计于表7中,表中的统计误差的数值精度满足：

(48)

值得指出的是，Henrick等[14]的数值实验中得出的统计误差的数值精度与式(48)所述一致，只是Henrick等[14]只给出了m=1,2,∞这3种情况下的精度，且没有对统计误差的L1、L2和L∞范数所表现的不同的数值精度给出解释与说明。

由统计误差的数值精度式(22)和式(23)易知R1=3，即统计区域中的最低数值精度为三阶。对于周期性初值式(47)而言，在极值点附近u′=0、u″≠0、u‴≠0，Henrick等[14]的分析指出此时非线性加权的WENO格式只有三阶精度，则采用类似非线性加权方式构造的WCNS格式在此类极值点附近也应仅具有三阶精度[18]，与表7中的最低数值精度为三阶相吻合。

由于本算例基于二维网格，故dN=2，由式(22)易知d1=1，即在极值点附近WCNS格式[3]的降阶以条带状形式出现，条带的宽度与格式的计算模板有关，但条带宽度不影响统计误差的数值精度。

表7 统计误差的数值精度

3 结 论

针对数值误差统计时采用的不同赋范方式可能表现出不同数值精度的问题，从理论上详细分析了统计误差赋范方式对其数值精度的影响，并给出了一般规律：当且仅当整个计算区域内的数值精度完全一致时，统计误差各范数的数值精度才相等，此时统计误差的各范数在数值精度上才具有等价性；而当全场的数值精度不完全一致时，统计误差各范数的数值精度并不一致，其具体精度与赋范方式有关，且一般情况下表现为L1范数的数值精度最高，L2范数的数值精度次之，以此类推，L∞范数的数值精度最低。

基于统计误差数值精度的理论分析结论，本文有针对性地设计了相应的数值试验予以验证，对于全场精度完全一致和不完全一致的情况，数值试验均给出了与理论分析完全吻合的结果，表明本文的研究结论能够为CFD中算法的精度验证提供理论依据。