Matlab中求解结构动力响应问题

bestingmoon

现有一简支梁，在外界加速度作用下振动，需要求解在某一加速度下梁沿横向变形情况，程序如下，不知道哪里有问题，求解后的值与理论值相差很大。望各位高手赐教！
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')

sogooda

大伙没有你说的文献，就算有也可能和你专业背景不同帮不上你。
最可行的办法就是你自己在程序中设置断点，跟踪程序的运行情况，随时检查程序运行结果是否正确，直到找到错误。

ChaChing

非常同意sogooda主任的看法!
或者楼主需说明更详细背景资料, 譬如解什麽? 过程? 理论值? 楼主解的值? ...
个人认为适当的转换描述, 或许有机会有人帮忙!

科技在线

基于MATLAB和传递函数的结构动力响应计算

在中国期刊数据库里面有这方面的文章，本人跟你的专业背景不一样，你自己看看。题目如上