基于自適應小波配點法的水擊數值計算

許志宇，譚永華，李小明

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

0 引言

高速流動的液體因閥門突然關閉或泵驟停會產生水擊現象。如果水擊壓力過高、劇烈，可能造成液壓元件或管路破壞、液體泄漏等問題。

水擊的數值計算方法包括有限元法、有限差分法、有限體積法和譜方法等[1-4]。其中，譜方法選擇適當的基函數逼近計算域內的函數，然后結合加權余量法求解偏微分方程，具有精度高、收斂性好的優點。但不論是Fourier級數、Chebyshev或Legendre多項式，在計算域內均不具有緊支撐特性，當求解水擊這一類具有局部奇異性問題時，分辨率越高，產生的數值振蕩越嚴重，需要配合使用濾波才能消除[1]。

小波數值方法是基于計算流體力學中精確計算湍流和捕捉激波發展起來的高精度高分辨率方法，屬于譜方法范疇，主要有小波-迦遼金法(Wavelet-Galerkin Method)和小波配點法(Wavelet Collocation Method)。由于小波函數具有緊支撐特性，因此相對于Fourier等基函數，能有效抑制數值振蕩；此外，根據多分辨分析，可以采用自適應方法自動提高或降低分辨率，在保障精度的前提下可減少計算量，適合求解局部奇異性問題。文獻[5]詳細介紹了利用3次樣條小波有限元法求解水擊問題的理論，但并未開展數值計算。文獻[6]將區間樣條小波應用于太陽能帆板的軌跡優化。文獻[7]詳細討論了利用樣條小波求解初邊值問題，文獻[8]利用小波數值方法求解多維偏微分方程，文獻[9-14]利用自適應小波配點法求解雙曲守恒律方程，用于計算激波、爆轟波等間斷問題，精度和分辨率均較高，文獻[15-16]將小波方法應用于湍流數值模擬，可有效降低計算量。

本文采用自適應小波配點法(Adaptive Wavelet Collocation Method, AWCM)對一維管路水擊進行數值計算，給出了小波數值方法的原理、計算過程和結果，分析了小波數值方法計算水擊問題的能力和特點。

1 自適應小波配點法

根據多分辨分析(MultiResolution Analysis, MRA)理論可知，尺度函數φ(x)伸縮平移變換張成的閉子空間Vj構成一組多分辨分析，尺度函數對應的小波函數ψ(x)通過伸縮平移變換張成的小波空間Wj構成尺度空間Vj的正交補空間。可表示為

φj,k(x)=2j/2φ(2jx-k)

(1)

ψj,k(x)=2j/2(2jx-k)

(2)

Vj=Vj-1⊕Wj-1

(3)

設函數f(x)∈L2(R)，根據式(3)可得f(x)的多尺度分解

(4)

式中s，d分別為尺度函數系數和小波系數。小波系數的絕對值越大表明f(x)在相應尺度的分量越大。

實際上不可能對f(x)作無限分解，因此尺度因子j只能取至某個上限J-1，使分辨率達到2-J，式(4)則近似表示為

(5)

對于水擊這一類問題，解通常僅在局部區域不連續或梯度較大，而在其他區域相對平緩，相應的小波系數均非常小。因此式(5)可以根據小波系數相對于閾值ε的大小分為兩部分[7]

f(x)=f≥(x)+f<(x)

(6)

其中

(7)

(8)

Donoho(1992年)證明，對于正則方程，舍去式(8)系數較小的小波，逼近的誤差上限為

(9)

根據小波系數的大小判斷變化劇烈的局部位置，并選擇適當的閾值對局部區域增加或降低分辨率的方法稱為自適應算法。Mallat分解首先從最細尺度向尺度函數空間分解，得到尺度函數系數，然后根據雙尺度方程

(10)

(11)

利用快速算法，逐步計算上一尺度的尺度函數系數和小波系數，如果小波系數的絕對值大于設定的閾值，則保留該小波，如果小于設定的閾值，則刪除該小波[5-8]。這種自適應算法效率高，但對于求解發展方程仍有不足。由于隨著時間推進，奇異位置發生移動，為了更好捕捉奇異位置和減小誤差，應適當保留奇異位置附近的配點，使得奇異位置移動后仍能保持分辨率[7]。

小波配點法以小波分解為基礎，通過式(5)計算偏微分方程中的空間導數項，一階導數

(12)

從而化偏微分方程為關于時間的常微分方程[5]。然后利用Euler法、R-K法等常微分方程的數值方法求解。

2 一維流動方程和計算格式

2.1 方程半離散格式

假設流體在管路內作一維絕熱流動，考慮壁面摩擦，無量綱的一維流動方程為[1]

(13)

其中

P=p/pT

U=u/a

X=x/L

τ=ta/L

式中：P，U，X及τ分別為無量綱的壓力、速度、位置及時間;系數ε=p0/(ρa)，μ=fL/(2D);pT和p0分別為入口壓力和管中初始壓力;L和D分別為管長和管徑;ρ和a分別為流體密度和聲速。

時間離散采用Euler法，因此式(13)的半離散格式為

(14)

式中壓力和速度采用小波分解表達式

(15)

(16)

根據式(12)計算出壓力和速度的空間一階導數，一并代回式(14)中，選取合適的時間步長Δτ，求解下一時刻的壓力和速度分布。重復上述過程便可求解各時刻的物理場。

2.2 邊界條件設置

小波變換定義在整個實數域(-∞，+∞)，而實際問題均屬于有限域[L,R]。截斷(-∞，L)和(R，+∞)范圍內的小波，會在邊界處造成較大誤差，使得邊界處小波系數產生異常、精度降低，即邊界效應。本文采用邊界延拓的方法保障邊界精度，使邊界效應不在真實的物理邊界處發生。

邊界延拓的關鍵在于給延拓的部分賦值。對于實際問題，根據邊界條件給定，使延拓部分具有相應邊界條件的物理意義。閥門關閉形成剛性壁面條件，即壁面處的速度為0，壓力通過式(14)計算，而延拓部分與關于壁面對稱部分的速度大小相等、方向相反，壓力和溫度等標量大小相等。對于壓力入口，延拓區域的壓力設定為與入口壓力相同，速度則通過式(14)計算。

尺度函數和小波函數的支集一般為整數域，為了方便進行小波分解和建立統一的計算格式，計算域[L,R]設定為正整數作為基礎域，對于實際問題進行坐標變換。本文選用的小波為支集為4的四階三次樣條小波，基礎域為[0,16]。

3 計算結果和分析

選擇水平直管路閥門瞬間關閉產生水擊的模型(圖1)進行計算。管路參數和初邊值如表1所示[1]，無量綱的初始流速U0=5.67×10-3，當t>0，閥端流速設置為0，以模擬閥門瞬間關閉。

圖1 一維水擊模型Fig.1 One-dimensional water-hammer model

3.1 數值振蕩和誤差分析

首先針對無損水擊進行計算。由于不考慮管路沿程損失，管路中產生的水擊為周期T=4L/a的方波，根據如科夫斯基公式，可求得完全水擊壓力的理論值

Δp=Δuρa

(17)

由于無壓力損失，因此管中穩態壓力等于入口壓力pT，初始流速假設與摩擦系數f=0.018時相同，即U0=5.67×10-3，根據式(17)可得Δp≈0.544，則水擊的波峰和波谷分別為1.544和0.456。

調整閾值ε時，由于壓力和速度的值以及變化范圍差別很大，因此設定εU=εU0，εP=εP0。固定時間步長后，計算精度依賴于最大尺度因子和閾值，當尺度因子足夠大，則精度主要受閾值大小影響。數值試驗表明，同一時刻速度場的配點比壓力場配點多，并且速度場配點基本包含壓力場所有配點。如果尺度因子較小，可能會使速度場的計算精度達不到要求。

當取閾值ε=5×10-5(εU=εU0，εP=εP0)，尺度因子J=6，時間步長Δτ=2×10-7，計算無摩擦損失的水擊，可以得到較好的結果，如圖2(a)～圖2(d)所示。壓力穩定后，統計前3個周期(記為T1，T2和T3)的水擊波峰和波谷的平均值，如表2所示，水擊壓力平均值與理論值的誤差小于0.9%。時間和空間的大梯度處分辨率較高，且無明顯數值振蕩。由于時間積分采用的是只有一階精度的Euler格式，計算精度和穩定性對閾值、尺度因子以及時間步長要求嚴格，數值試驗表明，為了使計算穩定，參數應滿足：ε<0.001，J>3，Δτ=1×10-6。

圖2(a)表明，在第一個T/2內，壓力的數值振蕩較小，隨后數值振蕩相對明顯。由于閥前速度恒為0，因此圖2(b)表示的閥前速度計算值同時代表速度的計算誤差。閥前速度的計算誤差同時顯示，在第一個周期內計算誤差較小，約為1×10-8，在隨后的周期內最大誤差保持在4×10-5以內。

表2 波峰和波谷平均值

真實的流體管路有一定摩擦損失，造成水擊壓力逐漸下降。假定摩擦系數f=0.018，并取J=6，ε=5×10-5，Δτ=2×10-7，計算結果如圖3所示。

圖3 考慮摩擦(f=0.018)的水擊波傳播過程Fig.3 Progagation process of water-hammer with friction (f=0.018)

3.2 自適應算法和計算效率

自適應小波配點法能夠自動判斷梯度的大小和位置，自動在大梯度處增加配點，并刪除平穩區域內對計算精度影響較小的配點，在保證整個計算域內均達到設定精度的前提下減少配點數，從而減小計算量。因此在梯度較大的位置配點密集，而在變化平緩的位置配點稀疏。圖4給出了不同時刻壓力場配點分布和配點數(f=0.018，J=6，ε=5×10-5，Δτ=2×10-7)，為了和小波空間W0(下標為0)區分，尺度空間V0以-1表示。對比圖2(c)，可以看出配點加密的位置對應水擊波的波陣面，表明依據小波系數能夠準確判斷局部奇異性的位置。

圖4(c)表明，根據初始條件，壓力初值處處為1，數值隨空間無變化，因此僅需25個配點，水擊發展初期間斷較強，所需配點數較多，約為125，隨著水擊波的發展，梯度略有減小，不同時刻所需配點數在105附近，而不采用自適應算法所需均勻配點數為513，因此壓縮比約為4.89，節省了約80%的迭代計算量。

圖4 壓力場自適應配點Fig.4 Collocation points of pressure fields

4 結論

采用自適應小波配點法對一維管路水擊問題進行了數值計算，結果表明：

1)小波配點法能夠準確計算水擊波的傳播過程，在水擊波波陣面處，不會因高分辨率出現劇烈的數值振蕩現象，精度和分辨率可以同時提高。

2)自適應算法能夠自動捕捉水擊波位置并采用高分辨率計算，而在梯度較小區域降低分辨率，在保障精度和分辨率的前提下，顯著節約計算量。

本文的小波分解和自適應方法均針對空間，無法控制時間誤差。后續可采用時間-空間同步自適應方法，在局部變化劇烈的位置采用更精細的時間積分，從而控制時間積分過程產生的誤差。