環形激波聚焦流場特性的數值研究＊

2012-06-20 08:22:26陳二云趙改平楊愛玲

爆炸與沖擊 2012年3期

陳二云，趙改平，楊愛玲

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

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

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

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

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

1 間斷有限元方法

1.1 控制方程

忽略粘性和熱傳導效應的假設下，軸對稱Euler方程組表示為

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

1.2 數值離散

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

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

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

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

式（3）中的積分用數值積分代替

最后得到弱表達式

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

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

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

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

2 計算結果與分析

計算區域如圖1所示，其中L表示圓柱形激波管半徑，D與d分別表示環形管道的內外半徑。初始時刻，馬赫數為3.0的環形激波從環形管道向同軸圓柱形管道傳播，波前氣體處于靜止狀態，網格步長Δx＝Δy＝1／500。

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

圖1 計算區域示意圖Fig.1Computation domain

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

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

圖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，環形激波從環形管道向同軸圓柱形管道傳播過程中的壓力等值線分布。從圖中可以看出，3種不同計算工況下捕捉到的流場結構特點（包括漩渦、二次激波、三波交點和球面雙馬赫反射等）基本是相似的，但也存在一定的差異。比較圖3（a）和圖3（b）可以發現，外側臺階后面漩渦的位置隨著環形管道外徑的增大沿徑向外移，靠近中心臺階后面的流場結構非常相似，當x＞0.5時，圖3（b）中高壓區域的覆蓋面積明顯增大。比較圖3（b）和圖3（c）可以發現，中心臺階后面的流場結構具有較大差別，隨著環形激波管內徑變大，環形激波在中心軸線的聚焦點向下游偏移，高壓區域也相應地向下游偏移。

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

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

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

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

3 結論

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

［1］董剛，葉經方，范寶春.激波聚焦反射的實驗和數值研究［J］.高壓物理學報，2006，20（4）：359－365.

DONG Gang，YE Jing－fang，FAN 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－hui，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］滕宏輝，姜宗林，韓肇元.環形激波繞射、反射和聚焦的數值模擬研究［J］.力學學報，2004，36（1）：9－15.

TENG Hong－hui，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］滕宏輝，張德良，李輝煌，等.用環形激波聚焦實現爆轟波直接起爆的數值模擬［J］.爆炸與沖擊，2005，25（6）：512－518.

TENG Hong－hui，ZHANG De－liang，LI Hui－huang，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.