薄膜皺折之有限元素分析

國 立 交 通 大 學

機械工程學系碩士班

碩

士 論 文

薄膜皺折之有限元素分析

Finite Element Analysis of Thin Membrane

Wrinkling

研 究 生：黃孝衡

指導教授：蕭國模 博士

薄膜皺折之有限元素分析

Finite Element Analysis of Thin Membrane Wrinkling

研究生：黃孝衡 Student：Hsiao Heng Huang 指導教授：蕭國模 博士 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 2007

Hsinchu, Taiwan, Republic of China

薄膜皺折之有限元素分析

研究生：黃孝衡 指導教授：蕭國模 博士 國立交通大學機械工程學系碩士班

摘 要

本研究主要目的為以共旋轉有限元素法及殼元素(shell elememt) 對膜結構皺折後的幾何形狀作初步探討。 本研究採用文獻[19]的平面三角殼元素分析膜結構，並採用牛頓-拉福森(Newton-Raphson)法和弧長控制(arc length control)法的增量迭 代法來解受位移負荷之結構的非線性平衡方程式。 皺折為側方向的變形，由於平面膜結構受到平面上的剪位移作用並 沒有產生側方向變形的機制，本研究利用偵測結構平路徑的分歧點， 計算挫屈模態，利用挫屈模態當擾動位移，使平面膜結構進入次要平 衡路徑，以產生皺折。 本研究分析在邊界受剪位移作用的平面膜結構的皺折，並將分析結 果與文獻上的實驗結果及數值結果比較。

Finite Element Analysis of Thin Membrane Wrinkling

Student：Hsiao Heng Huang Advisor：Dr. Kuo Mo Hsiao

Department of Mechanical Engineering National Chiao Tung University

ABSTRACT

The formation of wrinkling deformation for thin membrane is simulated by using the co-rotational finite element formulation and shell element.

The shell element employed here is the flat three-node triangular shell element proposed by Bathe and Ho’s [19]. An incremental-iterative method based on the Newton-Raphson method and constant arc length method is used for solving nonlinear equilibrium equations with displacement loading.

In order to initiate the out-of-plane buckled deformation for planar membranes subjected to purely in-plane displacement loading at the bifurcation point, a perturbation displacement proportional to the first buckling mode is added to the equilibrium configuration, and then equilibrium iteration is carried out until a new equilibrium state is achieved.

Wrinkling of a planar membrane under uniform shear displacement loading along the edge of the membrane is analyzed here. The present results are compared with experimental data and numerical results in the literature. Good agreement between the present numerical results and experimental data is observed.

誌 謝

衷心感謝指導教授 蕭國模博士兩年來的指導與教誨，使得本論文 得以順利完成。在這兩年碩士班期間，讓我印象最深的是老師在研究 上嚴謹的態度及對我們這群學生的日常生活上的關懷，這都使我受益 良多，在此致上最高的謝意。同時也感謝尹慶中老師以及蔡佳霖老師 撥冗擔任口試委員，為本文提出寶貴的意見。 感謝陳弘虎、蔡明旭、劉峰成、劉宗帆、楊禮龍學長的引領照顧， 以及同學楊水勝的扶持與幫助，還有學弟林育丞、顏宏儒在學業上的 砥礪與成長。 感謝我的父母親、妹妹對我的關愛與照顧，家是我最溫暖的避風 港，感謝奶奶的祝福以及大伯的鼓勵，感謝所有關心我的朋友對我的 支持與協助，感謝學校大門口土地公的保佑。僅以此成果與榮耀，獻 給我親愛的父母以及所有關心我的人。

目 錄

中文摘要 ……….. ….. I 英文摘要 ……….. ….. II 致謝 ……….. ….. III 目錄 ……….. IV 表目錄 ……… VI 圖目錄 ………. VII 第一章 緒論 ……… 1 第二章 理論推導 ……….. 4 2.1 基本假設 ………. 4 2.2 座標系統 ………. 4 2.3 旋轉向量 ………….………..……….. 5 2.4 殼元素變形的描述 .………. 5 2.4.1 常應變三角元素(CST)的變形描述 ………. 6 2.4.2 DKT 元素的變形描述 .……… 8 2.5 元素內力與元素剛度矩陣 ……… 11 2.5.1 CST 元素之節點內力與剛度矩陣 ………. 12 2.5.2 DKT 元素的節點內力及剛度矩陣 ……… 13 2.6 元素幾何剛度矩陣 ………... 14 2.7 元素變形角的描述 ………... 15 2.8 系統的平衡方程式與收斂準則 ……….. 17 第三章 數值計算方法與程序 ……… 19 3.1 本研究使用的數值方法 ……… 19 3.2 增量迭代法 ……….. 19

3.3 弧長控制法 ……….. 22 3.4 二分法……….. 23 3.5 逆冪法 ………..………. 24 3.6 數值程序 ……..………..……….. 24 第四章 數值例題與結果 .……….….. 26 第五章 結論與展望 ……...……… 31 參考文獻 .……….. 32 附錄A DKT 元素的形狀函數 ……….. 86 附錄B CST 元素的剛度矩陣 ……….. 92 附錄C 逆冪法解廣義特徵值問題之數值方法 ………. 93

表 目 錄

圖 目 錄

圖 1.1 皺折的示意圖…...……….………... 36 圖 1.2 波長的示意圖…...……….……... 37 圖 1.3 文獻[11]例題之結構示意圖………...……….. 38 圖 2.1 三角殼元素的示意圖及節點自由度…...………... 39 圖 2.2 總體座標與三種元素座標……… 40 圖 2.3 旋轉向量…..……….………. 41 圖 2.4 DKT 元素的節點及其三邊上的局部座標示意圖…..……… 42 圖 2.5 CST 元素在元素座標上的變形位移………... 43 圖 2.6 變形前板元素中心面之單位法向量n 受旋轉向量θ 作用的 示意圖……….………. 44

圖 2.7 元素座標的剛體旋轉 (a)面外旋轉(out-of plane rotation), (b)面內旋轉(in-plane rotation)………. 45 圖 2.8 決定板元素節點變形轉角的第 3 個步驟的示意圖...…… 46 圖4.1 圓柱殼結構示意圖及其受之位移負荷圖（例題一）……… 47 圖4.2 M 點之反力-負荷參數曲線圖(例題一) ……….. 48 圖4.3 M 點之反力-負荷參數曲線與文獻[18]比較圖(例題一) …… 49 圖4.4 正方形膜結構示意圖及其受之位移負荷圖（例題二）……… 50 圖4.5 挫屈模態(例題二，Mesh 40×40) ……… 51 圖4.6 挫屈模態(例題二，Mesh 6060× ) ……… 52 圖4.7 挫屈模態(例題二，Mesh 8080× ) ……… 53 圖4.8 挫屈模態的剖面圖……… 54 圖4.9 I 點之位移w_I與E點之位移w -_E 負荷參數λ曲線圖 (例題二，Mesh 4040× )……….. 55

圖4.10 I點之位移w_I與E點之位移w -_E 負荷參數λ曲線圖， (例題二，Mesh 6060× )……….……… 56 圖4.11 I點之位移w_I與E點之位移w -_E 負荷參數λ曲線圖， (例題二，Mesh 8080× )……….………. ₅₇ 圖 4.12 不同增量位移下，I點之位移w_I與E點之位移w -_E 負荷參 數λ曲線圖，(例題二，Mesh 6060× )………. 58 圖 4.13 不同元素網格下，I點之位移w_I與E點之位移w -_E 負荷參 數λ曲線比較圖，(例題二)………..………. 59 圖4.14 例題二位移邊界AD上(a)X方向節點力合力

∑

R_Xi–位移負 荷參數λ 曲線圖(b) Y方向節點力合力

∑

R_Yi–位移負荷參 數λ曲線圖……… 60 圖4.15 (a)位移邊界 AD上X 方向的節點反力圖(b)位移邊界 AD上 Y 方向的節點反力圖，(例題二，Mesh 4040× )…………. 61 圖4.16 (a)位移邊界 AD上X 方向的節點反力圖(b)位移邊界 AD上 Y 方向的節點反力圖，(例題二，Mesh 6060× )………62 圖4.17 (a)位移邊界 AD上X 方向的節點反力圖(b)位移邊界 AD上 Y 方向的節點反力圖，(例題二，Mesh 8080× )………63 圖 4.18膜結構變形圖，例題二(a) Mesh 40×40, )λ=0.0507(mm , (b)Mesh60×60, )λ=0.0493(mm , (c)Mesh80×80, ) ( 0507 . 0 mm = λ ………. 64 圖 4.19膜結構變形圖，例題二(a)Mesh 4040× , )λ =0.401(mm , (b)Mesh 6060× , )λ=0.404(mm , (c)Mesh 80×80, ) ( 406 . 0 mm = λ ……….. 65 圖 4.20膜結構變形圖，例題二(a)Mesh 4040× , )λ =0.601(mm ,

(b)Mesh 6060× , )λ=0.604(mm , (c)Mesh 80×80, ) ( 606 . 0 mm = λ ………66 圖 4.21膜結構變形圖，例題二(a)Mesh 4040× , )λ =0.801(mm , (b)Mesh 6060× , )λ=0.804(mm , (c)Mesh 80×80, ) ( 806 . 0 mm = λ ………67 圖 4.22膜結構變形圖，例題二(a)Mesh 4040× , )λ =1(mm , (b)Mesh 60×60,λ=1(mm),(c)Mesh80×80,λ=1(mm)…… 68 圖 4.23軸線 AB之Z方向位移分佈圖，例題二， (a)Mesh 4040× ，(b) Mesh 60×60……….. 69 圖 4.24軸線 AB之Z方向位移分佈圖，例題二， (a)Mesh 6060× ，(b) Mesh 80×80………70 圖 4.25軸線 FH之Z方向位移分佈圖，例題二， (a)Mesh 4040× ，(b) Mesh 60×60……….. 71 圖 4.26軸線 FH之Z方向位移分佈圖，例題二， (a)Mesh 6060× ，(b) Mesh 80×80………72 圖 4.27軸線 CD之Z方向位移分佈圖，例題二， (a)Mesh 4040× ，(b) Mesh 60×60………73 圖 4.28軸線 CD之Z方向位移分佈圖，例題二， (a)Mesh 6060× ，(b) Mesh 80×80……… 74 圖 4.29軸線 JK之 Z方向位移分佈圖，例題二， (a)Mesh 4040× ，(b) Mesh 60×60……… 75 圖 4.30軸線 JK之 Z方向位移分佈圖，例題二， (a)Mesh 6060× ，(b) Mesh 80×80……… 76 圖 4.31軸線 EG之Z方向位移分佈圖，例題二，

(a)Mesh 4040× ，(b) Mesh 60×60……… 77 圖 4.32軸線 EG之Z方向位移分佈圖，例題二， (a)Mesh 6060× ，(b) Mesh 80×80……… 78 圖 4.33軸線 MN之Z方向位移分佈圖，例題二， (a)Mesh 4040× ，(b) Mesh 60×60……… 79 圖 4.34軸線 MN之Z方向位移分佈圖，例題二， (a)Mesh 6060× ，(b) Mesh 80×80……… 80 圖 4.35對角軸線 BD之Z方向位移分佈圖，例題二， (a)Mesh 4040× ，(b) Mesh 60×60……… 81 圖 4.36對角軸線 BD之Z方向位移分佈圖，例題二， (a)Mesh 6060× ，(b) Mesh 80×80……… 82 圖 4.37膜結構皺折與文獻比較圖，(例題二，Mesh 4040× )……… 83 圖 4.38膜結構皺折與文獻比較圖，(例題二，Mesh 6060× )……… 84 圖 4.39膜結構皺折與文獻比較圖，(例題二，Mesh 8080× )……… 85

第一章 緒論

膜結構的歷史由來已久，在自然界中存在許多膜結構形態的東西， 許多動植物細胞以及人體器官都是由膜結構的系統所組成，比如說心 臟，血管，以及皮膚。 膜結構(Membrane Structure) 因為其材料重量輕，可以有大跨度的 延伸，方便折疊以及容易做成各種美觀的形狀而廣泛地應用在各類結 構物上，如航太科技的太陽能板、天線，以及可透光的大跨度輕質屋 頂，展覽場的暫時性結構……等等[1]。 膜結構是依靠膜材自身的張拉力和特殊的幾何形狀而構成的穩定 的承力系統，相對於膜結構承受拉力的能力，其承受壓力和承受彎矩 的能力極弱，故文獻上一般假設膜結構不能承受壓力和彎矩。膜結構 的變形依其受力狀態可以分類成三種型態[2,3,4]： 1. 緊繃(taut)：σ₁ >0, σ₂ >0 膜結構承受雙向(biaxial)的拉應力。 2. 鬆弛(slack)：σ₁ =0, σ₂ =0 膜結構沒有承受任何方向的拉力。 3. 皺折(wrinkle)：σ₁ >0, 0σ₂ = 膜結構既不是緊繃的狀態也不是鬆弛的狀態時，膜結構承受單軸向 的拉應力，另一個主應力恰欲變成壓應力，此時膜結構會產生面外的 (out of plane)的變形，變形如圖1.1所示。 皺折可分為材料性皺折與結構性皺折兩種[5]，材料性的皺折是在外 力卸載後皺折依舊存在，材料性的皺折發生的部位與本來的膜結構會 有不同的結構行為。結構性皺折是由於壓應力的產生所造成的暫時性 挫屈，當外力卸載後皺折也就消失了。 皺折對膜結構的外型，品質，可靠度會有重大的影響，所以研究皺

折發生的方向以及瞭解皺折區域的應力型態是重要的。文獻上討論膜 結構皺折的論文很多[2-17]。文獻上一般以有限元素法分析膜結構， 所採用之元素大都為膜元素(membrane element)[2-7]，因為一般討論 皺折的文章大多從結構的觀點在對膜結構皺折區域的應力與應變進 行探討。 雖然膜元素較簡單，但因膜元素不能承受壓力及彎矩，如果用膜元 素來進行皺折模擬，可以得到皺折的形態(pattern)，卻無法得到皺折 的 細 節 資 訊 ， 比 如 說 振 幅 A (amplitude) 以 及 波 長 λ (wavelength)[5,8-9]，在此振幅的定義為從基準面起算到皺折的最高或 最低點。波長的定義是兩個連續皺折的最高點或兩個連續皺折的最低 點之間的距離[9]，如圖 1.2 所示。而在航太科技的太陽帆船(solar sail) 中，膜結構的反射性(reflectivity)是皺折波長與振幅的函數，因此在太 陽帆船設計上必須要知道皺折的細節[8]。所以亦有少數的論文 [5,8-11,13]利用殼元素(shell elememt)來分析膜結構，在此所定義的殼 元素是考慮抗彎矩能力很弱的薄殼元素。膜結構皺折細節的關鍵在於 膜結構抵抗彎矩的能力(bending rigidity)相對於抗拉應力來說很弱，但 並不是零[3]，故選用殼元素就可以研究這些面外的變形(out of plane deformation)的細節。文獻[11]分析膜結構邊界受剪位移負荷時產生的 變形，圖 1.3 為文獻[11]分析的例題結構示意圖，由於平面膜結構單 純受到平面上的位移負荷，故沒有產生側方向變形的機制，因此文獻 [11]對結構側向加了擾動，使其幾何形狀不完美(Imperfection)，這個 擾動的大小是利用膜結構厚度的函數來計算的。文獻[11]使用10 個4 ABAQUS 軟體的 S4R5 四邊型元素分析膜結構邊界受剪位移後皺折 的振幅，並與文獻[12]的實驗結果比較。 由於文獻上利用殼元素來分析膜結構的論文甚少，且多是利用商用

軟體ABAQUS的殼元素來分析膜結構的皺折[5, 8-11, 13]，甚少提及分 析的細節，故本研究擬採用殼元素探討膜結構受力後完整的行為，以 及膜結構皺折後幾何形狀的細節。文獻[19]中將CST(constant strain triangle)平面元素[20] 與DKT(discrete Kirchhoff theory)三角板元素 [21]疊 加 成 一 平 面 三 角 殼 元 素 ， 並 使 用 更 新 拉 格 蘭 日 法 (updated Lagrangian formulation)將該元素用在具有大位移及大旋轉的薄殼結 構分析，由[19]的數值例題可見該元素相當的簡單及精確。文獻[19] 中在計算對應於DKT元素部分的內力時是將對應於每一増量位移的 増量內力相加而成，但[19]中在計算増量內力時將旋轉視為向量，且 沒有扣除剛體旋轉的部分。因三維有限旋轉並非向量，所以[19]的方 法在兩増量間的增量旋轉必須是小角度，才能得到精確的答案。為了 克服此困難，文獻[22]採用共旋轉法來描述殼元素的變形，並提出一 運動過程來決定平面三角殼元素的節點總變形角，再用總變形角來計 算DKT元素[21]的節點內力。因[22]的方法可以有效的解決兩增量間 增量轉角的限制，故本研究擬採用[22]中的共旋轉法及[19]的平面三 角殼元素。 本研究僅考慮彈性材料的膜結構，及結構性的皺折。本研究在第二 章介紹本文所使用的平面三角殼元素及座標系統。在第三章中介紹本 文的數值計算方法及程序。在第四章中以數值例題驗證本研究所使用 的程式可以偵測結構平衡路徑的分歧點，然後再以文獻上的例題探討 膜結構邊界受到剪位移時，膜結構皺折的情形。

第二章 理論推導

本研究使用文獻[22]的共旋轉法(co-rotational formulation)，及文獻 [19]中提出的平面三角殼元素，如圖 2.1 所示，文獻[19]中的殼元素是 由CST(constant strain triangle)常應變三角形元素[20]及在文獻[21]中 提出的DKT(discrete Kirchhoff theory)三角殼元素所疊加而成。為了本 文的完整性，在本章將簡單描述文獻[19]中元素變形的假設與內力、 剛度矩陣的推導，及本研究所修改的內力及剛度矩陣。文章中亦將簡 單描述文獻[22]中決定元素變形角的方法。 2.1 基本假設 在文獻[19]中對於非線性平面三角殼元素的推導，做以下的假設 (1) 薄膜變形(membrane deformation)以及彎曲變形(bending deformation)之間無耦合作用。 (2) 殼元素的變形為小變形。 (3) 在元素變形前，垂直於中心面的法向線段，在元素變形後，依然 保持直線，且沒有伸長及縮短，除了在元素三個頂點以及三個邊 的中央點外，該線段不必垂直於元素變形後的中心面。 2.2 座標系統 本文採用共旋轉法，為了描述系統的運動以及元素的變形，本文定 義了兩組座標系統： (1) 固定總體座標系統(global coordinate):X_i(i = 1, 2, 3)；如圖 2.2，結 構體節點的座標在此座標系統中定義。 (2) 元素座標系統(element coordinate):x (i = 1, 2, 3)；如圖 2.2，此座_i 標系統是建立在每一殼元素變形後的最新位置上，其座標原點是

在元素節點 1，x₁軸即是元素節點1 與元素節點 2 在元素平面上 的連線，x₂軸是在元素平面上垂直於x₁軸且朝著元素節點3 的方 向，x₃軸則是由x₁軸及x₂軸外積而得。元素變形、元素內力與元 素剛度矩陣是在此座標系統中定義，然後經由標準的座標轉換， 將其轉換至總體座標系統。本文中0x_i及Ix_i (i=1,2,3)表示初始 未變形元素座標及在第 I 個增量疊代收斂後的元素座標。 2.3 旋轉向量 本文中使用旋轉向量來表示一個有限旋轉，如圖2.3 所示，一向量 R受到一旋轉向量

φ

n的作用而轉到一個新位置R′，R 與R′之間的關 係可表示成[23] ) ( sin ) )( cos 1 ( cos R n R n n R R′=

φ

+ −

φ

⋅ +

φ

× (2.1) 其中φ表繞旋轉軸的旋轉角，n 表旋轉軸的單位向量。 2.4 殼元素變形的描述 如圖2.1 所示殼元素中心面上有三個節點，每個節點有6個自由 度，分別是x 、₁ x 、₂ x 軸方向的位移₃ u 、_j v 、_j w 以及繞_j x 、₁ x 、₂ x₃ 軸方向的位移轉角

θ

_xj、

θ

_yj、

θ

_zj。本殼元素假設元素的薄膜變形與彎 曲變形之間無耦合作用，所以元素的變形可由薄膜變形及彎曲變形疊 加而成。本元素的變形是由CST (constant strain triangle)常應變三角形 元素[20]的薄膜變形及文獻[21]中的DKT(discrete Kirchhoff theory)三 角形殼元素的彎曲變形疊加而成。

在圖 2.1 中的元素節點位移u 與_j v_j (j =1,2,3)是 CST 元素節點位 移，而

θ

_xj、

θ

_yj以及w_j(j =1, 2,3)為在[21]中的 DKT元素節點位移。

zj

θ

是為了不使元素剛度內的面內旋轉剛度(in-plane rotational stiffness)為 0，而人為加上去的自由度。圖2.4中的元素是在[21]中所 提出的DKT元素，其節點 1、節點 2、節點 3是元素的三個頂點，節 點 4、節點 5、節點 6為元素的三個邊的中央點，這三個中央點的自 由度僅在元素推導的過程中暫時使用，在最後不會出現在元素的節點 自由度。在本文以下的推導中，元素變形、元素內力以及元素剛度矩 陣都是定義在元素座標系統上。 2.4.1 常應變三角元素(CST)的變形描述 因為CST 元素內的應變為常數，所以其位移場為線性位移場，其 位移場可表示成 y a x a a u = ₁+ ₂ + ₃ (2.2) y a x a a v = ₄ + ₅ + ₆ (2.3) 其中u 跟 v 為在x 軸與₁ x 軸方向的位移， x 與 y 是元素內任意點變形₂ 前的座標值，a_i(i=1−6)是未定常數。 由元素座標的定義可知，在元素座標中圖2.1的節點 1、節點2、 節點3之變形前座標值可表示成(0,0)、(x₂,0)、(x₃,y₃)，令節點 j 在 1 x 軸、x₂軸的位移分別是u_j ,v_j (j=1,2,3)。 將三個節點的座標值及節點位移u_j ,v_j (j=1,2,3)代入(2.2)、(2.3) 可以得到 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ 3 2 1 3 3 2 3 2 1 1 0 1 0 0 1 a a a y x x u u u (2.4)

⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ 6 5 4 3 3 2 3 2 1 1 0 1 0 0 1 a a a y x x v v v (2.5) 由(2.4)、(2.5)式可將a_i(i=1−6)表示成u 及_i v 的函數，所以_i (2.2)、(2.3) 式可改寫成 u=Nu_m (2.6) u={u ,v} (2.7) _⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 3 2 1 3 2 1 0 0 0 0 0 0 N N N N N N N (2.8) u_m =

{

u₁ v₁ u₂ v₂ u₃ v₃

}

(2.9) 1 ( _{2 3 3 2 3} ) 3 2 1 x y xy x y x y y x N = − − + 1 ( _{3 3} ) 3 2 2 xy x y y x N = − 3 3 y y N = (2.10) 依據本文中對元素座標的假設，元素座標原點是定義在節點1，元 素座標的x₁軸是定義在節點1 與節點 2 的連線，如圖 2.5 所示，所以 (2.9)式的u_m中的u₁=v₁ =v₂ =0，故(2.9)式可表示成 } 0 0 0 { u₂ u₃ v₃ m = u (2.11) 其中u₂、u₃、v₃可由圖 2.5 中元素節點在變形前後的座標決定。 CST 元素的應變包含在x₁軸與x₂軸方向的應變

ε

_x與

ε

_y以及剪應變 xy

γ

，因本文中假設元素的變形為小變形所以

ε

_x、

ε

_y、

γ

_xy可表示成

x u x _∂ ∂ =

ε

, y v y _∂ ∂ =

ε

, x v y u xy _∂ ∂ + ∂ ∂ =

γ

(2.12) 將(2.6)式代入(2.12)式得到 m m m B u ε = (2.13) ε_m ={

ε

_x ,

ε

_y ,

γ

_xy} (2.14) ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = x y x y x y y y y x x x m N N N N N N N N N N N N , 3 , 3 , 2 , 2 , 1 , 1 , 3 , 2 , 1 , 3 , 2 , 1 0 0 0 0 0 0 B ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − − + − − + − − = 0 0 0 0 0 0 0 0 1 2 3 3 3 3 2 2 3 3 2 3 3 3 2 x y x y x x x x x x y y y x (2.15) 其中B_m稱為 CST元素的位移-應變轉換矩陣。 2.4.2 DKT元素的變形描述

本文中是採用在文獻[21]中提出的DKT(discrete Kirchhoff theory)三 角板元素，在圖2.6 中，n 為殼元素中心面變形前的單位法線向量，n_d 為n 在元素變形後的新位置，圖2.6 中θ為一在x₁−x₂平面上的旋轉 向量，將θ作用在n 可將n 轉到n_d。由2.1 節所提的假設(3)可知垂直 於變形前的元素中心面法線向量變形後仍為直線且長度不變，由假設 (2)知變形角為小角度，所以文獻[21]之DKT 元素的位移場可表示成 ) , ( yx z u=

θ

_y v=−z

θ

_x( yx, ) w=w( yx, ) (2.16) 其中x、y、 z 為元素上任一點分別在x₁、x₂、x₃軸的座標值，

θ

_y是 θ在x₂軸方向的分量，

θ

_x是θ在x₁軸方向的分量， u是在x₁軸方向的 位移，v是在x₂軸方向的位移，w是在x₃軸方向位移。當 θ <<1時，

y

θ

與

θ

_x可視為法向量n 繞x₂軸及x₁軸的轉角。

DKT元素變形包含彎曲應變(bending strain)

ε

_x、

ε

_y與剪應變

γ

_xy以 及橫向剪應變(transverse shear strain)

γ

_xz、

γ

_yz。

因本文假設元素的變形為小變形，

ε

_x、

ε

_y和

γ

_xy可表示成(2.12)式， xz

γ

、

γ

_yz可表示成 z u x w xz _∂ ∂ + ∂ ∂ =

γ

, z v y w yz _∂ ∂ + ∂ ∂ =

γ

(2.17) 將(2.16)式代入(2.12)式可得 κ ε_b =z (2.18) } , , { _{x y xy} b =

ε

ε

γ

ε (2.19) } , , {

θ

_y,x −

θ

_x,y

θ

_y,y −

θ

_x,x = κ (2.20) 將(2.16)式代入(2.17)式可得 γ={

γ

_xz ,

γ

_yz}={w_,x +

θ

_y , w_,y −

θ

_x} (2.21) 由(2.16)可知w、

θ

_y、

θ

_x與x₃無關，所以可由(2.21)式可知橫向剪 應變在厚度方向為常數。 本文中稱圖2.4 中沿著元素邊緣方向s為切線方向，而垂直於元素 邊緣方向n為法線方向。 在[21]中對於其所提出的DKT元素做了下列的假設 (1)

θ

_y、

θ

_x在元素內為二次變化，也就是

∑

= = 6 1 i yi i y N

θ

θ

,

∑

= − = 6 1 i xi i x N

θ

θ

(2.22) 其中

θ

_yi、

θ

_xi是

θ

_y、

θ

_x在圖 2.4中節點i的節點值，N_i(i =1−6)為形 狀函數[21]，其表示式詳見附錄A。 (2) 元素的頂點以及三個邊的中央點可滿足克希霍夫板理(Kirchhoff plate theory)的假設，即

(a) 在頂點

γ

_xzi =w_,xi +

θ

_yi =0 i=1, 2,3 (2.23a)

γ

_yzi =w_,yi −

θ

_xi =0 i =1,2,3 (2.23b) 其中

γ

_xzi、

γ

_yzi、w 、_xi w 、_yi

θ

_yi、

θ

_xi分別是

γ

_xz、

γ

_yz、 ( ) x w w_x ∂ ∂ = 、 ) ( y w w_y ∂ ∂ = 、

θ

_y、

θ

_x在節點i 的值。 (b) 在三個邊的中央點 0 , = + −

θ

_nk w_sk k =4,5,6 (2.24a) 0 , = + _nk sk w

θ

k =4,5,6 (2.24b) 其中

θ

_nk、

θ

_sk分別是

θ

_n、

θ

_s在節點 k 的值，

θ

_n與

θ

_s分別是θ在 n 與 s 方 向的分量，w_,sk 、w_,nk分別是 _, ( ) s w w_s ∂ ∂ = 、 _, ( ) n w w_n ∂ ∂ = 在節點 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.25) 其中w_i、w 是 w 在節點 i 及 j 的值，_j w_,sk 是w_,s在節點 k 的值， 6 , 5 , 4 = k 分別為邊23、邊31、邊12的中點，ij 邊為節點i 與節點 j 之 間的邊(見圖2.4)，其中i =1−3, j =1−3且i ≠ 。j (4)

θ

_s在元素邊緣是呈現線性變化，即 ) ( 2 1 sj si sk

θ

θ

θ

= + (2.26) 其中

θ

_sk、

θ

_si、

θ

_sj分別是

θ

_s在節點 k 、i、 j 之值，

θ

_s是θ在 s 方向的 分量，在圖2.4 中節點k =4,5,6分別為邊 23、邊 31、邊 12 的中點。 在圖2.4 中元素三個邊上的

θ

_y、

θ

_x與

θ

_s、

θ

_n之幾何轉換關係可表 示成

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ s n y x c s s c

θ

θ

θ

θ

(2.27) x w_, 、w_,y與w_,s、w_,n的幾何轉換關係為 ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ s n y x w w c s s c w w , , , , _(2.28)

其中c =cos

α

_ij、s =sin

α

_ij，

α

_ij為元素的邊ij 上的法線n 與_ij x 軸的夾₁ 角，見圖2.4。 由(2.23)-(2.28)式可以把(2.22)式表示成[21] b T x y H (

ξ

,

η

)u

θ

= b T y x H (

ξ

,

η

)u

θ

=− (2.29)

{

1 x1 y1 2 x2 y2 3 x3 y3

}

b = w

θ

θ

w

θ

θ

w

θ

θ

u (2.30) 其中u_b為 DKT 元素的節點位移，H_x與H_y是對應於元素節點位移的 新形狀函數，其表示式詳見附錄A，

ξ

與η是元素內任一點在元素自 然座標[21]的座標值，其中1≤

ξ

≤0、1≤

η

≤0。 將(2.29)式代入(2.20)式可以得到 κ =B_bu_b (2.31) 其中B_b為 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.32) 其中 2 3 2y x A= ， A為三角形面積。 2.5 元素內力與元素剛度矩陣

本文中殼元素的節點內力是由CST 及 DKT 元素的節點內力組成而 成，元素剛度矩陣是由CST 元素剛度矩陣k 、DKT 元素剛度矩陣_m k_b 以及面內旋轉剛度k_θz所疊加而成。本文中k_θz的值是取k 之對角線元_b 素的絕對值中的最小值。本節中將用虛功原理推導CST 元素及 DKT 元素的節點內力及剛度矩陣。 在平面應力狀態，等向性線彈性材料的應變與應力關係為 σ=Eε , (2.33) } , , {

σ

_x

σ

_y

τ

_xy = σ } , , {

ε

_x

ε

_y

γ

_xy = ε (2.34) ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − = 2 1 0 0 0 1 0 1 1 2

ν

ν

ν

ν

E E (2.35) yz yz G

γ

τ

= xz xz G

γ

τ

= (2.36) 其中 E 是楊氏模數(Young's module)，ν 是蒲松比(Poisson ratio)，G 為 剪力模數，(2.33)式之ε 可以是(2.14)式中的ε 及(2.19)式中的_m ε 。 _b 2.5.1 CST 元素之節點內力與剛度矩陣 將(2.33)式代入(2.13)式可得 σ_m =EB_mu_m (2.37) 由虛功原理可得 =

_∫

V m t m m t mf ε σ dV u

δ

δ

(2.38) } { _{m1 m2 m3} m f f f f =

f_mj ={f_xj f_yj} j =1,2,3 (2.39) 其中f 是 CST 元素對應於_m

δ

u_m的節點內力，V 是元素的體積，u 定_m 義於(2.11)式。 將(2.13)式、(2.37)式代入(2.38)式可得 =

_∫

V m m t m t m m t mf u B EB u dV u

δ

δ

(2.40) 由(2.40)式可得 f_m =k_mu_m (2.41) =

∫

V m t m m B EB dV k (2.42) 其中k_m是CST 元素的剛度矩陣，k_m的表示式詳見附錄B。 2.5.2 DKT 元素的節點內力及剛度矩陣 將(2.18)式、(2.31)式代入(2.33)式可得 σ_b = zEB_bu_b (2.43) 在薄殼中剪應力

τ

_yz與

τ

_xz所做的虛功可以忽略，所以本文中用虛功 原理推導DKT 元素的節點內力時僅考慮

σ

_x、

σ

_y及

τ

_xy所做的虛功。 由虛功原理可得 =

∫

V b t b b t bf ε σ dV u

δ

δ

(2.44) f_b ={f_b1 f_b2 f_b3} f_bj ={f_zj m_xj m_yj} j =1, 2,3 (2.45) 其中f 是 DKT 元素對應於_b

δ

u_b的節點內力，u 定義於(2.30)式，V 為_b DKT 元素的體積。 將(2.18)式、(2.31)式、(2.43)式代入(2.44)式可得 =

∫ ∫

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.46) 其中 ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − = =

_∫

₋ 2 1 0 0 0 1 0 1 ) 1 ( 12 2 3 5 . 0 5 . 0 2

ν

ν

ν

ν

Eh dz z h h b E D (2.47) 其中 h 為板元素厚度。 由(2.46)式可得[21] f_b =k_bu_b (2.48)

∫ ∫

− =2A₀1 ₀1 η t_{b b b} d

ξ

d

η

b B D B k (2.49) 其中k 是 DKT 元素剛度矩陣。 _b 2.6 元素幾何剛度矩陣 為了改善收斂速度，本文中在平衡迭代的過程中將元素幾何剛度矩 陣加入元素剛度矩陣中。本文中採用文獻[19]中的元素幾何剛度，其 表示式為 dA g A t g g B NB k =

∫

(2.50) 其中 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − − = 2 3 32 3 3 2 3 32 3 3 2 3 32 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 x x x y y x x x y y x x x y y A g B

⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = * * * N N N N 其中 x₃₂ =x₃ −x_{2 ⎥} ⎦ ⎤ ⎢ ⎣ ⎡ = 2 12 12 1 * N N N N N N =hEB_mu_m 其中 h 為板元素厚度， E 定義於(2.36)式，B_m定義於(2.15)式，u_m定 義於(2.9)式。 (2.50)式為一近似的幾何剛度矩陣，僅考慮面內應力

σ

_x、

σ

_y以及

τ

_xy 在剛體運動時的效應。文獻[24]提到將元素幾何剛度矩陣加入元素剛 度矩陣可以有效提高平衡迭代的收斂速率。最後，將剛度k 、_m k 、_b g k 疊加後，旋轉至總體座標系統上。 2.7 元素變形角的描述 本文中採用文獻[22]中的兩階段旋轉法來決定元素的剛體旋轉及節 點在元素座標的變形位移與轉角。假設第 I 個位置已知，此處所謂的 第 I 個位置，是指第 I 個增量疊代收斂後的平衡位置。Ix_i(i=1,2,3)為 元素在第 I 個位置元素座標，IX 、_j

∆

U_j及

∆

Φ_j (j =1,2,3)分別是元 素節點 j 在總體座標中第 I 個位置的位置向量、增量位移向量，及增量 旋轉向量，節點 j 當前的座標X_j (j =1, 2,3)可由IX 加上_j

∆

U_j得到。 由X 可以利用元素座標定義求出當前的元素座標_j x_i (i =1,2,3)，本文 中稱元素從Ix 到_i x 的運動為對應於_i

∆

U_j、

∆

Φ_j (j =1,2,3)之剛體運 動，本文中假設該剛體運動是由以下三個步驟達成的。

(1)

∆

U₁造成的位移：元素座標Ix 的原點受到_i

∆

U₁的作用移動到元 素座標x_i的原點，如圖2.7(a)，其中

∆

U₁是作用在元素節點 1的元素 節點增量位移向量。

(2) 旋轉向量α造成的面外旋轉(out-of plane rotation)：在圖 2.7(a)中可 以看到元素的Ix_i (i =1, 2,3)軸因為受到在x₁−x₂平面的旋轉向量α 作用而旋轉，其中Ix₃軸旋轉到x 軸處並且與之重疊， α 的表示式為₃ 3 3 3 3 3 3 1_{( )} cos e e e e e e α × × ⋅ = − I I I _(2.51) 其中Ie_3與 3 e 分別是Ix₃軸與x 軸上的單位向量。₃ 因為旋轉向量α造成的運動是朝遠離x₁ − 平面方向，所以x₂ α所造 成的旋轉稱為面外旋轉。 (3) 旋轉向量β造成的面內旋轉(in-plane rotation)：元素的Ix 軸因為受_i 到在x₁−x₂平面的旋轉向量α作用而旋轉到x′_i (i =1, 2,3)軸，由圖 2.7(b)可知x′ 軸與₃ x 軸重疊，而₃ x′ 軸受到旋轉向量_i β作用轉到x 軸，_i 旋轉向量β的表示式為 β=cos−1(e₁′ ⋅e₁)e₃ (2.52) 其中e₁′ 、e₁、e₃分別是x′ 、₁ x 、₁ x 軸方向的單位向量。₃ 因為旋轉向量β造成的運動的方向是平行於x₁−x₂平面，所以β所 造成的旋轉稱為面內旋轉。 本文使用文獻[22]提出的直接法(direct method)來計算對應於

∆

U_j 及

∆

Φ_j (j=1,2,3)的元素變形角。 在文獻[22]所提出的直接法可以分成以下四個步驟。 (1)

∆

U₁造成的剛體移動： 整個元素因為

∆

U₁而移動，其中

∆

U₁是 作用在元素節點1的元素節點增量位移向量。

(2) 旋轉向量α造成的剛體旋轉：旋轉向量α (2.51 式)通過元素節點1 使得整個元素除了In_dj之外由Ix_i (i=1,2,3)座標轉到x′_i (i =1,2,3) 座標，其中In_dj是在垂直於在第 I 次疊代收斂後位在元素節點 j 的變 形後元素中心面的法線向量。 (3) ∆Φ_tj造成I_{ndj的有限旋轉：法線向量} dj I_{n 受到} tj Φ

∆

作用而旋轉 到新位置n′ ，如圖_dj 2.8所示。其中

∆

Φ_tj是增量旋轉位移向量

∆

Φ_j 在x₁′ −x₂′平面的投影向量。 (4) 旋轉向量β造成的剛體旋轉：旋轉向量β (2.52式)通過元素節點 1 使得整個元素包括n′ 及_u n′ 旋轉到最後的位置_dj n_u及n_dj，其中n_u′ 是 u I_{n 受到α作用後的新位置，} u I_{n 是在第 I 次疊代收斂後，垂直於未變} 形的元素中心面的法線向量，n_u是當前垂直於未變形元素中心面的 法線向量。 令n_u及n_dj為將n_u與n_dj表示成在當前的元素座標x_i (i=1,2,3)軸 分量的向量，則元素變形角θ_j可表示成 dj u dj u dj u j j j n n n n n n θ × × ⋅ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ = cos− ( ) 0 1 2 1

θ

θ

(2.53) 2.8 系統的平衡方程式與收斂準則 結構系統受位移負荷時，其平衡方程式可以表示為

(

Q Q

)

0 F Ψ= ,

λ

_P =

ٛ

(2.54) 其中 Ψ 為系統不平衡力向量，F為系統節點內力， Q 為系統位移向 量，

_λ

為負荷參數，_{Q 為參考位移負荷向量。P F}可由(2.39)與(2.45)

式的元素節點力，用標準的座標轉換，轉換到對應總體座標後再組合 而成。

本文以不平衡力 Ψ 的weighted Euclidean norm 作為迭代時的誤差 度量，而且收斂準則表示為 tol e e e= ≤ F Ψ (2.55) 其中_{F 為對應於位移負荷的系統節點反力向量。e etol}為一給定之容許 誤差值。

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

本文解非線性平衡方程式(2.54)式的數值計算方法是採用文獻[25] 中所提出基於牛頓-拉福森(Newton-Raphson)法和弧長控制(arc length control)法的增量迭代法。為了文章的完整性，本章將簡單介紹文獻[25] 中提出的數值計算的方法與程序。 3.1 本研究使用的數值計算方法 本研究使用牛頓-拉福森法和弧長控制法[25]的增量迭代，本研究利 用判斷系統切線剛度的行列式值是否變號並使用二分法來偵測分歧 點，找到第一個分歧點後，使用逆冪法計算挫屈模態，以第一個挫屈 模態當擾動位移，使平面膜結構進入次要平衡路徑，以產生皺折。本 研究在結構進入次要平衡路徑之後，使用弧長控制法容易出現找不到 實根的情形，故改以標準牛頓法進行迭代。 3.2 增量迭代法 在增量迭代法中，若第 I 個增量的平衡位置為已知，令其位移向量 為Q 、負荷參數為_I

λ

_I，則第I+1 個增量的初始增量位移向量 Q

∆

，可 利用尤拉預測值( Euler predictor )求得[26]

∆

Q =

∆

λ

r_T ， (3.1) I _P T T K F r =−( )−1 ， (3.2)

∂λ

∂

F F_P = ， (3.3) 其中K_TI 為在第 I 個平衡位置的系統切線剛度矩陣，

∆

λ

為初始增量負 荷參數，F 為(2.54)式中的系統節點內力，F 之計算方法可說明如下_P : 因式(3.3)中的系統節點內力 F 是由在總體座標中的元素節點內力

G f 組合而成，因此，F_P亦可由 ∂λ ∂ = G G P f f ， (3.4) 組合而成；(3.4)式中之 ∂λ ∂_fG 可表示成 G P G G G G G r k r r f f ₌ ∂λ ∂ ∂ ∂ = ∂λ ∂ (3.5) G G G r f k ∂ ∂ = (3.6) ∂λ ∂ = G G P r r (3.7) 其中，k 為在總體座標上的元素剛度矩陣，G _{k 是由在元素座標中的}G 元素剛度矩陣經標準轉換得到，r 為元素的節點參考位移負荷，可_PG 以由系統參考位移負荷Q_P萃取而得到。 (3.1)式中的

∆

λ

可利用下式求出[27]

( )

12 T t T l r r

∆

λ

∆

= ± ， (3.8) 其中正負符號之決定方法如下：若第 I-1 與 I 個增量收斂時，其系統 切線剛度矩陣之行列式值同號，則

∆

λ

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

∆

表示第 I+1 個增量的增量弧長，其值可以 如下決定[27]

(

)

12 I D I J J l l

∆

∆

= ， (3.9) 其中J_D為給定的期望迭代次數，J_I為第 I 個增量迭代至平衡所使用 的迭代次數，

∆

l_I 為第I 個增量的增量弧長。

本文中第一個增量的增量弧長

∆

l₁是由下式決定 c l r I R max 0 max 1 R =

∆

， (3.10) 上式中R_max為給定的參考自由度之最大位移量，R₀為

λ

=0時(3.2)式 中之r ，_T R 為₀ R 的₀ Euclidean norm，I_max為給定之最大增量次數，

c r 為R₀在參考自由度的分量之絕對值。 當增量位移向量

∆

Q及增量位移負荷參數

∆

λ

已知時，則可進行以下 的平衡迭代：由Q_{I 1+} =Q_I +

∆

Q，

λ

_{I 1+} =

λ

_I +

∆

λ

，參考位移負荷向量 P Q 及2.5 節與 2.7 節的方法，則可以求得系統中各殼元素在當前的元 素座標、節點變形位移、節點變形角。再利用2.6 節的方法，求得在 當前元素座標上的節點內力及剛度矩陣，再轉換到總體座標的節點內 力及剛度矩陣。由(2.54)式此時之不平衡力可表示成 ΨٛI+1 =F

(

QI+1,λI+1QP

)

， (3.11) 將(2.54)式在

λ

=

λ

_I+1，Q =Q_I+1時用泰勒展開式展開為 = ₊ + + δλ+ λ ∂ ∂ δ ∂ ∂ F Q Q F Ψ Ψ ٛ I ₁ (二次以上高次項) =Ψ_I+1+K_Tr+

δλ

F_P +(二次以上高次項) (3.12) 其中， Q F K

∂

∂

= T 為系統切線剛度矩陣，r =

δ

Q為增量位移修正量，FP 在(3.3)式中已定義，在此我們必須注意F_P在每一次迭代後都要重新 計算。 依牛頓法，忽略(3.12)式中二次以上高次項，可得到位移修正量 r =−K_T−1(Ψ_I+1+

δλ

F_P)=r₀ +

δλ

r_T， (3.13)

其中，r₀ = −K_T−1

ٛ

Ψ_I+1，r 在_T ( 3.2 )式中已定義。 (3.13)式中的位移負荷參數修正量δλ有很多不同的決定方法[28]， 如Crisfield 在文獻[27]所提出的定弧長控制法及 Fried 在文獻[29]提出 的正交軌線法，本文使用Crisfield 在文獻[27]所提出的定弧長控制法 決定δλ，並將於3.3 節中介紹。 將求得之位移負荷參數修正量

δλ

及增量位移修正量r 加入上次迭 代之

∆

λ

與

∆

Q中，可得到一新的增量位移向量與增量位移負荷參數， 再進行下一次迭代。以上之迭代過程一直重複至(3.11)式中的不平衡 力滿足(2.55)式的收斂準則為止。 3.3 弧長控制法 (3.13)式中的位移負荷參數修正量δλ可利用文獻[27]中 Crisfield 所 提出的定弧長控制法決定，其方法使在每一個增量中固定其增量位移 的Euclidean norm 為一定值 l

∆

。由新的增量位移向量為

(

∆

Q+r

)

可得

(

) (

t _T

)

T l

∆

Q r

δλ

r

∆

Q r

δλ

r

∆

2 = + ₀ + + ₀ + _{， (3.14)} 上式經過整理後可以得到δλ的二次方程式 a₁

δλ

2 +a₂

δλ

+a₃ =0， (3.15) 其中 t _T T a₁=r r ，

(

)

t _T a₂ = 2

∆

Q+r₀ r ， a₃ =

(

∆

Q+r₀

) (

t

∆

Q+r₀

)

−

∆

l2， 當(3.15)式的解為兩實根

δλ

₁，

δλ

₂時，本文中取兩根中能使新的增 量位移

(

∆

Q+r

)

與前次迭代增量位移

∆

Q間的內積(inner product)較大

者為δλ。當(3.15)式無實根時，則把(3.9)式之 l

∆

減小，重做第I+1 個 增量。 3.4 二分法 利用3.2 節的增量迭代法可以求得結構之主要平衡路徑。在每個增 量的迭代收斂時，可以得到該增量在其平衡位置的負荷參數λ 及結構 切線剛度矩陣的行列式值D(λ 。令) λ 及_I D(λ 分別表示第 I 個增量在_I) 其平衡位置的λ 及D(λ 之值。) λ_I+1及D(λ_I+1)分別表示第 I+1 個增量在 其平衡位置的λ 及D(λ 之值。) ∆_lI+1表示第I+1 個增量的增量位移向 量之弧長。若D(λ 大於零且_I) D(λ_I+1)小於零則可利用以下二分法求得 挫屈負荷參數λ ： _CR (1) 令_∆lL =0,∆_lR =∆_lI+1,λ_L =λ_I,λ_R =λ_I+1,其中下標 L 及 R 表示左界 及右界。 (2) 取 2 R L 1 I l l l = ∆ + ∆ ∆ ₊ ，重作第I+1 個增量迭代，並求得新的λ 及_I+1 ) D(λ_I+1 。 (3) 若D(λ_I+1)大於零，則令λ_L =λ_I+1,∆_lL = ∆_lI+1，若D(λ_I+1)大於零， 則令λ_R =λ_I+1,∆_lR =∆_lI+1。 (4) 若下列兩式挫屈誤差準則同時滿足 I 1 e_D ) 0 ( D ) ( D < λ ₊ (3.16) λ + < λ λ − λ e 1 I L R _(3.17) 其中e_D，e_λ為給定的容許誤差值，則λ_I+1即為系統挫屈負荷，否則回 到步驟(2)，重新展開下一次二分法迭代。

3.5逆冪法( Inverse Power Method ) 由3.4 節的二分法僅能求得挫屈負荷，無法得知其挫屈模態為了求 得挫屈模態，本文中採用了逆冪法，並說明如下： 令矩陣A 與矩陣 B 分別為 A=K_T(λ_CR), B=I (3.18) 則對應於挫屈負荷參數λ_CR的挫屈模態為對應於特徵值為λ_CR之廣義 特徵值問題的特徵向量 AX=λBX (3.19) 本文用逆冪法解式(3.19)的數值程序放在附錄C中。 3.6 數值程序 本文使用的增量迭代法之數值程序可以分為三個主要部分 1. 輸入與計算開始分析所需的資料 (a) 輸入結構與負荷資料。 (b) 選擇一個參考自由度，並設定期望此自由度達到的最大位移。 (c) 給定最大增量數、每個增量期望的迭代數與最大的迭代數、收 斂時的容許誤差。 (d) 形成剛度矩陣並求得(3.10)式中的R 。₀ (e) 利用(3.1)-(3.3)式、(3.9)式與(3.10)式計算第一次增量的增量弧 長、增量負荷參數與初始增量位移。 2. 使用迭代法求增量的收斂解 (a) 利 用 已 知 的 增 量 位 移 求 得 目 前 元 素 的 變 形 ， 並 以(2.39)、 (2.45)式計算元素節點內力。 (b) 計算(2.54)式的不平衡力Ψ 。 (c) 檢查(2.55)式的收斂準則，若滿足則進行(e)；否則檢查迭代數， 如果小於給定之最大迭代數，則進行(d)；否則減少增量弧長並以

(3.1)-(3.3)式與(3.9)式計算新的增量位移與增量負荷，回到步驟(a)重 新此增量。 (d) 以(3.13)式與3.2節的方法計算增量位移修正量與增量負荷參數 修正量，然後回到步驟(a)。 (e)判斷系統切線剛度的行列式值是否變號，並使用二分法來偵測分 歧點，若找到第一個分歧點，則進行(f)，否則進入第三部分。 (f)使用逆冪法計算挫屈模態，進入第3部分，且此後第二部分改用 標準牛頓法進行迭代。 3. 計算下一次增量所需的資料 (a) 檢查參考自由度的位移或已進行的增量次數是否已達給定值。 若已達到給定值則停止分析工作；否則進行下一步工作。 (b) 計算(3.12)式中的切線剛度矩陣與F_P。 (c) 若第二部分有進行(f)步驟，則調整所算出的挫屈模態向量的長 度，並以此挫屈模態作為下一增量的增量位移，擾動進入次要路徑。 若第二部分沒有進行(f)步驟，則以(3.1)-(3.3)式與(3.8)-(3.9)式計算下 一次增量的增量位移、增量弧長與增量負荷參數。 (d) 回到第二部份執行迭代工作。

第四章 數值分析與結果

在本章中，將先以文獻上的數值例題驗證本文第二章之板元素及第 三章所提的幾何非線性分析之數值程序與方法的正確性及可行性，再 以文獻上的例題探討膜結構受到剪位移負荷作用下的非線性行為。 例題一：圓柱殼受到位移負荷作用 圖4.1 為圓柱殼結構示意及其受之位移負荷圖，如圖4.1所示，AD， BC 兩邊為鉸接，AB，CD 兩邊為自由邊，本例題在圖 4.1 的中點 M 施加一個向下的位移負荷λ_M ，其中λ_M =λ(mm)，λ 是負荷參數 (Loading parameter)，由於結構為雙對稱結構，本例題分別取1/4個結 構，與 1/2個結構，還有完整結構進行分析。分析結構平衡路徑的主 要路徑時，本例題使用10×10元素網格，容許誤差值取1×10−4，取 1/4個結構分析時，使用 20 個增量，平均迭代次數約為4，取 1/2個 結構分析時，使用 20 個增量，平均迭代次數約為 4，取完整結構分 析時，使用 23 個增量，平均迭代次數約為 4。分析結構平衡路徑的 次要路徑時，本例題使用10×10元素網格，容許誤差值取1×10−4，取 1/2個結構分析時，使用 27個增量，平均迭代次數約為4，取完整結 構分析時，使用33個增量，平均迭代次數約為4。圖4.2為本例題所 得到 M 點的位移-負荷曲線圖的結果。由圖 4.2 可以看出本例題在取 1/2 個結構或整個結構作分析時，都可以偵測到平衡路徑上的分歧 點，並可以挫屈模態擾動進入結構的次要平衡路徑。但在取1/4個結 構分析時，無法找到平衡路徑上的分歧點，也無法找到結構的次要平 衡路徑，這應是因為若僅取1/4結構進行分析，在設定邊界條件的時 候，將圖4.1 上 M點的X 與 Y方向位移鎖住了，所以不會偵測出結 構平衡路徑上的分歧點。圖 4.3 為本例題與文獻的比較，其中R_M 為

在位移負荷方向的反力，R_M及 λ相當於文獻之力負荷及力負荷方向 的位移，由圖4.3 可以看出本例題的結果與文獻[18]的結果相當重合。 例題二：正方形膜結構邊界受剪位移負荷 圖4.4 所示為本例題分析之正方形膜結構與受負荷圖，本例題是參 考文獻[11]提出的例題，正方形膜結構邊長L=229 (mm)，厚度 0762 . 0 h= (mm)，楊氏係數E=3790 (MPa)，蒲松比ν=0.38。如圖 4.4，本例題AB，CD邊為自由邊，BC邊為固接(Fixed End)，AD邊 為固接，但可在X軸方向上移動，本例題在 AD邊上施加均勻的水平 位移λ (mm) (圖4.4)。λ稱為負荷參數(loading parameter)，在施加剪位 移λ的過程中，AD 與BC兩邊之間的距離是維持固定不變的。 本例題的挫屈負荷為切線剛度行列式值的正負號恰欲由正變負時 的位移負荷，40×40元素網格挫屈負荷值為6.56×10−4 (mm)，60×60 元素網格挫屈負荷值為6.55×10−4 (mm)，80×80元素網格挫屈負荷 值為6.549×10−4 (mm)，可以看出使用80×80元素網格挫屈負荷的值 收斂了，本例題使用80×80元素網格對膜結構進行分析，並與使用 40 40× 以及60×60元素網格的結果以及文獻[11,12]作比較。使用 40 40× 元素網格例題在結構挫屈後使用標準牛頓法進行迭代，平衡迭 代容許誤差為5×10−5，使用 110 個增量，每個增量平均迭代次數為 約為 16 次，使用60×60元素網格例題結構挫屈後使用標準牛頓法進 行迭代，平衡迭代容許誤差為5×10−5，使用 202 個增量，每個增量 平均迭代次數為6，使用80×80元素網格例題結構挫屈後使用標準牛 頓法進行迭代，平衡迭代容許誤差為1×10−4，使用 203 個增量，每 個增量平均迭代次數為5。圖 4.5-4-7為膜結構的第一個挫屈模態，可 以看出使用不同元素網格計算出來的挫屈模態型態非常相似。圖 4.8

為挫屈模態的剖面圖。圖4.9-4.11為膜結構正中央點I 點之位移w_I與 AB 邊中點 E 點之位移w -_E 負荷參數 λ 曲線圖，I 點與 E 點位置如圖 4.4所示，由圖 4.9-4.11可以看出w_I隨著負荷參數λ增加一直在增加， 且使用40×40元素網格與60×60元素網格以及80×80元素網格w_I的 方向相反。w_E在初期會有不穩定的跳動現象，之後則趨於穩定，圖 4.12是使用60×60元素網格，w_E在使用不同的增量負荷參數∆ (mm)λ 進行迭代的結果，可以看出w_I曲線是趨於穩定的，但w_E會有不同的 結果，這可能是膜結構自由邊邊界鬆弛了，所以使的結構受到差異極 小的負荷條件，但結構的變形差異卻很大，圖 4.13 是使用60×60元 素網格以及80×80元素網格的比較圖，可以觀察到w_I的曲線是接近 的，但w_E的曲線仍有差異。圖 4.14(a)為位移邊界 AD 之節點合力

∑

RXi –位移負荷參數λ 曲線圖，圖4.14(b)為位移邊界 AD之節點合 力

∑

R_Yi–位移負荷參數λ 曲線圖，圖4.14中之反力

∑

R_Xi，

∑

R_Yi 各為圖4.4 中的位移邊界AD之節點 X與Y 方向反力的和，由圖4.14 可看出不同網格下，位移邊界AD 之節點反力的合力

∑

R_Xi，

∑

R_Yi– 位移負荷參數λ 的曲線相當接近。圖 4.15-4.17 為位移邊界 AD 在 X 與Y軸方向上的節點反力。由圖 4.15-4.17 可以看出越靠近A點的地 方節點反力越大，圖4.18-4.22是膜結構變形圖，由圖4.18-4.22可以 看出皺折的走向是呈AC的方向，膜結構最大拉應力的分佈的應是沿 著對角線AC的走向，隨著位移邊界 AD上的點越靠近 D點，拉應力 逐漸變小，故位移邊界上之節點力的分佈是合理的。 圖4.18-4.22為等高線圖及透視圖，透視圖是使用視角53o繪製。膜 結構側方向最大正負位移，以及最大正負位移所在位置的初始座標如 表4.1 所示。由圖4.18-4.22可以看出本例題使用60×60與80×80元素 網格，膜結構的皺折數目較使用40×40元素網格多。圖 4.23-4.36 為

結構重要軸線在 Z 方向上的變形，由圖 4.23-4.36 可以看出隨著負荷 參數λ 越來越大，膜結構在Z方向的位移也越來越大，且結構變形相 當對稱。由圖 4.31與 4.32 軸線 EG 在 Z方向的位移分佈圖可以看出 使用60×60與80×80元素網格，膜結構的皺折數目較使用40×40元素 網格膜結構的皺折數目多，且膜結構正中央點的位移方向相反，故分 析本題目時若使用40×40元素網格可能還不夠密，由圖 4.31 與 4.32 也可以看出使用60×60與80×80元素網格膜結構整體變形趨勢是很 接近的。 圖 4.22 是本例題在負荷參數λ=1(mm)時的膜結構變形圖，可以看 出使用40×40與60×60元素網格及80×80元素網格結果是有差異 的，使用60×60與80×80元素網格結果較接近。使用40×40元素網 格 ， 主 要 有 三 條 突 起 的 皺 折 ， 膜 結 構 側 方 向 最 大 正 負 位 移 為 = ±w_max 1.666， -2.805(mm)， 最 大 正 位 移 與 最 大 負 位 移 的 差 ) ( 471 . 4 ) ( _max max w mm w − − = + ，使用60×60元素網格，主要有四條突 起 的 皺 折 ， 膜 結 構 側 方 向 最 大 正 負 位 移 為 ±w_max = 1.024， -2.376(mm) ， 最 大 正 位 移 與 最 大 負 位 移 的 差 ) ( 400 . 3 ) ( _max max w mm w − − = + ，使用80×80元素網格，主要有四條突 起 的 皺 折 ， 膜 結 構 側 方 向 最 大 正 負 位 移 為 ±w_max = 0.930， -2.153(mm) ， 最 大 正 位 移 與 最 大 負 位 移 的 差 ) ( 083 . 3 ) ( _max max w mm w − − = + 。圖 4.37-4.39 的(b)圖為文獻[12]實驗的 結果，可以看出有四條主要的皺折，膜結構側方向最大正負位移為 = ±w_max 0.67, -3.11(mm)，最大位移與最小位移大約相差3.78(mm)。 圖 4.37-4.39的(c)圖為文獻[11] 以 ABAQUS 以104個 S4R5殼元素數 值模擬的結果，可以看出主要是有三條突起的皺折，膜結構側方向最 大正負位移為±w_max =1.18, -3.15(mm)，膜結構側方向最大位移與最

小位移大約相差4.33(mm)，以本例題與文獻[11,12]比較，可以看出本 例題以40×40元素網格分析的結果與文獻[11]數值例題的結果較接 近，都有三條突起的皺折，膜結構側方向最大正負位移的差也接近， 但是本例題以40×40元素網格分析的膜結構側方向最大正位移比起 文獻[11]的結果大了許多。本例題使用60×60與80×80元素網格分析 的結果與文獻[12]實驗結果較接近。主要有四條主要突起的皺折，膜 結構側方向最大正負位移的差也接近，但是本例題以80×80元素網格 分析的側方向最大正位移比起文獻[12]的結果大，最大負位移比起文 獻[12]的結果小，比較本例題使用80×80元素網格分析的結果與文獻 [12]實驗結果，可以發現兩者的皺折型態有相同的特徵，但細部仍有 差異。造成差異的原因有可能是本例題使用80×80元素網格還不夠 密，也有可能是文獻[12]在做此實驗時，有先將膜結構稍微彎曲一下， 才開始施加剪位移，這部分在本例題的數值計算中沒有被考慮到，也 有可能文獻[12]實驗真正的邊界條件與本例題所設定的邊界條件有差 異，也有可能是膜結構某些地方很容易鬆弛，差異極小的負荷條件也 會造成不同的變形結果。

第五章 結論與展望

本文以共旋轉有限元素推導法，殼元素來分析膜結構的皺折，由本 文分析之數值例題的結果，可得以下的結論： (1)膜結構邊界受剪力側向挫屈後有很多不同平衡路徑，且膜結構受 力後的變形很敏感，差異很小的負荷條件，有可能會造成差異很大的 變形結果。 (2)由本文例題結果發現膜結構受剪位移時，若元素網格不夠密，則 膜結構皺折的數目太少。若欲得到精確的結果，應使用較密元素網格。 (3)本文對於平面膜結構產生側方向位移的機制為在分歧點加上挫 屈模態擾動，使膜結構有側方向的位移，本研究發現挫屈模態與膜結 構挫屈後皺折的的形狀差異很大。 (4)由本研究分析的結果與文獻上的實驗結果比較，可以發現兩者的 皺折型態有相同的特徵，但細部仍有差異，其原因可能是本研究分析 時，元素的數目不夠，文獻上的結構有初始不完美，且實驗操作上真 實的邊界條件可能與分析的邊界條件不同。 本文例題發現使用太密的元素網格，會出現不易收斂的情形，這可 能是因為本研究以平板三角殼元素來模擬薄膜結構的行為，可能因平 板殼元素沒有考慮撓曲變形與膜變形間的耦合及採用近似的幾何剛 度，造成收斂的困難。也可能是個數值方法上的問題，未來可考慮採 用其他元素或數值方法，若能解決這些困難，使用較密的元素網格， 有可能可以得到更精確的結果。