1平面桿件結構的有限單元法

第一篇

桿件結構的有限單元法及程序設計

第1章

平面桿件結構的有限單元法1.1有限單元位移法的基本概念

所謂桿件是指從構造上來說其長度遠大于其截面尺寸的一維構件。在結構力學上我們通常將承受軸力或扭矩的桿件稱為桿，而將承受橫向力和彎矩的桿件稱為梁。在有限單元法中這兩種情況的單元分別稱為桿單元和梁單元。為方便起見，本書都稱之為桿單元。桿系結構是最簡單的一類結構，也是我們在工程上最常見的一類結構。一、有限單元法的基本思路

先分后合，既先將結構劃分成各個單元，進行單元分析，然后再將各單元集合成結構整體，進行整體分析。二、舉例說明

以簡單的結構連續(xù)梁為例，介紹有限元位移法的基本概念和求解方法

圖示簡單的結構連續(xù)梁

1、劃分單元取每一跨為一獨立單元，單元編號為①、②，單元的線剛度為i1,i2

，結點編號為：1，2，3，（稱為整體號）結點角位移為：

1

，

2

，3結點力偶：M1，M2，M3。

建立梁的整體坐標系xoy

根據(jù)右手坐標系規(guī)定，結點角位移和力偶以順時針為正。

每個單元也建立坐標系

，稱為局部坐標系其方向與整體坐標系一致。每一單元的始、末端分別記為i,j端，稱為局部號。

每個單元也建立坐標系

，稱為局部坐標系其方向與整體坐標系一致。每一單元的始、末端分別記為i,j端，稱為局部號。

離散化后的桿件單元和結點，其內(nèi)力只畫出桿端彎矩，各桿端彎矩以順時針為正。

2、建立數(shù)學模型

將結點1，2，3處的結點角位移為基本未知量。得出單元①，②的桿端彎矩和桿端轉角關系式：

單元①

單元②根據(jù)結點1、2、3處位移連續(xù)條件，有代入式（a)、(b)：為本例的位移法方程根據(jù)結點1、2、3處平衡條件有式（d)、(e)代入(f)得：

3、采用矩陣形式

以表示單元序號，則上面的（a）、（b）式可統(tǒng)一寫為：

或簡寫為：

式中：單元剛度矩陣單元桿端位移列陣單元桿端力列陣

或簡寫為：用矩陣形式表示得：

式中：整體剛度矩陣結點位移列陣結點荷載列三、幾點說明

（1）、剛度集成法的應用

建立整體剛度矩陣，常用的方法：剛度集成法（直接剛度法），它是直接利用單元剛度矩陣的“疊加”來形成整體剛度矩陣。

①、將式（1－3）單元剛度矩陣擴階單元①121233單元貢獻矩陣

單元②121233

單元貢獻矩陣相疊加，形成整體剛度矩陣121233121233對于具有n個結點和（n-1）個單元的多跨連續(xù)梁，其整體剛度矩陣如下：

（2）、兩端支承條件的引入

以上在推導連續(xù)梁整體剛度矩陣時，沒有涉及連續(xù)梁兩端有固定端支座的情況。圖示右端為固定端的兩跨連續(xù)梁，先不考慮右端的約束條件，得出整體剛度矩陣與式（1－6）相同。

考慮右端轉角為零的支承條件，求解基本未知量的基本方程為：

為了便于編寫程序，希望引入支承條件后，矩陣的階數(shù)和排列次序不變，，而又達到修正整體剛度方程的目的，將（k）式修改為如下形式：等效于基本方程式（1）和支承條件

連續(xù)梁兩端支承條件的引入方法：將整體矩陣的主對角線元素改為1，第i行、i列的其余元素改為0，對應的荷載元素也改為0，其中i=1或n。

i＝1對應左端為固定支座；i＝n對應右端為固定支座。左右端同時為固定端支座時，應同時進行修改。

（3）、非結點荷載的處理首先在各結點加約束，阻止結點轉動，各桿獨立承擔所受的荷載，桿端產(chǎn)生固端彎矩，記為：

各結點的約束力矩分別為交于該結點的各相關單元的固端力矩之和，以順時針為正。

然后去掉這些附加的約束，這相當于在各結點施加一外力荷載，其大小與約束力矩相同，但方向相反。原非結點荷載的等效結點荷載

疊加既得非結點荷載作用下得各桿桿端彎矩。

非結點荷載作用下的桿端彎矩為兩部分：1、在結點加阻止轉動的約束條件下的固端彎矩；2、在等效結點荷載作用下的桿端彎矩。

例1.1應用有限元位移法求圖示連續(xù)梁的內(nèi)力

解（1）結構離散化，單元及結點編號如圖所示（2）求固端力矩及等效結點荷載三個單元固端力矩為：

等效結點荷載列陣為：（3）求單元剛度矩陣及相應的單元貢獻矩陣

（4）求整體剛度矩陣（5）引入支承條件本例支承條件為：得基本方程：（6）求解基本方程（7）計算各桿端彎矩連續(xù)梁的彎矩圖1.2局部坐標系中的單元剛度矩陣單元剛度方程是單元的桿端力與桿端位移之間的關系式，單元剛度矩陣是單元的桿端位移－桿端力變換矩陣。1.2.1、一般單元

桿單元的長度為，截面面積為，截面慣性矩為，彈性模量為，單元的、端各有三個桿端力為、、，其對應的位移為、、單元桿端力列陣為：單元桿端位移列陣為：為了導出一般單元桿端力與桿端位移的關系，考慮以下兩種情況。

考慮桿端彎矩、和桿端剪力、與桿端轉角位移、和桿端的橫向位移、的關系，可得：將式（a)(b)合在一起，得：簡寫為：其中單元剛度矩陣為：②單元剛度矩陣為對稱矩陣，其元素

①

單元剛度矩陣中的每個元素代表單位桿端位移引起的桿端力。其中的任意元素的物理意義是第個桿端位移分量等于1（其余位移分量等于0）時，所引起的第個桿端力的分量值。1.2.2、單元剛度矩陣的性質(zhì)

③一般單元的單元剛度矩陣是奇異矩陣，它的元素組成的行列式等于零，即。根據(jù)奇異矩陣的性質(zhì)，沒有逆矩陣。也就是說，如果給定桿端位移，根據(jù)（1-18）式可以求出桿端力的惟一解，但反過來，如果已知桿端力，則不能根據(jù)來確定桿端位移的惟一解。因為即使在桿端力已知的情況下，由于單元兩端無任何約束，因此除出桿端自身變形外，還可以發(fā)生任意的剛體位移。

④

單元剛度矩陣具有分塊的性質(zhì)。用虛線把分為四個子矩陣，把和各分為兩個子矩陣，因此，單元剛度矩陣方程（1－17）寫為：

式中：

其任意子矩陣表示桿端力和桿端位移之間的關系。1.2.3、軸力單元只考慮軸向桿端位移和軸向桿端力的單元，稱為軸力單元。為了便于坐標變化，上式寫為：式中：1.3單元剛度矩陣的坐標變換

局部坐標系中的桿端力分量整體坐標系中的桿端力分量兩個坐標系中，力偶分量不變，即：對于單元j端的桿端力：將(a)(b)(c)式子用矩陣形式可表示為：上式可以簡寫成：兩種坐標系中單元桿端力的變換式局部坐標系中的單元桿端力列陣：單元坐標轉換矩陣為：整體坐標系中的單元桿端力列陣：

為正交矩陣，其逆矩陣等于其轉置矩陣，即有：

同樣也適用于兩種坐標系下的桿端位移之間的變換：

將式（1－25）、（1－28）代入式（1－18），得：上式兩邊同乘以，可以得到：令則得：兩種坐標系中單元剛度矩陣的變換公式寫成展開形式為：式中：以上推導方法和步驟完全適用于軸力單元。軸力單元在整體坐標系下的桿端力列陣和桿端位移列陣為：軸力單元中不需要考慮桿端轉角

和桿端力M，坐標變換轉換矩陣為：在整體坐標系下軸力單元中的單元剛度矩陣展開形式為：式中：

單元剛度矩陣中的每個元素代表單位桿端位移引起的桿端力。其中的任意元素的物理意義是第個桿端位移分量等于1（其余位移分量等于0）時，所引起的第個桿端力的分量值。為對稱矩陣，可以用分塊子矩陣表示。在整體坐標系下的單元剛度矩陣與局部坐標系下的單元剛度矩陣有類似的性質(zhì)：對于一般單元，具有奇異性。1.4單元未知量編碼建立結構整體剛度方程求解位置結點位移的方式有：“先處理法”和“后處理法”。1.4.1后處理法由單元剛度矩陣形成整體剛度矩陣，建立剛度方程后再引入支承條件，進而求解結點的未知位移的方法稱為“后處理法”。圖示一平面剛架，設所有結點位移都是未知量。結點位移列陣為：求出各單元剛度方程后，根據(jù)平衡條件和位移連續(xù)條件，可以建立整個結構的位移法方程：式中：結構的整體剛度矩陣應該注意到，在建立整體剛度矩陣方程式（1－38）時，我們假設所有結點都可能位移，相當于整體結構并無支座。因此在外力作用下，除了發(fā)生彈性變形外，還可能發(fā)生剛體位移，這樣各結點位移不能唯一確定。這說明整體剛度矩陣式（1－39）為一奇異矩陣，不能求逆矩陣，故利用整體剛度方程式（1－38）也不能求出結點位移。

實際上，在圖示剛架中，結點1和結點4均為固定端，因此結點位移是已知的。支承條件為：

這樣，將上述支承條件引入到方程中，對整體剛度方程進行修改，可得：

對上述方程進行化簡，分為兩個方程組：

和

這樣，利用第（1－41）式可以求得結點位移和，再根據(jù)第（1－42）2式可以計算未知的支座反力。對于一般桿件結構，都可以按上述步驟進行分析。無論結構具有多少個結點位移分量，經(jīng)過調(diào)整其排列次序，總可以將他們分成兩組：一組包括所有未知結點位移分量（或稱自由結點位移分量），以表示；另一組為支座結點位移分量，以表示。全部結點力分量也分為兩組，與相應者為已知的結點力列陣，以表示；與相應者為支座節(jié)點力列陣，以表示。則整體剛度方程可寫成下列形式:展開上式得:當已知及時，上式可用以求自由結點位移分量；后式可用來計算支座反力。當無支座移動，既時，以上兩式為：修正的整體剛度方程后處理法是將單元剛度矩陣集成整體剛度矩陣后，再引入支承條件予以修改，最后建立用以求解自由結點位移的修正的整體剛度方程。后處理法在對單元、結點的位移、力等未知量進行編號時比較方便。但當結構復雜，例如結構為復式框架、組合結構且支承較為復雜時，應用后處理法解題有諸多不便，此時常采用“先處理法”。

采用“先處理法”解題時，先考慮支承條件，根據(jù)支承條件僅對未知的自由結點位移分量編號，得到的結點位移矩陣中不包含已知的約束結點位移分量，這樣得到的整體剛度方程就是式（1－44）所表示的“修正的整體剛度方程”。用它可以直接求解自由結點的位移分量。1.4.2先處理法具有組合結點的剛架劃分為3個單元，其編號為①、②、③，各桿之上的箭頭表示局部坐標系的軸正方向?？紤]各單元結點的位移分量編號。采用“先處理法”需做如下規(guī)定：①僅對獨立的位移分量按自然順序編號，稱為位移號。若某些位移分量由于聯(lián)結條件或直桿軸向剛性條件（忽略軸向變形）的限制彼此相等，則將他們編為同一位移號。②在支座處，由于剛性約束而使某些位移分量為0，此位移分量編號為0。桿段單元編號單元結點編號單元位移分量編號始端末端AB①120,0,01,2,3BC②231,2,34,5,6DC③540,0,04,5,7在計算機程序中，單元兩端的結點號可采用二維數(shù)組JE（i，e）表示，稱為“單元兩端結點號數(shù)組”JE（i，1）＝單元e始端的結點號JE（i，2）＝單元e末端的結點號本例中：JE（1，1）＝1JE（2，1）＝2JE（1，2）＝2JE（2，2）＝3JE（1，3）＝5JE（2，3）＝4任意結點位移分量的位移號可采用二維數(shù)組JN（i，j）表示，稱為“結點位移數(shù)組”JN（1，j）＝結點j沿x方向的位移號JN（2，j）＝結點j沿y方向的位移號JN（3，j）＝結點j角位移的位移號本例中，對第3結點而言：JN（1，3）＝4JN（2，3）＝5JN（3，3）＝6將單元的始端及末端的位移編碼排成一行（始端在前)，此數(shù)碼稱為“單元定位數(shù)組”。本例中例1.2試對圖示結構的結點和單元進行編號，并用單元兩端的結點號數(shù)組表表示EF桿單元兩端結點號；用結點位移號數(shù)表示F結點的位移分量編號；寫出CD桿單元定位數(shù)組。分別考慮以下兩種情況：（1）考慮各桿的軸向變形；（2）忽略各桿的軸向變形。解：建立整體坐標系，將結構劃分為6個單元，用箭頭表示各單元的局部坐標系軸的正方向，根據(jù)結構特征編為9個結點號。其中聯(lián)結①、③、④單元的鉸接處編有3個結點號，原因是三個單元在鉸接點處具有相同的線位移，但角位移不相同。單元②、③、⑤相交點編有兩個結點號，因為單元②、③在此處為剛性聯(lián)結，具有相同的角位移和線位移。（1）考慮各桿的軸向變形時，各結點位移分量編號如圖所示。EF桿為第⑥單元，單元兩端結點號數(shù)組表示為：JE（1，6）＝8JE（2，6）＝9F點的結點號為9，結點位移號數(shù)組表示為：JN（1，9）＝14JN（2，9）＝15JN（3，9）＝16CD桿為第③單元，單元定位數(shù)組表示為：（1）忽略各桿的軸向變形時，各結點位移分量編號如圖所示。EF桿為第⑥單元，單元兩端結點號數(shù)組不變，仍為：JE（1，6）＝8JE（2，6）＝9F點的結點號為9，結點位移號數(shù)組表示為：JN（1，9）＝8JN（2，9）＝0JN（3，9）＝10CD桿為第③單元，單元定位數(shù)組表示為：1.5平面結構的整體剛度矩陣在進行了單元分析得出單元剛度矩陣之后，需要進行整體分析，以“先處理法”為例，將離散單元重新組合成原結構，使其滿足結構結點的位移連續(xù)條件和力的平衡條件，從而得到修正的結構整體剛度方程，既前面的式（1－44）。修正的結構整體剛度矩陣自由結點位移分量列陣自由結點荷載分量列陣由單元剛度矩陣集成整體剛度矩陣通常采用“直接剛度法”。其計算過程可以分為兩步：首先求出各單元的貢獻矩陣，然后將它們疊加起來，得出整體剛度矩陣。然而在實際電算中，這種做法是不便采用的，因為在計算中需先將所有單元的貢獻矩陣儲存起來，而各單元貢獻矩陣的階數(shù)與整體剛度矩陣的階數(shù)相同，因這樣就要占用大量的計算機存儲容量。故在實際中并不是采用貢獻矩陣法，而是利用各單元的定位數(shù)組，采用“邊定位，邊累加”的方法。

如圖所示的結構有3個單元，5個結點，7個獨立的位移分量，這樣其整體剛度矩陣應為階矩陣，即：

各元素在整體剛度矩陣中的位置為：單元單元剛度矩陣：單元②的剛度矩陣為：

中各元素在整體剛度矩陣中的位置為：

單元③的剛度矩陣為：

中各元素在整體剛度矩陣中的位置為：

這樣，按照以上所講的定位方法，將、、中的相關元素累加到整體剛度矩陣對應的元素上，可以得到整體剛度矩陣為：

在實際電算程中，采用“邊定位，邊累加”的方法，過程如下：

例1.3求圖示剛架的整體剛度矩陣，設各桿截面尺寸相同。

解（1）整理數(shù)據(jù)并進行編號（2）求局部坐標系中單元剛度矩陣（3）求整體坐標系中單元剛度矩陣123000123004（4）形成整體剛度矩陣1.6非結點荷載的處理為分析平面結構而建立的整體剛度方程，反映了結構的結點荷載與結點位移之間的關系。作用在結構上的荷載除了直接作用在結點上的荷載之外，還有作用在桿件上的分布荷載、集中荷載等。這些非結點荷載應轉換成等效結點荷載。將與相疊加，可得綜合結點荷載。綜合結點荷載亦稱“總結點荷載”其下標c通?？陕匀ゲ粚懀矗?/p>

等效結點荷載的計算步驟如下：第一步，在局部坐標系下求單元的固端力。簡圖剪力彎矩QABQBAMABMBA表1－2

兩端固定梁的固端力

第二步，求單元的等效結點荷載。

單元桿端力的變換式：固端力在兩種坐標系下的變換式：有：等效結點荷載列陣可以由下式寫出：展開，得到：

第三步，求整體結構的等效結點荷載。

求得單元等效結點荷載之后，利用單元結點位移分量編號，將中的各分量疊加到結構等效荷載列陣中去。因為中的各元素也是按結點位移分量編號排列的，中的6個元素也與結點位移分量編號一一對應，所以也按對號入座方法，將其逐一累加到相應的位置上去。

例1.4求圖示剛架的等效結點荷載和綜合結點荷載。解：單元局部坐標系及結點位移分量編號如圖所示。（1）：求各單元在局部坐標系中的固端內(nèi)力單元①：由表1－2可得：單元②：因此有：（2）求整體坐標系下各單元的等效結點荷載：單元①，＝900單元②，＝00利用式（1－49）或（1－50）可得：（3）求剛架等效結點荷載列陣，結果為：（4）求直接作用在結點上的荷載，（5）求總結點荷載列陣，應用有限元位移法對桿系結構進行分析，如果按照“先處理法”的思想，其步驟可以歸結為如下：（1）整理原始數(shù)據(jù)，對單元和結點進行編號，確定每個單元的局部坐標系整個結構的整體坐標系。（2）計算局部坐標系中的單元剛度矩陣。（3）按結點位移分量編號，“對號入座，同號疊加”，形成結構整體剛度矩陣。

1.7平面結構分析算例（5）求總的結點荷載列陣。在這一步里，首先必須將非結點荷載轉換成等效結點荷載，再與對應的結點荷載疊加，形成總的結點荷載矩陣。