相比于各向同性结构，各向异性结构^[1-6]由于拥有高刚度比和强韧性等特点而被广泛应用于各种工程领域中。三维各向异性曲梁和直梁的振动特性分析一直是振动噪声控制领域的热门课题^[7-21]。学者利用各种数值计算方法对三维直梁和曲梁进行结构建模和动力学特性分析，如传统有限元法，动刚度法、傅里叶法等。然而，大部分的数值方法在计算曲梁时很难保证结构几何的精确性和高阶函数连续等问题。等几何方法^[22]是一种能够实现计算机辅助设计(computer aided design，CAD)与计算机辅助工程(computer aided engineering，CAE)的无缝连接，并具有很高的精确性的数值方法。由于该方法具有高精度，高收敛，网格细化方便与高阶函数连续性等优点而被应用于求解各复合的梁、板和壳的力学问题^[23-31]。本文针各向异性的直梁和曲梁结构，并基于三维弹性理论，采用非均匀有理B样条(non-uniform rational B-spline，NURBS)函数，同时结合有限元思想进行几何建模和振动分析。通过一系列的数值算例验证该方法计算三维各向异性直梁和曲梁振动特性的快速收敛性与良好的精确性，同时还分析了不同边界下的参数影响。

1 三维各向异性梁

如图 1所示为三维各向异性梁结构，其中a、b和h分别为梁的长、宽和高。笛卡尔坐标选择建立在三维梁的轴向起始面。

基于三维弹性线性方程，各向异性梁的正应力和切应力可表示为:

$ \left[ {\begin{array}{*{20}{l}} {{\sigma _x}}\\ {{\sigma _y}}\\ {{\sigma _z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{C_{11}}}&{{C_{12}}}&{{C_{13}}}\\ {{C_{12}}}&{{C_{22}}}&{{C_{23}}}\\ {{C_{13}}}&{{C_{23}}}&{{C_{33}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{u_{,x}}}\\ {{v_{,y}}}\\ {{w_{,z}}} \end{array}} \right] $

(1)

$ \left[ {\begin{array}{*{20}{l}} {{\tau _{yz}}}\\ {{\tau _{xz}}}\\ {{\tau _{xy}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{C_{44}}}&0&0\\ 0&{{C_{55}}}&0\\ 0&0&{{C_{66}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{v_{,z}} + {w_{,y}}}\\ {{u_{,z}} + {w_{,x}}}\\ {{u_{,y}} + {v_{,x}}} \end{array}} \right] $

(2)

式中：带逗号的下角标x、y、z表示相应的位移分量的一阶偏导数；u、v和w分别为x、y和z方向上的位移分量；C_ij(i, j=1, 2，…，6)为弹性材料系数，其具体表达式为：

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

(3)

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

(4)

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

(5)

$ {\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}}}} $

(6)

$ {{C_{44}} = {G_{23}},{C_{55}} = {G_{13}},{C_{66}} = {G_{12}}} $

(7)

式中：E₁、E₂、E₃为3个方向上的杨氏模量；υ_ij为泊松比，其满足以下方程：

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

(8)

2 三维等几何有限元 2.1 非均匀有理B样条基函数

设节点矢量E=(ξ₁, ξ₂, ξ₃, …, ξ_n+p, ξ_n+p+1)，ξ_i(i=1, 2，…，n+p+1)为非负不减参数。则B样条的表达式^[32]为:

1) p=0:

$ {B_{i,0}}(\xi ) = \left\{ {\begin{array}{*{20}{l}} {1,{\xi _i} \le \xi \le {\xi _{i + 1}}}\\ {0,{\rm{ 其他 }}} \end{array}} \right. $

(9)

2) p>0:

$ {B_{i,p}}(\xi ) = \frac{{\xi - {\xi _i}}}{{{\xi _{i + p}} - {\xi _i}}}{B_{i,p - 1}}(\xi ) + \frac{{{\xi _{i + p + 1}} - \xi }}{{{\xi _{i + p + 1}} - {\xi _{i + 1}}}}{B_{i + 1,p - 1}}(\xi ) $

(10)

式中：p为基函数的阶数；n为基函数总数。基于节点矢量E=(0, 0, 0, 0.5, 1, 1, 1)的2阶的B样条函数如图 2所示。

在实际的建模仿真中，B样条函数并不能很好的满足工程设计师的要求。通过对每个B样条函数引入相应的权值，可以构建非均匀有理B样条函数(NURBS), 其一维函数表达式为：

$ {N_{i,p}}(\xi ) = \frac{{{w_{(i)}}{B_{i,p}}(\xi )}}{{\sum\limits_{i = 1}^n {{w_{(i)}}} {B_{i,p}}(\xi )}} $

(11)

式中：B_{i, p}(ξ)为B样条函数；w_i为相对应的权值。三维的NURBS为：

$ N_{i,j,k}^{p,q,r}(\xi ,\eta ,\zeta ) = \frac{{{w_{(i,j,k)}}{B_{i,p}}(\xi ){B_{j,q}}(\eta ){B_{r,k}}(\zeta )}}{{\sum\limits_{k = 1}^l {\sum\limits_{j = 1}^m {\sum\limits_{i = 1}^n {{w_{(i,j,k)}}} } } {B_{i,p}}(\xi ){B_{j,q}}(\eta ){B_{r,k}}(\zeta )}} $

(12)

式中：B_{i, p}(ξ)、B_{j, q}(η)、B_{r, k}(ζ)为3个参数空间上的B样条函数；w_{(i, j, k)}为对应的权值。

2.2 三维等几何振动控制方程

在三维等几何分析中，采用有限元思想，将作为形函数对未知位移分量的NURBS样条函数：

$ {\mathit{\boldsymbol{u}}^h}(\xi ,\eta ,\zeta ) = \sum\limits_{a = 1}^{n \times m \times l} {{N_a}} (\xi ,\eta ,\zeta ){\mathit{\boldsymbol{u}}_a} $

(13)

式中$a=i+(j-1) m+(k-1) m l ; \boldsymbol{u}_{a}=\left(u_{a}, v_{a}, w_{a}\right)^{\mathrm{T}} $为节点a的位移矢量；N_a为控制点a的NURBS函数。

三维各向异性梁结构的势能U和动能T为：

$ \begin{align} & U=\frac{1}{2}\iint\limits_{\varOmega }{({{\sigma }_{x}}{{\varepsilon }_{x}}+{{\sigma }_{y}}{{\varepsilon }_{y}}+{{\sigma }_{z}}{{\varepsilon }_{z}}+{{\tau }_{yz}}{{\gamma }_{yz}}+{{\tau }_{xz}}{{\gamma }_{xz}}+} \\ & \ \ \ \ \ \ {{\tau }_{xy}}{{\gamma }_{xy}})\text{d}x\text{d}y\text{d}z \\ \end{align} $

(14)

$ T=\frac{1}{2}\iint\limits_{\varOmega }{\left\{ \rho {{\left( \frac{\partial u}{\partial t} \right)}^{2}}+{{\left( \frac{\partial v}{\partial t} \right)}^{2}}+{{\left( \frac{\partial w}{\partial t} \right)}^{2}} \right\}\text{d}x\text{d}y\text{d}z} $

(15)

基于哈密尔顿原理，三维各向异性梁结构的虚功可表示为：

$ \delta \int_{{{t}_{1}}}^{{{t}_{2}}}{(T-U)}\text{d}t=0 $

(16)

将式(1)~(6)代入式(16)，并结合有限元方法，三维各向异性梁结构的振动方程为：

$ [\mathit{\boldsymbol{K}} - {\omega ^2}\mathit{\boldsymbol{M}}]\mathit{\boldsymbol{u}} = 0 $

(17)

式中：K和M分别为整体刚度矩阵和质量矩阵；ω为固有角频率；u为梁的控制点位移矢量。

3 三维梁等几何数值算例

通过前部分的理论推导和求解方法，本文将对各向异性的直梁和曲梁进行验证和数值结果的分析。并对异性的直梁和曲梁的参数变化对结构的振特性影响进行研究。

3.1 各向异性直梁数值结果

考虑到计算结果的精确性和计算机内存与计算速度之间的复杂关系。几何结构的控制点数选取对三维各向异性梁的振动是必不可少的。在本次收敛性分析中，三维各向异性直梁的边界考虑为C-C，梁的几何参数为a=1 m, b=h=0.1 m, 材料参数为E₁=40E₂，E₂=E₃=2 GPa, G₂₃=0.5E₂, G₁₃=G₁₂=0.6E₂, υ₁₂=υ₁₃=υ₂₃=0.25, ρ=1 500 kg/m³。无量纲频率$ \omega^{*}=\omega a \sqrt{\rho / E_{1}}$。图 3为三维各向异性梁结构的前4阶无量纲频率收敛性。在宽度和厚度方向上的NURBS样条函数的阶数取q=3, r=3。同时在宽度和厚度上单元数取4×4。在长度方向上的阶数p分别取2、3、4;单元取3~10个。表 1为三维各向异性直梁的前4阶振动频率f=ω/(2π)的对比。从表 1可以看出该方法在计算三维各向异性梁结构时拥有很好的收敛性；同时，样条函数的阶数越高，其收敛速度越快。

在本次对比中，梁结构的几何参数和材料属性与表 1中相同。在ANSYS计算中节点的数目大于一定数值时，其计算结果收敛。通过仿真计算可知传统有限元计算该结构的前4阶固有频率在节点大于2 793才趋于收敛。通过图 3可知，该方法计算结果在p=2时，控制点大于600收敛；在p=3时，控制点大于500收敛；在p=4时，控制点大于400收敛。可见该方法可通过提高函数的阶数，从而使计算效率优于传统有限元法。在本次对比分析中，ANSYS的数值结果采用SOLID45划分了24 696个节点计算而得。而本方法的数值结果采用p=q=r=3的NURBS并结合10×4×4单元(控制点=637)计算而得。通过数据对比可以看出，该方法计算出的振动频率与ANSYS的计算结果能很好的吻合。从而验证了该方法具有很好的精确性。

图 4为边界为C-S的各向异性直梁的前4阶模态图，进一步验证该方法计算三维各向异性直梁的精确性。该模态图计算所用的几何参数和材料参数与表 1中的一样。从图 4和表 1可以看出直梁结构的第1阶频率与第2阶频率相等。这是因梁的宽和高是相等的，且该结构的宽和高2个方向的材料属性也是相等的。因而，在该参数下三维梁的第1阶模态图与第2阶一样，只是运动方向不一样。通过模态图的对比，可以验证该方法不仅能够精确计算三维各向异性梁结构，同时也能精确的描述梁的振动模态。

为了更好地分析几何参数对各向异性梁的频率影响，在本次分析宽度b=0.3 m，长度a=1 m，厚度h为变量。材料参数保持与收敛性分析一样。图 5给出了长厚比对三维各向异性梁结构的前4阶无量纲频率的影响。通过图 5可以看出，不管在什么边界下，三维各向异性直梁结构的固有频率随厚度的增加而增加。当直梁的厚度增加时，会导致其在长度方向的抗弯刚度急剧增加。从而三维直梁的固有频率增加。当厚度h从0.05 m增加到0.125 m时，三维直梁的第3阶固有频率并没有变化，其原因有可能是该振动模态主要沿宽度方向振动。

3.2 各向异性曲梁数值结果

通过前文分析，可以看出该方法计算直梁频率的可靠性和精确性，然而在实际工程中，曲梁的应用比直梁更加广泛。曲梁的材料参数保持不变：E₁=40E₂，E₂=E₃=2 GPa, G₂₃=0.5E₂, G₁₃=G₁₂=0.6E₂, υ₁₂=υ₁₃=υ₂₃=0.25, ρ=1 500 kg/m³。曲梁的结构模型如图 6所示。

图 7为三维各向异性曲梁无量纲振动频率的收敛性，该收收敛性分析结构的边界、NURBS的阶数和几何划分的单元与直梁分析一样。曲梁的几何参数为：a=1 m，b=h=0.1 m, R=1 m。从图 7可以看出该方法计算曲梁时也拥有快速收敛的特性。表 2给出了三维各向异性曲梁的频率对比。ANSYS的数值结果采用SOLID45划分了25 872个节点计算而得。本方法的数值结果采用p=q=r=3的NURBS并结合10×4×4单元(控制点=637)计算而得。通过表 2可以看出该方法计算曲梁与ANSYS计算的曲梁频率结果也能吻合。表 2中的曲梁频率的误差比表 1中直梁的误差大，是因为ANSYS对曲面划分网格与本方法差距所造成的。图 8为三维各向异性曲梁在C-S边界下的前4阶模态对比图。从图中可以看出曲梁的模态对比图与ANSYS计算的模态图能很好地吻合。

图 9为厚度变化对三维各向异性曲梁无量纲固有频率的影响。在该分析中，材料参数保持不变。几何参数设置为a=1 m; b=0.1 m; R=0.1;从图 9可以看出，三维各向异性曲梁的无量纲频率在整体上随厚度的增加而增加；无量纲曲线的斜率随厚度的增加而减小。部分无量纲频率在一定范围内随着厚度的增加有可能保持不变，或者缓慢减少，如图 9中的C-C边界下的第2阶模态；C-C、C-S、S-S边界下的第3阶模态等。图 10为曲梁半径变化对三维各向异性曲梁的无量纲固有频率的影响。在该分析中，除了厚度h和半径R外，三维各向异性曲梁的材料参数和几何参数设置与图 9一样。曲梁的厚度h=0.1 m。从图 10可以看出，半径R对三维各向异性曲梁的无量纲固有频率的影响与厚度相似。即除了三维各向异性曲梁的部分频率在一定范围内不一定随半径的增加增加；但在整体上，三维各向异性曲梁的无量纲固有频率随半径增加而增加。

4 结论

1) 该方法在计算三维各向异性梁时具有快速收敛性，且基函数的阶数越高，计算结果收敛速度越快。

2) 该方法具有良好的精确性；其计算曲梁的固有频率与传统有限元结果吻合度低于该方法计算直梁固有频率和传统有限元结果的吻合度。

3) 边界约束越强，固有频率越高；随着半径的增加，曲梁的前四阶固有频率增加；整体上固有频率随着厚度的增加而增加，但部分阶次的固有频率有可能不变。

4) 借助NURBS强大的几何建模能力，该方法还适用于求解更为复杂的三维曲梁结构。