微型激波管內部激波特性的數值模擬

張 光，崔寶玲，金英子，金羲東

（1.浙江理工大學浙江省流體傳輸技術研究重點實驗室，杭州310018；2.安東國立大學機械工程學院，韓國安東760749）

微型激波管內部激波特性的數值模擬

張 光1，崔寶玲1，金英子1，金羲東2

（1.浙江理工大學浙江省流體傳輸技術研究重點實驗室，杭州310018；2.安東國立大學機械工程學院，韓國安東760749）

為了研究微型激波管內部的不穩定流動和激波運動特性，采用數值模擬的方法對微型激波管內部流動進行分析。對比分析不同隔膜壓力比（高壓腔與低壓腔的初始壓力之比）和激波管直徑對微型激波管內激波、交接面以及流體運動特性的影響，并與實驗數據進行比較。結果表明：隨著隔膜壓力比的增大，激波和交接面的運動速度逐漸增大；激波在微型激波管內運動時，其強度逐漸減弱；在低壓被驅動腔內壓強較低時，觀察到厚度較大的邊界層，這說明低壓影響對微型激波管內的激波和流體運動會產生一定的能量損失；激波前后的壓力梯度隨著激波運動逐漸減??；交接面在微型激波管內運動時，運動速度逐漸增大；參數S值可以反映微型激波管內的低壓和小尺寸影響。

微型激波管；激波；交接面；非定常流動；數值模擬

0　引 言

近年來，微型激波管已廣泛用于機械和生物工程設備中，如微型發動機、微型燃燒機、無針注射等。與傳統的激波管類似，微型激波管由隔膜隔開的高壓腔和低壓腔組成，高壓腔也稱為驅動腔，低壓腔稱為被驅動腔。由于高壓腔與低壓腔之間的壓力差，隔膜的瞬間破裂會使激波管內產生不穩定的激波流動。如果壓力足夠大，隔膜將會自行破裂，否則需要用手動的方式使隔膜破裂。隔膜破裂產生的激波和交接面會向低壓腔運動，而產生的膨脹波向高壓腔運動。在微型激波管內，含有運動激波的可壓縮流動需要更多地考慮耗散帶來的影響，而在大型激波管內可以忽略這種影響。

在過去的一個世紀里，盡管許多學者不斷努力研究微型激波管內部流動，但仍然存在一些沒有探究和無法解釋的現象。在微型激波管內，由低壓和小尺寸產生的稀疏影響表現得尤為明顯，該影響使激波強度逐漸衰弱，而在大型激波管內可以忽略這種影響；流體的粘性影響和激波與壁面之間的摩擦也使激波產生更多的能量損失。因此，關于產生激波強度損失的原因一直是國內外學者研究的重點［1-6］。文獻［1-3］采用試驗的方法探究了微型激波管內激波強度衰減的原因，結果表明：激波管內粘性影響、邊界層的發展以及壁面摩擦是導致試驗過程中激波強度損失的主要原因。微型激波管內邊界層的發展使激波和交接面的運動與理論預測產生一定的偏差。交接面是假想隔開高壓腔氣體與被激波加熱氣體的一個面。在交接面兩邊，氣體壓強和速度大小相等，而溫度和密度不同。在邊界層的作用下，激波運動速度逐漸減小，而交接面速度逐漸增大。Sturtevant等［7］、Tanaki等［8］運用數值模擬的方法探究了不同邊界條件下邊界層對激波運動的影響；Mirels［9］和Roshko［10］采用理論分析的方法研究了微型激波管內的邊界層效應，并提出了能夠預測邊界層效應的理論模型。

與傳統的大型激波管相比，尺寸對微型激波管內激波和流體的運動特性影響極大。Brouillette［11］采用試驗方法探究了微小尺寸激波管內可壓縮的湍流，引入控制體法，并提出表明小尺寸影響的一個參數S。在S較小的條件下，分析了微型激波管內摩擦和熱傳遞產生的擴散影響，結果表明：控制體模型能夠有效預測小尺寸對激波強度的影響。但作者沒有給出S值隨時間的變化情況。Park等［12］和Sun等［13］運用試驗的方法研究了不同直徑的微型激波管內激波運動，所得到的實驗結果與理論預測基本一致。Zeitoun等［14-15］采用數值模擬的方法探究了包含與速度滑移和溫度跳躍的邊界條件相耦合的湍流納維-斯托克斯方程在微型激波管內部流動方面的應用。Arun等［16-17］對不同邊界條件下的微型激波管進行數值模擬，并與理論預測進行比較，數值模擬結果與理論結果相吻合。

本文采用數值模擬的方法對微型激波管內激波和交接面的運動特性進行分析，并揭示不同邊界條件下激波強度損失的原因。采用CFD軟件Fluent進行數值模擬。保持高壓腔的壓力為大氣壓強不變，改變低壓腔的壓力，從而探究不同隔膜壓力比對激波運動的影響；采用不同直徑的激波管來研究尺寸對激波管內部流動特性的影響；將數值模擬結果與實驗結果進行對比，從而驗證數值模擬方法的正確性；通過計算得到激波與交接面之間的距離，進而得到反應尺寸影響的參數S的變化情況。

1　計算模型及數值方法

1.1計算模型

本文中數值模擬采用的計算模型如圖1所示。采用2維軸對稱模型，高壓腔和低壓腔都采用理想氣體作為工作流體，高壓腔氣體壓強為P4，低壓腔氣體壓強為P1。定義隔膜壓力比（diaphragm pressure ratio，DPR）為P4/P1。管內截面為圓形，直徑用D來表示，高壓腔長度L4為100 mm，低壓腔長度L1為200 mm。X表示軸線方向的位移，向右為正方向。

圖1　數值模擬的計算模型

2.2數值方法

為了研究微型激波管內部詳細的非定常流動，本文采用結構化四邊形網格，并在壁面附近創建邊界層網格。合適的邊界層網格有助于在數值模擬過程中觀察邊界層形成對激波運動的影響。假設高壓腔和低壓腔的工作流體為理想氣體；采用非定常的雷諾平均納維-斯托克斯方程來分析計算流體運動特性，湍流模型選用SST k-ω模型，粘度模型采用包含流體粘度隨溫度變化的薩瑟蘭粘度模型；采用AUSM格式作為通量格式，時間離散選用二階隱性格式，空間離散采用二階迎風格式。

2.3邊界條件

在不同的隔膜壓力比和激波管直徑下，采用數值模擬的方法分析了微型激波管內部流動。邊界條件如表1所示，高壓腔和低壓腔的溫度都初始化為300 K。模型網格數量為163000，非定常迭代計算的時間步長由最小的計算網格尺寸決定，設定為10-8s。隔膜的初始邊界條件為壁面，在改變其邊界條件為交互面的過程中，這一過程可視為隔膜瞬間破裂的過程。管壁采用絕熱無滑移壁面，并且管壁溫度保持300K不變。因此，在數值模擬過程中，工作流體與激波管管壁之間不會產生熱傳遞。數值模擬編號表示在不同邊界條件下進行的數值模擬，本文采用四種不同的邊界條件來模擬微型激波管內部流體的運動特性。

表1　不同數值模擬條件下的參數值

2.4方法驗證

為了驗證數值模擬方法的正確性，將數值模擬的結果與實驗數據進行對比。文獻［12］采用直徑為6mm的激波管，在隔膜壓力比為6、低壓腔初始壓強為大氣壓強的條件下進行試驗。本文數值模擬中采用相同的邊界條件，結果對比曲線如圖2所示。EH代表膨脹波頭部，SW代表激波。橫坐標代表激波和膨脹波頭部運動的位移，X/D表示激波和膨脹波頭部沿軸向運動的無因次位移，縱坐標表示運動時間。數值模擬結果與試驗結果存在一定的偏差，這種偏差是由在試驗過程中實際流體的粘性以及管內高溫氣體與激波管壁面之間的熱傳遞導致的。在試驗過程中，激波會加熱其運動過后區域內的氣體，高溫氣體會與激波管管壁之間產生熱傳遞，但在數值模擬中，這種熱傳遞不會發生，故激波和膨脹波在試驗過程中產生更多的強度損失。在運動的激波后面會產生湍流邊界層，這也是其強度產生衰減的一個主要原因。以上原因說明數值模擬結果與實驗結果之間的偏差是合理的，本文采用的數值方法是可靠的。因此，目前的數值模擬方法能用來預測微型激波管內激波運動特性。

圖2　數值模擬與試驗結果隨時間變化的對比曲線

2　理論分析

2.1激波管理論

當微型激波管內的薄膜破裂，會產生激波、交接面以及膨脹波。由無粘激波管理論可知，微型激波管內流體粘度可以忽略不計；工作流體為理想氣體；氣體絕熱指數為常數（γ=1.4）。根據高壓腔與低壓腔的初始壓力之比，激波和交接面向低壓腔運動的馬赫數可用式（1）和式（2）計算得到［11］。

其中：MS為激波運動馬赫數，MC為交接面運動馬赫數；γ1和γ4分別表示低壓腔和高壓腔內氣體的絕熱指數；a1和a4分別代表低壓腔和高壓腔氣體聲速大小；P2是激波與交接面之間區域的壓強值。

式（1）和（2）表明：隨著隔膜壓力比的增大，激波和交接面運動的馬赫數也逐漸增大。但由于流體的粘性影響和激波管內存在的稀疏影響，數值模擬得到的結果與理論預測的結果往往存在一定的偏差。

2.2尺寸影響參數S

Brouillette［11］發現，包含熱傳導和剪切力等擴散傳遞現象是導致數值模擬流體特性與理論預測有明顯偏差的主要原因。他提出了控制體分析法，控制體是一個取自激波與交接面之間區域的單元體，通過對控制體內流體特性進行分析，以說明激波管內部的尺寸影響和擴散傳遞現象。控制體的選取如圖3所示，EW表示膨脹波，CS表示交接面，SW表示激波。根據控制體分析法，表示尺寸影響的參數S取自于式（3），參數S的具體表達如式（4）所示。由式（3）可知，較小的S值會對計算激波前后的密度比會產生更大的影響。激波與交接面之間區域雷諾數Re和兩者之間的距離L是計算S必需的兩個參數，隨著Re的減小和L的增大，S值會減?。蝗绻鸖變得無限大，那么可以忽略尺寸影響。

圖3　用于分析激波管內擴散影響的控制體模型

其中：ρ1與ρ2分別表示激波前后的流體密度；Re和Pr分別代表雷諾數和普朗特數；D代表激波管直徑；L表示激波與交接面之間的距離；u2和μ2分別表示激波與交接面之間區域流體的速度和動力粘度。Re由式（5）計算得到。理論上，激波與交接面之間流體速度是相同的，但由于邊界層的存在，靠近壁面的流體速度相對較小，在本文Re的計算中，取交接面的運動速度為u2。

2.3激波前后的壓力梯度

隨著激波在激波管內運動，其前后的壓力梯度也發生變化。當隔膜破裂0.12ms后，激波前后的局部壓力分布如圖4所示。激波前后的壓力梯度與其前后的壓強差以及壓強變化的距離有關，如式（6）所示。式（7）中C表示激波前后壓力梯度的無量綱值，稱其為標準化壓力梯度（normalized pressure gradient，NPG）。由于激波強度隨時間變化逐漸衰弱，激波前后壓強差會逐漸減小，故其前后壓力梯度也不斷減小。

圖4　微型激波管內沿軸線方向的局部壓強變化

3　結果與分析

如表1所示，在隔膜壓力比分別為10、20和100的條件下，探究直徑為3 mm的激波管內激波運動的特征。當膨脹波頭部沒有運動到驅動腔的進口壁面時，激波管內流體沿軸線方向的壓力變化曲線如圖5所示。圖5（a）-（c）分別表示在DPR為10、20和100時，微型激波管流體沿軸線方向的壓力變化。縱坐標表示靜壓比，靜壓比是流體靜壓與低壓腔內流體的初始靜壓之比；橫坐標表示軸向位移。隨著激波向前運動，受激波管內流體粘度以及激波前端與壁面之間的摩擦的影響，其強度逐漸減弱。當DPR為100時，激波在同一時刻運動的位移比DPR為10和20時要大，這說明，隨著隔膜壓力比的增大，激波運動的速度也會增大，這與式（1）的預測相吻合。

圖5　不同隔膜壓力比下激波管內沿軸線方向的壓力變化曲線

激波前后的無量綱壓力梯度由式（7）計算得到，標準化壓力梯度的變化曲線如圖6所示。隨著激波向低壓腔運動，NPG逐漸減小，這主要是由于激波在運動過程中，其強度不斷減弱，導致其前后的壓強差逐漸減小。激波運動的初始階段，NPG衰減得更快，這是由于在初始階段激波前后壓強變化的距離|X2-X1|顯著增大。當激波運動一段位移后（X/D＞30），其前后壓力梯度變化逐漸減小，而這種變化可以視作線性變化。

圖6　不同隔膜壓力比下激波前后壓力梯度的變化曲線

在不同DPR的條件下，由式（5）計算得到的激波與交接面之間區域的雷諾數變化如圖7所示。在這個區域內，隨著激波和交接面向前運動，流體密度、速度與動力粘度發生改變。在三種不同DPR的條件下，Re都逐漸增大，這是由于激波與交接面之間區域流體速度逐漸增大，流體密度也逐漸增大。由于該區域流體溫度會升高，動力粘度也會變大，但流體粘度變大的趨勢沒有速度和密度變大的趨勢明顯，故Re總體變化還是逐漸增大的。當DPR為100時，Re比DPR為10和20時要小。這是由于在大DPR下，P4保持大氣壓強不變，較小的P1導致較小的流體密度，而較小的密度導致較小的雷諾數。

圖7　不同隔膜壓力比下激波與交接面之間區域雷諾數的變化曲線

在不同DPR的條件下，由式（4）計算得到的S值變化曲線如圖8所示。在三種情況下，S值都逐漸減小，這是由于激波和交接面之間存在速度差，導致兩者之間的距離L逐漸增大。當DPR為10時，S值比DPR為20和100時要大。當DPR較小時，雷諾數較大，激波和交接面之間的速度差較小，這導致L也較小。較小的S值能夠表明微型激波管內的尺寸影響，而這種尺寸影響會使激波強度衰減更快。在被驅動腔壓強較小的情況下，計算得到較小的S值，這說明S值同樣也可以用來衡量激波管內低壓對激波運動的影響。

圖8　不同的隔膜壓力比下S值的變化曲線

在不同的DPR下，當隔膜破裂0.03ms后，交接面前后流體的溫度分布以及交接面的位移如圖9所示。等溫線圖能夠準確地反映流體內部的邊界層厚度。如圖9所示，BL表示邊界層。邊界層是由于流體粘度而產生的；是流體在靠近壁面形成的薄流層。在邊界層內，靠近壁面的流體速度為0，沿遠離壁面的垂直方向，流體速度逐漸增大，并達到一恒定值。如圖9所示，在溫度云圖中，邊界層的厚度定義為從壁面到流體溫度達到恒定值的距離。當DPR較大時，交接面運動的位移也相對較大，即交接面運動的速度較大，這與式（2）的預測是相吻合的。與DPR為10和20相比，當DPR為100時，激波管內邊界層的厚度更大，這是因為在DPR較大時，P1較小，由低壓導致的稀疏影響使微型激波管內邊界層變厚。厚度較大的邊界層使激波強度衰減更快，同時使激波后面的流體產生更多的能量耗散和粘性損失。

圖9　不同隔膜壓力比下交接面前后的溫度分布

在不同直徑的激波管內，激波與交接面的速度變化曲線如圖10所示。激波的運動速度在兩種激波管內都逐漸減小，這是由激波前端與激波管壁面之間的摩擦、流體內部粘性以及激波后面邊界層的形成導致的。在小直徑的激波管內，激波速度衰減得更快，這是由微型激波管內的小尺寸影響導致的。交接面的運動速度在兩種激波管內都逐漸增大。在微型激波管內，邊界層隨著激波的運動在其后面逐漸形成，當交接面以馬赫數大于1的速度在激波管內部運動時，邊界層使交接面在一個超聲速噴管內運動，故速度會增加。在小直徑的激波管內，交接面的運動速度增加得更快，這是由于小直徑激波管內形成厚度更大的邊界層。

圖10　不同激波管直徑下激波與交接面速度的變化曲線

在不同直徑的激波管內，S值的變化曲線如圖11所示。在兩種激波管內，S值都隨著時間變化而減小，這是由于隨著激波與交接面向前運動，L會逐漸增大。由式（3）可知，較小的S值對計算激波前后的密度之比是有更明顯影響的。在直徑較小的激波管內，S值明顯較小，激波強度損失也更大，這與預測結果是一致的。這說明用S值來衡量微型激波管內的尺寸影響是可行的。

圖11　不同激波管直徑下S值的變化曲線

4　結 論

在不同隔膜壓力比和激波管直徑的條件下，采用數值模擬的方法來研究微型激波管內的激波、交接面以及流體運動特性，結果表明：隨著隔膜壓力比增大，激波強度逐漸增大，激波在激波管內的運動速度也逐漸增大，交接面的運動速度也逐漸增大；激波在激波管內運動時，其強度逐漸減弱；在被驅動腔壓強較低的情況下，觀察到厚度較大的邊界層，這說明低壓影響對微型激波管內的流體以及激波運動會產生更多能量損失；隨著激波在激波管內運動，其強度逐漸衰弱以及其兩邊壓強變化的距離逐漸變大，這使激波前后的壓力梯度隨著時間變化逐漸減?。唤唤用嬖谖⑿图げü軆冗\動時，其運動速度逐漸增大，這是由運動激波后面邊界層的形成導致的；激波在低壓和小直徑的激波管內運動時，激波強度衰減得更快，并且計算得到較小的S值，這表明S值可以反映微型激波管內的低壓和小尺寸影響。

［1］Duff R E.Shock tube performance at initial low pressure［J］.Phys Fluids，1959，2（1）：207-216.

［2］Ngomo D，Chaudhuri A，Chinnayya A，et al.Numerical study of shock propagation and attenuation in narrow tubes including friction and heat losses［J］.Computers &Fluids，2010，39（9）：1711-1721.

［3］Xiao H，Toshiyuki A，Naoya T.The feature of weak shock wave propagated in an overlong tunnel［J］.Open Journal of Fluid Dynamics，2012，2（4）：285-289.

［4］Mirshekari G，Brouillette M.One-dimensional model for microscale shock tube flow［J］.Shock Waves，2009，19（1）：25-38.

［5］Watvisave D S，Bhandarkar U V，Puranik B P.Effects of wall conduction and rarefaction on shock propagation in a micro-channel［J］.Shock Waves，2014，24（3）：295-306.

［6］Huynh D.The shock tube problem from a combined experimental and computational perspective［C］∥21th AIAA Computational Fluid Dynamics Conference.San Diego，California，USA，2013，3：2442-2458.

［7］Sturtevant B，Okamura T T.Dependence of shock tube boundary layers on shock strength［J］.Phys Fluids，1969，12（8）：1723-1725.

［8］Tanaki K，Inaba K，Yamamoto M.Numerical investigation on transition of shock induced boundary layer［C］∥47th AIAA Aerospace Science Meeting Including the New Horizons Forum and Aerospace Exposition.Orlando，Florida，2009：5-8.

［9］Mirels H.Test time in low pressure shock tubes［J］. Phys Fluids，1963，6（9）：1201-1214.

［10］Roshko A.On flow duration in low pressure shock tubes［J］.Phys Fluids，1960，3（6）：835-842.

［11］Brouillete M.Shock waves at microscales［J］.Shock Waves，2003，13（1）：3-12.

［12］Park JO，Kim G W，Kim H D.Experimental study of the shock wave dynamics in micro shock tube［J］. Journal of the Korean Society of Propulsion Engineers，2014，17（5）：54-59.

［13］Sun M，Ogawa T，Takayama K.Shock propagation in narrow channels［C］∥Processing of 24th International Symposium on Shock Waves.Tohoku University，Katahira，Japan，2001：1321-1327.

［14］Zeitoun D E.Microsize and initial pressure effects on shock wave propagation in a tube［J］.Shock Waves，2014，24（5）：515-520.

［15］Zeitoun D E，Burtschell Y.Navier-Stokes computations in micro shock tubes［J］.Shock Waves，2006，15（3）：241-246.

［16］Arun K R，Kim H D.Computational study of the unsteady flow characteristics of a micro shock tube［J］.Journal of Mechanical Science and Technology，2013，27（2）：451-459.

［17］Arun K R，Kim H D.Numerical visualization of the unsteady shock wave flow field in micro shock tube［J］.Journal of the Korean Society of Visualization，2012，10（1）：40-46.

Numerical Simulation of Shock Wave Characteristics in Micro Shock Tubes

ZHANGGuang1，CUI Bao-ling1，JIN Ying-zi1，KIM Heuy-dong2
（1.The Zhejiang Provincial Key Laboratory of Fluid Transmission Technology Research，Zhejiang Sci-Tech University，Hangzhou 310018，China；2.Department of Mechanical Engineering，Andong National University，Andong 760749，Republic of Korea）

Numerical simulation was applied to study internal flow of micro shock tube in order to study unsteady flow and shock wave motion characteristics in micro shock tube.Contrastive analysis was conducted for effects of different diaphragm pressure ratio（initial pressure ratio of high-pressure cavity and low-pressure cavity），and shock tube diameter on shock wave in micro shock tube，contact surface and fluid motion characteristics.Besides，experimental data were compared.The results show motion speed of shock wave and contact surface increases gradually with the increase in the diaphragm pressure ratio；when shock wave moves in micro shock tube，shock wave strength gradually weakens；a thick boundary layer is observed when pressure intensity in low-pressure driven cavity，which indicates that low pressure effect makes shock wave and flow more energy lose；pressure gradient in the front of and after shock wave gradually decreases as shock wave motion decreases gradually；the motion speed gradually increases when the contact surface moves in micro shock tube；in addition，S value can indicate effects of low pressure and small scale in micro shock tube.

micro shock tube；shock wave；contact surface；unsteady flow；numerical simulation

TH47

A

1673-3851（2015）06-0805-07

（責任編輯：康 鋒）

2015-01-13

國家自然科學基金項目（51406184）

張光（1988-），男，湖北孝昌人，碩士研究生，主要從事微型激波管方面的研究。

崔寶玲，E-mail：blcui@zstu.edu.cn