using System; using System.Collections.Generic; using FluidSim.Interfaces; namespace FluidSim.Components { public class Cylinder : IComponent { public Port IntakePort { get; } public Port ExhaustPort { get; } public Crankshaft Crankshaft { get; } private readonly Port[] _ports; IReadOnlyList IComponent.Ports => _ports; public float Bore { get; } public float Stroke { get; } public float ConRodLength { get; } public float CompressionRatio { get; } public float IVO, IVC, EVO, EVC; // degrees public float IntakeValveDiameter = 0.03f; public float ExhaustValveDiameter = 0.028f; public float IntakeValveLift = 0.005f; public float ExhaustValveLift = 0.005f; public float IntakeValveMaxArea => MathF.PI * IntakeValveDiameter * IntakeValveLift; public float ExhaustValveMaxArea => MathF.PI * ExhaustValveDiameter * ExhaustValveLift; public float SparkAdvance = 20f; public float WiebeA = 5f, WiebeM = 2f, WiebeDuration = 60f, WiebeStart = 5f; public float StoichiometricAFR = 14.7f; public float FuelLowerHeatingValue = 44e6f; public float EnergyVariationFraction = 0.05f; public float MisfireProbability = 0.0f; public float CylinderWallArea = 0.02f; public float HeatTransferCoefficient = 100f; public float AmbientTemperature = 300f; public float PhaseOffset; // rad public float Volume => cylinderVolume; public float Pressure => (Gamma - 1f) * cylinderEnergy / MathF.Max(cylinderVolume, 1e-12f); public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f); public float Density => Mass / MathF.Max(cylinderVolume, 1e-12f); public float Mass => _airMass + _exhaustMass; public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f); public float PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume; private float cylinderVolume, cylinderEnergy; private float _airMass, _exhaustMass; private float trappedAirMass, fuelMass, burnFraction; private bool combustionActive, fuelInjected; private float _energyFactor = 1f; private readonly Random _random = new Random(); private const float Gamma = 1.4f; private const float GasConstant = 287f; private const float MaxPressurePa = 200e5f; private const float MaxTemperatureK = 3500f; public Cylinder(float bore, float stroke, float conRodLength, float compressionRatio, float ivo, float ivc, float evo, float evc, Crankshaft crankshaft) { Bore = bore; Stroke = stroke; ConRodLength = conRodLength; CompressionRatio = compressionRatio; IVO = ivo; IVC = ivc; EVO = evo; EVC = evc; Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft)); cylinderVolume = clearanceVolume; float initRho = 1.225f; _airMass = initRho * clearanceVolume; _exhaustMass = 0f; cylinderEnergy = 101325f * clearanceVolume / (Gamma - 1f); IntakePort = new Port { Owner = this }; ExhaustPort = new Port { Owner = this }; _ports = new[] { IntakePort, ExhaustPort }; } private float SweptVolume => MathF.PI * 0.25f * Bore * Bore * Stroke; private float clearanceVolume => SweptVolume / (CompressionRatio - 1f); private float CrankRadius => Stroke * 0.5f; private float Obliquity => CrankRadius / ConRodLength; private float CrankDeg => ((Crankshaft.CrankAngle + PhaseOffset) % (4f * MathF.PI)) * 180f / MathF.PI % 720f; public float ComputeVolume(float thetaRad) { float r = CrankRadius, l = ConRodLength; float cosTh = MathF.Cos(thetaRad), sinTh = MathF.Sin(thetaRad); float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh); float x = r * (1f - cosTh) + l * (1f - term); float area = MathF.PI * 0.25f * Bore * Bore; return clearanceVolume + area * x; } private float ValveLift(float thetaDeg, float opens, float closes, float peakLift) { float deg = thetaDeg % 720f; if (deg < 0f) deg += 720f; float duration; float effectiveOpen = opens; float effectiveClose = closes; if (closes < opens) { effectiveClose += 720f; } duration = effectiveClose - effectiveOpen; if (duration <= 0f) return 0f; float mapped = deg; if (mapped < opens) mapped += 720f; if (mapped < opens || mapped > effectiveClose) return 0f; float rampDur = duration * 0.25f; float holdDur = duration - 2f * rampDur; if (mapped >= opens && mapped < opens + rampDur) { float t = (mapped - opens) / rampDur; return peakLift * t * t * (3f - 2f * t); } else if (mapped >= opens + rampDur && mapped < opens + rampDur + holdDur) { return peakLift; } else if (mapped >= opens + rampDur + holdDur && mapped <= effectiveClose) { float t = (mapped - (opens + rampDur + holdDur)) / rampDur; return peakLift * (1f - t) * (1f - t) * (1f + 2f * t); } return 0f; } public float IntakeValveArea => MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift); public float ExhaustValveArea => MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift); private float Wiebe(float angleSinceSpark) { if (angleSinceSpark < WiebeStart) return 0f; float phi = (angleSinceSpark - WiebeStart) / WiebeDuration; if (phi <= 0f) return 0f; return 1f - MathF.Exp(-WiebeA * MathF.Pow(phi, WiebeM + 1f)); } public void PreStep(float dt) { // Speed‑dependent spark advance (simple linear) float rpm = Crankshaft.AngularVelocity * 60f / (2f * MathF.PI); SparkAdvance = Math.Clamp(10f + rpm * 0.002f, 5f, 40f); // 10° at idle, ~30° at 10k rpm float prevVolume = cylinderVolume; float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset; cylinderVolume = ComputeVolume(crankAngleRad); float dV = cylinderVolume - prevVolume; float pRel = Pressure - 101325f; float sinTh = MathF.Sin(crankAngleRad), cosTh = MathF.Cos(crankAngleRad); float term = MathF.Sqrt(1f - Obliquity * Obliquity * sinTh * sinTh); float dxdtheta = CrankRadius * sinTh * (1f + Obliquity * cosTh / term); float pistonArea = MathF.PI * 0.25f * Bore * Bore; Crankshaft.AddTorque(pRel * pistonArea * dxdtheta); cylinderEnergy -= Pressure * dV; float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % 720f; float currDeg = crankAngleRad * 180f / MathF.PI % 720f; // Intake closing – triggers fuel injection if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC) { trappedAirMass = _airMass; fuelMass = trappedAirMass / StoichiometricAFR; fuelInjected = true; } // Spark float sparkAngle = 0f - SparkAdvance; if (sparkAngle < 0f) sparkAngle += 720f; bool crossedSpark = (prevDeg < sparkAngle && currDeg >= sparkAngle) || (prevDeg > sparkAngle + 360f && currDeg < sparkAngle); if (crossedSpark && !combustionActive && fuelInjected) { if (_random.NextDouble() < MisfireProbability) { combustionActive = false; } else { combustionActive = true; burnFraction = 0f; float range = EnergyVariationFraction; _energyFactor = 1f + range * (2f * (float)_random.NextDouble() - 1f); } } // Combustion if (combustionActive) { float angleSinceSpark = currDeg - sparkAngle; if (angleSinceSpark < 0f) angleSinceSpark += 720f; float newFraction = Wiebe(angleSinceSpark); if (newFraction >= 1f || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance)) { newFraction = 1f; combustionActive = false; float totalMass = _airMass + _exhaustMass; _airMass = 0f; _exhaustMass = totalMass; } fuelInjected = false; float dFraction = newFraction - burnFraction; if (dFraction > 0f) { float dQ = fuelMass * FuelLowerHeatingValue * _energyFactor * dFraction; cylinderEnergy += dQ; _exhaustMass += fuelMass * dFraction; burnFraction = newFraction; } } // Heat loss float dQ_loss = HeatTransferCoefficient * CylinderWallArea * (Temperature - AmbientTemperature) * dt; cylinderEnergy -= dQ_loss; // Update port states float p = Pressure, rho = Density, T = Temperature; float h = Gamma / (Gamma - 1f) * p / MathF.Max(rho, 1e-12f); float af = AirFraction; IntakePort.Pressure = p; IntakePort.Density = rho; IntakePort.Temperature = T; IntakePort.SpecificEnthalpy = h; IntakePort.AirFraction = af; ExhaustPort.Pressure = p; ExhaustPort.Density = rho; ExhaustPort.Temperature = T; ExhaustPort.SpecificEnthalpy = h; ExhaustPort.AirFraction = af; } public void UpdateState(float dt) { float dmAir = 0f, dmExhaust = 0f, dE = 0f; foreach (var port in _ports) { float mdot = port.MassFlowRate; float af = mdot >= 0f ? port.AirFraction : AirFraction; dmAir += mdot * af * dt; dmExhaust += mdot * (1f - af) * dt; dE += mdot * port.SpecificEnthalpy * dt; } _airMass += dmAir; _exhaustMass += dmExhaust; cylinderEnergy += dE; float V = MathF.Max(cylinderVolume, 1e-12f); float currentP = (Gamma - 1f) * cylinderEnergy / V; if (currentP > MaxPressurePa) cylinderEnergy = MaxPressurePa * V / (Gamma - 1f); float currentRho = (_airMass + _exhaustMass) / V; float currentT = currentP / MathF.Max(currentRho * GasConstant, 1e-12f); if (currentT > MaxTemperatureK) { float pAtTlimit = currentRho * GasConstant * MaxTemperatureK; cylinderEnergy = pAtTlimit * V / (Gamma - 1f); } float totalMass = _airMass + _exhaustMass; if (totalMass < 1e-9f) { _airMass = 1e-9f; _exhaustMass = 0f; cylinderEnergy = 101325f * V / (Gamma - 1f); } else if (cylinderEnergy < 0f) { cylinderEnergy = 101325f * V / (Gamma - 1f); } if (_airMass < 0f) _airMass = 0f; if (_exhaustMass < 0f) _exhaustMass = 0f; } } }