1、第三章杆件系统的有限单元法有限元法与结构力学中刚架计算的矩阵位移法有密切关系，因此我们就从杆件系统的有限元法入手，来了解有限元法的基本概念和过程。杆单元(Truss or Bar)杆的定义：两端铰接，只受轴向力作用的基本构件 i j xui(FNi)uj(FNj)杆单元的结点位移和结点力()()NiijNjjiEAEAFuuFuullEAl根据材料力学的有关知识，我们可以立刻写出杆单元的结点位移与结点力之间的关系为为材料的弹性模量，为杆单元的长度 杆单元的刚度矩阵为杆单元的横截面面积1111NiieeeNjjFuEAFulfK 写成矩阵形式就是杆单元的刚度矩阵单元的刚度矩阵eTdKB DB平面
2、梁单元(Beam)平面梁单元的位移模式待定参数的确定MATLAB不仅可以进行数值运算，也能进行符号运算。如式不仅可以进行数值运算，也能进行符号运算。如式(3.20)中的矩中的矩阵阵Au和和Av的求逆运算，我们可以在的求逆运算，我们可以在MATLAB的命令窗口下输入的命令窗口下输入 syms L Au=1 0 1 L ;Av=1 0 0 0 0 1 0 0 1 L L2 L3 0 1 2*L 3*L2;第一句是定义符号变量第一句是定义符号变量L，后面定义两个矩阵，后面定义两个矩阵Au和和Av。然后我们再输入下。然后我们再输入下面求逆的命令面求逆的命令 inv(Au)ans=1,0 -1/L,1/
3、L inv(Av)ans=1,0,0,0 0,1,0,0 -3/L2,-2/L,3/L2,-1/L 2/L3,1/L2,-2/L3,1/L211011ullA12232321000010032312121vllllllllA位移模式应变矩阵B单元刚度矩阵 上面刚度矩阵的推导也可以应用上面刚度矩阵的推导也可以应用MATLAB的符号运算功能，计算的符号运算功能，计算的程序段如下的程序段如下syms E L x y A=1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 L 0 0 0 1 L 0 L2 L3 0 0 1 0 2*L 3*L2;Hu=1 0 0 x 0
4、 0;Hv=0 1 x 0 x2 x3;B=diff(Hu,x,1);-y*diff(Hv,x,2)*inv(A);Ke=int(E*transpose(B)*B,x,0,L)上述命令执行后，得到结果如下：上述命令执行后，得到结果如下：Ke=E/L,0,0,-E/L,0,0 0,12*E*y2/L3,6*E*y2/L2,0,-12*E*y2/L3,6*E*y2/L2 0,6*E*y2/L2,4*E*y2/L,0,-6*E*y2/L2,2*E*y2/L-E/L,0,0,E/L,0,0 0,-12*E*y2/L3,-6*E*y2/L2,0,12*E*y2/L3,-6*E*y2/L2 0,6*E*y
5、2/L2,2*E*y2/L,0,-6*E*y2/L2,4*E*y2/L03232322012064000126012062064leTAEdAdxEA lEI lEI lEI lEA lEA lEI lEI lEI lEI lEI lEI lEI lKB B对称2Iy dA等效结点力计算 2eTssBdBfN p所谓等效结点力，是指非结点载荷按照虚功相等的原则分配到单元结点上的力。等效表面力的普遍公式对于作用在杆件单元上的分布力，我们可以简写为 eTssldxfN p分布轴向力 推导线性分布轴力的等效结点力公式的推导线性分布轴力的等效结点力公式的MATLAB程序如下程序如下syms L x x
6、1 x2 p1 p2Au=1 0 1 L;H1=1 x;Nu=H1*inv(Au);x1=0;x2=L;L1=(x-x2)/(x1-x2);L2=(x-x1)/(x2-x1);p=p1*L1+p2*L2;fe=int(transpose(Nu)*p,x,0,L);fe=simple(fe);MATLAB给出的等效结点力公式是给出的等效结点力公式是fe=1/6*L*(p2+2*p1)1/6*L*(2*p2+p1)12(2)6NilFpp12(2)6NjlFpp0()leTsup xdxfN如果想推导其他分布形式的等效结点力，只要修改前面程序中的如果想推导其他分布形式的等效结点力，只要修改前面程序
7、中的p的形式。比的形式。比如我们要推导二次抛物线分布轴力的等效结点力，我们只要把如我们要推导二次抛物线分布轴力的等效结点力，我们只要把p写成二次的写成二次的Lagrange插值形式，即插值形式，即syms L x x1 x2 x3 p1 p2 p3Au=1 0 1 L;H1=1 x;Nu=H1*inv(Au);x1=0;x2=L/2;x3=L;L1=(x-x2)*(x-x3)/(x1-x2)/(x1-x3);L2=(x-x1)*(x-x3)/(x2-x1)/(x2-x3);L3=(x-x1)*(x-x2)/(x3-x1)/(x3-x2);p=p1*L1+p2*L2+p3*L3;fe=int(t
8、ranspose(Nu)*p,x,0,L);fe=simple(fe);执行后得到的等效结点力公式是执行后得到的等效结点力公式是fe=1/6*L*(p1+2*p2)1/6*L*(2*p2+p3)按这个思路可以按这个思路可以推导更高阶次分推导更高阶次分布轴力的等效结布轴力的等效结点力。点力。分布横向力 推导线性分布横向力等效结点力的推导线性分布横向力等效结点力的MATLAB程序为程序为syms L x x1 x2 p1 p2Av=1 0 0 0 0 1 0 0 1 L L2 L3 0 1 2*L 3*L2;H2=1 x x2 x3;Nv=H2*inv(Av);x1=0;x2=L;L1=(x-x2
9、)/(x1-x2);L2=(x-x1)/(x2-x1);p=p1*L1+p2*L2;fe=int(transpose(Nv)*p,x,0,L);fe=simple(fe);执行后，得到等效结点力为执行后，得到等效结点力为fe=1/20*L*(7*p1+3*p2)1/60*L2*(3*p1+2*p2)1/20*L*(3*p1+7*p2)-1/60*L2*(2*p1+3*p2)1221212212(73)20(32)60(37)20(23)60epp lpp lpp lpp lf0()leTsvp xdxfN分布弯曲力矩 syms L x x1 x2 m1 m2Av=1 0 0 0 0 1 0 0
10、 1 L L2 L3 0 1 2*L 3*L2;H2=1 x x2 x3;Nv=H2*inv(Av);x1=0;x2=L;L1=(x-x2)/(x1-x2);L2=(x-x1)/(x2-x1);m=m1*L1+m2*L2;fe=int(transpose(diff(Nv,x,1)*m,x,0,L);fe=simple(fe);0()leTszvm xdxfN单元的坐标变换 在前面几节中，推导单元刚度矩阵时采用的全部是局局部坐标系部坐标系，它的坐标轴方向是由单元即杆或梁的截面主方向确定的。采用这样的坐标系，可以得到具有统一形式的单元刚度矩阵。但是，实际的结构一般由具有不同方向和处于不同位置的杆或
11、梁所构成的。由不同方向单元组成的结构，它的整体刚度矩阵并不能由局部坐标下的单元刚度矩阵简单地叠加而成，所以必须建立一个统一的整体坐标系整体坐标系。计算时先将单元上的结点力和位移转换到整体坐标系，单元刚度矩阵亦作坐标变换，才可按叠加规则直接相加组成整体刚度矩阵。坐标变换平面杆单元的坐标转换矩阵 cos,siniiiiuuvucos0sin00cos0siniiijjjuvuuuvcos0sin00cos0sinTcosjixxlsinjiyyl22()()jijilxxyyiuivxyxyiu整体坐标系下的单元刚度矩阵22222222cos0sin011cossin000cos1100cossi
12、n0sincossineeTEAlccsccsccsscssEAslccsccscsscss KTK T平面梁单元的转换矩阵iuivxyxyiuiviiiiiiiivuvvuucossinsincosiiiiiivuvu1000cossin0sincos平面梁单元的转换矩阵cossin0000sincos0000001000000cossin0000sincos0000001T整体刚度矩阵13241234111121314151617181122232425262728123334353637382244454647482355565758336667683477784884xyxyxyxyf
13、kkkkkkkkufkkkkkkkvfkkkkkkufkkkkkvfkkkkufkkkvfkkufkv4?ijk 1324123413241234(1)(1)(1)111121(1)(2)(4)(2)(3)(1)(2)(4)22222111212222(2)(4)(4)(2)(4)311111233(3)(4)(3)(4)4222244kk00fkkkkkfffkkkffkkff集成后的整体刚度矩阵(e)ij()()()1112()()22eeeiieejjfkkfk()()()1112()()()1222eeeiieeejjijijkkfkkfeKef1eneeKKenee1fffK 整体刚
14、度矩阵整体结点力向量整体平衡方程整体刚度矩阵的特性边界约束条件的处理划行划列法111213141516112122232425262231323334353633414243444546445152535455565561626364656666kkkkkkfkkkkkkfkkkkkkfkkkkkkfkkkkkkfkkkkkkf 1144设有方程约束条件1122232526222112443233353633311344445253555655511544626365666661164410000000000001000000kkkkfkkkkkkfkkkkkkfkkkkkkfkk划行划列的程
15、序乘大数法11121314151611 11212223242526223132333435363341424344454642445152535455565561626364656666kkkkkkkkkkkkkfkkkkkkfkkkkkkkkkkkkkfkkkkkkf1111221331441551661 11kkkkkkk111111kk1111(2,3,4,5,6)jjkkj乘大数的程序第二类约束条件：例子1第二类约束条件：例子2第二类约束条件：例子3ij第二类约束条件的处理办法111314151216123133343532363414344454246425112533254425552225662561636465626661001000002(1)0kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkf12332442525222662kfkfkffkkfk第二类约束条件程序实现单元结点反力TeeexiyizixjyjzjFFMFFMfK 1、结点力中包含了分布力的等效结点力2、工程中经常要用的剪力、轴力是在单元的局部坐标系下定义的单元结点反力的程序温度应力温度变化结构变形结构受到约束产生温度应力00有限元中的温度应力处理梁单元的温度荷载温度荷载的程序算例一有限元程序框架算例二有限元离散模型
