计算气动声学近年来随着航空工业需求而快速发展，主要目的是研究和模拟气动噪声的产生和传播过程，从而为降噪方案设计提供指导.

格子Boltzmann方法（lattice Boltzmann method，LBM）的低耗散和低色散误差优点^[1]使其在计算气动声学领域逐步得到应用. Kam等^[2]采用改进的LBM方法在均匀网格中模拟声波的扰动及散射问题，Viggen^[3-4]应用LBM在正交均匀网格中对周期性声源和平面波的传播进行数值模拟，计算结果与解析解较吻合. 王勇等^[5]应用LBM方法在均匀网格中模拟平面波的衰减过程，司海青等^[6-7]在均匀等距网格应用LBM模拟高斯脉冲传播和方形圆柱噪声问题. 目前LBM方法模拟气动噪声问题主要采用均匀网格，而在非均匀网格和贴体网格中对气动噪声模拟的研究还较少. 邵卫东等^[8]发展伽辽金玻尔兹曼方法在非均匀网格中对声传播问题进行模拟，李凯^[9]采用有限体积可压缩LBM在贴体网格下对圆柱和翼型绕流的气动噪声进行模拟，此外还有学者^[10-12]采用浸入边界法及大涡LBM耦合方法处理曲面边界的流动和噪声问题. 近年来有学者^[13-14]提出通用的插值补充格子Boltzmann方法（generalized form of interpolation supplemented lattice Boltzmann method，GILBM）可以在贴体网格下对圆柱这类复杂外形的流动问题进行模拟. 相比浸入边界法，该方法可以处理贴体网格，从而能够更好模拟附面层内的流动信息，在流动问题模拟中可以应用，然而尚未有学者将该方法应用于非均匀网格下的气动声学问题模拟.

本研究通过顶盖驱动方腔流算例验证GILBM方法的正确性，在此基础上应用该方法对非均匀网格下的二维高斯脉冲、周期性声源传播以及低雷诺数下圆柱绕流涡脱落产生的气动噪声问题进行模拟，研究GILBM方法在非均匀网格和贴体网格中计算复杂外形噪声产生和传播的能力.

基于Bhatnagar–Gross–Krook（BGK）模型的LBM控制方程为

(1) ${f_i}({{x}} + {{{c}}_i}\Delta t,t + \Delta t) - {f_i}({{x}},t) = - \frac{1}{\tau }\left[ {{f_i}({{x}},t) - f_i^{(0)}({{x}},t)} \right].$

式中： ${f_i}({{x}},t)$为密度分布函数方程， ${{x}}$为粒子所在位置， $t$为时间，i为粒子碰撞与迁移的方向， $i = 0,\cdots,8$； ${{{c}}_i}$为离散速度； $\Delta t$为时间步长， $f_i^{(0)}({{x}},t)$为平衡分布函数； $\tau $为松弛时间，取决于运动黏性 $\upsilon $. $\tau $的定义为

(2) $\tau = \frac{{3\upsilon }}{{{c^2}\Delta t}} + \frac{1}{2}.$

式中： $c$为粒子运动速度，LBM-BGK控制方程在时间及空间上具有二阶计算精度^[15].

离散速度模型采用D2Q9模型， ${{{c}}_i}$定义为

(3) ${ c}_i=c\left[ {\begin{array}{*{20}{c}} 0&1&0&{ - 1}&0&1&{ - 1}&{ - 1}&1 \\ 0&0&1&0&{ - 1}&1&1&{ - 1}&{ - 1} \end{array}} \right].$

LBM控制方程求解主要包括碰撞和迁移过程，即将式（1）分解为

(4) $f_i^*({{x}},t) = {f_i}({{x}},t) - \frac{1}{\tau }\left[ {{f_i}({{x}},t) - f_i^{(0)}({{x}},t)} \right],$

(5) ${f_i}({{x}},t + \Delta t) = f_i^*({{x}} - {{{c}}_i}\Delta t,t).$

其中平衡分布函数 $f_i^{(0)}$表达式为

(6) $f_i^{(0)}({{x}},t) = {w_i}\rho \left[ {1 + \frac{{3\left( {{{{c}}_i}\cdot{{u}}} \right)}}{{{c^2}}} + \frac{{9{{\left( {{{{c}}_i}\cdot{{u}}} \right)}^2}}}{{2{c^4}}} - \frac{{3{\left({{u}}\cdot{{u}}\right)}}}{{2{c^2}}}} \right].$

式中：ρ为宏观密度；u为宏观速度； ${w_i}$为权重系数，取决于离散速度模型.

GILBM方法从补充差值格子Boltzmann方法（interpolation-supplemented lattice Boltzmann method，ISLBM）^[16]发展而来，以二维为例，在计算过程中坐标系的转换由物理平面 ${{x}} \equiv \left( {x_1,x_2} \right)$转化为计算平面 ${{\xi }} \equiv \left( {\xi_1 ,\xi_2 } \right)$. 式（4）的 ${{x}}$被替换成 ${{\xi }}$：

(7) $f_i^*({{\xi }},t) = {f_i}({{\xi }},t) - \frac{1}{\tau }\left[ {{f_i}({{\xi }},t) - f_i^{(0)}({{\xi }},t)} \right].$

在计算平面 ${{\xi }} \equiv \left( {\xi_1 ,\xi_2 } \right)$中，离散速度定义为 ${{{\tilde c}}_i} \equiv \left( {{{\tilde c}_{i,1}},{{\tilde c}_{i,2}}} \right)$，有

(8) ${\tilde c_{i,\alpha }} = { c_{i,\,\beta }}{{\partial {\xi _\alpha }} / {\partial {x_\beta }}};\quad\alpha ,\beta = 1,2.$

式中：c_i,1、c_i,2为c_i的分量.

式（5）的物理平面中的 ${{x}} - {{{c}}_i}\Delta t$，对应于计算平面的位置是通过对离散速度 ${{{\tilde c}}_i}$积分得到的，迁移公式转换为

(9) ${f_i}({{\xi }},t + \Delta t) = f_i^*({{\xi }} - \Delta {{{\xi }}_{{\rm{up}},i}},t).$

式中： $\Delta {{{\xi }}_{{\rm{up}},i}}$为对流距离，通过二步Runge-Kutta方法获得：

(10) ${\text{第1步：}}\Delta {{\xi }}_{{\rm{up}},i}^{\left( 1 \right)} = \frac{1}{2}\Delta t{{{\tilde c}}_i}\left( {{\xi }} \right),$

(11) ${\text{第2步：}}\Delta {{{\xi }}_{{\rm{up}},i}} \approx \Delta t{{{\tilde c}}_i}\left( {{{\xi }} - \Delta {{\xi }}_{{\rm{up}},i}^{\left( 1 \right)}} \right).$

在Runge-Kutta方法第2步中位于网格点之间位置的值是通过相邻网格点得到的. 为了保证LBM的二阶计算精度，在内部网格点处采用二阶迎风差值，与边界相邻的网格点处采用二阶中心插值. 二阶迎风差值表达式为

(12) ${\tilde { c}_i}\left( {{{\xi }} - \Delta {{\xi }}_{{\rm{up}},i}^{\left( 1 \right)}} \right) = \sum\limits_{m = 0}^2 {\sum\limits_{l = 0}^2 {{a_{i,\,m,2}}} } {a_{i,l,1}}{\tilde{{c}}_{i,m,l}}.$

式中： $\tilde { c}_{i,m,l} $为(m, l)网格点处的离散速度； ${a_{i,0,\alpha }}$、 ${a_{i,1,\alpha }}$、 ${a_{i,2,\alpha }}$为插值系数^[16]，定义为

(13) $ \left. {\begin{array}{*{20}{c}} {{a_{i,0,\alpha }} = {{\left( {\left| {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}} \right| - 1} \right)\left( {\left| {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}} \right| - 2} \right)} / {2,}}}\\ {{a_{i,1,\alpha }} = - \left| {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}} \right|\left( {\left| {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}} \right| - 2} \right),}\\ {{a_{i,2,\alpha }} = {{\left| {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}} \right|\left( {\left| {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}} \right| - 1} \right)} / {2.}}} \end{array}} \right\} $

式中： ${\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}} $为通过Runge-kutta方法第1步得到的对流距离.

二阶中心插值表达式为

(14) ${\tilde {{c}}_i}\left( {{{\xi }} - \Delta {{\xi }}_{{\rm{up}},i}^{\left( 1 \right)}} \right) = \sum\limits_{m = - 1}^1 {\sum\limits_{l = - 1}^1 {{a_{i,m,2}}} } {a_{i,l,1}}{\tilde { c}_{i,m,l}}.$

其中，插值系数^[16]定义为

(15) $\left. {\begin{array}{*{20}{l}} {{a_{i,-1,\alpha }} = \Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}\left( {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)} - 1} \right)/2,}\\ {{a_{i,0,\alpha }}{\rm{ = }}1 - {{\left( {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}} \right)}^2},}\\ {{a_{i,1,\alpha }}{\rm{ = }}\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)}\left( {\Delta \xi _{{\rm{up}},i,\alpha }^{\left( 1 \right)} + 1} \right)/2.} \end{array}} \right\}$

$f_i^*({\xi}-\Delta{\xi}_{{\rm{up}},i},t) $采用 $\tilde { c}_i({\xi} - \Delta{\xi}_{{\rm{up}},i}^{(1)}) $相同的插值方法进行求解.

为了保证GILBM数值计算的稳定性，时间步长 $\Delta t$通过CFL（Courant–Friedrichs–Lewy）条件获得：

(16) $\Delta t = \lambda \min \;\left| {{1 / {{{\left. {{{\tilde c}_{i,\alpha}}} \right|}_\xi }}}} \right|.$

式中： $\lambda $为CFL数， $\lambda \in \left( {{\rm{0,1}}{\rm{.0}}} \right]$； $\left| {{1 / {{{\left. {{{\tilde c}_{i,a}}} \right|}_\xi }}}} \right|$为计算平面中所有网格点的最小值.

为了检验GILBM方法在非均匀网格计算中的正确性，采用顶盖驱动方腔流动作为验证算例，分别对雷诺数Re = 100、1000的方腔进行数值模拟. 计算网格采用非均匀网格，如图1所示. 模型长度为 $100 \times 100$，网格数目为 $256 \times 256$，在4个边界进行加密，边界处第1层网格高度为0.1，方腔几何中心处网格高度为0.7. 顶盖驱动的速度 $U = 0.1c$（ $c = 1$），方腔长度为 $L$，雷诺数定义为 $Re = {{UL} / \upsilon }$，流场初始密度为1.0.

如图2（a）、（b）所示分别为沿着水平几何中心线的无量纲垂直速度分布曲线和沿着垂直几何中心线的无量纲水平速度分布曲线. 图中，u、v分别为水平、垂直速度. 可以看出，采用GILBM计算得到的速度分布曲线与Ghia等^[17]的计算结果较吻合，从定量上验证了GILBM方法应用在非均匀网格数值模拟的正确性.

对Re=10、20、40的圆柱流动进行模拟，在笛卡尔坐标系 $\left( {x,y} \right)$和极坐标系 $\left( {r,\theta } \right)$下的圆柱流动模型示意图如图3（a）所示. 图中， $\theta $为方位角，沿顺时针方向增加；均匀来流平行于 $x$轴，来流马赫数 $Ma = {{{U_\infty }} / {{c_{\rm{s}}}}}$（ ${U_\infty }$为远场来流速度， ${c_{\rm{s}}}$为音速）；雷诺数 $Re = {{{U_\infty }d} / \upsilon }$， $d$为圆柱直径， $d = 1.0$. 圆柱绕流计算网格如图3（b）所示，周向网格点为360，径向网格点为631，第1层网格高度为 $0.005$，圆柱壁面采用非平衡外推边界条件，远场采用无反射边界条件，时间步长为 $3.47 \times {10^{ - 3}}$.

如图4所示为不同雷诺数Re=10、20、40的圆柱表面压力系数分布曲线图. 图中， $C_{\rm{P}} = {\left( {p - {p_\infty }} \right)} /$ ${\left( {0.5\rho U_\infty ^2} \right)}$， ${{p_\infty }}$、 ${{U_\infty }}$分别为远场来流压力、速度. 可以看出，计算得到的表面压力系数分布与He等^[16]采用ISLBM计算的结果较吻合，并且圆柱上、下表面的压力系数呈现对称分布，最大压力出现在圆柱前缘的驻点处，验证了GILBM方法应用于非均匀贴体网格复杂流场计算的正确性.

采用GILBM方法在非均匀网格中计算二维高斯脉冲传播，本研究在高斯声源施加处进行网格加密，最小网格高度为0.00125，远离声源点处网格逐渐稀疏，最大的网格高度为0.00550，网格数目为 $401 \times 401$，计算域长度为 $1 \times 1$，计算时间步长 $\Delta t \approx 0.001\;25$，计算模型与计算网格如图5所示.

当 $t = 0$时，在 $\left( {{x_0} = 0.5,{y_0} = 0.5} \right)$位置处施加高斯脉冲声源^[18]：

(17) $\rho = {\rho _0} + \varepsilon \exp \;\left( { - \alpha \left( {{{\left( {x - 0.5} \right)}^2} + {{\left( {y - 0.5} \right)}^2}} \right)} \right).$

式中： ${\rho _0}$为流场初始密度，设置为1.0； $\varepsilon $为高斯脉冲扰动， $\varepsilon = 0.01$； $\alpha = (\ln\; 2)0.001\;6$；u=0、0.1， $v = 0$. 在采用GILBM方法计算时，为了将黏性对声波的影响降到最小，将式（3）中的 $\tau $设置为0.5，即为无黏状态.

如图6（a）、（c）所示分别为不同的水平来流速度下， $t = 0.4$时刻的密度云图. 可以看出，当水平来流速度 $u = 0$时，声波呈现出与声源施加点位置对称的结果，随着水平速度的增大，声波出现多普勒效应，左侧波向前传播的距离变小，右侧波向前传播的距离变大，声波整体向右平移. 如图6（b）、（d）所示分别为不同的水平速度条件下， $t = 0.4$时刻沿着 $y = 0.5$水平线的密度分布曲线，采用GILBM方法在非均匀网格下计算得到的密度分布曲线与文献[18]给出的解析解较吻合，验证了在非均匀网格中GILBM方法能较好地模拟高斯脉冲传播.

采用GILBM方法在非均匀网格上对周期性声源的传播算例进行数值模拟. 计算网格在如图5（b）所示的施加高斯声源的网格基础上放大100倍，即最小网格高度为0.125，最大网格高度为0.550，计算时间步长 $\Delta t \approx 0.125$. 周期性声源施加与文献[3]一致，即

(18) $\rho \left( {{{x}},t} \right) = {\rho _0} + {\rho _{\rm{s}}}\sin \;(2{\text{π}} {t / T}).$

式中：初始流场设置为静止流场； ${\rho _0}$为初始密度， ${\rho _0} = 1.0$； ${\rho _{\rm{s}}}$为周期性声源幅值， ${\rho _{\rm{s}}} = 0.01$； $T$为振荡周期， $T = 20$；取 $\tau = 0.6$.

如图7（a）所示为 $t = 75$时刻的扰动密度云图. 图中， $\Delta \rho = \rho - {\rho _0}$. 可以看出，在非均匀网格下，周期性声源以圆形波的形式向外传播，随着传播距离增加，密度扰动逐渐减小. 如图7（b）所示为 $t = 75$时沿着 $y = 50$的扰动密度曲线图. 可以看出，GILBM计算得到的扰动密度曲线与Viggen^[3]给出的解析解较吻合，验证了GILBM方法在非均匀网格下可以较好地模拟周期性声源的传播过程.

进一步对低雷诺数下圆柱绕流由于涡脱落产生的气动噪声进行数值模拟研究. 计算模型与网格如图3所示. 计算条件与Inoue等^[19]的计算条件一致，其中 $Re = 150$， $Ma = 0.2$.

如图8所示为应用GILBM方法计算得到的时均表面压力系数图. 与Inoue等^[19]采用二维非定常可压缩Navier-Stokes进行高精度直接数值模拟（direct numerical simulation，DNS）计算得到的结果和Lafitte等^[20]应用PowerFLOW软件基于LBM计算得到的结果进行对比，发现本研究计算得到的表面压力系数趋势与Inoue等^[19]的一致，在θ=0°和60°<θ≤180°处存在误差，在滞止点处采用GILBM计算得到的表面压力系数最大值 $C_{{\rm p}{\max }}$= 1.05，而Lafitte等^[20]计算得到的 $C_{{\rm p}{\max }}$=1.07，Inoue等^[19]计算得到的滞止点 $C_{\rm p}$=1.00，结果表明GILBM的误差比Lafitte等^[20]的误差小，优于Lafitte等^[20]所得结果.

如图9（a）所示为圆柱升力系数随时间的变化曲线. 图中， ${C_{\rm{L}}}$为升力系数， ${C_{\rm{L}}} = {{{F_{\rm{L}}}} / {\left( {{{{\rho _\infty }U_\infty ^2d} / 2}} \right)}}$； ${F_{\rm{L}}}$为圆柱升力. 如图9（b）所示为监测点r=75， $\theta = {90^\circ }$处的扰动声压 $\Delta p$随时间的变化曲线. 可以看出，升力系数与扰动声压的变化呈现周期振荡，周期均为 $T = 27.306$，通过升力系数曲线计算得到的斯特劳哈尔数（ $St$， $St = fd/{U_\infty }$，f为升力系数的振荡频率）为0.183，与Inoue等^[19]得到的 $St = 0.183$完全一致，并且接近于Williamson^[21]通过实验得到的 $St = 0.18$. 另外， ${C_{\rm{L}}}$、 $\Delta p$随时间变化的曲线与Inoue等^[19]的结果在频率和幅值上都较吻合，证明采用GILBM方法不仅可以准确计算出圆柱绕流旋涡脱落的情况，即模拟近场的声源情况，同时也可以准确模拟远场噪声的传播情况.

如图10所示为圆柱绕流气动噪声特性图. 如图10（a）所示为瞬时声压分布云图. 图中，白色为正声压，黑色为负声压. 可以看出，负声压的稀疏波与正声压的压缩波在圆柱周围交替产生并向上下游传播，随着传播距离的增大，声波的强度逐渐减弱，然而即使传播到远场，声波条纹依然较清晰，说明GILBM方法可以较好地模拟声波的传播过程. 如图10（b）所示为极坐标系下 $r = 75$处的声压波动均方根值，极坐标最外圈的声压均方根为 $1 \times {10^{ - 4}}$， $r = 75$处声压波动指向性呈现出明显的偶极子现象，由于来流速度造成的多普勒效应的影响，声压指向性向来流上游方向偏移.

（1）GILBM在非均匀网格下对顶盖驱动方腔流、圆柱绕流进行模拟，计算结果与参考值较吻合，验证了在非均匀网格下采用该方法模拟的正确性.

（2）GILBM可以在非均匀网格中较好地模拟气动噪声的传播过程，计算结果与解析解较吻合，验证了该方法模拟声传播问题的正确性.

（3）GILBM可以较好地模拟绕圆柱气动噪声的产生和传播过程，计算得到的斯特劳哈尔数、瞬时声压波动云图、声压波动曲线等均与参考值较吻合，瞬时声压分布云图和声压波动指向性图呈现出明显的偶极子指向性，验证了该方法模拟曲面物体气动噪声的可行性.

（4）本研究目前只开展了二维的数值模拟，后续可以将该方法应用于三维较复杂外形的气动噪声计算.