第5章 桿單元和梁單元

1、第五章桿單元和梁單元 本章主要介紹利用桿單元及梁單元進行結(jié)構(gòu)靜力學的有限利用桿單元及梁單元進行結(jié)構(gòu)靜力學的有限元分析原理元分析原理。首先介紹了桿單元的分析方法介紹了桿單元的分析方法，詳細給出了采用桿單元進行有限元分析的整個過程；緊接著介紹了平面梁單元介紹了平面梁單元，以一個平面懸臂梁力學模型為分析實例一個平面懸臂梁力學模型為分析實例，分別采用材料力學材料力學、彈性力學解析計算彈性力學解析計算以及有限元法有限元法進行了分析與求解，以加深讀者對有限元法的理解。桿單元桿單元-桁架結(jié)構(gòu)桁架結(jié)構(gòu)梁單元梁單元-軸系，轉(zhuǎn)子動力學軸系，轉(zhuǎn)子動力學 一般情況下，認為桿件只承受軸向力，只有一個方向的受力和相應(yīng)的變

2、形。本節(jié)將采用有限元法來分析桿件系統(tǒng)，以下給出規(guī)范的有限元法中關(guān)于桿單元的推導過程，以及整個桿系的求解過程。 如圖5-1所示的桿件結(jié)構(gòu)，左端鉸支，右端作用一個集中力，相關(guān)參數(shù)如圖。具體求解過程如下：圖 5-1 桿件結(jié)構(gòu) 待求解的問題 （1）確定坐標系、單元離散，確定位移變量確定坐標系、單元離散，確定位移變量, 外載荷及邊界外載荷及邊界條件。條件。 5.1.1. 一維桿單元一維桿單元 材料力學可輕易求解材料力學可輕易求解 要建立兩種坐標系：單元坐標系（局部坐標系）、整體坐要建立兩種坐標系：單元坐標系（局部坐標系）、整體坐標系。標系。根據(jù)自然離散根據(jù)自然離散, , 坐標系建立成一維坐標系建立成一維

3、, , 單元劃分為兩個單元劃分為兩個, , 給出相應(yīng)的節(jié)點給出相應(yīng)的節(jié)點1 1、2 2、3 3以及相應(yīng)的坐標值（見圖以及相應(yīng)的坐標值（見圖5-15-1）。在局）。在局部坐標系中，取桿單元的左端點為坐標原點，圖部坐標系中，取桿單元的左端點為坐標原點，圖5-25-2為任取的為任取的一個桿單元。一個桿單元。 圖圖 5-2 桿單元桿單元 對于兩個節(jié)點的桿單元，存在如下節(jié)點力和節(jié)點位移的關(guān)對于兩個節(jié)點的桿單元，存在如下節(jié)點力和節(jié)點位移的關(guān)系式系式 1122ePuPuk (5.1) 其中，其中， 稱為稱為單元剛度矩陣單元剛度矩陣ek （2）確定位移模式確定位移模式 假設(shè)單元位移場：2123( )u xaa

4、 xa x 取其線性部分，系數(shù) 、 可由節(jié)點位移 、 確定，稱為位移插值模式（interpolation model）.1a2a1u2u12( )u xaa x(5.2) （3）形函數(shù)矩陣的推導形函數(shù)矩陣的推導 由單元的節(jié)點條件, 兩個節(jié)點坐標為x1、x2,兩個節(jié)點位移為 ， ，代入上式插值模式公式得：11( )|x xu xu22( )|x xu xu12111222aa xuaa xu 求解得到求解得到111121221212()/()()/()aux uuxxauuxx 這樣，這樣， 可以寫成如下矩陣形式可以寫成如下矩陣形式12( )u xaa x12( )1au xxa11122211

5、uxauxa111122211axuaxu 導出導出11112221( )111axuu xxxaxu12( )uxuN=(5.3) 得到得到形函數(shù)矩陣形函數(shù)矩陣（shape function matrix）(5.4) 記節(jié)點位移矢量記節(jié)點位移矢量 (nodal displacement vector)是是12euu(5.5) 11221211( )1(1)1xxxxxxxxxxN 因此，用形函數(shù)矩陣表達的單元內(nèi)任一點的位移函數(shù)是因此，用形函數(shù)矩陣表達的單元內(nèi)任一點的位移函數(shù)是( )( )eu xx N(5.6) （4）應(yīng)變應(yīng)變 由彈性力學的幾何方程知由彈性力學的幾何方程知1維桿單元滿足維桿單

6、元滿足111222d( )11( )deeuuuuxxuuuxxll NB(5.7) （5）應(yīng)力應(yīng)力 由彈性力學的物理方程知：由彈性力學的物理方程知： 12( )( )( )( )eeeeeeeeuEExDxExxull BBS(5.8) （6）利用最小勢能原理導出單元剛度矩陣利用最小勢能原理導出單元剛度矩陣 單元的勢能表達式：單元的勢能表達式：B為為應(yīng)變矩陣（常應(yīng)變）。應(yīng)變矩陣（常應(yīng)變）。 S為為應(yīng)力矩陣（常應(yīng)力）。應(yīng)力矩陣（常應(yīng)力）。 112211202011d2211()2211d22eeeeeeleeeleTTeeeeTeUWuPPuuA dxPPuEAx SBBBP 上式記作如下矩陣

7、形式：上式記作如下矩陣形式：1122eeTeeeTe K P (5.9) 0e eeek P根據(jù)最小勢能原理，根據(jù)最小勢能原理，可以得到可以得到， (5.10) （7）把所有單元按結(jié)構(gòu)形狀進行組集（）把所有單元按結(jié)構(gòu)形狀進行組集（assembly of discrete elements） 對于圖對于圖5.1所示結(jié)構(gòu)所示結(jié)構(gòu) 第一個單元：第一個單元： 1(1)2uu(1)(1)(1)(1)1111EAlK1(1)2RRP其其中，中，單元剛度矩陣單元剛度矩陣（element stiffness matrix），或稱單元特），或稱單元特性矩陣性矩陣(element characteristic m

8、atrix)011d11eeeleTeeeE AKEA xlBB (5.11) 整體結(jié)構(gòu)的總勢能是所有單元的勢能的和，即整體結(jié)構(gòu)的總勢能是所有單元的勢能的和，即 第二個單元：第二個單元： 2(2)3uu(2)(2)(2)(2)1111EAlK2(2)3RFP(1)(1)111(1)(2)12(1)222(2)(2)22223(2)3331111112211111122TTuuuEARRuuuluuuEARFuuul 在這里，把表達成整體位移矢量在這里，把表達成整體位移矢量 的函數(shù)，如下：的函數(shù)，如下：123uuu(1)(1)(1)(1)(1)(1)1111(1)(1)(1)(1)(2)(2)(

9、2)(2)222(1)(1)(2)(2)(2)(2)(2)(2)3333(2)(2)0110220TTEAEAlluuRuEAEAEAEAuuulllluuFuEAEAll 可記作可記作1122TT KP (5.12) 上式的即為上式的即為整體剛度矩陣整體剛度矩陣。即根據(jù)最小勢能原理，由各單元。即根據(jù)最小勢能原理，由各單元剛度矩陣求出的整體剛度矩陣。下式是由整體剛度矩陣表達的系剛度矩陣求出的整體剛度矩陣。下式是由整體剛度矩陣表達的系統(tǒng)方程：統(tǒng)方程：(1)(1)(1)(1)(1)(1)11(1)(1)(1)(1)(2)(2)(2)(2)22(1)(1)(2)(2)(2)(2)(2)(2)33(2

10、)(2)00EAEAllRuEAEAEAEARullllFuEAEAll (5.13) （8）引入邊界條件（引入邊界條件（Treatment of boundary conditions） 為獲取許可位移場，需引入邊界條件為獲取許可位移場，需引入邊界條件1( ): 0BC uu (5.14) 由于由于 ，可劃去它所對應(yīng)的行和列，這樣基于許可位移，可劃去它所對應(yīng)的行和列，這樣基于許可位移場的系統(tǒng)總勢能為場的系統(tǒng)總勢能為10u (1)(1)(2)(2)(2)(2)(1)(2)(2)2223(2)(2)(2)(2)333(2)(2)11022TEAEAEAuuulllFuuuEAEAll 由由最小勢

11、能原理，勢能函數(shù)對未知位移最小勢能原理，勢能函數(shù)對未知位移 求求變分，滿足變分，滿足 的的條件是條件是 ，得如下方程式，得如下方程式23uu230, 0uu= （9）求解節(jié)點位移求解節(jié)點位移 由上式方程可以直接求解得到由上式方程可以直接求解得到 , 注意到注意到R2是內(nèi)是內(nèi)力，不做功。在求解過程中，可以視為力，不做功。在求解過程中，可以視為0。也就是。也就是23uu (5.15)(1)(1)(2)(2)(2)(2)(1)(2)(2)2(2)(2)(2)(2)3(2)(2)EAEAEAullluEAEAll30F （10）求單元應(yīng)變求單元應(yīng)變（5.16） （5.17） （11）各單元應(yīng)力各單元應(yīng)

12、力 利用物理方程，求單元的應(yīng)力利用物理方程，求單元的應(yīng)力(1)(1)(1)(1)( )EExB（5.18） -1(1)(1)(2)(2)(2)(2)(1)(2)(2)2(2)(2)(2)(2)33(2)(2)0=E AEAEAullluFEAEAll1(1)(1)(1)(1)(1)211uull B 2(2)(2)(2)(2)(2)311uull B(2)(2)(2)(2)EEB （12）各支點反力各支點反力 各支反力公式是由單元最小勢能原理得到的，即各支反力公式是由單元最小勢能原理得到的，即 (1)(1)(1)KP(1)(1)11(1)221111uREAuPl（5.19） 為了清楚起見為了

13、清楚起見, , 將上述兩桿結(jié)構(gòu)代入具體數(shù)將上述兩桿結(jié)構(gòu)代入具體數(shù)值：值： ， ， ， ，進行相應(yīng)的單元應(yīng)力計算。得到的結(jié)果如下：，進行相應(yīng)的單元應(yīng)力計算。得到的結(jié)果如下：(1)(2)72 10EEPa(1)(2)222AAcm(1)(2)10cmll1424302.5 107.5 10uumum(1)32.5 10(2)3( )5 10 x (1)0.05MPa(2)110RN 0.1MPa=310NF 1u2u12T12uvuv 1122,x yxy 這里從這里從坐標變換的角度坐標變換的角度出發(fā)來說明平面桿單元的建立出發(fā)來說明平面桿單元的建立方法。從勢能的角度，桿單元的勢能不會因坐標系變換而

14、方法。從勢能的角度，桿單元的勢能不會因坐標系變換而 產(chǎn)生能量的變化。如圖產(chǎn)生能量的變化。如圖5-35-3所示，局部坐標系中桿單元所示，局部坐標系中桿單元1 1維位移維位移和和，可以投影到整體坐標系中變成，可以投影到整體坐標系中變成，兩個節(jié)點的坐標變?yōu)?，兩個節(jié)點的坐標變?yōu)?圖圖 5-2 局部坐標與整體坐標局部坐標與整體坐標的變換的變換 5.1.2 平面桿單元平面桿單元整體坐標系下的位移和局部坐標系下的位移的變換關(guān)整體坐標系下的位移和局部坐標系下的位移的變換關(guān)系為系為 （5.205.20） 式中，坐標變換矩陣為式中，坐標變換矩陣為 1111122222cossin0000cossineuuuvvu

15、uuvvTcossin0000cossine（5.215.21）因此，平面桿單元節(jié)點位移矢量的變換關(guān)系記為因此，平面桿單元節(jié)點位移矢量的變換關(guān)系記為 （5.225.22）單元勢能是一個標量，不會因坐標系的不同而改變。導單元勢能是一個標量，不會因坐標系的不同而改變。導出出整體坐標系下的單元勢能函數(shù)：整體坐標系下的單元勢能函數(shù)： （5.235.23）從上式我們可以導出在整體坐標下平面桿單元的剛度矩從上式我們可以導出在整體坐標下平面桿單元的剛度矩陣陣 eeeT TTTTTTTT1111222211()()221122eeeeeeeeeeeeTeeeeeeeeTeeeeeeeT K P TK T T

16、PTK TTPK PTeeeekTk T（5.245.24）具本而言，在轉(zhuǎn)換矩陣中，有具本而言，在轉(zhuǎn)換矩陣中，有 （5.255.25）式中式中 為單元的長度，為單元的長度， 。這里。這里令令 ， ，則轉(zhuǎn)換矩陣，則轉(zhuǎn)換矩陣T T可表示為可表示為因此，因此，整體坐標系下的單元剛度矩陣為整體坐標系下的單元剛度矩陣為12cos( )exxl12sin( )eyyleL221212()()elxxyy12exxll12eyyml0000elmlm（5.265.26）2222T2222eeeeeeellmllmlmmlmmE AlllmllmlmmlmmkTk T（5.275.27）例5-1，如圖5-4所

17、示的平面桁架結(jié)構(gòu)，在點1處施加有10000N向下的力。材料的彈性模量為E=30106Pa。所有桿的橫截面面積A=2m2，點2、3和點4完全約束。試求解各點的位移、受力以及各個單元的應(yīng)力情況。clear all % clear memory% E; modulus of elasticity % A: area of cross section % L: length of barE=30e6; A=2; EA=E*A;% generation of coordinates and connectivitiesnumberElements=3;numberNodes=4;elementNodes

18、=1 2;1 3;1 4;nodeCoordinates= 0 0;0 120;120 120;120 0;xx=nodeCoordinates(:,1);yy=nodeCoordinates(:,2);GDof=2*numberNodes; % GDof: total number of degrees of freedomdisplacements=zeros(GDof,1);force=zeros(GDof,1);% applied load at node 2force(2)=-10000.0;120120Ge=zeros(4,GDof,numberElements);K=zeros(

19、GDof, GDof);Te=;for i=1:numberElementspos_1=elementNodes(i,1); pos_2=elementNodes(i,2); %節(jié)點編號Le(i)=sqrt(xx(pos_1)-xx(pos_2)2+(yy(pos_1)-yy(pos_2)2);Ke=E*A/Le(i)*1 -1; -1 1; %一維桿單元Te(:,:,i)=(xx(pos_2)-xx(pos_1)/Le(i) (yy(pos_2)-yy(pos_1)/Le(i) 0 0; 0, 0, (xx(pos_2)-xx(pos_1)/Le(i) , (yy(pos_2)-yy(pos

20、_1)/Le(i);Ke_b(:, :, i)=Te(:,:,i)*Ke*Te(:,:,i);Ge(1,2*pos_1-1,i)=1 ; %組集矩陣Ge(2,2*pos_1,i)=1 ; %組集矩陣Ge(3,2*pos_2-1,i)=1 ; %組集矩陣Ge(4,2*pos_2,i)=1 ; %組集矩陣K=K+Ge(:,:,i)*Ke_b(:, :, i)*Ge(:,:,i);endK_s=K(1:2,1:2); %添加邊界后的剛度矩陣force=force(1:2); %外載荷化簡x=inv(K_s)*force %得到位移結(jié)果X=x 0,0,0,0,0,0 %擴展為完整的節(jié)點位移for e=

21、1:numberElementspos_1=elementNodes(e,1); pos_2=elementNodes(e,2); sigma(e)=E*-1 1*Te(:,:,e)*X(2*pos_1-1); X(2*pos_1); X(2*pos_2-1); X(2*pos_2)/Le(e)end利用該程序求得該桁架在各節(jié)點處的位移為：利用該程序求得該桁架在各節(jié)點處的位移為：U =0.0041 -0.0159 0 0 0 0 0 0T 各單元應(yīng)力：1.0e03*3.9645 1.4645 -1.0355 T節(jié)點支反力？在在5.15.1節(jié)及節(jié)及5.25.2節(jié)中，考慮了一維及平面桿單元的情節(jié)中

22、，考慮了一維及平面桿單元的情況，接下來考慮空間桿單元的問題。如圖所示，局部況，接下來考慮空間桿單元的問題。如圖所示，局部坐標系中桿單元坐標系中桿單元1 1維位移維位移 ， ，可以投影到三維的整，可以投影到三維的整體坐標系中，變成體坐標系中，變成 ，兩個節(jié)點的坐，兩個節(jié)點的坐標變?yōu)闃俗優(yōu)?1u2uT112222uvwuvw 111222,x y zxyz5.1.3 空間桿單元空間桿單元兩者之間存在的關(guān)系是兩者之間存在的關(guān)系是 （5.285.28）式中，分別為桿單元在整體坐標系中與各軸的夾角式中，分別為桿單元在整體坐標系中與各軸的夾角 圖圖 5-5 局部局部坐標與整體坐標的變換坐標與整體坐標的變換

23、 11112222coscoscos000000coscoscosuvuwuuvw 在整體坐標系下，空間桿單元剛度矩陣為在整體坐標系下，空間桿單元剛度矩陣為 12cosexxlL12coseyymL12cosezznL000000elmnlmnT根據(jù)上一節(jié)內(nèi)容，令根據(jù)上一節(jié)內(nèi)容，令，則則 222222222222TeeeeeeellmlnllmlnlmmmnlmmmnlnmnnlnmnnE ALllmlnllmlnlmmmnlmmmnlnmnnlnmnnkTk Tcoscoscos000000coscoscoseT（5.29）（5.30）平面懸臂梁問題的解析分析平面懸臂梁問題的解析分析 平面梁

24、單元的分析與求解平面梁單元的分析與求解 新的一類單元，簡化單元，有著廣泛的應(yīng)用新的一類單元，簡化單元，有著廣泛的應(yīng)用拉刀桿鎖緊螺母主軸定位檢出塊主軸定位輪主軸皮帶輪環(huán)圈主軸尾端蓋主軸套管筒主軸鼻端套筒后支撐軸承后支撐軸承前軸承主軸端蓋主軸梁單元的應(yīng)用梁單元的應(yīng)用梁單元的應(yīng)用梁單元的應(yīng)用 作為對照，先用經(jīng)典材料力學法和彈性力學法對平面懸臂作為對照，先用經(jīng)典材料力學法和彈性力學法對平面懸臂梁進行分析求解。梁進行分析求解。 （1） 平面懸臂梁的材料力學求解：平面懸臂梁的材料力學求解： 一端受載荷作用的懸臂梁如圖一端受載荷作用的懸臂梁如圖5-6（a）所示，選取坐標系如）所示，選取坐標系如圖圖5-6 (

25、b)，任意橫截面上的彎矩為，任意橫截面上的彎矩為 c c W y R M (a) 結(jié)構(gòu)示意圖 (b) 力學模型 圖5-6 平面懸臂梁力學模型xLWM（5.31） 受載荷作用后梁產(chǎn)生變形，在受載荷作用后梁產(chǎn)生變形，在xy平面內(nèi)梁的軸線將變成一條平面內(nèi)梁的軸線將變成一條曲線，即撓曲線。根據(jù)材料力學有關(guān)假設(shè)，梁彎曲的撓曲線的近曲線，即撓曲線。根據(jù)材料力學有關(guān)假設(shè)，梁彎曲的撓曲線的近似微分方程為似微分方程為EIMdxvd22（5.32） 由這兩個公式可得撓曲線的微分方程為由這兩個公式可得撓曲線的微分方程為xLWMvEI 積分得積分得CWxxWEIv22DCxxWLxWEIv2326（5.33） 引入邊

26、界條件，左側(cè)固定端引入邊界條件，左側(cè)固定端A處的轉(zhuǎn)角和撓度均等于零，即處的轉(zhuǎn)角和撓度均等于零，即當當x=0時，時， 0AAv0Av （5.34） 把邊界條件式代入式（把邊界條件式代入式（4.22），得），得00AAEIvDEIC 再將所得積分常數(shù)再將所得積分常數(shù)C和和D代回前式，得轉(zhuǎn)角方程和撓曲線方代回前式，得轉(zhuǎn)角方程和撓曲線方程分別為程分別為WLxxWEIv222326xWLxWEIv （5.35） 將懸臂梁的右端受載荷將懸臂梁的右端受載荷W處的橫坐標處的橫坐標x=l代入以上兩式，得代入以上兩式，得右端受載荷截面的轉(zhuǎn)角和撓度分別為右端受載荷截面的轉(zhuǎn)角和撓度分別為 （2）平面懸臂梁的彈性力學求

27、解）平面懸臂梁的彈性力學求解EIWLvBB22EIWLvfBB33（5.36） 末端受集中載荷作用的平面懸臂梁的位移場可以用以下多項末端受集中載荷作用的平面懸臂梁的位移場可以用以下多項式表示式表示 x方向：方向： y方向：方向： 22366),(yxLxEIWyyxu232( , )33()6Wv x yLxxyLxEI 梁的中性面（梁的中性面（y=0的面）上的撓曲為的面）上的撓曲為 受載荷作用的懸臂梁上任何位置處的轉(zhuǎn)角為受載荷作用的懸臂梁上任何位置處的轉(zhuǎn)角為3236)0 ,(xLxEIWxv （5.37） 左側(cè)懸臂處（左側(cè)懸臂處（x=0）的撓曲為）的撓曲為 ，右端處（，右端處（x=L）受到集

28、）受到集中載荷作用，撓曲為中載荷作用，撓曲為 ，該結(jié)果與材料力學中的撓曲線，該結(jié)果與材料力學中的撓曲線公式相同。公式相同。0vEIWLv3322336621yxLxEIWyuxvxy （5.38） 梁中性面（梁中性面（y=0）上的轉(zhuǎn)角為）上的轉(zhuǎn)角為2366)0 ,(xLxEIWxxy 左端點（左端點（x=0）為懸臂點，轉(zhuǎn)角為）為懸臂點，轉(zhuǎn)角為0 xy 受載荷作用的懸臂梁的應(yīng)力場可在應(yīng)變場的基礎(chǔ)上，由彈受載荷作用的懸臂梁的應(yīng)力場可在應(yīng)變場的基礎(chǔ)上，由彈性力學物理方程直接求出性力學物理方程直接求出 右端點（右端點（x=L）為受集中載荷點，轉(zhuǎn)角為）為受集中載荷點，轉(zhuǎn)角為EIWLxy22 受載荷作用的

29、懸臂梁的應(yīng)變場可由彈性力學幾何方程求出受載荷作用的懸臂梁的應(yīng)變場可由彈性力學幾何方程求出xLEIWyxuxyxLEIWyvy)(0 xvyuxy （5.39） yxxE21yxLEIWxLEIWyE221yxLIW012xyyE0 xyxyG （5.40） 注意平面梁僅存注意平面梁僅存x方向應(yīng)力方向應(yīng)力 （1） 建立坐標系，進行單元離散。坐標系包括結(jié)構(gòu)的整體建立坐標系，進行單元離散。坐標系包括結(jié)構(gòu)的整體坐標系和坐標系和單元局部坐標系（利用局部坐標系進行單元分析）。單元局部坐標系（利用局部坐標系進行單元分析）。 （2） 建立平面梁單元的位移模式。建立平面梁單元的位移模式。 設(shè)一個平面梁單元有兩個

30、節(jié)點，如圖設(shè)一個平面梁單元有兩個節(jié)點，如圖5-7所示。在局部坐標所示。在局部坐標系內(nèi)，平面梁單元定義有系內(nèi)，平面梁單元定義有6個自由度個自由度 圖圖5-7平面梁單元模型平面梁單元模型 ,ijeTiizjjzu vu v (5.41)x1=0,x2=L節(jié)點坐標節(jié)點坐標 略去軸向位移，可以設(shè)平面梁單元有如下略去軸向位移，可以設(shè)平面梁單元有如下4個自由度個自由度 ,ijeTizjzvv (5.42) 對于平面梁單元，其彎曲變形的位移場對于平面梁單元，其彎曲變形的位移場 可以設(shè)為下式可以設(shè)為下式)(xv342321)(xxxxv (5.43) 因此，梁的斜率是（因此，梁的斜率是（Hermite型型）2

31、23423zdvxxdx (5.44) 位移模式寫成矩陣形式位移模式寫成矩陣形式1232234( )1 ( )01 23zv xx xxxxx (5.45) （3）推導形推導形函數(shù)矩陣函數(shù)矩陣-直接利用節(jié)點坐標使推導簡單直接利用節(jié)點坐標使推導簡單 代入節(jié)點位移和節(jié)點坐標，有代入節(jié)點位移和節(jié)點坐標，有1(0)vv01|xdvdx2( )v Lv2|x Ldvdx 其中，其中，L梁單元的長度。得到梁單元的長度。得到111223212342223423vvLLLLL 前兩個方程直接解出前兩個方程直接解出 和和 ，代入后兩個方程，解出，代入后兩個方程，解出 和和 ，具體如下，具體如下2341 上面的推

32、導可以寫成如下矩陣形式上面的推導可以寫成如下矩陣形式1121321122412123231()(2)21()()vvvLLvvLL23231111112221112323332222222442211000 0101230011 0101 2323zzxxxvxxvlxxxllxxll 求得求得1231122123322241 01 231 2301iiiiizjjjzjjxxxvxxvxxxxx (5.46) 將式將式(5.46)代入式代入式(5.45)， , 用節(jié)點用節(jié)點的位移形式重新整理，得的位移形式重新整理，得342321)(xxxxv23231122323222( )1 32232x

33、xxxv xvxLLLLxxxxvLLLL 得到的用形函數(shù)矩陣表達的單元內(nèi)任一點的位移是得到的用形函數(shù)矩陣表達的單元內(nèi)任一點的位移是12311111223232111222332222224221 01 23( )11 ( )( )0 10 123231 01 23zezzx xxvxxv xxxxxxxxxvxxxxxxxxx N (5.47) 其中，其中，N(x) 平面梁單元的形函數(shù)。平面梁單元的形函數(shù)。 節(jié)點位移向節(jié)點位移向量，量， 。對于。對于eTeiijjvv( )()() ()() eviivjjxNNNNv 對應(yīng)于撓度的形函數(shù)對應(yīng)于撓度的形函數(shù)N (x) 的的具體表達式是具體表達

34、式是2323223232)( ,23)(2)( ,231)(LxLxNLxLxNLxLxxNLxLxNjjviiv (5.48) (4) 推導應(yīng)變、應(yīng)力，根據(jù)最小勢能原理導出單元剛度矩陣推導應(yīng)變、應(yīng)力，根據(jù)最小勢能原理導出單元剛度矩陣。在這里直接根據(jù)瑞利法，也可以導出以節(jié)點位移的形式來表達梁在這里直接根據(jù)瑞利法，也可以導出以節(jié)點位移的形式來表達梁單元的應(yīng)變能。彎曲梁的應(yīng)變能是單元的應(yīng)變能。彎曲梁的應(yīng)變能是 (5.49) 二階導數(shù)可由方程二階導數(shù)可由方程(5.47)決定決定（利用形函數(shù)位移插值公利用形函數(shù)位移插值公式）式），表示為表示為2222)(dxNddxvdiv22)(dxNdi22)(d

35、xNdjvjjiijvvdxNd)(221234iiejjvBBBBvB (5.50)2222111dddd222dLLLMvUMxEIxEIx 其中其中1234B B B BB22243222322223222162)( ,126)(64)( ,126)(LxLdxNdBLxLdxNdBLxLdxNdBLxLdxNdBjjviiv (5.51) 代入梁單元應(yīng)變能公式，同時假設(shè)對于該單元而言是常量，代入梁單元應(yīng)變能公式，同時假設(shè)對于該單元而言是常量，得單元應(yīng)變能得單元應(yīng)變能1()d2eTTeLUEIx B B (5.52) 節(jié)點位移向量節(jié)點位移向量 不是不是x的的函數(shù)函數(shù), ,上式可以寫成上式

36、可以寫成e1()d2eTTeLUEIxB B 應(yīng)變能的一般形式可以表達成應(yīng)變能的一般形式可以表達成12eTeeU k (5.53) 其中，其中， 平面梁單元的單元剛度矩陣，即平面梁單元的單元剛度矩陣，即ek()eTLEIdxkB B (5.54) 考慮到考慮到B B是是x x的函數(shù)的函數(shù), , 上式所有項積分后得上式所有項積分后得局部坐標系下（代局部坐標系下（代入了坐標值）入了坐標值）的的平面梁單元的單元剛度矩陣平面梁單元的單元剛度矩陣223221261266462()1261266264eLLLLLLEILLLLLLLk (5.55) （5）整體剛度矩陣的組集與坐標變換整體剛度矩陣的組集與坐

37、標變換 a) 局部坐標系向整體坐標系的轉(zhuǎn)換局部坐標系向整體坐標系的轉(zhuǎn)換 局部坐標系，整體坐標系，兩種坐標系下的節(jié)點載荷、節(jié)點局部坐標系，整體坐標系，兩種坐標系下的節(jié)點載荷、節(jié)點位移和單元剛度矩陣的變換關(guān)系為位移和單元剛度矩陣的變換關(guān)系為eeRTReeT1eekT k T 其中坐標變換矩陣為其中坐標變換矩陣為cossin00sincos0000cossin00sincosT(5.57) 式中，式中， 是是 軸軸相對于相對于x軸的夾角?？梢宰C明，轉(zhuǎn)換矩陣軸的夾角?？梢宰C明，轉(zhuǎn)換矩陣T的逆的逆矩陣等于它的轉(zhuǎn)置矩陣，所以，在整體坐標系下的單元剛度矩陣矩陣等于它的轉(zhuǎn)置矩陣，所以，在整體坐標系下的單元剛度

38、矩陣為為eTekT k T(5.58)(5.56)x22ddvyx ey B b) 進行整體剛度矩陣的組集。可以采用直接剛度法進行整體剛度矩陣的組集?？梢圆捎弥苯觿偠确?。 （6）引入約束條件。引入約束條件。 （7）求解系統(tǒng)方程，得到所有的節(jié)點位移。求解系統(tǒng)方程，得到所有的節(jié)點位移。 （8）進而再求出單元的應(yīng)力應(yīng)變等。進而再求出單元的應(yīng)力應(yīng)變等。eyE B例例5-2 平面梁單元應(yīng)用舉例。設(shè)一方形截面的懸臂梁，截面每邊長平面梁單元應(yīng)用舉例。設(shè)一方形截面的懸臂梁，截面每邊長為為5cm，長度為，長度為10m，在左端約束固定，在右端施以一個沿，在左端約束固定，在右端施以一個沿y軸負軸負方向的集中力方向

39、的集中力w=100N，求其撓度與轉(zhuǎn)角。，求其撓度與轉(zhuǎn)角。圖5-8 平面梁單元實例圖 將整個梁分成兩個平面梁單元，求出每個單元的剛度矩陣，然后將兩個單元組集成總體剛度矩陣，引入邊界條件后，再求解出各節(jié)點的撓度和轉(zhuǎn)角。 具體的計算過程可參見以下MATLAB程序。注意到其中的有關(guān)坐標變換部分利用了直接給出的轉(zhuǎn)換矩陣G，其具體含義還可參考第三章中的有關(guān)內(nèi)容。clearx1=0;x2=sym(L);x=sym(x);j=0:3;v=x.jm=. 1 x1 x12 x13 0 1 2*x1 3*x12 1 x2 x22 x230 1 2*x2 3*x22mm=inv(m);N=v*mm%N=1 x x2

40、x3*(inv(m) B=diff(N,2)k=transpose(B)*B;Ke=int(k,0,L)% Element 1: E=4.0e11, I =bh3/12=5.2e-7EI=4.0e11*5.2e-7Ke1=EI*subs(Ke,L,5)Ke2=Ke1T=eye(4,4)Ke1=T*Ke1*T;Ke2=T*Ke2*T;K1=G1*Ke1*G1K2=G2*Ke2*G2K=K1+K2F=0, 0,0,0,-100,0 %u=F*inv(K) u=v1,xta1,v2,xta2,v3,xta3,v1=0 xta1=0%K(1,:)=0;K(:,1)=0;K(2,:)=0;K(:,2)=

41、0;KX=K(3:6,3:6)F(1,1)=0;F(2,1)=0;FX=F(3:6,1)u=inv(KX)*FX利用該程序求得懸臂梁節(jié)點1(左端點)、節(jié)點2（中間點）和節(jié)點3（右端點）處的撓曲和轉(zhuǎn)角為0 0 -0.0501 -0.0180 -0.1603 -0.0240TANSYS的命令流如下：/CONFIG,NRES,1E5/title,analysis of the beam element/prep7!選單元ET,1,beam3!每個節(jié)點 有3個自由度的粱單元!定義材料特性 MP,EX,1,4.0e11 MP,DENS,1,7850MP,PRXY,1,0.3 type, 1b=5e-2!

42、截面的尺寸參數(shù)h=bs=b*hI=(b*h*h*h)/12!截面的慣性矩r,1,s,I,h!繪制模型 用的是4個節(jié)點的板單元n,1,0,0n,2,5,0n,3,10,0e,1,2!自動分配單元號e,2,3d,1,all!加約束f,3,fy,-100!外力/soluSolve/post1!進入后處理PRNSOL, dof !列表顯示節(jié)點位移、應(yīng)力、應(yīng)變PRNSOL, s,compPRNSOL,epel,compANSYS 左端點沿左端點沿y方向位移方向位移(撓曲撓曲)：0 左端點繞左端點繞z軸的轉(zhuǎn)角：軸的轉(zhuǎn)角：0 中間點沿中間點沿y方向位移方向位移(撓曲撓曲)：-0.05 中間點繞中間點繞z軸的

43、轉(zhuǎn)角：軸的轉(zhuǎn)角：-0.018 右端節(jié)點沿右端節(jié)點沿y方向位移方向位移(撓曲撓曲)：-0.16 右端節(jié)點繞右端節(jié)點繞z軸的轉(zhuǎn)角：軸的轉(zhuǎn)角：-0.024利用利用matlab和和ansys兩種方法求得的結(jié)果基本一致兩種方法求得的結(jié)果基本一致:利用前面提到的材料力學公式求得右端點處的撓度值為3117100 1000-0.160256433 4 105.2 10WLyEI 可見，用有限單元法、材料力學方法和用ANSYS計算得到的結(jié)果基本一致。 對于具有兩個節(jié)點的空間梁單元，設(shè)其節(jié)點坐標和相應(yīng)的對于具有兩個節(jié)點的空間梁單元，設(shè)其節(jié)點坐標和相應(yīng)的節(jié)點力如下節(jié)點力如下 （未完）（未完） 節(jié)點（節(jié)點（1）：）：

44、11111111111111TexyzTexyzxyzuvwFNQQMMM (5.59) 節(jié)點（節(jié)點（2）：）：22222222222222TexyzTexyzxyzuvwFNQQMMM (5.60) 5.3.1 空間梁單元的節(jié)點坐標空間梁單元的節(jié)點坐標 整體坐標系記為整體坐標系記為OXYZ，梁單元的局部坐標系記為，梁單元的局部坐標系記為oxyz，其，其中中ox軸正方向由軸正方向由i端截面形心指向端截面形心指向j端面形心，端面形心，y軸和軸和z軸是梁截面軸是梁截面的兩個相互垂直的形心主軸，見的兩個相互垂直的形心主軸，見圖圖5-9。坐標變換公式具有如下坐標變換公式具有如下形式：形式：123123

45、1233 3lllmmmnnnt12 12ttttT (5.61)圖圖5-9 空間梁單元的坐標變換空間梁單元的坐標變換 由局部坐標向整體坐標的位移變換公式是由局部坐標向整體坐標的位移變換公式是TeT (5.62) 節(jié)點力的變換公式是節(jié)點力的變換公式是 TeFT F (5.63) 單元剛度矩陣變換公式是單元剛度矩陣變換公式是TeKT K T (5.64) 在三維空間中，設(shè)在三維空間中，設(shè)x,y,z是局部坐標系，是局部坐標系，X，Y，Z是整體坐標是整體坐標系系, 單元局部坐標系的三個坐標軸的方向余弦分別如下式：單元局部坐標系的三個坐標軸的方向余弦分別如下式：111cos( ,)cos( , )(

46、, )lxXmxYncos xZ222cos( ,)cos( , )cos( , )lyXmyYnyZ333cos( ,)cos( , )cos( , )lzXmzYnzZ (5.65) 坐標變換矩陣的具體求算方法包括如下步驟。坐標變換矩陣的具體求算方法包括如下步驟。 (1) 局部坐標系局部坐標系x軸在整體坐標系中的方向余弦：軸在整體坐標系中的方向余弦：211222212121211222212121211222212121cos( ,)()()()cos( , )()()()( , )()()()xxlxXxxyyzzyymxYxxyyzzzzncos xZxxyyzz (5.66) （2）局部坐標系局部坐標系y軸在整體坐標系中的方向余弦軸在整體坐標系中的方向余弦 現(xiàn)在討論具有任意方向的空間梁單元。首先，由節(jié)點現(xiàn)在討論具有任意方向的空間梁單元。首先，由節(jié)點i、j 在整體坐標系下的坐標即可確定在整體坐標系下的坐標即可確定e1在整體坐標系中的三個方在整體坐標系中的三個方向余弦，即向余弦，即LXXij11LYYij12LZZij13 (5.67) 其中其中2/1222ijijijZZYYXXL (5.68) 下面計算下面計算e2和和e3。 在單元的主慣性平面在單元的主慣性平面oxy上任取一點上任取一點k（但（但k點不能取在點不能取在ox軸軸上），點在整體坐標系中的坐標記為（上），點在