應(yīng)用維數(shù)分裂方法推廣MUSCL和WENO格式的若干問(wèn)題

劉君，韓芳，魏雁昕

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

2019年，在國(guó)家數(shù)值風(fēng)洞工程基礎(chǔ)研究課題的支持下，作者與國(guó)內(nèi)研究高超聲速邊界層轉(zhuǎn)捩問(wèn)題的專(zhuān)家學(xué)者合作開(kāi)展激波裝配法在邊界層轉(zhuǎn)捩感受性問(wèn)題中的應(yīng)用研究，期望利用裝配法消除采用捕捉法數(shù)值模擬感受性問(wèn)題時(shí)在激波附近出現(xiàn)的虛假數(shù)值噪聲，以發(fā)展全場(chǎng)一致的高精度數(shù)值算法。

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

根據(jù)針對(duì)對(duì)象的不同，幾何守恒律又可分為體積守恒律和面積守恒律，體積守恒律對(duì)應(yīng)于運(yùn)動(dòng)及變形網(wǎng)格，面積守恒律對(duì)應(yīng)于靜止網(wǎng)格，本文所討論問(wèn)題皆基于靜止網(wǎng)格。在有限體積法(Finite Volume Method，F(xiàn)VM)中，與網(wǎng)格相關(guān)的幾何參數(shù)可以精確地求得，幾何守恒律自動(dòng)滿(mǎn)足，而有限差分法(Finite Difference Method，F(xiàn)DM)則存在幾何守恒律不能得到滿(mǎn)足的情況。因此，關(guān)于WENO格式守恒性的討論也可以變成對(duì)WENO格式是否屬于有限體積法的討論。“WENO格式是否屬于有限體積法”是本文提出的第1個(gè)問(wèn)題。

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

在數(shù)值實(shí)驗(yàn)及結(jié)果分析中發(fā)現(xiàn)從一維Euler方程出發(fā)介紹WENO格式時(shí)，如果模板作用于對(duì)流項(xiàng)通量，即便在非均勻網(wǎng)格上也不存在GCL問(wèn)題。但是對(duì)于多維空間，即使是均勻正方形網(wǎng)格也會(huì)出現(xiàn)與幾何參數(shù)相關(guān)的均值計(jì)算問(wèn)題。由于WENO格式模板函數(shù)比較復(fù)雜，推導(dǎo)過(guò)程選擇形式較為簡(jiǎn)單的MUSCL格式。

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

1 G-FVM和FDM是否等價(jià)

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

(1)

式中：為守恒變量；、、分別為、、方向的對(duì)流通量，其表達(dá)式為

=[,,,,]

=[,+,,,(+)]

=[,,+,,(+)]

=[,,,+,(+)]

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

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

(2)

當(dāng)網(wǎng)格靜止時(shí)，

=

式中：為Jacobian行列式，有

、、、、、、、、為度量系數(shù)，可通過(guò)以下公式求得：

與FDM不同，G-FVM可以在曲線(xiàn)坐標(biāo)系下直接應(yīng)用，采用高斯公式把體積分變成面積分：

(3)

式中：是控制體體積；是控制體外表面；d表示微元體積；d表示外表面微元面積；=++表示外表面微元外法向單位矢量。對(duì)流項(xiàng)通量積分的函數(shù)可寫(xiě)為

()=++

(4)

為便于幾何表達(dá)和簡(jiǎn)潔符號(hào)，采用圖1所示二維空間模型，其可看做是三維網(wǎng)格在方向的等值面，灰色曲線(xiàn)表示建立FDM的網(wǎng)格點(diǎn)，黑色曲線(xiàn)構(gòu)成G-FVM的控制體。

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

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

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

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

除此之外，本文進(jìn)一步補(bǔ)充如下論據(jù)：

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

對(duì)于二階精度的G-FVM，控制體梯度是利用相鄰網(wǎng)格點(diǎn)參數(shù)重構(gòu)出來(lái)的，不能保證為0，難以滿(mǎn)足絕熱壁條件。根據(jù)格心值和格心到壁面矢量重構(gòu)得到的壁面溫度分布不是常數(shù)，最多只能在一個(gè)點(diǎn)滿(mǎn)足等溫壁條件。

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

通過(guò)以上論述可知：G-FVM和FDM存在本質(zhì)差異。

圖1 二維曲線(xiàn)網(wǎng)格Fig.1 Two-dimensional curvilinear grid

圖2 物面附近網(wǎng)格示意圖Fig.2 Schematic diagram of grid near surface

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

按照G-FVM和FDM的特點(diǎn)，本文認(rèn)為MUSCL類(lèi)和WENO格式不屬于G-FVM，下面通過(guò)G-FVM界面通量計(jì)算過(guò)程提供更多的論據(jù)。

以圖3所示直角網(wǎng)格為研究對(duì)象，其中陰影部分為控制體，其流場(chǎng)參數(shù)的平均值為

(5)

圖3 均勻直角網(wǎng)格Fig.3 Uniform rectangular grid

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

(6)

目前看到MUSCL和WENO格式的構(gòu)造過(guò)程中只有線(xiàn)積分，從一維Euler方程出發(fā)，在半點(diǎn)界面處通量的通用形式為

(7)

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

1) 沿著方向，在參數(shù)重構(gòu)中套用格式，因?yàn)樵诟袷綐?gòu)造過(guò)程中沒(méi)有考慮方向的變化，嚴(yán)格的說(shuō)計(jì)算模板式(7)中的整點(diǎn)通量時(shí)，平均值應(yīng)該采用式(8)計(jì)算：

(8)

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

(9)

以上是在控制體左右界面的計(jì)算過(guò)程，對(duì)于上下界面還應(yīng)有重新定義的平均值。通過(guò)以上分析可知：盡管式(8)和式(9)的界面都在-12,處，形式相似，但是數(shù)學(xué)本質(zhì)不同：在MUSCL和WENO格式中，-12,是自變量的當(dāng)?shù)厝≈?，在G-FVM中，-12,是參數(shù)，積分函數(shù)的自變量是d。這種數(shù)學(xué)差異導(dǎo)致MUSCL和WENO格式用于計(jì)算G-FVM界面通量均不夠嚴(yán)謹(jǐn)。

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

圖4 非結(jié)構(gòu)網(wǎng)格中邊界面通量計(jì)算示意圖Fig.4 Schematic diagram of boundary surface flux calculation in unstructured grid

在非結(jié)構(gòu)網(wǎng)格中應(yīng)用MUSCL類(lèi)格式，如圖4 所示，在已知格心位置的守恒變量和梯度后，計(jì)算邊界面通量時(shí)采用式(10)計(jì)算邊界面上的守恒變量平均值L、R：

(10)

3 如何讓MUSCL和WENO格式守恒

目前大部分高精度格式建立在均勻結(jié)構(gòu)網(wǎng)格基礎(chǔ)上，其中的模板函數(shù)公式在非均勻條件下需要增加幾何權(quán)系數(shù)進(jìn)行修正，以四階中心差分格式為例進(jìn)行說(shuō)明。

在空間一維均勻網(wǎng)格上，方向?qū)?shù)對(duì)應(yīng)的差分格式為

(11)

將守恒變量換成原始變量(,,)、對(duì)流項(xiàng)通量，式(11)同樣成立，文獻(xiàn)[16]關(guān)于限制器函數(shù)相同的結(jié)論成立。

(12)

對(duì)式(12)的離散是在計(jì)算空間上進(jìn)行的，理論上差分格式中的變量應(yīng)該采用曲線(xiàn)坐標(biāo)系下的相關(guān)參數(shù)。采用四階中心差分格式離散對(duì)流通量項(xiàng)，得到

(13)

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

(14)

(15)

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

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

4 課題組研究成果

基于以上對(duì)G-FVM和FDM存在本質(zhì)差異的認(rèn)識(shí)，課題組開(kāi)展如下研究，初步取得了一些成果：

1) 利用FDM不需要幾何清理的特點(diǎn)，結(jié)合非結(jié)構(gòu)網(wǎng)格有限差分法，建立了不生成表面網(wǎng)格就能進(jìn)行流場(chǎng)計(jì)算的“扎染算法”，如圖5所示。傳統(tǒng)的笛卡爾有限體積法需要物面網(wǎng)格，即關(guān)聯(lián)點(diǎn)和點(diǎn)形成通量積分，如果有微小縫隙就無(wú)法計(jì)算，但是FDM本身不關(guān)心點(diǎn)和點(diǎn)之間的關(guān)系，只要尋找實(shí)體模型表面點(diǎn)就能計(jì)算。

2) 認(rèn)識(shí)到G-FVM構(gòu)造高精度邊界算法的困難。除了對(duì)現(xiàn)有邊界算法進(jìn)行改進(jìn)外，本文課題組還在研究直接在邊界點(diǎn)求解控制方程的全場(chǎng)統(tǒng)一算法，如圖6所示。把無(wú)滑移物面條件=0(物面法向速度為0)作為約束加入到控制方程中，利用非結(jié)構(gòu)網(wǎng)格有限差分法直接計(jì)算物面網(wǎng)格點(diǎn)參數(shù)，而傳統(tǒng)CFD的內(nèi)點(diǎn)計(jì)算只能計(jì)算

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

圖6 全場(chǎng)統(tǒng)一求解算法示意圖Fig.6 Schematic diagram of whole-field unified algorithm

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

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

5 結(jié) 論

本文以3個(gè)問(wèn)題為出發(fā)點(diǎn)，在對(duì)問(wèn)題進(jìn)行分析和討論的過(guò)程中得出以下結(jié)論：

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

2) 有限差分法和高斯積分型有限體積法在邊界條件處理方面存在差異。有限差分法在網(wǎng)格點(diǎn)上進(jìn)行離散求解，第一類(lèi)邊界條件可以準(zhǔn)確地添加到計(jì)算表達(dá)式中。高斯積分型有限體積法建立在多個(gè)積分面圍成的控制體上，盡管在通量積分中可以嚴(yán)格滿(mǎn)足第一類(lèi)邊界條件，但是梯度重構(gòu)等計(jì)算過(guò)程會(huì)導(dǎo)致數(shù)值解不能滿(mǎn)足施加邊界條件的現(xiàn)象。

3) 與在多維非結(jié)構(gòu)網(wǎng)格上直接構(gòu)造的WENO格式不同，結(jié)構(gòu)網(wǎng)格上的MUSCL格式和WENO格式存在應(yīng)用于空間多維時(shí)出現(xiàn)的平均值多定義問(wèn)題，不屬于高斯積分型有限體積法。

4) 結(jié)構(gòu)網(wǎng)格上的MUSCL格式和WENO格式不守恒。

致 謝

作者在構(gòu)思過(guò)程中多次向南京航空航天大學(xué)朱君和中國(guó)空氣動(dòng)力研究與發(fā)展中心朱華君請(qǐng)教WENO格式方面的問(wèn)題，在此表示感謝。