When Update the wheel and engine speeds - 5x5 matrix why A.set(i,i,...) in this way?

after compute the tire force/torque in processSuspTireWheels

the tire torque = -fz*wheelRadius,fz is the longitudinal tire force

update the wheel and engine angular speed
in func solveDrive4WInternaDynamicsEnginePlusDrivenWheels

everything is clear except A.set(i,i,…)

//express ai as 
//ai = [wi(t+dt) - wi(t)]/dt
//and rearrange
//wi(t+dt) - wi(t)] = dt*G*gammai*K*{G*[alpha0*w0(t+dt) + alpha1*w1(t+dt) + alpha2*w2(t+dt) + ..... alpha(N-1)*w(N-1)(t+dt)] - wEng(t+dt)}/Ii + dt*(bt(i) + tt(i))/Ii

A.set(i,i,1.0f+dtKGGR*aveWheelSpeedContributions[i]+dt*wheels4SimData.getWheelData(i).mDampingRate);

in the equation ,there is no wheelData.mDampingRate, Only 1.0+dtKGGR*aveWheelSpeedContributions[i] is clear
so how to make the wheelData.mDampingRate into the equation ?

You just need to add an extra damping term to ai. The extra damping term is D_i * alpha_i/Ii, where D_i is the wheel damping rate.

The correct equation will be:

//acceleration applied to ith wheel is 
//ai = G*gammai*K*{G*[alpha0*w0 + alpha1*w1 alpha2*w2 + ..... alpha(N-1)*w(N-1)] - wEng}/Ii + (bt(i) + tt(i))/Ii - <b>D_i*alpha_i//Ii</b>

If you work that through it isn’t hard to show that the diagonal term has an extra term:

  • dt*D_i/Ii

In code we’ve already divided dt by Ii so the extram term is just

+ dt*D_i/Ii

Hope this helps.

Thanks !

There is still a question:
After create the suspension limit constraint,
the constraint is add to Sc::Scene
But where the constraint is solved?
There are so many threads, I can’t find the function to solve the constraint

Also what’s the difference between a hard constraint and spring constraint in PhysX?

Here is the Entire Derivation of soft constraint in Bullet post by Erin Catto.
softness and penetration(bias) → (CFM and ERP)

Newton’s Law:
M * (v2 - v1) = JT * lambda

Velocity constraint with softness and Baumgarte:
J * v2 + softness * lambda + bias = 0

where bias = biasFactor * C / dt

We know everything except v2 and lambda.

First solve Newton’s law for v2 in terms of lambda:

v2 = v1 + invM * JT * lambda

Substitute this expression into the velocity constraint:

J * (v1 + invM * JT * lambda) + softness * lambda + bias = 0

Now collect coefficients of lambda:

(J * invM * JT + softness) * lambda = - J * v1 - bias

Now we define:

K = J * invM * JT + softness

and Meff = invK is the effective mass.

Now keep in mind that v1 contains the velocities for both bodies:

v1 = column_vector(vb1, omegab1, vb2, omegab2)

Now here is where things get a bit hairy. Each iteration we are not computing the whole impulse, we are computing an increment to the impulse and we are updating the velocity. Also, as we solve each constraint we get a perfect v2, but then some other constraint will come along and mess it up. So we want to patch up the constraint while acknowledging the accumulated impulse and the damaged velocity. To help with that we use P for the accumulated impulse and dP as the update. Mathematically we have:

M * (v2new - v2damaged) = JT * dP
J * v2new + softness * (P + dP) + bias = 0

If we solve this we get:

v2new = v2damaged + invM * JT * dP
J * (v2damaged + invM * JT * dP) + softness * P + softness * dP + bias = 0

(J * invM * JT + softness) * dP = -(J * v2damaged + softness * P + bias)

The constraints are solved in physx. If you have 3.4 then you will need to look in LowLevelDynamics to see the solver code. Functions like solve1D in DySolverConstraints.cpp are the place to start looking. Be aware, this is highly optimised code so following the computation will be harder than the vehicle code.

I recommend reading the User’s Guide on constraints and joints to learn more about the difference between hard and soft constraints in physx.