矩形凹槽管道中激波传播的数值研究＊

沙 莎，陈志华，韩珺礼，2

（1.南京理工大学瞬态物理重点实验室，江苏 南京 210094；2.北京机电研究所，北京 100083）

激波在管道内传播过程中经常与不同几何形状的障碍物相互作用，并发生反射、绕射、激波诱导涡、激波与涡作用等系列流体物理现象，因而一直受到关注。另外，研究中发现，若障碍物形状合适，则可使激波强度加速衰减，因而可广泛应用于爆炸产生的激波灾害防治。

D.S.Dosanjh［1］提出利用激波与矩形障碍物作用对激波进行衰减，并进行了实验研究，详细分析了激波与矩形障碍物碰撞产生的复杂波系结构。A.Sasoh等［2］利用全息干涉测量对弱激波通过穿孔壁面的过程进行了实验研究，实验结果清晰地反映了激波与穿孔壁面作用的详细过程。另外，他们还基于Euler方程和TVD（total variation diminishing）格式对此进行了数值模拟。A.Britan等［3］则对不同几何形状和孔隙率的穿孔壁面中的激波衰减进行了实验及数值研究。他们采用一维无粘方程对此进行了数值模拟，发现管道出口中心处的压力峰值与障碍物和管道出口的距离成线性反比。N.Gongora－Orozco等［4］对激波通过5种不同角度弯折的矩形凹槽阵列的衰减进行了实验研究，发现弯折角越小，出口激波压力下降越大。D.B.Epstein等［5］基于二维和三维Euler方程以及TVD格式分别对激波在球形障碍物阵列中的衰减情况进行了数值研究，并分析了入射激波马赫数以及障碍物的孔隙率和密度对衰减的影响。许应清等［6］基于Euler方程及GRP（generalized Riemann problem）有限差分法对激波穿过孔板系统的传播过程进行了数值模拟。其他对激波进行衰减的方法还包括利用主动式水雾［7］、惰性颗粒［8－9］等。激波与多个障碍物相互作用发生反射、绕射以及激波诱导涡、激波与涡相互作用导致激波衰减的详细过程至今未见报道。高阶精度格式及相关计算技术的发展使准确捕捉激波流场得以实现。为了完全揭示激波与多个障碍物作用导致其波阵面强度衰减的详细过程，本文中基于高精度WENO格式，结合AMR技术，对激波在矩形凹槽管道中的传播进行数值模拟，并将计算结果与实验结果［4］进行比较，分析激波与矩形凹槽相互作用导致激波衰减的机理，以期为相关工程应用提供依据。

1 计算方法与模型

本文中基于二维Euler方程，通过选用五阶精度WENO格式，并结合AMR技术，对激波传播及其与障碍物的相互作用进行数值模拟。

黏性对激波传播及其与障碍物相互作用过程的影响较小，因而常被忽略。而WENO格式具有对激波的高精度分辨率［10］，自提出以来，在流体运动尤其是包含激波和复杂流动现象的数值模拟中得到了广泛应用，本文中选用五阶精度的WENO格式。另外，结合AMR技术，根据流场压力变化情况自动进行网格加密，以尽量少的网格和数值计算量来获取较高的计算精度。

本文的计算模型如图1所示，为了能与实验结果［4］进行比较，本文中取与实验中相同的尺寸以及初始边界条件。管尺寸为125mm×30mm，设壁面为光滑壁面，且上下壁面均为固壁反射条件，下壁面有15个矩形凹槽，凹槽高度为7.5mm，宽度为2.5mm，相邻2个凹槽之间的距离为5mm。激波管内介质为空气。左侧来流马赫数Ma＝1.66，波阵面压力p＝304kPa，右侧为开口条件。

图1 计算模型Fig.1The model for numerical computation

2 结果与讨论

不同时刻，本文中模拟得到的激波与矩形凹槽管道的作用结果与实验阴影［4］的比较如图2所示。由图2可知，计算结果与实验结果［4］吻合较好，且对作用过程中激波流场细节变化揭示得很清晰。

图2 激波通过管道时，本文计算结果与文献［4］的实验结果的比较Fig.2Computational shadowgraphs for the shock wave propagation through the tube compared with experimental results［4］

2.1 激波与单个凹槽的相互作用

以激波绕射第1个凹槽为例，对激波与单个矩形凹槽作用过程进行讨论。图3为激波绕过第1个凹槽的纹影时序图。可知，入射激波分别在凹槽的左侧和右侧发生绕射和反射，反射激波与绕射诱导的旋涡发生作用并在凹槽内部不断发生反射，形成复杂波系结构。首先，当入射激波I到达第1个凹槽左上角时，下端发生绕射。此时激波阵面绕射端开始弯曲（D），同时产生一道与此对应的向上游传播的扇形膨胀波E。在凹槽左侧凸角处，由于斜压效应，滑移层开始卷起形成一个顺时针旋涡V，如图3（a）所示。随着入射激波I的继续传播，凹槽左侧凸角附近形成一个典型的激波绕射诱导涡的波系结构，如图3（b）所示：分离线H、二次激波S和接触面G［11］。与此同时，入射激波I与凹槽右侧凸角发生碰撞，产生一道扇形反射激波R，波后压力急剧升高，在凹槽右侧凸角附近形成一个高压区。反射激波R一端与入射激波I相连，并随I向右传播的同时向上壁面运动，另一端与绕射激波D相连，并与D一起在凹槽右侧壁面形成规则反射，向凹槽底部传播。随后，规则反射演变成马赫反射，马赫杆为M，如图3（b）～（c）所示。反射激波R向左膨胀，与旋涡V相互作用，并发生变形，进而断裂成2段，R上侧靠近旋涡段因运动方向与旋涡旋转运动方向相反而被减速，形成减速激波R1，而下侧因相同而被加速，形成加速激波R2，如图3（c）所示。减速激波R1向左传播，前端与驻定二次激波S合并。加速激波R2在凹槽内部左右两侧壁上来回反射，并与旋涡V作用，在凹槽左侧形成一道旋涡诱导波J，如图3（d）所示。同时，马赫杆M不断变长，随后在凹槽底端发生平面反射并向上传播而与涡相互作用，最终离开凹槽向管内传播，使管内波系结构变得更复杂。而分离线H则在整个过程中基本不变。

图3 激波与第1个凹槽相互作用的过程Fig.3The interaction process of the incident shock wave with the first groove

2.2 激波与凹槽阵列的相互作用

为了方便，本文中以激波绕射前3个凹槽为例对激波与矩形凹槽阵列的作用过程进行讨论。图4为入射激波绕过前3个凹槽的时序阴影图。首先，如前所述，入射激波I与第1个凹槽相互作用，下端发生绕射弯曲（D1），同时产生一道与此对应的向上游传播的扇形膨胀波E1，如图4（a）所示。随后入射激波I与凹槽右侧凸角发生碰撞，产生一道扇形反射激波R1，如图4（b）所示。入射激波I继续传播，与第2个凹槽发生作用，该过程和入射激波I与第1个凹槽的作用过程相同。入射激波I在第2个凹槽左上角发生绕射（D2），同时产生膨胀波E2。随后入射激波I在第2个凹槽右上角发生反射，产生扇形反射激波R2。此时入射激波与第1个凹槽作用产生的膨胀波E1以及反射激波R1已向凹槽上方传播，其右端始终与入射激波I相连，如图4（c）～（d）所示。随后入射激波I与第3个凹槽作用，作用过程和其与前2个凹槽的作用过程类似，产生第3对膨胀波－反射波E3－R3，如图4（e）～（f）所示。产生的每对膨胀波Ei（下标i代表第几个凹槽）与反射激波Ri均呈圆弧形向管内传播，它们的一端均与入射激波相连并向管上方传播，从而影响入射激波波阵面的强度，如图4（f）所示。

图4 激波与前3个凹槽的作用过程Fig.4The interaction processes of incident shock wave with the first three grooves

2.3 入射激波波阵面的强度变化

图5 t＝46μs时，入射激波附近的计算阴影以及对应的波阵面压力变化曲线Fig.5The computational shadowgraph around incident shock wave at t＝46μs and the corresponding pressure distributions along the incident shock wave front

为了清晰展示激波传播时的波阵面强度变化，以入射激波绕第3个凹槽时波阵面压力的变化为例来讨论。图5为t＝46μs时，入射激波附近阴影及波阵面上不同位置所对应的压力变化曲线图。由图5可知，此时入射激波已绕过前3个凹槽，并产生3对膨胀波－反射激波（E1－R1、E2－R2、E3－R3），它们在入射激波波阵面的位置则分别对应于纵坐标y1～y6。入射激波波阵面在大于y1的区域，其压力稍有下降，但仍可视为初始强度；而在y1～y2、y3～y4、y5～y6区间内，压力下降较快；相反，在y2～y3、y4～y5以及＜y6的区间内，压力升高。这主要是因为区间y1～y2、y3～y4、y5～y6分别对应膨胀波与压缩波组成的区域：E1－R1、E2－R2和E3－R3。因膨胀波的波后压力降低，当其与入射激波波阵面相交时，入射激波波阵面压力受其影响明显下降，因而出现以上3个区间的压力下降。相反，反射波则使入射激波波阵面压力进一步升高，因此当3个反射激波R1、R2和R3的前端与入射激波波阵面相交时，导致入射激波波阵面压力开始升高，直至与下一个膨胀波相遇，入射激波波阵面压力重新开始下降，周而复始。另外，膨胀波E1导致入射激波波阵面压力降低的值Δp1显然大于反射激波R1导致波阵面压力升高的值Δp2，即Δp1＞Δp2。因而每对膨胀波－反射波与入射激波相交后，入射激波波阵面压力峰值都有所下降，从而导致入射激波在与凹槽序列作用过程中，入射激波波阵面压力呈周期性振荡下降趋势。值得一提的是，管壁附近压力下降幅度不大。

图6为管内入射激波波阵面压力和最大压力随时间的变化曲线。图6（a）为入射激波波阵面3个不同高度（对应图1中的1、2、3，且y1＜y2＜y3）上的压力随时间的变化曲线。初始时，3个高度的入射激波波阵面压力相同，均为初始入射压力，分别在t1、t2、t3时刻开始下降，且t1＜t2＜t3。可知，距凹槽越近，入射激波波阵面压力越先受到影响而衰减。3个不同高度处入射激波波阵面压力均呈周期性振荡下降，振荡幅度平均为Δpy1、Δpy2和Δpy3，显然Δpy1＞Δpy2＞Δpy3。可见，距离凹槽越近，入射激波波阵面压力振荡幅度越大，即受膨胀波和反射激波的影响越大。然而，各点的入射激波波阵面平均压力下降相差无几。图6（b）为整个管道中最大压力随时间的变化曲线。初始时刻，流场最大压力即为初始入射激波压力。随着入射激波与凹槽右上角发生碰撞以及形成激波反射，波后压力急剧升高，因而在凹槽右上角附近形成管内高压点。随后，入射激波在凹槽内向下运动与槽底碰撞，并来回反射，从而在左右底角附近来回产生高压区。另外，槽底反射激波与旋涡作用产生2道反射激波，这2道反射激波尾端扫过凹槽左右顶角时发生反射，因而产生对应的高压区。5个高压区的压力受流场变化影响，形成压力摆动，组成了图6（b）所示的管内不同时刻最高压力变化曲线。

图6 入射激波波阵面3个不同高度上的压力以及流场中最大压力随时间的变化Fig.6Pressure histories of incident shock fronts at three different heights as well as the pressure peak in the flow field

3 结 论

结合高精度WENO格式和AMR技术，对激波在矩形凹槽管道内的传播过程进行了数值模拟。计算结果与已有实验结果［4］吻合较好，准确描述了入射激波与矩形凹槽阵列相互作用的过程，并得到以下主要结论：（1）入射激波与单个凹槽作用时，会在凹槽左侧发生绕射而产生膨胀波，入射激波与凹槽右侧碰撞产生反射激波。入射激波与每个凹槽的相互作用过程相同，均产生一对膨胀波－反射激波传入管内。（2）产生的膨胀波导致入射激波压力降低，而反射激波则导致入射激波压力升高，但其升高值小于膨胀波导致的下降值，因此入射激波波阵面压力出现振荡下降趋势。（3）入射激波距离凹槽越近的部位受膨胀波和反射激波的影响越早，但最终入射激波波阵面上各点压力的平均值相差不大。

［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］许应清，徐燕侯，吴清松.激波穿过孔板系统的传播及其诱导的流场［J］.气动实验测量与控制，1996，10（2）：69－73.

Xu Ying－qing，Xu Yan－hou，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］谢波，范宝春，夏自柱，等.大型通道中主动式水雾抑爆现象的实验研究［J］.爆炸与冲击，2003，23（2）：208－214.

Xie Bo，Fan 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］范宝春，李鸿志.惰性颗粒抑爆过程的数值模拟［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－hua，Fan Bao－chun，Jiang Xiao－hai.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.