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; // Boundary fluxes (set by solver before each step) private double _fxL_mass, _fxL_mom, _fxL_ener; private double _fxR_mass, _fxR_mom, _fxR_ener; private bool _leftSet, _rightSet; 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); public Pipe1D(double length, double area, int nCells, int sampleRate) { _n = nCells; _dx = length / nCells; _dt = 1.0 / sampleRate; _area = area; _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 SetLeftBoundaryFlux(double m, double p, double e) { _fxL_mass = m; _fxL_mom = p; _fxL_ener = e; _leftSet = true; } public void SetRightBoundaryFlux(double m, double p, double e) { _fxR_mass = m; _fxR_mom = p; _fxR_ener = e; _rightSet = true; } public void Simulate() { int n = _n; double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1]; // Left face if (_leftSet) { Fm[0] = _fxL_mass; Fp[0] = _fxL_mom; Fe[0] = _fxL_ener; } else { Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0; } // reflective wall // Internal faces (HLLC) 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 face if (_rightSet) { Fm[n] = _fxR_mass; Fp[n] = _fxR_mom; Fe[n] = _fxR_ener; } else { Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0; } // Update cells 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; // Clamp to physical if (_rho[i] < 1e-12) _rho[i] = 1e-12; double u = _rhou[i] / _rho[i]; double kinetic = 0.5 * _rho[i] * u * u; if (_E[i] < kinetic) _E[i] = kinetic; } // Friction (energy‑conserving) if (FrictionFactor > 0) { double D = 2.0 * Math.Sqrt(_area / Math.PI); for (int i = 0; i < _n; i++) { double u = _rhou[i] / Math.Max(_rho[i], 1e-12); double f = FrictionFactor / (2.0 * D) * _rho[i] * u * Math.Abs(u); //_rhou[i] -= _dt * f; FRICTIN DISABLED!!! } } // Write port flows for the solver PortA.MassFlowRate = _leftSet ? _fxL_mass * _area : 0.0; PortB.MassFlowRate = _rightSet ? -_fxR_mass * _area : 0.0; // Enthalpy for upwinding PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0); PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1); // Reset for next step _leftSet = _rightSet = false; } 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; } 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) { 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; } } } }