張 亮，艾邦成，陳 智

(中國航天空氣動力技術研究院，北京 100074)

在計算流體力學領域，非結構有限體積方法由于其優(yōu)異的復雜外形適應性在各類武器型號設計中得到了廣泛的應用。在非結構數(shù)值模擬技術中，對流格式、限制器以及隱式迭代方法等均得到了廣泛的研究，但針對黏性項離散的耗散格式長期以來并未得到足夠的重視。近年來，隨著精細化數(shù)值模擬需求(氣動加熱、摩阻以及分離流高精度模擬)的不斷提升，非結構有限體積耗散格式逐漸成為重要的研究方向。

在有限體積方法框架內，構造耗散格式的核心在于確定界面處的流動梯度。由于黏性項天然的橢圓特性，耗散格式構造需要滿足離散極值原理以避免“奇偶失聯(lián)”及數(shù)值穩(wěn)定性問題[1-2]。對于結構網(wǎng)格，由于其規(guī)則的方向特性，一般可基于當?shù)赜嬎阕鴺讼抵苯硬捎弥行母袷接嬎憬缑嫣荻萚3]。而對于非結構網(wǎng)格，由于計算網(wǎng)格的無序性，中心型耗散格式難以應用，常用的處理方法是借鑒結構網(wǎng)格中的“薄層”簡化思想[4]，將界面梯度沿主方向和次方向進行分解，對不同方向的梯度采用不同方法進行計算[5-10]。在實際應用中，Edge-Normal格式[7]和Face-Tangent[8]格式得到了更為廣泛的采用。相比較而言，對于存在較大扭曲的計算網(wǎng)格，F(xiàn)ace-Tangent格式具有更好的數(shù)值穩(wěn)定性[11]。

與傳統(tǒng)耗散格式構造方法不同，Nishikawa基于一維雙曲系統(tǒng)提出了一類新型alpha-damping耗散格式[12]。該格式由相容項和阻尼項兩部分組成，其中相容項用以保證格式的數(shù)值相容性，阻尼項用于抑制高頻誤差的發(fā)展。Nishikawa還進一步證明了包括Edge-Normal格式和Face-Tangent格式在內的多種耗散格式均可等價或近似等價于具有不同高頻阻尼的alpha-damping格式。Jalali等針對大量非結構有限體積耗散格式開展了計算精度及數(shù)值穩(wěn)定性研究，結果表明采用優(yōu)化高頻阻尼系數(shù)的alpha-damping格式具有最優(yōu)的數(shù)值表現(xiàn)[13]。

Nishikawa提出的alpha-damping格式為不同耗散格式建立了一個統(tǒng)一框架，因此對于非結構有限體積耗散格式的構造具有重要意義。但其理論分析是基于均勻網(wǎng)格和固定中心梯度格式，非結構網(wǎng)格以及不同梯度格式條件下高頻阻尼對計算精度的影響并不明確。本文借鑒alpha-damping格式相容項和高頻阻尼項結合的構造形式，采用一維非均勻網(wǎng)格模擬非結構網(wǎng)格幾何特征，通過對模型方程的理論分析和數(shù)值試驗，研究高頻阻尼系數(shù)在不同梯度格式條件下對耗散格式精度的影響規(guī)律。

1 耗散格式理論分析

1.1 格式精度分析

考慮一維耗散模型方程：

ut=uxx

(1)

(2)

(3)

對于線性重構有：

(4)

不同于原始alpha-damping格式，為表征網(wǎng)格非均勻特性的影響，這里采用加權平均方法計算相容項部分，同時高頻阻尼系數(shù)也定義為當?shù)刂怠?/p>

將式(3)和式(4)代入式(2)并整理可得：

(5)

在實際計算中，需要采用梯度格式近似計算單元梯度，這里采用2階精度梯度格式：

(6)

將式(6)代入式(5)并在xj處進行泰勒展開可得修正方程：

(7)

其中：

c0=0

c1=0

需要強調的是，這里的分析結果僅嚴格適用于2階及高于2階精度的梯度格式。對于1階精度梯度格式，相容性分析結果則與梯度格式具體的截斷誤差形式相關。

1.2 梯度格式影響分析

Nishikawa指出，當α=4/3時，alpha-damping格式在均勻網(wǎng)格條件下可達到4階精度[12]。雖然該優(yōu)化系數(shù)由2階中心梯度格式推導得到，但其在后續(xù)的復雜多維問題及NS方程擴展中均得到了廣泛的應用[14-16]。為確定α=4/3取值的普適性，這里針對一般梯度格式給出高頻阻尼系數(shù)最優(yōu)值與梯度格式之間的定量關系。

考慮一般梯度格式：

在均勻網(wǎng)格條件下，式(7)各項的系數(shù)可做進一步簡化：

c0=0,c1=0,c2=1,c3=(1-α)a1Δx

c4=[(1-α)(a2+1/6)+α/12]Δx2

顯然，對于均勻網(wǎng)格，格式(3)自動滿足相容性條件。當α=1時，修正方程截斷誤差可至少保證2階精度。當α≠1時，1階精度梯度格式僅能保證1階精度截斷誤差，2階精度梯度格式至少可保證2階精度截斷誤差。為實現(xiàn)2階精度以上的截斷誤差，需要滿足c4=0，即：

(8)

顯然，可實現(xiàn)更高精度的高頻阻尼系數(shù)最優(yōu)值取決于梯度格式的2階截斷誤差系數(shù)。

對于2階精度中心梯度格式(6)，有a2=1/6，因此高頻阻尼系數(shù)最優(yōu)值為α=4/3。由于該格式同時滿足a3=0，因此格式最終可實現(xiàn)4階精度。

對于a2≠1/6的2階及高于2階精度的梯度格式，高頻阻尼系數(shù)最優(yōu)值需由式(8)確定。特別的，對于高于2階精度的梯度格式有a2=0，因此高頻阻尼系數(shù)最優(yōu)值為α=2。對于一般情況a2∈[0,+∞)，有α∈(1,2]，見圖1。

2 算例測試

2.1 一維算例

針對一維泊松方程uxx=S開展算例測試。選擇解析解u=sin(πx)，對應源項S=-π2u。計算域取為x∈[0,1]，函數(shù)分布曲線如圖2所示。

計算網(wǎng)格采用均勻網(wǎng)格、隨機網(wǎng)格和拉伸網(wǎng)格3種類型，每種網(wǎng)格類型的單元數(shù)分別取10，20，40，80。

采用3種梯度計算格式：

對于均勻網(wǎng)格，g1和g2等價，均為2階精度。對于非均勻網(wǎng)格，g1為1階精度，g2為2階精度。g3為解析解，可認為是無窮高階精度。

對于相容項的加權系數(shù)，大量數(shù)值試驗表明采用代數(shù)平均和基于距離的線性加權平均均可實現(xiàn)較高的計算精度，且兩者差異并不顯著[13-17]。結合前文關于非均勻網(wǎng)格的相容性分析，這里選擇代數(shù)平均進行相容項梯度加權。

2.1.1 均勻網(wǎng)格

均勻網(wǎng)格的單元尺寸分布特征如圖3所示。

圖4～圖5給出了采用不同梯度格式時，高頻阻尼系數(shù)對均勻網(wǎng)格數(shù)值精度的影響。由于g1和g2在均勻網(wǎng)格下等價，這里僅給出g1的計算結果。

對于g1，高頻阻尼系數(shù)α=4/3時，截斷誤差和離散誤差實現(xiàn)了4階精度，其他高頻阻尼系數(shù)僅能實現(xiàn)2階精度。而對于g3，高頻阻尼系數(shù)α=2時，截斷誤差和離散誤差實現(xiàn)了4階精度，其他高頻阻尼系數(shù)僅能實現(xiàn)2階精度。這一結果與本文的理論分析一致。

2.1.2 隨機網(wǎng)格

隨機網(wǎng)格單元i的網(wǎng)格尺度ΔxR,i=ZFi·ΔxU，其中ΔxU=δmax/N為等效網(wǎng)格尺度，N為網(wǎng)格單元數(shù)，δmax為計算域長度，ZFi為單元i處的隨機縮放因子，本文限制ZFi∈[0.8,1.2]。網(wǎng)格尺度分布特征如圖6所示。

圖7～圖9給出了采用不同梯度格式時，高頻阻尼系數(shù)對隨機網(wǎng)格數(shù)值精度的影響。

對于截斷誤差，當α=1時，3種梯度格式均近似實現(xiàn)了1階計算精度。其中，g1由于梯度計算精度低于2階，在較密網(wǎng)格時截斷誤差精度有所損失；g2和g3由于嚴格滿足2階梯度計算精度，截斷誤差精度至少達到了1階。而對于其他高頻阻尼系數(shù)，格式無法實現(xiàn)相容性。這一結果與前文的理論分析是一致的。

對于離散誤差，不同梯度格式以及不同高頻阻尼系數(shù)均實現(xiàn)了近似2階計算精度，但高頻阻尼系數(shù)對離散誤差絕對值的影響在不同梯度格式條件下表現(xiàn)不同。其中g1和g2在α=4/3時具有最小的離散誤差，而g3在α=2時具有最小的離散誤差。這說明，非均勻網(wǎng)格條件下產生最小離散誤差的高頻阻尼系數(shù)與式(8)確定的均勻網(wǎng)格最優(yōu)值一致。

2.1.3 拉伸網(wǎng)格

圖11～圖13給出了采用不同梯度格式時，高頻阻尼系數(shù)對拉伸網(wǎng)格數(shù)值精度的影響。

對于截斷誤差，不同梯度格式以及不同高頻阻尼系數(shù)均實現(xiàn)了格式的相容性。為解釋這一現(xiàn)象，可考慮式(7)中的相容性條件。在拉伸網(wǎng)格條件下，采用代數(shù)梯度加權的相容性條件可簡化為

顯然，當ΔxU趨近于0時，β趨近于1，因此c2趨近于1，格式相容性得到間接滿足。由于拉伸網(wǎng)格光滑性較好，截斷誤差的高階項可實現(xiàn)較好的相互抵消，因此截斷誤差精度優(yōu)于理論分析結果。

對于離散誤差，不同梯度格式以及不同高頻阻尼系數(shù)均實現(xiàn)了近似2階計算精度。與隨機網(wǎng)格相同，高頻阻尼系數(shù)對離散誤差絕對值的影響在不同梯度格式條件下表現(xiàn)不同，產生最小離散誤差的高頻阻尼系數(shù)與式(8)確定的均勻網(wǎng)格最優(yōu)值一致。

2.2 多維算例

采用3種梯度計算格式：g1為高斯梯度格式；g2為距離加權最小二乘梯度格式；g3為梯度解析解。相容項仍然采用代數(shù)平均方法進行加權。

1) 均勻四邊形網(wǎng)格

計算域為x∈[0,1]，y∈[0,1]，x和y方向的波數(shù)為nx=1，ny=1，網(wǎng)格尺度分別為1/10，1/20，1/40和1/80。典型計算網(wǎng)格及計算結果如圖14所示。

圖15～圖16給出了采用不同梯度格式時，高頻阻尼系數(shù)對數(shù)值精度的影響。由于在均勻四邊形網(wǎng)格條件下，g1與g2均等價于二階中心格式，因此這里僅給出g1的計算結果。

可以看到，與一維均勻網(wǎng)格計算結果相同，g1和g3分別在其對應的高頻阻尼系數(shù)最優(yōu)值時實現(xiàn)了4階計算精度，而其他高頻阻尼系數(shù)僅能實現(xiàn)2階計算精度。

2) 均勻三角形網(wǎng)格

計算域為x∈[0,1]，y∈[0,1]，x和y方向的波數(shù)為nx=1，ny=1，計算網(wǎng)格采用對角化技術生成，平均網(wǎng)格尺度分別為1/10，1/20，1/40和1/80。典型計算網(wǎng)格及計算結果如圖17所示。

圖18～圖20給出了采用不同梯度格式時，高頻阻尼系數(shù)對數(shù)值精度的影響。

對于截斷誤差，由于g1和g2均無法嚴格滿足2階梯度計算精度，因此截斷誤差無法實現(xiàn)相容性。尤其對于g1，截斷誤差與網(wǎng)格尺度出現(xiàn)了負相關，格式不相容性更為嚴重。g3由于嚴格滿足2階梯度計算精度，截斷誤差在α=1時達到了1階精度。

對于離散誤差，g1仍然無法滿足格式相容性。g2和g3則均實現(xiàn)了2階計算精度，且其產生最小離散誤差的高頻阻尼系數(shù)與式(8)確定的均勻網(wǎng)格最優(yōu)值一致。

3) 拉伸四邊形網(wǎng)格

計算域為x∈[0,1]，y∈[0,0.05]，x和y方向的波數(shù)為nx=1，ny=20，計算網(wǎng)格在x方向等距分布，y方向由2.1.3節(jié)的拉伸方法得到，壓縮因子CF取為100。x和y方向的網(wǎng)格單元數(shù)分別為10，20，40，80。典型計算網(wǎng)格及計算結果如圖21所示。

圖22～圖24給出了采用不同梯度格式時，高頻阻尼系數(shù)對數(shù)值精度的影響。

對于截斷誤差，與一維拉伸網(wǎng)格算例相同，由于網(wǎng)格的光滑性，不同梯度格式以及不同高頻阻尼系數(shù)均實現(xiàn)了2階計算精度。

對于離散誤差，與前文算例相同，不同梯度格式產生最小離散誤差的高頻阻尼系數(shù)仍然與式(8)確定的均勻網(wǎng)格最優(yōu)值一致。

3 結論

1) 對于非均勻網(wǎng)格，在采用2階精度梯度格式條件下，采用代數(shù)平均加權的相容項和α=1的高頻阻尼項是保證耗散格式嚴格相容的必要條件。對于光滑的非均勻拉伸網(wǎng)格，由于局部誤差的相互抵消，格式相容性有較大改善。采用低于2階精度的梯度格式會降低格式的相容性。

2) 均勻網(wǎng)格條件下高頻阻尼系數(shù)最優(yōu)值與梯度格式的2階截斷誤差系數(shù)相關，不同梯度格式應對應不同的高頻阻尼系數(shù)最優(yōu)值。在非均勻網(wǎng)格條件下，采用均勻網(wǎng)格最優(yōu)值的高頻阻尼系數(shù)可得到更小的離散誤差。

本文研究工作雖然基于一維理論模型，但研究結論對于多維問題仍具有一定指導意義。但由于耗散格式精度與梯度格式緊密相關，而多維梯度重構不同于一維梯度重構，其受重構算法及具體網(wǎng)格特征的影響巨大，因此后續(xù)研究重點將集中于多維情況下不同梯度構造格式的精度特性分析，并以此為基礎深入開展多維耗散格式精度及相容性研究。