一種新的可計(jì)算可壓縮流動(dòng)的預(yù)處理方法*

劉博 邢樸 丁松 謝明軍 馮林 時(shí)曉天?

1) (中國(guó)航天空氣動(dòng)力技術(shù)研究院,北京 100074)

2) (北京航空航天大學(xué)航空科學(xué)與工程學(xué)院,北京 100083)

3) (北京航空航天大學(xué)大型飛機(jī)高級(jí)人才培訓(xùn)班,北京 100083)

4) (南開大學(xué)數(shù)學(xué)科學(xué)學(xué)院,天津 300071)

針對(duì)使用可壓縮流動(dòng)數(shù)值方法求解不可壓縮流動(dòng)存在的剛性問題,基于虛擬壓縮法思想,構(gòu)造了一種以Mach 數(shù)、速度、密度、溫度等變量為元素的預(yù)處理矩陣,改變了控制方程組的特征根并使其量級(jí)更接近.通過理論推導(dǎo)與分析,證明新方法相比Weiss,Pletcher,Dailey 和Choi 的方法而言,不僅能降低方程組的剛性,提高了數(shù)值求解效率,而且擁有更好的穩(wěn)定性,此外還能實(shí)現(xiàn)低速流動(dòng)和高速流動(dòng)之間的光滑過渡.采用有限差分格式進(jìn)行離散,對(duì)流項(xiàng)的Roe 格式作為基本加權(quán)無振蕩(WENO)格式的求解器,黏性項(xiàng)則使用中心型緊致差分格式來計(jì)算,與預(yù)處理矩陣相結(jié)合展開數(shù)值實(shí)驗(yàn),結(jié)果表明新預(yù)處理方法可以實(shí)現(xiàn)對(duì)無黏和有黏不可壓縮流動(dòng)問題的高精度模擬,且擁有比Weiss 和Pletcher 等提出的方法更好的收斂性和穩(wěn)定性.

1 引言

計(jì)算流體力學(xué)(computational fluid dynamics,CFD)是流體力學(xué)、計(jì)算數(shù)學(xué)與計(jì)算機(jī)技術(shù)相結(jié)合而發(fā)展起來的新興學(xué)科,它通過數(shù)值仿真和可視化處理,對(duì)流體流動(dòng)和熱傳導(dǎo)等物理現(xiàn)象進(jìn)行計(jì)算機(jī)數(shù)值分析和研究.其中有限差分法[1?3]是應(yīng)用最廣泛、最成的方法之一.

隨著CFD 的不斷發(fā)展,人們所求解的問題也日益復(fù)雜化,這些問題通常歸結(jié)為高速和低速流動(dòng)問題,高速流動(dòng)需采用可壓縮流動(dòng)計(jì)算方法[4,5],而低速流動(dòng)采用不可壓縮流動(dòng)計(jì)算方法[6?8].然而,這兩種方法的特性存在很大差別,例如,高超聲速往往伴隨著如激波一類的強(qiáng)間斷存在[9],亞聲速流動(dòng)的流場(chǎng)中任何微小的擾動(dòng)都將影響整個(gè)流場(chǎng)[10],因此具體數(shù)值方法要在特定環(huán)境下使用.針對(duì)高速的可壓縮流動(dòng)數(shù)值求解問題,學(xué)者們已研究出較成熟的理論,而面對(duì)低速流動(dòng)時(shí),無論采取顯格式還是隱格式計(jì)算,速度越小對(duì)計(jì)算誤差和收斂速度影響越大,這也被稱為剛性問題[11,12],這種剛性給數(shù)值方法帶來了很大的困難.

同時(shí),不可壓縮流動(dòng)數(shù)值方法在實(shí)際工程應(yīng)用又有很大的需求,如翼型繞流,邊界層內(nèi)流動(dòng)[13?15].在某些問題中,還會(huì)出現(xiàn)高低速流動(dòng)并存的現(xiàn)象[16?18].為了求解此類問題,需要發(fā)展能夠有效求解低速流動(dòng)的數(shù)值方法.為了求解有很強(qiáng)不可壓縮性的低速流動(dòng)問題,Fasel[19]提出渦函數(shù)-流函數(shù)法,Aziz[20]提出勢(shì)函數(shù)-流函數(shù)法,其本質(zhì)是引入流函數(shù)和渦函數(shù)并對(duì)控制方程進(jìn)行變換,但是壁面處渦量難以處理,而且推廣到三維問題后過于復(fù)雜.Halrow 和Wetch[21]提出格子中標(biāo)記點(diǎn)算法(marker in cell,MAC),Hirt 等[22]則在MAC 算法上加以改進(jìn)提出SOLA 算法,此類方法可以有效地處理含有復(fù)雜的自由面或界面的流動(dòng),但是對(duì)時(shí)間步限制非常嚴(yán)格,所以求解效率較低.Spalding 和Patankar[23]提出壓力聯(lián)系方程的半隱式法(semi implicit method for pressure linked equations,SIMPLE),該方法核心在于不斷修正壓力項(xiàng),通過修正后的壓力求解速度并滿足連續(xù)方程,直到收斂.為了提高收斂速度,Patankar[24]又提出SIMPLER (SIMPLE revised)算法,van Doormaal 和Raithby[25]提出SIM PLEC (SIMPLE consistent)算法,Minkowycz 等[26]提出SIMPLEX (SIMPLE extraplation)算法,Issa等[27]提出PISO (pressure implicit with splitting of operators)算法,他們都通過加強(qiáng)速度與壓力的耦合來達(dá)到加速收斂.但是這類方法占用內(nèi)存很大,插值過程繁瑣,大量時(shí)間都消耗在修正壓力項(xiàng)上,且編程復(fù)雜.

實(shí)際上,不可壓縮流體是為了簡(jiǎn)化某些問題提出的理想模型,任何流體都是可壓縮的.而且很多時(shí)候熱量的影響是不可忽略的[28],質(zhì)量、動(dòng)量和能量方程需要耦合在一起求解.鑒于可壓縮流動(dòng)的數(shù)值方法相對(duì)成熟,因此有研究者提出使用計(jì)算可壓縮流動(dòng)的數(shù)值方法來計(jì)算不可壓縮流動(dòng)問題.Chorin[29]于1967 年提出虛擬壓縮法就是基于這種思想,這也成為了預(yù)處理方法的雛形,該方法另一優(yōu)勢(shì)是可將高速與低速流動(dòng)統(tǒng)一在同一算法框架下求解,避免了切換數(shù)值方法的繁瑣和交換信息而產(chǎn)生的精度損失.最早的預(yù)處理矩陣由Turkle[30]于1986 年提出,其核心作用在于改變時(shí)間導(dǎo)數(shù)項(xiàng),減小方程組特征速度的量級(jí)差異,1993 年Choi 和Merkle[31]提出與Chorin-Turkel 矩陣相似的適用于Euler 方程的預(yù)處理矩陣.最早適用于全速域的預(yù)處理矩陣是由Viviand[32]設(shè)計(jì)的四參數(shù)矩陣,但在駐點(diǎn)附近剛性依然很強(qiáng).針對(duì)Navier-Stokes 方程,Choi 與Merkle[33]以及Buelow 等[34]提出了考慮Reynolds 數(shù)的影響的預(yù)處理矩陣,Godfrey 等[35]提出的矩陣將Euler 處理矩陣結(jié)合了離散方程的黏性/無黏性的塊Jacobi 矩陣,但仍然難以降低駐點(diǎn)附近的剛性.Weiss 等[36],Pletcher 和Chen[37]提出的矩陣形式簡(jiǎn)單,且降低了駐點(diǎn)附近剛性,被廣泛的應(yīng)用于FLUENT,OVERFLOW 等商業(yè)軟件.Dailey 和Pletcher[38]提另一種矩陣并將其與多重網(wǎng)格法結(jié)合.且預(yù)處理方法不僅可以與有限差分法搭配,也可以與有限體積法[39]、間斷有限元法[40]搭配,此外也可以與雙時(shí)間步技術(shù)結(jié)合用于計(jì)算非定常問題[41].

本文提出一種新的預(yù)處理矩陣,該矩陣引入了兩個(gè)關(guān)于Mach 數(shù)的函數(shù)b和ML,使得方程組的剛性進(jìn)一步降低,在確保(甚至能提高)精度和穩(wěn)定性前提下加速了收斂速度,并且在高速與低速流動(dòng)之間過渡更加光滑(其中ML在有黏和無黏情況下?lián)碛胁煌男问?.通過數(shù)值實(shí)驗(yàn),初步證明了本文方法的可靠性,并與其他方法相對(duì)比,體現(xiàn)本文方法的優(yōu)勢(shì).

2 基本概念

2.1 物理量的無量綱化

為了方便實(shí)現(xiàn)數(shù)值實(shí)驗(yàn)的相似模擬,本文對(duì)變量采取了如下無量綱化處理:

其中,x,y,z代表空間直角坐標(biāo)系的三個(gè)坐標(biāo)分量;L代表參考長(zhǎng)度;t代表時(shí)間;u,v,w分別代表三個(gè)坐標(biāo)軸方向上的速度;g代表重力加速度;ρ代表密度;E為單位體積流體的總能量;c代表聲速;T代表溫度;μ為黏性系數(shù);上標(biāo)“ * ”代表有量綱量,下標(biāo)“ ∞ ”代表參考值;下文2.2 節(jié)中的Re為來流Reynolds 數(shù),Reρ∞c∞L/μ∞.

2.2 流體力學(xué)控制方程組

Navier-Stokes 方程組是流體動(dòng)力學(xué)中最重要的基本方程,它是黏性流體微團(tuán)應(yīng)用牛頓第二定律得到的運(yùn)動(dòng)微分方程.以三維可壓縮流動(dòng)的情形為例,忽略外加熱源和徹體力的無量綱化方程:

其中

黏性應(yīng)力項(xiàng)與熱傳導(dǎo)流通量分別為

其中當(dāng)i=j時(shí),δij=1,如果i≠j,δij=0;Pr是Prandtl 數(shù);U代表速度矢量,對(duì)于無黏流動(dòng),(1)式的右端為零.在實(shí)際的應(yīng)用中,經(jīng)常會(huì)用到曲線坐標(biāo)系形式的方程組,具體如下:

這里 (ξ,η,ζ) 為曲線坐標(biāo)系.

本文采用的是方程組(2),并取τ=t,后文將省略該方程組中變量的上標(biāo)“?”.

3 方程組的剛性問題

3.1 預(yù)處理方法

使用與時(shí)間相關(guān)的可壓縮流動(dòng)數(shù)值方法來求解不可壓縮流動(dòng)時(shí),方程組的特征根量級(jí)差異過大(也稱剛性過大),這會(huì)嚴(yán)重影響計(jì)算的收斂性.基于虛擬壓縮法原理的預(yù)處理法可以有效地解決這一問題,即在時(shí)間導(dǎo)數(shù)項(xiàng)前乘上個(gè)預(yù)處理矩陣,改變方程組的特征系統(tǒng).雖然改變了方程組形式,但對(duì)于定常問題而言,經(jīng)過足夠長(zhǎng)的時(shí)間后取得的收斂解不會(huì)受此影響.為了更好地計(jì)算黏性問題,通常在帶預(yù)處理矩陣的控制方程中使用原始變量QT(p,u,v,T)T,下面以二維帶預(yù)處理矩陣的有黏性控制方程為例:

其中各變量與通量意義同前文所述.Γ是預(yù)處理矩陣,其形式有多種,常見的形式有Weiss-Smith 矩陣[36]:

其中

k0通常取1,如果是有黏方程,還要限制:

Weiss 與Smith[36]構(gòu)造參考速度Ur,是為了既能在Ma較大時(shí),讓預(yù)處理矩陣還原成非預(yù)處理狀態(tài),又能Ma較小時(shí),減小方程組剛性.此外,對(duì)于黏性流動(dòng)而言,還需要限制Ur不小于局部擴(kuò)散速度μ/(ρΔL),這種限制在擴(kuò)散作用為主導(dǎo)且網(wǎng)格間距較小的情況下是有必要的.對(duì)預(yù)處理后的特征根的討論詳見3.2 節(jié).

Pletcher-Chen 矩陣[37]:

其中

Choi-Merkle 矩陣[33]:

其中

Ma0是參考Mach 數(shù),k2是常數(shù),通過大量實(shí)驗(yàn)發(fā)現(xiàn),Ma0取0.5,k2在Euler 方程時(shí)取0.5,在NS 方程時(shí)取1,這樣有利于長(zhǎng)時(shí)間運(yùn)算的穩(wěn)定性.該方案是Choi 與Turkle[31]提出的,其作用是為避免駐點(diǎn)附近預(yù)處理矩陣的奇異性,對(duì)βMar2實(shí)施一定的人為控制.

3.2 特征系統(tǒng)分析

如果不預(yù)處理,則方程(3)應(yīng)該為

將(11)式等號(hào)左端全部替換成以QT為變量的形式:

即

其中守恒變量Q到原始變量QT的Jacobi 變換矩陣為

下面以直角坐標(biāo)系的ξ方向?yàn)槔治鲱A(yù)處理前后方程組的特征系統(tǒng):

其中

特別地,當(dāng)ξx,ηy時(shí)有

條件數(shù)通常用來衡量方程組的剛性大小,當(dāng)u→ 0 時(shí):

假如快波傳輸了1 個(gè)網(wǎng)格的長(zhǎng)度,則有

可見,當(dāng)條件數(shù)非常大時(shí),相同時(shí)間內(nèi)快慢波傳播的距離會(huì)相差很多數(shù)量級(jí),這給方程的求解造成很大困難.若使用預(yù)處理矩陣替換P0,則由Choi[31],Pletcher[37]和Weiss[36]等提出的預(yù)處理法得到的特征根都有相同的形式:

其中f是關(guān)于Ma的函數(shù),雖然每種預(yù)處理法都不一樣,但該函數(shù)都是改變特征根的關(guān)鍵所在.以Weiss-Smith 形式為例:

其中Ur同(5)式:

當(dāng)Mach 數(shù)很小時(shí),特征根依然能保持在相近的量級(jí)范圍內(nèi),剛性得到降低,條件數(shù)的取值范圍被限制在1

對(duì)此,本文設(shè)計(jì)了一種新的預(yù)處理矩陣,使得方程組在預(yù)處理與非預(yù)處理之間的過度更加光滑,并進(jìn)一步降低了特征系統(tǒng)的條件數(shù),從而增加收斂效率,并且擁有良好的精度和穩(wěn)定性.

3.3 新預(yù)處理矩陣的提出

Turkle[30]曾在Wiess-Smith 矩陣[36]基礎(chǔ)上加以改進(jìn),增加一個(gè)參數(shù)從而提高預(yù)處理矩陣的靈活性,其矩陣具體如下:

式中

其中δ為根據(jù)求解需要選定的0—1 之間的常數(shù),βMar2同(10)式.本文受Turkel 引入調(diào)節(jié)參數(shù)δ的啟發(fā),引入調(diào)節(jié)函數(shù)b1,b2和ΩL,這些函數(shù)會(huì)隨著流場(chǎng)中某些物理量的變化而變化,從而可實(shí)現(xiàn)對(duì)預(yù)處理矩陣動(dòng)態(tài)調(diào)控,具體形式如下:

其中

εm是防止ML=0 的一個(gè)小量,本文取εm=10–6.若Reynolds 數(shù)較小,則需限制其不小于局部擴(kuò)散速度,同時(shí)既要避免駐點(diǎn)附近的奇異性,又要保證ML函數(shù)的光滑性,還要保證后文中偽聲速的光滑性,轉(zhuǎn)化為數(shù)學(xué)分析語言:

1)ML函數(shù)在(0,1)上為單調(diào)遞增的至少C1連續(xù)函數(shù),且(1)=0.

2) 當(dāng)x→ 0 時(shí),ML(x) → 0;ML(1)=1.

3) 若ML≤μ/(ρc△L),則ML=μ/(ρc△L);若ML≥ 1 時(shí),則ML=1.

根據(jù)以上原則構(gòu)造出如(24c)式的ML函數(shù),其中εm同(24b)式:

b1,b2是與Ma相關(guān)的函數(shù),可以相等也可以不相等,在本文中取b1=b2=b,具體的有

為了在低速時(shí)進(jìn)一步降低方程組剛性,在高速時(shí)還原成原方程,需使得

1) 當(dāng)Ma=0 時(shí),B=0;當(dāng)Ma≥ma0時(shí),B=1.

2) 當(dāng)0

3)B'(0)=B''(0)=B'(ma0)=B''(ma0)=0.

4) 在全Mach 數(shù)范圍內(nèi),需要滿足λ3>λ1,2>|λ4|.

5) 在Ma

6) 在Ma >ma0時(shí),要求λ3/λ1,2,|λ4|/λ1,2分別逼近未經(jīng)預(yù)處理時(shí)的λ3/λ1,2,|λ4|/λ1,2.

7) 使B函數(shù)在左右端點(diǎn)Ma=0 和Ma=ma0附近應(yīng)保持變化緩慢,從而抑制速度u在低Mach數(shù)和臨界點(diǎn)處由于預(yù)處理作用而出現(xiàn)的劇變.

8) 與(24c)式的ML匹配,使得后文中偽聲速具備良好的光滑性.

按照以上8 條要求,本文構(gòu)造了如(25b)式的B函數(shù):

其中ma0是事先給定的常數(shù);n1,n2是正整數(shù);CL是常數(shù).為了使B滿足上述5)—7)性質(zhì),且盡可能使特征根差異最小,并使經(jīng)與處理后的特征根滿足λ3>λ1,2>|λ4|,計(jì)算無黏和黏性很小的流動(dòng)時(shí),選擇:

如果黏性不可忽略,可以令

函數(shù)B的圖像可見圖1(d).通過簡(jiǎn)單運(yùn)算可知該矩陣的行列式為

圖1 經(jīng)不同預(yù)處理方法后的特征系統(tǒng)(對(duì)于每種預(yù)處理后的速度有,U'=(λ3+λ4)/2,C'=|λ3–λ4|/2,均為無量綱化速) (a) 條件數(shù)CN;(b) 偽速度U';(c) 偽聲速C';(d) B 函數(shù);(e) 本文預(yù)處理法得到的特征根Fig.1.Characteristic system of the governing equations obtained after different preconditioning:(a) Condition number;(b) pseudospeed;(c) pseudo-sound speed;(d) function B;(e) the eigenvalues obtained after preconditioning in this paper.

可見,該矩陣恒為可逆矩陣.此時(shí)的預(yù)處理矩陣在計(jì)算低速問題時(shí)效果較好,如果用于跨聲速和全速域問題的計(jì)算,可以選取ma0=Mamax(此最大Mach 數(shù)是實(shí)驗(yàn)前的估計(jì)值).

3.4 經(jīng)新預(yù)處理法后的特征系統(tǒng)

3.4.1 特征值和方程組剛性

下面來定量分析該處理法得到的特征系統(tǒng).在二維曲線坐標(biāo)系 (ξ,η)中,以ξ方向?yàn)槔?求出經(jīng)ΓLiu預(yù)處理后Euler 方程的特征根:

其中偽速度和偽聲速分別為

U,Cξ意義同(16) 式,預(yù)處理參數(shù)選用(26a) 式,且U在(0,ma0U)上時(shí),λ3為關(guān)于Ma的單調(diào)遞增函數(shù),與未經(jīng)預(yù)處理時(shí)特征根隨速度U增加而變化的趨勢(shì)一致.經(jīng)過ΓLiu預(yù)處理后的特征根受兩個(gè)函數(shù)B和ML調(diào)控,形式上比Weiss,Pletcher和Choi 等得到的特征根(19a)式擁有更多自由度.當(dāng)Ma→ 0 時(shí)有

此時(shí)所有特征根在不考慮正負(fù)號(hào)的情況下為等價(jià)無窮小,即在Ma很小的時(shí)候,條件數(shù)趨于1.先證明對(duì)于任意Ma >εm,有λ3>λ1,2>|λ4|.

當(dāng)Ma≤1 時(shí):

當(dāng)1

當(dāng)Ma>ma0時(shí):

由此可知對(duì)于任意Ma >εm有λ3>λ1,2,同理可證λ1,2>|λ4|,這也與未經(jīng)預(yù)處理的特征根的大小關(guān)系一致,圖1(e)更直觀地展示了預(yù)處理后的特征根大小關(guān)系.如果ma0>則會(huì)出現(xiàn)Ma0>εm,使得λ3<λ1,2,因此需要限制ma0≤,故上文選取ma0=1.73.λ3,λ1,2恒為正值.當(dāng)Ma<1時(shí),λ4<0;當(dāng)Ma>1 時(shí),λ4>0,與未經(jīng)預(yù)處理時(shí)特征根的符號(hào)相同.特殊地,選取直角坐標(biāo)系,可以求出其特征系統(tǒng)的條件數(shù)滿足:

當(dāng)Ma >ma0時(shí),預(yù)處理矩陣退化成P0,方程恢復(fù)成原控制方程.圖1 為控制方程經(jīng)ΓWS,ΓPC,ΓCM,ΓD以及ΓLiu預(yù)處理后系統(tǒng)的條件數(shù)隨Ma的發(fā)展趨勢(shì).從圖1 可見,ΓD在Ma∞附近效果良好,在駐點(diǎn)處剛性很大,經(jīng)ΓLiu處理后的系統(tǒng)剛性小于其他方法,且特征根隨著Mach 數(shù)從附近增大到Ma>ma0時(shí)的變化過程是光滑的,當(dāng)Ma超過ma0后,偽速度和偽聲速都還原成真實(shí)速度和聲速.由于

其中CFL 數(shù)是差分法中判斷收斂的條件數(shù),當(dāng)CFL 數(shù)固定時(shí),|λmax|減小可以增大時(shí)間步長(zhǎng),從而達(dá)到加速效果.

3.4.2 特征根的連續(xù)性

下面證明特征根λi關(guān)于 (ξ,η) 的連續(xù)性,以直角坐標(biāo)系的x方向?yàn)槔?注意到(28)式中U′,是關(guān)于Ma以1 和ma0為分界點(diǎn)的分段函數(shù),只需要證明分界點(diǎn)處連續(xù)性.經(jīng)無量綱化后,不考慮y方向的速度影響,則有Ma=U=u,從而U′,均可視為Ma(或u)的函數(shù),那么有

當(dāng)u存在關(guān)于x的二階連續(xù)偏導(dǎo)數(shù)時(shí),只需驗(yàn)證U′關(guān)于u的可導(dǎo)性.對(duì)(29a)式和(29b)式分別在(0,1),(1,ma0),(ma0,+∞)上關(guān)于u求導(dǎo),可得

而B在(0,+∞)上存在關(guān)于u的二階連續(xù)導(dǎo)函數(shù),那么易證:

即U′存在關(guān)于x的二階連續(xù)偏導(dǎo)數(shù),同理可證也存在關(guān)于x的二階連續(xù)偏導(dǎo)數(shù).即經(jīng)ΓLiu預(yù)處理后得到的,只要u存在關(guān)于x的二階連續(xù)導(dǎo)數(shù),λi亦存在關(guān)于x的二階連續(xù)導(dǎo)數(shù).從而克服了ΓWS,ΓPC和ΓCM預(yù)處理后的特征根,在預(yù)處理和非預(yù)處理間過渡不夠光滑的缺陷,圖1(a)—(c)也可反映出這點(diǎn).

4 數(shù)值求解

4.1 空間離散

采用Roe 的通量差分格式作為Riemann 求解器[42],下面以ξ方向介紹該方法,其他方向與之同理.帶預(yù)處理的Roe 格式可以表示為

其中

ΛΓ是ΓLiuA的特征根,即(28)式組成的對(duì)角矩陣,A的意義同(14)式,MΓ是將其化為對(duì)角矩陣的變換矩陣,即

此處|ΛΓ|表示將該矩陣所有元素取絕對(duì)值后組成的矩陣.高精度的通量E(Q)L和E(Q)R采用5 階WENO-JS 格式來求得,此過程可詳見文獻(xiàn)[43].對(duì)于黏性項(xiàng)的離散,本文采取一種4 階緊致的中心差分格式,具體可見文獻(xiàn)[44].

4.2 時(shí)間推進(jìn)

為了與空間離散的精度相匹配,采取三階顯式Runge-Kutta 法做時(shí)間推進(jìn)[45]:

4.3 非定常問題的時(shí)間推進(jìn)

本文僅是對(duì)新預(yù)處理法的合理性和可靠性進(jìn)行檢驗(yàn),故選取較為簡(jiǎn)單的定常流動(dòng)來驗(yàn)證.如用于非定常流動(dòng)計(jì)算,預(yù)處理會(huì)破壞時(shí)間精度,因此需要做特殊處理.這里介紹一種應(yīng)用較廣泛的雙時(shí)間步法,其本質(zhì)是在控制方程中填入一項(xiàng)關(guān)于偽時(shí)間τ的導(dǎo)數(shù)項(xiàng),并對(duì)其做預(yù)處理,即

其中τ和t分別是偽時(shí)間和真實(shí)物理時(shí)間,其余物理量同(3)式所述.當(dāng)τ→ ∞時(shí),上式因左邊的第一項(xiàng)消失而復(fù)原為方程.分別對(duì)τ和t離散,該方法包含兩個(gè)循環(huán),即偽時(shí)間τ的內(nèi)循環(huán)和物理時(shí)間t的外循環(huán),在每個(gè)物理時(shí)間步當(dāng)作定常問題用偽時(shí)間推進(jìn)求解,無需時(shí)間精確的偽時(shí)間內(nèi)迭代可以采用預(yù)處理加速收斂.當(dāng)解的殘差變化穩(wěn)定后,即達(dá)到當(dāng)前物理時(shí)間步的真實(shí)解,隨即物理時(shí)間t外循環(huán)推進(jìn)至下一步,然后進(jìn)行偽時(shí)間的內(nèi)迭代并重復(fù)上述步驟,最終得到欲求時(shí)刻的非定常解.

5 數(shù)值模擬驗(yàn)證

下面選取一維有精確解的方程組,平板層流邊界層,方腔頂蓋驅(qū)動(dòng),凸鼓包管道流動(dòng)等定常問題來驗(yàn)證本文方法的可靠性、收斂性和穩(wěn)定性.由于ΓD在遠(yuǎn)離Ma∞的效果不如其他方法,文獻(xiàn)[46]的實(shí)驗(yàn)說明了使用Weiss &Smith 和Choi &Merkle矩陣的效果接近,故本文選取Weiss &Smith,Pletcher &Chen 和本文提出的預(yù)處理矩陣進(jìn)行實(shí)驗(yàn),并對(duì)比得到的結(jié)果.算例5.1 采用matlab2019 編程,算例5.2—5.4 采用fortran90 編程.

5.1 精度測(cè)試

選取有定常精確解的一維Euler 方程組,以QT(p,u,T)T為變量:

其計(jì)算域?yàn)?

由表1 可見,當(dāng)k=0.5 時(shí),速度u會(huì)經(jīng)歷從小于1 (無量綱化)到大于ma0的跨越,使用本文的預(yù)處理法可以在保證精度的同時(shí)實(shí)現(xiàn)從預(yù)處理還原到不處理的連續(xù)過渡;在相同的計(jì)算時(shí)間t=100 下,使用預(yù)處理方法比未處理得到的數(shù)值解精度高出約1—2 個(gè)量級(jí).增加一倍的計(jì)算時(shí)間,到t=200 時(shí)刻未用預(yù)處理的解的精度才勉強(qiáng)達(dá)到t=100 時(shí)刻使用預(yù)處理的精度.當(dāng)k進(jìn)一步減小時(shí),預(yù)處理的效果更加明顯,從圖2(c)—(f)可以清晰地看出,預(yù)處理比不處理收斂速度更快.這表明本文的預(yù)處理方法能夠在保證精度的同時(shí),提高計(jì)算的效率.

圖2 數(shù)值方法得到的解曲線與精確解曲線 (a) k=0.5;(b) k=0.04;(c) k=0.01;(d) k=5×10–3;(e) k=1×10–3;(f) k=2×10–4Fig.2.Solution curves obtained by numerical solution and exact solution curves:(a) k=0.5;(b) k=0.04;(c) k=0.01;(d) k=5×10–3;(e) k=1×10–3;(f) k=2×10–4.

表1 數(shù)值解與精確解之間的誤差Table 1. Errors between numerical solutions and exact solutions.

5.2 平板層流邊界層

一長(zhǎng)為L(zhǎng)的平板靜置在水平面上,從左側(cè)的與平板平行的均勻來流通過其上表面,以板長(zhǎng)為參考長(zhǎng)度,以來流處的聲速為參考速度進(jìn)行無量綱化,給定來流速度u∞=0.1,基于板長(zhǎng)的Reynolds數(shù)Re∞=106,設(shè)定來流處壓強(qiáng)和溫度以及出口的壓強(qiáng),計(jì)算網(wǎng)格為180×50,其中x方向?yàn)榫鶆蚓W(wǎng)格,y方向隨著y減小而加密,貼近平板上表面的第1 層網(wǎng)格厚度為0.002,CFL=5.計(jì)算到收斂后,選取x=0.8 處的速度u隨y的分布,并與Bulsius 準(zhǔn)解析解對(duì)比,對(duì)比結(jié)果如圖3 所示,其中

圖3 Ma=0.1 時(shí)邊界層流動(dòng)的速度剖面Fig.3.Velocity profile of the boundary layer flow when Ma=0.1.

本文預(yù)處理法的收斂速度同Weiss &Smith和Pletcher &Chen 方法的對(duì)比如圖4 所示.可見,本文的預(yù)處理方法得到的數(shù)值解與Blasius 解基本一致.由圖4 可見,本文方法的收斂效率略好于Weiss &Smith 和Pletcher &Chen 方法.

圖4 邊界層流動(dòng)的收斂歷程Fig.4.Convergence histories of the boundary layer flow when Ma=0.1.

5.3 空腔頂蓋驅(qū)動(dòng)

此算例描述的是一邊長(zhǎng)為1 (此算例中變量均無量綱化)的正方形空腔,初始Mach 數(shù)Ma0=1,頂蓋以u(píng)0=1 速度向右運(yùn)動(dòng),Reynolds 數(shù)分別為100,1000,10000.計(jì)算網(wǎng)格為60×60 成中心對(duì)稱的非均勻網(wǎng)格,并在壁面附近加密,CFL=100,分別使用本文預(yù)處理法和Weiss &Smith 和Pletcher &Chen 方法計(jì)算此問題.圖5—圖7 是計(jì)算到收斂后得到的空腔內(nèi)部流場(chǎng)線圖,可以觀察到,在空腔的幾何中心附近產(chǎn)生一個(gè)大渦,在底部?jī)山歉浇a(chǎn)生兩個(gè)較小的渦結(jié)構(gòu).隨著Reynolds 數(shù)的增大,大渦中心逐漸向空腔幾何中心移動(dòng),底部小渦逐漸增大,當(dāng)Reynolds 數(shù)為10000 時(shí),左上角附近也產(chǎn)生一個(gè)渦結(jié)構(gòu).圖8 是本文水平與垂直中心線上速度與使用多重網(wǎng)格法的計(jì)算結(jié)果[47]的對(duì)比,圖9 是三種預(yù)處理法的收斂速度對(duì)比.

圖5 Re=100 時(shí)的空腔流線圖 (a) 本文方法;(b) Weiss &Smith 方法;(c) Pletcher &Chen 方法Fig.5.Velocity streamline diagram of cavity,Re=100:(a) The new method;(b) method of Weiss &Smith;(c) method of Pletcher &Chen.

圖6 Re=1000 時(shí)的空腔流線圖 (a) 本文方法;(b) Weiss &Smith 方法;(c) Pletcher &Chen 方法Fig.6.Velocity streamline diagram of cavity,Re=1000:(a) The new method;(b) method of Weiss &Smith;(c) method of Pletcher &Chen.

圖7 Re=10000 時(shí)的空腔流線圖 (a) 本文方法;(b) Weiss &Smith 方法;(c) Pletcher &Chen 方法Fig.7.Velocity streamline diagram of cavity,Re=10000:(a) The new method;(b) method of Weiss &Smith;(c) method of Pletcher &Chen.

通過圖8 和圖9 可知,本文提出的方法精確度略高于Weiss &Smith 和Pletcher &Chen 方法(這兩種方法的精度對(duì)比可參考文獻(xiàn)[46]).本文方法的收斂速度也更快一些,且隨著黏性效應(yīng)增加(Reynolds 數(shù)減小)效果越明顯,這表明本文格式可以實(shí)現(xiàn)對(duì)低速有黏問題的高精度高效率的模擬.

圖8 本文預(yù)處理后的中心線速度分布 (a)水平中心線上的速度v;(b)垂直中心線上的速度uFig.8.Velocity along lines through geometric center caculated by the new method:(a) v-velocity along horizontal center line;(b) vvelocity along horizontal center line.

圖9 方腔流動(dòng)的收斂歷程對(duì)比 (a) Re=100;(b) Re=1000;(c) Re=10000Fig.9.Convergence histories of the the laminar flow in a square cavity:(a) Re=100;(b) Re=1000;(c) Re=10000.

5.4 帶凸鼓包的二維管道流動(dòng)

此算例[48]目的在于驗(yàn)證本文提出的預(yù)處理方法在曲線坐標(biāo)系下的計(jì)算性能,故選擇下表面帶鼓包的二維管道流動(dòng)問題,該算例模型和流場(chǎng)相對(duì)簡(jiǎn)單,便于對(duì)結(jié)果進(jìn)行分析.具體模型如下:管道的長(zhǎng)度為3 (單位均無量綱化),高度為1,凸鼓包位于下表面的中心位置,長(zhǎng)度為1,其尺寸為y=0.1sin2πx(1

圖11 為三種方法的殘差收斂曲線對(duì)比.Ma0=0.05 時(shí)三種方法得到的流場(chǎng)Mach 數(shù)等值線分布見圖12(a)—(c).圖12(d)為在加密網(wǎng)格(600×200)下使用Weiss &Smith 法計(jì)算得到的流場(chǎng),并將其作為準(zhǔn)精確解.在低速流動(dòng)下,該流場(chǎng)的Mach 等值線接近軸對(duì)稱.可見本文方法的計(jì)算結(jié)果與準(zhǔn)精確解基本一致,相比另外兩種方法精度更高,尤其第3—5 條等值線的位置.另外本文方法的收斂效率也略高于其他方法,雖然Ma0=0.2 時(shí)由于剛性并不是很大,加速收斂的效果比Ma0更低時(shí)明顯,但最終的殘差均小于其他預(yù)處理法,這說明本文格式在保證精度的同時(shí),擁有更良好的穩(wěn)定性.

圖12 Ma0=0.05 時(shí)流場(chǎng)的等Mach 線 (a) 精確解;(b) Weiss &Smith;(c) Pletcher &Chen;(d) 本文方法Fig.12.Mach contour of the two-dimensional bump flow problem,Ma0=0.05:(a) Exact solution;(b) Weiss &Smith;(c) Pletcher &Chen;(d) the new method in this paper.

6 結(jié)論

針對(duì)低速流動(dòng)問題,本文基于虛擬壓縮法的思想,提出一種以當(dāng)?shù)豈ach 數(shù)、來流Mach 數(shù)、密度、溫度、每個(gè)維度的分速度和Reynolds 數(shù)為變量的預(yù)處理矩陣,相比Weiss,Choi,Dailey,Pletcher和Turkle 等提出的矩陣,把控制方程的剛性在Ma<0.5 時(shí)進(jìn)一步降低到1 左右,使得特征根之間的量級(jí)差異進(jìn)一步減小,從而提升了求解效率,同時(shí)改進(jìn)了從預(yù)處理到關(guān)閉預(yù)處理的過渡不夠光滑的不足.

展開數(shù)值實(shí)驗(yàn)驗(yàn)證了新方法的可靠性.通過一維算例驗(yàn)證了新方法的準(zhǔn)確性,將預(yù)處理和不處理的結(jié)果加以對(duì)比,證明了預(yù)處理方法可加速計(jì)算低速流動(dòng)問題的收斂.通過平板層流邊界層和黏性方腔流動(dòng)問題進(jìn)一步驗(yàn)證了本文方法在二維低速流動(dòng)問題中可以獲得高精度數(shù)值解,并且擁有更高效的收斂性,隨著黏性作用越強(qiáng),這種高效性越明顯.通過二維鼓包管道流動(dòng)問題,驗(yàn)證了本文方法在曲線坐標(biāo)系下仍能保持良好的精度和收斂性,而且具備更好的穩(wěn)定性.目前本文重點(diǎn)在于改進(jìn)算法,故選取算例的結(jié)構(gòu)比較簡(jiǎn)單,未來將會(huì)嘗試模擬復(fù)雜網(wǎng)格,非定常流動(dòng)或引入湍流模型.

感謝中國(guó)科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室的李詩(shī)堯博士,天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室的陳嘉禹碩士為本文提供的計(jì)算資源;感謝天津大學(xué)數(shù)學(xué)學(xué)院2018 級(jí)程啟豪同學(xué),南開大學(xué)數(shù)學(xué)學(xué)院學(xué)習(xí)部2018 級(jí)左晨琛、2019 級(jí)的戚紅利、2020 級(jí)趙振岐等同學(xué)為本文的矩陣特征系統(tǒng)在曲線坐標(biāo)系下數(shù)學(xué)性質(zhì)以及特定幾何結(jié)構(gòu)從笛卡爾坐標(biāo)系下到曲線坐標(biāo)下變換涉及的復(fù)分析和拓?fù)鋯栴}的討論;感謝國(guó)防科技大學(xué)空天科學(xué)學(xué)院田正雨教授,南京航空航天大學(xué)航空學(xué)院陳紅全教授,南京航空航天大學(xué)數(shù)學(xué)系張燕博士為本文提出建議.