传统的冲击波捕捉方法为了获得高分辨率、高精度解，需要非常细密的网格，以便能捕捉清晰的冲击波结构。但是，精度和分辨率越高，计算量则越大，特别是高维问题，除非使用非规则(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和s_j0分别为尺度函数和尺度函数系数，ψ_j和d_j分别为小波函数和小波系数。

尺度空间V和小波空间W互为正交补空间，即${{{V}}_{j + 1}} = {{{V}}_j} + {{{W}}_j}$，因此，当小波函数的尺度因子j取至j_max时，尺度空间V_J的尺度因子J=j_max+1。为了获得式(1)中的系数，需要对函数f(x)进行采样，获得（2^J+1）个配点的值，然后进行离散小波变换。

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

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

其中

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

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

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

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

2.   冲击波捕捉方法

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

2.1   冲击波定位函数

一维守恒方程：

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

式中：$d_k^{{j_{\max }}}$为最细尺度的小波系数，并利用对应的守恒量的范数对$d_{k,\max }^{{j_{\max }}}$归一化。Euler方程组包含多个守恒变量，理论上需要对每个变量进行小波变换，并分别计算冲击波定位函数，并选择最大值作为最终的定位函数值^[1]。这种方法适用范围广，但计算复杂，且增加了计算量，并不利于在工程应用中推广。考虑到冲击波问题和一般Riemann问题的共性，选择守恒变量密度ρ进行小波变换、生成适应性网格并构造冲击波定位函数。这有利于简化计算、节约计算量。但对于强冲击波问题，由于人工黏性不足，并不能有效抑制数值振荡。由于式(7)给出的定位函数为不大于1的正数，因此，采用幂函数形式的定位函数能够控制黏性分布的宽度：

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

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

图  1  冲击波与对应的小波系数和定位函数

Figure  1.  Shock and the corresponding wavelet coefficients and shock locator functions

下载: 全尺寸图片 幻灯片

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或更大值时可以得到较好的数值解。以N₀表示有效配点数量，N表示使用的配点数量，具体数值如表1所示。

表  1  弱冲击波管计算参数

Table  1.  Computational parameters for weak shock tube

J	N	α	ε/10−4
9	513	1.5	1.0
10	1 025	1.5	1.0
11	2 049	1.5	1.0

下载: 导出CSV 

| 显示表格

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

图  2  弱冲击波t=0.24时的气体密度分布

Figure  2.  Density distribution of weak shock at t=0.24

下载: 全尺寸图片 幻灯片

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

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

图  3  计算用配点分布，t=0.24

Figure  3.  Spatial distribution of used collocations, t=0.24

下载: 全尺寸图片 幻灯片

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

图  4  冲击波波阵面位置随时间变化

Figure  4.  Positions of the shock front with time

下载: 全尺寸图片 幻灯片

3.2   强冲击波问题

初始条件设为

当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（N₀=4 097）时，又增加了60个配点，压缩比（N₀/N）分别为3.61、6.00和10.20，分辨率越高，压缩效率越高。

表  2  强冲击波管计算参数

Table  2.  Computational parameters for strong shock tube

J	N	α	ε/10−5
10	1 025	0.3	1.0
11	2 049	0.3	1.0
12	4 097	0.3	1.0

下载: 导出CSV 

| 显示表格

图  5  强冲击波t=0.05时的气体密度分布

Figure  5.  Density distribution of strong shock at t=0.05

下载: 全尺寸图片 幻灯片

3.3   计算效率

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

尺度因子选择为9～14，对应6种分辨率，计算网格数范围为512～16 384。自适应小波配点法的相对计算时间如图6所示。若CFL条件与自适应小波配点法相同，采用WENO格式时冲击波后出现数值振荡，需要通过减小时间步长消除，取$\Delta t_{\rm{WENO}}=0.75\Delta t_{\rm{awcm}}$时，基本无数值振荡，但导致更长的计算时间。

图  6  自适应小波配点法相对迎风、ENO、WENO格式的计算时间

Figure  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时三种格式密度分布对比

Figure  7.  Spatial distributions of density through three schemes with J=12

下载: 全尺寸图片 幻灯片

4.   结　论

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

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

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

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