0.   引言

圆柱壳—流场耦合振动及声辐射特性的研究工作很多^[1-4]，但是这些工作主要针对的是无限流域这类工况。在实际工程问题中，自由液面是常见的边界类型，因此需要对考虑自由液面影响的壳—液耦合振动及声辐射问题进行研究。

自由液面对壳—液耦合振动影响的研究工作较少，Ergin等^[5]基于实验和三维水弹性软件对有限浸没深度下圆柱壳振动特性进行分析，发现结构离自由液面越近，同阶次固有频率越大。Amabili^[6]对部分充液圆柱壳进行研究，提出了用扇形边界替代自由液面的近似方法。之后，Amabili^[7]将该方法拓展到处理部分浸没问题。随后，Ergin等^[8]利用边界积分法和镜像原理对部分充液（浸没）圆柱壳振动特性进行了研究，结果与实验数据符合良好。王鹏等^[9-10]基于波传播方法和虚源法，分析了浅水域圆柱壳的自由振动特性，但是该方法假设振型不随浸没深度变化并且忽略了虚源声波对结构振速的影响。Guo等^[11]建立了有限浸没深度下圆柱壳流固耦合物理模型，考虑了虚源速度势对结构振速的影响，并采用边界元方法研究了远场声辐射特性。王斌等^[12]从圆柱壳表面的均方振速和辐射声功率的角度，对半浸状态和全浸状态下圆柱壳在无限长线激励作用下的声振特性进行比较分析，指出了二者之间的差别与联系。刘佩等^[13]采用有限元软件ANSYS对有限深度浸没圆柱壳进行仿真，得到了与文献[5]类似的结论，并指出自由液面对圆柱壳自由振动的影响在浸没深度大于4倍半径时可以忽略不计。

半空间声学问题的解析研究工作较多，Huang^[14]研究了平面波入射下二维圆柱壳的散射声场，Hasheminejad等^[15-16]开展了有限空间谐振动二维圆柱的声场研究工作。白振国等^[17]采用镜像法建立了有限水深环境中二维圆柱壳的振动声辐射数学物理模型，初步计算了浅水对圆柱壳振动声辐射的影响规律及水深、潜深对声场分布和衰减特性的影响规律。Li等^[18]基于镜像原理进行了自由液面下有限潜深无限长圆柱壳结构的声辐射性能研究，基于稳相法，最终得到了有限浸没深度下圆柱壳结构远场辐射声压的计算表达式。但是，这些工作的研究对象皆是无限长圆柱壳，当圆柱壳长度为有限时，半空间声学问题的求解难度较大，目前尚缺乏解析研究。

圆柱壳振动时会向外辐射声波，由于自由液面的反射作用，反射声波会在结构表面发生刚性散射，而散射声触及自由液面又会形成回波，继而在结构表面产生多次散射和在界面产生多次反射（互散射），并最终形成稳态声场。由于辐射声和散射声的相似性，可将其统一表示。基于镜像原理，反射声均可认为由虚源发出，故将所有声波分为2类，即实源声和虚源声。若流体假设为不可压缩，虽然流体中不存在声波，但是也可通过镜像原理，将速度势设为实源速度势和虚源速度势来分析此类问题。

本文将基于此类思路分析计及自由液面影响的有限长圆柱壳振动响应及声辐射问题，实际上考虑了互散射效应。采用镜像原理来处理自由液面处声压为零的边界条件，再利用Graf加法定理对实源和虚源这2种坐标系进行转换得到速度势在流场的分布，然后再结合能量泛函变分方法得到壳—液耦合振动方程，进而可以求解其振动特性。之后，利用傅里叶变换构建波数域壳体表面连续条件，再利用傅里叶逆变换和稳相法求得其远场声压。

1.   理论分析

圆柱壳长度为L，厚度为h，中面半径为R，浸没深度为H，u, v, w分别表示轴向、周向和径向的中面位移，壳体材料的密度为ρ，弹性模量为E，泊松比为μ，流体密度为ρ_f。取圆柱壳左端面中心为坐标原点o，对应直角坐标(x, y, z)。实际分析中选择柱坐标系(x, r, φ)，其中x表示轴向，r表示径向，φ为周向角（与y轴的夹角），如图 1所示。

图  1  模型及坐标系

Figure  1.  Model and coordinate system

下载: 全尺寸图片 幻灯片

1.1   振动分析

为便于研究，本文取两端简支边界条件，因此位移场如式（1）所示。

式中：m为轴向半波数；n为周向波数；U_mn，V_mn，W_mn为三向位移幅值；k_m=mπ/L；ω为角频率；t为时间；i为虚数单位。

本文采用能量泛函变分的方法研究有限浸没深度圆柱壳的振动特性，故首先应得到各部分能量的表达式。壳体应变能（基于Love壳体理论^[19]）可表示为

式中：ε为应变向量；σ为应力向量；V为圆柱壳体积分域。

根据位移函数的正交性，积分后，式（2）可写成

式中：ξ_{m, n}=[U_mn, V_mn, W_mn]；刚度矩阵K_mn为三阶Hermite矩阵。

壳体动能可表示为

同理，根据位移函数正交性，积分后可表示为

式中，质量矩阵M_mn为三阶对角矩阵。

为求解流体做功，首先需要得到速度势函数的解析表达式。本文基于势流理论，将流体视为不可压缩、无旋、无粘性的理想流体，因此速度势函数ϕ(r, x, φ, t)满足柱坐标系下的Laplace方程：

对于水下圆柱壳，满足无穷远处速度势为零的条件：

由于自由液面的存在，可以借鉴镜像原理进行分析，认为速度势可由结构振动直接引起的实源速度势和自由液面反射的虚源速度势叠加组成。虚源坐标系(x′, r′, φ′)与实源坐标系关于自由液面对称，如图 2所示。

图  2  镜像原理示意图

Figure  2.  Schematic diagram of image method

下载: 全尺寸图片 幻灯片

设流域中任意一点为点P，其速度势函数可以表示为

式中：ϕ^r(r, x, φ, t)为实源流体速度势；ϕⁱ(r′, x′, φ′, t)为虚源流体速度势。

满足式（6）和式（7）的速度势函数有以下形式：

式中，K_n（）为第2类修正贝赛尔函数。

由于自由液面处速度势为零，故有

当P点位于自由液面上时，满足如下位置关系：

将式（9）和式（11）代入式（10），正交化处理后得到

即可以得到速度势函数的解析表达式：

根据Graf加法定理^[20]：

式中，I_a（）为第a阶第1类修正贝赛尔函数。

对于有限深度浸没，结构表面处半径r≈R < 2H，因此速度势解析表达式

因为系数a和n地位等价，交换顺序级数求和后上式改写为

根据圆柱壳外壁面处速度连续条件，有

将式（16）代入式（17）中，正交化处理后可以得到速度势幅值向量φ_n与位移幅值向量ζ_n的关系：

式中：Q为速度势幅值向量和位移幅值向量的关系矩阵；φ_n={ϕ_{m, -N}, ϕ_{m, -N+1}, …, ϕ_{m, N}}^T；ζ_n={W_{m, -N}, W_{m, -N+1}, …, W_{m, N}}^T，即可将速度势幅值向量用位移幅值向量表示。

由伯努利方程可以得到壁面处的流体动压力

流体做功为

点激励力做功为

式中，F₀为激励力幅值，激励力位置为(x₀, φ₀)处，激励力频率为f。

由上述各能量分量可得到能量泛函表达式：

根据变分原理，满足

由对幅值U_mn，V_mn的偏导为0可以得到其与幅值W_mn的线性关系，简写为如下所示的形式：

式中，a₁，b₁，c₁和a₂，b₂，c₂都是关于角频率ω、刚度矩阵K_mn和质量矩阵M_mn的系数，即U_mn，V_mn可由W_mn进行代换。

由于轴向函数在域内正交，当n的截断数为-N~N时，最终根据$\partial \prod /\partial {W_{mn}} = 0$可以得到轴向波数m取任意值时的方程：

式中：γ_n表示激励力展开后的幅值向量，γ_n=F₀sin(k_mx₀){exp(iNθ₀), exp[i(N-1)θ₀], …, exp(-iNθ₀)}^T；T_m为2N+1阶矩阵。

求解自由振动时，γ_n为零向量，因为ζ_n中元素不全为0，所以T_m必然不是满秩矩阵，即det(T_m)=0。由此可以求解出轴向波数m取任意值时各阶角频率ω，从而可以得到固有频率值。

求解受迫振动时，给定激励力频率f，可求出对应的ζ_n=T_m^-1γ_n，从而得到结构任意位置的振动响应。

1.2   远场声辐射分析

本文采用傅里叶变换方法，将轴向x变换到波数域，构建新的壁面连续条件（假设圆柱两端有两个半无限长声障柱），再进行逆变换即可求得声压表达式。

定义傅里叶变换及逆变换形式如下：

式中，k为波数，由此可将轴向坐标转化到波数域。

声压可划分为实源声压P_r和虚源声压P_i，根据分离变量法，实源声压可以表示为

式中，p_n为与周向角无关的声压物理量。

波数域流体声压满足Helmholtz方程：

式中，k_f=ω/c_f，为压缩波数，其中c_f为流体声速，角频率ω=2πf_h，f_h为谐振动频率。

由式（28）可以得到实源傅氏声压如下形式的解：

式中：H_n⁽¹⁾()为第n阶第1类Hankel函数；${k_r} = \sqrt {k_{\rm{f}}^2 - {k^2}}$，为径向波数；A_n(k)为对应的波数域幅值。

同理，虚源傅氏声压可以表示为

式中，B_n(k)为虚源声波数域幅值。

考虑到空气中的波阻抗远小于水中的波阻抗，声波由水中射向空气，自由液面处可以看成是绝对软的边界，自由液面处的声压可以近似认为是零，则自由液面处某点（r=r′，φ=π-φ′）的傅氏声压满足${\tilde P_{\rm{r}}} + {\tilde P_{\rm{i}}} = 0$，可以推导出

根据柱贝塞尔函数的Graf加法定理可以实现坐标迁移：

式中，J_n()为第n阶第1类贝塞尔函数。

同理，将结构径向位移w变换到波数域，有

式中，${\mathit{\boldsymbol{\tilde w}}_n}\left( k \right) = \frac{1}{{2{\rm{ \mathsf{ π} }}}}\sum\limits_{m = 1}^M {\frac{{{W_{mn}}{k_m}}}{{k_m^2 - {k^2}}}} \left[ {1 - {{\left( { - 1} \right)}^m}\exp \left( { - {\rm{i}}kL} \right)} \right]$。

根据壁面处连续条件，变换到波数域，有

将虚源傅氏声压迁移到实源坐标系下，因r≈R < 2H，有

交换积分顺序后，有

将式（36）代入式（34）后，级数做有限截断，均从-N取到N，可以得到波数域声压幅值向量${\mathit{\boldsymbol{\tilde A}}_n}\left( k \right)$与位移幅值向量${\mathit{\boldsymbol{\tilde w}}_n}\left( k \right)$的关系：

式中：${\mathit{\boldsymbol{\tilde A}}_n}\left( k \right) = {\left\{ {{{\tilde A}_{ - N}}\left( k \right), {{\tilde A}_{ - N + 1}}\left( k \right), \cdots, {{\tilde A}_N}\left( k \right)} \right\}^{\rm{T}}}$；T_ran为波数域声压幅值到位移幅值的迁移矩阵；${\mathit{\boldsymbol{\tilde w}}_n}\left( k \right) = {\left\{ {{{\tilde w}_{ - N}}\left( k \right), {{\tilde w}_{ - N + 1}}\left( k \right), \cdots, {{\tilde w}_N}\left( k \right)} \right\}^{\rm{T}}}$。

通过傅里叶逆变换，即可求出任意场点的声压。

对于远场声压，可以采用稳相法^[21]求解。限于篇幅，本文略去详细推导过程，直接给出球坐标系下声压表达式：

式中：A_a(k_fcosθ)可由式（37）求得；R₀为场点到原点的距离；θ为观测角（与轴向的夹角）。此处是将柱坐标系（x, r, φ）转换到球坐标系(R₀, θ, φ)下，有r=R₀sinθ，x=R₀cosθ。

2.   数值分析

2.1   振动分析

计算模型参数如下：壳长L=1.284 m，半径R=0.18 m，厚度h=0.003 m，壳体密度ρ=7 850 kg/m³，杨氏模量E=206 GPa，泊松比μ=0.3，流体密度ρ_f=1 025 kg/m³，流体声速c_f=1 500 m/s。点谐激励力位置坐标（R，L/2，0），激励力幅值F₀=1。

2.1.1   收敛性分析

为了说明方法收敛性，本文取H=0.2和2 m，分别计算在激励力频率为200和400 Hz时径向均方振速V_m随截断数M，N的变化规律，为便于分析收敛性，假设M=N。定义${V_m} = \frac{1}{{2s}}\int_s {{{\left| {{V_n}} \right|}^2}{\rm{d}}s}$，其中，s为表面积，V_n为径向速度。

从图 3可以看出，随着截断项数的增加，径向均方振速V_m值很快趋于稳定，当M=N=10时已经足够收敛，因此本文算例截断项数取值均为10。

图  3  径向均方振速收敛性分析

Figure  3.  The convergence curves of the quadratic velocity

下载: 全尺寸图片 幻灯片

2.1.2   准确性验证

为验证本文方法计算振动响应的准确性，取计算模型不变，激励力位置及幅值不变，激励力频率范围为2~500 Hz，扫频间隔为2 Hz，在壳体表面取一个测点，测点坐标x₁=L/4，φ₁=π/2。定义径向位移响应级RDL=20*lg（|w|/w₀），单位为dB，其中|w|为测点径向位移绝对值，w₀=10^-12 m。绘制浸没深度H=0.2 m和H=0.5 m时的径向位移响应级频谱曲线，并与有限元软件Nastran仿真（虚拟质量法）结果进行对比，如图 4所示。

图  4  本文方法和数值仿真径向位移级对比

Figure  4.  Comparison of the radial displacement levels between the present method and the numerical simulation

下载: 全尺寸图片 幻灯片

从图 4可以看出，本文方法计算结果和有限元软件Nastran仿真结果吻合良好，验证了本文方法的准确性。此外，本文方法的计算效率远高于有限元方法，以浸没深度为0.5 m工况下计算250个频率点振动响应为例，计算机配置均为4 GHz CPU与24 GB RAM，本文方法编程计算（Matlab2012）耗时约4 min，有限元软件Nastran（2 400个单元）耗时约2.2 h，本文方法的计算效率远高于数值仿真的计算效率。

2.1.3   浸没深度对受迫振动的影响

当圆柱壳靠近自由液面时，反射波会改变速度分布，因此需要研究浸没深度对圆柱壳速度分布的影响。计算模型及激励力位置和幅值不变，取浸没深度分别为半径的1.25，2.5，5，10和20倍，对比分析各工况下的径向均方振速，如图 5所示。

图  5  不同浸没深度下径向均方振速频谱

Figure  5.  The quadratic radial velocity versus frequency at different submerged depth

下载: 全尺寸图片 幻灯片

从图 5可以看出，越靠近自由液面，浸没深度对结构振动的影响越明显，并且随着浸没深度逐渐增大，振动响应趋于稳定（无限域工况）。此外，越靠近自由液面，曲线的峰值越往右移。换句话说，自由液面的存在会使共振频率增大，这与文献[5]的实验结论一致。

2.2   声辐射分析

2.2.1   声压指向性分析

得到壳体振动响应后，可以利用稳相法求解其远场辐射声压。取浸没深度为H=1，5和10 m这3组工况，分别计算激励力频率为200，300和400 Hz时的场点声压级。定义声压级SPL=20lg(|P_F|/P₀)，单位为dB，其中|P_F|为场点声压绝对值，P₀=10^-6 Pa。场点柱坐标为（1 000，1 000，φ），其中φ的取值范围为-π/2~π/2，取值间隔为π/180，流体声速c_f=1 500 m/s，从而可以绘制不同激励频率下各工况声压指向性图，并与边界元方法（BEM）计算结果进行对比，如图 6所示。

图  6  稳相法和边界元法声压级对比图

Figure  6.  Comparison of sound pressure level(SPL)between the stationary phase method and the boundary element method

下载: 全尺寸图片 幻灯片

从图 6可以发现，稳相法和边界元法计算得到的远场声压级指向性图吻合良好，说明稳相法准确可靠，而且相比于边界元方法，稳相法只是简单的赋值运算，计算效率显然更高。

从图 6还可以看出：远场声压指向性与频率和浸没深度相关，频率相同时，浸没深度越大，指向性曲线分瓣特性越明显；浸没深度相同时，频率越大，分瓣特性越明显。这是由于自由液面会带来类偶极子效应，若将结构认为是点源，则远场声压指向性与频率和浸没深度的乘积相关，乘积越大，分瓣越明显。但是因为结构并非极子，而且振动时表面速度并不均一，实际上指向性不能简单地用偶极子解释。

此外，值得注意的是，当浸没深度大于5倍半径以后，尽管振速会趋于稳定，但远场声辐射却不会趋于稳定。这是由于研究振动问题时，关注的是结构表面的速度势，虚源原点到结构表面的距离远大于实源原点到结构表面的距离，虚源速度势可以忽略；但是对于远场声辐射问题，实源原点和虚源原点到场点的距离几乎相同，虚源声压不可忽略。

2.2.2   声压波动性分析

为进一步研究远场声辐射随浸没深度变化的规律，选择观测点柱坐标（1 000，0，π/3），浸没深度H从0.2 m取到20 m，取值间隔0.01 m，激励力频率f分别取250和400 Hz，计算声场场点声压随浸没深度的变化规律，结果如图 7所示。

图  7  远场声压级随浸没深度变化曲线

Figure  7.  The curve of far-field sound pressure level with different immersion depth

下载: 全尺寸图片 幻灯片

从图 7可以看出，远场声压级随浸没深度变化呈现周期性波动，峰峰间距D随频率的增大而减小。但是在靠近自由液面时，声压级随浸没深度的变化比较剧烈，互散射效应比较强烈。

前文提及远场声压指向性规律可借鉴偶极子模型，为进一步探究水下圆柱壳受激振动远场声压级的波动规律，本文先以偶极子模型作为类比对象。

以两振幅相同而振动相位相反的脉动小球源为例，假设实源到场点的距离为r₁，虚源到场点的距离为r₂，两球相距2H，则场点声压可以表示为

式中：A为声压幅值；k_f为压缩波数。

假设平均半径r₀=(r₁+r₂)/2，则有r₁=r₀+Δr，r₂=r₀-Δr，其中Δr=(r₁-r₂)/2。由于观测点为远场，Δr远小于r₀，因此易于得到远场声压的近似表达式：

当Δr的增加量为S时，若满足k_f和S的乘积为π的整数倍，声压的绝对值保持不变，即声压呈现周期性波动，则波动一个周期时S=π/k_f。又因为声波波长λ=2π/k_f，所以S和λ满足如下数学关系：

若声场坐标系选择球坐标系（R₀, θ, φ），易于得到Δr与浸没深度H的几何关系：Δr=Hcosφsinθ，因此Δr的增量S与浸没深度H的增量D₀也满足相应的几何关系：S=D₀cosφsinθ，代入式（42）得

由式（43）可知，当浸没深度变化D₀=λ/（2cosφsinθ）时，声压级刚好变化一个周期，即偶极子模型声压级随浸没深度变化曲线的峰峰值为D₀。

在图 7（a）中激励力频率为250 Hz，因此声波波长λ=6 m，观测点位于柱坐标（1 000，0，π/3），转换到球坐标为（1 000，π/2，π/3），即θ=π/2，因此可用式（43）求出偶极子模型的峰峰值D₀=6 m，而由圆柱壳模型计算得到的远场声压级峰峰值D=5.85 m；同理，图 7（b）中圆柱壳模型的远场声压级峰峰值D=3.69 m，而偶极子模型的峰峰值D₀=3.75 m。从上述两例不难看出，偶极子模型计算的峰峰值与圆柱壳模型的计算值符合良好，这也从侧面说明本文方法是准确有效的。

为进一步说明该问题，选取激励力频率为500 Hz，以及球坐标系下4个观测点（1 000，π/6，π/6），（1 000，π/6，π/4），（1 000，π/4，π/6）和（1 000，π/4，π/4）。分别计算圆柱壳模型下远场声压峰峰值D和偶极子计算值D₀，如表 1所示。

表  1  不同观测点远场声压峰峰值及偶极子计算值对

Table  1.  peak-to-peak values of far field SPL and calculated values of dipole at different observation points

	φ=π/6			φ=π/4
	D/m	D0/m		D/m	D0/m
θ=π/6	3.40	3.46		4.15	4.24
θ=π/4	2.39	2.45		2.94	3.00

下载: 导出CSV 

| 显示表格

从表 1可以看出，对于不同观测点，圆柱壳模型下远场声压峰峰值D和偶极子计算值D₀符合良好，这也从场点声压波动性角度说明自由液面的存在会使得水下圆柱壳声辐射产生类偶极子效应。

3.   结论

本文提出了一种求解有限浸没深度下有限长圆柱壳振动特性及远场声的解析方法，相比于有限元或者边界元等数值方法，本文方法更加简便、计算效率更高，实现了外力激励下考虑自由液面影响的水下圆柱壳远场声辐射快速预报，对于半空间结构声振问题的研究具有一定的指导意义，并且可进一步拓展到研究刚性壁面或者组合边界问题。

具体结论如下：

1）自由液面的存在会增大同阶次共振频率，且越靠近自由液面，共振频率越大。

2）自由液面对声波的反射使得声场发生了菲涅耳干涉，远场声辐射的指向性和波动性都有类偶极子效应。