Structural-Mechanical
/
FEM-Course-Matlab
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
FEM-Course-Matlab/14.弹性地基梁的matlab有限元编程/main_frame.m

148 lines
4.3 KiB
Mathematica
Raw Permalink Normal View History

init project
2024-01-28 16:46:36 +00:00
clc,clear,close all
global m b D
height=20;%׮<EFBFBD>ܸ<EFBFBD>
deep=9;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
width=2;%׮<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
D=0.8;%ֱ<EFBFBD><EFBFBD>
E1=3e10;%׮<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD>
E2=3e10;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD>
b=0.8;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
h=0.6;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
m=4e6;%<EFBFBD><EFBFBD><EFBFBD>ɷ<EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
%%%%%%%%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>q1 q2%%%%%%%%%%%%
phi=25;%<EFBFBD><EFBFBD>Ħ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
q=10;%qΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
y=19.2;%<EFBFBD><EFBFBD><EFBFBD>ض<EFBFBD>
c=12;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ճ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
ka=tand(45-phi/2)^2;
sigma_a=(q+y*deep)*ka-2*c*sqrt(ka);
L0=deep*tand(45-phi/2);
beta=2*width/L0-(width/L0)^2;
Pa1=beta*sigma_a;%KPa
Pa2=(1-beta)*sigma_a;%KPa
p1=-Pa1*1e3;%Pa
p2=-Pa2*1e3;%Pa
%define nodes
for i=1:height+1
nodeCoord(i,:)=[0,i-1];%1mһ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ
end
nodeCoord(height+2,:)=[width,height];
for i=1:height
nodeCoord(height+2+i,:)=[width,height-i];%1mһ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ
end
node_num=length(nodeCoord(:,1));
x=nodeCoord(:,1)'; % <EFBFBD>ڵ<EFBFBD>x<EFBFBD>᷽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
y=nodeCoord(:,2)'; % <EFBFBD>ڵ<EFBFBD>y<EFBFBD>᷽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
A(1:height)=pi*D^2/4;
A(height+1)=b*h;
A(height+2:2*height+1)=pi*D^2/4;
I(1:height)=pi*D^4/36;
I(height+1)=b*h^3/12;
I(height+2:2*height+1)=pi*D^4/36;
E(1:height)=E1;
E(height+1)=E2;
E(height+2:2*height+1)=E1;
%<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>Ϣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ţ<EFBFBD><EFBFBD>ڵ<EFBFBD>1<EFBFBD><EFBFBD><EFBFBD>ţ<EFBFBD><EFBFBD>ڵ<EFBFBD>2<EFBFBD><EFBFBD><EFBFBD>ţ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD>
%define ele
for i=1:node_num-1
ele(i,:)=[i,i,i+1,A(i),E(i),I(i)];
end
q1=D*p1;%׮1<EFBFBD><EFBFBD><EFBFBD>ߺ<EFBFBD><EFBFBD><EFBFBD>
q2=D*p2;%׮2<EFBFBD><EFBFBD><EFBFBD>ߺ<EFBFBD><EFBFBD><EFBFBD>
%<EFBFBD>غ<EFBFBD><EFBFBD><EFBFBD>Ϣ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ţ<EFBFBD>x<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>y<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>z<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
for i=1:node_num
if y(i)<=deep && x(i)==min(x)
Load(i,:)=[i q1 0 0];
elseif y(i)>deep && x(i)==min(x)
Load(i,:)=[i q1/(height-deep)*(height-y(i)) 0 0];
elseif y(i)<=deep && x(i)==max(x)
Load(i,:)=[i q2 0 0];
elseif y(i)>deep && x(i)==max(x)
Load(i,:)=[i q2/(height-deep)*(height-y(i)) 0 0];
end
end
%Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ţ<EFBFBD>x<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>y<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>z<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD>
Constr=[1 NaN 0 NaN;node_num NaN 0 NaN];
Dofs=3*size(x,2); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD><EFBFBD><EFBFBD>
EleCount=size(ele,1); %<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
K=zeros(Dofs,Dofs); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
F=zeros(Dofs,1); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>غ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
U=zeros(Dofs,1); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
BarLength=BarsLength(x,y,ele);
figure('Name','Undeformed Truss')
RenderFrame(ele,Load,Constr,x,y,U,1,'-k',1) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
figure('Name','Deformed Truss')
RenderFrame(ele,Load,Constr,x,y,U,0.5,'-k',1) %<EFBFBD><EFBFBD><EFBFBD>Ʊ<EFBFBD><EFBFBD>κ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
hold on
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>е<EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD>ն<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֿ<EFBFBD><EFBFBD><EFBFBD>װ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ն<EFBFBD><EFBFBD><EFBFBD>
for iEle =1:EleCount
%<EFBFBD>õ<EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD>ı<EFBFBD><EFBFBD><EFBFBD>
n1=ele(iEle,2);n2=ele(iEle,3);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>任<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],BarLength(iEle));
%<EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD> Ke=R'*ke*R;<EFBFBD>ֲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD>µĵ<EFBFBD>Ԫ<EFBFBD>ն<EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD>Ϊȫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>µĵ<EFBFBD>Ԫ<EFBFBD>ն<EFBFBD><EFBFBD><EFBFBD>
if y(n1)<=deep
ke= FrameElementKe2(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
else
ke= FrameElementKe1(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
end
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD>նȷֿ<EFBFBD><EFBFBD><EFBFBD>װ<EFBFBD><EFBFBD><EFBFBD>ܸ<EFBFBD><EFBFBD><EFBFBD>Ӧλ<EFBFBD><EFBFBD>
eleDof=[n1*3-2:n1*3,n2*3-2:n2*3];
K(eleDof,eleDof)=K(eleDof,eleDof)+ke;
end
%<EFBFBD>γ<EFBFBD><EFBFBD>غ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>--1<EFBFBD><EFBFBD>2<EFBFBD><EFBFBD>3<EFBFBD><EFBFBD>4<EFBFBD><EFBFBD><EFBFBD>ɶȸ<EFBFBD><EFBFBD>غ<EFBFBD>ֵ
for LoadNum=1:size(Load,1)
for i=1:3
F(3*Load(LoadNum)+i-3,1)=Load(LoadNum,i+1);
end
end
%ʩ<EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD>--<EFBFBD>˴<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
for iConstr=1:size(Constr,1)
for j=2:4
if ~isnan(Constr(iConstr,j))
K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1 )+j-4)=...
1e12*K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1 )+j-4);
F(3*Constr(iConstr,1)+j-4)=1e12*Constr(iConstr,j)*...
K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1)+j-4);
end
end
end
U=K\F; %ȫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>,<EFBFBD>˴<EFBFBD>ԭ<EFBFBD><EFBFBD>Ϊ<EFBFBD>浫<EFBFBD>Ǿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>죬<EFBFBD>ʸij<EFBFBD>α<EFBFBD><EFBFBD>
for iEle =1:EleCount
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>˾ֲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>µ<EFBFBD>λ<EFBFBD><EFBFBD>
n1=ele(iEle,2);n2=ele(iEle,3);
R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],BarLength(iEle));
localU = R*[U(3*n1-2:3*n1,1);U(3*n2-2:3*n2,1)];
if y(n1)<=deep
[Ke_local] = FrameElementKeLocal1(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
else
[Ke_local] = FrameElementKeLocal2(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
end
EleForce(:, iEle)=Ke_local*localU; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
end
MomentForce=EleForce([3,6],:);
figure
MomentForcePlot1=[-MomentForce(1,1:20)];
plot(MomentForcePlot1,[-19:0])
hold on;
MomentForcePlot2=[MomentForce(2,22:41)];
plot(-MomentForcePlot2,[0:-1:-19])
legend('ǰ','<EFBFBD><EFBFBD>')
title('Moment');
figure
for i=1:node_num
Ux(i)=U((i-1)*3+1);
end
plot(Ux(1:21),[-20:0])
hold on;
plot(Ux(22:42),[0:-1:-20])
title('Displacement');
legend('ǰ','<EFBFBD><EFBFBD>')
figure
MomentForce(1,:)=-MomentForce(1,:);
ForcePlot(ele,Load,Constr,x,y,U,3,'-b',MomentForce)%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ
title('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ');