數值計算方法與程序

第三章 第三章

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

本文解(2.106)式的非線性平衡方程式所使用平衡迭代的數值計算方法 是採用文獻[45]中所提出基於牛頓-拉福森(Newton-Raphson)法和弧長控制 (arc length control)法的增量迭代法。本文中為了求得分歧點，以系統切線剛 度矩陣之行列式值為零來判斷。本文採用二分法決定增量位移向量的長 度，以求得系統切線剛度矩陣之行列式值為零的平衡位置。為了求得次要 平衡路徑，本文在平衡路徑的第一個挫屈負荷分歧點加入一個與第一挫屈 模態向量方向一致的擾動位移。為了本文的完整性，以下將簡單介紹文獻[45]

中提出的增量迭代數值計算方法與程序。

3.1 增量迭代法增量迭代法增量迭代法增量迭代法

若第I 個增量的平衡位置為已知，令其位移向量為Q_I、負荷參數為λ_I^， 則第I +1個增量的初始增量位移向量 Q∆ ，可以利用尤拉預測值(Euler predictor)求得[46]

rT

Q= ∆λ

∆ (3.1)

( )

^{K P}

rT = T^{I −1} (3.2)

其中∆λ為初始增量負荷參數，K_T^I 為第I 個平衡位置的系統切線剛度矩陣。

(3.1)式中的∆λ可利用下式求出[45]

( )

^{Tt T 1/2}

l r r

∆

±

=

∆λ (3.3)

其中正負符號決定方法為，當第I −1與I 個增量收歛時，其系統切線剛度 矩陣之行列式值同號，則∆λ的正負符號和第I 個增量時相同；若異號則符 號相反。∆l表示第I +1個增量的增量弧長，其值可以如下決定[45]

(

_{D I}

)

¹²

I J J

l l =∆

∆ (3.4)

其中J_D為給定的期望迭代次數，J_I 為第I 個增量迭代至平衡所使用的迭代 次數，∆l_I為第I 個增量的增量弧長。本文第一個增量的增量弧長∆l₁是由下

式決定

c max

r I l R

max 0

∆ 1 R

= (3.5)

上式中R_max為給定的參考自由度之最大位移量，R₀ 為參考外力負荷向量P 作用下的系統線性解之 Euclidean norm，I_max為給定之最大增量次數，r_c 為

R0在參考自由度的分量之絕對值。

在平衡迭代時若第I 個位置的增量位移向量 Q∆ 及增量位移負荷參數 λ

∆ 已知時，由Q_I+1 =Q_I +∆Q^，λ_{I 1+} =λ_I +∆λ，及 2.4 節、2.6 節與 2.7 節 的方法，則可以求得系統中各元素在當前的元素座標、節點位移、節點變 形角。再利用 2.9 節和 2.10 節的方法求得當前元素座標上的節點內力及剛 度矩陣，再由 2.12 節中提到的座標轉換轉到固定總體座標與節點基礎座標 的節點內力及剛度矩陣。當系統內力及外力求得後，系統不平衡力向量可 以由(2.106)式求得：

(

^Q

) (

^{F Q}

)

^P

Ψ

Ψ_I+1 = _I+1,λ_I+1 = _I+1 −λ_I+1 (3.6) 將(2.106)式在λ =λ_I+1^，Q=Q_I+1時用泰勒展開式展開為

(

1 1

)

⁰

2

1 + =

∂ + ∂

∂ + ∂

= _I+ δλ _I+ ,λ_I+

δ ^{Q Ψ}λ ^{O Q} Q

Ψ Ψ

Ψ (3.7)

其中

(

1 1

)

2

+

+ I

I ,λ Q

O 為二次以上的高次項，δ Q^{為增量位移修正量，}δλ^為負 荷參數修正量。根據(2. 106)式，(3.7)式中的∂Ψ ∂Q^與∂Ψ ∂λ^可表示為

KT

Q F Q

Ψ =

∂

= ∂

∂

∂ (3.8)

Ψ = −P

∂

∂

λ (3.9)

依牛頓法，忽略(3.7)式中二次以上的高次項，並將(3.8)式與(3.9)式代入，可 以得出

1+ − = 0

+ K r P

ΨI T δλ (3.10)

其中r=δ Q為增量位移修正量，整理(3.10)式，可得到位移修正量

(

I

)

T

T Ψ P r r

K

r= − ⁻¹ ₊₁ −δλ = ₀ +δλ (3.11)

其中r₀ = −K_T⁻¹Ψ_I+1為標準牛頓法的位移修正量，r_T在(3.2)式中已定義。

3.2 弧長控制弧長控制弧長控制弧長控制法法法 法

(3.11)式中的負荷參數修正量δλ可利用文獻[45]中所提出的定弧長控制

法決定，其方法在每一個增量中固定其增量位移向量的 Euclidean norm 為一 定值∆l^{，由新的增量位移向量}

(

∆^Q+r

)

^可以得到

(

^Q +^r

) (

^Q+^r

)

= ∆ ∆

∆l^{2 t} (3.12)

將(3.11)式代入(3.12)式可以得到

( ) (

T

)

t

l² = ∆Q+r₀ +δλrT ∆Q+r₀ +δλr

∆ (3.13)

(3.13)式經過整理後可以得到δλ^{的二次方程式}

3 0

2 2

1 +a +a =

aδλ δλ (3.14)

其中

T t

a₁ =rTr (3.15)

( )

T

a₂ = 2 ∆Q+r₀ tr (3.16)

(

0

) (

0

)

²

3 ∆ ∆ ∆l

a = Q+r ^t Q+r − (3.17)

將(3.11)式求得之增量位移修正量r以及(3.14)式求得之負荷參數修正量δλ 加入上次迭代之∆Q^與∆λ中，可以得到新的增量位移向量與增量負荷參 數，再進行下一次迭代，迭代過程將一直重複至(3.6)式中的系統不平衡力 向量滿足(2.107)式的收斂準則為止。

3.3 二分法二分法二分法二分法

利用 3.1 節的增量迭代法可以求得結構之主要平衡路徑。在每個增量的 迭代收歛時，可以得到該增量在其平衡位置的負荷參數λ^{及結構剛度矩陣的} 行列式值^D(λ)^。令λI ^及D⁽λI^{)分別表示第 I 個增量在其平衡位置的}λ^及

) (λ

D 值。λ_I+1^及D(λ_I+1^{)分別表示第 I+1 個增量在其平衡位置的}λ^及D(λ⁾

值。_∆l_I+1表示第 I+1 個增量的增量位移向量之弧長。若D(λ_I)^大於零且 )

( _I+1

D λ 小於零則可利用以下二分法求得挫屈負荷參數λ_cr^：

(1)令∆l_L =0，_∆l_{R =∆}l_I+1，λ_L =λ_I，λ_{R =}λ_I+1^{。其中下標 L 及 R 表示左界} 及右界。

(2)取

1 2

R L I

l l = ∆l +∆

∆ ₊ ，重做第 I+1 個增量迭代，並求得新的λ_I+1^及D⁽λI+1^)。

(3) 若D(λ_I+1)^{大於零，則令}λ_L =λ_I+1^，∆l_L =∆l_I+1。若D(λ_I+1)^{小於零，則令}

+1

= _I

R λ

λ ^，∆l_R =∆l_I+1。

(4) 若下列二式挫屈誤差準則同時滿足，則λ_I+1即為系統挫屈負荷，否則回 到步驟(2)，重新展開下一次二分法迭代。

^I e_D D

D ⁺ <

) 0 (

) (λ ₁

_λ

λ λ

λ e

I L R − <

+1

其中e_D及e_λ為給定的容許誤差值，本文例題之計算給定e_D =10⁻⁸及

10⁻5 λ =

e 。

3.4 數值程序數值程序數值程序數值程序

本文使用的增量迭代法之數值程序可以分為三個主要部分：

1. 輸入並計算開始分析所需的資料

(a)輸入結構資料、邊界條件與負荷參數。

(b)選擇一個參考自由度，並給定期望此自由度應達到的位移。

(c)給定最大增量次數、最大迭代次數、期望迭代次數與容許誤差值。

(d)形成系統剛度矩陣並求得(3.5)式中的R₀。

(e)利用(3.1)式、(3.3)式與(3.5)式計算初始增量位移向量、初始增量負荷

參數與第一次增量弧長。

2. 使用迭代法求增量的收歛解

(a)由 2.4 節的方法，利用已知的增量位移求得當前元素的變形向量q_θ和 ub，並計算(2.98)式中元素節點內力f_E，接著由 2.12 節座標轉換將f_E利 用(2.101)式轉換至f_B，然後將元素節點內力f_B組合成結構系統節點內力

F。

(b)計算(2.106)式的不平衡力_Ψ。

(c)檢查(2.107)式的收斂準則，若滿足則進入第 3 部分；否則檢查迭代數，

如果小於給定之最大迭代次數，則進行步驟(d)；否則減少增量弧長並以 (3.1)式與(3.3)式計算新的增量位移向量與增量負荷參數，回到步驟(a)重 新計算。

(d)利用(3.11)式與(3.14)式計算增量位移修正量與增量負荷參數修正量，

回到步驟(a)重新計算。

3. 計算下一次增量所需的資料

(a)檢查參考自由度的位移及進行的增量次數是否已達給定值，若已達到 給定值則停止分析工作；否則進行步驟(b)。

(b)計算(3.11)式中的切線剛度矩陣K_T。本文中K_T的計算方法是利用當

前變形位置的元素座標重新計算元素剛度矩陣k，再加上元素幾何剛度 矩陣k_σ，本文中忽略元素其他的幾何變形與內部應力對元素剛度矩陣造 成的影響，再利用 2.12 節座標轉換將k_E利用(2.105)式轉換至k_B，然後 將元素剛度矩陣k_B疊加成結構系統剛度矩陣K_T。

(c)利用(3.1)式、(3.3)式與(3.4)式計算下一次增量的增量位移向量、增量 負荷參數與增量弧長。

(d)回到第 2 部分進行迭代工作。

在文檔中 平面三角形殼元素之改善研究 (頁 41-46)