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

Energy conservation fixed

Browse Source
This commit is contained in:
max
2026-05-02 18:05:14 +02:00
parent 650b993c2b
commit ab50b808fb
2 changed files with 107 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

169
Components/Pipe1D.cs
View File

@@ -13,10 +13,9 @@ namespace FluidSim.Components
private double _dx, _dt, _gamma = 1.4, _area; private double _dx, _dt, _gamma = 1.4, _area;
private double[] _rho, _rhou, _E; private double[] _rho, _rhou, _E;
// Boundary fluxes (set by solver before each step) // Volume states at boundaries
private double _fxL_mass, _fxL_mom, _fxL_ener; private double _rhoLeft, _pLeft, _rhoRight, _pRight;
private double _fxR_mass, _fxR_mom, _fxR_ener; private bool _leftBCSet, _rightBCSet;
private bool _leftSet, _rightSet;
public double FrictionFactor { get; set; } = 0.02; public double FrictionFactor { get; set; } = 0.02;
@@ -57,79 +56,19 @@ namespace FluidSim.Components
public double GetLeftDensity() => _rho[0]; public double GetLeftDensity() => _rho[0];
public double GetRightDensity() => _rho[_n - 1]; public double GetRightDensity() => _rho[_n - 1];
public void SetLeftBoundaryFlux(double m, double p, double e) // ★ New: pass both density and pressure from the volume
public void SetLeftVolumeState(double rhoVol, double pVol)
{ {
_fxL_mass = m; _fxL_mom = p; _fxL_ener = e; _leftSet = true; _rhoLeft = rhoVol;
_pLeft = pVol;
_leftBCSet = true;
} }
public void SetRightBoundaryFlux(double m, double p, double e) public void SetRightVolumeState(double rhoVol, double pVol)
{ {
_fxR_mass = m; _fxR_mom = p; _fxR_ener = e; _rightSet = true; _rhoRight = rhoVol;
} _pRight = pVol;
_rightBCSet = true;
public void Simulate()
{
int n = _n;
double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1];
// Left face
if (_leftSet) { Fm[0] = _fxL_mass; Fp[0] = _fxL_mom; Fe[0] = _fxL_ener; }
else { Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0; } // reflective wall
// Internal faces (HLLC)
for (int i = 0; i < n - 1; i++)
{
double uL = _rhou[i] / Math.Max(_rho[i], 1e-12);
double uR = _rhou[i + 1] / Math.Max(_rho[i + 1], 1e-12);
HLLCFlux(_rho[i], uL, Pressure(i), _rho[i + 1], uR, Pressure(i + 1),
out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
}
// Right face
if (_rightSet) { Fm[n] = _fxR_mass; Fp[n] = _fxR_mom; Fe[n] = _fxR_ener; }
else { Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0; }
// Update cells
for (int i = 0; i < n; i++)
{
double dM = (Fm[i + 1] - Fm[i]) / _dx;
double dP = (Fp[i + 1] - Fp[i]) / _dx;
double dE = (Fe[i + 1] - Fe[i]) / _dx;
_rho[i] -= _dt * dM;
_rhou[i] -= _dt * dP;
_E[i] -= _dt * dE;
// Clamp to physical
if (_rho[i] < 1e-12) _rho[i] = 1e-12;
double u = _rhou[i] / _rho[i];
double kinetic = 0.5 * _rho[i] * u * u;
if (_E[i] < kinetic) _E[i] = kinetic;
}
// Friction (energy‑conserving)
if (FrictionFactor > 0)
{
double D = 2.0 * Math.Sqrt(_area / Math.PI);
for (int i = 0; i < _n; i++)
{
double u = _rhou[i] / Math.Max(_rho[i], 1e-12);
double f = FrictionFactor / (2.0 * D) * _rho[i] * u * Math.Abs(u);
//_rhou[i] -= _dt * f; FRICTIN DISABLED!!!
}
}
// Write port flows for the solver
PortA.MassFlowRate = _leftSet ? _fxL_mass * _area : 0.0;
PortB.MassFlowRate = _rightSet ? -_fxR_mass * _area : 0.0;
// Enthalpy for upwinding
PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0);
PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1);
// Reset for next step
_leftSet = _rightSet = false;
} }
private double GetCellTotalSpecificEnthalpy(int i) private double GetCellTotalSpecificEnthalpy(int i)
@@ -141,6 +80,90 @@ namespace FluidSim.Components
return h + 0.5 * u * u; return h + 0.5 * u * u;
} }
public void Simulate()
{
int n = _n;
double[] Fm = new double[n + 1], Fp = new double[n + 1], Fe = new double[n + 1];
// --- Left boundary (face 0) ---
if (_leftBCSet)
{
// Ghost = actual volume state (ρ_vol, u=0, p_vol)
double rhoL = _rhoLeft;
double uL = 0.0;
double pL = _pLeft;
double rhoR = _rho[0];
double uR = _rhou[0] / Math.Max(rhoR, 1e-12);
double pR = Pressure(0);
HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[0], out Fp[0], out Fe[0]);
}
else
{
Fm[0] = 0;
Fp[0] = Pressure(0);
Fe[0] = 0;
}
// --- Internal faces ---
for (int i = 0; i < n - 1; i++)
{
double uL = _rhou[i] / Math.Max(_rho[i], 1e-12);
double uR = _rhou[i + 1] / Math.Max(_rho[i + 1], 1e-12);
HLLCFlux(_rho[i], uL, Pressure(i), _rho[i + 1], uR, Pressure(i + 1),
out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
}
// --- Right boundary (face n) ---
if (_rightBCSet)
{
double rhoL = _rho[n - 1];
double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
double pL = Pressure(n - 1);
// Ghost = actual volume state (ρ_vol, u=0, p_vol)
double rhoR = _rhoRight;
double uR = 0.0;
double pR = _pRight;
HLLCFlux(rhoL, uL, pL, rhoR, uR, pR, out Fm[n], out Fp[n], out Fe[n]);
}
else
{
Fm[n] = 0;
Fp[n] = Pressure(n - 1);
Fe[n] = 0;
}
// --- Cell update ---
for (int i = 0; i < n; i++)
{
double dM = (Fm[i + 1] - Fm[i]) / _dx;
double dP = (Fp[i + 1] - Fp[i]) / _dx;
double dE = (Fe[i + 1] - Fe[i]) / _dx;
_rho[i] -= _dt * dM;
_rhou[i] -= _dt * dP;
_E[i] -= _dt * dE;
if (_rho[i] < 1e-12) _rho[i] = 1e-12;
double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
if (_E[i] < kinetic) _E[i] = kinetic;
}
// --- Friction disabled ---
// if (FrictionFactor > 0) { … }
// --- Port flows ---
PortA.MassFlowRate = _leftBCSet ? Fm[0] * _area : 0.0;
PortB.MassFlowRate = _rightBCSet ? -Fm[n] * _area : 0.0;
PortA.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(0);
PortB.SpecificEnthalpy = GetCellTotalSpecificEnthalpy(_n - 1);
_leftBCSet = _rightBCSet = false;
}
double Pressure(int i) => double Pressure(int i) =>
(_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12)); (_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));

59
Core/Solver.cs
View File

@@ -16,26 +16,24 @@ namespace FluidSim.Core
public void Step() public void Step()
{ {
// 1. Volumes publish state // 1. Volumes publish state to their ports
foreach (var v in _volumes) foreach (var v in _volumes)
v.PushStateToPort(); v.PushStateToPort();
// 2. Compute boundary fluxes (orifice model) // 2. Set volume states as boundary conditions on pipes
foreach (var conn in _connections) foreach (var conn in _connections)
{ {
if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
ApplyOrifice(conn, conn.PortA, conn.PortB); SetVolumeBC(conn.PortA, conn.PortB);
else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB)) else if (IsVolumePort(conn.PortA) && IsPipePort(conn.PortB))
ApplyOrifice(conn, conn.PortB, conn.PortA); SetVolumeBC(conn.PortB, conn.PortA);
else if (IsVolumePort(conn.PortA) && IsVolumePort(conn.PortB))
VolumeToVolume(conn);
} }
// 3. Pipe simulation step // 3. Run pipe simulations
foreach (var p in _pipes) foreach (var p in _pipes)
p.Simulate(); p.Simulate();
// 4. Transfer pipe‑port data to volumes // 4. Transfer pipe‑port flows to volume ports
foreach (var conn in _connections) foreach (var conn in _connections)
{ {
if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB)) if (IsPipePort(conn.PortA) && IsVolumePort(conn.PortB))
@@ -53,64 +51,29 @@ namespace FluidSim.Core
bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p); bool IsPipePort(Port p) => _pipes.Exists(pp => pp.PortA == p || pp.PortB == p);
Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p); Pipe1D GetPipe(Port p) => _pipes.Find(pp => pp.PortA == p || pp.PortB == p);
void ApplyOrifice(Connection conn, Port pipePort, Port volPort) void SetVolumeBC(Port pipePort, Port volPort)
{ {
Pipe1D pipe = GetPipe(pipePort); Pipe1D pipe = GetPipe(pipePort);
if (pipe == null) return; if (pipe == null) return;
bool isLeft = pipe.PortA == pipePort; bool isLeft = pipe.PortA == pipePort;
double pP = isLeft ? pipe.GetLeftPressure() : pipe.GetRightPressure();
double rhoP = isLeft ? pipe.GetLeftDensity() : pipe.GetRightDensity();
double uP = isLeft ? pipe.GetCellVelocity(0)
: pipe.GetCellVelocity(pipe.GetCellCount() - 1);
double pV = volPort.Pressure;
double rhoV = volPort.Density;
double uV = 0.0; // volume has zero organized velocity
OrificeBoundary.PipeVolumeFlux(
pP, rhoP, uP,
pV, rhoV, uV,
conn, pipe.Area, isLeft,
out double massFlux, out double momFlux, out double energyFlux);
if (isLeft) if (isLeft)
pipe.SetLeftBoundaryFlux(massFlux, momFlux, energyFlux); pipe.SetLeftVolumeState(volPort.Density, volPort.Pressure);
else else
pipe.SetRightBoundaryFlux(massFlux, momFlux, energyFlux); pipe.SetRightVolumeState(volPort.Density, volPort.Pressure);
}
void VolumeToVolume(Connection conn)
{
double mdot = OrificeBoundary.MassFlow(
conn.PortA.Pressure, conn.PortA.Density,
conn.PortB.Pressure, conn.PortB.Density, conn);
conn.PortA.MassFlowRate = -mdot;
conn.PortB.MassFlowRate = mdot;
if (mdot > 0)
conn.PortB.SpecificEnthalpy = conn.PortA.SpecificEnthalpy;
else if (mdot < 0)
conn.PortA.SpecificEnthalpy = conn.PortB.SpecificEnthalpy;
} }
void TransferPipeToVolume(Port pipePort, Port volPort) void TransferPipeToVolume(Port pipePort, Port volPort)
{ {
double mdot = pipePort.MassFlowRate; double mdot = pipePort.MassFlowRate;
// mdot > 0 → fluid enters pipe from volume
// mdot < 0 → fluid leaves pipe and enters volume
// Volume mass flow sign is opposite (positive into volume)
volPort.MassFlowRate = -mdot; volPort.MassFlowRate = -mdot;
if (mdot < 0) // pipe → volume if (mdot < 0) // pipe → volume
{ {
// ★ pipePort.SpecificEnthalpy now contains TOTAL enthalpy // pipePort.SpecificEnthalpy is already total (h + ½u²)
volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy; volPort.SpecificEnthalpy = pipePort.SpecificEnthalpy;
} }
// else: fluid goes volume → pipe → volume owns its own (static) enthalpy, // else: volume → pipe, volume’s own static enthalpy is used (already set)
// which is already correct.
} }
} }
} }