高階平面三角形殼元素之研究

國

立

交

通

大

學

機械工程學系

碩 士 論 文

高階平面三角形殼元素之研究

A study on high order flat triangular shell element

研 究 生：盧致群

指導教授：蕭國模 博士

高階平面三角形殼元素之研究

A study on high order flat triangular shell ekement

研 究 生： 盧致群 Student： Chih-Chun Lu

指導教授： 蕭國模 博士 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 2011

Hsinchu, Taiwan, Republic of China

高階平面三角形殼元素之研究

研究生：盧致群 指導教授：蕭國模博士

國立交通大學機械工程學系碩士班

摘要

本文主要目的為探討一個高階平面三角形殼元素在薄殼結構之幾何非 線性分析的精確性。本文以共旋轉(co-rotational formulation)有限元素法及增 量迭代法來探討薄殼的幾何非線性行為。本文將採用文獻上一個具旋轉自 由度的三角形平面元素與一個 連續的高階三角形板元素疊加成一個 3 節 點的高階平面三角殼元素，元素的節點自由度為 3 個位移、3 個旋轉、3 個 平面應變及 3 個側向位移二次微分。本文以殼結構之切線剛度矩陣的行列 式值來偵測平衡路徑上的分歧點及極限點。 1 C 本文採用牛頓-拉福森(Newton-Raphson)法和弧長控制(arc-length control)法的增量疊代法來解結構的非線性平衡方程式並以數值例題測試該 高階平面三角形殼元素的性能。

A Study on High Order Flat Triangular Shell Element

Student：Chih-Chun Lu Advisor：Dr. Kuo-Mo Hsiao

Department of Mechanical Engineering

National Chiao Tung University

ABSTRACT

A new high order flat triangular shell element for the geometrically

nonlinear analysis are presented. In this paper, co-rotational finite element

formulation and incremental-iterative method are employed. The new shell

element is the superposition of a triangular membrane element with drilling

degree of freedom and a continuous high order triangular plate element.

The element has 3 nodes and 12 degrees of freedom per node. The degrees of

freedom at each node are 3 translations, 3 rotations, 3 membrane strains and 3

second derivative of lateral displacement. The zero value of the tangent stiffness

matrix determinant of the structure is used as the criterion to detect the buckling

state.

1

C

An incremental-iterative method based on the Newton-Raphson method

combined with constant arc length of incremental displacement vector is used

for solving nonlinear equilibrium equations. Several numerical examples are

誌謝

衷心感謝指導教授 蕭國模博士在這兩年碩士班期間的指導與教誨，使本 論文得以順利的完成。老師在研究上嚴謹的態度及日常生活上的關懷，使我受 益良多，在此致上最高的謝意與敬意。感謝蔡佳霖老師及鄭文雅老師撥冗擔任 口試委員並對本論文所提出的指正與建議，使本論文能夠更臻完善。 感謝蔡明旭、劉峰成、林寬正、蔡秉宏、林運融學長的照顧，同學翁林甫 在學業上的砥礪與成長。感謝這一年來與我共處一室的室友在生活上的照顧， 以及留在交清兩校的同學和朋友，感謝你們陪伴我度過這兩年。 感謝父母親和關心我的親戚們對我的支持與鼓勵，僅以此成果與榮耀，獻 給我親愛的父母以及所有關心我的人。

目 錄

中文摘要 ... I 英文摘要 ... II 誌謝 ... III 目錄 ... IV 表目錄 ... VI 圖目錄 ... VII 第一章 緒論 ... 1 第二章 理論推導... 8 2.1 基本假設 ... 8 2.2 座標系統 ... 8 2.3 殼元素的變形描述... 10 2.4 元素内力與元素剛度矩陣... 17 2.5 節點參數在各座標系統之間的轉換... 21 2.6 元素節點變形參數的決定方法... 25 2.7 旋轉向量 ... 28 2.8 系統平衡方程式與收斂準則... 28 第三章 數值計算方法與程序... 30 3.1 增量迭代法 ... 30 3.2 弧長控制法 ... 32 3.3 二分法 ... 32 3.4 數值程序 ... 33 第四章 數值分析與結果... 36

4.2 RQT+QST 元素的幾何非線性與挫屈分析... 38 例題一：圓柱殼片段受集中力作用... 38 例題二：槽型梁受扭矩... 38 例題三：受均勻壓縮位移之簡支板... 39 例題四：槽型梁受均勻位移負荷之側向扭矩挫屈... 40 例題五：T 形斷面梁受集中側力 ... 41 例題六：角型梁受軸向集中力... 42 例題七：懸臂圓柱殼受一對端點集中力... 42 例題八：半球殼受二對集中力... 43 例題九：聚酯圓柱薄殼受兩階段負荷作用... 44 第五章 結論與展望... 46 5.1 結論 ... 46 5.2 未來研究方向 ... 47 參考文獻 ... 48 附錄 A ... 85 附錄 B ... 88 附錄 C ... 91 附錄 D ... 95

表目錄 表 4.1 簡支板受中點集中力 P 之無因次化中點側向位移 2 2 10 PL D W_C ... 52 表 4.2 簡支板受均勻分佈力q₀之無因次化中點側向位移 4 0 3 10 L q D W_C ... 52 表 4.3 固定板受中點集中力 P 之無因次化中點側向位移 2 3 10 PL D W_C ... 53 表 4.4 固定板受均勻分佈力q₀之無因次化中點側向位移 4 0 3 10 L q D W_C ... 53 表4.5 槽型梁受扭矩的極限點... 54 表4.6受均勻壓縮位移之簡支板的挫屈負荷... 54 表4.7槽型梁受均勻位移負荷之側向扭轉的挫屈負荷... 54 表4.8 T形斷面梁受集中側力的挫屈負荷... 54

圖目錄 圖2.1三角元素的示意圖及節點自由度... 55 圖2.2邊上切線和法線方向示意圖... 56 圖2.3元素變形示意圖... 57 圖2.4元素節點

j

中心面之 B軸受旋轉向量 ij

x

θ_nj作用 的示意圖... 58 圖2.5元素節點

j

中心面之 E為將 軸旋轉到 軸旋轉向量 zj 3e



B j x₁ 0 B j x₁ 的示意圖... 59 圖2.6旋轉向量... 60 圖4.1固定板 (a)結構示意圖 (b)網格種類及切割示意圖... 61 圖4.2圓柱殼片段受到單點集中力作用 (a)結構尺寸示意圖 (b)力負荷圖 (c)網格10×10示意圖... 62 圖4.3圓柱殼片段受集中力作用的位移-負荷曲線圖... 63 圖4.4槽型梁受扭矩之 (a)結構尺寸示意圖 (b)網格(1+2+1)×2 示意圖... 64 圖4.5槽型梁受扭矩的位移-負荷曲線圖... 65 圖4.6受均勻壓縮位移之簡支板的 (a)結構尺寸示意圖 (b)網格4×6示意圖... 66 圖4.7受均勻壓縮位移負荷之簡支板的位移-負荷曲線圖... 67 圖4.8 (a)圖4.7中 I 點(



0.596 W_E 0.223) (b)圖4.7 中 J 點(



3.069 W_E 5.142) (c)圖4.7 中 K 點 (



10 W_E 9.017) 側向位移場的等高線圖... 68 圖4.9槽型梁 (a)結構尺寸示意圖

(b)網格(1+2+1)×2 示意圖... 69 圖4.10 T形斷面梁 (a)結構尺寸示意圖 (b)網格(2+2+3)×4示意圖... 70 圖4.11 T形斷面梁受側力之位移-負荷曲線圖... 71 圖4.12角型梁 (a)結構尺寸示意圖 (b)網格(2+3)×2示意圖... 72 圖4.13角型梁受軸向集中力之位移負荷曲線圖... 73 圖4.14懸臂圓柱殼 (a)結構尺寸示意圖 (b)網格16×16示意圖... 74 圖4.15懸臂圓柱殼受一對端點集中力之位移負-荷曲線圖... 75 圖4.16半球殼 (a)結構尺寸示意圖 (b)網格12×12示意圖... 76 圖4.17半球殼受到二對集中力之位移-負荷曲線圖... 77 圖4.18圓柱薄殼 (a)結構示意圖 (b)俯視圖 (c)前視圖... 78 圖4.19圓柱薄殼第一階段受力 E 點之負荷參數-位移曲線圖... 79 圖4.20圓柱薄殼第二階段受力 E 點之位移-負荷曲線圖... 80 圖4.21圓柱薄殼第二階段受力 E 點位移-負荷曲線 之主要路徑與次要路徑... 81 圖4.22 (a)圖4.21中O 點(W_E 0



0) (b)圖4.21中 I 點(W_E 7.5309



6.5386) 側向位移場的等高線圖... 82 圖4.23 (a)圖4.21中 J 點(W_E 9.1726



9.1452) (b)圖4.21中 K 點(W_E 11.5455



14.5136) 側向位移場的等高線圖... 83 圖4.24 (a)圖4.21中 M 點(W_E 9.4736



9.4156) (b)圖4.21中 N 點(W_E 14.5441



9.14933)

第一章 緒論

在近代工程設計的發展上，對材料的要求與結構的表現趨向於高強度 與輕量化，不論機械、建築結構、航太設備及運輸工業等，設計上必須考 量以最小成本及重量來達到所要求的機能與強度，同時兼顧外型的美觀。 由於薄殼在承受彎曲應力與拉伸應力時有良好的表現，又擁有極佳的重量 強度比，因此薄殼結構常被廣泛應用在工程及生活上。目前常見的薄殼結 構有建築屋頂、飛機蒙皮、液體儲存槽、人造衛星、火箭、船體結構等。 薄殼結構受到外力作用經常會產生大位移和大旋轉，在大位移和大旋 轉的問題中，位移和外力往往不是線性關係，因此需要使用非線性分析的 方法來探討由幾何形狀改變所造成的非線性行為。常見的幾何非線性分析 的推導方法有三種：全拉格朗日法、更新拉格朗日法和共旋轉法。全拉格 朗日法是用初始狀態為參考位置來表示總位移和旋轉；更新拉格朗日法是 以結構上一個平衡狀態為參考位置來表示增量位移和旋轉；共旋轉法是利 用建立在元素當前變形位置的元素座標將剛體位移及旋轉從總位移及旋轉 中扣除，剩下的位移和旋轉即為小位移和小旋轉，因此若使用共旋轉法， 原本在線性分析的元素也可以應用在大位移、小應變的幾何非線性分析， 共旋轉法在梁與殼結構的幾何非線性分析已經被廣泛的使用[1-6]。 使用有限元素法分析薄殼結構是一個有系統的方法，但在殼元素的發

用在各種例題的完美殼元素，因此有很多文獻在推導新的殼元素的或是改 良殼元素的性能。殼元素大致分為三類：平面殼元素[1-6]、曲面殼元素和 等參數元素。平面殼元素是由一個平面板元素和一個平面元素疊加而成， 此種元素的推導方式簡單，而且在數值計算上比曲面殼元素更有效率[7]， 已分別應用在殼結構的線性和幾何非線性分析[1-6]。 從形狀分類，殼元素大致分為兩種：四邊形殼元素和三角形殼元素。 三角形殼元素具有容易切割複雜幾何形狀的優點，因此三角形的元素被廣 泛的探討及使用。

CST(constant strain triangle)元素和 LST(linear strain triangle)元素都是最

簡單的平面元素，在薄殼分析上常常使用它們與合適的板元素疊加，因為 這種平面三角殼元素缺少旋轉自由度(drilling degree of freedom)，所以其元 素剛度的面內旋轉剛度(in-plane rotational stiffness)為零，為了避免系統剛度 矩陣因奇異性(singularity)造成分析的困難，常見的解決方法有兩種：(1) 加 上一個人工的面內旋轉自由度[11-13,15]。(2) 採用具有旋轉自由度的平面 元素[8-11]。文獻[10]提到三角平面元素加入旋轉自由度的優點為：改善三 角平面元素的性能並避免使用到三角形邊上的節點，因為邊上的節點會影 響到網格生成，而且在模擬非線性分析與動態分析時較為困難；當三角平 面元素與三角板元素疊加時，能滿足物理上一個節點有 3 個旋轉自由度的 要求；三角形元素與殼元素、板元素或是梁元素同時使用時，能使接合簡

單化。1964 至 1983 年期間，許多人在研究如何在 3 節點的三角形平面元素 上加入節點旋轉自由度，希望能得到一個 3 節點 9 個自由度且具節點旋轉 自由度的三角形平面元素，但是都沒有得出可用的元素。1984 年文獻 [8]Allman 提出第一個成功帶有旋轉自由度、3 個節點 9 個自由度的三角形 平面元素。2002 年文獻[9]提出 DLST 元素是一個具有 12 個自由度，三角 形頂點自由度為 2 個位移、1 個旋轉並且邊上中點自由度為 1 個位移。2003 年文獻[10]提出的 OPT 元素具有 3 個節點，每個節點 2 個位移、1 個旋轉的 自由度。2008 年文獻[11] 成功的將一個具有旋轉自由度的 QST 元素[10]應 用在平面應變問題的共旋轉法幾何非線性分析上，該元素具有 3 個節點、 18 個自由度，每個節點有 2 個位移、1 個旋轉及 3 個應變自由度。由文獻[11] 的例題可以發現在使用同樣數目的元素網格時，使用 QST 元素分析的結果 比其他具有旋轉自由度之 3 個節點、9 個自由度的 Allman、OPT 元素好， 這個結果推測是正確的，因為 QST 元素的自由度比 Allman、OPT 元素多。 從 1960 年起有許多的三角形板元素被提出[5,6,7,8]，其中具有 3 個節 點，每個節點有 1 個位移、兩個旋轉自由度，總共 9 個自由度之三角形板 元素的研究發展最為迅速。 文獻[12]中比較 DKT、HSM、BCIZ、HCT、A-9、STRUDL 等 9 個自 由度三角形板元素的線性分析和振動分析後，認為 DKT 元素是這些三角形 板元素中最有效率的元素，在靜態和模態分析中均可以收斂到準確的答案。

文獻上有許多使用 DKT 板元素加上一個平面元素形成的殼元素進行幾 何非線性分析，1981 年文獻[1]使用 CST 平面元素和 DKT 板元素疊加成一 個殼元素，並用更新拉格朗日法分析大位移及大旋轉的薄殼結構，但更新 拉格朗日法的增量旋轉必須是小角度，因此 1987 年文獻[2]使用相同的殼元 素，搭配共旋轉法解決增量旋轉大小限制的問題。2006 年文獻[3]採用文獻 [2]的殼元素和共旋轉法，以數值例題探討殼結構受到各種位移負荷之幾何 非線性問題。2006 年文獻[4]將兩種平面元素(Allman、OPT)與多種板元素 (DKT、TRIC、AQR)互相搭配疊加，用共旋轉法測試挫屈例題，從文獻[4] 結論可以得知 DKT 元素和 OPT 元素搭配的殼元素在例題裡有比較好的結 果。 由以上文獻得知，DKT 元素在幾何非線性分析中是一個高精度、收歛 快的元素並且已被廣泛的使用，但 DKT 元素還是有些缺陷：文獻[10]中指 出該元素內部不能滿足 _y x w _

_

  、 _x y w_

_

  且沒有定義側向位移場 ，因為計 算板元素之質量矩陣及幾何剛度矩陣需要該元素的側向位移場 ，所以 DKT 元素的質量矩陣和幾何剛度矩陣皆使用其他元素的側向位移場來推 導，因此在薄板彎曲振動及挫屈分析的準確性不能令人滿意。 w w 2007 年文獻[5]沿用文獻[3]中 CST 平面元素和 DKT 板元素疊加成一個 殼元素，配合共旋轉法和位移負荷模擬文獻[13]中的實驗(詳見附錄 A)，但 無法模擬出實驗中 d-cone 旋轉的現象。2010 年文獻[6]使用 DKT 板元素和

具有平面旋轉自由度的 QST 平面元素[11]疊加成一個殼元素配合共旋轉法 分析幾何非線性和挫屈例題，雖然在大部份的例題中有很好的結果，但仍 然無法模擬出文獻[13]中的 d-cone 旋轉的現象，本人推測原因為 DKT 元素 不具有 連續的特性，且 DKT 元素內部沒有定義側向位移場，因此必須使 用其他元素的位移場推導幾何剛度矩陣。文獻[13]使用文獻[14]的元素模擬 實驗的圓柱薄殼，在薄殼變成波浪狀之前，數值模擬和實驗得出的受力-位 移曲線相當吻合，但數值模擬無法得出 d-cone 的旋轉的型態。文獻[14]使用 有限元素法，考慮的薄膜應變為不完整非線性項，並且假設 1 C x、y、 方向 的位移是完整五次多項式，計算 Von Karman 板殼理論推導出來的能量式， 利用找出該能量的最小值求得薄殼結構的變形。文獻[16]推測該數值模擬無 法得出 d-cone 的旋轉的型態，原因可能為離散化的不足。就本人所知，目 前還沒有文獻使用數值模擬得出文獻[13]的 d-cone 的旋轉的型態，因此該實 驗的數值模擬可以用來檢驗殼元素性能優劣的基準問題(Benchmark)。 z 在板元素的發展中，三角形板元素不全是只有 9 個自由度，一個 連 續且位移場是五次高階多項式的三角形板元素在 1965 年以後陸續被提出。 1968 年文獻[15]把位移假設為完整五次多項式的元素稱為 TUBA 6，除了在 三頂點有 1 個側向位移 、2個側向位移一次微分 、 ，3 個側向位 移兩次微分 、 、 之外，還有 3 個邊上中點側向位移對邊上法 線方向微分 ，共計 6 個節點、21 個自由度。1969 年文獻[16]中也假設 1 C i w, w xi w, w,_yi xxi w, nk w, xyi w, ,_yyi

側向位移是五次多項式，並且利用邊上的側向位移對邊上法線方向微分 是三次變化三個限制條件將三個邊上的自由度 、 、 去掉，因此 6 個節點減為 3 個節點、21 個自由度減為 18 個自由度。文獻[16]中有做靜 態和動態的例題測試，從得到的數據可以發現不但收斂速度很快，跟解析 解相比也非常的相近，因此文獻[16]作者 Cowper 認為 18 個自由度的高階板 元素是非常優秀的元素。文獻[17]比較了 6 節點、21 個自由度的板元素和 簡化過後的 3 節點、18 個自由度的板元素在線性例題上的準確度和收斂性， 發現簡化過後的 18 個自由度的板元素在精度和收斂性上相較於 21 度板元 素並沒有喪失很多，又由於只有 3 個頂點為節點的三角形板元素在網格切 割上會比 6 個節點、邊界上有節點的三角形板元素更加方便，所以 18 個自 由度的板元素會比 21 個自由度的板元素實用。文獻[17]最後的結論提到此 18 個自由度的高階板元素除了性能良好，此元素還有其他特點，其一：因 為有曲率當作節點自由度，因此可以直接計算出 internal moment。其二，有 曲率 、 、 當作節點自由度，在結構邊界條件上的設定可以更 加的滿足實際情況。 n w, 4 ,_n w 1 C 5 ,_n w w,_n6 xxi w, w,_xyi w,_yyi 雖然這種 18 個自由度 連續的高階板元素很早就被提出，也被證實是 一個性能優良的元素，但因為計算效率的問題，三角形板元素的發展著重 於 9 個自由度的板元素，因此在板殼的應用上鮮少有文獻使用。1990 年文 獻[18]使用文獻[19]的推導法，將 18 個自由度 連續的高階板元素剛度矩 1 C

陣拆成多個矩陣相乘，此推導法不需要使用數值積分，因此提昇了計算剛 度矩陣的效率。2004 年文獻[20]初次把此元素的形狀函數應用到 fusion MHD(magnetohydrodynamic)的問題上，並且在結論裡提到此 連續的元素 比 連續的元素更有效率、具有更多的優點。綜合此元素和 DKT 元素相 比的優點：擁有更高的精確度和收斂速度、 連續、具有定義的側向位移 場可以推導幾何剛度和質量矩陣。就本人所知，薄殼結構分析中還沒有文 獻使用此 18 個自由度的高階板元素和平面元素疊加形成的殼元素配合共旋 轉法分析幾何非線性和挫屈的問題。因此本文將使用文獻[18]推導的 18 個 自由度 連續的高階板元素和文獻[11]具有旋轉自由度的 QST 平面元素疊 加成一個殼元素，配合共旋轉法和推導準確的幾何剛度，試圖在分析薄殼 結構挫屈問題中能夠更加的精準偵測到平衡路徑上的分歧點和挫屈模態， 進一步模擬得到文獻[16]實驗中觀察到 d-cone 旋轉的現象。為了方便稱呼， 本文將 18 個自由度 連續的高階板元素簡稱為 RQT(reduced quintic triangular)元素 1 C 0 C 1 C 1 C 1 C 本文在第二章介紹本研究使用的平面殼元素以及決定變形參數的方 法。在第三章介紹本文的數值計算方法及程序。在第四章將使用線性及非 線性例題測試本文的平面三角形殼元素的性能，以及說明本文提出決定元 素變形參數的方法是可行的。

第二章

理論推導

本文將文獻[11]上具旋轉自由度的 QST 平面元素(quadratic strain triangle element)與 RQT 板元素(reduced quintic triangular element) [18]疊加成一個 3

節點 36 個自由度的平面三角殼元素，元素的節點自由度為 3 個位移、3 個 旋轉、3 個平面應變和 3 個位移二次微分。本章將描述本文中殼元素的基本 假設，並且介紹本文中所使用的座標系統，接下來分別推導平面元素、板 元素的剛度矩陣及幾何剛度以及節點參數在各座標系統之間的轉換關係， 最後提出一個決定元素節點變形參數的方法。 2.1 基本假設 本文中對非線性平面三角殼元素的推導，做以下假設： (1) 薄膜變形及彎曲變形之間無耦合作用，因此平面元素的變形和板元素的 變形可以各自描述。 (2)元素的變形位移及旋轉為小位移及小旋轉。 (3)在元素變形前，垂直於元素中心面的法向線段，在元素變形後，依然保 持直線，且沒有伸長或縮短及垂直於變形後的中心面。 2.2 座標系統 為了描述系統的運動及元素的變形，本文定義了三組座標系統：

(a) 固定總體座標系統(global coordinate system)： G i X



i1,2,3



結構體所有節點的座標、系統的邊界條件與其他座標系統的基底，以及結 構的平衡方程式，均在此座標系統中定義。在 座標系統內之座標值以 表示。 G i X



X,Y,Z



(b) 元素座標系統(element coordinate system)： E

i x



i1,2,3



此座標系統是建立在每一個元素變形後的最新位置上，其座標原點為元素 節點1， 軸為元素節點 1 與元素節點 2 在元素平面上的連線， 軸是在 元素平面上垂直於 軸，且朝著元素節點 3 的方向， 軸則是由 軸與 軸外積而得。元素的位移、元素變形、元素節點內力與元素剛度矩陣是在 此座標系統中定義，然後經由座標轉換，將其轉換至總體座標系統及基礎 座標系統。本文中 表示初始未變形時的元素座標， 表示第 E x₁ x₂E E x₁ E x₁ E i x 0 E x₃ x₂E E i I x I 個增量 迭代收斂後的元素座標， 表示當前變形位置的元素座標。在 座標系統 內之座標值以



表示。 E i x x_iE



z , y x,

(c) 節點基礎座標系統(base coordinate system)： B

ij x



i1,2,3





j1,2,3



此座標系統的原點是剛接在結構離散後的每一個節點，並與對應的節點一 起移動及旋轉。本文中節點 j 在其初始位置之 軸的方向為曲面在該節點 的法線方向， 軸、 軸的方向為曲面在該節點互相垂直的切線方向， 節點的應變、側向位移二次微分自由度是在此座標系統中定義。本文中 表示元素節點 B j x₃ B j x₁ x₂B_j B ij x 0 j 在初始未變形時的節點基礎座標，Ix_ijB表示元素節點

j

在第

I

個增量迭代收斂後的節點基礎座標，x_ijB表示元素節點 j 在當前變形位置的 節點基礎座標。 2.3 殼元素的變形描述 如圖2.1 所示，殼元素中心面上有三個節點，每個節點的自由度分別是 、 、 軸方向的位移 、 、 ( j=1,2,3)，繞 、 、 軸方 向的位移轉角 E x₁ x₂E x₃E u_j v_j w_j x₁E x₂E x₃E xj



、



_yj、



_zj (j 1,2,3)，應變分量



_xj、



_yj、



_xyj (j 1,2,3) 以及側向位移的二次微分w,_xxj、w,_xyj、w,_yyj (j 1,2,3)。 在圖 2.1 中的殼元素節點自由度u_j、v_j及



_zj (j1,2,3)是 QST 平面元 素[4]節點自由度，而w_j、



_xj、



_yj (j 1,2,3)為 RQT 板元素[16-18]的節點 自由度。本文以下敘述的元素變形、元素內力及元素剛度都是定義在元素 座標下。 2.3.1 QST 平面元素的變形描述 本文採用的三角平面元素其位移場為三次變化，應變場則為二次變化。 此元素有 3 個節點，每個節點有 6 個自由度，此元素可使用彼此間能互相 轉換的不同節點參數。本文中採用的節點參數為節點

j

在元素座標 、 軸 的位移分量 、 ，應變分量 E x₁ x₂E j u v_j



_xj、



_yj、



_xyj及逆時鐘方向的旋轉



_zj。但為 了方便推導，本文在元素推導時使用的節點參數是節點

j

在 、 軸的位 移 分 量 、 以 及 、 軸 的 位 移 分 量 分 別 對 面 積 座 標(area E x₁ x₂E j u v_j x₁E x₂E

coordinates)



、



的微分u,_j、u _j、 、 ，面積座標的介紹詳見附錄B。 3 N 6 N 2 N 0 N , _j u_  , 0 0 0 3 N 5 N , _j v_ j v,_ 7 N 0 } ,_j j v,_ 0 7 N 此元素的位移場可表示為： (2.1)  q  q 1  1 N N 0 j , Nt u u  (2.2) Nt v v  (2.3) } { _3  q q q  q_2 0 1 N j { N₂ u  N (2.4) } 0 0 0 0 _{5 8 9} 4 N N N 0 0 0 { N v  N (2.5) } 0 0 0 _{6 8 9} 4 N N (2.6) , { _j j u v u v   q 其中uu(





)、vv(



,



) v N 分別為在 、 軸的位移分量， 為對應於

u

的 形 狀 函 數 ， 為 對 應 於 的 形 狀 函 數 ， 為 節 點 參 數 向 量 ， 為面積座標的三次函數，其表示式詳見附錄 C，另外本文 中



表行矩陣。 E x₁ x₂E N_u

v

q_ i

N

i1,2,3 ,9)



( ,_ x 由小變形的假設，元素內任意點的正應變、剪應變及逆時鐘方向的剛體旋 轉可表示為： x u,



(2.7) y y v,



(2.8) x , y ,  xy u v



(2.9)

) , , ( 2 1 y x u v  



(2.10) 因

u

、

v

都是



、



的顯函數，所以 、 與 、 及 、 與 、 的關係可表示成(詳見附錄 C)： x u, u,_y u,_ u,_ v,_x v,_y v,_ v,_ (2.11)              y x x u u u u , , , ,    A (2.12)              y x x v v v v , , , ,    A (2.13)               31 31 21 21 , , , , y x y x y x y x x      A (2.14) j i ij x x x   (2.15) j i ij y y y   其中x_j、y_j分別是元素節點 在元素座標系統的座標值。

j

令節點參數向量q_x與q_m分別為： (2.16) } { _{x1 x2 x3} x q q q q  (2.17) } { _{m1 m2 m3} m q q q q  (2.18) } , , , , { _{j j xj yj xj yj} xj  u v u u v v q (2.19) } { _{j j xj yj xyj zj} mj  u v     q 其中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.11)式與(2.12) 式，可求得(2.3)式中q_與(2.16)式中q_x的關係： (2.20) x xq T q_  _

(2.21)            x x x x     R 0 0 0 R 0 0 0 R T (2.22)                               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.7)至(2.9)式可求得(2.16)式中q_x與(2.17)式中q_m的關係： (2.23) m xm x T q q  (2.24)            xm xm xm xm R 0 0 0 R 0 0 0 R T                     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 xm R (2.25) 將(2.20)式與(2.23)式代入(2.1)式與(2.2)式可得： (2.26) m t um m m t u uN T_ q N q (2.27) m t vm m m t v v N T_ q N q                        m m m xm x xm x xm x xm x m         R 0 0 0 R 0 0 0 R R R 0 0 0 R R 0 0 0 R R T T T (2.28)

                      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 xm x m R R R_{ } (2.29) 2.3.2 RQT 板元素的變形描述 本文採用的三角形板元素[16-18]有 3 個節點、每個節點有 6 個自由度。 節點自由度分別為節點 j 在元素座標下 軸的側向位移w_j、繞著 軸的旋 轉角 E x₃ x₁E xj



、繞著x₂E軸的旋轉角



_yj和側向位移二次微分 、 、 。文 獻[18]中推導的 RQT 元素自由度為 、 、 、 、 、 xx xxj w, w, w,_xy w yy w, xyj , j w w,_xj w,_yj w,_yyj ) (j 1,2,3 ，本文為了跟具有旋轉自由度的平面元素疊加，必須將 、 轉換為 xj w, w,_yj xj



、



_yj，因此本文先採用文獻[18]形狀函數的推導，其後轉換元素 自由度並且推導元素的內力及剛度矩陣。 文獻[18]中，假設元素的側向位移場是一個完整五次多項式，用面積座標表 示為： 2 4 3 9 1 4 3 8 3 4 2 7 1 4 2 6 3 4 1 5 2 4 1 4 5 3 3 5 2 2 5 1 1L L L L L L L L L L L L L L L w



































3 2 3 1 16 2 2 3 3 15 2 1 3 3 14 2 3 3 2 13 2 1 3 2 12 2 3 3 1 11 2 2 3 1 10L L



L L



L L



L L



L L



L L



L L L



       2 2 2 1 3 21 2 3 2 1 2 20 2 3 2 2 1 19 2 1 3 3 18 3 1 3 2 17L L L



L L L



LL L



L L L



L L L



     (2.30) 其中



₁~



₂₁為未定係數。 18 1 ~





可以很輕易的用 18 個節點自由度表示，再利用位移在邊上法線方

向 微 分w,_n為 三 次 函 數 的 三 個 限 制 條 件 ， 可 以 把



₁₉ ~



₂₁各 自 表 示 成 18 1 ~





的函數，邊上的切線和法線方向示意如圖2.2。最後把用



₁~



₂₁代 入(2.30)式，整理各節點自由度的係數，即可得到形狀函數N_y1 ~ N_y18。 側向位移場可以由形狀函數和節點自由度表示成[18]： y q 1 y 1 y q j t y wN t y  N N t y q q  t yj  w q y q (2.31)



Ny2 Ny3 ... Ny18





(2.32)



y2 qy3 (2.33)



w,xj w,yj w,xxj w,xyj w,yyj



(2.34) 其中 為節點自由度向量 根據克西荷夫的板理論，元素的位移場及旋轉場可以表示為為： x w  z u  (2.35) y w  z v  (2.36) y w x _  



(2.37) x w   b ybq yb 0 0 R (2.38) y 



由(2.37)和(2.38)式可以推導出不同節點自由度q_b、q 的轉換矩陣_y T_by： y T q  (2.39)            yb yb yb R 0 0 R 0 0 T (2.40)

                    1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 yb R (2.41) 其中



_{b1 b2 b3} b q q q q 



(2.42)



j xj yj xx xy yy



bj  w





w, w, w, q (2.43) 將(2.39)式代入(2.31)式推得： b t b b yb t y y t y wN q N T q N q (2.44) y t yb b T N N  (2.45) 其中N_b為本文RQT 板元素的形狀函數，其表示式詳見附錄 C 假設板的變形為小變形，其應變可以表示成： x u x   



(2.46) y v y   



(2.47) x v y u xy      



(2.48) 將(2.35)、(2.36)式是代入(2.46)至(2.48)式推得： 2 2 x w z x    



(2.49) 2 2 y w z y    



(2.50)

) 2 ( ) ( 2 2 2 y x w z y x w y x w z xy _{ }            



(2.51) 將(2.44)式代入(2.49)至(2.51)式可以將應變寫成形狀函數與節點自由度的函 數： b b b t b t b t b xy y x b z y x y x z q B q N N N ε                                       2 2 2 2 2 2







(2.52) 其中B 為應變-位移的關係矩陣： 18 3 2 2 2 2 2 2                           y x y x t b t b t b b N N N B (2.53) 2.4 元素內力與元素剛度矩陣 本文中殼元素的節點內力是由QST 平面元素及 RQT 板元素的節點內力 、 組合而成，元素剛度矩陣 m f f_b k 是由 QST 平面元素剛度矩陣 、RQT 板元素剛度矩陣 所疊加而成。本節中將用虛功原理推導QST 平面元素及 RQT 板元素的節點內力及剛度矩陣。 m k b k 在平面應力狀態，等向性線彈性材料的應變與應力關係為： ε E σ  (2.54)

} , , {



_x



_y



_xy  σ (2.55) } , , {



_x



_y



_xy  ε (2.56)                2 1 0 0 0 1 0 1 1 2









E E (2.57) yz yz G





 ,



_xz G



_xz (2.58)

其中E 是楊氏模數(Young’s module)，



是蒲松比(Poisson ratio)。

2.4.1 QST 平面元素之節點內力與剛度矩陣 將(2.26)式、(2.27)式代入(2.7)至(2.8)式可得： (2.59)  q BT ε_m  (2.60) } { _{x y xy} m 







ε (2.61)             x t v y t u y t v x t u , , , , N N N N B 其中 t _x、 、 、 的表示式詳見附錄B。 u

,

N

t _y u

,

N

t _x v

,

N

t _y v

,

N

將(2.59)式代入(2.54)式可得： (2.62) m m m EBT q σ  _ 由虛功原理可得： (2.63)



 V m t m m t mf ε σ dV q





} (2.64) { _{m1 m2 m3} m f f f f 

(2.65) } { _{xj yj xj yj j zj} mj  f f m m m m f 其中 為對應於 的元素節點內力， 、 分別為對應於 、 的節點 力， 、 、 分別為對應於 m f xj m_ m q j m_ xj f f_yj

u

_j

v

_j yj m_



_xj、



_yj、



_xyj的廣義節點力矩， 為對 應於 zj m zj



的節點傳統力矩，V 為元素體積。 將(2.59)式、(2.62)式代入(2.63)式可得 (2.66) m m V t t m t m m t mf q T B EBdVT q q



_{ }





_

由(2.66)式可得： (2.67) m m m k q f  (2.68) m T m m T kT k  (2.69)



 V T dV EB B k_ 其中k_、k_m分別為應於q_、q_m的QST 元素剛度矩陣。 2.4.2 RQT 板元素的節點內力及剛度矩陣 將(2.52)式代入(2.54)式可得： b t by b zEBT q σ  (2.70) 由虛功原理



W_I 



W_E可得： dV v b t b b t bf 



ε σ q





(2.71) } { _{b1 b2 b3} b f f f f  (2.72) } { _{zj xj yj xj xyj yj} bj  f m m b b b f (2.73)

其中 為對應於 的元素節點內力， 為對應於 的節點力， 、 為 對應於 b f q_b f_zj w_j m_xj m_yj xj



、



_yj的節點力矩， 、b 、 為對應於 、 、 w 的節 點雙力矩， xj b _xyj b_yj w,_xx w,_xy ,_yy V 為元素體積。 將(2.52)、(2.70)式代入(2.71)式推得： dV z V b t by t by t b b t bf  q



T B EBT q q



2



(2.74) b t by A t by b T B DBdAT q f 



(2.75)                

_

_ 2 1 0 0 0 1 0 1 ) 1 ( 12 2 3 2 2 2









Eh dz z h h E D (2.76) h 為元素厚度 b b b k q f  (2.77) t bx A t bx b T B DBdAT k 



(2.78) 其中k_b為RQT 板元素線性剛度矩陣。 2.4.3 元素幾何剛度矩陣 本文在幾何非線性分析中，切線剛度為線性剛度加上幾何剛度。文獻[26]

中以完整的Green’s strain 及 total Lagrangian 推導法推導出三維元素之切線

剛度矩陣的通式。文獻[6]修改、簡化文獻[26]推導元素幾何剛度矩陣的方 法，使其適用於推導殼元素的幾何剛度。本文採用文獻[6]推導殼元素幾何 剛度的方法，元素幾何剛度矩陣可表示式成

dV V t



 G MG k_ (2.79) } {G₁ G₂ G₃ G  I (2.80)                      y wj y vj y uj x wj x vj x uj j , , , , , , N 0 0 N 0 N N 0 0 N 0 N G        I I I I M y xy xy x









為 3×3 單位矩陣 (2.81) b m σ σ σ   σ {



_x ,



_y ,



_xy} (2.82) 其中 j = 1,2,3 為第 j 個節點，I 為 3×3 的單位矩陣、0 為 1×6 的零矩陣， 、 為(2.66)式、(2.74)式之 QST 平面元素及 RQT 板元素的應力場。 m σ b σ 為了簡化幾何剛度的計算，本文中 、 、 皆為線性位移場，其線 性位移場與 CST 元素[12]的形狀函數相同，因此本文在非線性迭代中使用 的切線剛度為近似切線剛度。 uj N N_vj N_wj 2.5 節點參數在各座標系統之間的轉換 為了建立結構的平衡方程式，必須將元素節點參數q 中的 、u_j v 、_j w_j和 xj



、



_yj、



_zj轉換成固定總體座標的分量，



_xj、



_yj、



_xyj、w 、 、 轉換成對應於節點基礎座標的分量，才能將各元素節點內力組合成結構系 統節點內力以及將元素剛度矩陣疊加成結構系統剛度矩陣。在元素座標的 節點位移 、 、 和節點旋轉 xx , w,_xy w,_yy j u v_j w_j



_xj、



_yj、



_zj與在固定總體座標的節點位

移和節點旋轉之關係式為： (2.83)                      j j j t GE j j j W V U w v u A                      G Zj G Yj G Xj t GE zj yj xj













A (2.84) (2.85) ] [ ₁E ₂E ₃E GE e e e A  其中 、 、 分別為在 、 、 軸的位移分量， 、 、 分 別為在 、X 、 軸的旋轉分量，e 為元素座標 方 向的單位向量。 j U X j V G 1 j W G 2 G X₁ X₂G X₃G (i  E i G Xj



x G Yj



 i G Yj





2,3



G X₃ 1,2,3) _iE 1, 其中節點應變



_xj、



_yj、



_xyj和 、 、 在元素座標與節點基礎座 標之間的關係式為： xx w, w,_xy w,_yy B j EBj j R ε ε   (2.86)



xj yj xyj



j     ε (2.87) (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 B j EBj j R κ κ   (2.90)



xxj xyj yyj



j  w, w, w, κ (2.91) (2.92)



B yyj B xyj B xxj B j  w, w, w, κ



(2.93)               2 2 2 2 2 2 2 2 Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj Bj EBj c s c s s c s c s c s s c c  R 其中 、 分別為 在元素座標 、 軸的分量， 、 、 、 、 、 為對應於節點基礎座標 、 軸的分量， Bj c w Bj s B yyj B j 1 e _xE 1 B j x₁ E x₂ B x₂ B xj



e B yj



, ( B xyj



2 , 1  B xxj w, 3 B xyj w, , _{j ij}B i j , )為節 點基礎座標x_ijB



i, j1,2,3



方向的單位向量。 B q 為位移、旋轉定義在總體座標，應變、側向位移二次微分自由度定義在 基礎座標的節點自由度向量。 為位移、旋轉、應變、側向位移二次微分 自由度皆定義在元素座標的節點自由度向量，其兩者轉換關係 為： E q EB T (2.94)



B1 B2 B3



B q q q q  B xyj B yj B xj G Zj G Yj G Xj j j j Bj {U V W













q , , ,B } (2.95) yyj B xyj B xxj w w w q_E T_EBq_B (2.96) q_Ej {u_j v_j w_j



_Xj



_Yj



_Zj



_xj



_yj



_xyj w,_xxj w,_xyj w,_yyj} (2.97)            3 2 1 EB EB EB EB T 0 0 0 T 0 0 0 T T (2.98)                  EBj EBj t GE t GE EBj R 0 0 0 0 R 0 0 0 0 A 0 0 0 0 A T (2.99)

殼元素在元素座標下的節點內力可以表示成： b m E f f f   (2.100) } { _{E1 E2 E3} E f f f f  (2.101) j yj xj zj yj xj zj yj xj Ej {f f f m m m m m m f b_xj b_xyj b_yj} (2.102) 由反梯度法則(contragradient law)及(2.96)式可得： (2.103) E t EB B T f f  } (2.104) { _{B1 B2 B3} B f f f f  B j B yj B xj G Zj G Yj G Xj Zj Yj Xj Bj {F F F m m m m m m f B } (2.105) yj B xyj B xj b b b_{  } 其中 為對應於 的元素節點內力， 、 、 分別為對應於 、 、 的節點力，m 、 、 為對應於 、 、 的節點傳統力矩， 、 、 分別為對應於 、 、 的廣義節點力矩， 、 b 、 分 別為對應於 、 、w 的廣義節點力矩。將(2.96)式、(2.100)式代入 (2.103)式可得： B f q_B G Xj w Xj F G Xj



B xyj Yj F G Yj



Zj F G Zj



j U B yj  j V B xj m_ B j b_ j W B yj m_ G Yj m B xyj , G Zj m B xj



B , B j m_



_yjB yyj



B xj b_ B xxj , w (2.106) B B B k q f  (2.107) EB E t EB B T k T k  其中 是使用直接剛度法將 和 結合的殼元素在元素座標下的線性剛 度。 E k k_m k_b

2.6 元素節點變形參數的決定方法 本文中採用增量迭代法解非線性平衡方程式，假設第 個位置已知，此 處的第 個位置，是指第 個增量迭代收斂後的平衡位置。 為元素節點

I

I

) 3 , 2

I

x_j

j

, 1 (j  在當前元素座標上的位置向量，I _j、 、 、I 、 _j、 、 以及 分別為元素節點 在固定總體座標中第 個位置的位置向 量、元素節點 在第 個位置的節點基礎座標、元素節點 在第 個位置的 應變向量、元素節點 在第 個位置的側向位移二次微分向量、元素節點 的 位移增量向量、元素節點 的應變增量向量、元素節點 的側向位移二次微 分增量向量以及元素節點 的旋轉增量。本文中假設元素節點 受 X B ij I_x B j I_ε

I

j

B j



U

I

j

B j ε Δ

j

j U B j κ Δ _j

j

j

I

j

j

I

j

j

 、 、 以及 B j ε Δ B j κ Δ _j作用後的變形過程如下： 1. 節點 受 的作用由第 個位置 平移到節點當前變形位置 ，在 移動過程中，節點 無剛體旋轉，即元素各邊在節點的切線方向維持不 變，剛接在其上的節點基礎座標的方位亦不變。

j

U_j

I

IX_j X_j

j

2. 節點 及剛接在其上的節點基礎座標軸

j

Ix_ijB受到_j作用，旋轉到當前變 形位置的 B。 ij x 3. 節點 在當前變形位置的節點基礎座標上之應變分量增加一個增量 。

j

B j ε Δ 4. 節點 在當前變形位置的節點基礎座標上之側向位移二次微分分量增加 一個增量 。

j

B j κ Δ

由上述的變形過程可知節點 當前變形位置的固定總體座標 可以由 加上 得到，由 可以利用元素座標定義求出當前變形位置的元素座標 。如圖2.3 所示，將元素變形後當前變形位置的元素座標 與初始未變 形時的元素座標0 重疊，並且將當前變形位置的元素座標 旋轉一個角 度

j

X_j E i x E i x j I_X j U  X_j E i x E i x



，



決定方式詳見附錄D。元素節點變形參數可決定如下： (1) 節點位移u_j 節點位移 可由元素節點 在當前變形位置的元素座標上的位置向量 與 元素節點 在初始未變形時的元素座標上位置向量的 的差得出： j u

j

j

x

_j j

x

0 u_j x_j0x_j (2.108) (2.109)



j j j j  u v w u





j j j



j  x y z x (2.110) (2.111)



j j j j x y z 0 0 0 0_{x }



由元素座標定義的方式可知0z_j  z_j w_j 0。 (2) 節點應變 ε_j 利用第 個位置在節點基礎座標上的節點應變向量 加上增量應變向量 可以求出當前變形位置的節點基礎座標上的節點應變向量：

I

IεB_j B j ε Δ B j B j I B j ε ε ε  Δ (2.112) (2.113)



B



xyj B yj B xj B j Δ Δ Δ Δε  其中 B為當前變形位置的基礎節點座標上的節點應變向量。 j

ε

(3) 側向位移二次微分κ_j 利用第 個位置在節點基礎座標上的節點應變向量 加上增量應變向量 可以求出當前變形位置的節點基礎座標上的節點應變向量：

I

IκB_j B j κ Δ B j B j I B j κ κ κ  Δ (2.114) (2.115)



B yyj B xyj B xxj B j Δw, Δw, Δw, Δκ 



(3) 節點變形角



_ij 如圖 2.4 所示，將旋轉向量θ_nj作用在 軸上，使其旋轉到 軸的新位 置， 可表示為： B ij x x_ijB 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.116) 其中 及 分別為當前變形位置的元素座標之 軸的單位向量，及當前 變形位置的節點基礎座標之 軸的單位向量。此時， 軸與 軸重合， 且 軸和 軸與 軸垂直，與 軸、 軸共平面，但不重合。 E 3 e B j  B j 3 e B j x₂ E x₃ B j x₃ x₃B_j x₃E x₁ x₃E x₁E x₂E 同理，在初始未變形時將旋轉向量0θ_nj作用在 B軸上，使其旋轉到 ij x 0 B ij x  0 _軸 的新位置，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.117) 其中 、 分別為初始未變形時元素座標之 軸的單位向量，以及初 始未變形時基礎節點座標之 軸的單位向量。若結構初始未變形時為曲 面，則 軸和 軸不平行，所以 ；若結構初始未變形時為平面， E 3 0_e x 0 B j 3 0_{e x}E 3 0 B j x₃ 0 E 3 0x3Bj 0θnj 0

則 _xE軸和 軸平行，所以 。 3 0 B j x₃ 0 _{θ 0} nj 0 2,3) 如圖 2.5 所示，



_zje₃E(j 1 , 為將0x₁B_j軸旋轉到x₁B_j軸的旋轉向量，



_zj可 表示為： E B j B j zj 1 1 3 0 ₎ ( e e e e e e e      _{ } B j B j B j B j 1 1 0 1 1 0 1 sin      B j  1 e



0_e (2.118) 其中 ₁B_j、 分別為初始未變形時之0x₁B_j軸的單位向量，以及當前變形位 置之_xB 軸的單位向量。 1 E nj θ 0 j  令 及 為將 與 分別以初始未變形時的元素座標及當前變形位 置的元素座標為基底，節點變形角 E nj θ 0θ_nj θ_nj ij



(i x,y) (j 1,2,3)可表示為： E nj yj xj θ 0 0       E nj θ      





(2.119) 2.7 旋轉向量 本文中使用旋轉向量來表示一個有限旋轉，如圖2.6 所示，一向量R受到一 旋轉向量



n的作用而轉到一個新位置R，R 與R之間的關係可表示成[39] cos R R



(1-cos



)(nR)nsin



(nR) (2.120) 2.8 系統平衡方程式與收斂準則 結構系統受外力負荷時，其平衡方程式可表示為



Q



 F

 

Q

 P



0

Ψ

,





(2.121) 其中Ψ 為系統不平衡力向量，系統節點內力 F 可由(2.114)式的元素節點內力

B

f 疊加得出，Q 為系統位移向量，



為負荷參數，

P

為參考外力負荷向量。

若外力為與變形位置相關(configuration dependent)的外力，則 在每一個變

形位置都須更新。本文以不平衡力 的weighted Euclidean norm 做為平衡迭

代時的誤差度量，且收歛準則表示為

P

Ψ tol e N e Ψ 



P (2.122) 其中

N

代表離散系統的自由度總數，e_tol為一給定的容許誤差值。

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

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

I

Q_I、負荷參數為



_I， 則第

I



1

個增量的初始增量位移向量 ，可以利用尤拉預測值(Euler Q predictor)求得[28] T r Q 



 (3.1) (3.2)

 

K P r  I 1 T T 其中 為初始增量負荷參數，



為第 個平衡位置的系統切線剛度矩陣。 (3.1)式中的 I T K

I



 可利用下式求出[27]

 

1/2 T t T l r r    



(3.3) 其中正負符號決定方法為，當第

I



1

與 個增量收歛時，其系統切線剛度 矩陣之行列式值同號，則

I



 的正負符號和第 個增量時相同；若異號則符 號相反。

I

l  表示第

I



1

個增量的增量弧長，其值可以如下決定[27]





12 I D I J J l l   (3.4) 其中J_D為給定的期望迭代次數，J_I為第 個增量迭代至平衡所使用的迭代

I