Computational methods for acoustic propagation based on linearized equations
-
摘要:
复杂流场如剪切层、旋涡等会改变气动噪声的传播特性，引起折射、反射和散射等现象，对声源识别、测量产生影响。基于线化方程的声传播计算方法是研究声波在复杂流场中传播的重要手段。本文针对非均匀流中声传播计算存在的问题，介绍了近年来课题组开展的研究工作：提出了改进梯度项抑制方法，以抑制基于线化方程模拟声波穿过剪切层时出现的数值不稳定波；发展了适用于线化方程的基于Boltzmann模型的有限体积法通量计算格式，以模拟声波在包含复杂外形流场中的传播问题；发展了简化的线化格子Boltzmann方法，改善了格子Boltzmann方法模拟声传播时内存占用大的问题；研究了剪切层对声源定位的影响规律，建立了定位误差随射流马赫数、斯特劳哈尔数变化的数学模型，为提高实验测量精度提供参考。
-
关键词:
- 线化方程 /
- 声传播 /
- 格子Boltzmann模型 /
- 计算气动声学
Abstract:Complex flows such as shear layers and vortices can change the propagation characteristics of aerodynamic noise, resulting in phenomena of refraction, reflection, and scattering, which have an impact on acoustic source measurement and identification. The linearized equation-based algorithm is an important tool for simulating the propagation of acoustic waves in complex flows. This paper briefly outlines the research carried out by our group in recent years to address the challenges associated with acoustic propagation in non-uniform flow. Specifically, we propose an improved gradient term suppression method to mitigate numerical instability when calculating acoustic wave propagating through shear layers using the linearized equations; develop a Boltzmann model-based flux solver under the finite volume method for the linearized equations to simulate acoustic wave propagation in flow fields containing complex geometries; establish a simplified linearized lattice Boltzmann method to reduce the computational memory requirements when using the lattice Boltzmann method; and evaluate the influences of shear layer characteristics on sound source localization, and establish a mathematical model of localization error with jet Mach number and Strouhal number to provide reference for improving the experimental measurement accuracy.
-
0 引 言
在实际工程问题中，大量存在着声波在旋涡、剪切层、激波等复杂流场中传播的现象[1-6]。例如在声学风洞中，试验件产生的噪声需穿过风洞出口射流与周围静止空气形成的剪切层，才能到达远场的麦克风阵列。研究表明：声波穿过旋涡、剪切层时会发生折射、反射、散射等现象，使麦克风接收信号产生幅值、频率和相位变化，进而影响测量及声源定位的精度[2-6]。航空发动机噪声中也有类似情形：涵道旁路冷流与喷口热射流、自由来流之间形成热、冷剪切层，极大改变了风扇和涡轮噪声向下游的传播特性。因此，开展声波在复杂流场中的传播研究具有重要意义。
数值模拟是研究声波在流场中传播的重要手段。随着计算流体力学（Computational Fluid Dynamics, CFD）的发展，学者们引入了CFD技术以研究气动声学问题。声波的本质是小扰动，而声传播模拟就是研究小扰动在流场中的传播过程。目前采用CFD模拟声波在流场中传播的方法主要有2种。其一，采用直接数值模拟对全流场进行计算，得到包含声扰动在内的流场物理量。由于声扰动通常比主流小数个量级，采用该方法模拟时需要很密的计算网格及很小的时间步长，使得声传播模拟的计算量巨大。其二，采用求解线化方程如线化Euler方程[7-10]、线化Navier‒Stokes（N−S）方程[11]的方法。该方法由于直接对扰动量进行计算，可以避免主流物理量将声扰动掩盖，在声传播模拟中得到广泛应用。但该方法仍然面临一些问题，具体表现在：
1）采用线化Euler方程模拟声波穿过梯度较强的剪切层时，方程中的速度扰动项与梯度项相互作用，容易产生Kelvin‒Helmholtz不稳定性，生成不稳定涡波，影响计算精度甚至导致计算发散。时域的不稳定波抑制方法主要分为梯度项抑制[12]和梯度项滤波[13]。梯度项抑制的概念最早由Bogey等[12]提出，他们通过直接移除动量方程和能量方程中的梯度项，有效减少了不稳定波。Ewert等[14-16]借助等熵假设将连续方程并入能量方程，推导出了声扰动方程。Seo等[17]参考移除梯度项的做法，削弱了涡波干扰。Hu[18]和李旦望[19]等采用在动量方程右端添加源项的方式来修正梯度项，并引入了波包概念，实现了多个频率的同时计算和滤波。Prax等[20]发现梯度项抑制方法并非完全无旋，需引入右端源项来削弱不稳定波。研究表明，梯度项抑制方法可减少不稳定波对计算的影响，但由于梯度项对声波的折射与反射有贡献，移除后计算误差会随着梯度的增强而增大[21]。对此，Deng等[22]将线化方程中的梯度项乘以系数，使其既能削弱不稳定波，又考虑了背景流场的影响，但系数需人为调节。Zhang等[13]提出了梯度项滤波的概念，认为声波对应Helmholtz分解出的无旋速度扰动，仅保留线化方程中的无旋速度扰动，能在一定程度上削弱不稳定波的影响。由此可见，声波穿过剪切层传播模拟引起的数值不稳定波问题仍未被完全解决，是一个亟需研究的重要问题。
2）求解线化Euler方程或线化N−S方程时，常采用有限差分法在结构网格上进行离散。有限体积法由于能应用于非结构网格，且对复杂外形具有较好的适应性，在流场求解中得到广泛应用。为了更好地模拟声波绕复杂外形的传播过程，采用有限体积法离散线化Euler方程或线化N−S方程是一种较优选择。通量计算是有限体积法的重要过程。传统流体力学控制方程的无黏通量计算，将单元交界面视作黎曼间断面，基于一维黎曼模型计算交界面的法向通量，采用被动标量法近似计算切向通量；而黏性通量计算则采用不同于无黏通量计算的方法进行，增加了算法的复杂性。基于Boltzmann模型的通量计算格式如气体动理学格式[23-25]和基于格子、连续Boltzmann模型的通量求解器[26-28]等可以推广至多维模型，无需分开计算无黏通量和黏性通量，算法一致性更好，近年来已用于各种流动问题的求解[29-38]。但这些格式并不适用于线化方程的通量计算。因此，如何发展基于Boltzmann模型的线化方程通量计算格式是另一个亟需解决的重要问题。
3）格子Boltzmann方法（Lattice Boltzmann Method, LBM）具有无需求解微分方程、在相同精度下数值耗散更低[39]、容易并行[40]等优点，近年来在气动声学计算中得到应用[41-43]。为了更好地模拟声扰动在流场中的传播过程，Vergnault等[44]将粒子分布函数分解为主流部分和扰动部分，推导出了线化格子Boltzmann方程。Pérez等[45]提出了碰撞算子与局部概率分布函数的线性化方法，推导了扰动粒子分布函数与宏观量之间的矩关系，证明了线化格子Boltzmann方程可以还原为线化N‒S方程以及利用迁移碰撞可以实现扰动粒子分布函数的演化。然而无论是LBM还是线化LBM，计算中都需存储网格点不同粒子速度方向上的分布函数，内存占用较大[46]；同时，由于LBM演化的是粒子分布函数，边界条件需转化为粒子分布函数才能施加[44, 47]。因此，线化LBM模拟中存在的内存占用大和边界处理难的问题也是声传播计算需要解决的重要问题。
下面针对以上基于线化方程的声传播计算方法存在的问题，介绍近年来课题组开展的研究工作。
1 线化Euler方程不稳定波抑制研究
1.1 线化Euler方程的数值不稳定波问题
线化Euler方程为：
$$ \dfrac{{ \partial \rho '}}{{ \partial t}} + \nabla \cdot \left( {\rho '{\bar {\boldsymbol{u}} }+ \bar \rho {{{\boldsymbol{u'}}}}} \right)= \dfrac{{ \partial {{\rho '}_{\rm{ S}}}}}{{ \partial t}} $$ (1) $$ \dfrac{{ \partial (\bar \rho {\boldsymbol{u'}})}}{{ \partial t}} + \nabla \cdot \left( {\bar \rho {\boldsymbol{u'}}\bar {\boldsymbol{u}}} \right) + \nabla p' + \boxed{\bar \rho {{{\boldsymbol{u'}}}_{}} \cdot \nabla {\bar {\boldsymbol{ u}}}} + \rho '{\bar {\boldsymbol{ u}}} \cdot \nabla {\bar {\boldsymbol{ u}}} = {\mathbf{0}} $$ (2) $$ \dfrac{{ \partial p {'}}}{{ \partial t}} + \nabla \cdot \left( { p {'} {\bar {\boldsymbol{ u}}}}+{\text{γ}}\,\bar p {\boldsymbol{ u{'}}} \right) + ({\text{γ}}-1)[p{'}\nabla \cdot {\bar {\boldsymbol{ u}}}- \boxed{{\boldsymbol{ u{'}}} \cdot \nabla \bar p }\,]= \dfrac{{ \partial p{ {'}}_{{\mathrm{ S}}}}}{{ \partial t}} $$ (3) 式中：$ \gamma $、$t$、$\rho$、$p$和$ \boldsymbol{u}\text{ = }\left(u,v,w\right) $分别为比热比、时间、密度、压力和速度；上标“”表示扰动量，上标“$\bar \;$”表示平均量，下标“s”表示源项；方框“$ \boxed{\mathop { }\limits^{} } $”标记了与不稳定波相关的项。$”表示扰动量，上标“$\bar \;$”表示平均量，下标“s”表示源项；方框“$ \boxed{\mathop { }\limits^{} } $”标记了与不稳定波相关的项。
方框内的速度扰动${\boldsymbol{u'}}$遇到强梯度项$ \nabla {\bar {\boldsymbol{ u}}} $与$ \nabla \bar p $时，容易引发数值不稳定波。Bailly等[7]将流场简化为平行剪切层，用稳定性方法分析了不稳定波问题。Zhang等[13]对剪切层失稳产生涡波的特征频率进行了理论预测。虽然稳定性分析方法难以用于实际剪切层，但速度扰动与速度/压力梯度的相互作用已被认为是产生不稳定波的原因。
1.2 改进梯度项抑制方法
直接移除梯度项是处理数值不稳定波最简单的方法，但随着剪切层速度梯度的增大，梯度引起的声折射效应不可忽略。为了去除不稳定波，同时考虑剪切层流场梯度对声折射的影响，课题组提出了改进梯度项抑制方法。
改进梯度项抑制方法受Zhang等[13]的梯度项滤波方法启发，共分为2步：第一步，用传统梯度项抑制方法获得近似无旋的速度扰动；第二步，用该速度扰动过滤完整线化Euler方程的梯度项。每个时间步分别采用有、无梯度项抑制的线化Euler方程求解，并用梯度项抑制方法得到的速度扰动单向替换线化Euler方程中梯度项的速度扰动，以达到滤波目的。
借助Seo等[17]的量纲分析方法，可说明高马赫数Ma下使用改进梯度项抑制方法的必要性。传统梯度项抑制方法的控制方程为：
$$ \dfrac{{ \partial \rho '}}{{ \partial t}} + \nabla \cdot \left( {\rho '{\bar {\boldsymbol{ u}}} + \bar \rho {{\boldsymbol{u}_{\mathrm{f}}^{\prime}}}} \right){\text{ }} = \dfrac{{ \partial {{\rho '}_{ {\mathrm{S}}}}}}{{ \partial t}} $$ (4) $$ \dfrac{{ \partial \bar \rho {\boldsymbol{u}_{\mathrm{f}}^{\prime}}}}{{ \partial t}} + \nabla \cdot \left( {\bar \rho {\boldsymbol{u}_{\mathrm{f}}^{\prime}} {\bar {\boldsymbol{ u}}}} \right) + \nabla p' = {{0}} $$ (5) $$ \dfrac{{ \partial p'}}{{ \partial t}} + \nabla \cdot \left( {p'{\bar {\boldsymbol{ u}}} + \gamma \bar p{\boldsymbol{u}_{\mathrm{f}}^{\prime}}} \right) = \dfrac{{ \partial {{p'}_{ {\rm{S}}}}}}{{ \partial t}} $$ (6) 引入滤波后的线化Euler方程为：
$$ \dfrac{{ \partial \rho '}}{{ \partial t}} + \nabla \cdot \left( {\rho '{\bar {\boldsymbol{ u}}} + \bar \rho {\boldsymbol{u'}}} \right){\text{ }} = \dfrac{{ \partial {{\rho '}_{ {\rm{S}}}}}}{{ \partial t}} $$ (7) $$ \dfrac{{ \partial (\bar \rho {\boldsymbol{u'}})}}{{ \partial t}} + \nabla \cdot \left( {\bar \rho {\boldsymbol{u'}} {\bar {\boldsymbol{ u}}}} \right) + \nabla p' + \boxed{\bar \rho {{\boldsymbol{u}_{\mathrm{f}}^{\prime}}} \cdot \nabla {\bar {\boldsymbol{ u}}}} + \rho '{\bar {\boldsymbol{ u}}} \cdot \nabla {\bar {\boldsymbol{ u}}} = {{0}} $$ (8) $$ \dfrac{ \partial p^{\prime}}{ \partial t} + \nabla \cdot \left(p^{\prime} \overline{\boldsymbol{u}} + {\gamma} \,\bar{p} \,\boldsymbol{u}^{\prime}\right) + ({\gamma}-1)\left(p^{\prime} \nabla \cdot \overline{\boldsymbol{u}} -\boxed{{\boldsymbol{u}_{\mathrm{f}}^{\prime}}\cdot\nabla \bar{p}}\,\right)=\dfrac{ \partial p_{{\mathrm{ S}}}^{\prime}}{ \partial t} $$ (9) 式中：${\boldsymbol{u}_{\mathrm{f}}^{\prime}}$为用于滤波或需要滤波的速度扰动。控制方程中仅方框内的${\boldsymbol{u}_{\mathrm{f}}^{\prime}}$需要滤波。
假设背景流场为平行剪切层，可知$ \nabla\cdot\overline{\boldsymbol{u}}=0 $。进一步用$\bar \rho \sim{\rho _{{\infty}} }$，${\bar u}\sim{u_{{\infty}} }$，$\bar p\sim{\rho _{{\infty}} }u_{{\infty}} ^2$，$t\sim l/{c_{{\infty}} }$，$ u'_{ }\sim M a\cdot u_{ }^{(\mathbf{1})} $，$(p',{p'_{{\mathrm{ S}}}})\sim{\rho _{{\infty}} }{c_{{\infty}} }u'$，$(\rho ',{\rho '_{{\mathrm{ S}}}})\sim({\rho _{{\infty}} }/{c_{{\infty}} }){u'}$，$ u_{ }= \overline{u}_{ }+M a\cdot u_{ }^{(\mathbf{1})}+\cdots $等关系进行量纲分析，（其中下标${{\infty}} $表示远场参数，l为特征长度，c为声速）可知改进梯度项抑制方法的马赫数相关性为：
$$ \begin{gathered} \underbrace {\dfrac{{{ \partial ^2}p'}}{{ \partial {t^2}}}}_{O\left( M a \right)} + \underbrace {{\bar {\boldsymbol{ u}}} \cdot \nabla \dfrac{{ \partial p'}}{{ \partial t}}}_{O\left( {{M a^2}} \right)} - \dfrac{{ \gamma \bar p}}{{\bar \rho }}\nabla \\ \cdot\left[ {\underbrace {\nabla \cdot \left( {\bar \rho {\boldsymbol{u'}} {\bar {\boldsymbol{ u}}}} \right) + \nabla p'}_{O\left( {{M a^3}} \right)} + \overbrace {\underbrace {\bar \rho \left( {{{{\boldsymbol{u'}}}_{\mathbf{ f}}} \cdot \nabla } \right){\bar {\boldsymbol{ u}}}}_{O\left( {{M a^4}} \right)} + \underbrace {\rho '\left( {{\bar {\boldsymbol{ u}}} \cdot \nabla } \right){\bar {\boldsymbol{ u}}}}_{O\left( {{M a^5}} \right)}}^{{\text{gradient terms of momentum equation}}}} \right] \\ + \overbrace {\dfrac{{\left( { \gamma - 1} \right)}}{{\bar \rho }}\left[ {\underbrace {\nabla \cdot \left( {{\bar\rho }{\boldsymbol{u}_{\mathrm{f}}^{\prime}} {\bar {\boldsymbol{ u}}}} \right)}_{O\left( {{M a^4}} \right)} + \underbrace {\nabla {{p'}}}_{O\left( {{M a^3}} \right)}} \right] \cdot \nabla \bar p}^{{\text{gradient terms of energy equation}}} = \underbrace {\dfrac{{{ \partial ^2}{{p'}_{ {\rm{S}}}}}}{{ \partial {t^2}}}}_{O\left( M a \right)} \\ \end{gathered} $$ (10) 式中：各项皆为${\rho _\infty }c_\infty ^3u^{({\mathbf{1}})}/{l^2} \cdot O\left( {{M a^n}} \right)$，其中上标“n”为阶数。动量方程与能量方程的梯度项是Ma的3次方以上，而等式左端第一项和右端源项都是Ma的1次方。低马赫数的指数次方量级较小，因此Seo等[17]将梯度项移除；Ma增大后，梯度项的折射效应不可忽略，故在改进梯度项抑制方法中保留了这些梯度项。
线化Euler方程采用有限差分法求解，空间离散使用7点4阶的色散关系保持格式，并借助11点10阶滤波器消除高频振荡。时间离散采用低耗散低色散的优化5步龙格库塔显式时间推进格式。
1.3 方法验证
采用声波在高速高温剪切层中传播的经典算例进行验证[18]。该算例以平行剪切层为背景流场，内侧为高温、高速气流，外侧为静止大气。当采用线化Euler方程进行计算时，较强的温度和速度梯度会导致下游声场出现明显的不稳定波。算例的平均流场如下：
$$ \left\{\begin{gathered} \dfrac{\overline{\rho}_{\text{jet}}}{\overline{\rho}} = \dfrac{\overline{T}_{\text{air}}}{\overline{T}_{\text{jet}}} - (\dfrac{\overline{T}_{\text{air}}}{\overline{T}_{\text{jet}}} - 1 )\dfrac{\overline{u}}{\overline{u}_{\text{jet}}} + \dfrac{\gamma - 1}{2}M a_{\text{jet}}^2\dfrac{\overline{u}}{\overline{u}_{\text{jet}}}( 1 - \dfrac{\overline{u}}{\overline{u}_{\text{jet}}} ) \\ \overline{u}=\overline{u}_{\text{jet}}\exp\left[-\mathrm{ln}2(\textit{y}/l_{\mathrm{ref}})^2\right] \\ \overline{v}=0 \\ \overline{p}=103\, 330\; \text{Pa} \\ \end{gathered}\right. $$ (11) 式中：下标“jet”和“air”分别表示射流和远场空气参数；$ M a = \overline{u}_{\text{jet}}/c_{\text{jet}} = 0.756 $，${\bar \rho _{\text{jet}}} = 0.6\;{\mathrm{kg}}/{\mathrm{m}}^3$，${\bar T_{\text{jet}}} = 600\;K$，${\bar T_{\text{air}}} $ = 300 K，${c_{\text{jet}}} = 491$ m/s，${l_{{\mathrm{ref}}}} = 1.3$ m。源项为：
$$ \left\{\begin{gathered}\dfrac{\partial\rho'_{ \rm{S}}}{\partial t}=0 \\ \dfrac{\partial p'_{ \rm{S}}}{\partial t}=A'\exp\left(-\vartheta_ax_{ }^2-\vartheta_b{\textit{y}} ^2\right)\cos(\omega t) \\ \end{gathered}\right. $$ (12) 式中：x、${\textit{y}} $为横、纵坐标；$ A'=0.001\; \rm{kg}\cdot\mathrm{m}^{-1}\cdot\mathrm{s}^{-3} $；${\vartheta _a}$ = 0.04ln2，${\vartheta _b} $ = 0.32ln2；圆频率$\omega $ = 76 rad/s。
图1为利用不同方法模拟声波穿过高速高温剪切层得到的瞬时脉动压力场。由图1（a）中线化Euler方程的结果可知：声源下游存在沿流向演化的不稳定波涡波，这些涡波与边界的反射波叠加后，会逐渐“污染”整个声场。为削弱涡波的影响，传统梯度项抑制方法直接移除梯度项，得到如图1（b）所示的声场。虽然图1（b）不存在涡波，但是绿色区域的范围与图1（a）存在区别。在图1（c）的改进梯度项抑制方法所得结果中，既不存在涡波，又获得了与图1（a）接近的声场特征。
![]()
提取第15个声源周期的瞬时脉动压力分布，并将传统梯度项抑制方法（GTS‒1）和改进梯度项抑制方法（GTSF）的计算结果与解析解进行对比，如图2所示。由图可知，传统梯度项抑制方法在声源上下游区域受剪切层折射影响而出现较大误差，而改进梯度项抑制方法的计算结果与解析解吻合更好。
![]()
2 线化方程的通量计算格式研究
2.1 基于Boltzmann模型的线化方程通量计算格式
针对采用有限体积法求解线化Euler方程或线化N−S方程的通量计算问题，提出了基于Boltzmann模型的线化方程通量计算格式[50-51]。该格式是基于线化格子Boltzmann方程推导得到，离散后的线化格子Boltzmann方程可以写为：
$$ \begin{aligned} & {f'}_{ \alpha}\left(\boldsymbol{X}+{{ξ}}_{\alpha}\delta_t,\; t+\delta_t\right)-{f'}_{ \alpha}(\boldsymbol{X},t) \\ & =-\dfrac{1}{\tau_0}\left[{f'}_{ \alpha}(\boldsymbol{X},t)-{f'}_{ \alpha}^{\rm{eq}}(\boldsymbol{X},t)\right]\end{aligned} $$ (13) 式中：$f'$为粒子分布函数的扰动分量；二维粒子速度模型采用的是经典D2Q9模型[52]，${{{ξ}}_\alpha }$为第α个粒子的速度；$\delta_t $为粒子迁移的时间步长；$\tau $0为无量纲松弛时间；${\boldsymbol{X}} $为空间坐标(x，${\textit{y}} $)。线化格子Boltzmann方程对应的平衡分布函数可写为：
$$ \begin{aligned} & {f'}_{ \alpha}^{\rm{eq}}=\dfrac{\rho'}{\overline{\rho}}\overline{f}_{\alpha}^{\rm{eq}} \\ & +\overline{\rho}w_{\alpha}\left[\dfrac{{ξ}_{\alpha}\cdot\boldsymbol{u'}}{c_{\rm{S}}^2}+\dfrac{\left({ξ}_{\alpha}\cdot\boldsymbol{u'}\right)\left({ξ}_{\alpha}\cdot\overline{\boldsymbol{u}}\right)}{c_{\rm{S}}^4}-\dfrac{\boldsymbol{u'}\cdot\overline{\boldsymbol{u}}}{c_{\rm{S}}^2}\right]\end{aligned} $$ (14) 式中：上标“eq”表示粒子分布函数的平衡部分；${w_\alpha }$为第α个粒子的权重；cS为无量纲声速，其取值可参考文献[52]。通过建立线化Boltzmann方程与宏观线化方程间的联系，可利用介观变量计算线化方程的通量。对平衡分布函数分别求零阶、一阶和二阶矩，对非平衡分布函数求二阶矩，得到下列矩关系式：
$$ \begin{aligned} &\sum_\alpha f_\alpha^{\prime \mathrm{eq}}=\rho^{\prime} \\ &\sum_\alpha\left(\xi_a\right)_a f_\alpha^{\prime \mathrm{eq}}=\bar{\rho} u_a^{\prime}+\rho^{\prime} \bar{u}_a \\ & \sum_\alpha\left(\xi_\alpha\right)_a\left(\xi_\alpha\right)_b f_\alpha^{\prime \text { eq }}=\bar{\rho} u_a^{\prime} \bar{u}_b+\rho^{\prime} \bar{u}_a \bar{u}_b+\bar{\rho} \bar{u}_a u_b^{\prime}+\rho^{\prime} c_{\mathrm{S}}^2 \delta_{a b} \\ & \sum_\alpha\left(\xi_\alpha\right)_a\left(\xi_\alpha\right)_b f_\alpha^{\prime \text { neq }}=-\tau _{0}\delta_t \bar{\rho} c_{\mathrm{S}}^2\left(\nabla_a u_b^{\prime}+\nabla_b u_a^{\prime}\right) \\ & -\tau_{0} \delta_t \rho^{\prime} c_{\mathrm{S}}^2\left(\nabla_a \bar{u}_b+\nabla_b \bar{u}_a\right) \end{aligned} $$ (15) 式中：上标“neq”为粒子分布函数的非平衡部分，与平衡分布函数有关；下标“a”和“b”均可分别取x或${\textit{y}} $，表示粒子速度在x或${\textit{y}} $方向的分量。基于式（15），通过Chapman‒Enskog展开分析，能够证明线化格子Boltzmann方程可以还原为线化N‒S方程，矩关系式建立了介观变量与宏观变量的关系，这也是构造基于Boltzmann模型的通量计算格式的依据。
通量计算首先需要通过插值计算得到网格交界面的宏观变量值，但来自交界面左右单元的插值往往不相等。传统通量计算格式基于黎曼问题获取界面通量，格子Boltzmann通量求解器则根据粒子速度方向判断相应平衡分布函数由左值还是右值计算。将相应位置的宏观量代入式（14）可获取每个粒子对应的平衡分布函数。非平衡部分则采用下式近似：
$$ f_{\alpha}^{'\text{ neq }}=\tau_0\left[f_{\alpha}^{'\rm{eq}}\left(-{{ξ}}_{\alpha}\delta_t,t-\delta_t\right)-f_{\alpha}^{'\rm{eq}}(0,t)\right] $$ (16) 即非平衡分布函数可由界面周围的平衡分布函数确定。最后，结合式（15）即可利用介观变量计算宏观方程的界面通量。
上述通量计算格式基于经典的D2Q9模型，对应的线化N‒S方程不包含能量方程。为了使该通量计算格式更具通用性，进一步发展了D1Q4粒子速度模型[51]，使线化格子Boltzmann方程可以还原到包含能量方程的线化Euler方程。通过Chapman‒Enskog展开分析发现，格子Boltzmann模型除了需要满足式（15）的前3项以外，还需要满足下列矩关系式（$e' $为内能，$ c'=\sqrt{({\text{γ}}−1)e'} $）：
$$ \begin{gathered} \sum\limits_\alpha {{f_\alpha }^{'{\rm{eq}}}{\xi _\alpha }{\xi _\alpha }{\xi _\alpha }} = {{\bar u}_n}\left( {2\bar \rho {{\bar u}_n}{u_n^{\prime}} + \rho '{{\bar u}_n}{{\bar u}_n} + 3p'} \right) \\ + {{u'}_{ n}}\left( {\bar \rho {{\bar u}_n}{{\bar u}_n} + \dfrac{{2\gamma }}{{\gamma - 1}}\bar p - 2\bar \rho e' + \bar \rho {{c'}^2}} \right) \\ \end{gathered} $$ (17) 关于可压缩Boltzmann模型，Qu等[53]基于矩关系反向推导了平衡分布函数表达式，并发展了D1Q4L2、D2Q13L2等粒子速度模型，为推导适用于线化Euler方程通量计算的粒子速度模型提供了很好的借鉴。设D1Q4模型中粒子的平衡分布函数分别为${g'_1}$、${g'_2} $、${g'_3} $、${g'_4} $，粒子速度d1、d2按下式确定：
$$ \left\{ \begin{gathered} {d_2} = \sqrt E \\ {d_1} = {d_2}/4 \\ \end{gathered} \right. $$ (18) 式中，E为总能。将4个平衡分布函数作为未知数代入式（15），可以得到方程组，求解该方程组即可得到粒子对应的平衡分布函数。
获得粒子速度模型后，可根据宏观变量计算交界面的平衡分布函数与粒子速度，并结合矩关系式得到交界面通量，具体过程可参考文献[51]。
2.2 方法验证
通过数个典型声传播算例验证该通量计算格式的正确性与鲁棒性。
2.2.1 声波绕过单个圆柱
算例描述了初始时刻高斯脉冲扰动传播及与壁面的干扰过程，圆柱直径为1，在其周围设置3个监测点A(0, 5)、B(5cos135°, 5sin135°)、C(5, 0)。初始状态下，在原点布置如下高斯脉冲扰动：
$$ \dfrac{{ \partial {{\rho '}_{ {\mathrm{S}}}}}}{{ \partial t}} = \dfrac{{ \partial {{p'}_{ {\mathrm{S}}}}}}{{ \partial t}} = \exp \left[ { - \ln 2\dfrac{{{{\left( {x - 4} \right)}^2} + {{\textit{y}} ^2}}}{{0.04}}} \right] $$ (19) 式中：x、${\textit{y}} $为空间位置坐标，变量均为无量纲量。
图3为无量纲时间t = 4和6时的声波与圆柱干扰的脉动压力云图。图4为3个监测点接收的脉动压力信号及其与解析结果的对比。由图可见，二者吻合很好，验证了算法模拟的声波在含曲面边界传播的正确性。
![]()
图 4 监测点接收的脉动压力信号与解析解的对比Fig. 4 Comparison between analytical solutions and perturbed pressure signals received from monitoring points2.2.2 周期性声源与单个圆柱相互作用
算例模拟高斯周期性声源在圆柱表面发生反射并相互干扰的复杂过程，如图5所示。圆柱直径为1，初始时刻在原点处布置一个周期性声源，描述如下：
$$ \dfrac{{ \partial {{\rho '}_{\rm{ S}}}}}{{ \partial t}} = \dfrac{{ \partial {{p'}_{\rm{ S}}}}}{{ \partial t}} = \exp \left[ { - \ln 2\dfrac{{{{\left( {x - 4} \right)}^2} + {{\textit{y}} ^2}}}{{0.04}}} \right]\sin \left( {\omega t} \right) $$ (20) 将圆柱表面的压力脉动均方根值${p'_{{\rm{rms}}}}$与精确解进行对比，结果如图6所示，二者吻合很好，验证了声波与曲面几何作用的正确性和鲁棒性。
![]()
图 6 圆柱表面脉动压力均方根值与精确解的对比Fig. 6 Comparison of the RMS and exact value of perturbed pressures near the cylinder3 简化的线化格子Boltzmann方法研究
近年来，Chen等[46, 54]引入分步求解技术，提出了一种简化的格子Boltzmann方法（Simplified Lattice Boltzmann Method, SLBM）。该方法通过宏观变量直接计算粒子分布函数的平衡部分，又通过其他位置和时刻的平衡分布函数进行二阶插值近似得到粒子分布函数的非平衡部分，使得粒子分布函数的演化可以转化为宏观量的演化。因此，SLBM不需要储存介观变量，能够节省大量内存，并且可以通过宏观变量施加边界条件。SLBM已经被应用于等熵/等温不可压缩流动[55]、多相流[56]、热流动[57]等问题的模拟中。为了模拟声传播问题，课题组从线化格子Boltzmann方程出发，通过分步求解推导了简化线化格子Boltzmann方法（Simplified Linearized Lattice Boltzmann Method, SLLBM）[58]。
3.1 简化的线化格子Boltzmann方程与算法
对于格子Boltzmann方程，粒子分布函数$f$可分解为平均量$\bar f$与扰动量$f'$，即$f = \bar f + f'$，得到BGK近似的线化格子Boltzmann方程：
$$\dfrac{ \partial f_\alpha^{\prime}}{ \partial t}+{{ξ}}_\alpha \cdot \nabla f_\alpha^{\prime}=-\dfrac{1}{\tau}\left(f_\alpha^{\prime}-f_\alpha^{\prime {\rm{eq}}}\right) $$ (21) 式中：$\tau = \mu / {p} $，$\mu $为动力黏度；$ {f'_\alpha } $为粒子分布函数$f'$沿$\alpha $方向上的分量；$f_\alpha^{\prime {\rm{eq}}}$为其平衡分布函数，表达式如式（14）。对于二维问题，格子Boltzmann模型采用D2Q9模型，其权重系数$ {w_\alpha } $和格子速度$ {{{ξ}}_\alpha } $由下式给出：
$$ \begin{gathered} \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\left| {{{{ξ}}_0}} \right| = 0,}&{\left| {{{{ξ}}_{1 - 4}}} \right| = 1,} \end{array}}&{\left| {{{{ξ}}_{5 - 8}}} \right| = \sqrt 2 } \end{array} \\ \begin{array}{*{20}{c}} {{w_0} = \dfrac{4}{9},}&{{w_{1 - 4}} = \dfrac{1}{9},}&{{w_{5 - 8}} = \dfrac{1}{{36}}} \end{array} \\ \end{gathered} $$ (22) 根据Chapman‒Enskog展开分析，粒子分布函数、时间导数和空间导数分别展开为：
$$ {f'_\alpha } = {f'_\alpha }^{\left( 0 \right)} + K n{f'_\alpha }^{\left( 1 \right)} + K {n^2}{f'_\alpha }^{\left( 2 \right)} $$ (23) $$ \dfrac{ \partial }{{ \partial t}} = K n\dfrac{ \partial }{{ \partial {t_0}}} + K {n^2}\dfrac{ \partial }{{ \partial {t_1}}} $$ (24) $$ \nabla = K n{\nabla _1} $$ (25) 式中：$K n$为Knudsen数。将上式代入式（21）中，可以还原线化Euler方程：
$$ \dfrac{{ \partial \rho '}}{{ \partial t}} + \nabla \cdot \left( {\sum\limits_\alpha {{{{ξ}}_\alpha }{{f'_{ \alpha}}^{{\rm{eq}}}} } } \right) = 0 $$ (26) $$ \dfrac{\partial \left(\rho^{\prime} \overline{\boldsymbol{u}} + \bar{\rho} \boldsymbol{u}^{\prime}\right)}{\partial t}+\nabla \cdot \left[ \sum_\alpha {ξ}_\alpha {ξ}_\alpha f_\alpha^{\prime {\rm{eq}}} + \left( 1 - \dfrac{1}{2 \tau} \right) \sum_\alpha {ξ}_\alpha {ξ}_\alpha f_\alpha^{\prime {{\mathrm{neq}} }} \right] = 0 $$ (27) 式中：非平衡分布函数$f_\alpha ^{{\rm{neq}}} = {f_\alpha } - f_\alpha ^{{\rm{eq}}} = K nf_\alpha ^{\left( 1 \right)}$，且满足以下矩关系 ：
$$ \sum_\alpha f_\alpha^{\prime \text { neq }}=0, \sum_\alpha {ξ}_\alpha {f'}_{ \alpha}^{{{\mathrm{neq}} }}=0$$ (28) 根据分步计算法，求解可以分成预估和校正2步。预估步为：
$$ \dfrac{{ \partial \rho '}}{{ \partial t}} + \nabla \cdot \left( {\sum\limits_\alpha {{{{ξ}}_\alpha }{{f'_\alpha }^{{\rm{eq}}}}} } \right) = 0 $$ (29) $$ \dfrac{ \partial\left(\rho^{\prime} \overline{\boldsymbol{u}}+\bar{\rho} \boldsymbol{u}^{\prime}\right)}{ \partial t}+\nabla \cdot\left(\sum_\alpha {ξ}_\alpha {ξ}_\alpha f_\alpha^{\prime {\rm{eq}}}+\dfrac{1}{2 \tau} \sum_\alpha {ξ}_\alpha {ξ}_\alpha {f'}_{ \alpha}^{ {{\mathrm{neq}} }} \right) = 0 $$ (30) 校正步为：
$$ \dfrac{{ \partial \rho '}}{{ \partial t}} = 0 $$ (31) $$ \dfrac{ \partial\left(\rho^{\prime} \overline{\boldsymbol{u}}+\bar{\rho} \boldsymbol{u}^{\prime}\right)}{ \partial t}+\nabla \cdot\left[\left(1-\dfrac{1}{\tau}\right) \sum_\alpha {ξ}_\alpha {ξ}_\alpha f_\alpha^{\prime { {\mathrm{neq}} }}\right]=0 $$ (32) 预估步可以通过下式求解：
$$ \rho^{\prime *} =\sum_\alpha f_\alpha^{\prime {\rm{eq}}}(\boldsymbol{X}-{ξ}_\alpha \delta_t, t-\delta_t) $$ (33) $$ \bar{\rho} \boldsymbol{u}^{\prime *}+\rho^{\prime *} \overline{\boldsymbol{u}} =\sum_\alpha {ξ}_\alpha f_\alpha^{\prime {\rm{eq}}}(\boldsymbol{X}-{ξ}_\alpha \delta_t, t-\delta_t) $$ (34) 式中：${\delta_t}$为计算过程中的时间步长；上标“$*$”表示通过求解预估步得到的宏观变量的中间值。校正步与预估步的求解类似，将矩关系带入非平衡分布函数$f_\alpha^{\prime {\rm{neq}}}(\boldsymbol{X}-{ξ}_\alpha \delta_t, t)$的泰勒展开中可得：
$$ \begin{aligned} &-\nabla \cdot\left[\left(1-\dfrac{1}{\tau}\right) \sum_\alpha {ξ}_\alpha {ξ}_\alpha {f'_\alpha}^{ {{\mathrm{neq }}}}(\boldsymbol{X}, t)\right]\\ &=\left(1-\dfrac{1}{\tau}\right) \dfrac{1}{\delta_t} \sum_\alpha {ξ}_\alpha {f'_\alpha}^{{{\mathrm{neq}} }}\left(\boldsymbol{X}-{ξ}_\alpha \delta_t, t\right) \end{aligned}$$ (35) 将上式代入式（32）中，可通过下式求解校正步：
$$ {\rho '^{n + 1}} = {\rho '^{\text{*}}} $$ (36) $$ \begin{aligned} &\bar{\rho} \boldsymbol{u}^{\prime n+1}+\rho^{\prime n+1} \overline{\boldsymbol{u}}=\bar{\rho} \boldsymbol{u}^{\prime *}+\rho^{\prime *} \overline{\boldsymbol{u}}+\\ &\left(1-\dfrac{1}{\tau}\right) \sum_\alpha {ξ}_\alpha f_\alpha^{\prime {\rm{neq}}}\left(\boldsymbol{X}-{ξ}_\alpha \delta_t, t\right) \end{aligned}$$ (37) 至此得到简化的线化格子Boltzmann方法。
3.2 方法验证
利用上述方法对均匀介质中的高斯脉冲扰动传播及其与壁面的反射过程进行了模拟，结果如图7所示。高斯脉冲扰动传播的计算域为$\left[ { - 100,100} \right] \times \left[ { - 100,100} \right]$，计算域内均匀布置网格点，网格尺度${{\text{δ}}_x} = {{\text{δ}}_{\textit{y}}} = 1.0$，时间步长${\delta_t} = 1.0$。初始时刻在计算域中心施加的高斯脉冲扰动由下式给出：
$$ \left\{ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\rho ' = \varepsilon \,\exp \left[ { - \beta \left( {{x^2} + {{\textit{y}} ^2}} \right)} \right]},\;{\bar \rho = {\rho _0}} \end{array}} \\ {\begin{array}{*{20}{c}} { u' = 0},{\bar u = {u_0}} \end{array}} \\ {\begin{array}{*{20}{c}} { v' = 0},{\bar v = {v_0}} \end{array}} \end{array}} \right. $$ (38) 式中：$ {\rho _0} = 1.0 $，$\varepsilon = 0.01$，$\beta = \left( {\ln 2} \right)/{8^2}$。图7为$t = 80$时刻${u_0} = 0$和$0.3$的均匀流场的瞬时扰动密度云图，图8为计算域中心线${\textit{y}} = 0$处的扰动密度曲线与精确解的对比。可以看出，SLLBM的计算结果均与精确解吻合很好。
进一步模拟壁面对声波的反射作用，声源位置为(0, −75)，$\beta = \left( {\ln 2} \right)/{3^2}$，${u_0} = 0$，网格尺度${{\text{δ}}_x} = {{\text{δ}}_{\textit{y}} } = 0.5$，时间步长${\delta_t} = 0.5$。图9为$ t=26\text{、}120\text{、} 160 $时刻的瞬时扰动密度云图，图10为相同时刻下${\textit{y}} = - 100$处的瞬时扰动密度曲线与精确解的对比。从图中可以看出，SLLBM能够很好地模拟壁面对声波的反射作用，且计算结果与精确解吻合很好。
![]()
图 9 t = 26、120、160时刻考虑壁面反射时的瞬时脉动密度云图Fig. 9 Instantaneous perturbed density contours with wall reflection at t = 26、120 and 160![]()
图 10 不同时刻下y = −100处的脉动密度分布对比Fig. 10 Instantaneous pertrubed density distribution at y = −100 in different moments4 基于线化方程的声传播计算应用
剪切层普遍存在于开口式、闭口式与凯夫拉式声学风洞中，其折射效应会改变声传播轨迹，影响声源定位的精度[59-60]。基于线化方程的声传播计算方法不仅可用于剪切层对声传播特性的影响分析，还可用于检验剪切层折射理论的准确性，此外，还能够计算获得剪切层引起的定位误差并构建数学模型，为提高声源定位精度提供参考。
当剪切层厚度不可忽略时，计算发现：厚度对声波垂直入射区域附近的声场影响较小，对垂直入射区两侧的影响较大[61]。声波垂直入射位置可用${x_{{\text{nom}}}} = {x_{{\text{S}}}} + {M a_{{{\infty}} }} h$进行估计，其中$ x_{\text{S}} $为声源横坐标，${M a_{{{\infty}} }}$为射流马赫数，$ h $为声源与剪切层的距离。在风洞实验中，可将阵列布置于声源正上方，以减弱剪切层对测量结果的影响。
为考虑实际剪切层对声传播的影响规律，将包含可压缩性与温度修正的Görtler自相似速度分布作为背景流场[48-49, 61]：
$$ \left\{ \begin{array}{l} \dfrac{{{{\bar u}} - {{\bar u}_{{\text{air}}}}}}{{{{\bar u}_{\text{jet}}} - {{\bar u}_{\text{air}}}}} = \dfrac{1}{2}\left[ {1 + {f _{{\rm{er}}}} \left( {{\text{β}} + {\text{β} _0}} \right)} \right] \\ \dfrac{{{\bar v} - {{\bar v}_{\text{air}}}}}{{{{\bar u}_{\text{jet}}} - {{\bar u}_{\text{air}}}}} = \dfrac{1}{{2\sigma }}\Big\{ {\text{β} _0}{f _{{\rm{er}}}} \left( {\text{β} + {\text{β} _0}} \right) + \Big.\\ \qquad\qquad\qquad\Big.\dfrac{1}{{\sqrt \pi }}\exp \left[ { - {{\left( {\text{β} + {\text{β} _0}} \right)}^2}} \right]\Big\} - \dfrac{{{\text{β} _0}}}{{2\sigma }} \end{array} \right. $$ (39) 式中：fer为误差函数；$\sigma $为扩张率；下标“jet”和“air”分别表示射流和空气。其余变量为：
$$ \text{β}=-\dfrac{\sigma\left(y -y _{\rm{S}h}\right)}{x_{ }-x_{\rm{S}h}} $$ (40) $$ \sigma = \dfrac{{{\sigma _0}}}{{\varPhi {\lambda _{\rm{S}}}}} $$ (41) $$ \lambda_{\rm{S}}\text{ }=\dfrac{\left(1-\overline{u}_{\text{air}}/\overline{u}_{\text{jet}}\right)\left(1+\sqrt{\overline{\rho}_{\text{air}}/\overline{\rho}_{\text{jet}}}\right)}{2\left(1+\overline{u}_{\text{air}}/\overline{u}_{\text{jet}}\sqrt{\overline{\rho}_{\text{air}}/\overline{\rho}_{\text{jet}}}\right)}=\dfrac{1+\mu'}{2} $$ (42) $$ \varPhi = 1 + 0.546\,8\left( {\dfrac{1}{{1 + 35.53M a_c^{8.614}}} - 1} \right) $$ (43) $$ M a_c=\dfrac{\overline{u}_{\text{jet}}-\overline{u}_{\text{air}}}{c_{\text{jet}}+c_{\text{air}}}=\dfrac{M a_{\text{jet}}}{\mu'+1} $$ (44) 式中：$ \overline{u}_{\text{air}}=0 $，$ (x_{\rm{S}h},\, y_{\rm{S}h}) $为剪切层起点；声速比$\mu'$ = cjet/cair；不可压剪切层的扩张率为9～15，在此基础上考虑可压缩性与温度后的扩张率为$\sigma$；${M a_c}$为对流马赫数；$\varPhi$为可压缩性修正系数；${\lambda _{\rm{S}}}$表征速度和密度对扩张率的影响。假设空气为理想气体、压力为常数、普朗特常数为1，则由Crocco‒Busemann关系可获得剪切层密度。
如图11所示，计算域尺寸为1.4 m × 1.4 m，各方向的网格点数为501，四周设置无反射边界条件。声源源项为：
$$ \left\{ \begin{gathered} \dfrac{{ \partial {{\rho '}_{ {\rm{S}}}}}}{{ \partial t}} = 0 \\ \dfrac{{ \partial {{p'}_{ {\rm{S}}}}}}{{ \partial t}} = \varepsilon\, \exp \left\{ { - \dfrac{{\ln 2}}{{{{(3{\text{Δ}})}^2}}} \left[ {{{({x} - {x_{{\mathrm{S}}}})}^2} + {{({y} - {y_{{\mathrm{S}}}})}^2}} \right]} \right\}\sin (\omega t) \\ \end{gathered} \right. $$ (45) 式中：${\text{Δ}}= 0.005$，$\varepsilon = 0.01$，声源坐标$ (x_{\rm{S}},y_{\mathrm{S}})= (0,-0.3) $。
使用改进梯度项抑制方法计算不同剪切流场下的声传播特性，结果如图12所示。由图12（a）可知，低速射流中声源产生的声波分别以$1 \mp {M a_{\text{jet}}}$的速度向上、下游传播，在接触剪切层后，发生折射与反射。当射流马赫数增至0.8后，图12（b）中出现类似“马赫锥”的现象。此时声源上游的声波波长较短，衰减较快，且剪切层上方的折射效应更明显。保持射流速度不变，升高温度则会增大剪切层扩张角，同时削弱剪切层强度的影响，如图12（c）所示。
![]()
基于线化方程的声传播计算方法可以用来检验理论预测无限薄剪切流中声波折射角与幅值的准确性。由图13可知，修正方法和改进梯度项抑制方法的结果接近。折射角从上游至下游逐渐减小，表明剪切层上方的声线都往左偏折，且互不相交。马赫数越大，剪切层强度越强，故折射角越大。射流速度和声源频率固定时，高温射流会增大声速，从而间接减小对流马赫数和折射角。在不同温度比下，修正方法和改进梯度项抑制方法的幅值比也基本相同，如图14所示。幅值比在高马赫数和上游位置存在最大值，在低马赫数和下游位置存在最小值。这是因为剪切流上方声线都向上游偏折，且随着射流温度升高，折射效果削弱。总体来看，修正方法对无限薄剪切流是有效的。
![]()
![]()
以基于线化方程的声传播计算方法获得的远场声信号作为延迟求和波束成形方法的输入，可以定位出声源的位置。对比定位声源与实际声源，得到定位误差${x_{{\rm{shift}}}}$。处理和分析定位误差数据，可以建立如图15和式（46）所示的线性拟合模型：
$$ {x_{{\text{shift}}}} = C{M a_{{\text{jet}}}}Sr $$ (46) 式中：C为拟合得到的常数，改进梯度项抑制方法的C =
0.2368 ，传统梯度项抑制方法的C =0.2049 ；${M a_{\text{jet}}}$为射流马赫数；斯特劳哈尔数$Sr = f{l_{{\text{mean}}}}/ {c_{{\text{air}}}} = {l_{{\text{mean}}}}/\lambda$，其中${l_{{\rm{mean}}}}$为剪切层的平均厚度，$f$为声源频率，${c_{{\rm{air}}}}$为远场声速，$\lambda $为声波波长。这一拟合公式定量反映了声源定位误差随射流马赫数和斯特劳哈尔数的变化关系。5 结 论
本文主要介绍了课题组在基于线化方程的声传播计算方面的研究工作，总结如下：
1）针对基于线化方程模拟声波在剪切层中传播时出现的数值不稳定波问题，提出了改进梯度项抑制方法。该方法可以有效过滤不稳定波，保证数值计算的稳定性，而且能够较准确地预测远场声压。该方法求解思路简单，没有调节参数，具有很好的鲁棒性。
2）针对采用有限体积法求解线化方程模拟声传播面临的通量计算问题，发展了基于格子Boltzmann模型的线化方程通量计算格式。该算法继承了传统基于Boltzmann模型通量计算格式的优点，从介观角度构造通量计算格式，通过局部求解线化格子Boltzmann方程，使得通量计算更符合物理。
3）针对线化格子Boltzmann方程模拟声传播时存在的内存占用大与边界条件处理难的问题，提出了一种简化的线化格子Boltzmann方法。该方法通过直接对宏观量进行演化，能够节省计算内存，也便于施加边界条件。
4）研究了声波穿过剪切层传播的规律，建立了定位误差随射流马赫数与斯特劳哈尔数变化的数学模型，为风洞实验声源定位修正提供了参考。
-
![]()
![]()
![]()
图 4 监测点接收的脉动压力信号与解析解的对比
Fig. 4 Comparison between analytical solutions and perturbed pressure signals received from monitoring points
![]()
图 6 圆柱表面脉动压力均方根值与精确解的对比
Fig. 6 Comparison of the RMS and exact value of perturbed pressures near the cylinder
![]()
图 9 t = 26、120、160时刻考虑壁面反射时的瞬时脉动密度云图
Fig. 9 Instantaneous perturbed density contours with wall reflection at t = 26、120 and 160
![]()
图 10 不同时刻下y = −100处的脉动密度分布对比
Fig. 10 Instantaneous pertrubed density distribution at y = −100 in different moments
![]()
![]()
![]()
-
[1] 张树海, 李虎, 王益民. 含激波、旋涡和声波的复杂多尺度流动数值模拟研究[J]. 空气动力学学报, 2018, 36(3): 449–462. ZHANG S H, LI H, WANG Y M. Numerical study for complex multi-scale flow with shock, vortex and sound wave[J]. Acta Aerodynamica Sinica, 2018, 36(3): 449–462.
[2] COLONIUS T, LELE S K, MOIN P. The scattering of sound waves by a vortex: numerical simulations and analytical solutions[J]. Journal of Fluid Mechanics, 1994, 260: 271–298. doi: 10.1017/s0022112094003514
[3] HONG Z L, WANG X Y, JING X D, et al. Vortex sound interaction in acoustic resonance of a flow duct containing a plate[J]. Journal of Sound and Vibration, 2020: 115482. doi: 10.1016/j.jsv.2020.115482
[4] CLAIR V, GABARD G. Spectral broadening of acoustic waves by convected vortices[J]. Journal of Fluid Mechanics, 2018, 841: 50–80. doi: 10.1017/jfm.2018.94
[5] DU Y L, YU H W, LIU Y C, et al. Hydrodynamic sources of the vortex sound in a two-dimensional shear layer[J]. International Journal of Aeroacoustics, 2023, 22: 41–59. doi: 10.1177/1475472X221150177
[6] MA R X, ZHANG S H, LUO Y, et al. Numerical study on multiple acoustic scattering by a vortex array[J]. Journal of Sound and Vibration, 2022, 527: 116815. doi: 10.1016/j.jsv.2022.116815
[7] BAILLY C, JUVÉ D. Numerical solution of acoustic propagation problems using linearized Euler equations[J]. AIAA Journal, 2000, 38(1): 22–29. doi: 10.2514/2.949
[8] 李晓东, 江旻, 高军辉, 等. 计算气动声学进展与展望[J]. 中国科学: 物理学 力学 天文学, 2014, 44(3): 234-248. LI X D, JIANG M, GAO J H, et al. Progress and prospective of computational aeroacoustics[J]. Scientia Sinica (Physica, Mechanica & Astronomica), 2014, 44(3): 234-248.
[9] 吕宏强, 朱国祥, 宋江勇, 等. 线化欧拉方程的高阶间断有限元数值解法研究[J]. 力学学报, 2011, 43(3): 621–624. DOI: 10.6052/0459-1879-2011-3-lxxb2010-077 LYU H Q, ZHU G X, SONG J Y, et al. High-order discontinuous Galerkin solution of linearized Euler equations[J]. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(3): 621–624. doi: 10.6052/0459-1879-2011-3-lxxb2010-077
[10] 余培汛, 白俊强, 杨海, 等. 耦合RANS/LES模型与LEE方程的气动噪声混合预测方法[J]. 西北工业大学学报, 2017, 35(6): 990–997. DOI: 10.3969/j.issn.1000-2758.2017.06.009 YU P X, BAI J Q, YANG H, et al. Hybrid aero-acoustics prediction method coupled with RANS/LES model and LEE equation[J]. Journal of Northwestern Polytechnical University, 2017, 35(6): 990–997. doi: 10.3969/j.issn.1000-2758.2017.06.009
[11] XIN B, SUN D K, JING X D, et al. Numerical study of acoustic instability in a partly lined flow duct using the full linearized Navier–Stokes equations[J]. Journal of Sound and Vibration, 2016, 373: 132–146. doi: 10.1016/j.jsv.2016.02.042
[12] BOGEY C, BAILLY C, JUVÉ D. Computation of flow noise using source terms in linearized Euler’s equations[J]. AIAA Journal, 2002, 40(2): 235–243. doi: 10.2514/2.1665
[13] ZHANG X, CHEN X X, GILL J R. Gradient term filtering for stable sound propagation with linearized Euler equations[C]//Proc of the 20th AIAA/CEAS Aeroacoustics Conference. 2014: 3306. doi: 10.2514/6.2014-3306
[14] EWERT R, SCHRÖDER W. Acoustic perturbation equations based on flow decomposition via source filtering[J]. Journal of Computational Physics, 2003, 188(2): 365–398. doi: 10.1016/S0021-9991(03)00168-2
[15] BAUER M, DIERKE J, EWERT R. Application of a discontinuous Galerkin method to discretize acoustic perturbation equations[J]. AIAA Journal, 2011, 49(5): 898–908. doi: 10.2514/1.J050333
[16] EWERT R, SCHRÖDER W. On the simulation of trailing edge noise with a hybrid LES/APE method[J]. Journal of Sound and Vibration, 2004, 270(3): 509–524. doi: 10.1016/j.jsv.2003.09.047
[17] SEO J H, MOON Y J. Linearized perturbed compressible equations for low Mach number aeroacoustics[J]. Journal of Computational Physics, 2006, 218(2): 702–719. doi: 10.1016/j.jcp.2006.03.003
[18] HU F Q, LI X D, LI X Y, et al. Time domain wave packet method and suppression of instability waves in aeroacoustic computations[J]. Journal of Fluids Engineering, 2014, 136(6): 060905. doi: 10.1115/1.4025866
[19] 李旦望, 李晓东, 李小艳, 等. 管道后传声数值模拟的不稳定波抑制[J]. 工程热物理学报, 2012, 33(4): 587–590. LI D W, LI X D, LI X Y, et al. Suppression of instability waves in numerical simulations of sound propagation from aft ducts[J]. Journal of Engineering Thermophysics, 2012, 33(4): 587–590.
[20] PRAX C, GOLANSKI F, NADAL L. Control of the vorticity mode in the linearized Euler equations for hybrid aeroacoustic prediction[J]. Journal of Computational Physics, 2008, 227(12): 6044–6057. doi: 10.1016/j.jcp.2008.02.022
[21] TESTER B, GABARD G, ÖZYÖRÜK Y. Influence of mean flow gradients on fan exhaust noise predictions[C]//Proc of the 14th AIAA/CEAS Aeroacoustics Conference (29th AIAA Aeroacoustics Conference). 2008: 2825. doi: 10.2514/6.2008-2825
[22] DENG Y Y, ALOMAR A, DRAGNA D, et al. Characterization and suppression of the hydrodynamic instability in the time domain for acoustic propagation in a lined flow duct[J]. Journal of Sound and Vibration, 2021, 500: 115999. doi: 10.1016/j.jsv.2021.115999
[23] XU K. A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method[J]. Journal of Computational Physics, 2001, 171(1): 289–335. doi: 10.1006/jcph.2001.6790
[24] Guo Z L, Wang L P, Qi Y M. Discrete unified gas kinetic scheme for continuum compressible flows[J]. Physical Review E, 2023, 107(2): 025304. doi: 10.1103/PhysRevE.107.025304
[25] Zhong M L, Zou S, Pan D X, et al. A simplified discrete unified gas-kinetic scheme for compressible flow[J]. Physics of Fluids, 2021, 33(3): 036103. doi: 10.1063/5.0033911
[26] SHU C, WANG Y, TEO C J, et al. Development of lattice boltzmann flux solver for simulation of incompressible flows[J]. Advances in Applied Mathematics and Mechanics, 2014, 6(4): 436–460. doi: 10.4208/aamm.2014.4.s2
[27] WANG Y, YANG L M, SHU C. From lattice boltzmann method to lattice boltzmann flux solver[J]. Entropy, 2015, 17(12): 7713–7735. doi: 10.3390/e17117713
[28] 杨鲤铭. 基于格子和连续Boltzmann模型的通量求解器研究及其应用[D]. 南京: 南京航空航天大学, 2016. YANG L M. Research and application of flux solver based on lattice and continuous Boltzmann model[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2016.
[29] ZHU Y J, ZHONG C W, XU K. GKS and UGKS for high-speed flows[J]. Aerospace, 2021, 8(5): 141. doi: 10.3390/aerospace8050141
[30] PAN D X, ZHONG C W, ZHUO C S. An implicit gas-kinetic scheme for turbulent flow on unstructured hybrid mesh[J]. Computers & Mathematics with Applications, 2018, 75(11): 3825–3848. doi: 10.1016/j.camwa.2018.02.032
[31] LI J, ZHONG C W, WANG Y, et al. Implementation of dual time-stepping strategy of the gas-kinetic scheme for unsteady flow simulations[J]. Physical Review E, 2017, 95(5): 053307. doi: 10.1103/physreve.95.053307
[32] ZHANG C, LI Q B, FU S, et al. A third-order gas-kinetic CPR method for the Euler and Navier–Stokes equations on triangular meshes[J]. Journal of Computational Physics, 2018, 363: 329–353. doi: 10.1016/j.jcp.2018.02.040
[33] CHEN S Z, GUO Z L, XU K. Simplification of the unified gas kinetic scheme[J]. Physical Review E, 2016, 94(2): 023313. doi: 10.1103/physreve.94.023313
[34] YANG L M, SHU C, WU J. Extension of lattice Boltzmann flux solver for simulation of 3D viscous compressible flows[J]. Computers & Mathematics with Applications, 2016, 71(10): 2069–2081. doi: 10.1016/j.camwa.2016.03.027
[35] YANG Y Q, PAN L, XU K. High-order gas-kinetic scheme on three-dimensional unstructured meshes for compressible flows[J]. Physics of Fluids, 2021, 33(9): 09610. doi: 10.1063/5.0062368
[36] ZHANG C, LI Q B. A third-order subcell finite volume gas-kinetic scheme for the Euler and Navier-Stokes equations on triangular meshes[J]. Journal of Computational Physics, 2021, 436: 110245. doi: 10.1016/j.jcp.2021.110245
[37] ZHAN N Y, CHEN R Q, YOU Y C. Meshfree method based on discrete gas-kinetic scheme to simulate incompressible/compressible flows[J]. Physics of Fluids, 2021, 33(1): 017112. doi: 10.1063/5.0033770
[38] LIU Y Y, SHU C, ZHANG H W, et al. A high order least square-based finite difference-finite volume method with lattice Boltzmann flux solver for simulation of incompressible flows on unstructured grids[J]. Journal of Computational Physics, 2020, 401: 109019. doi: 10.1016/j.jcp.2019.109019
[39] MARIÉ S, RICOT D, SAGAUT P. Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[J]. Journal of Computational Physics, 2009, 228(4): 1056–1070. doi: 10.1016/j.jcp.2008.10.021
[40] LATT J, MALASPINAS O, KONTAXAKIS D, et al. Palabos: parallel lattice boltzmann solver[J]. Computers & Mathematics with Applications, 2021, 81: 334–350. doi: 10.1016/j.camwa.2020.03.022
[41] CASALINO D, HAZIR A, MANN A. Turbofan broadband noise prediction using the lattice boltzmann method[J]. AIAA Journal, 2017, 56(2): 609–628. doi: 10.2514/1.J055674
[42] SCHÄFER R, BÖHLE M. Validation of the lattice boltzmann method for simulation of aerodynamics and aeroacoustics in a centrifugal fan[J]. Acoustics, 2020, 2(4): 735–752. doi: 10.3390/acoustics2040040
[43] 司海青, 朱卫军. 气动噪声计算方法及其应用[M]. 北京: 科学出版社, 2017. SI H Q, ZHU W J. Calculation method of aerodynamic noise and its application[M]. Beijing: Science Press, 2017.
[44] VERGNAULT E, MALASPINAS O, SAGAUT P. A lattice Boltzmann method for nonlinear disturbances around an arbitrary base flow[J]. Journal of Computational Physics, 2012, 231(24): 8070–8082. doi: 10.1016/j.jcp.2012.07.021
[45] PÉREZ J M, AGUILAR A, THEOFILIS V. Lattice Boltzmann methods for global linear instability analysis[J]. Theoretical and Computational Fluid Dynamics, 2017, 31(5): 643–664. doi: 10.1007/s00162-016-0416-7
[46] CHEN Z, SHU C, WANG Y, et al. A simplified lattice boltzmann method without evolution of distribution function[J]. Advances in Applied Mathematics and Mechanics, 2017, 9(1): 1–22. doi: 10.4208/aamm.oa-2016-0029
[47] VIGGEN E M. The lattice Boltzmann method with applications in acoustics[D]. Trondheim: Norwegian University of Science and Technology, 2009.
[48] 王李璨. 剪切流中声波传播问题的计算方法改进与影响规律研究[D]. 厦门: 厦门大学, 2021. WANG L C. Computational method improvement and influential assessment about sound propagation problems in A shear layer[D]. Xiamen: Xiamen University, 2021.
[49] WANG L C, CHEN R Q, YOU Y C, et al. Investigation of acoustic propagation and source localization in a hot jet flow[J]. Journal of Sound and Vibration, 2021, 492: 115801. doi: 10.1016/j.jsv.2020.115801
[50] ZHAN N Y, CHEN R Q, YOU Y C. Linear lattice Boltzmann flux solver for simulating acoustic propagation[J]. Computers & Mathematics with Applications, 2022, 114: 21–40. doi: 10.1016/j.camwa.2022.03.034
[51] ZHAN N Y, CHEN R Q, SONG Q C, et al. Linear discrete velocity model-based lattice Boltzmann flux solver for simulating acoustic propagation in fluids[J]. Physical Review E, 2022, 105(6): 065303. doi: 10.1103/physreve.105.065303
[52] QIAN Y H, D’HUMIÈRES D, LALLEMAND P. Lattice BGK models for Navier-Stokes equation[J]. Europhysics Letters (EPL), 1992, 17(6): 479–484. doi: 10.1209/0295-5075/17/6/001
[53] QU K, SHU C, CHEW Y T. Simulation of shock-wave propagation with finite volume lattice boltzmann method[J]. International Journal of Modern Physics C, 2007, 18(4): 447–454. doi: 10.1142/S012918310701067X
[54] CHEN Z, SHU C, TAN D. Highly accurate simplified lattice Boltzmann method[J]. Physics of Fluids, 2018, 30(10): 103605. doi: 10.1063/1.5050185
[55] CHEN Z, SHU C, TAN D. Immersed boundary-simplified lattice Boltzmann method for incompressible viscous flows[J]. Physics of Fluids, 2018, 30(5): 053601. doi: 10.1063/1.5028353
[56] YANG L M, SHU C, CHEN Z, et al. An improved multiphase lattice Boltzmann flux solver for the simulation of incompressible flow with large density ratio and complex interface[J]. Physics of Fluids, 2021, 33(3): 033306. doi: 10.1063/5.0038617
[57] LU J H, DAI C S, YU P. Analysis and reconstruction of the simplified thermal lattice Boltzmann method[J]. International Journal of Heat and Mass Transfer, 2022, 187: 122576. doi: 10.1016/j.ijheatmasstransfer.2022.122576
[58] SONG Q C, CHEN R Q, CAO S Q, et al. A simplified linearized lattice boltzmann method for acoustic propagation simulation[J]. Entropy, 2022, 24(11): 1622. doi: 10.3390/e24111622
[59] 张军, 王勋年, 张俊龙, 等. 开口风洞声阵列测量的剪切层修正方法[J]. 实验流体力学, 2018, 32(4): 39–46. DOI: 10.11729/syltlx20180013 ZHANG J, WANG X N, ZHANG J L, et al. Shear layer correction methods for open-jet wind tunnel phased array test[J]. Journal of Experiments in Fluid Mechanics, 2018, 32(4): 39–46. doi: 10.11729/syltlx20180013
[60] 张雪, 陈宝, 卢清华. Amiet剪切层理论的角度折射验证研究[J]. 应用声学, 2014, 33(5): 433–438. DOI: 10.11684/j.issn.1000-310X.2014.05.009 ZHANG X, CHEN B, LU Q H. Verification of angle refraction correction based on Amiet’s shear layer theory[J]. Journal of Applied Acoustics, 2014, 33(5): 433–438. doi: 10.11684/j.issn.1000-310X.2014.05.009
[61] 王李璨, 陈荣钱, 尤延铖, 等. 剪切层形态对声波传播与声源定位的影响研究[J]. 西北工业大学学报, 2019, 37(6): 1148–1157. DOI: 10.3969/j.issn.1000-2758.2019.06.008 WANG L C, CHEN R Q, YOU Y C, et al. Effects of shear layer characteristics on acoustic propagation and source localization[J]. Journal of Northwestern Polytechnical University, 2019, 37(6): 1148–1157. doi: 10.3969/j.issn.1000-2758.2019.06.008
下载:
计量
- 文章访问数: 354
- HTML全文浏览量: 44
- PDF下载量: 44
