How to solve system of second order differential equations

Hi there everyone.

I am trying to solve a system of second order differential equations for a mass spring damper as shown in the attached picture using ODE45. The data etc is below;

top mass (ms) = 100
lower mass (mu) = 18
spring rate sprung (Ks) = 20,000
spring rate unsprung (Ku) = 167,000
Damping ratio sprung (Cs) = 1200
Damping ratio unsprung (Cu) = 6
sprung displacement (xs) = unknown, variable trying to calculate
unsprung displacement (xu) = unknown, variable trying to calculate
road vertical displacement (y) = 1-cos^2(theta) bump, 260mm long, 50mm high at peak

I have been struggling to get any data other than a straight line when it should show something like in the graph in the second picture. I know that i'll have to modify the output to give me the positions of the sprung and unsprung masses (m1 & m2). However i don't understand why i'm getting these results. The only output i have managed to get is attached and looks like a typical damping graph for a motorcycle which is what i am modelling. However i'm looking to get positions.

I know i have to turn the equations into a first order system (for some weird reason when mathematica can do it as is?!) so following some of the posts on here i have done that, but what's confusing me most is that i have to define xu and xs, which is what i am trying to find? Basically i'm just trying to bodge it and could use some guidance and an explanation past the documentation as it from what i've found it is just talking about a system of equations to be solved, or solving a single second order differential, not a system of them.

I have put my code below which gives the damping type output, however i don't want to have to define xs and xu as that's what i'm trying to determine. the data in the mat file is also included as an attachment So any help would be greatly appreciated. Thanks, Ross

       function dydt = MSS(t,y)

      load('Standard_European_Bump_Data.mat','length','height')
      height1 = height/1000
      length1 = length/1000

      h_vs_l = plot(length1, height1)

      ms = 100
      mu = 18
      Ks = 20000
      Ku = 167000
      Cs = 1200
      Cu = 6

      y = height1

      xu = zeros(360,1)
      xs = zeros(360,1)

      %   y2' = (Ks/ms)*y1 + (Cs/ms)*y2 == 0
      %   y4' = (-Ks/ms)*y1 + (-Cs/ms)*y + (Ku/mu)*y3 + (Cu/mu)*y4 == 0

      y1 = (xs - xu)
      y2 = resample((diff(y1)),361,360) 
      y3 = (xu - y)
      y4 = resample((diff(y3)),361,360)

      dy1 = y2
      dy2 = y4
      dy3 = (Ks/ms)*y1 + (Cs/ms)*y2
      dy4 = (-Ks/ms)*y1 + (-Cs/ms)*y2 + (Ku/mu)*y3 + (Cu/mu)*y4

      dydt = [dy1;dy2;dy3;dy4]

      end

and the solving part of the code in a separate file

      y0 = zeros(1440,1)
      tspan = [0 10]

      [t,y]=ode45(@Quarter_Car, tspan, y0) 

      y1 = (1:1:77)

      % plot(t,y(:,1))
      % hold on
      % xlabel('time (s)')
      % ylabel('Displacement??')
      % axis([0,100,-120,120])

      plot(t,y1)

[Jan]

I've formatted your code. Please read http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup.

[Ross Hanna]

Hi there. My apologies for the bad formatting. I see in the code you have added another end and an if true.

May i ask why the second 'end', as i'm not sure what it is referring to? And also, why the 'if true' part? Are these just for formatting reasons? as if i put it into matlab it throws an error at them?

thanks Ross

[Jan]

@Ross: If you hit the "{} Code" button and no code was selected before, the forum's interface inserts "if true; end" for unknown reasons - a kind of default code. Either delete it or simply select the code before hitting the button.

[Ross Hanna]

ah i see, highlight the code, then press "{} Code". Thanks

Answer 1

[Jan]

Your function calculates:

xu = zeros(360,1)
xs = zeros(360,1)
y1 = (xs - xu)
y2 = resample((diff(y1)),361,360) 
y3 = (xu - y)
y4 = resample((diff(y3)),361,360)

If xu and xs contains zeros only, y1 and y2 must be zero also. The diff command is very strange here and it is meaningless to resample the data afterwards. I assume you try to create a system of 1st order is completely wrong. Please take a look into the example of the documentation, e.g. https://www.mathworks.com/help/matlab/ref/ode45.html#bu3uj8b.

This does not look meaningful also:

%   y2' = (Ks/ms)*y1 + (Cs/ms)*y2 == 0
%   y4' = (-Ks/ms)*y1 + (-Cs/ms)*y + (Ku/mu)*y3 + (Cu/mu)*y4 == 0

If y2' and y4' equals 0, you can omit the computations.

[Ross Hanna]

Hi Jan, thanks for responding.

the xu and xs being zeros was me trying to get the matrices the same as the "height" size. My issue here is that xu and xs are what i am trying to find, so i couldnt see why i had to define them. I have since found that by putting y(1) instead of y1 eliminates that but i cant find why?.

As for the second point of the equations equalling zero, that was me not having deleted the "== 0" after transposing it so y2' and y4' were the subject instead of all of the variables on one side equalling zero, such as in the equation

   F - Kx - Cv - Ma = 0

what i am trying to do is model the quarter car model of the mass spring dampers by finding the dynamic equations. Once i have done that i wanted to solve them in matlab, to get what the positions of xu and xs would be over the time period specified. Below are my original equations and what i have done to reduce them to first order.

    m2.a = Ks*(xs*xu) + Cs*(xs'-xu')
    m1.a = -Ks(xs-xu) + Ku(xu - y) - Cs*(xs' - xu') + Cu*(xu' - y')

therefore the equations can be written as;

    m2.a = Ks*y1 + Cs*y2
    m1.a = Ks*y1 + Ku*y3 - Cs*y2 + Cu*y4

where;

    y1 = (xs - xu)
    y2 = y1' = (xs' - xu')
    y3 = (xu - y)
    y4 = y3' = (xu' - y')

I get to this point and then i'm unsure how to go about it as there are two equations instead of one, so i have tried putting them in the brackets like in the equation. I have also changed y = height1 to Z = height1 in case Matlab was getting confused with the y's.

I am now getting graphs like i have attached, which is better than straight lines, but still nto correct. As i have y1 = (xs - xu), will i have to rearrange the equation to get the output like in the "Suspension response" picture i originally attached?

Below is my new code and i have attached a picture of the graphs i am getting out.

    function dydt = MSS(t,y)
    load('Standard_European_Bump_Data.mat','length','height')
    height1 = height/1000
    length1 = length/1000
    %h_vs_l = plot(length1, height1)
    ms = 100
    mu = 18
    Ks = 20000
    Ku = 167000
    Cs = 1200
    Cu = 6
    Z = height1
    % y(1) = (xs - xu)
    % y(3) = (xu - Z)
    % y1' = y(2)              % = (xs' - xu')
    % y2' = (Ks/ms)*y(1) + (Cs/ms)*y(2)
    % y3' = y(4)              % = (xu - Z)
    % y(4)' = -(Ks/ms)*y(1) - (Cs/ms) * y(2) + (Ku/mu)*y(3) + (Cu/mu) * y(4)
    dydt = [y(2); (Ks/ms)*y(1) + (Cs/ms)*y(2); y(4); -(Ks/ms)*y(1) - (Cs/ms) * y(2) + (Ku/mu)*y(3) + (Cu/mu) * y(4)]
    end

and what i have changed the ODE file to;

    y0 = [0 0.7 0 0]
    tspan = [0 10]
    [t,y]=ode45(@Quarter_Car_3, tspan, y0)
    hold on
    plot(t,y(:,1))
    plot(t,y(:,2))
    plot(t,y(:,3))
    plot(t,y(:,4))
    xlabel('time (s)')
    ylabel('Displacement??')
    axis([0,10,-120,120])