一種在網格內部捕捉間斷的Walsh 函數有限體積方法1)

任炯 王剛

(西北工業大學航空學院,西安 710072)

引言

在計算流體力學(computational fluid dynamics,CFD)領域中,對超聲速或高超聲速流場進行數值模擬時,經常會出現激波或剪切層(接觸面)等較為復雜的間斷流場結構.對于這些間斷結構,往往需要采用特殊的數值算法處理才能得到滿足精度或分辨率需求的模擬結果.目前通用的間斷處理方法主要有兩類:激波裝配法[1-6]和激波捕捉法[6-12].

激波裝配法的基本思想是:將激波作為流場區域里的運動邊界進行處理,利用Rankine Hugoniot(RH)條件精確地確定激波的運動速度和位置.該方法描述激波具有高精度和高分辨率的優點,缺點是算法復雜,尤其是流場中激波相互干擾作用,演化形成多尺度激波的復雜結構時,一般需要采用動網格技術才能完成[6]計算.而基于有限體積[6]或有限元[13-14]的激波捕捉方法在守恒型控制方程的框架下不需要對激波做過多特殊的處理,只需假設激波間斷位于流場網格離散后的控制體交界面處,根據交界面兩側的流場狀態量近似求解Riemann 問題,便可自動模擬出激波的位置和強度.近似求解Riemann 問題時通常采用的方法有通量矢量分裂格式[15-17]或通量差分分裂格式[18-22].相對于激波裝配方法,基于有限體積[6]或有限元[13-14]的激波捕捉方法算法簡單,而且可以在流場的光滑區域和間斷區域采用統一的格式計算,大大提高了計算效率.另外,該類方法還具有良好的幾何外形適應性.目前,它們已經發展成為CFD 領域中有效且主流的間斷數值模擬技術之一.然而,要達到與激波裝配方法相同的流場模擬精度和分辨率,其必須使用巨大的網格量或采用高階的流場重構技術[23-33].

通常,“間斷位于控制體交界面處”的假設使激波捕捉方法所引入的激波描述誤差較大,尤其在質量較差的粗網格上更為明顯,一般表現出鋸齒狀的激波面.為了突破以上假設的限制,能夠在“控制體內部捕捉到激波”的算法研究已經成為目前CFD 技術發展的目標之一.

國際上,已有一些學者進行了相關的研究,如Persson 和Peraire[34]在高階間斷Galerkin(discontinuous Galerkin,DG)方法的背景下,通過給原始控制方程增加人工黏性項將激波(間斷)尺度拉伸,使得可以利用高階連續近似的方法模擬黏性激波中包含的大梯度區域.他們的研究發現,p階多項式的連續近似只需要O(h/p)(h為網格尺度)的黏性就能分辨出激波,得到比控制體單元尺度更小的激波結構,由此達到子單元尺度分辨率的激波捕捉效果(即在控制體內捕捉激波的分辨率).不過該方法在引入數值或物理黏性模型的同時會增加需要經驗調節的參數,且對于不同的計算問題需要不同的黏性模型才能得到穩定的結果.此外,該方法只有在近似多項式的階次很高(≥8 階)時,才具有網格內部的間斷捕捉能力,對于簡單的常數/線性近似,其結果甚至不及原始方法的分辨率.

美國學者Gnoffo[35-36]將一組具有正交性的間斷Walsh 基函數[37-38]引入到CFD 領域,利用Walsh 基函數求解具有間斷解的非線性微分方程[35].其基本方法是根據Walsh 基函數的乘法封閉性和一些基本變換算子(如倒數變換、積分變換、微分變換等)將原始的微分方程變換到波數空間的多項式方程系統,在波數空間上求解原微分方程的全局解,以期得到比傳統有限體積方法更好的間斷捕捉效果.然而,在求解一維Burgers 方程和Euler 方程的算例中,該方法并沒有達到預期的理想效果.主要表現在兩個方面:其一是求解非線性偏微分方程時,各個原始變量的Walsh 基函數級數近似會引入很多復雜的多重交叉相乘項,導致最終得到的代數方程系統過于復雜,且邊界條件的處理也比較困難;其二是該全局解方法完全脫離了有限體積離散的框架,沒有明確的網格劃分過程,只針對最終的代數方程系統進行數值迭代求解,且采用數百個Walsh 基函數得到的數值結果與傳統采用迎風格式的有限體積方法的結果相比,也并未表現出明顯的優勢.

針對Persson 和Peraire[34]與Gnoffo[35]所提出的方法中存在的問題,本文提出了一種新型有限體積方法,該方法在保持有限體積離散優勢的基礎上,充分利用了Walsh 基函數的間斷特點,且能達到以下兩點目的:(1)避免Gnoffo[35]所提出的全局解方法的復雜性和數值模擬的局限性;(2)僅采用低階近似(如常數近似或線性近似)便可在網格內部捕捉到激波間斷,達到子單元尺度的激波分辨率.本文方法的基本原理是:摒棄傳統有限體積和有限元方法中流場變量在網格內部連續的假設,將守恒變量表示為截斷Walsh 基函數的級數形式.借助Walsh 基函數本身所具有的間斷方波屬性,在控制體內部配置間斷.所配置的間斷自然地將控制體虛分成若干個連續的子單元,將得到的關于Walsh 基函數級數的守恒型控制方程系統在這些子單元上積分離散,得到了一種新的求解方法.該方法對變量的基函數表示與有限元思想類似,且最終得到的離散系統也與有限元方法相似.但在子單元上的積分并沒有變分過程,而與有限體積方法完全相同.因此該方法雖然在思想上介于有限元和有限體積方法之間,但本質上屬于有限體積的離散框架.本文稱該方法為“Walsh 函數的有限體積方法(finite volume method with Walsh basis functions,FVM-WBF)”.僅依賴Walsh 基函數組離散的FVM-WBF 方法在子單元尺度上只有1 階計算精度.為了進一步提高對光滑解的分辨率,采用子單元上的流場變量均值進行流場的線性重構,將FVMWBF 方法的離散精度提升到2 階.在數值算例部分,運用1 階和2 階FVM-WBF 方法求解了無黏Burgers方程和Euler 方程的多個定解問題,以期驗證FVMWBF 方法的計算精度、計算效率、激波捕捉能力和魯棒性.

1 Walsh 基函數

Walsh 基函數由一類具有間斷特點的分段常數函數構成.在給定區域上,每個Walsh 基函數的圖像都由若干個相等幅度的方波構成.

若采用gn(x)表示閉區間[xa,xb]上具有n個波段的Walsh 基函數,則有如下計算公式[35]

其中

則稱gn(x)屬于p組.一般gn(x)的最小波段長度定義為

且p組Walsh 基函數的波段長度最多包含兩種,即dxp或2dxp.因此式(2)中(i)的取值為1 或2.另外,Walsh 基函數還具有正交性等優良數學特征,詳細描述可參考文獻[35],此處不贅述.

圖1 展示的是前4 組(0≤p≤3)Walsh 基函數(g1(x)～g8(x))的函數圖像,可以看到,各函數均由若干個等幅方波構成.

圖1 閉區間0≤x≤1 上, p=0,1,2 和3 組的Walsh 函數示意圖: g1(x)～g8(x)Fig.1 Walsh functions for p=0,1,2 and 3 in the interval 0≤x≤1: g1(x)～g8(x)

Gnoffo[35]通過研究Walsh 基函數級數對若干連續和間斷函數的擬合特點提出了求解非線性偏微分方程的Walsh 基函數全局解方法.該方法采用時空耦合算法,將微分方程的各變量表示成關于時間和空間的二元Walsh 基函數級數形式,并直接對微分方程做Walsh 基函數級數的積分或微分變換處理,最終得到關于基函數系數的復雜非線性代數方程系統,需要通過Newton 迭代[35]等數值手段才能完成計算.雖然Walsh 基函數全局解方法思路簡單,但存在一定的局限性,如:(1)求解不同的微分方程需要重新推導代數方程組,因此不具有普適性;(2)當微分方程的初邊值條件復雜時,其基函數級數展開可能導致最終的代數方程系統求解時收斂困難;(3)文獻[35]的測試算例顯示要得到與傳統有限體積方法類似的模擬效果,微分方程變量的每個維度采用的基函數個數一般要達到28量級左右,這將會給最終非線性代數方程系統的求解帶來較高的計算量負擔.

受Gnoffo[35]研究的啟發,本文將具有間斷特性的Walsh 基函數與傳統有限體積方法結合,提出了新的FVM-WBF 方法,以期在一定程度上改善Walsh 基函數全局解方法的局限性.

2 FVM-WBF 方法

根據Walsh 基函數級數對任意函數擬合的合理性[35],本文考慮采用截斷Walsh 基函數級數近似雙曲守恒律方程中的守恒變量,下面針對一維和二維情況分別進行討論.

2.1 一維FVM-WBF 方法

為描述簡便,我們先以標量守恒律方程為例

其中u=u(x,t)為守恒變量,f(u)=f(u(x,t))為通量函數.若采用傳統有限體積方法對守恒方程(5)進行離散,需將其在給定控制體?i上進行積分,并假設控制體單元內的守恒變量為連續分布函數,從而得到關于守恒變量單元平均值的半離散方程

現對守恒變量u進行截斷Walsh 基函數級數近似如下

式中ck(t)是隨時間變化的Walsh 基函數系數,且n=2p(p=0,1,2,···).將式(7)代入式(5),則得到關于系數組和Walsh 基函數組的守恒型方程組

圖2 控制體?i 內的變量分布示意圖Fig.2 Variable distribution in the control volume ?i

經整理,該式可寫成

根據式(7),不難發現n(n=2p,p=0,1,2,···)個Walsh 基函數的級數表達形式將會給控制體?i內的守恒變量近似分布中均勻引入n?1 個間斷.若取n=4,那么控制體?i內的3 個間斷位置如圖2 中的灰色豎直虛線所示.沿著n?1 個間斷所在位置進行劃分,控制體?i將被均勻地分割為n個大小相同的子單元,用?i,j(j=1,2,···,n)表示.對于一維情況,?i,j=[xi,j?1,xi,j],將關于Walsh 基函數系數組的守恒型方程(9)在子單元?i,j上積分,有

若定義子單元積分平均為

其中dxp=xi,j?xi,j?1,則方程(10)可寫為如下半離散形式

另外,由式(7)還可知子單元內守恒變量為常數分布近似,所以在整個控制體?i內,守恒變量的分布具有分段常數的特點,如圖2(c)所示.

根據式(11),對方程(10)進一步整理,可得到方程(12)的矢量形式

若將守恒型方程(9)在所有子單元上積分離散,則可在控制體?i上得到關于Walsh 基函數系數矢量C的緊致矩陣方程

該矩陣中元素的計算公式為矩陣A(n)表示n個Walsh 基函數對應得到的矩陣,其中所有元素的公因子為,所以該矩陣可歸一化為只包含±1 的常數矩陣.例如,當n=2,4,8 時有

顯然,矩陣A(n)對稱可逆,因此方程(15)的“解的存在性”得到保證.一個有趣的特點是矩陣A(n)的逆A(n)?1可由A(n)直接計算得到,即

根據該特性,在求解時可以省略矩陣求逆的過程,直接由預先儲存的A(n)進行計算,這在一定程度上對計算效率是有益的.

以上介紹的求解標量方程的FVM-WBF 方法同樣適用于求解守恒型方程系統,例如,對于一維Euler方程

其中ρ,u,p,E分別為密度、速度、壓強和單位體積總能,m=ρu,狀態方程為

γ=1.4 為理想氣體常數.

類似于標量情況,對3 個守恒變量(密度ρ,動量m和能量E)進行截斷Walsh 基函數級數近似

2.2 二維FVM-WBF 方法

對于二維情況,同樣先考慮標量守恒律方程

其中u=u(x,y,t)是守恒變量,f(u)和g(u)分別是x和y方向上的通量函數,(x,y)∈?(? 為求解域).

在控制體?i上,對守恒變量u=u(x,y,t)進行如下Walsh 基函數級數近似

其中,n1和n2分別為x和y方向上各自使用的Walsh基函數個數,ck1,k2(t)表示對應的基函數系數.

顯然,在控制體?i內,式(21)會給守恒變量的近似分布函數中引入n1+n2?2 個間斷(其中x方向為n1?1 個間斷,y方向為n2?1 個間斷).沿著間斷,可將原始控制體?i分割為n1·n2個大小完全相同的子單元,同樣用?i,j(j=1,2,···,n1·n2)表示.與一維方法類似,將式(21)代入守恒型方程(20)后,在每個子單元上進行積分離散.最終得到關于Walsh 基函數系數組

的緊致型矩陣方程

其中C=[c1,1(t)···c1,n2(t)c2,1(t)···cn1,n2(t)]T,F=[Fi,1Fi,2···Fi,n1·n2]T,(Fi,j(j=1,2,···,n1·n2)為?i,j上x方向的通量),G=[Gi,1Gi,2···Gi,n1·n2]T,(Gi,j(j=1,2,···,n1·n2)為?i,j上y方向的通量).對于二維情況,矩陣A(n1,n2)可由相應的2 個一維矩陣的張量積得到,即A(n1,n2)=A(n1)?A(n2)(?表示張量積),也可由相應的元素公式計算得到

該矩陣同樣對稱可逆,所以對于方程(22),“解的存在性”可以得到保證.而且A(n1,n2)的逆也可由其本身計算得到,即

(p1和p2分別是n1和n2對應的Walsh 基函數組數)或者由對應的兩個一維逆矩陣的張量積得到,即

實際求解時,矩陣A(n1,n2)及其逆都可預先計算并儲存,以減少計算過程中矩陣求逆帶來的計算量.

對于二維守恒型方程系統,將守恒變量按分量進行Walsh 基函數級數近似,可得到與方程(19)類似的塊對角緊致矩陣方程,對角線上的每個矩陣均為A(n1,n2),其推導過程類似一維情況,這里不再贅述.

由上述描述可知,FVM-WBF 方法采用時空分離的思想在子單元框架下對微分方程進行有限體積離散,具有更好的普適性.且FVM-WBF 方法是根據Walsh 基函數引入的間斷位置對原始控制體單元進行虛子單元劃分,并不改變原始網格的幾何拓撲結構,保持了原始網格良好的幾何適應性.另外,對于邊界條件的處理,FVM-WBF 方法與原始有限體積方法類似,為了保持與內場一致的數值通量計算,可以在邊界處延拓一層輔助網格單元,用與邊界單元相同的Walsh 基函數表達輔助網格上的流場變量,可以在邊界處達到子單元尺度的誤差精度.相對于Gnoffo[35]提出的Walsh 基函數全局解方法,FVM-WBF 方法并未在離散和求解方面引入新的困難.

2.3 誤差和精度分析

傳統有限體積方法假設控制體內的守恒變量為常數函數分布,如圖2(a)所示.那么相對于每個控制體?i的尺度,計算精度為1 階,其誤差為O(?x),即

而對于FVM-WBF 方法,守恒變量采用Walsh 基函數級數近似(如式(7)),控制體內守恒變量的分布為分段常數函數形式(見圖2(c)).相對于每個子單元?i,j的尺度,計算精度也為1 階,子單元上的誤差為O(dxi,p),即

而對于整個控制體?i的尺度,Walsh 基函數級數近似的積分平均如下

當n=2 時,

由子單元積分平均式(11),有

再由式(26)和式(27),可得

所以

同理,

當n=4 時,有

當n=8 時,有從以上推導可知,對于FVM-WBF 方法,當采用n=2p(p≥1)個Walsh 基函數對守恒變量近似時,其近似誤差為O(?x/2p).因此理論上,相對于傳統有限體積方法,FVM-WBF (n=2p)方法的誤差約以因子1/2p縮減,而收斂精度并沒有提高.

通常,為提高精度,傳統有限體積方法是在當前控制體和周圍多個控制體構成的模板上,采用守恒變量均值進行多項式重構,如線性重構后,變量的分布如圖2(b)所示.而本文根據FVM-WBF 方法的特點,采用控制體內的子單元均值(如式(11))進行線性重構實現2 階精度的計算,具體公式如下

其中

線性重構后,控制體?i內的變量分布如圖2(d)所示,為分段線性分布.傳統2 階有限體積方法與本文重構的2 階FVM-WBF 方法在計算精度和計算效率方面的比較將在下小節的數值測試中進行討論.

3 數值測試

本節采用FVM-WBF 方法對一/二維無黏Burgers方程和Euler 方程的幾個非定常問題進行求解,并通過與傳統有限體積方法的數值結果進行對比,測試和驗證該方法的計算精度、計算效率、激波捕捉能力和魯棒性.其中數值通量的計算均采用熵相容格式[21-22],時間推進均采用滿足TVD 條件的3 階顯式Runge-Kutta[39]方法.

為便于表述,算例中采用記號FVM-WBF-O 表示初始沒有線性重構的1 階FVM-WBF 方法的數值結果,而用FVM-WBF-R2 表示經線性重構實現的2 階FVM-WBF 方法的數值結果,作為對照,記號OFVM-1st 和OFVM-2nd 分別用來表示傳統有限體積方法(用OFVM 方法表示)的1 階和2 階計算結果.一維算例中的Exact 表示解析解.

3.1 一維無黏Burgers 方程連續初值問題

Burgers 方程的第 1 個測試算例是一維區域[?2,2]上的連續初值問題,主要用以測試FVM-WBF方法的計算精度和計算效率,其初始條件為

解析解為

取u0=0,u1=0.5,采用周期邊界條件.針對t=0.32和t=0.96 這兩個時刻,分別在表1 到表4 中給出了不同網格上的OFVM 和FVM-WBF (n=2,4,8)方法1 階與2 階數值計算的L1誤差與相應的收斂精度.可以發現,對于兩種計算精度,OFVM 與FVMWBF 方法在兩個時間點的收斂精度基本相同,但在相同網格上,相對于OFVM 方法,FVM-WBF 方法的誤差約以因子1/2p縮減(這與第2.3 節分析的結論一致).若OFVM 方法采用的網格單元數目等于FVMWBF 方法的網格單元數與基函數個數的乘積時,在整個計算域上兩種方法具有相同的自由度.對比表1 ～表4 中兩種方法自由度相同時的誤差,可以發現,1 階精度時兩種方法的誤差幾乎相同,這說明FVM-WBF 方法表現出了與OFVM 方法網格加密等價的誤差效應.而2 階精度時,FVM-WBF 比OFVM方法的誤差略小,這可能是由于子單元重構模板與原始控制體重構模板的差異造成的.另外,表2 和表4 給出的t=0.96 時刻的誤差和收斂精度相對于t=0.32 時刻都有所減小,這是因為t=0.96 時刻在x=0位置處形成的壓縮激波對精度造成了影響.圖3(a)和圖3(b)分別給出了網格量為M=20 時,兩個時間點上兩種方法的2 階精度結果與精確解的對比,可以看到FVM-WBF 方法取n=2,4,8,3 種基函數級數近似時,在光滑解的極值附近和間斷解處都表現出比OFVM 方法更好的分辨率.圖3(c)展示的是兩種方法取相同的自由度時(本算例的自由度為1280),其誤差精度與所消耗的CPU 時間之間的關系,通過比較發現,當誤差和收斂精度相同時,FVM-WBF 方法所消耗的CPU 時間都比OFVM 方法短,且隨著Walsh 基函數數目的增大,FVM-WBF 方法的計算效率也逐步提高.

表1 OFVM 方法與FVM-WBF-O 的L1 誤差和收斂精度(t=0.32)Table 1 L1-error and L1-order of OFVM and FVM-WBF-O method for t=0.32

表2 OFVM 方法與FVM-WBF-O 的L1 誤差和收斂精度(t=0.96)Table 2 L1-error and L1-order of OFVM and FVM-WBF-O method for t=0.96

表3 OFVM-2nd 方法與FVM-WBF-R2 的L1 誤差和收斂精度(t=0.32)Table 3 L1-error andL1-order of OFVM-2nd and FVM-WBF-R2 method for t=0.32

表4 OFVM-2nd 方法與FVM-WBF-R2 的L1 誤差和收斂精度(t=0.96)Table 4 L1-error and L1-order of OFVM-2nd and FVM-WBF-R2 method for t=0.96

3.2 一維Euler 方程激波管問題

本算例求解了一維Euler 方程的兩種激波管問題(SOD 激波管和LAX 激波管).在相同網格上比較了FVM-WBF 方法和OFVM 方法對弱間斷和強間斷的捕捉效果.兩個激波管問題的計算域均取為[?0.5,0.5],初始條件分別為:

(1)SOD

圖3 一維Burgers 方程連續初值問題數值解Fig.3 Numerical results of continuous initial value problem of 1D Burgers equation

(2)LAX

均采用齊次Neumann 邊界條件,網格量均為M=100.SOD 激波管問題計算到時間t=0.1,LAX 激波管問題計算到時間t=0.16.OFVM 方法與FVM-WBF方法對兩個激波管問題的1 階和2 階計算所得數值密度的分布分別如圖4 和圖5 所示.通過對比激波、接觸間斷和稀疏波拐角等處的數值結果,可以發現,網格相同的情況下,FVM-WBF 方法對弱間斷和強間斷均有更好的捕捉效果,且對光滑解的模擬也有一定的改善.圖4(c)和圖5(c)給出的是兩種方法的2 階計算捕捉到的激波與網格尺度對比的放大圖,其中灰色虛直線代表網格邊界.可以看到一維FVM-WBF方法在n=4 和8 時均能達到小于網格尺度的激波捕捉效果.

3.3 二維無黏Burgers 方程間斷初值問題

圖4 一維Euler 方程SOD 激波管問題數值解Fig.4 Numerical results of SOD-shock tube problem of 1D Euler equations

圖5 一維Euler 方程LAX 激波管問題數值解Fig.5 Numerical results of LAX-shock tube problem of 1D Euler equations

本算例求解了二維無黏Burgers 方程,比較了二維FVM-WBF 方法與OFVM 方法在相同網格上模模擬標量間斷初值問題的效果.計算域為[?1,1]×[?1,1],初始條件為采用周期邊界條件,計算到時間t=0.3,網格量為M=20×20.圖6 展示的是兩種方法2 階精度的數值結果.其中,圖6(a)和圖6(b)是等值線圖.在這兩幅圖中,同時給出了網格線以清楚地反映所得激波尺度與網格尺度的對比效果.圖6(c)進一步給出了x=y截線上的數值解分布(圖中灰色虛直線代表網格邊界),可以發現FVM-WBF 方法中的n1和n2取2,4,8 時捕捉到的激波尺度都小于網格尺度,說明FVM-WBF方法在二維標量情況下能夠達到網格內部捕捉間斷的分辨率.

圖6 二維Burgers 方程間斷初值問題數值結果Fig.6 Numerical results of the discontinuously initial value problem of 2D Burgers equation

3.4 二維Euler 方程黎曼問題

二維Euler 方程的黎曼問題反映了多變氣體的多尺度流場結構,一般由3 種常見的一維波(稀疏波,激波和接觸間斷波)相互作用形成.本算例在計算域[0,1]×[0,1]上考慮如下初始條件的黎曼問題[40]

采用齊次Neumann 邊界條件,在相同的網格M=40 × 40 上比較了OFVM 方法與FVM-WBF 方法計算到t=0.3 時刻對激波的捕捉效果.圖7(a)和圖7(b)分別展示的是兩種方法數值結果的密度等值線圖.可以看到,相對于OFVM 方法,FVM-WBF 方法顯著提高了激波間斷的分辨率.取y=0.9 的截線,如圖7(c)所示,可以發現FVM-WBF-R2 在n1,n2分別均取4 或8 時都能達到小于網格尺度的激波捕捉分辨率,而OFVM 方法捕捉的激波約為2 個網格尺度的寬度.該算例進一步說明FVM-WBF 方法在二維Euler 方程系統下也可以實現網格內部捕捉間斷的效果.

圖7 OFVM 和FVM-WBF 方法對二維Euler 方程黎曼問題的2 階精度數值密度結果Fig.7 Numerical density obtained by 2nd-order OFVM and FVM-WBF method for Riemann problem of 2D Euler equations

3.5 Rayleigh-Tayler 不穩定性問題

Rayleigh-Tayler (R-T)不穩定性問題[41]描述的是,對于上下兩層密度不同的流體,當密度大的流體受到重力影響流向密度小的流體時所產生的流動現象.經過一定時間的演化,輕流體的空泡和重流體的“尖頂”相互干擾,最終形成非常復雜的流場結構.

該算例的計算域為[0,0.25]×[0,1],初始條件為

另外,該算例還考察了兩種方法計算相同時間步數時的計算效率.因為網格量相同時,FVM-WBF方法的自由度是OFVM 方法的n1·n2倍,所以理論上,計算相同的時間步數,FVM-WBF 方法消耗的CPU 時間應該是OFVM 方法的n1·n2倍,而自由度數目相同時,兩者所耗費的時間應該相同,但測試結果顯示并非如此.表5 給出了1000 個時間步時兩種方法的CPU 時間對比.對于3 種基函數個數,OFVM 方法采用了與FVM-WBF 方法相同自由度數目的網格量.可以看到,自由度相同時,平均每個時間步FVMWBF 方法比OFVM 方法提高了約20%的計算效率,這可能是由于FVM-WBF 方法程序實現時數據結構存儲方面帶來的效益.未來可以從自適應算法的角度進一步提高FVM-WBF 方法的計算效率,相關的研究工作正在進行中.

圖8 OFVM 和FVM-WBF 方法對R-T 不穩定性問題的2 階精度數值密度等值線圖Fig.8 Density contours obtained by 2nd-order OFVM and FVM-WBF method for R-T problem

表5 FVM-WBF 方法與OFVM 方法計算效率比較Table 5 Comparison of computational costs between FVM-WBF and OFVM methods

4 結論

本文借助Walsh 基函數發展了一種可以在網格內部捕捉間斷的新型有限體積方法,簡稱FVM-WBF方法.其基本原理是利用具有方波特性的Walsh 基函數在控制體內離散表示守恒變量,從而將間斷配置在控制體內部.沿著間斷位置將控制體劃分為若干個連續的子單元,通過數值積分在子單元上對守恒型控制方程進行離散,求解以Walsh 基函數系數為解變量的離散控制方程,從而得到Walsh 基函數級數形式描述的數值解.文中理論推導和算例測試表明了FVM-WBF 方法具有高分辨率、高效率和良好的激波捕捉能力等優點,主要表現在:

(1)當FVM-WBF 方法與傳統有限體積方法具有相同的收斂精度時,其誤差縮減為傳統有限體積方法的1/2p倍(p為所采用的Walsh 基函數級數截斷對應的組數).因此,隨著Walsh 基函數組數p的增大,FVM-WBF 方法的誤差將以1/2 的p次方倍迅速下降,相應數值結果的分辨率也顯著提高.

(2)當FVM-WBF 方法與傳統有限體積方法具有相同的自由度時,收斂到相同的誤差精度或計算相同的時間步數,FVM-WBF 方法都表現出更高的計算效率.

(3)FVM-WBF 方法在控制體內的子單元劃分為有限體積高階重構提供了新的可能的模板選擇.本文采用子單元均值進行線性重構,得到了2 階收斂精度的FVM-WBF 方法.與傳統有限體積方法的2 階計算結果相比,FVM-WBF 方法的2 階計算依然表現出更高的激波捕捉能力,采用較少的Walsh 基函數(如2個或4 個)就能達到小于網格尺度的激波分辨率.

(4)FVM-WBF 方法在子單元上的積分離散與傳統有限體積方法相同,所以傳統有限體積方法中發展的各種成熟算法都可應用于FVM-WBF 方法.

另一方面,從文中的推導可以看到,將Walsh 基函數系數的守恒型控制方程在子單元上積分離散后,會得到關于Walsh 基函數系數的緊致型矩陣方程系統.因此,與通過擴展模板進行高階重構實現高分辨的傳統有限體積方法相比,FVM-WBF 方法更適用于并行方法的計算.而且隨著基函數組數的增大,控制體內劃分的子單元數目相應地按倍數增多,這也使得FVM-WBF 方法天然具有網格自適應和多重網格的潛力.這些特點應得到更為深入的研究,以期為CFD算法的探索和應用提供新的思路.