一種分裂形式CPR 格式在欠解析流動(dòng)中的穩(wěn)定性研究

賈斐然，朱華君，嚴(yán) 紅，燕振國(guó)，*，石京昶

(1.西北工業(yè)大學(xué) 動(dòng)力與能源學(xué)院，西安 710000；2.空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000)

0 引言

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

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

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

Gassner 等[13]研究了高階DG 離散在欠解析湍流模擬中的精度，發(fā)現(xiàn)低階近似顯示出難以接受的數(shù)值離散誤差，而高階離散受混淆誤差的影響，往往是不穩(wěn)定的。對(duì)流項(xiàng)中的非線性項(xiàng)使得相對(duì)較低波數(shù)的模態(tài)相互非線性疊加，產(chǎn)生較高波數(shù)的模態(tài)，甚至是當(dāng)前基函數(shù)無(wú)法描述的模態(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）框架，并對(duì)無(wú)黏Taylor-Green 渦（TGV）算例進(jìn)行模擬，結(jié)果表明，與標(biāo)準(zhǔn)格式相比，這種新的分裂形式在求解欠解析湍流渦主導(dǎo)的流動(dòng)時(shí)增強(qiáng)了模擬的穩(wěn)定性。Winters 等[15]研究高階DG 方法在欠解析湍流計(jì)算中的精度和穩(wěn)定性，考慮無(wú)黏TGV 來分析DG 方法在極高雷諾數(shù)下進(jìn)行隱式大渦模擬（implicit large eddy simulation,ILES）的能力。為了抑制混淆誤差，采用過積分[13]和分裂形式兩種方法對(duì)控制方程進(jìn)行離散，并比較了這兩種方法在無(wú)黏TGV 算例中的表現(xiàn)，結(jié)果表明，分裂形式具有更好的穩(wěn)定性。Chan 等[16]將高階熵穩(wěn)定DG 格式用于求解欠解析流動(dòng)，這種采用LG（Legendre-Gauss）點(diǎn)并經(jīng)過熵投影的格式具有很好的魯棒性，與采用LGL（Legendre-Gauss-Lobatto）點(diǎn)的高階熵穩(wěn)定DGSEM 格式[12]相比，能夠分辨尺度更小的流動(dòng)結(jié)構(gòu)。

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

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

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

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

1.1 控制方程

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

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

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

式中，μ為動(dòng)力黏度。

1.2 高階CPR 格式

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 分裂形式CPR 格式

同樣以一維Euler 方程為例（計(jì)算坐標(biāo)下），分裂形式[10]是將方程中的對(duì)流項(xiàng)分裂成守恒形式和非守恒形式的組合，其一般形式為：

式中，P為壓力項(xiàng)，將單獨(dú)處理，不會(huì)被分裂。

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

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

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

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

式（23）的第一項(xiàng)可寫為：

因此，解點(diǎn)處分裂形式的重構(gòu)通量散度項(xiàng)計(jì)算如下：

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

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

1.4 子單元限制和激波偵測(cè)器

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

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

2 數(shù)值測(cè)試

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

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

2.1 離散守恒律及精度測(cè)試

2.1.1 Sod 激波管

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

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

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

圖1 T=0.2時(shí)Sod 激波管問題的密度分布.修正采用LG 點(diǎn)、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 對(duì)流等熵渦

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

計(jì)算域 [-10,10]×[-10,10]，邊界均為周期邊界，CFL=0.1，計(jì)算時(shí)間T=0.1。

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

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

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

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

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

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

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

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

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

2.2 穩(wěn)定性測(cè)試

2.2.1 欠解析旋渦流動(dòng)

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

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

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

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

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

圖2 基于Roe 通量，不同對(duì)流形式和解點(diǎn)下的渦量場(chǎng)比較：LGL 點(diǎn)，散度形式；LGL 點(diǎn)，F(xiàn)e 分裂形式；LG 點(diǎn)，散度形式；LG 點(diǎn)，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 通量下的渦量計(jì)算結(jié)果。值得注意的是，散度形式下LGL 點(diǎn)和LG 點(diǎn)分別在T=15.66和T=56.84時(shí)崩潰，因此沒有畫出其渦量分布。這是由于LLF 通量在靠近出口邊界處會(huì)產(chǎn)生偽反射（圖中被紅圈圈出的），這些偽反射與傳入的旋渦非線性相互作用，導(dǎo)致小尺度噪聲，使得計(jì)算不穩(wěn)定。此外，即使存在明顯的小尺度噪聲，分裂形式均能穩(wěn)定計(jì)算到T=150。

圖3 基于LLF 通量和Fe 分裂形式、LGL 和LG 解點(diǎn)下的渦量場(chǎng)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 通量下的渦量計(jì)算結(jié)果。值得注意的是，散度形式下LGL 點(diǎn)和LG 點(diǎn)分別在T=15.56和T=50.23時(shí)崩潰，因此沒有畫出其渦量分布。兩種分裂形式均能穩(wěn)定計(jì)算到T=150。

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

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

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

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

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

圖5 不同時(shí)刻黏性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 多項(xiàng)式階數(shù)的影響

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

圖6 總自由度約為 643時(shí)，采用Div-LG-g2-Roe 方案的動(dòng)能、動(dòng)能耗散率、擬渦能計(jì)算結(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 對(duì)流項(xiàng)形式的影響

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

圖7 總自由度約為 643時(shí)，采用SplitFe-LG-g2-Roe 方案的動(dòng)能、動(dòng)能耗散率、擬渦能計(jì)算結(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 點(diǎn)的散度形式與分裂形式的對(duì)比，與LG 點(diǎn)所得結(jié)論一致，分裂形式極大地提高了格式的穩(wěn)定性。

圖8 總自由度約為 643時(shí)，采用Div-LGL-g2-Roe 方案的動(dòng)能、動(dòng)能耗散率、擬渦能計(jì)算結(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時(shí)，采用SplitFe-LGL-g2-Roe 方案的動(dòng)能、動(dòng)能耗散率、擬渦能計(jì)算結(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 解點(diǎn)類型的影響

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

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

2.3 激波捕捉能力測(cè)試

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

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

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

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

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

3 結(jié)論

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

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

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

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

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

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

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