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

2022-03-29 07:49:24劉君韓芳魏雁昕

航空學報 2022年2期

劉君,韓芳,魏雁昕

大連理工大學航天航空學院,大連 116024

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

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

1 間斷問題

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

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

1.1 激波

計算區域設為[0,1]×[0,1],網格尺度Δ=Δ=0.01,選取初始位于=01處的正激波作為初始計算條件,激波運動馬赫數為=30,波前波后參數為

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

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

這種波動現象在文獻[9]中已經提及,本文除了把不同精度格式的結果放在一起比較外,還給出了激波數值過渡區的內部特征結構。

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

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

提取流場內沿方向流動參數,對相鄰網格點進行Reimann分解,圖3 是WENO-Z 格式捕捉到的激波過渡區在=02和=025兩個時刻的內部結構分布。操作中發現即使相鄰點有非常微小的變化也會分解出Reimann結構,因此設置壓差Δ=001作為閾值,小于該值認為是光滑流場區域,用線段“”表示。

圖3 不同時刻激波過渡區內的Reimann結構分布Fig.3 Distribution of Reimann structures within excitation transition zone at different time points

從圖3可以看出,在激波的數值過渡區內主要存在兩種Reimann分解結構,圖中“”線段表示得到的是“激波-激波”結構,“”線段表示得到的是“稀疏波-激波”結構,而在不同時刻,激波過渡區內的Reimann分解結構分布并不相同。此外,作者還考察了不同格式和不同運動馬赫數的激波過渡區,發現在過渡區內兩種結構之間的變化復雜,其分布不具有規律性。

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

圖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

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

1.2 接觸間斷

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

圖7 二維接觸間斷：計算區域和初始條件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

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

2 坐標變換問題

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

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

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

圖10 均勻流場計算網格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

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

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

如果直接套用均勻網格的形式：

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

3 邊界條件問題

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

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

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

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

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

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

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

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

4 構造高精度格式的邏輯問題

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

4.1 守恒型方程

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

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

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

4.2 通量分裂格式

根據前面的數值模擬結果,如果掌握了激波過渡區內參數分布的規律性,就可以構建出沒有虛假波動的初始激波模型。由于理論建模較為困難,首先采用精確Riemann解對數值過渡區內部波系結構進行分析,包括不同的通量分裂格式和差分格式,至今未發現激波過渡區內波系結構隨激波參數的變化規律。

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

4.3 多點格式時空特性

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

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

文獻[19]采用大的時間步長來關聯空間多點的研究思路值得進一步探索。

4.4 特征線和依賴域

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