任意边界条件下变厚度中厚矩形板的静动态特性分析

吴子奇，王 治，史冬岩，王青山，姚熊亮

（1.哈尔滨工程大学 a.船舶工程学院；b.机电工程学院，哈尔滨150001；2.中南大学 机电工程学院，长沙410083）

矩形板结构作为主要的结构构件，已经被广泛应用于航空航天、机械、建筑、造船、桥梁等领域，在实际工程中，除了等厚度矩形板以外，变厚度板也有着广泛的应用价值，特别是在一些特殊的土建水利等方面。因此研究变厚度板的静动态特性具有较为重要的工程价值，使设计者在进行变厚度结构设计时能够了解其静动态特性。

从现有研究文献来看，在以往的研究中，对于等厚度结构的研究成果较为丰硕，其中主要的研究方法包括辛本征值展开法、无网格数值解法、离散积分法等，但是上述方法中的推导过程与求解过程都较为繁琐，边界条件较为单一，对于一般边界条件下的静动态特性无法计算；对于变厚度中厚板静动态特性的研究，从现有文献来看其主要的研究方法包括双重傅立叶级数[9]、半解析解法[10-13]、二次样条函数[14]、广义微分求积法[15]、辛几何法[16]和能量法[17]，对这些方法进行总结不难发现，大部分的方法都只局限于简单的经典边界条件，对于一般边界条件并不能进行求解，且求解过程较为繁琐。文献[18]提出了一种改进傅立叶级数方法，建立了一般边界条件下矩形薄板的静动态法分析模型，与传统的傅里叶级数相比收敛性有较大改善，且精度较高。

本文基于现有文献研究的不足，利用改进傅立叶级数研究方法，建立了一般边界条件下等厚度与变厚度矩形板结构的静动态特性分析模型。首先将结构的横向位移与旋转角度函数用改进傅立叶级数进行展开，并且对其进行能量描述。其次将函数展开的未知傅立叶级数系数作为广义向量，利用Rayleigh-Ritz方法进行求解，得到一个标准特征值问题。通过这个特征方程组，就能简便得到结构的静动态特性。本文建立的模型能够对于等厚度与变厚度中厚板结构在一般边界条件下的静动态特性进行求解。最后通过大量的数值算例对本文方法进行了验证，证明了本文方法的合理性，而且具有良好的计算精度和收敛速度。

1 矩形厚板弯曲振动理论模型

1.1 结构物理模型

本文所研究的变厚度Mindlin矩形板模型如图1所示，在直角坐标系（x，y，z）中，长为a，宽为b，厚度为h（x）的变厚度中厚矩形板。对于板结构中的任一质点，其横向振动位移函数采用w来描述，其单位为m；转角位移函数采用ψx和ψy来表示，其单位为m。板的一般边界条件可以采用沿各边均匀分布的线性位移（kil，i取值为 x、y，l的取值为0、a、b，在后面不再进行描述）和旋转约束弹簧（Kil）来模拟。由于结构在边界处还存在扭矩，因此为了更加精确地对其进行求解，在边界处还引入扭转弹簧类型（Kijl）来对其扭矩效果进行模拟，各种经典的边界条件可以通过改变弹簧刚度值的大小而简单实现。当沿边界分布的两类约束弹簧刚度值均为无穷大时，则为固支（C）边界条件。通过将两类约束弹簧刚度值设定为0时，则简单地实现自由（F）边界条件。简支（S）边界条件则可以通过将边界上线位移和旋转约束弹簧刚度值分别设置为无穷大和极小数而简便实现。对于边界上的线位移和旋转约束弹簧刚度的单位分别为N/m2和Nm/rad。板结构的边界条件采用逆指针顺序来进行描述，顺序如下：y=0、x=a、y=b、x=0。在本文计算中，无穷大的刚度值取为1×1014，而极小刚度值取为零。

1.2 结构运动关系与本构关系

根据线性小弹性变形理论，结构应变分量与结构位移变量之间的关系可以定义为：

根据胡克定律，应力与应变之间的关系定义为：

其中：

E为结构弹性杨氏模量，单位为Pa，μ为材料的泊松比，κ为剪切修正因子。

在物理模型的描述中可知，对于结构的一般边界条件模拟是依据设置在边界处的三类弹簧刚度来进行模拟，因此一般弹性边界可以表示为

其中：kx0（kxa）、ky0（kyb）分别代表 x=0（a）和 y=0（b）边上横向线性位移弹簧刚度系数，单位 N/m2；Kx0（Kxa）、Ky0（Kyb）分别代表 x=0（a）和 y=0（b）边上旋转约束弹簧刚度系数，单位 N/rad；Kxy0（Kxya）、Kyx0（Kyxb）分别代表 x=0（a）和 y=0（b）扭转约束弹簧刚度系数，单位 Nm/rad。

1.3 位移容许函数

在本文方法中，构造合适的位移容许函数对本文的研究有着很重要的关系。在以往的求解方法中，结构容许函数往往是根据相同边界下的梁函数进行设置，因此，不同的边界条件，需要自定义不同边界下的梁函数。这就会导致，当边界条件发生改变时需要设置不同类型的梁函数，当函数发生改变，整个求解过程需要重新推导与编程计算，不具有通用性，工作量繁琐。结构的容许函数除了设置成梁函数以外，还可以用简单、正交多项式和三角函数进行表示。当容许函数为多项式时，将容许函数展开为低阶多项式时，不能满足结构在高阶次的振动求解；展开为高阶多项式时，由于数值计算截断误差的原因会引起数值不稳定。当容许函数为传统傅立叶三角级数时，级数取无穷项时能够构成一个完整的无限维度向量空间，通过数值截断进行计算时，该容许函数具有良好的数值稳定性和一定的计算精度，也许可以克服多项式或者梁函数存在的问题。然而，传统的傅立叶级数在边界处会存在收敛性问题，除了简单的经典边界以外。因为，当容许函数展开为传统傅立叶三角级数时，其位移导数在边界处可能存在不连续现象，因此，可能导致不连续或者收敛速度极慢，不能满足结构的振动控制微分方程。为了克服这些难点，一种改进的傅立叶级数表达已经被提出，并且被应用到一般弹性边界下正交各向异性矩形薄板[18]的弯曲自由振动和环扇形板[19]的面内振动分析。在本文中，进一步将该方法扩展到变厚度中厚板结构的静动态特性分析中。因此，为了满足任意边界条件下的厚板结构，结构的位移函数与旋转函数被包含正弦三角级数的改进傅立叶级数。

在Mindlin板理论，位移分量假设由下式给出：

对于一般边界条件下中厚矩形板结构来说，其横向振动和转角位移函数可以采用二维改进傅里叶级数方法来进行表示：

1.4 能量泛函与求解步骤

在采用改进傅里叶级数建立了中厚矩形板的位移函数之后，需要对位移函数的未知傅里叶展开系数进行求解。由于本文所建立的位移函数足够光滑，位移函数的弱解（近似解）和强解（精确解）在数学意义上是等效的。因此，本文采用基于能量原理的Rayleigh-Ritz法来求解未知傅里叶展开系数。变厚度板的应变势能为：

将（1）、（2）式代入（15）式，结构的应变能可以表示为：

储存在边界的弹簧势能可以表示为：

当为一般弹性边界条件时，将弹簧考虑为无质量质点，此时结构的整体动能为

在理论模型的推导中，板的厚度h（ x，y）一般地表示为一个空间坐标函数。为了描述的统一和计算的方便，将空间坐标函数都采用一维或二维傅里叶级数来表示。由于板厚度的空间坐标函数为h（ x，y），板结构的厚度函数可以表示为二维傅里叶余弦级数

其中，傅里叶展开系数可以从下式进行计算：

为了简便起见，本文以沿x方向线性变化的中厚矩形板作为研究对象，假定厚度变化函数为

从而，板结构的厚度函数通过（19）式可以表示为一维的傅里叶余弦级数

其中，傅里叶展开系数根据（20）式可以计算得到

对施加在Mindlin矩形板上的外载荷所作功Wext可以表示为如下形式：

本文中将在两种不同类型的载荷下对Mindlin矩形板进行静态特性分析，这两种载荷形式分别为均布压强和线性压强。

其中：x1≤x≤x2，y1≤y≤y2为压强载荷作用有效区域。

为了得到结构位移函数的未知展开傅立叶系数，本文采用Rayleigh-Ritz方法进行求解。因为Rayleigh-Ritz方法作为一种数值近似解法，具有求解原理明晰、操作简单、计算效率高等优点。通过结构的应变势能、弹簧势能、动能以及外力功，结构的拉格朗日能量泛函可以表示为

将（10）-（12）式代入（26）式，对未知的傅立叶展开系数求极小值

结构的振动问题将转换为一个求特征值与特征向量的简单数学问题，将其矩阵化可以表示为

进行静态分析时，令（28）式中ω=0，则矩形板结构的横向弯曲振动位移的傅里叶系数向量可表示为：

得到位移的傅立叶系数后，代入（12）式即可得到矩形板的横向振动位移，通过对矩形板位移进行一系列数学操作即可得到其它感兴趣的物理量，诸如对位移求导能得到此时矩形板的速度和加速度等。

2 数值结果与分析

在本节中，将采用上一小节的理论模型，对不同边界下变厚度中厚板的静动态特性进行求解，将本文方法得到的计算结果与文献解及有限元计算结果进行对比验证，以验证本文方法的合理性。由于本文方法在计算过程中对级数进行阶段计算，所以首先对本文方法的收敛性进行分析；其次，对等厚度与变厚度下的静动态特性进行了研究分析。

2.1 本文方法收敛性研究

表1 不同截断数下CCCC匀厚中厚矩形板前8阶频率参数ΩTab.1 Frequency parameters Ω for CCCC moderately thick plates with uniform thickness and different truncate numbers

2.2 等厚度Mindlin矩形板静动态特性分析

本小节将对经典边界条件与弹性边界下的等厚度Mindlin矩形板进行静动态特性分析，静态特性主要分析Mindlin矩形板受局部压强下的位移变形，动态特性包括经典边界与一般弹性边界下的模态分析。在本小节所有算例计算中，矩形板的材料参数E=2×1011Pa，泊松比μ=0.3，密度为7800 kg/m3，矩形板结构参数设定为a=4 m，b=2 m；h/b=0.1，0.15，0.2。此外，由于缺少相关文献数据，为检验本方法的正确性，本小节的结果同有限元软件ABAQUS计算所得数据进行了比较，采用ABAQUS求解时矩形板所使用的单元类型为S4R单元，物理结构参数与边界条件与本文方法保持一致，网格尺寸为0.01 m×0.01 m。 其中表2-3 中（x1，x2，y1，y2）表示外界载荷作用在矩形板的范围，例如（2，3，1，2）表示外界载荷作用于矩形板的范围为x=2，x=3以及y=1，y=2所围成的矩形区域。

表2 均布压强P=1000 Pa作用于Mindlin矩形板下的最大变形 单位：mTab.2 The max structural deformation for moderately thick plates with uniform thickness and different boundary condition under load P=1000 Pa

表3 Mindlin矩形板在线性压强P（x）=100x+1 Pa作用下的最大变形 单位：mTab.3 The max structural deformation for moderately thick plates with uniform thickness and different boundary condition under load P(x)=100x+1 Pa

表2列出了均布压强在各种经典边界条件任意组合下在不同跨度比h/b、不同作用区域下结构的最大变形。表3列出了线性压强在各种经典边界条件任意组合下在不同跨度比h/b、不同作用区域下结构的最大变形。通过表2～3比较发现本文方法同有限元软件ABAQUS所得计算结果吻合良好，说明了本方法的正确性。同时也可以看出边界条件的变化对矩形板的变形影响很大，此外相同的边界条件和载荷形式及其大小，由于载荷作用的位置变化也会引起板结构最大变形的不同。图2与图3给出了 Mindlin 矩形板在 FCSC 边界条件下，均布压强载荷作用范围分别为（0，4，0，3）和（2，3，1，2）时Mindlin矩形板的最大变形图。由两组图形对比可知，两组图形中的矩形板受外界载荷后变形一致，进一步说明了本文方法的正确性。

接下来对等厚度中厚板的动态特性进行研究，由（28）式可知，模态特性可以通过求解一个标准的线性方程组简单得到。表4给出了不同经典边界条件及长宽比下结构前8阶固有频率参数，文献解及有限元分析结果作为参照值也在表4中列出，通过表4可知，本文方法的计算结果与文献解及有限元分析结果吻合良好。表5给出了弹性边界条件方板结构前6阶频率参数，边界处的弹性参数为kx0=kxa=ky0=kyb=Γw与Kx0=Kxa=Ky0=Kyb=Γx，扭转弹簧刚度系数与旋转弹簧保持一致，有限元分析结果作为参考值也列在表6中。通过表6可知，本文的计算结果与有限元计算结果吻合良好，因此本文方法对于任意边界条件下的等厚度Mindlin求解是可行的。

表4 不同边界条件下匀厚中厚矩形板前8阶频率参数ΩTab.4 Frequency parameters Ω for moderately thick plates with uniform thickness and different boundary conditions

表5 弹性边界条件下方板结构前6阶无量刚化固有频率参数ΩTab.5 Frequency parameters Ω for moderately square thick plates with uniform thickness and elastic rotation support

续表5

2.3 变厚度Mindlin矩形板静动态特性分析

表6 均布压强P=1000 Pa作用于Mindlin矩形板下的最大变形单位：mTab.6 The max structural deformation for moderately thick plates with variable thickness and different boundary condition under load P=1000 Pa unit：m

表7 Mindlin矩形板在线性压强P（x）=100x+1 Pa作用下的最大变形 单位：mTab.7 The max structural deformation for moderately thick plates with variable thickness and different boundary condition under load P(x)=100x+1 Pa unit：m

表8 不同边界条件下变厚度中厚矩形板前8阶频率参数ΩTab.8 Frequency parameters Ω for variable thickness square plates with different boundary condition and length-with ratio

表9 弹性边界条件下方板结构前6阶无量刚化固有频率参数ΩTab.9 Frequency parameters Ω for moderately square thick plates with variable thickness and elastic rotation support

3 结 论

本文基于傅立叶级数方法建立了任意边界条件下变厚度Mindlin板的静动态特性分析模型。该方法主要是将结构的位移容许函数利用改进傅立叶级数形式进行展开。在结构的边界处采用三类弹簧均匀布置来模拟任意边界支撑条件。将未知系数作为广义变量，结合Rayleigh-Ritz法对未知傅立叶展开系数求极值，将矩形厚板的弯曲振动问题转换为一个简单求解标准特征值问题。通过大量的数值算例，验证了本文方法的正确性和可靠性，并得出以下结论：

（1）任意边界条件下的变厚度Mindlin板位移容许函数可表示为一种通用的改进傅立叶三角级数形式；

（2）对级数进行截断后，随着截断项数的增加，计算结果快速收敛，并且数值稳定性很好；

（3）对于不同边界下的变厚度Mindlin板的静动态特性可以实现快速分析，无需重新推导与编程；

（4）除了单向线性变厚度，本文方法还可以很容易求解双向非线性厚度变化的结构的静动态特性，具有良好的扩展性。