面向三維激波問題的裝配方法

鄒東陽，林敬周，黃潔，劉君

1. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所，綿陽 621000 2. 大連理工大學(xué) 航空航天學(xué)院，大連 116024

強(qiáng)擾動在空氣中以超過聲速的速度進(jìn)行傳播，表現(xiàn)為流體中出現(xiàn)一個參數(shù)的間斷面，稱為激波[1]。激波間斷的存在給描述無黏流動的Euler方程的求解帶來了很大挑戰(zhàn)。過去半個多世紀(jì)里，為了發(fā)展出能夠精確地描述激波的數(shù)值方法，科研人員做了大量的工作。如今這些工作是計算流體力學(xué)(Computational Fluid Dynamics, CFD)的重要組成部分。根據(jù)對激波處理方式的不同，CFD可以分為激波捕捉方法[2-8]和激波裝配方法[9-11]2種。激波捕捉方法是以弱解理論為基礎(chǔ)，全場采用統(tǒng)一的計算方法，無需特殊處理，具有通用性的特點。激波裝配方法則是利用激波關(guān)系式(Rankine-Hugoniot，R-H)對激波進(jìn)行描述，理論上是精確的，但是需要找到匹配激波間斷解與其他光滑區(qū)域數(shù)值解的方法，存在兼容性的問題。對于復(fù)雜激波流動問題，特別是三維問題，激波裝配方法存在應(yīng)用上的困難。隨著激波捕捉方法的成熟，特別是TVD(Total Variation Diminishing)格式的誕生，對于激波的計算逐步滿足了一定的要求。過去30多年里，激波捕捉方法得到了迅猛發(fā)展，在航空航天、汽車建筑、天文地理等領(lǐng)域發(fā)揮了巨大的作用。

雖然激波裝配方法的研究在相當(dāng)長一段時間內(nèi)寥若晨星，但是其在計算激波時表現(xiàn)出來的優(yōu)勢還是給許多激波捕捉方法研究提供了思路。像激波裝配方法那樣，將計算網(wǎng)格面與激波面相匹配來減少激波計算誤差的處理手段在很多方法中都得到了應(yīng)用。在20世紀(jì)90年代中期，Parpia和Parikh[12]、van Rosendale[13]以及Camarero[14]等發(fā)展了不同的方法，全都是將非結(jié)構(gòu)自適應(yīng)三角形網(wǎng)格的棱邊與流場中出現(xiàn)的間斷進(jìn)行匹配。類似的思想在新近發(fā)展的間斷伽遼金有限元方法(DG-FEM)中也有一定的體現(xiàn)。兩個完全獨(dú)立的研究團(tuán)隊同時關(guān)注于解決網(wǎng)格節(jié)點(不包括流場參數(shù))與間斷面相匹配的問題[15-20]。

另外一方面，在過去的十幾年里，CFD工作者對于激波裝配技術(shù)的興趣得到了急劇的提升，將目光從激波捕捉方法重新轉(zhuǎn)向激波裝配方法。造成這一現(xiàn)象的原因主要有2個方面：一是由于CFD工作者在解決激波捕捉方法本質(zhì)性局限所遇到的困境；二是由于一些研究人員將激波裝配與非結(jié)構(gòu)網(wǎng)格相結(jié)合取得了顯著成果，給解決激波捕捉方法所遇到的問題帶來了希望。

意大利學(xué)者Paciorri和Bonfiglioli自2006年開始著力于研究基于浸沒邊界法的激波裝配技術(shù)。他們方法有效地結(jié)合了浮動激波裝配和邊界激波裝配的特點，激波的運(yùn)動過程中允許其在背景網(wǎng)格上進(jìn)行(該特點類似浮動激波裝配)，計算過程中通過挖洞構(gòu)造以激波為內(nèi)邊界的新的計算網(wǎng)格進(jìn)行求解(該特點類似邊界激波裝配方法)。Paciorri和Bonfiglioli的方法較傳統(tǒng)的激波裝配方法在靈活性上有了顯著的提高，利用發(fā)展的裝配方法成功地解決了若干含有激波構(gòu)型的問題[21-23]。在國內(nèi)，劉君等基于非結(jié)構(gòu)動網(wǎng)格技術(shù)對激波裝配方法進(jìn)行了發(fā)展，提出了新型非結(jié)構(gòu)激波裝配方法—嵌入式激波裝配方法。嵌入式激波裝配方法中，激波屬于網(wǎng)格的一部分，激波節(jié)點的移動帶動其他網(wǎng)格節(jié)點的運(yùn)動，靈活性有了很大提高。該方法被成功用于求解許多定常/非定常流動問題[24-29]。

基于非結(jié)構(gòu)網(wǎng)格的裝配方法在諸多二維問題中的成功應(yīng)用展示了其相比傳統(tǒng)激波裝配方法在靈活性上的優(yōu)勢。相比二維問題，三維問題雖然只是增加了一個維度，裝配難度卻是有了大幅度提升。三維空間的激波辨識難度大，給初始激波位置的選擇帶來了困難。描述激波面的網(wǎng)格點在運(yùn)動過程中彼此要相互匹配，以保證激波收斂運(yùn)動過程不發(fā)散。在三維空間中控制這些點的匹配關(guān)系要比二維更為復(fù)雜，對算法要求也較高。此外，三維激波相交等問題缺乏理論指導(dǎo)，使得裝配方法在復(fù)雜問題中的應(yīng)用難度較大。諸多問題的存在使得針對三維問題的裝配方法應(yīng)用研究還較少，目前只是在文獻(xiàn)[30]中看到了相關(guān)報道。Bonfiglioli等[30]應(yīng)用其發(fā)展的激波裝配方法研究了一些超聲速流動問題，證實了其方法在解決三維問題中的有效性。盡管如此，從文獻(xiàn)給出的例子還是可以看出，裝配的激波相對較為簡單，流場內(nèi)激波都是相互獨(dú)立的，規(guī)避了三維激波相交帶來的問題。本文基于前期關(guān)于嵌入式激波裝配方法的研究工作，進(jìn)一步從二維拓展應(yīng)用到對三維問題的模擬中。與Bonfiglioli等[30]給出的三維裝配不同，本文在對三維激波問題進(jìn)行裝配時，不只是針對獨(dú)立激波，還考慮了更為復(fù)雜并且具有代表性的三維激波相交、反射等問題，結(jié)合給出的激波交點的計算方法，完成了對三類超聲速/高超聲速三維激波問題的裝配。

1 三維激波裝配方法

激波裝配方法由2個要素組成：計算和追蹤。所謂的“追蹤”指的是對流場中出現(xiàn)各類間斷位置進(jìn)行追蹤。所謂的“計算”指的是對間斷邊界包絡(luò)計算域內(nèi)流動利用CFD求解器對偏微分控制方程進(jìn)行求解的過程。

無論使用什么樣的CFD求解器，在非結(jié)構(gòu)網(wǎng)格范圍內(nèi)裝配激波是由一系列激波面所表示，這些激波面是由激波點自動連接而成。裝配間斷的位置是通過在激波點上求解R-H關(guān)系式獲得。

根據(jù)求解器是格心型或者格點型的不同，追蹤激波的方式也會有很大的差異，需要相應(yīng)地去修改CFD求解器。本文的激波裝配方使用的是格心型有限體積求解器。追蹤激波時，在激波點上求解R-H關(guān)系式，以此來獲得激波運(yùn)動速度。通過激波面元的通量則是使用上游參數(shù)以及激波速度得到。一旦獲得通過激波面元的通量，包含激波面元的單元參數(shù)的更新方式就和其他單元一樣了。裝配的激波按照求解的激波點速度在拉格朗日體系下移動。采用三維非結(jié)構(gòu)激波裝配方法進(jìn)行數(shù)值分析的過程涉及4個主要內(nèi)容，下面分別進(jìn)行詳細(xì)介紹。

1) 激波間斷標(biāo)識

激波裝配方法不像激波捕捉方法那樣對全場采用同樣的處理方式，而是將激波間斷單獨(dú)標(biāo)識出來，用一種解析的有別于光滑區(qū)域的方法進(jìn)行求解。傳統(tǒng)上的邊界激波裝配方法是將激波看作一種特殊的邊界進(jìn)行標(biāo)識和計算的。這種方法需要將計算區(qū)域劃分成若干個子區(qū)域，對于含有懸空激波以及拓?fù)浣Y(jié)構(gòu)變化的問題不易操作。本文中使用的激波裝配方法是通過節(jié)點屬性定義來對激波進(jìn)行標(biāo)識的。

如圖1所示，采用網(wǎng)格生成程序(Tetgen)或軟件(Pointwise)生成空間四面體網(wǎng)格。對離散后的網(wǎng)格節(jié)點進(jìn)行屬性定義，劃分為激波節(jié)點(Shock nodes, S)和普通網(wǎng)格節(jié)點(Common nodes, C)2種。激波面元是由激波節(jié)點連接構(gòu)成的，規(guī)定如果一個網(wǎng)格面元包含的節(jié)點全部為激波節(jié)點則該面元為激波面元。在計算過程中可以根據(jù)流場變化對網(wǎng)格節(jié)點屬性進(jìn)行調(diào)整，這就極大地提高了激波裝配的靈活性。

圖1 激波節(jié)點標(biāo)識Fig.1 Labels of shock nodes

2) 網(wǎng)格節(jié)點參數(shù)確定

采用格心型有限體積方法時，流場參數(shù)是儲存在單元格心的。網(wǎng)格節(jié)點參數(shù)V=[ρ,ux,uy,uz,p]T需要利用格心的流場參數(shù)U=[ρ,ux,uy,uz,p]T通過重構(gòu)獲得。其中，ρ、p分別表示流場密度和壓力；ux、uy、uz為對應(yīng)于3個方向的流動速度。按照步驟1)中的分類，對于激波節(jié)點和普通節(jié)點2種節(jié)點類型需要2種不同的計算方法。

普通網(wǎng)格節(jié)點上的參數(shù)主要通過格心參數(shù)加權(quán)平均得到

(1)

式中：Vi表示節(jié)點i上流場參數(shù)；N為與網(wǎng)格節(jié)點i相鄰的單元數(shù)；αik和Uik分別為與之相鄰的第k個單元影響的權(quán)系數(shù)和格心參數(shù)。αik采用節(jié)點與單元距離|Rik|的反比進(jìn)行確定：

(2)

式中：xi和xk為節(jié)點i和與之相鄰的第k個單元的格心坐標(biāo);Rik為第k個單元與節(jié)點i的距離矢量。

基于單元格心參數(shù)，利用式(1)、式(2)確定節(jié)點參數(shù)的方式可以看出，對于普通節(jié)點而言，其影響域包含了四周全部的單元。

使用裝配方法對激波間斷處理時需要用到上下游的兩組參數(shù)。這個特點決定了被標(biāo)識的激波節(jié)點與普通節(jié)點的一個重要區(qū)別：包含2組參數(shù)，一組代表上游低壓區(qū)的影響，一組代表下游高壓區(qū)的影響，如圖2所示。對應(yīng)于上下游，確定其激波節(jié)點2組參數(shù)所使用的網(wǎng)格單元或影響域和普通節(jié)點有所不同，上游節(jié)點參數(shù)Vu只使用上游低壓區(qū)單元，下游節(jié)點參數(shù)Vd只使用下游高壓區(qū)單元

(3)

式中：下標(biāo)“u”“d”分別表示激波上下游；N1、N2為以激波面為分界包含激波節(jié)點的上下游單元個數(shù)。與式(1)中N不同的是，在求解式(3)時并非將所有與激波節(jié)點相鄰的單元都納入考慮范圍，而是要求影響單元需要包含激波面元。此外，加權(quán)系數(shù)αik的確定也有所不同，采用式(4)確定:

圖2 激波節(jié)點參數(shù)確定Fig.2 Determination of states on shock nodes

(4)

(5)

式中：mτk為第k個單元沿所包含激波面切向的馬赫數(shù)，即

(6)

式中：定義u=[ux,uy,uz]T為流場速度矢量，則uk為第k單元的格心流動速度；xc為其所包含的激波面元運(yùn)動速度；R′為對應(yīng)激波面元中心到激波節(jié)點的空間位置矢量；定義a為聲速，則ak為單元k個單元的格心聲速。

利用光滑區(qū)域內(nèi)單元格心參數(shù)確定激波節(jié)點上下游參數(shù)的過程建立了解析解和數(shù)值解之間的聯(lián)系。對于激波間斷而言，利用激波節(jié)點上下游參數(shù)通過求解R-H關(guān)系式可以得到激波節(jié)點的運(yùn)動速度。由于R-H關(guān)系式是一維的，所以三維激波裝配還涉及到三維問題和一維方程之間的匹配問題。通過建立基于激波節(jié)點法向的激波坐標(biāo)系，可以將三維問題降為一維進(jìn)行求解。

第i個激波節(jié)點法向ni是通過與之相鄰的激波面元法向確定：

(7)

式中：Ns為與激波節(jié)點i相鄰的激波面元的個數(shù)；nj為第j個激波面元的法向，并且定義激波方向由上游指向下游為正，加權(quán)系數(shù)α′ij的確定與式(4)類似：

(8)

式中：|R′ij|表示第j個激波面元中心到激波節(jié)點i之間的距離大小。

設(shè)定激波節(jié)點運(yùn)動速度大小為w，方向為n。在激波坐標(biāo)系下，激波節(jié)點上游流動的相對馬赫數(shù)為

(9)

根據(jù)上游相對馬赫數(shù)，利用R-H關(guān)系式并補(bǔ)充反向傳播黎曼變量保持不變的條件可以確定下游參數(shù)修正值V′d以及激波運(yùn)動速度

(10)

式中:前4式為補(bǔ)充了切向速度不變條件的R-H關(guān)系式，第5式是根據(jù)反向傳播黎曼變量保持不變條件給出，γ表示比熱比。求解上述非線性方程組可以得到激波節(jié)點運(yùn)動速度w以及下游參數(shù)修正值V′d。一般情況下，V′d與Vd不相同，修正值V′d主要用于下游單元格心梯度的求解。

3) 激波動態(tài)跟蹤

本文使用基于彈簧模型的非結(jié)構(gòu)動網(wǎng)格技術(shù)實現(xiàn)對激波的顯式跟蹤。根據(jù)式(10)求解的激波節(jié)點速度w以及運(yùn)動方向n，根據(jù)t時刻激波激波節(jié)點位置xt計算t+1時刻的激波節(jié)點位置

xt+1=xt+wnΔt

(11)

式中：Δt為時間步長。根據(jù)激波節(jié)點位移，利用動網(wǎng)格技術(shù)移動其他可移動的網(wǎng)格節(jié)點，以保證激波運(yùn)動變形過程中計算網(wǎng)格的質(zhì)量。圖3給出了從初始時刻到收斂時刻激波面的運(yùn)動示意圖。運(yùn)動過程中，網(wǎng)格質(zhì)量較差時還需要配合網(wǎng)格重構(gòu)的策略。

(4) 激波通量計算

采用格心型有限體積求解器對控制方程進(jìn)行求解時需要計算流經(jīng)單元面元的通量值。對于普通面元，可以采用各種成熟的迎風(fēng)格式或中心差分格式進(jìn)行計算。而對于激波面元來說，通量計算反而變得更加簡單和清晰。在激波坐標(biāo)系下，激波面元上游為超聲速流動，按照超聲速流動特征，面元通量Φ只和上游流動有關(guān)，因此可以根據(jù)式(12)進(jìn)行確定：

圖3 激波節(jié)點運(yùn)動示意圖Fig.3 Sketch of shock node motion

(12)

式中：φ1=ρ(u-xc)·n，φ2=xc·n，E=0.5·u·u+p/[ρ(γ-1)]。n=[nx,ny,nz]T為激波面元法向，激波面元運(yùn)動速度xc可以根據(jù)激波節(jié)點運(yùn)動速度wi和法向ni得到

(13)

式中：N為三維激波面元所包含的激波節(jié)點個數(shù)。

2 格心型有限體積方法

三維Euler方程的積分形式可以寫為

(14)

式中：Q為守恒變量；Fc(Q)為無粘對流通量；Ω和?Ω分別為控制體單元及其外邊界。

在選定的控制體上，將式(14)中的對流項進(jìn)行離散可以得到

(15)

式中：Fck、nk、Sk分別為控制單元的第k個表面邊界的通量、外法向單位矢量和面積；Nf為控制體的表面?zhèn)€數(shù)。

采用格心型有限體積方法求解上述方程的過程涉及到空間物理量的重構(gòu)?？臻g二階格式通常假設(shè)物理量在控制單元內(nèi)為線性分布，即控制單元內(nèi)梯度為常數(shù)?？刂茊卧砻嫖锢砹坎捎貌煌瑔卧裥闹抵貥?gòu)得到的結(jié)果在表面邊界處形成間斷，采用van Leer等格式進(jìn)行通量計算，可以獲得通過單元表面的對流通量值，對各個控制單元邊界遍歷求和從而獲得單元內(nèi)參數(shù)變化量。然后對時間導(dǎo)數(shù)項進(jìn)行離散，就可獲得下一時間步的格心處流動參數(shù)。關(guān)于求解器的詳細(xì)介紹可以參考文獻(xiàn)[31]，在此不再贅述。

3 三維算例驗證

3.1 三維弓形激波

高超聲速球柱體繞流是考核計算方法對強(qiáng)間斷問題計算能力的一個重要算例。本文首先使用三維激波裝配方法對來流馬赫數(shù)10條件下的該問題進(jìn)行分析。球柱模型設(shè)計如圖4所示，球頭半徑為L，柱體與圓柱相接，長度也為L。計算域選為6L×3L×6L的平行六面體。

不考慮攻角和側(cè)滑角情況下，球頭繞流為軸對稱的，流動結(jié)構(gòu)相對簡單，有很多文獻(xiàn)和數(shù)據(jù)可以用來對本文計算結(jié)果進(jìn)行考核。

根據(jù)研究內(nèi)容需要，設(shè)計了3套用于球柱體計算的網(wǎng)格，分別為細(xì)網(wǎng)格C-1、粗網(wǎng)格F-1和細(xì)網(wǎng)格F-2，網(wǎng)格節(jié)點和單元數(shù)如表1所示。細(xì)網(wǎng)格C-1使用均勻劃分的方式，首先對物面和計算域外邊界使用三角形進(jìn)行離散，控制參數(shù)設(shè)為Δ=0.05L，根據(jù)邊界網(wǎng)格生成內(nèi)部計算網(wǎng)格，計算域內(nèi)全部采用四面體進(jìn)行劃分。細(xì)網(wǎng)格C-1用于激波捕捉方法的計算。為了對激波點進(jìn)行標(biāo)記定義，在激波裝配計算網(wǎng)格F-1和F-2中預(yù)先設(shè)置一個空間面，對該面進(jìn)行離散后作為內(nèi)部邊界控制區(qū)域網(wǎng)格生成，并在計算中將面上的節(jié)點定義為激波點。采用與C-1同樣的網(wǎng)格生成方式，生成粗網(wǎng)格F-1和細(xì)網(wǎng)格F-2。其中，粗網(wǎng)格F-1中物面和激波面上的網(wǎng)格離散控制參數(shù)設(shè)為Δ=0.1L，外邊界控制參數(shù)設(shè)為Δ=0.3L；細(xì)網(wǎng)格F-2中對物面以及激波層內(nèi)進(jìn)行加密，網(wǎng)格離散控制參數(shù)設(shè)為Δ=0.05L，外邊界控制參數(shù)同樣設(shè)為Δ=0.3L。用于激波捕捉計算的網(wǎng)格C-1在計算過程中不發(fā)生網(wǎng)格變形，而網(wǎng)格F-1和F-2會隨著激波的收斂過程發(fā)生變形重構(gòu)。下面對使用3套網(wǎng)格和2種方法生成的計算結(jié)果進(jìn)行對比分析。

圖4 球柱體: 計算區(qū)域及網(wǎng)格Fig.4 Hemisphere-cylinder: computational domain and mesh

1) 激波捕捉和激波裝配對比

本文對激波裝配和激波捕捉2種結(jié)果對比時都是采用同樣的求解器和同樣的輸入?yún)?shù)進(jìn)行計算，不同之處在于激波裝配方法中存在激波節(jié)點的定義，網(wǎng)格允許運(yùn)動。

圖5給出了對稱面上的壓力分布的計算結(jié)果。其中上側(cè)為采用裝配方法得到的結(jié)果，下側(cè)為采用捕捉方法計算的結(jié)果，黑色線條為文獻(xiàn)[30]中給出的壓力分布。從對比結(jié)果可以明顯看出，采用裝配方法得到的結(jié)果和文獻(xiàn)符合得非常好，壓力分布完全一致；而采用捕捉方法的計算結(jié)果存在一個較厚的激波帶，激波帶對激波層內(nèi)壓力分布存一定影響，得到的壓力分布明顯前移。進(jìn)一步，圖6中對比了2種方法計算結(jié)果在物面上的壓力分布。從對比可以明顯看出，激波捕捉方法計算結(jié)果存在壓力分布前移的現(xiàn)象，對應(yīng)位置壓力略低于激波裝配的結(jié)果。此外，使用激波捕捉得到的壓力分布較不規(guī)則，這主要是由于采用非結(jié)構(gòu)網(wǎng)格離散時激波和網(wǎng)格面的斜交程度不同導(dǎo)致環(huán)向存在偏差。

表1 球柱體計算網(wǎng)格Table 1 Computational meshes for hemisphere-cylinder

圖5 無量綱壓力等值線云圖Fig.5 Normalized pressure iso-contour levels

圖6 球頭表面的壓力等值線圖Fig.6 Pressure iso-contour levels on body

2) 激波裝配中粗細(xì)網(wǎng)格結(jié)果對比

使用激波裝配方法分別基于網(wǎng)格F-1和F-2進(jìn)行了計算。圖7給出了對稱面上馬赫數(shù)計算結(jié)果的對比。為了更加直觀地顯示F-1和F-2空間網(wǎng)格的分布情況，圖7還給出了相應(yīng)的背景網(wǎng)格，上側(cè)為粗網(wǎng)格F-1，下側(cè)為細(xì)網(wǎng)格F-2?？梢钥闯?，在激波層內(nèi)細(xì)網(wǎng)格F-2要比粗網(wǎng)格F-1密得多，但基于2套網(wǎng)格得到的激波位置完全吻合，馬赫數(shù)分布無較大差異，特別是在球頭區(qū)域。對于馬赫數(shù)大于2區(qū)域內(nèi)計算結(jié)果的偏差主要是由于粗細(xì)網(wǎng)格對于膨脹波波系結(jié)構(gòu)的描述以及流場顯示時參數(shù)擬合的控制點存在差異造成的。

圖7 激波裝配馬赫數(shù)等值線圖Fig.7 Mach iso-contours for shock-fitting

對比網(wǎng)格F-1和F-2的結(jié)果可以得出結(jié)論，激波裝配方法計算結(jié)果對網(wǎng)格粗細(xì)并不敏感，使用較少網(wǎng)格量同樣可以獲得足夠高質(zhì)量的計算結(jié)果。從這點來看，使用激波裝配方法模擬時可以適當(dāng)減少網(wǎng)格量，從而提高計算效率。

3.2 三維激波反射

球柱體繞流只包含一個弓形激波，裝配難度相對較低。進(jìn)一步，我們考慮一個較為復(fù)雜的包括一個正規(guī)反射的管內(nèi)激波反射問題。設(shè)計模型如圖8所示，在圓型通道內(nèi)放置一個圓錐模型，以圓筒直徑L為參考，圓筒半徑為R=0.5L。圓錐模型底部到理論尖點的長度為L，圓錐半錐角δ=10°。為了更加貼近實際飛行器模型，圓錐頭部選為球形，半徑為0.02L。

考慮圓筒內(nèi)來流馬赫數(shù)Ma∞=2。超聲速氣流遇圓錐產(chǎn)生脫體圓錐激波，圓錐激波遇筒壁影響會發(fā)生反射。由于流動具有軸對稱的特點，所以計算區(qū)域只選擇1/4，即圖8中實線框內(nèi)區(qū)域。計算區(qū)域采用四面體單元進(jìn)行離散，四面體單元數(shù)約為35萬。

圖8 激波反射：計算域和網(wǎng)格Fig.8 Shock reflection: computational domain and mesh

采用三維激波裝配方法對入射激波和反射激波進(jìn)行裝配。由于激波在壁面發(fā)生反射，涉及到激波壁面交點運(yùn)動參數(shù)的確定。三維激波交點運(yùn)動參數(shù)的確定要比二維復(fù)雜得多，其運(yùn)動速度求解不當(dāng)會出現(xiàn)不相容問題，導(dǎo)致計算失敗。所以激波交點運(yùn)動參數(shù)的確定對于三維裝配而言至關(guān)重要。

根據(jù)正規(guī)反射的特點，從幾何構(gòu)型的角度出發(fā)認(rèn)為激波在筒壁上反射點的運(yùn)動速度主要受入射激波影響，應(yīng)與入射激波相匹配，即與壁面交點的運(yùn)動速度wR應(yīng)由入射激波決定，可以利用式(16)進(jìn)行確定：

wR=(wI·τw)τw

(16)

式中：wI為按照入射激波確定交點參數(shù)得到的激波節(jié)點速度;τw為壁面的切向向量。

裝配方法得到的激波型面如圖9所示。為了更加直觀地顯示激波面形狀將計算結(jié)果進(jìn)行對稱，可以看出激波面光滑合理，定性說明裝配的正確性。對計算結(jié)果進(jìn)行一步分析，給出了裝配激波的角度。根據(jù)氣體動力學(xué)相關(guān)知識可知，超聲速軸對稱繞流形成的附體圓錐激波從每個子午面來看都是斜激波。在馬赫數(shù)Ma∞=2、半錐角δ=10°條件下，激波角β=31.2°。從圖9的對比可以看出，由于圓錐頭部為球形，激波會脫體，對下游流動造成影響。隨著距離增加，激波角減小，最后收斂至錐形流理論激波角31.2°。入射激波壁面反射為正規(guī)反射，反射激波角度可以根據(jù)入射激波角31.2°和來流馬赫數(shù)Ma∞=2進(jìn)行確定。經(jīng)計算反射激波角為β′=30.6°，與裝配激波角度吻合較好。值得說明的是，遠(yuǎn)離洞壁時反射激波上游并非勻直流，反射激波角度會出現(xiàn)一定偏差。

圖9 激波反射：裝配得到的激波Fig.9 Shock reflection: fitted shock waves

3.3 三維激波-激波相交

最后考慮一個更為復(fù)雜的含有激波相交的流動問題—球柱錐組合體超聲速繞流問題。針對該問題Houtmann[32]、Bonfiglioli[30]分別進(jìn)行了實驗和數(shù)值研究。選用他們研究的模型和參數(shù)，如圖10所示，使用三維激波裝配方法進(jìn)行模擬計算。具體模型尺寸如圖10所示，來流馬赫數(shù)Ma∞=4.04，攻角為α=20°。該問題涉及到球頭和拐角處形成的脫體激波和附體激波的裝配以及對兩個激波形成的三維激波相交的裝配，裝配難度較大。

根據(jù)流動特點，計算區(qū)域選擇如圖11所示。采用完全的四面體網(wǎng)格對計算區(qū)域進(jìn)行離散，離散控制參數(shù)設(shè)為Δ=2 mm，離散網(wǎng)格量約為94萬，用于激波捕捉方法使用的計算網(wǎng)格。如前所述，使用激波裝配方法計算時為了對激波點進(jìn)行標(biāo)記定義，預(yù)先設(shè)置激波面，通過激波面和邊界對計算區(qū)域進(jìn)行網(wǎng)格劃分，網(wǎng)格量約為110萬。

圖10 模型尺寸圖Fig.10 Body geometry and dimension

圖11 激波相交：計算域及網(wǎng)格Fig.11 Shock interaction: Computational domain and mesh

在來流條件作用下最外層形成的脫體弓形激波與錐裙和圓柱的結(jié)合處形成的內(nèi)嵌激波相交形成VI類激波相交構(gòu)型，如圖12所示。如前文所述，對于該問題進(jìn)行裝配的難點和重點依然是如何確定激波交點運(yùn)動速度。根據(jù)流動特點，本文給出了一種適用于此類三維激波交點運(yùn)動速度的確定方法。通過對該類型激波結(jié)構(gòu)的特點進(jìn)行分析可知，激波交點T的運(yùn)動速度由處于上游的2個激波I1和I2的運(yùn)動決定。利用此關(guān)系，假設(shè)I1和I2的激波運(yùn)動速度分別為w1和w2，運(yùn)動方向分別為n1和n2，激波交點T的運(yùn)動速度w可以通過求解式(17)確定：

(17)

式中：n=n1×n2為平面Π的法向向量。

利用式(17)對交點運(yùn)動進(jìn)行計算，得到收斂的激波裝配結(jié)果，如圖13所示。從結(jié)果來看，激波型面光滑合理，根據(jù)控制方程解唯一性特點可以在一定程度上反映裝配結(jié)果的正確性。利用三維激波裝配方法得到的激波型面使得對三維激波的認(rèn)識更加直觀。

圖14給出了激波裝配方法和激波捕捉方法在對稱面上的馬赫數(shù)云圖對比，從對比結(jié)果可以看出：激波裝配和激波捕捉得到的流場結(jié)果分布基本一致，也反映了裝配結(jié)果的合理性，從而驗證了三維激波裝配方法適用于對復(fù)雜激波構(gòu)型問題的裝配。進(jìn)一步對比可以看出，捕捉方法得到的激波存在一個較厚的激波帶，隨著激波強(qiáng)度的變?nèi)鹾穸冗€有增大的趨勢，而激波裝配得到的激波為理論上的厚度為零。

圖12 激波結(jié)構(gòu)示意圖Fig.12 Sketch of shock structure

圖14 對稱面上馬赫數(shù)云圖Fig.14 Mach number iso-contours on symmetric plane

4 結(jié) 論

本文將基于非結(jié)構(gòu)動網(wǎng)格的激波裝配法拓展應(yīng)用于對三維問題的研究，數(shù)值模擬了若干含有不同激波個數(shù)以及不同激波結(jié)構(gòu)的超聲速/高超聲速流動問題。計算結(jié)果顯示該方法具有適用于三維問題的能力。通過對高超聲速球柱繞流問題的模擬驗證了該裝配方法在三維問題中的合理性。在激波壁面反射和激波相交問題中，通過對激波結(jié)構(gòu)的幾何關(guān)系的分析發(fā)展了適用于描述這兩類三維激波交點運(yùn)動的方法。理論上，在其他類型的三維激波相交問題中利用激波構(gòu)型的幾何特征確定激波交點速度的方法也是適用的。通過本文工作，對傳統(tǒng)裝配方法進(jìn)行了改進(jìn)，拓展了裝配方法的應(yīng)用范圍。三維激波裝配可以獲得合理的激波型面，增加了對三維激波形狀的直觀理解，對于相關(guān)理論方法的研究也具有一定參考意義。

激波裝配方法發(fā)展的一個重要目的就是解決其通用性問題，實現(xiàn)其在更為復(fù)雜的工程問題中的應(yīng)用價值。針對三維問題進(jìn)行研究，會發(fā)現(xiàn)二維裝配所遇不到的問題，以解決相關(guān)問題為導(dǎo)引有利于加速發(fā)展出更為通用的激波裝配方法。本文關(guān)于三維激波裝配方法的研究正是沿著這一思路而來。但是從目前來看，關(guān)于三維激波裝配仍然還有很多工作要做。本文提出的激波交點的計算方法是利用激波結(jié)構(gòu)的幾何關(guān)系獲得，背后的物理含義仍然需要作大量的探索研究。對于基于非結(jié)構(gòu)動網(wǎng)格的激波裝配方法，給定的初始激波和收斂時的激波位置距離較遠(yuǎn)時，激波運(yùn)動幅度較大，需要不斷進(jìn)行網(wǎng)格重構(gòu)，嚴(yán)重影響計算效率。精確的激波辨識可以提供相對準(zhǔn)確的激波面，有助于提高計算效率。關(guān)于激波面的辨識依然是三維激波問題的一個重點。諸如此類的問題在一定程度上都還制約著激波裝配方法的實際應(yīng)用，有待后續(xù)深入研究。