郭 嘉，朱華君，石國權(quán)，燕振國，*，宋松和，唐玲艷

(1. 中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000；2. 國防科技大學 理學院 數(shù)學系，長沙 410073)

0 引 言

高階精度格式相比于低階格式在精細結(jié)構(gòu)模擬中更具優(yōu)勢，隨著工程問題對于精細結(jié)構(gòu)模擬要求的逐步提升，高階精度格式成為計算流體力學中的一個重點研究方向[1-3]。在高階精度格式中，重構(gòu)修正方法（correction procedure via reconstruction, CPR）具有格式緊致、高效、可應(yīng)用于復雜非結(jié)構(gòu)網(wǎng)格的優(yōu)點。它最初由Huynh[4]提出，并由王志堅等[5]推廣至非結(jié)構(gòu)網(wǎng)格上， 現(xiàn)被較多應(yīng)用到科學計算和工程預測之中。但CPR 在捕捉激波時容易產(chǎn)生虛假的數(shù)值振蕩[1]，故CPR 的激波捕捉技術(shù)成為研究焦點[1-2,4-5]。

舒其望等[6]在CPR 方法中引入WENO 限制器，保持了光滑區(qū)域上的高階精度特性，并在間斷附近有效抑制了數(shù)值振蕩。李萬愛等[7]通過提出一個p加權(quán)的過程，在間斷Galerkin 方法（discontinuous Galerkin,DG）上構(gòu)造了新的光滑因子，提高了激波捕捉性能，然而間斷附近的振蕩現(xiàn)象未被完全消除。具有子單元分辨率的人工黏性法[8]也被提出，但是它依賴于參數(shù)的人工選擇。

一個可行的解決辦法是構(gòu)造有限元方法（如CPR 或DG）和有限差分法（ finite difference method,FD）/有限體積法（finite volume method, FV）的混合格式，使得在間斷附近具有更好的局部激波捕捉能力。這種方法的思路是，在光滑區(qū)域使用高階有限元方法保持緊致和高分辨率的性質(zhì)，在間斷附近采用FD 或FV 提供魯棒的激波捕捉能力。Dumbser 等提出DGFV 混合格式[9]并使用二階FV 來捕捉激波。劉鐵鋼等[10]提出具有良好幾何靈活性和高計算效率的Runge-Kutta 間 斷 Galerkin 方 法（ the Runge-Kutta discontinuous Galerkin method, RKDG）和加權(quán)本質(zhì)無振蕩格式（weighted essentially non-oscillatory schemes,WENO）的區(qū)域混合格式。這些格式避免了高階有限元方法直接處理激波，而是使用具有良好激波捕捉能力的格式進行激波捕捉，從而使得混合格式具有良好的激波捕捉特性。

朱華君等[11]構(gòu)造了基于區(qū)域混合思路的混合WCNS-CPR 格式，解決了混合格式的空間精度和幾何守恒律等基礎(chǔ)問題，該格式計算效率較高且易于應(yīng)用到復雜曲網(wǎng)格中。最近，基于單元動態(tài)混合思想，郭嘉和朱華君等[12]構(gòu)造了WCNS-CPR 格式的激波捕捉技術(shù)。本文，我們在混合WCNS-CPR 格式基礎(chǔ)上，針對高階CPR 方法，通過引入激波偵測器和基于高階WCNS 插值的二階格式來模擬帶有激波的流動問題，進一步提高混合算法的計算效率并增強混合算法的激波捕捉能力。

1 數(shù)值方法

以一維守恒律方程（1）為例，對三階CPR 方法和高分辨率二階格式作簡要的介紹。

1.1 三階CPR 方法

1.2 高分辨率二階格式

WCNS 格式[13-14]的構(gòu)造過程主要包含三部分：非線性插值、數(shù)值通量以及通量導數(shù)的求解。文獻[15]中指出，基于對稱守恒型網(wǎng)格導數(shù)算法（ symmetrical conservative metric method， SCMM） 的WCNS 可分解為多個二階格式。本文采用了基于高階WCNS 插值的二階格式，此格式為三階WCNS 的分解格式之一，即非線性插值為高階WCNS 插值，而通量導數(shù)使用二階差分算子的格式，這里簡稱為高分辨率二階格式（ high-resolution second order finite difference scheme，HRFD2）。

1.2.1 三階WCNS 非線性插值

1.3 混合格式

本文采用與高階WCNS-CPR 混合算法[12]中相同的混合策略，關(guān)鍵在于將WCNS 格式替換為基于WCNS 插值的二階格式HRFD2，進一步提高計算效率并增強混合算法的激波捕捉能力。下面對混合格式中的格式切換策略和界面處理方法進行介紹。

1.3.1 格式切換策略

混合格式在不光滑區(qū)域使用具有激波捕捉能力的HRFD2，而在光滑區(qū)域使用具有較高計算效率和較高分辨率的CPR 格式。因此需要激波偵測器來偵測激波所在位置。本文使用NI 激波偵測器[18]來判斷激波所在單元，并將其標記為問題單元（Troubled cell）。NI的公式為：

其中：r+1 是非線性加權(quán)的候選模板個數(shù)，對于HRFD2，r=1。為了保證激波不會移動到光滑單元（Smooth cell），在問題單元附近引入緩沖單元（Buffer cell），如圖1 所示。

圖1 三種單元類型示意圖Fig. 1 Illustration of three different types of cells

所采用的格式切換策略如下：當激波偵測器偵測到激波時，在激波附近的單元格式切換為HRFD2。當激波離開后，這些單元的格式由HRFD2 切換為CPR。綜上，即根據(jù)激波偵測器的值，將所有單元分為不光滑單元（問題單元和緩沖單元）和光滑單元，在不光滑單元采用HRFD2，而在光滑單元采用CPR。

1.3.2 界面處理

為了避免兩種格式間出現(xiàn)錯誤的信息交換，我們需要仔細研究界面處理，主要的困難在于HRFD2-CPR 界面的處理。下面以一維情況為例，給出界面處理方案。首先，如圖2 所示，第i+1 單元的HRFD2 需要其相鄰的第i單元（CPR 單元）的虛擬點信息{gi,k,k=1,2,3}。HRFD2 單元在虛擬點的函數(shù)值是通過使用CPR 的信息 {ci,k,k=1,2,3}插值得到的，即：

圖2 界面處理Fig. 2 Schematic of interface treatment

此外，我們也要考慮格式切換時的信息傳遞。當?shù)趇個單元在第tn時刻是光滑單元，而在下一時刻為不光滑單元時，HRFD2 的求解點值可以通過CPR 的拉格朗日多項式插值得到。而對于相反的情形，即第tn時刻為不光滑單元，下一時刻切換為光滑單元時，信息應(yīng)該由HRFD2 傳遞給CPR。這個過程涉及到了跨界面非線性插值，我們以第i個單元的第1 個求解點為例：

2 數(shù)值算例

本節(jié)中我們將會測試CPR 與HRFD2 的混合格式（簡記為CPR-HRFD2）的空間精度、激波捕捉能力和計算效率。同時會進行CPR-HRFD2 與HRFD2 的對比。

我們首先在二維等熵渦問題中測試CPRHRFD2 的精度，簡單對比了HRFD2 和CPR 的計算效率，并在一維激波管問題中驗證格式的激波捕捉性能，在二維黎曼問題和雙馬赫反射問題測試激波捕捉能力以及對比計算效率，隨后用高馬赫數(shù)問題來測試格式的強激波捕捉能力，最后在激波-旋渦干擾問題中研究對于復雜流動問題的模擬。時間推進格式為顯式三階龍格庫塔格式。L∞和L2誤差分別用來度量誤差的局部和全局值。L2誤差計算公式為：

2.1 精度測試

我們在二維等熵渦問題中對構(gòu)造的混合格式進行精度測試。下面考慮二維歐拉方程：

圖3、圖4、圖5、圖6 分別展示了CPR-HRFD2 的不同單元類型分布、x方向NI值、y方向NI值以及誤差分布圖。從圖4、圖5 可看出，當NI值大于 δNI時，單元會被判定為問題單元，單元分布圖中這個單元顯示為紅色，表示此單元由HRFD2 計算。緩沖單元也顯示為紅色，由HRFD2 計算。反之，當NI值小于δNI時，單元會被判定為光滑單元，單元分布圖中這個單元則顯示為藍色，此單元由CPR 計算。從表1可知，CPR-HRFD2 具有二階精度。

表1 二維精度測試Table 1 2D accuracy test

圖3 不同單元分布(自由度為120)Fig. 3 Distribution of different cells, DOFs = 120

圖4 x 方向NI 偵測結(jié)果分布(自由度為120)Fig. 4 NI detection result in x-direction, DOFs = 120

圖5 y 方向NI 偵測結(jié)果分布(自由度為120)Fig. 5 NI detection result in y-direction, DOFs = 120

圖6 誤差分布(自由度為120)Fig. 6 Error distribution, DOFs = 120

此外，針對二維等熵渦問題，進行了HRFD2 與CPR 的計算時間對比測試，取相同的自由度，時間步長ΔT=0.000 1， 計算到T=1。由表1 誤差分布和表2計算時間可知，在相同計算誤差下CPR 計算效率比HRFD2 高。

表2 二維等熵渦問題HRFD2 與CPR 計算時間對比Table 2 CPU time of HRFD2 and CPR in isentropic vortex problem

2.2 激波捕捉能力測試

2.2.1 激波管問題

考慮一維歐拉方程：

這一算例，采用基于解析解的Dirichlet 邊界條件。δNI=0.2 ， 時間步長為 ΔT=0.000 1， 計算到T=3.0。參考解由三階WCNS 在2100 自由度網(wǎng)格上計算得到。

由圖7 可知，激波和接觸間斷都可以被高分辨率地捕捉到。在接觸間斷和激波附近，沒有明顯的振蕩，這展現(xiàn)了CPR-HRFD2 良好的激波捕捉能力。此外，大部分單元由CPR 計算，保證具有更高的計算效率。

圖7 Sod 激波管問題，自由度為120Fig. 7 Sod shock tube problem, DOFs = 120

2.2.1.2 Shu-Osher 激波管問題[20]

初值為這一算例，采用Dirichlet 邊界條件。時間步長為ΔT=0.000 1， 參數(shù) δNI=0.2， 計算到T=1.8。參考解由三階WCNS 在9 000 自由度網(wǎng)格上計算得到。

由圖8 可知，在激波偵測器判別的充分光滑區(qū)域內(nèi)使用了CPR。此外，CPR-HRFD2 可以很好地捕捉激波和間斷。密度的數(shù)值解貼近參考解，展現(xiàn)了高分辨率特點。

圖8 Shu-Osher 激波管問題，自由度為600Fig. 8 Shu-Osher shock tube problem, DOFs = 600

2.2.1.3 Lax 激波管問題[21]

由圖9 可知，在x≈2.3和x≈3.8處，激波和間斷被高分辨率地捕捉到，整個流場沒有觀察到明顯振蕩。

圖9 Lax 激波管問題，自由度為300Fig. 9 Lax shock tube problem, DOFs = 300

2.2.2 二維黎曼問題

計算區(qū)域(x,y)∈[0,1]×[0,1]，初值為

計算到T=0.8， 取參數(shù) δNI=0.1。

圖10 為HRFD2 結(jié)果。由圖11、圖12 和圖13 可知，在光滑區(qū)域，流場由CPR 精確計算，且無明顯振蕩。此外，小尺度結(jié)構(gòu)也被CPR-HRFD2 高分辨率地捕捉到了，結(jié)果圖（圖11）和HRFD2 的結(jié)果圖（圖10）相近。大部分區(qū)域由CPR 計算，減少了計算代價。

圖10 HRFD2 密度計算結(jié)果，自由度為960×960Fig. 10 Density result of HRFD2, DOFs = 960×960

圖11 CPR-HRFD2 的計算密度結(jié)果，自由度為960×960Fig. 11 Density contour calculated by CPR-HRFD2 method，DOFs = 960×960

圖12 CPR-HRFD2 的x 方向NI 偵測Fig. 12 NI detection contour by CPR-HRFD2 method in x direction

圖13 CPR-HRFD2 的y 方向NI 偵測Fig. 13 NI detection contour by CPR-HRFD2 method in y direction

2.2.3 雙馬赫反射問題[22]

由圖14、圖15 和圖16 知，CPR-HRFD2 的結(jié)果圖和HRFD2 以及三階WCNS 的結(jié)果圖基本相同，此外，重要的小尺度和大尺度結(jié)構(gòu)都被CPR-HRFD2 捕捉到了，其他區(qū)域也由CPR 精確計算。由圖17 和圖18知，整個區(qū)域CPR 單元個數(shù)遠超過HRFD2 個數(shù)。另外，由表3 知，CPR-HRFD2 計算效率高于HRFD2（計算時間為18 480 s）和三階WCNS（計算時間為18 840 s），當 δNI=0.9時，相比于HRFD2 提高了34.7%的計算效率，相比于三階WCNS 提高了36.0%的計算效率。

表3 雙馬赫反射問題的CPU 時間對比（自由度為1 200×300）Table 3 CPU time comparison of double Mach problem, DOFs = 1 200×300

圖14 三階WCNS 格式密度計算結(jié)果，自由度為1 200×300Fig. 14 Density contour by WCNS3 method, DOFs = 1 200×300

圖15 HRFD2 密度計算結(jié)果，自由度為1 200×300Fig. 15 Density contour by HRFD2 method, DOFs = 1 200×300

圖16 CPR-HRFD2 密度計算結(jié)果，自由度為1 200×300Fig. 16 Density contour by CPR-HRFD2 method,DOFs = 1 200×300

圖17 CPR-HRFD2 的x 方向NI 偵測Fig. 17 NI detection contour by CPR-HRFD2 method in x direction

圖18 CPR-HRFD2 的y 方向NI 偵測Fig. 18 NI detection contour by CPR-HRFD2 method in y direction

2.2.4 馬赫數(shù) 2000 射流問題

我們研究馬赫數(shù)2 000 射流問題，計算域為[0,1]×[-0.25,0.25]， 參 數(shù) γ=5/3。 初 值 為(ρ,u,v,p)=(0.5,0,0,0.412 7)。右邊界、上邊界和下邊界是出流條件，左邊 界 的 入 流 條 件 為：若y∈[-0.05,0.05]，(ρ,u,v,p)=(5,800,0,0.412 7)； 否則 (ρ,u,v,p)=(5,0,0,0.412 7)。對于這個問題，噴流速度為800，相對于射流氣體中的聲速來說，大約是馬赫數(shù)2 100。我們用自由度1 710×855，計算到T=0.001， 時間步長 ΔT=5.0×10-4Δx,畫出密度和壓力的對數(shù)結(jié)果圖。文獻[23]中使用帶保正限制器的五階WCNS 格式對這個問題進行了模擬。如圖19 和圖20 所示，數(shù)值模擬結(jié)果與文獻[23]、文獻[24]相符，展現(xiàn)了CPR-HRFD2 格式的強激波捕捉能力。此外由圖21 和圖22 可知，大部分區(qū)域使用CPR 計算，保持了高計算效率。

圖19 CPR-HRFD2 密度計算結(jié)果，自由度1 710×855Fig. 19 Density contour by CPR-HRFD2 method, DOFs = 1 710×855

圖20 CPR-HRFD2 壓力計算結(jié)果，自由度1 710×855Fig. 20 Pressure contour by CPR-HRFD2 method,DOFs = 1 710×855

圖21 CPR-HRFD2 的x 方向NI 偵測Fig. 21 NI detection contour by CPR-HRFD2 method in x direction

圖22 CPR-HRFD2 的y 方向NI 偵測Fig. 22 NI detection contour by CPR-HRFD2 method in y direction

2.2.5 激波-旋渦干擾問題

在最后，我們測試了激波-旋渦干擾問題[25]。該問題包含了一個運動的旋渦和一道靜止正激波，可以發(fā)展出具有平滑特征和間斷的復雜流動結(jié)構(gòu)。如圖23 所示，該問題的計算域 (x,y)∈[0,2]×[0,1]。初始條 件 由 位 于 (0.5,y)的 靜 態(tài) 激 波Mas和 中 心 位 于(xc,yc)=(0.25,0.5)的 復合旋渦Mav給出，初始條件的具體設(shè)置參考文獻[25]。

圖23 激波-旋渦干擾的初始條件和邊界條件Fig. 23 Initial and boundary conditions of shock-vortex

左右邊界分別設(shè)為流入、流出邊界條件，上下邊界設(shè)為滑移固壁邊界條件。時間步長 ΔT=0.000 1，計算停止時間T=0.7， 參數(shù) δNI=0.9。

如圖24、圖25 所示，CPR-HRFD2 在保持激波捕捉能力的同時，對平滑旋渦產(chǎn)生的小尺度旋渦結(jié)構(gòu)也有很好的捕捉能力，計算結(jié)果與文獻[25]中的計算結(jié)果相符合。由圖26 和圖27 可知，大部分區(qū)域使用CPR 計算，保持了較高的計算效率。

圖24 HRFD2 密度計算結(jié)果，自由度960×480Fig. 24 Density contour by HRFD2 method, DOFs = 960×480

圖25 CPR-HRFD2 密度計算結(jié)果，自由度960×480Fig. 25 Density contour by CPR-HRFD2 method,DOFs = 960×480

圖26 CPR-HRFD2 的x 方向NI 偵測Fig. 26 NI detection contour by CPR-HRFD2method in x direction

激波-旋渦干擾問題包含復雜的流動結(jié)構(gòu)。按照如圖28 所示的切線位置，可以對不同結(jié)構(gòu)特征進行定量分析[26]。不同參考線位置對應(yīng)著不同的主要特征，見表4。

圖28 參考線位置Fig. 28 Illustration of the position of the reference lines

表4 參考線Table 4 Reference lines

為對比CPR-HRFD2 和HRFD2 的激波捕捉能力和分辨率特性，在此選?、佗冖蹍⒖季€進行分析。

從圖29 中①可看出，CPR-HRFD2 和HRFD2 都有較強的激波捕捉能力，激波附近和波后區(qū)域流場光滑，無明顯的振蕩。從圖29 中①③則可以明顯看出，對于旋渦的捕捉，HRFD2 的分辨率不如CPR-HRFD2。從圖26 和圖27 偵測的結(jié)果來看，CPR-HRFD2 在旋渦處也有使用HRFD2 計算，但是與完全使用HRFD2相比，旋渦結(jié)構(gòu)附近更多地使用CPR 求解，且相比于HRFD2，CPR 更高的精度避免了更大的耗散。如果增加 δNI，減少或避免偵測到旋渦結(jié)構(gòu)，分辨率特性差異將會更加明顯。

圖29 HRFD2 和CPR-HRFD2 的參考線結(jié)果Fig. 29 The reference lines calculated by HRFD2 and CPR-HRFD2

另外，文獻[26]中采用基于無參數(shù)梯度限制器的三階和四階CPR 方法對該問題進行了模擬，數(shù)值結(jié)果中觀察到數(shù)值振蕩，在定激波的右側(cè)產(chǎn)生許多高頻數(shù)值振蕩。而在本文的數(shù)值結(jié)果中不管是密度圖還是切線圖都沒有看到明顯的數(shù)值振蕩現(xiàn)象，說明了本文對高階CPR 方法所發(fā)展的基于混合算法的激波捕捉格式對激波旋渦問題具有更好的模擬能力。

3 結(jié) 論

本文發(fā)展了三階CPR 與HRFD2 的混合算法，使用基于非線性權(quán)的激波偵測器來判斷問題單元和緩沖單元，提出了格式切換和界面處理策略。在問題單元和緩沖單元內(nèi)，HRFD2 用來捕捉激波，而在光滑單元內(nèi)使用CPR，用于提高計算效率。文中進行了精度測試，激波捕捉能力測試，以及計算效率測試，得到了如下結(jié)論：

1）二維的精度測試結(jié)果表明CPR-HRFD2 保持了精度；

2）通過對一維激波管問題，二維黎曼問題，雙馬赫反射問題，高馬赫數(shù)射流問題和激波-旋渦干擾問題進行數(shù)值模擬，取得了比較好的結(jié)果，表明該混合格式具有強激波捕捉能力，可以應(yīng)用于高超聲速的高效數(shù)值模擬；

3）CPR-HRFD2 是一種高分辨率激波捕捉格式，同時具有較高計算效率。相比于HRFD2，CPR-HRFD2計算量更小，計算效率更高。

本文的構(gòu)造思想可以直接推廣到更高階的CPR方法與基于更高階WCNS 插值的二階格式的混合算法中。