周 帥 肖周芳 付 琳 汪丁順

* (中國航空發(fā)動(dòng)機(jī)研究院,北京 101300)

? (杭州電子科技大學(xué)計(jì)算機(jī)學(xué)院,杭州 310018)

引言

隨著計(jì)算機(jī)軟硬件技術(shù)的發(fā)展,以計(jì)算流體力學(xué)(computational fluid dynamics,CFD)為核心的數(shù)值模擬技術(shù)已成為空氣動(dòng)力學(xué)領(lǐng)域工程與科學(xué)問題分析中一項(xiàng)不可或缺的工具[1].然而,針對(duì)復(fù)雜幾何外形和復(fù)雜流動(dòng)問題,CFD 計(jì)算仍然面臨求解精度和效率的挑戰(zhàn)[2].網(wǎng)格自適應(yīng)技術(shù)[2-8]和高階精度數(shù)值方法[9-14]的結(jié)合有望提升CFD 復(fù)雜問題適應(yīng)能力,還可在更少的網(wǎng)格上獲得更快的漸進(jìn)收斂速度[2,10].然而,這兩項(xiàng)技術(shù)卻并未在主流商業(yè)CFD 軟件中使用.這主要是由于網(wǎng)格自適應(yīng)技術(shù)涉及諸多環(huán)節(jié)的協(xié)同工作,其還存在魯棒性問題,并且要將這兩項(xiàng)技術(shù)結(jié)合還需要解決一系列技術(shù)難題[2].其中一個(gè)難題便是如何將前一自適應(yīng)迭代步中的計(jì)算結(jié)果高精度地插值到后一迭代步的自適應(yīng)更新網(wǎng)格中,這一個(gè)過程稱為物理量插值.

物理量插值又被稱為物理量轉(zhuǎn)移或物理量重映[15-25],是將物理量從一網(wǎng)格上轉(zhuǎn)移到另一網(wǎng)格上,其常被應(yīng)用于網(wǎng)格變形[16-19]和網(wǎng)格自適應(yīng)更新[15]等計(jì)算問題中.在自適應(yīng)流場(chǎng)計(jì)算中,通過物理量插值可使當(dāng)前流場(chǎng)計(jì)算延續(xù)上一迭代步中的計(jì)算狀態(tài),從而加速流場(chǎng)計(jì)算的收斂速度[20-21].該過程中需關(guān)注物理量守恒和插值精度等特性.前者將物理量在從背景網(wǎng)格單元轉(zhuǎn)移到新網(wǎng)格單元上后,需保持物理量的守恒;對(duì)于后者,物理量的插值精度應(yīng)和當(dāng)前采用的數(shù)值方法精度一致.

以上物理量插值的兩個(gè)特性對(duì)于提升整個(gè)自適應(yīng)流場(chǎng)計(jì)算流程的穩(wěn)定性具有重要的作用[21],尤其在非定常流動(dòng)計(jì)算中,由物理量插值導(dǎo)致的數(shù)值誤差會(huì)累積到之后的流動(dòng)計(jì)算中.簡(jiǎn)單的線性插值算法無法滿足物理量守恒特性,且其只能達(dá)到一階插值精度.文獻(xiàn)[15]提出了針對(duì)二維非結(jié)構(gòu)網(wǎng)格的滿足物理量守恒及二階精度的物理量插值算法,物理守恒通過局部網(wǎng)格相交運(yùn)算實(shí)現(xiàn),二階精度通過重構(gòu)物理量梯度實(shí)現(xiàn).隨后,Alauzet[22]繼續(xù)將該物理量插值算法拓展到了三維非結(jié)構(gòu)網(wǎng)格上.文獻(xiàn)[17,23]采用基于ENO 方法實(shí)現(xiàn)物理守恒的物理插值方法,文獻(xiàn)[23]在二維網(wǎng)格上實(shí)現(xiàn)了三階精度的物理量插值.針對(duì)二維四邊形網(wǎng)格物理量插值問題,Lei 等[24]發(fā)展了基于WENO 的高精度插值方法.Zhang 等[25]針對(duì)高階間斷伽遼金(DG)數(shù)值解在變形網(wǎng)格和移動(dòng)網(wǎng)格應(yīng)用中的物理量插值需求,設(shè)計(jì)了DG 插值方法,實(shí)現(xiàn)了守恒和高階的物理量插值.

本文提出一類滿足物理量守恒和高精度的物理量插值方法,以適應(yīng)高階精度自適應(yīng)流動(dòng)計(jì)算.為實(shí)現(xiàn)流場(chǎng)插值過程中物理量守恒,該方法先計(jì)算新舊網(wǎng)格的重疊區(qū)域,然后將物理量從重疊區(qū)域的舊網(wǎng)格轉(zhuǎn)移到新網(wǎng)格中.為滿足高階精度要求,先采用kexact 最小二乘方法[26]對(duì)舊網(wǎng)格上的數(shù)值解進(jìn)行重構(gòu)，得到滿足精度要求的物理量分布多項(xiàng)式,隨后采用高斯數(shù)值積分實(shí)現(xiàn)物理量精確地轉(zhuǎn)移到新網(wǎng)格的各單元上.本文算法有助于推進(jìn)網(wǎng)格自適應(yīng)技術(shù)和高階精度數(shù)值方法的結(jié)合,從而提升CFD 復(fù)雜問題適應(yīng)能力.

1 算法流程

考慮二維自適應(yīng)流動(dòng)計(jì)算,記Mnew為前迭代步網(wǎng)格(也稱為新網(wǎng)格),Mback為上一迭代步中的網(wǎng)格(也稱為背景網(wǎng)格),本文算法目的為將分布在Mback上的物理量插值到Mnew中,使得插值后的解具有(r+1)階精度并滿足在插值過程中物理量守恒.為實(shí)現(xiàn)上述目的,本文算法步驟如下.

(1) 物理量重構(gòu):針對(duì)Mback中的每個(gè)網(wǎng)格單元,構(gòu)建描述流場(chǎng)物理量在該單元中分布的r次多項(xiàng)式表達(dá)式U(x,y).

(2) 網(wǎng)格相交計(jì)算:利用Mback中的網(wǎng)格邊將Mnew中的每個(gè)網(wǎng)格單元分割得到多個(gè)子單元,這些子單元為Mnew中對(duì)應(yīng)單元與Mback中不同網(wǎng)格單元的重疊部分.

(3) 物理量轉(zhuǎn)移:針對(duì)Mnew中的每個(gè)網(wǎng)格單元,根據(jù)分布于Mback上的流場(chǎng)物理量U(x,y)累積分布于該單元對(duì)應(yīng)子單元上的物理量,累積結(jié)果作為該單元上的物理量總量,并以該單元單位面積上的物理量為其參考點(diǎn)上的物理量值.

上述步驟中,步驟1 和步驟3 是實(shí)現(xiàn)高階精度流場(chǎng)物理量插值的關(guān)鍵.其中,步驟1 用于將背景網(wǎng)格上的物理量重構(gòu)為r階精度,步驟3 用于確保物理量轉(zhuǎn)移到新網(wǎng)格后保持r階精度,將在第2 節(jié)和第3 節(jié)分別討論這兩部分內(nèi)容.步驟2 中分割后的子單元用于實(shí)現(xiàn)插值過程中物理量守恒,其中涉及網(wǎng)格求交計(jì)算較為成熟[27-29],本文不再贅述.特別地,本文對(duì)在二維非結(jié)構(gòu)三角形網(wǎng)格間插值進(jìn)行闡述,其思想同樣適用于三維非結(jié)構(gòu)四面體網(wǎng)格.

2 物理量重構(gòu)

本文采用應(yīng)用于有限體積法中的k-exact 最小二乘法[26]進(jìn)行舊網(wǎng)格上的物理量重構(gòu).為了完整性,將結(jié)合本文中的流場(chǎng)插值需求,簡(jiǎn)述該算法過程.物理量重構(gòu)的目的為針對(duì)背景網(wǎng)格Mback上的任一網(wǎng)格單元i,構(gòu)建關(guān)于該單元參考點(diǎn)的物理量分布r次多項(xiàng)式

該式的計(jì)算僅跟網(wǎng)格幾何位置相關(guān),關(guān)于其計(jì)算方法請(qǐng)參考文獻(xiàn)[26].

從式(1)和式(2)可知,在二維情形中,次數(shù)為r的多項(xiàng)式含有(r+1)(r+2)/2 個(gè)未知項(xiàng)系數(shù).為求解這些未知項(xiàng),需選擇(r+1)(r+2)/2 個(gè)單元i的鄰接單元并針對(duì)每個(gè)單元建立一個(gè)如式(2)的等式.在實(shí)踐中,通常選擇多于(r+1)(r+2)/2 個(gè)鄰接單元構(gòu)建方程數(shù)多于未知項(xiàng)數(shù)的方程組[24],這些鄰接單元被稱為單元i的重構(gòu)模板.隨后,采用最小二乘法求解這些未知項(xiàng),即式(1)中的多項(xiàng)式系數(shù).

3 流場(chǎng)插值

3.1 滿足物理量守恒的插值

在自適應(yīng)流動(dòng)計(jì)算中,不同迭代步的網(wǎng)格通常相互獨(dú)立.不失一般性,本文考慮相互獨(dú)立的網(wǎng)格間的物理量插值情形.如圖1(a)所示,黑色邊表示的網(wǎng)格為Mback,紅邊表示的網(wǎng)格為Mnew(為便于展示,此處只以部分網(wǎng)格單元示意這兩套網(wǎng)格).圖1(b)展示網(wǎng)格相交計(jì)算后獲得Mnew中網(wǎng)格單元△ABC對(duì)應(yīng)的子單元結(jié)果,共有4 個(gè)子單元,分別為與Mback中三個(gè)單元重疊的區(qū)域.

圖1 二維物理量守恒插值示意圖,背景網(wǎng)格和當(dāng)前網(wǎng)格分別用黑邊和紅邊表示Fig.1 Schematic diagram of two-dimensional conservation of physical quantity interpolation,the background mesh and the current mes are represented by red and black edges respectively

3.2 高階精度插值

從式(4)可知,為求解Mnew中網(wǎng)格單元內(nèi)的物理量需進(jìn)行體積分計(jì)算,該計(jì)算的精度將影響最終的流場(chǎng)插值精度.為獲得高精度積分結(jié)果,本文采用高斯積分求解式(4)右端的積分項(xiàng),即

式中,Nq為高斯積分點(diǎn)個(gè)數(shù),(xq,yq)為積分點(diǎn)坐標(biāo),wq為對(duì)應(yīng)該積分點(diǎn)的權(quán)重.注意,高斯數(shù)值積分的精度和采用的高斯積分點(diǎn)個(gè)數(shù)相關(guān).為此,針對(duì)不同的精度要求,需采用不同的積分點(diǎn)數(shù)及相應(yīng)的積分點(diǎn)坐標(biāo)和積分點(diǎn)權(quán)重.

為支持高階精度的流場(chǎng)插值,針對(duì)積分區(qū)域?yàn)槿切吻樾?采用文獻(xiàn)[26,30]中推薦的高斯積分點(diǎn)信息.表1 展示了針對(duì)2～ 4 階精度被積函數(shù)所采用的高斯積分點(diǎn)數(shù)、坐標(biāo)和相應(yīng)權(quán)重.注意,表1 中所列的積分點(diǎn)坐標(biāo)為其重心坐標(biāo).當(dāng)所采用的流場(chǎng)數(shù)值計(jì)算方法精度高于4 階時(shí),需采用更多的高斯積分點(diǎn)數(shù).

表1 不同階數(shù)精度對(duì)應(yīng)的積分點(diǎn)信息Table 1 Information of integration points corresponding to different order

4 數(shù)值實(shí)驗(yàn)

選用兩個(gè)例子測(cè)試本文提出的面向高階自適應(yīng)流動(dòng)計(jì)算的流場(chǎng)插值方法,第一個(gè)例子通過具有解析解的流場(chǎng)驗(yàn)證該算法的插值誤差,第二個(gè)例子則將該算法應(yīng)用于高階自適應(yīng)流動(dòng)計(jì)算中.文中所有測(cè)試?yán)泳谛⌒凸ぷ髡旧线\(yùn)行,計(jì)算機(jī)配置為,內(nèi)存容量:16 GB;CPU 型號(hào):六核Inter Core i7-8700,主頻:3.2 GHz.

4.1 圓環(huán)模型

該模型為四分之一圓環(huán)模型[26],圓環(huán)內(nèi)圓半徑為2、外圓半徑為3 (如圖2 和圖3 所示).考慮非黏性超音速渦流流經(jīng)該模型,流場(chǎng)邊界條件為:內(nèi)圓邊界為入口邊界,在此處的馬赫數(shù)為Mi=2、質(zhì)量密度ρi=1.區(qū)域其他位置的流場(chǎng)變量(密度ρ、速度分量u、速度分量v和壓力P)值由以下式子定義

圖2 第一組網(wǎng)格Fig.2 The first group of meshes

圖3 第二組網(wǎng)格Fig.3 The second group of meshes

在自適應(yīng)流場(chǎng)計(jì)算中,為延續(xù)上一迭代步的計(jì)算狀態(tài),需要在每一迭代步開始前執(zhí)行流場(chǎng)插值過程,以獲得該迭代步的初始解.本文模仿自適應(yīng)流場(chǎng)計(jì)算過程中的插值過程,即將定義在一套網(wǎng)格(背景網(wǎng)格)中的流場(chǎng)插值到另一套網(wǎng)格(當(dāng)前網(wǎng)格)中,不同的是此處定義在背景網(wǎng)格上的流場(chǎng)通過流場(chǎng)變量的解析式(6)～ 式(9)求出.為驗(yàn)證本文插值算法的精度,除了從背景網(wǎng)格上獲得的插值解之外,在新網(wǎng)格上也通過流場(chǎng)變量的解析式計(jì)算精確解.隨后,針對(duì)每一個(gè)新網(wǎng)格上的網(wǎng)格單元,根據(jù)以下式子計(jì)算插值誤差

式中,φExact和 φi分別表示在該單元上基于精確值和流場(chǎng)插值得到的單位面積上的物理量值.

為實(shí)現(xiàn)上述插值誤差計(jì)算,設(shè)置兩組不同規(guī)模的網(wǎng)格,每組網(wǎng)格由兩套網(wǎng)格組成,分別為背景網(wǎng)格和當(dāng)前網(wǎng)格,且背景網(wǎng)格的網(wǎng)格單元數(shù)比當(dāng)前網(wǎng)格單元數(shù)少.圖2 和圖3 分別展示了這兩組網(wǎng)格,表2顯示了這兩組不同規(guī)模網(wǎng)格的網(wǎng)格單元數(shù)和網(wǎng)格點(diǎn)數(shù),第二組網(wǎng)格的規(guī)模大于第一組網(wǎng)格的規(guī)模.表3展示了在不同插值精度情形下在這兩組不同規(guī)模的網(wǎng)格上得到的插值誤差的L1范數(shù)、L2范數(shù)和L∞范數(shù).從表中數(shù)據(jù)可知,當(dāng)網(wǎng)格規(guī)模相同時(shí),插值精度階數(shù)越高,插值誤差越小;當(dāng)插值精度相同時(shí),網(wǎng)格規(guī)模越大,插值誤差越小.

表2 兩組不同規(guī)模的網(wǎng)格數(shù)據(jù)Table 2 Two groups of meshes with different scales

表3 在不同網(wǎng)格規(guī)模和不同插值精度下的流場(chǎng)插值誤差Table 3 Interpolation errors of flow field under different mesh scales and different orders of interpolation accuracy

4.2 NACA 0012 模型

以NACA 0012 翼型外流場(chǎng)數(shù)值計(jì)算為例驗(yàn)證本文插值方法在自適應(yīng)流動(dòng)計(jì)算中的有效性.該翼型外流場(chǎng)的計(jì)算條件為:雷諾數(shù)Re=5000,馬赫數(shù)Ma=0.5 和攻角a=0.為獲得該流場(chǎng)數(shù)值解,采用3 階精度的高階流場(chǎng)求解器進(jìn)行自適應(yīng)迭代求解.圖4(a)中的網(wǎng)格為初始計(jì)算網(wǎng)格,由7790 個(gè)網(wǎng)格單元和3966 個(gè)網(wǎng)格點(diǎn)組成.該初始網(wǎng)格在翼型附近區(qū)域分布較稀疏(見網(wǎng)格放大圖),難以捕捉翼型附近區(qū)域的精細(xì)流場(chǎng)特征.為此,采用自適應(yīng)技術(shù)根據(jù)流場(chǎng)解更新計(jì)算域網(wǎng)格,并基于新的網(wǎng)格重新計(jì)算流場(chǎng)數(shù)值解.

本文采用文獻(xiàn)[31-32]中的方法根據(jù)流場(chǎng)解計(jì)算用于更新網(wǎng)格的單元尺寸場(chǎng),并采用局部操作技術(shù)更新網(wǎng)格.最終,通過6 次自適應(yīng)迭代獲得收斂的流場(chǎng)數(shù)值解.圖4(b)～圖4(d)展示了前四次自適應(yīng)迭代計(jì)算的網(wǎng)格,可以發(fā)現(xiàn)隨著自適應(yīng)迭代次數(shù)的增加,翼型附近區(qū)域和尾跡區(qū)的網(wǎng)格越來越密,因此網(wǎng)格單元數(shù)和網(wǎng)格點(diǎn)數(shù)也不斷增加(見表3).圖5 展示了自適應(yīng)計(jì)算收斂過程中氣動(dòng)系數(shù)(升力系數(shù)和阻力系數(shù))隨網(wǎng)格規(guī)模變化的變化.最終,升力系數(shù)值收斂于0.000572(理想值為0);阻力系數(shù)收斂的值為0.0555,這和文獻(xiàn)[33]中展示的基于多種網(wǎng)格計(jì)算得到的阻力系數(shù)值接近.圖6 展示了自適應(yīng)迭代收斂后的馬赫數(shù)分布圖,可以發(fā)現(xiàn)在翼型附近區(qū)域和尾跡區(qū)具有清晰的流場(chǎng)特征.

圖4 翼型計(jì)算域在不同自適應(yīng)計(jì)算迭代步中的網(wǎng)格Fig.4 Meshes of the airfoil computational domain in different adaptive computation iteration steps

圖5 自適應(yīng)計(jì)算收斂過程中氣動(dòng)系數(shù)隨網(wǎng)格規(guī)模變化的變化Fig.5 Convergence of lift and drag coefficients against degrees of freedom (vertices)

圖6 自適應(yīng)迭代收斂后的馬赫數(shù)分布圖Fig.6 The distribution of Mach number after adaptive solution convergences

本文流場(chǎng)插值方法已集成于上述所采用的高階流場(chǎng)求解器中,為說明該插值方法在自適應(yīng)流場(chǎng)計(jì)算中的作用,對(duì)比有流場(chǎng)插值功能和無流場(chǎng)插值功能時(shí),求解器在每一自適應(yīng)迭代步中的求解收斂情況.在每一次自適應(yīng)計(jì)算中,需要在新網(wǎng)格上重新求解流場(chǎng)控制方程,該求解過程最終轉(zhuǎn)化為線性方程組的迭代求解,求解收斂條件為殘差小于給定門限值(本文中采用的值為1.0 × 10?9).在無流場(chǎng)插值步驟時(shí),每次自適應(yīng)計(jì)算從給定初始值開始進(jìn)行流場(chǎng)控制方程求解;而在有流場(chǎng)插值步驟時(shí),自適應(yīng)計(jì)算則以上一次計(jì)算結(jié)果插值到當(dāng)前網(wǎng)格后的解作為初始計(jì)算條件從而延續(xù)上一次自適應(yīng)計(jì)算狀態(tài).

表4 展示了在分別有流場(chǎng)插值和無流場(chǎng)插值步驟時(shí),每一次自適應(yīng)計(jì)算收斂迭代步數(shù)和收斂時(shí)間.因?yàn)榭梢匝永m(xù)上一次計(jì)算狀態(tài),流場(chǎng)插值功能可以有效縮短在新網(wǎng)格上的收斂迭代步數(shù),從而縮短收斂時(shí)間.如在第一次自適應(yīng)計(jì)算時(shí),無流場(chǎng)插值步驟時(shí),求解器需要15 迭代步來求解流場(chǎng)控制方程;而在有流場(chǎng)插值步驟時(shí)求解迭代步縮短到11 步,收斂時(shí)間從40.1 s 降到33 s,加速比為1.22.隨著自適應(yīng)次數(shù)的增加,網(wǎng)格規(guī)模越來越大,所需的數(shù)值求解時(shí)間也增多.從表中數(shù)據(jù)可知,當(dāng)網(wǎng)格規(guī)模越大時(shí),流場(chǎng)插值功能在提高收斂速度方面效果越明顯.如在第6 次自適應(yīng)迭代計(jì)算中,當(dāng)添加流場(chǎng)插值功能后,收斂迭代步數(shù)從29 步降低到10 步,收斂時(shí)間從294.6 s 降低到125.3 s,收斂加速比為2.35.需要說明的是,在本算例中,不管是否有流場(chǎng)插值功能,每一次自適應(yīng)迭代計(jì)算都收斂到相同的結(jié)果.

表4 有無流場(chǎng)插值功能時(shí)求解收斂情況Table 4 Convergence of solution with or without the flow field interpolation

5 結(jié)論

本文提出了一類高精度流場(chǎng)插值方法,實(shí)現(xiàn)將前一迭代步網(wǎng)格中流場(chǎng)數(shù)值解插值到當(dāng)前迭代步網(wǎng)格中,以延續(xù)前一迭代步中的計(jì)算狀態(tài),該插值方法可應(yīng)用于高階精度自適應(yīng)流動(dòng)計(jì)算中.為實(shí)現(xiàn)流場(chǎng)插值過程中物理量守恒,該方法先計(jì)算新舊網(wǎng)格的重疊區(qū)域,然后將重疊區(qū)域的物理量從舊網(wǎng)格轉(zhuǎn)移到新網(wǎng)格中.為滿足高階精度要求,先采用kexact 最小二乘方法對(duì)舊網(wǎng)格上的數(shù)值解進(jìn)行重構(gòu),得到滿足精度要求的物理量分布多項(xiàng)式,隨后采用高斯數(shù)值積分實(shí)現(xiàn)物理量精確地轉(zhuǎn)移到新網(wǎng)格的各單元上.最后,通過兩個(gè)算例驗(yàn)證了本文算法的有效性.

本文雖然以二維三角形網(wǎng)格為例闡述高階插值方法,但其思想同樣適應(yīng)于三維四面體網(wǎng)格情形.不同之處在于,三維情形需要對(duì)四面體網(wǎng)格進(jìn)行求交運(yùn)算以實(shí)現(xiàn)插值過程中物理量守恒,同時(shí)插值過程中需要進(jìn)行體積分計(jì)算.下一步的工作,將考慮將本文方法拓展到三維四面體網(wǎng)格的高階流場(chǎng)插值中.