A quarter car model with additional massless degree of freedom for the top mount

Hello,

I'd like to set up and solve a set of differential equations for a quarter car model but with an additional spring in series with the damper to represent top mount stiffness, as per the following image:

I already have a standard quarter car model running so have just modified these existing equations to include an extra massless degree of freedom, , between the suspension damper and the top mount spring. However, something doesn't seem right with the equations, not least because the whole lot disappears towards infinity very quickly. This might suggest I've made errors with force polarities but I can't find them.

My equations are:

The first two are rearranged for the sprung mass () and unsprung mass () accelerations by dividing through by the masses. The last equation is rearranged for the velocity of the massless element ():

As well as being one of the differential equations, this is also inserted into the equation for the unsprung mass (), but this seems to set equal to zero, also suggesting something isn't right.

My ODE code is:

dx(1,:) = x(3,:); % sprung mass velocity
dx(2,:) = x(4,:); % unsprung mass velocity
dx(3,:) = 1./ms.*(ks.*(x(2,:)-x(1,:)) + kp.*(x(5,:)-x(1,:))).*1000; % sprung mass acceleration
dx(4,:) = 1./mu.*((kt.*(i(1,:)-x(2,:))) + (ct.*(idot(1,:) - x(4,:))) + (ks.*(x(1,:)-x(2,:))) + (cs.*((x(4,:)-(kp.*(x(1,:)-x(5,:))./cs))-x(4,:)))).*1000; % unsprung mass acceleration
dx(5,:) = x(4,:)-(kp.*(x(1,:)-x(5,:))./cs); % top mount velocity
dx(6,:) = 0;

Is my approach to the problem correct?

Happy to post the full code if that helps.

Thanks.

[Sam Chak]

Hi @Simon Aldworth

Using the free-body diagram, you can apply Lagrangian mechanics to derive the Euler–Lagrange equations and compare the simplified version with the ODEs. If your modeling is accurate, the results should be consistent. Best to use Symbolic Math Toolbox for this task. You might consider using the functionalDerivative() function, and see relevant examples in the documentation.

Post the MATLAB code and results here so that other users can assist.

Answer 1

[Paul]

Hi @Simon Aldworth

Seems like the equation for the massless element would be:

m5*x5ddot = kp * (x1 - x5) + cs * (x2dot - x5dot)

m5 = 0 -> 0 = kp * (x1 - x5) + cs * (x2dot - x5dot) -> kp * (x5 - x1) = cs * (x2dot - x5dot) (note polarity of term on LHS)

x5dot = x2dot - kp / cs * (x5 - x1)

m2 * x2ddot = kt * (xi - x2) + ct * (xidot - x2dot) + cs * (x5dot - x2dot)

= kt * (xi - x2) + ct * (xidot - x2dot) - kp * (x5 - x1)

m1 * x1ddot = ks * (x2 - x1) + kp * (x1 - x5)

So we should have a system a five differential equations with state variables: z1 = x1, z2 = x1dot, z3 = x2, z4 = x2dot, and z5 = x5, with inputs u1 = xidot and u2 = xi.

If simulating this system one must maintain the mathematical relationship between u1 and u2. One way to do that would be to define xidot as the input to they system and add an additional state variable z6 = xi. Kind of depends on how xi(t) is specified.

[Simon Aldworth]

Hi @Paul,

Many thanks for spotting the sign error - the simulation runs fine now.

Note that you have missed out the component relating to ks for m2, but no matter to the outcome.

Interested to know what you meant by maintaining the mathematical relationship between the input displacement and veloctiy. Could you expand on that please?

My input is a chirp function and the velocity is calculated directly from it:

t0 = 0; % s. Input start time
t1 = 81.9; % s. Input end time
dt = 0.002; % s. Time resolution
tvec = t0:dt:t1; % horizontal vector between t0 and t1 incremented by dt
A = 10; % mm. Input amplitude
f0 = 0.01; % Hz. Input start frequency
f1 = 20; % Hz. Input end frequency
g = (f1./f0).^(1./(t1-t0)); % Exponential growth of chirp frequency
i = A.*sin(f0.*(((g.^tvec)-1)./log(g)).*2.*pi); % Ground input displacement - exponential chirp signal
idot = A.*cos(f0.*(((g.^tvec)-1)/log(g)).*2.*pi).*(2.*pi.*f0.*g.^tvec); % Ground input velocity - exponential chirp signal

plot(tvec, i)
xlabel('Time')
ylabel('Amplitude')
title('Ground input displacement - exponential chirp signal')

Thanks again.

Simon.

[Paul]

Yep, I missed the ks term on m2. Good catch.

From the original question I didn't know how you were getting i(t) and idot(t). Because you have closed form expressions for both, that should be fine. If you didn't have a closed form expression for i(t), e.g., a table lookup or something, you'd have to be careful to ensure that idot(t), however obtained, is really the derivative of i(t) when both are presented in the model.

[Simon Aldworth]

Ah, I see. I did wonder if that was what you meant because I have in the passed used a step input consisting of a simple vector of values in which a jump occurs over one time step. In this case, what do you propose I do to calculate idot? I have used the gradient function most of the time and accepted any inaccuracies.

[Sam Chak]

@Simon Aldworth, Since you are modeling a quarter car, the signal xi should represent the road surface profile, which acts as an "input disturbance" to the suspension dynamics that is generally not measured. For simulation purposes, it is acceptable to use closed-form mathematical formulas to describe the input disturbance. For example, a speed bump can be characterized using the Bump function.

Image source: https://www.sinoconcept.co.uk/car-park-safety-management/speed-bumps-and-speed-humps/one-part-speed-bump/

hw = 150;                               % half-width of the bump (mm)
x  = linspace(-2*hw, 2*hw, 60001);
y1 = exp(-1./(1 - (x/hw).^2));          % component 1
y2 = (sign(x + hw) - sign(x - hw))/2;   % component 2
H  = 50;                                % max height of the bump (mm)
A  = H/0.367879441171442;
y  = A*y1.*y2;                          % Speed Bump function 
plot(x, y, 'linewidth', 1.5), grid on, ylim([0 600])
title('5-cm height Speed Bump')
xlabel({'$x$/mm'}, 'interpreter', 'latex')
ylabel({'$y$/mm'}, 'interpreter', 'latex')

[Paul]

If i(t) is a step input, say at t = t0, and if the differential equation is affine in i(t) and idot(t), then the best way to solve the problem is to numerically integrate to t = t0 with i(t) = idot = 0, then integrate the solution "by-hand" a t = t0 to deal with idot = delta(t-t0), then continue the integration from t = t0 with the updated solution as the intitial condition.

If we have:

xdot = f(x,t) + g1(x,t)*i + g2(t)*idot w/ i = a*u(t) (step)

Then integrate:

xdot = f(x) from t = 0 to t0. (1)

then update x (assuming g2(t) is continuous at t = t0)

x = x + g2(t0)*a

Then continue the integration

xdot = f(x,t) + g1(x,t)*a (2)

Then append the solutions from (1) and (2).

I'd have to think a bit more about how to do this if g2 is also a function of x (I think I've thought about this in the past), though I'd be curious what type of physical system would have that structure.

See also impulse response ode45 help - MATLAB Answers - MATLAB Central for a similar discussion for the case where the system is LTI.

[Paul]

Hi Sam,

Question about how handle the bump function in this application.

First, I'm going to recast that function in the notation of this question where xi(t) is the height of the road (what you called y) and d defines the horizonal distance along the road (what you called x)

syms d h_w A real
assumeAlso(h_w,'positive')
x1 = exp(-1./(1 - (d/h_w).^2));          % component 1
x2 = (sign(d + h_w) - sign(d - h_w))/2;   % component 2
xi  = A*x1.*x2

xi = 

But d is really a function of t, so for simplicity I'll assume the that the car is traveling at a constant horizontal speed

syms v t  % assuming v ~= 0
xi(t) = subs(xi,d,v*t)

xi(t) = 

Now we have the expression for xi(t) as would be needed for simulation. But we also need the expression for its derivative

dxi = expand(diff(xi,t))

dxi(t) = 

Focusing on the second term (which is synonymous to the first)

c = children(dxi(t));
c2 = simplify(c{2})

c2 = 

I'm not sure how to handle that expression in simulation. We can further simplify the Dirac term.

c2 = subs(c2,dirac(h_w-t*v),dirac(t-h_w/v)/abs(v))

c2 = 

My concern is that the exponential term is discontinuous at t = hw/v. In that situation, one might argue that we could replace that term with the average of its limit on either side of the discontinuity. But in this case the limit is undefined as t approaches hw/v from the right.

[Sam Chak]

Hi @Paul

Thank you for your input. You are correct. I overlooked the issue that the horizontal distance traveled by the car must be a function of time t, and the bump function complicates the simulation due to discontinuities.

In vehicle suspension design, an engineer typically selects the spring stiffness (K) and damping coefficient (C) to optimize ride comfort, handling/safety, and durability for intended use on smooth asphalt with occasional speed bumps, in line with the majority of customer preferences.

Since the comfort goal is to minimize vertical acceleration and vibration transmitted to the occupants, it may be acceptable to use a simplified, infinitely differentiable analytic bump shape that provides a reasonably accurate representation of the bump profile for the quarter car simulation purposes.

hw  = 150;                              % half-width of the bump (mm)
d   = linspace(-2*hw, 2*hw, 60001);     % horizontal distance travelled by the car

% Original Bump function
x1  = exp(-1./(1 - (d/hw).^2));         % component 1
x2  = (sign(d + hw) - sign(d - hw))/2;  % component 2
H   = 50;                               % max height of the bump (mm)
A   = H/0.367879441171442;              % scaled amplitude of the Bump function 
xi  = A*x1.*x2;                         % complete bump profile

% Simplified, infinitely differentiable everywhere, analytic bump shape
sig = 85;                               % sigma             84.7743762163065
pow = 4;                                % power exponent    3.558511911363521
xis = H*exp(- 1/2*(d/sig).^pow);        % accurate version: H*exp(- 1/2*abs(d/sig).^pow)

plot(d, [xi; xis]), grid on, ylim([-300 300])
title('50-mm height Speed Bump profile')
xlabel({'$d$ / mm'}, 'interpreter', 'latex')
ylabel({'$x_{i}$ / mm'}, 'interpreter', 'latex')
legend('Original bump function', 'Simplified bump function')

However, it must be assumed that the car is traveling at very low speeds, allowing the tire to remain in contact with the road surface and the contact point (xi) to closely follow the bump profile.

At higher speeds, the contact can break, causing the tire's contact path to deviate from the static bump profile. If the wheel lifts off and then re-establishes contact, the tire contact path often resembles a lowercase "" (the wheel follows the bump up, then experiences free flight over the crest, and finally impacts back onto the downslope). Therefore, using the bump surface alone is inappropriate in that scenario.