羅　帥　劉　偉

(1.紹興文理學院　土木工程學院,浙江　紹興312000；2.廣東交通職業(yè)技術(shù)學院,廣東　廣州510650)

?

框架結(jié)構(gòu)有限元模型的坐標變換方法

羅帥1劉偉2

(1.紹興文理學院土木工程學院,浙江紹興312000；2.廣東交通職業(yè)技術(shù)學院,廣東廣州510650)

針對框架結(jié)構(gòu)有限元分析過程中的坐標轉(zhuǎn)換關(guān)系問題，基于三角函數(shù)和應變能關(guān)系推導了二維坐標系中框架單元局部剛度矩陣和整體剛度矩陣之間的轉(zhuǎn)換關(guān)系，利用功能原理推導荷載向量的局部坐標和整體坐標轉(zhuǎn)換方法，利用所推導的轉(zhuǎn)換矩陣，根據(jù)已知的框架結(jié)構(gòu)局部剛度矩陣變換得到整體剛度矩陣.與常用的基于幾何投影關(guān)系推導出來的單元坐標轉(zhuǎn)換公式相比，新方法具有推導過程邏輯嚴謹，物理意義明確等優(yōu)點，算例分析結(jié)果表明該方法適合編程應用.

框架結(jié)構(gòu)；有限元模型；坐標變換矩陣；局部剛度矩陣；整體剛度矩陣

0　引言

框架結(jié)構(gòu)在我國是應用較多的結(jié)構(gòu)形式.框架結(jié)構(gòu)又稱構(gòu)架式結(jié)構(gòu)，是由梁、板、柱構(gòu)成的承重體系：樓板搭在梁上，梁支撐在兩邊的柱子上，由此把荷載傳遞給柱子，最后沿著柱子在高度方向傳給結(jié)構(gòu)基礎.房屋的框架按所用材料分為鋼筋混凝土框架、木框架、鋼框架等，其中最常見的是鋼筋混凝土框架和鋼框架，鋼筋混凝土框架主要應用于學校、住宅、辦公樓等；鋼框架則多用于機場、火車站、停車場、體育館等一些有特殊用途的建筑結(jié)構(gòu)中.

目前框架結(jié)構(gòu)的力學分析主要利用有限元方法進行[1]，有限元方法分析框架結(jié)構(gòu)的基本原理是利用局部坐標系下各種規(guī)則的基本單元來模擬實體結(jié)構(gòu)中遇到的各種復雜幾何形狀[2].進行二維或三維結(jié)構(gòu)分析過程中，在局部坐標系下建立單元的剛度矩陣后，隨后的重要步驟就是通過局部坐標與整體坐標之間的映射關(guān)系組裝整體坐標系下的結(jié)構(gòu)剛度矩陣[3].已有的文獻中通常都是直接給出結(jié)構(gòu)整體矩陣和單元矩陣之間的變換關(guān)系[4]，這給基本方法的準確應用帶來了理解上的障礙，容易造成分析錯誤.

本文針對有限元模型中整體矩陣和局部矩陣之間的相互關(guān)系進行研究，以幾何關(guān)系及功能原理為分析基礎，推導了二維坐標系下兩者之間的轉(zhuǎn)換關(guān)系，建立了新的坐標變換方法，所得結(jié)果對采用有限元方法進行框架結(jié)構(gòu)的設計分析具有指導意義.

1　坐標轉(zhuǎn)換方法

1.1抗彎剛度矩陣坐標對應關(guān)系

如圖1所示框架結(jié)構(gòu)桿件單元由兩個節(jié)點組成，每節(jié)點具有二自由度，即豎向位移和轉(zhuǎn)角，不妨以變形之前的局部單元起點a(設此節(jié)點在整體坐標系中的節(jié)點編號為i)為原點，與終點b(在整體坐標系中的節(jié)點編號為j)的連線為坐標軸x,垂直于連線方向的坐標軸為y軸，建立框架構(gòu)件單元的局部坐標系.為方便起見，建議在整體坐標系中結(jié)構(gòu)構(gòu)件各單元的起始編號小于構(gòu)件的終點編號.

圖1　局部坐標系下單元節(jié)點位移

這里只考慮了彎曲，因此單元變形在局部坐標系中表現(xiàn)為節(jié)點沿豎向坐標軸上的位移和繞節(jié)點本身的轉(zhuǎn)角，在局部坐標系中將外力作用下的單元節(jié)點位移矢量記為d，節(jié)點的位移在局部坐標系中可表示為：

(1)

式(1)中，va，vb分別為點a，點b的豎向位移，φa，φb分別為單元在點a，點b處的轉(zhuǎn)角，本文中矢量上標T表示的是向量或者矩陣的轉(zhuǎn)置(下同).

如圖2所示，按右手螺旋定則確定此框架構(gòu)件單元在整體坐標系下的傾角為θ，則單元的方向余弦符合以下關(guān)系：

(2)

式(2)中X,Y是單元節(jié)點在整體坐標系中的位置坐標值，下標為其節(jié)點編號,L是單元長度：

(3)

將單元節(jié)點位移在整體坐標系中記為矢量Dij，則節(jié)點在整體系下的位移為：

(4)

式(4)中U，V為整體坐標系中節(jié)點i和節(jié)點j的位移在兩個坐標軸上的投影，不妨規(guī)定位移的方向與坐標軸正方向相一致取正號，反之則取負號，Φi，Φj分別為構(gòu)件在點a，點b處的轉(zhuǎn)角.

圖2　整體坐標系下單元的節(jié)點位移

由圖2可知，構(gòu)件的變形在局部坐標系和整體坐標系中保持一致，因此有：

(5)

此外有如下三角函數(shù)關(guān)系恒成立：

cos2θ+sin2θ=1

(6)

利用三角函數(shù)代數(shù)關(guān)系式(6)可將局部位移關(guān)系式(1)重寫為如下形式：

(7)

以代數(shù)的形式將整體坐標和局部坐標之間的對應關(guān)系式(5)代入到三角函數(shù)關(guān)系式(7)中(這一步為本文關(guān)鍵)，化簡后可得：

(8)

進一步將關(guān)系式(8)用矩陣的形式表示為：

(9)

根據(jù)式(9)可以將整體位移和局部位移間的對應關(guān)系矩陣S表示為：

(10)

此時可得變換關(guān)系：

d=SDij

(11)

假設框架單元構(gòu)件彈性模量為E，截面慣性矩為I，單元的長度為L，若已知該單元的局部剛度矩陣(是一個對稱矩陣)如下[5]：

(12)

則該單元的彈性應變能(標量形式)在局部坐標系中可寫成內(nèi)積形式：

(13)

此時該單元的彈性應變能在整體坐標系中還可以寫成內(nèi)積形式為：

(14)

將式(12)代入到式(13)中，化簡可得：

(15)

綜合式(13),(14)和(15)可得：

K=STkS

(16)

式(16)即為本文所推導的框架結(jié)構(gòu)有限元模型的整體坐標和局部坐標轉(zhuǎn)換公式，進一步將式(2)和式(10)代入到式(16)中可推導出框架單元在整體坐標系下的剛度矩陣：

(17)

1.2考慮軸向剛度的坐標轉(zhuǎn)換關(guān)系

在上述推導過程中僅考慮了單元的彎曲，即剛度矩陣僅考慮了豎向位移和轉(zhuǎn)角，為了更合理地表示每個節(jié)點自由度對單元矩陣的貢獻，本文擬以上述推導為基礎將軸向荷載作用下構(gòu)件的剛度矩陣和考慮彎曲的剛度矩陣聯(lián)合起來以建立更符合實際情況的坐標變化方法.已知考慮軸向變形的框架單元在局部坐標系下的剛度矩陣為[4]：

(18)

式(18)與式(12)相似，不同之處在于考慮了構(gòu)件軸向剛度的影響，由于同時考慮了框架單元的軸向變形和彎曲，此時構(gòu)件的變形在局部坐標系下表現(xiàn)為節(jié)點沿坐標軸軸向和垂直方向的位移及繞節(jié)點的轉(zhuǎn)角，將外力作用下的節(jié)點位移矢量在局部坐標系中表示為dn，此時在局部坐標系中，單元節(jié)點的局部位移為：

(19)

式(19)中ua，ub分別為點a，點b的豎向位移量，中va，vb分別為點a，點b處的豎向位移量，φa，φb分別為構(gòu)件繞點a和點b的轉(zhuǎn)角.此時節(jié)點位移矢量在整體坐標系下的表達式同式(4)，單元在整體坐標系下的傾角的方向余弦表達式同式(2).

根據(jù)疊加原理，可以將整體坐標系下節(jié)點的位移表示為單元節(jié)點沿豎向坐標軸的位移和轉(zhuǎn)角與單元軸向變形之和：

(20)

式(20)與式(5)相比，疊加了節(jié)點沿單元軸向位移的部分，化簡可得：

(21)

對投影矩陣求逆，進一步整理可得：

(22)

根據(jù)式(22)可得整體位移和局部位移間的變換關(guān)系矩陣Sn為：

(23)

式(23)與式(10)相似，此時由式(22)可得變換關(guān)系：

dn=SnDij

(24)

將單元的彈性應變能在局部坐標系下表示成內(nèi)積形式：

(25)

則該單元的彈性應變能在整體坐標系下表示成內(nèi)積形式如下：

(26)

將式(24)代入式(25)中，化簡得：

(27)

綜合式(25),(26),(27)得：

(28)

式(28)即為考慮軸向變形的框架結(jié)構(gòu)單元剛度矩陣的整體坐標和局部坐標轉(zhuǎn)換方程，進一步將式(18)和式(23)代入到式(28)中推導出整體坐標系下的框架結(jié)構(gòu)單元剛度矩陣如下：

(29)

盡管建立的整體坐標系下的剛度矩陣表達式(29)形式復雜，但是其中的局部坐標和整體坐標變換方程式(23)簡單合理，單元剛度矩陣的局部坐標和整體坐標轉(zhuǎn)換方程式(28)便于理解，適合編程實現(xiàn).

1.3二維荷載向量轉(zhuǎn)換關(guān)系

將節(jié)點的荷載向量在局部坐標系下表示為fn，則外力作功在局部坐標系中的內(nèi)積形式為：

(30)

由功能原理表達式Wn=Pn，式(30)中對位移求導可得局部坐標系中節(jié)點處力的平衡方程：

kndn=fn

(31)

將節(jié)點的荷載向量在整體坐標系下表示為Fn，則與其對應的外力作功在整體坐標系中寫成內(nèi)積的形式為：

(32)

由功能原理表達式Wn=Pn，將式(32)對位移矢量求導可得局部坐標系中單元節(jié)點處力的平衡方程如下：

KnDij=Fn

(33)

將式(24)代入式(30)中化簡得：

(34)

綜合式(24),(30),(34)得框架單元荷載向量的對應關(guān)系：

(35)

式(35)即為本文所推導的框架結(jié)構(gòu)單元荷載向量的整體坐標和局部坐標轉(zhuǎn)換關(guān)系式.根據(jù)式(33)可知整體坐標系中的節(jié)點處力的平衡方程：

(36)

式(36)即為本文建立的整體坐標系下框架結(jié)構(gòu)有限元模型力平衡方程.

由于單元的剛度矩陣和荷載向量是有限元建模中的重要環(huán)節(jié)，因此其推導過程對應用有限元方法求得正確分析結(jié)果非常重要.

2　算例分析

如圖3所示考慮軸向變形平面剛架[6]，其所受到的外荷載已經(jīng)在圖中標出，已知各桿件材料屬性為：EI=2.16×105KN·m2，EA=7.2×106KN，試求節(jié)點位移及支座反力.

圖3　結(jié)構(gòu)示意及位移編碼

第一步，根據(jù)式(18)推導單元剛度矩陣，對于單元①：單元長度L=5m，可得單元①的單元剛度矩陣：

(36)

對于單元③：單元長度L=4m，根據(jù)已知條件代入式(18)，可得單元③的單元剛度矩陣：

(37)

第二步，推導坐標轉(zhuǎn)換矩陣，根據(jù)式(2)求解單元在整體坐標系下的傾角，對于單元①：cosθ1=0.6，sinθ1=0.8(θ1=53.13°)；單元②整體坐標與局部坐標重合，不用進行坐標變換；單元③：cosθ3=0，sinθ3=1(θ3=90°)；根據(jù)式(23)得到單元①的變換矩陣為：

(38)

單元③的變換矩陣為：

(39)

將式(37)和式(39)代入式(28)中求得整體坐標系下單元單元③的剛度矩陣K3.至此，本文通過所推導的框架結(jié)構(gòu)有限元模型的坐標變換方法建立了算例的整體坐標下各單元的剛度矩陣.

第三步，單元組合，將求得的單元剛度矩陣K1，K2和K3組合形成整體剛度矩陣，由于算例結(jié)構(gòu)形式簡單，這里不妨使用Boole連通矩陣[7]或者更為簡便的對號入座法[8]來構(gòu)造結(jié)構(gòu)整體剛度矩陣KG，由于沒有施加邊界條件，此矩陣是一個12階的對稱奇異矩陣.

第四步，應用邊界條件，由于節(jié)點1和節(jié)點2是固定的，即在這些點上的節(jié)點位移和轉(zhuǎn)角為0，將這些支承條件應用到建立的整體剛度矩陣KG中，得到修改后的整體剛度矩陣K如下：

圖4　算例節(jié)點位移

(40)

第六步，施加荷載并建立結(jié)構(gòu)剛度方程，由于荷載均施加在單元②上，單元②整體坐標與局部坐標重合，不需要應用式(36)進行坐標變換，因此由圖3直接可得：

(41)

由單元②的等效節(jié)點荷載可得：

fE=

[0 -45KN -47.5KN·m 0 -45KN 37.5KN·m]T

(42)

從而求出節(jié)點荷載向量為：

f=fd+fE

(43)

由此建立結(jié)構(gòu)剛度方程為：

KD=f

(44)

第七步，求解結(jié)構(gòu)剛度方程得到節(jié)點3和4的變形向量為：

(45)

由式(45)得框架結(jié)構(gòu)的節(jié)點位移如圖4所示：

第八步，求解結(jié)構(gòu)支座反力，由變形向量表達式得到總的位移向量DG，支座反力可以從下列方程求出：

R=KGDG-FG

(43)

具體的表達式如下：

(44)

利用所得到的節(jié)點位移還可以進一步分析框架單元的內(nèi)力和正應力，且在計算過程中單元的整體位移和局部位移依然要通過坐標變換矩陣聯(lián)系起來，由于本文主要討論有限元建模過程中的坐標變換方法，這里不再贅述.

4　結(jié)論

本文圍繞框架結(jié)構(gòu)有限元分析中的建模問題，推導了框架單元局部剛度矩陣和整體剛度矩陣之間的轉(zhuǎn)換方法，結(jié)果表明：

1)利用三角恒等式和疊加原理推導的考慮軸向變形的框架單元的整體剛度矩陣合理且可行；

2)運用功能原理建立的荷載向量的局部坐標和整體坐標轉(zhuǎn)換方程式物理意義明確，便于理解；

3)新的框架結(jié)構(gòu)剛度矩陣模型變換方法簡便快捷，適合編程實現(xiàn).

[1]J N REDDY. An Introduction to the Finite Element Method [M].2nd ed, New York：McGraw-Hill， 1993.

[2]K J BATHE. Finite Element Procedures[M].Englewood Cliffs, NJ：Prentice Hall，1996.

[3]羅帥, 顏全勝. 桁架結(jié)構(gòu)有限元模型的坐標轉(zhuǎn)換方法[J]. 四川建筑科學研究. 2015, 41(3): 10-13.

[4]R D COOK, D S MALKUS,M E PLESHA, et al. Concepts and Applications of Finite Element Analysis [M]. 4th ed.New York：John Wiley & Sons, Inc., 2002: 61-63.

[5]S MOAVENI. Finite Element Analysis - Theory and Application with ANSYS [M]. 2nd ed. Upper Saddle River, NJ：Prentice-Hall, 2002.

[6]王煥定，陳少峰，邊文鳳.有限單元法基礎及MATLAB編程[M].北京：高等教育出版社，2012.

[7]A P BORESI, K P CHONG, S SAIGAL,Approximate solution methods in engineering mechanics [M]. 2nd ed. New York: John Wiley & Sons, Inc., 2003.

[8]朱慈勉，吳宇清，計算結(jié)構(gòu)力學[M].北京：科學出版社，2012.

(責任編輯王海雷)

Coordinate Transformation Method of Finite Element Model for Frame Structures

Luo Shuai1Liu Wei2

(1. School of Civil Engineering, Shaoxing University, Shaoxing, Zhejiang 312000;2. Guangdong Communication Polytechnic, Guangzhou, Guangdong 510650)

Coordinate transformation is an important step in frame structure analysis using Finite Element Method (FEM). In this study, the transformation matrix was deduced by the relationship of trigonometric functions and strain energy, and the load vector was transformed by work-energy theorem. The element stiffness matrices for frame structure were calculated in local coordinate systems and then transformed into global coordinate systems with the new method. Compared with the general derivation method, the new coordinate method results in rigorous logic reasoning and clear physical significance, and it is simple for comprehension and calculation. The example analysis shows that the method is suitable for computer programming.

frame structure model; finite element method; transformation matrix of coordinates; local stiffness matrix; global stiffness matrix

2016-04-15

廣東省自然科學基金資助項目(編號：2015A030310168)，紹興文理學院科研項目(編號：20155028).

羅帥(1981-),男,湖北天門人,助理研究員，博士，主要從事土木工程結(jié)構(gòu)動力學分析及振動控制研究工作.

10.16169/j.issn.1008-293x.k.2016.08.04

O302;TU12

A

1008-293X(2016)08-0017-07