張 亮，艾邦成，陳 智

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

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

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

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

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

1 耗散格式理論分析

1.1 格式精度分析

考慮一維耗散模型方程：

ut=uxx

(1)

(2)

(3)

對于線性重構(gòu)有：

(4)

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

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

(5)

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

(6)

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

(7)

其中：

c0=0

c1=0

需要強調(diào)的是，這里的分析結(jié)果僅嚴(yán)格適用于2階及高于2階精度的梯度格式。對于1階精度梯度格式，相容性分析結(jié)果則與梯度格式具體的截斷誤差形式相關(guān)。

1.2 梯度格式影響分析

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

考慮一般梯度格式：

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

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

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

顯然，對于均勻網(wǎng)格，格式(3)自動滿足相容性條件。當(dāng)α=1時，修正方程截斷誤差可至少保證2階精度。當(dāng)α≠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)，對應(yīng)源項S=-π2u。計算域取為x∈[0,1]，函數(shù)分布曲線如圖2所示。

計算網(wǎng)格采用均勻網(wǎng)格、隨機(jī)網(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為解析解，可認(rèn)為是無窮高階精度。

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

2.1.1 均勻網(wǎng)格

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

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

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

2.1.2 隨機(jī)網(wǎng)格

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

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

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

對于離散誤差，不同梯度格式以及不同高頻阻尼系數(shù)均實現(xiàn)了近似2階計算精度，但高頻阻尼系數(shù)對離散誤差絕對值的影響在不同梯度格式條件下表現(xiàn)不同。其中g(shù)1和g2在α=4/3時具有最小的離散誤差，而g3在α=2時具有最小的離散誤差。這說明，非均勻網(wǎng)格條件下產(chǎn)生最小離散誤差的高頻阻尼系數(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ù)梯度加權(quán)的相容性條件可簡化為

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

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

2.2 多維算例

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

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)格及計算結(jié)果如圖14所示。

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

可以看到，與一維均勻網(wǎng)格計算結(jié)果相同，g1和g3分別在其對應(yīng)的高頻阻尼系數(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)格采用對角化技術(shù)生成，平均網(wǎng)格尺度分別為1/10，1/20，1/40和1/80。典型計算網(wǎng)格及計算結(jié)果如圖17所示。

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

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

對于離散誤差，g1仍然無法滿足格式相容性。g2和g3則均實現(xiàn)了2階計算精度，且其產(chǎn)生最小離散誤差的高頻阻尼系數(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)格及計算結(jié)果如圖21所示。

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

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

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

3 結(jié)論

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

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

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