國
立
交
通
大
學
機械工程學系碩士班
碩 士 論 文
拘限彈性薄板受側向位移負荷的行為研究
A Study on the Behavior of Constrained Thin
Elastic Plate under Lateral Displacement
Loading
研 究 生：楊水勝
指導教授：蕭國模 博士
拘限彈性薄板受側向位移負荷的行為研究
A Study on the Behavior of Constrained Thin Elastic Plate
under Lateral Displacement loading
研 究 生：楊水勝 Student：Shui-Shen Yang
指導教授：蕭國模 博士 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
拘限彈性薄板受側向位移負荷的行為研究
研究生：楊水勝 指導教授：蕭國模博士國立交通大學機械工程學系碩士班
摘要
本文主要目的為探討拘限彈性薄板受側向力作用下的行為。本文以共 旋轉(co-rotational formulation)有限元素法及增量疊代法來探討薄殼的幾何 非線性行為。 本文採用基於牛頓-拉福森(Newton-Raphson)法和弧長控制(arc-length control)法的增量疊代法來解結構的非線性平衡方程式。本文以數值分析模 擬文獻上之拘限彈性薄板在受側向集中位移負荷後的實驗，本文求出完整 的平衡路徑，分布反力的變化，探討結構變形的轉換行為，並與文獻上的 實驗結果及分析結果比較。本文亦探討側向集中位移負荷偏移及結構邊界 偏移對結構變形行為的影響。A Study on the Behavior of Constrained Thin Elastic Plate under
Lateral Displacement loading
Student：Shui-Shen Yang Advisor：Dr. Kuo-Mo Hsiao
Department of Mechanical Engineering National Chiao Tung University Hsinchu, Taiwan
ABSTRACT
The objective of this paper is to investigate the behavior of constrained thin elastic plate under lateral displacement loading. The geometrically nonlinear behaviors of thin shell are investigated using the co-rotational finite element formulation and a incremental-iterative method.
An incremental-iterative method based on the Newton-Raphson method and constant arc-length method is used for solving nonlinear equilibrium equations. Numerical analysis is carried out to simulate an experiment of the constrained thin elastic plate under lateral displacement loading reported in the literature. The complete equilibrium path and the variation of distributed reaction are obtained. The transition of the deformation patterns is investigated and compared with the experimental and analytical result in the literature. The influence of the variance of lateral displacement loading point and that of the structural boundaries conditions on the transition of the deformation patterns is also investigated.
誌謝
衷心感謝指導教授 蕭國模博士在這兩年碩士班期間的指導與教誨，使本 論文得以順利的完成。老師在研究上嚴謹的態度及日常生活上的關懷，使我受 益良多，在此致上最高的謝意與敬意。感謝尹慶中老師及蔡佳霖老師撥冗擔任 口試委員並對本論文所提出的指正與建議，使本論文能夠更臻完善。 感謝陳弘虎、蔡明旭、劉峰成、劉宗帆、楊禮龍學長的照顧，同學黃孝衡 以及學弟顏宏儒、林育丞在學業上的砥礪與成長。 感謝父母親、兩位姊姊以及鄭念慈、侯玉真等關心我的朋友對我的支持與 鼓勵，僅以此成果與榮耀，獻給我親愛的父母以及所有關心我的人。目 錄 中文摘要 ... Ⅰ 英文摘要 ... Ⅱ 致謝 ... Ⅲ 目錄 ... Ⅳ 圖目錄 ... Ⅵ 第一章 緒論 ... 1 第二章 理論推導... 4 2.1 問題的描述 ... 4 2.2 座標系統 ... 4 2.3 平面三角殼元素... 5 2.3.1 殼的基本假設... 5 2.3.2 殼元素變形的描述... 5 2.3.2.1 常應變三角元素(CST)... 6 2.3.2.2 DKT 元素的變形描述 ... 8 2.4 元素內力與元素剛度矩陣... 12 2.4.1 CST 元素之節點內力與剛度矩陣 ... 13 2.4.2 DKT 元素之節點內力與剛度矩陣 ... 13 2.5 元素幾何剛度矩陣... 14 2.6 元素變形角的描述... 15 2.7 旋轉向量 ... 17 2.8 系統的平衡方程式與收斂準則... 17 第三章 數值計算方法與程序... 19 3.1 增量迭代法 ... 19 3.2 二分法 ... 21
3.3 數值程序 ... 22 第四章 數值分析與結果... 24 4.1 元素的收斂性... 24 4.2 完美結構的分析... 26 4.3 位移負荷λE偏移... 28 4.4 不完美結構的分析... 29 第五章 結論與未來研究方向... 31 5.1 結論 ... 31 5.2 未來研究方向... 31 參考文獻 ... 33 附圖 ... 38 附錄 A DKT 的形狀函數... 85 附錄 B CST 的剛度矩陣... 91
圖目錄 圖 1.1 文獻[33]實驗所觀察到四種變形轉換圖(a-d)及 對應於 a-c 圖結構的上視圖(e-g)... 38 圖 2.1 柱狀薄板結構示意圖... 39 圖 2.2 三角元素的示意圖及節點自由度... 40 圖 2.3 CST 元素在元素座標上的變形位移 ... 41 圖 2.4 DKT 元素的節點及其三邊上的局部座標示意圖 ... 42 圖 2.5 變形前板元素中心面之單位法向量 n 受旋轉向量 θ 作用的示意圖... 43
圖 2.6 元素座標的剛體旋轉 (a)面外旋轉(out-of plane rotation) (b)面內旋轉(in plane rotation) ... 44
圖 2.7 決定板元素節點變形轉角的第 3 個步驟之示意圖... 45 圖 2.8 旋轉向量... 46 圖 3.1 (a)第一階段位移負荷圖 (b)第二階段位移負荷圖 ... 47 圖 4.1 (a)圓柱薄板結構示意圖 (b)俯視圖 (c)前視圖 ... 48 圖 4.2 (a)結構示意圖 (b)位移負荷圖... 49 圖 4.3 結構中點的反力-位移負荷參數曲線圖... 50 圖 4.4 第一階段 E 點之位移-負荷參數圖 ... 51 圖 4.5 第二階段 E 點之反力-負荷參數圖 ... 51 圖 4.6 E 點之反力-負荷參數曲線圖 ... 52 圖 4.7 第一階段不同位移負荷λθ下 AB 線段在 Z 方向的 變形圖 ... 53 圖 4.8 第一階段不同位移負荷λθ下 CD 線段在 Z 方向 對 E 點的相對位移圖 ... 53
圖 4.9 第一階段邊界 HI 在不同位移負荷λθ下 X 方向的 反力分布圖... 54 圖 4.10 第一階段邊界 HI 在不同位移負荷λθ下 Y 方向的 反力分布圖... 55 圖 4.11 第一階段邊界 HI 在不同位移負荷λθ下 Z 方向的 反力分布圖... 56 圖 4.12 第二階段邊界 HI 在不同位移負荷λE下 X 方向的 反力分布圖... 57 圖 4.13 第二階段邊界 HI 在不同位移負荷λE下 Y 方向的 反力分布圖... 58 圖 4.14 第二階段邊界 HI 在不同位移負荷λE下 Z 方向的 反力分布圖... 59 圖 4.15 第二階段不同位移負荷下 CD 線段在 Z 方向的變形圖 .... 60 圖 4.16 第二階段不同位移負荷下 AB 線段在 Z 方向的變形圖 .... 61 圖 4.17 第二階段不同位移負荷λE下 CD 線段之位移分布圖... 62 圖 4.18 第二階段不同位移負荷λE下 AB 線段之位移分布圖... 63 圖 4.19 第二階段 E 點 X 向反力-負荷參數圖 ... 64 圖 4.20 第二階段 E 點 Y 向反力-負荷參數圖 ... 65 圖 4.21 (a)第二階段下λE =4.339mm 的結構變形圖... 66 (b)第二階段下λE =10.05mm 的結構變形圖... 66 (c)第二階段下λE =13.67mm 的結構變形圖 ... 67 (d)第二階段下λE =15.08mm 的結構變形圖... 67 圖 4.22 結構變形圖上視圖(a)λE =4.339mm (b)λE =10.05mm (c)λE =13.67mm (d)λE =15.08mm ... 68 圖 4.23 (a)位移負荷偏移不同距離後的位置圖 ... 69
(b)位移負荷向右偏移不同距離下中點的 Z 向反力- 位移負荷參數曲線圖... 69 圖 4.24 位移負荷向右偏移不同距離下中點的 X 向反力- 位移負荷參數曲線圖... 70 圖 4.25 位移負荷向右偏移不同距離下中點的 Y 向反力- 位移負荷參數曲線圖... 71 圖 4.26 位移負荷向右偏移 0.1mm 的結構變形上視圖 (a)λE=4.339mm (b)λE=10.05mm(c)λE=13.06mm ... 72 圖 4.27 E 點向右偏移 0.1mm 時在不同位移負荷λE ′下 AB 線段之位移分布圖 ... 73 圖 4.28 E 點向右偏移 0.1mm 時在不同位移負荷λE′下 CD 線段之位移分布圖 ... 74 圖 4.29 (a)位移負荷偏移後的位置圖(b)位移負荷偏移後的位置 為E4時在不同位移負荷下的結構變形上視圖... 75 圖 4.30 (a)位移負荷偏移後的位置圖(b)位移負荷偏移後的位置 為E5時在不同位移負荷下的結構變形上視圖... 76 圖 4.31 邊界不完美之結構示意圖... 77 圖 4.32 不完美結構之 E 點反力-位移負荷參數曲線圖 ... 78 圖 4.33 不完美結構在第二階段不同位移負荷λE下邊界H ′′I 在 X 方向的反力分布圖 ... 79 圖 4.34 不完美結構在第二階段不同位移負荷λE下邊界H ′′I 在 Y 方向的反力分布圖 ... 80 圖 4.35 不完美結構在第二階段不同位移負荷λE下邊界H ′′I 在 Z 方向的反力分布圖 ... 81 圖 4.36 (a)不完美結構在第二階段下λE =4.325mm 的結構變形圖 (b)不完美結構在第二階段下λE =12.31mm 的結構變形圖
... 82 圖 4.36 (c)不完美結構在第二階段下λE =13.68mm 的結構變形圖 (d)不完美結構在第二階段下λE =17.8mm 的結構變形圖 ... 83 圖 4.37 不完美結構變形圖的上視圖 (a)λE = 4.325mm (b)λE =12.31mm (c)λE =13.68mm (d)λE =17.8mm ... 84
第一章 緒論 在近代工程設計的發展上，對材料的要求與結構的表現趨向於高 強度與輕量化，不論機械、建築結構、航太設備及運輸工業等，設計 上必須考量以最小成本及重量來達到所要求的機能與強度，同時兼顧 外型的美觀。由於薄殼在承受彎曲應力與拉伸應力時有良好的表現， 又擁有極佳的重量強度比，因此薄殼結構常被廣泛應用在工程及生活 上。目前常見的薄殼結構有建築屋頂、飛機蒙皮、液體儲存槽、人造 衛星、火箭、船體結構等。 由於薄殼結構常用來承受大旋轉和大位移，即使在彈性範圍內， 其位移與受力的關係通常呈非線性變化，因此必須使用非線性分析來 探討由幾何形狀改變所造成的非線性行為。殼的研究從最早的線性分 析進入到較實際的材料非線性及大變形問題；後來殼的不穩定性、挫 屈及挫屈後行為和幾何非線性等問題則陸續被探討。殼的大位移分析 中，經常有snap-through和snap-back的現象，平衡路徑上亦常出現極 限點(limit point)或分歧點(bifurcation point)，為了求得完整的平衡路 徑，有很多方法被提出[1-6]。為測試殼元素的精確性及數值方法的有 效性，文獻中經常分析不同厚度及邊界條件的圓柱殼在側力下的行為 [7-15]，在對稱的情況下，為了節省計算時間，文獻上通常僅分析四 分之ㄧ的殼結構，文獻上發現殼的厚度減少時會出現snap-through、 snap-back及limit point，但沒有發現分歧點，文獻[2]分析二分之ㄧ的 殼結構發現在極限點前有一分歧點及次要平衡路徑。所以為了得到完 整的平衡路徑，應分析整個結構。 探討殼結構在軸力或側向力作用之挫屈行為的文獻很多[16-22]， 文獻[23]以實驗和數值方法探討一聚酯圓柱薄殼受位移負荷作用後的
非線性行為，文中將一薄片的兩個長邊以夾鉗固定，再將其彎成柱狀 殼，兩邊夾持端相距一固定距離，並與水平面維持一固定角度，然後 在薄殼中心施加一向下的位移負荷。在其實驗中隨著位移負荷的增 加，結構連續產生四個特殊的幾何變形，如圖1.1所示。第一個變形 是在薄殼中心附近出現兩個對稱X、Y軸的d-cone (developable cone) (圖1.1a)。第二個變形中出現兩個新的d-cone，而四個d-cone圍成一個 對薄殼中心轉了一個角度的菱形(圖1.1b)。第三個變形為四個d-cone 的連線形成一個梯型(圖1.1c)。第四個變形為梯形底邊兩個d-cone移到 薄殼自由端的邊界時，產生一個不連續的變化，使薄殼變成波浪狀的 圓柱殼(圖1.1d)。文獻[23]提到殼的長度若不夠長，則無法觀察到菱形 及後來的梯形圖形這兩種型態，會提早產生波浪狀的柱狀殼。文獻[23] 顯示負載-位移曲線的數值解與實驗結果相當吻合，文獻[23]並提到有 些階段能量及變形型態無法由數值結果看出，可能是因為離散不夠的 關係，文獻[23]是依據文獻[24]的數值方法來計算薄殼的能量分布， 文獻[24]是以Von Karman板殼理論及最小能量法來計算彈性板的變 形，並且在有限元素分析中，以一個五次多項式的函數來描述三個方 向u, v, w的位移。 據本文所知，在研究圓柱殼的非線性行為的文獻中，鮮少有人以 數值方法分析類似文獻[23]的結構及其實驗所觀察到的一連串現象， 而且文獻[1,2,10-12]分析的殼結構幾乎都是正方形，沒有文獻[23]所使 用的長方形。另外，文獻[23]實驗中殼結構的厚度與結構跨距的比值 43 . 471 1 ≈ d h 相對於文獻[1,2,10-12]的 80 1 ≈ d h 要來的小，文獻[23]殼結構 的曲率則比文獻[1,2,10-12]的大。雖然文獻[23]的負載-位移曲線圖顯 示其數值解與實驗的曲線相當接近，但是由其數值結果卻無法觀察到 實驗中出現的菱形旋轉及梯形的變形轉換。因此本文主要目的為使用
有限元素法並取整體結構來分析，找出柱狀薄殼受側向集中位移負荷 後的可能平衡路徑，觀察變形過程中結構的幾何形狀，探討此圓柱殼 的非線性行為和變形轉換可能的原因。
分析殼結構最常用的方法為有限元素法，文獻中有很多殼元素被 提出[7-10,25-39]，大致上可分為平面元素、曲面元素和等參數元素這 三類。文獻[30]中將CST(constant strain triangle)平面元素[26]與DKT (discrete Kirchhoff theory)三角板元素[29]疊加成一平面三角殼元素， 並使用更新拉格蘭日法(updated Lagrangian formulation)，將該元素用 在具大位移及大旋轉的薄殼結構分析中，由其例題的數值解可看出該 元素相當簡單及精確，故本研究擬採用文獻[30]的平面三角殼元素。 結構在受到力負荷及位移負荷的分析中，位移-負荷的平衡路徑 有時會出現snap-back與snap-through的情況[40]，文獻[41]提出基於牛 頓-拉福森(Newton-Raphson)法的弧長控制(arc-length control)法求解 力負荷作用下的非線性平衡方程式，並克服snap-back與snap- through 的情況，求得完整的平衡路徑。不過文獻[41]提出的方法，不能直接 使用在位移負荷的問題上，因此文獻[42]修改了文獻[41]的弧長控制 法，提出一個牛頓-拉福森(Newton-Raphson)法和弧長控制(arc-length control)法的增量疊代法來解決平面樑結構在位移負荷作用下的幾何 非線性問題，由於文獻[42]的方法非常有效，故本文將採用文獻[42] 的方法來解非線性平衡方程式。 本文將在第二章略為介紹本研究所使用的平面三角殼元素。在第 三章說明本研究所採用的數值計算方法以及分析時的數值程序。在第 四章探討柱狀薄殼的非線性行為，以及邊界不完美和位移負荷偏離柱 狀殼中心點時，對結構變形轉換的影響。
第二章 理論推導 2.1 問題的描述 考慮文獻[23]實驗所用的柱狀薄殼，如圖 2.1 所示。柱狀薄殼的 兩個短邊為自由端，兩長邊則以夾鉗固定，兩夾持端與 X-Y 平面(水 平面)為固定角度
α
= 20°，夾持端之間的跨距 d 固定為 16.5cm，薄板 的中心受到一集中位移負荷作用，觀察結構變形隨著位移負荷增加的 變化情形。 本文為了模擬此柱狀薄殼受集中負載後的行為，將採兩階段的位 移負荷來分析此問題。第一階段為形成柱狀薄殼受負荷前的初始狀 態，因此要先將矩形薄板彎成一個軸向與Y向平行的柱狀薄板，其執 行步驟為將薄板的兩個長邊以螺絲鉗夾持，螺絲鉗使薄板的兩個邊界 與水平方向呈一相同的角度；同時，右夾持端固定不動，左夾持端平 行向右推擠，使兩端的距離為期望之跨距。第二階段則是模擬柱狀薄 板受負荷後的狀態。其執行步驟為先重新設定兩夾持端的邊界條件， 將左夾持端的X向位移自由度 與兩夾持端的u Y向轉角自由度θ
y固 定後，再施加一個 向的集中位移負荷於柱狀板的中心。本文將以 有限元素法探討此柱狀薄板受側向集中負載下的變形，以及位移負載 偏移、結構的不完美對結構變形的影響。 Z − 2.2 座標系統 為了描述系統的運動及元素的變形，本文定義了兩組座標系統： (a)固定總體座標系統：XiG(i=1,2,3) 結構所有節點的座標、系統的邊界條件與其他座標系統的基底， 以及結構的平衡方程式，均在此座標系統中定義。(b)元素座標系統：xiE(i=1,2,3) 此座標系統建立在每一個殼元素變形後的最新位置上，其定義方 式將於 2.3 節中介紹。元素變形、元素內力及元素剛度矩陣在此座 標系統中定義，在經由標準的座標轉換，轉換至總體座標系統。 2.3 平面三角殼元素 本文使用文獻[10]的共旋轉法(co-rotational formulation)，本文與 文獻[10]一樣都採用文獻[30]中提出的非線性平面三角殼元素，如圖 2.2 所示，文獻[30]中的元素由 CST(constant strain triangle)三角形元素 [26]及文獻[29]中提出的 DKT(discrete Kirchhoff theory)三角殼元素所 疊加而成。為了本文完整性，在本節將簡單描述文獻[30]之殼元素變 形的假設、內力和剛度矩陣的推導，以及文獻[10]中元素變形角的決 定方法。 2.3.1 殼的基本假設 文獻[30]中對非線性平面三角殼元素的推導，做以下假設 (1) 薄膜變形及彎曲變形之間無偶合作用。 (2)殼元素的變形為小變形。 (3)在元素變形前，垂直於元素中心面的法向線段，在元素變形後，依 然保持直線，且沒有伸長或縮短，除了在元素三個頂點及三個邊的中 央點外，該線段不必垂直於變形後的中心面。 2.3.2 殼元素變形的描述 如圖 2.2 所示之殼元素中心面上有三個節點，其元素座標原點在 元素節點 1， E軸為元素節點 1 與元素節點 2 的連線， 軸是在元 x1 x2E
素平面上垂直於 軸，且朝元素節點 3 的方向， 軸則是由 軸與 軸外積而得。此殼元素每個節點有 6 個自由度，分別是 、 、 軸方向的位移 、 、 ( j=1,2,3)，以及繞 、 、 軸方向 的位移轉角 E x1 x3E x1E E x2 x1E x2E E x3 uj vj wj x1E x2E x3E xj
θ
、θ
yj、θ
zj( j=1,2,3)。此殼元素的變形可分為薄膜變形 (membrane deformation)與彎曲變形(bending deformation)兩部份，薄膜 變形是採用常應變三角形元素 (CST, constant strain triangle)[26]的薄 膜變形，彎曲變形是採用文獻[29]中的 DKT(discrete Kirchhoff theory) 三角殼元素的彎曲變形，由基本假設(1)可知此元素的變形可由薄膜變 形與彎曲變形疊加而成。 在圖 2.2 中的元素節點位移 與 ( j=1,2,3)是 CST 元素節點位 移，而 j u vj xjθ
、θ
yj及wj( j=1,2,3)為[29]中 DKT 元素節點位移，θ
zj是為了 不使元素剛度內的面內旋轉剛度(in-plane rotational stiffness)為零，而 人為加上的自由度。本文以下敘述的元素變形、元素內力以及元素剛 度矩陣都是在元素座標定義。 2.3.2.1 常應變三角元素(CST) 因為 CST 元素內的應變為常數，所以其位移場為線性位移場， 並可表示成： u =a1+a2x+a3y (2.1) v =a4+a5x+a6y (2.2) 其中 跟u v 為在 E軸與 軸方向的位移， x1 x2E x 與 是元素內任意點的座 標值， (i=1,2,…,6)是未定常數。 y i a 由元素座標的定義可知，在元素座標中圖 2.2 的節點 1、節點 2、 節點 3 之座標值可表示成(0,0)、(x2,0)、(x3,y3)，令節點 j在x1E軸、E x2 軸的位移分別是uj ,vj (j =1,2,3)。 本文中{}代表行矩陣。將三個節點的座標值及節點位移 代入(2.1)式、(2.2)式可以得到 j j v u , ) 3 , 2 , 1 (j= (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 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ 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.4) 由(2.3)、(2.4)式可將ai(i=1, 2 ....,, 6 表示成) 及 的函數，所以 (2.1)、(2.2)式可改寫成： i u vi u=Num (2.5) u={u,v} (2.6) N= ⎥ (2.7) ⎦ ⎤ ⎢ ⎣ ⎡ 3 2 1 3 2 1 0 0 0 0 0 0 N N N N N N um=
{
u1 v1 u2 v2 u3 v3}
(2.8) 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.9) 根據本文對元素座標的假設，元素座標原點定義在節點 1，元素 座標的 E軸定義在節點 1 與節點 2 的連線，如圖 2.3 所示，所以(2.8) x1式中um的u1=v1 =v2 = 0，故(2.8)式可表示成 =
{
(2.10) 其中 、 、 可由圖 2.3 中元素節點變形前後的座標決定出。 m u 0 0 u2 0 u3 v3}
2 u u3 v3 CST 元素的應變包含在x1E軸與x2E軸方向的應變εx與ε
y以及剪 應變γ
xy，因本文中假設元素的變形為小變形，所以εx、ε
y、γ
xy可 表示成： x u x ∂ ∂ =ε
, y v y ∂ ∂ =ε
, x v y u xy ∂ ∂ + ∂ ∂ =γ
(2.11) 將(2.5)式代入(2.11)式得到： εm =Bmum (2.12) εm={ε
x ,ε
y ,γ
xy} (2.13) Βm= ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ x y x y x y y y y x x x 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 ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − − + − − + − − = 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.14) 其中Bm稱為 CST 元素的位移-應變轉換矩陣。 2.3.2.2 DKT 元素的變形描述 圖 2.4 所示為文獻[29]中所提出的 DKT 板元素，節點 1、2、3 是 三角形的三個頂點，節點 4、5、6 為三角形三個邊的中點，這三個中的節點自由度。在圖 2.5 中，n為殼元素中心面變形前的單位法線向 量，n 為n在元素變形後的新位置，圖 2.5 中θ 為一在 平面上 的旋轉向量，將 作用在n可將n轉到 。由 2.3.1 的假設(3)可知垂 直於變形前的元素中心面法線向量變形後仍為直線且長度不變，所以 當 ′ E E x x1 − 2 θ nd θ <<1 時，DKT 元素的位移場可表示成： u =z
θ
y(x,y) v =−zθx(x,y) w=w(x,y) (2.15) 其中 x、 、y z為元素上任一點分別在x1E、x2E、x3E軸的座標值，θ
y是 在 軸方向的分量， θ E x1 θx是θ 在 軸方向的分量， 是在 軸方向 的位移， E x2 u x1E v 是在x2E軸方向的位移， 是在w x3E軸方向位移。當 θ <<1 時，θ
y與θx可視為法向量 繞n x1E軸及x2E軸的轉角。 DKT 元素的應變變形包含面內(in plane)正應變εx、ε
y與剪應變 xyγ
以及橫向剪應變(transverse shear strain)γ
yz、γxz。因本文假設元素的變形為小變形，εx、
ε
y和γ
xy可表示成(2.11) 式，γ
yz、γxz可表示成： z u x w xz ∂ ∂ + ∂ ∂ =γ
z v y w yz ∂ ∂ + ∂ ∂ =γ
(2.16) 將(2.15)式代入(2.11)式可得： κ εb = z (2.17) εb ={ε
x ,ε
y ,γ
xy} (2.18) } , , {θ
y,x −θ
x,yθ
y,y −θ
x,x = κ (2.19) 將(2.15)式代入(2.16)式可得： γ ={γ
xz ,γ
yz}={w,x +θ
y , w,y −θ
x} (2.20) 由(2.15)式可知 、wθ
y、θx與 無關，所以可由(2.20)式知橫向剪應 變在厚度方向為常數。 E x3本文中稱圖 2.4 中沿著元素邊緣方向 為切線方向，而垂直於元 素邊緣方向 為法線方向。 s n 在[29]中對於其所提出的 DKT 元素做了下列的假設： (1)
θ
y、θx在元素內為二次變化，也就是：∑
； (2.21) = = 6 1 i yi i y Nθ
θ
∑
= − = 6 1 i xi i x Nθ
θ
其中θyi、θxi是θ
y、θx在圖 2.4 中節點i的節點值， 為形狀函數[12]，其表示式詳見附錄 A。 ) 6 , , 1 (i = L Ni (2) 元素的三個頂點以及三個邊的中點滿足克希霍夫板理論 (Kirchhoff plate theory)的假設，即(a) 在三個頂點
γ
xzi =w,xi +θ
yi =0 i=1,2,3 (2.22a)γ
yzi =w,yi −θ
xi =0 i =1,2,3 (2.22b) 其中wxi、wyi、θ
yi、θxi分別是 ( ) x w wx ∂ ∂ = 、 ( ) y w wy ∂ ∂ = 、θ
y、 x θ 在節點i的值。 (b) 在三個邊的中點 0 , = + −θ
nk wsk k =4,5,6 (2.23a) 0 , = + nk sk wθ
k =4,5,6 (2.23b) 其中θnk、θsk分別是θn、θs在節點 的值，k θn與θs分別是θ 在 與 n s 方向的分量，w,sk、w,nk分別是 , ( ) s w ws ∂ ∂ = 、 ) n w ( ,n ∂ ∂ = w 在 節點 k 的值。 (3) w在元素邊緣的方向上是呈現三次變化，也就是： sj j si i sk w w w w w, =− 3 −1 , + 3 −1 , (2.24)其中wi、wj是 在節點w i及 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.25) 其中θsk、θsi、θ
sj分別是θs在節點 k 、i、 j之值，θs是θ 在 s 方 向的分量，在圖 2.4 中節點k =4,5,6分別為邊 23、邊 31、邊 12 的中點。 在圖 2.4 中元素三個邊上的θ
y、θx與θs、θn之幾何轉換關係可 表示成： ⎥ (2.26) ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ s n y x c s s cθ
θ
θ
θ
x w, 、w,y與w,s、w,n的幾何轉換關係為： ⎥ (2.27) ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ s n y x w w c s s c w w , , , , 其中c =cosα
ij、s=sinα
ij，α
ij為元素的邊ij上的法線 與 軸的夾 角，見圖 2.4。 ij n x1E 由(2.22)-(2.27)式可以把(2.21)式表示成[29]： b T x y H (ξ
,η
)uθ
= b T y x H (ξ
,η
)uθ
=− (2.28) = b u [w1θ
x1θ
y1 w2θ
x2θ
y2 w3θ
x3θ
y3] (2.29) 其中 為 DKT 元素的節點位移， 與 是對應於元素節點位移的 新形狀函數，其表示式詳見附錄 A， b u Hx Hyξ
與η
是元素內任一點在元素自 然座標[29]的座標值，其中1≤ξ
≤0、1≤η
≤0。將(2.28)式代入(2.19)式可以得到： κ =Bbub (2.30) 其中Bb為 DKT 元素的位移-應變轉換矩陣，表示式為： ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ + + − + − = T y T x T x T y T y T x b y x x x x y A ξ η ξ η ξ ξ
η
ξ
, 3 , 2 , 3 , 2 , 3 , 3 2 1 ) , ( H H H H H H B (2.31) 其中 2 3 2y x A= 為三角形面積。 2.4 元素內力與元素剛度矩陣 本文中殼元素的節點內力是由 CST 及 DKT 元素的節點內力組合 而成，元素剛度矩陣是由 CST 元素剛度矩陣 、DKT 元素剛度矩陣 以及面內旋轉剛度 所疊加而成， 為一 m k b k kθz kθz 3× 的對角矩陣，其3 對角線元素的值是取 之對角線元素的絕對值中的最小值。本節中 將用虛功原理推導 CST 元素及 DKT 元素的節點內力及剛度矩陣。 b k 在平面應力狀態，等向性線彈性材料的應變與應力關係為 ε E σ= (2.32) } , , {σ
xσ
yτ
xy = σ } , , {ε
xε
yγ
xy = ε (2.33) ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − = 2 1 0 0 0 1 0 1 1 2ν
ν
ν
ν
E E (2.34)τ
yz =Gγ
yz, τxz =Gγxz (2.35) 其中 E 是楊氏模數(Young’s module)，ν
是蒲松比(Poisson ratio)，G 為 剪力模數，(2.32)式之 可以是(2.13)式中的ε εm及(2.18)式中的 εb。2.4.1 CST 元素之節點內力與剛度矩陣 將(2.32)式代入(2.12)式可得 σm = EBmum (2.36) 由虛功原理可得 dV (2.37) V m t m m t mf =
∫
ε σ u δ δ } { m1 m2 m3 m f f f f = fmi ={fxj fyj} j=1,2,3 (2.38) 其中fm是 CST 元素對應於δum的節點內力，V 是元素的體積， 定 義於(2.10)式。 m u 將(2.12)式、(2.36)式代入(2.37)式可得 dV (2.39) V m m t m t m m t mf = u∫
B EB u u δ δ 由(2.39)式可得 fm =kmum (2.40) dV (2.41) V m t m m =∫
B EB k 其中km是 CST 元素的剛度矩陣，而Km的表示式詳見附錄 B。 2.4.2 DKT 元素的節點內力及剛度矩陣 將(2.17)式、(2.30)式代入(2.32)式可得 σb = zEBub (2.42) 假設在薄殼中剪應力τxz與τxz所做的虛功可以忽略，所以本文中用虛 功原理推導 DKT 元素的節點內力時僅考慮σx、σ
y及τ
xy所做的虛功。 由虛功原理可得 dV (2.43) V b t b b t bf =∫
ε σ u δ δ f ={f f f }fbi ={fzj mxj myj} j =1,2,3 (2.44) 其中fb是 DKT 元素對應於δub的節點內力，ub定義於(2.29)式，V 為 DKT 元素的體積。 將(2.17)式、(2.30)式、(2.42)式代入(2.43)式可得 =
∫∫
A b b t b t b b t bf u zB zEB u dzdA u δ δ =∫
(2.45) 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 2 2 2ν
ν
ν
ν
Eh dz z h h b E D (2.46) 其中 為元素厚度。 h 由(2.45)式可得[29] fb =kbub (2.47) η ξ η d d A tb b b b∫ ∫
− = 1 0 1 0 2 B D B k (2.48) 其中kb是 DKT 元素剛度矩陣。 2.5 元素幾何剛度矩陣 元素的幾何剛度矩陣對平衡迭代時的收斂速度有很大影響。本文 在平衡迭代的過程中將元素幾何剛度矩陣加入元素剛度中，以提高收 斂速度。本文採文獻[30]中的元素幾何剛度，其表示式為 dA A g t g g =∫
B NB k (2.49) 其中⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − − = 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 其中 x32 = x3 − x2 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ′ 2 12 12 1 N N N N N , N =hEBmum (2.49)式為一近似的幾何剛度矩陣，僅考慮面內應力σx、
σ
y及τ
xy 在剛體運動時的效應。 2.6 元素變形角的描述 本文中採用文獻[10]中的兩階段旋轉法來決定元素的剛體旋轉及 節點在元素座標的變形位移與轉角。假設第 I 個位置已知，此處第 I 個位置是指第 I 個增量迭代收斂後的平衡位置。 (i=1,2,3)為元素在 第 I 個位置的元素座標， 、 E i I x G i I X ∆ 及Ui ∆ (i=1,2,3)分別是元素節點 iΦi 在總體座標中第 I 個位置的位置向量、增量位移向量，及增量旋轉向 量，節點 i 當前的座標XiG(i=1,2,3)可由IXiG加上∆ 得到。由Ui 可 以利用元素座標定義求出當前的元素座標 (i=1,2,3)，本文稱元素從 到 的運動為對應於 G i X E i x E i I x xiE ∆Uj及∆Φj(j=1,2,3)之剛體運動，本文中假 設該剛體運動是由以下三個步驟達成的。(1) ∆U1造成的位移：元素座標 E的原點受到 i I x ∆ 的作用移動到元素U1 座標 的原點，其中 是作用在元素節點 1 的元素節點增量位移向 量。 E i x ∆U1
(2) 旋轉向量 α 造成的面外旋轉(out-of plane rotation)： 在圖 2.6(a)中可以看到元素的 (i=1,2,3)軸因為受到在 平面的旋轉向量 作用而旋轉，其中 軸旋轉到 軸處並與之重 疊， 的表示式為 E i I x x1E −x2E α I E x3 x3E α 3 3 3 3 3 3 1 ) ( cos e e e e e e α × × ⋅ = − I I I (2.50) 其中Ie3與 分別是 軸與 軸上的單位向量。 3 e I E x3 x3E (3) 旋轉向量 β 造成的面內旋轉(in-plane rotation)： 元素的 (i=1,2,3)軸因為受到在 平面的旋轉向量 作用 而旋轉到 (i=1,2,3)軸，由圖 2.6(b)可知 軸與 軸重疊，而 軸 受到旋轉向量 作用轉到 ，旋轉向量β 的表示式為 E i I x x1E −x2E α E i x′ x3′E x3′E x′iE β E i x β=cos−1(e1′⋅e1)e3 (2.51) 其中 、 與e′1 e1 e3分別是 E、 、 軸方向的單位向量。 x1′ x1E x3E 本文使用文獻[10]提出的直接法(direct method)來計算對應於 及 (j=1,2,3)的元素變形角。文獻[10]所提的直接法可分成以下 四個步驟。 i U ∆ ∆Φi (1) ∆U1造成的剛體移動：整個元素因為∆ 而移動，其中U1 是作 用在元素節點 1 的元素節點增量位移向量。 1 U ∆ (2) 旋轉向量 α 造成的剛體旋轉：旋轉向量 α (2.50 式)通過元素節點 1 使得整個元素除了Indj之外由IxiE(i=1,2,3)座標轉到x′iE(i=1,2,3)座
標，其中 是第 I 次迭代收斂後位在元素節點 j 上垂直於元素中心 面的法線向量。 dj In (3)
∆
Φ
tj造成Indj的有限旋轉：法線向量 受到 dj In tjΦ
∆
的作用而旋轉 到新位置n′dj，如圖 2.7 所示。其中∆
Φ
tj是增量旋轉位移向量 在 平面的投影向量。 j Φ ∆ E E x x1′ − 2′ (4) 旋轉向量 β 造成的剛體旋轉：旋轉向量 β (2.51 式)通過元素節點 1 使得整個元素包括 及 旋轉到最後的位置 及 ，其中 是 受到 作用後的新位置， 是第 I 次迭代收斂後，垂直於未變形的 元素中心面之法線向量， 是當前垂直於未變形的元素中心面之法 線向量。 u n′ n′dj nu ndj n′u Inu α Inu u n 令 及 為將 與 表示成在當前元素座標 (i=1,2,3)軸分 量的向量，則元素變形角可表示成 E u n E dj n nu ndj E i x E dj E u E dj E u E dj E u E j E j E j n n n n n n θ θ θ × × ⋅ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ = − ) ( cos 0 1 2 1 (2.52) 2.7 旋轉向量 本文中使用旋轉向量來表示一個有限旋轉，如圖 2.8 所示，一向量 受到一旋轉向量 Rφ
n的作用而轉到一個新位置R′，R與R′之間的關 係可表示成[43]R′=cosφR+(1-cosφ)(n⋅R)n+sinφ(n×R) (2.53)
2.8 系統的平衡方程式與收斂準則
Ψ =F(Q,
λ
Qp)=0 (2.54) 其中 為系統不平衡力向量， 為系統節點內力， 為系統位移向 量， Ψ F Qλ
為負荷參數， 為參考位移負荷向量。 可由(2.38)與(2.44) 式的元素節點力，以標準的座標轉換，轉換到對應總體座標後再組合 而成。 p Q F本文以不平衡力 的 weighted Euclidean norm 做為平衡迭代時的 誤差度量，且收歛準則表示為 Ψ tol e e e= ≤ F Ψ (2.55) 其中 為對應於位移負荷的系統節點反力向量。 為一給定的容許 誤差值。 e F etol
第三章 數值計算方法與程序 本文為了模擬 2.1 節所敘述之問題，將採兩階段的位移負荷來分 析。在第一階段的模擬中，同時施加三組位移負荷，這三組位移負荷 分別為施加 向的位移X λd於X =0的左夾持端，和施加 向的轉角Y λθ 於X=W的右夾持端，以及施加 向的轉角-Y λθ於X=0的左夾持端， 如圖 3.1(a)所示。為了使施加的位移與轉角這三組位移負荷同時到達 期望值，位移與轉角三組負荷參數間必須維持一固定比率。待位移與 轉角達到給定期望值，則形成圓柱薄板受側向集中負荷之前的初始狀 態。 第一階段進行平衡迭代時的數值方法大部分採用文獻[42]中所提 出基於牛頓-拉福森(Newton-Raphson)法和弧長控制(arc length control) 法的增量迭代法，而為了使第一階段的邊界位移與轉角精確地達到期 望值，將於第一階段的最後一步增量使用標準牛頓法來達成。第二階 段一開始先重新設定兩夾持端的邊界條件，再施加一個 向的集中 位移負荷 Z − E λ 於圓柱板的中心(E 點)，如圖 3.1(b)所示。第二階段則全 程採用文獻[42]中所提的增量迭代法。 本文中為了求得分歧點，以系統切線剛度矩陣之行列式值為零來 判斷。本文採用二分法決定增量位移向量的長度，以求得系統切線剛 度矩陣之行列式值為零的平衡位置。為了求得次要平衡路徑，本文在 平衡路徑的第一個挫曲負荷分歧點加入一個與第一挫曲模態向量方 向一致的擾動位移。 3.1 增量迭代法 若第 I 個增量的平衡位置為已知，令其位移向量為QI、負荷參
數為λI，則在此位置的系統切線剛度矩陣 可以求得，且第 I+1 個 增量的初始增量位移向量 T K Q ∆ ，可以利用尤拉預測值(Euler predictor) 求得[40] ∆Q=∆λrT rT =−(KTI )−1FP (3.1)
λ
∂ ∂ = F FP 其中∆λ
為初始增量負荷參數， 為參考負荷向量 的切線解， 為 (2.54)式中的系統節點內力。 T r FP Fλ
∆ 可利用下式求出[41] ∆λ
=±∆l( )
rTtrT 1/2 (3.2) 其中正負號決定方法如下：若第 I-1 與 I 個增量收歛時，其系統切線 剛度矩陣之行列式值同號，則∆λ
的正負號和第 I 個增量時相同；若 異號則符號相反。∆l表示第 I+1 個增量的增量弧長，其值可以如下決 定[41] ∆l=∆lI(
JD JI)
12 (3.3) 其中 為給定的期望迭代次數， 為第 I 個增量迭代至平衡所使用 的迭代次數， D J JI I l ∆ 為第 I 個增量的增量弧長。 本文第一個增量的增量弧長∆ 是由下式決定 l1 c r l max 0 max 1 I R R = ∆ (3.4) 上式中Rmax為給定的參考自由度之最大位移量， R0 為參考負荷向 量 作用下的系統線性解之 Euclidean norm， 為給定之最大增量 次數， P F Imax c r 為R0在參考自由度的分量之絕對值。 在平衡迭代時增量位移向量 Q∆ 及增量位移負荷參數∆λ
已知，由 ， = +1 Q Q +∆Q λ + =λ +∆λ，參考位移負荷向量Q 及 2.3.2 節與2.5 節的方法，則可求得系統中各殼元素在當前的元素座標、節點位 移、節點變形角。再利用 2.4 節的方法求得當前元素座標上的節點內 力及剛度矩陣，再由座標轉換轉到總體座標的節點內力及剛度矩陣。 當系統內力及外力求得後，不平衡力向量 可由(2.54)求得。若(2.55) 的收斂準則不能滿足，則利用定弧長控制法[41]求得一位移修正量 Ψ Q
δ
與負荷參數修正量δλ
，並加入前一次迭代的∆ 與Q ∆λ
中，而得到 一新的增量位移向量與增量負荷參數，再進行下一次迭代。δ
Q與δλ
可由下列兩式決定δ
Q=−KT−1(ΨI+1+δλ
FP) (3.5)δ
l2 =(∆Q+δ
Q)t(∆Q+δ
Q) (3.6) 以上迭代過程一直重複至不平衡力滿足(2.55)式的收斂準則。 3.2二分法 利用 3.1 節的增量迭代法可以求得結構之主要平衡路徑。在每個 增量的迭代收歛時，可以得到該增量在其平衡位置的負荷參數λ
及結 構剛度矩陣的行列式值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 1 R L I l l l = ∆ +∆ ∆ + ，重做第 I+1 個增量迭代，並求得新的λI+1及 ) ( I+1 D λ 。(3) 若D(λI+1)大於零，則令λL =λI+1，∆lL =∆lI+1。若D(λI+1)小於零， 則令λR =λI+1，∆lR =∆lI+1。 (4) 若下列二式挫曲誤差準則同時滿足，則λI+1即為系統挫曲負荷，否 則回到步驟(2)，重新展開下一次二分法迭代。 I eD D D < + ) 0 ( ) (
λ
1 λλ
λ
λ
e I L R − < +1 其中 及 為給定的容許誤差值，本文例題之計算給定 及 。 D e eλ eD =10−8 5 10− = λ e 3.3數值程序 本文針對此問題之模擬與分析，採用的數值程序說明如下： 1. 輸入初始的結構資料、材料常數，與分析所需的階段數目、每個 階段的位移負荷數，以及轉換至下一階段所需變換的自由度數目。 2. 設定該階段結構的邊界狀態，選擇該階段的參考自由度，並給定 期望此自由度於該階段應達到的位移，以及設定該階段收斂時的 容許誤差。 3. 設定該階段的位移負荷，若該階段有兩個以上之位移負荷同時作 用，則輸入位移負荷間的比例大小。 4. 使用迭代法求增量的收歛解 (a)利用已知的增量位移求當前元素變形並計算元素節點內力。 (b)計算不平衡力Ψ，檢查是否滿足收斂準則，若滿足，則進行 5.； 否則減少增量弧長並計算新的增量位移與增量負荷，回到(a)重 新此增量。5. 計算下一次增量所需的資料 (a)檢查該階段之參考自由度的位移是否已達的給定值。若已達到 給定值則進行6.；否則進行下一步工作。 (b)計算切線剛度矩陣 與參考負荷向量 ，並計算下一個增量 的增量位移、增量弧長與增量負荷參數。 T K FP (c)回到2.執行迭代工作。 6. 更換邊界狀態，並把邊界條件更動的節點位移儲存，供下一階段 的分析工作所使用。
第四章 數值分析與結果 在本章中，將以第二章之殼元素及第三章所提的數值程序，分析 圖 4.1 所示之圓柱薄殼受側向集中負載下的變形，以及討論位移負載 偏移、結構的不完美對結構變形的影響。該分析為文獻[23]之實驗的 數值模擬，文獻[23]的實驗中將一長寬比為 2:1 的薄板用夾鉗夾住薄 板的兩個長邊，夾鉗將薄板彎成一圓柱殼，兩夾持邊與水平面夾 20 度角(約 0.3491rad)，然後施加一位移負荷於結構中心，觀察結構變形 的過程。本文分析時採用兩階段的位移負荷來進行數值模擬，在第一 階段中使用三組位移負荷分別為λd、λθ、-λθ，如圖 3.1(a)所示。經 由第一階段的三組位移負荷將一薄板彎成如圖 4.1 的結構。第二階段 則使用一位移負荷λE(如圖 3.1(b))施加於結構中心，模擬如圖 4.1 的 圓柱薄殼中心受到一側向位移負荷後的行為並與文獻[23]的實驗結果 比較。在本章所考慮之結構的幾何及材料性質皆採用文獻[23]中的資 料，其數值分別為長度 L=35cm、寬度 W=17.5 cm、厚度 h=0.35mm， 蒲松比
ν
=0.4，楊氏係數 9 2 Nm 10 3.8 E= × − ，並稱其為標準尺寸。本 章除特別聲明外，皆使用標準尺寸並取整體結構來進行分析。 4.1 元素的收斂性與準確性分析 為了解本文中所採用之數值程序的可行性及三角殼元素的準確 性，本節首先分析如圖 4.2(a)所示之圓柱殼，其結構中心受到一位移 負荷λE作用，如圖 4.2(b)。結構的兩個直邊為鉸接，兩個曲邊為自由 邊。本例題分別取(a)整個結構與(b)二分之一的結構進行分析。以整 個結構分析時，將其離散為 200 個三角殼元素(10×10 之網格)，使用 18 個增量，每個增量的平均迭代次數約為 4，平衡迭代的容許誤差值取 4，挫曲負荷參數 10 1× − λcr=9.481mm；以二分之一的結構分析時， 則將結構離散成 100 個三角殼元素(5×10 之網格)，使用 20 個增量， 每個增量的平均迭代次數約為 4，平衡迭代的容許誤差值取 ， 挫曲負荷參數 4 10 1× − cr λ =9.48mm。圖 4.3 是反力-位移負荷參數圖，圖 4.3 中 的 Present (a)是以整個結構分析的結果，Present (b)是以二分之一結構 分析的結果。圖 4.3 中RE為在位移負荷方向之反力，RE及λE相當於 文獻[2]中的力負荷及力負荷方向的位移，由圖 4.3 可見本文的結果與 文獻[2]的結果相當吻合，故本文的數值程序可以準確找出主要平衡路 徑、分歧點及次要路徑。 為了解本文所採用之三角殼元素的收斂性，本節用不同的元素網 格分析如圖 4.1(a)之標準尺寸的完美圓柱薄殼，其中心受到一集中位 移負荷作用。本例題分析時，將一平板分割成 24×48、30×60、40×80 的網格，並以兩階段位移負荷的步驟來分析。圖 3.1(a)(b)為本章例題 中薄殼受到位移負荷的種類，圖 3.1 (a)為第一階段的三組位移負荷參 數λd、λθ與-λθ，其中λθ=λ1(rad)、λd=28.6479λ1(mm)；圖 3.1 (b)為 第二階段的位移負荷參數λE，λE=λ2(mm)。本文中以標準尺寸的薄 板用兩階段位移負荷來分析，並將其視為與文獻[23]中的實驗等效。 文獻[23]實驗中所施加的位移負荷固定在圓柱薄殼正中心的位置，其 容許偏移誤差為 0.1mm，因此本文將第二階段的位移負荷設定在薄殼 正中心的節點(如圖 4.1(b)之 E 點)上，並固定 E 點在 X-Y 平面上位置， 使作用在 E 點的位移負荷不隨著結構變形而偏移。圖 4.4 及圖 4.5 則 分別為用 24×48、30×60、40×80 的網格分析下，所得到第一階段結 構中心(E 點)在 Z 向的位移與負荷參數λθ的關係曲線圖，以及第二階 段結構中心(E 點)在 Z 向的反力與負荷參數λE的關係曲線圖。由圖 4.4 及圖 4.5 的結果可知 24×48 的網格結果已足夠準確。為了讓結構的受
力與變形能有更細部的呈現，本文中以後的例題將一律使用 30×60 的網格來分析。圖 4.6 為第二階段負荷所得到的 E 點之反力-負荷參 數曲線圖(30×60 之網格)與文獻[23]的結果比較，可見本文在極限點 (λE=13.6mm)前的曲線與文獻[23]曲線的趨勢相當接近。 4.2 完美結構的分析 本節分析標準尺寸的殼結構，其尺寸與邊界條件皆為對稱的狀 態，所以為一完美結構。本例題分析過程中第一階段使用 19 個增量， 每個增量的平均迭代次數約為 3，平衡迭代的容許誤差值取 ； 第二階段使用 45 個增量，每個增量的平均迭代次數約為 4，平衡迭 代的容許誤差值取 。 5 10 1× − 4 10 1× − 圖 4.7-圖 4.22 為本例題之結果，圖 4.7 為第一階段下薄板受到三 組位移負荷作用下，結構上 AB 線段(如圖 4.1(b))隨著位移負荷λθ的 變形過程，圖 4.8 為第一階段下結構 CD 線段(如圖 4.1(b))在不同移負 荷λθ下對 E 點的相對位移；圖 4.9-圖 4.11 為第一階段下邊界反力分 布圖，由於第一階段結束後的結構邊界與水平面夾 20 度角，因此定 義邊界座標Y軸方向與總體座標之 Y 軸重合，而右邊邊界(如圖 4.1(b) 之 HI 邊)座標的X、Z軸為總體座標的 X、Z 軸繞 Y 軸順時鐘方向旋 轉 20 度後的位置，左邊邊界(如圖 4.1(b)之 FG 邊)座標的X、Z軸為 總體座標的 X、Z 軸繞 Y 軸逆時鐘方向旋轉 20 度後的位置，如圖 4.1(c)。本文以下之邊界反力皆以結構 HI 邊(如圖 4.1(b))的邊界反力 表示。由圖 4.8 中可發現 CD 線段的兩端產生較大的相對位移，CD 線段並非完全的水平，因此邊界兩端的反力比邊界其他地方明顯，由 圖 4.9 可看到邊界的兩端附近受到拉力外，其餘皆受到壓力。圖 4.12-圖 4.14 為第二階段不同位移負荷λ 下，在邊界所受到的反力圖，由
圖 4.12-圖 4.14 的結果可見完美結構在邊界X、Y、Z三個方向的反 力對邊界中點對稱。邊界反力FX在靠近結構四個端點(如圖 4.1(b)之 F、G、H、I)附近為拉力區且有最大值，邊界其他區域則為壓力區； 若將邊界反力除以結構邊界的斷面積，則邊界端點的X向最大拉應力 可達 4.025 MPa，X向最大壓應力約 0.675MPa，Y向最大應力可達 0.845 MPa。邊界反力FX在邊界中點附近(如圖4.1(b)之A、B點附近) 在λ =4.339mmE 時本來為拉力區，在λ =10.05mmE 之後則轉變為壓力 區。由圖4.15的CD線段變形圖中可觀察到施加於中點的位移負荷λE 造成結構兩個短邊(如圖4.1 (b)之FH、HI邊)翹起，為了抵抗在結構 短邊的翹曲，夾持端在端點(如圖4.1(b)之F、G、H、I)附近產生拉力 來阻止結構從夾鉗中抽離。圖4.16 則為第二階段AB線段在不同位移 負荷下 Z 方向的變形圖。結構之 AB 與 CD 線段(如圖 4.1(b))在 Z 方 向的變形皆對稱結構中點。 圖4.17為第二階段不同位移負荷λE下CD線段(如圖4.1(b))之X 向位移 變化圖，圖u 4.18則為第二階段不同位移負荷λE下AB線段(如 圖4.1(b))之Y向位移 變化圖。由圖v 4.17可發現雖然在λ =12.47mmE 時，接近 C、D 點(如圖 4.1(b))附近的 X 向位移 較其他位移負荷下 相對較大，但是基本上 X向位移 幾乎接近於零，而且位移分布對稱 於E點。圖4.18的結果亦呈現出AB線段之Y向位移 對稱於E點(如 圖4.1(b))，在位移負荷 u u v E λ =13.67mm時有最大值，AB線段之 Y向位 移 大約介於v −4×10−5到3×10−5mm之間；λ =15.08mmE 時位移 分布 則又趨於平緩，此時結構為一波浪狀的薄殼。當初在第二階段分析 中，為了使位移負荷精準的施加於結構中點(E點)，而固定E點在X、 Y上的位置，因此由圖 4.17與圖4.18中可見E點的位移 、 為零。 圖4.19、圖4.20分別為 E點在X與Y 向反力-負荷參數 v E u vE λ 的曲線圖，
從圖中可發現在λ =12.47mmE 到極限點(λ =13.64mm)之間有較大的E 反力產生，由圖4.17、圖4.18可知在λ =12.47mmE 時AB與 CD線段 產生相對較大的位移，而E點的 X、Y 向自由度又被固定，因此在E 點在位移負荷λ =12.47mmE 時產生反力。 圖 4.21(a)-(d)則為第二階段四種不同位移負荷參數下結構整體的 變形圖，圖 4.22 則為不同位移負荷下結構變形的上視圖。由圖 4.22 的(a)-(c)可觀察到 D-cones 在 AB 線段上產生，並且凹陷部分的脊線 圍成一個菱形，菱形區域則隨著位移負荷增加而變大，D-cones 往兩 邊界的中點(如圖 4.1(b)之 C 點、D 點)靠近，直到位移負荷超過極限 點後，結構變成一波浪狀的薄殼，如圖4.22(d)。此結果雖然可觀察到 D-cones 的產生以及結構轉換成波浪狀的變形，但是與文獻[23]的數 值結果一樣，無法觀察到實驗中產生的菱形區域對 E 點旋轉一個角 度，以及菱形區域變成梯型區域這兩種變形現象。 4.3 位移負荷λE偏移 本例題將使用標準尺寸之薄板來測試位移負荷偏離 E 點(如圖 4.1(b))對結構變形的影響，觀察是否產生 4.2 節中提到的兩個未出現 的變形轉換。由於文獻[23]提到其位移負荷的容許偏移誤差在0.1mm 以內，因此本例題測試採用三組偏移距離 分別為偏離 E 點右邊 0.01mm、0.05mm 及 0.1mm，使第二階段的位移負荷 R d E λ 作用的位置 產生偏移，其中 E 點的 X 座標原為 92.5mm，偏移後的 X 座標則為 92.51、92.55以及92.6mm，並假設偏移後的位置分別為 、 、 點(如圖4.23(a))，偏移後的位移負荷稱為 1 E E2 E3 E ′ λ 。 圖4.23(b)為 、 及 點反力-負荷參數曲線圖，結果顯示三種 位移負荷偏移的情況與位移負荷不偏移時( 1 E E2 E3 0 = R d )的情況幾乎相
同。圖 4.24 為三種偏移距離下位移負荷施加點的 X 向反力-位移負 荷參數圖，圖 4.25則為三種偏移距離下位移負荷施加點的Y向反力 -位移負荷參數圖，由圖 4.24 及圖 4.25 可知 X 與 Y 向的反力皆在 E ′ λ =12.47mm附近產生較大的變化。而偏右距離dR ≤0.1mm的整體 結構變形圖在變形為波浪狀之前，亦無觀察到文獻[23]實驗中產生的 菱形區域對 E 點旋轉一個角度，以及菱形區域變成梯型區域這兩種 現象。本文僅顯示 mm 時的結構變形上視圖(圖 4.26)，以及 mm 結構中線 AB、CD 線段(如圖 4.1(b))的位移分布圖(圖 4.27、4.28)。由圖4.27、4.28可知因為位移負荷偏移的關係，CD中 線的 X 向位移較位移負荷不偏移的時候還要來的大，而AB 中線 Y 向位移分布在 1 . 0 = R d 1 . 0 = R d E ′ λ =12.37mm 之後產生不對稱的分布，但對整體結構 變形的轉換仍無影響。 本文亦測試兩種位移負荷偏移的情況：(1)在 X-Y 平面上偏 E 點 上方 =0.1mm，偏移後的位置稱為 (如圖4.29(a))；(2)偏E點右上 方 =0.1mm，偏移後的位置稱為 (如圖4.30(a))，本文中僅顯示兩 種位移負荷偏移情況的結構在不同位移負荷下的變形上視圖，如圖 4.29(a)、圖4.30(b)，由圖4.30(b)可見結構在位移負荷 U d E4 RU d E5 E λ =13.18mm的 時候，結構變形凹陷的部份近似一個三角形，但還是與文獻[23]實驗 觀察的結果不同。因此第二階段施加於結構中心的位移負荷λE，其 位置在偏移誤差 0.1mm 的範圍內皆不會產生文獻[23]實驗中觀察到 的菱形區域對中點旋轉一個角度，以及菱形區域變成梯型區域這兩個 變形轉換，但是對於D-cones移動及結構變形還是有影響。 4.4 不完美結構的分析 本節將採用圖 4.31 所示的薄板以兩階段位移負荷來分析，兩個
階段所施加的位移負荷與4.1節中相同。本例題之薄板在第一階段兩 長邊 FG與 HI(如圖 4.31)用夾鉗夾持時，夾鉗分別以 A 點與B 點(如 圖4.31)為基準，在X-Y平面上旋轉一個角度，此時兩邊界與Y軸夾 一個小角度θc，且兩夾持端互相平行，形成新的邊界為 和 (如 圖 4.31)。此時薄板可視為一近似矩形的平行四邊形薄板。本例題考 慮 ， 、 、 G F′ ′ H ′′I ° = 32740. c
θ
F ′F G ′G H ′H 及I ′I 線段(如圖4.28)皆為1mm時的 情況，討論邊界不完美對結構變形的影響。 圖4.32為本例題之結構與完美結構之E點(如圖 4.1(b)、圖4.31) 反力-位移負荷參數的比較圖，由圖中可發現不完美結構在位移負荷 E λ 約在12.47mm與 13.61mm之間的E點反力較完美結構小，但位移 負荷到了λ =13.61mmE 之後不完美結構 E 點的反力則繼續增加。圖 4.33-圖 4.35 為不完美結構在第二階段分析的邊界反力分布圖，由圖 4.33-圖 4.35 可發現由於邊界歪斜一個小角度的關係，造成邊界反力 分布不對稱，其中邊界反力在λ =11.2mmE 之後G′端的反力比 端 大、 端的反力比 端大(參考圖 4.31)。此不完美結構在四個不同位 移負荷下的整體變形圖顯示在圖 4.36(a)-(d)中，在圖 4.36(c)-(d)中可 看到明顯的四個D-cones正往結構的兩個短邊發展，在圖 4.36(b)中可 看到D-cones所圍成的菱形區域對 Z軸轉了一個角度，可以由圖4.37 F′ H′ I′ (a)-(d)看到不同位移負荷下結構變形的上視圖中，獲得較清楚的圖形 。由圖 4.37(b)中雖可看到菱形區域的旋轉，但還是無法觀察到轉換成 梯型區域的圖形，只能由圖4.37(d)看到一平行四邊形區域的轉換，其 結果與文獻[23]實驗亦不相同，但可以知道結構邊界的不完美可能會 造成結構產生不同的變形。第五章 結論與未來研究方向 5.1 結論 本文為了模擬一矩形薄板在長邊部分以夾鉗挾持，夾鉗先將薄板 彎成一圓柱狀，再施加一集中位移負荷於結構中心的情況，採用兩階 段的位移負荷分析。本文在第四章模擬完美結構的非線性行為，得到 的分析結果並未觀察到文獻實驗中某些位移負荷下的變形轉換，因此 本文測試不完美結構以及第二階段位移負荷偏移的情形，並討論文獻 實驗中結構變形轉換的原因。由本文的數值模擬結果，可以得到以下 結論 1. 結構上 D-cones 的移動與變形轉換可能是因為結構受力不均的緣 故，因此完美結構的受力與變形皆對稱於結構中心，所以無法觀 察到菱形區域的旋轉與梯形區域的轉換這兩種不對稱的變形，若 將結構離散為更多的元素亦不會產生不對稱變形。 2. 本文經由測試結果發現 D-cones 的移動會因為位移負荷施加的位 置而有改變，若位移負荷施加的位置造成結構受力不對稱，對結 構變形的影響比較明顯。 3. 本文由模擬結果推測 D-cones 的移動與結構的變形，可能與結構受 到夾持的邊界狀態有關，若邊界不完美，則可能會造成結構變形 的不對稱。並且由完美結構的模擬結果得知邊界X向最大拉應力 可達 4.025 MPa 左右，X向最大壓應力約 0.675MPa，Y向最大應 力可達 0.845 MPa 左右，若實驗中夾具與結構接觸面間無法提供足 夠的摩擦力，則可能造成結構在邊界無法完全固定而有些微滑動。 5.2 未來研究方向 本文推測 D-cones 的移動和結構變形可能與結構邊界狀態有關，
![圖 1.1 文獻[23]實驗所觀察到四種變形轉換(a-d)及 對應於 a-c 圖結構的上視圖(e-g)](https://thumb-ap.123doks.com/thumbv2/9libinfo/8486469.184467/49.892.242.662.239.901/圖11文獻23實驗所觀察到四種變形轉換ad及對應於ac圖結構的上視.webp)
