How to model a system made by a vehicle suspension and tire?

1. Apr 21, 2012

newlife

Hello everybody,
Im trying to model a system made by a suspension and a tire for a vehicle simulator written in c#.
I modeled both the suspension and the tire as a spring and a damper.
So the model is formed by a spring and a damper in parallel (the suspension) in series with another spring and a damper (the tire).
This is a simple scheme:

http://unitypackages.net/spring_damper.png [Broken]

I modeled the system in this way:
I throw a ray from the point where the suspension is connected to the vehicle towards the ground with length = suspension travel + wheel radius.
If the ray touches the ground, suspension compression is calculated in this way:

suspension force is:
springForce = suspensionRate*compression;

suspension damper force is:
damperForce = (damperVelocity - tireDeflectionVelocity)*damperRate;

with damperVelocity:
damperVelocity=(compression-lastCompression)/deltaTime;

total suspension force is springForce+damperForce.

For what concern the tire:
tireDeflection = (suspensionRate*compression)/verticalTireStiffness;

tireForce=verticalTireStiffness*tireDeflection;
damperForce = tireDeflectionVelocity * tireDamping;

with tireDeflectionVelocity:
tireDeflectionVelocity=(tireDeflection - lastTireDeflection)/deltaTime;

total tire force is tireForce+damperForce.

Everything works like it should, apart when suspensionRate is > verticalTireStiffness (it may happen when suspension is very stiff and tire is low inflated). In this case the system explode, in the sense that wheels shakes violently.

Is there a way to avoid this issue? I tested even with a very small timestep but it seems its not a problem of large timestep.

Sources:

Last edited by a moderator: May 5, 2017
2. Apr 21, 2012

chiro

Hey newlife and welcome to the forums.

If you have not taken a course on differential equations or numerical calculus then you may not be aware that a numerical calculation may not be stable.

The area of stability is one where you consider when a particular differential equation, a numerical evaluator and the parameter values of the DE and numerical evaluator work in a way that you get answers that 'make sense'.

If you get a situation where stability is not guaranteed then you get the kinds of problems that you are seeing.

What I would do is investigate a few different numerical schemes. A popular one is known as the Runge-Kutta method which is going to be a lot more accurate.

The basic idea with numerical integrators is that you have a system that models the changes of all the variables and then you give this information to a numerical integrator which calculates the outputs at certain points and then you will take these outputs at each computation step and then update the changes visually.

3. Apr 21, 2012

newlife

Hello chiro,
thank you for your reply. Im aware about differential equations and in particular 4th order Runge-Kutta method. I have used it for slip ratio calculation, which is very prone to instability.
The problem is that, as i wrote in the message, even if i choose a very small time step (0.001, 1000 iterations per sec) the problem is still present and, strangely enough, it seems even worst.
Btw, my standard time step is 0.02, 50 iterations per sec.

Some results of real cases test:

Susp. Rate: 50000 N/m
Tire Rate: 100000 N/m
OK

Susp. Rate: 100000 N/m
Tire Rate: 200000 N/m
OK

Susp. Rate: 100000 N/m
Tire Rate: 100000 N/m
KO-> oscillations

Susp. Rate: 80000 N/m
Tire Rate: 80000 N/m
KO-> oscillations

Susp. Rate: 70000 N/m
Tire Rate: 70000 N/m
KO-> oscillations

Susp. Rate: 60000 N/m
Tire Rate: 60000 N/m
KO-> oscillations

Susp. Rate: 50000 N/m
Tire Rate: 50000 N/m
OK

sidewall height 0.068 meters
weight: 1300 kg

4. Apr 21, 2012

chiro

I know you have written this implicitly in the code, but could you please post the DE relationships for your system just to clarify what you are calculating.

Again I'm aware that the code implicitly says this but it's just for clarity and context.

5. Apr 22, 2012

newlife

Sorry chiro what do you mean with DE relationship?

6. Apr 22, 2012

chiro

I just meant all the DE equations for all of the variables.

7. Apr 23, 2012

newlife

ok.

Code (Text):

// Spring calc
compression = suspensionTravel - springLen;

// Tire calc
rate=suspensionRate/verticalTireStiffness;
lastTireDeflection=tireDeflection;
tireDeflection = springForce/verticalTireStiffness;
tireDeflection=Mathf.Clamp(tireDeflection, 0, sidewallHeight);

float forces=verticalTireStiffness*tireDeflection;

oldDeflectionVelocity=deflectionVelocity;
float delta=tireDeflection - lastTireDeflection;
if (Mathf.Abs(delta)<0.00001f) delta=0;
deflectionVelocity=delta/Time.deltaTime;

damperForce = (deflectionVelocity) * tireDamping;

forces+=damperForce;

// Suspension Calc
springForce = suspensionRate*compression;

float damperForce;
float forces=springForce;

normalVelocity = normalVelocityR - deflectionVelocity; // normalVelocityR is the dot product between wheel velocity and wheel normal
damperForce = normalVelocity*damperRate;

forces += antiRollBarForce;
forces += damperForce;

Thats all. Tell me if you need something else

8. Apr 24, 2012

rvangaal

RK4 isn't going to do much for you, since all the parameters involved in a car's system make it too complex to just hop around in time trying to get valid states for each subtime.

A step of 0.001 is not bad; just 50Hz will get you nowhere. Try sticking to 0.001 step for a while (I use that in Racer; we also use that upto high-end users, although some even move towards higher values). Above 0.0001 you'll lose accuracy due to floating point accuracy.

I'd get rid of statements such as 'if (Mathf.Abs(delta)<0.00001f) delta=0;'. These normally do more harm than good and cause problems that are hard to fix. Just let the model jitter; it should be invisible visually anyway.

Also, you use:

float delta=tireDeflection - lastTireDeflection;
deflectionVelocity=delta/Time.deltaTime;

This means you're integrating velocity; you're losing a timestep here, which may cause the instability problems. It might be an idea to make the tire and wheel points actual point masses that have velocity, and update velocity through accelerations, rather than integrating directly (phase issues will probably occur).

Make sure your damping rate isn't too high; check how near your damping is to optimal damping (a unit from 0..1 rather than a damperRate). Although more damping normally just means overdamping, not quickly exploding.

Also try and only do suspension forces (no damping/ARB).

I'm doing explicit tire physics (implicit is highly complex) with just Euler integration for tire rates upto 300,000N/m and suspension rates of upto 150,000N/m which works. Good enough for Formula-style cars.

9. Apr 24, 2012

Bob Smith

Similar to what Ruud suggested, I would either model the wheel as a separate mass in the system, or just assume the wheel will always be in the resting position each frame, between the two springs (i.e. use combined stiffness to work out overall deflection, then use the individual stiffnesses to work out the individual spring deflections).

10. Apr 30, 2012

newlife

Hi Ruud, first of all i would like to thank you for the invaluable help you gave me (and all the sim developers) with your excellent web site, really packed with very useful and precise information for anyone interested in the argument.

Well im glad to tell you that I succeeded in making everything work even at 0.02!

This is what helped me to fix my problem!! thank you very much Ruud.
now deflectionVelocity is calculated like this:

accel = (normalSuspensionForce - normalTireForce)/mass;
deflectionVelocity += accel*Time.deltaTime; (in my case, i had to modify it with deflectionVelocity = accel*Time.deltaTime; it works much better)

More, I had to substitute this:

with this:
springLen += (deflectionVelocity - normalVelocityR) * Time.deltaTime;

cause springLen is influenced by the relative velocity of the system.

which is a realistic tire damping value? I found online a value of about 0.07 of critical damping (http://answers.yahoo.com/question/index?qid=20080422084944AA55GXu) but it seems a little too low for a correct behaviour of the tire in the sim. I choose a value of 0.3.

I succeded to obtain the more or less the same result (suspension rates of upto 150,000 N/m, tire rates upto 350,000 N/m, bump and rebound rate upto 12000 Ns/m) with 0.02 time step, with the trick of multiplying the mass value in the acceleration calc by a value, in order to damp the oscilation. Original mass value is 20 kg, now is 140 kg. Should affect much the simulation, it just damps wheel acceleration. What do you think?
Again, thanks for the help.
Michele

11. Apr 30, 2012

newlife

Hi Bob can you explain me better what you mean with "assume the wheel will always be in the resting position each frame"? it seems interesting but i cant get it..

12. May 4, 2012

rvangaal

You're welcome. :-) And good to hear you get it working at 50Hz, although with that high tire mass, hm.

Do realize here that you implicitly simulate the wheel moving up/down by integrating 'springLen'. Instead of the instantaneous you integrate. Quite ok if it works for you; there are dozens of paths to Rome.

For a tire rate of 300,000N/m, try 3000 N/m/s. Real tires damp more as velocity is low; at 0 m/s use 5000 for example, then go back to 2500 N/m/s when contact patch velocity reaches say 5 m/s. A stiffer rotating tire damps less, and we needed this for Formula style cars to behave stable. But that depends a lot on the other math. :-)

Doing this makes it hard to tune though; I'd try and split the 2 effects, using dampening on your springVelocity. Use F=m*a like with deflectionVelocity but separate the forces; this also makes it easier to cut out effects if you want to test things.
So add a 'springDampRate', then forceSpringDamp=springVelocity*springDampRate and add all forces to move the wheel center.

Just multiplying mass by 7 makes a tire stable, sure, but not for the right reason.