薄殼結構在位移負荷作用下之幾何非線性分析

(1)

國立交通大學

機械工程學系碩士班

碩士論文

薄殼結構在位移負荷作用之下的幾何非

線性分析

Geometrically Nonlinear Analysis of Thin

Shell Structures under Displacement Type

Loading

研究生：楊禮龍

指導教授：蕭國模博士

(2)

薄殼結構在位移負荷作用下之幾何非線性分析

研究生：楊禮龍指導教授：蕭國模博士國立交通大學機械工程學系碩士班

摘要

本文以共旋轉(co-rotational formulation)有限元素推導法及增量迭代法來探討薄殼結構在位移負荷作用下的幾何非線性行為。本文採用文獻[24]提出的平面三角殼元素，此元素是由 CST(constant strain triangle) 平面元素與 DKT(discrete Kirchhoff theory)三角板元素疊加而成。此元素有三個節點，每一節點具有六個自由度。結構的節點座標、增量位移與增量旋轉、以及結構的平衡方程式都是定義在一組固定的總體座標系統上；而殼元素的應變、節點內力以及元素剛度矩陣，則是在當前元素位置上建立的元素座標上定義。本文採用牛頓-拉福森(Newton-Raphson)法和弧長控制(are length control)法的增量迭代法來解受位移負荷之結構的非線性平衡方程式。最後，本研究以四個數值例題探討殼結構受到各種位移負荷之幾何非線性行為。例題一為鉸接球殼受到不同的側向位移負荷，例題二為圓柱殼受到不同的側向位移負荷，例題三為含缺口之簡支圓柱殼之邊界受到均勻位移及均勻力負荷，例題四為懸臂板受到不同的側向位移負荷。由例題結果可以知道結構受單一位移負荷及力負荷時，其平衡路徑是一樣的，但受多個位移負荷作用與受多個力負荷作用時，結構的行為有很大的差異。

(3)

Geometrically Nonlinear Analysis of Thin Shell Structures under

Displacement Type Loading

Student：Li Lung Yang Advisor：Dr. Kuo Mo Hsiao

Department of Mechanical Engineering National Chiao Tung University

ABSTRACT

The geometrically nonlinear behavior of thin shell structure under displacement loading are investigated using the co-rotational finite element formulation and a incremental-iterative method.

The shell element employed here is the flat three-node triangular shell element with six degrees-of-freedom per node proposed by Bathe and Ho’s [24] .The shell element is obtained by superimposing CST(constant strain triangle)element and DKT(discrete Kirchhoff theory) triangular

plate element. The nodal coordinates, displacements, rotations, and the

equilibrium equations of the structure are defined in a fixed global set of

coordinates. The strains of shell element, the element internal nodal forces,

the element stiffness matrix are defined in terms of element coordinates, which are constructed at the current configuration of the shell element. A incremental-iterative method based on the Newton-Raphson method and constant arc length method is used for solving nonlinear equilibrium equations with displacement loading.

(4)

nonlinear behavior of thin shell structures under different proportional displacement loadings. Example 1 is a hinged spherical shell under different lateral displacement loadings, Example 2 is a cylindrical shell under different lateral displacement loadings, Example 3 is a simply supported cylindrical shell with cutout under uniform in plane displacement loadings and force loadings, Example 4 is a cantilever plate under different lateral displacement loadings. It is found that a single concentrated displacement loading is equivalent to a single concentrated force loading as expected. However, the difference between the structure behaviors correspond to multiple displacement loading and force loading is remarked.

(5)

誌謝

衷心感謝指導教授蕭國模博士兩年來的指導與教誨，使得本論文得以順利完成。在這兩年碩士班期間，讓我印象最深的是老師在研究上嚴謹的態度及對我們這群學生的日常生活上的關懷，都使我受益良多，在此致上最高的謝意。同時也感謝葉孟考老師以及尹慶中老師撥冗擔任口試委員，為本文提出寶貴的意見。感謝陳弘虎、蔡明旭、劉峰成、陳致中、劉宗帆學長的引領照顧，以及學弟黃孝衡、楊水勝在學業上的砥礪與成長。感謝我的父母親、妹妹及所有關心我的朋友對我的支持與鼓勵，僅以此成果與榮耀，獻給我親愛的父母以及所有關心我的人。

(6)

目錄

中文摘要 ……….. I 英文摘要 ……….. II 致謝 ………. IV 目錄 ……….. V 圖目錄 ………. VI 第一章緒論 ……… 1 第二章理論推導 ……….. 3 2.1 基本假設 ………. 3 2.2 座標系統 ………. 3 2.3 旋轉向量 ………….………..……….. 4 2.4 殼元素變形的描述 .………. 4 2.4.1 常應變三角元素(CST)的變形描述 ………. 5 2.4.2 DKT 元素的變形描述 .……… 7 2.5 元素內力與元素剛度矩陣 ……… 10 2.5.1 CST 元素之節點內力與剛度矩陣 ………. 11 2.5.2 DKT 元素的節點內力及剛度矩陣 ……… 12 2.6 元素幾何剛度矩陣 ………... 13 2.7 元素變形角的描述 ………... 14 2.8 系統的平衡方程式與收斂準則 ……….. 16 第三章數值計算方法與程序 ……… 18 3.1 增量迭代法 ……….. 18 3.2 弧長控制法 ……….. 21 3.3 數值程序 ……….. 21 第四章數值例題與結果 .……….….. 23 第五章結論 ….……...……… 30 參考文獻 .……….. 31 附錄 A DKT 元素的形狀函數 ……….. 35 附錄 B CST 元素的剛度矩陣 ……….. 41 附錄 C 分佈反力之節點值與節點反力的關係 ………. 42

(7)

圖目錄

圖 2.1 三角殼元素的示意圖及節點自由度…...………... 44 圖 2.2 總體座標與三種元素座標………. ₄₅ 圖 2.3 旋轉向量…..……….……….. ₄₆ 圖 2.4 DKT 元素的節點及其三邊上的局部座標示意圖…..……… 47 圖 2.5 CST 元素在元素座標上的變形位移………... 48 圖 2.6 變形前板元素中心面之單位法向量受旋轉向量作用的示意圖……….………. n θ 49 圖 2.7 元素座標的剛體旋轉 (a)面外旋轉(out-of plane rotation),

(b)面內旋轉(in-plane rotation)………. 50 圖 2.8 決定板元素節點變形轉角的第 3 個步驟的示意圖....……… 51 圖 4.1 (a)例題一之 Hinged spherical shell 結構示意圖 (b)俯視圖

(c)位移負荷圖……….. 52 圖 4.2 E 點之反力–負荷參數曲線圖(例題一, Case 1)………... 53 圖 4.3 A 點之位移–負荷參數曲線圖(例題一, Case 1)………….... 54 圖 4.4 B 點與 E 點之反力–負荷參數曲線圖(例題一, Case 2)……. 55 圖 4.5 B 點與 E 點之位移–負荷參數曲線圖(例題一, Case 2)……. 56 圖 4.6 B 點與 E 點之反力–負荷參數曲線圖(例題一, Case 3)……. 57 圖 4.7 A 點、B 點、C 點之反力–負荷參數曲線圖(例題一, Case 4) 58 圖 4.8 A 點、B 點、F 點之反力–負荷參數曲線圖(例題一, Case 5) 59 圖 4.9 (a)例題二之 Hinged cylindrical shell 結構示意圖 (b)俯視圖

(c)位移負荷圖……….. 60 圖 4.10 E 點之反力–負荷參數曲線圖(例題二, Case 1)….………. 61 圖 4.11 A 點之位移–負荷參數曲線圖(例題二, Case 1)………….. 62 圖 4.12 B 點與 E 點之反力–負荷參數曲線圖(例題二, Case 2)….. 63 圖 4.13 B 點與 E 點之位移–負荷參數曲線圖(例題二, Case 2)….. 64

(8)

圖 4.14 B 點與 E 點之反力–負荷參數曲線圖(例題二, Case 3)….. 65 圖 4.15 A 點、B 點、C 點之反力–負荷參數曲線圖 (例題二, Case 4)……….. 66 圖 4.16 A 點、B 點、F 點之反力–負荷參數曲線圖 (例題二, Case 5)……….. 67 圖 4.17 B 點、E 點之位移–負荷參數曲線圖(例題二, Case 5)…… 68 圖 4.18 (a)例題三之 Hinged cylindrical shell 結構示意圖

(b)位移負荷之俯視圖………... 69 圖 4.18 (c)位移負荷之前視圖 (d)力負荷之俯視圖 (e)力負荷之前視圖... 70 圖 4.19 A 點之位移W_A–負荷參數λ_D曲線圖 (例題三, Case 1)……... 71 圖 4.20 反力參數λ_R–位移負荷參數λ_D曲線圖 (例題三, Case 1)……….. 72 圖 4.21 位移邊界上之無因次分佈反力圖 (例題三 Case 1 , mesh 10×10)………... 73 圖 4.22 位移邊界上之無因次分佈反力圖 (例題三 Case 1 , mesh 20×20)………... 74 圖 4.23 位移邊界上之無因次分佈反力圖 (例題三 Case 1 , mesh 25×25)……….. 75 圖 4.24 位移邊界上之無因次分佈反力圖 (例題三 Case 1 , mesh 30×30)……….. 76 圖 4.25 不同 Mesh 的位移邊界上之無因次分佈反力圖 (例題三 Case 1 , λ_D =0.266mm_{)……….. 77} 圖 4.26 W(0)在不同 Mesh 的數值比較圖 (例題三 Case 1 , 位移控制)………. 78 圖 4.27 位移邊界上之無因次分佈反力圖 (例題三 Case 2 , mesh 10×10)………... 79 圖 4.28 位移邊界上之無因次分佈反力圖

(9)

(例題三 Case 3 , mesh 10×10)………... 80 圖 4.29 位移邊界上之無因次分佈反力圖 (例題三 Case 4 , mesh 10×10)………... 81 圖 4.30 A 點之位移W_A–位移負荷參數λ_D曲線圖 (例題三, 位移負荷)……… 82 圖 4.31 F 點之位移W_F–位移負荷參數λ_D曲線圖 (例題三, 位移負荷)……… 83 圖 4.32 A 點之位移W_A–反力參數λ_R曲線圖 (例題三, 位移負荷)……… 84 圖 4.33 A 點之位移W_A–力負荷參數λ_F曲線圖 (例題三, 力負荷)……… 85 圖 4.34 A 點之位移W_A–力邊界上軸方向位移平均數X U_avg 曲線圖(例題三, 力負荷)……… 86 圖 4.35 力邊界上軸方向之無因次位移分佈圖 X (例題三 Case 1 , mesh 10×10)………... 87 圖 4.36 力邊界上軸方向之無因次位移分佈圖 X (例題三 Case 2 , mesh 10×10)……….. 88 圖 4.37 力邊界上軸方向之無因次位移分佈圖 X (例題三 Case 3 , mesh 10×10)……….. 89 圖 4.38 力邊界上軸方向之無因次位移分佈圖 X (例題三 Case 4 , mesh 10×10)……….. 90 圖 4.39 受位移負荷之λ_R–λ_D曲線圖及受力負荷之λ_F–U_avg 曲線圖(例題三)………... 91 圖 4.40 (a)Cantilever plate 結構示意圖 (b)位移負荷示意圖……... 92 圖 4.41 A 點之反力–位移曲線圖(例題四, Case 1)……….. 93 圖 4.42 A 點與 B 點之反力–負荷參數曲線圖(例題四, Case 2)….. 94 圖 4.43 A 點與 B 點之位移–負荷參數曲線圖(例題四, Case 2)... 95 圖 4.44 A 點與 B 點之位移–負荷參數曲線圖(例題四, Case 2)... 96 圖 4.45 C 點之位移–負荷參數曲線圖(例題四, Case 2)... 97

(10)

圖 4.46 A 點與 B 點之反力–負荷參數曲線圖(例題四, Case 3)….. 98 圖 4.47 A 點與 B 點之位移–負荷參數曲線圖(例題四, Case 3)…. 99 圖 4.48 C 點之位移–負荷參數曲線圖(例題四, Case 3)...100

(11)

第一章緒論

近代科學的發展，對於材料的要求趨向於高強度以及低重量，然而任何結構物的設計，不管是建築結構物、壓力容器、航空設備等等，設計者必須完成以最小成本達到其基本機能與美學要求為目的，由於薄殼是最經濟的結構之一，同時亦可承受大型彎曲應力與拉伸應力而可獲得最佳結構效果。常見的薄殼結構包括：薄殼屋頂、建築圓頂、飛機機身蒙皮結構、引擎葉片、壓力油管或水中潛體之外殼等。薄殼結構常需承受大位移及大旋轉，即使在彈性範圍內，其受力與變形的關係常常呈現非線性變化，所以必須使用非線性分析來探討由幾何形狀改變所造成的非線性行為，以求得較可靠的的結構行為。結構所受的負荷一般可以分為兩種，一種為力負荷(force-loading)，另一種為位移負荷(displacement-loading)[1]。文獻上有很多探討殼結構在力負荷作用下的幾何非線性行為的論文[2-24]，但甚少有探討結構在位移負荷作用下的幾何非線性行為的論文[25]。據作者所知，僅有文獻[25]探討過平面梁結構受位移負荷的幾何非線性行為，本文的主要目的為探討殼元素構成的結構在位移負荷下的幾何非線性行為，以補足文獻上的不足。文獻中有很多殼元素被提出[2-7,9,12-20,22-24,26-30,32]，文獻[24] 中將 CST(constant strain triangle)平面元素[31] 與 DKT(discrete

Kirchhoff theory)三角板元素[32]疊加成一平面三角殼元素，並使用更新拉格蘭日法(updated Lagrangian formulation)將該元素用在具大位移及大旋轉的薄殼結構分析，由[24]的數值例題可見該元素相當的簡單及精確。文獻[24]中在計算對應於 DKT 元素部分的內力時是將對應於每一増量位移的増量內力相加而成，但[24]中在計算増量內力時將旋轉視為向量，且沒有扣除剛體旋轉的部分。因三維有限旋轉並非向量，

(12)

所以[24]的方法在兩増量間的增量旋轉必須是小角度，才能得到精確的答案。為了克服此困難，文獻[33]採用共旋轉法來描述殼元素的變形，並提出一運動過程來決定平面三角殼元素的節點總變形角，再用總變形角來計算 DKT 元素的節點內力。因[33]的方法可以有效的解決兩增量間增量轉角的限制，故本研究擬採用[33]中的共旋轉法及[24] 的平面三角殼元素。結構在受到力負荷及位移負荷時位移-負荷的平衡路徑有時會出現 snap-back 與 snap-thoughs 的情況[34]，文獻[35]Crisfield 提出基於牛頓 -拉福森(Newton-Raphson)法的弧長控制(are length control)法求解在力負荷作用下的非線性平衡方程式，可以克服 snap-thoughs 與 snap-back 的情況，求得完整的平衡路徑。但是文獻[35]提出的方法有不能直接用在位移負荷的問題，所以文獻[25]修改了[35]的弧長控制法且提出了一個牛頓-拉福森(Newton-Raphson)法和弧長控制(are length control)法的增量疊代法來解受位移負荷的非線性平衡方程式，並將其用在平面梁結構受位移負荷的幾何非線性分析，因為文獻[25]提出的方法非常有效，故本文採用文獻[25]的方法來解非線性平衡方程式。最後，本研究以四個數值例題探討殼結構受到各種位移負荷之幾何非線性行為。例題一為鉸接球殼受到不同位移負荷作用的幾何非線性行為；例題二為圓柱殼受到不同位移負荷作用；例題三為含缺口的簡支圓柱殼之邊界受到均勻位移及均勻力負荷作用的幾何非線性行為；例題四為懸臂板受到不同的位移負荷作用。由例題結果可以知道結構受單一位移負荷及力負荷時，其平衡路徑是一樣的，但受多個位移負荷作用與受多個力負荷作用時，結構的行為有很大的差異。

(13)

第二章理論推導

本文使用文獻[33]的共旋轉法(co-rotational formulation)，與文獻[33] 一樣，本文章採用文獻[24]中提出的平面三角殼元素，如圖 2.1 所示，文獻[24]中的殼元素是由 CST(constant strain triangle)常應變三角形元素[31]及在文獻[32]中提出的 DKT(discrete Kirchhoff theory)三角殼元素所疊加而成。為了本文的完整性，在本章將簡單描述文獻[24]中元素變形的假設及內力、剛度矩陣的推導。文章中亦將簡單描述文獻[33] 中決定元素變形角的方法。 2.1 基本假設在文獻[24]中對於非線性平面三角殼元素的推導，做以下的假設 (1) 薄膜變形(membrane deformation)以及彎曲變形(bending deformation)之間無耦合作用。 (2) 殼元素的變形為小變形。 (3) 在元素變形前，垂直於中心面的法向線段，在元素變形後，依然保持直線，且沒有伸長及縮短，除了在元素三個頂點以及三個邊的中央點外，該線段不必垂直於元素變形後的中心面。 2.2 座標系統本文採用共旋轉法，為了描述系統的運動以及元素的變形，在圖 2.2 中本文定義了兩組座標系統 (1) 固定總體座標系統(global coordinate):X_i(i = 1,2,3)；結構體所有結點的座標、增量位移、增量旋轉、結構的平衡方程式均在此座標系統中定義。 (2) 元素座標系統(element coordinate):x_i (i = 1,2,3)；此座標系統是

(14)

建立在每一殼元素變形後的最新位置上，其座標原點是在元素節點 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′之間的關係可表示成[36] ) ( 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。本殼元素的變形分為薄膜變形(membrane

deformation)與彎曲變形(bending deformation)兩種，假設元素的薄膜變

形與彎曲變形之間無耦合作用，所以元素的變形可由薄膜變形及彎曲

變形疊加而成。本元素的變形是由CST (constant strain triangle)常應變

三角形元素[31]的薄膜變形及文獻[32]中的DKT(discrete Kirchhoff

theory)三角形殼元素的彎曲變形疊加而成。

(15)

移，而

θ

_xj、

θ

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

zj

θ

是為了不使元素剛度內的面內旋轉剛度(in-plane rotational stiffness)

為 0，而人為加上去的自由度。圖 2.4 中的元素是在[32]中所提出的 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 ) 6 1 (i = − a_i 由元素座標的定義可知，在元素座標中圖2.1 的節點1、節點2、節點3之變形前座標值可表示成(0,0)、(x₂,0)、(x₃,y₃)，令節點 j在x₁ 軸、x 軸的位移分別是₂ u_j ,v_j (j =1,2,3)。將三個節點的座標值及節點位移u_j ,v_j (j =1,2,3)代入(2.2)、(2.3) 可以得到 (2.4) ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ 3 2 1 3 3 2 3 2 1 1 0 1 0 0 1 a a a y x x u u u

(16)

⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ 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)表示成及的函數，所以(2.2)、(2.3) 式可改寫成 i u v_i u= Nu_m (2.6) u={u ,v} (2.7) _⎥ (2.8) ⎦ ⎤ ⎢ ⎣ ⎡ = 3 2 1 3 2 1 0 0 0 0 0 0 N N N N N N N u_m =

{

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

}

(2.9) 1 ( ₂ ₃ ₃ ₂ ₃ ) 3 2 1 x y xy x y x y y x N = − − + 1 ( ₃ ₃ ) 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可表示成

(17)

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 元素的變形描述

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

θ

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

θ

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

θ

_y

(18)

與θ_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為法線方向。在[32]中對於其所提出的DKT元素做了下列的假設 (1)

θ

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

∑

, (2.22) = = 6 1 i yi i y N

θ

_∑

= − = 6 1 i xi i x N

θ

其中

θ

_yi、θ_xi是

θ

_y、θ_x在圖2.4 中節點i的節點值，為形狀函數[32]，其表示式詳見附錄A。 ) 6 1 (i = − N_i (2) 元素的頂點以及三個邊的中央點可滿足克希霍夫板理(Kirchhoff plate theory)的假設，即

(19)

(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_j是在節點w i及 j的值，w_,_sk是w_,_s在節點k 的值，k =4,5,6 分別為邊 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之幾何轉換關係可表示成

(20)

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ 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)式表示成[32] 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) 其中為DKT元素的節點位移，與是對應於元素節點位移的新形狀函數，其表示式詳見附錄A， b u H_x H_y

ξ

與

η

是元素內任一點在元素自然座標[32]的座標值，其中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 元素內力與元素剛度矩陣

(21)

本文中殼元素的節點內力是由 CST 及 DKT 元素的節點內力組成而成，元素剛度矩陣是由 CST 元素剛度矩陣、DKT 元素剛度矩陣以及面內旋轉剛度所疊加而成。本文中的值是取之對角線元素的絕對值中的最小值。本節中將用虛功原理推導 CST 元素及 DKT 元素的節點內力及剛度矩陣。 m k k_b z k_θ k_θ_z k_b 在平面應力狀態，等向性線彈性材料的應變與應力關係為 σ=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)式中的ε ε_m及(2.19)式中的ε_b。 2.5.1 CST 元素之節點內力與剛度矩陣將(2.33)式代入(2.13)式可得 σ_m = EB_mu_m (2.37) 由虛功原理可得 =

∫

(2.38) V m t m m t mf ε σ dV u δ δ } { _m₁ _m₂ _m₃ m f f f f =

(22)

f_mj ={f_xj f_yj} j =1,2,3 (2.39) 其中f_m是 CST 元素對應於δu_m的節點內力，V 是元素的體積，定義於(2.11)式。 m u 將(2.13)式、(2.37)式代入(2.38)式可得 =

∫

(2.40) V m m t m t m m t mf u B EB u dV u δ δ 由(2.40)式可得 f_m =k_mu_m (2.41) =

∫

(2.42) V m t m m B EB dV k 其中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所做的虛功。由虛功原理可得 (2.44)

∫

= V b t b b t bf ε σ dV u δ δ f_b ={f_b₁ f_b₂ f_b₃} f_bj ={f_zj m_xj m_yj} j =1,2,3 (2.45) 其中f_b是 DKT 元素對應於δu_b的節點內力，u_b定義於(2.30)式，V 為 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 δ δ

(23)

=

∫

(2.46) A b b b t b t b B D B u dA u δ 其中 ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − = =

_∫

₋ 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)式可得[32] f_b =k_bu_b (2.48)

∫ ∫

− = 1 0 1 0 2A η t_b _b _b dξdη b B D B k (2.49) 其中k_b是 DKT 元素剛度矩陣。 2.6 元素幾何剛度矩陣元素的幾何剛度對平衡迭代時的收斂速度有很大的影響，為了改善收斂速度，本文中在平衡迭代的過程中將元素幾何剛度矩陣加入元素剛度矩陣中。本文中採用文獻[24]中的元素幾何剛度，其表示式為 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

(24)

其中 x₃₂ = x₃ − x₂ _⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 2 12 12 1 * N N N N N N =hEB_mu_m 其中 h 為板元素厚度，E定義於(2.36)式，定義於(2.15)式，定義於(2.9)式。 m B u_m (2.50)式為一近似的幾何剛度矩陣，僅考慮面內應力σ_x、

σ

_y以及

τ

_xy 在剛體運動時的效應。由本文分析的例題可發現將加入剛度矩陣可以加速平衡迭代的收斂速率。 g k 2.7 元素變形角的描述本文中採用文獻[33]中的兩階段旋轉法來決定元素的剛體旋轉及節點在元素座標的變形位移與轉角。假設第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得到。由可以利用元素座標定義求出當前的元素座標 j X x_i (i =1, 2,3)，本文中稱元素從I _i x 到x 的運動為對應於_i

∆

U_j、

∆

Φ_j (j =1,2,3)之剛體運動，本文中假設該剛體運動是由以下三個步驟達成的。 (1) ∆U₁造成的位移：元素座標Ix 的原點受到_i ∆U₁的作用移動到元素座標x 的原點，如圖_i 2.7(a)，其中∆U₁是作用在元素節點 1的元素節點增量位移向量。

(25)

(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 e Ix₃軸與x 軸上的單位向量。 ₃ 因為旋轉向量造成的運動是朝遠離α x₁ −x₂平面方向，所以α 所造 成的旋轉稱為面外旋轉。 (3) 旋轉向量 β 造成的面內旋轉(in-plane roration)：元素的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₂平面，所以所造成的旋轉稱為面內旋轉。 β 本文使用文獻[33]提出的直接法(direct method)來計算對應於

∆

U_j及 j Φ

∆

(j =1, 2,3)的元素變形角。在文獻[33]所提出的直接法可以分成以下四個步驟。 (1) ∆U₁造成的剛體移動：整個元素因為∆ U₁而移動，其中∆ U₁是作用在元素節點1的元素節點增量位移向量。 (2) 旋轉向量α造成的剛體旋轉：旋轉向量α (2.51式)通過元素節點1 使得整個元素除了In_dj之外由I _i x (i =1, 2,3)座標轉到x′_i (i =1, 2,3) 座標，其中In_dj是在垂直於在第I次疊代收斂後位在元素節點 j的變形

(26)

後元素中心面的法線向量。 (3)

∆

Φ_tj造成In_dj的有限旋轉：法線向量In_dj受到

∆

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

∆

Φ_tj是增量旋轉位移向量

∆

Φ_j 在x₁′− x₂′平面的投影向量。 (4) 旋轉向量 β 造成的剛體旋轉：旋轉向量 (2.52 式)通過元素節點 1 使得整個元素包括及旋轉到最後的位置及，其中是受到α 作用後的新位置， 是在第 β u n′ n′_dj n_u n_dj n′_u I_n_u u I_n I次疊代收斂後，垂直於未變形的元素中心面的法線向量，是當前垂直於未變形元素中心面的法線向量。 u n 令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為系統位移向量，

λ

為負荷參數，為參考位移負荷向量。_F可由(2.39)與(2.45)式的元素節點力，用標準的座標轉換，轉換到對應總體座標後再組合而成。 P Q

本文以不平衡力的weighted Euclidean norm 作為迭代時的誤差度

量，而且收斂準則表示為

(27)

tol e e e= ≤ F Ψ (2.55) 其中為對應於位移負荷的系統節點反力向量。為一給定之容許誤差值。 e F e_tol

(28)

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

本文解非線性平衡方程式(2.54)式的數值計算方法是採用文獻[25] 中所提出基於牛頓-拉福森(Newton-Raphson)法和弧長控制(arc length control)法的增量疊代法。為了文章的完整性，本章中將簡單介紹文獻 [25]中提出的數值計算的方法與程序。 3.1 增量迭代法在增量迭代法中，若第 I 個增量的平衡位置為已知，令其位移向量為Q_I、負荷參數為λ_I，則第 I+1 個增量的初始增量位移向量

∆

Q，可利用尤拉預測值( Euler predictor )求得[34] ∆Q = ∆λr_T ， (3.1) r_T K_TI 1F_P， (3.2) ) ( − − =

∂λ

∂

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

∆

λ

為初始增量負荷參數，為(2.54)式中的系統節點內力，F 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)

(29)

_G G G r f k

∂

= (3.6)

∂λ

∂

G G P r r = (3.7) 其中，k 為在總體座標上的元素剛度矩陣，G _{k 是由在元素座標中的}G 元素剛度矩陣經標準轉換得到，為元素的節點參考位移負荷，可以由系統參考位移負荷萃取而得到。 G P r P Q (3.1)式中的

∆

λ

可利用下式求出[35]

( )

12 T t T l r r

∆

λ

∆

=± ， (3.8) 其中正負符號之決定方法如下：若第 I-1 與 I 個增量收斂時，其系統切線剛度矩陣之行列式值同號，則

∆

λ

的正負符號和第 I 個增量時相同； 若異號則符號相反。 l

∆

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

(

)

12 I D I J J l l

∆

= ， (3.9) 其中為給定的期望迭代次數，為第 I 個增量迭代至平衡所使用的迭代次數， D J J_I I l ∆ 為第 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₀在參考自由度的分量之絕對值。

(30)

的平衡迭代：由Q_{I 1}₊ =Q_I +∆Q，λ_{I 1}₊ =λ_I +∆λ，參考位移負荷向量及 2.4 節與 2.6 節的方法，則可以求得系統中各殼元素在當前的元素座標、節點變形位移、節點變形角。再利用 2.5 節的方法，求得在當前元素座標上的節點內力及剛度矩陣再由標準轉換轉到總體座標的節點內力及剛度矩陣。由(2.54)式此時之不平衡力可表示成 P Q _I

ٛ

₊₁ =F

(

Q_I₊₁,λ_I₊₁Q_P

)

， (3.11) 將(2.54)式在λ =λ_I₊₁，Q=Q_I₊₁時用泰勒展開式展開為 = ₊ + +

δλ

+

λ

∂

δ

∂

F Q Q F Ψ Ψ

ٛ

_I ₁ (二次以上高次項) =Ψ_I₊₁ +K_Tr+δλF_P +(二次以上高次項) (3.12) 其中， Q F K

∂

= T 為系統切線剛度矩陣，r =

δ

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

δλ

F_P)=r₀ +

δλ

r_T ， (3.13) 其中，r₀ =−K_T−1

ٛ

Ψ_I₊₁，r_T在( 3.2 )式中已定義。 (3.13)式中的位移負荷參數修正量

δλ

有很多不同的決定方法[37]，如 Crisfield 在文獻[35]所提出的定弧長控制法及 Fried 在文獻[38]提出的正交軌線法，本文是使用 Crisfield[35]提出的定弧長控制法決定

δλ

，並將在 3.2 節中介紹。將求得之位移負荷參數修正量

δλ

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

∆

λ

與

∆

Q中，可得到一新的增量位移向量與增量位移負荷參數，再進行下一次迭代。以上之迭代過程一直重複至(3.11)式中的不平衡力滿足(2.55)式的收斂準則為止。

(31)

3.2 弧長控制法 (3.13)式中的位移負荷參數修正量

δλ

可利用文獻[35]中 Crisfield 所提出的定弧長控制法決定，其方法使在每一個增量中固定其增量位移 的 Euclidean norm 為一定值 l

∆

。由新的增量位移向量為

(

∆

Q+r

)

可得

∆

l2 =

(

∆

Q+r₀ +

δλ

r_T

) (

t

∆

Q+r₀ +

δλ

r_T

)

， (3.14) 上式經過整理後可以得到δλ的二次方程式 a₁

δλ

2 +a₂

δλ

+a₃ =0， (3.15) 其中 a₁ =r_Ttr_T， a₂ =2

(

∆

Q+r₀

)

tr_T ， a₃ =

(

∆

Q+r₀

) (

t

∆

Q+r₀

)

−

∆

l2，當(3.15)式的解為兩實根δλ₁，δλ₂時，本文中取兩根中能使新的增量位移

(

∆

Q+r

)

與前次迭代增量位移

∆

Q間的內積(inner product)較大者為

δλ

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

∆

減小，重做第 I+1 個增量。 3.3 數值程序本文使用的增量迭代法之數值程序可以分為三個主要部分 1. 輸入與計算開始分析所需的資料 (a) 輸入結構與負荷資料。 (b) 選擇一個參考自由度，並給定期望此自由度達到的最大位移。 (c) 給定最大增量數、每個增量期望的迭代數與最大的迭代數、收斂時的容許誤差。

(32)

(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.1 節的方法計算增量位移修正量與增量負荷參數修正量，然後回到步驟(a)。 3. 計算下一次增量所需的資料 (a) 檢查參考自由度的位移或已進行的增量次數是否已達給定值。若已達到給定值則停止分析工作；否則進行下一步工作。 (b) 計算(3.12)式中的切線剛度矩陣與F_P。 (c) 以(3.1)-(3.3)式與(3.8)-(3.9)式計算下一次增量的增量位移、增量弧長與增量負荷參數。 (d) 回到第二部份執行迭代工作。

(33)

第四章數值例題與結果

在本章中，將以第二章之梁元素及第三章所提的幾何非線性分析之數值程序與方法，探討不同的板殼結構在位移負荷作用下的非線性行為，並說明本文數值方法與程式之可行性與正確性。

例題一：鉸接球殼受到不同位移負荷作用

圖 4.1 所示為本例題分析之 Hinged spherical shell 的結構圖及其所受之位移負荷的示意圖，該結構四周為鉸接，圖 4.1(b) 中的 mm 75 . 159 d = 。本例題考慮的位移負荷有下列五種情形：(Case 1) E 點受一向下(-Z方向)的位移負荷λ_E作用，其中λ_E =λ (mm)，λ 是負

荷參數(loading parameter)，(Case 2) B 點、E 點受一向下( 方向)的位

移負荷 Z -E B ,λ λ 作用，其中λ_E =λ (mm) λ_B =0.5λ (mm)，(Case 3) B 點、 E 點受一向下 ( -Z 方向 ) 的位移負荷 λ_B ,λ_E 作用，其中 ) mm ( B E λ λ λ = = ，(Case 4) A 點、B 點、C 點受一向下( 方向)的位移負荷 Z -C B A ,λ ,λ λ 作用，其中λ_A =λ_B =λ_C =λ (mm)，(Case 5) A 點、 B 點、F 點受一向下(-Z方向)的位移負荷λ_A ,λ_B ,λ_F作用，其中 ) mm ( F B A λ λ λ λ = = = 。本例題的 Case 1 與文獻[33]中在相同殼結構之 E 點施加向下的力負荷是等效的。Case 1 由於雙對稱之故，本文僅取分析 1/4 的結構並將其離散成 50 個元素(5

×

5 網格)。除了 Case 1 外，本例題以整個結構進行分析，並將此結構離散成 200 個元素(10

×

10 網格)，本例題之平衡迭代的容許誤差值取1×10−5，Case 1 的分析過程中使用了 8 個增量，每個增量的平均迭代次數約為 4；Case 2 的分析過程中使用了 18 個增量，每個增量的平均迭代次數約為 4；Case 3 的分析過程中使用了 18 個增量，每個增量的平均迭代次數約為 5；Case 4 的分析過程中使用了 19 個增量，每個增量的平均迭代次數約為 5；

(34)

Case 5 的分析過程中使用了 18 個增量，每個增量的平均迭代次數約為 5。圖 4.2-圖 4.8 為本例題的結果，圖 4.2 為 Case 1 的結果與文獻[33] 比較，圖中R_E為在位移負荷方向的反力，R_E及

λ

相當於文獻[33]之力負荷及力負荷方向的位移，由圖 4.2 可見本文的結果與文獻[33]的結果完全重合，由圖 4.2 的結果可以得知施加單一位移負荷與施加單一力負荷是等效的。圖 4.3 為 Case 1 中 A 點的位移–負荷參數曲線圖，圖 4.4 為 Case 2 中 B 點與 E 點之反力–負荷參數曲線圖，由圖 4.4 可以發現 B 點的反力為負值，即反力的方向與設點位移的方向相反。由圖 4.3 可知當 Case 1 中 E 點受位移負荷

λ

時 A 點的位移大於0.5

λ

，因為對稱的關係所以 B 點的位移等於 A 點也大於0.5

λ

，在 Case 2 中 λ λ λ λ_E = , _B =0.5 ，所以需有一與λ_B 方向相反的反力，才能使 λ λ_B =0.5 。圖 4.5 為 B 點與 E 點的位移–負荷參數曲線圖，因 Case 2 的位移負荷對 E 點不對稱，所以不為。圖 4.6-4.8 為 Case 3-5 之反力–負荷參數曲線圖，因 Case 3 的 E U 0 λ λ λ_E = _B = ，所以圖 4.6 中的值是正的，即 B 點的反力與其位移同向，由圖 4.7 中可以發現與 –負荷參數的曲線重疊，這是因為 A 和 C 點的位移負荷對通過 EB 之大圓(Great Circle)的平面是對稱的。由圖 4.2、4.4、4.6-4.8 可以發現不同的位移負荷造成的反力有很大的差異。 B R A R C R 例題二：圓柱殼受到不同位移負荷作用

圖 4.9 中所示為本例題分析之 Hinged cylindrical shell 之結構圖及其

所受之位移負荷圖，此結構與軸平行的邊界為鉸接，而與軸平行的邊界為自由邊，圖 4.9(b)中的 Y X mm b mm a=50.71 , =50.8 。與例題一一樣，本例題考慮的位移負荷有下列五種情形：(Case 1) E 點受一個向下(-Z方向)的位移負荷λ_E作用，其中λ_E =λ (mm)，λ 是負荷參數

(35)

(loading parameter)。(Case 2) B 點、E 點受一向下(-Z方向)的位移負荷 E B ,λ λ 作用，其中λ_E = λ (mm)、λ_B =0.5λ (mm)。(Case 3) B 點、E 點受一向下 ( -Z 方向 ) 的位移負荷 λ_B ,λ_E 作用，其中 ) mm ( B E λ λ λ = = 。(Case 4)分別在 A 點、B 點、C 點受一向下( 方向)位移負荷 Z -C B A ,λ ,λ λ 作用，其中λ_A =λ_B = λ_C =λ (mm)。(Case 5) A 點、B 點、F 點受一向下(-Z方向)的位移負荷λ_A ,λ_B ,λ_F作用，其中 ) mm ( F B A λ λ λ λ = = = 。本例題的 Case 1 與文獻[33]中在相同殼結構之 E 點施加向下的力負荷是等效的。Case 1 由於雙對稱之故，本文僅取分析 1/4 的結構並將其離散成 50 個元素(5

×

5 網格)。除了 Case 1 外，本例題以整個結構進行分析，並將此結構離散成 200 個元素(10

×

10 網格)，本例題之平衡迭代的容許誤差值取1 4 10− × 。Case 1 的分析過程中使用了 7 個增量，每個增量的平均迭代次數約為 6；Case 2 的分析過程中使用了 16 個增量，每個增量的平均迭代次數約為 4；Case 3 的分析過程中使用了 12 個增量，每個增量的平均迭代次數約為 5；Case 4 的分析過程中使用了 14 個增量，每個增量的平均迭代次數約為 5； Case 5 的分析過程中使用了 17 個增量，每個增量的平均迭代次數約為 6。圖 4.10-4.11 為本例題的結果，圖 4.10-4.11 為 Case 1 中 E 點反力– 負荷參數曲線圖與 A 點之位移–負荷參數曲線圖，圖 4.10 中為在位移負荷方向的反力，及 E R E R

λ

相當於文獻[33]之力負荷及力負荷方向的位移，由圖 4.10 可見本文的結果與文獻[33]的結果完全重合。圖 4.12-4.13 為 Case 2 中 B 點與 E 點之反力及位移–負荷參數曲線圖，圖 4.14-4.16 為 Case 3-5 反力–負荷參數曲線圖，圖 4.17 為 Case 5 中 B 點與 E 點之位移–負荷參數曲線圖。例題三：含缺口的簡支圓柱殼邊界受到均勻位移及力負荷作用

(36)

圖 4.18 中所示為本例題分析之 Cylindrical shell 結構與受負荷圖，如圖 4.18 所示，本例題的圓柱殼在中央有一缺口 BCDE，本例題考慮了 4 種缺口的大小， Case 1 ： a =0,b=0 ，即無缺口。 Case 2 ：，Case 3： mm b mm a=25.39 , =25.4 a=25.39mm,b=50.8mm，Case 4：。圖 4.18 與軸平行之右邊邊界為鉸接，僅能在Y軸方向旋轉，左邊邊界為滾子，僅能在方向平移及在軸方向旋轉，而與軸平行的邊界為自由邊。本例題考慮了兩種負荷情形，第一種為在左邊邊界加一均勻的位移負荷 mm b mm a =50.78 , =50.8 Y X Y X D λ (圖 4.18(b)(c))，第二種為左邊邊界加一均勻的力負荷λ_F(圖 4.18(d)(e))。本例題平衡迭代的容許誤差值取1×10−4，用 10

×

10 的網格分析時，Case 1 在分析過程中使用 11 個增量，每個增量的平均迭代次數為 3；Case 2 在分析過程中使用 11 個增量，每個增量的平均迭代次數為 3；Case 3 在分析過程中使用 19 個增量，每個增量的平均迭代次數為 4；Case 4 在分析過程中使用 19 個增量，每個增量的平均迭代次數為 4。圖 4.19-4.38 為本例題的結果，圖 4.19 為 Case 1 中 A 點之位移

W

_A–負荷參數λ_D曲線圖，圖 4.20 為反力參數λ_R–位移負荷參數λ_D曲線圖，圖 4.20 中之反力參數λ_R =

∑

R_i/2L，

∑

為圖 4.18 中邊界 I J 之節點反力的和，為 I J 的邊長。由圖 4.19 及圖 4.20 可發現用 10 i R 2L

×

10 的網格即可求得準確的及 A W λ_R。在位移邊界上的反力為分佈反力，令表示分佈反力，若假設兩相鄰節點間的分佈反力為線性變化，則利用邊界節點的反力可以求出分佈反力的節點值(詳見附錄 C)。圖 4.21-4.24 為 Case 1 中使用不同元素網格時位移邊界上之無因次分佈反力圖，圖 4.25 為負荷參數 ) Y ( Wp mm 266 . 0 D = λ 不同元素網格之無因次分佈反力圖，圖 4.26 為在不同元素網格時的比較圖。由圖 4.21-4.25 可以發現在均勻位 ) 0 ( Wp

(37)

移邊界上的分佈反力並不均勻，Wp 在 0 L Y = 附近有應力集中的現象，但即使使用 100

×

100 的網格，仍未收斂到令人滿意的結果，這可能是因為在附近之分佈反力的變化率太大造成的。當 ) 0 ( Wp 0 X= λ_D較小時，Wp在 1 L Y = 有局部的最大值，隨著 L Y 減少，亦慢慢減小，但當 Wp D λ 逐漸增加時，Wp 的局部最大值逐漸移到 0.5 L Y = 附近。由圖中亦可發現用 10

×

10 的網格求得之Wp(Y)在 1 L Y 2 . 0 ≤ ≤ 時已足夠準確，為了節省計算時間，本研究中 Case 2-4 皆使用 10×10 的網格來分析。圖 4.27-4.29 為 Case 2-4 的位移邊界上之無因次分佈反力圖，由圖中可以發現缺口愈大Wp(L)的值愈小，且其最大值發生在 0.5 L Y = 附近。圖 4.30 為結構受到位移負荷作用時 A 點的位移W_A–位移負荷參數λ_D曲線圖，圖 4.31 為結構受到位移負荷作用時 F 點的側向位移 –位移負荷參數 F W D λ 曲線圖，其中 F 點是指 Case 4 中的 F 點(如圖 4.18(b)(d) 所示)。圖 4.32 為 A 點的位移W_A–反力參數λ_R曲線圖，由圖 4.30 可以發現 Case 1-4 之W_A–λ_D的曲線都重合，即 A 點的側向位移似乎與缺口的大小無關，會造成此一結果的原因應為缺口最大的長與寬也僅佔整個 mesh 的 51/ ，而 A 點距缺口的最小距離為缺口最大長度的 4 倍，故缺口的大小對 A 點的側向位移的影響是很小的。由圖 4.31 可發現當邊界位移相同時，缺口愈大則側向位移愈大；由圖 4.32 可發現當邊界位移相同時，缺口愈大則對應的反力愈小。本例題受力負荷的分析皆是使用 10×10 mesh 來作分析，圖 4.33-4.34 為結構受到力負荷作用時 A 點的位移 –力負荷參數 F W A W λ_F曲線圖以及力邊界上軸方向的節點的位移平均值曲線圖，由圖 4.33 可以發現缺口愈大時， X avg U

(38)

相同的λ_F會造成較大的，由圖 4.34 可發現相同時，的大小與缺口的大小無關。圖 4.35-4.38 為結構受力負荷時 Case 1-4 之力邊界上X軸方向之無因次位移分佈圖，由圖中可以發現對應於均勻分佈外力的位移並不均勻。當沒有缺口時，方向的位移在 A W U_avg W_A X 1 L Y = 時有最小值，當缺口愈大時，在 1 L Y = 的位移逐漸增加。隨著λ_F的增加，的分佈形態也逐漸改變，之最大值及最小值的位置也隨著改變。由圖中亦可發現和之差的最大值約為的，見圖 4.38。圖 4.39 為 Case 1-4 中受位移負荷之 i U i U i U U_avg U_avg 1.5% R λ –λ_D曲線圖及受力負荷之λ_F– 曲線圖，由圖中可以發現在相同的 Case 的 avg U λ_R–λ_D及λ_F– 的曲線幾乎重合。 avg U 例題四：懸臂板受到不同的位移負荷作用在圖 4.40 中所示為本例題分析之 Cantilever plate 的結構與負荷圖，本例題是用 6×8 的網格共 96 個元素來作分析。本例題考慮了三種位移負荷：(Case 1) 只在 A 點加一個向下( 方向)位移負荷，(Case 2) 分別在 A 點及 B 點加一個向下 ( 方向 ) 的位移負荷，其中 Z -Z -λ λ λ λ_A = , _B =− ，(Case 3) 分別在 A 點及 B 點加一個向下( 方向) 的位移負荷，其中 Z -λ λ λ λ_A = , _B = 0.5 。本例題平衡迭代的容許誤差值取，Case 1 在分析過程中使用 8 個增量，每個增量的平均迭代次數為 5；Case 2 在分析過程中使用 19 個增量，每個增量的平均迭代次數為 4；Case 3 在分析過程中使用 13 個增量，每個增量的平均迭代次數為 5。圖 4.41-4.48 為本例題的結果。圖 4.41 為 Case 1 的反力– 位移負荷參數曲線圖。因 4 10 1× − 2 b Y= 為本例題結構的對稱面，亦為 Case 2 的位移負荷之反對稱面，所以可以預期在 Case 2 的位移負荷下，有一

(39)

平衡路徑之 A 點反力和 B 點反力的大小相同方向相反，A 點和 B 點在方向的位移相同，在方向的位移大小相同方向相反，在 A R R_B X Y 0 Z , 2 b Y= = 的直線上，Y和方向的位移為 0。圖 4.42-4.44 為 Case 2 中 A 點與 B 點的反力及位移–負荷參數曲線圖，圖 4.45 為 Case 2 中 C 點之位移–負荷參數曲線圖。由圖 4.42-4.45 可以發現，當負荷參數 Z cm 5 . 3 <

λ

時 A、B 兩點的位移、反力及 C 點的位移應為上述預期的平衡路徑之解，但當

λ

>3.5cm時 A、B 點的位移、反力 C 點的位移顯然為另一平衡路徑的解，這可能是在

λ

=3.5cm附近有分歧點 (Bifurcation point)，

λ

<3.5cm時為主平衡路徑的解，

λ

>3.5cm時為次要平衡路徑的解，但因本文的元素切線剛度僅為一近似剛度，故不能用系統切線剛度的行列式值為零來偵測分歧點。為了求得分歧點後不穩定的主要平衡路徑，本文將 Case 2 在 ,Z 0 2 b Y = = 之直線上在方向的位移令為，並稱其為 Case 2a，由圖 4.42-4.45 可以發現當 Y 0 cm 5 . 3 <

λ

時，Case 2 和 Case 2a 的結果重合。圖 4.46-4.47 為 Case 3

中 A 點與 B 點的反力–負荷參數及位移–負荷參數曲線圖，圖 4.48 為 Case 3 中 C 點之位移–負荷參數曲線圖。由圖 4.46 可發現 B 點的

反力和位移的方向相反。由圖 4.41 可知當只有 A 點受位移負荷

λ

時，

B 點的位移大於0.5

λ

，所以當 B 點向位移負荷為0.5

λ

時，B 點的反力

(40)

第五章結論

本文以共旋轉(co-rotational formulation)有限元素推導法及增量疊代法來探討薄殼結構在位移負荷作用下的幾何非線性行為。由本文分析之數值例題的結果，可得以下的結論

(1)

本文所使用的數值計算方法可以分析結構受位移負荷作用時的行為，並有正確的結果。

(2)

在分析過程中亦發現本文在增量平衡迭代時將幾何剛度矩陣加入元素切線剛度矩陣中可以有效地改善收斂速度。

(3)

由本文例題結果可以知道結構受多個位移負荷作用與受多個力負荷作用時，結構的行為有很大的差異。

(4)

本文所使用的元素切線剛度矩陣僅為一近似剛度矩陣，故不能使用系統切線剛度的行列式值偵測平衡路徑的分歧點。

(41)

參考文獻

1. Croll JG, Walker AC. Elements of Structural Stability. New York: John Wiley and Sons; 1972.

2. Noor Ak, Hartley SJ. Nonlinear shell analysis via mixed Isoparametric elements. Comput Struct 1977; 7: 615-626.

3. Horrigmoe G, Bergan PG. Nonlinear analysis of free-form shells by flat finite element. Comput Meth Appl Mech Engng 1978;16:11-35. 4. Pica A, Wood RD, Hinton E. Finite element analysis of geometrically

nonlinear plate behavior using a Minflin formulation. Comput Struct 1980; 11: 203-215.

5. Hughes TJR, Lui WK. Nonlinear finite element analysis of shell: Part1. Three-dimensional shells. Comput Meth Appl Mech Engng 1981; 26: 331-362.

6. Surana KS. Geometrically nonlinear formulation for the axisymmetric Shell elements. Int J Numer Meth Engng. 1982; 18: 477-502.

7. Surana KS. Geometrically nonlinear formulation for the curved shell element. Int J Numer Meth Engng 1983; 19: 581-615.

8. Barut A, Madenci E, Tessler A, Stames JH. A new stiffened shell element for geometrically nonlinear analysis of composite laminates. Comput Struct 2000; 77: 11-40.

9. Levy R, Gal E. Geometrically Nonlinear three-noded flat triangular shell element. Comput Struct 2001; 79: 2349-2355.

10. Maccarini RR, Saetta A, Vitaliani R. A non-linear finite element formulation for shells of arbitrary geometry. Comput Meth Appl Mech Engng 2001; 190: 4967-86.

11. Kuznetsov VV, Levyakov SV. Nonlinear pure bending of toroidal shells of arbitrary cross-section. Int J Solids Struct 2001; 38: 7343-54.

薄殼結構在位移負荷作用下之幾何非線性分析

國 立 交 通 大 學

機械工程學系碩士班

碩 士 論 文

薄殼結構在位移負荷作用之下的幾何非

線性分析

Geometrically Nonlinear Analysis of Thin

Shell Structures under Displacement Type

Loading

研 究 生：楊禮龍

指導教授：蕭國模 博士