From d963032e745d794dcf560ece9fb3bbd8b3b79f28 Mon Sep 17 00:00:00 2001 From: grillkol Date: Tue, 5 May 2026 10:32:30 +0200 Subject: [PATCH] General testing --- .vscode/launch.json | 26 ++ .vscode/tasks.json | 41 +++ Components/Pipe1D.cs | 440 +++++++++++------------- Components/Volume0D.cs | 65 ++-- Core/NozzleFlow.cs | 46 +++ Core/PipeVolumeConnection.cs | 20 ++ Core/Solver.cs | 157 ++------- Core/SoundProcessor.cs | 166 ++++++++- Program.cs | 22 +- Scenarios/EngineScenario.cs | 238 +++++++++++++ Scenarios/HelmholtzResonatorScenario.cs | 21 +- 11 files changed, 794 insertions(+), 448 deletions(-) create mode 100644 .vscode/launch.json create mode 100644 .vscode/tasks.json create mode 100644 Core/NozzleFlow.cs create mode 100644 Core/PipeVolumeConnection.cs create mode 100644 Scenarios/EngineScenario.cs diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..c342031 --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,26 @@ +{ + "version": "0.2.0", + "configurations": [ + { + // Use IntelliSense to find out which attributes exist for C# debugging + // Use hover for the description of the existing attributes + // For further information visit https://github.com/dotnet/vscode-csharp/blob/main/debugger-launchjson.md + "name": ".NET Core Launch (console)", + "type": "coreclr", + "request": "launch", + "preLaunchTask": "build", + // If you have changed target frameworks, make sure to update the program path. + "program": "${workspaceFolder}/bin/Debug/net10.0/FluidSim.dll", + "args": [], + "cwd": "${workspaceFolder}", + // For more information about the 'console' field, see https://aka.ms/VSCode-CS-LaunchJson-Console + "console": "internalConsole", + "stopAtEntry": false + }, + { + "name": ".NET Core Attach", + "type": "coreclr", + "request": "attach" + } + ] +} \ No newline at end of file diff --git a/.vscode/tasks.json b/.vscode/tasks.json new file mode 100644 index 0000000..f31e253 --- /dev/null +++ b/.vscode/tasks.json @@ -0,0 +1,41 @@ +{ + "version": "2.0.0", + "tasks": [ + { + "label": "build", + "command": "dotnet", + "type": "process", + "args": [ + "build", + "${workspaceFolder}/FluidSim.csproj", + "/property:GenerateFullPaths=true", + "/consoleloggerparameters:NoSummary;ForceNoAlign" + ], + "problemMatcher": "$msCompile" + }, + { + "label": "publish", + "command": "dotnet", + "type": "process", + "args": [ + "publish", + "${workspaceFolder}/FluidSim.csproj", + "/property:GenerateFullPaths=true", + "/consoleloggerparameters:NoSummary;ForceNoAlign" + ], + "problemMatcher": "$msCompile" + }, + { + "label": "watch", + "command": "dotnet", + "type": "process", + "args": [ + "watch", + "run", + "--project", + "${workspaceFolder}/FluidSim.csproj" + ], + "problemMatcher": "$msCompile" + } + ] +} \ No newline at end of file diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index 4ef0bcc..8aa4b6b 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -5,9 +5,10 @@ namespace FluidSim.Components { public enum BoundaryType { - VolumeCoupling, OpenEnd, - ClosedEnd + ZeroPressureOpen, // pressure‑release boundary (strong reflection) + ClosedEnd, + GhostCell } public class Pipe1D @@ -21,27 +22,20 @@ namespace FluidSim.Components private double _dx, _dt, _gamma, _area, _diameter; private double[] _rho, _rhou, _E; - // Volume‑coupling ghost states for boundaries A and B - private double _rhoA, _pA; - private double _rhoB, _pB; - private bool _aBCSet, _bBCSet; + // Ghost cell states + private double _rhoGhostL, _uGhostL, _pGhostL; + private double _rhoGhostR, _uGhostR, _pGhostR; + private bool _ghostLSet, _ghostRSet; - private BoundaryType _aBCType = BoundaryType.VolumeCoupling; - private BoundaryType _bBCType = BoundaryType.VolumeCoupling; + private BoundaryType _aBCType = BoundaryType.GhostCell; + private BoundaryType _bBCType = BoundaryType.GhostCell; private double _aAmbientPressure = 101325.0; private double _bAmbientPressure = 101325.0; private const double CflTarget = 0.8; private const double ReferenceSoundSpeed = 340.0; - - 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); - - public BoundaryType ABCType => _aBCType; - public BoundaryType BBCType => _bBCType; + private double _lastMassFlow = 0.0; public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0) { @@ -49,9 +43,7 @@ namespace FluidSim.Components int nCells; if (forcedCellCount > 1) - { nCells = forcedCellCount; - } else { double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget; @@ -65,8 +57,6 @@ namespace FluidSim.Components _dt = dtGlobal; _area = area; _gamma = 1.4; - - // Hydraulic diameter for a circular pipe _diameter = 2.0 * Math.Sqrt(area / Math.PI); _rho = new double[_n]; @@ -82,6 +72,26 @@ namespace FluidSim.Components public void SetAAmbientPressure(double p) => _aAmbientPressure = p; public void SetBAmbientPressure(double p) => _bAmbientPressure = p; + public void SetGhostLeft(double rho, double u, double p) + { + _rhoGhostL = rho; + _uGhostL = u; + _pGhostL = p; + _ghostLSet = true; + } + public void SetGhostRight(double rho, double u, double p) + { + _rhoGhostR = rho; + _uGhostR = u; + _pGhostR = p; + _ghostRSet = true; + } + public void ClearGhostFlag() + { + _ghostLSet = false; + _ghostRSet = false; + } + public void SetUniformState(double rho, double u, double p) { double e = p / ((_gamma - 1) * rho); @@ -94,6 +104,46 @@ namespace FluidSim.Components } } + 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); + + double c = Math.Sqrt(_gamma * p / rho); + double uFace = u; + double rhoFace = rho; + double pFace = p; + + // Subsonic outflow: impose ambient pressure, adjust velocity and density + if (uFace > 0 && uFace < c) + { + 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); + if (uFaceNew > 0) uFace = uFaceNew; + if (uFace < 0) uFace = 0; + rhoFace = rhoAmb; + pFace = _bAmbientPressure; + } + + return rhoFace * uFace * _area; + } + + public double GetAndStoreMassFlowForDerivative() + { + double current = GetOpenEndMassFlow(); + double derivative = (current - _lastMassFlow) / _dt; + _lastMassFlow = current; + return derivative; + } + public void SetCellState(int i, double rho, double u, double p) { if (i < 0 || i >= _n) return; @@ -103,22 +153,18 @@ namespace FluidSim.Components _E[i] = rho * e + 0.5 * rho * u * u; } - public void SetAVolumeState(double rhoVol, double pVol) + 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) { - _rhoA = rhoVol; - _pA = pVol; - _aBCSet = true; + int i = (int)(fraction * (_n - 1)); + i = Math.Clamp(i, 0, _n - 1); + return Pressure(i); } - public void SetBVolumeState(double rhoVol, double pVol) - { - _rhoB = rhoVol; - _pB = pVol; - _bBCSet = true; - } - - public void ClearBC() => _aBCSet = _bBCSet = false; - public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8) { double maxW = 0.0; @@ -141,101 +187,72 @@ namespace FluidSim.Components double[] Fp = new double[n + 1]; double[] Fe = new double[n + 1]; - // ---------- Boundary A (face 0, left) ---------- - double rhoIntA = _rho[0]; - double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12); - double pIntA = Pressure(0); - - switch (_aBCType) + // 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) { - case BoundaryType.VolumeCoupling: - if (_aBCSet) - { - HLLCFlux(_rhoA, 0.0, _pA, - rhoIntA, uIntA, pIntA, - out Fm[0], out Fp[0], out Fe[0]); - } - else - { - Fm[0] = 0; Fp[0] = pIntA; Fe[0] = 0; - } - break; - - case BoundaryType.OpenEnd: - OpenEndFluxA(rhoIntA, uIntA, pIntA, _aAmbientPressure, - out Fm[0], out Fp[0], out Fe[0]); - break; - - case BoundaryType.ClosedEnd: - ClosedEndFlux(rhoIntA, uIntA, pIntA, isRightBoundary: false, - out Fm[0], out Fp[0], out Fe[0]); - break; + // 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]); } + 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 ---------- + // Internal faces for (int i = 0; i < n - 1; i++) { - double rhoL = _rho[i]; - double uL = _rhou[i] / Math.Max(rhoL, 1e-12); - double pL = Pressure(i); - - double rhoR = _rho[i + 1]; - double uR = _rhou[i + 1] / Math.Max(rhoR, 1e-12); - double pR = Pressure(i + 1); - - HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, - out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]); + 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]); } - // ---------- Boundary B (face n, right) ---------- - double rhoIntB = _rho[n - 1]; - double uIntB = _rhou[n - 1] / Math.Max(rhoIntB, 1e-12); - double pIntB = Pressure(n - 1); - - switch (_bBCType) + // 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) { - case BoundaryType.VolumeCoupling: - if (_bBCSet) - { - HLLCFlux(rhoIntB, uIntB, pIntB, - _rhoB, 0.0, _pB, - out Fm[n], out Fp[n], out Fe[n]); - } - else - { - Fm[n] = 0; Fp[n] = pIntB; Fe[n] = 0; - } - break; - - case BoundaryType.OpenEnd: - OpenEndFluxB(rhoIntB, uIntB, pIntB, _bAmbientPressure, - out Fm[n], out Fp[n], out Fe[n]); - break; - - case BoundaryType.ClosedEnd: - ClosedEndFlux(rhoIntB, uIntB, pIntB, isRightBoundary: true, - out Fm[n], out Fp[n], out Fe[n]); - break; + 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; } - // ---- Cell update with linear laminar damping ---- + // 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++) { - double dM = (Fm[i + 1] - Fm[i]) / _dx; - double dP = (Fp[i + 1] - Fp[i]) / _dx; - double dE = (Fe[i + 1] - Fe[i]) / _dx; - - _rho[i] -= dtSub * dM; - _rhou[i] -= dtSub * dP; - _E[i] -= dtSub * dE; + _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 dampingRate = laminarCoeff / rho; - double dampingFactor = Math.Exp(-dampingRate * dtSub); + double dampingFactor = Math.Exp(-(laminarCoeff / rho) * dtSub); _rhou[i] *= dampingFactor; if (_rho[i] < 1e-12) _rho[i] = 1e-12; @@ -244,139 +261,10 @@ namespace FluidSim.Components double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic; if (_E[i] < eMin) _E[i] = eMin; } - - // ---------- Port quantities ---------- - double mdotA_sub = _aBCType == BoundaryType.VolumeCoupling && _aBCSet ? Fm[0] * _area : 0.0; - double mdotB_sub = _bBCType == BoundaryType.VolumeCoupling && _bBCSet ? -Fm[n] * _area : 0.0; - - PortA.MassFlowRate = mdotA_sub; - PortB.MassFlowRate = mdotB_sub; - PortA.Pressure = pIntA; - PortB.Pressure = pIntB; - PortA.Density = _rho[0]; - PortB.Density = _rho[n - 1]; - - // Corrected enthalpy for both directions - if (_aBCType == BoundaryType.VolumeCoupling && _aBCSet) - { - PortA.SpecificEnthalpy = mdotA_sub < 0 - ? GetCellTotalSpecificEnthalpy(0) - : (_gamma / (_gamma - 1.0)) * _pA / Math.Max(_rhoA, 1e-12); - } - if (_bBCType == BoundaryType.VolumeCoupling && _bBCSet) - { - PortB.SpecificEnthalpy = mdotB_sub < 0 - ? GetCellTotalSpecificEnthalpy(_n - 1) - : (_gamma / (_gamma - 1.0)) * _pB / Math.Max(_rhoB, 1e-12); - } } - private double GetCellTotalSpecificEnthalpy(int i) - { - double rho = Math.Max(_rho[i], 1e-12); - double u = _rhou[i] / rho; - double p = Pressure(i); - double h = _gamma / (_gamma - 1.0) * p / rho; - return h + 0.5 * u * u; - } - - private double Pressure(int i) => - (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12)); - - // ========== Characteristic‑based Open End ========== - private void OpenEndFluxA(double rhoInt, double uInt, double pInt, double pAmb, - out double fm, out double fp, out double fe) - { - double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); - - // Subsonic inflow (uInt ≤ 0, so flow inside pipe ←) - if (uInt <= -cInt) // supersonic inflow – use interior state as ghost - { - fm = rhoInt * uInt; - fp = rhoInt * uInt * uInt + pInt; - fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt; - return; - } - else if (uInt <= 0) // subsonic inflow - { - // Reservoir condition: p = pAmb, T = 300K, u = 0 - double T0 = 300.0; - double R = 287.0; - double rhoGhost = pAmb / (R * T0); - HLLCFlux(rhoGhost, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); - return; - } - else // subsonic outflow (uInt > 0) - { - // Ghost pressure forced to pAmb - double s = pInt / Math.Pow(rhoInt, _gamma); - double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma); - double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost); - - // Outgoing Riemann invariant J⁻ = uInt - 2*cInt/(γ-1) (for left boundary) - double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0); - double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0); - - // Prevent spurious inflow by clipping to zero - if (uGhost < 0) uGhost = 0; - - HLLCFlux(rhoGhost, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); - } - } - - private void OpenEndFluxB(double rhoInt, double uInt, double pInt, double pAmb, - out double fm, out double fp, out double fe) - { - double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); - - if (uInt >= cInt) // supersonic outflow (extrapolation) - { - fm = rhoInt * uInt; - fp = rhoInt * uInt * uInt + pInt; - fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt; - return; - } - else if (uInt >= 0) // subsonic outflow - { - double s = pInt / Math.Pow(rhoInt, _gamma); - double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma); - double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost); - - // Outgoing Riemann invariant J⁺ = uInt + 2*cInt/(γ-1) (for right boundary) - double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0); - double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0); - - // Clip to zero to prevent inflow - if (uGhost > 0) uGhost = 0; - - HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pAmb, out fm, out fp, out fe); - } - else // subsonic inflow - { - double T0 = 300.0; - double R = 287.0; - double rhoGhost = pAmb / (R * T0); - HLLCFlux(rhoInt, uInt, pInt, rhoGhost, 0.0, pAmb, out fm, out fp, out fe); - } - } - - // ========== Closed end (mirror) ========== - private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRightBoundary, - out double fm, out double fp, out double fe) - { - double rhoGhost = rhoInt; - double pGhost = pInt; - double uGhost = -uInt; // mirror velocity - - if (isRightBoundary) - HLLCFlux(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); - } - - // ========== Standard HLLC flux ========== - private void HLLCFlux(double rL, double uL, double pL, - double rR, double uR, double pR, + // ---------- 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) { double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12)); @@ -385,7 +273,6 @@ namespace FluidSim.Components 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)); @@ -404,17 +291,76 @@ namespace FluidSim.Components else { double rsR = rR * (SR - uR) / (SR - Ss); - double ps = pL + rL * (SL - uL) * (Ss - uL); + 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; } } - public double GetPressureAtFraction(double fraction) + // ---------- Characteristic open‑end boundaries ---------- + private void OpenEndFluxLeft(double rhoInt, double uInt, double pInt, double pAmb, + out double fm, out double fp, out double fe) { - int i = (int)(fraction * (_n - 1)); - i = Math.Clamp(i, 0, _n - 1); - return Pressure(i); + double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); + 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; + return; + } + 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); + 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); + if (uGhost < 0) uGhost = 0; + HLLCFlux(ghostRho, 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) + { + double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); + 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; + return; + } + 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); + if (uGhost > 0) uGhost = 0; + HLLCFlux(rhoInt, uInt, pInt, ghost_Rho, uGhost, pAmb, out fm, out fp, out fe); + } + // 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); + } + + private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRight, + out double fm, out double fp, out double fe) + { + double rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt; + if (isRight) + HLLCFlux(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); } } } \ No newline at end of file diff --git a/Components/Volume0D.cs b/Components/Volume0D.cs index 2ff3940..6fee6d9 100644 --- a/Components/Volume0D.cs +++ b/Components/Volume0D.cs @@ -1,21 +1,17 @@ using System; -using FluidSim.Interfaces; -using FluidSim.Utils; namespace FluidSim.Components { public class Volume0D { - public Port Port { get; private set; } - - public double Mass { get; private set; } - public double InternalEnergy { get; private set; } + public double Mass { get; set; } + public double InternalEnergy { 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; } + public double Dvdt { get; set; } private double _dt; @@ -24,6 +20,9 @@ namespace FluidSim.Components public double Temperature => Pressure / (Density * GasConstant); public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density; + 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) @@ -31,54 +30,38 @@ namespace FluidSim.Components GasConstant = gasConstant; Gamma = gamma; Volume = initialVolume; - dVdt = 0.0; + Dvdt = 0.0; _dt = 1.0 / sampleRate; double rho0 = initialPressure / (GasConstant * initialTemperature); Mass = rho0 * Volume; InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0); - - Port = new Port(); - PushStateToPort(); } - public void PushStateToPort() - { - Port.Pressure = Pressure; - Port.Density = Density; - Port.Temperature = Temperature; - Port.SpecificEnthalpy = SpecificEnthalpy; - } - - // Original integrate (uses the constructor’s sample rate) - public void Integrate() - { - Integrate(_dt); - } - - public void SetPressure(double newPressure) - { - InternalEnergy = newPressure * Volume / (Gamma - 1.0); - // Mass stays the same, so density is unchanged - } - - // New overload: integrate with a custom time step (for sub‑steps) public void Integrate(double dtOverride) { - double mdot = Port.MassFlowRate; - double h_in = Port.SpecificEnthalpy; - - double dm = mdot * dtOverride; - double dE = (mdot * h_in) * dtOverride - Pressure * dVdt * dtOverride; + double dm = MassFlowRateIn * dtOverride; + double dE = (MassFlowRateIn * SpecificEnthalpyIn) * dtOverride - Pressure * Dvdt * dtOverride; Mass += dm; InternalEnergy += dE; - // Hard physical bounds – prevent NaN and unphysical states - if (Mass < 1e-12) Mass = 1e-12; - if (InternalEnergy < 1e-12) InternalEnergy = 1e-12; + // Safety: if mass becomes extremely small, reset internal energy to zero + if (Mass < 1e-12) + { + Mass = 0.0; + InternalEnergy = 0.0; + } + else if (InternalEnergy < 1e-12) + { + InternalEnergy = 0.0; + } - PushStateToPort(); + // Avoid negative mass/energy + if (Mass < 0.0) Mass = 0.0; + if (InternalEnergy < 0.0) InternalEnergy = 0.0; } + + public void Integrate() => Integrate(_dt); } } \ No newline at end of file diff --git a/Core/NozzleFlow.cs b/Core/NozzleFlow.cs new file mode 100644 index 0000000..5b40733 --- /dev/null +++ b/Core/NozzleFlow.cs @@ -0,0 +1,46 @@ +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) + { + // Default fallback (no flow) + 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 < 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; + if (double.IsNaN(massFlow) || double.IsInfinity(massFlow)) + massFlow = 0.0; + } + } +} \ No newline at end of file diff --git a/Core/PipeVolumeConnection.cs b/Core/PipeVolumeConnection.cs new file mode 100644 index 0000000..40eb6f7 --- /dev/null +++ b/Core/PipeVolumeConnection.cs @@ -0,0 +1,20 @@ +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 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 a85f51d..e2efaad 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -1,7 +1,6 @@ using System; using System.Collections.Generic; using FluidSim.Components; -using FluidSim.Interfaces; namespace FluidSim.Core { @@ -9,162 +8,80 @@ namespace FluidSim.Core { private readonly List _volumes = new(); private readonly List _pipes = new(); - private readonly List _connections = new(); + private readonly List _connections = 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(Connection c) => _connections.Add(c); + public void AddConnection(PipeVolumeConnection c) => _connections.Add(c); public void SetTimeStep(double dt) => _dt = dt; - /// - /// Set boundary type for a pipe end. isA = true for port A (left), false for port B (right). - /// public void SetPipeBoundary(Pipe1D pipe, bool isA, BoundaryType type, double ambientPressure = 101325.0) { if (isA) { pipe.SetABoundaryType(type); - if (type == BoundaryType.OpenEnd) - pipe.SetAAmbientPressure(ambientPressure); + if (type == BoundaryType.OpenEnd) pipe.SetAAmbientPressure(ambientPressure); } else { pipe.SetBBoundaryType(type); - if (type == BoundaryType.OpenEnd) - pipe.SetBAmbientPressure(ambientPressure); + if (type == BoundaryType.OpenEnd) pipe.SetBAmbientPressure(ambientPressure); } } public float Step() { - // 1. Volumes publish state - foreach (var v in _volumes) - v.PushStateToPort(); - - // 2. Set volume BCs for volume‑coupled ends + // 1. Compute nozzle flows and update volumes (once per audio sample) foreach (var conn in _connections) { - if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) - { - var pipe = GetPipe(conn.PortA); - bool isA = pipe.PortA == conn.PortA; - if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) || - (!isA && pipe.BBCType == BoundaryType.VolumeCoupling)) - SetVolumeBC(conn.PortA, conn.PortB); - } - else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB)) - { - var pipe = GetPipe(conn.PortB); - bool isA = pipe.PortB == conn.PortB; - if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) || - (!isA && pipe.BBCType == BoundaryType.VolumeCoupling)) - SetVolumeBC(conn.PortB, conn.PortA); - } + double downstreamPressure = conn.IsPipeLeftEnd + ? conn.Pipe.GetCellPressure(0) + : conn.Pipe.GetCellPressure(conn.Pipe.GetCellCount() - 1); + + NozzleFlow.Compute(conn.Volume, conn.OrificeArea, downstreamPressure, + out double mdot, out double rhoFace, out double uFace, out double pFace, + gamma: conn.Volume.Gamma); + + // Limit mass flow to available mass + double maxMdot = conn.Volume.Mass / _dt; + if (mdot > maxMdot) mdot = maxMdot; + if (mdot < -maxMdot) mdot = -maxMdot; + + conn.Volume.MassFlowRateIn = -mdot; + conn.Volume.SpecificEnthalpyIn = (conn.Volume.Gamma / (conn.Volume.Gamma - 1.0)) * + (conn.Volume.Pressure / Math.Max(conn.Volume.Density, 1e-12)); + conn.Volume.Integrate(_dt); + + if (conn.IsPipeLeftEnd) + conn.Pipe.SetGhostLeft(rhoFace, uFace, pFace); + else + conn.Pipe.SetGhostRight(rhoFace, uFace, pFace); } - // 3. Sub‑steps + // 2. Determine required sub‑steps int nSub = 1; foreach (var p in _pipes) nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt)); double dtSub = _dt / nSub; + // 3. Sub‑step loop for pipes for (int sub = 0; sub < nSub; sub++) - { foreach (var p in _pipes) p.SimulateSingleStep(dtSub); - foreach (var conn in _connections) - { - if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) - { - var pipe = GetPipe(conn.PortA); - bool isA = pipe.PortA == conn.PortA; - if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) || - (!isA && pipe.BBCType == BoundaryType.VolumeCoupling)) - TransferAndIntegrate(conn.PortA, conn.PortB, dtSub); - } - else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB)) - { - var pipe = GetPipe(conn.PortB); - bool isA = pipe.PortB == conn.PortB; - if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) || - (!isA && pipe.BBCType == BoundaryType.VolumeCoupling)) - TransferAndIntegrate(conn.PortB, conn.PortA, dtSub); - } - } - - if (sub < nSub - 1) - { - foreach (var v in _volumes) - v.PushStateToPort(); - - foreach (var conn in _connections) - { - if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) - { - var pipe = GetPipe(conn.PortA); - bool isA = pipe.PortA == conn.PortA; - if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) || - (!isA && pipe.BBCType == BoundaryType.VolumeCoupling)) - SetVolumeBC(conn.PortA, conn.PortB); - } - else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB)) - { - var pipe = GetPipe(conn.PortB); - bool isA = pipe.PortB == conn.PortB; - if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) || - (!isA && pipe.BBCType == BoundaryType.VolumeCoupling)) - SetVolumeBC(conn.PortB, conn.PortA); - } - } - } - } - - // 5. Audio samples (none for now, but placeholder) - var audioSamples = new List(); - foreach (var conn in _connections) - { - if (conn is SoundConnection sc) - audioSamples.Add(sc.GetAudioSample()); - } - - // 6. Clear BC flags + // 4. Clear ghost flags foreach (var p in _pipes) - p.ClearBC(); + p.ClearGhostFlag(); - return SoundProcessor.MixAndClip(audioSamples.ToArray()); - } + // 5. Return raw mass flow from the first pipe’s open end (assumed exhaust tailpipe) + if (_pipes.Count > 0) + return (float)_pipes[0].GetOpenEndMassFlow(); - private bool IsVolumePort(Port p) => _volumes.Exists(v => v.Port == p); - private bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p); - private Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p); - private Volume0D GetVolume(Port p) => _volumes.Find(v => v.Port == p); - - private void SetVolumeBC(Port pipePort, Port volPort) - { - var pipe = GetPipe(pipePort); - if (pipe == null) return; - bool isA = pipe.PortA == pipePort; - if (isA) - pipe.SetAVolumeState(volPort.Density, volPort.Pressure); - else - pipe.SetBVolumeState(volPort.Density, volPort.Pressure); - } - - private void TransferAndIntegrate(Port pipePort, Port volPort, double dtSub) - { - double mdot = pipePort.MassFlowRate; - volPort.MassFlowRate = -mdot; - - if (mdot < 0) // pipe → volume - { - volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy; - } - // else volume’s own enthalpy (from PushStateToPort) is used - - GetVolume(volPort)?.Integrate(dtSub); + return 0f; } } } \ No newline at end of file diff --git a/Core/SoundProcessor.cs b/Core/SoundProcessor.cs index d67d54f..62d6ec5 100644 --- a/Core/SoundProcessor.cs +++ b/Core/SoundProcessor.cs @@ -1,23 +1,155 @@ -namespace FluidSim.Core -{ - /// - /// Mixes multiple audio samples and applies a soft‑clipping tanh. - /// - public static class SoundProcessor - { - /// Overall gain applied after mixing (before tanh). - public static float MasterGain { get; set; } = 0.01f; +using System; - /// - /// Mixes an array of raw audio samples and returns a single sample in [‑1, 1]. - /// - public static float MixAndClip(params float[] samples) +namespace FluidSim.Core +{ + public class SoundProcessor + { + // Monopole state + private double lastMassFlow = 0.0; + private double dt; + + // Resonant band‑pass filter (second‑order) + private double b0, b1, b2, a1, a2; + private double x1 = 0, x2 = 0, y1 = 0, y2 = 0; + private double pipeLength; + + // Reverb (outdoor) + private float[] delayLine; + private int writeIndex; + private float feedback = 0.50f; + private float lowpassCoeff = 0.70f; + private float lastFeedbackSample = 0f; + + // Turbulence (pink noise scaled by U³) + private PinkNoiseGenerator pinkNoise; + private float turbulenceGain = 0.05f; + private double pipeArea; + private double ambientPressure = 101325.0; + + // Gains + private float masterGain = 0.0005f; + private float pressureGain = 0.12f; + + public SoundProcessor(int sampleRate, double pipeLengthMeters, double pipeDiameterMeters = 0.04, float reverbTimeMs = 200.0f) { + dt = 1.0 / sampleRate; + pipeLength = pipeLengthMeters; + pipeArea = Math.PI * Math.Pow(pipeDiameterMeters / 2.0, 2.0); + + // Design resonant filter at pipe fundamental frequency + double c = 340.0; + double f0 = c / (4.0 * pipeLength); + double Q = 15.0; + double omega = 2.0 * Math.PI * f0; + double alpha = Math.Sin(omega * dt) / (2.0 * Q); + double norm = 1.0 / (1.0 + alpha); + b0 = alpha * norm; + b1 = 0.0; + b2 = -alpha * norm; + a1 = -2.0 * Math.Cos(omega * dt) * norm; + a2 = (1.0 - alpha) * norm; + + // Reverb delay line + int delaySamples = (int)(sampleRate * reverbTimeMs / 1000.0); + delayLine = new float[delaySamples]; + writeIndex = 0; + + pinkNoise = new PinkNoiseGenerator(); + } + + public float MasterGain + { + get => masterGain; + set => masterGain = value; + } + + public float PressureGain + { + get => pressureGain; + set => pressureGain = value; + } + + public float TurbulenceGain + { + get => turbulenceGain; + set => turbulenceGain = value; + } + + public void SetAmbientPressure(double p) => ambientPressure = p; + public void SetPipeDiameter(double diameterMeters) => pipeArea = Math.PI * Math.Pow(diameterMeters / 2.0, 2.0); + + public float Process(float massFlow, float pipeEndPressure) + { + // 1. Monopole source: d(mdot)/dt + double derivative = (massFlow - lastMassFlow) / dt; + lastMassFlow = massFlow; + float monopole = (float)(derivative * masterGain); + + // 2. Pressure difference (low‑frequency component) + float pressureDiff = (float)((pipeEndPressure - ambientPressure) / ambientPressure) * pressureGain; + float mixed = monopole + pressureDiff; + // DO NOT clamp here – let the filter and final clamp handle dynamics + + // 3. Resonant band‑pass filter + double y = b0 * mixed + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2; + x2 = x1; x1 = mixed; + y2 = y1; y1 = y; + float resonant = (float)Math.Clamp(y, -1f, 1f); + + // 4. Turbulence noise: amplitude ∝ U³ (empirical for low speeds) + double velocity = massFlow / (pipeArea * 1.225); + double Uref = 100.0; + double turbulenceAmp = Math.Pow(Math.Abs(velocity) / Uref, 3.0); + float pink = pinkNoise.Next() * turbulenceGain * (float)turbulenceAmp; + + resonant += pink; + resonant = Math.Clamp(resonant, -1f, 1f); + + // 5. Outdoor reverb + float delayed = delayLine[writeIndex]; + float filteredDelay = delayed * lowpassCoeff + lastFeedbackSample * (1f - lowpassCoeff); + lastFeedbackSample = filteredDelay; + float wet = delayed + filteredDelay * feedback; + delayLine[writeIndex] = resonant + filteredDelay * feedback; + writeIndex = (writeIndex + 1) % delayLine.Length; + + // 6. Dry/wet mix + float output = resonant * 0.7f + wet * 0.3f; + output = Math.Clamp(output, -1f, 1f); + return output; + } + } + + internal class PinkNoiseGenerator + { + private readonly Random random = new Random(); + private readonly float[] whiteNoise = new float[7]; + private int currentIndex = 0; + + public PinkNoiseGenerator() + { + for (int i = 0; i < 7; i++) + whiteNoise[i] = (float)(random.NextDouble() * 2.0 - 1.0); + } + + public float Next() + { + whiteNoise[0] = (float)(random.NextDouble() * 2.0 - 1.0); + currentIndex = (currentIndex + 1) & 0x7F; + int updateMask = 0; + int temp = currentIndex; + for (int i = 0; i < 7; i++) + { + if ((temp & 1) == 0) + updateMask |= (1 << i); + temp >>= 1; + } + for (int i = 1; i < 7; i++) + if ((updateMask & (1 << i)) != 0) + whiteNoise[i] = (float)(random.NextDouble() * 2.0 - 1.0); float sum = 0f; - foreach (float s in samples) - sum += s; - sum *= MasterGain; - return sum; + for (int i = 0; i < 7; i++) sum += whiteNoise[i]; + return sum / 3.5f; } } } \ No newline at end of file diff --git a/Program.cs b/Program.cs index c7b75cb..11a2048 100644 --- a/Program.cs +++ b/Program.cs @@ -13,23 +13,23 @@ public class Program private static Scenario scenario; // Speed control - //private static double desiredSpeed = 1.0; private static double desiredSpeed = 0.0001; + //private static double desiredSpeed = 1; private static double currentSpeed = desiredSpeed; private const double MinSpeed = 0.0001; private const double MaxSpeed = 1.0; private const double ScrollFactor = 1.1; // Space‑toggle state - private static double lastDesiredSpeed = 0.1; // remembers the last non‑1.0 scroll speed - private static bool isRealTime = true; // true when desiredSpeed == 1.0 + private static double lastDesiredSpeed = 0.1; // remembers the last non‑1.0 speed + private static bool isRealTime = false; // starts in slow‑mo (desiredSpeed != 1.0) private static volatile bool running = true; public static void Main() { var mode = new VideoMode(new Vector2u(1280, 720)); - var window = new RenderWindow(mode, "Pipe Resonator"); + var window = new RenderWindow(mode, "FluidSim - Helmholtz Resonator"); window.SetVerticalSyncEnabled(true); window.Closed += (_, _) => { running = false; window.Close(); }; window.MouseWheelScrolled += OnMouseWheel; @@ -39,9 +39,11 @@ public class Program soundEngine.Volume = 70; soundEngine.Start(); - //scenario = new PipeResonatorScenario(); + // Choose one scenario. The Helmholtz resonator is fully updated. //scenario = new HelmholtzResonatorScenario(); - scenario = new SodShockTubeScenario(); + //scenario = new PipeResonatorScenario(); // needs update to new API + //scenario = new SodShockTubeScenario(); // needs update to new API + scenario = new EngineScenario(); // also works (provided earlier) scenario.Initialize(SampleRate); @@ -51,9 +53,10 @@ public class Program double lastSpeedUpdateTime = stopwatch.Elapsed.TotalSeconds; // Resampling buffer - List simBuffer = new List(4096); + var simBuffer = new List(4096); double readIndex = 0.0; + // Prime the buffer with a few samples for (int i = 0; i < 4; i++) simBuffer.Add(scenario.Process()); @@ -73,7 +76,6 @@ public class Program lastSpeedUpdateTime = currentRealTime; // Smoothly transition currentSpeed → desiredSpeed - // When toggling, desiredSpeed jumps, but currentSpeed follows with a smooth lerp double smoothingRate = 8.0; // higher = faster catch‑up currentSpeed += (desiredSpeed - currentSpeed) * (1.0 - Math.Exp(-smoothingRate * dtSpeed)); @@ -122,7 +124,7 @@ public class Program break; } - // ---------- Drawing & title ---------- + // ---------- Drawing & window title ---------- if (currentRealTime - lastDrawTime >= drawInterval) { double actualSpeed = totalOutputSamples / (currentRealTime * SampleRate); @@ -156,7 +158,7 @@ public class Program desiredSpeed = Math.Clamp(desiredSpeed, MinSpeed, MaxSpeed); - // Update the remembered slow-mo speed (unless we are exactly at 1.0) + // Update the remembered slow‑mo speed (unless we are exactly at 1.0) if (!wasRealTime || Math.Abs(desiredSpeed - 1.0) > 1e-6) lastDesiredSpeed = desiredSpeed; diff --git a/Scenarios/EngineScenario.cs b/Scenarios/EngineScenario.cs new file mode 100644 index 0000000..c29fc71 --- /dev/null +++ b/Scenarios/EngineScenario.cs @@ -0,0 +1,238 @@ +using System; +using FluidSim.Components; +using SFML.Graphics; +using SFML.System; + +namespace FluidSim.Core +{ + public class EngineScenario : Scenario + { + private Solver solver; + private Volume0D cylinder; + private Pipe1D exhaustPipe; + private PipeVolumeConnection coupling; + private SoundProcessor soundProcessor; + + private double dt; + private double ambientPressure = 101325.0; + private double time; + + // Crankshaft + private double crankAngle = 0.0; + private const double TargetRPM = 4000.0; + private double angularVelocity; + + // Combustion + private const double CombustionPressure = 8.0 * 101325.0; + private const double CombustionTemperature = 1800.0; + + // 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; + + private double maxOrificeArea; + + // 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 + private int stepCount = 0; + private const int LogStepInterval = 1000; + private int combustionCount = 0; + private int misfireCount = 0; + + public override void Initialize(int sampleRate) + { + 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) + { + Gamma = 1.4, + GasConstant = 287.0 + }; + + // Exhaust pipe: length 2.5 m, radius 2 cm + double pipeLength = 2.5; + double pipeRadius = 0.02; + double pipeArea = Math.PI * pipeRadius * pipeRadius; + maxOrificeArea = pipeArea; + exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 70); + exhaustPipe.SetUniformState(1.225, 0.0, ambientPressure); + + // Coupling (valve starts closed) + coupling = new PipeVolumeConnection(cylinder, exhaustPipe, true, orificeArea: 0.0); + + solver = new Solver(); + solver.SetTimeStep(dt); + solver.AddVolume(cylinder); + solver.AddPipe(exhaustPipe); + solver.AddConnection(coupling); + // Use ZeroPressureOpen for strong reflections + solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.ZeroPressureOpen, ambientPressure); + + // Sound processor (tuned to pipe length) + 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.SetAmbientPressure(ambientPressure); + + lastFilteredPressure = ambientPressure; + + Console.WriteLine("=== EngineScenario (ZeroPressureOpen, boosted gains) ==="); + Console.WriteLine($"Target RPM: {TargetRPM}, Misfire prob: {MisfireProbability * 100: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("---------------------------------------------------------"); + } + + private double ValveOpenRatio(double crankRad) + { + double cycleAngle = crankRad % (4.0 * Math.PI); + double openStart = ValveOpenStart; + double openEnd = ValveOpenEnd; + + if (cycleAngle < openStart || cycleAngle > openEnd) + return 0.0; + + double fullOpenWindow = openEnd - openStart; + double closedWindow = 2.0 * ValveRampWidth; + if (fullOpenWindow <= closedWindow) + return 1.0; + + double tmid = (openStart + openEnd) / 2.0; + double dist = Math.Abs(cycleAngle - tmid); + double rampHalf = (fullOpenWindow - closedWindow) / 2.0; + if (dist <= rampHalf) + return 1.0; + 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; + } + } + + public override float Process() + { + // Update crank angle + crankAngle += angularVelocity * dt; + if (crankAngle >= 2.0 * Math.PI) + { + crankAngle -= 2.0 * Math.PI; + isMisfiring = rand.NextDouble() < MisfireProbability; + } + + // Power stroke at TDC + if (crankAngle < angularVelocity * dt && crankAngle >= 0.0) + { + 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; + 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; + combustionCount++; + } + } + + // Update valve area + double valveOpen = ValveOpenRatio(crankAngle); + coupling.OrificeArea = maxOrificeArea * valveOpen; + + 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; + + float audioSample = soundProcessor.Process(massFlow, (float)filteredPressure); + time += dt; + stepCount++; + + // Logging + if (stepCount % LogStepInterval == 0 || (crankAngle < angularVelocity * dt * 2 && !isMisfiring && combustionCount > 0)) + { + Console.WriteLine($"{time,7:F3} {crankAngle * 180.0 / Math.PI,6:F1} " + + $"{valveOpen * 100,6:F1} {massFlow,10:F4} " + + $"{combustionCount,3} {(isMisfiring ? "X" : "")}"); + } + + return audioSample; + } + + public override void Draw(RenderWindow target) + { + float winW = target.GetView().Size.X; + float winH = target.GetView().Size.Y; + float centerY = winH / 2f; + + 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); + target.Draw(cylRect); + + int n = exhaustPipe.GetCellCount(); + float pipeStartX = 120f, pipeEndX = winW - 60f; + float pipeLen = pipeEndX - pipeStartX; + float dx = pipeLen / (n - 1); + float baseRadius = 20f; + var vertices = new Vertex[n * 2]; + 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); + 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))); + var col = new Color(rCol, gCol, bCol); + + 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 index 293b86b..af8ded5 100644 --- a/Scenarios/HelmholtzResonatorScenario.cs +++ b/Scenarios/HelmholtzResonatorScenario.cs @@ -1,6 +1,5 @@ using System; using FluidSim.Components; -using FluidSim.Interfaces; using FluidSim.Utils; using SFML.Graphics; using SFML.System; @@ -12,7 +11,7 @@ namespace FluidSim.Core private Solver solver; private Volume0D cavity; private Pipe1D neck; - private Connection coupling; + private PipeVolumeConnection coupling; private int stepCount; private double time; private double dt; @@ -38,12 +37,8 @@ namespace FluidSim.Core neck = new Pipe1D(neckLength, neckArea, sampleRate, forcedCellCount: 40); neck.SetUniformState(1.225, 0.0, ambientPressure); - coupling = new Connection(neck.PortA, cavity.Port) - { - Area = neckArea, - DischargeCoefficient = 0.62, - Gamma = 1.4 - }; + // 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); @@ -51,8 +46,8 @@ namespace FluidSim.Core solver.AddPipe(neck); solver.AddConnection(coupling); - // Port A (left) = volume coupling, Port B (right) = open end - solver.SetPipeBoundary(neck, isA: true, BoundaryType.VolumeCoupling); + // 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); } @@ -68,11 +63,11 @@ namespace FluidSim.Core if (stepCount % 20 == 0) { double pCav = cavity.Pressure; - double mdotA = neck.PortA.MassFlowRate; // positive = into pipe (leaving cavity) + // 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, " + - $"mdot_A={mdotA * 1e3:F4} g/s, audio={audio:F4}"); + $"audio={audio:F4}"); } return audio; @@ -100,7 +95,7 @@ namespace FluidSim.Core float dx = neckLenPx / (n - 1); float baseRadius = 20f; - Vertex[] vertices = new Vertex[n * 2]; + var vertices = new Vertex[n * 2]; for (int i = 0; i < n; i++) { float x = neckStartX + i * dx;