振动工程学报

阶数	干频率/Hz	振型(m,n)	偏差/%
1阶	59.49	60.70	(1,3)	2.03
2阶	75.44	77.41	(1,4)	2.61
3阶	86.44	85.75	(1,2)	-0.80
4阶	112.47	114.29	(1,5)	1.62
5阶	123.53	121.90	(2,4)	-1.32
6阶	126.88	125.94	(2,5)	0.74

阶数	干频率/Hz	振型(m,n)	偏差/%
1阶	59.49	60.70	(1,3)	2.03
2阶	75.44	77.41	(1,4)	2.61
3阶	86.44	85.75	(1,2)	-0.80
4阶	112.47	114.29	(1,5)	1.62
5阶	123.53	121.90	(2,4)	-1.32
6阶	126.88	125.94	(2,5)	0.74

阶数	湿频率/Hz	振型(m,n)	偏差/%
1阶	34.33	34.71	(1,3)	3.18
2阶	44.67	45.29	(1,2)	1.38
3阶	46.07	47.31	(1,4)	2.69
4阶	71.60	71.79	(1,5)	0.26
5阶	71.70	71.25	(2,4)	-0.62
6阶	77.20	76.45	(2,5)	0.97

阶数	湿频率/Hz	振型(m,n)	偏差/%
1阶	34.33	34.71	(1,3)	3.18
2阶	44.67	45.29	(1,2)	1.38
3阶	46.07	47.31	(1,4)	2.69
4阶	71.60	71.79	(1,5)	0.26
5阶	71.70	71.25	(2,4)	-0.62
6阶	77.20	76.45	(2,5)	0.97

动水中带封板圆柱壳振动的轴对称有限元-边界元耦合计算方法

PDF下载

张德春 ¹ , 陈浩然 ¹ , 李鹏 ¹ , 刘桂祥 ² , 杨翊仁 ¹

振动工程学报 | 2025,38(4): 687-696

收起

振动工程学报 | 2025, 38(4): 687-696

动水中带封板圆柱壳振动的轴对称有限元-边界元耦合计算方法

全屏

张德春¹, 陈浩然¹, 李鹏¹, 刘桂祥², 杨翊仁¹

作者信息

^1.西南交通大学力学与航空航天学院，四川　成都　611756

^2.中国核动力研究设计院，四川　成都　610213

张德春（1994—），男，博士研究生。E-mail：zhangdc_vib@qq.com

通讯作者:

李鹏（1983—），男，博士，教授。E-mail：lp_vib@126.com

Vibration analysis of cylindrical shells with endplates in axial flow：an axisymmetric FEM-BEM method

Dechun ZHANG¹, Haoran CHEN¹, Peng LI¹, Guixiang LIU², Yiren YANG¹

Affiliations

^1.School of Mechanics and Aerospace Engineering，Southwest Jiaotong University，Chengdu 611756，China

^2.Nuclear Power Institute of China，Chengdu 610213，China

出版时间: 2025-04-10 doi: 10.16385/j.cnki.issn.1004-4523.2025.04.003

文章导航

摘要

收起

针对水环境中带封板圆柱壳的动力学计算提出了一种轴对称有限元-边界元耦合计算方法。基于结构的周期性特征，圆柱壳和封板的运动可通过其母线的运动唯一确定。将母线离散为若干个单元以进行有限元求解。通过能量法获得单元矩阵并组装成总体矩阵，建立了结构运动的有限元计算格式；考虑流体的轴对称特征，在结构母线上建立轴对称边界元并构建了流弹压力的边界元计算格式；获得流体压力后将其转换为有限元节点的等效节点力，其表现为流体附加质量、阻尼和刚度；结合结构运动方程和流弹压力计算格式，构造了求解该类流-固耦合问题的轴对称有限元-边界元耦合计算方法。本文计算结果与已有理论解和商用分析软件计算结果吻合较好，表明了本文方法的准确性。相比于商用分析软件，本文方法具有更高的计算效率，可直接获得圆柱壳在流体中的质量和刚度矩阵，有利于受迫振动的快速计算。基于本文方法分析了端部含封板复杂圆柱壳在动水中的湿频率和稳定性特性，结果表明：封板会改变壳体表面流动特性，更易诱发圆柱壳的流弹失稳。

关键词

圆柱壳 / 流-固耦合 / 封板 / 轴向流 / 有限元-边界元耦合

Abstract

收起

An axisymmetric finite element-boundary element coupled computational method is proposed for the dynamics of cylindrical shells with endplates in water. Due to the periodic characteristics of the structure，the movement of the cylindrical shell and end plate can be determined by examining their meridional movement，and the meridians are discretized into several elements for further analysis. By using the energy method，the elemental matrices are obtained and assembled into the global matrices. Then the discrete schemes for structural motion with the finite element is established. Considering the axisymmetric characteristic of the fluid，an axisymmetric boundary element on the meridians is established and the calculation format for fluid-elastic forces on the finite elements is given. The fluid forces are then converted into an equivalent nodal force at the nodes and the effects of the fluid are evaluated by the added mass，damping，and stiffness matrices. By combining the structural motion and the schemes of fluid-elastic forces，an axisymmetric finite element-boundary element coupling method is developed to solve this fluid-structure coupling problem. The calculated results presented in this paper demonstrate good agreement with existing theoretical solutions and commercial analysis software，validating the accuracy of present method. Moreover，this method exhibits higher computational efficiency compared to commercial analysis software. It allows for direct determination of the mass and stiffness matrices of cylindrical shells in fluid，facilitating fast calculation of forced vibration. Based on this method，the wet frequencies and stability characteristics of complex cylindrical shells with end plates in axial flow are analyzed. The results show that the end plates will change the flow velocities on the cylindrical shell and make it more easily to instability.

Key words

cylindrical shell / fluid-structure coupling / endplates / axial flow / FEM-BEM coupling

引用本文

张德春, 陈浩然, 李鹏, 刘桂祥, 杨翊仁. 动水中带封板圆柱壳振动的轴对称有限元-边界元耦合计算方法. 振动工程学报, 2025 , 38 (4) : 687 -696 . DOI: 10.16385/j.cnki.issn.1004-4523.2025.04.003

Dechun ZHANG, Haoran CHEN, Peng LI, Guixiang LIU, Yiren YANG. Vibration analysis of cylindrical shells with endplates in axial flow：an axisymmetric FEM-BEM method[J]. Journal of Vibration Engineering, 2025 , 38 (4) : 687 -696 . DOI: 10.16385/j.cnki.issn.1004-4523.2025.04.003

正文

收起

薄壁圆柱壳是工程中广泛应用的重要结构，如导弹、储液罐、反应堆吊篮以及潜艇的外壳均可经简化而视为典型的薄壁圆柱壳结构。这些结构往往工作在流体环境中，其在水中的频率（以下称为湿频率）与结构的响应密切相关，对工程应用具有重要意义。针对圆柱壳湿频率计算最早提出的研究模型是轴向流（流动方向平行于圆柱壳轴线）中无封板的均质圆柱壳模型。LIN等^[1]采用行波法计算了静水中圆柱壳的湿频率，该方法也被推广至动水中圆柱壳的湿频率和稳定性计算^[2]。但需要指出的是，行波法只对长径比较大或两端简支的圆柱壳才有较好的适用性^[3]。PAIDOUSSIS等^[4-5]采用伽辽金法建立了轴向流中不同边界的圆柱壳的振动方程，并通过傅里叶变换计算流体力进而建立了流-固耦合方程，求解了其湿频率和临界流速。为考虑更一般的边界条件，DAI等^[6]通过在圆柱壳边界上布置弹簧模拟任意的弹性支承，结合假设振型和能量方法，分析了弹性支承对空气中圆柱壳频率（以下称为干频率）的影响规律，这本质上是一种半解析的方法^[7]。本文作者也提出了一种基于微分求积法（DQM）的计算方法^[8]，可有效地求解复杂边界圆柱壳湿频率。此外，有限单元法也是计算圆柱壳湿频率的一种重要方法^[9]。对于工程中常见的变截面（非均匀）圆柱壳，李海超等^[7]、CHEN等^[10]分别采用Jacobi多项式和行波法对圆柱壳的每段位移进行离散，通过力耦合和布置弹簧的方式考虑每段间的连续性条件，计算了变截面圆柱壳的干频率。然而据作者所知，暂无变截面圆柱壳湿频率计算的相关研究报道，因此有必要对此开展相关研究以解决该重要问题。

流体中的圆柱壳结构如潜航器^[11]及其舱段^[12]的端部通常带有封板。封板与圆柱壳的耦合建模是计算这类封闭结构频率的关键技术。石先杰等^[13]采用半解析法，分别假设圆柱壳和封板的振型，通过在连接处布置弹簧以模拟连接条件，通过能量方法计算了结构的干频率，与有限元法吻合很好。为考虑带封板圆柱壳与流体的耦合作用，CHENG^[14]联合使用半解析法和格林函数法计算了圆柱腔内的声辐射；CARESTA等^[11]采用行波法计算了含封板和肋板的圆柱腔声辐射问题；基于结构完整边界的广义谱-配点边界元法，张艳涛等^[15]求解了球-柱-球复合壳体的声辐射问题；XIE等^[16]结合半解析法和一维声学边界元法，计算了含肋板、锥壳和球壳封头及舱段间隔板的组合圆柱壳的响应和声辐射问题，该方法具有很好的适用性。但声学边界元方法需要通过特定处理消除矩阵的奇异性，计算较为复杂。陈美霞等^[17]也通过行波法计算了带封板圆柱壳的频响，但该研究没有考虑端板和水的耦合效应。此外，目前已有研究还主要是针对静水中的圆柱壳结构，暂未发现有文献考虑动水环境对含封板圆柱壳结构湿频率的影响，这不利于潜航器等在水中行进的结构的流-固耦合分析。

综上所述，目前仍缺乏针对动水环境中带封板的非均匀（如变材料、局部加固和变截面等）圆柱壳的求解方法和技术。值得指出的是，传统半解析法通常需使用弹簧模拟结构间的连接，当弹簧刚度过低时难以达到耦合效果，而过高时则会导致刚度矩阵病态从而影响求解精度，分析弹簧刚度的影响也会带来额外的计算量。正如庞福振等^[18]所指出的，为兼顾效率，结构壳段数划分不宜过大，因此该方法在任意非均匀圆柱壳的干频率和湿频率的求解中并未展示出足够的优势。有限元方法可通过共节点方式简单地实现不同结构间的耦合，但传统的有限元方法需要大量的单元以对圆柱壳结构进行离散^[9]，这会导致计算效率偏低。这一缺点是由于在建模中未关注到圆柱壳的轴对称特性及结构变形的周期性特征而造成的。事实上，基于这种周期性特征，圆柱壳上的运动可通过其母线上的运动进行表示^[8]，从而会很大程度上简化结构方程的建模并提高计算效率。TO等^[19]、SENJANOVIĆ等^[20]分别使用轴对称有限元计算了截锥壳和旋转圆柱壳的干频率，与理论解吻合很好，然而这种简便的建模和计算方法尚未被拓展到更复杂壳体组合结构的流-固耦合计算中。另外，注意到空气动力学中的边界元方法（面元法）对轴向流中薄壁类结构气动力求解表现出很好的适用性^[21]，因此本文拟在已有研究方法的基础上，发展轴对称有限元-边界元耦合计算方法以求解复杂的非均匀圆柱壳组合结构的流-固耦合振动问题，这也是本文主要的工作和创新点。

本文针对带封板非均匀圆柱壳的流-固耦合提出了轴对称有限元-边界元耦合计算方法。首先，对于圆柱壳沿其长度方向上截取的一段壳和在封板上截取的圆环板，考虑结构轴对称特征，这些结构的运动可通过母线上的运动表示。因此将母线离散为若干个有限单元。通过能量法获得不同周向波数下的单元矩阵，进而组装成结构总体质量和刚度矩阵，建立结构离散方程；然后，在结构的有限元离散网格上，扩展周期性条件至流体建模，构建一种轴对称边界元方法；计算各个边界元上的流体压力并将其转换为节点等效载荷，获得了流体影响的离散方程（表现为附加质量、阻尼和刚度矩阵）；最后，耦合结构和流体的离散方程，提出了动水环境中带封板圆柱壳的有限元-边界元耦合的计算方法。本文计算结果与已有理论解和商用分析软件计算结果吻合较好，且具有更高的计算效率。其基于轴对称有限元-边界元进行耦合求解，相比于传统的半解析法具有更强的适用性，可拓展到含环肋和多舱段复杂圆柱壳模型，也有望推广到变截面、变曲率等结构的动力学分析中。

1　结构模型及轴对称有限元方法

收起

如图1所示，本文研究的力学模型为一浸没在无黏无旋不可压缩的理想流体中的两端带封板的圆柱壳。外部流体密度为ρ_f，沿圆柱壳轴向以流速U_f流动。圆柱壳的厚度是可变的。圆柱壳及封板的材料相同，其密度、弹性模量和泊松比分别为ρ、E和ν。圆柱壳半径为R。坐标原点选取在圆柱壳轴线的中心，在原点分别建立直角坐标系和柱坐标系以描述位置和结构位移。圆柱壳轴向、环向和径向位移分别为U_s (x，θ，t)，V_s (x，θ，t)和W_s (x，θ，t)；封板径向、环向和法向的位移分别为U_p (x，θ，t)，V_p (x，θ，t)和W_p (x，θ，t)。考虑图1中蓝色粗线所围成的位于坐标轴x-y（x-r）轴形成的上半平面的阴影部分，该平面绕圆柱壳轴线旋转一周则可生成整个结构（轴对称结构），因此圆柱壳上与轴向平行的直线和封板上沿半径的径线在本文被定义为母线。圆柱壳和封板母线上的位移分别记作u_s、v_s、w_s和u_p、v_p、w_p。考虑这种周期性特征，圆柱壳和封板表面上任一点的位移均可分别表示为^[4-5]：

（1）

（2）

式中，n为周向波数。

上式表明，圆柱壳和封板的位移均可通过母线位移和周向的周期函数得到。传统的半解析方法常假设母线位移为三角函数或多项式函数。对于变截面等非均匀圆柱壳则需要对其每一段都进行假设，这会导致复杂的计算。事实上，考虑轴对称有限元方法^[19]，可将母线离散为多个单元（本文称之为轴对称有限元单元），每个单元的位移由单元节点位移进行插值计算获得，进而可建立结构轴对称有限元离散格式，这也是本文结构建模的基本思路。下面分别对圆柱壳和圆板结构进行离散，建立相应的有限元求解格式。

1.1　圆柱壳的离散格式

圆柱壳考虑为薄壳。如图2所示，沿轴线将其分割为若干个圆柱壳段。每段的位移可通过母线上（蓝线）的位移表示，将这段母线定义为一个轴对称单元^[19]，其长度为l，以两个节点的连线定义局部坐标s。值得注意的是，虽然轴对称有限单元定义为一维的，但其代表一段圆柱壳，可依据方程（1）描述该圆柱壳段上任意点的位移。根据有限元方法，母线单元上的任意一点位移可由节点位移经如下插值计算而获得：

（3）

式中，β_s为w_s关于x的1阶导数；ξ=s/l。插值函数选取为：

（4）

得到该段母线的位移后，由式（1）计算壳体位移，进而获得应变和应力。对于薄圆柱壳，根据Donell壳体理论，圆柱壳轴向、环向正应变和切应变与位移的关系为^[22]：

（5）

圆柱壳轴向、环向正应力和切应力通过其应变计算为：

（6）

式中，G为剪切模量。

由此可得该段圆柱壳的动能和势能分别为：

（7）

式中，dV=Rdzdθdx；h为该段壳体的厚度。

将式（1）和（3）代入式（7）即可获得系统能量的矩阵表达式：

（8）

式中，g=[u_s1 v_s1 w_s1 β_s1 u_s2 v_s2 w_s2 β_s2]^T为由单元节点位移构成的广义位移向量。圆柱壳轴对称单元的质量矩阵和刚度矩阵分别计算为：

（9）

式中，ρ_s为圆柱壳材料密度；(B_i_s)_3×8为与圆柱壳轴对称单元第i阶导数有关的矩阵（见附录A）；D为材料参数矩阵，

。

使用高斯求积法计算式（9）中的积分。需要注意的是，周向波数n=0表示轴对称的情况，此时节点的周向位移为0，可通过划去单元矩阵第2、6行和列上的元素而去掉该自由度。由式（7）中的积分可知，这种情况需将式（9）中的系数π变为2π。

1.2　封板的离散格式

考虑到封板为圆板，将其离散为若干个圆环，同样将母线离散为若干轴对称单元，如图3所示。单元内的位移可由节点上的位移插值计算获得：

（10）

式中，β_p为w_p2关于r的1阶导数。

圆环板径向、环向正应变和切应变与位移的关系为：

（11）

相应地，该单元的质量矩阵和刚度矩阵分别计算为：

（12）

式中，ρ_p为封板材料密度；r为单元中点处的半径；B_i_p为与单元第i阶导数有关的矩阵（见附录A）。对于n=0的情况，单元矩阵也要采取和圆柱壳单元同样的操作。

上述推导中，封板和圆柱壳轴对称单元的节点位移是定义在各自的局部坐标中的，在有限元法中需要计算总体坐标系下的单元矩阵以组装成总体矩阵。以图1中定义的柱坐标系为总体坐标系，母线上的点在柱坐标系下的x、θ和r方向上的位移分别为u、v和w，如图4所示，图中，τ和n分别表示切向量和法向量。注意到总体坐标系与封板局部坐标系的关系，最终得到转换后的封板轴对称单元矩阵为：

（13）

其中：

（14）

1.3　总体矩阵的组装

带封板圆柱壳的单元划分和节点编号顺序如图4所示，两侧封板均划分为N_p个单元，圆柱壳划分为N_s个单元，总单元数为N_m=2N_p+N_s，总节点数为N_n=N_m+1。圆柱壳和封板单元的单元矩阵分别通过式（9）和（13）进行计算（第i个单元的单元质量和刚度矩阵分别记为M_i和K_i）。在获得各个单元的单元矩阵后，系统的总体质量和总体刚度矩阵分别计算为^[23]：

（15）

式中，G_i为第i个单元节点自由度转换矩阵，它是行数和列数分别为8和4N_n的稀疏矩阵，矩阵中各非零元素定义为：

（16）

式中，G_i(j,4(i-1)+j)表示G_i的第j行、第4（i-1）+j列元素。即单元节点自由度转换矩阵非零元素的位置与单元和节点编号是相关的。由各节点自由度组成的未知位移矢量设为V。当不存在约束时，V的维度为4N_n；当某个自由度被约束时，需要划去总体矩阵中该自由度对应的行和列^[23]。

2　流-固耦合及轴对称边界元方法

收起

上一节获得了带封板圆柱壳的结构有限元离散方程，为考虑其与水体的耦合作用，下面建立流体作用力的离散计算格式。当圆柱壳端部存在封板时，流体的边界变得较为复杂，传统的一些解析方法很难适用^[4-5]。注意到空气动力学中广泛使用面元法求解具有复杂外形结构表面的气动力^[24]，但对于三维结构其需要较多的网格。考虑到流场也存在周期性，流体速度势和压力可通过其母线上的值获得。在圆柱壳和圆板母线上构建流体轴对称边界元，如图4所示，边界元的中点为控制点，在每个控制点处都需要满足无穿透条件^[24]：

（17）

式中，φ_i为第i个边界元控制点处的流体速度势；n为边界元的法线；v_N_i为第i个控制点即轴对称单元中点的位移，撇号和点号分别表示其对单元切向和时间的导数；U_i为第i个边界元控制点的切向流速；w_a_i为第i个控制点处的下洗速度。

为求解出满足边界条件的势函数，假设第j个边界元上均匀分布着源σ_j，对于任意边界元所代表的结构表面，源沿周向也是周期分布的。如图5所示，圆柱壳表面上一点P(x，y=R cos θ，z=R sin θ)处的源强度为σ_j cos (nθ)，则点P在第i个边界元的控制点Q (x_i，y_i，0)处产生的流体速度势为^[24]：

（18）

式中，r_PQ为P与Q两点之间的距离。

因此当边界元j为圆柱壳上的边界元时，其在Q点产生的速度势由以下的面积分得到：

（19）

式中，

，为第j边界元在第i个控制点处的势函数影响系数，式中积分采用高斯求积法计算。对于圆柱壳上的边界元，

，其中，x_i和y_i为第i个控制点的坐标。

与上式类似，当边界元j为封板上的边界元时，^pI_ij计算为：

（20）

此时，

。

由此获得各边界元之间的势函数影响系数。采用数值微分方法计算边界元在控制点处的方向导数，即可确定各边界元之间法向ⁿI_ij和切向诱导速度的影响系数^tI_ij。通过边界条件（17）确定各边界元的源强度后即可计算各个控制点的势函数，并由非定常伯努利方程计算各边界元控制点处流体压强，写成向量形式为：

（21）

式中，W_a为各边界元控制点处法向速度组成的向量；I_T、I_P和I_N分别为控制点的切向诱导速度、诱导势函数和法向诱导速度矩阵，其第i行第j列分别等于^tI_ij、^p I_ij和ⁿI_ij；上标“∘”表示矩阵Hadamard积，对于二维矩阵B，有：

（22）

式中，U为各控制点处的切向流速向量，其第i个元素表示第i个控制点的流速。U的计算式为：

（23）

式中，I为元素全为1的N_m维向量；I_T0和I_N0分别为周向波数为0时的切向和法向诱导速度矩阵；n₁为各控制点法线在x上的分量组成的向量。首先通过式（23）计算恒定流速时各个控制点的速度分布，再通过式（21）得到各控制点处的压强。

在面元法中，边界元上的压强被假设为均匀分布且其值等于控制点处压强^[24]。注意到流体压力只作用在壳体的法线方向，因此不会在轴向和周向的自由度上产生节点力。对于任一有限单元j，其上的压强为P_j。若其第i个自由度与单元法向位移相关，其等效节点力计算为：

（24）

式中，r(s)为边界元上各点半径的函数；N_i(s)为插值函数。

最终计算出等效节点广义力矢量为：

（25）

此时下洗速度向量W_a是未知的，为计算等效节点力矢量P，需要确定W_a和节点自由度向量V的关系，由式（17），有：

（26）

式中，T_p、T_d和T_v分别为控制点压强转换为等效节点力、节点位移转换为控制点处转角以及节点速度转换为控制点速度的转换矩阵（见附录B）。

利用上述结果，建立系统的结构-流体耦合振动方程为：

（27）

将压强的表达式（21）代入式（27），并移项得到：

（28）

式中，

和分别为流体附加质量、附加刚度和附加阻尼矩阵，这表明上述流体作用力表现为典型的流弹作用力。

式（28）给出了带封板圆柱壳与动水耦合振动问题的轴对称有限元-边界元耦合计算格式。此外，当母线上的节点受到外力激励时，只需要将式（28）中等号右端的零向量替换为等效节点载荷向量（详见参考文献[23]的2.2.5节）即可得到系统的受迫振动方程。

3　结果和讨论

收起

作为一种新的计算方法，首先需要评估该方法的收敛性和准确性。为了便于后续讨论，定义半波数m，对于圆柱壳其表示沿轴向的波峰和波谷总数；对于圆板其表示沿其直径的波峰和波谷总数。

3.1　收敛性分析

首先通过计算带封板圆柱壳的干频率讨论本文轴对称有限元方法的收敛性。选择与文献[8]一致的参数，圆柱壳长度为4 m，半径为1 m，均匀圆柱壳和封板的厚度分别为0.01 m和0.03 m。材料为钢，ρ=7850 kg/m³，E=2×10¹¹ Pa和ν=0.3。按照图4的方式划分单元进行求解，本文假设图4中的各个轴对称有限单元的长度相同。同时在ANSYS中使用SHELL181单元建立相应的有限元计算模型进行对比验证。

由于结构无约束，其前6阶频率都为0。工程上更关注非零频率。改变有限单元数目考察非零频率的变化。通过本文方法计算其前3阶非零频率和模态，并与ANSYS进行对比。由于圆柱壳的对称性，对于同一圆柱壳的振型，ANSYS会计算出两个相近的频率，文中取第一个进行对比。ANSYS中计算模型的单元数为N_a，计算结果随N_a的变化如图6中的虚线所示。当N_a=4876时ANSYS计算结果几乎不再随N_a变化，认为此时结果是收敛的。两种方法结果对比如图6所示，随着单元数N_m的增加，本文方法结果收敛到特定值；本文方法与ANSYS计算的频率偏差非常小（图中最大的偏差小于2%）。偏差是由于不同的壳理论以及数值计算方法导致的^[8]。

从图6中可以看出，本文方法只需要48个单元便可获得收敛的结果，而ANSYS则需4876个单元，本文方法计算需要的单元数远小于ANSYS。为保证计算结果的准确性，后文计算中选取N_m=60。此时本文方法计算200阶频率耗费的时间为0.24 s（Intel Xeon Gold 3.0 GHz处理器），而ANSYS则需要13.69 s，这也表明了本文方法的高效性。

3.2　空气中工况对比验证

为进一步讨论本文方法对非均匀结构问题的适用性，将圆柱壳中间2 m部分的厚度设为0.02 m（变截面情况），计算圆柱壳前6阶非零频率，与ANSYS对比结果如表1所示。本文方法与ANSYS最大偏差不超过3%，这表明了本文方法的准确性，本文方法可快速地计算变截面圆柱壳的频率，对该类结构的优化问题具有良好的应用潜力。此外，通过改变单元矩阵的计算方法，本文方法也有望拓展到包含球壳、锥壳等结构的组合壳的频率计算中。

3.3　水中工况对比验证

结构干频率对比结果验证了本文方法的有限元部分的准确性，下面通过水中湿频率计算验证本文方法的边界元和耦合计算的准确性。针对已验证过的变截面圆柱壳模型，考虑其置于静水中的情况。在ANSYS中使用声学单元FLUID30建立足够大的水体，并通过共节点实现壳和水体的耦合。计算前6阶壳体频率并与本文方法对比，如表2所示。结果表明本文方法与ANSYS偏差很小，表明了本文轴对称边界元方法对流弹力附加质量部分计算的准确性。值得注意的是，本文方法计算200阶湿频率所耗费的时间为1.21 s，而ANSYS则需要268 s，本文方法可有效地提高带封板圆柱壳湿频率的计算效率。

3.4　强迫振动验证

下面对空气和静水中含封板圆柱壳的强迫振动进行求解，以进一步考察本文方法的适用性。相比于传统的半解析法及商业软件，本文方法可直接获得结构的受迫振动方程，即方程（28）（右端项不为0时）。对其进行求解即可得到给定周向波数下的各个节点的稳态响应^[25]；对不同周向波数下的响应进行叠加可得到节点的位移响应。考虑单点激励的情况，如图1所示，单位法向力F_e加载在点P_e (0，R，0)（圆柱壳的中点）处。该点的法向位移的稳态解计算为：

（29）

式中，d_n (P_e)为周向波数为n时式（28）中该节点法向位移的稳态解；N_C为周向波数的截断数。选取表1对应的计算模型并设置结构损耗因子为0.02^[10]，N_C=20^[12]。

本文方法的计算结果与ANSYS的计算结果对比如图7所示。本文方法的计算结果与ANSYS的计算结果表现出了良好的一致性。值得注意的是，本文方法计算水中的100个激励频率下的响应只需要34 s，而ANSYS则需要耗费1806 s，表明本文方法具有更高的计算效率。

由图7可知，空气和静水中圆柱壳响应的峰值分别与其干频率和湿频率（分别如表1和表2所示）密切相关，这也表明了圆柱壳湿频率的计算对水下结构振动响应分析具有重要意义。

3.5　动水对湿频率的影响规律

以上算例表明了本文方法对圆柱壳在空气和静水中自由和受迫振动计算的准确性，同时也验证了本文中结构质量、刚度以及流体附加质量矩阵计算的准确性。下面通过讨论动水中结构湿频率对流体附加阻尼和附加刚度两部分的计算进行验证。由于尚未有文献考虑带封板及变截面圆柱壳在动水中的频率计算问题，为了与文献[8]进行对比，使用图6对应的均质圆柱壳模型，不考虑封板，对圆柱壳两端边界施加固支约束。

需要注意的是文献[8]考虑了圆柱壳内外流动，对于本文模型需要将文献[8]中的外部流体区域设置为无限大，且将内部流体密度设置为0。两种方法的计算结果对比如图8所示。由图8可知，对于不同模态，当流速较大时，两种方法计算的湿频率具有较小的偏差。这可能是由于本文计算考虑边界元上所有点的压强均等于控制点压强，不满足文献[8]中壳端部压强为零这一假设而导致的。工程中一般要求圆柱壳工作流速小于临界流速（各频率曲线与横轴的第一个交点对应的流速^[8]），此时两种方法的计算结果展现了较好的一致性。图8表明了本文方法对动水中圆柱壳频率计算的有效性和对多种边界条件的适用性。

下面通过本文方法计算和讨论动水对带封板圆柱壳结构湿频率的影响规律。为便于对比，仍使用图6中的圆柱壳模型。由于封板的作用，圆柱壳每个边界元控制点处的流速不再相同，此时需要先通过式（23）计算各个控制点处的切向流速（圆柱壳表面不同位置处流速U (x)与U_f比值如图9中带圆点实线所示），之后再通过本文方法计算结构频率随流速的变化，结果如图9所示。当不考虑封板对流场的影响，即假设壳表面液体的流动是均匀的且都等于U_f，此时的计算结果也展示在图9中。和未带封板的两端固支的圆柱壳一样，流速增加会导致湿频率的降低；相比于图8中的两端固支约束，封板对圆柱壳的约束更弱，降低了湿频率和临界流速。对于不同振型，封板对流场的影响都会使湿频率随流速衰减得更快，这是封板使流体产生了绕角流动而诱发圆柱壳表面的平均流速增加而导致的。

4　结论

收起

针对水中带封板圆柱壳的湿频率计算问题，本文提出了一种基于轴对称有限元-边界元耦合的计算方法，同时求解了动水环境中带封板圆柱壳的湿频率和稳定性问题。在结构建模中，采用轴对称有限元方法对结构进行离散，有效地降低了计算量；在流弹力建模中，发展了相应的轴对称边界元方法，计算了有限单元上的流体压力并转换为结构上的等效节点作用力（其表现为流体附加质量、阻尼和刚度），从而建立了流-固耦合有限元-边界元耦合离散方程。通过多种工况验证了本文方法的收敛性和准确性。本文研究可得出以下结论：

（1）与传统计算方法相比，本文方法表现出了良好的收敛性，仅需较少的单元即可保证计算结果的精度，在对非均匀复杂圆柱壳结构（变截面等）的分析和计算方面具有显著的优势。

（2）本文方法适用于流体中圆柱壳湿频率和受迫振动的响应的计算，具有较好的精度和较高的计算效率，有利于结构受迫振动响应的快速计算。

（3）本文方法可用于评估动水环境（流体附加刚度与阻尼）对圆柱壳湿频率的影响，相较于现有基于有限元软件的流-固耦合分析而言（仅能评估流体附加质量的影响），本文方法在潜航器的隔振和水声辐射研究方面具有应用潜力。

基金

收起

国家自然科学基金资助项目(12072298)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

LIN

T C

，MORGAN

G W

. Wave propagation through fluid contained in a cylindrical，elastic shell[J]. The Journal of the Acoustical Society of America，1956,28(6):1165-1176.

[2]

KOZAROV

，MLADENOV

. Hydroelastic stability of coaxial cylindrical shells[J]. Soviet Applied Mechanics，1981,17(5):449-456.

[3]

XING

Y F

，LIU

，XU

T F

. Exact solutions for free vibration of circular cylindrical shells with classical boundary conditions[J]. International Journal of Mechanical Sciences，2013,75:178-188.

[4]

PAIDOUSSIS

M P

，CHAN

S P

，MISRA

A K

. Dynamics and stability of coaxial cylindrical shells containing flowing fluid[J]. Journal of Sound and Vibration，1984,97(2):201-235.

[5]

PAIDOUSSIS

M P

，NGUYEN

V B

，MISRA

A K

. A theoretical study of the stability of cantilevered coaxial cylindrical shells conveying fluid[J]. Journal of Fluids and Structures，1991,5(2):127-164.

[6]

DAI

，YANG

T J

，LI

W L

，et al. Dynamic analysis of circular cylindrical shells with general boundary conditions using modified Fourier series method[J]. Journal of Vibration and Acoustics，2012,134(4):041004.

[7]

李海超，庞福振，张航，等．阶梯厚度圆柱壳自由振动特性分析[J].振动工程学报，2020,33(6):1226-1233.

H C

，PANG

F Z

，ZHANG

，et al. Free vibration of stepped cylindrical shell[J]. Journal of Vibration Engineering，2020,33(6):1226-1233.

[8]

ZHANG

D C

，LI

，LU

，et al. Stability analysis of cylindrical shell in axial flow：a DQ-based approach and an instability prediction formula[J]. Ocean Engineering，2023,267:113198.

[9]

董宇，杨翊仁，鲁丽.同心圆柱壳之间含静水的流固耦合振动分析[J].核动力工程，2015,36(5):80-82.

DONG

，YANG

Y R

，LU

. Vibration analysis of coupled concentric cylindrical shells in contact with an annular fluid region[J]. Nuclear Power Engineering，2015,36(5):80-82.

[10]

CHEN

M X

，XIE

，XU

，et al. Wave based method for free and forced vibration analysis of cylindrical shells with discontinuity in thickness[J]. Journal of Vibration and Acoustics，2015,137(5):051004.

[11]

CARESTA

，KESSISSOGLOU

N J

. Structural and acoustic responses of a fluid-loaded cylindrical hull with structural discontinuities[J]. Applied Acoustics，2009,70(7):954-963.

[12]

D T

，SU

J P

，HUA

H X

. Effects of circular plates on forced vibration of orthogonally stiffened cylindrical shell in wavenumber-frequency domain[J]. Ocean Engineering，2023,268:113439.

[13]

石先杰，李春丽，蒋华兵.复杂边界条件下圆柱壳-环板耦合结构振动特性分析[J].振动工程学报，2018,31(1):118-124.

SHI

X J

，LI

C L

，JIANG

H B

. Vibration characteristics of an annular plate-cylindrical shell structure with complex boundary conditions[J]. Journal of Vibration Engineering，2018,31(1):118-124.

[14]

CHENG

. Fluid-structural coupling of a plate-ended cylindrical shell：vibration and internal sound field[J]. Journal of Sound and Vibration，1994,174(5):641-654.

[15]

张艳涛，靳国永，杨宇航，等．球-柱-球复合壳体结构声辐射建模及特性分析[J].哈尔滨工程大学学报，2023,44(1):138-145.

ZHANG

Y T

，JIN

G Y

，YANG

Y H

，et al. Acoustic radiation modeling and analysis of sphere-cylindersphere composite shells[J]. Journal of Harbin Engineering University，2023,44(1):138-145.

[16]

XIE

，CHEN

M X

. Effects of end plates on vibro-acoustic responses of conical and cylindrical shells[J]. Applied Mathematical Modelling，2022,105:454-474.

[17]

陈美霞，张聪，邓乃旗，等．波传播法求解低频激励下水中加端板圆柱壳的振动[J].振动工程学报，2014,27(6):842-851.

CHEN

M X

，ZHANG

，DENG

N Q

，et al. Analysis of the low frequency vibration of a submerged cylindrical shell with endplates based on wave propagation approach[J]. Journal of Vibration Engineering，2014,27(6):842-851.

[18]

庞福振，李海超，霍瑞东，等．基于Jacobi-Ritz法的旋转组合结构自由振动特性分析[J].振动工程学报，2018,31(5):827-836.

PANG

F Z

，LI

H C

，HUO

R D

，et al. Jacobi-Ritz method for analysis of free vibration characteristics of combined shells in revolution[J]. Journal of Vibration Engineering，2018,31(5):827-836.

[19]

C W S

，WANG

. An axisymmetric thin shell finite element for vibration analysis[J]. Computers and Structures，1991,40(3):555-568.

[20]

SENJANOVIĆ

，ĆATIPOVIĆ

，ALUJEVIĆ

，et al. A finite strip for the vibration analysis of rotating cylindrical shells[J]. Thin-walled Structures，2018,122:158-172.

[21]

ZHANG

D C

，LIANG

，LI

，et al. A numerical and experimental study on the divergence instability of an inverted cantilevered plate in wall effect[J]. Archive of Applied Mechanics，2020,90(7):1509-1528.

[22]

曹志远．板壳振动理论[M]．北京：中国铁道出版社，1989.

[23]

王勖成．有限单元法[M]．北京：清华大学出版社，2003.

WANG

X C

. Finite Element Method[M]. Beijing：Tsinghua University Press，2003.

[24]

KATZ

，PLOTKIN

. Low-speed Aerodynamics[M]. Cambridge：Cambridge University Press，2001.

[25]

GAO

，PANG

F Z

，CUI

，et al. Free and forced vibration analysis of uniform and stepped combined conical-cylindrical-spherical shells：a unified formulation[J]. Ocean Engineering，2022,260:111842.

2025年第38卷第4期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.16385/j.cnki.issn.1004-4523.2025.04.003

接收时间：2023-04-17
首发时间：2026-02-12
出版时间：2025-04-10

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2023-04-17
修回日期：2023-06-28

基金

国家自然科学基金资助项目(12072298)

作者信息

^1.西南交通大学力学与航空航天学院，四川　成都　611756

^2.中国核动力研究设计院，四川　成都　610213

通讯作者:

李鹏（1983—），男，博士，教授。E-mail：lp_vib@126.com

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/zdgcxb/CN/10.16385/j.cnki.issn.1004-4523.2025.04.003

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

阶数	干频率/Hz		振型(m,n)	偏差/%
阶数	ANSYS	本文方法	振型(m,n)	偏差/%
1阶	59.49	60.70	(1,3)	2.03
2阶	75.44	77.41	(1,4)	2.61
3阶	86.44	85.75	(1,2)	-0.80
4阶	112.47	114.29	(1,5)	1.62
5阶	123.53	121.90	(2,4)	-1.32
6阶	126.88	125.94	(2,5)	0.74