diff --git a/Components/Atmosphere.cs b/Components/Atmosphere.cs index 8f07602..81b3c90 100644 --- a/Components/Atmosphere.cs +++ b/Components/Atmosphere.cs @@ -2,18 +2,15 @@ using FluidSim.Interfaces; namespace FluidSim.Components { - /// - /// Represents the ambient atmosphere – constant pressure/temperature reservoir. - /// public class Atmosphere : IComponent { - public double Pressure { get; set; } = 101325.0; - public double Temperature { get; set; } = 300.0; - public double GasConstant { get; set; } = 287.0; - public double Gamma => 1.4; + public float Pressure { get; set; } = 101325f; + public float Temperature { get; set; } = 300f; + public float GasConstant { get; set; } = 287f; + public float Gamma => 1.4f; - public double Density => Pressure / (GasConstant * Temperature); - public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density; + public float Density => Pressure / (GasConstant * Temperature); + public float SpecificEnthalpy => Gamma / (Gamma - 1f) * Pressure / Density; public Port Port { get; } @@ -25,9 +22,8 @@ namespace FluidSim.Components public IReadOnlyList Ports => new[] { Port }; - public void UpdateState(double dt) + public void UpdateState(float dt) { - // Atmosphere is static – just ensure the port reflects current values UpdatePort(); } @@ -37,7 +33,7 @@ namespace FluidSim.Components Port.Density = Density; Port.Temperature = Temperature; Port.SpecificEnthalpy = SpecificEnthalpy; - Port.AirFraction = 1.0; + Port.AirFraction = 1f; } } } \ No newline at end of file diff --git a/Components/Crankshaft.cs b/Components/Crankshaft.cs index 94fff97..cd4090e 100644 --- a/Components/Crankshaft.cs +++ b/Components/Crankshaft.cs @@ -1,54 +1,52 @@ -// Components/Crankshaft.cs using System; namespace FluidSim.Components { public class Crankshaft { - public double AngularVelocity { get; set; } // rad/s - public double CrankAngle { get; set; } // rad, 0 … 4π (four‑stroke cycle) - public double PreviousAngle { get; set; } // ← now has public setter + public float AngularVelocity; // rad/s + public float CrankAngle; // rad, 0 … 4π + public float PreviousAngle; - public double Inertia { get; set; } = 0.2; - public double FrictionConstant { get; set; } = 0.0; // N·m - public double FrictionViscous { get; set; } = 0.000; // N·m per rad/s + public float Inertia = 0.2f; + public float FrictionConstant; // N·m + public float FrictionViscous; // N·m per rad/s - private double externalTorque; + private float externalTorque; - public Crankshaft(double initialRPM = 400.0) + public Crankshaft(float initialRPM = 400f) { - AngularVelocity = initialRPM * 2.0 * Math.PI / 60.0; - CrankAngle = 0.0; - PreviousAngle = 0.0; + AngularVelocity = initialRPM * 2f * MathF.PI / 60f; + CrankAngle = 0f; + PreviousAngle = 0f; } - public void AddTorque(double torque) => externalTorque += torque; + public void AddTorque(float torque) => externalTorque += torque; - public void Step(double dt) + public void Step(float dt) { - // Catch NaN before it propagates - if (double.IsNaN(AngularVelocity) || double.IsInfinity(AngularVelocity)) - AngularVelocity = 0.0; - if (double.IsNaN(externalTorque) || double.IsInfinity(externalTorque)) - externalTorque = 0.0; + if (float.IsNaN(AngularVelocity) || float.IsInfinity(AngularVelocity)) + AngularVelocity = 0f; + if (float.IsNaN(externalTorque) || float.IsInfinity(externalTorque)) + externalTorque = 0f; PreviousAngle = CrankAngle; - double friction = FrictionConstant * Math.Sign(AngularVelocity) + FrictionViscous * AngularVelocity; - double netTorque = externalTorque - friction; - double alpha = netTorque / Inertia; + float friction = FrictionConstant * MathF.Sign(AngularVelocity) + + FrictionViscous * AngularVelocity; + float netTorque = externalTorque - friction; + float alpha = netTorque / Inertia; AngularVelocity += alpha * dt; - if (AngularVelocity < 0) AngularVelocity = 0; + if (AngularVelocity < 0f) AngularVelocity = 0f; CrankAngle += AngularVelocity * dt; + if (CrankAngle >= 4f * MathF.PI) + CrankAngle -= 4f * MathF.PI; + else if (CrankAngle < 0f) + CrankAngle += 4f * MathF.PI; - if (CrankAngle >= 4.0 * Math.PI) - CrankAngle -= 4.0 * Math.PI; - else if (CrankAngle < 0) - CrankAngle += 4.0 * Math.PI; - - externalTorque = 0.0; + externalTorque = 0f; } } } \ No newline at end of file diff --git a/Components/Cylinder.cs b/Components/Cylinder.cs index 4dabab5..3f14d2f 100644 --- a/Components/Cylinder.cs +++ b/Components/Cylinder.cs @@ -13,144 +13,103 @@ namespace FluidSim.Components private readonly Port[] _ports; IReadOnlyList IComponent.Ports => _ports; - // Geometry - public double Bore { get; } - public double Stroke { get; } - public double ConRodLength { get; } - public double CompressionRatio { get; } + public float Bore { get; } + public float Stroke { get; } + public float ConRodLength { get; } + public float CompressionRatio { get; } - // Valve timings (degrees, 0 = TDC compression, 720° full cycle) - public double IVO { get; } - public double IVC { get; } - public double EVO { get; } - public double EVC { 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; - // Valve geometry - public double IntakeValveDiameter { get; set; } = 0.030; - public double ExhaustValveDiameter { get; set; } = 0.028; - public double IntakeValveLift { get; set; } = 0.005; - public double ExhaustValveLift { get; set; } = 0.005; + public float IntakeValveMaxArea => MathF.PI * IntakeValveDiameter * IntakeValveLift; + public float ExhaustValveMaxArea => MathF.PI * ExhaustValveDiameter * ExhaustValveLift; - public double IntakeValveMaxArea => Math.PI * IntakeValveDiameter * IntakeValveLift; - public double ExhaustValveMaxArea => Math.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.01f; + public float CylinderWallArea = 0.02f; + public float HeatTransferCoefficient = 100f; + public float AmbientTemperature = 300f; - // Ignition and combustion - public double SparkAdvance { get; set; } = 20.0; - public double WiebeA { get; set; } = 5.0; - public double WiebeM { get; set; } = 2.0; - public double WiebeDuration { get; set; } = 60.0; - public double WiebeStart { get; set; } = 5.0; + public float PhaseOffset; // rad - // Fuel - public double StoichiometricAFR { get; set; } = 14.7; - public double FuelLowerHeatingValue { get; set; } = 44e6; + 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; - // Cycle‑to‑cycle randomness - public double EnergyVariationFraction { get; set; } = 0.05; - public double MisfireProbability { get; set; } = 0.01; - - // Heat loss - public double CylinderWallArea { get; set; } = 0.02; - public double HeatTransferCoefficient { get; set; } = 100.0; - public double AmbientTemperature { get; set; } = 300.0; - - // ---- Multi‑cylinder support ---- - /// - /// Phase offset (radians) added to the crankshaft angle for this cylinder. - /// Used for multi‑cylinder engines; set to 0 for single‑cylinder. - /// - public double PhaseOffset { get; set; } = 0.0; - - // State (public for drawing) - public double Volume => cylinderVolume; - public double Pressure => (Gamma - 1.0) * cylinderEnergy / Math.Max(cylinderVolume, 1e-12); - public double Temperature => Pressure / Math.Max(Density * GasConstant, 1e-12); - public double Density => Mass / Math.Max(cylinderVolume, 1e-12); - public double Mass => _airMass + _exhaustMass; - public double AirFraction => _airMass / Math.Max(Mass, 1e-12); - public double PistonFraction => (cylinderVolume - clearanceVolume) / SweptVolume; - - private double cylinderVolume; - private double cylinderEnergy; - private double _airMass; - private double _exhaustMass; - private double trappedAirMass; - private double fuelMass; - private double burnFraction; - private bool combustionActive; - private bool fuelInjected; - - private double _energyFactor = 1.0; + 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 double Gamma = 1.4; - private const double GasConstant = 287.0; + private const float Gamma = 1.4f; + private const float GasConstant = 287f; + private const float MaxPressurePa = 200e5f; + private const float MaxTemperatureK = 3500f; - private const double MaxPressurePa = 200e5; - private const double MaxTemperatureK = 3500.0; - - public Cylinder(double bore, double stroke, double conRodLength, double compressionRatio, - double ivo, double ivc, double evo, double evc, Crankshaft crankshaft) + 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; + Bore = bore; Stroke = stroke; ConRodLength = conRodLength; CompressionRatio = compressionRatio; - IVO = ivo; - IVC = ivc; - EVO = evo; - EVC = evc; - + IVO = ivo; IVC = ivc; EVO = evo; EVC = evc; Crankshaft = crankshaft ?? throw new ArgumentNullException(nameof(crankshaft)); cylinderVolume = clearanceVolume; - double initRho = 1.225; + float initRho = 1.225f; _airMass = initRho * clearanceVolume; - _exhaustMass = 0.0; - cylinderEnergy = 101325.0 * clearanceVolume / (Gamma - 1.0); + _exhaustMass = 0f; + cylinderEnergy = 101325f * clearanceVolume / (Gamma - 1f); IntakePort = new Port { Owner = this }; ExhaustPort = new Port { Owner = this }; _ports = new[] { IntakePort, ExhaustPort }; } - // Derived volumes - private double SweptVolume => Math.PI * 0.25 * Bore * Bore * Stroke; - private double clearanceVolume => SweptVolume / (CompressionRatio - 1.0); - private double CrankRadius => Stroke / 2.0; - private double Obliquity => CrankRadius / ConRodLength; + 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; - // Offset-aware crank angle in degrees - private double CrankDeg => - ((Crankshaft.CrankAngle + PhaseOffset) % (4.0 * Math.PI)) * 180.0 / Math.PI % 720.0; + private float CrankDeg => + ((Crankshaft.CrankAngle + PhaseOffset) % (4f * MathF.PI)) * 180f / MathF.PI % 720f; - public double ComputeVolume(double thetaRad) + public float ComputeVolume(float thetaRad) { - double r = CrankRadius; - double l = ConRodLength; - double cosTh = Math.Cos(thetaRad); - double sinTh = Math.Sin(thetaRad); - double term = Math.Sqrt(1.0 - Obliquity * Obliquity * sinTh * sinTh); - double x = r * (1.0 - cosTh) + l * (1.0 - term); - double area = Math.PI * 0.25 * Bore * Bore; + 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 double ValveLift(double thetaDeg, double opens, double closes, double peakLift) + private float ValveLift(float thetaDeg, float opens, float closes, float peakLift) { - double deg = thetaDeg % 720.0; - if (deg < 0) deg += 720.0; + float deg = thetaDeg % 720f; + if (deg < 0f) deg += 720f; + float duration = closes - opens; + if (duration <= 0f) return 0f; - double duration = closes - opens; - if (duration <= 0) return 0.0; - - double rampDur = duration * 0.25; - double holdDur = duration - 2.0 * rampDur; + float rampDur = duration * 0.25f; + float holdDur = duration - 2f * rampDur; if (deg >= opens && deg < opens + rampDur) { - double t = (deg - opens) / rampDur; - return peakLift * t * t * (3.0 - 2.0 * t); + float t = (deg - opens) / rampDur; + return peakLift * t * t * (3f - 2f * t); } else if (deg >= opens + rampDur && deg < opens + rampDur + holdDur) { @@ -158,54 +117,45 @@ namespace FluidSim.Components } else if (deg >= opens + rampDur + holdDur && deg <= closes) { - double t = (deg - (opens + rampDur + holdDur)) / rampDur; - return peakLift * (1.0 - t) * (1.0 - t) * (1.0 + 2.0 * t); + float t = (deg - (opens + rampDur + holdDur)) / rampDur; + return peakLift * (1f - t) * (1f - t) * (1f + 2f * t); } - return 0.0; + return 0f; } - public double IntakeValveArea => - Math.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift); + public float IntakeValveArea => + MathF.PI * IntakeValveDiameter * ValveLift(CrankDeg, IVO, IVC, IntakeValveLift); + public float ExhaustValveArea => + MathF.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift); - public double ExhaustValveArea => - Math.PI * ExhaustValveDiameter * ValveLift(CrankDeg, EVO, EVC, ExhaustValveLift); - - private double Wiebe(double angleSinceSpark) + private float Wiebe(float angleSinceSpark) { - if (angleSinceSpark < WiebeStart) return 0.0; - double phi = (angleSinceSpark - WiebeStart) / WiebeDuration; - if (phi <= 0) return 0.0; - return 1.0 - Math.Exp(-WiebeA * Math.Pow(phi, WiebeM + 1)); + 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(double dt) + public void PreStep(float dt) { - double prevVolume = cylinderVolume; - - // ----- Use phase‑offset crank angle for this cylinder ----- - double crankAngleRad = Crankshaft.CrankAngle + PhaseOffset; + float prevVolume = cylinderVolume; + float crankAngleRad = Crankshaft.CrankAngle + PhaseOffset; cylinderVolume = ComputeVolume(crankAngleRad); - double dV = cylinderVolume - prevVolume; - - // Piston torque - double pRel = Pressure - 101325.0; - double sinTh = Math.Sin(crankAngleRad); - double cosTh = Math.Cos(crankAngleRad); - double term = Math.Sqrt(1.0 - Obliquity * Obliquity * sinTh * sinTh); - double dxdtheta = CrankRadius * sinTh * (1.0 + Obliquity * cosTh / term); - double pistonArea = Math.PI * 0.25 * Bore * Bore; - double torque = pRel * pistonArea * dxdtheta; - Crankshaft.AddTorque(torque); + 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; - // Also use offset angle for event detection - double crankshaftPrevAngle = Crankshaft.PreviousAngle; - double prevDeg = (crankshaftPrevAngle + PhaseOffset) * 180.0 / Math.PI % 720.0; - double currDeg = crankAngleRad * 180.0 / Math.PI % 720.0; + float prevDeg = (Crankshaft.PreviousAngle + PhaseOffset) * 180f / MathF.PI % 720f; + float currDeg = crankAngleRad * 180f / MathF.PI % 720f; - // ----- Intake closing: capture trapped air mass and compute fuel ----- + // Intake closing if (prevDeg >= IVO && prevDeg < IVC && currDeg >= IVC) { trappedAirMass = _airMass; @@ -213,122 +163,103 @@ namespace FluidSim.Components fuelInjected = true; } - // ----- Spark ignition ----- - double sparkAngle = 0.0 - SparkAdvance; - if (sparkAngle < 0) sparkAngle += 720.0; - + // Spark + float sparkAngle = 0f - SparkAdvance; + if (sparkAngle < 0f) sparkAngle += 720f; bool crossedSpark = (prevDeg < sparkAngle && currDeg >= sparkAngle) || - (prevDeg > sparkAngle + 360.0 && currDeg < sparkAngle); + (prevDeg > sparkAngle + 360f && currDeg < sparkAngle); if (crossedSpark && !combustionActive && fuelInjected) { - bool misfire = _random.NextDouble() < MisfireProbability; - if (misfire) + if (_random.NextDouble() < MisfireProbability) { combustionActive = false; } else { - combustionActive = true; - burnFraction = 0.0; - double range = EnergyVariationFraction; - _energyFactor = 1.0 + range * (2.0 * _random.NextDouble() - 1.0); + combustionActive = true; burnFraction = 0f; + float range = EnergyVariationFraction; + _energyFactor = 1f + range * (2f * (float)_random.NextDouble() - 1f); } } - // ----- Combustion progress ----- + // Combustion if (combustionActive) { - double angleSinceSpark = currDeg - sparkAngle; - if (angleSinceSpark < 0) angleSinceSpark += 720.0; - double newFraction = Wiebe(angleSinceSpark); - - if (newFraction >= 1.0 || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance)) + float angleSinceSpark = currDeg - sparkAngle; + if (angleSinceSpark < 0f) angleSinceSpark += 720f; + float newFraction = Wiebe(angleSinceSpark); + if (newFraction >= 1f || angleSinceSpark > (WiebeDuration + WiebeStart + SparkAdvance)) { - newFraction = 1.0; - combustionActive = false; - double totalMass = _airMass + _exhaustMass; - _airMass = 0.0; - _exhaustMass = totalMass; + newFraction = 1f; combustionActive = false; + float totalMass = _airMass + _exhaustMass; + _airMass = 0f; _exhaustMass = totalMass; } - double dFraction = newFraction - burnFraction; - if (dFraction > 0) + float dFraction = newFraction - burnFraction; + if (dFraction > 0f) { - double dQ = fuelMass * FuelLowerHeatingValue * _energyFactor * dFraction; + float dQ = fuelMass * FuelLowerHeatingValue * _energyFactor * dFraction; cylinderEnergy += dQ; _exhaustMass += fuelMass * dFraction; burnFraction = newFraction; } } - // ----- Heat loss ----- - double dQ_loss = HeatTransferCoefficient * CylinderWallArea * - (Temperature - AmbientTemperature) * dt; + // Heat loss + float dQ_loss = HeatTransferCoefficient * CylinderWallArea * + (Temperature - AmbientTemperature) * dt; cylinderEnergy -= dQ_loss; // Update port states - double p = Pressure, rho = Density, T = Temperature; - double h = Gamma / (Gamma - 1.0) * p / Math.Max(rho, 1e-12); - double 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; + 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(double dt) + public void UpdateState(float dt) { - double dmAir = 0.0, dmExhaust = 0.0, dE = 0.0; - + float dmAir = 0f, dmExhaust = 0f, dE = 0f; foreach (var port in _ports) { - double mdot = port.MassFlowRate; - double af = mdot >= 0 ? port.AirFraction : AirFraction; + float mdot = port.MassFlowRate; + float af = mdot >= 0f ? port.AirFraction : AirFraction; dmAir += mdot * af * dt; - dmExhaust += mdot * (1.0 - af) * dt; + dmExhaust += mdot * (1f - af) * dt; dE += mdot * port.SpecificEnthalpy * dt; } - _airMass += dmAir; - _exhaustMass += dmExhaust; + _airMass += dmAir; _exhaustMass += dmExhaust; cylinderEnergy += dE; - double V = Math.Max(cylinderVolume, 1e-12); + float V = MathF.Max(cylinderVolume, 1e-12f); + float currentP = (Gamma - 1f) * cylinderEnergy / V; + if (currentP > MaxPressurePa) cylinderEnergy = MaxPressurePa * V / (Gamma - 1f); - double currentP = (Gamma - 1.0) * cylinderEnergy / V; - if (currentP > MaxPressurePa) - cylinderEnergy = MaxPressurePa * V / (Gamma - 1.0); - - double currentRho = (_airMass + _exhaustMass) / V; - double currentT = currentP / Math.Max(currentRho * GasConstant, 1e-12); + float currentRho = (_airMass + _exhaustMass) / V; + float currentT = currentP / MathF.Max(currentRho * GasConstant, 1e-12f); if (currentT > MaxTemperatureK) { - double pAtTlimit = currentRho * GasConstant * MaxTemperatureK; - cylinderEnergy = pAtTlimit * V / (Gamma - 1.0); + float pAtTlimit = currentRho * GasConstant * MaxTemperatureK; + cylinderEnergy = pAtTlimit * V / (Gamma - 1f); } - double totalMass = _airMass + _exhaustMass; - if (totalMass < 1e-9) + float totalMass = _airMass + _exhaustMass; + if (totalMass < 1e-9f) { - _airMass = 1e-9; - _exhaustMass = 0.0; - cylinderEnergy = 101325.0 * V / (Gamma - 1.0); + _airMass = 1e-9f; _exhaustMass = 0f; + cylinderEnergy = 101325f * V / (Gamma - 1f); } - else if (cylinderEnergy < 0.0) + else if (cylinderEnergy < 0f) { - cylinderEnergy = 101325.0 * V / (Gamma - 1.0); + cylinderEnergy = 101325f * V / (Gamma - 1f); } - if (_airMass < 0.0) _airMass = 0.0; - if (_exhaustMass < 0.0) _exhaustMass = 0.0; + if (_airMass < 0f) _airMass = 0f; + if (_exhaustMass < 0f) _exhaustMass = 0f; } } } \ No newline at end of file diff --git a/Components/Nozzleflow.cs b/Components/Nozzleflow.cs deleted file mode 100644 index a365f71..0000000 --- a/Components/Nozzleflow.cs +++ /dev/null @@ -1,46 +0,0 @@ -using System; - -namespace FluidSim.Components -{ - public static class NozzleFlow - { - /// - /// Computes the nozzle‑exit primitive state and mass flow rate from a - /// volume to a pipe, using isentropic relations. Follows ensim4's flow() logic. - /// - public static void Compute(double Pt_high, double Tt_high, - double P_low, double gamma, double R, double area, - out double rhoExit, out double uExit, - out double pExit, out double mdot) - { - double gm1 = gamma - 1.0; - double Pt_over_Ps = Pt_high / P_low; - - // Mach number (subsonic, clamped to 1) - double M = Math.Sqrt(Math.Max(0.0, - (2.0 / gm1) * (Math.Pow(Pt_over_Ps, gm1 / gamma) - 1.0))); - if (M > 1.0) M = 1.0; - - double T_star = Tt_high / (1.0 + 0.5 * gm1 * M * M); - double a_star = Math.Sqrt(gamma * R * T_star); - double u_star = M * a_star; - pExit = Pt_high * Math.Pow(1.0 + 0.5 * gm1 * M * M, -gamma / gm1); - rhoExit = pExit / (R * T_star); - uExit = u_star; // positive away from high‑pressure side - mdot = rhoExit * uExit * area; - } - - /// - /// Ambient cell for non‑reflecting open end (ensim4 calc_ambient_cell). - /// - public static void ComputeAmbientCell(double rhoInt, double uInt, double pInt, - double pAmbient, double gamma, - out double rhoAmb, out double uAmb, - out double pAmb) - { - pAmb = pAmbient; - uAmb = uInt; - rhoAmb = rhoInt * Math.Pow(pAmb / pInt, 1.0 / gamma); - } - } -} \ No newline at end of file diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs deleted file mode 100644 index a95301d..0000000 --- a/Components/Pipe1D.cs +++ /dev/null @@ -1,451 +0,0 @@ -using System; -using System.Diagnostics; -using FluidSim.Interfaces; - -namespace FluidSim.Components -{ - /// - /// 1‑D compressible Euler pipe with Lax‑Friedrichs finite‑volume scheme. - /// Ghost states are set externally via SetGhostLeft/Right; they are always required. - /// Now includes a passive scalar for air mass fraction. - /// - public class Pipe1D : IComponent - { - public const bool EnableDetailedProfiling = false; // set to false in release builds - - public Port PortA { get; } - public Port PortB { get; } - public double Area { get; } - public double DampingMultiplier { get; set; } = 10.0; - public double EnergyRelaxationRate { get; set; } = 5.0; // 1/s - public string Name = "Pipe"; - - private double _ambientPressure = 101325.0; - public double AmbientPressure - { - get => _ambientPressure; - set - { - _ambientPressure = value; - _ambientEnergyReference = value / (_gamma - 1.0); - } - } - - private readonly int _n; - private readonly double _dx; - private readonly double _diameter; - private readonly double _gamma = 1.4; - - private double[] _rho, _rhou, _E; - private double[] _Y; // air mass fraction - private double[] _fluxM, _fluxP, _fluxE; - - private double _rhoGhostL, _uGhostL, _pGhostL, _YGhostL; - private double _rhoGhostR, _uGhostR, _pGhostR, _YGhostR; - private bool _ghostLValid, _ghostRValid; - - private double _laminarCoeff; - private double _ambientEnergyReference; - - // ---------- Profiling accumulators ---------- - private long _profPrecomputeTicks; - private long _profLeftFluxTicks; - private long _profInteriorLoopTicks; - private long _profRightFluxTicks; - private long _profPortUpdateTicks; - private long _profCallCount; - - public Pipe1D(double length, double area, int cellCount) - { - if (cellCount < 4) throw new ArgumentException("cellCount must be at least 4"); - - _n = cellCount; - _dx = length / _n; - Area = area; - _diameter = 2.0 * Math.Sqrt(area / Math.PI); - - _rho = new double[_n]; - _rhou = new double[_n]; - _E = new double[_n]; - _Y = new double[_n]; - _fluxM = new double[_n + 1]; - _fluxP = new double[_n + 1]; - _fluxE = new double[_n + 1]; - - double mu_air = 1.8e-5; - double radius = _diameter * 0.5; - _laminarCoeff = 8.0 * mu_air / (radius * radius); - - _ambientEnergyReference = 101325.0 / (_gamma - 1.0); - - PortA = new Port { Owner = this }; - PortB = new Port { Owner = this }; - - SetUniformState(1.225, 0.0, 101325.0); - } - - IReadOnlyList IComponent.Ports => new[] { PortA, PortB }; - - public void UpdateState(double dt) { } - - // ---------- Ghost interface ---------- - public void SetGhostLeft(double rho, double u, double p, double airFraction) - { - _rhoGhostL = rho; _uGhostL = u; _pGhostL = p; _YGhostL = airFraction; _ghostLValid = true; - } - public void SetGhostRight(double rho, double u, double p, double airFraction) - { - _rhoGhostR = rho; _uGhostR = u; _pGhostR = p; _YGhostR = airFraction; _ghostRValid = true; - } - public void ClearGhostFlags() { _ghostLValid = false; _ghostRValid = false; } - - public double GetInteriorAirFractionLeft() => _Y[0]; - public double GetInteriorAirFractionRight() => _Y[_n - 1]; - - public (double rho, double u, double p) GetInteriorStateLeft() - { - double rho = Math.Max(_rho[0], 1e-12); - double u = _rhou[0] / rho; - double p = PressureScalar(0); - return (rho, u, p); - } - public (double rho, double u, double p) GetInteriorStateRight() - { - double rho = Math.Max(_rho[_n - 1], 1e-12); - double u = _rhou[_n - 1] / rho; - double p = PressureScalar(_n - 1); - return (rho, u, p); - } - public int CellCount => _n; - public double GetCellDensity(int i) => _rho[i]; - public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12); - public double GetCellPressure(int i) => PressureScalar(i); - - public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8) - { - double maxW = 0.0; - for (int i = 0; i < _n; i++) - { - double rho = Math.Max(_rho[i], 1e-12); - double u = Math.Abs(_rhou[i] / rho); - double p = PressureScalar(i); - double c = Math.Sqrt(_gamma * p / rho); - double local = u + c; - if (local > maxW) maxW = local; - } - maxW = Math.Max(maxW, 1e-8); - return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx))); - } - - // ---------- Main step (per sub‑step) ---------- - public void SimulateSingleStep(double dtSub) - { - if (!_ghostLValid || !_ghostRValid) - throw new InvalidOperationException("Ghost cells not set before SimulateSingleStep."); - - double dt = dtSub; - int n = _n; - double dt_dx = dt / _dx; - double coeff = _laminarCoeff * DampingMultiplier; - double relaxRate = EnergyRelaxationRate; - double gamma = _gamma; - double gm1 = gamma - 1.0; - - // ---------- Profiling start ---------- - long t0 = 0, t1 = 0; - if (EnableDetailedProfiling) - { - t0 = Stopwatch.GetTimestamp(); - _profCallCount++; - } - - // ---------- Phase 1: Pre‑compute pressure and speed of sound ---------- - double[] p = new double[n]; - double[] c = new double[n]; - for (int i = 0; i < n; i++) - { - double rho = Math.Max(_rho[i], 1e-12); - double u = _rhou[i] / rho; - p[i] = gm1 * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho); - c[i] = Math.Sqrt(gamma * p[i] / rho); - } - - if (EnableDetailedProfiling) - { - t1 = Stopwatch.GetTimestamp(); - _profPrecomputeTicks += (t1 - t0); - t0 = t1; - } - - // ---------- Local flux functions ---------- - void LaxFlux(double rL, double uL, double pL, double cL, - double rR, double uR, double pR, double cR, - out double fm, out double fp, out double fe) - { - double EL = pL / (gm1 * rL) + 0.5 * uL * uL; - double ER = pR / (gm1 * rR) + 0.5 * uR * uR; - double Fm_L = rL * uL; - double Fp_L = rL * uL * uL + pL; - double Fe_L = (rL * EL + pL) * uL; - double Fm_R = rR * uR; - double Fp_R = rR * uR * uR + pR; - double Fe_R = (rR * ER + pR) * uR; - double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR); - fm = 0.5 * (Fm_L + Fm_R) - 0.5 * alpha * (rR - rL); - fp = 0.5 * (Fp_L + Fp_R) - 0.5 * alpha * (rR * uR - rL * uL); - fe = 0.5 * (Fe_L + Fe_R) - 0.5 * alpha * (rR * ER - rL * EL); - } - - void ScalarFlux(double rL, double uL, double cL, double YL, - double rR, double uR, double cR, double YR, - double alpha, out double fy) - { - double Fm_L = rL * uL; - double Fm_R = rR * uR; - fy = 0.5 * (Fm_L * YL + Fm_R * YR) - 0.5 * alpha * (rR * YR - rL * YL); - } - - // ---------- Phase 2: Left face flux (ghostL – cell 0) ---------- - double rL_ghost = Math.Max(_rhoGhostL, 1e-12); - double pL_ghost = _pGhostL; - double uL_ghost = _uGhostL; - double cL_ghost = Math.Sqrt(gamma * pL_ghost / rL_ghost); - - LaxFlux(rL_ghost, uL_ghost, pL_ghost, cL_ghost, - _rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), p[0], c[0], - out double fluxM_left, out double fluxP_left, out double fluxE_left); - - double alphaLeft = Math.Max(Math.Abs(uL_ghost) + cL_ghost, - Math.Abs(_rhou[0] / Math.Max(_rho[0], 1e-12)) + c[0]); - ScalarFlux(rL_ghost, uL_ghost, cL_ghost, _YGhostL, - _rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), c[0], _Y[0], - alphaLeft, out double fluxY_left); - - if (EnableDetailedProfiling) - { - t1 = Stopwatch.GetTimestamp(); - _profLeftFluxTicks += (t1 - t0); - t0 = t1; - } - - // ---------- Phase 3: Interior loop (fluxes + cell updates) ---------- - double fluxM_prev = fluxM_left; - double fluxP_prev = fluxP_left; - double fluxE_prev = fluxE_left; - double fluxY_prev = fluxY_left; - - for (int i = 0; i < n - 1; i++) - { - int iL = i; - int iR = i + 1; - - double rL = Math.Max(_rho[iL], 1e-12); - double uL = _rhou[iL] / rL; - double pL = p[iL]; - double cL = c[iL]; - double YL = _Y[iL]; - - double rR = Math.Max(_rho[iR], 1e-12); - double uR = _rhou[iR] / rR; - double pR = p[iR]; - double cR = c[iR]; - double YR = _Y[iR]; - - LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, - out double fluxM_right, out double fluxP_right, out double fluxE_right); - - double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR); - ScalarFlux(rL, uL, cL, YL, rR, uR, cR, YR, alpha, out double fluxY_right); - - // Update cell i - double r = _rho[i]; - double ru = _rhou[i]; - double E = _E[i]; - double Y = _Y[i]; - - double newR = r - dt_dx * (fluxM_right - fluxM_prev); - double newRu = ru - dt_dx * (fluxP_right - fluxP_prev); - double newE = E - dt_dx * (fluxE_right - fluxE_prev); - double oldRhoY = r * Y; - double newRhoY = oldRhoY - dt_dx * (fluxY_right - fluxY_prev); - - double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt); - newRu *= dampingFactor; - double relaxFactor = Math.Exp(-relaxRate * dt); - newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor; - - newR = Math.Max(newR, 1e-12); - double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12); - double eMin = 100.0 / gm1 + kin; - newE = Math.Max(newE, eMin); - - _rho[i] = newR; - _rhou[i] = newRu; - _E[i] = newE; - _Y[i] = Math.Clamp(newRhoY / newR, 0.0, 1.0); - - fluxM_prev = fluxM_right; - fluxP_prev = fluxP_right; - fluxE_prev = fluxE_right; - fluxY_prev = fluxY_right; - } - - if (EnableDetailedProfiling) - { - t1 = Stopwatch.GetTimestamp(); - _profInteriorLoopTicks += (t1 - t0); - t0 = t1; - } - - // ---------- Phase 4: Right face flux (cell n‑1 – ghostR) ---------- - double rR_ghost = Math.Max(_rhoGhostR, 1e-12); - double pR_ghost = _pGhostR; - double uR_ghost = _uGhostR; - double cR_ghost = Math.Sqrt(gamma * pR_ghost / rR_ghost); - - double rInt = _rho[n - 1]; - double uInt = _rhou[n - 1] / Math.Max(rInt, 1e-12); - - LaxFlux(rInt, uInt, p[n - 1], c[n - 1], - rR_ghost, uR_ghost, pR_ghost, cR_ghost, - out double fluxM_right_final, out double fluxP_right_final, out double fluxE_right_final); - - double alphaRight = Math.Max(Math.Abs(uInt) + c[n - 1], Math.Abs(uR_ghost) + cR_ghost); - ScalarFlux(rInt, uInt, c[n - 1], _Y[n - 1], - rR_ghost, uR_ghost, cR_ghost, _YGhostR, - alphaRight, out double fluxY_right_final); - - // Update last cell - { - int i = n - 1; - double r = _rho[i]; - double ru = _rhou[i]; - double E = _E[i]; - double Y = _Y[i]; - - double newR = r - dt_dx * (fluxM_right_final - fluxM_prev); - double newRu = ru - dt_dx * (fluxP_right_final - fluxP_prev); - double newE = E - dt_dx * (fluxE_right_final - fluxE_prev); - double oldRhoY = r * Y; - double newRhoY = oldRhoY - dt_dx * (fluxY_right_final - fluxY_prev); - - double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt); - newRu *= dampingFactor; - double relaxFactor = Math.Exp(-relaxRate * dt); - newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor; - - newR = Math.Max(newR, 1e-12); - double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12); - double eMin = 100.0 / gm1 + kin; - newE = Math.Max(newE, eMin); - - _rho[i] = newR; - _rhou[i] = newRu; - _E[i] = newE; - _Y[i] = Math.Clamp(newRhoY / newR, 0.0, 1.0); - } - - if (EnableDetailedProfiling) - { - t1 = Stopwatch.GetTimestamp(); - _profRightFluxTicks += (t1 - t0); - t0 = t1; - } - - // ---------- Phase 5: Update port states ---------- - (double rhoA, double uA, double pA) = GetInteriorStateLeft(); - PortA.Pressure = pA; PortA.Density = rhoA; - PortA.Temperature = pA / (rhoA * 287.0); - PortA.SpecificEnthalpy = gm1 / (gamma - 1.0) * pA / rhoA; - PortA.AirFraction = _Y[0]; - - (double rhoB, double uB, double pB) = GetInteriorStateRight(); - PortB.Pressure = pB; PortB.Density = rhoB; - PortB.Temperature = pB / (rhoB * 287.0); - PortB.SpecificEnthalpy = gm1 / (gamma - 1.0) * pB / rhoB; - PortB.AirFraction = _Y[_n - 1]; - - if (EnableDetailedProfiling) - { - t1 = Stopwatch.GetTimestamp(); - _profPortUpdateTicks += (t1 - t0); - } - } - - private double PressureScalar(int i) - { - double rho = Math.Max(_rho[i], 1e-12); - return (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho); - } - - public void SetUniformState(double rho, double u, double p) - { - double e = p / ((_gamma - 1.0) * rho); - double E = rho * e + 0.5 * rho * u * u; - for (int i = 0; i < _n; i++) - { - _rho[i] = rho; - _rhou[i] = rho * u; - _E[i] = E; - _Y[i] = 1.0; // initially pure air - } - } - - public void SetCellState(int i, double rho, double u, double p) - { - if (i < 0 || i >= _n) return; - double e = p / ((_gamma - 1.0) * rho); - double E = rho * e + 0.5 * rho * u * u; - _rho[i] = rho; - _rhou[i] = rho * u; - _E[i] = E; - _Y[i] = 1.0; - } - - public void SetCellPressure(int i, double p) - { - if (i < 0 || i >= _n) return; - double rho = _rho[i]; - double u = _rhou[i] / rho; - double e = p / ((_gamma - 1.0) * rho); - _E[i] = rho * e + 0.5 * rho * u * u; - } - - // ---------- Public profiling interface ---------- - public void ResetDetailCounters() - { - _profPrecomputeTicks = 0; - _profLeftFluxTicks = 0; - _profInteriorLoopTicks = 0; - _profRightFluxTicks = 0; - _profPortUpdateTicks = 0; - _profCallCount = 0; - } - - public string GetDetailProfileReport() - { - if (!EnableDetailedProfiling) - return "Detailed profiling disabled."; - - double freq = Stopwatch.Frequency; - long totalTicks = _profPrecomputeTicks + _profLeftFluxTicks + - _profInteriorLoopTicks + _profRightFluxTicks + - _profPortUpdateTicks; - - if (totalTicks == 0) return "No profiling data."; - - double totalSec = totalTicks / freq; - double avgCallSec = totalSec / _profCallCount; - double avgCallUs = avgCallSec * 1e6; - - string report = $" Pipe detailed (over {_profCallCount} calls, total {totalSec * 1000:F2} ms):\n"; - report += $" Avg per call: {avgCallUs:F2} µs\n"; - report += $" Precompute p,c: {_profPrecomputeTicks * 100.0 / totalTicks:F1} % ({_profPrecomputeTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n"; - report += $" Left face flux: {_profLeftFluxTicks * 100.0 / totalTicks:F1} % ({_profLeftFluxTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n"; - report += $" Interior loop: {_profInteriorLoopTicks * 100.0 / totalTicks:F1} % ({_profInteriorLoopTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n"; - report += $" Right face flux: {_profRightFluxTicks * 100.0 / totalTicks:F1} % ({_profRightFluxTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n"; - report += $" Port update: {_profPortUpdateTicks * 100.0 / totalTicks:F1} % ({_profPortUpdateTicks / freq * 1e6 / _profCallCount:F2} µs/call)\n"; - return report; - } - } -} \ No newline at end of file diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs index 35b517f..1fc1087 100644 --- a/Components/Volume0D.cs +++ b/Components/Volume0D.cs @@ -8,36 +8,40 @@ namespace FluidSim.Components { public List Ports { get; } = new List(); - private double _airMass; - private double _exhaustMass; - public double InternalEnergy { get; set; } - public double Volume { get; set; } - public double Dvdt { get; set; } - public double Gamma { get; set; } = 1.4; - public double GasConstant { get; set; } = 287.0; + private float _airMass; + private float _exhaustMass; + public float InternalEnergy; + public float Volume; + public float Dvdt; + public float Gamma { get; set; } = 1.4f; + public float GasConstant { get; set; } = 287f; + public float AmbientPressure { get; set; } = 101325f; - public double AmbientPressure { get; set; } = 101325.0; + // ---------- Thermal relaxation to environment ---------- + /// Rate of heat transfer to the surroundings (1/s). 0 = adiabatic. + public float EnergyRelaxationRate { get; set; } = 10f; + /// Temperature to relax toward (K). Default is room temperature. + public float AmbientTemperature { get; set; } = 300f; - // Derived quantities - public double Mass => _airMass + _exhaustMass; - public double AirFraction => _airMass / Math.Max(Mass, 1e-12); - public double Density => Mass / Math.Max(Volume, 1e-12); - public double Pressure => (Gamma - 1.0) * InternalEnergy / Math.Max(Volume, 1e-12); - public double Temperature => Pressure / Math.Max(Density * GasConstant, 1e-12); - public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Math.Max(Density, 1e-12); + public float Mass => _airMass + _exhaustMass; + public float AirFraction => _airMass / MathF.Max(Mass, 1e-12f); + public float Density => Mass / MathF.Max(Volume, 1e-12f); + public float Pressure => (Gamma - 1f) * InternalEnergy / MathF.Max(Volume, 1e-12f); + public float Temperature => Pressure / MathF.Max(Density * GasConstant, 1e-12f); + public float SpecificEnthalpy => Gamma / (Gamma - 1f) * Pressure / MathF.Max(Density, 1e-12f); - public Volume0D(double initialVolume, double initialPressure, - double initialTemperature, double gasConstant = 287.0, double gamma = 1.4) + public Volume0D(float initialVolume, float initialPressure, + float initialTemperature, float gasConstant = 287f, float gamma = 1.4f) { GasConstant = gasConstant; Gamma = gamma; Volume = initialVolume; - Dvdt = 0.0; + Dvdt = 0f; - double rho0 = initialPressure / (GasConstant * initialTemperature); - _airMass = rho0 * Volume; // starts with all air - _exhaustMass = 0.0; - InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0); + float rho0 = initialPressure / (GasConstant * initialTemperature); + _airMass = rho0 * Volume; + _exhaustMass = 0f; + InternalEnergy = (initialPressure * Volume) / (Gamma - 1f); } public Port CreatePort() @@ -52,66 +56,75 @@ namespace FluidSim.Components return port; } - public void SetPressure(double pressure, double? temperature = null) + public void SetPressure(float pressure, float? temperature = null) { - double V = Math.Max(Volume, 1e-12); - double T = temperature ?? Temperature; - double rho = pressure / (GasConstant * T); - double totalMass = rho * V; - // Keep current air fraction when setting pressure? - double af = AirFraction; + float V = MathF.Max(Volume, 1e-12f); + float T = temperature ?? Temperature; + float rho = pressure / (GasConstant * T); + float totalMass = rho * V; + float af = AirFraction; _airMass = totalMass * af; - _exhaustMass = totalMass * (1.0 - af); - InternalEnergy = pressure * V / (Gamma - 1.0); + _exhaustMass = totalMass * (1f - af); + InternalEnergy = pressure * V / (Gamma - 1f); } - public void UpdateState(double dt) + public void UpdateState(float dt) { - double totalMdotAir = 0.0; - double totalMdotExhaust = 0.0; - double totalEdot = 0.0; - + float totalMdotAir = 0f, totalMdotExhaust = 0f, totalEdot = 0f; foreach (var port in Ports) { - double mdot = port.MassFlowRate; // positive INTO volume - double af = mdot >= 0 ? port.AirFraction : AirFraction; // inflow: use port's fraction; outflow: well-mixed + float mdot = port.MassFlowRate; + float af = mdot >= 0f ? port.AirFraction : AirFraction; totalMdotAir += mdot * af; - totalMdotExhaust += mdot * (1.0 - af); + totalMdotExhaust += mdot * (1f - af); totalEdot += mdot * port.SpecificEnthalpy; } - double dAir = totalMdotAir * dt; - double dExhaust = totalMdotExhaust * dt; - double dE = totalEdot * dt - Pressure * Dvdt * dt; + float dAir = totalMdotAir * dt; + float dExhaust = totalMdotExhaust * dt; + float dE = totalEdot * dt - Pressure * Dvdt * dt; _airMass += dAir; _exhaustMass += dExhaust; InternalEnergy += dE; - double V = Math.Max(Volume, 1e-12); - double totalMass = _airMass + _exhaustMass; - if (totalMass < 1e-9) + // ---- Thermal relaxation ---- + if (EnergyRelaxationRate > 0f) { - _airMass = 1e-9; - _exhaustMass = 0.0; - InternalEnergy = AmbientPressure * V / (Gamma - 1.0); - } - else if (InternalEnergy < 0.0) - { - InternalEnergy = AmbientPressure * V / (Gamma - 1.0); + float currentMass = Mass; + if (currentMass > 1e-12f) + { + // Target internal energy: current mass at ambient temperature + float targetE = currentMass * GasConstant * AmbientTemperature / (Gamma - 1f); + float relaxFactor = MathF.Exp(-EnergyRelaxationRate * dt); + InternalEnergy = targetE + (InternalEnergy - targetE) * relaxFactor; + } } - if (_airMass < 0.0) _airMass = 0.0; - if (_exhaustMass < 0.0) _exhaustMass = 0.0; + float V = MathF.Max(Volume, 1e-12f); + float totalMass = _airMass + _exhaustMass; + if (totalMass < 1e-9f) + { + _airMass = 1e-9f; + _exhaustMass = 0f; + InternalEnergy = AmbientPressure * V / (Gamma - 1f); + } + else if (InternalEnergy < 0f) + { + InternalEnergy = AmbientPressure * V / (Gamma - 1f); + } - double p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy, afrac = AirFraction; + if (_airMass < 0f) _airMass = 0f; + if (_exhaustMass < 0f) _exhaustMass = 0f; + + float p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy, afr = AirFraction; foreach (var port in Ports) { port.Pressure = p; port.Density = rho; port.Temperature = T; port.SpecificEnthalpy = h; - port.AirFraction = afrac; + port.AirFraction = afr; } } diff --git a/Core/BoundarySystem.cs b/Core/BoundarySystem.cs new file mode 100644 index 0000000..715d464 --- /dev/null +++ b/Core/BoundarySystem.cs @@ -0,0 +1,330 @@ +using FluidSim.Components; +using FluidSim.Interfaces; +using System; + +namespace FluidSim.Core +{ + public class BoundarySystem + { + public struct OrificeDesc + { + public Port VolumePort; + public int PipeIndex; + public bool IsLeftEnd; + public int AreaIndex; + public float DischargeCoeff; + + // --- Inertance support --- + public bool UseInertance; + public float EffectiveLength; + public float CurrentMdot; // kg/s, positive = volume → pipe + + // --- Dissipative loss --- + public float LossCoefficient; // K factor for pressure drop = K * 0.5*rho*u^2 + } + + public struct OpenEndDesc + { + public int PipeIndex; + public bool IsLeftEnd; + public float AmbientPressure; + public float Gamma; + public float PipeArea; + public float LastMassFlowRate; + public float LastFacePressure; + } + + private OrificeDesc[] _orifices; + private OpenEndDesc[] _openEnds; + private float[] _orificeAreas; + private PipeSystem _pipeSystem; + + public BoundarySystem(PipeSystem pipeSystem, int maxOrifices, int maxOpenEnds) + { + _pipeSystem = pipeSystem; + _orifices = new OrificeDesc[maxOrifices]; + _openEnds = new OpenEndDesc[maxOpenEnds]; + _orificeAreas = new float[maxOrifices]; + } + + public int OrificeCount { get; private set; } + public int OpenEndCount { get; private set; } + + public void AddOrifice(Port volumePort, int pipeIndex, bool isLeftEnd, + int areaIndex, float dischargeCoeff = 1f, + float lossCoefficient = 0f) + { + _orifices[OrificeCount] = new OrificeDesc + { + VolumePort = volumePort, + PipeIndex = pipeIndex, + IsLeftEnd = isLeftEnd, + AreaIndex = areaIndex, + DischargeCoeff = dischargeCoeff, + UseInertance = false, + EffectiveLength = 0f, + CurrentMdot = 0f, + LossCoefficient = lossCoefficient + }; + OrificeCount++; + } + + public void AddOrificeWithInertance(Port volumePort, int pipeIndex, bool isLeftEnd, + int areaIndex, float dischargeCoeff, + float effectiveLength, float lossCoefficient = 0f) + { + AddOrifice(volumePort, pipeIndex, isLeftEnd, areaIndex, dischargeCoeff, lossCoefficient); + ref var d = ref _orifices[OrificeCount - 1]; + d.UseInertance = true; + d.EffectiveLength = effectiveLength; + } + + public void AddOpenEnd(int pipeIndex, bool isLeftEnd, + float ambientPressure, float pipeArea, float gamma = 1.4f) + { + int idx = OpenEndCount; + _openEnds[idx] = new OpenEndDesc + { + PipeIndex = pipeIndex, + IsLeftEnd = isLeftEnd, + AmbientPressure = ambientPressure, + Gamma = gamma, + PipeArea = pipeArea + }; + OpenEndCount++; + } + + public void SetOrificeAreas(float[] areas) + { + for (int i = 0; i < OrificeCount; i++) + _orificeAreas[i] = areas[i]; + } + + public float GetOpenEndMassFlow(int openEndIndex) + { + if (openEndIndex < 0 || openEndIndex >= OpenEndCount) return 0f; + return _openEnds[openEndIndex].LastMassFlowRate; + } + + public float GetOpenEndPressure(int openEndIndex) + { + if (openEndIndex < 0 || openEndIndex >= OpenEndCount) return 101325f; + return _openEnds[openEndIndex].LastFacePressure; + } + + public void ResolveOrifices(float dt) + { + for (int i = 0; i < OrificeCount; i++) + { + ref var d = ref _orifices[i]; + float area = _orificeAreas[d.AreaIndex]; + if (area < 1e-12f || d.VolumePort == null) + { + // Closed wall – reflect interior state + var (rInt, uInt, pInt) = d.IsLeftEnd + ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex) + : _pipeSystem.GetInteriorStateRight(d.PipeIndex); + float afInt = d.IsLeftEnd + ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex) + : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex); + + if (d.IsLeftEnd) + _pipeSystem.SetGhostLeft(d.PipeIndex, rInt, -uInt, pInt, afInt); + else + _pipeSystem.SetGhostRight(d.PipeIndex, rInt, -uInt, pInt, afInt); + + if (d.VolumePort != null) d.VolumePort.MassFlowRate = 0f; + continue; + } + + // Gather states + float volP = d.VolumePort.Pressure; + float volRho = d.VolumePort.Density; + float volT = d.VolumePort.Temperature; + float volH = d.VolumePort.SpecificEnthalpy; + float volAF = d.VolumePort.AirFraction; + + var (pipeRho, pipeU, pipeP) = d.IsLeftEnd + ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex) + : _pipeSystem.GetInteriorStateRight(d.PipeIndex); + float pipeT = pipeP / MathF.Max(pipeRho * 287f, 1e-12f); + float pipeAF = d.IsLeftEnd + ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex) + : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex); + + float gamma = 1.4f, R = 287f, Cd = d.DischargeCoeff; + + // --- Preliminary nozzle solution (no loss) to estimate flow direction and velocity --- + float mdotEst, rhoFaceEst, uFaceEst, pFaceEst; + if (volP >= pipeP) + { + IsentropicOrifice.Compute(volP, volRho, volT, pipeP, gamma, R, area, Cd, + out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst); + } + else + { + IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, gamma, R, area, Cd, + out mdotEst, out rhoFaceEst, out uFaceEst, out pFaceEst); + mdotEst = -mdotEst; + } + + // --- Apply symmetric loss if LossCoefficient > 0 --- + float volP_eff = volP; + float pipeP_eff = pipeP; + if (d.LossCoefficient > 0f && MathF.Abs(mdotEst) > 1e-12f) + { + float rhoRef = mdotEst >= 0 ? rhoFaceEst : rhoFaceEst; // rhoFaceEst already reflects the correct side + float uRef = uFaceEst; + float dynP = 0.5f * rhoRef * uRef * uRef * d.LossCoefficient; + + // Clamp the loss to avoid overshoot (max 80% of pressure difference) + float dp = MathF.Abs(volP - pipeP); + dynP = MathF.Min(dynP, 0.8f * dp); + + // Apply symmetrically: loss reduces the higher pressure and increases the lower pressure + if (mdotEst >= 0) // volume → pipe + { + volP_eff -= dynP; + pipeP_eff += dynP; + } + else // pipe → volume + { + pipeP_eff -= dynP; + volP_eff += dynP; + } + } + + // --- Final nozzle solution with corrected pressures --- + float mdotSS, rhoFace0, uFace0, pFace0; + if (volP_eff >= pipeP_eff) + { + IsentropicOrifice.Compute(volP_eff, volRho, volT, pipeP_eff, gamma, R, area, Cd, + out float mUp, out rhoFace0, out uFace0, out pFace0); + mdotSS = mUp; + } + else + { + IsentropicOrifice.Compute(pipeP_eff, pipeRho, pipeT, volP_eff, gamma, R, area, Cd, + out float mUp, out rhoFace0, out uFace0, out pFace0); + mdotSS = -mUp; + } + + float mdot; + if (d.UseInertance) + { + float rhoUp = d.CurrentMdot >= 0 ? volRho : pipeRho; + float inertance = rhoUp * d.EffectiveLength / area; + float dp = volP_eff - pipeP_eff; + float resistance = MathF.Abs(dp) / MathF.Max(MathF.Abs(mdotSS), 1e-12f); + float dmdot_dt = (dp - resistance * d.CurrentMdot) / inertance; + mdot = d.CurrentMdot + dmdot_dt * dt; + + if (d.VolumePort.Owner is Volume0D vol0) + { + float maxOut = vol0.Mass / dt; + if (mdot > maxOut) mdot = maxOut; + } + if (float.IsNaN(mdot)) mdot = 0f; + } + else + { + mdot = mdotSS; + if (d.VolumePort.Owner is Volume0D vol0) + { + float maxOut = vol0.Mass / dt; + if (mdot > maxOut) mdot = maxOut; + } + } + + d.CurrentMdot = mdot; // stored for future steps (inertance or loss) + + // Ghost state construction + float rhoFace = mdot >= 0 ? volRho : pipeRho; + float pFace = pFace0; + float uFace = MathF.Abs(mdot) / MathF.Max(rhoFace * area, 1e-12f); + float airFracGhost; + if (mdot >= 0) + airFracGhost = volAF; + else + { + airFracGhost = pipeAF; + d.VolumePort.AirFraction = pipeAF; + } + + if (mdot >= 0 && d.IsLeftEnd) uFace = +uFace; + else if (mdot >= 0 && !d.IsLeftEnd) uFace = -uFace; + else if (mdot < 0 && d.IsLeftEnd) uFace = -uFace; + else if (mdot < 0 && !d.IsLeftEnd) uFace = +uFace; + + if (d.IsLeftEnd) + _pipeSystem.SetGhostLeft(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost); + else + _pipeSystem.SetGhostRight(d.PipeIndex, rhoFace, uFace, pFace, airFracGhost); + + d.VolumePort.MassFlowRate = -mdot; + if (-mdot >= 0) + { + float pipeH = gamma / (gamma - 1f) * pipeP / MathF.Max(pipeRho, 1e-12f); + d.VolumePort.SpecificEnthalpy = pipeH; + } + else + { + d.VolumePort.SpecificEnthalpy = volH; + } + } + } + + public void ResolveOpenEnds(float dt) + { + for (int i = 0; i < OpenEndCount; i++) + { + ref var d = ref _openEnds[i]; + + var (rhoInt, uInt, pInt) = d.IsLeftEnd + ? _pipeSystem.GetInteriorStateLeft(d.PipeIndex) + : _pipeSystem.GetInteriorStateRight(d.PipeIndex); + float afInt = d.IsLeftEnd + ? _pipeSystem.GetInteriorAirFractionLeft(d.PipeIndex) + : _pipeSystem.GetInteriorAirFractionRight(d.PipeIndex); + + float gamma = d.Gamma; + float gm1 = gamma - 1f; + float cInt = MathF.Sqrt(gamma * pInt / MathF.Max(rhoInt, 1e-12f)); + float pAmb = d.AmbientPressure; + + float Jplus = uInt + 2f * cInt / gm1; + float Jminus = uInt - 2f * cInt / gm1; + float s = pInt / MathF.Pow(rhoInt, gamma); + float rhoIso = MathF.Pow(pAmb / s, 1f / gamma); + float cIso = MathF.Sqrt(gamma * pAmb / MathF.Max(rhoIso, 1e-12f)); + float uIso = d.IsLeftEnd + ? (Jminus + 2f * cIso / gm1) + : (Jplus - 2f * cIso / gm1); + + bool supersonic = d.IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt); + float rhoGhost, uGhost, pGhost, afGhost; + if (supersonic) + { + rhoGhost = rhoInt; uGhost = uInt; pGhost = pInt; afGhost = afInt; + } + else + { + rhoGhost = rhoIso; uGhost = uIso; pGhost = pAmb; + bool inflow = d.IsLeftEnd ? (uIso >= 0f) : (uIso <= 0f); + afGhost = inflow ? 1f : afInt; + } + + if (d.IsLeftEnd) + _pipeSystem.SetGhostLeft(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost); + else + _pipeSystem.SetGhostRight(d.PipeIndex, rhoGhost, uGhost, pGhost, afGhost); + + float area = d.PipeArea; + float mdot = rhoGhost * uGhost * area; + if (d.IsLeftEnd) mdot = -mdot; + d.LastMassFlowRate = mdot; + d.LastFacePressure = pGhost; + } + } + } +} \ No newline at end of file diff --git a/Core/Constants.cs b/Core/Constants.cs index ca5d277..a0d27ff 100644 --- a/Core/Constants.cs +++ b/Core/Constants.cs @@ -2,10 +2,10 @@ namespace FluidSim.Core { public static class Constants { - public const double Gamma = 1.4; - public const double R_gas = 287.0; // J/(kg·K) - public const double P_amb = 101325.0; // Pa - public const double T_amb = 300.0; // K - public static readonly double Rho_amb = P_amb / (R_gas * T_amb); // ≈ 1.177 kg/m³ + public const float Gamma = 1.4f; + public const float R_gas = 287f; + public const float P_amb = 101325f; + public const float T_amb = 300f; + public static readonly float Rho_amb = P_amb / (R_gas * T_amb); } } \ No newline at end of file diff --git a/Core/GhostBuffer.cs b/Core/GhostBuffer.cs new file mode 100644 index 0000000..cdb55d9 --- /dev/null +++ b/Core/GhostBuffer.cs @@ -0,0 +1,27 @@ +namespace FluidSim.Core +{ + public class GhostBuffer + { + public float[] Rho, U, P, Y; + public int PipeCount { get; } + + public GhostBuffer(int pipeCount) + { + PipeCount = pipeCount; + int size = pipeCount * 2; + Rho = new float[size]; + U = new float[size]; + P = new float[size]; + Y = new float[size]; + } + + public void Set(int pipeIndex, bool isLeftEnd, float rho, float u, float p, float y) + { + int idx = pipeIndex * 2 + (isLeftEnd ? 0 : 1); + Rho[idx] = rho; + U[idx] = u; + P[idx] = p; + Y[idx] = y; + } + } +} \ No newline at end of file diff --git a/Core/IsentropicOrifice.cs b/Core/IsentropicOrifice.cs index ce23e83..b14e584 100644 --- a/Core/IsentropicOrifice.cs +++ b/Core/IsentropicOrifice.cs @@ -2,40 +2,30 @@ using System; namespace FluidSim.Core { - /// - /// Compressible flow through an orifice, modelled as an isentropic nozzle. - /// The caller provides the upstream stagnation state (pUp, rhoUp, TUp), - /// downstream pressure, orifice area, discharge coefficient, and gas properties. - /// Returns the face state and mass flow from upstream to downstream. - /// public static class IsentropicOrifice { public static void Compute( - double pUp, double rhoUp, double TUp, // upstream stagnation - double pDown, // downstream back pressure - double gamma, double R, double area, double Cd, - out double mdot, out double rhoFace, out double uFace, out double pFace) + float pUp, float rhoUp, float TUp, + float pDown, float gamma, float R, float area, float Cd, + out float mdot, out float rhoFace, out float uFace, out float pFace) { - mdot = 0; rhoFace = rhoUp; uFace = 0; pFace = pUp; + mdot = 0f; rhoFace = rhoUp; uFace = 0f; pFace = pUp; + if (area <= 0f || pUp <= 0f || rhoUp <= 0f || TUp <= 0f) return; - if (area <= 0 || pUp <= 0 || rhoUp <= 0 || TUp <= 0) - return; + float pr = MathF.Min(pDown / pUp, 1f); + if (pr < 1e-6f) pr = 1e-6f; + float prCrit = MathF.Pow(2f / (gamma + 1f), gamma / (gamma - 1f)); + if (pr < prCrit) pr = prCrit; - double pr = pDown / pUp; - if (pr < 1e-6) pr = 1e-6; + float exponent = (gamma - 1f) / gamma; + float M = MathF.Sqrt((2f / (gamma - 1f)) * (MathF.Pow(pr, -exponent) - 1f)); + if (float.IsNaN(M)) M = 0f; - double prCrit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0)); - if (pr < prCrit) pr = prCrit; // choked flow - - double exponent = (gamma - 1.0) / gamma; - double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -exponent) - 1.0)); - if (double.IsNaN(M)) M = 0; - - double aUp = Math.Sqrt(gamma * R * TUp); + float aUp = MathF.Sqrt(gamma * R * TUp); uFace = M * aUp; - rhoFace = rhoUp * Math.Pow(pr, 1.0 / gamma); + rhoFace = rhoUp * MathF.Pow(pr, 1f / gamma); pFace = pUp * pr; - mdot = rhoFace * uFace * area * Cd; // positive from upstream to downstream + mdot = rhoFace * uFace * area * Cd; } } } \ No newline at end of file diff --git a/Core/OpenEndLink.cs b/Core/OpenEndLink.cs deleted file mode 100644 index ba57d13..0000000 --- a/Core/OpenEndLink.cs +++ /dev/null @@ -1,120 +0,0 @@ -using System; -using FluidSim.Components; - -namespace FluidSim.Core -{ - /// - /// Characteristic open‑end boundary condition after Jones (1978). - /// For all subsonic flow (outflow and inflow), the ghost state is derived - /// from the isentropic expansion to ambient pressure, using the pipe's entropy, - /// and the outgoing Riemann invariant. This avoids a density jump at flow reversal. - /// Supersonic outflow extrapolates the interior state. - /// Now includes air fraction tracking: incoming air is fresh (AF=1), outgoing uses interior pipe AF. - /// - public class OpenEndLink - { - public Pipe1D Pipe { get; } - public bool IsLeftEnd { get; } - public double AmbientPressure { get; set; } = 101325.0; - public double Gamma { get; set; } = 1.4; - - public double LastMassFlowRate { get; private set; } - public double LastFaceDensity { get; private set; } - public double LastFaceVelocity { get; private set; } - public double LastFacePressure { get; private set; } - - public OpenEndLink(Pipe1D pipe, bool isLeftEnd) - { - Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe)); - IsLeftEnd = isLeftEnd; - } - - public void Resolve(double dtSub) - { - (double rhoInt, double uInt, double pInt) = IsLeftEnd - ? Pipe.GetInteriorStateLeft() - : Pipe.GetInteriorStateRight(); - - double airFracInt = IsLeftEnd - ? Pipe.GetInteriorAirFractionLeft() - : Pipe.GetInteriorAirFractionRight(); - - double gamma = Gamma; - double gm1 = gamma - 1.0; - double cInt = Math.Sqrt(gamma * pInt / Math.Max(rhoInt, 1e-12)); - double pAmb = AmbientPressure; - - // Riemann invariants - double J_plus = uInt + 2.0 * cInt / gm1; - double J_minus = uInt - 2.0 * cInt / gm1; - - double rhoGhost, uGhost, pGhost, airFracGhost; - - // ---- Subsonic branch (used for both outflow and inflow) ---- - double s = pInt / Math.Pow(rhoInt, gamma); // entropy constant - double rhoIso = Math.Pow(pAmb / s, 1.0 / gamma); - double cIso = Math.Sqrt(gamma * pAmb / Math.Max(rhoIso, 1e-12)); - double uIso = IsLeftEnd - ? (J_minus + 2.0 * cIso / gm1) - : (J_plus - 2.0 * cIso / gm1); - - // Check for supersonic outflow - bool supersonic = IsLeftEnd - ? (uInt <= -cInt) - : (uInt >= cInt); - - if (!supersonic) - { - if (IsLeftEnd) - supersonic = uIso <= -cIso; - else - supersonic = uIso >= cIso; - } - - if (supersonic) - { - // Supersonic outflow – extrapolate interior - rhoGhost = rhoInt; - uGhost = uInt; - pGhost = pInt; - airFracGhost = airFracInt; // whatever is leaving - } - else - { - // Subsonic flow – use isentropic state to ambient - rhoGhost = rhoIso; - uGhost = uIso; - pGhost = pAmb; - - // Determine if inflow or outflow - bool isInflow = IsLeftEnd ? (uIso >= 0) : (uIso <= 0); // positive u means into pipe from left end? Wait: left end: u>0 means flow to the right, into pipe. Right end: u>0 means flow to the right, out of pipe. Let's use mass flow sign later. - // More straightforward: if using the isentropic state, the ghost velocity direction indicates flow. For inflow (ambient to pipe), airFraction = 1.0; for outflow, airFraction = interior's AF. - if ((IsLeftEnd && uIso >= 0) || (!IsLeftEnd && uIso <= 0)) - { - // Inflow (ambient enters pipe) - airFracGhost = 1.0; - } - else - { - // Outflow (pipe exits to ambient) - airFracGhost = airFracInt; - } - } - - // Apply ghost to pipe - if (IsLeftEnd) - Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost, airFracGhost); - else - Pipe.SetGhostRight(rhoGhost, uGhost, pGhost, airFracGhost); - - // Mass flow out of the pipe (positive = leaving) - double area = Pipe.Area; - double mdot = rhoGhost * uGhost * area; - if (IsLeftEnd) mdot = -mdot; // left end: positive u is into pipe, outward flow is -u - LastMassFlowRate = mdot; - LastFaceDensity = rhoGhost; - LastFaceVelocity = uGhost; - LastFacePressure = pGhost; - } - } -} \ No newline at end of file diff --git a/Core/OrificeLink.cs b/Core/OrificeLink.cs deleted file mode 100644 index a416e6f..0000000 --- a/Core/OrificeLink.cs +++ /dev/null @@ -1,173 +0,0 @@ -using System; -using FluidSim.Components; -using FluidSim.Interfaces; - -namespace FluidSim.Core -{ - public class OrificeLink - { - public Port? VolumePort { get; } - public Pipe1D Pipe { get; } - public bool IsPipeLeftEnd { get; } - public Func AreaProvider { get; set; } - public double DischargeCoefficient { get; set; } = 0.62; - - public double EffectiveLength { get; set; } = 0.001; - public bool UseInertance { get; set; } = false; - - // Current mass flow (kg/s, positive = volume → pipe) - private double _mdot; - - public double LastMassFlowRate { get; private set; } // positive = into volume - public double LastFaceDensity { get; private set; } - public double LastFaceVelocity { get; private set; } - public double LastFacePressure { get; private set; } - - public OrificeLink(Port? volumePort, Pipe1D pipe, bool isPipeLeftEnd, Func areaProvider) - { - VolumePort = volumePort; - Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe)); - IsPipeLeftEnd = isPipeLeftEnd; - AreaProvider = areaProvider ?? throw new ArgumentNullException(nameof(areaProvider)); - _mdot = 0.0; - } - - public void Resolve(double dtSub) - { - double area = AreaProvider(); - if (area < 1e-12 || VolumePort == null) - { - SetClosedWall(); - return; - } - - // Gather states - double volP = VolumePort.Pressure; - double volRho = VolumePort.Density; - double volT = VolumePort.Temperature; - double volH = VolumePort.SpecificEnthalpy; - double volAF = VolumePort.AirFraction; - - (double pipeRho, double pipeU, double pipeP) = IsPipeLeftEnd - ? Pipe.GetInteriorStateLeft() - : Pipe.GetInteriorStateRight(); - - double pipeT = pipeP / Math.Max(pipeRho * 287.0, 1e-12); - double pipeAF = IsPipeLeftEnd - ? Pipe.GetInteriorAirFractionLeft() - : Pipe.GetInteriorAirFractionRight(); - - double gamma = 1.4; - double R = 287.0; - - // ---- Steady‑state nozzle solution ---- - double mdotSS; // positive = volume → pipe - double rhoFace0, uFace0, pFace0; - if (volP >= pipeP) - { - IsentropicOrifice.Compute(volP, volRho, volT, pipeP, gamma, R, area, DischargeCoefficient, - out double mdotUpToDown, out rhoFace0, out uFace0, out pFace0); - mdotSS = mdotUpToDown; - } - else - { - IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, gamma, R, area, DischargeCoefficient, - out double mdotUpToDown, out rhoFace0, out uFace0, out pFace0); - mdotSS = -mdotUpToDown; - } - - // ---- Dynamic update ---- - if (UseInertance) - { - double rhoUp = _mdot >= 0 ? volRho : pipeRho; - double inertance = rhoUp * EffectiveLength / area; - double dp = volP - pipeP; - double resistance = Math.Abs(dp) / Math.Max(Math.Abs(mdotSS), 1e-12); - double dmdot_dt = (dp - resistance * _mdot) / inertance; - _mdot += dmdot_dt * dtSub; - } - else - { - _mdot = mdotSS; - } - - // Clamp outflow to available mass (if finite volume) - if (VolumePort.Owner is Volume0D vol) - { - double maxOut = vol.Mass / dtSub; - if (_mdot > maxOut) _mdot = maxOut; - } - - // ---- Ghost state with air fraction ---- - double rhoFace = _mdot >= 0 ? volRho : pipeRho; - double pFace = pFace0; - double mdotMag = Math.Abs(_mdot); - double uFace = mdotMag / (rhoFace * area); - - // Determine air fraction for ghost and for volume port - double airFracGhost; // air fraction of ghost cell (at pipe end) - double airFracForVolume; // if flow reverses into volume, this is the air fraction entering volume - - if (_mdot >= 0) // volume → pipe - { - airFracGhost = volAF; - // Flow enters pipe; no need to set volume's air fraction (port already has its own) - airFracForVolume = volAF; // unused - } - else // pipe → volume - { - airFracGhost = pipeAF; - airFracForVolume = pipeAF; - VolumePort.AirFraction = airFracForVolume; - } - - if (IsPipeLeftEnd) - uFace = _mdot >= 0 ? uFace : -uFace; - else - uFace = _mdot >= 0 ? -uFace : uFace; - - if (IsPipeLeftEnd) - Pipe.SetGhostLeft(rhoFace, uFace, pFace, airFracGhost); - else - Pipe.SetGhostRight(rhoFace, uFace, pFace, airFracGhost); - - // Store results (positive = into volume) - LastMassFlowRate = -_mdot; - LastFaceDensity = rhoFace; - LastFaceVelocity = uFace; - LastFacePressure = pFace; - - VolumePort.MassFlowRate = -_mdot; - - // Enthalpy transport - if (-_mdot >= 0) // inflow → pipe enthalpy - { - double hPipe = gamma / (gamma - 1.0) * pipeP / Math.Max(pipeRho, 1e-12); - VolumePort.SpecificEnthalpy = hPipe; - } - else - { - VolumePort.SpecificEnthalpy = volH; - } - } - - private void SetClosedWall() - { - var (rInt, uInt, pInt) = IsPipeLeftEnd - ? Pipe.GetInteriorStateLeft() - : Pipe.GetInteriorStateRight(); - - if (IsPipeLeftEnd) - Pipe.SetGhostLeft(rInt, -uInt, pInt, IsPipeLeftEnd ? Pipe.GetInteriorAirFractionLeft() : Pipe.GetInteriorAirFractionRight()); - else - Pipe.SetGhostRight(rInt, -uInt, pInt, IsPipeLeftEnd ? Pipe.GetInteriorAirFractionLeft() : Pipe.GetInteriorAirFractionRight()); - - LastMassFlowRate = 0.0; - LastFaceDensity = rInt; - LastFaceVelocity = 0.0; - LastFacePressure = pInt; - if (VolumePort != null) - VolumePort.MassFlowRate = 0.0; - } - } -} \ No newline at end of file diff --git a/Core/Pipesystem.cs b/Core/Pipesystem.cs new file mode 100644 index 0000000..349ebec --- /dev/null +++ b/Core/Pipesystem.cs @@ -0,0 +1,560 @@ +using System; +using System.Diagnostics; +using System.Numerics; + +namespace FluidSim.Core +{ + public class PipeSystem + { + // ---------- Master arrays ---------- + private float[] _rho, _rhou, _E, _Y; + private readonly float[] _area; + private readonly float[] _dx; + private readonly int[] _pipeStart; + private readonly int[] _pipeEnd; + private readonly int _totalCells; // original cell count (visible) + private readonly int _allCells; // total allocated (padded to Vector.Count) + private readonly int _pipeCount; + + // Derived state – _p is kept for visualization, _c is gone + private float[] _p; + + // Flux arrays (size = _allCells + 1) + private float[] _fluxM, _fluxP, _fluxE, _fluxY; + + // Damping and relaxation (computed on‑the‑fly only if used) + private float[] _dampingFactors; + private float[] _relaxFactors; + private bool _applyDamping; + private bool _applyRelax; + + // Ghost buffer + private readonly GhostBuffer _ghost; + + // Wall mask – precomputed once + private readonly bool[] _isWallFace; + + // ---------- Physical constants ---------- + private const float Gamma = 1.4f; + private const float Gm1 = 0.4f; + private const float Gm1Inv = 1f / Gm1; // 2.5 + private const float GammaOverGm1 = Gamma / Gm1; // 3.5 + private float _coeffBase; + private float _relaxRate; + private float _ambientPressure = 101325f; + private float _ambientEnergyRef; + + public float DampingMultiplier + { + set + { + _coeffBase = 0.1f * value; + _applyDamping = _coeffBase != 0f; + } + } + public float EnergyRelaxationRate + { + set + { + _relaxRate = value; + _applyRelax = _relaxRate != 0f; + } + } + public float AmbientPressure + { + set + { + _ambientPressure = value; + _ambientEnergyRef = value * Gm1Inv; + } + } + + // ---------- Profiling ---------- + public bool EnableProfiling { get; set; } + private long _profFluxTicks; + private long _profUpdateTicks; + private long _profCallCount; + + // ---------- Construction ---------- + public PipeSystem(int totalCells, int[] pipeStart, int[] pipeEnd, + float[] area, float[] dx, + float initialRho, float initialU, float initialP) + { + _pipeStart = pipeStart; + _pipeEnd = pipeEnd; + _pipeCount = pipeStart.Length; + _totalCells = totalCells; + _area = area; + _dx = dx; + + // Pad to SIMD width so all vectorized loops cover the whole data + int vecSize = Vector.Count; + _allCells = totalCells % vecSize == 0 ? totalCells : totalCells + vecSize - (totalCells % vecSize); + + _rho = new float[_allCells]; + _rhou = new float[_allCells]; + _E = new float[_allCells]; + _Y = new float[_allCells]; + _p = new float[_allCells]; // pressure for drawing + int faceCount = _allCells + 1; + _fluxM = new float[faceCount]; + _fluxP = new float[faceCount]; + _fluxE = new float[faceCount]; + _fluxY = new float[faceCount]; + + _dampingFactors = new float[_allCells]; + _relaxFactors = new float[_allCells]; + _applyDamping = _coeffBase != 0f; + _applyRelax = _relaxRate != 0f; + + _ghost = new GhostBuffer(_pipeCount); + _ambientEnergyRef = initialP * Gm1Inv; + + // Pre‑compute wall face flags: each face that sits between two different pipes is a wall + _isWallFace = new bool[faceCount]; + for (int f = 1; f < _totalCells; f++) + { + for (int p = 0; p < _pipeCount; p++) + { + if (f == _pipeEnd[p] && f < _totalCells) + { + _isWallFace[f] = true; + break; + } + } + } + + // Initialize uniform state + float initE = initialP / (Gm1 * initialRho); + float rhoE = initialRho * initE + 0.5f * initialRho * initialU * initialU; + for (int i = 0; i < totalCells; i++) + { + _rho[i] = initialRho; + _rhou[i] = initialRho * initialU; + _E[i] = rhoE; + _Y[i] = 1f; + } + } + + // ---------- Ghost setters (for BoundarySystem) ---------- + public void SetGhostLeft(int pipeIndex, float rho, float u, float p, float y) + => _ghost.Set(pipeIndex, true, rho, u, p, y); + public void SetGhostRight(int pipeIndex, float rho, float u, float p, float y) + => _ghost.Set(pipeIndex, false, rho, u, p, y); + + // ---------- Public read methods ---------- + public int TotalCells => _totalCells; + public int PipeCount => _pipeCount; + public int GetPipeStart(int pipeIdx) => _pipeStart[pipeIdx]; + public int GetPipeEnd(int pipeIdx) => _pipeEnd[pipeIdx]; + public float GetCellPressure(int i) => _p[i]; + public float GetCellDensity(int i) => _rho[i]; + public float GetCellVelocity(int i) + { + float rho = _rho[i]; + return rho > 1e-12f ? _rhou[i] / rho : 0f; + } + public float GetCellAirFraction(int i) => _Y[i]; + + public (float rho, float u, float p) GetInteriorStateLeft(int pipeIdx) + { + int i = _pipeStart[pipeIdx]; + float rho = _rho[i]; + float rhou = _rhou[i]; + float u = rhou / MathF.Max(rho, 1e-12f); + float p = Gm1 * (_E[i] - 0.5f * rhou * u); + return (rho, u, p); + } + public (float rho, float u, float p) GetInteriorStateRight(int pipeIdx) + { + int i = _pipeEnd[pipeIdx] - 1; + float rho = _rho[i]; + float rhou = _rhou[i]; + float u = rhou / MathF.Max(rho, 1e-12f); + float p = Gm1 * (_E[i] - 0.5f * rhou * u); + return (rho, u, p); + } + public float GetInteriorAirFractionLeft(int pipeIdx) => _Y[_pipeStart[pipeIdx]]; + public float GetInteriorAirFractionRight(int pipeIdx) => _Y[_pipeEnd[pipeIdx] - 1]; + + public void SetCellState(int i, float rho, float u, float p, float y = 1f) + { + if (i < 0 || i >= _totalCells) return; + _rho[i] = rho; + _rhou[i] = rho * u; + _E[i] = p * Gm1Inv + 0.5f * rho * u * u; + _Y[i] = y; + } + + // ---------- Main step ---------- + public void SimulateStep(float dt) + { + long t0 = 0, t1 = 0; + if (EnableProfiling) + { + _profCallCount++; + t0 = Stopwatch.GetTimestamp(); + } + + ComputeFluxes(dt); + + if (EnableProfiling) + { + t1 = Stopwatch.GetTimestamp(); + _profFluxTicks += (t1 - t0); + t0 = t1; + } + + UpdateCells(dt); + + if (EnableProfiling) + { + t1 = Stopwatch.GetTimestamp(); + _profUpdateTicks += (t1 - t0); + } + } + + // ---------- Flux computation: fuses primitive calculation and flux evaluation ---------- + private void ComputeFluxes(float dt) + { + float fm, fp, fe; + int vecSize = Vector.Count; + + // ---- 1. Left ghost boundaries ---- + for (int p = 0; p < _pipeCount; p++) + { + int idx = _pipeStart[p]; + int ghostIdx = p * 2; + float rL = _ghost.Rho[ghostIdx]; + float uL = _ghost.U[ghostIdx]; + float pL = _ghost.P[ghostIdx]; + float YL = _ghost.Y[ghostIdx]; + float cL = MathF.Sqrt(Gamma * pL / MathF.Max(rL, 1e-12f)); + + float rR = _rho[idx], rhouR = _rhou[idx]; + float invRhoR = MathF.ReciprocalEstimate(MathF.Max(rR, 1e-12f)); + float uR = rhouR * invRhoR; + float pR = Gm1 * (_E[idx] - 0.5f * rhouR * uR); + float cR = MathF.Sqrt(Gamma * pR * invRhoR); + float YR = _Y[idx]; + + // store pressure for cell idx + _p[idx] = pR; + + LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe); + _fluxM[idx] = fm; _fluxP[idx] = fp; _fluxE[idx] = fe; + + float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR); + ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy); + _fluxY[idx] = fy; + } + + // ---- 2. Right ghost boundaries ---- + for (int p = 0; p < _pipeCount; p++) + { + int idx = _pipeEnd[p] - 1; + int face = idx + 1; + int ghostIdx = p * 2 + 1; + float rR = _ghost.Rho[ghostIdx]; + float uR = _ghost.U[ghostIdx]; + float pR = _ghost.P[ghostIdx]; + float YR = _ghost.Y[ghostIdx]; + float cR = MathF.Sqrt(Gamma * pR / MathF.Max(rR, 1e-12f)); + + float rL = _rho[idx], rhouL = _rhou[idx]; + float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f)); + float uL = rhouL * invRhoL; + float pL = Gm1 * (_E[idx] - 0.5f * rhouL * uL); + float cL = MathF.Sqrt(Gamma * pL * invRhoL); + float YL = _Y[idx]; + + // store pressure for cell idx + _p[idx] = pL; + + LaxFlux(rL, uL, pL, cL, rR, uR, pR, cR, out fm, out fp, out fe); + _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe; + + float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR); + ScalarFlux(rL, uL, YL, rR, uR, YR, alpha, out float fy); + _fluxY[face] = fy; + } + + // ---- 3. Interior faces – vectorised SIMD ---- + for (int face = 1; face < _totalCells; face++) + { + // Handle walls (rare) with scalar code + if (_isWallFace[face]) + { + int iL = face - 1; + float rL = _rho[iL], rhouL = _rhou[iL]; + float invRhoL = MathF.ReciprocalEstimate(MathF.Max(rL, 1e-12f)); + float uL = rhouL * invRhoL; + float pL = Gm1 * (_E[iL] - 0.5f * rhouL * uL); + float cL = MathF.Sqrt(Gamma * pL * invRhoL); + _p[iL] = pL; + + LaxFlux(rL, uL, pL, cL, rL, -uL, pL, cL, out fm, out fp, out fe); + _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe; + _fluxY[face] = 0f; + continue; + } + + // If the next vecSize faces contain a wall, fall back to scalar for this block + if (face + vecSize - 1 < _totalCells) + { + bool hasWall = false; + for (int f = face; f < face + vecSize; f++) + if (_isWallFace[f]) { hasWall = true; break; } + + if (!hasWall) + { + // --- Vectorised block --- + var rhoL = new Vector(_rho, face - 1); + var rhouL = new Vector(_rhou, face - 1); + var EL = new Vector(_E, face - 1); + var YL = new Vector(_Y, face - 1); + var rhoR = new Vector(_rho, face); + var rhouR = new Vector(_rhou, face); + var ER = new Vector(_E, face); + var YR = new Vector(_Y, face); + + var invRhoL = Vector.One / Vector.Max(rhoL, new Vector(1e-12f)); + var invRhoR = Vector.One / Vector.Max(rhoR, new Vector(1e-12f)); + var uL = rhouL * invRhoL; + var uR = rhouR * invRhoR; + var kinL = 0.5f * rhouL * uL; + var kinR = 0.5f * rhouR * uR; + var pL = Gm1 * (EL - kinL); + var pR = Gm1 * (ER - kinR); + var cL = Vector.SquareRoot(Gamma * pL * invRhoL); + var cR = Vector.SquareRoot(Gamma * pR * invRhoR); + + // Store pressures for visualisation (left cell of each face) + pL.CopyTo(_p, face - 1); + + // Lax‑Friedrichs fluxes + var ELs = pL * Gm1Inv * invRhoL + 0.5f * uL * uL; // energy per mass + var ERs = pR * Gm1Inv * invRhoR + 0.5f * uR * uR; + + var FmL = rhoL * uL; + var FpL = rhoL * uL * uL + pL; + var FeL = (rhoL * ELs + pL) * uL; + + var FmR = rhoR * uR; + var FpR = rhoR * uR * uR + pR; + var FeR = (rhoR * ERs + pR) * uR; + + var absUL = Vector.Abs(uL); + var absUR = Vector.Abs(uR); + var alpha = Vector.Max(absUL + cL, absUR + cR); + + var fmVec = 0.5f * (FmL + FmR) - 0.5f * alpha * (rhoR - rhoL); + var fpVec = 0.5f * (FpL + FpR) - 0.5f * alpha * (rhouR - rhouL); + var feVec = 0.5f * (FeL + FeR) - 0.5f * alpha * (rhoR * ERs - rhoL * ELs); + + var fyL = FmL * YL; + var fyR = FmR * YR; + var fyVec = 0.5f * (fyL + fyR) - 0.5f * alpha * (rhoR * YR - rhoL * YL); + + fmVec.CopyTo(_fluxM, face); + fpVec.CopyTo(_fluxP, face); + feVec.CopyTo(_fluxE, face); + fyVec.CopyTo(_fluxY, face); + + face += vecSize - 1; // loop increment will add 1, so we advance vecSize faces + continue; + } + } + + // --- Scalar interior face (fallback) --- + { + int iLf = face - 1, iRf = face; + float rLf = _rho[iLf], rhouLf = _rhou[iLf]; + float invRhoLf = MathF.ReciprocalEstimate(MathF.Max(rLf, 1e-12f)); + float uLf = rhouLf * invRhoLf; + float pLf = Gm1 * (_E[iLf] - 0.5f * rhouLf * uLf); + float cLf = MathF.Sqrt(Gamma * pLf * invRhoLf); + float YLf = _Y[iLf]; + _p[iLf] = pLf; + + float rRf = _rho[iRf], rhouRf = _rhou[iRf]; + float invRhoRf = MathF.ReciprocalEstimate(MathF.Max(rRf, 1e-12f)); + float uRf = rhouRf * invRhoRf; + float pRf = Gm1 * (_E[iRf] - 0.5f * rhouRf * uRf); + float cRf = MathF.Sqrt(Gamma * pRf * invRhoRf); + float YRf = _Y[iRf]; + + LaxFlux(rLf, uLf, pLf, cLf, rRf, uRf, pRf, cRf, out fm, out fp, out fe); + _fluxM[face] = fm; _fluxP[face] = fp; _fluxE[face] = fe; + + float alpha = MathF.Max(MathF.Abs(uLf) + cLf, MathF.Abs(uRf) + cRf); + ScalarFlux(rLf, uLf, YLf, rRf, uRf, YRf, alpha, out float fy); + _fluxY[face] = fy; + } + } + + // If damping/relaxation are active, compute the factors here (re-uses _dampingFactors/_relaxFactors arrays, + // but we no longer have a separate precompute pass). We compute them on demand in UpdateCells anyway? + // Actually UpdateCells multiplies by these factors; we can compute them there if needed. + } + + // ---------- Cell update (unchanged core, but skips relaxation/damping when not needed) ---------- + private void UpdateCells(float dt) + { + int vecSize = Vector.Count; + float dtRelax = -_relaxRate * dt; + + // Compute damping and relaxation factors if needed + if (_applyDamping) + { + for (int i = 0; i < _totalCells; i++) + { + float rho = _rho[i]; + _dampingFactors[i] = rho > 1e-12f + ? MathF.Exp(-_coeffBase * dt / rho) + : 1f; + } + } + if (_applyRelax) + { + var relaxVal = MathF.Exp(dtRelax); + for (int i = 0; i < _totalCells; i++) + _relaxFactors[i] = relaxVal; + } + + int iCell = 0; + for (; iCell <= _totalCells - vecSize; iCell += vecSize) + { + var rhoOld = new Vector(_rho, iCell); + var rhouOld = new Vector(_rhou, iCell); + var EOld = new Vector(_E, iCell); + var YOld = new Vector(_Y, iCell); + + var fluxM_L = new Vector(_fluxM, iCell); + var fluxP_L = new Vector(_fluxP, iCell); + var fluxE_L = new Vector(_fluxE, iCell); + var fluxY_L = new Vector(_fluxY, iCell); + + var fluxM_R = new Vector(_fluxM, iCell + 1); + var fluxP_R = new Vector(_fluxP, iCell + 1); + var fluxE_R = new Vector(_fluxE, iCell + 1); + var fluxY_R = new Vector(_fluxY, iCell + 1); + + var dtdx = new Vector(dt) / new Vector(_dx, iCell); + + var rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L); + var rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L); + var ENew = EOld - dtdx * (fluxE_R - fluxE_L); + var rhoYOld = rhoOld * YOld; + var rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L); + + if (_applyDamping) + rhouNew *= new Vector(_dampingFactors, iCell); + if (_applyRelax) + { + var ambRef = new Vector(_ambientEnergyRef); + var relax = new Vector(_relaxFactors, iCell); + ENew = ambRef + (ENew - ambRef) * relax; + } + + rhoNew = Vector.Max(rhoNew, new Vector(1e-12f)); + var kinNew = 0.5f * rhouNew * rhouNew / rhoNew; + var eMin = new Vector(100f * Gm1Inv) + kinNew; + ENew = Vector.Max(ENew, eMin); + + rhoNew.CopyTo(_rho, iCell); + rhouNew.CopyTo(_rhou, iCell); + ENew.CopyTo(_E, iCell); + var yNew = rhoYNew / rhoNew; + yNew = Vector.Min(Vector.Max(yNew, Vector.Zero), Vector.One); + yNew.CopyTo(_Y, iCell); + } + + // Scalar remainder (only a few cells) + for (; iCell < _totalCells; iCell++) + { + float rhoOld = _rho[iCell], rhouOld = _rhou[iCell], EOld = _E[iCell], YOld = _Y[iCell]; + float fluxM_L = _fluxM[iCell], fluxP_L = _fluxP[iCell], fluxE_L = _fluxE[iCell], fluxY_L = _fluxY[iCell]; + float fluxM_R = _fluxM[iCell + 1], fluxP_R = _fluxP[iCell + 1], fluxE_R = _fluxE[iCell + 1], fluxY_R = _fluxY[iCell + 1]; + float dtdx = dt / _dx[iCell]; + + float rhoNew = rhoOld - dtdx * (fluxM_R - fluxM_L); + float rhouNew = rhouOld - dtdx * (fluxP_R - fluxP_L); + float ENew = EOld - dtdx * (fluxE_R - fluxE_L); + float rhoYOld = rhoOld * YOld; + float rhoYNew = rhoYOld - dtdx * (fluxY_R - fluxY_L); + + if (_applyDamping) rhouNew *= _dampingFactors[iCell]; + if (_applyRelax) ENew = _ambientEnergyRef + (ENew - _ambientEnergyRef) * _relaxFactors[iCell]; + + rhoNew = MathF.Max(rhoNew, 1e-12f); + float kin = 0.5f * rhouNew * rhouNew / rhoNew; + float eMin = 100f * Gm1Inv + kin; + ENew = MathF.Max(ENew, eMin); + + _rho[iCell] = rhoNew; + _rhou[iCell] = rhouNew; + _E[iCell] = ENew; + _Y[iCell] = Math.Clamp(rhoYNew / rhoNew, 0f, 1f); + } + } + + // ---------- Scalar flux helpers (used in boundaries and scalar fallback) ---------- + private static void LaxFlux(float rL, float uL, float pL, float cL, + float rR, float uR, float pR, float cR, + out float fm, out float fp, out float fe) + { + float EL = pL * Gm1Inv / rL + 0.5f * uL * uL; + float ER = pR * Gm1Inv / rR + 0.5f * uR * uR; + float FmL = rL * uL; + float FpL = rL * uL * uL + pL; + float FeL = (rL * EL + pL) * uL; + float FmR = rR * uR; + float FpR = rR * uR * uR + pR; + float FeR = (rR * ER + pR) * uR; + float alpha = MathF.Max(MathF.Abs(uL) + cL, MathF.Abs(uR) + cR); + fm = 0.5f * (FmL + FmR) - 0.5f * alpha * (rR - rL); + fp = 0.5f * (FpL + FpR) - 0.5f * alpha * (rR * uR - rL * uL); + fe = 0.5f * (FeL + FeR) - 0.5f * alpha * (rR * ER - rL * EL); + } + + private static void ScalarFlux(float rL, float uL, float YL, + float rR, float uR, float YR, + float alpha, out float fy) + { + float FyL = rL * uL * YL; + float FyR = rR * uR * YR; + fy = 0.5f * (FyL + FyR) - 0.5f * alpha * (rR * YR - rL * YL); + } + + // ---------- Profiling report ---------- + public string GetProfileReport() + { + if (!EnableProfiling || _profCallCount == 0) + return "Pipe profiling disabled or no data."; + + double freq = Stopwatch.Frequency; + long totalTicks = _profFluxTicks + _profUpdateTicks; + if (totalTicks == 0) return "No pipe profile data collected."; + + double totalMs = totalTicks * 1000.0 / freq; + double avgCallUs = totalMs * 1000.0 / _profCallCount; + + double fluxMs = _profFluxTicks * 1000.0 / freq; + double updateMs = _profUpdateTicks * 1000.0 / freq; + + double fluxAvgUs = fluxMs * 1000.0 / _profCallCount; + double updateAvgUs = updateMs * 1000.0 / _profCallCount; + + string report = $" Pipe kernel (over {_profCallCount} calls, total {totalMs:F2} ms, avg {avgCallUs:F2} µs/call):\n"; + report += $" Fluxes (incl. primitives): {fluxMs:F2} ms ({_profFluxTicks * 100.0 / totalTicks:F1}%), avg {fluxAvgUs:F2} µs/call\n"; + report += $" Update cells: {updateMs:F2} ms ({_profUpdateTicks * 100.0 / totalTicks:F1}%), avg {updateAvgUs:F2} µs/call\n"; + + _profFluxTicks = 0; + _profUpdateTicks = 0; + _profCallCount = 0; + + return report; + } + } +} \ No newline at end of file diff --git a/Core/Solver.cs b/Core/Solver.cs index e6a42a4..1559fa2 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -10,136 +10,91 @@ namespace FluidSim.Core public class Solver { private readonly List _components = new(); - private readonly List _orificeLinks = new(); - private readonly List _openEndLinks = new(); - + private PipeSystem _pipeSystem; + private BoundarySystem _boundarySystem; private double _dt; - /// CFL target for sub‑stepping (0.3‑0.8). Lower values are safer for shocks. - public double CflTarget { get; set; } = 0.9; + public int SubStepCount { get; set; } = 4; + public bool EnableProfiling { get; set; } = false; - // ---------- Timing accumulators (reset every LogInterval steps) ---------- private long _stepCount; - private double _timeTotal, _timeCFL, _timeOrifice, _timeOpenEnd, - _timePipe, _timeClearGhosts, _timeUpdateState; - - private const int LogInterval = 5000; - private const bool EnableLogging = false; // temporarily ON for debugging + private long _ticksOrifice, _ticksOpenEnd, _ticksPipe, _ticksUpdate; public void SetTimeStep(double dt) => _dt = dt; - public void AddComponent(IComponent component) => _components.Add(component); - public void AddOrificeLink(OrificeLink link) => _orificeLinks.Add(link); - public void AddOpenEndLink(OpenEndLink link) => _openEndLinks.Add(link); + + public void SetPipeSystem(PipeSystem pipeSystem) + { + _pipeSystem = pipeSystem; + } + public void SetBoundarySystem(BoundarySystem boundarySystem) + { + _boundarySystem = boundarySystem; + } public void Step() { - var pipes = _components.OfType().ToList(); - if (pipes.Count == 0) return; + if (_pipeSystem == null || _boundarySystem == null) return; - var sw = Stopwatch.StartNew(); - - // CFL count – track which pipe demands the most sub‑steps - int nSub = 1; - Pipe1D worstPipe = pipes[0]; - foreach (var p in pipes) - { - int n = p.GetRequiredSubSteps(_dt, CflTarget); - if (n > nSub) - { - nSub = n; - worstPipe = p; - } - } - double dtSub = _dt / nSub; - - // ----- Diagnostic: warn if nSub is high ----- - if (nSub > 50) - { - double maxW = 0; - for (int i = 0; i < worstPipe.CellCount; i++) - { - double rho = worstPipe.GetCellDensity(i); - double u = Math.Abs(worstPipe.GetCellVelocity(i)); - double p = worstPipe.GetCellPressure(i); - double c = Math.Sqrt(1.4 * p / Math.Max(rho, 1e-12)); - if (u + c > maxW) maxW = u + c; - } - Console.WriteLine($"nSub = {nSub} (worst pipe: {worstPipe.Name}, maxW = {maxW:F0} m/s)"); - } - - _timeCFL += sw.Elapsed.TotalSeconds; - - // ----- Safety cap – prevent the solver from hanging ----- - const int maxSubSteps = 10000; - const int hardLimit = 500; // temporary low cap for debugging - - if (nSub > hardLimit) - { - Console.WriteLine($"nSub ({nSub}) exceeds hard limit {hardLimit}. Simulation step skipped."); - return; - } + int nSub = SubStepCount; + float dtSub = (float)(_dt / nSub); for (int sub = 0; sub < nSub; sub++) { - double t0; + long t0; - t0 = sw.Elapsed.TotalSeconds; - foreach (var link in _orificeLinks) - link.Resolve(dtSub); - _timeOrifice += sw.Elapsed.TotalSeconds - t0; + t0 = Stopwatch.GetTimestamp(); + _boundarySystem.ResolveOrifices(dtSub); + _ticksOrifice += Stopwatch.GetTimestamp() - t0; - t0 = sw.Elapsed.TotalSeconds; - foreach (var link in _openEndLinks) - link.Resolve(dtSub); - _timeOpenEnd += sw.Elapsed.TotalSeconds - t0; + t0 = Stopwatch.GetTimestamp(); + _boundarySystem.ResolveOpenEnds(dtSub); + _ticksOpenEnd += Stopwatch.GetTimestamp() - t0; - t0 = sw.Elapsed.TotalSeconds; - foreach (var p in pipes) - p.SimulateSingleStep(dtSub); - _timePipe += sw.Elapsed.TotalSeconds - t0; + t0 = Stopwatch.GetTimestamp(); + _pipeSystem.SimulateStep(dtSub); + _ticksPipe += Stopwatch.GetTimestamp() - t0; } - double tCG = sw.Elapsed.TotalSeconds; - foreach (var p in pipes) - p.ClearGhostFlags(); - _timeClearGhosts += sw.Elapsed.TotalSeconds - tCG; - - double tUS = sw.Elapsed.TotalSeconds; + long tUS = Stopwatch.GetTimestamp(); foreach (var comp in _components) - comp.UpdateState(_dt); - _timeUpdateState += sw.Elapsed.TotalSeconds - tUS; - - _timeTotal += sw.Elapsed.TotalSeconds; + comp.UpdateState((float)_dt); + _ticksUpdate += Stopwatch.GetTimestamp() - tUS; _stepCount++; - if (_stepCount % LogInterval == 0 && EnableLogging) + if (_stepCount % 5000 == 0 && EnableProfiling) { - if (_timeTotal > 0) - { - double stepsPerSec = LogInterval / _timeTotal; - double avgUs = (_timeTotal / LogInterval) * 1e6; + double freq = Stopwatch.Frequency; + double total = _ticksOrifice + _ticksOpenEnd + _ticksPipe + _ticksUpdate; + double avgStepUs = (total / freq) * 1e6 / 5000.0; - Console.WriteLine($"--- Solver timing ({LogInterval} steps) ---"); - Console.WriteLine($" Steps per second: {stepsPerSec:F1}"); - Console.WriteLine($" Avg step time: {avgUs:F1} µs (last nSub = {nSub})"); - Console.WriteLine($" CFL calc: {_timeCFL / _timeTotal * 100:F1} %"); - Console.WriteLine($" Sub‑step loop:"); - Console.WriteLine($" Orifice: {_timeOrifice / _timeTotal * 100:F1} %"); - Console.WriteLine($" OpenEnd: {_timeOpenEnd / _timeTotal * 100:F1} %"); - Console.WriteLine($" Pipe steps: {_timePipe / _timeTotal * 100:F1} %"); - Console.WriteLine($" Clear ghosts: {_timeClearGhosts / _timeTotal * 100:F1} %"); - Console.WriteLine($" Update state: {_timeUpdateState / _timeTotal * 100:F1} %"); - Console.WriteLine(); + int orificeCalls = 5000 * nSub; + int updateCalls = 5000; + + double orificeMs = _ticksOrifice * 1000.0 / freq; + double openEndMs = _ticksOpenEnd * 1000.0 / freq; + double pipeMs = _ticksPipe * 1000.0 / freq; + double updateMs = _ticksUpdate * 1000.0 / freq; + + double orificeAvgUs = orificeMs * 1000.0 / orificeCalls; + double openEndAvgUs = openEndMs * 1000.0 / orificeCalls; + double pipeAvgUs = pipeMs * 1000.0 / orificeCalls; + double updateAvgUs = updateMs * 1000.0 / updateCalls; + + Console.WriteLine($"--- Solver ({5000} steps, nSub={nSub}) ---"); + Console.WriteLine($" Average step: {avgStepUs:F2} µs"); + Console.WriteLine($" Orifice: {orificeMs:F2} ms ({(double)_ticksOrifice / total * 100:F1}%), avg {orificeAvgUs:F2} µs/call"); + Console.WriteLine($" OpenEnd: {openEndMs:F2} ms ({(double)_ticksOpenEnd / total * 100:F1}%), avg {openEndAvgUs:F2} µs/call"); + Console.WriteLine($" Pipe: {pipeMs:F2} ms ({(double)_ticksPipe / total * 100:F1}%), avg {pipeAvgUs:F2} µs/call"); + Console.WriteLine($" Update: {updateMs:F2} ms ({(double)_ticksUpdate / total * 100:F1}%), avg {updateAvgUs:F2} µs/call"); + + // Pipe internal breakdown (with per-phase averages) + if (_pipeSystem.EnableProfiling) + { + Console.WriteLine(_pipeSystem.GetProfileReport()); } - _timeTotal = 0; - _timeCFL = 0; - _timeOrifice = 0; - _timeOpenEnd = 0; - _timePipe = 0; - _timeClearGhosts = 0; - _timeUpdateState = 0; + _ticksOrifice = _ticksOpenEnd = _ticksPipe = _ticksUpdate = 0; } } } diff --git a/Core/SoundProcessor.cs b/Core/SoundProcessor.cs index a597d5b..978a50c 100644 --- a/Core/SoundProcessor.cs +++ b/Core/SoundProcessor.cs @@ -1,76 +1,34 @@ using System; -using FluidSim.Core; namespace FluidSim.Core { - /// - /// Synthesises far‑field exhaust sound using the monopole model - /// of Jones (1978). The radiated pressure is proportional to the - /// time derivative of the mass flow at the pipe exit. - /// - /// Reference: - /// Jones, A.D. (1978) "Noise characteristics and exhaust process - /// gas dynamics of a small 2-stroke engine", PhD thesis, Univ. Adelaide. - /// public class SoundProcessor { - private readonly double dt; - private readonly double r; // listener distance (m) - private readonly double scaleFactor; // 1 / (4π r) (free-field monopole) + private readonly float dt; + private readonly float scaleFactor; // 1 / (4π r) + private float flowLP, prevMassFlowOut, smoothDMdt; + private readonly float lpAlpha, alpha; - // ---------- Mass‑flow derivative (identical to original) ---------- - private double flowLP; - private readonly double lpAlpha; - private double prevMassFlowOut; - private double smoothDMdt; - private readonly double alpha; + public float Gain = 1f; - public float Gain { get; set; } = 1.0f; - - /// - /// - /// Audio sample rate (Hz). - /// Listener distance (m). - /// Ignored in this model; kept for compatibility. - public SoundProcessor(int sampleRate, - double listenerDistanceMeters = 1.0, - double pipeDiameterMeters = 0.0217) + public SoundProcessor(int sampleRate, float listenerDistance = 1f) { - dt = 1.0 / sampleRate; - r = listenerDistanceMeters; - scaleFactor = 1.0 / (4.0 * Math.PI * r); // free‑field monopole - - // ---- Smoothing time constants (unchanged) ---- - double tau = 0.02; // 2 ms for derivative - alpha = Math.Exp(-dt / tau); - - double tauLP = 0.00001; // 5 ms low‑pass on mass flow - lpAlpha = Math.Exp(-dt / tauLP); + dt = 1f / sampleRate; + scaleFactor = 1f / (4f * MathF.PI * listenerDistance); + float tau = 0.02f; + alpha = MathF.Exp(-dt / tau); + float tauLP = 0.005f; + lpAlpha = MathF.Exp(-dt / tauLP); } - /// - /// Process one sample. The OpenEndLink provides the instantaneous - /// exit‑plane mass flow. - /// - public float Process(OpenEndLink openEnd) + public float Process(float massFlowOut) { - double flowOut = openEnd.LastMassFlowRate; // kg/s, positive = leaving pipe - - // Low‑pass the mass flow signal - flowLP = lpAlpha * flowLP + (1.0 - lpAlpha) * flowOut; - - // Derivative of the smoothed mass flow - double rawDerivative = (flowLP - prevMassFlowOut) / dt; + flowLP = lpAlpha * flowLP + (1f - lpAlpha) * massFlowOut; + float rawDerivative = (flowLP - prevMassFlowOut) / dt; prevMassFlowOut = flowLP; - - // Smooth the derivative - smoothDMdt = alpha * smoothDMdt + (1.0 - alpha) * rawDerivative; - - // Far‑field monopole pressure (free‑field, Jones eq. 2.15 adapted) - double pressure = smoothDMdt * scaleFactor * Gain; - - // Soft clip to ±1 - return (float)pressure; + smoothDMdt = alpha * smoothDMdt + (1f - alpha) * rawDerivative; + float pressure = smoothDMdt * scaleFactor * Gain; + return MathF.Tanh(pressure); } } } \ No newline at end of file diff --git a/Interfaces/IComponent.cs b/Interfaces/IComponent.cs index b23be76..6f1ffed 100644 --- a/Interfaces/IComponent.cs +++ b/Interfaces/IComponent.cs @@ -2,18 +2,9 @@ using System.Collections.Generic; namespace FluidSim.Interfaces { - /// - /// Minimal interface for all simulation components that have ports. - /// public interface IComponent { - /// All ports exposed by this component. IReadOnlyList Ports { get; } - - /// - /// Called once per global time step to update the component's internal state - /// using the port flow data accumulated during sub‑steps. - /// - void UpdateState(double dt); + void UpdateState(float dt); } } \ No newline at end of file diff --git a/Interfaces/Port.cs b/Interfaces/Port.cs index 413c6b4..57f2d30 100644 --- a/Interfaces/Port.cs +++ b/Interfaces/Port.cs @@ -2,23 +2,23 @@ { public class Port { - public double MassFlowRate; // kg/s, positive INTO the component that owns this port - public double SpecificEnthalpy; // J/kg - public double Pressure; // Pa - public double Density; // kg/m³ - public double Temperature; // K - public double AirFraction; // mass fraction of air (0 = exhaust, 1 = air) + public float MassFlowRate; // kg/s, positive INTO owning component + public float SpecificEnthalpy; // J/kg + public float Pressure; // Pa + public float Density; // kg/m³ + public float Temperature; // K + public float AirFraction; // mass fraction (0 = exhaust, 1 = air) public object? Owner { get; set; } public Port() { - MassFlowRate = 0.0; - SpecificEnthalpy = 0.0; - Pressure = 101325.0; - Density = 1.225; - Temperature = 300.0; - AirFraction = 1.0; // default fresh air + MassFlowRate = 0f; + SpecificEnthalpy = 0f; + Pressure = 101325f; + Density = 1.225f; + Temperature = 300f; + AirFraction = 1f; } } } \ No newline at end of file diff --git a/Program.cs b/Program.cs index 95e639d..784dd3b 100644 --- a/Program.cs +++ b/Program.cs @@ -17,7 +17,7 @@ public class Program private const double DrawFrequency = 60.0; // Playback speed - private static double _desiredSpeed = 0.01; + private static double _desiredSpeed = 0.001; private static double _currentDisplaySpeed = _desiredSpeed; private const double MinSpeed = 0.0001; private const double MaxSpeed = 1.0; @@ -38,11 +38,11 @@ public class Program private static Text? _overlayText; // Throttle control - private static double _throttleTarget = 1.0; // 0‑1, set by arrow keys - private static double _throttleCurrent = 0.0; // actual current fraction (lerped) - private const double ThrottleLerpRate = 10.0; // times per second (speed of movement) + private static float _throttleTarget = 1.0f; // 0‑1, set by arrow keys + private static float _throttleCurrent = 0.0f; // actual current fraction (lerped) + private const float ThrottleLerpRate = 10.0f; // times per second (speed of movement) private static bool _wKeyHeld = false; - private static double _lastThrottleUpdateTime; + private static float _lastThrottleUpdateTime; private const int TargetMaxFill = (int)(SampleRate * 0.2); @@ -50,11 +50,11 @@ public class Program { var window = CreateWindow(); LoadFont(); - _scenario = new Inline4Scenario(); + _scenario = new HelmholtzScenario(); _scenario.Initialize(SampleRate); - _lastThrottleUpdateTime = 0.0; + _lastThrottleUpdateTime = 0.0f; - _simRingBuffer = new SimulationRingBuffer(131072); + _simRingBuffer = new SimulationRingBuffer(8192); _soundEngine = new SoundEngine(_simRingBuffer) { Volume = 100 }; _soundEngine.Start(); @@ -77,19 +77,19 @@ public class Program _soundEngine.Speed = _currentDisplaySpeed; // ---- Throttle update ---- - double dtThrottle = now - _lastThrottleUpdateTime; - _lastThrottleUpdateTime = now; + float dtThrottle = (float)now - _lastThrottleUpdateTime; + _lastThrottleUpdateTime = (float)now; - double throttleDesiredFraction = _wKeyHeld ? _throttleTarget : 0.0; + float throttleDesiredFraction = _wKeyHeld ? _throttleTarget : 0.0f; // Snap to zero instantly when target is zero (key released) if (throttleDesiredFraction == 0.0) { - _throttleCurrent = 0.0; + _throttleCurrent = 0.0f; } else { - double smoothing = 1.0 - Math.Exp(-ThrottleLerpRate * dtThrottle); + float smoothing = 1.0f - MathF.Exp(-ThrottleLerpRate * dtThrottle); _throttleCurrent += (throttleDesiredFraction - _throttleCurrent) * smoothing; } @@ -199,11 +199,11 @@ public class Program break; case Keyboard.Key.Up: - _throttleTarget = Math.Min(1.0, _throttleTarget + 0.05); + _throttleTarget = MathF.Min(1.0f, _throttleTarget + 0.05f); break; case Keyboard.Key.Down: - _throttleTarget = Math.Max(0.0, _throttleTarget - 0.05); + _throttleTarget = MathF.Max(0.0f, _throttleTarget - 0.05f); break; } } diff --git a/Scenarios/HelmholtzScenario.cs b/Scenarios/HelmholtzScenario.cs new file mode 100644 index 0000000..b62d2c7 --- /dev/null +++ b/Scenarios/HelmholtzScenario.cs @@ -0,0 +1,152 @@ +using FluidSim.Components; +using FluidSim.Core; +using FluidSim.Interfaces; +using SFML.Graphics; +using SFML.System; +using System; + +namespace FluidSim.Tests +{ + public class HelmholtzScenario : Scenario + { + private Volume0D cavity; + private Port cavityPort; + + private PipeSystem pipeSystem; + private int[] pipeStart = { 0 }; + private int[] pipeEnd; + + private BoundarySystem boundaries; + private int cavityOrificeIdx = 0; + private int openEndIdx = 0; + + private Solver solver; + private double dt; + private int stepCount; + + private SoundProcessor soundProcessor; + + public override void Initialize(int sampleRate) + { + dt = 1.0 / sampleRate; + + // --- Realistic Helmholtz resonator dimensions --- + float cavityVolume = 1e-3f; // 1 liter + float neckLength = 0.05f; // 5 cm + float neckDiameter = 0.02f; // 2 cm diameter + float neckArea = MathF.PI * 0.25f * neckDiameter * neckDiameter; + int neckCells = 20; + + // --- Volume (cavity) --- + float initialPressure = 1.2f * 101325f; // slight overpressure + float initialTemperature = 300f; + cavity = new Volume0D(cavityVolume, initialPressure, initialTemperature); + cavityPort = cavity.CreatePort(); + + // --- Pipe (neck) --- + float[] areas = new float[neckCells]; + float[] dxs = new float[neckCells]; + float dx = neckLength / neckCells; + for (int i = 0; i < neckCells; i++) + { + areas[i] = neckArea; + dxs[i] = dx; + } + pipeEnd = new[] { neckCells }; + + float rho0 = 101325f / (287f * 300f); + pipeSystem = new PipeSystem(neckCells, pipeStart, pipeEnd, areas, dxs, rho0, 0f, 101325f); + + + // Energy loss + cavity.EnergyRelaxationRate = 80f; + pipeSystem.EnergyRelaxationRate = 0f; + pipeSystem.DampingMultiplier = 2000f; + + // --- Boundary system --- + boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1); + + // Use steady orifice – the pipe already provides the inertia + boundaries.AddOrifice(cavityPort, pipeIndex: 0, isLeftEnd: true, areaIndex: cavityOrificeIdx, dischargeCoeff: 1f, lossCoefficient: 0.1f); + // LOSS COEFFICIENT BREAKS THE SYSTEM AT ~0.55, AT VALUES LOWER THAN THAT, IT SEEMS TO ONLY AFFECT VOLUME, NOT COMPOUND + + // Open end at right side of pipe + boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: false, 101325f, neckArea); + + float[] orificeAreas = new float[1] { neckArea }; + boundaries.SetOrificeAreas(orificeAreas); + + // --- Solver --- + // Slightly higher sub‑step count to ensure stability of the resonant oscillation + solver = new Solver { SubStepCount = 6, EnableProfiling = true }; + solver.SetTimeStep(dt); + solver.SetPipeSystem(pipeSystem); + solver.SetBoundarySystem(boundaries); + solver.AddComponent(cavity); + + // --- Sound --- + soundProcessor = new SoundProcessor(sampleRate, 1f) { Gain = 2f }; + + Console.WriteLine("Helmholtz resonator ready."); + stepCount = 0; + } + + public override float Process() + { + stepCount++; + if (stepCount <= 8192) return 0f; // let buffer pre‑fill + + solver.Step(); + + float flow = boundaries.GetOpenEndMassFlow(openEndIdx); + float sample = soundProcessor.Process(flow); + + if (stepCount % 10000 == 0) + { + float cavityP = cavity.Pressure; + float cavityT = cavity.Temperature; + float cavityRho = cavity.Density; + float cCavity = MathF.Sqrt(1.4f * cavityP / MathF.Max(cavityRho, 1e-12f)); + + // Temperature in the middle of the neck + int midCell = 10; + float pMid = pipeSystem.GetCellPressure(midCell); + float rhoMid = pipeSystem.GetCellDensity(midCell); + float tMid = pMid / MathF.Max(rhoMid * 287f, 1e-12f); + + // Neck effective length (physical + end correction) + float neckLen = 0.05f; // physical + float neckDia = 0.02f; + float neckArea = MathF.PI * 0.25f * neckDia * neckDia; + float endCorr = 0.85f * neckDia; // unflanged end + float L_eff = neckLen + endCorr; + + // Theoretical Helmholtz frequency from current cavity sound speed + float fHelmholtz = cCavity / (2f * MathF.PI) * + MathF.Sqrt(neckArea / (cavity.Volume * L_eff)); + + Console.WriteLine( + $"Step {stepCount}: cav P={cavityP / 1e5f:F4} bar, T={cavityT:F1} K, " + + $"pipeMid T={tMid:F1} K, est f={fHelmholtz:F1} Hz"); + } + + return sample; + } + + public override void Draw(RenderWindow target) + { + float winW = target.GetView().Size.X; + float winH = target.GetView().Size.Y; + + float cavityCenterX = 100f; + float cavityWidth = 80f, cavityHeight = 100f; + float cavityTopY = winH / 2f - cavityHeight / 2f; + DrawVolume(target, cavity, cavityCenterX, cavityTopY - 40f, cavityWidth, cavityHeight); + + float pipeStartX = cavityCenterX + cavityWidth / 2f + 10f; + float pipeEndX = winW - 50f; + float pipeCenterY = winH / 2f; + DrawPipe(target, pipeSystem, 0, pipeCenterY, pipeStartX, pipeEndX); + } + } +} \ No newline at end of file diff --git a/Scenarios/Inline4Scenario.cs b/Scenarios/Inline4Scenario.cs deleted file mode 100644 index 58531b9..0000000 --- a/Scenarios/Inline4Scenario.cs +++ /dev/null @@ -1,472 +0,0 @@ -using System; -using SFML.Graphics; -using SFML.System; -using FluidSim.Components; -using FluidSim.Core; -using FluidSim.Utils; - -namespace FluidSim.Tests -{ - public class Inline4Scenario : Scenario - { - // Crankshaft - private Crankshaft crankshaft; - - // Cylinders - private Cylinder cyl1, cyl2, cyl3, cyl4; - - // Intake - private Pipe1D intakePipeBeforeThrottle; - private Volume0D intakePlenum; - - // Runners (shorter, fewer cells) - private Pipe1D runner1, runner2, runner3, runner4; - - // Exhaust collector + tailpipe - private Volume0D exhaustCollector; - private Pipe1D tailPipe; - - // Exhaust stubs (short pipes between cylinders and collector) - private Pipe1D exhStub1, exhStub2, exhStub3, exhStub4; - - // Links – intake - private OpenEndLink intakeOpenEnd; - private OrificeLink throttleOrifice; - - // Plenum‑to‑runner orifices - private OrificeLink plenumToRunner1, plenumToRunner2, plenumToRunner3, plenumToRunner4; - - // Intake valves - private OrificeLink intakeValve1, intakeValve2, intakeValve3, intakeValve4; - - // Exhaust valves (cylinder → stub) - private OrificeLink exhaustValve1, exhaustValve2, exhaustValve3, exhaustValve4; - - // Stub‑to‑collector orifices - private OrificeLink stubToCollector1, stubToCollector2, stubToCollector3, stubToCollector4; - - // Collector‑to‑tailpipe orifice - private OrificeLink collectorToTailpipe; - - // Exhaust open end (tailpipe exit) - private OpenEndLink exhaustOpenEnd; - - private Solver solver; - private SoundProcessor exhaustSoundProcessor; - private SoundProcessor intakeSoundProcessor; - private OutdoorExhaustReverb reverb; - private double dt; - private int stepCount; - - public double MaxThrottleArea { get; set; } = 10 * Units.cm2; - - public override void Initialize(int sampleRate) - { - dt = 1.0 / sampleRate; - - solver = new Solver(); - solver.SetTimeStep(dt); - solver.CflTarget = 1; - - // ---- Shared crankshaft ---- - crankshaft = new Crankshaft(800); - crankshaft.Inertia = 1; - crankshaft.FrictionConstant = 3; - crankshaft.FrictionViscous = 0.2; - - // ---- Cylinder geometry ---- - double bore = 0.056, stroke = 0.057, conRod = 0.110, compRatio = 10; - double ivo = 350.0, ivc = 580.0, evo = 120.0, evc = 370.0; - - // Firing order 1-3-4-2 with 180° intervals (0°, 180°, 360°, 540°) - double phaseCyl1 = 0.0; - double phaseCyl3 = Math.PI; // 180° - double phaseCyl4 = 2.0 * Math.PI; // 360° - double phaseCyl2 = 3.0 * Math.PI; // 540° - - cyl1 = new Cylinder(bore, stroke, conRod, compRatio, ivo, ivc, evo, evc, crankshaft) - { - IntakeValveDiameter = 30 * Units.mm, - IntakeValveLift = 5 * Units.mm, - ExhaustValveDiameter = 28 * Units.mm, - ExhaustValveLift = 5 * Units.mm, - PhaseOffset = phaseCyl1, - EnergyVariationFraction = 0.03, - MisfireProbability = 0.0 - }; - cyl2 = new Cylinder(bore, stroke, conRod, compRatio, ivo, ivc, evo, evc, crankshaft) - { - IntakeValveDiameter = 30 * Units.mm, - IntakeValveLift = 5 * Units.mm, - ExhaustValveDiameter = 28 * Units.mm, - ExhaustValveLift = 5 * Units.mm, - PhaseOffset = phaseCyl2, - EnergyVariationFraction = 0.03, - MisfireProbability = 0.0 - }; - cyl3 = new Cylinder(bore, stroke, conRod, compRatio, ivo, ivc, evo, evc, crankshaft) - { - IntakeValveDiameter = 30 * Units.mm, - IntakeValveLift = 5 * Units.mm, - ExhaustValveDiameter = 28 * Units.mm, - ExhaustValveLift = 5 * Units.mm, - PhaseOffset = phaseCyl3, - EnergyVariationFraction = 0.03, - MisfireProbability = 0.0 - }; - cyl4 = new Cylinder(bore, stroke, conRod, compRatio, ivo, ivc, evo, evc, crankshaft) - { - IntakeValveDiameter = 30 * Units.mm, - IntakeValveLift = 5 * Units.mm, - ExhaustValveDiameter = 28 * Units.mm, - ExhaustValveLift = 5 * Units.mm, - PhaseOffset = phaseCyl4, - EnergyVariationFraction = 0.03, - MisfireProbability = 0.0 - }; - solver.AddComponent(cyl1); - solver.AddComponent(cyl2); - solver.AddComponent(cyl3); - solver.AddComponent(cyl4); - - double pipeDiameter = 4 * Units.cm; - double pipeArea = Units.AreaFromDiameter(pipeDiameter); - - // Sound processors (only one exhaust source now) - exhaustSoundProcessor = new SoundProcessor(sampleRate, 1, pipeDiameter) { Gain = 0.2f }; - intakeSoundProcessor = new SoundProcessor(sampleRate, 1, pipeDiameter) { Gain = 0.2f }; - reverb = new OutdoorExhaustReverb(sampleRate); - - // ---- Intake pipe before throttle (shorter, fewer cells) ---- - intakePipeBeforeThrottle = new Pipe1D(0.1, pipeArea, 10); - intakePipeBeforeThrottle.Name = "Intake pipe"; - solver.AddComponent(intakePipeBeforeThrottle); - - // ---- Plenum ---- - intakePlenum = new Volume0D(100 * Units.mL, 101325.0, 300.0); - var plenumInlet = intakePlenum.CreatePort(); // port 0 - var plenumOut1 = intakePlenum.CreatePort(); // port 1 - var plenumOut2 = intakePlenum.CreatePort(); // port 2 - var plenumOut3 = intakePlenum.CreatePort(); // port 3 - var plenumOut4 = intakePlenum.CreatePort(); // port 4 - solver.AddComponent(intakePlenum); - - // ---- Intake runners (shorter, fewer cells) ---- - runner1 = new Pipe1D(0.1, pipeArea, 5) { Name = "Runner 1" }; - runner2 = new Pipe1D(0.1, pipeArea, 5) { Name = "Runner 2" }; - runner3 = new Pipe1D(0.1, pipeArea, 5) { Name = "Runner 3" }; - runner4 = new Pipe1D(0.1, pipeArea, 5) { Name = "Runner 4" }; - solver.AddComponent(runner1); - solver.AddComponent(runner2); - solver.AddComponent(runner3); - solver.AddComponent(runner4); - - // ---- Exhaust collector volume ---- - exhaustCollector = new Volume0D(200 * Units.mL, 101325.0, 800.0); - var colIn1 = exhaustCollector.CreatePort(); // cylinder 1 stub - var colIn2 = exhaustCollector.CreatePort(); // cylinder 2 stub - var colIn3 = exhaustCollector.CreatePort(); // cylinder 3 stub - var colIn4 = exhaustCollector.CreatePort(); // cylinder 4 stub - var colOut = exhaustCollector.CreatePort(); // to tailpipe - solver.AddComponent(exhaustCollector); - - // ---- Exhaust stub pipes (short connection cylinder → collector) ---- - exhStub1 = new Pipe1D(0.1, pipeArea, 5) { Name = "ExhStub 1" }; - exhStub2 = new Pipe1D(0.1, pipeArea, 5) { Name = "ExhStub 2" }; - exhStub3 = new Pipe1D(0.1, pipeArea, 5) { Name = "ExhStub 3" }; - exhStub4 = new Pipe1D(0.1, pipeArea, 5) { Name = "ExhStub 4" }; - solver.AddComponent(exhStub1); - solver.AddComponent(exhStub2); - solver.AddComponent(exhStub3); - solver.AddComponent(exhStub4); - - foreach (var p in new[] { runner1, runner2, runner3, runner4, exhStub1, exhStub2, exhStub3, exhStub4, intakePipeBeforeThrottle }) - { - p.DampingMultiplier = 0.5; - p.EnergyRelaxationRate = 0.0; - } - - // ---- Tailpipe (single exhaust pipe) ---- - tailPipe = new Pipe1D(0.5, pipeArea, 20) - { - Name = "Tailpipe", - DampingMultiplier = 0.5, - EnergyRelaxationRate = 0.0 - }; - solver.AddComponent(tailPipe); - - // ---- Plenum → runner orifices (volume port to pipe left end) ---- - plenumToRunner1 = new OrificeLink(plenumOut1, runner1, isPipeLeftEnd: true, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - plenumToRunner2 = new OrificeLink(plenumOut2, runner2, isPipeLeftEnd: true, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - plenumToRunner3 = new OrificeLink(plenumOut3, runner3, isPipeLeftEnd: true, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - plenumToRunner4 = new OrificeLink(plenumOut4, runner4, isPipeLeftEnd: true, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - solver.AddOrificeLink(plenumToRunner1); - solver.AddOrificeLink(plenumToRunner2); - solver.AddOrificeLink(plenumToRunner3); - solver.AddOrificeLink(plenumToRunner4); - - // ---- Intake valves (cylinder port to runner right end) ---- - intakeValve1 = new OrificeLink(cyl1.IntakePort, runner1, isPipeLeftEnd: false, areaProvider: () => cyl1.IntakeValveArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - intakeValve2 = new OrificeLink(cyl2.IntakePort, runner2, isPipeLeftEnd: false, areaProvider: () => cyl2.IntakeValveArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - intakeValve3 = new OrificeLink(cyl3.IntakePort, runner3, isPipeLeftEnd: false, areaProvider: () => cyl3.IntakeValveArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - intakeValve4 = new OrificeLink(cyl4.IntakePort, runner4, isPipeLeftEnd: false, areaProvider: () => cyl4.IntakeValveArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - solver.AddOrificeLink(intakeValve1); - solver.AddOrificeLink(intakeValve2); - solver.AddOrificeLink(intakeValve3); - solver.AddOrificeLink(intakeValve4); - - // ---- Exhaust valves (cylinder port to stub left end) ---- - exhaustValve1 = new OrificeLink(cyl1.ExhaustPort, exhStub1, isPipeLeftEnd: true, areaProvider: () => cyl1.ExhaustValveArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - exhaustValve2 = new OrificeLink(cyl2.ExhaustPort, exhStub2, isPipeLeftEnd: true, areaProvider: () => cyl2.ExhaustValveArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - exhaustValve3 = new OrificeLink(cyl3.ExhaustPort, exhStub3, isPipeLeftEnd: true, areaProvider: () => cyl3.ExhaustValveArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - exhaustValve4 = new OrificeLink(cyl4.ExhaustPort, exhStub4, isPipeLeftEnd: true, areaProvider: () => cyl4.ExhaustValveArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - solver.AddOrificeLink(exhaustValve1); - solver.AddOrificeLink(exhaustValve2); - solver.AddOrificeLink(exhaustValve3); - solver.AddOrificeLink(exhaustValve4); - - // ---- Stub → collector orifices (collector port to stub right end) ---- - stubToCollector1 = new OrificeLink(colIn1, exhStub1, isPipeLeftEnd: false, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - stubToCollector2 = new OrificeLink(colIn2, exhStub2, isPipeLeftEnd: false, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - stubToCollector3 = new OrificeLink(colIn3, exhStub3, isPipeLeftEnd: false, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - stubToCollector4 = new OrificeLink(colIn4, exhStub4, isPipeLeftEnd: false, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - solver.AddOrificeLink(stubToCollector1); - solver.AddOrificeLink(stubToCollector2); - solver.AddOrificeLink(stubToCollector3); - solver.AddOrificeLink(stubToCollector4); - - // ---- Collector → tailpipe (collector port to tailpipe left end) ---- - collectorToTailpipe = new OrificeLink(colOut, tailPipe, isPipeLeftEnd: true, areaProvider: () => pipeArea) - { DischargeCoefficient = 1.0, UseInertance = false }; - solver.AddOrificeLink(collectorToTailpipe); - - // ---- Exhaust open end (tailpipe exit) ---- - exhaustOpenEnd = new OpenEndLink(tailPipe, isLeftEnd: false) - { - AmbientPressure = 101325.0, - Gamma = 1.4 - }; - solver.AddOpenEndLink(exhaustOpenEnd); - - // ---- Intake open end ---- - intakeOpenEnd = new OpenEndLink(intakePipeBeforeThrottle, isLeftEnd: true) - { - AmbientPressure = 101325.0, - Gamma = 1.4 - }; - solver.AddOpenEndLink(intakeOpenEnd); - - // ---- Throttle ---- - throttleOrifice = new OrificeLink(plenumInlet, intakePipeBeforeThrottle, isPipeLeftEnd: false, - areaProvider: () => MaxThrottleArea * Math.Clamp(Throttle, 0.0005, 1.0)) - { - DischargeCoefficient = 0.9, - UseInertance = false - }; - solver.AddOrificeLink(throttleOrifice); - - stepCount = 0; - Console.WriteLine("Inline-4 engine test"); - Console.WriteLine($"Bore {bore * 1000:F0}mm, Stroke {stroke * 1000:F0}mm, CR {compRatio}"); - Console.WriteLine("Firing order 1-3-4-2, 180° intervals"); - } - - public override float Process() - { - crankshaft.Step(dt); - - cyl1.PreStep(dt); - cyl2.PreStep(dt); - cyl3.PreStep(dt); - cyl4.PreStep(dt); - - solver.Step(); - stepCount++; - - if (stepCount % 10000 == 0) - { - double rpm = crankshaft.AngularVelocity * 60.0 / (2.0 * Math.PI); - Console.WriteLine($"Step {stepCount}, RPM = {rpm:F0}, " + - $"cyl1 P = {cyl1.Pressure / 1e5:F2} bar, " + - $"plenum P = {intakePlenum.Pressure / 1e5:F2} bar"); - } - - // Sound: only one exhaust source now - float exhaustSound = exhaustSoundProcessor.Process(exhaustOpenEnd); - float intakeSound = intakeSoundProcessor.Process(intakeOpenEnd); - return reverb.Process(exhaustSound * 0.25f + intakeSound); - } - - public override void Draw(RenderWindow target) - { - float winW = target.GetView().Size.X; - float winH = target.GetView().Size.Y; - - // --- Layout constants --- - float leftMargin = 40f; - float plenumW = 50f, plenumH = 120f; - float cylinderWidth = 60f, cylinderMaxHeight = 180f; - float cylinderSpacing = 90f; - float cylinderTopY = winH * 0.25f; - - // Plenum position - float plenumCenterX = leftMargin + plenumW / 2f; - float plenumTopY = cylinderTopY - 20f; - DrawVolume(target, intakePlenum, plenumCenterX, plenumTopY, plenumW, plenumH); - - // Throttle symbol (yellow rectangle) left of plenum - float throttleWidth = 8f, throttleHeight = 30f; - float throttleCenterX = leftMargin - 10f; - var throttleRect = new RectangleShape(new Vector2f(throttleWidth, throttleHeight)) - { - FillColor = Color.Yellow, - Position = new Vector2f(throttleCenterX - throttleWidth / 2f, plenumTopY + plenumH / 2f - throttleHeight / 2f) - }; - target.Draw(throttleRect); - - // Intake pipe before throttle (left of throttle) - float intakePipeEndX = throttleCenterX - throttleWidth / 2f; - float intakePipeStartX = intakePipeEndX - 100f; - float intakePipeY = plenumTopY + plenumH / 2f; - DrawPipe(target, intakePipeBeforeThrottle, intakePipeY, intakePipeStartX, intakePipeEndX); - - // Intake open end marker - var intakeMark = new CircleShape(4f) { FillColor = Color.Magenta }; - intakeMark.Position = new Vector2f(intakePipeStartX - 4f, intakePipeY - 4f); - target.Draw(intakeMark); - - // Cylinders and runners - float runnerStartX = leftMargin + plenumW; - Cylinder[] cyls = { cyl1, cyl2, cyl3, cyl4 }; - Pipe1D[] runners = { runner1, runner2, runner3, runner4 }; - - for (int i = 0; i < 4; i++) - { - float cylCenterX = runnerStartX + 40f + i * cylinderSpacing; - float runnerEndX = cylCenterX; - DrawPipe(target, runners[i], plenumTopY + plenumH / 2f, runnerStartX, runnerEndX); - DrawCylinder(target, cyls[i], cylCenterX, cylinderTopY, cylinderWidth, cylinderMaxHeight); - } - - // Exhaust collector below cylinders - float collectorLeftX = runnerStartX + 40f - cylinderWidth / 2f; - float collectorWidth = 3 * cylinderSpacing + cylinderWidth; - float collectorTopY = cylinderTopY + cylinderMaxHeight + 40f; - float collectorHeight = 50f; - float collectorCenterX = collectorLeftX + collectorWidth / 2f; - DrawVolume(target, exhaustCollector, collectorCenterX, collectorTopY, collectorWidth, collectorHeight); - - // Tailpipe from right edge of collector - float tailStartX = collectorLeftX + collectorWidth; - float tailEndX = tailStartX + 150f; - float tailCenterY = collectorTopY + collectorHeight / 2f; - DrawPipe(target, tailPipe, tailCenterY, tailStartX, tailEndX); - - // Exhaust open end marker - var exhaustMark = new CircleShape(4f) { FillColor = Color.Magenta }; - exhaustMark.Position = new Vector2f(tailEndX - 4f, tailCenterY - 4f); - target.Draw(exhaustMark); - - // Exhaust stubs (vertical connections from cylinder bottom to collector) - Pipe1D[] stubs = { exhStub1, exhStub2, exhStub3, exhStub4 }; - for (int i = 0; i < 4; i++) - { - float cylCenterX = runnerStartX + 40f + i * cylinderSpacing; - float vertStartY = cylinderTopY + cylinderMaxHeight; - float vertEndY = collectorTopY; - // Draw stub as a vertical pipe - DrawPipeVertical(target, stubs[i], cylCenterX, vertStartY, vertEndY); - } - } - - // Helper to draw a pipe vertically (reuse temperature coloring) - private void DrawPipeVertical(RenderWindow target, Pipe1D pipe, float centerX, float topY, float bottomY) - { - int n = pipe.CellCount; - if (n < 2) return; - - float pipeLengthPx = bottomY - topY; - float dy = pipeLengthPx / (n - 1); - - float baseRadius = 25f; - float rangeFactor = 2f; - float scaleFactor = 2f; - - static float SmoothStep(float edge0, float edge1, float x) - { - float t = Math.Clamp((x - edge0) / (edge1 - edge0), 0f, 1f); - return t * t * (3f - 2f * t); - } - - var centersY = new float[n]; - var radii = new float[n]; - var temperatures = new double[n]; - double R_gas = 287.0; - - for (int i = 0; i < n; i++) - { - double p = pipe.GetCellPressure(i); - double rho = pipe.GetCellDensity(i); - double T = p / Math.Max(rho * R_gas, 1e-12); - temperatures[i] = T; - - float deviation = (float)Math.Tanh((p - AmbientPressure) / AmbientPressure / rangeFactor); - radii[i] = baseRadius * (1f + deviation * scaleFactor); - if (radii[i] < 2f) radii[i] = 2f; - centersY[i] = topY + i * dy; - } - - int segmentsPerCell = 8; - int totalPoints = n + (n - 1) * segmentsPerCell; - Vertex[] stripVertices = new Vertex[totalPoints * 2]; - int idx = 0; - - for (int i = 0; i < n; i++) - { - float y = centersY[i]; - float r = radii[i]; - Color col = TemperatureColor(temperatures[i]); - - stripVertices[idx++] = new Vertex(new Vector2f(centerX - r, y), col); - stripVertices[idx++] = new Vertex(new Vector2f(centerX + r, y), col); - - if (i < n - 1) - { - for (int s = 1; s <= segmentsPerCell; s++) - { - float t = s / (float)segmentsPerCell; - float st = SmoothStep(0f, 1f, t); - float yi = centersY[i] + (centersY[i + 1] - centersY[i]) * t; - float ri = radii[i] + (radii[i + 1] - radii[i]) * st; - double Ti = temperatures[i] + (temperatures[i + 1] - temperatures[i]) * st; - Color coli = TemperatureColor(Ti); - - stripVertices[idx++] = new Vertex(new Vector2f(centerX - ri, yi), coli); - stripVertices[idx++] = new Vertex(new Vector2f(centerX + ri, yi), coli); - } - } - } - - var pipeMesh = new VertexArray(PrimitiveType.TriangleStrip, (uint)stripVertices.Length); - for (int i = 0; i < stripVertices.Length; i++) - pipeMesh[(uint)i] = stripVertices[i]; - target.Draw(pipeMesh); - } - } -} \ No newline at end of file diff --git a/Scenarios/Scenario.cs b/Scenarios/Scenario.cs index d519112..c8a2cde 100644 --- a/Scenarios/Scenario.cs +++ b/Scenarios/Scenario.cs @@ -1,72 +1,60 @@ -using System; -using SFML.Graphics; +using SFML.Graphics; using SFML.System; +using FluidSim.Core; using FluidSim.Components; namespace FluidSim.Tests { public abstract class Scenario { + protected const float AmbientPressure = 101325f; + protected const float AmbientTemperature = 300f; + public float Throttle { get; set; } + public abstract void Initialize(int sampleRate); public abstract float Process(); public abstract void Draw(RenderWindow target); - protected const double AmbientPressure = 101325.0; - protected const double AmbientTemperature = 300.0; - public double Throttle { get; set; } = 0.0; - - // ---------- Color from pressure (volumes) ---------- - protected Color PressureColor(double pressurePa) + protected Color PressureColor(float pressurePa) { - double bar = pressurePa / 1e5; // convert to bar for easier mapping + float bar = pressurePa / 1e5f; byte r, g, b; - - if (bar < 1.0) // vacuum → blue to green + if (bar < 1f) { - double factor = Math.Clamp(bar, 0.0, 1.0); - r = 0; - g = (byte)(255 * factor); - b = (byte)(255 * (1.0 - factor)); - } - else // above ambient → green to red - { - double factor = Math.Min((bar - 1.0) / 9.0, 1.0); // 1→10 bar maps to 0→1 - r = (byte)(255 * factor); - g = (byte)(255 * (1.0 - factor)); - b = 0; - } - return new Color(r, g, b); - } - - // ---------- Color from temperature (pipes) ---------- - protected Color TemperatureColor(double temperature) - { - double t = Math.Clamp(temperature, 0.0, 2000.0); - byte r, g, b; - if (t < AmbientTemperature) - { - double factor = t / AmbientTemperature; - r = 0; - g = (byte)(255 * factor); - b = (byte)(255 * (1.0 - factor)); + float f = Math.Clamp(bar, 0f, 1f); + r = 0; g = (byte)(255 * f); b = (byte)(255 * (1 - f)); } else { - double factor = (t - AmbientTemperature) / (2000.0 - AmbientTemperature); - r = (byte)(255 * factor); - g = (byte)(255 * (1.0 - factor)); - b = 0; + float f = Math.Min((bar - 1f) / 9f, 1f); + r = (byte)(255 * f); g = (byte)(255 * (1 - f)); b = 0; + } + return new Color(r, g, b); + } + + protected Color TemperatureColor(float t) + { + t = Math.Clamp(t, 0f, 2000f); + byte r, g, b; + if (t < AmbientTemperature) + { + float f = t / AmbientTemperature; + r = 0; g = (byte)(255 * f); b = (byte)(255 * (1 - f)); + } + else + { + float f = (t - AmbientTemperature) / (2000f - AmbientTemperature); + r = (byte)(255 * f); g = (byte)(255 * (1 - f)); b = 0; } return new Color(r, g, b); } - // ---------- Draw a generic volume (e.g. plenum) ---------- protected void DrawVolume(RenderWindow target, Volume0D volume, float centerX, float topY, float width, float height) { var rect = new RectangleShape(new Vector2f(width, height)) { - FillColor = PressureColor(volume.Pressure), // ← pressure‑based + FillColor = PressureColor(volume.Pressure), Position = new Vector2f(centerX - width / 2f, topY) }; target.Draw(rect); @@ -75,122 +63,99 @@ namespace FluidSim.Tests FillColor = Color.Transparent, OutlineColor = Color.White, OutlineThickness = 1f, - Position = new Vector2f(centerX - width / 2f, topY) + Position = rect.Position }; target.Draw(border); } - // ---------- Draw an engine cylinder ---------- protected void DrawCylinder(RenderWindow target, Cylinder cylinder, float centerX, float topY, float width, float maxHeight) { - double fraction = cylinder.PistonFraction; - float currentHeight = (float)(maxHeight * fraction); - - // Walls - var wall = new RectangleShape(new Vector2f(width, maxHeight)); - wall.FillColor = new Color(60, 60, 60); - wall.Position = new Vector2f(centerX - width / 2f, topY); + float fraction = cylinder.PistonFraction; + float currentHeight = maxHeight * fraction; + var wall = new RectangleShape(new Vector2f(width, maxHeight)) + { + FillColor = new Color(60, 60, 60), + Position = new Vector2f(centerX - width / 2f, topY) + }; target.Draw(wall); - - // Gas – colored by pressure now - float gasTop = topY; - var gasRect = new RectangleShape(new Vector2f(width, currentHeight)); - gasRect.FillColor = PressureColor(cylinder.Pressure); // ← pressure‑based - gasRect.Position = new Vector2f(centerX - width / 2f, gasTop); - target.Draw(gasRect); - - // Piston line - var pistonLine = new RectangleShape(new Vector2f(width, 4f)); - pistonLine.FillColor = Color.White; - pistonLine.Position = new Vector2f(centerX - width / 2f, topY + currentHeight); - target.Draw(pistonLine); - - // Valve indicators + var gas = new RectangleShape(new Vector2f(width, currentHeight)) + { + FillColor = PressureColor(cylinder.Pressure), + Position = new Vector2f(centerX - width / 2f, topY) + }; + target.Draw(gas); + var piston = new RectangleShape(new Vector2f(width, 4f)) + { + FillColor = Color.White, + Position = new Vector2f(centerX - width / 2f, topY + currentHeight) + }; + target.Draw(piston); float valveW = 6f, valveH = 10f, valveY = topY + 4f; - var intakeValve = new RectangleShape(new Vector2f(valveW, valveH)); - intakeValve.FillColor = cylinder.IntakeValveArea > 0 ? Color.Green : Color.Red; - intakeValve.Position = new Vector2f(centerX - width / 2f - valveW - 2f, valveY); - target.Draw(intakeValve); - - var exhaustValve = new RectangleShape(new Vector2f(valveW, valveH)); - exhaustValve.FillColor = cylinder.ExhaustValveArea > 0 ? Color.Green : Color.Red; - exhaustValve.Position = new Vector2f(centerX + width / 2f + 2f, valveY); - target.Draw(exhaustValve); + var iv = new RectangleShape(new Vector2f(valveW, valveH)) + { + FillColor = cylinder.IntakeValveArea > 0f ? Color.Green : Color.Red, + Position = new Vector2f(centerX - width / 2f - valveW - 2f, valveY) + }; + target.Draw(iv); + var ev = new RectangleShape(new Vector2f(valveW, valveH)) + { + FillColor = cylinder.ExhaustValveArea > 0f ? Color.Green : Color.Red, + Position = new Vector2f(centerX + width / 2f + 2f, valveY) + }; + target.Draw(ev); } - // ---------- Draw a pipe (unchanged) ---------- - protected void DrawPipe(RenderWindow target, Pipe1D pipe, float pipeCenterY, float pipeStartX, float pipeEndX) + protected void DrawPipe(RenderWindow target, PipeSystem pipeSystem, int pipeIndex, + float pipeCenterY, float pipeStartX, float pipeEndX) { - int n = pipe.CellCount; + int start = pipeSystem.GetPipeStart(pipeIndex); + int end = pipeSystem.GetPipeEnd(pipeIndex); + int n = end - start; if (n < 2) return; - float pipeLengthPx = pipeEndX - pipeStartX; - float dx = pipeLengthPx / (n - 1); - + float pipeLen = pipeEndX - pipeStartX; + float dx = pipeLen / (n - 1); float baseRadius = 25f; - float rangeFactor = 2f; - float scaleFactor = 2f; - - static float SmoothStep(float edge0, float edge1, float x) - { - float t = Math.Clamp((x - edge0) / (edge1 - edge0), 0f, 1f); - return t * t * (3f - 2f * t); - } var centers = new float[n]; var radii = new float[n]; - var temperatures = new double[n]; - double R_gas = 287.0; - + var temps = new float[n]; for (int i = 0; i < n; i++) { - double p = pipe.GetCellPressure(i); - double rho = pipe.GetCellDensity(i); - double T = p / Math.Max(rho * R_gas, 1e-12); - temperatures[i] = T; - - float deviation = (float)Math.Tanh((p - AmbientPressure) / AmbientPressure / rangeFactor); - radii[i] = baseRadius * (1f + deviation * scaleFactor); + int cell = start + i; + float p = pipeSystem.GetCellPressure(cell); + float rho = pipeSystem.GetCellDensity(cell); + temps[i] = p / MathF.Max(rho * 287f, 1e-12f); + float dev = MathF.Tanh((p - AmbientPressure) / AmbientPressure * 0.5f); + radii[i] = baseRadius * (1f + dev * 2f); if (radii[i] < 2f) radii[i] = 2f; centers[i] = pipeStartX + i * dx; } - int segmentsPerCell = 8; - int totalPoints = n + (n - 1) * segmentsPerCell; - Vertex[] stripVertices = new Vertex[totalPoints * 2]; - int idx = 0; - + int segments = 8; + var va = new VertexArray(PrimitiveType.TriangleStrip); for (int i = 0; i < n; i++) { - float x = centers[i]; - float r = radii[i]; - Color col = TemperatureColor(temperatures[i]); // pipes still use temperature - - stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY - r), col); - stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY + r), col); - + float x = centers[i], r = radii[i]; + Color col = TemperatureColor(temps[i]); + va.Append(new Vertex(new Vector2f(x, pipeCenterY - r), col)); + va.Append(new Vertex(new Vector2f(x, pipeCenterY + r), col)); if (i < n - 1) { - for (int s = 1; s <= segmentsPerCell; s++) + for (int s = 1; s <= segments; s++) { - float t = s / (float)segmentsPerCell; - float st = SmoothStep(0f, 1f, t); + float t = s / (float)segments; float xi = centers[i] + (centers[i + 1] - centers[i]) * t; - float ri = radii[i] + (radii[i + 1] - radii[i]) * st; - double Ti = temperatures[i] + (temperatures[i + 1] - temperatures[i]) * st; - Color coli = TemperatureColor(Ti); - - stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY - ri), coli); - stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY + ri), coli); + float ri = radii[i] + (radii[i + 1] - radii[i]) * t; + float Ti = temps[i] + (temps[i + 1] - temps[i]) * t; + Color colS = TemperatureColor(Ti); + va.Append(new Vertex(new Vector2f(xi, pipeCenterY - ri), colS)); + va.Append(new Vertex(new Vector2f(xi, pipeCenterY + ri), colS)); } } } - - var pipeMesh = new VertexArray(PrimitiveType.TriangleStrip, (uint)stripVertices.Length); - for (int i = 0; i < stripVertices.Length; i++) - pipeMesh[(uint)i] = stripVertices[i]; - target.Draw(pipeMesh); + target.Draw(va); } } } \ No newline at end of file diff --git a/Scenarios/SingleCylScenario.cs b/Scenarios/SingleCylScenario.cs new file mode 100644 index 0000000..f149fac --- /dev/null +++ b/Scenarios/SingleCylScenario.cs @@ -0,0 +1,220 @@ +using FluidSim.Components; +using FluidSim.Core; +using FluidSim.Interfaces; +using SFML.Graphics; +using SFML.System; +using System; + +namespace FluidSim.Tests +{ + public class SingleCylScenario : Scenario + { + private Crankshaft crankshaft; + private Cylinder cylinder; + + private PipeSystem pipeSystem; + private BoundarySystem boundaries; + private Solver solver; + + private Volume0D intakePlenum; + private Port plenumInlet, plenumOutlet; + private Volume0D exhaustCollector; + private Port colIn, colOut; + + private int throttleAreaIdx, plenumRunnerAreaIdx, intakeValveIdx, exhaustValveIdx; + private float[] orificeAreas; + private int intakeOpenIdx, exhaustOpenIdx; + + private SoundProcessor exhaustSound, intakeSound; + private OutdoorExhaustReverb reverb; + + private double dt; + private int stepCount; + public float MaxThrottleArea = 1e-4f; // 1 cm² + + // pipe area for open end calculations + private float pipeArea; + + public override void Initialize(int sampleRate) + { + dt = 1.0 / sampleRate; + + // ---- Crankshaft ---- + crankshaft = new Crankshaft(600); + crankshaft.Inertia = 0.2f; + crankshaft.FrictionConstant = 2f; + crankshaft.FrictionViscous = 0.04f; + + // ---- Cylinder ---- + float bore = 0.056f, stroke = 0.057f, conRod = 0.110f, compRatio = 9.2f; + float ivo = 350f, ivc = 580f, evo = 120f, evc = 370f; + cylinder = new Cylinder(bore, stroke, conRod, compRatio, + ivo, ivc, evo, evc, crankshaft) + { + IntakeValveDiameter = 0.03f, + IntakeValveLift = 0.005f, + ExhaustValveDiameter = 0.028f, + ExhaustValveLift = 0.005f + }; + + // ---- Pipe system ---- + int totalCells = 10 + 10 + 50; + int[] pipeStart = { 0, 10, 20 }; + int[] pipeEnd = { 10, 20, 70 }; + float[] area = new float[totalCells]; + float[] dx = new float[totalCells]; + float pipeDiameter = 0.02f; // 2 cm + pipeArea = MathF.PI * 0.25f * pipeDiameter * pipeDiameter; + float areaVal = pipeArea; + float intakeLenBefore = 0.2f, intakeLenRunner = 0.2f, exhaustLen = 0.5f; + for (int i = 0; i < totalCells; i++) + { + area[i] = areaVal; + if (i < 10) dx[i] = intakeLenBefore / 10f; + else if (i < 20) dx[i] = intakeLenRunner / 10f; + else dx[i] = exhaustLen / 50f; + } + + pipeSystem = new PipeSystem(totalCells, pipeStart, pipeEnd, area, dx, + 1.225f, 0f, 101325f); + pipeSystem.DampingMultiplier = 0.5f; + pipeSystem.EnergyRelaxationRate = 0f; + pipeSystem.AmbientPressure = 101325f; + + // ---- Volumes ---- + intakePlenum = new Volume0D(5e-6f, 101325f, 300f); // 5 mL + plenumInlet = intakePlenum.CreatePort(); + plenumOutlet = intakePlenum.CreatePort(); + exhaustCollector = new Volume0D(10e-6f, 101325f, 800f); // 10 mL (unused but present) + colIn = exhaustCollector.CreatePort(); + colOut = exhaustCollector.CreatePort(); + + // ---- Boundary system ---- + boundaries = new BoundarySystem(pipeSystem, maxOrifices: 4, maxOpenEnds: 2); + + throttleAreaIdx = 0; + plenumRunnerAreaIdx = 1; + intakeValveIdx = 2; + exhaustValveIdx = 3; + + // Intake open end (pipe0 left) + boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: true, 101325f, pipeArea); + intakeOpenIdx = 0; + + // Throttle orifice (plenum inlet to pipe0 right) + boundaries.AddOrifice(plenumInlet, pipeIndex: 0, isLeftEnd: false, throttleAreaIdx, 0.2f); + + // Plenum to runner (plenum outlet to pipe1 left) + boundaries.AddOrifice(plenumOutlet, pipeIndex: 1, isLeftEnd: true, plenumRunnerAreaIdx, 1f); + + // Intake valve (cylinder intake to pipe1 right) + boundaries.AddOrifice(cylinder.IntakePort, pipeIndex: 1, isLeftEnd: false, intakeValveIdx, 1f); + + // Exhaust valve (cylinder exhaust to pipe2 left) + boundaries.AddOrifice(cylinder.ExhaustPort, pipeIndex: 2, isLeftEnd: true, exhaustValveIdx, 1f); + + // Exhaust open end (pipe2 right) + boundaries.AddOpenEnd(pipeIndex: 2, isLeftEnd: false, 101325f, pipeArea); + exhaustOpenIdx = 1; + + orificeAreas = new float[4]; + orificeAreas[plenumRunnerAreaIdx] = areaVal; // fixed plenum->runner area + + // ---- Solver ---- + solver = new Solver { SubStepCount = 4, EnableProfiling = false }; + solver.SetTimeStep(dt); + solver.SetPipeSystem(pipeSystem); + solver.SetBoundarySystem(boundaries); + solver.AddComponent(cylinder); + solver.AddComponent(intakePlenum); + solver.AddComponent(exhaustCollector); + + // ---- Sound ---- + exhaustSound = new SoundProcessor(sampleRate, 1f) { Gain = 1f }; + intakeSound = new SoundProcessor(sampleRate, 1f) { Gain = 1f }; + reverb = new OutdoorExhaustReverb(sampleRate); + + stepCount = 0; + Console.WriteLine("TestScenario ready."); + } + + public override float Process() + { + crankshaft.Step((float)dt); + cylinder.PreStep((float)dt); + + // Update variable orifice areas + float throttledArea = MaxThrottleArea * Math.Clamp(Throttle, 0.0001f, 1f); + orificeAreas[throttleAreaIdx] = throttledArea; + orificeAreas[intakeValveIdx] = cylinder.IntakeValveArea; + orificeAreas[exhaustValveIdx] = cylinder.ExhaustValveArea; + boundaries.SetOrificeAreas(orificeAreas); + + solver.Step(); + stepCount++; + + // Retrieve open‑end mass flows for sound synthesis + float exhaustFlow = boundaries.GetOpenEndMassFlow(exhaustOpenIdx); + float intakeFlow = boundaries.GetOpenEndMassFlow(intakeOpenIdx); + + float exhaustDry = exhaustSound.Process(exhaustFlow); + float intakeDry = intakeSound.Process(intakeFlow); + + if (stepCount % 1000 == 0) + { + float rpm = crankshaft.AngularVelocity * 60f / (2f * MathF.PI); + Console.WriteLine($"Step {stepCount}, RPM={rpm:F0}, CylP={cylinder.Pressure / 1e5f:F2} bar"); + Console.WriteLine($"intake flow: {intakeFlow:F12}, exhaust flow: {exhaustFlow:F16}"); + } + + return reverb.Process(intakeDry); + } + + public override void Draw(RenderWindow target) + { + float winW = target.GetView().Size.X; + float winH = target.GetView().Size.Y; + + float intakeY = winH / 2f - 40f; + float exhaustY = winH / 2f + 80f; + float openEndX = 40f; + + // Intake pipe before throttle (pipe 0) + float pipe1StartX = openEndX; + float pipe1EndX = pipe1StartX + 120f; + DrawPipe(target, pipeSystem, 0, intakeY, pipe1StartX, pipe1EndX); + + // Throttle symbol + float throttleX = pipe1EndX + 5f; + var throttleRect = new RectangleShape(new Vector2f(8f, 30f)) + { + FillColor = Color.Yellow, + Position = new Vector2f(throttleX, intakeY - 15f) + }; + target.Draw(throttleRect); + + // Plenum + float plenW = 60f, plenH = 80f; + float plenLeftX = throttleX + 10f; + float plenCenterX = plenLeftX + plenW / 2f; + float plenTopY = intakeY - plenH / 2f; + DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH); + + // Runner pipe (pipe 1) + float runnerStartX = plenLeftX + plenW + 5f; + float runnerEndX = runnerStartX + 100f; + DrawPipe(target, pipeSystem, 1, intakeY, runnerStartX, runnerEndX); + + // Cylinder + float cylCX = runnerEndX + 50f; + float cylTopY = intakeY - 120f; + float cylW = 80f, cylMaxH = 240f; + DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH); + + // Exhaust pipe (pipe 2) + float exhStartX = cylCX + cylW / 2f + 20f; + float exhEndX = winW - 60f; + DrawPipe(target, pipeSystem, 2, exhaustY, exhStartX, exhEndX); + } + } +} \ No newline at end of file diff --git a/Scenarios/TestScenario.cs b/Scenarios/TestScenario.cs index 6996863..9e00f9f 100644 --- a/Scenarios/TestScenario.cs +++ b/Scenarios/TestScenario.cs @@ -1,176 +1,91 @@ -using System; +using System; using SFML.Graphics; using SFML.System; -using FluidSim.Components; using FluidSim.Core; -using FluidSim.Utils; namespace FluidSim.Tests { public class TestScenario : Scenario { - // Engine - private Cylinder cylinder; - private Crankshaft crankshaft; - - // Intake side - private Pipe1D intakePipeBeforeThrottle; - private Volume0D intakePlenum; // 5 mL - private Pipe1D intakeRunner; - - // Exhaust side - private Pipe1D exhaustPipe; - - // Links - private OpenEndLink intakeOpenEnd; - private OrificeLink throttleOrifice; - private OrificeLink plenumToRunner; - private OrificeLink intakeValve; - private OrificeLink exhaustValve; - private OpenEndLink exhaustOpenEnd; - + private PipeSystem pipeSystem; + private BoundarySystem boundaries; private Solver solver; - private SoundProcessor exhaustSoundProcessor; - private SoundProcessor intakeSoundProcessor; - private OutdoorExhaustReverb reverb; + + private int[] pipeStart = { 0 }; + private int[] pipeEnd; + private double dt; private int stepCount; - // ---------- Throttle control ---------- - public double MaxThrottleArea { get; set; } = 1 * Units.cm2; // 2 cm² + // Sound output: use pressure at open end + private SoundProcessor openEndSound; + private int openEndIdx = 0; // index of the open end in BoundarySystem (we added only one) public override void Initialize(int sampleRate) { dt = 1.0 / sampleRate; - solver = new Solver(); + const int cellCount = 200; + float length = 2f; + float dia = 0.02f; + float area = MathF.PI * 0.25f * dia * dia; + + float[] areas = new float[cellCount]; + float[] dxs = new float[cellCount]; + float dx = length / cellCount; + for (int i = 0; i < cellCount; i++) + { + areas[i] = area; + dxs[i] = dx; + } + + pipeEnd = new[] { cellCount }; + + float rho0 = 101325f / (287f * 300f); + pipeSystem = new PipeSystem(cellCount, pipeStart, pipeEnd, areas, dxs, + rho0, 0f, 101325f); + pipeSystem.DampingMultiplier = 0f; + pipeSystem.EnergyRelaxationRate = 0f; + pipeSystem.AmbientPressure = 101325f; + + // Pressure bubble near right end + float pBubble = 10f * 101325f; + float TBubble = 2000f; + float rhoBubble = pBubble / (287f * TBubble); + for (int i = 0; i <= 10; i++) + pipeSystem.SetCellState(i, rhoBubble, 0f, pBubble); + + // Boundaries: left closed, right open + boundaries = new BoundarySystem(pipeSystem, maxOrifices: 1, maxOpenEnds: 1); + boundaries.AddOrifice(null, pipeIndex: 0, isLeftEnd: true, areaIndex: 0, 1f); + boundaries.AddOpenEnd(pipeIndex: 0, isLeftEnd: false, 101325f, area); + float[] orificeAreas = new float[1] { 0f }; + boundaries.SetOrificeAreas(orificeAreas); + + solver = new Solver { SubStepCount = 3}; solver.SetTimeStep(dt); - solver.CflTarget = 0.9; + solver.SetPipeSystem(pipeSystem); + solver.SetBoundarySystem(boundaries); - // ---- Crankshaft (external, passed to cylinder) ---- - crankshaft = new Crankshaft(600); - crankshaft.Inertia = 0.2; - crankshaft.FrictionConstant = 2; - crankshaft.FrictionViscous = 0.04; + solver.EnableProfiling = true; + pipeSystem.EnableProfiling = true; - // ---- Cylinder ---- - double bore = 0.056, stroke = 0.057, conRod = 0.110, compRatio = 9.2; - double ivo = 350.0, ivc = 580.0, evo = 120.0, evc = 370.0; - cylinder = new Cylinder(bore, stroke, conRod, compRatio, ivo, ivc, evo, evc, crankshaft) - { - IntakeValveDiameter = 30 * Units.mm, // 30 mm - IntakeValveLift = 5 * Units.mm, // 5 mm - ExhaustValveDiameter = 28 * Units.mm, // 28 mm - ExhaustValveLift = 5 * Units.mm // 5 mm - }; - solver.AddComponent(cylinder); - - double pipeDiameter = 2 * Units.cm; - double pipeArea = Units.AreaFromDiameter(pipeDiameter); - - exhaustSoundProcessor = new SoundProcessor(sampleRate, 1, pipeDiameter) { Gain = 0.1f }; - intakeSoundProcessor = new SoundProcessor(sampleRate, 1, pipeDiameter) { Gain = 0.1f }; - reverb = new OutdoorExhaustReverb(sampleRate); - - // ---- Pipes ---- - intakePipeBeforeThrottle = new Pipe1D(0.2, pipeArea, 10); - intakeRunner = new Pipe1D(0.2, pipeArea, 10); - exhaustPipe = new Pipe1D(0.5, pipeArea, 50); - solver.AddComponent(intakePipeBeforeThrottle); - solver.AddComponent(intakeRunner); - solver.AddComponent(exhaustPipe); - - intakePlenum = new Volume0D(5 * Units.mL, 101325.0, 300.0); - var plenumInlet = intakePlenum.CreatePort(); - var plenumOutlet = intakePlenum.CreatePort(); - solver.AddComponent(intakePlenum); - - // ---- Intake open end ---- - intakeOpenEnd = new OpenEndLink(intakePipeBeforeThrottle, isLeftEnd: true) - { - AmbientPressure = 101325.0, - Gamma = 1.4 - }; - solver.AddOpenEndLink(intakeOpenEnd); - - // ---- Throttle orifice (variable area) ---- - throttleOrifice = new OrificeLink(plenumInlet, intakePipeBeforeThrottle, isPipeLeftEnd: false, - areaProvider: () => MaxThrottleArea * Math.Clamp(Throttle, 0.0001, 1)) - { - DischargeCoefficient = 0.2, - UseInertance = false - }; - solver.AddOrificeLink(throttleOrifice); - - // ---- Plenum to runner (fixed area) ---- - plenumToRunner = new OrificeLink(plenumOutlet, intakeRunner, isPipeLeftEnd: true, - areaProvider: () => pipeArea) - { - DischargeCoefficient = 1.0, - UseInertance = false - }; - solver.AddOrificeLink(plenumToRunner); - - // ---- Intake valve ---- - intakeValve = new OrificeLink(cylinder.IntakePort, intakeRunner, isPipeLeftEnd: false, - areaProvider: () => cylinder.IntakeValveArea) - { - DischargeCoefficient = 1.0, - UseInertance = false - }; - solver.AddOrificeLink(intakeValve); - - // ---- Exhaust valve ---- - exhaustValve = new OrificeLink(cylinder.ExhaustPort, exhaustPipe, isPipeLeftEnd: true, - areaProvider: () => cylinder.ExhaustValveArea) - { - DischargeCoefficient = 1.0, - UseInertance = false - }; - solver.AddOrificeLink(exhaustValve); - - // ---- Exhaust open end ---- - exhaustOpenEnd = new OpenEndLink(exhaustPipe, isLeftEnd: false) - { - AmbientPressure = 101325.0, - Gamma = 1.4 - }; - solver.AddOpenEndLink(exhaustOpenEnd); + // Simple sound processor: convert mass flow rate to audio + openEndSound = new SoundProcessor(sampleRate, 1f) { Gain = 2f }; + Console.WriteLine("Pulse test ready."); stepCount = 0; - Console.WriteLine("4‑Stroke engine test (plenum + two pipes)"); - Console.WriteLine($"Bore {bore * 1000:F0}mm, Stroke {stroke * 1000:F0}mm, CR {compRatio}"); - Console.WriteLine($"IVO {ivo}°, IVC {ivc}°, EVO {evo}°, EVC {evc}° (no overlap)"); } public override float Process() { - cylinder.Crankshaft.Step(dt); - cylinder.PreStep(dt); solver.Step(); stepCount++; - if (stepCount % 10000 == 0) - { - double crankDeg = cylinder.Crankshaft.CrankAngle * 180.0 / Math.PI % 720.0; - double cylP = cylinder.Pressure / 1e5; - double cylT = cylinder.Temperature; - double cylMass = cylinder.Mass * 1e6; - double mdotI = intakeValve.LastMassFlowRate; - double mdotE = exhaustValve.LastMassFlowRate; - double pipeR = exhaustPipe.GetCellPressure(exhaustPipe.CellCount - 1) / 1e5; - double plenumP = intakePlenum.Pressure / 1e5; - double actualArea = MaxThrottleArea * Throttle; + float flow = boundaries.GetOpenEndMassFlow(openEndIdx); + float sample = openEndSound.Process(flow); - Console.WriteLine($"Step {stepCount}: Angle={crankDeg:F1}°, " + - $"CylP={cylP:F2} bar, T={cylT:F0} K, mass={cylMass:F1} mg, " + - $"mdotI={mdotI:E4} kg/s, mdotE={mdotE:E4} kg/s, PipeR={pipeR:F2} bar"); - Console.WriteLine($"Throttle = {Throttle * 100:F0}% area = {actualArea * 1e6:F2} mm², Plenum P = {plenumP:F3} bar"); - } - - float exhaustDry = exhaustSoundProcessor.Process(exhaustOpenEnd); - float intakeDry = intakeSoundProcessor.Process(intakeOpenEnd); - return reverb.Process(exhaustDry + intakeDry); + return sample; } public override void Draw(RenderWindow target) @@ -178,56 +93,10 @@ namespace FluidSim.Tests float winW = target.GetView().Size.X; float winH = target.GetView().Size.Y; - float intakeY = winH / 2f - 40f; - float exhaustY = winH / 2f + 80f; - - // Open end marker - float openEndX = 40f; - var openEndMark = new CircleShape(5f) { FillColor = Color.Cyan }; - openEndMark.Position = new Vector2f(openEndX - 5f, intakeY - 5f); - target.Draw(openEndMark); - - // First intake pipe - float pipe1StartX = openEndX; - float pipe1EndX = pipe1StartX + 120f; - DrawPipe(target, intakePipeBeforeThrottle, intakeY, pipe1StartX, pipe1EndX); - - // Throttle symbol - float throttleX = pipe1EndX + 5f; - var throttleRect = new RectangleShape(new Vector2f(8f, 30f)) - { - FillColor = Color.Yellow, - Position = new Vector2f(throttleX, intakeY - 15f) - }; - target.Draw(throttleRect); - - // Plenum - float plenW = 60f, plenH = 80f; - float plenLeftX = throttleX + 10f; - float plenCenterX = plenLeftX + plenW / 2f; - float plenTopY = intakeY - plenH / 2f; - DrawVolume(target, intakePlenum, plenCenterX, plenTopY, plenW, plenH); - - // Runner pipe - float runnerStartX = plenLeftX + plenW + 5f; - float runnerEndX = runnerStartX + 100f; - DrawPipe(target, intakeRunner, intakeY, runnerStartX, runnerEndX); - - // Cylinder - float cylCX = runnerEndX + 50f; - float cylTopY = intakeY - 120f; - float cylW = 80f, cylMaxH = 240f; - DrawCylinder(target, cylinder, cylCX, cylTopY, cylW, cylMaxH); - - // Exhaust pipe - float exhStartX = cylCX + cylW / 2f + 20f; - float exhEndX = winW - 60f; - DrawPipe(target, exhaustPipe, exhaustY, exhStartX, exhEndX); - - // Exhaust open end marker - var exhOpenEndMark = new CircleShape(5f) { FillColor = Color.Magenta }; - exhOpenEndMark.Position = new Vector2f(exhEndX - 5f, exhaustY - 5f); - target.Draw(exhOpenEndMark); + float startX = 50f; + float endX = winW - 50f; + float y = winH / 2f; + DrawPipe(target, pipeSystem, 0, y, startX, endX); } } } \ No newline at end of file