xd=input(''); xn=size(xd,1); AA=1.0e30; % AA为大数 for i=1:xn K1(2*xd(i,1)-1,2*xd(i,1)-1)=K1(2*xd(i,1)-1,2*xd(i,1)-1)*AA; F(2*xd(i,1)-1)=xd(i,2)*K1(2*xd(i,1)-1,2*xd(i,1)-1); end fprintf('请输入y方向受到位移约束的结点编号及其约束位移值(形式为[node1 v1;node2 v2;……;noden vn]) \n '); yd=input(''); yn=size(yd,1); for i=1:yn K1(2*yd(i,1),2*yd(i,1))=K1(2*yd(i,1),2*yd(i,1))*AA; F(2*yd(i,1))=yd(i,2)*K1(2*yd(i,1),2*yd(i,1)); end u=K1\F; % 解整体刚度方程KU=P fprintf('边界修正后的整体刚度矩阵为K1=\n'); disp(K1); fprintf('边界修正后的外荷载列为F=\n'); disp(F); fprintf('在外荷载作用下各结点的位移为u=\n'); disp(u); 〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓 function k=k(co,E,NU,D,hou) % co为单元的结点坐标,为8行2列矩阵 % 函数返回单元的刚度矩阵 syms s t; % s,t表示局部坐标 N=sym(zeros(8,1)); N(1)=(1-s)*(1-t)*(-s-t-1)/4; N(2)=(1+s)*(1-t)*(s-t-1)/4; N(3)=(1+s)*(1+t)*(s+t-1)/4; N(4)=(1-s)*(1+t)*(-s+t-1)/4; N(5)=(1-t)*(1+s)*(1-s)/2; N(6)=(1+s)*(1+t)*(1-t)/2; N(7)=(1+t)*(1+s)*(1-s)/2; N(8)=(1-s)*(1+t)*(1-t)/2; x=0; y=0; for i=1:8 x=x+N(i)*co(i,1); y=y+N(i)*co(i,2); end (责任编辑:泉水) |