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

[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

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:
compression = suspensionTravel - (rayLength - radiusLoaded);

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.

radiusLoaded=maxRadius - tireDeflection;

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:
https://docs.google.com/viewer?a=v&...jK471r&sig=AHIEtbSih5uKFNJzLgXZjBH1wFwQ6L93qQ
https://www.physicsforums.com/showthread.php?t=225036

 

[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.

 

[newlife]

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.

Hello chiro,
thank you for your reply. I am 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

wheel radius: 0.34 meters
sidewall height 0.068 meters
weight: 1300 kg

 

[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.

 

[newlife]

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.

Sorry chiro what do you mean with DE relationship?

 

[chiro]

Sorry chiro what do you mean with DE relationship?

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

 

[newlife]

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

ok.

// Spring calc
SpringLen=(rayLength - radiusLoaded); 
compression = suspensionTravel - springLen;


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

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

 

[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.

 

[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).

 

[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.

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.

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

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).

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:
springLen=(rayLength - radiusLoaded);

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

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

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.

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.

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.

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

 

[newlife]

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).

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 can't get it..

 

[rvangaal]

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 I am glad to tell you that I succeeded in making everything work even at 0.02!

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

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:
springLen=(rayLength - radiusLoaded);

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

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

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.

which is a realistic tire damping value?

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. :-)

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?

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.