基于二代小波的電爆管沖擊波數(shù)值計(jì)算

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

(1.西安航天動(dòng)力研究所，陜西 西安 710100; 2.液體火箭發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710100;3.航天推進(jìn)技術(shù)研究院，陜西 西安 710100)

0 引言

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

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

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

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

1 自適應(yīng)小波配點(diǎn)法

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

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

(1)

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

(2)

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

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

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

預(yù)測(cè)：利用偶序列中相鄰的2N個(gè)點(diǎn)預(yù)測(cè)奇序列，設(shè)P為預(yù)測(cè)算子，則有

(3)

(4)

逆變換過(guò)程為

(5)

(6)

1.2 自適應(yīng)網(wǎng)格

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

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

(7)

其中

(8)

(9)

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

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

(10)

1.3 沖擊波捕捉技術(shù)

添加人工黏性項(xiàng)的一維Euler方程

(11)

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

(12)

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

(13)

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

(14)

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

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

(15)

具體步驟為：

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

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

3)在自適應(yīng)網(wǎng)格上計(jì)算空間導(dǎo)數(shù)；

4)選擇時(shí)間步長(zhǎng)Δt，計(jì)算(t+Δt)時(shí)刻的值，并返回1)繼續(xù)計(jì)算。

這種算法可以適當(dāng)調(diào)整網(wǎng)格，尺度因子越大，則分辨率越高；選擇的閾值越小，精度越高。

2 電爆管傳火通道模型

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

(16)

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

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

假設(shè)炸藥在電爆管爆腔內(nèi)瞬時(shí)完成定容反應(yīng)，并完全達(dá)到化學(xué)平衡，初始條件如圖3所示，根據(jù)爆炸腔的容積，炸藥及其大量爆炸產(chǎn)生物的性質(zhì)計(jì)算初始參數(shù)如表1所示。狀態(tài)方程采用理想氣體狀態(tài)方程，燃?xì)獗葻岜圈?1.25。電爆管封閉端設(shè)為剛性壁面條件，出口可使用壓力出口或壓力遠(yuǎn)場(chǎng)條件，在沖擊波到達(dá)出口前，均不影響上游沖擊波的計(jì)算。

圖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

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

3 計(jì)算結(jié)果和分析

3.1 模型和方法驗(yàn)證

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

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

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

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

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

(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

為了說(shuō)明相似律表達(dá)式的適用性，選擇規(guī)格不同的同類電爆管進(jìn)行驗(yàn)證。爆腔初始?jí)毫υO(shè)為p0=187.68 MPa，電爆管爆腔長(zhǎng)L=5 mm，通過(guò)數(shù)值計(jì)算和式(17)預(yù)估的壓力如表2所示。不同位置的估算誤差一般小于5%，表明式(17)能夠準(zhǔn)確預(yù)估同一類電爆管爆炸沖擊波在等截面?zhèn)骰鹜ǖ纼?nèi)的峰值壓力。

表2 等截面通道峰值壓力估算驗(yàn)證

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

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

對(duì)于拉瓦爾噴管型結(jié)構(gòu)，喉部面積減小導(dǎo)致輸出壓力降低，且對(duì)爆炸近區(qū)的影響更顯著，η=0.75時(shí)，x/L=1和2處的峰值壓力分別降低18%和9%。而對(duì)于純收縮或擴(kuò)張型結(jié)構(gòu)，傳火通道截面積越小，輸出壓力越高，在x/L=2的位置，η=0.75和0.5時(shí)，峰值壓力對(duì)應(yīng)增大27%和78%，而面積增加25%和50%時(shí)，峰值壓力相應(yīng)減小15%和28%?？梢钥闯?，在相同的初始條件下，變截面?zhèn)骰鹜ǖ纼?nèi)不同位置的峰值壓力不僅取決于傳播距離x,而且受傳火通道結(jié)構(gòu)形式和參數(shù)η的影響，傳播規(guī)律更為復(fù)雜。構(gòu)建無(wú)量綱的峰值壓力pm/p0、傳播距離x/L以及局部截面面積比η，利用數(shù)值計(jì)算的結(jié)果擬合得到不同位置無(wú)量綱的峰值壓力pm/p0與傳播距離x/L和局部截面積比η之間的關(guān)系如下

(18)

(19)

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

圖7 不同結(jié)構(gòu)傳火通道的壓力特性Fig.7 Pressure characteristics of different passages

表3 局部拉瓦爾結(jié)構(gòu)傳火通道峰值壓力估算驗(yàn)證

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 局部收縮或擴(kuò)張傳火通道峰值壓力估算驗(yàn)證

4 結(jié)論

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

1)自適應(yīng)小波配點(diǎn)法結(jié)合人工黏性技術(shù)可以準(zhǔn)確計(jì)算電爆管強(qiáng)沖擊波傳播問(wèn)題。

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

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