250 lines
9.5 KiB
C#
250 lines
9.5 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;
|
||
private double _dx, _dt, _gamma = 1.4, _area;
|
||
private double[] _rho, _rhou, _E;
|
||
private double _hydraulicDiameter;
|
||
|
||
private double _rhoLeft, _pLeft, _rhoRight, _pRight;
|
||
private bool _leftBCSet, _rightBCSet;
|
||
|
||
public double FrictionFactor { get; set; }
|
||
|
||
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>
|
||
/// Create a pipe with CFL‑stable automatic cell count.
|
||
/// </summary>
|
||
/// <param name="length">Pipe length [m].</param>
|
||
/// <param name="area">Cross‑sectional area [m²].</param>
|
||
/// <param name="sampleRate">Simulation step rate [Hz].</param>
|
||
/// <param name="c0">Speed of sound [m/s] (default 343).</param>
|
||
/// <param name="frictionFactor">Darcy friction factor (0 = inviscid).</param>
|
||
/// <param name="cflSafety">CFL safety factor ≤ 1 (0.8 recommended).</param>
|
||
public Pipe1D(double length, double area, int sampleRate,
|
||
double c0 = 343.0, double frictionFactor = 0.02,
|
||
double cflSafety = 0.8)
|
||
{
|
||
if (area <= 0) throw new ArgumentException("Pipe area must be > 0");
|
||
_area = area;
|
||
_dt = 1.0 / sampleRate;
|
||
FrictionFactor = frictionFactor;
|
||
|
||
// Nyquist‑based cell count (wave resolution)
|
||
double nNyquist = Math.Ceiling(length * sampleRate / c0);
|
||
// CFL‑stable cell count: dx ≥ maxSpeed·dt / cflSafety, maxSpeed = 2·c0 (supersonic safe)
|
||
double maxSpeed = 2.0 * c0;
|
||
double dxMinStable = maxSpeed * _dt / cflSafety;
|
||
double nStable = Math.Floor(length / dxMinStable);
|
||
|
||
_n = Math.Max(2, (int)Math.Min(nNyquist, nStable));
|
||
_dx = length / _n;
|
||
|
||
_rho = new double[_n];
|
||
_rhou = new double[_n];
|
||
_E = new double[_n];
|
||
|
||
_hydraulicDiameter = Math.Max(2.0 * Math.Sqrt(_area / Math.PI), 1e-9);
|
||
|
||
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;
|
||
}
|
||
|
||
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)
|
||
{
|
||
double rhoL = _rhoLeft, uL = 0.0, pL = _pLeft;
|
||
double rhoR = _rho[0], uR = _rhou[0] / Math.Max(rhoR, 1e-12), 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], uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12), pL = Pressure(n - 1);
|
||
double rhoR = _rhoRight, uR = 0.0, 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 (inviscid fluxes) ---
|
||
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;
|
||
|
||
// Emergency reset if NaN
|
||
if (double.IsNaN(_rho[i]) || double.IsNaN(_rhou[i]) || double.IsNaN(_E[i]))
|
||
{
|
||
_rho[i] = 1.225; // reset to atmospheric air at 300 K
|
||
_rhou[i] = 0.0;
|
||
_E[i] = 101325.0 / (_gamma - 1.0); // internal energy at 1 atm
|
||
}
|
||
}
|
||
|
||
// --- Friction (Darcy–Weisbach, energy‑conserving) ---
|
||
if (FrictionFactor > 0)
|
||
{
|
||
double D = _hydraulicDiameter;
|
||
double twoD = 2.0 * D;
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
double rho = _rho[i];
|
||
double u = _rhou[i] / rho;
|
||
double absU = Math.Abs(u);
|
||
double src = FrictionFactor * rho * absU * u / twoD;
|
||
double kinOld = 0.5 * rho * u * u;
|
||
_rhou[i] -= _dt * src;
|
||
double uNew = _rhou[i] / rho;
|
||
double kinNew = 0.5 * rho * uNew * uNew;
|
||
_E[i] += (kinOld - kinNew);
|
||
}
|
||
}
|
||
|
||
// --- Publish to ports ---
|
||
PortA.Pressure = Pressure(0);
|
||
PortA.Density = _rho[0];
|
||
PortB.Pressure = Pressure(_n - 1);
|
||
PortB.Density = _rho[_n - 1];
|
||
|
||
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));
|
||
|
||
void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR,
|
||
out double fm, out double fp, out double fe)
|
||
{
|
||
const double eps = 1e-12;
|
||
pL = Math.Max(pL, eps);
|
||
pR = Math.Max(pR, eps);
|
||
|
||
double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, eps));
|
||
double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, eps));
|
||
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 denom = rL * (SL - uL) - rR * (SR - uR);
|
||
double Ss;
|
||
if (Math.Abs(denom) < eps)
|
||
Ss = 0.5 * (uL + uR);
|
||
else
|
||
Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;
|
||
|
||
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 diffSL = SL - uL;
|
||
if (Math.Abs(diffSL) < eps) diffSL = eps;
|
||
double rsL = rL * diffSL / (SL - Ss);
|
||
double ps = pL + rL * diffSL * (Ss - uL);
|
||
double EsL = EL + (Ss - uL) * (Ss + pL / (rL * diffSL));
|
||
fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
|
||
}
|
||
else
|
||
{
|
||
double diffSR = SR - uR;
|
||
if (Math.Abs(diffSR) < eps) diffSR = eps;
|
||
double rsR = rR * diffSR / (SR - Ss);
|
||
double ps = pL + rL * (SL - uL) * (Ss - uL);
|
||
double EsR = ER + (Ss - uR) * (Ss + pR / (rR * diffSR));
|
||
fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
|
||
}
|
||
}
|
||
}
|
||
} |