0 引　言

圆柱壳是水下航行器的重要组成结构，其振动声辐射特性的研究具有重要的工程意义。Wang^[1]研究了有限长圆柱壳的振动特性。内部含有肋骨^[2]、铺板^[3]、舱壁^[4]等结构的圆柱壳也得到了学者的关注。

Pan等^[5]采用波动法研究了加筋圆柱壳的自由振动特性，环肋采用平摊法处理，该方法局限于均匀分布的小尺寸环向加筋。庞福振等^[6]采用精细传递矩阵法评估水下加筋圆柱壳的辐射声压，应用叠加原理、增量存储以及改进非齐次项求解方法，提出了改进的精细传递矩阵法，并进一步分析了边界条件、流体介质以及圆柱壳厚度对辐射声压的影响。肖邵予^[7]基于Hamilton变分原理研究了两端用端板封口的加筋圆柱壳声辐射特性，以及圆柱壳内声场对其声振特性的影响。陈美霞等^[8]采用波动法研究了多舱段加筋圆柱壳的振动特性，加筋采用圆环板模型，舱壁采用圆板模型，给出了舱段截断时边界条件的选取原则。Wei等^[9]采用波动法研究了非均匀加筋圆柱壳的自由振动特性。

然而，在实际工程中，圆柱壳内部的结构是不规则的，解析方法不再适用。有限元法和边界元法可以用来分析复杂结构的振动声辐射问题，但是随着分析频率增大，网格数量急剧增多，计算效率会降低。

Maxit^[10]提出一种子结构法CAA（Circumferential Admittance Approach）来分析水下含轴对称内部结构（环肋、舱壁、端板等）圆柱壳的振动声辐射，圆柱壳在波数域中求解，内部结构利用轴对称有限元法处理。由于轴对称有限元法只能分析轴对称结构，这在一定程度上限制了CAA的应用范围。

本文结合解析法计算效率高和有限元法适用范围广的优点，将水下含内部结构圆柱壳分为水下圆柱壳和内部结构两部分。首先，采用轴对称边界元法^[11-12]计算圆柱壳表面声阻抗；其次，结合Flügge壳体理论^[13]以及三角函数，建立水下圆柱壳的数理模型；再次，采用有限元法得到内部结构的动力学方程，并对方程进行动力学缩聚，以缩减计算规模。最后，根据圆柱壳与内部结构接触点的位移协调和内力平衡条件，建立整个结构的动力学方程，求解方程得到壳体表面振动速度，再采用边界元法计算壳体辐射声压。

1 数理模型

整个模型可分解为圆柱壳和内部结构，由于采用有限元法对内部结构进行动力学建模，内部结构可以是梁、板以及其他复杂结构，为方便陈述，以内部基座为代表进行本文方法的介绍，如图1所示。柱坐标系 $ (r,\theta ,x) $ 中，两端简支圆柱壳长度为 $ L $ ，半径为 $ R $ ，厚度为 $ h $ ，密度为 $ \ \rho $ ，杨氏模量为 $ E $ ，泊松比为 $\ \mu $ ；水的密度为 $ {\ \rho _f} $ ，水中声速为 $ {c_f} $ ；基座的长、宽、高分别为 $ {L_b} $ ， $ {W_b} $ ， $ {H_b} $ 。

1.1 水下圆柱壳建模分析

基于Flügge壳体理论，水下圆柱壳运动微分方程可写为：

$ \begin{split}\dfrac{{Eh}}{{1 - {\mu ^2}}}\left[ {\boldsymbol{L}} \right]\left\{ {\begin{array}{*{20}{c}} u \\ v \\ w \end{array}} \right\} =& \left\{ {\begin{array}{*{20}{c}} {{F_x} + \dfrac{1}{{2R}}\dfrac{{\partial {M_r}}}{{\partial \theta }}} \\ {{F_\theta } + \dfrac{1}{R}{M_x} - \dfrac{1}{2}\dfrac{{\partial {M_r}}}{{\partial x}}} \\ {{F_r} - \dfrac{1}{R}\dfrac{{\partial {M_x}}}{{\partial \theta }} + \dfrac{{\partial {M_\theta }}}{{\partial x}}} \end{array}} \right\}\times\\ &\delta (x - {x_0})\delta (\theta - {\theta _0}) - \left\{ {\begin{array}{*{20}{c}} 0 \\ 0 \\ p \end{array}} \right\} \text{。}\end{split}$

(1)

式中， $ {\boldsymbol{L}} $ 是Flügge壳体微分算子^[13]； $ u $ ， $ v $ 和 $ w $ 分别是壳体轴向，周向和径向位移； $ {F_x} $ ， $ {F_\theta } $ ， $ {F_r} $ ， $ {M_x} $ ， $ {M_\theta } $ ， $ {M_r} $ 分别是作用在壳体 $ ({x_0},{\theta _0}) $ 位置的轴向力，周向力，径向力，轴向弯矩，周向弯矩，径向弯矩。

对于两端简支圆柱壳， $ u $ ， $ v $ ， $ w $ 可以表示为：

$ \left\{ \begin{gathered} u = \sum\limits_{m = 1}^M {\sum\limits_{n = - N}^N {{U_{mn}}\cos ({k_m}x)\exp (jn\theta )\exp (j\omega t)} }\text{，} \hfill \\ v = \sum\limits_{m = 1}^M {\sum\limits_{n = - N}^N {{V_{mn}}\sin ({k_m}x)\exp (jn\theta )\exp (j\omega t)} }\text{，} \hfill \\ w = \sum\limits_{m = 1}^M {\sum\limits_{n = - N}^N {{W_{mn}}\sin ({k_m}x)\exp (jn\theta )\exp (j\omega t)} }\text{。} \hfill \\ \end{gathered} \right. $

(2)

式中： $ {k_m} = {{m\text{π} } \mathord{\left/ {\vphantom {{m\pi } L}} \right. } L} $ 为壳体轴向波数； $ \omega $ 为角频率； $ j = \sqrt { - 1} $ 。 $ M $ 和 $ N $ 为模态截断数。

对一封闭结构， $ Q $ 点位于结构表面， $ P $ 是结构外部声场中某一场点的坐标，基于直接边界元积分方程， $ P $ 点的声压可表示为^[12]：

$ c(P)p(P) = \int\limits_S \left[ {p(Q)\frac{{\partial G(P,Q)}}{{\partial {n_q}}} + j{\rho _0}\omega \bar v(Q)G(P,Q)} \right]{\text{d}}{S_q}\text{。} $

(3)

式中： $ p(Q) $ 是结构表面 $ Q $ 点的声压； $ \bar v(Q) $ 是结构表面 $ Q $ 点的法向速度。当 $ P $ 点位于结构外部时， $ c(P) = 1 $ ；当 $ P $ 点位于结构表面时， $ c(P) = 1 $ 可通过下式计算。

$ c(P) = 1 + \int\limits_S {\frac{\partial }{{\partial {n_q}}}\left[ {\frac{1}{{4\text{π} \bar R(P,Q)}}} \right]{\text{d}}{S_q}}\text{。} $

(4)

对于轴对称封闭结构（见图2），式（3）中的物理量可展开为Fourier级数：

$ \begin{split} & p(P) = \sum\limits_{n = 0}^\infty {\left[ {p_n^{c * }\cos (n{\theta _p}) + p_n^{s * }\sin (n{\theta _p})} \right]} \text{，} \\ & p(Q) = \sum\limits_{n = 0}^\infty {\left[ {p_n^c\cos (n{\theta _q}) + p_n^s\sin (n{\theta _q})} \right]} \text{，} \\ & \bar v = \sum\limits_{n = 0}^\infty {\left[ {\bar v_n^c\cos (n{\theta _q}) + \bar v_n^s\sin (n{\theta _q})} \right]} \text{，} \end{split} $

(5)

式中， $ * $ 用于区分结构表面的点 $ Q $ 和结构外的点 $ P $ 。

经过数学推导^[11-12]，式（3）可转化为柱坐标系中的线积分方程：

$ \begin{split} & c(P)p_n^{c * }(P) = \int_\Gamma {\left[ {p_n^cH_n^\prime + j{\rho _0}\omega \bar v_n^c{H_n}} \right]{r_q}{\text{d}}{\varGamma _q}} \text{，} \\ & c(P)p_n^{s * }(P) = \int_\Gamma {\left[ {p_n^sH_n^\prime + j{\rho _0}\omega \bar v_n^s{H_n}} \right]{r_q}{\text{d}}{\varGamma _q}} \text{。} \end{split} $

(6)

其中：

$ {H_n} = \int_0^{2\text{π} } {\frac{{{e^{ - jk\bar R}}}}{{4\text{π} \bar R}}\cos (n\theta ){\text{d}}\theta }\text{。} $

(7)

$ {H_n} $ 可通过高斯积分计算。

针对式（6）采用高斯求积公式，将母线 $ \Gamma $ 离散为一系列线单元。为保证计算精度，本文采用二次单元，每个单元有3个节点，单元几何特征可以通过对节点插值进行描述：

$ r = \sum\limits_{i = 1}^3 {{r_i}{N_i}(\xi )} ,\;\;z = \sum\limits_{i = 1}^3 {{z_i}{N_i}(\xi )} \text{。}$

(8)

式中： ${N_1}{\text{ = }}\xi (\xi - 1)/2,\;{N_2} = (1 - \xi )(1 + \xi ),\;{N_3} = \xi (\xi + 1)/2$ 为形函数， $ - 1 \leqslant \xi \leqslant 1 $ 。

母线 $ \Gamma $ 上的微元 $ {\text{d}}\Gamma $ 可表示为：

$\begin{split} {\text{d}}\Gamma =& \sqrt {{{{\text{(d}}r)}^2} + {{({\text{d}}z)}^2}} = \sqrt {{{\left(\frac{{{\text{d}}r}}{{{\text{d}}\xi }}\right)}^2} + {{\left(\frac{{{\text{d}}z}}{{{\text{d}}\xi }}\right)}^2}} {\text{d}}\xi = \\ &\sqrt {{{\left(\sum\limits_{i = 1}^3 {{r_i}\frac{{{\text{d}}{N_i}}}{{{\text{d}}\xi }}} \right)}^2} + {{\left(\sum\limits_{i = 1}^3 {{z_i}\frac{{{\text{d}}{N_i}}}{{{\text{d}}\xi }}} \right)}^2}} {\text{d}}\xi = J{\text{d}}\xi\text{，}\end{split} $

(9)

采用节点插值，单元内的声压和法向速度可以表示为：

$ \begin{split} & p_n^c = \sum\limits_{i = 1}^3 {p_{ni}^c{N_i}(\xi )} ,\;\;\;\bar v_n^c = \sum\limits_{i = 1}^3 {\bar v_{ni}^c{N_i}(\xi )} \text{，} \\ & p_n^s = \sum\limits_{i = 1}^3 {p_{ni}^s{N_i}(\xi )} ,\;\;\;\bar v_n^s = \sum\limits_{i = 1}^3 {\bar v_{ni}^s{N_i}(\xi )} \text{。} \end{split} $

(10)

将式（8）～式（10）代入式（6），可将母线 $ \Gamma $ 上的线积分写一系列单元积分之和，即为

$ \begin{split} c(P)p_n^{c * }(P) =& \sum\limits_{e = 1}^{{N_e}} \left[ \int_{ - 1}^1 {p_{en}^cH_{en}^\prime {r_{eq}}{J_e}{\text{d}}\xi } + \right.\\ &\left.j{\rho _0}\omega \int_{ - 1}^1 {\bar v_{en}^c{H_{en}}{r_{eq}}{J_e}{\text{d}}\xi } \right]\text{，} \\ c(P)p_n^{s * }(P) = &\sum\limits_{e = 1}^{{N_e}} \left[ \int_{ - 1}^1 {p_{en}^sH_{en}^\prime {r_{eq}}{J_e}{\text{d}}\xi } +\right.\\ &\left.j{\rho _0}\omega \int_{ - 1}^1 {\bar v_{en}^s{H_{en}}{r_{eq}}{J_e}{\text{d}}\xi } \right]\text{，} \end{split} $

(11)

式中， $ {N_e} $ 为高斯单元的数量。单元内的积分采用高斯求积公式，式（11）可表示为：

$ \begin{split} c(P)p_n^{c * }(P) =& \sum\limits_{e = 1}^{{N_e}} \left[ \sum\limits_{g = 1}^{{N_g}} {{\kappa _g}p_{en}^c({\zeta _g})H_{en}^\prime ({\zeta _g}){r_{eq}}({\zeta _g}){J_e}{\text{(}}{\zeta _g}{\text{)}}} + \right.\\ &\left.j{\rho _0}\omega \sum\limits_{g = 1}^{{N_g}} {{\kappa _g}\bar v_{en}^c({\zeta _g}){H_{en}}({\zeta _g}){r_{eq}}({\zeta _g}){J_e}{\text{(}}{\zeta _g}{\text{)}}} \right] \text{，} \\ c(P)p_n^{s * }(P) =& \sum\limits_{e = 1}^{{N_e}} \left[ \sum\limits_{g = 1}^{{N_g}} {{\kappa _g}p_{en}^s({\zeta _g})H_{en}^\prime ({\zeta _g}){r_{eq}}({\zeta _g}){J_e}{\text{(}}{\zeta _g}{\text{)}}} +\right.\\ &\left.j{\rho _0}\omega \sum\limits_{g = 1}^{{N_g}} {{\kappa _g}\bar v_{en}^s({\zeta _g}){H_{en}}({\zeta _g}){r_{eq}}({\zeta _g}){J_e}{\text{(}}{\zeta _g}{\text{)}}} \right] \text{。} \end{split} $

(12)

其中： $ {N_g} $ 为高斯点的数量； $ {\zeta _g} $ 为高斯点坐标； $ {\kappa _g} $ 为高斯点对应的权系数。

采用“配点法”，将式（12）中的 $ P $ 点依次置于单元节点上，构造方程组，建立节点声压与节点速度的关联矩阵为：

$ \begin{split} &{\boldsymbol{{\rm A}}}_n^c \cdot {\boldsymbol{p}}_n^c = {\boldsymbol{B}}_n^c \cdot {\boldsymbol{\bar v}}_n^c \text{，} \\ &{\boldsymbol{{\rm A}}}_n^s \cdot {\boldsymbol{p}}_n^s = {\boldsymbol{B}}_n^s \cdot {\boldsymbol{\bar v}}_n^s \text{。} \end{split} $

(13)

式中，若单元 $ {\boldsymbol{\bar v}}_n^c $ 和 $ {\boldsymbol{\bar v}}_n^s $ 已知，可计算 $ {\boldsymbol{p}}_n^c $ 和 $ {\boldsymbol{p}}_n^s $ 。

根据式（2）和式（13），可计算圆柱壳表面的声阻抗，即圆柱壳表面声压与法向速度的关系，可表示为：

$ \begin{split}p =& C\bar v = j\omega Cw = \\ &j\omega C\sum\limits_{m = 1}^M {\sum\limits_{n = - N}^N {{W_{mn}}\sin ({k_m}x)\exp (jn\theta )\exp (j\omega t)} }\text{，} \end{split}$

(14)

将式（2）和式（14）代入式（1），可得

$ \left[ {{\boldsymbol{K}} + j\omega {\boldsymbol{C}}} \right]{\left\{ {{U_{mn}}\;\;{V_{mn}}\;\;{W_{mn}}} \right\}^{\rm T}} = {\boldsymbol{F}} \text{。}$

(15)

求解上式，可得 $ {U_{mn}} $ ， $ {V_{mn}} $ ， $ {W_{mn}} $ ，将其代入式（2），可得圆柱壳与内部结构接触节点的位移响应函数为：

$ \begin{split}{{\boldsymbol{U}}_0} = &\sum\limits_{m = 1}^M \sum\limits_{n = - N}^N {\left[ {{{\boldsymbol{R}}_{mn}}} \right]{{\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{mn}}}&{{{\boldsymbol{V}}_{mn}}}&{{{\boldsymbol{W}}_{mn}}} \end{array}} \right\}}^{\rm T}}} = \\ &\sum\limits_{m = 1}^M {\sum\limits_{n = - N}^N {\left[ {{{\boldsymbol{R}}_{mn}}} \right]{{\left[ {{\boldsymbol{K}} + j\omega {\boldsymbol{C}}} \right]}^{ - 1}}} {\boldsymbol{F}}} =\\ &\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{D}}_1}}&{{{\boldsymbol{D}}_2}} \end{array}} \right]{\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_1}}&{{{\boldsymbol{F}}_2}} \end{array}} \right\}^{\rm T}}\text{。}\end{split} $

(16)

其中： $ {{\boldsymbol{R}}_{mn}} $ 是表示圆柱壳位移响应与位移函数系数关系的矩阵； $ {{\boldsymbol{F}}_1} $ 是圆柱壳上的外载荷； $ {{\boldsymbol{F}}_2} $ 表示圆柱壳与内部接触点的耦合载荷。

1.2 内部结构建模分析

内部结构采用有限元法进行分析，材料损耗因子以复杨氏模量的形式考虑，内部结构的动力学方程可写为：

$ \left( {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{{ K}}}_{mm}}}&{{{\boldsymbol{{ K}}}_{ms}}} \\ {{{\boldsymbol{{ K}}}_{sm}}}&{{{\boldsymbol{{ K}}}_{ss}}} \end{array}} \right] - {\omega ^2}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{{ M}}}_{mm}}}&{{{\boldsymbol{{ M}}}_{ms}}} \\ {{{\boldsymbol{{ M}}}_{sm}}}&{{{\boldsymbol{{ M}}}_{ss}}} \end{array}} \right]} \right)\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_m}} \\ {{{\boldsymbol{U}}_s}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_m}} \\ {\mathbf{0}} \end{array}} \right\} \text{。}$

(17)

式中： $ {{\boldsymbol{U}}_m} $ 为主自由度，代表有外部载荷作用的节点自由度以及与圆柱壳接触的节点自由度； $ {{\boldsymbol{U}}_s} $ 为从自由度。基于式（17）中的第2组方程，从自由度 $ {{\boldsymbol{U}}_s} $ 可用主自由度 $ {{\boldsymbol{U}}_m} $ 表示为：

$ {{\boldsymbol{U}}_s} = - {\left( {{{\boldsymbol{{ K}}}_{ss}} - {\omega ^2}{{\boldsymbol{{\rm M}}}_{ss}}} \right)^{ - 1}}\left( {{{\boldsymbol{{ K}}}_{sm}} - {\omega ^2}{{\boldsymbol{{\rm M}}}_{sm}}} \right){{\boldsymbol{U}}_m}\text{，} $

(18)

将式（18）代入式（17）的第1组方程，可得：

$ \left[ {{{\boldsymbol{Z}}_b}} \right]\left\{ {{{\boldsymbol{U}}_m}} \right\} = \left\{ {{{\boldsymbol{F}}_m}} \right\}\text{。} $

(19)

其中：

$\begin{split} {{\mathbf{Z}}_b} =& \left( {{{\boldsymbol{K}}_{mm}} - {\omega ^2}{{\boldsymbol{M}}_{mm}}} \right) - \left( {{{\boldsymbol{K}}_{ms}} - {\omega ^2}{{\boldsymbol{M}}_{ms}}} \right)\times\\ &{\left( {{{\boldsymbol{K}}_{ss}} - {\omega ^2}{{\boldsymbol{M}}_{ss}}} \right)^{ - 1}}\left( {{{\boldsymbol{K}}_{sm}} - {\omega ^2}{{\boldsymbol{M}}_{sm}}} \right)\text{，}\end{split} $

(20)

进一步对主自由度 $ {{\boldsymbol{U}}_m} $ 进行分类，可表示为：

$ {{\boldsymbol{U}}_m} = {\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{b1}}}&{{{\boldsymbol{U}}_{b2}}} \end{array}} \right\}^{\rm T}}\text{，} $

(21)

式中， $ {{\boldsymbol{U}}_{b1}} $ 为内部结构上外部载荷作用的节点自由度， $ {{\boldsymbol{U}}_{b2}} $ 为内部结构与圆柱壳接触的节点自由度。将式（21）代入式（19），式（19）可改写为：

$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{Z}}_{b11}}}&{{{\boldsymbol{Z}}_{b12}}} \\ {{{\boldsymbol{Z}}_{b21}}}&{{{\boldsymbol{Z}}_{b22}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{b1}}} \\ {{{\boldsymbol{U}}_{b2}}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_{b1}}} \\ {{{\boldsymbol{F}}_{b2}}} \end{array}} \right\} \text{。}$

(22)

其中： $ {{\boldsymbol{F}}_{b1}} $ 为作用在内部结构上的外载荷； $ {{\boldsymbol{F}}_{b2}} $ 为内部结构与圆柱壳接触节点的耦合载荷。

1.3 圆柱壳与内部结构的耦合

在圆柱壳与内部结构的接触节点，位移和耦合载荷应满足：

$ {{\boldsymbol{U}}_{b2}} = {{\boldsymbol{U}}_0} \text{，}$

(23)

$ {{\boldsymbol{F}}_{b2}} = - {{\boldsymbol{F}}_2}\text{。} $

(24)

综合圆柱壳位移响应函数式（16）、内部结构缩聚动力学方程式（22），以及式（23）和式（24），经过简单的数学推导，整个结果的动力学方程可表示为：

$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{Z}}_{b11}}}&{{{\boldsymbol{Z}}_{b12}}{{\boldsymbol{D}}_2}} \\ {{{\boldsymbol{Z}}_{b21}}}&{{{\boldsymbol{Z}}_{b22}}{{\boldsymbol{D}}_2} + {\boldsymbol{I}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{b1}}} \\ {{{\boldsymbol{F}}_2}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_{b1}} - {{\boldsymbol{Z}}_{b12}}{{\boldsymbol{D}}_1}{{\boldsymbol{F}}_1}} \\ { - {{\boldsymbol{Z}}_{b22}}{{\boldsymbol{D}}_1}{{\boldsymbol{F}}_1}} \end{array}} \right\}\text{。} $

(25)

求解上式，可以得到圆柱壳与内部结构接触节点的耦合力 $ {{\boldsymbol{F}}_2} $ ，将其代入圆柱壳的控制方程，可求得圆柱壳上的振动响应，进一步采用轴对称边界元法计算圆柱壳的辐射声压。

2 结果讨论 2.1 方法验证

以水下含内部基座圆柱壳为例验证本文方法的正确性，圆柱壳长度0.8 m，半径0.2 m，厚度0.003 m；内部基座长度0.6 m（由 $ x = 0.1\;{\text{m}} $ 到 $ x = 0.7\;{\text{m}} $ ），宽度0.14 m，高度0.08 m，面板和2个腹板的厚度均为0.003 m。圆柱壳和基座由钢材制成，材料参数为：杨氏模量 $ {\text{2}}{\text{.1}} \times {\text{1}}{{\text{0}}^{11}}{\text{Pa}} $ ，泊松比0.3，密度 $ {{7\;850\;{\rm{kg}}/}}{{\rm{m}}^3} $ ，损耗因子0.01。圆柱壳的边界条件为两端简支。激励力F₁作用于基座面板中心点且垂直于面板，幅值为1N，分析频率为1～1 000 Hz，间隔1 Hz。选取场点P₁ $ (1,{0^{\circ }},0.4) $ 和场点P₂ $ (1,{90^{{\circ }}},0.4) $ 进行声辐射分析。

将本文方法与数值计算方法结果进行对比，以验证本文方法的正确性。在保证本文方法收敛的前提下，模态截断数M和N的值均取15。数值计算的过程为：在Ansys16.0中建立水下含内部基座圆柱壳的有限元模型（见图3），计算结构在水中振动响应；基于圆柱壳表面振动，在Virtual.Lab13.6中采用直接边界元法计算辐射声压。

本文方法与有限元法计算得到的基座面板中心点速度和圆柱壳均方速度（MQV）如图4所示。可以看出，除了部分峰值频率略有偏移，本文方法和有限元法的计算结果吻合良好，验证了本文方法计算水中含内部结构圆柱壳振动响应的正确性。

图5给出了2种方法在场点P₁和P₂的辐射声压对比曲线，本文方法与有限元耦合边界元方法计算结果吻合较好。为进一步对本文方法辐射声压预报结果的进行验证，图6给出了243 Hz和427 Hz的声压指向性曲线，2种方法的计算结果有较高的吻合度，验证了本文方法预报含内部结构圆柱壳辐射声压的正确性。

与有限元耦合边界元法相比，当内部结构变化时，本文方法不需要重新进行声-固耦合计算，只需刷新内部结构的动力学方程，可快速预报结构的辐射声压，这在内部结构变化的多方案对比分析时具有显著优势。

2.2 周向模态对辐射声压贡献量分析

图7给出了不同周向模态对应的场点P₁辐射声压对比曲线。可以看出，n=0∶3与n=0∶15的辐射声压曲线基本重合。这说明针对本文的研究模型，在1～1 000 Hz的频段内，结构向外场的辐射声压主要由 $ n \leqslant 3 $ 的周向模态贡献， $ n \geqslant 4 $ 的周向模态对声辐射的贡献可忽略不计。

3 结　语

本文以含内部结构圆柱壳为研究对象，结合解析法效率高和有限元法适用性广的优点，将整个结构分为水下圆柱壳和内部结构两部分。针对水下圆柱壳，采用轴对称边界元法计算壳体表面声阻抗，再结合Flügge壳体理论，建立水下圆柱壳的数理方程；针对内部结构，采用有限元法进行动力学建模。最后，根据圆柱壳与内部结构在连接点的内力平衡及位移连续条件，建立整个结构的控制方程。该方法可较为准确地预报含内部结构的水下圆柱壳的振动和辐射声压，且当内部结构发生变化时，不需要重新进行声-固耦合计算，只需更新内部结构的动力学方程，即可快速预报结构的声振特性，这在改变内部结构的多方案对比分析时具有显著的计算速度优势。