第三章 數值計算方法與程序
3.1 靜態分析
本文解非線性靜平衡方程式(2.157)式(不包含系統顯節點慣性力F )I 的數值計算方法是基於牛頓－拉福森(Newton-Raphson)法配合弧長控制(arc length control)法的增量迭代法[81]。為了求得挫屈負荷，本文採用文獻[3]
中所提出的二分法，決定增量位移向量的長度，以求得系統切線剛度矩陣 之行列式值為零的平衡位置。為了改善平衡迭代的收斂情況，本文中採用 了N循環迭代法[56]。以下對本文使用的增量迭代法、二分法及N循環迭 代法加以說明。
3.1.1 增量迭代法
若第I 個增量的平衡位置為已知，則在此位置的系統切線剛度矩陣 KT可以求得，且第 I+1個增量的初始增量位移向量 Q∆ ，可利用尤拉預測 值(Euler predictor)求得
T, Q Q=∆λ
∆ (3.1)
其中∆λ為初始增量負荷參數，QT =KT−1Pref 為參考負荷向量Pref 的切線 解。∆λ可利用下式求出[81]
, ) (QtTQT 1/2
∆l
±
=
∆λ (3.2)
其中正負符號之決定方法如下：若第 I 與 I-1個增量收斂時，系統切線剛度 矩陣之行列式值同號，則∆λ的正負符號和第 I 個增量時相同；若異號則符 號相反。 l∆ 表第 I+1個增量的增量弧長，其值可以如下決定
, )
(JD JI 1/2 lI
l= ∆
∆ (3.3)
其中JD為給定的期望迭代次數，JI為第I 個增量，迭代至平衡所使用的 迭代次數，∆lI為第I個增量的增量弧長。
本文中第一個增量的增量弧長∆l1是由下式決定
rc
I R
max 0 max 1
= R
∆l (3.4)
上式中Rmax為給定的參考自由度之最大位移量， R0 為參考負荷向量Pref
作用下的系統線性解R0的 Euclidean norm，Imax為給定之最大增量次數，
rc 為R0在參考自由度的分量的絕對值。
在平衡迭代時，若增量位移向量 Q∆ 及增量負荷∆λ已知，則利用2.6 節的方法，可求得梁結構新的變形位置，再利用2.9與2.10節的方法，可求 得 元 素 座 標 上 的 節 點 內 力 及 剛 度 矩 陣 。 對 應 此 位 置 的 負 荷 參 數 為
λ λ
λ = I +∆ ，其中λI 為第I個增量達平衡時的負荷參數，∆λ即增量負荷 參數。當系統內力及外力求得後，不平衡力量 ψ 向量可由(2.157)式求得。
若(2.158)式的收斂準則不能滿足，則利用定弧長控制法[81]，求得一位移 修正量 Qδ 與負荷參數修正量δλ，並加入前一次迭代的 Q∆ 與∆λ中，而得 一新的增量位移向量與增量負荷參數，再進行下一次的迭代。 Qδ 與δλ可 由下列二式決定
)
1( P
K
Q δλ
δ = T− −ψ+ (3.5)
) (
)
2 =(∆Q+δQ ∆Q +δQ
∆l t (3.6)
以上之迭代計算過程一直重覆至滿足(2.158)式的收斂準則為止。
3.1.2 二分法
利用3.1節的增量迭代法可以求得結構之主要平衡路徑。在每個增量的 迭代收斂時，可以得到該增量在其平衡位置的負荷參數λ及結構切線剛度 矩陣的行列式值D(λ)。令λI 及D(λI)分別表示第I個增量在其平衡位置的 λ及D(λ)之值。λI+1及D(λI+1)分別表示第I +1的增量在其平衡位置的λ及
) (λ
D 之值。∆lI+1表示第I +1個增量的增量位移向量之弧長。若D(λI)大 於零且D(λI+1)小於零則可利用以下二分法求得挫屈負荷參數λNB：
(1) 令
∆
lL =0,∆
lR =∆lI+1,λL =λI ,λR =λI+1，其中下標L及R表示左界 及右界。(2) 取
1 L 2 R
I l l
l = ∆ +∆
∆ + ，重作第I +1個增量迭代，並求得新的λI+1及 )
( I+1
D λ 。
(3) 若D(λI+1)大於零，則令λL =λI+1,∆lL =∆lI+1 若D(λI+1)小於零，則令λR =λI+1,∆lR =∆lI+1 (4) 若下列二式挫屈誤差準則同時滿足
D
I e
D
D + <
) 0 (
) (λ 1
(3.7)
λ λ
λ
λ e
I L R − <
+1 (3.8)
其中eD及eλ為給定的容許誤差值
則λI+1即為系統挫屈負荷，否則回到步驟(2)，重新展開下一次二分法迭 代。
經由二分法求得挫屈負荷λ，再利用系統切線剛度K(λ)計算挫屈模 態，以下將說明挫屈模態的計算程序
(1) 將K(λ)分解成一下三角矩陣L及上三角矩陣U使K =LU，其中L的對