# Undesired Graph for Transmissibility vs Frequency.

7 views (last 30 days)
Ankit Joshi on 29 Sep 2021
Greetings All,
I've been wokring on a 4DOF model problem and I am not able to get the desired result for Transmissibility graph of the system. As it is a 4DOF system, the graph should have 4 peaks, but I'm only getting one peak.
The transmissibility if of the 4th system.
Following is the code:
clear all;clc;
% Displacement and Velocity Calculation
tspan = [0 4]; % Time (sec)
ic = [0.04 -0.02 0.05 -0.04 0 0 0 0]; % Initial conditions
[t,x]=ode45(@fourforceddamped,tspan,ic);
x_ = x(:,1:4); % All Displacement
xdot_ = x(:,5:8); % All Velocity
% Acceleration Calculation
m1 = 95; m2 = 945; m3 = 10; m4 = 80; % Mass(kg)
c1 = 800; c2 = 2000; c3 = 1200; c4 = 1200; % Damping Coefficient
k1 = 660000; k2 = 86000; k3 = 1000; k4 = 10000; % Stiffness
F = (m1+m2+m3+m4)*30; %Input force
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];
[e,d]=eig(K,M);
w=sqrt(abs(d));
a1 = (F)/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) ;
z=c4/(2*sqrt(m4*k4));
om=0:1000;
fr1=om/w(1);
fr2=om/w(6);
fr3=om/w(11);
fr4=om/w(16);
T1=sqrt((1+(2.*fr1.*z).^2)./((1-fr1.^2).^2+(2.*fr1.*z).^2));
T2=sqrt((1+(2.*fr2.*z).^2)./((1-fr2.^2).^2+(2.*fr2.*z).^2));
T3=sqrt((1+(2.*fr3.*z).^2)./((1-fr3.^2).^2+(2.*fr3.*z).^2));
T4=sqrt((1+(2.*fr4.*z).^2)./((1-fr4.^2).^2+(2.*fr4.*z).^2));
plot(om,T1);
plot(om,T2);
plot(om,T3);
plot(om,T4);
hold on;
grid on;
% ODE Function
function dxdt = fourforceddamped(t,x)
m1 = 95; m2 = 945; m3 = 10; m4 = 80; % Mass(kg)
c1 = 800; c2 = 2000; c3 = 1200; c4 = 1200; % Damping Coefficient
k1 = 660000; k2 = 86000; k3 = 1000; k4 = 10000; % Stiffness
F = (m1+m2+m3+m4)*30;
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];
[e,d]=eig(K,M);
w=sqrt(abs(d));
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)/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

Sulaymon Eshkabilov on 30 Sep 2021
This is how you should plot it in this way:
...
T4=sqrt((1+(2.*fr4.*z).^2)./((1-fr4.^2).^2+(2.*fr4.*z).^2));
loglog(om, T1+T2+T3+T4, 'b', 'linewidth', 2) % This is logarithmic scaled plot
grid on
xlabel('Frequency, Hz')
ylabel('|T|')
...
##### 2 CommentsShowHide 1 older comment
Sulaymon Eshkabilov on 30 Sep 2021

Sulaymon Eshkabilov on 29 Sep 2021
You'd need to plot in this way:
plot(om,T1); hold on
plot(om,T2);
plot(om,T3);
plot(om,T4);
grid on;
hold off
OR
plot(om, [T1; T2; T3; T4])
##### 2 CommentsShowHide 1 older comment
Ankit Joshi on 29 Sep 2021
Something Like this