基于二次特征值法的矩形薄板的动力稳定性研究*

钟子林，刘爱荣，卢汉文，黄友钦

(广州大学-淡江大学工程结构灾害与控制联合研究中心，广东 广州 510006)

基于二次特征值法的矩形薄板的动力稳定性研究*

钟子林，刘爱荣，卢汉文，黄友钦

(广州大学-淡江大学工程结构灾害与控制联合研究中心，广东 广州 510006)

基于Von-Karman薄板大挠度理论，利用Galerkin法得到面内周期荷载作用下四边简支矩形薄板的二阶常微Mathieu-Hill型参数振动方程；运用二次特征值法分别求出矩形薄板线性参数振动方程周期为2T和T时的主要与次要动力不稳定域，并用有限元数值分析方法验证了二次特征值法的精确性，同时定性地分析了主要参数共振下非线性弹性对系统定态振幅的影响。分析结果表明：① 当激发力频率近薄板两倍自振频率时，薄板发生强烈的横向参数共振；② 二次特征值法可精确计算矩形薄板发生动力不稳定时对应的频率和激发系数；③ 随着薄板振幅的增长，非线性的存在抑制了定态振动幅值的无限增长，牵引系统向大频率方向振动，导致振幅稳定增加或迅速增大的复杂振动状态。

矩形薄板；二次特征值法；参数共振；动力不稳定域；周期荷载

板壳结构在航空、建筑、桥梁、机械等结构中有着广泛的应用，因此对板壳结构力学行为的研究受到学者的普遍关注。结构在压力作用下处于平衡状态，当压力值达到临界荷载时，微小的外部干扰将会使其发生失稳现象[1]。

例如，2004年，温州大学体育场看台的板壳结构在台风“云娜”的袭击下发生整体破坏的事故；2005年，宁波北仑体艺中心板壳结构屋顶在台风“麦莎”的作用下出现撕毁破坏事故[2]。所以揭示板结构在地震、风荷或车辆等动力荷载作用下的稳定性能亦是当前研究的重要课题。

关于板结构动力稳定性能的研究已有广泛涉及[3-13]。例如，鲍洛金[8]在其著作中探讨了面内受周期荷载作用的板在一定条件下会发生剧烈的参数共振问题；Konishi等[9]基于薄板小挠度理论，研究了不同边界条件下面内受周期变化荷载作用的矩形薄板的动力稳定性；Ganapathi等[10]利用Newmark积分和迭代法研究了面内受周期荷载作用的各向同性和各向异性弹性板的非线性动力不稳定特性；Dey等[11]研究了面内周期荷载作用下复合材料层合板的动力稳定特性，并利用Bolotin法获得不同设计参数薄板的主要和次要不稳定域；Amabili[12]分别对不同边界条件下受到面内周期荷载的矩形薄板的非线性振动进行理论与实验的研究分析；傅衣铭等[13]运用谐波增量法研究了四边简支矩形板在纵横谐载共同作用下非线性动力稳定性问题；彭凡[14]采用Boltzman积分型本构关系，运用谐波增量法获得粘弹性板在不同设计参数下的主要不稳定域。

但目前关于板的动力稳定性的研究主要是采用近似法来获得动力系统的不稳定域，当动力系统的激发系数较大以及迭代次数较多时，此方法的计算结果误差较大[14]；此外，许多文献中并未考虑阻尼对矩形薄板动力稳定性的影响，事实上，研究表明阻尼对结构动力稳定域的结果影响较大，尤其体现在次要不稳定域的结果中[15]。本文针对前人研究的不足，采用二次特征值法精确求出矩形薄板线性参数振动方程周期为2T和T时的主要与次要不稳定域，并用有限元数值法验证了二次特征值法的精确性，同时定性地分析了主要参数共振下非线性弹性对系统定态振幅的影响。

1 动力系统模型

为探讨薄板在周期荷载作用下的动力稳定性能，本文选择四边简支、长为a，宽为b，厚度为h的各向同性矩形薄板为研究对象(图1)。并假设：矩形薄板的加载边在面内可以移动和转动，非加载边在面内不可动；加载边受周期简谐荷载Px=Po+Ptcos(t)作用；板的四边无剪切荷载，忽略板的面内惯性力的作用。

基于Von-Karman薄板大挠度理论，考虑面外惯性作用以及阻尼的影响，得到板在外部激励作用下系统的动力方程组为［16］:

(1)

式中，ω为板的横向位移函数；φ为板应力函数；ω为板的密度；ε为阻尼比，ε=ΩΔ/2π，Δ、Ω分别为承受纵向力不变分量的板的固有振动阻尼缩减率和固有振动频率；D=Eh3/12(1-ν2)，E和ν为材料的弹性模量和泊松比。

图1 动力系统模型Fig.1 Dynamic system model

板的非线性振动的边界条件包括侧向边界条件和面内边界条件。

侧向边界条件(简支)

(2)

面内边界条件

(3)

(4)

由薄板振动理论可知，板在振动时沿板长、宽方向形成屈曲半波，根据屈曲半波数选择其横向振动挠度函数［17］:

(5)

本文采用一阶模态Galerkin逼近进行分析，选取横向位移函数展开式中的第一项，即

(6)

式中，ω11和f11分别为基本模态的横向位移和振动幅值，下标表示板振动时沿板长、宽方向各形成一个屈曲半波，为简单起见，在下文表述中省略其下标。

将式(5)和式(6)代入方程组(1)的第二式，可得

(7)

根据应力边界条件(3)和(4)求得满足式(6)的应力函数

(8)

将应力函数φ代入式(1)的第一式，并进行伽辽金积分，得

(9)

于是，方程(8)可表示成以下含非线性项的Mathieu-Hill方程

(10)

(11)

设f具有周期为2T的周期解展开的傅里叶级数形式为［8］:

(12)

同理，可设f具有周期为T的周期解展开的傅里叶级数形式为式(14)[8]。

为使方程组(13)、 (15)具有非零解，则其系数行列式应等于零，从其系数行列式中可知周期为2T和周期为T时Mathieu-Hill 方程的激发系数μ和外荷载频率与2倍结构固有频率比θ/2Ω之间存在一定的关系，通过进一步的分析可得到周期为2T和T时的动力不稳定域[18]。

鲍洛金[7]应用连分式的数值方法，对方程组(13)和(15)的系数行列式进行逐次迭代求得周期为2T和T时的动力不稳定边界。但是，这种方法在求解过程中需要较大的计算量，每次迭代只能得出一个不稳定区域的边界线，且在激发系数较大时，误差较大。为较为准确地得出振动系统的不稳定域，本文采用二次特征值法求解周期为2T和T时振动系统的动力不稳定边界。

2 二次特征值法求解方程组

观察周期为2T时振动系统的方程组，对其系数行列式进行分析。

(13)

(14)

(15)

(16)

(17)

(18)

方程组(14)可利用上面的矩阵表示为：

(19)

(20)

运用特征值法便可求得周期为2T时振动系统的不稳定域边界，同时观察到周期为T时的系统方程不能简单地按上述方法建立求解特征值的相应方程。由于方程组(15)的系数矩阵右下部分B0的奇异性大，影响矩阵的特征值求解，可对该矩阵进行一次特征值计算，然后将求解得到的特征值重新加入到总矩阵中并运用上述求解周期为2T时动力不稳定边界的方法求解得振动系统周期为T时的动力不稳定边界。

对矩阵B0进行一次特征值求解，由于该矩阵对角元素含有实数1，求解时会出现无穷大的特征值，将无穷大的特征值舍去，得到含有特征值解的矩阵B：其中kn为特征值，将化简后的矩阵B重新加入到总矩阵中，运用二次特征值法求解系数行列式(23)便可得到周期为T时的动力不稳定边界。

(21)

(22)

(23)

3 动力不稳定性分析

基于以上理论分析，利用MATLAB编制程序，计算图1所示薄板模型的动力不稳定域。本文选取如下计算参数，即四边简支板薄板的长a=1 m，宽b=0.8 m，厚度h=5 mm;材质为铝材：弹性模量E=6.9×1010Pa，泊松比v=0.33，密度ρ=2 700 kg/m3；阻尼比ε=0.01。有限元模型见图2。

图3、4为矩形薄板一对边在面内受到周期荷载作用并考虑阻尼影响，周期解分别为2T(θ=2Ω/k，k=1，3，5，7，9)和T(θ=2Ω/k，k=2，4，6，8，10)时前五阶动力不稳定域(阴影部分)，其中横坐标表示外荷载作用频率θ与2倍薄板固有频率Ω比值的变化，纵坐标表示外荷载pt与薄板屈曲荷载Pcr比值的变化，虚线是利用近似迭代法得到的主要不稳定域解[7]。由图3和图4可知，面内受周期荷载作用的矩形薄板除了会发生周期为T的强迫共振外，还会发生周期为2T的强烈横向参数共振；此外，主要不稳定域和次要不稳定域的激发系数起始值大于零，这是由于考虑阻尼影响导致的，通常，若不考虑阻尼影响不稳定域激发系数的起始值为零[15]；此外，从图3和图4中可以看出，在阻尼相同时，主要不稳定域的起始值与激发系数原点的距离远小于次要不稳定与激发系数原点的距离，即阻尼对次要不稳定域的影响较主要不稳定域明显；图3和图4亦表明随着频率比的降低，不稳定域逐渐减小，当阻尼一定时，较小的激发系数就可使薄板出现主要不稳定域，而不会出现次要不稳定域，即只有当激发系数较大时才会出现次要不稳定域；当激发系数较大时，由二次特征值法求解得到的解与近似法求解结果差异较大，在求解过程中迭代法用先前计算得到的结果代入到求解下一个不稳定区域的计算当中，使得误差逐级累加并且所需计算量较大，而二次特征值法的精度随着所选取的谐波数的增加而增加。

为了验证二次特征值法的正确性，利用有限元软件ANSYS对该简支薄板模型进行瞬态仿真模拟分析。根据图1所示薄板模型的四边简支边界条件，ANSYS建模时可约束四边面外的自由度Uz，并分别约束平行于X轴的对称轴面内平动自由度Uy和平行于Y轴的对称轴面内平动自由度Ux。结合ANSYS在求解过程中的计算量与精确度，将薄板单元网格的大小设为25mm。

考虑到在次要不稳定域中的薄板振动幅值较小，不能凸显薄板在动力荷载作用下的不稳定性，故在图3、4中的主要不稳定域中选取点A(0.65,0.6)、B(1,26,0.6)、C(0.36,0.6)、D(0.52,0.6)来验证本文方法的准确性。通过计算得到该矩形薄板的自振频率为3.11Hz，以上所

选取点对应的外荷载作用频率分别为θA=4.04 Hz，θB=7.84 Hz，θC=2.24 Hz，θD=3.23 Hz，根据各点对应的激发系数取周期外荷载幅值PtA=PtB=PtC=PtD=34 N，令外荷载恒值P0等于零，分别编写相应的ANSYS命令流，对受到简谐荷载作用下的矩形薄板进行定频模拟分析，得到点A、B、C、D的时域图5-图8。

图5 点B的瞬态响应Fig.5 Transient response of point B

图6 点A的瞬态响应Fig.6 Transient response of point A

图7 点D的瞬态响应Fig.7 Transient response of point D

图8 点C的瞬态响应Fig.8 Transient response of point C

由图3和图4可知：点B和点D位于本文方法计算得到的不稳定域之内而位于近似法求解得到的不稳定域之外，结合此两点的响应图可知，在与点B对应的激发系数和频率比相同的动力荷载作用下，薄板处于不稳定的跳动状态，当该动力荷载持续作用于薄板时将会发生动力失稳[8]，同样，D点也是如此；同时观察到点A和点C位于本文方法计算得到的不稳定域之外而位于近似法求解得到的不稳定域之内，结合此两点的响应图可知，在与点A对应的激发系数和频率比相同的动力作用下，薄板发生相对于板厚为微小振幅的稳定状态，没有发生失稳现象，C点亦然。综上可知，二次特征值法比近似法更能准确地反映动力系统发生失稳时激发系数与外荷载作用频率的关系。

(24)

式中，A=a2+b2，表示定态振动的振幅；Δ、Ω同式(1)，γ同式(9)。通过计算可得：λ=6.22×108；Ω=2π×3.11=19.54;μ=pt/2Pcr=0.6，取阻尼比为0.01，将相关系数代入式(24)中，可得含定态振幅与频率比的非线性变化趋势。

在周期荷载作用下，由线性系统所求的振幅是无限增长的，但这与实际的情况不相符合。由点B的响应图可知，动力系统发生参数共振时薄板的最大的位移为1.225mm，根据点B的定态振幅与频率比的关系曲线(图9)，可知点B位于稳定解的上方；同理可知处于平稳状态的点A则位于稳定解的下方，表明当薄板发生参数共振时动力系统表现出了强烈的非线性特征，有效抑制了定态振动幅值的无限增长，牵引系统向大频率方向振动，导致系统处于振幅稳定增加或迅速增大的复杂振动状态。

图9 非线性系统下定态振幅与频率比的关系Fig.9 The relation between steady amplitude and frequency ratio in nonlinear system

4 结 论

本文基于Von-Karman薄板大挠度理论，利用Galerkin法得到面内周期荷载作用下四边简支矩形薄板的二阶常微分Mathieu-Hill型参数振动方程，运用二次特征值法和ANSYS有限元数值方法分析了矩形薄板的动力不稳定性，可以得出以下结论：

1)二次特征值法可精确计算矩形薄板发生动力失稳时的主要不稳定域和次要不稳定域，较其他计算方法精度高。

2)当外荷载作用的频率与薄板的固有频率的比值在两倍附近时，面内受到周期荷载作用的矩形薄板除了会发生强迫共振外，还可能引发强烈的横向参数共振现象。

3)随着频率比的降低，矩形薄板的不稳定域逐渐减小；阻尼对次要不稳定域的影响比较明显，当阻尼一定时，较小的激发系数就可使薄板出现主要不稳定域，而不会出现次要不稳定域，即只有当激发系数较大时才会出现次要不稳定域。

4)在薄板振动过程中，非线性因素随着振动系统振幅的增长开始发生作用，其主要的作用是抑制定态振动幅值的无限增长，牵引系统向大频率方向振动，逐渐远离发生参数共振的频率段，导致系统处于振幅逐渐减小或增大的复杂振动状态。

[1] 禹奇才，刘爱荣，肖才涛，等. 基于尖点突变理论的平行组拼双肋拱侧倾失稳临界荷载计算新方法[J]. 中山大学学报 (自然科学版)，2014，53 (2)：6-11.

[2] 沈世钊, 武岳. 大跨空间结构抗风研究新进展与展望[C]∥ 第十三届全国结构风工程学术会议.大连，2007.

[3]BIRMANV.Dynamicstabilityofunsymmetricallylaminatedrectangularplates[J].MechanicsResearchCommunications, 1985, 12(2):81-86.

[4]LEEHP,NGTY.Dynamicstabilityofamovingrectangularplatesubjecttoin-planeaccelerationandforceperturbations[J].AppliedAcoustics, 1995, 45(94):47-59.

[5]SINHASK.Commentson``dynamicstabilityofarectangularplatesubjectedtodistributedin-planedynamicforce''[J].JournalofSound&Vibration, 1989, 128(3):505-507.

[6]CHIACY.Largeamplitudevibrationsoflaminatedrectangularplates[J].FibreScience&Technology, 1982, 17(2):123-131.

[7]RAOGV,RAJUIS,RAJUKK,etal.Afiniteelementformulationforlargeamplitudeflexuralvibrationsofthinrectangularplates[J].Computers&Structures, 1976, 6(6):163-167.

[8] 鲍洛金. 弹性体系的动力稳定性[M]. 北京：高等教育出版社, 1960.

[9]TAKAHASHIK,KONISHIY.Dynamicstabilityofarectangularplatesubjectedtodistributedin-planedynamicforce[J].JournalofSound&Vibration, 1988, 123(1):115-127.

[10]GANAPATHIM,PATELBP,BOISSEP,etal.Non-lineardynamicstabilitycharacteristicsofelasticplatessubjectedtoperiodicin-planeload[J].InternationalJournalofNon-LinearMechanics, 2000, 35(3):467-480.

[11]DEYP,SINGHAMK.Dynamicstabilityanalysisofcompositeskewplatessubjectedtoperiodicin-planeload[J].Thin-WalledStructures, 2006, 44(9):937-942.

[12]AMABILIM.Nonlinearvibrationsofrectangularplateswithdifferentboundaryconditions:theoryandexperiments[J].Computers&Structures, 2004, 82(s31/32):2587-2605.

[13] 傅衣铭, 张思进. 中厚矩形板的非线性动力稳定性分析[J]. 工程力学, 1995(3):97-106.

[14] 彭凡, 傅衣铭. 粘弹性板的非线性动力稳定特性分析[J]. 固体力学学报, 2004, 25(1):115-118.

[15]HUANGYQ,LUHW,FUJY,etal.Dynamicstabilityofeulerbeamsunderaxialunsteadywindforce[J].MathematicalProblemsinEngineering, 2014, 2014(2):1-12.

[16] 袁尚平, 张建武, 王庆宇. 薄板3倍超谐振动的分析与试验[J]. 上海交通大学学报, 2001, 35(7):955-961.

[17]AMABILIM.Theoryandexperimentsforlarge-amplitudevibrationsofrectangularplateswithgeometricimperfections[J].JournalofSound&Vibration, 2006, 291(3):539-565.

[18] 李康杰, 刘爱荣, 卢汉文,等. 圆弧浅拱面外动力稳定性实验研究[J]. 中山大学学报(自然科学版), 2015,54(3):30-35.

Dynamic stability research on rectangular thin plate based on the quadratic eigenvalue method

ZHONGZilin,LIUAirong,LUHanwen,HUANGYouqin

(Guangzhou University-Tamkang University Engineering Structure Disaster and Control Joint Research Center，Guangzhou University，Guangzhou 510006，China)

Based on the Von-Karman large deflection theory of thin plates, the second-order ordinary differential Mathieu-hill Parametric vibration equation was established for of simply supported rectangular thin plate under in-plane periodic loading using the Galerkin method. The quadratic eigenvalue method was adopted to obtain the main and secondary instability domains of the rectangular thin plate with the periodic solutions of 2T and T, then the accuracy of the method was verified by finite element method, meanwhile, the influence of nonlinear elasticity on stationary amplitude caused by the main Parametric resonance was qualitatively analyzed. The analytic results showed that strong transverse parametric resonance will occur when the frequency of excitation force is about twice as large as the natural vibration frequency of the thin plate. The quadratic eigenvalue method can be used to calculate accurately those parameters, such as the frequency and the excitation coefficient of the dynamic instability region of the rectangular thin plate. As the increasing of the amplitude, the infinite increase of vibration amplitude was controlled owing to the nonlinear factor which tows system to the direction of large vibration frequency. Furthermore, the nonlinear factor can bring the system into a complex vibration state in which the vibration amplitude increased stably or rapidly.

rectangular thin plate; quadratic eigenvalue method; parametric resonance; dynamic instability domain；periodic load

10.13471/j.cnki.acta.snus.2016.03.012

2015-12-30

国家自然科学基金资助项目(51578166) ;广州市羊城学者资助项目(1201541551)；广东省教育厅高层次人才资助项目

钟子林( 1990年生)，男；研究方向：新型拱桥的静、动力稳定性；通讯作者： 刘爱荣;E-mail:liu-a-r@163.com

TU

A

0529-6579(2016)03-0068-09