Helmholtz test, sod shock tube
This commit is contained in:
@@ -21,15 +21,16 @@ namespace FluidSim.Components
|
||||
private double _dx, _dt, _gamma, _area, _diameter;
|
||||
private double[] _rho, _rhou, _E;
|
||||
|
||||
private double _rhoLeft, _pLeft;
|
||||
private double _rhoRight, _pRight;
|
||||
private bool _leftBCSet, _rightBCSet;
|
||||
// Volume‑coupling ghost states for boundaries A and B
|
||||
private double _rhoA, _pA;
|
||||
private double _rhoB, _pB;
|
||||
private bool _aBCSet, _bBCSet;
|
||||
|
||||
private BoundaryType _leftBCType = BoundaryType.VolumeCoupling;
|
||||
private BoundaryType _rightBCType = BoundaryType.VolumeCoupling;
|
||||
private BoundaryType _aBCType = BoundaryType.VolumeCoupling;
|
||||
private BoundaryType _bBCType = BoundaryType.VolumeCoupling;
|
||||
|
||||
private double _leftAmbientPressure = 101325.0;
|
||||
private double _rightAmbientPressure = 101325.0;
|
||||
private double _aAmbientPressure = 101325.0;
|
||||
private double _bAmbientPressure = 101325.0;
|
||||
|
||||
private const double CflTarget = 0.8;
|
||||
private const double ReferenceSoundSpeed = 340.0;
|
||||
@@ -39,8 +40,8 @@ namespace FluidSim.Components
|
||||
public double GetCellPressure(int i) => Pressure(i);
|
||||
public double GetCellVelocity(int i) => _rhou[i] / Math.Max(_rho[i], 1e-12);
|
||||
|
||||
public BoundaryType LeftBCType => _leftBCType;
|
||||
public BoundaryType RightBCType => _rightBCType;
|
||||
public BoundaryType ABCType => _aBCType;
|
||||
public BoundaryType BBCType => _bBCType;
|
||||
|
||||
public Pipe1D(double length, double area, int sampleRate, int forcedCellCount = 0)
|
||||
{
|
||||
@@ -76,10 +77,10 @@ namespace FluidSim.Components
|
||||
PortB = new Port();
|
||||
}
|
||||
|
||||
public void SetLeftBoundaryType(BoundaryType type) => _leftBCType = type;
|
||||
public void SetRightBoundaryType(BoundaryType type) => _rightBCType = type;
|
||||
public void SetLeftAmbientPressure(double p) => _leftAmbientPressure = p;
|
||||
public void SetRightAmbientPressure(double p) => _rightAmbientPressure = p;
|
||||
public void SetABoundaryType(BoundaryType type) => _aBCType = type;
|
||||
public void SetBBoundaryType(BoundaryType type) => _bBCType = type;
|
||||
public void SetAAmbientPressure(double p) => _aAmbientPressure = p;
|
||||
public void SetBAmbientPressure(double p) => _bAmbientPressure = p;
|
||||
|
||||
public void SetUniformState(double rho, double u, double p)
|
||||
{
|
||||
@@ -102,21 +103,21 @@ namespace FluidSim.Components
|
||||
_E[i] = rho * e + 0.5 * rho * u * u;
|
||||
}
|
||||
|
||||
public void SetLeftVolumeState(double rhoVol, double pVol)
|
||||
public void SetAVolumeState(double rhoVol, double pVol)
|
||||
{
|
||||
_rhoLeft = rhoVol;
|
||||
_pLeft = pVol;
|
||||
_leftBCSet = true;
|
||||
_rhoA = rhoVol;
|
||||
_pA = pVol;
|
||||
_aBCSet = true;
|
||||
}
|
||||
|
||||
public void SetRightVolumeState(double rhoVol, double pVol)
|
||||
public void SetBVolumeState(double rhoVol, double pVol)
|
||||
{
|
||||
_rhoRight = rhoVol;
|
||||
_pRight = pVol;
|
||||
_rightBCSet = true;
|
||||
_rhoB = rhoVol;
|
||||
_pB = pVol;
|
||||
_bBCSet = true;
|
||||
}
|
||||
|
||||
public void ClearBC() => _leftBCSet = _rightBCSet = false;
|
||||
public void ClearBC() => _aBCSet = _bBCSet = false;
|
||||
|
||||
public int GetRequiredSubSteps(double dtGlobal, double cflTarget = 0.8)
|
||||
{
|
||||
@@ -140,105 +141,90 @@ namespace FluidSim.Components
|
||||
double[] Fp = new double[n + 1];
|
||||
double[] Fe = new double[n + 1];
|
||||
|
||||
// Left boundary (face 0)
|
||||
switch (_leftBCType)
|
||||
// ---------- Boundary A (face 0, left) ----------
|
||||
double rhoIntA = _rho[0];
|
||||
double uIntA = _rhou[0] / Math.Max(rhoIntA, 1e-12);
|
||||
double pIntA = Pressure(0);
|
||||
|
||||
switch (_aBCType)
|
||||
{
|
||||
case BoundaryType.VolumeCoupling:
|
||||
if (_leftBCSet)
|
||||
if (_aBCSet)
|
||||
{
|
||||
HLLCFlux(_rhoLeft, 0.0, _pLeft,
|
||||
_rho[0], _rhou[0] / Math.Max(_rho[0], 1e-12), Pressure(0),
|
||||
HLLCFlux(_rhoA, 0.0, _pA,
|
||||
rhoIntA, uIntA, pIntA,
|
||||
out Fm[0], out Fp[0], out Fe[0]);
|
||||
}
|
||||
else
|
||||
{
|
||||
Fm[0] = 0; Fp[0] = Pressure(0); Fe[0] = 0;
|
||||
Fm[0] = 0; Fp[0] = pIntA; Fe[0] = 0;
|
||||
}
|
||||
break;
|
||||
|
||||
case BoundaryType.OpenEnd:
|
||||
{
|
||||
double rhoR = _rho[0];
|
||||
double uR = _rhou[0] / Math.Max(rhoR, 1e-12);
|
||||
double pR = Pressure(0);
|
||||
HLLCFlux(rhoR, uR, _leftAmbientPressure,
|
||||
rhoR, uR, pR,
|
||||
OpenEndFluxA(rhoIntA, uIntA, pIntA, _aAmbientPressure,
|
||||
out Fm[0], out Fp[0], out Fe[0]);
|
||||
}
|
||||
break;
|
||||
|
||||
case BoundaryType.ClosedEnd:
|
||||
{
|
||||
double rhoR = _rho[0];
|
||||
double uR = _rhou[0] / Math.Max(rhoR, 1e-12);
|
||||
double pR = Pressure(0);
|
||||
HLLCFlux(rhoR, -uR, pR,
|
||||
rhoR, uR, pR,
|
||||
out Fm[0], out Fp[0], out Fe[0]);
|
||||
}
|
||||
ClosedEndFlux(rhoIntA, uIntA, pIntA, isRightBoundary: false,
|
||||
out Fm[0], out Fp[0], out Fe[0]);
|
||||
break;
|
||||
}
|
||||
|
||||
// Internal faces
|
||||
// ---------- 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),
|
||||
double rhoL = _rho[i];
|
||||
double uL = _rhou[i] / Math.Max(rhoL, 1e-12);
|
||||
double pL = Pressure(i);
|
||||
|
||||
double rhoR = _rho[i + 1];
|
||||
double uR = _rhou[i + 1] / Math.Max(rhoR, 1e-12);
|
||||
double pR = Pressure(i + 1);
|
||||
|
||||
HLLCFlux(rhoL, uL, pL, rhoR, uR, pR,
|
||||
out Fm[i + 1], out Fp[i + 1], out Fe[i + 1]);
|
||||
}
|
||||
|
||||
// Right boundary (face n)
|
||||
switch (_rightBCType)
|
||||
// ---------- Boundary B (face n, right) ----------
|
||||
double rhoIntB = _rho[n - 1];
|
||||
double uIntB = _rhou[n - 1] / Math.Max(rhoIntB, 1e-12);
|
||||
double pIntB = Pressure(n - 1);
|
||||
|
||||
switch (_bBCType)
|
||||
{
|
||||
case BoundaryType.VolumeCoupling:
|
||||
if (_rightBCSet)
|
||||
if (_bBCSet)
|
||||
{
|
||||
double rhoL = _rho[n - 1];
|
||||
double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
|
||||
double pL = Pressure(n - 1);
|
||||
HLLCFlux(rhoL, uL, pL,
|
||||
_rhoRight, 0.0, _pRight,
|
||||
HLLCFlux(rhoIntB, uIntB, pIntB,
|
||||
_rhoB, 0.0, _pB,
|
||||
out Fm[n], out Fp[n], out Fe[n]);
|
||||
}
|
||||
else
|
||||
{
|
||||
Fm[n] = 0; Fp[n] = Pressure(n - 1); Fe[n] = 0;
|
||||
Fm[n] = 0; Fp[n] = pIntB; Fe[n] = 0;
|
||||
}
|
||||
break;
|
||||
|
||||
case BoundaryType.OpenEnd:
|
||||
{
|
||||
double rhoL = _rho[n - 1];
|
||||
double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
|
||||
double pL = Pressure(n - 1);
|
||||
HLLCFlux(rhoL, uL, pL,
|
||||
rhoL, uL, _rightAmbientPressure,
|
||||
OpenEndFluxB(rhoIntB, uIntB, pIntB, _bAmbientPressure,
|
||||
out Fm[n], out Fp[n], out Fe[n]);
|
||||
}
|
||||
break;
|
||||
|
||||
case BoundaryType.ClosedEnd:
|
||||
{
|
||||
double rhoL = _rho[n - 1];
|
||||
double uL = _rhou[n - 1] / Math.Max(rhoL, 1e-12);
|
||||
double pL = Pressure(n - 1);
|
||||
HLLCFlux(rhoL, uL, pL,
|
||||
rhoL, -uL, pL,
|
||||
out Fm[n], out Fp[n], out Fe[n]);
|
||||
}
|
||||
ClosedEndFlux(rhoIntB, uIntB, pIntB, isRightBoundary: true,
|
||||
out Fm[n], out Fp[n], out Fe[n]);
|
||||
break;
|
||||
}
|
||||
|
||||
// ---- Cell update with linear laminar damping ----
|
||||
double radius = _diameter / 2.0;
|
||||
double mu_air = 1.8e-5; // dynamic viscosity of air (Pa·s)
|
||||
double mu_air = 1.8e-5;
|
||||
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;
|
||||
@@ -247,14 +233,11 @@ namespace FluidSim.Components
|
||||
_rhou[i] -= dtSub * dP;
|
||||
_E[i] -= dtSub * dE;
|
||||
|
||||
// Laminar viscous damping on momentum (implicit exponential decay)
|
||||
double rho = Math.Max(_rho[i], 1e-12);
|
||||
double dampingRate = laminarCoeff / rho; // 1/s
|
||||
double dampingRate = laminarCoeff / rho;
|
||||
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;
|
||||
@@ -262,28 +245,29 @@ namespace FluidSim.Components
|
||||
if (_E[i] < eMin) _E[i] = eMin;
|
||||
}
|
||||
|
||||
// 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;
|
||||
// ---------- Port quantities ----------
|
||||
double mdotA_sub = _aBCType == BoundaryType.VolumeCoupling && _aBCSet ? Fm[0] * _area : 0.0;
|
||||
double mdotB_sub = _bBCType == BoundaryType.VolumeCoupling && _bBCSet ? -Fm[n] * _area : 0.0;
|
||||
|
||||
PortA.MassFlowRate = mdotA_sub;
|
||||
PortB.MassFlowRate = mdotB_sub;
|
||||
PortA.Pressure = Pressure(0);
|
||||
PortB.Pressure = Pressure(_n - 1);
|
||||
PortA.Pressure = pIntA;
|
||||
PortB.Pressure = pIntB;
|
||||
PortA.Density = _rho[0];
|
||||
PortB.Density = _rho[_n - 1];
|
||||
PortB.Density = _rho[n - 1];
|
||||
|
||||
if (_leftBCType == BoundaryType.VolumeCoupling && _leftBCSet)
|
||||
// Corrected enthalpy for both directions
|
||||
if (_aBCType == BoundaryType.VolumeCoupling && _aBCSet)
|
||||
{
|
||||
PortA.SpecificEnthalpy = mdotA_sub < 0
|
||||
? GetCellTotalSpecificEnthalpy(0)
|
||||
: 0.0;
|
||||
: (_gamma / (_gamma - 1.0)) * _pA / Math.Max(_rhoA, 1e-12);
|
||||
}
|
||||
if (_rightBCType == BoundaryType.VolumeCoupling && _rightBCSet)
|
||||
if (_bBCType == BoundaryType.VolumeCoupling && _bBCSet)
|
||||
{
|
||||
PortB.SpecificEnthalpy = mdotB_sub < 0
|
||||
? GetCellTotalSpecificEnthalpy(_n - 1)
|
||||
: 0.0;
|
||||
: (_gamma / (_gamma - 1.0)) * _pB / Math.Max(_rhoB, 1e-12);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -299,6 +283,98 @@ namespace FluidSim.Components
|
||||
private double Pressure(int i) =>
|
||||
(_gamma - 1.0) * (_E[i] - 0.5 * _rhou[i] * _rhou[i] / Math.Max(_rho[i], 1e-12));
|
||||
|
||||
// ========== Characteristic‑based Open End ==========
|
||||
private void OpenEndFluxA(double rhoInt, double uInt, double pInt, double pAmb,
|
||||
out double fm, out double fp, out double fe)
|
||||
{
|
||||
double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
|
||||
|
||||
// Subsonic inflow (uInt ≤ 0, so flow inside pipe ←)
|
||||
if (uInt <= -cInt) // supersonic inflow – use interior state as ghost
|
||||
{
|
||||
fm = rhoInt * uInt;
|
||||
fp = rhoInt * uInt * uInt + pInt;
|
||||
fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
|
||||
return;
|
||||
}
|
||||
else if (uInt <= 0) // subsonic inflow
|
||||
{
|
||||
// Reservoir condition: p = pAmb, T = 300K, u = 0
|
||||
double T0 = 300.0;
|
||||
double R = 287.0;
|
||||
double rhoGhost = pAmb / (R * T0);
|
||||
HLLCFlux(rhoGhost, 0.0, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||||
return;
|
||||
}
|
||||
else // subsonic outflow (uInt > 0)
|
||||
{
|
||||
// Ghost pressure forced to pAmb
|
||||
double s = pInt / Math.Pow(rhoInt, _gamma);
|
||||
double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
|
||||
double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
|
||||
|
||||
// Outgoing Riemann invariant J⁻ = uInt - 2*cInt/(γ-1) (for left boundary)
|
||||
double J_minus = uInt - 2.0 * cInt / (_gamma - 1.0);
|
||||
double uGhost = J_minus + 2.0 * cGhost / (_gamma - 1.0);
|
||||
|
||||
// Prevent spurious inflow by clipping to zero
|
||||
if (uGhost < 0) uGhost = 0;
|
||||
|
||||
HLLCFlux(rhoGhost, uGhost, pAmb, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||||
}
|
||||
}
|
||||
|
||||
private void OpenEndFluxB(double rhoInt, double uInt, double pInt, double pAmb,
|
||||
out double fm, out double fp, out double fe)
|
||||
{
|
||||
double cInt = Math.Sqrt(_gamma * pInt / Math.Max(rhoInt, 1e-12));
|
||||
|
||||
if (uInt >= cInt) // supersonic outflow (extrapolation)
|
||||
{
|
||||
fm = rhoInt * uInt;
|
||||
fp = rhoInt * uInt * uInt + pInt;
|
||||
fe = (rhoInt * (pInt / ((_gamma - 1) * rhoInt) + 0.5 * uInt * uInt) + pInt) * uInt;
|
||||
return;
|
||||
}
|
||||
else if (uInt >= 0) // subsonic outflow
|
||||
{
|
||||
double s = pInt / Math.Pow(rhoInt, _gamma);
|
||||
double rhoGhost = Math.Pow(pAmb / s, 1.0 / _gamma);
|
||||
double cGhost = Math.Sqrt(_gamma * pAmb / rhoGhost);
|
||||
|
||||
// Outgoing Riemann invariant J⁺ = uInt + 2*cInt/(γ-1) (for right boundary)
|
||||
double J_plus = uInt + 2.0 * cInt / (_gamma - 1.0);
|
||||
double uGhost = J_plus - 2.0 * cGhost / (_gamma - 1.0);
|
||||
|
||||
// Clip to zero to prevent inflow
|
||||
if (uGhost > 0) uGhost = 0;
|
||||
|
||||
HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pAmb, out fm, out fp, out fe);
|
||||
}
|
||||
else // subsonic inflow
|
||||
{
|
||||
double T0 = 300.0;
|
||||
double R = 287.0;
|
||||
double rhoGhost = pAmb / (R * T0);
|
||||
HLLCFlux(rhoInt, uInt, pInt, rhoGhost, 0.0, pAmb, out fm, out fp, out fe);
|
||||
}
|
||||
}
|
||||
|
||||
// ========== Closed end (mirror) ==========
|
||||
private void ClosedEndFlux(double rhoInt, double uInt, double pInt, bool isRightBoundary,
|
||||
out double fm, out double fp, out double fe)
|
||||
{
|
||||
double rhoGhost = rhoInt;
|
||||
double pGhost = pInt;
|
||||
double uGhost = -uInt; // mirror velocity
|
||||
|
||||
if (isRightBoundary)
|
||||
HLLCFlux(rhoInt, uInt, pInt, rhoGhost, uGhost, pGhost, out fm, out fp, out fe);
|
||||
else
|
||||
HLLCFlux(rhoGhost, uGhost, pGhost, rhoInt, uInt, pInt, out fm, out fp, out fe);
|
||||
}
|
||||
|
||||
// ========== Standard HLLC flux ==========
|
||||
private void HLLCFlux(double rL, double uL, double pL,
|
||||
double rR, double uR, double pR,
|
||||
out double fm, out double fp, out double fe)
|
||||
|
||||
Reference in New Issue
Block a user