Matlab有限元结构动力学分析（源程序）

%-------------------------------------------------------------------------% % Example 5.4.7 % To solve transient response of a 2-d frame structure. % The solution methods are: 1) central difference scheme. 3) Houbolt integration scheme. % 4) Wilson ?? integration scheme. 5) Newmark integration scheme % nodal dof: {u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z2} % Problem description % Find the response of a frame structure which is made of three beams of lengths of 4 m, % 3 m and 4 m, respectively. All beams have cross- section of 0.10 m height by 0.05 m % width. The elastic modulus is 2.10x10^11 Pa. The frame is subjected to an impulse load % of amplitude 500 N in the middle of the top beam. One end of the each vertical beam is % fixed. (see Fig. 5-9 for the element discretization) % Variable descriptions % k, m - element stiffness matrix and mass matrix % kk, mm - system stiffness matrix and mass matrix % ff - system force vector % index - a vector containing system dofs associated with each element % bcdof - a vector containing dofs associated with boundary conditions % bcval - a vector containing boundary condition values associated with the dofs in 'bcdof' % dsp - displacement matrix % vel - velocity matrix % acc - acceleartion matrix %-------------------------------------------------------------------------- % (0) input control data %-------------------------------------------------------------------------- clear; clc; Beam_InputData547; % import the input data for the information of % nodes, elements, loads, constraints and materials Opt_beam=1; % option for type of the beam % =1 Euler Bernoulli beam % =2 Timoshenko beam Opt_mass=2; % option for mass matrix % =1 consistent mass matrix % =2 lumped mass matrix Opt_section=1; % option for type of cross-section % = 1 rectangular cross-section % = 2 circular cross-section TypeMethod=1; % option for selecting the solution method % = 1 central difference scheme % = 3 Houbolt integration scheme % = 4 Wilson integration scheme % = 5 Newmark integration scheme Typeload=1; % option for selecting the load type % = 1 impulse load % = 2 step load % = 3 Harmonic load dt=0.00001; % time step size ti=0; % initial time tf=0.200; % final time nt=fix((tf-ti)/dt); % number of time steps tt=ti:dt:ti+nt*dt; % generate time samples vector ac=0.00002; bc=0.00008; % Parameters for proportional damping al=0; % angle between the reference coordinate system and % the local coordinate system for the space element %-------------------------------------------------------------------------- % (1) initialization of matrices and vectors to zero %-------------------------------------------------------------------------- k=zeros(No_nel*No_dof,No_nel*No_dof); % element stiffness matrix m=zeros(No_nel*No_dof,No_nel*No_dof); % element mass matrix kk=zeros(Sys_dof,Sys_dof); % initialization of system stiffness matrix mm=zeros(Sys_dof,Sys_dof); % initialization of system mass matrix ff=zeros(Sys_dof,1); % initialization of system force vector bcdof=zeros(Sys_dof,1); % initializing the vector bcdof bcval=zeros(Sys_dof,1); % initializing the vector bcval index=zeros(No_nel*No_dof,1); % initialization of index vector %-------------------------------------------------------------------------- % (2) calculation of constraints %-------------------------------------------------------------------------- [n1,n2]=size(ConNode); % calculate the constrained dofs for ni=1:n1 ki=ConNode(ni,1); bcdof((ki-1)*No_dof+1:ki*No_dof)=ConNode(ni,2:No_dof+1); % the code for constrained dofs bcval((ki-1)*No_dof+1:ki*No_dof)=ConVal(ni,2:No_dof+1); % the value at constrained dofs end %-------------------------------------------------------------------------- % (3) applied nodal loads %-------------------------------------------------------------------------- [n1,n2]=size(P); for ni=1:n1 ff(No_dof*(P(ni,2)-1)+P(ni,3))=P(ni,1); end %-------------------------------------------------------------------------- % (4) calculate the element matrices and assembling %-------------------------------------------------------------------------- for iel=1:No_el % loop for the total number of elements nd(1)=iel; % starting node number for element 'iel' nd(2)=iel+1; % ending node number for element 'iel' x(1)=gcoord(nd(1),1); y(1)=gcoord(nd(1),2); z(1)=gcoord(nd(1),3); % coordinate of 1st node x(2)=gcoord(nd(2),1); y(2)=gcoord(nd(2),2); z(2)=gcoord(nd(2),3); % coordinate of 2nd node leng=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2+(z(2)-z(1))^2); xi(1)=(x(2)-x(1))/leng; % cos of the local x-axis and system x-axis xi(2)=(y(2)-y(1))/leng; % cos of the local x-axis and system y-axis xi(3)=(z(2)-z(1))/leng; % cos of the local x-axis and system y-axis if Opt_beam==1 [k,m]=FrameElement31(prop,leng,xi,al,Opt_section,Opt_mass); % compute element matrix for Euler-Bernoulli beam elseif Opt_beam==2 [k,m]=FrameElement32(prop,leng,xi,al,Opt_section,Opt_mass); % compute element matrix for Timoshenko beam end index=femEldof(nd,No_nel,No_dof); % extract system dofs associated with element kk=femAssemble1(kk,k,index); % assemble system stiffness matrix mm=femAssemble1(mm,m,index); % assemble system mass matrix end %-------------------------------------------------------------------------- % (5) apply constraints and solve the matrix equation %-------------------------------------------------------------------------- [kk0,mm0,ff0,bcdof0,bcval0,sdof0]=kkCheck1(kk,mm,ff,bcdof,bcval); % check the zero main elements in kk and eliminate the rows % and columns in equation associated with the zero main elements [kk1,mm1,ff1,sdof1]=femApplybc1(kk0,mm0,ff0,bcdof0,bcval0); % apply boundary conditions to the system equation [kk2,mm2,ff2,sdof2]=bcCheck1(kk0,mm0,ff0,bcdof0,bcval0);