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

97 次查看（过去 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 中查找有关 Programming 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by