随着柔性梁在航空、航天和风力机等领域的广泛应用，研究者越来越重视其高精度的动力特性分析。如直升机螺旋桨^[1]、风力发电机叶片^[2]等,这类梁结构具有明显的非线性变形特点，工作时存在不同的回转速度。精确控制柔性旋转梁，需准确分析柔性梁的时域和频域特性。工程中为提高分析效率，通常采用简化梁模型，主要有线性有限元法、模态坐标近似法和柔性多体系统法，但这些方法存在模型精度较低、动力刚化模拟不明显、转动模态仿真不准确等问题。直接采用有限元模型，虽然能精确描述梁结构，但由于自由度显著增加，难以将其用于后续的多场耦合分析。现在计算大转动、非线性大变形梁的动力学问题的主流方法有GEBF(Geometrically Exact Beam Formula)和ANCF(Absolute Node Coordinate Formula)^[3]。前者在机翼气弹耦合分析领域得到广泛应用，但存在梁截面姿态插值问题，且广义坐标包含角度和位移，给频域分析的线性化处理带来较大困难。绝对节点坐标法由于质量矩阵为常数矩阵，且直接基于连续介质理论，能方便描述梁的几何非线性特征，并支持高效的时域求解。此外，在构建频域线性化动力学方程时，仅需对刚度矩阵线性化，因此适用于柔性旋转梁模型的时域和频域分析建模。目前，国内外已有专家开展了广泛研究。

Ahmed^[4]介绍了绝对节点坐标公式的计算机实现过程，为后续的数值仿真奠定了基础。Zhang等 ^[5-6]使用不同插值函数的ANCF对平面梁进行了动力特性分析，研究了弯曲和拉伸振动之间的模态耦合现象，但缺少弯扭耦合的相关研究。郑彤等^[7]对三维大变形柔性梁系统的动力学建模和仿真进行了研究。建立三维大变形柔性梁系统的动力学模型，解决了ADAMS在计算大变形物体动力学时的局限性问题，但并未涉及动力特性分析。赵春璋等^[8-12]建立动力学模型,研究了平面柔性梁大范围转动的变形、仿真精度的影响以及不同梁的理论适用范围，但未对频域响应进行相关的研究。姚震等^[13-14]对ANCF的精确性和灵敏度展开了研究，认为不同材料和截面特征对系统末端的运动规律有一定影响，同样未涉及动力特性分析方面的工作。

综上，目前研究主要集中于平面梁的动力学仿真和弯曲梁小变形情况下固有频率的计算，对基于连续介质理论和非线性变形假设的三维柔性旋转梁动力刚化及其动力特性的研究相对较少。

为此，本文基于绝对节点坐标法，结合连续介质理论和几何非线性假设，建立三维柔性梁的动力学模型。通过非线性梁自由单摆仿真实验，验证本文算法的精确性；然后对三维梁的动力学方程进行线性化，构建柔性旋转梁的频域振动方程，计算不同转速下梁的动力刚化现象，并研究柔性旋转梁在不同转速下的模态振型变化。为后续风力机或飞机叶片等的气弹阻尼分析提供指导。

1 绝对节点坐标柔性梁建模 1.1 梁单元位置插值

图1为三维二节点梁单元。其中， $ OXYZ $ 为绝对坐标系， $ (x,y,z) $ 为梁单元中任意一节点在单元坐标系中的坐标。 ${\boldsymbol{r}}^{p}$ 为梁单元中线上任意一点 p在绝对坐标系中的位置矢量；r₁，r₂，r₃，X，Y，Z为该点三维坐标参数； ${\boldsymbol{e}}\left(t\right)$ 为绝对坐标系中的单元节点坐标； ${{\boldsymbol{S}}}(x,y,z)$ 为单元形函数。

$ {{{\boldsymbol{r}}}}^{p}=\left[\begin{array}{c}{r}_{1}\\ {r}_{2}\\ {r}_{3}\end{array}\right]=\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]={{\boldsymbol{S}}}(x,y,z){{\boldsymbol{e}}}\left(t\right) $

(1)

本文采用三次插值多项式来描述单元中任意一点的绝对位置。

$ \left\{\begin{array}{l}{r}_{1}={a}_{0}+{a}_{1}x+{a}_{2}y+{a}_{3}{\textit{z}}+{a}_{4}xy+{a}_{5}x{\textit{z}}+{a}_{6}{x}^{2}+{a}_{7}{x}^{3}\\ {r}_{2}={b}_{0}+{b}_{1}x+{b}_{2}y+{b}_{3}{\textit{z}}+{b}_{4}xy+{b}_{5}x{\textit{z}}+{b}_{6}{x}^{2}+{b}_{7}{x}^{3}\\ {r}_{3}={c}_{0}+{c}_{1}x+{c}_{2}y+{c}_{3}{\textit{z}}+{c}_{4}xy+{c}_{5}x{\textit{z}}+{c}_{6}{x}^{2}+{c}_{7}{x}^{3}\end{array}\right. $

式中： $ {a}_{i} $ 、 $ {b}_{i} $ 、 $ {c}_{i} $ 为待定系数， $ x、y、z $ 为梁单元未变形时该点在单元坐标系下的坐标。

$ {{\boldsymbol{e}}}\left(t\right)={\left[\begin{array}{cc}{\boldsymbol{e}}_{i}& {\boldsymbol{e}}_{j}\end{array}\right]}^{\mathrm{T}} = {\left[\begin{array}{cccccccc}{\boldsymbol{r}}^{i}& {\boldsymbol{r}}_{x}^{i}& {\boldsymbol{r}}_{y}^{i}& {\boldsymbol{r}}_{z}^{i}& {\boldsymbol{r}}^{j}& {\boldsymbol{r}}_{x}^{j}& {\boldsymbol{r}}_{y}^{j}& {\boldsymbol{r}}_{z}^{j}\end{array}\right]}^{\mathrm{T}} $

(2)

式中： ${\boldsymbol{r}}^{i}$ 、 ${\boldsymbol{r}}^{j}$ 分别表示单元首尾节点在 $ OXYZ $ 下的位置矢量； ${\boldsymbol{r}}_{x}$ 表示位置矢量对梁单元坐标 $ x $ 的偏导数矢量， ${\boldsymbol{r}}_y$ 、 ${\boldsymbol{r}}_{\textit{z}}$ 同理。

将单元坐标式(2)与三次插值多项式联立，求解待定系数，可得到单元形函数^[15]。

$ \boldsymbol{S}\left(x,y,z\right)= \left[\begin{array}{cccccccc}{\boldsymbol{S}}_{1}\boldsymbol{I}& {\boldsymbol{S}}_{2}\boldsymbol{I}& {\boldsymbol{S}}_{3}\boldsymbol{I}& {\boldsymbol{S}}_{4}\boldsymbol{I}& {\boldsymbol{S}}_{5}\boldsymbol{I}& {\boldsymbol{S}}_{6}\boldsymbol{I}& {\boldsymbol{S}}_{7}\boldsymbol{I}& {\boldsymbol{S}}_{8}\boldsymbol{I}\end{array}\right] $

(3)

式中： $\boldsymbol{I}$ 为 $ 3\times 3 $ 的单位矩阵。

$ \left\{\begin{array}{ll}{\boldsymbol{S}}_{1}=1-3{\xi }^{2}+2{\xi }^{3};& \boldsymbol{S}_{2}=l\left(\xi -2{\xi }^{2}+{\xi }^{3}\right);\\ \boldsymbol{S}_{3}=l\left(\eta -\xi \eta \right);& \boldsymbol{S}_{4}=l\left(\zeta -\xi \zeta \right);\\ \boldsymbol{S}_{5}=3{\xi }^{2}-2{\xi }^{3};& \boldsymbol{S}_{6}=l\left({\xi }^{3}-{\xi }^{2}\right);\\ \boldsymbol{S}_{7}=l\left(\xi \eta \right);& \boldsymbol{S}_{8}=l\left(\xi \zeta \right)\end{array}\right. $

式中： $ \xi =x/l $ ， $ \eta =y/l $ ， $\zeta ={\textit{z}}/l$ ， $ l $ 为未变形时梁的单元长度。

1.2 单元质量矩阵

对式(1)进行时间求导可得到任意一点的绝对速度 $\dot{\boldsymbol{r}}$ 。

$ \dot{\boldsymbol{r}}=\boldsymbol{S}\dot{\boldsymbol{e}} $

(4)

梁单元的动能可表示为

$ {E}_{\mathrm{k}}=\frac{1}{2}\rho \iiint {\dot{\boldsymbol{r}}}^{\mathrm{T}}\dot{\boldsymbol{r}}\mathrm{d}V=\\ \frac{1}{2}{\dot{\boldsymbol{e}}}^{\mathrm{T}}\left(\rho \iiint {\boldsymbol{S}}^{\mathrm{T}}\boldsymbol{S}\mathrm{d}V\right)\dot{\boldsymbol{e}}=\frac{1}{2}{\dot{\boldsymbol{e}}}^{\mathrm{T}}\boldsymbol{M}\dot{\boldsymbol{e}} $

(5)

式中： $\boldsymbol{M}$ 为单元质量矩阵，该矩阵与广义坐标无关，为常数矩阵。

$ \boldsymbol{M}=\rho {\int }_{0}^{l}\left(\iint {\boldsymbol{S}}^{\mathrm{T}}\boldsymbol{S}\mathrm{d}A\right)\mathrm{d}x $

(6)

式中： $ \rho $ 、 $ \mathrm{d}A $ 、 $ l $ 分别为梁单元的密度、截面面积微元、梁未变形时的单元长度。

1.3 单元刚度矩阵

利用连续介质理论和非线性格林应变张量假设，计算梁单元的位移应变关系。梁单元的变形梯度 $\boldsymbol{H}$ 为

$ \boldsymbol{H}=\left[\begin{array}{ccc}\dfrac{{{\partial}} {r}_{1}}{ \partial x}& \dfrac{\partial {r}_{1}}{ \partial y}& \dfrac{\partial {r}_{1}}{ \partial z}\\ \dfrac{\partial {r}_{2}}{ \partial x}& \dfrac{\partial {r}_{2}}{ \partial y}& \dfrac{\partial {r}_{2}}{ \partial z}\\ \dfrac{\partial {r}_{3}}{ \partial x}& \dfrac{\partial {r}_{3}}{ \partial y}& \dfrac{\partial {r}_{3}}{ \partial z}\end{array}\right] $

(7)

非线性格林应变张量为

$ {{\boldsymbol{G}}}=\frac{1}{2}({\boldsymbol{H}}^{\mathrm{T}}\boldsymbol{H}-\boldsymbol{I}) $

(8)

将式(1)和式(7)代入式(8)得到

$ \boldsymbol{G}=\dfrac{1}{2}\left(\left[\begin{array}{lll}{\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{x}\boldsymbol{e}& {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{xy}\boldsymbol{e}& {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{{\textit{z}}x}\boldsymbol{e}\\ {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{xy}\boldsymbol{e}& {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{y}\boldsymbol{e}& {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{y{\textit{z}}}\boldsymbol{e}\\ {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{{\textit{z}}x}\boldsymbol{e}& {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{y{\textit{z}}}\boldsymbol{e}& {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{{\textit{z}}}\boldsymbol{e}\end{array}\right]-\boldsymbol{I}\right) $

(9)

式中：

$ \left\{\begin{array}{c}{\boldsymbol{S}}_{x}={\left(\dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial x}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial x}+{\left(\dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial x}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial x}+{\left(\dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial x}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial x}\\ {\boldsymbol{S}}_{y}={\left(\dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial y}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial y}+{\left(\dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial y}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial y}+{\left(\dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial y}\right)}^{{{\rm{T}}}} \dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial y}\\ {\boldsymbol{S}}_{{\textit{z}}}={\left(\dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial {\textit{z}}}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial {\textit{z}}}+{\left(\dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial {\textit{z}}}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial {\textit{z}}}+{\left(\dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial {\textit{z}}}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial {\textit{z}}}\\ {\boldsymbol{S}}_{xy}={\left(\dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial x}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial y}+{\left(\dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial x}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial y}+{\left(\dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial x}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial y}\\ {\boldsymbol{S}}_{y{\textit{z}}}={\left(\dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial y}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial {\textit{z}}}+{\left(\dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial y}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial {\textit{z}}}+{\left(\dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial y}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial {\textit{z}}}\\ {\boldsymbol{S}}_{{\textit{z}}x}={\left(\dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial {\textit{z}}}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{1}}{ \partial x}+{\left(\dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial {\textit{z}}}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{2}}{ \partial x}+{\left(\dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial {\textit{z}}}\right)}^{\mathrm{T}} \dfrac{\partial {\boldsymbol{S}}_{3}}{ \partial x}\end{array}\right. $

将对称的非线性格林应变张量表示为列矢量形式。

$ \boldsymbol{\varepsilon }={\left[\begin{array}{cccccc}{\boldsymbol{G}}_{11}& {\boldsymbol{G}}_{22}& {\boldsymbol{G}}_{33}& {2\boldsymbol{G}}_{12}& {2\boldsymbol{G}}_{13}& {2\boldsymbol{G}}_{23}\end{array}\right]}^{\mathrm{T}} $

(10)

在线弹性本构关系下，应力 $\boldsymbol{\sigma }$ 与应变 $\boldsymbol{\varepsilon }$ 的关系为

$ \boldsymbol{\sigma }=\boldsymbol{D}\boldsymbol{\varepsilon } $

(11)

式中： $\boldsymbol{D}$ 为材料弹性系数张量。

根据虚功原理弹性力所作虚功为

$ {\text{δ}} {V}_{{\rm{w}}}={\int }_{V}{\left(\boldsymbol{D}\boldsymbol{\varepsilon }\right)}^{\mathrm{T}}{\text{δ}} \boldsymbol{\varepsilon }\mathrm{d}{V}={\boldsymbol{e}}^{\mathrm{T}}\boldsymbol{K}{\text{δ}} \boldsymbol{e} $

(12)

式中： $\boldsymbol{K}$ 为单元刚度矩阵，V为梁单元的体积，可将单元刚度矩阵整理得到式(13)。

$ \boldsymbol{K}={\boldsymbol{K}}_{\mathrm{l}}+{\boldsymbol{K}}_{\mathrm{n}} $

(13)

式中： ${\boldsymbol{K}}_{\mathrm{l}}$ 为刚度矩阵的线性部分，与绝对节点坐标无关。

$ {\boldsymbol{K}}_{\mathrm{l}}=-\frac{\left(3\lambda +2\mu \right)}{2}\iiint ({\boldsymbol{S}}_{x}+{\boldsymbol{S}}_{y}+{\boldsymbol{S}}_{z})\mathrm{d}V $

(14)

${\boldsymbol{K}}_{\mathrm{n}}$ 为非线性部分， ${\boldsymbol{K}}_{\mathrm{n}}={\boldsymbol{K}}_{\mathrm{n},1}+{\boldsymbol{K}}_{\mathrm{n},2}+{\boldsymbol{K}}_{\mathrm{n},3}$ ，其中 $ \lambda $ 、 $ \mu $ 为lame常数，具体计算为

$ \left\{\begin{array}{l}{\boldsymbol{K}}_{\mathrm{n},1}=\dfrac{1}{2}(\lambda +2\mu )\displaystyle\iiint (({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{x}\boldsymbol{e}){\boldsymbol{S}}_{x}+\\ \qquad \;\;({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{y}\boldsymbol{e}){\boldsymbol{S}}_{y}+({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{{\textit{z}}}\boldsymbol{e}){\boldsymbol{S}}_{{\textit{z}}})\mathrm{d}V\\ {\boldsymbol{K}}_{\mathrm{n},2}=\dfrac{1}{2}\lambda \displaystyle\iiint (({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{x}\boldsymbol{e}){(\boldsymbol{S}}_{y}+{\boldsymbol{S}}_{{\textit{z}}})+\\ \qquad \;\;({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{y}\boldsymbol{e})({\boldsymbol{S}}_{x}+{\boldsymbol{S}}_{{\textit{z}}})+({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{{\textit{z}}}\boldsymbol{e}){(\boldsymbol{S}}_{x}+{\boldsymbol{S}}_{y}\left)\right)\mathrm{d}V\\ {{\boldsymbol{K}}_{\mathrm{n},3}=\mu \displaystyle\iiint (({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{xy}\boldsymbol{e})}({\boldsymbol{S}_{xy}}+{\boldsymbol{S}}_{xy}^{\mathrm{T}})+\\ \qquad \;\;({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{y{\textit{z}}}\boldsymbol{e}){(\boldsymbol{S}}_{y{\textit{z}}}+{\boldsymbol{S}}_{y{\textit{z}}}^{\mathrm{T}})+({\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{S}}_{{\textit{z}}x}\boldsymbol{e}){(\boldsymbol{S}}_{{\textit{z}}x}+{\boldsymbol{S}}_{{\textit{z}}x}^{\mathrm{T}}))\mathrm{d}V\end{array}\right. $

$ \begin{aligned}[b] &{\lambda }=\frac{v}{1+v} \cdot \frac{E}{1-2v} \\ &\mu =\frac{E}{2(1+v)} \end{aligned} $

(15)

式(15)中： $ E $ 为材料的弹性模量， $ v $ 为材料泊松比。

1.4 单元广义力矩阵

根据虚功原理以及式(1)，得到外力虚功为

$ {\text{δ}} {W}_{\mathrm{e}\mathrm{l}\mathrm{e}}={\boldsymbol{F}}^{\mathrm{T}}{\text{δ}} \boldsymbol{r}={\boldsymbol{F}}^{\mathrm{T}}\boldsymbol{S}{\text{δ}} \boldsymbol{e}={\boldsymbol{Q}}_{\mathrm{e}\mathrm{l}\mathrm{e}}^{\mathrm{T}}{\text{δ}} \boldsymbol{e} $

(16)

式中： $\boldsymbol{r}$ 是节点的位置矢量, $\boldsymbol{F}$ 为梁单元上任意一点所受的外力， ${\boldsymbol{Q}}_{\mathrm{e}\mathrm{l}\mathrm{e}}^{\mathrm{T}}$ 是与单元节点坐标相关的广义力的矢量。

$ {\boldsymbol{Q}}_{\mathrm{e}\mathrm{l}\mathrm{e}}^{\mathrm{T}}={\boldsymbol{F}}^{\mathrm{T}}\boldsymbol{S} $

(17)

2 柔性梁系统动力学方程 2.1 系统的约束表达

单摆仿真时，需要在单摆被约束的根部施加一转动副约束，使单摆绕着被约束点旋转，此外该转动副亦可施加于梁与梁之间，形成两梁之间的相对转动；多个梁单元时，单元与单元之间不会发生相对移动、扭转和旋转，单元与单元需固结在一起。故本文特介绍2种约束，一种为绕z轴旋转的自由单摆的转动副约束，一种为梁单元之间的固定铰约束。

转动副约束定义如图2(a)所示，点 $ p $ 是两单元连接处的节点，故两单元在点 $ p $ 处的绝对位置坐标和z方向的偏导坐标始终相同，x和y方向的偏导坐标可以不同。其约束方程为

$ {\boldsymbol{r}}^{ip}={\boldsymbol{r}}^{jp} $

(18)

$ {\boldsymbol{r}}_{{\textit{z}}}^{ip}={\boldsymbol{r}}_{{\textit{z}}}^{jp} $

(19)

固定铰约束如图2(b)所示，点 $ p $ 在 $ i $ 单元和 $ j $ 单元中的绝对位置坐标和偏导坐标始终相等。其约束方程为

$ {\boldsymbol{r}}^{ip}={\boldsymbol{r}}^{jp} \text{；} {\boldsymbol{r}}_{\alpha }^{ip}={\boldsymbol{r}}_{\alpha }^{jp},\;\alpha =x,y,z $

(20)

式中： ${\boldsymbol{r}}^{ip}$ 为 $ i $ 单元中 $ p $ 点的绝对坐标矢量， ${\boldsymbol{r}}_{\alpha }^{ip}$ 为 $ i $ 单元中 $ p $ 点关于单元坐标的偏导数。

2.2 系统动力学方程

由单元质量矩阵、单元刚度矩阵、单元广义力矩阵和组装布尔矩阵 $\boldsymbol{B}$ ，可获得梁模型的质量矩阵 ${\boldsymbol{M}}_{\mathrm{t}\mathrm{m}}$ 、刚度矩阵 ${\boldsymbol{K}}_{\mathrm{t}\mathrm{m}}$ 、广义力矩阵 ${\boldsymbol{Q}}_{\mathrm{t}\mathrm{m}}$ 。

$ \left\{\begin{array}{l}{\boldsymbol{M}}_{\mathrm{t}\mathrm{m}}={\boldsymbol{B}}^{\mathrm{T}}\boldsymbol{M}\boldsymbol{B}\\ {\boldsymbol{K}}_{\mathrm{t}\mathrm{m}}={\boldsymbol{B}}^{\mathrm{T}}\boldsymbol{K}\boldsymbol{B}\\ {\boldsymbol{Q}}_{\mathrm{t}\mathrm{m}}={\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{Q}}_{\mathrm{e}\mathrm{l}\mathrm{e}}\end{array}\right. $

(21)

$\boldsymbol{B}$ 是 $ m\times n $ 维的组装布尔矩阵， $ m $ 为梁单元坐标数与单元数目的乘积， $ n $ 为柔性梁系统节点坐标数^[16]。

考虑梁铰约束，将约束式(18)~(20)与式(21)联立。约束组装矩阵由 ${\boldsymbol{B}}_{\mathrm{c}}$ 表示， $\boldsymbol{\ddot e} $ 为加速度，约束柔性梁的动力学方程可表示为

$ \begin{array}{c}{\boldsymbol{M}}_{\mathrm{c}}\left({\boldsymbol{B}}_{\mathrm{c}}\ddot{\boldsymbol{e}}\right)+{\boldsymbol{K}}_{\mathrm{c}}\left({\boldsymbol{B}}_{\mathrm{c}}\boldsymbol{e}\right)={\boldsymbol{Q}}_{\mathrm{c}}\end{array} $

(22)

式中： ${\boldsymbol{M}}_{\mathrm{c}}$ 为约束后的质量矩阵， ${\boldsymbol{K}}_{\mathrm{c}}$ 为约束后的刚度矩阵， ${\boldsymbol{Q}}_{\mathrm{c}}$ 为约束后的广义力矩阵。

$ \left\{\begin{array}{l}{\boldsymbol{M}}_{\mathrm{c}}={{{\boldsymbol{B}}_{\rm{c}}^{\rm{T}}}\boldsymbol{M}}_{\mathrm{t}\mathrm{m}}{\boldsymbol{B}}_{\mathrm{c}}\\ {\boldsymbol{K}}_{\mathrm{c}}={{{\boldsymbol{B}}_{\rm{c}}^{\rm{T}}}\mathit{K}}_{\mathrm{t}\mathrm{m}}{\boldsymbol{B}}_{\mathrm{c}}\\ {\boldsymbol{Q}}_{\mathrm{c}}={{\boldsymbol{B}}_{\mathrm{c}}\boldsymbol{Q}}_{\mathrm{t}\mathrm{m}}\end{array}\right. $

2.3 旋转梁动力学方程

用绝对节点坐标方法建模时，节点坐标相对绝对坐标进行定义。当梁绕定轴转动时，为方便描述刚性大范围定轴转动和梁本身的变形运动，引入浮点坐标系，建立旋转梁动力学模型。

在三维空间坐标中，三维柔性梁绕着Z轴旋转，将其投影于旋转平面如图3所示。

图3中 ${r}_{1} $ 、 ${r}_{2} $ 、 ${\eta }_{1}$ 、 ${\eta }_{2} $ 分别表示点在绝对坐标系下的坐标和浮动坐标系下的坐标， $ \theta $ 为浮动坐标系相对于绝对坐标系的转角。该点在两坐标系中的关系如式(23)所示。

$ \left[\begin{array}{c}{r}_{1}\\ {r}_{2}\\ {r}_{3}\end{array}\right]=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\;\theta & -\mathrm{s}\mathrm{i}\mathrm{n}\;\theta & 0\\ \mathrm{s}\mathrm{i}\mathrm{n}\;\theta &\;\;\; \mathrm{c}\mathrm{o}\mathrm{s}\;\theta & 0\\ 0& \;\;\;0& 1\end{array}\right]\left[\begin{array}{c}{\eta }_{1}\\ {\eta }_{2}\\ {\eta }_{3}\end{array}\right]+R\left[\begin{array}{c}\mathrm{c}\mathrm{o}\mathrm{s}\;\theta \\ \mathrm{s}\mathrm{i}\mathrm{n}\;\theta \\ 0\end{array}\right] $

(23)

根据单元节点坐标 $\boldsymbol{e}\left(t\right)$ 的表示，结合式(23)，单元浮动坐标 ${\tilde {\boldsymbol{e}}_{{\rm{ele}}}} $ 与单元节点坐标 $\boldsymbol{e}\left(t\right)$ 关系为

$ \boldsymbol{e}\left(t\right)={\boldsymbol{U}}_{\mathrm{C}\mathrm{T}}{\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{l}\mathrm{e}}+{\boldsymbol{e}}_{\mathrm{e}\mathrm{l}\mathrm{e},R} $

(24)

式中： ${\boldsymbol{e}}_{\mathrm{e}\mathrm{l}\mathrm{e},R}$ 为轮毂半径坐标，R 为轮毂半径。

$ \begin{aligned}[b] &\qquad\qquad\qquad {\boldsymbol{e}}_{\mathrm{e}\mathrm{l}\mathrm{e},R}=R{\boldsymbol{U}}_{\mathrm{C}\mathrm{T}}\boldsymbol{V}\\ &\boldsymbol{V}={\left[\begin{array}{cccccccc}1& 0& 0& \cdots & 1& 0& 0& \cdots \end{array}\right]}^{\mathrm{T}} \end{aligned} $

(25)

$ {\boldsymbol{U}}_{\mathrm{C}\mathrm{T}}=\left[\begin{array}{cc}{\boldsymbol{N}}_{\mathrm{C}\mathrm{T}}& 0\\ 0& {\boldsymbol{N}}_{\mathrm{C}\mathrm{T}}\end{array}\right] $

(26)

式(25)~(26)中： $\boldsymbol{V}$ 为单元坐标组装的布尔向量， ${\boldsymbol{N}}_{\mathrm{C}\mathrm{T}}$ 为单元节点坐标变化矩阵。

$ {\boldsymbol{N}}_{\mathrm{C}\mathrm{T}}=\left[\begin{array}{cccc}{\boldsymbol{P}}_{0}& 0& 0& 0\\ 0& {\boldsymbol{P}}_{0}& 0& 0\\ 0& 0& {\boldsymbol{P}}_{0}& 0\\ 0& 0& 0& {\boldsymbol{P}}_{0}\end{array}\right] $

(27)

$ {\boldsymbol{P}}_{0}=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\;\theta & -\mathrm{s}\mathrm{i}\mathrm{n}\;\theta & 0\\ \mathrm{s}\mathrm{i}\mathrm{n}\;\theta & \mathrm{c}\mathrm{o}\mathrm{s}\;\theta & 0\\ 0& 0& 1\end{array}\right] $

(28)

将式(24)代入式(22)中，旋转梁动力学方程为

$ {\boldsymbol{M}}_{\mathrm{c}}{\ddot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}+2\dot{\theta }{\boldsymbol{M}}_{1\mathrm{c}}{\dot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}-{\dot{\theta }}^{2}{\boldsymbol{M}}_{2\mathrm{c}}{\tilde {\boldsymbol{e}}}_{\mathrm{c}} - R{\dot{\theta }}^{2}{\boldsymbol{M}}_{2\mathrm{c}}\boldsymbol{V}+{\boldsymbol{K}}_{\mathrm{c}}\left(\tilde {\boldsymbol{e}}\right){\tilde {\boldsymbol{e}}}_{\mathrm{c}}={\boldsymbol{Q}}_{\mathrm{c}} $

(29)

式中： ${\boldsymbol{M}}_{1\mathrm{c}}$ 、 ${\boldsymbol{M}}_{2\mathrm{c}}$ 为做了不同变化后的质量矩阵， $ \dot{\theta } $ 为旋转角速度， $\tilde {\boldsymbol{e}}$ 为浮动坐标向量， ${\tilde {\boldsymbol{e}}}_{\mathrm{c}}={\boldsymbol{B}}_{\mathrm{c}}\tilde {\boldsymbol{e}}$ 为系统约束后的浮动坐标向量， ${\ddot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}$ 为加速度， ${\dot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}$ 为速度。

$ {\boldsymbol{M}}_{1\mathrm{c}}=\rho {{\boldsymbol{B}}_{\rm{c}}^{\rm{T}}}\left({\int }_{0}^{l}\left(\iint {\boldsymbol{S}}_{\theta }\mathrm{d}A\right)\mathrm{d}x\right){\boldsymbol{B}}_{\mathrm{c}} $

(30)

$ {\boldsymbol{S}}_{\theta }={\left({\boldsymbol{U}}_{\mathrm{C}\mathrm{T}}\right)}^{\mathrm{T}}{\boldsymbol{S}}^{\mathrm{T}}\boldsymbol{S}{{\boldsymbol{U}}_{{\mathrm{C}\mathrm{T}},{\theta}}} \qquad \qquad\;\;$

(31)

$\;{{\boldsymbol{M}}}_{2{\rm{c}}}=\rho {{\boldsymbol{B}}_{\rm{c}}^{\rm{T}}}\left({\int }_{0}^{l}\left(\iint {\boldsymbol{S}}_{\theta \theta }\mathrm{d}A\right)\mathrm{d}x\right){\boldsymbol{B}}_{\mathrm{c}} $

(32)

$ {\boldsymbol{S}}_{\theta \theta }={\left({\boldsymbol{U}}_{\mathrm{C}\mathrm{T}}\right)}^{\mathrm{T}}{\boldsymbol{S}}^{\mathrm{T}}\boldsymbol{S}{{\boldsymbol{U}}_{\mathrm{C}\mathrm{T},\theta \theta }} \qquad \qquad $

(33)

${\boldsymbol{U}}_{\mathrm{C}\mathrm{T}}$ 为 $ 24\times 24 $ 的单元坐标变化矩阵， ${{\boldsymbol{U}}_{\mathrm{C}\mathrm{T},\theta }}_{}$ 和 ${{\boldsymbol{U}}_{\mathrm{C}\mathrm{T},\theta \theta }}_{}$ 分别为其对旋转角度的一阶导数和二阶导数。联合式(24)~(33)，计算得 ${\boldsymbol{M}}_{1\rm{c}}$ 和 $ {\mathit{M}}_{2\mathrm{c}} $ 为常数矩阵^[17]。

2.4 旋转梁动力学方程线性化

为了分析柔性旋转梁的频率特性及其模态变化，采用摄动理论对旋转动力学方程式(29)分别进行 ${\ddot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}$ 、 ${\dot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}$ 和 ${\tilde {\boldsymbol{e}}}_{\mathrm{c}}$ 的线性化。获得旋转梁振动方程为

$ {\boldsymbol{M}}_{\mathrm{c}}{\text{δ}} {\ddot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}+2\dot{\theta }{\boldsymbol{M}}_{1\mathrm{c}}{\text{δ}} {\dot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}+ ({\mathit{K}}_{\mathrm{T}1\mathrm{c}}({\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}})-{\dot{\theta }}^{2}{\boldsymbol{M}}_{2\mathrm{c}}){\text{δ}} {\tilde {\boldsymbol{e}}}_{\mathrm{c}}=0 $

(34)

式中： ${\boldsymbol{K}}_{\mathrm{T}}$ 为切线刚度矩阵, ${\boldsymbol{K}}_{\mathrm{T}1}$ 为切线刚度矩阵的线性部分， ${\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}$ 为当前旋转振动梁的平衡位置， ${\boldsymbol{K}}_{\mathrm{T}1\mathrm{c}}({\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}})$ 表示约束后的 ${\boldsymbol{K}}_{\mathrm{T}1}({\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}})$ ，约束方式与前一节相同。

$ {\boldsymbol{K}}_{\mathrm{T}1}\left({\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}\right)=\frac{\partial ({\boldsymbol{K}}_{\mathrm{t}\mathrm{m}}({\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}))}{\partial {\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}}{\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}+{\boldsymbol{K}}_{\mathrm{t}\mathrm{m}}({\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}) $

(35)

将 ${\ddot{\tilde {\boldsymbol{e}}}}_{\mathrm{c}}=0$ ， ${\tilde {\boldsymbol{e}}}_{\mathrm{c}}=0$ 代入式(29)，通过牛顿迭代法求解式(36)得到 ${\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}$ 。

$ -{\dot{\theta }}^{2}{\boldsymbol{M}}_{2\mathrm{c}}{\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}+{\boldsymbol{K}}_{\mathrm{c}}({\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}){\tilde {\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}-R{\dot{\theta }}^{2}{\boldsymbol{M}}_{2\mathrm{c}}\boldsymbol{C}=0 $

(36)

将 ${\text{δ}} {\tilde {\boldsymbol{e}}}_{\mathrm{c}}={{\boldsymbol{e}}}^{\mathrm{j}\omega t}\boldsymbol{Z}$ 代入式(34)，得到特征值方程

$ \left(\mathrm{j}\omega \left[\begin{array}{cc}{\boldsymbol{K}}_{\mathrm{T}}& 0\\ 0& {\boldsymbol{M}}_{\mathrm{c}}\end{array}\right]+\left[\begin{array}{cc}0& -{\boldsymbol{K}}_{\mathrm{T}}\\ {\boldsymbol{K}}_{\mathrm{T}}& 2\dot{\theta }{\boldsymbol{M}}_{1\mathrm{c}}\end{array}\right]\right)\boldsymbol{Z}=0 $

(37)

式中： $ \mathrm{j} $ 为虚数，t为摆动时间， $ \omega $ 为固有频率， $\boldsymbol{Z}$ 为复杂列向量常量， ${\boldsymbol{K}}_{\mathrm{T}}$ 为切线刚度矩阵。

$ {\boldsymbol{K}}_{\mathrm{T}}={\boldsymbol{K}}_{\mathrm{T}1\mathrm{c}}\left({\tilde{\boldsymbol{e}}}_{\mathrm{e}\mathrm{p}}\right)-{\dot{\theta }}^{2}{\boldsymbol{M}}_{2\mathrm{c}} $

(38)

3 仿真案例 3.1 自由单摆动力学仿真

以自由单摆仅受重力作用下的动力学仿真为例，验证本文方法所建梁模型的精确性。此时外力 $\boldsymbol{F}$ 为重力，是梁所受的体积力。

$ \boldsymbol{F}=\left[\begin{array}{ccc}0& -\rho {g}& 0\end{array}\right] $

(39)

根据形函数 $\boldsymbol{S}$ 和式(17)、式(39)可得到单元广义分布重力的表达式

$ {\boldsymbol{Q}}_{\mathrm{e}\mathrm{l}\mathrm{e}}^{\mathrm{T}}={\int }_{{V}}\left[\begin{array}{ccc}0& -\rho {g}& 0\end{array}\right]\boldsymbol{S}\mathrm{d}V $

(40)

式(39)~(40)中： $ \rho $ 为柔性梁密度， $ {g} $ 为重力加速度。

自由单摆如图4所示，柔性梁截面高为 $ h $ ，截面宽为 $ w $ ，总长为 $ L $ ，质量为m。

本文根据文献[12]中的算例进行了仿真实验1：弹性模量 $ E $ 为0.8 MPa，泊松比 $ \nu $ 为0.3，截面高 $ h $ 为0.02 m，截面宽 $ w $ 为0.08 m，总长 $ L $ 为1.2 m，密度 $ \rho $ 为5 540 $ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $ ，将单摆分成12单元。图5为仿真实验1中柔性单摆不同时间下的位置及形状，表1为仿真实验1与文献[12]的算例在t=0.3 s时的实验结果对比，证明了本文所采用的模型具有较高的准确性。

本文也对文献[18]中的算例进行了仿真实验2，其弹性模量 $ E $ 为2 MPa，泊松比 $ \nu $ 为0.3，截面高 $ h $ 为0.01 m，截面宽 $ w $ 为0.01 m，总长 $ L $ 为1.0 m，密度 $ \rho $ 为7200 $ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $ 。图6(a)为仿真实验2中柔性单摆不同时间下的位置及形状，表2为仿真实验2和文献[18]的算例在t=0.4 s时的实验结果对比，两者结果基本一致，进一步验证了本文方法的准确性。图6(b)为单摆能量曲线图，其中纵坐标为能量E_whole，横坐标为时间t，E_k为动能，E_{p, ela}为弹性变形能，E_{p, g}为重力势能，E_total为总能量。从图中可以看出整个仿真过程中机械能始终守恒。

3.2 非线性柔性梁固有频率的计算 3.2.1 柔性旋转梁固有频率的计算

为了验证旋转梁动力特性仿真精度，本文计算不同转速下的无量纲固有频率，并与文献[19]的计算结果进行对比。 $ {\omega }_{i} $ 是指不同阶数的固有频率，无量纲固有频率用不同阶数的固有频率与0阶的固有频率的比值表示，即 $ {\omega }_{i}/{\omega }_{0} $ 。

$ {\omega }_{0}=\sqrt{\frac{E\varGamma }{m{L}^{3}}} $

(41)

$ \varOmega =\gamma {\omega }_{0} $

(42)

式(41)~(42)中： $ E $ 为柔性梁的弹性模量， $ \varGamma $ 为截面惯性矩， $ m $ 为柔性梁的质量， $ L $ 为柔性梁的总长度， $ {\varOmega } $ 为无量纲旋转速度， $ \mathrm{\gamma } $ 为转速。

不同转速下无量纲固有频率前三阶计算结果如表3所示，两者结果基本一致，验证了旋转梁动力特性仿真精度。

3.2.2 不同转速下的各阶振型

梁模型的空间位置如图4所示，设置弹性模量 $ E $ 为 $ 5\times {10}^{8} $ Pa，泊松比 $ \nu $ 为0.3，截面高 $ h $ 为0.02 m，截面宽 $ w $ 为0.2 m，总长 $ L $ 为8.0 m，密度 $ \rho $ 为2677.67 $ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $ ，将单摆分成10单元进行仿真。绘制了不同转速下的前6阶模态振型，并将各阶振型投影到旋转平面(xy平面)或垂直于旋转平面的平面(xz平面)。

如图7~9所示，z方向模态振型从静止时的第3阶振型，变为了图8中转速为10 r/min的第2、6阶振型以及图9中转速为100 r/min的第2、4、6阶振型。随着转速的增加，前6阶中垂直于旋转平面(z方向)的振型所对应的频率会逐渐转换为低阶频率。

图10中，转速为140 r/min，y方向出现了3阶的振型，z方向出现了2阶振型，第6阶为扭转振型；图11中，转速为300 r/min，y方向出现了2阶的振型，z方向出现了2阶，第4阶为扭转振型，第6阶为拉伸振型。从此可见，随着转速的增加，扭转和拉伸模态振型所对应的频率会逐渐转换为低阶频率。

4 结论

本文在绝对节点坐标法基础上，运用连续介质理论推导了三维柔性梁的动力学方程，并对不同文献中的单摆案例进行了仿真，验证了本文方法的有效性。然后对动力学方程进行了线性化处理，计算了不同转速下柔性旋转梁的无量纲固有频率，并分析了非线性柔性旋转梁在不同转速下的模态振型变化特性。通过三维柔性梁动力学建模与动力特性的研究可得到以下结论。(1) 本文所建非线性梁模型能精确分析三维柔性梁的时域动力学响应；(2) 带浮动坐标的绝对节点坐标建模方法，能有效分析柔性旋转柔性梁的频域振动特性；(3) 三维柔性旋转梁随着转速的增加出现动力刚化现象，且垂直旋转平面的模态振型以及扭转和拉伸模态振型所对应的频率会转化为低阶频率。