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

108 lines
3.8 KiB
Matlab
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

clc,clear,close all
tic
nodes=zeros(3,244); %创建节点
ele0=zeros(2,5);ele1=zeros(24,5);ele2=zeros(60,5);ele3=zeros(20,5);
ele4=zeros(64,5);ele5=zeros(168,5);ele6=zeros(28,5);ele7=zeros(144,5);
ele8=zeros(396,5);ele9=zeros(36,5);
d1=4270/2;d2=5335;r2=4265;d3=7470;r3=6400;d4=9605;r4=8535;
l3=43890;l2=l3+29260;l1=l2+21950;
dl3=43890/12;dl2=29260/8;dl1=21950/6;
coordinate2=[d2,r2,0,-r2,-d2,-r2,0,r2;0,r2,d2,r2,0,-r2,-d2,-r2];%定义坐标
coordinate3=[d3,r3,d1,-d1,-r3,-d3,-d3,-r3,-d1,d1,r3,d3;...
d1,r3,d3,d3,r3,d1,-d1,-r3,-d3,-d3,-r3,-d1];
coordinate4=[d4,r4,d1,-d1,-r4,-d4,-d4,-r4,-d1,d1,r4,d4;...
d1,r4,d4,d4,r4,d1,-d1,-r4,-d4,-d4,-r4,-d1];
%节点坐标
for i=1:6
nodes(:,4*i-3:4*i)=[d1,-d1,-d1,d1;d1,d1,-d1,-d1;l1-(i-1)*dl1*ones(1,4)];
end
for i=1:8
nodes(:,24+8*(i-1)+1:24+8*i)=[coordinate2;l2-(i-1)*dl2*ones(1,8)];
end
for i=1:12
nodes(:,88+12*(i-1)+1:88+12*i)=[coordinate3;l3-(i-1)*dl3*ones(1,12)];
end
nodes(:,233:244)=[coordinate4;zeros(1,12)];
%单元信息:单元编号,单元节点,单元节点
ele0(:,2:3)=[1,3;2,4];
for i=1:6
ele1(4*i-3:4*i,2:3)=[4*i-3,4*i-2;4*i-2,4*i-1;4*i-1,4*i;4*i,4*i-3];%1-24
if i==6
ele3(1:5,2:3)=[21,32;21,25;21,26;21,27;21,28];
ele3(6:10,2:3)=[22,26;22,27;22,28;22,29;22,30];
ele3(11:15,2:3)=[23,28;23,29;23,30;23,31;23,32];
ele3(16:20,2:3)=[24,30;24,31;24,32;24,25;24,26];
else
ele2(12*(i-1)+3-2:12*(i-1)+3,2:3)=...
[4*(i-1)+1,4*i+4; 4*(i-1)+1,4*i+1; 4*(i-1)+1,4*i+2];
for j=2:3
ele2(12*(i-1)+3*j-2:12*(i-1)+3*j,2:3)=...
[4*(i-1)+j,4*i+j-1; 4*(i-1)+j,4*i+j; 4*(i-1)+j,4*i+j+1];
end
ele2(12*(i-1)+12-2:12*(i-1)+12,2:3)=[4*(i-1)+4,4*i+4-1;...
4*(i-1)+4,4*i+4; 4*(i-1)+4,4*i+1];
end
end
for i=1:8
ele4(8*i-7:8*i,2:3)=[8*i-7,8*i-6;8*i-6,8*i-5;8*i-5,8*i-4;...
8*i-4,8*i-3;8*i-3,8*i-2;8*i-2,8*i-1;8*i-1,8*i;8*i,8*i-7;];
if i==8
ele6(1:4,2:3)=[81,99;81,100;81,89;81,90];
for j=1:3
ele6(4*(j-1)+5:4*(j-1)+8,2:3)=[81+2*j,88+3*j-1;...
81+2*j,88+3*j;81+2*j,88+3*j+1;81+2*j,88+3*j+2];
end
for j=1:4
ele6(3*(j-1)+17:3*(j-1)+19,2:3)=[80+2*j,89+3*(j-1);...
80+2*j,89+3*(j-1)+1;80+2*j,89+3*(j-1)+2];
end
else
ele5(24*(i-1)+1:24*(i-1)+3,2:3)=...
[8*(i-1)+1,8*i+8;8*(i-1)+1,8*i+1; 8*(i-1)+1,8*i+2];
for j=2:7
ele5(24*(i-1)+3*j-2:24*(i-1)+3*j,2:3)=...
[8*(i-1)+j,8*i+j-1; 8*(i-1)+j,8*i+j; 8*(i-1)+j,8*i+j+1];
end
ele5(24*(i-1)+24-2:24*(i-1)+24,2:3)=...
[8*(i-1)+8,8*i+8-1; 8*(i-1)+8,8*i+8; 8*(i-1)+8,8*i+1];
end
end
ele4=ele4+24;
ele5=ele5+24;
for i=1:12
ele7(12*i-11:12*i,2:3)=[12*i-11,12*i-10;12*i-10,12*i-9;12*i-9,...
12*i-8;12*i-8,12*i-7;12*i-7,12*i-6;12*i-6,12*i-5;12*i-5,12*i-4;...
12*i-4,12*i-3;12*i-3,12*i-2;12*i-2,12*i-1;12*i-1,12*i;12*i,12*i-11;];%1-24
if i==12
ele9(1:3,2:3)=[233,232;233,221;233,222];
for j=2:11
ele9(3*(j-1)+1:3*(j-1)+3,2:3)=[232+j,220+j-1;232+j,220+j;232+j,220+j+1];
end
ele9(34:36,2:3)=[244,231;244,232;244,221];
else
ele8(36*(i-1)+1:36*(i-1)+3,2:3)=[12*(i-1)+1,...
12*i+12;12*(i-1)+1,12*i+1; 12*(i-1)+1,12*i+2];
for j=2:11
ele8(36*(i-1)+3*(j-1)+1:36*(i-1)+3*j,2:3)=[12*(i-1)+j,...
12*i+j-1; 12*(i-1)+j,12*i+j; 12*(i-1)+j,12*i+j+1];
end
ele8(36*(i-1)+36-2:36*(i-1)+36,2:3)=[12*(i-1)+12,12*i+12-1; ...
12*(i-1)+12,12*i+12; 12*(i-1)+12,12*i+1];
end
end
ele7=ele7+88;
ele8=ele8+88;
ele=[ele0;ele1;ele2;ele3;ele4;ele5;ele6;ele7;ele8;ele9];
for i=1:length(ele(:,1))
ele(i,1)=i;
end
A=4;E=2.1E005; % 定义各单元截面面积和材料弹性模量
ele(:,4:5)=ones(length(ele(:,1)),2)*[A 0;0 E]; %定义单元信息: 截面积,弹性模量
Load=[1 0 400 -100;2 0 400 -100]; %载荷信息单元编号x方向力y方向力z方向力
x=nodes(1,:);y=nodes(2,:);z=nodes(3,:);
Constr=zeros(12,4);%定义约束
for i=1:12
Constr(i,1)=i+232;
end
[U, Strain, Stress, AxialForce]=Truss3DStaticsFEA(x,y,z,ele,Load,Constr,10)
toc