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");
+ }
+ }
}
}
}