📄 104776-beam3d.m
字号:
function beam3d
% 空间弹性梁单元结构－静力分析
global node element material K knode
file_in=input('input file name: ','s');
file_out=input('output file name: ','s');
fid_in=fopen(file_in,'r');
node_number=fscanf(fid_in,'%d',1);
node=zeros(node_number,3);
for i=1:node_number
nn=fscanf(fid_in,'%d',1);
node(i,:)=fscanf(fid_in,'%f',[1,3]);
end
knode_number=fscanf(fid_in,'%d',1);
knode=zeros(knode_number,3);
for i=1:knode_number
nn=fscanf(fid_in,'%d',1);
knode(i,:)=fscanf(fid_in,'%f',[1,3]);
end
element_number=fscanf(fid_in,'%d',1);
element=zeros(element_number,4);
for i=1:element_number
ne=fscanf(fid_in,'%d',1);
element(i,:)=fscanf(fid_in,'%d',[1,4]);
end
material_number=fscanf(fid_in,'%d',1);
material=zeros(material_number,7);
for i=1:material_number
nm=fscanf(fid_in,'%d',1);
material(i,:)=fscanf(fid_in,'%f',[1,7]);
end
bc_number=fscanf(fid_in,'%d',1);
bc=zeros(bc_number,3);
for i=1:bc_number
nb=fscanf(fid_in,'%d',1);
bc(i,1)=fscanf(fid_in,'%d',1);
bc(i,2)=fscanf(fid_in,'%d',1);
bc(i,3)=fscanf(fid_in,'%f',1);
end
nf_number=fscanf(fid_in,'%d',1);
nf=zeros(nf_number,3);
for i=1:nf_number
nnf=fscanf(fid_in,'%d',1);
nf(i,1)=fscanf(fid_in,'%d',1);
nf(i,2)=fscanf(fid_in,'%d',1);
nf(i,3)=fscanf(fid_in,'%f',1);
end
fclose(fid_in);
K=sparse(node_number*6,node_number*6);
f=sparse(node_number*6,1);
for ie=1:element_number
k=StiffnessMatrix(ie);
AssembleStiffnessMatrix(ie,k);
end
for ibc=1:bc_number
n=bc(ibc,1);
d=bc(ibc,2);
m=(n-1)*6+d;
f(m)=bc(ibc,3)*K(m,m)*1E15;
K(m,m)=K(m,m)*1E15;
end
for inf=1:nf_number
n=nf(inf,1);
d=nf(inf,2);
f((n-1)*6+d)=nf(inf,3);
end
delta=K\f;
%disp(delta);
fid_out=fopen(file_out,'w');
for i=1:node_number
fprintf(fid_out,'node=%d -- x= %.5f y= %.5f z= %.5f ax= %.5f ay= %.5f az= %.5f\r\n',i,delta((i-1)*6+1),delta((i-1)*6+2),delta((i-1)*6+3),delta((i-1)*6+4),delta((i-1)*6+5),delta((i-1)*6+6));
end
for i=1:element_number
fprintf(fid_out,'\r\nElement %d\r\n\r\n',i);
for j=1:2
fprintf(fid_out,' Node %d\r\n',j);
for p=1:6
fprintf(fid_out,' F%d= %.5f\r\n',p,getforce(delta,i,j,p));
end
end
end
fclose(fid_out);
return;
function k=StiffnessMatrix(ie)
global node element material K knode
k=zeros(12,12);
E=material(element(ie,3),1);
G=material(element(ie,3),2);
Iy=material(element(ie,3),3);
Iz=material(element(ie,3),4);
Jx=material(element(ie,3),5);
A=material(element(ie,3),6);
xi=node(element(ie,1),1);
yi=node(element(ie,1),2);
zi=node(element(ie,1),3);
xj=node(element(ie,2),1);
yj=node(element(ie,2),2);
zj=node(element(ie,2),3);
xk=knode(element(ie,4),1);
yk=knode(element(ie,4),2);
zk=knode(element(ie,4),3);
L=((xj-xi)^2+(yj-yi)^2+(zj-zi)^2)^(1/2);
Lk=((xk-xi)^2+(yk-yi)^2+(zk-zi)^2)^(1/2);
l1=(xj-xi)/L;
m1=(yj-yi)/L;
n1=(zj-zi)/L;
g1=(xk-xi)/Lk;
g2=(yk-yi)/Lk;
g3=(zk-zi)/Lk;
s=((m1*g3-n1*g2)^2+(n1*g1-l1*g3)^2+(l1*g2-m1*g1)^2)^(1/2);
k=[E*A/L 0 0 0 0 0 -E*A/L 0 0 0 0 0
0 12*E*Iz/L^3 0 0 0 6*E*Iz/L^2 0 -12*E*Iz/L^3 0 0 0 6*E*Iz/L^2
0 0 12*E*Iy/L^3 0 -6*E*Iy/L^2 0 0 0 -12*E*Iy/L^3 0 -6*E*Iy/L^2 0
0 0 0 G*Jx/L 0 0 0 0 0 -G*Jx/L 0 0
0 0 -6*E*Iy/L^2 0 4*E*Iy/L 0 0 0 6*E*Iy/L^2 0 2*E*Iy/L 0
0 6*E*Iz/L^2 0 0 0 4*E*Iz/L 0 -6*E*Iz/L^2 0 0 0 2*E*Iz/L
-E*A/L 0 0 0 0 0 E*A/L 0 0 0 0 0
0 -12*E*Iz/L^3 0 0 0 -6*E*Iz/L^2 0 12*E*Iz/L^3 0 0 0 -6*E*Iz/L^2
0 0 -12*E*Iy/L^3 0 6*E*Iy/L^2 0 0 0 12*E*Iy/L^3 0 6*E*Iy/L^2 0
0 0 0 -G*Jx/L 0 0 0 0 0 G*Jx/L 0 0
0 0 -6*E*Iy/L^2 0 2*E*Iy/L 0 0 0 6*E*Iy/L^2 0 4*E*Iy/L 0
0 6*E*Iz/L^2 0 0 0 2*E*Iz/L 0 -6*E*Iz/L^2 0 0 0 4*E*Iz/L];
% k点取自X-Y平面内
t11=l1;
t12=(g1-l1*(l1*g1+m1*g2+n1*g3))/s;
t13=(m1*g3-n1*g2)/s;
t21=m1;
t22=(g2-m1*(l1*g1+m1*g2+n1*g3))/s;
t23=(n1*g1-l1*g3)/s;
t31=n1;
t32=(g3-n1*(l1*g1+m1*g2+n1*g3))/s;
t33=(l1*g2-m1*g1)/s;
t=[t11 t12 t13 0 0 0 0 0 0 0 0 0
t21 t22 t23 0 0 0 0 0 0 0 0 0
t31 t32 t33 0 0 0 0 0 0 0 0 0
0 0 0 t11 t12 t13 0 0 0 0 0 0
0 0 0 t21 t22 t23 0 0 0 0 0 0
0 0 0 t31 t32 t33 0 0 0 0 0 0
0 0 0 0 0 0 t11 t12 t13 0 0 0
0 0 0 0 0 0 t21 t22 t23 0 0 0
0 0 0 0 0 0 t31 t32 t33 0 0 0
0 0 0 0 0 0 0 0 0 t11 t12 t13
0 0 0 0 0 0 0 0 0 t21 t22 t23
0 0 0 0 0 0 0 0 0 t31 t32 t33];
k=t*k*transpose(t);
return
function AssembleStiffnessMatrix(ie,k)
global element K
for i=1:2
for j=1:2
for p=1:6
for q=1:6
m=(i-1)*6+p;
n=(j-1)*6+q;
M=(element(ie,i)-1)*6+p;
N=(element(ie,j)-1)*6+q;
K(M,N)=K(M,N)+k(m,n);
end
end
end
end
return
function shearf=getforce(b,ie,i,d) %获得杆端内力 位移 单元 杆端 自由度
global node element material K knode
k=zeros(12,12);
E=material(element(ie,3),1);
G=material(element(ie,3),2);
Iy=material(element(ie,3),3);
Iz=material(element(ie,3),4);
Jx=material(element(ie,3),5);
A=material(element(ie,3),6);
xi=node(element(ie,1),1);
yi=node(element(ie,1),2);
zi=node(element(ie,1),3);
xj=node(element(ie,2),1);
yj=node(element(ie,2),2);
zj=node(element(ie,2),3);
xk=knode(element(ie,4),1);
yk=knode(element(ie,4),2);
zk=knode(element(ie,4),3);
L=((xj-xi)^2+(yj-yi)^2+(zj-zi)^2)^(1/2);
Lk=((xk-xi)^2+(yk-yi)^2+(zk-zi)^2)^(1/2);
l1=(xj-xi)/L;
m1=(yj-yi)/L;
n1=(zj-zi)/L;
g1=(xk-xi)/Lk;
g2=(yk-yi)/Lk;
g3=(zk-zi)/Lk;
s=((m1*g3-n1*g2)^2+(n1*g1-l1*g3)^2+(l1*g2-m1*g1)^2)^(1/2);
k=[E*A/L 0 0 0 0 0 -E*A/L 0 0 0 0 0
0 12*E*Iz/L^3 0 0 0 6*E*Iz/L^2 0 -12*E*Iz/L^3 0 0 0 6*E*Iz/L^2
0 0 12*E*Iy/L^3 0 -6*E*Iy/L^2 0 0 0 -12*E*Iy/L^3 0 -6*E*Iy/L^2 0
0 0 0 G*Jx/L 0 0 0 0 0 -G*Jx/L 0 0
0 0 -6*E*Iy/L^2 0 4*E*Iy/L 0 0 0 6*E*Iy/L^2 0 2*E*Iy/L 0
0 6*E*Iz/L^2 0 0 0 4*E*Iz/L 0 -6*E*Iz/L^2 0 0 0 2*E*Iz/L
-E*A/L 0 0 0 0 0 E*A/L 0 0 0 0 0
0 -12*E*Iz/L^3 0 0 0 -6*E*Iz/L^2 0 12*E*Iz/L^3 0 0 0 -6*E*Iz/L^2
0 0 -12*E*Iy/L^3 0 6*E*Iy/L^2 0 0 0 12*E*Iy/L^3 0 6*E*Iy/L^2 0
0 0 0 -G*Jx/L 0 0 0 0 0 G*Jx/L 0 0
0 0 -6*E*Iy/L^2 0 2*E*Iy/L 0 0 0 6*E*Iy/L^2 0 4*E*Iy/L 0
0 6*E*Iz/L^2 0 0 0 2*E*Iz/L 0 -6*E*Iz/L^2 0 0 0 4*E*Iz/L];
% k点取自X-Y平面内
t11=l1;
t12=(g1-l1*(l1*g1+m1*g2+n1*g3))/s;
t13=(m1*g3-n1*g2)/s;
t21=m1;
t22=(g2-m1*(l1*g1+m1*g2+n1*g3))/s;
t23=(n1*g1-l1*g3)/s;
t31=n1;
t32=(g3-n1*(l1*g1+m1*g2+n1*g3))/s;
t33=(l1*g2-m1*g1)/s;
t=[t11 t12 t13 0 0 0 0 0 0 0 0 0
t21 t22 t23 0 0 0 0 0 0 0 0 0
t31 t32 t33 0 0 0 0 0 0 0 0 0
0 0 0 t11 t12 t13 0 0 0 0 0 0
0 0 0 t21 t22 t23 0 0 0 0 0 0
0 0 0 t31 t32 t33 0 0 0 0 0 0
0 0 0 0 0 0 t11 t12 t13 0 0 0
0 0 0 0 0 0 t21 t22 t23 0 0 0
0 0 0 0 0 0 t31 t32 t33 0 0 0
0 0 0 0 0 0 0 0 0 t11 t12 t13
0 0 0 0 0 0 0 0 0 t21 t22 t23
0 0 0 0 0 0 0 0 0 t31 t32 t33];
globdelta=zeros(12,1);
for ij=1:2
nodenum=element(ie,ij);
for p=1:6
globdelta((ij-1)*6+p,1)=b((nodenum-1)*6+p,1);
end
end
delta=transpose(t)*globdelta;
f=k*delta;
shearf=f((i-1)*6+d);
return
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -