基于混合重構高階間斷伽遼金方法的二維層流和湍流數值模擬

熊 為,張衛國,丁 玨,楊小權,翁培奮

(1.上海大學力學與工程科學學院,上海 200444;2.中國空氣動力研究與發展中心旋翼空氣動力學重點實驗室,四川 綿陽 621000)

目前,基于二階精度數值格式求解雷諾平均Navier-Stokes(Reynolds-averaged Navier-Stokes,RANS)方程的方法得到了廣泛應用,定常湍流流動問題已經得到了很好的數值解決.但是,基于高精度精細化的湍流數值模擬依然發展緩慢.高精度有限元方法以間斷伽遼金(discontinuous Galerkin,DG)有限元方法為代表,構造單元間分段的多項式解函數,僅增加解函數的自由度就能提高求解精度,因此具有較大的發展前景.相對于有限體積方法,基于有限元方法的間斷伽遼金方法本身存在不足,即處理相同的網格數量時,間斷伽遼金方法需要求解一個含有更多變量的矩陣方程組,從而導致間斷伽遼金方法在進行流場計算時通常需要消耗大量的計算資源和內存資源.

已有很多研究致力于降低間斷伽遼金方法在計算資源方面的消耗.Dumbser等[1-2]引進了一系列重構間斷伽遼金(reconstructed discontinuous Galerkin,rDG)方法的思路,以PnPm格式為標志提出了第一個重構間斷伽遼金方法,即最小二乘復原方法(least-square recovery method),其中Pn代表用于估計間斷伽遼金方法所求守恒量的解時,使用分段多項式的階數為n;而Pm則代表重構的分段多項式的階數為m.Dumbser等[1-2]的工作是依據文獻[3-4]提出的在單元內進行復原(recovery)的思路展開的.在使用最小二乘復原方法獲得多項式解時,認為重構多項式與相關控制體單元里間斷伽遼金方法直接獲得的解在弱形式上相同.最小二乘復原方法不僅對于單元內均值,而且對于單元內所有的高階矩都保證了較高程度的守恒性.然而,高守恒性的最小二乘復原方法在計算資源方面消耗比較高,同時重構過程依賴相鄰網格信息,在處理邊界附近網格單元時,會由于相鄰單元數目不足而降低重構的精度.Luo等[5]和Xia等[6]提出了基于泰勒基函數,可在任意網格上用于求解可壓縮Euler方程和Navier-Stokes方程的最小二乘重構方法(least-square reconstructed method).這種重構方法對相鄰單元中心點的值和導數進行插值,從而得到新的多項式方程.這類方法只在von Neumann相鄰單元上進行重構,因此格式相對而言具有緊致性、簡單性、魯棒性和可塑性,并且能夠保證較為精確的變量守恒性,包括單元內均值和泰勒基函數導數的守恒性.然而,此類重構方法無法保證2-exact的特性.為此,Cheng等[7-8]提出了在任意網格上針對可壓縮流動的混合最小二乘重構方法(hybrid least-squares reconstructed methods).該方法綜合了最小二乘點值重構和最小二乘積分重構的優點,確保了計算精度滿足2-exact的特性.在高階重構方法中,黏性項的離散大多采用第二Bassi-Rebay(second Bassi-Rebay,BR2)格式,但是傳統方法沒有考慮BR2格式局部和全局提升算子高階重構量的貢獻,這勢必會對計算精度造成一定的損失.

基于上述討論,本工作在混合最小二乘重構方法的基礎上,發展三階重構的BR2格式黏性數值通量的重構間斷伽遼金方法,用于求解可壓縮Navier-Stokes(NS)方程和RANS-SA(Reynolds averaged Navier Stokes-Spalart Allmaras)系統方程.通過典型算例,著重分析混合重構的rDGP1P2方法的計算效率、計算精度和收斂性.

1 控制方程

為了提高計算精度以及數值穩定性,本工作根據文獻[9-10]將RANS方程和一方程Negative Spalart-Allmaras湍流模型方程耦合成為系統方程,

式中:守恒向量U、對流項通量F、黏性項通量G和源項S的表達式分別為

其中壓力p可以由狀態方程得到,

Γ取1.4.應力分量為

式中:δij是Kronecker符號;k是流體的熱傳導系數.

湍流模型中的生成項P和破壞項D分別為

渦黏性系數μT的定義為

修正的渦量~S為

其他參量請參考文獻[9-10].對于層流算例只求解NS方程,不考慮SA湍流模型.

2 重構的間斷伽遼金方法

2.1 間斷伽遼金方法

方程(1)弱積分下的DG方法可表示為

黏性通量離散采用BR2格式[11],該格式的形式如下:

為了保證格式的穩定性,η取1～6.r為局部算子,R為全局算子,定義

兩種算子的離散形式為

其中

式(12)中Uh可近似為

基函數Bi(x)采用Xia等[6]提出的Taylor基函數.

2.2 混合高階重構方法

在傳統DG方法中,方程(5)的各階導數是通過求解方程組(4)得到的.三階數值模擬由于矩陣規模巨大,導致計算時間長,且占用大量內存.為了減少計算量,進一步提高計算效率,本工作采用重構方法,即自由度中的平均值和一階導數通過傳統的DG方法求得,而二階導數則采用重構的方法獲得,同樣可以達到三階精度.

文獻[7-8,12]發展的混合最小二乘重構方法兼容了積分重構方法和點值重構的優點,即平均值的計算采用積分方程,其余采用點值導數.該方法的優勢是計算效率高于傳統DG方法且滿足2-exact特性,使得非均勻網格能夠達到三階精度.具體做法如下:

在混合最小二乘重構方法中,重構以后的UR不是間斷伽遼金方法中單元中心點的值,而是將單元內平均值作為重構間斷伽遼金方法解向量多項式的第一個自由度,而解向量根據相鄰單元內平均值進行重構.利用與控制體單元i擁有公共面的單元j,建立如下關系式:

式中:R為重構方程組右端項.方程(19)左乘系數矩陣的轉置,得到線性方程組

其中

線性方程組(20)直接求解即可,無需進行迭代求解.對于間斷伽遼金方法,如果進行最小二乘重構,需要對中心點進行插值.而混合最小二乘重構方法用單元內守恒量的平均值進行多項式重構,滿足2-exact特性,保證了求解守恒量方程的精度.

當采用BR2格式進行黏性項離散時,在傳統rDGP1P2的BR2計算中,局部提升算子和全局提升算子均采用P1分布.為了進一步提高計算精度,本工作對兩類提升算子均采用P2多項式分布.具體表達式如下:

結合對守恒量進行重構的方式,對BR2格式的提升算子進行重構,用低階導數項求解高階導數項.

2.3 隱式時間推進格式

方程(8)寫成半離散形式為

式中:M是質量矩陣;U是解矢量;R是右端殘值.R包括無黏通量、黏性通量和它們對應的域積分項,這里無黏通量的求解采用Harten-Lax-van Leer接觸(Harten-Lax-van Leer contact)格式[13],黏性通量的求解采用BR2格式[11].

方程(25)線化后可以表示為

式中:?Rn/?U是Jacobian矩陣;Δt是當地時間步長.當Δt趨于無窮大時,系統方程的求解就變成了典型的牛頓迭代.

在通常情況下,Jacobian矩陣的重構方法模板比較寬泛,很難精確求解.主要存在的問題是,重構的高階導數對應的貢獻很難準確得到.為了盡可能地提高Jacobian矩陣的精度,本工作在處理Jacobian矩陣時考慮BR2格式局部提升算子和全局提升算子中高階導數的貢獻部分.除此之外,對平均值和一階導數的Jacobian可精確求解(即P1的精確Jacobian和局部提升算子、全局提升算子對二階導數的貢獻).另外,壁面邊界、遠場邊界等邊界上對應無黏和黏性數值通量的Jacobian矩陣依然給出了P1的精確求解.

對于定常問題,方程(21)的求解采用基于上-下對稱高斯賽德爾格式(lower-upper symmetric Gauss-Seidel scheme,LU-SGS)預處理的廣義極小剩余(generalized minimal residual,GMRES)方法;對于非定常計算采用四步隱式Runge-Kutta方法[14]求解,其中隱式子迭代同樣采用基于LU-SGS預處理的GMRES方法進行求解,子迭代的設置為10步或殘值下降3個量級.

3 算例分析

3.1 二維平板層流流動

二維絕熱平板層流的自由來流馬赫數為Ma∞=0.2,雷諾數為Re=1×105.計算域為[-0.5,1.0]×[0,1.0],其中底部[0,1.0]為無滑移邊界條件,底部其余部分為滑移邊界條件;其他邊界為特征邊界.分別采用三角形(2 400個單元)、四邊形(1 200個單元)和三角形/四邊形混合(1 800個單元)網格進行計算,如圖1所示.

圖1 平板層流計算網格Fig.1 Grid of laminar flow over a flat plate

圖2給出了平板層流流動表面摩擦力的常用對數在x方向的分布.基于PnPm格式重構間斷伽遼金(rDG)方法,將rDGP1P2和DGP1方法的計算結果與Blasius近似解進行比較.表面摩擦力分布顯示,rDGP1P2方法的求解結果與Blasius解吻合很好;與DGP1方法相比,rDGP1P2方法的計算結果與Blasius解吻合得更好.圖3給出了平板層流流動在混合網格上的密度殘差收斂曲線.由圖可知,DGP1方法的收斂需要迭代25步,rDGP1P2方法的收斂需要迭代28步,可見rDGP1P2方法在提升求解的精度基礎上,迭代步數和DGP1相當.

圖2 平板層流流動計算表面摩擦力對數分布Fig.2 Logarithmic skin friction coefficient of laminar flow over a flat plate

圖3 平板層流流動計算密度殘差曲線(混合網格)Fig.3 Convergence histories for density of laminar flow over a flat plate(hybrid mesh)

與二階精度的DGP1方法相比,三階DGP2方法提高了精度階數,守恒量多項式的自由度從3個增加到6個,對應求解時每個單元的自由度由原來的3×3(一階三角形網格)或4×3(一階四邊形網格)分別變為3×6和4×6.理論上,應用傳統DG方法時,三階精度模擬比二階精度模擬的計算量增大一倍,每個時間步長消耗的平均CPU時間也增加一倍.在DGP1(二階)基礎上重構的P2格式,即DGP1P2具有三階精度,但是每個單元上求解的自由度只有3個,因此計算量大幅度減小.表1給出了二維平板層流數值算例下rDGP1P2方法和DGP1方法的單步長平均CPU時間.由表可知,rDGP1P2方法每個時間步長消耗的平均CPU時間與DGP1方法相當,且與網格數、網格類型或控制方程類型均無關.通過混合重構,具有三階精度的rDGP1P2額外的3個自由度可利用二階(P1)精度通過求解一個較小的矩陣得到,在提高計算精度的同時,提高了高精度數值方法的計算效率.

表1 二維平板層流流動單步長平均CPU時間(單位:秒)Table 1 Average CPU time of each step of laminar flow over the plate(Unit:s)

3.2 Bump湍流流動

二維Bump湍流的自由來流馬赫數為Ma∞=0.5,雷諾數為Re=3×106.計算域為[-25.0,25.0]×[0,5.0]的矩形通道,底部[0,1.5]為帶有對稱凸起的無滑移壁面邊界條件,上下面其余部分為滑移壁面邊界條件,入流和出流為特征邊界條件.計算中分別采用2種四邊形網格進行數值計算,粗網格和細網格的節點數分別為11 495和53 361,網格單元數分別為52×54和110×120,其中粗網格以及對應的Bump區域局部放大網格如圖4所示,細網格僅對該結構按比例加密.

圖4 Bump湍流流動網格Fig.4 Grid of turbulent flow over the Bump

圖5顯示2種四邊形網格下,rDGP1P2和DGP1方法計算的Bump湍流流動渦黏度等值線.由圖可知,隨著網格密度的增加,rDGP1P2方法和DGP1方法所得的渦黏度計算結果均更加穩定和光滑.圖6給出Bump表面摩擦阻力系數曲線.由圖可見,Bump頂端的摩擦阻力系數有略微的數值振蕩現象,且隨著網格加密振蕩消失.圖7給出rDGP1P2方法和DGP1方法求解Bump算例得到的密度殘差收斂曲線.當殘差下降6個數量級時,rDGP1P2方法需要迭代249時間步,DGP1方法需要迭代187時間步,兩種方法的迭代時間步數略有差別.對比計算效率與精度,rDGP1P2方法單步的計算時間與DGP1方法相當,但是rDGP1P2方法在粗網格上的收斂解與DGP1在細網格上的收斂解基本一致,這說明rDGP1P2方法的計算精度高于DGP1方法.

圖5 Bump湍流流動渦黏度云圖Fig.5 Contours of vortex viscosity of turbulent flow over the Bump

圖6 Bump湍流流動表面摩擦阻力系數Fig.6 Skin friction coefficient of turbulent flow over the Bump

圖7 Bump湍流流動計算密度殘差曲線(粗網格)Fig.7 Convergence histories for density of turbulent flow over the Bump(coarse mesh)

3.3 非定常圓柱繞流層流流動

非定常圓柱繞流層流流動的自由來流馬赫數為Ma∞=0.2,雷諾數為Re=200.特征長度選取圓柱的直徑,圓柱壁面采用無滑移壁面邊界條件,遠場采用特征邊界條件.計算網格采用二階全單元網格,包含8 406個三角形單元和4 568個四邊形單元,如圖8所示.計算采用四階顯式單對角隱式Runge-Kutta(explicit singly-diagonal implicit Runge-Kutta,ESDIRK)[14]進行時間推進,Δt=0.5,共計算2 000個時間步長.

圖8 二維非定常圓柱繞流層流計算網格Fig.8 Grid of 2-dimensional transient laminar cylindrical winding flow

圖9給出了rDGP1P2重構方法計算t=1 000時刻的馬赫數云圖和熵云圖.由圖可見,對于非定常圓柱繞流層流,用基于BR2格式黏性通量離散的rDGP1P2方法能夠準確地捕捉層流的渦街結構.圖10給出rDGP1P2方法和DGP1方法計算所得的升力系數和阻力系數,橫坐標為無量綱時間.由圖可見,隨著時間的推進,非定常渦結構體現為某一定頻率的升阻力變化周期現象,與實際情況相符;同時,rDGP1P2方法和DGP1方法求解的升阻力系數曲線幅度和頻率吻合程度較高,驗證了混合重構BR2間斷伽遼金方法在非定常數值計算中的算法穩定性和程序可靠性.

圖9 二維非定常圓柱繞流層流計算馬赫數和熵分布云圖Fig.9 Contours of Mach number and entropy of 2-dimensional transient laminar cylindrical winding flow

圖10 二維非定常圓柱繞流層流計算升力系數和阻力系數Fig.10 Lift and drag coefficient of 2-dimensional transient laminar cylindrical winding flow

4 結論

針對三階精度傳統間斷伽遼金方法求解NS方程和RANS方程自由度多、運算量大等特點,本工作綜合利用點值重構和積分least-squares重構方法的優勢,并結合BR2黏性通量求解格式,發展了基于BR2格式的三階混合重構間斷伽遼金方法rDGP1P2.通過層流和湍流典型算例驗證了發展的rDGP1P2方法的計算精度、計算效率和收斂性,所得結論如下.

(1)rDGP1P2方法在P1自由度的基礎上通過運算量相對少的混合重構使計算精度由二階提高到三階,從而大幅度減少運算量,提高了計算效率.

(2)發展的rDGP1P2方法具有良好的收斂性,能夠滿足層流、湍流和非定常流動數值模擬的要求.

(3)發展的混合重構間斷伽遼金方法不僅能夠用于層流流場的計算,而且能夠用于湍流流場的數值模擬,數值計算結果和實驗結果吻合較好.