FluidSim/Components/Crankcase.cs

// ============================================================
// File: Crankcase.cs
// ============================================================
using System.Collections.Generic;
using FluidSim.Interfaces;

namespace FluidSim.Components
{
    public class Crankcase : IComponent
    {
        private readonly Crankshaft _crankshaft;
        private readonly float _crankRadius, _conrodLength, _pistonArea;
        private readonly float _clearanceVolume, _obliquity;

        private float _mass, _internalEnergy, _airFraction;
        public float Pressure { get; private set; }
        public float Temperature { get; private set; }
        public float Density => _mass / MathF.Max(Volume, 1e-12f);
        public float Volume { get; private set; }

        public Port IntakePort { get; }
        public Port TransferPort { get; }
        private readonly List<Port> _ports;
        public IReadOnlyList<Port> Ports => _ports;

        // FIX: store previous volume to calculate PdV work
        private float _prevVolume;

        private const float Rgas = 287.0f;
        private const float Gamma = 1.4f;
        private const float Cv = Rgas / (Gamma - 1.0f);

        public Crankcase(Crankshaft crankshaft,
                         float crankRadius, float conrodLength, float bore,
                         float clearanceVolume,
                         float initialPressure, float initialTemperature)
        {
            _crankshaft = crankshaft;
            _crankRadius = crankRadius;
            _conrodLength = conrodLength;
            _pistonArea = MathF.PI * 0.25f * bore * bore;
            _clearanceVolume = clearanceVolume;
            _obliquity = crankRadius / conrodLength;

            Pressure = initialPressure;
            Temperature = initialTemperature;
            float rho = initialPressure / (Rgas * initialTemperature);
            _mass = rho * clearanceVolume;
            _internalEnergy = _mass * Cv * initialTemperature;
            _airFraction = 1.0f;
            Volume = clearanceVolume;
            _prevVolume = Volume;

            IntakePort = new Port { Owner = this };
            TransferPort = new Port { Owner = this };
            _ports = new List<Port> { IntakePort, TransferPort };
        }

        public void PreStep(float dt)
        {
            // Save previous volume before updating
            _prevVolume = Volume;

            float theta = _crankshaft.CrankAngleRad % (2f * MathF.PI);
            float cosTh = MathF.Cos(theta);
            float sinTh = MathF.Sin(theta);
            float term = MathF.Sqrt(1f - _obliquity * _obliquity * sinTh * sinTh);

            // FIX: correct piston displacement: downstroke reduces crankcase volume
            float x = _crankRadius * (1f - cosTh) + _conrodLength * (1f - term);
            // Maximum volume at TDC (x = 0), minimum at BDC (x = stroke)
            float maxVolume = _clearanceVolume + _pistonArea * 2f * _crankRadius; // stroke = 2 * crankRadius
            Volume = maxVolume - _pistonArea * x;

            // Update thermodynamic state using the new volume (before mass transfer)
            if (_mass > 1e-12f && Volume > 1e-12f)
            {
                Temperature = _internalEnergy / (_mass * Cv);
                Pressure = _mass * Rgas * Temperature / Volume;
            }
        }

        public void UpdateState(float dt)
        {
            // ---- Mass and energy transport (identical to original) ----
            float mdotIn = IntakePort.MassFlowRate;
            float mdotOut = TransferPort.MassFlowRate;

            float dm = (mdotIn - mdotOut) * dt;
            float dE = (mdotIn * IntakePort.SpecificEnthalpy
                      - mdotOut * (Cv * Temperature + Pressure / MathF.Max(Density, 1e-12f))) * dt;
            float dY = (mdotIn * IntakePort.AirFraction - mdotOut * _airFraction) * dt;

            _mass += dm;
            _internalEnergy += dE;
            if (_mass > 1e-12f)
                _airFraction = Math.Clamp((_airFraction * (_mass - dm) + dY) / _mass, 0f, 1f);

            // ---- FIX: add mechanical work done BY the gas ON the piston ----
            // During a step, volume changed from _prevVolume to current Volume.
            // Work done BY gas = P * dV (if dV > 0, gas expands and does work, losing energy)
            float dV = Volume - _prevVolume;
            // Use average pressure during the step (approximate with current pressure)
            _internalEnergy -= Pressure * dV;   // removes energy when volume increases

            // Safety floors
            if (_mass < 1e-9f)
            {
                _mass = 1e-9f;
                _internalEnergy = _mass * Cv * 300f;
                _airFraction = 1f;
            }
            if (_internalEnergy < 0f)
                _internalEnergy = _mass * Cv * 300f;

            // Final state update
            if (_mass > 1e-12f && Volume > 1e-12f)
            {
                Temperature = _internalEnergy / (_mass * Cv);
                Pressure = _mass * Rgas * Temperature / Volume;
            }
            else
            {
                Temperature = 300f;
                Pressure = 101325f;
            }

            // Safety limits (unchanged, but now rarely triggered)
            const float safetyPressure = 1.0f;
            if (Pressure < safetyPressure && _mass > 1e-12f && Volume > 1e-12f)
            {
                Temperature = safetyPressure * Volume / (_mass * Rgas);
                _internalEnergy = _mass * Cv * Temperature;
                Pressure = safetyPressure;
            }

            const float maxPressure = 5e5f;
            if (Pressure > maxPressure && _mass > 1e-12f && Volume > 1e-12f)
            {
                float targetMass = maxPressure * Volume / (Rgas * Temperature);
                if (_mass > targetMass)
                {
                    _mass = targetMass;
                    _internalEnergy = _mass * Cv * Temperature;
                    Pressure = maxPressure;
                }
            }
        }
    }
}