<P>function beam3d<BR>% 空間彈性梁單元結構－靜力分析<BR>global node element material K knode<BR>file_in=input('input file name: ','s');<BR>file_out=input('output file name: ','s');<BR>fid_in=fopen(file_in,'r');<BR>node_number=fscanf(fid_in,'%d',1);<BR>node=zeros(node_number,3);<BR>for i=1:node_number<BR> nn=fscanf(fid_in,'%d',1);<BR> node(i,:)=fscanf(fid_in,'%f',[1,3]);<BR>end<BR>knode_number=fscanf(fid_in,'%d',1);<BR>knode=zeros(knode_number,3);<BR>for i=1:knode_number<BR> nn=fscanf(fid_in,'%d',1);<BR> knode(i,:)=fscanf(fid_in,'%f',[1,3]);<BR>end<BR>element_number=fscanf(fid_in,'%d',1);<BR>element=zeros(element_number,4);<BR>for i=1:element_number<BR> ne=fscanf(fid_in,'%d',1);<BR> element(i,:)=fscanf(fid_in,'%d',[1,4]);<BR>end<BR>material_number=fscanf(fid_in,'%d',1);<BR>material=zeros(material_number,7);<BR>for i=1:material_number<BR> nm=fscanf(fid_in,'%d',1);<BR> material(i,:)=fscanf(fid_in,'%f',[1,7]);<BR>end<BR>bc_number=fscanf(fid_in,'%d',1);<BR>bc=zeros(bc_number,3);<BR>for i=1:bc_number<BR> nb=fscanf(fid_in,'%d',1);<BR> bc(i,1)=fscanf(fid_in,'%d',1);<BR> bc(i,2)=fscanf(fid_in,'%d',1);<BR> bc(i,3)=fscanf(fid_in,'%f',1);<BR>end<BR>nf_number=fscanf(fid_in,'%d',1);<BR>nf=zeros(nf_number,3);<BR>for i=1:nf_number<BR> nnf=fscanf(fid_in,'%d',1);<BR> nf(i,1)=fscanf(fid_in,'%d',1);<BR> nf(i,2)=fscanf(fid_in,'%d',1);<BR> nf(i,3)=fscanf(fid_in,'%f',1);<BR>end<BR>fclose(fid_in);</P>
<P>K=sparse(node_number*6,node_number*6);<BR>f=sparse(node_number*6,1);<BR>for ie=1:element_number<BR> k=StiffnessMatrix(ie);<BR> AssembleStiffnessMatrix(ie,k);<BR>end</P>
<P>for ibc=1:bc_number<BR> n=bc(ibc,1);<BR> d=bc(ibc,2);<BR> m=(n-1)*6+d;<BR> f(m)=bc(ibc,3)*K(m,m)*1E15;<BR> K(m,m)=K(m,m)*1E15;<BR>end</P>
<P>for inf=1:nf_number<BR> n=nf(inf,1);<BR> d=nf(inf,2);<BR> f((n-1)*6+d)=nf(inf,3);<BR>end<BR>delta=K\f;<BR>%disp(delta);<BR>fid_out=fopen(file_out,'w');<BR>for i=1:node_number<BR> 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));<BR>end<BR>for i=1:element_number<BR> fprintf(fid_out,'\r\nElement %d\r\n\r\n',i);<BR> for j=1:2<BR> fprintf(fid_out,' Node %d\r\n',j);<BR> for p=1:6<BR> fprintf(fid_out,' F%d= %.5f\r\n',p,getforce(delta,i,j,p));<BR> end<BR> end<BR>end</P>
<P>fclose(fid_out);<BR>return;</P>
<P>function k=StiffnessMatrix(ie)<BR>global node element material K knode<BR>k=zeros(12,12);<BR>E=material(element(ie,3),1);<BR>G=material(element(ie,3),2);<BR>Iy=material(element(ie,3),3);<BR>Iz=material(element(ie,3),4);<BR>Jx=material(element(ie,3),5);<BR>A=material(element(ie,3),6);<BR>xi=node(element(ie,1),1);<BR>yi=node(element(ie,1),2);<BR>zi=node(element(ie,1),3);<BR>xj=node(element(ie,2),1);<BR>yj=node(element(ie,2),2);<BR>zj=node(element(ie,2),3);<BR>xk=knode(element(ie,4),1);<BR>yk=knode(element(ie,4),2);<BR>zk=knode(element(ie,4),3);<BR>L=((xj-xi)^2+(yj-yi)^2+(zj-zi)^2)^(1/2);<BR>Lk=((xk-xi)^2+(yk-yi)^2+(zk-zi)^2)^(1/2);<BR>l1=(xj-xi)/L;<BR>m1=(yj-yi)/L;<BR>n1=(zj-zi)/L;<BR>g1=(xk-xi)/Lk;<BR>g2=(yk-yi)/Lk;<BR>g3=(zk-zi)/Lk;<BR>s=((m1*g3-n1*g2)^2+(n1*g1-l1*g3)^2+(l1*g2-m1*g1)^2)^(1/2);</P>
<P>k=[E*A/L 0 0 0 0 0 -E*A/L 0 0 0 0 0<BR> 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<BR> 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<BR> 0 0 0 G*Jx/L 0 0 0 0 0 -G*Jx/L 0 0<BR> 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<BR> 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<BR> -E*A/L 0 0 0 0 0 E*A/L 0 0 0 0 0<BR> 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 <BR> 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<BR> 0 0 0 -G*Jx/L 0 0 0 0 0 G*Jx/L 0 0<BR> 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<BR> 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];<BR> <BR>% k點取自X-Y平面內<BR>t11=l1;<BR>t12=(g1-l1*(l1*g1+m1*g2+n1*g3))/s;<BR>t13=(m1*g3-n1*g2)/s;<BR>t21=m1;<BR>t22=(g2-m1*(l1*g1+m1*g2+n1*g3))/s;<BR>t23=(n1*g1-l1*g3)/s;<BR>t31=n1;<BR>t32=(g3-n1*(l1*g1+m1*g2+n1*g3))/s;<BR>t33=(l1*g2-m1*g1)/s;</P>
<P>t=[t11 t12 t13 0 0 0 0 0 0 0 0 0<BR> t21 t22 t23 0 0 0 0 0 0 0 0 0<BR> t31 t32 t33 0 0 0 0 0 0 0 0 0<BR> 0 0 0 t11 t12 t13 0 0 0 0 0 0<BR> 0 0 0 t21 t22 t23 0 0 0 0 0 0<BR> 0 0 0 t31 t32 t33 0 0 0 0 0 0<BR> 0 0 0 0 0 0 t11 t12 t13 0 0 0<BR> 0 0 0 0 0 0 t21 t22 t23 0 0 0<BR> 0 0 0 0 0 0 t31 t32 t33 0 0 0<BR> 0 0 0 0 0 0 0 0 0 t11 t12 t13<BR> 0 0 0 0 0 0 0 0 0 t21 t22 t23<BR> 0 0 0 0 0 0 0 0 0 t31 t32 t33];<BR>k=t*k*transpose(t);<BR>return</P>
<P>function AssembleStiffnessMatrix(ie,k)<BR>global element K<BR>for i=1:2<BR> for j=1:2<BR> for p=1:6<BR> for q=1:6<BR> m=(i-1)*6+p;<BR> n=(j-1)*6+q;<BR> M=(element(ie,i)-1)*6+p;<BR> N=(element(ie,j)-1)*6+q;<BR> K(M,N)=K(M,N)+k(m,n);<BR> end<BR> end<BR> end<BR>end<BR>return</P>
<P><BR>function shearf=getforce(b,ie,i,d) %獲得桿端內力 位移 單元 桿端 自由度<BR>global node element material K knode<BR>k=zeros(12,12);<BR>E=material(element(ie,3),1);<BR>G=material(element(ie,3),2);<BR>Iy=material(element(ie,3),3);<BR>Iz=material(element(ie,3),4);<BR>Jx=material(element(ie,3),5);<BR>A=material(element(ie,3),6);<BR>xi=node(element(ie,1),1);<BR>yi=node(element(ie,1),2);<BR>zi=node(element(ie,1),3);<BR>xj=node(element(ie,2),1);<BR>yj=node(element(ie,2),2);<BR>zj=node(element(ie,2),3);<BR>xk=knode(element(ie,4),1);<BR>yk=knode(element(ie,4),2);<BR>zk=knode(element(ie,4),3);<BR>L=((xj-xi)^2+(yj-yi)^2+(zj-zi)^2)^(1/2);<BR>Lk=((xk-xi)^2+(yk-yi)^2+(zk-zi)^2)^(1/2);<BR>l1=(xj-xi)/L;<BR>m1=(yj-yi)/L;<BR>n1=(zj-zi)/L;<BR>g1=(xk-xi)/Lk;<BR>g2=(yk-yi)/Lk;<BR>g3=(zk-zi)/Lk;<BR>s=((m1*g3-n1*g2)^2+(n1*g1-l1*g3)^2+(l1*g2-m1*g1)^2)^(1/2);</P>
<P>k=[E*A/L 0 0 0 0 0 -E*A/L 0 0 0 0 0<BR> 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<BR> 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<BR> 0 0 0 G*Jx/L 0 0 0 0 0 -G*Jx/L 0 0<BR> 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<BR> 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<BR> -E*A/L 0 0 0 0 0 E*A/L 0 0 0 0 0<BR> 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 <BR> 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<BR> 0 0 0 -G*Jx/L 0 0 0 0 0 G*Jx/L 0 0<BR> 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<BR> 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];<BR> <BR>% k點取自X-Y平面內<BR>t11=l1;<BR>t12=(g1-l1*(l1*g1+m1*g2+n1*g3))/s;<BR>t13=(m1*g3-n1*g2)/s;<BR>t21=m1;<BR>t22=(g2-m1*(l1*g1+m1*g2+n1*g3))/s;<BR>t23=(n1*g1-l1*g3)/s;<BR>t31=n1;<BR>t32=(g3-n1*(l1*g1+m1*g2+n1*g3))/s;<BR>t33=(l1*g2-m1*g1)/s;</P>
<P>t=[t11 t12 t13 0 0 0 0 0 0 0 0 0<BR> t21 t22 t23 0 0 0 0 0 0 0 0 0<BR> t31 t32 t33 0 0 0 0 0 0 0 0 0<BR> 0 0 0 t11 t12 t13 0 0 0 0 0 0<BR> 0 0 0 t21 t22 t23 0 0 0 0 0 0<BR> 0 0 0 t31 t32 t33 0 0 0 0 0 0<BR> 0 0 0 0 0 0 t11 t12 t13 0 0 0<BR> 0 0 0 0 0 0 t21 t22 t23 0 0 0<BR> 0 0 0 0 0 0 t31 t32 t33 0 0 0<BR> 0 0 0 0 0 0 0 0 0 t11 t12 t13<BR> 0 0 0 0 0 0 0 0 0 t21 t22 t23<BR> 0 0 0 0 0 0 0 0 0 t31 t32 t33];</P>
<P>globdelta=zeros(12,1);</P>
<P>for ij=1:2<BR> nodenum=element(ie,ij);<BR> for p=1:6<BR> globdelta((ij-1)*6+p,1)=b((nodenum-1)*6+p,1); <BR> end<BR>end</P>
<P>delta=transpose(t)*globdelta;</P>
<P>f=k*delta;<BR>shearf=f((i-1)*6+d);<BR>return</P>
<P> </P>