Mike Deffenbaugh

Graphing MDOF Systems


In the Vibrations section of my website I have a graph of a MDOF (Multiple Degree of Freedom) system's frequency response plot and an explanation of how to graph a SDOF in MATLAB. Well someone recently sent me an email asking if they could help them graph a MDOF system. Here is the one on SDOF and below I will explain how to graph a MDOF system.

The system that I will be using for this example will be the serpentine belt from an old Ford engine. Solving this system was a part of one of my assignments in Vibrations class. I won't go into the details of this specific system, I'm assuming that you already have your system set up (skip ahead if you do). If you don't have an example of a MDOF system then you can use my code to set up a system. If you don't understand how a MDOF system works you should read more about MDOF systems first because I will not be explaining their basic concepts in this guide.

This system can be defined using the following script file and function file. The script file begins with the following system definitions.


%7DOF pulley system
%input parameters
mass = .2663; %kg, mass of pulley on tensioner spring
J = [0 41.5 13.1 2.63 42.1 17.6 2.08]/10000;  %kg m^2, rotational inertias
r = [81.25 64.5 70.60 -41.15 30 67.5 -38.1]/1000; %m, radii of pulleys
x = [0 261.5 252 90.3 86 0 128]/1000; %m, x location of pulleys
y = [0 60 234 251.1 354 167.5 126]/1000; %m, y location of pulleys

EA = 65000; %N
kT = 9340; %N/m, stiffness of tensioner spring, assumed to be linear

%frequency parameters for graphing
Ldiv = 10; %increment for omega in rpm
Lom = 3000; %max omega in rpm

Then all of this information gets put into the MDOFpulley() function like this:


%response plot
plotOn = 1; %shows plot
[wnrpm,V,Amat,list] = MDOFpulley(J,r,x,y,mass,EA,kT,Ldiv,Lom,plotOn);

This function takes the system information and calculates the lengths of the pulley sections and then creates the equivalent mass, equivalent damping, and equivalent stiffness matrices. 


function [wnrpm,V,Amat,list] = MDOFpulley(J,r,x,y,mass,EA,kT,Ldiv,Lom,ploton)
%list = {[lo1,lo2],[EA,kT],J,r,x,y} this gives all the important parameters
%in a large list

    %solve for Lengths
    L = [0,0,0,0,0,0,0];

    n = 1;
    while n <= 7
        if n==7;
            L(n) = ((x(1)-x(n))^2 + (y(1)-y(n))^2 - (r(1)-r(n))^2)^(1/2);
            n = 8;
        else
            L(n) = ((x(n+1)-x(n))^2 + (y(n+1)-y(n))^2 - (r(n+1)-r(n))^2)^(1/2);
            n = n+1;
        end
    end

    k = EA./L;

    %create stiffness matrix
    keq1 = [(k(1)+k(2))*r(2)^2, -k(2)*r(2)*r(3), 0,0,0,0,0];
    keq2 = [-k(2)*r(2)*r(3), (k(2)+k(3))*r(3)^2, -k(3)*r(3)*(-r(4)),0,0,0,0];
    keq3 = [0, -k(3)*r(3)*(-r(4)), (k(3)+k(4))*(-r(4))^2, -k(4)*(-r(4))*r(5),0,0,0];
    keq4 = [0,0,-k(4)*(-r(4))*r(5), (k(4)+k(5))*r(5)^2, -k(5)*r(5)*r(6),0,0];
    keq5 = [0,0,0, -k(5)*r(5)*r(6),(k(5)+k(6))*r(6)^2, -k(6)*r(6)*(-r(7)),...
		(k(6)*r(6))*(x(6)-x(7))/L(6)];
    keq6 = [0,0,0,0,-k(6)*r(6)*(-r(7)),(k(6)+k(7))*(-r(7))^2 ,...
		-(k(6)*(-r(7))*(x(6)-x(7))/L(6))-(k(7)*(-r(7))*(x(7)-x(1))/L(7))];
    keq7 = [0,0,0,0, k(6)*(x(6)-x(7))*r(6)/L(6), ...
		-(k(7)*(x(7)-x(1))*(-r(7))/L(7)) - (k(6)*(x(6)-x(7))*(-r(7))/L(6)),...
		(k(7)*(x(7)-x(1))*(x(7)-x(1))/(L(7)^2)) + ((k(6)*(x(6)-x(7))^2)/(L(6)^2)) + kT];

    keq = [keq1;keq2;keq3;keq4;keq5;keq6;keq7];

    %adjust intertias, J = (mr^2)/2
    rorig = [81.25, 64.5, 70.60, -41.15, 30, 67.5, -38.1]/1000;
    for i = 1:6
        J(i) = (J(i)/((rorig(i)^4)))*(r(i)^4);
    end
    
    %create mass matrix
    meq = [];
    for i = 1:6
        meq(i,i) = J(i+1);
    end
    meq(7,7) = mass;

    %create damping "matrix", no damping for this problem
    ceq = 0;

Now the function takes these mass, damping, and stiffness matrices and calculates the natural frequencies and mode shapes of the systems. Most of the calculations are done with MATLABs built-in polyeig function. The rest of the code just extracts the important information. 


    [Q,lambda] = polyeig(keq,ceq,meq);
    %natural frequency = imaginary component
    lambda2 = imag(lambda);

    %extract only positive lambdas
    lambda3 = [];
    n = 1;
    %zeros = 0
    for i = 1:length(lambda2)
        if lambda2(i) > 0
            lambda3(n) = lambda2(i);
            n = n+1;

        end
    end
    %natural freq in RPM
    wnrpm = lambda3./(2*pi)*60*(1/3);

    %extract mode shapes from eigenvectors
    Q2 = [];
    n = 1;
    for i = 1:(length(Q)/2)
        Q2(:,n) = Q(:,2*i-1);
        n = n+1;
    end
    Q3 = abs(Q2);
    V = Q3;

    for i = 1:length(Q3)
        maxQ3 = max(Q3(:,i));
        Q3(:,i) = Q3(:,i)./maxQ3;
    end
    %mode shapes for each pulley at each frequency
    V = Q3;

This is the important part if you already have your own system. First we make the forces matrix that just states which pulleys have a force applied to them and what are the strength of each force relative to the others. In my case I just said that pulley 1, 6, and 7 all had an equal force of magnitude 1. Then make two empty matrices that will store all of the magnitude responses for each frequency that you want to graph. Then run it through a for-loop that finds these response magnitudes. 


    %MDOF frequency response
    forces = [1;0;0;0;0;1;1];
    Amat = zeros(7,Lom/Ldiv+1);
    omegavector = zeros(1,Lom/Ldiv+1);

    n = 0;
    for omega = 0:Ldiv:Lom
        A = ((-(omega/60*(2*pi)*3).^2.*(meq) + (sqrt(-1))*(omega/60*(2*pi)*3).*(ceq) + (keq))^-1)*(forces);
        n = n+1;
        omegavector(n) = omega;
        A = abs(A);
        Amat(:,n) = A;
    end

My test frequencies (omega) happen to be RPM values, but if you were using normal rad/s instead then your magnitudes should be written like this:


%if the input frequencies were in rad/s...
for rads = 0:Ldiv:Lom
    A = (  (  -rads.^2.*(meq) + (sqrt(-1))*(rads).*(ceq) + (keq)  )^-1  )*(forces);

Now we just need to plot our magnitudes stored in 'Amat' vs. our test frequencies stored in 'omegavector'


    figure()
    clf
    hold on
    axis([0, Lom, 0, .025])
    if ploton >= 1
        plot(omegavector,Amat(1,:),'r','Linewidth',2)
        plot(omegavector,Amat(2,:),'g','Linewidth',2)
        plot(omegavector,Amat(3,:),'b','Linewidth',2)
        plot(omegavector,Amat(4,:),'c','Linewidth',2)
        plot(omegavector,Amat(5,:),'y','Linewidth',2)
        plot(omegavector,Amat(6,:),'k','Linewidth',2)
        plot(omegavector,Amat(7,:),'m','Linewidth',2)
    end

Here is the final graph zoomed in on the first two natural frequencies. Also, remember that we used abs() to calculate the magnitude of our frequency response. Besides magnitudes there are also phase shifts to take into consideration. I left damping in this system as ceq = 0 but we could introduce damping into the system. If we did this the problem becomes more complicated, so for now I'm going to leave out a discussion of damping.