基于波传播方法和多元分析的正交各向异性圆柱壳振动特性研究

汪志强， 李学斌， 黄利华

(武汉第二船舶设计研究所，武汉 430205)

圆柱壳结构在工业上得到了广泛的应用。随着复合材料的使用，很多圆柱壳结构形式都能够采用正交各向异性圆柱壳的模型进行讨论[1-2]。这种圆柱壳在轴向和周向具有不同的杨氏弹性模量和泊松比，并满足一定的约束。很多学者对这种圆柱壳结构形式的动态特性进行了研究，如自由振动[3]，冲击响应研究[4]等。

在对圆柱壳的研究中，为了探索物理参数以及几何参数对于动态特性(如频率)，通常的做法是，先设定物理参数以及几何参数的范围，确定需要研究的因素(如厚度和半径比值)，通过变化这些因素来求解它对频率的影响。这样往往会固定某些参数以突出该因素的作用进行分析。这样的方法效率比较低，难以顾及其他参数的同时变化，因此对于振动特性的全局性分析存在一定的局限性。

本文基于波传播方法[5]讨论了正交各向异性圆柱壳在多种边界下的自由振动特性。给定一定的物理和几何参数之后，通过试验设计[6]的方法获得很多频率解，这样就得到了比较大的设计空间。为了讨论该设计空间的属性，本文引入了多元变量分析的概念[7]，讨论了这些几何参数、材料参数与频率的关联性。通过方差分析(Analysis of Variance, ANOVA)判断出这些参数对频率影响的重要程度，可以为圆柱壳的动态特性设计提供参考。自组织映射[8]是神经网络技术的一种，它在数据分析上具有重要的地位。本文还采用该技术对这个空间进行了数据分析以及可视化研究。

1 理论分析

1.1 频率求解

正交各向异性圆柱壳和坐标系示意图，如图1所示。圆柱壳半径R，长度L，厚度h。假定材料是正交各向异性，剪切弹性模量G。圆柱壳中面点在坐标轴方向的位移分别为u，v和w(向内为正)。

图1 正交各向异性圆柱壳及振型示意图

设圆柱壳的材料主轴和坐标系主轴重合，则满足如下关系式

vxEθ=vθEx

(1)

式中：Ex，vx和Eθ，vθ分别是轴向和切向的杨氏弹性模量与泊松比。应用Flügge经典壳体理论[9]，轴向压力P作用下正交各向异性圆柱壳的振动方程为

(2)

Dx,Dθ,Dxθ为拉伸刚度；Kx,Kθ,Kxθ为弯曲刚度

(3)

根据波传播方法，可假设位移函数形式为

u(x,θ,t)=U0e-iknsxcos(nθ)eiωt，

v(x,θ,t)=V0e-iknsxsin(nθ)eiωt，

w(x,θ,t)=W0e-iknsxcos(nθ)eiωt

(4)

式中：U0、V0和W0为振型幅值；kns为轴向波数；n为环向波数；t为时间；ω为圆频率(1/s)；i2=-1。将这个位移函数代入式(2)，可以得到如下矩阵形式的方程

(5)

矩阵A的元素如下

圆柱壳的振动频率就是这个齐次方程组的非零解，即要求系数行列式等于零。将行列式展开，得到如下关于圆频率ω2的3次方程

g6(ω2)3+g4(ω2)2+g2ω2+g0=0

(6)

式(6)可通过解析公式准确求解。通过式(6)求解圆柱壳的频率，还需要考虑圆柱壳两端的边界条件。本文采用波传播方法，即利用对应梁弯曲振动的边界条件[10]近似为圆柱壳的边界条件。这里考虑了3种边界情况：两端薄膜简支，两端固定以及简支-固定边界。这些边界条件对应的梁特征函数以及解的形式在表1中给出。表中m是轴向形成的振型波数。

表1 梁弯曲振动特征函数和波数解

1.2 多变量分析

在圆柱壳振动的研究和设计过程中，研究者往往都会通过变换不同的参数，获得它们对于频率的影响，并研究各参数的重要性。或者说，频率对哪个参数的变化更加敏感一些。通常还可以设计成图谱等形式。对于正交各向异性圆柱壳而言，结构参数比较多，特性变化也比较复杂。本文拟采用试验设计的方法，研究圆柱壳的材料特性以及几何参数对于自振频率的影响。本文主要关注两个事项，材料特性以及几何参数中，参数和频率的相关性如何，以及参数改变对频率变化的影响分析。

从前面的分析可知，可认为圆柱壳频率ω实际上是材料特性、几何参数、边界条件和外力的函数，可以表示为

(7)

多元统计分析就是以p个变量，N次观察数据所形成的数据矩阵X=[xij]，i=1,2,…,N；j=1,2,…,p为依据。矩阵X的第i行表示第i个样品的观察值，它是一个p维的向量。矩阵X的第j列表示对第j个变量的n次观测值，它是一个N维变量。

变量之间的相互依赖关系多元分析的一个重要内容，它研究1个或几个变量的变化是否依赖一些变量的变化。对于本文，就是研究频率变量和几何参数以及材料特性参数的变化依赖性。X矩阵中，两个变量(N维)a和b之间的相关性可以由如下公式获得[11]

(8)

式中，Cov(a,b)是向量a和b的协方差，Var(a)和Var(b)分别表示a和b的方差。ρab也称为Pearson相关系数。

了解输入参数和频率的相关性能之后，作为圆柱壳动态特性的设计者，还希望知道设计参数对于频率的影响程度，以更好帮助设计。ANOVA是数据分析的一种方法，它解决的基本问题是，通过数据分析弄清楚和研究对象有关的各因素以及各因素之间相互作用对于该对象的影响。ANOVA的研究可以放在一般线性模型的框架内，它给出的假设是：观测值独立并且成正态分布，所有这些观测值都有同样的方差，观测值的均值可以表示成为某些参数的线性组合[12]。多元方差分析是ANOVA的拓展，它的主要步骤包括：问题确定，选择变量，假设，模型估算，结果分析和验证等环节。这种技术已经是多元统计分析的重要部分，有很多统计类软件提供这个功能[13]。本文拟采用该技术研究几何参数以及材料参数各因素对频率的影响程度。

从上述过程可以看到，经过试验设计得到的解空间和性能空间是多维的。当我们考虑的系统超过3维之后，其可视化就会遇到很多困难。对于这个多维数据构成的设计空间，本文引入自组织神经网络(Self-Organizing Facture Map, SOM)技术进行了研究。SOM是一种无导师的神经网络学习方法。它最重要的特点是通过自动寻找样本中的内在规律和本质属性，自组织、自适应地改变网络参数与结构，广泛用于数据挖掘以及可视化研究。

SOM网络共有两层，即输入层和输出层。输入层的神经元和样本的维数相等。输出层也就是竞争层，通常的形式是二维平面阵。输出神经元常常采用矩形或者六角形侧向连接。在网络训练过程中，竞争获胜神经元和其周围一定半径内的神经元可以依据特定的形式调整权向量。SOM的算法过程和简要描述如下：

步骤1 初始化 权向量和学习率初始化，建立初始优胜邻域；

步骤2 竞争 样本输入以及计算获胜节点；

步骤3 合作 竞争优胜神经元确定获胜区域，确定邻域节点；

步骤4 调整 对优胜神经元获胜区域内的神经元进行节点权值调整；

步骤5 检查判断 检查学习率是否达到预设要求。否则转到步骤2。

SOM技术通过这两层网络形式能够把多维数据(输入层)转换到二维平面图(输出层)，它既是一种数据分析工具，也是一种可视化技术[14]。通过形成的二维图形，能够获得更多数据之间的内在关系。本文拟采用这个技术分析结构几何参数以及材料参数对于频率的影响，以及输入参数之间的关联关系。

2 数值算例和讨论

2.1 不同边界下圆柱壳的自振频率

这里首先使用波传播方法计算了不同边界和材料特性情况下，受不同轴向力作用时圆柱壳的频率，见表2。 对于正交各向异性圆柱壳，引入频率因子Ω2=ρR2ω2[15]。

表2 不同边界条件下正交各向异性圆柱壳的频率Ω，Ω×10-5

对于简支边界条件，表2还给出使用经典的解析解方法的结果[16]。对于这种边界，波传播方法和解析解的结果非常接近。对于这两个方法的比较，Gan等[17]针对环肋圆柱壳结构还给出理论分析，证明了两者之间的等价性。

2.2 多元统计和数据分析

这里讨论一个数值算例，计算圆柱壳在一定的材料特性和几何参数范围的频率特性。引入试验设计技术研究这些参数对于自振频率的影响。参数范围以及试验设计的水平列在表3。 假设不考虑外力，这样共有7个参数，每个参数选取3个水平，采用全因子设计方式。

表3 材料特性、几何参数范围和水平取值

通过这些参数的组合，并计及式(1)以及vθ<0.2的约束，可以得到1 782个数据。这样X矩阵的大小就是1 782×10。如果仅考虑最低频率，则X矩阵的大小就是1 782×8。

根据式(8)可以得到最低频率和这些参数的相关系数ρ以及p值见表4。计算中采用显著性水平α=0.05，p值是否定假设检验中原假设是否适当的参数。

从表4可知,L/R和最低频率的相关系数为负值，这说明壳体越长，则其最低频率越低。壳体越厚(h/R增加)以及轴向模态越高，频率也越大。两个坐标轴方向的杨氏弹性模量和频率的相关性相差不大。注意到表中vx和n对应的p值并不为零，这说明它们在置信度为0.95时并不是显著相关的。

表4 最低频率Ω1和参数的相关性指标以及p值

使用一般线性模型[18]，对于最低频率进行多元方差分析的结果见表5。表5中自由度等于变量的水平-1；Seq-SS(连续平方和)为组间平方和(因子)以及组内平方和(误差)，MS-平方和是连续平和和除以自由度得到的均方。f是将因子MS除以误差得到。p值是用以确定某个因子是否显著，通常和α=0.05进行比较。如果p<0.05，这表示该因子是显著的。

表5 最低频率多元方差分析结果

从表5可知，可以得到7个输入参数对于最低频率的重要程度。L/R的重要性达到5 055.44/17 916.62×100%=28.22%. 这说明圆柱壳的长度/半径比值对频率的影响最大。其次是轴向模态m以及厚度因素，分别达到了16.34%和7.05%。其他参数对于频率的影响重要性均比较小。

针对表3中的数据，因子范围保持不变，加大因子的水平到5，采用全因子设计共得到59 375个试验样本。以圆柱壳两端简支，轴向压力P=0为例，将这个59 375×10维的数据(含中间频率Ω2和最高频率Ω3)通过SOM技术映射到二维空间，可以得到图2。本文利用SOM Toolbox软件[19]进行分析求解。该映射图能够将正交各向异性圆柱壳的7个输入参数和3个频率(即共计10维的高维空间数据)用二维的形式展示出来，可以方便研究者从整体上了解振动特性。和传统用曲线表示两两变量间的变化方法相比，这里能更加直观给出高维空间里变量之间的关系。

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

在图2中，神经元以六角形侧向紧密连接，映射图中网格大小为47(纵向)×26(横向)，即59 375个数据共投影到这1 222个神经元内。当考察某属性映射图里(比如L/R图)单个神经元特性的时候，需要同时计及该神经元在其它映射图中(对应坐标网格)的属性(比如h/R等)。

在SOM映射过程中的权向量是多维的，有多种表示权向量差别的度量方法，主要有基于距离的统一距离矩阵(U-matrix)和基于密度估计的密度度量矩阵(P-matrix)方式。图2中给出U-matrix结果以及各个变量的变化情况。U-matrix矩阵中存储的是对应位置节点和相邻节点权向量之间的距离。从这个图形能够依据距离对原始空间的数据点进行聚类分析。本文形成SOM输入数据的是通过全因子分析得到的，每个参数给出的水平值相同，在U-matrix图像上显示的聚类信息并不显著。

观察7个输入变量以及3个频率参数的图像，可以分析如下：

1)最低频率Ω1(图2(i))在整个范围内变化均比较小，Ω1远小于1。在右下角有比较集中的大值，这对应于轴向波数m较大、h/R较大以及壳体较短(L/R较小)的情况。此区域的Ex值要远大于Eθ。结合表5的数据可知，Ex的重要程度也大于Eθ。比较L/R和Ω1的分布图，壳体越短，在其他参数不变的情况下，频率就越高。分析这两个图从左上角变化到右下角的趋势可以看出，这和表4中L/R与Ω1的相关系数-0.48(负值)是对应的(大小变化顺序相反)。轴向模态m对最低频率Ω1的影响，与L/R的影响类似，不过是正向的，即m越大，则Ω1也越大。这与m和Ω1的相关系数0.40(正值)也是符合的。

2)周向波数n对应Ω1的影响并不显著，除了从相关系数可以看出外，对比图2(i)可知，n值的取值范围内可以得到比较小最低频率。在最左下角处，即使是n值很大，结合h/R，L/R值，频率也比较低。

3)对于各向异性材料，为了获得较高的基础频率，在几何参数确定之后，比较图2(i)的左半幅和图2(f)与图2(g)，Ex和Eθ的增加导致频率增加恰恰是比较小的。Ex和Eθ的配合需要协调，并不是越高越好。这对于材料设计和制造经济性而言是有启示意义的。

4)圆柱壳的最高频率Ω3和最低频率Ω1有相似的变化趋势。最大值出现在右下角，壳体越短，m越大则越高。不过，在壳体比较薄(h/R小)以及低周波数n时，也会出现比较大的最高Ω3频率。对应于轴向位移为主的中间频率，其变化规律和Ω3相近。

3 结 论

本文根据Flügge经典壳体理论以及波传播方法讨论了圆柱壳在不同边界下的自振频率特性。引入试验设计的思想，结合多元统计分析方法以及自组织映射技术，深入分析了几何参数以及物理特性参数对于受到轴向压力作用正交各向异性圆柱壳的频率影响。本文给出的方法能够有效探知设计参数对于频率的影响，可以为圆柱壳的动力设计提供有效手段。从上述过程可以得到如下结论：

(1) 壳体的长度/半径比值和最低频率的相关系数为负值。壳体越长则其最低频率越低。壳体越厚(h/R增加)以及轴向模态越高，频率也越大。

(2) 长度/半径比值对频率的影响最大。其次是轴向模态m以及厚度因素，其他参数对于频率的影响重要性均比较小。

(3) 通过全因子分析得到的数据，在U-matrix图像上显示的聚类信息并不显著。

(4) 最低频率Ω1在整个范围内变化均比较小且远小于1。比较集中的大值对应于m较大、h/R较大以及壳体较短(L/R较小)的情况。Ex的重要程度也大于Eθ。周向波数n对应Ω1的影响并不显著，即使是n值很大，结合h/R，L/R值，频率也可能比较低。

(5) 对应各向异性材料，对于动态特性设计，Ex和Eθ的配合需要协调，可提高制造经济性。

[1] LEISSA A W. Vibration of shells (NASA SP 288), NASA SP 288[R]. Washington DC, US: Government Printing Office, 1993.

[2] 赵鑫, 张博, 李跃明. 不同边界条件下截顶正交各向异性圆锥壳的振动和声辐射研究[J]. 振动与冲击, 2016,35(3):99-106.

ZHAO Xin,ZHANG Bo,LI Yueming.Influences of boundary conditions on the vibration and sound radiation of a truncated orthotropic conical shell[J].Journal of Vibration and Shock,2016,35(3):99-106.

[3] LIU B, XING Y F, QATU M S, et al. Exact characteristic equations for free vibrations of thin orthotropic circular cylindrical shells[J]. Composite Structures, 2012,94:484-493.

[4] 李学斌. 正交各向异性圆柱壳的稳态动力响应分析[J]. 船舶力学, 2007,11(1):79-87.

LI Xuebin. Hairmonic response analysis of orthotropic cylindrical shells[J]. Journal of Ship Mechanics, 2007,11(1):79-87.

[5] LI Xuebin. Study on free vibration analysis of circular cylindrical shells using wave propagation[J]. Journal of Sound and Vibration, 2008,311:667-682.

[6] 马成良, 张海军, 李素平. 现代试验设计优化方法及应用[M]. 郑州: 郑州大学出版社, 2007.

[7] JOHNSON R A, WICHERN D W. Applied multivariate statistical analysis[M]. 6th. New Jersey: Prentice Hall, 2007.

[8] KOHONEN T. Self-organizing maps[M]. 3rd. Berlin: Springer-Verlag, 2001.

[9] FLÜGGE. Stresses in shells[M]. Berlin: Springer-Verlag, 1973.

[10] 季文美, 方同, 陈松淇. 机械振动[M]. 北京: 科学出版社, 1985.

[12] DEGROOT M H. Probability and Statistics[M]. Wokingham, England: Addison-Wesley Publishing Company, 1989.

[13] 高慧璇. 应用多元统计分析[M]. 北京: 北京大学出版社, 2005.

[14] MWASIAGI J I. Self organizing maps-applications and novel algorithm design[M]. Rijeka, Croatia: InTech, 2011.

[15] WARBURTON G B, SONI S R. Resonant response of orthotropic cylindrical shells[J]. Journal of Sound and Vibration, 1977,53(1):1-23.

[16] 李学斌. 正交各向异性圆柱壳静动态特性分析及比较研究[D]. 武汉: 华中科技大学, 2004.

[17] GAN L, LI X, ZHANG Z. Free vibration analysis of ring-stiffened cylindrical shells using wave propagation approach[J]. Journal of Sound and Vibration, 2009,326:633-646.

[18] TIMM N H. Applied multivariate analysis[M]. New York: Springer-Verlag, 2001.

[19] VESANTO J, HIMBERG J, ALHONIEMI E, et al. Self-organizing map in MATLAB:the som toolbox:[C]//Proceedings of the Matlab DSP Conference. Espoo, Finland, 1999.