1、 nnnxaxaxaaxPY2210)(nmiiniyxPD)(),()()(21012101212nmiininimiiinmiiaaaaFyxaxaayxPD 0)(1110jimiininjxyxaxaaaFmijiiminjinmijimijixyxaxaxa1111110,;11kkimiikmikiTxySx), 2 , 1 , 0(0njTSajnijii 01100TaSaSaSnnnnnnnnTaSaSaSaS222110 32dNcNbNaT 床身及导轨进给箱刀架车床属座溜板箱主轴箱离合器主轴组件 中间变速机构主轴主轴齿轮主轴轴承 naaaa321,),(321naaaa
2、461235 )(MPab)(MPas)(1MPa)(1MPa wnslnsyyysxx规格长像素宽像素规格 yxo(x,y)xyo),(nnyxxyo),(ssyxa用户坐标系b规格化坐标系c设备坐标系 三种坐标系的关系 EF0000GKLK0110CD1010100100010101IJIHG 0101010144max33min22max11minccyyccyyccxxccxx否则当否则当否则当否则当1234,cccc sfeqdcpbaAecyaxufdybxvsqypxw xTyTXY1010001yxTT平移变换矩阵 XY比例平移向比例系数向比例系数ySxSyx1000000yx
3、SS XY0100010001关于坐标轴反射对称轴为Y轴变换矩阵 XY0关于原点反射变换矩阵100010001 XY0绕原点旋转变换矩阵1000cossin0sincos XY0变换矩阵错切变换10001001a MNABCyxyx11),(00yx 101000100yxA1000cossin0sincosB 101000100yxC1sin)cos1 (sin)cos1 (1cossin0sincos0000 xyyxABCT 1p2p3p6p14p16p13p15p 全，边界，初始条件不健材料的非线性和不均匀数值）得到通解，如弹性力学微分方程（边界，初值体连续体中无限小的微分解析 ）（边
4、界元法边界有限差分有限元区域BEMFDMFEM Lhyxyzb YXYZh yzxyzxyyx PZOYX yzxzz, xyyx0zz0z 1MM1MM1),(MzyxMZYXyxwvM1Mu ZYXMABCDdxdzdy111.CBA dyyvvdydxxvvu+1B2BBAC1C2CV1A12XYudxxuuUBdxxvvVBdyyuuUCdyyvvVC 2211BABA211BAB dxxuudxxu为：xxxudxdxxuABABB11AyyxvdyvdyyvVACACCA)(11 A1B1-AB= A1B2-AB =uB-uA=( )-u= 即线段AB的正应变= 同理，线段AC的正
5、应变为： = zwzxy12同样的方法研究了微分体在XOZ或YOZ坐标面上的投影 ABCD内，AB与AC所夹直角的变化即投影面的剪应变两部分：一部分与X轴平行 AB向Y的转角一部分与Y轴平行 AC向X的转角xuxvudxxuudxvdxxvvBABBtg1)()(211211 1xxu1xv 故上式可简化为： =2同理: =yu xy12xvyuyuxvyvxuxyyx +=+ 平面应变几何方程，位移应变 当物体的位移分量确定时，应变分量完全确定，反过来，当应变分量确定时，位移分量不完全确定。 xyxyxuyvxvyu说明：以ABCD为例令其应变分量为零，即=0 =0 =0 即:=0 =0 +
6、=0 （C） （1）将（1）式分别对X和Y积分得：u=f1(y) , v=f2(x)代入（C）得:dxxdfdyydf)()(21 dyydf)(1dxxdf)(2)(1yf)(2xf左边是Y的任意函数，右边是X的函数，故两边等于同一常数w，故：= w = w=-wy+u0 =wx+v0 U=u0-wy v=v0+wx平面内的应变全为0时，面内的各点位移，与变形无关，是刚体位移。U0,v0代表沿x轴及y轴的刚体平移，w代表绕z轴的转动（证明略）。 oyzxABCPZZXZYYYXYZzxy（1）空间问题每个面上有三个应力分量：一个正应力，两个剪应力， 垂直作用于表示，前表示垂直于垂直于Z轴的面
7、上沿Z向作用。剪应力用垂直于哪一个轴，后表示沿哪个坐标轴，X轴沿Y轴方向移动。2, ,物理方程应变与应力之间的关系。 zyx线应变zxyzxy角应变（剪应变）x,y两线段间的角应变如果某一个截面上的外法线方向是沿坐标轴的正方向，这个截面称一个正面，这个面上的应力沿正向为正，负方向为负。相反如果某截面上的外法线是沿坐标轴的负方向，截面为负面，这个面上的应力以沿坐标轴负向为正，正向为负。空间问题有九个应力分量。三个正应力六个剪应力三个独立应力。xy=yx yz=zy zx=xz空间应力状态有6个独立应力分量:对应6个应变分量： xxEE为拉伸弹性模量zxzxyzyzxyxyyxzzzxyyzyxx
8、GGGuEuEuE111)(1)(1)(1)1 (uEG完全弹性,各向同性物体,应变与应力关系：（胡克定律导出）胡克定律：在单向应力状态下，处于弹性阶段的物体应力与应变是线性关系，即：（1）E弹性模量剪切弹性模量泊松比G （2）平面应力问题0)()1 (2)(1)(10 xzyzyxzxyxyxyyyxxxzyzzEuEuuEuE可得出应力方程： 210001011.21.1)1 (2)(1)(12222uuuuEDDuuEuEuuEuuExyyxxyyxxyxyxyxyyyxxD平面应力问题的弹性矩阵，对称与E，有关 （3）平面应变问题的物理方程xyxyxyxyyyxxyxzzxyzzxzy
9、zuEuuEuuuEuuuEuu2221)11 (2)1 (2)1(1)1(1100）式中代入（平面应变的物理方程。 ）（换成换成2112uuuuEE )1 (22100011011)21)(1 ()1 (DuuuuuuuuuED平面应变弹性矩阵两类平面问题类似：平面应力问题 得到平面应变问题的方程式。 （4）平衡方程应力与外力间的关系：作用于物体的外力分为体积力和表面力，简称体力和面力。体力是分布在物体体积中的力，如重力和惯性力。面力是分布于物体表面上的力，如接触力。（1）应力与体力之间的平衡方程 XYCxyxxyyxyyyydyxyxyxdxxxyxydxxxxZ向原度=1 平衡微分方程0
10、000yxyxyxFyFxxyyyxx （1）应力和面力关系： DABYX1S2Sn dsxyBDAnYXxxyyxy 静力边界条件sincossincos00yxyyxxYXFyFx 三虚功方程 弹性力学在能量法求解中，特别是在用有限单元法求解中，要用到变形体的虚位移原理和虚功方程，这里简介一下。 所谓虚位移，是一种假想加到系统上，为系统的约束条件所允许的任意微小位移。而元功是指实际的力在虚位移上所作的功，元功也称为虚功。变形体虚功原理表述如下： 设变形体在力的作用下处于平衡状态，又设变形体由于别的原因产生符合约束条件的微小的连续变形（虚位移，虚应变），则外力在虚位移上所做的虚功W，恒等于应
11、力在虚应变上所做的虚变形功V。 即： W=V 其中：tdsvYuXtdxdyYvuXWSA1)()(*虚位移为板厚ttdxdyUAxyxyyyxx)(* V如果只有刚体的虚位移而没有虚应变，虚变形功V=0.故W=0，变成刚体的虚功方程PxiuD+Pyivi + TP=TdV(xx+yy+zz+xyzy+yzyz+ zxzx)dxdydz如果变形体处于平衡状态，当给予虚位移时，外力在虚位移上的虚功等于整个变形体内应力在虚应变上的虚功。在弹性力学与有限元法中，通常用虚位移方程代替平衡微分方程与应力边界条件。 注意：1.本节用“变形体”而没用“弹性体”，是因为 这里分析结构并不局限于弹性范围；2.虚
12、功方程的两种运用；（1）虚设位移求未知力。（虚位移原理）求位移。（虚功原理）给定力系.虚设位移.求未知力.（2）虚设力状态给定位移虚设力系.求未知位移. 4.2 4.2 弹性力学平面问题的有限单元法弹性力学平面问题的有限单元法一、弹性力学的有限元分析计算可分为三个步骤: 1、结构离散化 这是有限元法的基础，用由有限个方位不同但几何性质及物理性质均相似的单元组成的集合体来代替原来的连续体或结构。每个单元仅在节点处和其他单元及外部有联系。对于不同的问题，根据自身的特点，可选用不同类型的单元。对同一问题也可以分别或同时选用多种单元。 1 2 3 X2 Y2 1 2 12xF 12yF 11xF 11
13、yF 22yF 23yF 2 3 22xF 23xF 节点载荷节点载荷节点力节点力 不同材料节点不合法 2P2P2mLA1m 例：图示一悬臂梁，梁的厚度为t，设泊松比= ，弹性模量为E，试用三节点三角形单元进行离散。31 2p2pXYijij4132 2.单元分析主要内容：由节点位移求内部任意一点的位移，由节点位移求单元应变、应力和节点力。3.整体分析（1）由节点平衡方程，建立以整体刚度矩阵K为系数的，整体节点位移和外载R的关系式整体平衡方程。（2）考虑几何边界条件，修改总体刚度矩阵，求解全部未知位移分量。 P P 力学模型力学模型(平面应力问题平面应力问题)有限元模型有限元模型有限元法是一种
14、数值计算方法。可广泛应用于各种微分方程描述的场问题的求解。有限元法是一种数值计算方法。可广泛应用于各种微分方程描述的场问题的求解。 二二. . 有限元求解例（一维单元）有限元求解例（一维单元）引例:用有限元法求图1所示受拉阶梯杆的位移和应力。已知杆截面面积A（1）=210-4m2，A（2）=110-4m2，各段杆长L（1）=L（2）=0.1m；材料弹性模量E（1）=E（2）=2105MPa，作用于杆端的拉力F3=100N。x)1(A,)1(E)2(A)1(E)1(L)2(L3F0 a) 13212e)(xiFix)(ei)(eL)(ejijjFjxx b) c) 图1 受拉阶梯杆a)示意图 b
15、)有限元模型 c)单元图 1.单元划分 根据材料力学的平面假设，等截面受拉杆的同一截面的不同点上可认为具有相同的位移和应力，即位移只与截面的轴向坐标（图1中为x）有关，所以可将阶梯杆看作由两个“一维单元”组成，同一个单元内截面积及材料特性不变，并用线段表示一维单元。最简单的情况是，每一个单元有两个结点，他们分别位于单元两端。相邻两单元靠公共结点联结。这样，图1a所视的受拉阶梯杆就简化为由两个一维单元和三个结点构成的有限单元模型（图1b）。图中和是单元号，1、2、3是结点号。取结点位移作为基本未知量，应力由求得的结点位移算出。 )()(xe)( ei)(ej2.确定单元插值函数（形函数） 有限元
16、法将整个求解域离散为一系列仅靠公共结点联结的单元，而每一个单元本身却视为光滑连续体。单元内任一点的场变量（如本例中的位移）可由本单元的结点值根据场变量在单元中的假定分布规律（插值函数）插值求得。 本例中，每单元有两个结点，采用线性插值方式是适宜的。图1c是一典型单元图。两结点分别为i和j。设单元中坐标为x处的场变量为结点场变量值分别记为和。 )()()()()()()()()()()()()()()()()(eejeijiejeieiejejijieiijjiijeiejeiexxLxxLxxxxxxxxxxxxxxx)()(xxji根据线性插值关系得： （1） 是单元自由度列阵；形函数矩阵，
17、因为它与单元的结点坐标（即单元形状）和单元插值形状有关。 式中，L（e）=xj-xi是单元长度；)(eTejei)()( 称为单元 形函数矩阵的分量数目应与单元自由度数目相等。 对于两自由度线性插值单元由式（1）可知形函数矩阵的两分量为： )()()()(eijejiLxxxLxxx （2） )()()()()()()()()()()()()()(ejeiejeeeeiejeieeeFLEAFLEA)(eiF)(ejF3.单元方程（单元结点位移与结点力的关系）由等截面杆变形与拉力的关系（虎克定律）得到 式中，和分别为作用于单元e的结点i和（3）结点j的结点力。 )()()()()()()(11
18、11ejeiejeieeeFFLEA)()()(eeeFK)(eKTejeieFFF)()()(式（3）写成矩阵形式为：或简记为： 式中，称为单元特性矩阵，在力学问题称为单元结点力列阵。（4）（5）中常称为单元刚度矩阵，简称单元刚阵；式（5）称为单元方程。 )(eF)(eF)(eiF)(ejFj123 到目前为止，单元方程（4）或（5）尚不能求解，因为结点力列阵尚属未知。的分量和 外载荷之间的关系。 4.单元组集建立总体方程组 为获得总体方程 和总体自由度（、）的对应关系进行扩展。 是相邻单元作用于单元e的结点i和j的力，即属于单元之间的作用力。只有将具有公共结点的单元“组集”在一起才能确定上
19、述结点力和结点和组，必须先将单元方程按照局部自由度 （ ）i和 具体来说，单元1的扩展方程为：0000011011)2(2)1(1321)1()1()1(FFLEA 式中，各项上角码表示单元序号；下角码表示自由 度总体序号。（6） 单元2的扩展方程为：)2(3)2(2321)2()2()2(0110110000FFLEA （7）由于相邻两单元公共结点上的基本场变量（位移）相同，所以可将扩展后的各单元方程相加。 将式（6）和式（7）相加得：321321)2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2()2(0
20、0FFFLEALEALEALEALEALEALEALEA （8）上述组集过程可记为：NEeeNEeeFK1)(1)()( 式中，NE代表有限元模型的单元总数。（9） FK组集后的结果简记为： 式中，K称为总体特性矩阵（力学中常称为总体刚度矩阵和总刚阵），F称为总体结点载荷列阵。需指出的是，对单元的一个公共结点而言，除了有相邻单元作用于该结点的力之外，还可能有做用于该结点的外载荷（包括以后要讲到的当量结点载荷）。若一结点上无外载荷作用（如本例中结点2），则说明各相邻单元作用于该结点的力是平衡的，即该结点的结点合力为零。 0)(33)2()2()2(2)2()2()2(FLEALEA1022026
21、40441013216F 若某结点上有外载荷作用（如本例中结点3），则各单元作用于该结点的内力和（即方程（8）中第3式左端项的负值）与该结点的外载荷（F3）相平衡，即：这就是说，列阵F各分量的含义是作用于相应自由度（结点位移）上的结点外载荷。将相应数据代入式（8）得： （11）（10） 上式即为本题的总体线性代数方程组，但不能获得唯一解，因为上式中的矩阵是奇异的。这种奇异性不是因数据巧合造成的，而是有其必然性。原因在于总体方程组式（8）只考虑了力平衡条件，而只根据力平衡不能唯一地确定系统的位移，因为系统在有任意刚性位移的情况下仍可处于力平衡状态。为获得各结点位移的唯一解，必须消除可能产生的刚体
22、位移，即必须计入位移边界条件。 5.计入边界条件，解方程组 0123011022261032623本题的位移边界条件为那么，式（11）和说，可从式（11）中消去一个方程。譬如，代入后得： 解得：=0.2510-6 m,=0.7510-6m。，中只剩下两个待求的自由度去第一个方程并将。也就是舍（12）这与材料力学求得的结果相同。 6.计算单元应变和应力 )(e)(edxxdxee)()()()()()()()()(111)(eeeeeLBdxdx)()()()()()()(eeeeeBExEx由材料力学得知，单元中任一点的应变（x）（x） 将式（1）及式（2）代入上式得 式中，B称为单元应变一结
23、点位移转换矩阵。应力应变关系为: （15）与位移的关系为:（13）（14） MPaLE5 . 01 . 0)1025. 00(102)(65) 1 (21) 1 () 1 (MPaLE11 . 010)75. 025. 0(102)(55)2(32)2()2(对于单元1对于单元2 yxohvsjimiviujvjumvmu三有限元求解例（有限元求解例（线性三角形单元）（1）单元位移模式三角形三节点单元如上图。 126) 1 (),(),(654321yxyxvyxyxu ),(10000001),(),(),(654321yxmyxyxyxvyxuyxfi、j、m为节点。六个位移分量需六个待定
24、参数、设单元位移分量是坐标x.y的线性函数，即：写成矩阵的形式为： e )(321321321ayxuyxuyxummmjjjiiiAA11AA22AA33（2）由单元节点位移求位移参数设节点i,j,m坐标分别是xi,yi;xj,yj;xm,ym。把三个节点的坐标及其水平位移代入式（1）中得：解得：， 45616mjikkmmjjiimjikkmmjjiimjikkmmjjiimjikkmmjjiimjikkmmjjiimjikkmmjjiivcvcvcvcvbvbvbvbvavavavaucucucucububububuauauaua,6,5,4,3,2,121)(2121)(2121)(2
25、121)(2121)(2121)(21对v同理可列出、解出的结果如下：方程。 jmimiijmmjixxcyybyxyxa mji轮换)(2111121ijmijmiiimmjmmijiiyxyxyxyxyxyxyxyxyx式中 为三角形单元面积。 eA mjimjimjimjimjimjicccbbbaaacccbbbaaaA00000000000000000021 e将写成矩阵形式，有由单元节点位移求单元内部任一点位移f(x,y)=m(x,y)A ef(x,y)=m(x,y) ),(),(yxvyxummjjiimmjjiivuvuvuyxNyxNyxNyxNyxNyxN.),(00),(
26、),(00),(),(00),( emjiyxN),(1001f(x,y)= =Ni(x,y) Nj(x,y) Nm(x,y)=二阶单位阵 ),(yxNi)(21ycxbaiii mji轮换形函数物理意义：，i节点单位位移，其他 Ni，Nj，Nm是坐标的连续函数，它反映单元内位移的分布状态，称为位移的形状函数，简称形函数。矩阵N称为形函数矩阵。节点位移分量为0，单元内部产生位移分布形状Ni(x,y)= 形函数的性质： 1、在单元任一点上，三个形函数之和等 于1. 2、形函数Ni在i点的函数值为1，在j点及m点的函数值为零。3、三角形单元i,j,m在ij边上的形函数与第三个顶点的坐标无关。 例：
27、求图示单元和单元的形函数矩阵2p2pXYij4132mmij如无特殊说明：约定节点编号逆时针方向如无特殊说明：约定节点编号逆时针方向 yjxabi)0 , 0(m), 0(a)0 ,(byxab(b,a)mji(b,0) （a）单元、分别如上图所示(b) abyxyxyxyxyxyximmimiimijii21)(211.单元如图a所示。设a=1m,b=2m.(或直接由图形可知其面积) (2)求系数ai，aj，am，bi，bj，bm，ci，cj，cm bxxcbxxcxxcayybyybayybabyxyxayxyxayxyxaijmmijjmijimimjmjiijjimmiimjjmmji
28、0000 bxycxbayxNiiii)(21),(ayycxbayxNjjjj)(21),(aybxycxbayxNmmmm1)(21),( (3)求形函数矩阵 代入相关常数： aybxaybxbyaybxbxNNNNNNNmmjjii10010000000000aybxyxyyxxN10021002002ab21将a=1,b=2代入得：单元如图b所示， =。 ai=xjym-xmyj=ab aj=xmyi-xiym=ab am=xiyj-xjyi=-ab bi=yj-ym=-abj=ym-yi=0bm=yi-yj=aci=xm-xj=0 cj=xi-xm=-bcm=xj-xi=b（1）求常
29、数 bxycxbayxNiiii1)(21),(ayycxbayxNjjjj1)(21),(aybxycxbayxNmmmm1)(21),(2)求形函数矩阵 )1(00)1(10011001aybxaybxayaybxbxyxyxyyxx2100211001210021N=将a=1,b=2m代入上式得： N= mmjjiimmmmjjjjiiiixyyxvuvuvubccbbccbbccbxvyuyvxu.00000021四四 1.由节点位移求单元的应变单元应变矩阵 = xNiyNi iiiibccb0021简记为=Be B可写成分块的形式：B=Bi Bj BmBi= B称为应变矩阵，它的元素
30、都只与单元的几何性单元应变与单元结点位移关系(i, j, m )质有关的常量。这种单元称为平面问题的常应变三角形单元。 2.单元应力与结点位移关系=D e (求应变的表达式)e 应力矩阵=DB=Si Sj Sm=DB=S 1)平面应力问题： 代入D及B 得:S=Si Sj Sm. 对于平面应力： Si=iiiiiibccubucbsE2121)1 (22 (i, j, m )* *注：注：D D 的表达式见的表达式见P72 4-8P72 4-8，BB的表达式见的表达式见P86 4-34P86 4-34 21E1iiiiiibccbcbE)1 (221)1 (22111)21)(1 (2)1 (
31、2) 平面应变问题. 将上式中以代E，以代则子矩阵： (i, j, m )Si= 3.单元刚度矩阵yxyyx0mVmViViVjVjVij(a)mUiUjU *mV*mV*iV*iV*jV*jVij(b)m*mu*iu*ju*mv*iv*jv mmjjiiVUVUVUxyyxFe= = *mmjjiivuvuvu*xyyx节点虚位移列阵及虚应变：*e= *= AA(*e)TFe=*T tdxdy=B e知*=B*e*e)TFe=(*e)TBTtdxdy由(AA由于*e中的元素为常量，提至前，故： BTtdxdy （B,为常量，dxdy为面积）Fe= Ammmjmijmjjjiimijiikkk
32、kkkkkkFe=BT tFe=BT t=BTDt =BTDBet =KeKe=BTDBt或Ke=BTDBtdxdy单元刚度矩阵=e srsrsrsrsrsrsrsrbbcccbbucbccubccbbEt21212121)1 (42对于平面应力问题krs=BrTDBst = (r=i,j,m; s=i,j,m)21E1对于平面应变问题：将上式中的E换成换成可得k，式略。, mmmmjjjjiiiibccbbccbbccb0000002121baabaabb111001011000100001求例4.2（p84）单元的单元刚度矩阵解：（1）求矩阵BB=ab 和bi,bj,bm及 ci,cj,c
33、m得:B= 2100010112EbaababaaabbbE21211102110021000112（2）求矩阵S D=S=DB= baababaaabbbbaabaabbEabt212111021100210001110101010100100001122（3）求矩阵keke=BTDBt=BTSt= )211()21(12121)21()211(212111100212102121021210212101001)1 (222222222222222222baababaabbabababababaabbaabaabbaaababababbabbabbabEt= 代入a=1 b=2m 得： )8
34、11 ()412(141812)2141(221414110022141081041121uuuuuuuuuuuuuuuEtR对称可算出，当a=b时单元刚度矩阵与尺寸a,b无关。 单元刚体矩阵的特性：单元刚度矩阵的物理意义： 单元刚度矩阵ke表示了单元抵抗变形的能力，即表示了节点位移e与节点力Fe之间 krs表示点s发生单位位移时，在节点r上的关系。产生的节点力。 例：证明右图所示中单元刚度矩阵.k=kayx(0,a)(a,a)(2a,a)(0,0)(a,0)(2a,0) mmjjiimjimjibcbcbccccbbb00000021证明：由于单元刚度矩阵 ke=BTDBt 可知：当两个三角
35、形单元几何尺寸相同时，t值和单元面积值均相同；当两个单元的材料不难验证,.单元的上述br.cr (i,j,m )性质相同时，弹性矩阵D也时相同的。故ke是否相同，取决于矩阵B 值均相等。B= 结论：两个单元刚度矩阵ke相等的条件为：只要两单元的形状、大小、方向和单元弹性常数均相同，并且编号的方式也相同，如按逆时针方向编号为 i,j,m，直角顶点编号为m，则两个单元的刚度矩阵时相等的。刚度矩阵的一些重要性质：（1）对称性 单元刚度矩阵是一个对称矩阵（2）奇异性 单元刚度矩阵是一个奇异矩阵0ek，表明其逆矩阵不存在。也就是说，如果给定了单元节点位移可以得出惟一的节点力。反之，如果给出节点力却无法求
36、出确定的节点位移。这是因为单一单元来考虑所受约束时，可能存在不引起单元应力和节点力的刚性位移。（3）分块性质 单元刚度矩阵可以分块运算。 五.单元载荷的移置（离散时每个单元受载作用于节点上）1.原则：将单元载荷向节点处移置，按照虚功等效的原则进行。对于变形体，虚功等效是指原载荷与节点载荷在任何虚位移上做的虚功相等。当位移模式确定后，载荷移置是唯一的。静力等效是虚功等效的特例。 mmjjiiYXYXYX*2.载荷移置公式（1）集中力 设单元i,j,m中任一点M(x,y)处受有集中力P=PxPyT移置到该单元各节点处载荷列阵为Re= 假设该单元发生一微小虚位移，M点相应的虚位移为f*该单元各节点处
37、相应虚位移为 由静力等效原理，载荷与节点等效载荷载虚位移上所作虚功相等：(e)TRe=f*TpT Re=NTpRe=Xi Yi Xj Yj Xm YmT*由f*=N e代入：=Nipx Nipy Njpx Njpy Nmpx NmpyT yxppp psTdspNt（2）面力 设单元i,j,m的一边受有分布的面力将微分面积tds上的面力合力tds当作集中载荷，可得面力的移置公式:RRe e=t=t（3）体力 设单元i,j,m受有分布体力p=px pyT将微分体积tdxdy上的体力合力ptdxdy当作集中载荷dp同理可得： dxdypNtRTe sTdspNiiqjq),0(jy),(mmyxm
38、)0 ,0( jxy例：设三角形单元i,j,m的ij边作用有线性分布的法向载荷，i和j两点的压力集度分别为qi和qj，试用公式 Re =求其等效节点载荷。单元厚度为t， 节点坐标如下图所示。 解：（1）计算常数 ai=xjym-xmyj=0 bi=yj-ym=-ym ci=xm-xj=xm aj=xmyi-xiym=xmyj bj=ym-yi=ym cj=xi-xm=-xm am=xiyj-xjyi=0 bm=yi-yj=yi cm=xj-xi=0 212121212121(2)计算形函数 Ni=(ai+bix+ciy)= (- ymx +xmy)(aj+bjx+cjy)= (xmyi+(ym
39、-yi)x-xmy(am+bmx+cmy)= yixNj=Nm= tdsPNRTe tdsPNtdsPNtdsPNRmijmijLTLTLTe（3）计算等效节点载荷 由 得： 在边界jm和mi上的面力为零，故上式积分中后两项为0(a) 0_jjiiqqqyyyxPtdsqqqyyINININRjjiimjiLije0在ij边上的面力分量可表示为：代入（a）式中得： 积分沿逆时针方向（ij)，有ds=-dy 所以 000)(000_0_0_00dyxNdyxNdyxNtytdxNxNxNRmyjyiyymjieiiii Tjijiieqqqqt yR000)2(0)2 (6 TietqyR000
40、1012在ij边上x=0，代入 Ni Nj Nm中：若在ij边上有沿x方向均布的压力q，则有qi=qj=qm，代入得： 4.3 整体分析 步骤： 建立整体刚度矩阵引入支承条件解方程求位移任务：建立整个离散体系的平衡方程，即结构的总刚度方程。求单元应力 4.3.1 建 立 整 体 刚 度 矩 阵1234mmjmiji2p2p 2Y 2Y yp2X2 2 X2 2 xp2U2 2 U2 2 V2 2 （作用于单元的节点力）环绕节点移置节点载荷（节点力单元I施加）一单元加节点力V2 2 i j m (4 2 3)作用于节点的集中力作用于节点的集中力 0 xF2222XPXUexeee0yF2222Y
41、PYVeyeeee eRF22由节点的平衡条件，有平衡方程： 环绕节点2的所有单元求和。故： TeeVUF222 TYXR222mjimmmjmijmjjjiimijiimjikkkkkkkkkFFF单元e作用在2节点上的节点力绕节点2的各单元作用在2节点上的等效节点载荷之和及直接作用于2的集中力。每一个单元可用节点位移表示节点力，采用分块形式P89 4.46 2*1 2*2 2*1iiYX nmjininmimjijiiiikkkkF, 22RFe 2,2Rknemjinn 节点2对应单元的局部编号为i即I单元i节点，故单元i节点上的节点力：代入得 每个节点均可得到类似的方程，每个节点可写出
42、两个平衡方程，按节点的序号排列，用矩阵表示： Rk节点位移列阵结构的节点载荷列阵整体的刚度矩阵 TYXYXR2211 4321444342413433323124232221141312114321KKKKKKKKKKKKKKKKRRRR故 iiiYXR iiivujijijijiijKKKKK2,212,22, 1212, 12 i，j=14 i=14 整体刚度矩阵的集成规则：（1）先求出每个单元的刚度矩阵ke（2）将ke的每个子块kije进行 换号，换成对应的整体编号。（3）将换号后的子块填入整体刚度矩（4）若在同一位置上有几个单元的相应阵上对应的位置。子块填入同一位置，则进行叠加。 11
43、1412414442212422142kkkkkkkkkmji333234232224434244324kkkkkkkkkmji前述刚度矩阵 kkkkkkkkkkkkkkkkkkk4444434242413433322424232222211412110022 总体刚度矩阵的特性：（1）对称性 总体刚度矩阵是一个对称矩阵。因单元刚度矩阵升阶后对称性不变，由之合成的总体刚度矩阵自然是对称矩阵。 （2）奇异性 总体刚度矩阵行列式的值 0K （3）稀疏性 总体刚度矩阵是一个稀疏矩阵；即矩阵中的绝大多数元素为0，非0元素只占元素总数的很小的一部分。因为只有当节点ij相关时Kij才不是0。与一个节点相关
44、的节点数很少，故其非零元素少，绝大部分是0。 （4）带状分布规律 分布在以主对角线为中心的带状区域内。 讲例：接前述例题。解：（1）求各单元刚度矩阵（P90 4-50)求得单元的刚度矩阵为： 8941141812414232214141121002412102141081410418102412004112EtkI k=k （单元刚度矩阵相等） 111412414442212422kkkkkkkkk333234232224434244kkkkkkkkk(2)整体刚度矩阵k 换号的单元刚度矩阵： k=k= 将 的单元刚度矩阵填入 （P101 4-62） 444342413433323124232
45、2211413121188kkkkkkkkkkkkkkkkk如果都有就相加 33211226661443633161100000yxyxQQYQQvukkkkkkk123222444334330Yvukkkk22442432342330Yvkukvkuk其中第3行及第4行表达式为： 0000010000001000000000000000010000001233221144433433Yvuvuvukkkk等效：u1=v1=u3=v3=0 与支承条件及前述矩阵一致。划掉0对应的行及同行号的列，划掉相应元素可得，便于计算机运算。 多个单元：节点n的水平位移un=0 则k =R改为： k中的2n-
46、1行与列中主对角元素为1其他为0 载荷向量R中2n-1元素置0 若Vn=0 则2n行作修改。 4433221166656463565554534645444336353433100000000100000000000000000000000000001000000001vuvuvuvukkkkkkkkkkkkkkkk00202000pp仍以P68（课本P98 4.24）u1=v1=u4=v4=0= （4）解方程 求节点位移 k=R 将整体刚度矩阵K代入支承条件202089411414142322112890412104 23133222ppvuvuEt 31332213412272413073
47、23vuvuEt2020pp设得：故： 99. 888. 142. 85 . 13322EtPvuvu TEtP0099. 888. 142. 85 . 10031（5）求单元应力S 单元1： a=1 b=2 = 000042. 85 . 1613103161016110061312131002189EEtP580. 1281. 0844. 0tpxyyx =Se =（）= xxyxyy1.580.8440.281tptptp1.58tp同理可求单元23. 254. 15 . 4163p418. 0289. 0844. 0tP= （6）求节点力及支承反力Fe=kee单元的节点力 a=1,b=2
48、, =31k = 121331161121613112761316141161100616131031610121610611210614161004189Et 单元的节点力 F= PPPPPP070. 10 . 2281. 058. 1789. 0422. 000yxFFpRpRyx07. 1211pRpRyx071. 0244同理可得F 也进行验算由受力图可得 验算 总结总结1 1：有限元求解弹性力学平面问题步骤如下：有限元求解弹性力学平面问题步骤如下：1.整理原始数据，结构离散化，对单元和节点编号。2.求单元刚刚度矩阵ke3.用刚度集成法，形成刚度矩阵k4.求节点等效载荷，写出载荷列阵
49、k5.引入支承条件6.解k=R 求出节点位移 7.求单元应力e=Se8.求单元节点力 Fe=ke9.整理结果，作节点位移图及应力图。 作业P151 4.4 P152 4.10ee 4.4 线性三角形元解弹性力学平面问题的matlab程序：一、matlab 基础在 matlab 的提示符 符“”下输入命令. 3*4+5 ans= 17 Cos（30*pi/180） ans= 0.8660 x=4 x=4 2/sqrt(3+x) ans= 0.7559 如果不让 matlab 输出运算数据，在命令行的结尾输入分号： y=32; z=5; x=2*y-z; w=3*y+4*z W=116 mable
50、t区分大小写 x=1 x=1 x=2 x=2 使用hecp命令可以获得所有matlab命令的详细用法 hecp inv 下面的例子显示如何输入矩阵并实现简单的矩阵运算。 x=1 2 3;4 5 6;7 8 1 2 3x=5 6 8 9 47 y=2 ; 0 ; -3 y= 2 0 -3 w=x*y w= -7 -10 -13 求解下面的联立方程组采用高斯消去法解方程组:A= 2 -1 3 0 ;1 5 -2 4;2 0 3 -2; 1 2 3 4 b= 3; 1; -2 ; 2 b= 3 1 -2 2 x=Ab. x= 1.9259 -1.8148 -0.8889 1.5926A= 2 -1
