引　言

空腔外形广泛应用于航空飞行器中, 如飞机起落架、武器舱等. 空腔流动包含剪切层、旋涡、复杂波系等典型流场结构, 在超声速流动中还包含激波. 这些复杂波系之间以及与空腔壁面的相互作用会产生自持振荡而引起强烈的噪声, 噪声声压级可达到170 dB以上^[1], 强噪声不仅会导致环境污染, 还会引起飞行器结构疲劳, 影响飞行安全^[2]. 准确认识空腔流动中复杂波系相互作用产生噪声的物理机理可为研究工程上可行的降噪方法提供理论基础. 为了掌握空腔内部流动结构和振荡频谱特性, 大量学者通过实验和数值模拟研究了不同类型的空腔流动, 相关研究工作可参考相关综述文献[3, 4].

空腔流动的分类与来流条件和空腔长宽比直接相关, 如超声速流动长深比在10以下为开式流动, 长深比大于13为闭式流动, 长深比介于二者之间为过渡式流动^[4]. 开式与闭式流动的主要区别在于, 对于开式流动, 空腔上方的剪切层会流过腔口而撞击到空腔后壁; 而对于闭式流动, 剪切层会撞击空腔底部, 类似于后台阶流动^[4]. 在工程实际中, 空腔的长深比大约在4—10之间, 表现为开式流动^[5]. 开式空腔会因流动的自持振荡产生强噪声, 主要过程为^[6]: 上游边界层脱离空腔前缘, 由于流体的Kelvin-Helmholtz不稳定性, 剪切层逐渐失稳并与空腔后壁面发生碰撞形成反射声波, 反射声波向上游传播到达空腔前缘进一步激发剪切层失稳, 形成反馈回路. 这种声反馈回路机制最早可见于Powell^[7]关于边缘音的研究中. 通过一系列实验, Krishnamurty^[8]发现亚声速和超声速空腔流动会产生强烈的具有离散频率的噪声. Rossiter^[6]进一步拓展了实验参数范围并给出了预测空腔振荡离散频率的半经验公式, 该公式推导基于声速为常数, 因而在马赫数较高时(M > 2)与实验测量偏差明显. Heller等^[9]考虑空腔中的声速变化, 给出了一个改进的半经验预测公式. 部分学者试图给出无经验参数的空腔振荡离散频率预测公式^[10-13], 其中Kerschen等^[12,13]通过将剪切层、空腔中向上下游传播的波系以及空腔两端的散射情况考虑在内提出了一个理论模型, 该模型不包含经验参数且可以预测共振模态的增长与衰减情况.

上述半经验公式和模型基于二维波系结构建立, 然而在实际空腔流动中, 空腔还存在展向流动, 此外空腔侧壁也会对腔内振荡产生影响^[14,15]. Kegerise等^[16]在长深比为2、马赫数为0.4的空腔流动实验中, 观察到腔内压力信号频谱中存在一个未知的低频成分, 通过小波分析发现, 该低频成分会对Rossiter主模态产生调制作用. Brès和Colonius ^[17]基于周期性的展向边界条件数值模拟了三维空腔流动, 他们在速度脉动中同样发现了低频成分, 并认为这种低频成分为三维模态; 同时发现三维模态与Rossiter振荡模态可同时存在, 二者间的相互作用可能会影响空腔流动自持振荡主模态的阶数; 三维模态的波长与空腔深度接近, 比Rossiter振荡模态在幅值上小一个量级. Brès和Colonius ^[17]认为, 三维模态的产生机制与空腔下游壁面附近再循环涡流中的闭合涡流相关的离心不稳定性有关. 因此, 理论模型还须考虑三维模态的影响.

空腔流动产生的噪声在频谱上可能包含多个具有离散频率的声模态, 设计高效空腔噪声控制方案需要厘清这些模态的演化规律. 对于二维空腔流动, Rowley等^[18]采用直接数值模拟研究了亚声速情形, 结果表明一阶、二阶Rossiter模态是共存的. 我们前期针对二维空腔流动的直接数值结果表明^[19—21], 在亚声速情况下, 低阶模态会随着流动稳定逐渐消失; 在超声速情况下, 主模态会产生谐振模态, 主模态与其是共存的. 上述模态演化规律与二维层流这种特定条件有关, 在三维流动中, 声模态的演化特征与上述情形有明显区别. 1998年, Cattafesta等^[22]开展的马赫数0.6的空腔流动实验结果表明, 主要的Rossiter模态之间存在模态切换现象. Gloerfelt等^[23]采用LES数值模拟了展向为周期边界、马赫数0.8的三维空腔流动, 在流动结构中发现了模态切换的证据. Thangamani^[24]实验研究了马赫数1.58的三维超声速空腔流动, 同样发现一阶和二阶主模态之间存在模态切换现象. 我们前期针对三维空腔流动的数值模拟结果同样表明^[25], 空腔振荡的主导模态之间存在切换现象. 但上述研究对模态切换性质的认识仍有不足, 对于空腔流动中各阶模态在流动振荡过程中的非线性相互作用仍需进一步研究.

受Kerschen等^[12,13]工作的启发, 本文通过研究空腔中包含展向波系在内的主要波系结构, 建立相关波系模型. 各部分内容如下: 在第1节中, 区别于以往的二维模型^[6,9], 通过在空腔两端引入展向波系, 同时考虑亚声速与超声速的区别, 分析空腔两端波系的散射过程, 对亚声速与超声速两种情况下的三维空腔波系分别进行建模, 同时分析波系之间的耦合作用; 在第2节中, 介绍马赫数0.9和马赫数1.5空腔流动的实验情况; 在第3节中, 基于马赫数0.9和马赫数1.5空腔流动实验数据, 采用信号分析工具, 结合所建模型分析空腔振荡模态的非线性耦合作用, 同时研究空腔流动中各阶模态的演化规律; 第4节为结论部分.

1.   散射波系模型

1.1   波系模型

1.1.1   亚声速情形

在亚声速情况下, 三维空腔流动中的波系结构示意图如图1所示. 空腔的长、宽、深分别设为L, W, D. 空腔内部包含向上游传播的声波${w_{\text{b}}}$, 其主要是由剪切层撞击后壁导致; 空腔内部向下游传播的声波${w_{\text{f}}}$, 其主要由${w_{\text{b}}}$传播至空腔前端时由前壁面反射形成; 剪切层中的不稳定波${w_{\text{s}}}$, 这是由于流体介质的Kelvin-Helmholtz不稳定性导致的; 空腔口沿剪切层向下游传播的声波${w_{\text{d}}}$以及向上游传播的声波${w_{\text{u}}}$, ${w_{\text{d}}}$的产生与向上游传播的波与空腔前缘相互作用有关, ${w_{\text{u}}}$与剪切层撞击后壁过程直接关联; ${w_{\text{t}}}$表示流动在展向上的波动. 下文中, 在空腔前端的上述波系及相关参数用$\widehat {[\cdot]}$表示, 以区别于在空腔后端的波系及相关参数.

参考Kerschen等^[12,13]的做法, 接下来考虑上述波系在空腔两端的散射情况. 在空腔前端, 考虑波系之间的相互关系有:

$$\left\{ \begin{aligned} &{{{\hat w}_{\text{s}}} ={{\hat C}_{{\text{sb}}}}{{\hat w}_{\text{b}}} + {{\hat C}_{{\text{su}}}}{{\hat w}_{\text{u}}} + {{\hat C}_{{\text{st}}}}{{\hat w}_{\text{t}}},} \\ &{{{\hat w}_{\text{f}}} = {{\hat C}_{{\text{fb}}}}{{\hat w}_{\text{b}}} + {{\hat C}_{{\text{fu}}}}{{\hat w}_{\text{u}}} + {{\hat C}_{{\text{ft}}}}{{\hat w}_{\text{t}}},} \\ & {{{\hat w}_{\text{d}}} ={{\hat C}_{{\text{db}}}}{{\hat w}_b} + {{\hat C}_{{\text{du}}}}{{\hat w}_{\text{u}}} + {{\hat C}_{{\text{dt}}}}{{\hat w}_{\text{t}}}. } \end{aligned} \right.$$

(1)

其中, ${\hat C_{[ \cdot]}}$表示波系之间的散射系数. 类似地, 在空腔后端同样有:

$$\left\{\begin{aligned} &{{w_{\text{u}}} = {C_{{\text{us}}}}{w_{\text{s}}} + {C_{{\text{ud}}}}{w_{\text{d}}} + {C_{{\text{uf}}}}{w_{\text{f}}} + {C_{{\text{ut}}}}{w_{\text{t}}},} \\ &{{w_{\text{b}}} = {C_{{\text{bs}}}}{w_{\text{s}}} + {C_{{\text{bd}}}}{w_{\text{d}}} + {C_{{\text{bf}}}}{w_{\text{f}}} + {C_{{\text{bt}}}}{w_{\text{t}}}. } \end{aligned} \right.$$

(2)

对于空腔前后同一个波, 根据波的增长关系, 可以将波的发展写成复数形式, 对于空腔内的波和剪切层不稳定波有:

$${\hat w_{\text{b}}} = {w_{\text{b}}}{{\text{e}}^{{i} {\alpha _{\text{b}}}L}},\;{w_{\text{s}}} = {\hat w_{\text{s}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{s}}}L}},\;{w_{\text{f}}} = {\hat w_{\text{f}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{f}}}L}}.$$

(3)

对于三维展向波有:

$${\hat w_t} = {\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}},\;{w_{\text{t}}} = {\sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{t}}}W}}.$$

(4)

对于空腔外部的波, 其沿空腔流向以${l^{ - 3/2}}$衰减^{[12, 13]}, l表示离空腔前缘的距离, 因此其形式为

$${\hat w_{\text{u}}} = {w_{\text{u}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{u}}}L}}/{L^{3/2}},\;\;{w_{\text{d}}} = {\hat w_{\text{d}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{d}}}L}}/{L^{3/2}}.$$

(5)

定义变量${\boldsymbol{X}}$为

$${\boldsymbol{X}} = {({\hat w_{\text{s}}},{\hat w_{\text{f}}},{\hat w_{\text{d}}},{w_{\text{b}}},{w_{\text{u}}})^{\text{T}}}.$$

(6)

将(3)式—(5)式代入(1)式和(2)式并整理, 即可得到矩阵形式的方程组: ${\boldsymbol{A}}{\mkern 1 mu} {\boldsymbol{X}} = {\boldsymbol{b}}$, 其中,

$${\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}} { - 1}&0&0&{{{\hat C}_{{\text{sb}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{b}}}L}}}&{{{\hat C}_{{\text{su}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{u}}}L}}/{L^{3/2}}} \\ 0&{ - 1}&0&{{{\hat C}_{{\text{fb}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{b}}}L}}}&{{{\hat C}_{{\text{fu}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{u}}}L}}/{L^{3/2}}} \\ 0&0&{ - 1}&{{{\hat C}_{{\text{db}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{b}}}L}}}&{{{\hat C}_{{\text{du}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{u}}}L}}/{L^{3/2}}} \\ {{C_{{\text{bs}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{s}}}L}}}&{{C_{{\text{bf}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{f}}}L}}}&{{C_{{\text{bd}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{d}}}L}}/{L^{3/2}}}&{ - 1}&0 \\ {{C_{{\text{us}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{s}}}L}}}&{{C_{{\text{uf}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{f}}}L}}}&{{C_{{\text{ud}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{d}}}L}}/{L^{3/2}}}&0&{ - 1} \end{array}} \right) ,$$

(7)

$${\boldsymbol{b}} = - {\Big({\hat C_{{\text{st}}}}{\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}},{\hat C_{{\text{ft}}}}{\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}},{\hat C_{{\text{dt}}}}{\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}},{C_{{\text{bt}}}}{\sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{t}}}W}},{C_{{\text{ut}}}}{\sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{t}}}W}}\Big)^{\text{T}}} .$$

(8)

于是, $\left| {\boldsymbol{A}} \right|{\boldsymbol{X}} = {{\boldsymbol{A}}^*}{\boldsymbol{b}}$, 其中$\left| {\boldsymbol{A}} \right|$表示矩阵${\boldsymbol{A}}$的行列式, ${{\boldsymbol{A}}^*}$是矩阵${\boldsymbol{A}}$的伴随矩阵. 当忽略展向波时, 上述波系为二维情形, 此时${\boldsymbol{b}} = 0$, 此时方程退化为Kerschen等^[12]描述的情形.

显然, 上述方程组解的属性与$\left| {\boldsymbol{A}} \right|$的性质密切相关, 其表达式为

$$\begin{split} \left| {\boldsymbol{A}} \right| =\;& - 1 + {C_{{\text{bs}}}}{{\hat C}_{{\text{sb}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{b}}}{\text{ + }}{\alpha _{\text{s}}})L}} + {C_{{\text{bf}}}}{{\hat C}_{{\text{fb}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{b}}}{\text{ + }}{\alpha _{\text{f}}})L}} + {L^{ - 3/2}}{C_{{\text{bd}}}}{{\hat C}_{{\text{db}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{b}}}{\text{ + }}{\alpha _{\text{d}}})L}} \\& + {L^{ - 3/2}}{C_{{\text{us}}}}{{\hat C}_{{\text{su}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{u}}}{\text{ + }}{\alpha _{\text{s}}})L}} - {L^{ - 3/2}}{C_{{\text{uf}}}}{{\hat C}_{{\text{fu}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{f}}}{\text{ + }}{\alpha _{\text{u}}})L}} + {L^{ - 3}}{C_{{\text{ud}}}}{{\hat C}_{{\text{du}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{u}}}{\text{ + }}{\alpha _{\text{d}}})L}} + {\varGamma _4}, \end{split}$$

(9)

其中, ${\varGamma _4}$表示行列式值的剩余的四阶组合项:

$$\begin{split} {\varGamma _4} =\;& {L^{ - 3}}{{\text{e}}^{{i} ({\alpha _{\text{b}}} + {\alpha _{\text{d}}} + {\alpha _{\text{f}}} + {\alpha _{\text{u}}})L}} {\text{ }}( - {C_{{\text{bf}}}}{C_{{\text{ud}}}}{{\hat C}_{{\text{du}}}}{{\hat C}_{{\text{fb}}}} + {C_{{\text{bd}}}}{C_{{\text{uf}}}}{{\hat C}_{{\text{du}}}}{{\hat C}_{{\text{fb}}}} + {\text{ }}{C_{{\text{bf}}}}{C_{{\text{ud}}}}{{\hat C}_{{\text{db}}}}{{\hat C}_{{\text{fu}}}} - {C_{{\text{bd}}}}{C_{{\text{uf}}}}{{\hat C}_{{\text{db}}}}{{\hat C}_{{\text{fu}}}}) \\& + {L^{ - 3}}{{\text{e}}^{i({\alpha _{\text{b}}} + {\alpha _{\text{d}}} + {\alpha _{\text{s}}} + {\alpha _{\text{u}}})L}} {\text{ }}( - {C_{{\text{bs}}}}{C_{{\text{ud}}}}{{\hat C}_{{\text{du}}}}{{\hat C}_{{\text{sb}}}} + {C_{{\text{bd}}}}{C_{{\text{us}}}}{{\hat C}_{{\text{du}}}}{{\hat C}_{{\text{sb}}}} + {\text{ }}{C_{{\text{bs}}}}{C_{{\text{ud}}}}{{\hat C}_{{\text{us}}}}{{\hat C}_{{\text{su}}}} - {C_{{\text{bd}}}}{C_{{\text{us}}}}{{\hat C}_{{\text{db}}}}{{\hat C}_{{\text{su}}}}) \\& + {L^{ - 3/2}}{{\text{e}}^{{\text{i}}({\alpha _{\text{b}}} + {\alpha _{\text{f}}} + {\alpha _{\text{s}}} + {\alpha _{\text{u}}}{\text{)}}L}} {\text{ }}( - {C_{{\text{bs}}}}{C_{{\text{uf}}}}{{\hat C}_{{\text{fu}}}}{{\hat C}_{{\text{sb}}}} + {C_{{\text{bf}}}}{C_{{\text{us}}}}{{\hat C}_{{\text{fu}}}}{{\hat C}_{{\text{sb}}}} + {\text{ }}{C_{{\text{bs}}}}{C_{{\text{uf}}}}{{\hat C}_{{\text{fb}}}}{{\hat C}_{{\text{su}}}} - {C_{{\text{bf}}}}{C_{{\text{us}}}}{{\hat C}_{{\text{fb}}}}{{\hat C}_{{\text{su}}}}). \end{split}$$

(10)

1.1.2   超声速情形

对于超声速空腔流动, 不考虑激波结构, 其波系结构示意图如图2所示. 在超声速情况下, 空腔中的波系与亚声速情况略有不同, 在空腔上方, 由于流动是超声速的, 声波不能向上游传播, 但向下游传播的声波包含慢声波和快声波, 这里分别记为${w_{{\text{ds}}}}$和${w_{{\text{df}}}}$. 整体波系结构中, 只有空腔中声波${w_{\text{b}}}$向上游传播.

类似于亚声速情况, 分别考虑空腔前后两端的波系散射情况. 在空腔前端有:

$$\left\{\begin{aligned} &{{{\hat w}_{\text{s}}} = {{\hat C}_{{\text{sb}}}}{{\hat w}_{\text{b}}} + {{\hat C}_{{\text{st}}}}{{\hat w}_{\text{t}}}{, }} \\ &{{{\hat w}_{f} } = {{\hat C}_{{\text{fb}}}}{{\hat w}_{\text{b}}} + {{\hat C}_{{\text{ft}}}}{{\hat w}_{\text{t}}}{, }} \\ &{{{\hat w}_{{\text{ds}}}} = {{\hat C}_{{\text{dsb}}}}{{\hat w}_{\text{b}}} + {{\hat C}_{{{\rm{dst}}} }}{{\hat w}_{\text{t}}},} \\ &{{{\hat w}_{{\text{df}}}} = {{\hat C}_{{\text{dfb}}}}{{\hat w}_{\text{b}}} + {{\hat C}_{{\text{dft}}}}{{\hat w}_{\text{t}}}. } \end{aligned} \right.$$

(11)

在空腔后端有:

$${w_{\text{b}}} = {C_{{\text{bs}}}}{w_{\text{s}}} + {C_{{\text{bf}}}}{w_{\text{f}}} + {C_{{\text{bds}}}}{w_{{\text{ds}}}} + {C_{{\text{bdf}}}}{w_{{\text{df}}}} + {C_{{\text{bt}}}}{w_{\text{t}}},$$

(12)

其中, 慢、快声波的演化衰减关系为^[12]

$${w_{{\text{ds}}}} = {\hat w_{{\text{ds}}}}{{\text{e}}^{{\text{i}}{\alpha _{{\text{ds}}}}L}}/{L^{3/2}},\;{w_{{\text{df}}}} = {\hat w_{{{\rm{df}}} }}{{\text{e}}^{{\text{i}}{\alpha _{{\text{df}}}}L}}/{L^{3/2}}.$$

(13)

定义变量${\boldsymbol{X}}$为

$${\boldsymbol{X}} = {({\hat w_{\text{s}}},{\hat w_{\text{f}}},{\hat w_{{\text{ds}}}},{\hat w_{{\text{df}}}},{w_{\text{b}}})^{\text{T}}},$$

(14)

将(3)式、(4)式和(13)式代入(11)式和(12)式并整理, 即可得到矩阵形式的方程组: ${\boldsymbol{A}}{\mkern 1 mu} {\boldsymbol{X}} = {\boldsymbol{b}}$, 其中,

$$\begin{split} {\boldsymbol{A}} =\left( {\begin{array}{*{20}{c}} { - 1}&0&0&0&{{{\hat C}_{{\text{sb}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{b}}}L}}} \\ 0&{ - 1}&0&0&{{{\hat C}_{{\text{fb}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{b}}}L}}} \\ 0&0&{ - 1}&0&{{{\hat C}_{{\text{dsb}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{b}}}L}}} \\ 0&0&0&{ - 1}&{{{\hat C}_{{\text{dfb}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{b}}}L}}} \\ {{C_{{\text{bs}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{s}}}L}}}&{{C_{{\text{bf}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{f}}}L}}}&{\dfrac{{{C_{{{\rm{bds}}} }}{{\text{e}}^{{\text{i}}{\alpha _{{\text{ds}}}}L}}}}{{{L^{3/2}}}}}&{\dfrac{{{C_{{\text{bdf}}}}{{\text{e}}^{{\text{i}}{\alpha _{{df} }}L}}}}{{{L^{3/2}}}}}&{ - 1} \end{array}} \right) , \end{split}$$

(15)

$${\boldsymbol{b}} = - \left({\hat C_{{\text{st}}}}{\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}},{\hat C_{{\text{ft}}}}{\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}}, {\hat C_{{\text{dst}}}}{\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}},{\hat C_{{\text{dft}}}}{\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}},{C_{{\text{bt}}}}{\sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}}\right)^{\text{T}} ,$$

(16)

行列式值$\left| {\boldsymbol{A}} \right|$为

$$\left| {\boldsymbol{A}} \right| = - 1 + {C_{{\text{bs}}}}{{\hat C}_{{\text{sb}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{b}}}{\text{ + }}{\alpha _s})L}} + {C_{{\text{bf}}}}{{\hat C}_{{\text{fb}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{b}}}{\text{ + }}{\alpha _{\text{f}}})L}} + {L^{ - 3/2}}{C_{{\text{bds}}}}{{\hat C}_{{\text{dsb}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{b}}}{\text{ + }}{\alpha _{{\text{ds}}}})L}} + {L^{ - 3/2}}{C_{{\text{bdf}}}}{{\hat C}_{{\text{dfb}}}}{{\text{e}}^{{\text{i}}({\alpha _{\text{b}}}{\text{ + }}{\alpha _{{\text{df}}}})L}}.$$

(17)

不难发现, 由于超声速情形在空腔两端的散射特征较亚声速情形简单, 因此, 方程组的系数矩阵和对应的行列式值也更为简洁.

1.2   模型分析

在Rossiter^[6]得到的半经验公式中, 图1和图2所示的波系只考虑了${w_{\text{s}}}$和${w_{\text{b}}}$, 方程组变为${\boldsymbol{A}}{\mkern 1 mu} {\boldsymbol{X}} = 0$. 其中, 矩阵${\boldsymbol{A}}$简化为

$${\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}} { - 1}&{{{\hat C}_{{\text{sb}}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{b}}}L}}} \\ {{C_{{\text{bs}}}}{{\text{e}}^{{\text{i}}{{\text{α }}_{\text{s}}}L}}}&{ - 1} \end{array}} \right) .$$

(18)

因为方程组有非零解, 于是必有$\left| {\boldsymbol{A}} \right| = 0$, 即

$$\left| {\boldsymbol{A}} \right| = - 1 + {C_{{\text{bs}}}}{\hat C_{{\text{sb}}}}{{\text{e}}^{{\text{i(}}{\alpha _{\text{b}}}{\text{ + }}{\alpha _{\text{s}}})L}} = 0.$$

(19)

利用复数欧拉公式, 可以得到:

$${\text{Arg}}[{C_{{\text{bs}}}}] + {\text{Arg}}[{\hat C_{{\text{sb}}}}] + {Re} [{\alpha _{\text{b}}} + {\alpha _{\text{s}}}]L = 2n{\text{π}} ,$$

(20)

其中, Arg[▪]表示复数的幅角, Re[▪]表示复数的实部. (20)式可以进一步表示为

$$L = 2{\text{π}}\frac{{m - \dfrac{{{\text{Arg}}[{C_{{\text{bs}}}}] + {\text{Arg}}[{{\hat C}_{{\text{sb}}}}]}}{{2{\text{π}}}}}}{{{Re} [{\alpha _{\text{b}}} + {\alpha _{\rm{s}}}]}} ,$$

(21)

其中${\alpha _{\text{s}}}$是剪切层不稳定波的波数, 设剪切层的对流速度为${U_c}$, 即有${\alpha _{\text{s}}} = \omega /{U_{\text{c}}}$. ${\alpha _{\text{b}}}$是空腔后壁面向上游传播的声波的波数, 其传播速度近似为声速, 因此, ${\alpha _{\text{b}}} = \omega /a = \omega \cdot M/U$. 将${\alpha _{\text{s}}}$和${\alpha _{\text{b}}}$代入(21)式并化简可以得到:

$${\text{S}}{{\text{t}}_n} = \frac{{\omega L}}{{2{\text{π}}U}} = \frac{{{f_n}L}}{U} = \frac{{n - \hat \gamma }}{{M + 1/\hat k}},\;n = 1,2,\cdots ,$$

(22)

其中,

$$\left\{ {\begin{aligned} \hat \gamma =\;& {{\left( {{\text{Arg}}[{C_{{\text{bs}}}}] + {\text{Arg}}[{{\hat C}_{{\text{sb}}}}]} \right)} \big/{(2{\text{π}}),}} \\ \hat k =\;& {U_{\text{c}}}/U. \end{aligned}} \right.$$

(23)

(22)式和(23)式是Kerschen等^[12]通过空腔两端的散射模型获得的空腔振荡频率预测公式, 可以看出该公式与Rossiter半经验公式^[6]形式上完全一致. 在(23)式中, 函数$\hat \gamma$和$\hat k$与具体的空腔尺寸和来流条件有关, 在实际中一般变化较小, 因此可以近似取为常数. 在Rossiter半经验公式中, 一般取$\hat \gamma = 0.25$, $\hat k$为剪切层相对对流速度, 一般取$\hat k = 0.57$. 在我们的二维层流空腔流动高精度数值模拟中^[19-21], 获得的$\hat k = 0.55$, 与常用值较为接近. 在实际三维空腔流动中, 由于湍流、侧壁等的影响, 其值可能会有所区别, 可以根据实际测量值进行估计, 如在Ahuja和Mendoza^[14]以及Kegerise^[26]的实验中, 得到的估计值为$\hat k = 0.66$.

当进一步考虑三维展向波系时, 右端项${\boldsymbol{b}}$简化为:

$${\boldsymbol{b}} = - {({\hat C_{{\text{st}}}}{\hat \sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{{\hat \alpha }_{\text{t}}}W}},{C_{bt}}{\sigma _{\text{t}}}{{\text{e}}^{{\text{i}}{\alpha _{\text{t}}}W}})^{\text{T}}}.$$

(24)

即有:

$${\boldsymbol{X}} = \left( {\begin{aligned} {\frac{{{{\hat C}_{{\text{bt}}}}{{\hat C}_{{\text{sb}}}}{\sigma _{\text{t}}}{{\text{e}}^{{\text{i}}\left( {L{\alpha _{\text{b}}} + W{\alpha _{\text{t}}}} \right)}} + {{\hat C}_{{\text{st}}}}{{\hat \sigma }_{\text{t}}}{{\text{e}}^{{\text{i}}W{{\hat \sigma }_{\text{t}}}}}}}{{ - 1 + {C_{{\text{bs}}}}{{\hat C}_{{\text{sb}}}}{{\text{e}}^{{\text{i}}L\left( {{\alpha _{\text{b}}} + {\alpha _{\text{s}}}} \right)}}}}} \\ {\frac{{{{\text{e}}^{{\text{i}}\left( {L{\alpha _{\text{s}}} + W{{\hat \sigma }_t}} \right)}}{C_{{\text{bs}}}}{x_{{\text{st}}}}{{\hat \sigma }_{\text{t}}} + {{\text{e}}^{{\text{i}}W{\alpha _{\text{t}}}}}{{\hat C}_{{\text{bt}}}}{\sigma _{\text{t}}}}}{{ - 1 + {{\text{e}}^{{\text{i}}L\left( {{\alpha _{\text{b}}} + {\alpha _{\text{s}}}} \right)}}{C_{{\text{bs}}}}{{\hat C}_{{\text{sb}}}}}}} \end{aligned}} \right) .$$

(25)

此时方程组${\boldsymbol{A}}{\mkern 1 mu} {\boldsymbol{X}} = {\boldsymbol{b}}$决定了三维展向波系将与Rossiter模态产生调制. 尽管三维展向模态频率和幅值较低, 但在特定情况下, 三维模态仍然会与主模态产生非线性作用. 在Neary和Stephanoff^[27]基于水介质的层流空腔实验中, 他们发现了三维模态与主模态调制产生的谐频. 在Kegerise等^[16]马赫数为0.4的空腔流动实验中, 三维模态对Rossiter主模态产生了调制作用. 同时, 空腔宽度参数$W$也会影响不同离散频率组分的幅值^[3].

当考虑空腔中的其余波系时, 方程组${\boldsymbol{A}}{\mkern 1 mu} {\boldsymbol{X}} = {\boldsymbol{b}}$的形式变得较为复杂, 矩阵${\boldsymbol{A}}$的行列式值如(9)式和(17)式所示. 此时, 不同波系之间将产生非线性的相互作用, 尽管部分散射过程波系相互作用的系数值较小, 但仍将产生相互调制作用, 这种调制过程可能产生其余频率成分. Kegerise等^[16]通过对空腔中测得的压力信号采用高阶谱分析发现, 在特定的流动条件下, 不同Rossiter模态之间会产生非线性作用, 这种非线性作用将导致产生不同于Rossiter模态的谐频成分, 频率大小与两个Rossiter模态频率的差值相关.

2.   三维空腔实验

2.1   实验设置

空腔模型C201由中国空气动力研究与发展中心高速空气动力研究所设计^[28], 模型示意图如图3所示. 空腔的长(L)、宽(W)、深(D)分别为200, 66.7, 33.3 mm, 即满足L∶W∶D = 6∶2∶1, 模型的其余尺寸见图3标注. 在实验中, 亚声速和超声速的前缘设计略有区别, 具体可参考引用文献中的说明^[28].

实验在中国空气动力研究与发展中心高速空气动力研究所亚跨超声速风洞中完成, 实验包含亚声速和超声速两种来流条件, 马赫数分别为0.9和1.5, 其余参数条件如表1所列.

2.2   实验验证

实验中的脉动压力测量点设置在空腔底部中截面上, 下面通过重复性实验, 根据壁面压力信号总声压级验证实验测量的可靠性. 总声压级(overall sound pressure level, OASPL)计算公式如下:

$${\text{OASPL}} = 20{\text{lo}}{{\text{g}}_{10}}\left( {\frac{{{p_{\text{e}}}}}{{{P_{{\text{ref}}}}}}} \right),$$

(26)

其中,

$${p_{\text{e}}} = \sqrt{\dfrac1 N \displaystyle\sum\limits_{i=1}^N (p-\bar p)^2},$$

(27)

$\bar p$为压力平均值, 空气中参考声压一般取${p_{{\text{ref}}}} = 2.0 \times {10^{ - 5}}{\text{Pa}}$, 对应正常人耳对 1 kHz声音的可听阈声压. 对于马赫数1.5, 两次实验的对比结果如图4所示, 两次测压最大差别为1.1 dB, 表明实验重复性较好, 实验系统可靠.

3.   模态演化特性分析

3.1   频谱调制特性

声压级的频域分布通过压力脉动的声功率谱密度(pressure spectral density, PSD)进行计算, 公式如下,

$${\text{SPL}} = 10{\log _{10}}\left(\frac{{{\text{PSD}}}}{{P_{{\text{ref}}}^2}}\right) .$$

(28)

对于当前两个马赫数下, 空腔前壁测点的压力信号声压级频谱计算结果如图5所示. 两种马赫数的主导频率如表2所列. 在(23)式中, 参数$\hat \gamma$和$\hat k$与具体的空腔尺寸和流动条件相关, 这里可以通过得到的表2中的频谱数据拟合得到. 对于马赫数0.9, 通过线性拟合, 可以得到$\hat r = 0.45$, $\hat k = 0.68$; 类似地, 对于马赫数1.5, 可以得到$\hat r = 0.40$, $\hat k = 0.86$. 将获得的参数代入(22)式即可得到离散频率预测公式, 公式预测值见图5中灰色竖线所示. 可以看到, 两种马赫数下, 空腔流动的振荡主导模态都为Rossiter二阶模态, 同时预测的各阶模态频率位置与压力信号频谱峰值符合良好.

从图5可以观察到, 两个马赫数下的频谱中都存在一个低频成分, 对马赫数0.9, 这个低频频率为$f = 134.3{\text{ Hz}}$, 对马赫数1.5, 这个低频频率为$f = 122.2{\text{ Hz}}$. 根据Brès和Colonius^[17]的研究, 这种低频成分与三维模态有关, 其波长与空腔深度接近, 因此建议采用空腔深度作为其特征长度. 对于马赫数0.9, 低频成分的无量纲频率为${\text{S}}{{\text{t}}_D} = {{f_l}D}/ U = 0.016$, 对于马赫数1.5, 低频成分的无量纲频率为${\text{S}}{{\text{t}}_D} = {{{f_l}D} /U} = 0.01$. 在Kegerise等^[16]马赫数0.4的空腔流动实验中, 观察到的低频成分的无量纲频率为${\text{S}}{{\text{t}}_D} = {{{f_l}D} / U}= 0.011$, 与当前两个马赫数下观测到的低频成分无量纲频率较为接近.

在第2节的分析中, 空腔中出现的主要波系的主导频率可以通过方程组确定, 不同波系之间可能产生相互调制, 这种调制过程可能产生区别于主导Rossiter模态之外的其余频率成分. 这种现象可以在本文两个实验的频谱结果中观察到. 在图5(a)所示的马赫数0.9的频谱结果中, 在圆圈圈出的位置可以观察到一个显著的区别于主频的频率: ${f_{\text{s}}} = {\text{3613 Hz}}$. 根据表2中马赫数0.9对应的结果可以发现${f_{\text{s}}} \approx {f_3} + {f_4} = 3612{\text{ Hz}}$. 即频率${f_{\text{s}}}$是由Rossiter三阶模态与Rossiter四阶模态调制产生的. 在图5(b)所示的马赫数1.5的频谱结果中, 在圆圈圈出的位置可以观察到有两个显著的区别于主频的频率${f_{{\text{s1}}}} = {\text{2526 Hz}}$和${f_{{\text{s2}}}} = {\text{3314 Hz}}$. 根据表2中马赫数1.5对应的结果可以发现${f_{{\text{s1}}}} \approx {f_2} + {f_2} = {\text{2528 Hz}}$, ${f_{{\text{s2}}}} \approx {f_2} + {f_3} = 3306{\text{ Hz}}$. 即频率${f_{{\text{s1}}}}$是由Rossiter二阶模态与自身调制产生的, 频率${f_{{\text{s2}}}}$是由Rossiter二阶模态与Rossiter三阶模态调制产生的. 这里仅列出了个别调制产生的幅值显著的谐频, 接下来采用高阶谱分析进一步分析这种调制作用.

从上述分析可以看出, 对于马赫数1.5, 其Rossiter二阶主模态的幅值显著高于其余振荡模态, 并可与多个其余模态调制, 产生区别于主模态的谐频. 这里采用双谱分析(bispectral analysis)来分析调制过程. 对于一个离散时间信号$x(t)$, 其对应的频域设为${X_k}$, 即满足:

$$x(t) = \mathop \sum \limits_{k = - \infty }^\infty {X_k}{{\text{e}}^{ - {\text{i}}{\omega _k}t}}$$

(29)

和

$${X_k} = \mathop {\lim }\limits_{k \to \infty } \frac{1}{T}\mathop \int \nolimits_{ - T/2}^{T/2} {{\text{e}}^{{\text{i}}{\omega _k}t}}x(t){\mkern 1mu} {\text{d}}t ,$$

(30)

其中, ${\omega _k} = 2{\text{π}}k/T$, $T$是信号$x(t)$的长度. 相应的三阶累计谱为

$$B(k,l) = E[{X_k}{X_l}X_{k + l}^*] ,$$

(31)

其中, E[▪]表示期望值. 双谱分析相关系数bic定义为^[29]

$${\text{bi}}{{\text{c}}^2}(k,l) = \frac{{|B(k,l){|^2}}}{{E[|{X_k}{X_l}{|^2}]E[|{X_{k + l}}{|^2}]}} .$$

(32)

利用Cauchy-Schwarz不等式容易证明上述相关系数满足$0 \leqslant {\text{bic}} \leqslant 1$. 上述定义中的双谱相关系数具有对称性质. 特别是, 当频率${\omega _k}$和频率${\omega _l}$存在非线性耦合时, ${\text{bic}} \to 1$; 而当频率${\omega _k}$和频率${\omega _l}$相互独立时, ${\text{bic}} \to 0$; 当bic为介于0—1之间的其余数值时, 其值的大小表征频率${\omega _k}$和频率${\omega _l}$非线性作用产生频率${\omega _{k + l}}$占实际频率${\omega _{k + l}}$幅值的比例. 于是可以利用上述性质来分析空腔振荡中不同模态的非线性耦合情况.

马赫数1.5空腔前面测量信号的双谱分析系数计算结果如图6所示, 图6(a)为整体结果, 图6(b)为局部区域放大情况, 结果关于斜中心线是对称的. 由于二阶Rossiter主模态较强, 这里主要关注与其相关的非线性作用情况. 从图6(a)可以看出, 与主频${f_2}$相关的非线性作用十字区域相较于其余区域十分明显. 在坐标$({f_2}, {f_2})$处, bic = 0.94, 即主模态自相关性很强, 由此会产生相关频率$2{f_2}$, 这与上文中结论一致, 从该频率成分在图5(b)中的频谱结果可以看到, 幅值十分显著. 从图6(b)还可以看到, 在坐标$({f_2}, f = 1361{\text{ Hz}})$处双谱相关系数较大: bic = 0.85, 意味着将产生频率为${f_2} + f = 2625{\text{ Hz}}$的成分, 该频率成分见图5(b)中的标注.

3.2   模态演化规律

在3.1节中, 使用快速傅里叶变换(FFT)等获得了主要振荡模态的频谱特征, 但FFT缺乏时频定位功能, 且无法获得非平稳信号的频率演化信息. 从前文论述中可知, 在三维空腔流动中, 流动振荡的各阶模态之间存在模态切换现象, 因此在接下来的分析中, 采用时频分析方法来分析其压力脉动特征. 采用连续小波变换(continuous wavelet transform, CWT)来分析压力脉动信号, 小波变换相较于短时傅里叶变换等时频分析方法, 可以实现随分辨率变化而自动调节分析带宽, 且计算过程存储量小, 适合分析长时间尺度的信号, 其具体变换公式如下^[30]:

$$\begin{split} {W_x}(a,b) = \;&\frac{1}{{\sqrt {|a|} }}\int_R {x(t){\mkern 1mu} {\psi ^*}\left( {\frac{{t - b}}{a}} \right)} {\mkern 1mu} {\text{dt}} \\ = \;&\int_R {x(t){\mkern 1mu} \psi _{a,b}^*(t)} {\mkern 1mu} {\text{d}} t \\ =\;& \left\langle {x(t),\;{\psi _{a,b}}(t)} \right\rangle , \end{split}$$

(33)

其中, $\varPsi$为小波函数. 本文选取Morlet小波, 形式为

$$\varPsi (t) = {{\text{e}}^{ - {t^2}/2}}{{\text{e}}^{{\text{i}}{\varOmega _0}t}} .$$

(34)

小波变换计算结果如图7所示, 图中反映了压力脉动中不同频率成分随时间的演化情况. 在前面关于二维层流空腔流动高精度数值模拟的腔内压力信号的时频分析结果中显示, 主导模态随时间演化过程中幅值稳定, 当前的时频结果与二维简单情形有明显的不同. 根据表2中的数值, 从图7可以看出, 在两个马赫数下, 主要振荡模态的幅值随时间有明显的改变, 存在模态切换现象, 但Rossiter 2模态占主导时幅值要更高, 与前面的频谱分析结果一致. 为了进一步分析主要振荡模态随时间的演化特征, 以马赫数0.9为例, 将一阶模态${f_1} = 330{\text{ Hz}}$和二阶模态${f_2} = 940{\text{ Hz}}$从小波分析结果中单独截取出来, 结果如图8所示. 需要说明的是, 从前面的理论模型分析中可知, 由于不同波系之间存在非线性作用, 主要振荡模态的频率在随时间演化过程中可能会发生微小的改变, 由于变化很小, 这里将其忽略. 从图8可以看出, 模态${f_1}$和${f_2}$的幅值随时间变化剧烈, 对截取结果进行FFT分析结果如图9所示. 从图9可以看出, 模态成分${f_1}$和${f_2}$的幅值随时间变化的波动频率较低, 但没有显著的主频. 进一步地, 使用概率密度估计来计算模态${f_1}$和${f_2}$随时间演化信息的概率密度函数分布, 结果如图10所示. 从图10可以看出, 两个模态随时间演化的概率密度分布形态与正态分布形态较为接近, 即模态演化整体表现出随机性.

4.   结　论

本文通过研究空腔中包含展向波系在内的主要波系结构, 分析不同波系在空腔两端的散射情况, 同时考虑亚声速和超声速流动下波系传播的差异, 分别建立了针对亚声速和超声速空腔流动的三维波系模型. 三维波系模型包含了空腔中不同波系之间的非线性相互作用, 这种非线性作用可导致产生不同于Rossiter模态的其余频率成分. 当仅保留剪切层不稳定波及其撞击后壁引起的反射声波时, 模型退化为二维简单情况, 此时方程表征空腔流动主要振荡模态—Rossiter模态.

基于三维空腔流动实验测量的马赫数0.9和马赫数1.5的腔内压力信号数据, 通过计算其声压级频谱, 获得了两种典型流动条件下, 空腔自持振荡的主要频谱特征. 根据获得的主要振荡频率, 对所建模型中的参数进行了线性估计, 模型预测的各阶模态频率位置与压力信号频谱峰值符合良好. 对压力信号采用双谱分析, 结果表明, 振荡主模态之间会产生非线性作用, 这种非线性作用产生了幅值较高的谐频, 与频谱分析结果一致. 采用连续小波变换方法对压力信号进行了时频分析, 结果表明, 主要的振荡模态之间存在模态切换现象, 且模态切换呈低频特征, 整体表现出随机性特性.

需要注意的是本文的空腔波系结构模型未考虑有舱门的情况, 舱门的存在会对空腔口的波系产生反射作用, 且舱门在不同角度时, 这种反射作用会有所差别. 针对上述情形, 仍需结合仿真和实验开展进一步研究.