第 28 卷 第 3 期 振 动 工 程 学 报 Vol.28 No.3 2015 年 6 月 Journal of Vibration Enginering Jun.2015
一种基于应变插值大变形梁单元的刚-柔耦合 动力学分析*
张 志 刚 1 ， 齐 朝 晖 1 ， 吴 志 刚 1 ，2 (1.大连理工大学工业装备结构分析国家重点实验室，辽宁 大连116024;
2.大连理工大学航空航天学院，辽宁 大连116024)
摘要:柔性多体系统动力学中的“动力刚化”现象起因于变形间的耦合。一次近似模型成功地解决了小变形情况下 的刚柔耦合建模问题，但在大变形情况下则需要考虑更多的耦合效应。选取表征梁弯曲应变的曲率和轴向应变作 为单元参数进行离散;在大变形大转动基础上得到了单元两端节点运动学参数的递推关系，构造出了能够自动计 及“动力刚化项”且适用于大变形刚柔耦合动力学分析的平面梁单元。最后采用所提应变插值单元求解了包含大 变形和刚柔耦合动力学柔性梁的数值算例，验证了算法的正确性和有效性。
关键词:刚柔耦合;梁单元;应变插值;动力刚化
中 图 分 类 号 : O 3 1 3 .7 文 献 标 志 码 : A 文 章 编 号 : 1 0 0 4 -4 5 2 3 (2 0 1 5 )0 3 -0 3 3 7 -0 8
DOI:10.16385/j.cnki.isn.1004-4523.2015.03.001
引 言
梁的刚柔耦合动力学分析是近年来力学、机械 及 航 空 航 天 领 域 研 究 的 一 个 热 点 和 难 点 [1-12]。 柔 性 梁的运动包含大范围刚体运动与变形运动的耦合， 且描述梁变形所需的形心位移和截面角之间也存在 耦合关系。这些耦合的作用机理非常复杂，给做大 范围运动梁的刚柔耦合动力学的建模和求解带来了 困难。
传统的柔性多体系统动力学混合坐标法直接套 用 了 结 构 力 学 小 变 形 假 设 ，Kane[1]在 用 其 对 旋 转 柔 性梁进行求解时发现高速转动时数值结果发散，表 现为随着转速增加梁变得更加柔软，这与事实相违 背，并首次提出了“动力刚化”的概念。这一现象引 发了国内外广大学者的关注，此后大量研究围绕刚 体运动与弹性运动的耦合及“动力刚化”问题展开: 洪 嘉 振 [2-6]指 出 ，传 统 混 合 坐 标 法 在 小 变 形 基 础 上 忽 略了纵向位移对轴向位移的影响是一种零次近似模 型，低速转动忽略掉变形耦合对仿真结果影响不大， 但在高速转动情况下这种近似将会产生完全错误的 结果;其提出小变形情况下在轴向位移表达式中补 充纵向位移二阶项的所谓一次近似模型，用该模型
*收 稿 日 期 :2013-12-26;修 订 日 期 :2014-06-12
基 金 项 目 :国 家 自 然 科 学 基 金 资 助 项 目 (11372057)
进行了中心刚体旋转柔性梁及末端带集中质量的旋 转柔性梁动力学特性，并进一步通过实验验证了一 次 近 似 模 型 能 够 捕 捉 “动 力 刚 化 项 ”。 章 定 国 [7 ] 在 一 次近似模型基础上对平面运动柔性梁进行了进一步 研究，指出大范围运动与弹性运动的耦合不仅会产 生 “动 力 刚 化 ”现 象 ，还 会 产 生 “动 力 柔 化 ”现 象 ，并 对 旋转内接悬臂梁的“动力柔化”效应进行了分析。和 兴 锁 等 [8]在 Shi P 等 [9]研 究 基 础 上 提 出 一 次 近 似 模 型未能充分考虑几何非线性变形对横向变形的影 响，因此在纵向与横向位移中增加了新的变形二次 耦合项。
上述研究大都局限于单元本身为小变形情况， 针对包含大变形的柔性梁刚柔耦合问题，刘锦阳[10] 采用绝对坐标法对旋转柔性梁进行了研究，并与一 次近似模型进行了对比，发现随着梁刚度降低一次 近似模型结果发散，而绝对坐标梁单元能够得到收 敛的结果。陈思佳[11]在一次近似模型基础上采用 模态坐标对梁进行离散，通过保留了动力学中方程 中的变形耦合量相关的高阶项，对考虑大变形的平 面梁刚柔耦合动力学分析做了初步研究，但是通过 模态坐标离散的求解精度和效率取决于模态集的选 取方式。
现有研究大都选取节点位移作为有限元离散参
338 振 动 工 程 学 报 第 28 卷
量，在小变形情况下根据伯努利梁变形耦合关系添 加位移高阶耦合项来构造了一次近似模型或者高次 近似模型。而直接构造满足大变形大转动的变形耦 合位移梁单元仍然十分困难和复杂。本文在对曲率 插值梁单元研究[12]的基础上，选取梁的弯曲应变 (曲率)和轴向应变作为单元参数进行离散，从而得 到了自动计及位移耦合关系的梁单元。采用本单元 对平面梁大范围运动刚柔耦合问题进行了研究，所 得 结 果 与 相 关 工 作 [1 3 -1 6 ] 吻 合 ， 从 而 验 证 了 所 构 造 单 元的正确性和有效性。
1 大范围运动梁单元虚功率方程 1 .1 梁 单 元 运 动 学 描 述
其中(·)′=d(·)/ds，θ为截面相对于单元坐标 系的转角。
单元左端截面相对于整体坐标系的转角φ1 表 征了单元大范围刚体转动，截面相对于单元坐标系 的转角θ代表了单元的弹性变形，因此单元域内截 面相对于整体坐标系的转角为
平面梁单元如图1所示，整体惯性坐标系为 δpe = (EAεδε+EIθ′δθ′)ds (6)∫
{g 1 g 2 } ， 在 梁 左 端 截 面 建 立 固 结 在 截 面 上 的 单 元 坐 标系 {e1e2}，第一根轴e1 指向梁长方向，初始时刻 与整体坐标系第一根轴夹角φ0 。在每个梁截面建 立 随 截 面 转 动 的 截 面 坐 标 系 {t 1t 2 } ， 变 形 前 截 面 坐 标系与单元坐标系方向保持一致。
图1 大范围运动平面梁单元
Fig.1A plane beam element undergoes large range of
motion 如图1所示，弧长s处截面形心矢径为
由上式可知，梁单元应变虚功率只与曲率和轴向应 变有关，而不受单元位移和转动大小的影响，因此选 取梁曲率和轴向应变进行离散可以方便地得到单元 应变虚功率和简洁的单元刚度矩阵。
相比于弯曲变形，梁的轴向变形一般为小量，因 此可进一步认为单元域内轴向应变ε为常量，曲率 θ′按线性变化，可由单元两端曲率θ′1，θ′2 做线性插
φ = φ 1 + θ 1 .2 梁 单 元 应 变 场 的 离 散
(4 )
对于不超过弹性范围的平面梁，其本构关系可 用截面坐标系中力的分量描述为
F = EAε;M = EIzθ′ (5) 式 中F 为 截 面 轴 向 力 ;M 为 弯 矩 ;EA ，EIz 分 别 为
截面拉压刚度和弯曲刚度。梁的应变虚功率 l
0
值得到 其中
() 珡 θ′s =Nγ
()
T ;珡
γ=[εθ′1 θ′2]N=[0 1-ξξ]8
( ) 珡 ()
其 中 ξ = s/l ，梁 曲 率 速 率 虚 变 分 由 式 (7 )可 得
δθ′ = Nδγ 9 将 式 (7 ) ， (9 ) 及 轴 向 应 变 ε 是 常 量 的 假 设 代 入 单 元
7
(1 ) 式中 r1 为单元坐标系原点相对整体坐标系原点 矢径，ρ为弧长s处截面形心相对于单元坐标系原 点矢径。取出原始弧长为ds的一段梁微元，变形后
应 变 虚 功 率 式 (6 )
δpe = (δγ)TKγ (10)
式中K 为单元刚度阵
EAl 0 0 熿燄
EIl/3EIl/6 (11) //
r = r 1 + ρ
K = 0 燀 0
几何关系
式 中 ε 为 轴 向 应 变 ， 即 : ε = d s珋 / d s – 1 ; d s珋 为 变 形
后梁微元长度。单元坐标系下截面形心坐标的弧长 导数为
单元域内截面的转角可由曲率插值函数式(7) 积分得到
dρ/ds= (1+ε)t1 (2)
x′=e1·′=(1+ε)cosθ
ρ (3) [ξ-ξξ]
EIl 6EIl 3燅 1 .3 单 元 内 运 动 学 递 推
其中
(13) w′=e2 ρ′= 1+εsinθ 转角θ为截面相对于单元转角 它完全由单元的弯
N= 0 (2 2)l/2 2l/2 {·() ，
s
θs= Nγds=Nγ () ∫珡
0
( )
12
烄
0l
[]
即 :ε = T εγ ;矩 阵 函 数 g (A ，B ) = A T B + B T A ;与 截
第3期 张志刚，等:一种基于应变插值大变形梁单元的刚-柔耦合动力学分析 339
曲 变 形 确 定 ，将 其 代 入 式 (12)可 得 截 面 相 对 整 体 坐 标系转角为
φ = φ 1 + N γ (1 4 ) 由 梁 截 面 形 心 坐 标 弧 长 导 数 式 (3 ) ， 对 弧 长 s 积 分 并
注意到轴向应变ε 为常量，可以得到截面形心坐标 s
( ∫) 烄x= 1+ε cosθds
0
烅s(15) ()
^^
H=熿N1 N2 N2 0 0燄 (21)
l2 -1，N2=2 3-2。 ξ(ξ) ξ(ξ)
cos 1 + 0 – sin 1 + 0
A= (φφ) (φφ)(23)
ss
εds- (w x )ds
10 到形心加速度
烆
18 其中sθ =sinθ，cθ =cosθ，λ=1+ε，单元坐标系下 左端面形心坐标及截面转角为0，即:x0 =0，w0 = 0，θ0 =0;右端面坐标和转角可由式(14)和(15)得 到。因此可由赫米特插值重新构造单元域内形心坐 标 x ， w ， 用 ρ珔 代 表 单 元 坐 标 系 下 截 面 形 心 坐 标 ， 则
有 式中
1∫
/2
烅s φρρφρφρ
0l ()
ρ珔=(x w)T =HΨl (ll)
(19) Ψl=λxlλcθwlλsθ (20)
φ
T
̈ ̈ ̈珘̈
() ()
Dl= (θl
燀0 0 0 N2 ^N2燅 ^2^
其中赫米特插值多项式N l1 ，N
1 =( – ) 2 =
弧长s处截面形心坐标为
r(s) r A珔 (22)
[ ]
sθTε +λcθNl 燀ll燅
熿 01×3 燄 T(，)
ξξ
烆
0
w=1+ε sinθds
从这两个式子可以看出截面形心坐标与转角和轴向 应变高度耦合。需要说明的是传统混合坐标零次近 似模型及一次近似模型都假设转角为小量，即:
θ 0 ，从 而 将 三 角 函 数 近 似 为 :sinθ ≈ θ ，cosθ ≈ 1 – θ 2 / 2 ， 并 由 式 (3 ) 得
////() w x = (w s ) (x s )= tanθ ≈ θ 16
[]
将 其 带 入 式 (15)得 烄x = s +∫
φρρ
其中I珘= 0 -1 。对上式再求一次时间导数，得
∫
=1+ρ 其中单元坐标系与整体坐标系转换矩阵
sin 1 + 0 cos 1 + 0 (φφ) (φφ)
式中 φ0 为单元初始方位角。
对 式 (22)求 时 间 导 数 得 到 截 面 形 心 速 度
r(s)=r1+ 1AI珘珔+A (24)
02 w =∫ w/x ds
0
2
珔 珓 ( )
(17)
̈()̈̈ 珓珔 rs=r1+1AI+A-1A+21AI 25
0
()
根 据 单 元 坐 标 系 下 截 面 形 心 坐 标 式 (19)，可 得 ̈
烆 上式推导过程中忽略掉了转角θ的二次以上项
及θ和轴向应变ε的耦合项，式中加下划线部分为 一次近似模型在零次模型基础上增加的横向位移对 轴线变形的二次耦合项。由此可见混合坐标一次近 似模型不适用于包含大转动的大变形情况。
本文不作小转角假设，基于此推导的单元完整 地保留了变形间的耦合关系，不仅适用于小转角情 况 ，同 样 也 能 满 足 大 转 角 情 况 。 由 式 (12)可 知 截 面 转 角 θ 为 弧 坐 标 s 的 二 次 式 ， 因 此 式 (1 5 ) 两 个 积 分 不存在有理形式解，本文由高斯积分得到数值解。
直 接 由 式 (15)就 可 得 到 截 面 形 心 坐 标 场 ，但 是 在推导质量阵和广义力时将会出现二重数值积分， 从而会引起公式推导和数值求解上的不便。再次考 察单元坐标系下截面形心坐标x，w 在端部满足
x(0)=x0，x′(0)=λcθ ;x(l)=xl =x′(l)=λcθ
ρ=HΨl=HDlγ;ρ=HDlγ+HDlγ (26) 其中
熿 Tε 燄 ClTε -λSnl
cT sN
θ ε-λθ l ;
Dl =
SlTε +λCnl
ll
-γ(gSnlTε +λCnl) TT
-γ s (N ，T )+λc N N
g l ε θl l l) (27)
T(，) γ(gCnlTε -λSnl)
T(，) T
燀γ (cθg Nl Tε -λsθNlNl)燅
ll
式中 Tε = 1 0 0 为单元轴线应变定位矩阵，
烅 () ， () ;() ， () w0=w0w′0=λsθ wl=wlw′l=λsθ
面转角θ相关积分系数 lll
∫ ; ∫ ; ∫ T烄Sl = sθdsSnl = sθNdsSnl = sθN Nds
000
烅lll
Cl = cθdsCnl = cθNdsCnl = cθN Nds
(28)
将 式 (26)代 入 得 到 截 面 形 心 速 度
()珘 () rs =r1 +AIHΨl 1 +AHDlγ 29
对速度取变分
()珘 () δrs =δr1+AIHΨlδ1+AHDlδγ 30
烆000 ∫;∫;∫T
将 式 (26)代 入 式 (25)得 到 截 面 形 心 加 速 度
rs =r1+ 1AIHΨl+AHDlγ+Aar 31 φ
φ
340 振 动 工 程 学 报 第 28 卷
其中加速度余量
珘() ar =HDlγ- 21HΨl+21IHDlγ 32
质量阵和广义力中常系数矩阵
1 1
烄NH = Hdξ;NHH = HTHdξ;
φφ
∫
∫
对 截 面 转 角 式 (14)求 导 可 得 截 面 转 动 角 速 度 及 虚 角 速度
角加速度
0
0
∫l T ̈
() () δpa = [δr r+δωIzω]ds 39
ρ 其中ρ代表梁的线密度，Iz 代表梁截面惯性矩。
fqg= I mg (50)
将 式 (30)~ (34)代 入 上 式
AINHΨl
(AN D )Tm (51)
1T q1 qγ q γg ̈ ̈ f=Hlg
δp a = (δq ) (M q + M γ + F ) +
另外，还可以将结构阻尼当做外力处理，本文采 用 文 献 [13]阻 尼 模 型 ，单 元 结 构 阻 尼 力 虚 功 率 为
T γq1 γ γ ̈ ̈
(δγ ) (M q + M γ + F ) 质量矩阵和广义力为
Mr Mrφ M q q = [
l
Mqγ = [ φγ
Fq = () φ
其中
烄
φr φ M M
0
式中为阻尼系数。对比单元应变虚功率式 η
Mrγ
;Mγq = (Mqγ)T (41) ]
(10)，可将机构阻尼力虚功率表示为
M Fr
(42)
δpd = (δγ)Tfγd (53) 其中广义结构阻尼力列阵为
F
fγd =ηKγ (54) 2 平面梁刚柔耦合动力学方程
2 .1 单 元 间 运 动 学 递 推
r ;φ T
M = mI M = mΨ l N H H Ψ l + mIz
(40 )
; δ p d = η E A ε δ ε + η E I zθ ′ δ θ ′ d s 5 2
()() ]∫
11
∫;∫T;() NN = NdξNNN = NNdξ 45
烅
烆
1 .5 梁 单 元 外 力 虚 功 率
梁单元上作用的外力虚功率为
; () ω = 1 + Nγ δω = δ 1 + Nδγ 33
0
1
0
φφ
̈ ̈() ω = φ1 + Nγ 34
为了方便推导，本文将整体坐标系坐标下截面形心 坐标r和转角φ 组成列阵为q ，两端面运动学参数
r1 r2
1烄烌2烄烌 2
∫珘
TT NHIH = HIHd
0
ξ
其中
其中
ll
IAIHlΨl 1 .4 梁 单 元 惯 性 力 虚 功 率
(珘) 珦 m1 +m2 +Ψl AIHl f2 +ΨlΦp
B =
燀燅
q = ，q = (35) φφ
∫l
( )T ( )T ()
烆 1烎 烆 2烎
由 式 (14)，(29)，(31)，(33)，(34)，单 元 两 端 截 面 运
δpf =∑[δri fi+δωimi]+ δr psds i=1 0
动学递推关系可以用矩阵表示为
; ; ̈̈̈
q2 =q1+σq2 =Bq1+Qγq2 =Bq1+Qγ+ν
(46) 式中 p(s)为分布载荷;fi，mi 分别为单元节点处
所 受 集 中 力 和 力 矩 。 将 式 (30)，(33)代 入 上 式
AHlΨl Aalr
σ=，ν=烄烌(37)
∫l
f1+f2+ pds 燄
N lγ 珘
烆 0 烎
fq =
γ T T()T T ()
(36)
δpf =(δq1)Tfq+(δγ)Tfγ (47)
[]
熿
AHlDl [0 1 ] [Nl ]
根据伯努利梁理论，平面梁可以看作由无数做 平面运动的刚性截面组成，其惯性力虚功率为
0
珦 ∫(珘 ) ; ∫( ) 。 式中 Φp = AIH TpdsΦp = AH Tpds
;Q =
(38)
TTT
(48) f=m2Nl+DlAHl f2+DlΦp 49
0
对于考虑重力的情况，可将重力看做一种特殊 的分布力，其对应广义力为
[珘 ] ()
0
T
0
; 烅rφ珘;rγ ;
烆
γ T
M = m D l N H H D l + m I z N N N (4 3 )
M = mAINHΨl M = mANHDl
φγ T
M = mΨlNHIHDl +mIzNN .
r2
( 珘)
由于选取单元应变进行离散，截面形心坐标及
转动角度由几何方程积分得到，并给出了单元两端
运动学递推关系式( )。因此单元两端运动学参数
F = m A N H D lγ – N H Ψ l 1 + 2 I N H D lγ 1
烄
φ φ
φ T( 2 ) F = m Ψ l N H I H D lγ – N H I H Ψ l 1 + 2 N H H D lγ 1
烅
烆 F γ = m D lT N H H D lγ – N H H Ψ lφ 21 + 2 N TH I H D lγ φ 1
φφ
( ) (44)
36
并不独立，如果知道了单元一端运动学参数，另一端
第3期 张志刚，等:一种基于应变插值大变形梁单元的刚-柔耦合动力学分析 341
参数就可以递推得到。对于由n个单元组成的梁系
统，可以选定一个节点为运动学递推起始点，其他节
将 梁 单 元 应 变 虚 功 率 式 (10)、惯 性 力 虚 功 率 式 (40)、
外 力 虚 功 率 式 (47)及 单 元 间 运 动 学 递 推 关 系 式 (55)
代入上式得到
点的参数便可以通过递推关系得到。本文把选定的 1
起始节点编号规定为1，其运动学参数记为 。为 q0
1T 001 0Γ 0 ()̈̈
了便于统一进行运动学递推，这里约定单元坐标系 的选取满足下面规则:
(1)单元坐标系原点选在靠近起始点一端节点， 第一根轴方向背离起始节点指向第二个节点;
(2)沿 背 离 起 始 节 点 方 向 单 元 编 号 递 增 。 单元连接拓扑关系通过定义邻接单元数组L ∑
描 述 :L 为 1 × n 行 向 量 ， 其 中 L (i ) 代 表 与 第 i 号 单 元相邻并靠近起始节点1的那一侧单元编号;规定 与起始节点相连单元的邻接单元为0。
如图2所示，下面单元标号满足规则标号，箭头 标出了单元坐标系原点所在节点和第一根坐标轴方 向。
图2 梁单元规则标号 Fig.2 Regularly labeling of beam elements
按照上述约定，图示由三个单元组成的梁的邻 接 单 元 数 组 为 :L = [0 1 2 ]。
利用邻接单元数组L，根据单元内运动学递推
关 系 式 (3 6 ) ， 可 将 i 号 单 元 坐 标 系 所 在 节 点 运 动 参
量与起始节点关系表示为
; ; ̈̈̈
q i1 = q 10 + κ i q i1 = Ω iq 10 + Λ iΓ q i1 = Ω iq 10 + Λ iΓ + τ i
δq0 Mq0+MΓ+F+ ()
T Γ01 ΓΓ Γ () ̈ ̈ ()
δΓ (M q0 + M Γ + F )= 0 59 其中质量矩阵
M
Ωi Mi Ωi ∑
n
00 Tq
，
=
M = Ωi MiΛi+MiTi M = M
i=1 n
qγ ; Γ0 ()
( 0Γ )T ， ΓΓ [Tq qγ T(γq
0Γ
T q q
i=1 n
M = ∑ Λi MiΛi +MiTi +Ti ()
MiΛi +
Mi Ti 广义力列阵
i=1
γ )]
(60)
F=∑Λi Miτi+Fi-fi + [()
n
0 Tqqq，
F= Ωi Miτi+Fi-fi
∑
() n()
i=1
Γ Tqqq 61
i=1
Tγqγ γ
T i ( M i τ i + F i + K iγ i – f i ) ]
由 梁 系 统 虚 功 率 方 程 式 (59)中 虚 变 分 独 立 ，其 系 数 为零可得平面梁动力学方程
M00 M0Γ̈q10 Fi0
烄烌烄烌 ()
其中
(55)
3 数值算例
在这一节中采用本文单元编写了动力学分析程 序，分别选取了悬臂梁纯弯曲、大范围运动已知的旋 转柔性梁及包含大变形的自由下落柔性单摆3个算 例 进 行 了 求 解 ，并 将 所 得 结 果 与 相 关 文 献 [10-11，14-15] 进行了分析比较。
3 .1 悬 臂 梁 受 梁 端 集 中 弯 矩 作 用
悬 臂 梁 端 受 集 中 弯 矩 M 作 用 ，梁 长 度 L = 10 m ，梁 端 线 位 移 u ，v 如 图 3 所 示 。
由于存在结构阻尼，梁最终会处于静平衡状态，
这里采用本文算法将整个梁划分为5个单元，为了
减少悬臂梁达到平衡状态所用时间，这里选取阻尼
比 =0.1进行仿真。在不同弯矩作用下梁最终的 η
变形状态如图4所示。
κi =κL(i)+σL(i);Λi =BL(i)ΛL(i)+QL(i)TL(i) {;
Ωi =BL(i)ΩL(i)τi =BL(i)τL(i)+νL(i)
单元应变参数列阵为: T T … T T ; Γ=γ1 γ2 γn
() Ti 为i号单元应变参量定位列阵，即:γi =TiΓ ;上
式中相关量初值为
κ0 =0;σ0 =0;τ0 =0;ν0 =0; {;;;;
(57) Ω0 =IB0 =IΛ0 =0Q0 =0T0 =0
2 .2 梁 系 统 刚 柔 耦 合 动 力 学 方 程
弹性体虚功率方程表述为:弹性体所受外力虚 功率等于惯性力虚功率和应变虚功率之和。划分为
n个单元柔性梁的虚功率方程为 nnn
δpia + δpie = δpif ∑∑∑ i=1 i=1 i=1
(58)
(56)
[ + =0 62
] ̈ Γ 烆烎烆烎
Γ0 ΓΓ
M M Γ Fi
这是一个包含大范围刚体运动参数 ̈q10 和弹性 运动参数 ̈Γ 的刚柔耦合动力学方程，可以采用适用 于刚性微分方程的数值方法求解。
342
振 动 工 程 学 报 第 28 卷 相 同 :梁 长 L = 10 m 截 面 面 积 A = 1 m 2 ，截 面 惯 性
图3 悬臂梁受集中弯矩作用 Fig.3 Cantilever beam undergoes bending load
矩I =5×10-4 m4，弹性模量E=2.8×107 Pa，密度 z
ρ = 1 .2 k g / m 3 。
图6 旋转柔性梁 Fig.6 A flexible beam on a rotating base
柔 性 梁 绕 转 动 中 心 的 转 角 φ (t ) 规 律 为 6 t2 15 2 2πt
cos – 10 ≤ t ≤ 15 s [(2)(15 )]
烄 +
φ(t)=152 π 烅
烆6t – 45t > 15 s 这里忽略梁结构阻尼，将梁划分为5个单元，仿真时
(63) 间为30s，得到梁自由端位移响应曲线如图7所示。
图4 不同弯矩作用下悬臂梁的大变形 Fig.4 Large deformation of the cantilever under bending
moments
下面给出悬臂梁在弯矩 M =2πEIz/L 作用下
梁末端位移随时间变化曲线如图5所示。
图5 悬臂梁末端受集中弯矩的位移
Fig.5The tip displacements of the cantilever with an
end bending load
从位移时间曲线图可以看出，在较大结构阻尼 作用下悬臂梁经过两个周期就达到了最终的平衡状 态，且仅用5个单元就得到了与解析解高度吻合的 结果，由此表明本文单元在处理包含大变形问题的 优越性。
3 .2 旋 转 柔 性 梁 大 范 围 运 动 已 知
大范围运动已知的旋转柔性梁已被很多学 者 [14-15]用 来 考 察 所 提 方 法 补 充 “动 力 刚 化 项 ”的 能 力。如图6所示，柔性梁物理参数选取与文献[14]
图7 自由端位移曲线 Fig.7 Displacement curves at the fre end
从图7可以看出:梁端部位移随着转动速度增 加而增大，当转速增加到一定数值维持不变，轴向位 移 逐 渐 回 落 并 稳 定 在 5 .1 4 × 1 0 – 4 m ， 横 向 位 移 在 0 附近小幅度波动。采用本文单元的仿真结果与文 献 [14-15]完 全 相 符 ，可 以 验 证 算 法 对 动 力 刚 化 问 题 的
第3期 张志刚，等:一种基于应变插值大变形梁单元的刚-柔耦合动力学分析 343
适用性。
3 .3 自 由 下 落 柔 性 单 摆
最后一个算例研究柔性单摆在重力作用下自由 下落问题，单摆初始时刻水平放置如图8所示。
图8 柔性单摆 Fig.8 The flexible pendulum
型仿真结果发散，表明其不适用与包含大变形的梁 的刚柔耦合动力学分析。而本文采用应变插值所得 结 果 与 文 献 [10 -11 ]相 吻 合 ，具 有 较 好 的 计 算 精 度 。
从两种方法所用 CPU 时间对比可以看出，本 文单元求解效率优于文献[10]绝对节点坐标单元， 也具有很好的计算效率。
4 结 论
本文通过直接选取单元应变进行插值构造了一 种不受单元变形大小限制，适用于刚柔耦合动力学 分析的平面梁单元。梁变形的位移及转角之间的耦 合关系通过几何方程积分得到，并由此给出了节点 运动学递推关系。本文所提梁单元有如下特点:
(1)选取单元应变进行离散不受单元刚体运动 的影响，得到了简洁的单元刚度矩阵及节点力;
(2)根据应变场由几何方程积分得到梁截面形 心坐标和转角耦合关系，自动计及了“动力刚化项”;
(3)单 元 不 受 变 形 大 小 的 限 制 ，适 用 于 包 含 大 变 形梁的刚柔耦合动力学分析。
参考文献:
[1] KaneT R，Ryan R，Banerje A K.Dynamics of a can- tilever beam atached to a moving base[J].Journal of Guidance，Control，and Dynamics，1987，10 (2):
139—151. [2] 杨辉，洪嘉振，余征跃.动力刚化问题的实验研究
[J].力 学 学 报 ，2004，36(1):118—124. YANG Hui，HONG Jiazhen，YU Zhengyue.Experi- mental investigation on dynamic stifening phenomenon [J].Acta Mechanica Sinica，2004，36(1):119—124.
[3] 杨辉，洪嘉振，余征跃.刚柔耦合建模理论的实验验 证 [J].力 学 学 报 ，2003，35(2):253—256. YANG Hui，HONG Jia-zhen，YU Zheng-yue.Experi- ment validation of modeling theory for rigid-flexible
coupling systems[J].Acta Mechanica Sinica，2003，35 (2):253—256.
[4] Cai G P，Hong J Z，Yang S X.Dynamic analysis of a flexible hub-beam system with tip mas[J].Mechanics
Research Communications，2005，32(2):173—190. [5] Cai G P，Lim C W.Dynamics studies of a flexible hub- beam system with significant damping efect[J].Jour-
nal of Sound and Vibration，2008，318(1):1—17. [6] 刘锦阳，洪嘉振.柔性梁的刚-柔耦合动力学特性研究
这个算例的大范围刚体运动是未知的，是个典 型的大范围运动与变形运动耦合的问题。针对这个 算 例 刘 锦 阳 [1 0 ] 和 陈 思 佳 [1 1 ] 分 别 采 用 不 同 建 模 方 法 进行仿真求解，并与混合坐标法一次近似模型进行 了对比。通过选取较低的材料弹性模量使得弹性变 形不再是小量，以此来验证方法对大变形问题的适 用 性 。 柔 性 摆 的 长 度 L = 1 .8 m ，E = 6 .895 × 109
P a ，A = 2 .5 c m 2 ，I z = 0 .1 3 c m 4 ， 密 度 ρ = 27 6 6 .6 7 kg/m3 。
这个算例中将单摆划分为5个单元，以本文方 法和绝对节点坐标方法[10]在 Matlab环境中选用刚 性 微 分 方 程 组 求 解 器 ode23tb，数 值 积 分 参 数 设 定 为:绝对误差限 AbsT=0.1，相对误差限 RelT=
0 .0 1 ， 在 同 一 台 微 机 上 进 行 仿 真 2 .5 s ， 所 用 C P U 时 间如表1所示。
表1 绝对坐标方法与本文方法 CPU 时间对比 Tab.1 The CPU time comparison betwen absolute cordinate
method and the proposed method
方 法
C P U 时 间 两种方法所得单摆末端位移时间曲线如图9所示。
根 据 文 献 [10 -11 ]结 果 ，混 合 坐 标 法 一 次 近 似 模
图9 柔性单摆的末端位移 Fig.9 The tip deformation of the flexible pendulum
文 献 [10]绝 对 坐 标 法 2 2 5 .0 5 s
本 文 方 法
6 .2 9 s
344 振 动 工 程 学 报 第 28 卷
[J].振 动 工 程 学 报 ，2002，15(2):194—198. Liu Jinyang， Hong Jiazhen.Study on rigid-flexible coupling dynamic behaviour of flexible beam[J].Jour-
nal of Vibration Enginering，2002，15(2):194—198. [7] 方建士，章定国.旋转内接悬臂梁的刚柔耦合动力学
特 性 分 析 [J].物 理 学 报 ，2013，62(4):305—311. Fang Jianshi，Zhang Dingguo.Analyses of rigid-flexi- ble coupling dynamic properties of a rotating internal cantilever beam[J].Acta Physica Sinica，2013，62(4): 305—311.
[8] 邓峰岩，和兴锁，杨永锋，等.计及非线性变形的刚- 柔耦合 动 力 学 建 模 [J].机 械 强 度，2006，28(6):
800—804. Deng Fengyan，He Xingsuo，Yang Yongfeng，et al. Dynamics modeling for a rigid-flexible coupling system with onlinear deformation field[J].Journal of Mechan-
ical Strength，2006，28(6):800—804 [9] Shi P，McPhe J，Heppler G R.A deformation field
for Euler-Bernouli beams with applications to flexible multibody dynamics[J].Multibody System Dynamics，
2001，5(1):79—104. [10]李彬，刘锦阳.大变形柔性梁系统的绝对坐标方法
[J].上 海 交 通 大 学 学 报 ，2005，39(5):827—831. Li Bin，Liu Jinyang.Application of absolute nodal co- ordination formulation in flexible beams with large de- formation[J].Journal of Shanghai Jiaotong Universi-
ty，2005，39(5):827—831. [11]陈思佳，章定国，洪嘉振.大变形旋转柔性梁的一种
高 次 刚 柔 耦 合 动 力 学 模 型 [J]. 力 学 学 报 ，2013，45 (2):251—256.
Chen Sijia，Zhang Dingguo，Hong Jiazhen.A high-or- der rigid-flexible coupling model of a rotating flexible beam under lagre deformation [J]. Acta Mechanica Sinica，2013，45(2):251—256.
[12]张志刚，齐朝晖，吴志刚.基于曲率插值的大变形梁 单 元 [J].应 用 数 学 和 力 学 ，2013 34(6):620—629.
Zhang Zhigang，Qi Zhaohui，Wu Zhigang.A large de- formation beam element based on curvature interpola- tion[J].Applied Mathematics and Mechanics，2013， 34(6):620—629.
[13]Garcíd D， Valverde J， Dominguez J. An internal damping model for the absolute nodal coordinate for- mulation [J]. Nonlinear Dynamics， 2005， 42 (4 ): 347—369.
[14]Simo J C，Vu-Quoc L.On the dynamics in space of rods undergoing large motions—ageometricaly exact approach[J].Computer Methods in Applied Mechan- ics and Enginering，1988，66(2):125—161.
[15 ] Zhao Z ， Ren G . A quaternion-based formulation of Euler-Bernouli beam without singularity[J].Nonlin-
ear Dynamics，2012，67(3):1 825—1 835.
Rigid-flexible dynamics analysis of a large deformation beam element based on interpolation of strains
ZHANG Zhi-gang1 ，QI Zhao-hui1 ，WU Zhi-gang1，2 (1.State Key Laboratory of Structural Analysis for Industrial Equipment，Dalian University of Technology，
Dalian 116024，China; 2.School of Aeronautics and Astronautics，Dalian University of Technology，Dalian 116024，China)
Abstract:The “dynamic stifening”phenomenon in flexible multi-body system dynamics is due to the deformation coupling. The first-order approximation model has ben sucesfuly applied in the rigid-flexible coupling modeling of smal deformation. However，it is found necesary to consider more deformation coupling efects in lager deformation cases.In this paper，the beam bending curvature and axial strain are selected as the element parameters.Then the recursion formulations are obtained
for the kinematic parameters of two end nodes of an element based on theories of large deformation and finite rotation.A planar beam element used for large deformation rigid-flexible dynamics analysis is proposed，which can automaticaly take into acount the “dynamic stif fening terms”.Final ly，the validity and ef fectivenes s of the proposed algorithm are verified through some nu- merical examples which involve the large deformations and rigid-flexible dynamics of beams.
Key words:rigid-flexible coupling;beam element;strain interpolation;dynamic stifening
作 者 简 介 :张 志 刚 (1984— )，男 ，博 士 研 究 生 。 电 话 :15326175369;E-mail:[email protected] 通 讯 作 者 :齐 朝 晖 (1964— )，男 ，教 授 ，博 士 生 导 师 。E-mail:[email protected]
Reviews
There are no reviews yet.