沖擊波的小波數(shù)值計算方法*

許志宇，譚永華，李小明

（1. 西安航天動力研究所，陜西 西安 710100；2. 西安航天動力研究所液體火箭發(fā)動機技術(shù)重點實驗室，陜西 西安 710100；3. 航天推進技術(shù)研究院，陜西 西安 710100）

傳統(tǒng)的沖擊波捕捉方法為了獲得高分辨率、高精度解，需要非常細密的網(wǎng)格，以便能捕捉清晰的沖擊波結(jié)構(gòu)。但是，精度和分辨率越高，計算量則越大，特別是高維問題，除非使用非規(guī)則(irregular)網(wǎng)格，否則因計算量巨大的原因?qū)⒑茈y實現(xiàn)[1]。小波數(shù)值方法是基于多分辨分析(multi resolution analysis,MRA)發(fā)展的新方法，由于小波函數(shù)具有緊支撐特性，因此能夠?qū)α鲌鰯?shù)據(jù)進行壓縮，生成捕捉流場局部結(jié)構(gòu)的適應(yīng)性網(wǎng)格，適合描述局部流動特征顯著的問題[2]。

小波數(shù)值方法主要有兩類，即小波-迦遼金(wavelet-Galerkin)[3]和小波配點(wavelet collocation)法[4-6]。小波-迦遼金法不適合處理非線性算子和任意邊界條件問題，而自適應(yīng)小波配點法在這兩方面均具有優(yōu)勢，特別是二代小波，在真實物理域中進行變換，可以方便處理任意邊界條件，因此發(fā)展迅速[7-8]。

Harten[9]最早將小波用于壓縮算法，從而減少高精度計算格式所需的網(wǎng)格數(shù)；Bürger 等[10]、王昱[11]和孫陽等[12]采用小波壓縮聯(lián)合WENO 或TVD 等格式，各自實現(xiàn)了傳統(tǒng)高精度格式的高分辨率計算，表明自適應(yīng)多分辨率格式為傳統(tǒng)的沖擊波捕捉格式節(jié)約了計算量。

由于小波數(shù)值方法基于嵌套的動態(tài)自適應(yīng)網(wǎng)格，而在動態(tài)網(wǎng)格上使用迎風(fēng)格式或其他高精度格式非常復(fù)雜[1]。趙勇等[13-14]提出了雙重小波收縮法和第二黏性法，避免了在動態(tài)網(wǎng)格上使用迎風(fēng)格式。Regele 等[1,15-17]、Kassoy 等[18]和Schneider 等[19]根據(jù)傳統(tǒng)的構(gòu)造沖擊波定位函數(shù)和人工黏性的方法，利用最細層小波系數(shù)的絕對值和守恒變量構(gòu)造沖擊波定位函數(shù)，然后利用沖擊波定位函數(shù)控制人工黏性的大小和分布，成功求解了沖擊波、爆轟波等問題。構(gòu)造的沖擊波定位函數(shù)能夠準確判斷沖擊波位置，但由于需要對所有守恒量進行小波分解和計算范數(shù)，因此計算稍顯復(fù)雜，不利于解決復(fù)雜的工程問題。

本文給出一種使用簡單、計算穩(wěn)定的沖擊波定位函數(shù)和控制人工黏性的方法；詳細介紹利用小波多尺度分解生成自適應(yīng)網(wǎng)格、構(gòu)造沖擊波定位函數(shù)、控制人工黏性的原理和過程，并對強、弱沖擊波管問題進行計算驗證，分析計算沖擊波的能力和特點。

1 自適應(yīng)小波配點法

由多分辨分析理論可知，任意函數(shù)f(x)可以多尺度分解為如下形式：式中：φj0和sj0分別為尺度函數(shù)和尺度函數(shù)系數(shù)，ψj和dj分別為小波函數(shù)和小波系數(shù)。

小波配點法以小波分解為基礎(chǔ)，通過式(1)計算偏微分方程中的空間導(dǎo)數(shù)項，從而將偏微分方程化為關(guān)于時間的常微分方程；然后利用Euler 法、R-K 法等常微分方程的數(shù)值方法步進求解。

沖擊波等具有局部特性的問題，小波系數(shù)的絕對值僅在變化急劇的局部區(qū)域較大，而在連續(xù)區(qū)域很小。根據(jù)小波系數(shù)相對于某個閾值ε（ε＞0）的大小，將式(1)分為兩部分：

其中

對于正則方程，舍去式（4）中的小波，誤差的上限滿足[7]

因此能夠保障計算精度。由于小波和配點一一對應(yīng)，因此舍去小波的同時刪除了對應(yīng)位置的配點，從而生成自適應(yīng)網(wǎng)格。

生成自適應(yīng)網(wǎng)格的過程需要注意兩個細節(jié)：(1) 由于局部流動結(jié)構(gòu)發(fā)生變化和移動，因此，根據(jù)當(dāng)前流場生成的網(wǎng)格應(yīng)具有預(yù)測下一時刻流場特性的能力，所以需要適當(dāng)保留局部特征附近的網(wǎng)格點[7]；(2) 為了能夠在自適應(yīng)網(wǎng)格上進行小波變換，需要保留小波變換所需的上一層的網(wǎng)格點。

為了方便處理邊界和減少離散小波變換的計算量，本文采用一種二代小波—提升插值小波進行小波變換，預(yù)測和更新分別利用兩側(cè)各兩個點，具體可參考文獻[7]。

2 沖擊波捕捉方法

利用小波對流場進行多尺度分解，對于最細尺度，系數(shù)最大的小波對應(yīng)沖擊波位置，并且系數(shù)越大表明梯度越大。因此通過最細尺度的小波系數(shù)構(gòu)造沖擊波定位函數(shù)。

2.1 沖擊波定位函數(shù)

一維守恒方程：式中：U 為守恒變量，F(xiàn) 為通量，υ 為人工黏性，Φ 為沖擊波定位函數(shù)（或通量限制器），等號右側(cè)為人工黏性項。最細尺度小波系數(shù)最大的區(qū)域能反映沖擊波位置，位置誤差為最小網(wǎng)格單元。而對于j＜jmax，沖擊波定位函數(shù)不能為判斷沖擊波位置提供準確信息。因此選擇最細尺度jmax對應(yīng)的小波系數(shù)構(gòu)造沖擊波定位函數(shù)：[1]

當(dāng)α＜1 時，沖擊波波陣面兩側(cè)的定位函數(shù)值增大，數(shù)值控制的沖擊波區(qū)域變寬；當(dāng)α＞1 時，波陣面兩側(cè)定位函數(shù)值減小，區(qū)域變窄。

圖 1 沖擊波與對應(yīng)的小波系數(shù)和定位函數(shù)Fig. 1 Shock and the corresponding wavelet coefficients and shock locator functions

式(7)和式(8)均未給出尺度j＜jmax對應(yīng)配點的Φ 值。由于j＜jmax和jmax對應(yīng)的配點滿足二分性，因此j＜jmax對應(yīng)配點的Φ 值可通過j=jmax對應(yīng)的Φ 值插值給定。對于如圖1(a)所示的沖擊波，參照閾值ε 刪除部分配點后，保留的小波系數(shù)分布如圖1(b)所示，箭頭代表小波系數(shù)，箭頭尾部縱坐標為尺度因子，箭頭長短代表小波系數(shù)的大小。從圖1(b)中可以看出，在沖擊波波陣面附近小波系數(shù)較大，因此保留了大量配點，而在遠離波陣面的區(qū)域僅保留了少量的配點。按式(8)計算沖擊波定位函數(shù)，并取α=1/3，如圖1(c)所示，可見集中在波陣面附近的7 個點處Φ 值較大，遠離波陣面處，Φ 迅速趨于0。Φ 值較大的這7 個點構(gòu)成數(shù)值計算的沖擊波區(qū)，人工黏性將在該區(qū)域內(nèi)起作用。

2.2 計算格式

通量的空間導(dǎo)數(shù)選擇簡單的中心差分格式，二階差分格式可以計算沖擊波，并作為基本計算格式，但其精度較低。當(dāng)采用四階或更高階格式時，可以使接觸間斷的精度更高。時間積分均采用一階步進格式。

人工黏性項的離散采用中心差分格式，因此式(6)的二階計算格式為[20]

式中：Δx 為空間步長，由尺度因子J 決定，U 為守恒變量。根據(jù)非線性穩(wěn)定性條件確定的最小黏性，式(9)標準的顯示格式為

式中：a 為對流速度。對于Euler 方程組，為了計算簡便，a 取Jacobian 矩陣的最大特征值c+|u|，其中c 為當(dāng)?shù)芈曀?，c=(γp/ρ)1/2，γ 為比熱比，p 為壓力，ρ 為密度；u 為速度[20]。這種方法的數(shù)值穩(wěn)定性好，但對接觸間斷有一定抹平[1]。比較式(9)和(10)的人工黏性項，可得人工黏性表達式為

為了保證人工黏性僅在大梯度區(qū)起作用，因此將沖擊波定位函數(shù)Φ 與人工黏性υ 聯(lián)合，由Φ 控制υ 的大小分布，形成對整個計算域通用的人工黏性項：

半點i±1/2 的值取相鄰整點i 和i±1 的平均值：

3 計算結(jié)果和分析

為了驗證計算沖擊波的能力和特點，選擇可得精確解的一維沖擊波管問題進行計算，控制方程為一維Euler 方程，狀態(tài)方程采用理想氣體狀態(tài)方程，氣體比熱比γ=1.4。

3.1 弱沖擊波問題

沖擊波管初始條件設(shè)為

選取不同的尺度因子進行計算，當(dāng)J＜9 時，沖擊波抹平嚴重，J 分別取9、10、11 或更大值時可以得到較好的數(shù)值解。以N0表示有效配點數(shù)量，N 表示使用的配點數(shù)量，具體數(shù)值如表1所示。

表 1 弱沖擊波管計算參數(shù)Table 1 Computational parameters for weak shock tube

當(dāng)t=0.24 時，氣體密度分布分別如圖2 所示。圖2 表明，分辨率越高，計算的沖擊波越接近Riemann 解。在相同分辨率條件下，接觸間斷處的抹平較為明顯，精度較低，采用四階或更高階差分格式可以在一定程度上提高接觸間斷處的精度。

圖 2 弱沖擊波t=0.24 時的氣體密度分布Fig. 2 Density distribution of weak shock at t=0.24

基礎(chǔ)網(wǎng)格點數(shù)量N0=2J+1，J 增加1，分辨率和基礎(chǔ)網(wǎng)格數(shù)量增加一倍，但所用配點數(shù)量并未相應(yīng)增加一倍。J=9（N0=513）時，計算中使用202 個配點，當(dāng)J=10（N0=1 025）時，增加了91 個配點，而J=11（N0=2 049）時又增加了81 個配點，壓縮比（N0/N）分別為2.54、3.50 和5.51，說明分辨率越高，數(shù)據(jù)壓縮程度和相對計算效率越高。

圖3 為不同分辨率條件下小波函數(shù)對應(yīng)的配點分布，尺度函數(shù)空間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

在不同分辨率條件下，沖擊波定位函數(shù)顯示的沖擊波位置如圖4 所示。可以看出，沖擊波位置隨時間的變化與Riemann 解一致。由二代小波分解的特點可知，最細尺度的小波系數(shù)判別沖擊波位置，產(chǎn)生的最大位置誤差為一個網(wǎng)格寬度。因此，分辨率越高，位置誤差越小。

3.2 強沖擊波問題

圖 4 沖擊波波陣面位置隨時間變化Fig. 4 Positions of the shock front with time

初始條件設(shè)為

當(dāng)J＜10，沖擊波抹平嚴重，J≥10 時，可以得到較好的結(jié)果，因此分別取J=10、11 和12 進行計算。對于強沖擊波問題，需令參數(shù)α＜1，使得沖擊波區(qū)的寬度適當(dāng)增加，從而能夠有效抑制振蕩。表2 給出了J=10, 11, 12 時使用的計算參數(shù)，圖5 對應(yīng)給出了t=0.05 時刻的密度分布。J=10 時，計算使用284 個配點；當(dāng)J=11 時，分辨率提升一倍，使用的配點增加了58 個；而J=12（N0=4 097）時，又增加了60 個配點，壓縮比（N0/N）分別為3.61、6.00 和10.20，分辨率越高，壓縮效率越高。

表 2 強沖擊波管計算參數(shù)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 計算效率

為了說明自適應(yīng)小波配點法在計算效率上的優(yōu)勢，比較幾種傳統(tǒng)的沖擊波捕捉格式在相同分辨率和CFL 條件下的運算時間。分別選取基于Steger-Warming 矢通量分裂的一階迎風(fēng)（Up-wind）格式、三階ENO 格式和五階WENO 格式對3.1 節(jié)和3.2 節(jié)中的沖擊波管問題進行計算。四種格式的計算時間記為tAWCM、tUp-wind、tENO和tWENO，自適應(yīng)小波配點法相對于迎風(fēng)格式、ENO 格式和WENO 格式的計算時間分別記為τA_U、τA_E和τA_W，定義為

圖 6 自適應(yīng)小波配點法相對迎風(fēng)、ENO、WENO 格式的計算時間Fig. 6 Relative computational time costs of AWCM to Up-wind, ENO and WENO schemes

相對于一階迎風(fēng)格式，自適應(yīng)小波配點法節(jié)省40%～70% 的計算時間；相對于ENO 和WENO 高精度格式，能節(jié)約超過70%的計算時間。分辨率提升一倍時，相對計算效率更高。

需要指出的是，本文的自適應(yīng)配點法沒有考慮流動特性，采用了簡單的中心差分格式，對于局部間斷問題，不利于獲得高精度解。提升分辨率后，計算結(jié)果的質(zhì)量低于ENO 和WENO 格式，特別是在接觸間斷附近，如圖7 所示。因此，后續(xù)可對考慮流動特性并利用小波構(gòu)造高精度格式進行研究。

圖 7 J=12 時三種格式密度分布對比Fig. 7 Spatial distributions of density through three schemes with J=12

4 結(jié) 論

基于自適應(yīng)小波配點法和人工黏性技術(shù)構(gòu)造了沖擊波數(shù)值計算格式；利用小波系數(shù)構(gòu)建了指數(shù)形式的沖擊波定位函數(shù)，并用以控制人工黏性的分布；對強/弱沖擊波管問題進行計算，得出如下結(jié)論：

（1）利用小波系數(shù)構(gòu)建的指數(shù)形式的沖擊波定位函數(shù)能夠準確捕捉?jīng)_擊波位置，最大位置誤差為一個網(wǎng)格寬度；

（2）沖擊波定位函數(shù)的指數(shù)可以控制沖擊波區(qū)域的寬度，指數(shù)越小，沖擊波區(qū)域越寬，計算越穩(wěn)定；

（3）自適應(yīng)配點法利用小波閾值濾波刪除大量網(wǎng)格點，比傳統(tǒng)方法計算高效，且分辨率越高，相對計算效率越高。