FEM-Course-Matlab/6.桁架结构matlab有限元编程/Truss3D/Truss3DStaticsFEA.m

76 lines
2.5 KiB
Mathematica
Raw Permalink Normal View History

2024-01-28 16:46:36 +00:00
%**************************************************************************
function [U, Strain, Stress, AxialForce]=TrussStaticsFEA(x,y,z,ele,Load,Constr,Scale)
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,z,ele);
figure('Name','Undeformed Truss')
RenderTruss(ele,Load,Constr,x,y,z,U,1,1,'-k',1) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
axis equal
figure('Name','Deformed Truss')
RenderTruss(ele,Load,Constr,x,y,z,U,0.5,0.5,'-k',1) %<EFBFBD><EFBFBD><EFBFBD>Ʊ<EFBFBD><EFBFBD>κ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
axis equal
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)],[z(n1) z(n2)],BarLength(iEle));
%<EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
ke= BarElementKe(ele(iEle,4),ele(iEle,5),R,BarLength(iEle));
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD>նȷֿ<EFBFBD><EFBFBD><EFBFBD>װ<EFBFBD><EFBFBD><EFBFBD>ܸ<EFBFBD><EFBFBD><EFBFBD>Ӧλ<EFBFBD><EFBFBD>
K(3*n1-2:3*n1,3*n1-2:3*n1)=K(3*n1-2:3*n1,3*n1-2:3*n1)+ke(1:3,1:3);
K(3*n1-2:3*n1,3*n2-2:3*n2)=K(3*n1-2:3*n1,3*n2-2:3*n2)+ke(1:3,4:6);
K(3*n2-2:3*n2,3*n1-2:3*n1)=K(3*n2-2:3*n2,3*n1-2:3*n1)+ke(4:6,1:3);
K(3*n2-2:3*n2,3*n2-2:3*n2)=K(3*n2-2:3*n2,3*n2-2:3*n2)+ke(4:6,4:6);
end
%<EFBFBD>γ<EFBFBD><EFBFBD>غ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>--1<EFBFBD><EFBFBD>2<EFBFBD><EFBFBD>3<EFBFBD><EFBFBD>4<EFBFBD><EFBFBD>5<EFBFBD><EFBFBD>6<EFBFBD><EFBFBD><EFBFBD>ɶȸ<EFBFBD><EFBFBD>غ<EFBFBD>ֵ
for LoadNum=1:size(Load,1)
for i=2:4
F(3*Load(LoadNum)+i-4,1)=Load(LoadNum,i);
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>
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)],[z(n1) z(n2)],BarLength(iEle));
localU = R*[U(3*n1-2:3*n1,1);U(3*n2-2:3*n2,1)];
Strain(1, iEle)=[-1/BarLength(iEle) 1/BarLength(iEle)]*localU; %Ӧ<EFBFBD><EFBFBD>
Stress(1, iEle)=ele(iEle,5)* Strain(1, iEle); %Ӧ<EFBFBD><EFBFBD>
AxialForce(1, iEle)=ele(iEle,4)* Stress(1, iEle); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
end
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD>ƣ<EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ı<EFBFBD><EFBFBD>ļ<EFBFBD>
fp=fopen('Result.txt','a');
str = [char(13,10)','U',' ',num2str(U'),char(13,10)','Stress',' ',...
num2str(Stress),char(13,10)','AxialForce',' ',num2str(AxialForce)];
fprintf(fp,str);
fclose(fp);
RenderTruss(ele,Load,Constr,x,y,z,U,1,1,'-.b',Scale) %<EFBFBD><EFBFBD><EFBFBD>Ʊ<EFBFBD><EFBFBD>κ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
figure
AxisForcePlot(ele,Load,Constr,x,y,z,U,1,2.5,'-b',Scale,AxialForce)
axis equal
end %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>