From f16a1aa7635a191cd056da2fe9f0e2e703ede05d Mon Sep 17 00:00:00 2001 From: grillkol Date: Tue, 5 May 2026 11:24:32 +0200 Subject: [PATCH] insane engine sound --- Components/Pipe1D.cs | 651 +++++++++++++++++++++++++----------- Core/SoundProcessor.cs | 3 +- FluidSim.csproj | 2 +- Scenarios/EngineScenario.cs | 326 ++++++++++++------ 4 files changed, 666 insertions(+), 316 deletions(-) diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index 8aa4b6b..ac444a7 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -1,4 +1,5 @@ using System; +using System.Numerics; using FluidSim.Interfaces; namespace FluidSim.Components @@ -6,7 +7,7 @@ namespace FluidSim.Components public enum BoundaryType { OpenEnd, - ZeroPressureOpen, // pressure‑release boundary (strong reflection) + ZeroPressureOpen, ClosedEnd, GhostCell } @@ -18,72 +19,99 @@ namespace FluidSim.Components public double Area => _area; public double DampingMultiplier { get; set; } = 1.0; - private int _n; - private double _dx, _dt, _gamma, _area, _diameter; - private double[] _rho, _rhou, _E; + 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 double _rhoGhostL, _uGhostL, _pGhostL; - private double _rhoGhostR, _uGhostR, _pGhostR; + 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 double _aAmbientPressure = 101325.0; - private double _bAmbientPressure = 101325.0; + private float _aAmbientPressure = 101325f; + private float _bAmbientPressure = 101325f; - private const double CflTarget = 0.8; - private const double ReferenceSoundSpeed = 340.0; - private double _lastMassFlow = 0.0; + // 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 public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0) { - double dtGlobal = 1.0 / sampleRate; + float dtGlobal = 1f / sampleRate; int nCells; + float dxTarget = ReferenceSoundSpeed * dtGlobal / CflTarget; if (forcedCellCount > 1) nCells = forcedCellCount; else { - double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget; - nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero)); - while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1) + nCells = Math.Max(2, (int)Math.Round((float)length / dxTarget, MidpointRounding.AwayFromZero)); + while (length / nCells > dxTarget * 1.01f && nCells < int.MaxValue - 1) nCells++; } _n = nCells; - _dx = length / _n; + _dx = (float)(length / nCells); _dt = dtGlobal; - _area = area; - _gamma = 1.4; - _diameter = 2.0 * Math.Sqrt(area / Math.PI); + _area = (float)area; + _diameter = (float)(2.0 * Math.Sqrt(area / Math.PI)); + _gamma = 1.4f; - _rho = new double[_n]; - _rhou = new double[_n]; - _E = new double[_n]; + _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 PortA = new Port(); PortB = new Port(); } + // ==================== PUBLIC API (unchanged) ============================ + public void SetABoundaryType(BoundaryType type) => _aBCType = type; public void SetBBoundaryType(BoundaryType type) => _bBCType = type; - public void SetAAmbientPressure(double p) => _aAmbientPressure = p; - public void SetBAmbientPressure(double p) => _bAmbientPressure = p; + public void SetAAmbientPressure(double p) => _aAmbientPressure = (float)p; + public void SetBAmbientPressure(double p) => _bAmbientPressure = (float)p; public void SetGhostLeft(double rho, double u, double p) { - _rhoGhostL = rho; - _uGhostL = u; - _pGhostL = p; + _rhoGhostL = (float)rho; + _uGhostL = (float)u; + _pGhostL = (float)p; _ghostLSet = true; } public void SetGhostRight(double rho, double u, double p) { - _rhoGhostR = rho; - _uGhostR = u; - _pGhostR = p; + _rhoGhostR = (float)rho; + _uGhostR = (float)u; + _pGhostR = (float)p; _ghostRSet = true; } public void ClearGhostFlag() @@ -94,39 +122,73 @@ namespace FluidSim.Components public void SetUniformState(double rho, double u, double p) { - double e = p / ((_gamma - 1) * rho); - double Etot = rho * e + 0.5 * rho * u * u; + 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] = rho; - _rhou[i] = rho * u; + _rho[i] = r; + _rhou[i] = r * vel; _E[i] = Etot; } } + public int GetCellCount() => _n; + public double GetCellDensity(int i) => _rho[i]; + public double GetCellVelocity(int i) + { + float rho = Math.Max(_rho[i], 1e-12f); + return _rhou[i] / rho; + } + public double GetCellPressure(int i) + { + 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; - double rho = _rho[lastCell]; - double u = _rhou[lastCell] / Math.Max(rho, 1e-12); - double p = Pressure(lastCell); + float rho = _rho[lastCell]; + float u = _rhou[lastCell] / Math.Max(rho, 1e-12f); + float p = PressureScalar(lastCell); - double c = Math.Sqrt(_gamma * p / rho); - double uFace = u; - double rhoFace = rho; - double pFace = p; + float c = MathF.Sqrt(_gamma * p / rho); + float uFace = u; + float rhoFace = rho; + float pFace = p; - // Subsonic outflow: impose ambient pressure, adjust velocity and density - if (uFace > 0 && uFace < c) + if (uFace > 0 && uFace < c) // subsonic outflow { - double s = p / Math.Pow(rho, _gamma); - double rhoAmb = Math.Pow(_bAmbientPressure / s, 1.0 / _gamma); - double cAmb = Math.Sqrt(_gamma * _bAmbientPressure / rhoAmb); - double J_plus = u + 2.0 * c / (_gamma - 1.0); - double uFaceNew = J_plus - 2.0 * cAmb / (_gamma - 1.0); + 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; @@ -138,229 +200,410 @@ namespace FluidSim.Components public double GetAndStoreMassFlowForDerivative() { - double current = GetOpenEndMassFlow(); + float current = (float)GetOpenEndMassFlow(); double derivative = (current - _lastMassFlow) / _dt; _lastMassFlow = current; return derivative; } - public void SetCellState(int i, double rho, double u, double p) + public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8f) { - if (i < 0 || i >= _n) return; - _rho[i] = rho; - _rhou[i] = rho * u; - double e = p / ((_gamma - 1) * rho); - _E[i] = rho * e + 0.5 * rho * u * u; - } - - public int GetCellCount() => _n; - public double GetCellDensity(int i) => _rho[i]; - public double GetCellPressure(int i) => Pressure(i); - public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12); - private double Pressure(int i) => (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12)); - public double GetPressureAtFraction(double fraction) - { - int i = (int)(fraction * (_n - 1)); - i = Math.Clamp(i, 0, _n - 1); - return Pressure(i); - } - - public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8) - { - double maxW = 0.0; + float maxW = 0f; for (int i = 0; i < _n; i++) { - double rho = _rho[i]; - double u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12)); - double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12)); - double local = u + c; + 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; if (local > maxW) maxW = local; } - maxW = Math.Max(maxW, 1e-8); - return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx))); + maxW = Math.Max(maxW, 1e-8f); + return Math.Max(1, (int)Math.Ceiling((float)dtGlobal * maxW / ((float)cflTarget * _dx))); } + // ==================== MAIN SIMULATION ================================== public void SimulateSingleStep(double dtSub) { + float dt = (float)dtSub; int n = _n; - double[] Fm = new double[n + 1]; - double[] Fp = new double[n + 1]; - double[] Fe = new double[n + 1]; - // Left boundary (face 0) - double rhoL = _rho[0]; - double uL = _rhou[0] / Math.Max(rhoL, 1e-12); - double pL = Pressure(0); - if (_aBCType == BoundaryType.GhostCell && _ghostLSet) - HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, rhoL, uL, pL, out Fm[0], out Fp[0], out Fe[0]); - else if (_aBCType == BoundaryType.OpenEnd) - OpenEndFluxLeft(rhoL, uL, pL, _aAmbientPressure, out Fm[0], out Fp[0], out Fe[0]); - else if (_aBCType == BoundaryType.ZeroPressureOpen) + // --- 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]); + + // --- 2. Internal faces (1 .. n-1) – SIMD --------------------------- + int vectorSize = Vector.Count; + int lastSimdFace = n - vectorSize; // highest face index that starts a full vector block + // Face index f is between cell f-1 (left) and cell f (right). + // We want to cover faces 1..n-1. + for (int f = 1; f <= lastSimdFace; f += vectorSize) { - // Strong reflection: force pressure to ambient, extrapolate density and velocity - double rhoFace = rhoL; - double uFace = uL; - double pFace = _aAmbientPressure; - HLLCFlux(rhoFace, uFace, pFace, rhoL, uL, pL, out Fm[0], out Fp[0], out Fe[0]); + SimdInternalFluxBlock(f, vectorSize); } - else if (_aBCType == BoundaryType.ClosedEnd) - ClosedEndFlux(rhoL, uL, pL, false, out Fm[0], out Fp[0], out Fe[0]); - else - { Fm[0] = 0; Fp[0] = pL; Fe[0] = 0; } - - // Internal faces - for (int i = 0; i < n - 1; i++) + // Scalar remainder for faces f .. n-1 + for (int f = Math.Max(1, lastSimdFace + 1); f < n; f++) { - double rhoLi = _rho[i]; - double uLi = _rhou[i] / Math.Max(rhoLi, 1e-12); - double pLi = Pressure(i); - double rhoRi = _rho[i + 1]; - double uRi = _rhou[i + 1] / Math.Max(rhoRi, 1e-12); - double pRi = Pressure(i + 1); - HLLCFlux(rhoLi, uLi, pLi, rhoRi, uRi, pRi, out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]); + 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]); } - // Right boundary (face n) - double rhoR = _rho[n - 1]; - double uR = _rhou[n - 1] / Math.Max(rhoR, 1e-12); - double pR = Pressure(n - 1); - if (_bBCType == BoundaryType.GhostCell && _ghostRSet) - HLLCFlux(rhoR, uR, pR, _rhoGhostR, _uGhostR, _pGhostR, out Fm[n], out Fp[n], out Fe[n]); - else if (_bBCType == BoundaryType.OpenEnd) - OpenEndFluxRight(rhoR, uR, pR, _bAmbientPressure, out Fm[n], out Fp[n], out Fe[n]); - else if (_bBCType == BoundaryType.ZeroPressureOpen) - { - double rhoFace = rhoR; - double uFace = uR; - double pFace = _bAmbientPressure; - HLLCFlux(rhoR, uR, pR, rhoFace, uFace, pFace, out Fm[n], out Fp[n], out Fe[n]); - } - else if (_bBCType == BoundaryType.ClosedEnd) - ClosedEndFlux(rhoR, uR, pR, true, out Fm[n], out Fp[n], out Fe[n]); - else - { Fm[n] = 0; Fp[n] = pR; Fe[n] = 0; } + // --- 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]); - // Cell update (linear damping) - double radius = _diameter / 2.0; - double mu_air = 1.8e-5; - double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius); - - for (int i = 0; i < n; i++) - { - _rho[i] -= dtSub * (Fm[i + 1] - Fm[i]) / _dx; - _rhou[i] -= dtSub * (Fp[i + 1] - Fp[i]) / _dx; - _E[i] -= dtSub * (Fe[i + 1] - Fe[i]) / _dx; - - double rho = Math.Max(_rho[i], 1e-12); - double dampingFactor = Math.Exp(-(laminarCoeff / rho) * dtSub); - _rhou[i] *= dampingFactor; - - if (_rho[i] < 1e-12) _rho[i] = 1e-12; - double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i]; - double pMin = 100.0; - double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic; - if (_E[i] < eMin) _E[i] = eMin; - } + // --- 4. Cell update + damping – SIMD ------------------------------ + SimdCellUpdate(dt); } - // ---------- HLLC Riemann solver ---------- - private void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR, - out double fm, out double fp, out double fe) + // ==================== PRIVATE SCALAR HELPERS =========================== + private float PressureScalar(int i) { - double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12)); - double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12)); - double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL; - double ER = pR / ((_gamma - 1) * rR) + 0.5 * uR * uR; - double SL = Math.Min(uL - cL, uR - cR); - double SR = Math.Max(uL + cL, uR + cR); - double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) - / (rL * (SL - uL) - rR * (SR - uR)); + float rho = Math.Max(_rho[i], 1e-12f); + return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho); + } - double FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL; - double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR; + 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) + { + 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; } + } + + 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; } + } + + // ==================== 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); + + 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) { - 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))); + 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 { - 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))); + 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; } } - // ---------- Characteristic open‑end boundaries ---------- - private void OpenEndFluxLeft(double rhoInt, double uInt, double pInt, double pAmb, - out double fm, out double fp, out double fe) + private void OpenEndFluxLeft(float rhoInt, float uInt, float pInt, float pAmb, + out float fm, out float fp, out float fe) { - double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); - if (uInt <= -cInt) // supersonic inflow + 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 - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt; + fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt; return; } - if (uInt <= 0) // subsonic inflow + if (uInt <= 0) // subsonic inflow { - double T0 = 300.0, R = 287.0; - double ghost_Rho = pAmb / (R * T0); - HLLCFlux(ghost_Rho, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); + 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 - double s = pInt / Math.Pow(rhoInt, _gamma); - double ghostRho = Math.Pow(pAmb / s, 1.0 / _gamma); - double cGhost = Math.Sqrt(_gamma * pAmb / ghostRho); - double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0); - double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0); + 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; - HLLCFlux(ghostRho, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); + HLLCFluxScalar(ghostRho2, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); } - private void OpenEndFluxRight(double rhoInt, double uInt, double pInt, double pAmb, - out double fm, out double fp, out double fe) + private void OpenEndFluxRight(float rhoInt, float uInt, float pInt, float pAmb, + out float fm, out float fp, out float fe) { - double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); - if (uInt >= cInt) // supersonic outflow + 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 - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt; + fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt; return; } - if (uInt >= 0) // subsonic outflow + if (uInt >= 0) // subsonic outflow { - double s = pInt / Math.Pow(rhoInt, _gamma); - double ghost_Rho = Math.Pow(pAmb / s, 1.0 / _gamma); - double cGhost = Math.Sqrt(_gamma * pAmb / ghost_Rho); - double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0); - double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0); + 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; - HLLCFlux(rhoInt, uInt, pInt, ghost_Rho, uGhost, pAmb, out fm, out fp, out fe); + HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho, uGhost, pAmb, out fm, out fp, out fe); + return; } // subsonic inflow - double T0 = 300.0, R = 287.0; - double ghostRho = pAmb / (R * T0); - HLLCFlux(rhoInt, uInt, pInt, ghostRho, 0.0, pAmb, out fm, out fp, out fe); + 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(double rhoInt, double uInt, double pInt, bool isRight, - out double fm, out double fp, out double fe) + private void ClosedEndFlux(float rhoInt, float uInt, float pInt, bool isRight, + out float fm, out float fp, out float fe) { - double rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt; + float rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt; if (isRight) - HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe); + HLLCFluxScalar(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe); else - HLLCFlux(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe); + HLLCFluxScalar(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe); + } + + // ==================== SIMD INTERNAL FACE ROUTINE ======================== + private void SimdInternalFluxBlock(int startFace, int count) + { + // startFace is the first face index; we process 'count' consecutive faces. + // For face f, left cell = f-1, right cell = f. + // We load left and right states for faces [startFace .. startFace+count-1]. + + int leftIdx = startFace - 1; + int rightIdx = startFace; + + // Load conserved variables for left cells and right cells as vectors. + 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); + + // Derived quantities: u = ru / r, p = (gamma-1)*(E - 0.5*ru^2 / r) + 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); + + // Sound speeds + Vector cL = Vector.SquareRoot(gammaVec * pL / rL); + Vector cR = Vector.SquareRoot(gammaVec * pR / rR); + + // Wave speeds + Vector SL = Vector.Min(uL - cL, uR - cR); + Vector SR = Vector.Max(uL + cL, uR + cR); + + // Star speed + Vector num = (pR - pL) + rL * uL * (SL - uL) - rR * uR * (SR - uR); + Vector den = rL * (SL - uL) - rR * (SR - uR); + Vector Ss = num / den; + + // Total energy per unit mass (E/rho) for left/right (needed for star region) + Vector eL = EL / rL; + Vector eR = ER / rR; + + // --- Compute all four possible flux vectors --- + // Left flux + Vector Fm_L = ruL; + Vector Fp_L = ruL * uL + pL; + Vector Fe_L = (EL + pL) * uL; // (r*E + p)*u + + // Right flux + Vector Fm_R = ruR; + Vector Fp_R = ruR * uR + pR; + Vector Fe_R = (ER + pR) * uR; + + // Star‑left fluxes (when SL < 0 < Ss) + 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 (when Ss < 0 < SR) + 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; + + // --- Select the correct flux based on wave signs --- + var maskSLge0 = Vector.GreaterThanOrEqual(SL, Vector.Zero); + var maskSRle0 = Vector.LessThanOrEqual(SR, Vector.Zero); + var maskMiddle = ~(maskSLge0 | maskSRle0); // SL<0 && SR>0 + var maskStarL = maskMiddle & Vector.GreaterThanOrEqual(Ss, Vector.Zero); + var maskStarR = maskMiddle & Vector.LessThan(Ss, Vector.Zero); + + // Start with left flux, override with right/star as needed + 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)))); + + // Store to flux arrays at indices startFace .. startFace+count-1 + fm.CopyTo(_fluxM, startFace); + fp.CopyTo(_fluxP, startFace); + fe.CopyTo(_fluxE, startFace); + } + + // ==================== SIMD CELL UPDATE + DAMPING ======================== + 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‑define constants used in clamping + Vector half = new Vector(0.5f); + Vector gammaMinus1 = new Vector(_gamma - 1f); + + 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 factor = Vector.Exp(-vCoeff / r * vDt); + newRu *= factor; + + // 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 + 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 = _fluxE[i + 1] - _fluxE[i]; + + float newR = r - dt_dx * dM; + float newRu = ru - dt_dx * dP; + float newE = E - dt_dx * dE; + + // Damping + float factor = MathF.Exp(-coeff / Math.Max(r, 1e-12f) * dt); + newRu *= factor; + + // 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; + + _rho[i] = newR; + _rhou[i] = newRu; + _E[i] = newE; + } } } } \ No newline at end of file diff --git a/Core/SoundProcessor.cs b/Core/SoundProcessor.cs index 62d6ec5..9c48efb 100644 --- a/Core/SoundProcessor.cs +++ b/Core/SoundProcessor.cs @@ -82,6 +82,7 @@ namespace FluidSim.Core { // 1. Monopole source: d(mdot)/dt double derivative = (massFlow - lastMassFlow) / dt; + derivative = Math.Clamp(derivative, -500, 500); lastMassFlow = massFlow; float monopole = (float)(derivative * masterGain); @@ -115,7 +116,7 @@ namespace FluidSim.Core // 6. Dry/wet mix float output = resonant * 0.7f + wet * 0.3f; - output = Math.Clamp(output, -1f, 1f); + output = MathF.Tanh(output); return output; } } diff --git a/FluidSim.csproj b/FluidSim.csproj index 3fda142..785e8db 100644 --- a/FluidSim.csproj +++ b/FluidSim.csproj @@ -5,7 +5,7 @@ net10.0 enable enable - true + false true diff --git a/Scenarios/EngineScenario.cs b/Scenarios/EngineScenario.cs index c29fc71..01cd681 100644 --- a/Scenarios/EngineScenario.cs +++ b/Scenarios/EngineScenario.cs @@ -17,34 +17,48 @@ namespace FluidSim.Core private double ambientPressure = 101325.0; private double time; - // Crankshaft - private double crankAngle = 0.0; - private const double TargetRPM = 4000.0; - private double angularVelocity; + // ---- 4‑stroke cycle angle (0 … 4π) ---- + private double cycleCrankAngle = 0.0; // 0 to 4π, then resets + private const double TargetRPM = 1000.0; + private double angularVelocity; // rad/s of crankshaft - // Combustion - private const double CombustionPressure = 8.0 * 101325.0; - private const double CombustionTemperature = 1800.0; + // ---- Engine geometry ---- + private double bore = 0.065; // 65 mm + private double stroke = 0.0565; // 56.5 mm → 250 cc + private double conRodLength = 0.113; // roughly 2 * stroke + private double compressionRatio = 10.0; + private double V_disp; // displacement volume + private double V_clear; // clearance volume - // Valve timing - private const double ValveOpenStart = 120.0 * Math.PI / 180.0; - private const double ValveOpenEnd = 480.0 * Math.PI / 180.0; - private const double ValveRampWidth = 30.0 * Math.PI / 180.0; + // ---- Combustion ---- + private const double CombustionPressure = 50.0 * 101325.0; + private const double CombustionTemperature = 2500.0; + private bool burnInProgress = false; + private double burnStartAngle; // cycle angle when ignition began + private const double BurnDurationDeg = 40.0; + private const double BurnDurationRad = BurnDurationDeg * Math.PI / 180.0; + private double targetBurnEnergy; + private double totalBurnMass; + + // Pre‑ignition state (compressed fresh charge) for misfire restoration + private double preIgnitionMass; + private double preIgnitionInternalEnergy; + + // ---- Valve timing ---- + private const double ValveOpenStart = 120.0 * Math.PI / 180.0; // 120° after TDC power + private const double ValveOpenEnd = 480.0 * Math.PI / 180.0; // 480° ≈ 120° after TDC exhaust + private const double ValveRampWidth = 30.0 * Math.PI / 180.0; // 30° ramps private double maxOrificeArea; - // Misfire + // ---- Misfire ---- private Random rand = new Random(); private const double MisfireProbability = 0.02; private bool isMisfiring = false; - // Low‑pass filter for pressure - private double lastFilteredPressure; - private const double PressureCutoffHz = 50.0; - - // Logging + // ---- Logging ---- private int stepCount = 0; - private const int LogStepInterval = 1000; + private const int LogStepInterval = 10000; private int combustionCount = 0; private int misfireCount = 0; @@ -53,25 +67,29 @@ namespace FluidSim.Core dt = 1.0 / sampleRate; angularVelocity = TargetRPM * 2.0 * Math.PI / 60.0; - // Cylinder: 0.5 litre, initially at ambient - double cylVolume = 0.5e-3; - double initialPressure = ambientPressure; - double initialTemperature = 300.0; - cylinder = new Volume0D(cylVolume, initialPressure, initialTemperature, sampleRate) + // Displacement volume + V_disp = (Math.PI / 4.0) * bore * bore * stroke; + V_clear = V_disp / (compressionRatio - 1.0); + + // Cylinder (starts at TDC clearance volume with compressed ambient charge) + double initialPressure = ambientPressure * Math.Pow(compressionRatio, 1.4); // isentropic compression + double initialTemperature = 300.0 * Math.Pow(compressionRatio, 1.4 - 1.0); + double initialVolume = V_clear; + cylinder = new Volume0D(initialVolume, initialPressure, initialTemperature, sampleRate) { Gamma = 1.4, GasConstant = 287.0 }; - // Exhaust pipe: length 2.5 m, radius 2 cm + // Exhaust pipe (2.5 m long, 3 cm radius) double pipeLength = 2.5; - double pipeRadius = 0.02; + double pipeRadius = 0.03; double pipeArea = Math.PI * pipeRadius * pipeRadius; maxOrificeArea = pipeArea; - exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 70); + exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 100); exhaustPipe.SetUniformState(1.225, 0.0, ambientPressure); - // Coupling (valve starts closed) + // Coupling (valve initially closed) coupling = new PipeVolumeConnection(cylinder, exhaustPipe, true, orificeArea: 0.0); solver = new Solver(); @@ -79,134 +97,218 @@ namespace FluidSim.Core solver.AddVolume(cylinder); solver.AddPipe(exhaustPipe); solver.AddConnection(coupling); - // Use ZeroPressureOpen for strong reflections - solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.ZeroPressureOpen, ambientPressure); + // Open end with characteristic radiation (softer reflections) + solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.OpenEnd, ambientPressure); - // Sound processor (tuned to pipe length) + // Sound processor (keep your carefully tuned gains) soundProcessor = new SoundProcessor(sampleRate, pipeLength, pipeRadius * 2, reverbTimeMs: 200.0f); - soundProcessor.MasterGain = 0.02f; // boosted from 0.0008 - soundProcessor.PressureGain = 4.0f; // boosted from6 0.12 - soundProcessor.TurbulenceGain = 0.0002f; // reduced from 0.02 + soundProcessor.MasterGain = 0.0002f; + soundProcessor.PressureGain = 10.0f; + soundProcessor.TurbulenceGain = 0.00005f; soundProcessor.SetAmbientPressure(ambientPressure); - lastFilteredPressure = ambientPressure; - - Console.WriteLine("=== EngineScenario (ZeroPressureOpen, boosted gains) ==="); + // Log startup info + Console.WriteLine("=== EngineScenario (improved physics) ==="); Console.WriteLine($"Target RPM: {TargetRPM}, Misfire prob: {MisfireProbability * 100:F1}%"); + Console.WriteLine($"Bore x Stroke: {bore*1000:F0} x {stroke*1000:F0} mm, CR: {compressionRatio:F1}"); Console.WriteLine($"Pipe length: {pipeLength} m, fundamental: {340/(4*pipeLength):F1} Hz"); - Console.WriteLine($"Valve opens at {ValveOpenStart*180/Math.PI:F0}°, closes at {ValveOpenEnd*180/Math.PI:F0}°, ramp {ValveRampWidth*180/Math.PI:F0}°"); - Console.WriteLine($"Sample rate: {sampleRate} Hz, dt = {dt*1000:F3} ms"); - Console.WriteLine("Time[s] Crank[°] Valve[%] MassFlow[kg/s] Comb# Misfire"); - Console.WriteLine("---------------------------------------------------------"); + Console.WriteLine($"Combustion: {CombustionPressure/101325:F0} bar, {CombustionTemperature} K"); + Console.WriteLine($"Valve opens at {ValveOpenStart*180/Math.PI:F0}°, closes at {ValveOpenEnd*180/Math.PI:F0}° (ramp {ValveRampWidth*180/Math.PI:F0}°)"); + Console.WriteLine($"Burn duration: {BurnDurationDeg}°"); + Console.WriteLine("Time[s] Crank[°] Vol[cc] MassFlow[kg/s] Comb# Misfire"); + Console.WriteLine("-------------------------------------------------------------"); } - private double ValveOpenRatio(double crankRad) + // ---- Piston volume & derivative ---- + private (double volume, double dvdt) PistonKinematics(double theta) { - double cycleAngle = crankRad % (4.0 * Math.PI); - double openStart = ValveOpenStart; - double openEnd = ValveOpenEnd; + // theta = crankshaft angle (0 at TDC of power stroke) + double R = stroke / 2.0; + double cosT = Math.Cos(theta); + double sinT = Math.Sin(theta); + double L = conRodLength; - if (cycleAngle < openStart || cycleAngle > openEnd) + // Slider‑crank position relative to TDC + double s = R * (1 - cosT) + L - Math.Sqrt(L * L - R * R * sinT * sinT); + double V = V_clear + (Math.PI / 4.0) * bore * bore * s; + + // Derivative dV/dθ + double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT); + double dVdθ = (Math.PI / 4.0) * bore * bore * (R * sinT + (R * R * sinT * cosT) / sqrtTerm); + + double dvdt = dVdθ * angularVelocity; // rad/s → volume change rate + return (V, dvdt); + } + + // ---- Valve lift (trapezoidal) ---- + private double ValveOpenRatio(double cycleRad) + { + // cycleRad: 0 … 4π + if (cycleRad < ValveOpenStart || cycleRad > ValveOpenEnd) return 0.0; - double fullOpenWindow = openEnd - openStart; - double closedWindow = 2.0 * ValveRampWidth; - if (fullOpenWindow <= closedWindow) - return 1.0; + double duration = ValveOpenEnd - ValveOpenStart; + double ramp = ValveRampWidth; - double tmid = (openStart + openEnd) / 2.0; - double dist = Math.Abs(cycleAngle - tmid); - double rampHalf = (fullOpenWindow - closedWindow) / 2.0; - if (dist <= rampHalf) - return 1.0; + double t = (cycleRad - ValveOpenStart) / duration; + + if (t < ramp / duration) + return t / (ramp / duration); + else if (t > 1.0 - ramp / duration) + return (1.0 - t) / (ramp / duration); else - { - double frac = (dist - rampHalf) / ValveRampWidth; - frac = Math.Clamp(frac, 0.0, 1.0); - double lift = Math.Cos(frac * Math.PI / 2.0); - return lift * lift; - } + return 1.0; + } + + // ---- Wiebe burn fraction ---- + private double WiebeFraction(double angleFromIgnition) + { + if (angleFromIgnition >= BurnDurationRad) return 1.0; + double a = 5.0, m = 2.0; + double x = angleFromIgnition / BurnDurationRad; + return 1.0 - Math.Exp(-a * Math.Pow(x, m + 1)); } public override float Process() { - // Update crank angle - crankAngle += angularVelocity * dt; - if (crankAngle >= 2.0 * Math.PI) + // Advance cycle crank angle + cycleCrankAngle += angularVelocity * dt; + if (cycleCrankAngle >= 4.0 * Math.PI) // 720° { - crankAngle -= 2.0 * Math.PI; + cycleCrankAngle -= 4.0 * Math.PI; isMisfiring = rand.NextDouble() < MisfireProbability; - } - // Power stroke at TDC - if (crankAngle < angularVelocity * dt && crankAngle >= 0.0) - { + // ---- Prepare cylinder for new power stroke ---- + // Fill cylinder with fresh charge at BDC, then compress isentropically to TDC clearance volume. + double T_bdc = 300.0; // intake temperature + double p_bdc = ambientPressure; // intake pressure + double V_bdc = V_clear + V_disp; // volume at BDC (intake valve closing) + double freshMass = p_bdc * V_bdc / (287.0 * T_bdc); + double freshInternalEnergy = p_bdc * V_bdc / (1.4 - 1.0); + + // Compress isentropically to V_clear + double V1 = V_bdc, V2 = V_clear; + double gamma = 1.4; + double p2 = p_bdc * Math.Pow(V1 / V2, gamma); + double T2 = T_bdc * Math.Pow(V1 / V2, gamma - 1); + + // Set compressed state + cylinder.Volume = V_clear; + cylinder.Mass = freshMass; + cylinder.InternalEnergy = p2 * V_clear / (gamma - 1.0); // consistent with pressure/temperature + + // Store pre‑ignition state for misfire recovery + preIgnitionMass = cylinder.Mass; + preIgnitionInternalEnergy = cylinder.InternalEnergy; + if (isMisfiring) { - double vol = cylinder.Volume; - double R = cylinder.GasConstant; - double T0 = 300.0; - double newMass = ambientPressure * vol / (R * T0); - double newInternalEnergy = ambientPressure * vol / (cylinder.Gamma - 1.0); - cylinder.Mass = newMass; - cylinder.InternalEnergy = newInternalEnergy; + // No combustion – just expand from compressed state misfireCount++; } else { - double volume = cylinder.Volume; - double gamma = cylinder.Gamma; - double newInternalEnergy = CombustionPressure * volume / (gamma - 1.0); - double R = cylinder.GasConstant; - double newMass = CombustionPressure * volume / (R * CombustionTemperature); - cylinder.InternalEnergy = newInternalEnergy; - cylinder.Mass = newMass; + // Start Wiebe burn + double V = V_clear; + targetBurnEnergy = CombustionPressure * V / (gamma - 1.0); + double R = 287.0; + totalBurnMass = CombustionPressure * V / (R * CombustionTemperature); + burnInProgress = true; + burnStartAngle = cycleCrankAngle; // now = 0 combustionCount++; } } - // Update valve area - double valveOpen = ValveOpenRatio(crankAngle); - coupling.OrificeArea = maxOrificeArea * valveOpen; + // ---- Combustion progress (if active) ---- + if (burnInProgress) + { + double angleFromIgnition = cycleCrankAngle - burnStartAngle; + if (angleFromIgnition >= BurnDurationRad) + { + // Burn complete + cylinder.Mass = totalBurnMass; + cylinder.InternalEnergy = targetBurnEnergy; + burnInProgress = false; + } + else + { + double fraction = WiebeFraction(angleFromIgnition); + // Interpolate between pre‑ignition (compressed charge) and final burned state + double gamma = 1.4; + double V = cylinder.Volume; // still near clearance + double baseEnergy = preIgnitionInternalEnergy; + double baseMass = preIgnitionMass; + cylinder.InternalEnergy = baseEnergy * (1.0 - fraction) + targetBurnEnergy * fraction; + cylinder.Mass = baseMass * (1.0 - fraction) + totalBurnMass * fraction; + } + } + + // ---- Update cylinder volume from piston kinematics ---- + double theta = cycleCrankAngle % (2.0 * Math.PI); // crank angle for piston position + var (vol, dvdt) = PistonKinematics(theta); + cylinder.Volume = vol; + cylinder.Dvdt = dvdt; + + // ---- Valve lift & orifice area ---- + double lift = ValveOpenRatio(cycleCrankAngle); + coupling.OrificeArea = maxOrificeArea * lift; + + // ---- Solver step ---- float massFlow = solver.Step(); float endPressure = (float)exhaustPipe.GetCellPressure(exhaustPipe.GetCellCount() - 1); - // Low‑pass filter the pressure (emphasise low frequencies) - double rc = 1.0 / (2.0 * Math.PI * PressureCutoffHz); - double alpha = dt / (rc + dt); - double filteredPressure = alpha * endPressure + (1.0 - alpha) * lastFilteredPressure; - lastFilteredPressure = filteredPressure; + // ---- Audio (no filter, feed raw pressure) ---- + float audioSample = soundProcessor.Process(massFlow, endPressure); - float audioSample = soundProcessor.Process(massFlow, (float)filteredPressure); + // Log occasionally time += dt; stepCount++; - - // Logging - if (stepCount % LogStepInterval == 0 || (crankAngle < angularVelocity * dt * 2 && !isMisfiring && combustionCount > 0)) + if (stepCount % LogStepInterval == 0) { - Console.WriteLine($"{time,7:F3} {crankAngle * 180.0 / Math.PI,6:F1} " + - $"{valveOpen * 100,6:F1} {massFlow,10:F4} " + - $"{combustionCount,3} {(isMisfiring ? "X" : "")}"); + double crankDeg = cycleCrankAngle * 180.0 / Math.PI; + double volCC = cylinder.Volume * 1e6; // cc + Console.WriteLine($"{time,5:F3} {crankDeg,7:F1}° {volCC,5:F1} {massFlow,14:E4} {combustionCount,4} {misfireCount,4}"); } return audioSample; } + // ---- Drawing (unchanged) ---- 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); + } + float cylW = 80f, cylH = 150f; var cylRect = new RectangleShape(new Vector2f(cylW, cylH)); cylRect.Position = new Vector2f(40f, centerY - cylH / 2f); - double pNorm = (cylinder.Pressure - ambientPressure) / ambientPressure; - if (double.IsNaN(pNorm)) pNorm = 0; - byte red = (byte)(Math.Clamp(pNorm * 128, 0, 255)); - byte blue = (byte)(Math.Clamp(-pNorm * 128, 0, 255)); - cylRect.FillColor = new Color(red, 0, blue); + + double tempCyl = cylinder.Temperature; // Volume0D now has Temperature + float tnCyl = NormaliseTemperature(tempCyl); + byte redCyl = (byte)(tnCyl > 0 ? 255 * tnCyl : 0); + byte blueCyl = (byte)(tnCyl < 0 ? -255 * tnCyl : 0); + byte greenCyl = (byte)(255 * (1 - Math.Abs(tnCyl))); + cylRect.FillColor = new Color(redCyl, greenCyl, blueCyl); target.Draw(cylRect); int n = exhaustPipe.GetCellCount(); @@ -215,21 +317,25 @@ namespace FluidSim.Core float dx = pipeLen / (n - 1); float baseRadius = 20f; var vertices = new Vertex[n * 2]; + float ambientPressure = 101325f; + for (int i = 0; i < n; i++) { float x = pipeStartX + i * dx; - double p = exhaustPipe.GetCellPressure(i); - float r = baseRadius * (float)(1.0 + (p - ambientPressure) / ambientPressure); + double p = exhaustPipe.GetCellPressure(i); + double rho = exhaustPipe.GetCellDensity(i); + double T = p / (rho * R); + + float r = baseRadius * 0.1f * (float)(1.0 + (p - ambientPressure) / ambientPressure); if (r < 2f) r = 2f; - double t = (p - ambientPressure) / ambientPressure; - t = Math.Clamp(t, -1.0, 1.0); - byte rCol = (byte)(t > 0 ? 255 * t : 0); - byte bCol = (byte)(t < 0 ? -255 * t : 0); - byte gCol = (byte)(255 * (1 - Math.Abs(t))); + float tn = NormaliseTemperature(T); + byte rCol = (byte)(tn > 0 ? 255 * tn : 0); + byte bCol = (byte)(tn < 0 ? -255 * tn : 0); + byte gCol = (byte)(255 * (1 - Math.Abs(tn))); var col = new Color(rCol, gCol, bCol); - vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col); + 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);