有限空间中结构的声辐射特性往往不可避免地受到边界影响，其中半空间声辐射问题的研究较为成熟。黎胜等利用半空间声学边界元法研究了自由液面和刚性壁面对脉动球声辐射的影响，结果表明刚性壁面对脉动球辐射声功率的影响远小于自由液面对脉动球辐射声功率的影响^[1-2]。邹元杰等建立了半无限流体域中结构流固耦合振动方程，探讨了自由液面和刚性壁面对结构的固有频率、振动响应和有关声学物理量的影响，结果表明边界影响明显存在并随结构离开边界而减弱^[3]。另外还有一些与半空间声学边界元法相关的文献，其重点在于流固耦合算法和半空间快速多极声学边界元法的研究^[4-9]。陈炉云等对垂直和水平壁面同时存在的四分之一和八分之一空间声辐射问题展开研究，建立了多虚源的边界元模型，拓展了半空间声学边界元法的应用范围^[10-12]。对类似平行波导的有限空间中结构声辐射问题的研究并不多见。邹元杰等建立浅水域声学边界元方程和相应的有限元/边界元流固耦合振动方程，研究平行边界对结构振动和声辐射的影响，其水底边界近似为刚性水底^[13]。Chen等利用波叠加法研究了半空间和平行波导中结构的声辐射问题^[14]；Zou等采用三维声弹性理论分析了浅海环境对结构声辐射的影响，其中水底边界的反射系数为掠射角的函数^[15-16]。白振国等采用解析法建立了浅水环境中二维圆柱壳的振动声辐射数学物理模型，研究了水深和浸深对圆柱壳声振特性的影响，结果表明圆柱壳离开水面一段小距离后其振动响应和表面声压与在无限水域中一致，边界对辐射声场的影响十分显著^[17]。在浅海环境中(水深一般为几十到几百米)，海底吸声作用不可忽略，本文首先对波导中格林函数进行修正，提出可考虑边界吸声作用的平行波导空间声学边界元法，然后研究吸声和刚性海底边界条件下浅海波导域中圆柱壳的声辐射特性及其规律，最后运用声学点源的辐射叠加原理解释辐射规律的产生机理。

1 海底吸声边界条件下的波导声学边界元法

浅海中结构声辐射和坐标系如图 1所示，设海面和海底构成介质均匀的平行波导空间，x轴垂直于海面。声波在海面和海底之间多次反射，任意场点声压可看作是实源和多个虚源在该点的声压叠加，求解其声场时采用的Helmholtz边界积分方程为^[18]

$ C\left( P \right)p\left( P \right) = \int_S {\left[{G\left( {P, Q} \right)\frac{{\partial p\left( Q \right)}}{{\partial {\boldsymbol{n}_s}}}-p\left( P \right)\frac{{\partial G\left( {P, Q} \right)}}{{\partial {\boldsymbol{n}_s}}}} \right]} \;{\rm{d}}S $

(1)

式中：P为声场中任意点，Q为结构边界上任意积分点，Q_Un和Q_Ln分别为Q在海面以上和海底以下的虚源，n为保证计算收敛的反射次数，也表示虚源序号，θ_n为Q_Ln与海底的掠射角，S表示结构边界，C(P)为P点声压系数，光滑连续边界上C(P)=0.5。p(P)和p(Q)分别表示P点和Q点声压，∂p(Q)/∂n_s=-iρ₀ωv_n，n_s为该边界点处法向单位向量，指向背离声域，在外部辐射问题中，指向结构内部。i表示单位虚数，v_n是单元法向速度，ω是圆频率，单元声压和法向速度通过插值获得。

平行波导空间中三维声学Helmholtz方程的基本解(格林函数)由实源项和虚源项叠加而成^[13]，可表示为

$ \begin{array}{*{20}{l}} {G\left( {P,Q} \right) = \frac{{{{\rm{e}}^{-{\rm{i}}{k_0}{R_0}}}}}{{4\pi {R_0}}}}\\ { + \sum\limits_{j = 1}^n {\left( {F_{Un}^mF_{Ln}^{j-m}\frac{{{{\rm{e}}^{{\rm{-i}}{k_0}{R_{Un}}}}}}{{4\pi {R_{Un}}}} + F_{Un}^{j - m}F_{Ln}^m\frac{{{{\rm{e}}^{{\rm{ - i}}{k_0}{R_{Ln}}}}}}{{4\pi {R_{Ln}}}}} \right)} } \end{array} $

(2)

式中：第一项代表实源，第二项代表虚源，R₀为P和实源Q的距离，R_Un为P和虚源Q_Un的距离，R_Ln为P和虚源Q_Ln的距离。F_Un和F_Ln分别为Q_Un和Q_Ln在海面和海底的反射系数，对于自由液面F_Un=-1，对于刚性海底F_Ln=1。k₀表示海水中的波数，m是与反射次数n相关的系数，m=j/2+[1+(-1)^j+1]/4。

由于实源和虚源在一条直线上，Q_Un和Q_Ln与Q的x坐标关系如下

$ \left\{ {\begin{array}{*{20}{c}} {{x_{Un}} = \left( {{h_U}-{h_L}-2{x_0}} \right)\frac{{1 + {{\left( {-1} \right)}^{n + 1}}}}{2} + nh + {x_0}}\\ {{x_{Ln}} = \left( {{h_U} - {h_L} - 2{x_0}} \right)\frac{{1 + {{\left( { - 1} \right)}^{n + 1}}}}{2} + nh + {x_0}} \end{array}} \right. $

(3)

式中：h_U和h_L分别为海面和海底到坐标原点o的距离，x₀为Q的x坐标，h=h_U+h_L表示浅海深度。

实际上海底并非刚性，其对声波的吸收是很明显的。由于基本解中的三项e^-ik₀R₀/(4πR₀)、e^-ik₀R_Un/(4πR_Un)和e^-ik₀R_Ln/4πR_Ln分别表示位于Q、Q_Un和Q_Ln的点源作用于场点P的声压，故式(1) 表示分布于结构表面不同强度点源及其镜像在场点P的声压叠加^[14-16]。在吸声边界条件下，点源格林函数中海底反射系数可表示为^[19]

$ {F_{Ln}} = F + B $

(4)

其中，

$ \begin{array}{l} F = \frac{{M{\rm{sin}}{\theta _n}-\sqrt {{N^2}-{\rm{co}}{{\rm{s}}^2}{\theta _n}} }}{{M{\rm{sin}}{\theta _n} + \sqrt {{N^2}-{\rm{co}}{{\rm{s}}^2}{\theta _n}} }}{\rm{, }}M = {\rho _1}/\rho, \\ N = {c_0}/{c_{10}}, B = \frac{1}{{2{\rm{i}}{k_0}{R_{Ln}}}}\left( {\frac{{{\partial ^2}F}}{{\partial \theta _n^2}} - \frac{{\partial F}}{{\partial {\theta _n}}}\tan {\theta _n}} \right)。 \end{array} $

海水密度和声速为ρ₀=1 025 kg/m³和c₀=1 500 m/s，海底沉积层密度ρ₁=1 400 kg/m³，考虑沉积层对声波吸收的声速^[20] c₁=1 530(1-iη)，其中η=0.021。

在结构表面划分边界单元，通过插值求得x₀并代入式(3) 得到x_Ln，则掠射角θ_n为

$ {\theta _n} = \arcsin \left( {\frac{{{x_0} + {x_P}}}{{{R_{Ln}}}}} \right) $

(5)

式中x_p是P的x坐标。

由式(4) 和(5) 求得反射系数F_Ln，再联合式(1)~(3) 便可求得吸声海底边界条件下的浅海平行波导空间声场。在F_Ln中，当k₀R_Ln足够大时可将B忽略，则F_Ln≈F。因为R_Ln随n的增大而增大，当h_L=50 m时，R_Ln≥50 m。令k₀=1 rad/m，R_Ln分别取50、500、1 000 m，其他参数同上。反射系数随掠射角的变化曲线如图 2所示，其中反射系数由F_Ln和F的绝对值表示。

图 2中F_Ln曲线随着R_Ln的增大而靠近F曲线，当R_Ln=500 m时，两条曲线已经十分接近，R_Ln=1 000 m时，两条曲线基本重合。同时可看出当掠射角超过30°后几条曲线几乎重合，所以在大角度掠射情况下k₀R_Ln对反射系数的影响较小。由于在40°~90°区间内，反射系数下降非常缓慢，为方便收敛性分析及声场计算，此区间内反射系数可近似看作常数。

2 收敛性分析和算例验证

在Matlab中编程分析格林函数G的收敛性，设频率f=200 Hz和800 Hz，h_U=10 m，h_L=50 m，Q点坐标为(0，0，0)，P点坐标为(0，50 m，0)，格林函数幅值|G|随反射次数n的变化曲线如图 3所示。

由图 3可看出吸声海底条件下的格林函数很快收敛，而刚性海底条件下的格林函数波动明显，收敛缓慢，这是因为两种海底边界的反射能力有很大差别，导致叠加效果不同。为保证收敛，在下文的计算中两种海底条件下的n分别取5和20。

在海底吸声边界条件下，以浅海平行波导空间中脉动球辐射解析解验证提出的声学边界元法的正确性。在脉动球表面划分若干个四边形单元，球心位于坐标原点o，f=200 Hz，场点分布为y轴上10~50 m的一条直线。此时F_Ln近似为常数，根据叠加原理^[1]可知场点声压解析解为

$ \begin{array}{l} p = \frac{{{\rm{i}}{\rho _0}\omega {a^2}{v_n}}}{{{r_0}\left( {1 + i{k_0}a} \right)}}{{\rm{e}}^{ - {\rm{i}}{k_0}({r_0} - a)}}\\ + \sum\limits_{j = 1}^n {\left[ {F_{Un}^mF_{Ln}^{j - m}\frac{{{\rm{i}}{\rho _0}\omega {a^2}{v_n}}}{{{r_{Un}}\left( {1 + {\rm{i}}{k_0}a} \right)}}{{\rm{e}}^{ - {\rm{i}}{k_0} ]\kern-0.15em] \left. {{r_{Un}} - a} \right)}} + } \right.} \\ \;\;\;\;\;\left. {F_{Un}^{j - m}F_{Ln}^m\frac{{{\rm{i}}{\rho _0}\omega {a^2}{v_n}}}{{{r_{Ln}}\left( {1 + {\rm{i}}{k_0}a} \right)}}{{\rm{e}}^{ - {\rm{i}}{k_0}\left( {{r_{Ln}} - a} \right)}}} \right] \end{array} $

(6)

式中：v_n=1×10^-6 m/s，脉动球半径a=1 m，r₀为球心到场点的距离，r_Un和r_Ln分别是球心在海面以上和海底以下的虚源到场点的距离。

图 4展示了吸声海底条件下声学边界元法与解析法计算的声压级(SPL)曲线，可看出两种方法的结果十分吻合，证明所提方法是正确的。

3 浅海中圆柱壳的声辐射特性

简支在半无限障板上的圆柱壳模型如图 5所示，圆柱壳长L_c=1.284 m，半径R_c=0.18 m，厚度为0.003 m，材料的杨氏模量为2.1×10¹¹ Pa，泊松比为0.3，密度为7 850 kg/m³。圆柱壳轴线与z轴重合，圆柱壳中心位于坐标系原点o，场点P在平面xy上，P到o的距离为r，γ为周向角。在圆柱壳下表面沿x轴施加频率为f，幅值为1 N的法向简谐激励力，激励力指向o。利用有限元法计算圆柱壳流固耦合模型的表面振速，当h=50 m时其均方振速在圆柱壳离开边界五倍半径距离后不再变化，这种情况与文献[17]研究结果一致。当h增大，边界影响将更小，故在h≥50 m情况下，圆柱壳离开边界五倍半径距离后可采用无限流域中圆柱壳的振速输入声学边界元Matlab程序中计算不同水深情况下圆柱壳的声辐射特性。

设f=200 Hz，h_U=10 m，h_L=50 m，不同半径下的声压级指向性曲线如图 6所示，其中SPLu、SPLa和SPLr分别表示忽略海底只考虑海面边界、海面+吸声海底和海面+刚性海底边界条件下的声压级。从图 6可看出两种海底边界条件下的声压级指向性曲线总体上存在明显差别，SPLa曲线比SPLr曲线更接近SPLu曲线，这是因为吸声海底对声波的反射能力较弱。当h_U和h_L不变，r越大，SPLa和SPLr曲线与SPLu曲线的差别越大，这是因为不同场点处的叠加效果不同。

令f=200 Hz，h_L=100 m，r=15 m，浸深h_U分别为50 m和200 m，且考虑静水压力作用，不同浸深下的声压级指向性曲线如图 7所示，其中SPLi为无限流域的声压级。从图 7可看出随着圆柱壳离开海面，SPLa曲线越来越接近SPLi曲线，这是因为圆柱壳离海底足够远，其反射作用可忽略，而海面的反射作用随浸深增加而减弱。再令h_U=h_L=50 m，声压级沿x轴的变化如图 8所示。从图 8看出，在圆柱壳下方三条曲线差别显著，SPLu曲线单调下降，而SPLa和SPLr曲线波动下降，并且SPLr曲线波动的剧烈程度远高于SPLa曲线；在圆柱壳上方三条曲线差别较小，波动情况相同。进一步发现，声压在任何情况下都表现出一种周期波动性，表明边界对声波的反射造成深度方向上出现驻波。波峰位置和波峰间距L并不因边界条件改变而改变，当f=200 Hz时波长λ=7.5 m，L=3.79 m，λ/L=1.98；当f=800 Hz时波长λ=1.875 m，L=0.94 m，λ/L=1.99，波长与波峰间距之比近似为一常数。

当h_U=10 m，h_L=50 m时，两种海底边界条件下的声压级分布云图如图 9所示。从图 9可见两种海底边界条件下的声压级分布有明显不同，吸声海底边界条件下的声压变化较小。在整个空间都有驻波分布，频率越高，驻波越密集，叠加声场中波峰波谷交错出现。

4 圆柱壳声辐射特性的产生机理

圆柱壳辐射声场中P点声压可以看作是无数个不同强度点源的辐射声压和反射声压在P点的叠加，式(2) 中格林函数第一项幅值为G₁，代表辐射声，第二项幅值为G₂，代表反射声。当f=200 Hz时，设h_U=10 m，h_L=50 m，Q点坐标为(0，0，0)，P点坐标为(0，y，0)；再设h_L=100 m，P点坐标为(0，15 m，0)，G₁和G₂随P点位置和浸深变化如图 10所示。从图 10看到G₂/ G₁随着P点离开圆柱壳而增大，随着圆柱壳离开海面而减小，其根本原因是反射声与辐射声的成分比例随场点到虚源距离的变化而变化。

在结构声振问题中，半无限域和无限域指在某些条件下可以忽略某些边界的声波反射作用的情况。在浅海波导域中，当P点离圆柱壳较近和圆柱壳距离海面或海底足够远时，有可能忽略边界反射作用，在此情况下布置测点得到圆柱壳的辐射声压等同于半无限域或无限域测试结果。为寻找满足声学半无限域或无限域条件的P点，可采用点源格林函数进行预测，不同参数下的格林函数幅值如图 11所示，其中G_a、G_u和G_i表示海面+吸声海底、忽略海底只考虑海面和无限流域中的格林函数幅值。从图 11可见格林函数幅值的匹配情况与图 6和图 7中声压级的匹配情况完全一致，这点进一步解释了图 6和图 7所展示现象的产生原因，说明采用点源格林函数预测浅海中圆柱壳声辐射特性是可行的。

声压在深度方向上的周期性波动现象同样源自点源的声压叠加性。图 1中实源和无限个虚源组成一条点源链，其中点源相位周期性变化。以构成一个周期的四个点源Q_U1、Q、Q_L1、Q_L2为例，设点源强度为A，则该周期的点源在P点的叠加声压为

$ \begin{array}{l} {\rm{ }}{p_0} =-\frac{A}{{{R_{U1}}}}{{\rm{e}}^{-{\rm{i}}{k_0}{R_{U1}}}} + \frac{A}{{{R_0}}}{{\rm{e}}^{-{\rm{i}}{k_0}{R_0}}}\\ + F\frac{A}{{{R_{L1}}}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{L1}}}} - F\frac{A}{{{R_{L2}}}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{L2}}}} \end{array} $

(7)

当P点下降ΔR=λ/2时，P点的叠加声压变为

$ \begin{array}{l} {p_1} =-\frac{A}{{{R_{U1}} + \Delta R}}{{\rm{e}}^{-{\rm{i}}{k_0}\left( {{R_{U1}} + \Delta R} \right)}} + \frac{A}{{{R_0} + \Delta R}}{{\rm{e}}^{-{\rm{i}}{k_0}({R_0} + \Delta R)}}\\ + F\frac{A}{{{R_{L1}} - \Delta R}}{{\rm{e}}^{{\rm{ - i}}{k_0}\left( {{R_{L1}} - \Delta R} \right)}} - F\frac{A}{{{R_{L2}} - \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}\left( {{R_{L2}} - \Delta R} \right)}}\\ = - \frac{A}{{R_{U1} + \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{U1}}}}{{\rm{e}}^{ - {\rm{i\pi }}}} + \frac{A}{{{R_0} + \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_0}}}{{\rm{e}}^{{\rm{ - i\pi }}}}\\ + F\frac{A}{{{R_{L1}} - \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{L1}}}}{{\rm{e}}^{{\rm{i\pi }}}} - F\frac{A}{{{R_{L2}} - \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{L2}}}}{{\rm{e}}^{{\rm{i\pi }}}} = \\ \frac{A}{{{R_{U1}} + \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{U1}}}} - \frac{A}{{{R_0} + \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_0}}}\\ - F\frac{A}{{{R_{L1}} - \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{L1}}}} + F\frac{A}{{{R_{L2}} - \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{L2}}}} \end{array} $

(8)

由式(7)、(8) 和指数函数的周期性质可知，p₀和p₁具有相同的周期性，当| p₀|为极大值时| p₁|也为极大值，距离和边界条件的变化只改变幅值大小而不改变周期性，这便是λ/L=2的原因。

当只考虑海面反射作用时，实源和虚源组成相位相反的偶极子，P点的叠加声压为

$ {p_2} =-\frac{A}{{{R_{U1}}}}{{\rm{e}}^{-{\rm{i}}{k_0}{R_{U1}}}} + \frac{A}{{{R_0}}}{{\rm{e}}^{-{\rm{i}}{k_0}{R_0}}} $

(9)

当P点在实源和海面之间下降ΔR=λ/2时，P点的叠加声压变为

$ \begin{array}{l} {p_3} =-\frac{A}{{{R_{U1}} + \Delta R}}{{\rm{e}}^{-{\rm{i}}{k_0}\left( {{R_{U1}} + \Delta R} \right)}} + \frac{A}{{{R_0}-\Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}({R_0} - \Delta R)}} = \\ - \frac{A}{{{R_{U1}} + \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{U1}}}}{{\rm{e}}^{{\rm{ - i\pi }}}} + \frac{A}{{{R_0} - \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_0}}}{{\rm{e}}^{{\rm{i\pi }}}} = \\ \frac{A}{{{R_{U1}} + \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_{U1}}}} - \frac{A}{{{R_0} - \Delta R}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_0}}} \end{array} $

(10)

同样，p₂和p₃具有相同的周期性，当时| p₂|为极大值时| p₃|也为极大值，在实源和海面之间λ/L=2。

当P点在实源下方下降时，λ/L=2规律的产生条件便不再满足。此时R_U1=R₀+2h_U恒成立，P点的叠加声压为：

$ \begin{array}{l} {p_4} = -\frac{A}{{{R_{U1}}}}{{\rm{e}}^{- {\rm{i}}{k_0}{R_{U1}}}} + \frac{A}{{{R_0}}}{{\rm{e}}^{- {\rm{i}}{k_0}{R_0}}} = \\ - \frac{A}{{{R_0} + 2{h_U}}}{{\rm{e}}^{ - {\rm{i}}{k_0}({R_0} + 2{h_U})}} + \frac{A}{{{R_0}}}{{\rm{e}}^{ - {\rm{i}}{k_0}{R_0}}} = \\ A{{\rm{e}}^{ - {\rm{i}}{k_0}({R_0} + {h_U})}}\left( {\frac{{{{\rm{e}}^{{\rm{i}}{k_0}{h_U}}}}}{{{R_0}}} - \frac{{{{\rm{e}}^{ - {\rm{i}}{k_0}{h_U}}}}}{{{R_0} + 2{h_U}}}} \right) = \\ A{{\rm{e}}^{ - {\rm{i}}{k_0}({R_0} + {h_U})}}\left[{\left( {\frac{1}{{{R_0}}}-\frac{1}{{{R_0} + 2{h_U}}}} \right)\cos \left( {{k_0}{h_U}} \right) + {\rm{i}}\left( {\frac{1}{{{R_0}}} + \frac{1}{{{R_0} + 2{h_U}}}} \right)\sin \left( {{k_0}{h_U}} \right)} \right] \end{array} $

(11)

可得到：

$ \begin{array}{l} \left| {{p_4}} \right| = A\sqrt {{{\left( {\frac{1}{{{R_0}}}-\frac{1}{{{R_0} + 2{h_U}}}} \right)}^2} + \frac{{4{{\sin }^2}\left( {{k_0}{h_U}} \right)}}{{{R_0}\left( {{R_0} + 2{h_U}} \right)}}} \le \\ \frac{A}{{{R_0}}} + \frac{A}{{{R_0} + 2{h_U}}} \end{array} $

(12)

由式(12) 可知|p₄|随着P点下降而减小，故图 8中SPLu在圆柱壳下方单调衰减。对应图 8的格林函数幅值如图 12所示(G_r表示海面+刚性海底边界条件下的格林函数幅值)，从图 12可看出格林函数幅值表现出与声压级相同的变化规律。以上现象是由流域边界对声波的反射作用产生的，而这点是早期圆柱壳声辐射问题研究文献并未关注的^[21]。

5 结论

1) 刚性和吸声海底边界条件下圆柱壳的辐射声压幅值有明显差异，吸声海底的声场叠加效果总体弱于刚性海底的声场叠加效果。

2) 可采用点源格林函数预测浅海波导域中圆柱壳的声辐射特性和寻找满足声学半无限域或无限域条件的场点，边界反射作用随着圆柱壳离开海面或海底而减弱，随着场点离开圆柱壳而增强。当圆柱壳与海面和海底距离满足一定条件时，其近场声辐射特性与半无限域或无限域中情况相同。

3) 海面和海底对声波的多次反射在浅海波导域中产生复杂的叠加声场，波峰波谷交错出现，频率越高，驻波越密集，其中深度方向上的驻波遵循波峰间距为波长一半的规律，此规律不受边界类型影响。