李小明，許志宇,2，譚永華，胡 攀

(1.西安航天動力研究所，陜西 西安 710100; 2.液體火箭發(fā)動機技術重點實驗室，陜西 西安 710100;3.航天推進技術研究院，陜西 西安 710100)

0 引言

電爆管是電爆閥、拔銷器等裝置的驅動能量源，具有能量密度高、體積小和工作可靠等優(yōu)點。電爆管起爆后產生高溫高壓燃氣，在傳火通道中形成強烈沖擊波，沖擊波沿傳火通道傳播和驅動做功的同時，伴隨著接觸間斷、稀疏波的傳播、反射及其相互作用等復雜過程[1-2]。電爆管的壓力輸出特性決定著電爆閥、拔銷器等裝置的工作特性和可靠性。但是，由于爆腔和傳火通道十分狹小，壓力和溫度極高，很難通過試驗準確獲取傳火通道中壓力分布和變化規(guī)律。目前利用容積數(shù)倍于真實爆腔的標準密閉爆發(fā)器測試電爆管的壓力曲線，并利用峰值壓力通過靜態(tài)等溫過程換算真實容腔的壓力和驅動力[3]。由于忽略了沖擊波傳播的動態(tài)過程和空間分布特性，估算誤差可能較大，而且不能評估傳火通道結構對沖擊波傳播規(guī)律的影響。因此建立考慮空間分布的一維模型，分析電爆管沖擊波傳播過程和傳播規(guī)律。

傳統(tǒng)的沖擊波計算方法如ENO，WENO，MUSCL及AUSM等為了獲得高分辨率、高精度解，需要非常細密的網格，以便能捕捉清晰的沖擊波結構。但是，精度和分辨率越高，計算量越大，特別是高維問題，除非使用非規(guī)則網格，否則因計算量巨大的原因將很難實現(xiàn)[4]。小波數(shù)值方法是基于多分辨分析(MultiResolution Analysis, MRA)發(fā)展的新方法，由于小波函數(shù)具有緊支撐特性，因此能夠對流場數(shù)據進行壓縮，生成捕捉流場局部結構的適應性網格，適合描述局部流動特征顯著的問題[4-6]，但仍處于發(fā)展階段，并未廣泛應用。

小波數(shù)值方法主要有兩類，即小波-迦遼金和小波配點法[5]。小波-迦遼金法不適合處理非線性算子和任意邊界條件，而自適應小波配點法在這兩方面均具有優(yōu)勢，特別是二代小波，在真實物理域中進行變換，可以方便處理任意邊界條件，因此發(fā)展迅速。文獻[6]利用二代小波，采用自適應配點法計算了Burgers方程和含化學反應的無黏Euler方程；文獻[7]利用小波配點方法模擬了地震波在各向異性介質中的傳播；文獻[8-14]研究了小波數(shù)值方法在激波、湍流、燃燒、爆轟波、水擊等方面的應用。研究表明，小波數(shù)值方法在捕捉水擊波、激波波、爆轟波以及火焰局部結構等方面具有顯著優(yōu)勢，分辨率和精度高，計算量小。

本文基于一維Euler方程，采用自適應小波配點法(Adaptive Wavelet Collocation Methods,AWCM)對電爆管沖擊波的傳播過程進行數(shù)值計算。詳細介紹自適應二代小波配點方法的原理、沖擊波捕捉方法、計算過程，初步分析電爆管沖擊波沿傳火通道的傳播規(guī)律。

1 自適應小波配點法

1.1 多分辨分析和二代小波變換

根據多分辨分析理論可知，任意函數(shù)f(x)∈L2(R)的多尺度分解可近似為

(1)

式中：sj0和φj0分別為尺度系數(shù)和尺度函數(shù)；dj和ψj分別為小波系數(shù)和小波函數(shù)；J為尺度因子。

(2)

分解過程分為3步：分裂(Split)、預測(Predict)和更新(Update)，如圖1所示。

圖1 二代小波變換示意圖Fig.1 Forward and inverse transformation of second-generation wavelet

分裂：按尺度因子J進行規(guī)則采樣S={ck,k=0,…,2J}，然后將S按順序號奇偶分為Se和So兩部分。

預測：利用偶序列中相鄰的2N個點預測奇序列，設P為預測算子，則有

(3)

(4)

逆變換過程為

(5)

(6)

1.2 自適應網格

對于局部包含小尺度特征的函數(shù)，僅在相應的局部位置小波系數(shù)較大，而大部分小波系數(shù)將會很小，即便舍去這些系數(shù)很小的小波函數(shù)，也能很好地逼近原函數(shù)?；谶@種思想，人工選定或通過某種規(guī)則計算得到一個閾值ε(ε>0)，根據小波系數(shù)相對于閾值的大小分為兩部分

f(x)=f≥(x)+f<(x)

(7)

其中

(8)

(9)

由于每個小波函數(shù)(也包括尺度函數(shù))和配點一一對應，因此舍去小波的同時刪除了對應的配點，從而生成自適應網格。對于正則方程，舍去式(9)中系數(shù)較小的小波，逼近的誤差上限為[6]

‖f(x)-f≥(x)‖≥C1ε‖f(x)‖

(10)

1.3 沖擊波捕捉技術

添加人工黏性項的一維Euler方程

(11)

式中：U為守恒變量；F為通量。等號右側為人工黏性項，人工黏性υ為沖擊波定位函數(shù)Φ(或稱通量限制器)的函數(shù)，而Φ隨小波系數(shù)的大小變化。最細尺度的小波系數(shù)能反應沖擊波位置，而對于j

(12)

Euler方程組包含多個守恒變量，理論上需要對每個變量進行小波變換，并分別計算沖擊波定位函數(shù)?？紤]到沖擊波問題和一般黎曼問題的共性，選擇守恒變量ρ進行小波變換、生成適應性網格以及構造沖擊波定位函數(shù)，有利于簡化計算和節(jié)約計算量。但對于強沖擊波問題，可能會由于人工黏性不足，并不能有效抑制數(shù)值振蕩?？紤]到Φ在[0,1]區(qū)間內變化，采用冪函數(shù)形式的定位函數(shù)能夠控制黏性分布的寬度

(13)

對于強沖擊波問題，一般指數(shù)因子α<1，以便有效抑制數(shù)值振蕩。利用沖擊波定位函數(shù)控制的人工粘性[4]

(14)

式中：c為當?shù)芈曀?，c=(γp/ρ)1/2；u為速度。

通量的空間導數(shù)采用二階中心差分格式，守恒量的時間導數(shù)采用一階步進格式，則方程(11)的離散格式為

(15)

具體步驟為：

1)由t時刻的值通過二代小波變換獲得各尺度的小波系數(shù)；

2)根據小波系數(shù)相對于閾值ε的大小刪除相應的小波函數(shù)，從而生成自適應網格；

3)在自適應網格上計算空間導數(shù)；

4)選擇時間步長Δt，計算(t+Δt)時刻的值，并返回1)繼續(xù)計算。

這種算法可以適當調整網格，尺度因子越大，則分辨率越高；選擇的閾值越小，精度越高。

2 電爆管傳火通道模型

電爆管一般屬于標準火工品元件，匹配不同的閥門或者拔銷器的結構，可能在連接部位形成變截面的傳火通道。根據不同的匹配形式，將電爆管傳火通道簡化為4類一維流道，作為不同結構型面的近似，即等截面、錐形擴張、錐形收縮以及拉瓦爾噴管型結構。局部截面積線性增大或減小，截面變化的過渡距離與爆腔長度相同，因此變截面特征可以僅通過一個參數(shù)——截面面積比

(16)

式中：A0為爆腔橫截面積，當為拉瓦爾噴管結構時，傳火通道與爆腔橫截面積相同;A為喉部截面積，當為錐形結構時，A為擴張或收縮后的傳火通道截面積；η<1為收縮型；η>1為擴張型。如圖2所示，截面變化的過渡距離與爆腔長度L均為10 mm。特別地，η=1為理想的等截面模型，對于變截面通道，定義Δη=|1-η|，用于表征相對于等截面通道的變化程度。

圖2 傳火通道橫截面積變化示意圖Fig.2 Cross-sectional area for different passage types

假設炸藥在電爆管爆腔內瞬時完成定容反應，并完全達到化學平衡，初始條件如圖3所示，根據爆炸腔的容積，炸藥及其大量爆炸產生物的性質計算初始參數(shù)如表1所示。狀態(tài)方程采用理想氣體狀態(tài)方程，燃氣比熱比γ=1.25。電爆管封閉端設為剛性壁面條件，出口可使用壓力出口或壓力遠場條件，在沖擊波到達出口前，均不影響上游沖擊波的計算。

圖3 一維等截面?zhèn)骰鹜ǖ滥Ｐ虵ig.3 One-dimensional physical model of the cartridge with an equal section passage

表1 電爆管模型初始條件

Tab.1 Initial conditions of the cartridge model

計算域壓力/ MPa密度/(kg·m-3)Ⅰ375.36254.65Ⅱ0.101 3251.00

3 計算結果和分析

3.1 模型和方法驗證

為了驗證計算方法在強沖擊波條件下的準確性，利用可求得精確解的沖擊波管問題進行驗證，計算域為[-5,5]，初始條件與電爆管模型完全相同，即表1所示?；A網格點為1 025個(J=10)，濾波閾值ε=2×10-5。圖4所示為5 μs和10 μs時刻速度分布，可以看出自適應小波配點法(AWCM)的計算結果與黎曼解符合得很好，并且無明顯的數(shù)值振蕩，表明沖擊波的小波捕捉方法能夠準確計算強沖擊波問題，適用于電爆管沖擊波傳播特性的預測。采用自適應算法，刪除了變化平緩區(qū)域的大量網格點，使網格點主要集中在沖擊波波陣面、稀疏波以及接觸間斷區(qū)域。

3.2 等截面通道沖擊波傳播規(guī)律

取爆腔長度L=10 mm、傳火通道長度l=10L，電爆管沖擊波傳播過程的計算參數(shù)與沖擊波管相同。計算的沖擊波傳播速度為3.62 km/s，經過約27.6 μs傳播至傳火通道的右端出口，與黎曼解完全一致，表明計算準確。圖5為傳火通道中不同位置的壓力歷程和不同時刻壓力分布，可以看出，爆腔內壓力衰減迅速，10 μs時，爆腔內最大壓力衰減為200 MPa，20 μs時，最大壓力約為40 MPa；壓力隨傳播距離迅速衰減，距離爆腔越近，峰值壓力衰減率越大，如圖6所示。

圖4 沖擊波管的黎曼解與小波數(shù)值解Fig.4 Riemann and wavelet numerical solutions of the shock tube

對于等截面的傳火通道，不同位置的峰值壓力p與爆腔長度L、初始壓力p0和傳播距離x有關，構建無量綱的傳播距離x/L和峰值壓力pm/p0，分兩個區(qū)間擬合指數(shù)形式的無量綱函數(shù)，結果如式(17)所示，其中x/L>1時的函數(shù)曲線如圖6所示。對于同一類電爆管，爆炸沖擊波傳播相似律如下

(17)

圖5 等截面?zhèn)骰鹜ǖ乐袎毫ψ兓头植糉ig.5 Pressure variation and pressure distribution in the equal-sectional passage

圖6 等截面?zhèn)骰鹜ǖ啦煌恢玫撵o壓峰值Fig.6 Static pressure peaks at different locations of equal-section passage

為了說明相似律表達式的適用性，選擇規(guī)格不同的同類電爆管進行驗證。爆腔初始壓力設為p0=187.68 MPa，電爆管爆腔長L=5 mm，通過數(shù)值計算和式(17)預估的壓力如表2所示。不同位置的估算誤差一般小于5%，表明式(17)能夠準確預估同一類電爆管爆炸沖擊波在等截面?zhèn)骰鹜ǖ纼鹊姆逯祲毫Α?/p>

表2 等截面通道峰值壓力估算驗證

3.3 變截面通道沖擊波傳播規(guī)律

在相同的初始條件(表1)下，計算局部變截面?zhèn)骰鹜ǖ缐毫υ?～25 μs內的發(fā)展變化，結果如圖7所示，其中η=1表示等截面?zhèn)骰鹜ǖ?，作為不同結構比較的基準。

對于拉瓦爾噴管型結構，喉部面積減小導致輸出壓力降低，且對爆炸近區(qū)的影響更顯著，η=0.75時，x/L=1和2處的峰值壓力分別降低18%和9%。而對于純收縮或擴張型結構，傳火通道截面積越小，輸出壓力越高，在x/L=2的位置，η=0.75和0.5時，峰值壓力對應增大27%和78%，而面積增加25%和50%時，峰值壓力相應減小15%和28%?？梢钥闯觯谙嗤某跏紬l件下，變截面?zhèn)骰鹜ǖ纼炔煌恢玫姆逯祲毫Σ粌H取決于傳播距離x,而且受傳火通道結構形式和參數(shù)η的影響，傳播規(guī)律更為復雜。構建無量綱的峰值壓力pm/p0、傳播距離x/L以及局部截面面積比η，利用數(shù)值計算的結果擬合得到不同位置無量綱的峰值壓力pm/p0與傳播距離x/L和局部截面積比η之間的關系如下

(18)

(19)

利用3.2節(jié)中假設的電爆管模型(p0=187.68 MPa，L=5 mm)進行驗證，不同位置峰值壓力和估算誤差如表3和表4所示。當1

圖7 不同結構傳火通道的壓力特性Fig.7 Pressure characteristics of different passages

表3 局部拉瓦爾結構傳火通道峰值壓力估算驗證

Tab.3Verification of the estimated peak pressures forthe passages with a Laval structure

η位置/mm峰值壓力/MPa數(shù)值解相似律誤差/%0.6620.521.96.81116.214.7-9.21612.111.1-8.3219.68.9-7.30.8625.126.24.41118.116.7-7.71613.412.3-8.22110.19.6-5.0

表4 局部收縮或擴張傳火通道峰值壓力估算驗證

4 結論

基于二代小波配點法和人工黏性技術，構造了沖擊波的小波數(shù)值計算方法，并應用于電爆管一維沖擊波數(shù)值計算，根據計算結果初步分析了電爆管爆炸沖擊波的傳播規(guī)律，分析表明：

1)自適應小波配點法結合人工黏性技術可以準確計算電爆管強沖擊波傳播問題。

2)電爆管爆炸沖擊波沿一維等截面?zhèn)骰鹜ǖ纻鞑M足一定相似律，建立的相似律表達式能夠準確估算同一類電爆管沖擊波在爆炸近區(qū)的峰值壓力。

3)傳火通道結構形式和參數(shù)變化對沖擊波傳播規(guī)律影響顯著，針對同一類型的電爆管和傳火通道建立的相似律表達式在1