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

許志宇，譚永華，李小明

（1. 西安航天動(dòng)力研究所，陜西 西安 710100；2. 西安航天動(dòng)力研究所液體火箭發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710100；3. 航天推進(jìn)技術(shù)研究院，陜西 西安 710100）

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

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

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

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

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

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

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

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

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

其中

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

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

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

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

2 沖擊波捕捉方法

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

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

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

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

2.2 計(jì)算格式

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

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

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

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

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

半點(diǎn)i±1/2 的值取相鄰整點(diǎn)i 和i±1 的平均值：

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

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

3.1 弱沖擊波問題

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

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

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

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

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

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

圖3 為不同分辨率條件下小波函數(shù)對應(yīng)的配點(diǎn)分布，尺度函數(shù)空間V2的5 個(gè)配點(diǎn){?0.5, ?0.25, 0,0.25, 0.5}均不刪除，因此未在圖中展示。從圖3 可以看出，隨著分辨率的提高，增加的配點(diǎn)集中在沖擊波波陣面、接觸間斷、稀疏波波頭和波尾4 個(gè)梯度較大的位置。

圖 3 計(jì)算用配點(diǎn)分布，t=0.24Fig. 3 Spatial distribution of used collocations, t=0.24

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

3.2 強(qiáng)沖擊波問題

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

初始條件設(shè)為

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

表 2 強(qiáng)沖擊波管計(jì)算參數(shù)Table 2 Computational parameters for strong shock tube

圖 5 強(qiáng)沖擊波t=0.05 時(shí)的氣體密度分布Fig. 5 Density distribution of strong shock at t=0.05

3.3 計(jì)算效率

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

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

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

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

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

4 結(jié) 論

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

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

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

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