From f4c0b9133b6256925a4f9d624f233fc4342c20a7 Mon Sep 17 00:00:00 2001 From: peternewell Date: Tue, 26 Sep 2023 15:56:12 +1000 Subject: [PATCH 1/2] Implement Polach Adhesion --- .../RollingStocks/MSTSLocomotive.cs | 9 +- .../RollingStocks/MSTSSteamLocomotive.cs | 6 +- .../SubSystems/PowerTransmissions/Axle.cs | 375 ++++++++++++------ .../PowerTransmissions/ElectricMotor.cs | 2 +- .../PowerTransmissions/InductionMotor.cs | 4 +- .../PowerTransmissions/SeriesMotor.cs | 4 +- .../RunActivity/Viewer3D/Popups/HUDWindow.cs | 22 +- .../Viewer3D/RollingStock/MSTSWagonViewer.cs | 4 +- 8 files changed, 290 insertions(+), 136 deletions(-) diff --git a/Source/Orts.Simulation/Simulation/RollingStocks/MSTSLocomotive.cs b/Source/Orts.Simulation/Simulation/RollingStocks/MSTSLocomotive.cs index 5e0ec98b07..84affd9aa8 100644 --- a/Source/Orts.Simulation/Simulation/RollingStocks/MSTSLocomotive.cs +++ b/Source/Orts.Simulation/Simulation/RollingStocks/MSTSLocomotive.cs @@ -1469,7 +1469,7 @@ public override void Initialize() // Ensure Drive Axles is set with a default value if user doesn't supply an OR value in ENG file if (LocoNumDrvAxles == 0) { - if (MSTSLocoNumDrvWheels != 0 && MSTSLocoNumDrvWheels < 6) + if (MSTSLocoNumDrvWheels != 0 && MSTSLocoNumDrvWheels < 7) { LocoNumDrvAxles = (int) MSTSLocoNumDrvWheels; } @@ -2800,6 +2800,9 @@ public virtual void AdvancedAdhesion(float elapsedClockSeconds) axle.BrakeRetardForceN = BrakeRetardForceN/LocomotiveAxles.Count; axle.TrainSpeedMpS = SpeedMpS; //Set the train speed of the axle mod axle.WheelRadiusM = DriverWheelRadiusM; + axle.WheelDistanceGaugeM = TrackGaugeM; + axle.CurrentCurveRadiusM = CurrentCurveRadiusM; + axle.BogieRigidWheelBaseM = RigidWheelBaseM; } LocomotiveAxles.Update(elapsedClockSeconds); MotiveForceN = LocomotiveAxles.CompensatedForceN; @@ -2813,10 +2816,10 @@ public virtual void AdvancedAdhesion(float elapsedClockSeconds) // This enables steam locomotives to have different speeds for driven and non-driven wheels. if (EngineType == EngineTypes.Steam && SteamEngineType != MSTSSteamLocomotive.SteamEngineTypes.Geared) { - WheelSpeedSlipMpS = LocomotiveAxles[0].AxleSpeedMpS; + WheelSpeedSlipMpS = (float)LocomotiveAxles[0].AxleSpeedMpS; WheelSpeedMpS = SpeedMpS; } - else WheelSpeedMpS = LocomotiveAxles[0].AxleSpeedMpS; + else WheelSpeedMpS = (float)LocomotiveAxles[0].AxleSpeedMpS; } diff --git a/Source/Orts.Simulation/Simulation/RollingStocks/MSTSSteamLocomotive.cs b/Source/Orts.Simulation/Simulation/RollingStocks/MSTSSteamLocomotive.cs index 88106c0ab3..e5359cc21f 100644 --- a/Source/Orts.Simulation/Simulation/RollingStocks/MSTSSteamLocomotive.cs +++ b/Source/Orts.Simulation/Simulation/RollingStocks/MSTSSteamLocomotive.cs @@ -5045,7 +5045,7 @@ protected override void UpdateTractiveForce(float elapsedClockSeconds, float t, // Crank Angle - converts the above range to 0 - 180 - 0 - this is the principle reference used so that it lines up with reference // tables used to buold this function // Normalised Crank Angle - converts the above to a 0 - 360 range, this is used for triggering special steam effects, etc. - float axlePostionRad = LocomotiveAxles[0].AxlePositionRad; + float axlePostionRad = (float)LocomotiveAxles[0].AxlePositionRad; float crankAngleRad = (float)(axlePostionRad + WheelCrankAngleDiffRad[i]); crankAngleRad = (float)(MathHelper.WrapAngle(crankAngleRad)); // Ensures that crank angle is in the range 0 - 180 - 0 @@ -5299,8 +5299,8 @@ protected override void UpdateTractiveForce(float elapsedClockSeconds, float t, /// private float NormalisedCrankAngle(int cylinderNumber) { - float normalisedCrankAngleRad = (float)MathHelper.WrapAngle(LocomotiveAxles[0].AxlePositionRad + WheelCrankAngleDiffRad[cylinderNumber]); - + float normalisedCrankAngleRad = (float)MathHelper.WrapAngle((float)LocomotiveAxles[0].AxlePositionRad + WheelCrankAngleDiffRad[cylinderNumber]); + if (normalisedCrankAngleRad < 0) { normalisedCrankAngleRad += (float)(2 * Math.PI); diff --git a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs index b7540a132e..01692743df 100644 --- a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs +++ b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs @@ -25,6 +25,7 @@ using Orts.Parsers.Msts; using Orts.Simulation.RollingStocks.SubSystems.PowerTransmissions; using SharpDX.Direct2D1; +using SharpDX.Direct3D9; namespace Orts.Simulation.RollingStocks.SubSystems.PowerTransmissions { @@ -321,6 +322,8 @@ public class Axle : ISubSystem /// public float DampingNs { set { dampingNs = Math.Abs(value); } get { return dampingNs; } } + int count; + protected float frictionN; public float FrictionN { set { frictionN = Math.Abs(value); } get { return frictionN; } } @@ -417,6 +420,11 @@ public float InertiaKgm2 /// float forceToAccelerationFactor; + /// + /// Pre-calculation of slip characteristics at 0 slip speed + /// + double axleStaticForceN; + /// /// Transmission ratio on gearbox covered by TransmissionRatio interface /// @@ -470,11 +478,36 @@ public float TransmissionEfficiency /// public float WheelRadiusM; + /// + /// Gauge of Track + /// + public float WheelDistanceGaugeM; + + /// + /// Radius of Track Curve + /// + public float CurrentCurveRadiusM; + + /// + /// Bogie Rigid Wheel Base - distance between wheel in the bogie + /// + public float BogieRigidWheelBaseM; + + /// + /// Axles in group of wheels + /// + public float NumAxles = 1; + /// /// Static adhesion coefficient, as given by Curtius-Kniffler formula /// public float AdhesionLimit; + /// + /// Wheel adhesion as calculated by Polach + /// + public float WheelAdhesion; + /// /// Correction parameter of adhesion, it has proportional impact on adhesion limit /// Should be set to 1.0 for most cases @@ -484,11 +517,13 @@ public float TransmissionEfficiency /// /// Axle speed value, in metric meters per second /// - public float AxleSpeedMpS { get; private set; } + public double AxleSpeedMpS { get; private set; } + /// /// Axle angular position in radians /// - public float AxlePositionRad { get; private set; } + public double AxlePositionRad { get; private set; } + /// /// Read only axle force value, in Newtons /// @@ -517,19 +552,44 @@ public float TransmissionEfficiency float WheelSlipTimeS; /// - /// Read only wheelslip threshold value used to indicate maximal effective slip - /// - its value is computed as a maximum of slip function: - /// 2*K*umax^2 * dV - /// f(dV) = u = --------------------- - /// umax^2*dV^2 + K^2 + /// Wheelslip threshold value used to indicate maximal effective slip + /// - its value is computed as a maximum of slip function /// maximum can be found as a derivation f'(dV) = 0 /// - public float WheelSlipThresholdMpS + public float WheelSlipThresholdMpS; + public void ComputeWheelSlipThresholdMpS() { - get + // Bisection algorithm. We assume adhesion maximum is between 0 and 4 m/s + double a = 0; + double b = 4; + // We have to find the zero of the derivative of adhesion curve + // i.e. the point where slope changes from positive (adhesion region) + // to negative (slip region) + double dx = 0.001; + double fa = SlipCharacteristics(a + dx) - SlipCharacteristics(a); + double fb = SlipCharacteristics(b + dx) - SlipCharacteristics(b); + if (fa * fb > 0) + { + // If sign does not change, bisection fails + WheelSlipThresholdMpS = MpS.FromKpH(0.2f); + return; + } + while (Math.Abs(b - a) > MpS.FromKpH(0.05f)) { - return MpS.FromKpH(AdhesionK / AdhesionLimit); + double c = (a + b) / 2; + double fc = SlipCharacteristics(c + dx) - SlipCharacteristics(c); + if (fa * fc > 0) + { + a = c; + fa = fc; + } + else + { + b = c; + fb = fc; + } } + WheelSlipThresholdMpS = (float)Math.Max((a + b) / 2, MpS.FromKpH(0.2f)); } /// @@ -549,7 +609,7 @@ public float SlipSpeedMpS { get { - return (AxleSpeedMpS - TrainSpeedMpS); + return (float)(AxleSpeedMpS - TrainSpeedMpS); } } @@ -605,7 +665,7 @@ public float SlipDerivationPercentpS } } - float integratorError; + double integratorError; int waitBeforeSpeedingUp; /// @@ -613,6 +673,8 @@ public float SlipDerivationPercentpS /// public float SlipWarningTresholdPercent { set; get; } + PolachCalculator Polach; + public List AnimatedParts = new List(); /// @@ -629,6 +691,7 @@ public Axle() SlipWarningTresholdPercent = 70.0f; DriveType = AxleDriveType.ForceDriven; totalInertiaKgm2 = inertiaKgm2; + Polach = new PolachCalculator(this); } public void Initialize() { @@ -649,6 +712,9 @@ public void Parse(STFReader stf) case "ortsradius": WheelRadiusM = stf.ReadFloatBlock(STFReader.UNITS.Distance, null); break; + case "ortsnumberwheelaxles": + NumAxles = stf.ReadFloatBlock(STFReader.UNITS.Distance, null); + break; case "ortsinertia": InertiaKgm2 = stf.ReadFloatBlock(STFReader.UNITS.RotationalInertia, null); break; @@ -670,6 +736,7 @@ public void Parse(STFReader stf) public void Copy(Axle other) { WheelRadiusM = other.WheelRadiusM; + NumAxles = other.NumAxles; InertiaKgm2 = other.InertiaKgm2; AxleWeightN = other.AxleWeightN; AnimatedParts.Clear(); @@ -707,19 +774,31 @@ public void Save(BinaryWriter outf) /// /// Compute variation in axle dynamics. Calculates axle speed, axle angular position and in/out forces. /// - public (float, float, float, float) GetAxleMotionVariation(float axleSpeedMpS) + public (double, double, double, double) GetAxleMotionVariation(double axleSpeedMpS, double elapsedClockSeconds) { - float axleOutForceN = AxleWeightN * SlipCharacteristics(axleSpeedMpS - TrainSpeedMpS, TrainSpeedMpS, AdhesionK, AdhesionLimit); - - float axleInForceN = 0; + double slipSpeedMpS = axleSpeedMpS - TrainSpeedMpS; + double axleOutForceN = Math.Sign(slipSpeedMpS) * AxleWeightN * SlipCharacteristics(slipSpeedMpS); + double axleInForceN = 0; if (DriveType == AxleDriveType.ForceDriven) axleInForceN = DriveForceN * transmissionEfficiency; else if (DriveType == AxleDriveType.MotorDriven) axleInForceN = motor.GetDevelopedTorqueNm(axleSpeedMpS * transmissionRatio / WheelRadiusM) * transmissionEfficiency / WheelRadiusM; - float totalAxleForceN = axleInForceN - axleOutForceN - dampingNs * (axleSpeedMpS - TrainSpeedMpS); // Force transmitted to rail + heat losses - - totalAxleForceN -= Math.Sign(axleSpeedMpS) * (BrakeRetardForceN + frictionN); // Dissipative forces: they will never increase wheel speed + double motionForceN = axleInForceN - dampingNs * (axleSpeedMpS - TrainSpeedMpS); // Drive force + heat losses + double frictionForceN = BrakeRetardForceN + frictionN; // Dissipative forces: they will never increase wheel speed + double totalAxleForceN = motionForceN - Math.Sign(axleSpeedMpS) * frictionForceN; + if (Math.Abs(TrainSpeedMpS) < 0.001f && Math.Abs(slipSpeedMpS) < 0.001f && Math.Abs(motionForceN) < frictionForceN) + { + return (-slipSpeedMpS / elapsedClockSeconds, axleSpeedMpS / WheelRadiusM, 0, axleInForceN); + } + if (Math.Abs(totalAxleForceN) < axleStaticForceN) + { + if (Math.Abs(slipSpeedMpS) < 0.001f || Math.Sign(slipSpeedMpS) != Math.Sign(slipSpeedMpS + (totalAxleForceN - axleOutForceN) * forceToAccelerationFactor * elapsedClockSeconds)) + { + axleOutForceN = slipSpeedMpS / elapsedClockSeconds / forceToAccelerationFactor + totalAxleForceN; + } + } + totalAxleForceN -= axleOutForceN; return (totalAxleForceN * forceToAccelerationFactor, axleSpeedMpS / WheelRadiusM, axleOutForceN, axleInForceN); } @@ -732,9 +811,9 @@ public void Save(BinaryWriter outf) void Integrate(float elapsedClockSeconds) { if (elapsedClockSeconds <= 0) return; - float prevSpeedMpS = AxleSpeedMpS; + double prevSpeedMpS = AxleSpeedMpS; - if (Math.Abs(integratorError) > Math.Max((Math.Abs(SlipSpeedMpS)-1)*0.01f, 0.001f)) + if (Math.Abs(integratorError) > Math.Max((Math.Abs(SlipSpeedMpS) - 1) * 0.01, 0.001)) { ++NumOfSubstepsPS; waitBeforeSpeedingUp = 100; @@ -748,15 +827,15 @@ void Integrate(float elapsedClockSeconds) } } - NumOfSubstepsPS = Math.Max(Math.Min(NumOfSubstepsPS, 50), 1); - - float dt = elapsedClockSeconds / NumOfSubstepsPS; - float hdt = dt / 2.0f; - float axleInForceSumN = 0; - float axleOutForceSumN = 0; - for (int i=0; i Math.Max((Math.Abs(SlipSpeedMpS) - 1) * 10, 1) / 100) @@ -764,60 +843,20 @@ void Integrate(float elapsedClockSeconds) NumOfSubstepsPS = Math.Min(NumOfSubstepsPS + 5, 50); dt = elapsedClockSeconds / NumOfSubstepsPS; hdt = dt / 2; - } - - if (Math.Sign(AxleSpeedMpS + k1.Item1 * dt) != Math.Sign(AxleSpeedMpS) && BrakeRetardForceN + frictionN > Math.Abs(driveForceN - k1.Item3)) - { - AxlePositionRad += AxleSpeedMpS * hdt; - AxlePositionRad = MathHelper.WrapAngle(AxlePositionRad); - AxleSpeedMpS = 0; - AxleForceN = 0; - DriveForceN = k1.Item4; - return; - } + } } - var k2 = GetAxleMotionVariation(AxleSpeedMpS + k1.Item1 * hdt); - var k3 = GetAxleMotionVariation(AxleSpeedMpS + k2.Item1 * hdt); - var k4 = GetAxleMotionVariation(AxleSpeedMpS + k3.Item1 * dt); - AxleSpeedMpS += (integratorError = (k1.Item1 + 2 * (k2.Item1 + k3.Item1) + k4.Item1) * dt / 6.0f); - AxlePositionRad += (k1.Item2 + 2 * (k2.Item2 + k3.Item2) + k4.Item2) * dt / 6.0f; + var k2 = GetAxleMotionVariation(AxleSpeedMpS + k1.Item1 * hdt, hdt); + var k3 = GetAxleMotionVariation(AxleSpeedMpS + k2.Item1 * hdt, hdt); + var k4 = GetAxleMotionVariation(AxleSpeedMpS + k3.Item1 * dt, dt); + + AxleSpeedMpS += (integratorError = (k1.Item1 + 2 * (k2.Item1 + k3.Item1) + k4.Item1) * dt / 6); + AxlePositionRad += (k1.Item2 + 2 * (k2.Item2 + k3.Item2) + k4.Item2) * dt / 6; axleOutForceSumN += (k1.Item3 + 2 * (k2.Item3 + k3.Item3) + k4.Item3); axleInForceSumN += (k1.Item4 + 2 * (k2.Item4 + k3.Item4) + k4.Item4); } - AxleForceN = axleOutForceSumN / (NumOfSubstepsPS * 6); - DriveForceN = axleInForceSumN / (NumOfSubstepsPS * 6); - AxlePositionRad = MathHelper.WrapAngle(AxlePositionRad); - } - - /// - /// Work in progress. Calculates wheel creep assuming that wheel inertia - /// is low, removing the need of an integrator - /// Useful for slow CPUs - /// - 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; + AxleForceN = (float)(axleOutForceSumN / (NumOfSubstepsPS * 6)); + DriveForceN = (float)(axleInForceSumN / (NumOfSubstepsPS * 6)); + AxlePositionRad = MathHelper.WrapAngle((float)AxlePositionRad); } /// @@ -831,6 +870,40 @@ public virtual void Update(float timeSpan) { forceToAccelerationFactor = WheelRadiusM * WheelRadiusM / totalInertiaKgm2; + Polach.Update(); + axleStaticForceN = AxleWeightN * SlipCharacteristics(0); + ComputeWheelSlipThresholdMpS(); + + if (count < 6 && count++ == 5) + { + TrainSpeedMpS = 10 / 3.6f; + Polach.Update(); + axleStaticForceN = AxleWeightN * SlipCharacteristics(0); + + /* + double[] spd = new double[50]; + double[] adh = new double[50]; + for (int i = 0; i < spd.Length; i++) + { + spd[i] = i / (float)spd.Length; + adh[i] = SlipCharacteristics(spd[i]); + } + for (int i = 0; i < spd.Length; i++) + { + Console.Write(spd[i]); + Console.Write(" "); + } + Console.WriteLine(""); + Console.WriteLine(""); + for (int i = 0; i < spd.Length; i++) + { + Console.Write(adh[i]); + Console.Write(" "); + } + Console.WriteLine(""); + */ + } + motor?.Update(timeSpan); Integrate(timeSpan); @@ -840,10 +913,10 @@ public virtual void Update(float timeSpan) // The Axle module subtracts brake force from the motive force for calculation purposes. However brake force is already taken into account in the braking module. // And thus there is a duplication of the braking effect in OR. To compensate for this, after the slip characteristics have been calculated, the output of the axle module // has the brake force "added" back in to give the appropriate motive force output for the locomotive. Braking force is handled separately. - // Hence CompensatedAxleForce is the actual output force on the axle. - CompensatedAxleForceN = AxleForceN + Math.Sign(TrainSpeedMpS) * BrakeRetardForceN; - if (AxleForceN == 0) CompensatedAxleForceN = 0; - + // Hence CompensatedAxleForce is the actual output force on the axle. + if (Math.Abs(TrainSpeedMpS) < 0.001f && AxleForceN == 0) CompensatedAxleForceN = DriveForceN; + else if (TrainSpeedMpS < 0) CompensatedAxleForceN = AxleForceN - BrakeRetardForceN; + else CompensatedAxleForceN = AxleForceN + BrakeRetardForceN; if (Math.Abs(SlipSpeedMpS) > WheelSlipThresholdMpS) { // Wait some time before indicating wheelslip to avoid false triggers @@ -877,43 +950,121 @@ public virtual void Update(float timeSpan) } } + class PolachCalculator + { + Axle Axle; + double StiffnessPreFactors; + double polach_A; + double polach_B; + double polach_Ka; + double polach_Ks; + double umax; + double speedMpS; + double spinM1; + double Sy = 0; + double Sy2 = 0; + double Syc2 = 0; + double kalker_C11; + double kalker_C22; + public PolachCalculator(Axle axle) + { + Axle = axle; + } + public void Update() + { + umax = Axle.AdhesionLimit; + speedMpS = Math.Abs(Axle.TrainSpeedMpS); + + var wheelRadiusMM = Axle.WheelRadiusM * 1000; + var wheelDistanceGaugeMM = Axle.WheelDistanceGaugeM * 1000; + var GNm2 = 8.40E+10; + var wheelLoadN = Axle.AxleWeightN / (Axle.NumAxles * 2); // Assume two wheels per axle, and thus wheel weight will be have the value - multiple axles???? + var wheelLoadkN = Axle.AxleWeightN / (Axle.NumAxles * 2 * 1000); // Assume two wheels per axle, and thus wheel weight will be have the value - multiple axles???? + var Young_ModulusMPa = 207000; + + // Calculate Hertzian values - assume 2b = 12mm. + var b_HertzianMM = 6.0; + var a_HertzianMM = (3.04f / 2) * Math.Sqrt(((wheelLoadkN * wheelRadiusMM) / (2 * b_HertzianMM * Young_ModulusMPa))); + + var hertzianMM = a_HertzianMM / b_HertzianMM; + var hertzianMM2 = hertzianMM * hertzianMM; + // Calculate Kalker values + kalker_C11 = 0.32955 * hertzianMM2 + 0.48538 * hertzianMM + 3.4992; + kalker_C22 = 0.90909 * hertzianMM2 + 1.2594 * hertzianMM + 2.3853; + + StiffnessPreFactors = (GNm2 * Math.PI * a_HertzianMM * b_HertzianMM) / (4.0 * wheelLoadN); + + // Calculate slip and creep values + var wheelProfileConicityRad = 0.5; + var wheelCentreDeviationMM = 3.0; + double YawAngleRad = 0; + if (Axle.CurrentCurveRadiusM > 0) + { + YawAngleRad = Math.Asin(2.0f * Axle.BogieRigidWheelBaseM / Axle.CurrentCurveRadiusM); + } + var lateralSlipVelocityMpS = speedMpS * (Axle.CurrentCurveRadiusM * wheelProfileConicityRad * wheelCentreDeviationMM * YawAngleRad) / (wheelDistanceGaugeMM * wheelRadiusMM); + spinM1 = Math.Sin(wheelProfileConicityRad) / wheelRadiusMM; + Sy = lateralSlipVelocityMpS / speedMpS; + Sy2 = Sy * Sy; + var Syc = Sy + spinM1 * a_HertzianMM; + Syc2 = Syc * Syc; + + // calculate "standard" Polach adhesion parameters as straight line approximations as u varies + polach_A = 0.4; + polach_B = (1.6 * umax) - 0.28; + polach_Ka = (2.8 * umax) - 0.54; + polach_Ks = (1.2 * umax) - 0.26; + } + public double SlipCharacteristics(double slipSpeedMpS) + { + var polach_uadhesion = umax * ((1 - polach_A) * Math.Exp(-polach_B * slipSpeedMpS) + polach_A); + + if (speedMpS < 0.05f) + return polach_uadhesion; + + var Sx = slipSpeedMpS / speedMpS; + var Sx2 = Sx * Sx; + var S = Math.Sqrt(Sx2 + Sy2); + var Sc = Math.Sqrt(Sx2 + Syc2); + + var kalker_Cjj = 1.0; + if (S != 0) // prevents NaN calculation if slip values are zero + { + var coef1 = kalker_C11 * Sx / S; + var coef2 = kalker_C22 * Sy / S; + kalker_Cjj = Math.Sqrt(coef1 * coef1 + coef2 * coef2); + } + + var Stiffness = StiffnessPreFactors * kalker_Cjj * Sc / polach_uadhesion; + var KaStiffness = polach_Ka * Stiffness; + var adhesionComponent = KaStiffness / (1 + KaStiffness * KaStiffness); + var slipComponent = Math.Atan(polach_Ks * Stiffness); + var f = polach_uadhesion / (Math.PI / 2) * (adhesionComponent + slipComponent); + var fx = f * Sx / Sc; + return fx; + } + } + /// - /// Slip characteristics computation - /// - Uses adhesion limit calculated by Curtius-Kniffler formula: - /// 7.5 - /// umax = --------------------- + 0.161 - /// speed * 3.6 + 44.0 - /// - Computes slip speed - /// - Computes relative adhesion force as a result of slip characteristics: - /// 2*K*umax^2*dV - /// u = --------------------- - /// umax^2*dv^2 + K^2 + /// Uses the Polach creep force curves calculation described in the following document + /// "Creep forces in simulations of traction vehicles running on adhesion limit" by O. Polach 2005 Wear - http://www.sze.hu/~szenasy/VILLVONT/polachslipvizsg.pdf + /// + /// To calculate the Polach curves the Hertzian, and Kalker values also need to be calculated as inputs into the Polach formula. + /// The following paper describes a methodology for calculating these values - /// - /// 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) + /// "Determination of Wheel-Rail Contact Characteristics by Creating a Special Program for Calculation" by Ioan Sebesan and Yahia Zakaria + /// - https://www.researchgate.net/publication/273303594_Determination_of_Wheel-Rail_Contact_Characteristics_by_Creating_a_Special_Program_for_Calculation /// /// /// Difference between train speed and wheel speed /// Current speed - /// Slip speed correction - /// Relative weather conditions, usually from 0.2 to 1.0 /// Relative force transmitted to the rail - public float SlipCharacteristics(float slipSpeedMpS, float speedMpS, float K, float umax) + public double SlipCharacteristics(double slipSpeedMpS) { - var slipSpeedKpH = MpS.ToKpH(slipSpeedMpS); - float x = slipSpeedKpH * umax / K; // Slip percentage - 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; - return Math.Sign(slipSpeedKpH) * umax * ((sqrt3 / 2 - inftyFactor) * (float)Math.Exp((sqrt3 - absx) / (2 * sqrt3 - 4 * inftyFactor)) + inftyFactor); - } - return 2.0f * umax * x / (1 + x * x); + slipSpeedMpS = Math.Abs(slipSpeedMpS); + double fx = Polach.SlipCharacteristics(slipSpeedMpS); + WheelAdhesion = (float)fx; + return fx; } } } diff --git a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/ElectricMotor.cs b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/ElectricMotor.cs index 594e43a1e9..65266a3e3f 100644 --- a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/ElectricMotor.cs +++ b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/ElectricMotor.cs @@ -61,7 +61,7 @@ public virtual void UpdateTractiveForce(float elapsedClockSeconds, float t) } - public virtual float GetDevelopedTorqueNm(float motorSpeedRadpS) + public virtual double GetDevelopedTorqueNm(double motorSpeedRadpS) { return 0; } diff --git a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/InductionMotor.cs b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/InductionMotor.cs index f24e84cb09..4377d81d3d 100644 --- a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/InductionMotor.cs +++ b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/InductionMotor.cs @@ -48,9 +48,9 @@ public override void UpdateTractiveForce(float elapsedClockSeconds, float t) { TargetForceN = Locomotive.TractiveForceN / Locomotive.LocomotiveAxles.Count; } - public override float GetDevelopedTorqueNm(float motorSpeedRadpS) + public override double GetDevelopedTorqueNm(double motorSpeedRadpS) { - return requiredTorqueNm * MathHelper.Clamp((DriveSpeedRadpS - motorSpeedRadpS) / OptimalAsyncSpeedRadpS, -1, 1); + return requiredTorqueNm * MathHelper.Clamp((float)(DriveSpeedRadpS - motorSpeedRadpS) / OptimalAsyncSpeedRadpS, -1, 1); } public override void Update(float timeSpan) { diff --git a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/SeriesMotor.cs b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/SeriesMotor.cs index f3ce96b850..6cc0f1e002 100644 --- a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/SeriesMotor.cs +++ b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/SeriesMotor.cs @@ -136,9 +136,9 @@ public SeriesMotor(float nomCurrentA, float nomVoltageV, float nomRevolutionsRad NominalRevolutionsRad = nomRevolutionsRad; } - public override float GetDevelopedTorqueNm(float revolutionsRad) + public override double GetDevelopedTorqueNm(double revolutionsRad) { - BackEMFvoltageV = revolutionsRad * fieldWb; + BackEMFvoltageV = (float)revolutionsRad * fieldWb; return fieldWb * armatureCurrentA/* - (frictionTorqueNm * revolutionsRad / NominalRevolutionsRad * revolutionsRad / NominalRevolutionsRad)*/; } diff --git a/Source/RunActivity/Viewer3D/Popups/HUDWindow.cs b/Source/RunActivity/Viewer3D/Popups/HUDWindow.cs index 53094d1588..6ada2f75c5 100644 --- a/Source/RunActivity/Viewer3D/Popups/HUDWindow.cs +++ b/Source/RunActivity/Viewer3D/Popups/HUDWindow.cs @@ -1085,26 +1085,26 @@ void TextPageForceInfo(TableData table) TableSetCell(table, table.CurrentRow++, table.CurrentLabelColumn, Viewer.Catalog.GetString("Axle drive force")); TableSetCell(table, table.CurrentRow++, table.CurrentLabelColumn, Viewer.Catalog.GetString("Axle brake force")); TableSetCell(table, table.CurrentRow++, table.CurrentLabelColumn, Viewer.Catalog.GetString("Number of substeps")); + TableSetCell(table, table.CurrentRow++, table.CurrentLabelColumn, Viewer.Catalog.GetString("Wheel Adhesion")); TableSetCell(table, table.CurrentRow++, table.CurrentLabelColumn, Viewer.Catalog.GetString("Axle out force")); TableSetCell(table, table.CurrentRow++, table.CurrentLabelColumn, Viewer.Catalog.GetString("Comp Axle out force")); TableSetCell(table, table.CurrentRow++, table.CurrentLabelColumn, Viewer.Catalog.GetString("Wheel speed")); - for (int i=0; i= 0 && kvp.Key < loco.LocomotiveAxles.Count ? loco.LocomotiveAxles[kvp.Key] : null; if (axle != null) //TODO: next code line has been modified to flip trainset physics in order to get viewing direction coincident with loco direction when using rear cab. - kvp.Value.UpdateLoop(((MSTSWagon.Train != null && MSTSWagon.Train.IsPlayerDriven && loco.UsingRearCab) ? -1 : 1) * axle.AxleSpeedMpS * elapsedTime.ClockSeconds / MathHelper.TwoPi / axle.WheelRadiusM); + kvp.Value.UpdateLoop(((MSTSWagon.Train != null && MSTSWagon.Train.IsPlayerDriven && loco.UsingRearCab) ? -1 : 1) * (float)axle.AxleSpeedMpS * elapsedTime.ClockSeconds / MathHelper.TwoPi / axle.WheelRadiusM); else if (AnimationDriveWheelRadiusM > 0.001) kvp.Value.UpdateLoop(distanceTravelledM / MathHelper.TwoPi / AnimationDriveWheelRadiusM); } @@ -794,7 +794,7 @@ private void UpdateAnimation(RenderFrame frame, ElapsedTime elapsedTime) if (axle != null) { //TODO: next code line has been modified to flip trainset physics in order to get viewing direction coincident with loco direction when using rear cab. - wheelRotationMatrix = Matrix.CreateRotationX(((MSTSWagon.Train != null && MSTSWagon.Train.IsPlayerDriven && loco.UsingRearCab) ? -1 : 1) * -axle.AxlePositionRad); + wheelRotationMatrix = Matrix.CreateRotationX(((MSTSWagon.Train != null && MSTSWagon.Train.IsPlayerDriven && loco.UsingRearCab) ? -1 : 1) * -(float)axle.AxlePositionRad); } else { From 0b44e34871dd1d7b2b742a2749af7bf66814f1d1 Mon Sep 17 00:00:00 2001 From: peternewell Date: Wed, 27 Sep 2023 15:33:19 +1000 Subject: [PATCH 2/2] Correct -ve occuring issues --- .../SubSystems/PowerTransmissions/Axle.cs | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs index 01692743df..4ec3dab792 100644 --- a/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs +++ b/Source/Orts.Simulation/Simulation/RollingStocks/SubSystems/PowerTransmissions/Axle.cs @@ -623,7 +623,7 @@ public float SlipSpeedPercent { var temp = SlipSpeedMpS / WheelSlipThresholdMpS * 100.0f; if (float.IsNaN(temp)) temp = 0;//avoid NaN on HuD display when first starting OR - return temp; + return Math.Abs(temp); } } @@ -1009,11 +1009,17 @@ public void Update() var Syc = Sy + spinM1 * a_HertzianMM; Syc2 = Syc * Syc; - // calculate "standard" Polach adhesion parameters as straight line approximations as u varies + // calculate "standard" Polach adhesion parameters as straight line approximations as u varies - these values are capped at the moment at the u=0.3 level + // Taking them lower may reduce the stability of the calculations polach_A = 0.4; polach_B = (1.6 * umax) - 0.28; + if (polach_B < 0.2) polach_B = 0.2f; + polach_Ka = (2.8 * umax) - 0.54; + if (polach_Ka < 0.3) polach_Ka = 0.3f; + polach_Ks = (1.2 * umax) - 0.26; + if (polach_Ks < 0.1) polach_Ks = 0.1f; } public double SlipCharacteristics(double slipSpeedMpS) { @@ -1041,6 +1047,14 @@ public double SlipCharacteristics(double slipSpeedMpS) var slipComponent = Math.Atan(polach_Ks * Stiffness); var f = polach_uadhesion / (Math.PI / 2) * (adhesionComponent + slipComponent); var fx = f * Sx / Sc; +/* + if (fx < 0) + { + Trace.TraceInformation("Negative Fx - Fx {0} Speed {1} SlipSpeed {2} Sx {3} PolachAdhesion {4} adhesionComponent {5} slipComponent {6} Polach_Ks {7}", fx, MpS.ToMpH((float)speedMpS), MpS.ToMpH((float)slipSpeedMpS), MpS.ToMpH((float)Sx), polach_uadhesion, adhesionComponent, slipComponent, polach_Ks); + + } +*/ + return fx; } }