應用維數分裂方法推廣MUSCL和WENO格式的若干問題

劉君，韓芳，魏雁昕

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

2019年，在國家數值風洞工程基礎研究課題的支持下，作者與國內研究高超聲速邊界層轉捩問題的專家學者合作開展激波裝配法在邊界層轉捩感受性問題中的應用研究，期望利用裝配法消除采用捕捉法數值模擬感受性問題時在激波附近出現的虛假數值噪聲，以發展全場一致的高精度數值算法。

為分析數值噪聲產生的機理，選擇應用廣泛的五階WENO格式作為高精度格式在有限差分框架下進行數值實驗。在對數值實驗結果進行分析時，發現在流場中除了激波過渡區容易出現非物理振蕩外，高精度格式也可能放大因坐標變換而引起的幾何誘導誤差，至少在給出的曲線坐標系下的均勻流算例中，五階WENO格式模擬的流場參數誤差大于一階迎風格式的數值模擬結果。文獻調研發現，在曲線坐標系下的WENO格式確實存在會產生幾何誘導誤差的問題，需要構造與之相對應的消除算法，文獻[8-9]的算例也驗證了這一觀點。幾何誘導誤差問題屬于幾何守恒律(Geometric Conservation Law，GCL)研究領域，有幾何誘導誤差產生即證明算法不滿足幾何守恒律。但是，本文作者于2020年1月參加“第二屆計算流體力學中高精度方法及其應用研討會”提出這一觀點時卻受到了質疑，“定性分析，構造WENO格式時使用的是積分運算，模板作用于控制體內物理量的均值，計算又是從守恒型Euler方程出發，怎么會不守恒？”

根據針對對象的不同，幾何守恒律又可分為體積守恒律和面積守恒律，體積守恒律對應于運動及變形網格，面積守恒律對應于靜止網格，本文所討論問題皆基于靜止網格。在有限體積法(Finite Volume Method，FVM)中，與網格相關的幾何參數可以精確地求得，幾何守恒律自動滿足，而有限差分法(Finite Difference Method，FDM)則存在幾何守恒律不能得到滿足的情況。因此，關于WENO格式守恒性的討論也可以變成對WENO格式是否屬于有限體積法的討論。“WENO格式是否屬于有限體積法”是本文提出的第1個問題。

由于國內外計算流體力學(CFD)相關文獻及教材對FVM的定義不同，因此本文討論的特指對控制體采用高斯公式將體積分轉換成面積通量積分的有限體積法，可稱為“高斯積分型有限體積法”，為表述方便，以G-FVM表示。早期國內教材認為G-FVM和FDM存在明顯差異，“對于較均勻的網格，如果坐標變換函數的計算能具有高精度，那么有限差分法可以通過多點格式或緊致格式獲得高精度，但對于有限體積法很難獲得二階以上精度”。后來的教材采用了國外文獻觀點，認為“有限差分法和有限體積法的不同主要是對網格的幾何處理方法不同，二者沒有本質區別”；“在一般曲線坐標系中，對幾何度量系數做適當處理后，守恒型流動方程組有限差分算法和物理空間里有限體積算法是等價的”。綜合以上觀點，在回答WENO格式是否屬于有限體積法之前，需要先弄清楚第2個問題：G-FVM和FDM是否等價？

在數值實驗及結果分析中發現從一維Euler方程出發介紹WENO格式時，如果模板作用于對流項通量，即便在非均勻網格上也不存在GCL問題。但是對于多維空間，即使是均勻正方形網格也會出現與幾何參數相關的均值計算問題。由于WENO格式模板函數比較復雜，推導過程選擇形式較為簡單的MUSCL格式。

最早構造MUSCL格式時采用“Specific Volume”來體現網格的非均勻特性，但是國內CFD教材在介紹MUSCL格式時大部分采用直角坐標系下的均勻網格，這種情況下幾何度量沒有變化，“限制器的變量可以是守恒變量、原始變量、通量、特征變量等，但限制器的函數都是相同的”。即使空間一維，在非均勻網格上插值也要考慮距離的影響，定性分析空間多維曲線坐標系下的方程表達式，守恒變量和對流項通量中的幾何度量系數存在量綱上的差異，由此引出第3個問題：如何使MUSCL格式和WENO格式守恒？

1 G-FVM和FDM是否等價

首先給出FDM和G-FVM的控制方程，FDM的控制方程在笛卡爾坐標系(,,,)下得到，以三維守恒型Euler方程為例：

(1)

式中：為守恒變量；、、分別為、、方向的對流通量，其表達式為

=[,,,,]

=[,+,,,(+)]

=[,,+,,(+)]

=[,,,+,(+)]

其中：為密度；、、分別為、、方向的速度；為壓力；為單位質量流體的總能，其計算式為

在計算具有復雜邊界或外形的流場時，需要將控制方程由物理空間(,,,)變換到計算空間(,,,)，變換后的守恒型控制方程為

(2)

當網格靜止時，

=

式中：為Jacobian行列式，有

、、、、、、、、為度量系數，可通過以下公式求得：

與FDM不同，G-FVM可以在曲線坐標系下直接應用，采用高斯公式把體積分變成面積分：

(3)

式中：是控制體體積；是控制體外表面；d表示微元體積；d表示外表面微元面積；=++表示外表面微元外法向單位矢量。對流項通量積分的函數可寫為

()=++

(4)

為便于幾何表達和簡潔符號，采用圖1所示二維空間模型，其可看做是三維網格在方向的等值面，灰色曲線表示建立FDM的網格點，黑色曲線構成G-FVM的控制體。

文獻[12-13]討論了G-FVM和FDM的不同，這里引用其主要觀點：

1) FDM的數值解是網格節點處的物理量，可以采用Taylor級數寫出修正方程，從而直接得到精度、相容性等格式特性；G-FVM得到的是控制體的平均值，分析格式特性需要采用其他數學手段。

2) 即使兩種方法采用數值上相同的幾何度量系數，其精度也有差異，控制體和積分面的幾何參數對于G-FVM是精確的，用在FDM中必然是近似，根據物理空間(,,,)變換到計算空間(,,,)上的函數解析表達式計算出來的幾何度量系數不能直接用于G-FVM。

3) G-FVM保持守恒性。對于均勻流動參數，采用各種精度的格式和任意形狀的網格，圍成封閉控制體所有面的通量積分之和為0，自動滿足GCL。

除此之外，本文進一步補充如下論據：

4) FDM不需要封閉的控制體或積分面(對封閉性的理解見附錄A)。在實際CFD數值模擬中，一般需要先通過CAD實體模型生成物面網格，進而形成流場計算所需的空間網格。對G-FVM來說，要求控制體的每個積分面滿足封閉性，否則無法進行積分計算。例如，對圖2所示物體表面存在細縫的情況，便需要通過幾何修形將細縫“抹平”。而課題組近期研究工作表明FDM本質上不需要表面網格，只要找到物面點的位置信息就可以進行計算，不需要關心不同物面網格點之間的關系。

對于二階精度的G-FVM，控制體梯度是利用相鄰網格點參數重構出來的，不能保證為0，難以滿足絕熱壁條件。根據格心值和格心到壁面矢量重構得到的壁面溫度分布不是常數，最多只能在一個點滿足等溫壁條件。

6) 高斯公式只能用到空間二維，把面積分轉換成線積分進行計算，在空間一維中無法使用，因此，G-FVM不能退化到空間一維。

通過以上論述可知：G-FVM和FDM存在本質差異。

圖1 二維曲線網格Fig.1 Two-dimensional curvilinear grid

圖2 物面附近網格示意圖Fig.2 Schematic diagram of grid near surface

2 MUSCL和WENO格式是否屬于G-FVM

按照G-FVM和FDM的特點，本文認為MUSCL類和WENO格式不屬于G-FVM，下面通過G-FVM界面通量計算過程提供更多的論據。

以圖3所示直角網格為研究對象，其中陰影部分為控制體，其流場參數的平均值為

(5)

圖3 均勻直角網格Fig.3 Uniform rectangular grid

應用高斯公式后通量積分變成圍繞控制體的線積分，以圖3中紅色線段-12,為例，法向矢量=-，對流項通量為

(6)

目前看到MUSCL和WENO格式的構造過程中只有線積分，從一維Euler方程出發，在半點界面處通量的通用形式為

(7)

如果采用MUSCL和WENO格式計算式(6) 中G-FVM的界面通量，有兩種途徑：

1) 沿著方向，在參數重構中套用格式，因為在格式構造過程中沒有考慮方向的變化，嚴格的說計算模板式(7)中的整點通量時，平均值應該采用式(8)計算：

(8)

2) 同樣沿著方向的積分中使用格式，計算整點通量也只能采用儲存在-12,面心處的平均值：

(9)

以上是在控制體左右界面的計算過程，對于上下界面還應有重新定義的平均值。通過以上分析可知：盡管式(8)和式(9)的界面都在-12,處，形式相似，但是數學本質不同：在MUSCL和WENO格式中，-12,是自變量的當地取值，在G-FVM中，-12,是參數，積分函數的自變量是d。這種數學差異導致MUSCL和WENO格式用于計算G-FVM界面通量均不夠嚴謹。

Roe格式計算界面通量也使用了平均概念，最近Roe在文獻[19]中也提到這一格式不屬于FVM。因此，基于平均值構造的格式不能簡單的歸類于FVM。這類格式不能像FDM那樣直接采用Taylor級數分析，屬于FDM也不合適。MUSCL格式的原始文獻采用“守恒差分格式(Conservative Difference Scheme)”和“積分格式(Integration Scheme)”表示，本文認為“積分格式”這一定義更能準確反映這類格式的特點。

圖4 非結構網格中邊界面通量計算示意圖Fig.4 Schematic diagram of boundary surface flux calculation in unstructured grid

在非結構網格中應用MUSCL類格式，如圖4 所示，在已知格心位置的守恒變量和梯度后，計算邊界面通量時采用式(10)計算邊界面上的守恒變量平均值L、R：

(10)

3 如何讓MUSCL和WENO格式守恒

目前大部分高精度格式建立在均勻結構網格基礎上，其中的模板函數公式在非均勻條件下需要增加幾何權系數進行修正，以四階中心差分格式為例進行說明。

在空間一維均勻網格上，方向導數對應的差分格式為

(11)

將守恒變量換成原始變量(,,)、對流項通量，式(11)同樣成立，文獻[16]關于限制器函數相同的結論成立。

(12)

對式(12)的離散是在計算空間上進行的，理論上差分格式中的變量應該采用曲線坐標系下的相關參數。采用四階中心差分格式離散對流通量項，得到

(13)

此時，式(13)仍保持四階精度。但是如果采用式(11)直接離散守恒變量，得到

(14)

(15)

用GCL消除坐標變換引起的幾何誘導誤差，原則上不同的格式需要構建不同的GCL算法。根據相關研究文獻，在高精度格式方面，SCMM(Symmetrical Conservative Metric Method)理論完善，算例驗證能夠消除WCNS類格式的幾何誘導誤差，從提供的驗證算例看，借鑒SCMM構建的WENO類差分格式的GCL算法可以有效減少誤差，但是誤差曲線還是有明顯變化。

另外一種消除幾何誘導誤差的算法是從離散等價方程出發應用差分格式，具體理論推導和算例驗證可以參考文獻[8-9]，這里不再詳述。

4 課題組研究成果

基于以上對G-FVM和FDM存在本質差異的認識，課題組開展如下研究，初步取得了一些成果：

1) 利用FDM不需要幾何清理的特點，結合非結構網格有限差分法，建立了不生成表面網格就能進行流場計算的“扎染算法”，如圖5所示。傳統的笛卡爾有限體積法需要物面網格，即關聯點和點形成通量積分，如果有微小縫隙就無法計算，但是FDM本身不關心點和點之間的關系，只要尋找實體模型表面點就能計算。

2) 認識到G-FVM構造高精度邊界算法的困難。除了對現有邊界算法進行改進外，本文課題組還在研究直接在邊界點求解控制方程的全場統一算法，如圖6所示。把無滑移物面條件=0(物面法向速度為0)作為約束加入到控制方程中，利用非結構網格有限差分法直接計算物面網格點參數，而傳統CFD的內點計算只能計算

圖5 扎染算法示意圖Fig.5 Schematic diagram for tie-dye algorithm

圖6 全場統一求解算法示意圖Fig.6 Schematic diagram of whole-field unified algorithm

到第2層網格，需采用邊界條件假設得到物面參數。圖7為分別采用兩種方法模擬Prandtl-Myer流動得到的物面附近密度等值線圖，可以看出，全場統一求解算法的等值線更加光滑，定性分析更為合理。

圖7 不同方法得到的密度等值線Fig.7 Density contours for different methods

5 結 論

本文以3個問題為出發點，在對問題進行分析和討論的過程中得出以下結論：

1) 有限差分法和高斯積分型有限體積法在網格需求方面存在差異。有限差分法計算不需要封閉的控制體或積分面，而高斯積分型有限體積法無法在不封閉的控制體上進行積分運算。

2) 有限差分法和高斯積分型有限體積法在邊界條件處理方面存在差異。有限差分法在網格點上進行離散求解，第一類邊界條件可以準確地添加到計算表達式中。高斯積分型有限體積法建立在多個積分面圍成的控制體上，盡管在通量積分中可以嚴格滿足第一類邊界條件，但是梯度重構等計算過程會導致數值解不能滿足施加邊界條件的現象。

3) 與在多維非結構網格上直接構造的WENO格式不同，結構網格上的MUSCL格式和WENO格式存在應用于空間多維時出現的平均值多定義問題，不屬于高斯積分型有限體積法。

4) 結構網格上的MUSCL格式和WENO格式不守恒。

致 謝

作者在構思過程中多次向南京航空航天大學朱君和中國空氣動力研究與發展中心朱華君請教WENO格式方面的問題，在此表示感謝。