本文解非線性平衡方程式(2.54)式的數值計算方法是採用文獻[25]
中所提出基於牛頓-拉福森(Newton-Raphson)法和弧長控制(arc length control)法的增量疊代法。為了文章的完整性,本章中將簡單介紹文獻 [25]中提出的數值計算的方法與程序。
3.1 增量迭代法
在增量迭代法中,若第 I 個增量的平衡位置為已知,令其位移向量 為QI、負荷參數為λI,則第 I+1 個增量的初始增量位移向量∆Q,可 利用尤拉預測值( Euler predictor )求得[34]
∆Q = ∆λrT , (3.1)
G
的平衡迭代:由QI 1+ =QI +∆Q,λI 1+ =λI +∆λ,參考位移負荷向量 及 2.4 節與 2.6 節的方法,則可以求得系統中各殼元素在當前的元 素座標、節點變形位移、節點變形角。再利用 2.5 節的方法,求得在 當前元素座標上的節點內力及剛度矩陣再由標準轉換轉到總體座標的 節點內力及剛度矩陣。由(2.54)式此時之不平衡力可表示成
QP
Iٛ+1 =F
(
QI+1,λI+1QP)
, (3.11) 將(2.54)式在λ =λI+1,Q=QI+1時用泰勒展開式展開為= + + + δλ+ λ
∂ δ ∂
∂
∂ F
Q Q Ψ F
Ψ ٛ I 1 (二次以上高次項)
=ΨI+1 +KTr+δλFP +(二次以上高次項) (3.12)
其中, Q
K F
∂
= ∂
T 為系統切線剛度矩陣,r =δQ為增量位移修正量,
在(3.3)式中已定義,在此我們必須注意 在每一次迭代後都要重新計 算。
FP
FP
依牛頓法,忽略(3.12)式中二次以上高次項,可得到位移修正量 r = −KT−1(ΨI+1 +δλFP)=r0 +δλrT , (3.13) 其中,r0 =−KT−1ٛΨI+1,rT在( 3.2 )式中已定義。
(3.13)式中的位移負荷參數修正量δλ有很多不同的決定方法[37],如 Crisfield 在文獻[35]所提出的定弧長控制法及 Fried 在文獻[38]提出的 正交軌線法,本文是使用 Crisfield[35]提出的定弧長控制法決定δλ, 並將在 3.2 節中介紹。
將求得之位移負荷參數修正量δλ及增量位移修正量r加入上次迭代 之∆λ與∆Q中,可得到一新的增量位移向量與增量位移負荷參數,再 進行下一次迭代。以上之迭代過程一直重複至(3.11)式中的不平衡力滿 足(2.55)式的收斂準則為止。
3.2 弧長控制法
(3.13)式中的位移負荷參數修正量δλ可利用文獻[35]中 Crisfield 所 提出的定弧長控制法決定,其方法使在每一個增量中固定其增量位移 的 Euclidean norm 為一定值 l∆ 。由新的增量位移向量為
(
∆Q+r)
可得 ∆l2 =(
∆Q+r0 +δλrT) (
t ∆Q+r0 +δλrT)
, (3.14) 上式經過整理後可以得到δλ的二次方程式a1δλ2 +a2δλ+a3 =0, (3.15) 其中
a1 =rTtrT,
a2 =2
(
∆Q+r0)
trT ,a3 =
(
∆Q+r0) (
t ∆Q+r0)
−∆l2,當(3.15)式的解為兩實根δλ1,δλ2時,本文中取兩根中能使新的增 量位移
(
∆Q+r)
與前次迭代增量位移∆Q間的內積(inner product)較大 者為δλ。當(3.15)式無實根時,則把(3.9)式之 l∆ 減小,重做第 I+1 個 增量。3.3 數值程序
本文使用的增量迭代法之數值程序可以分為三個主要部分 1. 輸入與計算開始分析所需的資料
(a) 輸入結構與負荷資料。
(b) 選擇一個參考自由度,並給定期望此自由度達到的最大位移。
(c) 給定最大增量數、每個增量期望的迭代數與最大的迭代數、收 斂時的容許誤差。
(d) 形成剛度矩陣並求得(3.10)式中的R0。
(e) 利用(3.1)-(3.3)式、(3.9)式與(3.10)式計算第一次增量的增量弧 長、增量負荷參數與初始增量位移。
2. 使用迭代法求增量的收斂解
(a) 利 用 已 知 的 增 量 位 移 求 得 目 前 元 素 的 變 形 , 並 以 (2.39) 、 (2.45)式計算元素節點內力。
(b) 計算(2.54)式的不平衡力Ψ。
(c) 檢查(2.55)式的收斂準則,若滿足,則進行(e);否則檢查迭代數,
如果小於給定之最大迭代數,則進行(d);否則減少增量弧長並以 (3.1)-(3.3)式與(3.9)式計算新的增量位移與增量負荷,回到步驟(a)重新 此增量。
(d) 以(3.13)式與 3.1 節的方法計算增量位移修正量與增量負荷參數 修正量,然後回到步驟(a)。
3. 計算下一次增量所需的資料
(a) 檢查參考自由度的位移或已進行的增量次數是否已達給定值。
若已達到給定值則停止分析工作;否則進行下一步工作。
(b) 計算(3.12)式中的切線剛度矩陣與FP。
(c) 以(3.1)-(3.3)式與(3.8)-(3.9)式計算下一次增量的增量位移、增量 弧長與增量負荷參數。
(d) 回到第二部份執行迭代工作。