元素節點慣性力向量之推導

+ + +

+

= ^{t t D} _R^{t t D}

D Tθφ Tθφ fθ T TϕφTθφ fθ

f ( ^{0 1} ) ( ^{1 0} ) (2.112)

其中f_θ^D+是將f_θ^D在(2.108)及(2.109)式中 f₁₂^θG_iL (i=b,c)改成 f₁₂G_iL。

2.9.2 元素節點慣性力向量之推導

為推導方便，本文將(2.87)式之δq 拆成四個廣義位移向量ϕ

δuϕ_i (i=a,b,c,d)

δu^ϕ_a ={δu₁,δu₂}， δu^ϕ_b ={δv₁,δϕ₃₁,δv₂,δϕ₃₂}

δu^ϕ_c ={δw₁,δϕ₂₁,δw₂,δϕ₂₂} ，δu^ϕ_d ={δϕ₁₁,δβ₁,δϕ₁₂,δβ₂} (2.113) 對應於δu^ϕ_i ⁽i=a^,b^,c^,d⁾的廣義節點慣性力向量f_i^{I (}i=a^,b^,c^,d⁾可表示如 下

} , { ₁₁^I ₁₂^I

I

a = f f

f ， f_b^I ={f₂₁^I ,m₃₁^I , f₂₂^I ,m₃₂^I }

}

利用(2.113)及(2.114)式，(2.93)式可改寫成

∫

]

( && && && && && && & &

&& N u^θ N u^θ θ θ

dx

+ρAy²_p

∫

N_d(α_x +θ&&₁−Ω_yv&_,x +Ω_zw&_,x +2Ω_yN^t_eu&^θ_b)dx

∫

⁺

∫

^{+ −}

∫

⁺

dx

+

} AEy AEy

L

AEz AEy AEy

∫

AEz AEy

L

∫

ϕ && && && && && &&

&& u ϕ u ϕ

q = (2.139)

其中f 定義在(2.90)式，可以由(2.129)式之^I f_i^I (i=a,b,c,d)組合而成，

u&& 、j ϕ&&_j (j=1,2)分別為節點 j 的絕對加速度及絕對角加速度。 m 為一對稱 矩陣，利用直接勁度法[76]，m 可由下列的次矩陣組合而成

_ϕ

j I i

ij u

m f

&&

∂

= ∂ (2.140)

其 中 i, j =a,b,c,d ， f_i^I 為 (2.129) 式 之 節 點 慣 性 力 ， u&& =^ϕ_a {u&&₁,u&&₂} ， }

, , ,

{ ₁ ϕ_{31 2} ϕ₃₂

ϕ && && && &&

&&_b = v v

u ，u&&^ϕ_c ={w&&₁,ϕ&&₂₁,w&&₂,ϕ_&&22} ， u&& =^ϕ_d {ϕ_&&11,β^&&₁,ϕ_&&12,β^&&₂}。 將 (2.23)－(2.27)、(2.55)、(2.58)及(2.129)式代入(2.139)式可得

dx A _a ^t_a

aa ⁼

∫

^{N N}

m ρ

∫

= Ay_{p a} ^t_edx

ab N N

m ρ

∫

= Az_{p a} ^t_fdx

ac N N

m ρ 0 m_ad =

∫

⁺

∫

⁺

∫

^{′ ′}

= Ay_{p e} ^t_edx A _b ^t_bdx I_{z b b}^tdx

bb N N N N N N

m ρ ² ρ ρ

∫

= Ay_pz_{p e} ^t_fdx

bc N N

m ρ

∫

= Az_{p b} ^t_ddx

bd N N

m ρ

∫

⁺

∫

⁺

∫

^{′ ′}

= Az_{p f} ^t_fdx A _c ^t_cdx I_{y c c}^tdx

cc N N N N N N

m ρ ² ρ ρ

∫

−

= Ay_{p c} ^t_ddx

cd N N

m ρ

∫ ∫

∫

^{+ ′ ′ + +}

= I_{p d} ^t_ddx I _{d d}^tdx A y_p z_{p d} ^t_ddx

dd N N N N N N

m ρ ρ _ω ρ ( ^{2 2})

(2.141)

因為系統平衡方程式是定義於總體座標系統中，所以(2.131)式的元 素顯切線剛度矩陣及(2.138)式質量矩陣需經下列的座標轉換，方可組合成 系統矩陣。

GE,

GE G T kTt

k =

GE,

GE

G T mTt

m =

(2.142)

⎥⎥

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎣

⎡

=

I2

0 0

0 0

0 A

0 0

0

0 0 A

0 0

0 0 0

A 0

0 0 0

0 A

T

t t

t t

GE GE

GE GE

GE (2.143)

其中0是一3× 3階的零矩陣，0是一3× 2階的零矩陣。及A_{GE為(2.1)式中總} 體座標與元素座標間之轉換矩陣。

2.11 與變形位置相關之節點作用力與負荷剛度矩陣

本文中僅考慮由固定方向之保守力造成的與變形位置相關之節點力 矩。由保守力造成的節點力矩是保守力矩[77, 78]；但是其大小與方向通常 會隨著變形而改變[77,79]。與變形位置相關的節點作用力，在結構變形時 對系統剛度矩陣的貢獻稱為負荷剛度矩陣(load stiffness matrix)。本文中與 變形位置相關的節點作用力(configuration dependent load)、負荷向量(load vector)、以及負荷剛度矩陣都是在一個負荷基底座標X_i^P( i = 1, 2, 3)上定 義。如2.2節所述，X_i^P的原點是剛接在與變形位置相關的外力作用的節點 上。在本文中上標 P 表示該量是在負荷基底座標上定義。在本節之推導，

除有特別之聲明以外，各向量都是表示成負荷基底座標的分量。為了表達 之簡潔，在不會造成混淆時，本節中有時省略了向量之上標 P 。本文中定 義了兩種與變形位置相關的節點作用力機制，並分別說明如下：

(1)第一型節點作用力機制。如圖2.7所示，一個半徑 R 的剛性圓盤，

其圓心 O 剛接於節點上，在盤緣上纏繞著二條繩子，二條繩子的自由端承 受了一對大小相等方向相反的保守力。假設在梁變形前，繩子和圓盤位於 同一平面上；梁變形後，圓盤和梁的節點一起平移和旋轉，但繩子的方向 仍維持不變，所以在梁變形後，繩子和圓盤通常不在同一平面上。 A 點為

繩子與圓盤邊緣的接觸點，若假設剛性圓盤是半徑為 R 球心為 O 之球上的 一個大圓，則OA 和繩子互相垂直。

為了描述上的方便，如圖2.7所示，本文將X₁^P軸定在圓盤的法線方 向，X₂^P軸定在盤面上的任意方向，X₃^P軸則由右手定則決定。

在圖2.7中，因 OA 和繩子及X₁^P軸垂直，所以 OA 線上的單位向量e_A 可表示成

2 /

)1

(a a a

e_A = ^t (2.144)

} , , 0

{ _{3 2}

1 = l −l

×

=e^P_P e^P

a (2.145)

其 中 e^P_P ={l₁,l₂,l₃}為 繩 子 方 向 ( 或 保 守 力 P 方 向 ) 上 的 單 位 向 量 ， 0}

0, , 1

1^P ={

e 為X₁^P軸上的單位向量。雖然繩子在固定座標中的方向是固定 的，但因負荷基底座標隨梁變形改變而改變，所以e^P_P的分量l_i( i = 1, 2, 3) 亦隨變形改變而改變。

圖2.7中之二平行力對O 點造成的力矩M可以表示成

P P

RPeA e

M= 2 × (2.146)

當 與 圓 盤 剛 接 的 節 點 受 到 一 個 微 小 的 旋 轉 向 量 ∆ϕ^P = }

, ,

{∆ϕ₁ ∆ϕ₂ ∆ϕ₃ 擾動時，負荷基底座標會有一微小的旋轉，所以e₁^P會有 一微小的變量∆e₁^P。當∆ϕ_i(i=1,2,3)足夠小時，∆e₁^P可以表示成

} ,

, 0

{ _{3 2}

1

1 =∆ × = ∆ϕ −∆ϕ

∆e^P ϕ^P e^P (2.147)

因(2.146)式中之e_A為e₁^P的函數，所以(2.146)式中的M亦有一變量∆M^如下 式所示

)

⎥⎥ 力矩、第二型準切線(QT2, quasitangential of the second type)力矩、及半準 切線(ST, semitangential)力矩。從圖2.8中，可以發現QT1力矩中外力 P 和

由於系統平衡方程式是定義在固定總體座標系統中，負荷剛度矩陣 KP需先經以下的座標轉換，方可疊加到系統剛度矩陣中，亦即

t GP P GP G

P A K A

K = (2.156)

其中A_GP為(2.1)式中，總體座標與負荷基底座標間的轉換矩陣。

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

系統離散化後的動平衡方程式，是由在固定元素座標上算出的元素節 點變形力、慣性力經標準的座標轉換後，在總體座標上所合成之系統內力 向量及系統外力向量所組成，可表示如下

0 P F F

Ψ = ^D + ^I − = (2.157)

其中Ψ為不平衡力向量，F^D為系統變形力內力向量，F^I是系統慣性力內 力向量，P表系統外力向量。F^D、F^I可分別由(2.112)式之元素節點變形力

向量及(2.90)式之元素節點慣性力向量經由元素座標轉換到固定總體座標

上組合而成，P可由與變形位置相關的節點作用力及與變形無關之節點作 用力組合而成。

本文以不平衡力向量Ψ 的 weighted Euclidean norm 作為迭代時的誤 差度量，而且收斂準則表示為

etol

e= N ≤ λ*

ψ (2.158)

其中N表離散系統的自由度數，在靜態分析時λ^{*取為 P ，動態分析} 時λ^*取為1，etol是一給定的容許誤差值。

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

本章將在3.1節及3.2節分別說明開口薄壁梁之幾何非線性靜態分析及 動態分析的數值計算方法與程序。本研究分析薄壁梁結構的自然振動及其 在負載下之靜態平衡位置的自然振動是採用次空間法(subspace iteration method) [75]解廣義的標準的特徵值問題。

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)

其中∆λ為初始增量負荷參數，Q_T =K_T⁻¹P_ref 為參考負荷向量P_{ref 的切線} 解。∆λ可利用下式求出[81]

, ) (Q^t_TQ_T ^1/2

∆l

±

=

∆λ ^(3.2)

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

, )

(J_D J_I ^1/2 _lI

l= ∆

∆ (3.3)

其中J_D為給定的期望迭代次數，J_I為第I 個增量，迭代至平衡所使用的 迭代次數，∆l_I為第I個增量的增量弧長。

本文中第一個增量的增量弧長∆l₁是由下式決定

rc

I R

max 0 max 1

= R

∆_l (3.4)

上式中R_max為給定的參考自由度之最大位移量， R₀ 為參考負荷向量Pref

作用下的系統線性解R₀的 Euclidean norm，I_max為給定之最大增量次數，

rc 為R₀在參考自由度的分量的絕對值。

在平衡迭代時，若增量位移向量 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,∆l_L =∆l_I+1 若D(λ_I+1)小於零，則令λ_R =λ_I+1,∆l_R =∆l_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的對

3.2 動態分析

本文解非線性動平衡方程式(2.157)的數值計算方法是採用基於牛頓-拉福 森(Newton-Raphson)法配合 Newmark 直接積分法[84]的增量迭代法 。本文 中使用 Newmark trapezoidal rule為數值積分方法，然而此法應用在受拘束 的系統時，會有不穩定的現象發生，若加入 Numerical damping 則可獲得 改善。近年來，有許多具有 Numerical damping 的數值方法被提出，如 Nemark-β damping scheme及 Hilber-Hughes-Taylor scheme [85]。故本文除 了 Newmark trapezoidal rule外，亦探討了上述兩種具 damping 效果的數值 方法的準確性與收斂性。

3.2.1 Newmark直接積分法

Newmark 直接積分法乃假設在時間t_n時，滿足動平衡方程式(2.157)的平 衡位置、速度U& 和加速度_n U&& 為已知，則在_n t_n+1 =t_n +∆t時刻之平衡位 置、速度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.11)

U&⁰_n+1=U& _n+∆t[(1−γ)U&&_n+γU&&_n+1] (3.12) )2

1 4(

1 α

β = − 、 (1 2 ) 2

1 α

γ = − (3.13)

其中α ≤0。

2. 令U&&_n+1 =U&&⁰_n+1，U&_n+1 =U&⁰_n+1

3. 由∆ (U 或δU)及上次迭代後的變形位置得到此迭代的變形位置，再算出 系統的節點變形力F_n^D₊₁(詳見附錄 F)。由U&_n+1、U&& _n+1及最近的變形位置 算出系統的節點慣性力F_n^I₊₁(詳見附錄 F)。由在t_n+1時刻的外力P_n+1、

D n+1

F 、F_{n 1}^I₊ 及系統平衡方程式(2.157)算出系統不平衡力

Ψ_n+1 =F_n^I₊₁ +F_n^D₊₁ −P_n+1 (3.14)

4. 檢查 Ψ_n+1 是否滿足收斂準則，若滿足則迭代停止。

5. 若不滿足，則由

¹ ] 1

[ˆ ⁻ ₊

−

= K Ψ_n

δU (3.15)

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

1 [ ] [ ] ˆ]

[K 2 M + K

= ∆

β t (3.16)

為系統有效剛度矩陣(effective stiffness matrix) ，[M 、] [K ] 表系統之質量與剛度矩陣。

6. 令

U U δU β ²

1 0

1

1

n t

n₊ = && ₊ + ∆

&& (3.17)

U U δU β

γ

n t

n₊ = ₊₁ + ∆

0

1 &

& (3.18)

7. 回到 2.

若上述迭代程序經過一預先設定的迭代次數還不收斂時，則將時間增量

∆t減半，令

1 2 t t

t_n+ = _n + ∆ ，然後再重複上述迭代程序。但在下一時刻，則

將時間增量改為

2 t+ ∆t

∆ 。

Newmark 直接積分法中，(3.13)式中若取α =0，即β =1 4、γ =1 2 則 稱為 Newmark trapezoidal rule [84]，本文的數值例題中除非另有說明都採 用 Newmark trapezoidal rule。因 Newmark trapezoidal rule 不會造成振幅衰 減，無法除掉高頻的振動，為了使 Newmark 直接積分法能造成振幅衰

Newmark 直接積分法中，(3.13)式中若取α =0，即β =1 4、γ =1 2 則 稱為 Newmark trapezoidal rule [84]，本文的數值例題中除非另有說明都採 用 Newmark trapezoidal rule。因 Newmark trapezoidal rule 不會造成振幅衰 減，無法除掉高頻的振動，為了使 Newmark 直接積分法能造成振幅衰

在文檔中 薄壁開口梁之自由振動分析及幾何非線性動態反應研究 (頁 60-0)