板壳结构广泛应用于各种工程领域，该种结构通常非常薄，其刚度、强度可能无法满足要求。由于增加板壳厚度会显著地增加结构的重量，此种方法并不可取。为了提高板的刚度和强度，加筋结构以其重量轻、刚度大等优势被广泛应用于工程结构中，于是其振动特性也非常值得关注。目前国内外对加筋板的振动问题的研究成果已经较为丰富，但是还并不完备，对于带有开口的加筋板结构及曲线型加筋结构的振动特性的分析相对还很少。

已有大量应用数值方法对加筋结构的振动性能的分析成果，如应用有限元法^[1]、微分求积法^[2]、光滑有限元法^[3]、无网格伽辽金法^[4]等。赵芝梅等^[5]应用ANSYS软件，对基于改进的子结构线导纳方法建立的多点激励下加筋板壳耦合结构的振动模型进行验证。Asokendu等^[6]采用有限元法对加筋壳体的自振特性进行分析，提出了一种新的壳单元，具有应用范围更广，更加灵活的特点。在分析较为复杂的结构问题时，数值方法具有较强的优势，但是解析方法作为一种能够探求结构内部机理的方法，对其进行研究具有重要的意义。

解析或半解析方法的应用甚至更早。王宏伟等^[7]应用了拉普拉斯变换技术对L型加筋板的自由振动进行研究，同时考虑了加强筋和板的材料阻尼影响。李凯等^[8]基于能量泛函变分的方法，研究了附加多个集中质量纵横加筋板的自由振动特性，引入拉格朗日乘子把板、梁组合振动分析问题转化为处理一类无约束泛函变分问题。杜菲等^[9]基于瑞利李兹法对四边固支的加筋板进行研究，通过与实验及ANSYS结果进行对比说明计算的准确性。应用解析方法对加筋结构进行分析具有良好的适用性^[10-11]，但是以上研究多是对简单矩形板的直线型加筋的分析，对带有复杂开口加筋矩形板，开口曲线加强筋的分析还较少。

本文应用瑞利李兹法对带有复杂开口的矩形加筋板结构的自振特性进行分析，引入弹簧模型模拟各种复杂的边界条件，并对结构的模态振型进行分析，为实际的工程问题中加筋结构的振动性能的分析提供参考。

1 理论分析

本文的研究对象为带有开口加强筋的加筋矩形板，示意图如图 1，板面为一矩形板，长为a，宽为b，中心为一开口，板面带有横向和纵向的加筋，开口处存在加强筋，四周为沿矩形板边界线性分布的位移约束线弹簧和转角约束线弹簧，用以模拟复杂的外边界条件。

边界处存在2种线性分布的弹簧，模拟边界条件时改变两类弹簧刚度的取值。表 1给出了几种边界条件的弹簧取值。

应用瑞利李兹法时，位移试函数对精度影响较大。改进的傅里叶级数可以克服边界处不连续现象，且不满足特定的边界条件，边界条件只与弹簧的刚度系数有关，适用于模拟各种复杂的边界条件。本文方法选用改进的傅里叶级数作为试函数^[12-13]，具体形式为：

$ \left\{ {\begin{array}{*{20}{l}} {{f_m}(x) = \sin \left( {m{\rm{ \mathsf{ π} }}{x^*}} \right),\quad 0 < m < 5}\\ {{f_m}(x) = \cos \left[ {(m - 5){\rm{ \mathsf{ π} }}{x^*}} \right],\quad m \ge 5}\\ {{g_n}(y) = \sin \left( {n{\rm{ \mathsf{ π} }}{y^*}} \right),\quad 0 < n < 5}\\ {{g_n}(y) = \cos \left[ {(n - 5){\rm{ \mathsf{ π} }}{y^*}} \right],\quad n \ge 5} \end{array}} \right. $

(1)

式中：m=1，2，3，…, M；n=1，2，3，…, N。

求解整体的应变能、动能、弹性势能时，将整个结构视为开口矩形板和加筋的组合，分别求解相关的能量，最后再进行求和。

在计算开口矩形板时，应用瑞利李兹法，首先假设一个位移试函数^[14]：

$ w(x,y,t) = \left( {\sum\limits_{m = 1}^M {\sum\limits_{n = 1}^N {{A_{mn}}} } {f_m}(x){g_n}(y)} \right){{\rm{e}}^{{\rm{i}}\omega t}} $

(2)

式中：A_mn为傅里叶级数展开系数；简谐时间因子e^iωt表示垂向位移与时间相关的项；M、N为截断项数；f_m(x)、g_n(y)为x、y方向的容许梁函数，其具体表达式见式(1)。

无开口的矩形板的弯曲应变能表示为:

$ \begin{array}{l} {V_p} = \frac{D}{2}\int_0^a {\int_0^b {\left[ {{{\left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}}} \right)}^2} + {{\left( {\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right)}^2} + 2\mu \left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}}} \right)\left( {\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right) + } \right.} } \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {2(1 - \mu ){{\left( {\frac{{{\partial ^2}w}}{{\partial x\partial y}}} \right)}^2}} \right]{\rm{d}}x{\rm{d}}y \end{array} $

(3)

式中：$D=\frac{E h^{3}}{12\left(1-\mu^{2}\right)}$，E、μ分别为板材料的杨氏模量和泊松比。

无开口的矩形板的动能为：

$ T = \frac{{\rho h}}{2}\int_0^a {\int_0^b {{{\left( {\frac{{\partial w}}{{\partial t}}} \right)}^2}} } {\rm{d}}x{\rm{d}}y $

(4)

式中ρ为材料的密度。

开口部分结构的弯曲应变能为：

$ \begin{array}{l} {V_{po}} = \frac{D}{2} {\iint_s {\left[ {{{\left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}}} \right)}^2} + {{\left( {\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right)}^2} + 2\mu \left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}}} \right)\left( {\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right) + } \right.} } \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {2(1 - \mu ){{\left( {\frac{{{\partial ^2}w}}{{\partial x\partial y}}} \right)}^2}} \right]{\rm{d}}x{\rm{d}}y \end{array} $

(5)

式中s为开口部分形状。

开口部分动能可以表示为：

$ {T_o} = \frac{{\rho h}}{2} {\iint_s {{{\left( {\frac{{\partial w}}{{\partial t}}} \right)}^2}} } {\rm{d}}x{\rm{d}}y $

(6)

储存在边界约束弹簧中的弹性势能为：

$ \begin{array}{l} {V_s} = \frac{1}{2}\int_0^a {{{\left[ {{k_{y0}}{w^2} + {K_{y0}}{{\left( {\frac{{\partial w}}{{\partial y}}} \right)}^2}} \right]}_{y = 0}}} {\rm{d}}x + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}\int_0^a {{{\left[ {{k_{yb}}{w^2} + {K_{yb}}{{\left( {\frac{{\partial w}}{{\partial y}}} \right)}^2}} \right]}_{y = b}}} {\rm{d}}x + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}\int_0^b {{{\left[ {{k_{x0}}{w^2} + {K_{x0}}{{\left( {\frac{{\partial w}}{{\partial x}}} \right)}^2}} \right]}_{x = 0}}} {\rm{d}}y + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}\int_0^b {{{\left[ {{k_{xa}}{w^2} + {K_{xa}}{{\left( {\frac{{\partial w}}{{\partial x}}} \right)}^2}} \right]}_{x = a}}} {\rm{d}}y \end{array} $

(7)

式中：k_x0、k_y0、k_xb、k_yb分别为x=0、y=0、x=a、y=b处的位移约束弹簧的刚度值。K_x0、K_y0、K_xb、K_yb分别为x=0、y=0、x=a、y=b处的转角约束弹簧的刚度值。

计算筋条的弯曲应变能及动能时，选用与矩形板相同的位移试函数，通过使筋条的横向位移与筋条所在位置板的横向位移相等来保证筋条与板的位移连续性，将加筋视为欧拉梁模型。则横向、纵向加筋筋条的弯曲应变能、动能分别为：

$ \left\{ {\begin{array}{*{20}{l}} {{V_{bh}} = \frac{{{{(EI)}_{bh}}}}{2}\int_l {{{\left( {\frac{{{\partial ^2}{w_{bh}}}}{{\partial {x^2}}}} \right)}^2}} {\rm{d}}x}\\ {{V_{bv}} = \frac{{{{(EI)}_{bv}}}}{2}\int_l {{{\left( {\frac{{{\partial ^2}{w_{bv}}}}{{\partial {y^2}}}} \right)}^2}} {\rm{dy}}}\\ {{T_{bh}} = \frac{{{m_{bh}}}}{2}\int_l {{{\left( {\frac{{{\partial ^2}{w_{bh}}}}{{\partial {t^2}}}} \right)}^2}} {\rm{d}}x}\\ {{T_{bv}} = \frac{{{m_{bv}}}}{2}\int_l {{{\left( {\frac{{{\partial ^2}{w_{bv}}}}{{\partial {t^2}}}} \right)}^2}} {\rm{dy}}} \end{array}} \right. $

(8)

式中：E为梁的弹性模量；I为梁的惯性矩；l为梁长；m_bh和m_bv分别为横向加筋和纵向加筋的质量密度。

由于加筋的存在，使结构整体的中性轴与光板中面位置相比下移，在计算梁的惯性矩时，应当对其进行修正，修正后的偏心距为^[15]:

$ {e^*} = \frac{{{S_1}}}{{{S_1} + {S_2}}}e $

(9)

式中：S₁、S₂分别为板、梁横截面面积；e为梁的中性轴到板中面的距离；e^*为梁的中性轴到组合截面中性轴的距离。

在对曲线加筋进行计算时，有时形状较为不规则，无法直接积分求得弯曲应变能及动能，引入数值方法对积分进行计算，在计算积分时将梁沿长度进行离散，得到小单元的中心点坐标，再利用中点坐标得到微段的应变能和动能，再将所有微段的应变能和动能求和，得到整体梁的应变能和动能。

离散后的梁的弯曲应变能为^[16]：

$ {V_b} = \frac{{{{(EI)}_b}}}{2}\sum\limits_{i = 1}^Q {{{\left( {{{\left. {\frac{{{\partial ^2}w(x,y)}}{{\partial {\mathit{\boldsymbol{n}}^2}}}} \right|}_{x = {x_i},y = {y_i}}}} \right)}^2}} \Delta {l_i} $

(10)

离散后的梁的动能为：

$ {T_b} = \frac{{{m_b}}}{2}\sum\limits_{i = 1}^Q {{{\left( {{{\left. {\frac{{{\partial ^2}w(x,y)}}{{\partial {t^2}}}} \right|}_{x = {x_i},y = {y_i}}}} \right)}^2}} \Delta {l_i} $

(11)

式中：n为边界的法线方向；l为边界的总长度。如图 2^[16]。

式(11)中偏导数为：

$ \frac{{\partial w(x,y)}}{{\partial \mathit{\boldsymbol{n}}}} = \frac{{\partial w(x,y)}}{{\partial x}}\cos \theta + \frac{{\partial w(x,y)}}{{\partial y}}\sin \theta $

(12)

离散筋条如图 3^[16]。

每个微段的长度为^[16]：

$ \Delta {l_j} = \sqrt {\Delta {y^2} + \Delta x_j^2} $

(13)

于是系统的能量泛函可以表示为：

$ \begin{array}{*{20}{c}} {\varPi = {V_p} - {V_{po}} + {V_{bh}} + {V_{bv}} + {V_s}}\\ {T + {T_o} - {T_{bh}} - {T_{bv}}} \end{array} $

(14)

将式(3)~(8)代入式(14)，并求极值：

$ \frac{{\partial \varPi }}{{\partial {A_{mn}}}} = 0 $

(15)

式中A_mn是用来描述薄板弯曲振动的未知系数。于是可以表示为：

$ \left( {\mathit{\boldsymbol{K}} - {\omega ^2}\mathit{\boldsymbol{M}}} \right)\mathit{\boldsymbol{A}} = {\bf{0}} $

(16)

式中：K为弹性势能与整体结构应变能之和；M为质量矩阵；A为未知系数向量；ω为圆频率。

2 收敛性分析 2.1 收敛性

截断项数M、N对能否取得较为准确的结果影响较大，选取矩形加筋板进行收敛性分析，边界条件为四边自由，参数如下：矩形板长a=10 m，宽b=10 m，厚度h=0.05 m。材料为钢材：杨氏模量E=2.1×10¹¹ Pa，泊松比μ=0.3，密度ρ=7 850 kg/m³。加筋的截面为矩形截面，宽B=0.1 m，高H=0.1 m，加筋位置为横向加筋y=5处，纵向加筋x=4、6处。

选取自由边界对截断项数进行收敛性分析，与有限元模型进行对比。表 2给出加筋板在不同截断项数M=N=q时对应的前6阶固有频率。有限元网格数为36 860结果收敛。由表 2可知，当M=N=12时，矩形加筋板的固有频率结果已经收敛，后续算例中取M=N=12。

2.2 实例分析

以上述纵横加筋板为例，计算固有频率，与有限元结果进行对比，计算两者相对误差。

本文将边界条件进行简写，其中C表示固支；S表示简支；F表示自由；E表示弹性边界。

表 3中弹性边界，y=0边弹簧刚度系数为k=10⁵ N/m²，K=0 N/rad。

选取带有组合形状开口的矩形加筋板，对其自振特性进行分析，计算固有频率，得到振型图，与有限元结果进行对比，说明分析的准确性，示意图如图 1。

模型的主要参数为长A=10 m，宽B=10 m，板厚H=0.05 m，中心开口半圆半径为r=1，连接直边为l=2 m，筋条截面为矩形，宽b₁=0.1 m，高h₁=0.1 m，加筋位置为横向加筋y=5处，纵向加筋x=4，x=6处，开口处加强筋宽b₂=0.1 m，高h₂=0.1 m。有限元建模如图 4所示，计算时对网格划分进行收敛性分析，当网格划分为4 425时结果已收敛。

根据上面的分析，应用Matlab编程，选取不同边界条件对上述模型进行计算，并与有限元模型结果对比，结果如表 4所示。给出四边固支前六阶模态振型图的对比图如图 5，说明方法的准确性。

表 4据表明，本文方法计算结果与有限元仿真分析计算结果误差很小，图 5中2种方法计算所得模态振型图吻合良好。说明本方法在计算带有开口的加筋板自由振动的准确性。

3 结论

1) 传统瑞利里兹法结合整体能量减去开口部分能量、弹簧模型模拟边界条件等方法能够非常灵活地对开口加筋等振动问题进行计算。

2) 结合离散的方法能够非常灵活地实现对曲线型加筋的振动问题的计算，可以模拟围壁加强等情况。