How to solve Multiple DOF Mass Spring Damper system and find/plot acceleration, velocity & displacement vs time

64 次查看（过去 30 天）
AJMIT KUMAR
AJMIT KUMAR 2021-1-24
回答： AJMIT KUMAR 2021-1-25
I know to solve 1 DOF system but i dont know to solve multiple DOF (Matrix will come instead of single value);
Also how to find acceleration?
I am getting absurd/wrong values. Kindly correct the code
Function Code :-
% ODE Function
function dxdt = Ransom(t, x, M, K, C)
m1 = 25.81;m2 = 15.36; m3 = 40.35; m4 = 10.12;
c1 = 3564; c2 = 5685; c3 = 8550; c4 = 620;
k1 = 950000; k2 = 262800; k3 = 213000; k4 = 450000;
w = 50; F = 5;
M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];
C = [c1+c2 -c2 0 0; -c2 c2+c3 -c3 0; 0 -c3 c3+c4 -c4; 0 0 -c4 c4];
K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];
dxdt = zeros(8, 1) ;
dxdt(1) = x(5) ;
dxdt(2) = x(6) ;
dxdt(3) = x(7) ;
dxdt(4) = x(8) ;
dxdt(5) = 5*sin(w*t) -K(1)/M(1) * x(1) - K(2)/M(1) * x(2) -C(1)/M(1) * x(5) - C(2)/M(1) * x(6) ;
dxdt(6) = -K(5)/M(6) * x(1) - K(6)/M(6) * x(2) - K(7)/M(6) * x(3) -C(5)/M(6) * x(5) - C(6)/M(6) * x(6) - C(7)/M(6) * x(7) ;
dxdt(7) = -K(10)/M(11) * x(1) - K(11)/M(11) * x(2) - K(12)/M(11) * x(3) -C(10)/M(11) * x(5) - C(11)/M(11) * x(6) - C(12)/M(11) * x(7) ;
dxdt(8) = -K(15)/M(16) * x(3) - K(16)/M(16) * x(4) -C(15)/M(16) * x(7) - C(16)/M(16) * x(8) ;
end
Calling Code: -
clear all;clc;
m1 = 25.81;m2 = 15.36; m3 = 40.35; m4 = 10.12;
c1 = 3564; c2 = 5685; c3 = 8550; c4 = 620;
k1 = 950000; k2 = 262800; k3 = 213000; k4 = 450000;
w = 50; F = 5;
M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];
C = [c1+c2 -c2 0 0; -c2 c2+c3 -c3 0; 0 -c3 c3+c4 -c4; 0 0 -c4 c4];
K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];
tspan = [0 1000] ;
y0 = [0 0 0 0 0 0 0 0] ;
[t,x]=ode45('Ransom',tspan,y0);
x_ = x(:, 1:4) ;
xdot_ = x(:, 5:8) ;
subplot(2,1,1)
plot(t, x_);
subplot(2,1,2)
plot(t, xdot_);

请先登录，再回答此问题。

采纳的回答

AJMIT KUMAR
AJMIT KUMAR 2021-1-25
Finally I did it myself
Function Code: -
% ODE Function
function dxdt = Ransom(t,x)
m1 = 25.81; m2 = 15.36; m3 = 40.35; m4 = 10.12; % Mass(kg)
c1 = 3564; c2 = 5685; c3 = 8550; c4 = 620; % Damping Coefficient
k1 = 950000; k2 = 262800; k3 = 213000; k4 = 450000; % Stiffness
w = 100; % Angular Velocity (rad/s)
F = 5;
M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];
C = [c1+c2 -c2 0 0; -c2 c2+c3 -c3 0; 0 -c3 c3+c4 -c4; 0 0 -c4 c4];
K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];
dxdt = zeros(8,1) ;
% x(1),x(2),x(3)and x(4) are displacement of mass m1, m2, m3 & m4 respectively
dxdt(1) = x(5) ; % velocity of mass m1
dxdt(2) = x(6) ; % velocity of mass m2
dxdt(3) = x(7) ; % velocity of mass m3
dxdt(4) = x(8) ; % velocity of mass m4
dxdt(5) = (F*sin(w*t))/M(1) -K(1)/M(1) * x(1) - K(2)/M(1) * x(2) -C(1)/M(1) * x(5) - C(2)/M(1) * x(6) ;
dxdt(6) = -K(5)/M(6) * x(1) - K(6)/M(6) * x(2) - K(7)/M(6) * x(3) -C(5)/M(6) * x(5) - C(6)/M(6) * x(6) - C(7)/M(6) * x(7) ;
dxdt(7) = -K(10)/M(11) * x(2) - K(11)/M(11) * x(3) - K(12)/M(11) * x(4) -C(10)/M(11) * x(6) - C(11)/M(11) * x(7) - C(12)/M(11) * x(8) ;
dxdt(8) = -K(15)/M(16) * x(3) - K(16)/M(16) * x(4) -C(15)/M(16) * x(7) - C(16)/M(16) * x(8) ;
end
Calling Code: -
clear all;clc;
% Displacement and Velocity Calculation
tspan = [0 1]; % Time (sec)
y0 = [1 1 1 1 0 0 0 0]; % Initial conditions
[t,x]=ode45('Ransom',tspan,y0);
x_ = x(:,1:4); % All Displacement
xdot_ = x(:,5:8); % All Velocity
% Plot Displacement wrt Time
subplot(3,1,1)
plot(t, x_ ,'color','b','LineWidth',2);
grid on
xlabel('Time (sec)')
ylabel('Displacement(m)')
title('Displacement Vs Time')
% Plot Velocity wrt Time
subplot(3,1,2)
plot(t, xdot_ ,'color','g','LineWidth',2);
grid on
xlabel('Time (sec)')
ylabel('Velocity (m/s)')
title('Velocity Vs Time')
% Acceleration Calculation
m1 = 25.81; m2 = 15.36; m3 = 40.35; m4 = 10.12; % Mass(kg)
c1 = 3564; c2 = 5685; c3 = 8550; c4 = 620; % Damping Coefficient
k1 = 950000; k2 = 262800; k3 = 213000; k4 = 450000; % Stiffness
w = 100; % Angular Velocity (rad/s)
F = 5;
M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];
C = [c1+c2 -c2 0 0; -c2 c2+c3 -c3 0; 0 -c3 c3+c4 -c4; 0 0 -c4 c4];
K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];
a1 = (F*sin(w*t))/M(1) -K(1)/M(1) * x(:,1) - K(2)/M(1) * x(:,2) -C(1)/M(1) * x(:,5) - C(2)/M(1) * x(:,6) ;
a2 = -K(5)/M(6) * x(:,1) - K(6)/M(6) * x(:,2) - K(7)/M(6) * x(:,3) -C(5)/M(6) * x(:,5) - C(6)/M(6) * x(:,6) - C(7)/M(6) * x(:,7) ;
a3 = -K(10)/M(11) * x(:,2) - K(11)/M(11) * x(:,3) - K(12)/M(11) * x(:,4) -C(10)/M(11) * x(:,6) - C(11)/M(11) * x(:,7) - C(12)/M(11) * x(:,8) ;
a4 = -K(15)/M(16) * x(:,3) - K(16)/M(16) * x(:,4) -C(15)/M(16) * x(:,7) - C(16)/M(16) * x(:,8) ;
a = [a1 a2 a3 a4];
% Plot Acceleration wrt Time
subplot(3,1,3)
plot(t, a ,'color','r','LineWidth',2);
grid on
xlabel('Time (sec)')
ylabel('Acceleration (m/s^2)')
title('Acceleration Vs Time')

更多回答（1 个）

Stephan
Stephan 2021-1-24
tspan=[0 10]; % Simulation time
z0=[0;0]; % Initial conditions for forced vibration
[t,z]=ode45(@forced_vibration,tspan,z0);
subplot(2,1,1)
plot(t,z(:,1),'color','b','LineWidth',2);
grid on
xlabel('Time (t) (sec)')
ylabel('Displacement (x) (m)')
title('Displacement Vs Time (Forced Vibration)')
subplot(2,1,2)
plot(t,z(:,2),'color','b','LineWidth',2);
grid on
xlabel('Time (t) (sec)')
ylabel('Velocity (xdot) (m/s)')
title('Velocity Vs Time (Forced Vibration)')
function f = forced_vibration(t,z)
m = 10; % Mass(kg)
c = 10; % Damping coefficient (Ns/m)
k = 1000; % Stiffness coefficient (N/m)
F = 1; % Force (N)
w = 50; % Excitation frequency(rad/s)
f = [z(2);((F/m)*sin(w*t)-(c/m)*z(2)-(k/m)*z(1))];
end

类别

在 Help Center 和 File Exchange 中查找有关 Acoustics, Noise and Vibration 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by