特定條件下高階WENO 格式計算結(jié)果誤差

劉君,韓芳,魏雁昕

大連理工大學(xué) 航天航空學(xué)院,大連 116024

計算流體力學(xué)(CFD)應(yīng)用于工程項目時,和所有工業(yè)軟件一樣,希望它“又快又準(zhǔn)”?；仡機(jī)FD 的發(fā)展歷史,40年前需要花費(fèi)幾年時間才能完成的計算工作量現(xiàn)在僅需幾分鐘,CFD 依靠計算機(jī)技術(shù)在“快”方面取得了巨大進(jìn)展。如果把增加網(wǎng)格量對提高計算結(jié)果精度的貢獻(xiàn)去掉,相比日新月異的“快”而言,CFD 在“準(zhǔn)”方面的發(fā)展可以用“停滯不前”來形容。盡管高精度格式作為研究熱點已持續(xù)了20多年,目前商業(yè)軟件廣泛采用的還是20世紀(jì)80年代提出的二階精度格式。

近年來,作者在對激波裝配(Shock-fitting)法和激波捕捉(Shock-capture)法的對比研究中發(fā)現(xiàn),對捕捉得到的計算結(jié)果,有時候使用高階WENO 格式的計算誤差明顯比一階迎風(fēng)格式的計算誤差大。針對這一現(xiàn)象,本文通過定性分析及算例測試等手段進(jìn)行討論,根據(jù)類型的不同,可分為以下4個問題。

1 間斷問題

間斷問題是超/高超聲速可壓縮流中不可避免的問題,本文對間斷問題的討論主要是因激波或接觸間斷運(yùn)動而引起的非物理波動問題。

以笛卡爾直角坐標(biāo)系下的二維守恒型Euler方程作為控制方程,在均勻正交網(wǎng)格上進(jìn)行離散求解。通量分裂采用van Leer格式,空間離散采用一階迎風(fēng)格式及五階WENO格式,包括WENO-JS、WENO-M、WENO-Z等。對應(yīng)的,時間離散采用一階顯示格式或具有TVD 性質(zhì)的三階Runge-Kutta格式,時間推進(jìn)步長按CFL=0.5計算。

1.1 激 波

計算區(qū)域設(shè)為[0,1]×[0,1],網(wǎng)格尺度Δ=Δ=0.01,選取初始位于=01處的正激波作為初始計算條件,激波運(yùn)動馬赫數(shù)為=30,波前波后參數(shù)為

選取無量綱時間=025時的計算結(jié)果進(jìn)行討論,此時激波運(yùn)動至=085的位置。圖1給出了流場密度、壓力及熵沿=05的分布曲線,由于3種WENO 格式的分布曲線幾乎一致,因此本文只給出了WENO-Z 格式的計算結(jié)果。從圖中可以看出,在激波波后,存在一等熵一非等熵兩個非物理波動,且與一階迎風(fēng)格式的結(jié)果相比,WENO-Z 格式捕捉到的波動寬度更窄,幅值更高。

圖1 y=0.5線上不同格式的流場參數(shù)分布曲線Fig.1 Flow field parameter curves along line y=0.5

這種波動現(xiàn)象在文獻(xiàn)[9]中已經(jīng)提及,本文除了把不同精度格式的結(jié)果放在一起比較外,還給出了激波數(shù)值過渡區(qū)的內(nèi)部特征結(jié)構(gòu)。

在激波坐標(biāo)系下,激波波前波后的參數(shù)應(yīng)符合R-H 關(guān)系式,由式(1)所決定的激波質(zhì)量通量、動量通量和能量通量分別為(,,)=(42,136,294)。圖2給出了=05線上WENO-Z格式和一階迎風(fēng)格式所模擬激波守恒通量的分布曲線,從圖中可以看出,盡管計算是以守恒型Euler方程作為控制方程,但在捕捉到的數(shù)值激波的過渡區(qū)內(nèi),其物理量的守恒特性不能保持,且與流場參數(shù)分布曲線相似,WENO-Z格式計算的守恒通量的局部最大誤差要大于一階迎風(fēng)格式的計算結(jié)果。

圖2 y=0.5線上不同格式的守恒通量分布曲線Fig.2 Conserved flux curves in different schemes along line y=0.5

提取流場內(nèi)沿方向流動參數(shù),對相鄰網(wǎng)格點進(jìn)行Reimann分解,圖3 是WENO-Z 格式捕捉到的激波過渡區(qū)在=02和=025兩個時刻的內(nèi)部結(jié)構(gòu)分布。操作中發(fā)現(xiàn)即使相鄰點有非常微小的變化也會分解出Reimann結(jié)構(gòu),因此設(shè)置壓差Δ=001作為閾值,小于該值認(rèn)為是光滑流場區(qū)域,用線段“”表示。

圖3 不同時刻激波過渡區(qū)內(nèi)的Reimann結(jié)構(gòu)分布Fig.3 Distribution of Reimann structures within excitation transition zone at different time points

從圖3可以看出,在激波的數(shù)值過渡區(qū)內(nèi)主要存在兩種Reimann分解結(jié)構(gòu),圖中“”線段表示得到的是“激波-激波”結(jié)構(gòu),“”線段表示得到的是“稀疏波-激波”結(jié)構(gòu),而在不同時刻,激波過渡區(qū)內(nèi)的Reimann分解結(jié)構(gòu)分布并不相同。此外,作者還考察了不同格式和不同運(yùn)動馬赫數(shù)的激波過渡區(qū),發(fā)現(xiàn)在過渡區(qū)內(nèi)兩種結(jié)構(gòu)之間的變化復(fù)雜,其分布不具有規(guī)律性。

保持計算網(wǎng)格不變,將以上正激波旋轉(zhuǎn)一定角度,激波與軸的夾角記為,初始激波放置在網(wǎng)格對角線下游以保證超聲速出口邊界不影響內(nèi)部流場,如圖4所示。文獻(xiàn)[9]給出了不同旋轉(zhuǎn)角度下激波運(yùn)動所產(chǎn)生的密度誤差云圖和渦量云圖,認(rèn)為在網(wǎng)格點與初始激波不匹配的條件下,高精度格式的高分辨率特性很容易誘導(dǎo)出類似“湍流”的虛假結(jié)構(gòu)。圖5給出了=43.5°時一階迎風(fēng)格式及五階WENO-Z格式所計算的流場渦量云圖。圖6是=01時刻,從(0.4,0)點出發(fā)與初始激波平行的直線上的渦量值分布曲線。由圖5的渦量分布云圖及圖6的渦量值分布曲線可以看出,WENO-Z格式得到的激波波后流場結(jié)構(gòu)分布更復(fù)雜,流場參數(shù)誤差的最大幅值也更大。

圖4 計算模型示意圖Fig.4 Schematic diagram of calculation model

圖5 不同精度格式得到的渦量云圖分布Fig.5 Vorticity contours obtained by different precision formats

圖6 不同精度格式得到的渦量分布曲線Fig.6 Vorticity distribution curves obtained by different precision formats

在文獻(xiàn)[9]中,作者通過理論分析及數(shù)值實驗得出以下結(jié)論：激波在計算過程中由數(shù)學(xué)間斷發(fā)展成數(shù)值過渡區(qū)的過程中必然會出現(xiàn)以上非物理現(xiàn)象。結(jié)合本文結(jié)果可進(jìn)一步發(fā)現(xiàn),高精度格式的高分辨率特性會放大這種非物理的虛假波動,導(dǎo)致波后流場局部出現(xiàn)比一階迎風(fēng)格式結(jié)果更大的數(shù)值誤差。

1.2 接觸間斷

接觸間斷算例的計算區(qū)域及初始條件如圖7所示,網(wǎng)格尺度Δ=Δ=0.01。在此算例中,流場全場壓力相等、速度相等,因此渦量=05×(v -u )的理論值為0。由于流場的壓力及渦量分布在接觸間斷上下兩側(cè)非常相似,因此分別取接觸間斷上側(cè)WENO-Z格式的結(jié)果和下側(cè)一階迎風(fēng)格式的結(jié)果進(jìn)行定性比較,如圖8所示。由于WENO-Z格式的結(jié)果存在周期波動,圖中取一階迎風(fēng)格式滿足收斂條件2倍時間的計算結(jié)果進(jìn)行比較。同時為了考察流場參數(shù)的波動特性,在流場下游設(shè)置了3個觀測點(,)=(1.5,-0.5),(1.5,0),(1.5,0.5),觀測點上壓力隨時間的變化曲線如圖9所示。

圖7 二維接觸間斷：計算區(qū)域和初始條件Fig.7 2D contact discontinuity flow：computational domain and initial conditions

圖8 不同精度格式的壓力和渦量云圖Fig.8 Pressure and vorticity contours obtained by different precision formats

圖9 接觸間斷測點處不同格式壓力值隨時間變化曲線Fig.9 Pressure change curves with time in different schemes at measurement points

在超聲速流場求解過程中,接觸間斷從初始數(shù)學(xué)上的間斷變成有厚度的數(shù)值剪切層,必然會誘導(dǎo)出旋渦、激波、膨脹波及其相互干擾的非物理結(jié)構(gòu)。從本文算例看,不同精度格式的數(shù)值誤差量級接近,但是,高精度格式出現(xiàn)了和物理流動非常相似的剪切層失穩(wěn)現(xiàn)象。目前高精度格式常用來模擬轉(zhuǎn)捩、DNS等非定常問題,如何從數(shù)值解中剔除這種虛假波動值得關(guān)注。

2 坐標(biāo)變換問題

從笛卡爾直角坐標(biāo)系(,)變換到柱坐標(biāo)系(,)的變換函數(shù)可寫成如下形式：

為了避免圓心處的奇性,取半徑∈[1,2]的圓環(huán)作為計算區(qū)域,網(wǎng)格沿徑向和周向均勻分布,這是貼體坐標(biāo)系中嚴(yán)格正交、充分光滑的理想網(wǎng)格。采用與間斷問題相同的計算格式數(shù)值模擬馬赫數(shù)=3的均勻流,無量綱流動參數(shù)為

在文獻(xiàn)[10]的計算中,作者發(fā)現(xiàn)在＞0 區(qū)域內(nèi)環(huán)進(jìn)口邊界的處理算法會影響流場參數(shù)誤差的分布特性,因此,本文直接選擇＜0的區(qū)域計算,計算網(wǎng)格如圖10所示。圖11(a)是由一階迎風(fēng)格式計算得到的流場壓力云圖分布,圖11(b)是由五階WENO-Z格式計算得到的流場壓力云圖分布。圖12是兩種精度格式所得壓力沿=15的周向分布曲線,其中橫坐標(biāo)為角度,以逆時針旋轉(zhuǎn)為正。

圖10 均勻流場計算網(wǎng)格Fig.10 Computational grid for uniform flow

圖11 流場壓力云圖Fig.11 Contour of flow field pressure

圖12 不同格式得到的沿r=1.5的周向壓力分布曲線Fig.12 Pressure distribution curves along line r=1.5 for different schemes

在貼體變標(biāo)系下應(yīng)用WENO 格式計算均勻流場出現(xiàn)誘導(dǎo)誤差的現(xiàn)象很早就被注意到,CFD工作者們對其消除方法進(jìn)行了研究,本文通過定量比較,為“高精度格式的幾何誘導(dǎo)誤差更大”這一論點提供了論據(jù)。下面根據(jù)WENO 格式構(gòu)造過程對壓力沿周向變化這個新現(xiàn)象進(jìn)行定性解釋。

盡管計算空間網(wǎng)格總是均勻的,但是構(gòu)造模板函數(shù)時不能直接采用直角坐標(biāo)系下的原始變量(,,,)和守恒變量=(,,,),因 為這些量沒有體現(xiàn)坐標(biāo)變換過程中網(wǎng)格距離和網(wǎng)格面積等幾何參數(shù)的變化。以二階精度中心格式為例簡單說明,在非均勻網(wǎng)格Δ≠Δ條件下,速度導(dǎo)數(shù)的離散形式為

如果直接套用均勻網(wǎng)格的形式：

近年來國內(nèi)外研究幾何守恒律的學(xué)者認(rèn)為,在有限差分法中坐標(biāo)變換誘導(dǎo)誤差出現(xiàn)的原因是采用軟件生成的網(wǎng)格不能保證其正交性和光滑性,坐標(biāo)變換系數(shù)的算法精度不能匹配有限差分格式的高精度,上述算例不支持這種觀點。本文作者對其產(chǎn)生機(jī)理進(jìn)行理論分析并且提出新的消除算法。

3 邊界條件問題

CFD 在數(shù)學(xué)領(lǐng)域又稱為偏微分方程的初邊值問題,其中對于定常問題,初值僅影響收斂過程,由邊界條件決定流場特性,因此理論上應(yīng)該是將邊界條件作為約束代入到流體控制方程建立邊界方程。但是,由于邊界方程與內(nèi)部方程形式不同,再加上基于內(nèi)點構(gòu)造的差分格式很難直接用于邊界點計算,因此,在實際應(yīng)用中常采用簡化模型。例如,無黏流動在固壁處滿足不穿透條件,很少按照早期嚴(yán)謹(jǐn)?shù)淖龇?把V =0代入Euler方程的特征線形式推導(dǎo)出邊界方程后離散求解,目前常規(guī)的處理是在假設(shè)?/?=0的基礎(chǔ)上,將其他物面參數(shù)根據(jù)等熵假設(shè)或者直接再進(jìn)一步補(bǔ)充條件?/?=0和?V/?=0,然后離散以上空間導(dǎo)數(shù),形成邊界點和內(nèi)點的關(guān)系式,在時間推進(jìn)過程中根據(jù)內(nèi)點信息預(yù)測更新邊界值。

邊界問題的算例計算模型如圖13所示,模擬全部區(qū)域是激波對稱相交流動,僅模擬上半部分就是激波在固體壁面的正規(guī)反射。按激波在固體壁面做正規(guī)反射求解時,由于固體壁面是平面,流場上下對稱,可以通過在固體內(nèi)部設(shè)置虛擬網(wǎng)格點直接求解邊界值,計算得到的壁面參數(shù)和全流場中間對稱流線上的流場參數(shù)完全相等。對于存在曲率的固體壁面,對稱原理明顯不成立,物體內(nèi)部虛擬點如何取值是個問題,因此這種“鏡像法”不具有普適性。在物體內(nèi)部沒有網(wǎng)格點的情況下,面臨如何離散?/?=0等空間導(dǎo)數(shù)的問題,應(yīng)用中發(fā)現(xiàn)單邊高階格式離散容易引起穩(wěn)定性問題,因此,很多時候盡管內(nèi)部網(wǎng)格點采用高精度格式,邊界網(wǎng)格點也是采用一階外推得到。但是在圖13所示算例中,在中心對稱流線處?/?=0成立,可以用來評價固體壁面的邊界算法。

圖13 激波正規(guī)反射計算模型示意圖Fig.13 Sketch map for shock regular reflection

在全流場計算中,中間對稱流線的壓力記為,上下相鄰網(wǎng)格線上的流場壓力分別記為和,壓力差Δ=-是數(shù)值模擬計算結(jié)果。如果模擬固體壁面正規(guī)反射,一階外推得到邊界值表示為：=,這種簡化處理等效于在邊界處引入Δ誤差。

取圖13虛線框內(nèi)區(qū)域為計算區(qū)域,使用笛卡爾直角坐標(biāo)系下的均勻正交網(wǎng)格,網(wǎng)格初始尺度為Δ=Δ=0.01。同時,為考察網(wǎng)格尺度對計算結(jié)果的影響,對網(wǎng)格進(jìn)行2次加密,加密方式為一次加密即將初始網(wǎng)格或上一次加密網(wǎng)格的一個單元均分為4個網(wǎng)格單元,將初始網(wǎng)格及2次加密網(wǎng)格分別記為網(wǎng)格1、網(wǎng)格2、網(wǎng)格3。一階迎風(fēng)格式和WENO-JS 格式的Δ誤差分布曲線見圖14,可以看出WENO-JS格式得到的誤差曲線幅值更大。圖15給出了網(wǎng)格3上Δ=-的分布曲線,發(fā)現(xiàn)一階迎風(fēng)格式得到的結(jié)果為0,WENO-JS格式得到的曲線顯示流場參數(shù)在中心對稱線上下有明顯的不對稱現(xiàn)象。

圖14 不同格式得到的邊界誤差分布Fig.14 Boundary error distributions for different schemes

圖15 不同格式中心對稱線上下壓力差分布曲線Fig.15 Pressure difference distributions curves for different schemes above and below center line

采用以上設(shè)計模型,通過全流場Δ評估格式在邊界降階算法的誤差,結(jié)論沒有普適性,但至少表明高精度格式存在放大邊界處算法誤差的風(fēng)險。

4 構(gòu)造高精度格式的邏輯問題

從2013 年開始,NASA 聯(lián)合航空航天領(lǐng)域眾多高校、研究機(jī)構(gòu)和企業(yè)的知名專家,面向美國已經(jīng)開展及正在規(guī)劃的重大工程項目需求,歷時兩年形成了一份預(yù)測未來CFD 軟件發(fā)展前景的咨詢報告。報告認(rèn)為CFD 技術(shù)的發(fā)展目前已進(jìn)入瓶頸階段,并列舉了CFD 要得到進(jìn)一步發(fā)展需要攻克的六大關(guān)鍵技術(shù),其中在“數(shù)值算法(Numerical Algorithms)”一項中提到了應(yīng)用高精度格式遇到的問題。在以上具體算例的基礎(chǔ)上,本文作者提出如下問題與同行交流討論。

4.1 守恒型方程

通過前面的算例和定性分析可以發(fā)現(xiàn),從守恒型方程出發(fā)的捕捉法,得到的激波數(shù)值過渡區(qū)內(nèi)必然會出現(xiàn)質(zhì)量通量、動量通量和能量通量不守恒的現(xiàn)象。既然“守恒”方程捕捉到的激波過渡區(qū)并不守恒,是否可以突破守恒方程的限制,采用新的方法數(shù)值模擬激波呢?

受肖鋒等構(gòu)造THINC 格式時用雙曲正切函數(shù)描述流場接觸間斷的啟發(fā),根據(jù)激波前后參數(shù)在5個網(wǎng)格點之間擬合人工過渡區(qū)作為激波初始條件,顯然內(nèi)部3點不滿足守恒性,但圖16所示的計算結(jié)果表明,以人工過渡區(qū)作為初始條件的正激波在運(yùn)動中不會產(chǎn)生1.1節(jié)中的非物理波動,采用函數(shù)模型來模擬激波也是可行的。

圖16 不同初始條件下的正激波密度曲線(t=0.25)Fig.16 Density curves under different initial conditions for positive shock waves(t=0.25)

4.2 通量分裂格式

根據(jù)前面的數(shù)值模擬結(jié)果,如果掌握了激波過渡區(qū)內(nèi)參數(shù)分布的規(guī)律性,就可以構(gòu)建出沒有虛假波動的初始激波模型。由于理論建模較為困難,首先采用精確Riemann解對數(shù)值過渡區(qū)內(nèi)部波系結(jié)構(gòu)進(jìn)行分析,包括不同的通量分裂格式和差分格式,至今未發(fā)現(xiàn)激波過渡區(qū)內(nèi)波系結(jié)構(gòu)隨激波參數(shù)的變化規(guī)律。

通量分裂格式既不能避免間斷非物理,又不能準(zhǔn)確反映雙曲型方程守恒量沿特征線傳播的空間特性,那么研究高精度格式是否固守這個技術(shù)途徑值得探討。

4.3 多點格式時空特性

考慮一維Euler方程,擾動的傳播如圖17所示。對一初始為靜止的流場,在+2點處給予一擾動,擾動以聲速傳播。對于本文所用的顯式格式,在滿足穩(wěn)定性要求的Δ時間內(nèi),向上游傳播的影響范圍不會超過+1點,因此,不可能影響到過點的流場參數(shù)。其次,即使不考慮時間因素,+2點的擾動能影響到點的流場,但是擾動經(jīng)過+1點時要發(fā)生波的相互干擾,根據(jù)1.1節(jié)的Riemann分解結(jié)果,流場結(jié)構(gòu)復(fù)雜,相互干擾效應(yīng)不是高精度格式的有限模板及其固定系數(shù)能夠描述的。

圖17 擾動傳播示意圖Fig.17 Schematic diagram of disturbance propagation

文獻(xiàn)[19]采用大的時間步長來關(guān)聯(lián)空間多點的研究思路值得進(jìn)一步探索。

4.4 特征線和依賴域

圖18表示的是一個脈沖小擾動在超聲速均勻流場中的傳播,隨著速度向下游傳播過程中影響范圍以聲速的倍數(shù)逐漸擴(kuò)大,流體微團(tuán)到達(dá)點時至少點以上的區(qū)域已經(jīng)感受不到曾經(jīng)的擾動了。如果是持續(xù)擾動,在達(dá)到定常以后,點物理量的變化沿特征線保持,換言之,馬赫錐內(nèi)流場不再感受其擾動。由此看來,擾動的影響域在定常和非定常兩種條件下是有差異的。同樣,對于依賴域也存在類似的問題。圖19是在超聲速流場設(shè)計中得到廣泛應(yīng)用的二維特征線法,對于定常二維超聲速等熵流動,只要知道2條特征線和流線的上游參數(shù)就可以確定點的流場,如果整個依賴域內(nèi)還有其他點對其有影響,則變成了超定問題,與理論存在矛盾。此外,在構(gòu)造點處的高階格式時,必然會用到點在方向的相鄰多點作為構(gòu)造模板,但根據(jù)圖19 可以看出,這些點對點處的流場并無影響。

圖18 小擾動在超聲速流場中的傳播[20]Fig.18 Propagation of weak perturbation in supersonic flow fields[20]

圖19 定常二維超聲速流場中的特征線與流線[21]Fig.19 Characteristic lines and streamlines in constant two-dimensional supersonic flows[21]

由此看來,目前采用多點擬合方法構(gòu)建高精度格式的思路和雙曲型方程描述波的傳播過程在理論上存在矛盾之處。

5 結(jié) 論

有限差分的高階WENO 格式因其強(qiáng)大的捕捉能力及辨識能力,使其能夠捕捉到流場中的細(xì)微結(jié)構(gòu),對提高CFD 對流場的數(shù)值模擬精度具有重大意義。但以上算例表明,在某些特定情況下,高階WENO 格式存在放大流場計算結(jié)果誤差的風(fēng)險,在使用時需要特別注意。通過對以上算例結(jié)果的討論和分析,得到如下結(jié)論：

1)通過對數(shù)值激波附近守恒通量的比較,可以看出高階WENO 格式得到的通量數(shù)值誤差大于一階迎風(fēng)格式,且高階格式捕捉到的流場參數(shù)在激波和接觸間斷處更容易出現(xiàn)非物理振蕩。

2)如果沒有配套的幾何守恒律,在同一曲線網(wǎng)格下,高階WENO 格式得到的幾何誘導(dǎo)誤差大于低階格式。

3)高階WENO 格式存在放大邊界處算法誤差的風(fēng)險。

4)高階WENO 格式采用多點模板的思路本質(zhì)上不符合雙曲型方程的特性。

綜上,或許正如文獻(xiàn)[22]所說：也許到了突破守恒方程、通量分裂、多點格式這些傳統(tǒng)理論的時候了。