系統平衡方程式

本研究考慮的負荷有力負荷及位移負荷兩種，所以非線性系統的運動

方程式可表示成

Ψ(Q,λ)F(Q,λ_dQ_d)_fP0 (2.7.1) 其中Ψ 為系統的不平衡力，F 為系統節點內力，Q 為系統的節點位移，_d為 位移負荷參數，Q_d為參考位移負荷向量，_f 為力負荷參數，P為參考力負 荷向量，其中F由系統節點變形力F^D、系統節點慣性力F^I組成。F^D及F^I

可由(2.4.40)–(2.4.44)式中的元素節點變形力、元素節點慣性力從當前的元素

座標轉換到總體座標後組合而成。本研究在每一增量平衡迭代的過程中，

將_d或是_f 的值固定，(2.7.1)式中之在迭代中代表_d或是_f 中可變動的 負荷參數。

若不考慮慣性力，則(2.7.1)式為非線性平衡方程式。若將(2.7.1)式在靜 態平衡的位置，用泰勒級數展開到一次項，則可得到系統在靜態平衡點的 線性振動方程式

0 Q K Q

M    (2.7.2) 其中M為系統的質量矩陣，Q _{為系統的節點加速度，}K為系統的切線剛度 矩陣，Q為靜態平衡點算起的系統節點位移。M、K 是由元素的質量矩陣、

切線剛度矩陣從當前的元素座標轉換到總體座標後組合而成。

若振動方程式(2.7.2)式存在自然振動頻率，則其解的形式可以表示如下

t

ei^

Φ

Q  (2.7.3) 其中Φ是振態，是自然振動頻率。將(2.7.3)式代入(2.7.2)式可得

MΦ

KΦ² (2.7.4)

(2.7.4)式為一廣義特徵值問題，振動模態Φ 是特徵向量，特徵值是²。

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

本章將說明挫屈梁在不同位移負荷及/或力負荷下的幾何非線性分析 及線性自由振動分析所需要用到的數值計算方法。本文的挫屈梁是指如圖 1-1 所示之直梁，受軸向壓縮 引起側向挫屈的梁，本文以梁中點拱起的高 度 作為挫屈梁挫屈的程度。壓縮量 與梁中點拱起的高度 有一非線性關 係，本文中視壓縮位移 為位移負荷求得中點拱起高度為 的挫屈梁，並探 討挫屈梁的的的自然振動頻率及振態，接著在挫屈梁的中點施加側向力負 荷 P，探討挫屈梁的幾何非線性行為，挫屈梁受側力的平衡路徑上可能會有 極限點、分歧點，在極限點有不穩定平衡路徑，會產生跳躍現象，在分歧 點有穩定或不穩定的次要平衡路徑。本文採用牛頓—拉福森（Newton - Raphson）法配合弧長控制（arc - length control）法[9]的增量迭代法，來求 得完整的平衡路徑。

在平衡路徑上的臨界點(極限點、分歧點)，系統切線剛度之行列式值 為零，可用來決定挫屈梁之臨界點的位置。因本文要追蹤不同（或 ）之 挫屈梁的臨界點，即臨界點的摺線[12]，如用上述方法，需要大量的計算，

為節省計算量，本文中使用 Eriksson[12-14]中延伸系統的觀念，即將挫屈時 切線剛度矩陣絕對值最小的特徵值為零，當作約束條件，得到一平衡方程 式的延伸系統，可以直接求得平衡路徑上的臨界點。在追蹤不同 之挫屈 梁的臨界點時，為節省計算量，不再重新計算其平衡路徑，而直接由目前 之臨界點對應的力/位移負荷及平衡位置，將一增量力/位移負荷加到原來 的臨界力/位移負荷，當作對應於某一未知 /P_CR之挫屈梁的臨界力/位移 負荷，將 /P當作位移/力負荷參數，使用正交軌線法[16]得到一個初始增 量負荷參數與初始增量位移，然後再回到延伸系統求得新的臨界力/位移負

荷 對 應 的 / P_CR 。 延 伸 系 統 中 附 加 的 條 件 可 以 根 據 計 算 的 目 的 更 動 [12-14]，用於弧長控制法與計算特徵值與特徵向量。

本文使用子空間法[2]，求出廣義特徵值問題(2.7.4)式的特徵值²及對 應的特徵向量Φ 。在直接計算挫屈負荷時，若將質量矩陣取為單位矩陣，

則(2.7.4)式可以寫成KΦ g IΦ，該標準特徵值問題，在平衡路徑上的臨界 點，有一特徵值為零，與其對應之特徵向量是挫屈模態。 

本章將在 3.1 節中詳述本文需要的各種數值方法，在 3.2 節說明整個數 值程序迭代的流程。

3.1 數值方法

本研究考慮的負荷有力負荷及位移負荷兩種，在每一增量平衡迭代的 過程中，將(2.7.1)式之平衡方程式中之_d或是_f 的值固定。

3.1.2 增量迭代法

增量迭代法通常用來分析結構的幾何非線性分析，若第I 個增量的平衡 位置為已知，其位移向量為Q_I，負荷參數為_I，則可以求得這個平衡位置 的系統切線剛度K_T ，且第 I+1 個增量的初始增量位移可以用尤拉預測值

（Euler predictor）求得 rT

Q

 (3.1.2) 其中r_T是第I 次增量的切線解，在力負荷時，

P r

K_{T T}  (3.1.3) 在位移負荷時，

P T

Tr F

K  (3.1.3) 其中 

F

F_P ，是第 I+1 增量初始增量負荷參數，本文配合弧長控制法

[9]，選用

t T Tr r







 (3.1.5)

其中 Q 是初始增量位移的弧長， 之正負號則與 r_T^tQ_I的正負號一 致，Q_I^為上一增量中收斂的增量位移。

迭代過程中對增量的修正，即位移改正量和負荷改正量，詳見 3.1.6 節 之弧長控制法。

3.1.3 收斂準則 

本文以系統的不平衡力的 weighted Euclidean norm 作為迭代時的誤差 度量，當負荷為力負荷時，收斂準則表示為

e_tol e N 

P Ψ

 (3.1.6) 當負荷為位移負荷時，收斂準則表示為

_tol

R

e e 

F

Ψ (3.1.7)

其中N表離散系統的自由度數，F_R為系統節點反力向量，e_tol為一給定之容 許誤差值。

3.1.4 延伸系統（Extended system）

考慮一個靜態力學平衡的離散化模型，可定義一個延伸系統的控制方 程式[12-14]

















0 0 Z

G Z Ψ

) (

)

( (3.1.8)

其中Ψ是N個平衡方程式的集合，G 是M個輔助條件的集合，這些條件在

3.1.5 延伸系統延遲分解法（Postponed factorization of an extended system） 如果一個矩陣可以表示成延伸系統的模式：

3.1.6 弧長控制法（Arc-length method） (3.1.18)式中，L 通常可以忽略不計，所以(3.1.18)式之負荷參數改正量可以 近似成

3.1.7 特徵值直接計算法（Direct computation of eigenvalues）

標準特徵值問題可以表示成

到g  g^(i1)、ΦΦ^(i1)。

3.1.8 臨界值直接計算法（Direct computation of critical points）

平衡路徑的臨界點通常指極限點或分歧點，在臨界點系統切線剛度矩

Qr₀ r_T (3.1.34)

r0是(3.1.19)式中標準牛頓法中，不平衡力的位移改正量。

3.1.9 正交軌線法（Orthogonal trajectory accession）

增量迭代過程中的負荷參數修正量，也可以使用文獻[16]所提的正 交軌線法得到，其方法是在由N 維的位移向量與負荷參數所組成的N 1維 空間中，令

t

 

r_T,1^t (3.1.35) n^t 



Q,



(3.1.36) n^tt0 (3.1.37) 其中n 為包含負荷參數改正量 Q^t  與增量位移改正量的N 1維修正量，

t 是將參考負荷 P 或F 的切線解向量_P r 放大一維，將_T (3.1.34)式帶入(3.1.37) 式整理，可以得到負荷參數改正量

t T T t T

r r

r r

 

 1

 0 (3.1.38)

將再帶回(3.1.34)式則可以得到增量位移改正量Qr₀ r_T。

3.1.10 免特徵值分析法（Eigenanalys-free method）

文獻[17,29]中提出一個的方法，可以不需要進行特徵值的計算，直接 在剛度矩陣的LDL^t分解（LDL^t-decomposed stiffness matrix）運算中得到 一個最接近 0 的近似特徵值，當特徵值靠近 0



g 0.0,g 0.0



時，它對應 的近似特徵向量

 

t j

j L e

S  ^1 (3.1.39) 其中 j 是對角矩陣D中絕對值最小的位置，e 是第_j j 維純量等於1，其他部

分為 0的N 維單位向量。

這個方法同時可以得到近似特徵值

2 j

Dj

g S

 (3.1.40)

其中D_j是對角矩陣D 中第 j 個分量，當我們的平衡位置越靠近臨界負荷，

即特徵值越靠近 0，所求得的近似特徵值與特徵向量也越靠近真正的值。

這個方法所花費的計算時間非常少，因為剛度矩陣的分解K LDL^t在 進行平衡路徑的迭代的過程中會一直被重複計算，L^t矩陣有上三角結構特 性，所以解L^tS_j  得到e_j S 。 _j

3.1.11 子空間法（Subspace method）

分析結構振動的自然振動頻率，要求解方程組KΦ²MΦ，其中K 是系統剛度矩陣，M 是系統質量矩陣，可以把他當作廣義特徵值問題

BΦ

AΦ g (3.1.41) 來求解，本文採用子空間法 [2]，求出特徵值g及對應的特徵向量Φ ，而在 程式中當我們需要計算標準特徵值（standard eigenvalue problem）問題

Φ

AΦ g 時，令B 單位矩陣（II dentity matrix）。

子空間法以將N維的大型特徵值問題近似的轉化為一個P維向量空間 的小型特徵值問題，從而使求解的問題縮小，節省計算時間，但是當遇到 特徵值為負的，無法求出答案，所以本文在計算之前首先利用 3.1.10 的免 特徵值分析法[17,29]，這個方法會提供一個切線剛度的近似特徵值(3.1.40) 式，雖然這個近似值在它距離 0 越遠時越不準確，但是它的正負號與真正 的特徵值解是相同的，所以一旦發現它是負號時，在運算中給一常數c，讓 方程組(3.1.41)成為

其中g' gc0，如此可以順利迭代求出特徵值 g。

1

STOP

1

 _C

C N

N

YES

NO 數值程序流程圖：

 

輸入初始值 結構幾何、材料性質 負載條件、求解類型

增量數、最大迭代數及收斂時的容許誤差

階段數NC_max、每一階段增量數I_max

計算線性剛度、線性解

0 NC

開始第NC 階段分析

計算第一個增量代所需資料

0 I

初始增量、增量位移 、弧長Q _

A

CMAX

C N

N 

B

新增量_I _I1

二分法求 turning points

更改力負荷_f 或

下一增量所需資料

記錄原來負荷與位移

給一固定力負荷_f

或位移負荷_d

目前平衡位置的r 、₀ r_T

正交軌線法 新增量、 、Q _

計算新增量

新增量位移 Q

新弧長_

1

 I

I 2

3

C

NO 選擇求解

種類 平衡路徑

Qmax

Q_C 

YES 臨界點

 

   

     

       

     

                               

A

B B

1

3

其中 是第 I 個增量所使的弧長，_I K_max是給定的每一增量 最大迭代數，K 是第 I 個增量迭代至平衡所使用的迭代數。 _I  正負號之決定方法如下：若 Qr_T 0，則 的正負號與

第 I 個增量相同，反之則相反。

        因為參考自由度改變，利用正交軌線法求新增量猜值^{ 、}  、 Q 

t T T t T

r r

r r

 





1

 0 (3.1.48)

Qr_T r₀  Q ²

       

第四章 挫屈梁的初始變形及自由振動分析

圖 1-1（a）所示為一長方形斷面的細長直梁，其長為L_T，寬為 b，高 為 h，該梁兩端受到側向夾持，一端受軸向壓縮 ，當壓縮量到達挫屈臨界 值，會造成側向挫屈，如繼續施加軸向壓縮，則該梁會如圖 1-1（b）所示 朝側向拱起，文獻上通常稱其為挫屈梁，稱該變形為挫屈梁的初始變形，

並用 /h表示挫屈的程度，其中 梁中點的側向位移，文獻上通常假設其形 狀為該梁的第一挫屈模態[1,4-6,28,34,36-38,40,43]。文獻上有不少挫屈梁振 動的研究[1,4-6,11,18,19,27,28,31,34-38,40,41,43-45]，其中大都是 /h較小 (例如 /h10)的淺挫屈梁，文獻[27,28]用淺拱（shallow arch）理論探討淺 挫屈梁之自然振動頻率及模態的解析解，並以實驗驗證其結果。文獻[43]

中，假設挫屈梁之側向振動位移為該挫屈梁對應之直梁之前兩個挫屈模態 的組合，用淺拱理論、Galerkin’s method 及假設模態法，探討淺挫屈梁

（ /h1.5）之自然振動頻率及受基礎振動的非線性振動，並且用實驗驗 證其分析結果。由文獻[11,27,28,43]的分析及實驗可以發現，當 /h很小時，

不 同 細 長 比 的 挫 屈 梁 ， 若 /h 相 同 ， 則 有 相 同 的 無 因 次 自 然 頻 率

2 1 4) /(EIAL_T

 ，其中為自然頻率、E 為彈性模數，I 為梁斷面的二次矩、

為密度、A 為斷面積，L_T為長度。當 /h很大（例如 /h100）時，文

為密度、A 為斷面積，L_T為長度。當 /h很大（例如 /h100）時，文

在文檔中 挫屈梁在側向負荷下的幾何非線性分析 (頁 36-0)