座標系統描述

本 研 究 是 使 用 共 旋 轉 有 限 元 素 法 (co-rotational finite element formulation)，將梁分割成若干個兩個節點的梁元素。為了描述系統的運動、

元素的變形、邊界條件，本文中共定義了三套直角座標系統：

(1) 固定總體座標系統(圖二)，

X

_i^G

( i  1 , 2, 3)

系統的節點座標及方位，其他座標系統之座標軸的方向餘弦，皆 在此座標系統中定義。

(2)

元素座標系統

(

圖二

)

，

x

_i

( i  1 , 2, 3)

此座標系統附屬在每一梁元素上，其原點位於該元素的節點

1

上，

軸通過該元素的兩端節點

(1

，

2)

， 軸與 軸在元素變形前與斷 面的主軸方向一致，而元素變形後的 軸與 軸，可由該元素未 翹曲的兩端斷面的方位來決定，是分別將位於節點

1

^，

2

^變形後 的斷面繞一個與該斷面之法線及 軸垂直的旋轉軸旋轉一角度使

斷面之法線方向與 軸方向一致(此時並不考慮斷面之翹曲變

形，否則斷面的法線方向無法定義)，然後再以兩斷面主軸方向的 角平分線作為 軸及 軸的方向。梁元素的變形、節點內力、以 及元素剛度矩陣，皆在此座標系統中建立。

x

1

x

₂

x x

3

2

x

3

x

1

x

1

2

x

x

₃

(3)

元素斷面座標系統

(

圖二

)

，

x

_i^S

( i  1 , 2, 3)

此座標系統與元素的斷面一起平移和旋轉，其原點剛接於未翹曲 斷面的剪力中心， 軸為未翹曲斷面的法線方向， 軸與 軸分 別與未翹曲斷面的主軸重合。元素的變形是由斷面座標相對於元 素座標的旋轉來決定。

x

₁S

x

₂^S

x

₃^S

本文中以符號

 

代表行矩陣。總體座標系統X^G 

S 3

{ }與元

素座標x =



,元素斷面座標 ={ }的關係可表示如下:

G G

G

X X

X

₁

,

₂

,

₃

3



2 1

, x , x

x

x^S

x

₁^S

, x

₂^S

, x X

^G

 A

_GE

x

X

^G

 A

_GS

x

^S

(2.3.1)

其中 、 分別代表元素座標、元素斷面座標對於固定總體座標系統

的方向餘弦矩陣。

AGE A_GS

2.

量

4

旋轉向量

本文中使用旋轉向量來表示一個有限旋轉，如圖三所示，一向量b受 到一旋轉向



a的作用而轉到一個新的位置b，向量 與 之間的關係可

成 :

b b

表示

b   cos  b  ( 1  cos  )( a  b ) a  sin  ( a  b )

(2.4.1)

其中符號．與分別代

表向量的內積與外積，



表逆時鐘方向旋轉角， 表

轉軸的單位向量。

形可以由其形心軸的位移、截面的翹曲

(warping)

及其截面

之形心軸上的對應點。在旋轉座標上

Q

點的變

a 旋

2.5 Euler

梁的變形描述

本文在旋轉座標上描述梁元素的變形，由

(2.2)

節中的基本假設可知

Euler

梁元素的變 的旋轉來描述。

本研究採用

[23]

中之梁的變型機制，如圖二中

Q

點為梁元素中的任意

點，

P

點為

Q

點在同一斷面 形前後位置可以表示如下

:

r

₀

 x e

₁

 y e

₂

 z e

₃

(2.5.1)

r  x

_p

e  v e

₂

 w e

₃



_,s

ω e

₁^S

 y e

^S₂

 z

^S₃

(2.5.2)

其中

θ

₁

e

1

x 、 y 、 z 為變形前 Q 點在元素座標 x (

_i

i

1 3,2, )上的座標， x 亦為 點 變形前 軸的座標，

y 、

P

x

1

z 亦同時是 Q 點在 x

₂^S與

x

₃^S軸的座標。

x

_p(

x

)、 (x

v

) 以及

w

(x)分別是變形後 點在元素座標P

x (

_i

i

1,2,3)上的座標，

v

(x)、

w

(x)亦 為 P 點在

x 及

₂

x 軸方向的位移，

₃

θ s



θ

,s

  ¹

1 是梁斷面沿變形後形心軸的軸向

扭轉率，

θ

₁(

x

)為形心軸的扭轉角，s 為變形後形心軸的弧長，







(

y

,

z

)代 表等斷面梁的聖維南翹曲函數，

e

及

e

^S(

i

1,2,3)分別為

x

與

x 軸的單位向

^S

量，梁變形後形心軸的單位切線向量可表示為

旋轉向量連續作用於單位向量

e

_i (

i

1,2,3)來決定[18]、[23]：

cos }

₃



₂



₃



₁

( 

₂²



₁²

)] 

_1,



₂



數(shape function)。

N

i

將(2.5.19)至(2.5.21)式代入(2.5.18)式整理可得



₀ ) 分。本文中採用了 Green strain 及工程應變(Engineering strain)來描述梁的變 形，但本章中的推導僅列出用 Green strain 的推導，工程應變的推導在[17]

中有詳細的推導，兩種應變推導的節點內力中的軸力及扭矩有些微的差

x y  z

)

)

元素對應於虛位移

 q

_的虛應變，

E 為楊氏係數，G 為剪力模數，V 為梁元

27.1 等效節點變形內力 的推導

f

_

由(2.7.1)、(2.7.2)、(2.7.6)、(2.5.20)－(2.5.23)式及(2.5.29)式可得



^{ }







0 P Q

Q F Q

Ψ( ,



) ( ,

λ

_{d d})

λ

_f  (2.9.1)

其中Ψ 為系統的不平衡力，F 為系統節點內力，Q 為系統的節點位移，

為參考位移負荷向量，P 為參考力負荷向量，

Q

d



f 為力負荷參數，



_d 為位移 負荷參數，其中 可由(2.7.13)元素節點變形力求得。 F

本文以不平衡力向量 的 weighted Euclidean norm 作為迭代時的誤 差度量，而且收斂準則表示為

tol r

N e

e  

P



(2.9.2)

其中

N

表離散系統的自由度數， 為端點反力， 是一給定的容許誤差 值。

P

r

e

_tol

第三章 數值計算方法與程序

本文解非線性平衡方程式(2.9.1)式的數值計算方法是基於牛頓－拉福 森(Newton-Raphson)法配合弧長控制法(arc length control)的增量迭代法 [20]。本文中以系統切線剛度的行列式值為零當作挫屈的準則，為了求得挫 屈負荷，本文採用二分法[21]，決定增量位移的長度，以求得系統切線剛度 矩陣之行列式值為零的平衡位置。為了求得次要平衡路徑，本文中在平衡 路徑的第一個挫屈負荷分歧點加入一個與第一挫屈模態向量方向一致的擾 動位移。

3.1 增量迭代法

若第 I 個增量的平衡位置為已知，則在此位置的系統切線剛度矩陣 可以求得，且第

K

T

I +1個增量的初始增量位移向量 Q

 ，可利用尤拉預測值 (Euler predictor)求得

T

, Q Q   

 (3.1)

其中 為初始增量負荷參數，



Q_T K_T^1P_ref 為參考負荷向量 的切線 解。

Pref

 可利用下式求出

 [20]

, ) (Q_T^t Q_T ^1/2











^(3.2)

其中正負符號之決定方法如下：若第 I 與 I -1個增量收斂時，系統切線剛度 矩陣之行列式值同號，則 的正負符號和第 I 個增量時相同；若異號則符



號相反。表第

I +1個增量的增量弧長，其值可以如下決定

, )

( J

_D

J

_I ^1/2



_I

  

 (3.3)

其中

J

_D為給定的期望迭代次數，

J

_I為第

I

個增量，迭代至平衡所使用的迭 代次數，

 

_I為第

I

個增量的增量弧長。

本文中第一個增量的增量弧長

 

₁是由下式決定

r

c

I R

max 0 max 1



R

 (3.4)

上式中

R

_max為給定的參考自由度之最大位移量， ₀ 為參考負荷向量 _ref 作用下的系統線性解 ₀的

Euclidean norm

_max為給定之最大增量次數

R

，

P

R

I

，

r

c 為R 在參考自由度的分量的絕對值。₀

在平衡迭代時增量位移向量 Q 及增量負荷

 

已知，由 可求得梁 結構新的變形位置。再利用

2.7

與

2.8

節的方法，求得元素座標上的節點內力 及剛度矩陣。而對應此位置的負荷參數為

Q





 

_I

 

，其中



_{I 為第 I 個增} 量達平衡時的負荷參數，

 

即增量負荷參數。當系統內力及外力求得後，

不平衡力量  向量可由

(2.9.1)

式求得。若

(2.9.2)

式的收斂準則不能滿足，則 利用定弧長控制法，求得一位移修正量 Q



與負荷參數修正量



，並加入前 一次迭代的Q與

 

Q

中，而得一新的增量位移向量與增量負荷參數，再進 行下一次的迭代。



與



可由下列二式決定

(3.5) )

1

( P

K

Q 

 

_T^

  

(3.6) )

( )

2

 (  Q   Q  Q   Q

 

^t

以上之迭代計算過程一直重覆至滿足

(2.9.2)

式的收斂準則為止。

3.2

二分法

利用

3.1

節的增量迭代法可以求得結構之主要平衡路徑。在每個增量 的迭代收斂時，可以得到該增量在其平衡位置的負荷參數



及結構切線剛度 矩陣的行列式值

D

(



)。令I及

D

(



_I)分別表示第 I 個增量在其平衡位置的



及

D

(



)之值。



_I1及

D

(



_I1)分別表示第 I +1的增量在其平衡位置的



及

) (



D

之值。

 

_I1表示第

I +1個增量的增量位移向量之弧長。若 D

(



_I)大於 零且

D

(



_I1)小於零則可利用以下二分法求得挫屈負荷參數



_NB：

(1) 令

 _

_L



0,

 _

_R

  _

_I1,

^

_L

^{ }

_I ,



_R

 

_I₁，其中下標 L 及 R 表示左界 及右界。

(2) 取

1 2

R L

I



   

 _ ，重作第 I +1個增量迭代，並求得新的



_I1及 )

( _I1

D 

。

(3) 若

D

(



_I1)大於零，則令



_L 



_I1,_L _I1 若

D

(



_I1)小於零，則令



_R 



_I1,_R _I1 (4) 若下列二式挫屈誤差準則同時滿足

D

I

e

D

D

^  ) 0 (

) (



₁

(3.7)







 e

I L

R  

1 (3.8)

其中

e

_D及

e

_為給定的容許誤差值

則



_I1即為系統挫屈負荷，否則回到步驟(2)，重新展開下一次二分法迭代。

經由二分法求得挫屈負荷



，再利用系統切線剛度

K

(



)_{計算挫屈模} 態，以下將說明挫屈模態的計算程序

(1) 將

K

(



)分解成一下三角矩陣L 及上三角矩陣 使U K LU，其中L 的對角 線元素值皆為1。

(2) 找出矩陣 U 主對角線元素之絕對值有最小值的行令為第 I 行。

(3) 令模態第 I 個分量



_I

 1

，再將KΘ 改寫成下式 0

第 四 章 數值例題與結果

本章中將在 4.1 節驗證理論之準確性以及程式的正確性，在 4.2 節為探 討文獻[1,2]實驗，並與其實驗結果與理論結果比較。本章中

2

4

2 T

y

cr

L

P  EI



為

固端梁受軸向力的挫屈負荷，

T y

cr

L

M

2.861

 EI

 為固端梁受軸向扭矩的挫

屈扭矩。

4.1

本文理論之準確性以及程式的正確性

為了確定本文使用之理論的準確性以及程式是否正確，本節中分析了

兩個例題並與文獻的結果比較。

例題一：固端梁受端點軸向力及扭矩的挫屈分析

如圖四所示，一長度為 的圓形斷面梁，A 端為固定端，B 端只允許

軸向位移及轉角，本例題有兩階段的力負荷，第一階段：在 B 端施加一固 定軸向外力

L

T

P ，第二階段：在 B 端不斷增加軸向扭矩 M ，此細長梁其楊氏

係數

E

5710³

N

/

mm

²，剪力係數

G

2010³

N

/

mm

²，元素的數目使用 80 個元素，

L

_T

 2000 mm

，斷面半徑

r

0.5

mm

，本例題圓形斷面之斷面性 質列於表一，本例題(2.9.2)式之容許誤差

e

_tol

 1  10

^7，表二為在不同

P

下 所對應的挫屈扭矩 ，圖五為挫屈扭矩

-

軸力關係圖，本文的結果和文獻

[1]

的結果幾乎重合，由圖五可看出在

M

nb

 0

P

時 會近似於 ，而在

時

M

nb

M

_cr

P

cr

P  M 會約等於 0

。

例題二：細長梁受端點軸向扭矩之次要平衡路徑

內，在軸向扭轉越大的時候，將使用越長的梁進行實驗，但大多數的梁都 沒有超過

500mm

。文獻

[2]

實驗首先以鑽機

(drill)

夾鑽頭的方式將細長梁兩 端盡量保持軸方向的方式固定於

AB

兩端，在

B

端施加強制轉角



，再施 加軸向強制位移 ，由兩端之感測器量測端點軸向力

P 與軸向扭矩 M 。而

文獻

[1,2]

實驗的結果與文獻

[1,2]

理論的結果有所差異，文獻

[2]

認為此差異 可能是沒有考慮重力影響、沒有考慮摩擦力、沒有考慮剪變形、沒有考慮 材料非線性、及實際的邊界條件與分析的邊界條件不一致等原因造成的，

此 外 當 應 變 不 是 很 小 時 ， 由 工 程 應 變

(Engineering strain)

及 工 程 應 力

(Engineering stress) [17]

與使用

Green strain

及

Second Piolla Kirchhoff stress

[11, 14]

推導出的梁之軸向力的二次項有些差異，此差異亦可能影響分析的

結果，所以本章將分析梁之自重、實際的邊界條件與分析的邊界條件不一 致、梁結構之初始不完美

(initial imperfection)

以及採用

Engineering strain

與

Green strain

對分析結果的差異。本章例題三將探討梁受重力以及採用

Engineering strain

與

Green strain

的差異。

例題三：固端梁受自重與端點扭角

如圖十一所示，一長度為 的圓形斷面梁，

A 、 B

端皆為固定端，其

中

P

及

M

為

B

端對梁之反力，

L

T



為端點強制扭角， 為中點側向位移，

本例題有兩階段的負荷，第一階段：考慮梁的自重，第二階段：在

B

端施 加強制轉角

V

m



。本例題分為四種狀況：不考慮自重採用

Engineering strain

分析、不考慮自重採用

Green strain

分析、考慮自重採用

Engineering strain

分析及考慮自重採用

Green strain

分析，梁的長度考慮 、

400

、

500

、

、 五 種 長 度 ， 斷 面 半 徑

 300 L

T

600 700 mm r



0 . 5 mm

， 此 細 長 梁 其 楊 氏 係 數

質與文獻

[2]

使用之細長梁相同，其值列於表三，在考慮自重時將自重視

Engineering strain

分析及採用

Green strain

分析的兩種狀況在不同長度 及 扭角

L

T



下對應之 、 和 ，表五為不考慮自重下採用

Engineering strain

分析及採用

Green strain

分析兩種狀況在不同長度 及扭 角

考慮自重以及採用

Engineering strain

與

Green strain

對軸向扭矩影響十分的 小。圖十二、十三分別為不考慮自重及考慮自重在

Engineering strain

分析的軸力差異不大，在考慮自重下 越大，則 會

越大。圖十七為考慮自重在不同轉角

L

T

P / P

_cr



下

-

之關係圖，因採用

Green strain

分析的軸力較大，因此以採用

Green strain

分析的中點軸向位移應較 小，而圖十七也符合預期的結果。。

L

T

V

_m

因

Green strain

的物理意義不明確，而

Engineering strain

有明確的物理 意義，所以本章中將以

Engineering strain

的結果與文獻的結果比較。

在文檔中 三維挫屈梁之非線性分析 (頁 14-0)