结构力学大作业-matlab-矩阵位移(共11页).doc
精选优质文档-倾情为你奉上结构力学大作业一、编号二、输出结果的过程（编写完程序后）1.输入A=N×6的矩阵（单元数N不限） 2.通过函数k(A)得到结点位移列矩阵W4.通过函数li(W,A)得到各单元的局部坐标下的杆端力矩阵L。（该题的操作过程详见附录一，详细程序见附录二）三、程序设计过程根据上述的输出结果过程，可得，后期工作量较小，前期编程工作中，得到整体坐标下的单元刚度矩阵以及单元等效结点荷载都仅仅是码字工作，唯一有些难度的是定位向量的利用，即将整体坐标下的单元刚度矩阵以及单元等效结点荷载变换成该单元对应的整体刚度矩阵以及整体等效荷载。设计思路：定位向量的作用是定位，则可以试想根据定位向量得到行列变换矩阵，从而得到该单元对应的整体刚度矩阵。假设共有N个位移编号，则建立r=(N+6)×(N+6)的单位矩阵。kk=(N+6)×(N+6)的空矩阵，单元定位向量为rr(1,2,5,13,14,0)，整体坐标系下单元刚度矩阵为ke。其中kk左上角6×6置放该刚度矩阵ke；目标：通过矩阵行列变换，将kk的右下角N×N变成该单元对应的整体刚度矩阵，即得到相应的行变换矩阵(RR)以及列变换矩阵(RR)。（具体过程见附件 rr.m）四、计算结果（与求解器对比）整体刚度矩阵K结点等效荷载P得出：与求解器结果的比较（上述方法的竖直位移方向与求解器的方向相反 ）求解器建模过程：结点,1,0,0 结点,2,0.6,3 结点,3,1,5 结点,4,4,5 结点,5,4,3 结点,6,4,0单元,1,2,1,1,1,1,1,0 单元,2,3,1,1,0,1,1,1 单元,3,4,1,1,1,1,1,1单元,4,5,1,1,1,1,1,1 单元,5,6,1,1,1,1,1,1 单元,2,5,1,1,0,1,1,1单元,2,6,1,1,0,1,1,1单元荷载,3,3,1,0,1,90 单元荷载2,4,5,3,0.5,0,1,90 单元荷载,6,1,4,1/2,90结点支承,1,6,0,0,0,0 结点支承,6,6,0,0,0,0 单元材料性质,1,7,0.5,1,0,0,-1 杆端 1 杆端 2 - - 单元码 u -水平位移 v -竖直位移 ?-转角 u -水平位移 v -竖直位移 ?-转角- 1 0. 0. 0. -5. -17. 1. 2 -5. -17. -0. -6. -24. -0. 3 -6. -24. -0. -9. -24. 1. 4 -9. -24. 1. -6. -19. 2. 5 -6. -19. 2. 0. 0. 0. 6 -5. -17. -3. -6. -19. 2. 7 -5. -17. 5. 0. 0. 0.即，整体得到结构力学求解器 MATLAB计算易看出，除了精度稍微低了一些外，计算结果完全一致。五、内力MATLAB计算的单元杆端力下面是求解器中结果 杆端 1 杆端 2单元码 轴力 剪力 弯矩 轴力 剪力 弯矩 1 -2. -0. 0. -2. -0. 0. 2 -1. -0. 0. -1. -0. -0. 3 -0. 1. -0. -0. -1. 0. 4 -1. 0. 0. -1. -0. 0. 5 -3. -0. 0. -3. -1. -2. 6 -0. 2. -0. -0. -1. 0. 7 -0. -0. 0. -0. -0. -2.易得，除了精度以外，数据完全相同（部分杆件编号不同）。 弯矩图 剪力图轴力图六、总结虽然结果与求解器完全一致，但肯定的说编写的程序远远不及求解器的，在编写的程序中单元的定位向量需要手动输入，而求解器是自动完成的，方便了很多。但除此以外，在求解杆端力，结点位移上，应该是比较接近求解器的，同样对结构的单元数没有限制。编写程序中，因诸多因素以及时间安排不合理，导致没能花太多的时间在这上面，也使得可以臻于完美的功能没能实现。通过编写程序，也让我极大的理解了矩阵位移法，也让我意识到土木人也是应该学好程序的。附录一1.输入原始数据矩阵ac2.运行函数k(ac)，得到结点位移de3.运行函数li(de,ac)，得到杆端力L结束！易看出，工作量主要集中在输入原始数据的过程中，且容易输错。附录二rr.m （该函数得到的是根据定位向量，得到使局部刚度矩阵转换成整体刚度矩阵的矩阵。）function R=rr(a,b,c,d,e,f,g)a=a+6;b=b+6;c=c+6;d=d+6;e=e+6;f=f+6;g=g+6;R=eye(g);if a>6 r=eye(g); r(1,a)=1; r(a,1)=1; r(1,1)=0; r(a,a)=0; R=r*R;endif b>6 r=eye(g); r(2,b)=1; r(b,2)=1; r(2,2)=0; r(b,b)=0; R=r*R;endif c>6 r=eye(g); r(3,c)=1; r(c,3)=1; r(3,3)=0; r(c,c)=0; R=r*R;endif d>6 r=eye(g); r(4,d)=1; r(d,4)=1; r(4,4)=0; r(d,d)=0; R=r*R; endif e>6 r=eye(g); r(5,e)=1; r(e,5)=1; r(5,5)=0; r(e,e)=0; R=r*R; endif f>6 r=eye(g); r(6,f)=1; r(f,6)=1; r(6,6)=0; r(f,f)=0; R=r*R;endk.m（得到结点位移）function KKs=k(X)h=size(X,1);E=1;I=1;A=0.5;dw=zeros(h,6);dw=X(1:h,7:12);zs=max(max(dw);zzs=zs+6;KKK=zeros(zs);FF=zeros(zs,1);for di=1:h x1=X(di,1); y1=X(di,2); x2=X(di,3); y2=X(di,4);cc=(x2-x1)/sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);ss=(y1-y2)/sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);l=sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);p=X(di,5);q=X(di,6);RR=rr(X(di,7),X(di,8),X(di,9),X(di,10),X(di,11),X(di,12),zs);rrk=zeros(1,6);rrk=X(di,7:12);m1=-q*l*l/12;m2=q*l*l/12;f2y=-q*l/2;f1y=-q*l/2;f1x=0;f2x=0;m3=-p*l/8;m4=p*l/8;f3y=-p/2;f4y=-p/2;f1x=0;f2x=0;fe=zeros(1,6);fe=f1x,f1y,m1,f2x,f2y,m2+f1x,f3y,m3,f2x,f4y,m4;T=cc,ss,0,0,0,0;-ss,cc,0,0,0,0;0,0,1,0,0,0;0,0,0,cc,ss,0;0,0,0,-ss,cc,0;0,0,0,0,0,1;i=E*I/l;kee=E*A/l,0,0,-E*A/l,0,0;0,12*i/(l*l),6*i/l,0,-12*i/(l*l),6*i/l;0,6*i/l,4*i,0,-6*i/l,2*i;-E*A/l,0,0,E*A/l,0,0;0,-12*i/(l*l),-6*i/l,0,12*i/(l*l),-6*i/l;0,6*i/l,2*i,0,-6*i/l,4*i;ke=inv(T)*kee*T;kk=zeros(zzs);kk(1:6,1:6)=ke;kk=RR*kk*RR'kkk=zeros(zs);kkk=kk(7:zzs,7:zzs);KKK=KKK+kkk;fee=zeros(6,1);ff=zeros(zzs,1);fee=inv(T)*fe'ff(1:6,1)=fee;ff=RR*ff;f=zeros(zs,1);f=ff(7:zzs,1);FF=FF+f;endKKs=zeros(zs,1);KKs=-inv(KKK)*FF; li.m（得到杆端力）function L=li(de,X)h=size(X,1);E=1;I=1;A=0.5;dw=zeros(h,6);dw=X(1:h,7:12);zs=max(max(dw);zzs=zs+6;KKK=zeros(zs);FF=zeros(zs,1);L=zeros(6,h);for di=1:h x1=X(di,1); y1=X(di,2); x2=X(di,3); y2=X(di,4);cc=(x2-x1)/sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);ss=(y1-y2)/sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);l=sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);p=X(di,5);q=X(di,6);RR=rr(X(di,7),X(di,8),X(di,9),X(di,10),X(di,11),X(di,12),zs);rrk=zeros(1,6);rrk=X(di,7:12);m1=-q*l*l/12;m2=q*l*l/12;f2y=-q*l/2;f1y=-q*l/2;f1x=0;f2x=0;m3=-p*l/8;m4=p*l/8;f3y=-p/2;f4y=-p/2;f1x=0;f2x=0;fe=zeros(1,6);fe=f1x,f1y,m1,f2x,f2y,m2+f1x,f3y,m3,f2x,f4y,m4;T=cc,ss,0,0,0,0;-ss,cc,0,0,0,0;0,0,1,0,0,0;0,0,0,cc,ss,0;0,0,0,-ss,cc,0;0,0,0,0,0,1;i=E*I/l;kee=E*A/l,0,0,-E*A/l,0,0;0,12*i/(l*l),6*i/l,0,-12*i/(l*l),6*i/l;0,6*i/l,4*i,0,-6*i/l,2*i;-E*A/l,0,0,E*A/l,0,0;0,-12*i/(l*l),-6*i/l,0,12*i/(l*l),-6*i/l;0,6*i/l,2*i,0,-6*i/l,4*i;ke=inv(T)*kee*T;kk=zeros(zzs);kk(1:6,1:6)=ke;kk=RR*kk*RR'kkk=zeros(zs);kkk=kk(7:zzs,7:zzs);KKK=KKK+kkk;fee=zeros(6,1);ff=zeros(zzs,1);fee=inv(T)*fe'deb=zeros(6,1);for ii=1:6 if rrk(ii)=0; deb(ii,1)=0; else deb(ii,1)=de(rrk(ii),1); endendfb=zeros(6,1);fb=ke*deb-fee;fbz=T*fb;L(1:6,di)=-fbz;endKKs=zeros(zs,1);KKs=-inv(KKK)*FF; 专心-专注-专业
