using System; 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; public double DampingMultiplier { get; set; } = 1.0; private int _n; private double _dx, _dt, _gamma, _area, _diameter; private double[] _rho, _rhou, _E; // Volume‑coupling ghost states for boundaries A and B private double _rhoA, _pA; private double _rhoB, _pB; private bool _aBCSet, _bBCSet; private BoundaryType _aBCType = BoundaryType.VolumeCoupling; private BoundaryType _bBCType = BoundaryType.VolumeCoupling; private double _aAmbientPressure = 101325.0; private double _bAmbientPressure = 101325.0; private const double CflTarget = 0.8; private const double ReferenceSoundSpeed = 340.0; 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 BoundaryType ABCType => _aBCType; public BoundaryType BBCType => _bBCType; 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; _area = area; _gamma = 1.4; // Hydraulic diameter for a circular pipe _diameter = 2.0 * Math.Sqrt(area / Math.PI); _rho = new double[_n]; _rhou = new double[_n]; _E = new double[_n]; PortA = new Port(); PortB = new Port(); } public void SetABoundaryType(BoundaryType type) => _aBCType = type; public void SetBBoundaryType(BoundaryType type) => _bBCType = type; public void SetAAmbientPressure(double p) => _aAmbientPressure = p; public void SetBAmbientPressure(double p) => _bAmbientPressure = p; 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 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 SetAVolumeState(double rhoVol, double pVol) { _rhoA = rhoVol; _pA = pVol; _aBCSet = true; } public void SetBVolumeState(double rhoVol, double pVol) { _rhoB = rhoVol; _pB = pVol; _bBCSet = true; } public void ClearBC() => _aBCSet = _bBCSet = 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]; // ---------- Boundary A (face 0, left) ---------- double rhoIntA = _rho[0]; double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12); double pIntA = Pressure(0); switch (_aBCType) { case BoundaryType.VolumeCoupling: if (_aBCSet) { HLLCFlux(_rhoA, 0.0, _pA, rhoIntA, uIntA, pIntA, out Fm[0], out Fp[0], out Fe[0]); } else { Fm[0] = 0; Fp[0] = pIntA; Fe[0] = 0; } break; case BoundaryType.OpenEnd: OpenEndFluxA(rhoIntA, uIntA, pIntA, _aAmbientPressure, out Fm[0], out Fp[0], out Fe[0]); break; case BoundaryType.ClosedEnd: ClosedEndFlux(rhoIntA, uIntA, pIntA, isRightBoundary: false, out Fm[0], out Fp[0], out Fe[0]); break; } // ---------- Internal faces ---------- for (int i = 0; i < n - 1; i++) { double rhoL = _rho[i]; double uL = _rhou[i] / Math.Max(rhoL, 1e-12); double pL = Pressure(i); double rhoR = _rho[i + 1]; double uR = _rhou[i + 1] / Math.Max(rhoR, 1e-12); double pR = Pressure(i + 1); HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]); } // ---------- Boundary B (face n, right) ---------- double rhoIntB = _rho[n - 1]; double uIntB = _rhou[n - 1] / Math.Max(rhoIntB, 1e-12); double pIntB = Pressure(n - 1); switch (_bBCType) { case BoundaryType.VolumeCoupling: if (_bBCSet) { HLLCFlux(rhoIntB, uIntB, pIntB, _rhoB, 0.0, _pB, out Fm[n], out Fp[n], out Fe[n]); } else { Fm[n] = 0; Fp[n] = pIntB; Fe[n] = 0; } break; case BoundaryType.OpenEnd: OpenEndFluxB(rhoIntB, uIntB, pIntB, _bAmbientPressure, out Fm[n], out Fp[n], out Fe[n]); break; case BoundaryType.ClosedEnd: ClosedEndFlux(rhoIntB, uIntB, pIntB, isRightBoundary: true, out Fm[n], out Fp[n], out Fe[n]); break; } // ---- Cell update with linear laminar damping ---- double radius = _diameter / 2.0; double mu_air = 1.8e-5; double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius); 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; double rho = Math.Max(_rho[i], 1e-12); double dampingRate = laminarCoeff / rho; double dampingFactor = Math.Exp(-dampingRate * dtSub); _rhou[i] *= dampingFactor; 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 ---------- double mdotA_sub = _aBCType == BoundaryType.VolumeCoupling && _aBCSet ? Fm[0] * _area : 0.0; double mdotB_sub = _bBCType == BoundaryType.VolumeCoupling && _bBCSet ? -Fm[n] * _area : 0.0; PortA.MassFlowRate = mdotA_sub; PortB.MassFlowRate = mdotB_sub; PortA.Pressure = pIntA; PortB.Pressure = pIntB; PortA.Density = _rho[0]; PortB.Density = _rho[n - 1]; // Corrected enthalpy for both directions if (_aBCType == BoundaryType.VolumeCoupling && _aBCSet) { PortA.SpecificEnthalpy = mdotA_sub < 0 ? GetCellTotalSpecificEnthalpy(0) : (_gamma / (_gamma - 1.0)) * _pA / Math.Max(_rhoA, 1e-12); } if (_bBCType == BoundaryType.VolumeCoupling && _bBCSet) { PortB.SpecificEnthalpy = mdotB_sub < 0 ? GetCellTotalSpecificEnthalpy(_n - 1) : (_gamma / (_gamma - 1.0)) * _pB / Math.Max(_rhoB, 1e-12); } } 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; } private double Pressure(int i) => (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12)); // ========== Characteristic‑based Open End ========== private void OpenEndFluxA(double rhoInt, double uInt, double pInt, double pAmb, out double fm, out double fp, out double fe) { double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); // Subsonic inflow (uInt ≤ 0, so flow inside pipe ←) if (uInt <= -cInt) // supersonic inflow – use interior state as ghost { fm = rhoInt * uInt; fp = rhoInt * uInt * uInt + pInt; fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt; return; } else if (uInt <= 0) // subsonic inflow { // Reservoir condition: p = pAmb, T = 300K, u = 0 double T0 = 300.0; double R = 287.0; double rhoGhost = pAmb / (R * T0); HLLCFlux(rhoGhost, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); return; } else // subsonic outflow (uInt > 0) { // Ghost pressure forced to pAmb double s = pInt / Math.Pow(rhoInt, _gamma); double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma); double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost); // Outgoing Riemann invariant J⁻ = uInt - 2*cInt/(γ-1) (for left boundary) double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0); double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0); // Prevent spurious inflow by clipping to zero if (uGhost < 0) uGhost = 0; HLLCFlux(rhoGhost, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); } } private void OpenEndFluxB(double rhoInt, double uInt, double pInt, double pAmb, out double fm, out double fp, out double fe) { double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); if (uInt >= cInt) // supersonic outflow (extrapolation) { fm = rhoInt * uInt; fp = rhoInt * uInt * uInt + pInt; fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt; return; } else if (uInt >= 0) // subsonic outflow { double s = pInt / Math.Pow(rhoInt, _gamma); double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma); double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost); // Outgoing Riemann invariant J⁺ = uInt + 2*cInt/(γ-1) (for right boundary) double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0); double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0); // Clip to zero to prevent inflow if (uGhost > 0) uGhost = 0; HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pAmb, out fm, out fp, out fe); } else // subsonic inflow { double T0 = 300.0; double R = 287.0; double rhoGhost = pAmb / (R * T0); HLLCFlux(rhoInt, uInt, pInt, rhoGhost, 0.0, pAmb, out fm, out fp, out fe); } } // ========== Closed end (mirror) ========== private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRightBoundary, out double fm, out double fp, out double fe) { double rhoGhost = rhoInt; double pGhost = pInt; double uGhost = -uInt; // mirror velocity if (isRightBoundary) HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe); else HLLCFlux(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe); } // ========== Standard HLLC flux ========== 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; } } public double GetPressureAtFraction(double fraction) { int i = (int)(fraction * (_n - 1)); i = Math.Clamp(i, 0, _n - 1); return Pressure(i); } } }