開口薄壁梁之自由振動及幾何非線性動態分析(I)

(1)

行政院國家科學委員會專題研究計畫期中進度報告

開口薄壁梁之自由振動及幾何非線性動態分析(1/3)

計畫類別：個別型計畫計畫編號： NSC94-2211-E-009-027- 執行期間： 94 年 08 月 01 日至 95 年 10 月 31 日執行單位：國立交通大學機械工程學系(所) 計畫主持人：蕭國模報告類型：精簡報告報告附件：出席國際會議研究心得報告及發表論文處理方式：本計畫可公開查詢

中華民國 95 年 7 月 6 日

(2)

開口薄壁梁之自由振動及幾何非線性動態分析

(1/3)

On free vibration analysis and geometrically nonlinear dynamic analysis of

thin-walled open-section beams (1/3)

計畫編號：

94－2211－E－009－027－

執行期限：2005 年 08 月 1 日至 2006 年 05 月 31 日

主持人：蕭國模國立交通大學機械工程學系

計畫參與人員：陳弘虎、劉宗帆

中文摘要：本研究為三年期計畫的第一年，主要目的是以一致性共旋轉法推導一個二節點十四個自由度的薄壁開口梁元素，並建立開口薄壁梁的幾何非線性運動方程式及振動方程式。本文中推導的梁元素有兩個節點，每個節點有七個自由度。本文中將元素節點定在斷面剪心，並取剪心軸當作梁元素變形的參考軸。本研究在當前梁元素變形的位置上建立兩個重合的元素座標，一個為固定在元素當前的變形位置之固定元素座標，一個為與元素一起剛體運動但不一起變形的移動元素座標，本研究在當前的固定元素座標定義元素的節點位移向量、旋轉向量、虛位移向量、虛旋轉向量、速度、加速度、角速度、角加速度、節點力、剛度矩陣、質量矩陣，但在移動元素座標上描述元素的變形及定義節點變形參數。移動元素座標的速度、加速度、角速度、角加速度是由元素當前的節點速度、加速度、角速度、角加速度決定，元素節點變形參數在移動元素座標的虛擾動量是由元素節點之虛位移向量、虛旋轉向量及當前的節點位移向量、旋轉向量、變形參數的決定。本研究利用虛功原理和 D’Alembert 原理，以及完整的幾何非線性梁理論的一致性二次線性化在固定元素座標元素座標上推導元素節點變形力及慣性力。本研究中保留了變形力中撓曲、扭曲及軸向變形間之耦合項和慣性力中速度間的耦合項。軸向扭轉率的三階項是所有三階項中的支配項，而且是反映梁受到純扭矩時產生非線性行為重要的內力項，因此在元素節點內力中必須予以考慮。本研究推導的元素節點慣性力時不考慮變形與速度及加速度間的耦合項，因其在元素增加時會趨近於零。本研究在推導元素的節點變形力時，保留了當前的元素節點位移向量、旋轉向量，所以由元素節點變形力對節點參數微分可求得切線剛度矩陣。元素的一致性質量矩陣 (consistent mass matrix)可由節點慣性力對節點加速度微分求得。本研究中建立開口薄壁梁的幾何非線性運動方程式

(3)

Abstract.

This study is the first year part of a three-year project. A two-node displacement-based thin-walled open-section beam element with seven degrees of freedom per node is developed by using consistent corotational formulation. The vibration equation and the geometrically nonlinear dynamic equation are constructed for thin-walled beam with open section.

In this study, element nodes are chosen to be the shear centers of end sections of the element. The shear center axis is employed as the reference axis of the beam element. Two coincident element coordinate systems are constructed at the current configuration of the beam element. One is fixed element coordinates and the other is element moving coordinates. The element nodal displacement vector, rotation vector, virtual nodal displacement vector, virtual rotation vector, nodal velocity, nodal acceleration, nodal angular velocity, nodal angular acceleration, internal nodal force, stiffness matrices and mass matrix are defined in the current fixed element coordinates. The element deformations and nodal rotation parameters are referred to the initial undeformed geometry of the beam element and described in the current moving element coordinates. The velocity, acceleration, angular velocity and angular acceleration of the moving element coordinates are determined by the current element nodal velocity, nodal acceleration, nodal angular velocity and nodal angular acceleration. Both the element deformation nodal forces and inertia nodal forces are systematically derived by consistent linearization of the fully geometrically nonlinear beam theory using the D’Alembert principle and the virtual work principle. All the deformational coupling effect and velocity coupling terms are retained. The third-order term of twist rate is the dominant term for the third-order terms and may be a very important term to reflect the nonlinear behavior of the beam subjected to a pure torque. Hence, the third-order term of twist rate is also considered in the nodal forces. The terms relevant to the twist angle, slopes and the length of the beam element, which will converge to zero with the decrease of element size, are removed from the element nodal forces and the element matrices.

Due to the way used to define the element coordinates, except the current axial displacement at element node 2, the value of the current element nodal displacements is zero. The value of the element nodal rotation vector is reset to zero at the current configuration of the beam element. However, all the current element nodal displacement vectors and rotation vectors are retained in the derivation of the element nodal force. Thus, the element tangent stiffness matrix may be obtained by differentiating the element nodal force vector with respect to the element nodal parameters. The element consistent mass matrix may be obtained by differentiating the element nodal inertial force vector with respect to the second time derivative of the nodal parameters. In this study, the geometrically nonlinear equations of motion of thin-walled beam with open section are constructed.

Keywords: Thin-walled beam, Geometrical nonlinearity, Co-rotational formulation, Finite element method, Equations of motion

(4)

1 導論 梁在結構工程系統中，長久以來一直扮演著非常重要的角色，在機械、航空太空、建築、車輛及土木工程中皆有很廣泛的應用。有些結構如飛機、太空船、船舶等，為了減輕重量而使用了高強度材料及薄壁斷面。開口薄壁梁結構運動時的耦合效應包括變形間的耦合、慣性間的耦合、及慣性與變形間的耦合，變形間的耦合效應包括撓曲、扭曲、扭轉翹曲及軸向變形間的耦合，慣性間的耦合效應包括撓曲、扭曲、扭轉翹曲及軸向速度間的耦合，慣性與變形間的耦合包括撓曲、扭曲、扭轉翹曲及軸向速度與變形間的耦合。為了正確地考慮開口薄壁梁各種變形及慣性間的耦合效應，在推導其運動方程式時必須使用開口薄壁梁正確的變形機制，其中包括了正確地考慮梁斷面的有限旋轉及有限旋轉的非向量特性，開口薄壁梁的運動方程式為一高度的非線性方程式，但該非線性為幾何非線性，幾何非線性主要是由剛體旋轉造成，若將剛體旋轉從總位移中去掉，則剩下的位移仍為小變形及小旋轉。為了去掉運動方程式的高次非線性項但仍保持其精度，所以本研究擬採用共旋轉法[1-4]除掉剛體旋轉，即將梁分割成數個元素，然後在每一梁元素當前的變形位置上建立一元素座標。每一梁元素的變形、節點內力、運動方程式都是建立在該元素座標上。文獻[4]中用 corotational formulation 推導一個含翹曲剛度之不對稱薄壁梁元素，[2-4]中提出一開口薄壁梁正確的變形機制及建立一套描述梁斷面的有限旋轉的方法以決定變形後斷面座標之方位，該方法能夠正確地考慮撓曲、扭曲與軸向變形的耦合效應，[2-4]中探討薄壁梁的幾何非線性靜態反應及挫屈負荷，其結果相當精確，但文獻[2-4]沒有考慮動態效應，文獻[1]採用 corotational formulation 及完整之幾何非線性梁理論，利用虛功原理、d’Alembert 原理推導梁元素的運動方程式，並探討三維梁的幾何非線性動態反應，由[1]之例題可發現其結果相當精確，文獻[1]的方法能夠正確地考慮撓曲、扭曲與軸向變形及慣性間之耦合的效應，但文獻[1]沒有考慮梁斷面的翹曲剛度及剪力中心與形心不一致的問題，故本研究擬採用文獻[2-4]中提出的薄壁梁變形機制及以文獻[1] 推導梁元素的運動方程式的方法為基礎，並將其做必要的修改，推導薄壁梁元素的非線性運動方程式。為考慮各種變形及慣性間完整的非線性耦合，本研究將利用虛功原理、d’Alembert 原理與完整的非線性梁理論之一致性二階線性化來推導梁元素的節點變形內力及節點慣性內力。文獻[1-7]用虛功原理推導元素節點內力時，在梁元素當前的元素座標上定義元素節點旋轉參數的擾動量，在推導上很方便，這種方法可稱為廣義的共旋轉推導法，本研究將採用一致性的共旋轉推導法，即在當前梁元素變形的位置上建立兩個重合的元素座標，一個為固定在元素當前的變形位置之固定元素座標，一個為與元素一起剛體運動但不一起變形的移動元素座標，本研究在當前的固定元素座標定義元素的節點位移向量、旋轉向量、虛位移向量、虛旋轉向量、速度、加速度、角速度、角加速度、節點力、剛度矩陣、質量矩陣，但在移動元素座標上描述元素的變形及定義節點變形參數。移動元素座標的速度、加速度、角速度、角加速度是由元素當前的節點速度、加速度、角速度、角加速度決定，元素節點變形參數在移動元素座標的虛擾動量是由元素節點之虛位移向量、虛旋轉向量及當前的節點位移向量、旋轉向量、變形參數的決定。本研究在推導元素的節點變形力時，保留了當前的元素節點位移向量、旋轉向量，所以由元素節點變形力對節點參數微分可求得切線剛度矩陣。元素的一致性質量矩陣 (consistent mass matrix) 可由節點慣性力對節點加速度微分求得。梁在空間中高速運動時，速度間的耦合效應不能忽略，所以本研究將保留全部速度間的耦合項，本研究將忽略變形和慣性間的耦合效應趨近於零的項，故慣性力中僅部份保留至節點參數的一次項。本研究推導開口薄壁梁元素時，為了保留各種變形間的耦合，須先將節點變形參數保留到適當的階數，然後在適當的時候將節點變形力及慣性力作一致性二階線性化[1-9]，以求得正確之節點變形

(5)

力、節點慣性力、切線剛度矩陣及質量矩陣。本研究用一致性一次線性化去掉梁的非線性運動方程式之非線性項求得其線性振動方程式。 2 理論推導 本研究梁元素之基本假設、座標系統、梁的變形機制及其描述方法和文獻[1-4]大致一樣，本文採用一致性共旋轉法推導梁元素節點內力。 2.1 基本假設 本文對非線性梁元素的推導，作如下的假設： (1) 梁為細長的等斷面梁，且Euler-Bernoulli 假說成立。 (2) 梁斷面為不對稱斷面。 (3) 當去除剪心扭轉造成正應變，梁元素的形心軸剩餘之縱向正應變(longitudinal normal strian)為一常數。 (4) 梁元素在斷面上沒有變形，且斷面內的應變可以忽略。 (5) 梁元素斷面的翹曲為梁元素的軸向扭轉率與該梁的聖維南(Saint Venant)翹曲函數的乘積。 (6) 梁元素的變形為小變形。 2.2 座標系統

本文採用共旋轉之全拉格蘭日推導法(co-rotational total Lagrangian formulation)。為了描述系統的運動、元素的變形、邊界條件、以及與結構變形位置相關的外力 (configuration dependent load)，本文中共定義了四套直角座標系統：

(1) 固定總體座標系統(圖.1)，X_iG(i=1 ,2,3)；系統的節點座標及方位，其他座標系統之座標軸的方向餘弦，皆在此座標系統中定義。 (2) 元素座標系統(圖.1，圖.2)，x_i(i=1 ,2,3)、x_i(i=1 ,2,3)；此座標系統建立在當前梁元素變形的位置上，本研究採用兩個元素座標，一個為固定在元素當前的變形位置之固定元素座標xi，一個為與元素一起剛體運動，但不一起變形的移動元素座標xi。兩個元素座標在當前梁元素變形的位置上是重合的。此座標系統是由[7] 中提出的方法決定。本研究在當前的固定元素座標定義元素的節點位移向量、旋轉向量、虛位移向量、虛旋轉向量、速度、加速度、角速度、角加速度、節點力、剛度矩陣、質量矩陣，但在移動元素座標上描述元素的變形及定義節點變形參數。移動元素座標的速度、加速度、角速度、角加速度是由元素當前的節點速度、加速度、角速度、角加速度決定，元素節點變形參數在移動元素座標的虛擾動量是由元素節點之虛位移向量、虛旋轉向量及當前的節點位移向量、旋轉向量、變形參數的決定。 (3) 元素斷面座標系統(圖.1)，x_iS(i=1 ,2,3)；此座標系統與元素斷面一起剛體運動，其原點剛接於未翹曲斷面的形心上，x₁S軸為未翹曲斷面的法線方向，x₂S 軸與x₃S 軸分別與未翹曲斷面的主軸重合。元素的變形是由斷面座標相對於元素座標的旋轉來決定。 (4) 負荷基底座標系統，X_iP(i=1 ,2,3)；此座標系統是用來描述與結構變形位置相關的作用力機制。

(6)

2.3 旋轉向量及其時間微分 本文中使用旋轉向量來表示一個有限旋轉，一向量b受到一旋轉向量a=φe的作用而轉到一個新的位置b′，向量b′與b之間的關係可表示成[10] ) ( sin ) )( cos 1 ( cos e b e e b b b′= φ+ − φ ⋅ + φ × b Ι a a I a I sin ( ) 1 cos ( )] [ 2 × × − + × + = φ φ φ φ =Rb (1) 其中I 為3×3的單位矩陣，×表向量外積，．表向量內積，φ表逆時鐘方向旋轉角，e表旋轉軸的單位向量，R 為旋轉矩陣。 由(1)式 b′對時間的微分可表示成 b R R Rb b′ ₌ ₌ _t _′ dt d _& (2) 因 b′的長度固定，所以其對時間的微分可表示成 b ω b _′ × = ′ dt d (3) 其中ω為角速度向量。因R&Rt為一反對稱矩陣，從(2)式和(3)式可得 _ω_×_I₌_R_&_Rt₍₄₎ 由(1)及(4)式可得ω和旋轉向量 a&對時間的一次微分有以下關係[11]： ω=Γ(a)a_& )] ( ) ( [ ) (a = I+ ₁ a×I + ₁a× a×Ι Γ a b (5) 其中 2 1 1 cos φ φ − = a , 1 (1 sin ) 2 1 _φφ φ − = b 。當旋轉向量a有一微小變量δ 時會使向量 a b′繞x_i (i=1,2,3)軸做微小旋轉 i δϕ ,δ 與a δϕ={δϕ₁δϕ₂δϕ₃}的關係和_a&與ω的關係相同[11]，即 δϕ=Γ(a)δa (6) 將(5)式對時間微分可得 ω_& =Γ&(a)a_&+Γ(a)a_&&

) ( ) ( ) ( ) (a =a₁ a×I + ₁a× a×I + ₁ a×I Γ& & b& a &

) ( ) ( ₁ 1a× a×I + a× a×I +b & b & (7) 其中_ω&為角加速度,a&&為旋轉向量對時間的兩次微分。由(6)式可知當a=0時, Γ(0)=I所以 ω =a& (8) ω_{& =}a_&& (9)

(7)

即在旋轉向量a=0時，其對時間的一次，二次微分之值，等於角速度、角加速度，本研究中，在任何時刻及位置都將節點的旋轉向量重新設定為零，所以旋轉向量對時間的一次、二次微分之值與角速度、角加速度之值相同。 2.4 元素座標與元素斷面座標之關係 令e_i與e_iS_{(i = 1, 2, 3) 分別代表固定元素座標 x} i方向的單位向量與元素斷面座標 x_iS軸方向的單位向量。由座標系統的定義可知，在變形前x 軸與_i x 軸是一致的，而且_iS 變形後e 與(10)式的 t 重合。在本文中假設變形後的單位向量₁S e_iS_{(i = 1, 2, 3)，是由以下} 兩個旋轉向量連續作用於單位向量e (i = 1, 2, 3)來決定 [7]: _i n θ_n =θ_n _， θ_t =θ₁t, (10) 其中 } , , 0 { } ) ( , ) ( , 0 { 3 2 2 1 2 3 2 2 2 2 1 2 3 2 2 2 n n = + + = θ θ θ θ θ θ n ₍₁₁₎ t={cosθ_n,θ₂,θ₃} (12) cosθ_n = −(1 θ₂2 −θ₃2 1 2) (13) ds x dw )( 2=− θ , ds x dv )( 3= θ (14,15) n 為垂直 e₁ 和 t 的單位向量。θ_n為cosθ_n的反函數。 v s( ) 與 w s( ) 為剪心軸的側向位 移。 s 為剪心軸在變形後的弧長。 由 (1)、(10)-(15)式可得e_i 與 e_iS_{(i = 1, 2, 3)的關係可表示如下} e_iS t R R e Re i i =[ , ₁, ₂] = (16) 其中R稱為旋轉矩陣。因R為θi (i = 1, 2, 3)的函數，所以本文中稱θi為旋轉參數。當θ_i(i=1,2,3)分別有一微小變化δθ_i時，斷面座標會旋轉到一個新的位置，此一新的位置可由元素座標繞x (_i i = 1, 2, 3)軸分別作微小旋轉δϕ_i(i=1,2,3)而得[7] ϕ ϕ = θ δ δ θ θ θ θ δ 1 2 3 2 3 1 0 0 1 2 2 1 − = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − − T (17) 2.5 梁元素之位移、應變與絕對加速度 圖.1中Q點為梁元素中的任意點，P點為Q點在剪心軸上的對應點，即P點與Q點位於同一斷面上。在固定元素座標上，Q點在梁元素變形前後的位置向量可分別表示如下 3 2 1 0 e ( )e ( )e r =x + y−y_p + z−z_p ₍₁₈₎ S x S p S p p x t v x t w x t y y z z x t 1 , 1 3 2 3 2 1 ( , ) ( , ) ( ) ( ) ) , ( ) ( e e e e e e r ω θ + − + − + + + = (19)

(8)

其中x_p( tx, )、v( tx, )、以及w( tx, )分別是P點在x_i( i = 1, 2, 3)軸上的座標，y_p、z_p及 y 、 z分別是P和Q在 S x₂ 軸與x₃S軸的座標。式(19)中首三項是P點的位置向量，四、五項是Q 相對於P的位置向量，最末項表示Q點扭轉翹曲的位移。 x x _d d ₁ 1, θ θ = 是沿變形後的剪心軸的斷面之軸向扭轉率，ω =ω( zy, )代表等斷面梁的聖維南(Saint Venant)翹曲函數。本文中若( 表示定義於固定元素座標中的變數，則) (_)表示定義於移動元素座標中的對應變數，例如r 是移動元素座標x 上的位置向量。如圖.2所示，在時刻t，移動元素_i 座標x 與固定元素座標_i x 是重合的，所以_i r(t)的值與(19)式中r(t)的值相同，但其變分 並不相同。當元素節點 j ( j = 1, 2) 受到一定義於固定元素座標的擾動位移向量δ 及u_j 擾動旋轉向量δφ_j作用後，在 t+δt時刻，Q點位置向量可表示成 1 1 ) ) ( ( ) ( A r r u u r t+δt = _x_x t +δ + +δ (20) 其中A_x_x為受擾動後的移動元素座標與固定元素座標的轉換矩陣，u₁為元素節點1的位移向量，δu₁元素節點1的擾動位移向量。位置向量r 的擾動量δr可表示成 ) ( ) ( ) ) ( ( ) ( ) (t t r t A_x_x r t r u₁ t u₁ r t r r= +δ − = +δ + +δ − δ (21) 移動元素座標與固定元素座標的轉換矩陣移動A_x_x是由元素節點之虛位移向量、虛旋轉向量及當前的節點位移向量、旋轉向量、變形參數的決定。

本文中應變的度量是採用Green strain。本文中以ε_ij(i=1 ,2 ,3 ,j=1 ,2 ,3)表示Green strain。由於基本假設(4)，所以本文只考慮ε₁₁、ε₁₂與ε₁₃，因本文在移動元素座標上描述元素的變形及定義節點變形參數，所以應變表示如下 ₁₁ ₁t ₁ ₁₂ ₁t ₂ ₁₃ ₁t ₃ 2 1 , 2 1 , ) 1 ( 2 1 g g g g g g − = = = ε ε ε (22) z y x

∂

r g r g r g₁= ∂ , ₂= ∂ , ₃ = ∂ r δ 為r 對固定元素座標的擾動量，所以將(20)式δ 除以一微小時間r δ ，令t δt→0，可得絕對速度_{r& ，再將r& 對時間微分便可得絕對加速度r&&之完整表示式}

1 2

)

(Ω ΩΩ r Ωr r u

r_&&= & + _c + &+ &&+ _&& (23) 其中Ω& 與 Ω 是將擾動後元素座標x_i(i=1,2,3)的角加速度與角速度以矩陣形式表示。 2.6 元素顯節點參數的擾動與隱節點參數的擾動關係 本研究與文獻[4]一樣採用顯節點參數及隱節點參數，但本研究採用兩套顯節點參數，令δqϕ及δqφ為定義於固定元素座標中的第一及第二顯節點參數的擾動量並可表示成 } , , , , {δ ₁ δϕ₁ δ ₂ δϕ₂ δβ δq_ϕ = u u ₍₂₄₎ } , , , , {δ ₁ δφ₁ δ ₂ δφ₂ δβ δq_φ = u u ₍₂₅₎

(9)

其中δuj ={δuj,δvj,δwj}，δϕj ={δϕ1j,δϕ2j,δϕ3j}，δβ={δβ₁,δβ₂}， } , , { ₁_j ₂_j ₃_j j δφ δφ δφ δφ = ， j=1, 2，δuj的分量為節點 j 在當前變形位置的元素座標xi( i =1, 2, 3)軸方向的擾動位移，δϕj的分量為節點 j 在當前變形位置的元素座標繞xi( i = 1, 2, 3) 軸的擾動旋轉，δβj為扭轉率在節點 j 的擾動量，δφij( i = 1, 2, 3)為節點 j 的擾動旋轉向 量δφj在xi( i =1, 2, 3)軸方向的的分量。對應於δuij的廣義節點力 fij，為在xi軸方向的力；對應於δϕij的廣義節點力mij，為繞xi軸的傳統力矩。對應於δβj的廣義節點力Bj為廣義雙力矩，對應於δφij的廣義節點力m 為廣義力矩。φ_ij 依據(6)式，顯節點參數的擾動量δqϕ與旋轉向量的擾動量δq 有以下關係 φ φ ϕφ ϕ δ δq =(I+T ) q (26) 若梁元素在當前的變形位置受到一虛位移δq 的擾動，可以由文獻[4]的方法決定擾_φ 動後的移動元素座標及定義在該元素座標的節點旋轉參數θ_ij( i =1, 2, 3, j =1, 2)及元素 的弦長 l [4]，節點旋轉參數及元素的弦長的擾動量可表示成 ij ij ij θ θ θ δ = − ，δ_l= _l−_l (27) 本研究與文獻[4]採用一樣的隱節點參數，令δq 為定義於擾動後的移動元素座標的θ 隱節點參數的擾動量並可表示成 θ δq _{ _, _, _, *_, _} 2 2 * 1 1 δθ δ δθ δβ δu u = (28) 其中δuj ={δuj,δvj,δwj}，δθ*_j ={δθ₁*_j,δθ₂*_j,δθ₃*_j}={δθ₁_j,−δw′_j,δv′_j}， δβ={δβ₁,δβ₂} j=1, 2，δ 的分量為節點 j 在擾動後的移動元素座標uj xi( i =1, 2, 3)軸方向的擾動位 移，δθ1j定義於(27)式，δ 、w′j δ 分別為v′j x w ∂ ∂ 、 x v ∂ ∂ 在節點 j 的擾動量，δβ 為扭轉率j 在節點 j 的擾動量。對應於δ 的廣義節點力為uij fij，對應於δθij的廣義節點力為m ，ijθ 對應於δβ 的廣義節點力j Bj為廣義雙力矩。因為θ 在變形後不為零，所以其變分並不*_j 是繞x_i軸的無限小旋轉，所以廣義力矩m 並非繞θ_ij xi軸的傳統力矩。由(14)及(15)式節點參數的定義及(26)-(28)式，隱節點參數的擾動量δq 、第二顯節點參_θ 數的擾動量δq 及第一顯節點參數的擾動量_φ δq_ϕ有以下關係 ϕ ϕφ φ θ φ φ θ θ δ δ δq (T T ) q (T T )(I Tt ) q R R = + + + = (29) 2.7 梁元素之節點力 為了推導方便，本文將元素隱節點參數與顯節點參數分為四個廣義節點位移向量 i u 、u 與i uϕ_i ( i =a, b, c, d) } , {u₁ u₂ a = u ，u_b ={v₁,v₁′,v₂,v₂′} ，u_c ={w₁,w₁′,w₂,w₂′}，u_d ={θ₁₁,β₁,θ₁₂,β₂}

(10)

} , {u₁ u₂ a = ϕ u ，uϕ_b ={v₁,ϕ₃₁,v₂,ϕ₃₂}，uϕ_c ={w₁,ϕ₂₁,w₂,ϕ₂₂}， } , , , {ϕ₁₁ β₁ ϕ₁₂ β₂ ϕ ₌ d u (30) 應用D’Alembert 原理及虛功原理，可知 d t d c t c b t b a t a f u f u f u f uϕ δ ϕ δ ϕ δ ϕ δ + + +

∫

+ + + = V t V(σ11δε11 2σ12δε12 2σ13δε13)dV (ρδr &&r)dV (31) I d t d I c t c I b t b I a t a D d t d D c t c D b t b D a t af u f u f u f u f u f u f u f u δ δ δ δ ϕ δ ϕ δ ϕ δ ϕ δ + + + + + + + = 其中u_i( i =a,b,c,d)的定義方式與(30)式之u 相同，表示定義在_i x 座標的節點位移向_i 量。f_i =f_iD+f_iI，f 與_iD f_iI( i =a,b,c,d)與為對應於δu 的顯節點變形力與慣性力，ϕ_i D i f ( i =a,b,c,d)與為對應於δu_i( i =a,b,c,d)的隱節點變形力。V 為未變形梁元素的體積， ρ為梁元素之密度， ₂ 2 dt d r r_&&= 為絕對加速度。δε₁_j( j =1, 2, 3)是(21)式中ε₁_j的變分，σ₁_j( j =1, 2, 3)是對應於ε₁_j的第二類 Piola-Kirchhoff 應力。本文中假設材料為線彈性材料，其應力－應變關係為 , 2 , 2 , ₁₂ ₁₂ ₁₃ ₁₃ 11 11 ε σ ε σ ε σ =E = G = G ₍₃₂₎

其中 E 是楊氏係數(Young's modulus)，G是剪力係數(Shear modulus)。

將(21)、(23)及(32)式代入(31)式中，並保留元素隱節點參數至二次項，由(31)式的等號兩邊的δu_i (i =a,b,c,d)之係數對等之關係，則可得隱節點變形力f_iD (i=a,b,c,d)與顯節點慣性力f_iI (i=a,b,c,d)。利用(28)式與反梯度法則(contagradient law)[12]，顯節點變形力向量可以表示成 D t R t D θ φ θ ϕφ T T f T I f =( + )( + ) (33) 其中 } , , , , {f₁ m₁ f₂ m₂ B fD = D D D D ， { ₁ , ₂ , ₃D} j D j D j D j = f f f f ， { ₁ , ₂ , ₃D} j D j D j D j = m m m m } , , , , {f₁ m₁ f₂ m₂ B f_θD = θ θ θ θ ，fθ_j ={f₁θ_j, f₂θ_j, f₃θ_j}，mθ_j ={m₁θ_j,m₂θ_j,mθ₃_j}( j = 1, 2) ，B={B1,B2}₍₃₄₎ 2.8 元素剛度矩陣及質量矩陣 對應於元素顯節點參數的元素切線剛度矩陣k ，稱為顯切線剛度矩陣。依切線剛度矩 陣的定義可知df =kdq_ϕ，所以由(33)式的元素節點變形力向量f 對顯節點參數D q_ϕ微分求得顯切線剛度矩陣。依切線剛度矩陣的定義可得 ϕ ϕ ϕ _qf q q k f d d d D ∂ ∂ = = ， ϕ q f k ∂ ∂ = D (35) 元素的質量矩陣可由節點慣性力f_iI (i=a,b,c,d)式對節點加速度微分而得

(11)

_ϕ j I i ij u f M && ∂ ∂ = ₍i_,j=a_,b_,c_,d₎ (36) 2.10 系統運動方程式 系統離散化後的動平衡方程式, 是由在元素座標上算出的元素節點變形力、慣性力經標準的座標轉換後, 在總體座標上所合成之系統內力向量及系統外力向量所組成, 可表示如下： _ϕ=FI +FD −P=0 ₍₃₈₎ 其中ϕ 是不平衡力,F 是系統慣性力內力向量I ,_{F 是系統變形力內力向量}D _{, P 表系統} 外力向量。 2.9 梁的振動程式 本研究中開口薄壁梁的振動是以穩態為平衡點的振動，一致性一次線性化去掉梁的非線性運動方程式之非線性項求得其線性振動方程式，振動方程式可表示成 0 Θ M K − ) = ( G _λ2 G ₍₃₇₎ 其中λ為自然頻率，Θ為對應於λ振動模態，系統的切線剛度矩陣K 和質量矩陣G G M 是由梁元素的切線剛度矩陣 k 與質量矩陣 M 經座標轉換後對應而成。 3. 結論 本研究採用一致性共旋轉法推導一個二節點十四個自由度的薄壁開口梁元素，並建立開口薄壁梁的幾何非線性運動方程式及振動方程式。本文中推導的梁元素有兩個節點，每個節點有七個自由度。本文中將元素節點定在斷面剪心，並取剪心軸當作梁元素變形的參考軸。本研究在當前梁元素變形的位置上建立兩個重合的元素座標，一個為固定在元素當前的變形位置之固定元素座標，一個為與元素一起剛體運動但不一起變形的移動元素座標，本研究在當前的固定元素座標定義元素的節點位移向量、旋轉向量、虛位移向量、虛旋轉向量、速度、加速度、角速度、角加速度、節點力、剛度矩陣、質量矩陣，但在移動元素座標上描述元素的變形及定義節點變形參數。移動元素座標的速度、加速度、角速度、角加速度是由元素當前的節點速度、加速度、角速度、角加速度決定，元素節點變形參數在移動元素座標的虛擾動量是由元素節點之虛位移向量、虛旋轉向量及當前的節點位移向量、旋轉向量、變形參數的決定。本研究利用虛功原理和 D’Alembert 原理，以及完整的幾何非線性梁理論的一致性二次線性化在固定元素座標元素座標上推導元素節點變形力及慣性力。本研究中保留了變形力中撓曲、扭曲及軸向變形間之耦合項和慣性力中速度間的耦合項。軸向扭轉率的三階項是所有三階項中的支配項，而且是反映梁受到純扭矩時產生非線性行為重要的內力項，因此在元素節點內力中必須予以考慮。本研究推導的元素節點慣性力不考慮變形與速度及加速度間的耦合項，因其在元素增加時會趨近於零。本研究在推導元素的節點變形力時，保留了當前的元素節點位移向量、旋轉向量，所以由元素節點變形力對節點參數微分可求得切線剛度矩陣。元素的一致性質量矩陣 (consistent mass matrix)可由節點慣性力對節點加速度微分求得。本研究用一致性一次線性化去掉梁的非線性運動方程式之非線性項求得其線性振動方程式。

REFERENCES

(12)

Formulation for Geometrically Nonlinear Dynamic Analysis of 3-D Beams, Computer Methods in Applied Mechanics and Engineering 169 (1999) 1-18.

[2] K.M. Hsiao, and Wen Yi Lin, A Co-Rotational Formulation for Thin-Walled Beams with Monosymmetric Open Section, Computer Methods in Applied Mechanics and Engineering 190 (2000) 1163-1185.

[3] W. Y. Lin and K.M. Hsiao, Co-Rotational Formulation for Geometric Nonlinear analysis of Doubly Symmetric Thin-Walled Beams, Computer Methods in Applied Mechanics and Engineering 190 (2001) 6023-6052.

[4] W. Y. Lin and K.M. Hsiao, Corotational Finite Element for Thin-Walled Beams with Generic Open Cross section, Proceedings of The Fifth International Conference on Space Structures, Guildford, UK 19-21 August, 2002.

[5] K.M. Hsiao and W.Y. Lin, A co-rotational finite element formulation for buckling and postbuckling analysis of spatial beams, Computer Methods in Applied Mechanics and Engineering 188 (2000) 567-594.

[6] W. Y. Lin and K.M. Hsiao, A Buckling and Postbuckling Analysis of rods Under End Torque and Compressive Load, Computer Modeling in Engineering & Sciences 4 (2003) 259-271.

[7] K.M. Hsiao, Corotational total Lagrangian formulation for three- dimensional beam element, AIAA J. 30(3) (1992) 797-804.

[8] 蕭國模，91 年"二階梁理論及其在側向-扭轉挫屈的應用(I)"國科會研究報告， NSC91－2211－E009－041

[9] 蕭國模，92 年"二階梁理論及其在側向-扭轉挫屈的應用(1/2)"國科會研究報告， NSC91－2211－E009－046.

[10] H. Goldstein, Classical Mechanics, Addision-Wesley Publishing Company, 1980. [11] M. Borri, F. Mello and S.N. Atluri, Variational Approaches for Dynamics and

Time-Finite-Element: Numerical Studies, Computational Mechamics 7 49-76 (1990). [12] D.J. Dawe, Matrix and Finite Element Displacement Analysis of Structures, Oxford