第二部分 平面杆系结构的线弹性有限元法

§3 平面杆系结构的线弹性有限元法§3.1 概论在有限元法中，可以采用位移法，也可以采用力法或混合法。

其中提出最早并且应用最广的是位移法。

对于平面杆系结构来说，位移法实际上就是结构力学中的矩阵位移法（也称刚度法），在计算时以结点位移作为基本未知量。

杆系结构的矩阵分析实际上就是有限元法。

其基本思路是：先把结构离散成有限个数目的单元，然后再考虑某些条件，将这些离散的单元重新组合在一起进行分析计算。

这样使一个复杂的计算问题转化为简单的单元分析和集合问题。

根据这个思路，杆系结构的有限元法可分为两大步骤：（1）单元分析。

研究单元的受力与变形之间的关系；（2）整体分析。

研究如何将这些离散的单元重新组合得到与实际问题相符合的（如边界条件、外界荷载等等）的计算模型—整体刚度方程。

在有限元中，一般采用矩阵形式进行分析求解，因为矩阵运算不仅使公式非常紧骤，而且形式统一，易于编程，适合在电子计算机上进行自动求解。

因此，在有限元法的一般格式中，应尽量采用矩阵形式进行运算。

§3.2 局部坐标系下的单元刚度矩阵1 单元的划分。

在杆系结构的有限元法中，一般将由相同材料、具有相同横截面的一根杆件（即等截面直杆）当成一个单元，整个结构就是由有限个杆件单元组成的集合体。

杆件单元具有2个结点，即首结点和末结点，但一般是先确定结点的位置，结点一旦确定，则结点之间的单元也就确定了。

在进行杆系结构的单元划分时，应注意如下事项：○1结点位置的确定。

结点一般选在杆件的如下位置：杆件的转折点、杆件汇交点、支承点、截面或材料的突变点，这些点都是结构的构造点，有时为了使结构只承受结点荷载，在集中荷载的作用处也设置一个结点。

○2结点的编号。

为了使集合以后的总刚的带宽最小，一般应遵循尽量使相关结点（有单元相连的结点）编号差值的最大值最小的原则进行。

2 单元刚度矩阵考虑一等截面的平面梁单元，单元首末结点分别为j i ,，单元长为l ，单元抗弯刚度为EI ，E 为材料的弹性模量，I 是截面的抗弯惯矩，取x 轴为沿梁单元中心轴，y 轴与x 轴成90o，如图1所示。

○1位移模式和形函数 如果不考虑杆件的轴向变形与横向弯曲变形的互相影响，且设x 轴向的位移u （即单元轴向位移）取为x 的线性函数，而对于y 轴向的位移v （即单元横向位移，亦即梁单元挠度）取为x 的三次多项式表示。

于是有：⎩⎨⎧+++=+=332211xb x b x b b v xa a u o o (1) 不考虑剪切变形的影响，即假定0=γ，可推导出位移模式又可表示为{}{}ej j j i i i v vvvu uN v u v u N N N N N N v u r δθθ][00000432111=⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧= （2） 其中{}TTj Tie][δδδ=，其中{}{}Tj j jjTi iii v uv u ][,][θδθδ==；lx N lx N u u =-=21,1232433223232233221232231lx lx Nl x lx N l xl x x Nl x lx N vv vv+-=-=+-=+-=式中1u N 、2u N 、1v N 、2v N 、3v N 、4v N 反映了单元的位移形态，称为形态函数，简称形函数，][N 称为形函数矩阵。

○2用结点位移表示的应变和应力 如单元产生拉压变形和弯曲变形，则其纵向纤维的线应变可分成两个部分： （1）0ε－称为拉压应变，也称为轴向应变，整个梁截面都相同；（2）b ε－弯曲应变，沿梁高y 因存在曲率κ（称为广义弯曲应变）而不同。

假设不考虑剪切变形，由梁单元的几何方程可知⎪⎪⎩⎪⎪⎨⎧-=-==曲应变） （曲率，广义弯（轴向应变）220)( dx v d dx dv dx d dx du κε （3）在杆件截面高度y h =处，有曲率κ引起的应变为22)(dxv d yy y b -=⋅=κε （4）由（3）和（4）式可写成 {}{}e v v v v u u b x yN x yN x yN x yN x N x N dx v d y dx du δεεε⎥⎦⎤⎢⎣⎡''-''-''-''-''=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧-=⎭⎬⎫⎩⎨⎧=)()(0)()(00)(00)(432121220（5）式中)(1'x N u表示)(1x N u 对x 求一阶导，)(1''x N v 表示)(1''x N v 对x 求二阶导，其余说明类似。

上式可写成{}{}e B δε][= （6）由虎克定律，就可以得到利用结点位移表示的应力表达式{}{}εεεεεσσσ][][000D D E E b b b =⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧=（7） 式中][D 称为弹性矩阵，且⎥⎦⎤⎢⎣⎡⋅=1001][E D 由式（6）可得{}{}{}{}eeS B D D δδεσ][][][][=⋅⋅== （8）式中][][][B D S ⋅=，称为应力矩阵。

○3由虚位移原理导出梁单元的刚度方程 应用虚功原理于梁单元上，可得到梁单元保持平衡的刚度方程为{}{}[]{}[][][]{}eVT lTed e dv B D B dx q N F F δ⎰⎰⎰⎰=++ （9）令{}{}{}{}{}ed eq ed lTeP F F F dx q N F +=+=⎰][ （10）⎰⎰⎰Ω=υd B D B k Te]][[][][ （11）式中e k ][称为单元刚度矩阵，eq F ][为由于分布荷载而移置的等效结点荷载，ed F ][为直接作用在结点上的荷载。

于是{}{}{}eeeP ek F F δ][=+ （12）令{}{}{}eP eeeqF F F += （13）则{}{}eeeeqk F δ][= （14）上式即为局部坐标系中的单元刚度方程。

将][B 的具体的表达式代入（12）式，经积分和矩阵运算可得到平面梁单元的单元刚度矩阵为e k ][,具体表达式如下:⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=666563625655535244413635333226252322141100000000000000][k k k k k k k k k k k k k k k k k k k k k e（15） 式(13)中元素ij k 为平面杆系单元刚度矩阵元素，即：显然单元刚度矩阵ek ][为对称矩阵。

⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧=====-=-=-=====-=-===-=-==lEI k k l EI k k lEI k k k k k k k lEI k k k k l EA k k k k 246126336663326556356226322335225552241144411§3.3 整体坐标系下的单元刚度矩阵在前面的分析之中，单元刚度矩阵是在单元的局部坐标系y o x 中形成的，由于各个单元的局部坐标系不同，因此必须将每个单元的单刚转换到同一个公共的坐标系下，这个公共坐标系就是整体坐标系xoy 。

为了区别起见，在局部坐标系下的杆端分量符号顶上加“－”。

下面首先介绍转换矩阵的概念，然后据之建立整体坐标系下的单元刚度矩阵。

一 转换矩阵如图3－2所示，任一单元○e 的首端结点力在两种坐标系中的分量。

其中图(a )表示在局部坐标系y o x 中的三个分量i X 、i Y 和i M ，图(b )表示在整体坐标系xoy 中的三个分量i X 、i Y 和i M 。

为了导出i X 、i Y 、i M 与i X 、i Y 、i M 之间的关系式，在图中将两个力i X 、i Y 分别投影在x 和y 轴上，可得出下式中的前两式：⎪⎪⎭⎪⎪⎬⎫=+-=+=i i i i i i i iM M Y X Y Y X Xααααcos sin sin cos其中第三式表明，在两个坐标系中的力偶仅是仍彼此相等。

α表示由x 轴转到x 轴的角度，以逆时针方向为正。

同理对单元○e 的另一端力也可得出类似的关系：⎪⎪⎭⎪⎪⎬⎫=+-=+=j j j j j j j jM M Y X Y Y X Xααααcos sin sin cos把以上两个方程组组合成一个矩阵方程得⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--=⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧j j j i i i j j j i i i M Y X M Y X M Y X M Y X 1000cos sin 0000sin cos 0000001000000cos sin 0000sin cos αααααααα 或简写成eeF T F }]{[}{=eeeF T F }{][}{=式中{}eF是局部坐标系中的单元杆端力列阵，{}eF 是整体坐标系中的单元杆端力列阵，[T ]称为单元的坐标转换矩阵。

可以证明，[T ]矩阵是正交矩阵，其逆矩阵等于它的转置矩阵，即1][][-=T T T对于单元杆端结点位移，也可以同样进行转换，即有eeT }]{[}{δδ= eeeT }{][}{δδ=二 整体坐标系中的单元刚度矩阵 1 推导过程。

现在研究整体坐标系中的单元等效杆端力eeq F }{与单元结点位移e}{δ之间的关系式。

因为{}{}ee eeeeqT k k F }]{[][][δδ==eeq eeq F T F }]{[}{=所以有{}{}{}ee T eeq Teeq eeqT k T F T F T T F }]{[][][][][][1δ===-令][][][][T k T k eTe=则{}ee eeqk F }{][δ=上式即为整体坐标系下的单元刚度方程，简称单元整体刚度方程。

因此e k ][相应地称为单元整体刚度矩阵。

2 单元刚度矩阵的特点。

（1）由][][][][)]([][)]([])[]([])[][]([)]([T k T T k T T T k T k T k e T T e T T T T e T e T T e ====，可知单元刚度矩阵e k ][为对称矩阵。

（2）单元刚度的分块由于在以后的整体分析之中，是对结构的每个结点建立平衡方程，为了以后讨论方便，可把上式按单元的首未结点编号i 、j 进行分块，写成如下形式：式中：分别为首端i 和末端j 的杆端力和杆端位移列阵。

ejj eji eij eii k k k k 、、、为单元刚度矩阵e k ][的四个子块，每个子块都是33⨯阶的方阵，并且有eji eij k k ＝，两个子块的元素相同，ejj eii k k 、子块为对称子块。