产品
学校上课的作业，之前给过不少同学，明天结课考试也不知道还能复习啥，整理一下代码吧。非常按部就搬刚度矩阵叠加出来的……个人认为写得还可以的地方： 1. 利用了Matlab里稀疏矩阵sparse()的特性构造整体刚度矩阵，可以少写几个字 2. 查书，用边界条件解决了奇异的问题，也是这招帮几个同学de了bug 3. 能解决大部分平面刚架有限元问题——本来是梁单元改出来的代码，改回去也可以吧！ 写的不太好的地方有： 1. 不怎么写代码了，Matlab很多特性不清楚，感觉分块叠加整体刚度矩阵的部分可以写得更简洁，用不着四个循环== 2. 没有画图，只是单纯算了结果。毕竟老师不要图嘛。
代码用了两个例子测试，分别是《有限元分析——ANSYS理论与应用（第四版）》例4.5（例1）和北航材料力学第二册刚度矩阵章节的例题（例2.1与2.2），自己多分了几个节点，电子版有学号水印就不发了。下面是代码：
% bilibili没有Matlab代码块，晕！
F=[0 0 0 0 0 0 0 -4000 80000 0 0 0 0 0 0 0 0 0]' % 节点载荷
la=[10*12/2 10*12/2 3*12 3*12 3*12]; % 单元长度
Ea=[30000000 30000000 30000000 30000000 30000000]; % 单元弹性模量
Ia=[204 204 204 204 204];% 单元惯性矩
Aa=[7.65 7.65 7.65 7.65 7.65];% 单元面积
sitaa=[0 0 3*pi/2 3*pi/2 3*pi/2]% 单元转角
flag=[1 2 3 16 17 18];% 0位移标志
n=6;% 单元刚度矩阵大小
node=[1 2
2 3
3 4
4 5
5 6]% 哪两个节点间有杆件
a=asseml(Aa,la,Ea,Ia,flag,sitaa,F,n,node)
F=[0 0 0 20000 -20000 -40000000 0 0 0 0 0 0]'
la=[1000 1000 1000]
Ea=[210000 210000 210000]
Aa=[2000 2000 2000]
Ia=[3000800 3000000 3000000]
sitaa=[0 0 3*pi/2]
flag=[1 2 3 7 8 9 10 11 12]
node=[1 2
2 3
2 4]
a=asseml(Aa,la,Ea,Ia,flag,sitaa,F,n,node)
F=[0 0 0 0 0 0 20000 -20000 -40000000 0 0 0 0 0 0]'
la=[500 500 1000 1000]
Ea=[210000 210000 210000 210000]
Aa=[2000 2000 2000 2000]
Ia=[3000800 3000000 3000000 3000000]
sitaa=[0 0 0 3*pi/2]
flag=[1 2 3 10 11 12 13 14 15]
node=[1 2
2 3
3 4
3 5]
a=asseml(Aa,la,Ea,Ia,flag,sitaa,F,n,node)
function kt=matx(A,l,E,I,sita)
kk=[A*E/l 0 0 -A*E/l 0 0
0 12*E*I/(l^3) 6*E*I/(l^2) 0 -12*E*I/(l^3) 6*E*I/(l^2)
0 6*E*I/(l^2) 4*E*I/l 0 -6*E*I/(l^2) 2*E*I/l
-A*E/l 0 0 A*E/l 0 0
0 -12*E*I/(l^3) -6*E*I/(l^2) 0 12*E*I/(l^3) -6*E*I/(l^2)
0 6*E*I/(l^2) 2*E*I/l 0 -6*E*I/(l^2) 4*E*I/l];
T=[cos(sita) sin(sita) 0 0 0 0
-sin(sita) cos(sita) 0 0 0 0
0 0 1 0 0 0
0 0 0 cos(sita) sin(sita) 0
0 0 0 -sin(sita) cos(sita) 0
0 0 0 0 0 1];
kt=(T')*kk*(T);
end
function ans=asseml(Aa,la,Ea,Ia,flag,sitaa,F,n,node)
is=[];
js=[];
vs=[];
for t=1:length(la)
x=node(t,1);
y=node(t,2);
kt=matx(Aa(t),la(t),Ea(t),Ia(t),sitaa(t));
for i=1:3
for j=1:3
is=[is 3*(x-1)+i];
js=[js 3*(x-1)+j];
vs=[vs kt(i,j)];
end
for j=4:6
is=[is 3*(x-1)+i];
js=[js 3*(y-1)+j-3];
vs=[vs kt(i,j)];
end
end
for i=4:6
for j=1:3
is=[is 3*(y-1)+i-3];
js=[js 3*(x-1)+j];
vs=[vs kt(i,j)];
end
for j=4:6
is=[is 3*(y-1)+i-3];
js=[js 3*(y-1)+j-3];
vs=[vs kt(i,j)];
end
end
end
kk=full(sparse(is,js,vs))
for i=flag
for j=1:(3*(length(la)+1))
kk(i,j)=0;
end
kk(i,i)=1;
end
ans=kk\F;
end运行结果：
例1
例2.1
例2.2
附上一份没有用应变能的手写推导，字很丑，有空拿tex整理一下：
matepad paper挺好用的，已经一两个月没开过ipad了……
155-2731-8020
