行政院國家科學委員會專題研究計畫 期中進度報告
二階梁理論及其在側向-扭轉挫屈的應用(1/2)
計畫類別： 個別型計畫 計畫編號： NSC92-2211-E-009-046- 執行期間： 92 年 08 月 01 日至 93 年 07 月 31 日 執行單位： 國立交通大學機械工程學系 計畫主持人： 蕭國模 計畫參與人員： 陳弘虎、李政道、陳致中 報告類型： 精簡報告 報告附件： 出席國際會議研究心得報告及發表論文 處理方式： 本計畫可公開查詢中 華 民 國 93 年 5 月 27 日
二階梁理論及其在側向-扭轉挫屈的應用
(1/2)
A Second Order Beam Theory and Its Application In Lateral-Torsional
Buckling (1/2)
計畫編號：NSC 92-2211-E-009-046
執行期限：92 年 08 月 1 日至 93 年 07 月 31 日
主持人：蕭國模 國立交通大學機械工程學系
計畫參與人員：陳弘虎、
李政道、陳致中
中文摘要：本研究利用二階梁理論(Second order beam theory)之平衡方程式及構成方程式，對 一梁在扭矩及軸力作用下的幾何非線性側向挫屈分析。本研究採用共旋轉法(corotational formulation)來描述結構的變形及建立平衡方程式，即在分析時將梁結構分成很多小段， 每一段稱為梁元素。然後在每個梁元素變形後的最新位置上建立一個元素座標系統。梁 元素的變形及平衡方程式都是在元素座標上定義。本研究考慮了彎曲、扭轉、及軸向變 形間完整的非線性耦合。本研究先求出梁在軸力及扭矩作用下之主要平衡路徑的解析 解，然後在主要平衡路徑上加一擾動位移，並以一致線性化推導出次要平衡路徑之擾動 位移的統御方程式。本研究提出一級數解法求得擾動位移的正確一般解，再由邊界條件 及內部節點的平衡求得其挫屈狀態的統御齊次方程組。本研究採用二分法求挫屈扭矩。 本研究以數值例題來驗證本文中所提出的方法的有效性。 關鍵詞: 梁, 幾何非線性, 共旋轉法, 虛功原理, 側向-扭轉挫屈 Abstract.
In this paper, the geometrical nonlinear lateral buckling analysis of a beam subjected to torsional moment and axial force is studied using the equilibrium equations and constitutive equations of second order beam theory. A corotational formulation is used to describe the deformation of beam structure and establish equilibrium equations in this paper. The beam structure is divided into many small parts, named beam element. A set of element coordinate system is constructed at the current configuration of the deformed beam element. The deformation and equilibrium equations are defined in terms of the element coordinates. The complete nonlinear coupling among bending, twisting, and stretching deformations is considered here. The analytical solution of the primary path is solved first. Then perturbation displacements are added to the obtained solution of primary path. The new obtained, displacements are substituted into equilibrium differential equations and constitutive equations of the secondary beam theory and using first order linearization to derive the governing equations for the perturbation displacements of secondary equilibrium path. A power series solution method is proposed to solve the perturbation displacements. From the boundary conditions and the equilibrium condition at the internal nodes, the governing homogeneous linear algebra equations for the buckling state can be obtained. A bisection method is used to solve the buckling moment. Numerical examples are presented to investigate the effectiveness of proposed approach
.
Keywords: Beam, Geometrical Nonlinearity, Co-rotational Formulation, Virtual Work Principle, Lateral-Torsional Buckling
1. 緒論 當梁斷面為開口薄壁斷面時，因其撓曲、扭曲及軸向變形間的耦合效應容易引起挫 屈。因線性梁理論無法考慮這些耦合效應，薄壁梁的挫屈的統御方程式都需使用非線性 梁理論推導[1-9]，即使是線性挫屈分析的統御方程式[6]，也必需是由非線性梁理論的一 致線性化推導[4,6,7]才能得到正確的統御方程式及結果。為考慮各種變形間的耦合效應 的耦合效應，在非線性梁理論中最少要保留到變形參數的二次項。文獻中一般的非線性 梁理論[1-9]都僅保留到變形參數的二次項，故本文稱之為二階梁理論(second order beam theory)。但是文獻上有關變形間的耦合效應大都不是以梁之正確的變形機制(Kinematics) 及二階的一致線性化(consistent second order linearzation)推導[4,6,7]，文獻上有些二階梁 理論僅就幾何變形上所能觀察到會影響應變或平衡方程式的部分加入原先傳統之線性 推導中，或是在推導的過程中太早對梁的變形機制作線性化，故文獻中梁的平衡微分方 程式有很多都無法考慮到完整的耦合效應，僅適用於特定的問題，而且文獻上對扭矩與 扭率的關係也都採用線性分析所採用的構成方程式(constitutive equation)而不是由二階 梁理論推導出來的，所以文獻上的方法有很多無法求得含完整耦合效應之統御方程式， 故其分析的結果和結論有很多是不合理的或是錯誤的。 文獻上似乎仍沒有一個以完整的二階梁理論推導的平衡微分方程式及統御方程 式，故本人在[10]中用共旋轉法推導出一個二階梁理論之平衡微分方程式及構成方程 式，以添補文獻上的不足，並提供一個分析三維尤拉梁結構之各種非線性挫屈行為正確 及有效的理論工具。 梁同時受軸力與扭矩的情況在各種鑽探時的轉軸經常會發生，其側向-扭轉挫屈分 析很重要，但文獻上除了文獻[5]外皆為線性挫屈分析的結果，文獻[5]中雖然考慮了挫 屈前的變形，但採用線性分析所採用的扭矩與扭率的構成方程式，而不是由二階梁理論 推導的構成方程式，且文獻[5]沒有考慮扭轉率的三次向，故其結果的準確性需要以正 確的二階梁理論分析的結果來檢驗。故本研究的目的是將利用本人提出之二階梁理論的 平衡微分方程式及斷面合力與變形間的構成方程式探討三維梁同時受軸力與扭矩的側 向-扭轉挫屈分析，並以數值例題說明本研究的實用性及正確性，同時驗證文獻上之挫 屈扭矩之正確性。 2. 理論推導 本研究採用了[10]中提出的二階梁理論，為了完整性，本節將簡單的描述二階梁理 論。 2.1 基本假設 本文對非線性梁元素的推導，作如下的假設： (1) 為細長的（slender）等斷面（prismatic）梁，且Euler-Bernoulli 假說成立。 (2) 梁斷面為雙對稱斷面。 (3) 梁元素在斷面上沒有變形，且斷面內的應變可以忽略。 (4) 梁元素斷面的翹曲為梁元素的軸向扭轉率與該梁的聖維南（Saint Venant）翹曲函數 的乘積。 (5) 梁元素的變形為小變形。 由假設(1)與假設(4)可知梁元素的變形可由其形心軸的位移、斷面的方向(即斷面座 標)及其翹曲來決定。 2.2 座標系統
P X3 2 X X1 x3 2 x 1 x xc s 3 S y P x z Q xS 2 1 x u P v w ( ,0,λ 0) 2 x xS3 xS2 S 1
b
a
b
為了描述整個系統的變形，本文中定義了以下三個座標系統： (1) 固定總體座標系統， Xi (i = 1, 2, 3) (如圖.1)；系統的節點座標及方位，其他座標 系統之座標軸的方向餘弦，皆在此座標系統中定義。 (2) 元素座標系統，x
i (i = 1, 2, 3) (如圖.1)。 (3) 元素斷面座標系統，x
iS (i = 1, 2, 3) (如圖.1)；x
1S軸為未翹曲斷面的法線方向，x
2S 軸與x
3S軸分別與未翹曲斷面的主軸重合。元素的變形是由斷面座標相對於元素座標 的旋轉來決定。 圖.1 元素座標與元素截面座標 圖.2 旋轉向量 2.3 旋轉向量 本文中使用旋轉向量來表示一個有限旋轉，如圖 2.2 所示，一向量 b 受到一旋轉向 量φa的作用而轉到一個新的位置b ，向量 b 與b之間的關係可表示成[11] ) ( sin ) )( cos 1 ( cos b a b a b b= φ + − φ ⋅ + φ × (1) 其中符號．與×分別代表向量的內積與外積，φ表逆時鐘方向旋轉角，a表旋轉軸的單 位向量。 令ei與eiS (i = 1, 2, 3) 分別代表元素座標 xi方向的單位向量與元素斷面座標 xiS軸方 向的單位向量。由座標系統的定義可知，在變形前x 軸與i xiS軸是一致的，而且變形後e1S 與(2.3)式的t重合。在本文中假設變形後的單位向量eiS (i = 1, 2, 3)，是由以下兩個旋轉 向量連續作用於單位向量ei (i = 1, 2, 3)來決定 [4]: n θn =θn ， θt =θ1t, (2,3) 其中} , , 0 { } ) ( , ) ( , 0 { 3 2 2 1 2 3 2 2 2 2 1 2 3 2 2 2 n n = + + = θ θ θ θ θ θ n (4) } , , {cosθnθ2θ3 = t (5) cosθn = −(1 θ22 −θ32 1 2) (6) θ2 = − dw s ds ( ) , ds s dv )( 3 = θ (7,8) n為垂直 e1 和 t 的單位向量。θn為cosθn的反函數。v s( )與w s( )為形心軸的側向位 移。s 為形心軸在變形後的弧長。 由 (2)-(8)式可得ei 與 eiS (i = 1, 2, 3)的關係可表示如下 eiS t R R e Re i i =[ , 1, 2] = (9) 其中R稱為旋轉矩陣。因R為θi (i = 1, 2, 3)的函數，所以本文中稱θi為旋轉參數。 當θi(i=1,2,3)分別有一微小變化δθi時，斷面座標會旋轉到一個新的位置，此一新 的位置可由元素座標繞x ( i = 1, 2, 3)i 軸分別作微小旋轉δφi(i=1,2,3)而得[4] δ θ θ θ θ δ δ θ = φ φ. 1 2 2 1 0 0 1 3 2 3 2 1 − − =T− (10) 2.4 梁元素之位移與應變 圖. 1中Q點為梁元素中的任意點，P點為Q點在形心軸上的對應點，即P點與Q點位 於同一斷面上。在元素座標上，Q點在梁元素變形前後的位置向量可分別表示如下 r0 = xe1 + ye2 +ze3 (11) 和 r=xc( )s e1 +v s( )e2 +w s( )e3 +θ1,se1s +ye2s +ze3s (12) 其中 xc( )s = +x u s( ) (13) 其中xc(x)、v(x)、以及w(x)分別是P點在xi( i = 1, 2, 3)軸上的座標(見圖. 1)， x x d d 1 1, θ θ = 是沿變形後的形心軸的斷面之軸向扭轉率，ω =ω(y,z)代表等斷面梁的聖維南(Saint Venant)翹曲函數，y與z分別是Q點在x2S軸與x3S軸的座標。
本文中應變的度量是採用 Green strain。本文中以εij(i=1 ,2 ,3 ,j=1 ,2 ,3)表示 Green
ε11 1 2 1 = (r r,tx ,x − ) ， ε12 1 2 = r r,tx ,y ， ε13 1 2 = r r,tx ,z. (14) M 1 F 1 F F2 1 2 x x x3 2 M M1 3 M F3 a Bθ b θ B s 3 M F3 F2 2 M 1 圖.3 作用於梁元素中任一小段的端點負荷 2.5 平衡方程式及構成方程式 本文利用虛功原理在元素座標上推導梁元素的統御方程式及構成方程式。圖. 3 所 示為梁元素中的一小段，點a、b分別表示這一小段梁元素的端點。本文採用形心軸在
x
i (i = 1, 2, 3)軸方向的位移量ui (i = 1, 2, 3) ( u1=u, u2 =v, u3 =w)，旋轉參數θi (i = 1, 2, 3), 及 θ1,s作為廣義位移。在斷面上與廣義位移對應之廣義合力為在xi(i = 1, 2, 3) 軸方 向的力Fi (i = 1, 2, 3)，廣義力矩 Miθ (i = 1, 2, 3)及廣義雙力矩(Bimoment)Bθ。 廣義 力矩與傳統力矩Mi (i = 1, 2, 3)間的關係可以由(10)式及 Controgradient Law [13]求得 並可表示如下 M T M= −t θ (15) 其中 M={
M1 M2 M3}
,Mθ ={
M1θ M2θ M3θ}
,及 T−t 為(10)式中 T−1 的轉 置矩陣。 若給予端點a及b(見圖. 3)一個廣義位移(δui)j, (δθi)j , δθ1,s (i = 1, 2, 3, j = a,b)， 其中 j ) ( 表( )在端點j之值，則由虛功原理可知對應於虛位移，內力所做的虛功δWint等 於外力所做的虛功 ext W δ δWint =δWext (16) 因本文假設梁元素間無外加負荷，所以δWext即為端點a、b上內力之廣義合力所做的 虛功，δWext可以表示成( )
δWext tδ θ tδ Bθδθ s a b = + + F u M θ 1, (17) 其中 F={F F1, 2,F3}δu={δ δu1, u2,δu3}={δ δ δu v, , w} δθ ={δθ δθ δθ1, 2, 3} (18) (17)式表示 [ ]ab 內各項在端點b之值減掉各項在端點a之值。 由基本假設(3)，本文僅考慮ε11,ε12,ε13，則(16)式中內力所作的虛功可表示成
(
)
δWint σ δε σ δε σ δε dV a b =∫
11 11+ 12 12 + 13 13 (19) 其中δε1j (j = 1, 2, 3) 為(14)式中 ε1j 的變分。σ1j (j = 1, 2, 3)為工程應力，V為端點a , b間梁元素變形前的體積dV =dAds/(1+ε0)。又本文假設材料為線彈性材料，所以應力 與應變有如下的關係 σ11 =Eε11, σ12 =2Gε12, σ13 =2Gε13 (20) 其中 E 為楊氏模數，G為剪力模數。 將 (17), (19) 及 (20) 代入式 (16)，可得到下列方程式：(
)
M1θ,s = E Iy−Iz θ θ2,s 3,s (21) M2 s F3 F1 2 1C 1s 3s 2 , , , θ − = θ − θ θ (22) M3 s F2 F1 3 1C 1s 2 s 2 , , , θ + = θ + θ θ (23) F1,s =0, F2,s =0, F3,s =0 (24-26) M1 Bs E 0Ip 1s C 1 0 1s 1 2s 3 2 3s EI4 13s 2 1 2 1 2 θ + θ = ε θ + +ε θ + θ θ − θ θ + θ , , [( ) , , , ] , (27)(
)
M2 E 1 4 0 Iy 2s Iy Iz 1 3s 3 yz 1ss 3s 1C 1s 3 2 θ = [ + ε θ +( − )θ θ − θ θ ]+ θ θ , , Ω , , , (28)(
)
M3 E 1 4 0 Iz 3s Iy Iz 1 2s 3 yz 1ss 2s 1C 1s 2 2 θ = [ + ε θ +( − )θ θ − θ θ ]− θ θ , , Ω , , , (29)(
)
Bθ =C1 1+4ε θ0 1,ss −3EΩyzθ θ2,s 3,s (30) F1 E A 0 1 5 0 Ip 12s Iy 22s Iz 32s 12ssA C 12s 2 1 2 5 2 = [ ε ( + ε )+ θ , + ( θ , + θ , +θ , ωω)]+ θ , (31)其 中 C =GJ 為 扭 轉 剛 度 (torsional rigidity) 及 C1 = EAωω 為 翹 曲 剛 度 (warping rigidity)。 式 (21)-(26) 可視為平衡方程式，(27)-(31)式可視為構成方程式。在(21)-(31)式中 保留變形參數的全部一次項、二次項及部分三次項(含底線項)。在[10]中發現某些三次 項亦不能忽略，否則在某些情況中會造成誤差。所以本文亦保留[10]中的三次項。 2.6 主要平衡路徑之統御方程式 梁承受一定值之壓負荷P及保守扭矩M ，所以F1 =−P、θ2 =θ3 =0。將近似式 1 1+ε0 = 代入(31)式可得 2 1 1 2 1 0 2 1 ss , s , p p EA C I E C A θ θ ε ε − + − − = (32)
其中 AE P p = ε 。 由(21)式知 θ 1 M 為一常數(因θ2 =θ3 =0)，所以θ1,s為一常數，可得θ1,ss =θ1,sss =0。由 (15)式知M1θ =M 。因θ1,ss =θ2 =θ3 =0，由(30)式可得Bθ =0。將(32)式代入(27)式可 得梁挫屈前的統御方程式如下 3 1 2 4 1 ( ) 2 1 ) ( p p ,s p Ip ,s A C I A E EI EI C M = − ε θ + − − θ (33) 2.7 次要平衡路徑之統御方程式 若已經求得梁在主要平衡路徑的變形，則每一元素之θ10(s)及ε0為已知。若在 該平衡位置上給梁元素一x 及2 x3軸方向的擾動位移 v 、 w ，由(22)式、(23)式分別對 s 微分及(28)式、(29)式分別對 s 微分兩次的關係，可以得到以下之方程組
(
2 +)
=0 + + − − ζζζ ζζζζ ζζ ζζζζ αµ κζ κ γ αη , , , y , w v v I l w l(
2 −)
=0 − − + ζζ ζζζζ ζζζ ζζζζ µ κζ κ γ η , , , z , v w w I l v l (34) 其中 2 1 − = l s ζ ， (−0.5≤ζ ≤0.5) 0 1 ε η= +，
0 1 1) ,s (α θ κ = −，
γ θ0 β 1 s, =，
µ =εpA(35) 本節中v,w,θ2及θ3皆代表擾動量。假設(34)式之解可表示成下式
n n n n B A v w ζ ∑ = ∞ =0
(36) 將(36)式代入(34)式，可歸納出下列關係
Ntq v v= ， Nt q w w= (37) 本文中稱N ,v Nw 為形狀函數。其中q為未定係數組成的向量，可以由該元素兩端的擾 動位移及其微分決定，由(7)、(8)式與(10)、(37)式之關係可得 q T u= uq (38) 其中
{
v1 w1 φ21 φ31 v2 w2 φ22 φ32}
= u{
A0 B0 A1 B1 A2 B2 A3 B3}
= q(39) 2.8梁元素在擾動後的節點內力 將(28)、(29)式代入(15)式以及與(7)、(8)、(37)之關係，則可將梁元素節點上的內力 表示成 q N F= F (40) 其中
{
F21 F31 M21 M31 F22 F32 M22 M32}
= F 利用元素座標系統與總體座標系統之間的轉換關係矩陣 ge T 與(38)、(40)式之關係，梁元 素在總體座標上的節點內力向量可以表示成Fg =TgeNFTuq−1Tget ug (41) 其中ug元素在總體座標的節點擾動位移向量。 2.9 挫屈負荷 將相鄰元素在共同節點上之位移一致關係代入(41)式，將(41)式代入受力端之邊界條 件與元素在共同節點上之合力(矩)相等的關係式，則可以表示成下式 0 U K = (42) 其中N代表元素的數目，K為系統的剛度矩陣，U為系統的位移向量。 當U有非零解時，表示在原平衡位置附近還有其他的平衡位置，也就是發生挫屈。
因為(42)式為一齊項式(homogeneous equation)所以僅有當剛度矩陣K的行列式值det K
為零時才有非零的解。 3. 數值例題與結果 本章將探討在一定值之軸向保守力P及保守扭矩M 作用下，梁的側向挫屈。 如圖.4 所示，考慮一矩形斷面的懸臂梁，梁之長度為L = 480 mm，本文中所使用的 材料之楊氏模數E=1.931×105N mm2，剪力模數G=7.45×104N mm2的不銹鋼
。
其 A 端為固定端，B端為受力端。其受力端為自由端(如圖.4)，AB兩端都允許自由翹 曲(warping free)。扭矩的施加方式為 QT-1 型扭矩（如圖.4），其中M =2rF，r
為圓盤之半徑，F為保守力之大小，上標g代表總體座標系統。本文在進行挫屈分 析時，所施加的軸向負荷分別為P=0.0,0.4,0.6,0.8及0.9 Pcr，其中Pcr為梁僅受軸向壓 力時之挫屈負荷。在梁受力端為自由端時 2 2 4L EI Pcr =π z 。Mcr為梁在軸向保守扭矩作用 下 之 線 性 挫 屈 負 荷 的 經 典 解 。 在 梁 受 力 端 為 自 由 端 中 受 QT-1 型 扭 矩 作 用 下 之 z y cr I I L E M 2 π = 。 圖. 4 梁承受端點扭矩及軸力之結構圖 令Mnl為利用本文所提之分析方法而求得之挫屈負荷。由圖.5 中得知當軸向力增加 時，分析時用較少的元素即可得到相同的精度，在分析時為了確保精度皆取軸向力為零時所需的梁元素數目。本文中發現當長寬比(b/a)愈大時，分析時所需的元素數目亦需較 多。 No. of elements 30 50 100 150 200 0.0 2.0 4.0 0.4 0.0 0.9 0.8 0.6 = M cr /M nl (%) P/Pcr 圖. 5 梁受QT-1形扭矩及軸力作用下之元素收斂分析(CASE 6). 3 X a b X2 1 15 30 b/a QT-1 0 50 100 P/P 0.4 0.6 0.9 0.8 0.0 cr /M M nl cr (% ) 圖. 6 梁受QT-1型扭矩及軸力作用時之挫屈負荷圖 圖.6 為在不同之軸向保守力作用的情形下，不同長寬比(b/a)之斷面對應其受扭矩作 用時挫屈負荷在其經典解之百分比值之圖。由圖中可知，當軸向力為零時，扭矩作用在 梁之斷面的長寬比為一時(b/a = 1/1)，可得到幾乎與經典解有一致的結果，這是因(23)-(33) 式中許多耦合項會自動消失。而軸向力增加時，造成梁側向挫屈的扭矩負荷就變小。當
梁斷面的長寬比相差愈大時，則其挫屈負荷與經典解之差異也愈大。而軸向作用力的施 加，對長寬比(b/a)愈大之梁所造成之挫屈負荷影響愈小。 4 結論 本文中利用二階梁理論探討梁在軸力及扭矩作用下之側向-扭轉挫屈。本文中提出一 個以擾動位移的級數解求取挫屈負荷的方法，即先求出滿足主要平衡路徑之統御方程式 的解析解及滿足次要平衡路徑之統御方程式的側向擾動位移之級數解，再由相鄰元素在 共同節點上需滿足變形一致和合力一致的條件，及結構兩端點的位移和合力需滿足的邊 界條件，可以求得一齊項式之聯立方程式。使該聯立方程式之剛度矩陣的行列式值等於 零時的扭矩即為挫屈負荷。由數值結果我們可以觀察到以下的結論：當軸向力為零時， 扭矩作用在梁之斷面的長寬比為一時(b/a = 1/1)，可得到幾乎與經典解有一致的結果。當 梁斷面的長寬比(b/a)相差愈大時，則其挫屈負荷與經典解之比值愈小。軸向力增加時， 梁側向挫屈的扭矩負荷就變小，且長寬比(b/a)愈大其對挫屈負荷之影響愈小。 5 參考文獻
[1] M.M. Attard, Lateral buckling analysis of beams by the FEM, Comput. Struct. 23 (1985) 217-231.
[2] Y.L. Pi and N.S. Trahair, Prebuckling deflections and lateral buckling: part I- theory, J. Struct. Engrg. ASCE 118(11) (1992) 2949-2966.
[3] Y.L. Pi and N.S. Trahair, Prebuckling deflections and lateral buckling: part II- application, J. Struct. Engrg. ASCE 118(11) (1992) 2967-2985.
[4] K.M. Hsiao, Corotational total Lagrangian formulation for three- dimensional beam element, AIAA J. 30(3) (1992) 797-804.
[5] Y. Goto, X. S. Li and T. Kasugal, Buckling analysis of elastic space rods under torsional moment, J. of Engrg. Mech. ASCE 122 (1996) 826-833.
[6] K.M. Hsiao, R.T. Yang and W.Y. Lin, A consistent finite element formulation for linear buckling analysis of spatial beams, Comput. Methods Appl. Mech. Engrg. 156 (1998) 259-276.
[7] J.C. Simo and L. Vu-Quoc, The role of non-linear theories in transient dynamic analysis of flexible structures, J. Sound Vib. 119 (1987) 487-508.
[8] J.M. Anderson and N.S. Trahair, Stability of monosymmetric beams and cantilevers, J. Struct. Div. ASCE 98(1) (1972) 269-286.
[9] N.S. Trahair and S.T. Woolcock , Effect of major axis curvature on I-beam stability, J. Engrg. Mech. Div. ASCE 99(EM1) (1974) 85-98.
[10] 蕭國模, “ 二階梁理論及其在側向-扭轉挫屈的應用（I）” 國科會專題研究計畫成
果報告, NSC91-2212-E009-041, 新竹 (2003)
[11] H. Goldstein, Classical Mechanics, Addision-Wesley, Reading, MA, (1980).
[12] T.J. Chung, Continuum Mechanics, Prentice-Hall N.J., (1988).
[13] D.J. Dawe, Matrix and Finite Element Displacement Analysis of Structures, Oxford
