1. 引言
弹性矩形薄板是工程中常见的结构形式，如国际空间站的太阳能电池帆板 [1] ，电子设备中的芯片，高速铁路中的轨道板 [2] 等，弹性矩形薄板的自由振动特性受到了各国学者的广泛关注，并且提出了相关的模型和计算方法。Ritz [3] 基于变分原理，假设薄板的形函数为三角函数与双曲函数的乘积，分析了自由边界条件下方形薄板的横向自由振动，得到了自由边界条件下方形薄板横向自由振动的特征值及节线图。Mindlin [4] 考虑了转动惯量和剪切变形，导出了各向同性弹性板的应力应变关系、运动微分方程以及能量函数。Gorman [5] 提出了求解四边固支矩形薄板自由振动解析解的叠加方法，该方法将振型函数表示为两个函数的乘积，并将振型函数分为双对称、双反对称、对称反对称三种形式，计算过程比较复杂。钟阳 [6] 从四边固支矩形薄板自由振动的基本方程出发，在哈密顿体系中使用辛几何方法得到了四边固支矩形薄板自由振动的精确解。Liew [7] 提出了pb-2瑞利里兹方法，该方法以边界函数和正交多项式构造出里兹试函数。Zhou和Cheung [8] 在pb-2瑞利里兹方法的基础上提出了切比雪夫里兹方法，分析了矩形板的三维振动，该方法以一组边界函数与切比雪夫多项式的乘积作为里兹试函数。陈林 [9] 等基于改进的傅里叶级数法研究了任意边界条件下矩形薄板的振动特性。薄板自由振动的控制方程为四阶椭圆型偏微分方程，上述方法一般都需要假设满足位移边界的位移试函数，根据变分原理得到薄板的特征频率方程。
近年来，一些新的理论和方法也得到了发展。Silling [10] 提出了一套非局部连续力学体系——近场动力学理论。同时，有关谱方法的研究也逐渐增多 [11] - [17] 。Orszag等根据不同的权函数对谱方法进行了分类，系统地总结了谱方法数值稳定性、收敛性的分析理论。Canuto等在前人研究基础上阐明了谱方法在流体力学应用中的关键内容，包括谱逼近的正交多项式、二维映射方法等。Shen等应用谱方法研究了高阶微分方程的边值问题。Shen等将谱方法与有限元方法结合，采用谱元方法对薄壁梁前4阶自由振动的频率和振型进行了分析。Trefethen提出了一种切比雪夫谱方法，该方法采用切比雪夫插值点构造了任意区间上任意连续函数的一阶求导矩阵。
从目前国内外的研究现状来看，采用切比雪夫谱方法研究板壳振动的文献不多。本文在Trefethen的研究基础上，利用一阶求导矩阵构造出满足四边固支矩形薄板自由振动微分方程的高阶求导矩阵，得到了四边固支矩形薄板自由振动的特征频率矩阵方程。与传统的Ritz法、Galerkin法等不同的是，由切比雪夫谱方法得到的特征频率矩阵方程不含有刚度矩阵和质量矩阵。通过算例和有限元方法比较，验证了切比雪夫谱方法的正确性。
2. 位移方程和边界条件
分析如图1所示弹性矩形薄板的横向自由振动问题，一般采用基尔霍夫假定，忽略剪切变形和转动惯量对薄板弯曲变形的影响，薄板横向自由振动的位移方程可以写成式(1) [18] ：
(1)
(1)式中：
为薄板的抗弯刚度，E为薄板材料的弹性模量，
为薄板材料的泊松比，
为薄板材料的密度，h为薄板的厚度。双调和算子的运算满足式(2)：
(2)
假设薄板在时域做简谐运动，将时间变量与空间变量分离，则横向自由振动位移w可以写成式(3)：
(3)
将(3)式带入(1)式，得到特征频率方程式(4)：
(4)
(4)式中：
为波数，
为薄板的自由振动频率。对于四边固支边界，薄板在边界上不发生横向位移，也不产生由位移变化引起的转角。因此，四边固支薄板的边界关系满足式(5)。根据文献 [14] 的研究，薄板的固支边界条件为微分方程的第一类边界条件和第二类边界条件的组合。
(5)
3. 切比雪夫谱方法
切比雪夫多项式是以俄罗斯数学家Pafnuty Lvovich Chebyshev命名的一系列多项式序列。切比雪夫多项式是一类完备的正交多项式，在数值计算中切比雪夫多项式的收敛性和数值稳定性优于一般的多项式，如泰勒级数 [19] 。第一类切比雪夫多项式的表达式如式(6)所示：
(6)
(6)式中：n为非负整数。第一类切比雪夫多项式存在交错点组，使得
，在数值分析中被称为切比雪夫插值点，是在区间[−1,1]选取的N + 1个不同的点，切比雪夫插值点的定义如式(7)所示：
(7)
Figure 2. Distribution of Chebyshev points (N = 8)
图2. N = 8时切比雪夫插值点的分布
当N = 8时，切比雪夫插值点的分布如图2所示。对于任意区间
，切比雪夫插值点满足式(8)的映射关系；同理，对于任意矩形区域
，
，切比雪夫插值点满足式(9)的映射关系：
(8)
(9)
根据文献 [18] 所提出的切比雪夫谱方法，区间[−1,1]上的1阶切比雪夫谱求导矩阵
满足式(10)：
(10)
(10)式中：i、j分别代表1阶切比雪夫谱求导矩阵中第i + 1行、第j + 1列，
满足式(11)：
(11)
采用切比雪夫谱方法求解高阶微分方程时，高阶切比雪夫谱求导矩阵可以表达为1阶切比雪夫谱求导矩阵的高次幂，如式(12)所示：
(12)
将对应的各阶切比雪夫谱求导矩阵带入特征频率方程(4)，可以得到如下的特征频率矩阵方程：
(13)
(13)式中：
代表克罗内克积 [20] 。
4. 数值算例与结果分析
4.1. 算例
算例为一四边固支方形薄板，如图3所示。边长为L = 2 m，弹性模量E = 2.10 × 109 Pa，薄板厚度h = 0.01 m，泊松比ν = 0.3，薄板的密度ρ = 7800 kg/m3。
4.2. 频率分析
有限元模型是在ANSYS 15.0中建立的，采用shell181单元，方形网格边长为0.01 m，网格数量40,000，节点数量40,401。在切比雪夫谱方法中，切比雪夫插值点数取25，求解算法在MATLAB R2017b中实现，矩阵规模25 × 25。四边固支方形薄板横向自由振动频率的计算结果如表1和图4所示。
分别采用切比雪夫谱方法和有限元方法计算了算例的前25阶频率。由表1可知，随着频率的增大，切比雪夫谱方法与有限元方法计算结果之间的误差越来越大，最大值出现在25阶频率，但是最大误差只有0.31%。同时在前25阶频率中，出现了6对重复的频率值，分别为2阶和3阶、7阶和8阶、9阶和10阶、14阶和15阶、18阶和19阶、23阶和24阶。与有限元方法所需的网格数量相比，切比雪夫谱方法选取的切比雪夫插值点个数很少，矩阵规模远远小于有限元方法，计算机内存占用更少，计算效率更高。
在上述研究的基础上，根据1~25阶频率的特点，选取了两组频率，分别得到了两组频率值随切比雪夫插值点数变化的规律。如图4所示，第一组为9阶和10阶，频率值相同；第二组为24和25阶，频率
Table 1. Frequencies and errors calculated by spectral method and FEM
表1. 谱方法和有限元方法计算的频率及误差
值不同。由图4可知，当插值点数大于12时，9阶和10阶频率值快速收敛，基本不再变化；当插值点数大于14时，24阶和25阶频率值快速收敛，基本不再变化。计算结果表明，采用切比雪夫谱方法研究算例的前25阶频率时，采用15~25个切比雪夫插值点可以保证计算结果的收敛性。
4.3. 振型分析
在频率分析的基础上，采用切比雪夫谱方法和有限元方法，研究了算例前25阶频率对应的振型。如图5(1)~图5(25)所示，各子图左侧为切比雪夫谱方法的计算结果，右侧为有限元方法的计算结果。算例1~25阶振型表现出以下规律：1) 切比雪夫谱方法和有限元方法计算的1~25阶振型结果基本规律一致，但不完全相同。两种方法计算的1~7、11~13、16、17、20~22、25阶振型基本相同，8~10、14、15、18、19、23、24振型存在一定差别。2) 差别振型来自重复频率值对应的振型。7阶和8阶、9阶和10阶、14阶和15阶、18阶和19阶、23阶和24阶振型都是重复频率值对应的振型，而且重复频率值对应的振型都是中心对称图形(包括2阶和3阶频率值对应的图形)。采用切比雪夫谱方法和有限元方法计算振型，最明显的差别存在于8阶振型。8阶和7阶振型对应的频率值相同，两种方法计算7阶振型的结果一致，但是两种方法计算得到的8阶振型却不同。这种现象也存在于18阶和19阶振型的计算结果中。实际上，7阶和8阶、18和19阶两对重复频率值对应的振型若采用有限元方法计算，结果是一致的，只是振型图绕中心点旋转了90度。但是7阶和8阶、18和19阶两对重复频率值对应的振型若采用切比雪夫谱方法计算，结果是不一致的。综上所述，在计算四边固支矩形薄板振型时，采用切比雪夫谱方法和有限元方法出现了一些不同的结果，有待进一步研究。
Figure 4. The relationship between frequencies and the number of interpolation points
图4. 频率与插值点数的关系
5. 结论
本文基于切比雪夫谱方法分析了四边固支矩形薄板的自由振动特性。基于薄板自由振动的位移方程，采用分离变量法推导了薄板自由振动的特征频率方程，建立了薄板各物理参数与自由振动频率之间的关系式。建立了任意区间、任意矩形区域与[−1, 1]区间内切比雪夫插值点的映射关系，构造出满足四边固支矩形薄板自由振动微分方程的高阶求导矩阵，得到了薄板自由振动的特征频率矩阵方程，该方程不含刚度矩阵和质量矩阵。算例薄板前25阶自由振动频率的计算结果表明，切比雪夫谱方法选取25个插值点的计算结果与有限元的计算结果相比，两者最大误差仅有0.31%，切比雪夫谱方法在计算9阶和10阶、24阶和25阶频率时，选取15~25个插值点，计算结果趋于稳定，收敛速度很快。因此本文提出的切比雪夫谱方法具有选取插值点少、收敛速度快、计算精度高的特点。在分析算例薄板的振型时，发现有限元方法与切比雪夫谱方法的计算结果存在某些差别，有待进一步研究。此外，基于切比雪夫谱方法的中厚板、各向异性材料板、不规则形状板以及其他边界条件下板的自由振动特性分析，有待深入研究。
基金项目
国家自然科学基金(11472057)，北京信息科技大学“勤信拔尖人才”培育计划项目(QXTCPB201701)，北京市自然科学基金(1182010)。