國
立
交
通
大
學
機械工程學系
碩 士 論 文
以有限元素法分析旋轉傾斜尤拉梁的穩態變形與
自由振動
Steady state and free vibration analysis of a rotating inclined
Euler beam by finite element method
研 究 生：周裕淳
指導教授：蕭國模 博士
以有限元素法分析旋轉傾斜尤拉梁的穩態變形與自由振動
Steady state and free vibration analysis of a rotating inclined
Euler beam by finite element method
研 究 生： 周裕淳
Student： Yu-Chun Zhou
指導教授： 蕭國模 博士 Advisor： Dr. Kuo-Mo Hsiao
國 立 交 通 大 學
機 械 工 程 學 系
碩 士 論 文
A Thesis
Submitted to Department of Mechanical Engineering
College of Engineering
National Chiao Tung University
in Partial Fulfillment of the Requirements
for the Degree of
Master of Science
in
Mechanical Engineering
September 2009
Hsinchu, Taiwan, Republic of China
以有限元素法分析旋轉傾斜尤拉梁的穩態變形與自由振動
Steady state and free vibration analysis of a rotating inclined Euler beam by
finite element method
研究生：周裕淳 指導教授：蕭國模博士
國立交通大學機械工程學系碩士班
摘要
本研究主要利用共旋轉有限元素法推導旋轉傾斜尤拉梁的運動方程
式，探討設定角為
具不同傾斜角之等速旋轉傾斜尤拉梁的穩態變形及以
穩態變形為平衡點的自然振動頻率。
D0
本文利用非線性梁理論的一致線性化、d’Alembert 原理和虛功原理在當
前的旋轉元素座標上推導梁元素的節點變形力、節點慣性力、元素剛度矩
陣、元素向心力剛度矩陣(
centripetal stiffness matrix)，
元素質量矩陣(mass
matrix)，元素陀螺矩陣(
gyroscopic matrix
)。本研究中將變形參數對時間的微
分視為擾動量，故僅取到一次項，元素節點變形內力取到變形參數的二次
項，元素節點慣性內力僅取到變形參數的一次項且忽略變形參數與變形參
數對時間的微分的耦合項。將系統的非線性運動方程式中對時間的微分的
項去掉即為系統的穩態平衡方程式，將系統運動方程式用泰勒級數在穩態
變形的位置展開，取到一次項，即為旋轉梁微小振動的運動方程式。
本文利用基於牛頓法的增量迭代法求出軸向位移及側向位移的穩態
解。旋轉傾斜梁的振動方程式中存在陀螺矩陣，所以其自然振動頻率對應
的振動模態為複變數，其頻率方程式(frequency equations)為一組代數齊次方
程式，該組齊次方程式為一個二次特徵值問題，其係數形成之矩陣的行列
式值為零時的根，即為自然振動頻率。本文以二分法來求行列式值為零時
的根。
本研究以無因次化的數值例題探討傾斜角、無因次轉速、無因次轉軸
半徑及細長比對旋轉傾斜梁無因次側向穩態變形和無因次自然頻率的影
響，本研究還探討旋轉傾斜梁的兩個振動頻率接近時，對應之振動模態的
耦合現象、特徵值曲線轉向(eigenvalue curve veering)及振態交換的現象。
Steady state and free vibration analysis of a rotating inclined Euler beam by
finite element method
Student：Yu-Chun Zhou Advisor：Dr. Kuo-Mo Hsiao
Department of Mechanical Engineering
National Chiao Tung University
Abstract
In this paper a co-rotational finite element formulation is proposed to derive
the equations of motion for a rotating inclined Euler beam with constant angular
velocity. The steady state deformation and natural frequency of the
infinitesimal free vibration measured from the position of the corresponding
steady state deformation are investigated for rotating inclined Euler beams with
zero setting angle.
The element deformation nodal forces, inertia nodal forces, stiffness matrix,
centripetal stiffness matrix, mass matrix and gyroscopic matrix are
systematically derived by consistent linearization of the fully geometrically
non-linear beam theory using the d'Alembert principle and the virtual work
principle in the current rotating element coordinates. In this paper the terms up
to the second order of deformation parameters and their spatial derivatives
corresponding to the steady state element deformations are retained. However,
only the terms up to the first order of deformation parameters, their spatial
derivatives, and time derivatives corresponding to the free vibration of the beam
time derivatives are neglected. The steady state equilibrium equations may be
obtained by dropping the terms of the time derivatives in the equation of motion.
The governing equations for linear vibration may be obtained by the first order
power series expansion of the equation of motion at the position of the
corresponding steady state deformation.
An incremental-iterative method based on the Newton-Raphson method is
employed to solve the steady state deformation. The frequency equations for
free vibration of rotating inclined beam are a set of homogeneous equations. The
natural frequencies may be determined by solving the homogeneous equations
using the bisection method.
Dimensionless numerical examples are studied to investigate the
dimensionless steady lateral deformation and the dimensionless natural
frequency of rotating inclined beams with different inclined angle,
dimensionless angular velocities, dimensionless radius of the hub, and
slenderness ratios. The phenomenon of the coupling of mode shapes, the
phenomenon of eigenvalue curve veering and mode shape changes are also
investigated for rotating inclined beams that has two modes with closely spaced
誌謝
衷心感謝指導教授 蕭國模博士在這兩年期間的指導與教誨，使本論文得
以順利完成，蕭老師在研究上嚴謹的態度，使我受益良多，在此致上最高的敬
意及謝意。也感謝葉孟考老師、蔣長榮及尹慶中老師撥冗擔任口試委員並對本
論文所提出的指正與建議，使本論文能夠更臻完善。
感謝蔡明旭、林育丞、顏宏儒學長們在研究上的協助與照顧，以及學弟蔡
秉宏、林運融、林寬政在各方面的幫助。
感謝父母親、姊姊及弟弟以及所有關心我的親人、朋友對我的支持與鼓
勵，僅以此成果與榮耀，獻給所有關心我的人。
目錄
中文摘要 ... I
英文摘要 ... III
致謝 ... V
目錄 ... VI
表目錄 ... VIII
圖目錄 ... XII
第一章 導論 ... 1
第二章 理論推導 ... 5
2.1 問題描述 ... 5
2.2 基本假設 ... 5
2.3 座標系統描述 ... 5
2.4 Euler 梁的變形描述 ... 7
2.5 梁的應變及其變分、速度、加速度... 11
2.6 元素節點內力之推導 ... 18
2.7 元素剛度矩陣及慣性矩陣之推導... 23
2.8 系統的運動方程式 ... 26
第三章 數值方法及程序 ... 29
3.1 穩態解 ... 30
3.2 振動分析 ... 32
第四章 數值例題 ... 35
4.1 收斂分析 ... 36
4.2 準確性分析 ... 36
4.3 個案分析 ... 37
4.4 旋轉梁之特徵值曲線轉向(Eigenvalue curve veering)分析 ... 40
第五章 結論與展望... 42
參考文獻 ... 44
附表 ... 47
附圖 ... 92
附錄 A ... 135
附錄 B ... 136
表目錄
表一 不同傾斜角與不同轉速的旋轉梁之振動頻率的收斂分析
(
η
=
10
，r
=
1
.
5
)... 47
表二 不同傾斜角與不同轉速的旋轉梁之振動頻率的收斂分析
(
η
=
20
，
r
=
1
.
5
) ... 48
表三 不同傾斜角與不同轉速的旋轉梁之振動頻率的收斂分析
(
η
=
50
，
r
=
1
.
5
) ... 49
表四 不同傾斜角與不同轉速的旋轉梁之振動頻率的收斂分析
(
η
=
100
，
r
=
1
.
5
) ... 50
表五 不同傾斜角與不同轉速的旋轉梁之振動頻率的收斂分析
(
η
=
500
，
r
=
1
.
5
) ... 51
表六 不同傾斜角與不同轉速的旋轉梁之振動頻率的收斂分析
(
η
=
1000
，
r
=
1
.
5
) ... 52
表七 不同傾斜角的旋轉梁之振動頻率的準確性分析
(
70
5
=
k
，
r
=
1
，
η
=
70
)... 53
表八 不同傾斜角與不同細長比的旋轉梁之端點位移的準確性分析
(
r
=
1
) ... 54
表九 旋轉傾斜梁在不同細長比下的振動頻率
(
k
=
0 )... 55
表十 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0 )... 56
表十一 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
0
D) ... 57
表十二 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
5
D) ... 58
表十三 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
10
D) ... 59
表十四 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
15
D) ... 60
表十五 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
30
D) ... 61
表十六 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
45
D) ... 62
表十七 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
60
D) ... 63
表十八 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
75
D) ... 64
表十九 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
0
.
5
，
α
=
90
D) ... 65
表二十 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
0
D) ... 66
表二十一 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
5
D) ... 67
表二十二 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
10
D)... 68
表二十三 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
15
D) ... 69
表二十四 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
30
D)... 70
表二十五 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
45
D)... 71
表二十六 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
60
D) ... 72
表二十七 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
75
D)... 73
表二十八 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
，
α
=
90
D)... 74
表二十九 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
0
D) ... 75
表三十 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
5
D)... 76
表三十一 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
10
D) ... 77
表三十二 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
15
D) ... 78
表三十三 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
30
D) ... 79
表三十四 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
45
D) ... 80
表三十五 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
60
D) ... 81
表三十六 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
75
D) ... 82
表三十七 旋轉傾斜梁在不同細長比與不同轉速下的振動頻率
(
r
=
1
.
5
，
α
=
90
D) ... 83
表三十八 旋轉梁在不同傾斜角與不同轉速下的振動頻率
(
r
=
1
，
η
=
38
) ... 84
表三十九 旋轉梁在不同傾斜角與不同轉速下的振動頻率
(
r
=
1
，
η
=
38
.
5
) ... 85
表四十 旋轉梁在不同傾斜角與不同轉速下的振動頻率
(
r
=
1
，
η
=
39
) ... 86
表四十一 旋轉梁在不同傾斜角與不同轉速下的振動頻率
(
r
=
1
，
η
=
40
) ... 87
表四十二 旋轉梁在不同傾斜角與不同轉速下的振動頻率
(
r
=
1
，
η
=
77
) ... 88
表四十三 旋轉梁在不同傾斜角與不同轉速下的振動頻率
(
r
=
1
，
η
=
78
) ... 89
表四十四 旋轉梁在不同傾斜角與不同轉速下的振動頻率
(
r
=
1
，
η
=
79
) ... 90
表四十五 旋轉梁在不同傾斜角與不同轉速下的振動頻率
(
r
=
1
，
η
=
80
) ... 91
圖目錄
圖一 無傾斜角的旋轉梁結構 ... 92
圖二 具傾斜角的旋轉梁結構 ... 92
圖三 旋轉傾斜梁的上視圖 ... 93
圖四 旋轉傾斜梁的側視圖 ... 93
圖五 元素座標及總體座標關係圖... 94
圖六 具傾斜角的旋轉梁結構(
β
=
0
D) ... 95
圖七 梁的變形圖 ... 96
圖八 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
80
，
r
=
1
，
α
=
5
D) ... 97
圖九 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
80
，
r
=
1
，
α
=
30
D) ... 98
圖十 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
80
，
r
=
1
，
α
=
90
D) ... 99
圖十一 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
100
，
r
=
1
，
α
=
5
D) ... 100
圖十二 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
100
，
r
=
1
，
α
=
30
D)... 101
圖十三 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
100
，
r
=
1
，
α
=
90
D)... 102
圖十四 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
1000
，
r
=
1
，
α
=
5
D) ... 103
圖十五 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
1000
，
r
=
1
，
α
=
30
D) ... 104
圖十六 傾斜旋轉梁的
(a
)
穩態變形
(b
)
軸向位移
(c
)
側向位移
(
η
=
1000
，
r
=
1
，
α
=
90
D) ... 105
圖十七 不同轉速下的第一至第六振動模態
(
α
=
0
D，
r
=
1
，
η
=
50
) ... 106
圖十八 不同轉速下的第一至第六振動模態
(
α
=
5
D，
r
=
1
，
η
=
50
)... 107
圖十九 不同轉速下的第一至第六振動模態
(
α
=
30
D，
r
=
1
，
η
=
50
) ... 108
圖二十 不同轉速下的第一至第六振動模態
(
α
=
90
D，
r
=
1
，
η
=
50
) ... 109
圖二十一 不同轉速下的第一至第六振動模態
(
α
=
0
D，
r
=
1
，
η
=
100
) ... 110
圖二十二 不同轉速下的第一至第六振動模態
(
α
=
5
D，
r
=
1
，
η
=
100
)... 111
圖二十三 不同轉速下的第一至第六振動模態
(
α
=
30
D，
r
=
1
，
η
=
100
) ... 112
圖二十四 不同轉速下的第一至第六振動模態
(
α
=
90
D，
r
=
1
，
η
=
100
) ... 113
圖二十五 不同轉速下的第一至第六振動模態
(
α
=
0
D，
r
=
1
，
η
=
1000
) ... 114
圖二十六 不同轉速下的第一至第六振動模態
(
α
=
5
D，
r
=
1
，
η
=
1000
) ... 115
圖二十七 不同轉速下的第一至第六振動模態
(
α
=
30
D，
r
=
1
，
η
=
1000
) ... 116
圖二十八 不同轉速下的第一至第六振動模態
(
α
=
90
D，
r
=
1
，
η
=
1000
) ... 117
圖二十九 無因次振動頻率-無因次轉速曲線
(
r
=
1
，α
=
0
D) ... 118
圖三十 不同轉速下的第一至第六振動模態
(
α
=
0
D，
r
=
1
，
η
=
38
) ... 119
圖三十一 不同轉速下的第一至第六振動模態
(
α
=
5
D，
r
=
1
，
η
=
38
) ... 120
圖三十二 不同轉速下的第一至第六振動模態
(
α
=
30
D，
r
=
1
，
η
=
38
) ... 121
圖三十三 不同轉速下的第一至第六振動模態
(
α
=
90
D，
r
=
1
，
η
=
38
) ... 122
圖三十四 不同轉速下的第一至第六振動模態
(
α
=
0
D，
r
=
1
，
η
=
39
) ... 123
圖三十五 不同轉速下的第一至第六振動模態
(
α
=
5
D，
r
=
1
，
η
=
39
)... 124
圖三十六 不同轉速下的第一至第六振動模態
(
α
=
30
D，
r
=
1
，
η
=
39
) ... 125
圖三十七 不同轉速下的第一至第六振動模態
(
α
=
90
D，
r
=
1
，
η
=
39
) ... 126
圖三十八 不同轉速下的第一至第六振動模態
(
α
=
0
D，
r
=
1
，
η
=
40
) ... 127
圖三十九 不同轉速下的第一至第六振動模態
(
α
=
5
D，
r
=
1
，
η
=
40
)... 128
圖四十 不同轉速下的第一至第六振動模態
(
α
=
30
D，
r
=
1
，
η
=
40
) ... 129
圖四十一 不同轉速下的第一至第六振動模態
(
α
=
90
D，
r
=
1
，
η
=
40
) ... 130
圖四十二 不同轉速下的第一至第六振動模態
(
α
=
0
D，
r
=
1
，
η
=
77
)... 131
圖四十三 不同轉速下的第一至第六振動模態
(
α
=
5
D，
r
=
1
，
η
=
77
)... 132
圖四十四 不同轉速下的第一至第六振動模態
(
α
=
30
D，
r
=
1
，
η
=
77
) ... 133
圖四十五 不同轉速下的第一至第六振動模態
(
α
=
90
D，
r
=
1
，
η
=
77
) ... 134
第一章
導論
旋轉梁結構在實際上的應用是很重要的，像吊扇、渦輪的葉片、直升
機的旋轉翼、風力發電機的葉片、衛星的支臂、飛機的螺旋槳和機械手臂。
振動分析在旋轉梁結構的設計和分析上扮演很重要的角色，文獻上在這方
面已有很多的研究。[1-16]
關於旋轉梁結構的振動分析可從文獻[1,2]有詳細的探討與回顧，
Schilhansl [3]在考慮離心力，但忽略科氏力的情況下，導出了如圖一所示之
等速旋轉梁振動的微分方程式。Lee 與 Kuo[4]探討了如圖一所示之旋轉
Euler 梁，對其旋轉軸的中心半徑、設定角及轉速對旋轉梁彎矩振動自然頻
率的影響。Yokoyama[5]將旋轉慣量及剪變形、旋轉軸的中心半徑和設定角
合併到有限元素的模式中，探討其對自然頻率的影響。Lee and Lin[6]用線
性梁理論去推導旋轉 Timoshenko 梁之運動方程式，並探討了旋轉速度和質
量慣性矩(mass moment of inertia)的耦合效應、設定角和旋轉速度對彎矩自
然頻率的影響。Eick and Mignolet [7]探討旋轉梁在不同旋轉軸中心半徑與旋
轉梁長度之比值下，其受軸向壓應力挫屈時之臨界轉速。
文獻[3-7]均用線性梁理論推導旋轉梁的運動方程式，且在作其振動分
析時都不考慮科氏力，但均無討論其適當性或影響，在文獻[8]Simo and
Vu-Quoc 提到在分析旋轉結構需要用幾何非線性梁理論(至少取到二次項)
才能適當的計算離心力對彎矩剛度的影響，若用線性梁理論(只取到一次項)
將會產生虛假的彎矩剛度流失，所以文獻[3-7]中推導的旋轉梁之運動方程
式及所求得之振動的自然頻率應是不正確的。
文獻[9,10]利用非線性梁理論的一致線性化、虛功原理和 d’Alembert 原
理在旋轉座標上推導旋轉 Timoshenko 梁正確的線性運動方程式，文獻[9,10]
在分析時考慮了軸向變形及科氏力。旋轉梁的自然振動是指以其穩態解為
平衡點的微小振動，故須先求出其穩態解，除了設定角為
或
外，旋轉
梁之穩態變形是三維的變形，且其自然振動是軸向、側向與扭轉耦合的三
維振動，文獻[9,10]僅分析設定角為
或
之旋轉梁，並僅考慮軸向變形及
一個側方向的位移與旋轉的二維振動，文獻[9]提出一套旋轉梁之自然頻率
的級數解法及計算其自然頻率的數值計算程序，並探討科氏力對旋轉梁之
自 然 頻 率 的 影 響 。 文 獻 [10] 以 文 獻 [9] 提 出 的 方 法 及 細 長 比 很 大 的
Timoshenko 梁模擬旋轉 Euler 梁，文獻[10]發現在低轉速時，科氏力對細長
比很大的旋轉梁的自然頻率影響不大，但文獻[9,10]中並無高轉速的結果，
因在高轉速時，文獻[9,10]的數值方法對細長比很大的旋轉梁無法收斂。
D0
90
D D 0 90D文獻[11]分析如圖一所示，設定角為
或
的旋轉 Euler 梁，利用虛
功原理與 d’Alembert 原理，配合非線性梁理論的一致線性化，在旋轉座標
上推導旋轉 Euler 梁正確的線性運動方程式，並僅考慮軸向變形及一個側方
向的位移與旋轉的二維振動，文獻[11]將旋轉梁分成數段，每段稱為一個元
素，每個元素用一個級數解來表示其自由振動，文獻[11]發現當細長比很大
時，在高轉速下僅用一個元素無法求得正確的自然頻率，需將旋轉梁分成
兩個以上的元素，才能求得精確的自然頻率，但文獻[11]並未探討其原因。
文獻[12]考慮一具雙軸對稱之三維旋轉 Timoshenko 梁，利用共旋轉有限元
素法(Co-rotational finite element formulation)和虛功原理配合非線性梁理論
的一致線性化，推導梁元素節點慣性力與節點變形力。具雙軸對稱之三維
旋轉梁的穩態解包含軸向和扭轉變形，文獻[12]保留軸向和扭轉變形的穩態
解到二次項及扭轉率的三次項，文獻[12]將旋轉梁的運動方程式中的時間函
數去掉求得系統穩態平衡式，再用牛頓法求得穩態解，文獻[12]用泰勒級數
在穩態平衡點將運動方程式一致線性化，求得旋轉梁的振動方程式，文獻[12]
D0
90
D探討旋轉速度和設定角對三維旋轉梁之穩態變形及自然頻率的影響。
文獻[13-16]探討，如圖二所示，具有設定角與傾斜角的旋轉梁之自由
振動。當傾斜角不為零時，若設定角不為
，旋轉梁之軸向位移及側向位
移的穩態解都不為零，但文獻[13]僅考慮軸向的穩態變形對自然頻率的影
響，並未考慮側向穩態位移的影響，因文獻[13-15]忽略了側向位移的穩態
解，故其旋轉梁之自然頻率可能不準確。當傾斜角不為零時，僅有設定角
為
或
時，旋轉傾斜尤拉梁的穩態解及自然振動是二維運動，文獻[16]
用虛功原理與 d’Alembert 原理及幾何非線性梁理論的一致線性化，推導設
定角為
或
之旋轉 Euler 梁的二維運動方程式。當設定角為
時，旋轉
傾斜尤拉梁的側向穩態解為零，且因為科氏加速度的值為零，旋轉梁的軸
向與側向振動不耦合，所以文獻[16]用與文獻[11]相同的方法，求得軸向穩
態變形與自然振動頻率及振態，並探討傾斜角、轉速、轉軸半徑及細長比
對等速旋轉傾斜尤拉梁之自然振動頻率的影響；當設定角為
時，旋轉傾
斜尤拉梁之側向穩態變形不為零，且因為科氏加速度的值不為零，其軸向
與側向振動耦合。文獻[16]假設側向穩態變形為小變形，文獻[16]先求出其
軸向穩態解析解，再用級數解法求得旋轉梁之穩態解。文獻[16]將穩態解代
入運動方程式，用一致線性化求得振動的統御方程式，但文獻[16]並未求其
自然振動頻率及振態。
D90
D0
90
D D0
90
D90
D D0
文獻[16]僅保留到變形的二次項，所以側向穩態變形太大時可能不準
確，且文獻[16]以旋轉梁變形前所受的離心力求其軸向及側向穩態變形，但
旋轉梁所受的離心力為與結構變形位置相關的外力(configuration dependent
load)，當旋轉梁的側向穩態變形不是很小時，必須考慮幾何非線性，才能
得到可靠的側向穩態變形。當旋轉傾斜梁的細長比很大，在高轉速時，其
側向位移的穩態解可能很大，其軸向和側向穩態變形會互相耦合，為高度
的幾何非線性問題，若將旋轉梁分成數段(元素)，用共旋轉法描述旋轉梁的
變形，則可除去梁元素的剛體旋轉，故可以僅保留到變形的二次項，但仍
維持解的精度，不過須以迭代的方式求得正確的穩態變形。如果採用文獻[16]
的級數解，求穩態變形的過程將會很複雜，文獻[12]利用共旋轉有限元素法
成功的求出旋轉梁之軸向及扭轉耦合的穩態變形，故本研究擬採用共旋轉
有限元素法求軸向和側向耦合的穩態變形。本研究的主要目的為以共旋轉
有限元素法探討設定角為
之旋轉傾斜梁的穩態變形及自然振動頻率。本
文在第二章中先以梁變形前的形心軸之長度為獨立變數推導梁元素的變
形，本研究利用文獻[16]中旋轉梁的變形機制，由梁元素在當前元素座標之
位置向量及轉速求得梁元素的加速度，再以虛功原理、有限元素法及非線
性梁理論的一致線性化，在梁元素當前之元素座標上推導節點慣性力和節
點變形力，將元素的節點力轉到總體座標後組合成系統的非線性運動方程
式。本研究將旋轉梁的運動方程式的時間函數項去掉求得系統穩態平衡方
程式，再用基於牛頓法的增量迭代法求出軸向位移及側向位移的穩態解，
將運動方程式在穩態平衡位置用泰勒級數展開，取到一次項，求得旋轉傾
斜梁的振動方程式。假設自然振動頻率存在，可獲得一組代數齊次方程式，
該組齊次方程式係數形成之矩陣的行列式為零時，即可求得旋轉梁以穩態
解為平衡點的自然振動的頻率及其對應的振態。本研究擬探討傾斜角、旋
轉速度、轉軸半徑及細長比等對旋轉 Euler 梁之穩態變形、自然頻率及振態
的影響。
D0
第二章 理論推導
2.1 問題描述
如圖二所示，本文考慮一長度為
具均勻斷面且雙軸對稱之尤拉梁，
其支承端以設定角(setting angle)
TL
β
與傾斜角(inclination angle)
α
剛接在一半
徑為
R
剛性圓柱上，該圓柱以等角速率
Ω
繞其軸心旋轉。本文中所有梁的
位移、變形和振動指的是在一個以等角速率
Ω
繞圓柱中心軸旋轉的旋轉座
標上描述的位移、變形和振動。本文中僅考慮梁的軸向位移，單一個側向
位移及旋轉。以等角速率旋轉的傾斜梁存在著一個含軸向及側向的穩態變
形。本文中所有的振動都是指以該穩態變形為平衡點的振動。本文中考慮
的振動是線性振動，所以由振動造成的位移、速度和加速度都視為是一微
小量(infinitesimal quantity)。
2.2 基本假設
本文對梁元素的推導，做如下的假設：
(1)Euler-Bernoulli 假說成立。
(2)梁元素的形心軸之單位長度伸長量(unit extension)為均勻的伸長。
(3)梁元素的變形與應變皆為小變形與小應變。
2.3 座標系統描述
本研究是使用共旋轉有限元素法(co-rotational finite element formulation
)，將梁分割成若干個兩個節點的梁元素，節點1及節點2為元素的兩個端點。
為了描述旋轉梁系統的運動，本文中使用以下二個座標系統：
(1) 總體座標系統(Global coordinates)
X
1、
X
2、
X
3總體座標系統是以等角速率
Ω
繞圓柱中心軸旋轉，如圖三與圖四所
示，總體座標系統的原點是取在旋轉梁斷面的形心軸與旋轉圓柱的交點(即
O點
)上，其
軸和梁變形前的斷面形心軸一致，其
和
軸是取旋轉梁
變形前的斷面主軸方向，將圓柱的轉軸方向繞
軸逆時鐘方向轉
1X
X
2X
3 1X
β
角即為
和
軸的方向。本文中假設梁只有在
、
方向的變形，本文中旋轉梁
的節點座標、節點位移、節點速度、節點加速度及整個系統的運動方程式
均在此座標系統中定義。
3X
X
1X
2(2)元素座標系統(Element coordinates)
x
1、
x
2、
x
3元素座標系統是建立在每個元素當前的位置上，且以一個等角速率
Ω
繞圓柱中心軸旋轉，如圖五所示，元素座標系統的原點是定義在節點1(即o
點)上，令o點當前的總體座標為(
,
, 0)，
軸的方向為梁元素兩節點連
線的方向，令
軸與總體座標的
軸間的夾角為
oX
Y
ox
1 1x
X
1θ
e，因本文中假設只有
平面的變形，所以
軸與總體座標
軸的方向一致，
軸的方向由右
手定則決定，本文中梁元素的位移、變形、速度、加速度及運動方程式，
均在此座標系統定義。
1x
x
2x
3X
3x
2元素座標系統與總體座標系統關係可表示成
X
=
A
GEx
(2.3.1)
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
−
=
1
0
0
0
cos
sin
0
sin
cos
e e e e GEθ
θ
θ
θ
A
(2.3.2)
其中
x
=
{
x
1,
x
2,
x
3}
，
X
=
{
X
1,
X
2,
X
3}
。本文中以{ }代表行矩陣。
令
Ω
X為旋轉梁的角速度向量在總體座標上的表示式，其分量可表示如下：
Ω
X=
Ω
{
0
,
sin
β
,
cos
β
}
(2.3.3)
由(2.3.1)、(2.3.3)式可得，令
Ω
為旋轉梁的角速度向量在元素座標上的表示
式如下：
Ω
=
{
Ω
1,
Ω
2,
Ω
3}
=
A
GEtΩ
X=
Ω
n
(2.3.4)
n
=
{
n
1,
n
2,
n
3}
=
{sin
β
sin
θ
e,
sin
β
cos
θ
e,
cos
β
}
(2.3.5)
其中
n為旋轉軸的單位向量， (i=1,2,3)為其在元素座標軸 的分量，
n
ix
iβ
為
梁的設定角，
θ
e為梁元素變形後元素座標之
x
1軸與總體座標之
X
1軸的夾角。
因除了
及
外，傾斜旋轉梁的穩態解是三維的，故本文假設僅
適用於
及
，為了推導上的方便，本文在推導過程中仍視為
D0
=
β
D90
D0
=
β
D90
β
變數。
當
時，穩態解僅有軸向位移，[16]用級數解詳盡的探討了其自然振動
頻率。當
時，穩態解有側向和軸向位移，故其所受離心力為變形位置
的函數，文獻上仍沒有其側向位移的非線性穩態解及其自然振動頻率，故
本研究主要探討如圖六所示設定角
之旋轉傾斜梁的軸向、側向位移的
穩態解及自然振動頻率，當
時(2.3.4)式的角速度向量
退化成
D90
=
β
D0
=
β
D0
=
β
D0
=
β
Ω
Ω
=
{
Ω
1,
Ω
2,
Ω
3}
=
Ω
n
(2.3.6)
n
=
{
n
1,
n
2,
n
3}
=
{
0
,
0
,
1
}
(2.3.7)
2.4 Euler 梁的變形描述
本文是在旋轉元素座標上描述梁元素的變形。由2.2節中的基本假設可
知，梁元素的變形可由其形心軸在元素座標上的位移及其斷面的旋轉決
定。本文採用梁變形前形心軸的長度為獨立變數。
2.4.1梁元素之位移
圖七中的P點為梁中的任意點，Q點為P點在形心軸上的對應點，即P點
與Q點位於梁的同一斷面上。在元素座標上，Q點在梁變形前後的位置向量
可分別表示為
{x
,
0
,
0
}
與
{
x
p(
x
,
t
),
v
(
x
,
t
),
0
}
。其中t為時間，
x
p( t
x
,
)
及
v
( t
x
,
)
分別是Q點在
x
1與
x
2軸方向的座標。
P點在梁變形前後的位置向量可分別表示如下
r
0=
{
x
,
y
,
z
}
(2.4.1)
r
=
{
r
1,
r
2,
r
3}
=
{
x
p−
y
sin
θ
,
y
cos
θ
+
v
,
z
}
(2.4.2)
其中
x、
y與
z分別為梁變形前P點在 (i = 1,2,3)軸的座標，
x
iθ
為
軸和形心軸
的切線向量的夾角。
1x
(2.4.2)式之
sin 及
θ
cos 可表示成如下
θ
, 0 0
1
sin
1
1
xv
v
x v
v
s
s x
x
θ
ε
ε
∂
∂ ∂
∂
=
=
=
=
∂
∂ ∂
+
∂
+
(2.4.3)
2 / 1 2)
sin
1
(
cos
θ
=
−
θ
∂
∂
∂
∂
=
∂
∂
=
s
x
x
x
s
x
p p(2.4.4)
0∂
−
1
∂
=
∂
∂
−
∂
=
x
s
x
x
s
ε
(2.4.5)
其中s為圖七中o點到Q點間形心軸在變形後的弧長，
ε
0為形心軸的單位伸長
量。
由(2.4.3)-(2.4.5)式可得
x
px
t
=
u
+
∫
x+
−
v
xdx
(2.4.6)
0 2 / 1 2 , 2 0 1[(
1
)
]
)
,
(
ε
其中
為節點1在
方向上的位移，由元素座標系統的定義，其值為零，但
其變分及對時間的微分並不為零。
1u
x
1由小變形的假設，利用近似式
)
2
1
1
(
]
)
1
[(
+
ε
0 2−
v
,2x 1/2≈
+
ε
0−
v
,2x，(2.4.6)式
可表示成
x
px
t
=
u
+
∫
x+
−
v
xdx
0 2 , 0 1)
2
1
1
(
)
,
(
ε
(2.4.7)
由(2.4.7)式及梁元素的形心軸之單位伸長量為均勻的假設，可以得到形心軸
的單位伸長量
ε
0=
−
+
∫
Lv
xdx
L
L
L
l
0 2 , 02
1
ε
(2.4.8)
l
=
L
+
u
2−
u
1(2.4.9)
其中L為梁元素變形前的長度， l 為梁元素之形心軸變形後的弦長，
為節
點2在
方向的位移。
2u
1x
本文中假設梁元素變形後的形心軸的側向位移
v
( t
x
,
)
為 x 的三次 Hermitian
多項式。因此
v
( t
x
,
)
可表示成
v
(
x
,
t
)
=
{
N
1,
N
2,
N
3,
N
4}
t{
v
1,
v
1′
,
v
2,
v
2′
}
=
N
btu
b(2.4.10)
(
1
)
(
2
)
4
1
2 1=
−
ξ
+
ξ
N
，
(
1
)(
1
)
8
2 2=
−
ξ
−
ξ
L
N
， (2.4.11)
(
1
)
(
2
)
4
1
2 3=
+
ξ
−
ξ
N
，
(
1
)(
1
)
8
2 4=
−
+
ξ
+
ξ
L
N
L
x
2
1+
−
=
ξ
(2.4.12)
u
b= u
b(
t
)
=
{
v
1,
v
1′
,
v
2,
v
′
2}
(2.4.13)
其中
v
j(
j
=
1
,
2
)
是
v 在節點
j
的節點值
，v′ 則是
jx
v
v
∂
∂
=
′
在節點
j
(
j
=
1
,
2
)
之節
點值，
N
i(
i
=
1
−
4
)
代表形狀函數(shape function)。
在小變形的假設下，由(2.4.2)式及近似式
sin
θ
≈
θ
，
22
1
1
cos
θ
≈
−
θ
，可以將
位置向量
r
重新寫成如下
r
=
{
r
1,
r
2,
r
3}
=
)
,
}
2
1
1
(
,
{
x
p−
y
θ
y
−
θ
2+
v
z
(2.4.14)
將(2.4.10)式代入(2.4.8)式整理可得
ε
0)
2
1
(
1
0 2 ,∫
+
=
t a L x av
dx
L
G
u
2
)
1
(
1
b t b a t aL
G
u
+
G
u
=
(2.4.15)
其中
u
a=
{
u
1,
u
2}
(2.4.16)
G
a=
{−
1
,
1
}
(2.4.17)
G
b=
{
G
b1,
G
b2,
G
b3,
G
b4}
=
∫
N
′
bv
,xdx
(2.4.18)
將(2.4.15)式代入(2.4.7)式整理可得
∫
−
+
+
=
t b x x b a t a pv
dx
L
x
x
t
x
x
0 2 ,2
1
2
)
,
(
N
u
G
u
(2.4.19)
⎭
⎬
⎫
⎩
⎨
⎧
−
+
=
2
1
,
2
1
ξ
ξ
aN
(2.4.20)
2.4.2 梁之位置向量的變分
將(2.4.14)式的位置向量
r
變分可表示成
r
δ
=
{
δ
r
1,
δ
r
2,
0
}
=
{
−
y
δθ
+
δ
x
p,
−
y
θδθ
+
δ
v
,
0
}
(2.4.21)
在小變形的假設下，利用近似式
sin
θ
≈
θ
、
(
1
)
)
1
(
1
0 0ε
ε
≈
−
+
，將(2.4.3)式及
變分可得
δθ
=
−
δε
0v
,x+
(1
−
ε
0)
δ
v
,x(2.4.22)
由(2.4.15)式的變分可得
δε
01
(
ta a tb b)
L
δ
u
G
+
δ
u
G
=
(2.4.23)
x 的一次微分可表示成
將(2.4.10)式對
v
,x=
N′
btu
b(2.4.24)
由(2.4.24)式的變分可得
δ
v
,x=
δ
u
tbN
′
b(2.4.25)
由(2.4.19)式的變分可得
px
δ
v
v
dx
L
x
x x x b t b a t a , 0 ,∫
−
+
=
δ
u
N
δ
u
G
δ
(2.4.26)
將(2.4.23)式代入(2.4.22)式整理可得
)
(
)
1
(
)
1
(
0 , 0 0 , , b t b a t a x x xL
v
v
v
ε
δθ
δε
ε
δθ
δ
u
G
δ
u
G
δ
≈
+
+
≈
+
+
+
(2.4.27)
將(2.4.22)及(2.4.26)式代入(2.4.21)式，位置向量的變分 r
δ
可以寫成
1r
δ
=
yv
,xδε
0−
y
(
1
−
ε
0)
δ
v
,x+
δ
u
atN
a tb bL
x
G
u
δ
+
−
∫
xv
xv
xdx
0 ,δ
,(2.4.28)
2r
δ
=
y
(
1
−
ε
0)
v
2,xδε
0−
y
(
1
−
ε
0)
2v
,xδ
v
,x+
δ
v
因為梁元素為小變形，所以
在元素較多時都將趨近於零，故在計算慣性
力時，上式中畫底線的項可以忽略。
xv
,2.5 梁的應變及其變分、速度、加速度
2.5.1梁的應變及其變分
本文中的應變採用工程應變。為了推導上的方便，本文中先推導出
Green strain
ε
ij，再由Green strain求得與其對應之工程應變。Euler梁的Green
strain非為零的應變只有
ε11，可表示成
(
1
)
2
1
1 1 11=
g
tg
−
ε
(2.5.1)
其中
x
∂
∂
=
r
g
1(2.5.2)
將(2.4.7)、(2.4.14)式代入(2.5.2)式，可得
g
1的分量
g
11和
g
12如下
0 , 2 , 0 11
1
2
1
1
ε
ε
+
−
−
+
=
xx xv
y
v
g
2 0 , , , 12
)
1
(
+
ε
−
=
x xx xv
v
y
v
g
(2.5.3)
由小變形的假設，將(2.5.3)式及近似式
(
1
)
)
1
(
1
0 0ε
ε
≈
−
+
，代入(2.5.1)式，且
保留變形參數及其微分到二次項，可得：
11 0 02 , 2 ,22
1
2
1
xx xxy
v
yv
+
−
+
=
ε
ε
ε
(2.5.4)
Green strain
ε
11與對應之工程應變
e 的關係可表示成[17]：
11e
11=
(
1
+
2
ε
11)
1/2−
1
(2.5.5)
將(2.5.4)式代入(2.5.5)式，且保留變形參數及其微分到二次項可得：
e
11=
ε
0−
(
1
−
ε
0)
yv
,xx(2.5.6)
由(2.5.6)式的變分可表示成
δ
e
11=
δε
0+
yv
,xxδε
0−
(
1
−
ε
0)
y
δ
v
,xx(2.5.7)
x 的二次微分可表示成
將(2.4.10)式對
v
,xx=
N ′′
btu
b由上式的變分可得
δ
v
,xx=
δ
u
tbN
b′′
(2.5.8)
將(2.4.23)、(2.5.8)式，代入(2.5.7)式，可得
xx ta a xx tb b
y
tb bL
yv
L
yv
e
=
+
δ
u
G
+
+
δ
u
G
−
−
ε
δ
u
N
′′
δ
11(
1
,)
(
1
,)
(
1
0)
(2.5.9)
2.5.2 梁的速度及加速度
因梁的位置向量是在旋轉元素座標上描述，所以P點的絕對速度在當前
元素座標的分量可表示成
r
r
Ω
v
v
=
{
v
1,
v
2,
v
3}
=
o+
×
+
(2.5.10)
o o o o ov
v
v
Ω
r
v
=
Ω
{
1,
2,
3}
=
×
(2.5.11)
oG t GE o o o or
r
r
A
r
r
=
{
1,
2,
3}
=
(2.5.12)
}
sin
sin
,
cos
sin
,
cos
{
R
α
X
oR
α
β
Y
oR
α
β
oG=
+
−
+
r
(2.5.13)
其中
為 o 點的絕對速度，r
為 P 點對元素座標原點 o 的速度，
為元素
座標原點 o 在總體座標的位置向量的表示式，
為原點 o 在當前元素座標
ov
r
oGr
的表示式。
將(2.3.2)、(2.3.4)式代入(2.5.11)、(2.5.12)式可得：
=
1 ov
n
2r
o3−
n
3r
o2(2.5.14)
=
2 ov
n
3r
o1−
n
1r
o3=
3 ov
n
1r
o2−
n
2r
o1=
1 or
cos
θ
e(
R
cos
α
+
X
o)
−
sin
θ
e(
R
sin
α
cos
β
−
Y
o)
(2.5.15)
=
2
o
r
−
sin
θ
e(
R
cos
α
+
X
o)
−
cos
θ
e(
R
sin
α
−
Y
o)
β
α
sin
sin
3R
r
o=
由(2.3.6)式可知，當
β
=
0
D時，
n
1= n
2=
0
，
n
3=
1
，故(2.5.14)、(2.5.15)
式會退化成
=
1 ov
−
r
o2，
v
o2=
r
o1，
v
o3=
0
(2.5.16)
=
1 or
cos
θ
e(
R
cos
α
+
X
o)
−
sin
θ
e(
R
sin
α
−
Y
o)
(2.5.17)
=
2
o
