节点、单元编号
只计算到节点位移,应力没计算
%% 《有限单元法》--王勖成,习题2.12 MATLAB 程序
% 仅供参考
function fem()
P=[-10 0 10 0 10 0 -10 0]';%结构节点载荷列向量
[K,P]=xQiYi(P,make_K());
solu(K,P);
%% 存储各单元的节点 x、y 值
function [a,b]=get_xy(m)
x=[0 12 0;
0 12 12;];
y=[4 0 0;
4 4 0];
a=x(m,:);
b=y(m,:);
end
%% 求解并显示
function solu(K,P)
n=1000;
JieDianWeiYi=inv(K)*P
% YingLi=WeiYi(3,1)/12*2e5 %x向拉、压时可用
plot([0 12 12 0],[4 4 0 0],'b');hold on;xlim([-1 14]);ylim([-4 8]);
plot(JieDianWeiYi([1 3 5 7],1)'*n+[0 12 12 0],JieDianWeiYi([2 4 6 8],1)'*n+[4 4 0 0],'r');
end
%% 用“对角元素改 1 法” ,引入位移约束并使 总的K 不奇异
function [b,a]=xQiYi(P,K)
%---------------------------------------------------------------------------------%
d=[0 0 1 1 1 1 0 0]';%左侧两结点铰支
% 共 4 个节点,顺序为 :左上角开始1 ,顺时针排序
% 共 8个位移分量,[dx1 dy1 dx2 dy2 ... ]
% 位移为 0的分量,对应项填 0,其他填非 0 值
% 例:右侧两点铰支,则为 d=[ 1 1 0 0 0 0 1 1]'
for i=1:1:8
if d(i,1)==0
for j=1:1:8
K(i,j)=0;
K(j,i)=0;
end
P(i,1)=0;
K(i,i)=1;
end
end
a=P;
b=K;
end
%% 组装总的结构刚度矩阵
function K=make_K()
K=[su([1 1],[1 1]) su(0,[1 2]) su([1 2],[1 3]) su([1 3],0);
su(0,[2 1]) su(0,[2 2]) su(0,[2 3]) su(0,0);
su([2 1],[3 1]) su(0,[3 2]) su([2 2],[3 3]) su([2 3],0);
su([3 1],0) su(0,0) su([3 2],0) su([3 3],0)];
end
%% 分块计算矩阵的和
function b=su(a,b)
K1=get_K(1);
K2=get_K(2);
if a==0
K1=zeros(6,6);
a=[1,1];
end
if b==0
K2=zeros(6,6);
b=[1,1];
end
h11=a(1)*2-1; h12=a(1)*2; l11=a(2)*2-1; l12=a(2)*2;
h21=b(1)*2-1; h22=b(1)*2; l21=b(2)*2-1; l22=b(2)*2;
b=[K1(h11,l11)+K2(h21,l21) K1(h11,l12)+K2(h21,l22);
K1(h12,l11)+K2(h22,l21) K1(h12,l12)+K2(h22,l22)];
end
%% 计算 m 单元的 刚度矩阵
function K=get_K(m)
%-------------------------------------------------------%
v=0.3;
D=2e5/(1-v^2)*[ 1 ,v ,0 ;
v ,1 , 0 ;
0 ,0 ,(1-v)/2];
%弹性矩阵
[x,y]=get_xy(m);
B=get_B(x,y);
A=get_A(x,y);
K=B'*D*B*A;
end
%% 计算单元应变矩阵 B
function B=get_B(x,y)
[~,b,c]=get_abc(x,y);
B=zeros(3,6);
[x,y]=get_xy(1);
A=get_A(x,y);
for i=1:1:3
B(:,i*2-1:i*2)=[b(i) 0;
0 c(i);
c(i) b(i)];
end
B=B/(2*A);
end
%% 计算 ai bi ci
function [a,b,c]=get_abc(x , y)
a=zeros(2,3);
b=a;
c=a;
for l=0:1:2
i=l+1;
j=rem(l+1,3)+1;
m=rem(l+2,3)+1;
a(i)=det([x(j) y(j);
x(m) y(m)]);
b(i)= y(j) - y(m);
c(i)= -x(j) + x(m);
end
end
%% 计算单元面积
function A=get_A(x,y)
A=-det([1 x(1) y(1); 1 x(2) y(2) ; 1 x(3) y(3)])/2;
end
end