📄 104777-beam3dm.m
字号:
function beam3dm
% 空间弹性梁单元结构－模态分析
global node element material K knode M
file_in=input('input file name: ','s');
%file_out=input('output file name: ','s');
tm=input('Modal type No.: ','s');
typem=str2num(tm);
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
fclose(fid_in);
K=sparse(node_number*6,node_number*6);
M=sparse(node_number*6,node_number*6);
for ie=1:element_number
k=StiffnessMatrix(ie);
AssembleStiffnessMatrix(ie,k);
me=MassMatrix(ie);
AssembleMassMatrix(ie,me);
end
for ibc=1:bc_number
n=bc(ibc,1);
d=bc(ibc,2);
m=(n-1)*6+d;
K(m,m)=K(m,m)*1E15;
M(m,m)=M(m,m)*1E15;
end
m=typem;
% H=(M^-1)*K;
% [u1,ev]=eigs(H,m,'SM');
% 利用子空间迭代法求H的m个最小特征值
n=node_number*6;
u1=2*rand(n,m)-1; % u0为迭代初始值
er=zeros(m,m);
H2=(K^-1)*M;
while(1)
u1=H2*u1;
for i=1:m
s=max(abs(u1(:,i)));
u1(:,i)=u1(:,i)/s;
end
Km=u1'*K*u1;
Mm=u1'*M*u1;
%Hm=(Mm^-1)*Km;
[A,ev]=eig(Km,Mm);
ev
e=ev-er;
for i=1:m
s=max(abs(A(:,i)));
A(:,i)=A(:,i)/s;
end
u=u1*A;
% eu=u-u1;
if sum(diag(abs(e)))<0.000000001
break;
else
er=ev;
u1=u;
end
end
pi=3.1415926535897932384626;
f=1000*sqrt(diag(ev))/pi/2
% u1=A;
% figure(1);
% for i=1:element_number
% xe(1,1)=node(element(i,1),1);
% xe(2,1)=node(element(i,2),1);
% ye(1,1)=node(element(i,1),2);
% ye(2,1)=node(element(i,2),2);
% ze(1,1)=node(element(i,1),3);
% ze(2,1)=node(element(i,2),3);
% line(xe,ye,ze);
% end
nodem=node;
mj=typem;
for i=1:node_number
nodem(i,1)=nodem(i,1)+u1((i-1)*6+1,mj)*20;
nodem(i,2)=nodem(i,2)+u1((i-1)*6+2,mj)*20;
nodem(i,3)=nodem(i,3)+u1((i-1)*6+3,mj)*20;
end
for i=1:element_number
xe(1,1)=nodem(element(i,1),1);
xe(2,1)=nodem(element(i,2),1);
ye(1,1)=nodem(element(i,1),2);
ye(2,1)=nodem(element(i,2),2);
ze(1,1)=nodem(element(i,1),3);
ze(2,1)=nodem(element(i,2),3);
line(xe,ye,ze);
end
return;
function k=StiffnessMatrix(ie)
global node element material K knode
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 me=MassMatrix(ie)
global node element material K knode M
%me=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);
m=material(element(ie,3),7); %线重度
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);
% 计入转动惯量影响的一致质量矩阵
me=[ 1/3*m*L, 0, 0, 0, 0, 0, 1/6*m*L, 0, 0, 0, 0, 0
0, 13/35*m*L+1/10*m/A/L*12*Iz, 0, 0, 0, 11/210*m*L^2+1/120*m/A*12*Iz, 0, 9/70*m*L-1/10*m/A/L*12*Iz, 0, 0, 0, -13/420*m*L^2+1/120*m/A*12*Iz
0, 0, 13/35*m*L+1/10*m/A/L*12*Iy, 0, -11/210*m*L^2-1/120*m/A*12*Iy, 0, 0, 0, 9/70*m*L-1/10*m/A/L*12*Iy, 0, 13/420*m*L^2-1/120*m/A*12*Iy, 0
0, 0, 0, 1/3*m/A*Jx*L, 0, 0, 0, 0, 0, 1/6*m/A*Jx*L, 0, 0
0, 0, -11/210*m*L^2-1/120*m/A*12*Iy, 0, 1/105*m*L^3+1/90*m/A*L*12*Iy, 0, 0, 0, -13/420*m*L^2+1/120*m/A*12*Iy, 0, -1/140*m*L^3-1/360*m/A*L*12*Iy, 0
0, 11/210*m*L^2+1/120*m/A*12*Iz, 0, 0, 0, 1/105*m*L^3+1/90*m/A*L*12*Iz, 0, 13/420*m*L^2-1/120*m/A*12*Iz, 0, 0, 0, -1/140*m*L^3-1/360*m/A*L*12*Iz
1/6*m*L, 0, 0, 0, 0, 0, 1/3*m*L, 0, 0, 0, 0, 0
0, 9/70*m*L-1/10*m/A/L*12*Iz, 0, 0, 0, 13/420*m*L^2-1/120*m/A*12*Iz, 0, 13/35*m*L+1/10*m/A/L*12*Iz, 0, 0, 0, -11/210*m*L^2-1/120*m/A*12*Iz
0, 0, 9/70*m*L-1/10*m/A/L*12*Iy, 0, -13/420*m*L^2+1/120*m/A*12*Iy, 0, 0, 0, 13/35*m*L+1/10*m/A/L*12*Iy, 0, 11/210*m*L^2+1/120*m/A*12*Iy, 0
0, 0, 0, 1/6*m/A*Jx*L, 0, 0, 0, 0, 0, 1/3*m/A*Jx*L, 0, 0
0, 0, 13/420*m*L^2-1/120*m/A*12*Iy, 0, -1/140*m*L^3-1/360*m/A*L*12*Iy, 0, 0, 0, 11/210*m*L^2+1/120*m/A*12*Iy, 0, 1/105*m*L^3+1/90*m/A*L*12*Iy, 0
0, -13/420*m*L^2+1/120*m/A*12*Iz, 0, 0, 0, -1/140*m*L^3-1/360*m/A*L*12*Iz, 0, -11/210*m*L^2-1/120*m/A*12*Iz, 0, 0, 0, 1/105*m*L^3+1/90*m/A*L*12*Iz];
% 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];
me=t*me*transpose(t);
return
function AssembleMassMatrix(ie,me)
global element M
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;
Q=(element(ie,i)-1)*6+p;
N=(element(ie,j)-1)*6+q;
M(Q,N)=M(Q,N)+me(m,n);
end
end
end
end
return
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -