beam.m

matlab有限元程序

function beam
global node element material K
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,2);
for i=1:node_number
    nn=fscanf(fid_in,'%d',1);
    node(i,:)=fscanf(fid_in,'%f',[1,2]);
end
element_number=fscanf(fid_in,'%d',1);
element=zeros(element_number,3);
for i=1:element_number
    ne=fscanf(fid_in,'%d',1);
    element(i,:)=fscanf(fid_in,'%d',[1,3]);
end
material_number=fscanf(fid_in,'%d',1);
material=zeros(material_number,3);
for i=1:material_number
    nm=fscanf(fid_in,'%d',1);
    material(i,:)=fscanf(fid_in,'%f',[1,3]);
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*3,node_number*3);
f=sparse(node_number*3,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)*3+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)*3+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 -- u= %.4e  v= %.4e  sita= %.4e\r\n',i,delta((i-1)*3+1),delta((i-1)*3+2),delta((i-1)*3+3));
end
fclose(fid_out);
return;

function k=StiffnessMatrix(ie)
global node element material K
k=zeros(6,6);
E=material(element(ie,3),1);
I=material(element(ie,3),2);
A=material(element(ie,3),3);
xi=node(element(ie,1),1);
yi=node(element(ie,1),2);
xj=node(element(ie,2),1);
yj=node(element(ie,2),2);
L=((xj-xi)^2+(yj-yi)^2)^(1/2);
l1=(xj-xi)/L;
m1=(yj-yi)/L;

k=[E*A/L             0                    0                 -E*A/L             0                     0
     0                   12*E*I/L^3       6*E*I/L^2      0                   -12*E*I/L^3       6*E*I/L^2
     0                   6*E*I/L^2         4*E*I/L         0                   -6*E*I/L^2         2*E*I/L
    -E*A/L            0                     0                 E*A/L              0                    0
     0                  -12*E*I/L^3       -6*E*I/L^2    0                    12*E*I/L^3      -6*E*I/L^2
     0                   6*E*I/L^2         2*E*I/L         0                   -6*E*I/L^2         4*E*I/L];
t=[l1    -m1     0        0        0       0
    m1   l1       0        0        0       0
    0      0        1       0        0       0
    0      0        0       l1      -m1    0
    0      0        0       m1     l1      0
    0      0        0       0       0        1];
k=t*k*transpose(t);
return

function AssembleStiffnessMatrix(ie,k)
global element K
for i=1:2
	for j=1:2
		for p=1:3
			for q=1:3
				m=(i-1)*3+p;
				n=(j-1)*3+q;
				M=(element(ie,i)-1)*3+p;
				N=(element(ie,j)-1)*3+q;
				K(M,N)=K(M,N)+k(m,n);
			end
		end
	end
end
return