grillkol/FluidSim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added boundary states for correct resonances

Browse Source
This commit is contained in:
max
2026-05-03 01:52:55 +02:00
parent 3926ed7ef9
commit a006a07049
9 changed files with 432 additions and 244 deletions
Download Patch File Download Diff File Expand all files Collapse all files

379
Components/Pipe1D.cs
View File

@@ -3,24 +3,35 @@ using FluidSim.Interfaces;
namespace FluidSim.Components
{
public enum BoundaryType
{
VolumeCoupling,
OpenEnd,
ClosedEnd
}
public class Pipe1D
{
public Port PortA { get; }
public Port PortB { get; }
public double Area => _area;
private int _n; // number of cells
private int _n;
private double _dx, _dt, _gamma, _area;
private double[] _rho, _rhou, _E;
// Volume boundary states, constant during sub‑steps
private double _rhoLeft, _pLeft;
private double _rhoRight, _pRight;
private bool _leftBCSet, _rightBCSet;
// CFL control
private BoundaryType _leftBCType = BoundaryType.VolumeCoupling;
private BoundaryType _rightBCType = BoundaryType.VolumeCoupling;
private double _leftAmbientPressure = 101325.0;
private double _rightAmbientPressure = 101325.0;
private const double CflTarget = 0.8;
private const double ReferenceSoundSpeed = 340.0; // m/s, standard air
private const double ReferenceSoundSpeed = 340.0;
public double FrictionFactor { get; set; } = 0.02;
@@ -29,28 +40,29 @@ namespace FluidSim.Components
public double GetCellPressure(int i) => Pressure(i);
public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);
/// <summary>
/// Creates a 1D pipe.
/// Cell count is automatically determined to satisfy CFL in still air.
/// </summary>
/// <param name="length">Pipe length in metres.</param>
/// <param name="area">Cross‑sectional area in m².</param>
/// <param name="sampleRate">Global simulation sample rate (Hz).</param>
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;
public BoundaryType LeftBCType => _leftBCType;
public BoundaryType RightBCType => _rightBCType;
// 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++;
public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
{
double dtGlobal = 1.0 / sampleRate;
int nCells;
if (forcedCellCount > 1)
{
nCells = forcedCellCount;
}
else
{
double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget;
nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero));
while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1)
nCells++;
}
_n = nCells;
_dx = length / _n;
_dt = dtGlobal; // global (audio) time step
_dt = dtGlobal;
_area = area;
_gamma = 1.4;
@@ -62,6 +74,11 @@ 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 SetUniformState(double rho, double u, double p)
{
double e = p / ((_gamma - 1) * rho);
@@ -74,10 +91,14 @@ namespace FluidSim.Components
}
}
public double GetLeftPressure() => Pressure(0);
public double GetRightPressure() => Pressure(_n - 1);
public double GetLeftDensity() => _rho[0];
public double GetRightDensity() => _rho[_n - 1];
public void SetCellState(int i, double rho, double u, double p)
{
if (i < 0 || i >= _n) return;
_rho[i] = rho;
_rhou[i] = rho * u;
double e = p / ((_gamma - 1) * rho);
_E[i] = rho * e + 0.5 * rho * u * u;
}
public void SetLeftVolumeState(double rhoVol, double pVol)
{
@@ -93,6 +114,165 @@ namespace FluidSim.Components
_rightBCSet = true;
}
public void ClearBC() => _leftBCSet = _rightBCSet = false;
public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
{
double maxW = 0.0;
for (int i = 0; i < _n; i++)
{
double rho = _rho[i];
double u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12));
double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12));
double local = u + c;
if (local > maxW) maxW = local;
}
maxW = Math.Max(maxW, 1e-8);
return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
}
public void SimulateSingleStep(double dtSub)
{
int n = _n;
double[] Fm = new double[n + 1];
double[] Fp = new double[n + 1];
double[] Fe = new double[n + 1];
// Left boundary (face 0)
switch (_leftBCType)
{
case BoundaryType.VolumeCoupling:
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;
}
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,
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]);
}
break;
}
// 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)
switch (_rightBCType)
{
case BoundaryType.VolumeCoupling:
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;
}
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,
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]);
}
break;
}
// 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];
double pMin = 100.0;
double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
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;
PortA.MassFlowRate = mdotA_sub;
PortB.MassFlowRate = mdotB_sub;
PortA.Pressure = Pressure(0);
PortB.Pressure = Pressure(_n - 1);
PortA.Density = _rho[0];
PortB.Density = _rho[_n - 1];
if (_leftBCType == BoundaryType.VolumeCoupling && _leftBCSet)
{
PortA.SpecificEnthalpy = mdotA_sub < 0
? GetCellTotalSpecificEnthalpy(0)
: 0.0;
}
if (_rightBCType == BoundaryType.VolumeCoupling && _rightBCSet)
{
PortB.SpecificEnthalpy = mdotB_sub < 0
? GetCellTotalSpecificEnthalpy(_n - 1)
: 0.0;
}
}
private double GetCellTotalSpecificEnthalpy(int i)
{
double rho = Math.Max(_rho[i], 1e-12);
@@ -102,146 +282,6 @@ namespace FluidSim.Components
return h + 0.5 * u * u;
}
/// <summary>
/// Advance the pipe over one global time step using sub‑stepping.
/// Must be called once per global simulation cycle.
/// </summary>
public void Simulate()
{
int n = _n;
// --- Determine maximum wave speed in the pipe ---
double maxWaveSpeed = 0.0;
for (int i = 0; i < n; i++)
{
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;
}
if (maxWaveSpeed < 1e-8) maxWaveSpeed = 1e-8;
int nSub = Math.Max(1, (int)Math.Ceiling(_dt * maxWaveSpeed / (CflTarget * _dx)));
double dtSub = _dt / nSub;
// 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;
}
// 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));
@@ -279,5 +319,12 @@ namespace FluidSim.Components
fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
}
}
public double GetPressureAtFraction(double fraction)
{
int i = (int)(fraction * (_n - 1));
i = Math.Clamp(i, 0, _n - 1);
return Pressure(i);
}
}
}