FluidSim/Components/EngineCylinder.cs

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;

        public double V_disp { get; private set; }
        public double V_clear { get; private set; }
        public bool ignition = false;

        // ---- Exhaust valve ----
        private double exhMaxOrificeArea;
        private double exhValveOpenStart = 120.0 * Math.PI / 180.0;   // 120° (EVO)
        private double exhValveOpenEnd   = 480.0 * Math.PI / 180.0;   // 480° (EVC)
        private double exhValveRampWidth = 30.0  * Math.PI / 180.0;
        public double ExhaustOrificeArea => ExhaustValveLift() * exhMaxOrificeArea;

        // ---- Intake valve ----
        private double intMaxOrificeArea;
        private double intValveOpenStart = 380.0 * Math.PI / 180.0;   // 380° (IVO)
        private double intValveOpenEnd   = 560.0 * Math.PI / 180.0;   // 560° (IVC)
        private double intValveRampWidth = 30.0  * Math.PI / 180.0;
        public double IntakeOrificeArea => IntakeValveLift() * intMaxOrificeArea;

        // ---- Combustion ----
        public double TargetPeakPressure { get; set; } = 50.0 * 101325.0;
        private const double PeakTemperature = 2500.0;
        private bool burnInProgress = false;
        private double burnStartAngle;         // cycle angle (0–4π)
        private double burnDuration = 40.0 * Math.PI / 180.0;
        private double targetBurnEnergy;
        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; }

        // Cycle‑aware angle (0 – 4π)
        public double CycleAngle => crankshaft.CrankAngle % (4.0 * Math.PI);
        private double prevCycleAngle;

        // Piston position fraction (0 = TDC, 1 = BDC)
        public double PistonPositionFraction
        {
            get
            {
                double currentVol = Cylinder.Volume;
                if (currentVol <= V_clear) return 0.0;
                if (currentVol >= V_clear + V_disp) return 1.0;
                return (currentVol - V_clear) / V_disp;
            }
        }

        public EngineCylinder(Crankshaft crankshaft,
                              double bore, double stroke, double compressionRatio,
                              double exhPipeArea, double intPipeArea, int sampleRate)
        {
            this.crankshaft = crankshaft;
            this.bore = bore;
            this.stroke = stroke;
            conRodLength = 2.0 * stroke;
            this.compressionRatio = compressionRatio;
            exhMaxOrificeArea = exhPipeArea;
            intMaxOrificeArea = intPipeArea;
            pistonArea = Math.PI / 4.0 * bore * bore;

            V_disp = pistonArea * stroke;
            V_clear = V_disp / (compressionRatio - 1.0);

            // Start at BDC with fresh ambient charge
            double V_bdc = V_clear + V_disp;
            double p_amb = 101325.0;
            double T_amb = 300.0;
            double rho0 = p_amb / (287.0 * T_amb);
            double mass0 = rho0 * V_bdc;
            double energy0 = p_amb * V_bdc / (1.4 - 1.0);

            Cylinder = new Volume0D(V_bdc, p_amb, T_amb, sampleRate)
            {
                Gamma = 1.4,
                GasConstant = 287.0
            };
            Cylinder.Volume = V_bdc;
            Cylinder.Mass = mass0;
            Cylinder.InternalEnergy = energy0;

            prevCycleAngle = CycleAngle;

            preIgnitionMass = Cylinder.Mass;
            preIgnitionInternalEnergy = Cylinder.InternalEnergy;
        }

        // ---- Piston kinematics ----
        private (double volume, double dvdt) PistonKinematics(double cycleAngle)
        {
            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 lifts (cycle‑aware) ----
        private double ExhaustValveLift()
        {
            double a = CycleAngle;
            if (a < exhValveOpenStart || a > exhValveOpenEnd) return 0.0;
            double duration = exhValveOpenEnd - exhValveOpenStart;
            double ramp = exhValveRampWidth;
            double t = (a - exhValveOpenStart) / duration;
            double rampFrac = ramp / duration;
            if (t < rampFrac) return t / rampFrac;
            if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac;
            return 1.0;
        }

        private double IntakeValveLift()
        {
            double a = CycleAngle;
            if (a < intValveOpenStart || a > intValveOpenEnd) return 0.0;
            double duration = intValveOpenEnd - intValveOpenStart;
            double ramp = intValveRampWidth;
            double t = (a - intValveOpenStart) / duration;
            double rampFrac = ramp / duration;
            if (t < rampFrac) return t / rampFrac;
            if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac;
            return 1.0;
        }

        // ---- Wiebe burn fraction ----
        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;
        }

        // ---- TDC detection (power stroke, at angle 0 mod 4π) ----
        private bool DetectTDCPowerStroke()
        {
            double current = CycleAngle;
            double previous = prevCycleAngle;
            prevCycleAngle = current;
            return (previous > 3.8 * Math.PI && current < 0.2 * Math.PI);
        }

        public void Step(double dt)
        {
            bool crossingTDC = DetectTDCPowerStroke();

            if (crossingTDC)
            {
                misfireCurrent = rand.NextDouble() < MisfireProbability;

                // *** Always capture the state at TDC, whether we burn or not ***
                preIgnitionMass = Cylinder.Mass;
                preIgnitionInternalEnergy = Cylinder.InternalEnergy;

                if (misfireCurrent)
                {
                    MisfireCount++;
                }
                else if (ignition)
                {
                    double V = Cylinder.Volume;
                    targetBurnEnergy = TargetPeakPressure * V / (Cylinder.Gamma - 1.0);
                    if (double.IsNaN(targetBurnEnergy))
                        targetBurnEnergy = 101325.0 * V / (Cylinder.Gamma - 1.0);
                    burnInProgress = true;
                    burnStartAngle = CycleAngle;
                    CombustionCount++;
                }
            }

            if (burnInProgress)
            {
                double angleFromIgnition = CycleAngle - burnStartAngle;
                if (angleFromIgnition < 0) angleFromIgnition += 4.0 * Math.PI;

                if (angleFromIgnition >= burnDuration)
                {
                    Cylinder.InternalEnergy = targetBurnEnergy;
                    burnInProgress = false;
                }
                else
                {
                    double fraction = WiebeFraction(angleFromIgnition);
                    Cylinder.InternalEnergy = preIgnitionInternalEnergy * (1.0 - fraction)
                                              + targetBurnEnergy * fraction;
                    Cylinder.Mass = preIgnitionMass;
                }
            }

            var (vol, dvdt) = PistonKinematics(CycleAngle);
            Cylinder.Volume = vol;
            Cylinder.Dvdt = dvdt;

            if (double.IsNaN(Cylinder.Pressure) || double.IsNaN(Cylinder.Temperature) || Cylinder.Mass < 1e-9)
            {
                double V = Math.Max(vol, V_clear);
                Cylinder.Mass = 1.225 * V;
                Cylinder.InternalEnergy = 101325.0 * V / (1.4 - 1.0);
            }

            double torque = ComputeTorque();
            crankshaft.AddTorque(torque);
        }
    }
}