環(huán)形激波聚焦流場(chǎng)特性的數(shù)值研究＊

陳二云，趙改平，楊?lèi)?ài)玲

（1.上海理工大學(xué)能源與動(dòng)力工程學(xué)院，上海 200093；2.上海理工大學(xué)醫(yī)療器械學(xué)院，上海 200093）

激波聚焦［1］是指在一定條件下激波在其傳播方向上發(fā)生收斂的行為。激波聚焦產(chǎn)生的瞬態(tài)高溫高壓脈沖在工程技術(shù)領(lǐng)域具有重要的應(yīng)用價(jià)值，如第一臺(tái)體外激波粉碎腎結(jié)石的新型醫(yī)療機(jī)械就利用了激波聚焦過(guò)程產(chǎn)生高能密度的力學(xué)特性。激波聚焦的一種簡(jiǎn)單實(shí)現(xiàn)方法是讓平面入射激波在凹形壁面反射實(shí)現(xiàn)聚焦。而環(huán)形激波聚焦不需要反射物面，具有三維會(huì)聚的特點(diǎn)，能獲得更好的聚焦效果［2］。

激波聚焦的研究方法主要有理論分析、實(shí)驗(yàn)測(cè)試和數(shù)值模擬。T.Satio［3］從理論、計(jì)算和實(shí)驗(yàn)等3個(gè)不同的角度對(duì)半球形弱沖擊波的聚焦問(wèn)題進(jìn)行了研究。S.Hamid等［4］利用無(wú)膜激波管對(duì)低馬赫數(shù)環(huán)形激波在激波管內(nèi)的聚焦情況進(jìn)行了實(shí)驗(yàn)研究，獲得了流場(chǎng)在不同時(shí)刻的全息干涉條紋和壁面各點(diǎn)壓力時(shí)程曲線(xiàn)。對(duì)于高馬赫數(shù)下環(huán)形激波聚焦實(shí)驗(yàn)測(cè)試技術(shù)比較困難，因此數(shù)值模擬是一種較好的研究手段。目前，對(duì)于此類(lèi)流動(dòng)的數(shù)值模擬方法主要有DCD［5］、FVM［6］和 WENO［7］格式，能夠較好的模擬出激波聚焦過(guò)程中的復(fù)雜流場(chǎng)結(jié)構(gòu)。但這些方法都必須通過(guò)擴(kuò)大節(jié)點(diǎn)模板來(lái)提高解的精度，因此在求解包含復(fù)雜幾何外形的流動(dòng)問(wèn)題時(shí)，不利于程序的編制和邊界條件的處理。

B.Cockburn等［8］、J.Qiu等［9］、N.Chevaugeon［10］提出了一種間斷有限元方法。該方法吸收了有限元方法和有限體積方法的優(yōu)點(diǎn)，具有一致高階精度、網(wǎng)格適應(yīng)性強(qiáng)、結(jié)構(gòu)守恒性和容易向高維推廣等特點(diǎn)。與其它方法相比，該方法最大的優(yōu)點(diǎn)在于解的精度是通過(guò)提高插值多項(xiàng)式的階數(shù)來(lái)實(shí)現(xiàn)，無(wú)需擴(kuò)大節(jié)點(diǎn)模板，便于程序的編制和邊界條件的處理。

本文中，在已有工作的基礎(chǔ)上，采用間斷有限元方法模擬了環(huán)形激波在同軸圓柱形激波管內(nèi)的聚焦流場(chǎng)特性，分析環(huán)形管道內(nèi)外半徑對(duì)聚焦峰值壓力的影響作用。計(jì)算的結(jié)果可以為實(shí)際的工程應(yīng)用提供一定的理論指導(dǎo)。

1 間斷有限元方法

1.1 控制方程

忽略粘性和熱傳導(dǎo)效應(yīng)的假設(shè)下，軸對(duì)稱(chēng)Euler方程組表示為

式中：U＝，t是時(shí)間變量，ρ是密度，u和v分別為對(duì)應(yīng)于x和y方向的速度分量，p是壓強(qiáng)，E＝ρe＋ρ（u2＋v2）／2，表示單位體積的總能量，e是單位質(zhì)量?jī)?nèi)能，對(duì)于完全氣體，狀態(tài)方程p＝（γ－1）ρe，γ是比熱容比。

1.2 數(shù)值離散

在數(shù)值模擬中，對(duì)式（1）的空間離散形式采用間斷有限元方法。設(shè)Γh是區(qū)域Ω的一個(gè)有限剖分，單元K∈Γh，?K表示單元K的邊界，n表示單元邊界的外法向。將間斷有限元的局部解空間定義為單元上次數(shù)不超過(guò)2的多項(xiàng)式集合，記為V（K）。對(duì)于任意時(shí)刻t，在間斷有限元空間Vh＝中尋找近似解Uh（X，t），其中X＝（x，y）。

首先在單元K上用連續(xù)函數(shù)vh∈Vh乘式（1）的兩端，并用他的近似解Uh（X，t）∈Vh代替式（1）中的精確解U（X，t），由Green公式，得

式中：F＝，用數(shù)值通量H（X，t）代替通量F（Uh（X，t））·n，得

數(shù)值通量定義為H（X，t）＝H（Uh（Xint（K），t），Uh（Xext（K），t）），體現(xiàn)了與相鄰單元之間的信息傳遞，同時(shí)也保證了離散格式的守恒性。

式（3）中的積分用數(shù)值積分代替

最后得到弱表達(dá)式

為了計(jì)算的方便，在單元K中取正交基函數(shù)（如勒讓得多項(xiàng)式）｛φ1，φ2，φ3，…，φj｝，則質(zhì)量矩陣成為分塊對(duì)角矩陣，故有限元解可以表示為

令式（7）中vh（X）＝φi（X），得

設(shè)MK為單元K的質(zhì)量矩陣，則vh∈V（K），?K∈Γh，式（9）可寫(xiě)成常微分方程形式

式中：算子Lh表示對(duì)（－?·F＋S）的近似離散，γh為跡函數(shù)。在時(shí)間推進(jìn)方面，為了與空間離散保持一致高階精度，采用三階Runge－Kutta時(shí)間離散形式［8］。

2 計(jì)算結(jié)果與分析

計(jì)算區(qū)域如圖1所示，其中L表示圓柱形激波管半徑，D與d分別表示環(huán)形管道的內(nèi)外半徑。初始時(shí)刻，馬赫數(shù)為3.0的環(huán)形激波從環(huán)形管道向同軸圓柱形管道傳播，波前氣體處于靜止?fàn)顟B(tài)，網(wǎng)格步長(zhǎng)Δx＝Δy＝1／500。

圖2給出了當(dāng)L∶D∶d＝1∶0.4∶0.2時(shí)，環(huán)形激波從環(huán)形管道向同軸圓柱形管道傳播過(guò)程的壓力等值線(xiàn)分布。

圖1 計(jì)算區(qū)域示意圖Fig.1Computation domain

圖2 壓力等值線(xiàn)分布 （Ma＝3.0）Fig.2Distributions of pressure contours（Ma＝3.0）

從圖中可以看出，當(dāng)激波進(jìn)入圓柱形管道后，由于傳播截面突然增大，激波發(fā)生繞射，繞射波在臺(tái)階拐角處產(chǎn)生2個(gè)很強(qiáng)的漩渦，如圖2（a）所示。當(dāng)繞射激波到達(dá)對(duì)稱(chēng)軸時(shí)，激波聚焦在對(duì)稱(chēng)軸上，在聚焦點(diǎn)附近形成局部的高溫高壓區(qū)域，且在漩渦附近出現(xiàn)了二次激波，同時(shí)激波在對(duì)稱(chēng)軸上發(fā)生由規(guī)則反射向馬赫反射的轉(zhuǎn)變，如圖2（b）所示。激波在對(duì)稱(chēng)軸上發(fā)生馬赫反射之后，馬赫桿受到激波聚焦誘發(fā)的強(qiáng)射流沖擊作用，再次發(fā)生馬赫反射，使得馬赫桿向前凸起呈半球面形狀，并且與原馬赫桿相交位置形成三波交點(diǎn)，如圖2（c）所示。隨著時(shí)間的進(jìn)一步發(fā)展，對(duì)稱(chēng)軸附近的球面狀馬赫桿繼續(xù)凸起，前端已經(jīng)超越繞射激波的波陣面，后方有1個(gè)明顯的漩渦，如圖2（d）所示，這種激波結(jié)構(gòu)又被稱(chēng)為球面雙馬赫反射，是環(huán)形激波聚焦過(guò)程中產(chǎn)生的典型流場(chǎng)結(jié)構(gòu)。

圖3（a）～（c）分別給出了L∶D∶d＝1∶0.4∶0.2，L∶D∶d＝1∶0.5∶0.2和L∶D∶d＝1∶0.5∶0.3，環(huán)形激波從環(huán)形管道向同軸圓柱形管道傳播過(guò)程中的壓力等值線(xiàn)分布。從圖中可以看出，3種不同計(jì)算工況下捕捉到的流場(chǎng)結(jié)構(gòu)特點(diǎn)（包括漩渦、二次激波、三波交點(diǎn)和球面雙馬赫反射等）基本是相似的，但也存在一定的差異。比較圖3（a）和圖3（b）可以發(fā)現(xiàn)，外側(cè)臺(tái)階后面漩渦的位置隨著環(huán)形管道外徑的增大沿徑向外移，靠近中心臺(tái)階后面的流場(chǎng)結(jié)構(gòu)非常相似，當(dāng)x＞0.5時(shí)，圖3（b）中高壓區(qū)域的覆蓋面積明顯增大。比較圖3（b）和圖3（c）可以發(fā)現(xiàn)，中心臺(tái)階后面的流場(chǎng)結(jié)構(gòu)具有較大差別，隨著環(huán)形激波管內(nèi)徑變大，環(huán)形激波在中心軸線(xiàn)的聚焦點(diǎn)向下游偏移，高壓區(qū)域也相應(yīng)地向下游偏移。

圖3 壓力等值線(xiàn)分布 （Ma＝3.0，t＝0.35）Fig.3Distributions of pressure contours（Ma＝3.0，t＝0.35）

圖4（a）給出了L∶D∶d＝1∶0.4∶0.2和L∶D∶d＝1∶0.5∶0.2的情況下，環(huán)形激波在整個(gè)傳播過(guò)程中軸線(xiàn)各點(diǎn)所能達(dá)到的最大壓力分布曲線(xiàn)對(duì)比圖，其中圓圈代表L∶D∶d＝1∶0.5∶0.2時(shí)軸線(xiàn)各點(diǎn)最大壓力分布曲線(xiàn)，實(shí)線(xiàn)代表L∶D∶d＝1∶0.4∶0.2時(shí)軸線(xiàn)各點(diǎn)最大壓力分布曲線(xiàn)，pmax／p0表示中心軸線(xiàn)上各點(diǎn)最大壓力與環(huán)境靜壓之比。從圖4（a）中可以看出，其它流動(dòng)條件不變時(shí)，環(huán)形激波管外徑對(duì)中心軸線(xiàn)各點(diǎn)最大壓力值幾乎沒(méi)有影響，僅在下游較遠(yuǎn)位置二者有一定的差別。

圖4（b）給出了L∶D∶d＝1∶0.5∶0.2和L∶D∶d＝1∶0.5∶0.3情況下，環(huán)形激波在整個(gè)流動(dòng)過(guò)程中軸線(xiàn)各點(diǎn)所能達(dá)到的最大壓力分布曲線(xiàn)對(duì)比圖，其中圓圈代表L∶D∶d＝1∶0.5∶0.2時(shí)軸線(xiàn)各點(diǎn)最大壓力分布曲線(xiàn)，實(shí)線(xiàn)代表L∶D∶d＝1∶0.5∶0.3時(shí)軸線(xiàn)各點(diǎn)最大壓力分布曲線(xiàn)。從圖4（b）中可以看出，在其他流動(dòng)條件不變時(shí)，環(huán)形管道內(nèi)徑的增大使環(huán)形激波管中心軸線(xiàn)上絕大部分位置的最大壓力都有一定增加，最大峰值壓力pmax為約105p0，增加了約1.25倍，峰值壓力位置也向下游偏移。

圖4（c）給出了L∶D∶d＝1∶0.4∶0.2和L∶D∶d＝1∶0.5∶0.3情況下，環(huán)形激波在整個(gè)流動(dòng)過(guò)程中軸線(xiàn)上各點(diǎn)所能達(dá)到的最大壓力分布曲線(xiàn)對(duì)比圖，其中圓圈代表L∶D∶d＝1∶0.4∶0.2時(shí)軸線(xiàn)上各點(diǎn)最大壓力分布曲線(xiàn)，實(shí)線(xiàn)代表L∶D∶d＝1∶0.5∶0.3時(shí)軸線(xiàn)上各點(diǎn)最大壓力分布曲線(xiàn)。從圖4（c）中可以看出，中心軸線(xiàn)上各點(diǎn)最大壓力的分布曲線(xiàn)與圖4（b）相似，由此可知，環(huán)形激波在圓柱形管道內(nèi)傳播時(shí)，中心軸線(xiàn)上各點(diǎn)最大壓力值受環(huán)形管道外徑的影響較小，受環(huán)形管道內(nèi)徑的影響較大。

圖4 中心軸線(xiàn)最大壓力分布曲線(xiàn) （Ma＝3.0）Fig.4Distributions of maximum pressure along the axis（Ma＝3.0）

3 結(jié) 論

采用三階間斷有限元方法對(duì)環(huán)形激波在圓柱形激波管內(nèi)聚焦流場(chǎng)進(jìn)行數(shù)值模擬。結(jié)果表明，三階間斷有限元方法能夠較準(zhǔn)確的捕捉到激波聚焦過(guò)程中復(fù)雜流場(chǎng)波系結(jié)構(gòu)，通過(guò)改變環(huán)形管道內(nèi)外半徑對(duì)聚焦流場(chǎng)進(jìn)行研究發(fā)現(xiàn)，環(huán)形激波在圓柱形管道內(nèi)傳播時(shí)，中心軸線(xiàn)上各點(diǎn)處最大壓力值受環(huán)形管道外徑的影響較小，受環(huán)形管道內(nèi)徑的影響較大。此外，環(huán)形激波在中心軸線(xiàn)上的聚焦位置隨環(huán)形管道內(nèi)徑的增加向下游移動(dòng)。計(jì)算結(jié)果可以為實(shí)際的工程應(yīng)用提供一定的理論指導(dǎo)。

［1］董剛，葉經(jīng)方，范寶春.激波聚焦反射的實(shí)驗(yàn)和數(shù)值研究［J］.高壓物理學(xué)報(bào)，2006，20（4）：359－365.

DONG Gang，YE Jing－fang，F(xiàn)AN Bao－chun.Experimental and numerical investigation of shock wave focusing and reflection［J］.Chinese Journal of High Pressure Physics，2006，20（4）：359－365.

［2］TENG Hong－h(huán)ui，JIANG Zong－lin.Gasdynamic characteristics of toroidal shock and detonation wave converging［J］.Science in China Series G：Physics Mechanics and Astronomy，2005，48（2）：739－749.

［3］Satio T.An experimental analytical and numerical study of temperatures near hemispherical implosion foci［R］.UTIAS Report No.260，1982.

［4］Hamid S，Hosseini R，Takayama K.Study of shock wave focusing and reflection over symmetrical axis of a compact vertical co－axial diaphragmless shock tube［C］∥Proceedings of ISSW23，2001：1550－1557.

［5］滕宏輝，姜宗林，韓肇元.環(huán)形激波繞射、反射和聚焦的數(shù)值模擬研究［J］.力學(xué)學(xué)報(bào)，2004，36（1）：9－15.

TENG Hong－h(huán)ui，JIANG Zong－lin，HAN Zhao－yuan.Numerical investigation of diffraction，focusing and reflection of toroidal shock waves［J］.Acta Mechanica Sinica，2004，36（1）：9－15.

［6］滕宏輝，張德良，李輝煌，等.用環(huán)形激波聚焦實(shí)現(xiàn)爆轟波直接起爆的數(shù)值模擬［J］.爆炸與沖擊，2005，25（6）：512－518.

TENG Hong－h(huán)ui，ZHANG De－liang，LI Hui－h(huán)uang，et al.Numerical investigation of detonation direct initiation induced by toroidal shock wave focusing［J］.Explosion and Shock Waves，2005，25（6）：512－518.

［7］Liang S M，Wu L N，Hsu R L.Numerical investigation of axisymmetric shock wave focusing over paraboloidial reflectors［J］.Shock Waves，1999，9（6）：367－379.

［8］Cockburn B，Shu C W.Runge－Kutta discontinuous galerkin methods for convection－dominated problems［J］.Journal of Scientific Computing，2001，16（3）：173－261.

［9］Qiu J X，Khoo B C，Shu C W.A numerical study for the performance of the Runge－Kutta discontinuous galerkin method based on different numerical fluxes［J］.Journal of Computational Physics，2006，212（2）：540－565.

［10］Chevaugeon N，Xin J，Hu P，et al.Discontinuous galerkin methods applied to shock and blast problems［J］.Journal of Scientific Computing，2005，22－23（1／2／3）：227－243.