diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index bd0ba4d..39171ef 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -4,10 +4,8 @@ using FluidSim.Interfaces; namespace FluidSim.Components { /// - /// 1‑D compressible Euler pipe (finite‑volume, HLLC flux). - /// Boundary conditions are set externally via SetGhostLeft/Right. - /// Enforces that ghosts are always valid before stepping. - /// Uses exponential damping and Newtonian energy relaxation. + /// 1‑D compressible Euler pipe with Lax‑Friedrichs finite‑volume scheme. + /// Ghost states are set externally via SetGhostLeft/Right; they are always required. /// public class Pipe1D : IComponent { @@ -17,7 +15,6 @@ namespace FluidSim.Components public double DampingMultiplier { get; set; } = 1.0; public double EnergyRelaxationRate { get; set; } = 0.0; // 1/s - // Ambient pressure for the energy relaxation term (default 101325 Pa) private double _ambientPressure = 101325.0; public double AmbientPressure { @@ -29,37 +26,21 @@ namespace FluidSim.Components } } - // Geometry - private readonly int _n; // number of real cells - private readonly double _dx; // cell size (m) - private readonly double _diameter; // m + private readonly int _n; + private readonly double _dx; + private readonly double _diameter; private readonly double _gamma = 1.4; - // Conserved variables [0 .. _n-1] - private double[] _rho; - private double[] _rhou; - private double[] _E; + private double[] _rho, _rhou, _E; + private double[] _fluxM, _fluxP, _fluxE; // flux at cell faces (0.._n) - // Face fluxes [0 .. _n] - private double[] _fluxM; - private double[] _fluxP; - private double[] _fluxE; - - // Ghost cells (set externally) private double _rhoGhostL, _uGhostL, _pGhostL; private double _rhoGhostR, _uGhostR, _pGhostR; private bool _ghostLValid, _ghostRValid; - // Pre‑computed damping coefficient private double _laminarCoeff; - private double _ambientEnergyReference; // internal energy density at ambient pressure + private double _ambientEnergyReference; - /// - /// Initialise a pipe with a given cell count. - /// - /// Pipe length (m). - /// Cross‑sectional area (m²). - /// Number of finite‑volume cells (≥ 4). public Pipe1D(double length, double area, int cellCount) { if (cellCount < 4) throw new ArgumentException("cellCount must be at least 4"); @@ -77,48 +58,32 @@ namespace FluidSim.Components _fluxP = new double[_n + 1]; _fluxE = new double[_n + 1]; - // Laminar damping coefficient for air at 20°C (multiplied by DampingMultiplier each step) double mu_air = 1.8e-5; double radius = _diameter * 0.5; _laminarCoeff = 8.0 * mu_air / (radius * radius); - // Ambient energy reference (internal energy per unit volume at 101325 Pa) _ambientEnergyReference = 101325.0 / (_gamma - 1.0); PortA = new Port { Owner = this }; PortB = new Port { Owner = this }; - // Initial state = still air at ambient conditions SetUniformState(1.225, 0.0, 101325.0); } IReadOnlyList IComponent.Ports => new[] { PortA, PortB }; - // No integration needed for pipes – state is advanced via sub‑steps public void UpdateState(double dt) { } - // ---------- Ghost cell interface ---------- + // ---------- Ghost interface ---------- public void SetGhostLeft(double rho, double u, double p) { - _rhoGhostL = rho; - _uGhostL = u; - _pGhostL = p; - _ghostLValid = true; + _rhoGhostL = rho; _uGhostL = u; _pGhostL = p; _ghostLValid = true; } - public void SetGhostRight(double rho, double u, double p) { - _rhoGhostR = rho; - _uGhostR = u; - _pGhostR = p; - _ghostRValid = true; - } - - public void ClearGhostFlags() - { - _ghostLValid = false; - _ghostRValid = false; + _rhoGhostR = rho; _uGhostR = u; _pGhostR = p; _ghostRValid = true; } + public void ClearGhostFlags() { _ghostLValid = false; _ghostRValid = false; } public (double rho, double u, double p) GetInteriorStateLeft() { @@ -127,7 +92,6 @@ namespace FluidSim.Components double p = PressureScalar(0); return (rho, u, p); } - public (double rho, double u, double p) GetInteriorStateRight() { double rho = Math.Max(_rho[_n - 1], 1e-12); @@ -135,18 +99,11 @@ namespace FluidSim.Components double p = PressureScalar(_n - 1); return (rho, u, p); } - public int CellCount => _n; - public double GetCellDensity(int i) => _rho[i]; - public double GetCellVelocity(int i) - { - double rho = Math.Max(_rho[i], 1e-12); - return _rhou[i] / rho; - } + public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12); public double GetCellPressure(int i) => PressureScalar(i); - // ---------- Sub‑stepping ---------- public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8) { double maxW = 0.0; @@ -163,38 +120,65 @@ namespace FluidSim.Components return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx))); } - // ---------- Main simulation step (per sub‑step) ---------- + // ---------- Main step (per sub‑step) ---------- public void SimulateSingleStep(double dtSub) { - // Enforce that both ends have been provided with ghost states if (!_ghostLValid || !_ghostRValid) - throw new InvalidOperationException("Pipe boundary ghosts not set before SimulateSingleStep."); + throw new InvalidOperationException("Ghost cells not set before SimulateSingleStep."); double dt = dtSub; int n = _n; - // Left boundary face (index 0) - HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, _rho[0], _rhou[0] / _rho[0], PressureScalar(0), - out _fluxM[0], out _fluxP[0], out _fluxE[0]); + // ---- Compute fluxes at all faces using Lax‑Friedrichs ---- + // Left face (0): between ghostL and cell 0 + double rL = Math.Max(_rhoGhostL, 1e-12); + double pL = _pGhostL; + double uL = _uGhostL; + double eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL; + + double rR = Math.Max(_rho[0], 1e-12); + double pR = PressureScalar(0); + double uR = _rhou[0] / rR; + double eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR; + + LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR, + out _fluxM[0], out _fluxP[0], out _fluxE[0]); // Internal faces (1 .. n-1) for (int f = 1; f < n; f++) { - double rhoL = Math.Max(_rho[f - 1], 1e-12); - double uL = _rhou[f - 1] / rhoL; - double pL = PressureScalar(f - 1); - double rhoR = Math.Max(_rho[f], 1e-12); - double uR = _rhou[f] / rhoR; - double pR = PressureScalar(f); - HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out _fluxM[f], out _fluxP[f], out _fluxE[f]); + int iL = f - 1; + int iR = f; + + rL = Math.Max(_rho[iL], 1e-12); + pL = PressureScalar(iL); + uL = _rhou[iL] / rL; + eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL; + + rR = Math.Max(_rho[iR], 1e-12); + pR = PressureScalar(iR); + uR = _rhou[iR] / rR; + eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR; + + LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR, + out _fluxM[f], out _fluxP[f], out _fluxE[f]); } - // Right boundary face (index n) - HLLCFlux(_rho[_n - 1], _rhou[_n - 1] / _rho[_n - 1], PressureScalar(_n - 1), - _rhoGhostR, _uGhostR, _pGhostR, - out _fluxM[n], out _fluxP[n], out _fluxE[n]); + // Right face (n): between cell n-1 and ghostR + rL = Math.Max(_rho[n - 1], 1e-12); + pL = PressureScalar(n - 1); + uL = _rhou[n - 1] / rL; + eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL; - // Cell update + rR = Math.Max(_rhoGhostR, 1e-12); + pR = _pGhostR; + uR = _uGhostR; + eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR; + + LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR, + out _fluxM[n], out _fluxP[n], out _fluxE[n]); + + // ---- Cell update ---- double dt_dx = dt / _dx; double coeff = _laminarCoeff * DampingMultiplier; double relaxRate = EnergyRelaxationRate; @@ -213,15 +197,12 @@ namespace FluidSim.Components double newRu = ru - dt_dx * dP; double newE = E - dt_dx * dE_flux; - // Wall friction damping (laminar) double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt); newRu *= dampingFactor; - // Newtonian cooling toward ambient energy double relaxFactor = Math.Exp(-relaxRate * dt); newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor; - // Clamps – minimum density 1e-12, minimum pressure 100 Pa newR = Math.Max(newR, 1e-12); double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12); double eMin = 100.0 / (_gamma - 1.0) + kin; @@ -232,77 +213,54 @@ namespace FluidSim.Components _E[i] = newE; } - // Update port states to reflect the current interior state (for audio / monitoring) + // Update port states (double rhoA, double uA, double pA) = GetInteriorStateLeft(); - PortA.Pressure = pA; - PortA.Density = rhoA; + PortA.Pressure = pA; PortA.Density = rhoA; PortA.Temperature = pA / (rhoA * 287.0); PortA.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pA / rhoA; (double rhoB, double uB, double pB) = GetInteriorStateRight(); - PortB.Pressure = pB; - PortB.Density = rhoB; + PortB.Pressure = pB; PortB.Density = rhoB; PortB.Temperature = pB / (rhoB * 287.0); PortB.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pB / rhoB; } - // ---------- Private helpers ---------- + // ---------- Lax‑Friedrichs flux ---------- + private void LaxFriedrichsFlux(double rL, double uL, double pL, double eL, + double rR, double uR, double pR, double eR, + out double fm, out double fp, out double fe) + { + // Primitive states + double rhoL = rL, rhoR = rR; + double EL = rhoL * eL; // total energy per volume = rho * (specific total energy) + double ER = rhoR * eR; + + // Conserved vectors U = (ρ, ρu, E) + // Flux F = (ρu, ρu²+p, (E+p)u) + double Fm_L = rhoL * uL; + double Fp_L = rhoL * uL * uL + pL; + double Fe_L = (EL + pL) * uL; + + double Fm_R = rhoR * uR; + double Fp_R = rhoR * uR * uR + pR; + double Fe_R = (ER + pR) * uR; + + // Lax‑Friedrichs dissipation coefficient α = max(|u|+c) over whole domain, but here we use local max to be simple: + double cL = Math.Sqrt(_gamma * pL / rL); + double cR = Math.Sqrt(_gamma * pR / rR); + double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR); + + fm = 0.5 * (Fm_L + Fm_R) - 0.5 * alpha * (rhoR - rhoL); + fp = 0.5 * (Fp_L + Fp_R) - 0.5 * alpha * (rhoR * uR - rhoL * uL); + fe = 0.5 * (Fe_L + Fe_R) - 0.5 * alpha * (ER - EL); + } + private double PressureScalar(int i) { double rho = Math.Max(_rho[i], 1e-12); return (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho); } - /// - /// HLLC approximate Riemann solver (Toro, 1997). - /// Computes the numerical flux at a face given left and right states. - /// - private void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR, - out double fm, out double fp, out double fe) - { - double cL = Math.Sqrt(_gamma * pL / rL); - double cR = Math.Sqrt(_gamma * pR / rR); - double EL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL; // specific total energy - double ER = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR; - - // Wave speed estimates (Davis, 1988) - double SL = Math.Min(uL - cL, uR - cR); - double SR = Math.Max(uL + cL, uR + cR); - - double denom = rL * (SL - uL) - rR * (SR - uR); - double Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom; - - double Fm_L = rL * uL; - double Fp_L = rL * uL * uL + pL; - double Fe_L = (rL * EL + pL) * uL; - - double Fm_R = rR * uR; - double Fp_R = rR * uR * uR + pR; - double Fe_R = (rR * ER + pR) * uR; - - if (SL >= 0) { fm = Fm_L; fp = Fp_L; fe = Fe_L; } - else if (SR <= 0) { fm = Fm_R; fp = Fp_R; fe = Fe_R; } - else if (Ss >= 0) - { - double rsL = rL * (SL - uL) / (SL - Ss); - double ps = pL + rL * (SL - uL) * (Ss - uL); - double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL))); - fm = rsL * Ss; - fp = rsL * Ss * Ss + ps; - fe = (rsL * EsL + ps) * Ss; - } - else - { - double rsR = rR * (SR - uR) / (SR - Ss); - double ps = pR + rR * (SR - uR) * (Ss - uR); - double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR))); - fm = rsR * Ss; - fp = rsR * Ss * Ss + ps; - fe = (rsR * EsR + ps) * Ss; - } - } - - /// Initialise all cells to a uniform state (rho, u, p). public void SetUniformState(double rho, double u, double p) { double e = p / ((_gamma - 1.0) * rho); @@ -314,5 +272,24 @@ namespace FluidSim.Components _E[i] = E; } } + + public void SetCellState(int i, double rho, double u, double p) + { + if (i < 0 || i >= _n) return; + double e = p / ((_gamma - 1.0) * rho); + double E = rho * e + 0.5 * rho * u * u; + _rho[i] = rho; + _rhou[i] = rho * u; + _E[i] = E; + } + + public void SetCellPressure(int i, double p) + { + if (i < 0 || i >= _n) return; + double rho = _rho[i]; + double u = _rhou[i] / rho; + double e = p / ((_gamma - 1.0) * rho); + _E[i] = rho * e + 0.5 * rho * u * u; + } } } \ No newline at end of file diff --git a/Core/IsentropicOrifice.cs b/Core/IsentropicOrifice.cs index effe861..ce23e83 100644 --- a/Core/IsentropicOrifice.cs +++ b/Core/IsentropicOrifice.cs @@ -4,38 +4,38 @@ namespace FluidSim.Core { /// /// Compressible flow through an orifice, modelled as an isentropic nozzle. - /// Supports choked and unchoked flow, forward and reverse. + /// The caller provides the upstream stagnation state (pUp, rhoUp, TUp), + /// downstream pressure, orifice area, discharge coefficient, and gas properties. + /// Returns the face state and mass flow from upstream to downstream. /// public static class IsentropicOrifice { - /// - /// Compute mass flow and face primitive state for an orifice. - /// - /// Upstream stagnation pressure (Pa). - /// Upstream stagnation density (kg/m³). - /// Ratio of specific heats. - /// Specific gas constant (J/kg·K). - /// Downstream static pressure (Pa). - /// Effective orifice area (m²). - /// Discharge coefficient (default 0.62). - /// Mass flow rate (kg/s), positive from upstream to downstream. - /// Face density (kg/m³). - /// Face velocity (m/s). - /// Face pressure (Pa). - public static void Compute(double pUp, double rhoUp, double TUp, double gamma, double R, - double pDown, double area, double Cd, - out double mdot, out double rhoFace, out double uFace, out double pFace) + public static void Compute( + double pUp, double rhoUp, double TUp, // upstream stagnation + double pDown, // downstream back pressure + double gamma, double R, double area, double Cd, + out double mdot, out double rhoFace, out double uFace, out double pFace) { - // mdot is positive from upstream to downstream. - double pr = Math.Max(pDown / pUp, 1e-6); - double prCrit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0)); - if (pr < prCrit) pr = prCrit; + mdot = 0; rhoFace = rhoUp; uFace = 0; pFace = pUp; - double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0)); - uFace = M * Math.Sqrt(gamma * R * TUp); + if (area <= 0 || pUp <= 0 || rhoUp <= 0 || TUp <= 0) + return; + + double pr = pDown / pUp; + if (pr < 1e-6) pr = 1e-6; + + double prCrit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0)); + if (pr < prCrit) pr = prCrit; // choked flow + + double exponent = (gamma - 1.0) / gamma; + double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -exponent) - 1.0)); + if (double.IsNaN(M)) M = 0; + + double aUp = Math.Sqrt(gamma * R * TUp); + uFace = M * aUp; rhoFace = rhoUp * Math.Pow(pr, 1.0 / gamma); pFace = pUp * pr; - mdot = rhoFace * uFace * area * Cd; // mass flow from upstream to downstream + mdot = rhoFace * uFace * area * Cd; // positive from upstream to downstream } } } \ No newline at end of file diff --git a/Core/OpenEndLink.cs b/Core/OpenEndLink.cs index 939c1fe..a6efce3 100644 --- a/Core/OpenEndLink.cs +++ b/Core/OpenEndLink.cs @@ -3,19 +3,15 @@ using FluidSim.Components; namespace FluidSim.Core { - /// - /// Characteristic open‑end boundary condition. - /// For subsonic outflow the outgoing Riemann invariant is conserved, - /// and the ghost pressure is set to the prescribed ambient value. - /// public class OpenEndLink { public Pipe1D Pipe { get; } public bool IsLeftEnd { get; } public double AmbientPressure { get; set; } = 101325.0; public double Gamma { get; set; } = 1.4; + public double GasConstant { get; set; } = 287.0; + public double AmbientTemperature { get; set; } = 300.0; - // Last resolved state (for audio / monitoring) public double LastMassFlowRate { get; private set; } public double LastFaceDensity { get; private set; } public double LastFaceVelocity { get; private set; } @@ -27,9 +23,6 @@ namespace FluidSim.Core IsLeftEnd = isLeftEnd; } - /// - /// Compute the ghost state and mass flow for one sub‑step. - /// public void Resolve(double dtSub) { (double rhoInt, double uInt, double pInt) = IsLeftEnd @@ -40,80 +33,61 @@ namespace FluidSim.Core double gm1 = gamma - 1.0; double cInt = Math.Sqrt(gamma * pInt / Math.Max(rhoInt, 1e-12)); double pAmb = AmbientPressure; + double rhoAmb = pAmb / (GasConstant * AmbientTemperature); + double aAmb = Math.Sqrt(gamma * pAmb / rhoAmb); double rhoGhost, uGhost, pGhost; - double mdot; - if (IsLeftEnd) + // ----- Supersonic outflow: extrapolate interior ----- + bool supersonicOut = IsLeftEnd ? (uInt <= -cInt) : (uInt >= cInt); + if (supersonicOut) { - // Left end: outgoing invariant is J- = u - 2c/(γ-1) + rhoGhost = rhoInt; + uGhost = uInt; + pGhost = pInt; + } + else + { + // Riemann invariants + double J_plus = uInt + 2.0 * cInt / gm1; double J_minus = uInt - 2.0 * cInt / gm1; - if (uInt <= -cInt) // supersonic inflow (all info from outside) - { - // Simple reservoir model – use ambient density and temperature 300 K - rhoGhost = pAmb / (287.0 * 300.0); - uGhost = uInt; // keep interior velocity (should be supersonic inward) - pGhost = pAmb; - } - else if (uInt < 0) // subsonic inflow - { - double rhoAmb = pAmb / (287.0 * 300.0); - double cAmb = Math.Sqrt(gamma * pAmb / rhoAmb); - uGhost = J_minus + 2.0 * cAmb / gm1; - rhoGhost = rhoAmb; - pGhost = pAmb; - } - else // subsonic outflow (uInt >= 0) - { - double s = pInt / Math.Pow(rhoInt, gamma); - rhoGhost = Math.Pow(pAmb / s, 1.0 / gamma); - double cGhost = Math.Sqrt(gamma * pAmb / rhoGhost); - uGhost = J_minus + 2.0 * cGhost / gm1; - if (uGhost < 0) uGhost = 0; - pGhost = pAmb; - } - } - else // Right end - { - // Right end: outgoing invariant is J+ = u + 2c/(γ-1) - double J_plus = uInt + 2.0 * cInt / gm1; + // Trial subsonic outflow ghost state + double s = pInt / Math.Pow(rhoInt, gamma); + double rhoOut = Math.Pow(pAmb / s, 1.0 / gamma); + double cOut = Math.Sqrt(gamma * pAmb / rhoOut); + double uOut = IsLeftEnd + ? (J_minus + 2.0 * cOut / gm1) + : (J_plus - 2.0 * cOut / gm1); - if (uInt >= cInt) // supersonic outflow + bool outflowPossible = IsLeftEnd ? (uOut <= 0) : (uOut >= 0); + + if (outflowPossible) { - rhoGhost = rhoInt; - uGhost = uInt; - pGhost = pInt; - } - else if (uInt >= 0) // subsonic outflow - { - double s = pInt / Math.Pow(rhoInt, gamma); - rhoGhost = Math.Pow(pAmb / s, 1.0 / gamma); - double cGhost = Math.Sqrt(gamma * pAmb / rhoGhost); - uGhost = J_plus - 2.0 * cGhost / gm1; - if (uGhost < 0) uGhost = 0; + // Subsonic outflow pGhost = pAmb; + rhoGhost = rhoOut; + uGhost = uOut; } - else // subsonic inflow (uInt < 0) + else { - double rhoAmb = pAmb / (287.0 * 300.0); - double cAmb = Math.Sqrt(gamma * pAmb / rhoAmb); - uGhost = J_plus - 2.0 * cAmb / gm1; + // Subsonic inflow (ambient reservoir model) + pGhost = pAmb; rhoGhost = rhoAmb; - pGhost = pAmb; + uGhost = IsLeftEnd + ? (J_minus + 2.0 * aAmb / gm1) + : (J_plus - 2.0 * aAmb / gm1); } } - // Apply ghost to pipe if (IsLeftEnd) Pipe.SetGhostLeft(rhoGhost, uGhost, pGhost); else Pipe.SetGhostRight(rhoGhost, uGhost, pGhost); - // Mass flow (positive = out of pipe) double area = Pipe.Area; - mdot = rhoGhost * uGhost * area; - if (IsLeftEnd) mdot = -mdot; // positive u into pipe, so out of pipe is negative u + double mdot = rhoGhost * uGhost * area; + if (IsLeftEnd) mdot = -mdot; LastMassFlowRate = mdot; LastFaceDensity = rhoGhost; LastFaceVelocity = uGhost; diff --git a/Core/OrificeLink.cs b/Core/OrificeLink.cs index 156590b..2843589 100644 --- a/Core/OrificeLink.cs +++ b/Core/OrificeLink.cs @@ -6,7 +6,8 @@ namespace FluidSim.Core { /// /// Connects a port (volume or atmosphere) to one end of a pipe via an orifice. - /// The area can be dynamic (Func). + /// Uses the isentropic nozzle model for the steady‑state relationship, + /// and includes acoustic inertance for dynamic (Helmholtz) behaviour. /// public class OrificeLink { @@ -15,105 +16,131 @@ namespace FluidSim.Core public bool IsPipeLeftEnd { get; } public Func AreaProvider { get; set; } public double DischargeCoefficient { get; set; } = 0.62; - public double Gamma { get; set; } = 1.4; - public double GasConstant { get; set; } = 287.0; - // Last resolved state (for audio/monitoring) + // Acoustic length (wall thickness + end correction) – controls the resonance frequency + public double EffectiveLength { get; set; } = 0.001; // 1 mm + + // Whether to include inertance; if false, uses the steady‑state nozzle model directly + public bool UseInertance { get; set; } = true; + + // Current mass flow (kg/s, positive = volume → pipe) + private double _mdot; + public double LastMassFlowRate { get; private set; } public double LastFaceDensity { get; private set; } public double LastFaceVelocity { get; private set; } public double LastFacePressure { get; private set; } - public OrificeLink(Port volumePort, Pipe1D pipe, bool isPipeLeftEnd, Func areaProvider) + public OrificeLink(Port? volumePort, Pipe1D pipe, bool isPipeLeftEnd, Func areaProvider) { - VolumePort = volumePort ?? throw new ArgumentNullException(nameof(volumePort)); + VolumePort = volumePort; // null is allowed Pipe = pipe ?? throw new ArgumentNullException(nameof(pipe)); IsPipeLeftEnd = isPipeLeftEnd; AreaProvider = areaProvider ?? throw new ArgumentNullException(nameof(areaProvider)); + _mdot = 0.0; } - /// - /// Resolve the coupling for one sub‑step. Computes nozzle flow (isentropic) - /// and sets the pipe ghost cell and the port flow rates. - /// public void Resolve(double dtSub) { double area = AreaProvider(); - if (area < 1e-12) + // Closed wall or missing volume port => reflective boundary + if (area < 1e-12 || VolumePort == null) { SetClosedWall(); return; } - // Retrieve volume state - double volP = VolumePort.Pressure; + // Gather volume state + double volP = VolumePort.Pressure; double volRho = VolumePort.Density; - double volT = VolumePort.Temperature; - double volH = VolumePort.SpecificEnthalpy; + double volT = VolumePort.Temperature; + double volH = VolumePort.SpecificEnthalpy; - // Retrieve pipe interior state at the connected end + // Gather pipe interior state at the connected end (double pipeRho, double pipeU, double pipeP) = IsPipeLeftEnd ? Pipe.GetInteriorStateLeft() : Pipe.GetInteriorStateRight(); - // Determine upstream/downstream: if volume pressure > pipe pressure, flow is out of volume (negative into volume). - bool flowOutOfVolume = volP > pipeP; - double pUp, rhoUp, TUp, pDown; - if (flowOutOfVolume) + double pipeT = pipeP / Math.Max(pipeRho * 287.0, 1e-12); + double gamma = 1.4; + double R = 287.0; + + // ---- Steady‑state mass flow from isentropic nozzle ---- + double mdotSS; // positive = volume → pipe + double rhoFace, uFace, pFace; + + if (volP >= pipeP) { - pUp = volP; rhoUp = volRho; TUp = volT; pDown = pipeP; + IsentropicOrifice.Compute(volP, volRho, volT, pipeP, gamma, R, area, DischargeCoefficient, + out double mdotUpToDown, out rhoFace, out uFace, out pFace); + mdotSS = mdotUpToDown; // volume → pipe } else { - // Pipe is upstream - pUp = pipeP; rhoUp = pipeRho; TUp = pipeP / (pipeRho * GasConstant); // temperature from pipe - pDown = volP; + IsentropicOrifice.Compute(pipeP, pipeRho, pipeT, volP, gamma, R, area, DischargeCoefficient, + out double mdotUpToDown, out rhoFace, out uFace, out pFace); + mdotSS = -mdotUpToDown; // pipe → volume → negative for volume→pipe convention } - // Compute isentropic nozzle flow - IsentropicOrifice.Compute(pUp, rhoUp, TUp, Gamma, GasConstant, pDown, area, DischargeCoefficient, - out double mdotUpstreamToDown, out double rhoFace, out double uFace, out double pFace); - - // mdotUpstreamToDown is positive from upstream to downstream. - // Convert to mass flow into volume (positive mdot = into volume). - double mdotVolume; - if (flowOutOfVolume) - mdotVolume = -mdotUpstreamToDown; // out of volume is negative + // ---- Inertance ODE (optional) ---- + if (UseInertance) + { + double rhoUp = _mdot >= 0 ? volRho : pipeRho; + double inertance = rhoUp * EffectiveLength / area; + double dp = volP - pipeP; + double resistance = Math.Abs(dp) / Math.Max(Math.Abs(mdotSS), 1e-12); + double dmdot_dt = (dp - resistance * _mdot) / inertance; + _mdot += dmdot_dt * dtSub; + } else - mdotVolume = mdotUpstreamToDown; // into volume is positive + { + _mdot = mdotSS; + } - // Clamp mass flow to available mass in volume (if it is a Volume0D) + // Clamp outflow to available mass (if finite volume) if (VolumePort.Owner is Volume0D vol) { - double maxMdot = vol.Mass / dtSub; - if (mdotVolume > maxMdot) mdotVolume = maxMdot; - if (mdotVolume < -maxMdot) mdotVolume = -maxMdot; + double maxOut = vol.Mass / dtSub; + if (_mdot > maxOut) _mdot = maxOut; } - // Apply ghost state to pipe + // ---- Ghost state ---- + // Density = upstream density (consistent with current flow direction) + rhoFace = _mdot >= 0 ? volRho : pipeRho; + // Pressure = downstream pressure (consistent with nozzle exit) + pFace = _mdot >= 0 ? pipeP : volP; + // Velocity magnitude derived from actual mass flow + double mdotMag = Math.Abs(_mdot); + uFace = mdotMag / (rhoFace * area); + + if (IsPipeLeftEnd) + uFace = _mdot >= 0 ? uFace : -uFace; // left end: positive u = into pipe + else + uFace = _mdot >= 0 ? -uFace : uFace; // right end: positive u = out of pipe + + // Apply ghost to pipe if (IsPipeLeftEnd) Pipe.SetGhostLeft(rhoFace, uFace, pFace); else Pipe.SetGhostRight(rhoFace, uFace, pFace); - // Store results - LastMassFlowRate = mdotVolume; - LastFaceDensity = rhoFace; + // ---- Store results ---- + double mdotIntoVolume = -_mdot; // positive = into volume + LastMassFlowRate = mdotIntoVolume; + LastFaceDensity = rhoFace; LastFaceVelocity = uFace; LastFacePressure = pFace; - // Set port flow rates for volume integration - VolumePort.MassFlowRate = mdotVolume; - if (mdotVolume >= 0) + VolumePort.MassFlowRate = mdotIntoVolume; + + // Enthalpy for volume integration + if (mdotIntoVolume >= 0) // inflow → pipe enthalpy { - // Inflow: enthalpy comes from upstream (pipe) - double pPipe = pipeP; - double rhoPipe = pipeRho; - VolumePort.SpecificEnthalpy = Gamma / (Gamma - 1.0) * pPipe / rhoPipe; + double hPipe = gamma / (gamma - 1.0) * pipeP / Math.Max(pipeRho, 1e-12); + VolumePort.SpecificEnthalpy = hPipe; } - else + else // outflow → volume's own enthalpy { - // Outflow: volume's own specific enthalpy VolumePort.SpecificEnthalpy = volH; } } @@ -130,11 +157,12 @@ namespace FluidSim.Core Pipe.SetGhostRight(rInt, -uInt, pInt); LastMassFlowRate = 0.0; - LastFaceDensity = rInt; + LastFaceDensity = rInt; LastFaceVelocity = 0.0; LastFacePressure = pInt; - VolumePort.MassFlowRate = 0.0; - // Keep specific enthalpy as is (not used) + // Don't touch VolumePort if it's null + if (VolumePort != null) + VolumePort.MassFlowRate = 0.0; } } } \ No newline at end of file diff --git a/Core/SoundProcessor.cs b/Core/SoundProcessor.cs index 9229e9f..59d9fba 100644 --- a/Core/SoundProcessor.cs +++ b/Core/SoundProcessor.cs @@ -1,18 +1,20 @@ using System; -using FluidSim.Interfaces; +using FluidSim.Core; namespace FluidSim.Core { public class SoundProcessor { private readonly double dt; - private readonly double scaleFactor; // 1 / (4π r) and a user gain + private readonly double scaleFactor; // 1 / (4π r) private double prevMassFlowOut; - - // Simple low‑pass for derivative smoothing (≈ 2‑3 ms) private double smoothDMdt; private readonly double alpha; + // New: low‑pass the mass flow signal before derivative + private double flowLP; + private readonly double lpAlpha; + public float Gain { get; set; } = 1.0f; public SoundProcessor(int sampleRate, double listenerDistanceMeters = 1.0) @@ -20,29 +22,34 @@ namespace FluidSim.Core dt = 1.0 / sampleRate; scaleFactor = 1.0 / (4.0 * Math.PI * listenerDistanceMeters); - // Smoothing time constant ~ 2 ms, blocks single‑sample spikes - double tau = 0.002; + // Smoothing time constant for the derivative: 10 ms (much smoother) + double tau = 0.010; // 10 ms alpha = Math.Exp(-dt / tau); + + // Low‑pass time constant for the mass flow: 5 ms (kneecap high‑freq directly) + double tauLP = 0.005; + lpAlpha = Math.Exp(-dt / tauLP); } - public float Process(Port port) + public float Process(OpenEndLink openEnd) { - // Outflow mass flow (positive = leaving pipe) - double flowOut = -port.MassFlowRate; + double flowOut = openEnd.LastMassFlowRate; - // Derivative - double rawDerivative = (flowOut - prevMassFlowOut) / dt; - prevMassFlowOut = flowOut; + // Low‑pass the mass flow signal + flowLP = lpAlpha * flowLP + (1.0 - lpAlpha) * flowOut; - // Smooth the derivative to kill isolated spikes + // Derivative of the smoothed mass flow + double rawDerivative = (flowLP - prevMassFlowOut) / dt; + prevMassFlowOut = flowLP; + + // Smooth the derivative smoothDMdt = alpha * smoothDMdt + (1.0 - alpha) * rawDerivative; // Far‑field monopole pressure double pressure = smoothDMdt * scaleFactor * Gain; - // Soft clip to ±1 for audio output (safe limit) - float sample = (float)Math.Tanh(pressure); - return sample; + // Soft clip to ±1 (should rarely trigger now) + return (float)Math.Tanh(pressure); } } } \ No newline at end of file diff --git a/Program.cs b/Program.cs index 330b59b..f4cd113 100644 --- a/Program.cs +++ b/Program.cs @@ -14,7 +14,8 @@ public class Program private static Scenario scenario; // Speed control (existing + new throttle) - private static double desiredSpeed = 0.01; + private static double desiredSpeed = 0.001; + //private static double desiredSpeed = 1.0; private static double currentSpeed = desiredSpeed; private const double MinSpeed = 0.0001; private const double MaxSpeed = 1.0; diff --git a/Scenarios/TestScenario.cs b/Scenarios/TestScenario.cs index 6fead4d..6178091 100644 --- a/Scenarios/TestScenario.cs +++ b/Scenarios/TestScenario.cs @@ -9,40 +9,41 @@ namespace FluidSim.Tests public class TestScenario : Scenario { private Solver solver; - private Volume0D volume; private Pipe1D pipe; - private Atmosphere atmosphere; - private OrificeLink orificeLink; - private OpenEndLink openEndLink; + private OrificeLink closedEnd; // left end – closed wall + private OpenEndLink openEndLink; // right end – atmosphere + private SoundProcessor soundProcessor; + private OutdoorExhaustReverb reverb; private int stepCount; + private double simTime; // elapsed simulation time (seconds) + private double pulseInterval = 0.1; // seconds between pulses + private double nextPulseTime; + private double dt; public override void Initialize(int sampleRate) { - double dt = 1.0 / sampleRate; + dt = 1.0 / sampleRate; + + soundProcessor = new SoundProcessor(sampleRate, 1); + soundProcessor.Gain = 10; + + reverb = new OutdoorExhaustReverb(sampleRate); solver = new Solver(); solver.SetTimeStep(dt); - volume = new Volume0D(1e-3, 150000.0, 300.0); - solver.AddComponent(volume); - - pipe = new Pipe1D(2.0, 1e-4, 20); + // Pipe: 2 m long, 1 cm² area, 200 cells + pipe = new Pipe1D(length: 2, area: 1e-4, cellCount: 100); solver.AddComponent(pipe); - atmosphere = new Atmosphere(); - solver.AddComponent(atmosphere); + // Initially pipe at ambient conditions + pipe.SetUniformState(1.225, 0.0, 101325.0); - // Volume → left pipe end (orifice) - var volPort = volume.CreatePort(); - orificeLink = new OrificeLink(volPort, pipe, isPipeLeftEnd: true, areaProvider: () => 1e-5) - { - DischargeCoefficient = 0.62, - Gamma = volume.Gamma, - GasConstant = volume.GasConstant - }; - solver.AddOrificeLink(orificeLink); + // Left end: closed wall (area = 0 → reflective) + closedEnd = new OrificeLink(null, pipe, isPipeLeftEnd: true, areaProvider: () => 0.0); + solver.AddOrificeLink(closedEnd); - // Right pipe end → atmosphere (characteristic open‑end) + // Right end: open to atmosphere openEndLink = new OpenEndLink(pipe, isLeftEnd: false) { AmbientPressure = 101325.0, @@ -51,45 +52,82 @@ namespace FluidSim.Tests solver.AddOpenEndLink(openEndLink); stepCount = 0; - Console.WriteLine("TestScenario initialized with sampleRate = " + sampleRate); + simTime = 0.0; + nextPulseTime = pulseInterval; // first pulse at 100 ms + + Console.WriteLine("Pulse reflection test – closed left, open right"); + Console.WriteLine("Pulse injected every 100 ms at left end (cell 0)"); } public override float Process() { solver.Step(); stepCount++; + simTime += dt; - if (stepCount % 100 == 0) + // ---- Inject a pressure pulse at the closed end every 100 ms ---- + if (simTime >= nextPulseTime) { - double volPressure = volume.Pressure; - double volMass = volume.Mass; - double pipeLeftPressure = pipe.GetCellPressure(0); - double pipeRightPressure = pipe.GetCellPressure(pipe.CellCount - 1); - double mdotOrifice = orificeLink.LastMassFlowRate; - double mdotOpen = openEndLink.LastMassFlowRate; + // Apply a Gaussian pulse to the first few cells + double ambientPressure = 101325.0; + double pulseAmplitude = 20 * ambientPressure; // 0.5 atm overpressure + double pulseWidth = 0.05; // m (spread over a few cells) + int n = pipe.CellCount; + double dx = 2.0 / n; - Console.WriteLine($"Step {stepCount}:"); - Console.WriteLine($" Vol Pressure = {volPressure:F1} Pa, Mass = {volMass:E4} kg"); - Console.WriteLine($" Pipe left P = {pipeLeftPressure:F1} Pa, right P = {pipeRightPressure:F1} Pa"); - Console.WriteLine($" Orifice mdot = {mdotOrifice:E4} kg/s, Open‑end mdot = {mdotOpen:E4} kg/s"); - Console.WriteLine(); + // Only modify cells within 2*pulseWidth from the left end + int maxCell = Math.Min(5, n - 1); // at most the first 5 cells + for (int i = 0; i <= maxCell; i++) + { + double x = (i + 0.5) * dx; + double P = pulseAmplitude * Math.Exp(-x * x / (pulseWidth * pulseWidth)); + double currentP = pipe.GetCellPressure(i); + double newP = P; + // Update pressure, keeping density and velocity unchanged + // We recompute total energy accordingly + double rho = pipe.GetCellDensity(i); + double u = pipe.GetCellVelocity(i); + double e = newP / ((1.4 - 1.0) * rho); + double E = rho * e + 0.5 * rho * u * u; + pipe.SetCellState(i, rho, u, newP); + } + Console.WriteLine($"Pulse injected at t = {simTime:F3} s"); + nextPulseTime += pulseInterval; } - // Audio sample from the open‑end mass flow - return (float)openEndLink.LastMassFlowRate; + // Audio from open‑end mass flow + float sample = soundProcessor.Process(openEndLink); + + // Log every 200 steps + if (stepCount % 1000 == 0) + { + int leftIdx = 0; + int midIdx = pipe.CellCount / 2; + int rightIdx = pipe.CellCount - 1; + double pL = pipe.GetCellPressure(leftIdx); + double pM = pipe.GetCellPressure(midIdx); + double pR = pipe.GetCellPressure(rightIdx); + + Console.WriteLine($"Step {stepCount}: P_left={pL:F1} Pa, P_mid={pM:F1} Pa, P_right={pR:F1} Pa"); + } + + if (double.IsNaN(pipe.GetCellPressure(0))) + { + Console.WriteLine("NaN detected – stopping simulation."); + return 0f; + } + + return reverb.Process(sample); } public override void Draw(RenderWindow target) { float winWidth = target.GetView().Size.X; float winHeight = target.GetView().Size.Y; - float pipeCenterY = winHeight / 2f; float margin = 60f; float pipeStartX = margin; float pipeEndX = winWidth - margin; - - // Use the shared pipe drawing from the base class DrawPipe(target, pipe, pipeCenterY, pipeStartX, pipeEndX); } }