1.   引　言

声波在流场中传播时, 流场会与声波发生流声耦合作用, 使声波幅值和相位等产生明显变化^[1,2]. 涡旋是流体流动基本结构之一, 涡流对声波的散射问题是声波在复杂流场中传播的基本问题^[3-5], 在水声学、气动声学中具有广泛应用. 目前涡流对声波的散射问题的相关研究大多是定常涡流对声波的散射^[6-8], 定常涡流场并不随时间发生变化, 只具有空间分布, 而由于环境扰动, 随时间周期性变化的振荡涡流场在自然界和实际工程中广泛存在, 如水下运动目标激发的振荡伴流场^[9]、钝体结构物在绕流场中产生的振荡尾涡^[10,11]. 研究水下振荡涡流场声散射调制产生机理及相关特性, 对水下目标探测和流场声成像等领域具有重要意义.

流声耦合问题引起学者的广泛关注, 进行了大量相关研究, 建立了理论分析模型. Lighthill^[12]针对流声耦合问题, 采用Lighthill 声比拟方法发展了最早的理论分析模型, 该研究工作集中于声场的空间散射, Kraichnan^[13]基于Lighthill 声比拟方法推导了各向同性湍流的声散射表达式. Bogey等 ^[14,15]通过对欧拉方程进行线性化推导得到了线性欧拉方程(LEE), 在气动声学领域具有广泛应用, 但由于有关气体常数的存在, 无法直接用于水声学领域. Pierce ^[16,17]基于线性化理论, 在可压缩流体力学方程基础上推导了运动介质的波动方程, 能描述缓变非均匀运动流体介质对声传播的影响, 可以求解水下振荡涡流场中的声传播过程, 其进一步对运动介质的波动方程进行一阶近似, 通过引入势函数推导了单一变量的近似波动方程. 同时, 数值模拟计算是目前研究流声耦合问题的主要手段, Colonius等^[5]基于线性纳维-斯托克斯方程(LNS), 选取了多个基本涡旋结构, 采用直接数值模拟(DNS)开展了声散射研究, 并将结果与基于WKB 近似和 Born 近似方法的理论分析结果进行了比较, 为后来模型验证提供了数值参考. Karabasov等^[6]和Cheinet等^[18]分别基于时域有限差分(FDTD)方法和加权本质无震荡(WENO)方法对线性欧拉方程进行数值求解, 研究了定常涡旋对声场的散射作用. 对时空变化涡流与声场耦合作用进行数值模拟, 对于数值方法精度和稳定性具有一定要求, 间断伽辽金(DG)方法能够实现高阶数值格式, 数值精度高且有良好的收敛性和稳定性. Cockburn和Shu^[19,20]基于间断伽辽金方法和显式时间积分方法求解了非线性双曲型问题, 建立了龙格-库塔间断伽辽金方法, Lee和Kwon^[21]、Williamschen和Gabard^[22]基于间断伽辽金方法的高阶格式计算了背景流场对声场的影响.

尽管涡流的声散射问题已经研究了很长时间, 但目前鲜有关于水下振荡涡流场声散射调制的相关研究. 针对水下低频振荡涡流场声散射调制问题, 本文建立了理论分析模型与数值计算方法, 研究了声散射调制声场产生机理与时空频特性. 首先, 基于运动介质的波动方程, 通过引入势函数将波动方程分解为流声耦合项和非耦合项, 通过对流声耦合项的频域分析, 理论揭示了振荡涡流场的声散射调制机理; 其次, 采用间断伽辽金数值计算方法对水下低频振荡涡流场中声传播过程进行了数值模拟, 研究了不同入射声波频率、涡流场的振荡频率、马赫数和涡核尺度条件下振荡涡流场声散射调制声场的时空频特性.

2.   物理模型

声波在具有时空变化的涡流场中传播时, 声场空间分布和时间特性会发生改变, 发生流声耦合作用. 为了对流声耦合作用进行数学描述, 对流体力学方程组进行线性化近似后, 得到缓变非均匀运动介质的波动方程^[16,17,23]:

$$\left\{ \begin{gathered} \frac{1}{{{\rho _0}{c^2}}}\frac{{{\text{D}}{p^\prime }}}{{{\text{D}}t}} + \nabla \cdot {{\boldsymbol{u}}^\prime } + \frac{1}{{{\rho _0}{c^2}}}{{\boldsymbol{u}}^\prime } \cdot \nabla {p_0} + \frac{1}{{{\rho _0}{c^2}}}{p^\prime }\nabla \cdot {{\boldsymbol{u}}_0} = 0, \\ \frac{{{\text{D}}{{\boldsymbol{u}}^\prime }}}{{{\text{D}}t}} + \frac{1}{{{\rho _0}}}\nabla {p^\prime } + \left( {{{\boldsymbol{u}}^\prime } \cdot \nabla } \right){{\boldsymbol{u}}_0} - \frac{{{\rho ^\prime }}}{{\rho _0^2{c^2}}}\nabla {p_0} = 0, \\ \end{gathered} \right.$$

(1)

其中, ${\text{D}}/{\text{D}}t = \partial /\partial t + {{\boldsymbol{u}}_0} \cdot \nabla$, ${\rho }_{0}(r, t), {u}_{0}(r, t), {p}_{0}(r, t)$分别为不考虑声扰动时的流体密度、速度矢量、压强, ${\rho }'(r, t), {u}'(r, t), {p}'(r, t)$为对应的声场扰动量, 分别为密度逾量、振速矢量、声压. 方程组(1)能够理论分析和数值求解水下振荡流场与声场发生的耦合作用.

为了理论揭示振荡涡流场声散射调制机理, 将含有多个变量的耦合方程组(1)化为单一变量的近似波动方程. 对方程组(1)第一个式子求物质导数${\text{D}}/{\text{D}}t$, 对第二个式子求散度$\nabla \cdot$, 相减得到:

$$\begin{split} & \nabla \cdot \left( {\frac{1}{{{\rho _0}}}\nabla {p^\prime }} \right) - \frac{1}{{{\rho _0}{c^2}}}\frac{{{{\text{D}}^2}{p^\prime }}}{{{\text{D}}{t^2}}} - \frac{1}{{{\rho _0}{c^2}}}\frac{{\text{D}}}{{{\text{D}}t}}\left( {{p^\prime }\nabla \cdot {{\boldsymbol{u}}_0}} \right) \\ & + \nabla \cdot \left[ {\left( {{{\boldsymbol{u}}_0} \cdot \nabla } \right){{\boldsymbol{u}}^\prime }} \right] - {{\boldsymbol{u}}_0} \cdot \nabla \left( {\nabla \cdot {{\boldsymbol{u}}^\prime }} \right)\\ & + \nabla \cdot \left[ {\left( {{{\boldsymbol{u}}^\prime } \cdot \nabla } \right){{\boldsymbol{u}}_0}} \right] \\ =\; & \frac{1}{{{c^2}\rho _0^2}}\nabla \cdot \left( {{p^\prime }\nabla {p_0}} \right) + \frac{1}{{{\rho _0}{c^2}}}\frac{{\text{D}}}{{{\text{D}}t}}\left( {{{\boldsymbol{u}}^\prime } \cdot \nabla {p_0}} \right). \\[-15pt] \end{split}$$

(2)

基于声学中经典研究思想, 对方程(2)进行一阶近似, 得到^[17]:

$$\begin{split} & - \nabla \cdot \left[ {\frac{1}{{{\rho _0}}}\nabla {p^\prime }} \right] + \frac{{{{\text{D}}^2}}}{{{\text{D}}{t^2}}}\left( {\frac{1}{{{\rho _0}{c^2}}}{p^\prime }} \right) + 2( \nabla \cdot [ {\left( {{{\boldsymbol{u}}_0} \cdot \nabla } \right){{\boldsymbol{u}}^\prime }} ] \\ & - {{\boldsymbol{u}}_0} \cdot \nabla ( {\nabla \cdot {{\boldsymbol{u}}^\prime }} t) ) \approx 0.\\[-10pt] \end{split}$$

(3)

引入声压${p^\prime }({\boldsymbol{r}}, t)$、振速${{\boldsymbol{u}}^\prime }({\boldsymbol{r}}, t)$与速度势函数$\psi ({\boldsymbol{r}}, t)$的关系:

$${p^\prime }({\boldsymbol{r}},t) = - {\rho _0}\frac{{{\text{D}}\psi ({\boldsymbol{r}},t)}}{{{\text{D}}t}},$$

(4)

$${{\boldsymbol{u}}^\prime }({\boldsymbol{r}},t) = \nabla \psi ({\boldsymbol{r}},t) + O\left( {\frac{1}{R}} \right) + O\left( {\frac{1}{T}} \right) \approx \nabla \psi ({\boldsymbol{r}},t),$$

(5)

其中, $R, T$分别为介质非均匀变化的空间尺度和时间尺度, 对于缓变非均匀介质, $O\left( {1/R} \right)$, $O\left( {1/T} \right)$为高阶小量.

将势函数$\psi ({\boldsymbol{r}}, t)$与声压${p^\prime }({\boldsymbol{r}}, t)$、振速${{\boldsymbol{u}}^\prime }({\boldsymbol{r}}, t)$的关系代入方程(3)中, 经过化简后得到单一变量的波动方程, 展开后^[17]:

$$\begin{split} & {\nabla ^2}\psi - \frac{1}{{c_{}^2}}\frac{{{\partial ^2}\psi }}{{\partial {t^2}}} - \frac{1}{{c_{}^2}}\left[\frac{{\partial {{\boldsymbol{u}}_0}}}{{\partial t}}\nabla \psi + 2{{\boldsymbol{u}}_0} \cdot \frac{\partial }{{\partial t}}\nabla \psi \right.\\ & + {{\boldsymbol{u}}_0}\nabla \cdot {{\boldsymbol{u}}_0} \cdot \nabla \psi + {{\boldsymbol{u}}_0} \cdot {{\boldsymbol{u}}_0} \cdot {\nabla ^2}\psi ] = 0. \end{split}$$

(6)

为方便描述, 将方程(6)称为近似波动方程, 将方程(6)展开并分解为非耦合项$I({\boldsymbol{r}}, t)$和流声耦合项$R({\boldsymbol{r}}, t)$两部分:

$$I({\boldsymbol{r}},t) - \frac{1}{{c_{}^2}}R({\boldsymbol{r}},t) = 0,$$

(7)

$$I({\boldsymbol{r}},t) = {\nabla ^2}\psi ({\boldsymbol{r}},t) - \frac{1}{{c_{}^2}}\frac{{{\partial ^2}\psi ({\boldsymbol{r}},t)}}{{\partial {t^2}}},$$

(8)

$$\begin{split} R({\boldsymbol{r}},t) =& \frac{{\partial {{\boldsymbol{u}}_0}({\boldsymbol{r}},t)}}{{\partial t}} \cdot \nabla \psi ({\boldsymbol{r}},t) + 2{{\boldsymbol{u}}_0}({\boldsymbol{r}},t) \frac{\partial }{{\partial t}}\nabla \psi ({\boldsymbol{r}},t) \\ &+ {{\boldsymbol{u}}_0}({\boldsymbol{r}},t)\nabla \cdot {{\boldsymbol{u}}_0}({\boldsymbol{r}},t) \nabla \psi ({\boldsymbol{r}},t) \\ &+ {{\boldsymbol{u}}_0}({\boldsymbol{r}},t) \cdot {{\boldsymbol{u}}_0}({\boldsymbol{r}},t){\nabla ^2}\psi ({\boldsymbol{r}},t), \\[-12pt]\end{split}$$

(9)

其中, $I({\boldsymbol{r}}, t)$是与流场速度${{\boldsymbol{u}}_0}({\boldsymbol{r}}, t)$无关的项, 该项表示声波不会与流场发生耦合作用, 称为非耦合项; $R({\boldsymbol{r}}, t)$是与流场速度${{\boldsymbol{u}}_0}({\boldsymbol{r}}, t)$有关的项, 该项表示流场与声场发生耦合作用, 称为流声耦合项.

实际的振荡流场是一种复杂现象, 其时间变化特征并不一定具有严格的周期性, 但作为一种简化, 研究相关问题时可以近似为周期性变化^[9,24] , 涡流场速度函数${{\boldsymbol{u}}_0}({\boldsymbol{r}}, t)$为：

$${{\boldsymbol{u}}_0}({\boldsymbol{r}},t) = {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) + {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos ({\varOmega }t),$$

(10)

其中, ${\varOmega } = 2{\text{π }}F$为涡流场振荡角频率, ${{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}})$代表涡流场定常速度分量, ${{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos ({\varOmega }t)$代表涡流场振荡速度分量, 将(10)式代入(9)式, 此时流声耦合项$R({\boldsymbol{r}}, t)$为：

$$\begin{split} R({\boldsymbol{r}},t) =\;& 2{{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot \frac{\partial }{{\partial t}}\nabla \psi ({\boldsymbol{r}},t) + 2{{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t) \cdot \frac{\partial }{{\partial t}}\nabla \psi ({\boldsymbol{r}},t) - {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\varOmega \sin (\varOmega t) \cdot \nabla \psi ({\boldsymbol{r}},t) \\ &+ {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}})\nabla {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot \nabla \psi ({\boldsymbol{r}},t) + {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}})\nabla {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t) \cdot \nabla \psi ({\boldsymbol{r}},t) \\ &+ {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t)\nabla {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot \nabla \psi ({\boldsymbol{r}},t) + {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t)\nabla {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t) \cdot \nabla \psi ({\boldsymbol{r}},t) \\ &+ {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}){\nabla ^2}\psi ({\boldsymbol{r}},t) + {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t){\nabla ^2}\psi ({\boldsymbol{r}},t) \\ & + {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t) \cdot {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}){\nabla ^2}\psi ({\boldsymbol{r}},t) + {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t) \cdot {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\cos (\varOmega t){\nabla ^2}\psi ({\boldsymbol{r}},t). \end{split}$$

(11)

将流声耦合项从时域变换至频域进行分析, 用$\psi ({\boldsymbol{r}}, \omega )$表示$\psi ({\boldsymbol{r}}, t)$的频域形式, 根据Fourier变换性质和欧拉公式有:

$$\begin{split} & \mathcal{F}\left\{ {\psi ({\boldsymbol{r}},t)\cos (\varOmega t)} \right\} = \frac{1}{2}\psi ({\boldsymbol{r}},\omega \pm \varOmega ),~~ \mathcal{F}\left\{ {\psi ({\boldsymbol{r}},t)\sin (\varOmega t)} \right\} = \frac{1}{{2{\text{i}}}}[\psi ({\boldsymbol{r}},\omega - \varOmega ) - \psi ({\boldsymbol{r}},\omega + \varOmega )]. \end{split}$$

(12)

经过Fourier变换, $R({\boldsymbol{r}}, \omega)$为：

$$\begin{split} R({\boldsymbol{r}},\omega ) =\;& \mathcal{F}\left\{ {R({\boldsymbol{r}},t)} \right\}=\Big\{ 2{{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot \omega \nabla \psi ({\boldsymbol{r}},\omega ) + {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}})\nabla {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot \nabla \psi ({\boldsymbol{r}},\omega ) + {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}})\cdot {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}){\nabla ^2}\psi ({\boldsymbol{r}},\omega ) \\ & + \frac{1}{4}{{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\nabla {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}) \cdot \nabla \psi ({\boldsymbol{r}},\omega ) + \frac{1}{4}{{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}) \cdot {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\nabla \psi ({\boldsymbol{r}},\omega ) \Big\} \\ & +\Big\{ {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})(\omega \pm \varOmega )\nabla \psi ({\boldsymbol{r}},\omega \pm \varOmega ) \pm \frac{\varOmega }{{2i}}{{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\nabla \psi ({\boldsymbol{r}},\omega \pm \varOmega ) \\ &+ \frac{1}{2}{{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}})\nabla {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}) \cdot \nabla \psi ({\boldsymbol{r}},\omega \pm \varOmega ) + \frac{1}{2}{{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\nabla {{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot \nabla \psi ({\boldsymbol{r}},\omega \pm \varOmega ) \\ & + \frac{1}{2}{{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}){\nabla ^2}\psi ({\boldsymbol{r}},\omega \pm \varOmega ) + \frac{1}{2}{{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}){{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}}) \cdot {\nabla ^2}\psi ({\boldsymbol{r}},\omega \pm \varOmega ) \Big\} \\ & +\left\{ {\frac{1}{4}{{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\nabla {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}) \cdot \nabla \psi ({\boldsymbol{r}},\omega \pm 2\varOmega ) } + {\frac{1}{4}{{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}) \cdot {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})\nabla \psi ({\boldsymbol{r}},\omega \pm 2\varOmega )} \right\}. \\[-10pt] \end{split}$$

(13)

为了便于分析, 将频域流声耦合项$R({\boldsymbol{r}}, \omega )$写为：

$$R({\boldsymbol{r}},\omega ) = \left\{ {S({\boldsymbol{r}},\omega )} \right\} + \left\{ {{M_1}({\boldsymbol{r}},\omega \pm \varOmega )} \right\} + \left\{ {{M_2}({\boldsymbol{r}},\omega \pm 2\varOmega )} \right\},$$

(14)

$$S({\boldsymbol{r}},\omega ) = \sum \left\{ {\omega ^{{m_0}}}{\boldsymbol{u}}^{{i_0}}_{{\text{ste}}}({\boldsymbol{r}})[\nabla {\boldsymbol{u}}^{{j_0}}_{{\text{ste}}}({\boldsymbol{r}})]{\boldsymbol{u}}^{{k_0}}_{{\text{osc}}}({\boldsymbol{r}}) {{[\nabla {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})]}^{{l_0}}}{\nabla ^{{n_0}}}\psi ({\boldsymbol{r}},\omega ) \right\} ,$$

(15)

$${M_1}({\boldsymbol{r}},\omega \pm \varOmega ) = \sum \left\{ {{(\omega \pm \varOmega )}^{{m_1}}}{{\boldsymbol{u}}^{{i_1}}}_{{\text{ste}}}({\boldsymbol{r}})[\nabla {{\boldsymbol{u}}^{{j_1}}}_{{\text{ste}}}({\boldsymbol{r}})]{\boldsymbol{u}}^{{k_1}}_{{\text{osc}}}({\boldsymbol{r}}){{[\nabla {\boldsymbol{u}}_{{\text{osc}}}({\boldsymbol{r}})]}^{{l_1}}}{\nabla ^{{n_1}}}\psi ({\boldsymbol{r}},\omega \pm \varOmega ) \right\} ,$$

(16)

$${M_2}({\boldsymbol{r}},\omega \pm 2\varOmega ) = \sum \left\{ {\boldsymbol{u}}^{{k_2}}_{{\text{osc}}}({\boldsymbol{r}}){{[\nabla {{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})]}^{{l_2}}} \nabla \psi ({\boldsymbol{r}},\omega \pm 2\varOmega ) \right\} .$$

(17)

其中, ${m_{0, 1}}, {i_{0, 1}}, {j_{0, 1}}, {k_{0, 1, 2}}, {l_{0, 1, 2}} \in \left( {0, 1} \right)$, ${n_{0, 1}} \in (1, 2)$, 其特定组合能够使(14)式与忽略常系数后的(13)式相等.

下面基于频域流声耦合项(13)—(17)式对振荡涡流场与声场发生的耦合作用进行理论分析. 当涡流场处于振荡状态时, 即${{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}) \ne 0$, 振荡涡流场具有空间分布和时间变化特征, (14)式中$S({\boldsymbol{r}}, \omega )$表明涡流场能够使声场空间分布发生变化, 但不会改变声波频率, 产生频率为入射声波频率的散射声场. 为了便于区分和描述, 将$S({\boldsymbol{r}}, \omega )$称为散射项, 将不改变声场频率的耦合作用称为声散射作用, 散射声场强弱取决于$S({\boldsymbol{r}}, \omega )$大小. (14)式中${M_1}({\boldsymbol{r}}, \omega \pm \varOmega )$, ${M_2}({\boldsymbol{r}}, \omega \pm 2\varOmega )$表明振荡涡流场不仅能使声场空间分布发生改变, 振荡涡流场对声波散射可产生包含涡流场振荡频率双边带调制谐波的散射调制声场, 在频域中这种调制现象表现为入射声波频率两侧出现包含涡流场振荡频率分量的边频带$\omega \pm \varOmega$, $\omega \pm 2\varOmega$, 属于双边带调制, 其中, 边带包含一阶边频$\omega \pm \varOmega$, 二阶边频$\omega \pm 2\varOmega$. 将${M_1}({\boldsymbol{r}}, \omega \pm \varOmega ), {M_2}({\boldsymbol{r}}, \omega \pm 2\varOmega )$称为散射调制项, 把这种使声场频率发生改变的现象称为声散射调制作用, 将声场中产生的新频率称为调制频率, 散射调制声场强弱取决于散射调制项${M_1}({\boldsymbol{r}}, \omega \pm \varOmega )$, ${M_2}({\boldsymbol{r}}, \omega \pm 2\varOmega )$大小. 综上所述, 振荡涡流场对声波的散射会同时产生声散射作用和声散射调制作用, 对应的声场称为总散射声场. 当涡流场处于定常状态时, 即${{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}) = 0$, (14)式中散射调制项${M_1}({\boldsymbol{r}}, \omega \pm \varOmega )$, ${M_2}({\boldsymbol{r}}, \omega \pm 2\varOmega )$消失, 此时涡流场只会对声场产生散射作用.

3.   数值计算模型及验证

3.1   计算模型

基于间断伽辽金数值计算方法(DG)对运动介质的波动方程(1)进行求解, 数值模拟涡流场中的声传播过程, 时间步进采用四阶Runge-Kutta格式. 计算模型示意图如图1所示, 数值计算在直角坐标系$(x, y)$下进行, 原点定义在涡核中心. Gaussian涡是经典二维涡旋模型, 能够反应真实流场速度分布的部分特征^[3,6,25] , 基于对Gaussian涡模型的简单假设, 建立具有周期性变化的振荡涡流场模型, 其速度函数为：

$$\begin{split} &{u_x}(x,y,t) = \left[ {{\varGamma _{{\text{ste}}}} + {\varGamma _{{\text{osc}}}}\cos (2{\text{π }}Ft)} \right]\\ & \times\frac{{\left[ {1 - \exp \left( { - \alpha \left(\dfrac{{{x^2} + {y^2}}}{{{L^2}}}\right)} \right)} \right]\cos \left[ {\arctan \left( {\dfrac{y}{x}} \right)} \right]}}{{2{\pi}\sqrt {{x^2} + {y^2}} }}, \\ &{u_y}(x,y,t) = \left[ {{\varGamma _{{\text{ste}}}} + {\varGamma _{{\text{osc}}}}\cos (2{\text{π }}Ft)} \right]\\ &\times \frac{{\left[ {1 - \exp \left( { - \alpha \left(\dfrac{{{x^2} + {y^2}}}{{{L^2}}}\right)} \right)} \right]\sin \left[ {\arctan \left( {\dfrac{y}{x}} \right)} \right]}}{{2{\text{π }}\sqrt {{x^2} + {y^2}} }}, \end{split}$$

(18)

其中, $\alpha$为常数$1.256431$; $r = \sqrt {{x^2} + {y^2}}$为观测点到原点距离; ${\varGamma _{{\text{ste}}}} = 2.8{\text{π }}LcM{a_{{\text{ste}}}}$, $c$为介质声速, $M{a_{{\text{ste}}}}$表示涡流场定常速度分量的马赫数, 定义为涡流场最大定常速度与介质声速之比, $L$为涡核尺度, 涡流场定常速度分量在$r = L$处达到最大值; ${\varGamma _{{\text{osc}}}} = 2.8{\text{π }}LcM{a_{{\text{osc}}}}$, $M{a_{{\text{osc}}}}$表示涡流场振荡速度分量的马赫数, 定义为涡流场最大振荡速度与介质声速之比, 涡流场振荡速度分量在$r = L$处达到最大值, 当$M{a_{{\text{osc}}}} = 0$时, 涡流场处于定常状态, 涡流场速度不随时间变化; $F$为涡流场振荡频率, 振荡周期为${T_{{\text{osc}}}}$. 计算域是边长为$40 L \times 40 L$的矩形区域, 左边界为入口边界, 频率为$f$的平面入射声波沿$x$轴正向传播, 表达式为${p_{\text{i}}} = \sin (2{\text{π }}ft - kx)$, 声波反射主要发生在出口边界处, 而入口边界和上、下边界处基本无反射, 在出口边界处设置吸收层^[7,8], 以避免反射声波的影响. Gaussian涡为圆形涡流, 关于原点具有对称性, 不同时刻$y = 0$轴线上的涡流场速度随$x$坐标的变化曲线如图2所示, 涡流场速度分布随时间变化.

3.2   模型验证

为了考察本文数值方法的可行性和计算结果的准确性, 选取文献[5]标准算例进行对比, 文献计算了定常Gaussian涡流对声波的散射作用. 参数设置为: 马赫数$M{a_{{\text{ste}}}} = 0.125$, $M{a_{{\text{osc}}}} = 0$, 入射声波频率$f = c/(4 L)$, 涡核尺度$L = 1{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{m}}$, 观测曲线半径$r = 10 L$. 由${p_{{\text{rms}}}}/{p_{\text{i}}}$表示散射指向性, ${p_{{\text{rms}}}}$为散射有效声压:

$${p_{{\text{rms}}}} = \sqrt {\int p_{\text{s}}^2\dfrac{{\text{d}}t}{T} } = \sqrt {\int {{{({p_{\text{t}}} - {p_{\text{i}}})}^2}\dfrac{{\text{d}}t}{T}} } ,$$

(19)

其中, $p_{\text{s}}^{}$为散射声压, 由声波在涡流场中传播的总声压${p_{\text{t}}}$减去入射声压${p_{\text{i}}}$得到. 本文计算结果与文献结果比较结果如图3所示, 吻合程度很高, 证明了本文数值方法的有效性.

对于水声学问题, 通常涡流场马赫数较小, 且振荡频率远小于入射声波频率^[24] , 水下目标探测和海洋环流监测等问题多以低频声波为主^[9,26]. 参数设置为: 入射声波频率$f = 375{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{Hz}}$, 涡流场的振荡频率$F = 3.75{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{Hz}}$, 定义入射声波频率与涡流场振荡频率之比为频率比, 频率比$f/F = 100$; 涡核尺度$L = 1\;{\text{m}}$, 观测曲线半径$r = 10 L$; 涡流场马赫数为$M{a_{{\text{ste}}}} = 0.005$, $M{a_{{\text{osc}}}} = 0.0025$. 选择每波长网格数^[27]N分别为10, 15, 20, 30个和${\text{CFL}}$数^[19]分别为0.2, 0.1, 0.05, 0.025, 对振荡涡流场中的声传播过程进行数值求解, 散射指向性图如图4所示, 在$N \geqslant 15$, ${\text{CFL}} \leqslant {\text{0}}{\text{.1}}$时, 不同网格数和时间步长计算结果具有一致性, 无关性验证结果良好, 本文数值计算选择的每波长网格数$N$为15, ${\text{CFL}}$数为0.1.

4.   数值计算结果与分析

采用间断伽辽金数值计算方法(DG)对水下定常涡流场和低频振荡涡流场中的声传播过程进行数值求解. 研究了不同入射声波频率、涡流场的振荡频率、马赫数和涡核尺度条件下振荡涡流场声散射调制声场时空频特性, 并结合第2节中流声耦合项对其相关特性进行理论分析.

4.1   振荡涡流场声散射调制作用

振荡涡流场马赫数分别设置为: $M{a_{{\text{ste}}}} = 0.005$, $M{a_{{\text{osc}}}} = 0.0025$; $M{a_{{\text{ste}}}} = 0.005$, $M{a_{{\text{osc}}}} = 0.000625$; $M{a_{{\text{ste}}}} = 0.005$, $M{a_{{\text{osc}}}} = 0$. 其余参数与上文无关性验证部分的参数设置相同. 其中, $M{a_{{\text{osc}}}} \ne 0$时, 涡流场为振荡状态, $M{a_{{\text{osc}}}} = 0$时, 涡流场为定常状态. 散射声压周期性变化达到稳定后, 定常涡流场和振荡涡流场对声波散射得到的总散射声场指向性图如图5所示, 其中虚线示意主瓣方向. 不同时刻散射声压云图如图6—图8所示, 圆圈示意涡核位置. 散射指向性关于声波入射方向对称, 这与基于低马赫数大波长的理论分析结果一致^[28,29], 且有明显主瓣和旁瓣. 由于频率比远大于1, 定常涡流场和振荡涡流场对声波散射产生的散射声波频率相近, 因此其主瓣方位角一致, 散射指向性曲线形状基本相似. 当$M{a_{{\text{ste}}}} = 0.005$, $M{a_{{\text{osc}}}} = 0$时, 涡流场为定常状态, 散射声压幅值不随时间发生变化, 仅相位发生变化, 如图6所示. 当$M{a_{{\text{osc}}}} = 0.0025$, $M{a_{{\text{osc}}}} = 0.000625$时, 涡流场为振荡状态, 散射声压相位和幅值均随时间发生变化, 如图7和图8所示. 由于涡流场速度具有逆时针旋转方向, 为了方便讨论, 以$x$轴为分界线, 定义总散射声场上方区域为迎流区, 下方区域为顺流区^[25]. 在迎流区内, 涡流场速度水平分量方向与声波入射方向相反, 沿$x$轴负向. 在顺流区, 涡流场速度水平分量方向与声波入射方相同, 沿$x$轴正向. 可以看到迎流区与顺流区散射声压幅值呈现对称性, 而相位相反, 前向散射均远强于后向散射. 迎流区和顺流区主瓣方位角${\varTheta _ \pm }$位于约$\pm {29^ \circ }$处.

对振荡涡流场声散射调制时频特性进行分析, 选择观测曲线上迎流区主瓣方位角${\varTheta _ + }$处为观测点. 观测点散射声压时域图如图9所示, 为了便于观察, 通过Hilbert变换提取散射声压包络, 如时域图中红色线所示, 对散射声压进行频谱分析如图10所示. 当$M{a_{{\text{ste}}}} = 0.005$, $M{a_{{\text{osc}}}} = 0$时, 涡流场为定常状态, (13)式中${{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}}) = 0$, 此时流声耦合项(14)式中只有散射项$S({\boldsymbol{r}}, \omega )$, 涡流场只会对声场产生散射作用, 散射声压幅值不随时间发生变化, 如图9(a)所示, 散射声波频率仍为入射声波频率$f$. 当$M{a_{{\text{osc}}}} = 0.0025$, $M{a_{{\text{osc}}}} = 0.000625$时, 涡流场为振荡状态, 散射声压幅值随时间发生周期性变化, 如图9(b)和(c)所示. 振荡涡流场对声波的散射可产生散射调制声场, 散射声波频率变为$f$, $f \pm F$, $f \pm 2 F$. 其中, 一阶边频调制远强于二阶边频调制, 且由于二阶边频调制, 导致散射声压相邻两波包幅值略有差别. 关于马赫数变化对涡流场声散射调制作用的影响规律, 将在4.4节中进一步讨论.

4.2   不同入射声波频率条件下声散射调制特性

本节研究不同入射声波频率$f$条件下振荡涡流场声散射调制特性. 涡流场马赫数为$M{a_{{\text{ste}}}} = 0.005 {\text{ }}M{a_{{\text{osc}}}} = 0.0025$, 假设涡流场振荡频率$F = 3.75{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{Hz}}$不变, 改变入射声波频率$f$, 频率比$f/F$计算范围为25—400. 不同入射声波频率条件下散射指向性图如图11所示, 散射声压云图如图12所示. 随着入射声波频率增大, 不同角度处散射有效声压${p_{{\text{rms}}}}$明显增大, 散射指向性主瓣逐渐尖锐, 主瓣宽度减小, 主瓣方位角逐渐趋近于入射声波传播方向, 散射指向性仍关于声波入射方向对称.

在频域上对散射声场和散射调制声场强度变化规律进行量化分析, 将频率$f$在频域上对应的频谱幅值定义为${P_f}$, $f \pm F$在频域上对应的频谱幅值定义为${P_{f \pm F}}$, 以${P_f}$衡量散射声场强度, 以$({P_{f + F}} + {P_{f - F}})$衡量散射调制声场强度.

对不同入射声波频率$f$条件下声散射调制时频特性进行分析, 主要研究一阶边频. 观测点散射声压时域图如图13所示, 以最大频谱幅值进行归一化, 绘制归一化频域图如图14所示, 散射调制声场强度$({P_{f + F}} + {P_{f - F}})$和散射声场强度${P_f}$随入射声波频率变化规律如图15所示. 随着入射声波频率增大, 散射声波频率发生改变, 散射声压幅值明显增大; 散射声场强度${P_f}$增大, 观测点处散射声场强度与入射声波频率$f$具有接近正比关系; 散射调制声场强度$({P_{f + F}} + {P_{f - F}})$增大, 观测点处散射调制声场强度与入射声波频率$f$具有接近正比关系. 可以解释为(14)式中散射项$S({\boldsymbol{r}}, \omega )$和散射调制项${M_1}({\boldsymbol{r}}, \omega \pm \varOmega )$中包含有系数$f$、$f \pm F$, 随着入射声波频率$f$增大, 系数$f$, $f \pm F$增大, 散射项和散射调制项均增大, 使得振荡涡流场对声场的散射作用和散射调制作用均增强.

4.3   不同涡流场振荡频率条件下声散射调制特性

本小节研究不同涡流场振荡频率$F$条件下声散射调制特性. 涡流场马赫数为$M{a_{{\text{ste}}}} = 0.005 {\text{ }}M{a_{{\text{osc}}}} = 0.0025$, 假设入射声波频率$f = 375{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{Hz}}$不变, 改变涡流场振荡频率$F$, 频率比$f/F$计算范围为25—400. 不同涡流场振荡频率条件下散射指向性图如图16所示, 散射声压云图如图17所示, 在频率比$f/F \gg 1$条件下, 随着涡流场振荡频率的变化, 散射指向性和声场空间分布几乎没有变化.

不同涡流场振荡频率条件下观测点散射声压时域图及归一化频域图如图18、19所示, 频率比远大于1时, 随着涡流场振荡频率$F$增大, 散射声压频率发生明显变化, 但散射声压幅值几乎没有变化. 可以解释为(14)式中散射调制项${M_1}({\boldsymbol{r}}, \omega \pm \varOmega )$中包含有系数$f \pm F$, 由于频率比$f/F \gg 1$, 涡流场振荡频率$F$的改变对$f \pm F$几乎没有影响, 因此振荡频率变化对散射调制声场强度几乎没有影响.

4.4   不同涡流场马赫数条件下声散射调制特性

本小节研究不同涡流场马赫数条件下声散射调制特性. 研究低马赫数条件下, 水下结构物振荡尾涡、水下目标振荡伴流场马赫数范围内的声散射调制问题, 分别计算 $M{a_{{\text{osc}}}} = 0.0025$, $M{a_{{\text{ste}}}}$为$0.005—0.04$; $M{a_{{\text{ste}}}} = 0.005$, $M{a_{{\text{osc}}}}$为$0.0003125— 0.0025$. 不同涡流场马赫数条件下散射指向性图如图20所示, 散射声压云图如图21所示. 在低马赫数条件下, 随着马赫数增大, 不同角度处散射有效声压${p_{{\text{rms}}}}$增大, 由于涡流场速度具有逆时针旋转方向, 迎流区主瓣幅值略大于顺流区主瓣幅值, 但散射指向性主瓣方位角基本不变, 散射指向性曲线形状仍然相近.

不同涡流场马赫数条件下观测点散射声压时域图如图22和24所示, 绘制归一化频域图如图23和25所示, 散射调制声场强度$({P_{f + F}} + {P_{f - F}})$和散射声场强度${P_f}$随马赫数变化规律如图26所示. 散射声压频率取决于入射声波频率与涡流场振荡频率, 马赫数变化不会发生改变散射声波频率. 随着涡流场定常速度分量的马赫数$M{a_{{\text{ste}}}}$增大, (13)式中${{\boldsymbol{u}}_{{\text{ste}}}}({\boldsymbol{r}})$增大, 散射项$S({\boldsymbol{r}}, \omega )$增大, 使得散射声场强度${P_f}$增大, 而散射调制声场强度几乎没有变化, 观测点处散射声场强度与$M{a_{{\text{ste}}}}$具有接近正比关系. 同理, 随着涡流场振荡速度分量的马赫数$M{a_{{\text{osc}}}}$增大, 式(13)中${{\boldsymbol{u}}_{{\text{osc}}}}({\boldsymbol{r}})$增大, 散射调制项${M_1}({\boldsymbol{r}}, \omega \pm \varOmega )$增大, 使得散射调制声场强度$({P_{f + F}} + {P_{f - F}})$增大, 而散射声场强度几乎没有变化, 散射调制声场强度与$M{a_{{\text{osc}}}}$具有接近正比关系. 当$M{a_{{\text{ste}}}} \gg M{a_{{\text{osc}}}}$时, 散射声压时域图变为图9(a).

4.5   不同涡核尺度条件下声散射调制特性

本小节研究不同涡核尺度$L$条件下振荡涡流场声散射调制特性, 涡核尺度$L$计算范围为0.25—4 m. 不同涡核尺度$L$条件下散射声压指向性图如图27所示, 散射声压云图如图28所示. 随着涡核尺度$L$变化, 散射指向性变化规律与入射声波频率变化相同. 随着涡核尺度$L$增大, 指向性主瓣逐渐尖锐, 主瓣宽度减小, 主瓣方位角逐渐趋近于入射波传播方向.

不同涡核尺度$L$条件下观测点散射声压时域图及归一化频域图如图29和30所示, 调制声场强度$({P_{f + F}} + {P_{f - F}})$和散射声场强度${P_f}$随涡核尺度变化规律如图31所示. 随着涡核尺度$L$增大, 散射声压幅值明显增大, 散射声场强度${P_f}$增大, 观测点处散射声场强度${P_f}$与涡核尺度$L$具有接近正比关系. 散射调制声场强度$({P_{f + F}} + {P_{f - F}})$增大, 散射调制声场强度与涡核尺度$L$具有接近正比关系. 涡核尺度$L$变化对散射调制声场强度的影响等价于入射声波频率$f$的变化, 这与经典刚性球散射规律类似, 不同之处在于调制频率$f \pm F$会随着入射声波频率$f$变化, 而涡核尺度的变化不会影响调制频率.

5.   结　论

本文对水下低频振荡涡流场与声场耦合作用进行理论建模与数值模拟, 通过理论模型分析揭示了振荡涡流场的声散射调制机理, 采用间断伽辽金数值方法对水下低频振荡涡流场的声散射调制特性进行了计算分析. 探究了在低马赫数条件下不同入射声波频率、涡流场的振荡频率和涡核尺度对涡流场的声散射调制特性的影响规律, 并结合理论分析模型对其特性进行了解释, 得到以下结论:

1)振荡涡流场对声波的散射可产生包含涡流场振荡频率双边带调制谐波的散射调制声场, 散射声压幅值随时间发生周期性变化, 前向散射远强于后向散射, 一阶边频调制远强于二阶边频调制.

2)随着入射声波频率、涡流场振荡速度分量马赫数和涡核尺度的增加, 散射声压幅值增大, 散射调制声场强度增强, 在频率比远大于1条件下, 涡流场的振荡频率变化对散射调制声场强度影响较小.

3)在低马赫数和频率比远大于1条件下, 总散射声场空间分布具有对称性和明显的主瓣; 随着入射声波频率和涡核尺度增大, 指向性主瓣逐渐尖锐, 且主瓣方位角趋近于入射波传播方向; 涡流场的振荡频率和马赫数变化对散射指向性影响较小.

关于水下低频振荡涡流场声散射调制问题, 将在后面理论和实验中进一步深化研究.