How quarter car modeling on MATLAB?

I am trying to model the quarter suspension model in Inman's engineering vibration book in MATLAB, but I have not made any progress. The vehicle will move on the sinusoidal path as in the figure and the height of the road is 5 cm.

I want to plot the largest vibration amplitude of the vehicle for the speed range v=[0 200] km/h and the total vibration for r=0.3 for the time interval t=[0 20] s. The average vehicle weight will be M=1000 kg. I predict the spring stiffness to be K=140 kN/m, damping coefficient c=900 Ns/m.

[Image Analyst]

Do you have any formulas at all? The Crystal Ball Toolbox is still under development so I cannot currently see your engineering vibration book.

[Mustafa Furkan]

@Image Analyst Although the values are different, if I want to solve this problem analytically, it will be a solution like the image above.

Answer 1

[Sam Chak]

Hi @Mustafa Furkan,

I prefer the simulation-based numerical solution over the formula-based analytical solution in the example because there is no graph provided in your image.

By Newton's 2nd law, we should get

where the disturbance force caused by the wavy road surface is given by

.

Thus, the model can be rewritten as

,

where the road surface is given by . I believe you can use calculus to find .

The following is based on the info in Example 2.4.2 (provided by you).

tspan  = linspace(0, 20, 20001);
x0     = [0 0];                     % initial rest condition
[t, x] = ode45(@quarterCar, tspan, x0);
plot(t, x(:,1)), grid on, xlabel('Time, [seconds]'), ylabel('Deflection, [meter]')
title('Deflection experienced by the car')

% Maximum deflection experienced by the car (same as in Example 2.4.2)
dmax = max(x(:,1)) - min(x(:,1))

dmax = 0.0320

% Quarter Car suspension model
function xdot = quarterCar(t, x)
    xdot = zeros(2, 1);	
    % Example 2.4.2
    m    = 1007;                % mass of car
    c    = 2000;                % damping coefficient
    k    = 4e4;                 % stiffness coefficient
    h    = 1;                   % road height in cm
    v    = 20;                  % car velocity (fixed speed)
    
    % Original Problem
%     m    = 1000;                % mass of car
%     c    = 900;                 % damping coefficient 
%     k    = 140e3;               % stiffness coefficient 
%     h    = 5;                   % road height in cm
%     v    = (200/20)*t;          % car velocity ramp up from 0 to 200 km/h in 20 seconds

    % Standard formulas and state-space model (1st-order ODE form)
    wb   = 0.2909*v;            % formula given in Example 2.4.2
    a    = (h/2)/100;           % sinusoidal amplitude in m
    d    = k*a*sin(wb*t) + c*a*wb*cos(wb*t);    % disturbance due to wavy road surface
    A    = [0 1; -k/m -c/m];    % state matrix
    B    = [0; 1/m];            % input matrix
    xdot = A*x + B*d;
end

[Mustafa Furkan]

@Sam Chak thank you your code really helped me. If I want to add the following value to the code and have an r-X plot drawn, how can I do that?

r = wb/6.303; %frequency ratio

[Sam Chak]

HI @Mustafa Furkan

The r-X plot is not shown in your image, so I don't know how it looks like.

Can you type out the formulas {, ζ, X} as shown in the image (Example 2.4.2) so that they are easy to be copied and tested later?

[Mustafa Furkan]

@Sam Chak hi,

You can see an example of the r-x plot above.The formula required to obtain the r values here is used at the bottom of the image of example 2.4.2. The formula is as follows:

r=wb/wn (frequency ratio)

[Sam Chak]

Yes I see it. Can you paste the same equation/formula here in the comment section?

It will be convenient for copying your original error-free equation.

[Mustafa Furkan]

r=wb/wn; % this is the basic formula of frequency response

% wn=6.303 (constant)

% wb=0.2909*v (from previous code written by you)

[Sam Chak]

I'm outside on this festive day. If I return, I may show the plot properly. Here is how you plot it:

r = linspace(0, 3, 3001);
X = 0.01*sqrt((1 + (2*0.158*r).^2)./((1 - r.^2).^2 + (2*0.158*r).^2));
plot(r, X), grid on
xlabel('r'), ylabel('X')

[Mustafa Furkan]

Thank you @Sam Chak, waiting for your return. I tried something different from the code you wrote, but like the one below, but the result was not what I wanted.

t=0:0.01:20;
v=(200/20)*t;
wb=0.2909*v;
r = wb/6.303;
ksi= 0.038;
X = 0.025*sqrt((1 + (2*ksi*r).^2)./((1 - r.^2).^2 + (2*ksi*r).^2));
plot(r, X)

[Sam Chak]

Hi @Mustafa Furkan

Unless you implied the author, Prof. Daniel J. Inman made an error in his book, I don't see any problem with the original formula that you supplied (in the image). I have checked the formula a few times. THe code looks the same as in the image.

Perhaps, you can be a little specific about the result or plot that you are expecting.

% t   = 0:0.01:20;
% v   = (200/20)*t;
% wb  = 0.2909*v;
% r   = wb/6.303;
r   = 0:0.01:10;
ksi = 0.038;
X   = 0.025*sqrt((1 + (2*ksi*r).^2)./((1 - r.^2).^2 + (2*ksi*r).^2));
plot(r, X), grid on
xlabel('r'), ylabel('X')

[KARTHIK]

when i try this code, automatically more that 20 plots are genreating. because of that I am getting pop-up window, which I could not able to close. please guide.

For ref;

tspan  = linspace(0, 20, 20001);
x0     = [0 0];                     % initial rest condition
[t, x] = ode45(@quarterCar, tspan, x0);
plot(t, x(:,1)), grid on, xlabel('Time, [seconds]'), ylabel('Deflection, [meter]')
title('Deflection experienced by the car')

% Maximum deflection experienced by the car (same as in Example 2.4.2)
dmax = max(x(:,1)) - min(x(:,1))

% Quarter Car suspension model
function xdot = quarterCar(t, x)
    xdot = zeros(2, 1);	
    % Example 2.4.2
    m    = 1007;                % mass of car
    c    = 2000;                % damping coefficient
    k    = 4e4;                 % stiffness coefficient
    h    = 1;                   % road height in cm
    v    = 20;                  % car velocity (fixed speed)
    
    % Original Problem
    %     m    = 1000;                % mass of car
    %     c    = 900;                 % damping coefficient 
    %     k    = 140e3;               % stiffness coefficient 
    %     h    = 5;                   % road height in cm
    %     v    = (200/20)*t;          % car velocity ramp up from 0 to 200 km/h in 20 seconds
    
    % Standard formulas and state-space model (1st-order ODE form)
    wb   = 0.2909*v;            % formula given in Example 2.4.2
    a    = (h/2)/100;           % sinusoidal amplitude in m
    d    = k*a*sin(wb*t) + c*a*wb*cos(wb*t);    % disturbance due to wavy road surface
    A    = [0 1; -k/m -c/m];    % state matrix
    B    = [0; 1/m];            % input matrix
    xdot = A*x + B*d;
    
    %% ----- I inserted these inside the function quarterCar() -----
    r   = 0:0.01:10;
    ksi = 0.038;
    X   = 0.025*sqrt((1 + (2*ksi*r).^2)./((1 - r.^2).^2 + (2*ksi*r).^2));
    plot(r, X), grid on
    xlabel('r'), ylabel('X')
    %% -------------------------------------------------------------
end

[Sam Chak]

Hi @KARTHIK

Your code generated multiple plots because you incorrectly inserted the plotting lines into the quarterCar() function and then called the ode45() solver to perform internal integration in the quarterCar(). The following code fixes your issue. Follow the steps as guided.

%% Step 2 of 5: Call ode45 to solve quarterCar()
tspan  = linspace(0, 20, 20001);
x0     = [0 0];                     % initial rest condition
[t, x] = ode45(@quarterCar, tspan, x0);

%% Step 3 of 5: Plot deflection in Figure 1
figure(1)
plot(t, x(:,1)), grid on, xlabel('Time, [seconds]'), ylabel('Deflection, [meter]')
title('Deflection experienced by the car')

%% Step 4 of 5: Find the Maximum deflection experienced by the car
dmax = max(x(:,1)) - min(x(:,1))

dmax = 0.0320

%% Step 5 of 5: Plot the Frequency response in Figure 2
r   = 0:0.01:10;
ksi = 0.038;
X   = 0.025*sqrt((1 + (2*ksi*r).^2)./((1 - r.^2).^2 + (2*ksi*r).^2));
figure(2)
plot(r, X), grid on
xlabel('r'), ylabel('X')
title('Frequency response by the car')

%% Step 1 of 4: Write the code for the Quarter Car suspension model
function xdot = quarterCar(t, x)
    xdot = zeros(2, 1);	
    % Example 2.4.2
    m    = 1007;                % mass of car
    c    = 2000;                % damping coefficient
    k    = 4e4;                 % stiffness coefficient
    h    = 1;                   % road height in cm
    v    = 20;                  % car velocity (fixed speed)
    
    % Original Problem
    %     m    = 1000;                % mass of car
    %     c    = 900;                 % damping coefficient 
    %     k    = 140e3;               % stiffness coefficient 
    %     h    = 5;                   % road height in cm
    %     v    = (200/20)*t;          % car velocity ramp up from 0 to 200 km/h in 20 seconds
    
    % Standard formulas and state-space model (1st-order ODE form)
    wb   = 0.2909*v;            % formula given in Example 2.4.2
    a    = (h/2)/100;           % sinusoidal amplitude in m
    d    = k*a*sin(wb*t) + c*a*wb*cos(wb*t);    % disturbance due to wavy road surface
    A    = [0 1; -k/m -c/m];    % state matrix
    B    = [0; 1/m];            % input matrix
    xdot = A*x + B*d;
end

[KARTHIK]

Thank you @Sam Chak. Code is working fine

[KARTHIK]

@Sam Chak Sir, Please guide me for this problem by ode45

I tried with above provided code, I couldnot able to feed this function.

I need to plot amplitude vs velocity

Answer 2

[Sulaymon Eshkabilov]

Hi,

Your system equation is: M*ddx+C*dx+K*x = F(t) and you are studying a forced vibration.

There are a few different ways how to model and solve this spring-mass-damper system.

(1) To obtain an analytical solution, use a symbolic math function: dsolve()

(2) To obtain an analytical solution, use symbolic math functions: laplace(), ilaplace()

(3) To obtain a numerical solution, use ode solvers: ode45, ode23, ode113, etc

(4) To obtain a numerical solution, use Simulink modelling: ode45, ode23, ode113, etc

(5) To obtain a numerical solution, write a code using Euler, Runge-Kutta and other methods

(6) To obtain a numerical solution, use control toolbox functions: tf(), lsim()

There are a few nice sources how to model and solve this spring-mass-damper system.

(4) https://www.mathworks.com/matlabcentral/fileexchange/94585-mass-spring-damper-systems

(6) https://www.mathworks.com/matlabcentral/fileexchange/95288-mass-spring-damper-1-dof

(6) https://www.instructables.com/How-to-Model-a-Simple-Spring-mass-damper-Dynamic-S/

Nice sim demo: https://www.mathworks.com/matlabcentral/fileexchange/94585-mass-spring-damper-systems

[Mustafa Furkan]

Thank you. I've reviewed all of your examples before. I'm new to MATLAB so I'm having trouble implementing it. I don't know how to apply range values or draw related graphs.