295 lines
11 KiB
C#
295 lines
11 KiB
C#
using System;
|
||
using FluidSim.Interfaces;
|
||
|
||
namespace FluidSim.Components
|
||
{
|
||
/// <summary>
|
||
/// 1‑D compressible Euler pipe with Lax‑Friedrichs finite‑volume scheme.
|
||
/// Ghost states are set externally via SetGhostLeft/Right; they are always required.
|
||
/// </summary>
|
||
public class Pipe1D : IComponent
|
||
{
|
||
public Port PortA { get; }
|
||
public Port PortB { get; }
|
||
public double Area { get; }
|
||
public double DampingMultiplier { get; set; } = 1.0;
|
||
public double EnergyRelaxationRate { get; set; } = 0.0; // 1/s
|
||
|
||
private double _ambientPressure = 101325.0;
|
||
public double AmbientPressure
|
||
{
|
||
get => _ambientPressure;
|
||
set
|
||
{
|
||
_ambientPressure = value;
|
||
_ambientEnergyReference = value / (_gamma - 1.0);
|
||
}
|
||
}
|
||
|
||
private readonly int _n;
|
||
private readonly double _dx;
|
||
private readonly double _diameter;
|
||
private readonly double _gamma = 1.4;
|
||
|
||
private double[] _rho, _rhou, _E;
|
||
private double[] _fluxM, _fluxP, _fluxE; // flux at cell faces (0.._n)
|
||
|
||
private double _rhoGhostL, _uGhostL, _pGhostL;
|
||
private double _rhoGhostR, _uGhostR, _pGhostR;
|
||
private bool _ghostLValid, _ghostRValid;
|
||
|
||
private double _laminarCoeff;
|
||
private double _ambientEnergyReference;
|
||
|
||
public Pipe1D(double length, double area, int cellCount)
|
||
{
|
||
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];
|
||
|
||
double mu_air = 1.8e-5;
|
||
double radius = _diameter * 0.5;
|
||
_laminarCoeff = 8.0 * mu_air / (radius * radius);
|
||
|
||
_ambientEnergyReference = 101325.0 / (_gamma - 1.0);
|
||
|
||
PortA = new Port { Owner = this };
|
||
PortB = new Port { Owner = this };
|
||
|
||
SetUniformState(1.225, 0.0, 101325.0);
|
||
}
|
||
|
||
IReadOnlyList<Port> IComponent.Ports => new[] { PortA, PortB };
|
||
|
||
public void UpdateState(double dt) { }
|
||
|
||
// ---------- Ghost interface ----------
|
||
public void SetGhostLeft(double rho, double u, double p)
|
||
{
|
||
_rhoGhostL = rho; _uGhostL = u; _pGhostL = p; _ghostLValid = true;
|
||
}
|
||
public void SetGhostRight(double rho, double u, double p)
|
||
{
|
||
_rhoGhostR = rho; _uGhostR = u; _pGhostR = p; _ghostRValid = true;
|
||
}
|
||
public void ClearGhostFlags() { _ghostLValid = false; _ghostRValid = false; }
|
||
|
||
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) => _rhou[i] / Math.Max(_rho[i], 1e-12);
|
||
public double GetCellPressure(int i) => PressureScalar(i);
|
||
|
||
public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
|
||
{
|
||
double maxW = 0.0;
|
||
for (int i = 0; i < _n; i++)
|
||
{
|
||
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-8);
|
||
return Math.Max(1, (int)Math.Ceiling(dtGlobal * maxW / (cflTarget * _dx)));
|
||
}
|
||
|
||
// ---------- Main step (per sub‑step) ----------
|
||
public void SimulateSingleStep(double dtSub)
|
||
{
|
||
if (!_ghostLValid || !_ghostRValid)
|
||
throw new InvalidOperationException("Ghost cells not set before SimulateSingleStep.");
|
||
|
||
double dt = dtSub;
|
||
int n = _n;
|
||
|
||
// ---- Compute fluxes at all faces using Lax‑Friedrichs ----
|
||
// Left face (0): between ghostL and cell 0
|
||
double rL = Math.Max(_rhoGhostL, 1e-12);
|
||
double pL = _pGhostL;
|
||
double uL = _uGhostL;
|
||
double eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;
|
||
|
||
double rR = Math.Max(_rho[0], 1e-12);
|
||
double pR = PressureScalar(0);
|
||
double uR = _rhou[0] / rR;
|
||
double eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
|
||
|
||
LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR,
|
||
out _fluxM[0], out _fluxP[0], out _fluxE[0]);
|
||
|
||
// Internal faces (1 .. n-1)
|
||
for (int f = 1; f < n; f++)
|
||
{
|
||
int iL = f - 1;
|
||
int iR = f;
|
||
|
||
rL = Math.Max(_rho[iL], 1e-12);
|
||
pL = PressureScalar(iL);
|
||
uL = _rhou[iL] / rL;
|
||
eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;
|
||
|
||
rR = Math.Max(_rho[iR], 1e-12);
|
||
pR = PressureScalar(iR);
|
||
uR = _rhou[iR] / rR;
|
||
eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
|
||
|
||
LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR,
|
||
out _fluxM[f], out _fluxP[f], out _fluxE[f]);
|
||
}
|
||
|
||
// Right face (n): between cell n-1 and ghostR
|
||
rL = Math.Max(_rho[n - 1], 1e-12);
|
||
pL = PressureScalar(n - 1);
|
||
uL = _rhou[n - 1] / rL;
|
||
eL = pL / ((_gamma - 1.0) * rL) + 0.5 * uL * uL;
|
||
|
||
rR = Math.Max(_rhoGhostR, 1e-12);
|
||
pR = _pGhostR;
|
||
uR = _uGhostR;
|
||
eR = pR / ((_gamma - 1.0) * rR) + 0.5 * uR * uR;
|
||
|
||
LaxFriedrichsFlux(rL, uL, pL, eL, rR, uR, pR, eR,
|
||
out _fluxM[n], out _fluxP[n], out _fluxE[n]);
|
||
|
||
// ---- Cell update ----
|
||
double dt_dx = dt / _dx;
|
||
double coeff = _laminarCoeff * DampingMultiplier;
|
||
double relaxRate = EnergyRelaxationRate;
|
||
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
double r = _rho[i];
|
||
double ru = _rhou[i];
|
||
double E = _E[i];
|
||
|
||
double dM = _fluxM[i + 1] - _fluxM[i];
|
||
double dP = _fluxP[i + 1] - _fluxP[i];
|
||
double dE_flux = _fluxE[i + 1] - _fluxE[i];
|
||
|
||
double newR = r - dt_dx * dM;
|
||
double newRu = ru - dt_dx * dP;
|
||
double newE = E - dt_dx * dE_flux;
|
||
|
||
double dampingFactor = Math.Exp(-coeff / Math.Max(r, 1e-12) * dt);
|
||
newRu *= dampingFactor;
|
||
|
||
double relaxFactor = Math.Exp(-relaxRate * dt);
|
||
newE = _ambientEnergyReference + (newE - _ambientEnergyReference) * relaxFactor;
|
||
|
||
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
|
||
(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;
|
||
}
|
||
|
||
// ---------- Lax‑Friedrichs flux ----------
|
||
private void LaxFriedrichsFlux(double rL, double uL, double pL, double eL,
|
||
double rR, double uR, double pR, double eR,
|
||
out double fm, out double fp, out double fe)
|
||
{
|
||
// Primitive states
|
||
double rhoL = rL, rhoR = rR;
|
||
double EL = rhoL * eL; // total energy per volume = rho * (specific total energy)
|
||
double ER = rhoR * eR;
|
||
|
||
// Conserved vectors U = (ρ, ρu, E)
|
||
// Flux F = (ρu, ρu²+p, (E+p)u)
|
||
double Fm_L = rhoL * uL;
|
||
double Fp_L = rhoL * uL * uL + pL;
|
||
double Fe_L = (EL + pL) * uL;
|
||
|
||
double Fm_R = rhoR * uR;
|
||
double Fp_R = rhoR * uR * uR + pR;
|
||
double Fe_R = (ER + pR) * uR;
|
||
|
||
// Lax‑Friedrichs dissipation coefficient α = max(|u|+c) over whole domain, but here we use local max to be simple:
|
||
double cL = Math.Sqrt(_gamma * pL / rL);
|
||
double cR = Math.Sqrt(_gamma * pR / rR);
|
||
double alpha = Math.Max(Math.Abs(uL) + cL, Math.Abs(uR) + cR);
|
||
|
||
fm = 0.5 * (Fm_L + Fm_R) - 0.5 * alpha * (rhoR - rhoL);
|
||
fp = 0.5 * (Fp_L + Fp_R) - 0.5 * alpha * (rhoR * uR - rhoL * uL);
|
||
fe = 0.5 * (Fe_L + Fe_R) - 0.5 * alpha * (ER - EL);
|
||
}
|
||
|
||
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);
|
||
}
|
||
|
||
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;
|
||
}
|
||
}
|
||
|
||
public void SetCellState(int i, double rho, double u, double p)
|
||
{
|
||
if (i < 0 || i >= _n) return;
|
||
double e = p / ((_gamma - 1.0) * rho);
|
||
double E = rho * e + 0.5 * rho * u * u;
|
||
_rho[i] = rho;
|
||
_rhou[i] = rho * u;
|
||
_E[i] = E;
|
||
}
|
||
|
||
public void SetCellPressure(int i, double p)
|
||
{
|
||
if (i < 0 || i >= _n) return;
|
||
double rho = _rho[i];
|
||
double u = _rhou[i] / rho;
|
||
double e = p / ((_gamma - 1.0) * rho);
|
||
_E[i] = rho * e + 0.5 * rho * u * u;
|
||
}
|
||
}
|
||
} |