76 lines
2.5 KiB
Matlab
76 lines
2.5 KiB
Matlab
%**************************************************************************
|
|
function [U, Strain, Stress, AxialForce]=TrussStaticsFEA(x,y,z,ele,Load,Constr,Scale)
|
|
Dofs=3*size(x,2); %总自由度数
|
|
EleCount=size(ele,1); %单元总数
|
|
K=zeros(Dofs,Dofs); %总体刚度矩阵
|
|
F=zeros(Dofs,1); %总体载荷列阵
|
|
U=zeros(Dofs,1); %总体位移列阵
|
|
BarLength=BarsLength(x,y,z,ele);
|
|
figure('Name','Undeformed Truss')
|
|
RenderTruss(ele,Load,Constr,x,y,z,U,1,1,'-k',1) %绘制桁架
|
|
axis equal
|
|
figure('Name','Deformed Truss')
|
|
RenderTruss(ele,Load,Constr,x,y,z,U,0.5,0.5,'-k',1) %绘制变形后桁架
|
|
axis equal
|
|
hold on
|
|
%遍历所有单元，将各单元刚度阵分块组装到总体刚度阵
|
|
for iEle =1:EleCount
|
|
%该单元的两个节点的编号
|
|
n1=ele(iEle,2);n2=ele(iEle,3);
|
|
%计算坐标变换矩阵
|
|
R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],[z(n1) z(n2)],BarLength(iEle));
|
|
%计算单元刚度矩阵
|
|
ke= BarElementKe(ele(iEle,4),ele(iEle,5),R,BarLength(iEle));
|
|
%将各单元刚度分块组装到总刚相应位置
|
|
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
|
|
%形成载荷列阵--1、2、3、4、5、6自由度赋载荷值
|
|
for LoadNum=1:size(Load,1)
|
|
for i=2:4
|
|
F(3*Load(LoadNum)+i-4,1)=Load(LoadNum,i);
|
|
end
|
|
end
|
|
%施加约束--乘大数法
|
|
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; %全局坐标系下位移
|
|
for iEle =1:EleCount
|
|
%计算杆局部坐标下的位移
|
|
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; %应变
|
|
Stress(1, iEle)=ele(iEle,5)* Strain(1, iEle); %应力
|
|
AxialForce(1, iEle)=ele(iEle,4)* Stress(1, iEle); %轴力
|
|
end
|
|
%保存位移，应力与轴力到文本文件
|
|
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) %绘制变形后桁架
|
|
figure
|
|
AxisForcePlot(ele,Load,Constr,x,y,z,U,1,2.5,'-b',Scale,AxialForce)
|
|
axis equal
|
|
end %主程序结束
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|