From 650b993c2b618646c3bf43f43be7ea755e4a9420 Mon Sep 17 00:00:00 2001 From: max Date: Sat, 2 May 2026 17:40:43 +0200 Subject: [PATCH] Known bad point (unphysical energy loss) --- Components/Pipe1D.cs | 22 +++++++++---- Core/OrificeBoundary.cs | 13 +++++--- Core/Simulation.cs | 71 +++++++++++++++-------------------------- Core/Solver.cs | 54 +++++++++++++++++++++---------- Interfaces/Junction0.cs | 10 ------ Interfaces/Junction1.cs | 10 ------ 6 files changed, 89 insertions(+), 91 deletions(-) delete mode 100644 Interfaces/Junction0.cs delete mode 100644 Interfaces/Junction1.cs diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index 027af15..53d45cf 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -67,6 +67,8 @@ namespace FluidSim.Components _fxR_mass = m; _fxR_mom = p; _fxR_ener = e; _rightSet = true; } + + public void Simulate() { int n = _n; @@ -106,16 +108,15 @@ namespace FluidSim.Components if (_E[i] < kinetic) _E[i] = kinetic; } - // Friction + // Friction (energy‑conserving) if (FrictionFactor > 0) { double D = 2.0 * Math.Sqrt(_area / Math.PI); - for (int i = 0; i < n; i++) + for (int i = 0; i < _n; i++) { double u = _rhou[i] / Math.Max(_rho[i], 1e-12); double f = FrictionFactor / (2.0 * D) * _rho[i] * u * Math.Abs(u); - _rhou[i] -= _dt * f; - if (_E[i] > _dt * f * u) _E[i] -= _dt * f * u; + //_rhou[i] -= _dt * f; FRICTIN DISABLED!!! } } @@ -124,13 +125,22 @@ namespace FluidSim.Components PortB.MassFlowRate = _rightSet ? -_fxR_mass * _area : 0.0; // Enthalpy for upwinding - PortA.SpecificEnthalpy = _gamma / (_gamma - 1.0) * Pressure(0) / Math.Max(_rho[0], 1e-12); - PortB.SpecificEnthalpy = _gamma / (_gamma - 1.0) * Pressure(_n - 1) / Math.Max(_rho[_n - 1], 1e-12); + PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0); + PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1); // Reset for next step _leftSet = _rightSet = false; } + 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; + } + double Pressure(int i) => (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12)); diff --git a/Core/OrificeBoundary.cs b/Core/OrificeBoundary.cs index c4c22c2..733eba2 100644 --- a/Core/OrificeBoundary.cs +++ b/Core/OrificeBoundary.cs @@ -55,14 +55,19 @@ namespace FluidSim.Core bool isLeftBoundary, out double massFlux, out double momFlux, out double energyFlux) { - // mass flow from pipe to volume (positive = pipe → volume) - double mdot = MassFlow(pPipe, rhoPipe, pVol, rhoVol, conn); + // ----- Compute STAGNATION pressures ----- + double pStagPipe = pPipe + 0.5 * rhoPipe * uPipe * uPipe; + double pStagVol = pVol + 0.5 * rhoVol * uVol * uVol; // uVol is always 0 for your volumes + + // Mass flow driven by stagnation pressure difference (positive = pipe→volume) + double mdot = MassFlow(pStagPipe, rhoPipe, pStagVol, rhoVol, conn); // Limit mass flow to the amount that can leave/enter the pipe cell double maxMdot = rhoPipe * pipeArea * 343.0; if (Math.Abs(mdot) > maxMdot) mdot = Math.Sign(mdot) * maxMdot; - bool flowLeavesPipe = mdot > 0; + bool flowLeavesPipe = mdot > 0; // pipe → volume + double uFace, pFace, rhoFace; double massFluxPerArea; @@ -83,7 +88,7 @@ namespace FluidSim.Core { uFace = uVol; pFace = pVol; rhoFace = rhoVol; } } - // Total enthalpy of the injected fluid (corrected: mass flux × total enthalpy) + // Total enthalpy of the injected fluid double specificEnthalpy = (1.4 / (1.4 - 1.0)) * pFace / Math.Max(rhoFace, 1e-12); double totalEnthalpy = specificEnthalpy + 0.5 * uFace * uFace; diff --git a/Core/Simulation.cs b/Core/Simulation.cs index 45a9a3c..2b310c5 100644 --- a/Core/Simulation.cs +++ b/Core/Simulation.cs @@ -6,58 +6,46 @@ namespace FluidSim.Core { public static class Simulation { + private static Solver solver; + private static Volume0D volA, volB; private static Pipe1D pipe; - private static Connection leftConn, rightConn; // dummy connections for orifice params + private static Connection connA, connB; + private static int stepCount; private static double time; private static double dt; - private static int stepCount; public static void Initialize(int sampleRate) { dt = 1.0 / sampleRate; + double V = 5.0 * Units.L; + volA = new Volume0D(V, 1.1 * Units.atm, Units.Celsius(20), sampleRate); + volB = new Volume0D(V, 1.0 * Units.atm, Units.Celsius(20), sampleRate); + double length = 150 * Units.mm; double diameter = 25 * Units.mm; double area = Units.AreaFromDiameter(25, Units.mm); int nCells = 10; pipe = new Pipe1D(length, area, nCells, sampleRate); - pipe.SetUniformState(1.2, 0.0, 1.0 * Units.atm); // start at 1 atm + pipe.SetUniformState(volA.Density, 0.0, volA.Pressure); pipe.FrictionFactor = 0.02; - // Dummy connections – only used for orifice parameters - leftConn = new Connection(null, null) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 }; - rightConn = new Connection(null, null) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 }; + // Connections with orifice area equal to pipe area (flange joint) + connA = new Connection(volA.Port, pipe.PortA) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 }; + connB = new Connection(pipe.PortB, volB.Port) { Area = area, DischargeCoefficient = 1.0, Gamma = 1.4 }; + + solver = new Solver(); + solver.AddVolume(volA); + solver.AddVolume(volB); + solver.AddPipe(pipe); + solver.AddConnection(connA); + solver.AddConnection(connB); } public static float Process() { - // Fixed boundary reservoirs - double pLeft = 1.1 * Units.atm; - double rhoLeft = 1.2; - double uLeft = 0.0; - - double pRight = 1.0 * Units.atm; - double rhoRight = 1.2; - double uRight = 0.0; - - // Compute boundary fluxes via orifice model - OrificeBoundary.PipeVolumeFlux( - pipe.GetLeftPressure(), pipe.GetLeftDensity(), 0.0, - pLeft, rhoLeft, uLeft, - leftConn, pipe.Area, true, - out double leftMassFlux, out double leftMomFlux, out double leftEnergyFlux); - - OrificeBoundary.PipeVolumeFlux( - pipe.GetRightPressure(), pipe.GetRightDensity(), 0.0, - pRight, rhoRight, uRight, - rightConn, pipe.Area, false, - out double rightMassFlux, out double rightMomFlux, out double rightEnergyFlux); - - pipe.SetLeftBoundaryFlux(leftMassFlux, leftMomFlux, leftEnergyFlux); - pipe.SetRightBoundaryFlux(rightMassFlux, rightMomFlux, rightEnergyFlux); - pipe.Simulate(); - + solver.Step(); time += dt; stepCount++; Log(); @@ -66,20 +54,13 @@ namespace FluidSim.Core public static void Log() { - if (stepCount <= 20 || stepCount % 50 == 0) + if (stepCount <= 50 || stepCount % 200 == 0) { - Console.WriteLine($"Step {stepCount:D4} t = {time * 1e3:F3} ms"); - for (int i = 0; i < pipe.GetCellCount(); i++) - { - double rho = pipe.GetCellDensity(i); - double p = pipe.GetCellPressure(i); - double u = pipe.GetCellVelocity(i); - Console.WriteLine($" Cell {i}: ρ={rho:F4} kg/m³ p={p / 1e5:F6} bar u={u:F3} m/s"); - } - double leftFlow = pipe.PortA.MassFlowRate; - double rightFlow = pipe.PortB.MassFlowRate; - Console.WriteLine($" Left flow = {leftFlow * 1e3:F4} g/s Right flow = {rightFlow * 1e3:F4} g/s"); - Console.WriteLine(); + Console.WriteLine( + $"t = {time * 1e3:F3} ms Step {stepCount:D4}: " + + $"PA = {volA.Pressure / 1e5:F6} bar, " + + $"PB = {volB.Pressure / 1e5:F6} bar, " + + $"FlowA = {pipe.PortA.MassFlowRate * 1e3:F2} g/s"); } } } diff --git a/Core/Solver.cs b/Core/Solver.cs index 58c2cf6..f1e5e1f 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -20,7 +20,7 @@ namespace FluidSim.Core foreach (var v in _volumes) v.PushStateToPort(); - // 2. Apply orifice boundaries to pipes + // 2. Compute boundary fluxes (orifice model) foreach (var conn in _connections) { if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) @@ -31,20 +31,20 @@ namespace FluidSim.Core VolumeToVolume(conn); } - // 3. Pipes simulate + // 3. Pipe simulation step foreach (var p in _pipes) p.Simulate(); - // 4. Transfer pipe flows to connected volumes + // 4. Transfer pipe‑port data to volumes foreach (var conn in _connections) { if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) - Transfer(conn.PortA, conn.PortB); + TransferPipeToVolume(conn.PortA, conn.PortB); else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB)) - Transfer(conn.PortB, conn.PortA); + TransferPipeToVolume(conn.PortB, conn.PortA); } - // 5. Volumes integrate + // 5. Integrate volumes foreach (var v in _volumes) v.Integrate(); } @@ -61,34 +61,56 @@ namespace FluidSim.Core double pP = isLeft ? pipe.GetLeftPressure() : pipe.GetRightPressure(); double rhoP = isLeft ? pipe.GetLeftDensity() : pipe.GetRightDensity(); - double uP = 0.0; - double pV = volPort.Pressure, rhoV = volPort.Density, uV = 0.0; + double uP = isLeft ? pipe.GetCellVelocity(0) + : pipe.GetCellVelocity(pipe.GetCellCount() - 1); + + double pV = volPort.Pressure; + double rhoV = volPort.Density; + double uV = 0.0; // volume has zero organized velocity + + OrificeBoundary.PipeVolumeFlux( + pP, rhoP, uP, + pV, rhoV, uV, + conn, pipe.Area, isLeft, + out double massFlux, out double momFlux, out double energyFlux); - OrificeBoundary.PipeVolumeFlux(pP, rhoP, uP, pV, rhoV, uV, conn, pipe.Area, - isLeft, out double mf, out double pf, out double ef); if (isLeft) - pipe.SetLeftBoundaryFlux(mf, pf, ef); + pipe.SetLeftBoundaryFlux(massFlux, momFlux, energyFlux); else - pipe.SetRightBoundaryFlux(mf, pf, ef); + pipe.SetRightBoundaryFlux(massFlux, momFlux, energyFlux); } void VolumeToVolume(Connection conn) { - double mdot = OrificeBoundary.MassFlow(conn.PortA.Pressure, conn.PortA.Density, - conn.PortB.Pressure, conn.PortB.Density, conn); + double mdot = OrificeBoundary.MassFlow( + conn.PortA.Pressure, conn.PortA.Density, + conn.PortB.Pressure, conn.PortB.Density, conn); + conn.PortA.MassFlowRate = -mdot; conn.PortB.MassFlowRate = mdot; + if (mdot > 0) conn.PortB.SpecificEnthalpy = conn.PortA.SpecificEnthalpy; else if (mdot < 0) conn.PortA.SpecificEnthalpy = conn.PortB.SpecificEnthalpy; } - void Transfer(Port pipePort, Port volPort) + void TransferPipeToVolume(Port pipePort, Port volPort) { double mdot = pipePort.MassFlowRate; + // mdot > 0 → fluid enters pipe from volume + // mdot < 0 → fluid leaves pipe and enters volume + + // Volume mass flow sign is opposite (positive into volume) volPort.MassFlowRate = -mdot; - volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy; + + if (mdot < 0) // pipe → volume + { + // ★ pipePort.SpecificEnthalpy now contains TOTAL enthalpy + volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy; + } + // else: fluid goes volume → pipe → volume owns its own (static) enthalpy, + // which is already correct. } } } \ No newline at end of file diff --git a/Interfaces/Junction0.cs b/Interfaces/Junction0.cs deleted file mode 100644 index 882a7d1..0000000 --- a/Interfaces/Junction0.cs +++ /dev/null @@ -1,10 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Text; - -namespace FluidSim.Interfaces -{ - internal class Junction0 - { - } -} diff --git a/Interfaces/Junction1.cs b/Interfaces/Junction1.cs deleted file mode 100644 index 5ecc76d..0000000 --- a/Interfaces/Junction1.cs +++ /dev/null @@ -1,10 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Text; - -namespace FluidSim.Interfaces -{ - internal class Junction1 - { - } -}