第十一章有限元分析法概述

1、第十一章有限元分析方法概述1、基本概念有限元分析方法是随着电子计算机的发展而迅速发展起来的一种现代没计计算方法。它是 20 世纪 50 年代首先在连续体力学领域飞机结构静、动态特性分析中应用的一种有效的数值分析方法，随后很快就广泛地应用于求解热传导、电磁场、流体力学等连续性问题。在工程分析和科学研究中，常常会遇到大量的由常微分方程、偏微分方程及相应的边界条件描述的场问题，如位移场、应力场和温度场等问题。求解这类场问题的方法主要有两种：用解析法求得精确解；用数值解法求其近似解。应该指出，能用解析法求出精确解的只是方程性质比较简单且几何边界相当规则的少数问题。而对于绝大多数问题，则很少能得出解析解

2、。这就需要研究它的数值解法，以求出近似解。目前工程中实用的数值解法主要有三种：有限差分法、有限元法和边界元法。其中，以有限元法通用性最好，解题效率高，目前在工程中的应用最为广泛。下面通过一个具体例子，分别采用解析法和数值解法进行求解，从而体会一下有限元分析方法的含义及其相关的一些基本概念。如下图所示为一变横截面杆，杆的一端固定，另一端承受负荷P ，试求杆沿长度方向任一截面的变形大小。其中，杆的上边宽度为w1 ，下边宽度为 w2 ，厚度为 t ，长度为 L ，杆的材料弹性模量为 。已知4450N ，12tLEEPw 50mm，w=25mm ， =3mm ， =250mm ，=72GPa。 采用解

3、析法精确求解假设杆任一横截面面积为A( y) ，其上平均应力为，应变为。根据静力平衡条件有：P A( y)0根据虎克定律有：E而任一横截面面积为：A( y) (w1w2w1 y)tLdu任一横截面产生的应变为：dy将上述方程代入静力平衡条件，进行变换后有：PdudyEA( y)沿杆的长度方向对上式两边进行积分，可得：uyPyPdydu0 EA ( y )dyw2w100y )Et ( w1L1将 A( y) 表达式代入上式，并对两边进行积分，得杆沿长度方向任一横截面的变形量：u( y)PLln( w1w2 w1 y ) ln w1 Et ( w2 w1 )L当 y 分别取 0、62.5 、12

4、5、 187.5 、250 值时，变截面杆相应横截面处的沿杆长方向的变形量分别为：u0m; u227.51 106 m; u 59.27 10 6 m; u496.83 10 6 m; u5142.80 10 6 m13 采用数值解法近似求解将变横截面杆沿长度方向分成独立的 4 小段，每一小段采用等截面直杆近似，等截面直杆的横截面面积为相应的变截面杆横截面面积的平均面积表示，每一小段称为一个单元，小段之间通过节点连接起来。这样，变横截面杆就可以用 5 个节点和 4 个单元组成的模型来近似表示，如右图所示。假设任一横截面面积为 A、长为 l 的等截面直杆，在轴向拉力 F 的作用下产生变形量 l

5、，则该直杆横截面上的应力和应变分别为：FlAl根据虎克定律：E可得：AEFll上述方程与线性弹簧的方程 F kx 极为相似，表明一个中心点集中受力且横截面相等的等截面直杆可以等效为一个弹簧，其等效刚度为 :AEkeql因此，变横截面杆可以看作由四个线性弹簧串联起来的模型来近似表示，如下图所示，每一个单元都可以视为一个线性弹簧，其弹性行为符合以下方程：f keq(ui 1 ui )Aavg E( Ai 1Ai ) Eui )(ui 1 ui )(ui 1l2l下面考虑每一个节点的受力，根据静力平衡条件，每一个节点上的受力总和为0，即：节点 1： R1k1 (u2u1 )0节点 2：k1 (u2u

6、1 )k2 (u3u2 )0节点 3： k2 (u3u2 )k3 (u4u3 )0节点 4： k3 (u4u3 )k4 (u5u4 )0节点 5： k4 (u5u4 )P02将反作用力 R1和外力 P 从内力中分离出来，重新对上述五个方程组成的方程组进行变换，得：节点 1： k1u1k1u2R1节点 2： k1u1(k1k2 )u2k2u30节点 3：k2u2(k2k3 )u3k3u40节点 4：k3u3(k3k4 )u4k4u50节点 5：k4 u4k5u5P将上述方程组写成矩阵形式，有：k1k1000u1R1k1k1 k2k200u200k2k2 k3k30u3000k3k3 k4k4u4

7、0000k4k4u5P将反作用力和外力分离出来，可以重组上述矩阵，得：R1k1k1000u100k1k1 k2k200u2000k2k2 k3k30u30000k3k3 k4k4u400000k4k4u5P写成一般形式，可得：RK U F即表示：反作用力矩阵总体刚度矩阵 位移矩阵 负荷矩阵 引入边界条件，根据本题的要求，节点1 的位移为 0，即u10 ，则有如下矩阵形式：10000u10k1 k1k2k200u200k2k2 k3k30u3000k3k3 k4k4u40000k4k4u5P通过求解上述矩阵方程，可得每个节点的位移，进而可以求得每个节点的反作用力，每一个单元的应力和应变。即：ui

8、 1uiE iiil根据变横截面杆结构的已知参数可得：lL62.5mmA( y) (w1w2w1 y)t 1500.3 y4L3当 y 0 时， A1 150mm2当 yL 时， A2 1500.3L1500.3 250131.25mm2444当 y2L 时， A31500.32L1500.32250112.5mm2444当 y3L时， A41500.33L1500.3325093.75mm2444当 yL时， A51500.3L1500.3 25075mm2每个单元的等效刚度系数keq( Ai1Ai ) E2lk1( A1A2)E(150131.25)72106162106N /M2l262

9、.5k2( A2A3 )E(131.25112.5) 7210 6140.4106N/M2l262.5k3( A3A4 )E(112.5 93.75)72106118.8106N/M2l262.5k4( A4A5 )E(93.7575) 7210697.2 106 N / M2l262.5总体刚度矩阵：162162000162162140.4140.400K(G)1060140.4140.4118.8118.8000118.8118.897.297.200097.297.2应用边界条件 u10和负荷 P4450N，可以得到：10000u10162162140.4140.400u20106014

10、0.4140.4118.8118.80u3000118.8118.897.297.2u4000097.297.2u54450求解该方程，可得：u0m;u227.47 10 6 m;u359.1610 6 m;u496.6210 6 m;u142.4010 6 m15而第一种精确求解方法求得的每个节点处的位移分别为：u10m;u227.51 10 6 m;u359.2710 6 m;u496.8310 6 m;u5142.8010 6 m比较两种结果表明：采用数值解法近似求解的结果与解析法精确求解的结果相当接近，如果将变横截面杆沿杆长方向分离成的单元越多，数值解法求解的结果将与精确解法求得的结果

11、误差将会越来越小。42、有限元法的分析过程2.1连续体离散化所谓连续体是指所求解的对象 （物体或结构），所谓离散化就是将所求解的对象划分为有限个具有规则形状的微小块体，每个微小块体称为单元，两相邻单元之间只通过若干点互相连接，每个连接点称为节点。因而，相邻单元只在节点处连接，载荷也只通过节点在各单元之间传递，用这些有限个单元构成的集合体来近似代替原来的连续体。这种由单元和节点构成的集合体称为有限元分析模型。离散化也称为划分网格或网格化。单元划分后，给每个单元及节点进行合理编号；选定坐标系，计算各个节点坐标；确定各个单元的形态和性态参数以及边界条件等。下图所示为将一悬臂梁建立有限元分析模型的例子

12、，图中将该悬臂梁划分为许多三角形单元，三角形单元的三个顶点都是节点。结构离散化后，单元与单元之间利用单元的节点相互连结起来，单元节点的设置、性质、数目等应视具体问题的性质、描述变形形态的需要和计算精度而定。所以有限元法中分析的结构已不是原有的物体或结构物，而是同样材料的由众多单元以一定方式连结成的离散物体。这样，用有限元分析计算所获得的结果只是近似的。如果划分单元数目非常多而又合理，则所获得的结果就与实际情况相接近。2.2单元特性分析连续体离散化后，即可对单元体进行特性分析，简称为单元特性分析。单元特性分析主要有两项：选择单元位移模式（位移函数）和分析单元的特性，即建立单元刚度矩阵。根据材料学

13、、工程力学原理可知，弹性连续体在载荷或其他因素作用下产生的应力、应变和位移， 都可以用位置函数来表示。 为了能用节点位移来表示单元体内任一点的位移、应变和应力，必须对各单元中位移的分布作出某种假设，也就是假定单元中任一点的位移是单元节点位移的某种简单的函数，以此模拟单元内位移的分布规律，这种函数就称为位移模式或位移函数。选择适当的位移函数是有限单元法分析与计算中的关键，通常采用多项式作为位移模式。因为多项式的数学运算比较方便，并且所有光滑函数的局部都可以用多项式逼近。至于多项式的项数和阶次的选择，则要考虑到单元的自由度和解的收敛性。一般来说，多项式的项数应等于单元的自由度数 （单元节点独立位移

14、的个数） ，多项式的阶次应包含常数项和线性项等。5选定好单元位移模式后，即可进行单元力学特性分析，将作用在单元上的所有力（表面力、体积力、集中力）等效地移置为节点载荷。根据单元的材料性质、形状、尺寸、节点数目、位置及其含义等，应用有关的力学原理建立单元内节点力与节点位移之间的方程式，从而导出单元刚度矩阵。2.3整体分析在对全部单元进行完单元分析之后，就要进行单元组集，即把各个单元的刚度矩阵集成为总体刚度矩阵，以及将各单元的节点力向量集成总的力向量，求得整体平衡方程。集成过程所依据的原理是节点变形协调条件和平衡条件。2.4确定约束条件由上述所形成的整体平衡方程是一组线性代数方程， 在求解之前，

15、必需根据具体情况，分析与确定求解对象问题的边界约束条件，并对这些方程进行适当修正。2.5有限元方程求解解方程，即可求得各节点的位移，进而根据位移计算单元的应力及应变。3、有限元法的理论基础有限元法是一种离散化的数值解法，对于结构力学特性的分析而言，其理论基础是能量原理。根据未知数的性质和分析方法的不同，有三种基本解法：1) 位移法。位移法采用最小势能原理或虚位移原理进行分析。 它以节点位移作为基本未知量，选择适当的位移函数，进行单元的力学特性分析，在节点处建立单元刚度方程，再合并组成整体刚度矩阵，解出节点位移后，由节点位移再求解出应力。位移法优点是比较简单，规律性强，易于编写计算机程序。所以得

16、到广泛应用，其缺点是精度稍低。2) 应力法。应力法常采用最小余能原理进行分析。它以节点力作为基本未知量，在节点处建立位移连续方程，求解出节点力后，再求解节点位移和单元应力。力法的特点是计算精度高。3) 混合法。这种方法常采用混合变分原理进行分析。 它取一部分节点位移和一部分节点力作为基本未知量，建立平衡方程进行求解。在进行结构静力学分析中，对大多数问题，位移法要比应力法简单得多，从而得到了最广泛的应用和发展，在本书中只讨论有限元位移法。4、有限元分析软件目前有限元分析软件可以分为三类：1）通用有限元分析软件：这类软件自成体系，侧重点有所不同，但解决工程问题的领域比较宽，适应性和通用性强，比较有

17、代表性的有：ABAQUS 、ADINA 、ANSYS 、MARC和 NASTRAN 等2）专用有限元分析软件：主要特点是在某一专门领域内，开发了专门的功能，强调专用性。比较有代表性的有： ADAMS 、DADS 、MSC/FATIGUE 等等。3）嵌套在 CAD/CAM 系统中的有限元分析模块： 这类分析模块与设计软件集成为一体，6有限元分析在工程师所熟悉的设计环境中进行，功能没有专用或通用有限元分析软件那么强大全面，但它们解决一般工程问题的能力也是很强的，比较有代表性的有：I-DEAS 、PRO/ENGINEER 和 UNIGRAPHICS 等 CAD/CAM/CAE 系统中的有限元分析模块

18、。采用有限元分析软件进行结构分析的基本步骤包括：1）预处理阶段，分析对象的有限元网格剖分与数据生成。主要包括以下几个方面：建立求解域并将之离散化成有限元，即将问题分解成节点和单元；假设代表单元物理行为的形函数，即假设代表单元解的近似连续函数；对单元建立方程；将单元组合成总体的问题，构造总体刚度矩阵；应用边界条件、初值条件和负荷。2）求解阶段。求解线性和非线性的微分方程组，得到节点的值。3）后处理阶段。根据工程和产品模型与设计要求，对有限元分析结果进行用户所要求的加工和检查，并已图形方式将结果提供给用户，辅助用户判定计算结果与设计方案的合理性。具体包括：有限元分析结果的数据平滑，各种物理量的加工

19、和检查，如：结构变形图、应力分布图和结构动力振型图等，针对工程和产品设计的要求与工程规范对结果进行校核，根据计算结果进行设计优化与模型修改，还包括计算结果的文档整理等。目前，应用比较广泛的有限元分析软件主要是ANSYS软件。ANSYS软件是融结构、传热、流体、电场、磁场、声场分析于一体的大型通用有限元分析软件，由美国有限元分析软件公司 ANSYS开发，并能与多数 CAD 软件接口，实现数据的共享和交换，如 Pro/Engineer ， NASTRAN，Alogor ，I-DEAS，AutoCAD 等，是现代产品设计中的高级 CAD 工具之一。软件主要包括 3 个部分：前处理模块： 提供强大的实

20、体建模及网格划分工具， 用户可以方便地构造有限元模型；分析计算模块：包括结构分析 ( 可进行线性分析、非线性分析和高度非线性分析 ) 、流体动力学分析、电磁场分析、声场分析、压电分析以及多物理场的耦合分析，可模拟多种物理介质的相互作用，具有灵敏度分析及优化分析能力；后处理模块：可将计算结果以彩色等值线显示、梯度显示、矢量显示、粒子流迹显示、立体切片显示、透明及半透明显示 ( 可看到结构内部 ) 等图形方式显示出来，也可将计算结果以图表、曲线形式显示或输出。软件提供了 100 种以上的单元类型，用来模拟工程中的各种结构和材料。分析类型主要包括以下几种：1) 结构静力分析。用来求解外载荷引起的位移

21、、应力和力。静力分析很适合求解惯性和阻尼对结构的影响并不显著的问题。 ANSYS程序中的静力分析不仅可以进行线性分析，而且也可以进行非线性分析，如塑性、蠕变、膨胀、大变形、大应变及接触分析。2) 结构动力学分析。结构动力学分析用来求解随时间变化的载荷对结构或部件的影响。与静力分析不同，动力分析要考虑随时间变化的力载荷以及它对阻尼和惯性的影响。ANSYS可进行的结构动力学分析类型包括：瞬态动力学分析、模态分析、谐波响应分析及随机振动响应分析。3) 结构非线性分析。结构非线性导致结构或部件的响应随外载荷不成比例变化。ANSYS7程序可求解静态和瞬态非线性问题，包括材料非线性、几何非线性和单元非线性

22、3 种。4) 动力学分析。 ANSYS程序可以分析大型三维柔体运动。 当运动的积累影响起主要作用时，可使用这些功能分析复杂结构在空间中的运动特性， 并确定结构中由此产生的应力、应变和变形。5) 热分析。程序可处理热传递的 3 种基本类型：传导、对流和辐射。热传递的 3 种类型均可进行稳态和瞬态、线性和非线性分析。热分析还具有可以模拟材料固化和熔解过程的相变分析能力以及模拟热与结构应力之间的热 - 结构耦合分析能力。6) 电磁场分析。主要用于电磁场问题的分析，如电感、电容、磁通量密度、涡流、电场分布、磁力线分布、力、运动效应、电路和能量损失等。还可用于螺线管、调节器、发电机、变换器、磁体、加速器

23、、电解槽及无损检测装置等的设计和分析领域。7) 流体动力学分析。 ANSYS流体单元能进行流体动力学分析， 分析类型可以为瞬态或稳态。分析结果可以是每个节点的压力和通过每个单元的流率。并且可以利用后处理功能产生压力、流率和温度分布的图形显示。另外，还可以使用三维表面效应单元和热流管单元模拟结构的流体绕流并包括对流换热效应。8) 声场分析。程序的声学功能用来研究在含有流体的介质中声波的传播， 或分析浸在流体中的固体结构的动态特性。这些功能可用来确定音响话筒的频率响应，研究音乐大厅的声场强度分布，或预测水对振动船体的阻尼效应。9) 压电分析。用于分析二维或三维结构对交流 (AC)、直流 (DC)或

24、任意随时间变化的电流或机械载荷的响应。这种分析类型可用于换热器、振荡器、谐振器、麦克风等部件及其他电子设备的结构动态性能分析。可进行 4 种类型的分析：静态分析、模态分析、谐波响应分析、瞬态响应分析。下面主要讨论一下总体刚度矩阵与单元刚度矩阵的关系我们取任一等截面直杆单元，用线性弹簧近似表示，如图所示。在节点处传递的力为：fiki (uiui 1 )fi 1ki (ui 1ui )将其表示成矩阵形式：f ikikiuifi 1kikiui 1则单元刚度矩阵为： K (i )kikikiki总体刚度矩阵为： K ( G ) K (1) K ( 2) K (3). K ( n )对于本题，其总体刚

25、度矩阵为：8k1k1000k1 k1k2k200 K (G)0k2k2 k3k3000k3k3 k4k4000k4k4下面采用最小势能原理来推导总体刚度矩阵物体受到外力的作用会产生变形，外力所做的功以弹性能的形式储存在物体中，即物体储存应变能。假设：有一等截面杆，长L ，在集中力 F 的作用下变形伸长 L根据虎克定律，有： F(AE) LkyL杆内储存的能量： dyy1 k y ) yFd yky dy (002将上式写成标准应力和应变的形式：d( 1y dxdz )dy1y y dV22则变形能为：( e )d(1E2y y) dVVV22y dV由 n 个单元和 m 个节点组成的物体的总势

26、能II 为总应变能和外力所做功的差:IInm( e )F i uie1i1根据最小总势能原理，对于一个稳定的系统，平衡位置会发生位移，并使系统的总势能最小。即：IIn( e)mF i u i0u iu i e 1u i i 1对于例 1，任意单元（ e）的应变能由应变能公式可得：( e )Aavg E(u i2u i22u i 1u i )21( e )A avg E则：u i1 )(u iu i( e )AavgE1u i )( uiu i1写成矩阵形式：9( e )u ik ek eui( e )k ek eu i1u i1而：( Fi ui ) Fi( Fi 1ui 1)Fi 1uiui

27、 1将总势能最小公式写成矩阵形式，得出k1k1000u1F1k1k1 k2k200u2F20k2k2 k3k30u3F300k3k3k4k4u4F4000k4k4u5F5根据力边界条件，可知：F1=-RF2=0F3=0F4=0F5=P则得：k1k1000u1Rk1 k1k2k200u200k2k2 k3k30u3000k3k3 k4k4u40000k4k4u5P根据位移边界条件： u1=0则有：10000u10k1 k1k2k200u200k2k2 k3k30u3000k3k3 k4k4u40000k4k4u5P求解上式矩阵，得出的结果与采用直接公式法得出的结果相同。10第十二章杆梁结构的有限

28、元分析第一节弹簧单元任取一弹簧单元，如图所示。两个节点 ij ，节点位移 ui u j ，节点力 fi f j，弹簧刚度 k ， 根据节点处力的平衡可知：fi k(uiu j )f jk(u jui )将其表示成矩阵形式：fikkuif jkku j则定义单元刚度矩阵为： K kkk。 很显然，单元刚度矩阵为一对称矩阵。k如果将两个弹簧串联组成一个弹簧系统，如下右图所示，则系统有两个单元，三个节点。对于单元 1，有矩阵方程：f11k1k1u1f 21k1k1u2对于单元 2，也有矩阵方程：f12k2k2u2f 22k2k2u3式中， fi m 表示第 m 单元的第 i 节点上作用的内力。根据节

29、点上作用力平衡条件，有：节点 1：1节点1f12F2节点：2F3f1 F12： f 23f2也就是下面方程组：k1u1k1u2F1k1u1( k1k2 )u2k2u3F2k2 u2k2u3F3写成矩阵形式，有：k1k10u1F1k1k1k2k2u2F20k2k2u3F3则串连弹簧系统的总体刚度矩阵为：11k1k10Kk1k1 k2k20k2k2u1F1位移矩阵为： Uu2节点负荷矩阵为： FF2则矩阵方程可简化为：u3F3KUF假设：位移边界条件为 u10 ； 力边界条件为 F2F3P ，则矩阵方程变为：k1k100F1k1 k1k2k2u2P0k2k2u3P将其进一步简化为： F1k1u2k

30、1k2k2u2Pk2k2u3Pu22P / k1F12P求解可得：2P / k12P / k2u3例 1：如图所示一弹簧系统，已知k1100N/mm ， k2200N/mm ， k3 100N/mm ，P 500N ， u1 u40 。 求： 1）系统的整体刚度矩阵；2）第 2 与第 3 节点的位移； 3）第1 与第 4 节点的反力； 4）中间弹簧的受力大小。解：根据上面可知，系统总体刚度矩阵为：k1k10010010000k1k1 k2k20100 3002000Kk2k2 k3k30200300100000k3k30010010010010000u1F1系统矩阵方程为：1003002000

31、u200200300100u3P00100100u4F4100100000F1施加边界条件 u11003002000u20u4 0 可得：0200300100u3P001001000F412300200u20将矩阵方程简化可得：300u3P200u2P/2502求解上述矩阵方程得：3P / 500(mm)u33则节点 1 与节点 4 的反力： F1100u2200NF4100u3300N200200u2f2单元 2 的矩阵方程为：200u3f3200则单元 2 的弹簧力为： F f3f2200(u3 u2 )200 N例 2：下图所示一弹簧系统，写出其整体刚度矩阵和边界条件。解： K1k1k1

32、K 2k2k2k1k1k2k2K 3k3k3K 4k4k4k3k3k4k4整体刚度矩阵为：k4k4000k4 k1k2 k4k2k10K0k2k2 k30k30k10k1000k30k3系统刚度方程为：k4k4000u1F1k4k1 k2 k4k2k10u2F2K0k2k2 k30k3u3F30k10k10u4F400k30k3u5F5边界条件为： uu0FPFP451132课后习题：考虑如图所示的弹簧单元系统，假定 k 120KN/m ， P 20KN 。求：1）结构的整体刚度矩阵； 2 ）节点 3 、4、5 的位移； 3）节点 1 、2 的支反力； 4 ）每个弹簧的内力。13第二节线性等截

33、面直杆单元1.1单元特性分析右下图所示为一线性等截面直杆单元，其中，杆的长度为L ，横截面积为 A ，材料弹性模量为 E 。根据材料力学可知，单元在节点轴向力的作用下，杆内应力和应变在轴线各点处均是恒定常数，因而杆内任一点位移沿杆轴线呈线性规律变化。建立杆单元的整体坐标系， 则杆单元的位移函数可以表示为：ua1a2 x写成矩阵形式：u1 xa1a2当 x xi 时 u ui；当 x x j时 u u j代入位移函数可得：uia1a2 xiu ja1a2 x jui1xia1写成矩阵形式：1x ja2u j1xi1a11xi1ui可得：将上式两边左乘x ja21x ju j1则杆单元中的任一点位

34、移可表示为：1xi1u1xui1x ju j简化可得：uxjxxxiuix jxix jxiujx jxxxi称为形函数，则 N N iN j 称为形函数矩阵，它将单元令： NiN jxixjxix j的节点位移与单元的内位移连接起来。则：uNU建立杆单元的局部坐标，则有：xi0x jL形函数矩阵可变为： NN iN j1xxLL14则位移函数为： u1xxui1 xuixu ju jLLLL令：x01， 称为自然坐标。则形函数可表示为：LN i( ) 1， N j ( )形函数矩阵变为 :N Ni ()N j () ， 则单元中任一点位移为：u( x) N i ( )uiN j ( )u j

35、于是，应变为 :dudNUBUdxdx其中： BdN称为单元应变矩阵，则：dxBdNd N i ( ) N j ( )d N i( ) N j ( ) d1 1dxdxddxL L因此，应力可表示为：EEBU单元应变能：U1 TdV1(U T BT EBU )dV1 U T(BT EB) dV U2 V2 V2V单元的外力势能：W( fi uif j u j )U T f则单元的总势能为：IIUW1 U T(BT EB)dV UU T f2V根据最小势能原理可知：II0 ，则有：U( BT EB )dV UfV令 K( BT EB)dV ，称为杆单元在局部坐标系中的单元刚度矩阵为。则杆单元刚度方程为：VK Uf将单元应变矩阵表达式代入单元刚度矩阵可得 :K(BTLEB) dVV01/ LEA11E1/ L 1/L Adx111/ LL这就是等截面直杆在局部坐标系中的单元刚度矩阵。15例 1：如图所示，求节点1 和节点 3 的反力。解：这是由两个等截面直杆