崔鵬程， 鄧有奇， 唐靜， 李彬

中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所， 綿陽 621000

基于伴隨方程的網(wǎng)格自適應(yīng)及誤差修正

崔鵬程， 鄧有奇， 唐靜*， 李彬

中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所， 綿陽 621000

基于流場方程的離散伴隨優(yōu)化理論和三維非結(jié)構(gòu)網(wǎng)格，建立了網(wǎng)格自適應(yīng)技術(shù)和目標(biāo)函數(shù)誤差修正方法。詳細(xì)研究了用流動(dòng)伴隨變量進(jìn)行目標(biāo)函數(shù)的誤差估計(jì)和修正技術(shù)，構(gòu)造了適用于格心格式有限體積法的流場變量插值技術(shù)和網(wǎng)格單元剖分判據(jù)，初步實(shí)現(xiàn)了網(wǎng)格物面投影和空間單元優(yōu)化，發(fā)展了適用于有限體積法的整套網(wǎng)格自適應(yīng)方法。對NACA0012翼型和ONERA-M6機(jī)翼繞流進(jìn)行了自適應(yīng)數(shù)值模擬，并對升、阻力等目標(biāo)函數(shù)進(jìn)行了誤差修正。數(shù)值結(jié)果表明，本文自適應(yīng)方法能正確地捕捉到影響目標(biāo)函數(shù)計(jì)算精度的敏感區(qū)域，網(wǎng)格自適應(yīng)和誤差修正兩項(xiàng)技術(shù)顯著提高了升、阻力等目標(biāo)函數(shù)的計(jì)算精度。

非結(jié)構(gòu)網(wǎng)格； 伴隨方程； 目標(biāo)函數(shù)； 誤差修正； 網(wǎng)格自適應(yīng)

飛行器氣動(dòng)力的高精度預(yù)測是計(jì)算流體力學(xué)一直以來追求的目標(biāo)。計(jì)算精度既取決于所使用的數(shù)值格式精度，也取決于計(jì)算網(wǎng)格的質(zhì)量。網(wǎng)格質(zhì)量不僅影響計(jì)算精度，還影響計(jì)算收斂性。人為地全局加密和預(yù)測流場特征并局部加密是兩種獲得足夠網(wǎng)格分辨率的常用方法。對于復(fù)雜流場或不熟悉的流場，往往通過全局加密或試算的方法來解決。全局加密不僅浪費(fèi)計(jì)算資源，還會(huì)延長計(jì)算時(shí)間。預(yù)測流場特征需要大量的工程實(shí)踐經(jīng)驗(yàn)，并且面對新型復(fù)雜飛行器仍需不斷的摸索和學(xué)習(xí)。網(wǎng)格自適應(yīng)技術(shù)是一種自動(dòng)優(yōu)化計(jì)算網(wǎng)格、提高流動(dòng)模擬精度的有效方法。

近年來局部網(wǎng)格加密的自適應(yīng)技術(shù)在非結(jié)構(gòu)網(wǎng)格中得到廣泛研究?，F(xiàn)有的大多數(shù)自適應(yīng)方法都基于流場特征進(jìn)行自適應(yīng)加密，如激波、邊界層、旋渦、尾流、滑移線或駐點(diǎn)等?；诹鲌鎏卣鞯木W(wǎng)格自適應(yīng)方法假設(shè)局部誤差與流動(dòng)特征量相關(guān)，利用流場特征探測器(如梯度探測器、熵增探測器、法向馬赫數(shù)探測器等)標(biāo)識需要加密網(wǎng)格的流場區(qū)域，通過局部網(wǎng)格加密減小流場局部誤差。該方法致力于減小流場特征處的誤差，忽略了影響目標(biāo)函數(shù)全局誤差的其他區(qū)域。該自適應(yīng)方法只捕捉流場特征，不一定能提高流場的全局計(jì)算精度[1]。

基于伴隨方程的網(wǎng)格自適應(yīng)方法[1-2]，更能有效地提高全機(jī)氣動(dòng)特性模擬精度。該方法通過伴隨方程直接建立局部誤差與目標(biāo)函數(shù)全局誤差之間的關(guān)系，可直接度量局部誤差對全局誤差的影響，對目標(biāo)函數(shù)的全局誤差進(jìn)行估計(jì)和修正。該方法在設(shè)定目標(biāo)函數(shù)全局誤差上限后，通過自適應(yīng)加密對全局誤差產(chǎn)生較大影響的網(wǎng)格區(qū)域，使計(jì)算誤差在網(wǎng)格單元間均勻分布并減小到給定的誤差上限以下。相比基于流動(dòng)特征的自適應(yīng)技術(shù)，達(dá)到相同的計(jì)算準(zhǔn)度，該方法增加的網(wǎng)格量更少，能構(gòu)造基于計(jì)算精度的自適應(yīng)終止判據(jù)，對飛行器的氣動(dòng)設(shè)計(jì)更有指導(dǎo)意義[3-4]。

國內(nèi)外很多學(xué)者對基于伴隨方程的網(wǎng)格自適應(yīng)和誤差修正技術(shù)進(jìn)行了研究。Park[5-8]發(fā)展了基于伴隨方程的三維各項(xiàng)同性網(wǎng)格自適應(yīng)方法，并探索了混合網(wǎng)格的自適應(yīng)和誤差修正技術(shù)；Park等[9]發(fā)展了結(jié)構(gòu)網(wǎng)格的伴隨自適應(yīng)與誤差估計(jì)理論，對音爆問題進(jìn)行了探討；Oliver和James等[10-11]將伴隨自適應(yīng)與高階算法相結(jié)合，大幅提高了目標(biāo)函數(shù)計(jì)算的準(zhǔn)確性；Yamaleev和Diskin[12-13]針對非定常問題研究了二維各向異性非結(jié)構(gòu)網(wǎng)格的伴隨自適應(yīng)和誤差修正；國內(nèi)楊振虎和周磊[14]探索了二維非結(jié)構(gòu)網(wǎng)格的自適應(yīng)技術(shù)，楊夏勰和周春華[15]發(fā)展了一種多目標(biāo)函數(shù)修正的二維非結(jié)構(gòu)網(wǎng)格自適應(yīng)方法。

但是現(xiàn)有的大部分文獻(xiàn)主要研究格點(diǎn)格式有限體積法的誤差修正及自適應(yīng)方法，對格心格式有限體積法的誤差修正和自適應(yīng)方法研究較少。格心格式有限體積法在邊界條件處理和基于面的循環(huán)遍歷等方面更有優(yōu)勢，在工程中應(yīng)用廣泛。國內(nèi)學(xué)者對基于伴隨方程的目標(biāo)函數(shù)誤差修正方法研究較少，也無人介紹基于伴隨方程的誤差修正方法在三維非結(jié)構(gòu)網(wǎng)格的應(yīng)用。

本文基于三維非結(jié)構(gòu)網(wǎng)格和離散伴隨方程理論，發(fā)展了一種適用于格心格式有限體積法的網(wǎng)格自適應(yīng)方法和升力、阻力等目標(biāo)函數(shù)誤差修正方法。經(jīng)過數(shù)值試驗(yàn)驗(yàn)證，該方法能明顯提高目標(biāo)函數(shù)的計(jì)算精度。

1 數(shù)值方法

1.1 流場方程

本文使用的流場解算器為自主研發(fā)的MFlow軟件[16]，時(shí)間離散使用LU-SGS (Lower-Upper Symmetric Gauss-Seidel)方法[17]，如無特殊說明，無黏通量離散使用Roe格式[18]。本文針對無黏流動(dòng)研究純四面體網(wǎng)格的自適應(yīng)和誤差修正方法，控制方程采用守恒型三維歐拉方程,即

(1)

式中：U為流動(dòng)解向量;F、G、H為通量項(xiàng)。在定常情況下?U/?t=0，歐拉方程簡化為

(2)

1.2 伴隨方程

得到流場方程的離散解之后，目標(biāo)函數(shù)f(如升力、阻力和力矩等)可以表示為流場離散解U的函數(shù)，f=f(U)，目標(biāo)函數(shù)對流場變量求導(dǎo)，由鏈?zhǔn)椒▌t得

(3)

式中：R為流場方程的殘差;定義伴隨變量Ψ為目標(biāo)函數(shù)對離散殘差的敏感度[6],即

(4)

綜合式(3)與式(4)，得到伴隨方程為

(5)

文中伴隨方程求解采用LU-SGS方法[19]和多重網(wǎng)格技術(shù)[20]，數(shù)值離散格式與流場方程離散格式保持一致。

1.3 誤差修正

假設(shè)存在兩套不同的網(wǎng)格，粗網(wǎng)格ΩH和細(xì)網(wǎng)格Ωh，粗網(wǎng)格和細(xì)網(wǎng)格上計(jì)算的目標(biāo)函數(shù)分別記為fH(UH)和fh(Uh)，H與h分別為粗網(wǎng)格和細(xì)網(wǎng)格的尺度。粗網(wǎng)格消耗的計(jì)算資源較少、計(jì)算時(shí)間也較短，但計(jì)算得到的目標(biāo)函數(shù)準(zhǔn)確度不高。細(xì)網(wǎng)格上計(jì)算的目標(biāo)函數(shù)準(zhǔn)確度很高，但是往往消耗大量的計(jì)算資源及計(jì)算時(shí)間。通過對fH(UH)進(jìn)行修正獲得fh(Uh)的近似值是可以期待的[5,21]。本文用伴隨方程修正目標(biāo)函數(shù)的方法來自文獻(xiàn)[5,21]，簡述如下。

(6)

(7)

把式(7)代入到式(6)中，得到細(xì)網(wǎng)格上目標(biāo)函數(shù)的估計(jì)值：

(8)

引入伴隨方程：

(9)

則式(8)可寫為

(10)

(11)

(12)

(13)

1.4 自適應(yīng)參數(shù)

假設(shè)插值足夠精確，則式(11)對線性目標(biāo)函數(shù)的修正是精確的，但在實(shí)際應(yīng)用中目標(biāo)函數(shù)往往是非線性的。對非線性目標(biāo)函數(shù)，式(11)還有一定剩余誤差，希望通過網(wǎng)格自適應(yīng)減少目標(biāo)函數(shù)修正后的剩余誤差，故將剩余誤差設(shè)為自適應(yīng)參數(shù)。自適應(yīng)參數(shù)的設(shè)置方法來自參考文獻(xiàn)[6,21]，簡述如下。將式(10)進(jìn)一步展開整理得

(14)

則目標(biāo)函數(shù)修正后的剩余誤差為

(15)

式(15)中剩余誤差是伴隨解的誤差和流場方程殘差的內(nèi)積。根據(jù)伴隨方程進(jìn)一步推導(dǎo)[6]可得另一種形式的剩余誤差，如式(16)可以表示成伴隨方程殘差和流場方程解誤差的內(nèi)積。

(16)

(17)

(18)

(19)

把這兩種形式的剩余誤差絕對值的均值設(shè)為自適應(yīng)參數(shù)，則自適應(yīng)參數(shù)ε表示為

(20)

2 自適應(yīng)策略

2.1 粗網(wǎng)格剖分方法

基于伴隨方程的誤差估計(jì)和網(wǎng)格自適應(yīng)方法需要構(gòu)造一套密網(wǎng)格，本文通過把粗網(wǎng)格全局剖分細(xì)化得到內(nèi)嵌細(xì)網(wǎng)格，從而修正全局目標(biāo)函數(shù)和計(jì)算自適應(yīng)參數(shù)。四面體網(wǎng)格單元的剖分方式有一分二、一分四、一分八等[22]。

如圖1所示，四面體一分二和一分四方法得到的小四面體與原四面體不相似，若原四面體質(zhì)量較差，剖分后得到的小四面體質(zhì)量將更差。為了使剖分后的網(wǎng)格盡可能地保持原四面體的形狀，本文中，統(tǒng)一采取一分八的網(wǎng)格剖分方法，即在每個(gè)四面體的邊上添加一個(gè)新的網(wǎng)格點(diǎn)，將原四面體剖分為8個(gè)小四面體。四面體一分八由于點(diǎn)的連接方式不同，有3種不同的剖分方式，本文選取使內(nèi)嵌細(xì)網(wǎng)格的最大二面角最小的一種剖分方法。

圖1 粗網(wǎng)格剖分方法Fig.1 Methods to divide coarse grid

2.2 插值技術(shù)

在修正目標(biāo)函數(shù)和計(jì)算自適應(yīng)參數(shù)時(shí)，需計(jì)算粗網(wǎng)格流場解和伴隨解在細(xì)網(wǎng)格上的插值解。因此，一種準(zhǔn)確的插值方法至關(guān)重要。很多文獻(xiàn)采用格點(diǎn)格式有限體積法，流場變量和伴隨變量存儲(chǔ)在網(wǎng)格節(jié)點(diǎn)處，因此大部分文獻(xiàn)采用節(jié)點(diǎn)插值法[5-6]。本文使用格心格式有限體積法，流場變量和伴隨變量存儲(chǔ)在網(wǎng)格形心處[23]，若參照節(jié)點(diǎn)插值法，則細(xì)網(wǎng)格上插值解的準(zhǔn)確度不高。因此，本文提出了一種多項(xiàng)式插值法，直接將流場變量和伴隨變量從粗網(wǎng)格體心插值到細(xì)網(wǎng)格體心處，其基本思路如下。

對于低階插值，假設(shè)流場變量和伴隨變量在網(wǎng)格單元k內(nèi)線性分布，滿足線性方程：

q=a+bx+cy+dz

(21)

式(21)有4個(gè)未知系數(shù)a、b、c和d，如圖2所示，假設(shè)在與該網(wǎng)格單元共點(diǎn)的所有m個(gè)網(wǎng)格單元內(nèi)，流場變量和伴隨變量也滿足該線性多項(xiàng)式，則會(huì)形成m個(gè)線性方程：

qi=a+bxi+cyi+dzii=1，2，…，m

(22)

再用最小二乘法求解該線性方程組，得到各項(xiàng)系數(shù)。由于距離的不同，網(wǎng)格單元k周邊的網(wǎng)格單元對其影響不同，故求解方程組時(shí)對各網(wǎng)格單元所對應(yīng)的方程設(shè)置不同的權(quán)重，離網(wǎng)格單元k越近，權(quán)重越大。設(shè)置權(quán)重wr為

wr=1/l2

(23)

式中：l為各網(wǎng)格單元和網(wǎng)格k形心間的距離。得到多項(xiàng)式各系數(shù)后，把各細(xì)網(wǎng)格形心坐標(biāo)代入線性方程，就可以得到細(xì)網(wǎng)格上的低階插值。

同理對于高階插值，假設(shè)流場變量和伴隨變量在網(wǎng)格單元k內(nèi)呈二次分布，滿足二次方程：

q=a+bx+cy+dz+ex2+fy2+

gz2+hxy+lxz+myz

(24)

式(24)有10個(gè)未知系數(shù)，與低階插值類似，假設(shè)在與該網(wǎng)格單元共點(diǎn)的所有m個(gè)網(wǎng)格單元內(nèi)，流場變量和伴隨變量也滿足該二次多項(xiàng)式，則會(huì)形成m個(gè)二次方程：

(25)

式中：j=1，2，…，m。

圖2 二維多項(xiàng)式插值法Fig.2 Two-dimensional interpolation method using polynomial

與線性插值類似，用最小二乘法求解帶距離權(quán)的二次方程組，可得細(xì)網(wǎng)格上的高階插值。

2.3 自適應(yīng)判據(jù)

自適應(yīng)的目標(biāo)是降低目標(biāo)函數(shù)修正后的剩余誤差，即減小計(jì)算域內(nèi)自適應(yīng)參數(shù)的值。內(nèi)嵌在粗網(wǎng)格內(nèi)的細(xì)網(wǎng)格構(gòu)造好以后，可以方便計(jì)算每個(gè)粗網(wǎng)格單元k的自適應(yīng)參數(shù)εk：

(26)

式中：εk為粗網(wǎng)格單元k內(nèi)的內(nèi)嵌細(xì)網(wǎng)格單元自適應(yīng)參數(shù)之和;lk為粗網(wǎng)格單元k的細(xì)網(wǎng)格單元集。定義局部自適應(yīng)準(zhǔn)則：

(27)

此外，定義全局自適應(yīng)參數(shù)ε為所有粗網(wǎng)格自適應(yīng)參數(shù)之和，即ε=∑εk。ε雖然并不是準(zhǔn)確的目標(biāo)函數(shù)的全局誤差，但是ε是全局誤差的一種度量[6]。由此，一種全局的自適應(yīng)判據(jù)可以定義為

(28)

當(dāng)ηz>1，即粗網(wǎng)格上的目標(biāo)函數(shù)全局誤差大于人為設(shè)定的誤差上限，網(wǎng)格需要進(jìn)一步細(xì)化。與文獻(xiàn)[6]類似，本文結(jié)合ηk和ηz，設(shè)定綜合自適應(yīng)判據(jù)，對每個(gè)粗網(wǎng)格單元k，判斷ηzηk，若ηzηk>1，則標(biāo)記該網(wǎng)格單元k需要細(xì)化，若ηzηk<1，則該網(wǎng)格單元k無需細(xì)化。

2.4 網(wǎng)格優(yōu)化

對于物面上的粗網(wǎng)格單元，剖分成細(xì)網(wǎng)格后新添加的網(wǎng)格點(diǎn)通常不在物面上，需要將新網(wǎng)格點(diǎn)投影到物面上以保持網(wǎng)格與物面的相容性。

如圖3所示，對于物面網(wǎng)格，本文利用兩節(jié)點(diǎn)X1、X2及其法向與切向信息N1、N2，r1、r2構(gòu)造出物面曲線的二次擬合曲線，將自適應(yīng)新加入的點(diǎn)XNEW移動(dòng)到擬合曲線的中點(diǎn)。物面網(wǎng)格投影之后，用距離權(quán)函數(shù)方法[24]移動(dòng)相應(yīng)的空間網(wǎng)格，保持空間網(wǎng)格與物面網(wǎng)格的相容性。

圖3 物面投影方法Fig.3 Projection method on surface geometry

圖4 Laplacian網(wǎng)格光滑方法Fig.4 Laplacian method to smooth grid

自適應(yīng)后的網(wǎng)格質(zhì)量檢測及網(wǎng)格光滑對于自適應(yīng)技術(shù)也起到重要的作用。通常，由自適應(yīng)得到的網(wǎng)格局部質(zhì)量不高，需要光滑處理。本文采用Laplacian光滑方法[22]，如圖4所示，Laplacian光滑方法是將網(wǎng)格點(diǎn)Vi移動(dòng)到相鄰網(wǎng)格點(diǎn)(V1,V2，…，VN)的形心Vc處，使網(wǎng)格布局均勻分布，通過網(wǎng)格光滑可以改善自適應(yīng)網(wǎng)格局部疏密不均的情況。

2.5 自適應(yīng)流程

基于伴隨方程的自適應(yīng)是一個(gè)迭代過程，圖5 為其流程圖。首先把初始粗網(wǎng)格全局剖分細(xì)化得到內(nèi)嵌的細(xì)網(wǎng)格，并構(gòu)造高階插值與低階插值。然后計(jì)算局部自適應(yīng)參數(shù)與全局自適應(yīng)參數(shù)，對目標(biāo)函數(shù)進(jìn)行修正。最后判斷是否滿足自適應(yīng)判據(jù)，若滿足則自適應(yīng)過程結(jié)束，若不滿足，則標(biāo)記需要加密的粗網(wǎng)格并加密，把自適應(yīng)網(wǎng)格優(yōu)化之后進(jìn)行下一輪自適應(yīng)過程，直到滿足自適應(yīng)判據(jù)為止。

圖5 自適應(yīng)迭代過程Fig.5 Iteration process of adaptation

3 算例分析

本文采用NACA0012翼型和ONERA-M6機(jī)翼對所建立的非結(jié)構(gòu)網(wǎng)格自適應(yīng)方法進(jìn)行了驗(yàn)證，控制方程為歐拉方程，網(wǎng)格為純四面體非結(jié)構(gòu)網(wǎng)格。為了校驗(yàn)網(wǎng)格自適應(yīng)和目標(biāo)函數(shù)誤差修正的效果，本文將一系列均勻加密的網(wǎng)格與自適應(yīng)網(wǎng)格作對比，把均勻加密網(wǎng)格計(jì)算結(jié)果的Richardson外插值或極其精密網(wǎng)格的計(jì)算結(jié)果當(dāng)做期望值。

3.1 NACA0012翼型

NACA0012翼型擾流流場來流條件為：來流馬赫數(shù)Ma∞=0.8,迎角α=1.25°，來流溫度T∞=279.3 K,雷諾數(shù)Re=9×106。目標(biāo)函數(shù)為阻力系數(shù)CD，每次自適應(yīng)過程將誤差上限設(shè)為上次自適應(yīng)阻力系數(shù)的10%。

圖6是NACA0012翼型原始網(wǎng)格與自適應(yīng)后的對稱面網(wǎng)格分布，自適應(yīng)網(wǎng)格本身就反應(yīng)了一定的流場特征，自適應(yīng)后的網(wǎng)格在機(jī)翼前緣、后緣和機(jī)翼上下表面的強(qiáng)弱激波處都有明顯的加密，因?yàn)檫@些網(wǎng)格區(qū)域?qū)ACA0012翼型的阻力系數(shù)影響較大。

圖6 NACA0012翼型自適應(yīng)前后網(wǎng)格Fig.6 NACA0012 airfoil grid before and after adaptation

圖7是NACA0012翼型原始網(wǎng)格與伴隨自適應(yīng)網(wǎng)格計(jì)算得到的壓力分布p，從原始網(wǎng)格的壓力分布可以看出，機(jī)翼上表面捕捉到的強(qiáng)激波范圍較大，機(jī)翼下表面基本捕捉不到弱激波。而自適應(yīng)后網(wǎng)格計(jì)算得出的上表面強(qiáng)激波與機(jī)翼下表面弱激波都很精細(xì)，壓力等值線非常光滑。伴隨自適應(yīng)網(wǎng)格精確地捕捉到了NACA0012機(jī)翼的流場特征。

圖8給出了NACA0012翼型的阻力系數(shù)隨自適應(yīng)次數(shù)與網(wǎng)格量的變化曲線。隨著自適應(yīng)次數(shù)的遞增，網(wǎng)格量不斷加大，阻力系數(shù)逐漸收斂到外插值。從圖8可以看出，阻力系數(shù)誤差修正效果明顯，相同網(wǎng)格量，經(jīng)過誤差修正后阻力系數(shù)的計(jì)算準(zhǔn)確度比全局加密網(wǎng)格的準(zhǔn)確度提高了15.50%。

圖7 NACA0012翼型自適應(yīng)前后壓力分布Fig.7 Pressure distribution of NACA0012 airfoil before and after adaptation

圖8 NACA0012翼型阻力系數(shù)自適應(yīng)Fig.8 Drag coefficient adaption of NACA0012 airfoil

3.2 ONERA-M6機(jī)翼

1) 目標(biāo)函數(shù)為阻力系數(shù)

目標(biāo)函數(shù)為阻力系數(shù)CD，每一次自適應(yīng)過程將誤差上限設(shè)為上一次自適應(yīng)過程阻力系數(shù)的10%。計(jì)算狀態(tài)選擇為：Ma∞=0.84,α=3.06°，T∞=255.56 K,Re=1.172×107。在該計(jì)算狀態(tài)下，機(jī)翼表面為附著流，機(jī)翼上表面有λ型激波。

圖9是ONERA-M6機(jī)翼自適應(yīng)前后的上表面網(wǎng)格分布與空間網(wǎng)格分布，可以看出，自適應(yīng)后的網(wǎng)格在機(jī)翼前緣、后緣和激波處有明顯的加密，在遠(yuǎn)離機(jī)翼表面的區(qū)域網(wǎng)格加密較少。離機(jī)翼表面較遠(yuǎn)的網(wǎng)格區(qū)域?qū)ψ枇ο禂?shù)的貢獻(xiàn)較小，機(jī)翼前后緣以及激波等網(wǎng)格區(qū)域?qū)ψ枇ο禂?shù)貢獻(xiàn)很大。

圖9 ONERA-M6機(jī)翼阻力系數(shù)自適應(yīng)前后網(wǎng)格Fig.9 ONERA-M6 wing grid before and after adaptation for drag coefficient

圖10是文獻(xiàn)[5]中的自適應(yīng)網(wǎng)格，與本文的自適應(yīng)結(jié)果基本一致。

圖11是用ONERA-M6機(jī)翼原始網(wǎng)格與自適應(yīng)網(wǎng)格計(jì)算得到的壓力分布，從原始網(wǎng)格的壓力分布可以看出，機(jī)翼上表面的激波范圍較大，而從自適應(yīng)后的壓力分布可以看到精細(xì)的激波分布，壓力等值線非常光滑，這說明伴隨自適應(yīng)在流場特征處有明顯的加密，能較好地捕捉到流場特征。

圖10 文獻(xiàn)[5]中ONERA-M6機(jī)翼阻力系數(shù)自適應(yīng)后網(wǎng)格 Fig.10 Adapted ONERA-M6 wing grid for drag coefficient in Ref.[5]

圖11 ONERA-M6機(jī)翼阻力系數(shù)自適應(yīng)前后壓力分布Fig.11 Pressure distribution of ONERA-M6 wing before and after adaptation for drag coefficient

圖12給出了ONERA-M6機(jī)翼阻力系數(shù)隨自適應(yīng)次數(shù)和網(wǎng)格量的變化曲線，可以看出，誤差修正可以使阻力系數(shù)迅速地降到人工設(shè)定的誤差范圍以內(nèi)。阻力系數(shù)的誤差修正效果明顯，相同網(wǎng)格量，經(jīng)過誤差修正后阻力系數(shù)的計(jì)算準(zhǔn)確度比全局加密網(wǎng)格的準(zhǔn)確度提高了17.74%，修正后的阻力系數(shù)與外插值很接近。同時(shí)數(shù)值試驗(yàn)發(fā)現(xiàn)，自適應(yīng)網(wǎng)格比均勻加密網(wǎng)格計(jì)算收斂速度更快。由于自適應(yīng)網(wǎng)格捕捉到了一定的流場特征，在對目標(biāo)函數(shù)影響較大的區(qū)域進(jìn)行大量加密，使殘差更容易收斂，從而使目標(biāo)函數(shù)的收斂速度更快。

圖12 ONERA-M6機(jī)翼阻力系數(shù)自適應(yīng)Fig.12 Drag coefficient adaption of ONERA-M6 wing

2) 目標(biāo)函數(shù)為升力系數(shù)

此處采用的計(jì)算條件與目標(biāo)函數(shù)為阻力系數(shù)時(shí)一致，但目標(biāo)函數(shù)為升力系數(shù)CL。經(jīng)過數(shù)值試驗(yàn)發(fā)現(xiàn)，升力系數(shù)對網(wǎng)格不是很敏感，經(jīng)反復(fù)探索和參考國內(nèi)外文獻(xiàn)，每一次自適應(yīng)過程將誤差上限設(shè)為上一次自適應(yīng)過程升力系數(shù)的5%。

圖13是ONERA-M6機(jī)翼原始網(wǎng)格與自適應(yīng)后的上表面網(wǎng)格分布與空間網(wǎng)格分布，可以看出，自適應(yīng)后的網(wǎng)格在機(jī)翼前緣、后緣和激波處有明顯的加密，因?yàn)檫@些網(wǎng)格區(qū)域?qū)ιο禂?shù)貢獻(xiàn)很大。

圖13 ONERA-M6機(jī)翼升力系數(shù)自適應(yīng)前后網(wǎng)格Fig.13 ONERA-M6 wing grid before and after adaptation for lift coefficient

圖14給出了用ONERA-M6機(jī)翼原始網(wǎng)格與自適應(yīng)網(wǎng)格計(jì)算得到的壓力分布，從原始網(wǎng)格的壓力分布可以看出，機(jī)翼上表面的激波范圍較大，壓力等值線較粗糙，而從自適應(yīng)網(wǎng)格的壓力分布可以看到有精細(xì)的激波分布，壓力等值線非常光滑。伴隨自適應(yīng)能較好地捕捉到流場特征。

圖14 ONERA-M6機(jī)翼升力系數(shù)自適應(yīng)前后壓力分布Fig.14 Pressure distribution of ONERA-M6 wing before and after adaptation for lift coefficient

圖15為升力系數(shù)隨自適應(yīng)次數(shù)和網(wǎng)格量的變化曲線，把一套極其精密網(wǎng)格的計(jì)算結(jié)果當(dāng)做期望值來檢驗(yàn)誤差修正效果。相同網(wǎng)格量時(shí)升力系數(shù)的誤差修正效果明顯，經(jīng)過誤差修正后升力系數(shù)的計(jì)算準(zhǔn)確度比全局加密網(wǎng)格的準(zhǔn)確度提高了1.05%，誤差修正能使計(jì)算的升力系數(shù)快速的接近期望值。同時(shí)數(shù)值試驗(yàn)也發(fā)現(xiàn)，由于自適應(yīng)網(wǎng)格在對升力系數(shù)誤差敏感的區(qū)域進(jìn)行了自適應(yīng)加密，自適應(yīng)網(wǎng)格比均勻加密網(wǎng)格計(jì)算收斂速度更快。

圖15 ONERA-M6機(jī)翼升力系數(shù)自適應(yīng)Fig.15 Lift ocefficient adaption of ONERA-M6 wing

從圖12可以發(fā)現(xiàn)，在第4次自適應(yīng)之后用自適應(yīng)網(wǎng)格計(jì)算得出的阻力系數(shù)沒有相同數(shù)量全局加密網(wǎng)格計(jì)算得出的阻力系數(shù)精確。圖3為本文采用的物面網(wǎng)格投影方法示意圖，用兩節(jié)點(diǎn)及其法向構(gòu)造出物面曲線的二次擬合曲線，將自適應(yīng)新加入的點(diǎn)移動(dòng)到擬合曲線的中點(diǎn)。這種投影方法在普通物面能較為準(zhǔn)確地投影新添加的網(wǎng)格點(diǎn)，但是如圖16所示，在機(jī)翼前后緣等物面曲率較大的區(qū)域，這種投影方法在局部區(qū)域不能與物面完全匹配，自適應(yīng)后期由于投影誤差累積使得自適應(yīng)網(wǎng)格計(jì)算的阻力系數(shù)精度略低。數(shù)值模擬結(jié)果表明，這種投影方法對自適應(yīng)過程和誤差修正結(jié)果影響較小。

圖16 擬合曲線投影方法Fig.16 Fitted curves projection method

4 結(jié) 論

1) 對比相同數(shù)量的全局加密網(wǎng)格，經(jīng)過自適應(yīng)迭代后的網(wǎng)格分布更為合理，氣動(dòng)力計(jì)算精度得到明顯提高，為誤差修正打下良好基礎(chǔ)。

2) 采用本文誤差修正方法對升力、阻力系數(shù)等目標(biāo)函數(shù)進(jìn)行修正，氣動(dòng)力計(jì)算精度進(jìn)一步提高，并且隨自適應(yīng)迭代過程快速收斂到期望值。

3) 針對格心格式的有限體積方法，本文提出的流場變量和伴隨變量多項(xiàng)式插值方法適應(yīng)性好、插值精度高。

本文實(shí)現(xiàn)了純四面體非結(jié)構(gòu)網(wǎng)格的自適應(yīng)和誤差修正方法，下一步將開展混合網(wǎng)格(包含四面體、棱柱、金字塔等單元)的自適應(yīng)和誤差修正研究，并探索新的物面網(wǎng)格投影方法。

[1] PARK M A. Adjoint-based, three-dimensional error prediction and grid adaptation[J]. AIAA Journal, 2004, 42(9): 1854-1862.

[2] BECKER R, RANNACHER R. An optimal control approach to a posteriori error estimation in finite element methods[M]. 1st ed. Cambridge: Cambridge University Press, 2001: 1-102.

[3] FIDKOWSKI K J, DARMOFAL D L. Review of output-based error estimation and mesh adaptation in computational fluid dynamics[J]. AIAA Journal, 2011, 49(4): 673-694.

[4] GILES M B, SüLI E. Adjoint methods for PDEs: A posteriori error analysis and postprocessing by duality[J]. Acta Numerica, 2002, 11: 145-236.

[5] PARK M A. Adjoint-based, three-dimensional error prediction and grid adaptation: AIAA-2002-3286[R]. Reston: AIAA, 2002.

[6] PARK M A. Three-dimensional turbulent RANS adjoint-based error correction: AIAA-2003-3849[R]. Reston: AIAA, 2003.

[7] PARK M A, DARMOFAL D L. Output-adaptive tetrahedral cut-cell validation for sonic boom prediction: AIAA-2008-6594[R]. Reston: AIAA, 2008.

[8] PARK M A. Low boom configuration analysis with FUN3D adjoint simulation framework: AIAA-2011-3337[R]. Reston: AIAA, 2011.

[9] PARK M A, AFTOSMIS M A, CAMPBELL R C. Summary of the 2008 NASA fundamental aeronautics program sonic boom prediction workshop: AIAA-2013-0649[R]. Reston: AIAA, 2013.

[10] OLIVER T A, DARMOFAL D L. An unsteady adaptation algorithm for discontinuous Galerkin discretizations of the RANS equations: AIAA-2007-3940[R]. Reston: AIAA, 2007.

[11] JAMES M Y, MODISETTE M, DARMOFAL D L. The importance of mesh adaptation for higher-order discretizations of aerodynamic flows: AIAA-2011-3852[R]. Reston: AIAA, 2011.

[12] YAMALEEV N K, DISKIN B, PATHAK K. Error minimization via adjoint-based anisotropic grid adaptation: AIAA-2010-4436[R]. Reston: AIAA, 2010.

[13] DISKIN B, YAMALEEV N K. Grid adaptation using adjoint-based error minimization: AIAA-2011-3986[R]. Reston: AIAA, 2011.

[14] 楊振虎, 周磊. 基于誤差估計(jì)的伴隨網(wǎng)格自適應(yīng)方法[J]. 航空計(jì)算技術(shù), 2011, 41(3): 14-16.

YANG Z H, ZHOU L. Output-based error estimation and grid adaptive mesh refinement[J]. Aeronautical Computing Technique, 2011, 41(3): 14-16 (in Chinese).

[15] 楊夏勰, 周春華. 目標(biāo)函數(shù)誤差估計(jì)及網(wǎng)格自適應(yīng)處理[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2014, 32(5): 688-693.

YANG X X, ZHOU C H. Output-based error estimation and grid adaptation[J]. Acta Aerodynamica Sinica, 2014, 32(5): 688-693 (in Chinese).

[16] 張耀冰, 周乃春, 陳江濤. 小展弦比飛翼標(biāo)模雷諾數(shù)影響數(shù)值模擬研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2015, 33(3): 279-288.

ZHANG Y B, ZHOU N C, CHENG J T. Numerical investigation of Reynolds number effects on a low-aspect-ratio flying-wing model[J]. Acta Aerodynamica Sinica, 2015, 33(3): 279-288 (in Chinese).

[17] KIM J S, KWON O J. Improvement on block LU-SGS scheme for unstructured mesh Navier-Stokes computations: AIAA-2002-1061[R]. Reston: AIAA, 2002.

[18] ROE P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372.

[19] 李彬, 鄧有奇, 唐靜, 等. 基于三維非結(jié)構(gòu)網(wǎng)格的離散型伴隨方法[J]. 航空學(xué)報(bào), 2014, 35 (3): 674-686.

LI B, DENG Y Q, TANG J, et al. Discrete adjoint method for 3D on unstructured grid[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(3): 674-686 (in Chinese).

[20] 李彬, 唐靜, 鄧有奇, 等. 并行的多重網(wǎng)格方法在離散伴隨優(yōu)化中的應(yīng)用[J]. 航空學(xué)報(bào), 2014, 35(8): 2091-2101.

LI B, TANG J, DENG Y Q, et al. Application of parallel multigrid algorithm to discrete adjoint optimization[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(8): 2091-2101 (in Chinese).

[21] VENDITTI D A, DARMOFAL D L. Grid adaptation for functional outputs: application to two-dimensional inviscid flows[J]. Journal of Computational Physics, 2002, 176(1): 40-69.

[22] 唐靜, 鄭鳴, 鄧有奇, 等. 網(wǎng)格自適應(yīng)技術(shù)在復(fù)雜外形流場模擬中的應(yīng)用[J]. 計(jì)算力學(xué)學(xué)報(bào), 2015, 32(6): 752-757.

TANG J, ZHENG M, DENG Y Q, et al. Grid adaptation for simulation of complicated configuration[J]. Chinese Journal of Computational Mechanics, 2015, 32(6): 752-757 (in Chinese).

[23] BLAZEK J. Computational fluid dynamics: principles and applications[M]. 1st ed. Oxford: Elsevier, 2001: 75-120.

[24] 唐靜， 鄧有奇， 馬明生， 等. 飛翼氣動(dòng)優(yōu)化中參數(shù)化和網(wǎng)格變形技術(shù)研究[J]. 航空學(xué)報(bào), 2015, 36(5): 1480-1490.

TANG J, DENG Y Q, MA M S, et al. Parametrization and grid deformation techniques for fly-wing shape optimization[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(5): 1480-1490 (in Chinese).

崔鵬程男， 碩士， 研究實(shí)習(xí)員。主要研究方向： 非結(jié)構(gòu)網(wǎng)格技術(shù)， 飛行器氣動(dòng)外形優(yōu)化設(shè)計(jì)。

Tel: 0816-2463270

E-mail: ficojustdoit@163.com

鄧有奇男， 博士， 研究員， 博士生導(dǎo)師。主要研究方向： 計(jì)算空氣動(dòng)力學(xué)， 飛行器優(yōu)化設(shè)計(jì)。

Tel: 0816-2463007

E-mail: cai@cardc.cn

唐靜男， 博士研究生， 助理研究員。主要研究方向： 飛行器氣動(dòng)外形優(yōu)化設(shè)計(jì)， 非結(jié)構(gòu)網(wǎng)格技術(shù)。

Tel: 0816-2463270

E-mail: tangjingn@foxmail.com

李彬男， 博士， 助理研究員。主要研究方向： 飛行器氣動(dòng)外形優(yōu)化設(shè)計(jì)， 外掛物分離投放數(shù)值模擬。

Tel: 0816-2463091

E-mail: leebin2008@hotmail.com

URL：www.cnki.net/kcms/detail/11.1929.V.20160321.1035.002.html

Adjointequations-basedgridadaptationanderrorcorrection

CUIPengcheng，DENGYouqi，TANGJing*，LIBin

ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

Basedonthediscreteadjointoptimizingtheoryandthree-dimensionalunstructuredgrid,agridadaptationtechnologyandanerrorcorrectionmethodforobjectivefunctionarebuilt.Amethodtopredictandcorrecttheerrorofobjectivefunctionusingadjointequationsispresented.Then,aninterpolationtechnologywhichsuitsforcentre-basedfinitevolumemethodisproposed,somemethodstodividetetrahedralgrids,projectsurfacegridsandoptimizespatialgridsarediscussed,andacompletegridadaptationsystemwhichsuitsforfinitevolumemethodisbuilt.Finally,thegridadaptationmethodisappliedtothesimulationofinviscidflowsaroundNACA0012airfoilandONERA-M6wing,andtheerrorofobjectivefunction,suchasthecoefficientofdragandlift,iscorrected.Numericalresultsshowthatthesensitivegridsforobjectivefunctionaredetectedandrefinedbythisgridadaptationmethod,andtheaccuracyofobjectivefunctionisobviouslyimprovedaftergridadaptationanderrorcorrection.

unstructuredgrid;adjointequation;objectivefunction;errorcorrection;gridadaptation

2015-11-24；Revised2016-02-17；Accepted2016-03-14；Publishedonline2016-03-211035

.Tel.：0816-2463270E-mailtangjingn@foxmail.com

2015-11-24；退修日期2016-02-17；錄用日期2016-03-14； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-03-211035

www.cnki.net/kcms/detail/11.1929.V.20160321.1035.002.html

.Tel.：0816-2463270E-mailtangjingn@foxmail.com

崔鵬程， 鄧有奇， 唐靜， 等． 基于伴隨方程的網(wǎng)格自適應(yīng)及誤差修正J． 航空學(xué)報(bào),2016,37(10):2992-3002.CUIPC,DENGYQ,TANGJ,etal.Adjointequations-basedgridadaptationanderrorcorrectionJ.ActaAeronauticaetAstronauticaSinica,2016,37(10):2992-3002．

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0079

V211.3

A

1000-6893(2016)10-2992-11