基于通用炸藥狀態(tài)方程分析飛板運(yùn)動(dòng)規(guī)律的特征線法＊

李曉杰，趙春風(fēng)，2

（1.大連理工大學(xué)運(yùn)載工程與力學(xué)學(xué)部工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.大連理工大學(xué)工程抗震研究所，遼寧 大連 116024）

由于在滑移爆轟驅(qū)動(dòng)下飛板的運(yùn)動(dòng)狀態(tài)直接影響爆炸焊接工藝中各種參數(shù)的選擇和焊接質(zhì)量，所以研究飛板的運(yùn)動(dòng)規(guī)律是爆炸焊接理論研究中很重要的內(nèi)容。對(duì)于穩(wěn)定的滑移驅(qū)動(dòng)飛板運(yùn)動(dòng)，將坐標(biāo)系置于爆轟波頭上，爆轟產(chǎn)物的流場(chǎng)就是定常的，可按二維可壓縮定常流處理爆轟產(chǎn)物的流動(dòng)，因此，可以采用特征線差分法對(duì)爆轟產(chǎn)物流場(chǎng)進(jìn)行數(shù)值計(jì)算。以往在對(duì)這種定常爆轟流場(chǎng)進(jìn)行計(jì)算時(shí)，一直用多方等熵方程來(lái)描述爆轟產(chǎn)物［1－5］，沒有使用JWL和簡(jiǎn)化JWL等描述炸藥爆轟的專用狀態(tài)方程。這不僅影響爆轟流場(chǎng)的計(jì)算精確度，而且影響對(duì)爆轟驅(qū)動(dòng)力學(xué)實(shí)質(zhì)的分析。為此，本文中力求在使用通用狀態(tài)方程的基礎(chǔ)上，由二維可壓縮定常流的基本理論導(dǎo)出特征線方程，并針對(duì)二維滑移爆轟飛板拋擲問(wèn)題重新編制計(jì)算程序，以常用的JWL狀態(tài)為例進(jìn)行數(shù)值計(jì)算。沿用文獻(xiàn)［1］中的爆轟波后采用的弱爆轟假設(shè)來(lái)處理聲速面問(wèn)題。另外，采用一種直接從物理概念出發(fā)獲得等熵?zé)o旋定常超聲速二維流動(dòng)特征線方程的方法，而沒有使用通常的由控制方程出發(fā)，采用不定線法或方向?qū)?shù)法來(lái)推導(dǎo)特征線方程的方法［6］，以期更明晰地描述特征線的物理意義。

1 二維定常超聲速流的特征線方程

在均勻的超聲速流場(chǎng)中，處于流場(chǎng)中某一固定點(diǎn)的小擾動(dòng)源，所產(chǎn)生的擾動(dòng)如圖1所示，球?qū)ΨQ地以聲速c相對(duì)流體向外傳播；流體本身以速度q運(yùn)動(dòng)，小擾動(dòng)傳播的絕對(duì)速度是聲速c與流速q的矢量迭加。對(duì)于超聲速流，有q＞c，擾動(dòng)面是一個(gè)錐面，且與從擾動(dòng)點(diǎn)出發(fā)的一串球面相切，該錐面為馬赫錐，馬赫錐的邊界面（線）稱為馬赫面（線）。半錐角α為［1，5］

通常把半錐角α稱為馬赫角，Ma稱為馬赫數(shù)。

圖1 超聲速流中的小擾動(dòng)Fig.1Aperturbation in supersonic flow

對(duì)于本文中要研究的平面問(wèn)題，擾動(dòng)源是無(wú)限長(zhǎng)的線源，沿?cái)_動(dòng)線一系列馬赫錐組成一個(gè)楔形。在二維平面上，如圖2所示，馬赫錐組成的楔形退化為2條馬赫線I和II，在Oxy坐標(biāo)系中，這2條馬赫線的方程為［1，5］

式中：θ是流動(dòng)方向角。如果馬赫波后的小擾動(dòng)使流速由q變?yōu)閝＋δq，流動(dòng)方向角由θ變?yōu)棣龋摩?，則流動(dòng)穿過(guò)馬赫波I時(shí)，沿馬赫波的平行方向的速度動(dòng)量方程有

圖2 超聲速流中的特征線Fig.2Characteristic curves in supersonic flow

式中：Qm＝ρqsinα為通過(guò)馬赫波上的質(zhì)量流量。展開式（3），并考慮到δq→dq，δθ→dθ，且cosδθ→1，sinδθ→dθ，并忽略二階小量，將式（1）代入式（3）可得

同理，對(duì)穿過(guò)馬赫線II的參數(shù)變化可以寫出

從數(shù)學(xué)意義上來(lái)講，馬赫線就是特征線，馬赫線I和II分別對(duì)應(yīng)第1族和第2族特征線，式（2）就是這2族特征線的方程，而式（4）和（5）分別是波陣面前、后的相容關(guān)系。任何復(fù)雜的超聲速流場(chǎng)均可以被看成是由一系列連續(xù)的空間小擾動(dòng)相互作用的結(jié)果，所以，常利用2族特征線方程和對(duì)應(yīng)的相容關(guān)系來(lái)求解超聲速流場(chǎng)。在進(jìn)行數(shù)值求解時(shí)，一般使用沿特征線的相容關(guān)系。由于在超聲速膨脹流中同族特征線不相交，一族特征線穿過(guò)另一族特征線，因此，一族特征線波陣面上的相容關(guān)系就是沿另一族特征線上的相容關(guān)系，所以可以獲得沿第1、2族特征線上的相容關(guān)系為［1，5］

從物理概念出發(fā)獲得的平面等熵定常超聲速流動(dòng)的特征線方程（2）和（6），由于方程式不依賴流體物態(tài)方程相關(guān)的參數(shù)，所以式（2）和（6）是通用物態(tài)方程的特征線關(guān)系。

然后，可以利用能量方程和具體的氣體等熵方程來(lái)建立馬赫數(shù)Ma和流速q與流動(dòng)方向θ之間的關(guān)系。對(duì)于定常可壓縮流，用熱焓i和流速q的關(guān)系表示其伯努利方程［5，7］

式中：熱焓i＝e＋p／ρ，e為比內(nèi)能，p為爆壓，ρ為炸藥密度；i0為滯止焓。微分式（7）可得

再利用Ma2＝ （q／c）2＝2 （i0－i）／c2，等熵條件di＝dp／ρ和聲速c2＝（?p／ ?ρ）S＝dp／dρ，由上面推導(dǎo)可以得到下面的關(guān)系

由于在等熵線上，流體的熱焓i、聲速c，包括馬赫數(shù)Ma都是密度ρ的單值函數(shù)，所以式（9）是可以積分的，故可以將式（9）的積分定義成一個(gè)函數(shù)

式中：函數(shù)ν（ρ）相當(dāng)于Prandtl－Meyer（P－M）函數(shù)［6］，ρCJ為炸藥CJ狀態(tài)時(shí)的密度。

將式（10）定義為不依賴流體狀態(tài)方程的P－M函數(shù)，或稱為通用狀態(tài)方程的P－M函數(shù)。在式（6）中引入該函數(shù)寫成新的特征線相容關(guān)系，再與式（2）合并，可以寫成如下超聲速定常流特征線方程組

式中：RⅠ和RⅡ是常數(shù)，分別為第1和第2黎曼不變量。顯然，只要根據(jù)具體的流體等熵狀態(tài)方程求出熱焓i、聲速c代入式（10）求得P－M函數(shù)ν（ρ），就可以和式（11）組成封閉的求解方程組了。如對(duì)于常用的JWL炸藥狀態(tài)方程，其等熵方程如下

式中：pS為等熵壓力，eS為等熵比內(nèi)能，ρ0為初始密度，A、B、C、R1、R2和ω為常數(shù)，V＝ρ0／ρ是相對(duì)比體積。根據(jù)相應(yīng)的定義可以推導(dǎo)出JWL狀態(tài)方程的聲速和熱焓表達(dá)式

由于特征線方程（11）為常微分方程組，所包含的P－M函數(shù)ν（ρ）可以用辛普森積分等高精度積分方法求解，所以使用特征線方程可以獲得高精度的差分解。

2 滑移爆轟驅(qū)動(dòng)飛板的特征線差分法

在確立了平面超聲速流的特征線方程后，可將特征線方程離散化，對(duì)圖3所示的滑移爆轟驅(qū)動(dòng)飛板問(wèn)題進(jìn)行數(shù)值計(jì)算。

2.1 初值選取

如圖3所示，炸藥向空氣飛散在O點(diǎn)處可以近似用中心稀疏波處理。由于在爆轟波CJ面上的馬赫數(shù)MaCJ＝1，對(duì)應(yīng)的密度ρ＝ρCJ，氣流轉(zhuǎn)角θ＝0，根據(jù)普朗特繞流理論的第2族特征線相容關(guān)系式［7］

圖3 滑移爆轟作用下的飛板運(yùn)動(dòng)姿態(tài)Fig.3Movement of flyer plate driven by glancing detonation

可以在ρCJ到某一較小ρn范圍內(nèi)給出一系列的離散點(diǎn)（ρi，θi），其中i＝0，1，2，…，n；離散點(diǎn)的（x，y）坐標(biāo)均在O點(diǎn)處，這些離散點(diǎn)就是差分的初值點(diǎn)。再利用下式確定滯止焓

2.2 聲速面處理

由于在爆轟產(chǎn)物流場(chǎng)中，爆轟波面后的CJ面是聲速面，從O點(diǎn)發(fā)出的馬赫線會(huì)在飛板壁面上垂直反射，導(dǎo)致無(wú)法計(jì)算差分值。具體的處理方法可以簡(jiǎn)單地令炸藥爆轟稍向弱爆轟方向偏離CJ點(diǎn)［5］，只對(duì)整個(gè)流場(chǎng)產(chǎn)生微小的影響。即令＝ （1－Δ）ρCJ，使爆轟波后完全變?yōu)槌曀倭鲌?chǎng)，就完全可使用特征線差分求解。這樣做引起的誤差也很小，把Δ取足夠小量，對(duì)計(jì)算結(jié)果影響不大。在實(shí)際計(jì)算中，Δ可取為10－2～10－3量級(jí)，對(duì)量綱一壓力的擾動(dòng)為10－4～10－6量級(jí)。

2.3 超聲速流場(chǎng)域內(nèi)差分

在爆轟產(chǎn)物的超聲速流場(chǎng)內(nèi)部，如圖4所示，當(dāng)已知點(diǎn)1、2的流動(dòng)參數(shù)時(shí)，可以用2族特征線相交求取出未知點(diǎn)3上的參數(shù)。由式（11）將第1族特征線RⅠ方程由已知點(diǎn)1至未知點(diǎn)3離散成差分形式，第2族特征線RⅡ方程由已知點(diǎn)2至未知點(diǎn)3離散，經(jīng)聯(lián)立整理可得

圖4 爆轟產(chǎn)物中差分示意圖Fig.4Scheme of differences in detonation products

式中：下標(biāo)1～3分別表示點(diǎn)1～3。按式（16）由點(diǎn)1、2的參數(shù)差分解出θ3，進(jìn)而用P－M反函數(shù)求得ρ3。然后，利用狀態(tài)方程，求解最后代入式（17）確定點(diǎn)3的坐標(biāo)。

2.4 飛板邊界差分

飛板在爆轟壓力驅(qū)動(dòng)下的運(yùn)動(dòng)微分方程為［6］

式中：R為炸藥質(zhì)量比，即炸藥的厚度和密度的乘積與飛板的厚度和密度乘積的比值；＝p／pCJ，pCJ為等熵狀態(tài)的爆轟壓力；D為炸藥爆速；pk為中間變量；x、y為藥厚量綱一化的坐標(biāo)。如圖5所示，當(dāng)已知飛板邊界點(diǎn)1和流場(chǎng)域內(nèi)點(diǎn)2時(shí)，可利用飛板運(yùn)動(dòng)方程（18），與通過(guò)點(diǎn)2的第1族特征線求出點(diǎn)3的參數(shù)。將上述邊界方程作如下離散

圖5 飛板邊界差分示意圖Fig.5Scheme of differences on the boundary of flyer plate

對(duì)通過(guò)點(diǎn)2、3的第1族特征線方程也離散成差分形式

將式（19）代入式（20），整理得

預(yù)估方程組（21）可以用預(yù)估－校正法獲得較高的差分精度。具體做法是，以k＊＝cot（θ2＋α2）和pk＝Rp1／（ρ0D2）作為預(yù)估初值代入式（21）解出θ3、ρ3，并用狀態(tài)方程解出對(duì)應(yīng)的點(diǎn)3的i3、c3、Ma3和α3；再用k＊＝ ［cot（θ2＋α2）＋cot（θ3＋α3）］／2和pk＝R（p1＋p3）／（2ρ0D2）代入式（21），重新求解一遍，可獲得經(jīng)過(guò)校正的點(diǎn)3的參數(shù)。最后，用式（19）求出點(diǎn)3的坐標(biāo)。

2.5 實(shí)例計(jì)算

用上述通用狀態(tài)方程的差分格式和方法，編制計(jì)算程序?qū)︼w板二維拋擲問(wèn)題進(jìn)行數(shù)值求解。計(jì)算中TNT炸藥的JWL狀態(tài)方程參數(shù)［8］：多方指數(shù)K，2.727 6；CJ狀態(tài)的密度，1.630g／cm3；CJ狀態(tài)的爆壓，21.0GPa；CJ狀態(tài)的爆速，6.93km／s；A，373.8GPa；B，2.747GPa；C，0.734GPa；R1，4.15；R2，0.90；ω，0.35。乳化炸藥的JWL狀態(tài)方程參數(shù)［9］：多方指數(shù)K，3.008；CJ狀態(tài)的密度，1.145g／cm3；CJ狀態(tài)的爆壓，7.62GPa；CJ狀態(tài)的爆速，5.141 4km／s；A，326.42GPa；B，5.808 9GPa；C，1.032 5GPa；R1，5.8；R2，1.56；ω，0.57。圖6為所計(jì)算的 TNT和乳化炸藥滑移爆轟驅(qū)動(dòng)飛板的拋擲曲線，圖7為飛板的動(dòng)態(tài)彎折角與量綱一飛行距離的關(guān)系，圖8為TNT、乳化炸藥的JWL和多方方程的等熵卸載線。

由圖6～7可以看出，對(duì)于TNT炸藥，用多方狀態(tài)方程描述的對(duì)飛板的加速能力較用JWL方程描述的大；取量綱一距離飛行距離y＝1處的彎折角進(jìn)行對(duì)比，分別為11.96°和10.57°，用JWL狀態(tài)方程的計(jì)算結(jié)果較用多方方程的計(jì)算結(jié)果小11.62%。對(duì)于乳化炸藥，用多方狀態(tài)方程描述的對(duì)飛板的加速能力較用JWL方程描述的??；取量綱一距離飛行距離y＝1處的彎折角對(duì)比，分別為11.20°和12.31°，用JWL狀態(tài)方程的計(jì)算結(jié)果較用多方方程的計(jì)算結(jié)果大9.91%。這與圖8等熵膨脹線所表示的炸藥膨脹做功能力完全一致，即TNT的JWL等熵膨脹線在多方方程等熵線的下方，膨脹做功能力較??；乳化炸藥的JWL等熵線居于多方方程等熵線的上方，做功驅(qū)動(dòng)能力較大。由此可見，通用狀態(tài)方程特征線差分方法能夠正確地反映飛板的爆轟驅(qū)動(dòng)過(guò)程。

圖6 爆轟驅(qū)動(dòng)飛板的拋擲曲線Fig.6Flying curves of flyer plate driven by detonation

圖7 飛板量綱一飛行距離與彎折角的關(guān)系Fig.7Relationship of dimensionless flying distance and bending angle for flyer plate

圖8 JWL等熵線與多方等熵線的對(duì)比Fig.8Comparision of JWL and polytropic isentropic curves

3 結(jié) 論

（1）從馬赫波的物理概念出發(fā)，推導(dǎo)出了平面二維超聲速定常流的特征線方程；并且重新定義了以流體密度為單一自變量的Prantl－Meyer函數(shù)ν（ρ），從而構(gòu)成了求解平面二維超聲速定常流的封閉方程組。由于所獲得的特征線微分方程組不依賴流體狀態(tài)方程的表達(dá)形式，所以可以用于流體各種形式狀態(tài)方程的求解，是一種通用物態(tài)方程的超聲速定常流特征線求解方法。

（2）為了檢驗(yàn)通用物態(tài)方程特征線解法，針對(duì)滑移爆轟驅(qū)動(dòng)飛板運(yùn)動(dòng)問(wèn)題構(gòu)建了爆轟產(chǎn)物流場(chǎng)內(nèi)部和飛板邊界特征線差分法格式，對(duì)爆轟波頭聲速面的處理采用了向超聲速的弱爆轟攝動(dòng)的處理方法。

（3）對(duì)TNT炸藥和乳化炸藥采用JWL狀態(tài)方程與多方方程進(jìn)行對(duì)比計(jì)算的結(jié)果表明，這種通用狀態(tài)方程特征線差分方法可以準(zhǔn)確地反映炸藥狀態(tài)方程對(duì)飛板拋擲的作用，飛板的爆轟驅(qū)動(dòng)過(guò)程與其狀態(tài)方程所表述的做功驅(qū)動(dòng)能力一致。

最后應(yīng)該說(shuō)明的是，特征線差分方法所求解的方程組本質(zhì)上是常微分方程組，所以完全可以參照常微分方程組的數(shù)值求解方法（如歐拉預(yù)估－校正法）獲得高精度特征線差分格式，從而提高定常爆轟驅(qū)動(dòng)問(wèn)題的求解精度。另外，由于特征線差分方法計(jì)算簡(jiǎn)單，并且可以很好地反映滑移爆轟驅(qū)動(dòng)問(wèn)題的實(shí)質(zhì)，如果將通用狀態(tài)方程特征線差分程序作為參數(shù)優(yōu)化程序的內(nèi)嵌，可以構(gòu)成精確簡(jiǎn)單的炸藥狀態(tài)方程參數(shù)擬合程序，簡(jiǎn)化炸藥狀態(tài)方程的參數(shù)擬合計(jì)算過(guò)程。

［1］李曉杰.滑移爆轟下帶覆蓋物的飛板拋擲及對(duì)切割的影響［D］.大連：大連理工大學(xué)，1987.

［2］Besshaposhnikov Y P.Motion of a plate thrown by aglancing detonation wave［J］.Combustion，Explosion，and Shock Waves，1999，35（1）：109－112.

［3］Ivanov A G，Kochkin L I，Ogorodnikov V A，et al.Characteristics of the acceleration of plates by aglancing detonation wave in the presence of an additional or concentrated mass［J］.Combustion，Explosion，and Shock Waves，1990，26（5）：612－614.

［4］邵丙璜，張凱.爆炸焊接原理及其工程應(yīng)用［M］.大連：大連工學(xué)院出版社，1987.

［5］呂洪生，曾新吾.連續(xù)介質(zhì)力學(xué)（中冊(cè)）：流體力學(xué)與爆炸力學(xué)［M］.長(zhǎng)沙：國(guó)防科技大學(xué)出版社，1999.

［6］Courant R，F(xiàn)riedrichs K O.Supersonic flow and shock waves［M］.5th Edition.New York：Springer－Verlag，1976：37－75.

［7］趙春風(fēng).基于通用狀態(tài)方程的爆炸焊接特征線法研究［D］.大連：大連理工大學(xué)，2010.

［8］徐錫申，張萬(wàn)箱.實(shí)用物態(tài)方程理論導(dǎo)引［M］.北京：科學(xué)出版社，1986：402－403.

［9］宋錦泉.乳化炸藥爆轟特性研究［D］.北京：北京科技大學(xué)，2000.