0 引言

柔性旋转梁是空间机械臂、涡轮叶片、直升机旋翼及卫星天线等诸多工程结构的理论模型，其动力学特性对系统精度、性能与安全可靠性具有重要影响^[1-4].相关问题一直是国内外学者的研究热点，在动力学模型和建模理论方法方面，发展出了动力刚化模型^[5]、一次近似耦合模型^[6-8]、一次精确耦合模型^[9-11]、高次耦合模型^[12]、几何精确模型^[13]、浮动坐标法^[14-15]、绝对节点坐标法^[16]、随转坐标法^[17]、有限段方法^[18]、一致线性化方法^[19]等；在变形场离散与插值方法方面发展了Bezier插值法^[20-21]、B样条插值法^[22]及无网格插值法^[23-24]等；动力学方程分析求解方面，发展了假设模态法^[25]、有限差分法^[26]、有限元法等方法^[27]；取得了诸多研究成果.

为了准确进行旋转结构的设计，需要正确估计旋转系统结构及状态参数对其固有振动特性的影响. 在以旋转叶片、旋翼为背景的平面内柔性旋转梁问题研究中, 有学者发现系统频率、模态振型随旋转速度变化出现频率转向与振型转换现象^[28].在以含梁式柔性附件的航天器为背景的中心刚体-柔性旋转梁系统动力学特性研究中，蒋丽忠和洪嘉振^[29]指出大范围运动转速对频率和模态均有影响，杨辉等^[30-31]通过理论分析和实验验证发现定边界的模态与动边界的模态存在较大差异，指出刚柔耦合系统模态频率与刚体-柔性梁转动惯量比值有一定关系，方建士等^[32]发现随着转速的增加，旋转悬臂梁将出现频率转向与振型转换现象，Zhao等^[33]也同样发现刚柔耦合系统动力学特性会随着中心刚体转动惯量及转速的变化而变化.但上述研究在动力学建模时多认为旋转梁根部与中心刚体为刚性固接，而实际工程中由于连接件弹性变形^[34]、根部裂纹^[35-37]等因素，根部边界处不是理想的刚性连接，而应视为具有一定扭转刚度的铰接.此外，现有研究在动力学特性分析时多集中于旋转速度对系统固有频率与振型的影响，而对系统径长比、梁的长细比等因素的影响分析较少. 最后，现有离散方法往往需要根据不同的边界条件构造不同的形函数，处理过程复杂^[38].Chebyshev谱方法作为具有指数收敛特性的高精度数值方法，能够很好地解决各类线性、非线性梁动力学问题，结合投影矩阵法能够"一致"地处理固支、铰支及弹性支撑等不同边界问题，为解决上述问题提供了一种新的途径^[38-39].

本文基于Gauss-Lobattoe节点与Chebyshev多项式插值对柔性梁变形场进行离散，在此基础上考虑根部连接弹性，对平面内旋转柔性梁动力学特性进行研究，通过投影矩阵法施加固定及弹性连接边界条件，利用Chebyshev谱方法获得系统固有频率及模态振型的数值解，将仿真结果与有限元法、加权残余法等方法进行对比，验证方法的有效性；并对弹性连接刚度、角速度比率、系统径长比及梁的长细比等参数对系统固有频率及模态振型的影响进行分析.

1 旋转柔性梁动力学模型 1.1 系统动力学方程

考虑图 1所示平面内旋转柔性梁系统，中心刚体半径为$R$，在外力驱动下作角速度为$\Omega $的定轴转动，中心刚体转动惯量远大于柔性梁，考虑二者连接处的弹性将其等效为铰接加一扭簧连接，扭转弹簧系数为$k_{\rm c}$，梁未变形时的轴线延长线过转轴. 柔性梁长为$L$，密度为$\rho $，截面面积为$A$,截面惯性矩为$I$，弹性模量为$E$.通过中心刚体转轴$O$建立惯性坐标系$OXYZ$，在柔性梁与中心刚体连接点建立浮动坐标系$oxyz$.

基于von-Kavman非线性应变理论，采用Euler-Bernoulli梁模型，平面内旋转柔性梁系统自由振动动力学方程可表示为^{[28, 40]}

$\left. {\matrix{ {\rho A\left( {\ddot u - 2\Omega \dot v - {\Omega ^2}u - \dot \Omega v} \right) - EAu'' = } \cr {\rho A{\Omega ^2}\left( {R + x} \right)} \cr {\rho A\left( {\ddot v + 2\Omega \dot u + \dot \Omega u - {\Omega ^2}v} \right) - EA{{\left( {u'v'} \right)}^\prime } + EIv'''' = } \cr { - \rho A\dot \Omega \left( {R + x} \right)} \cr } } \right\}$

(1)

式中，$u,v$分别为点$P$在浮动坐标系$oxyz$中$x,y$方向的位移，右上标符号" $'$ "表示对空间变量$x$的导数，顶标"$ \cdot $"表示对时间变量$t$的导数.从式中可以看出柔性梁横向弯曲变形与轴向拉伸变形相互耦合. 对于刚性连接，有边界条件

$ \left. \!\! u\left( {0,t} \right) = 0 ,v\left( {0,t} \right) = 0 ,\ {v}'\left( {0,t} \right) = 0 \\ {u}'\left( {L,t} \right) = 0 ,{v}"\left( {L,t} \right) = 0 ,\ {v}"'\left( {L,t} \right) = 0 \!\!\right\} $

(2)

对于弹性连接，有边界条件

$ \left. \!\! u\left( {0,t} \right) = 0 ,v\left( {0,t} \right) = 0 ,\ \left. {EI{v}" - k_{\rm c} {v}'} \right|_{x = 0} = 0 \\ {u}'\left( {L,t} \right) = 0 ,\ {v}"\left( {L,t} \right) = 0 ,\ {v}"'\left( {L,t} \right) = 0 \!\!\right\}$

(3)

考虑大范围匀速运动情况，此时$\dot {\Omega } =0$，基于微扰法，柔性梁轴向变形可看作平衡状态变形量与扰动变形量之和，即

$ u = u_0 + \Delta u $

(4)

平衡状态变形量$u_0 $对时间的导数为零，将式(4)代入式(1)，可解得稳态轴向变形量

$ u_0 = \dfrac{R\beta \cos \left[{\beta \left( {L - x} \right)} \right] + \sin \left( {\beta x} \right)}{\beta \cos \left( {\beta L} \right)} - \left( {R + x} \right) $

(5)

式中$\beta = \Omega \sqrt {\rho/E} $，将式(5)代入式(1)，得到系统在平衡位置附近自由振动运动方程

$ \left. \begin{array}{l} \ddot {u} - 2\Omega \dot {v} - \Omega ^2u - \chi {u}" = 0 \\ \ddot {v} + 2\Omega \dot {u} - \Omega ^2v - \chi \left( {{u}"_0 {v}' + {u}'_0 {v}"} \right) + \tau {{v}"}" = 0 \end{array} \right\}$

(6)

式中$\chi = E/ \rho $，$\tau = {EI}/ (\rho A)$.

1.2 Chebyshev谱方法

Chebyshev谱方法以Chebyshev多项式插值为基础的方法.本文采用第一类Chebyshev多项式，它是一类递归正交多项式，其第$k$项可表示为^[38]

$ T_k (x) = \cos (k\arccos x) ,\quad k = 1,2,3,\cdots$

(7)

且满足递推关系

$ T_k (x) = 2xT_{k - 1} (x) - T_{k - 2} (x)$

(8)

式中$x \in [- 1,1]$. 对于$x \in \left[{0,L} \right]$的情况可采取变换 $\xi = {(2x - L)} / L$，此时式中采用 $T_k \left( \xi \right) = T_k \left( {\xi \left( x \right)} \right)$替代即可. 对于任意在$x \in [-1,1]$上连续可导且平方可积的函数$y(x)$采用$N -1$阶Chebyshev多项式逼近

$ y_N (x) = \sum_{k = 0}^{N - 1} a_k T_k (x)$

(9)

若采用 Gauss-Lobatto 节点$x_k = \cos \left[{{i\pi }/{\left( {N - 1} \right)}} \right]$，$i = 0,1,\cdots,N - 1$，节点函数值为$y_k \left( {x_k } \right)$，则有节点值列向量${ y}$与Chebyshev多项式系数列向量${ a}$关系

${ y} = { \Gamma }_{\rm B} { a}$

(10)

$ { \Gamma }_{\rm B} = \left[\!\! \begin{array}{cccc} {T_0 (x_0 )} & {T_1 (x_0 )} & \cdots & {T_{N - 1} (x_0 )} \\ {T_0 (x_1 )} & {T_1 (x_1 )} & \cdots & {T_{N - 1} (x_1 )} \\ \vdots & \vdots & \ddots & \vdots \\ {T_0 (x_{N - 1} )} & {T_0 (x_{N - 1} )} & \cdots &{T_0 (x_{N - 1} )} \end{array}\!\! \right]$

(11)

并可得函数$y ( x )$导数的多项式逼近

$ {y}'_N (x) = \sum_{k = 0}^{N - 1} a_k {T}'_k (x) = \sum_{k = 0}^{N - 1} b_k T_k (x) $

(12)

多项式系数$b_k $组成$N$维系数列向量${ b} $，并有

$ { b} = { D a}$

(13)

式中${ D }$为Chebyshev微分矩阵. 同理可得函数$y(x)$的$n$阶导数逼近

$ { y }^{(n)} = { \Gamma }_{\rm B} { D }^n{ \Gamma }_{\rm B}^{ - 1} { y } = { Q }_n { y }$

(14)

${ Q }_n $为$n$阶Chebyshev微分矩阵. 函数$y ( x )$的积分可表示为

$ \int_{l_1 }^{l_2 } {y_N (x)} d x = \sum_{k = 0}^N {a_k } \int_{l_1 }^{l_2 } {T_k (x)} d x = { v }^{\rm T}{ a }_N$

(15)

式中${ v }$为$N$维列向量，称作Chebyshev积分向量. 在区间$\left[{l_1 ,l_2 } \right]$平方可积的函数$f ( x )$和 $g ( x )$的内积可表示为

$ \int_{l_1 }^{l_2 } {f(x)g(x)d x = { f}^{\rm T}} { V g}$

(16)

${ V}$为$N\times N$ Chebyshev内积矩阵，${ f}$和${g}$分别为函数$f\left( x \right)$和$g\left( x \right)$的Gauss-Lobattoa节点值列向量. 同理，函数$f\left( x \right)$、$g\left( x \right)$关于函数$h\left( x \right)$的加权内积可表示为

$ \int_{l_1 }^{l_2 } {f\left( x \right)} g\left( x \right)h\left( x \right)d x = { f }^{\rm T}{ V }_h { g }$

(17)

式中${ V}_h $是关于函数$h\left( x \right)$的$N\times N$ Chebyshev加权内积矩阵.

由Chebyshev多项式理论可知，在已知函数$N$个Gauss-Lobatto 节点值时，可以用Chebyshev多项式逼近函数及其各阶导数、积分及函数间的内积.

1.3 动力学方程离散

采用Chebyshev多项式对变形场进行逼近，则变形量$u\left( {x,t}\right)$和$v\left( {x,t} \right)$可用$N$阶Chebyshev多项式逼近

$u\left( {x,t} \right) \approx \sum\limits_{k = 0}^{N - 1} {{a_{uk}}\left( t \right)} {T_k}\left( x \right)$

(18)

$ v\left( {x,t} \right) \approx \sum_{k = 0}^{N - 1} {a_{vk} \left( t \right)} T_k (x) $

(19)

由于采用有限项多项式逼近，将式(18)和式(19)代入式(6)后，等式右端并不为0，将存在残差，根据伽辽金方法，选择可用Chebyshev多项式表示的函数$\theta_u (x) $,$\theta _v (x) $作为权函数，则应有

${\matrix{ {\int_0^L {{\theta _u}\left( {\ddot u - 2\Omega \dot v - {\Omega ^2}u - \chi u''} \right)} dx = 0} \cr {\int_0^L {{\theta _v}\left[ {\ddot v + 2\Omega \dot u - {\Omega ^2}v - \chi \left( {{{u''}_0}v' + {{u'}_0}v''} \right) + \tau v''''} \right]dx} } \cr } }$

(20)

采用Gauss-Lobatto节点，节点值构成$N$维列向量${ u }$和${ v }$，根据式(14)和式(16) 得

$\left. {\matrix{ {\theta _u^{\rm{T}}V\left( {\ddot u - 2\Omega \dot v - {\Omega ^2}u - \chi {Q_2}u} \right) = {\bf{0}}} \cr {\theta _v^{\rm{T}}\left[ {V\left( {\ddot v + 2\Omega \dot u - {\Omega ^2}v + \tau {Q_4}v} \right) - } \right.} \cr {\left. {\chi \left( {{V_{\matrix{ {u''} \cr {_0} \cr } }}{Q_1}v + {V_{{{u'}_0}}}{Q_2}v} \right)} \right] = {\bf{0}}} \cr } } \right\}$

(21)

表示为矩阵形式为

$ { \theta }^{\rm T}\left( {{ M}\ddot{ q } + { G}\dot { q} + { K}{ q} } \right) = {\bf 0}$

(22)

式中

${ q } = \left\{ \!\! \begin{array}{c} { u } \\ { v } \end{array} \!\! \right\} , \quad { \theta } = \left[\!\! \begin{array}{c} {{ \theta }_u } \\ {{ \theta }_v } \end{array} \!\! \right]$

(23)

${ M } = \left[\!\! \begin{array}{cc} { V } & {\bf 0} \\ { 0 } & { V} \end{array} \!\! \right] , \quad { G } = \left[\!\! \begin{array}{cc} {\bf 0 } & { - 2\Omega { V }} \\ {2\Omega { V }} & {\bf 0} \end{array} \!\! \right]$

(24)

$\eqalign{ & K = \cr & \left[ {\matrix{ { - {\Omega ^2}V - \chi V{Q_2}} & {\bf{0}} \cr {\bf{0}} & { - {\Omega ^2}V + \tau V{Q_4} - \chi \left( {{V_{{{u''}_0}}}{Q_1} + {V_{{{u'}_0}}}{Q_2}} \right)} \cr } } \right] \cr} $

(25)

离散后，式(2)中刚性连接边界条件可表示为

${ B q } = { b}$

(26)

${ B } = \left[\!\! \begin{array}{cc} { e }_1 & {\bf 0} \\ {\bf 0 } & { e }_1 \\ {\bf 0 } & { e }_1 { Q }_1 \\ { e }_N { Q }_1 & {\bf 0 } \\ {\bf 0} & { e }_N { Q }_2 \\ {\bf 0} & { e }_N { Q }_3 \end{array} \!\! \right] , \quad { b } = \left[\!\! \begin{array}{c} { 0} \\ { 0} \\ { 0} \\ { 0} \\ { 0} \\ { 0} \end{array} \!\! \right]$

(27)

对于式(3) 中弹性边界条件，上式中${ B} $则变为

$ { B } = \left[\!\! \begin{array}{cc} { e }_1 & {\bf 0} \\ {\bf 0} & { e }_1 \\ {\bf 0} & EI{ e }_1 { Q }_2 - k_c { E }_1 { Q }_1 \\ { e }_N { Q }_1 & {\bf 0} \\ {\bf 0} & { e }_N { Q }_2 \\ {\bf 0} & { e }_N { Q }_3 \end{array} \!\! \right]$

(28)

上述两式中${ B }$为$6\times 2N$边界条件矩阵，${ e }_i $为第$i$个元素为1，其他元素均为0的$N$维行向量.对于满足式(26)的解${ q }$可表示为^[39]

$ { q} = { P}\tilde { q} + { R b}$

(29)

其中$\tilde { q}$为$2N - 6$维列向量，${ P}$为$2N\times \left( {2N - 6} \right)$矩阵，是从$2N$维${q}$空间到其$(2N-6)$维$\tilde{ q}$子空间的投影矩阵，${ P}$和${ R}$可以通过对${B}$的奇异值分解获得. 采用投影矩阵法在离散时无需考虑边界条件问题，不同边界条件下离散过程具有一致性. 由于式(21) 中试函数$\theta (x) $亦需满足边界条件，则有

$ { \theta } = { P}\tilde {\theta } + { R b}$

(30)

将式(29)和式 (30)代入式(22)，并考虑到${ b }$为零向量，得

$ \tilde{ \theta }^{\rm T}{ P }^{\rm T}\left( {{ M}{ P}\ddot {\tilde { q}} + { G P}\dot {\tilde { q}} + { K P}\tilde { q}} \right) = {\bf 0}$

(31)

由于试函数${ \theta }_u (x) $与 ${ \theta }_v (x) $具有任意性，可得到

$\tilde{ M}\ddot {\tilde { q}} + \tilde{ G}\dot {\tilde { q} } + \tilde { K}\tilde{ q} = {\bf 0}$

(32)

$ \tilde{ M} = { P }^{\rm T}{ M P } ,\quad \tilde{ G} = { P }^{\rm T}{ G P } ,\quad \tilde{ K } = { P }^{\rm T}{ K P }$

(33)

根据式(32)采用状态空间法可得到系统模态频率与振型，令

$ { z } = \left[ \tilde{ q }^{\rm T} \dot {\tilde { q} }^{\rm T} \right]^{\rm T}$

(34)

式(32)转换为状态方程

$ { C}\dot { z} + { D z} = {\bf 0}$

(35)

式中

$ { C} = \left[\!\! \begin{array}{cc} \tilde{ G} & \tilde{ M} \\ \tilde{ M} & {\bf 0} \end{array}\!\! \right] , \quad { D} = \left[\!\! \begin{array}{cc} \tilde{ K} & {\bf 0} \\ {\bf 0} & -\tilde{ M} \\ \end{array} \!\! \right]$

(36)

设方程的解为${ z} = { Z}{\rm e}^{st}$ 将其代入式(35)得复特征值问题

$ \left( {s{ C} + { D} } \right){ Z} = {\bf 0}$

(37)

式中，$s$为复特征值，${ Z}$为复模态振型.

2 动力学特性分析 2.1 弹性连接刚度影响

引入以下无量纲变量

$\eqalign{ & \cr & \left. {\matrix{ {\delta = {R \over L},\quad \alpha = \sqrt {{{A{L^2}} \over I}} ,\quad \gamma = \Omega {L^2}\sqrt {{{\rho A} \over {EI}}} } \cr {f = \omega {L^2}\sqrt {{{\rho A} \over {EI}}} ,\quad \kappa = {{{k_{\rm{c}}}L} \over {EI}}} \cr } } \right\} \cr} $

(38)

式中，$\delta $为系统径长比，$\alpha $为梁的长细比，$\gamma$为无量纲角速度比率，$f$为无量纲固有频率，$\omega $为系统固有频率，$\kappa $为无量纲连接刚度.

为验证Chebyshev谱方法的有效性，考虑刚性连接边界条件下$\delta = \gamma =0$,$\alpha = 70$时的情况，当$N = 12$时，前四阶无量纲固有频率分别收敛至3.156 0,22.034 5,61.697 1, 109.955 7,120.902 0，与文献 ^[41] 中给出的精确解3.516 0,22.034 5,61.697 1,109.955 7,120.901 9吻合很好，而有限元方法中，当单元个数为100时，其获得的解分别为3.516 0,22.034 5,61.697 2,109.956 9,120.901 9^[41]； 另当$\delta = 0$,$\gamma = 10$,$\alpha = 70$时，当$N = 13$时前四阶无量纲固有频率分别收敛至5.055 9, 32.124 1,74.008 0,111.316 1,134.573 4，如图 2(a)所示，文献 ^[28] 中前四阶固有频率分别收敛至5.056 2, 32.124 7,74.009 4,111.316 1,134.575 6，二者吻合较好；通过与文献 ^{[28, 41]} 的比较验证了Chebyshev方法的有效性. 同时从图 2可以看出，Chebyshev谱方法具有快速收敛性，低阶频率比高阶频率随Chebyshev多项式阶数收敛更快，采用较少阶数多项式便可获得高精度解. 弹性连接边界条件下，$\delta = 0$,$\gamma = 10$,$\alpha = 70$, $\kappa = 100$时，计算结果收敛性如图 2(b) 所示，当$N = 12$ 时，前四阶无量纲固有频率分别收敛至3.767 6,29.920 1, 68.629 8,111.314 9,125.249 5，可以看出计算结果仍然具有快速收敛特性；同时与刚性连接条件相比，各阶固有频率均下降.

图 3为$\delta = 0$,$\gamma = 2$,$\alpha =70$时，弹性连接刚度对系统模态频率与振型的影响，点$A_1 $和$A_2 $模态振型如图 3 (b) 所示，以轴向变形量$u$为主(以下简称拉伸模态)，随着无量纲连接刚度$\kappa$的增大拉伸模态频率保持不变，其轴向变形振型不变，但横向变形增大；点$B_1 $,$B_2 $振型如图 3 (c)所示，以横向变形$v$为主(以下简称弯曲模态)，显然两点为同一弯曲模态，但其固有频率随着连接刚度$\kappa$的增加而呈现阶梯上升状变化，增大并超过原高一阶拉伸模态频率，如图 3 (a)所示，即发生了频率转向与振型转换，此外说明，为保持系统动力学特性稳定，柔性梁根部连接刚度应保持图中阈值$\kappa _{\rm c} $之上，否则系统弯曲模态频率会随着连接刚度的下降而迅速下降，尤其是高阶频率.

2.2 系统参数影响

图 4为$\delta = 0$，$\alpha = 70$，$\kappa = 10^2$时，无量纲角速度比率对系统模态频率与振型的影响，图 4 (a) 结果与文献 ^[28] 一致，为分析图中出现传递性频率转向的原因，选择图中7个点分析其模态振型，图 4 (b)为$A_1 $和$A_2$模态振型，显然其为同一拉伸模态(第一阶拉伸模态)，随着角速度比率的增加其轴向变形振型几乎不变，其频率增加较弯曲模态慢，因此同样会出现低阶弯曲模态频率增大并超过上一阶拉伸模态频率；同样从图 4 (c) 中可以看出$B_1 $和$B_2$均为第二阶拉伸模态；从而说明出现频率转向的原因是弯曲模态频率对角速度比率"敏感"，随其增大而快速增大，而拉伸模态频率随其增大增加较小，因此低价弯曲模态频率能在增大的过程中超过高一阶拉伸模态频率，从而出现频率转向与振型转换现象.

图 5为$\gamma = 10$,$\alpha = 70$,$\kappa =10^3$时系统径长比对无量纲固有频率的影响，由于拉伸模态频率随径长比增大而不变，弯曲模态频率则随之增大，出现拉伸模态频率轨迹与弯曲模态频率轨迹的交叉，发生频率转向与振型转换.

图 6为$\gamma = 10$,$\delta = 0.1$,$\kappa = 10^3$梁的长细比对无量纲固有频率影响，从图中可以看出，弯曲模态频率不随长细比增大而变化，而拉伸模态频率则线性增加，从而出现了低阶拉伸模态频率增大并超过上一阶弯曲模态频率的现象，即频率转向与振型转换. 这说明系统结构参数的变化亦将导致频率转向与振型转换的发生.

3 结论

(1) Chebyshev谱方法应用于柔性旋转梁动力学分析结果与有限元及加权残余法结果基本一致，采用较少的多项式阶数即可获得与有限元方法相当的精度，采用投影矩阵方法可以方便地用同样形式处理固定及弹性边界条件.

(2) 出现频率转向与振型转换的原因是系统弯曲模态、拉伸模态随系统参数变化规律不一致，变化快的低价模态频率增大并超过原高阶模态频率.

(3) 系统弯曲模态频率随连接刚度呈阶梯状变化，拉伸模态频率不随连接刚度变化而变化，二者频率变化曲线交叉，出现频率转向与振型转换现象，为保证系统动力学特性稳定，连接刚度应大于关键刚度$\kappa _{\rm c}$.

(4) 系统弯曲模态频率随系统径长比增大而增大，拉伸模态频率不随径长比变化而变化，二者频率变化曲线交叉，出现频率转向与振型转换现象.

(5) 拉伸模态频率随梁的长细比的增大而增大，弯曲模态频率保持不变，二者频率变化曲线交叉，出现频率转向与振型转换现象.