From ab50b808fb1f01efd15f928c32a359ee09e05487 Mon Sep 17 00:00:00 2001 From: max Date: Sat, 2 May 2026 18:05:14 +0200 Subject: [PATCH] Energy conservation fixed --- Components/Pipe1D.cs | 169 ++++++++++++++++++++++++------------------- Core/Solver.cs | 59 +++------------ 2 files changed, 107 insertions(+), 121 deletions(-) diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index 53d45cf..5b6e2aa 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -13,10 +13,9 @@ namespace FluidSim.Components private double _dx, _dt, _gamma = 1.4, _area; private double[] _rho, _rhou, _E; - // Boundary fluxes (set by solver before each step) - private double _fxL_mass, _fxL_mom, _fxL_ener; - private double _fxR_mass, _fxR_mom, _fxR_ener; - private bool _leftSet, _rightSet; + // Volume states at boundaries + private double _rhoLeft, _pLeft, _rhoRight, _pRight; + private bool _leftBCSet, _rightBCSet; public double FrictionFactor { get; set; } = 0.02; @@ -57,79 +56,19 @@ namespace FluidSim.Components public double GetLeftDensity() => _rho[0]; public double GetRightDensity() => _rho[_n - 1]; - public void SetLeftBoundaryFlux(double m, double p, double e) + // ★ New: pass both density and pressure from the volume + public void SetLeftVolumeState(double rhoVol, double pVol) { - _fxL_mass = m; _fxL_mom = p; _fxL_ener = e; _leftSet = true; + _rhoLeft = rhoVol; + _pLeft = pVol; + _leftBCSet = true; } - public void SetRightBoundaryFlux(double m, double p, double e) + public void SetRightVolumeState(double rhoVol, double pVol) { - _fxR_mass = m; _fxR_mom = p; _fxR_ener = e; _rightSet = true; - } - - - - public void Simulate() - { - int n = _n; - double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1]; - - // Left face - if (_leftSet) { Fm[0] = _fxL_mass; Fp[0] = _fxL_mom; Fe[0] = _fxL_ener; } - else { Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0; } // reflective wall - - // Internal faces (HLLC) - 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), - out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]); - } - - // Right face - if (_rightSet) { Fm[n] = _fxR_mass; Fp[n] = _fxR_mom; Fe[n] = _fxR_ener; } - else { Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0; } - - // Update cells - 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] -= _dt * dM; - _rhou[i] -= _dt * dP; - _E[i] -= _dt * dE; - - // Clamp to physical - if (_rho[i] < 1e-12) _rho[i] = 1e-12; - double u = _rhou[i] / _rho[i]; - double kinetic = 0.5 * _rho[i] * u * u; - if (_E[i] < kinetic) _E[i] = kinetic; - } - - // Friction (energy‑conserving) - if (FrictionFactor > 0) - { - double D = 2.0 * Math.Sqrt(_area / Math.PI); - 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; FRICTIN DISABLED!!! - } - } - - // Write port flows for the solver - PortA.MassFlowRate = _leftSet ? _fxL_mass * _area : 0.0; - PortB.MassFlowRate = _rightSet ? -_fxR_mass * _area : 0.0; - - // Enthalpy for upwinding - PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0); - PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1); - - // Reset for next step - _leftSet = _rightSet = false; + _rhoRight = rhoVol; + _pRight = pVol; + _rightBCSet = true; } private double GetCellTotalSpecificEnthalpy(int i) @@ -141,6 +80,90 @@ namespace FluidSim.Components return h + 0.5 * u * u; } + public void Simulate() + { + int n = _n; + double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1]; + + // --- Left boundary (face 0) --- + if (_leftBCSet) + { + // Ghost = actual volume state (ρ_vol, u=0, p_vol) + double rhoL = _rhoLeft; + double uL = 0.0; + double pL = _pLeft; + + double rhoR = _rho[0]; + double uR = _rhou[0] / Math.Max(rhoR, 1e-12); + double pR = Pressure(0); + + HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[0], out Fp[0], out Fe[0]); + } + else + { + Fm[0] = 0; + Fp[0] = Pressure(0); + Fe[0] = 0; + } + + // --- 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), + out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]); + } + + // --- Right boundary (face n) --- + if (_rightBCSet) + { + double rhoL = _rho[n - 1]; + double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12); + double pL = Pressure(n - 1); + + // Ghost = actual volume state (ρ_vol, u=0, p_vol) + double rhoR = _rhoRight; + double uR = 0.0; + double pR = _pRight; + + HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[n], out Fp[n], out Fe[n]); + } + else + { + Fm[n] = 0; + Fp[n] = Pressure(n - 1); + Fe[n] = 0; + } + + // --- Cell update --- + 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] -= _dt * dM; + _rhou[i] -= _dt * dP; + _E[i] -= _dt * dE; + + if (_rho[i] < 1e-12) _rho[i] = 1e-12; + double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i]; + if (_E[i] < kinetic) _E[i] = kinetic; + } + + // --- Friction disabled --- + // if (FrictionFactor > 0) { … } + + // --- Port flows --- + PortA.MassFlowRate = _leftBCSet ? Fm[0] * _area : 0.0; + PortB.MassFlowRate = _rightBCSet ? -Fm[n] * _area : 0.0; + + PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0); + PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1); + + _leftBCSet = _rightBCSet = false; + } + 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/Solver.cs b/Core/Solver.cs index f1e5e1f..0c7807f 100644 --- a/Core/Solver.cs +++ b/Core/Solver.cs @@ -16,26 +16,24 @@ namespace FluidSim.Core public void Step() { - // 1. Volumes publish state + // 1. Volumes publish state to their ports foreach (var v in _volumes) v.PushStateToPort(); - // 2. Compute boundary fluxes (orifice model) + // 2. Set volume states as boundary conditions on pipes foreach (var conn in _connections) { if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) - ApplyOrifice(conn, conn.PortA, conn.PortB); + SetVolumeBC(conn.PortA, conn.PortB); else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB)) - ApplyOrifice(conn, conn.PortB, conn.PortA); - else if (IsVolumePort(conn.PortA) && IsVolumePort(conn.PortB)) - VolumeToVolume(conn); + SetVolumeBC(conn.PortB, conn.PortA); } - // 3. Pipe simulation step + // 3. Run pipe simulations foreach (var p in _pipes) p.Simulate(); - // 4. Transfer pipe‑port data to volumes + // 4. Transfer pipe‑port flows to volume ports foreach (var conn in _connections) { if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) @@ -53,64 +51,29 @@ namespace FluidSim.Core bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p); Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p); - void ApplyOrifice(Connection conn, Port pipePort, Port volPort) + void SetVolumeBC(Port pipePort, Port volPort) { Pipe1D pipe = GetPipe(pipePort); if (pipe == null) return; bool isLeft = pipe.PortA == pipePort; - double pP = isLeft ? pipe.GetLeftPressure() : pipe.GetRightPressure(); - double rhoP = isLeft ? pipe.GetLeftDensity() : pipe.GetRightDensity(); - 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); - if (isLeft) - pipe.SetLeftBoundaryFlux(massFlux, momFlux, energyFlux); + pipe.SetLeftVolumeState(volPort.Density, volPort.Pressure); else - 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); - - 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; + pipe.SetRightVolumeState(volPort.Density, volPort.Pressure); } 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; if (mdot < 0) // pipe → volume { - // ★ pipePort.SpecificEnthalpy now contains TOTAL enthalpy + // pipePort.SpecificEnthalpy is already total (h + ½u²) volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy; } - // else: fluid goes volume → pipe → volume owns its own (static) enthalpy, - // which is already correct. + // else: volume → pipe, volume’s own static enthalpy is used (already set) } } } \ No newline at end of file