0 引　言

自然界和工业生产过程中存在大量声波与剪切层流动相互作用的物理现象。1858年Leconte^[1]发现音乐会上的蜡烛会随着乐器的声音“随音起舞”，即声波对燃烧射流的影响。在声学风洞实验中，声波穿过开口声学风洞与周围流体形成的剪切层，会产生折射、反射现象，从而导致测量结果失真^[2]。因此，声波与剪切层的相互作用的机理研究可用于提升人们对燃烧、湍流涡环的控制能力，亦或是准确的声源定位技术，具有重要的科学与工程价值。

声激励对流动发展的影响存在于大量流动稳定性研究之中。1878年Rayleigh^[3]研究了射流的稳定性问题，探讨了流动失稳的数学原理。在大量关于射流或混合层稳定性研究的实验^[4-6]中，声波被用于充当初始小扰动。Sato和Sakao^[7]使用纵向声波激励平面射流，通过流场结构变化确定出二维平面射流失稳的临界雷诺数；Ho和Huang^[8]通过实验发现使用与混合层不稳定波频率呈一定数值关系的声波激励混合层，可以控制流场中涡的融合，并依此提出亚谐频反馈回路机制。在一些实验^[9-10]中为凸显出相干结构的演化规律，研究人员对流场施加弱声波激励使其结构变得清晰可见。

尽管声波影响流动的现象被广泛研究，但流动稳定性研究的关注重点是不稳定波的演化，声波只被简单当作初始扰动引入，而缺乏更精细的过程描述；声波与剪切层失稳波的耦合也在大部分工作中被忽略。Tam的报道^[11-12]指出关于声波如何激励起不稳定波的机理问题是感受性问题，通过建立波束与二维剪切层相互作用的数学模型，求解得到声波激励剪切层失稳的两个条件：1）声波与不稳定波频率需要匹配；2）二者在流向上的相速度需要匹配。吴介之^[13]亦在文中指出感受性与稳定性的不同之处，即感受性是激励声波在流场中的表现。目前关于声波在剪切层中的演化规律研究较少。

另一方面，流动对声传播的影响被广泛关注。许多学者针对声波穿过剪切层的物理过程开展了大量研究。从Miles^[14]研究穿过有相对运动速度流体的声波传播规律开始，研究者对波动方程做了诸多修正并提出多种模型，如Amiet^[15]创建的剪切层修正理论以及Padios等^[16]修正了对流引起的声源偏移。在数值模拟方面，张雪等^[17]通过求解线化Euler方程数值模拟了二维单极子在剪切层中的声传播，对比验证了Amiet理论修正剪切层声折射现象的准确性。尽管线性Euler方程可以捕捉声传播的过程，但由于其未考虑黏性和非线性项，预测的剪切层Kelvin-Helmholtz（K-H）不稳定波的发展并不符合物理规律，也就无法精确描述声波与剪切层之间的相互作用。更进一步地，关于K-H不稳定性现象的研究是流体力学中的一个经典问题^[18]，尽管近年来其在可压缩流^[19-20]、磁流体^[21]、多相流^[22]等环境下有长足的进展，但对于声激励条件下剪切流动失稳机理尚不十分明确。即，从剪切层对声传播干扰方面来看，现有理论只注重了剪切层外的传播规律及接收信号的变化，只有剪切层对声传播单方面的影响，而较少关注声波在剪切层内部是如何与流动相互作用（声波激励起剪切层中复杂的涡系结构）的。

格子Boltzmann方法（lattice Boltzmann method，LBM）^[23]是20世纪90年代发展起来的一种介观计算流体力学方法，备受关注。它具有编程简单、并行性好等特征，适用于大规模计算。研究表明，它能够很好地对具有复杂固壁边界、湍流等的流动现象进行准确预测。在1998年Buick等^[24]采用LBM成功地获得了声信号传播的线性与非线性行为，而关于LBM方法在声学领域适用性的论证源自于Lallemand和Luo^[25]对LBM数值耗散、色散特性的分析。Marié等^[26]进一步分析了多种LBM数值格式并指出：LBM可以近乎理想地模拟声学现象。Najafi-Yazdi和Mongeau^[27]将完美匹配层方法（perfect match layer，PML）添加到LBM框架中，这是一种效果较好的声学无反射边界条件。上述工作为LBM方法在声学模拟领域的发展奠定了坚实的基础，使其得以用于很多流场模拟领域。以LBM为核心求解方案的流体力学商业软件PowerFlow，已成功地应用于轿车气动声音的分析。Viggen^[28]研究了流体的剪切黏性、体黏性对声传播的影响，其研究进一步确定了LBM对声传播捕捉的精确性。对中高雷诺数流动^[29]、非线性声学现象^[30]以及涡运动噪声及其与壁面相互作用^[31]的研究也证实了LBM对流动噪声问题的模拟具有令人满意的效果。

综上所述，声波穿过剪切层的前期研究大都没有考虑双向耦合。剪切层稳定性受声场影响，并进一步给声传播提供复杂流动背景，改变声传播特征，其内在规律有待研究。因此本文采用LBM直接数值模拟了不同波数声波穿过剪切层的传播规律，以及在其影响下剪切层K-H波的激励与发展过程。本文后续部分中，将首先介绍所采用的数值方法和计算模型，然后结合剪切层稳定性理论，展示和讨论剪切层对声传播的影响以及声波激励剪切层流动失稳的机理，最后给出结论。本文相关研究成果可用于声学风洞声信号解读、指导流动控制等技术的发展，具有重要的学术意义和工程应用价值。

1 数值方法与计算模型 1.1 数值方法

LBM基于经典Boltzmann方程，是一种基于气体动理论的介观方法。采用指标法以及Einstein求和约定表达（下文相同）的控制方程为：

$ \frac{\partial f}{\partial t}+{\xi }_{i}{\partial }_{i}f = \mathrm{\varOmega }\left(f\right) $

(1)

其中， $ f(\boldsymbol{x},\boldsymbol{\xi },t) $ 是粒子速度分布函数，表征了空间为 $ \boldsymbol{x} $ 、时间为 $ t $ 、速度为 $ \boldsymbol{\xi } $ 的粒子的数密度； $\mathrm{\varOmega }\left(f\right)$ 表征了粒子碰撞影响，是一个关于速度空间的二次积分。Boltzmann方程刻画了粒子输运过程所遵循的一般规律，但该积分微分方程求解困难。基于弱非平衡假设、粒子碰撞的熵增结果，以及质量和动量守恒的约束条件，碰撞函数可以简化成一个关于 $ f $ 和 ${f}^{{\rm{eq}}}$ 的算子：

$ \mathrm{\varOmega }\left(f\right) = -\frac{f-{f}^{{\rm{eq}}}}{\tau } $

(2)

其中 ${f}^{{\rm{eq}}}$ 是平衡态分布函数：

$ {f}^{{\rm{eq}}} = \rho {\left(\frac{1}{2\text{π} RT}\right)}^{D/2}\mathrm{exp}\left[-{\left|\boldsymbol{v}\right|}^{2}/\left(2RT\right)\right] $

(3)

式中， $ \boldsymbol{v} $ 表征了粒子的无序热运动， $ \boldsymbol{v}\left(\boldsymbol{x},t\right) = \boldsymbol{\xi }\left(\boldsymbol{x},t\right)-\boldsymbol{u}(\boldsymbol{x},t) $ ， $ \boldsymbol{u}\left(\boldsymbol{x},t\right) = \int \boldsymbol{\xi }f{d}^{3}\boldsymbol{\xi } $ 为宏观流动速度； $ \;\rho $ 指粒子密度； $ T $ 为温度； $ R $ 为气体常数； $ D $ 为维度。碰撞算子中 $ \tau $ 为松弛系数，表征了 $ f $ 向平衡态演化的速率，其决定了气体黏性。这种碰撞被称为BGK（Bhatnagar-Gross-Krook）碰撞算子。它因形式简单、物理意义明朗而被广泛采用。多种数学分析^[32-33]显示Boltzmann方程是与Navier-Stokes方程相容的，并具有更高精度。

不同于传统计算流体力学方法，LBM首先对速度相空间进行离散（ $ \boldsymbol{\xi }\to {\boldsymbol{c}}_{i} $ ），得到速度离散Boltzmann方程（discrete velocity Boltzmann equation）：

$ {\partial }_{t}{f}_{i}+{c}_{i\alpha }{\partial }_{\alpha }{f}_{i} = -\frac{1}{\tau }({f}_{i}-{f}_{i}^{{\rm{eq}}}) $

(4)

文中，希腊字母（ $\alpha 、\beta$ 等）表示笛卡尔坐标分量，指标 $ i $ 表示离散速度的方向。 $ {f}_{i} $ 指经速度离散后的分布函数， $ {\boldsymbol{c}}_{i} $ 是离散后的粒子速度。Shan等^[33]将公式（1）在Hermit基上进行投影，推导获得任意阶平衡态分布函数，本文采用四阶精度形式^[34-35]：

$ {f}_{i}^{{\rm{eq}}}\approx {w}_{i}\rho \left[1+\frac{{c}_{i\alpha }{u}_{\alpha }}{{c}_{s}^{2}}+\frac{{A}^{\left(2\right)}}{2{c}_{s}^{4}}+\frac{{A}^{\left(3\right)}}{2{c}_{s}^{6}}++\frac{{A}^{\left(4\right)}}{4{c}_{s}^{8}}+O\left({u}^{5}\right)\right] $

(5)

$ {A}^{\left(2\right)} = {H}_{i\alpha \beta }^{\left(2\right)}{u}_{\alpha }{u}_{\beta } $

(6)

$\begin{split} {A}^{\left(3\right)} =& {H}_{ixxy}^{\left(3\right)}{u}_{x}^{2}{u}_{y}+{H}_{ixxz}^{\left(3\right)}{u}_{x}^{2}{u}_{z}+{H}_{ixyy}^{\left(3\right)}{u}_{y}^{2}{u}_{x}+{H}_{ixzz}^{\left(3\right)}{u}_{z}^{2}{u}_{x} +\\&{H}_{iyzz}^{\left(3\right)}{u}_{z}^{2}{u}_{y}+{H}_{iyyz}^{\left(3\right)}{u}_{y}^{2}{u}_{z}+{H}_{ixyz}^{\left(3\right)}{u}_{x}{u}_{y}{u}_{z} \end{split}$

(7)

$ \begin{split}&{A}^{\left(4\right)} = {H}_{ixxyy}^{\left(4\right)}{u}_{x}^{2}{u}_{y}^{2}+{H}_{ixxyy}^{\left(4\right)}{u}_{x}^{2}{u}_{y}^{2}+{H}_{ixxzz}^{\left(4\right)}{u}_{x}^{2}{u}_{z}^{2}+{H}_{iyyzz}^{\left(4\right)}{u}_{z}^{2}{u}_{y}^{2} +\\&2\left({H}_{ixyzz}^{\left(4\right)}{u}_{x}{u}_{y}{u}_{z}^{2}+{H}_{ixyyz}^{\left(4\right)}{u}_{x}{u}_{y}^{2}{u}_{z}+{H}_{ixxyy}^{\left(4\right)}{u}_{x}^{2}{u}_{y}^{2}+{H}_{ixxyz}^{\left(4\right)}{u}_{x}^{2}{u}_{y}{u}_{z}\right) \end{split}$

(8)

其中， $ {w}_{i} $ 是不同离散速度的权重系数， $ {c}_{s} $ 是格子声速，这两个量取决于离散速度 $ {\boldsymbol{c}}_{i} $ 的选取。公式中 $ {\boldsymbol{H}}^{\left(n\right)} $ 是 $ n $ 阶Hermit多项式，具有张量属性。采用高阶平衡态函数，LBM具有更好的稳定性，能有效避免相对马赫数（ $M {a}_{1}-M {a}_{2}$ ）过大造成的剪切流动模拟发散。

本文模拟的是二维问题（D = 2），选用D2Q9速度模型，其中离散粒子速度定义为：

$ {\boldsymbol{c}}_{i} = \left\{\begin{array}{ll}\left(\mathrm{0,0}\right)，& i = 0\\ \left[\mathrm{sin}\left(\dfrac{i-1}{2}\text{π} \right),\mathrm{cos}\left(\dfrac{i-1}{2}\text{π} \right)\right]，& i = 1～4\\ \sqrt{2}\left[\mathrm{sin}\left(\dfrac{2i-1}{4}\text{π} \right),\mathrm{cos}\left(\dfrac{2i-1}{4}\text{π} \right)\right]，& i = 5～8\end{array}\right. $

(9)

速度权重为：

$ {w}_{i} = \left\{\begin{array}{ll}4/9，& i = 0\\ 1/9，& i = 1～4\\ 1/36，& i = 5～8\end{array}\right. $

(10)

格子声速为 $ 1/\sqrt{3} $ 。对式（4）进行时间和空间离散可得格子Boltzmann方程：

$ {f}_{i}\left(\boldsymbol{x}+{\boldsymbol{c}}_{i}\mathrm{\Delta }t,t+\mathrm{\Delta }t\right) = {f}_{i}\left(\boldsymbol{x},t\right)-\frac{\mathrm{\Delta }t}{\tau }[{f}_{i}\left(\boldsymbol{x},t\right)-{f}_{i}^{{\rm{eq}}}(\boldsymbol{x},t\left)\right] $

(11)

其中， $ {\boldsymbol{c}}_{i}\mathrm{\Delta }t $ 表明方程中 $ {f}_{i} $ 用一个时间步长 $ \mathrm{\Delta }t $ 以速度 $ {\boldsymbol{c}}_{i} $ 从一个格子点运动到对应方向的相邻格点。这表明LBM中的空间离散尺度和时间尺度是耦合的： $ {\boldsymbol{c}}_{i} = \mathrm{\Delta }{\boldsymbol{x}}_{i}/\mathrm{\Delta }t $ 。上述方程可以用以模拟弱可压流动，流场压力满足理想气体假设 $ p = \rho {c}_{s}^{2} $ ，因而下文中以密度场表征压力分布。

本文涉及点声源在流场中的叠加，采用Viggen^[36]提出的多极子声源方法。该方法可以统一描述多极子声源作用，并以源项 $ {s}_{i}(\boldsymbol{x},t) $ 形式被施加于式（11）：

$ \left[\begin{array}{c}{s}_{0}\\[16.4pt] {s}_{1}\\[15.5pt] {s}_{2}\\[16.2pt] {s}_{3}\\ {s}_{4}\\ {s}_{5}\\[19.2pt] {s}_{6}\\ {s}_{7}\\ {s}_{8}\end{array}\right] = \left[\left.\begin{array}{ccccccccc}{w}_{0}& 0& 0& -1& -1& 0& 0& -\dfrac{1}{2}& -\dfrac{1}{2}\\ {w}_{1}& \dfrac{1}{2}& 0& \dfrac{1}{2}& 0& 0& 0& 0& 0\\ {w}_{2}& 0& \dfrac{1}{2}& 0& \dfrac{1}{2}& 0& 0& 0& 0\\ {w}_{3}& -\dfrac{1}{2}& 0& \dfrac{1}{2}& 0& 0& 0& 0& 0\\ {w}_{4}& 0& -\dfrac{1}{2}& 0& \dfrac{1}{2}& 0& 0& 0& 0\\ {w}_{5}& 0& 0& 0& 0& \dfrac{1}{\sqrt{8}}& 0& \dfrac{1}{4}& 0\\ {w}_{6}& 0& 0& 0& 0& 0& \dfrac{1}{\sqrt{8}}& 0& \dfrac{1}{4}\\ {w}_{7}& 0& 0& 0& 0& -\dfrac{1}{\sqrt{8}}& 0& \dfrac{1}{4}& 0\\ {w}_{8}& 0& 0& 0& 0& 0& -\dfrac{1}{\sqrt{8}}& 0& \dfrac{1}{4}\end{array}\right]\right.\left[\begin{array}{c}{M}_{0}\\ {M}_{x}\\ {M}_{y}\\ {M}_{xx}\\ {M}_{yy}\\ {M}_{{x}'}\\ {M}_{{y}'}\\ {M}_{{x}'{x}'}\\ {M}_{{y}'{y}'}\end{array}\right] $

(12)

其中， $ {M}_{0}(\boldsymbol{x},t) $ 表示单极子作用， $ {M}_{\alpha }(\boldsymbol{x},t) $ 表示偶极子作用， $ {M}_{\alpha \beta }(\boldsymbol{x},t) $ 表示四极子作用；向量 $ \boldsymbol{M} $ 前的矩阵为转换矩阵。Viggen基于Chapmann-Enskog展开，证明了 $ {s}_{i} $ 的形式与多级声源作用相容。本文仅研究单极子声源，其形式为： ${s}_{i}\left({\boldsymbol{x}},t\right) = {w}_{i}{M}_{0} = {w}_{i}\delta \mathrm{sin}\left(\omega t\right)，\omega = 2\text{π} {c}_{s}/\lambda$ 为声源频率（ $ \lambda $ 表征声波波长），与前人所采用的单极子声源数学形式^[37]一致；同时向量 $ \boldsymbol{M} $ 的其他分量为0。

1.2 计算模型

点声源声波影响下的剪切层发展计算模型如图1所示。计算域上下两个区域分别有水平速度 $ {u}_{1} $ 和 $ {u}_{2} $ ，上下两层流体由剪切层分开（图1中橘色区域，中心处 $ y = 0 $ ），声源坐标为 $ (0,h) $ 。剪切层内速度剖面分布为：

$ \left\{\begin{array}{l}u(x,y) = \dfrac{{u}_{1}+{u}_{2}}{2}+\dfrac{{u}_{1}-{u}_{2}}{2}\mathrm{tanh}\left(\dfrac{2y}{{\theta }_{m}}\right)\\ v(x,y) = 0\end{array}\right. $

(13)

式中 $ {\theta }_{m} $ 指剪切层厚度。

计算域上下两边设置声波吸收条带，采用PML技术实现^[27]，这主要为了吸收计算域中的传出声波，不使边界反射波污染流场。在PML层中，额外增加起吸波效果的碰撞项：

$ {\mathrm{\varOmega }}^{{\rm{PML}}}\left({f}_{i}\right) = -\left[\epsilon{c}_{i\alpha }{\partial }_{\alpha }{Q}_{i}+2\epsilon{\tilde f}_{i}^{{\rm{eq}}}+{\epsilon}^{2}{Q}_{i}\right] $

(14)

$ \begin{split}\\{\partial }_{t}{Q}_{i} = {\tilde f}_{i}^{{\rm{eq}}} = {f}_{i}^{{\rm{{\rm{eq}}}}}\left(\rho ,\boldsymbol{u},t\right)-{f}_{i}^{{\rm{eq}}}\left(\bar{\rho },\bar{\boldsymbol{u}},t\right) \end{split}$

(15)

其中， $\; \bar{\rho }、\bar{\boldsymbol{u}}$ 分别为假想的远场流动参数， $ \epsilon $ 为吸波强度。为了有效抑制压力波在边界上的反射， $ \epsilon$ 以一定的数学形式给出，本文取：

$ \epsilon = {\epsilon}_{0}\frac{{l}^{2}}{{L}_{{\rm{PML}}}^{2}} $

$ l $ 为吸收层内网格节点到流场-吸收层边界的距离， ${L}_{{\rm{PML}}}$ 为吸收层厚度， $ {\epsilon}_{0} $ 取0.08。

在本文计算中，计算域尺度为 $ 5\;000\times 800 $ ，坐标系原点在计算域中心。声波吸收边界层厚度 $ 100 $ 。计算域上下两侧采用滑移边界条件，左右两侧为周期边界条件。在模拟中，计算域的长度足够长，使得在研究关注的时间范围内声波不会到达计算域两侧边界。声源密度波动幅值为背景流场密度 $\; {\rho }_{0} $ 的 $ 1\% $ ，满足线性声波要求。为了较为精确地解析声场作用下的剪切层演化及内部流动细节，本文在关注的受干扰剪切层区域上采用网格加密技术（图1中蓝色虚线标记区域）。网格逐层加密两次，每层网格分辨率增加一倍。本研究中，测试所用声波波长在8～64（格子长度，下同），剪切层厚度为 $ 3.2 $ ，因此本文网格可以解析流动细节。另一方面，由于在网格加密过程中，会不可避免地采用插值、数据归并等操作（本文采用文献[38]所用方法），所以会带来一定的误差。本文前期针对平直声波跨加密层边界、不同加密区域下剪切层失稳等算例开展了验证，结果表明当前网格加密技术所引起的额外的声耗散和色散效应较原始激励声波总体可以忽略，对剪切层失稳的影响可被忽略。

1.3 算例验证

计算结果表明，在初始无扰动情况下，剪切层可以在声波穿越过程中维持一段时间的稳定。在这段时间中，声传播方向满足Amiet理论的描述^[15]，即背景流场具有稳定的速度剖面时，多普勒效应导致声传播在剪切层上发生偏折。Amiet公式为：

$ \mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{c} = \frac{\chi }{M a+(1-{M a}^{2})\mathrm{c}\mathrm{o}\mathrm{s}\theta } $

(16)

$ \chi = \sqrt{{\left(1-M a\;\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)}^{2}-{\mathrm{cos}}^{2}\theta } $

(17)

其中， $M a = u/{c}_{s}$ 为马赫数， $ {\theta }_{c} $ 为入射声波的偏角， $ \theta $ 为经边界层折射的出射声波的偏角。作为验证性典型算例，本文取剪切层上侧马赫数为 $ 0.2 $ ，下侧为 $ 0 $ ；剪切层厚度 $ {\theta }_{m} = 0.8 $ 或 $ 0.4 $ ； $ h = 80 $ ；声波角频率为 $ 0.045 $ 。

图2给出某一时刻流场压力脉动（压力与密度成正比）分布云图，红线标记了剪切层位置，可以看到声场在剪切层两侧的偏折现象明显。进一步在图中取三个测点：（-40，-40）、（80，-40）、（120，-40），进行声波波前到达时间以及密度脉动幅度的测量。结果显示，LBM的计算结果与Amiet理论预测吻合较好（见表1），证实了本文数值方法的准确性。

2 结果与分析 2.1 平行剪切流失稳及声场作用的一般描述

本文研究不同波数声波与有一定厚度的水平剪切层相互作用后的流场和声场的变化。根据流向速度可将流场划分为3个部分：剪切层上侧（声源一侧）、剪切层区域以及剪切层下侧（声波接收区）。受速度梯度的影响，声波传播会发生反射与折射，剪切层上侧的声场同时受入射波与反射波影响，剪切层下侧基本只有折射声波。剪切层中不仅存在3种声波的叠加，还包含由声波激励起的K-H波造成的流场、压力场变化。

计算结果显示，在声传播过程中流场声场演化基本分为三个阶段：初始阶段，平直剪切层稳定，声波穿越剪切层过程主要具有折射特征，满足Amiet理论预测；第二阶段，随着流动的发展，剪切层对声传播的影响进一步显现，声压分布在不同的区域表现出一些特殊结构，其局部的流动参量分布开始偏离稳定流状态，此时剪切流场仍然没有可见的变化；第三阶段，在声场分布特殊结构的作用下，剪切层开始失稳，并进一步迅速发展为不稳定流动，其中剪切层的波动分布与前述声波的特殊结构相匹配。研究同时表明，声场对剪切层不稳定性的发展起到诱发作用，在剪切层不稳定性进一步发展的过程中则摆脱了声场的影响。

数值模拟发现，在演化过程中激励声波波数 $ k $ 会影响剪切层失稳状态。据前人关于剪切层流动稳定性的研究可知，剪切层失稳波的发展与初始扰动波波数之间存在一定关系。线性假设下剪切层稳定性分析方法较为成熟，典型的如模态分析法。其通常以背景流动满足平行假设为前提，将流动变量分解为背景流与扰动量： $ q\left(\boldsymbol{x},t\right) = \bar{q}\left(\boldsymbol{x}\right)+{q}'(\boldsymbol{x},t) $ 。基于正则模态分析将小扰动量 ${q}'(\boldsymbol{x},t)$ 作Fourier展开为 $\widehat{q}{{\rm{e}}}^{{\rm{i}}\left(kx-\omega t\right)}$ ，其中 $ k $ 为波数， $ \omega $ 为频率。因本文算例中剪切层附近流动速度为0，且关注剪切层失稳波的形成过程，故采用时间模式的稳定性分析： $ k $ 为实数， $ \omega $ 为复数。 $ \omega $ 的实部和虚部分别为K-H波的频率和扰动增长率。将波形式的扰动带入线性化Navier-Stokes方程，可得稳定性方程： $ A\widehat{q} = B\omega \widehat{q} $ ， $ A $ 和 $ B $ 代表相应的系数矩阵，其与流动速度分布、波长等参数相关。本文采用文献[39]的程序计算，基于Chebyshev配置点谱方法，获得不同波数条件下扰动波的增长率（见图3），其中实心点表示本文分析的算例，空心点表示测试过的算例。流场演化满足正文所提的长波、短波发展规律。可以看到当前剪切流参数条件下，不稳定波数 $ k $ 范围为（0.07，0.61），最大增长率的不稳定波波数为 $ 0.3 $ 。

测算结果表明，以 $ 0.61 $ 为界，短波长和长波长声波对于剪切层失稳的影响有明显区别。本文围绕两个算例：声波波数 $ k = 1.59 $ （算例1）和 $ 0.265 $ （算例2）进行分析。上下两层流体中，流动速度为 $ {u}_{1} = -{u}_{2} = 0.0577 $ ，对应 $ M = \pm 0.1 $ 、 $ {\theta }_{m} = 3.2 $ ，声源位置（0，80）。

2.2 短波长声波作用下的剪切层失稳

本节主要针对短波长声波传播与剪切层发展耦合现象（算例1）开展分析。图4给出了 $ t = 184.7 $ 时， $ -400\leqslant x\leqslant 400 $ 、 $ -40\leqslant y\leqslant 40 $ 区域的密度脉动场，图中的密度不均匀分布包含了声波影响（细小条纹结构）与流动影响（较大尺度的密度变化）。其中 $ t $ 为无量纲时间（下同），特征时间为 $ 2L/\left|{u}_{1}-{u}_{2}\right| $ 。图中显现了在声场作用下剪切层失稳K-H波的典型特征：压力波峰和波谷交替出现。基于剪切流动的稳定性，整个流场分为三个区域。位于声源左侧的区域1中剪切流动最早失稳，且肉眼可见的K-H波动首先出现在该区域的右侧，即区域1与区域2的衔接部附近。模拟结果同时显示，在上下两层流体中，声波的传播相对独立。区域2位于声源的下方（声波垂直入射点两侧），声波穿透剪切层的现象明显，声传播相对连续。这个区域的流动失稳较区域1滞后。区域3处于声源的右侧（流动下游方向），该区域中流动可以长期维持稳定。穿越稳定剪切层的声传播过程可以用前文所提的Amiet公式加以描述。

上述描述表明，区域1中的流动具有较为丰富的物理内涵，而区域2中的声场-流动演化与2.3节所述内容相近。因此，下面的内容主要针对区域1的声场分布与不稳定性演化过程加以展开。

图5给出了 $ -400\leqslant x\leqslant 0 $ 、 $ -15\leqslant y\leqslant 15 $ 区域压力脉动场随时间的演化云图。在上层流体中可以发现，声场压力显现出明显的放射状的条带结构（图5（a）中红色虚线指示了自由传播声波与剪切层反射声波的叠加）。考虑到剪切层上下两侧流体的声传播特性有差异，剪切层对声波具有折射和反射作用。自由传播的声波与反射声波叠加，在上层流体区域构成干涉条带结构。更进一步，当入射声波与剪切层平面夹角超过一个临界值时，声波会在剪切层上发生全反射。尽管本文考虑的剪切层厚度有限，但该夹角仍可用Amiet理论进行估算：

$ \theta \left({M}_{1}\right) = \mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left[-\frac{1}{1+{M}_{1}}\right] $

当前算例中， $\theta = {146.44}^{\circ }$ ，对应 $x = -120.6$ （图5中竖直黄色虚线表示声波全反射区起始位置）。在该标记线左侧的下层流体中，没有穿透剪切层的声波，而只有上游透射声波的影响，因而声压条纹较为规整。在更为重要的剪切层区域，可以观察到压力脉动的波包结构。该波包结构一开始在区域1的右侧相对孤立，间距较大；随着时间的推进，该波包结构逐渐变强，并最终导致了流动的失稳。

为了进一步看清压力波包结构的发展，我们提取了剪切层附近一条水平线 $ y = -2 $ 上压力的演化过程（图6）。图6（a）显示，在剪切流较为稳定的阶段，压力波包结构（ $ x < -160 $ ）较为清晰、规整。在声波全反射点附近，波包结构平滑过渡到区域2。在区域2中由于声波入射角不同，导致提取的压力脉动波长向声源方向增加。区域1中的压力波包结构暗示了两道波长相近的声波在此区域叠加。针对图6（a）中 $ x\leqslant -200 $ 区域的压力分布进行频谱分析（图7），可以得到这两道波的波数分别为 $ {k}_{1} = 1.44 $ 和 $ {k}_{2} = 1.69 $ 。考虑背景流动带来的声传播多普勒效应，剪切层上下两部分的有效声传播速度分别为 $ {c}_{e} = 0.9{c}_{s} $ 和 $ 1.1{c}_{s} $ ，且声波频率 $ \omega = 0.181 $ ，对应的波数 $ k = \omega /{c}_{e} $ 分别为 $ 1.43 $ 和 $ 1.74 $ 。数值模拟结果与理论结果吻合。图7同时给出了区域1中不同时刻的压力波空间分布频谱，其显示声波分量在K-H波发展过程中维持稳定。

图6进一步显示，在时刻 $ t = 126.93 $ 以后，前述波包结构进一步发展为剪切层的K-H波。对比 $ t = 69.24 $ 和 $ t = 126.93 $ 时密度波的包络线（以彩色曲线标记），可以发现声波与流体压力波动的叠加，导致了包络线的畸变。剪切层不稳定性最终发展为如 $ t = 150.01 $ 时刻大尺度的压力波动，这也可以从图5中被看到。在频谱图7中，剪切层不稳定性波的波数为 $ 0.245 $ 。该数值与前述两个声波叠加结果吻合： ${k}_{_{{\rm{KH}}}} = {k}_{2}-{k}_{1} = 1.69-1.44 = 0.25$ 。由流动不稳定性的发展可知，该波数的波是随着时间增长的。图6中点划线所标记的密度过零点位置以及图7的压力脉动分量演化均显示，虽然波包结构是导致K-H不稳定波产生的直接原因，但在流动失稳被诱导出以后，其发展具有一定的独立性。这是因为失稳波扰动流动动能与声波能量在量级上存在差异。

最后需要指出两点：

1） 结合区域3的声场压力分布（图4）以及声波在剪切层全反射的条件判据可知，区域3中，声波以穿越剪切层为主，因而不存在区域1中的波叠加现象，也就无法诱导出类似于区域1中的波包结构；并且声波脉动本身的空间尺度又明显小于K-H最不稳定波的尺度。因而此区域的K-H波发展缓慢，流动相对稳定。

2） 尽管区域1剪切层中波包结构较为规整，但从图6可以看出靠近声源区域的波包波长（点划线指示部分）仍然存在变化，这主要是由声波、背景流动速度、剪切层的夹角在这个区域变化较大所造成的。

2.3 长波长声波作用下的剪切层失稳

图8给出了算例2中， $ -400\leqslant x\leqslant 200 $ 、 $ -15\leqslant y\leqslant 15 $ 区域压力脉动场随时间的演化云图，红线标记了声源水平位置。算例中采用的声波波数为0.265，根据2.1节的内容，该波长理论上处于剪切层失稳区域。声源左侧区域中剪切层上下侧声场分布仍然是不同的，但剪切层区域内未发现明显的如2.2节所述的密度波包的情形。其原因可能有三种：第一种，该算例中区域2的不稳定波发展较快，对压力场造成的影响较大，从而弱化了波包结构的显示；第二种，由于当前声波波长较大，由不同波数引起的波包结构波长较大，受计算域大小的限制，波包结构不能被观察出来；第三种，当前情况下，区域1波包结构对应的波数已经处于振幅负增长区域，这意味着流动是稳定的。

在当前算例中，区域2的剪切流动首先失稳，其表现为剪切层中，剪切层上声波垂直入射位置的两侧出现压力极小值区（图8（b）中的亮斑）。该压力极值区将进一步发展出剪切层变形，并最终失稳（图8（e））。对于其他大波长声波（ $\lambda = 16、 48$ ），其失稳过程与当前描述相符。从图8中可以看出区域1的失稳与区域2交界的地方没有出现明显的不稳定波，因而两者的流动失稳没有直接关联。观测图8（b～d）中剪切层波动几何特征，可以得到本算例中发展出的不稳定波波长约为42，对应波数0.15。图3表明其并不对应剪切层最不稳定波的波数。图9给出了算例2中 $ -400\leqslant x\leqslant 200 $ 、 $ -15\leqslant y\leqslant 15 $ 区域失稳后的涡量分布，其中红线依旧指示声源水平位置。

与2.2节中区域1的声波入射条件不同，区域2中声波入射角度受水平位置影响更大。二维平面内点声源的传播有理论解，其分布可以用Bessel函数加以描述（本文采用第一类0阶Bessel函数： $ {J}_{0}\left(r\right) $ ）。在此给出静止、均匀流体中点声源造成的压力波动沿剪切层直线的分布（图10（a）中红色虚线表示均匀流中密度波动分布，蓝色点线基于理论解 $ {J}_{0}\left(r\right) $ 获得），其与数值结果具有一定的相似性。但同时也发现，算例结果中剪切层上的密度波相对于垂直入射点不对称分布，并且随着时间的推移，垂直入射点两侧附近压力波动波长减小、幅值逐渐增加。另一点值得注意的是，区域1中的声波干涉条纹也是逐渐出现的，并且条纹结构总体指向区域2的不稳定区域（图8（a）中的虚线）。

上述现象表明，含剪切层和背景流空间中，声波的传播发展具有特殊性。为了区分两者的作用，计算了声波在均匀流场中（即 $ {u}_{2} = {u}_{1} $ ）的声场分布，并测量了与剪切层相对应的直线上的压力分布（图10（a））。可以看出：

1） 均匀流条件下声波波幅在垂直入射点两侧呈衰减趋势，且衰减速率与经典声学理论解（Bessel函数）预测一致；

2） 剪切层条件下，声场在垂直入射点右侧与均匀流一致，但左侧声脉动幅值明显提高，且在远场方向上的衰减微弱（或不衰减）；

3） 剪切层条件下，声波垂直入射点左侧剪切层中压力波动波长发生了变化。

上述现象是在剪切层没有明显失稳的前提下被观察到的，由此我们可以判断：尽管在区域2中声波的演化、传播仍然具有行波特征，但剪切层的反射作用加剧了区域2左半部分的声压脉动幅值；并且由于非线性效应，压力波动的波长逐步减小，幅值随时间增加。剪切层对声波的透射和反射作用具有指向的选择性，这进一步使得剪切层上压力波动分布左右不对称，加剧了流动的不均匀性。综合2.2节所得结果，我们可以得到声场波动幅值的空间不均匀性，促使剪切层流动的不稳定发展。

本节结果表明区域2中的声波传播受更多因素的影响，比如剪切层对不同角度入射声波的折射/反射规律、声-流耦合；而现有理论模型、数值和实验结果仍然不足以支撑对有关现象的精确描述，因此有待将来开展更深入的物理机理研究。

3 结　论

本文采用格子Boltzmann数值方法，对声波传播与平行剪切流稳定性发展的声作用机理开展了研究。探讨了声波在剪切流界面上的折射、反射现象，发现剪切层对声波的反射具有方向选择性。发现多普勒效应—声波在上下两层流体中的波长变化是导致流动稳定性的一个重要诱因。最终结合平行剪切流动稳定性分析，探讨了声压幅值分布对流动稳定性的影响；发现声场对剪切层稳定影响可以分为（声波）长波模式和短波模式。基于声源一侧流动方向，剪切层不稳定性发展可分为上游区、声波近垂直入射区和下游区。三个区域的失稳原因和失稳条件存在不同。研究结果表明：

1）在背景流动下游，声传播方向与流动方向夹角较小，声波透过剪切层成分较高，其传播规律服从经典Amiet理论，声场作用对剪切层稳定性影响较小。在背景流动上游区域，随着声波传播方向与来流方向夹角的增长，剪切层上的声反射作用增强，直至全反射出现，反射波与入射波在声源一侧区域形成干涉条纹结构。

2）在短波长声波情况下，剪切层上侧声波与下侧声波形成干涉。由于其波长不同，叠加声波构成稳定的波包结构，该结构的波长如果与剪切层失稳波长吻合，会诱导产生流动失稳。在流动失稳形成后，流动与声波发展可以解耦。

3）在长波长声波条件下，前述剪切层上游的干涉波包结构较弱，声源垂直入射区扰动增长占优。剪切层的反射作用使得该区域上游部分声压脉动增强，造成剪切层内声扰动分布对称性的缺失，同时压力波动波长随时间减小；当这样的波长落入失稳波长区域后，K-H波被诱导产生，并逐渐进入相对独立的发展阶段。

本研究揭示了不同波长声波作用下，平行剪切流动失稳的机理，可以为流动稳定性控制、声探测、流动噪声调制等研究提供借鉴。研究同时也显示，在声场作用下，剪切层内部流动、剪切层的声反射、声传播机制等方面的研究工作仍然亟待开展。