refactoring (broken right now)
This commit is contained in:
43
Components/Atmosphere.cs
Normal file
43
Components/Atmosphere.cs
Normal file
@@ -0,0 +1,43 @@
|
||||
using FluidSim.Interfaces;
|
||||
|
||||
namespace FluidSim.Components
|
||||
{
|
||||
/// <summary>
|
||||
/// Represents the ambient atmosphere – constant pressure/temperature reservoir.
|
||||
/// </summary>
|
||||
public class Atmosphere : IComponent
|
||||
{
|
||||
public double Pressure { get; set; } = 101325.0;
|
||||
public double Temperature { get; set; } = 300.0;
|
||||
public double GasConstant { get; set; } = 287.0;
|
||||
public double Gamma => 1.4;
|
||||
|
||||
public double Density => Pressure / (GasConstant * Temperature);
|
||||
public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Density;
|
||||
|
||||
public Port Port { get; }
|
||||
|
||||
public Atmosphere()
|
||||
{
|
||||
Port = new Port { Owner = this };
|
||||
UpdatePort();
|
||||
}
|
||||
|
||||
public IReadOnlyList<Port> Ports => new[] { Port };
|
||||
|
||||
public void UpdateState(double dt)
|
||||
{
|
||||
// Atmosphere is static – just ensure the port reflects current values
|
||||
UpdatePort();
|
||||
}
|
||||
|
||||
private void UpdatePort()
|
||||
{
|
||||
Port.Pressure = Pressure;
|
||||
Port.Density = Density;
|
||||
Port.Temperature = Temperature;
|
||||
Port.SpecificEnthalpy = SpecificEnthalpy;
|
||||
// MassFlowRate is set by the solver connector
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -1,249 +0,0 @@
|
||||
using System;
|
||||
using FluidSim.Components;
|
||||
|
||||
namespace FluidSim.Core
|
||||
{
|
||||
public class EngineCylinder
|
||||
{
|
||||
public Volume0D Cylinder { get; private set; }
|
||||
private Crankshaft crankshaft;
|
||||
|
||||
private double bore, stroke, conRodLength, compressionRatio;
|
||||
private double pistonArea;
|
||||
|
||||
public double V_disp { get; private set; }
|
||||
public double V_clear { get; private set; }
|
||||
public bool ignition = false;
|
||||
|
||||
// ---- Exhaust valve ----
|
||||
private double exhMaxOrificeArea;
|
||||
private double exhValveOpenStart = 130.0 * Math.PI / 180.0;
|
||||
private double exhValveOpenEnd = 390.0 * Math.PI / 180.0;
|
||||
private double exhValveRampWidth = 30.0 * Math.PI / 180.0;
|
||||
public double ExhaustOrificeArea => ExhaustValveLift() * exhMaxOrificeArea;
|
||||
public double ExhaustValveLiftCurrent => ExhaustValveLift();
|
||||
|
||||
// ---- Intake valve ----
|
||||
private double intMaxOrificeArea;
|
||||
private double intValveOpenStart = 340.0 * Math.PI / 180.0;
|
||||
private double intValveOpenEnd = 600.0 * Math.PI / 180.0;
|
||||
private double intValveRampWidth = 30.0 * Math.PI / 180.0;
|
||||
public double IntakeOrificeArea => IntakeValveLift() * intMaxOrificeArea;
|
||||
public double IntakeValveLiftCurrent => IntakeValveLift();
|
||||
|
||||
// ---- Combustion ----
|
||||
public double TargetPeakPressure { get; set; } = 50.0 * 101325.0;
|
||||
private const double PeakTemperature = 2500.0;
|
||||
private bool burnInProgress = false;
|
||||
private double burnStartAngle; // cycle angle (0–4π)
|
||||
private double burnDuration = 40.0 * Math.PI / 180.0;
|
||||
private double targetBurnEnergy;
|
||||
private double preIgnitionMass, preIgnitionInternalEnergy;
|
||||
|
||||
private Random rand = new Random();
|
||||
public double MisfireProbability { get; set; } = 0.02;
|
||||
private bool misfireCurrent = false;
|
||||
|
||||
public int CombustionCount { get; private set; }
|
||||
public int MisfireCount { get; private set; }
|
||||
|
||||
// Cycle‑aware angle (0 – 4π)
|
||||
public double CycleAngle => crankshaft.CrankAngle % (4.0 * Math.PI);
|
||||
private double prevCycleAngle;
|
||||
|
||||
// Piston position fraction (0 = TDC, 1 = BDC)
|
||||
public double PistonPositionFraction
|
||||
{
|
||||
get
|
||||
{
|
||||
double currentVol = Cylinder.Volume;
|
||||
if (currentVol <= V_clear) return 0.0;
|
||||
if (currentVol >= V_clear + V_disp) return 1.0;
|
||||
return (currentVol - V_clear) / V_disp;
|
||||
}
|
||||
}
|
||||
|
||||
public EngineCylinder(Crankshaft crankshaft,
|
||||
double bore, double stroke, double compressionRatio,
|
||||
double exhPipeArea, double intPipeArea, int sampleRate)
|
||||
{
|
||||
this.crankshaft = crankshaft;
|
||||
this.bore = bore;
|
||||
this.stroke = stroke;
|
||||
conRodLength = 2.0 * stroke;
|
||||
this.compressionRatio = compressionRatio;
|
||||
exhMaxOrificeArea = exhPipeArea * 0.5;
|
||||
intMaxOrificeArea = intPipeArea * 0.5;
|
||||
pistonArea = Math.PI / 4.0 * bore * bore;
|
||||
|
||||
V_disp = pistonArea * stroke;
|
||||
V_clear = V_disp / (compressionRatio - 1.0);
|
||||
|
||||
// Start at BDC with fresh ambient charge
|
||||
double V_bdc = V_clear + V_disp;
|
||||
double p_amb = 101325.0;
|
||||
double T_amb = 300.0;
|
||||
double rho0 = p_amb / (287.0 * T_amb);
|
||||
double mass0 = rho0 * V_bdc;
|
||||
double energy0 = p_amb * V_bdc / (1.4 - 1.0);
|
||||
|
||||
Cylinder = new Volume0D(V_bdc, p_amb, T_amb, sampleRate)
|
||||
{
|
||||
Gamma = 1.4,
|
||||
GasConstant = 287.0
|
||||
};
|
||||
Cylinder.Volume = V_bdc;
|
||||
Cylinder.Mass = mass0;
|
||||
Cylinder.InternalEnergy = energy0;
|
||||
|
||||
prevCycleAngle = CycleAngle;
|
||||
|
||||
preIgnitionMass = Cylinder.Mass;
|
||||
preIgnitionInternalEnergy = Cylinder.InternalEnergy;
|
||||
}
|
||||
|
||||
// ---- Piston kinematics ----
|
||||
private (double volume, double dvdt) PistonKinematics(double cycleAngle)
|
||||
{
|
||||
double theta = cycleAngle % (2.0 * Math.PI);
|
||||
double R = stroke / 2.0;
|
||||
double cosT = Math.Cos(theta);
|
||||
double sinT = Math.Sin(theta);
|
||||
double L = conRodLength;
|
||||
|
||||
double s = R * (1 - cosT) + L - Math.Sqrt(L * L - R * R * sinT * sinT);
|
||||
double V = V_clear + pistonArea * s;
|
||||
|
||||
double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT);
|
||||
double dVdθ = pistonArea * (R * sinT + (R * R * sinT * cosT) / sqrtTerm);
|
||||
double dvdt = dVdθ * crankshaft.AngularVelocity;
|
||||
return (V, dvdt);
|
||||
}
|
||||
|
||||
// ---- Valve lifts (cycle‑aware) ----
|
||||
private double ExhaustValveLift()
|
||||
{
|
||||
double a = CycleAngle;
|
||||
if (a < exhValveOpenStart || a > exhValveOpenEnd) return 0.0;
|
||||
double duration = exhValveOpenEnd - exhValveOpenStart;
|
||||
double ramp = exhValveRampWidth;
|
||||
double t = (a - exhValveOpenStart) / duration;
|
||||
double rampFrac = ramp / duration;
|
||||
if (t < rampFrac) return t / rampFrac;
|
||||
if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac;
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
private double IntakeValveLift()
|
||||
{
|
||||
double a = CycleAngle;
|
||||
if (a < intValveOpenStart || a > intValveOpenEnd) return 0.0;
|
||||
double duration = intValveOpenEnd - intValveOpenStart;
|
||||
double ramp = intValveRampWidth;
|
||||
double t = (a - intValveOpenStart) / duration;
|
||||
double rampFrac = ramp / duration;
|
||||
if (t < rampFrac) return t / rampFrac;
|
||||
if (t > 1.0 - rampFrac) return (1.0 - t) / rampFrac;
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
// ---- Wiebe burn fraction ----
|
||||
private double WiebeFraction(double angleFromIgnition)
|
||||
{
|
||||
if (angleFromIgnition >= burnDuration) return 1.0;
|
||||
double a = 5.0, m = 2.0;
|
||||
double x = angleFromIgnition / burnDuration;
|
||||
return 1.0 - Math.Exp(-a * Math.Pow(x, m + 1));
|
||||
}
|
||||
|
||||
// ---- Torque from pressure ----
|
||||
private double ComputeTorque()
|
||||
{
|
||||
double p = Cylinder.Pressure;
|
||||
double ambient = 101325.0;
|
||||
double force = (p - ambient) * pistonArea;
|
||||
if (force <= 0) return 0.0;
|
||||
|
||||
double theta = crankshaft.CrankAngle % (2.0 * Math.PI);
|
||||
double R = stroke / 2.0;
|
||||
double L = conRodLength;
|
||||
double sinT = Math.Sin(theta);
|
||||
double cosT = Math.Cos(theta);
|
||||
|
||||
double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT);
|
||||
double lever = R * sinT * (1.0 + (R * cosT) / sqrtTerm);
|
||||
return force * lever;
|
||||
}
|
||||
|
||||
// ---- TDC detection (power stroke, at angle 0 mod 4π) ----
|
||||
private bool DetectTDCPowerStroke()
|
||||
{
|
||||
double current = CycleAngle;
|
||||
double previous = prevCycleAngle;
|
||||
prevCycleAngle = current;
|
||||
return (previous > 3.8 * Math.PI && current < 0.2 * Math.PI);
|
||||
}
|
||||
|
||||
public void Step(double dt)
|
||||
{
|
||||
bool crossingTDC = DetectTDCPowerStroke();
|
||||
|
||||
if (crossingTDC)
|
||||
{
|
||||
misfireCurrent = rand.NextDouble() < MisfireProbability;
|
||||
|
||||
// *** Always capture the state at TDC, whether we burn or not ***
|
||||
preIgnitionMass = Cylinder.Mass;
|
||||
preIgnitionInternalEnergy = Cylinder.InternalEnergy;
|
||||
|
||||
if (misfireCurrent)
|
||||
{
|
||||
MisfireCount++;
|
||||
}
|
||||
else if (ignition)
|
||||
{
|
||||
double V = Cylinder.Volume;
|
||||
targetBurnEnergy = TargetPeakPressure * V / (Cylinder.Gamma - 1.0);
|
||||
if (double.IsNaN(targetBurnEnergy))
|
||||
targetBurnEnergy = 101325.0 * V / (Cylinder.Gamma - 1.0);
|
||||
burnInProgress = true;
|
||||
burnStartAngle = CycleAngle;
|
||||
CombustionCount++;
|
||||
}
|
||||
}
|
||||
|
||||
if (burnInProgress)
|
||||
{
|
||||
double angleFromIgnition = CycleAngle - burnStartAngle;
|
||||
if (angleFromIgnition < 0) angleFromIgnition += 4.0 * Math.PI;
|
||||
|
||||
if (angleFromIgnition >= burnDuration)
|
||||
{
|
||||
Cylinder.InternalEnergy = targetBurnEnergy;
|
||||
burnInProgress = false;
|
||||
}
|
||||
else
|
||||
{
|
||||
double fraction = WiebeFraction(angleFromIgnition);
|
||||
Cylinder.InternalEnergy = preIgnitionInternalEnergy * (1.0 - fraction)
|
||||
+ targetBurnEnergy * fraction;
|
||||
Cylinder.Mass = preIgnitionMass;
|
||||
}
|
||||
}
|
||||
|
||||
var (vol, dvdt) = PistonKinematics(CycleAngle);
|
||||
Cylinder.Volume = vol;
|
||||
Cylinder.Dvdt = dvdt;
|
||||
|
||||
if (double.IsNaN(Cylinder.Pressure) || double.IsNaN(Cylinder.Temperature) || Cylinder.Mass < 1e-9)
|
||||
{
|
||||
double V = Math.Max(vol, V_clear);
|
||||
Cylinder.Mass = 1.225 * V;
|
||||
Cylinder.InternalEnergy = 101325.0 * V / (1.4 - 1.0);
|
||||
}
|
||||
|
||||
double torque = ComputeTorque();
|
||||
crankshaft.AddTorque(torque);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -1,616 +1,318 @@
|
||||
using System;
|
||||
using System.Numerics;
|
||||
using FluidSim.Interfaces;
|
||||
|
||||
namespace FluidSim.Components
|
||||
{
|
||||
public enum BoundaryType
|
||||
{
|
||||
OpenEnd,
|
||||
ZeroPressureOpen,
|
||||
ClosedEnd,
|
||||
GhostCell
|
||||
}
|
||||
|
||||
public class Pipe1D
|
||||
/// <summary>
|
||||
/// 1‑D compressible Euler pipe (finite‑volume, HLLC flux).
|
||||
/// Boundary conditions are set externally via SetGhostLeft/Right.
|
||||
/// Enforces that ghosts are always valid before stepping.
|
||||
/// Uses exponential damping and Newtonian energy relaxation.
|
||||
/// </summary>
|
||||
public class Pipe1D : IComponent
|
||||
{
|
||||
public Port PortA { get; }
|
||||
public Port PortB { get; }
|
||||
public double Area => _area;
|
||||
public double Area { get; }
|
||||
public double DampingMultiplier { get; set; } = 1.0;
|
||||
public double EnergyRelaxationRate { get; set; } = 0.0; // 1/s
|
||||
|
||||
private int _n; // number of real cells
|
||||
private float _dx, _dt; // spatial step, global time step
|
||||
private float _area, _diameter; // cross‑section, diameter
|
||||
private float _gamma; // ratio of specific heats (1.4)
|
||||
|
||||
// Conserved variables – arrays sized [_n] (only real cells, ghosts handled externally)
|
||||
private float[] _rho;
|
||||
private float[] _rhou;
|
||||
private float[] _E;
|
||||
|
||||
// Flux arrays for faces 0 .. _n (face i is between cell i-1 and i)
|
||||
private float[] _fluxM; // mass flux
|
||||
private float[] _fluxP; // momentum flux
|
||||
private float[] _fluxE; // energy flux
|
||||
|
||||
// Ghost cell states
|
||||
private float _rhoGhostL, _uGhostL, _pGhostL;
|
||||
private float _rhoGhostR, _uGhostR, _pGhostR;
|
||||
private bool _ghostLSet, _ghostRSet;
|
||||
|
||||
private BoundaryType _aBCType = BoundaryType.GhostCell;
|
||||
private BoundaryType _bBCType = BoundaryType.GhostCell;
|
||||
|
||||
private float _aAmbientPressure = 101325f;
|
||||
private float _bAmbientPressure = 101325f;
|
||||
|
||||
// CFL / sub-stepping
|
||||
private const float CflTarget = 0.8f;
|
||||
private const float ReferenceSoundSpeed = 340f;
|
||||
private float _lastMassFlow = 0f;
|
||||
|
||||
// Pre‑computed for damping
|
||||
private float _laminarCoeff; // 8*mu / r^2, then multiplied by DampingMultiplier
|
||||
|
||||
// ---- Energy loss (Newton cooling) ----
|
||||
private float _ambientEnergyReference; // total energy density at ambient (Pamb / (γ-1))
|
||||
public float EnergyRelaxationRate { get; set; } = 0.0f; // 1/s
|
||||
|
||||
public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
|
||||
// Ambient pressure for the energy relaxation term (default 101325 Pa)
|
||||
private double _ambientPressure = 101325.0;
|
||||
public double AmbientPressure
|
||||
{
|
||||
float dtGlobal = 1f / sampleRate;
|
||||
int nCells;
|
||||
float dxTarget = ReferenceSoundSpeed * dtGlobal / CflTarget;
|
||||
|
||||
if (forcedCellCount > 1)
|
||||
nCells = forcedCellCount;
|
||||
else
|
||||
get => _ambientPressure;
|
||||
set
|
||||
{
|
||||
nCells = Math.Max(2, (int)Math.Round((float)length / dxTarget, MidpointRounding.AwayFromZero));
|
||||
while (length / nCells > dxTarget * 1.01f && nCells < int.MaxValue - 1)
|
||||
nCells++;
|
||||
_ambientPressure = value;
|
||||
_ambientEnergyReference = value / (_gamma - 1.0);
|
||||
}
|
||||
|
||||
_n = nCells;
|
||||
_dx = (float)(length / nCells);
|
||||
_dt = dtGlobal;
|
||||
_area = (float)area;
|
||||
_diameter = (float)(2.0 * Math.Sqrt(area / Math.PI));
|
||||
_gamma = 1.4f;
|
||||
|
||||
_rho = new float[_n];
|
||||
_rhou = new float[_n];
|
||||
_E = new float[_n];
|
||||
|
||||
// +1 because there are _n+1 faces
|
||||
_fluxM = new float[_n + 1];
|
||||
_fluxP = new float[_n + 1];
|
||||
_fluxE = new float[_n + 1];
|
||||
|
||||
// Pre‑compute laminar damping coefficient (using air at 20°C)
|
||||
float mu_air = 1.8e-5f;
|
||||
float radius = _diameter * 0.5f;
|
||||
_laminarCoeff = 8f * mu_air / (radius * radius); // will be multiplied by DampingMultiplier at each step
|
||||
|
||||
// Ambient reference energy (internal energy per unit volume at 101325 Pa)
|
||||
_ambientEnergyReference = 101325f / (_gamma - 1f); // ≈ 253312.5 J/m³
|
||||
|
||||
PortA = new Port();
|
||||
PortB = new Port();
|
||||
}
|
||||
|
||||
// ==================== PUBLIC API ============================
|
||||
public void SetABoundaryType(BoundaryType type) => _aBCType = type;
|
||||
public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
|
||||
public void SetAAmbientPressure(double p) => _aAmbientPressure = (float)p;
|
||||
public void SetBAmbientPressure(double p) => _bAmbientPressure = (float)p;
|
||||
// Geometry
|
||||
private readonly int _n; // number of real cells
|
||||
private readonly double _dx; // cell size (m)
|
||||
private readonly double _diameter; // m
|
||||
private readonly double _gamma = 1.4;
|
||||
|
||||
public float GetFaceMassFlux(int faceIndex)
|
||||
// Conserved variables [0 .. _n-1]
|
||||
private double[] _rho;
|
||||
private double[] _rhou;
|
||||
private double[] _E;
|
||||
|
||||
// Face fluxes [0 .. _n]
|
||||
private double[] _fluxM;
|
||||
private double[] _fluxP;
|
||||
private double[] _fluxE;
|
||||
|
||||
// Ghost cells (set externally)
|
||||
private double _rhoGhostL, _uGhostL, _pGhostL;
|
||||
private double _rhoGhostR, _uGhostR, _pGhostR;
|
||||
private bool _ghostLValid, _ghostRValid;
|
||||
|
||||
// Pre‑computed damping coefficient
|
||||
private double _laminarCoeff;
|
||||
private double _ambientEnergyReference; // internal energy density at ambient pressure
|
||||
|
||||
/// <summary>
|
||||
/// Initialise a pipe with a given cell count.
|
||||
/// </summary>
|
||||
/// <param name="length">Pipe length (m).</param>
|
||||
/// <param name="area">Cross‑sectional area (m²).</param>
|
||||
/// <param name="cellCount">Number of finite‑volume cells (≥ 4).</param>
|
||||
public Pipe1D(double length, double area, int cellCount)
|
||||
{
|
||||
if (faceIndex < 0 || faceIndex > _n) return 0f;
|
||||
return _fluxM[faceIndex];
|
||||
if (cellCount < 4) throw new ArgumentException("cellCount must be at least 4");
|
||||
|
||||
_n = cellCount;
|
||||
_dx = length / _n;
|
||||
Area = area;
|
||||
_diameter = 2.0 * Math.Sqrt(area / Math.PI);
|
||||
|
||||
_rho = new double[_n];
|
||||
_rhou = new double[_n];
|
||||
_E = new double[_n];
|
||||
|
||||
_fluxM = new double[_n + 1];
|
||||
_fluxP = new double[_n + 1];
|
||||
_fluxE = new double[_n + 1];
|
||||
|
||||
// Laminar damping coefficient for air at 20°C (multiplied by DampingMultiplier each step)
|
||||
double mu_air = 1.8e-5;
|
||||
double radius = _diameter * 0.5;
|
||||
_laminarCoeff = 8.0 * mu_air / (radius * radius);
|
||||
|
||||
// Ambient energy reference (internal energy per unit volume at 101325 Pa)
|
||||
_ambientEnergyReference = 101325.0 / (_gamma - 1.0);
|
||||
|
||||
PortA = new Port { Owner = this };
|
||||
PortB = new Port { Owner = this };
|
||||
|
||||
// Initial state = still air at ambient conditions
|
||||
SetUniformState(1.225, 0.0, 101325.0);
|
||||
}
|
||||
|
||||
IReadOnlyList<Port> IComponent.Ports => new[] { PortA, PortB };
|
||||
|
||||
// No integration needed for pipes – state is advanced via sub‑steps
|
||||
public void UpdateState(double dt) { }
|
||||
|
||||
// ---------- Ghost cell interface ----------
|
||||
public void SetGhostLeft(double rho, double u, double p)
|
||||
{
|
||||
_rhoGhostL = (float)rho;
|
||||
_uGhostL = (float)u;
|
||||
_pGhostL = (float)p;
|
||||
_ghostLSet = true;
|
||||
_rhoGhostL = rho;
|
||||
_uGhostL = u;
|
||||
_pGhostL = p;
|
||||
_ghostLValid = true;
|
||||
}
|
||||
|
||||
public void SetGhostRight(double rho, double u, double p)
|
||||
{
|
||||
_rhoGhostR = (float)rho;
|
||||
_uGhostR = (float)u;
|
||||
_pGhostR = (float)p;
|
||||
_ghostRSet = true;
|
||||
}
|
||||
public void ClearGhostFlag()
|
||||
{
|
||||
_ghostLSet = false;
|
||||
_ghostRSet = false;
|
||||
_rhoGhostR = rho;
|
||||
_uGhostR = u;
|
||||
_pGhostR = p;
|
||||
_ghostRValid = true;
|
||||
}
|
||||
|
||||
public void SetUniformState(double rho, double u, double p)
|
||||
public void ClearGhostFlags()
|
||||
{
|
||||
float r = (float)rho;
|
||||
float vel = (float)u;
|
||||
float pr = (float)p;
|
||||
float e = pr / ((_gamma - 1f) * r);
|
||||
float Etot = r * e + 0.5f * r * vel * vel;
|
||||
for (int i = 0; i < _n; i++)
|
||||
{
|
||||
_rho[i] = r;
|
||||
_rhou[i] = r * vel;
|
||||
_E[i] = Etot;
|
||||
}
|
||||
_ghostLValid = false;
|
||||
_ghostRValid = false;
|
||||
}
|
||||
|
||||
public int GetCellCount() => _n;
|
||||
public (double rho, double u, double p) GetInteriorStateLeft()
|
||||
{
|
||||
double rho = Math.Max(_rho[0], 1e-12);
|
||||
double u = _rhou[0] / rho;
|
||||
double p = PressureScalar(0);
|
||||
return (rho, u, p);
|
||||
}
|
||||
|
||||
public (double rho, double u, double p) GetInteriorStateRight()
|
||||
{
|
||||
double rho = Math.Max(_rho[_n - 1], 1e-12);
|
||||
double u = _rhou[_n - 1] / rho;
|
||||
double p = PressureScalar(_n - 1);
|
||||
return (rho, u, p);
|
||||
}
|
||||
|
||||
public int CellCount => _n;
|
||||
|
||||
public double GetCellDensity(int i) => _rho[i];
|
||||
public double GetCellVelocity(int i)
|
||||
{
|
||||
float rho = Math.Max(_rho[i], 1e-12f);
|
||||
double rho = Math.Max(_rho[i], 1e-12);
|
||||
return _rhou[i] / rho;
|
||||
}
|
||||
public double GetCellPressure(int i)
|
||||
public double GetCellPressure(int i) => PressureScalar(i);
|
||||
|
||||
// ---------- Sub‑stepping ----------
|
||||
public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
|
||||
{
|
||||
float rho = Math.Max(_rho[i], 1e-12f);
|
||||
return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
|
||||
}
|
||||
|
||||
public double GetPressureAtFraction(double fraction)
|
||||
{
|
||||
int i = (int)(fraction * (_n - 1));
|
||||
i = Math.Clamp(i, 0, _n - 1);
|
||||
return GetCellPressure(i);
|
||||
}
|
||||
|
||||
public void SetCellState(int i, double rho, double u, double p)
|
||||
{
|
||||
if (i < 0 || i >= _n) return;
|
||||
float r = (float)rho;
|
||||
float vel = (float)u;
|
||||
float pr = (float)p;
|
||||
_rho[i] = r;
|
||||
_rhou[i] = r * vel;
|
||||
float e = pr / ((_gamma - 1f) * r);
|
||||
_E[i] = r * e + 0.5f * r * vel * vel;
|
||||
}
|
||||
|
||||
public double GetOpenEndMassFlow()
|
||||
{
|
||||
if (_bBCType != BoundaryType.OpenEnd && _bBCType != BoundaryType.ZeroPressureOpen)
|
||||
return 0.0;
|
||||
|
||||
int lastCell = _n - 1;
|
||||
float rho = _rho[lastCell];
|
||||
float u = _rhou[lastCell] / Math.Max(rho, 1e-12f);
|
||||
float p = PressureScalar(lastCell);
|
||||
|
||||
float c = MathF.Sqrt(_gamma * p / rho);
|
||||
float uFace = u;
|
||||
float rhoFace = rho;
|
||||
float pFace = p;
|
||||
|
||||
if (uFace > 0 && uFace < c) // subsonic outflow
|
||||
{
|
||||
float s = p / MathF.Pow(rho, _gamma);
|
||||
float rhoAmb = MathF.Pow(_bAmbientPressure / s, 1f / _gamma);
|
||||
float cAmb = MathF.Sqrt(_gamma * _bAmbientPressure / rhoAmb);
|
||||
float J_plus = u + 2f * c / (_gamma - 1f);
|
||||
float uFaceNew = J_plus - 2f * cAmb / (_gamma - 1f);
|
||||
if (uFaceNew > 0) uFace = uFaceNew;
|
||||
if (uFace < 0) uFace = 0;
|
||||
rhoFace = rhoAmb;
|
||||
pFace = _bAmbientPressure;
|
||||
}
|
||||
|
||||
return rhoFace * uFace * _area;
|
||||
}
|
||||
|
||||
public double GetAndStoreMassFlowForDerivative()
|
||||
{
|
||||
float current = (float)GetOpenEndMassFlow();
|
||||
double derivative = (current - _lastMassFlow) / _dt;
|
||||
_lastMassFlow = current;
|
||||
return derivative;
|
||||
}
|
||||
|
||||
public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8f)
|
||||
{
|
||||
float maxW = 0f;
|
||||
double maxW = 0.0;
|
||||
for (int i = 0; i < _n; i++)
|
||||
{
|
||||
float rho = _rho[i];
|
||||
float u = MathF.Abs(_rhou[i] / Math.Max(rho, 1e-12f));
|
||||
float p = PressureScalar(i);
|
||||
float c = MathF.Sqrt(_gamma * p / Math.Max(rho, 1e-12f));
|
||||
float local = u + c;
|
||||
double rho = Math.Max(_rho[i], 1e-12);
|
||||
double u = Math.Abs(_rhou[i] / rho);
|
||||
double p = PressureScalar(i);
|
||||
double c = Math.Sqrt(_gamma * p / rho);
|
||||
double local = u + c;
|
||||
if (local > maxW) maxW = local;
|
||||
}
|
||||
maxW = Math.Max(maxW, 1e-8f);
|
||||
return Math.Max(1, (int)Math.Ceiling((float)dtGlobal * maxW / ((float)cflTarget * _dx)));
|
||||
maxW = Math.Max(maxW, 1e-8);
|
||||
return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
|
||||
}
|
||||
|
||||
// ==================== MAIN SIMULATION ==================================
|
||||
// ---------- Main simulation step (per sub‑step) ----------
|
||||
public void SimulateSingleStep(double dtSub)
|
||||
{
|
||||
float dt = (float)dtSub;
|
||||
// Enforce that both ends have been provided with ghost states
|
||||
if (!_ghostLValid || !_ghostRValid)
|
||||
throw new InvalidOperationException("Pipe boundary ghosts not set before SimulateSingleStep.");
|
||||
|
||||
double dt = dtSub;
|
||||
int n = _n;
|
||||
|
||||
// --- 1. Left boundary face (index 0) – scalar -----------------------
|
||||
float rhoL = _rho[0];
|
||||
float uL = _rhou[0] / Math.Max(rhoL, 1e-12f);
|
||||
float pL = PressureScalar(0);
|
||||
ComputeLeftBoundaryFlux(rhoL, uL, pL, out _fluxM[0], out _fluxP[0], out _fluxE[0]);
|
||||
// Left boundary face (index 0)
|
||||
HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, _rho[0], _rhou[0] / _rho[0], PressureScalar(0),
|
||||
out _fluxM[0], out _fluxP[0], out _fluxE[0]);
|
||||
|
||||
// --- 2. Internal faces (1 .. n-1) – SIMD ---------------------------
|
||||
int vectorSize = Vector<float>.Count;
|
||||
int lastSimdFace = n - vectorSize; // highest face index that starts a full vector block
|
||||
for (int f = 1; f <= lastSimdFace; f += vectorSize)
|
||||
// Internal faces (1 .. n-1)
|
||||
for (int f = 1; f < n; f++)
|
||||
{
|
||||
SimdInternalFluxBlock(f, vectorSize);
|
||||
}
|
||||
// Scalar remainder for faces f .. n-1
|
||||
for (int f = Math.Max(1, lastSimdFace + 1); f < n; f++)
|
||||
{
|
||||
float rhoLi = _rho[f - 1];
|
||||
float uLi = _rhou[f - 1] / Math.Max(rhoLi, 1e-12f);
|
||||
float pLi = PressureScalar(f - 1);
|
||||
float rhoRi = _rho[f];
|
||||
float uRi = _rhou[f] / Math.Max(rhoRi, 1e-12f);
|
||||
float pRi = PressureScalar(f);
|
||||
HLLCFluxScalar(rhoLi, uLi, pLi, rhoRi, uRi, pRi,
|
||||
out _fluxM[f], out _fluxP[f], out _fluxE[f]);
|
||||
double rhoL = Math.Max(_rho[f - 1], 1e-12);
|
||||
double uL = _rhou[f - 1] / rhoL;
|
||||
double pL = PressureScalar(f - 1);
|
||||
double rhoR = Math.Max(_rho[f], 1e-12);
|
||||
double uR = _rhou[f] / rhoR;
|
||||
double pR = PressureScalar(f);
|
||||
HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out _fluxM[f], out _fluxP[f], out _fluxE[f]);
|
||||
}
|
||||
|
||||
// --- 3. Right boundary face (index n) – scalar --------------------
|
||||
float rhoR = _rho[n - 1];
|
||||
float uR = _rhou[n - 1] / Math.Max(rhoR, 1e-12f);
|
||||
float pR = PressureScalar(n - 1);
|
||||
ComputeRightBoundaryFlux(rhoR, uR, pR, out _fluxM[n], out _fluxP[n], out _fluxE[n]);
|
||||
// Right boundary face (index n)
|
||||
HLLCFlux(_rho[_n - 1], _rhou[_n - 1] / _rho[_n - 1], PressureScalar(_n - 1),
|
||||
_rhoGhostR, _uGhostR, _pGhostR,
|
||||
out _fluxM[n], out _fluxP[n], out _fluxE[n]);
|
||||
|
||||
// --- 4. Cell update + damping + energy loss – SIMD -----------------
|
||||
SimdCellUpdate(dt);
|
||||
}
|
||||
// Cell update
|
||||
double dt_dx = dt / _dx;
|
||||
double coeff = _laminarCoeff * DampingMultiplier;
|
||||
double relaxRate = EnergyRelaxationRate;
|
||||
|
||||
// ==================== PRIVATE SCALAR HELPERS ===========================
|
||||
private float PressureScalar(int i)
|
||||
{
|
||||
float rho = Math.Max(_rho[i], 1e-12f);
|
||||
return (_gamma - 1f) * (_E[i] - 0.5f * _rhou[i] * _rhou[i] / rho);
|
||||
}
|
||||
|
||||
private void ComputeLeftBoundaryFlux(float rhoInt, float uInt, float pInt,
|
||||
out float fm, out float fp, out float fe)
|
||||
{
|
||||
if (_aBCType == BoundaryType.GhostCell && _ghostLSet)
|
||||
HLLCFluxScalar(_rhoGhostL, _uGhostL, _pGhostL, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||||
else if (_aBCType == BoundaryType.OpenEnd)
|
||||
OpenEndFluxLeft(rhoInt, uInt, pInt, _aAmbientPressure, out fm, out fp, out fe);
|
||||
else if (_aBCType == BoundaryType.ZeroPressureOpen)
|
||||
for (int i = 0; i < n; i++)
|
||||
{
|
||||
float rhoFace = rhoInt;
|
||||
float uFace = uInt;
|
||||
float pFace = _aAmbientPressure;
|
||||
HLLCFluxScalar(rhoFace, uFace, pFace, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||||
}
|
||||
else if (_aBCType == BoundaryType.ClosedEnd)
|
||||
ClosedEndFlux(rhoInt, uInt, pInt, false, out fm, out fp, out fe);
|
||||
else
|
||||
{ fm = 0; fp = pInt; fe = 0; }
|
||||
}
|
||||
double r = _rho[i];
|
||||
double ru = _rhou[i];
|
||||
double E = _E[i];
|
||||
|
||||
private void ComputeRightBoundaryFlux(float rhoInt, float uInt, float pInt,
|
||||
out float fm, out float fp, out float fe)
|
||||
{
|
||||
if (_bBCType == BoundaryType.GhostCell && _ghostRSet)
|
||||
HLLCFluxScalar(rhoInt, uInt, pInt, _rhoGhostR, _uGhostR, _pGhostR, out fm, out fp, out fe);
|
||||
else if (_bBCType == BoundaryType.OpenEnd)
|
||||
OpenEndFluxRight(rhoInt, uInt, pInt, _bAmbientPressure, out fm, out fp, out fe);
|
||||
else if (_bBCType == BoundaryType.ZeroPressureOpen)
|
||||
{
|
||||
float rhoFace = rhoInt;
|
||||
float uFace = uInt;
|
||||
float pFace = _bAmbientPressure;
|
||||
HLLCFluxScalar(rhoInt, uInt, pInt, rhoFace, uFace, pFace, out fm, out fp, out fe);
|
||||
}
|
||||
else if (_bBCType == BoundaryType.ClosedEnd)
|
||||
ClosedEndFlux(rhoInt, uInt, pInt, true, out fm, out fp, out fe);
|
||||
else
|
||||
{ fm = 0; fp = pInt; fe = 0; }
|
||||
}
|
||||
double dM = _fluxM[i + 1] - _fluxM[i];
|
||||
double dP = _fluxP[i + 1] - _fluxP[i];
|
||||
double dE_flux = _fluxE[i + 1] - _fluxE[i];
|
||||
|
||||
// ==================== SCALAR HLLC & BOUNDARY FLUX ======================
|
||||
private void HLLCFluxScalar(float rL, float uL, float pL, float rR, float uR, float pR,
|
||||
out float fm, out float fp, out float fe)
|
||||
{
|
||||
float cL = MathF.Sqrt(_gamma * pL / Math.Max(rL, 1e-12f));
|
||||
float cR = MathF.Sqrt(_gamma * pR / Math.Max(rR, 1e-12f));
|
||||
float EL = pL / ((_gamma - 1f) * rL) + 0.5f * uL * uL;
|
||||
float ER = pR / ((_gamma - 1f) * rR) + 0.5f * uR * uR;
|
||||
float SL = Math.Min(uL - cL, uR - cR);
|
||||
float SR = Math.Max(uL + cL, uR + cR);
|
||||
double newR = r - dt_dx * dM;
|
||||
double newRu = ru - dt_dx * dP;
|
||||
double newE = E - dt_dx * dE_flux;
|
||||
|
||||
float denom = rL * (SL - uL) - rR * (SR - uR);
|
||||
float Ss = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;
|
||||
|
||||
float FrL_m = rL * uL, FrL_p = rL * uL * uL + pL, FrL_e = (rL * EL + pL) * uL;
|
||||
float 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)
|
||||
{
|
||||
float rsL = rL * (SL - uL) / (SL - Ss);
|
||||
float ps = pL + rL * (SL - uL) * (Ss - uL);
|
||||
float EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
|
||||
fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
|
||||
}
|
||||
else
|
||||
{
|
||||
float rsR = rR * (SR - uR) / (SR - Ss);
|
||||
float ps = pR + rR * (SR - uR) * (Ss - uR);
|
||||
float EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
|
||||
fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
|
||||
}
|
||||
}
|
||||
|
||||
private void OpenEndFluxLeft(float rhoInt, float uInt, float pInt, float pAmb,
|
||||
out float fm, out float fp, out float fe)
|
||||
{
|
||||
float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f));
|
||||
if (uInt <= -cInt) // supersonic inflow
|
||||
{
|
||||
fm = rhoInt * uInt;
|
||||
fp = rhoInt * uInt * uInt + pInt;
|
||||
fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt;
|
||||
return;
|
||||
}
|
||||
if (uInt <= 0) // subsonic inflow
|
||||
{
|
||||
float T0 = 300f, R = 287f;
|
||||
float ghostRho = pAmb / (R * T0);
|
||||
HLLCFluxScalar(ghostRho, 0f, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||||
return;
|
||||
}
|
||||
// subsonic outflow
|
||||
float s = pInt / MathF.Pow(rhoInt, _gamma);
|
||||
float ghostRho2 = MathF.Pow(pAmb / s, 1f / _gamma);
|
||||
float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho2);
|
||||
float J_minus = uInt - 2f * cInt / (_gamma - 1f);
|
||||
float uGhost = J_minus + 2f * cGhost / (_gamma - 1f);
|
||||
if (uGhost < 0) uGhost = 0;
|
||||
HLLCFluxScalar(ghostRho2, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||||
}
|
||||
|
||||
private void OpenEndFluxRight(float rhoInt, float uInt, float pInt, float pAmb,
|
||||
out float fm, out float fp, out float fe)
|
||||
{
|
||||
float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f));
|
||||
if (uInt >= cInt) // supersonic outflow
|
||||
{
|
||||
fm = rhoInt * uInt;
|
||||
fp = rhoInt * uInt * uInt + pInt;
|
||||
fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt;
|
||||
return;
|
||||
}
|
||||
if (uInt >= 0) // subsonic outflow
|
||||
{
|
||||
float s = pInt / MathF.Pow(rhoInt, _gamma);
|
||||
float ghostRho = MathF.Pow(pAmb / s, 1f / _gamma);
|
||||
float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho);
|
||||
float J_plus = uInt + 2f * cInt / (_gamma - 1f);
|
||||
float uGhost = J_plus - 2f * cGhost / (_gamma - 1f);
|
||||
if (uGhost > 0) uGhost = 0;
|
||||
HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho, uGhost, pAmb, out fm, out fp, out fe);
|
||||
return;
|
||||
}
|
||||
// subsonic inflow
|
||||
float T0 = 300f, R = 287f;
|
||||
float ghostRho2 = pAmb / (R * T0);
|
||||
HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho2, 0f, pAmb, out fm, out fp, out fe);
|
||||
}
|
||||
|
||||
private void ClosedEndFlux(float rhoInt, float uInt, float pInt, bool isRight,
|
||||
out float fm, out float fp, out float fe)
|
||||
{
|
||||
float rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt;
|
||||
if (isRight)
|
||||
HLLCFluxScalar(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
|
||||
else
|
||||
HLLCFluxScalar(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||||
}
|
||||
|
||||
// ==================== SIMD INTERNAL FACE ROUTINE ========================
|
||||
private void SimdInternalFluxBlock(int startFace, int count)
|
||||
{
|
||||
int leftIdx = startFace - 1;
|
||||
int rightIdx = startFace;
|
||||
|
||||
Vector<float> rL = new Vector<float>(_rho, leftIdx);
|
||||
Vector<float> ruL = new Vector<float>(_rhou, leftIdx);
|
||||
Vector<float> EL = new Vector<float>(_E, leftIdx);
|
||||
|
||||
Vector<float> rR = new Vector<float>(_rho, rightIdx);
|
||||
Vector<float> ruR = new Vector<float>(_rhou, rightIdx);
|
||||
Vector<float> ER = new Vector<float>(_E, rightIdx);
|
||||
|
||||
Vector<float> uL = ruL / rL;
|
||||
Vector<float> uR = ruR / rR;
|
||||
|
||||
Vector<float> half = new Vector<float>(0.5f);
|
||||
Vector<float> gammaMinus1 = new Vector<float>(_gamma - 1f);
|
||||
Vector<float> gammaVec = new Vector<float>(_gamma);
|
||||
|
||||
Vector<float> pL = gammaMinus1 * (EL - half * ruL * ruL / rL);
|
||||
Vector<float> pR = gammaMinus1 * (ER - half * ruR * ruR / rR);
|
||||
|
||||
Vector<float> cL = Vector.SquareRoot(gammaVec * pL / rL);
|
||||
Vector<float> cR = Vector.SquareRoot(gammaVec * pR / rR);
|
||||
|
||||
Vector<float> SL = Vector.Min(uL - cL, uR - cR);
|
||||
Vector<float> SR = Vector.Max(uL + cL, uR + cR);
|
||||
|
||||
Vector<float> num = (pR - pL) + rL * uL * (SL - uL) - rR * uR * (SR - uR);
|
||||
Vector<float> den = rL * (SL - uL) - rR * (SR - uR);
|
||||
Vector<float> Ss = num / den;
|
||||
|
||||
Vector<float> eL = EL / rL;
|
||||
Vector<float> eR = ER / rR;
|
||||
|
||||
// Left flux
|
||||
Vector<float> Fm_L = ruL;
|
||||
Vector<float> Fp_L = ruL * uL + pL;
|
||||
Vector<float> Fe_L = (EL + pL) * uL;
|
||||
|
||||
// Right flux
|
||||
Vector<float> Fm_R = ruR;
|
||||
Vector<float> Fp_R = ruR * uR + pR;
|
||||
Vector<float> Fe_R = (ER + pR) * uR;
|
||||
|
||||
// Star‑left fluxes
|
||||
Vector<float> diffL = SL - uL;
|
||||
Vector<float> dL_den = SL - Ss;
|
||||
Vector<float> rsL = rL * diffL / dL_den;
|
||||
Vector<float> psSL = pL + rL * diffL * (Ss - uL);
|
||||
Vector<float> EsL = eL + (Ss - uL) * (Ss + pL / (rL * diffL));
|
||||
Vector<float> Fm_starL = rsL * Ss;
|
||||
Vector<float> Fp_starL = rsL * Ss * Ss + psSL;
|
||||
Vector<float> Fe_starL = (rsL * EsL + psSL) * Ss;
|
||||
|
||||
// Star‑right fluxes
|
||||
Vector<float> diffR = SR - uR;
|
||||
Vector<float> dR_den = SR - Ss;
|
||||
Vector<float> rsR = rR * diffR / dR_den;
|
||||
Vector<float> psSR = pR + rR * diffR * (Ss - uR);
|
||||
Vector<float> EsR = eR + (Ss - uR) * (Ss + pR / (rR * diffR));
|
||||
Vector<float> Fm_starR = rsR * Ss;
|
||||
Vector<float> Fp_starR = rsR * Ss * Ss + psSR;
|
||||
Vector<float> Fe_starR = (rsR * EsR + psSR) * Ss;
|
||||
|
||||
var maskSLge0 = Vector.GreaterThanOrEqual(SL, Vector<float>.Zero);
|
||||
var maskSRle0 = Vector.LessThanOrEqual(SR, Vector<float>.Zero);
|
||||
var maskMiddle = ~(maskSLge0 | maskSRle0);
|
||||
var maskStarL = maskMiddle & Vector.GreaterThanOrEqual(Ss, Vector<float>.Zero);
|
||||
var maskStarR = maskMiddle & Vector.LessThan(Ss, Vector<float>.Zero);
|
||||
|
||||
Vector<float> fm = Vector.ConditionalSelect(maskSLge0, Fm_L,
|
||||
Vector.ConditionalSelect(maskSRle0, Fm_R,
|
||||
Vector.ConditionalSelect(maskStarL, Fm_starL,
|
||||
Vector.ConditionalSelect(maskStarR, Fm_starR, Vector<float>.Zero))));
|
||||
|
||||
Vector<float> fp = Vector.ConditionalSelect(maskSLge0, Fp_L,
|
||||
Vector.ConditionalSelect(maskSRle0, Fp_R,
|
||||
Vector.ConditionalSelect(maskStarL, Fp_starL,
|
||||
Vector.ConditionalSelect(maskStarR, Fp_starR, Vector<float>.Zero))));
|
||||
|
||||
Vector<float> fe = Vector.ConditionalSelect(maskSLge0, Fe_L,
|
||||
Vector.ConditionalSelect(maskSRle0, Fe_R,
|
||||
Vector.ConditionalSelect(maskStarL, Fe_starL,
|
||||
Vector.ConditionalSelect(maskStarR, Fe_starR, Vector<float>.Zero))));
|
||||
|
||||
fm.CopyTo(_fluxM, startFace);
|
||||
fp.CopyTo(_fluxP, startFace);
|
||||
fe.CopyTo(_fluxE, startFace);
|
||||
}
|
||||
|
||||
// ==================== SIMD CELL UPDATE + DAMPING + ENERGY LOSS =========
|
||||
private void SimdCellUpdate(float dt)
|
||||
{
|
||||
float dt_dx = dt / _dx;
|
||||
Vector<float> vDtDx = new Vector<float>(dt_dx);
|
||||
float coeff = _laminarCoeff * (float)DampingMultiplier;
|
||||
Vector<float> vCoeff = new Vector<float>(coeff);
|
||||
Vector<float> vDt = new Vector<float>(dt);
|
||||
|
||||
int vectorSize = Vector<float>.Count;
|
||||
int n = _n;
|
||||
int lastSimdCell = n - vectorSize;
|
||||
|
||||
// Pre‑defined constants used in clamping
|
||||
Vector<float> half = new Vector<float>(0.5f);
|
||||
Vector<float> gammaMinus1 = new Vector<float>(_gamma - 1f);
|
||||
Vector<float> ambientEnergyVec = new Vector<float>(_ambientEnergyReference);
|
||||
Vector<float> energyRelaxRateVec = new Vector<float>(EnergyRelaxationRate);
|
||||
|
||||
for (int i = 0; i <= lastSimdCell; i += vectorSize)
|
||||
{
|
||||
// Load conserved
|
||||
Vector<float> r = new Vector<float>(_rho, i);
|
||||
Vector<float> ru = new Vector<float>(_rhou, i);
|
||||
Vector<float> E = new Vector<float>(_E, i);
|
||||
|
||||
// Load fluxes at faces i (left) and i+1 (right)
|
||||
Vector<float> flxM_L = new Vector<float>(_fluxM, i);
|
||||
Vector<float> flxM_R = new Vector<float>(_fluxM, i + 1);
|
||||
Vector<float> flxP_L = new Vector<float>(_fluxP, i);
|
||||
Vector<float> flxP_R = new Vector<float>(_fluxP, i + 1);
|
||||
Vector<float> flxE_L = new Vector<float>(_fluxE, i);
|
||||
Vector<float> flxE_R = new Vector<float>(_fluxE, i + 1);
|
||||
|
||||
// Update conserved: Q_new = Q - dt/dx * (flux_right - flux_left)
|
||||
Vector<float> newR = r - vDtDx * (flxM_R - flxM_L);
|
||||
Vector<float> newRu = ru - vDtDx * (flxP_R - flxP_L);
|
||||
Vector<float> newE = E - vDtDx * (flxE_R - flxE_L);
|
||||
|
||||
// Damping
|
||||
Vector<float> dampingFactor = Vector.Exp(-vCoeff / r * vDt);
|
||||
// Wall friction damping (laminar)
|
||||
double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
|
||||
newRu *= dampingFactor;
|
||||
|
||||
// Energy loss (Newton cooling toward ambient)
|
||||
Vector<float> relaxFactor = Vector.Exp(-energyRelaxRateVec * vDt);
|
||||
newE = ambientEnergyVec + (newE - ambientEnergyVec) * relaxFactor;
|
||||
|
||||
// Clamp density
|
||||
newR = Vector.Max(newR, new Vector<float>(1e-12f));
|
||||
|
||||
// Enforce minimal pressure (100 Pa) -> minimal energy
|
||||
Vector<float> kinE = half * newRu * newRu / newR;
|
||||
Vector<float> eMin = new Vector<float>(100f) / gammaMinus1 + kinE;
|
||||
newE = Vector.Max(newE, eMin);
|
||||
|
||||
newR.CopyTo(_rho, i);
|
||||
newRu.CopyTo(_rhou, i);
|
||||
newE.CopyTo(_E, i);
|
||||
}
|
||||
|
||||
// Scalar remainder
|
||||
float relaxRate = EnergyRelaxationRate;
|
||||
for (int i = Math.Max(0, lastSimdCell + 1); i < n; i++)
|
||||
{
|
||||
float r = _rho[i];
|
||||
float ru = _rhou[i];
|
||||
float E = _E[i];
|
||||
|
||||
float dM = _fluxM[i + 1] - _fluxM[i];
|
||||
float dP = _fluxP[i + 1] - _fluxP[i];
|
||||
float dE_flux = _fluxE[i + 1] - _fluxE[i];
|
||||
|
||||
float newR = r - dt_dx * dM;
|
||||
float newRu = ru - dt_dx * dP;
|
||||
float newE = E - dt_dx * dE_flux;
|
||||
|
||||
// Damping
|
||||
float dampingFactor = MathF.Exp(-coeff / Math.Max(r, 1e-12f) * dt);
|
||||
newRu *= dampingFactor;
|
||||
|
||||
// Energy loss
|
||||
float relaxFactor = MathF.Exp(-relaxRate * dt);
|
||||
// Newtonian cooling toward ambient energy
|
||||
double relaxFactor = Math.Exp(-relaxRate * dt);
|
||||
newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;
|
||||
|
||||
// Clamps
|
||||
if (newR < 1e-12f) newR = 1e-12f;
|
||||
float kin = 0.5f * newRu * newRu / newR;
|
||||
float emin = 100f / (_gamma - 1f) + kin;
|
||||
if (newE < emin) newE = emin;
|
||||
// Clamps – minimum density 1e-12, minimum pressure 100 Pa
|
||||
newR = Math.Max(newR, 1e-12);
|
||||
double kin = 0.5 * newRu * newRu / Math.Max(newR, 1e-12);
|
||||
double eMin = 100.0 / (_gamma - 1.0) + kin;
|
||||
newE = Math.Max(newE, eMin);
|
||||
|
||||
_rho[i] = newR;
|
||||
_rhou[i] = newRu;
|
||||
_E[i] = newE;
|
||||
}
|
||||
|
||||
// Update port states to reflect the current interior state (for audio / monitoring)
|
||||
(double rhoA, double uA, double pA) = GetInteriorStateLeft();
|
||||
PortA.Pressure = pA;
|
||||
PortA.Density = rhoA;
|
||||
PortA.Temperature = pA / (rhoA * 287.0);
|
||||
PortA.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pA / rhoA;
|
||||
|
||||
(double rhoB, double uB, double pB) = GetInteriorStateRight();
|
||||
PortB.Pressure = pB;
|
||||
PortB.Density = rhoB;
|
||||
PortB.Temperature = pB / (rhoB * 287.0);
|
||||
PortB.SpecificEnthalpy = _gamma / (_gamma - 1.0) * pB / rhoB;
|
||||
}
|
||||
|
||||
// ---------- Private helpers ----------
|
||||
private double PressureScalar(int i)
|
||||
{
|
||||
double rho = Math.Max(_rho[i], 1e-12);
|
||||
return (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / rho);
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// HLLC approximate Riemann solver (Toro, 1997).
|
||||
/// Computes the numerical flux at a face given left and right states.
|
||||
/// </summary>
|
||||
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 / rL);
|
||||
double cR = Math.Sqrt(_gamma * pR / rR);
|
||||
double EL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL; // specific total energy
|
||||
double ER = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
|
||||
|
||||
// Wave speed estimates (Davis, 1988)
|
||||
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 = (pR - pL + rL * uL * (SL - uL) - rR * uR * (SR - uR)) / denom;
|
||||
|
||||
double Fm_L = rL * uL;
|
||||
double Fp_L = rL * uL * uL + pL;
|
||||
double Fe_L = (rL * EL + pL) * uL;
|
||||
|
||||
double Fm_R = rR * uR;
|
||||
double Fp_R = rR * uR * uR + pR;
|
||||
double Fe_R = (rR * ER + pR) * uR;
|
||||
|
||||
if (SL >= 0) { fm = Fm_L; fp = Fp_L; fe = Fe_L; }
|
||||
else if (SR <= 0) { fm = Fm_R; fp = Fp_R; fe = Fe_R; }
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
/// <summary>Initialise all cells to a uniform state (rho, u, p).</summary>
|
||||
public void SetUniformState(double rho, double u, double p)
|
||||
{
|
||||
double e = p / ((_gamma - 1.0) * rho);
|
||||
double E = rho * e + 0.5 * rho * u * u;
|
||||
for (int i = 0; i < _n; i++)
|
||||
{
|
||||
_rho[i] = rho;
|
||||
_rhou[i] = rho * u;
|
||||
_E[i] = E;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -1,69 +1,109 @@
|
||||
using System;
|
||||
using System.Collections.Generic;
|
||||
using FluidSim.Interfaces;
|
||||
|
||||
namespace FluidSim.Components
|
||||
{
|
||||
public class Volume0D
|
||||
/// <summary>
|
||||
/// Zero‑dimensional control volume with arbitrary number of ports.
|
||||
/// Integrates mass and energy fluxes from all ports.
|
||||
/// Safeguards keep a tiny amount of gas to avoid negative states.
|
||||
/// </summary>
|
||||
public class Volume0D : IComponent
|
||||
{
|
||||
public double Mass { get; set; }
|
||||
public double InternalEnergy { get; set; }
|
||||
public List<Port> Ports { get; } = new List<Port>();
|
||||
|
||||
public double Mass { get; private set; }
|
||||
public double InternalEnergy { get; private set; }
|
||||
public double Volume { get; set; }
|
||||
public double Dvdt { get; set; }
|
||||
public double Gamma { get; set; } = 1.4;
|
||||
public double GasConstant { get; set; } = 287.0;
|
||||
|
||||
public double Volume { get; set; }
|
||||
public double Dvdt { get; set; }
|
||||
|
||||
private double _dt;
|
||||
// Ambient pressure used for emergency refill – default 101325 Pa
|
||||
public double AmbientPressure { get; set; } = 101325.0;
|
||||
|
||||
// Derived quantities
|
||||
public double Density => Mass / Math.Max(Volume, 1e-12);
|
||||
public double Pressure => (Gamma - 1.0) * InternalEnergy / Math.Max(Volume, 1e-12);
|
||||
public double Temperature => Pressure / Math.Max(Density * GasConstant, 1e-12);
|
||||
public double SpecificEnthalpy => Gamma / (Gamma - 1.0) * Pressure / Math.Max(Density, 1e-12);
|
||||
|
||||
public double MassFlowRateIn { get; set; }
|
||||
public double SpecificEnthalpyIn { get; set; }
|
||||
|
||||
public Volume0D(double initialVolume, double initialPressure,
|
||||
double initialTemperature, int sampleRate,
|
||||
double gasConstant = 287.0, double gamma = 1.4)
|
||||
double initialTemperature, double gasConstant = 287.0, double gamma = 1.4)
|
||||
{
|
||||
GasConstant = gasConstant;
|
||||
Gamma = gamma;
|
||||
Volume = initialVolume;
|
||||
Dvdt = 0.0;
|
||||
_dt = 1.0 / sampleRate;
|
||||
|
||||
double rho0 = initialPressure / (GasConstant * initialTemperature);
|
||||
Mass = rho0 * Volume;
|
||||
InternalEnergy = (initialPressure * Volume) / (Gamma - 1.0);
|
||||
}
|
||||
|
||||
public void Integrate(double dtOverride)
|
||||
/// <summary>Add a new port and initialise it to the volume's current state.</summary>
|
||||
public Port CreatePort()
|
||||
{
|
||||
double dm = MassFlowRateIn * dtOverride;
|
||||
double dE = (MassFlowRateIn * SpecificEnthalpyIn) * dtOverride - Pressure * Dvdt * dtOverride;
|
||||
var port = new Port { Owner = this };
|
||||
// Set the port state immediately to avoid a mismatch before the first integration
|
||||
port.Pressure = Pressure;
|
||||
port.Density = Density;
|
||||
port.Temperature = Temperature;
|
||||
port.SpecificEnthalpy = SpecificEnthalpy;
|
||||
Ports.Add(port);
|
||||
return port;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Integrate over dt using the MassFlowRate and SpecificEnthalpy
|
||||
/// that have been set on each port during the coupling resolution phase.
|
||||
/// </summary>
|
||||
public void UpdateState(double dt)
|
||||
{
|
||||
double totalMdot = 0.0;
|
||||
double totalEdot = 0.0;
|
||||
|
||||
foreach (var port in Ports)
|
||||
{
|
||||
totalMdot += port.MassFlowRate;
|
||||
// mdot * h gives energy flow: positive mdot = inflow, negative = outflow
|
||||
totalEdot += port.MassFlowRate * port.SpecificEnthalpy;
|
||||
}
|
||||
|
||||
double dm = totalMdot * dt;
|
||||
double dE = totalEdot * dt - Pressure * Dvdt * dt; // piston work
|
||||
|
||||
Mass += dm;
|
||||
InternalEnergy += dE;
|
||||
|
||||
// ---- ABSOLUTE SAFEGUARD: keep at least 1 µg of gas at ambient pressure ----
|
||||
double minMass = 1e-9;
|
||||
// Safeguards: keep at least 1 µg of gas at a small pressure
|
||||
double V = Math.Max(Volume, 1e-12);
|
||||
if (Mass < minMass)
|
||||
if (Mass < 1e-9)
|
||||
{
|
||||
Mass = minMass;
|
||||
InternalEnergy = 5000.0 * V / (Gamma - 1.0); // 0.05 bar, not ambient
|
||||
Mass = 1e-9;
|
||||
InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
|
||||
}
|
||||
else if (InternalEnergy < 0.0)
|
||||
{
|
||||
InternalEnergy = 101325.0 * V / (Gamma - 1.0);
|
||||
InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
|
||||
}
|
||||
|
||||
// Final non‑negative clamp
|
||||
if (Mass < 0.0) Mass = 0.0;
|
||||
if (InternalEnergy < 0.0) InternalEnergy = 0.0;
|
||||
// Final non‑negative clamps (should not be needed after above)
|
||||
if (Mass < 0.0) Mass = 1e-9;
|
||||
if (InternalEnergy < 0.0) InternalEnergy = AmbientPressure * V / (Gamma - 1.0);
|
||||
|
||||
// Push updated state back to all ports
|
||||
double p = Pressure, rho = Density, T = Temperature, h = SpecificEnthalpy;
|
||||
foreach (var port in Ports)
|
||||
{
|
||||
port.Pressure = p;
|
||||
port.Density = rho;
|
||||
port.Temperature = T;
|
||||
port.SpecificEnthalpy = h; // will be overwritten by couplings for inflow, but this is the default
|
||||
}
|
||||
}
|
||||
|
||||
public void Integrate() => Integrate(_dt);
|
||||
IReadOnlyList<Port> IComponent.Ports => Ports;
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user