输流圆柱壳固有特性分析的有限元传递矩阵法

陈爱志，王浩然，柳贡民

（1.海装沈阳局驻哈尔滨汽轮机厂有限责任公司 军事代表室，哈尔滨 150046；2.耐世特汽车系统有限公司，江苏 苏州 215000；3.哈尔滨工程大学 动力与能源工程学院，哈尔滨 150001）

圆柱壳结构常见于船舶、航空航天、水利水电等工程领域。结构振动直接影响圆柱壳结构的可靠性。圆柱壳结构固有特性的计算分析是结构优化、安全性设计中的重要环节。因而，准确、高效地完成圆柱壳结构的固有特性计算就成为人们关注的重要研究内容。

结构固有特性的分析方法有很多，包括当今应用广泛的有限元法（FEM）[1]，能高效计算链式结构的传递矩阵法（TMM）[2]，将二者结合起来的有限元传递矩阵法（FETMM）[3]以及常用于瞬态应力分析的特征线方法（MOC）[4]等。其中，FETMM是由Dokanish[5]在1972年提出的，该方法将FEM和TMM联合起来进行结构的分析，拓宽了传统TMM的研究范围并改善了FEM在使用过程中计算效率较低的问题。FETMM处理复杂结构时既具有FEM的优点，又兼备TMM的高效率[6]。徐铭陶[7]创建了一种新的始端线和终端线节点数目不等的传递元件，提出了广义传递矩阵和有限单元组合法，此法便于分析复杂形状平板的振动。何斌、芮筱亭等[8]基于欧拉梁理论，用FETMM对细长弹箭的振动特性进行了计算。当前，对FETMM的利用局限于对梁和板的分析，较少出现在圆柱壳的结构振动特性的分析中。不过根据Petyt[9]的研究，可以通过对圆柱壳周向模态的限制，利用类梁的壳单元来减少自由度，来达到简化计算的目的。张勇亮[10]利用此理论对内部含有流体的圆柱壳耦合振动波传播问题进行了研究。

本文利用圆柱壳的轴对称特性，基于Sanders薄壳理论[11]建立类梁的壳模型，运用MATLAB软件完成了有限元传递矩阵法程序的编写和圆柱壳固有特性的计算，验证了有限元传递矩阵法与类梁的壳模型结合运用于输流圆柱壳固有特性研究的准确性和高效性。

1 基于类梁壳模型的圆柱壳流固耦合方程

图1 含流圆柱壳模型

参考文献[12] ，将圆柱壳3个方向的位移进行傅里叶余弦展开,具体表达式如下:

其中，n为周向波数，ux为轴向位移，uθ为径向位移，ur为周向位移。轴向和径向位移的模态假设为cosnθ，周向假设为sinnθ[9]。将位移进行如此变换，则只需进行轴向网格划分而不必进行周向网格划分，如此，有限元方程的自由度得到降低，计算量和计算时间大大减少。这样划分所得单元称为类梁的壳单元。

图2分别为基于壳单元的一般圆柱壳网格划分模型与基于本文类梁壳单元的圆柱壳网格划分模型。

图2 两种壳体模型

在求解过程中，每个节点各方向位移的具体形式可以通过相应的形函数和节点位移呈现[12]

因此，单元位移表达式可以简单地写为

[Ns] 为形函数，{uˉ}为节点位移矢量。每个单元有2个节点，每个节点包含4个位移，即轴向uˉx、周向uˉr、径向uˉθ和一个旋转量φ=∂uˉr/∂x，节点位移矢量和形函数具体形式如下

其中：

ξ=x∕le，le是单元的长度。此处，圆柱壳的径向位移类似梁的横向位移。因此，方程式(4)中的形函数和梁横向振动的形函数相同。

参考文献[12] ，圆柱壳的单元质量、刚度矩阵和力矢量可以分别表示为

其中：ρs为壳体密度，h为壳体厚度，R为壳体中面半径，[B ] 为应变位移矩阵，[D] 为应力矩阵，pd为流体作用于壳壁的压力，ps为外部分布力。

另外，考虑管内流体为不可压缩、非黏性流体时，极坐标系下的拉普拉斯方程可以表示为[13]

对于管内以恒定速度流动的流体，其单元质量、刚度矩阵和耦合阻尼矩阵可表示为[12]

hf表示流体的有效厚度，它取决于周向模数n和特征值λ。

流体作用于圆柱壳内壁的压力可以用势函数来表达，其中U为流体速度。

2 有限元传递矩阵法

将输流圆柱壳整体沿轴向离散成若干单元，则单元的流固耦合振动微分方程可以表示为[10]

假设圆柱壳振动为简谐振动形式，则单元位移列向量以及单元力向量可表示为

ω为管道振动圆频率。因此，对式(12)进行变换得

式(15)可以简写成如下形式

其中

分块以后得

对式(16)进行恒等变换，写成如下形式

其中

由此得圆柱壳轴向划分为n个单元时，前后端的状态矢量传递关系满足

从输入端矢量到输出端矢量的整体传递关系可用传递矩阵Ut表示

值得注意的是，在运用有限元传递矩阵法求解圆柱壳振动问题时，由于矩阵的连乘，会出现计算不稳定的情况，导致结果出现较大误差。本文引入Riccati变换后，可以较好地解决这一问题。图3给出了引入Riccati变换前后的计算结果对比。

图3 Riccati传递矩阵法与传统传递矩阵法对比

其中，纵坐标Log(miu)为编程求解过程中所设的参数，其峰值所对应的横坐标频率即为计算所得的圆柱壳固有频率。

3 计算精度影响因素分析及算例验证

3.1 计算精度影响因素分析

如前所述，将类梁壳模型与Riccati有限元传递矩阵法结合进行圆柱壳固有频率求解。作为对比算例[13]的薄壁圆柱壳与流体具体参数如表1所示，边界条件为一端固支一端自由。

表1 输流圆柱壳与流体参数

在求解数学问题时，选择恰当的方法是顺利求解的关键。本文运用3、4、5次高斯积分分别进行计算，但是在求解过程中发现运用3次高斯积分所得结果并不理想，具体结果对比如图4和表2所示。

从图4和表2的计算结果很容易发现，相对于3次高斯积分，4次高级积分和5次高斯积分得到的计算结果更接近于实验值和商业软件仿真结果。因此，在接下来的计算分析过程中采用计算结果更精确的5次高斯积分进行求解。

表2 不同高斯积分阶数计算结果对比

简便和高效是本文方法的突出特点，但在求解过程中发现轴向划分的单元数目会影响算法的准确性和效率。图5展示了同一模态下划分不同轴向单元数时圆柱壳固有频率计算结果对比，表3列出不同单元数下固有频率的同时，也列出了计算所耗费的时间。

可以看出单元数目太低会导致计算结果准确性下降，网格数增加虽能提高计算精度，却又会导致运算时间增加。

表3 不同单元数下的固有频率计算结果与耗时对比

3.2 算例验证

在3.1所选定的高斯积分阶数(5阶)及轴向单元划分数(20个)的基础之上，利用本文方法对圆柱壳的固有特性进行计算。此外，采用商业软件ANSYS对同一圆柱壳振动固有频率进行同步求解，仿真计算中选取的单元类型为SHELL63。

图4 不同高斯积分阶数计算结果对比

图5 不同单元数下的固有频率计算结果对比

图6所示为计算得到的不考虑内含流体时不同周向模态下的圆柱壳振动固有频率；表4、图7表现的是实验值、有限元软件所得数据与本文方法计算结果的对比，从中可以看出本文模型和方法的准确性。

图6 不同周向模态下的圆柱壳固有频率

图7 固有频率对比

表4 不含流圆柱壳的固有频率

4 流体对圆柱壳固有特性的影响

分析圆柱壳固有特性与其内含流体是否有着密切的联系。

流体对圆柱壳的作用包括附加质量、阻尼和刚度。其中，作为附加质量的流体质量是变化的，这种变化和周向模态有关。

图8表示的是不同周向模态下流体作为附加质量作用于圆柱壳时的有效厚度值。

图8 流体在不同周向模态下的附加质量有效厚度

从图8中可以发现当n=1时，有效厚度值为1即流体完全作用于附加质量。随着n的增加,流体作用于圆柱壳的有效厚度逐渐变小，变化也是由快到慢的。

对同一圆柱壳模型，内含静止流体时固有频率相对不含流体时的固有频率值会降低，这是因为内含的流体相当于增加了圆柱壳的质量，而圆柱壳固有频率会随着其质量的增加而减小。含静止流体圆柱壳固有频率计算结果与实验值以及有限元值的对比见表5和图9。

表5 含静止流体圆柱壳固有频率

图9 含静水固有频率对比

图10显示了n=1时含流与不含流圆柱壳固有频率计算结果对比，可以看出含流时同一频率区间内固有频率(Log(miu)的峰值)个数明显增加。

图10 n=1时固有频率对比

这是因为这些固有频率的计算结果之中不仅包含了含流圆柱壳整体振动的固有频率，还包含了内含的水柱的固有频率。

接下来考虑圆柱壳内含的流动流体对其固有特性的影响。采用的算例来自于Weaver[12]，具体参数如下：ρf∕ρt=0.128，h∕R=0.01，L∕R=2。无量纲流速Uˉ =U∕U0，无量纲频率为ωˉ=ω∕ω0。

本文得到的结果和Weaver以及Selmane[13]得到的结果对比如图11所示。可以很明显看出圆柱壳固有频率随着流速的增加而降低。同时从对比的结果可以发现，流速较低时，本文方法和Weaver以及Selman得到的结果很相近，但随着流速的增加，对比周向模态数m=2时的结果可以发现，本文方法与Weaver出现较大不同，但与Selman文中的结果吻合相对较好。这可能是由于Weaver在应用伽辽金法时采用的项较少。总体来看，本文方法在高流速时计算结果的精确度没有流速较低时好，这可能与本文处理圆柱壳内流体的方式有关。

图11 无量纲频率计算结果对比

5 结语

本文结合类梁壳模型及有限元理论，考虑圆柱壳与流体的耦合作用，采用有限元传递矩阵法，对输流圆柱壳的固有振动特性进行了分析。扩展了传递矩阵法的研究对象，同时在保持有限元计算精度的前提下，提高了计算速度与效率。通过与实验和仿真数据的对比验证了本文模型和算法的准确性，同时研究了高斯积分次数以及划分的轴向单元个数对计算精度的影响，以及流体对圆柱壳振动固有频率的影响

计算结果表明随着高斯积分次数和轴向划分单元数的增加，计算精度也会随之增加，但增加的幅度会越来越小，且计算时间也会相应延长，因此在满足一定计算精度要求的前提，应使2个参数尽量的小以提高计算效率。

含流与不含流圆柱壳对比，同一频率区间内固有频率个数明显增加，同时圆柱壳振动固有频率随着流速增加而降低。