I’m developing a simulation software for flight simulation.
I have problems with solver stability.

I have a rigidbody subject to gravity and aerodynamic forces. My timestep is fixed to 0.01s and lift and drag forces are computed at timestep n using body linear and angular velocity. (I compute AoA and relative speed body/fluid in wing coordinate system) and apply to the actor with simply functions: AddForceAtLocalPosition for example.
after few timesteps (200/300) the simulation diverge and velocities goes to infinite. If I lower the timestep obviously the instability come up at a later moment.

I’ve tried to use function as:
ComputeVelocityDeltaFromImpulse and use velocity at timestep n+1 to compute forces in order to try an “implicit method” but this not improve the simulation stability.

I understood that Physx solver use semi implicit integrator method but i don’t if i can modify some parameters about the solver.

Any method to improve solver stability without reducing the timestep too much in this case? There a lot of simulator and i can’t figure out how they deal such stability problems for realtime simulations?

I suppose you’re using quadratic (or other non-linear) equations to calculate air drag. I also had problems with this and settled for an incorrect but believable linear model. Can you show your equations? So from rb.GetPointVelocity to AddForceAtLocalPosition. Does your system also become instable at low velocities?

my code is quite simple, i use local velocity in airfoil coordinate system in order to calculate angle of attack so interpolate Cl,Cd coeffiecients and calculate lift and drag as:
Lift = 0.5AreaCl*Vxz^2

obiously is a not linear model. I’m interested in make a realistic simulation so don’t want to change formulas.

Yes, after some seconds a quite slow speed (5/6 m/s) the simulation become unstable and diverge with a fixed timestep of 0.01s.

v *= -1; //Get flow velocity
v += Vinf; //Add flow velocity at infinite
v.Transform(trf); //Velocity in local BladeElement coord system
v_dir = new Vector3d(v);
v_dir.Unitize();
geomAoA = blElem.ComputeGeomAoA(v);
cl = blElem.fs.ComputeCl(geomAoA);
cd = blElem.fs.ComputeCd(geomAoA);
cm = blElem.fs.ComputeCm(geomAoA);
Debug.WriteLine("AoA = " + geomAoA*180/Math.PI);
Debug.WriteLine("Cl = " + cl);
Debug.WriteLine("Area = " + blElem.Area);
q = 0.5 * rho * blElem.Area;
lift = q * (v.X * v.X + v.Z * v.Z) * cl; //Neglect spanwise velocity vy
drag = q * (v.X * v.X + v.Y * v.Y + v.Z * v.Z) * cd; //Use Full velocity
lift_dir = Vector3d.CrossProduct(new Vector3d(0.0,1.0,0.0), v_dir); //Lift direction perpendicular to Velocity and local Y axis
lift_dir.Unitize();
Vector3 localElemtoCG = blElem.GetLocalCollocationPoint().ToSystemVector() - ibody.cgPosition;
Debug.WriteLine("local CG=" + ibody.cgPosition.ToRhinoPoint().ToString());
Debug.WriteLine("localElemtoCG"+localElemtoCG.ToString());
lift_vec = lift * lift_dir;
drag_vec = drag * v_dir;
force_i = lift_vec + drag_vec; //force in BladeElem coord system
//Rotate to World coord system
force_i.Transform(trf_inv);
moment_i = Vector3d.CrossProduct(blElemToCG, force_i);
//moment_i.Transform(trf_inv);
lift_vec.Transform(trf_inv);
drag_vec.Transform(trf_inv);
forceV += force_i;
momentV += moment_i;