李曉杰,張程嬌,閆鴻浩,王小紅,王宇新
(大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧 大連116024)
水下爆炸的過(guò)程分為炸藥的爆轟、沖擊波在水中的形成及傳播、氣體球膨脹收縮運(yùn)動(dòng)等3個(gè)部分[1]。以往對(duì)遠(yuǎn)場(chǎng)中低壓力區(qū)的爆炸沖擊波和氣泡脈動(dòng)的研究較多,其中遠(yuǎn)場(chǎng)沖擊波壓力與氣泡脈動(dòng)壓力的測(cè)試常被用于炸藥爆炸能量評(píng)定。對(duì)水下爆炸沖擊波傳播問(wèn)題的研究通常采用實(shí)驗(yàn)和數(shù)值模擬的方法。實(shí)驗(yàn)測(cè)試近場(chǎng)沖擊波時(shí),由于沖擊波峰值壓力可以達(dá)到數(shù)萬(wàn)至數(shù)十萬(wàn)個(gè)大氣壓,沖擊波對(duì)壓力傳感器性能要求很高并容易造成嚴(yán)重破壞使得沖擊波壓力難于測(cè)量,因此水下爆炸近場(chǎng)沖擊波的實(shí)驗(yàn)研究比較困難。近年來(lái)對(duì)水下爆炸近場(chǎng)問(wèn)題的研究已逐步采用數(shù)值計(jì)算的方法解決。目前采用的數(shù)值計(jì)算方法很多,但本質(zhì)上仍然是有限差分和有限元2類(lèi),采用有限元或高精差分方法模擬沖擊波間斷問(wèn)題時(shí),在計(jì)算中必須引入人工粘性,對(duì)沖擊波間斷做連續(xù)化處理,所計(jì)算的沖擊波壓力受人工粘性系數(shù)取值的影響[2-4]。盡管人們已經(jīng)通過(guò)將大量計(jì)算結(jié)果與實(shí)驗(yàn)比較,獲得了較合適的人工粘性系數(shù),但人工粘性系數(shù)的使用引入了人為干預(yù)因素。
特征線差分方法[5]是早期流場(chǎng)計(jì)算中常用的一種方法,特征線差分法具有物理概念明確、計(jì)算簡(jiǎn)單、計(jì)算精度高等特點(diǎn)。但由于特征線差分法是根據(jù)流場(chǎng)全域等熵假定(即均熵假定)獲得的,所以只能用于研究均熵流動(dòng)問(wèn)題。盡管早期學(xué)者也使用特征線法對(duì)水下爆炸問(wèn)題進(jìn)行過(guò)研究,但一直忽略了水下變強(qiáng)度沖擊波后的非均熵流特征,而是采用簡(jiǎn)單的均熵假設(shè)來(lái)近似處理沖擊波問(wèn)題,這使得特征線差分法計(jì)算的正確性受到了質(zhì)疑。但是特征線差分方法的優(yōu)點(diǎn)非常明顯:(1)無(wú)需對(duì)沖擊波做連續(xù)化處理,(2)可以直接將流體的偏微分方程轉(zhuǎn)化為計(jì)算精度較高的常微分方程來(lái)求解,(3)特征線本身對(duì)應(yīng)著聲波的傳輸軌跡具有明確的物理意義;所以該方法的研究一直受到重視。
本文中從可壓縮流的基本方程出發(fā),重新推導(dǎo)含有熵變修正項(xiàng)的2族特征線方程,然后將質(zhì)點(diǎn)的跡線方程及其能量方程補(bǔ)充為第3族特征線方程,使得特征線方法適用于非均熵流場(chǎng)的求解。而后利用所獲得的3族特征線方程組,對(duì)質(zhì)量分?jǐn)?shù)比w(TNT)∶w(RDX)=40∶60的TNT/RDX球形炸藥水下爆炸進(jìn)行數(shù)值計(jì)算,并將沖擊波近場(chǎng)壓力的計(jì)算結(jié)果與文獻(xiàn)[6]中的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較,驗(yàn)證這種改進(jìn)的特征線差分法的正確性。
水下爆炸作用時(shí)間極短,通常忽略傳熱效應(yīng),視為絕熱過(guò)程。忽略熱傳導(dǎo)效應(yīng)和粘性效應(yīng)的一維可壓縮流的連續(xù)方程和運(yùn)動(dòng)方程如下
式中:ρ、u、p、r和t分別為流體密度、速度、壓強(qiáng)、空間坐標(biāo)和時(shí)間坐標(biāo)。N=0,1,2分別表示一維平面、柱面和球面流動(dòng)問(wèn)題。以上2式是拉格朗日形式的連續(xù)方程和運(yùn)動(dòng)方程,也就是沿質(zhì)點(diǎn)跡線的流動(dòng)方程。流場(chǎng)內(nèi)質(zhì)點(diǎn)的壓力可以用熵S與密度ρ來(lái)表示p=p(ρ,S),可以寫(xiě)成
式中:下標(biāo)III表示物理量沿質(zhì)點(diǎn)跡線的變化關(guān)系,引入系數(shù)λ,聯(lián)立式(4)與式(2)得到
若令λ=±ρc,那么上式就可以寫(xiě)成
若時(shí)空平面內(nèi)的2族特征線方程為dr/dt=u±c,代入式(6)就可以獲得物理平面上的相容關(guān)系,其中右行特征線方程與相容關(guān)系如下
左行特征線方程與相容關(guān)系如下
式中:(dQ)III=(TdS)III,能量方程中顯含有熵變項(xiàng),對(duì)于一般問(wèn)題必須給出熵S的具體表達(dá)式。但對(duì)于水下爆炸這類(lèi)高速?zèng)_擊問(wèn)題,依然可以假設(shè)流體質(zhì)點(diǎn)之間是絕熱的,即每個(gè)質(zhì)點(diǎn)沿跡線是等熵的,但各質(zhì)點(diǎn)之間的熵變互不相同。因此不要求全流場(chǎng)具有相同的熵變,這就是熵變沿流體質(zhì)點(diǎn)跡線等熵的全場(chǎng)非均熵流,即(dS)III=0,因此能量方程就可以寫(xiě)成隱含熵S的形式,即
同樣地,由于(dS)III=0,式(7)、(8)中的項(xiàng)也為0。式(7)、(8)、(10)就構(gòu)成了可以用于求解非均熵可壓縮流的隱含熵S項(xiàng)3族特征線方程。應(yīng)該注意的是,盡管這3族特征線方程與以往均熵流的方程在形式上幾乎完全相同,但所代表的物理意義卻完全不同,可以用于處理初始熵變非均勻分布的絕熱非均熵流了。例如:對(duì)于球?qū)ΨQ一維水下爆炸這類(lèi)變沖擊波強(qiáng)度的問(wèn)題,由沖擊波強(qiáng)度確定了波陣面的熵變后,沖擊波后方的流場(chǎng)就變成這種沿質(zhì)點(diǎn)跡線等熵的全場(chǎng)非均熵流場(chǎng);另外,對(duì)于單個(gè)沖擊波追趕卸載問(wèn)題也具有同樣的特征,也可以化為這類(lèi)問(wèn)題用3族特征線方程組求解。
從上述方程可見(jiàn),在時(shí)空平面內(nèi),第I、II族特征線方程與標(biāo)準(zhǔn)特征線方程完全相同,均為dr/dt=u±c;在物理平面內(nèi),第I、II族特征線代表了聲波小擾動(dòng)的傳播線,推導(dǎo)的特征線方程和標(biāo)準(zhǔn)特征線方程具有相同的物理意義,不同在于推導(dǎo)的相容關(guān)系中包含了熵變影響項(xiàng)。這也說(shuō)明了在任何復(fù)雜的流場(chǎng)中,其中任意一點(diǎn)的擾動(dòng)均會(huì)沿著聲波小擾動(dòng)線傳播,即第I、II族特征線是物理存在的,而復(fù)雜問(wèn)題在于如何確定特征線上的相容關(guān)系。第3族特征線方程組(10)采用了絕熱等熵假設(shè),如果使用更復(fù)雜的能量方程,該方法可以用于求解帶化學(xué)反應(yīng)或帶耗散的復(fù)雜流動(dòng)。
在應(yīng)用標(biāo)準(zhǔn)特征線法求解均熵流場(chǎng)的問(wèn)題時(shí),只需要前2族特征線就可以求解,其求解方法是用一條右行特征線I與一條左行特征線II相交于某一點(diǎn),該點(diǎn)的各個(gè)參數(shù)由這2條特征線所通過(guò)的上一時(shí)刻的2個(gè)節(jié)點(diǎn)參數(shù)來(lái)確定。如圖1所示,右行特征線I與左行特征線II相交于A點(diǎn),A點(diǎn)的各個(gè)參數(shù)可以用這2條特征線所通過(guò)的前一時(shí)刻的B、C等2點(diǎn)來(lái)求解;如此在給定了計(jì)算域初值條件的情況下,就可以求解出域內(nèi)相關(guān)點(diǎn)的參數(shù)。標(biāo)準(zhǔn)特征線法求解過(guò)程簡(jiǎn)單并且物理意義明確,然而這種方法是以計(jì)算域內(nèi)流場(chǎng)均熵為基礎(chǔ)的,即對(duì)均熵流場(chǎng)求解,在計(jì)算過(guò)程中并未考慮相鄰質(zhì)點(diǎn)間熵值不相等對(duì)計(jì)算結(jié)果的影響。
圖1 標(biāo)準(zhǔn)特征線法在時(shí)空域內(nèi)的差分方法示意圖Fig.1 Standard difference method of characteristics in r-t plane
圖2 改進(jìn)的特征線法的差分格式Fig.2 Difference scheme for improved difference method of characteristics
在絕熱的流場(chǎng)中,各質(zhì)點(diǎn)的熵值沿跡線保持不變,本文的特征線方程中添加了熵修正項(xiàng),并補(bǔ)充加入了質(zhì)點(diǎn)的跡線方程作為第3族特征線方程,因此可以求解非均熵流場(chǎng)。由于引入了3族特征線也就再無(wú)法使用簡(jiǎn)單的差分方法進(jìn)行計(jì)算了,因此必須使用新的差分方法進(jìn)行求解。改進(jìn)特征線法的差分求解方法如下:如圖2所示,假定在t0時(shí)刻點(diǎn)A、B0、C的參數(shù)已知,經(jīng)過(guò)時(shí)間Δt,質(zhì)點(diǎn)由B0沿跡線III運(yùn)動(dòng)至B1處,在時(shí)空平面上可以沿跡線III求出B1點(diǎn)的坐標(biāo);通過(guò)B1點(diǎn)的2條特征線I和II又與t0時(shí)間線相交于點(diǎn)E和點(diǎn)F,這2點(diǎn)的坐標(biāo)及各參數(shù)可以通過(guò)在A、B0與B0、C之間分別插值求得;最后再用點(diǎn)E和點(diǎn)F就可以求出B1點(diǎn)處質(zhì)點(diǎn)的各個(gè)狀態(tài)參數(shù)。若給定初值條件和邊界條件,則如此逐層迭代計(jì)算就可得出整個(gè)域內(nèi)相關(guān)點(diǎn)的參數(shù)。
一維球面水下爆炸問(wèn)題是典型的非均熵流問(wèn)題。爆轟波入射到水中產(chǎn)生球面的水中沖擊波,而水中沖擊波強(qiáng)度隨傳播距離增加而不斷衰減。由于沖擊波引起的熵增與沖擊波強(qiáng)度有關(guān),因此水下爆炸沖擊波后的熵增也隨著傳播距離增加而不斷減小,所以說(shuō)球形炸藥水下爆炸是典型的非均熵流問(wèn)題。
將改進(jìn)的特征線法程序化,對(duì)質(zhì)量分?jǐn)?shù)比w(TNT)∶w(RDX)=40∶60的TNT/RDX炸藥的水下爆炸問(wèn)題進(jìn)行數(shù)值計(jì)算。選用JWL狀態(tài)方程
式中:e為單位質(zhì)量?jī)?nèi)能;V 為相對(duì)體積,V=v/v0,v0=1/ρ0,ρ0為炸藥初始密度,ρ0=1.717g/cm3;方程中的其他參數(shù)分別取值為:A=524.23GPa,B=7.678GPa,R1=4.2,R2=1.1,ω=0.34;炸藥爆速D=7.98km/s。
對(duì)水選用多項(xiàng)式形式的狀態(tài)方程,當(dāng)水受壓和膨脹時(shí),有
式中:μ 是壓縮度,μ=ρ/ρ0-1;e是單位質(zhì)量?jī)?nèi)能;ρ0是初始密度,ρ0=1.0g/cm3;方程中的其他參數(shù):A1=2.2GPa,A2=9.54GPa,A3=14.57GPa,T1=2.2GPa,T2=0,B0=0.28,B1=0.28。
計(jì)算中采用一維球?qū)ΨQ模型,在無(wú)限水域中放入1kg的質(zhì)量分?jǐn)?shù)比w(TNT)∶w(RDX)=40∶60的TNT/RDX炸藥球,藥球半徑為5.18cm。為了避開(kāi)藥球中心處的計(jì)算奇點(diǎn),起爆點(diǎn)設(shè)置在距離球心0.2cm的球面上,起爆后在該處設(shè)置固壁邊界條件;網(wǎng)格大小設(shè)置為0.249cm,炸藥內(nèi)劃分20個(gè)網(wǎng)格。
數(shù)值計(jì)算結(jié)果如圖3所示,圖中繪制出了1≤R/R0≤4范圍內(nèi)的峰值壓力pmax-比例距離R/R0的曲線。為檢驗(yàn)本文中提出的改進(jìn)的特征線法,將計(jì)算結(jié)果與池家春等[6]的實(shí)驗(yàn)結(jié)果進(jìn)行了比較,實(shí)驗(yàn)值公式如下
從圖3中可以看出,沖擊波峰值壓力的計(jì)算值與實(shí)驗(yàn)值符合較好。在比例距離1≤R/R0≤4的近場(chǎng)范圍內(nèi),兩者最大相對(duì)誤差為13.1%,這說(shuō)明本文的特征線法完全可以用于處理水下爆炸近場(chǎng)的非均熵流問(wèn)題。
圖3 近場(chǎng)沖擊波峰值壓力的比較Fig.3 Comparison of the peak pressures in the near-field region
本文中提出了一種改進(jìn)的特征線法,在特征線方程中添加了熵的修正項(xiàng),并將質(zhì)點(diǎn)的跡線作為第3族特征線加入到特征線方程組中,使它可以用于求解變強(qiáng)度沖擊波及波后的非均熵流場(chǎng)。由于時(shí)空平面內(nèi)的第I、II族特征線與標(biāo)準(zhǔn)特征線方程在形式上完全相同,因此非均熵流的前2族特征線方程同樣也代表了聲波小擾動(dòng)的傳播線,與標(biāo)準(zhǔn)特征線方程具有相同的物理意義。
將改進(jìn)的特征線法程序化,對(duì)質(zhì)量分?jǐn)?shù)比w(TNT)∶w(RDX)=40∶60的TNT/RDX炸藥球的水下爆炸問(wèn)題做了計(jì)算。在比例距離1≤R/R0≤4的近場(chǎng)范圍內(nèi),用本文的特征線法計(jì)算得到的沖擊波峰值壓力與實(shí)驗(yàn)值[6]最大誤差為13.1%,從而驗(yàn)證了非均熵特征線法的正確性。另外,這種方法無(wú)需采用人工粘性系數(shù)處理沖擊波間斷,只需給定計(jì)算所需的相關(guān)物理參數(shù)即可求解,因此減少了對(duì)沖擊波的人為干預(yù)因素。
[1]Cole R H.Underwater explosions[M].New Jersy:Princeton University Press,1948:3-13,110.
[2]劉科種,徐更光,辛春亮,等.AUTODYN水下爆炸數(shù)值模擬研究[J].爆破,2009,26(3):18-21.LIU Ke-zhong,XU Geng-guang,XIN Chun-liang,et al.Research on numerical simulation in underwater explosioln by AUTODYN[J].Blasting,2009,26(3):18-21.
[3]徐豫新,王樹(shù)山,李園.水下爆炸數(shù)值仿真研究[J].彈箭與制導(dǎo)學(xué)報(bào),2009,29(6):95-97,102.XU Yu-xin,WANG Shu-shan,LI Yuan.Study on numerical simulation of the underwater explosive[J].Journal of Projectiles,Rockets,Missiles and Guidance,2009,29(6):95-97,102.
[4]方斌,朱錫,張振華,等.水下爆炸沖擊波數(shù)值模擬中的參數(shù)影響[J].哈爾濱工程大學(xué)學(xué)報(bào),2005,26(4):419-424.FANG Bin,ZHU Xi,ZHANG Zhen-h(huán)ua,et al.Effect of parameters in numerical simulation of underwater shock wave[J].Journal of Harbin Engineering University,2005,26(4):419-424.
[5]Chou P C,Huang S L,Karpp R R.Numerical calculation of blast wave by the method of characteristics[J].AIAA Journal,1967,5(4):618-623.
[6]池家春,馬冰.TNT/RDX(40/60)炸藥球水中爆炸波研究[J].高壓物理學(xué)報(bào),1999,13(3):199-204.CHI Jia-chun,MA Bing.Underwater explosion wave by a spherical charge of composition B-3[J].Chinese Journal of High Pressure Physics,1999,13(3):199-244.