羅 帥，朱佳洋，梁超鋒，劉 偉

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

在結(jié)構(gòu)的有限元分析過程中，對于所劃分的單元體，在求得單元?jiǎng)偠确匠讨蟊仨毥⒄麄€(gè)計(jì)算模型的整體剛度方程[1-5]．完成這一步的關(guān)鍵在于怎樣將單元的剛度矩陣組裝成整體剛度矩陣[6-9]．羅少云等[10-11]根據(jù)結(jié)點(diǎn)平衡方程建立整體剛度矩陣，這些方法直觀易懂，但是不能提供編寫有限元分析程序的合理思路；BORESI等[12]使用Boole連通矩陣進(jìn)行剛度矩陣的推導(dǎo)，這種方法數(shù)學(xué)意義直觀，但計(jì)算效率較低；LI等[13]使用對號入座法來構(gòu)造結(jié)構(gòu)整體剛度矩陣，但是這種方法編程效率較低，且不適合處理單元結(jié)點(diǎn)耦合的問題，只適合分析自由度數(shù)較少的結(jié)構(gòu)體系．本研究對有限元模型分析中整體剛度矩陣的集成方法進(jìn)行了分析，采用關(guān)聯(lián)表法建立了有限元整體剛度矩陣的集成方法，所得結(jié)果對加強(qiáng)有限元建模過程的理解以及大規(guī)模有限元模型分析的開展具有指導(dǎo)意義．

1 利用能量原理建立整體矩陣

如圖1所示的平面桁架結(jié)構(gòu)，包含6個(gè)結(jié)點(diǎn)，每個(gè)結(jié)點(diǎn)有2個(gè)自由度，因此結(jié)構(gòu)的整體剛度矩陣為12階，組成桁架結(jié)構(gòu)的桿單元彈性模量為E， 單元截面積為A， 單元長度為L， 設(shè)單元在整體坐標(biāo)系下的傾角為θ， 則在整體坐標(biāo)系下的單元?jiǎng)偠染仃?為對稱矩陣)可表示為[14]

(1)

桿單元軸向變形在整體坐標(biāo)系下表現(xiàn)為結(jié)點(diǎn)沿水平坐標(biāo)軸的位移．將結(jié)點(diǎn)位移矢量在整體坐標(biāo)系下記為D， 則結(jié)點(diǎn)的位移可以表示為

D=[U1V1U2V2…UnVn]T

(2)

其中，U和V分別為桁架結(jié)點(diǎn)的位移在整體坐標(biāo)系兩個(gè)坐標(biāo)軸上的投影，下標(biāo)為單元結(jié)點(diǎn)編號，在此分析模型中n=12．

圖1 桁架結(jié)構(gòu)平面圖Fig.1 Plane truss structure

進(jìn)一步，將所分析的單元(編號為ij)結(jié)點(diǎn)位移矢量在整體坐標(biāo)系下記為Dij， 則該單元節(jié)點(diǎn)的整體位移表達(dá)式為

Dij=[LUiVi…UjVj…]T

(3)

其中，U和V分別為結(jié)點(diǎn)i和結(jié)點(diǎn)j的位移在整體坐標(biāo)系兩個(gè)坐標(biāo)軸上的投影，這里省略了其他結(jié)點(diǎn)的位移．因此有

Dij=Lij·D[UiViUjVj]T

(4)

式(4)中的Lij是一個(gè)構(gòu)造的矩陣，功能是從桁架結(jié)構(gòu)的整體位移向量中提取出與單元?jiǎng)偠染仃囅鄬?yīng)的單元位移向量，其表達(dá)式為

Lij=

(5)

其中，下標(biāo)為元素所在的行和列．

將單元ij的彈性應(yīng)變能(標(biāo)量形式)寫成內(nèi)積形式為

(6)

將式(4)代入式(6)中可得

(7)

已知此單元的彈性應(yīng)變能在整體坐標(biāo)系下保持不變，其內(nèi)積形式為

(8)

對比式(7)和式(8)，可得單元?jiǎng)偠染仃囋谡w坐標(biāo)系下可表示為

Ke=LTijkLij

(9)

由式(9)求得的整體剛度矩陣是一個(gè)過渡性的矩陣，由于剛度矩陣為對稱矩陣，因此可以略去上標(biāo)轉(zhuǎn)置符號．結(jié)構(gòu)的彈性應(yīng)變能在整體坐標(biāo)系下寫成內(nèi)積形式為

(10)

根據(jù)能量原理，結(jié)構(gòu)的彈性應(yīng)變能就是結(jié)構(gòu)內(nèi)部構(gòu)件的應(yīng)變能之和，即

P=∑Pe

(11)

將式(7)代入式(11)，可得

(12)

對比式(10)和式(12)，可得單元?jiǎng)偠染仃嚍?/p>

K=∑LTijkLij

(13)

式(13)即為基于能量原理，在已經(jīng)推導(dǎo)的單元?jiǎng)偠染仃噆的基礎(chǔ)上，通過構(gòu)造功能矩陣Lij推導(dǎo)的有限元模型的整體剛度矩陣建立的方法．該方法揭示了從單元?jiǎng)偠染仃嚨浇Y(jié)構(gòu)整體剛度矩陣的過程，運(yùn)算可知該方法的本質(zhì)就是將單元?jiǎng)偠染仃嚨脑?，按照結(jié)構(gòu)構(gòu)件的聯(lián)接規(guī)律置入整體剛度矩陣的相應(yīng)位置上．該方法具有物理意義明確、推導(dǎo)過程嚴(yán)謹(jǐn)、便于理解的優(yōu)點(diǎn)，但是轉(zhuǎn)換過程較為復(fù)雜，本研究引入關(guān)聯(lián)表以簡化整體矩陣的建立過程．

2 利用關(guān)聯(lián)表建立整體矩陣

由式(9)可知，集成總體剛度矩陣的過程就是將單元?jiǎng)偠染仃囍械脑?，按照?gòu)件結(jié)點(diǎn)的聯(lián)接規(guī)律置入整體剛度矩陣的相應(yīng)位置．

根據(jù)這一規(guī)律，對圖1所示的模型，將各單元的單元編號與整體結(jié)點(diǎn)編號的對應(yīng)關(guān)系列于表1．該結(jié)構(gòu)結(jié)點(diǎn)號碼表確定了每個(gè)單元的結(jié)點(diǎn)位置，這里稱之為關(guān)聯(lián)表，可以把它作為有限元編程過程中要求輸入的基本數(shù)據(jù)之一．隨后利用關(guān)聯(lián)表1，通過定位運(yùn)算得到各結(jié)構(gòu)單元的關(guān)聯(lián)向量(表2)可推導(dǎo)出單元?jiǎng)偠染仃嚫髟卦诮Y(jié)構(gòu)整體剛度矩陣中的對應(yīng)位置，不妨以圖1中的單元6為例來說明具體的推導(dǎo)過程．

表1 局部結(jié)點(diǎn)與整體結(jié)點(diǎn)的關(guān)聯(lián)表

首先，根據(jù)式(1)求得單元6的在整體坐標(biāo)系下的整體剛度矩陣，是一個(gè)4階的對稱矩陣．然后，構(gòu)造12階的零矩陣作為過渡矩陣．接著，根據(jù)式(9)將求得的單元?jiǎng)偠染仃嚫髟靥畛涞竭^渡矩陣的對應(yīng)位置中．由表2可知，其關(guān)聯(lián)向量(行向量)為[3 4 7 8]，構(gòu)造一個(gè)與關(guān)聯(lián)向量同階的全1行向量[1 1 1 1]，求這兩個(gè)行向量的直積[6]，得到所有16個(gè)元素的行位置[5]．交換兩個(gè)向量的順序，再次求得兩個(gè)行向量的直積，得到所有16個(gè)元素的列位置，根據(jù)行列指標(biāo)的定位，將單元?jiǎng)偠染仃嚨脑刂糜谶^渡矩陣的相應(yīng)位置中完成這一步．最后，應(yīng)用式(13)對各個(gè)單元所求得的過渡矩陣求和，得到結(jié)構(gòu)的整體剛度矩陣．

表2 關(guān)聯(lián)向量計(jì)算表

由式(1)可得單元6的單元?jiǎng)偠染仃?，它是一個(gè)4階的對稱矩陣．連接單元6的結(jié)點(diǎn)為結(jié)點(diǎn)2和結(jié)點(diǎn)4，每個(gè)結(jié)點(diǎn)的自由度為2．根據(jù)直積運(yùn)算求得的指標(biāo)向量可知：單元?jiǎng)偠染仃嚨?行、第1列的元素k11在整體剛度矩陣中的位置是第3行、第3列；第1行、第2列的元素k12在整體剛度矩陣中的位置是第3行、第4列；第1行、第3列的元素k13在整體剛度矩陣中的位置是第3行、第7列；第1行、第4列的元素k14在整體剛度矩陣中的位置是第3行、第8列；依此類推，可得第2、3、4行各個(gè)元素在整體矩陣中的位置，從而建立單元6的過渡矩陣kt6． 值得指出的是，當(dāng)以顯示的形式求出每一個(gè)單元的過渡矩陣時(shí)會消耗大量的計(jì)算資源，但是以隱式的形式求和直接得到整體剛度矩陣時(shí)可以顯著提高計(jì)算效率．

由于結(jié)構(gòu)的整體剛度矩陣集成是有限元編程中的重要環(huán)節(jié)，因此其推導(dǎo)過程對運(yùn)用有限元方法獲得正確分析結(jié)果至關(guān)重要，本研究采用的關(guān)聯(lián)表方法簡化了整體矩陣的建立過程，能有效節(jié)約計(jì)算資源，適合編程實(shí)現(xiàn),對大型有限元模型分析的開展具有指導(dǎo)意義．

3 算例分析

圖1所示的平面桁架，其所受到的外荷載已在圖中標(biāo)出，已知各桿件的材料物理屬性為EA=2×108N， 試求各桿件內(nèi)力．

1)建立單元?jiǎng)偠染仃嚕鶕?jù)圖1中各節(jié)點(diǎn)的幾何坐標(biāo)，采用式(1)推導(dǎo)各單元?jiǎng)偠染仃噆ei． 以單元6為例，由單元6的整體坐標(biāo)位置可得單元長度L=4 m，單元傾角θ=36.87°， 由式(1)得到單元6的單元?jiǎng)偠染仃?為簡化這里略去單位，下同)為

(14)

2)建立結(jié)構(gòu)各單元的過渡矩陣．建立與整體剛度矩陣同階的全零矩陣作為過渡矩陣(本例模型有6個(gè)結(jié)點(diǎn)，每結(jié)點(diǎn)2個(gè)自由度，共12階)，根據(jù)單元節(jié)點(diǎn)的連接關(guān)系建立結(jié)構(gòu)關(guān)聯(lián)向量(表2)．以結(jié)點(diǎn)6為例,其關(guān)聯(lián)向量(行向量)為[3 4 7 8]，構(gòu)造一個(gè)與關(guān)聯(lián)向量同階的全1行向量[1 1 1 1]，求這兩個(gè)行向量的直積[15](也叫Kronnecker張量積)，得到所有16個(gè)元素的行位置為

[3 3 3 3 4 4 4 4 7 7 7 7 8 8 8 8]．

交換兩個(gè)向量的順序，再次求得兩個(gè)行向量的直積，得到所有16個(gè)元素的列位置為

[3 4 7 8 3 4 7 8 3 4 7 8 3 4 7 8]．

把單元6的剛度矩陣各元素按行排成一個(gè)行向量,根據(jù)得到的行列指標(biāo)定位，將此行向量的各個(gè)元素置于過渡矩陣kt6的相應(yīng)位置,得到相應(yīng)單元的過渡矩陣為

(15)

3)建立結(jié)構(gòu)整體剛度矩陣．將所求得的各單元過渡矩陣加起來得到結(jié)構(gòu)整體剛度矩陣為

(16)

值得指出的是,因?yàn)檫€沒有考慮邊界約束，此步驟得到的整體剛度矩陣是一個(gè)奇異矩陣．

4)施加邊界條件,由于結(jié)點(diǎn)1在水平和豎直方向的位移被約束,即在結(jié)點(diǎn)1上水平和豎直方向的位移為0，結(jié)點(diǎn)4在豎直方向的位移被約束．即在結(jié)點(diǎn)4上豎直方向的位移為0,采用劃0置1法[16]將這些邊界條件應(yīng)用到建立的整體剛度矩陣KG中,得到修改后的整體剛度矩陣K為

(17)

5)施加荷載并建立結(jié)構(gòu)剛度方程．由于荷載均施加在結(jié)點(diǎn)6的豎直向下方向上，因此由圖1直接可得

f=[0 0 0 0 0 0 0 0 0 0 0 -10]T

(18)

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

KD=f

(19)

6)計(jì)算單元內(nèi)力．將式(17)和式(18)代入式(19)中求解結(jié)構(gòu)剛度方程，得到各結(jié)點(diǎn)的位移，圖2中的虛線即結(jié)點(diǎn)位移后的結(jié)構(gòu)變形圖．根據(jù)得到的結(jié)點(diǎn)位移可求得桁架桿單元的內(nèi)力，見圖2中數(shù)值．由于本研究主要討論有限元建模過程中的整體矩陣集成方法，這里不再贅述．

圖2 單元內(nèi)力(單位：N)Fig.2 Diagram of element internal force(unit：N)

結(jié) 語

本研究圍繞結(jié)構(gòu)有限元分析中的建模問題,推導(dǎo)了基于關(guān)聯(lián)表的有限元模型整體剛度矩陣集成方法，結(jié)果表明，運(yùn)用功能原理推導(dǎo)結(jié)構(gòu)的整體剛度矩陣物理意義明確，便于理解；運(yùn)用關(guān)聯(lián)表建立的結(jié)構(gòu)整體剛度矩陣有效可行；本研究基于關(guān)聯(lián)表方法建立的結(jié)構(gòu)整體剛度矩陣不僅適合桁架結(jié)構(gòu)，對其他結(jié)構(gòu)也適用．