一類面向高階精度自適應(yīng)流動計算的流場插值方法1)

周 帥 肖周芳 付 琳 汪丁順

* (中國航空發(fā)動機研究院,北京 101300)

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

引言

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

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

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

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

1 算法流程

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

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

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

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

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

2 物理量重構(gòu)

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

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

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

3 流場插值

3.1 滿足物理量守恒的插值

在自適應(yīng)流動計算中,不同迭代步的網(wǎng)格通常相互獨立.不失一般性,本文考慮相互獨立的網(wǎng)格間的物理量插值情形.如圖1(a)所示,黑色邊表示的網(wǎng)格為Mback,紅邊表示的網(wǎng)格為Mnew(為便于展示,此處只以部分網(wǎng)格單元示意這兩套網(wǎng)格).圖1(b)展示網(wǎng)格相交計算后獲得Mnew中網(wǎng)格單元△ABC對應(yīng)的子單元結(jié)果,共有4 個子單元,分別為與Mback中三個單元重疊的區(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)的物理量需進行體積分計算,該計算的精度將影響最終的流場插值精度.為獲得高精度積分結(jié)果,本文采用高斯積分求解式(4)右端的積分項,即

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

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

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

4 數(shù)值實驗

選用兩個例子測試本文提出的面向高階自適應(yīng)流動計算的流場插值方法,第一個例子通過具有解析解的流場驗證該算法的插值誤差,第二個例子則將該算法應(yīng)用于高階自適應(yīng)流動計算中.文中所有測試例子均在小型工作站上運行,計算機配置為,內(nèi)存容量:16 GB;CPU 型號:六核Inter Core i7-8700,主頻:3.2 GHz.

4.1 圓環(huán)模型

該模型為四分之一圓環(huán)模型[26],圓環(huán)內(nèi)圓半徑為2、外圓半徑為3 (如圖2 和圖3 所示).考慮非黏性超音速渦流流經(jīng)該模型,流場邊界條件為:內(nèi)圓邊界為入口邊界,在此處的馬赫數(shù)為Mi=2、質(zhì)量密度ρi=1.區(qū)域其他位置的流場變量(密度ρ、速度分量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)流場計算中,為延續(xù)上一迭代步的計算狀態(tài),需要在每一迭代步開始前執(zhí)行流場插值過程,以獲得該迭代步的初始解.本文模仿自適應(yīng)流場計算過程中的插值過程,即將定義在一套網(wǎng)格(背景網(wǎng)格)中的流場插值到另一套網(wǎng)格(當(dāng)前網(wǎng)格)中,不同的是此處定義在背景網(wǎng)格上的流場通過流場變量的解析式(6)～ 式(9)求出.為驗證本文插值算法的精度,除了從背景網(wǎng)格上獲得的插值解之外,在新網(wǎng)格上也通過流場變量的解析式計算精確解.隨后,針對每一個新網(wǎng)格上的網(wǎng)格單元,根據(jù)以下式子計算插值誤差

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

為實現(xiàn)上述插值誤差計算,設(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)格點數(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ù)越高,插值誤差越小;當(dāng)插值精度相同時,網(wǎng)格規(guī)模越大,插值誤差越小.

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

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

4.2 NACA 0012 模型

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

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

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

圖5 自適應(yī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

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

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

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

5 結(jié)論

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

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