Matlab中求解结构动力响应问题
现有一简支梁，在外界加速度作用下振动，需要求解在某一加速度下梁沿横向变形情况，程序如下，不知道哪里有问题，求解后的值与理论值相差很大。望各位高手赐教！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*ell4*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^20 -22*ell4*ell^2;]
me=T*me*transpose(T);
fqi=-mass/2;mzi=-mass*ell/12;
fe=[ 0 fqi mzi 0fqi -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 ];
% 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法
= 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=;
A=;
% 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')
回复 楼主 bestingmoon 的帖子
大伙没有你说的文献，就算有也可能和你专业背景不同帮不上你。最可行的办法就是你自己在程序中设置断点，跟踪程序的运行情况，随时检查程序运行结果是否正确，直到找到错误。
回复 楼主 bestingmoon 的帖子
非常同意sogooda主任的看法!或者楼主需说明更详细背景资料, 譬如解什麽? 过程? 理论值? 楼主解的值? ...
个人认为适当的转换描述, 或许有机会有人帮忙! 基于MATLAB和传递函数的结构动力响应计算
在中国期刊数据库里面有这方面的文章，本人跟你的专业背景不一样，你自己看看。题目如上
页:
[1]