一種改進(jìn)的非定常激波裝配算法

常思源，白曉征，崔小強(qiáng)，劉君,*

1. 大連理工大學(xué) 航空航天學(xué)院，大連 116024 2. 中國人民解放軍95791部隊(duì)，酒泉 735018

在可壓縮非定常流動(dòng)中，激波的產(chǎn)生、傳播及其與幾何邊界的相互干擾，往往支配著整個(gè)流動(dòng)的特性。在計(jì)算流體力學(xué)(CFD)領(lǐng)域中，一些學(xué)者采用激波捕捉(Shock-Capturing, S-C)法研究激波與湍流/聲波的干擾時(shí)，發(fā)現(xiàn)激波處產(chǎn)生的非物理波動(dòng)會(huì)顯著降低下游流場的精度，即無論格式的設(shè)計(jì)精度有多高，激波下游的流場難以恢復(fù)到應(yīng)有的格式精度[1-2]。

本質(zhì)上講，這些缺陷的根源在于激波捕捉格式的基本構(gòu)造，即捕捉的激波等強(qiáng)間斷往往會(huì)跨越若干網(wǎng)格點(diǎn)，而這些節(jié)點(diǎn)上的參數(shù)僅僅是由人為引進(jìn)而來的，絲毫不能反映真實(shí)間斷內(nèi)部的流動(dòng)結(jié)構(gòu)[3]。因此近年來有很多學(xué)者在認(rèn)識(shí)到捕捉法的先天不足后，堅(jiān)信只有類似激波裝配(Shock-Fitting, S-F)的算法才能從根本上解決這些問題[2,4-5]。

傳統(tǒng)激波裝配法主要基于結(jié)構(gòu)網(wǎng)格框架，取得了一系列卓有成效的研究成果[6-8]。然而，考慮到其深受網(wǎng)格拓?fù)湎拗?、算法邏輯比較復(fù)雜及程序移植性差等缺點(diǎn)，因此近年來激波裝配法的研究熱點(diǎn)逐漸轉(zhuǎn)移到較為靈活的非結(jié)構(gòu)網(wǎng)格算法框架。經(jīng)過十幾年的發(fā)展，基于非結(jié)構(gòu)網(wǎng)格的激波裝配法已成功應(yīng)用于一些超/高超聲速、二/三維定常流動(dòng)問題[9-11]，然而將其拓展到復(fù)雜的非定常流動(dòng)中仍然困難重重，目前僅有一些探索性的嘗試。

2016年，Bonfiglioli等成功將其基于非結(jié)構(gòu)網(wǎng)格和格點(diǎn)型有限體積法的裝配算法[9]拓展到非定常計(jì)算并模擬了激波與渦的相互作用[12]，隨后又進(jìn)一步模擬了運(yùn)動(dòng)激波正規(guī)反射問題[13]；劉君等模擬了運(yùn)動(dòng)激波與氦氣柱的相互作用，采用裝配方法對(duì)氣柱界面進(jìn)行裝配，而運(yùn)動(dòng)激波仍用捕捉法進(jìn)行模擬[14]；隨后，鄒東陽等對(duì)等截面通道內(nèi)的運(yùn)動(dòng)激波進(jìn)行裝配，與理論值比較考察了算法的計(jì)算精度[15]；此外，作者所在團(tuán)隊(duì)采用邊界裝配策略對(duì)超高速彈丸發(fā)射時(shí)頭部大范圍運(yùn)動(dòng)的強(qiáng)激波進(jìn)行處理，拓展了裝配法在此類問題中的應(yīng)用[16]。然而，從以上研究來看，裝配的激波只在直壁面上運(yùn)動(dòng)，而對(duì)于激波在彎曲壁面上的傳播并未涉及；在整個(gè)計(jì)算中，裝配的激波陣面大多不能自發(fā)地實(shí)現(xiàn)大變形，如拉伸、縮短和彎曲，往往需要人為調(diào)整改變裝配激波點(diǎn)的分布；且當(dāng)激波運(yùn)動(dòng)致使流場網(wǎng)格質(zhì)量變差時(shí)，經(jīng)常需要人為重新替換網(wǎng)格，算法的自動(dòng)化程度比較低，計(jì)算代價(jià)比較大；最重要的一點(diǎn)，由于缺少成熟有效的基于非結(jié)構(gòu)網(wǎng)格的運(yùn)動(dòng)激波探測(cè)(Shock Detection)算法，目前尚不能自動(dòng)處理激波發(fā)生拓?fù)渥兓@類情況[12]，如激波的正規(guī)反射向馬赫反射演化、激波的新生與弱化等。

針對(duì)上述前3個(gè)問題，本文在近年來發(fā)展的自適應(yīng)間斷裝配求解器(Adaptive Discontinuity Fitting solver，ADFs)[17]的基礎(chǔ)上進(jìn)行了若干算法改進(jìn)，實(shí)現(xiàn)了裝配激波點(diǎn)在曲壁面上的運(yùn)動(dòng)；通過自動(dòng)調(diào)整間斷節(jié)點(diǎn)的分布，保證了激波陣面在演化中不失真；局部自動(dòng)網(wǎng)格重構(gòu)策略的實(shí)現(xiàn)大大減少了人為干預(yù)，顯著提高了算法的計(jì)算效率；此外，結(jié)合數(shù)值算例，提出了一種針對(duì)激波干擾點(diǎn)運(yùn)動(dòng)的裝配策略。仿真結(jié)果表明了算法的合理性和有效性。

1 自適應(yīng)間斷裝配求解器

本節(jié)首先對(duì)激波捕捉算法和以此為基礎(chǔ)發(fā)展而來的ADFs進(jìn)行簡單介紹，然后重點(diǎn)說明針對(duì)非定常問題的一些改進(jìn)。

1.1 激波捕捉算法

考慮網(wǎng)格運(yùn)動(dòng)變形，控制方程采用任意拉格朗日-歐拉(Arbitrary Lagrange-Euler，ALE)法描述的二維可壓縮非定常Euler方程的積分形式：

(1)

式中：Q為守恒變量；Fc為對(duì)流項(xiàng)通量；xc為網(wǎng)格運(yùn)動(dòng)速度；n為控制體表面外法向單位向量；V和S分別表示控制體和控制體表面；dσ和dψ分別表示體積微元和面積微元；ρ為密度；ui為速度；E為單位質(zhì)量流體總能；p為壓強(qiáng)；δij為Kronecker函數(shù)；考慮量熱完全氣體，比熱比γ=1.4。

數(shù)值計(jì)算采用格心型有限體積法，流場變量僅存儲(chǔ)在網(wǎng)格單元中心。不特別說明，本文均采用AUSM+(Advection Upstream Splitting Method+)格式計(jì)算面元對(duì)流通量，采用顯式4步Runge-Kutta格式進(jìn)行時(shí)間推進(jìn)，以此保證在光滑流場中時(shí)空離散上均具有二階精度。

非結(jié)構(gòu)動(dòng)網(wǎng)格技術(shù)方面，本文利用改進(jìn)后的彈簧近似法[18]，避免網(wǎng)格在變形運(yùn)動(dòng)中出現(xiàn)失效單元，保持較好的網(wǎng)格質(zhì)量，且不會(huì)造成求解器計(jì)算效率的大幅度降低。此外，當(dāng)流場邊界因出現(xiàn)大變形、大位移而導(dǎo)致網(wǎng)格質(zhì)量急劇下降時(shí)，往往需要進(jìn)行網(wǎng)格重構(gòu)。這又涉及到新/舊網(wǎng)格間的流場信息傳遞，為了盡量減小插值引起的空間精度下降和非物理振蕩，本文采用基于單元相交的混合網(wǎng)格精確守恒插值方法[19]。同時(shí)，為了防止因網(wǎng)格運(yùn)動(dòng)引入非物理解，還需要根據(jù)“離散幾何守恒率”對(duì)計(jì)算方法進(jìn)行修正。由于篇幅所限，相關(guān)算法這里不再詳述，可參考文獻(xiàn)[18]。

1.2 自適應(yīng)間斷裝配算法

1.2.1 標(biāo)記間斷位置

在使用激波裝配法時(shí)，首先都需要給定擬裝配的初始間斷位置。一般來說，對(duì)于定常流動(dòng)，只需給定大概位置，不要求十分準(zhǔn)確，隨著流場收斂間斷最終便會(huì)移動(dòng)到準(zhǔn)確的位置[17]；然而，對(duì)于間斷處于不斷運(yùn)動(dòng)的非定常流動(dòng)，初始間斷位置的標(biāo)記比較重要，若偏差過大可能會(huì)嚴(yán)重影響后續(xù)流場的演化。

圖1 自適應(yīng)間斷求解器的關(guān)鍵步驟

如圖1(a)所示，ADFs通過指定網(wǎng)格節(jié)點(diǎn)的屬性(間斷或普通)來標(biāo)記間斷的位置，間斷陣面由一系列間斷面元組成，而間斷面元同時(shí)也作為流場網(wǎng)格單元的邊界；每個(gè)間斷陣面上的間斷節(jié)點(diǎn)均存儲(chǔ)兩組流場數(shù)據(jù)，分別對(duì)應(yīng)間斷兩側(cè)的狀態(tài)。特別地，對(duì)于由m條間斷匯聚的干擾點(diǎn)，則該點(diǎn)將存儲(chǔ)2m組流場參數(shù)，分別對(duì)應(yīng)m條間斷的兩側(cè)。

1.2.2 初始化間斷節(jié)點(diǎn)參數(shù)

對(duì)于格心型有限體積法，已知t時(shí)刻的流場數(shù)據(jù)僅存儲(chǔ)于所有網(wǎng)格單元中心?？紤]到ADFs的裝配計(jì)算是建立在網(wǎng)格節(jié)點(diǎn)上的，因此采用反距離加權(quán)(Inverse Distance Weighted)法，如圖1(b)所示，由間斷兩側(cè)貢獻(xiàn)單元的格心值(Uu和Ud)插值得到所有間斷節(jié)點(diǎn)兩側(cè)參數(shù)(Vu和Vd)，即

(2)

此外，間斷節(jié)點(diǎn)j的法向構(gòu)造如下：

(3)

1.2.3 計(jì)算并修正間斷節(jié)點(diǎn)參數(shù)

1.2.2節(jié)得到的間斷節(jié)點(diǎn)兩側(cè)參數(shù)并不滿足描述間斷的跳躍關(guān)系式，如Rankine-Hugoniot(R-H)關(guān)系式。對(duì)激波來講，在激波坐標(biāo)系下，由于上游流動(dòng)是超聲速的，即下游流動(dòng)信息理論上不會(huì)污染上游參數(shù)，因此上游流動(dòng)參數(shù)Vu無需修正。沿著間斷節(jié)點(diǎn)法向n，一共有4個(gè)參數(shù)需要計(jì)算，即修正后的激波下游流動(dòng)參數(shù)V′d=[ρ′d,u′d,p′d]和激波節(jié)點(diǎn)的運(yùn)動(dòng)速度ω=ω·n。R-H關(guān)系式和下游特征相容關(guān)系式分別為

(4)

(5)

聯(lián)立式(4)和式(5)即可求得V′d和ω。

此外，在求得激波節(jié)點(diǎn)運(yùn)動(dòng)速度后，便可以計(jì)算出間斷節(jié)點(diǎn)在激波坐標(biāo)系下的上游相對(duì)馬赫數(shù)為

(6)

1.2.4 網(wǎng)格節(jié)點(diǎn)運(yùn)動(dòng)

(7)

如圖1(d)所示，通過彈簧近似法便可求得流場其余所有普通節(jié)點(diǎn)在新時(shí)刻的位置，從而更新得到新的流場網(wǎng)格。值得注意的是，對(duì)于定常流動(dòng)，最終ω將趨于0，即間斷位置收斂，流場網(wǎng)格將不發(fā)生變形。

1.2.5 流場更新

基于新的流場網(wǎng)格，采用1.1節(jié)所述的激波捕捉求解器進(jìn)行流場更新，獲得t+Δt時(shí)刻各個(gè)網(wǎng)格單元的格心值，為此需要計(jì)算各個(gè)面元的對(duì)流通量。對(duì)于普通面元，仍可以采用各種空間離散方法進(jìn)行通量分解計(jì)算，無需任何修正；而對(duì)于間斷面元，其通量可以根據(jù)上游參數(shù)直接寫出：

(8)

最后，通過時(shí)間推進(jìn)算法便可得到t+Δt時(shí)刻的流場，重復(fù)1.2.1節(jié)～1.2.5節(jié)的步驟進(jìn)行下一時(shí)間步的裝配計(jì)算直到計(jì)算終止。

1.3 非定常激波裝配的一些改進(jìn)

非定常激波相比定常激波的一個(gè)顯著特點(diǎn)是其陣面的位置或形狀會(huì)隨時(shí)間發(fā)生變化，因此算法需要具備有效模擬激波形狀或長度發(fā)生劇烈變化的能力。本節(jié)針對(duì)運(yùn)動(dòng)激波傳播特性，對(duì)ADFs進(jìn)行相應(yīng)的改進(jìn)，使其能準(zhǔn)確高效地模擬激波運(yùn)動(dòng)問題。

1.3.1 間斷節(jié)點(diǎn)沿壁面運(yùn)動(dòng)

采用ADFs模擬管道中的運(yùn)動(dòng)激波時(shí)，需要對(duì)沿壁面運(yùn)動(dòng)的間斷節(jié)點(diǎn)進(jìn)行特殊處理。根據(jù)壁面的彎曲特性，一般可以將其分為直壁面和曲壁面兩大類，對(duì)兩者的處理有明顯區(qū)別。

對(duì)于直壁面，如圖1(d)所示，間斷節(jié)點(diǎn)可以沿壁面方向自由滑動(dòng)，不需要跨越任何壁面網(wǎng)格節(jié)點(diǎn)；而間斷附近的壁面節(jié)點(diǎn)根據(jù)非結(jié)構(gòu)動(dòng)網(wǎng)格技術(shù)調(diào)整位置，同樣沿壁面滑動(dòng)，不會(huì)破壞直壁面的形狀。

對(duì)于曲壁面，處理要復(fù)雜一點(diǎn)，由于其是由一系列非共線的小直線段(壁面網(wǎng)格面元)離散的，為了保證貼體特性，曲壁面上的網(wǎng)格節(jié)點(diǎn)在計(jì)算中通常固定不動(dòng)，因此間斷節(jié)點(diǎn)必然會(huì)跨越壁面網(wǎng)格節(jié)點(diǎn)。針對(duì)彎曲方向的不同，曲壁面又可分為凸曲壁和凹曲壁，兩者的處理方式略有差異，下面簡要介紹。

如圖2所示，用i、j、k、m、n表示各個(gè)網(wǎng)格節(jié)點(diǎn)。在t時(shí)刻，間斷節(jié)點(diǎn)i沿面元mk向壁面節(jié)點(diǎn)k滑動(dòng)，根據(jù)間斷節(jié)點(diǎn)沿壁面運(yùn)動(dòng)速度ωt,i，若計(jì)算有ωt,iΔt>lik，其中l(wèi)ik為節(jié)點(diǎn)i、k之間的距離，則說明下一時(shí)刻間斷節(jié)點(diǎn)i將跨越節(jié)點(diǎn)k。

圖2 間斷節(jié)點(diǎn)沿曲壁面的運(yùn)動(dòng)

在確定t+Δt時(shí)刻新的壁面間斷節(jié)點(diǎn)位置時(shí)，首先根據(jù)式(7)計(jì)算出間斷節(jié)點(diǎn)i的虛擬推進(jìn)位置k′；對(duì)于凹曲壁(見圖2(a2))，新間斷節(jié)點(diǎn)取直線jk′與面元in的交點(diǎn)；對(duì)于凸曲壁(見圖2(b2))，過點(diǎn)k′作面元in的垂線，新間斷節(jié)點(diǎn)取垂足。

實(shí)際上，從動(dòng)網(wǎng)格角度考慮，該過程相當(dāng)于間斷節(jié)點(diǎn)i移動(dòng)到了壁面網(wǎng)格節(jié)點(diǎn)處并退化成普通節(jié)點(diǎn)，而對(duì)應(yīng)壁面網(wǎng)格節(jié)點(diǎn)k移動(dòng)到新的壁面間斷位置并轉(zhuǎn)化成間斷節(jié)點(diǎn)。從間斷角度考慮，相當(dāng)于間斷節(jié)點(diǎn)跨越了運(yùn)動(dòng)方向上鄰近的曲壁面網(wǎng)格節(jié)點(diǎn)。值得注意的是，在跨越前后，由節(jié)點(diǎn)i、j、k組成的三角形網(wǎng)格單元①從激波上游移動(dòng)到了激波下游，因此需要對(duì)其格心參數(shù)重新進(jìn)行賦值，直接由間斷節(jié)點(diǎn)下游激波參數(shù)插值即可得到。

1.3.2 間斷節(jié)點(diǎn)分布的自動(dòng)重構(gòu)

對(duì)于定常間斷裝配，由于間斷移動(dòng)變形有限，通常間斷節(jié)點(diǎn)能維持較好的分布；在前期的一些裝配計(jì)算中，若有必要調(diào)整間斷節(jié)點(diǎn)分布時(shí)，通常要人為重新劃分網(wǎng)格[17]，費(fèi)時(shí)費(fèi)力。而間斷在移動(dòng)過程中往往會(huì)不斷變長或縮短，從而造成相鄰間斷節(jié)點(diǎn)的間距過大或過小，使得間斷陣面附近的網(wǎng)格質(zhì)量迅速下降，如圖3所示。當(dāng)這種情況發(fā)生時(shí)，需要沿間斷陣面重新調(diào)整間斷節(jié)點(diǎn)分布，使其間距近似等于周邊網(wǎng)格單元的尺寸。對(duì)于二維問題，具體步驟如下：

圖3 間斷節(jié)點(diǎn)分布重構(gòu)前后的網(wǎng)格

步驟1判斷間斷節(jié)點(diǎn)是否滿足分布重構(gòu)條件關(guān)系式：

lk/lave<δmin或lk/lave>δmax

(9)

式中:lk為第k個(gè)舊間斷面元的長度；lave為間斷鄰近網(wǎng)格單元的平均尺度；δmin和δmax分別為預(yù)設(shè)的最小間距和最大間距系數(shù)，通常取0.3和1.7。

步驟2對(duì)需要分布重構(gòu)的各組間斷散點(diǎn)進(jìn)行排序并曲線擬合?？紤]到每條間斷陣面的始點(diǎn)和終點(diǎn)在重構(gòu)前后的位置固定，本文采用多次Bézier曲線擬合法[20]得到各條間斷擬合曲線。

步驟3將擬合曲線按照尺寸βlave進(jìn)行平均分割獲得各個(gè)均分點(diǎn)。為了更好地描述曲線的彎曲程度，用β來控制均分點(diǎn)的疏密程度，對(duì)較為平直的曲線，β通常取1.2；而對(duì)于比較彎曲的曲線，β通常取0.7。

步驟4將均分點(diǎn)映射到舊間斷面元上從而得到分布重構(gòu)后的間斷節(jié)點(diǎn)。注意，如果不進(jìn)行映射，直接采用均分點(diǎn)作為新間斷節(jié)點(diǎn)，會(huì)在一定程度上損失位置精度，當(dāng)進(jìn)行多次分布重構(gòu)后位置誤差積累可能造成間斷位置形狀出現(xiàn)顯著偏差。

步驟5基于新的間斷節(jié)點(diǎn)，重新劃分間斷附近網(wǎng)格，通過插值算法獲得新網(wǎng)格單元格心的流場參數(shù)。

該方法不僅用于間斷節(jié)點(diǎn)的分布重構(gòu)，同樣，當(dāng)間斷節(jié)點(diǎn)在直壁面上滑動(dòng)時(shí)，若干壁面網(wǎng)格節(jié)點(diǎn)難免會(huì)過近或過遠(yuǎn)，此時(shí)便需進(jìn)行壁面網(wǎng)格節(jié)點(diǎn)的分布重構(gòu)。

1.3.3 流場網(wǎng)格的局部自動(dòng)重構(gòu)

對(duì)于定常激波裝配，由于間斷陣面的變形幅度較小，一般只需進(jìn)行幾次網(wǎng)格重構(gòu)即可，可以輸出網(wǎng)格并導(dǎo)入專業(yè)CFD網(wǎng)格生成軟件(如Pointwise)人工生成網(wǎng)格，然后再讀入網(wǎng)格完成重構(gòu)過程。而對(duì)于運(yùn)動(dòng)激波裝配，由于網(wǎng)格變形比較劇烈，需要頻繁重構(gòu)網(wǎng)格，若再進(jìn)行人工手動(dòng)重構(gòu)將付出極大的計(jì)算代價(jià)，嚴(yán)重影響計(jì)算效率，因此需要高效的自動(dòng)網(wǎng)格生成策略。

本文成功將Shewchuk研發(fā)并開源的二維Delaunay三角形網(wǎng)格生成程序“Triangle”[21]嵌入到ADFs中，實(shí)現(xiàn)了流場網(wǎng)格的自動(dòng)重構(gòu)。下面以一道右行激波在平直管道中運(yùn)動(dòng)為例，通過展示若干關(guān)鍵時(shí)刻的流場網(wǎng)格，來具體說明該重構(gòu)策略的特性：

1) 采用局部網(wǎng)格變形/重構(gòu)策略能顯著提高計(jì)算效率，同時(shí)重構(gòu)單元數(shù)量的縮減很大程度上可以減小流場信息插值引入的誤差。如圖4(a)所示，以間斷位置作為起點(diǎn)向兩側(cè)搜尋網(wǎng)格，標(biāo)記距離間斷最近的n層網(wǎng)格單元作為可變形單元且可以重構(gòu)，其他網(wǎng)格靜止且不參與重構(gòu)，即流場可分為“動(dòng)網(wǎng)格區(qū)”(圖4陰影區(qū))和“靜網(wǎng)格區(qū)”。當(dāng)激波在直壁面上運(yùn)動(dòng)時(shí)，n通常取10；而在曲壁面上運(yùn)動(dòng)時(shí)，由于網(wǎng)格重構(gòu)較為頻繁，動(dòng)網(wǎng)格區(qū)可以取小一點(diǎn)，如n=3。

2) 當(dāng)激波運(yùn)動(dòng)使得網(wǎng)格質(zhì)量較差時(shí)(如圖4(b))，需要進(jìn)行網(wǎng)格重構(gòu)?？紤]到計(jì)算初始網(wǎng)格通常用專業(yè)網(wǎng)格生成軟件劃分，網(wǎng)格質(zhì)量可以控制得比較好(如大部分為正三角形)，因此ADFs在網(wǎng)格重構(gòu)時(shí)利用了初始網(wǎng)格信息。如圖4(c)所示，提取位于重構(gòu)區(qū)內(nèi)的初始網(wǎng)格點(diǎn)，作為網(wǎng)格生成“控制點(diǎn)”供Triangle程序調(diào)用生成Delaunay三角形網(wǎng)格。注意，為了保證網(wǎng)格質(zhì)量，距離間斷過近(約一個(gè)網(wǎng)格尺寸內(nèi))的初始網(wǎng)格點(diǎn)不作為控制點(diǎn)。這樣，經(jīng)控制生成的大部分網(wǎng)格單元的尺寸和形狀與初始高質(zhì)量網(wǎng)格保持一致(對(duì)比圖4(a)與圖4(c))，一定程度上減弱了網(wǎng)格重構(gòu)對(duì)流場網(wǎng)格分布和拓?fù)涞挠绊憽?/p>

圖4 間斷附近局部網(wǎng)格自動(dòng)重構(gòu)(n=5)

3) 每次完成網(wǎng)格重構(gòu)后，由于舊網(wǎng)格變形區(qū)已經(jīng)不太適合新的間斷位置，因此需要重新標(biāo)記動(dòng)網(wǎng)格區(qū)和靜網(wǎng)格區(qū)，如圖4(d)所示，即計(jì)算過程中流場網(wǎng)格重構(gòu)區(qū)是動(dòng)態(tài)變化的，隨間斷移動(dòng)需要?jiǎng)討B(tài)調(diào)整。

2 驗(yàn)證算例

本節(jié)通過3個(gè)二維算例對(duì)若干改進(jìn)的準(zhǔn)確性、可行性和可靠性做了驗(yàn)證，體現(xiàn)了該非定常激波裝配算法在處理運(yùn)動(dòng)激波的潛力。值得強(qiáng)調(diào)的是，以下所有裝配算例整個(gè)計(jì)算完全自動(dòng)進(jìn)行，沒有經(jīng)過任何人為干預(yù)。

2.1 激波繞射90°拐角

平面運(yùn)動(dòng)激波繞射尖銳拐角是一類相當(dāng)常見的氣動(dòng)問題，許多學(xué)者從試驗(yàn)和數(shù)值的角度出發(fā)，對(duì)90°拐角的激波繞射問題做了大量研究[22-24]。根據(jù)入射運(yùn)動(dòng)激波的強(qiáng)弱，該問題主要分為兩類，即激波下游的流場是超聲速或亞聲速。圖5給出了強(qiáng)激波繞射下該流場的若干主要波系結(jié)構(gòu)，流場在演化中具有明顯的自相似性和準(zhǔn)定常特性。

圖5 強(qiáng)激波繞射90°拐角流場結(jié)構(gòu)[22]

本節(jié)考慮一道平直入射激波(運(yùn)動(dòng)激波馬赫數(shù)Mas=2.4)繞射90°拐角的情況，初始流場參數(shù)詳見表1，流場計(jì)算區(qū)域如圖6所示，初始激波位于拐角前0.5L處，網(wǎng)格尺度Δ=0.25L，全場初始共有28 762個(gè)均勻分布的三角形網(wǎng)格單元。本算例中的時(shí)間均為無量綱量，計(jì)算中只對(duì)入射強(qiáng)激波進(jìn)行裝配，流場演化產(chǎn)生的新激波仍然用捕捉法計(jì)算。

表1 激波繞射初始流場條件

圖6直觀地給出了一系列裝配得到的運(yùn)動(dòng)激波陣面，無量綱間隔時(shí)間dt=0.05。可以看出當(dāng)入射強(qiáng)激波繞過壁面拐角后，激波會(huì)逐漸向垂直壁面彎曲，這是由于拐角附近的流場在膨脹波干擾下，激波強(qiáng)度減弱，激波運(yùn)動(dòng)速度降低；此外，通過觀察可以發(fā)現(xiàn)，相比流場上方未受干擾的入射激波，彎向垂直壁面的激波間距明顯減小，經(jīng)推算此處的激波馬赫數(shù)約降為1.36。

圖6 不同時(shí)刻下繞射激波的裝配陣面

由于該算例的激波陣面在傳播過程發(fā)生了劇烈地彎曲和伸長，因此需要頻繁進(jìn)行間斷節(jié)點(diǎn)的分布重構(gòu)。從某一時(shí)刻流場出發(fā)，考察了一段時(shí)間后是否進(jìn)行間斷節(jié)點(diǎn)分布重構(gòu)的流場網(wǎng)格，結(jié)果對(duì)比如圖7所示。其中，左圖未進(jìn)行間斷節(jié)點(diǎn)重構(gòu)，且采用傳統(tǒng)的全局網(wǎng)格變形策略，可以發(fā)現(xiàn)局部間斷面元被過度拉長或擠短，不能較好地描述間斷曲率；同時(shí)激波下游網(wǎng)格尺度將逐漸增大，會(huì)降低此處復(fù)雜流場的解析精度。相反，算法改進(jìn)后，經(jīng)過若干次網(wǎng)格重構(gòu)，間斷陣面及其附近的網(wǎng)格尺度和鄰近網(wǎng)格始終保證接近，且間斷的形狀更加光滑合理，充分反映了1.3節(jié)改進(jìn)策略的必要性和有效性。

為了定性說明裝配結(jié)果的準(zhǔn)確性，基于相同網(wǎng)格尺度用激波捕捉法重新計(jì)算了該問題。圖8給出了2個(gè)無量綱時(shí)刻流場密度等值線的對(duì)比(ρ0為初始入射激波上游密度)，由于捕捉激波本質(zhì)上存在的數(shù)值誤差，捕捉法計(jì)算的密度等值線(白色實(shí)線)會(huì)出現(xiàn)兩條不合理的波帶，嚴(yán)重降低了流場等值線的光滑性；而裝配計(jì)算的等值線(黑色實(shí)線)光滑合理，沒有出現(xiàn)誤差帶，整體看來裝配法的流場結(jié)果更加準(zhǔn)確可信。

圖8 不同時(shí)刻流場密度云圖對(duì)比

圖9對(duì)比了裝配法和捕捉法計(jì)算終止時(shí)刻流場壓強(qiáng)云圖(p0為初始入射激波上游壓強(qiáng))，可以看到兩者的激波位置和大部分區(qū)域的等值線符合得較好。進(jìn)一步提取x=1.6L直線上的數(shù)據(jù)，定量比較兩者壓強(qiáng)分布，如圖10所示，裝配法避免了捕捉法處理激波產(chǎn)生的虛假數(shù)值過渡區(qū)，而在膨脹波區(qū)兩者吻合很好。以上說明了對(duì)于非定常流場，裝配法即使經(jīng)歷了多次網(wǎng)格重構(gòu)和流場信息插值，仍能保持較高的計(jì)算精度。

圖9 t=0.8時(shí)刻流場壓強(qiáng)云圖對(duì)比

圖10 x=1.6L處壓強(qiáng)分布對(duì)比

2.2 二維激波管內(nèi)激波增強(qiáng)

本節(jié)模擬二維激波管內(nèi)激波增強(qiáng)問題。計(jì)算模型如圖11所示，通過設(shè)計(jì)特定的上下壁面收縮型線[25]，初始平面運(yùn)動(dòng)激波依次經(jīng)過光滑的凹形曲線段、斜直線段和光滑的凸形曲線段，先后受到壁面產(chǎn)生的“激波-壓縮”和“激波-膨脹”擾動(dòng)，形狀不斷發(fā)生改變，最終在收縮段出口恢復(fù)為增強(qiáng)的平面激波，且激波波面上沒有明顯擾動(dòng)。

圖11 收縮壁面型線設(shè)計(jì)示意圖[25]

該激波管收縮段入口高度H0=70 mm，出口高度H1=8 mm，匯集角θ=20°；收縮段總長179.4 mm，其中凹曲壁、斜直壁和凸曲壁的水平長度分別為139.8、18.1、21.5 mm。仿真計(jì)算時(shí)，初始入射激波(Mas=3.2)位于收縮段入口前方8 mm處，激波上游為靜止空氣，且上游壓強(qiáng)p0和密度ρ0分別為6 kPa和0.072 855 kg/m3。采用激波裝配法計(jì)算時(shí)，全場網(wǎng)格尺寸ΔLSF=1 mm，一共有25 082個(gè)均勻分布的三角形網(wǎng)格單元。

在整個(gè)過程中，激波一共運(yùn)動(dòng)了170 μs。圖12依次給出了間隔10 μs裝配的激波陣面。可以清晰直觀地看到，平面激波在進(jìn)入收縮段后開始逐漸自底部向中心彎曲，在斜直壁段時(shí)曲率保持不變，當(dāng)傳播到凸曲壁時(shí)，激波陣面曲率逐漸減小，最終在出口位置再次轉(zhuǎn)變?yōu)槠矫婕げā氖湛s段入口運(yùn)動(dòng)到出口，激波長度明顯縮小，采用1.3.2節(jié)所述的間斷節(jié)點(diǎn)自動(dòng)分布重構(gòu)技術(shù)，激波陣面上的間斷節(jié)點(diǎn)從初始71個(gè)逐漸減少到12個(gè)。根據(jù)160～170 μs時(shí)段中激波的水平運(yùn)動(dòng)16.4 mm，近似求出收縮段出口位置的激波馬赫數(shù)為4.83，與試驗(yàn)紋影[25]推算出的激波馬赫數(shù)4.80吻合得非常好，相對(duì)偏差僅有0.6%。

圖12 激波管內(nèi)不同時(shí)刻裝配激波陣面

文獻(xiàn)[26]理論分析了激波運(yùn)動(dòng)到不同水平位置的強(qiáng)度，即運(yùn)動(dòng)激波馬赫數(shù)Mas，考慮到裝配法需要計(jì)算出間斷節(jié)點(diǎn)的運(yùn)動(dòng)速度，進(jìn)而很方便地推導(dǎo)出激波運(yùn)動(dòng)馬赫數(shù)，因此本文就該數(shù)值進(jìn)行對(duì)比，結(jié)果如圖13所示?？梢钥闯觯诎记诙?，裝配結(jié)果和理論預(yù)測(cè)符合得很好，而在斜直壁段和凸曲壁段，二者略有差別，有待進(jìn)一步校核；中軸線附近的激波強(qiáng)度在x<80 mm段基本沒變化，受到壓縮波作用之后迅速增強(qiáng)，且在斜直壁段與理論預(yù)測(cè)出的壁面激波強(qiáng)度較為吻合，同樣，在靠近出口處，受膨脹波的影響強(qiáng)度有所下降?？傊?，采用激波裝配法可以很容易得到激波在傳播過程中的強(qiáng)度變化，而激波捕捉法很難直接得到，這也是裝配法的一個(gè)優(yōu)勢(shì)所在。

圖14依次對(duì)比了激波捕捉法和裝配法模擬的密度ρ、壓強(qiáng)p和熵S=p/ργ的云圖，圖中S0為初始入射激波上游熵值。捕捉法的流場網(wǎng)格尺度ΔLSC=ΔLSF/5，全場共有626 128個(gè)三角形網(wǎng)格單元，即網(wǎng)格單元總量約是裝配法的25倍。對(duì)比密度和壓強(qiáng)云圖可以發(fā)現(xiàn)，捕捉法得到的等值線在x=30，130 mm左右存在明顯的振蕩，而裝配法得到的等值線光滑合理沒有波動(dòng)；從熵云圖可以大致看出造成捕捉法出現(xiàn)非物理波動(dòng)的原因，捕捉得到的激波下游會(huì)出現(xiàn)一條逐漸遠(yuǎn)離激波的熵帶，一定程度污染了下游流場，這可能是由捕捉激波從本質(zhì)上產(chǎn)生的數(shù)值誤差造成的，文獻(xiàn)[12,15]在計(jì)算靜止激波與渦的相互作用中也同樣發(fā)現(xiàn)這種現(xiàn)象，此處有待進(jìn)一步探討。相比之下，裝配法由于精確處理了激波，即使用很少的網(wǎng)格也能得到較好的流場。

圖13 激波在不同水平位置的強(qiáng)度對(duì)比

圖14 t=165 μs時(shí)刻裝配法和捕捉法模擬的流場對(duì)比

2.3 彎管內(nèi)激波傳播

圖15 激波在彎管中傳播的試驗(yàn)陰影圖[27]

爆震發(fā)動(dòng)機(jī)中經(jīng)常會(huì)出現(xiàn)激波在彎曲管道中的傳播現(xiàn)象，其典型波系結(jié)構(gòu)如圖15試驗(yàn)陰影圖[27]所示。本節(jié)采用激波裝配法模擬激波在“L”形彎管內(nèi)的傳播、加速過程，進(jìn)一步考核本方法的可靠性。

計(jì)算模型由直管段和“L”形彎管段組成，如圖16所示，內(nèi)徑均為L，且彎管段的內(nèi)壁面曲率半徑r=0.5L。初始入射激波馬赫數(shù)Mas=4.0，下游與上游的壓強(qiáng)比和密度比分別為18.5、4.571；且初始入射激波上游為靜止空氣，下游氣流馬赫數(shù)為1.553。全場網(wǎng)格尺度Δ=0.03L，初始共有11 654個(gè)均勻分布的三角形網(wǎng)格單元。本算例中的時(shí)間均為無量綱量。

圖16 計(jì)算區(qū)域、初始條件和網(wǎng)格

圖17 不同時(shí)刻下L形彎管內(nèi)裝配的激波陣面

圖17給出了間隔dt=0.04裝配的激波陣面，入射激波于t=0.22時(shí)刻到達(dá)彎曲段入口，此時(shí)激波陣面仍是平整的，激波進(jìn)入彎曲管道后同時(shí)受到內(nèi)壁面的膨脹干擾及外壁面的壓縮干擾，使得靠近內(nèi)壁面的激波逐漸減弱并發(fā)生彎曲以保持傳播方向與管道軸線一致，而靠近外壁面的激波受到壓縮波的作用而逐漸增強(qiáng)最終形成馬赫反射結(jié)構(gòu)。圖中可以清晰地看到馬赫桿隨著激波的傳播逐漸變長，即三波點(diǎn)逐漸向內(nèi)壁面移動(dòng)；且從間距可以看出馬赫桿的強(qiáng)度更高，運(yùn)動(dòng)速度更快。

非定常流動(dòng)中三波點(diǎn)等間斷相交點(diǎn)運(yùn)動(dòng)速度的求解一直以來是個(gè)難題，目前尚未有完善的理論模型，因此運(yùn)動(dòng)相交點(diǎn)的處理是裝配法的一個(gè)研究重點(diǎn)。圖18給出了某時(shí)刻間斷節(jié)點(diǎn)運(yùn)動(dòng)示意，其中三波點(diǎn)的運(yùn)動(dòng)速度沒有修正，造成鄰近的間斷面元斜率出現(xiàn)異常，隨著時(shí)間推進(jìn)，三波點(diǎn)位置會(huì)迅速發(fā)散引起計(jì)算失敗。

本文提出一種“預(yù)估位移，推算速度”的策略間接完成運(yùn)動(dòng)三波點(diǎn)的裝配，下面結(jié)合圖19介紹其關(guān)鍵步驟：

圖18 未修正的三波點(diǎn)運(yùn)動(dòng)

圖19 三波點(diǎn)運(yùn)動(dòng)速度計(jì)算模型

步驟1辨識(shí)三波點(diǎn)(TP)。分析可知三波點(diǎn)一般位于幾個(gè)不同斜率激波線的交匯處，如圖19所示，求解每個(gè)間斷節(jié)點(diǎn)左右兩個(gè)相鄰間斷線元的夾角ψ，若ψ高于某閾值ψthr，則標(biāo)記該節(jié)點(diǎn)為三波點(diǎn)。數(shù)值試驗(yàn)表明ψthr通常取15°即可。

步驟2求解相交點(diǎn)新時(shí)刻位置。首先根據(jù)1.2.3節(jié)求得的間斷節(jié)點(diǎn)運(yùn)動(dòng)速度，分別計(jì)算三波點(diǎn)TP左鄰居間斷節(jié)點(diǎn)(P1～P3)和右鄰居間斷節(jié)點(diǎn)(P4～P6)在t+Δt時(shí)刻的位置(圖19藍(lán)色實(shí)心方塊)，然后分別擬合出兩條直線(圖19黑色實(shí)線)，可以認(rèn)為這兩條直線與新時(shí)刻三波點(diǎn)附近的激波線斜率近似一致，最后求解這兩條直線的相交點(diǎn)(圖19紅色實(shí)心方塊)作為新時(shí)刻三波點(diǎn)的位置。

步驟3計(jì)算三波點(diǎn)速度。由舊三波點(diǎn)到新三波點(diǎn)的位移矢量為LTP，則三波點(diǎn)的運(yùn)動(dòng)速度矢量ωTP=LTP/Δt。

值得注意的是，雖然此處針對(duì)的是三波點(diǎn)結(jié)構(gòu)，但其他間斷相交點(diǎn)的運(yùn)動(dòng)也可以類似處理，這里不再詳述。

為了校核本算例三波點(diǎn)位置計(jì)算的準(zhǔn)確性，圖20對(duì)比了任意兩個(gè)時(shí)刻捕捉法和裝配法的結(jié)果，可見裝配的入射激波、馬赫桿和三波點(diǎn)位置和捕捉結(jié)果基本吻合。此外，對(duì)比圖20和圖16的流場網(wǎng)格，可見裝配激波走過的區(qū)域，經(jīng)1.3.3節(jié)所述的網(wǎng)格重構(gòu)策略，流場網(wǎng)格可以恢復(fù)到初始高質(zhì)量網(wǎng)格，很好地減弱了裝配陣面對(duì)流場網(wǎng)格的影響。

圖20 不同時(shí)刻流場密度云圖及網(wǎng)格對(duì)比

3 結(jié) 論

本文面向非定常流動(dòng)中激波傳播問題，對(duì)最近發(fā)展的自適應(yīng)間斷裝配求解器ADFs進(jìn)行了若干改進(jìn)，主要結(jié)論如下：

1) 當(dāng)激波在曲壁面上運(yùn)動(dòng)時(shí)，運(yùn)動(dòng)間斷節(jié)點(diǎn)可以在固定壁面網(wǎng)格節(jié)點(diǎn)的前提下實(shí)現(xiàn)裝配。

2) 針對(duì)激波產(chǎn)生大變形、大位移運(yùn)動(dòng)問題，通過間斷節(jié)點(diǎn)分布的自動(dòng)重構(gòu)保證激波陣面不失真，同時(shí)采用局部網(wǎng)格自動(dòng)重構(gòu)策略確保了網(wǎng)格的高質(zhì)量，并提高了程序計(jì)算效率。

3) 對(duì)于激波相交點(diǎn)的運(yùn)動(dòng)，設(shè)計(jì)了一種根據(jù)位移推算速度的方法進(jìn)行裝配，數(shù)值算例表明該方法行之有效。

總之，采用激波裝配法處理激波傳播問題，相比激波捕捉方法可以獲得流場間斷更加直觀清晰的圖譜，有利于進(jìn)一步深入認(rèn)識(shí)激波的演化過程。本文僅對(duì)單個(gè)運(yùn)動(dòng)激波進(jìn)行裝配，流場演化中的新生激波仍用捕捉法進(jìn)行處理，而激波裝配法的目標(biāo)是顯式地自動(dòng)模擬這些間斷拓?fù)渥兓瘑栴}，因此將激波辨識(shí)技術(shù)融合到激波裝配法中將是下一步工作的重點(diǎn)。

致 謝

中國科學(xué)技術(shù)大學(xué)的楊基明老師在激波增強(qiáng)管壁型線幾何模型方面給予了協(xié)助和寶貴的指導(dǎo)，在此致以衷心的感謝!