// EngineCylinder.cs (in Core namespace) using System; using FluidSim.Components; namespace FluidSim.Core { public class EngineCylinder { public Volume0D Cylinder { get; private set; } private Crankshaft crankshaft; private double bore, stroke, conRodLength, compressionRatio; private double pistonArea; private double V_disp, V_clear; private double maxOrificeArea; private double valveOpenStart = 120.0 * Math.PI / 180.0; private double valveOpenEnd = 480.0 * Math.PI / 180.0; private double valveRampWidth = 30.0 * Math.PI / 180.0; public double OrificeArea => ValveLift() * maxOrificeArea; public double TargetPeakPressure { get; set; } = 50.0 * 101325.0; private const double PeakTemperature = 2500.0; private bool burnInProgress = false; private double burnStartAngle; // full cycle angle when ignition began private double burnDuration = 40.0 * Math.PI / 180.0; private double targetBurnEnergy; private double totalBurnMass; private double preIgnitionMass, preIgnitionInternalEnergy; private Random rand = new Random(); public double MisfireProbability { get; set; } = 0.02; private bool misfireCurrent = false; public int CombustionCount { get; private set; } public int MisfireCount { get; private set; } public EngineCylinder(Crankshaft crankshaft, double bore, double stroke, double compressionRatio, double pipeArea, int sampleRate) { this.crankshaft = crankshaft; this.bore = bore; this.stroke = stroke; conRodLength = 2.0 * stroke; this.compressionRatio = compressionRatio; maxOrificeArea = pipeArea; pistonArea = Math.PI / 4.0 * bore * bore; V_disp = pistonArea * stroke; V_clear = V_disp / (compressionRatio - 1.0); // Initial compressed charge at TDC (no burn) double T_bdc = 300.0; double p_bdc = 101325.0; double V_bdc = V_clear + V_disp; double freshMass = p_bdc * V_bdc / (287.0 * T_bdc); double freshInternalEnergy = p_bdc * V_bdc / (1.4 - 1.0); double p_tdc = p_bdc * Math.Pow(V_bdc / V_clear, 1.4); Cylinder = new Volume0D(V_clear, p_tdc, T_bdc * Math.Pow(V_bdc / V_clear, 1.4 - 1.0), sampleRate) { Gamma = 1.4, GasConstant = 287.0 }; Cylinder.Volume = V_clear; Cylinder.Mass = freshMass; Cylinder.InternalEnergy = p_tdc * V_clear / (1.4 - 1.0); preIgnitionMass = Cylinder.Mass; preIgnitionInternalEnergy = Cylinder.InternalEnergy; } // ---- Piston kinematics (uses full cycle angle for position) ---- private (double volume, double dvdt) PistonKinematics(double cycleAngle) { // Slider-crank uses 0–2π, but we want the same motion for 0–2π (power/exhaust) and 2π–4π (intake/compression) double theta = cycleAngle % (2.0 * Math.PI); double R = stroke / 2.0; double cosT = Math.Cos(theta); double sinT = Math.Sin(theta); double L = conRodLength; double s = R * (1 - cosT) + L - Math.Sqrt(L * L - R * R * sinT * sinT); double V = V_clear + pistonArea * s; double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT); double dVdθ = pistonArea * (R * sinT + (R * R * sinT * cosT) / sqrtTerm); double dvdt = dVdθ * crankshaft.AngularVelocity; return (V, dvdt); } // ---- Valve lift ---- private double ValveLift() { double cycleRad = crankshaft.CrankAngle; if (cycleRad < valveOpenStart || cycleRad > valveOpenEnd) return 0.0; double duration = valveOpenEnd - valveOpenStart; double ramp = valveRampWidth; double t = (cycleRad - valveOpenStart) / duration; double rampFrac = ramp / duration; if (t < rampFrac) return t / rampFrac; else if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac; else return 1.0; } private double WiebeFraction(double angleFromIgnition) { if (angleFromIgnition >= burnDuration) return 1.0; double a = 5.0, m = 2.0; double x = angleFromIgnition / burnDuration; return 1.0 - Math.Exp(-a * Math.Pow(x, m + 1)); } // ---- Torque from pressure ---- private double ComputeTorque() { double p = Cylinder.Pressure; double ambient = 101325.0; double force = (p - ambient) * pistonArea; if (force <= 0) return 0.0; double theta = crankshaft.CrankAngle % (2.0 * Math.PI); double R = stroke / 2.0; double L = conRodLength; double sinT = Math.Sin(theta); double cosT = Math.Cos(theta); double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT); double lever = R * sinT * (1.0 + (R * cosT) / sqrtTerm); return force * lever; } public void Step(double dt) { double cycleAngle = crankshaft.CrankAngle; double prevAngle = crankshaft.PreviousAngle; // ----- TDC crossing detection (power stroke) ----- // Power stroke TDC occurs at angle 0 (mod 4π). We detect when PreviousAngle was near 4π and CrankAngle wraps to near 0. bool crossingTDC = (prevAngle > 3.8 * Math.PI && cycleAngle < 0.2 * Math.PI) // normal forward || (prevAngle < 0.2 * Math.PI && cycleAngle > 3.8 * Math.PI); // (rare backward, ignore) if (crossingTDC) { misfireCurrent = rand.NextDouble() < MisfireProbability; // Fresh charge: trapped at BDC, compressed isentropically to V_clear double T_bdc = 300.0; double p_bdc = 101325.0; double V_bdc = V_clear + V_disp; double freshMass = p_bdc * V_bdc / (287.0 * T_bdc); double freshInternalEnergy = p_bdc * V_bdc / (1.4 - 1.0); double gamma = 1.4; double p_tdc = p_bdc * Math.Pow(V_bdc / V_clear, gamma); Cylinder.Volume = V_clear; Cylinder.Mass = freshMass; Cylinder.InternalEnergy = p_tdc * V_clear / (gamma - 1.0); preIgnitionMass = Cylinder.Mass; preIgnitionInternalEnergy = Cylinder.InternalEnergy; if (misfireCurrent) { MisfireCount++; } else { double V = V_clear; targetBurnEnergy = TargetPeakPressure * V / (gamma - 1.0); totalBurnMass = TargetPeakPressure * V / (287.0 * PeakTemperature); burnInProgress = true; burnStartAngle = cycleAngle; CombustionCount++; } } // ----- Burn progress ----- if (burnInProgress) { double angleFromIgnition = cycleAngle - burnStartAngle; if (angleFromIgnition < 0) angleFromIgnition += 4.0 * Math.PI; // wrap if needed if (angleFromIgnition >= burnDuration) { Cylinder.Mass = totalBurnMass; Cylinder.InternalEnergy = targetBurnEnergy; burnInProgress = false; } else { double fraction = WiebeFraction(angleFromIgnition); Cylinder.InternalEnergy = preIgnitionInternalEnergy * (1.0 - fraction) + targetBurnEnergy * fraction; Cylinder.Mass = preIgnitionMass * (1.0 - fraction) + totalBurnMass * fraction; } } // ----- Piston motion ----- var (vol, dvdt) = PistonKinematics(cycleAngle); Cylinder.Volume = vol; Cylinder.Dvdt = dvdt; // ----- Torque contribution ----- double torque = ComputeTorque(); crankshaft.AddTorque(torque); } } }