fembeam2d.m

梁元有限元分析源码（用matlab语言编写而成）

% plane stress problem
clc;clear;
Lx=9;Ly=3;
%x=Lx;y=Ly;%单元范围[0,Lx]*[0,Ly]
Nx=12;Ny=6;%x,y方向划分单元数,Nx=Ny
Nodes=(Nx+1)*(Ny+1);%节点总数
Enum=Nx*Ny;
A=(Lx*Ly)/(2*Nx*Ny);%单元面积
E=2e10;
mu=0.167;
t=1;
q=10*Lx;
qe=q/Nx;

D=E*[1,mu,0;mu,1,0;0,0,(1-mu)/2]/(1-mu*mu);
xe1=[Lx/Nx;0;0];ye1=[Ly/Ny;Ly/Ny;0];%上三角，单数单元
xe2=[0;Lx/Nx;Lx/Nx];ye2=[0;0;Ly/Ny];%下三角，双数单元
b1=zeros(3,1);c1=zeros(3,1);
b2=zeros(3,1);c2=zeros(3,1);
for i=1:3
    if(i==1)	j=2;m=3;                                       %i
		else if(i==2)	j=3;m=1;                               %%  %
		    else if(i==3)	j=1;m=2;                           %%    %
                 end                                           %%%%%%%%%%% m
             end                                                %j
    end
    b1(i)=ye1(j)-ye1(m);
    c1(i)=-xe1(j)+xe1(m);
    b2(i)=ye2(j)-ye2(m);
    c2(i)=-xe2(j)+xe2(m);
end

B1=[b1(1),0,b1(2),0,b1(3),0;0,c1(1),0,c1(2),0,c1(3);c1(1),b1(1),c1(2),b1(2),c1(3),b1(3)]/(2*A);
B2=[b2(1),0,b2(2),0,b2(3),0;0,c2(1),0,c2(2),0,c2(3);c2(1),b2(1),c2(2),b2(2),c2(3),b2(3)]/(2*A);
S1=D*B1;
S2=D*B2;
k=(B1')*D*B1*t*A;
% kkk=(B2')*D*B2*t*A; % k=kkk;

K=zeros(Nodes*2,Nodes*2);
qNodes=zeros(Nodes*2,1);
U=zeros(Nodes*2,1);
index1=[];index2=[];
%left
for i=2:Ny
    K(2*i-1:2*i,2*i-1:2*i)=k(1:2,1:2)+k(3:4,3:4)+k(5:6,5:6);
    K(2*i-1:2*i,2*(i+1)-1:2*(i+1))=k(3:4,5:6);
    K(2*i-1:2*i,2*(i+Ny+1)-1:2*(i+Ny+1))=k(3:4,1:2)+k(1:2,3:4);
    K(2*i-1:2*i,2*(i+Ny)-1:2*(i+Ny))=k(1:2,5:6)+k(5:6,1:2);
    K(2*i-1:2*i,2*(i-1)-1:2*(i-1))=k(5:6,3:4);
%     i
    U(2*i-1)=0;
    index1=[index1;2*i-1];
end
%bottom
for i=2:Nx
    s=i*(Ny+1);
    K(2*s-1:2*s,2*s-1:2*s)=k(1:2,1:2)+k(3:4,3:4)+k(5:6,5:6);
    K(2*s-1:2*s,2*(s-Ny-1)-1:2*(s-Ny-1))=k(3:4,1:2);
    K(2*s-1:2*s,2*(s+Ny+1)-1:2*(s+Ny+1))=k(1:2,3:4);
    K(2*s-1:2*s,2*(s+Ny)-1:2*(s+Ny))=k(5:6,1:2)+k(1:2,5:6);
    K(2*s-1:2*s,2*(s-1)-1:2*(s-1))=k(5:6,3:4)+k(3:4,5:6);
end
%right
for i=Nodes-Ny+1:Nodes-1
%     i
    K(2*i-1:2*i,2*i-1:2*i)=k(1:2,1:2)+k(3:4,3:4)+k(5:6,5:6);
    K(2*i-1:2*i,2*(i+1)-1:2*(i+1))=k(5:6,3:4);
    K(2*i-1:2*i,2*(i-Ny-1)-1:2*(i-Ny-1))=k(3:4,1:2)+k(1:2,3:4);
    K(2*i-1:2*i,2*(i-Ny)-1:2*(i-Ny))=k(1:2,5:6)+k(5:6,1:2);
    K(2*i-1:2*i,2*(i-1)-1:2*(i-1))=k(3:4,5:6);
end
%top
for i=1:Nx-1
    r=1+i*(Ny+1);
    K(2*r-1:2*r,2*r-1:2*r)=k(1:2,1:2)+k(3:4,3:4)+k(5:6,5:6);
    K(2*r-1:2*r,2*(r-Ny-1)-1:2*(r-Ny-1))=k(1:2,3:4);
    K(2*r-1:2*r,2*(r-Ny)-1:2*(r-Ny))=k(5:6,1:2)+k(1:2,5:6);
    K(2*r-1:2*r,2*(r+1)-1:2*(r+1))=k(5:6,3:4)+k(3:4,5:6);
    K(2*r-1:2*r,2*(r+Ny+1)-1:2*(r+Ny+1))=k(3:4,1:2);
    qNodes(2*r)=-qe;
end
%corner
K(1:2,1:2)=k(3:4,3:4);
K(1:2,3:4)=k(3:4,5:6);
K(1:2,2*(1+Ny+1)-1:2*(1+Ny+1))=k(3:4,1:2);
U(1)=0;
index1=[index1;1];
qNodes(2)=-qe/2;

K(2*(1+Ny)-1:2*(1+Ny),2*(1+Ny)-1:2*(1+Ny))=k(1:2,1:2)+k(5:6,5:6);
K(2*(1+Ny)-1:2*(1+Ny),2*Ny-1:2*(Ny))=k(5:6,3:4);
K(2*(1+Ny)-1:2*(1+Ny),2*(1+Ny+1+Ny)-1:2*(1+Ny+1+Ny))=k(1:2,3:4);
K(2*(1+Ny)-1:2*(1+Ny),2*(1+Ny+Ny)-1:2*(1+Ny+Ny))=k(1:2,5:6)+k(5:6,1:2);
U(2*(1+Ny)-1)=0;
index1=[index1;2*(1+Ny)-1];

K(2*Nodes-1:2*Nodes,2*Nodes-1:2*Nodes)=k(3:4,3:4);
K(2*Nodes-1:2*Nodes,2*(Nodes-1)-1:2*(Nodes-1))=k(3:4,5:6);
K(2*Nodes-1:2*Nodes,2*(Nodes-Ny-1)-1:2*(Nodes-Ny-1))=k(3:4,1:2);
U(2*Nodes)=0;
index2=[index2;2*Nodes];

K(2*(Nodes-Ny)-1:2*(Nodes-Ny),2*(Nodes-Ny)-1:2*(Nodes-Ny))=k(1:2,1:2)+k(5:6,5:6);
K(2*(Nodes-Ny)-1:2*(Nodes-Ny),2*(Nodes-Ny+1)-1:2*(Nodes-Ny+1))=k(5:6,3:4);
K(2*(Nodes-Ny)-1:2*(Nodes-Ny),2*(Nodes-Ny-Ny-1)-1:2*(Nodes-Ny-Ny-1))=k(1:2,3:4);
K(2*(Nodes-Ny)-1:2*(Nodes-Ny),2*(Nodes-Ny-Ny)-1:2*(Nodes-Ny-Ny))=k(5:6,1:2)+k(1:2,5:6);
qNodes(2*(Nodes-Ny))=-qe/2;

%inner
for i=1:Nx-1
    for j=1:Ny-1
        l=1+i*(Ny+1)+j;
        K(2*l-1:2*l,2*l-1:2*l)=2*k(5:6,5:6)+2*k(3:4,3:4)+2*k(1:2,1:2);
        K(2*l-1:2*l,2*(l-Ny-1)-1:2*(l-Ny-1))=k(3:4,1:2)+k(1:2,3:4);
        K(2*l-1:2*l,2*(l-Ny)-1:2*(l-Ny))=k(1:2,5:6)+k(5:6,1:2);
        K(2*l-1:2*l,2*(l+1)-1:2*(l+1))=k(5:6,3:4)+k(3:4,5:6);
        K(2*l-1:2*l,2*(l+Ny+1)-1:2*(l+Ny+1))=k(3:4,1:2)+k(1:2,3:4);
        K(2*l-1:2*l,2*(l+Ny)-1:2*(l+Ny))=k(1:2,5:6)+k(5:6,1:2);
        K(2*l-1:2*l,2*(l-1)-1:2*(l-1))=k(5:6,3:4)+k(3:4,5:6);
    end
end

KTemp=K;% save K as KTemp
qNodesTemp=qNodes;

index_total_zero=[index1;index2];
number_zero=numel(index_total_zero);
%删去位移为零对应的行和列
for i=1:number_zero
    maxindex=max(index_total_zero);
    K(maxindex,:)=[];
    K(:,maxindex)=[];
    qNodes(maxindex)=[];
    index_total_zero=setdiff(index_total_zero,maxindex);
end
%零位移标号
index_total_zero=[index1;index2];
%位移标号
index_total=[];
for i=1:Nodes*2
    index_total=[index_total;i];
end
%计算非零位移对应的在原位移序列中的标号
index_total_nonzero=setdiff(index_total,index_total_zero);

number_nonzero=numel(index_total_nonzero);
%计算非零位移
U_temp=inv(K)*qNodes;
%把计算所的的非零位移放回原位移序列中对应的位置
for i=1:number_nonzero
    U(index_total_nonzero(i))=U_temp(i);
end

Ux=zeros(Nodes,1);
Uy=zeros(Nodes,1);
for i=1:Nodes
    Ux(i)=U(2*i-1);
    Uy(i)=U(2*i);
end


xgrid=0:(Lx/Nx):Lx;
ygrid=Ly:-(Ly/Ny):0;
[X,Y]=meshgrid(xgrid,ygrid);

% % mesh
% M=(reshape(Ux,Ny+1,Nx+1));
% figure(1);
% subplot(2,2,1)
% Mesh(X,Y,M)
% xlabel('X-axis'),ylabel('Y-axis'),zlabel('x向位移')
% title('有限元法解简支梁问题计算结果数据可视化 图 1')
% 
% subplot(2,2,2)
% Mesh(X,Y,M)
% xlabel('X-axis'),ylabel('Y-axis'),zlabel('x向位移')
% title('图 2')
% view(30,20)
% 
% subplot(2,2,3)
% Mesh(X,Y,M)
% xlabel('X-axis'),ylabel('Y-axis'),zlabel('x向位移')
% title('图 3')
% view(-30,45)
% 
% subplot(2,2,4)
% surfl(X,Y,M)
% view(-60,45)
% shading interp
% xlabel('X-axis'),ylabel('Y-axis'),zlabel('x向位移')
% title('图 4')

%绘制网格
figure(2)
xmin1=0;%区域左端的x值
xmax1=Lx;%区域右端的x值
ymin1=0;%区域下端的y值
ymax1=Ly;%区域上端的y值
n1=Nx+1;%区域划分的纵线数
m1=Ny+1;%区域划分的横线数

a1=zeros(m1*n1);b1=zeros(m1*n1);
xy1=zeros(m1*2,2);x1=zeros(1,m1*n1);y1=zeros(1,m1*n1);
%计算节点坐标
for kk=1:n1
    for i=1:m1
        x1(i+(kk-1)*m1)=xmin1+(kk-1)*(xmax1-xmin1)/(n1-1);
        y1(i+(kk-1)*m1)=ymin1+(i-1)*(ymax1-ymin1)/(m1-1);
    end
end
clear kk i
a1(m1,m1)=2;
a1((n1-1)*m1+1,(n1-1)*m1+1)=2;
a1(1,1)=3;
a1(n1*m1,n1*m1)=3;
for i=1:(n1-2)
    a1(i*m1+1,i*m1+1)=4;
    a1((i+1)*m1,(i+1)*m1)=4;
end
clear i
for j=2:(m1-1)
    a1(j,j)=4;
    a1((n1-1)*m1+j,(n1-1)*m1+j)=4;
end
clear j
for i=1:(n1-2)
    for j=2:(m1-1)
        a1(i*m1+j,i*m1+j)=6;
    end
end
clear i j
a1=sparse(a1);
clear i j
for i=1:n1*m1-1
    b1(i,i+1)=-1;
end
clear i
for i=1:(n1-1)*m1
    b1(i,i+m1)=-1;
end
clear i
for i=1:(n1-1)*m1-1
    b1(i,i+m1+1)=-1;
end
for kk=1:n1-1
    for h=1:n1-1
        b1(m1*kk,m1*h+1)=0;
    end
end
clear kk h
b1=sparse(b1);
b1=a1+b1'+b1;
clear a1
a1=sparse(b1);
clear b1
xy1=[x1',y1'];
subplot(2,2,1);
gplot(a1,xy1);%绘制网格图


%compute the stress of every element F
UElement=zeros(6,Nx*Ny*2); % element's displacement
thigma=zeros(3,Nx*Ny*2);
for i=1:Nx
    for j=(i-1)*Ny+1:i*Ny
%         2*j-1;
%         j;
        w=j+i-1;
        UElement(:,2*j-1)=[Ux(w+Ny+1);Uy(w+Ny+1);Ux(w);Uy(w);Ux(w+1);Uy(w+1)];
        UElement(:,2*j)=[Ux(w+1);Uy(w+1);Ux(w+1+Ny+1);Uy(w+1+Ny+1);Ux(w+Ny+1);Uy(w+Ny+1)];
        thigma(:,2*j-1)=S1*UElement(:,2*j-1);
        thigma(:,2*j)=S2*UElement(:,2*j);
    end
%     disp('con');
end

%calculate the stress of every nodes
thigma_nodes=zeros(3,Nodes);
%modify the strss of nodes
%left
for i=2:Ny
    thigma_nodes(:,i)=(thigma(:,2*i-3)+thigma(:,2*i-2)+thigma(:,2*i-1))/3;
end

%bottom
for i=2:Nx
%     i
    s=i*(Ny+1);
%     2*(i-1)*Ny
%     2*i*Ny
%     2*i*Ny-1
%     disp('con');
    thigma_nodes(:,s)=(thigma(:,2*(i-1)*Ny)+thigma(:,2*i*Ny)+thigma(:,2*i*Ny-1))/3;
end

%right
for i=Nodes-Ny+1:Nodes-1
%     i
    ii=2*(Nx-1)*Ny+2+2*(i-(Nodes-Ny+1));
%     disp('con');
    thigma_nodes(:,i)=(thigma(:,ii)+thigma(:,ii+1)+thigma(:,ii+2))/3;
end

%top
for i=1:Nx-1
    r=1+i*(Ny+1);
%     2*(i-1)*Ny+1
%     2*(i-1)*Ny+2
%     2*i*Ny+1
%     disp('con');
    thigma_nodes(:,r)=(thigma(:,2*(i-1)*Ny+1)+thigma(:,2*(i-1)*Ny+2)+thigma(:,2*i*Ny+1))/3;
end

%corner
thigma_nodes(:,1)=thigma(:,1);
thigma_nodes(:,Nodes)=thigma(:,2*Nx*Ny);
thigma_nodes(:,1+Ny)=(thigma(:,2*Ny-1)+thigma(:,2*Ny))/2;
thigma_nodes(:,Nodes-Ny)=(thigma(:,2*(Nx-1)*Ny+1)+thigma(:,2*(Nx-1)*Ny+2))/2;


%inner
for i=1:Nx-1
    for j=1:Ny-1
        l=1+i*(Ny+1)+j;
        ll=2*(i-1)*Ny+2*j;
        lll=2*i*Ny+2*j;
%         disp('con');
    thigma_nodes(:,l)=(thigma(:,ll)+thigma(:,ll+1)+thigma(:,ll+2)+thigma(:,lll)+thigma(:,lll-1)+thigma(:,ll+1))/6;
    end
end

% mesh
M1=thigma_nodes(1,:);
M11=reshape(M1,Ny+1,Nx+1);

figure(1);
subplot(2,2,1)
Mesh(X,Y,M11)
xlabel('X-axis'),ylabel('Y-axis'),zlabel('Stress1')
title('有限元法解简支梁问题计算结果数据(Stress)可视化 图 1')

subplot(2,2,2)
Mesh(X,Y,M11)
xlabel('X-axis'),ylabel('Y-axis'),zlabel('Stress1')
title('图 2')
view(30,20)

subplot(2,2,3)
Mesh(X,Y,M11)
xlabel('X-axis'),ylabel('Y-axis'),zlabel('Stress1')
title('图 3')
view(-30,45)

subplot(2,2,4)
surfl(X,Y,M11)
view(-60,45)
shading interp
xlabel('X-axis'),ylabel('Y-axis'),zlabel('Stress1')
title('图 4')