賈斐然，朱華君，嚴 紅，燕振國，*，石京昶

(1.西北工業(yè)大學(xué) 動力與能源學(xué)院，西安 710000；2.空氣動力學(xué)國家重點實驗室，綿陽 621000)

0 引言

計算流體力學(xué)（Computational Fluid Dynamics,CFD）是一種通過數(shù)值手段研究流體力學(xué)問題的方法，是對理論研究和實驗研究的補充。近些年來，為了對流動問題進行精細模擬，高階格式因其數(shù)值耗散誤差和色散誤差小的特性，逐漸成為CFD 中的一個重要研究方向[1]。

重構(gòu)修正過程（correction procedure via reconstruction,CPR）方法的思想最早由Huynh[2]提出，當(dāng)時稱之為通量重構(gòu)（flux reconstruction,FR）方法。Huynh 將其用于求解結(jié)構(gòu)網(wǎng)格上的雙曲守恒律方程，王志堅等[3-4]將此方法進行了推廣，提出適用于非結(jié)構(gòu)網(wǎng)格的提升配點懲罰法（lifting collocation penalty formulation,LCP），這兩種方法聯(lián)系緊密，被統(tǒng)稱為CPR 方法[5]。CPR 方法的優(yōu)勢主要在于三個方面：第一，在求解線性標量方程時，通過調(diào)整修正函數(shù)類型，CPR 方法可以恢復(fù)成間斷伽遼金（discontinuous Galerkin,DG）[6]、譜差分（spectral difference,SD）等方法[7]，但不同于DG 方法在求解過程中涉及非線性項的積分，CPR 方法是差分型方法，計算量小且計算復(fù)雜度低；第二，CPR 方法對復(fù)雜邊界有很好的適應(yīng)性，便于推廣到非結(jié)構(gòu)網(wǎng)格[8]；第三，該方法具有良好的緊致性，便于進行并行計算[9]。

CPR 方法或DG 方法在求解非線性問題時，即使流動中沒有出現(xiàn)像激波這樣的不連續(xù)性情況，也可能會由于離散非線性對流項引起的混淆誤差的累積而引發(fā)穩(wěn)定性問題[10]。這類問題在欠解析流動中尤為突出。為了增強欠解析流動模擬的穩(wěn)定性，目前主要存在的穩(wěn)定機制有濾波[6]、過積分[11]以及分裂形式[12]。

Gassner 等[13]研究了高階DG 離散在欠解析湍流模擬中的精度，發(fā)現(xiàn)低階近似顯示出難以接受的數(shù)值離散誤差，而高階離散受混淆誤差的影響，往往是不穩(wěn)定的。對流項中的非線性項使得相對較低波數(shù)的模態(tài)相互非線性疊加，產(chǎn)生較高波數(shù)的模態(tài)，甚至是當(dāng)前基函數(shù)無法描述的模態(tài)。高于Nyquist 波數(shù)的未分辨波數(shù)被“混淆”到分辨波數(shù)范圍。Blaisdell 等[14]在譜方法中提出一種分裂形式，發(fā)現(xiàn)其減少了高波數(shù)范圍內(nèi)的混淆誤差，Gassner 和Abe 等[10,12,15]將這種思路推廣到DG 和CPR 方法中。Gassner 等[12]構(gòu)造了一種高階分裂形式間斷伽遼金譜元法（discontinuous Galerkin spectral element method,DGSEM）框架，并對無黏Taylor-Green 渦（TGV）算例進行模擬，結(jié)果表明，與標準格式相比，這種新的分裂形式在求解欠解析湍流渦主導(dǎo)的流動時增強了模擬的穩(wěn)定性。Winters 等[15]研究高階DG 方法在欠解析湍流計算中的精度和穩(wěn)定性，考慮無黏TGV 來分析DG 方法在極高雷諾數(shù)下進行隱式大渦模擬（implicit large eddy simulation,ILES）的能力。為了抑制混淆誤差，采用過積分[13]和分裂形式兩種方法對控制方程進行離散，并比較了這兩種方法在無黏TGV 算例中的表現(xiàn)，結(jié)果表明，分裂形式具有更好的穩(wěn)定性。Chan 等[16]將高階熵穩(wěn)定DG 格式用于求解欠解析流動，這種采用LG（Legendre-Gauss）點并經(jīng)過熵投影的格式具有很好的魯棒性，與采用LGL（Legendre-Gauss-Lobatto）點的高階熵穩(wěn)定DGSEM 格式[12]相比，能夠分辨尺度更小的流動結(jié)構(gòu)。

CPR 方法在欠解析流動中的研究相比于DG 方法較少。Ball 和Jameson 研究了兩種新的高階FR 格式[17-18]即優(yōu)化能量穩(wěn)定通量重構(gòu)（optimized energy stable flux reconstruction,OESFR）和優(yōu)化通量重構(gòu)（optimized flux reconstruction,OFR）在粗網(wǎng)格上求解黏性TGV 的能力，并與能量穩(wěn)定通量重構(gòu)（energy stable flux reconstruction,ESFR）方法進行了比較。結(jié)果表明，OFR 格式精度最高，所計算的能譜與參考能譜吻合得最好，且在高波數(shù)下優(yōu)于ESFR 方法的能譜。Abe 等[10]針對可壓縮Euler 和Navier-Stokes（NS）方程，構(gòu)造并證明了一種穩(wěn)定且守恒的FR 格式。這種格式采用LGL 解點，對對流項采用分裂形式，基于多項式分析嚴格推導(dǎo)了滿足離散守恒和動能保持性質(zhì)的充分條件，通過黏性TGV 算例證明：基于這種分裂形式的FR 框架，使用無耗散動能保持（kinetic energy preservation,KEP）通量[19]，在非常高階（p15）和相對粗糙網(wǎng)格下的模擬仍然是穩(wěn)定的。Abe 提出了一種基于LG 點的滿足離散守恒律的分裂形式CPR 格式，但并未研究其求解欠解析流動的特性。

Ramírez 等[20]提出了一種通用的子單元限制策略來構(gòu)造魯棒的高精度節(jié)點DG 格式，主要思路是構(gòu)造兼容的低階有限體積（FV）離散，并與高階DG 進行凸混合。這種策略可以有效處理強激波，在KPP 問題[21]、湍流和高超聲速Euler 模擬，及以激波和湍流為特征的MHD 問題上有很好表現(xiàn)。Zhu 等[22–24]針對基于LG 點的CPR 方法，提出了一類先驗子單元限制格式。首先，利用激波偵測器檢測存在間斷的問題單元；然后，將問題單元分解為非均勻子單元，每個子單元包含一個解點；最后，對光滑單元使用CPR 方法計算，對問題單元使用緊致非均勻非線性加權(quán)（compact nonuniform nonlinear weighted,CNNW）插值格式計算。這種新策略在分辨率和激波捕捉魯棒性之間具有良好的平衡。

分裂形式CPR 格式在欠解析流動中的研究較少，且主要以LGL 點作為解點。本文研究了以LG 點為解點的分裂形式CPR 格式在欠解析流動中的應(yīng)用。首先，對比了基于LG 點的分裂形式和散度形式在求解欠解析流動時的穩(wěn)定性，展現(xiàn)了分裂形式在增強格式穩(wěn)定性方面的優(yōu)勢；其次，從分辨率和穩(wěn)定性的角度比較了LG 點和LGL 點分裂形式CPR 格式；最后，首次將基于LG 點的分裂形式CPR 格式與Zhu 等提出的CNNW 子單元限制格式相結(jié)合，求解含激波的Kelvin-Helmholtz 不穩(wěn)定性問題。

1 控制方程與數(shù)值方法

1.1 控制方程

考慮守恒形式的三維N-S 方程：

式中，U為守恒變量，F(xiàn)、G、H分別為x、y、z方向的無黏通量，F(xiàn)v、Gv、Hv分別為x、y、z方向的黏性通量，具體形式如下：

式中，ρ 為密度，u、v、w分別為x、y、z方向的速度，p為壓力，E為單位質(zhì)量流體的總內(nèi)能。本文只涉及理想氣體，比熱比 γ=1.4。黏性應(yīng)力滿足：

式中，μ為動力黏度。

1.2 高階CPR 格式

本小節(jié)以一維守恒形式Euler 方程為例，介紹CPR 格式的構(gòu)造方法。守恒方程形式為：

計算域 [a,b]，首先將其剖分為N個區(qū)間，將其中第j個區(qū)間設(shè)為Ij=[xj-1/2,xj+1/2],j=1,·,N。為 了便于討論，將每個區(qū)間映射到標準單元I=[-1,1]上。假設(shè)每個區(qū)間Ij上的點x到標準單元I的點 ξ存在一個線性映射關(guān)系：

式中，xj=(xj-1/2+xj+1/2)/2 是區(qū)間Ij的中心，hj=xj+1/2-xj-1/2為區(qū)間Ij的長度。

在區(qū)間Ij上定義K個解點xj,k，k=1,·,K，解點上對應(yīng)的守恒變量為Uj,k，k=1,·,K。為了便于分析，本文對所有區(qū)間均采用相同類型的解點，即LG 點或LGL 點。設(shè)標準單元上的解點位置為 ξk，k=1,·,K，則區(qū)間Ij上的解點位置由下式計算得到：

由解點處的守恒變量Uj,k，k=1,·,K，通過構(gòu)造K-1 次Lagrange 插值多項式來逼近Uj：

其中l(wèi)k(ξ)是Lagrange 插值基函數(shù)，形式如下：

同理，區(qū)間Ij上的通量函數(shù)由下式逼近：

通量多項式Fj(ξ)的構(gòu)造只依賴于區(qū)間內(nèi)部的解點通量值，在相鄰區(qū)間交界面處的通量值不連續(xù)，即Fj(1)≠Fj+1(-1)。因此，需要構(gòu)造一個連續(xù)通量多項式表達式為：

式中，gL、gR為邊界的左右修正函數(shù)，滿足：

修正函數(shù)的形式一般有以下三種：

式中，PK是K階Legendre 多項式。gDG、gGa和g2分別是可將CPR 格式恢復(fù)成DG 格式、SD 格式和Huynh 型格式的修正函數(shù)[2]。

為了保證相鄰區(qū)間交界面處的連續(xù)性，需要引入Riemann 通量常用的Riemann 通量有局部Lax-Friedrichs（LLF）、Roe 通量[25]等。

最后，得到Euler 方程的半離散形式：

針對上式，本文使用三階TVD Runge-Kutta 格式[26-27]進行時間離散。以上為一維CPR 格式的實現(xiàn)過程，二維及三維CPR 格式可以直接由一維形式擴展而來[10,22]。

1.3 分裂形式CPR 格式

同樣以一維Euler 方程為例（計算坐標下），分裂形式[10]是將方程中的對流項分裂成守恒形式和非守恒形式的組合，其一般形式為：

式中，P為壓力項，將單獨處理，不會被分裂。

常用的分裂形式類型有Fe 分裂[28]、Bl 分裂[14]、KG1 分裂[29]和KG2 分裂[29]等。本文只考慮第一種分裂形式，表達式為：

其中，?={1,u,H}T。式（21）中的 (Split)C可由式（22）代替。

(Split)C表示分裂形式下重構(gòu)通量（通量散度）的ξ方向?qū)?shù)。首先，使用解點處的守恒變量值計算解點處分裂形式的通量散度項。在第n個單元上，假設(shè)符號函數(shù) {fn,gn,hn} 代表 {ρ,u,?}，則式（22）中等號右端項在解點處可由以下形式表示：

式中，引入符號ISP;K[ψn] 表示一個任意的函數(shù) ψn的K階多項式近似：

式（23）的第一項可寫為：

因此，解點處分裂形式的重構(gòu)通量散度項計算如下：

Abe[10]證明了以LGL 點作為解點時，F(xiàn)e 分裂形式滿足離散守恒律。同時也指出：如果以LG 點作為解點，要使分裂形式滿足離散守恒律，需要對邊界通量項進行修正，形式如下：

式中I[?] 即為ISP;K[?]。

1.4 子單元限制和激波偵測器

本文所使用的子單元限制策略的主要思想是，利用激波偵測器偵測出流場中存在的大梯度或間斷的單元，并將其標記為問題單元，然后將這些問題單元劃分為非等距子單元，每個子單元包含一個解點，最后在子單元上添加限制手段，具體限制方式在文獻[24]中給出。光滑單元用分裂形式CPR 格式計算，而問題單元則使用二階CNNW 格式計算以捕捉激波。值得注意的是，基于LG 點的分裂形式與子單元限制結(jié)合時并不直接滿足離散守恒律，仍然需要對分裂形式使用邊界通量修正，在2.1.2 小節(jié)對此進行了數(shù)值驗證。

本文采用的激波偵測器為最高模態(tài)衰減（highest modal decay,MDH）。該方法利用解多項式的最高模態(tài)能量系數(shù)在光滑區(qū)域衰減較快、非光滑區(qū)域衰減較慢的特點[30]，通過設(shè)置閾值來判斷單元內(nèi)是否存在大梯度或間斷。為了避免奇偶效應(yīng)，Hennemann 等[31]使用最高和次高模態(tài)系數(shù)改進了此方法，本文采用這一改進方法。

2 數(shù)值測試

使用一維Sod 激波管、二維對流等熵渦[22]問題來測試格式的離散守恒律和數(shù)值精度；使用二維欠解析旋渦流動[32-33]和三維黏性TGV[10,18]算例來測試格式的穩(wěn)定性；使用二維無黏Kelvin-Helmholtz 不穩(wěn)定性算例來測試分裂形式結(jié)合子單元限制技術(shù)的激波捕捉能力。為了便于描述，所測試的計算方案按以下方式命名：對流項-解點類型-修正函數(shù)-無黏通量。例如，Div-LG-gDG-Roe 表示對流項采用散度形式、解點選擇LG 點、修正函數(shù)選擇gDG、無黏通量選擇Roe 通量。表1 列出了本文所采用的算例及其對應(yīng)的計算條件設(shè)置。此外，對于黏性通量，均使用BR2 格式[34]。

表1 算例類型及對應(yīng)的計算條件設(shè)置Table 1 Simulation cases and the corresponding calculation condition settings

2.1 離散守恒律及精度測試

2.1.1 Sod 激波管

Sod 激波管問題計算域為 0 ≤x≤1，被劃分成等距的200 個網(wǎng)格。初始條件如下：

該算例存在精確解[35]。左右邊界施加Dirichlet邊界條件，計算時間T=0.2，CFL=0.1。解多項式階為p2（即三階空間精度），總自由度為600。為了便于評估數(shù)值穩(wěn)定性，不采用任何激波捕捉方法。

圖1 是使用LG 點和修正函數(shù)g2的計算結(jié)果。圖1（a）表明：基于LG 點的分裂形式CPR 格式直接求解Sod 激波管問題時，并不能正確計算激波位置，即不滿足離散守恒律。圖1（b）是經(jīng)過邊界通量修正[10]后的計算結(jié)果，此時散度形式和Fe 分裂形式均能滿足離散守恒律。

圖1 T=0.2時Sod 激波管問題的密度分布.修正采用LG 點、g2 修正函數(shù)和Roe 通量，總自由度為600(p2)Fig.1 Density profiles of the Sod shock tube problem at T=0.2:(a) without boundary flux fix;(b) with boundary flux fix.LG points,g2 correction function and Roe flux are used,and the total number of DoFs is 600(p2)

2.1.2 對流等熵渦

對流等熵渦問題除了可以用來驗證格式的離散守恒律[22]外，還可以用于測試精度。該算例是在一個均勻流動的基礎(chǔ)上添加一個等熵旋渦擾動。均勻流動設(shè)置如下：

計算域 [-10,10]×[-10,10]，邊界均為周期邊界，CFL=0.1，計算時間T=0.1。

首先用該算例測試離散守恒律。全局離散守恒律誤差表達式如下：

其中，q可以用于表示 ρ、ρu等守恒變量，下標“0”表示初始時刻的守恒變量。

表2 給出了對流項的不同形式在使用LG 點時的離散守恒律誤差，其中Div-subcell 和SplitFe-subcell分別表示結(jié)合子單元限制的散度形式和分裂形式CPR 格式。當(dāng)離散守恒律誤差接近機器零時，認為格式滿足離散守恒律。由表2 可知，當(dāng)分裂形式與子單元限制相結(jié)合時，也需要經(jīng)過邊界通量修正才能滿足離散守恒律。

表2 對流等熵渦的全局離散守恒律誤差Table 2 Errors of the global discrete conservation law in the computation of convecting isentropic vortex

表3 和表4 分別是基于LG 點的散度形式和Fe 分裂形式的精度測試（p4）結(jié)果，這兩種對流項形式計算的誤差與收斂階基本相同。表5 和表6 則是基于LGL 點的散度形式和Fe 分裂形式的精度測試（p4）結(jié)果，F(xiàn)e 分裂形式的L1誤差小于散度形式，而L∞誤差大于散度形式。最后對比表4 和表6，發(fā)現(xiàn)基于LG 點的格式誤差明顯小于LGL 點（接近一個量級）。

表3 基于LG 點的散度形式精度測試(p4)Table 3 Accuracy test of divergence form based on LG points (p4)

表4 基于LG 點的Fe 分裂形式精度測試(p4)Table 4 Accuracy test of the Fe split form based on LG points (p4)

表5 基于LGL 點的散度形式精度測試(p4)Table 5 Accuracy test of the divergence form based on LGL points (p4)

表6 基于LGL 點的Fe 分裂形式精度測試(p4)Table 6 Accuracy test of the Fe split form based on LGL points (p4)

2.2 穩(wěn)定性測試

2.2.1 欠解析旋渦流動

為了評估格式對實際流動的欠解析模擬的影響，使用二維可壓縮Euler 方程和非定常流入邊界條件，模擬渦流沿下游傳播的被動發(fā)生器。該問題是無黏的，使用Euler 方程來研究數(shù)值方法在黏性消失極限（極高雷諾數(shù)）下的性質(zhì)和行為[15,36]。算例設(shè)置如下：

計算域為 [0,20π]×[-π,π]，計算網(wǎng)格數(shù) 120×12，CFL=0.5，T=150.0 。在y=±π位置施加壁面邊界條件，x=0位置施加流入邊界條件：

式中，ρ∞=1、u∞=1是自由流密度和均勻流速度，是自由流靜壓，用于通過聲速c∞=u∞Ma-1定義流動的參考馬赫數(shù)。此外，入口處擾動參數(shù)定義為A=1/2、K=5、?=1，出口邊界為初始條件與出口邊界相同。T=20π左右時，流場將達到漸近狀態(tài)。

本算例使用的多項式階為p5。主要考慮三種CPR格式中常見的Riemann 求解器，即HLL 通量[37-38]、Roe 通量和LLF 通量。在本算例中，在多項式高階和欠解析條件下，Riemann 求解器會有不同的特性。本算例中選擇馬赫數(shù)為0.03，因為在低馬赫數(shù)下，會放大HLL 和LLF 通量產(chǎn)生的偽反射和數(shù)值噪聲，從而更好地展現(xiàn)數(shù)值不穩(wěn)定性帶來的影響[32]。除了使用多種Riemann 求解器外，還使用兩種不同對流項形式的CPR 格式（散度形式和Fe 分裂）和兩種解點類型（LG 點和LGL 點）。為了評估這些方面的影響，考慮流場中的渦量分布

圖2 給出了在Roe 通量下，采用不同對流項形式和不同解點的渦量計算結(jié)果。所有情況均能穩(wěn)定計算至結(jié)束。比較圖2 中的第一和第三小圖，發(fā)現(xiàn)采用LG 點對流場中的微小結(jié)構(gòu)捕捉得更清楚。

圖2 基于Roe 通量，不同對流形式和解點下的渦量場比較：LGL 點，散度形式；LGL 點，F(xiàn)e 分裂形式；LG 點，散度形式；LG 點，F(xiàn)e 分裂形式Fig.2 Comparison of the vorticity fields computed with different convention from and solution point combination based on the Roe flux:LGL points,Div form;LGL points,SplitFe form;LG points,Div form;LG points,SplitFe form

圖3 是在LLF 通量下的渦量計算結(jié)果。值得注意的是，散度形式下LGL 點和LG 點分別在T=15.66和T=56.84時崩潰，因此沒有畫出其渦量分布。這是由于LLF 通量在靠近出口邊界處會產(chǎn)生偽反射（圖中被紅圈圈出的），這些偽反射與傳入的旋渦非線性相互作用，導(dǎo)致小尺度噪聲，使得計算不穩(wěn)定。此外，即使存在明顯的小尺度噪聲，分裂形式均能穩(wěn)定計算到T=150。

圖3 基于LLF 通量和Fe 分裂形式、LGL 和LG 解點下的渦量場Fig.3 Comparison of the vorticity fields computed with the LGL and LG solution points based on the LLF flux and Fe split form

圖4 是在HLL 通量下的渦量計算結(jié)果。值得注意的是，散度形式下LGL 點和LG 點分別在T=15.56和T=50.23時崩潰，因此沒有畫出其渦量分布。兩種分裂形式均能穩(wěn)定計算到T=150。

圖4 基于HLL 通量和Fe 分裂形式、LGL 和LG 解點下的渦量場Fig.4 Comparison of the vorticity fields computed with the LGL and LG solution points based on the HLL flux and Fe split form

2.2.2 黏性Taylor-Green 渦

在本小節(jié)中，我們模擬了黏性TGV 問題。該問題通常用于研究旋渦動力學(xué)和湍流轉(zhuǎn)捩與衰減[39]。該問題由最初包含光滑渦量分布的立方體流體組成，能夠反映數(shù)值格式的魯棒性以及對多尺度流動結(jié)構(gòu)的模擬能力。在美國航空航天學(xué)會航空航天科學(xué)會議[40]上舉行的計算流體動力學(xué)高階方法國際研討會上，該算例用于評估湍流模擬方法。許多作者已經(jīng)成功使用高階模擬方法預(yù)測該流場[13,18,41-42]。本文將使用TGV 來比較散度形式和分裂形式在湍流欠解析模擬中的精度和穩(wěn)定性。初始條件設(shè)置如下：

其中，馬赫數(shù)和雷諾數(shù)分別設(shè)為0.1 和1 600，L為參考長度。計算域是 -πL≤x,y,z≤πL，被劃分為互不重疊的均勻六面體單元。在本測試中，包括內(nèi)部解點在內(nèi)的自由度總數(shù)固定在 643左右，以便在不同階的解多項式之間進行比較。表7 列出了網(wǎng)格數(shù)與解多項式階的對應(yīng)關(guān)系。計算時間為 0 ≤T≤20，使用固定時間步長 ?t=3.14×10-4。本算例只考慮了Fe 分裂形式。為了更好地描述該問題隨時間的演化過程，需要定義以下幾個參數(shù)：

表7 總單元數(shù)和解多項式階數(shù)Table 7 Total number of cells and order of the solution polynomial

式中，u為速度矢量，ω為渦量矢量。

圖5 分別給出了不同時刻下，多項式階為p3 的方案SplitFe-LG-g2-Roe 的z方向渦量等值面圖。從T=0開始，初始流場存在明顯的大渦結(jié)構(gòu)，隨后不斷地拉伸和旋轉(zhuǎn)，逐漸破裂成較小的渦結(jié)構(gòu)[43]。

圖5 不同時刻黏性Taylor-Green 渦 z 方向渦量等值面圖Fig.5 Iso-surfaces of vorticity in z-direction for the viscous Taylor-Green vortex at different time instances

2.2.2.1 多項式階數(shù)的影響

為了節(jié)省計算資源和比較格式的穩(wěn)定性，本算例的計算均是在欠解析即網(wǎng)格數(shù)量不足條件下進行的。我們從p1 開始，逐漸提高多項式階數(shù)，直到出現(xiàn)不穩(wěn)定解。如圖6 所示，可以明顯看出：在保證總自由度不變時，提高多項式階數(shù)使得數(shù)值解逼近參考解，這與COX 之前的結(jié)論一致[44]。本文主要關(guān)注的是欠解析情況下格式的非線性穩(wěn)定性問題。當(dāng)多項式階數(shù)提高時，在T=5 時流場很容易由于混淆誤差而崩潰，這時就需要具有去混淆效果的分裂形式來增強格式的穩(wěn)定性。

圖6 總自由度約為 643時，采用Div-LG-g2-Roe 方案的動能、動能耗散率、擬渦能計算結(jié)果Fig.6 Computational results of the Div-LG-g2-Roe scheme with the total degree of freedom of 643: (a) kinetic energy;(b) kinetic energy dissipation rate;(c) enstrophy

2.2.2.2 對流項形式的影響

圖6 和圖7 對比了方案Div-LG-g2-Roe 和方案SplitFe-LG-g2-Roe。在圖6 中，多項式階從p1 到p4 的方案均能穩(wěn)定計算至T=20，而當(dāng)階數(shù)提高到很高階（p7）時，方案在T=5左右崩潰。而圖7 中，所有多項式階的方案均穩(wěn)定。這一結(jié)果表明，即使是很高階（p8）離散，分裂格式CPR 格式仍能保證數(shù)值穩(wěn)定性。

圖7 總自由度約為 643時，采用SplitFe-LG-g2-Roe 方案的動能、動能耗散率、擬渦能計算結(jié)果Fig.7 Computational results of the SplitFe-LG-g2-Roe scheme with the total degree of freedom of 643: (a) kinetic energy;(b) kinetic energy dissipation rate;(c) enstrophy

圖8 和圖9 是基于LGL 點的散度形式與分裂形式的對比，與LG 點所得結(jié)論一致，分裂形式極大地提高了格式的穩(wěn)定性。

圖8 總自由度約為 643時，采用Div-LGL-g2-Roe 方案的動能、動能耗散率、擬渦能計算結(jié)果Fig.8 Computational results of the Div-LGL-g2-Roe scheme with the total degree of freedom of 643: (a) kinetic energy;(b) kinetic energy dissipation rate;(c) enstrophy

圖9 總自由度約為 643時，采用SplitFe-LGL-g2-Roe 方案的動能、動能耗散率、擬渦能計算結(jié)果Fig.9 Computational results of the SplitFe-LGL-g2-Roe scheme with a total degree of freedom of 643: (a) kinetic energy;(b) kinetic energy dissipation rate;(c) enstrophy

2.2.2.3 解點類型的影響

如圖6 和圖8 所示，首先對比散度形式下LG 點和LGL 點的計算結(jié)果。LG 點在p1～p4 方案下穩(wěn)定，而LGL 點僅能在p1～p2 方案下穩(wěn)定。因此，這一結(jié)果表明，在多項式階不是很高的條件下（p5 以下），基于LG 點比基于LGL 點的散度形式CPR 格式穩(wěn)定性更優(yōu)。

圖7 和圖9 是分裂形式下LG 點和LGL 點的計算結(jié)果。兩種解點類型均能在所要求的多項式階下保證穩(wěn)定。從圖7（c）和圖9（c）可以看出，LG 點相較于LGL 點的優(yōu)勢在于：其計算精度更高，尤其是多項式階為p7 時，LG 點的擬渦能與參考解吻合得更好。

2.3 激波捕捉能力測試

本小節(jié)考慮二維無黏Kelvin-Helmholtz 不穩(wěn)定性問題[45]。該算例對于CPR 方法來說很具有挑戰(zhàn)性，因為它包含嚴重欠解析的渦結(jié)構(gòu)，馬赫數(shù)約等于0.6。本測試的目的是將結(jié)合子單元限制技術(shù)的分裂形式CPR 格式應(yīng)用于求解欠解析流動。初始條件由下式給出：

其中B=tanh(15y+7.5)-tanh(15y-7.5) 。計算域[-1,1]2且均為周期邊界。我們使用 64×64個四邊形單元對計算域進行均勻劃分，多項式階數(shù)為p7，計算時間T=10。

在圖10 中，我們分別繪制了T=3.7 和T=10時分裂形式結(jié)合子單元限制求解Kelvin-Helmholtz 不穩(wěn)定性問題的密度等值線。在使用相同激波偵測器和相同自由度的情況下，相比于Ramírez 等[20]提出的基于LGL 點、結(jié)合子單元限制技術(shù)的DGSEM 格式的計算結(jié)果，這種新策略能分辨更多的小尺度特征且振蕩更少。圖11 分別給出了兩個時刻下的問題單元分布。

圖10 不同時刻下Kelvin Helmholtz 不穩(wěn)定性模擬的密度云圖Fig.10 Density contours for the Kelvin Helmholtz instability simulation at different time instances

圖11 不同時刻下Kelvin Helmholtz 不穩(wěn)定性模擬的問題單元分布Fig.11 Distribution of problematic cells for the Kelvin Helmholtz instability simulation at different time instances

3 結(jié)論

本文在CPR 格式的基礎(chǔ)上，比較了采用不同解點類型的分裂形式對求解欠解析流動時的穩(wěn)定性和精度的影響，主要有以下結(jié)論：

1）Sod 激波管和等熵渦算例結(jié)果表明，以LG 點作為解點時，通過邊界通量修正的Fe 分裂形式滿足離散守恒律。

2）精度測試結(jié)果表明，使用LG 點的散度形式和分裂形式的L1、L∞誤差基本相同；在相同散度形式或分裂形式下，使用LG 點的誤差明顯小于LGL 點（約一個量級）。

3）二維欠解析旋渦流動算例用來考察數(shù)值通量的特性和格式的穩(wěn)定性。LLF 和HLL 通量在出口邊界處會產(chǎn)生偽反射和小尺度噪聲，導(dǎo)致散度形式CPR 格式在計算時崩潰，而Fe 分裂形式能夠保證穩(wěn)定計算。Roe 通量則能完全抑制這種偽反射。此外，使用LGL 點時放大了出口處產(chǎn)生的偽反射，而LG 點對流場中的微小結(jié)構(gòu)捕捉得更清楚。

4）三維黏性Taylor-Green 渦算例結(jié)果表明，無論采用LG 點或LGL 點，即使是很高階（p7）的離散，分裂形式CPR 格式對于保證計算穩(wěn)定依然是有效的，且LG 點精度更高，在p7 下與參考解匹配得更好。如果統(tǒng)一使用散度形式，在多項式階不是很高的條件下（p5 以下），采用LG 點的格式穩(wěn)定性優(yōu)于LGL 點。

5）二維無黏Kelvin-Helmholtz 不穩(wěn)定性算例考察了分裂形式結(jié)合子單元限制技術(shù)的激波捕捉能力，相比于基于LGL 點的間斷譜元法子單元限制策略，前者在分辨率和穩(wěn)定性方面均有優(yōu)勢。

總的來說，本文將基于LG 點的分裂形式CPR 格式用于求解欠解析流動問題，不僅能提高格式的穩(wěn)定性，而且相比于LGL 點，該方法對流場的分辨率更高；在結(jié)合子單元限制策略后，其對含激波的欠解析流動也能得到很好的解。本文基于CPR 方法、采用“分裂+LG”策略時需要采用邊界通量修正保證守恒性，如果想將這一策略應(yīng)用于其他高階譜元方法，在如何滿足格式的守恒性方面有待進一步探索。