基于各向異性非結(jié)構(gòu)網(wǎng)格的超聲速流動(dòng)自適應(yīng)計(jì)算

鄒建鋒，盛 東，邢 菲，張 帥，楊永健

（浙江大學(xué) 玉泉校區(qū)航空航天學(xué)院，浙江 杭州 310027）

0 引 言

對(duì)于具有各向異性特征的流場(chǎng)，采用各向異性網(wǎng)格自適應(yīng)技術(shù)，將會(huì)很大程度的降低計(jì)算量和提高流場(chǎng)結(jié)構(gòu)的捕捉精度。各向異性網(wǎng)格自適應(yīng)算法包含兩方面的內(nèi)容，即度量張量的離散定義和基于度量張量的網(wǎng)格生成。文獻(xiàn)［1－8］等對(duì)各向異性網(wǎng)格的創(chuàng)建技術(shù)以及度量張量的離散定義方法做了廣泛的研究和實(shí)踐；文獻(xiàn)［9－10］則將各向異性自適應(yīng)技術(shù)應(yīng)用于計(jì)算具有強(qiáng)間斷的多介質(zhì)界面流動(dòng)問(wèn)題，并對(duì)計(jì)算的可靠性和準(zhǔn)確性做了對(duì)比研究。國(guó)內(nèi)已有的相關(guān)研究主要關(guān)注各向異性網(wǎng)格的自動(dòng)生成工作［11－12］。

超聲速楔形體繞流及超聲速橫向噴流問(wèn)題是高超聲速飛行器流動(dòng)控制及超燃沖壓發(fā)動(dòng)機(jī)燃料噴射等實(shí)際問(wèn)題中的典型流動(dòng)模式，具有顯著的各向異性特征。如果能將各向異性網(wǎng)格自適應(yīng)求解技術(shù)應(yīng)用于這一類問(wèn)題的研究，則對(duì)于深入分析高速流動(dòng)問(wèn)題有著很大的推動(dòng)作用。本文的第1小節(jié)將對(duì)我們采用的各向異性網(wǎng)格自適應(yīng)技術(shù)進(jìn)行簡(jiǎn)要的介紹。第2小節(jié)給出了基于密度的可壓縮N－S方程的求解方法。第3和第4小節(jié)是本文的主要部分，給出了各向異性自適應(yīng)算法應(yīng)用于超聲速楔形體繞流及超聲速橫向噴流問(wèn)題計(jì)算的收斂性和準(zhǔn)確性分析。第5小節(jié)是全文的小結(jié)。

1 各向異性網(wǎng)格自適應(yīng)技術(shù)

基于各向異性非結(jié)構(gòu)網(wǎng)格的自適應(yīng)求解技術(shù)在國(guó)內(nèi)外相關(guān)文獻(xiàn)中有著廣泛的研究。本文主要參照文獻(xiàn)［13－14］中的算法介紹，開(kāi)展各向異性網(wǎng)格生成技術(shù)與自適應(yīng)求解過(guò)程的實(shí)現(xiàn)；主要研究定常流動(dòng)問(wèn)題，對(duì)應(yīng)的自適應(yīng)迭代求解過(guò)程描述如下。（注：自適應(yīng)求解之前，需要預(yù)先給定網(wǎng)格質(zhì)量期望值q與網(wǎng)格單元期望數(shù)目N＊）

1）計(jì)算域初始網(wǎng)格系統(tǒng)創(chuàng)建

應(yīng)用基于陣面推進(jìn)技術(shù)（AFT）的三角形網(wǎng)格生成程序，對(duì)計(jì)算域進(jìn)行三角形網(wǎng)格劃分Tk（k＝0），得到初始網(wǎng)格系統(tǒng)。

2）流場(chǎng)解獲取

應(yīng)用本文第2小節(jié)介紹的可壓縮問(wèn)題求解器，獲得Tk（k≥0）網(wǎng)格上的流場(chǎng)解φk。

3）度量張量矩陣計(jì)算

基于流場(chǎng)解φk，獲得描述流場(chǎng)誤差信息的Hessian矩陣Hk及其對(duì)應(yīng)的度量張量矩陣Mk，其數(shù)學(xué)表達(dá)式如下：

其中，｜Δ｜和p（Δ）分別是網(wǎng)格單元的面積與周長(zhǎng)，h＊＝h＊（N＊）是計(jì)算域網(wǎng)格單元的平均尺寸，所有幾何測(cè)度都是在度量空間Mk內(nèi)進(jìn)行。函數(shù)Ft（0≤F≤1）是極值點(diǎn)在t＝1處的凸函數(shù)。網(wǎng)格系統(tǒng)的質(zhì)量Qk定義為質(zhì)量最差的網(wǎng)格單元質(zhì)量值。

5）基于度量矩陣的網(wǎng)格重劃分

對(duì)計(jì)算域執(zhí)行新的三角劃分Tk＋1，以獲得在度量張量Mk定義的幾何空間內(nèi)為最優(yōu)的三角化系統(tǒng)，并返回第2步。

其中，φk是流場(chǎng)物理量的第k次迭代解，λi與V分別是Hessian矩陣Hk的特征值與特征矢量。本文應(yīng)用了基于Mach數(shù)場(chǎng)的度量張量矩陣。

4）網(wǎng)格系統(tǒng)質(zhì)量判斷

網(wǎng)格系統(tǒng)質(zhì)量達(dá)到預(yù)定值q，則自適應(yīng)循環(huán)結(jié)束。各向異性非結(jié)構(gòu)網(wǎng)格系統(tǒng)的質(zhì)量定義建立在度量張量Mk所定義的幾何空間之上。本文采用的單元質(zhì)量定義式為：

2 可壓縮流動(dòng)求解技術(shù)

需要提及的是，應(yīng)用各向異性網(wǎng)格系統(tǒng)不需要專門的流體力學(xué)求解器。本文采用的是作者所在研究小組開(kāi)發(fā)的并行流動(dòng)及燃燒求解器，它建立在格點(diǎn)型非結(jié)構(gòu)網(wǎng)格系統(tǒng)上，采用有限元形式的高階低耗散格式。

3 高速楔形體繞流問(wèn)題計(jì)算

美國(guó)國(guó)家應(yīng)用軟件研究中心CFD研究所（NPARC）的 CFD驗(yàn)證與確認(rèn)網(wǎng)站［15］，提供了高速楔形體繞流問(wèn)題（Ma＝2．5）的 Wind軟件計(jì)算結(jié)果。本文對(duì)該問(wèn)題進(jìn)行了各向異性網(wǎng)格自適應(yīng)求解，以驗(yàn)證該方法對(duì)高速流動(dòng)問(wèn)題的適應(yīng)性和優(yōu)越性。

圖1展示了自適應(yīng)求解過(guò)程中網(wǎng)格與馬赫數(shù)分布情況的演變與收斂歷程。初始網(wǎng)格的單元數(shù)為4866，節(jié)點(diǎn)數(shù)為2582。自適應(yīng)計(jì)算中，網(wǎng)格單元數(shù)N＊被設(shè)置為2000。

由圖1可以看出，自適應(yīng)求解的第一次迭代就產(chǎn)生顯著的效果，網(wǎng)格資源明顯向激波附近匯集，計(jì)算得到的激波變得尖銳。經(jīng)過(guò)四次迭代，自適應(yīng)求解過(guò)程趨于收斂，網(wǎng)格單元被有效分配至激波附近的狹窄區(qū)域。

圖2和圖3給出了楔形體繞流自適應(yīng)求解過(guò)程中馬赫數(shù)分布曲線的演變及其與分析解［16］和 Wind軟件數(shù)值解［15］的對(duì)比情況。Wind軟件在計(jì)算中設(shè)計(jì)了155×101的結(jié)構(gòu)網(wǎng)格系統(tǒng)。本文設(shè)計(jì)的各向異性網(wǎng)格系統(tǒng)的單元數(shù)為2000，約為Wind軟件采用的網(wǎng)格數(shù)目的13%，網(wǎng)格規(guī)模有了大幅度的降低。

圖1 楔形體繞流問(wèn)題的自適應(yīng)求解過(guò)程Fig．1 Adaptive solution for supersonic flow around a wedge

圖2 楔形體繞流自適應(yīng)求解過(guò)程馬赫數(shù)分布曲線Fig．2 Mach number along free－stream direction

圖3 馬赫數(shù)分布曲線的局部放大特征Fig．3 Zoom in on the flow around shock wave

可以看出，波前馬赫數(shù)在自適應(yīng)迭代中與分析解吻合的很好。經(jīng)過(guò)第三次迭代，波后馬赫數(shù)與分析解達(dá)到一致，捕捉到的激波區(qū)域?qū)挾扰cWind軟件基本一致。第四次迭代后，激波區(qū)域?qū)挾冗M(jìn)一步變窄并趨于收斂。Wind軟件的計(jì)算結(jié)果在激波附近出現(xiàn)了明顯的數(shù)值振蕩，而在各向異性網(wǎng)格自適應(yīng)求解中，數(shù)值振蕩被有效減弱。

表1是本文計(jì)算結(jié)果與理論分析解及Wind軟件模擬結(jié)果的定量對(duì)比?？梢钥闯?，本文的計(jì)算結(jié)果不但在定性上能捕捉到相對(duì)更加尖銳的激波效應(yīng)，與理論分析解的定量對(duì)比也吻合的很好，物理量誤差值與Wind軟件模擬結(jié)果相比有了明顯的下降。

表1 激波前后流場(chǎng)數(shù)據(jù)對(duì)比Table 1 Comparison of the present results with theory and Wind code

4 高速橫向噴流干擾問(wèn)題計(jì)算

為了研究各向異性非結(jié)構(gòu)網(wǎng)格自適應(yīng)算法在更加復(fù)雜的高速流動(dòng)問(wèn)題中的應(yīng)用，我們計(jì)算了超聲速橫向噴流問(wèn)題。

我們模擬的是聲速射流橫向噴入超聲速來(lái)流所形成的復(fù)雜流動(dòng)結(jié)構(gòu)。超聲速空氣來(lái)流壓力為3145Pa，來(lái)流馬赫數(shù)為3．5。橫向噴口位于超聲速來(lái)流進(jìn)口下游228．6mm處的位置，噴口尺寸為0．2667mm，射流壓力與空氣來(lái)流壓力之比為8．74。

本文創(chuàng)建的各向同性網(wǎng)格系統(tǒng)包含165226個(gè)三角形單元，合計(jì)84280個(gè)節(jié)點(diǎn)。而各向異性網(wǎng)格自適應(yīng)求解過(guò)程，我們把計(jì)算網(wǎng)格系統(tǒng)的單元規(guī)模N＊限定為5000。

圖4給出了自適應(yīng)求解過(guò)程中不同網(wǎng)格對(duì)應(yīng)的凈質(zhì)量流率的收斂曲線?？梢钥闯?，在各向同性網(wǎng)格系統(tǒng)上，即初始網(wǎng)格，凈流率的收斂速度較為緩慢，迭代4000余次凈質(zhì)量流率達(dá)到5．0×10－4的量級(jí)；而對(duì)于各向異性網(wǎng)格系統(tǒng)，迭代1000次后就能達(dá)到收斂，而且凈質(zhì)量流率為1．0×10－6的量級(jí)。

圖4 不同網(wǎng)格系統(tǒng)上的收斂曲線對(duì)比Fig．4 Mass flow rate vs．a(chǎn)daption process

圖5 是四次自適應(yīng)計(jì)算過(guò)程后得到的網(wǎng)格系統(tǒng)及噴口附近流動(dòng)情況?？梢园l(fā)現(xiàn)，在自適應(yīng)過(guò)程中，計(jì)算區(qū)域的網(wǎng)格分布逐漸向激波間斷等大梯度的局部區(qū)域集中，從而極大的提升了大梯度區(qū)域流場(chǎng)結(jié)構(gòu)的分辨率。因此，各向異性網(wǎng)格自適應(yīng)技術(shù)具有對(duì)有限的網(wǎng)格節(jié)點(diǎn)資源進(jìn)行合理有效布置的功能。

圖5 自適應(yīng)求解得到的網(wǎng)格與馬赫數(shù)分布Fig．5 Adaptive solution mesh and flow features

然而從圖5可以發(fā)現(xiàn)，湍流邊界層附近的流場(chǎng)沒(méi)能被準(zhǔn)確模擬和識(shí)別。也就是說(shuō)，純粹基于Mach數(shù)場(chǎng)的度量張量矩陣，不能有效定義和描述邊界層內(nèi)網(wǎng)格單元的尺寸分布及形狀朝向。因此，對(duì)基于Mach數(shù)場(chǎng)的度量張量矩陣進(jìn)行合理的改進(jìn)使之能準(zhǔn)確有效模擬邊界層內(nèi)的流動(dòng)，這是對(duì)超聲速橫向噴流問(wèn)題及其他更加復(fù)雜的高速流動(dòng)問(wèn)題進(jìn)行各向異性網(wǎng)格自適應(yīng)求解需要解決的關(guān)鍵技術(shù)難點(diǎn)。

5 存在的問(wèn)題和進(jìn)一步工作

本文對(duì)基于Hessian度量張量矩陣的各向異性非結(jié)構(gòu)網(wǎng)格的創(chuàng)建，以及對(duì)基于各向異性非結(jié)構(gòu)網(wǎng)格的自適應(yīng)求解算法的實(shí)現(xiàn)，做了詳細(xì)的介紹。并通過(guò)對(duì)高速楔形體繞流問(wèn)題和高速橫向噴流問(wèn)題的自適應(yīng)計(jì)算，展現(xiàn)各向異性非結(jié)構(gòu)網(wǎng)格技術(shù)在準(zhǔn)確模擬激波流場(chǎng)方面的獨(dú)特優(yōu)勢(shì)以及存在的困難和問(wèn)題。

本文研究表明，基于各向異性非結(jié)構(gòu)網(wǎng)格的自適應(yīng)求解具有如下幾個(gè)方面的優(yōu)點(diǎn)：

1）具有對(duì)網(wǎng)格單元資源進(jìn)行合理的分配和布置的能力，有效降低問(wèn)題的規(guī)模；

2）提升大梯度流動(dòng)區(qū)域的計(jì)算分辨率；

3）有效降低激波附近流場(chǎng)的數(shù)值振蕩；

4）提升激波前后流場(chǎng)物理量的計(jì)算精度；

5）具有良好的收斂性和收斂精度。

但在對(duì)高速橫向噴流問(wèn)題的計(jì)算中我們也發(fā)現(xiàn)，純粹基于Mach數(shù)分布場(chǎng)的度量矩陣，不能準(zhǔn)確描述湍流邊界層內(nèi)部的流動(dòng)信息，尤其是對(duì)橫向噴口附近的局部小渦結(jié)構(gòu)以及激波／邊界層的相互作用區(qū)域缺乏準(zhǔn)確模擬的能力。因此，還需要對(duì)原有的度量張量矩陣進(jìn)行合理的局部修改和完善，以滿足上述流場(chǎng)特性的捕捉和計(jì)算。

［1］ HABASHI W，DOMPIERRE J，BOURGAULT Y．Anisotropic mesh adaptation：towards user－independent，mesh－independent and solver－independent CFD［J］．International Journal for Numerical Methods in Fluids，2000，32（6）：725－744．

［2］ DOMPIERRE J，VALLET M，BOURGAULT Y．Anisotropic mesh adaptation：towards user－independent，mesh－independent and solver－independent CFD［J］．International Journal for Numerical Methods in Fluids，2002，39（8）：675－702．

［3］ TAM A，AIT－ALI－YAHIA D，ROBICHAUD M．Anisotropic mesh adaptation for 3Dflows on structured and unstructured grids［J］．Computer Methods in Applied Mechanics and Engineering，2000，189（4）：1205－1230．

［4］ DOLEJSI V．Anisotropic mesh adaption technique for viscous flow simulationp［J］．East West Journal of Numerical Mathematics，2001，9（1）：1－24．

［5］ BOTTASSO C．Anisotropic mesh adaption by metricdriven optimization［J］．International Journal for Numeri－cal Methods in Engineering，2004，60（3）：597－639．

［6］ DOLEJSI V，F(xiàn)ELCMAN J．Anisotropic mesh adaptation for numerical solution of boundary value problems［J］．Numerical Methods for Partial Differential E－quations，2004，20（4）：576－608．

［7］ FREY P，ALAUZET F．Anisotropic mesh adaptation for CFD computations［J］．Computer Methods in Applied Mechanics and Engineering，2005，194（48－49）：5068－5082．

［8］ REMACLE J，LI X，SHEPHARD M．Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods［J］．International Journal for Numerical Methods in Engineering，2005，62（7）：899－923．

［9］ BELHAMADIA Y，F(xiàn)ORTIN A，CHAMBERLAND E．Anisotropic mesh adaptation for the solution of the Stefan problem［J］．Journal of Computational Physics，2004，194（1）：233－255．

［10］BELHAMADIA Y，F(xiàn)ORTIN A，CHARNBERLAND E．Three－dimensional anisotropic mesh adaptation for phase change problems［J］．Journal of Computational Physics，2004，201（2）：753－770．

［11］徐明海，陶文銓．二維各向異性非結(jié)構(gòu)化網(wǎng)格的自動(dòng)生成與應(yīng)用［J］．西安交通大學(xué)學(xué)報(bào)，2002，3（3）：222－235．

［12］趙建軍，鐘毅芳，周濟(jì)．各向異性網(wǎng)格生成及其在曲面三角化中的應(yīng)用［J］．中國(guó)圖像圖形學(xué)報(bào)，2002，7（9）：950－955．

［13］VASSILEVSKI Y，DYADECHKO V，LIPNIKOV K．Hessian－based anisotropic mesh adaptation in domains with discrete boundaries［J］．Russian Journal of Numerical Analysis and Mathematical Modelling，2005，20（4）：391－402．

［14］LIPNIKOV K，VASSILEVSKI Y．Analysis of Hessian recovery methods for generating adaptive meshes［R］．In：IMR．Los Alamos National Laboratory Institute of Numerical Mathematics，2006，163－171．

［15］http：／／www．grc．nasa．gov／WWW／wind／valid／archive［EB］．html．

［16］ANDERSON J．Modern compressible flow［M］．New York：McGraw Hill Inc，1982．