%****************** 有 限 元 编 程 作 业 *****************

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('输入各杆件的长度,单位为:m\n')

E = input('弹性模量,单位为:KN/m^2\n')

A = input('单元面积，单位为:m^2\n')

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(e)^2,a(e)*b(e);a(e)*b(e),b(e)^2];

K=[k,-k;-k,k];%整体坐标中的单元刚度矩阵

K1=zeros(nn*ndof);

m=JM(e,1);n=JM(e,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))=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,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('输入各杆件的长度,单位：m\n')

E = input('弹性模量，单位KN/m^2\n')

I = input('输入惯性矩，单位：m^4\n')

A = input('单元面积,单位：m^2\n')

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,0,0,-b(e),a(e),0;0,0,0,0,0,1];

K1=E*I/L(e)^3*[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,