龔小權(quán) 吳曉軍 唐靜 李明 張健

(中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 綿陽 621000)

間斷Galerkin(discontinuous Galerkin,DG)有限元[1-2]方法作為一種適用于非結(jié)構(gòu)網(wǎng)格的高精度方法,具有計算精度高,計算模板緊致,適合處理復雜外形及邊界、初值問題,并行能力強,自適應計算方便等優(yōu)點。

與其他高精度格式類似,DG 方法在模擬包含強間斷流場時存在收斂性和魯棒性差的問題。目前,大部分學者通過添加人工黏性[3-4]或限制器[5-6]來抑制激波附近的震蕩。 人工黏性方法作為一種先驗方法,在控制方程中添加代表黏性的二階導數(shù)項,該項數(shù)值與計算網(wǎng)格、模擬問題相關(guān),具有一定經(jīng)驗性,其對數(shù)值解的干擾沒有統(tǒng)一的評估標準。 限制器作為一種后驗方法,通過限制或重構(gòu)多項式高階自由度抑制數(shù)值震蕩,其存在殘差收斂滯止、結(jié)果收斂性差的問題。 目前,DG 方法的激波捕捉技術(shù)已成為其發(fā)展的瓶頸之一,國內(nèi)外在限制器方面開展了大量研究,包括Barth-Jespersen 限制器[7]、Venkatakrishnan[8]限制器,以及Shu 等[9]按照WENO 重構(gòu)思想在三角形網(wǎng)格上提出的HWENO 限制器。

激波前后密度、壓力變化大,在數(shù)值模擬前由于無法預測激波位置,初始網(wǎng)格一般為均勻分布網(wǎng)格,激波附近網(wǎng)格分辨率不足進一步加劇DG方法在模擬強間斷問題時的收斂性和魯棒性問題。 基于以上原因,基于網(wǎng)格自適應和限制器技術(shù)開展DG 方法的間斷流場數(shù)值模擬是重要研究方向,能夠提高計算精度,改善激波模擬分辨率。

耦合流動特征、動態(tài)調(diào)整網(wǎng)格單元分布的r型、h 型網(wǎng)格自適應技術(shù)能夠改善氣動特性計算精度,提高復雜流動模擬精度[10-11],其在誤差修正[12]、網(wǎng)格無關(guān)性[13]、驗證與確認[14]、非定常流動模擬[15]等方面都開展了應用研究。 Karen等[16]在數(shù)值模擬飛行器再入過程中采用基于網(wǎng)格點溫度的r 型網(wǎng)格點移動方法,引入網(wǎng)格點周圍單元的最小內(nèi)切圓半徑及基于雷諾數(shù)的網(wǎng)格點目標長度作為控制參數(shù)。 從結(jié)果看到,該方法實現(xiàn)了r 型網(wǎng)格自適應,在強激波附近聚集2 ～3 層網(wǎng)格。 吳澤艷[17]、孫強[18]等將h 型網(wǎng)格自適應方法應用在二維DG 的激波捕捉方法中,采用HWENO 限制器捕捉激波,有效實現(xiàn)了激波區(qū)域的網(wǎng)格自適應加密。 王利和周偉江[19]基于伴隨方法開展了網(wǎng)格自適應技術(shù)的二維DG 方法研究,提高了計算精度。 Woopen[20]、Fidkowski[21-22]等在DG 方法上建立了基于伴隨方法發(fā)展網(wǎng)格自適應技術(shù),并將其用于阻力預測,取得了較好的效果。 r 型網(wǎng)格自適應方法[23-24]作為一種激波捕捉和誤差估計方法,通過移動網(wǎng)格點使網(wǎng)格聚集在流動變化劇烈區(qū)域,實現(xiàn)對流動特征區(qū)域網(wǎng)格加密,改善網(wǎng)格分辨率,提高流場模擬精度和氣動特性求解精度。 r 型網(wǎng)格自適應的最大優(yōu)點是網(wǎng)格單元的連接關(guān)系和網(wǎng)格的空間拓撲結(jié)構(gòu)在自適應迭代過程中始終保持不變,流場解算器在運行過程中不需要調(diào)整網(wǎng)格數(shù)據(jù)結(jié)構(gòu),不必考慮并行動態(tài)負載平衡問題。

本文針對DG 方法在數(shù)值模擬強間斷方面存在的計算收斂性差、魯棒性低的問題,建立了一種基于網(wǎng)格點壓力和網(wǎng)格點相對位移變化為控制參數(shù)的r 型網(wǎng)格自適應方法及Venkatakrishnan 限制器來提高DG 方法模擬強間斷流場的能力。 重點介紹了r 型網(wǎng)格自適應獨特的網(wǎng)格點移動驅(qū)動策略;Venkatakrishnan 限制器考慮了DG 方法高斯數(shù)值積分的特殊性。 首先,給出了DG 方法,包括控制方程、空間離散、時間離散和Venkatakrishnan限制器;其次,介紹了r 型網(wǎng)格自適應方法,詳細給出了網(wǎng)格移動策略;再次,采用NACA0012 翼型跨聲速繞流驗證DG 方法及限制器;然后,數(shù)值模擬了并列NACA0012 翼型超聲速流動,驗證了二維流動下本文r 型網(wǎng)格自適應和限制器技術(shù);最后,模擬了三維雙半球-圓柱流動干擾流動,驗證了基于r 型網(wǎng)格自適應的DG 方法在激波捕捉方面的能力。

1 間斷Galerkin 有限元方法

1.1 控制方程

1.2 DG 方法空間離散

空間離散前,整個求解區(qū)域Ω被劃分為互不重疊的非結(jié)構(gòu)單元Ωk,對于三維問題,Ωk包括三棱柱、四面體、六面體和金字塔。

將方程(1)兩端乘以測試函數(shù)?j(x,y,z),再在單元內(nèi)積分,并使用分部積分方法得到DG 的弱形式。

式中:n為單元邊界面的外法向;?Ω為單元Ω的邊界;Γ為單元Ω的邊界?Ω的微分符號。

假設(shè)變量在單元內(nèi)具有以下多項式分布:

1.3 DG 方法時間離散

對方程(4)采用一階歐拉后差:

本文采用LU-SGS 隱式離散[25]求解該線性系統(tǒng),求解過程中引入局部時間步長、變CFL 數(shù)等加速收斂技術(shù)。

1.4 Venkatakrishnan 限制器

為保證計算魯棒性,在間斷區(qū)域及網(wǎng)格質(zhì)量較差區(qū)域,采用限制器來抑制震蕩DG 方法的限制器可能帶來計算精度降低、殘值收斂滯止問題。DG 方法作為一種內(nèi)自由度方法,具有模板緊致特性,這要求限制過程也能保持緊致特性。 有限體積的斜率類限制器通過改進后能用于DG 方法,其核心思想是:限制單元變量梯度,使單元任何一點變量值不超過相鄰單元的最值,防止出現(xiàn)“過沖”現(xiàn)象。 Barth 和Jespersen[7]建立了非結(jié)構(gòu)網(wǎng)格中常用的一類Barth-Jespersen 限制器,該限制器在計算限制因子過程中引入不可微函數(shù)帶來殘值收斂性問題。 Venkatakrishnan[8]通過近似min(1,y) 函數(shù)避免不可微過程,建立了Venkatakrishnan 限 制 器。 Venkatakrishnan 限 制 器 由于其較好的殘值收斂性、魯棒性和計算精度而得到了廣泛應用。 本文參考有限體積法中這類限制器的構(gòu)造原理,得到適合DG 方法的斜率限制器。

判斷當前單元所有邊界面的高斯積分點變量值是否滿足以下條件:

積分點限制因子αj具體形式如下:

其中:K為量級為1 的自由可調(diào)參數(shù),其大小影響殘值收斂性,值越小代表限制越大;h為單元的參考長度。

韓國的風水說主要應用于墓地、部落等的選址以及都城、寺剎、住宅等建造方面。 就墓地而言，中國道教風水學說傳入朝鮮半島后，古代朝鮮人民也認為如果把祖先葬在藏風納氣的風水寶地上，這種“佳吉氣”就會感應到子孫后代身上，讓子孫繁衍昌盛； 如果祖先所葬之地沒有“氣”，則會出現(xiàn)天壤之別的情形。

2 r 型網(wǎng)格自適應方法

r 型網(wǎng)格自適應方法通過移動網(wǎng)格點使網(wǎng)格聚集在流動變化劇烈的區(qū)域,實現(xiàn)流動特征區(qū)域網(wǎng)格加密,改善網(wǎng)格分辨率,進而提高流場模擬精度。 網(wǎng)格自適應技術(shù)主要包括誤差估計、網(wǎng)格加密/稀疏、表面網(wǎng)格投影、空間體網(wǎng)格匹配及并行實現(xiàn)。

本文r 型網(wǎng)格自適應方法在自適應過程中限制表面網(wǎng)格點移動,因此不考慮表面網(wǎng)格投影。r 型網(wǎng)格自適應在自適應過程中網(wǎng)格單元的連接關(guān)系和網(wǎng)格空間拓撲結(jié)構(gòu)保持不變,也不考慮并行問題。 本文重點介紹r 型網(wǎng)格自適應網(wǎng)格點移動策略。

Karen 等[16]采用溫度作為網(wǎng)格點移動控制變量適合高超聲速流動,對于一般的跨聲速或超聲速流動基于壓力變化可能更具有普適性。 從文獻[16]算例看到,基于單一的溫度作為控制參數(shù),網(wǎng)格自適應區(qū)域的網(wǎng)格聚集數(shù)較少,經(jīng)過100 次自適應迭代后,只有2 ～3 排網(wǎng)格點移動,激波附近網(wǎng)格點移動信息向其周圍網(wǎng)格點傳遞不明顯。

本文r 型網(wǎng)格自適應網(wǎng)格點移動的驅(qū)動力由網(wǎng)格點的流場變量值確定,以歸一化網(wǎng)格點的變量差及網(wǎng)格點相對位移量這2 類參數(shù)作為網(wǎng)格點移動的權(quán)函數(shù)。 針對激波問題,本文以壓力P作為網(wǎng)格點移動控制變量。 網(wǎng)格點移動的計算流場如圖1 所示。

圖1 r 型網(wǎng)格自適應網(wǎng)格點移動位移計算流程Fig.1 Calculation process of grid point displacement in r-grid adaptive

以圖2 為例具體說明網(wǎng)格點移動控制過程。

圖2 網(wǎng)格點分布示意圖Fig.2 Distribution diagram of grid points

1) 采用距離權(quán)計算網(wǎng)格點0 及周圍網(wǎng)格點1、2、3、4、5 的壓力。

7) 限制特殊點的位移,如對稱面上網(wǎng)格點y方向位移為零(對稱面在y=0 平面)。

8) 移動網(wǎng)格后檢查網(wǎng)格質(zhì)量,如果出現(xiàn)負體積,則修改移動量Δr0,Δ?r0=C·Δr0,C為限制系數(shù),本文取0.5,并再次檢查網(wǎng)格質(zhì)量,如果出現(xiàn)負體積,再按Δ?r0=C·Δr0減小移動量,直到體積正常。

3 數(shù)值算例

本文所有算例都在中國空氣動力研究與發(fā)展中心的計算機群上完成。 首先,采用NACA0012翼型驗證本文DG 方法和Venkatakrishnan 限制器;然后,采用典型超聲速二維及三維流動算例展示r 型網(wǎng)格自適應方法在DG 方法的激波捕捉技術(shù)中的應用能力。 算例表明,本文限制器方法具有較好的魯棒性,能夠處理各向異性網(wǎng)格及網(wǎng)格尺度急劇變化的網(wǎng)格。

3.1 DG 方法及其限制器驗證

圖3 NACA0012 翼型計算網(wǎng)格Fig.3 Computational grid of NACA0012 airfoil

圖4 不同方法壓力分布比較Fig.4 Comparison of pressure distribution between different methods

3.2 并列兩段NACA0012 翼型超聲速繞流

本節(jié)數(shù)值模擬了2 個并列放置的NACA0012翼型的超聲速繞流。 計算的初始網(wǎng)格如圖5 所示,遠場為30 倍弦長,共25 982 個三棱柱單元,展向拉伸一排網(wǎng)格。 其中一個翼型相對另外一個翼型上移和后移0. 5 m,計算來流馬赫數(shù)Ma=2.0,來流迎角0°,來流靜壓P=26 500 Pa,來流靜溫T=223.25 K,自適應迭代11 次,在自適應迭代中保持翼型表面網(wǎng)格點不移動。

翼型錯位放置,翼型間產(chǎn)生復雜的激波干擾,上部翼型的脫體激波在下部翼型的上表面產(chǎn)生反射激波后反射到上部翼型的下表面。 圖5 ～圖9分別給出了初始網(wǎng)格、第1 次自適應、第3 次自適應、第5 次自適應、第11 次自適應后的網(wǎng)格分布及壓力系數(shù)空間分布。 可以看到,初始網(wǎng)格和自適應后的網(wǎng)格都捕捉到脫體激波及翼型間的反射激波;r 型網(wǎng)格自適應方法準確判別到流場劇烈變化區(qū)域,并將網(wǎng)格沿著激波方向各向異性自適應聚合,自適應后的激波捕捉更清晰和銳利;隨著自適應迭代次數(shù)的增加,激波附近的網(wǎng)格逐步被自適應加密,網(wǎng)格更緊密聚合于激波處,約有4 ～5排網(wǎng)格,捕捉激波的網(wǎng)格寬度更小,拉伸比更大。自適應方法同樣也準確識別翼型之間的反射激波,并沿激波方向進行了自適應加密。

圖7 NACA0012 翼型第3 次自適應后計算網(wǎng)格及流場Fig.7 Computational grid and flowfield of NACA0012 airfoil after third r-adaptation

圖8 NACA0012 翼型第5 次自適應后計算網(wǎng)格及流場Fig.8 Computational grid and flowfield of NACA0012 airfoil after fifth r-adaptation

圖9 NACA0012 翼型第11 次自適應后計算網(wǎng)格及流場Fig.9 Computational grid and flowfield of NACA0012 airfoil after eleventh r-adaptation

圖10 給出了第11 次網(wǎng)格自適應后翼型脫體激波附近網(wǎng)格放大圖。 自適應后,激波附近網(wǎng)格拉伸比變大,網(wǎng)格尺寸劇烈變化,網(wǎng)格質(zhì)量變差,本文的Venkatakrishnan 限制器能夠很好處理該類網(wǎng)格,保證DG 方法計算的魯棒性。 圖11 給出了空間截面位置,圖12 和圖13 分別給出了空間截面z=0. 25 m 和z=0. 068 m 處的壓力分布?？梢钥吹?自適應前后DG 方法都能夠捕捉激波,隨著r 型網(wǎng)格自適應迭代,網(wǎng)格逐步向激波位置聚集,DG 方法捕捉激波更加清晰,主要表現(xiàn)在跨過激波前后壓力峰值更大,激波寬度更窄,激波更銳利。 從圖12 看到,在自適應迭代11 次后,數(shù)值模擬的激波變化已經(jīng)很小。

圖10 自適應網(wǎng)格局部放大圖Fig.10 Local enlarged view of adaptation grid

圖11 空間截面位置Fig.11 Space section position

圖12 z =0.25 m 處壓力分布及激波處局部放大圖Fig.12 Pressure distribution at section z =0.25 m and local enlarged view of shock wave

圖13 z =0.068 m 處壓力分布及激波處局部放大圖Fig.13 Pressure distribution at section z =0.068 m and local enlarged view of shock wave

圖14 展示了自適應計算過程中翼型阻力系數(shù)CD隨計算步數(shù)變化曲線。 阻力系數(shù)在每次網(wǎng)格自適應后只需較小計算步就收斂,隨著自適應迭代增加,阻力系數(shù)增加量逐漸減小。

圖14 自適應過程阻力系數(shù)隨計算步數(shù)變化曲線Fig.14 Curve of drag coefficient varying with iteration steps in adaptive process

圖15 給出了阻力系數(shù)和Venkatakrishnan 限制器系數(shù)平均值隨自適應迭代步變化曲線。 阻力系數(shù)隨網(wǎng)格自適應迭代逐步增大,斜率逐步減小,限制系數(shù)平均值及其斜率隨網(wǎng)格自適應迭代步數(shù)都逐步減小。 網(wǎng)格自適應步數(shù)越大,網(wǎng)格拉伸比越大,網(wǎng)格分布呈現(xiàn)更加不均勻性,因此限制系數(shù)會逐漸減小。 但由于網(wǎng)格主要在激波附近發(fā)生移動,而數(shù)值模擬激波時計算自身就要降階,此處的網(wǎng)格移動對計算精度影響很小,算例驗證了本文r 網(wǎng)格自適應方法、DG 方法及Venkatakrishnan 限制器的魯棒性。

圖15 阻力系數(shù)及限制器系數(shù)平均值隨自適應步數(shù)變化Fig.15 Resistance and average value of limitation coefficients with adaptive steps

3.3 三維雙半球-圓柱流動干擾模擬

為進一步驗證本文DG 方法的限制器及網(wǎng)格自適應技術(shù),本節(jié)計算了三維雙半球-圓柱流動干擾模擬算例。 計算網(wǎng)格如圖16 所示,共527 612 個四面體單元。 計算來流馬赫數(shù)Ma=2.0,來流迎角0°,來流靜壓P= 26 500 Pa,來流靜溫T=223.25 K,自適應迭代4 次。 圖16 ～圖18 給出了雙半球-圓柱對稱面及表面初始網(wǎng)格和流場、第2次自適應后的網(wǎng)格及流場、第4 次自適應后的網(wǎng)格及流場。 隨著自適應迭代步的增加,網(wǎng)格逐步向激波附近聚集,用于捕捉激波的網(wǎng)格的寬度減小,拉伸比增大,激波分辨率增高。

圖16 三維雙半球-圓柱初始計算網(wǎng)格及流場Fig.16 Initial computational grid and flowfield of 3D double hemispherical cylinder

圖17 雙半球-圓柱第2 次自適應計算網(wǎng)格及流場Fig.17 Computational grid and flowfield of double hemispherical cylinder after second r-adaptation

圖18 雙半球-圓柱第4 次自適應計算網(wǎng)格及流場Fig.18 Computational grid and flowfield of double hemispherical cylinder after fourth r-adaptation

圖19 ～圖21 給出了x=0.22 m 處的初始網(wǎng)格及流場、第2 次自適應后的網(wǎng)格及流場、第4 次自適應后的網(wǎng)格及流場。 從圖21 看到,r 型網(wǎng)格自適應后空間網(wǎng)格向激波附近聚集,來自上部圓柱的脫體激波打到下部圓柱,其激波分辨率提高。

圖19 雙半球-圓柱初始x =0.22 m 處空間網(wǎng)格及流場Fig.19 Spatial grid and flowfield at x =0.22 m section of double hemispherial cylinder

圖22 給出了自適應前后對稱面上z=0.04 m截面壓力分布。 自適應后激波寬度減小,激波前后壓力峰值更大。 圖23 和圖24 分別給出了自適應前后下部圓柱和上部圓柱y=0 截面表面壓力分布。 自適應后影響上下圓柱壓力分布,由于本文的r 型網(wǎng)格自適應未移動表面網(wǎng)格,影響量較小。 圖25 給出了x=0.22 m 截面處下部圓柱的表面壓力分布。 上部圓柱的脫體激波在此處打到下部圓柱,看到自適應前后該處的壓力有明顯變化,自適應后的表面壓力更低。

圖22 自適應前后對稱面z =0.04 m 截面壓力分布Fig.22 Pressure distribution at z =0.04 m section of symmetry plane before and after adaptation

圖23 自適應前后下部圓柱y =0 截面表面壓力分布Fig.23 Surface pressure distribution of lower cylinder at y =0 section before and after adaptation

圖24 自適應前后上部圓柱y =0 截面表面壓力分布Fig.24 Surface pressure distribution of upper cylinder at y =0 section before and after adaptation

圖25 自適應前后下部圓柱x =0.22 m 截面表面壓力分布Fig.25 Surface pressure distribution of lower cylinder at x =0.22 m section before and after adaptotion

4 結(jié) 論

本文針對DG 方法在模擬強間斷流場面臨的問題,發(fā)展了Venkatakrishnan 限制器和r 型網(wǎng)格自適應方法,通過二維及三維算例分析不同自適應迭代步后的表面和空間壓力分布,展示了本文方法在激波精細捕捉方面的能力,并能得到如下結(jié)論:

1) r 型網(wǎng)格自適應方法具備根據(jù)流場特征開展網(wǎng)格自適應能力,能夠沿著激波方向?qū)崿F(xiàn)網(wǎng)格的各向異性自適應移動加密,具有較好的魯棒性,保證網(wǎng)格不出現(xiàn)負體積。

2) DG 方法及Venkatakrishnan 限制器針對超聲速問題具有較好的激波捕捉能力,能夠魯棒處理各向異性非結(jié)構(gòu)網(wǎng)格、針對自適應后的相鄰網(wǎng)格尺寸差異很大的網(wǎng)格具有較高的計算魯棒性。

3) 通過在DG 方法上發(fā)展r 型網(wǎng)格自適應方法能夠明顯改善數(shù)值模擬強激波能力,提高激波分辨率,進一步提高計算精度。