CPR 方法與高分辨率二階格式的混合算法

郭 嘉，朱華君，石國權，燕振國，*，宋松和，唐玲艷

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

0 引 言

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

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

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

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

1 數值方法

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

1.1 三階CPR 方法

1.2 高分辨率二階格式

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

1.2.1 三階WCNS 非線性插值

1.3 混合格式

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

1.3.1 格式切換策略

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

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

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

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

1.3.2 界面處理

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

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

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

2 數值算例

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

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

2.1 精度測試

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

圖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 偵測結果分布(自由度為120)Fig. 4 NI detection result in x-direction, DOFs = 120

圖5 y 方向NI 偵測結果分布(自由度為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 自由度網格上計算得到。

由圖7 可知，激波和接觸間斷都可以被高分辨率地捕捉到。在接觸間斷和激波附近，沒有明顯的振蕩，這展現了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， 參數 δNI=0.2， 計算到T=1.8。參考解由三階WCNS 在9 000 自由度網格上計算得到。

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

圖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 二維黎曼問題

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

計算到T=0.8， 取參數 δNI=0.1。

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

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

圖11 CPR-HRFD2 的計算密度結果，自由度為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 的結果圖和HRFD2 以及三階WCNS 的結果圖基本相同，此外，重要的小尺度和大尺度結構都被CPR-HRFD2 捕捉到了，其他區域也由CPR 精確計算。由圖17 和圖18知，整個區域CPR 單元個數遠超過HRFD2 個數。另外，由表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 格式密度計算結果，自由度為1 200×300Fig. 14 Density contour by WCNS3 method, DOFs = 1 200×300

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

圖16 CPR-HRFD2 密度計算結果，自由度為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 馬赫數 2000 射流問題

我們研究馬赫數2 000 射流問題，計算域為[0,1]×[-0.25,0.25]， 參 數 γ=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，相對于射流氣體中的聲速來說，大約是馬赫數2 100。我們用自由度1 710×855，計算到T=0.001， 時間步長 ΔT=5.0×10-4Δx,畫出密度和壓力的對數結果圖。文獻[23]中使用帶保正限制器的五階WCNS 格式對這個問題進行了模擬。如圖19 和圖20 所示，數值模擬結果與文獻[23]、文獻[24]相符，展現了CPR-HRFD2 格式的強激波捕捉能力。此外由圖21 和圖22 可知，大部分區域使用CPR 計算，保持了高計算效率。

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

圖20 CPR-HRFD2 壓力計算結果，自由度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]。該問題包含了一個運動的旋渦和一道靜止正激波，可以發展出具有平滑特征和間斷的復雜流動結構。如圖23 所示，該問題的計算域 (x,y)∈[0,2]×[0,1]。初始條 件 由 位 于 (0.5,y)的 靜 態 激 波Mas和 中 心 位 于(xc,yc)=(0.25,0.5)的 復合旋渦Mav給出，初始條件的具體設置參考文獻[25]。

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

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

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

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

圖25 CPR-HRFD2 密度計算結果，自由度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

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

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

表4 參考線Table 4 Reference lines

為對比CPR-HRFD2 和HRFD2 的激波捕捉能力和分辨率特性，在此選取①②③參考線進行分析。

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

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

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

3 結 論

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

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

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

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

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