Spring Mass system (displacement)

Question

0 投票

I have acceleration data, m,c,k and how to write ode45 to find displacement? It is a 3DOF system

The below is my matlab code

Mx"+cx'+kx=0

1 件のコメント
-1 件の古いコメントを表示 -1 件の古いコメントを非表示

Jason Hui 2018 年 3 月 7 日

編集済み: KSSV 2018 年 3 月 7 日

MATLAB Online で開く

function [x,v]=dynresp(t, x0,v0, M, C, K, F, fs, NMODE, DOF_Output)
% [t,x,v]=strdynrk(t,x0,v0,m,c,k,functim)
% This function uses ode45 to solve the differential equation
% modal superposition is applied to system equations: M*X"+C*X'+K*X=F(t)
% the first NMODE modes are used only
% t       - row vector of solution times (1*NTIME)
% x0,v0   - initial position and velocity vectors (NDOF*NTIME)
% M,C,K   - mass, damping and stiffness matrices (NDOF*NDOF)
% F       - driving force vector (NDOF*NTIME)
% fs      - sampling frequency
% x,v     - arrays containing solution values for position and velocity (NDOF*NTIME)
%DOF_Output: if available, only x and v at this point are output.
%
NDOF=length(M); 
% eigen-analysis
DISPY=struct('disp',0);
[MODE,EIGVALUE]=eigs(K,M,NMODE,0,DISPY);
EIGVALUE=diag(EIGVALUE);
[EIGVALUE,ORDER]=sort(EIGVALUE); MODE=MODE(:,ORDER);
MODE=mod_norm_mass(MODE, M);
% transform
F=MODE'*F; F=[F  zeros(size(F,1),1)]; 
x0=MODE\x0;
v0=MODE\v0;
%M=eye(NMODE);
%C=MODE'*C*MODE;
%K=diag(EIGVALUE);
M=ones(NMODE,1); C=diag(MODE'*C*MODE); K=EIGVALUE;
[t,z]=ode45(@fun_intg,t,[x0; v0],[], M,C,K,F, fs);
if nargin==10
    x=MODE(DOF_Output,:)*z(:,1:NMODE)'; v=MODE(DOF_Output,:)*z(:,NMODE+1:2*NMODE)';
else
    x=MODE*z(:,1:NMODE)'; v=MODE*z(:,NMODE+1:2*NMODE)';
end
return
%================================
function out=fun_intg(t,z, M,C,K,F, fs)
NDOF=length(M); 
%TIME=round(t*fs); f=F(:,TIME+1);
TIME=t*fs; TIME1=floor(TIME); TIME2=TIME1+1; f=F(:,TIME1+1)*(TIME2-TIME) + F(:,TIME2+1)*(TIME-TIME1); % interpolation
out=[z(NDOF+1:2*NDOF); (f - C.*z(NDOF+1:2*NDOF) - K.*z(1:NDOF))./M];
return
%================================
% A typical call to strdynrk function is:
% trial for integration
% for a 3dof system, it takes 7 seconds only to get responses of
% 1000 seconds with frequency of 100 Hz.
m=eye(3,3); k=[2,-1,0;-1,2,-1;0,-1,2];
c=.05*k; x0=zeros(3,1); v0=zeros(3,1);
% t=linspace(0,10,101);
% F=[0:10 10:-1:0 zeros(1,79);zeros(2,101)];
t=linspace(0,1000,100001);
F=[rand(1,length(t));zeros(2,length(t))];
[x,v]=dynresp(t,x0,v0,m,c,k,F,10,3);
plot(t,x), hold on, plot(t,F'), legend('1','2','3','F')