矩形凹槽管道中激波傳播的數(shù)值研究＊

沙 莎，陳志華，韓珺禮，2

（1.南京理工大學(xué)瞬態(tài)物理重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210094；2.北京機(jī)電研究所，北京 100083）

激波在管道內(nèi)傳播過程中經(jīng)常與不同幾何形狀的障礙物相互作用，并發(fā)生反射、繞射、激波誘導(dǎo)渦、激波與渦作用等系列流體物理現(xiàn)象，因而一直受到關(guān)注。另外，研究中發(fā)現(xiàn)，若障礙物形狀合適，則可使激波強(qiáng)度加速衰減，因而可廣泛應(yīng)用于爆炸產(chǎn)生的激波災(zāi)害防治。

D.S.Dosanjh［1］提出利用激波與矩形障礙物作用對(duì)激波進(jìn)行衰減，并進(jìn)行了實(shí)驗(yàn)研究，詳細(xì)分析了激波與矩形障礙物碰撞產(chǎn)生的復(fù)雜波系結(jié)構(gòu)。A.Sasoh等［2］利用全息干涉測(cè)量對(duì)弱激波通過穿孔壁面的過程進(jìn)行了實(shí)驗(yàn)研究，實(shí)驗(yàn)結(jié)果清晰地反映了激波與穿孔壁面作用的詳細(xì)過程。另外，他們還基于Euler方程和TVD（total variation diminishing）格式對(duì)此進(jìn)行了數(shù)值模擬。A.Britan等［3］則對(duì)不同幾何形狀和孔隙率的穿孔壁面中的激波衰減進(jìn)行了實(shí)驗(yàn)及數(shù)值研究。他們采用一維無粘方程對(duì)此進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)管道出口中心處的壓力峰值與障礙物和管道出口的距離成線性反比。N.Gongora－Orozco等［4］對(duì)激波通過5種不同角度彎折的矩形凹槽陣列的衰減進(jìn)行了實(shí)驗(yàn)研究，發(fā)現(xiàn)彎折角越小，出口激波壓力下降越大。D.B.Epstein等［5］基于二維和三維Euler方程以及TVD格式分別對(duì)激波在球形障礙物陣列中的衰減情況進(jìn)行了數(shù)值研究，并分析了入射激波馬赫數(shù)以及障礙物的孔隙率和密度對(duì)衰減的影響。許應(yīng)清等［6］基于Euler方程及GRP（generalized Riemann problem）有限差分法對(duì)激波穿過孔板系統(tǒng)的傳播過程進(jìn)行了數(shù)值模擬。其他對(duì)激波進(jìn)行衰減的方法還包括利用主動(dòng)式水霧［7］、惰性顆粒［8－9］等。激波與多個(gè)障礙物相互作用發(fā)生反射、繞射以及激波誘導(dǎo)渦、激波與渦相互作用導(dǎo)致激波衰減的詳細(xì)過程至今未見報(bào)道。高階精度格式及相關(guān)計(jì)算技術(shù)的發(fā)展使準(zhǔn)確捕捉激波流場(chǎng)得以實(shí)現(xiàn)。為了完全揭示激波與多個(gè)障礙物作用導(dǎo)致其波陣面強(qiáng)度衰減的詳細(xì)過程，本文中基于高精度WENO格式，結(jié)合AMR技術(shù)，對(duì)激波在矩形凹槽管道中的傳播進(jìn)行數(shù)值模擬，并將計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果［4］進(jìn)行比較，分析激波與矩形凹槽相互作用導(dǎo)致激波衰減的機(jī)理，以期為相關(guān)工程應(yīng)用提供依據(jù)。

1 計(jì)算方法與模型

本文中基于二維Euler方程，通過選用五階精度WENO格式，并結(jié)合AMR技術(shù)，對(duì)激波傳播及其與障礙物的相互作用進(jìn)行數(shù)值模擬。

黏性對(duì)激波傳播及其與障礙物相互作用過程的影響較小，因而常被忽略。而WENO格式具有對(duì)激波的高精度分辨率［10］，自提出以來，在流體運(yùn)動(dòng)尤其是包含激波和復(fù)雜流動(dòng)現(xiàn)象的數(shù)值模擬中得到了廣泛應(yīng)用，本文中選用五階精度的WENO格式。另外，結(jié)合AMR技術(shù)，根據(jù)流場(chǎng)壓力變化情況自動(dòng)進(jìn)行網(wǎng)格加密，以盡量少的網(wǎng)格和數(shù)值計(jì)算量來獲取較高的計(jì)算精度。

本文的計(jì)算模型如圖1所示，為了能與實(shí)驗(yàn)結(jié)果［4］進(jìn)行比較，本文中取與實(shí)驗(yàn)中相同的尺寸以及初始邊界條件。管尺寸為125mm×30mm，設(shè)壁面為光滑壁面，且上下壁面均為固壁反射條件，下壁面有15個(gè)矩形凹槽，凹槽高度為7.5mm，寬度為2.5mm，相鄰2個(gè)凹槽之間的距離為5mm。激波管內(nèi)介質(zhì)為空氣。左側(cè)來流馬赫數(shù)Ma＝1.66，波陣面壓力p＝304kPa，右側(cè)為開口條件。

圖1 計(jì)算模型Fig.1The model for numerical computation

2 結(jié)果與討論

不同時(shí)刻，本文中模擬得到的激波與矩形凹槽管道的作用結(jié)果與實(shí)驗(yàn)陰影［4］的比較如圖2所示。由圖2可知，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果［4］吻合較好，且對(duì)作用過程中激波流場(chǎng)細(xì)節(jié)變化揭示得很清晰。

圖2 激波通過管道時(shí)，本文計(jì)算結(jié)果與文獻(xiàn)［4］的實(shí)驗(yàn)結(jié)果的比較Fig.2Computational shadowgraphs for the shock wave propagation through the tube compared with experimental results［4］

2.1 激波與單個(gè)凹槽的相互作用

以激波繞射第1個(gè)凹槽為例，對(duì)激波與單個(gè)矩形凹槽作用過程進(jìn)行討論。圖3為激波繞過第1個(gè)凹槽的紋影時(shí)序圖?？芍肷浼げǚ謩e在凹槽的左側(cè)和右側(cè)發(fā)生繞射和反射，反射激波與繞射誘導(dǎo)的旋渦發(fā)生作用并在凹槽內(nèi)部不斷發(fā)生反射，形成復(fù)雜波系結(jié)構(gòu)。首先，當(dāng)入射激波I到達(dá)第1個(gè)凹槽左上角時(shí)，下端發(fā)生繞射。此時(shí)激波陣面繞射端開始彎曲（D），同時(shí)產(chǎn)生一道與此對(duì)應(yīng)的向上游傳播的扇形膨脹波E。在凹槽左側(cè)凸角處，由于斜壓效應(yīng)，滑移層開始卷起形成一個(gè)順時(shí)針旋渦V，如圖3（a）所示。隨著入射激波I的繼續(xù)傳播，凹槽左側(cè)凸角附近形成一個(gè)典型的激波繞射誘導(dǎo)渦的波系結(jié)構(gòu)，如圖3（b）所示：分離線H、二次激波S和接觸面G［11］。與此同時(shí)，入射激波I與凹槽右側(cè)凸角發(fā)生碰撞，產(chǎn)生一道扇形反射激波R，波后壓力急劇升高，在凹槽右側(cè)凸角附近形成一個(gè)高壓區(qū)。反射激波R一端與入射激波I相連，并隨I向右傳播的同時(shí)向上壁面運(yùn)動(dòng)，另一端與繞射激波D相連，并與D一起在凹槽右側(cè)壁面形成規(guī)則反射，向凹槽底部傳播。隨后，規(guī)則反射演變成馬赫反射，馬赫桿為M，如圖3（b）～（c）所示。反射激波R向左膨脹，與旋渦V相互作用，并發(fā)生變形，進(jìn)而斷裂成2段，R上側(cè)靠近旋渦段因運(yùn)動(dòng)方向與旋渦旋轉(zhuǎn)運(yùn)動(dòng)方向相反而被減速，形成減速激波R1，而下側(cè)因相同而被加速，形成加速激波R2，如圖3（c）所示。減速激波R1向左傳播，前端與駐定二次激波S合并。加速激波R2在凹槽內(nèi)部左右兩側(cè)壁上來回反射，并與旋渦V作用，在凹槽左側(cè)形成一道旋渦誘導(dǎo)波J，如圖3（d）所示。同時(shí)，馬赫桿M不斷變長，隨后在凹槽底端發(fā)生平面反射并向上傳播而與渦相互作用，最終離開凹槽向管內(nèi)傳播，使管內(nèi)波系結(jié)構(gòu)變得更復(fù)雜。而分離線H則在整個(gè)過程中基本不變。

圖3 激波與第1個(gè)凹槽相互作用的過程Fig.3The interaction process of the incident shock wave with the first groove

2.2 激波與凹槽陣列的相互作用

為了方便，本文中以激波繞射前3個(gè)凹槽為例對(duì)激波與矩形凹槽陣列的作用過程進(jìn)行討論。圖4為入射激波繞過前3個(gè)凹槽的時(shí)序陰影圖。首先，如前所述，入射激波I與第1個(gè)凹槽相互作用，下端發(fā)生繞射彎曲（D1），同時(shí)產(chǎn)生一道與此對(duì)應(yīng)的向上游傳播的扇形膨脹波E1，如圖4（a）所示。隨后入射激波I與凹槽右側(cè)凸角發(fā)生碰撞，產(chǎn)生一道扇形反射激波R1，如圖4（b）所示。入射激波I繼續(xù)傳播，與第2個(gè)凹槽發(fā)生作用，該過程和入射激波I與第1個(gè)凹槽的作用過程相同。入射激波I在第2個(gè)凹槽左上角發(fā)生繞射（D2），同時(shí)產(chǎn)生膨脹波E2。隨后入射激波I在第2個(gè)凹槽右上角發(fā)生反射，產(chǎn)生扇形反射激波R2。此時(shí)入射激波與第1個(gè)凹槽作用產(chǎn)生的膨脹波E1以及反射激波R1已向凹槽上方傳播，其右端始終與入射激波I相連，如圖4（c）～（d）所示。隨后入射激波I與第3個(gè)凹槽作用，作用過程和其與前2個(gè)凹槽的作用過程類似，產(chǎn)生第3對(duì)膨脹波－反射波E3－R3，如圖4（e）～（f）所示。產(chǎn)生的每對(duì)膨脹波Ei（下標(biāo)i代表第幾個(gè)凹槽）與反射激波Ri均呈圓弧形向管內(nèi)傳播，它們的一端均與入射激波相連并向管上方傳播，從而影響入射激波波陣面的強(qiáng)度，如圖4（f）所示。

圖4 激波與前3個(gè)凹槽的作用過程Fig.4The interaction processes of incident shock wave with the first three grooves

2.3 入射激波波陣面的強(qiáng)度變化

圖5 t＝46μs時(shí)，入射激波附近的計(jì)算陰影以及對(duì)應(yīng)的波陣面壓力變化曲線Fig.5The computational shadowgraph around incident shock wave at t＝46μs and the corresponding pressure distributions along the incident shock wave front

為了清晰展示激波傳播時(shí)的波陣面強(qiáng)度變化，以入射激波繞第3個(gè)凹槽時(shí)波陣面壓力的變化為例來討論。圖5為t＝46μs時(shí)，入射激波附近陰影及波陣面上不同位置所對(duì)應(yīng)的壓力變化曲線圖。由圖5可知，此時(shí)入射激波已繞過前3個(gè)凹槽，并產(chǎn)生3對(duì)膨脹波－反射激波（E1－R1、E2－R2、E3－R3），它們?cè)谌肷浼げú嚸娴奈恢脛t分別對(duì)應(yīng)于縱坐標(biāo)y1～y6。入射激波波陣面在大于y1的區(qū)域，其壓力稍有下降，但仍可視為初始強(qiáng)度；而在y1～y2、y3～y4、y5～y6區(qū)間內(nèi)，壓力下降較快；相反，在y2～y3、y4～y5以及＜y6的區(qū)間內(nèi)，壓力升高。這主要是因?yàn)閰^(qū)間y1～y2、y3～y4、y5～y6分別對(duì)應(yīng)膨脹波與壓縮波組成的區(qū)域：E1－R1、E2－R2和E3－R3。因膨脹波的波后壓力降低，當(dāng)其與入射激波波陣面相交時(shí)，入射激波波陣面壓力受其影響明顯下降，因而出現(xiàn)以上3個(gè)區(qū)間的壓力下降。相反，反射波則使入射激波波陣面壓力進(jìn)一步升高，因此當(dāng)3個(gè)反射激波R1、R2和R3的前端與入射激波波陣面相交時(shí)，導(dǎo)致入射激波波陣面壓力開始升高，直至與下一個(gè)膨脹波相遇，入射激波波陣面壓力重新開始下降，周而復(fù)始。另外，膨脹波E1導(dǎo)致入射激波波陣面壓力降低的值Δp1顯然大于反射激波R1導(dǎo)致波陣面壓力升高的值Δp2，即Δp1＞Δp2。因而每對(duì)膨脹波－反射波與入射激波相交后，入射激波波陣面壓力峰值都有所下降，從而導(dǎo)致入射激波在與凹槽序列作用過程中，入射激波波陣面壓力呈周期性振蕩下降趨勢(shì)。值得一提的是，管壁附近壓力下降幅度不大。

圖6為管內(nèi)入射激波波陣面壓力和最大壓力隨時(shí)間的變化曲線。圖6（a）為入射激波波陣面3個(gè)不同高度（對(duì)應(yīng)圖1中的1、2、3，且y1＜y2＜y3）上的壓力隨時(shí)間的變化曲線。初始時(shí)，3個(gè)高度的入射激波波陣面壓力相同，均為初始入射壓力，分別在t1、t2、t3時(shí)刻開始下降，且t1＜t2＜t3。可知，距凹槽越近，入射激波波陣面壓力越先受到影響而衰減。3個(gè)不同高度處入射激波波陣面壓力均呈周期性振蕩下降，振蕩幅度平均為Δpy1、Δpy2和Δpy3，顯然Δpy1＞Δpy2＞Δpy3?？梢姡嚯x凹槽越近，入射激波波陣面壓力振蕩幅度越大，即受膨脹波和反射激波的影響越大。然而，各點(diǎn)的入射激波波陣面平均壓力下降相差無幾。圖6（b）為整個(gè)管道中最大壓力隨時(shí)間的變化曲線。初始時(shí)刻，流場(chǎng)最大壓力即為初始入射激波壓力。隨著入射激波與凹槽右上角發(fā)生碰撞以及形成激波反射，波后壓力急劇升高，因而在凹槽右上角附近形成管內(nèi)高壓點(diǎn)。隨后，入射激波在凹槽內(nèi)向下運(yùn)動(dòng)與槽底碰撞，并來回反射，從而在左右底角附近來回產(chǎn)生高壓區(qū)。另外，槽底反射激波與旋渦作用產(chǎn)生2道反射激波，這2道反射激波尾端掃過凹槽左右頂角時(shí)發(fā)生反射，因而產(chǎn)生對(duì)應(yīng)的高壓區(qū)。5個(gè)高壓區(qū)的壓力受流場(chǎng)變化影響，形成壓力擺動(dòng)，組成了圖6（b）所示的管內(nèi)不同時(shí)刻最高壓力變化曲線。

圖6 入射激波波陣面3個(gè)不同高度上的壓力以及流場(chǎng)中最大壓力隨時(shí)間的變化Fig.6Pressure histories of incident shock fronts at three different heights as well as the pressure peak in the flow field

3 結(jié) 論

結(jié)合高精度WENO格式和AMR技術(shù)，對(duì)激波在矩形凹槽管道內(nèi)的傳播過程進(jìn)行了數(shù)值模擬。計(jì)算結(jié)果與已有實(shí)驗(yàn)結(jié)果［4］吻合較好，準(zhǔn)確描述了入射激波與矩形凹槽陣列相互作用的過程，并得到以下主要結(jié)論：（1）入射激波與單個(gè)凹槽作用時(shí)，會(huì)在凹槽左側(cè)發(fā)生繞射而產(chǎn)生膨脹波，入射激波與凹槽右側(cè)碰撞產(chǎn)生反射激波。入射激波與每個(gè)凹槽的相互作用過程相同，均產(chǎn)生一對(duì)膨脹波－反射激波傳入管內(nèi)。（2）產(chǎn)生的膨脹波導(dǎo)致入射激波壓力降低，而反射激波則導(dǎo)致入射激波壓力升高，但其升高值小于膨脹波導(dǎo)致的下降值，因此入射激波波陣面壓力出現(xiàn)振蕩下降趨勢(shì)。（3）入射激波距離凹槽越近的部位受膨脹波和反射激波的影響越早，但最終入射激波波陣面上各點(diǎn)壓力的平均值相差不大。

［1］Dosanjh D S.Interaction of grids with travelling shock waves［R］.NACA－TN－3680，1956.

［2］Sasoh A，Matsuoka K，Nakashio K，et al.Attenuation of weak shock waves along pseudo－perforated walls［J］.Shock Waves，1998，8（3）：149－159.

［3］Britan A，Igra O，Ben－Dor G，et al.Shock wave attenuation by grids and orifice plates［J］.Shock Waves，2006，16（1）：1－15.

［4］Gongora－Orozco N，Zare－Behtash H，Kontis K.Experimental studies on shock wave propagating through junction with grooves［C］∥Proceedings of the 47th AIAA Aerospace Sciences Meeting Including New Horizons Forum and Aerospace Exposition.Florida，USA：American Institute of Aeronautics and Astronautics，2009：5－8.

［5］Epstein D B，Kudryavtsev A N.Shock and blast wave propagation through a porous barrier［C］∥Proceedings of the 28th International Symposium on Shock Waves：Vol 2.Springer Berlin－Heidelberg，2012：537－542.

［6］許應(yīng)清，徐燕侯，吳清松.激波穿過孔板系統(tǒng)的傳播及其誘導(dǎo)的流場(chǎng)［J］.氣動(dòng)實(shí)驗(yàn)測(cè)量與控制，1996，10（2）：69－73.

Xu Ying－qing，Xu Yan－h(huán)ou，Wu Qing－song.The propagation of shock wave through a baffle system and induced flow field［J］.Aerodynamic Experiment and Measurement ＆ Control，1996，10（2）：69－73.

［7］謝波，范寶春，夏自柱，等.大型通道中主動(dòng)式水霧抑爆現(xiàn)象的實(shí)驗(yàn)研究［J］.爆炸與沖擊，2003，23（2）：208－214.

Xie Bo，F(xiàn)an Bao－chun，Xia Zi－zhu，et al.Experimental study of explosion suppression by active water sprays in large－scale duct［J］.Explosion and Shock Waves，2003，23（2）：208－214.

［8］范寶春，李鴻志.惰性顆粒抑爆過程的數(shù)值模擬［J］.爆炸與沖擊，2000，20（3）：208－214.

Fan Bao－chun，Li Hong－zhi.Numerical simulations of explosion suppression by inert particles［J］.Explosion and Shock Waves，2000，20（3）：208－214.

［9］Chen Zhi－h(huán)ua，F(xiàn)an Bao－chun，Jiang Xiao－h(huán)ai.Suppression effects of powder suppressants on the explosions of oxyhydrogen gas［J］.Journal of Loss Prevention in the Process Industries，2006，19（6）：648－655.

［10］Liu X D，Osher S，Chan T.Weighted essentially non－oscillatory schemes［J］.Journal of Computational Physics，1994，115（1）：200－212.

［11］Hillier R.Computation of shock wave diffraction at a ninety degrees convex edge［J］.Shock Waves，1991，1（2）：89－98.