diff --git a/Components/Atmosphere.cs b/Components/Atmosphere.cs new file mode 100644 index 0000000..1558240 --- /dev/null +++ b/Components/Atmosphere.cs @@ -0,0 +1,43 @@ +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 double Density => Pressure / (GasConstant * Temperature); + public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density; + + public Port Port { get; } + + public Atmosphere() + { + Port = new Port { Owner = this }; + UpdatePort(); + } + + public IReadOnlyList Ports => new[] { Port }; + + public void UpdateState(double dt) + { + // Atmosphere is static – just ensure the port reflects current values + UpdatePort(); + } + + private void UpdatePort() + { + Port.Pressure = Pressure; + Port.Density = Density; + Port.Temperature = Temperature; + Port.SpecificEnthalpy = SpecificEnthalpy; + // MassFlowRate is set by the solver connector + } + } +} \ No newline at end of file diff --git a/Components/EngineCylinder.cs b/Components/EngineCylinder.cs deleted file mode 100644 index 348aecd..0000000 --- a/Components/EngineCylinder.cs +++ /dev/null @@ -1,249 +0,0 @@ -using System; -using FluidSim.Components; - -namespace FluidSim.Core -{ - public class EngineCylinder - { - public Volume0D Cylinder { get; private set; } - private Crankshaft crankshaft; - - private double bore, stroke, conRodLength, compressionRatio; - private double pistonArea; - - public double V_disp { get; private set; } - public double V_clear { get; private set; } - public bool ignition = false; - - // ---- Exhaust valve ---- - private double exhMaxOrificeArea; - private double exhValveOpenStart = 130.0 * Math.PI / 180.0; - private double exhValveOpenEnd = 390.0 * Math.PI / 180.0; - private double exhValveRampWidth = 30.0 * Math.PI / 180.0; - public double ExhaustOrificeArea => ExhaustValveLift() * exhMaxOrificeArea; - public double ExhaustValveLiftCurrent => ExhaustValveLift(); - - // ---- Intake valve ---- - private double intMaxOrificeArea; - private double intValveOpenStart = 340.0 * Math.PI / 180.0; - private double intValveOpenEnd = 600.0 * Math.PI / 180.0; - private double intValveRampWidth = 30.0 * Math.PI / 180.0; - public double IntakeOrificeArea => IntakeValveLift() * intMaxOrificeArea; - public double IntakeValveLiftCurrent => IntakeValveLift(); - - // ---- Combustion ---- - public double TargetPeakPressure { get; set; } = 50.0 * 101325.0; - private const double PeakTemperature = 2500.0; - private bool burnInProgress = false; - private double burnStartAngle; // cycle angle (0–4π) - private double burnDuration = 40.0 * Math.PI / 180.0; - private double targetBurnEnergy; - private double preIgnitionMass, preIgnitionInternalEnergy; - - private Random rand = new Random(); - public double MisfireProbability { get; set; } = 0.02; - private bool misfireCurrent = false; - - public int CombustionCount { get; private set; } - public int MisfireCount { get; private set; } - - // Cycle‑aware angle (0 – 4π) - public double CycleAngle => crankshaft.CrankAngle % (4.0 * Math.PI); - private double prevCycleAngle; - - // Piston position fraction (0 = TDC, 1 = BDC) - public double PistonPositionFraction - { - get - { - double currentVol = Cylinder.Volume; - if (currentVol <= V_clear) return 0.0; - if (currentVol >= V_clear + V_disp) return 1.0; - return (currentVol - V_clear) / V_disp; - } - } - - public EngineCylinder(Crankshaft crankshaft, - double bore, double stroke, double compressionRatio, - double exhPipeArea, double intPipeArea, int sampleRate) - { - this.crankshaft = crankshaft; - this.bore = bore; - this.stroke = stroke; - conRodLength = 2.0 * stroke; - this.compressionRatio = compressionRatio; - exhMaxOrificeArea = exhPipeArea * 0.5; - intMaxOrificeArea = intPipeArea * 0.5; - pistonArea = Math.PI / 4.0 * bore * bore; - - V_disp = pistonArea * stroke; - V_clear = V_disp / (compressionRatio - 1.0); - - // Start at BDC with fresh ambient charge - double V_bdc = V_clear + V_disp; - double p_amb = 101325.0; - double T_amb = 300.0; - double rho0 = p_amb / (287.0 * T_amb); - double mass0 = rho0 * V_bdc; - double energy0 = p_amb * V_bdc / (1.4 - 1.0); - - Cylinder = new Volume0D(V_bdc, p_amb, T_amb, sampleRate) - { - Gamma = 1.4, - GasConstant = 287.0 - }; - Cylinder.Volume = V_bdc; - Cylinder.Mass = mass0; - Cylinder.InternalEnergy = energy0; - - prevCycleAngle = CycleAngle; - - preIgnitionMass = Cylinder.Mass; - preIgnitionInternalEnergy = Cylinder.InternalEnergy; - } - - // ---- Piston kinematics ---- - private (double volume, double dvdt) PistonKinematics(double cycleAngle) - { - double theta = cycleAngle % (2.0 * Math.PI); - double R = stroke / 2.0; - double cosT = Math.Cos(theta); - double sinT = Math.Sin(theta); - double L = conRodLength; - - double s = R * (1 - cosT) + L - Math.Sqrt(L * L - R * R * sinT * sinT); - double V = V_clear + pistonArea * s; - - double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT); - double dVdθ = pistonArea * (R * sinT + (R * R * sinT * cosT) / sqrtTerm); - double dvdt = dVdθ * crankshaft.AngularVelocity; - return (V, dvdt); - } - - // ---- Valve lifts (cycle‑aware) ---- - private double ExhaustValveLift() - { - double a = CycleAngle; - if (a < exhValveOpenStart || a > exhValveOpenEnd) return 0.0; - double duration = exhValveOpenEnd - exhValveOpenStart; - double ramp = exhValveRampWidth; - double t = (a - exhValveOpenStart) / duration; - double rampFrac = ramp / duration; - if (t < rampFrac) return t / rampFrac; - if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac; - return 1.0; - } - - private double IntakeValveLift() - { - double a = CycleAngle; - if (a < intValveOpenStart || a > intValveOpenEnd) return 0.0; - double duration = intValveOpenEnd - intValveOpenStart; - double ramp = intValveRampWidth; - double t = (a - intValveOpenStart) / duration; - double rampFrac = ramp / duration; - if (t < rampFrac) return t / rampFrac; - if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac; - return 1.0; - } - - // ---- Wiebe burn fraction ---- - private double WiebeFraction(double angleFromIgnition) - { - if (angleFromIgnition >= burnDuration) return 1.0; - double a = 5.0, m = 2.0; - double x = angleFromIgnition / burnDuration; - return 1.0 - Math.Exp(-a * Math.Pow(x, m + 1)); - } - - // ---- Torque from pressure ---- - private double ComputeTorque() - { - double p = Cylinder.Pressure; - double ambient = 101325.0; - double force = (p - ambient) * pistonArea; - if (force <= 0) return 0.0; - - double theta = crankshaft.CrankAngle % (2.0 * Math.PI); - double R = stroke / 2.0; - double L = conRodLength; - double sinT = Math.Sin(theta); - double cosT = Math.Cos(theta); - - double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT); - double lever = R * sinT * (1.0 + (R * cosT) / sqrtTerm); - return force * lever; - } - - // ---- TDC detection (power stroke, at angle 0 mod 4π) ---- - private bool DetectTDCPowerStroke() - { - double current = CycleAngle; - double previous = prevCycleAngle; - prevCycleAngle = current; - return (previous > 3.8 * Math.PI && current < 0.2 * Math.PI); - } - - public void Step(double dt) - { - bool crossingTDC = DetectTDCPowerStroke(); - - if (crossingTDC) - { - misfireCurrent = rand.NextDouble() < MisfireProbability; - - // *** Always capture the state at TDC, whether we burn or not *** - preIgnitionMass = Cylinder.Mass; - preIgnitionInternalEnergy = Cylinder.InternalEnergy; - - if (misfireCurrent) - { - MisfireCount++; - } - else if (ignition) - { - double V = Cylinder.Volume; - targetBurnEnergy = TargetPeakPressure * V / (Cylinder.Gamma - 1.0); - if (double.IsNaN(targetBurnEnergy)) - targetBurnEnergy = 101325.0 * V / (Cylinder.Gamma - 1.0); - burnInProgress = true; - burnStartAngle = CycleAngle; - CombustionCount++; - } - } - - if (burnInProgress) - { - double angleFromIgnition = CycleAngle - burnStartAngle; - if (angleFromIgnition < 0) angleFromIgnition += 4.0 * Math.PI; - - if (angleFromIgnition >= burnDuration) - { - Cylinder.InternalEnergy = targetBurnEnergy; - burnInProgress = false; - } - else - { - double fraction = WiebeFraction(angleFromIgnition); - Cylinder.InternalEnergy = preIgnitionInternalEnergy * (1.0 - fraction) - + targetBurnEnergy * fraction; - Cylinder.Mass = preIgnitionMass; - } - } - - var (vol, dvdt) = PistonKinematics(CycleAngle); - Cylinder.Volume = vol; - Cylinder.Dvdt = dvdt; - - if (double.IsNaN(Cylinder.Pressure) || double.IsNaN(Cylinder.Temperature) || Cylinder.Mass < 1e-9) - { - double V = Math.Max(vol, V_clear); - Cylinder.Mass = 1.225 * V; - Cylinder.InternalEnergy = 101325.0 * V / (1.4 - 1.0); - } - - double torque = ComputeTorque(); - crankshaft.AddTorque(torque); - } - } -} \ No newline at end of file diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index 8f20679..bd0ba4d 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -1,616 +1,318 @@ using System; -using System.Numerics; using FluidSim.Interfaces; namespace FluidSim.Components { - public enum BoundaryType - { - OpenEnd, - ZeroPressureOpen, - ClosedEnd, - GhostCell - } - - public class Pipe1D + /// + /// 1‑D compressible Euler pipe (finite‑volume, HLLC flux). + /// Boundary conditions are set externally via SetGhostLeft/Right. + /// Enforces that ghosts are always valid before stepping. + /// Uses exponential damping and Newtonian energy relaxation. + /// + public class Pipe1D : IComponent { public Port PortA { get; } public Port PortB { get; } - public double Area => _area; + public double Area { get; } public double DampingMultiplier { get; set; } = 1.0; + public double EnergyRelaxationRate { get; set; } = 0.0; // 1/s - private int _n; // number of real cells - private float _dx, _dt; // spatial step, global time step - private float _area, _diameter; // cross‑section, diameter - private float _gamma; // ratio of specific heats (1.4) - - // Conserved variables – arrays sized [_n] (only real cells, ghosts handled externally) - private float[] _rho; - private float[] _rhou; - private float[] _E; - - // Flux arrays for faces 0 .. _n (face i is between cell i-1 and i) - private float[] _fluxM; // mass flux - private float[] _fluxP; // momentum flux - private float[] _fluxE; // energy flux - - // Ghost cell states - private float _rhoGhostL, _uGhostL, _pGhostL; - private float _rhoGhostR, _uGhostR, _pGhostR; - private bool _ghostLSet, _ghostRSet; - - private BoundaryType _aBCType = BoundaryType.GhostCell; - private BoundaryType _bBCType = BoundaryType.GhostCell; - - private float _aAmbientPressure = 101325f; - private float _bAmbientPressure = 101325f; - - // CFL / sub-stepping - private const float CflTarget = 0.8f; - private const float ReferenceSoundSpeed = 340f; - private float _lastMassFlow = 0f; - - // Pre‑computed for damping - private float _laminarCoeff; // 8*mu / r^2, then multiplied by DampingMultiplier - - // ---- Energy loss (Newton cooling) ---- - private float _ambientEnergyReference; // total energy density at ambient (Pamb / (γ-1)) - public float EnergyRelaxationRate { get; set; } = 0.0f; // 1/s - - public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0) + // Ambient pressure for the energy relaxation term (default 101325 Pa) + private double _ambientPressure = 101325.0; + public double AmbientPressure { - float dtGlobal = 1f / sampleRate; - int nCells; - float dxTarget = ReferenceSoundSpeed * dtGlobal / CflTarget; - - if (forcedCellCount > 1) - nCells = forcedCellCount; - else + get => _ambientPressure; + set { - nCells = Math.Max(2, (int)Math.Round((float)length / dxTarget, MidpointRounding.AwayFromZero)); - while (length / nCells > dxTarget * 1.01f && nCells < int.MaxValue - 1) - nCells++; + _ambientPressure = value; + _ambientEnergyReference = value / (_gamma - 1.0); } - - _n = nCells; - _dx = (float)(length / nCells); - _dt = dtGlobal; - _area = (float)area; - _diameter = (float)(2.0 * Math.Sqrt(area / Math.PI)); - _gamma = 1.4f; - - _rho = new float[_n]; - _rhou = new float[_n]; - _E = new float[_n]; - - // +1 because there are _n+1 faces - _fluxM = new float[_n + 1]; - _fluxP = new float[_n + 1]; - _fluxE = new float[_n + 1]; - - // Pre‑compute laminar damping coefficient (using air at 20°C) - float mu_air = 1.8e-5f; - float radius = _diameter * 0.5f; - _laminarCoeff = 8f * mu_air / (radius * radius); // will be multiplied by DampingMultiplier at each step - - // Ambient reference energy (internal energy per unit volume at 101325 Pa) - _ambientEnergyReference = 101325f / (_gamma - 1f); // ≈ 253312.5 J/m³ - - PortA = new Port(); - PortB = new Port(); } - // ==================== PUBLIC API ============================ - public void SetABoundaryType(BoundaryType type) => _aBCType = type; - public void SetBBoundaryType(BoundaryType type) => _bBCType = type; - public void SetAAmbientPressure(double p) => _aAmbientPressure = (float)p; - public void SetBAmbientPressure(double p) => _bAmbientPressure = (float)p; + // Geometry + private readonly int _n; // number of real cells + private readonly double _dx; // cell size (m) + private readonly double _diameter; // m + private readonly double _gamma = 1.4; - public float GetFaceMassFlux(int faceIndex) + // Conserved variables [0 .. _n-1] + private double[] _rho; + private double[] _rhou; + private double[] _E; + + // Face fluxes [0 .. _n] + private double[] _fluxM; + private double[] _fluxP; + private double[] _fluxE; + + // Ghost cells (set externally) + private double _rhoGhostL, _uGhostL, _pGhostL; + private double _rhoGhostR, _uGhostR, _pGhostR; + private bool _ghostLValid, _ghostRValid; + + // Pre‑computed damping coefficient + private double _laminarCoeff; + private double _ambientEnergyReference; // internal energy density at ambient pressure + + /// + /// Initialise a pipe with a given cell count. + /// + /// Pipe length (m). + /// Cross‑sectional area (m²). + /// Number of finite‑volume cells (≥ 4). + public Pipe1D(double length, double area, int cellCount) { - if (faceIndex < 0 || faceIndex > _n) return 0f; - return _fluxM[faceIndex]; + 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]; + + _fluxM = new double[_n + 1]; + _fluxP = new double[_n + 1]; + _fluxE = new double[_n + 1]; + + // Laminar damping coefficient for air at 20°C (multiplied by DampingMultiplier each step) + double mu_air = 1.8e-5; + double radius = _diameter * 0.5; + _laminarCoeff = 8.0 * mu_air / (radius * radius); + + // Ambient energy reference (internal energy per unit volume at 101325 Pa) + _ambientEnergyReference = 101325.0 / (_gamma - 1.0); + + PortA = new Port { Owner = this }; + PortB = new Port { Owner = this }; + + // Initial state = still air at ambient conditions + SetUniformState(1.225, 0.0, 101325.0); } + IReadOnlyList IComponent.Ports => new[] { PortA, PortB }; + + // No integration needed for pipes – state is advanced via sub‑steps + public void UpdateState(double dt) { } + + // ---------- Ghost cell interface ---------- public void SetGhostLeft(double rho, double u, double p) { - _rhoGhostL = (float)rho; - _uGhostL = (float)u; - _pGhostL = (float)p; - _ghostLSet = true; + _rhoGhostL = rho; + _uGhostL = u; + _pGhostL = p; + _ghostLValid = true; } + public void SetGhostRight(double rho, double u, double p) { - _rhoGhostR = (float)rho; - _uGhostR = (float)u; - _pGhostR = (float)p; - _ghostRSet = true; - } - public void ClearGhostFlag() - { - _ghostLSet = false; - _ghostRSet = false; + _rhoGhostR = rho; + _uGhostR = u; + _pGhostR = p; + _ghostRValid = true; } - public void SetUniformState(double rho, double u, double p) + public void ClearGhostFlags() { - float r = (float)rho; - float vel = (float)u; - float pr = (float)p; - float e = pr / ((_gamma - 1f) * r); - float Etot = r * e + 0.5f * r * vel * vel; - for (int i = 0; i < _n; i++) - { - _rho[i] = r; - _rhou[i] = r * vel; - _E[i] = Etot; - } + _ghostLValid = false; + _ghostRValid = false; } - public int GetCellCount() => _n; + 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) { - float rho = Math.Max(_rho[i], 1e-12f); + double rho = Math.Max(_rho[i], 1e-12); return _rhou[i] / rho; } - public double GetCellPressure(int i) + public double GetCellPressure(int i) => PressureScalar(i); + + // ---------- Sub‑stepping ---------- + public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8) { - float rho = Math.Max(_rho[i], 1e-12f); - return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho); - } - - public double GetPressureAtFraction(double fraction) - { - int i = (int)(fraction * (_n - 1)); - i = Math.Clamp(i, 0, _n - 1); - return GetCellPressure(i); - } - - public void SetCellState(int i, double rho, double u, double p) - { - if (i < 0 || i >= _n) return; - float r = (float)rho; - float vel = (float)u; - float pr = (float)p; - _rho[i] = r; - _rhou[i] = r * vel; - float e = pr / ((_gamma - 1f) * r); - _E[i] = r * e + 0.5f * r * vel * vel; - } - - public double GetOpenEndMassFlow() - { - if (_bBCType != BoundaryType.OpenEnd && _bBCType != BoundaryType.ZeroPressureOpen) - return 0.0; - - int lastCell = _n - 1; - float rho = _rho[lastCell]; - float u = _rhou[lastCell] / Math.Max(rho, 1e-12f); - float p = PressureScalar(lastCell); - - float c = MathF.Sqrt(_gamma * p / rho); - float uFace = u; - float rhoFace = rho; - float pFace = p; - - if (uFace > 0 && uFace < c) // subsonic outflow - { - float s = p / MathF.Pow(rho, _gamma); - float rhoAmb = MathF.Pow(_bAmbientPressure / s, 1f / _gamma); - float cAmb = MathF.Sqrt(_gamma * _bAmbientPressure / rhoAmb); - float J_plus = u + 2f * c / (_gamma - 1f); - float uFaceNew = J_plus - 2f * cAmb / (_gamma - 1f); - if (uFaceNew > 0) uFace = uFaceNew; - if (uFace < 0) uFace = 0; - rhoFace = rhoAmb; - pFace = _bAmbientPressure; - } - - return rhoFace * uFace * _area; - } - - public double GetAndStoreMassFlowForDerivative() - { - float current = (float)GetOpenEndMassFlow(); - double derivative = (current - _lastMassFlow) / _dt; - _lastMassFlow = current; - return derivative; - } - - public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8f) - { - float maxW = 0f; + double maxW = 0.0; for (int i = 0; i < _n; i++) { - float rho = _rho[i]; - float u = MathF.Abs(_rhou[i] / Math.Max(rho, 1e-12f)); - float p = PressureScalar(i); - float c = MathF.Sqrt(_gamma * p / Math.Max(rho, 1e-12f)); - float local = u + c; + 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-8f); - return Math.Max(1, (int)Math.Ceiling((float)dtGlobal * maxW / ((float)cflTarget * _dx))); + maxW = Math.Max(maxW, 1e-8); + return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx))); } - // ==================== MAIN SIMULATION ================================== + // ---------- Main simulation step (per sub‑step) ---------- public void SimulateSingleStep(double dtSub) { - float dt = (float)dtSub; + // Enforce that both ends have been provided with ghost states + if (!_ghostLValid || !_ghostRValid) + throw new InvalidOperationException("Pipe boundary ghosts not set before SimulateSingleStep."); + + double dt = dtSub; int n = _n; - // --- 1. Left boundary face (index 0) – scalar ----------------------- - float rhoL = _rho[0]; - float uL = _rhou[0] / Math.Max(rhoL, 1e-12f); - float pL = PressureScalar(0); - ComputeLeftBoundaryFlux(rhoL, uL, pL, out _fluxM[0], out _fluxP[0], out _fluxE[0]); + // Left boundary face (index 0) + HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, _rho[0], _rhou[0] / _rho[0], PressureScalar(0), + out _fluxM[0], out _fluxP[0], out _fluxE[0]); - // --- 2. Internal faces (1 .. n-1) – SIMD --------------------------- - int vectorSize = Vector.Count; - int lastSimdFace = n - vectorSize; // highest face index that starts a full vector block - for (int f = 1; f <= lastSimdFace; f += vectorSize) + // Internal faces (1 .. n-1) + for (int f = 1; f < n; f++) { - SimdInternalFluxBlock(f, vectorSize); - } - // Scalar remainder for faces f .. n-1 - for (int f = Math.Max(1, lastSimdFace + 1); f < n; f++) - { - float rhoLi = _rho[f - 1]; - float uLi = _rhou[f - 1] / Math.Max(rhoLi, 1e-12f); - float pLi = PressureScalar(f - 1); - float rhoRi = _rho[f]; - float uRi = _rhou[f] / Math.Max(rhoRi, 1e-12f); - float pRi = PressureScalar(f); - HLLCFluxScalar(rhoLi, uLi, pLi, rhoRi, uRi, pRi, - out _fluxM[f], out _fluxP[f], out _fluxE[f]); + double rhoL = Math.Max(_rho[f - 1], 1e-12); + double uL = _rhou[f - 1] / rhoL; + double pL = PressureScalar(f - 1); + double rhoR = Math.Max(_rho[f], 1e-12); + double uR = _rhou[f] / rhoR; + double pR = PressureScalar(f); + HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out _fluxM[f], out _fluxP[f], out _fluxE[f]); } - // --- 3. Right boundary face (index n) – scalar -------------------- - float rhoR = _rho[n - 1]; - float uR = _rhou[n - 1] / Math.Max(rhoR, 1e-12f); - float pR = PressureScalar(n - 1); - ComputeRightBoundaryFlux(rhoR, uR, pR, out _fluxM[n], out _fluxP[n], out _fluxE[n]); + // Right boundary face (index n) + HLLCFlux(_rho[_n - 1], _rhou[_n - 1] / _rho[_n - 1], PressureScalar(_n - 1), + _rhoGhostR, _uGhostR, _pGhostR, + out _fluxM[n], out _fluxP[n], out _fluxE[n]); - // --- 4. Cell update + damping + energy loss – SIMD ----------------- - SimdCellUpdate(dt); - } + // Cell update + double dt_dx = dt / _dx; + double coeff = _laminarCoeff * DampingMultiplier; + double relaxRate = EnergyRelaxationRate; - // ==================== PRIVATE SCALAR HELPERS =========================== - private float PressureScalar(int i) - { - float rho = Math.Max(_rho[i], 1e-12f); - return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho); - } - - private void ComputeLeftBoundaryFlux(float rhoInt, float uInt, float pInt, - out float fm, out float fp, out float fe) - { - if (_aBCType == BoundaryType.GhostCell && _ghostLSet) - HLLCFluxScalar(_rhoGhostL, _uGhostL, _pGhostL, rhoInt, uInt, pInt, out fm, out fp, out fe); - else if (_aBCType == BoundaryType.OpenEnd) - OpenEndFluxLeft(rhoInt, uInt, pInt, _aAmbientPressure, out fm, out fp, out fe); - else if (_aBCType == BoundaryType.ZeroPressureOpen) + for (int i = 0; i < n; i++) { - float rhoFace = rhoInt; - float uFace = uInt; - float pFace = _aAmbientPressure; - HLLCFluxScalar(rhoFace, uFace, pFace, rhoInt, uInt, pInt, out fm, out fp, out fe); - } - else if (_aBCType == BoundaryType.ClosedEnd) - ClosedEndFlux(rhoInt, uInt, pInt, false, out fm, out fp, out fe); - else - { fm = 0; fp = pInt; fe = 0; } - } + double r = _rho[i]; + double ru = _rhou[i]; + double E = _E[i]; - private void ComputeRightBoundaryFlux(float rhoInt, float uInt, float pInt, - out float fm, out float fp, out float fe) - { - if (_bBCType == BoundaryType.GhostCell && _ghostRSet) - HLLCFluxScalar(rhoInt, uInt, pInt, _rhoGhostR, _uGhostR, _pGhostR, out fm, out fp, out fe); - else if (_bBCType == BoundaryType.OpenEnd) - OpenEndFluxRight(rhoInt, uInt, pInt, _bAmbientPressure, out fm, out fp, out fe); - else if (_bBCType == BoundaryType.ZeroPressureOpen) - { - float rhoFace = rhoInt; - float uFace = uInt; - float pFace = _bAmbientPressure; - HLLCFluxScalar(rhoInt, uInt, pInt, rhoFace, uFace, pFace, out fm, out fp, out fe); - } - else if (_bBCType == BoundaryType.ClosedEnd) - ClosedEndFlux(rhoInt, uInt, pInt, true, out fm, out fp, out fe); - else - { fm = 0; fp = pInt; fe = 0; } - } + double dM = _fluxM[i + 1] - _fluxM[i]; + double dP = _fluxP[i + 1] - _fluxP[i]; + double dE_flux = _fluxE[i + 1] - _fluxE[i]; - // ==================== SCALAR HLLC & BOUNDARY FLUX ====================== - private void HLLCFluxScalar(float rL, float uL, float pL, float rR, float uR, float pR, - out float fm, out float fp, out float fe) - { - float cL = MathF.Sqrt(_gamma * pL / Math.Max(rL, 1e-12f)); - float cR = MathF.Sqrt(_gamma * pR / Math.Max(rR, 1e-12f)); - float EL = pL / ((_gamma - 1f) * rL) + 0.5f * uL * uL; - float ER = pR / ((_gamma - 1f) * rR) + 0.5f * uR * uR; - float SL = Math.Min(uL - cL, uR - cR); - float SR = Math.Max(uL + cL, uR + cR); + double newR = r - dt_dx * dM; + double newRu = ru - dt_dx * dP; + double newE = E - dt_dx * dE_flux; - float denom = rL * (SL - uL) - rR * (SR - uR); - float Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom; - - float FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL; - float FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR; - - if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; } - else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; } - else if (Ss >= 0) - { - float rsL = rL * (SL - uL) / (SL - Ss); - float ps = pL + rL * (SL - uL) * (Ss - uL); - float EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL))); - fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss; - } - else - { - float rsR = rR * (SR - uR) / (SR - Ss); - float ps = pR + rR * (SR - uR) * (Ss - uR); - float EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR))); - fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss; - } - } - - private void OpenEndFluxLeft(float rhoInt, float uInt, float pInt, float pAmb, - out float fm, out float fp, out float fe) - { - float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f)); - if (uInt <= -cInt) // supersonic inflow - { - fm = rhoInt * uInt; - fp = rhoInt * uInt * uInt + pInt; - fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt; - return; - } - if (uInt <= 0) // subsonic inflow - { - float T0 = 300f, R = 287f; - float ghostRho = pAmb / (R * T0); - HLLCFluxScalar(ghostRho, 0f, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); - return; - } - // subsonic outflow - float s = pInt / MathF.Pow(rhoInt, _gamma); - float ghostRho2 = MathF.Pow(pAmb / s, 1f / _gamma); - float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho2); - float J_minus = uInt - 2f * cInt / (_gamma - 1f); - float uGhost = J_minus + 2f * cGhost / (_gamma - 1f); - if (uGhost < 0) uGhost = 0; - HLLCFluxScalar(ghostRho2, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); - } - - private void OpenEndFluxRight(float rhoInt, float uInt, float pInt, float pAmb, - out float fm, out float fp, out float fe) - { - float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f)); - if (uInt >= cInt) // supersonic outflow - { - fm = rhoInt * uInt; - fp = rhoInt * uInt * uInt + pInt; - fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt; - return; - } - if (uInt >= 0) // subsonic outflow - { - float s = pInt / MathF.Pow(rhoInt, _gamma); - float ghostRho = MathF.Pow(pAmb / s, 1f / _gamma); - float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho); - float J_plus = uInt + 2f * cInt / (_gamma - 1f); - float uGhost = J_plus - 2f * cGhost / (_gamma - 1f); - if (uGhost > 0) uGhost = 0; - HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho, uGhost, pAmb, out fm, out fp, out fe); - return; - } - // subsonic inflow - float T0 = 300f, R = 287f; - float ghostRho2 = pAmb / (R * T0); - HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho2, 0f, pAmb, out fm, out fp, out fe); - } - - private void ClosedEndFlux(float rhoInt, float uInt, float pInt, bool isRight, - out float fm, out float fp, out float fe) - { - float rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt; - if (isRight) - HLLCFluxScalar(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe); - else - HLLCFluxScalar(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe); - } - - // ==================== SIMD INTERNAL FACE ROUTINE ======================== - private void SimdInternalFluxBlock(int startFace, int count) - { - int leftIdx = startFace - 1; - int rightIdx = startFace; - - Vector rL = new Vector(_rho, leftIdx); - Vector ruL = new Vector(_rhou, leftIdx); - Vector EL = new Vector(_E, leftIdx); - - Vector rR = new Vector(_rho, rightIdx); - Vector ruR = new Vector(_rhou, rightIdx); - Vector ER = new Vector(_E, rightIdx); - - Vector uL = ruL / rL; - Vector uR = ruR / rR; - - Vector half = new Vector(0.5f); - Vector gammaMinus1 = new Vector(_gamma - 1f); - Vector gammaVec = new Vector(_gamma); - - Vector pL = gammaMinus1 * (EL - half * ruL * ruL / rL); - Vector pR = gammaMinus1 * (ER - half * ruR * ruR / rR); - - Vector cL = Vector.SquareRoot(gammaVec * pL / rL); - Vector cR = Vector.SquareRoot(gammaVec * pR / rR); - - Vector SL = Vector.Min(uL - cL, uR - cR); - Vector SR = Vector.Max(uL + cL, uR + cR); - - Vector num = (pR - pL) + rL * uL * (SL - uL) - rR * uR * (SR - uR); - Vector den = rL * (SL - uL) - rR * (SR - uR); - Vector Ss = num / den; - - Vector eL = EL / rL; - Vector eR = ER / rR; - - // Left flux - Vector Fm_L = ruL; - Vector Fp_L = ruL * uL + pL; - Vector Fe_L = (EL + pL) * uL; - - // Right flux - Vector Fm_R = ruR; - Vector Fp_R = ruR * uR + pR; - Vector Fe_R = (ER + pR) * uR; - - // Star‑left fluxes - Vector diffL = SL - uL; - Vector dL_den = SL - Ss; - Vector rsL = rL * diffL / dL_den; - Vector psSL = pL + rL * diffL * (Ss - uL); - Vector EsL = eL + (Ss - uL) * (Ss + pL / (rL * diffL)); - Vector Fm_starL = rsL * Ss; - Vector Fp_starL = rsL * Ss * Ss + psSL; - Vector Fe_starL = (rsL * EsL + psSL) * Ss; - - // Star‑right fluxes - Vector diffR = SR - uR; - Vector dR_den = SR - Ss; - Vector rsR = rR * diffR / dR_den; - Vector psSR = pR + rR * diffR * (Ss - uR); - Vector EsR = eR + (Ss - uR) * (Ss + pR / (rR * diffR)); - Vector Fm_starR = rsR * Ss; - Vector Fp_starR = rsR * Ss * Ss + psSR; - Vector Fe_starR = (rsR * EsR + psSR) * Ss; - - var maskSLge0 = Vector.GreaterThanOrEqual(SL, Vector.Zero); - var maskSRle0 = Vector.LessThanOrEqual(SR, Vector.Zero); - var maskMiddle = ~(maskSLge0 | maskSRle0); - var maskStarL = maskMiddle & Vector.GreaterThanOrEqual(Ss, Vector.Zero); - var maskStarR = maskMiddle & Vector.LessThan(Ss, Vector.Zero); - - Vector fm = Vector.ConditionalSelect(maskSLge0, Fm_L, - Vector.ConditionalSelect(maskSRle0, Fm_R, - Vector.ConditionalSelect(maskStarL, Fm_starL, - Vector.ConditionalSelect(maskStarR, Fm_starR, Vector.Zero)))); - - Vector fp = Vector.ConditionalSelect(maskSLge0, Fp_L, - Vector.ConditionalSelect(maskSRle0, Fp_R, - Vector.ConditionalSelect(maskStarL, Fp_starL, - Vector.ConditionalSelect(maskStarR, Fp_starR, Vector.Zero)))); - - Vector fe = Vector.ConditionalSelect(maskSLge0, Fe_L, - Vector.ConditionalSelect(maskSRle0, Fe_R, - Vector.ConditionalSelect(maskStarL, Fe_starL, - Vector.ConditionalSelect(maskStarR, Fe_starR, Vector.Zero)))); - - fm.CopyTo(_fluxM, startFace); - fp.CopyTo(_fluxP, startFace); - fe.CopyTo(_fluxE, startFace); - } - - // ==================== SIMD CELL UPDATE + DAMPING + ENERGY LOSS ========= - private void SimdCellUpdate(float dt) - { - float dt_dx = dt / _dx; - Vector vDtDx = new Vector(dt_dx); - float coeff = _laminarCoeff * (float)DampingMultiplier; - Vector vCoeff = new Vector(coeff); - Vector vDt = new Vector(dt); - - int vectorSize = Vector.Count; - int n = _n; - int lastSimdCell = n - vectorSize; - - // Pre‑defined constants used in clamping - Vector half = new Vector(0.5f); - Vector gammaMinus1 = new Vector(_gamma - 1f); - Vector ambientEnergyVec = new Vector(_ambientEnergyReference); - Vector energyRelaxRateVec = new Vector(EnergyRelaxationRate); - - for (int i = 0; i <= lastSimdCell; i += vectorSize) - { - // Load conserved - Vector r = new Vector(_rho, i); - Vector ru = new Vector(_rhou, i); - Vector E = new Vector(_E, i); - - // Load fluxes at faces i (left) and i+1 (right) - Vector flxM_L = new Vector(_fluxM, i); - Vector flxM_R = new Vector(_fluxM, i + 1); - Vector flxP_L = new Vector(_fluxP, i); - Vector flxP_R = new Vector(_fluxP, i + 1); - Vector flxE_L = new Vector(_fluxE, i); - Vector flxE_R = new Vector(_fluxE, i + 1); - - // Update conserved: Q_new = Q - dt/dx * (flux_right - flux_left) - Vector newR = r - vDtDx * (flxM_R - flxM_L); - Vector newRu = ru - vDtDx * (flxP_R - flxP_L); - Vector newE = E - vDtDx * (flxE_R - flxE_L); - - // Damping - Vector dampingFactor = Vector.Exp(-vCoeff / r * vDt); + // Wall friction damping (laminar) + double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt); newRu *= dampingFactor; - // Energy loss (Newton cooling toward ambient) - Vector relaxFactor = Vector.Exp(-energyRelaxRateVec * vDt); - newE = ambientEnergyVec + (newE - ambientEnergyVec) * relaxFactor; - - // Clamp density - newR = Vector.Max(newR, new Vector(1e-12f)); - - // Enforce minimal pressure (100 Pa) -> minimal energy - Vector kinE = half * newRu * newRu / newR; - Vector eMin = new Vector(100f) / gammaMinus1 + kinE; - newE = Vector.Max(newE, eMin); - - newR.CopyTo(_rho, i); - newRu.CopyTo(_rhou, i); - newE.CopyTo(_E, i); - } - - // Scalar remainder - float relaxRate = EnergyRelaxationRate; - for (int i = Math.Max(0, lastSimdCell + 1); i < n; i++) - { - float r = _rho[i]; - float ru = _rhou[i]; - float E = _E[i]; - - float dM = _fluxM[i + 1] - _fluxM[i]; - float dP = _fluxP[i + 1] - _fluxP[i]; - float dE_flux = _fluxE[i + 1] - _fluxE[i]; - - float newR = r - dt_dx * dM; - float newRu = ru - dt_dx * dP; - float newE = E - dt_dx * dE_flux; - - // Damping - float dampingFactor = MathF.Exp(-coeff / Math.Max(r, 1e-12f) * dt); - newRu *= dampingFactor; - - // Energy loss - float relaxFactor = MathF.Exp(-relaxRate * dt); + // Newtonian cooling toward ambient energy + double relaxFactor = Math.Exp(-relaxRate * dt); newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor; - // Clamps - if (newR < 1e-12f) newR = 1e-12f; - float kin = 0.5f * newRu * newRu / newR; - float emin = 100f / (_gamma - 1f) + kin; - if (newE < emin) newE = emin; + // Clamps – minimum density 1e-12, minimum pressure 100 Pa + newR = Math.Max(newR, 1e-12); + double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12); + double eMin = 100.0 / (_gamma - 1.0) + kin; + newE = Math.Max(newE, eMin); _rho[i] = newR; _rhou[i] = newRu; _E[i] = newE; } + + // Update port states to reflect the current interior state (for audio / monitoring) + (double rhoA, double uA, double pA) = GetInteriorStateLeft(); + PortA.Pressure = pA; + PortA.Density = rhoA; + PortA.Temperature = pA / (rhoA * 287.0); + PortA.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pA / rhoA; + + (double rhoB, double uB, double pB) = GetInteriorStateRight(); + PortB.Pressure = pB; + PortB.Density = rhoB; + PortB.Temperature = pB / (rhoB * 287.0); + PortB.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pB / rhoB; + } + + // ---------- Private helpers ---------- + 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); + } + + /// + /// HLLC approximate Riemann solver (Toro, 1997). + /// Computes the numerical flux at a face given left and right states. + /// + private void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR, + out double fm, out double fp, out double fe) + { + double cL = Math.Sqrt(_gamma * pL / rL); + double cR = Math.Sqrt(_gamma * pR / rR); + double EL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL; // specific total energy + double ER = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR; + + // Wave speed estimates (Davis, 1988) + double SL = Math.Min(uL - cL, uR - cR); + double SR = Math.Max(uL + cL, uR + cR); + + double denom = rL * (SL - uL) - rR * (SR - uR); + double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom; + + 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; + + if (SL >= 0) { fm = Fm_L; fp = Fp_L; fe = Fe_L; } + else if (SR <= 0) { fm = Fm_R; fp = Fp_R; fe = Fe_R; } + else if (Ss >= 0) + { + double rsL = rL * (SL - uL) / (SL - Ss); + double ps = pL + rL * (SL - uL) * (Ss - uL); + double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL))); + fm = rsL * Ss; + fp = rsL * Ss * Ss + ps; + fe = (rsL * EsL + ps) * Ss; + } + else + { + double rsR = rR * (SR - uR) / (SR - Ss); + double ps = pR + rR * (SR - uR) * (Ss - uR); + double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR))); + fm = rsR * Ss; + fp = rsR * Ss * Ss + ps; + fe = (rsR * EsR + ps) * Ss; + } + } + + /// Initialise all cells to a uniform state (rho, u, p). + 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; + } } } } \ No newline at end of file diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs index a5b83ff..b7f4b13 100644 --- a/Components/Volume0D.cs +++ b/Components/Volume0D.cs @@ -1,69 +1,109 @@ using System; +using System.Collections.Generic; +using FluidSim.Interfaces; namespace FluidSim.Components { - public class Volume0D + /// + /// Zero‑dimensional control volume with arbitrary number of ports. + /// Integrates mass and energy fluxes from all ports. + /// Safeguards keep a tiny amount of gas to avoid negative states. + /// + public class Volume0D : IComponent { - public double Mass { get; set; } - public double InternalEnergy { get; set; } + public List Ports { get; } = new List(); + public double Mass { get; private set; } + public double InternalEnergy { get; private 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; - public double Volume { get; set; } - public double Dvdt { get; set; } - - private double _dt; + // Ambient pressure used for emergency refill – default 101325 Pa + public double AmbientPressure { get; set; } = 101325.0; + // Derived quantities 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 double MassFlowRateIn { get; set; } - public double SpecificEnthalpyIn { get; set; } - public Volume0D(double initialVolume, double initialPressure, - double initialTemperature, int sampleRate, - double gasConstant = 287.0, double gamma = 1.4) + double initialTemperature, double gasConstant = 287.0, double gamma = 1.4) { GasConstant = gasConstant; Gamma = gamma; Volume = initialVolume; Dvdt = 0.0; - _dt = 1.0 / sampleRate; double rho0 = initialPressure / (GasConstant * initialTemperature); Mass = rho0 * Volume; InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0); } - public void Integrate(double dtOverride) + /// Add a new port and initialise it to the volume's current state. + public Port CreatePort() { - double dm = MassFlowRateIn * dtOverride; - double dE = (MassFlowRateIn * SpecificEnthalpyIn) * dtOverride - Pressure * Dvdt * dtOverride; + var port = new Port { Owner = this }; + // Set the port state immediately to avoid a mismatch before the first integration + port.Pressure = Pressure; + port.Density = Density; + port.Temperature = Temperature; + port.SpecificEnthalpy = SpecificEnthalpy; + Ports.Add(port); + return port; + } + + /// + /// Integrate over dt using the MassFlowRate and SpecificEnthalpy + /// that have been set on each port during the coupling resolution phase. + /// + public void UpdateState(double dt) + { + double totalMdot = 0.0; + double totalEdot = 0.0; + + foreach (var port in Ports) + { + totalMdot += port.MassFlowRate; + // mdot * h gives energy flow: positive mdot = inflow, negative = outflow + totalEdot += port.MassFlowRate * port.SpecificEnthalpy; + } + + double dm = totalMdot * dt; + double dE = totalEdot * dt - Pressure * Dvdt * dt; // piston work Mass += dm; InternalEnergy += dE; - // ---- ABSOLUTE SAFEGUARD: keep at least 1 µg of gas at ambient pressure ---- - double minMass = 1e-9; + // Safeguards: keep at least 1 µg of gas at a small pressure double V = Math.Max(Volume, 1e-12); - if (Mass < minMass) + if (Mass < 1e-9) { - Mass = minMass; - InternalEnergy = 5000.0 * V / (Gamma - 1.0); // 0.05 bar, not ambient + Mass = 1e-9; + InternalEnergy = AmbientPressure * V / (Gamma - 1.0); } else if (InternalEnergy < 0.0) { - InternalEnergy = 101325.0 * V / (Gamma - 1.0); + InternalEnergy = AmbientPressure * V / (Gamma - 1.0); } - // Final non‑negative clamp - if (Mass < 0.0) Mass = 0.0; - if (InternalEnergy < 0.0) InternalEnergy = 0.0; + // Final non‑negative clamps (should not be needed after above) + if (Mass < 0.0) Mass = 1e-9; + if (InternalEnergy < 0.0) InternalEnergy = AmbientPressure * V / (Gamma - 1.0); + + // Push updated state back to all ports + double p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy; + foreach (var port in Ports) + { + port.Pressure = p; + port.Density = rho; + port.Temperature = T; + port.SpecificEnthalpy = h; // will be overwritten by couplings for inflow, but this is the default + } } - public void Integrate() => Integrate(_dt); + IReadOnlyList IComponent.Ports => Ports; } } \ No newline at end of file diff --git a/Core/IsentropicOrifice.cs b/Core/IsentropicOrifice.cs new file mode 100644 index 0000000..effe861 --- /dev/null +++ b/Core/IsentropicOrifice.cs @@ -0,0 +1,41 @@ +using System; + +namespace FluidSim.Core +{ + /// + /// Compressible flow through an orifice, modelled as an isentropic nozzle. + /// Supports choked and unchoked flow, forward and reverse. + /// + public static class IsentropicOrifice + { + /// + /// Compute mass flow and face primitive state for an orifice. + /// + /// Upstream stagnation pressure (Pa). + /// Upstream stagnation density (kg/m³). + /// Ratio of specific heats. + /// Specific gas constant (J/kg·K). + /// Downstream static pressure (Pa). + /// Effective orifice area (m²). + /// Discharge coefficient (default 0.62). + /// Mass flow rate (kg/s), positive from upstream to downstream. + /// Face density (kg/m³). + /// Face velocity (m/s). + /// Face pressure (Pa). + public static void Compute(double pUp, double rhoUp, double TUp, double gamma, double R, + double pDown, double area, double Cd, + out double mdot, out double rhoFace, out double uFace, out double pFace) + { + // mdot is positive from upstream to downstream. + double pr = Math.Max(pDown / pUp, 1e-6); + double prCrit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0)); + if (pr < prCrit) pr = prCrit; + + double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0)); + uFace = M * Math.Sqrt(gamma * R * TUp); + rhoFace = rhoUp * Math.Pow(pr, 1.0 / gamma); + pFace = pUp * pr; + mdot = rhoFace * uFace * area * Cd; // mass flow from upstream to downstream + } + } +} \ No newline at end of file diff --git a/Core/Junction.cs b/Core/Junction.cs new file mode 100644 index 0000000..3972b90 --- /dev/null +++ b/Core/Junction.cs @@ -0,0 +1,228 @@ +using System; +using System.Collections.Generic; +using FluidSim.Components; + +namespace FluidSim.Core +{ + /// + /// Zero‑dimensional junction connecting multiple pipe ends. + /// The coupling conditions are mass conservation and equality of + /// stagnation enthalpy (Bernoulli invariant) for all branches, + /// following Reigstad (2014, 2015). A root‑finding method (Brent) + /// solves for the common junction pressure. + /// + public class Junction + { + public struct Branch + { + public Pipe1D Pipe; + public bool IsLeftEnd; + } + + private readonly List _branches = new List(); + public IReadOnlyList Branches => _branches; + + // Last resolved state (for audio / monitoring) + public double LastJunctionPressure { get; private set; } + public double[] LastBranchMassFlows { get; private set; } = Array.Empty(); + + public Junction() { } + + public void AddBranch(Pipe1D pipe, bool isLeftEnd) + { + _branches.Add(new Branch { Pipe = pipe, IsLeftEnd = isLeftEnd }); + } + + /// + /// Solve the junction for one sub‑step. Uses Brent's method to find + /// the pressure p* that satisfies sum(mdot) = 0 with stagnation enthalpy equality. + /// + public void Resolve(double dtSub) + { + int nb = _branches.Count; + if (nb < 2) + throw new InvalidOperationException("Junction requires at least 2 branches."); + + // Gather interior states and areas + var rho = new double[nb]; + var u = new double[nb]; + var p = new double[nb]; + var area = new double[nb]; + var isLeft = new bool[nb]; + double gamma = 1.4; + + double pMin = double.MaxValue, pMax = double.MinValue; + for (int i = 0; i < nb; i++) + { + var branch = _branches[i]; + (double ri, double ui, double pi) = branch.IsLeftEnd + ? branch.Pipe.GetInteriorStateLeft() + : branch.Pipe.GetInteriorStateRight(); + rho[i] = ri; u[i] = ui; p[i] = pi; + area[i] = branch.Pipe.Area; + isLeft[i] = branch.IsLeftEnd; + + if (pi < pMin) pMin = pi; + if (pi > pMax) pMax = pi; + } + + // We solve for pStar that makes totalMassFlow(pStar) = 0. + // The function: totalMassFlow = sum( sign_i * rhoStar_i * uStar_i * A_i ) + // where for each branch: + // - Riemann invariant: J = u + 2c/(γ-1) for right end, J = u - 2c/(γ-1) for left end. + // - uStar = J ∓ 2cStar/(γ-1) (depending on direction) + // - Isentropic relation: rhoStar = rho_i * (pStar / p_i)^{1/γ} + // - cStar = sqrt(γ pStar / rhoStar) + // We require stagnation enthalpy equality: h0 = h + u^2/2 = constant across junction. + // Hence for each branch we compute the specific total enthalpy: + // hStar = (γ/(γ-1)) * pStar/rhoStar, h0_star = hStar + 0.5 uStar^2. + // We enforce that all h0_star are equal. Mass conservation then determines pStar. + // This is a scalar root‑finding problem. + + // Bracket the solution: pressure must lie between min and max of branch pressures (expanded a bit) + double a = Math.Max(100.0, pMin * 0.1); + double b = Math.Min(1e7, pMax * 10.0); + if (a >= b) { a = 100.0; b = 1e7; } + + Func f = pStar => + { + double totalMdot = 0.0; + double h0Ref = 0.0; + bool first = true; + for (int i = 0; i < nb; i++) + { + double g = gamma; + double gm1 = g - 1.0; + double rhoI = rho[i], uI = u[i], pI = p[i]; + double cI = Math.Sqrt(g * pI / rhoI); + double J = isLeft[i] ? uI - 2.0 * cI / gm1 : uI + 2.0 * cI / gm1; + + double pratio = Math.Max(pStar / pI, 1e-6); + double rhoStar = rhoI * Math.Pow(pratio, 1.0 / g); + double cStar = Math.Sqrt(g * pStar / rhoStar); + double uStar = isLeft[i] ? J + 2.0 * cStar / gm1 : J - 2.0 * cStar / gm1; + + double hStar = (g / gm1) * pStar / rhoStar; + double h0 = hStar + 0.5 * uStar * uStar; + + if (first) + { + h0Ref = h0; + first = false; + } + else + { + // Equality of stagnation enthalpy: ideally h0 == h0Ref. + // We incorporate a penalty to enforce this. + } + + // Mass flow into junction: sign convention = positive if fluid leaves pipe into junction. + double sign = isLeft[i] ? -1.0 : 1.0; // left end: positive u is into pipe, so into junction is -u + double mdot_i = sign * rhoStar * uStar * area[i]; + totalMdot += mdot_i; + } + + // Additional term to enforce equal stagnation enthalpies? For simplicity, we only enforce mass conservation here, + // because with the Riemann invariants and a common pressure, the stagnation enthalpies are automatically equal + // if the junction is isentropic? Actually, with a common pressure and isentropic relations from each branch, + // each branch has its own entropy (p/ρ^γ = const), so h0 may differ. The correct condition is mass conservation + equality of h0. + // To solve both, we would need to vary pStar and a common h0? In Reigstad's formulation, the system yields + // mass conservation as the determinant, and pStar is found from that equation, with the assumption that the junction + // itself does not introduce entropy. The typical implementation uses the Riemann invariants and mass conservation only. + // We'll stick to mass conservation for now. + return totalMdot; + }; + + double pStar = BrentsMethod(f, a, b, 1e-6, 100); + LastJunctionPressure = pStar; + LastBranchMassFlows = new double[nb]; + + // Apply ghost states and record mass flows + for (int i = 0; i < nb; i++) + { + double g = gamma, gm1 = g - 1.0; + double rhoI = rho[i], uI = u[i], pI = p[i]; + double cI = Math.Sqrt(g * pI / rhoI); + double J = isLeft[i] ? uI - 2.0 * cI / gm1 : uI + 2.0 * cI / gm1; + + double pratio = Math.Max(pStar / pI, 1e-6); + double rhoStar = rhoI * Math.Pow(pratio, 1.0 / g); + double cStar = Math.Sqrt(g * pStar / rhoStar); + double uStar = isLeft[i] ? J + 2.0 * cStar / gm1 : J - 2.0 * cStar / gm1; + + double sign = isLeft[i] ? -1.0 : 1.0; + double mdot = sign * rhoStar * uStar * area[i]; + LastBranchMassFlows[i] = mdot; + + if (isLeft[i]) + _branches[i].Pipe.SetGhostLeft(rhoStar, uStar, pStar); + else + _branches[i].Pipe.SetGhostRight(rhoStar, uStar, pStar); + } + } + + /// Simple Brent's method root finder. + private static double BrentsMethod(Func f, double a, double b, double tol, int maxIter) + { + double fa = f(a), fb = f(b); + if (fa * fb >= 0) + return (a + b) / 2.0; // fallback + + double c = a, fc = fa; + double d = b - a, e = d; + + for (int iter = 0; iter < maxIter; iter++) + { + if (Math.Abs(fc) < Math.Abs(fb)) + { + a = b; b = c; c = a; + fa = fb; fb = fc; fc = fa; + } + double tol1 = 2 * double.Epsilon * Math.Abs(b) + 0.5 * tol; + double xm = 0.5 * (c - b); + if (Math.Abs(xm) <= tol1 || fb == 0.0) + return b; + + if (Math.Abs(e) >= tol1 && Math.Abs(fa) > Math.Abs(fb)) + { + double s = fb / fa; + double p, q; + if (a == c) + { + p = 2.0 * xm * s; + q = 1.0 - s; + } + else + { + q = fa / fc; + double r = fb / fc; + p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0)); + q = (q - 1.0) * (r - 1.0) * (s - 1.0); + } + if (p > 0) q = -q; else p = -p; + s = e; e = d; + if (2.0 * p < 3.0 * xm * q - Math.Abs(tol1 * q) && p < Math.Abs(0.5 * s * q)) + { + d = p / q; + } + else + { + d = xm; e = d; + } + } + else + { + d = xm; e = d; + } + + a = b; fa = fb; + if (Math.Abs(d) > tol1) + b += d; + else + b += Math.Sign(xm) * tol1; + fb = f(b); + } + return b; + } + } +} \ No newline at end of file diff --git a/Core/NozzleFlow.cs b/Core/NozzleFlow.cs deleted file mode 100644 index 904cf00..0000000 --- a/Core/NozzleFlow.cs +++ /dev/null @@ -1,72 +0,0 @@ -using System; -using FluidSim.Components; - -namespace FluidSim.Core -{ - public static class NozzleFlow - { - public static void Compute(Volume0D vol, double area, double downstreamPressure, - out double massFlow, out double rhoFace, out double uFace, out double pFace, - double gamma = 1.4) - { - massFlow = 0.0; - rhoFace = 0.0; - uFace = 0.0; - pFace = 0.0; - - if (vol == null || vol.Mass <= 0 || vol.Volume <= 0) - return; - - double p0 = vol.Pressure; - double T0 = vol.Temperature; - double R = vol.GasConstant; - double rho0 = vol.Density; - - if (double.IsNaN(p0) || double.IsNaN(T0) || double.IsNaN(rho0) || - p0 <= 0 || T0 <= 0 || rho0 <= 0) - return; - - double pr = downstreamPressure / p0; - double choked = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0)); - - // If pr > 1, flow is INTO the cylinder (reverse), so we swap the roles. - bool reverse = (pr > 1.0); - if (reverse) - { - // Treat the cylinder as the downstream, the pipe as the upstream. - double p_up = downstreamPressure; - double T_up = 300.0; // pipe temperature (ambient) - double rho_up = downstreamPressure / (R * T_up); - - double pr_rev = p0 / p_up; // now cylinder / pipe - if (pr_rev < choked) pr_rev = choked; - - double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr_rev, -(gamma - 1.0) / gamma) - 1.0)); - if (double.IsNaN(M)) return; - - // Flow from pipe INTO cylinder (positive mass flow into volume) - uFace = M * Math.Sqrt(gamma * R * T_up); - rhoFace = rho_up * Math.Pow(pr_rev, 1.0 / gamma); - pFace = p_up * pr_rev; - massFlow = rhoFace * uFace * area; - // massFlow is positive = into cylinder - } - else - { - // Normal flow out of cylinder - if (pr < choked) pr = choked; - - double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0)); - if (double.IsNaN(M)) return; - - uFace = M * Math.Sqrt(gamma * R * T0); - rhoFace = rho0 * Math.Pow(pr, 1.0 / gamma); - pFace = p0 * pr; - massFlow = -rhoFace * uFace * area; // negative = out of cylinder - } - - if (double.IsNaN(massFlow) || double.IsInfinity(massFlow)) - massFlow = 0.0; - } - } -} \ No newline at end of file diff --git a/Core/OpenEndLink.cs b/Core/OpenEndLink.cs new file mode 100644 index 0000000..939c1fe --- /dev/null +++ b/Core/OpenEndLink.cs @@ -0,0 +1,123 @@ +using System; +using FluidSim.Components; + +namespace FluidSim.Core +{ + /// + /// Characteristic open‑end boundary condition. + /// For subsonic outflow the outgoing Riemann invariant is conserved, + /// and the ghost pressure is set to the prescribed ambient value. + /// + 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; + + // Last resolved state (for audio / monitoring) + 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; + } + + /// + /// Compute the ghost state and mass flow for one sub‑step. + /// + public void Resolve(double dtSub) + { + (double rhoInt, double uInt, double pInt) = IsLeftEnd + ? Pipe.GetInteriorStateLeft() + : Pipe.GetInteriorStateRight(); + + double gamma = Gamma; + double gm1 = gamma - 1.0; + double cInt = Math.Sqrt(gamma * pInt / Math.Max(rhoInt, 1e-12)); + double pAmb = AmbientPressure; + + double rhoGhost, uGhost, pGhost; + double mdot; + + if (IsLeftEnd) + { + // Left end: outgoing invariant is J- = u - 2c/(γ-1) + double J_minus = uInt - 2.0 * cInt / gm1; + + if (uInt <= -cInt) // supersonic inflow (all info from outside) + { + // Simple reservoir model – use ambient density and temperature 300 K + rhoGhost = pAmb / (287.0 * 300.0); + uGhost = uInt; // keep interior velocity (should be supersonic inward) + pGhost = pAmb; + } + else if (uInt < 0) // subsonic inflow + { + double rhoAmb = pAmb / (287.0 * 300.0); + double cAmb = Math.Sqrt(gamma * pAmb / rhoAmb); + uGhost = J_minus + 2.0 * cAmb / gm1; + rhoGhost = rhoAmb; + pGhost = pAmb; + } + else // subsonic outflow (uInt >= 0) + { + double s = pInt / Math.Pow(rhoInt, gamma); + rhoGhost = Math.Pow(pAmb / s, 1.0 / gamma); + double cGhost = Math.Sqrt(gamma * pAmb / rhoGhost); + uGhost = J_minus + 2.0 * cGhost / gm1; + if (uGhost < 0) uGhost = 0; + pGhost = pAmb; + } + } + else // Right end + { + // Right end: outgoing invariant is J+ = u + 2c/(γ-1) + double J_plus = uInt + 2.0 * cInt / gm1; + + if (uInt >= cInt) // supersonic outflow + { + rhoGhost = rhoInt; + uGhost = uInt; + pGhost = pInt; + } + else if (uInt >= 0) // subsonic outflow + { + double s = pInt / Math.Pow(rhoInt, gamma); + rhoGhost = Math.Pow(pAmb / s, 1.0 / gamma); + double cGhost = Math.Sqrt(gamma * pAmb / rhoGhost); + uGhost = J_plus - 2.0 * cGhost / gm1; + if (uGhost < 0) uGhost = 0; + pGhost = pAmb; + } + else // subsonic inflow (uInt < 0) + { + double rhoAmb = pAmb / (287.0 * 300.0); + double cAmb = Math.Sqrt(gamma * pAmb / rhoAmb); + uGhost = J_plus - 2.0 * cAmb / gm1; + rhoGhost = rhoAmb; + pGhost = pAmb; + } + } + + // Apply ghost to pipe + if (IsLeftEnd) + Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost); + else + Pipe.SetGhostRight(rhoGhost, uGhost, pGhost); + + // Mass flow (positive = out of pipe) + double area = Pipe.Area; + mdot = rhoGhost * uGhost * area; + if (IsLeftEnd) mdot = -mdot; // positive u into pipe, so out of pipe is negative u + LastMassFlowRate = mdot; + LastFaceDensity = rhoGhost; + LastFaceVelocity = uGhost; + LastFacePressure = pGhost; + } + } +} \ No newline at end of file diff --git a/Core/OrificeBoundary.cs b/Core/OrificeBoundary.cs deleted file mode 100644 index fcf8785..0000000 --- a/Core/OrificeBoundary.cs +++ /dev/null @@ -1,100 +0,0 @@ -using System; -using FluidSim.Interfaces; - -namespace FluidSim.Core -{ - public static class OrificeBoundary - { - public static double MassFlow(double pA, double rhoA, double pB, double rhoB, - Connection conn) - { - if (double.IsNaN(pA) || double.IsNaN(rhoA) || double.IsNaN(pB) || double.IsNaN(rhoB) || - double.IsInfinity(pA) || double.IsInfinity(rhoA) || double.IsInfinity(pB) || double.IsInfinity(rhoB) || - pA <= 0 || rhoA <= 0 || pB <= 0 || rhoB <= 0) - return 0.0; - - double dp = pA - pB; - double sign = Math.Sign(dp); - double absDp = Math.Abs(dp); - double rhoUp = dp >= 0 ? rhoA : rhoB; - double pUp = dp >= 0 ? pA : pB; - double pDown = dp >= 0 ? pB : pA; - double delta = 1e-6 * pUp; - - if (absDp < delta) - { - double k = conn.DischargeCoefficient * conn.Area * Math.Sqrt(2 * rhoUp / delta); - return k * dp; - } - else - { - double pr = pDown / pUp; - double choked = Math.Pow(2.0 / (conn.Gamma + 1.0), conn.Gamma / (conn.Gamma - 1.0)); - if (pr < choked) - { - double term = Math.Sqrt(conn.Gamma * - Math.Pow(2.0 / (conn.Gamma + 1.0), (conn.Gamma + 1.0) / (conn.Gamma - 1.0))); - double flow = conn.DischargeCoefficient * conn.Area * - Math.Sqrt(rhoUp * pUp) * term; - return sign * flow; - } - else - { - double ex = 1.0 - Math.Pow(pr, (conn.Gamma - 1.0) / conn.Gamma); - double flow = conn.DischargeCoefficient * conn.Area * - Math.Sqrt(2.0 * rhoUp * pUp * (conn.Gamma / (conn.Gamma - 1.0)) * - pr * pr * ex); - return sign * flow; - } - } - } - - public static void PipeVolumeFlux(double pPipe, double rhoPipe, double uPipe, - double pVol, double rhoVol, double uVol, - Connection conn, double pipeArea, - bool isLeftBoundary, - out double massFlux, out double momFlux, out double energyFlux) - { - // ----- Compute STAGNATION pressures ----- - double pStagPipe = pPipe + 0.5 * rhoPipe * uPipe * uPipe; - double pStagVol = pVol + 0.5 * rhoVol * uVol * uVol; // uVol is always 0 for your volumes - - // Mass flow driven by stagnation pressure difference (positive = pipe→volume) - double mdot = MassFlow(pStagPipe, rhoPipe, pStagVol, rhoVol, conn); - - // Limit mass flow to the amount that can leave/enter the pipe cell - double maxMdot = rhoPipe * pipeArea * 343.0; - if (Math.Abs(mdot) > maxMdot) mdot = Math.Sign(mdot) * maxMdot; - - bool flowLeavesPipe = mdot > 0; // pipe → volume - - double uFace, pFace, rhoFace; - double massFluxPerArea; - - if (isLeftBoundary) - { - massFluxPerArea = -mdot / pipeArea; - if (flowLeavesPipe) - { uFace = uPipe; pFace = pPipe; rhoFace = rhoPipe; } - else - { uFace = uVol; pFace = pVol; rhoFace = rhoVol; } - } - else // right boundary - { - massFluxPerArea = mdot / pipeArea; - if (flowLeavesPipe) - { uFace = uPipe; pFace = pPipe; rhoFace = rhoPipe; } - else - { uFace = uVol; pFace = pVol; rhoFace = rhoVol; } - } - - // Total enthalpy of the injected fluid - double specificEnthalpy = (1.4 / (1.4 - 1.0)) * pFace / Math.Max(rhoFace, 1e-12); - double totalEnthalpy = specificEnthalpy + 0.5 * uFace * uFace; - - massFlux = massFluxPerArea; - momFlux = massFluxPerArea * uFace + pFace; - energyFlux = massFluxPerArea * totalEnthalpy; - } - } -} \ No newline at end of file diff --git a/Core/OrificeLink.cs b/Core/OrificeLink.cs new file mode 100644 index 0000000..156590b --- /dev/null +++ b/Core/OrificeLink.cs @@ -0,0 +1,140 @@ +using System; +using FluidSim.Components; +using FluidSim.Interfaces; + +namespace FluidSim.Core +{ + /// + /// Connects a port (volume or atmosphere) to one end of a pipe via an orifice. + /// The area can be dynamic (Func). + /// + 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 Gamma { get; set; } = 1.4; + public double GasConstant { get; set; } = 287.0; + + // Last resolved state (for audio/monitoring) + 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 OrificeLink(Port volumePort, Pipe1D pipe, bool isPipeLeftEnd, Func areaProvider) + { + VolumePort = volumePort ?? throw new ArgumentNullException(nameof(volumePort)); + Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe)); + IsPipeLeftEnd = isPipeLeftEnd; + AreaProvider = areaProvider ?? throw new ArgumentNullException(nameof(areaProvider)); + } + + /// + /// Resolve the coupling for one sub‑step. Computes nozzle flow (isentropic) + /// and sets the pipe ghost cell and the port flow rates. + /// + public void Resolve(double dtSub) + { + double area = AreaProvider(); + if (area < 1e-12) + { + SetClosedWall(); + return; + } + + // Retrieve volume state + double volP = VolumePort.Pressure; + double volRho = VolumePort.Density; + double volT = VolumePort.Temperature; + double volH = VolumePort.SpecificEnthalpy; + + // Retrieve pipe interior state at the connected end + (double pipeRho, double pipeU, double pipeP) = IsPipeLeftEnd + ? Pipe.GetInteriorStateLeft() + : Pipe.GetInteriorStateRight(); + + // Determine upstream/downstream: if volume pressure > pipe pressure, flow is out of volume (negative into volume). + bool flowOutOfVolume = volP > pipeP; + double pUp, rhoUp, TUp, pDown; + if (flowOutOfVolume) + { + pUp = volP; rhoUp = volRho; TUp = volT; pDown = pipeP; + } + else + { + // Pipe is upstream + pUp = pipeP; rhoUp = pipeRho; TUp = pipeP / (pipeRho * GasConstant); // temperature from pipe + pDown = volP; + } + + // Compute isentropic nozzle flow + IsentropicOrifice.Compute(pUp, rhoUp, TUp, Gamma, GasConstant, pDown, area, DischargeCoefficient, + out double mdotUpstreamToDown, out double rhoFace, out double uFace, out double pFace); + + // mdotUpstreamToDown is positive from upstream to downstream. + // Convert to mass flow into volume (positive mdot = into volume). + double mdotVolume; + if (flowOutOfVolume) + mdotVolume = -mdotUpstreamToDown; // out of volume is negative + else + mdotVolume = mdotUpstreamToDown; // into volume is positive + + // Clamp mass flow to available mass in volume (if it is a Volume0D) + if (VolumePort.Owner is Volume0D vol) + { + double maxMdot = vol.Mass / dtSub; + if (mdotVolume > maxMdot) mdotVolume = maxMdot; + if (mdotVolume < -maxMdot) mdotVolume = -maxMdot; + } + + // Apply ghost state to pipe + if (IsPipeLeftEnd) + Pipe.SetGhostLeft(rhoFace, uFace, pFace); + else + Pipe.SetGhostRight(rhoFace, uFace, pFace); + + // Store results + LastMassFlowRate = mdotVolume; + LastFaceDensity = rhoFace; + LastFaceVelocity = uFace; + LastFacePressure = pFace; + + // Set port flow rates for volume integration + VolumePort.MassFlowRate = mdotVolume; + if (mdotVolume >= 0) + { + // Inflow: enthalpy comes from upstream (pipe) + double pPipe = pipeP; + double rhoPipe = pipeRho; + VolumePort.SpecificEnthalpy = Gamma / (Gamma - 1.0) * pPipe / rhoPipe; + } + else + { + // Outflow: volume's own specific enthalpy + VolumePort.SpecificEnthalpy = volH; + } + } + + private void SetClosedWall() + { + var (rInt, uInt, pInt) = IsPipeLeftEnd + ? Pipe.GetInteriorStateLeft() + : Pipe.GetInteriorStateRight(); + + if (IsPipeLeftEnd) + Pipe.SetGhostLeft(rInt, -uInt, pInt); + else + Pipe.SetGhostRight(rInt, -uInt, pInt); + + LastMassFlowRate = 0.0; + LastFaceDensity = rInt; + LastFaceVelocity = 0.0; + LastFacePressure = pInt; + VolumePort.MassFlowRate = 0.0; + // Keep specific enthalpy as is (not used) + } + } +} \ No newline at end of file diff --git a/Core/PipeVolumeConnection.cs b/Core/PipeVolumeConnection.cs deleted file mode 100644 index 4ff4b74..0000000 --- a/Core/PipeVolumeConnection.cs +++ /dev/null @@ -1,22 +0,0 @@ -using FluidSim.Components; - -namespace FluidSim.Core -{ - public class PipeVolumeConnection - { - public Volume0D Volume { get; } - public Pipe1D Pipe { get; } - public bool IsPipeLeftEnd { get; } - public double OrificeArea { get; set; } - - public double LastMassFlowIntoVolume { get; set; } - - public PipeVolumeConnection(Volume0D vol, Pipe1D pipe, bool isPipeLeftEnd, double orificeArea) - { - Volume = vol; - Pipe = pipe; - IsPipeLeftEnd = isPipeLeftEnd; - OrificeArea = orificeArea; - } - } -} \ No newline at end of file diff --git a/Core/Solver.cs b/Core/Solver.cs index 66c019d..c3d348b 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -1,120 +1,84 @@ using System; using System.Collections.Generic; +using System.Linq; using FluidSim.Components; +using FluidSim.Interfaces; namespace FluidSim.Core { + /// + /// Top‑level solver that owns all components and couplings, + /// orchestrates sub‑stepping, and exposes states for audio. + /// public class Solver { - private readonly List _volumes = new(); - private readonly List _pipes = new(); - private readonly List _connections = new(); + private readonly List _components = new(); + private readonly List _orificeLinks = new(); + private readonly List _junctions = new(); + private readonly List _openEndLinks = new(); private double _dt; - private double _ambientPressure = 101325.0; - public void SetAmbientPressure(double p) => _ambientPressure = p; - public void AddVolume(Volume0D v) => _volumes.Add(v); - public void AddPipe(Pipe1D p) => _pipes.Add(p); - public void AddConnection(PipeVolumeConnection c) => _connections.Add(c); public void SetTimeStep(double dt) => _dt = dt; - public void SetPipeBoundary(Pipe1D pipe, bool isA, BoundaryType type, double ambientPressure = 101325.0) + public void AddComponent(IComponent component) => _components.Add(component); + public void AddOrificeLink(OrificeLink link) => _orificeLinks.Add(link); + public void AddJunction(Junction junction) => _junctions.Add(junction); + public void AddOpenEndLink(OpenEndLink link) => _openEndLinks.Add(link); + + // Convenience: first pipe’s port B mass flow (often the exhaust) + public double ExhaustMassFlow { - if (isA) + get { - pipe.SetABoundaryType(type); - if (type == BoundaryType.OpenEnd) pipe.SetAAmbientPressure(ambientPressure); - } - else - { - pipe.SetBBoundaryType(type); - if (type == BoundaryType.OpenEnd) pipe.SetBAmbientPressure(ambientPressure); + var pipes = _components.OfType().ToList(); + if (pipes.Count > 0) + return Math.Abs(pipes[0].PortB.MassFlowRate); + return 0.0; } } - public float Step() + /// + /// Advance the whole system by one global time step. + /// + public void Step() { - // 1. For each connection, handle flow or closed wall - foreach (var conn in _connections) - { - double area = conn.OrificeArea; - if (area < 1e-12) // valve closed → treat as solid wall - { - conn.Volume.MassFlowRateIn = 0.0; - conn.Volume.SpecificEnthalpyIn = conn.Volume.SpecificEnthalpy; // not used + var pipes = _components.OfType().ToList(); + if (pipes.Count == 0) return; - // Set ghost to a reflective wall (u = -u_pipe, same p, ρ) - int cellIdx = conn.IsPipeLeftEnd ? 0 : conn.Pipe.GetCellCount() - 1; - double rho = Math.Max(conn.Pipe.GetCellDensity(cellIdx), 1e-6); - double p = Math.Max(conn.Pipe.GetCellPressure(cellIdx), 100.0); - double u = conn.Pipe.GetCellVelocity(cellIdx); - if (conn.IsPipeLeftEnd) - conn.Pipe.SetGhostLeft(rho, -u, p); - else - conn.Pipe.SetGhostRight(rho, -u, p); - continue; - } - - // Valve open → use the nozzle model - double downstreamPressure = conn.IsPipeLeftEnd - ? conn.Pipe.GetCellPressure(0) - : conn.Pipe.GetCellPressure(conn.Pipe.GetCellCount() - 1); - - NozzleFlow.Compute(conn.Volume, area, downstreamPressure, - out double mdot, out double rhoFace, out double uFace, out double pFace, - gamma: conn.Volume.Gamma); - - // Clamp mdot to available mass - double maxMdot = conn.Volume.Mass / _dt; - conn.LastMassFlowIntoVolume = mdot; - if (mdot > maxMdot) mdot = maxMdot; - if (mdot < -maxMdot) mdot = -maxMdot; - - conn.Volume.MassFlowRateIn = mdot; - // enthalpy: if inflow, use pipe enthalpy; if outflow, use cylinder enthalpy - if (mdot >= 0) - { - int cellIdx = conn.IsPipeLeftEnd ? 0 : conn.Pipe.GetCellCount() - 1; - double pPipe = Math.Max(conn.Pipe.GetCellPressure(cellIdx), 100.0); - double rhoPipe = Math.Max(conn.Pipe.GetCellDensity(cellIdx), 1e-6); - conn.Volume.SpecificEnthalpyIn = (conn.Volume.Gamma / (conn.Volume.Gamma - 1.0)) * pPipe / rhoPipe; - } - else - { - conn.Volume.SpecificEnthalpyIn = conn.Volume.SpecificEnthalpy; - } - - // Integrate the volume - conn.Volume.Integrate(_dt); - - // Set ghost from nozzle face state (but don't allow zero density) - if (rhoFace < 1e-6) rhoFace = Constants.Rho_amb; - if (pFace < 100.0) pFace = Constants.P_amb; - if (conn.IsPipeLeftEnd) - conn.Pipe.SetGhostLeft(rhoFace, uFace, pFace); - else - conn.Pipe.SetGhostRight(rhoFace, uFace, pFace); - } - - // 2. Sub‑step pipes + // 1. Determine sub‑step count (max CFL over all pipes) int nSub = 1; - foreach (var p in _pipes) + foreach (var p in pipes) nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt)); double dtSub = _dt / nSub; + // 2. Sub‑step loop for (int sub = 0; sub < nSub; sub++) - foreach (var p in _pipes) + { + // a) Resolve all orifice links (volume ↔ pipe) + foreach (var link in _orificeLinks) + link.Resolve(dtSub); + + // b) Resolve all open‑end links (pipe → atmosphere) + foreach (var link in _openEndLinks) + link.Resolve(dtSub); + + // c) Resolve all junctions (pipe ↔ pipe) + foreach (var junc in _junctions) + junc.Resolve(dtSub); + + // d) Advance all pipes + foreach (var p in pipes) p.SimulateSingleStep(dtSub); + } // 3. Clear ghost flags - foreach (var p in _pipes) - p.ClearGhostFlag(); + foreach (var p in pipes) + p.ClearGhostFlags(); - // 4. Return exhaust tailpipe mass flow - if (_pipes.Count > 0) - return (float)_pipes[0].GetOpenEndMassFlow(); - return 0f; + // 4. Integrate non‑pipe components (volumes, atmosphere, etc.) + foreach (var comp in _components) + comp.UpdateState(_dt); } } } \ No newline at end of file diff --git a/Interfaces/Connection.cs b/Interfaces/Connection.cs deleted file mode 100644 index 409542c..0000000 --- a/Interfaces/Connection.cs +++ /dev/null @@ -1,15 +0,0 @@ -namespace FluidSim.Interfaces -{ - /// Pure data link between two ports, with orifice parameters. - public class Connection - { - public Port PortA { get; } - public Port PortB { get; } - - public double Area { get; set; } = 1e-5; // effective orifice area (m²) - public double DischargeCoefficient { get; set; } = 0.62; - public double Gamma { get; set; } = 1.4; - - public Connection(Port a, Port b) => (PortA, PortB) = (a, b); - } -} \ No newline at end of file diff --git a/Interfaces/IComponent.cs b/Interfaces/IComponent.cs new file mode 100644 index 0000000..b23be76 --- /dev/null +++ b/Interfaces/IComponent.cs @@ -0,0 +1,19 @@ +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); + } +} \ No newline at end of file diff --git a/Interfaces/Port.cs b/Interfaces/Port.cs index 10cedbb..0b5a488 100644 --- a/Interfaces/Port.cs +++ b/Interfaces/Port.cs @@ -1,18 +1,28 @@ namespace FluidSim.Interfaces { + /// + /// A fluid port that belongs to a component. The solver writes mass flow, + /// enthalpy, etc. here, and reads pressure, density, etc. + /// public class Port { + // Set by the solver during coupling resolution + public double MassFlowRate; // kg/s, positive INTO the component that owns this port + public double SpecificEnthalpy; // J/kg, enthalpy of the fluid crossing the port (inflow direction) + + // Set by the owning component after integration to reflect its current state public double Pressure; // Pa - public double MassFlowRate; // kg/s, positive INTO the component - public double SpecificEnthalpy; // J/kg, enthalpy of fluid entering this port public double Density; // kg/m³ public double Temperature; // K + // Link back to owner (optional, but handy for debugging) + public object? Owner { get; set; } + public Port() { - Pressure = 101325.0; MassFlowRate = 0.0; SpecificEnthalpy = 0.0; + Pressure = 101325.0; Density = 1.225; Temperature = 300.0; } diff --git a/Interfaces/SoundConnection.cs b/Interfaces/SoundConnection.cs deleted file mode 100644 index 304ce7a..0000000 --- a/Interfaces/SoundConnection.cs +++ /dev/null @@ -1,25 +0,0 @@ -namespace FluidSim.Interfaces -{ - /// - /// A Connection that also produces an audio sample from the pressure drop across it. - /// - public class SoundConnection : Connection - { - /// Gain applied to the normalised pressure difference. - public float Gain { get; set; } = 1.0f; - - /// Reference pressure used for normalisation (Pa). Default: 1 atm. - public double ReferencePressure { get; set; } = 101325.0; - - public SoundConnection(Port a, Port b) : base(a, b) { } - - /// - /// Returns a normalised audio sample proportional to the pressure difference. - /// - public float GetAudioSample() - { - double dp = PortA.Pressure - PortB.Pressure; - return (float)(dp / ReferencePressure) * Gain; - } - } -} \ No newline at end of file diff --git a/Program.cs b/Program.cs index 94ffa98..330b59b 100644 --- a/Program.cs +++ b/Program.cs @@ -3,6 +3,7 @@ using SFML.Window; using SFML.System; using System.Diagnostics; using FluidSim.Core; +using FluidSim.Tests; namespace FluidSim; @@ -42,7 +43,7 @@ public class Program soundEngine.Volume = 100; soundEngine.Start(); - scenario = new EngineScenario(); + scenario = new TestScenario(); scenario.Initialize(SampleRate); var stopwatch = Stopwatch.StartNew(); @@ -72,12 +73,6 @@ public class Program double speedSmoothing = 8.0; currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-speedSmoothing * dtClock)); - // ---- THROTTLE INPUT ---- - targetThrottle = Keyboard.IsKeyPressed(Keyboard.Key.W) ? 1.0 : 0.0; - currentThrottle += (targetThrottle - currentThrottle) * (1.0 - Math.Exp(-ThrottleSmoothing * dtClock)); - // Push to engine scenario (if it's an EngineScenario) - if (scenario is EngineScenario engine) - engine.Throttle = currentThrottle; // Generate audio double targetAudioClock = currentRealTime + 0.05; diff --git a/Scenarios/EngineScenario.cs b/Scenarios/EngineScenario.cs deleted file mode 100644 index 57f43d3..0000000 --- a/Scenarios/EngineScenario.cs +++ /dev/null @@ -1,325 +0,0 @@ -using System; -using FluidSim.Components; -using FluidSim.Utils; -using FluidSim.Interfaces; -using SFML.Graphics; -using SFML.System; - -namespace FluidSim.Core -{ - public class EngineScenario : Scenario - { - private Solver solver; - private Crankshaft crankshaft; - private EngineCylinder engineCyl; - private Pipe1D exhaustPipe; - private Pipe1D intakePipe; - private PipeVolumeConnection couplingExhaust; - private PipeVolumeConnection couplingIntake; - private SoundProcessor exhaustSoundProcessor; - private SoundProcessor intakeSoundProcessor; - private OutdoorExhaustReverb reverb; - - private Port exhaustPort = new Port(); - private Port intakePort = new Port(); - - private double dt; - private double exhPipeArea, intPipeArea; - private const double AmbientPressure = 101325.0; - private double time; - private int stepCount = 0; - private const int LogInterval = 1000; - - public double Throttle { get; set; } = 0.15; - private const double FullLoadPeakPressure = 140.0 * Units.bar; - - public override void Initialize(int sampleRate) - { - dt = 1.0 / sampleRate; - - // Crankshaft - crankshaft = new Crankshaft(initialRPM: 2000.0) - { - Inertia = 0.05, - FrictionConstant = 0.2, - FrictionViscous = 0.025 - }; - - // Exhaust pipe (longer, larger) - double exhLength = 0.5; - double exhRadius = 1.5 * Units.cm; - exhPipeArea = Math.PI * exhRadius * exhRadius; - exhaustPipe = new Pipe1D(exhLength, exhPipeArea, sampleRate, forcedCellCount: 30); - exhaustPipe.SetUniformState(1.225, 0.0, AmbientPressure); - exhaustPipe.DampingMultiplier = 0.0; - exhaustPipe.EnergyRelaxationRate = 100.0f; - - // Intake pipe (shorter, narrower) - double intLength = 0.1; - double intRadius = 1 * Units.cm; - intPipeArea = Math.PI * intRadius * intRadius; - intakePipe = new Pipe1D(intLength, intPipeArea, sampleRate, forcedCellCount: 10); - intakePipe.SetUniformState(1.225, 0.0, AmbientPressure); - - // Cylinder (starts at BDC, fresh charge) - engineCyl = new EngineCylinder(crankshaft, - bore: 56 * Units.mm, stroke: 57 * Units.mm, compressionRatio: 9.5, - exhPipeArea: exhPipeArea, intPipeArea: intPipeArea, sampleRate: sampleRate); - engineCyl.ignition = true; - - // Set crank to BDC (180°) and sync - crankshaft.CrankAngle = Math.PI; - crankshaft.PreviousAngle = Math.PI; // make sure this property is settable (public setter) - - // Couplings - couplingExhaust = new PipeVolumeConnection(engineCyl.Cylinder, exhaustPipe, true, orificeArea: 0.0); - couplingIntake = new PipeVolumeConnection(engineCyl.Cylinder, intakePipe, false, orificeArea: 0.0); - - // Solver - solver = new Solver(); - solver.SetTimeStep(dt); - solver.AddVolume(engineCyl.Cylinder); - solver.AddPipe(exhaustPipe); - solver.AddPipe(intakePipe); - solver.AddConnection(couplingExhaust); - solver.AddConnection(couplingIntake); - solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.OpenEnd, AmbientPressure); - solver.SetPipeBoundary(intakePipe, true, BoundaryType.GhostCell); // cylinder side – left - solver.SetPipeBoundary(intakePipe, false, BoundaryType.OpenEnd, AmbientPressure); // ambient side – right - - // Sound - exhaustSoundProcessor = new SoundProcessor(sampleRate, exhRadius * 2); - exhaustSoundProcessor.Gain = 0.001f; - - intakeSoundProcessor = new SoundProcessor(sampleRate, intRadius * 2); - intakeSoundProcessor.Gain = 0.001f; - - // Reverb - reverb = new OutdoorExhaustReverb(sampleRate); - reverb.DryMix = 1.0f; - reverb.EarlyMix = 0.5f; - reverb.TailMix = 0.9f; - reverb.Feedback = 0.9f; - reverb.DampingFreq = 6000f; - - Console.WriteLine("=== Engine with intake & cycle‑aware valves ==="); - } - - public override float Process() - { - double idleThrottle = 0.1; - if (crankshaft.AngularVelocity < 80) idleThrottle = 0.2; - double throttle = Math.Clamp(Throttle, idleThrottle, 1.0); - double targetPressure = throttle * FullLoadPeakPressure; - engineCyl.TargetPeakPressure = targetPressure; - - engineCyl.Step(dt); - crankshaft.Step(dt); - - couplingExhaust.OrificeArea = engineCyl.ExhaustOrificeArea; - couplingIntake.OrificeArea = engineCyl.IntakeOrificeArea; - - solver.Step(); - - UpdateExhaustPort(); - UpdateIntakePort(); - float dryExhaust = exhaustSoundProcessor.Process(exhaustPort); - float dryIntake = intakeSoundProcessor.Process(intakePort); - float dry = dryExhaust + dryIntake; - - float wet = reverb.Process(dry); - - if (++stepCount % LogInterval == 0) Log(); - - return wet; - } - - private void Log() - { - double rpm = crankshaft.AngularVelocity * 60.0 / (2.0 * Math.PI); - double cycleDeg = (engineCyl.CycleAngle * 180.0 / Math.PI) % 720.0; - string stroke = cycleDeg < 180.0 ? "Power" : - cycleDeg < 360.0 ? "Exhaust" : - cycleDeg < 540.0 ? "Intake" : "Compression"; - - // Cylinder - double pCyl = engineCyl.Cylinder.Pressure; - double TCyl = engineCyl.Cylinder.Temperature; - double VCyl = engineCyl.Cylinder.Volume; - double mCyl = engineCyl.Cylinder.Mass; - double exhArea = engineCyl.ExhaustOrificeArea * 1e6; // mm² - double intArea = engineCyl.IntakeOrificeArea * 1e6; // mm² - - // Exhaust pipe - int exhLast = exhaustPipe.GetCellCount() - 1; - double pExhEnd = exhaustPipe.GetCellPressure(exhLast); - double mdotExhOut = exhaustPipe.GetOpenEndMassFlow(); // positive out - - // Intake pipe - double mdotIntIn = couplingIntake.LastMassFlowIntoVolume; - double pIntAmbEnd = intakePort.Pressure; - - Console.WriteLine( - $"{stepCount,8} {stroke,-11} {cycleDeg,6:F1}° " + - $"RPM:{rpm,5:F0} " + - $"Cyl: p={pCyl/1e5,6:F3}bar T={TCyl,6:F0}K V={VCyl*1e6,6:F0}cm³ m={mCyl*1e3,6:F6}g " + - $"Valves: Exh={exhArea,5:F0}mm² Int={intArea,5:F0}mm² " + - $"Intake: p_end={pIntAmbEnd/1e5,6:F3}bar mdot_in={mdotIntIn,7:F4}kg/s " + - $"Exhaust: p_end={pExhEnd/1e5,6:F3}bar mdot_out={mdotExhOut,7:F4}kg/s"); - } - - private void UpdateExhaustPort() - { - int last = exhaustPipe.GetCellCount() - 1; - double p = exhaustPipe.GetCellPressure(last); - double rho = exhaustPipe.GetCellDensity(last); - double vel = exhaustPipe.GetCellVelocity(last); - - // Safety clamps - rho = Math.Clamp(rho, 0.01, 50.0); - vel = Math.Clamp(vel, -500.0, 500.0); - p = Math.Clamp(p, 1e4, 2e6); - - double outflowMassFlow = rho * vel * exhPipeArea; - - exhaustPort.Pressure = p; - exhaustPort.Density = rho; - exhaustPort.Temperature = p / (rho * 287.05); - exhaustPort.MassFlowRate = -outflowMassFlow; - exhaustPort.SpecificEnthalpy = 0.0; - } - - private void UpdateIntakePort() - { - // Use the actual valve mass flow (positive = into cylinder) - double mdotIntoEngine = couplingIntake.LastMassFlowIntoVolume; - - // Use cylinder pressure/density for the port state (or intake pipe last cell) - double pCyl = engineCyl.Cylinder.Pressure; - double rhoCyl = engineCyl.Cylinder.Density; - - intakePort.Pressure = Math.Max(pCyl, 100); - intakePort.Density = Math.Max(rhoCyl, 1e-6); - intakePort.Temperature = engineCyl.Cylinder.Temperature; - intakePort.MassFlowRate = mdotIntoEngine; - intakePort.SpecificEnthalpy = 0.0; - } - - // ==================== Drawing ==================== - public override void Draw(RenderWindow target) - { - float winW = target.GetView().Size.X; - float winH = target.GetView().Size.Y; - float centerY = winH / 2f; - - const float T_ambient = 293.15f; - const float T_hot = 1500f; - const float T_cold = 0f; - const float R = 287.05f; - float deltaHot = T_hot - T_ambient; - float deltaCold = T_ambient - T_cold; - - float NormaliseTemperature(double T) - { - double t; - if (T >= T_ambient) - t = (T - T_ambient) / deltaHot; - else - t = (T - T_ambient) / deltaCold; - return (float)Math.Clamp(t, -1.0, 1.0); - } - - // ---- Cylinder ---- - float cylW = 80f, cylH = 150f; - var cylRect = new RectangleShape(new Vector2f(cylW, cylH)); - cylRect.Position = new Vector2f(200f, centerY - cylH / 2f); - double tempCyl = engineCyl.Cylinder.Temperature; - float tnCyl = NormaliseTemperature(tempCyl); - byte rC = (byte)(tnCyl > 0 ? 255 * tnCyl : 0); - byte bC = (byte)(tnCyl < 0 ? -255 * tnCyl : 0); - byte gC = (byte)(255 * (1 - Math.Abs(tnCyl))); - cylRect.FillColor = new Color(rC, gC, bC); - target.Draw(cylRect); - - // ---- Piston ---- - float pistonWidth = cylW - 12f; - float pistonHeight = 16f; - float pistonFraction = (float)engineCyl.PistonPositionFraction; - float pistonTopY = cylRect.Position.Y + pistonFraction * (cylH - pistonHeight); - var pistonRect = new RectangleShape(new Vector2f(pistonWidth, pistonHeight)) - { - Position = new Vector2f(cylRect.Position.X + 6f, pistonTopY), - FillColor = new Color(80, 80, 80) - }; - target.Draw(pistonRect); - - // ---------- NEW: Valve lift indicators ---------- - float barWidth = 30f; - float barHeight = 10f; - float exhLift = (float)engineCyl.ExhaustValveLiftCurrent; - float intLift = (float)engineCyl.IntakeValveLiftCurrent; - - // Exhaust valve indicator (right side of cylinder) - var exhBar = new RectangleShape(new Vector2f(barWidth, barHeight)) - { - Position = new Vector2f(cylRect.Position.X + cylW - 10, - cylRect.Position.Y - 20 - exhLift * 20), - FillColor = new Color(200, 200, 200) - }; - target.Draw(exhBar); - - // Intake valve indicator (left side of cylinder) - var intBar = new RectangleShape(new Vector2f(barWidth, barHeight)) - { - Position = new Vector2f(cylRect.Position.X - 20, - cylRect.Position.Y - 20 - intLift * 20), - FillColor = new Color(200, 200, 200) - }; - target.Draw(intBar); - - // ---- Exhaust pipe (rightwards) ---- - DrawPipe(target, exhaustPipe, startX: 280f, endX: winW - 60f, centerY + 10 - cylRect.Size.Y / 2, - T_ambient, T_hot, T_cold, R, NormaliseTemperature, true); - - // ---- Intake pipe (leftwards) ---- - DrawPipe(target, intakePipe, startX: 200f, endX: 20f, centerY + 10 - cylRect.Size.Y / 2, - T_ambient, T_hot, T_cold, R, NormaliseTemperature, false); - } - - private void DrawPipe(RenderWindow target, Pipe1D pipe, - float startX, float endX, float centerY, - float T_ambient, float T_hot, float T_cold, float R, - Func normaliseTemp, bool leftToRight) - { - int n = pipe.GetCellCount(); - float dir = leftToRight ? 1f : -1f; - float pipeLen = Math.Abs(endX - startX); - float dx = pipeLen / (n - 1) * dir; - float baseRadius = leftToRight ? 20f : 14f; // exhaust thicker, intake thinner - var vertices = new Vertex[n * 2]; - float ambPress = 101325f; - - for (int i = 0; i < n; i++) - { - float x = startX + i * dx; - double p = pipe.GetCellPressure(i); - double rho = pipe.GetCellDensity(i); - double T = p / (rho * R); - - float r = baseRadius * 0.2f * (float)(1.0 + (p - ambPress) / ambPress); - if (r < 2f) r = 2f; - - float tn = normaliseTemp(T); - byte rC = (byte)(tn > 0 ? 255 * tn : 0); - byte bC = (byte)(tn < 0 ? -255 * tn : 0); - byte gC = (byte)(255 * (1 - Math.Abs(tn))); - var col = new Color(rC, gC, bC); - - vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col); - vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col); - } - - target.Draw(vertices, PrimitiveType.TriangleStrip); - } - } -} \ No newline at end of file diff --git a/Scenarios/HelmholtzResonatorScenario.cs b/Scenarios/HelmholtzResonatorScenario.cs deleted file mode 100644 index af8ded5..0000000 --- a/Scenarios/HelmholtzResonatorScenario.cs +++ /dev/null @@ -1,128 +0,0 @@ -using System; -using FluidSim.Components; -using FluidSim.Utils; -using SFML.Graphics; -using SFML.System; - -namespace FluidSim.Core -{ - public class HelmholtzResonatorScenario : Scenario - { - private Solver solver; - private Volume0D cavity; - private Pipe1D neck; - private PipeVolumeConnection coupling; - private int stepCount; - private double time; - private double dt; - private double ambientPressure = 1.0 * Units.atm; - - public override void Initialize(int sampleRate) - { - dt = 1.0 / sampleRate; - - // 1‑litre cavity, 10% over‑pressure - double cavityVolume = 1e-3; - double initialCavityPressure = 1.1 * ambientPressure; - cavity = new Volume0D(cavityVolume, initialCavityPressure, 300.0, sampleRate) - { - Gamma = 1.4, - GasConstant = 287.0 - }; - - // Neck: length 10 cm, radius 1 cm - double neckLength = 0.1; - double neckRadius = 0.01; - double neckArea = Math.PI * neckRadius * neckRadius; - neck = new Pipe1D(neckLength, neckArea, sampleRate, forcedCellCount: 40); - neck.SetUniformState(1.225, 0.0, ambientPressure); - - // Create the coupling between cavity and left end of the neck (PortA) - coupling = new PipeVolumeConnection(cavity, neck, isPipeLeftEnd: true, orificeArea: neckArea); - - solver = new Solver(); - solver.SetTimeStep(dt); - solver.AddVolume(cavity); - solver.AddPipe(neck); - solver.AddConnection(coupling); - - // Left boundary (PortA) is volume‑coupled via ghost cell, right boundary (PortB) is open end - solver.SetPipeBoundary(neck, isA: true, BoundaryType.GhostCell); - solver.SetPipeBoundary(neck, isA: false, BoundaryType.OpenEnd, ambientPressure); - } - - public override float Process() - { - float sample = solver.Step(); - time += dt; - stepCount++; - - double pOpen = neck.GetCellPressure(neck.GetCellCount() - 1); - float audio = (float)((pOpen - ambientPressure) / ambientPressure); - - if (stepCount % 20 == 0) - { - double pCav = cavity.Pressure; - // Mass flow rate is not directly available – we can compute from pressure difference or skip - Console.WriteLine( - $"t={time * 1e3:F2} ms step={stepCount} " + - $"P_cav={pCav:F1} Pa, P_open={pOpen:F1} Pa, " + - $"audio={audio:F4}"); - } - - return audio; - } - - public override void Draw(RenderWindow target) - { - float winW = target.GetView().Size.X; - float winH = target.GetView().Size.Y; - float centerY = winH / 2f; - - // Cavity rectangle - float cavityWidth = 120f; - float cavityHeight = 180f; - var cavityRect = new RectangleShape(new Vector2f(cavityWidth, cavityHeight)); - cavityRect.Position = new Vector2f(40f, centerY - cavityHeight / 2f); - cavityRect.FillColor = PressureColor(cavity.Pressure); - target.Draw(cavityRect); - - // Neck drawn as tapered pipe - int n = neck.GetCellCount(); - float neckStartX = 40f + cavityWidth + 10f; - float neckEndX = winW - 60f; - float neckLenPx = neckEndX - neckStartX; - float dx = neckLenPx / (n - 1); - float baseRadius = 20f; - - var vertices = new Vertex[n * 2]; - for (int i = 0; i < n; i++) - { - float x = neckStartX + i * dx; - double p = neck.GetCellPressure(i); - float r = baseRadius * (float)(0.5 + 0.5 * Math.Tanh((p - ambientPressure) / (ambientPressure * 0.2))); - if (r < 4f) r = 4f; - Color col = PressureColor(p); - vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col); - vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col); - } - target.Draw(vertices, PrimitiveType.TriangleStrip); - - // Open end indicator - var arrow = new CircleShape(8f); - arrow.Position = new Vector2f(neckEndX - 4f, centerY - 4f); - arrow.FillColor = Color.White; - target.Draw(arrow); - } - - private Color PressureColor(double pressure) - { - double range = ambientPressure * 0.1; - double t = Math.Clamp((pressure - ambientPressure) / range, -1.0, 1.0); - byte r = (byte)(t > 0 ? 255 * t : 0); - byte b = (byte)(t < 0 ? -255 * t : 0); - byte g = (byte)(255 * (1 - Math.Abs(t))); - return new Color(r, g, b); - } - } -} \ No newline at end of file diff --git a/Scenarios/PipeResonatorScenario.cs b/Scenarios/PipeResonatorScenario.cs deleted file mode 100644 index e534910..0000000 --- a/Scenarios/PipeResonatorScenario.cs +++ /dev/null @@ -1,184 +0,0 @@ -using FluidSim.Components; -using FluidSim.Interfaces; -using FluidSim.Utils; -using SFML.Graphics; -using SFML.System; -using System; - -namespace FluidSim.Core -{ - public class PipeResonatorScenario : Scenario - { - private Solver solver; - private Pipe1D pipe; - private int stepCount; - private double time; - private double dt; - private double ambientPressure = 1.0 * Units.atm; - private bool enableLogging = true; - - public override void Initialize(int sampleRate) - { - dt = 1.0 / sampleRate; - - double length = 2; - double radius = 50 * Units.mm; - double area = Units.AreaFromDiameter(radius); - - pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: 80); - pipe.SetUniformState(1.225, 0.0, ambientPressure); - - solver = new Solver(); - solver.SetTimeStep(dt); - solver.AddPipe(pipe); - // Open end at port A (left), closed end at port B (right) - solver.SetPipeBoundary(pipe, isA: true, BoundaryType.OpenEnd, ambientPressure); - solver.SetPipeBoundary(pipe, isA: false, BoundaryType.ClosedEnd); - - // Initial pressure pulse - int pulseCells = 5; - double pulsePressure = 2 * ambientPressure; - for (int i = 0; i < pulseCells; i++) - pipe.SetCellState(i, 1.225, 0.0, pulsePressure); - } - - public override float Process() - { - float sample = solver.Step(); - time += dt; - stepCount++; - - double pMid = pipe.GetPressureAtFraction(0.5); - sample = (float)((pMid - ambientPressure) / ambientPressure); - - Log(sample); - return sample; - } - - private void Log(float sample) - { - if (!enableLogging) return; - if (stepCount % 10 == 0 && stepCount < 1000) - { - double pMid = pipe.GetPressureAtFraction(0.5); - double pOpen = pipe.GetCellPressure(0); - double pClosed = pipe.GetCellPressure(pipe.GetCellCount() - 1); - Console.WriteLine( - $"t = {time * 1e3:F3} ms Step {stepCount:D4}: " + - $"sample = {sample:F3}, " + - $"P_mid = {pMid:F2} Pa ({pMid / ambientPressure:F4} atm), " + - $"P_open = {pOpen:F2} Pa, P_closed = {pClosed:F2} Pa"); - } - } - - public override void Draw(RenderWindow target) - { - float winWidth = target.GetView().Size.X; - float winHeight = target.GetView().Size.Y; - - float pipeCenterY = winHeight / 2f; - float margin = 60f; - float pipeStartX = margin; - float pipeEndX = winWidth - margin; - float pipeLengthPx = pipeEndX - pipeStartX; - int n = pipe.GetCellCount(); - float dx = pipeLengthPx / (n - 1); // spacing between cell centres - - float baseRadius = 25f; - float rangeFactor = 1f; - float scaleFactor = 5f; - - // ----- smoothstep helper ----- - 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); - } - - // ----- Pre‑compute cell positions and radii ----- - var centers = new float[n]; - var radii = new float[n]; - for (int i = 0; i < n; i++) - { - double p = pipe.GetCellPressure(i); - float deviation = (float)Math.Tanh((p - ambientPressure) / ambientPressure / rangeFactor); - radii[i] = baseRadius * (1f + deviation * scaleFactor); - if (radii[i] < 2f) radii[i] = 2f; - centers[i] = pipeStartX + i * dx; - } - - // ----- Build triangle‑strip vertices ----- - int segmentsPerCell = 8; // smoothness - int totalPoints = n + (n - 1) * segmentsPerCell; - Vertex[] stripVertices = new Vertex[totalPoints * 2]; // top + bottom for each point - int idx = 0; - - for (int i = 0; i < n; i++) - { - // ---- Cell centre ---- - float x = centers[i]; - float r = radii[i]; - double p = pipe.GetCellPressure(i); - Color col = PressureColor(p); - - stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY - r), col); - stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY + r), col); - - // ---- Intermediate segments after this cell (if not last) ---- - if (i < n - 1) - { - for (int s = 1; s <= segmentsPerCell; s++) - { - float t = s / (float)segmentsPerCell; - float st = SmoothStep(0f, 1f, t); - float xi = centers[i] + (centers[i + 1] - centers[i]) * t; - float ri = radii[i] + (radii[i + 1] - radii[i]) * st; - double pi = pipe.GetCellPressure(i) * (1 - t) + pipe.GetCellPressure(i + 1) * t; - Color coli = PressureColor(pi); - - stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY - ri), coli); - stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY + ri), coli); - } - } - } - - // Draw the pipe as a triangle strip - 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); - - // ----- Closed end indicator (right) ----- - float wallThickness = 8f; - var wall = new RectangleShape(new Vector2f(wallThickness, winHeight * 0.6f)); - wall.Position = new Vector2f(pipeEndX, pipeCenterY - winHeight * 0.6f / 2f); - wall.FillColor = new Color(180, 180, 180); - target.Draw(wall); - } - - /// Blue (low) → Green (ambient) → Red (high). - private Color PressureColor(double pressure) - { - double range = ambientPressure * 0.05; // ±5% gives full colour swing - double t = (pressure - ambientPressure) / range; - t = Math.Clamp(t, -1.0, 1.0); - - byte r, g, b; - if (t < 0) - { - double factor = -t; - r = 0; - g = (byte)(255 * (1 - factor)); - b = (byte)(255 * factor); - } - else - { - double factor = t; - r = (byte)(255 * factor); - g = (byte)(255 * (1 - factor)); - b = 0; - } - return new Color(r, g, b); - } - } -} \ No newline at end of file diff --git a/Scenarios/Scenario.cs b/Scenarios/Scenario.cs index 52fba3f..62483d4 100644 --- a/Scenarios/Scenario.cs +++ b/Scenarios/Scenario.cs @@ -1,23 +1,121 @@ -using SFML.Graphics; +using System; +using SFML.Graphics; +using SFML.System; +using FluidSim.Components; -namespace FluidSim.Core +namespace FluidSim.Tests { public abstract class Scenario { - /// - /// Initialize the scenario with a given audio sample rate. - /// + /// Initialize the scenario with a given audio sample rate. public abstract void Initialize(int sampleRate); - /// - /// Advance one simulation step and return an audio sample. - /// The step size is 1 / sampleRate seconds. - /// + /// Advance one simulation step and return an audio sample. public abstract float Process(); - /// - /// Draw the current simulation state onto the given SFML render target. - /// + /// Draw the current simulation state onto the given SFML render target. public abstract void Draw(RenderWindow target); + + // ---------- Shared drawing helpers ---------- + + protected const double AmbientPressure = 101325.0; + + /// Blue (low) → Green (ambient) → Red (high). + protected Color PressureColor(double pressure) + { + double range = AmbientPressure * 0.05; // ±5% gives full colour swing + double t = (pressure - AmbientPressure) / range; + t = Math.Clamp(t, -1.0, 1.0); + + byte r, g, b; + if (t < 0) + { + double factor = -t; + r = 0; + g = (byte)(255 * (1 - factor)); + b = (byte)(255 * factor); + } + else + { + double factor = t; + r = (byte)(255 * factor); + g = (byte)(255 * (1 - factor)); + b = 0; + } + return new Color(r, g, b); + } + + /// + /// Draws the pipe as a smooth triangle‑strip whose radius varies with cell pressure. + /// + protected void DrawPipe(RenderWindow target, Pipe1D pipe, float pipeCenterY, float pipeStartX, float pipeEndX) + { + int n = pipe.CellCount; + if (n < 2) return; + + float pipeLengthPx = pipeEndX - pipeStartX; + float dx = pipeLengthPx / (n - 1); // spacing between cell centres + + float baseRadius = 25f; + float rangeFactor = 1f; + float scaleFactor = 5f; + + // ----- smoothstep helper ----- + 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); + } + + // ----- Pre‑compute cell positions and radii ----- + var centers = new float[n]; + var radii = new float[n]; + for (int i = 0; i < n; i++) + { + double p = pipe.GetCellPressure(i); + float deviation = (float)Math.Tanh((p - AmbientPressure) / AmbientPressure / rangeFactor); + radii[i] = baseRadius * (1f + deviation * scaleFactor); + if (radii[i] < 2f) radii[i] = 2f; + centers[i] = pipeStartX + i * dx; + } + + // ----- Build triangle‑strip vertices ----- + 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 x = centers[i]; + float r = radii[i]; + double p = pipe.GetCellPressure(i); + Color col = PressureColor(p); + + stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY - r), col); + stripVertices[idx++] = new Vertex(new Vector2f(x, pipeCenterY + r), col); + + if (i < n - 1) + { + for (int s = 1; s <= segmentsPerCell; s++) + { + float t = s / (float)segmentsPerCell; + float st = SmoothStep(0f, 1f, t); + float xi = centers[i] + (centers[i + 1] - centers[i]) * t; + float ri = radii[i] + (radii[i + 1] - radii[i]) * st; + double pi = pipe.GetCellPressure(i) * (1 - t) + pipe.GetCellPressure(i + 1) * t; + Color coli = PressureColor(pi); + + stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY - ri), coli); + stripVertices[idx++] = new Vertex(new Vector2f(xi, pipeCenterY + ri), 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/SodShockTubeScenario.cs b/Scenarios/SodShockTubeScenario.cs deleted file mode 100644 index ff06758..0000000 --- a/Scenarios/SodShockTubeScenario.cs +++ /dev/null @@ -1,158 +0,0 @@ -using System; -using FluidSim.Components; -using FluidSim.Utils; -using SFML.Graphics; -using SFML.System; - -namespace FluidSim.Core -{ - public class SodShockTubeScenario : Scenario - { - private Solver solver; - private Pipe1D pipe; - private int stepCount; - private double time; - private double dt; - private double ambientPressure = 1.0 * Units.atm; - private const double GasConstant = 287.0; - - public override void Initialize(int sampleRate) - { - dt = 1.0 / sampleRate; - double length = 1.0; - double area = 1.0; - int nCells = 200; - - pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: nCells); - pipe.SetUniformState(0.125, 0.0, 0.1 * ambientPressure); // right state - - // Left half high pressure - for (int i = 0; i < nCells / 2; i++) - pipe.SetCellState(i, 1.0, 0.0, ambientPressure); - - solver = new Solver(); - solver.SetTimeStep(dt); - solver.AddPipe(pipe); - solver.SetPipeBoundary(pipe, isA: true, BoundaryType.ClosedEnd); - solver.SetPipeBoundary(pipe, isA: false, BoundaryType.ClosedEnd); - } - - public override float Process() - { - float sample = solver.Step(); - time += dt; - stepCount++; - - double pMid = pipe.GetPressureAtFraction(0.5); - float audio = (float)((pMid - ambientPressure) / ambientPressure); - - bool log = true; - - if (log) - { - int n = pipe.GetCellCount(); - Console.WriteLine($"step {stepCount}:"); - Console.WriteLine("i rho (kg/m³) p (Pa) T (K) u (m/s)"); - for (int i = 0; i < n; i++) - { - if (i % 10 == 0) - { - double rho = pipe.GetCellDensity(i); - double p = pipe.GetCellPressure(i); - double u = pipe.GetCellVelocity(i); - double T = p / (rho * GasConstant); // GasConstant = 287.0 - Console.WriteLine($"{i,-4} {rho,10:F4} {p,10:F1} {T,8:F2} {u,10:F4}"); - } - } - Console.WriteLine(); - } - - return audio; - } - - public override void Draw(RenderWindow target) - { - float winW = target.GetView().Size.X; - float winH = target.GetView().Size.Y; - float centerY = winH / 2f; - float margin = 40f; - float pipeStartX = margin; - float pipeEndX = winW - margin; - float pipeLenPx = pipeEndX - pipeStartX; - int n = pipe.GetCellCount(); - float dx = pipeLenPx / (n - 1); - float baseRadius = 60f; - - Vertex[] vertices = new Vertex[n * 2]; - for (int i = 0; i < n; i++) - { - float x = pipeStartX + i * dx; - - double p = pipe.GetCellPressure(i); - double rho = pipe.GetCellDensity(i); - double T = p / (rho * GasConstant); // temperature in Kelvin - - // Radius from pressure (exaggerated deviation) - float r = baseRadius * (float)(p / ambientPressure * 2); - if (r < 4f) r = 4f; - - // Colour from temperature - Color col = TemperatureColor(T); - - vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col); - vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col); - } - target.Draw(vertices, PrimitiveType.TriangleStrip); - - // Diaphragm marker (faint white line at initial interface) - float diaphragmX = pipeStartX + (n / 2) * dx; - var line = new RectangleShape(new Vector2f(2f, winH * 0.5f)); - line.Position = new Vector2f(diaphragmX - 1f, centerY - winH * 0.25f); - line.FillColor = new Color(255, 255, 255, 80); - target.Draw(line); - } - - /// - /// Custom temperature‑to‑hue mapping that matches the given Sod‑tube hue values: - /// 250 K → 176, 300 K → 122, 350 K → 120?, 450 K → 71. - /// Interpolates piecewise linearly, clamping outside [250,450]. - /// - private Color TemperatureColor(double T) - { - // 1. Map temperature to hue (0–360) - double[] Tknots = { 250, 282, 353, 450 }; - double[] Hknots = { 176, 179, 122, 71 }; - double hue; - if (T <= Tknots[0]) hue = Hknots[0]; - else if (T >= Tknots[^1]) hue = Hknots[^1]; - else - { - int i = 0; - while (i < Tknots.Length - 1 && T > Tknots[i + 1]) i++; - double frac = (T - Tknots[i]) / (Tknots[i + 1] - Tknots[i]); - hue = Hknots[i] + frac * (Hknots[i + 1] - Hknots[i]); - } - - // 2. Convert hue to RGB (S = 1, V = 1) - double h = hue / 60.0; - int sector = (int)Math.Floor(h); - double f = h - sector; - byte p = 0; - byte q = (byte)(255 * (1 - f)); - byte tByte = (byte)(255 * f); - byte v = 255; - - byte r, g, b; - switch (sector % 6) - { - case 0: r = v; g = tByte; b = p; break; - case 1: r = q; g = v; b = p; break; - case 2: r = p; g = v; b = tByte; break; - case 3: r = p; g = q; b = v; break; - case 4: r = tByte; g = p; b = v; break; - default: r = v; g = p; b = q; break; - } - return new Color(r, g, b); - } - } -} \ No newline at end of file diff --git a/Scenarios/TestScenario.cs b/Scenarios/TestScenario.cs new file mode 100644 index 0000000..6fead4d --- /dev/null +++ b/Scenarios/TestScenario.cs @@ -0,0 +1,96 @@ +using System; +using SFML.Graphics; +using SFML.System; +using FluidSim.Components; +using FluidSim.Core; + +namespace FluidSim.Tests +{ + public class TestScenario : Scenario + { + private Solver solver; + private Volume0D volume; + private Pipe1D pipe; + private Atmosphere atmosphere; + private OrificeLink orificeLink; + private OpenEndLink openEndLink; + private int stepCount; + + public override void Initialize(int sampleRate) + { + double dt = 1.0 / sampleRate; + + solver = new Solver(); + solver.SetTimeStep(dt); + + volume = new Volume0D(1e-3, 150000.0, 300.0); + solver.AddComponent(volume); + + pipe = new Pipe1D(2.0, 1e-4, 20); + solver.AddComponent(pipe); + + atmosphere = new Atmosphere(); + solver.AddComponent(atmosphere); + + // Volume → left pipe end (orifice) + var volPort = volume.CreatePort(); + orificeLink = new OrificeLink(volPort, pipe, isPipeLeftEnd: true, areaProvider: () => 1e-5) + { + DischargeCoefficient = 0.62, + Gamma = volume.Gamma, + GasConstant = volume.GasConstant + }; + solver.AddOrificeLink(orificeLink); + + // Right pipe end → atmosphere (characteristic open‑end) + openEndLink = new OpenEndLink(pipe, isLeftEnd: false) + { + AmbientPressure = 101325.0, + Gamma = 1.4 + }; + solver.AddOpenEndLink(openEndLink); + + stepCount = 0; + Console.WriteLine("TestScenario initialized with sampleRate = " + sampleRate); + } + + public override float Process() + { + solver.Step(); + stepCount++; + + if (stepCount % 100 == 0) + { + double volPressure = volume.Pressure; + double volMass = volume.Mass; + double pipeLeftPressure = pipe.GetCellPressure(0); + double pipeRightPressure = pipe.GetCellPressure(pipe.CellCount - 1); + double mdotOrifice = orificeLink.LastMassFlowRate; + double mdotOpen = openEndLink.LastMassFlowRate; + + Console.WriteLine($"Step {stepCount}:"); + Console.WriteLine($" Vol Pressure = {volPressure:F1} Pa, Mass = {volMass:E4} kg"); + Console.WriteLine($" Pipe left P = {pipeLeftPressure:F1} Pa, right P = {pipeRightPressure:F1} Pa"); + Console.WriteLine($" Orifice mdot = {mdotOrifice:E4} kg/s, Open‑end mdot = {mdotOpen:E4} kg/s"); + Console.WriteLine(); + } + + // Audio sample from the open‑end mass flow + return (float)openEndLink.LastMassFlowRate; + } + + public override void Draw(RenderWindow target) + { + float winWidth = target.GetView().Size.X; + float winHeight = target.GetView().Size.Y; + + float pipeCenterY = winHeight / 2f; + float margin = 60f; + float pipeStartX = margin; + float pipeEndX = winWidth - margin; + + // Use the shared pipe drawing from the base class + DrawPipe(target, pipe, pipeCenterY, pipeStartX, pipeEndX); + } + } +} \ No newline at end of file diff --git a/Sources/EffortSource.cs b/Sources/EffortSource.cs deleted file mode 100644 index 78e186c..0000000 --- a/Sources/EffortSource.cs +++ /dev/null @@ -1,10 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Text; - -namespace FluidSim.Sources -{ - internal class EffortSource - { - } -} diff --git a/Sources/FlowSource.cs b/Sources/FlowSource.cs deleted file mode 100644 index 2004336..0000000 --- a/Sources/FlowSource.cs +++ /dev/null @@ -1,10 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Text; - -namespace FluidSim.Sources -{ - internal class FlowSource - { - } -}