双曲守恒律方程有着广泛的应用，在内爆动力学，惯性约束聚变，计算天文学，磁流体力学，计算爆炸力学，计算生物学等领域都和它密切相关。它有一个明显的特征，无论初值条件多么的光滑，随着时间的演化，最终都可能会演化为带有强间断的解。爆炸与冲击问题是强非线性的双曲方程，这些方程的解通常带有奇异性，也就是在解的附近存在较大的变化，包括激波，稀疏波，接触间断等。为了克服上述难题，一些高精度的数值格式被提出，如MUSCL、WENO^[1-2]等。但采用高阶重构时，解在间断附近会产生虚假震荡。因此构造有效的高精度数值算法，不仅要求能捕捉到间断区域，同时又可以消除虚假震荡。针对虚假震荡的非物理现象，可以考虑将原始守恒方程映射到特征空间上^[3]，在特征空间上求得数值通量之后^[4]，再将通量映射回原始空间。这种通过对方程组多变量的耦合进行解耦的过程，可以巧妙地避开虚假震荡。

对于常规的有限体积方法，计算网格是均匀网格，若提高计算精度，那么相应的策略是增加网格数目。但是对于物理量梯度变化较小的区域，粗糙的网格便可以得到满意的结果，对光滑区域求解增加的计算量是没有必要的，因此这无疑在一定程度上降低了计算效率。基于此，根据解的性质，对网格进行合理的节点分布，对于高效、精确地计算双曲守恒方程有着极其重要的作用。伪弧长算法是解决这种问题的有效方法之一。在伪弧长算法中，网格节点会随着时间的变化而变化，在保证网格节点数目不变的情况下，在物理量梯度较大的区域自动加密，而在解较为平滑的区域自动稀疏。Tang等^[5]通过求解一个基于变分原理的网格方程得到了一种移动网格方法，2006年Tang^[6]将移动网格算法拓展到二维情况，近年来，自适应网格算法得到了较大的发展^[7-9]。

本文将高精度的WENO格式同伪弧长算法相结合，针对网格移动之后造成的网格非均匀和非正交使得通量的计算和网格物理量重构的过程较为复杂的现象，提出采用坐标变换，将坐标映射到均匀正交的计算坐标系下进行计算。结果表明这样可以提高数值精度。通过一维激波管问题结果的比较，可以看出伪弧长算法相较于传统的固定网格算法可以更好地捕捉到强间端，而高阶伪弧长算法对间断捕捉的分辨率最高。结果表明高阶伪弧长算法相较于有限体积以及低阶格式的伪弧长算法，它更适合处理处理爆炸与冲击强间断问题。此外，可以改变弧长参数使得网格移动更加明显，进而提高对间断捕捉的分辨率。

1.   伪弧长算法

本文采用的伪弧长算法包含3个部分：第1部分给出物理空间的控制方程在时间和空间上的离散格式；第2部分是坐标变换，将不均匀的物理空间变换到均匀的计算空间中进行求解；第3部分则是伪弧长算法。

1.1   离散格式

考虑如下的双曲守恒系统：

式中：${\boldsymbol{w}}$为质量、动量、能量组成的守恒变量向量，${\boldsymbol{F}}({\boldsymbol{w}})$、${\boldsymbol{G}}({\boldsymbol{w}})$、${\boldsymbol{S}}({\boldsymbol{w}})$为${\boldsymbol{w}}$的函数。

对式(1)在网格${K_{i,j}}$单元上进行积分，得到：

式中：${\text{d}}\sigma$为面积微元，${K_{i,j}}$为网格区域。

对上式采用Green公式进行变换，则有：

式中：${\overline {\boldsymbol{w}} _{i,j}}$为网格的物理量平均值，$\left| {{K_{i,j}}} \right|$为${K_{i,j}}$的面积，$\partial {K_{i,j}}$为网格的边界，ds为网格边界长度微元（取逆时针方向为正），$\delta ({\boldsymbol{w}}){\text{ = (}}{\boldsymbol{F}}({\boldsymbol{w}}){\text{,}}{\boldsymbol{G}}({\boldsymbol{w}}){\text{)}}$，${{\boldsymbol{n}}_{i,j}}$为边界$\partial {K_{i,j}}$的单位外法向量。

数值通量采用局部Lax-Friedrichs格式，详细定义如下：

式中：u和v为守恒变量的内外重构值，${\boldsymbol{n}}$为网格每一条边的为单位外法向量，$\ \mu {\text{ = max(}}{{\boldsymbol{\delta }}'_{\boldsymbol{u}}} \cdot {\boldsymbol{n}},{{\boldsymbol{\delta }}'_{\boldsymbol{v}}} \cdot {\boldsymbol{n}}{\text{)}}$。

式(3)可以写成半离散格式：

式中： ${\partial ^k}{K_{i,j}}$为网格${K_{i,j}}$的第$k$条边，$k = 1,2,3,4$；$\left| {{\partial ^k}{K_{i,j}}} \right|$为第k条边的长度；${\boldsymbol{w}}_{i,j}^{{\text{int}}(k)},{\boldsymbol{w}}_{i,j}^{{\text{ext}}(k)}$为边界${\partial ^k}{K_{i,j}}$的内外物理量，${\boldsymbol{w}}_{i,j}^{{\text{int}}(k)},{\boldsymbol{w}}_{i,j}^{{\text{ext}}(k)}$的二阶MUSCL形式详见文献[10]。下面给出WENO格式的内外物理量具体形式。

以下仅对x方向进行离散，y方向仅需要将变量i替换为j即可。

式中：${\boldsymbol{w}}_l^ + \left( {{x_{i + 1/2}}} \right)$和${\boldsymbol{w}}_l^ - \left( {{x_{i - 1/2}}} \right)$为网格右边界和左边界守恒变量的的线性插值多项式，${{\boldsymbol{w}}_{{\rm{R}},i}}$和${{\boldsymbol{w}}_{{\rm{L}},i}}$和网格${K_{i,j}}$的左右内边界守恒变量的非线性加权值，$\omega _l^ \pm$为加权因子，$d_l^ \pm$为理想光滑区域的最优化系数，$\ \beta _l^ \pm$为光滑因子函数，$ϵ$和$\lambda$为常数。

对于时间的离散，忽略掉网格的索引$i,j$，则采用如下的四阶TVD-RK格式：

式中：$\Delta t$为时间步长，${\overline {\boldsymbol{w}} ^{\textit z}}$为上步物理量，${\overline {\boldsymbol{w}} ^{{\textit z} + 1}}$为更新的物理量，${\overline {\boldsymbol{w}} ^{(k)}}$为求解过程的中间物理量值，${L_j}\left( {{{\bar u}^{(k)}}} \right)$为微分算子，$k = 0,1,2,3$。

式(8)和式(9)的重构中，重构的物理量可以选择原始变量，守恒变量以及特征变量。由于欧拉格式的非线性性质，采用高阶格式容易引起非物理震荡，因此可以考虑采用将变量映射到特征空间中，在特征空间中进行重构，重构完成之后再映射回原始空间，具体计算如下：

先计算$\dfrac{{\partial {\boldsymbol{F}}({\boldsymbol{w}})}}{{\partial {\boldsymbol{w}}}}$的左特征矢量矩阵${\boldsymbol{L}}({\boldsymbol{w}})$，元素可以采用Roe平均，也可以简单的取平均值，然后进行变量映射，${\boldsymbol{w'}} = {\boldsymbol{L}}({\boldsymbol{w}}){\boldsymbol{w}}$，将映射完的特征变量用WENO格式重构出网格边界处的值，最后用${\boldsymbol{w}} = {\boldsymbol{R}}({\boldsymbol{w}}){\boldsymbol{w'}}$ 返回到守恒变量或者原始变量。

1.2   坐标变换

在伪弧长算法中，移动之后的网格是非均匀网格，而传统的数值格式均是基于均匀网格给出的。通过限制网格的移动距离，可以近似用传统的格式进行插值运算，然而仅仅在一阶或者二阶情况下，这种插值引起的误差不大。若直接运用到高阶格式，将会引起较大的误差。因此，对于高阶伪弧长算法的通量计算，必须考虑网格尺度的影响。为了解决这种问题，通常有两种策略，第一种策略是将网格尺度引入WENO格式的非线性权重，然后直接进行构造高阶格式^[11-13]。这种构造清晰明了，但是格式较为复杂，而且不易推广到二维或者三维情况。另外一种策略是转换坐标系，将当前坐标系变换到曲线坐标系下，因为曲线坐标系中计算网格是均匀正交的，因此可以直接用传统的格式进行通量计算。详细变换如下。

引入曲线坐标

则式(1)可以由原始坐标转换到曲线坐标下进行计算，即

式中：ρ为密度；p为压力；u和v为速度；$U = {\xi _x}u + {\xi _y}v + {\xi _t}{\kern 1pt} ,{\text{ }}V = {\eta _x}u + {\eta _y}v + {\eta _t}{\kern 1pt} {\kern 1pt}$；ψ=w/J，J为Jacobian行列式，表达式为

通过式(12)在计算坐标下求解物理量，求解完成之后进行逆变换${\boldsymbol{w}} = J{\boldsymbol{\psi }}$，则可以得到原始物理量。式(12)的计算核心是计算Jacobian矩阵$J$和度量系数${x_\xi }\text{、}{x_\eta }\text{、}{y_\xi }\text{、}{y_\eta }\text{、}{x_\tau }\text{、}{y_\tau }$。而$J$的推进可以采用几何守恒条件进行计算，具体如下

因此只要给出$J$的初值，$J$的推进可以利用上式进行计算。其中：${x_\tau },{y_\tau }$为网格的运动速度。对于速度的表达式，详见下面的伪弧长算法。下面给出${x_\xi }$的计算过程，其他度量系数的计算则类似。

式中：为了提高数值精度，采用了六阶中心差分进行计算${x_{\xi + 1/2}}$，但采用其他计算格式如WENO也是可以的。因为计算空间是均匀分布的，所以${\varDelta _\xi } = 1$。

1.3   网格移动算法

为了保证网格向梯度大的地方移动，文献[5]基于变分方法，提出仅需求泛函$E(x,y)$的最小值即可，$E(x,y)$的表达式如下

式中：$\xi = \xi (x,y),\eta = \eta (x,y)$表示弧长空间坐标，${{\boldsymbol{G}}_1}\text{、}{{\boldsymbol{G}}_2}$是给定的正定对称矩阵，为了方便计算令${\boldsymbol{G}} = w{\boldsymbol{I}},$$\nabla = {({\partial _\xi },{\partial _\eta })^{\text{T}}},$${\partial _\xi } = \partial /\partial \xi ,{\partial _\eta } = \partial /\partial \eta$，根据变分原理知(15)的E-L方程为

伪弧长控制函数取为

式中：${a_1}\text{、}{a_2}$为可调参数，它的大小决定了网格的移动程度。

为了保证网格的光滑性，可以用低通滤波器进行光滑处理，一维情况表达式为

二维情况采用下述公式

式中：${\kappa _i}\text{、}{\kappa _{i,j}}$为中间控制函数，作为中间保存步备用。

具体光滑处理次数应该视情况而定，通常情况下处理2~3次即可，若非线性情况较严重，可以处理10次左右。

因此方程(16)改写为

通过对式(20)的计算，可以求得网格的分布。式(20)的计算可以采用Jacobian迭代进行计算，计算步骤如下

式中：${\boldsymbol{x}}_{_{i,j}}^{({\textit z} + 1)}$为下一时刻网格的坐标$(x,y)_{_{i,j}}^{({\textit z} + 1)}$，其中$0 \leqslant i \leqslant N,0 \leqslant j \leqslant M$，并且对于边界（i=0, j=0, i=N或j=M），网格保持不动。

因此根据式(21)可以得到式(13)的网格移动速度为

式中：$\Delta t$为时间步长。

网格移动完成之后，接下来的工作就是需要得到物理量在新网格上的值，这一重要步骤又称物理量重映。常用的物理量重映方法分为2类，一是基于通量的物理量重映方法，二是基于网格相交的精确积分守恒重映方法。考虑到在相邻的时间间隔要求，网格移动距离较小，而且网格始终处于保凸性，因此本文采用的映射方法为第一种。

如图1所示，${D_k}$为新旧四个节点组成的通量四边形，根据文献[5]的工作，可以得到新旧物理量之间的关系为：

图  1  新旧网格转化示意图

Figure  1.  Diagram of old grid transforming to new grid

下载: 全尺寸图片 幻灯片

式中：$\left| {K_{i,j}^{({\textit z} + 1)}} \right|$代表新网格的面积 ，${\boldsymbol{w}}_{i,j}^{({\textit z} + 1)}$代表新网格的物理量，$\displaystyle \sum\limits_{k = 1}^4 {{F_k}({\boldsymbol{w}}_{i,j}^{{{\rm{int}}} (k)},{\boldsymbol{w}}_{i,j}^{{\text{ext(}}k)})}$代表第k个四边形的通量，${\boldsymbol{w}}_{i,j}^{{{\rm{int}}} (k)},{\boldsymbol{w}}_{i,j}^{{\text{ext}}(k)}$的计算见式(16)。通量的计算过程如下：

式中：${D_k}$代表通量四边形的四个顶点按逆时针进行排序计算出的面积。

综上，伪弧长算法的详细步骤如下。

第1步，初始化，给定初始时刻网格节点的坐标$x_{i,j}^{(0)},y_{i,j}^{(0)}$，初始物理量${\boldsymbol{w}}_{i,j}^{(0)}$。

第2步，求解网格控制函数式(17)，并用式(18)和式(19)进行光滑打磨控制函数。

第3步，伪弧长网格移动，根据式(21)进行迭代计算，迭代的限制条件为$|{\boldsymbol{x}}_{i,j}^{(\textit z)} - {\boldsymbol{x}}_{i,j}^{({\textit z} + 1)}| < \varepsilon$。其中，$|{\boldsymbol{x}}_{i,j}^{(\textit z)} - {\boldsymbol{x}}_{i,j}^{({\textit z} + 1)}|$代表相邻时间步网格的距离差，$\varepsilon$为一个自定义的极小常数。最终求得网格的坐标${\boldsymbol{x}}_{i,j}^{({\textit z} + 1)}$。

第4步，物理量映射，根据式(23)更新新网格下的物理量。

第5步，求解物理量控制方程：(a)如果采用低阶伪弧长算法，则直接在物理空间进行求解。通过等式(10)更新${t_{n + 1}}$时刻的物理量${\boldsymbol{w}}_{_{i,j}}^{(n + 1)}$；(b)如果采用高阶伪弧长算法，则利用式(12)、式(13)和式(14)将物理空间的物理量转换到均匀计算坐标系下进行求解。通过式(10)更新${t_{n + 1}}$时刻的物理量${\boldsymbol{w}}_{_{i,j}}^{(n + 1)}$，求解完成之后根据逆变换再把物理量映射回物理空间。

第6步，如果求解达到终点时间，则程序终止，否则转到第2步。

2.   一维数值算例

下面基于伪弧长数值算法进行一些算例测试，由于上面的理论是基于二维情况推导的，对于一维情况直接进行降维处理，然后再计算即可。限于篇幅这里不再详细说明一维的理论。此外，以下所有数值算例均采用无量纲进行计算，伪弧长控制函数均采用式(17)计算。

2.1   一维Euler方程激波管问题

计算域为[−5，5]，终止时间为${t_{\max }} = 2$，伪弧长控制函数(17)取为$({a_1},{a_2}) = (1,{\text{ }}5)$，边界条件为自由输出边界条件，计算结果如图2所示。

图  2  一维Euler方程激波管图

Figure  2.  Results of shock tube for one dimensional Euler equation

下载: 全尺寸图片 幻灯片

上面的计算均是在150个网格数目下求得的，其中参照解是采用20 000个固定网格下结合五阶特征变换进行计算的。其中MUSCL,WENO-5,WENOCV-5分别代表的是固定网格下用二阶MUSCL，五阶WENO，五阶特征变量WENO格式计算；而PALM-2, PALM-5代表的是在采用二阶伪弧长和五阶伪弧长的格式计算。

从图2(a)中可以看出在固定网格下，高阶格式相对于低阶格式可以更加精确地捕捉到间断，但是从局部放大中可以看出，五阶固定网格在间断附近产生了数值震荡，这是高阶格式带来的弊端，因此可以通过特征变换进行求解，进行特征变换之后，高阶格式消除了局部震荡，整体解的光滑性较好，因此可以经过特征变换得到较为满意的结果。从图2(b)中可以看出，伪弧长算法可以很容易的捕捉到间断面，而且二阶伪弧长算法对间断的捕捉几乎接近于固定网格下五阶的捕捉程度，但是弱于伪弧长五阶算法。因此这说明了高阶伪弧长算法是有效的。图2(c)是网格的运动轨迹，说明伪弧长算法很好地控制了网格的自适应移动，使得网格在梯度较大的区域进行加密，进而捕捉到间断解。

2.2   一维爆炸问题

该算例的计算域为[0,1]，终止时间${t_{\max }} = 0.038$，两套伪弧长控制函数分别取为${({a_1},{a_2})_1} = (10,20)$，${({a_1},{a_2})_2} = (10,80)$，边界条件为固壁反射条件，计算结果如图3所示，其中2-1, 2-2, 5-1, 5-2分别代表二阶和五阶分别采用系数1（${({a_1},{a_2})_1} = (10,20)$）和系数2（${({a_1},{a_2})_2} = (10,80)$）时的计算结果，可以看出，当采用不同的系数时，逼近的分辨率也不同。

图  3  一维爆炸算例计算结果

Figure  3.  Results of one dimensional explosion example

下载: 全尺寸图片 幻灯片

该算例是一个爆炸问题的算例，是检验高精度格式的经典算例。普通低阶格式在网格数目较少的情况下很难捕捉到强间断，而且在计算的过程中，由于稀疏波的传播，压力接近于零，使得计算容易产生负值，导致程序终止。因此，在计算的时候均采用了保正性条件，详细的保正性理论可以参考文献[9]。图中的参照解是在固定网格下用五阶特征变换在网格数为25 000下进行计算得到的。

从上面的计算结果可以看出，二阶固定网格耗散较强，对最高点的捕捉能力较弱；五阶固定网格较二阶提升很大，可以较为清晰的捕捉到多个间断点，而且在间断点附近耗散性非常小；二阶伪弧长算法较固定网格算法有一定提升，固定网格对断点的捕捉能力很弱，耗散性较强，而采用伪弧长算法以后，耗散能力减弱，对间断点的捕捉能力增强；采用五阶伪弧长算法最为满意，各个间断点都可以较为准确的捕捉，而且在间断点邻域的结果和参照解也最为接近。

通过网格的轨迹线图，可以看出在t=0.027时刻附近，两个冲击波相遇，同时伴随着稀疏波的传播，正是因为稀疏波的传播，造成压力过小，接近于零。这样，采用普通的通量差分格式计算会因为计算机的浮点误差和计算格式的数值误差造成物理量为负，从而导致程序的终止。因此需要采用特殊的格式进行修正才能保证程序的进行，本文采用了保正性格式进行限制。

3.   数值精度验证

3.1   一维精度验证

初值条件：$\rho (x,0) = 1 + 0.2\sin (x),{\text{ }}u(x,0) = 0.5,{\text{ }}p(x,0) = 1.0$。

人为解：$\rho (x,t) = 1 + 0.2\sin (x - 0.5t),{\text{ }}u(x,t) = 0.5,{\text{ }}p(x,t) = 1.0$。

伪弧长控制函数取为$({a_1},{a_2}) = (15,15)$，边界条件采用周期性边界条件。

对于范数误差和收敛阶的公式参照文献[8]，其中范数误差在一维情况下取${L_1} =\displaystyle \sum\limits_{i = 1}^N {{\text{|}}{w_i} - w_i^{{\text{exact}}}|} {{\Delta }}x$，而在二维情况下取${L_1} =\displaystyle \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {|{{\text{w}}_{i,j}} - {\text{w}}_{i,j}^{{\text{exact}}}|} } {{\Delta }}x{{\Delta }}y$，${w^{{\text{exact}}}}$代表解析解，收敛阶采取$\dfrac{{\ln ({\varepsilon _{k - 1}}/{\varepsilon _k})}}{{\ln ({{\Delta }}{s_{k - 1}}/{{\Delta }}s{}_k)}}$ ，其中${\varepsilon _k}$表示范数误差，${{\Delta }}s{}_k$表示网格步长。验证结果如表1所示，可以看出，伪弧长算法的范数误差和收敛阶较为满意，采用高阶伪弧长算法，并没有损失太多的精度。表1比较了二阶固定网格FG-2 (Fixed-Grid-2)，五阶固定网格FG-5 (Fixed-Grid-5)，二阶伪弧长PALM-2，五阶伪弧长PALM-5算法下的误差和收敛阶。

表  1  有限体积法与伪弧长算法在不同网格数下的误差和精度

Table  1.  Numerical errors and precision of FVM and PALM changing with grid numbers (example 3)

网格数	L1					Order
	FG-2	PALM-2	FG-5	PALM-5		FG-2	PALM-2	FG-5	PALM-5
40	3.197×10−2	3.185×10−2	4.050×10−5	6.079×10−5
80	9.173×10−3	9.181×10−3	1.021×10−6	1.047×10−6		1.801	1.794	5.310	4.892
160	2.502×10−3	2.500×10−3	3.042×10−8	3.265×10−8		1.874	1.877	5.068	5.003
320	6.804×10−4	6.712×10−4	1.365×10−9	1.675×10−9		1.879	1.897	4.478	4.285

下载: 导出CSV 

| 显示表格

3.2   二维精度验证

采用二维两步化学反应流来测试伪弧长算法的范数误差和收敛阶。对应式(1)的详细表达式如下：

式中：$\ \rho$为密度，$p$为压力，$\gamma$为比热比，$u$、$v$为速度，$\alpha$、$\ \beta$为化学反应量，$Q$为热释放速率，$R$为理想气体常数，${\omega _\alpha }$、${\omega _\beta }$为化学反应变化速率，${k_\alpha }$、${k_\beta }$为反应率常数，${E_\alpha }$、${E_\beta }$为激活能。

采取构造一组人为解进行数值验证。相应的源项需要根据人为解进行修正。其中构造的人为解详细如下：

为了保证计算的完整性，取终止时间T=$2{\text{π}}$，计算域为$[0,2{\text{π]}} \times [0,2{\text{π}}]$，采用周期性边界条件。伪弧长控制参数为$({a_1},{a_2}) = (4,4)$。

表2给出了二维验证计算结果，可以看出，随着计算网格的加密，误差在逐渐减小，而且收敛阶和算法本身较为吻合。通过图4可以看出，高阶有限体积算法对等密度图的捕捉较低阶有限体积更为精细，而且采用伪弧长算法之后，网格自动在物理量梯度较大的区域进行自适应加密，同时也没有失去解的准确性。

表  2  有限体积法与伪弧长算法在不同网格数下的误差和精度

Table  2.  Numerical errors and precision of FVM and PALM changing with grid numbers (Example 4)

Mesh	L1					Order
	FG-2	PALM-2	FG-5	PALM-5		FG-2	PALM-2	FG-5	PALM-5
20×20	1.756	1.626	6.889×10−2	7.224×10−2
40×40	0.597	0.558	2.116×10−3	2.009×10−3		1.558	1.541	5.025	5.169
80×80	0.157	0.140	7.577×10−5	7.564×10−5		1.924	1.997	4.803	4.731
160×160	0.042	0.033	2.365×10−6	2.672×10−6		1.890	2.079	5.002	4.823

下载: 导出CSV 

| 显示表格

图  4  密度云图（$T = 2{\text{π}}$，网格$40 \times 40$）

Figure  4.  Density contours ($T = 2{\text{π}}$, grid $40 \times 40$)

下载: 全尺寸图片 幻灯片

4.   爆轰化学反应流的计算

4.1   一维化学反应流

一维爆炸冲击问题，控制方程取为如下：

伪弧长控制函数取为$({a_1},{a_2}) = (0.5,5)$，边界条件为反射边界，终止时间T=0.1。计算结果如图5所示。

图  5  一维化学反应流密度

Figure  5.  Density of one dimensional chemical reaction flow

下载: 全尺寸图片 幻灯片

图  7  一维化学反应流的网格轨迹

Figure  7.  Mesh trajectory of one dimensional chemical reaction flow

下载: 全尺寸图片 幻灯片

一维化学反应流的参照解是基于5 000个固定网格结合特征变换进行计算的，而二阶格式采用的是120个网格，五阶格式采用的是120个网格。可以看出，在固定网格下，五阶格式对间断的捕捉要优于二阶格式，而经过伪弧长的变化之后，二阶格式也可以较为准确是捕捉到强间断。对解捕捉的最好结果是采用了五阶伪弧长格式，它对网格捕捉的分辨率最高，图中可以看出高阶伪弧长和参照解几乎重合。这充分说明了高阶伪弧长格式的优越性。在本算例中，由于压力的初值变化跨越了10个数量级以上，计算时采用的变量均为原始守恒变量，并没有采用局部特征变换，从而导致在间断附近存在微小震荡。采用了伪弧长算法，并没有完全消除震荡，但是却很明显的减弱了数值震荡，这说明伪弧长算法对于减弱方程的奇异性有着很好的效果。需要注意的是，由于爆炸冲击问题的强烈的奇异性，常规的高阶格式会因为计算物理量为负，而导致程序终止。在本文中，修正式(7)的$\varepsilon = {10^{ - {\text{15}}}}$，并结合保正性算法^[9]，这样才能使得程序可以继续计算。此外可以看到，网格的的轨迹得到了很好的自适应，它在解梯度较大的区域自适应加密。

4.2   二维化学反应流前台阶爆轰问题

计算域为$[0,3] \times [0,1]$。障碍物的区域为$[1,3] \times$$[0,0.4]$，伪弧长控制函数取为$({a_1},{a_2}) = (2,0.2)$，初值条件取ZND^[10]的解析解作为初值。控制方程采用式(27)两步化学反应流。边界条件：除了左右两侧，其他均设为反射边界。左右两侧为流入流出边界。终止时间T=0.22。图6为各方法给出的T=0.22时刻的状态。

图  6  一维化学反应流压力

Figure  6.  Pressure of one dimensional chemical reaction flow

下载: 全尺寸图片 幻灯片

图  8  二维化学反应流图

Figure  8.  Results of two dimensional chemical reaction flow

下载: 全尺寸图片 幻灯片

图6的参照解是采用了$1\;200 \times 400$网格数结合五阶WENO有限体积进行计算的结果。其余计算均在$360 \times 90$网格下进行计算。通过对比可以发现伪弧长算法都较为准确的捕捉到了爆轰反应的细节，对爆炸冲击波阵面都得到了较为精确的捕捉，而且用伪弧长五阶算法的结果和参照解的结果更为接近，并且在局部拐角区域，高阶算法对细节的捕捉更为精确。这说明了高阶伪弧长算法收敛速度较快，因此可以更好地处理爆炸与冲击强间断的问题。

5.   结　语

本文将高精度数值算法和伪弧长算法进行结合，较为详细的介绍了高阶伪弧长算法从推导到应用的过程。针对高阶伪弧长算法在处理通量计算和网格变换之后的物理量重构中的难题，通过引入计算坐标，经过坐标变换将计算过程转换到曲线坐标系下进行计算。通过构造的人为解，验证了高阶伪弧长算法的数值精度。

各算例的计算验证表明：伪弧长算法相较于有限体积算法可以更成功地捕捉到强间断；而且高阶伪弧长算法对强间断的捕捉相对于传统高阶有限体积法和二阶伪弧长算法有更大的优势，它可以以更少的网格数目获得更高的分辨率。

本文采用伪弧长算法结合高阶WENO有限体积格式进行计算，有限体积格式在处理物理量重映过程中的过程较为繁琐，而且对网格移动的限制较为严格，要求每两步网格迭代变化的距离不能太大。因此为了提高计算效率，可以尝试采取有限差分进行计算，有限差分结合动网格和坐标映射可以避免物理量的重映过程，因此可以适当提高计算效率。为了提高计算精度获得高分辨率的结果，而且又不需要很高的网格数目，可以采用高阶伪弧长算法进行计算，它可以以较少的网格数目得到较为满意的结果，从而提高计算效率。