有限元軟件中通用梁單元的統(tǒng)一形成方法1)

陳 璞 杜 暉 孫樹立 傅向榮

?(北京大學(xué)工學(xué)院力學(xué)與工程科學(xué)系，北京100871)

?(中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，北京100083)

梁單元無疑是有限元分析中應(yīng)用最為廣泛的單元之一。為了滿足各種工程需求，梁單元常?？赡芤缶邆湟韵绿匦曰蚬δ埽?/p>

(1)剪切變形；

(2)端頭剛域，梁?柱交匯處的區(qū)域；

(3)部分剛接，梁?柱交匯等處的剛度損失；

(4)端頭鉸接，梁?柱交匯等處的構(gòu)造或剛度缺失；

(5)截面變化，例如橋梁等大跨結(jié)構(gòu)；

(6)部分偏心，例如與墻柱相連的外立面梁；

(7)幾何剛度，軸向力引起的附加剛度；

(8)軸線為曲線，例如：拱。

在中外文獻(xiàn)中關(guān)于梁單元的形成方法層出不窮，這里不一一列舉，但大部分文獻(xiàn)或只考慮某一特定的單元形式，或形成方案過于繁瑣，不宜在工程軟件中獲得應(yīng)用。本文從工程實(shí)用的角度出發(fā)，著重討論前5類功能在長為L的二維平面梁單元中的統(tǒng)一實(shí)現(xiàn)，讀者不難將其推廣至三維情形。

1 梁單元?jiǎng)偠染仃嚨慕y(tǒng)一形成方式

以Timoshenko梁(T?梁)討論剛度矩陣以及載荷向量的形成方案，Euler梁(E?梁)可以作為T?梁的一個(gè)特例。為了處理變截面梁、端頭剛域以及端頭部分剛接，我們以柔度法為討論重點(diǎn)。

1.1 梁單元與相應(yīng)的簡支梁單元

考慮一個(gè)置于Oxy平面內(nèi)長為L的梁單元，其主軸與x軸重合，左端與原點(diǎn)O重合。與標(biāo)準(zhǔn)的有限元教材[1-2]稍許不同，考慮到單元的通用性，工程軟件中梁單元的4×4剛度矩陣(這里只考慮彎曲部分)往往不是按相應(yīng)的公式直接形成的，而是通過簡支梁或懸臂梁的2×2剛度矩陣擴(kuò)張得到的。為此，對(duì)照梁單元與同跨度的簡支梁(圖1)，它們的撓度和轉(zhuǎn)角滿足變換

圖1 梁單元和簡支梁單元的位移向量

其中w(ξ),θ(ξ),ws(ξ),θs(ξ)分別以自然坐標(biāo)ξ=x/L表示梁單元和相應(yīng)簡支梁的撓度與轉(zhuǎn)角。w(ξ),θ(ξ)滿足T?梁的控制方程，當(dāng)且僅當(dāng)ws(ξ),θs(ξ)滿足該方程。記u(e)={w i,θi,w j,θj}為梁單元兩端的撓度與截面轉(zhuǎn)角，u(e)s={αi,αj}為簡支梁兩端截面轉(zhuǎn)角，則有

其中，S是兩組結(jié)點(diǎn)位移之間的幾何關(guān)系矩陣，其轉(zhuǎn)置是相應(yīng)的兩端外力關(guān)系矩陣，即

其中，f(e)={Q i M i Q j M j}是梁單元兩端的作用力，而f(e)s={M i M j}是相應(yīng)簡支梁兩端的作用力矩，Q i=?Q j=(M i+M j)/L反映了梁端作用力的平衡關(guān)系。由此可知，通過這兩個(gè)關(guān)系矩陣，任意梁單元的變形和梁端作用力均可由一個(gè)相應(yīng)的簡支梁表達(dá)出來，當(dāng)然我們也可以由簡支梁得到普通梁單元的剛度矩陣與等效載荷向量。

1.2 等截面簡支梁的柔度矩陣與剛度矩陣

記ξ=x/L為梁單元的自然坐標(biāo)，(?)′為(?)對(duì)ξ的導(dǎo)數(shù)。由于簡支梁是靜定結(jié)構(gòu)，我們可以通過考慮剪切效應(yīng)的最小余能原理求簡支T?梁的柔度矩陣

其中的彎矩取線性插值

滿足梁平衡關(guān)系。對(duì)于等截面梁，簡支T?梁的柔度矩陣可解析給出

其中無量綱參數(shù)η=12EI/(GAsL2)，表示彎曲與剪切剛度之比。其剛度矩陣可由式(4)求逆獲得

另一方面，我們還可以利用最小勢(shì)能原理獲得梁單元與相應(yīng)簡支梁單元的剛度矩陣變換關(guān)系以及簡支梁單元的剛度矩陣[3]。記H2(ξ)=ξ(ξ?1)2,H4(ξ)=ξ2(ξ?1)為[0,1]上標(biāo)準(zhǔn)的三次Hermite插值函數(shù)中關(guān)于導(dǎo)數(shù)的兩個(gè)，則滿足無分布力簡支T?梁方程的位移插值形函數(shù)可以寫為

其中

基于最小勢(shì)能原理的單元?jiǎng)偠染仃嚍?/p>

由式(7)中的積分，可以獲得與式(5)完全相同的剛度矩陣。式(4)和式(5)互為逆，表明兩端截面轉(zhuǎn)角、力矩和插值所得的跨內(nèi)撓度、彎矩之間的準(zhǔn)確關(guān)系。

當(dāng)η=0時(shí)，剪切剛度無窮大，T?梁回歸到E?梁。作為ξ的函數(shù)，Ns就是Hermite插值H，Ms是其微商H′。此時(shí)，式(5)也退化到E?梁的剛度矩陣。

1.3 變截面梁單元

當(dāng)截面沿中性軸變化時(shí)，位移的三次Hermite插值形函數(shù)不再滿足梁的平衡方程，因此最小勢(shì)能原理所導(dǎo)出的剛度矩陣不再是簡支梁的精確剛度矩陣。但作為靜定結(jié)構(gòu)，簡支梁的彎矩分布仍然是線性的，插值函數(shù)式(3)仍然給出結(jié)構(gòu)力學(xué)意義下的精確內(nèi)力，式(2)是精確的表達(dá)式。由此，最小余能原理所導(dǎo)出的柔度矩陣還是簡支梁的精確柔度矩陣。在這種意義下，變截面梁單元的剛度矩陣可以從柔度矩陣求逆而準(zhǔn)確得到。

式(2)一般可通過較高階的Gauss積分或Lobatto積分獲得。我們的經(jīng)驗(yàn)是5點(diǎn)Gauss積分的精度已足夠工程應(yīng)用。變截面梁的剛度矩陣如采用位移法得到，則可能隱含較大的誤差[4]。

1.4 梁端剛域

在結(jié)構(gòu)力學(xué)中梁?柱的計(jì)算是所謂的中心線計(jì)算。在實(shí)際工程中，梁在其與柱的交匯區(qū)域內(nèi)幾乎沒有彎曲變形，這一部分稱為剛域(圖2)，其長度一般可取為柱的中線到柱外緣距離的75%。

圖2 剛域的簡化

一旦計(jì)入不變形的剛域，梁的彈性變形部分縮短，剛度增大。如圖3，仍記梁柱中心線處簡支梁的轉(zhuǎn)角為又記剛域外緣彈性梁的轉(zhuǎn)角為那么由線性化的幾何關(guān)系得到

圖3 含剛域的梁模型

現(xiàn)在，長為L?c?d的彈性簡支梁兩端的截面轉(zhuǎn)角與外力矩滿足

綜合起來，中心線位置的簡支梁柔度矩陣為

1.5 梁端部分剛接

梁端有時(shí)與柱等其他構(gòu)件的連接不緊密，不能完全傳遞力矩，此時(shí)一般稱為部分剛接，它也可在彈塑性計(jì)算中模擬理想塑性鉸。在有限元計(jì)算中，部分剛接用梁端的附加扭轉(zhuǎn)彈簧來模擬。

較為方便地獲得剛度矩陣的方法仍然是利用最小余能原理，修改不含轉(zhuǎn)角彈簧的柔度矩陣。例如，用最小余能原理得到圖4所示含右端轉(zhuǎn)動(dòng)彈簧的柔度矩陣

圖4 部分剛接簡支梁模型

其逆即為相應(yīng)的剛度矩陣。

由以上討論可以得到梁單元?jiǎng)偠染仃囆纬煞绞剑菏紫扔?jì)算彈性簡支梁的柔度矩陣，如果存在剛域或部分剛接的情形，則需進(jìn)行端頭剛域與集中彈簧修正，獲得，最后通過變換S計(jì)算梁單元的剛度矩陣。

2 等效結(jié)點(diǎn)力

2.1 梁跨中載荷的等效結(jié)點(diǎn)力

梁載荷的等效結(jié)點(diǎn)力大約是工程有限元軟件中最為復(fù)雜的部分。由于工程實(shí)際需求，必須要考慮多種形式的載荷，例如：集中力、均布載荷、三角形分布載荷、梯形分布載荷等，且載荷作用在純彈性梁、組合剛域?彈性梁、組合彈簧?彈性梁。按梁的形式分工況形成載荷無疑是一種可能性，但我們?cè)谶@里再一次從簡支梁出發(fā)，給出一個(gè)簡潔的方案。限于篇幅，僅敘述集中力的作用下的等效載荷，其他分布載荷情況可考慮為單位集中力的等效結(jié)點(diǎn)力的積分。

假設(shè)集中力p作用在梁上，距i端的距離為a，記b=L?a，用力的平衡關(guān)系可得簡支梁跨內(nèi)集中力作用在兩端支座上的等效作用力

用柔度法可得兩端的轉(zhuǎn)角

其中

是簡支梁在集中力p作用下的彎矩。對(duì)于等截面簡支梁[5]可得到

如果是變截面梁，式(12)的數(shù)值積分需要按函數(shù)的光滑性分段進(jìn)行。

用式(13)跨內(nèi)載荷產(chǎn)生的轉(zhuǎn)角，可以得到跨內(nèi)載荷在兩端的等效作用力矩

同時(shí)按照平衡關(guān)系，這兩個(gè)等效力矩也會(huì)帶來在支座上的附加作用力，應(yīng)用式(1)中的轉(zhuǎn)換矩陣S，我們獲得最終等效結(jié)點(diǎn)載荷

對(duì)于等截面梁，等效結(jié)點(diǎn)力的表達(dá)式是

如不考慮η，即剪切剛度無窮大，式(17)退化為等截面E?梁的等效結(jié)點(diǎn)力。

上面的推導(dǎo)過程中，我們沒有區(qū)分是否是等截面梁。如果是變截面梁，封閉形式的等效結(jié)點(diǎn)力一般難于得到，在程序中建議用分段數(shù)值積分計(jì)算。

如果要考慮帶剛域梁的等效結(jié)點(diǎn)力，則需要分別處理作用在彈性部分與剛域部分的跨內(nèi)載荷。我們首先用前述的方法計(jì)算出作用在彈性部分載荷的等效結(jié)點(diǎn)力⌒p(e)，沿用圖3的記號(hào)，記中心線的等效結(jié)點(diǎn)力為p(e)，用平衡即可得到它們之間的轉(zhuǎn)換關(guān)系

如果集中力p作用在左、右端的剛域上，即a≤c或b≤d，那么等效集中力僅出現(xiàn)在一端，分別為

在實(shí)際計(jì)算中并不需要顯式地采用式(17)，而是通過計(jì)算簡支梁的支反力與轉(zhuǎn)角間接獲得等效結(jié)點(diǎn)力，然后再實(shí)施關(guān)于剛域的變換。

特別注意，部分剛接或桿端扭轉(zhuǎn)彈簧對(duì)簡支梁的反力與轉(zhuǎn)角計(jì)算在目前的方案中均無直接影響，它們對(duì)等效結(jié)點(diǎn)載荷的影響間接地通過剛度矩陣的修改計(jì)入。當(dāng)然我們假定，梁的一端不同時(shí)出現(xiàn)矛盾的部分剛接與剛域。

2.2 梁的內(nèi)力求解

如圖5所示，假設(shè)我們已經(jīng)通過求解有限元方程獲得了梁兩端的位移u(e)={w i,θi,w j,θj}，則兩端的真實(shí)作用力可按產(chǎn)生結(jié)點(diǎn)位移所必需的作用力與跨間載荷的等效結(jié)點(diǎn)力疊加計(jì)算，

圖5 梁的內(nèi)力計(jì)算

其中p(e)是等效結(jié)點(diǎn)力[6]。對(duì)于梁單元，這樣得到的端頭作用力與端頭內(nèi)力之間僅有符號(hào)的差異。內(nèi)力輸出時(shí)應(yīng)按工程習(xí)慣進(jìn)行正負(fù)號(hào)的調(diào)整。無論是否含端頭剛域與彈簧，跨內(nèi)分布橫向載荷q(x)的梁內(nèi)力用梁段的平衡關(guān)系給出

3 從簡支梁剛度到整體坐標(biāo)系梁的剛度

各方向定義在圖6中給出。

圖6 梁的方向定義

為了獲得整體坐標(biāo)系下的梁(彎?拉?扭)的12×12單元?jiǎng)偠染仃?，需要?jīng)過兩次變換，第一次是將整體的位移通過坐標(biāo)變換轉(zhuǎn)換到梁局部坐標(biāo)系中其中{i;n1,n2,n3}是梁的右手局部坐標(biāo)系，n1是從i到j(luò)的方向，n2是次彎曲軸方向，n3是主彎曲軸方向。第二次是將局部坐標(biāo)系的位移變換到“簡支”坐標(biāo)系的位移，它的形式是

4 總結(jié)

梁的剛度矩陣與等效結(jié)點(diǎn)載荷的形成方案要考慮多方面的因素，如圖7所示。

圖7 梁單元

梁單元?jiǎng)偠染仃嚨男纬煞桨赴匆韵马樞蜻M(jìn)行:

(1)由式(2)求簡支梁的柔度矩陣，如有剛域或部分剛接情形，需作修正；

(2)求逆獲得簡支梁的剛度矩陣；

(3)按式(7)求梁的剛度矩陣。

前兩步針對(duì)變截面梁(無顯式)，對(duì)于有顯式表達(dá)的等截面梁，可直接寫出Ks，再變換即可。

梁的剛度矩陣與等效結(jié)點(diǎn)載荷在有限元軟件中應(yīng)用非常廣泛，本文從有限元軟件實(shí)現(xiàn)的角度給出了一個(gè)通用梁單元的框架，便于工程軟件開發(fā)人員參考借鑒。以上僅是我們對(duì)有限元梁單元的一點(diǎn)思考，有不當(dāng)與謬誤之處，請(qǐng)同行批評(píng)指正。