穩態解

為了簡潔，在不造成混淆的情況下，本節中將(2.9.2)式之下標s及上標_D 省略，再令P²F_sref^I ，則旋轉梁系統的穩態平衡方程式可改寫成

ψ  F  k

²

P  0

(3.4) 其中P 稱為參考負荷，

k

稱為負荷參數。

本文以基於牛頓法的增量迭代法解非線性代數方程式(3.1)式，求得在 不同無因次轉速 下，旋轉梁的節點位移向量

k

Q。

3.1.1 增量迭代數值計算方法

本文中將選定之最大負荷參數，即最大無因次轉速 等分成數個增量 負荷參數

k

max

 k

，若第

I 個增量的平衡位置為已知，即其位移向量為

、負荷參數 為 為已知，則對應於第

Q

I

k

I

I

1個增量之負荷參數

k  k

_I

  k

的初始增量位移向 量Q，可利用尤拉預測值(Euler predictor)[31]求得

P K Q

( 2



)

 ^1



k

_I

k k

_T

(3.5)

I

T

Q

QQ

K

 

 

(3.6)

其中K_T K 

Ω

²K_Ω為第

I 個增量的平衡位置之系統切線剛度矩陣，

K 和 為系統的剛度矩陣及向心力剛度矩陣。

K

Ω

由 可求得每個元素當前的元素座標及節點變形位移，將其代入

(2.7.13)

及

(2.7.17)

式，可算出元素的節點變形力及慣性力，將元素的節點力轉換到

Q Q

Q 

_I

 

總體座標，可組合得到

(3.4)

式之不平衡力， 依牛頓法，可得位移修正量如下：

再

1









_T

δ Q K (3.7)

不平衡力，

正量 得新的節點位移向

其中 為

K

_T為當前的系統切線剛度矩陣。

將求得的位移修

δ ，加入上次迭代之 Q Q

中，可

量，再進行下一次迭代，此過程一直重複至

(3.4)

式中的不平衡力滿足斂準則為 止。本文以不平衡力 的

weighted Euclidean norm

做為平衡迭代時的誤差度 量，所使用的收斂準則為

P N k

²

 

e

_

(3.8)

其中

 為 的歐幾里德範數

(Euclidean norm)

， 為方程式的數目，

N e

_為一設 定的容許誤差，本文中取

e

_ 

10

^7。

3.1.2

數值程序

量迭代法之數值之數值程序可以分成三個部分：

、參考負荷 。 收斂

(3.8)

本文所使用的增

1.

輸入與計算開始分析所需要的資料

(a)

輸入結構資料及給定外力負荷參數的最大值。

(b)

給定增量數、最大迭代數及收斂時的容許誤差。

(c)

計算增量負荷參數、負荷參數、

Ω K

² Ω

 (d)

用

(3.6)

式計算系統切線剛度

K

_T

 K

解

P

2.

使用迭代法求在已知負荷參數的

Q

(a)

利用 式求初始增量位移向量 。

(b)

將前一個平衡位置的節點位移向量

Q

_I加上Q得到 。

的節點力，再 Q

(c)

由

Q

中萃取元素之節點位移，計算出當前的元素座標及元素 計算

(3.7)

式之不平衡力 。

(d)

檢查

(3.11)

式的收斂準則，若滿足則進行

(e)

；若不滿足，檢查迭代次數，

考負荷。

振動分析

。本文先用本文採用子 空間法

若迭代次數小於給定之最大迭代次數，則利用

(3.10)

式求得位移修正向量

Q

δ ，將當前的節點位移向量 Q

加上

δ 得到一個新的 Q Q

，再回

2(c)

進行 迭代；若迭代次數大於最大迭代次數則停止迭代並印出迭代相關資料。

(e)

檢查增量次數是否大於最大增量次數，若滿足，則完成增量迭代步驟；若 不滿足，則進行步驟

3

。

3.

計算下一次增量所需要的資料

(a)

計算

(3.9)

式中的切線剛度及參

(b)

計算下一次增量的負荷參數。

(c)

回到

2

執行迭代工作。

3.2

本節將說明求旋轉梁自然頻率及振動模態的計算程序

(subspace method)[32]

，求出當無因次轉速

k



0

時的無因次自然頻率

K 及特徵向量

Q，將由

3.1

節之增量迭代法求得對 無因次轉速

k 的穩態

入

(2.9.6)

式中，再以二分法

(bisection method)[33]

解

(2.9.5)

式。

(2.9.5)

式之

( )

應於 解代

K

H

可分解成

H ( K )  LDL

^t，其中L 為下三角矩陣，D為對角 線矩陣。

令

D ( K )



det H ( K )



0 (3.9)

其中

D ( K )

為

H ( K )

的行列式值，其值為D矩陣之對角線元素的乘積。

若 滿足 ，則 為旋轉梁之ㄧ無因次自然頻率。由於H 的 維數隨著元素數目的增加而變大，為了避免其行列式的數值過大，所以本文中將

K

B

D ( K

_B

)  0 K

_B

) ( K

D

做以下的標準化

(normalization)

處理：

) (

) ) (

(

K

0

D K K D

D



(3.10)

其中

K

₀為一參考值。

若 ，其中 為一無因次自然振動頻率，則 ，其中

、 分別為 、

D

中負的對角線元素的個數。如已知 及 ，

則可由二分法求得 。

R B

L

K K

K  

N

R

D ( K

_L

K

B

K

B R

) K

R

L

N

N  K

R

N

L

) ( K

_L

本文解

(3.9)

式所採用的計算程序如下：

設定需要的自然頻率的數目，先用子空間法

(subspace method)

求出當無

因次轉速

k



0

時的無因次自然頻率

K 及特徵向量 。若第 I 個增量的穩態平

衡位置已求出，即其位移向量為

Q

Q

I、負荷參數 為已知，先以前一個無因次 轉速 得到之無因次自然頻率為參考值，設定無因次自然頻率的起始值 及增量

k

I

1

k

I

K

₀

 ，計算

K D ( K

₀

)

中負的對角線元素的數目

N

₀。

(A)

(1)

令

K

_n

 K

₀

 ( n  K ( 1) n



1, 2, 3, ... )

，由

K

_n、

k

_I、K 及

Q

_I計算

D ( K

_n

)

，

一直到相鄰兩個

D ( K

_n

)

中，負的對角線元素的數目不一樣。

(5)將

Θ

除以其分量中絕對值最大的分量



_max，即

max



Θ

Θ

(3.12)

(6)以步驟(A)求得之

K

_B及前一個無因次轉速

k

_I1得到之無因次自然頻率為 參值，設定無因次自然頻率的起始值

K

₀及增量 ，回到步驟(A)繼續求

K

下一個自然頻率及振動模態。

第四章 數值例題之結果與討論

本章將分析不同斷面、細長比、無因次轉軸半徑 r、傾斜角



、設定角



及 預錐角



的三維尤拉梁在不同的無因次轉速 k 下之穩態解、無因次自然頻率

K 及振態。本章中並將比較本文的結果與文獻上的結果以說明本文方法的

正確性與有效性。本文的結果皆考慮抑制翹曲(Warping restraint)對自然頻率 的影響。本章所考慮的斷面如圖七所示之橢圓斷面、矩形斷面、I 型斷面及十字 斷面，各種斷面之斷面常數列於附錄 B 中。梁的長度 與斷面高度比，在橢圓 斷面及矩形斷面指的是 ； I 型斷面上指的是 ，其中 為 I 型斷 面的名義上的(nominal)高度，如

T nom

L d / a

L

_T

/ L

_T

d

_nom

100 10

W

的

d

_nom為10

in

；十字斷面

d



b



t

中 指的是

L

_T

d

，其中

d

為真正的高度 。在橢圓斷面及矩形斷面中

d a

/ 為梁斷面

b

的高度與寬度比。

本章的參數中，

n

_y 

AL

²

I

_z 為附錄 A 中斷面主軸

y

方向的細長比，

y

z

AL I

n 

² 為

z

方向的細長比，因本文中總體座標系統 、梁斷面座標系統 及元素座標系統 (i＝1,2,3）的座標軸方向，在梁未變形時，都是相同的，

方向為 的座標軸方向，

z

方向為 的座標軸方向，所以 及

G

X

i S

x

i

y

x

i

X

₂G

X

₃^G n_y n_z分別為

梁在

X

₁^G-

X

₂^G及

X

₁^G-

X

₃^G平面的細長比。

K

_i (

i 

1,2,3...)為第 個無因次自然頻 率，

i

U

為在

X

₁^G方向之無因次軸向位移，

V

、

W

為在 及 方向之無因次側 向位移， 為

X

方向之軸向扭轉角((2.10.1)式)。當無因次轉速

X

₂G

X

₃^G



1 ₁^G

k 

0時，無因

次振動頻率

K

僅與梁斷面以及細長比有關且其軸向和側向及扭轉的振態不 互相耦合；當

k 

0時由(2.3.3)-(2.3.6)、(2.6.7)、(2.6.8)、(2.6.15)、(2.7.17) 式可知旋轉梁之穩態變形的形態會受無因次轉軸半徑 r、傾斜角



、設定角



及預錐角



的影響，本章中將對應於不同 r 、



、



、



的穩態解型態列於 表一中，由表一可以發現當

r 

0時，



對穩態變形沒有影響；當



 0或90

及



 0時，各種型態的穩態解都存在。

本章中將旋轉梁的自然頻率依其在轉速為零時的振態分為(I = 1, 2, 3, 4)四類：

I = 1 - 轉速為零時，軸向振動

I = 2 - 轉速為零時，在斷面主軸

y

(

X

₂^G)方向的側向振動 I = 3 - 轉速為零時，在斷面主軸

z

(

X

₃^G)方向的側向振動 I = 4 - 轉速為零時，扭轉振動

當轉速不為零時，各種振動會耦合在一起，I (I = 1, 2, 3, 4)代表該振動中第 I 類的 振動有最大的分量。本文中的振動模態圖皆以轉速為零時的振態 I = 1, 2, 3, 4 表 示。

本文中假設穩態解的應變



1，為確保梁的應變在合理範圍內，本文中 將穩態解的允許最大應變定為10^2。本文將限制旋轉梁的最大無因次轉速 轉速

k

，使其最大膜應變(membrane strain)



_0max及兩個側方向的最大撓曲應 變(flexural strain)



_bmax、



_cmax的和儘量不超過10^2。由附錄 C 可知旋轉梁 的最大膜應變和最大撓曲應變為

k

、



、



、



、

r 、 I

_y、

I 的函數，所以

_z 具不同



、



、



、r 以及斷面的旋轉梁有不同的最大無因次轉速 k 。

為了探討穩態變形對旋轉梁之自然頻率的影響，並與[19]的結果比較，本章 中考慮了以下 3 種情況：

A : 考慮穩態變形及保留了(2.7.17)式中所有的項。

B : 不考慮穩態變形，穩態分析時僅考慮(2.7.17)式含 虛底線之項。

C : 穩 態 分 析 時 除 了 將 B 中 第 一 個 含 虛 底 線 項 改 成 外與

B

相同。振動分析時不考慮對應於科氏力的陀 螺矩陣

((2.8.13)

式)。



²

 A  N

_a

a

_Ao1

dx

  a dx

A cos  N

_{a Ao1}





²

N

因

[19]

沒有考慮科氏力且使用了錯誤的軸向離心力，即使用了相當於

C

中的軸向 離心力。為了與

[19]

的結果比較，所以本章中考慮了情況

C

。本章中

A N、 B N、

C

分別表示情況

A

、

B

、

C

，使用 N 個元素的結果。

4.1

個案分析

本節中將探討具不同無因次轉軸半徑r 、傾斜角



、設定角



、預錐角



之 三維旋轉尤拉梁在不同斷面、不同無因次轉速k 之下的穩態解、無因次振動頻率

K 及振態。

本節中首先分析文獻

[19]

的例題，本文的分析使用橢圓斷面；接著本文 分析

.1.1

橢圓斷面旋轉梁的穩態變形及自然頻率

文獻

[23,24]

的例題，其斷面為矩形斷面。本節中將使用的各種斷面梁的

細長比列於表二。

4

表三至表五為橢圓斷面

a / b  10

、

L

_T

/ a  50

、

r

1、

k  0 . 005

、設定角



 0



、

45 

、

90 

、預錐角



 0、

3 0 

、 態解及自 表 ₀



60

及不同元素數目 N 的穩

然頻率， 中



為穩態解之最大膜應變，



_b及



_c為 y 及 z 兩方向之最大撓曲應變，

U

tip為在

X

₁^G 向之無因次端點軸向位移，

V 、

_tip

W

_tip為在

X

₂^G及

X

₃^G方向之無因 點側向 移， ₁

( 1 )

方

次端 位



為梁自由端上的扭轉角，



₁

( 0 . 1 )

為梁在 個長度位置

的扭轉率，

K

_i

( i =1

表五中可以發現

C10

的結果與文獻

19]

的結果非常接近，這是可預期的，因情況

C

與文獻

[19]

用了相 同錯誤的離心力。在

 0 

0.1

~7)

為前七個無因次自然頻率。由表三至

[



時，

A20

、

A50

、

B10

及

C10

的自然頻率很接近，這 也是可預期的，因



 0時，穩態解僅有軸向變形，各種情況的離心力幾乎相同。

而當



增加時，各種情況之穩態解跟

I = 2

、

3

之自然頻率的差距也增加，這應是 離心力的影響。表三至表四中

A20

及

A50

之結果的差距在



增加時也跟著增加，

當



 90、



 60時，

A20

無法收斂，故表五中使用

A50 A8

0，由表五中可 發現兩者的答案接近。由表三可知在



及

  0

、

  60 

時，

A20

及

A50

都發生了 發散不穩定

(divergence instability)

的 ，而 及

C10

都還沒發生，這可 能是因為

B10

及

C10

沒有考慮側向位移及軸向的拉伸變形。

圖八至圖十三為橢圓斷面

10

情形

B10

/ b 

a

、

L

_T

/ a  50

、

r



1

，分別在設定角

0

、

45

^、

90

^及不同預錐角不同轉速下的第一自然頻率曲線圖與第二自然 線圖 圖 發現





0

^、

45

^、

90

^之曲線在

0





頻率曲 ， 中可

k

 時有共同的起點。由

圖八可知在





0

、





0

時，

K

₁隨著轉速增加而增加；在





30

時，頻率

K

₁ 曲線變得平緩；在





45

時，轉速增加到

k



0 . 004785

時

K

₁

0

；在





60

 時，轉速增加到

k



0 . 00 38 9

時

K

₁

 0

。由圖九可知在



 0、



 0 時，

著轉速增加而增加 

2隨

K

；在





30

時，自然頻率

K 曲線變得

₂ 平緩；在





45

時，自 然頻率

K 曲線變得有下滑的傾向。從圖九可以發現在

₂





0

時，當



越大，曲 線的斜率就越小。由圖十與圖十二得知在





45

、

K 都隨著轉速增

₁ 加而增加，且當同一轉速下，

90

時，



越大則

K

₁越小。從圖十 可一 以發現在



 45時，

當



越大，

K 就越快從

₂

I = 2

轉為

I = 4

，且當同一轉速下，



越大則 。

由圖十三可 

2越

K

小

知在





90

、





0

時，

K 隨著轉速增加而增加；在

₂





3 0

時，

K

₂ 隨著轉速增加到

k



0 . 00225

而

；在

增加， 轉速再增加，

K

₂會從

I = I =

且隨著轉速增加而 

但

2

轉為

4

，

減少



 45 時，

K

₂隨著轉速增加到

k



0 . 00175

而增 加，但轉速再增加，

K 會從

₂

I = 2

轉為

I = 4

，且隨著轉速增加而減少；在





60

 時，

K 隨著轉速增加到

₂

k



0 . 0015

而增加，但轉速再增加，

K 會從

₂

I =

I = 4

， 隨著轉速增加而減 十三可以發現在 

2

轉為 且 少。從圖



 90 時 當，



越大，

K

₂

就越快從

I = 2

轉為

I = 4

，且當同一轉速下，



越大則 。

圖十四至圖十六為橢圓斷面

10

K

2越小

/ b



a

、

L

_T

/ a



50

、

r



1

，在不同預錐角不





0



同轉速下的位移分布圖。由圖十四可知在 時，旋轉 僅有一個側向位移

同轉速下的位移分布圖。由圖十四可知在 時，旋轉 僅有一個側向位移

在文檔中 以有限元素法分析具預錐角之三維旋轉傾斜尤拉梁的穩態變形與自由振動 (頁 70-0)