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

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

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

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

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

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

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

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

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

1.1 激波捕捉算法

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

(1)

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

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

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

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

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

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

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

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

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

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

(2)

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

(3)

1.2.3 計算并修正間斷節(jié)點參數(shù)

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

(4)

(5)

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

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

(6)

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

(7)

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

1.2.5 流場更新

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

(8)

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

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

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

1.3.1 間斷節(jié)點沿壁面運(yùn)動

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

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

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

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

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

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

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

1.3.2 間斷節(jié)點分布的自動重構(gòu)

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

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

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

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

(9)

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

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

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

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

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

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

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

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

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

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

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

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

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

2 驗證算例

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

2.1 激波繞射90°拐角

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

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

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

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

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

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

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

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

圖8 不同時刻流場密度云圖對比

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

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

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

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

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

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

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

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

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

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

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

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

圖14 t=165 μs時刻裝配法和捕捉法模擬的流場對比

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

圖15 激波在彎管中傳播的試驗陰影圖[27]

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

計算模型由直管段和“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個均勻分布的三角形網(wǎng)格單元。本算例中的時間均為無量綱量。

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

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

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

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

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

圖18 未修正的三波點運(yùn)動

圖19 三波點運(yùn)動速度計算模型

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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

致 謝

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