第5章 杆单元和梁单元

1、第五章杆单元和梁单元 本章主要介绍利用杆单元及梁单元进行结构静力学的有限利用杆单元及梁单元进行结构静力学的有限元分析原理元分析原理。首先介绍了杆单元的分析方法介绍了杆单元的分析方法，详细给出了采用杆单元进行有限元分析的整个过程；紧接着介绍了平面梁单元介绍了平面梁单元，以一个平面悬臂梁力学模型为分析实例一个平面悬臂梁力学模型为分析实例，分别采用材料力学材料力学、弹性力学解析计算弹性力学解析计算以及有限元法有限元法进行了分析与求解，以加深读者对有限元法的理解。杆单元杆单元-桁架结构桁架结构梁单元梁单元-轴系，转子动力学轴系，转子动力学 一般情况下，认为杆件只承受轴向力，只有一个方向的受力和相应的变

2、形。本节将采用有限元法来分析杆件系统，以下给出规范的有限元法中关于杆单元的推导过程，以及整个杆系的求解过程。 如图5-1所示的杆件结构，左端铰支，右端作用一个集中力，相关参数如图。具体求解过程如下：图 5-1 杆件结构 待求解的问题 （1）确定坐标系、单元离散，确定位移变量确定坐标系、单元离散，确定位移变量, 外载荷及边界外载荷及边界条件。条件。 5.1.1. 一维杆单元一维杆单元 材料力学可轻易求解材料力学可轻易求解 要建立两种坐标系：单元坐标系（局部坐标系）、整体坐要建立两种坐标系：单元坐标系（局部坐标系）、整体坐标系。标系。根据自然离散根据自然离散, , 坐标系建立成一维坐标系建立成一维

3、, , 单元划分为两个单元划分为两个, , 给出相应的节点给出相应的节点1 1、2 2、3 3以及相应的坐标值（见图以及相应的坐标值（见图5-15-1）。在局）。在局部坐标系中，取杆单元的左端点为坐标原点，图部坐标系中，取杆单元的左端点为坐标原点，图5-25-2为任取的为任取的一个杆单元。一个杆单元。 图图 5-2 杆单元杆单元 对于两个节点的杆单元，存在如下节点力和节点位移的关对于两个节点的杆单元，存在如下节点力和节点位移的关系式系式 1122ePuPuk (5.1) 其中，其中， 称为称为单元刚度矩阵单元刚度矩阵ek （2）确定位移模式确定位移模式 假设单元位移场：2123( )u xaa

4、 xa x 取其线性部分，系数 、 可由节点位移 、 确定，称为位移插值模式（interpolation model）.1a2a1u2u12( )u xaa x(5.2) （3）形函数矩阵的推导形函数矩阵的推导 由单元的节点条件, 两个节点坐标为x1、x2,两个节点位移为 ， ，代入上式插值模式公式得：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) 得到得到形函数矩阵形函数矩阵（shape function matrix）(5.4) 记节点位移矢量记节点位移矢量 (nodal displacement vector)是是12euu(5.5) 11221211( )1(1)1xxxxxxxxxxN 因此，用形函数矩阵表达的单元内任一点的位移函数是因此，用形函数矩阵表达的单元内任一点的位移函数是( )( )eu xx N(5.6) （4）应变应变 由弹性力学的几何方程知由弹性力学的几何方程知1维杆单元满足维杆单

6、元满足111222d( )11( )deeuuuuxxuuuxxll NB(5.7) （5）应力应力 由弹性力学的物理方程知：由弹性力学的物理方程知： 12( )( )( )( )eeeeeeeeuEExDxExxull BBS(5.8) （6）利用最小势能原理导出单元刚度矩阵利用最小势能原理导出单元刚度矩阵 单元的势能表达式：单元的势能表达式：B为为应变矩阵（常应变）。应变矩阵（常应变）。 S为为应力矩阵（常应力）。应力矩阵（常应力）。 112211202011d2211()2211d22eeeeeeleeeleTTeeeeTeUWuPPuuA dxPPuEAx SBBBP 上式记作如下矩阵

7、形式：上式记作如下矩阵形式：1122eeTeeeTe K P (5.9) 0e eeek P根据最小势能原理，根据最小势能原理，可以得到可以得到， (5.10) （7）把所有单元按结构形状进行组集（）把所有单元按结构形状进行组集（assembly of discrete elements） 对于图对于图5.1所示结构所示结构 第一个单元：第一个单元： 1(1)2uu(1)(1)(1)(1)1111EAlK1(1)2RRP其其中，中，单元刚度矩阵单元刚度矩阵（element stiffness matrix），或称单元特），或称单元特性矩阵性矩阵(element characteristic m

8、atrix)011d11eeeleTeeeE AKEA xlBB (5.11) 整体结构的总势能是所有单元的势能的和，即整体结构的总势能是所有单元的势能的和，即 第二个单元：第二个单元： 2(2)3uu(2)(2)(2)(2)1111EAlK2(2)3RFP(1)(1)111(1)(2)12(1)222(2)(2)22223(2)3331111112211111122TTuuuEARRuuuluuuEARFuuul 在这里，把表达成整体位移矢量在这里，把表达成整体位移矢量 的函数，如下：的函数，如下：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) 上式的即为上式的即为整体刚度矩阵整体刚度矩阵。即根据最小势能原理，由各单元。即根据最小势能原理，由各单元刚度矩阵求出的整体刚度矩阵。下式是由整体刚度矩阵表达的系刚度矩阵求出的整体刚度矩阵。下式是由整体刚度矩阵表达的系统方程：统方程：(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) 由于由于 ，可划去它所对应的行和列，这样基于许可位移，可划去它所对应的行和列，这样基于许可位移场的系统总势能为场的系统总势能为10u (1)(1)(2)(2)(2)(2)(1)(2)(2)2223(2)(2)(2)(2)333(2)(2)11022TEAEAEAuuulllFuuuEAEAll 由由最小势

11、能原理，势能函数对未知位移最小势能原理，势能函数对未知位移 求求变分，满足变分，满足 的的条件是条件是 ，得如下方程式，得如下方程式23uu230, 0uu= （9）求解节点位移求解节点位移 由上式方程可以直接求解得到由上式方程可以直接求解得到 , 注意到注意到R2是内是内力，不做功。在求解过程中，可以视为力，不做功。在求解过程中，可以视为0。也就是。也就是23uu (5.15)(1)(1)(2)(2)(2)(2)(1)(2)(2)2(2)(2)(2)(2)3(2)(2)EAEAEAullluEAEAll30F （10）求单元应变求单元应变（5.16） （5.17） （11）各单元应力各单元应

12、力 利用物理方程，求单元的应力利用物理方程，求单元的应力(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、清楚起见, , 将上述两杆结构代入具体数将上述两杆结构代入具体数值：值： ， ， ， ，进行相应的单元应力计算。得到的结果如下：，进行相应的单元应力计算。得到的结果如下：(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 这里从这里从坐标变换的角度坐标变换的角度出发来说明平面杆单元的建立出发来说明平面杆单元的建立方法。从势能的角度，杆单元的势能不会因坐标系变换而

14、方法。从势能的角度，杆单元的势能不会因坐标系变换而 产生能量的变化。如图产生能量的变化。如图5-35-3所示，局部坐标系中杆单元所示，局部坐标系中杆单元1 1维位移维位移和和，可以投影到整体坐标系中变成，可以投影到整体坐标系中变成，两个节点的坐标变为，两个节点的坐标变为 图图 5-2 局部坐标与整体坐标局部坐标与整体坐标的变换的变换 5.1.2 平面杆单元平面杆单元整体坐标系下的位移和局部坐标系下的位移的变换关整体坐标系下的位移和局部坐标系下的位移的变换关系为系为 （5.205.20） 式中，坐标变换矩阵为式中，坐标变换矩阵为 1111122222cossin0000cossineuuuvvu

15、uuvvTcossin0000cossine（5.215.21）因此，平面杆单元节点位移矢量的变换关系记为因此，平面杆单元节点位移矢量的变换关系记为 （5.225.22）单元势能是一个标量，不会因坐标系的不同而改变。导单元势能是一个标量，不会因坐标系的不同而改变。导出出整体坐标系下的单元势能函数：整体坐标系下的单元势能函数： （5.235.23）从上式我们可以导出在整体坐标下平面杆单元的刚度矩从上式我们可以导出在整体坐标下平面杆单元的刚度矩阵阵 eeeT TTTTTTTT1111222211()()221122eeeeeeeeeeeeTeeeeeeeeTeeeeeeeT K P TK T T

16、PTK TTPK PTeeeekTk T（5.245.24）具本而言，在转换矩阵中，有具本而言，在转换矩阵中，有 （5.255.25）式中式中 为单元的长度，为单元的长度， 。这里。这里令令 ， ，则转换矩阵，则转换矩阵T T可表示为可表示为因此，因此，整体坐标系下的单元刚度矩阵为整体坐标系下的单元刚度矩阵为12cos( )exxl12sin( )eyyleL221212()()elxxyy12exxll12eyyml0000elmlm（5.265.26）2222T2222eeeeeeellmllmlmmlmmE AlllmllmlmmlmmkTk T（5.275.27）例5-1，如图5-4所

17、示的平面桁架结构，在点1处施加有10000N向下的力。材料的弹性模量为E=30106Pa。所有杆的横截面面积A=2m2，点2、3和点4完全约束。试求解各点的位移、受力以及各个单元的应力情况。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); %节点编号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 %得到位移结果X=x 0,0,0,0,0,0 %扩展为完整的节点位移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利用该程序求得该桁架在各节点处的位移为：利用该程序求得该桁架在各节点处的位移为：U =0.0041 -0.0159 0 0 0 0 0 0T 各单元应力：1.0e03*3.9645 1.4645 -1.0355 T节点支反力？在在5.15.1节及节及5.25.2节中，考虑了一维及平面杆单元的情节中

22、，考虑了一维及平面杆单元的情况，接下来考虑空间杆单元的问题。如图所示，局部况，接下来考虑空间杆单元的问题。如图所示，局部坐标系中杆单元坐标系中杆单元1 1维位移维位移 ， ，可以投影到三维的整，可以投影到三维的整体坐标系中，变成体坐标系中，变成 ，两个节点的坐，两个节点的坐标变为标变为 1u2uT112222uvwuvw 111222,x y zxyz5.1.3 空间杆单元空间杆单元两者之间存在的关系是两者之间存在的关系是 （5.285.28）式中，分别为杆单元在整体坐标系中与各轴的夹角式中，分别为杆单元在整体坐标系中与各轴的夹角 图图 5-5 局部局部坐标与整体坐标的变换坐标与整体坐标的变换

23、 11112222coscoscos000000coscoscosuvuwuuvw 在整体坐标系下，空间杆单元刚度矩阵为在整体坐标系下，空间杆单元刚度矩阵为 12cosexxlL12coseyymL12cosezznL000000elmnlmnT根据上一节内容，令根据上一节内容，令，则则 222222222222TeeeeeeellmlnllmlnlmmmnlmmmnlnmnnlnmnnE ALllmlnllmlnlmmmnlmmmnlnmnnlnmnnkTk Tcoscoscos000000coscoscoseT（5.29）（5.30）平面悬臂梁问题的解析分析平面悬臂梁问题的解析分析 平面梁

24、单元的分析与求解平面梁单元的分析与求解 新的一类单元，简化单元，有着广泛的应用新的一类单元，简化单元，有着广泛的应用拉刀杆锁紧螺母主轴定位检出块主轴定位轮主轴皮带轮环圈主轴尾端盖主轴套管筒主轴鼻端套筒后支撑轴承后支撑轴承前轴承主轴端盖主轴梁单元的应用梁单元的应用梁单元的应用梁单元的应用 作为对照，先用经典材料力学法和弹性力学法对平面悬臂作为对照，先用经典材料力学法和弹性力学法对平面悬臂梁进行分析求解。梁进行分析求解。 （1） 平面悬臂梁的材料力学求解：平面悬臂梁的材料力学求解： 一端受载荷作用的悬臂梁如图一端受载荷作用的悬臂梁如图5-6（a）所示，选取坐标系如）所示，选取坐标系如图图5-6 (

25、b)，任意横截面上的弯矩为，任意横截面上的弯矩为 c c W y R M (a) 结构示意图 (b) 力学模型 图5-6 平面悬臂梁力学模型xLWM（5.31） 受载荷作用后梁产生变形，在受载荷作用后梁产生变形，在xy平面内梁的轴线将变成一条平面内梁的轴线将变成一条曲线，即挠曲线。根据材料力学有关假设，梁弯曲的挠曲线的近曲线，即挠曲线。根据材料力学有关假设，梁弯曲的挠曲线的近似微分方程为似微分方程为EIMdxvd22（5.32） 由这两个公式可得挠曲线的微分方程为由这两个公式可得挠曲线的微分方程为xLWMvEI 积分得积分得CWxxWEIv22DCxxWLxWEIv2326（5.33） 引入边

26、界条件，左侧固定端引入边界条件，左侧固定端A处的转角和挠度均等于零，即处的转角和挠度均等于零，即当当x=0时，时， 0AAv0Av （5.34） 把边界条件式代入式（把边界条件式代入式（4.22），得），得00AAEIvDEIC 再将所得积分常数再将所得积分常数C和和D代回前式，得转角方程和挠曲线方代回前式，得转角方程和挠曲线方程分别为程分别为WLxxWEIv222326xWLxWEIv （5.35） 将悬臂梁的右端受载荷将悬臂梁的右端受载荷W处的横坐标处的横坐标x=l代入以上两式，得代入以上两式，得右端受载荷截面的转角和挠度分别为右端受载荷截面的转角和挠度分别为 （2）平面悬臂梁的弹性力学求

27、解）平面悬臂梁的弹性力学求解EIWLvBB22EIWLvfBB33（5.36） 末端受集中载荷作用的平面悬臂梁的位移场可以用以下多项末端受集中载荷作用的平面悬臂梁的位移场可以用以下多项式表示式表示 x方向：方向： y方向：方向： 22366),(yxLxEIWyyxu232( , )33()6Wv x yLxxyLxEI 梁的中性面（梁的中性面（y=0的面）上的挠曲为的面）上的挠曲为 受载荷作用的悬臂梁上任何位置处的转角为受载荷作用的悬臂梁上任何位置处的转角为3236)0 ,(xLxEIWxv （5.37） 左侧悬臂处（左侧悬臂处（x=0）的挠曲为）的挠曲为 ，右端处（，右端处（x=L）受到集

28、）受到集中载荷作用，挠曲为中载荷作用，挠曲为 ，该结果与材料力学中的挠曲线，该结果与材料力学中的挠曲线公式相同。公式相同。0vEIWLv3322336621yxLxEIWyuxvxy （5.38） 梁中性面（梁中性面（y=0）上的转角为）上的转角为2366)0 ,(xLxEIWxxy 左端点（左端点（x=0）为悬臂点，转角为）为悬臂点，转角为0 xy 受载荷作用的悬臂梁的应力场可在应变场的基础上，由弹受载荷作用的悬臂梁的应力场可在应变场的基础上，由弹性力学物理方程直接求出性力学物理方程直接求出 右端点（右端点（x=L）为受集中载荷点，转角为）为受集中载荷点，转角为EIWLxy22 受载荷作用的

29、悬臂梁的应变场可由弹性力学几何方程求出受载荷作用的悬臂梁的应变场可由弹性力学几何方程求出xLEIWyxuxyxLEIWyvy)(0 xvyuxy （5.39） yxxE21yxLEIWxLEIWyE221yxLIW012xyyE0 xyxyG （5.40） 注意平面梁仅存注意平面梁仅存x方向应力方向应力 （1） 建立坐标系，进行单元离散。坐标系包括结构的整体建立坐标系，进行单元离散。坐标系包括结构的整体坐标系和坐标系和单元局部坐标系（利用局部坐标系进行单元分析）。单元局部坐标系（利用局部坐标系进行单元分析）。 （2） 建立平面梁单元的位移模式。建立平面梁单元的位移模式。 设一个平面梁单元有两个

30、节点，如图设一个平面梁单元有两个节点，如图5-7所示。在局部坐标所示。在局部坐标系内，平面梁单元定义有系内，平面梁单元定义有6个自由度个自由度 图图5-7平面梁单元模型平面梁单元模型 ,ijeTiizjjzu vu v (5.41)x1=0,x2=L节点坐标节点坐标 略去轴向位移，可以设平面梁单元有如下略去轴向位移，可以设平面梁单元有如下4个自由度个自由度 ,ijeTizjzvv (5.42) 对于平面梁单元，其弯曲变形的位移场对于平面梁单元，其弯曲变形的位移场 可以设为下式可以设为下式)(xv342321)(xxxxv (5.43) 因此，梁的斜率是（因此，梁的斜率是（Hermite型型）2

31、23423zdvxxdx (5.44) 位移模式写成矩阵形式位移模式写成矩阵形式1232234( )1 ( )01 23zv xx xxxxx (5.45) （3）推导形推导形函数矩阵函数矩阵-直接利用节点坐标使推导简单直接利用节点坐标使推导简单 代入节点位移和节点坐标，有代入节点位移和节点坐标，有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)， , 用节点用节点的位移形式重新整理，得的位移形式重新整理，得342321)(xxxxv23231122323222( )1 32232x

33、xxxv xvxLLLLxxxxvLLLL 得到的用形函数矩阵表达的单元内任一点的位移是得到的用形函数矩阵表达的单元内任一点的位移是12311111223232111222332222224221 01 23( )11 ( )( )0 10 123231 01 23zezzx xxvxxv xxxxxxxxxvxxxxxxxxx N (5.47) 其中，其中，N(x) 平面梁单元的形函数。平面梁单元的形函数。 节点位移向节点位移向量，量， 。对于。对于eTeiijjvv( )()() ()() eviivjjxNNNNv 对应于挠度的形函数对应于挠度的形函数N (x) 的的具体表达式是具体表达

34、式是2323223232)( ,23)(2)( ,231)(LxLxNLxLxNLxLxxNLxLxNjjviiv (5.48) (4) 推导应变、应力，根据最小势能原理导出单元刚度矩阵推导应变、应力，根据最小势能原理导出单元刚度矩阵。在这里直接根据瑞利法，也可以导出以节点位移的形式来表达梁在这里直接根据瑞利法，也可以导出以节点位移的形式来表达梁单元的应变能。弯曲梁的应变能是单元的应变能。弯曲梁的应变能是 (5.49) 二阶导数可由方程二阶导数可由方程(5.47)决定决定（利用形函数位移插值公利用形函数位移插值公式）式），表示为表示为2222)(dxNddxvdiv22)(dxNdi22)(d

35、xNdjvjjiijvvdxNd)(221234iiejjvBBBBvB (5.50)2222111dddd222dLLLMvUMxEIxEIx 其中其中1234B B B BB22243222322223222162)( ,126)(64)( ,126)(LxLdxNdBLxLdxNdBLxLdxNdBLxLdxNdBjjviiv (5.51) 代入梁单元应变能公式，同时假设对于该单元而言是常量，代入梁单元应变能公式，同时假设对于该单元而言是常量，得单元应变能得单元应变能1()d2eTTeLUEIx B B (5.52) 节点位移向量节点位移向量 不是不是x的的函数函数, ,上式可以写成上式

36、可以写成e1()d2eTTeLUEIxB B 应变能的一般形式可以表达成应变能的一般形式可以表达成12eTeeU k (5.53) 其中，其中， 平面梁单元的单元刚度矩阵，即平面梁单元的单元刚度矩阵，即ek()eTLEIdxkB B (5.54) 考虑到考虑到B B是是x x的函数的函数, , 上式所有项积分后得上式所有项积分后得局部坐标系下（代局部坐标系下（代入了坐标值）入了坐标值）的的平面梁单元的单元刚度矩阵平面梁单元的单元刚度矩阵223221261266462()1261266264eLLLLLLEILLLLLLLk (5.55) （5）整体刚度矩阵的组集与坐标变换整体刚度矩阵的组集与坐

37、标变换 a) 局部坐标系向整体坐标系的转换局部坐标系向整体坐标系的转换 局部坐标系，整体坐标系，两种坐标系下的节点载荷、节点局部坐标系，整体坐标系，两种坐标系下的节点载荷、节点位移和单元刚度矩阵的变换关系为位移和单元刚度矩阵的变换关系为eeRTReeT1eekT k T 其中坐标变换矩阵为其中坐标变换矩阵为cossin00sincos0000cossin00sincosT(5.57) 式中，式中， 是是 轴轴相对于相对于x轴的夹角。可以证明，转换矩阵轴的夹角。可以证明，转换矩阵T的逆的逆矩阵等于它的转置矩阵，所以，在整体坐标系下的单元刚度矩阵矩阵等于它的转置矩阵，所以，在整体坐标系下的单元刚度

38、矩阵为为eTekT k T(5.58)(5.56)x22ddvyx ey B b) 进行整体刚度矩阵的组集。可以采用直接刚度法进行整体刚度矩阵的组集。可以采用直接刚度法 。 （6）引入约束条件。引入约束条件。 （7）求解系统方程，得到所有的节点位移。求解系统方程，得到所有的节点位移。 （8）进而再求出单元的应力应变等。进而再求出单元的应力应变等。eyE B例例5-2 平面梁单元应用举例。设一方形截面的悬臂梁，截面每边长平面梁单元应用举例。设一方形截面的悬臂梁，截面每边长为为5cm，长度为，长度为10m，在左端约束固定，在右端施以一个沿，在左端约束固定，在右端施以一个沿y轴负轴负方向的集中力方向

39、的集中力w=100N，求其挠度与转角。，求其挠度与转角。图5-8 平面梁单元实例图 将整个梁分成两个平面梁单元，求出每个单元的刚度矩阵，然后将两个单元组集成总体刚度矩阵，引入边界条件后，再求解出各节点的挠度和转角。 具体的计算过程可参见以下MATLAB程序。注意到其中的有关坐标变换部分利用了直接给出的转换矩阵G，其具体含义还可参考第三章中的有关内容。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利用该程序求得悬臂梁节点1(左端点)、节点2（中间点）和节点3（右端点）处的挠曲和转角为0 0 -0.0501 -0.0180 -0.1603 -0.0240TANSYS的命令流如下：/CONFIG,NRES,1E5/title,analysis of the beam element/prep7!选单元ET,1,beam3!每个节点 有3个自由度的粱单元!定义材料特性 MP,EX,1,4.0e11 MP,DENS,1,7850MP,PRXY,1,0.3 type, 1b=5e-2!

42、截面的尺寸参数h=bs=b*hI=(b*h*h*h)/12!截面的惯性矩r,1,s,I,h!绘制模型 用的是4个节点的板单元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 !列表显示节点位移、应力、应变PRNSOL, s,compPRNSOL,epel,compANSYS 左端点沿左端点沿y方向位移方向位移(挠曲挠曲)：0 左端点绕左端点绕z轴的转角：轴的转角：0 中间点沿中间点沿y方向位移方向位移(挠曲挠曲)：-0.05 中间点绕中间点绕z轴的

43、转角：轴的转角：-0.018 右端节点沿右端节点沿y方向位移方向位移(挠曲挠曲)：-0.16 右端节点绕右端节点绕z轴的转角：轴的转角：-0.024利用利用matlab和和ansys两种方法求得的结果基本一致两种方法求得的结果基本一致:利用前面提到的材料力学公式求得右端点处的挠度值为3117100 1000-0.160256433 4 105.2 10WLyEI 可见，用有限单元法、材料力学方法和用ANSYS计算得到的结果基本一致。 对于具有两个节点的空间梁单元，设其节点坐标和相应的对于具有两个节点的空间梁单元，设其节点坐标和相应的节点力如下节点力如下 （未完）（未完） 节点（节点（1）：）：

44、11111111111111TexyzTexyzxyzuvwFNQQMMM (5.59) 节点（节点（2）：）：22222222222222TexyzTexyzxyzuvwFNQQMMM (5.60) 5.3.1 空间梁单元的节点坐标空间梁单元的节点坐标 整体坐标系记为整体坐标系记为OXYZ，梁单元的局部坐标系记为，梁单元的局部坐标系记为oxyz，其，其中中ox轴正方向由轴正方向由i端截面形心指向端截面形心指向j端面形心，端面形心，y轴和轴和z轴是梁截面轴是梁截面的两个相互垂直的形心主轴，见的两个相互垂直的形心主轴，见图图5-9。坐标变换公式具有如下坐标变换公式具有如下形式：形式：123123

45、1233 3lllmmmnnnt12 12ttttT (5.61)图图5-9 空间梁单元的坐标变换空间梁单元的坐标变换 由局部坐标向整体坐标的位移变换公式是由局部坐标向整体坐标的位移变换公式是TeT (5.62) 节点力的变换公式是节点力的变换公式是 TeFT F (5.63) 单元刚度矩阵变换公式是单元刚度矩阵变换公式是TeKT K T (5.64) 在三维空间中，设在三维空间中，设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轴在整体坐标系中的方向余弦轴在整体坐标系中的方向余弦 现在讨论具有任意方向的空间梁单元。首先，由节点现在讨论具有任意方向的空间梁单元。首先，由节点i、j 在整体坐标系下的坐标即可确定在整体坐标系下的坐标即可确定e1在整体坐标系中的三个方在整体坐标系中的三个方向余弦，即向余弦，即LXXij11LYYij12LZZij13 (5.67) 其中其中2/1222ijijijZZYYXXL (5.68) 下面计算下面计算e2和和e3。 在单元的主惯性平面在单元的主惯性平面oxy上任取一点上任取一点k（但（但k点不能取在点不能取在ox轴轴上），点在整体坐标系中的坐标记为（上），点在