國 立 交 通 大 學
機械工程學系碩士班
碩士論文
雙對稱開口薄壁 Timoshenko 梁之非線性分析
Nonlinear analysis of bisymmetric thin-walled open-section Timoshenko beam
研 究 生：莊士緯
指導教授：蕭國模 博士
雙對稱開口薄壁 Timoshenko 梁之非線性分析
Nonlinear analysis of bisymmetric thin-walled open-section Timoshenko beam
研 究 生： 莊士緯 Student： Shih-Wei Chuang
指導教授： 蕭國模 博士 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
中華民國一百零二年九月
雙對稱開口薄壁 Timoshenko 梁之非線性分析
Nonlinear analysis of bisymmetric thin-walled open-sectionTimoshenko beam 研究生：莊士緯 指導教授：蕭國模博士
國立交通大學機械工程學系碩士班
摘要 本 研 究 的 主 要 目 的 是 以 一 致 性 共 旋 轉 法 探 討 雙 對 稱 開 口 薄 壁 Timoshenko 梁的幾何非線性挫屈及挫屈後的分析。 本文中推導的梁元素有兩個節點，每個節點有 7 個自由度。本文中將元素 節點定在斷面剪心，並取剪心軸當作梁元素變形的參考軸。本研究在當前 梁元素變形的位置上建立元素座標，並在其上描述元素的變形。為了準確 描述 Timoshenko 元素的變形，本研究利用虛功原理並考慮剪力修正係數推 導元素節點內力，元素節點內力所作的虛功是在元素受虛位移擾動前的元 素座標上推導，但元素應力所作的虛功是在元素受虛位移擾動後的元素座 標上推導，即將元素座標建立在元素受虛位移擾動後的位置，並在其上定義元素的變形及推導虛應變。本研究推導的元素節點內力能滿足靜力的平 衡。本研究由元素節點內力的改變與擾動位移的關係推導梁元素的切線剛 度矩陣。因本研究在推導元素的節點內力時，扣除了虛位移中剛體運動的 部分，而元素的節點內力與元素一起剛體運動，所以不能僅由元素節點內 力對節點參數微分求得，還要考慮元素節點內力在剛體運動時因方向改變 造成的元素節點內力的改變。線性剛度矩陣包含在元素切線剛度矩陣裡面。 本文求解非線性平衡方程式的數值計算方法是基於牛頓－拉福森 (Newton-Raphson)法配合弧長控制(arc length control)法的增量迭代法。本研 究中以系統切線剛度矩陣之行列式值為零當作挫屈準則，利用弧長的二分 法求得挫屈負荷。
為了驗證本研究提出的方法的準確性與有效性，本研究以不同數值例 題探討剪應變對雙對稱開口薄壁梁之負荷—位移曲線及挫屈負荷的影響， 並與 Euler–Bernoulli 梁的結果做比較。
Nonlinear analysis of bisymmetric thin-walled open-section
Timoshenko beam
Student：Shih-Wei Chuang Advisor：Kuo-Mo Hsiao
Department of Mechanical Engineering
National Chiao Tung University, Hsinchu, Taiwan, R.O.C.
ABSTRACT
A consistent co-rotational total Lagrangian finite element formulation for the geometric nonlinear buckling and postbuckling analysis of bisymmetric thin-walled Timoshenko beams is presented.
The element developed here has two nodes with seven degrees of freedom per node. The element nodes are chosen to be located at the centroid of the end cross-sections of the beam element and the axis of centroid is chosen to be the reference axis. The deformations of the beam element are described in the current element coordinate system constructed at the current configuration of the beam element. The exact kinematics of the Timoshenko beam is considered. The element nodal forces are derived using the virtual work principle with the consideration of the shear correction factor. The virtual rigid body motion corresponding to the virtual nodal displacements is excluded in the derivation of the element nodal forces. A procedure is proposed to determine the virtual rigid body motion. A consistent second-order linearization of the element nodal forces is used here. Thus, all coupling among bending, shearing, twisting, and stretching deformations of the beam element is retained. In the derivation of the element tangent stiffness matrix, the change of element nodal forces induced by the element rigid body rotations should be considered for the present method. Thus, a stability matrix is included in the element tangent stiffness matrix.
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 matrix determinant of the structure is used as the criterion of the buckling state. A bisection method of the arc length is used to find the buckling load.
Numerical examples are studied and compared with the results obtained by using Euler beam element to demonstrate the accuracy and efficiency of the proposed method and to investigate the effect of the shear deformation on the loading–deflection curves and buckling load of the bisymmetric thin-walled beams.
誌謝
衷心感謝指導教授 蕭國模博士在這兩年期間的指導與教誨，使本論文得 以順利完成，蕭老師在研究上嚴謹的態度以及對日常生活上的關懷，使我受益 良多，在此致上最高的敬意及謝意。也感謝鄭文雅老師及尹慶中老師撥冗擔任 口試委員並對本論文所提出的指正與建議，使本論文能夠更臻完善。 感謝蔡明旭學長、盧致群學長、翁林甫學長、黃楚璋學長與許彤羽學 姐在研究上的協助與照顧以及生活上的互相照應。感謝同學林琮琪、沈佳 鴻在課業及研究上的幫忙。最後再感謝學弟林群禮、高嘉鴻、金長虹、蔡 耀賢在學業以及各方面的砥礪與成長。 感謝父母親、關心我的親人、及朋友對我的支持與鼓勵，僅以此成果與榮 耀，獻給我親愛的父母、家人以及所有關心我的人。目錄
中文摘要……….………...Ⅰ 英文摘要……….………...………... VI 致謝……….………...V 目錄.……..……….VI 表目錄..……….……….. VIII 圖目錄..………..………...IX 第一章 緒論 ………..1 第二章 理論推導 ……….….3 2.1 基本假設 ……….3 2.2 座標系統描述………...3 2.3 旋轉向量 ……….….4 2.4 Timoshenko梁的變形描述……….5 2.5 Tmoshenko梁的應變………....12 2.5.1梁的應變………...12 2.5.2梁之位置向量的變分及應變的變分……….14 2.6 節點參數與節點力……….……….16 2.6.1 系統節點參數與節點力……….16 2.6.2 元素節點參數與節點力……….172.6.3 元素顯節點參數的擾動與隱節點參數的擾動關係…….18 2.7 元素節點內力f之推導………....22 27.1 元素顯節點內力與隱節點內力的關係……….23 27.2 隱節點內力f 的推導………24 2.8 元素剛度矩陣之推導………....27 2.9 系統平衡方程式與收斂準則 ……….….31 第三章 數值計算方法與程序 ………...33 3.1 增量迭代法 ………...33 3.2 二分法 ………...35 3.3 N循環迭代法 ……….….36 第四章 數值例題與結果 …....……….………....37 第五章 結論與展望 ….……...……….………....47 參考文獻 ………..….49 附錄A ……….………...….52 附錄B ……….……….………...54 附錄C ……….………....57 附表..………...61 附圖..………...….…70
表 目 錄 表 4.1 例題一：矩形懸臂梁末端受集中負荷之挫屈負荷………...61 表 4.2 例題二：兩端固接矩形梁中間承受鉛直力之挫屈負荷………...62 表 4.3 例題三：I 形懸臂梁末端受一集中載重之挫屈負荷……….…...64 表 4.4 例題四：兩端固接 I 形梁中間承受鉛直力之挫屈負荷…….…...65 表 4.5 例題五：兩端固接 I 形梁承受兩個鉛直力之挫屈負荷..…...66 表 4.6 例題九之十字斷面性質………..………..67 表 4.7 例題十：懸臂梁末端斷面不同位置受集中載重之挫屈負荷…....68 表 4.8 例題十、十一之 I 形斷面性質………...69
圖 目 錄 圖 2.1 梁元素之位移以及座標系統關係圖………70 圖2.2 旋轉向量 ………71 圖 4.1 矩形懸臂梁末端承受集中載重之結構圖(例題一)….…………72 圖 4.2 端點 B 在X2G,X3G方向之負荷-位移曲線圖(例題一)...73 圖 4.3 矩形梁中間承受集中載重之結構圖(例題二)……….74 圖 4.4梁端點C在XG XG方向的負荷-位移曲線圖(例題二)…………75 3 2 , 圖 4.5 I 形懸臂梁末端承受集中載重之結構圖(例題三)…………...…76 圖4.6自由端B的負荷-位移曲線圖(例題三)………..….……77 圖4.7 I形梁中間承受一集中載重之結構圖(例題四)……….…78 圖4.8梁端點C在X2G,X3G方向的負荷-位移曲線圖(例題四)…………79 圖4.9梁端點C在X2G,X3G方向的負荷-位移曲線圖(例題四)…………80 圖4.10梁端點C在X2G,X3G方向的負荷-位移曲線圖(例題四)…...…81 圖4.11梁端點C在X2G,X3G方向的負荷-位移曲線圖(例題四)……..…82 圖4.12 I形梁承受兩個集中載重之結構圖(例題五)………...……83 圖4.13梁端點C在X2G,X3G方向的負荷-位移曲線圖(例題五)….…..…84 圖4.14梁端點C在X2G,X3G方向的負荷-位移曲線圖(例題五)……...…85 圖4.15梁端點C在X2G,X3G方向的負荷-位移曲線圖(例題五)……...86 圖4.16梁端點C在X2G,X3G方向的負荷-位移曲線圖(例題五)….…..…87 圖4.17簡支梁兩端承受偏心軸力之結構圖(例題六)………..……88
圖4.18簡支梁受偏心軸力的負荷-位移曲線圖(例題六)...…89 圖4.19懸臂L形構架自由端承受水平力之結構圖(例題七)……….…90 圖4.20端點 C在X1G方向的負荷-位移曲線圖(例題七)……….…91 圖4.21 端點 C在X2G方向的負荷-位移曲線圖(例題七)…...92 圖4.22 端點 C在X3G方向的負荷-位移曲線圖(例題七)……….…93 圖4.23簡支梁兩端承受均勻彎矩之結構圖(例題八)...………..………94 圖4.24簡支梁兩端皆為自由翹曲受均勻彎矩的負荷-位移曲線圖 (例題八)………..…………....95 圖4.25簡支梁兩端皆為抑制翹曲受均勻彎矩的負荷-位移曲線圖 (例題八)………...96 圖4.26 懸臂梁末端承受一軸力之結構圖(例題九)…….………..…...97 圖4.27 懸臂梁受軸力的負荷-位移曲線圖(例題九)………...98 圖 4.28 懸臂梁承受均佈載重之結構圖(例題十)……….………....99 圖4.29偏心力施加位置示意圖(例題十)...………..…100 圖4.30懸臂梁頂部受偏心分佈載重負荷-位移曲線圖(例題十)……...….101 圖4.31懸臂梁底部受偏心分佈載重負荷-位移曲線圖(例題十)...……….102 圖4.32偏心力施加位置示意圖(例題十一)………..….……...103 圖4.33懸臂梁末端斷面頂部受集中載重負荷-位移曲線圖 (例題十一)………104 圖4.34懸臂梁末端斷面底部受集中載重之負荷-位移曲線圖 (例題十一)………105
圖4.35懸臂梁末端斷面中間受集中載重之負荷-位移曲線圖
第一章 緒論
長久以來薄壁梁結構被廣泛應用在工程上，隨著使用上的不同對 結構的選擇也有所不同，有些梁結構在使用中，需能承受大變形與大旋轉； 有些結構如飛機、太空梭、船舶等，為了減輕結構整體重量而使用高強度 的材料及薄壁斷面；有些特殊結構為了兼顧重量及盡量使受力後變形小， 也會使用薄壁梁來組成。 在結構力學上，Euler–Bernoulli梁的變形假設較Timoshenko梁單純，因 此被廣泛應用在建模上。尤拉梁與Timoshenko梁變形假設得最大差別在於 尤拉梁的斷面在變形後仍然垂直梁的中心軸，而Timoshenko梁考慮斷面剪 應變所以較複雜。當梁的細長比(Slenderness ratio) I AL2
較大時，梁的側 向位移主要是由撓曲正應變造成的，剪應變的貢獻相對地很小，因此細長 梁不論用Timoshenko梁模型還是Euler–Bernoulli梁模型分析，得到的側向變 形都差不多；而粗短梁的細長比很小，梁的側向位移就不能忽略剪應變的 影響。開口薄壁梁的扭轉剛度通常比撓曲剛度小很多，當兩個斷面主軸之 二次矩相差較大的開口薄壁梁，在垂直於強主軸的方向受側向力作用時， 即使該方向的細長比不大(如
10)，還是很容易引起側向－扭轉挫屈 (Lateral-torsional buckling)。在這種情況下，梁的挫屈負荷及在挫屈前後的 側向位移，也許不能忽略剪應變的影響。據本人所知，甚少有文獻探討剪 應變對薄壁梁之挫屈負荷及負荷－位移曲線的影響。所以本文希望探討Timoshenko梁與Euler–Bernoulli梁之幾何非線性行為及挫屈負荷的差異。文 獻[1]以一致性共旋轉法推導三維Euler–Bernoulli梁元素，考慮聖維南定理與 扭轉撓曲的影響，並以例題探討薄壁彈性梁側向扭轉挫屈的現象。文獻[2,3] 研究三維薄壁梁的材料及幾何非線性行為，文獻[2,3]採用Timoshenko梁與 Euler–Bernoulli梁兩種模型，其例題主要在探討彈塑性梁的幾何非線性行 為。但其例題之梁的細長比都很大，所以兩種梁模型的結果都差不多。文 獻[4]提出一三維Timoshenko梁的變形機制，推導一梁元素，並將其用在旋 轉梁的穩態變形及自然振動分析，但文獻[4]推導的梁元素並沒有考慮剪力 修正係數，且文獻[4]的穩態變形僅考慮軸向扭轉及位移。文獻[5]提出一個 決定剪切係數(shear coefficient)的方法，並與Timoshenko[6]與Cowper[7]的剪 切係數做比較。文獻[8]推導一三維Timoshenko梁元素的shape function。 本文將利用文獻[4]黃之三維Timoshenko梁的變形機制推導一兩個節點 十四個自由度的梁元素，但在推導梁元素時加入剪力修正係數。本文以 d'Alembert原理、虛功原理配合幾何非線性梁理論的一致線性化，推導出元 素之節點變形力。本文將利用文獻[9]提出的方法推導元素之切線剛度矩 陣。本研究將以數值例題探討三維薄壁梁結構在不同邊界條件、斷面下受 不同負荷的幾何非線性行為及側向挫屈的現象，以說明本研究中提出的方 法的正確性，並與Euler–Bernoulli梁的結果比較，以探討剪應變對薄壁梁之 挫屈負荷及負荷－位移曲線的影響。
第二章 理論推導
本章所推導的三維Timoshenko梁元素之基本假設、座標系統、梁的變形機 制、梁元素的推導方法與文獻[4]類似。 2.1 基本假設 本文對梁元素的推導，做如下的假設： (1) 梁為等斷面、雙對稱的Timoshenko梁。 (2) 梁元素的形心軸之單位長度伸長量(unit extension)為均勻的伸長。 (3) 梁元素的變形與應變皆為小變形與小應變。 (4) 梁元素斷面的翹曲為梁元素的軸向扭轉率與該梁的聖維南(Saint Venant) 翹曲函數的乘積。 2.2 座標系統描述本 研 究 是 使 用 共 旋 轉 有 限 元 素 法 (co-rotational finite element formulation)，將梁分割成若干個兩個節點的梁元素。為了描述梁系統的變 形，本文中使用三個右手座標系統： (1) 總體座標系統XiG,( ＝1,2,3） i 如圖四所示，系統的節點座標及方位，其他座標系統之座標軸的方 向餘弦，皆在此座標系統中定義。 (2) 梁斷面座標系統xiS,(i＝1,2,3）
該座標系統的原點是剛接在梁斷面的形心上，其x1S軸取在未翹曲斷 面的法線方向，x2S、x3S軸取在未翹曲斷面的主軸方向。 (3) 元素座標系統xi , ( ＝1,2,3） i 元素座標系統是建立在每個元素當前的位置上，如圖 2.1 所示，元 素座標系統的原點是定義在元素節點 1 上， 軸的方向為梁元素兩節點 連線的方向， 與 軸在元素變形前與斷面的主軸方向一致，而元素變 形後的 與 軸，可以由該元素未翹曲的兩端斷面的方位來決定[10]， 本文是分別將位於節點 1、2 後的斷面繞一個與該斷面之法線及與 軸 垂直的旋轉軸旋轉一角度使斷面之法線方向與 軸方向一致(此時不考 慮斷面之翹曲變形，否則斷面的法線方向將無法定義)，然後再以兩斷面 的主軸方向的角平分線作為 軸與 軸的方向。本文中梁元素的位移、 變形，均在此座標系統定義。 1 x 2 x 3 3 x 2 x x 1 x 1 x 2 x x3 本文中以符號
代表行矩陣。總體座標系統XG S 3 { ， 元素座標x=
，元素斷面座標 ={ }的關係可表示如下: G G G X X X1 , 2 , 3 }
3 2 1,x ,x x xS x1S,x2S,x XG AGEx XG AGSxS (2.2.1) 其中 、 分別代表元素座標、元素斷面座標對於固定總體座標系統 的方向餘弦矩陣。 GE A AGS 2.3 旋轉向量本文中使用旋轉向量來表示一個有限旋轉，如圖2.2 所示，一向量b受 到一旋轉向量ψ{
1
2
3}
e的作用而轉到一個新的位置 ，向量與b之間的關係可表示成[11]:
b b
bcos
b(1cos
)(eb)esin
(eb) (2.3.1)其中I為33的單位矩陣，符號．與分別代表向量的內積與外積，
表繞 旋轉軸的轉角，e表旋轉軸的單位向量。 當旋轉向量ψ有一微小變化
ψ時會使向量b繞 軸作微小 旋轉 i x (i1,2,3) i
，
ψ與
φ有以下的關係[11]
φΓ(ψ)
ψ (2.3.2) Γ(ψ)[Ia1(ψI)b1ψ(ψI)] (2.3.3) 其中 2 1 cos 1
a ， 1 (1 sin ) 2 1
b ，當
0， 2 1 1 a ， 6 1 1 b 當
0，將(2.3.3)式線性化，Γ(ψ)可近似成 0 0 0 2 1 1 2 1 3 2 3
I T (2.3.4) 2.4 Timoshenko 梁的變形描述 本文在元素座標上描述梁元素的變形，由(2.1)節中的基本假設可知 Timoshenko 梁元素的變形可以由其形心軸的位移、截面的翹曲(warping)及其截面的旋轉來描述。 本研究採用[4]中之梁的變形機制，如圖 2.1 中 Q 點為梁元素中的任意 點，P 點為 Q 點在同一斷面之形心軸上的對應點。在元素座標上 Q 點的變 形前後位置可以表示如下: r0 xe1 ye2 ze3 (2.4.1) r xpe1 ve2 we3
1,xωe1S ye2S ze3S (2.4.2) 其中x 、 、y z 為變形前 Q 點在元素座標 (xi i1,2,3)上的座標， x 亦為 P 點 變形前x1軸的座標， 、y z 亦同時是 點在Q x2S與x3S軸的座標。x (p x)、v(x) 以及w(x)分別是變形後 點在元素座標P xi(i1,2,3)上的座標，v(x)、w(x)亦 為 點在P x2及x3軸方向的位移， x ,x 1 1 是梁斷面沿變形後形心軸的軸向扭 轉率，
1(x)為形心軸的扭轉角，
(y, z)代表等斷面梁的聖維南翹曲函 數， 及ei eSi (i1,2,3 )分別為xi與xiS軸的單位向量。 令 s 為變形後形心軸的弧長， 為弧長的微小變形，由畢氏定理可以 知道： ds 2 2 2 2 dx dv dw ds p (2.4.3) 由(2.4.3)式可得 dx w v t u x xp
x x x 0 2 1 2 , 2 , 2 0 1( ) [(1 ) ( ) ( ) ] ) (
(2.4.4) 其中u1(t)為節點1在x1方向上的位移，由元素座標系統的定義，其值為零。 x v v x , ， x w wx , ， 0 dx 1 ds
(2.4.5)0
為形心軸的單位伸長量。 由 梁 元 素 的 變 形 為 小 變 形 的 假 設 ， 利 用 近 似 式 2 / 1 2 , 2 , 2 0) ] 1 [(
vxwx ) 2 1 2 1 1 ( 0 v,2x w,2x
，上式可表示成
x x x p x u t v w dx x 0 2 , 2 , 0 1 ) 2 1 2 1 1 ( ) ( ) (
(2.4.6) 由座標系統的定義可知，在變形前xi軸與xiS(i1,2,3)軸的方向是一致 的，即 與 的方向是一致的。在本文中假設變形後的單位向量 的 方 向 是 由 以 下 兩 個 旋 轉 向 量 連 續 作 用 於 單 位 向 量 來決定： i e 2, , 1 2, , 1 3) 2, , 1 (i S i e 3) 3) (i S i e i e (i 1
1e1 (2.4.7) m
mm
0
2
3
(2.4.8) {0, , } {0, 2, 3} m 3 2 m m m
m (2.4.9) 其中m 為垂直於 軸之單位向量，x1 2 1/2為 和 的夾角， 3 2 2 ) (
m e1 e1S
1為 斷面繞x1軸的轉角。 將旋轉向量1作用在 上，使其轉至一中繼位置ei e，再將i m作用在ei， 將其轉到eiS。若 、ei m、以及t已知，則元素斷面座標 就唯一決定；反 之，若 與 已知，則旋轉向量 S i e i e eiS m與1亦唯一決定。 eiS與 之關係可表示如下ei [9]、[10] (2.4.10) i i S i n R R e Re e [ 1 2] R1 cos
1r1sin
1r2 R2 sin
1r1cos
1r2 } } ) cos 1 ( ) cos 1 ( cos sin { 3 22 2 3 1 m
m
m
m m
m m m r ) cos 1 ( cos ) cos 1 ( sin { 2 2 3 32 2 m
m
m m m
m
m m r 其中R 稱為旋轉矩陣。因 R 為
i(i1,2,3)的函數，所以本文中稱
i為旋轉 參數，{
1,
2,
3} ，為旋轉參數向量。 當
i(i1,2,3)分別有一微小變化
i時，斷面座標會旋轉到一個新的 位置，此一新的位置可由元素座標繞xi(i 1,2,3)軸分別作微小旋轉
i (i1,2,3 )而得。 } , , {
1
2
3
與
φ{
1,
2,
3}之關係可表示如下[9]、[10]：
φ[n, t1an, t2 bn] T (2.4.11) } , , { 2 3 2 2 2 2 2 3 1 m m m
t } , , { 2 2 3 2 3 2 2 2 2 m m m
t 2 3 m a
22 m b
(2.4.11)式之反函數可表示如下:φ T φ T φ 1 2 2 3 2 2 2 3 2 2 2 3 2 2 2 3 2 2 3 2 3 ) sin( ) sin( ) cos( ) sin( )) sin( ) cos( ( ) sin( )) sin( ) cos( ( ) sin( ) cos( ) sin( ) 1 ) (cos( ) sin( ) 1 ) (cos( ) sin( 1 m m m m m m m m m m m m m m m m m m m m m m m m m m (2.4.12) 當旋轉參數向量{
1,
2,
3} 很小時，T矩陣可近似如下式 1 0 0 1 1 2 3 2 2 1 3 2 1
T (2.4.13) 將(2.4.10)式代入(2.4.2)式，利用近似式 22 32 2 1 2 1 1 cos
m
並保留 變形參數至二次項,則位置向量r 可以化簡成 3 3 2 2 1 1e e e r r r r (2.4.14) r1 xp y(
3
1
2)z(
2
1
3)
1,x
2 1 2 3 1, 3 3 2 1 2 ) 2 1 ( ) 2 1 2 1 1 ( z x y v r
2 1, 2 2 2 1 3 2 1 3 ) 2 1 2 1 1 ( ) 2 1 ( z x y w r 由梁之形心軸單位長度的伸長量為均勻的伸長量之假設及(2.4.6)式，可 以得到形心軸單位長度伸長量
0可表示如下
L v x wx dx L L L l 0 2 , 2 , 0 ( ) 2 1
(2.4.15) l Lu2u1 (2.4.16)其中L 為梁元素變形前的長度，l為梁元素變形後之形心軸的弦長， 、 分別為節點1與2在 方向的位移。 1 u u2 1 x 本文中假設梁元素形心軸的側向位移v (x)、w(x)及旋轉參數
1(x)， 2
(x)，
3(x)的形狀函數如下[8]:
t b b t b b b N N N v N x v( ) 1 2 b3 4}
31 2
32 N u c3 v1 w1 { ) (2.4.17)
t c c t c c c N N N w N x w( { 1 2 4}
21 2
22
N u (2.4.18)
t d d d d d N N N x 1 2 4 11 1 12 2 N u 1( { }
t
w d N 3 e N 3 ) ) (2.4.19)
t c e t e e e N N w N x 1 2 4 1 21 2 22 N u 2( { }
(2.4.20)
t b f f f f N N v v N x 1 2 3 4 1 31 2 32 N u 3( {
}t f N ) i u (2.4.21) 其中 (i bd)為梁在元素座標上的節點位移 ) 1
)) ( 2 ( 2 3 ( )) ( 2 2 ( )) ( 3 1 ( 2 2 3 2 3 2 3 2
z z z z z z z z d LD d D d LD d D N 2 3 2
b (2.4.22) ( 1 (
y y )) ( 2 ( ) 2 3 ( )) 2 2 ( )) 2 3 1 ( 2 3 2 3 2 3 3
y y y y y y d LD d D d LD d D N 2 2 2
) ) 3 3 3 3
c d (2.4.23) ( 2 3 2 ( 2 3 1 2 2 2 2
L L N (2.4.24) ) 3 2 ( ) ( 6 )) 1 ( 3 4 1 ( ) ( 6 2 2 2 2
y y y y y y e d D L D d D L D N (2.4.25) ) 3 2 ( ) ( 6 )) 1 ( 3 4 1 ( ) ( 6 2 2 2 2
z z z z z z f d D L D d D L D N (2.4.26) ) 1 ( 1 z z d D 12 2 L GA EI d sy y z z
(2.4.27) ) 1 ( 1 y y d D 2 12 L GA EI d sz z y y
L x
其中Ni(i a f )為Timoshenko梁的形狀函數(shape function) [8]。 與 ( j v j w j= 1, 2)分別是v與 在節點w j 的節點值。
ij( =1-3,i j=1, 2)是
i在節點 j 的節點值，
j(j 1,2)則是 x x 1 , 1
在節點 j 的節點值。
y、
z為 、 方向的剪力修正係數， 、 為 、 方向的有效剪力面積。 S x2 x3S sy A Asz x2S x3S 將(2.4.16)至(2.4.18)式代入(2.4.15)式整理可得 ) 2 1 2 1 ( 1 0 c t c b t b a t a L G u G u G u
(2.4.28) } 1 1 { a G (2.4.29)dx v G G G Gb b b b b x b { 1, 2, 3, 4}
N , G (2.4.30) dx w G G G Gc c c c c x c { 1, 2, 3, 4}
N , G (2.4.31) } , {u1 u2 a u (2.4.32) 將(2.4.28)代入(2.4.6)整理可得 p x Ntaua x L x 2 ( c) t c b t b u G u G
x vx wx dx 0 2 , 2 , ) ( 2 1 (2.4.33) } , 1 {
a N (2.4.34) 2.5 Tmoshenko梁的應變 為了推導上的方便，本節中我們將先推導出梁的應變以及梁的位置向 量變分、應變的變分。 2.5.1 梁的應變 假如將(2.5.1)式中的x、 、y z 視為拉格蘭日座標(Lagrange coordinates) ，則Green strains11、12、13可以表示成[20]: ( 1) 2 1 1 t 1 11 g g
， 12 1t 2 2 1 g g
， 13 1t 3 2 1 g g
(2.5.1)z
y
x
r
g
r
g
r
g
1,
2,
3 (2.5.2) 將(2.4.6)、(2.4.14)代回(2.5.2),可得gi的分量g (ij i, j1,2,3)如下 ) ( ) ( 2 1 1 ,2 ,2 11 o vx wx y 1 2,x 2 1,x 3,x g
(2.5.3) z(
1
3,x
3
1,x
2,x)
1,xx
x x xx x x x x x x y z v g , 3 , 1 3 , 1 , 3 2 3 , 2 , 1 , 3 3 , 1 1 , 12 ) 2 1 2 1 ( ) (
x x xx x x x x x x y z w g , 2 , 1 2 , 1 , 2 2 , 1 1 , 3 2 3 , 2 , 1 , 13 ) ( ) 2 1 2 1 (
g21
3
1
2
1,x
,y 22 12 32 1, , 3 2 1 2 1 1
x
y
g 2 , , 1 3 2 1 23 2 1
x y g g31
2
1
3
1,x
,z 32 1 2 3 1, , 3 2 1
x z g 2 , , 1 2 2 2 1 33 2 1 2 1 1
x
z
g 為了能清楚分辨出應變中那些項代表斷面剪應變，令 0 , 2 1
wx ds dw 0 , 3 1
vx ds dv (2.5.4) w,x
2(1
0) v,x
3(1
0) 再將(2.5.4)代入(2.5.2)式及(2.5.3)式再代回(2.5.1)式，保留變形參數及其微分 到二次項，並在剪應變加上修正係數
12
y，
13
z，其中
y、
z[18] 分別為(2.4.27)式中x2S、x3S方向的剪力修正係數，可得: 2 11 1 11 11
(2.5.5-a) xx x x y z 2, 3, 1, 0 1 11
2 , 1 2 , 3 1 , 2 , 1 2 , 2 2 , 3 , 2 , 3 , 1 , 2 1 2 , 3 2 , 1 12 3 3 13 2 2 2 , 1 2 2 0 , 3 , 2 , 1 2 0 2 11 2 1 ) ( 2 1 ) ( 2 1 ) ) ( ) (( ) ( 2 1 ) ( 2 1 xx x x xx x x x x xx x x x x x x xx z z y z y z y z y y z 2 12 1 12 12
(2.5.5-b) 12 3 3 , 1 , 1 12 ( ) 2 1 ) ( 2 1
y z x ) ( 4 1 ] ) ( ) ( ) ) ( [( 2 1 3 , 2 , 3 2 , 3 , 1 , , 2 , 1 , , 1 , 1 , 1 13 2 2 0 12 3 3 , , 1 2 12 x x x x y x x y x xx y y x z y z 2 13 1 13 13
(2.5.5-c) 13 2 2 , 1 , 1 13 ( ) 2 1 ) ( 2 1
z y x ) ( 4 1 ] ) ( ) ( ) ) ( [( 2 1 2 , 3 3 , 2 , 3 , 1 , , 2 , 1 , , 1 , 1 , 1 12 3 3 0 13 2 2 , , 1 2 13
x x x x z x x z x xx z z x y y z ) 2 , 1 , 3 , 2 , 1 ( 1 j k k j
代表 中之 次項。(2.5.5)式中加底線的項為梁的斷面 剪應變項，在這些項中都加入修正係數 k j 1
k 12
或
13。 2.5.2 梁之位置向量的變分及應變的變分 本文利用虛功原理及d’Alembert原理推導平衡方程式，所以需要求出 位置向量及應變的變分。由(2.4.14)及(2.4.33)式，我們可以求得位置向量變 分的分量 } , , {
r1
r2
r3
r (2.5.6)
r1 xp 1(z 3 y 2) 2(z y 1) 3(y z 1) 1,x v z y z y z r x x
) ( ) 2 1 ( ) 2 1 ( ) ( 3 , 1 , 1 2 3 3 3 2 1 1 2 w y z y z y r x x
) ( ) 2 1 ( ) 2 1 ( ) ( 2 , 1 2 3 , 1 2 3 2 1 1 3 其中 p x
v v w w dx L x x x x x x c t b t a t a (
b
c ) 0( ,
, ,
, )
u N u G u G
(2.5.7) 而由(2.5.5)及(2.4.28)式可以得到應變的變分 2 11 1 11 11
(2.5.8-a) xx x x y z 2, 3, 1, 0 1 11
) ( ) ( ) ( ) ) ( ) ( ( ) ( ) ) ( ( ) ( ) ) ( ) (( ) ( ) ( , 1 13 , , 1 12 , , 1 2 , 2 , 3 0 , 1 12 3 12 , , 1 2 13 , 13 2 2 , 1 , 1 , 3 , 2 1 1 , 1 , 2 0 , 3 2 , 3 , 1 12 3 , 1 0 , 2 2 , 3 1 , 2 , 1 13 2 , 1 , 3 , 2 0 0 2 11 x x x x xx x x xx x x x x x x x xx x x x x xx x x x x xx x x y w z v z y z v z y w y z y z y z y z z z y z y y z
2 12 1 12 12
(2.5.8-b) ) ) 1 ( 2 1 ( ) 2 1 2 1 ( ) 2 1 ( ) 2 1 ( 12 , 3 12 1, , , 0 12 0 1 12
vx x y z v x ) 2 1 ( ) 2 1 ( ) 2 1 ( ) 2 1 2 1 2 1 2 1 2 1 ( ) 2 1 2 1 ( ) 4 1 2 1 2 1 ( ) 4 1 2 1 ( ) 4 1 2 1 ( ) 4 1 2 1 ( ) 2 1 2 1 2 1 ( 1 13 , 12 0 , , 1 , , 1 , 3 , 2 , , 1 , 0 , , 3 , , 1 13 , 13 2 1 2 , 1 , , 1 , 3 , 2 12 0 3 3 , , 1 , 2 , 3 1 13 2 12 , 12 3 , , 1 0 2 12
x x x y xx x x y xx y y x y x x x y x x x y x x x x y x w v z y w z y z z z z v 2 13 1 13 13
(2.5.8-c) ) ) 1 ( 2 1 ( ) 2 1 2 1 ( ) 2 1 ( ) 2 1 ( , 13 2 13 1, , , 0 13 0 1 13
wx x z y wx ) 2 1 ( ) 2 1 ( ) 2 1 ( ) 2 1 2 1 2 1 2 1 2 1 ( ) 2 1 2 1 ( ) 4 1 2 1 ( ) 4 1 2 1 ( ) 2 1 2 1 4 1 ( ) 4 1 2 1 ( ) 2 1 2 1 2 1 ( 13 0 , 1 12 , , 1 , , 1 , 2 , 1 , 0 , , 2 , , 3 , , 1 12 , 12 3 1 2 , , 1 , 3 , 2 1 12 3 , 1 , , 1 3 , 2 , 3 13 0 2 13 , 13 2 , , 1 0 2 13
x x x z xx x xx z z x z x z x x z x x x x z x x x x z x w v z y v y y y z y y w 其中 ) ( 1 0 c t c b t b a t a L
u G
u G
u G
(2.5.9) 2.6 節點參數與節點力 2.6.1 系統節點參數與節點力 本文中在元素節點 j ( j =1, 2)上的系統節點參數為UijG與
ijG及
j，其中 是節點 G ij U j的平移向量UGj 在固定總體座標軸XiG(i1,2,3)方向的分量；G ij
是節點 j的旋轉向量ψGj 在固定總體座標軸XiG(i1,2,3)方向的分量； j
是節點 j的剪心軸扭轉率。本文中在任何變形位置的 都被重新設定 為零，所以旋轉向量分量的變分 與繞 軸的微小旋轉 的值相同[12]。 G j ψ ) 3 , 2 , 1 (i G ij
XiG
ijG 2.6.2 元素節點參數與節點力 本文中採用了兩組元素節點參數，分別稱為顯節點參數(explicit nodalparameter)與隱節點參數(implicit nodal parameter)。元素的顯節點參數是用
來將元素矩陣組合為系統矩陣，所以其定義必需與系統的節點參數一致。 本 文 中 在 元 素 節 點 j ( j1,2)的 顯 節 點 參 數 為uij 與
ij 及
j ， 其 中 j ij u (u1j uj,u2j vj,u3 j wj)是節點 的平移向量 在其當前固定元素座 標軸 ( i )方向的分量； j u i x 1,2,3 ij是節點 j的旋轉向量 在當前固定元素 座標軸 ( )方向的分量； j ψ j i x i1,2,3
j是節點 的剪心軸扭轉率。本文中在任 何變形位置的 都被重新設定為零。對應於 的廣義節點力 ，為在 軸 方向的力；對應於 j ψ uij fij xi ij 的廣義節點力 ，與繞 軸的傳統力矩的值相同。 對應於 ij m xi j
的廣義節點力Bj為廣義雙力矩 。 元素的隱節點參數是用以決定梁元素的變形。在元素節點 j ( j1,2) 的隱節點參數為uij與
ij及
j，其中uij (u1 j uj,u2j vj,u3j wj)是節點 j 的位移向量uj在其元素座標軸xi (i1,2,3)方向的分量；
ij是旋轉參數
i在節點 j的值，
j是節點 j的剪心軸扭轉率。對應於 的節點力為在 軸 方向的力 ；對應於 及 ij u xi ij f
ij
j的節點力分別為廣義力矩mij與廣義雙力矩 。因為 在變形後不為零，所以其變分並不是繞 軸的無限小旋轉， 所以廣義力矩 j B
ij xi ij m 3 , 2 , 1 並非繞xi軸的傳統力矩。 2.6.3 元素顯節點參數的擾動與隱節點參數的擾動關係 令 x ii ( )表示在當前變形位置的固定元素座標，ei 、
ij 、 ij
(i 1,2,3, j1,2)、 分別表示在 軸方向的單位向量、元素節點xi j的旋 旋轉向量及轉參數分量、及元素的弦長。此時元素的顯節點參數除了 外其餘的值皆為零。令 u2 L (2.6.1) } , β q { u1, φ1 u2,
φ2,
j
} j , j v , j { u
w
j
u , 2 j , {
1 j
3j}
φ , 2} {
1
其中
q 為一擾動節點位移向量，( j1,2)，
uj、
vj、
wj為節點 j在固 定元素座標xi(i1,2,3)軸方向的擾動位移，
ij為節點 j繞固定元素座標 ( i )軸的擾動旋轉， i x 1,2,3
為擾動扭轉率向量，
j為元素扭轉率在節 點 j的擾動量。 令
q 為對應於
q 之顯節點參數向量 q的擾動量，
q 可表示成 } , , , , { u1 ψ1 u2 ψ2 β q
(2.6.2) } , , { 1j 2j 3j j
ψ 其中
uj、
在(2.6.1)式中已有定義。 由(2.3.2)式可知
q 與
q 有以下的關係
q T q (2.6.3) 其中矩陣T的顯式列於附錄A。 因本文中在任何變形位置的ψj都被重新設定為零，所以T 的值等於I14。 若梁元素在當前的變形位置受到
q 的擾動 ，則可由文獻[13]的方法決定 受擾動後的元素座標x (i i1,2,3)、定義在擾動後的元素座標之元素節點 j 的旋轉參數
ij(i1,2,3, j1,2)及元素的弦長 。若取到擾動量的一次 項，則擾動前後的元素座標的關係可表示成 [13] 3 2 x x x x1 1 2 12 11 12 2 12 11 1 12 12 12 1 2 ) ( 2 ) ( 1 1 A A w A A v w v
x Axx x x 3 2 x 1 (2.6.4) 其中
v12
v2
v1,
w12
w2
w1, ( ) 2 1 12 12 11 11 1
A , ) ( 4 1 32 21 31 d
32 2 d A 22d
21
31d
22
, 12 2 2 w j d j
, 12 3 3 v j d j
. 若取到擾動量的一次項，則擾動後元素的弦長 可表示成12 u
(2.6.5) 1 2 12 u u u
擾動後元素節點 j 的旋轉參數
ij(i 1,2,3, j 1,2)在擾動後的元素座標x 可i 表示成 4 ) ( 4 ) ( 2 2 2 2 32 31 12 22 21 12 12 11 12 11 11 12
v w 4 4 4 4 32 22 31 21 22 32 21 31
(2.6.6) j j j j j j j w v 3 11 3 12 3 1 12 11 12 2 12 2 2 2 1 2 1 2 ) (
2 ) ( 11 12 3
j j j j j j j j w v 1 2 12 2 11 2 12 11 12 3 12 3 3 2 1 2 1 2 ) (
2 ) ( 11 12 2
j 其中
11
12，所以
11
12 0，但其變分不為零，且在推導元素剛度矩 陣時需用到該變分，故在上式中仍需保留該項。 當元素的顯節點參數受到
q 擾動時，其運動過程可視為先做一剛體 運動後再做一變形運動，在該剛體運動後，元素從當前的變形位置移動到 擾動後的元素座標x ，但仍保持其原來的變形，即定義在i x 座標的元素節i 點旋轉參數及元素的弦長和擾動前定義在xi座標的元素節點旋轉參數及元素的弦長有相同的值，仍為
ij和 。在該變形運動後，定義在x 座標的i 元素節點旋轉參數及元素的弦長就變成(2.6.5)式及(2.6.6)式之
ij及 。所 以由(2.6.5)式及(2.6.6)式可知在x 座標中元素節點旋轉參數及元素的弦長i 受到的擾動可表示成 ij ij ij
(2.6.7)
(2.6.8) 由元素座標的定義可知元素節點位移在x 方向的擾動可表示成i 0 1 vj wj u
(2.6.9)
u2
令 } , , , , {
1
1
2
2
q u u (2.6.10) } , , { j j j j
u
v
w
u } , , { 1j 2j 3j j
其中
在(2.6.1)式中已有定義。 由(2.4.15)式、(2.4.28)式-(2.4.31)式、(2.5.9)式、(2.6.5)式-(2.6.10)式可知
q 與在x 座標中隱節點參數的擾動量i
q 有以下的關係
q T q (2.6.11)其中矩陣T 的顯式列於附錄A 。 由(2.6.3)式及(2.6.11)式可知
q 與顯隱節點參數的擾動量
q 有以下的關 係
q T T q (2.6.12) 2.7 元素節點內力f之推導 對應於(2.6.11)式及(2.6.12)式之
q 、
q 與
q 的節點內力向量分別為 顯節點內力向量f {f1,m1,f2,m2,B} } , 2 B m 、傳統節點內力向量 , , , {f1 m1 f2 f 與隱節點內力向量f {f1,m1,f2 ,m2,B}，其中 } 3 , , { 1 2 j { 1 , 2 , j j j f f f f ，fj fj fj f3j}，fj {f1j, f2j, f3j}， } 3 , , { 1 2 j m { 1 , j j j m m m ，mj mj m2j,m3j}，mj {m1j,m2j,m3j}( j = 1, 2)， ， } , {B1 B2 B fij為作用在元素節點 j，xi方向的內力， f1j為作用在元素節 點 j ， 方向的內力，xi mij 為作用在元素節點 j ，繞 軸的力矩、xi mij 為作用 在元素節點 j 的廣義力矩之 分量，xi mij為作用在元素節點 j 的廣義力矩之xi 分量，Bj為作用在元素節點 j 的雙力矩(Bimoment)。 及 的分量是定義 在當前的元素座標上，而 f f f 的分量是定義在擾動後的元素座標上，雙力矩 和座標系統無關。 本文利用虛功原理推導元素節點內力，本文將元素的節點內力 及 視為作用元素節點的外力，應力所作的虛功視為是隱節點內力 f f f 所作的虛 功。在元素當前的變形位置給元素節點 j( j1,2)一個虛位移，則由虛功原理可知
