梁结构广泛应用于船舶与海洋工程、建筑工程和航空航天等领域中，对于此类结构的振动特性及其参数影响规律研究受到科研工作者的广泛关注^[1-11]。随着材料技术和加工工艺的进步，复合材料由于自身比刚度大、可设计性强、机械性能突出等优点，成为了工程应用中的首选材料。复合材料往往呈现出较为明显的正交各向异性，对各向异性梁结构动力学特性的研究具有十分重要的意义。经典的梁理论模型主要包含Euler-Bernoulli和Timoshenko 2种理论。然而，欧拉-伯努利梁理论忽略了结构横向拉伸和剪切变形以及转动惯量的影响，其计算结果与实际情况相比存在较大误差。铁木申科梁理论虽然考虑了横向剪切变形和转动惯量的影响，但其引入的剪切修正因子对梁结构的几何截面形状和材料参数均具有较高的依赖性，故而不适用于复合材料梁结构的研究。另外，经典梁理论更多地关注于梁结构沿着特定方向上的形变或振动问题。当梁结构的长径比较小或截面为不规则形状时，梁结构在外界载荷作用下不仅会出现弯曲变形，亦会产生扭转、拉伸或弯扭组合变形。此时需寻求具有较高计算精度和适用性的三维理论。Carrera等^[12]基于卡诺统一定理提出了一种改进的高阶梁理论模型。该理论将梁结构的位移容许函数分解为截面和轴向2部分，利用高阶泰勒公式拟合梁结构的截面位移，从而以泰勒展开的阶次控制理论的计算精度。该理论模型也被成功地应用于求解一些复杂问题，如薄壁结构、复合材料问题、旋转叶片等^[13-16]，并展现出较高的计算精度。

本文结合卡雷拉理论(Carrera unified formulation, CUF)和有限元理论建立正交各向异性梁结构的三维振动分析模型。基于梁结构的应变能和动能，由最小势能原理推导梁单元的质量矩阵和刚度矩阵。分别从质量和刚度矩阵中提取出3×3阶的质量、刚度核心矩阵，利用循环得到整个梁的质量矩阵和刚度矩阵。随后对多种边界条件下各向异性梁结构进行了自由振动分析，并将计算结果与有限元分析软件结果对比分析。

1 正交各向异性梁三维理论模型 1.1 梁结构位移场离散

如图 1所示，考虑一矩形正交各向异性梁结构，建立直角坐标系(x，y和z)。正切轴x与生长环相切，纵向轴y与纤维方向平行, 径向轴z与生长环垂直。梁形结构上任意一点沿坐标x，y和z 3个方向上的位移变形分别用U_x，U_y和U_z表示。梁结构的长度为L，宽度和高度分别为b和h。

为了简化三维实体模型，本文通过分离变量将梁形结构的位移函数分为截面位移函数和轴向位移函数2个部分。其中，利用卡诺高阶理论对梁结构截面面内位移进行拟合，卡诺高阶理论中用到的二维泰勒展开分量可以简写为统一的形式。泰勒展开的阶次可直接控制该方法的计算精度。梁结构的位移函数为：

$ \left\{ {\begin{array}{*{20}{l}} {{U_x} = {u_{x1}} + x{u_{x2}} + \cdots + {z^N}{u_{x(N + 1)(N + 2)/2}} = \sum\limits_{\tau = 1}^\varGamma {{F_\tau }} {u_{x\tau }}}\\ {{U_y} = {u_{y1}} + x{u_{y2}} + \cdots + {z^N}{u_{y(N + 1)(N + 2)/2}} = \sum\limits_{\tau = 1}^\varGamma {{F_\tau }} {u_{y\tau }}}\\ {{U_z} = {u_{z1}} + x{u_{z2}} + \cdots + {z^N}{u_{z(N + 1)(N + 2)/2}} = \sum\limits_{\tau = 1}^\varGamma {{F_\tau }} {u_{z\tau }}} \end{array}} \right. $

(1)

式中：N为二维泰勒展开的阶次；Γ为二维泰勒展开的总项数；F_τ 表示泰勒展开中的第τ项；变量u_xτ、u_yτ和u_zτ为梁形结构的轴向位移变量，是关于轴向y的函数。

轴向位移变量u_xτ、u_yτ和u_zτ利用有限元方法进行轴向解域离散，可表示为：

$ {u_{x\tau }} = \sum\limits_{i = 1}^M {{A_{\tau i}}} {N_i},{u_{y\tau }} = \sum\limits_{i = 1}^M {{B_{\tau i}}} {N_i},{u_{z\tau }} = \sum\limits_{i = 1}^M {{C_{\tau i}}} {N_i} $

(2)

以2节点线性梁单元(M=2)为例，函数为：

$ {N_1} = \frac{1}{2}(1 - r),{N_2} = \frac{1}{2}(1 + r),r \in [ - 1,1] $

(3)

将式(2)代入式(1)中，梁结构的位移函数为：

$ \left\{ {\begin{array}{*{20}{l}} {{U_x}(x,y,z) = \sum\limits_{\tau = 1}^\varGamma {\sum\limits_{i = 1}^M {{A_{\tau i}}} } {F_\tau }{N_i}}\\ {{U_y}(x,y,z) = \sum\limits_{\tau = 1}^\varGamma {\sum\limits_{i = 1}^M {{B_{\tau i}}} } {F_\tau }{N_i}}\\ {{U_z}(x,y,z) = \sum\limits_{\tau = 1}^\varGamma {\sum\limits_{i = 1}^M {{C_{\tau i}}} } {F_\tau }{N_i}} \end{array}} \right. $

(4)

通过式(4)可以看出，对于任意一梁单元，每一节点有3×Γ个未知变量，用向量表示为：

$ \mathit{\boldsymbol{q}} = {({A_{11}},{B_{11}},{C_{11}}, \cdots ,{A_{\varGamma 1}},{B_{\varGamma 1}},{C_{\varGamma 1}})^{\rm{T}}} $

(5)

值得注意的是，通过增加分析中所使用的线性单元的数目，或使用高阶的插值函数，都能够提高有限元结果的精度。

1.2 梁结构自由振动控制方程

考虑结构发生微小变形，根据应变位移关系可得：

$ \left\{ {\begin{array}{*{20}{l}} {{\varepsilon _x} = \frac{{\partial {U_x}}}{{\partial x}},{\varepsilon _y} = \frac{{\partial {U_y}}}{{\partial y}},{\varepsilon _z} = \frac{{\partial {U_z}}}{{\partial z}}}\\ {{\gamma _{xy}} = \frac{{\partial {U_x}}}{{\partial y}} + \frac{{\partial {U_y}}}{{\partial x}},{\gamma _{xz}} = \frac{{\partial {U_x}}}{{\partial z}} + \frac{{\partial {U_z}}}{{\partial x}}}\\ {{\gamma _{yz}} = \frac{{\partial {U_y}}}{{\partial z}} + \frac{{\partial {U_z}}}{{\partial y}}} \end{array}} \right. $

(6)

根据广义胡克定律，正交各向异性结构的应力应变关系可表示为：

$ \left[ {\begin{array}{*{20}{l}} {{\sigma _x}}\\ {{\sigma _y}}\\ {{\sigma _z}}\\ {{\tau _{xy}}}\\ {{\tau _{yz}}}\\ {{\tau _{xz}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{C_{11}}}&{{C_{12}}}&{{C_{13}}}&0&0&0\\ {{C_{12}}}&{{C_{22}}}&{{C_{23}}}&0&0&0\\ {{C_{13}}}&{{C_{23}}}&{{C_{33}}}&0&0&0\\ 0&0&0&{{C_{44}}}&0&0\\ 0&0&0&0&{{C_{55}}}&0\\ 0&0&0&0&0&{{C_{66}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\varepsilon _x}}\\ {{\varepsilon _y}}\\ {{\varepsilon _z}}\\ {{\gamma _{xy}}}\\ {{\gamma _{yz}}}\\ {{\gamma _{xz}}} \end{array}} \right] $

(7)

正交各向异性材料的弹性系数可表示为^[17]：

$ {C_{11}} = \frac{{1 - {v_{23}}{v_{32}}}}{{{E_2}{E_3}\varDelta }},{C_{22}} = \frac{{1 - {v_{13}}{v_{31}}}}{{{E_1}{E_3}\varDelta }},{C_{33}} = \frac{{1 - {v_{12}}{v_{21}}}}{{{E_1}{E_2}\varDelta }} $

(8)

$ {C_{12}} = \frac{{{v_{12}} + {v_{32}}{v_{13}}}}{{{E_1}{E_3}\varDelta }},{C_{13}} = \frac{{{v_{13}} + {v_{12}}{v_{23}}}}{{{E_1}{E_2}\varDelta }},{C_{12}} = \frac{{{v_{32}} + {v_{12}}{v_{31}}}}{{{E_1}{E_3}\varDelta }} $

(9)

$ \varDelta = \frac{{1 - {v_{12}}{v_{21}} - {v_{23}}{v_{32}} - {v_{31}}{v_{13}} - 2{v_{21}}{v_{32}}{v_{13}}}}{{{E_1}{E_2}{E_3}}} $

(10)

另外，υ_ij和υ_ji存在以下关系：

$ {v_{ij}}{E_j} = {v_{ji}}{E_i} $

(11)

如此，需要由9个独立参数来确定正交各向异性材料的物性参数：E₁、E₂、E₃、G₁₂、G₁₃、G₂₃、υ₁₂、υ₁₃、υ₂₃。

基于能量法推导梁单元的刚度矩阵梁单元的应变能：

$ {\varLambda ^{(e)}} = \int_V {\frac{1}{2}} {\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}\mathit{\boldsymbol{\sigma }}{\rm{d}}V = \frac{1}{2}{\mathit{\boldsymbol{q}}^{\rm{T}}}\mathit{\boldsymbol{Kq}} $

(12)

式中：ε为应变向量；σ为应力向量；K即为梁单元刚度矩阵。仔细观察方阵K的内部构造可以发现，该方阵是由一个3×3阶方阵K^ijτs堆积而成，其中的9个组成部分如下：

$ \left\{ {\begin{array}{*{20}{l}} {K_{xx}^{ij\tau s} = {C_{22}}{{\langle {F_{\tau ,x}}{F_{s,x}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l} + {C_{66}}{{\langle {F_{\tau ,z}}{F_{s,z}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l} + {C_{44}}{{\langle {F_\tau }{F_s}\rangle }_\varOmega }{{\langle {N_{i,y}}{N_{j,y}}\rangle }_l}}\\ {K_{xy}^{ij\tau s} = {C_{23}}{{\langle {F_{\tau ,x}}{F_s}\rangle }_\varOmega }{N_i}{N_{j,y}} + {C_{44}}{{\langle {F_\tau }{F_{s,x}}\rangle }_\varOmega }{{\langle {N_{i,y}}{N_j}\rangle }_l}}\\ {K_{xz}^{ij\tau s} = {C_{12}}{{\langle {F_{\tau ,x}}{F_{s,z}}\rangle }_\varOmega }{N_i}{N_j} + {C_{66}}{{\langle {F_{\tau ,z}}{F_{s,x}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l}} \end{array}} \right. $

(13)

$ \left\{ {\begin{array}{*{20}{l}} {K_{yx}^{ij\tau s} = {C_{44}}{{\langle {F_{\tau ,x}}{F_{s,z}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l} + {C_{23}}{{\langle {F_\tau }{F_{s,x}}\rangle }_\varOmega }{{\langle {N_{i,y}}{N_j}\rangle }_l}}\\ {K_{yy}^{ij\tau s} = {C_{55}}{{\langle {F_{\tau ,z}}{F_{s,z}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l} + {C_{44}}{{\langle {F_{\tau ,x}}{F_{s,x}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l} + {C_{33}}{{\langle {F_\tau }{F_s}\rangle }_\varOmega }{{\langle {N_{i,y}}{N_{j,y}}\rangle }_l}}\\ {K_{yz}^{ij\tau s} = {C_{55}}{{\langle {F_{\tau ,z}}{F_s}\rangle }_\varOmega }{{\langle {N_i}{N_{j,y}}\rangle }_l} + {C_{13}}{{\langle {F_\tau }{F_{s,z}}\rangle }_\varOmega }{{\langle {N_{i,y}}{N_j}\rangle }_l}} \end{array}} \right. $

(14)

$ \left\{ {\begin{array}{*{20}{l}} {K_{zx}^{ijs} = {C_{12}}{{\langle {F_{\tau ,z}}{F_{s,x}}\rangle }_\varOmega }\langle {N_i}{N_j}\rangle + {C_{66}}{{\langle {F_{\tau ,x}}{F_{s,z}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l}}\\ {K_{zy}^{ij\tau s} = {C_{13}}{{\langle {F_{\tau ,z}}{F_s}\rangle }_\varOmega }{{\langle {N_i}{N_{j,y}}\rangle }_l} + {C_{55}}{{\langle {F_\tau }{F_{s,z}}\rangle }_\varOmega }{{\langle {N_{i,y}}{N_j}\rangle }_l}}\\ {K_{zz}^{ij\tau s} = {C_{11}}{{\langle {F_{\tau ,z}}{F_{s,z}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l} + {C_{66}}{{\langle {F_{\tau ,x}}{F_{s,x}}\rangle }_\varOmega }{{\langle {N_i}{N_j}\rangle }_l} + {C_{55}}{{\langle {F_\tau }{F_s}\rangle }_\varOmega }{{\langle {N_{i,y}}{N_{j,y}}\rangle }_l}} \end{array}} \right. $

(15)

式中：${{\left\langle {{F}_{\tau }}{{F}_{s}} \right\rangle }_{\mathit{\Omega} }}\text{ } $表示截面内积分；$ {{\left\langle {{N}_{i}}{{N}_{j}} \right\rangle }_{l}}$表示轴向单元内积分。当梁结构的厚度较小时，此处的有限元方法会产生剪切自锁现象，本文用低阶的高斯-勒让德数值积分消除此处的剪切自锁现象^[18]。

梁单元的动能由各连续质点的动能累加而成：

$ T = \int_V {\frac{\rho }{2}} {\mathit{\boldsymbol{u}}^{\prime {\rm{T}}}}{\mathit{\boldsymbol{u}}^\prime }{\rm{d}}x $

(16)

式中：ρ为梁的密度；u′为梁上任意一点的速度，是关于时间t的函数。将位移函数代入上式中，则：

$ \begin{array}{*{20}{l}} {T = \frac{{\rho {\omega ^2}}}{2}\int_V {{N_i}} {F_\tau }\mathit{\boldsymbol{q}}_{\tau i}^{\rm{T}}{N_j}{F_s}{\mathit{\boldsymbol{q}}_{sj}}{\rm{d}}V = }\\ {{\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} \frac{{{\omega ^2}}}{2}\mathit{\boldsymbol{q}}_{\tau i}^{\rm{T}}{\mathit{\boldsymbol{M}}^{ij\tau s}}{\mathit{\boldsymbol{q}}_{sj}}} \end{array} $

(17)

与刚度矩阵相似，M^ijτs为质量矩阵的核心，质量矩阵可由3×3阶方阵M^ijτs堆积而成，其中的9个组成部分如下：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{M}}_{xx}^{ij\tau s} = \mathit{\boldsymbol{M}}_{yy}^{ij\tau s} = \mathit{\boldsymbol{M}}_{zz}^{ij\tau s} = \rho \int_\varOmega {{F_\tau }} {F_s}{\rm{d}}\varOmega \int_l {{N_i}} {N_j}{\rm{d}}y\\ \mathit{\boldsymbol{M}}_{xy}^{ij\tau s} = \mathit{\boldsymbol{M}}_{xz}^{ij\tau s} = \mathit{\boldsymbol{M}}_{yx}^{ij\tau s} = \mathit{\boldsymbol{M}}_{yz}^{ij\tau s} = \mathit{\boldsymbol{M}}_{zx}^{ij\tau s} = \mathit{\boldsymbol{M}}_{zy}^{ij\tau s} = 0 \end{array} \right. $

(18)

采用划行化列、赋大数等有限元边界施加方法处理正交各向异性梁结构的边界条件，正交各向异性梁自由振动控制方程可表示为：

$ ({\mathit{\boldsymbol{K}}_z} - {\omega ^2}{\mathit{\boldsymbol{M}}_z})\mathit{\boldsymbol{q}} = 0 $

(19)

式中：K_z为梁的总刚度矩阵；M_z为总质量矩阵。

2 正交各向异性梁自由振动分析

本文通过对一正交各向异性梁结构的自由振动分析来验证本文模型的合理性。考虑梁结构的材料为正交各向异性材料，其材料属性可表示为：E₁=145 GPa，E₂=9.6 GPa，G₁₂=G₁₃=4.1 GPa，G₂₃=3.4 GPa，υ₁₂=0.3，密度ρ=1 389 kg/m³。为了方便与文献[19]对比，纤维方向角取为90°，梁结构的长高比取为10，考虑两端固支边界条件，并引入无量纲频率参数：

$ {\omega ^*} = \omega {L^2}\sqrt {I/(QJ)} $

(20)

式中：$ Q=\frac{E_{1}^{2}}{E_{1}-v_{12}^{2} E_{2}} ; I=\rho b h ; J=b h^{3} / 12$。

表 1给出了正交各向异性梁结构的前6阶弯曲模态的无量纲频率参数对比结果。本文计算结果与文献[19]结果具有较好的一致性。值得注意的是，本文所建立的梁结构理论模型为三维模型，不仅可以预测梁结构的横向弯曲振动，也可以获得梁结构的扭转模态和纵振模态。

为了进一步验证本文三维振动分析结果的收敛性和正确性，本文将与有限元商业软件ANSYS的三维仿真结果进行对比。选取梁结构材料(graphite fabric-carbon matrix layers，GFCML)，其材料属性可表示为：E₁=173.1 GPa，E₂=33.1 GPa，E₃=5.17 GPa，G₁₂=9.38 GPa，G₁₃=8.27 GPa，G₂₃=3.24 GPa，υ₁₂=0.036，υ₁₃=0.25，υ₂₃=0.171，密度ρ=1 650 kg/m³。令梁结构的宽度和高度相等b=h=0.2 m，长度L为2 m。为了较为准确地预测梁结构的扭转模态，本文选取截面拟合多项式的阶次N为4。

表 2为正交各向异性梁结构在简支边界条件下，前10阶模态的收敛情况，每一行代表结构不同的模态阶次，每一列代表着不同的单元类型和单元数。B₂表示2节点线性单元，B₃表示3节点抛物线单元。从表中可以较为清楚地看到随着单元数量的增加，本文方法的计算结果快速收敛至稳定值。选用较高的单元类型会进一步加快计算方法的收敛速度。同时，表 2也列出了利用ANSYS软件计算出的数据，对于该仿真结果，本文选用了SOLID186单元类型，划分640个单元，3 665个节点。其有限元网格划分结果如图 2所示。另外，2种计算方法所需的自由度数(DOFs)均在表 2中给出。通过对比可以发现，本文方法所需较少的自由度数而获得与ANSYS相近的计算结果。

表 3为3种经典边界条件(简支、悬臂和固支边界)下正交各向异性梁结构的前10阶固有频率。通过与ANSYS仿真结果的对比，可以发现，对于不同的约束条件，本文方法均具有较高的计算精度。另外图 3给出了简支边界条件下正交各向异性梁结构的振型图，图 3(a)~(c)为ANSYS仿真结果，图 3(d)~(f)为本文方法绘制的结果。(a)和(d)为梁结构第1阶弯曲振型，(b)和(e)为结构1阶扭转振型，(c)和(f)为结构第1阶纵振振型。可以看出本文方法不仅可以用于预测正交各向异性梁结构的弯曲模态，也可以较为准确地预测梁结构的扭转模态和纵振模态。

3 结论

1) 本文方法不仅可以提供高精度的解，还可以较为全面地展现其振动特性(弯曲模态、扭转模态、纵振模态)。

2) 本文方法所需较少的自由度数而获得与商业有限元相同的计算精度。

3) 对于正交各向异性梁结构的几何参数设计及动力学特性优化具有指导意义。