matlab矩阵位移法编程结构力学

1、矩阵位移法编程大作业(091210211)一、编制原理本程序的原理是基于结构力学矩阵位移法原理，以结构结点位移作基本未知量，将要分析的结构拆成已知节点力一结点力位移关系的单跨梁集合，通过强令结构发生待定的基本未知位移，在各个单跨梁受力分析结果的基础上通过保证结构平衡建立位移法的线性方程组，从而求得基本未知量。二、程序说明本程序是计算 3 层 11 跨框架右侧结点的位移和弯矩的程序，编译过程是按照矩阵位移法的先处理法进行的。首先将结构杆件的交汇点作为结点，共有 36 个结点和 108个位移编号，然后根据梁、柱、斜杆的不同分别建立单元刚度矩阵，然后转换为整体坐标系下的刚度矩阵，然后将所有杆件的单元

2、刚度矩阵整合成为总体刚度矩阵，在进行整合时连续运用 for 函数，最终形成 108 阶的总体刚度矩阵。然后通过对荷载的分析自己确定出荷载矩阵，直接写进程序。这样就可以把36 个结点的 108 个位移求得，然后再利用各个单元的单元刚度矩阵和所得的位移求得单元杆件的内力。离散化编号如下图：2U7LB043寓LSI羽曲LE2站L50L63站L64囱47国U65咄L53酸37一任L6738牌u&a38用LE940国L54L叩二二二DL3731L3632L3933L40L4135L423GL4325IL+42GL452JL29L3I29一hF5由L2S出L刘1L33L3fl小i；G卜L14liL

3、15而L1B21L1722Lie23L192iL20S12111L221：L2315U9117二J上口16L7LOL11L12LI5三、算法流程四、源代码%吉构力学大作业 3 层 11 跨框架矩阵位移法编程王贝 091210211h=input(输入单层高 h：);L=input(输入单跨度 L：);EIc=input(输入柱子的抗弯刚度 EIc：);EAc=input(输入柱子的抗压刚度 EAc);EIb=input(输入梁的抗弯刚度 EIb:);EAb=input(输入梁的抗压刚度 EAb);EIo=input(输入斜杆的抗弯刚度 EIo：);EAo=input(输入斜杆的抗压刚度 EAa

4、);q=input(输入侧向均布荷载集度 q:);T1=1,0,0,0,0,0;0,1,0,0,0,0;0,0,1,0,0,0;0,0,0,1,0,0;0,0,0,0,1,0;0,0,0,0,0,1;%角度为 0的转换矩阵T2=0,1,0,0,0,0;-1,0,0,0,0,0;0,0,1,0,0,0;0,0,0,0,1,0;0,0,0,-1,0,0;0,0,0,0,0,1;新度为 90的转换矩阵x=atan(h/L);T=cos(x),sin(x),0,0,0,0;-sin(x),cos(x),0,0,0,0;0,0,1,0,0,0;0,0,0,cos(x),sin(x),0;0,0,0,-s

5、in(x),cos(x),0;0,0,0,0,0,1;脸斗杆的转换矩阵T3=T;%梁的单元刚度矩阵kb0=EAb/L00-EAb/L00;012*EIb/(L*L*L)6*EIb/(L*L)0-12*EIb/(L*L*L)6*EIb/(L*L);06*EIb/(L*L)4*EIb/L0-6*EIb/(L*L)2*EIb/L;- EAb/L00EAb/L00;0-12*EIb/(L*L*L)-6*EIb/(L*L)012*EIb/(L*L*L)-6*EIb/(L*L);06*EIb/(L*L)2*EIb/L0-6*EIb/(L*L)4*EIb/L;%柱子的单元刚度矩阵kc0=EAc/h00-EA

6、c/h00;012*EIc/(h*h*h)6*EIc/(h*h)0-12*EIc/(h*h*h)6*EIc/(h*h);06*EIc/(h*h)4*EIc/h0-6*EIc/(h*h)2*EIc/h;输出结果结束0- EAc/h00EAc/h00;0-12*EIc/(h*h*h)-6*EIc/(h*h)012*EIc/(h*h*h)-6*EIc/(h*h);06*EIc/(h*h)2*EIc/h0-6*EIc/(h*h)4*EIc/h;%斜杆的单元刚度矩阵H=sqrt(h*h+L*L);ko0=EAo/H00-EAo/H00;012*EIo/(H*H*H)6*EIo/(H*H)0-12*EIo

7、/(H*H*H)6*EIo/(H*H);06*EIo/(H*H)4*EIo/H0-6*EIo/(H*H)2*EIo/H;- EAo/H00EAo/H00;0-12*EIo/(H*H*H)-6*EIo/(H*H)012*EIo/(H*H*H)-6*EIo/(H*H);06*EIo/(H*H)2*EIo/H0-6*EIo/(H*H)4*EIo/H;kb=T1*kb0*T1;%、体坐标下梁的单元刚度矩阵kc=T2*kc0*T2;%、体坐标下柱子的单元刚度矩阵ko=T3*ko0*T3;%、体坐标斜杆的单元刚度矩阵X=zeros(108,108);Y=zeros(108,108);Z=zeros(108

8、,108);温义 108 阶 0 矩阵K1=zeros(108,108);K2=zeros(108,108);K3=zeros(108,108);K4=zeros(108,108);K5=zeros(108,108);K6=zeros(108,108);K7=zeros(108,108);K8=zeros(108,108);K9=zeros(108,108);%m 梁杆单元矩阵整合到总体刚度矩阵的循环语句forii=1:11X(3 洵-2:3 洵+3,3*ii-2:3*ii+3)=kb;K1=K1+X;X=zeros(108,108);endforii=13:23Y(3*ii-2:3*ii+3

9、,3*ii-2:3*ii+3)=kb;K1=K1+Y;Y=zeros(108,108);endforii=25:35Z(3*ii-2:3*ii+3,3*ii-2:3*ii+3)=kb;K1=K1+Z;Z=zeros(108,108);end%m 柱杆单元矩阵整合到总体刚度矩阵的循环语句forjj=1:36K2(3*jj-2:3*jj,3*jj-2:3*jj 尸 kc(4:6,4:6);endforjj=1:24K3(3*jj-2:3*jj,3*jj-2:3*jj)=kc(1:3,1:3);endforjj=1:24K4(3*jj-2:3*jj,3*jj+34:3*jj+36)=kc(1:3,4

10、:6);endforjj=1:24K5(3*jj+34:3*jj+36,3*jj-2:3*jj 尸 kc(4:6,1:3);end%e 斜杆杆单元矩阵整合到总体刚度矩阵的循环语句forgg=3:12:27K6(3*gg-2:3*gg,3*gg-2:3*gg)=ko(4:6,4:6);endforgg=2:12:14K7(3*gg-2:3*gg,3*gg-2:3*gg)=ko(1:3,1:3);endforgg=2:12:14K8(3*gg-2:3*gg,3*gg+37:3*gg+39)=ko(1:3,4:6);endforgg=2:12:14K9(3*gg+37:3*gg+39,3*gg-2:

11、3*gg)=ko(4:6,1:3);endK=K1+K2+K3+K4+K5+K6+K7+K8+%9;体刚度矩阵P=zeros(108,1);P(1,1)=h*q;P(37,1)=h*q;P(73,1)=h*q/2;P(75,1)=q*h*h/12;A=KP;%吉构位移B1=kb*A(103:108,1);B2=kb*A(67:72,1);B3=kb*A(31:36,1);D1=zeros(6,1);D1(1:3,1)=A(70:72,1);D1(4:6,1)=A(106:108,1);D2=zeros(6,1);D2(1:3,1)=A(34:36,1);D2(4:6,1)=A(70:72,1)

12、;C1=kc*D1;C2=kc*D2;C3=kc(4:6,4:6)*A(34:36,1);M1(1,1)=B1(6,1);M2(1,1)=C1(6,1);M3(1,1)=B2(6,1);M4(1,1)=C2(6,1);M5(1,1)=C1(3,1);M6(1,1)=B3(6,1);M7(1,1)=C3(3,1);M8(1,1)=C2(3,1);fori=1:3m=36*i-2;fprintf(第dg 最右侧节点的位移是dn,i,A(m,1)endfprintf(第 1 层最右侧节点左侧杆的弯矩是%fn,M6)fprintf(第 1 层最右侧节点下侧杆的弯矩是%fn,M7)fprintf(第 1

13、 层最右侧节点上侧杆的弯矩是%fn,M8)fprintf(第 2 层最右侧节点左侧杆的弯矩是%fn,M3)fprintf(第 2 层最右侧节点下侧杆的弯矩是%fn,M4)fprintf(第 2 层最右侧节点上侧杆的弯矩是%fn,M5)fprintf(第 3 层最右侧节点左侧杆的弯矩是%fn,M1)fprintf(第 3 层最右侧节点下侧杆的弯矩是%fn,M2)五、试算算例输入数据：卒&入单层高 h：1输入单跨度 L:1输入柱子的抗弯刚度 EIc：1输入柱子的抗压刚度 EAc：1输入梁的抗弯刚度 EIb:1输入梁的抗压刚度 EAb:1输入斜杆的抗弯刚度 EIo：1输入斜杆的抗压刚度 EAo：1输入侧向均布荷载集度 q：1计算结果：第 1 层最右侧节点的位移是-6.219850e-003第 2 层最右侧节点的位移是-2.152659e-002第 3 层最右侧节点的位移是-4.131873e-002第 1 层最右侧节点左侧杆的弯矩是