|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
现有一简支梁,在外界加速度作用下振动,需要求解在某一加速度下梁沿横向变形情况,程序如下,不知道哪里有问题,求解后的值与理论值相差很大。望各位高手赐教!
clc
clear
format long
%
% Input data
%
% default data!!
L = 100*10^-6; % beam span(meter)
b = 5*10^-7; % beam width(meter)
h = 15*10^-7; % beam thickness(meter)
E = 92.e+9; % young's modulus(Pascal) of multiple beam
Iz = b*h^3/12; % second area moment
EI = E*Iz; % bending rigidity
V=b*h*L; % volumn of beam
rho = 2500; % density (kg/meter^3) of multiple beam
m_total = rho*V; % total mass of beam
global gElement gK
Wmax = 2000;% Maximum frequency Wmax
I = 2; % Station # of Excitation
J = 5; % Station # of Response
disp(' ');
disp('Welcome to vibration analysis of a beam by the finite element method');
nelem =input('number of elements to model the beam = ');
% begin computations
%
node_g = nelem+1; % total number of global nodes
node_e = 2; % number of nodes in each element
dof_node = 3; % degrees of freedom per node
dof_elem = dof_node*node_e; % total degrees of freedom per elemental
dof_g = dof_node*node_g; % total global degree of freedom
%
% element stiffness matrix
%
ell = L/nelem; % elemental beam length
kr = EI/ell^3;
ke=[
E*b*h/L 0 0 E*b*h/L 0 0;
0 12*EI/ell^3 0 0 12*EI/ell^3 0;
0 6*EI/ell^2 4*EI/ell 0 6*EI/ell^2 4*EI/ell;
-E*b*h/L 0 0 E*b*h/L 0 0;
0 -12*EI/ell^3 -6*EI/ell^2 0 -12*EI/ell^3 0
0 6*EI/ell^2 2*EI/ell 0 -6*EI/ell^2 4*EI/ell ]
% count transform matrix
T=[
1 0 0 0 0 0;
0 1 0 0 0 0;
0 0 1 0 0 0;
0 0 0 1 0 0;
0 0 0 0 1 0;
0 0 0 0 0 1;]
ke=T*ke*transpose(T);
% lumped element mass matrix
mass=rho*b*h*ell*9.8; % mass of element
me=mass/420*[
140 0 0 70 0 0 ;
0 156 22*ell 0 54 -13*ell;
0 22*ell 4*ell^2 0 13*ell -3*ell^2;
70 0 0 140 0 0;
0 54 13*ell 0 156 -22*ell;
0 -13*ell -3*ell^2 0 -22*ell 4*ell^2;]
me=T*me*transpose(T);
fqi=-mass/2;mzi=-mass*ell/12;
fe=[ 0 fqi mzi 0 fqi -mzi]
fe=T*fe';
gF=zeros(node_g*3,1);
%count global mass and stiffness matrix
gK=zeros(node_g*3,node_g*3);
gM=zeros(node_g*3,node_g*3);
gElement=zeros(nelem,2)
for i_g=1:nelem
gElement(i_g,1)=gElement(i_g,1)+i_g;
gElement(i_g,2)=gElement(i_g,2)+i_g+1;
end
for ie=1:nelem
for i=1:1:2
for j=1:1:2
for p=1:1:3
for q=1:1:3
m=(i-1)*3+p;
n=(j-1)*3+q;
M=(gElement(ie,i)-1)*3+p;
N=(gElement(ie,j)-1)*3+q;
gK(M,N)=gK(M,N)+ke(m,n);
gM(M,N)=gM(M,N)+me(m,n);
end
end
end
end
end
for ie=1:nelem
i=gElement(ie,1)
j=gElement(ie,2)
gF( (i-1)*3+1 : (i-1)*3+3 ) = gF( (i-1)*3+1 : (i-1)*3+3 ) + fe( 1:3 ) ;
gF( (j-1)*3+1 : (j-1)*3+3 ) = gF( (j-1)*3+1 : (j-1)*3+3 ) + fe( 4:6 ) ;
end
% 约束条件 ,限制了三个自由度
% 节点号 自由度号 约束值
gBC1 = [ 1, 1, 0.0
1, 2, 0.0
1, 3, 0.0
node_g, 1, 0.0
node_g, 2, 0.0
node_g, 3, 0.0 ];
% 处理约束条件,修改刚度矩阵和节点力向量。采用乘大数法
[bc_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
gF = gF - gBC1(ibc,3) * gK(:,m) ;
gF(m) = gBC1(ibc,3) ;
gK(:,m) = zeros( node_g*3, 1 ) ;
gK(m,:) = zeros( 1, node_g*3 ) ;
gK(m,m) = 1.0 ;
end
%n = gBC1(ibc, 1 ) ;
%d = gBC1(ibc, 2 ) ;
%m = (n-1)*3 + d ;
%gF(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;
%gK(m,m) = gK(m,m) * 1e15 ;
end
% count damnp index
b=[0.1;0.1];
A=[1/(2*10*2*pi) 2*pi*10/2; 1/(2*5000*2*pi) 2*pi*5000/2];
% A*damp=b
damp=A\b;
alphadM=damp(1)
betadK=damp(2)
gC=alphadM*gM+betadK*gK
% count d v a
omiga=10*2*pi %frequency
KK=-omiga^2*gM+omiga*gC+gK
delta=KK\gF
for di=1:1:nelem+1
d1(di)=delta((di-1)*3+2);
d2(di)=di
end
plot(d2,d1,'r') |
|