283 lines
11 KiB
C#
283 lines
11 KiB
C#
using System;
|
||
using FluidSim.Interfaces;
|
||
|
||
namespace FluidSim.Components
|
||
{
|
||
public class Pipe1D
|
||
{
|
||
public Port PortA { get; }
|
||
public Port PortB { get; }
|
||
public double Area => _area;
|
||
|
||
private int _n; // number of cells
|
||
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 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;
|
||
public double GetCellDensity(int i) => _rho[i];
|
||
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;
|
||
|
||
// 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 / _n;
|
||
_dt = dtGlobal; // global (audio) time step
|
||
_area = area;
|
||
_gamma = 1.4;
|
||
|
||
_rho = new double[_n];
|
||
_rhou = new double[_n];
|
||
_E = new double[_n];
|
||
|
||
PortA = new Port();
|
||
PortB = new Port();
|
||
}
|
||
|
||
public void SetUniformState(double rho, double u, double p)
|
||
{
|
||
double e = p / ((_gamma - 1) * rho);
|
||
double Etot = rho * e + 0.5 * rho * u * u;
|
||
for (int i = 0; i < _n; i++)
|
||
{
|
||
_rho[i] = rho;
|
||
_rhou[i] = rho * u;
|
||
_E[i] = Etot;
|
||
}
|
||
}
|
||
|
||
public double GetLeftPressure() => Pressure(0);
|
||
public double GetRightPressure() => Pressure(_n - 1);
|
||
public double GetLeftDensity() => _rho[0];
|
||
public double GetRightDensity() => _rho[_n - 1];
|
||
|
||
public void SetLeftVolumeState(double rhoVol, double pVol)
|
||
{
|
||
_rhoLeft = rhoVol;
|
||
_pLeft = pVol;
|
||
_leftBCSet = true;
|
||
}
|
||
|
||
public void SetRightVolumeState(double rhoVol, double pVol)
|
||
{
|
||
_rhoRight = rhoVol;
|
||
_pRight = pVol;
|
||
_rightBCSet = true;
|
||
}
|
||
|
||
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;
|
||
}
|
||
|
||
/// <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));
|
||
|
||
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));
|
||
double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL;
|
||
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));
|
||
|
||
double FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL;
|
||
double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * uR;
|
||
|
||
if (SL >= 0) { fm = FrL_m; fp = FrL_p; fe = FrL_e; }
|
||
else if (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
|
||
else if (Ss >= 0)
|
||
{
|
||
double rsL = rL * (SL - uL) / (SL - Ss);
|
||
double ps = pL + rL * (SL - uL) * (Ss - uL);
|
||
double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
|
||
fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
|
||
}
|
||
else
|
||
{
|
||
double rsR = rR * (SR - uR) / (SR - Ss);
|
||
double ps = pL + rL * (SL - uL) * (Ss - uL);
|
||
double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
|
||
fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
|
||
}
|
||
}
|
||
}
|
||
} |