系統運動方程式

Ψ 是不平衡力，F 是系統慣性力內力向量，^I F 是系統變形力內力向量，P^D

是系統外力向量。由(2.6.19)式中可知，慣性力可拆成速度項及加速度項 0

P F F

F

Ψ  ^I1 ^I2  ^D   (2.9.2)

I1

F 為慣性力加速度項，F 為慣性力速度耦合項^I2

本文中平衡迭代所使用的收斂準則為

e_tol

e ΨN 

其中N為系統方程式的數目， Ψ 不平衡力的Euclidean Norm，e_tol是容許 誤差值。

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

本文使用顯積分法及隱積分法求解非線性運動方程式，並比較兩種積 分法解三維梁之非線性動態反應的效率及準確性，因此在本章中將介紹兩 種 不 同 的 數 值 計 算 之 方 法 ， 隱 積 分 法 使 用 Newmark 直 接 積 分 法 和

Newton-Raphson增量迭代法來求解非線性動態平衡方程式，然而此法應用

在受拘束的系統時，會有不穩定的現象發生，若加入 Numerical damping 則可穫得改善，如 Newmark-β damping scheme 及 Hilber-Hughes-Taylor

scheme [25]。故本文亦探討上述兩種具damping效果的數值方法。而本文

中顯積分法使用中央差分法。

本文中探討t_n時刻的動平衡方程式時，初始時間t₀ 0，而t_n nt， 其中 為時間增量，t n為增量次數

3.1 Newmark 直接積分法

3.1.1 Newmark 直接積分法

Newmark 直接積分法乃假設在時間t_n時(n0)，滿足動平衡方程式

(2.9.1)的平衡位置，速度U_{n和加速度U}_{n為已知，則在}t_n1時刻之平衡位置、

速度U _n1和加速度U_n1可由下述迭代過程得到。

1. 令在t_n1時刻的初始增量位移猜測值U0

n n

n t U tU U

U  1) 

2 ( 1 1

1

2

0 1  

 

 

     (3.1.1) ]

) 1

[( ₁

0_n1U_n t  U_n  U_n

U      (3.1.2)

2. 令U_n1U⁰_n1，U_n1U ⁰_n1

3. 由U(或U)及上次迭代後的變形位置得到此次迭代後的變形位置，再 算出系統的節點變形力F_n^D_1。由U_n1、U _n1及最新的變形位置算出系統 的節點慣性力F_n^I_1。由在t_n1時刻的外力P_n1、F_n^D_1、F_n^I_1、及系統動平 衡方程式(2.9.1)算出系統不平衡力

1 1

1

1   

  _n^I  _n^D  _n

n F F P

Ψ (3.1.3)

4. 檢查 Ψ_n1 是否滿足收斂準則，若滿足則迭代停止，本文採用的收斂準 則為

n etol

N ^ Ψ 1

etol為給定的容許誤差，N為方程式的數目

5. 若不滿足，則由

U[Kˆ]^1Ψ_n1 (3.1.4)

得一增量位移修正量U，其中

] [ ] [ ]

1 [ ˆ ]

[K ₂ M C  K

 

 

t  t



 (3.1.5) 為系統有效剛度矩陣(effective stiffness matrix)，[M]、[C]、[K ]

表系統之質量、陀螺與剛度矩陣(其推導詳見(附錄H))。

6. 令

U U

U 

 ²

0 1

1 1

n t

n_   _  

 (3.1.6)

U U

U 





n t

n_  _1 

0 1 

 (3.1.7) 7. 回到 2.

Newmark直接積分法中，若取 0.25 0.5則稱Newmark trapezoidal rule，因Newmark trapezoidal rule[25]不會造成振幅衰減無法除掉高頻的振 動，為了使 Newmark直接積分法能造成振幅衰減，以除掉不必要的高頻振 動，和必需滿足下式[26]：

)2

1 4(

1 

  

) 2 1 2(

1 

   (3.1.8) 上式中 0。

Newmark 直 接 積 分 法 中 若  和 滿 足(3.1.8)式 ， 則 稱 為 Newmark- β damping scheme[25]。本文中取 0.05、 0.275625、 0.55。

3.1.2 Hilber-Hughes-Taylor scheme

Hilber-Hughes-Taylor scheme(HHT)是一種Multistep的數值方法[25]，及 其在求取目前的收斂位置時，不僅考慮目前的內力、外力，同時亦考慮了 上一次收斂時刻的內力、外力。在應用 HHT scheme 時，前節中 Newmark 直接積分法之數值方法及方程式(3.1.1)式、(3.1.2)式、(3.1.6)式及(3.1.7)式仍 可適用，但(3.1.3)式之不平衡力需改為

n D n

D n I n

n

n F F F P P

Ψ _1 _1(1) _1 (1) _1 (3.1.9)

其中F_n^D及P_n為t_n時刻在平衡位置的節點變形力及外力， 0為一控制參 數，當 0時，(3.1.9)式與(3.1.3)式，及 Newmark trapezoidal rule相同。

HHT scheme中之、 亦需符合(3.1.8)式之要求。

3.2 中央差分法

中央差分法在計算二階運動方程中擁有最高的精確度及穩定度，其缺 點為時間增量 需要取非常小，其原理為t t_n時刻建立運動方程式，求下個 時刻t_n1的位移U_n1，其過程如下：

當n1時在時間t_n之位移U_n及之時間t_n1之位移U_n1已知，欲求t_n1時 刻的位移U_n1，由動平衡方程式(2.9.2)，並將慣性力加速度項F_n^I1，由質量 矩陣M以及t_n時刻的加速度U _n表示可得

) ( )

( _{n n}^I2 _n

nD n

n P F U F U

U

M  _{  }  (3.2.1)

nD

F 為t_n時刻的變形力，P_n為t_n時刻的外力，F_n^I2在計算時需用到t_n時刻的 速度U_n，但在顯積分法中在t_n時刻速度為未知，故本研究以t_n1時刻速度 代替以求得F_n^I2。

令時間t_n下速度及加速度(其推導詳見(附錄G))

t t

n n

n

n n 





 



 ^  ^{ }

2 2

1 1

1 U U U

U U (3.2.2)

2 1

0 2 0

0

1 U U 2 U

U  t 

t ^





 (3.2.9) 當n1時，在t_n時刻計算U_n之數值計算過程如下：

1. 令U _n U_n1

2. 利用當前的變形位置U_n及U_n1，計算系統的節點變形力F^D(U_n)、由 當前的速度U 及變形位置求出節點慣性力速度耦合項_n F^I2(U_n)

3. 由(3.2.6)式解出位移增量U_n代回(3.2.2)、(3.2.3)可得t_n時的速度U _n、 加速度U _n

4. U_n1 U_n U_n

第四章 數值例題與結果

為了測試本研究提出的方法的有效性、準確性及效率，本章中分析了

多種不同的數值例題。本章中為了比較不同之數值計算方法的準確性及效 率，例題一、例題七、例題八及例題十，除了 Newmark trapezoidal rule 外，

還使用中央差分法(CDM)，其餘的例題皆只使用Newmark trapezoidal rule， 其中第二題為了除掉虛假的高頻振動，另外使用了 Newmark-β damping

scheme 及HHT scheme。在比較數值方法不同時，是以計算次數位移增量及

改正量次數做為探討效率的比較，Newmark法的計算次數即為總迭代次 數，CDM法的計算次數則為增量次數。本章中所有例題之容許誤差皆設為

10^4 tol 

e 。

由第二章的說明可知元素節點內力及剛度矩陣中會有趨近於零的項，

為探討這些項的影響，本研究考慮三種類型的內力及剛度，並說明如下 第一類型內力: (2.6.15)式全部考慮

第二類型內力: (2.6.15)式中去掉加底線











 的項

第三類型內力：(2.6.15)式中去掉加底線











 及



















 的項

第一類型剛度: (2.7.7)式全部考慮 第二類型剛度: (2.7.7)式中去掉加底線











 的項

第三類型剛度：(2.7.7)式中去掉加底線











 及



















 的項

為探討慣性矩陣中不同項的影響，本研究考慮三種類型的慣性矩陣，並說 明如下

第一類型：(2.7.17)式，考慮所有慣性矩陣

第二類型:(2.7.18)式，考慮所有慣性矩陣，但質量矩陣及陀螺矩陣用近似式 m

m 及cc_

第三類型:(2.7.19)式，只考慮質量矩陣且用近似式mm_

第四類型：(2.7.20)式，只考慮質量矩陣及陀螺矩陣且用近似式mm_及 c

c 

本章中若無特別說明，在 Newmark積分法時，一律使用第二類型的內力及 剛度、第四類型的慣性矩陣；在 CDM時，一律使用第二類型的內力，質量 矩陣取mT_^t m_T_。

本章中例題一、例題八、例題十為實際並不存在的梁斷面，為考量本 研究的真實性，另外討論包括矩形斷面及 I型斷面(圖五)。本章中除了例題 十一有考慮抑制翹曲，其餘皆為自由翹曲。計算形心的軸向應變皆僅考慮 元素弦長的變化。

例題一 懸臂直角梁

如圖六所示，考慮一根在X₁X₂平面上成直角的懸臂梁，在其肘處受 X3

 方向上的外力作用，其作用力大小如圖六所示，在時間t 2及移除外 力，讓它自由振動[27,28,29]。此題之幾何及材料性質為：L10，EA10⁶，

103





EI GJ

EI_{y z} ，A 1，I_y I_z 10。數值方法使用Newmark積 分法及 CDM比較。本題結果如圖七至圖十二及表一、表二所示，除另有說 明外，本例題分析時使用八個元素，Newmark積分法時，使用時間增量

2 .

0

t 。圖七為 A、B兩點位移的歷時分析，圖七可發現四個元素與八個 元素的解仍有少許的差異，圖八及圖九中分別為A、B 兩點在X₃^G方向上位 移的歷時分析與文獻[27,28,29]的結果比較，由圖八及圖九可以發現本文的 結果與文獻的結果相當吻合。

圖十Case 1為考慮(2.6.19)式之慣性力所有的項，Case 2為不考慮(2.6.19) 式之加底線的項，即不考慮速度耦合項，因本例題之軸向角速度與側向角 速度的乘積相當大，所以圖十Case 1 及 Case 2的結果在最大位移處約有 10

％的誤差。

表一為 Newmark積分法以不同類型之剛度矩陣、慣性矩陣，元素數目

去及時間增量分析時，所需的迭代次數。由表一可發現分析時用第三類型 的慣性矩陣，即只考慮近似的質量矩陣，比其他類型的慣性矩陣分析時所 需的迭代次數約多了10%。因第三類型的慣性矩陣為傳統的質量矩陣，在 推導上最簡單，計算上所需的空間及時間最少，故本研究提出之三維梁共 旋轉推導法也可只考慮質量矩陣。第四類型的慣性矩陣在計算上所需的時 間及空間比第三類型的慣性矩陣增加不多，但需要的迭代次數可以減少

%

10 ，所以本章中除另有說明外都是用第四類型的慣性矩陣。

圖十一為 CDM與Newmark積分法的位移-時間曲線。兩種方法都是使 用八個元素的結果，由圖十一可發現 CDM與Newmark積分法數值結果非 常接近，但 CDM使用的時間增量t 10^3，共310⁴個增量，Newmark積 分法使用的迭代次數約 1000次(見表一)，所以本例題CDM使用的計算次

數約為 Newmark積分法的30倍。探討質量矩陣對 CDM的最大容許時間增

量及準確性分析之影響，本例題考慮了以下兩種情況，第一種情況B1：

 

m T T

m ^t 以及第二種情況B2：mm_，表二為 CDM在不同元素數目 的最大容許時間增量。由表二發現在 CDM中，若取四個元素，使用質量矩 陣mT_^t m_T_的最大容許時間較mm_的大，但取八個元素的兩種情況 的最大容許時間增量都是t 10^3。兩種情況所得到的結果都一樣。由圖十 二可知本例題四個元素與八個元素的位移時間曲線有些微的差異，故在本

例題使用 Newmark積分法不管是效率及準確性都比CDM好。

例題二：旋轉薄圓盤的運動

本題考慮如圖十三所示之旋轉薄圓盤的運動，此圓盤藉著一無質量的

剛性桿與一球接頭連接並受重力的作用，由於此系統之運動為剛體運動，

可以應用三維剛體動力學求出此系統之初始速度及加速度見圖十四。為了 模擬此無質量剛性桿的作用，假設此剛性桿之半徑r 0.01，A10^10，

1012



EA ，GJ EI_y EI_z 2.510⁷。本例題中僅使用一個元素來進行分

析，同時將此圓盤視為元素端點上的 lump mass，本例題時間增量 取為t 10^3



t 。

在受拘束的系統中若使用 Newmark trapezoidal rule( 0)，則其加速度 會有不穩定的現象[25]，本文使用了[25]中所建議的 Newmark- damping scheme )( 0.05 及HHT scheme ( 0.05)來改善其收斂性。圖十五為 本研究 0及 0.05時，圓盤中心X 方向之加速度的歷時分析的比較，

由圖中可以發現當 0，X方向加速度會有劇烈抖動的情形，若使用HHT scheme( 0.05)或 Newmark- damping scheme( 0.05)，則抖動情形 會獲得明顯改善。本例題使用 Newmark- damping scheme( 0.05)之結 果與 HHT( 0.05)的結果相同，但在總迭代次數HHT( 0.05)卻是 Newmark- damping scheme( 0.05)的兩倍，圖十六、圖十七及圖十八

皆為使用 Newmark- damping scheme的結果。圖十六為系統角速度的變

化過程，圖十七為圓盤中心在 XY平面上的運動軌跡。圖十八為本研究圓 盤中心 X方向之加速度的歷時分析與文獻結果的比較。由圖十六-十八可知 本文結果與文獻[25]的結果相當吻合。

例題三：一雙對稱I型斷面梁承受均佈側向載重

本題考慮如圖十九所示之一雙對稱 I型斷面梁承受均佈側向載重，此梁 之幾何及材料性質：L 10 m，b0.3m，t_f 0.03m，d 0.56m，

m

t_w 0.012 ，x_C 6.9048m，楊氏係數E 210 GPa，剪力係數 GPa

G 80.77 ，線密度m188.40kg/m。表三為其斷面性質。本例題考慮 A、B兩端的邊界情況如下：(a) u_A u_B 0,w_A w_B 0,w_A w_B 0, (b) 0u_A u_B  , 0w_A w_B  , 0w_A  , (c) u_A u_B 0, 0w_A w_B  , (d) 0u_A  , 0w_A w_B  , (e) u_A 0, 0w_A w_B  , 0w_A  , (f) u_A 0,

0

 _B

A w

w , w_A w_B 0, (g) u_A 0, w_A 0, w_A 0。情況(a)-(f)負載 m

kN

P₀ 3000 / ，情況(g)負載P₀ 300kN/m。本例題中皆使用二十個元素

P₀ 3000 / ，情況(g)負載P₀ 300kN/m。本例題中皆使用二十個元素

在文檔中 雙對稱斷面薄壁梁之幾何非線性動態分析 (頁 49-0)