计算力学(有限元)matlab编程大作业(空间网架)
matlab网架编程
计算力学小论文《网架matlab编程》姓名：学号：专业：土木工程比例：本程序为适应正放四角锥、正放抽放四角锥、棋盘星四角锥三种类型的网架，可根据输入不同的参数进行变化。
前处理比较简单，运用方便。
本程序节点编号采用横向进行，先上层再下层。
这三种都可以运行，并已经附了真实的算例。
最后尝试了非线性的求解，但是结果不收敛，但在这里附上程序。
E=input('输入弹性模量（单位为Pa）:');A=input('输入各个单元截面面积（单位为平方米）:');m=input('输入横向网格个数:');n=input('输入纵向网格个数:');Lx=input('输入横向网格尺寸:');Ly=input('输入纵向网格尺寸:');h=input('输入网架厚度:');p=input('向下的节点力大小:');flag=input('网架的边界约束，0代表铰接；1代表刚接:');t1=input('第一个节点编号:');t2=input('第二个节点编号:');M=input('1:正放四角锥；2：正放抽空四角锥；3：棋盘星四角锥:'); if M==1%定位坐标for i=1:m+1for j=1:n+1x(i+(m+1)*(j-1))=(i-1)*Lx;y(i+(m+1)*(j-1))=(j-1)*Ly;z(i+(m+1)*(j-1))=h;endendfor i=1:mfor j=1:nx((m+1)*(n+1)+i+m*(j-1))=Lx/2+(i-1)*Lx;y((m+1)*(n+1)+i+m*(j-1))=Ly/2+(j-1)*Ly;z((m+1)*(n+1)+i+m*(j-1))=0;endend%建立总刚矩阵N=((m+1)*(n+1)+m*n)*3;u=zeros(N,1);while e<=10^-6for i=1:N/3x(i)=x(i)+u([i],1);y(i)=y(i)+u([i],1);z(i)=z(i)+u([i],1);endKK=zeros(N,N);for j=1:n+1for i=1:ma=i+(m+1)*(j-1);b=i+(m+1)*(j-1)+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:nfor i=1:ma=i+(m+1)*(j-1);b=i+(m+1)*j;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+(m+1)*(j-1);b=(m+1)*(n+1)+i+m*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+1+(m+1)*(j-1);b=(m+1)*(n+1)+i+m*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+(m+1)*j;b=(m+1)*(n+1)+i+m*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+1+(m+1)*j;b=(m+1)*(n+1)+i+m*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:nfor i=1:m-1a=(m+1)*(n+1)+i+m*(j-1);b=(m+1)*(n+1)+i+m*(j-1)+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:n-1for i=1:ma=(m+1)*(n+1)+i+m*(j-1);b=(m+1)*(n+1)+i+m*j;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endK=KK;F=K*u;P=zeros(((m+1)*(n+1)+m*n)*3,1);for i=1:(m+1)*(n+1)P(3*i-1,1)=-p;endP=P-F;Restrain=zeros(1,((m+1)*(n+1)+m*n)*3);if flag==1Restrain(1,(m+1)*(n+1)*3+1:((m+1)*(n+1)+m)*3)=1;Restrain(1,((m+1)*(n+1)+m*n-m)*3+1:((m+1)*(n+1)+m*n)*3)=1;for j=1:nRestrain(1,(1+(m+1)*(n+1)+m*(j-1))*3-2:(1+(m+1)*(n+1)+m*(j-1))*3)=1; Restrain(1,(m+(m+1)*(n+1)+m*(j-1))*3-2:(m+(m+1)*(n+1)+m*(j-1))*3)=1; endendif flag==0for i=1:mRestrain(1,((m+1)*(n+1)+i)*3-2:((m+1)*(n+1)+i)*3-1)=1;Restrain(1,((m+1)*(n+1)+m*n-m+i)*3-2:((m+1)*(n+1)+m*n-m+i)*3-1)=1; endfor j=1:nRestrain(1,(1+(m+1)*(n+1)+m*(j-1))*3-2:(1+(m+1)*(n+1)+m*(j-1))*3-1)=1; Restrain(1,(m+(m+1)*(n+1)+m*(j-1))*3-2:(m+(m+1)*(n+1)+m*(j-1))*3-1)=1; endendendif M==2%定位坐标for i=1:m+1for j=1:n+1x(i+(m+1)*(j-1))=(i-1)*Lx;y(i+(m+1)*(j-1))=(j-1)*Ly;z(i+(m+1)*(j-1))=h;endendm1=(m-1)/2+1;n1=(n-1)/2+1;m2=(m+1)*(n+1);for j=1:n1for i=1:mx(m2+i+(m+m1)*(j-1))=Lx/2+(i-1)*Lx;y(m2+i+(m+m1)*(j-1))=Ly/2+(j-1)*2*Ly;z(m2+i+(m+m1)*(j-1))=0;endfor j=1:n1-1for i=1:m1x(m2+m+i+(m+m1)*(j-1))=Lx/2+(i-1)*2*Lx;y(m2+i+(m+m1)*(j-1))=Ly/2+j*2*Ly;z(m2+i+(m+m1)*(j-1))=0;endend%建立总刚矩阵N1=((m+1)*(n+1)+m+(m+m1)*(n1-1))*3;KK=zeros(N1,N1);for j=1:n+1for i=1:ma=i+(m+1)*(j-1);b=i+(m+1)*(j-1)+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:nfor i=1:ma=i+(m+1)*(j-1);b=i+(m+1)*j;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:n1for i=1:m-1a=m2+i+(m+m1)*(j-1);b=m2+i+(m+m1)*(j-1)+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor i=1:m1for j=1:n1-1a=m2+2*i-1+(m+m1)*(j-1);b=m2+m+i+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m2+m+i+(m+m1)*(j-1);b=m2+2*i-1+(m+m1)*j;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:n1for i=1:ma=i+2*(m+1)*(j-1);b=m2+i+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+1+2*(m+1)*(j-1);b=m2+i+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+m+1+2*(m+1)*(j-1);b=m2+i+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+m+2+2*(m+1)*(j-1);b=m2+i+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:n1-1for i=1:m1a=2*i-1+m+1+2*(m+1)*(j-1);b=m2+2*i-1+m+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=2*i-1+m+1+2*(m+1)*(j-1)+1;b=m2+2*i-1+m+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=2*i-1+2*(m+1)+2*(m+1)*(j-1);b=m2+2*i-1+m+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=2*i+2*(m+1)+2*(m+1)*(j-1);b=m2+2*i-1+m+(m+m1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendK=KK;P=zeros(N1,1);for i=1:(m+1)*(n+1)P(3*i-1,1)=-p;endRestrain=zeros(1,N1);if flag==1Restrain(1,m2*3+1:(m2+m)*3)=1;Restrain(1,N1-m*3+1:N1)=1;for j=1:n1Restrain(1,(1+m2+(m+m1)*(j-1))*3-2:(1+m2+(m+m1)*(j-1))*3)=1;Restrain(1,(m+m2+(m+m1)*(j-1))*3-2:(m+m2+(m+m1)*(j-1))*3)=1;endfor j=1:n1-1Restrain(1,(1+m2+m+(m+m1)*(j-1))*3-2:(1+m2+m+(m+m1)*(j-1))*3)=1; Restrain(1,(m+m2+m1+(m+m1)*(j-1))*3-2:(m+m2+m1+(m+m1)*(j-1))*3)=1; endendif flag==0for i=1:mRestrain(1,(m2+i)*3-2:(m2+i)*3-1)=1;Restrain(1,(N1/3-m+i)*3-2:(N1/3-m+i)*3-1)=1;endfor j=1:n1Restrain(1,(1+m2+(m+m1)*(j-1))*3-2:(1+m2+(m+m1)*(j-1))*3-1)=1; Restrain(1,(m+m2+(m+m1)*(j-1))*3-2:(m+m2+(m+m1)*(j-1))*3-1)=1;endfor j=1:n1-1Restrain(1,(1+m2+m+(m+m1)*(j-1))*3-2:(m+m2+1+(m+m1)*(j-1))*3-1)=1; Restrain(1,(m+m2+m1+(m+m1)*(j-1))*3-2:(m+m2+m1+(m+m1)*(j-1))*3-1)=1; endendendif M==3%定位坐标for i=1:m+1for j=1:n+1x(i+(m+1)*(j-1))=(i-1)*Lx;y(i+(m+1)*(j-1))=(j-1)*Ly;z(i+(m+1)*(j-1))=h;endendm3=(m-1)/2+1;n3=(n-1)/2+1;m2=(m+1)*(n+1);N1=(m-3)/2*(2*m3+1)+m3+2*m+(m+1)*(m+1);for i=1:mx(m2+i)=Lx/2+(i-1)*Lx;y(m2+i)=Ly/2;z(m2+i)=0;x(N1-m+i)=Lx/2+(i-1)*Lx;y(N1-m+i)=Ly/2+(n-1)*Ly;z(N1-m+i)=0;endfor i=1:m3for j=1:n3-1x(m2+m+i+(2*m3+1)*(j-1))=Lx/2+(i-1)*2*Lx;y(m2+m+i+(2*m3+1)*(j-1))=Ly/2*3+(j-1)*2*Ly;z(m2+m+i+(2*m3+1)*(j-1))=0;endendfor i=1:m3-1for j=1:n3-2x(m2+m+m3+1+i+(2*m3+1)*(j-1))=Lx/2+(2*i-1)*Lx;y(m2+m+m3+1+i+(2*m3+1)*(j-1))=Ly/2*5+(j-1)*2*Ly;z(m2+m+m3+1+i+(2*m3+1)*(j-1))=0;endendfor j=1:n3-2x(m2+m+m3+1+(2*m3+1)*(j-1))=Lx/2;y(m2+m+m3+1+(2*m3+1)*(j-1))=Ly/2*5+(j-1)*2*Ly;z(m2+m+m3+1+(2*m3+1)*(j-1))=0;end%建立总刚矩阵N2=N1*3;KK=zeros(N2,N2);for j=1:n+1for i=1:ma=i+(m+1)*(j-1);b=i+(m+1)*(j-1)+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:nfor i=1:ma=i+(m+1)*(j-1);b=i+(m+1)*j;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor i=1:m-1a=m2+i;b=a+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=N1-m+i;b=a+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);enda=m2+1;b=a+m;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m2+m;b=a+m3;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=N1-m-m3+1;b=a+m3;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=N1-m;b=N1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);for j=1:n3-2a=m2+m+1+(2*m3+1)*(j-1);b=a+m3;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m2+m+1+m3+(2*m3+1)*(j-1);b=a+m3+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m2+1+(2*m3+1)+(2*m3+1)*(j-1);b=a+m3+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m2+1+(2*m3+1)+m3+1+(2*m3+1)*(j-1);b=a+m3;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);enda=m2;b=m2+m;for i=1:m3-1a=a+2;b=b+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);enda=N1-m;b=N1-m-m3+1;for i=1:m3-1a=a+2;b=b+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);enda=m2;b=m2+m+1;for i=1:m3-1a=a+2;b=b+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);enda=N1-m;b=a-m3;for i=1:m3-1a=a+2;b=b+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);enda=m2+m+1;b=a+m3;for j=1:m3for i=1:m3-1a=a+1;b=b+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endxx=mod(j,2);if xx==1a=a+2;b=b+2;elsea=a+1;b=b+1;endenda=m2+m+1;b=a+m3-1;for j=1:m3for i=1:m3-1a=a+1;b=b+1;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endxx=mod(j,2);if xx==1a=a+1;b=b+1;elsea=a+2;b=b+2;endendfor i=1:ma=i;b=m2+i;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+1;b=m2+i;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+m+1;b=m2+i;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=i+m+2;b=m2+i;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endfor i=1:ma=m2-m-1+i;b=N1-m-1+i;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m2-m+i;b=N1-m-1+i;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m2-2*m-2+i;b=N1-m-1+i;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m2-2*m-1+i;b=N1-m-1+i;k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endfor i=1:m3for j=1:n3-1a=m+1+(2*i-1)+2*(m+1)*(j-1);b=m2+m+i+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m+1+1+(2*i-1)+2*(m+1)*(j-1);b=m2+m+i+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m+1+m+1+(2*i-1)+2*(m+1)*(j-1);b=m2+m+i+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=m+1+m+2+(2*i-1)+2*(m+1)*(j-1);b=m2+m+i+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);endendfor j=1:n3-1a=2*(m+1)+1+2*(m+1)*(j-1);b=m2+m+m3+1+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=2*(m+1)+2+2*(m+1)*(j-1);b=m2+m+m3+1+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=3*(m+1)+1+2*(m+1)*(j-1);b=m2+m+m3+1+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=3*(m+1)+2+2*(m+1)*(j-1);b=m2+m+m3+1+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=2*(m+1)+1+2*(m+1)*(j-1);b=m2+m+m3+1+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=3*(m+1)-1+2*(m+1)*(j-1);b=m2+m+(2*m3+1)+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b)); KK=SpaceTrussAssemble(KK,k,a,b);a=3*(m+1)+2*(m+1)*(j-1);b=m2+m+(2*m3+1)+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));KK=SpaceTrussAssemble(KK,k,a,b);a=4*(m+1)-1+2*(m+1)*(j-1);b=m2+m+(2*m3+1)+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));KK=SpaceTrussAssemble(KK,k,a,b);a=4*(m+1)+2*(m+1)*(j-1);b=m2+m+(2*m3+1)+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));KK=SpaceTrussAssemble(KK,k,a,b);endfor i=1:m3-1for j=1:n3-2a=2*(m+1)+1+(2*i-1)+2*(m+1)*(j-1);b=m2+m+m3+1+i+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));KK=SpaceTrussAssemble(KK,k,a,b);a=2*(m+1)+2+(2*i-1)+2*(m+1)*(j-1);b=m2+m+m3+1+i+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));KK=SpaceTrussAssemble(KK,k,a,b);a=3*(m+1)+1+(2*i-1)+2*(m+1)*(j-1);b=m2+m+m3+1+i+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));KK=SpaceTrussAssemble(KK,k,a,b);a=3*(m+1)+2+(2*i-1)+2*(m+1)*(j-1);b=m2+m+m3+1+i+(2*m3+1)*(j-1);k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));KK=SpaceTrussAssemble(KK,k,a,b);endendK=KK;P=zeros(N1*3,1);for i=1:(m+1)*(n+1)P(3*i-1,1)=-p;endRestrain=zeros(1,N1*3);if flag==1Restrain(1,m2*3+1:(m2+m)*3)=1;Restrain(1,(N1-m)*3+1:N1*3)=1;for j=1:n3-1Restrain(1,(1+m2+m+(2*m3+1)*(j-1))*3-2:(1+m2+m+(2*m3+1)*(j-1))*3)=1; Restrain(1,(m2+m+m3+(2*m3+1)*(j-1))*3-2:(m2+m+m3+(2*m3+1)*(j-1))*3)=1; endfor j=1:n3-2Restrain(1,(m+m2+m3+1+(2*m3+1)*(j-1))*3-2:(m+m2+m3+1+(2*m3+1)*(j-1))*3)=1; Restrain(1,(m+m2+m3+1+m3+(2*m3+1)*(j-1))*3-2:(m+m2+m3+1+m3+(2*m3+1)*(j-1))*3)=1; endendif flag==0for i=1:mRestrain(1,(m2+i)*3-2:(m2+i)*3-1)=1;Restrain(1,(N1-m+i)*3-2:(N1/3-m+i)*3-1)=1;endfor j=1:n3-1Restrain(1,(1+m2+m+(2*m3+1)*(j-1))*3-2:(1+m2+m+(2*m3+1)*(j-1))*3-1)=1;Restrain(1,(m2+m+m3+(2*m3+1)*(j-1))*3-2:(m2+m+m3+(2*m3+1)*(j-1))*3-1)=1;endfor j=1:n3-2Restrain(1,(m+m2+m3+1+(2*m3+1)*(j-1))*3-2:(m+m2+m3+1+(2*m3+1)*(j-1))*3-1)=1; Restrain(1,(m+m2+m3+1+m3+(2*m3+1)*(j-1))*3-2:(m+m2+m3+1+m3+(2*m3+1)*(j-1))*3-1)=1 ;endendendwz=find(Restrain); %找出边界条件为1的for i=1:length(wz)KK(wz(i),:)=0; %总刚行置零KK(:,wz(i))=0; %总刚列置零KK(wz(i),wz(i))=1; %对角线元素置1P(wz(i),1)=0; %対应载荷置零endu=KK\P ;u' %求出总体坐标系下节点位移值Forces=K*u ;Forces' %求出总体坐标系下节点力%求节点编号为t1，t2的杆件的单元应力u1=[u(t1*3-2);u(t1*3-1);u(t1*3);u(t2*3-2);u(t2*3-1);u(t1*3)];y=SpaceTrussElementStress(E,x(t1),y(t1),z(t1),x(t2),y(t2),z(t2),u1)\%主程序中用到的函数：function k= SpaceTrussElementStiffness(E,A,x1,y1,z1,x2,y2,z2)L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));thetax=acos((x2-x1)/L);thetay=acos((y2-y1)/L);thetaz=acos((z2-z1)/L);Cx = cos(thetax);Cy = cos(thetay);Cz = cos(thetaz);w = [Cx*Cx Cx*Cy Cx*Cz ; Cy*Cx Cy*Cy Cy*Cz ; Cz*Cx Cz*Cy Cz*Cz]; k= E*A/L*[w -w ; -w w];function z=SpaceTrussAssemble(KK,k,i,j)%输入单元刚度矩阵k，单元的节点编号为i，jDOF(1)=3*i-2;DOF(2)=3*i-1;DOF(3)=3*i;DOF(4)=3*j-2;DOF(5)=3*j-1;DOF(6)=3*j;for n1=1:6for n2=1:6KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;function y = SpaceTrussElementStress(E,A,x1,y1,z1,x2,y2,z2,u)L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));thetax=acos((x2-x1)/L);thetay=acos((y2-y1)/L);thetaz=acos((z2-z1)/L);Cx = cos(thetax);Cy = cos(thetay);Cz = cos(thetaz);y = E/L*[-Cx -Cy -Cz Cx Cy Cz]*u;算例1输入弹性模量（单位为Pa）:2.2e11输入各个单元截面面积（单位为平方米）:0.0004输入横向网格个数:5输入纵向网格个数:5输入横向网格尺寸:0.5输入纵向网格尺寸:0.4输入网架厚度:0.4向下的节点力大小:50网架的边界约束，0代表铰接；1代表刚接:0第一个节点编号:1第二个节点编号:21:正放四角锥；2：正放抽空四角锥；3：棋盘星四角锥:1Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.423002e-19.u =1.0e+10 *Columns 1 through 180.0079 0.0587 0.1613 0.0079 0.0587 0.1514 0.0079 0.0587 0.1415 0.0079 0.0587 0.1316 0.0079 0.0587 0.1217 0.0079 0.8052 0.4851Columns 19 through 360.0079 0.0587 0.1026 0.0079 0.0587 0.0927 0.0079 0.0587 0.0828 0.0079 0.0587 0.0729 0.0079 0.0587 0.0630 0.0079 0.0587 0.0531Columns 37 through 540.0079 0.0587 0.0440 0.0079 0.0587 0.0341 0.0079 0.0587 0.0242 0.0079 0.0587 0.0143 0.0079 0.0587 0.0044 0.0079 0.0587 -0.0055Columns 55 through 720.0079 0.0587 -0.0147 0.0079 0.0587 -0.0246 0.0079 0.0587 -0.0345 0.0079 0.0587 -0.0444 0.0079 0.0587 -0.0543 0.0079 0.0587 -0.0642Columns 73 through 900.0079 0.0587 -0.0734 0.0079 0.0587 -0.0833 0.0079 0.0587 -0.0932 0.0079 0.0587 -0.1031 0.0079 0.0587 -0.1130 0.0079 0.0587 -0.1229Columns 91 through 1080.0079 0.0587 -0.1321 0.0079 0.0587 -0.1420 0.0079 0.0587 -0.1519 0.0079 0.0587 -0.1618 0.0079 0.0587 -0.1717 0.0079 -2.1771 0.93630 0 0.1270 0 0 0.1171 0 0 0.1072 0 0 0.0973 0 0 0.0874 0 0 0.0683Columns 127 through 1440.0000 -0.0000 0.0584 -0.0000 -0.0000 0.0485 -0.0000 -0.0000 0.0386 0 0 0.0287 0 0 0.0097 -0.0000 -0.0000 -0.0002Columns 145 through 162-0.0000 -0.0000 -0.0101 -0.0000 -0.0000 -0.0200 0 0 -0.0299 0 0 -0.0490 -0.0000 -0.0000 -0.0589 -0.0000 -0.0000 -0.0688Columns 163 through 180-0.0000 -0.0000 -0.0787 0 0 -0.0886 0 0 -0.1077 0 0 -0.1176 0 0 -0.1275 0 0 -0.1374Columns 181 through 1830 0 -0.1473Forces =Columns 1 through 180 -64.0000 16.0000 16.0000 -56.0000 48.0000 24.0000 -64.0000 64.0000 0 -80.0000 64.0000 0 -80.0000 80.0000 -48.0000 -56.0000 16.0000Columns 19 through 3624.0000 -92.0000 -16.0000 8.0000 -28.0000 -40.0000 32.0000 -32.0000 56.0000 -32.0000 4.0000 72.0000 -12.0000 -82.0000 72.0000 -12.0000 -54.0000 0Columns 37 through 54-11.0000 -41.0000 18.0000 21.4063 -114.2500 20.5313 20.0000 -26.0000 38.0000 -14.0000 -40.0000 16.0000 8.0000 -68.0000 4.0000 -4.0000 -50.0000 -12.0000Columns 55 through 724.0000 -40.0000 0 8.0000 -16.0000 40.0000 8.0000 -116.0000 48.0000 8.0000 -80.0000 32.0000 8.0000 -152.0000 64.0000 -16.0000 -56.0000 0Columns 73 through 900 -32.0000 16.0000 -16.0000 -64.0000 32.0000 -16.0000 -64.0000 0 16.0000 -64.0000 80.0000 48.0000 0 96.0000 16.0000 -48.0000 0Columns 91 through 108-8.0000 -64.0000 0 -16.0000 -16.0000 0 -16.0000 -16.0000 0 16.0000 -56.0000 -48.0000 -16.0000 -80.0000 -96.0000 48.0000 -144.0000 96.0000-104.2941 87.1209 -64.0000 -136.7383 206.8394 192.0000 7.8845 197.9146 -64.0000 56.5073 207.8739 -128.0000 48.0630 108.8262 64.0000 -34.9133 111.4466 -0.0000Columns 127 through 144-0.0891 -37.2942 64.0000 4.6814 -31.8969 32.0000 -3.7053 -9.7667 0.0000 17.6676 156.1519 0.0000 3.0587 162.9703 -24.0000 -11.9842 -5.6752 6.6250Columns 145 through 162-14.9397 -8.4694 32.0000 2.8699 -17.1463 16.0000 20.9631 163.8546 -16.0000 52.8506 112.6978 -64.0000 -15.7860 -18.1998 -128.0000 -27.1637 -12.8956 64.0000Columns 163 through 180-21.9125 -27.9474 -128.0000 52.6079 122.7961 -128.0000 112.8397 8.0503 0.0000 72.3954 153.9700 -192.0000 33.0182 208.8754 -128.0000 -118.3590 198.1338 -256.0000Columns 181 through 18317.1968 226.1485 128.0000y =-12288算例2输入弹性模量（单位为Pa）:2.2e11输入各个单元截面面积（单位为平方米）:0.0004输入横向网格个数:5输入纵向网格个数:5输入横向网格尺寸:0.5输入纵向网格尺寸:0.5输入网架厚度:0.4向下的节点力大小:50网架的边界约束，0代表铰接；1代表刚接:0第一个节点编号:1第二个节点编号:21:正放四角锥；2：正放抽空四角锥；3：棋盘星四角锥:2Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.422492e-51.ans =1.0e+24 *Columns 1 through 180.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 19 through 360.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 37 through 540.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 55 through 720.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 73 through 900.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000Columns 91 through 1080.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000Columns 109 through 1260 0 0.0000 0 0 -0.0000 0 0 -0.0000 0 0 -0.0000 0 0 -0.0000 0 0 0.0000Columns 127 through 1440.0000 0.0000 0.0000 0 0 -0.0000 0 0 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000Columns 145 through 1620 0 -0.0000 0 0 0.0000 5.3002 0.0000 -5.3002 0 0 -0.0000 0 0 0.0000 0 0 -0.0000Columns 163 through 1710 0 -0.0000 0 0 -0.0000 0 0 -0.0000ans =1.0e+03 *Columns 1 through 180.0027 -0.0480 -0.0025 -0.0100 -0.0520 0 -0.0520 -0.0400 0.0020 -0.0960 0.0320 0 -0.0320 -0.0960 0 -0.0320 -0.0320 0Columns 19 through 360.0780 -0.0480 -0.0160 -0.0040 -0.0520 0.0080 -0.0320 -0.0520 -0.0120 0 -0.1120 0.0800 0.0780 -0.1160 0.1839 0.0260 -0.0585 0.1007Columns 37 through 540.0190 -0.0480 -0.0030 0.0635 -0.0620 -0.0005 -0.0360 -0.0640 0.0040 -0.0960 -0.0960 0.0640 -0.0320 0.1600 0 -0.0320 -0.0640 -0.0640Columns 55 through 72-0.0437 -0.0440 0.0032 0.0516 -0.0480 -0.0027 -0.0040 -0.0600 0.0060 -0.0520 0.0240 -0.0440 0.0910 -0.1240 0.1005 0.0220 0.0083 -0.1404Columns 73 through 90-0.0080 -0.0480 0 0.0280 -0.0640 -0.0160 -0.0160 -0.0480 -0.0640 0 0 0 0 0.0960 -0.0640 -0.0320 -0.1280 0Columns 91 through 108-0.0080 -0.0480 0 -0.0040 -0.0560 -0.0080 0.0640 -0.0480 0.0640 -0.0640 0.0640 0.0640 -0.0960 0.0960 -0.0640 0.2240 0.2240 0.0640Columns 109 through 126-0.0404 0.0338 0.0000 -0.2985 -0.0160 0.0240 -0.2423 -0.0535 0.0160 -0.6021 -0.1815 -0.2560 -0.1834 0.0403 0.0000 -0.0216 0.2704 -0.0000Columns 127 through 144-0.0000 -0.0052 -0.0000 -0.0640 0.3071 0.0000 0.8166 -0.0950 -0.0080 -0.0033 0.0800 -0.0280 -0.0033 0.1860 -0.0320 -0.0005 0.0224 0.2560Columns 145 through 1620.9555 0.1959 -0.2560 -0.0206 -0.2296 -0.0005 0.0000 0.0009 0.0000 -0.0320 0.2504 0.0000 -0.0450 -0.1532 0.0960 -1.1203 0.7380 -0.0320Columns 163 through 1710.1491 0.0725 -0.2560 0.1532 0.0917 0.0000 -0.2543 0.0117 0.2560y =算例3输入弹性模量（单位为Pa）:2.2e11输入各个单元截面面积（单位为平方米）:0.0004输入横向网格个数:5输入纵向网格个数:5输入横向网格尺寸:0.5输入纵向网格尺寸:0.5输入网架厚度:0.4向下的节点力大小:50网架的边界约束，0代表铰接；1代表刚接:0第一个节点编号:1第二个节点编号:21:正放四角锥；2：正放抽空四角锥；3：棋盘星四角锥:3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.303872e-49.ans =1.0e+23 *Columns 1 through 180.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 19 through 360.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 37 through 540.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 55 through 720.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 73 through 900.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000Columns 91 through 1080.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000Columns 109 through 1260 0 0.0000 0 0 0.0000 0 0 -0.0000 0 0 -0.0000 0 0 -0.0000 0 0 0.0000Columns 127 through 1440.0000 -0.0000 -0.0000 0 0 -0.0000 0 0 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 0 0 -0.0000Columns 145 through 1620 0 0.0000 0.0000 -0.0000 0.0000 0 0 -0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000Columns 163 through 1680.0000 0.0000 -0.0000 0.0000 0.0000 -8.4171ans =1.0e+03 *Columns 1 through 180.0160 0 0 -0.0160 0 -0.0200 -0.0320 -0.0160 -0.0160 0.0640 0 0 -0.0320 -0.0960 0 -0.2240 0.1600 -0.1920Columns 19 through 360.0480 -0.0800 0 0.0480 -0.1440 -0.0960 0.0400 -0.1040 0 -0.0480 0.0720 0.0160 0.1600 -0.1600 -0.3840 -0.1280 -0.1280 0Columns 37 through 540.0640 0 -0.1920 -0.0800 -0.0800 -0.0320 0.0640 -0.0640 0 0.0960 -0.1280 0 0.0440 0.0640 -0.0320 0.0400 -0.0800 -0.0880Columns 55 through 720.0960 -0.1280 0.1920 -0.0320 0.0320 0.0640 -0.1160 -0.0840 0.0160 -0.0920 0.0360 0 -0.0320 -0.0640 -0.0320 -0.0640 -0.0320 -0.0320Columns 73 through 90-0.0640 -0.1280 0.1280 0.2560 -0.1280 -0.1280 -0.0640 -0.1280 0.0320 -0.0400 -0.0320 0.0320 -0.0320 -0.0800 -0.0160 0.1920 -0.1000 0.0200Columns 91 through 108-0.0640 -0.0320 0.1280 0.1280 -0.1280 0.1280 -0.0640 0.0640 0.0320 -0.0400 -0.0720 -0.0160 0.0960 -0.0800 -0.0160 0.0480 -0.0620 -0.0280Columns 109 through 1260.2252 -0.1035 0.1280 0.1578 -0.0242 -0.0480 0.1963 0.0053 0.0000 -0.1540 0.2248 0.2560 -0.4782 0.6520 0.5120 -0.0426 0.1847 0.6400Columns 127 through 1440.0454 -0.0262 0.0000 -0.2212 0.5215 0.0000 -0.2880 0.5694 -0.2560 -0.0857 0.0617 0.2560 0.0473 0.0089 0.1280 -0.0595 0.0159 -0.0240Columns 145 through 162-1.0624 0.4815 -0.2560 0.0266 0.0185 0.0640 0.0194 0.5300 0.2560 0.0646 -0.0320 -0.5120 -0.0590 -0.0132 -0.1280 0.0916 0.0182 0.0320Columns 163 through 168-0.0660 0.0950 0.0640 0 0.0020 0y =-131072E=input('输入弹性模量（单位为Pa）:');A=input('输入各个单元截面面积（单位为平方米）:');m=input('输入横向网格个数:');n=input('输入纵向网格个数:');Lx=input('输入横向网格尺寸:');Ly=input('输入纵向网格尺寸:');h=input('输入网架厚度:');p=input('向下的节点力大小:');flag=input('网架的边界约束，0代表铰接；1代表刚接:');t1=input('第一个节点编号:');t2=input('第二个节点编号:');M=input('1:正放四角锥；2：正放抽空四角锥；3：棋盘星四角锥:'); e=1;if M==1%定位坐标for i=1:m+1for j=1:n+1x(i+(m+1)*(j-1))=(i-1)*Lx;y(i+(m+1)*(j-1))=(j-1)*Ly;z(i+(m+1)*(j-1))=h;endendfor i=1:mfor j=1:nx((m+1)*(n+1)+i+m*(j-1))=Lx/2+(i-1)*Lx;y((m+1)*(n+1)+i+m*(j-1))=Ly/2+(j-1)*Ly;z((m+1)*(n+1)+i+m*(j-1))=0;endend%建立总刚矩阵N=((m+1)*(n+1)+m*n)*3;u=zeros(N,1);U=zeros(N,1);while e>=10^-6for i=1:N/3x(i)=x(i)+u([i],1);y(i)=y(i)+u([i],1);z(i)=z(i)+u([i],1);endKK=zeros(N,N);for j=1:n+1for i=1:ma=i+(m+1)*(j-1);b=i+(m+1)*(j-1)+1;。
(完整版)有限元大作业matlab---课程设计例子
有限元大作业程序设计学校：天津大学院系：建筑工程与力学学院专业：01级工程力学姓名：刘秀学号：\\\\\\\\\\\指导老师：连续体平面问题的有限元程序分析[题目]：如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力，m kNp 1=，同时在沿对角线y 轴上受一对集中压力，载荷为2KN ，若取板厚1=t ，泊松比0=v 。
[分析过程]：由于连续平板的对称性，只需要取其在第一象限的四分之一部分参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。
采用将此模型化分为4个全等的直角三角型单元。
利用其对称性，四分之一部分的边界约束，载荷可等效如图所示。
[程序原理及实现]：用FORTRAN程序的实现。
由节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。
模型基本信息由文件为BASIC.IN生成。
该程序的特点如下：问题类型：可用于计算弹性力学平面问题和平面应变问题单元类型：采用常应变三角形单元位移模式：用用线性位移模式载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷材料性质：弹性体由单一的均匀材料组成约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束方程求解：针对半带宽刚度方程的Gauss消元法输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN结果文件：输出一般的结果文件DATA.OUT程序的原理如框图：（1）主要变量：ID：问题类型码，ID＝1时为平面应力问题，ID=2时为平面应变问题N_NODE：节点个数N_LOAD：节点载荷个数N_DOF：自由度，N_DOF=N_NODE*2（平面问题）N_ELE：单元个数N_BAND：矩阵半带宽N_BC：有约束的节点个数PE：弹性模量PR：泊松比PT：厚度LJK_ELE(I,3)：单元节点编号数组，LJK_ELE(I,1)，LJK_ELE(I,2)，LJK_ELE(I,3)分别放单元I的三个节点的整体编号X(N_NODE), Y(N_NODE)：节点坐标数组，X(I),Y(I)分别存放节点I的x，y 坐标值P_LJK(N_BC,3)：节点载荷数组，P_LJK(I,1)表示第I个作用有节点载荷的节点的编号，P_LJK(I,2)，P_LJK(I,3)分别为该节点沿x,y方向的节点载荷数值AK(N_DOF,N_BAND)：整体刚度矩阵AKE(6,6)：单元刚度矩阵BB(3,6)：位移……应变转换矩阵（三节点单元的几何矩阵）DD(3,3)：弹性矩阵SS(3,6)；应力矩阵RESULT_N(N_NOF)：节点载荷数组，存放节点载荷向量，解方程后该矩阵存放节点位移DISP_E(6):：单元的节点位移向量STS_ELE(N_ELE,3)：单元的应力分量STS_ND(N_NODE,3)：节点的应力分量（2）子程序说明：READ_IN：读入数据BAND_K：形成半带宽的整体刚度矩阵FORM_KE：计算单元刚度矩阵FORM_P：计算节点载荷CAL_AREA：计算单元面积DO_BC：处理边界条件CLA_DD：计算单元弹性矩阵SOLVE：计算节点位移CLA_BB：计算单元位移……应变关系矩阵CAL_STS：计算单元和节点应力（3）文件管理：源程序文件：chengxu.for程序需读入的数据文件：BASIC.IN，NODE.IN，ELEMENT.IN（需要手工生成）程序输出的数据文件：DATA.OUT（4）数据文件格式：需读入的模型基本信息文件BASIC.IN的格式如下表需读入的节点信息文件NODE.IN的格式如下表需读入的单元信息文件ELEMENT.IN的格式如下表输出结果文件DATA.OUT格式如下表[算例原始数据和程序分析]：（1）模型基本信息文件BASIC.IN的数据为1,4,6,5,31.,0.,1.1,1,0,2,1,0,4,1,1,5,0,1,6,0,11,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5（2）手工准备的节点信息文件NODE.IN的数据为1 0.0 2.02 0.0 1.03 1.0 1.04 0. 0.5 1.0 0.6 2.0 0.（3）手工准备的单元信息文件ELEMENT.IN的数据为1 2 3 3 0 0 0 0 1 1 1 1 0 12 4 5 5 0 0 0 0 1 1 1 1 0 25 3 2 2 0 0 0 0 1 1 1 1 0 33 5 6 6 0 0 0 0 1 1 1 1 04 （4）源程序文件chengxu.for为：PROGRAM FEM2DDIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),&RESULT_N(500),AK(500,100)D IMENSION STS_ELE(500,3),STS_ND(500,3)OPEN(4,FILE='BASIC.IN')OPEN(5,FILE='NODE.IN')OPEN(6,FILE='ELEMENT.IN')OPEN(8,FILE='DATA.OUT')OPEN(9,FILE='FOR_POST.DAT')READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOADIF(ID.EQ.1)WRITE(8,20)IF(ID.EQ.2)WRITE(8,25)20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========')25 FORMAT(/5X,'=========PLANE STRAIN PROBLEM========')CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT, & IJK_ELE,X,Y,IJK_U,P_IJK)CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,& IJK_ELE,X,Y,PE,PR,PT,AK)CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, & STS_ELE,STS_ND)c to putout a data fileWRITE(9,70)REAL(N_NODE),REAL(N_ELE)70 FORMAT(2f9.4)WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),& STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)71 FORMA T(7F9.4)WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)),&REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)),&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE)72 FORMAT(7f9.4)cCLOSE(4)CLOSE(5)CLOSE(6)CLOSE(8)CLOSE(9)E NDcc to get the original data in order to model the problemSUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR, &PT,IJK_ELE,X,Y,IJK_U,P_IJK)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3), & P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14)REAL ND_ANSYS(N_NODE,3)READ(4,*)PE,PR,PTREAD(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC)READ(4,*)((P_IJK(I,J),J=1,3),I=1,N_LOAD)READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_NODE)READ(6,*)((NE_ANSYS(I,J),J=1,14),I=1,N_ELE)DO 10 I=1,N_NODEX(I)=ND_ANSYS(I,2)Y(I)=ND_ANSYS(I,3)10 CONTINUEDO 11 I=1,N_ELEDO 11 J=1,3IJK_ELE(I,J)=NE_ANSYS(I,J)11 CONTINUEN_BAND=0DO 20 IE=1,N_ELEDO 20 I=1,3DO 20 J=1,3IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J))IF(N_BAND.LT.IW)N_BAND=IW20 CONTINUEN_BAND=(N_BAND+1)*2IF(ID.EQ.1) THENELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)END IFR ETURNENDcC to form the stiffness matrix of elementSUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3), & AKE(6,6), SS(6,6)CALL CAL_DD(PE,PR,DD)CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 10 I=1,3DO 10 J=1,6SS(I,J)=0.0DO 10 K=1,310 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 20 I=1,6DO 20 J=1,6AKE(I,J)=0.0DO 20 K=1,320 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PTRETURNENDcc to form banded global stiffness matrixSUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE, & PR,PT,AK)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)N_DOF=2*N_NODEDO 40 I=1,N_DOFDO 40 J=1,N_BAND40 AK(I,J)=0DO 50 IE=1,N_ELECALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DO 50 I=1,3DO 50 II=1,2IH=2*(I-1)+IIIDH=2*(IJK_ELE(IE,I)-1)+IIDO 50 J=1,3DO 50 JJ=1,2IL=2*(J-1)+JJIZL=2*(IJK_ELE(IE,J)-1)+JJIDL=IZL-IDH+1IF(IDL.LE.0) THENELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)END IF50 CONTINUERETURNENDcc to calculate the area of elementSUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=(XIJ*YIK-XIK*YIJ)/2.0RETURNENDcc to calculate the elastic matrix of elementSUBROUTINE CAL_DD(PE,PR,DD)DIMENSION DD(3,3)DO 10 I=1,3DO 10 J=1,310 DD(I,J)=0.0DD(1,1)=PE/(1.0-PR*PR)DD(1,2)=PE*PR/(1.0-PR*PR)DD(2,1)=DD(1,2)DD(2,2)=DD(1,1)DD(3,3)=PE/((1.0+PR)*2.0)RETURNENDcc to calculate the strain-displacement matrix of elementSUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)DO 10 II=1,3DO 10 JJ=1,310 BB(II,JJ)=0.0BB(1,1)=Y(J)-Y(K)BB(1,3)=Y(K)-Y(I)BB(1,5)=Y(I)-Y(J)BB(2,2)=X(K)-X(J)BB(2,4)=X(I)-X(K)BB(2,6)=X(J)-X(I)BB(3,1)=BB(2,2)BB(3,2)=BB(1,1)BB(3,3)=BB(2,4)BB(3,4)=BB(1,3)BB(3,5)=BB(2,6)BB(3,6)=BB(1,5)CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DO 20 I1=1,3DO 20 J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURNENDcc to form the global load matrixSUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3), & RESULT_N(N_DOF)DO 10 I=1,N_DOF10 RESULT_N(I)=0.0DO 20 I=1,N_LOADII=P_IJK(I,1)RESULT_N(2*II-1)=P_IJK(I,2)20 RESULT_N(2*II)=P_IJK(I,3)RETURNENDcc to deal with BC(u) (here only for fixed displacement) using "1-0" method SUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)DO 30 I=1,N_BCIR=IJK_U(I,1)DO 30 J=2,3IF(IJK_U(I,J).EQ.0)THENELSEII=2*IR+J-3AK(II,1)=1.0RESULT_N(II)=0.0DO 10 JJ=2,N_BAND10 AK(II,JJ)=0.0DO 20 JJ=2,II20 AK(II-JJ+1,JJ)=0.0END IF30 CONTINUERETURNENDcc to solve the banded FEM equation by GAUSS eliminationSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),AK(500,100)DO 20 K=1,N_DOF-1IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1IF(N_DOF.LE.K+N_BAND-1)IM=N_DOFDO 20 I=K+1,IML=I-K+1C=AK(K,L)/AK(K,1)IW=N_BAND-L+1DO 10 J=1,IWM=J+I-K10 AK(I,J)=AK(I,J)-C*AK(K,M)20 RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)DO 40 I1=1,N_DOF-1I=N_DOF-I1IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1IF(N_BAND.LE.N_DOF-I-1)JQ=N_BANDDO 30 J=2,JQK=J+I-130 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40 RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)50 FORMAT(/12X,'* * * * * RESULTS BY FEM2D * * * * *',//8X,&'--DISPLACEMENT OF NODE--'//5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP') DO 60 I=1,N_NODE60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)70 FORMAT(8X,I5,7X,2E15.6)RETURNENDcc calculate the stress components of element and nodeSUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, &STS_ELE,STS_ND)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6), &SS(3,6),RESULT_N(N_DOF),DISP_E(6)DIMENSION STS_ELE(500,3),STS_ND(500,3)WRITE(8,10)10 FORMAT(//8X,'--STRESSES OF ELEMENT--')CALL CAL_DD(PE,PR,DD)DO 50 IE=1,N_ELECALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 20 I=1,3DO 20 J=1,6SS(I,J)=0.0DO 20 K=1,320 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 30 I=1,3DO 30 J=1,2IH=2*(I-1)+JIW=2*(IJK_ELE(IE,I)-1)+J30 DISP_E(IH)=RESULT_N(IW)STX=0STY=0TXY=0DO 40 J=1,6STX=STX+SS(1,J)*DISP_E(J)STY=STY+SS(2,J)*DISP_E(J)40 TXY=TXY+SS(3,J)*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STYSTS_ELE(IE,3)=TXY50 WRITE(8,60)IE,STX,STY,TXY60 FORMAT(1X,'ELEMENT NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)c the following part is to calculate stress components of nodeWRITE(8,55)55 FORMAT(//8X,'--STRESSES OF NODE--')DO 90 I=1,N_NODEA=0.B=0.C=0.II=0DO 70 K=1,N_ELEDO 70 J=1,3IF(IJK_ELE(K,J).EQ.I) THENII=II+1A=A+STS_ELE(K,1)B=B+STS_ELE(K,2)C=C+STS_ELE(K,3)END IF70 CONTINUESTS_ND(I,1)=A/IISTS_ND(I,2)=B/IISTS_ND(I,3)=C/IIWRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)75 FORMAT(1X,'NODE NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)90 CONTINUERETURNENDc FEM2D programm end[算例结果]：chengxu.for所输出的数据文件DATA.OUT数据内容如下：=========PLANE STRESS PROBLEM========* * * * * RESULTS BY FEM2D * * * * *--DISPLACEMENT OF NODE--NODE NO X-DISP Y-DISP1 .000000E+00 -.525275E+012 .000000E+00 -.225275E+013 -.108791E+01 -.137363E+014 .000000E+00 .000000E+005 -.824176E+00 .000000E+006 -.182418E+01 .000000E+00--STRESSES OF ELEMENT--ELEMENT NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00ELEMENT NO.= 2STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00ELEMENT NO.= 3STX=-.108791E+01 STY=-.137363E+01 TXY= .307692E+00ELEMENT NO.= 4STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00--STRESSES OF NODE--NODE NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00NODE NO.= 2STX=-.100000E+01 STY=-.220879E+01 TXY= .249084E+00NODE NO.= 3STX=-.105861E+01 STY=-.191575E+01 TXY= .205128E+00NODE NO.= 4STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00NODE NO.= 5STX=-.970696E+00 STY=-.166667E+01 TXY= .586081E-01NODE NO.= 6STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00[结论与体会]:通过本次的课程设计，我对有限元的概念有了更加深刻的理解，同时也弥补了平时学习是疏忽的地方，充实了有限元知识。
matlab有限元三单元编程
matlab有限元三单元编程MATLAB是一种功能强大的编程语言和环境，广泛用于工程和科学领域的数值计算、数据分析和可视化。
在有限元分析中，MATLAB的强大功能和直观的语法使其成为一个理想的选择。
在本文中，我们将讨论MATLAB中有限元三单元的编程方法和实践。
有限元分析是一种数值方法，用于解决连续介质的力学问题。
它将一个复杂的结构分解成更简单的有限元单元，然后通过求解线性代数方程组来得到结构的应力和位移解。
在有限元分析中，三角形和四边形单元是最常用的有限元单元之一。
本文将重点讨论三角形单元的编程实现。
首先，我们需要定义一个三角形单元的几何信息。
在三角形单元中，我们有三个顶点，每个顶点有两个坐标。
我们可以使用一个3x2的矩阵来表示这些坐标。
例如，定义三角形ABC的顶点坐标矩阵为P：P = [x_A, y_A;x_B, y_B;x_C, y_C];接下来，我们需要定义三角形单元的连接性信息。
在MATLAB中，我们可以使用一个3x1的向量来表示三个顶点的连接性。
例如，定义三角形ABC的连接性向量为E：E = [1;2;3];然后，我们可以定义三角形单元材料属性和载荷信息。
这些信息包括杨氏模量、泊松比和外力向量。
我们可以将这些信息存储在一个结构体中，例如：properties.E = 210e9; % 杨氏模量properties.v = 0.3; % 泊松比properties.f = [0; -1000; 0]; % 外力向量接下来，我们可以使用这些信息来计算三角形单元的刚度矩阵和力向量。
刚度矩阵是一个3x3的矩阵，力向量是一个3x1的向量。
我们可以使用以下公式来计算它们：K = (E/(1-v^2)) * [1 v 0;v 1 0;0 0 (1-v)/2] * A;f = (A/3) * [1; 1; 1] * properties.f;其中，A是三角形的面积。
我们可以使用以下公式来计算它：A = abs(det([1, 1, 1;x_A, x_B, x_C;y_A, y_B, y_C])) / 2;最后，我们可以使用这些信息来求解三角形单元的位移解。
matlab课程设计大作业
matlab课程设计大作业一、教学目标本课程的教学目标是使学生掌握MATLAB基本语法、编程技巧以及MATLAB 在工程计算和数据分析中的应用。
通过本课程的学习，学生将能够熟练使用MATLAB进行简单数学计算、线性方程组求解、函数图像绘制等。
1.掌握MATLAB基本语法和编程结构。
2.了解MATLAB在工程计算和数据分析中的应用。
3.熟悉MATLAB的函数库和工具箱。
4.能够使用MATLAB进行简单数学计算。
5.能够使用MATLAB求解线性方程组。
6.能够使用MATLAB绘制函数图像。
7.能够利用MATLAB进行数据分析和处理。
情感态度价值观目标：1.培养学生对计算机辅助设计的兴趣和认识。
2.培养学生团队合作和自主学习的能力。
二、教学内容本课程的教学内容主要包括MATLAB基本语法、编程技巧以及MATLAB在工程计算和数据分析中的应用。
1.MATLAB基本语法：介绍MATLAB的工作环境、基本数据类型、运算符、编程结构等。
2.MATLAB编程技巧：讲解MATLAB的函数调用、脚本编写、函数文件编写等编程技巧。
3.MATLAB在工程计算中的应用：介绍MATLAB在数值计算、线性方程组求解、图像处理等方面的应用。
4.MATLAB在数据分析中的应用：讲解MATLAB在数据采集、数据分析、数据可视化等方面的应用。
三、教学方法本课程采用讲授法、案例分析法、实验法等多种教学方法相结合的方式进行教学。
1.讲授法：通过讲解MATLAB的基本语法、编程技巧以及应用案例，使学生掌握MATLAB的基本知识和技能。
2.案例分析法：通过分析实际工程案例，使学生了解MATLAB在工程计算和数据分析中的应用。
3.实验法：安排上机实验，使学生在实际操作中巩固所学知识，提高实际编程能力。
四、教学资源本课程的教学资源包括教材、实验设备、多媒体资料等。
1.教材：选用《MATLAB教程》作为主要教材，辅助以相关参考书籍。
2.实验设备：为学生提供计算机实验室，配备有MATLAB软件的计算机。
有限元-计算结构力学-大作业
有限元-计算结构力学-大作业本页仅作为文档页封面，使用时可以删除This document is for reference only-rar21year.MarchSHANGHAI JIAO TONG UNIVERSITY 平面应力问题解的Matlab实现姓名: heiya168 学号: 帆哥班级：指导老师：目录1绪论 (4)2平面问题的四节点四边形单元 (4)2.1单元的构造 (4)2.2等参变换 (5)2.3边界条件的处理——置“1”法 (7)3有限元分析流程 (8)3.1程序原理和流程 (8)3.2使用的函数 (9)3.3文件管理 (9)3.4数据文件格式 (10)4算例——开方孔的矩形板拉伸分析 (11)4.1问题的具体参数与载荷 (11)4.2Matlab程序计算 (11)4.3ANSYS建模计算 (13)4.4误差分析 (15)5总结 (15)参考文献 (16)附录 (17)1绪论有限元方法(finite element method)，是求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要基础性原理。
将它用于在科学研究中，可成为探究物质客观规律的先进手段。
将它应用于工程技术中，可成为工程设计和分析的可靠工具。
弹性体在载荷作用下，其基本方程可写成以下的三类方程和两种边界条件。
平衡方程——应力与外载荷的关系；几何方程——应变位移关系；物理方程——应力应变关系；力的边界条件；几何边界条件。
应用最小位能原理，并利用上述关系，最终建立由刚度方程，节点位移和等效节点载荷所构成的求解方程。
带入边界条件求解方程，就可以得出弹性力学问题的一般性解答。
本次大作业基于有限元方法的基本原理，使用Matlab这一平台，针对平面应力问题，采用四节点四边形单元编写了求解单元节点位移的程序。
主要内容包括：1）介绍有限元的基本原理；2）编程基本思路及流程介绍；3）程序原理及说明； 4）具体算例这四个部分。
matlab编译平面有限元计算
matlab编译平面有限元计算编译平面有限元计算是一种常用的数值计算方法，可以用于求解各种复杂的工程问题。
在本文中，我们将介绍如何使用MATLAB编写平面有限元计算程序，并通过一个实例来说明其应用。
让我们了解一下有限元方法的基本原理。
有限元方法是一种将连续体划分为有限个单元，根据物理方程和边界条件，在每个单元上建立离散方程组，最终求解得到整个连续体的近似解的方法。
在平面问题中，连续体被划分为三角形或四边形单元。
MATLAB是一种功能强大的数值计算软件，它提供了丰富的工具箱和函数，可以方便地进行有限元计算。
下面我们将以求解平面弹性力学问题为例，介绍如何使用MATLAB编写平面有限元计算程序。
我们需要定义问题的几何信息、边界条件和材料参数。
在MATLAB中，可以通过定义节点坐标、单元连接关系、边界条件和材料参数来描述问题。
节点坐标可以用一个矩阵表示，其中每一行代表一个节点的坐标。
单元连接关系可以用一个矩阵表示，其中每一行代表一个单元的节点编号。
边界条件可以用一个向量表示，其中每个元素代表一个节点的边界条件。
材料参数可以用一个矩阵表示，其中每一行代表一个单元的材料参数。
接下来，我们需要建立有限元离散方程组。
在平面弹性力学问题中，离散方程组可以通过组装单元刚度矩阵和外力向量得到。
单元刚度矩阵可以通过单元的几何信息和材料参数计算得到。
外力向量可以根据边界条件和载荷信息计算得到。
最终，我们可以通过求解离散方程组得到节点的位移解。
我们可以根据节点的位移解计算出单元的应力和应变。
在平面弹性力学问题中，应力可以通过单元的几何信息和位移解计算得到。
应变可以通过应力和材料参数计算得到。
通过计算应力和应变，我们可以评估结构的稳定性和性能。
MATLAB编译平面有限元计算程序是一种强大的工具，可以用于求解各种复杂的工程问题。
通过定义几何信息、边界条件和材料参数，建立有限元离散方程组，求解得到节点的位移解，最终计算出单元的应力和应变。
结构分析的有限元法与matlab程序设计
结构分析的有限元法与matlab程序设计有限元法是一种结构分析的数值分析方法，它基于有限元变分原理，通过有限个单元来拟合近似满足实际结构的力学模型，以解决扩展性非常大的形状复杂的结构问题。
在实际应用中，MATLAB程序设计可以极大地提高有限元法的实用性，因为它使算法的实现更加容易。
程序设计过程主要分为以下几个步骤：
首先，根据实际应用情况，建立结构的物理模型，第二，确定要求的结构的抗力参数，也就是确定边界条件；然后，通过建立模型，根据力学原理将原有结构分解成若干有限元，第四步，建立有限元内力函数和节点变量；再次，解决建立的有限元模型系数矩阵，解出系统未知形变量；最后，利用有限元的形变量计算正确的结构问题。
要完成上述算法步骤，MATLAB可以提供更加有效的工具。
MATLAB 有一个可以直接进行有限元分析的工具箱——FEA TOOLBOX，可以直接用它来计算形变量和结构参数。
另外，MATLAB中还有大量可以用来处理有限元问题的函数和工具，如系数矩阵求解函数等。
在开发有限元法的程序设计中，MATLAB是一个很方便的工具。
我们可以充分利用MATLAB中的基本功能，对有限元法进行灵活的操作。
同时，通过程序的可视化，可以更好地了解模型的状态，并降低有限元算法分析的误差和风险。
总之，有限元法与MATLAB程序设计是结构分析中同一不可分割的组合，它可以最大程度地利用计算机来准确解决结构分析中所遇到的复杂问题，使其得到高效的计算。
matlab有限元计算
matlab有限元计算有限元计算是一种常用的数值计算方法，广泛应用于工程领域。
而Matlab作为一种强大的数值计算软件，提供了丰富的工具和函数，使得有限元计算变得更加简单和高效。
有限元计算是一种将连续问题离散化为有限个简单子问题的方法。
它将复杂的连续问题转化为离散的有限元网格，然后通过求解每个单元上的方程，最终得到整个问题的解。
有限元计算可以用于求解结构力学、流体力学、热传导等各种物理问题。
在Matlab中进行有限元计算，首先需要构建有限元模型。
有限元模型由节点和单元组成，节点是问题的离散点，单元是连接节点的基本单元。
在Matlab中，可以使用函数如meshgrid、linspace等来生成节点坐标，使用函数如delaunay、trimesh等来生成单元。
然后，需要定义问题的边界条件和加载条件。
边界条件是指在问题的边界上给定的约束条件，加载条件是指在问题中施加的外部力或位移。
在Matlab中，可以使用函数如boundary、findEdges等来定义边界条件，使用函数如force、displacement等来定义加载条件。
接下来，需要定义问题的材料性质和单元特性。
材料性质是指问题中所使用的材料的力学性质，单元特性是指单元的几何形状和材料性质。
在Matlab中，可以使用函数如materialProperties、elementProperties等来定义材料性质和单元特性。
然后，需要建立有限元方程。
有限元方程是通过对每个单元上的方程进行组装得到的整体方程。
在Matlab中，可以使用函数如stiffnessMatrix、loadVector等来建立有限元方程。
最后，需要求解有限元方程。
在Matlab中，可以使用函数如solve、eigs等来求解有限元方程。
求解得到的结果可以用于分析问题的应力、位移、变形等。
除了上述基本步骤，Matlab还提供了丰富的后处理工具和函数，用于可视化和分析有限元计算的结果。
有限元课程设计matlab
有限元课程设计matlab一、课程目标知识目标：1. 学生能理解有限元分析的基本原理，掌握运用MATLAB进行有限元建模和求解的基本步骤。
2. 学生能够运用MATLAB软件进行简单物理场的有限元模拟，并解释模拟结果。
3. 学生掌握如何将实际问题抽象为有限元模型，并能够运用MATLAB进行模型参数的设定和调整。
技能目标：1. 学生能够独立操作MATLAB软件，进行有限元模型的构建和求解。
2. 学生能够通过MATLAB编程实现有限元模型的自动化处理，包括前处理、求解和后处理。
3. 学生通过解决实际问题，提高数值分析能力和计算机应用能力。
情感态度价值观目标：1. 学生培养对科学研究的兴趣，特别是在工程计算和仿真领域。
2. 学生通过解决实际问题，体会数学和工程结合的美，增强对工程问题的探究欲望。
3. 学生通过团队合作解决问题，培养协作精神和解决问题的能力。
本课程针对高年级本科生或研究生，他们具备一定的数学基础和编程能力。
课程性质偏重实践，旨在通过MATLAB这一工具将有限元理论应用于具体问题的求解。
课程目标旨在使学生不仅掌握理论知识，而且能够实际操作，将理论知识转化为解决实际问题的技能。
通过课程学习，学生应能够将所学知识应用于未来的学术研究或工程实践中。
二、教学内容1. 有限元方法基本原理回顾：包括有限元离散化、单元划分、形函数、刚度矩阵和载荷向量等概念。
- 教材章节：第二章 有限元方法基础2. MATLAB编程基础：介绍MATLAB的基本操作、数据结构、流程控制、函数编写等。
- 教材章节：第三章 MATLAB编程基础3. MATLAB中的有限元工具箱使用：学习如何使用MATLAB内置的有限元工具箱进行建模和求解。
- 教材章节：第四章 MATLAB有限元工具箱介绍4. 有限元模型构建与求解：结合实际问题，学习如何构建有限元模型，并进行求解。
- 教材章节：第五章 有限元模型构建与求解5. 实例分析与上机操作：通过案例分析，让学生实际操作MATLAB软件，解决具体的有限元问题。
matlab 有限元法
matlab 有限元法
Matlab中的有限元法（Finite Element Method，FEM）是一种常用的数值分析方法，用于模拟和解决包括结构力学、热传导、流体力学等问题。
它将连续介质划分为离散的有限单元，通过建立数学模型和使用近似解法来求解。
下面是一般步骤来使用Matlab进行有限元分析：
1. 剖分网格：将要模拟的连续介质划分为离散的有限单元（如三角形或四边形元素）。
2. 建立数学模型：根据具体问题的物理方程或导引方程，建立线性或非线性的方程模型。
3. 施加边界条件：确定并施加边界条件，如位移、载荷或约束等。
4. 组装刚度矩阵和载荷向量（Assembly）：通过元素刚度矩阵的组装，得到总系统的刚度矩阵和载荷向量。
5. 求解方程：通过求解总系统的线性方程组，得到未知位移或其他需要的结果。
6. 后处理结果：对求解结果进行可视化或分析，如绘制应力分布、位移云图、应变曲线等。
Matlab提供了丰富的工具箱和函数，用于各种结构和物理问题的有限元分析，例如Partial Differential Equation Toolbox（部分微分方程工具箱）和Structural Analysis T oolbox（结构分析工具箱），其中包含了常用的有限元分析函数和设置界面。
另外，Matlab还支持用户自定义编程，允许使用脚本或函
数来实现特定的有限元算法。
总之，通过Matlab的有限元分析工具和编程能力，可以方便地进行各种结构和物理问题的数值分析和模拟。
计算力学(有限元)matlab编程大作业(空间网架)
逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 输入结构约束总数：11 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：3 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：3 按照编号逐一输入结构约束编号（先输入节点号） ） ：4 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：4 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：3 输入施加在结构上的总的力的个数：3 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：-1 输入力的大小：1000 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：1
2012年-有限元作业-matlab编程实现有限元求解简单结构位移及应力
2012年有限元作业对于三角形单元有： 1、形函数⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎣⎡⎥⎦⎤=⎩⎨⎧⎭⎬⎫k k j j i i k jik j i k j i k j ik j ik j i v u v u v u c c c b b b a a a c c c b b b a a a y x y x A v u 00000000000000001000000121 其中：j m m i m m j j i y x y x y x y x a -==m j mji y y y y b -=-=11j m m j i x x x x c -==11 kkj ji iy x y x y x A 1112=相乘后得：()])()()[(21,k k k k j j j j i i i i u y c x b a u y c x b a u y c x b a A y x u ++++++++=()])()()[(21,k k k k j j j j i i i i v y c x b a v y c x b a v y c x b a Ay x v ++++++++=也可写成：()()kk j j i i k k j j i i v N v N v N y x v u N u N u N y x u ++=++=,,简写为：][}{q N v u =⎩⎨⎧⎭⎬⎫其中}}{{Tk kj j i iv u v u v uq =Ay c x b a N A y c x b a N A y c x b a N k k k i j j j i i i i i 2/)(2/)(2/)(++=++=++= 2、应变有弹性力学知道：xvy u y v x u xy y x ∂∂+∂∂=∂∂=∂∂==γεε,,，可得 }{}{}{q B v u v u v u c b c b c b c c c b b b A u b u b u b v c v c v c v c v c v c u b u b u b A v u x yy xy v x u y v x u k k j j i i k kjjiik j i k j ik k j j i i k k j j i i k k j j i i k k j j i i =⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤=⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤+++++++++=⎩⎨⎧⎭⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧⎪⎪⎪⎭⎪⎪⎪⎬⎫∂∂∂∂∂∂∂∂=⎪⎪⎪⎩⎪⎪⎪⎨⎧⎪⎪⎪⎪⎪⎪⎪⎭⎫∂∂+∂∂∂∂∂∂=000000212100ε3、应力}{][}{[][]{}q B D D ==εσ[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=2100010112μμμμE D 4、刚度矩阵[]()⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+-+-+-+-=s r s r sr s r s r s r s r s r rs b b c c cb bc b c b b c c b b A Et K 21212121142μμμμμμμ 题目：如图所示是一平面梁。
有限元的matlab编程
本程序的网架位移求解函数附在主程序后面，主程序运行时调用该函数。
完整ppt
14
几何建模 定义荷载
用自定义输入
加下划线的为单元编号
完整ppt
5
形成等效荷载列阵
f=[0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a];%每个节点两个自由度，a 为之前输入的节点力
集成总刚：
获取单元两端节点坐标
xi = Node( Element( ie, 1 ), 1 ) ;%ie为单元号，以下相同 yi = Node( Element( ie, 1 ), 2 ) ; xj = Node( Element( ie, 2 ), 1 ) ; yj = Node( Element( ie, 2 ), 2 ) ;
节点号：',num2str(Element(ie,2)),' 轴力:',num2str(nodef(1))] ) ;
完整ppt
12
end
例二：网架
完整ppt
13
思路分析
几点说明
网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成 网架的基本单元有三角锥，三棱体，正方体，截头四角锥等。鉴于网架的形 式较多，本程序提供一种通用的网架输入方法，但录入较为繁琐，同时提供 一种正放四角锥网架的简易输入方法作为典型。
有限元编程示例
完整ppt
1
例一：桁架
题目描述:
如下图所示的平面桁架，杆件长度、弹性模量、截 面积以及所受节点力P的大小可以自行定义。求节 点位移及杆件轴力。
matlab有限元编程
是的，MATLAB可以用于有限元编程，用于解决各种结构、固体力学、热传导、电磁场等物理问题的数值模拟。
有限元法是一种数值方法，用于求解偏微分方程，特别适用于复杂的几何形状和边界条件的问题。
以下是一些使用MATLAB进行有限元编程的基本步骤：
1. 建立几何模型：首先，需要定义模型的几何形状和边界条件。
这包括确定节点和单元（元素）的位置。
2. 网格生成：将模型划分为离散的节点和单元。
这些节点和单元构成了有限元网格，可以使用MATLAB中的函数或者专用的网格生成工具来实现。
3. 定义材料属性和载荷：为每个单元分配适当的材料属性，如弹性模量、密度等。
同时，还要定义在结构上施加的力或边界条件。
4. 组装刚度矩阵和载荷向量：根据有限元法原理，将每个单元的局部刚度矩阵组装成全局刚度矩阵，同时将载荷向量组装成全局载荷向量。
5. 施加边界条件：根据边界条件，在全局刚度矩阵和载荷向量中施加约束条件。
6. 求解方程：通过求解线性方程组，得到节点的位移和其他所需的结果。
7. 后处理：根据求解结果，进行结果的可视化和分析。
可以绘制应力、应变分布图，计算位移、反应力等。
MATLAB提供了丰富的工具箱和函数来进行有限元编程，如Partial Differential Equation Toolbox（偏微分方程工具箱）、Finite Element Analysis Toolbox（有限元分析工具箱）等，可以大大简化有限元编程的过程。
在使用这些工具时，可以参考MATLAB的官方文档和例子，以及相关的有限元理论和方法。
matlab 有限元 编程
MATLAB 是一种高级编程语言和交互式环境，常用于数值计算、数据分析、算法开发、数据可视化以及应用开发等。
在有限元分析（FEA）中，MATLAB 可以用来编写和运行有限元程序。
以下是一个简单的 MATLAB 有限元分析编程示例：这个例子假设你正在解决一个一维拉普拉斯方程，其形式为 -d^2u/dx^2 = f，在区间 [0, 1] 上。
```matlab% 参数定义L = 1; % 长度h = 0.01; % 步长N = L/h; % 元素数量x = linspace(0, L, N+1); % x 向量% 创建有限元网格nodes = x(1:N+1);elements = [nodes(2:end-1) nodes(1:end-1) nodes(2:end)]; % 定义载荷向量 ff = sin(pi*x);% 定义刚度矩阵 A 和载荷向量 FA = zeros(N, N);F = zeros(N, 1);for i = 1:Nfor j = 1:i-1A(i, j) = A(j, i) = h^2/(6*(x(i)-x(j))^2);endF(i) = -f(i)*h/2;end% 使用 MATLAB 的线性方程求解器u = A\F;% 绘制结果plot(x, u);xlabel('x');ylabel('u');title('有限元解');```这个代码首先定义了问题的参数，然后创建了一个有限元网格。
然后，它定义了刚度矩阵 A 和载荷向量 F。
最后，它使用 MATLAB 的线性方程求解器来求解方程，并绘制结果。
请注意，这只是一个非常简单的例子。
在实际的有限元分析中，问题可能会更复杂，并且需要更多的编程工作。
空间框架静力计算源程序(程序结构力学大作业)
主要为一些工具函数 变量定义，单元属性分析 矩阵求解 整体控制模块
1：Lxz_Tools.f90 module Lxz_Tools implicit none integer (kind(1)),parameter ::ikind=(kind(1)) integer (kind(1)),parameter ::rkind=(kind(0.D0)) real (rkind), parameter :: Zero=0.D0,One=1.D0,Two=2.D0,Three=3.D0, & & Four=4.D0,Five=5.D0,Six=6.D0,Seven=7.D0,Eight=8.D0,Nine=9.D0, & & Ten=10.D0 contains function matinv(A) result (B) real(rkind) ,intent (in)::A(:,:) !real(rkind) , allocatable::B(:,:) real(rkind) , pointer::B(:,:) integer(ikind):: N,I,J,K real(rkind)::D,T real(rkind), allocatable::IS(:),JS(:) N=size(A,dim=2) allocate(B(N,N)) allocate(IS(N));allocate(JS(N)) B=A do K=1,N D=0.0D0 do I=K,N do J=K,N if(abs(B(I,J))>D) then D=abs(B(I,J)) IS(K)=I JS(K)=J
计算力学大作业报告
计算力学基础大作业报告航21班苏浩尹思凡一、求解流程说明本次作业的求解过程主要分为三个步骤：1. 前处理网格划分——采用ABAQUS软件进行画网格，并输出前处理文件。
形成前处理文件——利用matlab编程，形成inmesh文件。
2. 有限元计算采用上机课上提供的FEATP程序计算，程序略有修改，使其能够计算单元和节点数目更多的算例。
本次大作业所有的图中str1~3,分别代表x（r）、y(theta)方向的正应力，剪力。
3. 后处理采用Tecplot软件画图，得到位移场和应力场的云图。
二、所编程序及功能简介二、具体算例1.悬臂梁问题1.1 问题描述图1所示悬臂梁，一端固支，承受集中力 P=1000N，梁的长度为 100mm，高度为10mm，厚度为1mm，材料弹性模量为2×1011Pa，泊松比为0.32，假定在变形过程中界面始终保持完整，利用有限单元法计算位移场分布，并利用求得的位移场计算结点应力，讨论与理论解的误差。
分别使用不同类型不同阶次单元进行计算，比较结果精确度，单元收敛性等。
1.2 解析计算采用材料力学方法进行计算，建立主轴坐标系，根据材料力学的推导，有如下公式：()0()z x zy y z xy z M x y I Q S y I tσστ=-==-（1.1）其中()z M x 是z 轴正方向的弯矩，y Q 是方向为y 轴正方向的剪力，z I 是梁对z 方向的惯性矩，()z S y 是横截面对z 轴的静面矩，t 为厚度。
由题目中的数据，可以得到：33422223221000()()1000(100)()10125012123111()()()(25)2242y z z yy z h h Q P NM x P L x x N mm h t I mm h S y ytdy ty t y y mm --=-=-=--=-⋅-⋅⨯======-=-⎰ （1.2）将（1.2）中的公式带入（1.1）中，可以得到:()212(100)06(25)()x y xyx y MPa y MPa σστ⎧=-⎪⎪=⎨⎪=-⎪⎩ (1.3) 材料力学中的挠度公式和转角公式：22()z z z M x d EI dx d dxννθ==(1.4) 而x 方向的位移：(,)()z u x y x y θ=- (1.5)将（1.2）中的计算结果和材料参数带入，并对式（1.4）进行积分，得到：521523121610(100)21610(50)6d x x c dx x x c x c νν--=-⨯⨯-+=-⨯⨯-++ （1.6）位移边界条件：0,0,0dvx dxν=== 得到位移场：525231610(100)21610(50)6u x x y x x ν--=-⨯⨯-=-⨯⨯- （1.7） 1.3 问题讨论（1）位移边界条件：固定位移：左边界——中点 u=0 v=0 其余点 u=0 （2）载荷条件：右上角点1.4 用不同类型的单元求解位移场和应力场 1.4.1 三角形三节点单元 （1）网格划分采用三种不同疏密度的网格计算：（2）计算结果这里仅选取网格密度最大的第三种情况的结果 位移场,u vσστ应力场,,x y xy1.4.2 四边形四节点单元（1）网格划分采用四种不同疏密度的网格计算：（2）计算结果位移场,u vσστ应力场,,x y xy1.4.3 四边形八节点单元（1）网格划分采用三种疏密度不同的网格计算（2）计算结果位移场,u vσστ应力场,,x y xy1.5 单元精度与单元收敛性1）考虑平面应力问题，理论上v 的挠度位移：3635PL PL v EI GA=+ （1.8）代入数据：P=1000N ，L=100mm E=2×105MPa G=E/(2（1+ν）)=75757.6MPa 得到理论挠度：3353610001006100010020.1584135575757.6110321011012PL PL v mm EI GA ⨯⨯⨯=+=+=⨯⨯⨯⨯⨯⨯⨯⨯以（2）再对应力场进行检验以（50,10）处的正应力值x σ作为收敛对象，理论()121005053000x MPa σ=⨯-⨯=上表及上图中，三角形单元T3的节点数分别为63、205、1111和4221，四边形单元Q4和Q8的单元数相同，节点数分别为Q4：63、205、1111、4221;Q8:165,569,3221。
有限元法MATLAB作业
例题：如图所示的受力均匀分布载荷作用的薄平板结构，将平板离散成两个线性三角元，如图所示，假定2210,0.3,0.025,3000/m E GPa v t m w KN ====。
求 （1）该结构的整体刚度矩阵。
（2）节点2和节点3的水平位移和垂直位移。
（3）节点1和节点4的支反力。
（4）每个单元的应力。
（5）每个单元的主应力和主应力方向角。
图2. 用双线性三角元离散化的薄木板 解：（1）离散化域我们将平板分为两个单元，4个节点，如图2所示，分布载荷的总作用力平均分给节点2和节点3，由于结构是薄平板，所以假定其属于平面应力情况。
MATLAB 中采用的计算单位是KN 和m 。
表1给出了该题的单元连通性。
（2）单元刚度矩阵通过调用MATLAB 的LinearTriangleElementStiffness 函数，得到两个单元矩阵k 1和k 2，每个矩阵都是66⨯的。
k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)的源程序：function k1= gangdujuzhen( f,E,NU,t,xi,yi,xj,yj,xm,ym,p ) %UNTITLED4 此处显示有关此函数的摘要% 此处显示详细说明k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)图1. 薄平板结构k1= f(E,NU,t,xi,yi,xj,yj,xm,ym,p);>>k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)k1 =1.0e+06 *2.0192 0 0 -1.0096 -2.0192 1.00960 5.7692 -0.8654 0 0.8654 -5.76920 -0.8654 1.4423 0 -1.4423 0.8654-1.0096 0 0 0.5048 1.0096 -0.5048-2.0192 0.8654 -1.4423 1.0096 3.4615 -1.87501.0096 -5.7692 0.8654 -0.5048 -1.8750 6.2740k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)的源程序：function k2= gangdujuzhen2( f,E,NU,t,xi,yi,xj,yj,xm,ym,p )%UNTITLED3 此处显示有关此函数的摘要% 此处显示详细说明k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)k2= f(E,NU,t,xi,yi,xj,yj,xm,ym,p);>>k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)k2 =1.0e+06 *1.4423 0 -1.4423 0.8654 0 -0.86540 0.5048 1.0096 -0.5048 -1.0096 0-1.4423 1.0096 3.4615 -1.8750 -2.0192 0.86540.8654 -0.5048 -1.8750 6.2740 1.0096 -5.76920 -1.0096 -2.0192 1.0096 2.0192 0-0.8654 0 0.8654 -5.7692 0 5.7692（3）集成整体刚度矩阵⨯的，因此为了得到整体刚度矩阵K，我们由于结构有4个节点，所以刚度矩阵是88⨯的零矩阵。
计算力学编程大作业
定位向量如下 namta{1}=[4 5 2 1]; namta{2}=[5 6 3 2]; namta{3}=[7 8 5 4]; namta{4}=[8 9 6 5]; namta{5}=[10 11 8 7]; namta{6}=[11 12 9 8];
三、 控制方程
图示平面问题的控制方程为 U - b - ST = 0 因为只是计算静力问题，所以不考虑惯性项 b ST = 0
1 0 0
D
E
1
2
0
0
1 0 0
0 0 0
0
0
1
2
应变位移矩阵
N a
x
B
a
0
0
N a y
0
N
a
y
0
N a x
计算单元刚度和单元荷载
单元刚度矩阵：
采用高斯积分，x、y 向各取两个点求得单元的刚度矩阵
1
1
i
3
3
i 1
1
边界荷载
1) 对于分布力荷载
3
如单元 1 和单元 2
因为 N1 N2 在上边界 1 的取值为零
形函数如下 N {1} =1/4*(1-x/(b/4)) * (1-y/(a/2)); N {2} =1/4*(1+x/(b/4)) * (1-y/(a/2)); N {3} =1/4*(1+x/(b/4)) * (1+y/(a/2)); N {4} =1/4*(1-x/(b/4)) * (1+y/(a/2)); 与材料性质有关的矩阵 D 如下（先考虑平面应力问题：薄板）
所以 f11
f21
0 0
只对节点 3、4 上的荷载进行考虑和计算
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
M=jiedianlianjie(n);
L=l(M,N,n);
%n*n 的矩阵，用于存储 n 个节点的连接关系，连接则对应元素为 1，否则为 0 %n*n 矩阵，用于存储 I,j 节点连接所得的杆件的长度
Y=ganjiangeshu(L,n); %数量值，表示杆件的个数 t=T(L,N,n); %3n*3n 矩阵，用于存储以上各杆件的坐标转换矩阵 S=dangangjihe(E,A,n,Y,L,t); %6*6y 矩阵，用于存储 y 根杆件的各自的单刚矩阵 K=zonggang(E,A,L,t,n); %3n*3n 矩阵，描述总刚 P=zhiling(K); %总刚矩阵主元根据约束条件置零 Q=jia1(P,n); %总刚矩阵置零主元加 1 z=lixiangliang(n); %由输入力组成的列向量 p=inv(Q)*z; %求解位移 F=neili(n,Y,S,L,p); %由位移及单刚矩阵求得各杆件内力 三、命令执行输入参数及计算结果： （5 节点九单元，实际模型见下图）
F = (KN)
0 0 -230 0 0 230 0 0 0 0 0 0 0 0 0 0 0 0 0 230 -230 0 -230 230 6.82E-13 0 -6.82E-13 -6.82E-13 0 6.82E-13 1000 1000 -1000 -1000 -1000 1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
计算力学编程大作业（网架结构）
一、 编程算法：
j节点
i节点
图1
1.总刚矩阵形成： 根据图 1，对于任意的 i,j 节点，不论原来的连接状况如何，只要 I,j 两节点有了联系， 则在总刚矩阵中对应的位置需要加上有 ij 杆的产生带来的刚度贡献，贡献具体可表现 为在 Kii 处需要加上有 ij 杆件贡献的 K11（经过坐标变换后） ，Kjj 处需要加上 ij 杆件贡 献的 K22（经过坐标变换后） ，Kij 处需要加上有 ij 杆件贡献的 K12（经过坐标变换后） ， Kji 处需要加上有 ij 杆件贡献的 K21（经过坐标变换后）,由此可以直接形成总刚矩阵。 2.非线性实现（分段力，存在系统误差） 假定杆件一端发生了总位移超过了其自身杆长的 1/1000（此参数可调）为大位移，需 要考虑由于其端点发生位移而导致结构刚度变化的影响，将由原来不考虑非线性的结 构各节点位移算出，求出各节点相对于自身杆件发生的位移，取其中最大值，将之设 定为 1/1000，由此可得出需要考虑非线性的临界荷载，将原荷载分成 n 份（n 取整） 临界荷载，逐级加载在结构上，每一次加载都考虑结构节点位移对结构刚度的影响。 由此计算出结构的最终位移。 3.牛顿拉普森迭代法实现非线性 利用牛顿拉普森迭代法，在原来线性的基础上，计算出在加上变形的基础上求出的结 构新的刚度矩阵K1 ， 由∆F = F − K1 ∗ p得到差量∆F， 再将结构在∆F作用下的位移求出， 并求出发生该位移后结构新的刚度，重复以上步骤，直到前后两次节点变形最大值小 于 1e-8，循环结束，输出结果。 二、命令窗口输入命令（不考虑非线性） ： n=input('请输入节点个数:'); E=input('请输入 E 的大小:'); A=input('请输入 A 的大小:'); N=dianzuobiao(n); %3*n 的矩阵，用于存储 n 个节点的坐标
逐一输入力对应的向量中的坐标值：0 输入力的大小：1230 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：1 逐一输入力000
四、计算结果： p = (m)
0.000000000000000 0.000000000000000 0.000000000000000 0.002272807092008 0.000440269119346 0.000115000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.002157807092008 0.000000000000000 0.000000000000000 0.000000000000000
请输入节点个数:5 请输入 E 的大小:2e8 请输入 A 的大小:0.01 从左往右从下往上编号输入第 i 节点 x 坐标：0 从左往右从下往上编号输入第 i 节点 y 坐标：0 从左往右从下往上编号输入第 i 节点 z 坐标：0 从左往右从下往上编号输入第 i 节点 x 坐标：0 从左往右从下往上编号输入第 i 节点 y 坐标：0 从左往右从下往上编号输入第 i 节点 z 坐标：1 从左往右从下往上编号输入第 i 节点 x 坐标：0 从左往右从下往上编号输入第 i 节点 y 坐标：1 从左往右从下往上编号输入第 i 节点 z 坐标：0 从左往右从下往上编号输入第 i 节点 x 坐标：1 从左往右从下往上编号输入第 i 节点 y 坐标：0 从左往右从下往上编号输入第 i 节点 z 坐标：0 从左往右从下往上编号输入第 i 节点 x 坐标：1 从左往右从下往上编号输入第 i 节点 y 坐标：1 从左往右从下往上编号输入第 i 节点 z 坐标：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0
五、有限元软件验算 模型建立：
模型计算结果：
软件内力输出： （此处内力是 x,y,z 方向的合力） 1 节点荷载 0.000000 0.000000 2 节点荷载 0.000000 0.000000 3 节点荷载 0.000000 0.000000 4 节点荷载 -325.269119 -325.269119 5 节点荷载 230.000000 230.000000 6 节点荷载 0.000000 0.000000 7 节点荷载 0.000000 0.000000 8 节点荷载 0.000000 0.000000 9 节点荷载 -1732.050808 -1732.050808 软件位移输出： 1 节点荷载 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2 节点荷载 0.002273 0.000440 0.000115 0.000000 0.000000 0.000000 3 节点荷载 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4 节点荷载 0.000000 0.000000 -0.002158 0.000000 0.000000 0.000000 5 节点荷载 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 完全一致。 六、窗口命令输入（分段力考虑非线性，在上文窗口命令的基础上输入） a=maxbianxing(L,p,n,Y); %找出所有杆件中（位移/杆长）最大值 [f,nmax]=fenduanli(z,a); %求出单元加载力向量以及需要加载的次数（nmax） [d,G,LL]=feixianxing(E,A,f,Q,N,n,M,P,nmax); %求出考虑非线性下的位移 d，并输出最终变形后的坐
逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 输入结构约束总数：11 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：3 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：3 按照编号逐一输入结构约束编号（先输入节点号） ） ：4 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：4 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：3 输入施加在结构上的总的力的个数：3 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：-1 输入力的大小：1000 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：1