國
立
交
通
大
學
機械工程學系碩士班
碩 士 論 文
平面三角形殼元素之改善研究
A study of improvement for the triangular plane
shell element
研 究 生：林寬政
指導教授：蕭國模 博士
平面三角形殼元素之改善研究
A study of improvement for the triangular plane shell element
研 究 生： 林寬政 Student： Kuan-Zheng Lin
指導教授： 蕭國模 博士 Advisor： Dr. Kuo-Mo Hsiao
國 立 交 通 大 學 機 械 工 程 學 系 碩 士 班
碩 士 論 文
A Thesis
Submitted to Department of Mechanical Engineering College of Engineering
National Chiao Tung University in Partial Fulfillment of the Requirements
for the Degree of Master of Science
in
Mechanical Engineering August 2010
Hsinchu, Taiwan, Republic of China
平面三角形殼元素之改善研究
研究生：林寬政 指導教授：蕭國模博士國立交通大學機械工程學系碩士班
摘要
本文主要目的為改善平面三角形殼元素在薄殼結構之幾何非線性分析 的精確性。本文以共旋轉(co-rotational formulation)有限元素法來探討薄殼的 幾何非線性行為。本文採用文獻上具旋轉自由度的三角形平面元素與三角 形板元素疊加成一個 3 節點的平面三角殼元素，該元素的節點自由度為 3 個位移、3 個旋轉、及 3 個平面應變。本文中提出一個決定三角殼元素之節 點變形參數的方法。本文以共旋轉 total Lagrangian 法推導平面三角殼元素 幾何剛度矩陣，並以殼結構之切線剛度矩陣的行列式值來偵測平衡路徑上 的分歧點及極限點。 本 文 採 用 牛 頓 - 拉 福 森 (Newton-Raphson) 法 和 弧 長 控 制 (arc-length control)法的增量疊代法來解結構的非線性平衡方程式。本研究以文獻上的 數值例題探討本研究採用之平面三角殼元素的性能，來說明本文提出的方 法的正確性及功效，同時探討不同的元素幾何剛度矩陣對平衡迭代和偵測 平衡路徑上分歧點及挫屈模態的影響。A study of improvement for the flat triangular shell element
Student：Kuan-Zheng Lin Advisor：Dr. Kuo-Mo Hsiao
Department of Mechanical Engineering
National Chiao Tung University
ABSTRACT
The objective of this paper is to improve the accuracy and efficiency of the flat triangular shell element for the geometric nonlinear analysis. In this paper, the co-rotational finite element formulation is employed. The 3-node triangular shell element employed here is the superposition of the triangular membrane element with drilling degrees of freedom and triangular plate element proposed in the literature. The element has nine degrees of freedom per node: three translations, three rotations, and three membrane strains. A motion process to determine the element deformation nodal rotations is proposed. A co-rotational total Lagrangian formulation is used to derive the geometric stiffness matrix of the triangular shell element. The zero value of the tangent stiffness matrix determinant of the structure is used as the criterion to detect the buckling state.
An incremental-iterative method based on the Newton-Raphson method combined with constant arc length of incremental displacement vector is employed for the solution of nonlinear equilibrium equations. Numerical examples are presented to investigate the accuracy and efficiency of the proposed method. The effect of different geometric stiffness matrices derived using different approximations on the convergence rate of equilibrium iteration
and the value of the buckling load are also investigated through numerical examples.
致謝
衷心感謝指導教授 蕭國模博士在這兩年碩士班期間的指導與教誨，使本 論文得以順利的完成。老師在研究上嚴謹的態度及日常生活上的關懷，使我受 益良多，在此致上最高的謝意與敬意。感謝尹慶中老師及蔡佳霖老師撥冗擔任 口試委員並對本論文所提出的指正與建議，使本論文能夠更臻完善。 感謝蔡明旭、周裕淳學長的照顧，同學林運融、蔡秉宏以及學弟盧致群、 翁林甫、黃楚璋在學業上的砥礪與成長。 感謝父母親等關心我的朋友對我的支持與鼓勵，僅以此成果與榮耀，獻給 我親愛的父母以及所有關心我的人。目 錄 中文摘要 ... i 英文摘要 ... ii 致謝 ... iv 目錄 ... v 表目錄 ... vii 圖目錄 ... viii 第一章 緒論 ... 1 第二章 理論推導... 9 2.1 基本假設... 9 2.2 座標系統... 10 2.3 旋轉向量... 11 2.4 元素節點變形參數的決定方法... 11 2.5 殼元素的變形描述... 14 2.6 QST 平面元素的變形描述 ... 14 2.7 DKT 板元素的變形描述... 18 2.8 元素內力與元素剛度矩陣... 22 2.9 QST 平面元素之節點內力與剛度矩陣 ... 23 2.10 DKT 板元素之節點內力與剛度矩陣 ... 24 2.11 元素幾何剛度矩陣... 25 2.12 座標系統轉換... 26 2.13 系統平衡方程式與收斂準則... 29 第三章 數值計算方法與程序... 30 3.1 增量迭代法... 30 3.2 弧長控制法... 32
3.3 二分法... 32 3.4 數值程序... 33 第四章 數值分析與結果... 35 4.1 半圓環受到單點集中力作用 ... 36 4.2 直角構架受到端點剪力作用... 37 4.3 圓柱殼片段受到單點集中力作用... 38
4.4 Lateral torsional buckling ... 39
4.5 Simply supported compressed plate... 40
4.6 Transversally loaded T profile ... 40
4.7 Channel section in torsion... 41
4.8 直角梁受到單點集中力作用... 42 4.9 懸臂圓柱殼受到單點集中力作用... 43 4.10 半球殼受到單點集中力作用... 44 第五章 結論與展望... 45 5.1 結論... 45 5.2 未來研究方向... 46 參考文獻 ... 48 附錄 A 面積座標 ... 82 附錄 B 不完整三階埃爾米特元素的形狀函數及其微分... 85 附錄 C DKT 元素的形狀函數... 87 附錄 D 元素幾何剛度矩陣... 93 附錄 E 等效節點外力 ... 100 附錄 F Battini 所使用的元素座標系統... 106
表目錄
表 4.1 圓柱殼片段受到單點集中力作用所偵測到的挫屈負荷
(例題 4.3) ... 54
表 4.2 Lateral torsional buckling 所偵測到的挫屈負荷(例題 4.4) ... 54 表 4.3 Simply supported compressed plate 所偵測到的挫屈負荷
(例題 4.5)... 55
表 4.4 Transversally loaded T profile 所偵測到的挫屈負荷
(例題 4.6)... 56
圖目錄 圖 1.1 文獻[37]實驗所觀察到四種變形轉換(a-d)及 對應於 a-c 圖結構的上視圖(e-g)... 57 圖 2.1 固定總體座標、元素座標與節點基礎座標... 58 圖 2.2 旋轉向量... 59 圖 2.3 元素變形示意圖... 60 圖 2.4 元素節點 j中心面之 B ij x 軸受旋轉向量−θnj作用的示意圖... 61 圖 2.5 元素節點 j中心面之 E zj 3e θ 為將 B j x1′ 0 軸旋轉到 B j x1′軸旋轉向量 的示意圖... 62 圖 2.6 三角元素的示意圖及節點自由度... 63 圖 2.7 DKT 元素的節點及其三邊上的局部座標示意圖 ... 64 圖 4.1 半圓環受到單點集中力作用 (a)結構尺寸示意圖 (b)網格 18×1 示意圖 (c)網格 18×2 示意圖 ... 65 圖 4.2 半圓環受到單點集中力作用，不同網格下之負荷位移 曲線圖。... 66 圖 4.3 直角構架受到端點剪力作用 (a)結構尺寸示意圖 (b)網格 M21 與網格 M22 示意圖... 67 圖 4.4 直角構架受到端點剪力作用，不同網格下之負荷位移曲 線圖。 ... 68 圖 4.5 圓柱殼片段受到單點集中力作用 (a)結構尺寸示意圖 (b)力負荷圖 (c)網格 10×10 示意圖 ... 69 圖 4.6 圓柱殼片段受到單點集中力作用之負荷位移曲線圖 ... 70
圖 4.7 Lateral torsional buckling (a)結構尺寸示意圖 (b)網格(1+2+1)×2 示意圖... 71
(b)網格 4×6 示意圖... 72
圖 4.9 Transversally loaded T profile (a)結構尺寸示意圖
(b)網格(2+2+3)×4 示意圖... 73
圖 4.10 Channel section in torsion (a)結構尺寸示意圖
(b)網格(1+2+1)×2 示意圖... 74
圖 4.11 Channel section in torsion，不同網格下之負荷位移
曲線圖。 ... 75 圖 4.12 懸臂圓柱殼受到單點集中力作用 (a)結構尺寸示意圖 (b)網格 16×16 示意圖... 76 圖 4.13 懸臂圓柱殼受到單點集中力作用之負荷位移曲線圖 ... 77 圖 4.14 半球殼受到單點集中力作用 (a)結構尺寸示意圖 (b)網格 12×12 示意圖... 78 圖 4.15 半球殼受到單點集中力作用之負荷位移曲線圖... 79 圖 4.16 直角梁受到單點集中力作用 (a)結構尺寸示意圖 (b)網格(2+2)×2 示意圖 ... 80 圖 4.17 直角梁受到單點集中力作用，不同網格下之負荷位移 曲線圖。... 81
第一章 第一章 第一章 第一章 緒論緒論緒論緒論 在近代工程設計的發展上，對材料的要求與結構的表現趨向於高強度 與輕量化，舉凡建築結構、航太設備、壓力容器、軍事載體及汽車工業等， 設計者要考量以最小成本來達到所需的機能與強度，並且兼顧外型的美 觀。由於薄殼在承受彎曲應力與拉伸應力的表現上有良好的表現，又能達 到經濟與輕量化的要求，因此薄殼為設計者最常使用的結構之ㄧ，而廣泛 應用在工程及生活上。常見的薄殼結構包括：建築屋頂、飛機蒙皮、液體 儲存槽、人造衛星、火箭、船體結構、水中潛體等。 由於薄殼結構常用來承受大旋轉和大位移，即使在彈性範圍內，其位 移與受力的關係通常呈非線性變化，因此必須使用非線性分析來探討由幾 何形狀改變所造成的非線性行為。殼的研究從最早的線性分析進入到較實 際的材料非線性及大變形問題；後來殼的不穩定性、挫屈及挫屈後行為和 幾何非線性等問題則陸續被探討。結構在受到力負荷及位移負荷的分析 中，負荷參數對位移曲線圖的平衡路徑有時會出現snap-back與snap-through 的情況[1]。結構在挫屈前僅有一平衡路徑，即主要平衡路徑；但如果出現 分歧點(bifurcation point)，則在挫屈後的平衡路徑除了主要平衡路徑外，還 有次要平衡路徑，主要和次要平衡路徑上的變形型態可能完全不同。 分析殼結構最常用的方法為有限元素法，世界各地的學者在薄板和薄 殼的有限元素模擬研究已超過五十年，會花這麼長時間在這方面研究主要
可以用兩個理由來解釋[2]，第一，在工業上固體力學問題的有限元素模擬 大約佔百分之七十都當成薄殼結構來處理，第二，一些挑戰性理論和數值 問題吸引很多研究者提出了「簡單、有效率、健全、可靠、容易收斂、精 確、最好…」的板殼元素互相競爭比較。 殼元素最常使用的形狀為三角形或四邊形，對任何不規則形狀的殼結 構，我們都可以輕易地將其切割成有限的三角形組合，但不一定適合將其 切割成四邊形的組合，故三角形元素在文獻上被廣泛的探討及使用。文獻 中有很多殼元素被提出[3-17]，大致上可分為平面元素、曲面元素和等參數 元 素 這 三 類 。 有 限 元 素 分 析 薄 殼 最 常 見 的 做 法 就 是 將 三 角 平 面 元 素
(membrane element)與三角板元素(plate bending element)疊加成一平面三角
殼元素(shell element)，這種殼元素推導簡單，並且數值計算比曲面元素更 有效率[14]，已分別應用在殼結構的線性分析[3-5,17]和幾何非線性分析
[7-17]上。
結構之幾何非線性分析常用的推導法有全拉格朗日法(total Lagrangian
formulation)[9]、更新拉格朗日法(updated Lagrangian formulation)[6,9,17]、
及共旋轉法(co-rotational formulation) [7,8,10-16]，全拉格朗日法是以結構初 始狀態為參考位置來表示總位移及旋轉；更新拉格朗日法是以結構上一個 平衡狀態為參考位置來表示增量位移及旋轉；共旋轉法是利用建立在元素 當前變形位置的元素座標將剛體位移及旋轉從總位移及旋轉中扣除，剩下 的總變形位移及旋轉仍為小位移及旋轉。所以若使用共旋轉法，則在線性
分析使用的元素可以用在大位移小應變問題之幾何非線性分析。共旋轉法 在梁與殼結構的幾何非線性分析已經被廣泛的使用[7,8,10-16]，但在平面應 變的問題，則很少見到被使用[18,26]。
CST(constant strain triangle)元素和LST(linear strain triangle)元素都是最
簡單的平面元素，它們常常跟合適的板元素疊加在薄殼分析上使用，因為 這種平面三角殼元素缺少旋轉自由度(drilling degree of freedom)，所以其元 素剛度的面內旋轉剛度(in-plane rotational stiffness)為零，為了避免系統剛度 矩陣因奇異性(singularity)造成分析的困難，文獻[6,7,11,12]是加上一人工 (Artifical)的面內旋轉剛度。 從1964至1983年期間，許多人在研究如何在3節點的平面三角形元素上 加入節點旋轉自由度，希望能得到一個3節點9個自由度的平面三角形元 素，節點自由度為2個位移、1個旋轉，但是都沒有得出可用的元素，一直 到在1984年Allman[19]才提出第一個成功帶有旋轉自由度、3個節點9個自由 度的平面三角形元素。三角平面元素加入旋轉自由度的原因為[24]：改善三 角平面元素的性能並避免使用到三角形邊上的節點，因為邊上的節點會影 響到網格生成，而且在模擬非線性分析與動態分析時較為困難；當三角平 面元素與三角板元素疊加時，能滿足物理上一個節點有3個旋轉自由度的要 求；三角形元素與殼元素、板元素或是梁元素同時使用時，能使接合簡單 化。 從 1984 年至今，有許多具旋轉自由度平面三角形元素被提出[19-22]，
其中除了文獻[22]的 DLST 元素是 12 個自由度，三角形頂點自由度為 2 個 位移、1 個旋轉並且邊上中點自由度為 1 個位移，其餘的元素都是 3 個節點 9 個自由度的元素。由文獻[22]的結果可以發現當元素網格很密時，所有元 素的結果都差不多，但是在使用同樣數目的元素網格時，DLST 元素的結果 比具旋轉自由度之 3 節點 9 個自由度之元素的結果好，這應是合理的，因 為 DLST 元素有 12 個自由度。1966 年 Felippa[23]提出具 3 節點 18 個自由 度的 QST 平面元素(quadratic strain triangle element)， QST 元素中有一種的 節點自由度為 2 個位移、1 個旋轉及 3 個應變，在文獻[24]中稱其為 QST-3/18RS 元素。文獻上使用 QST 元素分析例題的結果甚為少見，這可能 是因為現在的趨勢是使用低階元素，QST 元素不滿足簡單元素的要求，但 由文獻[25]的例題，可以發現 QST 元素有很好的性能，由文獻[26]的例題， 可以發現在使用同樣數目的元素網格時，QST 元素的結果比具旋轉自由度 之 3 節點 9 個自由度及 12 個自由度之元素的結果好，這應是合理的，因為 QST 元素有 18 個自由度。廣義的來說，應變也是物理上的自由度，所以該 QST 元素應可稱為高階的簡單元素。文獻[26]提出了具旋轉自由度之三角平 面元素的共旋轉推導法，將 QST-3/18RS 元素應用於平面應變問題的幾何非 線性分析，並提出一個決定元素節點變形參數的方法，由文獻[26]例題的結 果，可以發現其方法確實可應用在幾何非線性分析上，並有正確的結果。 但文獻[26]在平衡迭代時使用的剛度矩陣並沒有包含幾何剛度矩陣。 從1965年至今，有許多三角板元素被提出[27-36]，在1966年Bazeley[27]
提出了BCIZ元素，其側向位移場w僅滿足C0連續，BCIZ元素和HCT[28]、
DKT (Discrete Kirchoff Triangle) [29]、HSM (Hybrid Stress Model) [30,31] 元
素一樣都是3個節點9個自由度，其節點自由度為1個位移、2個旋轉，Batoz[32] 研究過幾個克希霍夫板元素BCIZ、DKT、HSM、HCT元素，認為DKT元素 是對薄板分析最可靠的元素，但是在1998年Y. K. Cheung[33]指出DKT元素 的缺陷，即該元素內部不能滿足 x y y x ∂ ∂ = ∂ ∂θ θ 的連續條件且沒有定義側向位移 場w，因為計算板元素之質量矩陣及幾何剛度矩陣需要該元素的側向位移場 w，所以其質量矩陣和幾何剛度矩陣皆使用其他元素之側向位移場來推導， 在薄板彎曲振動及挫屈分析的準確性不能令人滿意。文獻[6,7,11,12]在計算 該元素的幾何剛度矩陣時，假設側向位移場w為線性位移場，並且只考慮平 面元素產生的應力，但文獻[6,7,11,12]中僅有幾何非線性分析，並無挫屈分 析。Y. K. Cheung在1998年[33]提出RDKT元素(refined triangular discrete
Kirchhoff plate bending element)，其元素定義的側向位移場w能滿足C1連
續，RDKT元素的元素幾何剛度矩陣一部分由BCIZ元素的側向位移場w提 供，另一部分由私下假設 x x w =θ ∂ ∂ ， y y w =θ ∂ ∂ 提供，並且加入控制參數來分配， 由文獻[33]中薄板挫屈分析的結果可以發現RDKT元素比DKT元素有更好 的結果，但是無挫屈後分析。在2008年Kasparek[36]提出E30元素，E30元素 是利用拉格朗日插值法(Lagrange interpolation)將側向位移場w假設為三次 完整多項式，在線性分析的結果比HCT元素還要準確。在2007年Khosravi
等人在文獻[14]將OPT元素與DKT元素疊加成一平面三角殼元素，並且使用 LST-Ret平面元素[24]的位移場u、v及BCIZ板元素[27]的側向位移場w推導 該元素的幾何剛度矩陣，該元素在幾何非線性分析跟解析解相當接近。 Khosravi等人在文獻[15]用不同的平面元素Allman(3M)[21]、LST-Ret、OPT 與DKT元素疊加分析薄殼受均布力的例題，將結構分割成不同網格去測 試，文獻[15]從其負荷參數對位移曲線圖發現網格變多之下不同元素得到的 分歧點也會比較接近。在2007年Battini[13]將OPT元素與DKT元素疊加成一 平面三角殼元素，以共旋轉推導法來推導該元素的幾何剛度矩陣，使用較 少的元素在幾何非線性分析、挫屈分析以及挫屈後分析都有不錯的結果。 文獻[37]以實驗和數值方法[38]探討一聚酯圓柱薄殼受位移負荷作用後 的非線性行為，模擬一矩形薄板在長邊以夾鉗挾持，夾鉗先將薄板彎成一 圓柱狀，再施加一集中位移負荷於結構中心的情況，採用兩階段的位移負 荷分析。在其實驗中隨著位移負荷的增加，結構連續產生四個特殊的幾何 變形，如圖1.1所示。第一個變形是在薄殼中心附近出現兩個對稱X、Y軸的
d-cone (developable cone) (圖1.1a)。第二個變形中出現兩個新的d-cone，而
四個d-cone圍成一個對薄殼中心轉了一個角度的菱形(圖1.1b)。第三個變形 為四個d-cone的連線形成一個梯型(圖1.1c)。第四個變形為梯形底邊兩個
d-cone移到薄殼自由端的邊界時，產生一個不連續的變化，使薄殼變成波浪
狀的圓柱殼(圖1.1d)。雖然文獻[37]的負荷-位移曲線圖顯示其數值解與實驗 的曲線相當接近，但是由其數值結果卻無法觀察到實驗中出現的菱形旋轉
及梯形的變形轉換。文獻[12]將CST元素與DKT元素疊加成一平面三角殼元 素來分析，但是與文獻[37]的數值結果一樣，無法觀察實驗中出現的現象。 據本人所知，目前還沒有文獻用數值分析模擬出該實驗出現的各種現象， 所 以 該 實 驗 的 數 值 模 擬 可 以 用 來 檢 驗 殼 元 素 性 能 優 劣 之 基 準 問 題 (benchmark problem)。 一個好的殼元素應能精確的計算出殼結構之非線性平衡路徑及偵測其 上的分歧點及極限點。因一般共旋轉法採用的幾何剛度矩陣僅是近似的幾 何剛度矩陣，故本研究擬用精確的方法推導幾何剛度矩陣，並探討採用不 同位移場推導幾何剛度矩陣對平衡迭代及挫屈負荷的影響。文獻[12]所使用 的元素缺乏旋轉自由度，可能受其影響導致無法觀察實驗中出現的現象。 如前所述，QST元素具有旋轉自由度，也比其他的平面元素有更好的精度， 但文獻上尚未有人採用QST元素與板元素結合成殼元素，故本研究將採用文 獻上具旋轉自由度的QST三角形平面元素[26]與DKT三角板元素[32]疊加成 一平面三角殼元素，並使用共旋轉法，將該元素用在具大位移及大旋轉的 薄殼結構分析中，推導更精確的元素幾何剛度矩陣，並且用當前變形位置 重新計算QST元素的線性剛度矩陣，希望能有效地改善殼結構幾何非線性分 析中平衡迭代的收斂速度，並以精確的切線剛度矩陣來偵測平衡路徑上的 分歧點及挫屈模態。 本文將在第二章介紹本研究所使用的平面三角殼元素，以及推導元素 幾何剛度矩陣。在第三章說明本研究所採用的數值計算方法以及分析時的
數值程序。在第四章中將以非線性例題測試本研究所使用的平面三角殼元 素的性能，來說明本文提出的決定元素節點變形參數的方法是可行的。
第二章 第二章 第二章
第二章 理論推導理論推導理論推導理論推導
本文使用共旋轉推導法(co-rotational formulation)[7]，本文將文獻[26]
上具旋轉自由度的 QST 平面元素(quadratic strain triangle element)與 DKT 板 元素(discrete Kirchhoff theory element)[32]疊加成一個 3 節點 27 個自由度的 三角殼元素，元素的節點自由度為 3 個位移、3 個旋轉及 3 個平面應變。為 了文章的完整性，在本章將簡單描述殼元素變形的假設、內力和剛度矩陣 的推導。不同的元素若有相同的節點及節點自由度，在共旋轉推導法中， 可以用相同的方法決定其節點變形參數，在本章中將提出一個決定三角殼 元素之節點變形參數 (位移、旋轉及應變)的方法。本章中將以較精確的方 法推導平面三角殼元素的幾何剛度矩陣。 2.1 基本假設基本假設基本假設 基本假設 本文中對非線性平面三角殼元素的推導，做以下假設： (1) 薄膜變形及彎曲變形之間無耦合作用。 (2)元素的變形位移及旋轉為小位移及小旋轉。 (3)在元素變形前，垂直於元素中心面的法向線段，在元素變形後，依然保 持直線，且沒有伸長或縮短，該線段在元素的三個頂點垂直於變形後的中 心面，但其餘的點依不同元素可有不同的假設。
2.2 座標系統座標系統座標系統 座標系統
為了描述系統的運動及元素的變形，本文定義了三組座標系統：
(a) 固定總體座標系統(global coordinate system)：XiG
(
i=1,2,3)
如圖2.1所示，結構體所有節點的座標、系統的邊界條件與其他座標系統的
基底，以及結構的平衡方程式，均在此座標系統中定義。在XiG座標系統內
之座標值以
(
X,Y,Z)
表示。(b) 元素座標系統(element coordinate system)：xiE
(
i=1,2,3)
如圖2.1所示，此座標系統是建立在每一個元素變形後的最新位置上，其座 標原點為元素節點 1， E x1 軸為元素節點 1 與元素節點 2 在元素平面上的連 線， E x2 軸是在元素平面上垂直於 E x1 軸，且朝著元素節點 3 的方向， E x3 軸則 是由 E x1 軸與 E x2 軸外積而得。元素的位移、元素變形、元素節點內力與元素 剛度矩陣是在此座標系統中定義，然後經由座標轉換，將其轉換至總體座 標系統及基礎座標系統。本文中 E i x 0 表示初始未變形時的元素座標， E i I x 表 示第I 個增量迭代收斂後的元素座標， E i x 表示當前變形位置的元素座標。 在xiE座標系統內之座標值以
(
x,y,z)
表示。(c) 節點基礎座標系統(base coordinate system)：xijB
(
i=1,2,3)(
j=1,2,3)
如圖 2.1 所示，圖中 j 為節點在元素中的編號，此座標系統的原點是剛接在 結構離散後的每一個節點，並與對應的節點一起移動及旋轉。本文中節點 j 在其初始位置之 B j x3 軸的方向為曲面在該節點的法線方向，x1Bj軸、x2Bj軸的 方向為曲面在該節點互相垂直的切線方向，節點的應變自由度是在此座標
系統中定義。本文中 B ij x 0 表示元素節點 j 在初始未變形時的節點基礎座標， B ij I x 表示元素節點 j在第I 個增量迭代收斂後的節點基礎座標， B ij x 表示元素 節點 j在當前變形位置的節點基礎座標。 2.3 旋轉向量旋轉向量旋轉向量 旋轉向量 本文中使用旋轉向量來表示一個有限旋轉，如圖 2.2 所示，一向量R受到一 旋轉向量 nφ 的作用而轉到一個新位置R′，R與R′之間的關係可表示成[39] ) ( sin ) )( cos -(1 cos R n R n n R R′= φ + φ ⋅ + φ × (2.1) 2.4 元素節點變形參數的決定方法元素節點變形參數的決定方法 元素節點變形參數的決定方法元素節點變形參數的決定方法 本文中採用增量迭代法解非線性平衡方程式，假設第I 個位置已知，此 處的第I 個位置，是指第I 個增量迭代收斂後的平衡位置。xj 為元素節點 j (j =1,2,3)在當前元素座標上的位置向量，I j X 、 B ij I x 、∆Uj、 B j ε ∆ 以及∆φφφφj 分別為元素節點 j在固定總體座標中第I 個位置的位置向量、元素節點 j在 第I 個位置的節點基礎座標、元素節點 j的增量位移向量、元素節點 j的增 量應變向量以及元素節點 j的增量旋轉。本文中假設元素節點 j受∆Uj、 B j ε ∆ 以及∆φφφφj作用後的變形過程如下： 1. 節點 j受∆Uj的作用由第I 個位置I j X 平移到節點當前變形位置Xj，在 移動過程中，節點 j無剛體旋轉，即元素各邊在節點的切線方向維持不 變，剛接在其上的節點基礎座標的方位亦不變。 2. 節點 j及剛接在其上的節點基礎座標軸I B x 受到∆φφφφj作用，旋轉到當前變
形位置的 B ij x 。 3. 節點 j在當前變形位置的節點基礎座標上之應變分量增加一增量 B j ε ∆ 。 由上述的變形過程可知節點 j當前變形位置的固定總體座標Xj可以由I j X 加上∆Uj得到，由Xj可以利用元素座標定義求出當前變形位置的元素座標 E i x 。如圖 2.3 所示，將元素變形後當前變形位置的元素座標 E i x 與初始未變 形時的元素座標 E i x 0 重疊，則元素節點變形參數可決定如下： (1) 節點位移uj 節點位移uj可由元素節點 j在當前變形位置的元素座標上的位置向量xj與 元素節點 j在初始未變形時的元素座標上位置向量的0xj的差得出： j j j x x u = −0 (2.2)
{
j j j}
j = u v w u (2.3){
j j j}
j = x y z x (2.4){
j j j}
j 0x 0y 0z 0 = x (2.5) 由元素座標定義的方式可知0zj =zj =wj =0。 (2) 節點應變εj 利用第I 個位置在節點基礎座標上的節點應變向量 B j I ε 加上增量應變向量 B j ε ∆ 可以求出當前變形位置的節點基礎座標上的節點應變向量： B j B j I B j ε ε ε = +∆ (2.6){
B}
xyj B yj B xj B j ∆ε
∆ε
∆γ
∆ε = (2.7) 其中 B j ε 為當前變形位置的基礎節點座標上的節點應變向量。(3) 節點變形角θij 如圖 2.4 所示，將旋轉向量−θnj作用在xijB軸上，使其旋轉到xijB′軸的新位置， nj θ 可表示為： B j E B j E B j E nj 3 3 3 3 3 3 1 ) ( sin e e e e e e θ × × × = − (2.8) 其中 E 3 e 及e 分別為當前變形位置的元素座標之3Bj x3E軸的單位向量，及當前變 形位置的節點基礎座標之 B j x3 軸的單位向量。此時，x3Bj′軸與x3E軸重合，且x1Bj′ 軸和 B j x2′軸與x3E軸垂直，與x1E軸、x2E軸共平面，但不重合。 同理，在初始未變形時將旋轉向量−0θnj作用在0xijB軸上，使其旋轉到0xijB′軸 的新位置，0θnj可表示為： B j E B j E B j E nj 3 0 3 0 3 0 3 0 3 0 3 0 1 0 ) ( sin e e e e e e θ × × × = − (2.9) 其中 E 3 0 e 、0e 分別為初始未變形時元素座標之3Bj 0x3E軸的單位向量，以及初始 未變形時基礎節點座標之 B j x3 0 軸的單位向量。若結構初始未變形時為曲面， 則 E x3 0 軸和 B j x3 0 軸不平行，所以 0 θnj ≠ 0 ；若結構初始未變形時為平面，則 E x3 0 軸和 B j x3 0 軸平行，所以 0 θnj = 0 。 如圖2.5所示，θzje3E(j =1 ,2,3)為將0x1Bj′軸旋轉到x1Bj′軸的旋轉向量，θzj可表 示為： E B j B j B j B j B j B j zj 3 1 1 0 1 1 0 1 1 0 1 ) ( sin e e e e e e e ⋅ × × × = − ′ ′ ′′ ′′ θ (2.10) 其中 B j ′ 1 0 e 、e 分別為初始未變形時之1Bj′ 0x1Bj′軸的單位向量，以及當前變形位置
之 B j x1′軸的單位向量。 令 E nj θ 0 及 E nj θ 為將0θnj與θ 分別以初始未變形時的元素座標及當前變形位置nj 的元素座標為基底，節點變形角θij(i =x,y; j=1,2,3)可表示為： E nj E nj yj xj θ θ 0 0 − = θ θ (2.11) 2.5 殼元素的變形殼元素的變形殼元素的變形描述殼元素的變形描述描述描述 如圖 2.6 所示之殼元素中心面上有三個節點，每個節點的自由度至少分 別是 E x1 、 E x2 、 E x3 軸方向的位移uj、vj、w ( j j=1,2,3)，繞 E x1 、 E x2 、 E x3 軸方 向的位移轉角θxj、θyj、θzj( j=1,2,3)，以及應變分量εxj、εyj、γxyj。此殼元
素 的 變 形 可 分 為 薄 膜 變 形(membrane deformation)與 彎 曲 變 形(bending
deformation)兩部份，薄膜變形是來自平面元素，彎曲變形則是來自板元素， 由基本假設(1)可知此元素的變形可由薄膜變形與彎曲變形疊加而成。 在圖 2.6中的元素節點位移u 、j v 及j θzj ( j=1,2,3)是QST平面元素[26] 節點位移，而θxj、θyj及wj( j=1,2,3)為DKT 板元素[32]的節點位移。 因本文使用共旋轉推導法，且假設在元素座標上，元素的變形位移及 旋轉為小位移及小旋轉，所以本章中除了推導幾何剛度矩陣外，都僅考慮 線性的位移－應變關係。 2.6 QST 平面平面平面元素平面元素元素的變形描述元素的變形描述的變形描述的變形描述 本文採用的三角平面元素為文獻[26]中所提出的 QST元素，其位移場為
三次變化，應變場則為二次變化。此元素有 3 個節點，每個節點有 6 個自 由度，此元素可使用彼此間能互相轉換的不同節點參數。本文中採用的節 點參數為節點 j在元素座標 E x1、x2E軸的位移分量uj、vj，應變分量
ε
xj、ε
yj、 xyjγ
及逆時鐘方向的旋轉θzj。但為了方便推導，本文在元素推導時使用的節 點參數是節點 j在 E x1 、x2E軸的位移分量uj、vj以及x1E、x2E軸的位移分量 分別對面積座標(area coordinates)ξ
、η
的微分u,ξj、u,ηj、v,ξj、v,ηj，面積 座標的介紹詳見附錄 A。 此元素的位移場可表示為[40]： ξ q Ntu u= (2.12) ξ q Ntv v= (2.13) } { ξ1 ξ2 ξ3 ξ q q q q = (2.14) 0 0 0 {N1 N2 N3 u = N } 0 0 0 0 0 0 5 6 7 8 9 4 N N N N N N (2.15) 3 2 1 0 0 0 { N N N v = N } 0 0 0 0 0 0 N4 N5 N6 N7 N8 N9 (2.16) } , , , , { j j j j j j j u v uξ uη vξ vη ξ = q (2.17) 其中u =u(ξ,η)、v =v(ξ,η)分別為在 E x1 、 E x2 軸的位移分量，Nu為對應於u 的 形 狀 函 數 ， Nv 為 對 應 於 v 的 形 狀 函 數 ， qξ 為 節 點 參 數 向 量 ， i N (i=1,2,3,L,9)為面積座標的三次函數，其表示式詳見附錄 B，另外本文 中{ }
表行矩陣。由小變形的假設，元素內任意點的正應變、剪應變及逆時鐘方向的剛體旋 轉可表示為[41]： x x =u, ε (2.18) y y =v, ε (2.19) x y xy =u, +v, γ (2.20) ) , , ( 2 1 y x u v − = θ (2.21) 因u、v都是
ξ
、η
的顯函數，所以u,x、u,y與u,ξ、u,η及v,x、v,y與v,ξ、v,η 的關係可表示成(詳見附錄 A)： = y x x , u , u , u , u ξ η ξ A (2.22) = y x x , v , v , v , v ξ η ξ A (2.23) = = 31 31 21 21 y x y x , y , x , y , x x η η ξ ξ ξ A (2.24) j i ij x x x = − (2.25) j i ij y y y = − (2.26) 其中xj、yj分別是元素節點 j在元素座標系統的座標值。 令節點參數向量qx與qθ分別為： } { x1 x2 x3 x q q q q = (2.27) } { θ1 θ2 θ3 θ q q q q = (2.28) } , , , , { j j xj yj xj yj xj = u v u u v v q (2.29)} { j j xj yj xyj zj j u v ε ε γ θ θ = q (2.30) 其中u,xj、u,yj、v,xj、v,yj分別為u,x、u,y、v,x、v,y在元素節點 j之值，
ε
xj、 yjε
、γ
xyj、θzj分別為ε
x、ε
y、γ
xy、θ
在元素節點 j之值。由(2.22)式與(2.23) 式，可求得(2.14)式中qξ與(2.27)式中qx的關係： x xq T qξ = ξ (2.31) = x x x x ξ ξ ξ ξ R 0 0 0 R 0 0 0 R T (2.32) = = 31 31 21 21 31 31 21 21 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 y x y x y x y x x x x ξ ξ ξ A 0 0 0 A 0 0 0 I R (2.33) 由(2.18)至(2.21)式可求得(2.27)式中qx與(2.28)式中qθ的關係： θ θq T qx = x (2.34) = θ θ θ θ x x x x R 0 0 0 R 0 0 0 R T (2.35) − = 0 0 1 0 0 0 1 2 1 0 0 0 0 1 2 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 θ x R (2.36) 將(2.31)式與(2.34)式代入(2.12)式與(2.13)式可得：θ θ θ ξθq N q T Ntu u u= = (2.37) θ θ θ ξθq N q T Ntv v v= = (2.38) = = = ξθ ξθ ξθ θ ξ θ ξ θ ξ θ ξ ξθ R 0 0 0 R 0 0 0 R R R 0 0 0 R R 0 0 0 R R T T T x x x x x x x x (2.39) − − = = 31 31 31 21 21 21 31 31 31 21 21 21 2 0 0 0 2 0 0 0 2 0 0 0 2 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 x x y x x y y y x y y x x x θ ξ ξθ R R R (2.40) 2.7 DKT 板元素板元素板元素板元素的變形描述的變形描述的變形描述的變形描述 圖 2.7 所示為文獻[32]中所提出的 DKT 板元素，節點 1、2、3 是三角 形的三個頂點，節點 4、5、6 為三角形三個邊的中點，這三個中點的自由 度僅在元素推導的過程中暫時使用，在最後不會出現在元素的節點自由 度。在圖2.4中， B j x3 軸的方向為元素節點 j 在當前變形位置的法線方向， B j x3′ 軸的方向為元素中心面在當前變形位置的法線方向，θnj為一在 B j B j x x1 − 2 平面 上的旋轉向量，將−θnj作用在 B ij x 軸可將x3Bj軸轉到 B j x3′軸，同理， B j x3 0 軸的方 向為元素節點 j 在初始未變形時的法線方向， B j x3′ 0 軸的方向為元素中心面在 初始未變形時的法線方向，0θnj為一在 B j B j x x1 0 2 0 − 平面上的旋轉向量，將 nj θ 0 − 作用在 B ij x 0 軸可將 B j x3 0 軸轉到 B j x3′ 0 軸。由 2.1 的假設(3)可知垂直於變形前的 元 素 中 心 面 法 線 向 量 變 形 後 仍 為 直 線 且 長 度 不 變 ， 所 以 當 旋 轉 向 量
nj nj θ θ −0 <<1 時，DKT 元素的位移場可表示成： u =zθy(x,y) v=−zθx(x,y) w= w(x,y) (2.41) 其中x、y、z 為元素上任一點分別在x1E、x2E、x3E軸的座標值，θy是(θn−0θn) 在x1E軸方向的分量，θx是(θn−0θn)在x2E軸方向的分量，θn為當前變形位置 的元素變形角，0θn為初始未變形時的元素變形角，u是在x1E軸方向的位 移，v是在x2E軸方向的位移，w是在x3E軸方向位移。 DKT 元素的應變變形包含面內(in plane)正應變εx、εy與剪應變γxy以
及橫向剪應變(transverse shear strain)γyz、γxz。
因本文假設元素的變形為小變形，εx、εy和γxy可表示成(2.18)至(2.20) 式，γyz、γxz可表示成： z x xz =w, +u, γ γyz =w,y +v,z (2.42) 將(2.41)式代入(2.18)至(2.20)式可得： κ εb =z (2.43) εb ={εx ,εy ,γxy} (2.44) x x y y y x x y, , , , , , {θ −θ θ −θ = κ (2.45) 將(2.41)式代入(2.42)式可得： γ={γxz ,γyz}={w,x +θy , w,y −θx} (2.46) 由(2.41)式可知 w、θy、θx與x3E無關，所以可由(2.45)式知橫向剪應變在厚 度方向為常數。 本文中稱圖2.7中沿著元素邊緣方向s為切線方向，而垂直於元素邊緣
方向n為法線方向。 在文獻[32]中對於其所提出的 DKT 元素做了下列的假設： (1) θy、θx在元素內為二次變化，也就是：
∑
= = 6 1 i yi i y Nθ θ ；∑
= − = 6 1 i xi i x N θ θ (2.47) 其中θyi、θxi是θy、θx在圖 2.7 中節點i的節點值，Ni(i =1,L,6)為形 狀函數，其表示式詳見附錄 C。 (2) 元素的三個頂點以及三個邊的中點滿足克希霍夫板理論(Kirchhoff plate theory)的假設，即 (a) 在三個頂點 γxzi =w,xi +θyi =0 i =1,2,3 (2.48a) γyzi =w,yi −θxi =0 i =1,2,3 (2.48b) 其中wxi、wyi、θyi、θxi分別是 ( ) x w wx ∂ ∂ = 、 ( ) y w wy ∂ ∂ = 、θy、θx在 節點i的值。 (b) 在三個邊的中點 0 , = + −θnk wsk k =4,5,6 (2.49a) 0 , = + nk sk w θ k =4,5,6 (2.49b) 其中θnk、θsk分別是θn、θs在節點k的值，θn與θs分別是θ在 n 與 s 方向的分量，w,sk、w,nk分別是 , ( ) s w ws ∂ ∂ = 、 , ( ) n w wn ∂ ∂ = 在節點 k 的值。 (3) w 在元素邊緣的方向上是呈現三次變化，也就是：sj j ij si i ij sk w w l w w l w, , , 4 1 2 3 4 1 2 3 − + − − = (2.50) 其中wi、wj是w在節點i及 j 的值，w,sk是w,s在節點 k 的值，k =4,5,6 分別為邊 23、邊31、邊12的中點，ij 邊為節點i與節點 j 之間的邊(見 圖2.7)，其中i =1−3, j =1−3且i≠ j。 (4) θs在元素邊緣是呈現線性變化，即： ) ( 2 1 sj si sk θ θ θ = + (2.51) 其中θsk、θsi、θsj分別是θs在節點 k 、i、 j 之值，θs是θ在 s 方向的分 量，在圖 2.7 中節點k =4,5,6分別為邊 23、邊 31、邊 12 的中點。 在圖 2.7 中元素三個邊上的
θ
y、θx與θs、θn之幾何轉換關係可表示成： − = s n y x c s s c θ θ θ θ (2.52) x w, 、w,y與w,s、w,n的幾何轉換關係為： − = s n y x w w c s s c w w , , , , (2.53) 其中c =cosαij、s =sinαij，αij為元素的邊 ij 上的法線n 與ij x1E軸的夾角， 見圖 2.7。 由(2.48)至(2.53)式可以把(2.47)式表示成[32]： b T x y H (ξ ,η)u θ = b T y x H (ξ ,η)u θ =− (2.54) = b u [w1 θx1 θy1 w2 θx2 θy2 w3 θx3 θy3] (2.55)其中ub為 DKT 元素的節點位移，Hx與Hy是對應於元素節點位移的新形狀 函數，其表示式詳見附錄 C，ξ與η是元素內任一點在元素自然座標[]的座 標值，其中1≤ξ ≤0、1≤η ≤0。 將(2.54)式代入(2.45)式可以得到： κ =Bbub (2.56) 其中Bb為 DKT 元素的位移-應變轉換矩陣，表示式為： + + − + − = T y T x T x T y T y T x b y x x x x y A ξ η ξ η ξ ξ η ξ , 3 , 2 , 3 , 2 , 3 , 3 2 1 ) , ( H H H H H H B (2.57) 其中 2 3 2y x A= 為三角形面積。 2.8 元素內力與元素剛度矩陣元素內力與元素剛度矩陣元素內力與元素剛度矩陣 元素內力與元素剛度矩陣 本文中殼元素的節點內力是由 QST 平面元素及 DKT 板元素的節點內 力組合而成，元素剛度矩陣k是由 QST 平面元素剛度矩陣k 、DKT 板元素θ 剛度矩陣k 所疊加而成。本節中將用虛功原理推導 QST 平面元素及 DKTb 板元素的節點內力及剛度矩陣。 在平面應力狀態，等向性線彈性材料的應變與應力關係為 ε E σ= (2.58) } , , {σx σy τxy = σ (2.59) } , , {εx εy γxy = ε (2.60)
− − = 2 1 0 0 0 1 0 1 1 2 ν ν ν ν E E (2.61) yz yz Gγ τ = , τxz =Gγxz (2.62)
其中 E 是楊氏模數(Young’s module)，ν 是蒲松比(Poisson ratio)，G 為剪力
模數。 2.9 QST 平面平面平面元素之節點內力與剛度矩陣平面元素之節點內力與剛度矩陣元素之節點內力與剛度矩陣 元素之節點內力與剛度矩陣 將(2.37)式、(2.38)式代入(2.18)至(2.20)式可得： θ ξθq BT εm = (2.63) } { x y xy m = ε ε γ ε (2.64) + = x t v y t u y t v x t u , , , , N N N N B (2.65) 其中Ntu,x、Ntu,y、Ntv,x、Ntv,y的表示式詳見附錄 B。 將(2.63)式代入(2.58)式可得： θ ξθq EBT σm = (2.66) 由虛功原理可得[42]：
∫
= V m t m t dV σ ε f q δ δ θ θ (2.67) } {θ1 θ2 θ3 θ f f f f = (2.68) } { xj yj xj yj j zj j f f mε mε mγ m θ = f (2.69) 其中fθ為對應於qθ的元素節點內力， fxj、 fyj分別為對應於uj、vj的節點力，mεxj、mεyj、mγj分別為對應於
ε
xj、ε
yj、γ
xyj的廣義節點力矩，mzj為對 應於θzj的節點傳統力矩，V 為元素體積。由於應變的因次(dimension)與旋 轉相同，因此應變所對應的節點力mεxj、mεyj、mγj應與力矩有相同的因次為 一廣義節點力矩，而非應力。將(2.63)式、(2.66)式代入(2.67)式可得 θ ξθ ξθ θ θ θδ
δ
q f = q T∫
B EB T q V t t t t dV (2.70) 由(2.70)式可得： θ θ θ k q f = (2.71) ξθ ξ ξθ θ T k T k = T (2.72)∫
= V T dV EB B kξ (2.73) 其中kξ、kθ 分別為應於qξ、qθ的 QST 元素剛度矩陣。 2.10 DKT 板元素之節點內力與剛度矩陣板元素之節點內力與剛度矩陣板元素之節點內力與剛度矩陣 板元素之節點內力與剛度矩陣 將(2.43)式、(2.56)式代入(2.58)式可得： b b b zEB u σ = (2.74) 由虛功原理可得： dV V b t b b t bf =∫
ε σ u δ δ (2.75) 其中fb是板元素對應於δub的節點內力，V為元素的體積。 將(2.43)式、(2.56)式、(2.74)式代入(2.75)式可得：∫∫
= A b b t b t b b t bf u zB zEB u dzdA u δ δ =∫
A b b b t b t b B D B u dA u δ (2.76) 其中 − − = =
∫
− 2 1 0 0 0 1 0 1 ) 1 ( 12 2 3 2 2 2 ν ν ν ν Eh dz z h h b E D (2.77) 其中h為元素厚度。 由(2.76)式可得： b b b k u f = (2.78)∫
= A b b t b b B D B dA k (2.79) 其中kb是 DKT 元素剛度矩陣。 2.11 元素幾何剛度矩陣元素幾何剛度矩陣元素幾何剛度矩陣元素幾何剛度矩陣[43] 精確的元素幾何剛度矩陣，必須以非線性的位移－應變關係推導，文獻[43]中以完整的 Green’s strain 及 total Lagrangian 推導法推導出一三維元素 之切線剛度矩陣的通式，本節中將採用文獻[43]中推導元素幾何剛度矩陣的 方法，並加以修改、簡化，使其適用於本研究採用的元素，文獻[43]中的詳 細推導及本文之元素幾何剛度矩陣的詳細推導在附錄 D 中。本文的元素幾 何剛度矩陣可表示式成 dV V t
∫
= G MG kσ (2.80) } {G1 G2 G3 G= = y wj y yj y j v y xj y j u x wj x yj x j v x xj x j u j z z z z , , , , , , , , , , N 0 H N H N N 0 H N H N G θ θ θ θ (2.81) = I I I I M y xy xy x σ τ τ σ I為 3×3 單位矩陣 (2.82) b m σ σ σ = + σ ={σx ,σy ,τxy} (2.83) 其中 j = 1,2,3為第 j 個節點，Nuθj、Nvθj為 (2.37)式之Nuθ、(2.38)式之Nvθ對 應於節點 j 的部份，Hxj、Hyj為(2.54)式之Hx、Hy對應於節點 j 的部份，I 為3×3的單位矩陣、0為 1×6的零矩陣，σm、σb為(2.66)式、(2.74)式之QST 平面元素及DKT板元素的應力場。 因為DKT元素的內部沒有定義側向位移場w，所以本文假設w為線性位移 場w= Nwub，其中Nw與CST元素[12]的形狀函數相同，本文除了使用QST 元素本身的位移場u=Nuθ,xqθ、v=Nvθ,xqθ外，還用u、v為線性位移場推導 (2.81)式之元素幾何剛度矩陣，本文推導的兩種元素幾何剛度矩陣都僅是近 似的幾何剛度矩陣，但由本文的數值結果發現這兩種元素幾何剛度矩陣在 非線性分析的平衡迭代時都能提高收斂速度且性相能近，在挫屈分析時也 都能偵測到相當接近的挫屈負荷。 2.12 座標系統轉換座標系統轉換座標系統轉換座標系統轉換 為了建立結構的平衡方程式，必須將元素節點參數q 中的uj、vj、wj和 xj θ 、θyj、θzj轉換成固定總體座標的分量，
ε
xj、ε
yj、γ
xyj轉換成對應於節點 基礎座標的分量，才能將各元素節點內力組合成結構系統節點內力以及將 元素剛度矩陣疊加成結構系統剛度矩陣。在元素座標的節點位移uj、vj、wj 和節點旋轉θ 、θ 、θ 與在固定總體座標的節點位移和節點旋轉之關係式為： = j j j t GE j j j W V U w v u A (2.84) = G Zj G Yj G Xj t GE zj yj xj θ θ θ θ θ θ A (2.85) ] [ 1E 2E 3E GE e e e A = (2.86) 其中Uj、Vj、Wj分別為在X1G、 G X2 、X3G軸的位移分量，θXjG、θYjG、θYjG分 別為在 G X1 、 G X2 、 G X3 軸的旋轉分量，eEi (i=1,2,3)為元素座標xiE
(
i=1,2,3)
方 向的單位向量。 其中在元素座標的節點應變ε
xj、ε
yj、γ
xyj與在節點基礎座標的節點應 變之關係式為[41]： B j EBj j R ε ε = (2.87){
xj yj xyj}
j = ε ε γ ε (2.88){
B}
xyj B yj B xj B j = ε ε γ ε (2.89) − − − = 2 2 2 2 2 2 2 2 Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj EBj s c s c s c s c c s s c s c R (2.90) ] [ 1Bj 2Bj 3Bj GBj e e e A = (2.91) 其中cBj、sBj分別為e 在元素座標1Bj x1E、x2E軸的分量，ε
xjB、ε
yjB、γ
xyjB 為對應 於節點基礎座標 B j x1 、 x2Bj軸的應變分量，eijB(i, j=1,2,3)為節點基礎座標(
i, j=1,2,3)
xijB 方向的單位向量。} { j j j xj yj zj xj yj xyj Ej = u v w θ θ θ ε ε γ q (2.92) 將(2.80)式、(2.81)式、(2.86)式代入(2.90)式可得： B EB E T q q = (2.93) = 3 2 1 EB EB EB EB T 0 0 0 T 0 0 0 T T (2.94)
{
B1 B2 B3}
B q q q q = (2.95) = EBj t GE t GE EBj R 0 0 0 A 0 0 0 A T (2.96) } { j j j XjG YjG ZjG xjB yjB xyjB Bj = U V W θ θ θ ε ε γ q (2.97) 由(2.71)和(2.78)可得： b E f f f = θ + (2.98) } { E1 E2 E3 E f f f f = (2.99) } { xj yj zj xj yj zj xj yj j Ej = f f f m m m mε mε mγ f (2.100) 由反梯度法則(contragradient law)[44]及(2.93)式可得： E t EB B T f f = (2.101) } { B1 B2 B3 B f f f f = (2.102) } { Xj Yj Zj GXj YjG GZj Bxj Byj Bj Bj = F F F m m m mε mε mγ f (2.103) 其中fB為對應於qB的元素節點內力，FXj、FYj、FZj分別為對應於Uj、Vj、 j W 的節點力，mεBxj、 B yj mε 、 B j mγ 分別為對應於 B xjε
、 B yjε
、 B xyjγ
的廣義節點力矩， G Xj m 、mYjG、mGZj為對應於θXjG、θYjG、θZjG的節點傳統力矩。將(2.93)式、(2.98)式代入(2.101)式可得： B B B k q f = (2.104) EB E t EB B T k T k = (2.105) 其中kB是對應於qB的元素剛度矩陣。 2.13 系統平衡方程式與收斂準則系統平衡方程式與收斂準則系統平衡方程式與收斂準則系統平衡方程式與收斂準則 結構系統受外力負荷時，其平衡方程式可表示為
(
Q) ( )
=F Q − P = 0 Ψ ,λ
λ
(2.106) 其中Ψ為系統不平衡力向量，系統節點內力F可由(2.102)式的元素節點內力 B f 疊加得出，Q 為系統位移向量，λ為負荷參數，P為參考外力負荷向量。 若外力為與變形位置相關(configuration dependent)的外力，則P在每一個變形位置都須更新。本文以不平衡力Ψ的 weighted Euclidean norm 做為平衡迭
代時的誤差度量，且收歛準則表示為 tol e N e= ≤ P Ψ
λ
(2.107) 其中N代表離散系統的自由度總數，etol為一給定的容許誤差值。第三章 第三章 第三章 第三章 數值計算方法與程序數值計算方法與程序數值計算方法與程序數值計算方法與程序 本文解(2.106)式的非線性平衡方程式所使用平衡迭代的數值計算方法 是採用文獻[45]中所提出基於牛頓-拉福森(Newton-Raphson)法和弧長控制
(arc length control)法的增量迭代法。本文中為了求得分歧點，以系統切線剛
度矩陣之行列式值為零來判斷。本文採用二分法決定增量位移向量的長 度，以求得系統切線剛度矩陣之行列式值為零的平衡位置。為了求得次要 平衡路徑，本文在平衡路徑的第一個挫屈負荷分歧點加入一個與第一挫屈 模態向量方向一致的擾動位移。為了本文的完整性，以下將簡單介紹文獻[45] 中提出的增量迭代數值計算方法與程序。 3.1 增量迭代法增量迭代法增量迭代法增量迭代法 若第I 個增量的平衡位置為已知，令其位移向量為QI、負荷參數為λI， 則第I +1個增量的初始增量位移向量 Q∆ ，可以利用尤拉預測值(Euler predictor)求得[46] T r Q= ∆λ ∆ (3.1)
![圖 1.1 文獻[37]實驗所觀察到四種變形轉換(a-d)及 對應於 a-c 圖結構的上視圖(e-g)](https://thumb-ap.123doks.com/thumbv2/9libinfo/8433842.181387/68.892.235.655.238.902/圖11文獻37實驗所觀察到四種變形轉換ad及對應於ac圖結構的上視.webp)
![圖 4.2 半圓環受到單點集中力作用，不同網格下之負荷位移曲線圖。051015202530012345678910111213 (a) (b) ref.[18] 20x1 30x1 40x2 80x4 100x6Load F (103)](https://thumb-ap.123doks.com/thumbv2/9libinfo/8433842.181387/77.892.192.768.250.713/圖半圓環受到單點集中力作用不同網格下之負荷位移曲線82x13.webp)
![圖 4.4 直角構架受到端點剪力作用，不同網格下之負荷位移曲線圖。012345670.00.51.01.52.02.53.03.54.0 (a) (b) ref.[18] M21 M22Load F (104)](https://thumb-ap.123doks.com/thumbv2/9libinfo/8433842.181387/79.892.173.770.272.733/圖直角構架受到端點剪力作用不同網格下之負荷位移曲線ref1MMLoad.webp)
