grillkol/FluidSim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

insane engine sound

Browse Source
This commit is contained in:
grillkol
2026-05-05 11:24:32 +02:00
parent d963032e74
commit f16a1aa763
4 changed files with 666 additions and 316 deletions
Download Patch File Download Diff File Expand all files Collapse all files

623
Components/Pipe1D.cs
View File

@@ -1,4 +1,5 @@
using System; using System;
using System.Numerics;
using FluidSim.Interfaces; using FluidSim.Interfaces;
namespace FluidSim.Components namespace FluidSim.Components
@@ -6,7 +7,7 @@ namespace FluidSim.Components
public enum BoundaryType public enum BoundaryType
{ {
OpenEnd, OpenEnd,
ZeroPressureOpen, // pressure‑release boundary (strong reflection) ZeroPressureOpen,
ClosedEnd, ClosedEnd,
GhostCell GhostCell
} }
@@ -18,72 +19,99 @@ namespace FluidSim.Components
public double Area => _area; public double Area => _area;
public double DampingMultiplier { get; set; } = 1.0; public double DampingMultiplier { get; set; } = 1.0;
private int _n; private int _n; // number of real cells
private double _dx, _dt, _gamma, _area, _diameter; private float _dx, _dt; // spatial step, global time step
private double[] _rho, _rhou, _E; 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 // Ghost cell states
private double _rhoGhostL, _uGhostL, _pGhostL; private float _rhoGhostL, _uGhostL, _pGhostL;
private double _rhoGhostR, _uGhostR, _pGhostR; private float _rhoGhostR, _uGhostR, _pGhostR;
private bool _ghostLSet, _ghostRSet; private bool _ghostLSet, _ghostRSet;
private BoundaryType _aBCType = BoundaryType.GhostCell; private BoundaryType _aBCType = BoundaryType.GhostCell;
private BoundaryType _bBCType = BoundaryType.GhostCell; private BoundaryType _bBCType = BoundaryType.GhostCell;
private double _aAmbientPressure = 101325.0; private float _aAmbientPressure = 101325f;
private double _bAmbientPressure = 101325.0; private float _bAmbientPressure = 101325f;
private const double CflTarget = 0.8; // CFL / sub-stepping
private const double ReferenceSoundSpeed = 340.0; private const float CflTarget = 0.8f;
private double _lastMassFlow = 0.0; private const float ReferenceSoundSpeed = 340f;
private float _lastMassFlow = 0f;
// Pre‑computed for damping
private float _laminarCoeff; // 8*mu / r^2, then multiplied by DampingMultiplier
public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0) public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
{ {
double dtGlobal = 1.0 / sampleRate; float dtGlobal = 1f / sampleRate;
int nCells; int nCells;
float dxTarget = ReferenceSoundSpeed * dtGlobal / CflTarget;
if (forcedCellCount > 1) if (forcedCellCount > 1)
nCells = forcedCellCount; nCells = forcedCellCount;
else else
{ {
double dxTarget = ReferenceSoundSpeed * dtGlobal * CflTarget; nCells = Math.Max(2, (int)Math.Round((float)length / dxTarget, MidpointRounding.AwayFromZero));
nCells = Math.Max(2, (int)Math.Round(length / dxTarget, MidpointRounding.AwayFromZero)); while (length / nCells > dxTarget * 1.01f && nCells < int.MaxValue - 1)
while (length / nCells > dxTarget * 1.01 && nCells < int.MaxValue - 1)
nCells++; nCells++;
} }
_n = nCells; _n = nCells;
_dx = length / _n; _dx = (float)(length / nCells);
_dt = dtGlobal; _dt = dtGlobal;
_area = area; _area = (float)area;
_gamma = 1.4; _diameter = (float)(2.0 * Math.Sqrt(area / Math.PI));
_diameter = 2.0 * Math.Sqrt(area / Math.PI); _gamma = 1.4f;
_rho = new double[_n]; _rho = new float[_n];
_rhou = new double[_n]; _rhou = new float[_n];
_E = new double[_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
PortA = new Port(); PortA = new Port();
PortB = new Port(); PortB = new Port();
} }
// ==================== PUBLIC API (unchanged) ============================
public void SetABoundaryType(BoundaryType type) => _aBCType = type; public void SetABoundaryType(BoundaryType type) => _aBCType = type;
public void SetBBoundaryType(BoundaryType type) => _bBCType = type; public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
public void SetAAmbientPressure(double p) => _aAmbientPressure = p; public void SetAAmbientPressure(double p) => _aAmbientPressure = (float)p;
public void SetBAmbientPressure(double p) => _bAmbientPressure = p; public void SetBAmbientPressure(double p) => _bAmbientPressure = (float)p;
public void SetGhostLeft(double rho, double u, double p) public void SetGhostLeft(double rho, double u, double p)
{ {
_rhoGhostL = rho; _rhoGhostL = (float)rho;
_uGhostL = u; _uGhostL = (float)u;
_pGhostL = p; _pGhostL = (float)p;
_ghostLSet = true; _ghostLSet = true;
} }
public void SetGhostRight(double rho, double u, double p) public void SetGhostRight(double rho, double u, double p)
{ {
_rhoGhostR = rho; _rhoGhostR = (float)rho;
_uGhostR = u; _uGhostR = (float)u;
_pGhostR = p; _pGhostR = (float)p;
_ghostRSet = true; _ghostRSet = true;
} }
public void ClearGhostFlag() public void ClearGhostFlag()
@@ -94,39 +122,73 @@ namespace FluidSim.Components
public void SetUniformState(double rho, double u, double p) public void SetUniformState(double rho, double u, double p)
{ {
double e = p / ((_gamma - 1) * rho); float r = (float)rho;
double Etot = rho * e + 0.5 * rho * u * u; 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++) for (int i = 0; i < _n; i++)
{ {
_rho[i] = rho; _rho[i] = r;
_rhou[i] = rho * u; _rhou[i] = r * vel;
_E[i] = Etot; _E[i] = Etot;
} }
} }
public int GetCellCount() => _n;
public double GetCellDensity(int i) => _rho[i];
public double GetCellVelocity(int i)
{
float rho = Math.Max(_rho[i], 1e-12f);
return _rhou[i] / rho;
}
public double GetCellPressure(int i)
{
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() public double GetOpenEndMassFlow()
{ {
if (_bBCType != BoundaryType.OpenEnd && _bBCType != BoundaryType.ZeroPressureOpen) if (_bBCType != BoundaryType.OpenEnd && _bBCType != BoundaryType.ZeroPressureOpen)
return 0.0; return 0.0;
int lastCell = _n - 1; int lastCell = _n - 1;
double rho = _rho[lastCell]; float rho = _rho[lastCell];
double u = _rhou[lastCell] / Math.Max(rho, 1e-12); float u = _rhou[lastCell] / Math.Max(rho, 1e-12f);
double p = Pressure(lastCell); float p = PressureScalar(lastCell);
double c = Math.Sqrt(_gamma * p / rho); float c = MathF.Sqrt(_gamma * p / rho);
double uFace = u; float uFace = u;
double rhoFace = rho; float rhoFace = rho;
double pFace = p; float pFace = p;
// Subsonic outflow: impose ambient pressure, adjust velocity and density if (uFace > 0 && uFace < c) // subsonic outflow
if (uFace > 0 && uFace < c)
{ {
double s = p / Math.Pow(rho, _gamma); float s = p / MathF.Pow(rho, _gamma);
double rhoAmb = Math.Pow(_bAmbientPressure / s, 1.0 / _gamma); float rhoAmb = MathF.Pow(_bAmbientPressure / s, 1f / _gamma);
double cAmb = Math.Sqrt(_gamma * _bAmbientPressure / rhoAmb); float cAmb = MathF.Sqrt(_gamma * _bAmbientPressure / rhoAmb);
double J_plus = u + 2.0 * c / (_gamma - 1.0); float J_plus = u + 2f * c / (_gamma - 1f);
double uFaceNew = J_plus - 2.0 * cAmb / (_gamma - 1.0); float uFaceNew = J_plus - 2f * cAmb / (_gamma - 1f);
if (uFaceNew > 0) uFace = uFaceNew; if (uFaceNew > 0) uFace = uFaceNew;
if (uFace < 0) uFace = 0; if (uFace < 0) uFace = 0;
rhoFace = rhoAmb; rhoFace = rhoAmb;
@@ -138,229 +200,410 @@ namespace FluidSim.Components
public double GetAndStoreMassFlowForDerivative() public double GetAndStoreMassFlowForDerivative()
{ {
double current = GetOpenEndMassFlow(); float current = (float)GetOpenEndMassFlow();
double derivative = (current - _lastMassFlow) / _dt; double derivative = (current - _lastMassFlow) / _dt;
_lastMassFlow = current; _lastMassFlow = current;
return derivative; return derivative;
} }
public void SetCellState(int i, double rho, double u, double p) public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8f)
{ {
if (i < 0 || i >= _n) return; float maxW = 0f;
_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++) for (int i = 0; i < _n; i++)
{ {
double rho = _rho[i]; float rho = _rho[i];
double u = Math.Abs(_rhou[i] / Math.Max(rho, 1e-12)); float u = MathF.Abs(_rhou[i] / Math.Max(rho, 1e-12f));
double c = Math.Sqrt(_gamma * Pressure(i) / Math.Max(rho, 1e-12)); float p = PressureScalar(i);
double local = u + c; float c = MathF.Sqrt(_gamma * p / Math.Max(rho, 1e-12f));
float local = u + c;
if (local > maxW) maxW = local; if (local > maxW) maxW = local;
} }
maxW = Math.Max(maxW, 1e-8); maxW = Math.Max(maxW, 1e-8f);
return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx))); return Math.Max(1, (int)Math.Ceiling((float)dtGlobal * maxW / ((float)cflTarget * _dx)));
} }
// ==================== MAIN SIMULATION ==================================
public void SimulateSingleStep(double dtSub) public void SimulateSingleStep(double dtSub)
{ {
float dt = (float)dtSub;
int n = _n; 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) // --- 1. Left boundary face (index 0) – scalar -----------------------
double rhoL = _rho[0]; float rhoL = _rho[0];
double uL = _rhou[0] / Math.Max(rhoL, 1e-12); float uL = _rhou[0] / Math.Max(rhoL, 1e-12f);
double pL = Pressure(0); float pL = PressureScalar(0);
ComputeLeftBoundaryFlux(rhoL, uL, pL, 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
// Face index f is between cell f-1 (left) and cell f (right).
// We want to cover faces 1..n-1.
for (int f = 1; f <= lastSimdFace; f += vectorSize)
{
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]);
}
// --- 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]);
// --- 4. Cell update + damping – SIMD ------------------------------
SimdCellUpdate(dt);
}
// ==================== 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) if (_aBCType == BoundaryType.GhostCell && _ghostLSet)
HLLCFlux(_rhoGhostL, _uGhostL, _pGhostL, rhoL, uL, pL, out Fm[0], out Fp[0], out Fe[0]); HLLCFluxScalar(_rhoGhostL, _uGhostL, _pGhostL, rhoInt, uInt, pInt, out fm, out fp, out fe);
else if (_aBCType == BoundaryType.OpenEnd) else if (_aBCType == BoundaryType.OpenEnd)
OpenEndFluxLeft(rhoL, uL, pL, _aAmbientPressure, out Fm[0], out Fp[0], out Fe[0]); OpenEndFluxLeft(rhoInt, uInt, pInt, _aAmbientPressure, out fm, out fp, out fe);
else if (_aBCType == BoundaryType.ZeroPressureOpen) else if (_aBCType == BoundaryType.ZeroPressureOpen)
{ {
// Strong reflection: force pressure to ambient, extrapolate density and velocity float rhoFace = rhoInt;
double rhoFace = rhoL; float uFace = uInt;
double uFace = uL; float pFace = _aAmbientPressure;
double pFace = _aAmbientPressure; HLLCFluxScalar(rhoFace, uFace, pFace, rhoInt, uInt, pInt, out fm, out fp, out fe);
HLLCFlux(rhoFace, uFace, pFace, rhoL, uL, pL, out Fm[0], out Fp[0], out Fe[0]);
} }
else if (_aBCType == BoundaryType.ClosedEnd) else if (_aBCType == BoundaryType.ClosedEnd)
ClosedEndFlux(rhoL, uL, pL, false, out Fm[0], out Fp[0], out Fe[0]); ClosedEndFlux(rhoInt, uInt, pInt, false, out fm, out fp, out fe);
else else
{ Fm[0] = 0; Fp[0] = pL; Fe[0] = 0; } { fm = 0; fp = pInt; fe = 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) private void ComputeRightBoundaryFlux(float rhoInt, float uInt, float pInt,
double rhoR = _rho[n - 1]; out float fm, out float fp, out float fe)
double uR = _rhou[n - 1] / Math.Max(rhoR, 1e-12); {
double pR = Pressure(n - 1);
if (_bBCType == BoundaryType.GhostCell && _ghostRSet) if (_bBCType == BoundaryType.GhostCell && _ghostRSet)
HLLCFlux(rhoR, uR, pR, _rhoGhostR, _uGhostR, _pGhostR, out Fm[n], out Fp[n], out Fe[n]); HLLCFluxScalar(rhoInt, uInt, pInt, _rhoGhostR, _uGhostR, _pGhostR, out fm, out fp, out fe);
else if (_bBCType == BoundaryType.OpenEnd) else if (_bBCType == BoundaryType.OpenEnd)
OpenEndFluxRight(rhoR, uR, pR, _bAmbientPressure, out Fm[n], out Fp[n], out Fe[n]); OpenEndFluxRight(rhoInt, uInt, pInt, _bAmbientPressure, out fm, out fp, out fe);
else if (_bBCType == BoundaryType.ZeroPressureOpen) else if (_bBCType == BoundaryType.ZeroPressureOpen)
{ {
double rhoFace = rhoR; float rhoFace = rhoInt;
double uFace = uR; float uFace = uInt;
double pFace = _bAmbientPressure; float pFace = _bAmbientPressure;
HLLCFlux(rhoR, uR, pR, rhoFace, uFace, pFace, out Fm[n], out Fp[n], out Fe[n]); HLLCFluxScalar(rhoInt, uInt, pInt, rhoFace, uFace, pFace, out fm, out fp, out fe);
} }
else if (_bBCType == BoundaryType.ClosedEnd) else if (_bBCType == BoundaryType.ClosedEnd)
ClosedEndFlux(rhoR, uR, pR, true, out Fm[n], out Fp[n], out Fe[n]); ClosedEndFlux(rhoInt, uInt, pInt, true, out fm, out fp, out fe);
else else
{ Fm[n] = 0; Fp[n] = pR; Fe[n] = 0; } { fm = 0; fp = pInt; fe = 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 ---------- // ==================== SCALAR HLLC & BOUNDARY FLUX ======================
private void HLLCFlux(double rL, double uL, double pL, double rR, double uR, double pR, private void HLLCFluxScalar(float rL, float uL, float pL, float rR, float uR, float pR,
out double fm, out double fp, out double fe) out float fm, out float fp, out float fe)
{ {
double cL = Math.Sqrt(_gamma * pL / Math.Max(rL, 1e-12)); float cL = MathF.Sqrt(_gamma * pL / Math.Max(rL, 1e-12f));
double cR = Math.Sqrt(_gamma * pR / Math.Max(rR, 1e-12)); float cR = MathF.Sqrt(_gamma * pR / Math.Max(rR, 1e-12f));
double EL = pL / ((_gamma - 1) * rL) + 0.5 * uL * uL; float EL = pL / ((_gamma - 1f) * rL) + 0.5f * uL * uL;
double ER = pR / ((_gamma - 1) * rR) + 0.5 * uR * uR; float ER = pR / ((_gamma - 1f) * rR) + 0.5f * uR * uR;
double SL = Math.Min(uL - cL, uR - cR); float SL = Math.Min(uL - cL, uR - cR);
double SR = Math.Max(uL + cL, uR + cR); float 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; float denom = rL * (SL - uL) - rR * (SR - uR);
double FrR_m = rR * uR, FrR_p = rR * uR * uR + pR, FrR_e = (rR * ER + pR) * 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; } 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 (SR <= 0) { fm = FrR_m; fp = FrR_p; fe = FrR_e; }
else if (Ss >= 0) else if (Ss >= 0)
{ {
double rsL = rL * (SL - uL) / (SL - Ss); float rsL = rL * (SL - uL) / (SL - Ss);
double ps = pL + rL * (SL - uL) * (Ss - uL); float ps = pL + rL * (SL - uL) * (Ss - uL);
double EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL))); float EsL = EL + (Ss - uL) * (Ss + pL / (rL * (SL - uL)));
fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss; fm = rsL * Ss; fp = rsL * Ss * Ss + ps; fe = (rsL * EsL + ps) * Ss;
} }
else else
{ {
double rsR = rR * (SR - uR) / (SR - Ss); float rsR = rR * (SR - uR) / (SR - Ss);
double ps = pR + rR * (SR - uR) * (Ss - uR); float ps = pR + rR * (SR - uR) * (Ss - uR);
double EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR))); float EsR = ER + (Ss - uR) * (Ss + pR / (rR * (SR - uR)));
fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss; fm = rsR * Ss; fp = rsR * Ss * Ss + ps; fe = (rsR * EsR + ps) * Ss;
} }
} }
// ---------- Characteristic open‑end boundaries ---------- private void OpenEndFluxLeft(float rhoInt, float uInt, float pInt, float pAmb,
private void OpenEndFluxLeft(double rhoInt, double uInt, double pInt, double pAmb, out float fm, out float fp, out float fe)
out double fm, out double fp, out double fe)
{ {
double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f));
if (uInt <= -cInt) // supersonic inflow if (uInt <= -cInt) // supersonic inflow
{ {
fm = rhoInt * uInt; fm = rhoInt * uInt;
fp = rhoInt * uInt * uInt + pInt; fp = rhoInt * uInt * uInt + pInt;
fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt; fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt;
return; return;
} }
if (uInt <= 0) // subsonic inflow if (uInt <= 0) // subsonic inflow
{ {
double T0 = 300.0, R = 287.0; float T0 = 300f, R = 287f;
double ghost_Rho = pAmb / (R * T0); float ghostRho = pAmb / (R * T0);
HLLCFlux(ghost_Rho, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); HLLCFluxScalar(ghostRho, 0f, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
return; return;
} }
// subsonic outflow // subsonic outflow
double s = pInt / Math.Pow(rhoInt, _gamma); float s = pInt / MathF.Pow(rhoInt, _gamma);
double ghostRho = Math.Pow(pAmb / s, 1.0 / _gamma); float ghostRho2 = MathF.Pow(pAmb / s, 1f / _gamma);
double cGhost = Math.Sqrt(_gamma * pAmb / ghostRho); float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho2);
double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0); float J_minus = uInt - 2f * cInt / (_gamma - 1f);
double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0); float uGhost = J_minus + 2f * cGhost / (_gamma - 1f);
if (uGhost < 0) uGhost = 0; if (uGhost < 0) uGhost = 0;
HLLCFlux(ghostRho, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe); HLLCFluxScalar(ghostRho2, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
} }
private void OpenEndFluxRight(double rhoInt, double uInt, double pInt, double pAmb, private void OpenEndFluxRight(float rhoInt, float uInt, float pInt, float pAmb,
out double fm, out double fp, out double fe) out float fm, out float fp, out float fe)
{ {
double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12)); float cInt = MathF.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12f));
if (uInt >= cInt) // supersonic outflow if (uInt >= cInt) // supersonic outflow
{ {
fm = rhoInt * uInt; fm = rhoInt * uInt;
fp = rhoInt * uInt * uInt + pInt; fp = rhoInt * uInt * uInt + pInt;
fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt; fe = (rhoInt * (pInt / ((_gamma - 1f) * rhoInt) + 0.5f * uInt * uInt) + pInt) * uInt;
return; return;
} }
if (uInt >= 0) // subsonic outflow if (uInt >= 0) // subsonic outflow
{ {
double s = pInt / Math.Pow(rhoInt, _gamma); float s = pInt / MathF.Pow(rhoInt, _gamma);
double ghost_Rho = Math.Pow(pAmb / s, 1.0 / _gamma); float ghostRho = MathF.Pow(pAmb / s, 1f / _gamma);
double cGhost = Math.Sqrt(_gamma * pAmb / ghost_Rho); float cGhost = MathF.Sqrt(_gamma * pAmb / ghostRho);
double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0); float J_plus = uInt + 2f * cInt / (_gamma - 1f);
double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0); float uGhost = J_plus - 2f * cGhost / (_gamma - 1f);
if (uGhost > 0) uGhost = 0; if (uGhost > 0) uGhost = 0;
HLLCFlux(rhoInt, uInt, pInt, ghost_Rho, uGhost, pAmb, out fm, out fp, out fe); HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho, uGhost, pAmb, out fm, out fp, out fe);
return;
} }
// subsonic inflow // subsonic inflow
double T0 = 300.0, R = 287.0; float T0 = 300f, R = 287f;
double ghostRho = pAmb / (R * T0); float ghostRho2 = pAmb / (R * T0);
HLLCFlux(rhoInt, uInt, pInt, ghostRho, 0.0, pAmb, out fm, out fp, out fe); HLLCFluxScalar(rhoInt, uInt, pInt, ghostRho2, 0f, pAmb, out fm, out fp, out fe);
} }
private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRight, private void ClosedEndFlux(float rhoInt, float uInt, float pInt, bool isRight,
out double fm, out double fp, out double fe) out float fm, out float fp, out float fe)
{ {
double rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt; float rhoGhost = rhoInt, pGhost = pInt, uGhost = -uInt;
if (isRight) if (isRight)
HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe); HLLCFluxScalar(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
else else
HLLCFlux(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe); HLLCFluxScalar(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
}
// ==================== SIMD INTERNAL FACE ROUTINE ========================
private void SimdInternalFluxBlock(int startFace, int count)
{
// startFace is the first face index; we process 'count' consecutive faces.
// For face f, left cell = f-1, right cell = f.
// We load left and right states for faces [startFace .. startFace+count-1].
int leftIdx = startFace - 1;
int rightIdx = startFace;
// Load conserved variables for left cells and right cells as vectors.
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);
// Derived quantities: u = ru / r, p = (gamma-1)*(E - 0.5*ru^2 / r)
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);
// Sound speeds
Vector<float> cL = Vector.SquareRoot(gammaVec * pL / rL);
Vector<float> cR = Vector.SquareRoot(gammaVec * pR / rR);
// Wave speeds
Vector<float> SL = Vector.Min(uL - cL, uR - cR);
Vector<float> SR = Vector.Max(uL + cL, uR + cR);
// Star speed
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;
// Total energy per unit mass (E/rho) for left/right (needed for star region)
Vector<float> eL = EL / rL;
Vector<float> eR = ER / rR;
// --- Compute all four possible flux vectors ---
// Left flux
Vector<float> Fm_L = ruL;
Vector<float> Fp_L = ruL * uL + pL;
Vector<float> Fe_L = (EL + pL) * uL; // (r*E + p)*u
// Right flux
Vector<float> Fm_R = ruR;
Vector<float> Fp_R = ruR * uR + pR;
Vector<float> Fe_R = (ER + pR) * uR;
// Star‑left fluxes (when SL < 0 < Ss)
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 (when Ss < 0 < SR)
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;
// --- Select the correct flux based on wave signs ---
var maskSLge0 = Vector.GreaterThanOrEqual(SL, Vector<float>.Zero);
var maskSRle0 = Vector.LessThanOrEqual(SR, Vector<float>.Zero);
var maskMiddle = ~(maskSLge0 | maskSRle0); // SL<0 && SR>0
var maskStarL = maskMiddle & Vector.GreaterThanOrEqual(Ss, Vector<float>.Zero);
var maskStarR = maskMiddle & Vector.LessThan(Ss, Vector<float>.Zero);
// Start with left flux, override with right/star as needed
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))));
// Store to flux arrays at indices startFace .. startFace+count-1
fm.CopyTo(_fluxM, startFace);
fp.CopyTo(_fluxP, startFace);
fe.CopyTo(_fluxE, startFace);
}
// ==================== SIMD CELL UPDATE + DAMPING ========================
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‑define constants used in clamping
Vector<float> half = new Vector<float>(0.5f);
Vector<float> gammaMinus1 = new Vector<float>(_gamma - 1f);
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> factor = Vector.Exp(-vCoeff / r * vDt);
newRu *= factor;
// 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
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 = _fluxE[i + 1] - _fluxE[i];
float newR = r - dt_dx * dM;
float newRu = ru - dt_dx * dP;
float newE = E - dt_dx * dE;
// Damping
float factor = MathF.Exp(-coeff / Math.Max(r, 1e-12f) * dt);
newRu *= factor;
// 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;
_rho[i] = newR;
_rhou[i] = newRu;
_E[i] = newE;
}
} }
} }
} }

3
Core/SoundProcessor.cs
View File

@@ -82,6 +82,7 @@ namespace FluidSim.Core
{ {
// 1. Monopole source: d(mdot)/dt // 1. Monopole source: d(mdot)/dt
double derivative = (massFlow - lastMassFlow) / dt; double derivative = (massFlow - lastMassFlow) / dt;
derivative = Math.Clamp(derivative, -500, 500);
lastMassFlow = massFlow; lastMassFlow = massFlow;
float monopole = (float)(derivative * masterGain); float monopole = (float)(derivative * masterGain);
@@ -115,7 +116,7 @@ namespace FluidSim.Core
// 6. Dry/wet mix // 6. Dry/wet mix
float output = resonant * 0.7f + wet * 0.3f; float output = resonant * 0.7f + wet * 0.3f;
output = Math.Clamp(output, -1f, 1f); output = MathF.Tanh(output);
return output; return output;
} }
} }

2
FluidSim.csproj
View File

@@ -5,7 +5,7 @@
<TargetFramework>net10.0</TargetFramework> <TargetFramework>net10.0</TargetFramework>
<ImplicitUsings>enable</ImplicitUsings> <ImplicitUsings>enable</ImplicitUsings>
<Nullable>enable</Nullable> <Nullable>enable</Nullable>
<PublishAot>true</PublishAot> <PublishAot>false</PublishAot>
<InvariantGlobalization>true</InvariantGlobalization> <InvariantGlobalization>true</InvariantGlobalization>
</PropertyGroup> </PropertyGroup>

320
Scenarios/EngineScenario.cs
View File

@@ -17,34 +17,48 @@ namespace FluidSim.Core
private double ambientPressure = 101325.0; private double ambientPressure = 101325.0;
private double time; private double time;
// Crankshaft // ---- 4‑stroke cycle angle (0 … 4π) ----
private double crankAngle = 0.0; private double cycleCrankAngle = 0.0; // 0 to 4π, then resets
private const double TargetRPM = 4000.0; private const double TargetRPM = 1000.0;
private double angularVelocity; private double angularVelocity; // rad/s of crankshaft
// Combustion // ---- Engine geometry ----
private const double CombustionPressure = 8.0 * 101325.0; private double bore = 0.065; // 65 mm
private const double CombustionTemperature = 1800.0; private double stroke = 0.0565; // 56.5 mm → 250 cc
private double conRodLength = 0.113; // roughly 2 * stroke
private double compressionRatio = 10.0;
private double V_disp; // displacement volume
private double V_clear; // clearance volume
// Valve timing // ---- Combustion ----
private const double ValveOpenStart = 120.0 * Math.PI / 180.0; private const double CombustionPressure = 50.0 * 101325.0;
private const double ValveOpenEnd = 480.0 * Math.PI / 180.0; private const double CombustionTemperature = 2500.0;
private const double ValveRampWidth = 30.0 * Math.PI / 180.0; private bool burnInProgress = false;
private double burnStartAngle; // cycle angle when ignition began
private const double BurnDurationDeg = 40.0;
private const double BurnDurationRad = BurnDurationDeg * Math.PI / 180.0;
private double targetBurnEnergy;
private double totalBurnMass;
// Pre‑ignition state (compressed fresh charge) for misfire restoration
private double preIgnitionMass;
private double preIgnitionInternalEnergy;
// ---- Valve timing ----
private const double ValveOpenStart = 120.0 * Math.PI / 180.0; // 120° after TDC power
private const double ValveOpenEnd = 480.0 * Math.PI / 180.0; // 480° ≈ 120° after TDC exhaust
private const double ValveRampWidth = 30.0 * Math.PI / 180.0; // 30° ramps
private double maxOrificeArea; private double maxOrificeArea;
// Misfire // ---- Misfire ----
private Random rand = new Random(); private Random rand = new Random();
private const double MisfireProbability = 0.02; private const double MisfireProbability = 0.02;
private bool isMisfiring = false; private bool isMisfiring = false;
// Low‑pass filter for pressure // ---- Logging ----
private double lastFilteredPressure;
private const double PressureCutoffHz = 50.0;
// Logging
private int stepCount = 0; private int stepCount = 0;
private const int LogStepInterval = 1000; private const int LogStepInterval = 10000;
private int combustionCount = 0; private int combustionCount = 0;
private int misfireCount = 0; private int misfireCount = 0;
@@ -53,25 +67,29 @@ namespace FluidSim.Core
dt = 1.0 / sampleRate; dt = 1.0 / sampleRate;
angularVelocity = TargetRPM * 2.0 * Math.PI / 60.0; angularVelocity = TargetRPM * 2.0 * Math.PI / 60.0;
// Cylinder: 0.5 litre, initially at ambient // Displacement volume
double cylVolume = 0.5e-3; V_disp = (Math.PI / 4.0) * bore * bore * stroke;
double initialPressure = ambientPressure; V_clear = V_disp / (compressionRatio - 1.0);
double initialTemperature = 300.0;
cylinder = new Volume0D(cylVolume, initialPressure, initialTemperature, sampleRate) // Cylinder (starts at TDC clearance volume with compressed ambient charge)
double initialPressure = ambientPressure * Math.Pow(compressionRatio, 1.4); // isentropic compression
double initialTemperature = 300.0 * Math.Pow(compressionRatio, 1.4 - 1.0);
double initialVolume = V_clear;
cylinder = new Volume0D(initialVolume, initialPressure, initialTemperature, sampleRate)
{ {
Gamma = 1.4, Gamma = 1.4,
GasConstant = 287.0 GasConstant = 287.0
}; };
// Exhaust pipe: length 2.5 m, radius 2 cm // Exhaust pipe (2.5 m long, 3 cm radius)
double pipeLength = 2.5; double pipeLength = 2.5;
double pipeRadius = 0.02; double pipeRadius = 0.03;
double pipeArea = Math.PI * pipeRadius * pipeRadius; double pipeArea = Math.PI * pipeRadius * pipeRadius;
maxOrificeArea = pipeArea; maxOrificeArea = pipeArea;
exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 70); exhaustPipe = new Pipe1D(pipeLength, pipeArea, sampleRate, forcedCellCount: 100);
exhaustPipe.SetUniformState(1.225, 0.0, ambientPressure); exhaustPipe.SetUniformState(1.225, 0.0, ambientPressure);
// Coupling (valve starts closed) // Coupling (valve initially closed)
coupling = new PipeVolumeConnection(cylinder, exhaustPipe, true, orificeArea: 0.0); coupling = new PipeVolumeConnection(cylinder, exhaustPipe, true, orificeArea: 0.0);
solver = new Solver(); solver = new Solver();
@@ -79,134 +97,218 @@ namespace FluidSim.Core
solver.AddVolume(cylinder); solver.AddVolume(cylinder);
solver.AddPipe(exhaustPipe); solver.AddPipe(exhaustPipe);
solver.AddConnection(coupling); solver.AddConnection(coupling);
// Use ZeroPressureOpen for strong reflections // Open end with characteristic radiation (softer reflections)
solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.ZeroPressureOpen, ambientPressure); solver.SetPipeBoundary(exhaustPipe, false, BoundaryType.OpenEnd, ambientPressure);
// Sound processor (tuned to pipe length) // Sound processor (keep your carefully tuned gains)
soundProcessor = new SoundProcessor(sampleRate, pipeLength, pipeRadius * 2, reverbTimeMs: 200.0f); soundProcessor = new SoundProcessor(sampleRate, pipeLength, pipeRadius * 2, reverbTimeMs: 200.0f);
soundProcessor.MasterGain = 0.02f; // boosted from 0.0008 soundProcessor.MasterGain = 0.0002f;
soundProcessor.PressureGain = 4.0f; // boosted from6 0.12 soundProcessor.PressureGain = 10.0f;
soundProcessor.TurbulenceGain = 0.0002f; // reduced from 0.02 soundProcessor.TurbulenceGain = 0.00005f;
soundProcessor.SetAmbientPressure(ambientPressure); soundProcessor.SetAmbientPressure(ambientPressure);
lastFilteredPressure = ambientPressure; // Log startup info
Console.WriteLine("=== EngineScenario (improved physics) ===");
Console.WriteLine("=== EngineScenario (ZeroPressureOpen, boosted gains) ===");
Console.WriteLine($"Target RPM: {TargetRPM}, Misfire prob: {MisfireProbability * 100:F1}%"); Console.WriteLine($"Target RPM: {TargetRPM}, Misfire prob: {MisfireProbability * 100:F1}%");
Console.WriteLine($"Bore x Stroke: {bore*1000:F0} x {stroke*1000:F0} mm, CR: {compressionRatio:F1}");
Console.WriteLine($"Pipe length: {pipeLength} m, fundamental: {340/(4*pipeLength):F1} Hz"); Console.WriteLine($"Pipe length: {pipeLength} m, fundamental: {340/(4*pipeLength):F1} Hz");
Console.WriteLine($"Valve opens at {ValveOpenStart*180/Math.PI:F0}°, closes at {ValveOpenEnd*180/Math.PI:F0}°, ramp {ValveRampWidth*180/Math.PI:F0}°"); Console.WriteLine($"Combustion: {CombustionPressure/101325:F0} bar, {CombustionTemperature} K");
Console.WriteLine($"Sample rate: {sampleRate} Hz, dt = {dt*1000:F3} ms"); Console.WriteLine($"Valve opens at {ValveOpenStart*180/Math.PI:F0}°, closes at {ValveOpenEnd*180/Math.PI:F0}° (ramp {ValveRampWidth*180/Math.PI:F0}°)");
Console.WriteLine("Time[s] Crank[°] Valve[%] MassFlow[kg/s] Comb# Misfire"); Console.WriteLine($"Burn duration: {BurnDurationDeg}°");
Console.WriteLine("---------------------------------------------------------"); Console.WriteLine("Time[s] Crank[°] Vol[cc] MassFlow[kg/s] Comb# Misfire");
Console.WriteLine("-------------------------------------------------------------");
} }
private double ValveOpenRatio(double crankRad) // ---- Piston volume & derivative ----
private (double volume, double dvdt) PistonKinematics(double theta)
{ {
double cycleAngle = crankRad % (4.0 * Math.PI); // theta = crankshaft angle (0 at TDC of power stroke)
double openStart = ValveOpenStart; double R = stroke / 2.0;
double openEnd = ValveOpenEnd; double cosT = Math.Cos(theta);
double sinT = Math.Sin(theta);
double L = conRodLength;
if (cycleAngle < openStart || cycleAngle > openEnd) // Slider‑crank position relative to TDC
double s = R * (1 - cosT) + L - Math.Sqrt(L * L - R * R * sinT * sinT);
double V = V_clear + (Math.PI / 4.0) * bore * bore * s;
// Derivative dV/dθ
double sqrtTerm = Math.Sqrt(L * L - R * R * sinT * sinT);
double dVdθ = (Math.PI / 4.0) * bore * bore * (R * sinT + (R * R * sinT * cosT) / sqrtTerm);
double dvdt = dVdθ * angularVelocity; // rad/s → volume change rate
return (V, dvdt);
}
// ---- Valve lift (trapezoidal) ----
private double ValveOpenRatio(double cycleRad)
{
// cycleRad: 0 … 4π
if (cycleRad < ValveOpenStart || cycleRad > ValveOpenEnd)
return 0.0; return 0.0;
double fullOpenWindow = openEnd - openStart; double duration = ValveOpenEnd - ValveOpenStart;
double closedWindow = 2.0 * ValveRampWidth; double ramp = ValveRampWidth;
if (fullOpenWindow <= closedWindow)
return 1.0;
double tmid = (openStart + openEnd) / 2.0; double t = (cycleRad - ValveOpenStart) / duration;
double dist = Math.Abs(cycleAngle - tmid);
double rampHalf = (fullOpenWindow - closedWindow) / 2.0; if (t < ramp / duration)
if (dist <= rampHalf) return t / (ramp / duration);
return 1.0; else if (t > 1.0 - ramp / duration)
return (1.0 - t) / (ramp / duration);
else else
{ return 1.0;
double frac = (dist - rampHalf) / ValveRampWidth;
frac = Math.Clamp(frac, 0.0, 1.0);
double lift = Math.Cos(frac * Math.PI / 2.0);
return lift * lift;
} }
// ---- Wiebe burn fraction ----
private double WiebeFraction(double angleFromIgnition)
{
if (angleFromIgnition >= BurnDurationRad) return 1.0;
double a = 5.0, m = 2.0;
double x = angleFromIgnition / BurnDurationRad;
return 1.0 - Math.Exp(-a * Math.Pow(x, m + 1));
} }
public override float Process() public override float Process()
{ {
// Update crank angle // Advance cycle crank angle
crankAngle += angularVelocity * dt; cycleCrankAngle += angularVelocity * dt;
if (crankAngle >= 2.0 * Math.PI) if (cycleCrankAngle >= 4.0 * Math.PI) // 720°
{ {
crankAngle -= 2.0 * Math.PI; cycleCrankAngle -= 4.0 * Math.PI;
isMisfiring = rand.NextDouble() < MisfireProbability; isMisfiring = rand.NextDouble() < MisfireProbability;
}
// Power stroke at TDC // ---- Prepare cylinder for new power stroke ----
if (crankAngle < angularVelocity * dt && crankAngle >= 0.0) // Fill cylinder with fresh charge at BDC, then compress isentropically to TDC clearance volume.
{ double T_bdc = 300.0; // intake temperature
double p_bdc = ambientPressure; // intake pressure
double V_bdc = V_clear + V_disp; // volume at BDC (intake valve closing)
double freshMass = p_bdc * V_bdc / (287.0 * T_bdc);
double freshInternalEnergy = p_bdc * V_bdc / (1.4 - 1.0);
// Compress isentropically to V_clear
double V1 = V_bdc, V2 = V_clear;
double gamma = 1.4;
double p2 = p_bdc * Math.Pow(V1 / V2, gamma);
double T2 = T_bdc * Math.Pow(V1 / V2, gamma - 1);
// Set compressed state
cylinder.Volume = V_clear;
cylinder.Mass = freshMass;
cylinder.InternalEnergy = p2 * V_clear / (gamma - 1.0); // consistent with pressure/temperature
// Store pre‑ignition state for misfire recovery
preIgnitionMass = cylinder.Mass;
preIgnitionInternalEnergy = cylinder.InternalEnergy;
if (isMisfiring) if (isMisfiring)
{ {
double vol = cylinder.Volume; // No combustion – just expand from compressed state
double R = cylinder.GasConstant;
double T0 = 300.0;
double newMass = ambientPressure * vol / (R * T0);
double newInternalEnergy = ambientPressure * vol / (cylinder.Gamma - 1.0);
cylinder.Mass = newMass;
cylinder.InternalEnergy = newInternalEnergy;
misfireCount++; misfireCount++;
} }
else else
{ {
double volume = cylinder.Volume; // Start Wiebe burn
double gamma = cylinder.Gamma; double V = V_clear;
double newInternalEnergy = CombustionPressure * volume / (gamma - 1.0); targetBurnEnergy = CombustionPressure * V / (gamma - 1.0);
double R = cylinder.GasConstant; double R = 287.0;
double newMass = CombustionPressure * volume / (R * CombustionTemperature); totalBurnMass = CombustionPressure * V / (R * CombustionTemperature);
cylinder.InternalEnergy = newInternalEnergy; burnInProgress = true;
cylinder.Mass = newMass; burnStartAngle = cycleCrankAngle; // now = 0
combustionCount++; combustionCount++;
} }
} }
// Update valve area // ---- Combustion progress (if active) ----
double valveOpen = ValveOpenRatio(crankAngle); if (burnInProgress)
coupling.OrificeArea = maxOrificeArea * valveOpen; {
double angleFromIgnition = cycleCrankAngle - burnStartAngle;
if (angleFromIgnition >= BurnDurationRad)
{
// Burn complete
cylinder.Mass = totalBurnMass;
cylinder.InternalEnergy = targetBurnEnergy;
burnInProgress = false;
}
else
{
double fraction = WiebeFraction(angleFromIgnition);
// Interpolate between pre‑ignition (compressed charge) and final burned state
double gamma = 1.4;
double V = cylinder.Volume; // still near clearance
double baseEnergy = preIgnitionInternalEnergy;
double baseMass = preIgnitionMass;
cylinder.InternalEnergy = baseEnergy * (1.0 - fraction) + targetBurnEnergy * fraction;
cylinder.Mass = baseMass * (1.0 - fraction) + totalBurnMass * fraction;
}
}
// ---- Update cylinder volume from piston kinematics ----
double theta = cycleCrankAngle % (2.0 * Math.PI); // crank angle for piston position
var (vol, dvdt) = PistonKinematics(theta);
cylinder.Volume = vol;
cylinder.Dvdt = dvdt;
// ---- Valve lift & orifice area ----
double lift = ValveOpenRatio(cycleCrankAngle);
coupling.OrificeArea = maxOrificeArea * lift;
// ---- Solver step ----
float massFlow = solver.Step(); float massFlow = solver.Step();
float endPressure = (float)exhaustPipe.GetCellPressure(exhaustPipe.GetCellCount() - 1); float endPressure = (float)exhaustPipe.GetCellPressure(exhaustPipe.GetCellCount() - 1);
// Low‑pass filter the pressure (emphasise low frequencies) // ---- Audio (no filter, feed raw pressure) ----
double rc = 1.0 / (2.0 * Math.PI * PressureCutoffHz); float audioSample = soundProcessor.Process(massFlow, endPressure);
double alpha = dt / (rc + dt);
double filteredPressure = alpha * endPressure + (1.0 - alpha) * lastFilteredPressure;
lastFilteredPressure = filteredPressure;
float audioSample = soundProcessor.Process(massFlow, (float)filteredPressure); // Log occasionally
time += dt; time += dt;
stepCount++; stepCount++;
if (stepCount % LogStepInterval == 0)
// Logging
if (stepCount % LogStepInterval == 0 || (crankAngle < angularVelocity * dt * 2 && !isMisfiring && combustionCount > 0))
{ {
Console.WriteLine($"{time,7:F3} {crankAngle * 180.0 / Math.PI,6:F1} " + double crankDeg = cycleCrankAngle * 180.0 / Math.PI;
$"{valveOpen * 100,6:F1} {massFlow,10:F4} " + double volCC = cylinder.Volume * 1e6; // cc
$"{combustionCount,3} {(isMisfiring ? "X" : "")}"); Console.WriteLine($"{time,5:F3} {crankDeg,7:F1}° {volCC,5:F1} {massFlow,14:E4} {combustionCount,4} {misfireCount,4}");
} }
return audioSample; return audioSample;
} }
// ---- Drawing (unchanged) ----
public override void Draw(RenderWindow target) public override void Draw(RenderWindow target)
{ {
float winW = target.GetView().Size.X; float winW = target.GetView().Size.X;
float winH = target.GetView().Size.Y; float winH = target.GetView().Size.Y;
float centerY = winH / 2f; float centerY = winH / 2f;
const float T_ambient = 293.15f;
const float T_hot = 1500f;
const float T_cold = 0f;
const float R = 287.05f;
float deltaHot = T_hot - T_ambient;
float deltaCold = T_ambient - T_cold;
float NormaliseTemperature(double T)
{
double t;
if (T >= T_ambient)
t = (T - T_ambient) / deltaHot;
else
t = (T - T_ambient) / deltaCold;
return (float)Math.Clamp(t, -1.0, 1.0);
}
float cylW = 80f, cylH = 150f; float cylW = 80f, cylH = 150f;
var cylRect = new RectangleShape(new Vector2f(cylW, cylH)); var cylRect = new RectangleShape(new Vector2f(cylW, cylH));
cylRect.Position = new Vector2f(40f, centerY - cylH / 2f); cylRect.Position = new Vector2f(40f, centerY - cylH / 2f);
double pNorm = (cylinder.Pressure - ambientPressure) / ambientPressure;
if (double.IsNaN(pNorm)) pNorm = 0; double tempCyl = cylinder.Temperature; // Volume0D now has Temperature
byte red = (byte)(Math.Clamp(pNorm * 128, 0, 255)); float tnCyl = NormaliseTemperature(tempCyl);
byte blue = (byte)(Math.Clamp(-pNorm * 128, 0, 255)); byte redCyl = (byte)(tnCyl > 0 ? 255 * tnCyl : 0);
cylRect.FillColor = new Color(red, 0, blue); byte blueCyl = (byte)(tnCyl < 0 ? -255 * tnCyl : 0);
byte greenCyl = (byte)(255 * (1 - Math.Abs(tnCyl)));
cylRect.FillColor = new Color(redCyl, greenCyl, blueCyl);
target.Draw(cylRect); target.Draw(cylRect);
int n = exhaustPipe.GetCellCount(); int n = exhaustPipe.GetCellCount();
@@ -215,18 +317,22 @@ namespace FluidSim.Core
float dx = pipeLen / (n - 1); float dx = pipeLen / (n - 1);
float baseRadius = 20f; float baseRadius = 20f;
var vertices = new Vertex[n * 2]; var vertices = new Vertex[n * 2];
float ambientPressure = 101325f;
for (int i = 0; i < n; i++) for (int i = 0; i < n; i++)
{ {
float x = pipeStartX + i * dx; float x = pipeStartX + i * dx;
double p = exhaustPipe.GetCellPressure(i); double p = exhaustPipe.GetCellPressure(i);
float r = baseRadius * (float)(1.0 + (p - ambientPressure) / ambientPressure); double rho = exhaustPipe.GetCellDensity(i);
double T = p / (rho * R);
float r = baseRadius * 0.1f * (float)(1.0 + (p - ambientPressure) / ambientPressure);
if (r < 2f) r = 2f; if (r < 2f) r = 2f;
double t = (p - ambientPressure) / ambientPressure; float tn = NormaliseTemperature(T);
t = Math.Clamp(t, -1.0, 1.0); byte rCol = (byte)(tn > 0 ? 255 * tn : 0);
byte rCol = (byte)(t > 0 ? 255 * t : 0); byte bCol = (byte)(tn < 0 ? -255 * tn : 0);
byte bCol = (byte)(t < 0 ? -255 * t : 0); byte gCol = (byte)(255 * (1 - Math.Abs(tn)));
byte gCol = (byte)(255 * (1 - Math.Abs(t)));
var col = new Color(rCol, gCol, bCol); var col = new Color(rCol, gCol, bCol);
vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col); vertices[i * 2] = new Vertex(new Vector2f(x, centerY - r), col);