CN110457740B - 一种机械结构在基础激励下的响应特性分析方法 - Google Patents

Abstract

本发明属于旋转梁动力特性研究技术领域，尤其涉及一种机械结构在基础激励下的响应特性分析方法。该方法包括以下步骤：将机械结构的动力学分析退化成对受到基础激励的旋转悬臂梁模型的动态特性的分析；通过汉密尔顿原理和有限元方法将基础激励下旋转悬臂梁模型离散成一系列的铁木辛柯梁段；通过固有频率和基础激励响应的比较证明旋转悬臂梁模型的正确性。该方法采用了较为简化的数学模型，并且此分析方法具有较高的精确性。

Description

一种机械结构在基础激励下的响应特性分析方法

技术领域

本发明属于旋转梁动力特性研究技术领域，尤其涉及一种机械结构在基 础激励下的响应特性分析方法。

背景技术

很多学者对旋转梁的动力特性进行了大量的研究，通过对此研究能 够对机械设计过程的中选材、结构的设置具有指导意义，但就我们所知， 目前的研究大多局限于模态特性的分析，而对扭转振动的振动响应，特 别是在外力作用下的扭转振动的振动响应研究较少。另外，采用有限元 方法并考虑旋转效应(科氏力效应，刚化效应，软化效应和角加速度效应) 时，很少有人考虑梁的弯-摆-轴-扭耦合影响，特别是在模拟梁的扭转振 动时计算精度较差。

发明内容

(一)要解决的技术问题

针对现存的动力特征分析精确性差的技术问题，本发明提供一种机 械结构在基础激励下的响应特性分析方法。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

一种机械结构在基础激励下的响应特性分析方法，其包括以下步 骤：

将机械结构的动力学分析退化成对受到基础激励的旋转悬臂梁模型 的动态特性的分析；

通过汉密尔顿原理和有限元方法将基础激励下旋转悬臂梁模型离散 成一系列的铁木辛柯梁段；

通过固有频率和基础激励响应的比较证明旋转悬臂梁模型的正确 性。

(三)有益效果

本发明的有益效果是：本发明提供的一种机械结构在基础激励下的 响应特性分析方法，采用了较为简化的数学模型，并且此分析方法具有 较高的精确性。

附图说明

图1为旋转双锥度悬臂预扭梁的有限元离散：其中，图1(a)某航空 发动机一级风扇叶片，图1(b)梁的离散模型，图1(c)欧拉角的三个坐 标系图；

图2为γ(L)对扭转刚度和固有频率的影响：图2(a)对扭转刚度的影 响，图2(b)对固有频率的影响；

图3为旋转双锥度悬臂预扭梁有限元模型；

图4为宽度锥度τ_b值对固有频率的影响；

图5为在γ(L)＝45°，β₀＝0°，τ_b＝0.3和n＝10000r/min时的振动响应： (a)X向，(b)Y向，(c)Z向，(d)θ向；

图6为安装角对一阶动频的影响：(a)Y_I，(b)Z_I，(c)θ_I；

图7为预扭角对一阶动频的影响：(a)Y_I，(b)Z_I，(c)θ_I；

图8为幅频响应曲线：(a)X，(b)Y，(c)Z，(d)θ；

图9为一阶弯曲动频曲线；

图10为不同安装角和转速下的常值分量：(a)X方向，(b)θ方向；

图11θ向变形云图(n＝20000rev/min)：(a)β₀＝45°，γ(L)＝0°，(b)β₀＝-45°，γ(L)＝0°；

图12为幅频响应：(a)X，(b)Y，(c)Z，(d)θ；

图13为一阶弯曲动频曲线；

图14为不同安装角和转速下的常值分量：(a)X，(b)θ；

图15为θ向变形云图(n＝20000rev/min)：(a)β₀＝0°，γ(L)＝45°，(b)β₀＝0°，γ(L)＝-45°；

图16为A.1单元轴向变形；

图17为A.2单元扭转变形；

图18为A.3微元体在弯曲方向的受力及剪应变关系，其中，图18 (a)为微元体力平衡关系；图18(b)为剪应力关系；图18(c)为力矩 -正应力关系；

图19为A.4单元弯曲变形；

图20为A.5微元体在弯曲方向的受力及剪应变关系，其中，图20 (a)为微元体力平衡关系；图20(b)为剪应力关系；图20(c)为力矩 -正应力关系；

图21为A.6单元摆动变形；

图22为C1矩形截面尺寸。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实 施方式，对本发明作详细描述。

本发明公开了一种机械结构在基础激励下的响应特性分析方法，其 包括以下步骤：

将机械结构的动力学分析退化成对受到基础激励的旋转悬臂梁模型 的动态特性的分析；

通过汉密尔顿原理和有限元方法将基础激励下旋转悬臂梁模型离散 成一系列的铁木辛柯梁段；

通过固有频率和基础激励响应的比较证明旋转悬臂梁模型的正确 性。

为了对上述机械结构在基础激励下的响应特性分析方法进行进一步 的说明，本实施方式提供了上述分析方法的具体步骤：

基础激励下叶盘的工程实例如图1(a)所示，它由两个一级风扇叶 片，刚性盘和弹性支撑组成。在图1(a)中，松脱叶片和未松脱叶片之间 的转角是180°。一旦装置运转起来并达到某一转速时，在离心力的作 用下两个一级风扇叶片中的一个将从刚性盘中松脱，这将导致叶盘系统 的不平衡。而这个不平衡的叶盘系统将会产生涡动，同时，未松脱叶片 的运动可以看做是叶盘涡动和叶片转动的合成，也就是说，未松脱的旋 转叶片的动力学分析可以退化成对受到基础激励的旋转悬臂叶片的动态 特性的分析。

基于上述的工程实例，本文简化并采用了相似的数学模型，即通过 汉密尔顿原理和有限元方法将基础激励下旋转态双锥度悬臂预扭梁离散 成一系列的铁木辛柯梁段，通常我们称之为铁木辛柯梁单元(图1(b))。 在这一节中，重点是推导铁木辛柯梁的运动偏微分方程和单元矩阵。需 要提出的是，旋转引起的科氏力，离心刚化，旋转软化，角加速度效应 在弯曲方向，摆动方向，轴向和扭转方向都有体现。最后，通过现有技 术中的固有频率和基础激励响应的比较证明了本文所提出的旋转双锥度 悬臂预扭梁模型的正确性。ANSYS中采用Beam188和Shell181单元。

2.1偏微分方程的建立

以第k个铁木辛柯梁单元为例，建立相应的弯扭轴摆耦合运动方程 (图1(b))。在图1(b)中，OXYZ，OX^rY^rZ^r，o^bx^by^bz^b和o^ex^ey^ez^e分别是整体坐标系，旋 转坐标系，局部坐标系和单元局部坐标系，α为旋转角。在图1(c)中，三 个坐标系中的欧拉角θ，φ和

用来表示在单元局部坐标系任意梁截面处的 角位移，β₀为叶片安装角，γ_k为第k个梁单元处的预扭角。从而，第k个梁 单元处的扭角β_k＝β₀+γ_k。R_k＝R_d+x_k+X_u ^r，其中R_k是o^e在旋转坐标系下的x坐标，R_d是叶盘半径，x_k是o^e在局部坐标系中的x坐标，而X_u ^r是X^r的基础激励。梁 截面在整体坐标系下任意点P的坐标可以表示为：

其中(x,y,z)是P点在ox₃y₃z₃坐标系(图1(c))的坐标。u，v和w 分别是梁上任意截面形心'o'沿着x^e-，y^e-和z^e-轴正方向的线位移；A₁， A₂，A₃，A₄和A₅是从ox₃y₃z₃到OXYZ的转换矩阵，它们可以表示为：

在公式(2)中，基于小变形假设采用一阶近似，即：

cos(φ)≈1，sin(φ)≈φ，cos(θ)≈1和sin(θ)≈θ。

根据公式(1)和公式(2)，第k个铁木辛柯梁单元的动能T可以表示为：

其中，ρ，l_k和A分别为密度，第k个梁单元的单元长度和截面面 积。符号上点表示对时间的导数。

第k个梁单元的势能U可以表示为：

其中，E和G分别是弹性模量和剪切模量；I_x，I_y和I_z分别是关于 轴ox₃，oy₃和oz₃的惯性矩；J是扭转惯性矩，按2.4节方法计算；κ_y和κ_z分别是关于y-和z-方向的剪切系数；

是离心 力；/>

是基础激励在x-方向产生的力；

是一个常值。另外，公式(4)中的求导符号表示对x的偏导 数。需要指出的是，在公式(4)中，下划线的一项表示离心刚化效应对扭 转变形的影响，而双下划线的一项表示预扭角对势能的影响。

视u，v，w，θ，φ和

为独立变量，并应用汉密尔顿能量方程/>

从而获得六个偏微分方程。对于x^e-轴向的位移u，有：

y^e-弯曲方向的位移v为：

z^e-摆动方向的位移w为：

x^e方向的角位移θ为：

y^e方向的角位移φ为：

z^e方向的角位移

为：

2.2铁木辛柯梁单元的单元矩阵

为了求解公式(5)-(10)，u-，v-，w-，θ-，φ-和

-方向的假设位移场函数及 相关表达式见于附录A。将公式(A.1)-(A.4)代入公式(5)-(10)，可以得到如下 的二阶常微分方程：

其中，M^e和

分别表示单元质量矩阵和结构刚度矩阵；G^e，/>

和/>

分别 是当角加速度/>

时的科氏力矩阵，离心刚化矩阵和旋转软化矩阵； />

是与角加速度

有关的刚度矩阵；δ^e和F^e分别表示节点位移和外 力向量。单元矩阵的细节见于附录B。

2.3单元矩阵的组集

2.2节中所有的单元矩阵(M^e，G^e，

δ^e和F^e)都是 从单元局部坐标系o^ex^ey^ez^e中得到的。然而，在组集这些单元矩阵以建立 旋转双锥度悬臂预扭梁时，必须统一相应的单元坐标系。在本文中，所 有的单元矩阵都转化到旋转坐标系OX^rY^rZ^r下进行组集，转换矩阵公式如 下：

其中，T是从OX^rY^rZ^r到o^ex^ey^ez^e坐标系的转换矩阵。δ^e和δ_r ^e分别是节点在o^ex^ey^ez^e和OX^rY^rZ^r坐标系中的位移向量。

基于公式(11)和公式(12)，旋转坐标系下的单元矩阵如

和F_r ^e可以这样计算：

那么此悬臂梁总的矩阵χ(χ＝M,G,K_e,K_c,K_s,K_acc,F)可以通过

得到，其中N是梁单元的总数。

2.4扭形直梁的扭转刚度

考虑到宽厚比(h/b，附录C中图22)的影响，附录C给出了计算矩形直 梁扭转惯量(J_s)的详细推导过程。这里直接给出它的计算公式：

另外，预扭角γ(L)的存在会产生额外的扭转惯量J_a：

其中，γ(L)是在L位置处梁的扭角(见图1)。

从公式(14)和公式(15)可得，J＝J_s+J_a。在图2中，本文得出的扭转 刚度GJ和第一阶扭转固有频率θ_I分别与文献[1]吻合很好，从而证明了 上述公式的正确性。

文献[1]C.William,Vibrations of pre-twisted cantilever blading,Proceedings of the Institution of Mechanical Engineers 173(1)(1959)343-374.

上述文献的中文翻译为：卡内基威廉扭形悬臂叶片的振动分析， 机械工程师学会会议录173(1)(1959)343-374。

2.5模型验证

根据2.1节-2.4节，旋转双锥度悬臂预扭梁有限元模型如图3。此 外，应该采取一些方法来验证本文模型。针对侧边锥度和预扭角对梁结 构动力特性的影响，许多学者提出了各种梁模型。本文采用文献[2]的 扭型悬臂梁来佐证本文的模型。材料和几何参数定义如下： E＝2.0684×10¹¹Pa，＝7861.1kg/m³，泊松比υ＝0.3，κ_y＝κ_z＝5/6，R_d＝0m， 梁长度L＝0.1524m，叶根宽度b₀＝0.0254m，叶根厚度h₀＝0.0045974m， 厚度锥度τ_h＝0.436464。为了验证模型的正确性，比较模型之间的固有频 率和振动响应。ANSYS有限元模型采用Beam188和Shell181单元，如图 4和图5所示。

其中，文献[2]J.S.Rao,Flexural vibration of pretwisted taperedcantilever blades,Journal of Engineering for Industry 94(1)(1972)343-346.

上述文献[2]的中文翻译为：J.S.拉奥扭锥形悬臂叶片的弯曲振动 工业工程杂志94(1)(1972)343-346。

在图4中，Y_I,Z_I,Y_II和Z_II是从本文提出的模型获得的，不同宽度锥度τ_b值下用Beam188，Shell181单元和参考文献得到的结果相互吻合，这样就证实 了本文的模型。

基于这几个模型，梁顶部节点在基础激励

下X-，Y-，Z-和 θ-向的振动响应如图5所示。其中，A₀＝1×10^-4m，/>

n＝10000rev/min。 需要指出的是，计算了17个激励周期的振动响应，每个周期分为60个 时间步长。此外，考虑瑞利阻尼C，表达式如下：

其中，ω₁＝2πf₁和ω₂＝2πf₂。[f₁,f₂]是此模型中较接近一、二阶固有频率的频率范 围，而ξ₁和ξ₂是模型的模态阻尼比。在本文中，f₁＝195.56Hz，f₂＝733.05Hz，ξ₁＝ξ₂＝0.03。

总的来说，各个模型X-，Y-和Z-方向的响应吻合很好(图5(a)-5(c))，而 本文提出的模型在扭转方向的响应比Beam188单元更接近Shell181单元的响 应(图5(d)),这说明用Beam188单元分析旋转双锥度悬臂预扭梁的扭转振动 时可能有错误。

3.旋转扭锥悬臂梁的动力学特性

在这一节中，分析讨论了安装角β₀和预扭角γ(L)对旋转双锥度悬臂 预扭倾斜梁模态特性和振动影响特性的影响。相关仿真参数见于表1：

表1：模态和瞬态动力学分析的相关参数

3.1预应力模态分析

考虑旋转造成的科氏力效应，离心刚化效应和旋转软化效应，分别 研究安装角和预扭角对一阶弯曲固有频率(Y_I)，一阶摆动固有频率(Z_I)和 一阶扭转固有频率(θ_I)的影响(图6和图7)。

3.1.1工况1：安装角对Y_I，Z_I和θ_I的影响

在图6中，当n≠0rev/min时，Y_I和θ_I是一直增大的，而Z_I随着安装角的 增大减小。安装角越大，Y_I和θ_I增大量越大(图6(a)和图6(c))，Z_I的减小量 越大(图6(b))。这是由于旋转软化效应在Y-和θ-方向是随着安装角减弱的 (见于公式(6)中

和公式(8)中 />

然而在Z-方向是随安装角增强的(见于公 式(7)中/>

这就导致了上述现象。

3.1.2工况2：预扭角对Y_I，Z_I和θ_I的影响

图7中一些典型的特征总结如下：

(1)随着预扭角的增大，Y_I和θ_I增大，而Z_I减小。这是由于随着预扭 角的增大梁单元的结构刚度矩阵在Y-和θ-方向上减小但是在Z-方向增大。

(2)预扭角越大，Y_I和θ_I的增量也越大(图7(a)和图7(c))，同时，Z_I的增量就越小(图7(b))。这主要是由于预扭角的增加导致Y_-和θ_-向旋转软 化效应的减弱(见于公式(6)中

和公式(8)中 />

而Z-向旋转软化效应的增强(见于公式(7) 中/>

这就导致了Y_I和θ_I的增大，Z_I的减小(图7)。

3.2基础激励作用下旋转双锥度悬臂预扭梁在X-,Y-,Z-和θ-向的响应 特性

考虑X-向基础激励(图3)，对旋转双锥度悬臂预扭梁在各方向的振动 响应特性进行了详细分析讨论。相关仿真参数见于表1。

3.2.1工况3：安装角对基础激励下各方向振动响应的影响

X，Y，Z和θ向的幅频响应如图(8)，可以得到以下的几个结论：

(1)当β₀＝0°和γ(L)＝0°时，基础激励作用下系统的共振仅发生在X向和Y 向(图8(a)-图8(d))，而Z向和θ向的幅频响应为0(图8(c)和图8(d)所 示)。这主要是由于科氏力效应直接导致了X向和Y向间的振动耦合(见于 公式(5)中

和公式(6)中/>

而X，Y和 Z向间不存在耦合项(见于公式(5)中

公式(6)中 />

和公式(7)中

此外，由基础激励导致的Z- 和θ-向的附加外载荷为0(见于公式(7)中/>

和公式 (8)中

这也直接导致了Z向和θ向振动响应的消失 (图8(c)和图(d))。

(2)当β₀＝±20°且γ(L)＝0°时，X向,Y向,Z向和θ向间均存在振动耦合 项，因而导致各方向间均出现了共振峰。这些耦合项均体现在公式(1)中。

(3)当β₀＝±45°且γ(L)＝0°，X向,Y向,Z向和θ向均出现了耦合振动， 但无共振峰的出现，这主要是由于在β₀＝45°/-45°，γ(L)＝0°时激振频率1×n/60未 达到系统的一阶弯曲频率(见图9)。

另外，图(10)则给出了不同安装角和转速下X向和θ向常值分量的 变化规律，需要说明的是，Y和Z向因不存在常值外载荷作用而导致这两 个方向的常值分量为0(见于公式(6)，(7)，(8)，(9)，(10))，此处不作讨 论。

在图10(a)中，X向常值分量不受β₀的影响且随转速的增大而增大。 这主要是由于本文中离心载荷(见于公式(5)中

与β₀无关 且随转速/>

的增加而呈增大的趋势。

从图10(b)中可以看出，当β₀＝±45°时，θ向常值分量达到最大，而当β₀＝0° 时，θ向常值分量达到最小(0rad)。此外，正/负β₀会导致正/负θ向常值分量 (图10(b)和图11)。这可以通公式(8)中外力矩项

进行解释。当β₀分别为β₀＝±45°和0°时， 会产生相应的最大扭矩和最小扭矩，进而分别导致最大和最小θ向常值分 量的产生。

3.2.2工况4：γ(L)对X-,Y-,Z-和θ-向振动响应的影响

与3.2.1节中β₀对各方向振动响应的影响不同的是，γ(L)(γ(L)≠0°)会直接 导致X向，Y向，Z向和θ向间振动响应的耦合，进而导致各方向共振峰 的出现(图12和图13)。相关耦合项可参考3.2.1中的内容。

在图14(a)中，X向常值分量并不受γ(L)的影响，这主要是由于X向 离心力(见于公式(5)中

与γ(L)无关。

从图14(b)可以看出，当γ(L)(20°≤|γ(L)|≤45°)互为相反数时，θ向常 值分量关于θ＝0rad对称。此外，正/负γ(L)会导致负/正θ向常值分量 的产生(图14(b)和图15)。γ(L)越大，θ向常值分量越大。造成这种现象 的主要原因为：(1)θ向扭矩为

(见于公式(8)，I_y>I_z且 β₀＝0°)，且当20°≤|γ(L)|≤45°，式中带双下划线项要大于其前一项；(2)负/ 正γ(L)会导致正/负扭矩的产生，进而导致正/负θ向常值分量的产生，见 图5(d)，图14(b)和图15。

4.结论

采用Hamilton原理结合有限单元法建立了旋转双锥度悬臂预扭梁的 弯-摆-轴-扭耦合动力学模型，并通过与文献和ANSYS结果的对比验证 了所推导模型的有效性。基于所推导的模型，分析讨论了安装角、预扭 角和基础激励对系统模态特性和振动响应特性的影响。主要结论如下：

(1)安装角/预扭角的增加会导致一阶弯曲/扭转固有频率的增加，同 时也会造成一阶摆动固有频率的减小。值得说明的是，当转速为0 rev/min时，系统的一阶弯曲，摆动和扭转固有频率仅受扭角的影响， 与安装角无关。

(2)在轴向基础激励作用下，会出现系统的弯-摆耦合甚至弯-扭-轴- 摆耦合振动。此外，正/负的安装角会导致正/负的扭转常值变形，而正/ 负预扭角会导致负/正的扭转常值变形。特别地，当预扭角为0°且安装角 为±45°时，扭转方向的常值变形达到最大，而当安装角为0°且预扭角互 为相反数时，扭转方向的常值分量大小相等，方向相反。

(3)在线弹性理论范围内，相对于Beam188模型，所推导模型和 Shell181模型在模拟旋转双锥度悬臂预扭梁的动力学特性时，其数值结 果更为相近，尤其是对其扭转力学特性的模拟。

附录A：弯、摆、轴、扭方向的单元形函数

(1)轴向形函数：

假设梁单元轴向位移插值模式为线性模式，即

则根据位 移边界条件有u(0)＝u_i，u(l_k)＝u_j(见图16所示)。

式(A.2)给出了位移插值系数

与节点位移u_i，u_j的关系式：

由式(A.1)可进一步得到单元内任一点x处的轴向位移u(x)与节点位移 u_i，u_j的关系式如下：

式中，轴向形函数N_u1(x)和N_u2(x)的表达式分别为

记

则N_u1(ξ)＝1-ξ,N_u2(ξ)＝ξ。

(2)扭转形函数

扭转方向的形函数采用与轴向位移相似的线性位移插值模式，即

则根据位移边界条件有θ(0)＝θ_i，θ(l_k)＝θ_j(见图17所示)。

位移插值系数

与节点位移θ_i，θ_j的关系式如下：

基于式(A.3)，可进一步得到单元内任一点x处的扭转位移θ(x)与节点位 移θ_i，θ_j的关系式如下：

式中，轴向形函数N_θ1(x)和N_θ2(x)的表达式分别为

记

则N_θ1(ξ)＝1-ξ,N_θ2(ξ)＝ξ。

(3)弯曲方向形函数

1)内力矩，剪力和几何变形间的关系

对于y^eo^ex^e面内的平面梁弯曲，取一微段梁进行研究(见图18)。在图 18(a)中，对微段梁左端取矩可得力矩平衡方程：

略去式(A.5)中的二阶高阶小量，可进一步化简为如下表达式：

图18(b)表明y^eo^ex^e面内的剪应变满足如下关系式：

/>

需要说明的是，本文中剪应变γ_xy采用常值假设，即γ_xy＝γ₀。

微元体中的几何变形关系满足如下方程：

式中，u为由于梁的弯曲造成的在水平方向的位移。

微元体中力矩M_z与截面转角

满足如下关系式(见图18(c))：

剪应力Q_Y和剪应变γ₀满足如下关系式：

Q_y＝κ_yGAγ₀ (A.10)

2)位移插值模式

弯曲方向平动位移v采用3次多项式进行插值，即

此外，由式(A.7)可进一步求得弯曲方向转角位移/>

基 于式(A.6),(A.9)和式(A.10)，并代入平动和转角位移插值模式，可得出如 下关系式：

记

根据式(A.7)，可进一步得出弯曲转角位移

根据弯曲平动方向的位移边界条件：v(0)＝v_i，v(l_k)＝v_j，

(见图19)。写成矩阵形式，可得如下表达式：

记

基于式(A.12)可得/>

和/>

的计算表达式如下：

基于式A.13(a)～A.13(d)，并结合弯曲方向的位移插值模式，可得如 下表达式：

/>

且

和/>

的表达式如下：

式中

(4)摆动方向形函数

1)内力矩，剪力和几何变形间的关系

对于z^eo^ex^e面内的平面梁弯曲，同样取一微段梁进行研究(见图20)。在 图20(a)中，对微段梁左端取矩可得力矩平衡方程：

略去式(A.16)中的二阶高阶小量，可进一步化简为如下表达式：

图20(b)表明y^eo^ex^e面内的剪应变满足如下关系式：

γ_xy＝φ+w′ (A.18)

需要说明的是，本文中剪应变γ_xz采用常值假设，即γ_xz＝γ₀。

微元体中的几何变形关系满足如下方程：

u＝zφ (A.19(a))

γ_xy＝φ+w′ (A.19(b))

式中，u为由于梁的摆动造成的在水平方向的位移。

微元体中力矩M_y与截面转角φ满足如下关系式(见图20(c))：

剪应力Q_z和剪应变γ₀满足如下关系式：

Q_z＝κ_zGAγ₀ (A.21)

2)位移插值模式

摆动方向平动位移w采用3次多项式进行插值，即

此外，由式(A.18)可进一步求得弯曲方向转角位移 />

基于式(A.17),(A.20)和式(A.21)以及，代入摆动平 动和转角位移插值模式，可得出如下关系式：

记

根据式(A.18)可进一步得出摆动转角位移

根据摆动平动方向的位移边界条件：w(0)＝w_i，w(l_k)＝w_j，φ(0)＝φ_i， φ(l_k)＝φ_j(见图21)。写成矩阵形式，可得如下表达式：

记

基于式(A.23)可得/>

和/>

的计算表达式如 下：

基于式A.24(a)～A.24(d)，并结合摆动方向的位移插值模式，可得如 下表达式：

且

N_φ1(x)，N_φ2(x)，N_φ3(x)和N_φ4(x)的表达式如下：

/>

式中

附录B：梁单元的单元矩阵

微段直梁单元的质量矩阵M^e、科氏力矩阵G^e、结构刚度矩阵

离 心刚化矩阵/>

旋转软化矩阵/>

角加速度引起的刚度矩阵/>

和力向量 F^e的表达式如下：

式中：

/>

式中：

式中：

/>

/>

式中：

式中：

/>

F^e＝[F₁ F₂ F₃ F₄ F₅ F₆ F₇ F₈ F₉ F₁₀ F₁₁ F₁₂]^T (B.6)

式中：

/>

附录C：矩形截面梁的抗扭惯性矩

任意矩形截面的扭转应力函数满足如下关系式：

式中，h为矩形截面短边边长，τ₁(y,z)为考虑厚宽比h/b影响的应 力修正项。

此外,τ(y,z)应满足泊松方程：

将式(C.1)代入式(C.2)可进一步得出如下关系式：

在矩形截面的边界处，应力函数τ(y,z)满足：

将公式(C.1)代入公式(C.2)，得到：

设τ₁(y,z)＝f(y)·f(z)并带入公式(C.3)可得：

其中，λ为任意常数。

由公式(C.6)可得如下表达式：

公式(C.7)的通解有如下形式：

根据薄膜比拟法，τ₁(y,z)应为坐标y和z的偶函数，故：

τ₁(y,z)＝Acos(λy)cosh(λz) (C.9)

将公式(C.9)代入公式(C.5)的第二个式子可得：

公式(C.10)的通解具有以下形式：

将公式(C.11)代入(C.9)，τ₁(y,z)可以展开成如下形式：

将公式(C.12)代入(C.5)，得到：

将公式(C.13)两边同乘以cos(λ_ly)，并在区间[-b/2,b/2]进行积分 可得：

将公式(C.14)代入(C.12)，可得：

结合公式(C.1)和(C.15)，可得任意矩形截面的抗扭截面惯性矩J的 计算表达式如下所示：

其中

是J_s的抗扭截面惯性矩修正因子，在本文中

即取l＝0。

以上结合具体实施例描述了本发明的技术原理，这些描述只是为了 解释本发明的原理，不能以任何方式解释为对本发明保护范围的限制。 基于此处解释，本领域的技术人员不需要付出创造性的劳动即可联想到 本发明的其它具体实施方式，这些方式都将落入本发明的保护范围之 内。

Claims (4)

1.一种机械结构在基础激励下的响应特性分析方法，其特征在于：包括以下步骤：

将机械结构的动力学分析退化成对受到基础激励的旋转悬臂梁模型的动态特性的分析；

通过汉密尔顿原理和有限元方法将基础激励下旋转悬臂梁模型离散成一系列的铁木辛柯梁段；

通过固有频率和基础激励响应的比较证明旋转悬臂梁模型的正确性。

2.根据权利要求1所述的机械结构在基础激励下的响应特性分析方法，其特征在于，

通过汉密尔顿原理和有限元方法将基础激励下旋转悬臂梁模型离散成一系列的铁木辛柯梁段具体包括以下步骤：

偏微分方程的建立：

基于第k个铁木辛柯梁单元，建立相应的弯扭轴摆耦合运动方程，其中，OXYZ，OX^rY^rZ^r，o^bx^by^bz^b和o^ex^ey^ez^e分别是整体坐标系，旋转坐标系，局部坐标系和单元局部坐标系，α为旋转角；

三个坐标系中的欧拉角θ，φ和

用来表示在单元局部坐标系任意梁截面处的角位移，β₀为叶片安装角，γ_k为第k个梁单元处的预扭角；

从而，第k个梁单元处的扭角β_k＝β₀+γ_k，R_k＝R_d+x_k+X_u ^r，其中R_k是o^e在旋转坐标系下的x坐标，R_d是叶盘半径，x_k是o^e在局部坐标系中的x坐标，而X_u ^r是X^r的基础激励；

梁截面在整体坐标系下任意点P的坐标可以表示为：

其中(x,y,z)是P点在ox₃y₃z₃坐标系的坐标；u，v和w分别是梁上任意截面形心'o'沿着x^e-，y^e-和z^e-轴正方向的线位移；A₁，A₂，A₃，A₄和A₅是从ox₃y₃z₃到OXYZ的转换矩阵，它们可以表示为：

在公式(2)中，基于小变形假设采用一阶近似，即：

sin(φ)≈φ，cos(θ)≈1和sin(θ)≈θ；

根据公式(1)和公式(2)，第k个铁木辛柯梁单元的动能T可以表示为：

其中，ρ，l_k和A分别为密度，第k个梁单元的单元长度和截面面积，符号上点表示对时间的导数；

第k个梁单元的势能U可以表示为：

其中，E和G分别是弹性模量和剪切模量；I_x，I_y和I_z分别是关于轴ox₃，oy₃和oz₃的惯性矩；J是扭转惯性矩；κ_y和κ_z分别是关于y-和z-方向的剪切系数；

是离心力；

是基础激励在x-方向产生的力；

是一个常值；

另外，公式(4)中的求导符号表示对x的偏导数；在公式(4)中，下划线的一项表示离心刚化效应对扭转变形的影响，而双下划线的一项表示预扭角对势能的影响；

视u，v，w，θ，φ和

为独立变量，并应用汉密尔顿能量方程

从而获得六个偏微分方程；

对于x^e-轴向的位移u，有：

y^e-弯曲方向的位移v为：

z^e-摆动方向的位移w为：

x^e方向的角位移θ为：

y^e方向的角位移φ为：

z^e方向的角位移

为：

3.根据权利要求2所述的机械结构在基础激励下的响应特性分析方法，其特征在于，

铁木辛柯梁单元的单元矩阵

为了求解公式(5)-(10)，u-，v-，w-，θ-，φ-和

方向的假设位移场函数及相关表达式，将相关表达式代入公式(5)-(10)，可以得到如下的二阶常微分方程：

其中，M^e和

分别表示单元质量矩阵和结构刚度矩阵；Ge，

和

分别是当角加速度

时的科氏力矩阵，离心刚化矩阵和旋转软化矩阵；

是与角加速度

有关的刚度矩阵；δ^e和F^e分别表示节点位移和外力向量；

单元矩阵的组集

所有的单元矩阵(M^e，G^e，

δ^e和F^e)都是从单元局部坐标系o^ex^ey^ez^e中得到的；然而，在组集这些单元矩阵以建立旋转双锥度悬臂预扭梁时，必须统一相应的单元坐标系；

所有的单元矩阵都转化到旋转坐标系OX^rY^rZr下进行组集，转换矩阵公式如下：

其中，T是从OX^rY^rZ^r到o^ex^ey^ez^e坐标系的转换矩阵；δ^e和

分别是节点在o^ex^ey^ez^e和OX^rY^rZ^r坐标系中的位移向量；

基于公式(11)和公式(12)，旋转坐标系下的单元矩阵如

和

可以这样计算：

那么此悬臂梁梁总的矩阵χ(χ＝M,G,K_e,K_c,K_s,K_acc,F)可以通过

得到，其中N是梁单元的总数。

4.根据权利要求2所述的机械结构在基础激励下的响应特性分析方法，其特征在于，

扭形直梁的扭转刚度

计算矩形直梁扭转惯量(J_s)的计算公式：

另外，预扭角γ(L)的存在会产生额外的扭转惯量J_a：

其中，γ(L)是在L位置处梁的扭角；

从公式(14)和公式(15)可得，J＝J_s+J_a。

Images