求柱的固有频率和欧拉屈曲荷载

 

clear
clc

disp('please wait!!!!!!-The job is under run')

% Discretizing the Beam

nel=50;                 % number of elements
nnel=2;                 % number of nodes per element
ndof=2;                 % number of dofs per node
nnode=(nnel-1)*nel+1;   % total number of nodes in system
sdof=nnode*ndof;        % total system dofs 

% Material properties
E=2.1*10^11;            % Youngs modulus
I=2003.*10^-8;          % moment of inertia of cross-section
mass = 61.3;            % mass density of the beam
tleng = 7.;             % total length of the beam
leng = tleng/nel;       % uniform mesh (equal size of elements)
lengthvector = 0:leng:tleng ;
% Boundary Conditions
bc = 'c-f' ;             % clamped-free
%bc = 'c-c' ;            % clamped-clamped
%bc = 'c-s' ;            % clamped-supported
%bc = 's-s' ;            % supported-supported

kk=zeros(sdof,sdof);    % initialization of system stiffness matrix
kkg=zeros(sdof,sdof);   % initialization of system geomtric stiffness matrix 
mm=zeros(sdof,sdof);    % initialization of system mass matrix 
index=zeros(nel*ndof,1);  % initialization of index vector

for iel=1:nel           % loop for the total number of elements

index=elementdof(iel,nnel,ndof);  % extract system dofs associated with element

[k,kg,m]=beam(E,I,leng,mass); % compute element stiffness,geometric
                                   % stiffness & mass matrix
                                   
kk=assembel(kk,k,index); % assemble element stiffness matrices into system matrix
kkg=assembel(kkg,kg,index); % assemble geometric stiffness matrices into system matrix
mm=assembel(mm,m,index); % assemble element mass matrices into system matrix

end
%
% Applying the Boundary conditions
[nbcd,bcdof] = BoundaryConditions(sdof,bc); % Reducing the matrix size
[kk,mm] = constraints(kk,mm,bcdof) ;
[kk,kkg] = constraints(kk,kkg,bcdof) ;
%
% Natural frequencies and Buckling load
[vecfreq freq]=eig(kk,mm);   % solve the eigenvalue problem for Natural Frequencies
freq = diag(freq) ;
freq=sqrt(freq);   % UNITS :rad per sec
freqHz = freq/(2*pi) ; % UNITS : Hertz
%
[vecebl ebl] = eig(kk,kkg);  % solve the eigenvalue problem for Buckling Loads
ebl = diag(ebl) ;
%
% Plot Mode Shapes
h = figure ;
set(h,'name','Mode Shapes of Beam in rad/s','numbertitle','off')
PlotModeShapes(vecfreq,freq,lengthvector,nbcd)
h = figure ;
set(h,'name','Buckling Mode shape in N','numbertitle','off')
PlotModeShapes(vecebl,ebl,lengthvector,nbcd)
%
% Theoretical Natural Frequencies
 [thfreq,thfreqHz,pcr] =  theory(bc,E,I,mass,tleng) ;

% Code validitation
theory = thfreq(1:3) ;
fem = freq(nbcd+1:nbcd+3);
error = (fem-theory)./theory*100;
compare = [theory fem error] ;
disp('First three natural frequencies (rad/sec)')
disp('theory        fem       error%')
disp('---------------------------------' )
disp(compare)
%
theory = pcr ;
fem = ebl(1);
error = (fem-theory)./theory*100 ;
compare = [theory fem error] ;
disp('Euler Buckling load (N)')
disp('theory        fem       error%')
disp('---------------------------------' )
disp(compare)

d146