2018年启动的国家数值风洞（national numerical windtunnel, NNW）工程是我国工业软件里程碑事件，项目投资巨大，有100多家单位参与，未来影响深远。课题组2014年发表了基于非结构动网格技术的激波装配法后[1]，受到NNW的关注，2019年给予基础研究项目支持，研究目标是探索全场一致高精度计算方法，主要技术指标定为：“对于简单流场的超声速流动问题，全场流动计算精度不低于5阶；对于复杂流场的超声速流动问题，全场流动精度不低于3阶。”高精度格式研究文献中提到的的验证算例（Benchmark）并不多，常用的有双马赫反射、斜锲马赫反射和激波-自由界面相互干扰等几种，主要评价准则是定性考察流场中接触间断是否失稳出现旋涡结构。课题组以前也模拟过这些经典算例，项目申请时设想采用对比装配法和捕捉法的计算结果来落实技术考核指标，但是实施过程中发现装配法模拟这些流动，既无涡旋，也无失稳。虽然难以实现技术指标，但在努力查找原因过程中对计算流体力学（CFD）理解更全面，逐步形成如下观点：采用高精度格式模拟超声速流动得到的计算结果不全是高精度。近期在写文章和为学术期刊审稿中，新观点引起争议。作为投稿人，我写过7页文档回复审稿意见，作为审稿专家，接受过13页的回复。辩论受益匪浅，为使论证逻辑更严谨，现在整理出来供同行参考。根据前期讨论的经验，为关注中心议题，展开讨论之前作如下说明：（1）从1987年提出ENO算起，高精度格式经历30多年研究，文献资料汗牛充栋。因为本文是对所有高精度格式的质疑和批评，全面综述篇幅太多，仅提个别格式容易误解，所以参考文献只引用课题组自己的论文。（2）验证（verification）与确认（validation）是CFD可信度的两个方面。2011年以来课题组承担62个项目，大多来自工程部门的横向课题，让应用部门接受计算数据必须提供可信性论据，做过很多确认工作。本文只讨论验证问题，误差就是数值解和理论解的差异，不关心数值解和风洞试验或飞行试验数据相比精度达到多少个阻力单位。（3）作者从业CFD近40年，经历其迅速发展，享受其丰硕成果，深知每个高精度格式能够发表和应用，必然有其过人之处，本文指出的问题属于瑕不掩瑜，如果能改进就更完美。（4）高精度格式属于CFD专门研究领域，课题组介入并不多，编写求解器时可能出现曲解原始思想的误操作，作者愿意提供文中所有算例的Fortran源程序供检验。1　激波诱导的误差按照CFD理论，采用有限差分格式离散原始流动方程，计算得到的数值解满足修正方程，例如，一维Euler方程∂U∂t+∂F∂x=0 (1)有限差分格式的数值解对应的修正方程为∂U∂t+∂F∂x=∑n=2∞γn∂nU∂xn (2)满足出发方程(1)的初始激波位于两个网格点之间，捕捉法计算过程中逐步变成满足修正方程(2)的数值解。基于守恒型Euler方程的弱解理论认为激波是空间突跃的数学间断，同一空间位置可以有多组流动参数（图1正激波为两组，图 13激波相交有5组）。捕捉法无法描述这种间断特性，计算得到的过渡区本身不存在，用激波前或后的理论值定义数值解误差都不合适，因此，模拟激波时只能定性讨论结果相似性，不可能严谨讨论格式精度及其验证问题。准确模拟数学间断只能采用装配法，关于这个观点的论述参见近期文章[1-11]。10.12050/are20220101.F001图1正激波守恒量误差曲线Fig.1Shock wave conservation error curve10.12050/are20220101.F002图2斜激波示意图Fig.2Diagram of an oblique excitation wave10.12050/are20220101.F003图3一阶迎风涡量分布Fig.3Vorticity distribution of first order upwind scheme10.12050/are20220101.F004图4五阶格式涡量分布Fig.4Vorticity distribution of five order upwind scheme10.12050/are20220101.F005图5波的时空特性示意图Fig.5Diagram of the spatio-temporal properties of waves10.12050/are20220101.F006图6多维效应示意图Fig.6Schematic representation of the multidimensional effect10.12050/are20220101.F007图7一维压力误差曲线Fig.7One-dimensional pressure error curves10.12050/are20220101.F008图8二维压力误差云图Fig.8Two-dimensional pressure error contour10.12050/are20220101.F009图9超声速流动影响域Fig.9Region of influence of supersonic flow10.12050/are20220101.F010图10多维分解效应示意图Fig.10Schematic representation of the multidimensionaldecomposition effect10.12050/are20220101.F011图11捕捉法纹影图[13]Fig.11Capture method of shadow diagrams10.12050/are20220101.F012图12装配法纹影图[3]Fig.12Shock fitting method of shadow diagrams10.12050/are20220101.F013图13激波相交装配法压力云图Fig.13Shock wave intersection pressure contours of shockfitting method本文不讨论数值过渡区局部区域，关注激波之外的误差，因此，称为激波诱导误差。目前常用来验证格式的Sod、Shu-Osher等问题中包含多种Riemann结构，尽管有理论解，但是很难分辨出误差来自激波、接触间断或二者的相互干扰。为了定量确定误差，设计如下最简单的激波问题：流场中仅包含一道直线激波，流动参数转换到相对激波的静止坐标系下，整个计算区域上质量、动量和能量的守恒量理论值为常数，很容易确定数值解的误差，从而进行定量评价。1.1　一维正激波马赫数为Mas的正激波放在以速度uc运动的均匀流中，激波初始位置在xs处，速度采用波前声速无量纲（量纲一）化，初始流动参数及其分布为ρ, u, p=ρ2, u2- uc, p2,  x≤xs 1,- uc,1/1.4,   else          (3)其中流动参数根据激波R-H关系式确定ρ2=2.4Mas2/0.4Mas2+2p2=1.42.8Mas2- 0.4/2.4u2=Mas-Mawp2/ρ2Maw2=1+0.2Mas2/1.4Mas2- 0.2 (4)运动速度uc=Mas对应定常激波，可以得到整个计算域上保持常数的R-H守恒量f1, f2, f3=ρu, ρu2+p, h=Mas, Mas2+            1/1.4, 0.5Mas2+2.5 (5)采用相对误差δfi=Δfi/fi作为评价指标。选择uc=0和Mas≥3参数，处理边界条件没有困难，典型的守恒量误差分布如图 1所示，在过渡区外后产生以特征速度u2运动的波非等熵和以速度u2-a2运动的等熵波。计算结果表明高精度格式比一阶迎风格式的误差极值更大[12-13]，从后面定性分析符合格式特性。1.2　二维斜激波二维情况如图2所示，把以上正激波旋转成与x轴夹角为β的斜激波，为了避免边界条件影响内部流场，可以让β450的斜激波从x轴出发，同时法向马赫数Man≥3满足波后为超声速流动条件。激波前后理论值为常数，是涡量值全场为0的无旋流场，因此，误差评价指标除了采用R-H守恒量，还可以选择涡量ω=0.5(vx- uy)。比较高精度格式和一阶迎风[13]，守恒量误差评价指标的结论同前面一样，但是，涡量评价指标明显不同，如图3和图4所示，前者出现更多细小涡旋结构，不均匀特性更为显著。1.3　定性分析从满足式(1)理论激波的初始间断到满足式(2)数值激波的过渡区，本文认为所有格式均会出现这种诱导误差。现有计算结果表明，高精度格式误差极值比一阶迎风格式更大，多维空间产生的非物理波动现象更复杂，下面分析其原因。超声速流动具有双曲方程特性，CFD教科书建模分析常用波动方程：ut+aux=0。已知初始条件：u(0, x)=f(x)，按照特征线理论，不管初值函数f(x)是间断或连续，经过Δt后理论值为：u(t, x)=f(x-aΔt)，即初始波形位置移动到aΔt。显式格式时间步长Δt存在稳定性要求，一阶格式的CFL条件：aΔt/Δx≤1，在Δt内扰动传播范围不会超过一个网格：aΔt≤Δx，对应的物理机理是来自网格边界j-1/2的扰动信息不能影响到边界j+1/2，如图5所示。高阶格式满足CFL条件的时间步长更小，理论上波传播距离更短。目前构造高精度格式常用多点模板，例如，采用7点构造的五阶格式计算初始位于i和i+1之间的激波，无论Δt如何取值，时间推进一步后i+3点的流场参数发生变化，波未致解已变，传播的信息只能是误差。算法改变了波的传播速度，涉及点越多误差越大，出现高精度格式误差更大的计算结果也就不足为奇了。二维情况下，直线激波和网格不平行，初始激波形状就存在弯折，时间推进中扰动沿x方向和y方向传播速度不同，如图6所示，相同的扰动从A和B两点出发向C点传播，如果传播速度有差异到达C点有先后，就会引出各向异性现象，这是导致涡量非0的原因。高精度格式误差传播速度机理复杂，和维数组合产生结构丰富的非物理现象也合乎情理。2　接触间断诱导的误差考察以上一维激波算例数值过渡区内两个网格点流动参数的Riemann结构，发现随时间变化规律复杂[13]。由于激波前后密度、速度和压力均有变化，为剥离激波影响和便于建模分析，设计密度变化的接触间断验证算例。2.1　一维接触间断接触间断放在速度为u的流场中，采用波前声速无量纲化，初始位置在xd处的流动参数分布为ρ, u, p=ρ2, u, 1/1.4, x≤xd1, u, 1/1.4, else  (6)对应的间断前后马赫数为Ma1=u，Ma2=uρ2在计算区域内压力和速度为常数，均可作为误差评价指标。数值结果表明，在一侧或两侧亚声速情况下，通量分裂格式的影响比差分格式更为明显，通量差分分裂格式(flux difference splitting, FDS)能够维持压力和速度为常数，矢通量分裂 (flux vector splitting, FVS)出现了两道以u+a1和u-a2速度运动、压力和速度有变化的非物理波，典型结果如图7所示。2.2　二维接触间断为了考察定常流场中误差的空间传播特性，将上述算例的接触间断旋转90o后放置在超声速流场，提出如下二维验证模型ρ, u, v, p=ρ2, u, 0, 1/1.4,     y≤0    1, u, 0, 1/1.4, y0 Ma1=u1, Ma2=uρ21 (7)全场压力p=1/1.4为常数，y方向速度v=0，上下两层x方向速度相等。同前面的二维激波算例一样，误差评价指标可以选择压力、速度或涡量。为避免边界条件的影响，调整密度满足全场超声速流动的条件。实际计算中可以组合成多结构讨论边界层问题，图8为三层结构的压力误差云图，出现了6道误差波，相互干扰在下游形成非常复杂的流场结构。2.3　定性分析目前CFD领域教科书根据数值结果认为FVS比FDS的数值耗散更大。本文发现对于u≠0运动流场，从满足式(1)的初始间断开始计算中形成满足式(2)的接触间断过渡区，在形成计算初期表现为FVS密度过渡层厚度确实比FDS大，但是随着FVS的两道非物理波明显远离，再比较密度过渡层内的密度分布曲线，相同差分格式实际上相差不大，由此看来上述观点值得进一步探讨。由于篇幅所限，图7和图8计算模型的更多细节介绍，以及对接触间断诱导误差特性的讨论另文详述。根据CFD教科书截断误差定义，均匀流场满足任何格式修正方程，换而言之，如果从参数为常数的初场开始，计算过程应该维持不变，或者稳定以后流场维持常数。采用高精度格式计算流动参数空间线性分布的流场，理论上收敛解误差为0。保持其他参数为常数，采用FVS计算u=0的静止流场和u=cy的线性流场，无论一阶迎风还是五阶格式，都会出现误差。由此提出疑问：采用FVS的计算结果是否有一阶精度？图9在空气动力学教科书中常用来介绍超声速流体特性的示意图，扰动只能在特征线围成空间有限的影响域内传播，其他区域即使近在咫尺也无法感知。据此考察图 8，影响范围符合理论，但是误差分布在影响内变化明显，离开流线和特征线的大部分区域内流场恢复到初场均匀流参数。图 8的误差分布体现了定常流动特点，即扰动和影响域内大部分点之间没有关系。沿着特征线在上游形成依赖域，当地点不能感受到其依赖域之外的信息。如图10所示二维空间超声速流场，非定常流动只有上游依赖域内a- i对计算点0产生影响，如果只关心定常解，原则上应该消除其中有些点（如b和d点），否则传递过来的信息只能增加误差。目前很多时候采用维数分裂方法（dimension by dimension）把基于一维空间构造的高精度格式推广到多维空间，在y方向的差分格式把依赖域外1-4点牵扯进来，这些点给0点数值解带来的影响只能是误差。同样，y方向涉及的点越多引入的误差越大。讨论中同行提出“1-4点看似无关，但是属于a-i点的影响域内，因此，间接传递到0点的信息也是合理的”，这个看法定性合理，但是缺乏直接证据。3　接触间断诱导的误差验证高精度格式的经典问题，如斜锲马赫反射、双马赫反射、激波和自由界面相互干扰等问题，文献常把流场中是否出现接触间断失稳形成的涡街作为高精度格式标志。根据前期争论的经验，本文区分成两个问题：（1）接触间断失稳和涡街是不是非物理波动；（2）用这种现象评价高精度格式是否科学。3.1　失稳和涡街2.1节和2.2节中接触间断算例中出现的波动是非物理波，这个可以通过误差来论证，同行很少有争议。但是，经典验证算例中有些失稳和涡街不是初始间断产生的，而是计算过程中形成的，因此，如何证明是非物理波动？下面先给出数值结果，然后进行理论分析。从无黏流动Euler方程出发，模拟激波和自由界面相互干扰问题。图 11是2008年课题组采用捕捉法的计算结果，和试验图片非常相似，流场中清晰可见界面失稳和旋涡在超声速流动引起分布广泛的细小结构，详细比较不同时刻界面位值也符合很好，在此基础上称赞所采用的算法具有高分辨的优点[14]。近年课题组研究装配法，仅装配气泡界面，其他间断采用捕捉法，计算结果如图12所示，整个计算过程界面未出现失稳，流场中也没有旋涡及其诱导波系结构[3，14]。从精度角度看，装配法直接采用接触间断的数学表达式，理论上是准确的，根据图 12来评价图11，这些复杂流场结构是错误的计算结果。上述算例流动是非定常的，只能定性分析，设计如图 13所示的异侧斜激波相交定常流动模型进行定量考察。为消除边界影响，调整来流马赫数Ma∞、上下激波角β1和β2参数，整个流场的理论值均为超声速。尽管流场存在各种间断，但是其他区域流动参数为常数，采用装配法即使一阶迎风格式也能得到符合理论解的计算结果，图13是典型流场云图[16]。选取激波角相同条件下，采用高精度格式进行数值模拟，在流场中心对称线处设置监测点记录压力，位置如图14所示，计算中测点波动曲线如图15所示，出现明显的周期性振荡。图 16是交点下游涡量云图，从图16可以看出，类似于文献中经典算例的接触间断失稳和涡街现象。值得指出的是，在β1=β2条件下理论上没有接触间断，但是数值模拟结果表明，下游测点波动振幅明显变小，沿着对称线依然存在强度较弱的涡街。10.12050/are20220101.F014图14测点位置示意图Fig.14Diagram of the location of the monitoring points10.12050/are20220101.F015图15监测点压力时程曲线Fig.15Pressure variation curves over time at monitoring points10.12050/are20220101.F016图165阶WENO计算涡量云图Fig.16Contours of vorticity for five order WENO scheme激波正规相交形成定常流动，这些非定常波动属于算法引起的不稳定现象，这个算例可以作为“接触间断失稳和涡街”非物理波动的直接论据，此外还可以根据常微分方程定性理论进行理论分析。对于无黏定常流动，考虑重力影响，参数沿流线满足伯努利方程12u2+pρ+hg=const (8)式中：右侧常数与密度无关。接触间断满足流线定义，可以看成两侧介质共用的流线，各自满足伯努利方程。假如接触间断发展成涡街，二维情况下涡旋中心为螺点或中心点，两个涡旋之间鞍点连接，在以鞍点速度uc运动的相对坐标系中，接触间断两侧流线在鞍点相交，如图17所示。图中鞍点处压力和高度相等，伯努利方程变为12uc2+pρ1+gh=C1 (9)12uc2+pρ2+gh=C2 (10)10.12050/are20220101.F017图17鞍点结构示意图Fig.17Schematic diagram of the saddle point消除压力参数后，得到鞍点速度uc2=2ρ1C1- ρ2C2ρ1- ρ2-2gh (11)进一步推导出存在uc≠0鞍点条件ρ1C1- ρ2C2ghρ1- ρ2 (12)如果接触间断两侧的密度、速度和高度组合满足以上条件就可以出现鞍点，即使无黏流动方程也可能出现“波浪滚滚”的流场，这是模拟Rayleigh-Taylor不稳定现象的理论基础。在不考虑重力影响的条件下，式（12）变为ρ1C1ρ2C2 (13)无黏流动沿流线总焓相等，对于比热[容]比为常数的量热完全气体，不难证明式（13）结论同样成立。激波相交算例中流场满足C1=C2，存在鞍点条件变为ρ1ρ2 (14)接触间断两侧密度根据上下激波角β1和β2参数确定，如果β1和β2满足以上条件出现了涡街，那么将参数取β2和β1时就没有了。斜锲马赫反射、双马赫反射等非定常算例，接触间断也是均匀来流产生的，激波后流场满足C1=C2，按照鞍点条件同样推断出把接触间断参数互换后涡街变化的结论。这种不满足对称性的现象显然是不可能的，因此，流场中没有鞍点，也不会形成涡街。综上所述，在不考虑重力影响Euler方程的流场中，不管是初场固有的还是后来产生的接触间断，理论上难以出现失稳和涡街，目前高精度格式研究文献给出的验证算例中模拟出来的是非物理波动。3.2　评价标准满足式(1)的理论接触间断厚度为0，差分格式计算过程中数值黏性使其变成若干网格、满足方程(2)的数值剪切层。数值剪切层越薄精度越高，高精度格式模拟出现失稳和涡街后厚度明显增加，如果采用厚度作为评价指标误差更大，显然不符合逻辑。基于以上数值结果和定性分析，认为“出现失稳和涡旋是捕捉法过渡层内数值黏性引起的，按照错误原因不能产生正确结果的逻辑，用来证明高精度无异于说错误越离谱精度越高”。对此，同行提出如下反驳意见：(1)在高精度领域中，诸如双马赫反射，R-T不稳定性等无黏算例中的涡结构是由于数值黏性造成的，数值黏性越小，则修正方程的雷诺数越大，也就能搓出来更多的涡结构。因此求解这些经典算例的意义在于，可以通过流动细节的多少来推断修正方程的雷诺数，通过雷诺数来推断数值黏性，数值黏性越小则分辨率越高。(2)该问题的奇异性在于数值黏性越趋近于0，产生的涡越多，而当数值黏性等于0的时候，涡就“搓”不出来了。通过误差变化规律来揭示机理是科学研究常用的技术途径，因此，同行的观点有道理，采用“这些被广泛应用的经典算例”来论证高精度格式数值黏性大小符合逻辑。因为接触间断失稳和涡街现象定性是非物理波动，然后又用来定性验证数值黏性，论证过程很难严密，容易陷入逻辑困境。高精度格式研究中空间一维条件下，经常通过比较Fourier级数得到的波数和耗散特性曲线来讨论数值耗散特性。多维空间情况下理论分析困难，如果通过比较高精度格式计算得到的漩涡大小、数目合理和网格尺度关系，建立类似的定量评判科学标准，可以直接作为论证“数值黏性更小对大尺度流动细节分辨能力更强”的依据。4　膨胀波流场的精度讨论中同行认可在间断附近能够捕捉激波，保持鲁棒性，高精度格式也要降阶处理，但是不影响在光滑区域的使用精度。去掉超声速流场中激波和接触间断后只剩下膨胀波。一维空间的膨胀波难以稳定，非定常计算误差中辨识格式时空影响有难度，空气动力学教科书常见的Prandt-Meyer流动模型是最简单的二维膨胀波，有解析解的定常流动，可以计算误差，适合作为验证算例。由于Prandt-Meyer膨胀波是超声速流动适应弯曲物面形成的流动现象，计算涉及到曲线网格。在贴体坐标系下非均匀网格上应用差分格式，会出现几何诱导误差，常用自由流问题检验误差特性，即考察格式在不断迭代过程中能否保持初始均匀流场参数不变化或维持在计算机机器0。数值模拟Prandt-Meyer膨胀波之前，先介绍坐标变换诱导误差的消除算法。4.1　柱坐标模型从高精度格式应用早期就发现几何诱导误差比低阶格式更大，文献常用网格正交性和光滑性来解释，普遍观点是“高阶精度格式要求网格也是高阶光滑”。这个观点不完全准确，从笛卡尔坐标变换到柱坐标解析表达式r=x2+y2θ=arctan(x/y),x=rsinθy=rcosθ (15)为了消除边界处理引入的误差，建议采用超声速流动进出口边界条件，而且仅选择上游的半圆区域，即θ∈[π2,3π2]或θ∈[- π2,π2]。为避免网格奇性，取半径r0。坐标变换函数各阶导数存在且连续，生成具有完全正交充分光滑的均匀网格，计算Ma∞=3自由流问题。通量分裂格式采用VanLeer、AUSM、HLLC和Roe，坐标变换系数解析式和二阶差商数值解，比较了各种组合，验证了文献结论：“高阶格式误差比一阶迎风格式更大”，还发现高精度格式出现沿周向θ不均匀现象[17]。目前常用误差随网格变化斜率评价格式精度（精度测试），从初始径向21点、周向41的基础上采用不断均分的措施进行网格加密，流场内最大误差作为评价指标，精度测试曲线如图18所示，根据一阶精度（黑色虚线）和二阶精度（红色虚线）参考斜率评估误差下降速率，可以看出，五阶算例的收敛精度远小于其格式理论精度。10.12050/are20220101.F018图18全场最大压力误差随网格加密的变化Fig.18Maximum pressure errors with grid refinement过去消除坐标变换几何诱导误差的主要措施是针对每个差分格式建立相应的几何守恒律（geometric conservation law，GCL）算法。近年来本课题组提出离散等价方程及其离散准则（discrete equivalence equation and its discrete-rule，DEER），理论分析对坐标变换函数和差分格式没有限制，一阶迎风或五阶格式组合以上各种通量分裂格式模拟自由流问题，误差从始到终保持机器0，验证普适性。DEER具体理论推导过程、验证算例以及和GCL算法差异说明参阅参考文献[17]~文献[19]。4.2　Prandt-Meyer流动设计模型的计算区域和边界处理如下：根据流线确定固壁外形，采用延伸到固体内部的虚拟点处理边界，根据有限差分法理论，解析解代入差分格式的修正方程和流动方程相差为高阶小量，因此，在计算中这些虚拟点上固定为理论值不影响格式精度。在膨胀波上游设置多排均匀直角网格处理超声速入口边界，根据格式需要固定若干排为来流参数。采用简单外推处理出口边界和上边界对内部流场的影响范围有限，在考察精度时选择明显未被干扰的区域。图19为网格及测点示意图。10.12050/are20220101.F019图19网格及测点示意图Fig.19Schematic diagram of the grid and monitoring points入口马赫数Ma∞=1.067，对应的气流折转角度θ=20o，根据法向网格与壁面垂直，各向网格尺度Δx=Δy=0.09098。首先计算自由流问题，二阶MUSCL格式和五阶格式密度误差分布如图20所示，采用DEER算法可以消除几何诱导误差。在此基础上均分网格加密，图21是4套网格的局部示意图。10.12050/are20220101.F020图20几何诱导密度误差云图Fig.20Geometrically induced density error contour10.12050/are20220101.F021图21网格加密示意图Fig.21Schematic diagram of grid refinement定义压力误差绝对值参εp=p-pexct，图 22是该参数的全场云图（最粗网格），可以看出，膨胀波到达上边界后误差剧烈增加，这是边界点算法所致，因此，只考察未受到上边界影响区域的流场。10.12050/are20220101.F022图22压力误差云图Fig.22Pressure error contours在x≤1的有效区域内，最大误差在膨胀波起始位置，这是流场参数的导数不连续所致。误差等值线基本上沿特征线，符合超声速流动理论特征。选择局部伞形区域固定位置上的误差进行精度测试，如图23所示，根据一阶精度（黑色虚线）和二阶精度（红色虚线）参考斜率评估误差下降速率，这两个格式在误差带内只有一阶精度。测点2（蓝点）的精度测试曲线如图23(b)所示，也是同样结论。10.12050/are20220101.F023图23测点压力误差随网格加密的变化Fig.23Pressure errors with grid refinement of monitoring points从误差云图看，误差极值集中在膨胀波起始位置的局部伞形区域内，下游区域又恢复到高精度。上游一阶精度误差对下游影响有限，似乎违背超声速的扰动传播特性，对这种现象未见其他文献报道，机理认识有限。计算结果表明，高精度格式不但在激波和接触间断处发生降阶，对参数连续的膨胀波也不一定全是高精度。5　精度测试的合理性验证和确认是可行度的两个方面，在具体实施层面，确认工作较为简单，把给定的试验结果（数据或图像）作为标准，采用数理统计方法分析数值解误差和网格、差分格式、通量分裂格式、湍流模型和边界条件等因素的关联特性，只要数据量够大总能给出一些定量的评价指标和定性的影响结论。验证考察计算机代码正确求解方程的近似程度，以方程精确解作为评价标准，由于能够给出理论解的流动问题很少，因此，开展工作难度较大。通过确认过程给出CFD软件计算结果和试验的偏差，从而有效使用CFD数据对于工程应用意义重大，但是用于算法验证则存在风险。例如，前面已经论证接触间断失稳和涡街为非物理波动，数值流场图像和超声速剪切层流动非常相似，如果从确认角度看，FDS通量格式没有模拟出这种波动现象，计算结果比FVS较差，但是从验证角度看，对于无黏流动方程，维持上下两层参数不变化的FDS结果是精确的理论值。在没有解析解的情况下，常采用误差随着网格加密过程变化曲线的斜率来测评格式精度，简称精度测试。图18的计算中，除了格式就是坐标变换函数影响，相对比较简单，但是，在图23的计算中，误差还受流场梯度影响，无法分辨多种因素对误差的贡献，很难定性分析，因此，精度测试结果经常出现难以用CFD理论解释的非整数阶数。精度测试手段除了缺乏严谨的理论基础，也存在定性结论错误的风险。2008年前后，中国空气动力学学会计算流体力学专委会组织过“CFD验证与确认”专题研讨，会后本文作者有感而发，从Burgers方程出发建立了数值解精度预测模型，即给定网格以后就可以确定误差，2011年指导本科生毕业设计编程数值验证[20]。由于不易获取原始文献，下面对其中与本文相关内容进行介绍。在区间[0，1]两端给定边界条件，Burgers方程及其理论解dudx=1μd2udx2              μ0 u(0)=0 u(1)=1  (16)u(x)=exp(μx)-1exp(μ)-1 (17)在[0，1]区间均匀分布间距Δx的网格，计算中两侧边界直接固定为理论值。根据CFD理论，数值解u*满足的修正方程du*dx=1μd2u*dx2+∑k=2∞γkdku*dxk(Δx)k-1 (18)将真值u代入修正方程得到差分格式的局部截断误差R(x, Δx)=∑k=2∞γkμk(Δx)k-1exp(μx)exp(μ)-1=εμd2udx2 (19)式中：ε=∑k=2∞γk(μΔx)k-1，称为数值黏性系数。按照以上公式可以计算一阶、二迎风格式和二阶中心格式的数值黏性系数，如图 24所示，其中横坐标ReΔx=μΔx称为网格雷诺数。可以看出，尽管ε随ReΔx单调变化，但是，迎风和中心格式的符号相反。由此推测，用迎风和中心组成的混合格式，ε随ReΔx变化规律不一定保持单调，ε本身出现正负变号也是完全可能的。根据数值解误差和数值黏性的关系，如果ε不再单调或发生正负变号意味网格尺度减小引起误差增大。10.12050/are20220101.F024图 24数值黏性随网格变化Fig.24Numerical viscosity with grid refinement由于TVD、ENO等格式的限制器、WENO和WCNS的模板函数与流场参数相关，很难推导出ε表达式。但是，对于Jameson格式，只要固定二阶和四阶人工黏性项的系数，完全可以推导出ε随ReΔx变化规律。在理论指导下，调整人工黏性项系数，出现了ε正负变号情况，选择数值解最大误差Δumax作为精度评价参数，计算结果如图 25所示，验证了网格加密误差增大的理论预测。10.12050/are20220101.F025图25最大误差随网格变化Fig.25Maximum error with grid refinement证明定理需要进行各种逻辑推理和仔细推敲，而否定它仅需一个反例。图25给出的结果表明，采用精度测试手段评价格式精度不完全可靠。课题组目前推广DEER算法消除黏性项几何诱导误差取得成功，正在从N-S方程出发寻找证据验证以上结论。6　结束语数值结果和理论分析表明，没有计算超声速流动的高精度格式。CFD领域还有许多值得关注的问题[21]，需要开展大量的基础研究工作。