羅　帥　劉　偉

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

?

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

羅帥1劉偉2

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

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

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

0　引言

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

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

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

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

1.1抗彎剛度矩陣坐標(biāo)對(duì)應(yīng)關(guān)系

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

圖1　局部坐標(biāo)系下單元節(jié)點(diǎn)位移

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

(1)

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

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

(2)

式(2)中X,Y是單元節(jié)點(diǎn)在整體坐標(biāo)系中的位置坐標(biāo)值，下標(biāo)為其節(jié)點(diǎn)編號(hào),L是單元長(zhǎng)度：

(3)

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

(4)

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

圖2　整體坐標(biāo)系下單元的節(jié)點(diǎn)位移

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

(5)

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

cos2θ+sin2θ=1

(6)

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

(7)

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

(8)

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

(9)

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

(10)

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

d=SDij

(11)

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

(12)

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

(13)

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

(14)

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

(15)

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

K=STkS

(16)

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

(17)

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

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

(18)

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

(19)

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

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

(20)

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

(21)

對(duì)投影矩陣求逆，進(jìn)一步整理可得：

(22)

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

(23)

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

dn=SnDij

(24)

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

(25)

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

(26)

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

(27)

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

(28)

式(28)即為考慮軸向變形的框架結(jié)構(gòu)單元?jiǎng)偠染仃嚨恼w坐標(biāo)和局部坐標(biāo)轉(zhuǎn)換方程，進(jìn)一步將式(18)和式(23)代入到式(28)中推導(dǎo)出整體坐標(biāo)系下的框架結(jié)構(gòu)單元?jiǎng)偠染仃嚾缦拢?/p>

(29)

盡管建立的整體坐標(biāo)系下的剛度矩陣表達(dá)式(29)形式復(fù)雜，但是其中的局部坐標(biāo)和整體坐標(biāo)變換方程式(23)簡(jiǎn)單合理，單元?jiǎng)偠染仃嚨木植孔鴺?biāo)和整體坐標(biāo)轉(zhuǎn)換方程式(28)便于理解，適合編程實(shí)現(xiàn).

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

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

(30)

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

kndn=fn

(31)

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

(32)

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

KnDij=Fn

(33)

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

(34)

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

(35)

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

(36)

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

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

2　算例分析

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

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

第一步，根據(jù)式(18)推導(dǎo)單元?jiǎng)偠染仃?，?duì)于單元①：?jiǎn)卧L(zhǎng)度L=5m，可得單元①的單元?jiǎng)偠染仃嚕?/p>

(36)

對(duì)于單元③：?jiǎn)卧L(zhǎng)度L=4m，根據(jù)已知條件代入式(18)，可得單元③的單元?jiǎng)偠染仃嚕?/p>

(37)

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

(38)

單元③的變換矩陣為：

(39)

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

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

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

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

(40)

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

(41)

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

fE=

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

(42)

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

f=fd+fE

(43)

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

KD=f

(44)

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

(45)

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

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

R=KGDG-FG

(43)

具體的表達(dá)式如下：

(44)

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

4　結(jié)論

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

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

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

3)新的框架結(jié)構(gòu)剛度矩陣模型變換方法簡(jiǎn)便快捷，適合編程實(shí)現(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)有限元模型的坐標(biāo)轉(zhuǎn)換方法[J]. 四川建筑科學(xué)研究. 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]王煥定，陳少峰，邊文鳳.有限單元法基礎(chǔ)及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]朱慈勉，吳宇清，計(jì)算結(jié)構(gòu)力學(xué)[M].北京：科學(xué)出版社，2012.

(責(zé)任編輯王海雷)

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

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

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

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

O302;TU12

A

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