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

陳二云，趙改平，楊愛玲

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

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

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

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

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

1 間斷有限元方法

1.1 控制方程

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

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

1.2 數(shù)值離散

在數(shù)值模擬中，對式（1）的空間離散形式采用間斷有限元方法。設(shè)Γh是區(qū)域Ω的一個有限剖分，單元K∈Γh，?K表示單元K的邊界，n表示單元邊界的外法向。將間斷有限元的局部解空間定義為單元上次數(shù)不超過2的多項式集合，記為V（K）。對于任意時刻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)了與相鄰單元之間的信息傳遞，同時也保證了離散格式的守恒性。

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

最后得到弱表達式

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

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

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

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

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

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

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

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

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

從圖中可以看出，當(dāng)激波進入圓柱形管道后，由于傳播截面突然增大，激波發(fā)生繞射，繞射波在臺階拐角處產(chǎn)生2個很強的漩渦，如圖2（a）所示。當(dāng)繞射激波到達對稱軸時，激波聚焦在對稱軸上，在聚焦點附近形成局部的高溫高壓區(qū)域，且在漩渦附近出現(xiàn)了二次激波，同時激波在對稱軸上發(fā)生由規(guī)則反射向馬赫反射的轉(zhuǎn)變，如圖2（b）所示。激波在對稱軸上發(fā)生馬赫反射之后，馬赫桿受到激波聚焦誘發(fā)的強射流沖擊作用，再次發(fā)生馬赫反射，使得馬赫桿向前凸起呈半球面形狀，并且與原馬赫桿相交位置形成三波交點，如圖2（c）所示。隨著時間的進一步發(fā)展，對稱軸附近的球面狀馬赫桿繼續(xù)凸起，前端已經(jīng)超越繞射激波的波陣面，后方有1個明顯的漩渦，如圖2（d）所示，這種激波結(jié)構(gòu)又被稱為球面雙馬赫反射，是環(huán)形激波聚焦過程中產(chǎn)生的典型流場結(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)形管道向同軸圓柱形管道傳播過程中的壓力等值線分布。從圖中可以看出，3種不同計算工況下捕捉到的流場結(jié)構(gòu)特點（包括漩渦、二次激波、三波交點和球面雙馬赫反射等）基本是相似的，但也存在一定的差異。比較圖3（a）和圖3（b）可以發(fā)現(xiàn)，外側(cè)臺階后面漩渦的位置隨著環(huán)形管道外徑的增大沿徑向外移，靠近中心臺階后面的流場結(jié)構(gòu)非常相似，當(dāng)x＞0.5時，圖3（b）中高壓區(qū)域的覆蓋面積明顯增大。比較圖3（b）和圖3（c）可以發(fā)現(xiàn)，中心臺階后面的流場結(jié)構(gòu)具有較大差別，隨著環(huán)形激波管內(nèi)徑變大，環(huán)形激波在中心軸線的聚焦點向下游偏移，高壓區(qū)域也相應(yīng)地向下游偏移。

圖3 壓力等值線分布 （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)形激波在整個傳播過程中軸線各點所能達到的最大壓力分布曲線對比圖，其中圓圈代表L∶D∶d＝1∶0.5∶0.2時軸線各點最大壓力分布曲線，實線代表L∶D∶d＝1∶0.4∶0.2時軸線各點最大壓力分布曲線，pmax／p0表示中心軸線上各點最大壓力與環(huán)境靜壓之比。從圖4（a）中可以看出，其它流動條件不變時，環(huán)形激波管外徑對中心軸線各點最大壓力值幾乎沒有影響，僅在下游較遠位置二者有一定的差別。

圖4（b）給出了L∶D∶d＝1∶0.5∶0.2和L∶D∶d＝1∶0.5∶0.3情況下，環(huán)形激波在整個流動過程中軸線各點所能達到的最大壓力分布曲線對比圖，其中圓圈代表L∶D∶d＝1∶0.5∶0.2時軸線各點最大壓力分布曲線，實線代表L∶D∶d＝1∶0.5∶0.3時軸線各點最大壓力分布曲線。從圖4（b）中可以看出，在其他流動條件不變時，環(huán)形管道內(nèi)徑的增大使環(huá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)形激波在整個流動過程中軸線上各點所能達到的最大壓力分布曲線對比圖，其中圓圈代表L∶D∶d＝1∶0.4∶0.2時軸線上各點最大壓力分布曲線，實線代表L∶D∶d＝1∶0.5∶0.3時軸線上各點最大壓力分布曲線。從圖4（c）中可以看出，中心軸線上各點最大壓力的分布曲線與圖4（b）相似，由此可知，環(huán)形激波在圓柱形管道內(nèi)傳播時，中心軸線上各點最大壓力值受環(huán)形管道外徑的影響較小，受環(huán)形管道內(nèi)徑的影響較大。

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

3 結(jié) 論

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

［1］董剛，葉經(jīng)方，范寶春.激波聚焦反射的實驗和數(shù)值研究［J］.高壓物理學(xué)報，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é)報，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)形激波聚焦實現(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.