0316 收稿日期:20210316
修回日期:20210702
移动阅读
CHEN Changjiang,LIU Yong,WEI Jianping,et al. Flow field structure and pulsation frequency of air jet under different pressure ratios[J]. Journal of China Coal Society,2021,46(12):3883-3890.
高压磨料气体射流辅助钻井可有效避免水射流辅助钻进松软煤层、泥岩等水敏地层时，水分引入导致的矿物膨胀或井筒壁坍塌等问题[1-2]。另外，高压磨料气体射流以高压空气作为动力源，还可以减少水资源的污染和浪费，便于磨料粒子的二次回收利用[3]。然而，目前常用的磨料气体射流具有粒子加速效率低、冲击力小的缺点[4-5]，不能满足辅助钻井的要求。鉴于磨料气体射流的优、劣势，LIU等[6-7]提出了以Laval喷嘴为基础，利用超音速气流加速粒子的方式[6-7]。研究得出，超音速气流加速磨料粒子可有效提高粒子加速和射流破碎岩石效果，且随着气流压力或速度增加，磨料在时间和空间上出现了不均匀分布特征，该特征可有效降低磨料射流冲击靶体时的shield效应，提高冲蚀效率。根据气固两相流理论，气体流场的时空特征决定了磨料分布规律[8]。因此，为有效调控磨料在时间和空间上的分布规律，进一步提高磨料射流破岩效果，有必要对超音速气体射流的流场空间结构及脉动规律进行研究。国内外学者对气体射流流场结构及脉动频率开展了相关研究。其中，何枫、FRENDI等[9-11]研究表明，超音速气体射流冲击靶板过程中会存在离散的尖锐“啸声”，即气体射流流场出现了固定频率的振荡。流场振荡是由于流场内的激波导致了流场压力、速度、密度等参数出现了脉动。又由空气动力学可知，超音速气体流场内部的激波结构的形成和空间位置决定于射流压力和环境压力的比值，即膨胀比β[12]。ABEDI，ZHAO等[13-15]通过测量不同位置处气流的物性参数(压力、密度等)，研究了超音速燃料混和过程中流场的空间特征，得到相同膨胀比下流场具有相似的空间结构，但改变膨胀比、流场结构可能会存在突变。通过上述介绍，可知膨胀比是决定气体流场空间结构及脉动规律的重要因素，但关于膨胀比对超音速气体射流流场时空演化规律的研究较少。为此，笔者通过测定不同膨胀比β条件下，流体的物性参数随射流流场空间位置及时间变化规律，研究了膨胀比对超音速气体射流流场时空演化规律的影响。
对于流体参数测定可选择数值计算及试验测量。数值计算是通过分析流场流动状态，计算理论状态下变量对流场的影响[16-17]。由于高压磨料气体射流的流场中，气流物理性质复杂，选用的数学模型仅能等效于单方面的物性参数。此时，理论计算结果与流场的实际发展过程偏差较大。纹影方法是利用光学设备非接触式测量流场的内部结构，通过选择高分辨率、高灵敏度的纹影原件可获得毫秒量级上气体流场的灰度照片[18-19]。经标定、校准、编码计算，纹影照片可转换为气体射流过程中的密度场。
笔者采用纹影试验获得了不同膨胀比条件下超音速气体射流流场的灰度照片，并利用照片的灰度值计算流场的密度分布，又通过对不同时刻纹影照片的灰度值进行Fourier变换，计算流场的脉动频率。
纹影仪是根据沿程折射率的改变，光路会发生偏折的特性，通过在光路的传播方向上设置“刀口”对偏折后的光线遮挡，进而利用光强的变化，在高速摄像机镜头或光屏上形成清晰的灰度影像。目前，纹影试验主要被应用于风洞气流及火焰传播过程中边界层、激波结构、密度梯度等流场信息的定性研究[20]。通过采用高分辨率、高灵敏度的纹影原件及高拍摄帧率的高速像机，可得到流场实时、清晰的纹影照片。又根据空气动力学及相关研究，喷嘴膨胀比β决定了气流的压缩或膨胀程度，并影响射流流场结构及脉动规律[9-15]。因此，笔者采用纹影法采集了不同膨胀比下超音速气体射流流场的纹影图像，并经过再处理，定量研究了流场的密度分布特征及脉动规律。
超音速气体射流流场纹影试验研究的设备包括高压气体射流试验系统和纹影仪系统，如图1所示。
图1 气体射流及纹影仪系统
Fig.1 System diagram of airjet and schlieren
高压气体射流试验系统主要有高压空气压缩机、高压气瓶、数字压力表、调压阀、射流喷嘴构成。空气压缩机最高压力为40 MPa，最大吸气量2 m3/min；高压气瓶公称体积为100 L，设计压力为40 MPa，满足试验中气压和气量的要求；数字压力表量程为0～40 MPa，最小分度值为0.01 MPa；调压阀的进口压力调节范围为0～40 MPa，出口压力调节范围为0～30 MPa，出口压力可调精确度为0.1 MPa。射流喷嘴包括有出口膨胀比β=1.00，1.12，2.00，3.00的4种Laval喷嘴，喷嘴结构如图2所示，参数见表1。高压气体射流试验系统可保证射流出口压力的稳定时间大于1 min，满足射流试验的测试要求。
图2 喷嘴结构
Fig.2 Nozzle structure
表1 喷嘴参数
Table 1 Nozzle parameters
注：d1为入口直径；d2为收敛段直径；d3为出口直径；D为喷嘴直径；L1为收敛段长度；L2为扩张段长度。
βd1/mmd2/mmd3/mmD/mmL1/mmL2/mm1.004.82.06.614.01.821.31.124.82.06.214.01.819.52.004.82.05.214.01.815.23.004.82.04.614.01.812.8
纹影仪的元件包括有卤钨灯光源、纹影镜M1，M2、刀口装置、高速摄像机，如图1所示。为减少或消除流场成像时的相差，调整纹影仪元件的摆放位置，采用“Z”形光路测量流场信息，如图1所示。其中，纹影镜M1，M2是镜面直径D=300 mm的凹面镜，且为满足拍摄灵敏度，选择的凹面镜焦距f=2.6 m；试验中由于喷嘴出口处气流的马赫数为3～4，选择高速摄像机的拍摄帧率为3 000 Hz，像素为1 240×1 240；卤钨灯的电压U=24 V，功率P=300 W，曝光时间T=500 μs；纹影设备的“刀口”装置对焦完成后，可实现的拍摄范围直径为300 mm。
试验前，开启空气压缩机为高压气瓶充气，满足试验压力要求后关闭空气压缩机，待压力表示数稳定一段时间后，在射流系统出口处安装射流喷嘴；组装纹影仪，依照“Z”形光路摆放元件，打开光源及高速摄像机，调整高速摄像机像素为1 240×1 240，移动刀口位置进行对焦；为了计算密度场，试验前还需要调节，记录刀口的上、下进刀位置，并采集“刀口”处于不同位置时的灰度图像。试验时，保证射流入口压力为15 MPa，通过替换射流喷嘴，依次对β=1.00，1.12，2.00，3.00的流场进行拍摄取样，每次拍摄时间为4 s，拍摄流场长度为300 mm，每次取样3次。
在均匀介质中光线沿直线传播，而当介质中的密度不均匀时，会导致同一透明介质的局部折射率发生变化，使得光线传播发生扰动。扰动区域的图像会与周围的图像形成差异，即可见纹影[21]。但为实现流场空间结构特征的定量化研究，得到的纹影图像还需再处理，转化为流场的密度值。
纹影法是一种将光通过透明非均匀介质引起的相位分布转换为光强分布的方法。
根据Huygens原理，光线折射率n和偏转角θ的方程 [21] 为
(1)
式中，L为流场宽度；y为光线偏折的距离。
光线在刀口处的偏移量a与光线的偏转角θ之间的关系可表示为
(2)
式中，f2为纹影镜的焦距。
根据Gladstone-Dale(G-D)原理，折射率n与射流流场密度ρ的方程[22]可表示为
n=1+Kρ
(3)
式中，K为可见光的Gladstone-Dale常数[23]，K=2.22×10-4+2.97×10-19/λ2，λ为可见光的波长，波长范围390～780 nm。
结合式(1)～(3)可以得到流场的密度梯度公式为
(4)
由于Kρ的数量级为10-4，方程(4)可表示为
(5)
由于沿程折射率改变，光路发生偏折，经过“刀口”对偏折后的光线遮挡，纹影照片上出现明暗信息。当入射光强相同时，光线在“刀口”上的偏移量a，决定了照片的灰度值H。假设灰度值H与光线在刀口处的偏移量a满足:
a=P(H)
(6)
式中，P为函数表达式。
结合式(4)，(5)，可得流场密度与灰度值的关系为
(7)
式中，ξ0, ξ1为刀口在高速摄像机上的开始和结束位置；ρ0为环境密度。
采用MATLAB对式(7)的密度场计算模型进行编码，可实现照片灰度值转化为流场密度值。
为了提取流场图像的密度信息，需要对图像的灰度值和光照强度进行标定，获得式(6)的函数关系式。根据纹影照片的灰度值H决定于光线在“刀口”上的偏移量a，可通过研究纹影仪拍摄区域未加入流场时，“刀口”位置与照片平均灰度值之间的关系(图3)，得出式(6)具体表达式。即试验前，调节、记录“刀口”的进刀位置，并利用高速摄像机采集“刀口”处于不同位置时的灰度图像。通过MATLAB提取图像的平均灰度值，得出灰度值与“刀口”位置的关系。
图3 灰度值与“刀口”位置关系
Fig.3 Relationship between gray scale value and“knife edge” position
试验得到的“刀口”位置与图像灰度值之间的关系，且当灰度值在20～160时，灰度值H与“刀口”位置a满足一次函数。因此，式(6)可以表示为
H=117a-56
(8)
其拟合系数R2=0.99。
通过标定灰度值和“刀口”位置，获得了灰度值和光强的关系，即可以采用式(7)模型，由灰度值计算流场的密度值。
通过改变射流喷嘴结构，得到了喷嘴出口处膨胀比β=1.00，1.12，2.00，3.00的超音速气体射流流场的纹影图像，如图4(a)所示。又基于2.1节的流场密度计算模型，图像灰度值被转换成了流场密度，并通过灰度-颜色映射处理，得到了伪彩色的流场密度云图，如图4(b)所示。
如图4所示，不同膨胀比条件下射流形成的流场与环境之间均有明显自由界面，且穿过边界层，气流密度会迅速衰减为环境密度值。膨胀比β=1.00时，气体射流处于完全膨胀状态，流场中的密度变化平缓，压力变化波动较小，能量衰减缓慢，没有明显的波结构；β=1.12时，气体射流出现轻微膨胀，形成类圆筒状的波节结构，但由于膨胀-压缩效果不明显，维持的射流等速核较长；β=2.00，3.00时，由于超音速气流表现出明显的压缩-膨胀性，导致射流流场核心段沿射流轴向交替出现了高密度区和低密度区，核心段长度快速衰减。
膨胀比β>1时，如图4中β=1.12, 2.00，3.00所示。喷嘴处气流处于欠膨胀状态，静压力高于环境压力，压力差产生的扰动会沿喷嘴出口壁向流场下游形成圆锥状膨胀波系。气体射流穿过波系的波阵面，向自由界面偏折，流场径向中心出现低密度区，随着气流向下游运动，膨胀波系在射流自由界面反射形成压缩波系，气流穿过压缩波系的波阵面界面发生压缩叠加，动压力降低，静压力升高，温度升高，出现高密度区。随着波阵面沿射流轴向发展，压缩波系会再次反射为膨胀波系，流场内出现低密度区；同样，膨胀波系会再次反射为压缩波系，流场内静压力升高，出现高密度区。但由于膨胀-压缩波系循环，沿程射流压力势能和动能不断消散，气体可压缩性减弱，最终，类圆筒状的波节结构沿射流轴向收敛、消失。
图4 气体流场的灰度及密度云图
Fig.4 Gray scale and density contour image of air flow field
随着膨胀比β增加，气体射流经过膨胀波，气流压力会迅速衰减，且低于周围环境压力，即图4(b)中β=2.00,3.00的密度云图中的过膨胀区域。由于压力差的增加，膨胀波系在射流自由界面会反射形成激波。气体射流穿过激波后，气流的密度会急剧升高，如图4(b)中β=2.00所示，气流穿过过膨胀区域后，流场内出现了高密度区，密度最大值达到了5.6 kg/m3。同样，激波在射流自由界面会再次反射为膨胀波系，流场内会再次出现低密度区。密度的剧增或剧降导致了流场在空间上振荡性增加，且由于密度的快速变化，射流沿程能量损失严重，核心段长度会快速衰减。当膨胀比β继续增加，流场内膨胀-压缩效果继续增强，局部区域的激波会持续发生叠加，直至气流出现密度骤增、速度骤降等强间断扰动现象，形成马赫盘，如图4(b)中β=3.00所示。超音速气流穿过马赫盘等强间断扰动，能量损失严重，气流速度迅速衰减为亚音速，气流变成不可压缩流体，密度参数接近于周围环境密度值，流场空间上的振荡消失。
由于膨胀比β增加，流场出现波节结构，气流密度变化加剧，但β=1.12时，由于膨胀-压缩效果不明显，维持的射流等速核较长。随着膨胀比β增加，气流密度沿射流轴向交替出现剧增、剧降，流场空间结构的振荡性增加，核心段长度快速衰减。继续增加膨胀比β，气流密度会迅速衰减为环境密度值，振荡性减弱。水射流加速磨料的研究表明，等速核长度越长，越利于磨料粒子的加速[24]。因此，β=1.12时，超音速气体射流加速磨料的效果要优于β=1.00，2.00，3.00。
气体流场内复杂的膨胀-压缩波系，导致气体流场空间结构呈现出振荡分布特征，密度沿射流轴向振荡衰减。流场在空间振荡的同时，还会存在时间上脉动。为研究流场脉动频率，对纹影仪采集到的灰度图像进行了坐标处理，并利用MATLAB依次提取了射流轴线上距射流喷嘴出口50，100，150，200 mm位置不同时刻的灰度值。通过灰度值与时间的对应关系得出了气流密度参数的时域信息，再采用Fourier变换对时域数据进行了频谱分析。
对时域数据进行频域转换后，采用频谱分析可以研究不同膨胀比条件下自由流场的脉动频率的变化规律；其中，Fourier变换是进行时域到频域转化研究的主要方法[25]。
设f(t)为t时刻的函数值，对f(t)进行Fourier变换时，可得
F(w)=
f(t)e-iwtdt
(9)
式中，w为周期函数的频率；F(w)为频率w对应的振幅；i为虚部。
由于采集到的样本是离散的数据，需对连续傅里叶变换做离散化处理，式(9)可变为
(10)
式中:Xm为样本值；X(m·Δf)为样本函数；Δf为频率分辨率；φ为时间数；m为第m个采样点；M为采样点数量；Δt为采样时间间隔。
令频率f=0，即m=0时，有
(11)
由式(11)可以看出，频谱分析中频率分量为0时所对应的幅值表示流场中密度的平均值。密度的平均值远大于脉动值，直接对采集到的时域信号进行Fourier变换后得到的频谱图中，频率为0的位置幅值变化大，掩盖了射流的脉动信息。为研究不同膨胀比条件下流场脉动频率的变化规律，频域变换过程中选取灰度系数Cp代替所采集的灰度值。
(12)
式中，
为采集数据的平均灰度值。
通过提取不同时刻纹影图像上的灰度值，可得到流场时域变化信息。对流场时域信息采用3.1节的频谱计算方法进行频域转换，可得到流场的频域变化信息。
经过对得到的50，100，150，200 mm处的频域信息进行分析，得出气体流场的脉动频率决定于喷嘴膨胀比β，而与空间位置无关。因此，笔者以射流喷嘴出口150 mm处的流场时域及频域信息为例，分析了不同膨胀比条件下流场脉动频率的变化规律。图5为射流喷嘴出口150 mm处的流场时域及频域信息。
如图5所示，气体流场的脉动效果，决定于射流喷嘴膨胀比β。β=1.00时，气体流场的脉动频率集中于100 Hz；β=1.12时，自由流场中存在多个脉动主频率，频率值分别集中于100，550，1 100，1 400 Hz；β=2.00时，气体流场的脉动频率集中于50 Hz；β=3.00时，气体流场的脉动频率又变为100 Hz。
由2.3节可知，β=1.00时，流场内部没有波节结构，密度变化平缓，相对应的脉动频率会集中，即如图5(a)所示，流场脉动的主频率集中于100 Hz。随着膨胀比β增大，喷嘴出口处气流进入欠膨胀状态，气体流场内部开始出现膨胀-压缩波系，流场的脉动特性增强，即流场的脉动主频率增多，频率增大，如图5(b)所示，β=1.12时，自由流场中存在多个脉动主频率，频率值分别集中于100，550，1 100，1 400 Hz。随着膨胀比β继续增大，流场内部出现激波，气体射流的沿程压力势能和动能迅速衰减，流场的脉动强度受到抑制，脉动频率开始减小，直至β=2.00时，流场的脉动主频率减小到50 Hz，如图5(c)所示。膨胀比β继续增大，会导致流场内部的垂直激波发展为马赫盘，随着马赫盘的产生和移动，流场空间上的振荡性消失，流场的主频率也回归到了膨胀比β=1.00时的频率值，如图5(d)所示。
图5 气体流场的时域及频域信息
Fig.5 Time-domain and frequency-domain information of air flow field
气体射流随膨胀比β增加，流场脉动强度依次表现为增强、减弱、趋于稳定的变化趋势，且β=1.12时，气体流场会存在多个脉动主频率。研究表明，射流加速磨料粒子时，适当的增加脉动频率有利于磨料粒子向核心射流段掺混，有更好的加速磨料效果[26]。
(1)膨胀比不同，流场的空间结构不同。当膨胀比β=1.00时，流场内部不存在波节结构；随膨胀比增加，流场开始出现波节结构。当β=1.12时，由于膨胀-压缩效果不明显，维持的射流等速核较长；而随膨胀比继续增加，流场空间结构的振荡性增加，核心段长度会快速衰减。
(2)膨胀比不同，流场内气流密度的变化规律不同。当膨胀比β=1.00时，流场内气流的密度变化平缓；膨胀比增加，气流密度变化加剧；而随膨胀比继续增加，气流密度会迅速衰减为环境密度值。
(3)超音速气体射流流场的脉动强度决定于膨胀比β，随着随膨胀比β增加，流场脉动强度依次表现为增强、减弱、趋于稳定的变化趋势，且当β=1.12时，气体流场会存在多个脉动主频率分别为100，550，1 100，1 400 Hz。
[1] 刘勇,陈长江,魏建平,等. 磨料水射流与磨料气体射流破岩压力对比分析[J]. 煤炭学报,2018,43(9):2510-2517．
LIU Yong,CHEN Changjiang,WEI Jianping,et al. Comparison analysis on the rock breakage pressure induced by abrasive water jets and abrasive gas jets[J].Journal of China Coal Society,2018,43(9):2510-2517.
[2] ZHANG L,YAN X,YANG X,et al. An elastoplastic model of colla-pse pressure for deep coal seam drilling based on Hoek-Brown criterion related to drilling fluid loss to reservoir[J]. Journal of Petroleum Science and Engineering，2015，134:205-213.
[3] 吴彦君. 一种能将磨料回收循环再利用的喷砂设备[P]. 中国专利：CN111496689A,2020-08-07.
[4] 樊晶明,王成勇,王军. 微磨料空气射流加工技术的发展[J]. 金刚石与磨料磨具工程,2005(1):25-30
FAN Jingming,WANG Chengyong,WANG Jun. Development of micro abrasive jet machining technology[J]. Diamond and Abrasives engineering,2005(1):25-30.
[5] LI H,WANG J,KWOK N,et al. A study of the micro-hole geometry evolution on glass by abrasive air-jet micromachining[J]. Journal of Manufacturing Processes，2018，31:156-161.
[6] LIU Y,ZHANG J,WEI J,et al. Optimum structure of a Laval nozzle for an abrasive air jet based on nozzle pressure ratio[J]. Powder Technology,2020,364:343-362.
[7] LIU Y,CHEN C,WEI J,et al. Influence of abrasive hardness on erosion wear of abrasive air jets[J]. Journal of Central South University,2020,52(27):356-371.
[8] KAPIL A,PETER N L. The role of meso-scale structures in rapid gas-solid flows[J]. Journal of Fluid Mechanics,2001,445:151-185.
[9] 何枫,郝鹏飞,张锡文. 超音速射流啸叫模式切换不稳定性[J]. 声学学报,2003,17(2):88-92.
HE Feng,HAO Pengfei,ZHANG Xiwen. Instability of screech switch for under expanded free jet[J]. Acta Acustica,2003,17(2):88-92.
[10] FRENDI A,BROWN M. Flow structures and noise from a super-sonic impinging jet[J]. International Journal of Numerical Methods for Heat and Fluid Flow,2016,26(8):1-20.
[11] 李庠儒,刘年华,刘露菡,等. 冲击射流激波振荡与抑制[J/OL]. 气体物理:1-7[2021-12-28].DOI:10. 19527/j. cnki. 2096-1642. 0898.
LI Xiangru,LIU Nianhua,LIU Luhan,et al. Research and Suppression on shock oscillation of impinging jet[J/OL]. Physics of gases:1-7[2021-12-28].DOI:10. 19527/j. cnki. 2096-1642. 0898.
[12] 徐华舫. 空气动力学基础[M]. 北京：北京航空学院出版社,1987.
[13] ZHAO J,LU X,LIU W,et al. Microscopic characteristics of supersonic natural gas jet in a constant volume chamber under different pressure ratios[J]. Fuel,2020,277(4):118151.
[14] ABEDI M,ASKARI R,SEPAHI-YOUNSI J,et al. Axisymmetric and three-dimensional flow simulation of a mixed compression supersonic air inlet[J]. Propulsion and Power Research，2020，9(1):51-61.
[15] LEI Y,LIU J,QIU T,et al. Gas jet flow characteristic of high-pressure methane pulsed injection of single-hole cylindrical nozzle[J]. Fuel，2019，257:116081.
[16] JUNQUEIRA JUNIOR C,AZEVEDO J,PANETTA J,et al. On the sca-lability of CFD tool for supersonic jet flow configurations[J]. Parallel Computing，2020,93:102620.
[17] 刘海,姚朝晖. 超声速高度欠膨胀冲击射流的大涡模拟[J]. 空气动力学学报,2012,30(1):80-85.
LIU Hai,YAO Zhaohui. Large eddy simulation of superso-nic highly under-expanded impinging Jet [J]. Journal of Aerodynamics,2012,30(1):80-85.
[18] VOGEL A,APITZ I,FREIDANK S,et al. Sensitive high-resolution white-light Schlieren technique with a large dynamic range for the investigation of ablation dynamics[J]. Optics Letters,2006,31(12):1812-1814.
[19] KOLHE P,AGRAWAL A. Density measurements in a supersonic microjet using miniature rainbow schlieren reflectometry[J]. AIAA Journal,2009,47(4):830-838.
[20] 相龙凯,牙宇晨,聂晓康,等. 本生灯法结合纹影技术测量甲烷/空气层流燃烧速度及流场分析[J]. 实验流体力学,2020,34(1):27-34.
XIANG Longkai,YA Yuchen,NIE Xiaokang,et al. Measurement of methane/air laminar combustion velocity and analysis of flow field by Bunsen burner and schlieren technique[J]. Journal of Experiments in Fluid Mechanics,2020,34(1):27-34.
[21] 陈子阳,蒲继雄. 杨氏双缝干涉实验中光谱奇异现象的特性分析[J]. 光子学报,2007,36(4):733-737.
CHEN Ziyang,PU Jixiong. Characteristic analysis of spectral singularity in Yang′s double-slit interference experiment[J]. Acta Photonica Sinica,2007,36(4):733-737.
[22] MEHRA R,DAYANAND M. Application of refractive index mixing rules in binary systems of hexadecane and heptadecane with n -alkanols at different temperatures[J]. Journal of Chemical Sciences,2003,115(2):147-154.
[23] 李炜龙. 纹影仪流场密度定量测量技术研究[D]. 西安：西安工业大学,2017.
[24] 卢义玉,李晓红,廖勇,等. 脉冲磨料射流主要参数对切割性能的影响[J]. 重庆大学学报,2002,25(5):116-118,123.
LU Yiyu,LI Xiaohong,LIAO Yong,et al. Effect of main parameters on cutting of the pulsed abrasive water jet [J]. Journal of Chongqing University,2002,25(5):116-118,123.
[25] BRACEWELL R. The fourier transform and its applications[M]. New York:Mcgraw-hill book compay,1965.
[26] 李晓红,杨林,王建生,等. 自激振荡脉冲射流装置的固有频率特性[J]. 煤炭学报,2000,25(6):82-85.
LI Xiaohong,YANG Lin,WANG Jiansheng,et al. The natural frequency characteristic of the self-excited oscillation pulsed water jet device[J]. Journal of China Coal Society,2000,25(6):82-85.