文章快速检索  
  高级检索
CFD中统计误差的数值精度分析
闵耀兵, 马燕凯, 李松     
中国空气动力研究与发展中心, 绵阳 621000
摘要: 在计算流体力学(CFD)的算法研究中经常会对离散误差进行数值精度分析,通常以统计误差的各范数为研究对象,最常用的统计误差范数为L1范数、L2范数和L范数,一般认为各范数在数值精度上具有等价性。实际上,由于流场局部存在间断、网格局部不连续或者是在极值点附近采用非线性加权插值等可能使数值方法存在局部降阶问题,导致统计误差各范数所表达的数值精度并不一致。通过详细的理论分析,揭示了统计误差各范数所表达的数值精度之间的关系,并通过相应的数值试验予以验证。研究结果不仅能够指导CFD算法的数值精度验证工作,而且也可为更为复杂流动模拟的数值精度判定提供理论依据。
关键词: 统计误差    误差范数    数值精度    精度分析    降阶    
Accuracy analysis of numerical error with statistical forms in CFD
MIN Yaobing, MA Yankai, LI Song     
China Aerodynamics Research and Development Center, Mianyang 621000, China
Abstract: Generally, the accuracy of an algorithm should be validated numerically in Computational Fluid Dynamics(CFD), and the research object is composed by statistical norms of numerical error, which is usually represented by the L1 norm, the L2 norm, and the L norm based on the hypothesis that each norm is equivalent in accuracy order. In reality, there are locally discontinuities of flow variables, non-smoothness of mesh, and nonlinear interpolations near critical points that will result in the loss in accuracy of numerical algorithm, causing different numerical orders of accuracy of each norm of numerical error. By carrying out detailed theoretical analysis, the relationship of different error norms in accuracy order is exhibited in this paper and is soon validated by numerical experiments. The research results in this paper can not only serve as a guide to the validation of the accuracy order in a CFD algorithm, but also theoretically support the judgement of numerical simulation order with more complex flow.
Keywords: statistical error    norm of error    numerical accuracy    analysis of accuracy    loss in accuracy    

计算流体力学(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简单表示为

$ {e_i} = {c_i}{h^{{r_i}}}\;\;\;\;i = 1,2, \cdots ,{N_h} $ (1)
 

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

如果记误差向量eh=[e1, e2, …, eNh],则对整个计算区域的数值误差ei进行统计,相当于对误差向量eh取范数。对于向量eh定义Lmh

$ L_m^h \buildrel \Delta \over = {\left( {\frac{1}{{{N_h}}}\sum\limits_{i = 1}^{{N_h}} {{{\left| {{e_i}} \right|}^m}} } \right)^{\frac{1}{m}}}\;\;\;m \ge 1 $ (2)
 

式中:$ \buildrel \Delta \over = $为定义符号。根据赋范线性空间理论中关于向量范数的定义,易知式(2)中Lmh为误差向量eh的范数,不妨将Lmh记为误差向量ehLm范数,即

$ {\left\| {{\mathit{\boldsymbol{e}}_h}} \right\|_{{L_m}}} \buildrel \Delta \over = L_m^h $ (3)
 

随着计算网格的加密,网格尺度h逐步减小,而误差统计的总点数Nh逐渐增加。在计算网格的逐步加密过程中,不同网格尺度h下的误差向量eh的维数Nh也逐步增加,固定维数的赋范线性空间理论中向量范数的等价性定理不再成立,不能据此证明各范数Lmh之间(m取不同值时)的等价性。

根据误差向量ehLm范数Lmh的定义(3),在CFD中数值验证时经常采用的L1范数、L2范数和L范数可以分别表示为

$ \left\{ \begin{array}{l} L_1^h \buildrel \Delta \over = {\left\| {{\mathit{\boldsymbol{e}}_h}} \right\|_{{L_1}}} = \frac{1}{{{N_h}}}\sum\limits_{i = 1}^{{N_h}} {\left| {{e_i}} \right|} \\ L_2^h \buildrel \Delta \over = {\left\| {{\mathit{\boldsymbol{e}}_h}} \right\|_{{L_2}}} = \sqrt {\frac{1}{{{N_h}}}\sum\limits_{i = 1}^{{N_h}} {{{\left| {{e_i}} \right|}^2}} } \\ L_\infty ^h \buildrel \Delta \over = {\left\| {{\mathit{\boldsymbol{e}}_h}} \right\|_{{L_\infty }}} = \max \left( {\left| {{e_i}} \right|} \right) \end{array} \right. $ (4)
 

下面给出误差向量ehLm范数Lmh的数值精度Om的详细分析过程。

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

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

$ \min\left( {{r_i}} \right) < \max\left( {{x_i}} \right) $ (5)
 

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

$ {R_1} = \min \left( {{r_i}} \right) < {R_2} < \cdots < {R_n} = \max \left( {{r_i}} \right) $ (6)
 

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

$ \sum\limits_{i = 1}^n {{N_i}} = {N_h} $ (7)
 

则误差向量ehLm范数Lmh可以改写为

$ L_m^h = {\left( {\frac{1}{{{N_h}}}\sum\limits_{j = 1}^n {{h^{m{R_j}}}} \sum\limits_{i = 1}^{{N_j}} {{{\left| {{c_i}} \right|}^m}} } \right)^{\frac{1}{m}}} $ (8)
 

不妨记

$ {C_j} = {\left( {\frac{1}{{{N_j}}}\sum\limits_{i = 1}^{{N_j}} {{{\left| {{c_i}} \right|}^m}} } \right)^{\frac{1}{m}}}\;\;\;j = 1,2, \cdots ,n $ (9)
 

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

$ L_m^h = {\left( {\frac{1}{{{N_h}}}\sum\limits_{j = 1}^n {{h^{m{R_j}}}} {N_j}C_j^m} \right)^{\frac{1}{m}}} $ (10)
 

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

$ {r_{{\rm{fc}}}} = \frac{{{h_{\rm{c}}}}}{{{h_{\rm{f}}}}} $ (11)
 

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

$ {N_{{h_{\rm{f}}}}} = r_{{\rm{f}}{{\rm{c}}^N}}^d \cdot {N_{{h_{\rm{c}}}}} $ (12)
 

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

$ {N_{j,{\rm{f}}}} = r_{{\rm{f}}{{\rm{c}}^j}}^d \cdot {N_{j,{\rm{c}}}}\;\;\;\;j = 1,2, \cdots ,n $ (13)
 

由于式(7)及Nj < Nh对于任意网格尺度h均满足,容易得知在计算网格的加密过程中,Nj随网格尺度加密倍数rfc增长的速度必然不会大于总点数Nh增加的速度,即

$ {d_j} \le {d_N}\;\;\;j = 1,2, \cdots ,n $ (14)
 

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

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

$ {N_h} = \frac{{N_h^0}}{{{h^{{d_N}}}}},{N_j} = \frac{{N_j^0}}{{{h^{{d_j}}}}}\;\;\;\;j = 1,2, \cdots ,n $ (15)
 

式中:Nh0Nj0分别为最粗网格的网格总点数和精度为Rj的网格点数,与当前网格尺度h无关,则有

$ L_m^h = {\left( {{h^{{d_N}}}\sum\limits_{j = 1}^n {{h^{m{R_j} - {d_j}}}} \frac{{N_j^0}}{{N_h^0}}C_j^m} \right)^{\frac{1}{m}}} $ (16)
 

为便于分析,不妨记

$ {{\bar R}_{j,m}} = {R_j} - \frac{{{d_j}}}{m} $ (17)
 

进一步容易得到

$ \begin{array}{l} L_m^h = {\left( {{h^{{d_N}}}\sum\limits_{j = 1}^n {{h^{m{{\bar R}_{j,m}}}}} \frac{{N_j^0}}{{N_h^0}}C_j^m} \right)^{\frac{1}{m}}} = \\ \;\;\;\;\;{h^{\min\left( {{{\bar R}_{j,m}}} \right) + \frac{{{d_N}}}{m}}}{\left( {\sum\limits_{j = 1}^n {{h^{m\left( {{{\bar R}_{j,m}} - \min\left( {{{\bar R}_{j,m}}} \right)} \right)}}} \frac{{N_j^0}}{{N_h^0}}C_j^m} \right)^{\frac{1}{m}}} \end{array} $ (18)
 

由于Nh0Nj0Cj均为与网格尺度h无关的常数,且必然存在j使得Rj, m-min(Rj, m)=0,以及

$ {{\bar R}_{j,m}} = \min \left( {{{\bar R}_{j,m}}} \right) \ge 0 $ (19)
 

则对于j=1, 2, …, n$ {h^{m\left( {{{\bar R}_{j, m}} - \min \left( {{{\bar R}_{j, m}}} \right)} \right)}}\frac{{N_j^0}}{{N_h^0}}C_j^m $中必然包括与网格尺度h无关的常数项以及与h相关的小量。

根据数值精度的定义,误差向量ehLm范数Lmh的数值精度Om可以表示为

$ {O_m} = \text{lg}{_{\frac{{{h_{\rm{c}}}}}{{{h_{\rm{f}}}}}}}\frac{{L_m^{{h_{\rm{c}}}}}}{{L_m^{{h_{\rm{f}}}}}} = \frac{1}{{\text{ln}{r_{{\rm{fc}}}}}}\text{ln}\frac{{L_m^{{h_{\rm{c}}}}}}{{L_m^{{h_{\rm{f}}}}}} $ (20)
 

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

$ \begin{array}{l} {O_m} = \left( {\min\left( {{{\bar R}_{j,m}}} \right) + \frac{{{d_N}}}{m}} \right)\frac{1}{{\text{ln}{r_{{\rm{fc}}}}}}\text{ln}\frac{{{h_{\rm{c}}}}}{{{h_{\rm{f}}}}} + \\ \;\;\;\;\;\;\;\frac{1}{{m \cdot \text{ln}{r_{{\rm{fc}}}}}}\text{ln}\frac{{\sum\limits_{j = 1}^n {h_{\rm{c}}^{m\left( {{{\bar R}_{j,m}} - \min\left( {{{\bar R}_{j,m}}} \right)} \right)}} \frac{{N_j^0}}{{N_h^0}}C_j^m}}{{\sum\limits_{j = 1}^n {h_{\rm{f}}^{m\left( {{{\bar R}_{j,m}} - \min\left( {{{\bar R}_{j,m}}} \right)} \right)}} \frac{{N_j^0}}{{N_h^0}}C_j^m}} \end{array} $ (21)
 

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

$ {O_m} = \min \left( {{R_j} - \frac{{{d_j}}}{m}} \right) + \frac{{{d_N}}}{m} $ (22)
 

由式(22)可以看出,当整个区域中各点计算精度不完全相同时,随着计算网格的加密,其统计误差的Lm范数Lmh的数值精度Om不仅与全场的数值精度Rj及其误差点数Nj随网格加密倍数rfc增加的指数dj有关,还与误差统计的总点数Nh随网格加密倍数rfc增加的指数dN相关,甚至与误差统计时范数选取的方式m也相关。

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

$ {O_m} = {R_1} + \frac{{{d_N} - {d_1}}}{m} $ (23)
 

m趋于正无穷大时,有

$ {O_\infty } = \mathop {\min}\limits_{m \to + \infty } \left( {{O_m}} \right) = \min\left( {{R_j}} \right) = {R_1} $ (24)
 

而根据式(4)中Lh的定义,随着计算网格的进一步加密,其统计误差的L范数Lh应为计算区域中数值精度最低点的误差值,同时对应的数值精度O也应为全场最低精度R1=min(ri),与式(22)和式(23)在m趋于正无穷大时给出的精度值相同,故式(22)和式(23)中给出的统计误差的数值精度关系对于L范数也同样适用。

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

$ {O_1} \ge {O_2} \ge \cdots \ge {O_\infty } $ (25)
 

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

1.2 各点精度一致的情形

当全场精度均一致时,有

$ \min \left( {{r_i}} \right) = \max \left( {{r_i}} \right) = R $ (26)
 

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

$ {L_m} = {\left( {\frac{1}{{{N_h}}}{h^{mR}}{N_h}C_h^m} \right)^{\frac{1}{m}}} = {h^R}{C_h} $ (27)
 

式中:

$ {C_h} = {\left( {\frac{1}{{{N_h}}}\sum\limits_{i = 1}^{{N_h}} {{{\left| {{c_i}} \right|}^m}} } \right)^{\frac{1}{m}}} $ (28)
 

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

$ {O_m} = \frac{1}{{\text{ln}{r_{{\rm{fc}}}}}}\text{ln}\frac{{{L_{m,{\rm{c}}}}}}{{{L_{m,{\rm{f}}}}}} = \frac{R}{{\text{ln}{r_{{\rm{fc}}}}}}\text{ln}\frac{{{h_{\rm{c}}}}}{{{h_{\rm{f}}}}} $ (29)
 

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

$ {O_m} = R $ (30)
 

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

2 数值算例

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

2.1 网格生成

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

$ \begin{array}{l} {\rm{do}}\;\;\;\;j = 1,{N_j}\\ {\rm{do}}\;\;\;\;i = 1,{N_i}\\ \;\;\;\;\;\;\;\Delta x = \frac{1}{{{N_i} - 1}},\Delta y = \frac{1}{{{N_j} - 1}}\\ \;\;\;\;\;\;\;{x_0} = i \cdot \Delta x - \frac{1}{2},{y_0} = j \cdot \Delta y - \frac{1}{2}\\ \;\;\;\;\;\;\;{x_1} = {\rm{ \mathsf{ π} }}{y_0} + \frac{{\rm{ \mathsf{ π} }}}{2},{y_1} = {\sin ^2}\left( {{x_1}} \right)\\ \;\;\;\;\;\;\;{\rm{if}}\left( {\left| {{x_0}} \right| < 0.4} \right)\;\;\;\;{\rm{then}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{x_2} = \frac{{5{\rm{ \mathsf{ π} }}}}{4}{x_2} + \frac{{\rm{ \mathsf{ π} }}}{2}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{y_2} = {\sin ^\kappa }\left( {{x_2}} \right)\\ \;\;\;\;\;\;\;{\rm{else}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{y_2} = 0\\ \;\;\;\;\;\;\;{\rm{end}}\;{\rm{if}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{x_{i,j}} = {x_0}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{y_{i,j}} = {y_0} + A \cdot {y_1} \cdot {y_2}\\ \;\;\;\;\;\;\;{\rm{end}}\;{\rm{do}}\\ \;\;\;\;\;\;\;{\rm{end}}\;{\rm{do}} \end{array} $ (31)
 

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

图 1 函数y2分布 Fig. 1 Distribution of function y2
2.2 数值离散

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

$ {I_x} = \frac{{\partial {{\hat \xi }_x}}}{{\partial \xi }} + \frac{{\partial {{\hat \eta }_x}}}{{\partial \eta }} $ (32)
 

式中:

$ \left\{ {\begin{array}{*{20}{l}} {{{\hat \xi }_x} = {J^{ - 1}}{\xi _x} = {y_\eta }}\\ {{{\hat \eta }_x} = {J^{ - 1}}{\eta _x} = - {y_\xi }} \end{array}} \right. $ (33)
 

几何守恒律误差Ix在有限差分离散算子δ作用下的差分形式IxN可以表示为

$ I_x^N = {\rm{ \mathsf{ δ} }}_\xi ^1\left( {{{\hat \xi }_x}} \right) + \delta _\eta ^1\left( {{{\hat \eta }_x}} \right) = \delta _\xi ^1\left( {\delta _\eta ^2y} \right) - \delta _\eta ^1\left( {\delta _\xi ^2y} \right) $ (34)
 

式中:下标表示求导方向,上标表示求导顺序。根据Deng等[15]的分析,当δ1=δ2时,式(34)中IxN=0;而当δ1δ2时,IxN一般不为零,而是以一定的精度趋于零[16-17]

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

$ \begin{array}{*{20}{c}} {{{\left( {{\delta _\xi }E} \right)}_j} = \frac{\alpha }{{2\Delta \xi }}\left( {{E_{j + 1}} - {E_{j - 1}}} \right) + }\\ {\frac{{1 - \alpha }}{{4\Delta \xi }}\left( {{E_{j + 2}} - {E_{j - 2}}} \right)} \end{array} $ (35)
 

对式(35)进行Taylor展开易知:当α=4/3时,差分格式(35)具有四阶精度;当α≠4/3时格式仅为二阶精度。为便于计算几何守恒律误差IxN,在实际计算过程中取δ1δ2,且分别为

$ \left\{ {\begin{array}{*{20}{l}} {{\delta ^1}:{\alpha _1} = \frac{1}{3}}\\ {{\delta ^2}:{\alpha _2} = \frac{2}{3}} \end{array}} \right. $ (36)
 
2.3 结果分析

根据式(12)中的定义,对于本文采用的二维计算网格有dN=2。按照算法(31)生成的计算网格在x=±0.4处存在间断,导致IxNx=±0.4附近的数值精度最低,则根据定义式(13),容易得到d1=1。根据Mao等[13]的分析,在远离网格间断的光滑区域,几何守恒律误差IxN均应以二阶精度趋于零,即Rn=2;而在x=±0.4附近,由于网格的不连续导致IxN出现一定程度的降阶,具体表现为

$ \left\{ {\begin{array}{*{20}{l}} {\kappa = 1}&{{R_1} = 0}\\ {\kappa = 2}&{{R_1} = 1}\\ {\kappa = 3}&{{R_1} = 2} \end{array}} \right. $ (37)
 

κ=1时,几何守恒律误差IxN在整个计算区域上统计误差的Lm范数的数值精度Omκ=1

$ O_m^{\kappa = 1} = \frac{1}{m} $ (38)
 

同理有

$ O_m^{\kappa = 2} = 1 + \frac{1}{m} $ (39)
 

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

$ O_m^{\kappa = 3} = 2 $ (40)
 

数值试验结果见表 1~表 3,数值试验结果表明:当κ=1和κ=2时,几何守恒律误差IxN在整个计算区域统计误差的的各范数对应的数值精度分别与理论分析式(38)和式(39)中给出的完全相同;当κ=3时,统计误差的各范数在数值精度上具有等价性,均应为二阶精度,恰如式(40)所述,同时也与表 3中数值试验的统计误差精度完全相同。

表 1 几何守恒误差的统计精度(κ=1) Table 1 Statistical accuracy of geometry conservation error(κ=1)
NL1L2L3L4L5L
误差精度误差精度误差精度 误差精度误差精度误差精度
200 2.67E-03 2.04E-02 4.18E-02 6.03E-02 7.55E-02 2.06E-01
300 1.77E-03 1.01 1.67E-02 0.49 3.65E-02 0.33 5.46E-02 0.25 6.97E-02 0.20 2.06E-01 0
500 1.06E-03 1.01 1.30E-02 0.50 3.09E-02 0.33 4.81E-02 0.25 6.30E-02 0.20 2.06E-01 0
800 6.58E-04 1.01 1.03E-02 0.50 2.64E-02 0.33 4.28E-02 0.25 5.74E-02 0.20 2.06E-01 0
1 200 4.38E-04 1.00 8.39E-03 0.50 2.31E-02 0.33 3.86E-02 0.25 5.29E-02 0.20 2.06E-01 0
2 000 2.62E-04 1.00 6.50E-03 0.50 1.95E-02 0.33 3.40E-02 0.25 4.78E-02 0.20 2.06E-01 0
3 000 1.75E-04 1.00 5.31E-03 0.50 1.70E-02 0.33 3.07E-02 0.25 4.41E-02 0.20 2.06E-01 0
5 000 1.05E-04 1.00 4.11E-03 0.50 1.43E-02 0.33 2.71E-02 0.25 3.98E-02 0.20 2.06E-01 0
8 000 6.55E-05 1.00 3.25E-03 0.50 1.23E-02 0.33 2.41E-02 0.25 3.62E-02 0.20 2.06E-01 0
表 2 几何守恒误差的统计精度(κ=2) Table 2 Statistical accuracy of geometry conservation error(κ=2)
NL1L2L3L4L5L
误差精度误差精度误差精度误差精度误差精度误差精度
200 1.75E-04 7.03E-04 1.40E-03 2.05E-03 2.61E-03 8.07E-03
300 7.80E-05 1.99 3.82E-04 1.50 8.18E-04 1.33 1.24E-03 1.25 1.61E-03 1.20 5.38E-03 1.00
500 2.82E-05 1.99 1.77E-04 1.50 4.14E-04 1.33 6.54E-04 1.25 8.72E-04 1.20 3.23E-03 1.00
800 1.10E-05 2.00 8.75E-05 1.50 2.22E-04 1.33 3.64E-04 1.25 4.96E-04 1.20 2.02E-03 1.00
1 200 4.90E-06 2.00 4.76E-05 1.50 1.29E-04 1.33 2.19E-04 1.25 3.05E-04 1.20 1.35E-03 1.00
2 000 1.77E-06 2.00 2.21E-05 1.50 6.53E-05 1.33 1.16E-04 1.25 1.65E-04 1.20 8.07E-04 1.00
3 000 7.85E-07 2.00 1.20E-05 1.50 3.81E-05 1.33 6.97E-05 1.25 1.02E-04 1.20 5.38E-04 1.00
5 000 2.83E-07 2.00 5.60E-06 1.50 1.93E-05 1.33 3.68E-05 1.25 5.51E-05 1.20 3.23E-04 1.00
8 000 1.10E-07 2.00 2.76E-06 1.50 1.03E-05 1.33 2.05E-05 1.25 3.13E-05 1.20 2.02E-04 1.00
表 3 几何守恒误差的统计精度(κ=3) Table 3 Statistical accuracy of geometry conservation error(κ=3)
NL1L2L3 L4L5L
误差精度误差精度误差精度误差精度误差精度误差精度
200 2.49E-04 3.48E-04 4.13E-04 4.62E-04 5.02E-04 9.45E-04
300 1.11E-04 1.99 1.55E-04 1.99 1.85E-04 1.99 2.07E-04 1.99 2.24E-04 1.99 4.22E-04 1.99
500 4.01E-05 1.99 5.61E-05 1.99 6.67E-05 1.99 7.47E-05 1.99 8.11E-05 1.99 1.52E-04 2.00
800 1.57E-05 2.00 2.20E-05 2.00 2.61E-05 2.00 2.92E-05 2.00 3.18E-05 2.00 5.94E-05 2.00
1 200 6.98E-06 2.00 9.77E-06 2.00 1.16E-05 2.00 1.30E-05 2.00 1.41E-05 2.00 2.64E-05 2.00
2 000 2.52E-06 2.00 3.52E-06 2.00 4.19E-06 2.00 4.69E-06 2.00 5.09E-06 2.00 9.51E-06 2.00
3 000 1.12E-06 2.00 1.56E-06 2.00 1.86E-06 2.00 2.09E-06 2.00 2.27E-06 2.00 4.23E-06 2.00
5 000 4.03E-07 2.00 5.64E-07 2.00 6.70E-07 2.00 7.51E-07 2.00 8.16E-07 2.00 1.52E-06 2.00
8 000 1.57E-07 2.00 2.20E-07 2.00 2.62E-07 2.00 2.93E-07 2.00 3.19E-07 2.00 5.95E-07 2.00

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

$ \begin{array}{l} {\rm{if}}\left( {\left| {{x_0}} \right| = 0.4\& {y_0} = 0} \right)\;\;\;\;\;{\rm{then}}\\ \;\;\;\;\;\;{y_{i,j}} = B \cdot {\left( {\Delta x \cdot \Delta y} \right)^{\frac{\gamma }{2}}}\\ {\rm{end}}\;{\rm{if}} \end{array} $ (41)
 

式中:γ为网格间距指数;B为与网格间距无关的常数,其取值方式为

$ \left\{ {\begin{array}{*{20}{l}} {B = 0.3}&{\gamma = 1}\\ {B = 100}&{\gamma = 2}\\ {B = 20000}&{\gamma = 3} \end{array}} \right. $ (42)
 

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

图 2 局部网格示意图(γ=1) Fig. 2 Local enlarged image of mesh (γ=1)

算法(41)生成的网格仅在点x=±0.4、y=0处存在间断,导致IxNx=±0.4、y=0附近的数值精度降到最低,则根据定义式(13),容易得到d1=0。除了x=±0.4、y=0点及其附近区域外,IxN均以二阶精度趋于零。而IxN在间断点附近的精度,相当于在其光滑网格的分析基础上直接加上算法(41)中的关于网格间距的波动项,考虑到IxN的计算过程式(34)需要进行两次求导,则算法(41)中指数γx=±0.4、y=0点附近的最低精度R1之间存在如下关系:

$ \left\{ {\begin{array}{*{20}{l}} {\gamma = 1}&{{R_1} = - 1}\\ {\gamma = 2}&{{R_1} = 0}\\ {\gamma = 3}&{{R_1} = 1} \end{array}} \right. $ (43)
 

由上述分析以及dN=2易知,当γ=1时,式(22)和式(23)均给出:几何守恒律误差IxN在整个计算区域上统计误差的Lm范数的数值精度Omγ=1

$ O_m^{\gamma = 1} = \frac{2}{m} - 1 $ (44)
 

同理,当γ=2时,由式(22)和式(23)均给出统计误差的Lm范数的数值精度Omγ=2

$ O_m^{\gamma = 2} = \frac{2}{m} $ (45)
 

γ=3时,仅有m≥2时,由式(22)和式(23)给出的统计误差的Lm范数的数值精度Omγ=3相同,且均为

$ O_m^{\gamma = 3} = \frac{2}{m} + 1\;\;\;\;\;\left( {m \ge 2} \right) $ (46)
 

而当m=1时,式(22)给出O1γ=3=2,而式(23)却给出O1γ=3=3,此时式(22)和式(23)给出的统计误差的Lm范数的数值精度并不相同。尽管作为式(22)的简单估算式,式(23)在绝大多数情况下都是成立的,但是式(22)却总是严格成立的。

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

表 4 几何守恒误差的统计精度(γ=1) Table 4 Statistical accuracy of geometry conservation error(γ=1)
NL1L2L3L4L5L
误差精度误差精度误差精度误差精度误差精度误差精度
200 1.24E-03 4.98E-02 1.84E-01 3.53E-01 5.22E-01 2.50E+00
300 7.73E-04 1.16 4.98E-02 0 2.10E-01 -0.34 4.32E-01 -0.50 6.66E-01 -0.60 3.75E+00 -1.00
500 4.39E-04 1.11 4.99E-02 0 2.50E-01 -0.34 5.58E-01 -0.50 9.05E-01 -0.60 6.25E+00 -1.00
800 2.65E-04 1.07 4.99E-02 0 2.92E-01 -0.33 7.07E-01 -0.50 1.20E+00 -0.60 1.00E+01 -1.00
1 200 1.73E-04 1.05 5.00E-02 0 3.35E-01 -0.33 8.66E-01 -0.50 1.53E+00 -0.60 1.50E+01 -1.00
2 000 1.02E-04 1.03 5.00E-02 0 3.97E-01 -0.33 1.12E+00 -0.50 2.08E+00 -0.60 2.50E+01 -1.00
3 000 6.77E-05 1.02 5.00E-02 0 4.54E-01 -0.33 1.37E+00 -0.50 2.65E+00 -0.60 3.75E+01 -1.00
5 000 4.04E-05 1.01 5.00E-02 0 5.39E-01 -0.33 1.77E+00 -0.50 3.61E+00 -0.60 6.25E+01 -1.00
8 000 2.52E-05 1.01 5.00E-02 0 6.30E-01 -0.33 2.24E+00 -0.50 4.78E+00 -0.60 1.00E+02 -1.00
表 5 几何守恒误差的统计精度(γ=2) Table 5 Statistical accuracy of geometry conservation error(γ=2)
NL1L2L3L4L5L
误差精度误差精度误差精度误差精度误差精度误差精度
200 1.90E-03 8.29E-02 3.06E-01 5.88E-01 8.70E-01 4.17E+00
300 8.47E-04 1.99 5.54E-02 1.00 2.34E-01 0.66 4.80E-01 0.50 7.40E-01 0.40 4.17E+00 0
500 3.06E-04 1.99 3.33E-02 1.00 1.66E-01 0.67 3.72E-01 0.50 6.04E-01 0.40 4.17E+00 0
800 1.20E-04 2.00 2.08E-02 1.00 1.22E-01 0.67 2.94E-01 0.50 5.00E-01 0.40 4.17E+00 0
1 200 5.32E-05 2.00 1.39E-02 1.00 9.29E-02 0.67 2.40E-01 0.50 4.25E-01 0.40 4.17E+00 0
2 000 1.92E-05 2.00 8.33E-03 1.00 6.61E-02 0.67 1.86E-01 0.50 3.47E-01 0.40 4.17E+00 0
3 000 8.52E-06 2.00 5.55E-03 1.00 5.05E-02 0.67 1.52E-01 0.50 2.95E-01 0.40 4.17E+00 0
5 000 3.07E-06 2.00 3.33E-03 1.00 3.59E-02 0.67 1.18E-01 0.50 2.40E-01 0.40 4.17E+00 0
8 000 1.20E-06 2.00 2.08E-03 1.00 2.62E-02 0.67 9.32E-02 0.50 1.99E-01 0.40 4.17E+00 0
表 6 几何守恒误差的统计精度(γ=3) Table 6 Statistical accuracy of geometry conservation error(γ=3)
N L1L2L3L4L5L
误差精度误差精度误差精度误差精度误差精度误差精度
200 1.90E-03 8.29E-02 3.06E-01 5.88E-01 8.70E-01 4.17E+00
300 6.02E-04 2.83 3.69E-02 2.00 1.56E-01 1.66 3.20E-01 1.50 4.93E-01 1.40 2.78E+00 1.00
500 1.46E-04 2.77 1.33E-02 2.00 6.66E-02 1.66 1.49E-01 1.50 2.41E-01 1.40 1.67E+00 1.00
800 4.17E-05 2.67 5.20E-03 2.00 3.04E-02 1.67 7.36E-02 1.50 1.25E-01 1.40 1.04E+00 1.00
1 200 1.47E-05 2.57 2.31E-03 2.00 1.55E-02 1.67 4.01E-02 1.50 7.09E-02 1.40 6.94E-01 1.00
2 000 4.18E-06 2.46 8.33E-04 2.00 6.61E-03 1.67 1.86E-02 1.50 3.47E-02 1.40 4.17E-01 1.00
3 000 1.61E-06 2.35 3.70E-04 2.00 3.36E-03 1.67 1.01E-02 1.50 1.97E-02 1.40 2.78E-01 1.00
5 000 5.09E-07 2.26 1.33E-04 2.00 1.44E-03 1.67 4.71E-03 1.50 9.62E-03 1.40 1.67E-01 1.00
8 000 1.83E-07 2.17 5.21E-05 2.00 6.56E-04 1.67 2.33E-03 1.50 4.98E-03 1.40 1.04E-01 1.00
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]的测试方法给定周期性初值:

$ u\left( x \right) = \sin\left[ {2{\rm{ \mathsf{ π} }}x - \frac{{\sin\left( {2{\rm{ \mathsf{ π} }}x} \right)}}{{2{\rm{ \mathsf{ π} }}}}} \right] $ (47)
 

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

$ {O_m} = 3 + \frac{1}{m} $ (48)
 
表 7 统计误差的数值精度 Table 7 Numerical accuracy of statistical error
NL1L2L3L4L5L
误差精度误差精度误差精度误差精度误差精度误差精度
200 8.88E-08 3.01E-07 5.87E-07 8.46E-07 1.06E-06 2.67E-06
300 1.64E-08 4.17 7.66E-08 3.38 1.67E-07 3.10 2.52E-07 2.99 3.24E-07 2.93 8.84E-07 2.73
500 1.90E-09 4.22 1.26E-08 3.53 3.04E-08 3.33 4.80E-08 3.25 6.33E-08 3.20 1.91E-07 3.00
800 2.98E-10 3.94 2.12E-09 3.80 4.88E-09 3.89 7.57E-09 3.93 9.90E-09 3.95 3.07E-08 3.89
1 200 5.77E-11 4.05 5.11E-10 3.51 1.26E-09 3.35 2.01E-09 3.28 2.67E-09 3.23 8.61E-09 3.13
2 000 7.36E-12 4.03 8.54E-11 3.50 2.29E-10 3.33 3.81E-10 3.25 5.20E-10 3.20 1.85E-09 3.01
3 000 1.44E-12 4.02 2.07E-11 3.50 5.95E-11 3.32 1.03E-10 3.23 1.44E-10 3.17 5.78E-10 2.87
5 000 1.85E-13 4.02 3.49E-12 3.49 1.10E-11 3.30 2.01E-11 3.19 2.91E-11 3.13 1.35E-10 2.85
8 000 2.80E-14 4.02 6.85E-13 3.46 2.40E-12 3.25 4.64E-12 3.13 6.94E-12 3.05 3.61E-11 2.80

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

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

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

3 结论

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

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

参考文献
[1] DENG X G, MAEKAWA H. Compact high-order accurate nonlinear schemes[J]. Journal of Computational Physics, 1997, 130: 77-91.
Click to display the text
[2] DENG X G, MAO M L, TU G H, et al. Extending weighted compact nonlinear schemes to complex grids with characteristic-based interface conditions[J]. AIAA Journal, 2010, 48(12): 2840-2851.
Click to display the text
[3] DENG X G, ZHANG H X. Developing high-order weighted compact nonlinear schemes[J]. Journal of Computational Physics, 2000, 165: 22-44.
Click to display the text
[4] JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126: 202-228.
Click to display the text
[5] CASPER J, ATKINS H L. A finite-volume high order ENO scheme for two dimensional hyperbolic systems[J]. Journal of Computational Physics, 1993, 106: 62-76.
Click to display the text
[6] CASPER J, SHU C W, ATKINS H. Comparison of two formulations for high-order accurate essentially non-oscillatory schemes[J]. AIAA Journal, 1994, 32(10): 1970-1977.
Click to display the text
[7] 宫兆新, 鲁传敬, 黄华雄. 虚拟解法分析浸入边界法的精度[J]. 应用数学和力学, 2010, 31(10): 1141-1151.
GONG Z X, LU C J, HUANG H X. Accuracy analysis of the immersed boundary method using the method of manufactured solutions[J]. Applied Mathematics and Mechanics, 2010, 31(10): 1141-1151. (in Chinese)
Cited By in Cnki (7) | Click to display the text
[8] CARPENTER M H, CASPER J. The accuracy of shock capturing in two spatial dimensions: AIAA-1997-2107[R]. Reston, VA: AIAA, 1997.
[9] SVARD M, NORDSTROM J. Review of summation-by-parts schemes for initial-boundary-value problems[J]. Journal of Computational Physics, 2014, 268: 17-38.
Click to display the text
[10] ABARBANEL S, DITKOWSKI A, GUSTAFSSON B. On error bounds of finite difference approximations to partial differential equations-temporal behavior and rate of convergence[J]. Journal of Scientific Computing, 2000, 15(1): 79-116.
Click to display the text
[11] 赵秉新. 一维非定常对流扩散方程的高阶组合紧致迎风格式[J]. 数值计算与计算机应用, 2012, 33(2): 138-148.
ZHAO B X. A high-order combined compact upwind difference scheme for solving 1D unsteady convection-diffusion equation[J]. Journal on Numerical Methods and Computer Applications, 2012, 33(2): 138-148. (in Chinese)
Cited By in Cnki (6) | Click to display the text
[12] 涂国华, 邓小刚, 闵耀兵, 等. CFD空间精度分析方法及4种典型畸形网格中WCNS格式精度测试[J]. 空气动力学学报, 2014, 32(4): 425-432.
TU G H, DENG X G, MIN Y B, et al. Method for evaluating spatial accuracy order of CFD and applications to WCNS scheme on four typically distorted meshes[J]. Acta Aerodynamica Sinica, 2014, 32(4): 425-432. (in Chinese)
Cited By in Cnki (3) | Click to display the text
[13] MAO M L, ZHU H J, DENG X G, et al. Effect of geometric conservation law on improving spatial accuracy for finite difference schemes on two-dimensional nonsmooth grids[J]. Communications in Computational Physics, 2015, 18(3): 673-706.
Click to display the text
[14] HENRICK A K, ASLAM T D, POWERS J M. Mapped weighted essentially non-oscillatory schemes:Achieving optimal order near critical points[J]. Journal of Computational Physics, 2005, 207: 542-567.
Click to display the text
[15] DENG X G, MAO M L, TU G H, et al. Geometric conservation law and applications to high-order finite difference schemes with stationary grids[J]. Journal of Computational Physics, 2011, 230: 1100-1115.
Click to display the text
[16] STEGER J L. Implicit finite difference simulation of flow about arbitrary geometrics with application to airfoils: AIAA-1977-0665[R]. Reston, VA: AIAA, 1977.
[17] STEGER J L. Implicit finite difference simulation of flow about arbitrary two-dimensional geometries[J]. AIAA Journal, 1978, 16(7): 679-686.
Click to display the text
[18] YAN Z G, LIU H Y, MAO M L, et al. New nonlinear weights for improving accuracy and resolution of weighted compact nonlinear scheme[J]. Computers and Fluids, 2016, 127: 226-240.
Click to display the text
[19] YAN Z G, LIU H Y, MA Y K, et al. Further improvement of weighted compact nonlinear scheme using compact nonlinear interpolation[J]. Computers and Fluids, 2017, 156: 135-145.
Click to display the text
http://dx.doi.org/10.7527/S1000-6893.2019.23554
中国航空学会和北京航空航天大学主办。
0

文章信息

闵耀兵, 马燕凯, 李松
MIN Yaobing, MA Yankai, LI Song
CFD中统计误差的数值精度分析
Accuracy analysis of numerical error with statistical forms in CFD
航空学报, 2020, 41(4): 123554.
Acta Aeronautica et Astronautica Sinica, 2020, 41(4): 123554.
http://dx.doi.org/10.7527/S1000-6893.2019.23554

文章历史

收稿日期: 2019-10-09
退修日期: 2019-10-31
录用日期: 2019-11-18
网络出版时间: 2019-11-29 11:40

相关文章

工作空间