劉君，魏雁昕，韓芳

大連理工大學(xué) 航空航天學(xué)院，大連 116024

有限差分法應(yīng)用于具有復(fù)雜外形的網(wǎng)格時(shí)需要進(jìn)行坐標(biāo)變換，即采用數(shù)學(xué)方法將流體力學(xué)控制方程由笛卡爾坐標(biāo)系變換到與物體形狀相吻合的貼體坐標(biāo)系(Body-Fitted Coordinates，BFC)以后再進(jìn)行格式離散。在實(shí)際工程問題中，很少能直接給出坐標(biāo)變換的解析表達(dá)式，因此，網(wǎng)格生成技術(shù)成為計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)的一個(gè)重要分支。離散是對(duì)連續(xù)的近似表達(dá)，在生成的網(wǎng)格點(diǎn)基礎(chǔ)上計(jì)算坐標(biāo)變換系數(shù)存在誤差。例如，一維函數(shù)u=u(x)從物理空間變換到計(jì)算空間ξ=ξ(x)，導(dǎo)數(shù)ux變成ux=ξxuξ，可以看出，如果坐標(biāo)變換系數(shù)ξx采用數(shù)值算法得到，則必然也會(huì)引入誤差。算例表明，即使坐標(biāo)變換系數(shù)ξx采用解析值，誤差依然存在，本文把這種由坐標(biāo)變換引起的誤差稱為坐標(biāo)變換誘導(dǎo)誤差。

早在1961年Trulio和Trigger[1]就注意到了這個(gè)問題，提出了一維問題中的“體積守恒(conservation of volume)”概念。1978年Steger等[2-4]研究了從守恒型控制方程出發(fā)計(jì)算具有復(fù)雜幾何外形的流場(chǎng)時(shí)坐標(biāo)變換系數(shù)的影響，即使是均勻流場(chǎng)，采用傳統(tǒng)坐標(biāo)變換系數(shù)求解方法也會(huì)產(chǎn)生非零源項(xiàng)，該源項(xiàng)由坐標(biāo)變換引起，在流場(chǎng)中造成非物理解，嚴(yán)重時(shí)會(huì)引起數(shù)值振蕩，影響計(jì)算穩(wěn)定性[5-7]。1979年Thomas和Lombard[8-9]在前人研究基礎(chǔ)上提出了守恒型的坐標(biāo)變換系數(shù)計(jì)算方法，并正式提出了被后來廣泛引用的幾何守恒律(Geometric Conservation Law，GCL)概念。根據(jù)計(jì)算過程中網(wǎng)格是否隨時(shí)間變形，GCL被進(jìn)一步區(qū)分為“體積守恒律”(Volume Conservation Law，VCL)和“面積守恒律”(Surface Conservation Law，SCL)，前者解決“任意拉格朗日-歐拉”(Arbitrary Lagrange-Euler，ALE)方程動(dòng)網(wǎng)格問題，應(yīng)用范圍遠(yuǎn)沒有后者普遍，目前在高精度格式領(lǐng)域研究GCL常指后者。本文后續(xù)提到的GCL也指“面積守恒律”。

按照文獻(xiàn)[9]中網(wǎng)格導(dǎo)數(shù)變換理論，如果使用二階精度格式離散uξ，那么采用中心差分格式離散ξx不影響導(dǎo)數(shù)ux的精度?；谏鲜稣J(rèn)識(shí)，對(duì)GCL研究工作基本停頓下來，隨著高階精度格式的發(fā)展與應(yīng)用這個(gè)問題又引起重視。從1999年Gaitonde和Visbal[10]提出了緊致格式的GCL以來，近20年來對(duì)基于高精度有限差分的GCL的研究得到了充分發(fā)展，其中2011年鄧小剛等證明了滿足GCL的充分條件[11]并提出了坐標(biāo)變換系數(shù)的守恒計(jì)算方法——C MM(Conservative Metric Method)算法；隨后，在2013年進(jìn)一步發(fā)展出了坐標(biāo)變換系數(shù)計(jì)算結(jié)果唯一的對(duì)稱守恒計(jì)算方法——SC MM算法(Symmetrical Conservative Metric Method)[12]，該研究成果對(duì)于緊致類格式具有普適性指導(dǎo)意義。近年來科研工作者們對(duì)構(gòu)建WENO格式的GCL時(shí)常借鑒其思想，如將WENO格式分解成中心差分部分及數(shù)值耗散部分[13-15]，中心部分直接采用SC MM計(jì)算。

文獻(xiàn)[15]的觀點(diǎn)“面積守恒律的簡(jiǎn)化過程中，用到兩個(gè)數(shù)學(xué)前提，即網(wǎng)格物理坐標(biāo)對(duì)計(jì)算坐標(biāo)連續(xù)可導(dǎo)和混合導(dǎo)數(shù)都連續(xù)。但是在實(shí)際計(jì)算過程中，不能保證計(jì)算網(wǎng)格坐標(biāo)點(diǎn)等間距或二階可導(dǎo)，簡(jiǎn)化過程中所做假設(shè)不成立，并且網(wǎng)格導(dǎo)數(shù)由數(shù)值差分格式計(jì)算得到，因此引入數(shù)值誤差”在GCL研究領(lǐng)域具有普遍性。為了驗(yàn)證上述結(jié)論，針對(duì)坐標(biāo)變換函數(shù)的各階導(dǎo)數(shù)和混合導(dǎo)數(shù)均存在且連續(xù)的若干算例進(jìn)行考查。

1 數(shù)值現(xiàn)象

目前常見的高精度格式大多包含隨著流場(chǎng)變化的模板或限制器，對(duì)于參數(shù)變化的流場(chǎng)難以區(qū)分誤差來源是差分格式還是幾何誘導(dǎo)。因此，選擇超聲速均勻流動(dòng)，來流馬赫數(shù)M∞=2，無量綱流動(dòng)參數(shù)(ρ,u,v,p)=(1.4,2,0,1)，其中ρ、u、v、p分別為無量綱密度、X方向速度、Y方向速度和壓力。從二維Euler方程出發(fā)進(jìn)行數(shù)值模擬，首先給出柱坐標(biāo)算例。

1.1 柱坐標(biāo)均勻網(wǎng)格的坐標(biāo)變換誘導(dǎo)誤差

由柱坐標(biāo)和直角坐標(biāo)之間的變換函數(shù)及其逆函數(shù)可以寫出解析表達(dá)式：

(1)

式中：r為圓環(huán)區(qū)域半徑；x為X軸方向坐標(biāo)，x=rcosθ；y為Y軸方向坐標(biāo)，y=rcosθ；θ為徑向網(wǎng)格線與X軸正向夾角。

計(jì)算區(qū)域?yàn)閞∈[1,2]的圓環(huán)，沿徑向r和周向θ網(wǎng)格均勻分布，徑向網(wǎng)格點(diǎn)總數(shù)Nξ=51，周向網(wǎng)格點(diǎn)總數(shù)Nη=(360/7.2)+Nη0，其中Nη0為θ=0°網(wǎng)格線沿圓周方向編號(hào)。為便于處理周向邊界條件，進(jìn)行計(jì)算域的拓展，例如，一階迎風(fēng)格式沿周向拓展2個(gè)點(diǎn)，編碼η=Nη0=1和η=51 重合于x軸(即θ=0°網(wǎng)格線)，編碼η=0和η=50 對(duì)應(yīng)θ=-7.2°的網(wǎng)格線，編碼η=2和η=52均是θ=7.2°的網(wǎng)格線，對(duì)于七點(diǎn)WENO格式需拓展更多的點(diǎn)，此時(shí)Nη0=3，但是網(wǎng)格夾角保持不變。

基于以上設(shè)定條件，計(jì)算空間坐標(biāo)和柱坐標(biāo)之間的線性關(guān)系為

(2)

計(jì)算中用到的5個(gè)坐標(biāo)變換系數(shù)解析求出：

(3)

對(duì)流項(xiàng)計(jì)算采用一階迎風(fēng)離散格式及Van Leer通量分裂格式，計(jì)算收斂以后數(shù)值解和理論值的誤差就是坐標(biāo)變換誘導(dǎo)誤差，其流動(dòng)參數(shù)誤差的分布規(guī)律類似，此處以密度誤差為例進(jìn)行分析。當(dāng)流場(chǎng)徑向內(nèi)外圓邊界根據(jù)Riemann不變量處理時(shí)，其密度誤差δρ分布如圖1(a)所示，由于超聲速流動(dòng)擾動(dòng)不會(huì)向上游傳播，流場(chǎng)中x<0區(qū)域(左側(cè)半圓)的計(jì)算結(jié)果定性合理，但是在x>0 區(qū)域存在明顯的邊界影響。采用固定邊界值的方法處理徑向內(nèi)外圓邊界也可以得到收斂解，其密度誤差分布云圖如圖1(b)所示，可以看出，流場(chǎng)參數(shù)誤差分布左右兩側(cè)不對(duì)稱特性依然存在。當(dāng)徑向外圓邊界采用Riemann不變量處理時(shí)，內(nèi)圓邊界把x<0出口值作為在y相等高度上x>0對(duì)應(yīng)點(diǎn)的入口值，其計(jì)算結(jié)果如圖1(c)所示，可以看出流場(chǎng)參數(shù)誤差在x=0左右兩側(cè)呈現(xiàn)對(duì)稱分布。由此可以推斷：圖1(a)和圖1(b)中流場(chǎng)密度誤差分布出現(xiàn)左右兩側(cè)不對(duì)稱的原因是在x>0內(nèi)環(huán)進(jìn)口邊界處理中沒有考慮坐標(biāo)變換引起的誤差。

圖1 3種邊界條件下密度誤差云圖

定量比較圖1密度誤差特性，可以看出圓柱上游誤差分布相等且不受內(nèi)圓邊界處理影響，因此，僅對(duì)x<0區(qū)域的密度誤差分布進(jìn)行分析。為了后續(xù)研究，坐標(biāo)變換系數(shù)除了采用式(3)求導(dǎo)的理論值，還采用二階中心格式的數(shù)值解進(jìn)行計(jì)算，圖2為不同的系數(shù)計(jì)算方法和不同的差分格式下密度誤差δρ比較，為便于對(duì)比兩種格式，誤差云圖選取相同幅值。從結(jié)果看，WENO格式的誤差比一階迎風(fēng)格式更大。

圖2 柱坐標(biāo)系密度誤差云圖

同心圓網(wǎng)格算例表明進(jìn)出口邊界變換成曲線以后，坐標(biāo)變換也有影響，在超聲速流動(dòng)穩(wěn)定性、氣動(dòng)聲學(xué)、湍流轉(zhuǎn)捩等研究中，上游微小的擾動(dòng)會(huì)導(dǎo)致結(jié)果出現(xiàn)差異，數(shù)值模擬過程中應(yīng)該分辨和評(píng)估這種進(jìn)口邊界誤差的不均勻特性對(duì)結(jié)果造成的干擾。為了排除進(jìn)口邊界的影響，進(jìn)行1.2節(jié)的數(shù)值實(shí)驗(yàn)。

1.2 正弦函數(shù)網(wǎng)格的坐標(biāo)變換誘導(dǎo)誤差

設(shè)計(jì)了如圖3(a)所示的拼接網(wǎng)格，左側(cè)第1塊網(wǎng)格(x≤0.6)完全是直角坐標(biāo)系下的均勻網(wǎng)格：

(4)

式中：0≤ξ≤30, 0≤η≤100。

第2塊網(wǎng)格(x>0.6)按照如式(5)所示的變換函數(shù)生成：

(5)

式中：31≤ξ≤130, 0≤η≤100。

計(jì)算中用到的坐標(biāo)變換系數(shù)可以解析求出，5個(gè)系數(shù)中有2個(gè)為常數(shù)：

(6)

另外3個(gè)系數(shù)J、ξx、ξy的分布如圖3所示，JACB代表雅可比行列式，可以看出，在拼接線x=0.6處整體來說存在明顯不連續(xù)，但是在y=-1,0,1坐標(biāo)處系數(shù)基本上光滑過渡。

圖3 正弦網(wǎng)格及其變換系數(shù)

采用Steger通量分裂格式，不同離散格式和坐標(biāo)系數(shù)計(jì)算方法得到的密度誤差特性如圖4所示，可以看出邊界處理對(duì)計(jì)算結(jié)果沒有影響。同樣WENO格式的計(jì)算誤差比一階迎風(fēng)格式更大，其原因?qū)⒃诤笪闹杏懻?，這里主要關(guān)注坐標(biāo)系數(shù)精度的影響。

圖2和圖4的計(jì)算結(jié)果均表明坐標(biāo)系數(shù)計(jì)算方法影響坐標(biāo)變換誘導(dǎo)誤差，采用二階中心格式數(shù)值解的誤差比根據(jù)式(3)解析求得的誤差小兩個(gè)數(shù)量級(jí)。

圖4 不同坐標(biāo)變換系數(shù)下正弦網(wǎng)格密度誤差云圖

GCL研究的理論基礎(chǔ)如前所述，從Steger等研究GCL開始[2-4]，就是針對(duì)差分格式構(gòu)建更高精度的坐標(biāo)系數(shù)計(jì)算方法，從而實(shí)現(xiàn)不影響方程整體離散精度的目標(biāo)，例如，按照滿足SC MM算法計(jì)算ξx，對(duì)于五階精度WCNS-E-5格式，需要六階中心插值[16]，同樣對(duì)于WENO格式ξx也是形式復(fù)雜的高精度插值[13-15]。原則上坐標(biāo)系數(shù)計(jì)算方法的精度極限是準(zhǔn)確的理論值，圖2和圖4 密度誤差云圖中數(shù)值解的誤差比解析解的還要大，這種現(xiàn)象與GCL的研究目標(biāo)似乎相矛盾。

從同心圓算例和正弦網(wǎng)格算例計(jì)算結(jié)果可以看出，在“網(wǎng)格物理坐標(biāo)對(duì)計(jì)算坐標(biāo)連續(xù)可導(dǎo)和混合導(dǎo)數(shù)都連續(xù)”這兩個(gè)數(shù)學(xué)前提成立條件下，即使坐標(biāo)系數(shù)采用理論值誤差依然存在，由此看來，認(rèn)為坐標(biāo)變換誘導(dǎo)誤差來自變換導(dǎo)數(shù)不連續(xù)的觀點(diǎn)不準(zhǔn)確，不是坐標(biāo)系數(shù)計(jì)算方法及其精度可以解決的問題。

目前建立的GCL算法針對(duì)不同的對(duì)流項(xiàng)離散格式，還沒有考慮通量分裂格式這個(gè)因素。圖5 是一階迎風(fēng)格式條件下坐標(biāo)變換系數(shù)解析求解采用AUSM、HLLC、ROE、Van Leer 4種通量分裂格式的壓力誤差δp云圖。采用WENO格式得到的誤差更大，分布類似，這里不再顯示。

從圖5中可以看出，不同通量分裂格式的坐標(biāo)變換誘導(dǎo)誤差存在差異，研究GCL時(shí)應(yīng)該考慮這個(gè)因素。GCL把坐標(biāo)變換誘導(dǎo)誤差歸因于網(wǎng)格，但是調(diào)研未看到網(wǎng)格質(zhì)量因素如何影響的定量研究文獻(xiàn)。網(wǎng)格生成技術(shù)研究領(lǐng)域經(jīng)常采用正交性、光滑性和合理分布性評(píng)價(jià)網(wǎng)格質(zhì)量[17]，前述柱坐標(biāo)算例不但變換函數(shù)的導(dǎo)數(shù)連續(xù)，而且均勻網(wǎng)格完全正交且光滑，依然存在坐標(biāo)變換誘導(dǎo)誤差，說明這些網(wǎng)格質(zhì)量評(píng)價(jià)因素未能反映問題的本質(zhì)，設(shè)計(jì)1.3節(jié)的算例對(duì)這個(gè)問題進(jìn)行數(shù)值實(shí)驗(yàn)。

圖5 不同通量分裂格式的壓力誤差

1.3 網(wǎng)格正交性和光滑性的影響

在柱坐標(biāo)結(jié)果(圖1和圖2)分析中，注意到誤差分布具有軸對(duì)稱特征，考察計(jì)算中用到的5個(gè)坐標(biāo)變換系數(shù)分布，只有雅可比行列式J是軸對(duì)稱的。在正弦函數(shù)網(wǎng)格中沿著y=-1,0,1這3條線點(diǎn)J變化較小，對(duì)應(yīng)的誤差也較小。為考察J和坐標(biāo)變換誘導(dǎo)誤差之間的關(guān)系，設(shè)計(jì)和計(jì)算了正交網(wǎng)格和非正交網(wǎng)格。對(duì)于網(wǎng)格線與坐標(biāo)軸平行的完全正交網(wǎng)格，坐標(biāo)系數(shù)ξy=ηx=0，所有算例中的J、ξx和ηy同時(shí)不為常數(shù)，圖6(a)和圖6(b)為指數(shù)加密網(wǎng)格及J云圖。可以看出存在明顯的不均勻特性。除了常用的指數(shù)加密，還計(jì)算了Δxi+1=αΔxi和Δyi+1=βΔyi等比例加密，其中Δxi和Δyi為i點(diǎn)處網(wǎng)格尺度，α和β為參數(shù)，計(jì)算網(wǎng)格如圖6(c)所示，計(jì)算中坐標(biāo)變換系數(shù)采用二階中心差分的計(jì)算值，該網(wǎng)格下J云圖如圖6(d)所示。圖6(a)和圖6(c)兩種網(wǎng)格下的誤差云圖如圖7所示，說明單純J變化不一定產(chǎn)生坐標(biāo)變換誘導(dǎo)誤差。由于相鄰網(wǎng)格之間的網(wǎng)格尺度Δx和Δy均不相等，從網(wǎng)格質(zhì)量評(píng)價(jià)因素看光滑性不夠好，因此選取α=β=2.00得到的網(wǎng)格光滑性較差，且網(wǎng)格節(jié)點(diǎn)數(shù)減少為圖6(a)和圖6(c)中網(wǎng)格的1/10，該網(wǎng)格實(shí)際應(yīng)用很少采用，計(jì)算結(jié)果也沒有發(fā)現(xiàn)誤差，至少表明對(duì)均勻流動(dòng)來說，光滑性不會(huì)引起坐標(biāo)變換誘導(dǎo)誤差。

將圖6(a)中的直角網(wǎng)格變成如圖8(a)所示的傾斜網(wǎng)格，在不光滑網(wǎng)格基礎(chǔ)上考查正交性的影響。網(wǎng)格線傾斜以后除ηx以外其余4個(gè)坐標(biāo)系數(shù)均不為常數(shù)，定義網(wǎng)格y軸與x軸負(fù)向的夾角為θ′，例如θ′=45°工況下J、ξx和ξy分布如圖8(b)～圖8(d)所示。在0°<θ′<180°范圍內(nèi)，選擇夾角參數(shù)θ′=m×15°(m=1,2,…,11)進(jìn)行數(shù)值實(shí)驗(yàn)，圖9給出了夾角參數(shù)θ′=30°、θ′=45°和θ′=120° 的計(jì)算結(jié)果，所有算例的誤差云圖和圖7相同，均沒有坐標(biāo)變換誘導(dǎo)誤差，說明在光滑性較差的情況下，夾角保持常數(shù)時(shí)正交性較差的網(wǎng)格也不會(huì)產(chǎn)生誤差。

圖6 正交網(wǎng)格參數(shù)

圖7 3種網(wǎng)格誤差云圖

圖8 45°夾角傾斜網(wǎng)格及其坐標(biāo)系數(shù)分布

圖9 不同夾角傾斜網(wǎng)格密度誤差云圖

通過以上數(shù)值現(xiàn)象，可以得出：

1) 坐標(biāo)變換誘導(dǎo)誤差來自“網(wǎng)格物理坐標(biāo)對(duì)計(jì)算坐標(biāo)的變換函數(shù)的導(dǎo)數(shù)不連續(xù)”的觀點(diǎn)不準(zhǔn)確。

2) 不能僅采用網(wǎng)格的正交性和光滑性評(píng)估坐標(biāo)變換誘導(dǎo)誤差的影響特性。

2 坐標(biāo)變換誘導(dǎo)誤差的機(jī)理

將笛卡爾坐標(biāo)系下(t,x,y,z)的三維守恒型Euler方程變換到貼體坐標(biāo)系(τ,ξ,η,ζ)，其中t為時(shí)間，τ為貼體坐標(biāo)系下的時(shí)間，ζ為三維空間貼體坐標(biāo)下與ξ、η軸成右手系的坐標(biāo)。其中得到原始形式：

UIt+FIx+GIy+HIz=S

(7)

如果不考慮“體積守恒律”，那么It=0。“面積守恒律”方程為

(8)

式中：ξz、ηz、ζx、ζy和ζz均為坐標(biāo)變換度量系數(shù)。以Ix=0為例，根據(jù)逆變換關(guān)系式得

(9)

式中：yζ、zξ、zη和zζ均為逆變換度量系數(shù)。繼續(xù)求導(dǎo)，以式(9)中第1項(xiàng)為例，得

y″η ξzζ+yηz″ζ ξ-z″η ξyζ-zηy″ζ ξ

(10)

式中：上標(biāo)″表示求二階偏導(dǎo)數(shù)。

原則上在給定網(wǎng)格點(diǎn)(x,y,z)ijk位置以后，可以按照式(3)計(jì)算ξx，但是由于存在計(jì)算誤差，導(dǎo)致Ix=0不成立，即使對(duì)均勻來流也存在引起質(zhì)量不守恒的現(xiàn)象，這就是產(chǎn)生坐標(biāo)變換誘導(dǎo)誤差的機(jī)理。

根據(jù)式(10)可以看出，即使yη和ξx采用準(zhǔn)確的理論值，由于還有?/?ξ的離散誤差，恒等式也是不成立的，依然存在坐標(biāo)變換誘導(dǎo)誤差，可以合理解釋1.1節(jié)和1.2節(jié)算例采用理論值也出現(xiàn)坐標(biāo)變換誘導(dǎo)誤差的數(shù)值現(xiàn)象。坐標(biāo)變換系數(shù)采用數(shù)值解抹平了yη和ξx的變化斜率，對(duì)應(yīng)的?/?ξ也較小，形成的源項(xiàng)S也較小，由此可以解釋圖2和圖4中出現(xiàn)的二階中心格式數(shù)值解的誘導(dǎo)誤差比準(zhǔn)確的理論值還要小的數(shù)值現(xiàn)象。同樣，考查1.3節(jié)網(wǎng)格正交性和光滑性的算例，如果在離散以后恒等式還成立，那么也沒有誤差。

認(rèn)識(shí)到坐標(biāo)變換誘導(dǎo)誤差的機(jī)理和坐標(biāo)變換系數(shù)的理論值無法消除誤差，那么就通過構(gòu)建yη、ξx和?/?ξ相互匹配的算法使在離散空間上恒等式Ix=Iy=Iz=0重新成立，這就是GCL的出發(fā)點(diǎn)。由于yη和ξx只能計(jì)算得到，在GCL研究領(lǐng)域有個(gè)無法擺脫的邏輯困境：坐標(biāo)變換系數(shù)采用解析值是不滿足GCL算法的。

3 幾何守恒律算法

從量綱角度分析，yη等參數(shù)對(duì)應(yīng)矢量，ξx等參數(shù)是矢量面積，J-1對(duì)應(yīng)體積，如果采用二階中心差商計(jì)算矢量，涉及離散點(diǎn)(x,y,z)ijk及其相鄰網(wǎng)格點(diǎn)正好形成六面體，在網(wǎng)格分布非等距情況下，采用不同算法得到的yη結(jié)果不同，網(wǎng)格點(diǎn)不共面情況下矢量面積ξx也與算法相關(guān)，得到的體積也不唯一，因此，早期CFD應(yīng)用中Steger等解釋Ix≠0的原因就是矢量、面積和體積的算法不相容，采用一種加權(quán)平均算法計(jì)算坐標(biāo)變換系數(shù)，解決了二階中心格式離散流動(dòng)方程時(shí)遇到的問題[2-4]。后來研究者發(fā)現(xiàn)，如果流動(dòng)方程的離散格式非二階中心格式，即使采用Steger等幾何相容滿足Ix=0等面積守恒律，依然存在均勻來流質(zhì)量不守恒的現(xiàn)象，認(rèn)識(shí)到計(jì)算式(10)中y″ηξ等二階混合導(dǎo)數(shù)與流動(dòng)方程的離散格式相關(guān)，不同格式需要相適應(yīng)GCL算法，其中鄧小剛等[12]提出的SC MM實(shí)際上是一個(gè)普適性的準(zhǔn)則，應(yīng)用該準(zhǔn)則能夠滿足面積守恒律要求，即Ix=Iy=Iz=0。

按照SC MM準(zhǔn)則，如果計(jì)算過程中流動(dòng)方程算子發(fā)生變化，坐標(biāo)變換系數(shù)的算子也要做相應(yīng)調(diào)整。例如對(duì)于一階迎風(fēng)格式，在時(shí)間迭代過程中流動(dòng)算子隨著參數(shù)動(dòng)態(tài)變化，滿足SC MM要求的坐標(biāo)變換系數(shù)算子也要具有迎風(fēng)特性，且每步需要重新計(jì)算，這給具體應(yīng)用帶來困難。目前SC MM主要應(yīng)用于格式模板系數(shù)確定的格式，從給出的均勻流動(dòng)計(jì)算看，對(duì)于可以寫成若干中心格式組合的高階精度緊致類格式效果明顯，對(duì)于WENO格式不能完全消除坐標(biāo)變換誘導(dǎo)誤差[14]，推測(cè)與其中的迎風(fēng)型模板有關(guān)。

目前的GCL算法僅與格式相關(guān)，如果再將通量分裂格式、限制器等因素增加進(jìn)來，沿著這條技術(shù)路線消除坐標(biāo)變換誘導(dǎo)誤差構(gòu)建的算法會(huì)非常復(fù)雜，預(yù)測(cè)將給應(yīng)用帶來困難。

4 離散等價(jià)方程算法

除了第3節(jié)中所述GCL外，另一種思路是2018年劉君等提出離散等價(jià)方程算法[18]，放棄對(duì)恒等式Ix=Iy=Iz=0的約束條件，在離散對(duì)流項(xiàng)的同時(shí)也離散S≠0源項(xiàng)，以此消除坐標(biāo)變換誘導(dǎo)誤差。下面對(duì)該算法進(jìn)行簡(jiǎn)單介紹。

坐標(biāo)變換系數(shù)在以上三維守恒型Euler方程通量中是線性的，引入

(11)

對(duì)流項(xiàng)通量和源項(xiàng)表示為

(12)

(13)

(14)

根據(jù)Steger和Van Leer通量分裂格式的通用表達(dá)式，可以直接證明坐標(biāo)變換系數(shù)在通量分裂以后也是線性的，算子分別作用后

(15)

同樣使用算子δ1離散源項(xiàng)：

(16)

得到半離散形式：

(17)

(18)

對(duì)于Steger和Van Leer通量分裂格式可以推導(dǎo)出如式(18)所示的離散方程的算子形式，由于AUSM、HLLC和Roe通量分裂格式是在半點(diǎn)處重構(gòu)出來的，直接推導(dǎo)難度較大，采用數(shù)值實(shí)驗(yàn)驗(yàn)證，計(jì)算結(jié)果也沒有誤差。圖10和圖11是采用離散等價(jià)方程算法后第1節(jié)部分算例的計(jì)算結(jié)果，可以看出離散等價(jià)方程算法完全消除了坐標(biāo)變換誘導(dǎo)誤差。

圖10 柱坐標(biāo)系密度誤差云圖(一階迎風(fēng))

圖11 采用離散等價(jià)方程算法后的正弦網(wǎng)格密度誤差云圖

采用算子形式證明離散等價(jià)方程算法適用于任意格式，近期比較一階迎風(fēng)格式和五階ENO類格式[20]，計(jì)算結(jié)果和1.1節(jié)、1.2節(jié)定性一致，高精度格式的坐標(biāo)變換誘導(dǎo)誤差更大，通過分析誤差來源，得出的結(jié)論是ENO模板涉及的周圍相鄰網(wǎng)格點(diǎn)較多，在源項(xiàng)較大的局部區(qū)域引入的誤差較大，結(jié)論與本文一致。

式(11)～式(18)從帶有S≠0源項(xiàng)的離散等價(jià)方程算法推導(dǎo)中，算子作用在對(duì)流項(xiàng)守恒通量上，有些高精度格式采用原始變量重構(gòu)界面的形式，討論這種基于原始變量的離散等價(jià)方程算法實(shí)現(xiàn)。

5 非均勻網(wǎng)格界面通量重構(gòu)

以密度為例進(jìn)行討論，笛卡爾坐標(biāo)系下均勻網(wǎng)格的MUSCL格式為

(19)

(20)

對(duì)應(yīng)二階精度的3點(diǎn)中心格式：

(21)

根據(jù)Taylor級(jí)數(shù)：

(22)

(23)

式中：Δx1和Δx2分別為i點(diǎn)與i+1點(diǎn)網(wǎng)格間距和i點(diǎn)與i-1點(diǎn)網(wǎng)格間距；O(·)代表截?cái)嗾`差。

可以看出在非均勻網(wǎng)格Δx1≠Δx2條件下，應(yīng)用MUSCL的3點(diǎn)中心格式式(21)達(dá)不到二階計(jì)算精度，在該種網(wǎng)格下二階精度格式表示為

(24)

分析計(jì)算空間內(nèi)3點(diǎn)中心格式的精度：

(25)

不論坐標(biāo)變換函數(shù)如何變化，其實(shí)質(zhì)是物理空間到網(wǎng)格編碼的映射關(guān)系。以一維空間進(jìn)行簡(jiǎn)要說明，通過變換函數(shù)ξ=ξ(x)把物理空間x∈[a,b]映射到計(jì)算空間ξ∈[ξ(xa),ξ(xb)]，其中xa和xb分別為物理空間中a、b兩點(diǎn)坐標(biāo)，把ξ空間離散成N個(gè)網(wǎng)格點(diǎn)，格點(diǎn)間距為常數(shù)：

Δξ=[ξ(xa)-ξ(xb)]/(N-1)

(26)

計(jì)算中網(wǎng)格編碼i和ξ之間存在線性關(guān)系i=1+[ξ(xi)-ξ(xb)]/Δξ，常數(shù)Δξ合并至坐標(biāo)變換系數(shù)中，因此在離散后本質(zhì)還是物理空間x∈[a,b]到i∈[1,N]的映射關(guān)系。實(shí)際應(yīng)用中很少引入上述中間變換，而是直接從x∈[a,b]映射到ξ∈[1,N]，在計(jì)算空間離散點(diǎn)上滿足ξi=i。在離散后的計(jì)算空間上Δξ=1，決定計(jì)算空間大小的唯一參數(shù)是網(wǎng)格節(jié)點(diǎn)數(shù)N，如果網(wǎng)格加密過程中節(jié)點(diǎn)總數(shù)N保持不變，實(shí)際上只是調(diào)整離散點(diǎn)在物理空間[a,b]分布，與計(jì)算空間的映射關(guān)系沒有任何變化。如果通過增加網(wǎng)格總數(shù)N加密網(wǎng)格，即使物理區(qū)域[a,b]保持不變，其在計(jì)算空間的映射關(guān)系也發(fā)生相應(yīng)的變化。

根據(jù)以上分析可知通過式(25)難以直接得到截?cái)嗾`差量級(jí)，嚴(yán)格的精度分析還需要變換至物理空間，由Taylor級(jí)數(shù)可得

(27)

由于Ji+1同樣需要Taylor展開，具體推導(dǎo)過程非常復(fù)雜，因此針對(duì)一維空間變換函數(shù)ξ=ξ(x)進(jìn)行定性的量級(jí)分析。

根據(jù)微分關(guān)系式dξ=ξxdx可以得到坐標(biāo)變換系數(shù)的量級(jí)ξx～O(Δx-1)，進(jìn)一步推導(dǎo)出高階導(dǎo)數(shù)量級(jí)：

(28)

式中：f代表任意變量，用于任意物理量的高階導(dǎo)數(shù)量級(jí)分析。

在空間多維情況下，J至少與兩個(gè)變量有關(guān)，證明過程較為復(fù)雜，按照有限體積法的量級(jí)定性分析。MUSCL類型的格式應(yīng)用于有限體積法考慮了度量系數(shù)，常用形式為

(29)

基于以上分析，在坐標(biāo)變換后如果對(duì)原始變量ρ出發(fā)進(jìn)行重構(gòu)，應(yīng)該考慮相鄰網(wǎng)格尺度的變化，滿足精度要求的MUSCL格式為

(30)

對(duì)于均勻流場(chǎng)，在非均勻網(wǎng)格J不為常數(shù)的條件下，例如Ji≠Ji-1，其中Ji為節(jié)點(diǎn)i的雅可比行列式，即使?jié)M足ρi=ρi-1，也存在

(31)

由此推斷，重構(gòu)過程也會(huì)引起坐標(biāo)變換誘導(dǎo)誤差。

Euler方程是擬線性雙曲型方程組，對(duì)流項(xiàng)通量、守恒變量和系數(shù)矩陣之間具有性質(zhì)：

(32)

對(duì)物理空間自變量(x,y,z)求導(dǎo)：

(33)

同樣對(duì)計(jì)算空間變量(ξ、η、ζ)求導(dǎo)依然保持這種特性，以ξ為例：

(34)

對(duì)應(yīng)的算子形式有

(35)

代入第4節(jié)中帶有源項(xiàng)的離散等價(jià)方程出發(fā)得到算子形式中：

(36)

從式(32)～式(36)推導(dǎo)過程看，沒有任何精度損失。

為應(yīng)用通量分裂格式，同樣把系數(shù)矩陣作為常數(shù)移到空間算子內(nèi)，引入符號(hào)

(37)

表示離散點(diǎn)0及其相鄰網(wǎng)格點(diǎn)i的通量，用于編程的算子形式為

(38)

式(38)是在計(jì)算空間下的離散算子，在此基礎(chǔ)上針對(duì)U應(yīng)用MUSCL格式：

(39)

式(39)同直角坐標(biāo)系下均勻網(wǎng)格的重構(gòu)形式完全一致，但是二者的精度不同，式(39)達(dá)不到標(biāo)準(zhǔn)MUSCL格式的精度。根據(jù)上述分析，式(39)乘以坐標(biāo)變換系數(shù)Γ后轉(zhuǎn)換為計(jì)算空間的MUSCL格式才具有直角坐標(biāo)系下均勻網(wǎng)格的相應(yīng)精度，Γ本質(zhì)是非等距插值公式的影響系數(shù)。

對(duì)于均勻流場(chǎng)，不論參數(shù)κ如何取值，均滿足

(40)

因此必然不會(huì)產(chǎn)生坐標(biāo)變換誘導(dǎo)誤差。對(duì)于非均勻流場(chǎng)，計(jì)算誤差不僅與差分格式有關(guān)，還受到流場(chǎng)參數(shù)影響，如何驗(yàn)證還在探索之中。

6 結(jié) 論

通過針對(duì)不同網(wǎng)格的均勻流算例，展示了坐標(biāo)變換誘導(dǎo)誤差的普遍性，并對(duì)其產(chǎn)生機(jī)理進(jìn)行理論分析，討論了沿著GCL研究路線消除這類誤差存在的難度，提出基于離散等價(jià)方程的新思路，構(gòu)建離散對(duì)流項(xiàng)通量的計(jì)算方法，推導(dǎo)了針對(duì)MUSCL格式的應(yīng)用策略。提供的算例顯示離散等價(jià)方程方法能完全消除坐標(biāo)變換誘導(dǎo)誤差。