From 3926ed7ef9a526662d0bda656a746bfcce33a8b5 Mon Sep 17 00:00:00 2001 From: max Date: Sun, 3 May 2026 00:20:17 +0200 Subject: [PATCH] Introduced automatic sub stepping and pipe cell count --- Components/Pipe1D.cs | 242 ++++++++++++++++++++++++++++--------------- Core/Simulation.cs | 24 ++++- 2 files changed, 181 insertions(+), 85 deletions(-) diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs index 5b6e2aa..a77b885 100644 --- a/Components/Pipe1D.cs +++ b/Components/Pipe1D.cs @@ -9,14 +9,19 @@ namespace FluidSim.Components public Port PortB { get; } public double Area => _area; - private int _n; - private double _dx, _dt, _gamma = 1.4, _area; + private int _n; // number of cells + private double _dx, _dt, _gamma, _area; private double[] _rho, _rhou, _E; - // Volume states at boundaries - private double _rhoLeft, _pLeft, _rhoRight, _pRight; + // Volume boundary states, constant during sub‑steps + private double _rhoLeft, _pLeft; + private double _rhoRight, _pRight; private bool _leftBCSet, _rightBCSet; + // CFL control + private const double CflTarget = 0.8; + private const double ReferenceSoundSpeed = 340.0; // m/s, standard air + public double FrictionFactor { get; set; } = 0.02; public int GetCellCount() => _n; @@ -24,12 +29,30 @@ namespace FluidSim.Components public double GetCellPressure(int i) => Pressure(i); public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12); - public Pipe1D(double length, double area, int nCells, int sampleRate) + /// + /// Creates a 1D pipe. + /// Cell count is automatically determined to satisfy CFL in still air. + /// + /// Pipe length in metres. + /// Cross‑sectional area in m². + /// Global simulation sample rate (Hz). + public Pipe1D(double length, double area, int sampleRate) { + // Desired spatial step to keep CFL ≤ target for still air + double dtGlobal = 1.0 / sampleRate; + double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget; + + // Number of cells must be at least 2; try to hit dxTarget + int nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero)); + // Ensure we don't accidentally overshoot dxTarget by more than a factor + while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1) + nCells++; + _n = nCells; - _dx = length / nCells; - _dt = 1.0 / sampleRate; + _dx = length / _n; + _dt = dtGlobal; // global (audio) time step _area = area; + _gamma = 1.4; _rho = new double[_n]; _rhou = new double[_n]; @@ -56,7 +79,6 @@ namespace FluidSim.Components public double GetLeftDensity() => _rho[0]; public double GetRightDensity() => _rho[_n - 1]; - // ★ New: pass both density and pressure from the volume public void SetLeftVolumeState(double rhoVol, double pVol) { _rhoLeft = rhoVol; @@ -80,95 +102,152 @@ namespace FluidSim.Components return h + 0.5 * u * u; } + /// + /// Advance the pipe over one global time step using sub‑stepping. + /// Must be called once per global simulation cycle. + /// 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 --- + // --- Determine maximum wave speed in the pipe --- + double maxWaveSpeed = 0.0; 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; + double rho = Math.Max(_rho[i], 1e-12); + double u = Math.Abs(_rhou[i] / rho); + double c = Math.Sqrt(_gamma * Pressure(i) / rho); + double local = u + c; + if (local > maxWaveSpeed) maxWaveSpeed = local; } - // --- Friction disabled --- - // if (FrictionFactor > 0) { … } + if (maxWaveSpeed < 1e-8) maxWaveSpeed = 1e-8; - // --- Port flows --- - PortA.MassFlowRate = _leftBCSet ? Fm[0] * _area : 0.0; - PortB.MassFlowRate = _rightBCSet ? -Fm[n] * _area : 0.0; + int nSub = Math.Max(1, (int)Math.Ceiling(_dt * maxWaveSpeed / (CflTarget * _dx))); + double dtSub = _dt / nSub; - PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0); - PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1); + // Accumulators for net mass flows + double sumMdotA = 0.0, sumMdotB = 0.0; + + // Accumulators for fluid that ENTERS the volumes (pipe → volume) + double massInA = 0.0, energyInA = 0.0; + double massInB = 0.0, energyInB = 0.0; + + for (int step = 0; step < nSub; step++) + { + double[] Fm = new double[n + 1]; + double[] Fp = new double[n + 1]; + double[] Fe = new double[n + 1]; + + // Left boundary (face 0) + if (_leftBCSet) + { + HLLCFlux(_rhoLeft, 0.0, _pLeft, + _rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), Pressure(0), + 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); + HLLCFlux(rhoL, uL, pL, + _rhoRight, 0.0, _pRight, + 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] -= dtSub * dM; + _rhou[i] -= dtSub * dP; + _E[i] -= dtSub * 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; + } + + // Sub‑step mass flow rates (kg/s) + double mdotA_sub = _leftBCSet ? Fm[0] * _area : 0.0; // >0 = into pipe + double mdotB_sub = _rightBCSet ? -Fm[n] * _area : 0.0; // >0 = into pipe from right + + sumMdotA += mdotA_sub; + sumMdotB += mdotB_sub; + + // Flow FROM pipe INTO volume A: mdotA_sub < 0 + if (mdotA_sub < 0 && _leftBCSet) + { + double massRate = -mdotA_sub; // kg/s entering volume A + double h = GetCellTotalSpecificEnthalpy(0); + massInA += massRate * dtSub; + energyInA += massRate * dtSub * h; + } + + // Flow FROM pipe INTO volume B: mdotB_sub < 0 (because + // mdotB_sub = -Fm[n], and Fm[n] > 0 is flow to the right) + if (mdotB_sub < 0 && _rightBCSet) + { + double massRate = -mdotB_sub; // kg/s entering volume B + double h = GetCellTotalSpecificEnthalpy(_n - 1); + massInB += massRate * dtSub; + energyInB += massRate * dtSub * h; + } + } + + // Averaged net mass flows (sign: positive = into pipe) + PortA.MassFlowRate = sumMdotA / nSub; + PortB.MassFlowRate = sumMdotB / nSub; + + // Assign enthalpy ONLY for the fluid that physically entered the volume + if (massInA > 1e-12) + PortA.SpecificEnthalpy = energyInA / massInA; + + if (massInB > 1e-12) + PortB.SpecificEnthalpy = energyInB / massInB; + + // If no inflow occurred, leave the port’s enthalpy unchanged. + // (It will be set to the volume’s static enthalpy by PushStateToPort + // or overwritten by TransferPipeToVolume if flow reverses later.) _leftBCSet = _rightBCSet = false; } - double Pressure(int i) => + // Pressure and HLLC flux unchanged + private double Pressure(int i) => (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12)); - void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR, - out double fm, out double fp, out double fe) + 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)); double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12)); @@ -176,6 +255,7 @@ 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)); diff --git a/Core/Simulation.cs b/Core/Simulation.cs index 2b310c5..1e8d488 100644 --- a/Core/Simulation.cs +++ b/Core/Simulation.cs @@ -19,15 +19,14 @@ namespace FluidSim.Core dt = 1.0 / sampleRate; double V = 5.0 * Units.L; - volA = new Volume0D(V, 1.1 * Units.atm, Units.Celsius(20), sampleRate); + volA = new Volume0D(V, 2.0 * 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 = new Pipe1D(length, area, sampleRate); pipe.SetUniformState(volA.Density, 0.0, volA.Pressure); pipe.FrictionFactor = 0.02; @@ -54,13 +53,30 @@ namespace FluidSim.Core public static void Log() { - if (stepCount <= 50 || stepCount % 200 == 0) + bool logPipe = true; + + if ((stepCount <= 10 || (stepCount <= 1000 && stepCount % 100 == 0)) || stepCount % 1000 == 0 && stepCount < 10000) { + // Summary line 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"); + + // Per‑cell state + if (logPipe && stepCount <= 1000) + { + int n = pipe.GetCellCount(); + for (int i = 0; i < n; i++) + { + double rho = pipe.GetCellDensity(i); + double p = pipe.GetCellPressure(i); + double u = pipe.GetCellVelocity(i); + Console.WriteLine( + $" Cell {i,2}: ρ={rho,8:F4} kg/m³, p={p,10:F2} Pa, u={u,8:F3} m/s"); + } + } } } }