平面三角形薄殼元素之共旋轉推導法

(1)

國

立

交

通

大

學

機械工程學系碩士班

碩士論文

平面三角形薄殼元素之共旋轉推導法

Co-rotational formulation for the triangular thin

flat shell element

研究生：沈佳鴻

指導教授：蕭國模博士

(2)

平面三角形薄殼元素之共旋轉推導法

Co-rotational formulation for the triangular thin

flat shell element

研究生：沈佳鴻 Student： Jia-Hung Shen

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

Hsinchu, Taiwan, Republic of China

(3)

平面三角形薄殼元素之共旋轉推導法

研究生：沈佳鴻指導教授：蕭國模博士

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

摘要

本研究主要目的是以共旋轉全拉格朗日推導法(Co-rotational total Lagrangian formulation)、von Karman 板正確的變形機制、一致性二階線性

化(Consistent second order linearization)及虛功原理，推導一個具面內旋轉自由度平面三角殼元素，並將其應用在薄殼結構的幾何非線性及挫屈分析。本文中推導的平面三角殼元素有 3 節點、每個節點有 9 個自由度，元素的節點自由度為節點位移向量的 3 個分量、節點旋轉向量的 3 個分量及節點平面應變的 3 個分量。本研究將元素的三個節點建立在板的中平面上，在三個節點的當前位置建立一個元素座標，並用中平面變形後的位移及法向量描述殼元素的變形。本研究利用極分解定理(Polar decomposition theorem)將殼中平面的變形梯度(Deformation gradient)分解成一個旋轉矩陣

(Rotation matrix)和一個伸縮矩陣(Stretch matrix)的乘積，並用一剛接在元素

中平面的座標系統的旋轉表示變形梯度中的剛體旋轉。本研究用三個旋轉參數來描述該中平面座標系統的旋轉。本文採用牛頓-拉福森(Newton-Raphson)法和弧長控制(Arc-length control)法的增量疊代法來解結構的非線性平衡方程式。本研究分析文獻上常見的殼結構基準問題，並與文獻上的線性解、非線性解、挫屈負荷比較。本研究探討元素切線剛度矩陣中一些高次項對結構之非線性行為及挫屈負荷的影響。

(4)

Co-rotational formulation for the triangular thin flat shell element

Student：Jia-Hung Shen Advisor：Dr. Kuo-Mo Hsiao

Department of Mechanical Engineering

National Chiao Tung University

ABSTRACT

A facet triangular shell element with drilling degree of freedom is developed by using a corotational total Lagrangian formulation for the geometrically nonlinear analysis of thin shell structure with large rotation but small strain.

The element developed has three nodes with nine degrees of freedom per node. The element nodes are chosen to be located at the mid-plane of the plate element. The deformations of the shell element are described in a current element coordinate system constructed at the current configuration of the shell element. The element nodal forces are derived using the virtual work principle, the exact kinematics of the von Karman plate, and the consistent second order linearization. The element tangent stiffness matrix may be obtained by differencing the element nodal force with respect to nodal parameters. The deformation of the shell element is determined by the displacements of the mid-plane and the rotations of a mid-plane coordinate system associated with

(5)

each point of the mid-plane relative to the current element coordinate system. The origin of the mid-plane coordinates is rigidly tied to mid-plane. Three rotation parameters are defined to describe the rotation of the mid-plane coordinate system. For convenience, two set of nodal parameters are employed to determine the displacement fields of the element. The first set of nodal parameters is chosen to be three nodal displacements, three nodal rotation parameters, and three strains. The second set of nodal parameters is chosen to be three nodal displacement and six nodal values of the first spatial derivative of displacements. To determine the relationship between these two sets of nodal parameters, the deformation gradient at each element node is decomposed into the product of a proper orthogonal matrix and a right stretch matrix by using the polar decomposition theorem, and the rotation matrix corresponding to the rotation of the mid-plane coordinate system relative to the current element coordinate system is regarded as the proper orthogonal matrix. Two sets of nodal parameters are used for the assembly of the structural equations.

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. The zero value of the tangent stiffness determinant of the structure is used as the criterion of the buckling state. Benchmark problems for linear and geometric nonlinear analysis of shells given in the literature are studied to demonstrate the accuracy and efficiency of the proposed shell element. The effect of the first order terms of the transformation matrix between the variation of the two sets of nodal parameters on the equilibrium path and buckling load are also investigated.

(6)

致謝

衷心感謝指導教授蕭國模博士在這兩年碩士班期間的指導與教誨，使本論文得以順利的完成。老師在研究上嚴謹的態度及日常生活上的關懷，使我受益良多，在此致上最高的謝意與敬意。感謝尹慶中老師和鄭文雅老師撥冗擔任口試委員並對本論文所提出的指正與建議，使本論文能夠更臻完善。感謝黃楚璋學長、許彤羽學姊的照顧，同學林琮棋、莊士緯以及學弟林群禮、高嘉鴻在學業上的砥礪與成長。感謝父母親、哥哥、弟弟、女友等關心我的朋友對我的支持與鼓勵，僅以此成果與榮耀，獻給我親愛的父母以及所有關心我的人。

(7)

目錄

中文摘要 ... I 英文摘要 ... II 致謝 ... IV 目錄 ... V 表目錄 ... VIII 圖目錄 ... IX 第一章緒論 ... 1 第二章理論推導... 9 2.1 基本假設 ... 9 2.2 座標系統 ... 10 2.3 旋轉向量 ... 11 2.4 殼的變形描述... 12 2.5 元素節點參數及節點力向量... 18 2.6 元素的應變及其變分... 23 2.7 元素節點內力之推導... 28 2.8 元素剛度矩陣... 30 2.9 殼元素的位移場... 31 2.10 座標系統轉換... 34 2.11 系統平衡方程式與收斂準則... 37 2.12 元素節點變形參數的決定方法... 38 第三章數值計算方法與程序... 42 3.1 增量迭代法 ... 42 3.2 弧長控制法 ... 44

(8)

3.3 二分法 ... 45 3.4 數值程序 ... 46 第四章數值分析與結果... 49 4.1 半圓環受單點集中力作用... 50 4.2 直角構架受到端點剪力作用... 51 4.3 圓柱殼片段受到單點集中力作用... 51 4.4 槽型斷面梁之側向扭轉挫屈... 52 4.5 受壓之簡支承板... 54 4.6 T 型斷面梁受側向負荷 ... 55 4.7 槽型斷面梁受扭矩... 55 4.8 直角梁受到單點集中力作用... 56 4.9 懸臂圓柱殼受到單點集中力作用... 57 4.10 半球殼受到單點集中力作用... 58 4.11 圓柱殼受一對集中力... 59 4.12 圓柱殼受四個徑向集中力作用... 60 4.13 開口型半球殼受徑向集中力作用... 61 4.14 裂縫環形板受均勻力負荷... 62 4.15 聚酯圓柱薄殼受兩階段負荷作用... 63 第五章結論與展望... 65 5.1 結論 ... 65 5.2 未來研究方向... 66 參考文獻 ... 67 附錄 A 元素座標系統的決定 ... 113 附錄 B 切線剛度矩陣k ... 115 附錄 C 面積座標 ... 119

(9)

附錄 D QST 元素的形狀函數及其微分... 122

附錄 E DKT 元素的形狀函數 ... 123

附錄 F Projector matrix ... 127

(10)

表目錄

表 4.1 圓柱殼片受單點集中力作用的挫屈負荷 (例題 4.3，Mesh 10×10) ... 74 表 4.2 槽型斷面梁的挫屈負荷 (例題 4.4) ... 74 表 4.3 受壓之簡支承板的挫屈負荷 (例題 4.5) ... 75 表 4.4 T 型斷面梁的挫屈負荷 (例題 4.6) ... 76 表 4.5 槽型斷面梁的極限點 (例題 4.7) ... 77 表 4.6 圓柱殼受四個徑向集中力作用的挫屈負荷 (例題 4.12) ... 77 表 4.7 開口型半球殼 B 點位移之線性解 (例題 4.13，P = 1)... 78 表 4.8 開口型半球殼 B 點位移的非線性解 (例題 4.13，P = 400)... 78 表 4.9 聚酯圓柱薄殼受兩階段負荷作用的挫屈負荷 (例題 4.15) ... 79

(11)

圖目錄

圖 1.1 文獻[19]實驗所觀察到四種變形轉換(a-d)及 對應於a-c 圖結構的上視圖(e-g)... 80 圖 2.1 旋轉向量 ... 81 圖 2.2 薄殼中 P、Q 點之位移以及元素座標與中平面座標 的關係圖 ... 82 圖 2.3 元素節點 j 中心面之 B軸受旋轉向量 ij x θ_nj作用的示意圖 ... 83 圖 2.4 元素節點 j 中心面之0x₁B_j軸受旋轉向量



₃_{j 3}eE作用的示意圖 ... 84 圖 4.1 半圓環受到單點集中力作用 (a)結構尺寸示意圖 (b)網格 18×1 示意圖 (c)網格 18×2 示意圖 ... 85 圖 4.2 半圓環受到單點集中力作用之負荷－位移曲線圖... 86 圖 4.3 直角構架受到端點剪力作用 (a)結構尺寸示意圖 (b)網格 M21 與網格 M22 示意圖 ... 87 圖 4.4 直角構架受到端點剪力作用之負荷－位移曲線圖... 88 圖 4.5 圓柱殼片段受到單點集中力作用 (a)結構尺寸示意圖 (b)力負荷圖 (c)網格 10×10 示意圖 ... 89 圖 4.6 圓柱殼片段受到單點集中力作用之負荷－位移曲線圖... 90 圖 4.7 槽型斷面梁之側向扭轉挫屈(例題 4.4) (a)結構尺寸示意圖 (b)網格(1+2+1)×2 示意圖... 91 圖 4.8 受壓之簡支承板 (a)結構尺寸示意圖 (b)網格 4×6 示意圖 ... 92

(12)

圖 4.9 T 型斷面梁 (a)結構尺寸示意圖 (b)網格(2+2+3)×4 示意圖... 93 圖 4.10 槽型斷面梁(例題 4.7) (a)結構尺寸示意圖 (b)網格(1+2+1)×2 示意圖... 94 圖 4.11 槽型斷面梁之負荷－位移曲線圖(例題 4.7) ... 95 圖 4.12 直角梁受到單點集中力作用(例題 4.8) (a)結構尺寸示意圖 (b)網格(2+3)×2 示意圖 ... 96 圖 4.13 直角梁受到單點集中力作用之負荷－位移曲線圖 (例題 4.8) ... 97 圖 4.14 懸臂圓柱殼受到單點集中力作用(例題 4.9) (a)結構尺寸示意圖 (b)網格 16×16 示意圖 ... 98 圖 4.15 懸臂圓柱殼受到單點集中力作用之負荷－位移曲線圖 (例題 4.9) ... 99 圖 4.16 半球殼受到單點集中力作用(例題 4.10) (a)結構尺寸示意圖 (b)網格 12×12 示意圖 ... 100 圖 4.17 半球殼受到單點集中力作用之負荷－位移曲線圖 (例題 4.10) ... 101 圖 4.18 圓柱殼受單點集中力作用(例題 4.11) (a)結構尺寸示意圖 (b)網格 16×24 示意圖 ... 102 圖 4.19 圓柱殼受一對集中力作用之負荷－位移曲線圖 (例題 4.11) ... 103 圖 4.20 圓柱殼受四個徑向單點集中力作用(例題 4.12) (a)結構尺寸示意圖 (b)網格 8×48 示意圖 ... 104 圖 4.21 圓柱殼受四個徑向單點集中力作用之負荷－位移曲線圖 (例題 4.12) ... 105

(13)

圖 4.22 開口型半球殼受集中力作用(例題 4.13) (a)結構尺寸示意圖 (b)網格 12×12 示意圖 ... 106 圖 4.23 開口型半球殼受集中力作用之負荷－位移曲線圖 (例題 4.13) ... 107 圖 4.24 裂縫環形板受均勻力負荷作用(例題 4.14) (a)結構尺寸示意圖 (b)網格 6×30 示意圖 ... 108 圖 4.25 裂縫環形板受均勻力負荷作用之負荷－位移曲線圖 (例題 4.14) ... 109 圖 4.26 圓柱薄殼(例題 4.15) (a)結構示意圖 (b)俯視圖 (c)前視圖... 110 圖 4.27 圓柱薄殼第一階段 E 點之位移－負荷參數曲線圖 (例題 4.15) ... 111 圖 4.28 圓柱薄殼第二階段 E 點之反力－負荷參數曲線圖 (例題 4.15) ... 112

(14)

第一章緒論

在近代工程設計的發展上，對材料的要求與結構的表現趨向於高強度與輕量化，舉凡建築結構、航太設備、壓力容器、軍事載體及汽車工業等，設計者要考量以最小成本來達到所需的機能與強度，並且兼顧外型的美觀。由於薄殼在承受彎曲應力與拉伸應力的表現上有良好的表現，又能達到經濟與輕量化的要求，因此薄殼為設計者最常使用的結構之ㄧ，而廣泛應用在工程及生活上。常見的薄殼結構包括：建築屋頂、飛機蒙皮、液體儲存槽、人造衛星、火箭、船體結構、水中潛體等。薄殼結構受到外力作用經常會產生大位移和大旋轉，在大位移和大旋轉的問題中，位移和外力往往不是線性關係，因此需要使用非線性分析的方法來探討由幾何形狀改變所造成的非線性行為。常見的幾何非線性分析的推導方法有三種：全拉格朗日法(total

Lagrangian formulation)、更新拉格朗日法(updated Lagrangian formulation)和

共旋轉法(Co-rotational formulation)。全拉格朗日法是用初始狀態為參考位置來表示總位移和旋轉；更新拉格朗日法是以結構上一個平衡狀態為參考位置來表示增量位移和旋轉；共旋轉法是利用建立在元素當前變形位置的元素座標將剛體位移及旋轉從總位移及旋轉中扣除，剩下的位移和旋轉即為小位移和小旋轉，因此若使用共旋轉法，原本在線性分析的元素也可以

(15)

應用在大位移、小應變的幾何非線性分析，共旋轉法在梁與殼結構的幾何非線性分析已經被廣泛的使用[1-10]。分析殼結構最常用的方法為有限元素法。殼元素大致分為三類：平面殼元素、曲面殼元素和等參數元素。常見的平面殼元素是由一個平面板元素和一個平面元素疊加而成，此種元素的推導方式簡單，而且在數值計算上比曲面殼元素更有效率[11]，已分別應用在殼結構的線性[12-15]和幾何非線性分析[4-11,16]。平面殼元素最常使用的形狀為三角形或四邊形，對任何不規則形狀的殼結構，我們都可以輕易地將其切割成有限的三角形組合，但不一定適合將其切割成四邊形的組合，故三角形元素在文獻上被廣泛的探討及使用。

CST(Constant strain triangle)元素和 LST(Linear strain triangle)元素都是

最簡單的平面元素，在薄殼分析上常常使用它們與合適的板元素疊加，因為這種平面三角殼元素缺少旋轉自由度(Drilling degree of freedom)，所以其元素剛度的面內旋轉剛度(In-plane rotational stiffness)為零，為了避免系統剛度矩陣因奇異性(Singularity)造成分析的困難，常見的解決方法有兩種：(1) 加上一個人工的面內旋轉自由度[17-20]。(2) 採用具有旋轉自由度的平面元素[8,9,17,21-23]。文獻[23]提到三角平面元素加入旋轉自由度的優點為：改善三角平面元素的性能並避免使用到三角形邊上的節點，因為邊上的節點會影響到網格生成，而且在模擬非線性分析與動態分析時較為困難；當

(16)

三角平面元素與三角板元素疊加時，能滿足物理上一個節點有 3 個旋轉自由度的要求；三角形元素與殼元素、板元素或是梁元素同時使用時，能使接合簡單化。1964 至 1983 年期間，許多人在研究如何在 3 節點的三角形平面元素上加入節點旋轉自由度，希望能得到一個 3 節點 9 個自由度且具節點旋轉自由度的三角形平面元素，但是都沒有得出可用的元素。1984 年文獻[21]Allman 提出第一個成功帶有旋轉自由度、3 個節點 9 個自由度的三角形平面元素。2002 年文獻[22]提出 DLST 元素是一個具有 12 個自由度，三角形頂點自由度為 2 個位移、1 個旋轉並且邊上中點自由度為 1 個位移。2003 年文獻[23]提出的 OPT 元素具有 3 個節點，每個節點 2 個位移、1 個旋轉的自由度。2008 年文獻[17]成功的將一個具有旋轉自由度的 QST 元素[23]應用在平面應變問題的共旋轉法幾何非線性分析上，該元素具有 3 個節點、 18 個自由度，每個節點有 2 個位移、1 個旋轉及 3 個應變自由度。從 1960 年起有許多的三角形板元素被提出[24-28]，其中具有 3 個節點，每個節點有 1 個位移、兩個旋轉自由度，總共 9 個自由度之三角形板元素的研究發展最為迅速。文獻[18]中比較 DKT、HSM、BCIZ、HCT 等 9 個自由度三角形板元素的線性分析和振動分析後，認為 DKT 元素是這些三角形板元素中最有效率的元素，在靜態和模態分析中均可以收斂到準確的答案。文獻[23]中指出 DKT 元素內部不能滿足 x y y x     



的連續條件且沒有

(17)

定義側向位移場，因為計算板元素之質量矩陣及幾何剛度矩陣需要該元素的側向位移場，所以 DKT 元素的質量矩陣和幾何剛度矩陣皆使用其他元素的側向位移場來推導。1968 年文獻[20]把位移假設為完整五次多項式的元素稱為 TUBA 6，除了在三頂點有 1 個側向位移、2 個側向位移一次微分、，3 個側向位移兩次微分、、之外，還有 3 個邊上中點側向位移對邊上法線方向微分，共計 6 個節點、21 個自由度。 1969 年文獻[29]中也假設側向位移是五次多項式，並且利用邊上的側向位移對邊上法線方向微分是三次變化三個限制條件將三個邊上的自由度、、去掉，因此 6 個節點減為 3 個節點、21 個自由度減為 18 個自由度，文獻[9]將此元素稱為 RQT 元素。 w 6 ,_n w w i w, w xi w, 4 ,_n yi w, 5 ,_n xxi w, nk w, xyi w, ,_yyi n w, w w 1981 年文獻[1]中將 CST 平面元素與 DKT 三角板元素疊加成一平面三角殼元素，並使用更新拉格朗日法將該元素用在具大位移及大旋轉的薄殼結構分析，但更新拉格朗日法的增量旋轉必須是小角度，因此 1987 年文獻 [4]使用相同的殼元素，搭配共旋轉法解決增量旋轉大小限制的問題。2006 年文獻[6]採用文獻[4]的殼元素和共旋轉法，以數值例題探討殼結構受到各種位移負荷之非線性分析。一個好的殼元素應能分析大位移問題和偵測平衡路徑上的分歧點 (Bifurcation point)及極限點。為了要檢測元素的優劣性，文獻上通常藉由基準問題(Benchmark problem)來試驗推導出來的元素模型是否精準或收斂方

(18)

法是否有效率。文獻[30]整理出殼在幾何非線性分析中常見的基準問題，文獻[7-9,16,31-34]模擬殼在挫屈分析中常見的基準問題，但更好的基準問題是要兼具實驗數據和數值模擬的比較[19]。文獻[19]以實驗和數值方法[35]探討一聚酯圓柱薄殼受位移負荷作用後的非線性行為，模擬一矩形薄板在長邊以夾鉗挾持，夾鉗先將薄板彎成一圓柱狀，再施加一集中位移負荷於結構中心的情況，採用兩階段的位移負荷分析。在其實驗中隨著位移負荷的增加，結構連續產生四個特殊的幾何變形，如圖1.1所示。第一個變形是在薄殼中心附近出現兩個對稱X、Y軸的 d-cone (developable cone) (圖1.1a)。第二個變形中出現兩個新的d-cone，而

四個d-cone圍成一個對薄殼中心轉了一個角度的菱形(圖1.1b)。第三個變形為四個d-cone的連線形成一個梯型(圖1.1c)。第四個變形為梯形底邊兩個 d-cone移到薄殼自由端的邊界時，產生一個不連續的變化，使薄殼變成波浪狀的圓柱殼(圖1.1d)。文獻[35]所使用的數值方法為有限元素法，其元素為三角形殼元素，考慮的薄膜應變為不完整非線性項，並且假設位移場是完整五次多項式，多項式的21個係數由節點上的函數值、函數的一次微分、函數二次微分以及函數對邊上中點之法線方向微分計算而得，計算von Karman板殼理論推導出來的能量式，利用找出該能量的最小值求得薄殼結構的變形。雖然在薄殼變成波浪狀前，數值模擬和實驗得出的受力-位移曲線相當接近，但數值模擬無法得出第二個和第三個變形。文獻[7]採用文獻

(19)

[6]的數值方法及元素模擬文獻[19]的實驗，並討論位移負載偏移、結構的不完美對結構變形的影響，雖然與文獻[19]的挫屈例題之數值做比較，其數值的結果很相近，但無法觀察到實驗中出現的現象。文獻[8]採用QST平面元素和DKT板元素疊加成一三角平面殼元素，並探討採用不同位移場推導得的幾何剛度矩陣對平衡迭代及挫屈負荷的影響，與文獻[19]所偵測的平衡路徑和挫屈負荷相近，但無法觀察到實驗中d-cone旋轉的現象。文獻[9]採用 QST元素和RQT板元素疊加成一三角平面殼元素，此殼元素與文獻[8]相比的優點:擁有更高的精確度、元素內部滿足 x y y x     



的連續條件及具有定義的側向位移場可以推導幾何剛度矩陣。數值模擬的結果與大部分具挫屈分析例題相近，但無法觀察到實驗中d-cone旋轉的現象。本人猜測文獻[7-9]無法模擬出 d-cone 旋轉現象的原因在於其推導出的切線剛度矩陣並不完整和沒有考慮到無限小旋轉矩陣(Infinitesimal rotation matrix)的影響，故本篇研究將引用非線性的殼理論來描述殼的變形過程，將變形之間的耦合項考慮更加完善，試圖模擬出實驗中[19]d-cone 旋轉的現象。廣為人知的殼理論推導可分兩種：第一種是將殼結構視為三維物體 [36-39]，並設參考表面，通常取中平面，再對厚度方向積分。第二種是將殼結構視為 Cosserat surface[40]，在此表面上的任一向量皆可變形，且保有在剛體運動下，向量長度不變的性質。

(20)

本研究將三角殼元素的三個節點建立在板的中平面上，在三個節點當前的位置建立一個元素座標，並在當前的元素座標上用板中平面的位移及中平面變形後的法向量描述殼元素的變形。本研究利用極分解定理(Polar decomposition theorem)[41]將殼中平面的變形梯度(Deformation gradient)分

解成一個旋轉矩陣(Rotation matrix)，和一伸縮矩陣(Stretch matrix)的乘積，並用一剛接在元素中平面的座標系統的旋轉表示位移梯度中的剛體旋轉。本研究用三個旋轉參數來描述該中平面座標系統的旋轉。本研究採用 von Karman 板理論為基本假設[42,43]，即第一種殼理論推導，取中平面為參考平面對厚度方向積分的三維物體概念。位移場採用文獻[23,28]的元素位移場，推導出一個 3 節點、每個節點有 9 個自由度的平面三角殼元素，元素的節點自由度為節點位移向量的 3 個分量、節點旋轉向量的 3 個分量及節點平面應變的 3 個分量。廣義來說，應變也為物理上的自由度，且在文獻[9]的數值比較下，大部分的例題也有很好的結果。據本人所知，文獻上尚未有其他人用此方式來描述殼的非線性變形，故本研究將採用此方式推導一平面三角形殼元素，並使用共旋轉法將該元素用在具大位移及大旋轉的薄殼結構分析中，希望能夠有效地改善殼結構幾何非線性分析中平衡迭代的收斂速度，並精確偵測平衡路徑上的分歧點及挫屈模態。本文在第二章介紹本研究所使用的平面三角殼元素的變形機制以及推

(21)

導元素節點內力及剛度矩陣。在第三章說明本研究的數值計算方法和分析時的數值程序。在第四章中以非線性例題測試本研究所使用的平面三角殼元素的性能，以及說明本文提出決定元素節點變形參數的方法是可行的。

(22)

第二章

理論推導

本章將採用共旋轉全拉格蘭日推導法(Co-rotational total Lagrange

formulation)、von Karman 板理論[42,43]、虛功原理、一致性二階線性化 (Consistent second order linearization)推導出一個 3 節點、每個節點有 9 個自

由度的平面三角殼元素，元素的節點自由度為節點位移向量的 3 個分量、節點旋轉向量的 3 個分量及節點平面應變的 3 個分量。在本章將說明座標系統、殼元素變形之基本假設、變形描述、應變、元素節點內力及剛度矩陣的推導。 2.1 基本假設 本文中對非線性平面三角殼元素的推導，做以下假設： (1) 板為均勻厚度的薄板。 (2) 在元素變形前，垂直於元素中平面的線段，在元素變形後，依然保持直線、垂直於變形後之中平面且沒有伸長或縮短。 (3) 元素的變形位移及旋轉為小位移及小旋轉。 (4) 元素的應變為小應變且僅考慮面內的應變。假設(2)為 von Karman 板的變形假設。因採用共旋轉推導法，只要元素夠小，假設(3)一定可以滿足。

(23)

2.2 座標系統

為了描述系統的運動及元素的變形，本文定義了四組座標系統： (a) 固定總體座標系統(Global coordinate system)：X_iG



i1,2,3



結構體所有節點的座標、位移、旋轉、系統的邊界條件與其他座標系統的基底，以及結構的平衡方程式，均在此座標系統中定義。

(b) 元素座標系統(Element coordinate system)：x_iE



i1,2,3



此座標系統是建立在每一個元素變形後的最新位置上，先將座標原點定在元素節點1，x₁E軸為元素節點1與元素節點 2在元素平面上的連線， E x₂ 軸是在元素平面上垂直於x₁E軸，且朝著元素節點3的方向，x₃E軸則是由x₁E軸與x₂E軸外積而得，再將座標原點移至元素三個節點決定之三角形的形心並旋轉



角使得當前變形位置的元素與初始未變形時的元素在形心處無旋轉作用，詳細說明在附錄A。元素的位移、元素變形、元素節點內力與元素剛度矩陣是在此座標系統中定義，然後經由座標轉換，將其轉換至總體座標系統及基礎座標系統。 (c) 元素中平面座標系統：x_iS



i1,2,3



此座標系統的原點是剛接在元素的中平面上，與中平面一起平移及旋轉，在變形前，其座標軸的方向與元素座標的座標軸的方向一致，在變形後，其x₃S軸的方向與變形後之中平面法線方向一致。本文中將該座標系統對元素座標系統的旋轉視為中平面運動之變形梯度(Deformation

(24)

gradient)中剛體旋轉的部份。

(d) 節點基礎座標系統(Base coordinate system)：x_ijB



i, j 1,2,3



此座標系統的原點是剛接在結構離散後的每一個節點，並與對應的節點一起移動及旋轉， j 為節點在元素中的編號。本文中節點 j 在其初始位 置之x₃B_j軸的方向為曲面在該節點的法線方向，x₁B_j軸、x₂B_j軸的方向為曲面在該節點互相垂直的切線方向，節點的應變自由度是在此座標系統中定義。本文中0x_ijB表示元素節點 j 在初始未變形時的節點基礎座標， B ij I_{x 表示元素節點 j 在第 I 個增量迭代收斂後的節點基礎座標，} B ij x 表示 元素節點 j 在當前變形位置的節點基礎座標。 2.3 旋轉向量 本文中使用旋轉向量來表示一個有限旋轉，如圖 2.1 所示，一向量b 受 到一旋轉向量



a的作用而轉到一個新位置b，向量與b 之間的關係可表示成[45]： b  ) ( sin ) )( cos -(1 cos b a b a a b b







 



 (2.1) 其中符號．與分別代表向量的內積與外積，



表繞單位向量為之轉軸的旋轉角。 a 當旋轉向量{



₁



₂



₃}，



_i (i1,2,3)為 在某一直角座標的分量，有一微小變量



時，會繞該座標的座標軸做微小的旋轉b



_i (i 1,2,3)。當很小時 



{



₁



₂



₃}與



{





₁



₂



₃}之關係可表示成[46]：

(25)

 



T_ (2.2)                     1 2 2 2 1 2 2 2 1 1 2 1 3 2 3



 T (2.3) 2.4 殼的變形描述本文在當前的元素座標上描述殼元素的變形，由2.1節的基本假設(2) 可知von Karman板的變形可由其中平面的位移來描述。如圖2.2中，Q 點為殼中的任意點， P 點為Q 點在同一斷面之中平面上 的對應點。在元素座標上Q 點的變形前後位置可以表示如下： } { 0  x y z r (2.4) n r r{r₁r₂ r₃} _p z (2.5)



x₁ x₂ x₃

 

x u(x,y) y v(x,y) w(x,y) p    



r (2.6) 其中x、、為變形前y z Q點在元素座標x_iE(i1,2,3)上的座標，x、同時也是 y P 點變形前的x₁E、x₂E軸的座標。r 、r_p分別是變形後Q點、P 點在元 素座標x_iE上的位置向量，x₁、x₂以及x₃為P 點在元素座標x_iE的座標值， ) y , (x u 、v(x,y)以及w(x,y)為P 點在 、及軸方向的位移，n 是變形 後中平面在 E 1 x x₂E x₃E P 點的向外單位法線向量，可表示為：              



_n



_cos 1 1 ₁2 ₂2 1 2 2 2 1 2 n (2.7)

(26)

2 1 _x w   



(2.8) 1 2 _x w    



(2.9) 2 2 2 1 1 1 cos



   n (2.10)

本文中將x、y、z 視為拉格蘭日座標(Lagrange coordinates)。由(2.6)式可得：

                                                                             2 1 , , , , 2 1 2 1 2 1 1 1 x w x w v u v u x w x w y x y x x x x x y wx w y y x x _(2.11) 由(2.4)式，殼中任意點的變形梯度(deformation gradient)F [41]可表示如下：                                     z r y r x r z r y r x r z r y r x r z y x 3 3 3 2 2 2 1 1 1 ) , , ( F F (2.12)

由(2.11)式和極分解定理(Polar decomposition theorem) [41]可得：

RU F (2.13)               z zy zx yz y yx xz xy x



1 1 1 U (2.14)

其中R 為一為旋轉矩陣(Rotation matrix)，U為一伸縮矩陣(Stretch matrix)， i



、



_ij (i x,y,z, j x,y,z)為在任意點(x, y, z)的工程應變。由基本假設(3)

(27)

體旋轉，該旋轉可使在矩陣CFtF之主軸方向的向量轉到矩陣之主軸方向[41]。 1 1 1 ₍  ₎   _ _F _F B t i e eS_i S i x



i1,2,3



i e 3 3 3 e θ 令及分別為元素座標軸與剛接在P 點之中平面座標軸 的單位向量。由座標系統的定義可知，在變形前軸與軸的方向是一致的，即與的方向是一致的，而且變形後與(2.7)式的方向一樣。在本文中假設變形後的單位向量是由以下兩個旋轉向量連續作用於單位向量來決定： E i x i x x_iS i e , 1  S i e ) 3 S 3 e n S i e , 2 (i



 n n n e θ (2.15)



 (2.16) } 0 , , { 2 2 2 1 2 2 2 2 1 1



   n e _(2.17) 其中



₃為中平面繞e₃旋轉的角度，



_n為和n 的夾角，為垂直於與n 之 單位向量。 3 e e_n e₃ 將旋轉向量θ₃作用在上，將其轉至一中繼位置e_i e_i，再將θ_n作用在e_i 上，將其轉到。若、θ 以及 已知，則中平面座標就唯一決定；反之，若與已知，則旋轉向量與亦唯一決定。 S i e e_i ₃ θ_n 3 θ S i e i e eS_i θ_n S i e 與之關係可表示如下：e_i i SE i S i R R n e R e e [ ₁ ₂ ]  _(2.18) 2 3 1 3 1 cos r sin r R 







, R₂ sin



₃r₁cos



₃r₂ _(2.19) } sin ) cos 1 ( ) cos 1 ( cos { ₁2 ₁ ₂ ₂ 1 



n n 



n n n 



n n



n r

(28)

} sin ) cos 1 ( cos ) cos 1 ( { ₁ ₂ ₂2 ₁ 2  n n 



n



n n 



n n



n r 2 2 2 1 2 2 2 1 1 sin



    n 其中R_SE為一旋轉矩陣。因R_SE為



_i(i1,2,3)的函數，所以本文中稱



_i為旋轉參數。將R_SE保留到



_i (i1,2,3)的二次項之近似式可表示如下：                               2 2 2 1 3 2 1 3 1 2 1 2 3 2 1 2 1 3 2 2 1 3 2 3 2 2 2 1 2 1 1 2 1 2 1 1 2 1 2 1 2 1 2 1 1



SE R _(2.20) 當旋轉參數



_i(i1,2,3)分別有一微小變化



_i時，中平面座標的單位向量會旋轉到一個新的位置，此一新的位置可由繞元素座標軸軸分別作微小旋轉 S i e (i S i e ) 3 , 2 , 1 x_iE



_i(i1,2,3)而得。





₁



₂



₃





θ 與



{



₁



₂



₃}之關係可表示成：  



θ _ _(2.21)                1 1 1 1 2 1 2 2 2 1 2 2 1 2 1



 a a  _(2.22) 2 2 2 1 2 2 2 1 1 1



      a 將 保留到_



_i (i 1,2,3)的一次項之近似式可表示如下：              1 1 1 1 0 0 1 1 2



  (2.23)

(29)

由(2.21)式可知當



_i (i 1,2,3)不為零時，



_i並非繞軸的微小旋轉。因(2.18)式之隨 P 點一起剛體運動，所以(2.20)式之旋轉矩陣與(2.13)式中P 點的旋轉矩陣 ) 3 , 2 , 1 (i x_iE S i e (i1,2,3) SE R R 應為同一旋轉矩陣。在本文以後的內 容中，除另有說明外，(2.13)式之R 都是指在中平面上的旋轉矩陣R(x,y,0)。將(2.5)式代入(2.12)式，可得變形梯度，再將F z 0代入，可得在中平面之變形梯度 F ) 0 , , (x y F 。將F(x,y,0)之變形參數保留到二次項，可得：               n y x y x y x w w v v u u y x



cos 1 1 ) 0 , , ( , , 1 , , 2 , , F _(2.24) 將(2.24)式、(2.20)式和(2.14)式代入(2.13)式，並將變形參數保留到二次項可得： 3 2 3 2 2 , ₂ 1 2 1

_

_

_



_x _xy x u     (2.25) 2 1 3 , ₂ 1 ) 1 (







    _xy _y y u _(2.26) yx x x v_,  





₃ 



₁



₂ 



2 1 ) 1 ( _(2.27) 2 3 2 1 3 , ₂ 1 2 1

_







    _xy _y y v _(2.28) 1 3 2 ,x (1



x)



(





yx)



w     (2.29) 1 2 3 ,y (





xy)



(1



y)



w     (2.30) (2.25)式－(2.30)式為中平面之位移梯度(Displacement gradient)與旋轉參數 i



(i1,2,3)、中平面之工程應變



_x、



_y、



_xy 2



_xy 2



_yx間二次的非線性耦合關係。本文中將位移梯度u_,_ 、v_,_ 、 w_,_ (



 x, y) 及旋轉參數



_i

(30)

) 3 , 2 , 1 (i 、中平面之工程應變



_x、



_y、



_xy視為兩組獨立的變形參數。本研究用



_i (i1,2,3) ，



_x、



_y、



_xy在元素節點的值及(2.25)式－(2.30)式決定、 , y)  , u v 、_,_ w (_,_



x 在元素節點的值，用再 u 、v_,_ _,_、w (_,_



x, y)在元素節點的值及元素節點的節點位移決定元素的位移場。將(2.25) 式至 (2.28) 式的變形參數保留到一次項可得 u_,_x 



_x 、、和 3



v_, ,y 



xy  u _x 



₃ 



_yx v_,_y 



_y，該結果與線性分析的定義一致。    x y , 將(2.25)式－(2.28)式代入(2.29)式和(2.30)式，並將變形參數保留到二次項可得：        )  w            2 1 (1



 y u u_, 1  x y v v , , ) x , ( w , _(2.31) (2.31)式與(2.11)式之二次近似式的結果是一致的。將(2.31)式代入(2.5)式，並保留變形參數u_,_、v_,_、w_,_ (



x, y)到二次項，可得變形後Q 點在元素座標x_iE上的位置向量之二次近似式為：                (    2 1 x y v



   1 1 ( 1 z       ) , ) , ) , y y x y x         ( ( ( x w v y u x r x y w, 2 , y y y x x



2



, , , 2 1 ) ) x y x u v u



(2.32) (2.33)



y x w,



為了要推導出假設的節點參數之關係，故將(2.25)式至(2.30)式兩邊做變分可

(31)

xy xy x x u



















_,   ₂ ₂ (  ₃) ₃  ₃ (2.34) 1 2 2 1 3 3 , ₂ 1 2 1 ) 1 (





















u _y  _xy   _y  _y   (2.35) yx x x x v





















_,   ₃  ₃  ₁ ₂  ₂ ₁  2 1 2 1 ) 1 ( (2.36) xy y xy y v



















_, (  ₃) ₃   ₁ ₁  ₃ (2.37) yx yx x x x w

























_, (1 ) ₂  ₂ ( ₃  ) ₁  ₁ ₃  ₁ (2.38) y y xy xy y w

























_, ( ₃  ) ₂  ₂ ₃  ₂ (1 ) ₁  ₁ (2.39) 2.5元素節點參數及節點力向量本文推導的平面三角殼元素有3個位於中平面的節點，每個節點有9個自由度，為了簡潔，本文用

 

_j、

 

_j表示

 

、

 

中的每個元素都帶有下標。為了在推導上的需要，本文採用四類的元素節點參數向量，並分別表示如下： j } { _x₁ _x₂ _x₃ x q q q q  (2.40) } { _₁ _₂ _₃  q q q q  (2.41) } { _₁ _₂ _₃ 



q  q q q (2.42) } { _₁ _₂ _₃  q q q q  (2.43) 其中 j xj {u u v w} q     (2.44) j j {u θ ε} q_  (2.45)

(32)

j j { u ε} q



_   (2.46) j j {u ε} q_   (2.47) j (j1,2,3)表示元素局部節點 j ，u_j {u v w}_j j } ，、以及為元素節點位移向量在、及軸方向的位移分量，，，， j u v_j j   u j w x u j u j E x₁ j  w E x₂ x w { _, E 3 y , x w j y u } { _, _, y} x v v { _, _,  j  v  u_, _j ( u)_j



  _ 、 v,j (_

_

v)j 以及 j j w w_, ( )



  _ (



x,y)，θj {θ1 θ2 θ3}j，



ij為作用在節點 j的旋轉參數 i



(i 1,2,3)((2.8)、(2.9)、(2.15)式)，ε_j {



_x



_y



_xy}_j，



_xj、



_yj、



_xyj為中平面之工程應變((2.14)式)的節點值，



_j {



₁



₂



₃}_j，



_ij為作用在節點 j 的無限小旋轉參數，_j {



₁



₂



₃}_j，



_ij為作用在節點 j 的旋 轉向量((2.2)式)_j在x₁E、x₂E及x₃E軸方向的分量。本文中將作用在節點 j的旋轉向量在當前的變形位置重新設定為零，本文在平衡迭代的過程中，用節點旋轉向量的增量或擾動量作用在對應的節點基礎座標系統，決定該座標系統變形後的方位。 j  j  本研究用在元素節點 j (j1,2,3)的節點基礎座標系統當前的方位及當前的元素座標決定，然後用決定，再用決定元素中平面之位移場  q q_ q_x q_x ) , ( yx u 、v( yx, )以及w(x, y)。由(2.34)式－(2.39)式可得擾動節點參數向量



q_xj與



q__j有如下的關係： j j x xj 







q T q (2.48)

(33)

j xy y x yx xy x y xy j x                                                         2 1 2 3 1 2 1 3 3 3 1 3 1 2 3 1 2 3 3 2 2 1 0 ) ( ) 1 ( 0 0 0 2 1 0 ) 1 ( ) ( 0 0 0 2 1 1 0 ) ( 0 0 0 0 2 1 0 ) 1 ( 2 1 2 1 0 0 0 2 1 0 ) 1 ( 2 1 2 1 0 0 0 2 1 0 1 ) ( 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1































 T (2.49) 由(2.21)式可得擾動節點參數向量



q__j與



q__j有如下的關係： j j j   



q T q (2.50) j j                                  1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 2 1 2 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 2 1 2



 T (2.51) 由(2.2)式可得擾動節點參數向量



q__j與



q__j有如下的關係： j j j   



q T q (2.52)

(34)

j j                                     1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 2 1 2 1 0 0 0 0 0 0 2 1 1 2 1 0 0 0 0 0 0 2 1 2 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 2 1 3 2 3



 T (2.53) 由(2.48)式－(2.53)式可得



q_x、



q_、



q_、



q_間的轉換關係如下  



q_x T_x q (2.54)   



q T q (2.55)   



q T q (2.56)            3 2 1     x x x x T 0 0 0 T 0 0 0 T T (2.57)            3 2 1     T 0 0 0 T 0 0 0 T T (2.58)            3 2 1     T 0 0 0 T 0 0 0 T T (2.59) 對應於擾動節點參數向量



q_xj、



q__j、



q__j、



q__j的廣義節點力向量分別為： j w v u xj {f m  m  m } f (2.60) } {f m m f  (2.61)

(35)

j j {  }  f m m f  (2.62) j j {  }  f m m f  (2.63) 其中 j j {f1 f2 f3} f (2.64) j uy ux j u {m m } m ，m_v__j {m_vx m_vy}_j，m_w__j {m_wx m_wy}_j j j {m1 m2 m3}   m j j {m1 m2 m3}   m j j {m1 m2 m3}   m j y x j {m m m}   m 對應於



u_j、



v_j以及



w_j的廣義節點力為在x_iE (i1,2,3)軸方向的力；對應於 ij f j u_



_, 、



v_,__j、



w_,__j(



x,y, 3j1,2, )的廣義節點力分別為廣義力矩、、；對應於 j u m _ m_v__j m_w__j



_ij(i1,2,3, 3j1,2, )的廣義節點力為廣義力矩m__ij；對應於



_ij(i1,2,3, j 1,2,3)的廣義節點力為傳統力矩；對應於 ij m_ xj



、



_yj、



_xyj的廣義節點力分別為廣義力矩m__xj、m__yj、m__j；對應於 ij



(i1,2,3, j1,2,3)的廣義節點力為廣義力矩m__ij。因為



_ij在變形後不為零，所以其變分



_ij並不是繞軸的無限小旋轉(見(2.21)式)，所以廣義力矩並非繞軸的傳統力矩。如前所述，本文中將作用在節點 E i x ij  m x_iE j的旋轉向量((2.2)式)_j在當前的變形位置重新設定為零，所以其變分



_ij與繞軸的無限小旋轉 E i x ij



的值相同，所以廣義力矩m__ij與繞x_iE軸的傳統力矩的值相

(36)

同。不同元素在共同節點的廣義節點力需先轉換到共同座標系統才能組合成在該節點的系統節點內力。 j  f 對應於(2.40)式－(2.43)式之元素節點參數向量的元素節點內力向量可表示成： } { _x₁ _x₂ _x₃ x f f f f  (2.65) } {_₁ _₂ _₃  f f f f  (2.66) } { _₁ _₂ _₃  f f f f  (2.67) } {_₁ _₂ _₃  f f f f  (2.68) 其中f_xj、f__j、f__j、f__j定義於(2.60)－(2.63)式由(2.54)式－(2.56)式與反梯度法則(contragradient law)[47]，可得、、、間的轉換關係如下： x f f_ f_  f x t x f T f_  _ (2.69)    T f f  t (2.70)    T f f  t (2.71) 2.6元素的應變及其變分 2.6.1 應變本文中應變的度量是採用 Green strain，本文中以



_ij(i, j 1,2,3)表示在當前元素座標x_iE的Green strain。由基本假設(2)，本文只考慮



₁₁、



₁₂與



₂₂，並

(37)

表示如下[48]： ) 1 ( 2 1 1 t 1 11 g g 



(2.72) ) ( 2 1 2 t 1 12  g g



) 1 ( 2 1 2 t 2 22  g g 



x    r g₁ , y    r g₂ , z    r g₃ (2.73) 其中r 定義於(2.32)式。將(2.32)式代入(2.73)式，並保留變形參數及其微分到二次項，可得的分量如下： i g xE_j ) 3 , 2 , 1 , (i j g_ij ] ) 1 ( [ 1 _, _, _, _, _, _, _, 11 u x z u xx y u x y x v xx x v x x x g    



 











(2.74) ] ) 1 ( [ _, _, _, _, _, _, , 12 vx z vxy x v y x x u xy y u y y x g  



 











) ( _, _, , 13 wx z x x x y y x g   







] ) 1 ( [ _, _, _, _, _, _, , 21 u y z u xy y u x y y vxy x vx x y g   



 











(2.75) ] ) 1 ( [ 1 _, _, _, _, _, _, _, 22 v y z v yy x v y x y u yy y u y y y g   



 











) ( _, _, , 23 wy z x x y y y y g   







x x y x v u g₃₁ (1 _, )



 _,



(2.76) y y x y u v g₃₂ (1 _, )



 _,



2 2 33 ₂ 1 2 1 1 _x _y g  







(38)

將(2.74)式至(2.76)式代入(2.72)式，並保留參數及其微分到二次項，可得應變為： 2 11 1 11 11



  (2.77) x y x z u_, _, 1 11





  2 , 2 2 , 2 , , 2 , 2 , 2 , 2 11 ₂ 1 2 1 2 1 2 1 2 1 x y x x xx y x xx x x x v w zv z u z z u





       2 12 1 12 12      (2.78) y y x x x y v z z u_, _, _, _, 1 12 ₂ 1 2 1 2 1 2 1 _ _      y x x x y y x y xy x xy y y x y x y xu v v w w z u z v z z u_, _, _, _, _, _, _, _, 2 _, _, 2 _, _, 2 12 ₂ 1 2 1 2 1 2 1 2 1 _ _ _ _ _ _         2 22 1 22 22



  (2.79) y x y z v_, _, 1 22





  2 , 2 2 , 2 , , 2 , 2 , 2 , 2 22 ₂ 1 2 1 2 1 2 1 2 1 y y y x yy y yy x y y y v w z v z u z z u





       其中



_ijk(i, j,k 1,2)表示



_ij中之次項k 為了方便以後推導，本文中令



₁₂ 2



₁₂，並將(2.77)式至(2.79)式的應變表示為向量式： 2 2 2 1 2 1 b ab b a a ε zε zε z ε ε ε     (2.80)            12 22 11





ε (2.81)

(39)

            x y y x a v u v u , , , , 1 ε (2.82)                        y x y x y x y y y x x x a w w v v u u w v u w v u , , , , , , 2 , 2 , 2 , 2 , 2 , 2 , 2 2 1 2 1 2 1 2 1 2 1 2 1 ε (2.83)              y y x x y x x y b , , , , 1



ε (2.84)                xy x xy y yy y yy x xx y x xx ab v u u v u v , , , , , , 2 2 2



ε (2.85)                     y x x x y y x y y y y x x y x x b , , , , 2 , 2 , 2 , 2 , 2 2 1 2 1 2 1 2 1



ε (2.86) 2.6.2 應變的變分 本研究採用虛功原理推導元素節點內力，所以需要應變的變分。將(2.77) 式至(2.79)式變分可得： 2 11 1 11 11



  (2.87) ) ( , , 1 11



ux



y x z



  ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( , 2 , , 2 , , , , , , , , , , , 2 11 x y x y x x x x xx y xx x x xx y xx x x x x x x z z zu zv z v z u w w v v u u



















          

(40)

2 12 1 12 12      (2.88) ) 2 1 ( ) 2 1 ( ) 2 1 ( ) 2 1 ( _, _, _, , 1 12 u y vx x x z y y z       ) 2 1 ( ) 2 1 ( ) 2 1 ( ) 2 1 ( ) ( ) ( ) ( ) ( ) 2 1 ( ) 2 1 ( ) 2 1 ( ) 2 1 ( ) 2 1 ( ) 2 1 ( , 2 , , 2 , , 2 , , 2 , , , , , , , , , , , , , , , , , 2 12 x y y y y y x y x x y x y x x x xy y xy x x xy y xy x y y x x y y x x y y x z z z z zu zv z v z u w w w w v v v v u u u u                                      2 22 1 22 22



  (2.89) ) ( , , 1 22 



v y



x y z



) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( , 2 , , 2 , , , , , , , , , , , 2 22 y y y y y x y x yy y yy x x yy y yy y y y y y y z z zu zv z v z u w w v v u u



















           其中



_ijk(i, j,k 1,2)表示



_ij中之k次項的變分為了方便以後推導，本文將(2.87)式－(2.89)式的應變表示為向量式： 2 2 2 1 2 1 b ab b a a ε z ε z ε z ε ε ε



     (2.90)             x y y x a v u v u , , , , 1



ε (2.91)                     ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( , , , , , , , , , , , , , , , , , , , , , , , , 2 x y y x x y y x x y y x y y y y y y x x x x x x a w w w w v v v v u u u u w w v v u u w w v v u u

平面三角形薄殼元素之共旋轉推導法

國

立

交

通

大

學

機械工程學系碩士班

碩 士 論 文