杜 圆， 李海超， 庞福振， 缪旭弘,2

(1.哈尔滨工程大学 船舶工程学院，哈尔滨 150001； 2.海军研究院，北京 100161)

矩形板结构在船舶、海洋平台、潜艇及波浪能电站等结构中有着大量的应用，因而设计者需对其振动特性进行详细了解。目前常用的方法有能量法[1]、微分求积法[2]、辛本征展开法[3]。但使用能量法的前提是选取与边界条件相对应的挠度函数，微分求积法中加权系数的选择对结果影响很大，辛本展开法所建立的泛函缺乏与边界条件相关的表征项。Shen等[4]通过使用坐标函数的线性组合来表征弯曲位移函数得到了板结构的弯曲振动近似解；Cho等[5]基于假定模态法，通过对板的拉格朗日方程求解特征值求解板的自由振动问题；Reddy[6]基于Reissner-Midlin剪切变形理论与非线性应变位移关系对板的弯曲振动进行求解；Aksu等[7]运用有限差分方法对矩形板固有频率进行了计算；Gorman[8]提出用系统叠加法来求解各种边界条件下板的振动问题。上述研究求解过程复杂，且仅能对经典边界条件求解，无法对实际工程应用中常见的弹性边界条件进行求解。

本文将薄板振动的位移容许函数表示成二维傅里叶余弦级数和辅助傅里叶级数的线性组合，辅助傅里叶级数的引入使位移容许函数三阶连续且可偏导；从而满足薄板结构振动控制微分方程，解决了传统傅里叶级数法在边界处不连续的问题。基于改进后的位移容许函数列出矩形薄板拉格朗日能量泛函方程，然后通过Hamilton原理求解得到位移容许函数中的系数。通过改变横向位移约束弹簧刚度值k和旋转约束弹簧刚度值K模拟任意边界条件。

在与文献及数值仿真结果进行对比验证本文方法的有效性后，对弹性组合边界条件下矩形薄板的自由振动特性进行求解。并通过大量计算探讨边界条件与矩形薄板自由振动特性的关系，旨在为后续研究及工程应用提供参考。

1 矩形薄板振动理论模型

1.1 模型简介

本文所研究的任意边界条件下矩形薄板示意图如图1所示。矩形薄板的边界条件通过四条边上的横向位移约束弹簧与旋转约束弹簧来模拟。

图1 任意边界条件下的矩形薄板示意图

横向位移约束弹簧与旋转约束弹簧在四条边上均匀分布，其中kx0、ky0、kxa、kyb分别代表x=0,y=0,x=a,y=b边界上横向位移约束弹簧的刚度，Kx0、Ky0、Kxb、Kya分别代表x=0,y=0,x=a,y=b边界上旋转约束弹簧的刚度。当kx0、ky0、kxa、kyb、Kx0、Ky0、Kxa、Kyb均趋于无穷时，表征矩形薄板四边刚性固定；当kx0、ky0、kxa、kyb、Kx0、Ky0、Kxa、Kyb均设置为零时，表征矩形薄板四边自由；当kx0、ky0、kxa、kyb刚度设置为无穷,Kx0、Ky0、Kxa、Kyb设置为0时，表征矩形薄板四边简支。

1.2 位移容许函数

位移容许函数的选择将直接影响后续计算精度[9]，求解过程中需将式(1)代入结构控制方程式(3)。但传统傅里叶级数的导数在端点处仍会产生不连续或跳跃，如图2所示。通过将位移容许函数表征为式(1)的形式，辅助函数的引入使其满足求解域内三阶求导连续且四阶导数存在，有效克服边界不连续现象，如图3所示。

式中:λm=mπ/a;λn=nπ/b;Amn、Bmn、Cmn和Dmn为矩形薄板位移容许函数中的Fourier系数，简谐时间因子eiωt表示薄板位移随时间的变化。式(1)中第四项傅里叶正弦级数之和用来处理边界不连续，前三项中的m与n取值越大，位移容许函数越精确；在计算过程中展开级数项m与n的取值M、N称为截断值。

图2 传统傅里叶级数在端点位置可能产生不连续问题示意图

图3 改进傅里叶级数如何解决端点位置可能产生不连续问题示意图

1.3 结构控制方程

薄板仅承受法向载荷q(x,y,t)时，薄板运动微分方程如式(2)所示

(2)

本文仅研究矩形薄板的自由振动，薄板上并不受力，结构控制方程表示为

D▽4w(x,y)-ρhω2w(x,y)=0

(3)

式中:D为薄板的弯曲刚度，表达式为D=Eh3/(12(1-μ2));E为薄板的杨氏模量;μ为薄板的泊松比；ρ为薄板的密度;h为薄板的厚度;w(x,y)为“1.2”中的位移容许函数。

1.4 求解过程

矩形薄板弯曲产生的应变能如式(4)所示

(4)

矩形薄板边界弹簧所储存的弹性势能如式(5)所示

(5)

矩形薄板系统的动能T表示为

(6)

综上，将式(4)到式(6)代入可得矩形薄板拉格朗日函数为

L=Vp+Vs-T

(7)

由Hamilton原理可知

(8)

将式(7)代入式(8)求解可得矩阵方程如下

([K]-ω2[M]){B}=0

(9)

式中:[K]=[Kp]+[Ks];[Kp]为矩形薄板应变势能刚度矩阵;[Ks]为矩形薄板边界弹簧势能刚度矩阵;M为矩形薄板质量矩阵;B为未知的傅里叶级数矩阵；式(9)中各项可参照附录，对式(9)求解可得到矩形薄板的固有频率及振型。

2 方法有效性验证

为验证本文方法有效性，以下计算中矩形薄板材料参数参照文献[10]设置如下：密度ρ=7 850 kg/m3，泊松比μ=0.3，杨氏模量E=2.1×1011Pa；为了便于与已有文献[10]及有限元计算结果进行对比，对计算结果进行无量纲化处理，无量纲化公式为:Ω=ωa2(ρh/D)1/2。有限元计算通过ABAQUS实现，材料参数与文献[10]一致，单元类型为S4R，单元数为105。为便于表示，下文中边界条件：C为刚固，S为简支，F为自由，E为弹性。

2.1 弹簧刚度对收敛性影响分析

由上述分析可知，矩形薄板边界条件通过四条边上均布的横向位移约束弹簧与旋转约束弹簧来模拟；现讨论弹簧刚度值对矩形薄板固有频率的影响，矩形薄板边长a=1 m，宽b=1 m，厚度h=0.002 m。简支边界条件下逐渐增大旋转约束弹簧刚度，矩形薄板无量纲频率参数变化,如图4所示。旋转约束弹簧刚度大于1010N·m/rad后，矩形薄板无量纲频率参数基本不变；旋转约束弹簧刚度设为1010N·m/rad时，矩形薄板无量纲频率参数随横向位移约束弹簧刚度变化,如图4所示。

图4 矩形薄板无量纲频率参数随旋转约束弹簧刚度值变化关系

图5 矩形薄板无量纲频率参数随横向位移弹簧刚度值变化关系

由图4及图5可知，当旋转约束弹簧刚度与横向位移约束弹簧刚度大于1010时，矩形薄板无量纲频率参数已收敛，可作为刚性边界条件；约束弹簧刚度均设置为0时，可视为自由边界条件；横向位移约束弹簧与旋转约束弹簧刚度值分别取102～107N/m与102～105N·m/rad的某一值时，该边界条件介于自由边界条件与刚性固定边界条件之间，可视为某种弹性边界条件。综上所述，参照文献[11]本文任意边界弹簧刚度取值如表1所示。

表1 任意边界弹簧刚度值

2.2 截断数对收敛性影响分析

截断值M、N取值不同，会对计算精度与效率产生较大影响；因而有必要对截断值与本文计算方法收敛性之间关系进行探究。以S-S-S-S(四边简支)边界条件下的矩形薄板为研究对象，矩形薄板边长a=1 m，宽b=1 m，厚度h=0.002 m。简支边界条件中横向位移约束弹簧刚度k为1010N/m，旋转约束弹簧刚度K设置为0。

表2为截断值M=N取值从5变为14时，该矩形薄板的前5阶无量纲频率参数，此外有限元分析结果与参考文献[10]结果也列在表2中作为参考。由表2可知,截断数M=N=13后，算例中矩形薄板固有频率已基本不变，本文计算方法已收敛，在后面的计算中将选取截断数M=N=13。

表2不同的截断值M、N下S-S-S-S矩形薄板结构无量纲频率参数Ω

Tab.2 Dimensionless frequency parameters of S-S-S-S boundary condition with different truncated numbers

M=N模态阶次12345521.67251.30651.31080.919100.676621.67151.30551.30880.918100.674721.67151.30351.30580.918100.667821.67151.30251.30380.917100.666921.67051.30251.30380.917100.6641021.67051.30251.30280.917100.6631121.67051.30151.30180.917100.6621221.67051.30151.30180.917100.6611321.67051.30151.30180.917100.6611421.67051.30151.30180.917100.661FEM21.65051.28051.29080.890100.650文献[10]21.60051.28751.28780.886100.680

对应表2中任意给定阶次的模态频率，式(9)均可求解得到相对应的特征向量；特征向量对应该矩形板在给定阶次模态频率所对应振型的傅里叶展开系数，将系数代入式(1)所表征的位移容许函数，即可得到与给定阶次模态频率相对应的模态阵型。图6给出了S-S-S-S边界条件下采用本文计算方法得到的矩形薄板部分模态振型图。使用有限元分析软件(ABAQUS)计算，单元类型为S4R，单元数为105时得到的部分模态振型图如图7所示。

(a) 第1阶模态振型

(b) 第4阶模态振型

(a) 第1阶模态振型

(b) 第4阶模态振型

2.3 典型边界条件下矩形薄板振动特性计算

表3～表5列出S-S-F-F、S-S-S-S、C-F-F-F三种典型边界条件下，厚度h=0.002 m，不同长宽比矩形薄板使用本文方法与有限元法计算得到的无量纲频率参数及文献解。

表3 S-S-F-F边界条件下不同长宽比矩形薄板无量纲频率参数Ω

表4 S-S-S-S边界条件下不同长宽比矩形无量纲频率参数Ω

表5 C-F-F-F边界条件下不同长宽比矩形无量纲频率参数Ω

表3～表5中采用本文方法得到的解与有限元及文献解吻合良好。此外，图6～图7中采用本文方法与有限元得到的振型图趋势相仿；以上对比说明，本文方法准确可靠。此外，基于本文方法MATLAB求解时间仅需3 s(PC、3.4 GHz)，而有限元求解超过60 s；本方法对不同尺寸、边界条件的矩形薄板求解，仅需修改参数无需重新建模。

3 弹性组合边界条件矩形薄板自由振动特性分析

基于上述研究内容，对任意边界条件下不同长宽比，厚h=0.002 m的矩形薄板自由振动特性进行分析，旨在为后续研究及工程应用提供参考。从表6及表7可知，边界条件对矩形薄板无量纲频率参数有显著影响。固支边界条件与弹性边界条件组合中，随着固支边条界的增多，矩形薄板无量纲频率参数呈增大趋势。简支及自由边界条件与弹性边界条件组合中，随着弹性边条界的增多，矩形薄板无量纲频率参数呈增大趋势。

表6 不同边界矩形薄板自由振动无量纲频率参数Ω

Tab.6 Dimensionless frequency parameters of rectangular thin plates in different boundary conditions

长宽比r模态阶次边界条件C-C-C-EC-C-E-EC-E-E-EE-E-E-ES-S-S-E133.75432.18930.45729.16922.182259.26352.22550.41749.19248.237369.74767.68555.33349.19648.6041491.67877.00070.60564.53271.725596.48983.14374.25970.90382.6646127.065108.36290.44774.59995.727156.41453.15952.53651.98340.135287.37883.08178.90175.67664.5893118.712103.699102.989102.397100.2241.54141.000127.412119.885110.025109.3605144.928135.956122.755119.418121.1156193.532159.885156.656146.593161.383191.25485.86685.53385.24365.6622117.624110.700108.748107.02688.7293166.273157.983150.609145.346131.27224204.100178.161177.750177.380173.2815226.553197.675195.366193.435192.9486237.579228.259208.456194.921194.4181137.334129.299129.079128.89998.6302160.418150.510149.353148.317120.7813204.113192.221188.439185.107161.4212.54270.414256.904246.754238.450222.2555314.536274.658274.361274.119267.2846335.077291.831290.393289.136286.3271194.204182.963182.807182.686138.9972215.099201.769201.019200.286160.5283254.715238.789236.491234.081199.64334315.980297.729291.730286.417258.3575400.076379.768366.959354.642337.5436449.806392.928392.731392.503382.2521252.189246.701246.619246.488186.7062274.489263.862263.312262.773207.8513281.078297.287295.624294.112245.7633.54351.198350.970347.072343.349302.6085374.110427.171418.886411.329379.5556502.399526.706510.274495.878477.109

4 结 论

通过将薄板振动的位移函数表示成二维傅里叶余弦级数和辅助傅里叶级数的线性组合，列出矩形薄板的拉格朗日方程；然后通过Hamilton原理对拉格朗日方程求解得到矩形薄板自由振动无量纲频率参数。通过与文献及有限元计算结果对比，验证了本文方法的有效性。通过本文研究，可得如下主要结论：

表7 不同边界矩形薄板自由振动无量纲频率参数Ω

Tab.7 Dimensionless frequency parameters of rectangular thin plates in different boundary conditions

长宽比r模态阶次边界条件S-S-E-ES-E-E-EF-F-F-EF-F-E-EF-E-E-E124.71426.7568.18519.54320.739247.56148.3588.44522.17331.896348.00948.74418.95135.13643.9581467.19965.54027.87042.92653.154573.37372.93328.06845.91055.835693.92480.17846.19958.42868.620148.65750.0717.74744.08945.157267.42171.59914.06546.93255.3413101.073101.68631.41458.62080.8311.54108.026109.20045.88986.97597.6975116.744118.13653.86396.713106.5496151.198148.48868.45899.770120.851182.88483.89421.08278.48679.494298.861102.69842.93981.36688.9553134.042139.90063.96992.678112.25024176.394176.84973.902117.583153.3335190.238191.81291.444162.818173.0056192.035193.679113.479172.061181.5251127.117127.89129.437122.742123.6822141.459144.62447.347125.594132.6833172.362178.58883.824136.761154.0172.54224.676231.849111.930159.639192.4285273.338273.708136.638200.164249.5906286.278287.633147.401262.685269.8921181.266181.8959.283176.935177.5812194.628197.23339.298179.648186.3753222.427228.08060.702191.268205.98734269.682278.07695.207212.410241.4575339.301347.347158.304251.670296.2256391.881392.188191.591306.672369.0811245.340245.83550.661240.835241.7252258.040260.22974.098243.526250.2443283.652288.590108.545254.648268.9653.54326.761334.866149.489275.451302.1475390.935401.180232.355310.261353.0296477.963487.258247.703362.905423.776

(1) 本文算法收敛性与边界弹簧刚度及截断数相关。旋转约束弹簧刚度K>1010N·m/rad，横向位移约束弹簧刚度k>1010N/m时，矩形薄板无量纲频率参数已收敛，可作为刚性边界条件。横向位移约束弹簧刚度k为102N/m～107N/m中某一值时，可将其视为某种弹性边界条件。

(2) 典型边界条件下本文计算结果与文献及有限元解吻合良好，方法准确可靠。

(3) 本文方法与有限元相比具有更高的求解效率；对不同尺寸、边界条件的矩形薄板求解，仅需修改参数，无需重新建模。

(4) 固支边界条件与弹性边界条件组合中，随着固支边条界范围增大，矩形薄板无量纲频率参数呈增大趋势。简支及自由边界条件与弹性边界条件组合中，随着弹性边条界范围增大，矩形薄板无量纲频率参数呈增大趋势。

[cos(λm1x)cos(λn1y)]dxdy

[cos(λm1x)sin(λn1y)]dxdy

MBA=MABT，MCA=MACT

2(1-μ)λmλnλm1λn1sin(λmx)sin(λny)sin(λm1x)sin(λn1y)]dxdy

2(1-μ)λmλnλm1λn1sin(λmx)sin(λny)sin(λm1x)cos(λn1y)]dxdy