孔令發(fā)，劉 偉，董義道

(國(guó)防科技大學(xué) 空天科學(xué)學(xué)院，長(zhǎng)沙 410073)

0 引 言

相比多塊結(jié)構(gòu)化網(wǎng)格和笛卡爾網(wǎng)格，非結(jié)構(gòu)網(wǎng)格因其便捷靈活的生成方式，近年來逐漸成為CFD研究與應(yīng)用的熱點(diǎn)[1-4]，但伴隨著非結(jié)構(gòu)網(wǎng)格的廣泛使用，一個(gè)涉及網(wǎng)格生成和模擬精度的矛盾開始變得愈發(fā)突出[5-6]。因此，學(xué)者們不斷嘗試改進(jìn)基于非結(jié)構(gòu)網(wǎng)格的離散算法，來實(shí)現(xiàn)自動(dòng)化網(wǎng)格生成與精準(zhǔn)化數(shù)值模擬的統(tǒng)一。

目前，針對(duì)非結(jié)構(gòu)網(wǎng)格較常用的數(shù)值方法是具有二階精度的有限體積方法[7-10]，該方法因具有較好的數(shù)值表現(xiàn)而被應(yīng)用于很多知名商業(yè)軟件，如ANSYS公司開發(fā)的商業(yè)軟件Fluent，以及NASA開發(fā)的In-house軟件FUN3D等[11]。為了提高計(jì)算效率與流場(chǎng)分辨率，很多學(xué)者嘗試將此數(shù)值方法向更高精度推廣[12-16]。不同于二階精度，高精度非結(jié)構(gòu)有限體積方法的實(shí)現(xiàn)，除了需要引入更多的高斯積分點(diǎn)來實(shí)現(xiàn)高階精度通量積分外，還需要重構(gòu)控制體單元的高階導(dǎo)數(shù)。有關(guān)高階導(dǎo)數(shù)的重構(gòu)，Barth與Jespersen[17-18]首先提出了k-exact重構(gòu)算法，Ollivier-Gooch等[19-21]在此基礎(chǔ)上發(fā)展了具有高階精度的非結(jié)構(gòu)有限體積離散方法，并系統(tǒng)研究了有關(guān)高階精度非結(jié)構(gòu)FV離散的對(duì)流、粘性通量與源項(xiàng)積分過程，以及曲線邊界處理等一系列問題。

在高階離散中，不同模板對(duì)實(shí)際計(jì)算（尤其是梯度與高階導(dǎo)數(shù)重構(gòu)過程）具有更直接的影響，因此部分學(xué)者也開展了相應(yīng)的模板選擇研究。文獻(xiàn)[11]提出了一種“智能增加”(Smart Augmentation, SA)模板構(gòu)造方法。在共面模板的基礎(chǔ)上，從圍繞中心單元每個(gè)節(jié)點(diǎn)的共點(diǎn)模板中尋找了一個(gè)距離中心單元最近的單元。該方法完全基于距離信息，可能導(dǎo)致流場(chǎng)出現(xiàn)非物理解。為了克服這一問題，在SA模板的基礎(chǔ)上，Schwoppe等[11]繼續(xù)從共點(diǎn)模板中選擇單元，以改善最小二乘矩陣的條件數(shù)。Nishikawa在文獻(xiàn)[22]研究指出，盡管該方法帶來了一定程度的改進(jìn)，但在高度變形網(wǎng)格上其計(jì)算穩(wěn)定性較差。為克服這一問題，Nishikawa同樣基于共面模板增加額外單元，其中一個(gè)是通過增加額外單元改進(jìn)模板的對(duì)稱性；另外一種以減小最小二乘矩陣系統(tǒng)Frobenius范數(shù)的倒數(shù)為目標(biāo)，從而降低重構(gòu)梯度的數(shù)值，以提高數(shù)值模擬的穩(wěn)定性。但遺憾的是，這一模板構(gòu)造方法在一定程度上降低了數(shù)值模擬的離散精度。

不同于上述模板選擇方式，Xiong等[23-24]針對(duì)二階精度非結(jié)構(gòu)有限離散，提出了一種基于局部方向的模板選擇方式，使得模板選擇沿著兩個(gè)局部方向依次進(jìn)行。這種模板選擇方式在非結(jié)構(gòu)網(wǎng)格上初步體現(xiàn)出類似結(jié)構(gòu)網(wǎng)格的結(jié)構(gòu)化特征，易于并行，并且在各向同性的三角形網(wǎng)格上，其中一個(gè)局部方向接近壁面法向，另一個(gè)方向沿著流向，比較有效地反映了例如邊界層型流場(chǎng)的主要變化情況。但在高度各向異性三角形網(wǎng)格上，由于其中一個(gè)局部方向與偏離壁面法向間偏差相對(duì)較大，導(dǎo)致其穩(wěn)定性下降[24]。

在局部方向模板的基礎(chǔ)上，本文作者在前期工作中探索了一種基于計(jì)算域全局方向的模板選擇方法[25]。例如邊界層型流動(dòng)，由于沿著壁面法向的變化相對(duì)劇烈，于是希望重構(gòu)過程的模板單元能夠盡可能多地體現(xiàn)壁面法向的流動(dòng)信息。因此，其中一個(gè)全局方向可取當(dāng)?shù)乇诿娣ㄏ颍瑫r(shí)為了使模板具有更好的空間延展性，并保留結(jié)構(gòu)化特征，另一個(gè)方向可取流向。相比上述幾種重構(gòu)模板，全局方向的確定更加簡(jiǎn)便，擺脫了原始網(wǎng)格拓?fù)鋵?duì)模板構(gòu)型的影響，并且能夠準(zhǔn)確反映流場(chǎng)特征。在二階精度非結(jié)構(gòu)有限體積方法中，其計(jì)算誤差相比常用的共點(diǎn)與共面模板更低，并且所需的模板數(shù)量較少，在具有較高計(jì)算準(zhǔn)確性的同時(shí)，有效節(jié)約了計(jì)算開銷。

本文將全局方向模板推廣至具有三階精度非結(jié)構(gòu)FV求解器，以測(cè)試其在高階精度數(shù)值模擬中的計(jì)算效果。文章的總體結(jié)構(gòu)如下：第1節(jié)給出了高階精度非結(jié)構(gòu)有限體積方法的控制方程離散以及kexact重構(gòu)的基本過程；第2節(jié)主要討論了不同的模板選擇方法，并敘述了全局方向模板選擇方法的主要思想與基本過程，簡(jiǎn)要探討了在相對(duì)復(fù)雜的外形上確定全局方向的基本思路；第3節(jié)分別測(cè)試了全局方向模板在兩類數(shù)值算例中的表現(xiàn)情況，同時(shí)測(cè)試了局部方向模板在高階精度非結(jié)構(gòu)有限體積方法中的數(shù)值表現(xiàn)；第4節(jié)給出了本文結(jié)論與下一步工作展望。

1 非結(jié)構(gòu)有限體積方法

1.1 控制方法與離散

無粘流動(dòng)的積分型控制方程如下：

式中，u為守恒變量，F(xiàn)c(u)為對(duì)流通量，n為控制體邊界面外法向矢量，V與?V分別代表積分域與積分域邊界。在高階精度FV離散過程中，方程(1)通?？杀浑x散為：

圖1 非結(jié)構(gòu)有限體積高精度離散與通量構(gòu)造示意圖Fig.1 Sketch of high-order unstructured finite volume discretization and flux construction

1.2 k-exact重構(gòu)

在二階精度非結(jié)構(gòu)有限體積方法中，通常采用最小二乘方法重構(gòu)單元梯度，由于直接基于最小二乘方法進(jìn)行泰勒級(jí)數(shù)展開無法達(dá)到二階以上精度，因此本節(jié)采用k-exact方法來實(shí)現(xiàn)單元梯度與高階導(dǎo)數(shù)重構(gòu)。重構(gòu)函數(shù)的形式如下：

同時(shí)，其需要滿足的約束為：

因此，將重構(gòu)函數(shù)在單元j上積分可得：

式中nx為單元面單位外法矢量的x分量。

參照方程(6)，針對(duì)任意模板單元均可構(gòu)造一個(gè)重構(gòu)多項(xiàng)式，并最終組建矩陣方程組：

方程組的第一個(gè)方程是約束方程，除約束方程外，其余每個(gè)方程前均乘以權(quán)系數(shù) ωij：

該系數(shù)用來強(qiáng)調(diào)與中心控制體幾何上鄰近的模板單元對(duì)重構(gòu)結(jié)果的影響。在得到單元梯度與高階導(dǎo)數(shù)的基礎(chǔ)上，可通過約束方程得到單元變量點(diǎn)值：

該點(diǎn)值用來插值計(jì)算單元面高斯點(diǎn)處的左右狀態(tài)矢量，以構(gòu)造該點(diǎn)處的數(shù)值通量。

2 模板選擇方法

本節(jié)回顧了當(dāng)前應(yīng)用廣泛的共點(diǎn)、共面模板以及局部方向模板選擇方法的基本思路與存在的問題，并在此基礎(chǔ)上，給出了全局方向模板選擇方法的具體過程及其表現(xiàn)效果。

2.1 共點(diǎn)與共面模板

共點(diǎn)模板以及共面模板的應(yīng)用最為廣泛，其依賴于網(wǎng)格拓?fù)洌⑼ㄟ^相應(yīng)的模板層數(shù)來控制模板單元數(shù)量使其滿足重構(gòu)方程的要求，如圖2所示，共面模板由于網(wǎng)格單元面數(shù)量固定，因此隨著模板層數(shù)的遞增，模板單元數(shù)量增長(zhǎng)平緩，而共點(diǎn)模板的數(shù)量較難控制，并且隨著層數(shù)增加，單元數(shù)量快速上升。

圖2 共面模板與共點(diǎn)模板Fig.2 Face-neighbor and vertex-neighbor stencils

此外，兩種模板選擇方式均依賴于固定的網(wǎng)格拓?fù)潢P(guān)系，針對(duì)不同流動(dòng)無法有效捕捉到流場(chǎng)變化，并且兩種模板在滿足重構(gòu)要求的基礎(chǔ)上，包含過多冗余單元，尤其是進(jìn)行高階導(dǎo)數(shù)重構(gòu)時(shí)，所需要的模板數(shù)量較多，因此隨著層數(shù)的增加，冗余單元的數(shù)量也因此驟增，這將顯著增大計(jì)算開銷，降低求解效率。

2.2 局部方向模板

針對(duì)上述兩種模板存在的問題，Xiong等[23-24]構(gòu)建了基于兩個(gè)局部方向的模板選擇方法，這種模板選擇方法有效減少了用于重構(gòu)的模板單元數(shù)量，并且在各項(xiàng)同性或長(zhǎng)寬比較小的三角形網(wǎng)格上，兩個(gè)局部方向近似沿著壁面法向與流向，在二階精度非結(jié)構(gòu)有限體積求解器中取得了較好地?cái)?shù)值表現(xiàn)。

如圖3所示，對(duì)于四邊形單元，局部方向就是網(wǎng)格單元兩組對(duì)邊中點(diǎn)的連線；在處理三角形單元時(shí)，可將三角形單元的一個(gè)頂點(diǎn)視作由四邊形的一條邊退化而成，并將其稱為“退化點(diǎn)”，但由于存在三種可能的局部方向組合，因此需要采用陣面推進(jìn)方法[26]（Advancing Front Technique, AFT）來確定最終的局部方向組合方式。

圖3 四邊形與三角形單元的局部方向組合Fig.3 Local directions of quadrilateral and triangular grids

針對(duì)三角形網(wǎng)格，基于陣面推進(jìn)方法找到的局部方向的組合結(jié)果如圖4所示，從圖中可以看出，當(dāng)網(wǎng)格長(zhǎng)寬比較大時(shí)，其第一局部方向已偏離壁面法向。

圖4 三角形網(wǎng)格單元的局部方向組合Fig.4 Combination of local directions in a triangular grid

在確定局部方向的基礎(chǔ)上，可沿著兩個(gè)方向擴(kuò)充模板單元。如果中心單元為四邊形，沿著局部方向選擇的模板單元實(shí)際上是共面模板單元。

圖5展示了三角形單元的模板構(gòu)造方式。若局部方向通過“退化點(diǎn)”，例如圖中的P0點(diǎn)，此時(shí)選擇和局部方向夾角最大的單元；如果不通過“退化點(diǎn)”，則選擇相應(yīng)的共面模板。由于該模板選擇方法涉及陣面推進(jìn)與夾角判斷過程，相對(duì)增加了前處理過程的處理開銷。

圖5 三角形網(wǎng)格的模板單元確定Fig.5 Stencil selection on a triangular grid

此外，如圖6所示，當(dāng)網(wǎng)格長(zhǎng)寬比較小時(shí)，局部方向模板單元基本沿著壁面法向與流向，但在高度各向異性網(wǎng)格上，由于局部方向的偏離，導(dǎo)致沿兩個(gè)方向找到的模板空間延展性下降。經(jīng)測(cè)試，在此類網(wǎng)格上局部方向模板難以保證計(jì)算穩(wěn)定性[24]。因此，針對(duì)這一問題，迫切需要發(fā)展一種新的模板構(gòu)造方式，其能擺脫網(wǎng)格拓?fù)涞募s束，并準(zhǔn)確捕捉流動(dòng)的各向異性。

圖6 不同長(zhǎng)寬比網(wǎng)格上的局部方向模板示意圖Fig.6 Diagram of local-direction stencils on grids with different aspect ratios

2.3 全局方向模板

在局部方向模板的基礎(chǔ)上，我們針對(duì)二階精度非結(jié)構(gòu)有限體積方法發(fā)展了一種全局方向模板選擇方法，該方法不再依賴于每個(gè)網(wǎng)格單元內(nèi)部的局部方向，而是計(jì)算域的全局方向。因此如圖7所示，即使在大長(zhǎng)寬比三角形網(wǎng)格上，全局方向模板始終沿著壁面法向與流向，能夠準(zhǔn)確捕捉到流場(chǎng)的各項(xiàng)異性特征。

圖7 不同長(zhǎng)寬比網(wǎng)格上的全局方向模板示意圖Fig.7 Diagram of global-direction stencils on grids with different aspect ratios

具體而言，尋找全局方向模板單元可分為兩個(gè)步驟，首先需要找到與中心單元共點(diǎn)的模板單元集合，在此基礎(chǔ)上，過中心單元的幾何中心做兩條直線，找到集合中與兩條直線相交的模板單元即可。因此，模板單元數(shù)量可得到較為準(zhǔn)確地控制，并且整個(gè)過程并未引入任何復(fù)雜性。

針對(duì)簡(jiǎn)單外形，確定壁面法向的過程相對(duì)簡(jiǎn)便，例如圖7中給出的規(guī)則矩形計(jì)算域，其只包含一個(gè)壁面，確定全局方向時(shí)可直接取單元參考點(diǎn)處的壁面法向與切向。但對(duì)于更一般的外形，由于其物面不存在解析曲線，因此針對(duì)復(fù)雜外形確定計(jì)算域中單元參考點(diǎn)的全局方向時(shí)，可基于湍流模型中的壁函數(shù)計(jì)算首先獲得單元壁面距離[27-29]。在此基礎(chǔ)上，Jalali與Ollivier-Gooch[19]等提出了一種基于壁面距離來獲得壁面法向的有效手段，可采用該方法得到計(jì)算域中每一單元參考點(diǎn)處的壁面法向。此外，為使模板具有更好的空間延展性，另一個(gè)全局方向可取與壁面法向垂直的方向。

但在網(wǎng)格劃分好的前提下，對(duì)計(jì)算域中的每一個(gè)網(wǎng)格參考點(diǎn)計(jì)算一次壁面距離，并基于壁面距離確定壁面法向并不容易。這種處理方式會(huì)顯著增加算法復(fù)雜度與實(shí)現(xiàn)的困難程度，尤其在并行計(jì)算中。更重要的是，針對(duì)凹腔[30]這類局部包含多個(gè)壁面的情況，若在整個(gè)計(jì)算域內(nèi)均采用全局方向模板，會(huì)造成方向判斷混亂，導(dǎo)致方法無法實(shí)施。

因此針對(duì)復(fù)雜外形，較好的處理方法是在對(duì)物面附近網(wǎng)格單元進(jìn)行重構(gòu)時(shí)，采用全局方向模板，來捕捉邊界層內(nèi)部的流場(chǎng)變化，但在遠(yuǎn)離邊界層的各向同性區(qū)域仍采用依賴網(wǎng)格拓?fù)潢P(guān)系的共點(diǎn)或共面模板。這樣做的好處在于，在邊界層內(nèi)部可以有效捕捉流動(dòng)的各向異性特征，并減少重構(gòu)過程使用的模板單元數(shù)量；此外，針對(duì)凹腔等外形，考慮局域性后可保證方法的正確實(shí)施；同時(shí)，僅對(duì)物面附近的網(wǎng)格單元確定全局方向，可有效簡(jiǎn)化算法的實(shí)現(xiàn)過程，并且在并行計(jì)算中的處理也會(huì)更加簡(jiǎn)便。

本文首先在一些簡(jiǎn)單外形上驗(yàn)證了全局方向模板在高階精度非結(jié)構(gòu)有限體積方法中的數(shù)值表現(xiàn)，為下一步將方法向復(fù)雜外形推廣打下了基礎(chǔ)。

3 數(shù)值驗(yàn)證

為檢驗(yàn)不同重構(gòu)模板在高階精度非結(jié)構(gòu)有限體積求解器中的數(shù)值表現(xiàn)，本節(jié)通過基于制造解方法(Method of Manufactured Solutions, MMS)[31-33]的流動(dòng)以及超聲速渦流兩類數(shù)值算例，分別對(duì)比了共點(diǎn)、共面模板以及全局方向模板在三階精度求解器上的計(jì)算效果。由于2.2節(jié)分析了在各向異性的三角形網(wǎng)格上，基于局部方向找到的模板單元偏離壁面法向，導(dǎo)致其在二階精度求解器上，計(jì)算穩(wěn)定性下降[23]。本節(jié)為了進(jìn)一步檢驗(yàn)局部方向模模板在高階精度非結(jié)構(gòu)有限體積方法中應(yīng)用的可行性，在測(cè)試上述三種模板計(jì)算效果的基礎(chǔ)上，還測(cè)試了局部方向模板的數(shù)值表現(xiàn)。同時(shí)，為簡(jiǎn)化圖表中對(duì)不同模板的描述，我們分別使用 V-Stencil、F-Stencil、L-Stencil與 G-Stencil來表示共點(diǎn)、共面、局部方向模板與全局方向模板。

3.1 基于制造解方法的流動(dòng)

本節(jié)針對(duì)歐拉方程構(gòu)建了一個(gè)具有各向異性特征的MMS流動(dòng)[34]，并且解的形式如下：

因此，可首先根據(jù)解的形式得到源項(xiàng)s為：

由于本文探討三階精度數(shù)值模擬中的不同模板選擇方法，如果僅采用單元中心處的源項(xiàng)點(diǎn)值，無法達(dá)到二階以上精度。因此需要對(duì)源項(xiàng)在控制體單元上進(jìn)行數(shù)值積分，以保證其具有三階精度：

式中si與wi分別代表單元面中點(diǎn)處的源項(xiàng)與權(quán)系數(shù)，參照文獻(xiàn)[20]，wi可取由方程(10)得到的流場(chǎng)如圖8所示。

圖8 基于制造解的流場(chǎng)Fig.8 Flow fields with manufactured solutions

此外，我們分別測(cè)試了網(wǎng)格長(zhǎng)寬比AR= 5×102、1×103的規(guī)則與擾動(dòng)三角形網(wǎng)格，并且在每種長(zhǎng)寬比下，分別設(shè)置了五套由疏到密的網(wǎng)格以測(cè)試不同模板的計(jì)算精度。網(wǎng)格示意圖與背景四邊形網(wǎng)格量分布如圖9與表1所示。

表1 不同長(zhǎng)寬比下背景四邊形網(wǎng)格的分布量Table 1 The distribution of background quadrilateral grid cells with different aspect ratios

圖9 規(guī)則與擾動(dòng)三角形網(wǎng)格Fig.9 Regular and randomly perturbed triangular grids

為簡(jiǎn)化分析過程，并對(duì)比在高度各項(xiàng)異性網(wǎng)格上不同模板的表現(xiàn)，本節(jié)給出了長(zhǎng)寬比AR= 1×103時(shí)的計(jì)算結(jié)果。但在測(cè)試過程中，我們發(fā)現(xiàn)局部方向模板在長(zhǎng)寬比AR= 5×102與1×103的高度各向異性網(wǎng)格上計(jì)算發(fā)散。深入分析后發(fā)現(xiàn)，該方法在邊界附近找到的模板單元數(shù)量不足。

如圖10所示，對(duì)邊界單元C擴(kuò)充模板單元時(shí)，首先可沿第二局部方向，找到與單元C共面的模板單元C2，其次在沿第一局部方向?qū)ふ夷０鍐卧獣r(shí)，由于第一局部方向與CC2之間的夾角小于CC1，進(jìn)而沿第一局部方向找到的模板單元仍為C2，導(dǎo)致在邊界附近模板數(shù)量不足，需要不斷擴(kuò)充層數(shù)以保證正確求解邊界單元的重構(gòu)方程。因此后續(xù)的算例測(cè)試只給出了共點(diǎn)、共面模板與全局方向模板的計(jì)算結(jié)果。

圖10 邊界附近的局部方向模板單元擴(kuò)充Fig.10 The stencil augmentation of boundary-adjacent cells

從圖11與表2反映的計(jì)算結(jié)果來看，使用三種模板均能接近并達(dá)到求解器的設(shè)計(jì)精度，但在三種模板中，不論在規(guī)則網(wǎng)格還是添加隨機(jī)節(jié)點(diǎn)擾動(dòng)的網(wǎng)格上，全局方向模板得到的計(jì)算誤差更低，因此，該模板的使用能夠有效改善求解器的計(jì)算準(zhǔn)確性。

表2 規(guī)則密網(wǎng)格上的計(jì)算誤差以及最后兩套網(wǎng)格間的精度Table 2 Errors on the finest regular grid and computational accuracy between the last two sets of grids

圖11 規(guī)則網(wǎng)格與擾動(dòng)網(wǎng)格上的計(jì)算誤差Fig.11 Computational errors on the regular and perturbed grids

此外，由于三階精度非結(jié)構(gòu)FV方法的重構(gòu)過程需要求解6個(gè)未知量，因此，至少需要6個(gè)模板單元。但共點(diǎn)模板當(dāng)層數(shù)取1時(shí)，在邊界單元上無法滿足重構(gòu)方程組的求解，因此，需要將層數(shù)擴(kuò)充為2，而共面模板只有當(dāng)層數(shù)增加至4層時(shí)才可進(jìn)行流場(chǎng)計(jì)算。結(jié)合圖12可以明顯看出，相比這兩種模板，全局方向模板所需要的模板單元數(shù)量更少，因此在具有更高計(jì)算準(zhǔn)確性的同時(shí)，可相對(duì)降低重構(gòu)過程的計(jì)算開銷。

圖12 三種模板的平均模板單元數(shù)量Fig.12 Average stencil sizes of three different stencils

但在上述對(duì)比中，三種模板選擇方式均是基于邊界處的模板單元數(shù)量而給定整體的模板層數(shù)，從而造成共點(diǎn)與共面模板的單元數(shù)量在整個(gè)流場(chǎng)中明顯增加。在此基礎(chǔ)上，我們僅對(duì)邊界單元擴(kuò)充模板層數(shù)，而內(nèi)部單元的層數(shù)保持最低，這樣可有效縮減三種模板的規(guī)模。

在長(zhǎng)寬比AR= 1×103的網(wǎng)格上，通過上述方案得到三種模板的模板單元數(shù)量如表3所示。

為對(duì)比控制模板層數(shù)前后的計(jì)算誤差，兩種情況的誤差統(tǒng)計(jì)結(jié)果如圖13所示，并且修改層數(shù)后的模板選擇方法均添加標(biāo)識(shí)“improved”。

圖13 修改模板層數(shù)前后的計(jì)算誤差對(duì)比Fig.13 The comparison of computational errors before and after altering stencil layers

從圖中可以明顯看出，采用最小的模板層數(shù)后，三種模板的L2誤差均明顯降低，而L∞誤差相比縮減模板層數(shù)之前有所上升。但不論是否縮減模板層數(shù)，全局方向模板的誤差值在三種方法中均保持最低。同時(shí)結(jié)合表3中的數(shù)據(jù)可以看出，在對(duì)三種模板選擇方式縮減模板層數(shù)后，共點(diǎn)、共面模板，以及全局方向模板的模板單元數(shù)量均有所減少，但全局方向模板單元數(shù)量仍保持最少。因此，在模板數(shù)量差距較小的情況下，全局方向模板的計(jì)算誤差仍在三種方法中保持最低，有效證明了從完全多維的共點(diǎn)模板中提取具有分維特征的全局方向模板單元的必要性。

表3 三種模板選擇方式在不同網(wǎng)格上的平均模板單元數(shù)量Table 3 Average stencil sizes of three stencil selection methods on five sets of triangular grids

3.2 超聲速渦流動(dòng)

本節(jié)采用含曲線邊界的超聲速渦流動(dòng)，來進(jìn)一步檢驗(yàn)全局方向模板在真實(shí)流動(dòng)中的數(shù)值表現(xiàn)。該流動(dòng)具有解析解，其形式如下：

式中 γ為比熱比，‖v‖為速度的大小，并且其內(nèi)外徑分別為ri= 1以及r0= 1.384。式中在內(nèi)徑處的馬赫數(shù)Mi= 2.25，并且密度ρi= 1。此外聲速為：

該流動(dòng)的流場(chǎng)結(jié)構(gòu)如圖14所示。

圖14 超聲速渦流場(chǎng)Fig.14 Flow fields of the supersonic vortex flow

同樣，分別測(cè)試了兩種不同長(zhǎng)寬比網(wǎng)格下的計(jì)算結(jié)果，由于此算例包含曲線邊界，網(wǎng)格長(zhǎng)寬比不像直線邊界網(wǎng)格在全場(chǎng)統(tǒng)一。因此采用壁面第一層網(wǎng)格的長(zhǎng)寬比來近似代替。最疏一套網(wǎng)格示意圖，以及背景四邊形網(wǎng)格在徑向、周向的分布量如圖15與表4所示。該算例同樣測(cè)試了局部方向模板的計(jì)算效果，但與3.1節(jié)所反映出的計(jì)算結(jié)果一致，局部方向模板在該數(shù)值算例中出現(xiàn)計(jì)算發(fā)散的現(xiàn)象。因此，在統(tǒng)計(jì)誤差時(shí)只給出了共點(diǎn)、共面模板以及全局方向模板的計(jì)算結(jié)果。

圖15 規(guī)則與擾動(dòng)三角形網(wǎng)格Fig.15 Regular and randomly perturbed triangular grids

表4 不同長(zhǎng)寬比下的背景四邊形網(wǎng)格在徑向與周向的分布量Table 4 The distribution of background quadrilateral grid cells in the radial and circumferential directions

為簡(jiǎn)化結(jié)果分析，文中給出了長(zhǎng)寬比AR≈ 4.0時(shí)的計(jì)算結(jié)果。圖16與圖17分別為在規(guī)則網(wǎng)格和擾動(dòng)網(wǎng)格上，三種模板在全場(chǎng)以及壁面附近的誤差統(tǒng)計(jì)結(jié)果，表4中列舉的數(shù)據(jù)為規(guī)則密網(wǎng)格上的具體誤差值以及最后兩套網(wǎng)格間的計(jì)算精度。

從圖16中可以看出，不論在規(guī)則網(wǎng)格還是擾動(dòng)網(wǎng)格上，三種模板得到的計(jì)算結(jié)果均接近求解器的設(shè)計(jì)精度。從誤差值來看，全局方向模板的L2誤差在三種模板中最低，L∞誤差與共面模板接近，并稍高于共面模板，但低于共點(diǎn)模板。在此基礎(chǔ)上統(tǒng)計(jì)了壁面附近的誤差值，其中最大誤差通常位于壁面，而L∞誤差曲線已在全場(chǎng)中統(tǒng)計(jì)，因此壁面附近只給出了L2誤差曲線。

圖16 規(guī)則網(wǎng)格與擾動(dòng)網(wǎng)格上的全場(chǎng)壓力誤差Fig.16 Global pressure errors on the regular and perturbed grids

從圖17可以看出，三種模板在壁面附近的誤差接近，但結(jié)合表5的數(shù)據(jù)可以看出，其中共面模板的誤差最低，全局方向模板的誤差與共點(diǎn)模板接近并稍低于共點(diǎn)模板。

表5 規(guī)則密網(wǎng)格上的計(jì)算誤差以及最后兩套網(wǎng)格間的精度Table 5 Errors on the finest regular grid and computational accuracy between the last two sets of grids

圖17 規(guī)則網(wǎng)格與擾動(dòng)網(wǎng)格上的壁面壓力誤差Fig.17 Wall pressure errors on the regular and perturbed grids

此外，從圖18可以明顯看出，全局方向模板所需要的模板單元數(shù)量更少，相比共點(diǎn)模板與共面模板，模板數(shù)量分別減少了57.9%與49.6%，在具備較高計(jì)算準(zhǔn)確性的同時(shí)降低了重構(gòu)過程的計(jì)算開銷。

圖18 三種模板的平均模板單元數(shù)量Fig.18 Average stencil sizes of three different stencils

4 結(jié) 論

本文將全局方向模板推廣至高階精度非結(jié)構(gòu)有限體積方法，并通過基于制造解的數(shù)值流動(dòng)與真實(shí)超聲速渦流分別檢驗(yàn)了全局方向模板在具有三階精度非結(jié)構(gòu)有限體積求解器中的數(shù)值表現(xiàn)，得出以下結(jié)論：

1）首先從方法實(shí)現(xiàn)過程而言，全局方向模板選擇方法的過程簡(jiǎn)便，相比局部方向模板，全局方向模板單元的確定有效避免了陣面推進(jìn)與局部方向傳遞的繁瑣；相比常用的共點(diǎn)模板，新方法只需要在確定共點(diǎn)單元的的基礎(chǔ)上，找到與兩個(gè)全局方向相交的網(wǎng)格單元即可，該過程并未引入任何復(fù)雜性。

2）其次從計(jì)算結(jié)果來看，相比兩種常用的共點(diǎn)、共面模板，全局方向模板的使用可有效降低計(jì)算誤差，在測(cè)試的兩個(gè)數(shù)值算例中，全局方向模板均具有較好的數(shù)值表現(xiàn)。

3）此外，全局方向模板對(duì)不同長(zhǎng)寬比的網(wǎng)格具有較強(qiáng)的適應(yīng)性，擺脫了網(wǎng)格單元的拓?fù)浼s束，即使在大長(zhǎng)寬比三角形網(wǎng)格上，模板單元始終能夠保證沿著壁面法向與流向，具有較好的空間延展性。

4）最后，基于全局方向的模板選擇方法從完全多維模板中提取了具有分維特征的模板單元，保證空間延展性的同時(shí)，可有效減少重構(gòu)過程所需的模板單元數(shù)量。在具有三階精度的非結(jié)構(gòu)有限體積求解器中，相比共點(diǎn)與共面模板，分別減少模板總數(shù)的57.9%與49.6%，可節(jié)約計(jì)算開銷。

綜上，全局方向模板在高階精度非結(jié)構(gòu)有限體積方法中具有較好的數(shù)值表現(xiàn)，具備進(jìn)一步推廣與應(yīng)用的可行性。接下來的工作將從兩個(gè)方面分別開展：首先，從改善計(jì)算效率方面，針對(duì)計(jì)算域內(nèi)的不同網(wǎng)格單元，進(jìn)一步對(duì)比分析采用不同模板層數(shù)后幾種模板的計(jì)算效果，以控制重構(gòu)過程所需要的模板單元數(shù)量；其次，從方法推廣與實(shí)際應(yīng)用層面，應(yīng)推廣該重構(gòu)模板在二維、三維復(fù)雜外形以及粘性流動(dòng)中的使用。同時(shí)，粘性項(xiàng)的高精度離散同樣也是下一步研究與關(guān)注的重點(diǎn)。