Quater car modelling using Matlab

                                                      Quarter Car model 

Spring - Mass -Damper Dynamic System

Equations of motion of the mass damper system are

`mddotx+ c dotx+kx = 0 or F`

`omega_n = sqrt(k/m)``implies` Natural frequency

`C__critical = 2*sqrt(k*m)`

`zeta = ((C)/(C__critical))` `implies` = Damping ratio

 Quater Car Model: Equations of motion

`m_s` = sprung __mass

`m_(us)` = unsprung_mass

`k_s` = spring stiffness

`k_tr` = tyre stiffness

`C_sh` = Shock absorber damping coefficient

`C_t` = tyre damping coefficient

`Z_1` = Sprung mass displacement

`Z_2` = Unsprung mass displacement

`dotZ` = velocity

`ddotZ` = acceleration

Sprung mass equation of motion:

`m_s ddotZ_1 + C_sh(dotZ_1-dotZ_2) + k_s(dotZ_1 - dotZ_2) = 0`

Unsprung mass equation of motion

`m_(us) ddotZ_2 + C_sh(dotZ_2-dotZ_1) + k_s(dotZ_2 - dotZ_1) +C_t(dotZ_2-dotZ_0) + k_t(dotZ_2 - dotZ_0)  = 0`

Force on tyre

F(t) = `k_t*Z_0 + C_t*dotZ_0`

Free Responce / Undamped vibrations

The general solutions of differential equations are in the form of 

`Ae^(at)*Sin(omega*t)+Be^(at)*Sin(omega*t)`

`z1 = Z1sin(omega*t)`

`z2 = Z2sin(omega*t)`

`omega = omega_n`; For free responce natural frequency is taken into consideration

Damping terms are neglected

`C_sh` =0 ;  `C_t` =0 ;  `Z_1` =0 ; 

Substituting into sprung mass equation

`m_s d^2/dt^2(z_1sinomegat) + k_sz_1sinomegat-k_sz_2sinomegat =0`

Substituting into un_sprung mass equation

`m_(us) d^2/dt^2(z_2sinomegat) + k_sz_2sinomegat-k_sz_1sinomegat + k_tz_2sinomegat=0`

After simplification

`(-m_somega_n^2+k_s)z_1 - k_sz_2 =0`

`-k_sz1 + (-m_(us) omega_n^2 + k_s + k_t)z_2 =0`

writing the above equations in matrix form

`[[-m_somega_n^2+k_s,k_s],[-k_s,-m_(us) omega_n^2 + k_s + k_t]]((z1),(z2)) = 0`

If the determinant of the above matrix is zero; `abs(det)` = 0 then the equations are satisfied for any `z_1 and z_2` (characteristic equation)

`omega_n^4(m_sm_(us)) + omega_n^2(-m_sk_s-m_sk_t-m_(us)k_s) + k_sk_t = 0`

The above equation can be represented as `Ax^2+Bx+c`

`omega_n^2 to x`;

where  roots are `-B(+)/(-)sqrt(B^2-4AC)/(2A)`

`f_n = omega_n/(2*pi)`

`f_(n-s) = 1/2pi * sqrt((((kskt)/(ks+kt))/m_s)` `to` springs in series

`f_(n-us) = 1/2pi*sqrt((k_s+k_t)/m_(us))``to` springs in parallel

`star`Generally, unsprung mass is taken as 10% as that of the sprung mass.

 if `m_s` = 1800 kg ; `m_us` = 180 kg

    `k_s` = 88kN/m; `k_t` = 700 kN/m

then Natural frequency of sprung mass `f_(n-s)` = 1.04 Hz `to` body

       Natural frequency of unsprung mass `f_(n-us)` = 10 hz `to` wheel.

 Response to Road Excitation - Transmissibility Ratio

For external sources `omega_n` is replaced by `omega`

`[[-m_somega_n^2+c_somega_n+k_s,-c_somega_n-k_s],[-c_somega_n-k_s,-m_(us) omega_n^2 +(c_s+c_(us))omega_n+ k_s + k_(us)]]((z_s),(z_(us))) = ((0),(c_somega_n + k_(us))) Z_0`

A X = B

X = AinverseB (`A^(-1)B`)

`Z_s/Z_0  and Z_(us)/Z_0` are called Transmissibility ratio

 `Z_s/Z_0  = ((c_sc_(us)omega^2 + (k_sc_(us) + k_(us)c_s)omega + k_sk_(us))/( m_sm_(us)omega^4 + (m_sc_s+m_sc_(us)+m_u)omega^3+(k_sm_(us)+m_s(k_s+k_(us) +c_s(omega))omega^2)+(k_sc_(us)+k_(us)c_s)omega + k_sk_us))`

`Z_(us)/Z_0 = ((m_sc_(us)omega^3 + (c_sc_(us)+m_sk_(us))omega^2 + (k_sc_us+ k_w(s)omega)+k_sk_(us))/(( m_sm_(us)omega^4 + (m_sc_s+m_sc_(us)+m_u)omega^3+(k_sm_(us)+m_s(k_s+k_(us) +c_s(omega))omega^2)+(k_sc_(us)+k_(us)c_s)omega + k_sk_us))`

With regards to suspension performance, there are three aspects that need to be considered.

Suspension Performance Evaluation:

Vibration Isolation: 

This can be evaluated by the response of sprung mass (passenger area) to the excitation from the ground. 

Transmissibility ratio (`(Z_s/Z_0)`) can be used to assess vibration isolation.

Suspension travel:

It can be measured by spring deflection or in other terms it can be defined as the relative displacement between the sprung mass and the unsprung mass relative to road input         

                                       `(Z_(us) - Z_s)/(Z_0)`

It defines the space required to accommodate the suspension movement during jounce and rebound. Also known as 'Rattle space' 

Road holding:

It is necessary for the road to remain in contact with the ground. This is normally influenced by the force between the tyre and the road. This force can be represented by the dynamic tyre deflection relative to the road surface.

                                                `(Z_(us)-Z_o)/Z_o`

Due to the sprung mass on the tyre, there will be a static deflection on the tyre and due to this deflection, the tyre gets flattened at the contact patch and loses contact with the road. If the deflection is zero that is the tyre regains its original shape then the normal force will also be zero. Beyond this point, the tyre will be lifted from the ground. 

For a better design, there should be a tradeoff between these three transmissibility ratios so that the suspension will be optimized.

State Space Representation  :

 `dx_1/dt = a_(11) + a_(12)x_2+b_1v(t)`

  `dx_2/dt = a_(21) + a_(22)x_2+b_2v(t)`

 `y = c_1x_1 +c_2x_2+d_1v(t)`

`((dotx_1),(dotx_2)) = [[a_(11),a_(21)],[a_(21),a_(22)]] ((x_1),(x_2))+((b_1),(b_2))*v(t)`

`y = [[c_(1),c_2]][(x_1),(x_2)] + d_1*v(t)`

 `x_1` and `x_2` are called state variables

A = `[[a_(11),a_(21)],[a_(21),a_(22)]]` `to` System matrix

B, G = `((b_1),(b_2))``to` Input matrix

C =  `[[c_(1),c_2]]``to` Output matrix

D = `d_1` `to` feed forward matrix

Quater car model Equations of motion in State Space Representation :

`m_sddotZ_s+ c_s(dotZ_s - dot Z_(us)) + k_s(Z_s - Z_(us)) = - F`

`m_(us)ddot(Z_(us)) + c_s(dotZ_(us) - dotZ_s) + k_s(Z_(us) -Z_s)+c_(us)(dotZ_(us)-dotZ_o)+k_(us)(Z_(us) - Z_o) = F`

`d/dt[[Z_(us)-Z_o],[dotZ_(us)],[Z_(us)-Z_o],[dotZ]] =[[0,1,0,0],[-k_(us)/m_(us),-(c_s+c_(us))/m_(us),k_(us)/m_(us), c_s/m_(us)],[0,-1,0,1],[0,c_s/m_s,-k_s/m_s,-c_s/m_s]][[Z_(us)-Z_o],[dotZ_(us)],[Z_(us)-Z_o],[dotZ]]+ [(0),(m_s/m_(us)),(0),(-1)]F/m_s+[(-1),(c_(us)/m_(us)),(0),(0)]dotZ_o`

A `to` `[[0,1,0,0],[-k_(us)/m_(us),-(c_s+c_(us))/m_(us),k_(us)/m_(us), c_s/m_(us)],[0,-1,0,1],[0,c_s/m_s,-k_s/m_s,-c_s/m_s]]`

B `to` `[(0),(m_s/m_(us)),(0),(-1)]F/m_s`

G `to` `[(-1),(c_(us)/m_(us)),(0),(0)]dotZ_o`

Since `dotx` = Ax + Bu + Gw

            `y  = Cx + DG

Outputs `to` C

Tyre deflection = `[[1,0,1,0]]`

Suspension stroke = `[[0,0,1,0]]`

Sprung acceleration = `[[0 ,c_s/m_s, -k_s/m_s, -c_s/m_s]]` 

`rho = m_s/m_(us)

`omega_1 = sqrt(k_s/m_s)  `
`omega_2 = sqrt(k_(us)/m_(us))  `
`zeta_1 = c_s/( 2m_somega_1)`  
`zeta_2 = c_s/( 2m_(us)omega_2)`

A `to` `[[0,1,0,0],[-omega_(2)^2,-2(rhozeta_1omega_1+zeta_2omega_2),rhoomega_(1)^2,2rhozeta_1omega_1],[0,-1,0,1],[0,2zeta_1omega_1,-omega_(1)^2,-2]]`

B `to` `[[0,rho,0,-1]]`

G `to` `[[-1,2zeta_2omega_2,0,0]]`

clear all
close all
clc

% x(1) = tire deflection (zus-z0) 
% x(2) = velocity of unsprung mass (d(zus)/dt)
% x(3) = suspension stroke (zs-zus)
% x(4) = sprung mass speed (d(zs)/dt)


% parameters
w1 = 2*pi; % w1 = sqrt(ks/ms)
w2 = 20.0*pi; % w2 = sqrt(kus/mus)
z1 = 0.3; % z1 = cs/(2*ms*w1)
z2 = 0.0; % z2 = cus/(2*mus*w2)
rho = 10.; % rho = ms/mus

 % Open loop system equations:
 A = [0 1 0 0
     -w2^2 -2*(z1*w1*rho+z2*w2)  rho*w1^2     2*z1*w1*rho
     0 -1 0 1
     0 2*z1*w1 -w1^2 -2*z1*w1];
 B = [0 rho 0 -1]';
 G = [-1 2*z2*w2 0 0]';
 
 �fine outputs of interest:
 % output=tire deflection
 
C1=[1 0 0 0]; D1= 0.0; 
[num1, den1]=ss2tf(A,G,C1,D1,1);

% output=suspension stroke
C2=[0 0 1 0]; D2= 0.0; 
[num2, den2]=ss2tf(A,G,C2,D2,1);

% output=sprung mass acceleration.
C3=[A(4,:)]; D3= G(4); 
[num3, den3]=ss2tf(A,G,C3,D3,1);


%Generate the white noise input w(t)
t=[0:0.001:1]; 
w=0.1*randn(size(t));   % Assume road velocity is
          % white Gaussian with 
                              % standard deviations of 0.1 (m/sec)
                              
                              
                              
%% Simulate the response of interest:
y1=lsim(num1,den1,w,t);
y2=lsim(num2,den2,w,t);
y3=lsim(num3,den3,w,t);
clf;
figure(1), subplot, subplot(221)
plot(t,w,'r')
xlabel('Time (sec)'); ylabel('Road speed (m/sec)')
subplot(222), plot(t,y1,'r'); 
xlabel('Time (sec)'); ylabel('Tire Def (m)')
subplot(223), plot(t,y2,'r'); 
xlabel('Time (sec)'); 
ylabel('Susp Stroke (m)')
subplot(224), plot(t,y3,'r'); 
xlabel('Time (sec)'); 
ylabel('Sprung. mass accel (m/sec^2 )')       


figure(2)
freq=logspace(-1,2.7,200); ^-1 to 10^2.7 rad/sec, 200 points
[mag1, phase1]=bode(num1,den1,freq);
loglog(freq/(2*pi),mag1,'r');
title('Frequency Response Magnitude');
xlabel('Frequency (Hz)'); 
ylabel('Tire def (m/(m/sec))'); 
grid


figure(3)
[mag2, phase2]=bode(num2,den2,freq);
loglog(freq/(2*pi),mag2,'r'); 
title('Frequency Response Magnitude');
xlabel('Frequency (Hz)'); 
ylabel('Susp stroke (m/(m/sec))'); 
grid

figure(4)
[mag3, phase3]=bode(num3,den3,freq);
loglog(freq/(2*pi),mag3,'r'); 
title('Frequency Response Magnitude');
xlabel('Frequency (Hz)'); 
ylabel('Sprung mass accel (m/sec^2 /(m/sec))'); 
grid;