r 型網(wǎng)格自適應(yīng)在間斷Galerkin 有限元激波捕捉中的應(yīng)用

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

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

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

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

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

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

本文針對DG 方法在數(shù)值模擬強間斷方面存在的計算收斂性差、魯棒性低的問題,建立了一種基于網(wǎng)格點壓力和網(wǎng)格點相對位移變化為控制參數(shù)的r 型網(wǎng)格自適應(yīng)方法及Venkatakrishnan 限制器來提高DG 方法模擬強間斷流場的能力。 重點介紹了r 型網(wǎng)格自適應(yīng)獨特的網(wǎng)格點移動驅(qū)動策略;Venkatakrishnan 限制器考慮了DG 方法高斯數(shù)值積分的特殊性。 首先,給出了DG 方法,包括控制方程、空間離散、時間離散和Venkatakrishnan限制器;其次,介紹了r 型網(wǎng)格自適應(yīng)方法,詳細(xì)給出了網(wǎng)格移動策略;再次,采用NACA0012 翼型跨聲速繞流驗證DG 方法及限制器;然后,數(shù)值模擬了并列NACA0012 翼型超聲速流動,驗證了二維流動下本文r 型網(wǎng)格自適應(yīng)和限制器技術(shù);最后,模擬了三維雙半球-圓柱流動干擾流動,驗證了基于r 型網(wǎng)格自適應(yī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 限 制 器 由于其較好的殘值收斂性、魯棒性和計算精度而得到了廣泛應(yīng)用。 本文參考有限體積法中這類限制器的構(gòu)造原理,得到適合DG 方法的斜率限制器。

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

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

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

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

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

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

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

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

本文r 型網(wǎng)格自適應(yī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)格自適應(yī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)負(fù)體積,則修改移動量Δr0,Δ?r0=C·Δr0,C為限制系數(shù),本文取0.5,并再次檢查網(wǎng)格質(zhì)量,如果出現(xiàn)負(fù)體積,再按Δ?r0=C·Δr0減小移動量,直到體積正常。

3 數(shù)值算例

本文所有算例都在中國空氣動力研究與發(fā)展中心的計算機群上完成。 首先,采用NACA0012翼型驗證本文DG 方法和Venkatakrishnan 限制器;然后,采用典型超聲速二維及三維流動算例展示r 型網(wǎng)格自適應(yīng)方法在DG 方法的激波捕捉技術(shù)中的應(yīng)用能力。 算例表明,本文限制器方法具有較好的魯棒性,能夠處理各向異性網(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,自適應(yīng)迭代11 次,在自適應(yīng)迭代中保持翼型表面網(wǎng)格點不移動。

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

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

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

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

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

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

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

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

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

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

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

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

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

圖18 雙半球-圓柱第4 次自適應(yīng)計算網(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 次自適應(yīng)后的網(wǎng)格及流場、第4 次自適應(yīng)后的網(wǎng)格及流場。 從圖21 看到,r 型網(wǎng)格自適應(yī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 給出了自適應(yīng)前后對稱面上z=0.04 m截面壓力分布。 自適應(yīng)后激波寬度減小,激波前后壓力峰值更大。 圖23 和圖24 分別給出了自適應(yīng)前后下部圓柱和上部圓柱y=0 截面表面壓力分布。 自適應(yīng)后影響上下圓柱壓力分布,由于本文的r 型網(wǎng)格自適應(yīng)未移動表面網(wǎng)格,影響量較小。 圖25 給出了x=0.22 m 截面處下部圓柱的表面壓力分布。 上部圓柱的脫體激波在此處打到下部圓柱,看到自適應(yīng)前后該處的壓力有明顯變化,自適應(yīng)后的表面壓力更低。

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

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

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

圖25 自適應(yīng)前后下部圓柱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)格自適應(yīng)方法,通過二維及三維算例分析不同自適應(yīng)迭代步后的表面和空間壓力分布,展示了本文方法在激波精細(xì)捕捉方面的能力,并能得到如下結(jié)論:

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

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

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