[工学]第四章-杆单元和梁单元

1、 本章主要介绍利用杆单元及梁单元进行结构静力学的有限利用杆单元及梁单元进行结构静力学的有限元分析原理元分析原理。首先介绍了杆单元的分析方法介绍了杆单元的分析方法，详细给出了采用杆单元进行有限元分析的整个过程；紧接着介绍了平面梁单元介绍了平面梁单元，以一个平面悬臂梁力学模型为分析实例一个平面悬臂梁力学模型为分析实例，分别采用材料力学材料力学、弹性力学解析计算弹性力学解析计算以及有限元法有限元法进行了分析与求解，以加深读者对有限元法的理解。 杆件只承受轴向力，可以视为一种特殊的梁单元，本节将采用有限元法来分析杆件系统，以下给出规范的有限元法中关于杆单元的推导过程，以及整个杆系的求解过程。 如图4-

2、1所示的杆件结构，左端铰支，右端作用一个集中力，相关参数如图。具体求解过程如下：图 4-1 杆件结构 （1）确定坐标系、单元离散，确定位移变量）确定坐标系、单元离散，确定位移变量, 外载荷及边界外载荷及边界条件。条件。 要建立两种坐标系：单元坐标系（局部坐标系）、整体坐要建立两种坐标系：单元坐标系（局部坐标系）、整体坐标系。根据自然离散标系。根据自然离散, , 坐标系建立成一维坐标系建立成一维, , 单元划分为两个单元划分为两个, , 给出相应的节点给出相应的节点1 1、2 2、3 3以及相应的坐标值（见图以及相应的坐标值（见图4-14-1）。在局）。在局部坐标系中，取杆单元的左端点为坐标原点

3、，图部坐标系中，取杆单元的左端点为坐标原点，图4-24-2为任取的为任取的一个杆单元。一个杆单元。 图图 4-2 杆单元杆单元 对于两个节点的杆单元，存在如下节点力和节点位移的关对于两个节点的杆单元，存在如下节点力和节点位移的关系式系式 1122ePuPuk (4.1) 其中，其中， 称为单元刚度矩阵称为单元刚度矩阵ek （2）确定位移模式确定位移模式 假设单元位移场：2123( )u xaa xa x 取其线性部分，系数 、 可由节点位移 、 确定，称为位移插值模式（interpolation model）.1a2a1u2u12( )u xaa x(4.2) （3）形函数矩阵的推导形函数矩阵

4、的推导 由单元的节点条件, 两个节点坐标为x1、x2,两个节点位移为 ， ，代入上式插值模式公式得：11( )|x xu xu22( )|x xu xu12111222aa xuaa xu 求解得到求解得到111121221212()/()()/()aux uuxxauuxx 这样，这样， 可以写成如下矩阵形式可以写成如下矩阵形式12( )u xaa x12( )1au xxa11122211uxauxa111122211axuaxu 导出导出11112221( )111axuu xxxaxu12( )uxuN=(4.3) 得到形函数矩阵（得到形函数矩阵（shape function matr

5、ixshape function matrix）2121( )(1)xxxxxxxN(4.4) 记节点位移矢量记节点位移矢量 (nodal displacement vector)(nodal displacement vector)是是12euu(4.5) 因此，用形函数矩阵表达的单元内任一点的位移函数是因此，用形函数矩阵表达的单元内任一点的位移函数是( )( )eu xx N(4.6) （4）应变应变 由弹性力学的几何方程知由弹性力学的几何方程知1维杆单元满足维杆单元满足111222d( )11( )deeuuuuxxuuuxxll NB(4.7) （5）应力应力 由弹性力学的物理方程知：

6、由弹性力学的物理方程知： 12( )( )( )eeeeeeeuEExExxull BS(4.8) （6）利用最小势能原理导出单元刚度矩阵利用最小势能原理导出单元刚度矩阵 单元的势能表达式：单元的势能表达式：112211202011( ) ( )d2211( ( )( ) ( ( )( )2211d22eeeeeeleTeeleTTeeeeTeUWuxxPPuuxxxx A dxPPuEA x SBBBP 上式记作如下矩阵形式：上式记作如下矩阵形式：1122eeTeeeTe K P (4.9) 其中，单元刚度矩阵（其中，单元刚度矩阵（element stiffness matrix），或称单）

7、，或称单元特性矩阵元特性矩阵(element characteristic matrix)011d11eeeleTeeeE AKEA xlBB (4.10) 根据最小势能原理，根据最小势能原理， ，得，得0e eeeK P (4.11) 其中节点载荷矩阵为其中节点载荷矩阵为 12ePPP （7）把所有单元按结构形状进行组集（）把所有单元按结构形状进行组集（assembly of discrete elements） 对于图对于图4.1所示结构所示结构 第一个单元：第一个单元： 1(1)2uu(1)(1)(1)(1)1111EAlK1(1)2RRP 整体结构的总势能是所有单元的势能的和，即整体结

8、构的总势能是所有单元的势能的和，即 第二个单元：第二个单元： 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)(2)(2)2222(1)(1)(2)(2)(2)(2)(2)(2)3333(2)(2)011220TTEAEAlluuRuEA

9、EAEAEAuuRulllluuFuEAEAll 可记作可记作1122TT KP (4.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)(2)00EAEAllRuEAEAEAEARullllFuEAEAll (4.13) （8）引入边界条件（引入边界条件（T

10、reatment of boundary conditions） 为获取许可位移场，需引入边界条件为获取许可位移场，需引入边界条件1( ): 0BC uu (4.14) 由于由于 ，可划去它所对应的行和列，这样基于许可位移，可划去它所对应的行和列，这样基于许可位移场的系统总势能为场的系统总势能为10u (1)(1)(2)(2)(2)(2)(1)(2)(2)2223(2)(2)(2)(2)333(2)(2)11022TEAEAEAuuulllFuuuEAEAll （9）建立系统弹性方程建立系统弹性方程 由最小势能原理，势能函数对未知位移由最小势能原理，势能函数对未知位移 求变分，满足求变分，满足

11、 的条件是的条件是 ，得如下方程式，得如下方程式23uu23,minu u230, 0uu23RF(1)(1)(2)(2)(2)(2)(1)(2)(2)2(2)(2)(2)(2)3(2)(2)EAEAEAullluEAEAll= (4.15) （10）求解节点位移求解节点位移 由上式方程可以直接求解得到由上式方程可以直接求解得到 , 注意到注意到R2是内是内力，不做功。在求解过程中，可以视为力，不做功。在求解过程中，可以视为0。也就是。也就是23uu （11）求单元应变求单元应变(1)(1)(2)(2)(2)(2)(1)(2)(2)2(2)(2)(2)(2)3(2)(2)EAEAEAulllu

12、EAEAll30F（4.16） =1(1)(1)(1)(1)(1)211( )( )uxxull B2(2)(2)(2)(2)(2)311( )( )uxxull B（4.17） （12）各单元应力各单元应力 利用物理方程，求单元的应力利用物理方程，求单元的应力(1)(1)(1)(1)( )EExB（4.18） （13）各支点反力各支点反力 各支反力公式是由单元最小势能原理得到的，即各支反力公式是由单元最小势能原理得到的，即 (1)(1)(1)KP(1)(1)11(1)221111uREAuPl（4.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 作为对照，先用经典材料力学法和弹性力学法对平面悬臂作为对照，先用经典材料力学法和弹性力学法对平面悬臂梁进行分析求解。梁进行分析求解。 （1） 平面悬臂梁的材料力学求解：平面悬臂梁的材料力学求解： 一端受载荷作用的悬臂梁如图一端受载荷作用的悬臂梁如图4-6（a）所

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

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

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

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

18、受集中载荷点，转角为EIWLxy22 受载荷作用的悬臂梁的应变场可由弹性力学几何方程求出受载荷作用的悬臂梁的应变场可由弹性力学几何方程求出xLEIWyxuxyxLEIWyvy)(0 xvyuxy （4.28） yxxE21yxLEIWxLEIWyE221yxLIW012xyyE0 xyxyG （4.29） （1） 建立坐标系，进行单元离散。坐标系包括结构的整体建立坐标系，进行单元离散。坐标系包括结构的整体坐标系和单元局部坐标系。坐标系和单元局部坐标系。 （2） 建立平面梁单元的位移模式。建立平面梁单元的位移模式。 设一个平面梁单元有两个节点，如图设一个平面梁单元有两个节点，如图4-5所示。在局

19、部坐标所示。在局部坐标系内，平面梁单元定义有系内，平面梁单元定义有6个自由度个自由度 图图4-5平面梁单元模型平面梁单元模型 ,ijeTiizjjzu vu v (4.30) 略去轴向位移，可以设平面梁单元有如下略去轴向位移，可以设平面梁单元有如下4个自由度个自由度 ,ijeTizjzvv (4.31) 对于平面梁单元，其弯曲变形的位移场对于平面梁单元，其弯曲变形的位移场 可以设为下式可以设为下式)(xv342321)(xxxxv (4.32) 因此，梁的斜率是（因此，梁的斜率是（Hermite型）型）223423zdvxxdx (4.33) 位移模式写成矩阵形式位移模式写成矩阵形式12322

20、34( )1 ( )01 23zv xx xxxxx (4.34) （3）推导形函数矩阵）推导形函数矩阵 代入节点位移和节点坐标，有代入节点位移和节点坐标，有ivv)0(0|xidvdxjvLv)(|x Ljdvdx 其中，其中，L梁单元的长度。得到梁单元的长度。得到24323423212132LLLLLvvjjii 前两个方程直接解出前两个方程直接解出 和和 ，代入后两个方程，解出，代入后两个方程，解出 和和 ，具体如下，具体如下2341 上面的推导可以写成如下矩阵形式上面的推导可以写成如下矩阵形式123243231()(2)21()()iijiijijijvvvLLvvLL23231122

21、2232333224411000 01 23010011 01230123iiiiziiijjjjz jjjvxxxxxvlxxxllxxll 求得求得123122233241 01 231 2301iiiiziiijjjjz jjjvxxxxxvxxxxx (4.35) 将式将式(4.35)代入式代入式(4.34)， , 用节点用节点的位移形式重新整理，得的位移形式重新整理，得342321)(xxxxvjjiiLxLxvLxLxLxLxxvLxLxxv2323223232232231)( 得到的用形函数矩阵表达的单元内任一点的位移是得到的用形函数矩阵表达的单元内任一点的位移是12312232

22、3222233241 0 1 23( )11 ( )( )0 10 123231 2301iiiiziiiejzjjjzjjjvx xxxxv xxxxxxxxvxxxxxxxxxx N (4.36) 其中，其中，N(x) 平面梁单元的形函数。平面梁单元的形函数。 节点位移向节点位移向量，量， 。对于。对于eTeiijjvv( )()() ()() eviivjjxNNNNv 式中式中N(x)的具体表达式是的具体表达式是2323223232)( ,23)(2)( ,231)(LxLxNLxLxNLxLxxNLxLxNjjviiv (4.37) (4) 推导应变、应力，根据最小势能原理导出单元刚

23、度矩阵推导应变、应力，根据最小势能原理导出单元刚度矩阵。在这里直接根据瑞利法，也可以导出以节点位移的形式来表达梁在这里直接根据瑞利法，也可以导出以节点位移的形式来表达梁单元的应变能。弯曲梁的应变能是单元的应变能。弯曲梁的应变能是dxdxvdEIUL22221 (4.37) 二阶导数可由方程二阶导数可由方程(4.36)决定，表示为决定，表示为2222)(dxNddxvdiv22)(dxNdi22)(dxNdjvjjiijvvdxNd)(221234iiejjvBBBBvB (4.38) 其中其中1234B B B BB22243222322223222162)( ,126)(64)( ,126)

24、(LxLdxNdBLxLdxNdBLxLdxNdBLxLdxNdBjjviiv (4.39) 代入梁单元应变能公式，同时假设对于该单元而言是常量，代入梁单元应变能公式，同时假设对于该单元而言是常量，得单元应变能得单元应变能1()d2eTTeLUEIx B B (4.40) 节点位移向量节点位移向量 不是不是x x的函数的函数, ,上式可以写成上式可以写成e1()d2eTTeLUEIxB B 应变能的一般形式可以表达成应变能的一般形式可以表达成12eTeeU k (4.41) 其中，其中， 平面梁单元的单元刚度矩阵，即平面梁单元的单元刚度矩阵，即ek()eTLEIdxkB B (4.42) 考虑

25、到考虑到B B是是x x的函数的函数, , 上式所有项积分后得局部坐标系下的平上式所有项积分后得局部坐标系下的平面梁单元的单元刚度矩阵面梁单元的单元刚度矩阵223221261266462()1261266264eLLLLLLEILLLLLLLk (4.43) （5）整体刚度矩阵的组集与坐标变换整体刚度矩阵的组集与坐标变换 a) 局部坐标系向整体坐标系的转换局部坐标系向整体坐标系的转换 局部坐标系，整体坐标系，两种坐标系下的节点载荷、节点局部坐标系，整体坐标系，两种坐标系下的节点载荷、节点位移和单元刚度矩阵的变换关系为位移和单元刚度矩阵的变换关系为eeRTReeT1eekT k T 其中坐标变换

26、矩阵为其中坐标变换矩阵为cossin00sincos0000cossin00sincosT(4.45) 式中，式中， 是是x轴相对于轴相对于x轴的夹角。可以证明，转换矩阵轴的夹角。可以证明，转换矩阵T的的逆矩阵等于它的转置矩阵，所以，在整体坐标系下的单元刚度矩逆矩阵等于它的转置矩阵，所以，在整体坐标系下的单元刚度矩阵为阵为eeTkT k T(4.46) b) 进行整体刚度矩阵的组集。可以采用直接刚度法进行整体刚度矩阵的组集。可以采用直接刚度法 。 （6）引入约束条件。）引入约束条件。 （7）求解系统方程，得到所有的节点位移。）求解系统方程，得到所有的节点位移。 （8）进而再求出单元的应力应变等

27、。）进而再求出单元的应力应变等。 式中，式中， 是是x轴相对于轴相对于x轴的夹角。可以证明，转换矩阵轴的夹角。可以证明，转换矩阵T的的逆矩阵等于它的转置矩阵，所以，在整体坐标系下的单元刚度矩逆矩阵等于它的转置矩阵，所以，在整体坐标系下的单元刚度矩阵为阵为eeTkT k T(4.46)例例4-1 平面梁单元应用举例。设一方形截面的悬臂梁，截面每边长平面梁单元应用举例。设一方形截面的悬臂梁，截面每边长为为5cm，长度为，长度为10m，在左端约束固定，在右端施以一个沿，在左端约束固定，在右端施以一个沿y轴负轴负方向的集中力方向的集中力w=100N，求其挠度与转角。，求其挠度与转角。图4-6 平面梁单

28、元实例图利用利用matlab和和ansys两种方法求得的结果基本一致两种方法求得的结果基本一致: 左端点沿左端点沿y方向位移方向位移(挠曲挠曲)：0 左端点绕左端点绕z轴的转角：轴的转角：0 中间点沿中间点沿y方向位移方向位移(挠曲挠曲)：-0.05 中间点绕中间点绕z轴的转角：轴的转角：-0.018 右端节点沿右端节点沿y方向位移方向位移(挠曲挠曲)：-0.16 右端节点绕右端节点绕z轴的转角：轴的转角：-0.024 对于具有两个节点的空间梁单元，设其节点坐标和相应的对于具有两个节点的空间梁单元，设其节点坐标和相应的节点力如下节点力如下 节点（节点（1）：）：11111111111111Te

29、xyzTexyzxyzuvwFNQQMMM (4.47) 节点（节点（2）：）：22222222222222TexyzTexyzxyzuvwFNQQMMM (4.48) 整体坐标系记为整体坐标系记为OXYZ，梁单元的局部坐标系记为，梁单元的局部坐标系记为oxyz，其，其中中ox轴正方向由轴正方向由i端截面形心指向端截面形心指向j端面形心，端面形心，y轴和轴和z轴是梁截面轴是梁截面的两个相互垂直的形心主轴，见图的两个相互垂直的形心主轴，见图4-7。坐标变换公式具有如下。坐标变换公式具有如下形式：形式：1231231233 3lllmmmnnnt12 12ttttT (4.49)图图4-7 空间梁

30、单元的坐标变换空间梁单元的坐标变换 由局部坐标向整体坐标的位移变换公式是由局部坐标向整体坐标的位移变换公式是TeT (4.50) 节点力的变换公式是节点力的变换公式是 TeFT F (4.51) 单元刚度矩阵变换公式是单元刚度矩阵变换公式是TeKT K T (4.52) 在三维空间中，设在三维空间中，设x,y,z是局部坐标系，是局部坐标系，X，Y，Z是整体坐标是整体坐标系系, 单元局部坐标系的三个坐标轴的方向余弦分别如下式：单元局部坐标系的三个坐标轴的方向余弦分别如下式：111cos( ,)cos( , )( , )lxXmxYncos xZ222cos( ,)cos( , )cos( , )

31、lyXmyYnyZ333cos( ,)cos( , )cos( , )lzXmzYnzZ (4.53) 坐标变换矩阵的具体求算方法包括如下步骤。坐标变换矩阵的具体求算方法包括如下步骤。 (1) 局部坐标系局部坐标系x轴在整体坐标系中的方向余弦：轴在整体坐标系中的方向余弦：211222212121211222212121211222212121cos( ,)()()()cos( , )()()()( , )()()()xxlxXxxyyzzyymxYxxyyzzzzncos xZxxyyzz (4.54) （2）局部坐标系）局部坐标系y轴在整体坐标系中的方向余弦轴在整体坐标系中的方向余弦 现在讨

32、论具有任意方向的空间梁单元。首先，由节点现在讨论具有任意方向的空间梁单元。首先，由节点i、j 在整体坐标系下的坐标即可确定在整体坐标系下的坐标即可确定e1在整体坐标系中的三个方在整体坐标系中的三个方向余弦，即向余弦，即LXXij11LYYij12LZZij13 (4.55) 其中其中2/1222ijijijZZYYXXL (4.56) 下面计算下面计算e2和和e3。 在单元的主惯性平面在单元的主惯性平面oxy上任取一点上任取一点k（但（但k点不能取在点不能取在ox轴轴上），点在整体坐标系中的坐标记为（上），点在整体坐标系中的坐标记为（Xk，Yk,Zk）。）。 沿矢量沿矢量 方向取矢量方向取矢量

33、g， g在整体坐标系中的三个分量是在整体坐标系中的三个分量是1g2g3gKikDXX KikDYY KikDZZ (4.57) 2/1222)()()(ikikikkZZYYXXD (4.58) ik 因因z轴垂直于轴垂直于oxy平面，而平面，而e1和和g均在平面均在平面oxy上，故可取上，故可取3egege11 (4.59) e1和和e3既已确定，只需按右手直角坐标系条件确定既已确定，只需按右手直角坐标系条件确定e2即可。即可。因而可取因而可取 e2= e1+e3 (4.60) 由矢量叉乘法则：由矢量叉乘法则：1112131213131113211212131231321123EEEegggEggEggEggg 最后由式最后由式(4.61）有：）有：1121321311133112121/2222123AAAAAAA322gggggg (4.61) 则则3egege1131223AAAAAA1EEE 从而从而