420 lines
16 KiB
C#
420 lines
16 KiB
C#
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);
|
||
}
|
||
}
|
||
} |