khyiu0415 发表于 2013-1-26 18:44

finite element method using matlab 錯誤

本帖最后由 khyiu0415 于 2013-1-26 19:46 编辑

請問大家有看過這本書嗎,當中的m-file'fediresp'是否有錯誤嗎?請指教

function =fediresp(M,K,F,u,t,C,q0,dq0,a,b);

%----------------------------------------------------------------------------------
%Purpose:
%   The function subroutine fediresp.m calulates impulse response
%   for a damped structural system using modal analysis. It uses modal
%   coordinate equations to evaluate modal responses anaytically, then
%   convert modal coordinates into physical responses
%
%Synopsis:
%   =fediresp(M,K,F,u,t,C,q0,dq0,a,b)
%
%Variable Description:
%   Input parameters - M, K - Mass and stiffness matrices
%                        F - Input or forcing influence matrix
%                        t - Time of evaluation
%                        u - Index for excitation
%                        C - Output matrix
%                        q0, dq0 - Initial conditions
%                        a, b - Parameters for proportional damping =a+b
%      Outpur parameters - eta - modal coordinate response
%                        y - physical coordinate response
%---------------------------------------------------------------------------------

disp('')
disp('Please wait!! - The job is being performed.')
%---------------------------------------------------------------
%Solve the eigenvalue problem and normalize the eigenvectors
% --------------------------------------------------------------

=size(M);=size(F);
nstep=size(t');

=eig(K,M);

=sort(diag(D));% Sort the eigenvaules and eigenvectors

V=V(:,k);

Factor=diag(V'*M*V);

Vnorm=V*inv(sqrt(diag(Factor)));      %Eigenvectors are normailzed

omega=diag(sqrt(Vnorm'*K*Vnorm));% Natural frequencies
                  
Fnorm=Vnorm'*F;

%----------------------------------------------------------------------------
%Compute modal damping matrix from the proportional damping matrix
%----------------------------------------------------------------------------
Modamp=Vnorm'*(a*M+b*K)*Vnorm;
zeta=diag((1/2)*Modamp*inv(diag(omega)))

if (max(zeta) >= 1),
disp('Warning - Your maximu damping ratio is grater than or equal to 1')
disp('You have to reselect a and b ')
else
pause
disp('If you want to continue, type return key')
end
%----------------------------------------------------------------------
%Find out impulse response of each modal coordinate analytically
%-------------------------------------------------------------------

eta0=Vnorm'*M*q0;            % Initial conditions for modal coordinates
deta0=Vnorm'*M*dq0;          % - both displacement and velocity
eta=zeros(nstep,n);

for i=1:n               %Responses are obtained for n modes
omegad=omega(i)*sqrt(1-zeta(i)^2);
phase=omegad*t;
tcons=zeta(i)*omega(i)*t;
eta(:,i)=exp(-tcons)'.*(eta0(i)*(cos(phase')+zeta(i)/sqrt(1-zeta(i)^2)*...
sin(phase'))+deta0(i)*sin(phase')/omegad+sin(phase')*Fnorm(i,u)/omegad);
end



%-----------------------------------------------------------------------
%Convert modal coordinate responses to physical coordinate responses
%-----------------------------------------------------------------------

yim=C*Vnorm*eta';
%-----------------------------------------------------------------------

tenglang 发表于 2013-1-28 12:55

具体错误请给出。
页: [1]
查看完整版本: finite element method using matlab 錯誤