系統平衡方程式與收斂準則

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

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

例題四：固端梁受自重，端點受扭角及軸向壓縮

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

P

及

M

為

B

端

對梁之反力，此細長梁斷面性質與文獻

[2]

使用之細長梁相同，其值列於表 三，本例題採用

Engineering strain

分析，考慮三階段的負荷作用，第一階 段：兩端為固定端，考慮梁的自重，第二階段：允許

B

端的軸向轉角，在

L

T

B

端施加一固定軸向轉角



，第三階段：限制

B

端軸向轉角，允許

B

端軸 向 位 移 ， 在 B 端 施 加 軸 向 位 移  。 此 細 長 梁 其 楊 氏 係 數

2

3

/

10

57 N mm

E

  ，剪力係數

G



20



10

³

N / mm

²，在考慮自重時將自重

mm

N / 709

.

4

⁵

視為一個均佈負荷

q g



10

^ 作用在

X

₂^G的負方向，考慮軸 向轉角







 =0

、

1

、

2

、

3

、

4

、

5 

六種狀況，為了防止超出材料的彈性範圍，

在扭轉角大時，將使用較大的

L

_T ，由於本例題當



越大時越不易收斂，且 在考慮自重時，會更加不易收斂，因此隨著



的增加，將增加元素數及加 大

( 2.9 .2)

式中的容許誤差

e

_tol來進行分析。本例題在考慮自重

L

_T

 650 mm

，







5

，

e

_tol 

1



10

^6，元素數即使用到

600

個依然不易收斂，嘗試使用較 大的容許誤差

e

_tol，而經過分析

e

_tol 

1



10

^6與

5



10

^4所得的曲線幾乎完全 重疊，因此在此狀況下容許誤差將使用

e

_tol 

5



10

^4。以下為在不同軸向 轉角



所使用的

L

_T、元素數及

e

_tol：

0 

：

L

_T 

300 mm

，

100

個元素，

e

_tol 

5



10

^7

1 

：

L

_T 

400 mm

，

200

個元素，

e

_tol 

1



10

^6

2 

：

L

_T 

400 mm

，

200

個元素，

e

_tol 

1



10

^6

3 

：

L

_T 

500 mm

，

200

個元素，

e

_tol 

1



10

^6

4 

：

L

_T 

500 mm

，

600

個元素，

e

_tol 

1



10

^6

5 

：不考慮自重下

L

_T 

650 mm

，

600

個元素，

e

_tol 

1



10

^6



圖十九

-

二十九分別為

0-5 

之

 / L

_T

-

曲線圖及

-

曲線 圖，由各轉角之

-

圖中可以看出文獻

[2]

其理論之結果在受到一

極微小的軸向強制位移後，軸向力會瞬間大幅下降，而使

-

曲

線產生明顯的轉折處，而文獻

[2]

的實驗結果軸向力下降則較為緩慢，其

-

曲線相較文獻

[2]

之理論結果較為平滑，且文獻

[2]

實驗結果的 軸向力皆小於文獻

[1]

的理論結果，且隨著

P

cr

P / M / M

_cr

L

 /

L

T

 /

T

P / P L

T

 / P / P

_cr

cr

L

T

 / P / P

_cr



的增加兩者間的差距也越大。

本例題在

 = 0

不考慮自重下所得的結果皆與文獻

[1]

的理論解幾乎重 合，而在



0不考慮自重所得的

 / L

_T

-

曲線與文獻

[1]

的理論解幾乎 平行但兩者間會有差距，文獻

[2]

的理論結果會比較大且隨著

P

cr

P /



的增加兩者 間的差距也越大。由各轉角之

 / L

_T

-

曲線圖可看出，本例題在考慮 自重後

 -

的曲線會比沒有考慮自重時平滑，與文獻

[2]

實驗結果 較為相似，而本例題考慮自重與不考慮自重主要差異在



不大的時

後，當 增加後兩者會漸漸重合。

P

cr

/ P

L

T

/

L

T

/

P

cr

P /

L

T

/



本研究所得結果和文獻

[2]

的理論結果幾乎相同，都會與文獻

[2]

實驗結 果有差距，本文推測可能是實驗儀器造成的誤差，而文獻

[2]

中並未將其實 驗儀器介紹清楚，因此其誤差不知從何而來，本研究將曲線向下平移

(

即個 轉角



下

 / L

_T

- P / P

_cr圖中之實線，其平移的值列於圖中

)

，當



越大平移 量也越大，由圖二十八可看出在

 / L

_T

- P / P

_cr曲線平移後

 = 5 

可以幾乎重 合，而

 = 4 

之曲線也相當穩合

(

圖二十六

)

，但當軸力為壓力時

(  =0- 3 ) 

， 文獻

[2]

實驗解的線型與本題考慮自重的線型的差異隨著



的減少而逐漸增 加

(

圖十九、二十、二十二、二十四

)

，文獻

[2]

中提到實際的邊界條件與理

本研究所得結果和文獻

[2]

的理論結果幾乎相同，都會與文獻

[2]

實驗結 果有差距，本文推測可能是實驗儀器造成的誤差，而文獻

[2]

中並未將其實 驗儀器介紹清楚，因此其誤差不知從何而來，本研究將曲線向下平移

(

即個 轉角



下

 / L

_T

- P / P

_cr圖中之實線，其平移的值列於圖中

)

，當



越大平移 量也越大，由圖二十八可看出在

 / L

_T

- P / P

_cr曲線平移後

 = 5 

可以幾乎重 合，而

 = 4 

之曲線也相當穩合

(

圖二十六

)

，但當軸力為壓力時

(  =0- 3 ) 

， 文獻

[2]

實驗解的線型與本題考慮自重的線型的差異隨著



的減少而逐漸增 加

(

圖十九、二十、二十二、二十四

)

，文獻

[2]

中提到實際的邊界條件與理

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