虚拟弹簧边界Mindlin板的自由振动和屈曲分析*

2022-02-18单逸峰

起重运输机械 2022年1期

单逸峰向维

西南交通大学机械工程学院成都 610031

0 引言

起重机的箱形梁结构，由上下翼缘板、主副腹板以及横纵向加劲肋焊接组合而成[1]。准确求得板的自振频率和屈曲载荷对起重机的安全工作至关重要。以有限元法(以下简称FEM)为代表的数值解法易于求得高精度近似解，然而，FEM需要把研究对象划分成大量微小单元，使计算过程耗时更多。Xing Y等[2]提出的名为微分求积有限单元法(以下简称DQFEM)的全局数值方法在计算二维静态问题时表现出高精度和快速收敛性。

在应用数值算法时，Wang Q等[3，4]在对壳体的自由振动分析中多次应用到虚拟弹簧边界。Kim K等[5]在分析功能梯度矩形板的自由振动时，使用虚拟弹簧边界引入边界条件，所得数值结果准确有效。Du Y等[6]在分析球型盖的动态行为时指出，虚拟弹簧边界易于模拟不同的边界条件。由于采用虚拟弹簧边界不减少整体矩阵的阶数，故完整的模态向量也可直接求得。Guan X等[7]基于经典板理论(以下简称CPT)，将虚拟弹簧边界和DQFEM共同应用于求解薄板的自由振动问题。但是，CPT未考虑横向剪切变形，使求得的自然频率偏大[8]。

针对板的自由振动和特征屈曲问题，基于考虑了横向剪切变形的Mindlin板理论，采用DQFEM先后完成虚拟弹簧边界Mindlin板的自由振动和屈曲分析。数值结果表明，虚拟弹簧边界与DQFEM相结合的数值方法所需结点数更少，求得的自振频率精度更高、屈曲因子偏向于安全保守，且虚拟弹簧边界对DQFEM的计算效率有所提高。

1 虚拟弹簧边界的Mindlin板

Mindlin板理论能够较好地协调计算精度和效率，因而广泛应用于模拟中厚板的力学行为。其位移场为

式中：w为沿z轴的移动，φx为xOz平面内的转动，φy为yOz平面内的转动。

在板的边界配置3种虚拟弹簧，用以约束边界的每个广义自由度，如图1所示。弹簧标记kti、krxi和kryi的含义由表1列出。

图1 边界受虚拟弹簧约束的Mindlin板

表1 弹簧标记及其含义

Gauss-Lobatto结点已被证明是能有效保证计算精度的网格划分方法[9]，将Mindlin板分别沿x轴和y轴划分出M个和N个Gauss-Lobatto结点。虚拟弹簧边界离散在各个Gauss-Lobatto结点上，为展示清晰，图2仅完整表示边2和边3的全部边界弹簧。

图2 划分Gauss-Lobatto结点的虚拟弹簧边界Mindlin矩形板

1.1 自由振动特征方程

要推导自由振动特征方程，首先要列出系统的能量。Mindlin板的动能T和应变能Us分别为

式中：D为弯曲刚度；G为剪切模量；E为弹性模量，E=10 920 Pa；υ为泊松比，υ=0.3；ρ为密度，ρ=1；A为板中面的面积；I为转动惯量；к为剪切修正系数。

边界上，线弹簧的弹性势能Vt，扭转弹簧在xOz面内和yOz面内的弹性势能Vrx、Vry分别为

式中：i为板边编号，li为第i条边的长度。

由式(2)～式(4)可得Mindlin板的Lagrange函数为

势能U中含有充当边界条件的虚拟边界弹簧的弹性势能，由此推导出的运动方程中便已施加过边界条件，而无需修改整体矩阵，不同于传统数值方法。

基于DQFEM离散Lagrange函数，由Hamilton变分原理，令δL=0，可得Mindlin板的振动特征方程为

式中：K为系统的刚度矩阵，M为质量矩阵。

w是由Mindlin板上所有结点的位移所组成的列向量，即

刚度矩阵K中的Ks来自于板的应变能，Ksprings来自于虚拟弹簧边界的弹性势能。

式中：υ1=(1-υ)/2；υs=6κ(1-υ)/h2；分别为二元函数沿x轴和y轴方向的一阶权系数矩阵；C是二维Gauss-Lobatto积分系数矩阵；D1、D2、D3由式(12)求得，即

Ksprings中的各个分量分别为

为便于同前人文献的数据对比，将式(6)得出的自由振动频率ω代入式(16)，求得无量纲频率参数为

1.2 屈曲特征方程

虚拟弹簧边界Mindlin矩形板在图3所示的均匀初始应力作用下，系统的总能量为

图3 有面内应力作用的矩形板

2 数值方法参数的确定

为使数值方法得到准确的数值结果，需要配置足够的结点数以使数值解收敛；要使虚拟弹簧边界代替经典边界条件，需要确定各个虚拟弹簧的刚度取值。

2.1 结点数

为确定使数值解收敛的最少结点数，针对CCCC和SSSS边界，厚度比h/b为0.01和0.1的4种Mindlin方板，图4为结点数对第1、3、5阶频率参数的影响情况。

为便于计算，取长和宽方向结点数相同，即M=N。当模拟CCCC边界时，令所有弹簧刚度为1010N/m，κ=0.860 1；当模拟SSSS边界时，令kx1=kx2=ky1=ky2=0，其他弹簧刚度为1010N/m，κ=5/6。图中虚线分别为第1、3、5阶自振频率的收敛值，箭头指向前6阶频率收敛时的最少单边结点数。

由图4可得，使用虚拟弹簧边界代替经典边界条件，DQFEM仍有良好的收敛性；为保证前6阶无量纲频率参数收敛，应配置的结点数为12×12。

图4 结点数对频率参数的影响

2.2 虚拟边界弹簧的刚度取值

为分析不同类型弹簧在模拟经典边界条件时应取的刚度值，将弹簧分为3类：法向扭转弹簧krxi(i=1，2)与kryi(i=1，2)；切向扭转弹簧 krxi(i=3，4)与 kryi(i=3，4)；线弹簧kti(i=1，2，3，4)。针对线弹簧和扭转弹簧，引入无量纲刚度

用Sφ1和Sφ2分别表示法向和切向扭转弹簧的无量纲刚度。

图5为3类弹簧的刚度对厚度比h/b为0.01和0.1的Mindlin方板的第1、3、5阶无量纲频率参数的影响情况。由图5可得，不论哪种边界弹簧，当其刚度大于某一定值，频率参数不再随弹簧刚度的增大而变化。这可以解释为弹簧的刚度太大，以至于板边界被近似严格约束。对于法向、切向扭转弹簧以及线弹簧，实现这种严格约束的无量纲刚度分别为107N/m、108N/m和1010N/m。

图5 频率参数随3类弹簧刚度的变化情况

3 自由振动分析数值结果

针对Mindlin板自由振动问题，本小节首先验证数值结果的准确性，而后讨论虚拟弹簧边界对DQFEM计算效率的影响。采用12×12的结点划分Mindlin板，边界弹簧的刚度设置为：固支边(C)上，Sw=Sφ1=Sφ2=1010N/m；简支边(S)上，Sw=Sφ2=1010，Sφ1=0；自由边(F)上，Sw=Sφ1=Sφ2=0。

3.1 准确性

为验证本文数值解的准确性，以SSSS、SCSC和CCCC边界条件的厚度比h/b为0.01和0.1的Mindlin板为对象，计算频率参数。表2～表4将本文数值解同Hinton[10]的精确解、Ferreira[11]的有限元解、Dawe D J等[12]的Rayleygh-Ritz解进行对比。Ferreira的角标1和2分别为15×15和20×20的Q4有限元法。

由表2和表3可知，在SSSS和SCSC边界条件下得出的前5阶无量纲频率比FEM求得的数值解更接近精确解。由表4可知，在CCCC边界条件下求得的前5阶无量纲频率比FEM求得的数值解更接近Rayleygh-Ritz法的数值结果。

表2 SSSS边界，κ=0.833，Mindlin方板前5阶无量纲频率ω

表3 SCSC边界，κ=0.822，Mindlin方板前5阶无量纲频

h/b 计算方法频率阶次1 2 3 4 5 Present 0.141 1 0.266 8 0.337 6 0.460 4 0.497 7 Hinton 0.141 1 0.266 8 0.337 7 0.460 8 0.497 9 Ferreira1 0.142 4 0.271 0 0.348 4 0.472 2 0.519 1 0.01 Present 1.300 1 2.393 9 2.884 5 3.839 1 4.231 3 Hinton 1.302 2.398 2.888 3.852 4.237 Ferreira1 1.294 0 2.397 1 2.929 0 3.839 4 4.347 5 0.1

表4 CCCC边界，κ=0.860 1，Mindlin方板前5阶无量纲频

h/b 计算方法频率阶次1 2 3 4 5 Present 0.175 4 0.357 4 0.357 4 0.526 5 0.639 9 Ferreira2 0.175 0 0.363 5 0.363 5 0.535 8 0.663 4 Dawe 0.175 4 0.357 6 0.357 6 0.527 4 0.640 2 0.01 Present 1.591 0 3.038 9 3.038 9 4.262 5 5.024 7 Ferreira2 1.595 5 3.066 2 3.066 2 4.292 4 5.123 2 Dawe 1.594 0 3.039 0 3.039 0 4.265 0 5.035 0 0.1

由此，使用虚拟弹簧边界替代经典边界条件后，采用DQFEM求解Mindlin板的自由振动问题的准确性得以验证。值得注意的是，DQFEM仅需划分12×12的结点数，而求得的频率参数比20×20的Q4有限元法所得结果更加准确。

3.2 计算效率

由于边界条件施加方法的改变，DQFEM的计算效率可能改变。表5以厚度比h/b=0.1的Mindlin方板为研究对象，分别统计了对整体矩阵消元和使用虚拟弹簧施加SSSS和CCCC边界条件时DQFEM求得前6阶无量纲频率及完整的模态向量所花费的时间。

由表5可知，使用虚拟弹簧边界施加边界条件有助于减少DQFEM的计算时间。对于SSSS边界，消元法和虚拟弹簧边界的10个算例的平均时间分别为0.422 0 s和0.412 2 s，虚拟弹簧边界使DQFEM的平均计算时间缩短了2.32%；对于CCCC边界，消元法和虚拟弹簧边界的10个算例的平均时间分别为0.314 9 s和0.305 2 s，虚拟弹簧边界使DQFEM的平均计算时间缩短了3.08%。

表5 消元法和虚拟弹簧边界DQFEM的计算时间 s

4 屈曲分析数值结果

表6和表7列出SSSS、SCSF和SFSF边界厚度比h/b为0.05和0.1的Mindlin方板，在分别受x轴方向和受双轴向均匀载荷时的临界屈曲因子n*，同时列出Mizusawa T[13]的数值解、Hashemi S等[14]的精确解。剪切修正因子κ取π2/12。数值结果与Hashemi S等精确解的相对误差由式(26)求得，即

对比表6和表7数据可知，将DQFEM和虚拟弹簧边界共同用于Mindlin板的屈曲分析时，得出的屈曲因子与Mizusawa通过样条方法得出的数值解几乎完全吻合，而略小于Hashemi的精确解，但最大相对误差不超过1.5%。这说明使用虚拟弹簧边界模拟经典边界条件，并使用DQFEM分析Mindlin板的特征屈曲问题，具有一定准确性，求得的临界屈曲因子偏向于安全保守。