diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index dea0d24..4ef0bcc 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -21,15 +21,16 @@ namespace FluidSim.Components private double _dx, _dt, _gamma, _area, _diameter; private double[] _rho, _rhou, _E; - private double _rhoLeft, _pLeft; - private double _rhoRight, _pRight; - private bool _leftBCSet, _rightBCSet; + // Volume‑coupling ghost states for boundaries A and B + private double _rhoA, _pA; + private double _rhoB, _pB; + private bool _aBCSet, _bBCSet; - private BoundaryType _leftBCType = BoundaryType.VolumeCoupling; - private BoundaryType _rightBCType = BoundaryType.VolumeCoupling; + private BoundaryType _aBCType = BoundaryType.VolumeCoupling; + private BoundaryType _bBCType = BoundaryType.VolumeCoupling; - private double _leftAmbientPressure = 101325.0; - private double _rightAmbientPressure = 101325.0; + private double _aAmbientPressure = 101325.0; + private double _bAmbientPressure = 101325.0; private const double CflTarget = 0.8; private const double ReferenceSoundSpeed = 340.0; @@ -39,8 +40,8 @@ namespace FluidSim.Components public double GetCellPressure(int i) => Pressure(i); public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12); - public BoundaryType LeftBCType => _leftBCType; - public BoundaryType RightBCType => _rightBCType; + public BoundaryType ABCType => _aBCType; + public BoundaryType BBCType => _bBCType; public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0) { @@ -76,10 +77,10 @@ namespace FluidSim.Components PortB = new Port(); } - public void SetLeftBoundaryType(BoundaryType type) => _leftBCType = type; - public void SetRightBoundaryType(BoundaryType type) => _rightBCType = type; - public void SetLeftAmbientPressure(double p) => _leftAmbientPressure = p; - public void SetRightAmbientPressure(double p) => _rightAmbientPressure = p; + 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 SetUniformState(double rho, double u, double p) { @@ -102,21 +103,21 @@ namespace FluidSim.Components _E[i] = rho * e + 0.5 * rho * u * u; } - public void SetLeftVolumeState(double rhoVol, double pVol) + public void SetAVolumeState(double rhoVol, double pVol) { - _rhoLeft = rhoVol; - _pLeft = pVol; - _leftBCSet = true; + _rhoA = rhoVol; + _pA = pVol; + _aBCSet = true; } - public void SetRightVolumeState(double rhoVol, double pVol) + public void SetBVolumeState(double rhoVol, double pVol) { - _rhoRight = rhoVol; - _pRight = pVol; - _rightBCSet = true; + _rhoB = rhoVol; + _pB = pVol; + _bBCSet = true; } - public void ClearBC() => _leftBCSet = _rightBCSet = false; + public void ClearBC() => _aBCSet = _bBCSet = false; public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8) { @@ -140,105 +141,90 @@ namespace FluidSim.Components double[] Fp = new double[n + 1]; double[] Fe = new double[n + 1]; - // Left boundary (face 0) - switch (_leftBCType) + // ---------- Boundary A (face 0, left) ---------- + double rhoIntA = _rho[0]; + double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12); + double pIntA = Pressure(0); + + switch (_aBCType) { case BoundaryType.VolumeCoupling: - if (_leftBCSet) + if (_aBCSet) { - HLLCFlux(_rhoLeft, 0.0, _pLeft, - _rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), Pressure(0), + HLLCFlux(_rhoA, 0.0, _pA, + rhoIntA, uIntA, pIntA, out Fm[0], out Fp[0], out Fe[0]); } else { - Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0; + Fm[0] = 0; Fp[0] = pIntA; Fe[0] = 0; } break; case BoundaryType.OpenEnd: - { - double rhoR = _rho[0]; - double uR = _rhou[0] / Math.Max(rhoR, 1e-12); - double pR = Pressure(0); - HLLCFlux(rhoR, uR, _leftAmbientPressure, - rhoR, uR, pR, + OpenEndFluxA(rhoIntA, uIntA, pIntA, _aAmbientPressure, out Fm[0], out Fp[0], out Fe[0]); - } break; case BoundaryType.ClosedEnd: - { - double rhoR = _rho[0]; - double uR = _rhou[0] / Math.Max(rhoR, 1e-12); - double pR = Pressure(0); - HLLCFlux(rhoR, -uR, pR, - rhoR, uR, pR, - out Fm[0], out Fp[0], out Fe[0]); - } + ClosedEndFlux(rhoIntA, uIntA, pIntA, isRightBoundary: false, + out Fm[0], out Fp[0], out Fe[0]); break; } - // Internal faces + // ---------- Internal faces ---------- for (int i = 0; i < n - 1; i++) { - double uL = _rhou[i] / Math.Max(_rho[i], 1e-12); - double uR = _rhou[i + 1] / Math.Max(_rho[i + 1], 1e-12); - HLLCFlux(_rho[i], uL, Pressure(i), - _rho[i + 1], uR, Pressure(i + 1), + 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]); } - // Right boundary (face n) - switch (_rightBCType) + // ---------- 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) { case BoundaryType.VolumeCoupling: - if (_rightBCSet) + if (_bBCSet) { - double rhoL = _rho[n - 1]; - double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12); - double pL = Pressure(n - 1); - HLLCFlux(rhoL, uL, pL, - _rhoRight, 0.0, _pRight, + HLLCFlux(rhoIntB, uIntB, pIntB, + _rhoB, 0.0, _pB, out Fm[n], out Fp[n], out Fe[n]); } else { - Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0; + Fm[n] = 0; Fp[n] = pIntB; Fe[n] = 0; } break; case BoundaryType.OpenEnd: - { - double rhoL = _rho[n - 1]; - double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12); - double pL = Pressure(n - 1); - HLLCFlux(rhoL, uL, pL, - rhoL, uL, _rightAmbientPressure, + OpenEndFluxB(rhoIntB, uIntB, pIntB, _bAmbientPressure, out Fm[n], out Fp[n], out Fe[n]); - } break; case BoundaryType.ClosedEnd: - { - double rhoL = _rho[n - 1]; - double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12); - double pL = Pressure(n - 1); - HLLCFlux(rhoL, uL, pL, - rhoL, -uL, pL, - out Fm[n], out Fp[n], out Fe[n]); - } + ClosedEndFlux(rhoIntB, uIntB, pIntB, isRightBoundary: true, + out Fm[n], out Fp[n], out Fe[n]); break; } // ---- Cell update with linear laminar damping ---- double radius = _diameter / 2.0; - double mu_air = 1.8e-5; // dynamic viscosity of air (Pa·s) + double mu_air = 1.8e-5; double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius); for (int i = 0; i < n; i++) { - // Flux divergence double dM = (Fm[i + 1] - Fm[i]) / _dx; double dP = (Fp[i + 1] - Fp[i]) / _dx; double dE = (Fe[i + 1] - Fe[i]) / _dx; @@ -247,14 +233,11 @@ namespace FluidSim.Components _rhou[i] -= dtSub * dP; _E[i] -= dtSub * dE; - // Laminar viscous damping on momentum (implicit exponential decay) double rho = Math.Max(_rho[i], 1e-12); - double dampingRate = laminarCoeff / rho; // 1/s + double dampingRate = laminarCoeff / rho; double dampingFactor = Math.Exp(-dampingRate * dtSub); _rhou[i] *= dampingFactor; - // Note: total energy _E[i] is unchanged – kinetic energy loss becomes internal heat - // Physical bounds if (_rho[i] < 1e-12) _rho[i] = 1e-12; double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i]; double pMin = 100.0; @@ -262,28 +245,29 @@ namespace FluidSim.Components if (_E[i] < eMin) _E[i] = eMin; } - // Port quantities (only meaningful for volume‑coupled ends) - double mdotA_sub = _leftBCType == BoundaryType.VolumeCoupling && _leftBCSet ? Fm[0] * _area : 0.0; - double mdotB_sub = _rightBCType == BoundaryType.VolumeCoupling && _rightBCSet ? -Fm[n] * _area : 0.0; + // ---------- 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 = Pressure(0); - PortB.Pressure = Pressure(_n - 1); + PortA.Pressure = pIntA; + PortB.Pressure = pIntB; PortA.Density = _rho[0]; - PortB.Density = _rho[_n - 1]; + PortB.Density = _rho[n - 1]; - if (_leftBCType == BoundaryType.VolumeCoupling && _leftBCSet) + // Corrected enthalpy for both directions + if (_aBCType == BoundaryType.VolumeCoupling && _aBCSet) { PortA.SpecificEnthalpy = mdotA_sub < 0 ? GetCellTotalSpecificEnthalpy(0) - : 0.0; + : (_gamma / (_gamma - 1.0)) * _pA / Math.Max(_rhoA, 1e-12); } - if (_rightBCType == BoundaryType.VolumeCoupling && _rightBCSet) + if (_bBCType == BoundaryType.VolumeCoupling && _bBCSet) { PortB.SpecificEnthalpy = mdotB_sub < 0 ? GetCellTotalSpecificEnthalpy(_n - 1) - : 0.0; + : (_gamma / (_gamma - 1.0)) * _pB / Math.Max(_rhoB, 1e-12); } } @@ -299,6 +283,98 @@ namespace FluidSim.Components 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, out double fm, out double fp, out double fe) diff --git a/Core/Solver.cs b/Core/Solver.cs index 5b5b2c6..a85f51d 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -1,4 +1,5 @@ -using System.Collections.Generic; +using System; +using System.Collections.Generic; using FluidSim.Components; using FluidSim.Interfaces; @@ -10,93 +11,90 @@ namespace FluidSim.Core private readonly List _pipes = new(); private readonly List _connections = new(); - private double _dt; // global time step + private double _dt; public void AddVolume(Volume0D v) => _volumes.Add(v); public void AddPipe(Pipe1D p) => _pipes.Add(p); public void AddConnection(Connection c) => _connections.Add(c); - - /// Set the global time step (called from Simulation). public void SetTimeStep(double dt) => _dt = dt; /// - /// Convenient method to set the boundary type of a pipe end. + /// Set boundary type for a pipe end. isA = true for port A (left), false for port B (right). /// - public void SetPipeBoundary(Pipe1D pipe, bool isLeft, BoundaryType type, double ambientPressure = 101325.0) + public void SetPipeBoundary(Pipe1D pipe, bool isA, BoundaryType type, double ambientPressure = 101325.0) { - if (isLeft) + if (isA) { - pipe.SetLeftBoundaryType(type); + pipe.SetABoundaryType(type); if (type == BoundaryType.OpenEnd) - pipe.SetLeftAmbientPressure(ambientPressure); + pipe.SetAAmbientPressure(ambientPressure); } else { - pipe.SetRightBoundaryType(type); + pipe.SetBBoundaryType(type); if (type == BoundaryType.OpenEnd) - pipe.SetRightAmbientPressure(ambientPressure); + pipe.SetBAmbientPressure(ambientPressure); } } public float Step() { - // 1. Volumes publish state to ports (only needed if any volume exists) + // 1. Volumes publish state foreach (var v in _volumes) v.PushStateToPort(); - // 2. Set initial pipe boundary conditions ONLY for volume‑coupled ends + // 2. Set volume BCs for volume‑coupled ends foreach (var conn in _connections) { if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) { var pipe = GetPipe(conn.PortA); - bool isLeft = pipe.PortA == conn.PortA; - BoundaryType bc = isLeft ? pipe.LeftBCType : pipe.RightBCType; - if (bc == BoundaryType.VolumeCoupling) + 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 isLeft = pipe.PortB == conn.PortB; - BoundaryType bc = isLeft ? pipe.LeftBCType : pipe.RightBCType; - if (bc == BoundaryType.VolumeCoupling) + bool isA = pipe.PortB == conn.PortB; + if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) || + (!isA && pipe.BBCType == BoundaryType.VolumeCoupling)) SetVolumeBC(conn.PortB, conn.PortA); } } - // 3. Determine number of sub‑steps + // 3. Sub‑steps int nSub = 1; foreach (var p in _pipes) nSub = Math.Max(nSub, p.GetRequiredSubSteps(_dt)); double dtSub = _dt / nSub; - // 4. Sub‑step loop for (int sub = 0; sub < nSub; sub++) { foreach (var p in _pipes) p.SimulateSingleStep(dtSub); - // Transfer flows only for volume‑coupled connections foreach (var conn in _connections) { if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) { var pipe = GetPipe(conn.PortA); - bool isLeft = pipe.PortA == conn.PortA; - if (pipe.LeftBCType == BoundaryType.VolumeCoupling || pipe.RightBCType == BoundaryType.VolumeCoupling) + 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 isLeft = pipe.PortB == conn.PortB; - if (pipe.LeftBCType == BoundaryType.VolumeCoupling || pipe.RightBCType == BoundaryType.VolumeCoupling) + bool isA = pipe.PortB == conn.PortB; + if ((isA && pipe.ABCType == BoundaryType.VolumeCoupling) || + (!isA && pipe.BBCType == BoundaryType.VolumeCoupling)) TransferAndIntegrate(conn.PortB, conn.PortA, dtSub); } } - // Update BCs for volume‑coupled ends between sub‑steps if (sub < nSub - 1) { foreach (var v in _volumes) @@ -107,24 +105,24 @@ namespace FluidSim.Core if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) { var pipe = GetPipe(conn.PortA); - bool isLeft = pipe.PortA == conn.PortA; - if ((isLeft && pipe.LeftBCType == BoundaryType.VolumeCoupling) || - (!isLeft && pipe.RightBCType == BoundaryType.VolumeCoupling)) + 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 isLeft = pipe.PortB == conn.PortB; - if ((isLeft && pipe.LeftBCType == BoundaryType.VolumeCoupling) || - (!isLeft && pipe.RightBCType == BoundaryType.VolumeCoupling)) + 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 from SoundConnections (if any) + // 5. Audio samples (none for now, but placeholder) var audioSamples = new List(); foreach (var conn in _connections) { @@ -132,7 +130,7 @@ namespace FluidSim.Core audioSamples.Add(sc.GetAudioSample()); } - // 6. Clear volume BC flags + // 6. Clear BC flags foreach (var p in _pipes) p.ClearBC(); @@ -148,11 +146,11 @@ namespace FluidSim.Core { var pipe = GetPipe(pipePort); if (pipe == null) return; - bool isLeft = pipe.PortA == pipePort; - if (isLeft) - pipe.SetLeftVolumeState(volPort.Density, volPort.Pressure); + bool isA = pipe.PortA == pipePort; + if (isA) + pipe.SetAVolumeState(volPort.Density, volPort.Pressure); else - pipe.SetRightVolumeState(volPort.Density, volPort.Pressure); + pipe.SetBVolumeState(volPort.Density, volPort.Pressure); } private void TransferAndIntegrate(Port pipePort, Port volPort, double dtSub) @@ -164,7 +162,7 @@ namespace FluidSim.Core { volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy; } - // else: volume’s own enthalpy (set by PushStateToPort) is used + // else volume’s own enthalpy (from PushStateToPort) is used GetVolume(volPort)?.Integrate(dtSub); } diff --git a/Program.cs b/Program.cs index 3e7f273..c7b75cb 100644 --- a/Program.cs +++ b/Program.cs @@ -13,7 +13,8 @@ public class Program private static Scenario scenario; // Speed control - private static double desiredSpeed = 1.0; + //private static double desiredSpeed = 1.0; + private static double desiredSpeed = 0.0001; private static double currentSpeed = desiredSpeed; private const double MinSpeed = 0.0001; private const double MaxSpeed = 1.0; @@ -38,7 +39,10 @@ public class Program soundEngine.Volume = 70; soundEngine.Start(); - scenario = new PipeResonatorScenario(); + //scenario = new PipeResonatorScenario(); + //scenario = new HelmholtzResonatorScenario(); + scenario = new SodShockTubeScenario(); + scenario.Initialize(SampleRate); var stopwatch = Stopwatch.StartNew(); diff --git a/Scenarios/HelmholtzResonatorScenario.cs b/Scenarios/HelmholtzResonatorScenario.cs new file mode 100644 index 0000000..293b86b --- /dev/null +++ b/Scenarios/HelmholtzResonatorScenario.cs @@ -0,0 +1,133 @@ +using System; +using FluidSim.Components; +using FluidSim.Interfaces; +using FluidSim.Utils; +using SFML.Graphics; +using SFML.System; + +namespace FluidSim.Core +{ + public class HelmholtzResonatorScenario : Scenario + { + private Solver solver; + private Volume0D cavity; + private Pipe1D neck; + private Connection coupling; + private int stepCount; + private double time; + private double dt; + private double ambientPressure = 1.0 * Units.atm; + + public override void Initialize(int sampleRate) + { + dt = 1.0 / sampleRate; + + // 1‑litre cavity, 10% over‑pressure + double cavityVolume = 1e-3; + double initialCavityPressure = 1.1 * ambientPressure; + cavity = new Volume0D(cavityVolume, initialCavityPressure, 300.0, sampleRate) + { + Gamma = 1.4, + GasConstant = 287.0 + }; + + // Neck: length 10 cm, radius 1 cm + double neckLength = 0.1; + double neckRadius = 0.01; + double neckArea = Math.PI * neckRadius * neckRadius; + neck = new Pipe1D(neckLength, neckArea, sampleRate, forcedCellCount: 40); + neck.SetUniformState(1.225, 0.0, ambientPressure); + + coupling = new Connection(neck.PortA, cavity.Port) + { + Area = neckArea, + DischargeCoefficient = 0.62, + Gamma = 1.4 + }; + + solver = new Solver(); + solver.SetTimeStep(dt); + solver.AddVolume(cavity); + solver.AddPipe(neck); + solver.AddConnection(coupling); + + // Port A (left) = volume coupling, Port B (right) = open end + solver.SetPipeBoundary(neck, isA: true, BoundaryType.VolumeCoupling); + solver.SetPipeBoundary(neck, isA: false, BoundaryType.OpenEnd, ambientPressure); + } + + public override float Process() + { + float sample = solver.Step(); + time += dt; + stepCount++; + + double pOpen = neck.GetCellPressure(neck.GetCellCount() - 1); + float audio = (float)((pOpen - ambientPressure) / ambientPressure); + + if (stepCount % 20 == 0) + { + double pCav = cavity.Pressure; + double mdotA = neck.PortA.MassFlowRate; // positive = into pipe (leaving cavity) + 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}"); + } + + return audio; + } + + public override void Draw(RenderWindow target) + { + float winW = target.GetView().Size.X; + float winH = target.GetView().Size.Y; + float centerY = winH / 2f; + + // Cavity rectangle + float cavityWidth = 120f; + float cavityHeight = 180f; + var cavityRect = new RectangleShape(new Vector2f(cavityWidth, cavityHeight)); + cavityRect.Position = new Vector2f(40f, centerY - cavityHeight / 2f); + cavityRect.FillColor = PressureColor(cavity.Pressure); + target.Draw(cavityRect); + + // Neck drawn as tapered pipe + int n = neck.GetCellCount(); + float neckStartX = 40f + cavityWidth + 10f; + float neckEndX = winW - 60f; + float neckLenPx = neckEndX - neckStartX; + float dx = neckLenPx / (n - 1); + float baseRadius = 20f; + + Vertex[] vertices = new Vertex[n * 2]; + for (int i = 0; i < n; i++) + { + float x = neckStartX + i * dx; + double p = neck.GetCellPressure(i); + float r = baseRadius * (float)(0.5 + 0.5 * Math.Tanh((p - ambientPressure) / (ambientPressure * 0.2))); + if (r < 4f) r = 4f; + Color col = PressureColor(p); + vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col); + vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col); + } + target.Draw(vertices, PrimitiveType.TriangleStrip); + + // Open end indicator + var arrow = new CircleShape(8f); + arrow.Position = new Vector2f(neckEndX - 4f, centerY - 4f); + arrow.FillColor = Color.White; + target.Draw(arrow); + } + + private Color PressureColor(double pressure) + { + double range = ambientPressure * 0.1; + double t = Math.Clamp((pressure - ambientPressure) / range, -1.0, 1.0); + byte r = (byte)(t > 0 ? 255 * t : 0); + byte b = (byte)(t < 0 ? -255 * t : 0); + byte g = (byte)(255 * (1 - Math.Abs(t))); + return new Color(r, g, b); + } + } +} \ No newline at end of file diff --git a/Scenarios/PipeResonatorScenario.cs b/Scenarios/PipeResonatorScenario.cs index 5496c14..e534910 100644 --- a/Scenarios/PipeResonatorScenario.cs +++ b/Scenarios/PipeResonatorScenario.cs @@ -21,7 +21,7 @@ namespace FluidSim.Core { dt = 1.0 / sampleRate; - double length = 0.5; + double length = 2; double radius = 50 * Units.mm; double area = Units.AreaFromDiameter(radius); @@ -31,8 +31,9 @@ namespace FluidSim.Core solver = new Solver(); solver.SetTimeStep(dt); solver.AddPipe(pipe); - solver.SetPipeBoundary(pipe, isLeft: true, BoundaryType.OpenEnd, ambientPressure); - solver.SetPipeBoundary(pipe, isLeft: false, BoundaryType.ClosedEnd); + // Open end at port A (left), closed end at port B (right) + solver.SetPipeBoundary(pipe, isA: true, BoundaryType.OpenEnd, ambientPressure); + solver.SetPipeBoundary(pipe, isA: false, BoundaryType.ClosedEnd); // Initial pressure pulse int pulseCells = 5; diff --git a/Scenarios/SodShockTubeScenario.cs b/Scenarios/SodShockTubeScenario.cs new file mode 100644 index 0000000..ff06758 --- /dev/null +++ b/Scenarios/SodShockTubeScenario.cs @@ -0,0 +1,158 @@ +using System; +using FluidSim.Components; +using FluidSim.Utils; +using SFML.Graphics; +using SFML.System; + +namespace FluidSim.Core +{ + public class SodShockTubeScenario : Scenario + { + private Solver solver; + private Pipe1D pipe; + private int stepCount; + private double time; + private double dt; + private double ambientPressure = 1.0 * Units.atm; + private const double GasConstant = 287.0; + + public override void Initialize(int sampleRate) + { + dt = 1.0 / sampleRate; + double length = 1.0; + double area = 1.0; + int nCells = 200; + + pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: nCells); + pipe.SetUniformState(0.125, 0.0, 0.1 * ambientPressure); // right state + + // Left half high pressure + for (int i = 0; i < nCells / 2; i++) + pipe.SetCellState(i, 1.0, 0.0, ambientPressure); + + solver = new Solver(); + solver.SetTimeStep(dt); + solver.AddPipe(pipe); + solver.SetPipeBoundary(pipe, isA: true, BoundaryType.ClosedEnd); + solver.SetPipeBoundary(pipe, isA: false, BoundaryType.ClosedEnd); + } + + public override float Process() + { + float sample = solver.Step(); + time += dt; + stepCount++; + + double pMid = pipe.GetPressureAtFraction(0.5); + float audio = (float)((pMid - ambientPressure) / ambientPressure); + + bool log = true; + + if (log) + { + int n = pipe.GetCellCount(); + Console.WriteLine($"step {stepCount}:"); + Console.WriteLine("i rho (kg/m³) p (Pa) T (K) u (m/s)"); + for (int i = 0; i < n; i++) + { + if (i % 10 == 0) + { + double rho = pipe.GetCellDensity(i); + double p = pipe.GetCellPressure(i); + double u = pipe.GetCellVelocity(i); + double T = p / (rho * GasConstant); // GasConstant = 287.0 + Console.WriteLine($"{i,-4} {rho,10:F4} {p,10:F1} {T,8:F2} {u,10:F4}"); + } + } + Console.WriteLine(); + } + + return audio; + } + + public override void Draw(RenderWindow target) + { + float winW = target.GetView().Size.X; + float winH = target.GetView().Size.Y; + float centerY = winH / 2f; + float margin = 40f; + float pipeStartX = margin; + float pipeEndX = winW - margin; + float pipeLenPx = pipeEndX - pipeStartX; + int n = pipe.GetCellCount(); + float dx = pipeLenPx / (n - 1); + float baseRadius = 60f; + + Vertex[] vertices = new Vertex[n * 2]; + for (int i = 0; i < n; i++) + { + float x = pipeStartX + i * dx; + + double p = pipe.GetCellPressure(i); + double rho = pipe.GetCellDensity(i); + double T = p / (rho * GasConstant); // temperature in Kelvin + + // Radius from pressure (exaggerated deviation) + float r = baseRadius * (float)(p / ambientPressure * 2); + if (r < 4f) r = 4f; + + // Colour from temperature + Color col = TemperatureColor(T); + + vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col); + vertices[i * 2 + 1] = new Vertex(new Vector2f(x, centerY + r), col); + } + target.Draw(vertices, PrimitiveType.TriangleStrip); + + // Diaphragm marker (faint white line at initial interface) + float diaphragmX = pipeStartX + (n / 2) * dx; + var line = new RectangleShape(new Vector2f(2f, winH * 0.5f)); + line.Position = new Vector2f(diaphragmX - 1f, centerY - winH * 0.25f); + line.FillColor = new Color(255, 255, 255, 80); + target.Draw(line); + } + + /// + /// Custom temperature‑to‑hue mapping that matches the given Sod‑tube hue values: + /// 250 K → 176, 300 K → 122, 350 K → 120?, 450 K → 71. + /// Interpolates piecewise linearly, clamping outside [250,450]. + /// + private Color TemperatureColor(double T) + { + // 1. Map temperature to hue (0–360) + double[] Tknots = { 250, 282, 353, 450 }; + double[] Hknots = { 176, 179, 122, 71 }; + double hue; + if (T <= Tknots[0]) hue = Hknots[0]; + else if (T >= Tknots[^1]) hue = Hknots[^1]; + else + { + int i = 0; + while (i < Tknots.Length - 1 && T > Tknots[i + 1]) i++; + double frac = (T - Tknots[i]) / (Tknots[i + 1] - Tknots[i]); + hue = Hknots[i] + frac * (Hknots[i + 1] - Hknots[i]); + } + + // 2. Convert hue to RGB (S = 1, V = 1) + double h = hue / 60.0; + int sector = (int)Math.Floor(h); + double f = h - sector; + byte p = 0; + byte q = (byte)(255 * (1 - f)); + byte tByte = (byte)(255 * f); + byte v = 255; + + byte r, g, b; + switch (sector % 6) + { + case 0: r = v; g = tByte; b = p; break; + case 1: r = q; g = v; b = p; break; + case 2: r = p; g = v; b = tByte; break; + case 3: r = p; g = q; b = v; break; + case 4: r = tByte; g = p; b = v; break; + default: r = v; g = p; b = q; break; + } + return new Color(r, g, b); + } + } +} \ No newline at end of file