有關有限差分高精度格式兩個應用問題的討論

劉 君, 韓 芳

(大連理工大學 航空航天學院, 大連 116024)

0 引 言

模擬激波一直是計算流體力學(Computational Fluid Dynamics，CFD)的重要研究內容，Godunov、MUSCL、TVD、ENO、WENO、WCNS等差分格式，Steger、Van Leer、Roe、AUSM、HLLC等通量格式，這些被當做里程碑事件寫入CFD教科書的研究成果大多和激波有關。近些年本課題組提出了自適應間斷裝配求解器(Adaptive Discontinuity Fitting solver, ADFs)[1]，應用成熟的非結構動網格技術[2]很好地解決了因激波運動引起的計算域空間幾何變化的難題[3-5]，通過把R-H(Rankine-Hugoniot)關系式當作一種新的通量分裂形式在程序中實現了捕捉法和裝配法的自由轉換[6]，這種方法不僅編程方便，且應用簡單的結構辨識算法就能在計算得到的流場基礎上“逢山開路、遇水搭橋”很容易的完成復雜激波結構的自動裝配。

在高精度格式方面，基于體積平均思想的有限體積法(此處特指使用高斯定理把體積分轉換為面積分以后建立的算法)在構建高精度格式方面存在理論上的先天不足，國內外從事這方面算法和應用的研究者也較少[7-8]，而有限差分方法因理論完善、操作簡單而成為高精度格式發展的主要方向，這也是本文討論的范圍。搜索“高精度格式”文獻很多，遺憾的是至今還沒有哪種高精度格式能像當年TVD那樣在學術界得到普遍認可、在工程領域被廣泛使用，目前商業軟件的主流格式依然是二階精度。在本文中，作者針對有限差分的高精度格式在應用研究中所遇到的兩個問題進行分析討論，探索解決的途徑。每個問題論述結構相似，先介紹研究背景、數值現象和文獻綜述，然后分析討論和提出看法，最后是對解決問題的初步探索。

1 曲線坐標變換的誤差問題

裝配了激波以后流場空間形成互不重疊的相鄰區域，原則上可在多塊分區拼接網格上應用高精度有限差分法求解，但是由于激波在收斂過程中存在運動，因此，需要解決變形網格引起的幾何守恒律問題。

早期幾何守恒律用來消除因網格變形給流場引入的新誤差[9-10]，后來人們發現均勻流場在坐標變換后即使網格不發生變形也會出現誤差，為便于區分，把前者稱為體積守恒律(Volume Conservation Law，VCL)，后者稱為面積守恒律(Surface Conservation Law，SCL)。從文獻調研情況看，對體積守恒律的研究主要集中于有限體積法領域[2,10]；對有限差分法，日本的Abe、Nonomura[11-12]團隊分別于2012年、2013年提出了對稱守恒度量算法及非對稱守恒度量算法，對任意線性的高階有限差分格式具有很好的效果。近些年隨著高精度格式進入應用，對面積守恒律的研究越來越受關注[13-16]。基于變形網格的激波裝配法應用高精度格式，兩種幾何守恒律都需要注意，下面先討論靜止網格的面積守恒律。

關于曲線坐標變換以后出現的問題，文獻[17]在一階迎風格式的基礎上開展了研究，這里主要驗證其結論對高精度格式的有效性，因此，取同樣的算例進行對比，計算網格如圖1所示。

圖1 計算網格Fig.1 Computational grids

對任意坐標系下二維Euler方程：

(1)

通量分裂采用VanLeer格式，時間推進采用3階Runge-Kutta格式，按照CFL=0.9計算時間推進步長，計算時間t=2.0。通量導數計算分別選取五階WENO-JS[18]、WENO-M[19]、WENO-Z[20]、TENO[21]格式，首先計算進出口邊界條件確定的超聲速均勻流場，流場初場參數(ρ,u,v,p)=(1.4,2.0,0.0,1.0)，四種格式的最大殘差收斂曲線如圖2所示。圖3為采用不同格式計算得到的密度誤差分布云圖，從圖中可以看出密度誤差產生于網格非均勻區域，并沿著馬赫錐向下游傳播。圖4為流場參數(ρ,u,v,p)=(1.4,0.0,0.0,1.0)時的靜止流場的計算結果，可以看出，邊界附近并沒有更大的誤差。

圖2 最大殘差收斂曲線Fig.2 Maximum residual convergence curves for the uniform flow with different schemes

從幾何誘導均勻流場產生數值誤差的問題被發現開始，CFD研究者們對如何保持均勻流守恒性的研究便從未停止，其中發展較為成熟的是對幾何守恒律的研究。幾何守恒律認為原物理坐標系下的控制方程經坐標變換后在計算坐標系下的完整形式如式(2)所示，式(2)的右端項S=0在微分形式下自動成立，但經過有限差分離散之后一般不再成立；離散后的S不為0是導致均勻流不能保持守恒特性的主要原因。因此，幾何守恒律的研究主要是通過構造合適的坐標變換系數及坐標變換雅克比的離散格式使得S在離散形式下仍然能保持為0，以式(1)作為控制方程進行數值模擬，如現今被廣泛使用的由鄧小剛等研究學者提出的幾何守恒律的守恒計算方法(Conservative Metric Method, CMM)[13]及對稱守恒計算方法(Symmetrical Conservative Metric Method, SCMM)[14]。此外，Cai和Ladeinde[22]、Sun等[16]以式(2)作為控制方程在WENO格式上做了一系列工作，取得了很好的成果。本文的工作便是以式(2)作為控制方程進行開展的。為便于區分及表達，本文將式(1)稱為離散近似方程，將式(2)稱為離散等價方程。

(2)

式中：S=FIx+GIy。

作者在文獻[17]中提出了只要源項和對流項的離散格式相等就能滿足幾何守恒律的構造準則，得到方程(2)對應的半離散形式為：

(3)

其中：

由于采用離散等價方程和相容性準則消除坐標變換引起的誤差，和幾何守恒律采用構造坐標變換系數使得源項離散后仍為0的算法思路不同，文章發表后同行提出異議，因此在后續論文[23]中不再提幾何守恒律，而是采用“離散等價方程的相容性算法”來表述。文獻[23]采用算子形式證明了提出的離散等價方程相容性準則適用于任意格式，本文采用以上五階ENO類格式計算均勻流場，可以完全消除圖3和圖4中的誤差，進一步驗證了上述結論的正確性。

文獻[17]中給出了亞跨超聲速流場的算例，定性比較誤差云圖,和本文的高精度格式結果非常相似，因此本文不再給出更多馬赫數的計算結果，而是給出了采用一階迎風格式時的數值模擬結果(如圖5所示)作為對比。

對比圖3、圖4和圖5中相同來流條件下的計算結果，可發現在本文計算條件下，使用五階ENO類格式計算得到的誤差結果比使用一階迎風格式計算得到的誤差結果高一個量級，這是個還沒有看到文獻報道的新現象，下面進行初步分析。

(a) Ma=2.0

(b) Ma=0圖5 一階迎風格式密度誤差云圖Fig.5 Contours of density error for the uniform flow with the first-order upwind scheme

對于存在激波的流動，由于模板選擇與流動參數相關，給理論分析帶來難度，只能通過數值算例進行研究。對圖1所示計算網格，在x=0.1處放置一道Ma=3.0的正激波，給定初始計算條件為：

(ρ,u,v,p)=

(4)

圖6 通量計算選取計算點Fig.6 Calculation points selected for flux calculation

(a) 網格點 1

(b) 網格點 2圖7 不同計算方法得到的對流通量差值Fig.7 Errors of convex flux obtained by different methods

對流通量導數采用五階WENO-Z格式計算，時間步長參數選取CFL=0.5，經過時間t=0.25，激波運動至x=0.85的位置。比較數值解和理論值的計算結果，發現在激波附近密度誤差最大，因此，選擇0≤x≤0.5區域進行結果分析，密度誤差采用同樣幅度顯示，結果如圖8所示，其中圖8(a)是傳統CFD方法的計算結果，8(b)是采用離散等價方程的相容性算法得到的計算結果，圖8(c)為選擇同樣網格數目的正方形均勻網格的計算結果。

(a) 方法一 (b) 方法二 (c) 均勻網格圖8 x≤0.5區域的密度誤差云圖Fig.8 Contours of density error in the area of x≤0.5

對應云圖，圖9給出了t=0.25時刻沿著上下對稱線y=0.5的密度分布曲線。三種結果在激波下游都存在明顯的波動特性，從分布曲線及誤差云圖沿縱向的平行度來分析，采用離散等價方程的相容性算法后，誤差特性更接近沒有坐標變換的正方形均勻網格，表明這種算法對于存在激波的流場也是有效的。

圖9 y=0.5線上密度分布Fig.9 Density distributions along line y=0.5

圖10為計算區域中心點(x,y)=(0.5,0.5)處密度隨時間變化曲線，為了更好地分辨及分析計算結果，在其局部放大圖中給出一階迎風格式的計算結果。從圖中可以得到如下認識：

(1) 這種波動現象與坐標變換及格式精度無關，高精度格式的低耗散特性在維持波動向下游發展過程中比一階格式的幅值更高、寬度更窄；

(2) 波動從下游邊界離開后流場趨于穩定，收斂解存在可辨識的差異，齊次方程的高精度格式計算結果和一階迎風接近，曲線坐標離散等價方程和正方形網格的高精度格式計算結果接近。

經過以上數值試驗和理論分析，得出如下結論：對于采用模板擴展方式構造的高精度格式，在計算均勻流場時出現曲線坐標變換引起的誤差比低精度格式更大的現象，這是造成高精度格式在復雜工程應用中效果不佳的主要原因。從離散等價方程出發,采用相容性算法時使用高精度格式可以把曲線坐標變換引起的誤差降低到直角坐標系量級，有效解決了高精度格式的工程應用難題。

(a) 密度分布曲線

(b) 密度分布曲線的局部放大圖10 中心點處密度隨時間變化曲線及其局部放大Fig.10 Time-density curves at the central point and its partial enlarged drawing

2 激波捕捉法和Riemann分解的匹配問題

很多高精度格式的應用，例如旋渦主導的流動、湍流模擬、氣動聲學等，大多采用正方形均勻網格，以下算例均采用笛卡爾直角坐標系下Euler方程的直接離散，大部分算例是正方形網格，少量算例是長方形網格。

首先，采用正方形均勻網格考察式(4)所描述的流場，其對稱線上不同時刻密度和壓力的波動曲線如圖11和圖12所示。通過比較發現，密度曲線存在兩個逐步拉開的波動，而壓力曲線只有一個波動，且與密度曲線其中一個波動的位置相同。此外還給出一階迎風格式的計算結果，峰值基本一致，影響區域較寬，為便于分辨未標出。

選取時間t=0.1,0.2,0.4,0.6,0.8,1.0,1.25,1.5時的計算結果，畫出波1、波2的極小值位置隨時間變化曲線如圖13所示，其斜率正好是激波后流動的兩個特征速度u2和u2-a2，這個規律的發現有助于解釋波動現象的機理。

(a) t=0.25時波的位置示意圖

(b) 波1、波2速度分布圖13 波后擾動的傳播速度Fig.13 Propagation velocities of these two disturbance waves after the shock

計算網格采用正方形均勻網格時，由圖8(c)看出變量沿縱向沒有變化，因此對以上結果下面按照一維Euler方程進行理論分析。Euler方程不考慮黏性，激波是上下游參數滿足R-H關系式的數學間斷。20世紀50年代，Lax等提出了雙曲型守恒律方程的弱解理論，其中的Reimann分解理論能夠給出Euler方程在任意初始流場參數下的解析式。對CFD建立在離散空間網格上的求解模型，將相鄰網格之間物理量的差異當做間斷處理可以很好地利用弱解理論來構建算法。精確Reimann分解存在多種結構，計算較為耗時，因此后續發展出Steger、Van Leer、Roe、HLLC等通量格式來近似Reimann分解以提高計算效率。

激波捕捉法最早是von Neumann提出的，他認為激波具有有限的厚度，在離散空間可以由若干網格構成的過渡區描述，因此整個流場可以采用統一的計算格式。由于激波在理論上不存在過渡區，無法評價數值解的精度，只能是定性比較，認為過渡區涉及網格點越少或者曲線斜率較大的格式精度越高。而根據CFD理論，即使包含激波的流場，差分格式離散后的數值解也需要滿足如下修正方程：

(5)

以上兩個理論在激波處出現了矛盾：

(1) 由于每個網格點上只能存在一個數值解，利用捕捉法描述的間斷存在過渡區，而過渡區在數學上是不存在的，依然要進行Reimann分解，從而產生更多的非物理結構。高精度格式采用限制器抑制過渡區內的虛假結構不發散，因此，不同格式的效果不一樣(見圖3～圖5)。

(2) 初始時刻根據R-H關系式得到的激波本來是Reimann分解中的單一結構，由于不滿足以上修正方程在計算過程中也會被Reimann分解成多個結構。把初始激波分解出來的虛假結構當做輸入條件進行計算，不同格式精度的計算結果均會出現這種虛假波動現象(見圖10～圖12)。

很多時候網格和激波不能像以上算例保持完全匹配，為此，把以上正激波旋轉一定角度，旋轉后的激波與x軸的夾角記為α，計算初始條件為：

(ρ,u,v,p)=

(6)

對流通量導數仍然采用五階WENO-Z格式計算，通量分裂采用VanLeer格式，時間步長參數CFL=0.5，計算時間t=0.1。

選取的旋轉角度不同，計算結果也不相同，計算得到的密度誤差云圖和渦量云圖分別如圖14、圖15所示，可以看出小的角度變化會引起非常大的流場差異。同樣給出一階迎風格式的計算結果(圖16)作為對比，除了驗證以上分析結論，還發現高精度格式的高分辨率特性可以使這種虛假結構的結構更為清晰，在網格不匹配條件下相互干擾易出現虛假的“數值湍流”。

圖14 不同激波運動角度下的密度誤差云圖Fig.14 Contours of density error for the shock moving problem with different inclined angles

在激波捕捉法范圍內只能采用網格匹配來消除這種“數值湍流”，圖17是采用Δx≠Δy的長方形網格保證初始激波全部落在對角線網格點上的密度誤差云圖。不管定常或非定常問題，計算過程中激波都是運動的，因此要想實現激波與網格完全匹配只能采用激波裝配法。

關于裝配法的驗證和演示算例可以參考文獻[4-5]，本文給出以上正激波算例的裝配法計算結果來進行簡單比較。給定初始正激波位于x=0.4，時間步長參數CFL=0.5，運動時間t=0.1時計算停止，計算過程中網格變化如圖18所示。圖19為分別采用傳統CFD計算方法(方法1)和基于離散等價方程的相容性算法(方法2)得到的對稱線y=0.5上的密度分布曲線，因為裝配法直接采用R-H關系式計算激波，初始激波不會分解出虛假結構，激波后流動參數理論上為均勻流場的參數。圖19的計算結果也證明了這一觀點，采用方法1計算時，激波波后有因網格變換而產生的誤差，采用方法2后，誤差消除，激波波后為均勻流場。

圖15 不同激波運動角度下的渦量分布云圖Fig.15 Contours of vorticity for the shock moving problem with different inclined angles

圖16 一階迎風格式α=44.0°下時的密度誤差和渦量分布云圖Fig.16 Contours of density error and vorticity for the shock moving problem at α=44.0° with the first-order upwind scheme

圖17 網格匹配時計算得到的密度誤差云圖及 渦量分布云圖(α=44.0°)Fig.17 Contours of density error and vorticity for the shock matched grids (α=44.0°)

(a) 初始網格 (b) 計算終止時網格圖18 激波裝配計算過程中的網格變化Fig.18 Illustration diagrams of the grid changes during the shock-fitting process

圖19 y=0.5線上密度分布的局部放大圖Fig.19 Partial enlarged drawing of the density distribution curves alone the line y=0.5

3 結 論

(1) 有限差分格式應用于貼體坐標系時會因為不滿足幾何守恒律而產生數值誤差，本文通過數值實驗發現在同一計算條件下，對均勻流場來說，高階格式產生的誤差比一階迎風格式要大，經過初步分析認為因為高階格式選用的模板比一階格式的模板涉及到的網格點更多，每個網格點上的誤差疊加導致引入了更多的誤差。作者提出的基于離散等價方程的相容性算法在本文中被證明對高階格式同樣具有適用性，能夠消除高精度格式應用于貼體坐標系時產生的數值誤差。

(2) 本文在求解正方形均勻網格上的正激波運動問題時發現，在運動激波的波后會產生兩個新的波的結構，經過初步分析認為這是通量分裂格式在激波處產生的隨著特征線傳播的非物理波動。在二維問題中，高精度格式的高分辨特性導致在激波與網格不能完全匹配情況下多維波動相互干擾出現虛假的“數值湍流”現象，激波裝配法能夠很好地解決這個難題。