冲击波的小波数值计算方法*

许志宇，谭永华，李小明

（1. 西安航天动力研究所，陕西 西安 710100；2. 西安航天动力研究所液体火箭发动机技术重点实验室，陕西 西安 710100；3. 航天推进技术研究院，陕西 西安 710100）

传统的冲击波捕捉方法为了获得高分辨率、高精度解，需要非常细密的网格，以便能捕捉清晰的冲击波结构。但是，精度和分辨率越高，计算量则越大，特别是高维问题，除非使用非规则(irregular)网格，否则因计算量巨大的原因将很难实现[1]。小波数值方法是基于多分辨分析(multi resolution analysis,MRA)发展的新方法，由于小波函数具有紧支撑特性，因此能够对流场数据进行压缩，生成捕捉流场局部结构的适应性网格，适合描述局部流动特征显著的问题[2]。

小波数值方法主要有两类，即小波-迦辽金(wavelet-Galerkin)[3]和小波配点(wavelet collocation)法[4-6]。小波-迦辽金法不适合处理非线性算子和任意边界条件问题，而自适应小波配点法在这两方面均具有优势，特别是二代小波，在真实物理域中进行变换，可以方便处理任意边界条件，因此发展迅速[7-8]。

Harten[9]最早将小波用于压缩算法，从而减少高精度计算格式所需的网格数；Bürger 等[10]、王昱[11]和孙阳等[12]采用小波压缩联合WENO 或TVD 等格式，各自实现了传统高精度格式的高分辨率计算，表明自适应多分辨率格式为传统的冲击波捕捉格式节约了计算量。

由于小波数值方法基于嵌套的动态自适应网格，而在动态网格上使用迎风格式或其他高精度格式非常复杂[1]。赵勇等[13-14]提出了双重小波收缩法和第二黏性法，避免了在动态网格上使用迎风格式。Regele 等[1,15-17]、Kassoy 等[18]和Schneider 等[19]根据传统的构造冲击波定位函数和人工黏性的方法，利用最细层小波系数的绝对值和守恒变量构造冲击波定位函数，然后利用冲击波定位函数控制人工黏性的大小和分布，成功求解了冲击波、爆轰波等问题。构造的冲击波定位函数能够准确判断冲击波位置，但由于需要对所有守恒量进行小波分解和计算范数，因此计算稍显复杂，不利于解决复杂的工程问题。

本文给出一种使用简单、计算稳定的冲击波定位函数和控制人工黏性的方法；详细介绍利用小波多尺度分解生成自适应网格、构造冲击波定位函数、控制人工黏性的原理和过程，并对强、弱冲击波管问题进行计算验证，分析计算冲击波的能力和特点。

1 自适应小波配点法

由多分辨分析理论可知，任意函数f(x)可以多尺度分解为如下形式：式中：φj0和sj0分别为尺度函数和尺度函数系数，ψj和dj分别为小波函数和小波系数。

小波配点法以小波分解为基础，通过式(1)计算偏微分方程中的空间导数项，从而将偏微分方程化为关于时间的常微分方程；然后利用Euler 法、R-K 法等常微分方程的数值方法步进求解。

冲击波等具有局部特性的问题，小波系数的绝对值仅在变化急剧的局部区域较大，而在连续区域很小。根据小波系数相对于某个阈值ε（ε＞0）的大小，将式(1)分为两部分：

其中

对于正则方程，舍去式（4）中的小波，误差的上限满足[7]

因此能够保障计算精度。由于小波和配点一一对应，因此舍去小波的同时删除了对应位置的配点，从而生成自适应网格。

生成自适应网格的过程需要注意两个细节：(1) 由于局部流动结构发生变化和移动，因此，根据当前流场生成的网格应具有预测下一时刻流场特性的能力，所以需要适当保留局部特征附近的网格点[7]；(2) 为了能够在自适应网格上进行小波变换，需要保留小波变换所需的上一层的网格点。

为了方便处理边界和减少离散小波变换的计算量，本文采用一种二代小波—提升插值小波进行小波变换，预测和更新分别利用两侧各两个点，具体可参考文献[7]。

2 冲击波捕捉方法

利用小波对流场进行多尺度分解，对于最细尺度，系数最大的小波对应冲击波位置，并且系数越大表明梯度越大。因此通过最细尺度的小波系数构造冲击波定位函数。

2.1 冲击波定位函数

一维守恒方程：式中：U 为守恒变量，F 为通量，υ 为人工黏性，Φ 为冲击波定位函数（或通量限制器），等号右侧为人工黏性项。最细尺度小波系数最大的区域能反映冲击波位置，位置误差为最小网格单元。而对于j＜jmax，冲击波定位函数不能为判断冲击波位置提供准确信息。因此选择最细尺度jmax对应的小波系数构造冲击波定位函数：[1]

当α＜1 时，冲击波波阵面两侧的定位函数值增大，数值控制的冲击波区域变宽；当α＞1 时，波阵面两侧定位函数值减小，区域变窄。

图 1 冲击波与对应的小波系数和定位函数Fig. 1 Shock and the corresponding wavelet coefficients and shock locator functions

式(7)和式(8)均未给出尺度j＜jmax对应配点的Φ 值。由于j＜jmax和jmax对应的配点满足二分性，因此j＜jmax对应配点的Φ 值可通过j=jmax对应的Φ 值插值给定。对于如图1(a)所示的冲击波，参照阈值ε 删除部分配点后，保留的小波系数分布如图1(b)所示，箭头代表小波系数，箭头尾部纵坐标为尺度因子，箭头长短代表小波系数的大小。从图1(b)中可以看出，在冲击波波阵面附近小波系数较大，因此保留了大量配点，而在远离波阵面的区域仅保留了少量的配点。按式(8)计算冲击波定位函数，并取α=1/3，如图1(c)所示，可见集中在波阵面附近的7 个点处Φ 值较大，远离波阵面处，Φ 迅速趋于0。Φ 值较大的这7 个点构成数值计算的冲击波区，人工黏性将在该区域内起作用。

2.2 计算格式

通量的空间导数选择简单的中心差分格式，二阶差分格式可以计算冲击波，并作为基本计算格式，但其精度较低。当采用四阶或更高阶格式时，可以使接触间断的精度更高。时间积分均采用一阶步进格式。

人工黏性项的离散采用中心差分格式，因此式(6)的二阶计算格式为[20]

式中：Δx 为空间步长，由尺度因子J 决定，U 为守恒变量。根据非线性稳定性条件确定的最小黏性，式(9)标准的显示格式为

式中：a 为对流速度。对于Euler 方程组，为了计算简便，a 取Jacobian 矩阵的最大特征值c+|u|，其中c 为当地声速，c=(γp/ρ)1/2，γ 为比热比，p 为压力，ρ 为密度；u 为速度[20]。这种方法的数值稳定性好，但对接触间断有一定抹平[1]。比较式(9)和(10)的人工黏性项，可得人工黏性表达式为

为了保证人工黏性仅在大梯度区起作用，因此将冲击波定位函数Φ 与人工黏性υ 联合，由Φ 控制υ 的大小分布，形成对整个计算域通用的人工黏性项：

半点i±1/2 的值取相邻整点i 和i±1 的平均值：

3 计算结果和分析

为了验证计算冲击波的能力和特点，选择可得精确解的一维冲击波管问题进行计算，控制方程为一维Euler 方程，状态方程采用理想气体状态方程，气体比热比γ=1.4。

3.1 弱冲击波问题

冲击波管初始条件设为

选取不同的尺度因子进行计算，当J＜9 时，冲击波抹平严重，J 分别取9、10、11 或更大值时可以得到较好的数值解。以N0表示有效配点数量，N 表示使用的配点数量，具体数值如表1所示。

表 1 弱冲击波管计算参数Table 1 Computational parameters for weak shock tube

当t=0.24 时，气体密度分布分别如图2 所示。图2 表明，分辨率越高，计算的冲击波越接近Riemann 解。在相同分辨率条件下，接触间断处的抹平较为明显，精度较低，采用四阶或更高阶差分格式可以在一定程度上提高接触间断处的精度。

图 2 弱冲击波t=0.24 时的气体密度分布Fig. 2 Density distribution of weak shock at t=0.24

基础网格点数量N0=2J+1，J 增加1，分辨率和基础网格数量增加一倍，但所用配点数量并未相应增加一倍。J=9（N0=513）时，计算中使用202 个配点，当J=10（N0=1 025）时，增加了91 个配点，而J=11（N0=2 049）时又增加了81 个配点，压缩比（N0/N）分别为2.54、3.50 和5.51，说明分辨率越高，数据压缩程度和相对计算效率越高。

图3 为不同分辨率条件下小波函数对应的配点分布，尺度函数空间V2的5 个配点{−0.5, −0.25, 0,0.25, 0.5}均不删除，因此未在图中展示。从图3 可以看出，随着分辨率的提高，增加的配点集中在冲击波波阵面、接触间断、稀疏波波头和波尾4 个梯度较大的位置。

图 3 计算用配点分布，t=0.24Fig. 3 Spatial distribution of used collocations, t=0.24

在不同分辨率条件下，冲击波定位函数显示的冲击波位置如图4 所示。可以看出，冲击波位置随时间的变化与Riemann 解一致。由二代小波分解的特点可知，最细尺度的小波系数判别冲击波位置，产生的最大位置误差为一个网格宽度。因此，分辨率越高，位置误差越小。

3.2 强冲击波问题

图 4 冲击波波阵面位置随时间变化Fig. 4 Positions of the shock front with time

初始条件设为

当J＜10，冲击波抹平严重，J≥10 时，可以得到较好的结果，因此分别取J=10、11 和12 进行计算。对于强冲击波问题，需令参数α＜1，使得冲击波区的宽度适当增加，从而能够有效抑制振荡。表2 给出了J=10, 11, 12 时使用的计算参数，图5 对应给出了t=0.05 时刻的密度分布。J=10 时，计算使用284 个配点；当J=11 时，分辨率提升一倍，使用的配点增加了58 个；而J=12（N0=4 097）时，又增加了60 个配点，压缩比（N0/N）分别为3.61、6.00 和10.20，分辨率越高，压缩效率越高。

表 2 强冲击波管计算参数Table 2 Computational parameters for strong shock tube

图 5 强冲击波t=0.05 时的气体密度分布Fig. 5 Density distribution of strong shock at t=0.05

3.3 计算效率

为了说明自适应小波配点法在计算效率上的优势，比较几种传统的冲击波捕捉格式在相同分辨率和CFL 条件下的运算时间。分别选取基于Steger-Warming 矢通量分裂的一阶迎风（Up-wind）格式、三阶ENO 格式和五阶WENO 格式对3.1 节和3.2 节中的冲击波管问题进行计算。四种格式的计算时间记为tAWCM、tUp-wind、tENO和tWENO，自适应小波配点法相对于迎风格式、ENO 格式和WENO 格式的计算时间分别记为τA_U、τA_E和τA_W，定义为

图 6 自适应小波配点法相对迎风、ENO、WENO 格式的计算时间Fig. 6 Relative computational time costs of AWCM to Up-wind, ENO and WENO schemes

相对于一阶迎风格式，自适应小波配点法节省40%～70% 的计算时间；相对于ENO 和WENO 高精度格式，能节约超过70%的计算时间。分辨率提升一倍时，相对计算效率更高。

需要指出的是，本文的自适应配点法没有考虑流动特性，采用了简单的中心差分格式，对于局部间断问题，不利于获得高精度解。提升分辨率后，计算结果的质量低于ENO 和WENO 格式，特别是在接触间断附近，如图7 所示。因此，后续可对考虑流动特性并利用小波构造高精度格式进行研究。

图 7 J=12 时三种格式密度分布对比Fig. 7 Spatial distributions of density through three schemes with J=12

4 结 论

基于自适应小波配点法和人工黏性技术构造了冲击波数值计算格式；利用小波系数构建了指数形式的冲击波定位函数，并用以控制人工黏性的分布；对强/弱冲击波管问题进行计算，得出如下结论：

（1）利用小波系数构建的指数形式的冲击波定位函数能够准确捕捉冲击波位置，最大位置误差为一个网格宽度；

（2）冲击波定位函数的指数可以控制冲击波区域的宽度，指数越小，冲击波区域越宽，计算越稳定；

（3）自适应配点法利用小波阈值滤波删除大量网格点，比传统方法计算高效，且分辨率越高，相对计算效率越高。