366 lines
15 KiB
C#
366 lines
15 KiB
C#
using System;
|
||
using FluidSim.Interfaces;
|
||
|
||
namespace FluidSim.Components
|
||
{
|
||
public enum BoundaryType
|
||
{
|
||
OpenEnd,
|
||
ZeroPressureOpen, // pressure‑release boundary (strong reflection)
|
||
ClosedEnd,
|
||
GhostCell
|
||
}
|
||
|
||
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;
|
||
|
||
// Ghost cell states
|
||
private double _rhoGhostL, _uGhostL, _pGhostL;
|
||
private double _rhoGhostR, _uGhostR, _pGhostR;
|
||
private bool _ghostLSet, _ghostRSet;
|
||
|
||
private BoundaryType _aBCType = BoundaryType.GhostCell;
|
||
private BoundaryType _bBCType = BoundaryType.GhostCell;
|
||
|
||
private double _aAmbientPressure = 101325.0;
|
||
private double _bAmbientPressure = 101325.0;
|
||
|
||
private const double CflTarget = 0.8;
|
||
private const double ReferenceSoundSpeed = 340.0;
|
||
private double _lastMassFlow = 0.0;
|
||
|
||
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;
|
||
_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 SetGhostLeft(double rho, double u, double p)
|
||
{
|
||
_rhoGhostL = rho;
|
||
_uGhostL = u;
|
||
_pGhostL = p;
|
||
_ghostLSet = true;
|
||
}
|
||
public void SetGhostRight(double rho, double u, double p)
|
||
{
|
||
_rhoGhostR = rho;
|
||
_uGhostR = u;
|
||
_pGhostR = p;
|
||
_ghostRSet = true;
|
||
}
|
||
public void ClearGhostFlag()
|
||
{
|
||
_ghostLSet = false;
|
||
_ghostRSet = false;
|
||
}
|
||
|
||
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 GetOpenEndMassFlow()
|
||
{
|
||
if (_bBCType != BoundaryType.OpenEnd && _bBCType != BoundaryType.ZeroPressureOpen)
|
||
return 0.0;
|
||
|
||
int lastCell = _n - 1;
|
||
double rho = _rho[lastCell];
|
||
double u = _rhou[lastCell] / Math.Max(rho, 1e-12);
|
||
double p = Pressure(lastCell);
|
||
|
||
double c = Math.Sqrt(_gamma * p / rho);
|
||
double uFace = u;
|
||
double rhoFace = rho;
|
||
double pFace = p;
|
||
|
||
// Subsonic outflow: impose ambient pressure, adjust velocity and density
|
||
if (uFace > 0 && uFace < c)
|
||
{
|
||
double s = p / Math.Pow(rho, _gamma);
|
||
double rhoAmb = Math.Pow(_bAmbientPressure / s, 1.0 / _gamma);
|
||
double cAmb = Math.Sqrt(_gamma * _bAmbientPressure / rhoAmb);
|
||
double J_plus = u + 2.0 * c / (_gamma - 1.0);
|
||
double uFaceNew = J_plus - 2.0 * cAmb / (_gamma - 1.0);
|
||
if (uFaceNew > 0) uFace = uFaceNew;
|
||
if (uFace < 0) uFace = 0;
|
||
rhoFace = rhoAmb;
|
||
pFace = _bAmbientPressure;
|
||
}
|
||
|
||
return rhoFace * uFace * _area;
|
||
}
|
||
|
||
public double GetAndStoreMassFlowForDerivative()
|
||
{
|
||
double current = GetOpenEndMassFlow();
|
||
double derivative = (current - _lastMassFlow) / _dt;
|
||
_lastMassFlow = current;
|
||
return derivative;
|
||
}
|
||
|
||
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 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);
|
||
private double Pressure(int i) => (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
|
||
public double GetPressureAtFraction(double fraction)
|
||
{
|
||
int i = (int)(fraction * (_n - 1));
|
||
i = Math.Clamp(i, 0, _n - 1);
|
||
return Pressure(i);
|
||
}
|
||
|
||
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];
|
||
|
||
// Left boundary (face 0)
|
||
double rhoL = _rho[0];
|
||
double uL = _rhou[0] / Math.Max(rhoL, 1e-12);
|
||
double pL = Pressure(0);
|
||
if (_aBCType == BoundaryType.GhostCell && _ghostLSet)
|
||
HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, rhoL, uL, pL, out Fm[0], out Fp[0], out Fe[0]);
|
||
else if (_aBCType == BoundaryType.OpenEnd)
|
||
OpenEndFluxLeft(rhoL, uL, pL, _aAmbientPressure, out Fm[0], out Fp[0], out Fe[0]);
|
||
else if (_aBCType == BoundaryType.ZeroPressureOpen)
|
||
{
|
||
// Strong reflection: force pressure to ambient, extrapolate density and velocity
|
||
double rhoFace = rhoL;
|
||
double uFace = uL;
|
||
double pFace = _aAmbientPressure;
|
||
HLLCFlux(rhoFace, uFace, pFace, rhoL, uL, pL, out Fm[0], out Fp[0], out Fe[0]);
|
||
}
|
||
else if (_aBCType == BoundaryType.ClosedEnd)
|
||
ClosedEndFlux(rhoL, uL, pL, false, out Fm[0], out Fp[0], out Fe[0]);
|
||
else
|
||
{ Fm[0] = 0; Fp[0] = pL; Fe[0] = 0; }
|
||
|
||
// Internal faces
|
||
for (int i = 0; i < n - 1; i++)
|
||
{
|
||
double rhoLi = _rho[i];
|
||
double uLi = _rhou[i] / Math.Max(rhoLi, 1e-12);
|
||
double pLi = Pressure(i);
|
||
double rhoRi = _rho[i + 1];
|
||
double uRi = _rhou[i + 1] / Math.Max(rhoRi, 1e-12);
|
||
double pRi = Pressure(i + 1);
|
||
HLLCFlux(rhoLi, uLi, pLi, rhoRi, uRi, pRi, out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
|
||
}
|
||
|
||
// Right boundary (face n)
|
||
double rhoR = _rho[n - 1];
|
||
double uR = _rhou[n - 1] / Math.Max(rhoR, 1e-12);
|
||
double pR = Pressure(n - 1);
|
||
if (_bBCType == BoundaryType.GhostCell && _ghostRSet)
|
||
HLLCFlux(rhoR, uR, pR, _rhoGhostR, _uGhostR, _pGhostR, out Fm[n], out Fp[n], out Fe[n]);
|
||
else if (_bBCType == BoundaryType.OpenEnd)
|
||
OpenEndFluxRight(rhoR, uR, pR, _bAmbientPressure, out Fm[n], out Fp[n], out Fe[n]);
|
||
else if (_bBCType == BoundaryType.ZeroPressureOpen)
|
||
{
|
||
double rhoFace = rhoR;
|
||
double uFace = uR;
|
||
double pFace = _bAmbientPressure;
|
||
HLLCFlux(rhoR, uR, pR, rhoFace, uFace, pFace, out Fm[n], out Fp[n], out Fe[n]);
|
||
}
|
||
else if (_bBCType == BoundaryType.ClosedEnd)
|
||
ClosedEndFlux(rhoR, uR, pR, true, out Fm[n], out Fp[n], out Fe[n]);
|
||
else
|
||
{ Fm[n] = 0; Fp[n] = pR; Fe[n] = 0; }
|
||
|
||
// Cell update (linear 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++)
|
||
{
|
||
_rho[i] -= dtSub * (Fm[i + 1] - Fm[i]) / _dx;
|
||
_rhou[i] -= dtSub * (Fp[i + 1] - Fp[i]) / _dx;
|
||
_E[i] -= dtSub * (Fe[i + 1] - Fe[i]) / _dx;
|
||
|
||
double rho = Math.Max(_rho[i], 1e-12);
|
||
double dampingFactor = Math.Exp(-(laminarCoeff / rho) * 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;
|
||
}
|
||
}
|
||
|
||
// ---------- HLLC Riemann solver ----------
|
||
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 = pR + rR * (SR - uR) * (Ss - uR);
|
||
double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
|
||
fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
|
||
}
|
||
}
|
||
|
||
// ---------- Characteristic open‑end boundaries ----------
|
||
private void OpenEndFluxLeft(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 inflow
|
||
{
|
||
fm = rhoInt * uInt;
|
||
fp = rhoInt * uInt * uInt + pInt;
|
||
fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
|
||
return;
|
||
}
|
||
if (uInt <= 0) // subsonic inflow
|
||
{
|
||
double T0 = 300.0, R = 287.0;
|
||
double ghost_Rho = pAmb / (R * T0);
|
||
HLLCFlux(ghost_Rho, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||
return;
|
||
}
|
||
// subsonic outflow
|
||
double s = pInt / Math.Pow(rhoInt, _gamma);
|
||
double ghostRho = Math.Pow(pAmb / s, 1.0 / _gamma);
|
||
double cGhost = Math.Sqrt(_gamma * pAmb / ghostRho);
|
||
double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0);
|
||
double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0);
|
||
if (uGhost < 0) uGhost = 0;
|
||
HLLCFlux(ghostRho, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||
}
|
||
|
||
private void OpenEndFluxRight(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
|
||
{
|
||
fm = rhoInt * uInt;
|
||
fp = rhoInt * uInt * uInt + pInt;
|
||
fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
|
||
return;
|
||
}
|
||
if (uInt >= 0) // subsonic outflow
|
||
{
|
||
double s = pInt / Math.Pow(rhoInt, _gamma);
|
||
double ghost_Rho = Math.Pow(pAmb / s, 1.0 / _gamma);
|
||
double cGhost = Math.Sqrt(_gamma * pAmb / ghost_Rho);
|
||
double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0);
|
||
double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0);
|
||
if (uGhost > 0) uGhost = 0;
|
||
HLLCFlux(rhoInt, uInt, pInt, ghost_Rho, uGhost, pAmb, out fm, out fp, out fe);
|
||
}
|
||
// subsonic inflow
|
||
double T0 = 300.0, R = 287.0;
|
||
double ghostRho = pAmb / (R * T0);
|
||
HLLCFlux(rhoInt, uInt, pInt, ghostRho, 0.0, pAmb, out fm, out fp, out fe);
|
||
}
|
||
|
||
private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRight,
|
||
out double fm, out double fp, out double fe)
|
||
{
|
||
double rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt;
|
||
if (isRight)
|
||
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);
|
||
}
|
||
}
|
||
} |