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

Added pipe friction

Browse Source
This commit is contained in:
max
2026-05-03 02:10:28 +02:00
parent a006a07049
commit c427c1f7d3
2 changed files with 25 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
Components/Pipe1D.cs
View File

@@ -15,9 +15,10 @@ namespace FluidSim.Components
public Port PortA { get; }
public Port PortB { get; }
public double Area => _area;
public double DampingMultiplier { get; set; } = 1.0;
private int _n;
private double _dx, _dt, _gamma, _area;
private double _dx, _dt, _gamma, _area, _diameter;
private double[] _rho, _rhou, _E;
private double _rhoLeft, _pLeft;
@@ -33,8 +34,6 @@ namespace FluidSim.Components
private const double CflTarget = 0.8;
private const double ReferenceSoundSpeed = 340.0;
public double FrictionFactor { get; set; } = 0.02;
public int GetCellCount() => _n;
public double GetCellDensity(int i) => _rho[i];
public double GetCellPressure(int i) => Pressure(i);
@@ -66,6 +65,9 @@ namespace FluidSim.Components
_area = area;
_gamma = 1.4;
// Hydraulic diameter for a circular pipe
_diameter = 2.0 * Math.Sqrt(area / Math.PI);
_rho = new double[_n];
_rhou = new double[_n];
_E = new double[_n];
@@ -229,9 +231,14 @@ namespace FluidSim.Components
break;
}
// Cell update
// ---- Cell update with linear laminar damping ----
double radius = _diameter / 2.0;
double mu_air = 1.8e-5; // dynamic viscosity of air (Pa·s)
double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius);
for (int i = 0; i < n; i++)
{
// Flux divergence
double dM = (Fm[i + 1] - Fm[i]) / _dx;
double dP = (Fp[i + 1] - Fp[i]) / _dx;
double dE = (Fe[i + 1] - Fe[i]) / _dx;
@@ -240,15 +247,22 @@ namespace FluidSim.Components
_rhou[i] -= dtSub * dP;
_E[i] -= dtSub * dE;
if (_rho[i] < 1e-12) _rho[i] = 1e-12;
// Laminar viscous damping on momentum (implicit exponential decay)
double rho = Math.Max(_rho[i], 1e-12);
double dampingRate = laminarCoeff / rho; // 1/s
double dampingFactor = Math.Exp(-dampingRate * dtSub);
_rhou[i] *= dampingFactor;
// Note: total energy _E[i] is unchanged – kinetic energy loss becomes internal heat
// Physical bounds
if (_rho[i] < 1e-12) _rho[i] = 1e-12;
double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
double pMin = 100.0;
double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
if (_E[i] < eMin) _E[i] = eMin;
}
// Port quantities (only meaningful for volume coupled ends)
// Port quantities (only meaningful for volume‑coupled ends)
double mdotA_sub = _leftBCType == BoundaryType.VolumeCoupling && _leftBCSet ? Fm[0] * _area : 0.0;
double mdotB_sub = _rightBCType == BoundaryType.VolumeCoupling && _rightBCSet ? -Fm[n] * _area : 0.0;