结构动力学matlab代码

【核心代码】

%The author: wangjinghui
%Date:2021.04

clear all;
clc;
%Beam_InputData544;

Opt_beam = 1;
Opt_mass = 2;
Opt_section = 1;
Opt_graphics1 = 1;
Opt_graphics2 = 2;






Prob_title= 'static  analysisi  of  a  3-D  frame  structure';
  No_el=12;
  No_nel=2;
  No_dof=6;
  No_node=8;
  Sys_dof=No_node*No_dof;
  
  
  prop(1)=2.1e11;
  prop(2)=0.3;
  prop(3)=7860;
  prop(4)=0.01;
  prop(5)=0.01;
  prop(6)=0.01*0.01;
  prop(7)=0.01*0.01^3/12;
  prop(8)=0.01*0.01^3/12;
  
  prop(9)=1;
   
  prop(10)=1;
  
  Opt_section=1;
  
  
  gcoord=[1   0,0,0;
                 2   4,0,0;
                 3   4,0,0;
                 4   4,0,0;
                 5   4,0,0;
                 6   4,0,0;
                 7   4,0,0;
                 8   4,0,0];
             
             
gcoord=[gcoord(:,2:4)];
  
      al=0;
      node=[1   1 , 5;
                  2   2 , 6;
                  3   3 , 7;
                  4   4 , 8;
                  5   5 , 6;
                  6   6 , 7;
                  7   7 , 8;
                  8   8 , 5;
                  9   1 , 6;
                  10   2 , 7;
                  11   3 , 8;
                  12   4 , 5];
                 
nodes=[nodes(:,2:3)];
P=[-100000,7,2];

ConNode =[1 , 1 , 1 , 1 , 1 , 1 ;
                   2 , 1 , 1 , 1 , 1 , 1 ;
                   3 , 1 , 1 , 1 , 1 , 1 ;
                   4 , 1 , 1 , 1 , 1 , 1 ];
               
               
               ConVal =[1 , 0 , 0 , 0 , 0 , 0 ;
                   2 , 0 , 0 , 0 , 0 , 0 ;
                   3 , 0 , 0 , 0 , 0 , 0 ;
                   4 , 0 , 0 , 0 , 0 , 0 ];
               








k = zeros(No_nel*No_dof,No_nel*No_dof);
m = zeros(No_nel*No_dof,No_nel*No_dof);

kk = zeros(Sys_dof,Sys_dof);
mm = zeros(Sys_dof,Sys_dof);

ff = zeros(Sys_dof,1);
bcdof = zeros(Sys_dof,1);

index = zeros(No_nel*No_dof,1);



[n1,n2]=size(ConNode);

for ni = 1:n1
  ki=ConNode(ni,1);
  bcdof((ki-1)*No_dof:ki*No_dof)=ConNode(ni,2:No_dof 1);
  
  bcval((ki-1)*No_dof:ki*No_dof)=ConVal(ni,2:No_dof 1);
end

[n1,n2]=size(P);


for  ni=i:n1
  ff(No_dof*(P(ni,2)-1) P(ni,3))=P(ni,1);
end


for iel=1:No_el
  nd(1)=nodes(iel,1);
  nd(2)=nodes(iel,2);
  x(1)=gcoord(nd(1),1);
  y(1)=gcoord(nd(1),2);
  z(1)=gcoord(nd(1),3);
  

  x(2)=gcoord(nd(2),1);
  y(2)=gcoord(nd(2),2);
  z(2)=gcoord(nd(2),3);
  
  
  leng=sqrt((x(2)-x(1))^2 (y(2)-y(1))^2 (z(2)-z(1))^2);
  
  xi(1)=(x(2)-x(1))/leng;
  xi(2)=(y(2)-y(1))/leng;
  xi(3)=(z(2)-z(1))/leng;
  
  
  if Opt_beam==1
    [k,m]=FrameElement31(prop,leng,xi,al,Opt_section,Opt_mass);
  elseif Opt_beam==2
    [k,m]=FrameElement32(prop,leng,xi,al,Opt_section,Opt_mass);
  end
  
  index=femEldof(nd,No_nel,No_dof);
  
  kk=femAssemble(kk,k,index);
  mm=femAssemble(mm,m,index);
  end
  
  
  [kk0 , mm0 , ff0 ,bcdof0 , bcval0 , sdof0]=kkCheck1(kk , mm , ff , bcdof , bcval);
  
  [kk1 , mm1 , ff1 , sdof1]=femApplybc1(kk0 , mm0 , ff0 , bcdof0 , bcval0);
  
  displmt = kk1\ff1;
  No_dof1 = No_dof/(Sys_dof/sdof1);
  for ii=1:No_node
      for ij=1:No_dof1
          displmtnode(ii ,  ij)=displmt((ii-1)*No_dof1 ij,1);
      end
  end
  
  if Opt_graphics1==1
      for iel=1:No_el
          nd(1)=nodes(iel,1);
          nd(2)=nodes(iel,2);
          
          x(1)=gcoord(nd(1),1);
          y(1)=gcoord(nd(1),2);
          z(1)=gcoord(nd(1),3);
      
          
          x(2)=gcoord(nd(2),1);
          y(2)=gcoord(nd(2),2);
          z(2)=gcoord(nd(2),3);
          
          figure(1)
          plot3(x,y,z),xlabel('x'),ylabel('y'),zlabel('z'),hold on;
      end
      
      
      if  Opt_graphics2~=1 
          title('nodal connectivity of elementd'); hold off;
      end
  end
  
  if  Opt_graphics2==1
      gcoordA= gcoord 10*displmtnode( : ,1:3);
      for iel=1:No_el
          nd(1)=nodes(iel,1);
          nd(2)=nodes(iel,2);
          
          
          x(1)=gcoordA(nd(1),1);
          y(1)=gcoordA(nd(1),2);
          z(1)=gcoordA(nd(1),3);
      
          
          x(2)=gcoordA(nd(2),1);
          y(2)=gcoordA(nd(2),2);
          z(2)=gcoordA(nd(2),3);
          
          figure(1)
          plot3(x,y,z,'r'), title('Deflection of the structure'),  hold on;
      end
      
      hold off
  end
  
  disp('The calculation is use of')
  
  
  if Opt_beam==1
      disp('Euler-Bernoulli beam element')
  elseif Opt_beam==2
      diso('Timoshenko beam element')
  end
  disp('  ')
  disp('      numrical solution        ')
  disp('      node x  y  z  牟x   牟y   牟z    ')
  num  = 1:1:No_node;
  displmts=[num' displmtnode]