Add support for AC induction motor

Author: cesarBLG

Source/Orts.Simulation/Orts.Simulation.csproj

@@ -194,6 +194,7 @@
     <Compile Include="Simulation\Turntables.cs" />
     <Compile Include="Simulation\Weather.cs" />
     <Compile Include="Common\Scripting\EuropeanTrainControlSystem.cs" />
+    <Compile Include="Simulation\RollingStocks\SubSystems\PowerTransmissions\InductionMotor.cs" />
   </ItemGroup>
   <ItemGroup>
     <ProjectReference Include="..\ORTS.Common\ORTS.Common.csproj">

Source/Orts.Simulation/Simulation/RollingStocks/MSTSLocomotive.cs

@@ -442,10 +442,10 @@ public float OdometerM
 
         public enum TractionMotorTypes
         {
+            DC,
             AC,
-            DC
         }
-        public TractionMotorTypes TractionMotorType;
+        public TractionMotorTypes TractionMotorType = TractionMotorTypes.DC;
 
         public ILocomotivePowerSupply LocomotivePowerSupply => PowerSupply as ILocomotivePowerSupply;
         public ScriptedTrainControlSystem TrainControlSystem;
@@ -1460,7 +1460,11 @@ public override void Initialize()
                     Trace.TraceInformation("Number of Locomotive Drive Axles set to default value of {0}", LocoNumDrvAxles);
                 }
             }
-               
+            if (TractionMotorType == TractionMotorTypes.AC)
+            {
+                InductionMotor motor = new InductionMotor(LocomotiveAxle, this);
+            }
+
 
             // Calculate minimum speed to pickup water
             const float Aconst = 2;
@@ -1873,6 +1877,7 @@ public override void Update(float elapsedClockSeconds)
             }
             else
                 DynamicBrakeForceN = 0; // Set dynamic brake force to zero if in Notch 0 position
+                
 
             UpdateFrictionCoefficient(elapsedClockSeconds); // Find the current coefficient of friction depending upon the weather
 
@@ -2242,11 +2247,6 @@ protected virtual void UpdateTractiveForce(float elapsedClockSeconds, float t, f
                 if (TractionMotorType == TractionMotorTypes.AC)
                 {
                     AbsTractionSpeedMpS = AbsSpeedMpS;
-                    if (AbsWheelSpeedMpS > 1.1 * MaxSpeedMpS)
-                    {
-                        AverageForceN = TractiveForceN = 0;
-                        return;
-                    }
                 }
                 else
                 {
@@ -2644,21 +2644,30 @@ public void AdvancedAdhesion(float elapsedClockSeconds)
                 LocomotiveAxle.DampingNs = MassKG / 1000.0f;
                 LocomotiveAxle.FrictionN = MassKG / 1000.0f;
 
-                if (SlipControlSystem == SlipControlType.Full)
+                if (LocomotiveAxle.Motor is InductionMotor motor)
                 {
-                    // Simple slip control
-                    // Motive force is reduced to the maximum adhesive force
-                    // In wheelslip situations, motive force is set to zero
-                    MotiveForceN = Math.Sign(MotiveForceN) * Math.Min(LocomotiveAxle.AdhesionLimit * LocomotiveAxle.AxleWeightN, Math.Abs(MotiveForceN));
-                    if (LocomotiveAxle.IsWheelSlip) MotiveForceN = 0;
+                    motor.SlipControl = SlipControlSystem == SlipControlType.Full;
+                    motor.TargetForceN = MotiveForceN;
+                    motor.EngineMaxSpeedMpS = MaxSpeedMpS;
+                }
+                else
+                {
+                    if (SlipControlSystem == SlipControlType.Full)
+                    {
+                        // Simple slip control
+                        // Motive force is reduced to the maximum adhesive force
+                        // In wheelslip situations, motive force is set to zero
+                        MotiveForceN = Math.Sign(MotiveForceN) * Math.Min(LocomotiveAxle.AdhesionLimit * LocomotiveAxle.AxleWeightN, Math.Abs(MotiveForceN));
+                        if (LocomotiveAxle.IsWheelSlip) MotiveForceN = 0;
+                    }
+                    LocomotiveAxle.DriveForceN = MotiveForceN;              //Total force applied to wheels
                 }
 
                 //Set axle model parameters
 
                 // Inputs
                 LocomotiveAxle.BrakeRetardForceN = BrakeRetardForceN;
                 LocomotiveAxle.AxleWeightN = 9.81f * DrvWheelWeightKg;  //will be computed each time considering the tilting
-                LocomotiveAxle.DriveForceN = MotiveForceN;              //Total force applied to wheels
                 LocomotiveAxle.TrainSpeedMpS = SpeedMpS;                //Set the train speed of the axle mod
                 var watch = new Stopwatch();
                 watch.Start();

Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs

@@ -77,11 +77,7 @@ public class Axle
         /// <summary>
         /// Axle drive type covered by DriveType interface
         /// </summary>
-        protected AxleDriveType driveType;
-        /// <summary>
-        /// Read/Write Axle drive type flag
-        /// </summary>
-        public AxleDriveType DriveType { set { driveType = value; } get { return driveType; } }
+        public AxleDriveType DriveType;
 
         /// <summary>
         /// Axle drive represented by a motor, covered by ElectricMotor interface
@@ -96,17 +92,12 @@ public ElectricMotor Motor
             set
             {
                 motor = value;
-                switch(driveType)
+                DriveType = motor != null ? AxleDriveType.MotorDriven : AxleDriveType.ForceDriven;
+                switch(DriveType)
                 {
-                    case AxleDriveType.NotDriven:
-                        break;
                     case AxleDriveType.MotorDriven:
-                        //Total inertia considering gearbox
                         totalInertiaKgm2 = inertiaKgm2 + transmissionRatio * transmissionRatio * motor.InertiaKgm2;
                         break;
-                    case AxleDriveType.ForceDriven:
-                        totalInertiaKgm2 = inertiaKgm2;
-                        break;
                     default:
                         totalInertiaKgm2 = inertiaKgm2;
                         break;
@@ -149,7 +140,7 @@ public float InertiaKgm2
                 if (value <= 0.0)
                     throw new NotSupportedException("Inertia must be greater than zero");
                 inertiaKgm2 = value;
-                switch (driveType)
+                switch (DriveType)
                 {
                     case AxleDriveType.NotDriven:
                         break;
@@ -398,54 +389,8 @@ public Axle()
         {
             transmissionEfficiency = 0.99f;
             SlipWarningTresholdPercent = 70.0f;
-            driveType = AxleDriveType.ForceDriven;
-
-            switch (driveType)
-            {
-                case AxleDriveType.NotDriven:
-                    break;
-                case AxleDriveType.MotorDriven:
-                    totalInertiaKgm2 = inertiaKgm2 + transmissionRatio * transmissionRatio * motor.InertiaKgm2;
-                    break;
-                case AxleDriveType.ForceDriven:
-                    totalInertiaKgm2 = inertiaKgm2;
-                    break;
-                default:
-                    totalInertiaKgm2 = inertiaKgm2;
-                    break;
-            }
-        }
-
-        /// <summary>
-        /// Creates motor driven axle class instance
-        /// - sets TransmissionEfficiency to 0.99 (99%)
-        /// - sets SlipWarningThresholdPercent to 70%
-        /// - sets axle DriveType to MotorDriven
-        /// - updates totalInertiaKgm2 parameter
-        /// </summary>
-        /// <param name="electricMotor">Electric motor connected with the axle</param>
-        public Axle(ElectricMotor electricMotor)
-        {
-            motor = electricMotor;
-            motor.AxleConnected = this;
-            transmissionEfficiency = 0.99f;
-            driveType = AxleDriveType.MotorDriven;
-
-            switch (driveType)
-            {
-                case AxleDriveType.NotDriven:
-                    totalInertiaKgm2 = inertiaKgm2;
-                    break;
-                case AxleDriveType.MotorDriven:
-                    totalInertiaKgm2 = inertiaKgm2 + transmissionRatio * transmissionRatio * motor.InertiaKgm2;
-                    break;
-                case AxleDriveType.ForceDriven:
-                    totalInertiaKgm2 = inertiaKgm2;
-                    break;
-                default:
-                    totalInertiaKgm2 = inertiaKgm2;
-                    break;
-            }
+            DriveType = AxleDriveType.ForceDriven;
+            totalInertiaKgm2 = inertiaKgm2;
         }
 
         /// <summary>
@@ -478,10 +423,10 @@ public void Save(BinaryWriter outf)
             float axleForceN = AxleWeightN * SlipCharacteristics(axleSpeedMpS - TrainSpeedMpS, TrainSpeedMpS, AdhesionK, AdhesionLimit);
 
             float motiveAxleForceN = -axleForceN - dampingNs * (axleSpeedMpS - TrainSpeedMpS); // Force transmitted to rail + heat losses
-            if (driveType == AxleDriveType.ForceDriven)
-                motiveAxleForceN += driveForceN * transmissionEfficiency;
-            else if (driveType == AxleDriveType.MotorDriven)
-                motiveAxleForceN += motor.DevelopedTorqueNm * transmissionEfficiency * WheelRadiusM;
+            if (DriveType == AxleDriveType.ForceDriven)
+                motiveAxleForceN += DriveForceN * transmissionEfficiency;
+            else if (DriveType == AxleDriveType.MotorDriven)
+                motiveAxleForceN += motor.GetDevelopedTorqueNm(axleSpeedMpS * transmissionRatio / WheelRadiusM) * transmissionEfficiency / WheelRadiusM;
 
             // Dissipative forces: they will never increase wheel speed
             float frictionalForceN = BrakeRetardForceN + frictionN;
@@ -501,12 +446,17 @@ public void Save(BinaryWriter outf)
             return (totalAxleForceN * forceToAccelerationFactor, axleSpeedMpS / WheelRadiusM, axleForceN);
         }
 
+        /// <summary>
+        /// Integrates the wheel rotation movement using a RK4 method,
+        /// calculating the required number of substeps
+        /// Outputs: wheel speed, wheel angular position and motive force
+        /// </summary>
         void Integrate(float elapsedClockSeconds)
         {
             if (elapsedClockSeconds <= 0) return;
             float prevSpeedMpS = AxleSpeedMpS;
 
-            if (Math.Abs(integratorError) > Math.Max(SlipSpeedMpS, 1) / 1000)
+            if (Math.Abs(integratorError) > Math.Max((Math.Abs(SlipSpeedMpS)-1)*0.01f, 0.001f))
             {
                 ++NumOfSubstepsPS;
                 waitBeforeSpeedingUp = 100;
@@ -527,7 +477,7 @@ void Integrate(float elapsedClockSeconds)
             for (int i=0; i<NumOfSubstepsPS; i++)
             {
                 var k1 = GetAxleMotionVariation(AxleSpeedMpS);
-                if (i == 0 && k1.Item1 * dt > Math.Max(SlipSpeedMpS, 1) / 100)
+                if (i == 0 && k1.Item1 * dt > Math.Max((Math.Abs(SlipSpeedMpS) - 1) * 10, 1) / 100)
                 {
                     NumOfSubstepsPS = Math.Min(NumOfSubstepsPS + 5, 50);
                     dt = elapsedClockSeconds / NumOfSubstepsPS;
@@ -549,6 +499,37 @@ void Integrate(float elapsedClockSeconds)
             }
         }
 
+        /// <summary>
+        /// Work in progress. Calculates wheel creep assuming that wheel inertia
+        /// is low, removing the need of an integrator
+        /// Useful for slow CPUs
+        /// </summary>
+        void StationaryCalculation(float elapsedClockSeconds)
+        {
+            var res = GetAxleMotionVariation(AxleSpeedMpS);
+            float force = res.Item1 / forceToAccelerationFactor + res.Item3;
+            float maxAdhesiveForce = AxleWeightN * AdhesionLimit;
+            if (maxAdhesiveForce == 0) return;
+            float forceRatio = force / maxAdhesiveForce;
+            float absForceRatio = Math.Abs(forceRatio);
+            float characteristicTime = WheelSlipThresholdMpS / (maxAdhesiveForce * forceToAccelerationFactor);
+            if (absForceRatio > 1 || IsWheelSlip || Math.Abs(res.Item1 * elapsedClockSeconds) < WheelSlipThresholdMpS)
+            {
+                Integrate(elapsedClockSeconds);
+                return;
+            }
+            NumOfSubstepsPS = 1;
+            if (absForceRatio < 0.001f)
+            {
+                AxleSpeedMpS = TrainSpeedMpS;
+                AxleForceN = 0;
+                return;
+            }
+            float x = (float)((1 - Math.Sqrt(1 - forceRatio * forceRatio)) / forceRatio);
+            AxleSpeedMpS = TrainSpeedMpS + MpS.FromKpH(AdhesionK * x / AdhesionLimit);
+            AxleForceN = (force + res.Item3)/2;
+        }
+
         /// <summary>
         /// Main Update method
         /// - computes slip characteristics to get new axle force
@@ -570,11 +551,8 @@ public virtual void Update(float timeSpan)
             CompensatedAxleForceN = AxleForceN + Math.Sign(TrainSpeedMpS) * BrakeRetardForceN;
             if (AxleForceN == 0) CompensatedAxleForceN = 0;
 
-            if (driveType == AxleDriveType.MotorDriven)
-            {
-                motor.RevolutionsRad = AxleSpeedMpS * transmissionRatio / WheelRadiusM;
-                motor.Update(timeSpan);
-            }
+            motor?.Update(timeSpan);
+
             if (timeSpan > 0.0f)
             {
                 slipDerivationMpSS = (SlipSpeedMpS - previousSlipSpeedMpS) / timeSpan;
@@ -591,8 +569,7 @@ public virtual void Update(float timeSpan)
         public void Reset()
         {
             AxleSpeedMpS = 0;
-            if (motor != null)
-                motor.Reset();
+            motor?.Reset();
 
         }
 
@@ -604,8 +581,7 @@ public void Reset(double resetTime, float initValue)
         {
             AxleSpeedMpS = initValue;
             ResetTime = resetTime;
-            if (motor != null)
-                motor.Reset();
+            motor?.Reset();
         }
 
         /// <summary>
@@ -619,11 +595,12 @@ public void Reset(double resetTime, float initValue)
         ///             2*K*umax^2*dV
         ///     u = ---------------------
         ///           umax^2*dv^2 + K^2
-        /// For high slip speeds (after the inflexion point of u), the formula is
-        /// replaced with an exponentially decaying function (with smooth coupling)
-        /// reaching a 40% of maximum adhesion at infinity. Quick fix until
-        /// further investigation is done to get a single formula that provides
-        /// non zero adhesion at infinity.
+        /// 
+        /// For high slip speeds the formula is replaced with an exponentially 
+        /// decaying function (with smooth coupling) which reaches 40% of
+        /// maximum adhesion at infinity. The transition point between equations
+        /// is at dV = sqrt(3)*K/umax (inflection point)
+        /// 
         /// </summary>
         /// <param name="slipSpeedMpS">Difference between train speed and wheel speed</param>
         /// <param name="speedMpS">Current speed</param>
@@ -637,7 +614,7 @@ public float SlipCharacteristics(float slipSpeedMpS, float speedMpS, float K, fl
             float absx = Math.Abs(x);
             float sqrt3 = (float)Math.Sqrt(3);
             if (absx > sqrt3)
-            { 
+            {
                 // At infinity, adhesion is 40% of maximum (Polach, 2005)
                 // The value must be lower than 85% for the formula to work
                 float inftyFactor = 0.4f;