北航有限元编程大作业(Matlab)

1、%* 有 限 元 编 程 作 业 * y=1; while y=0 x=input('请选择单元类型：1 为二力杆单元；2 为平面刚架单元；3 为平面三角形单元n' %* 二 力 杆 单 元 * if x=1 disp('所选单元为：二力杆单元' %二力杆单元刚度矩阵的生成 %输入基本参数 ndof = 2;%节点自由度数 nn = input('输入节点总数n' nm = input('输入单元数n' L = input('输入各杆件的长度,单位为:mn' E = input('弹性模量,单位为:KN/m

2、2n' A = input('单元面积，单位为:m2n' Q = input('输入各单元倾角,单位为:°n' q = Q*pi/180;%将Q转化为弧度 a = cos(q; b = sin(q; JM = input('输入节点码数组'%将单元与节点进行对应 K2=zeros(nn*ndof; for e=1:nm%按单元数进行循环 k=A*E/L(e*a(e2,a(e*b(e;a(e*b(e,b(e2; K=k,-k;-k,k;%整体坐标中的单元刚度矩阵 K1=zeros(nn*ndof; m=JM(e,1;n=JM(e,

3、2;%确定单元刚度矩阵在总刚度矩阵中的位置 %组集总刚度矩阵 K1(2*m-1,(2*m-1=K(1,1; K1(2*m-1,(2*m=K(1,2; K1(2*m-1,(2*n-1=K(1,3; K1(2*m-1,(2*n=K(1,4; K1(2*m,(2*m-1=K(2,1; K1(2*m,(2*m=K(2,2; K1(2*m,(2*n-1=K(2,3; K1(2*m,(2*n=K(2,4; K1(2*n-1,(2*m-1=K(3,1; K1(2*n-1,(2*m=K(3,2; K1(2*n-1,(2*n-1=K(3,3; K1(2*n-1,(2*n=K(3,4; K1(2*n,(2*m-1

4、=K(4,1; K1(2*n,(2*m=K(4,2; K1(2*n,(2*n-1=K(4,3; K1(2*n,(2*n=K(4,4; K2=K2+K1; end disp('总刚度矩阵为:' K2 %- %组集载荷以及求解线性方程组 nf=input('输入节点载荷的个数n' NF=zeros(nf; F=input('输入已知节点载荷(单位为KN及对应行的数组(2维n' for i=1:nf for j=1:nf NF(i,j=K2(F(i,2,F(j,2; end end F1=zeros(nf,1; for i=1:nf F1(i=F(i,

5、1; end V=inv(NF*F1;%位移向量 U1=zeros(2*nn,1; for m1=1:nf U1(F(m1,2=V(m1; end disp('求得的位移为:' U1 end %*平 面 刚 架 单 元* if x=2 disp('所选单元为：平面刚架单元' %刚架单元总刚度矩阵的生成 %输入基本参数 ndof = 3;%节点自由度数 nn = input('输入节点总数n' nm = input('输入单元数n' L = input('输入各杆件的长度,单位：mn' E = input('

6、弹性模量，单位KN/m2n' I = input('输入惯性矩，单位：m4n' A = input('单元面积,单位：m2n' Q = input('输入各单元倾角，单位：°n' q = Q*pi/180;%将Q转化为弧度 a = cos(q; b = sin(q; JM = input('输入节点码数组n' K4=zeros(nn*ndof; for e=1:nm%按单元数进行循环 T=a(e,b(e,0,0,0,0;-b(e,a(e,0,0,0,0;0,0,1,0,0,0;0,0,0,a(e,b(e,0;0,

7、0,0,-b(e,a(e,0;0,0,0,0,0,1; K1=E*I/L(e3*A*L(e*L(e/I,0,0,-A*L(e*L(e/I,0,0;0,12,6*L(e,0,-12,6*L(e;0,6*L(e,4*L(e*L(e,0,-6*L(e,2*L(e*L(e;-A*L(e*L(e/I,0,0,A*L(e*L(e/I,0,0;0,-12,-6*L(e,0, 12,-6*L(e;0,6*L(e,2*L(e*L(e,0,-6*L(e,4*L(e*L(e; K2=T'*K1*T;%整体坐标中的单元刚度矩阵 K3=zeros(nn*ndof; m=JM(e,1; n=JM(e,2;%确定单

8、元刚度矩阵在总刚度矩阵中的位置 %组集总刚度矩阵 K3(3*m-2,(3*m-2=K2(1,1; K3(3*m-2,(3*m-1=K2(1,2; K3(3*m-2,(3*m=K2(1,3; K3(3*m-2,(3*n-2=K2(1,4; K3(3*m-2,(3*n-1=K2(1,5; K3(3*m-2,(3*n=K2(1,6; K3(3*m-1,(3*m-2=K2(2,1; K3(3*m-1,(3*m-1=K2(2,2; K3(3*m-1,(3*m=K2(2,3; K3(3*m-1,(3*n-2=K2(2,4; K3(3*m-1,(3*n-1=K2(2,5; K3(3*m-1,(3*n=K2(

9、2,6; K3(3*m,(3*m-2=K2(3,1; K3(3*m,(3*m-1=K2(3,2; K3(3*m,(3*m=K2(3,3; K3(3*m,(3*n-2=K2(3,4; K3(3*m,(3*n-1=K2(3,5; K3(3*m,(3*n=K2(3,6; K3(3*n-2,(3*m-2=K2(4,1; K3(3*n-2,(3*m-1=K2(4,2; K3(3*n-2,(3*m=K2(4,3; K3(3*n-2,(3*n-2=K2(4,4; K3(3*n-2,(3*n-1=K2(4,5; K3(3*n-2,(3*n=K2(4,6; K3(3*n-1,(3*m-2=K2(5,1; K3(

10、3*n-1,(3*m-1=K2(5,2; K3(3*n-1,(3*m=K2(5,3; K3(3*n-1,(3*n-2=K2(5,4; K3(3*n-1,(3*n-1=K2(5,5; K3(3*n-1,(3*n=K2(5,6; K3(3*n,(3*m-2=K2(6,1; K3(3*n,(3*m-1=K2(6,2; K3(3*n,(3*m=K2(6,3; K3(3*n,(3*n-2=K2(6,4; K3(3*n,(3*n-1=K2(6,5; K3(3*n,(3*n=K2(6,6; K4=K3+K4; end disp('总刚度矩阵为:' K4 %- nf=input('输入

11、节点载荷的个数n' NF=zeros(nf; F=input('输入已知节点载荷(单位为KN及对应行的数组(2维n' for i=1:nf for j=1:nf NF(i,j=K4(F(i,2,F(j,2; end end F1=zeros(nf,1; for i=1:nf F1(i=F(i,1; end V=inv(NF*F1;%位移向量 U2=zeros(3*nn,1; for m2=1:nf U2(F(m2,2=V(m2; end disp('求得的位移为:' U2 end %*三 角 形 单 元* if x=3 disp('所选单元为：三

12、角形单元' %三角形单元总刚度矩阵的生成 ndof =2;%节点自由度数 nn = input('输入节点数n' nm = input('输入单元数n' Z = input('输入节点坐标(单位为：mx= y=n' E = input('弹性模量值,单位为：108KN/m2n' u = input('泊松比n' t = input('输入单元的厚度，单位为mn' JM = input('输入单元对应节点编码的数组n' K3=zeros(ndof*nn;%总刚度矩阵初始值 fo

13、r e=1:nm%按单元编码循环 m=JM(e,1; n=JM(e,2; p=JM(e,3; D =E/(1-u2,u*E/(1-u2,0;u*E/(1-u2,E/(1-u2,0;0,0,0.5*E/(1+u; d1=1,Z(m,1,Z(m,2;1,Z(n,1,Z(n,2;1,Z(p,1,Z(p,2; d2=0.5*det(d1; B=0.5/d2*(Z(n,2-Z(p,2,0,(Z(p,1-Z(n,1;0,(Z(p,1-Z(n,1,(Z(n,2-Z(p,2;(Z(p,2-Z(m,2,0,(Z(m,1-Z(p,1;0,(Z(m,1-Z(p,1,(Z(p,2-Z(m,2;(Z(m,2-Z(n,2

14、,0,(Z(n,1-Z(m,1;0,(Z(n,1-Z(m,1,(Z(m,2-Z(n,2; K1=B*D*B'*d2;%整体坐标中的单元刚度矩阵 K2=zeros(ndof*nn;%扩展后的单元刚度矩阵 %组集总刚度矩阵 K2(2*m-1,(2*m-1=K1(1,1; K2(2*m-1,(2*m=K1(1,2; K2(2*m-1,(2*n-1=K1(1,3; K2(2*m-1,(2*n=K1(1,4; K2(2*m-1,(2*p-1=K1(1,5 ; K2(2*m-1,(2*p=K1(1,6; K2(2*m,(2*m-1=K1(2,1; K2(2*m,(2*m=K1(2,2; K2(2*

15、m,(2*n-1=K1(2,3; K2(2*m,(2*n=K1(2,4; K2(2*m,(2*p-1=K1(2,5; K2(2*m,(2*p=K1(2,6; K2(2*n-1,(2*m-1=K1(3,1; K2(2*n-1,(2*m=K1(3,2; K2(2*n-1,(2*n-1=K1(3,3; K2(2*n-1,(2*n=K1(3,4; K2(2*n-1,(2*p-1=K1(3,5; K2(2*n-1,(2*p=K1(3,6; K2(2*n,(2*m-1=K1(4,1; K2(2*n,(2*m=K1(4,2; K2(2*n,(2*n-1=K1(4,3; K2(2*n,(2*n=K1(4,4;

16、 K2(2*n,(2*p-1=K1(4,5; K2(2*n,(2*p=K1(4,6; K2(2*p-1,(2*m-1=K1(5,1; K2(2*p-1,(2*m=K1(5,2; K2(2*p-1,(2*n-1=K1(5,3; K2(2*p-1,(2*n=K1(5,4; K2(2*p-1,(2*p-1=K1(5,5; K2(2*p-1,(2*p=K1(5,6; K2(2*p,(2*m-1=K1(6,1; K2(2*p,(2*m=K1(6,2; K2(2*p,(2*n-1=K1(6,3; K2(2*p,(2*n=K1(6,4; K2(2*p,(2*p-1=K1(6,5; K2(2*p,(2*p=K1(6,6; K3=K3+K2;%总刚度矩阵 end disp('总刚度矩阵为:' K3 %- %组集载荷以及求解线性方程组 nf=input('输入节点载荷的个数n' NF=zeros(nf; F=input('输入已知节点载荷(单位为KN及对应行号的数组(2维n' for i=1:nf for j=1: