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

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

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

0 引 言

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

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

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

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

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

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

2）流場解獲取

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

3）度量張量矩陣計算

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

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

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

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

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

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 可壓縮流動求解技術(shù)

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

3 高速楔形體繞流問題計算

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

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

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

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

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

圖2 楔形體繞流自適應(yīng)求解過程馬赫數(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)過第三次迭代，波后馬赫數(shù)與分析解達(dá)到一致，捕捉到的激波區(qū)域?qū)挾扰cWind軟件基本一致。第四次迭代后，激波區(qū)域?qū)挾冗M(jìn)一步變窄并趨于收斂。Wind軟件的計算結(jié)果在激波附近出現(xiàn)了明顯的數(shù)值振蕩，而在各向異性網(wǎng)格自適應(yīng)求解中，數(shù)值振蕩被有效減弱。

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

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

4 高速橫向噴流干擾問題計算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4）提升激波前后流場物理量的計算精度；

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

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

［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)格的自動生成與應(yīng)用［J］．西安交通大學(xué)學(xué)報，2002，3（3）：222－235．

［12］趙建軍，鐘毅芳，周濟(jì)．各向異性網(wǎng)格生成及其在曲面三角化中的應(yīng)用［J］．中國圖像圖形學(xué)報，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．