629 lines
30 KiB
C++
629 lines
30 KiB
C++
#pragma once
|
|
|
|
#include "FluidSim/Engine.hpp"
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <deque>
|
|
#include <limits>
|
|
#include <stdexcept>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
namespace fluidsim::port {
|
|
|
|
class Volume0D : public Component {
|
|
public:
|
|
explicit Volume0D(float volume_m3, float pressure = 101325.0f, float temperature = 300.0f,
|
|
float gas_constant_value = 287.0f, float gamma_value = 1.4f)
|
|
: gamma(gamma_value), gas_constant(gas_constant_value), volume_(volume_m3),
|
|
air_mass_(pressure * volume_m3 / (gas_constant_value * temperature)),
|
|
energy_(pressure * volume_m3 / (gamma_value - 1.0f)) {
|
|
if (!std::isfinite(volume_m3) || volume_m3 <= 0 || !std::isfinite(pressure) || pressure <= 0 ||
|
|
!std::isfinite(temperature) || temperature <= 0 || !std::isfinite(gas_constant_value) ||
|
|
gas_constant_value <= 0 || !std::isfinite(gamma_value) || gamma_value <= 1)
|
|
throw std::invalid_argument("volume thermodynamic state must be finite and positive, gamma > 1");
|
|
}
|
|
|
|
float gamma{1.4f};
|
|
float gas_constant{287.0f};
|
|
float ambient_pressure{101325.0f};
|
|
float ambient_temperature{300.0f};
|
|
float energy_relaxation_rate{0.0f};
|
|
float volume_change_rate{};
|
|
|
|
Port &create_port() {
|
|
ports_.emplace_back();
|
|
auto &p = ports_.back();
|
|
p.owner_energy = &energy_;
|
|
p.owner_air_mass = &air_mass_;
|
|
p.owner_exhaust_mass = &exhaust_mass_;
|
|
sync_port(p);
|
|
return p;
|
|
}
|
|
float pressure() const { return (gamma - 1.0f) * energy_ / std::max(volume_, 1e-12f); }
|
|
float density() const { return mass() / std::max(volume_, 1e-12f); }
|
|
float temperature() const { return pressure() / std::max(density() * gas_constant, 1e-9f); }
|
|
float volume() const { return volume_; }
|
|
float mass() const { return air_mass_ + exhaust_mass_; }
|
|
float internal_energy() const { return energy_; }
|
|
float air_fraction() const { return air_mass_ / std::max(mass(), 1e-12f); }
|
|
void set_pressure(float value, float temp = -1.0f) {
|
|
if (!std::isfinite(value) || value <= 0 || !std::isfinite(temp))
|
|
throw std::invalid_argument("invalid volume pressure/temperature");
|
|
const float target_t = temp > 0.0f ? temp : temperature();
|
|
const float total = value * volume_ / (gas_constant * target_t);
|
|
const float fraction = air_fraction();
|
|
air_mass_ = total * fraction;
|
|
exhaust_mass_ = total * (1.0f - fraction);
|
|
energy_ = value * volume_ / (gamma - 1.0f);
|
|
air_roundoff_ = exhaust_roundoff_ = energy_roundoff_ = 0;
|
|
sync_ports();
|
|
}
|
|
void update(float dt) {
|
|
if (!std::isfinite(dt) || dt < 0)
|
|
throw std::invalid_argument("invalid volume timestep");
|
|
if (dt == 0)
|
|
return;
|
|
double d_air = 0, d_exhaust = 0, d_energy = 0;
|
|
const float fraction = air_fraction();
|
|
for (const Port &p : ports_) {
|
|
const float flow_fraction = p.mass_flow >= 0.0f ? p.air_fraction : fraction;
|
|
d_air += double(p.mass_flow) * flow_fraction;
|
|
d_exhaust += double(p.mass_flow) * (1.0f - flow_fraction);
|
|
d_energy += double(p.mass_flow) * p.enthalpy;
|
|
}
|
|
accumulate(air_mass_, air_roundoff_, d_air * dt);
|
|
accumulate(exhaust_mass_, exhaust_roundoff_, d_exhaust * dt);
|
|
accumulate(energy_, energy_roundoff_, d_energy * dt);
|
|
const float new_volume = volume_ + volume_change_rate * dt;
|
|
if (!(new_volume > 0) || !std::isfinite(new_volume))
|
|
throw std::runtime_error("invalid moving volume");
|
|
if (new_volume != volume_) {
|
|
const double scale = std::pow(double(volume_) / new_volume, gamma - 1.0f);
|
|
const double precise = (double(energy_)+energy_roundoff_)*scale;
|
|
energy_ = float(precise);
|
|
energy_roundoff_ = precise-energy_;
|
|
}
|
|
volume_ = new_volume;
|
|
if (energy_relaxation_rate > 0.0f && mass() > 1e-12f) {
|
|
const float target = mass() * gas_constant * ambient_temperature / (gamma - 1.0f);
|
|
energy_ = target + (energy_ - target) * std::exp(-energy_relaxation_rate * dt);
|
|
energy_roundoff_ = 0;
|
|
}
|
|
if (!(mass() > 0) || !(energy_ > 0) || !std::isfinite(mass()) || !std::isfinite(energy_))
|
|
throw std::runtime_error(
|
|
"nonphysical volume state; reduce timestep or check external port fluxes");
|
|
if (air_mass_ < -1e-12f || exhaust_mass_ < -1e-12f)
|
|
throw std::runtime_error("negative volume species mass");
|
|
air_mass_ = std::max(air_mass_, 0.0f);
|
|
exhaust_mass_ = std::max(exhaust_mass_, 0.0f);
|
|
sync_ports();
|
|
}
|
|
void update_state(float dt) override { update(dt); }
|
|
std::vector<Port *> ports() override {
|
|
std::vector<Port *> result;
|
|
result.reserve(ports_.size());
|
|
for (Port &port : ports_)
|
|
result.push_back(&port);
|
|
return result;
|
|
}
|
|
const std::deque<Port> &ports() const { return ports_; }
|
|
|
|
private:
|
|
float volume_, air_mass_{}, exhaust_mass_{}, energy_;
|
|
// Preserve sub-ulp transfers without changing the float port/budget ABI.
|
|
double air_roundoff_{}, exhaust_roundoff_{}, energy_roundoff_{};
|
|
static void accumulate(float& value, double& roundoff, double increment) {
|
|
const double precise = double(value)+roundoff+increment;
|
|
value = float(precise);
|
|
roundoff = precise-value;
|
|
}
|
|
std::deque<Port> ports_;
|
|
void sync_port(Port &p) const {
|
|
p.gamma_value = gamma;
|
|
p.gas_constant_value = gas_constant;
|
|
p.pressure = pressure();
|
|
p.density = density();
|
|
p.temperature = temperature();
|
|
p.enthalpy = gamma / (gamma - 1.0f) * p.pressure / std::max(p.density, 1e-12f);
|
|
p.air_fraction = air_fraction();
|
|
}
|
|
void sync_ports() {
|
|
for (Port &p : ports_)
|
|
sync_port(p);
|
|
}
|
|
};
|
|
|
|
struct GasCell {
|
|
float rho{1.225f}, momentum{};
|
|
float energy{1.225f * (101325.0f / (pipe_gm1 * 1.225f))};
|
|
float air{1.0f};
|
|
};
|
|
struct GasState {
|
|
float rho{1.225f}, velocity{}, pressure{101325.0f}, air{1.0f};
|
|
};
|
|
|
|
class Pipe1D {
|
|
public:
|
|
// Integrated face fluxes: kg/s, N, W, kg/s of fresh air, positive along +x.
|
|
struct Flux {
|
|
float m{}, p{}, e{}, y{};
|
|
};
|
|
float area{}, dx{};
|
|
float damping_multiplier{0.0f}, energy_relaxation_rate{0.0f}, ambient_pressure{101325.0f};
|
|
std::vector<float> cell_areas, cell_lengths;
|
|
std::vector<GasCell> cells;
|
|
GasState left_ghost{}, right_ghost{};
|
|
std::uint64_t positivity_repairs{};
|
|
Pipe1D(float a, float length, int count)
|
|
: Pipe1D(uniform(a, count), uniform(length / std::max(count, 1), count)) {}
|
|
Pipe1D(std::vector<float> areas, std::vector<float> lengths)
|
|
: area(areas.empty() ? 0 : areas.front()), dx(lengths.empty() ? 0 : lengths.front()),
|
|
cell_areas(std::move(areas)), cell_lengths(std::move(lengths)), cells(cell_areas.size()),
|
|
next_cells_(cells.size()), states_(cells.size()), waves_(cells.size()), faces_(cells.size() + 1) {
|
|
if (cells.empty() || cell_lengths.size() != cells.size())
|
|
throw std::invalid_argument("pipe profiles must be nonempty and have equal length");
|
|
for (std::size_t i = 0; i < cells.size(); ++i)
|
|
if (!std::isfinite(cell_areas[i]) || !std::isfinite(cell_lengths[i]) || cell_areas[i] <= 0 ||
|
|
cell_lengths[i] <= 0)
|
|
throw std::invalid_argument("pipe geometry must be finite and positive");
|
|
}
|
|
GasState state(int i) const { return primitive(cells.at(static_cast<std::size_t>(i))); }
|
|
float cell_mass(int i) const { return cells.at(i).rho * cell_areas.at(i) * cell_lengths.at(i); }
|
|
float cell_area(int i) const { return cell_areas.at(i); }
|
|
float cell_length(int i) const { return cell_lengths.at(i); }
|
|
float minimum_cell_length() const { return *std::min_element(cell_lengths.begin(), cell_lengths.end()); }
|
|
float stable_timestep(float cfl = 0.8f) const {
|
|
float dt = std::numeric_limits<float>::max();
|
|
for (std::size_t i = 0; i < cells.size(); ++i) {
|
|
auto s = primitive(cells[i]);
|
|
if (!(s.rho > 0) || !(s.pressure > 0) || !std::isfinite(s.velocity) || !std::isfinite(s.pressure))
|
|
throw std::runtime_error("nonphysical pipe state in CFL calculation");
|
|
dt = std::min(dt, cfl * cell_lengths[i] / wave(s));
|
|
}
|
|
return dt;
|
|
}
|
|
// Solver preparation is valid until advance; boundary resolution changes
|
|
// fluxes/ghosts, not cells. Standalone advance still prepares its own cache.
|
|
float prepare_step(float cfl = 0.8f) {
|
|
float dt = std::numeric_limits<float>::max();
|
|
for (std::size_t i = 0; i < cells.size(); ++i) {
|
|
states_[i] = primitive(cells[i]);
|
|
const auto &s = states_[i];
|
|
if (!(s.rho > 0) || !(s.pressure > 0) || !std::isfinite(s.velocity) || !std::isfinite(s.pressure))
|
|
throw std::runtime_error("nonphysical pipe state in preparation");
|
|
waves_[i] = wave(s);
|
|
dt = std::min(dt, cfl * cell_lengths[i] / waves_[i]);
|
|
}
|
|
prepared_ = true;
|
|
return dt;
|
|
}
|
|
int required_substeps(float dt, float cfl = 0.8f) const {
|
|
if (!std::isfinite(dt) || dt < 0 || !std::isfinite(cfl) || cfl <= 0 || cfl > 1)
|
|
throw std::invalid_argument("invalid timestep/CFL");
|
|
const double n = std::ceil(double(dt) / stable_timestep(cfl));
|
|
if (n > 1000000)
|
|
throw std::runtime_error("excessive pipe substeps");
|
|
return std::max(1, static_cast<int>(n));
|
|
}
|
|
void set_boundary_flux(bool left, Flux f) {
|
|
if (left) {
|
|
left_flux_ = f;
|
|
left_override_ = true;
|
|
} else {
|
|
right_flux_ = f;
|
|
right_override_ = true;
|
|
}
|
|
}
|
|
static GasState primitive(const GasCell &c) {
|
|
const float u = c.momentum / std::max(c.rho, 1e-12f);
|
|
return {c.rho, u, pipe_gm1 * (c.energy - .5f * c.momentum * u), std::clamp(c.air, 0.0f, 1.0f)};
|
|
}
|
|
static float wave(const GasState &s) {
|
|
return std::abs(s.velocity) + std::sqrt(gamma * s.pressure / s.rho);
|
|
}
|
|
static Flux riemann(const GasState &a, const GasState &b) {
|
|
return riemann(a, b, std::max(wave(a), wave(b)));
|
|
}
|
|
void advance(float dt) {
|
|
if (dt <= 0)
|
|
return;
|
|
if (!std::isfinite(dt))
|
|
throw std::invalid_argument("nonfinite pipe timestep");
|
|
const std::size_t n = cells.size();
|
|
if (!prepared_)
|
|
prepare_step();
|
|
faces_[0] = left_override_ ? left_flux_ : scaled(riemann(left_ghost, states_[0]), cell_areas[0]);
|
|
faces_[n] =
|
|
right_override_ ? right_flux_ : scaled(riemann(states_[n - 1], right_ghost), cell_areas[n - 1]);
|
|
for (std::size_t i = 1; i < n; ++i)
|
|
faces_[i] = scaled(riemann(states_[i - 1], states_[i], std::max(waves_[i - 1], waves_[i])),
|
|
.5f * (cell_areas[i - 1] + cell_areas[i]));
|
|
const float relaxation = energy_relaxation_rate > 0 ? std::exp(-energy_relaxation_rate * dt) : 1.0f;
|
|
for (std::size_t i = 0; i < n; ++i) {
|
|
const auto fl = faces_[i], fr = faces_[i + 1];
|
|
const float factor = dt / (cell_areas[i] * cell_lengths[i]);
|
|
const float al = i ? .5f * (cell_areas[i - 1] + cell_areas[i]) : cell_areas[i];
|
|
const float ar = i + 1 < n ? .5f * (cell_areas[i] + cell_areas[i + 1]) : cell_areas[i];
|
|
float rho = cells[i].rho - factor * (fr.m - fl.m);
|
|
float mom = cells[i].momentum - factor * (fr.p - fl.p - states_[i].pressure * (ar - al));
|
|
float en = cells[i].energy - factor * (fr.e - fl.e);
|
|
const float air = cells[i].rho * cells[i].air - factor * (fr.y - fl.y);
|
|
if (!std::isfinite(rho) || !std::isfinite(mom) || !std::isfinite(en) || !std::isfinite(air))
|
|
throw std::runtime_error("nonfinite finite-volume update");
|
|
if (rho < 1e-12f) {
|
|
rho = 1e-12f;
|
|
++positivity_repairs;
|
|
}
|
|
if (damping_multiplier > 0)
|
|
mom *= std::exp(-.1f * damping_multiplier * dt / cells[i].rho);
|
|
if (energy_relaxation_rate > 0)
|
|
en = ambient_pressure * pipe_gm1_inverse +
|
|
(en - ambient_pressure * pipe_gm1_inverse) * relaxation;
|
|
const float minimum = 1e-6f * pipe_gm1_inverse + .5f * mom * mom / rho;
|
|
if (en < minimum) {
|
|
en = minimum;
|
|
++positivity_repairs;
|
|
}
|
|
next_cells_[i] = {rho, mom, en, std::clamp(air / rho, 0.0f, 1.0f)};
|
|
}
|
|
cells.swap(next_cells_);
|
|
prepared_ = false;
|
|
left_override_ = right_override_ = false;
|
|
}
|
|
static Flux scaled(Flux f, float a) { return {f.m * a, f.p * a, f.e * a, f.y * a}; }
|
|
|
|
private:
|
|
std::vector<GasCell> next_cells_;
|
|
std::vector<GasState> states_;
|
|
std::vector<float> waves_;
|
|
std::vector<Flux> faces_;
|
|
Flux left_flux_{}, right_flux_{};
|
|
bool left_override_{}, right_override_{}, prepared_{};
|
|
static std::vector<float> uniform(float v, int n) {
|
|
if (n <= 0)
|
|
throw std::invalid_argument("pipe count must be positive");
|
|
return std::vector<float>(n, v);
|
|
}
|
|
static Flux riemann(const GasState &a, const GasState &b, float c) {
|
|
const float ma = a.rho * a.velocity, mb = b.rho * b.velocity;
|
|
const float ea = a.pressure * pipe_gm1_inverse + .5f * ma * a.velocity;
|
|
const float eb = b.pressure * pipe_gm1_inverse + .5f * mb * b.velocity;
|
|
return {.5f * (ma + mb - c * (b.rho - a.rho)),
|
|
.5f * (ma * a.velocity + a.pressure + mb * b.velocity + b.pressure - c * (mb - ma)),
|
|
.5f * ((ea + a.pressure) * a.velocity + (eb + b.pressure) * b.velocity - c * (eb - ea)),
|
|
.5f * (ma * a.air + mb * b.air - c * (b.rho * b.air - a.rho * a.air))};
|
|
}
|
|
};
|
|
|
|
class BoundarySystem {
|
|
public:
|
|
explicit BoundarySystem(std::vector<Pipe1D> &pipes) : pipes_(pipes) {}
|
|
void add_orifice(Port &port, int pipe, bool left, float &area, float cd = 1.0f) {
|
|
if (std::abs(port.gamma_value - gamma) > 1e-6f ||
|
|
std::abs(port.gas_constant_value - gas_constant) > .001f)
|
|
throw std::invalid_argument("pipe coupling requires the same calorically perfect gas");
|
|
if (!std::isfinite(cd) || cd < 0 || cd > 1)
|
|
throw std::invalid_argument("orifice Cd must be in [0,1]");
|
|
for (auto &o : orifices_)
|
|
if (o.port == &port)
|
|
throw std::invalid_argument("each connection requires a distinct component port");
|
|
claim(pipe, left);
|
|
Orifice o;
|
|
o.port = &port;
|
|
o.pipe = pipe;
|
|
o.left = left;
|
|
o.area = &area;
|
|
o.cd = cd;
|
|
if (port.owner_air_mass) {
|
|
auto it = std::find_if(budgets_.begin(), budgets_.end(),
|
|
[&](const Budget &b) { return b.air == port.owner_air_mass; });
|
|
if (it == budgets_.end()) {
|
|
o.budget = budgets_.size();
|
|
budgets_.push_back({port.owner_air_mass, port.owner_exhaust_mass, port.owner_energy});
|
|
} else
|
|
o.budget = static_cast<std::size_t>(it - budgets_.begin());
|
|
}
|
|
orifices_.push_back(o);
|
|
}
|
|
void add_orifice_with_inertance(Port &p, int pipe, bool left, float &a, float cd, float length,
|
|
float loss = 0) {
|
|
if (!std::isfinite(length) || length <= 0 || !std::isfinite(loss) || loss < 0)
|
|
throw std::invalid_argument("invalid inertance");
|
|
add_orifice(p, pipe, left, a, cd);
|
|
auto &o = orifices_.back();
|
|
o.length = length;
|
|
o.loss = loss;
|
|
}
|
|
void add_closed_end(int pipe, bool left) {
|
|
claim(pipe, left);
|
|
Orifice o;
|
|
o.pipe = pipe;
|
|
o.left = left;
|
|
orifices_.push_back(o);
|
|
}
|
|
void add_open_end(int pipe, bool left, float p, float a) {
|
|
if (!std::isfinite(p) || p <= 0 || !std::isfinite(a) || a <= 0)
|
|
throw std::invalid_argument("invalid open boundary");
|
|
claim(pipe, left);
|
|
opens_.push_back({pipe, left, p, a});
|
|
}
|
|
// An ideal, zero-volume junction between two duct faces. Orientation is explicit.
|
|
void add_connection(int a, bool left_a, int b, bool left_b) {
|
|
if (a < 0 || b < 0 || static_cast<std::size_t>(a) >= pipes_.size() ||
|
|
static_cast<std::size_t>(b) >= pipes_.size())
|
|
throw std::invalid_argument("invalid connection pipe index");
|
|
const float aa = pipes_[a].cell_area(end(pipes_[a], left_a));
|
|
const float ab = pipes_[b].cell_area(end(pipes_[b], left_b));
|
|
if (std::abs(aa - ab) > 1e-6f * std::max(aa, ab))
|
|
throw std::invalid_argument("unequal-area junction needs a profiled pipe transition");
|
|
if (a == b && left_a == left_b)
|
|
throw std::invalid_argument("cannot join a face to itself");
|
|
check_face(a, left_a);
|
|
check_face(b, left_b);
|
|
claim(a, left_a);
|
|
claim(b, left_b);
|
|
connections_.push_back({a, b, left_a, left_b});
|
|
}
|
|
float open_end_mass_flow(std::size_t i) const { return opens_.at(i).last_mass_flow; }
|
|
float open_end_pressure(std::size_t i) const { return opens_.at(i).last_face_pressure; }
|
|
std::uint64_t limited_fluxes{}, nonlinear_iterations{};
|
|
void validate_topology() const {
|
|
if (claims_.size() != 2 * pipes_.size())
|
|
throw std::runtime_error("every pipe endpoint needs an explicit boundary or connection");
|
|
}
|
|
void resolve(float dt) {
|
|
if (!std::isfinite(dt) || dt <= 0)
|
|
throw std::invalid_argument("boundary timestep must be positive");
|
|
for (auto &b : budgets_) {
|
|
b.mass = .25f * (*b.air + (b.exhaust ? *b.exhaust : 0));
|
|
b.heat = b.energy ? .25f * (*b.energy) : std::numeric_limits<float>::max();
|
|
}
|
|
for (auto &o : orifices_)
|
|
resolve_orifice(o, dt);
|
|
for (auto &o : opens_)
|
|
resolve_open(o);
|
|
for (auto &o : connections_)
|
|
resolve_connection(o);
|
|
}
|
|
|
|
private:
|
|
struct Budget {
|
|
float *air{}, *exhaust{}, *energy{};
|
|
float mass{}, heat{};
|
|
};
|
|
struct Orifice {
|
|
Port *port{};
|
|
int pipe{};
|
|
bool left{};
|
|
float *area{};
|
|
float cd{1}, length{}, loss{}, mdot{};
|
|
std::size_t budget{std::numeric_limits<std::size_t>::max()};
|
|
};
|
|
struct Open {
|
|
int pipe;
|
|
bool left;
|
|
float pressure, area, last_mass_flow{}, last_face_pressure{};
|
|
};
|
|
struct Connection {
|
|
int a, b;
|
|
bool left_a, left_b;
|
|
};
|
|
struct Claim {
|
|
int pipe;
|
|
bool left;
|
|
};
|
|
std::vector<Pipe1D> &pipes_;
|
|
std::vector<Orifice> orifices_;
|
|
std::vector<Open> opens_;
|
|
std::vector<Connection> connections_;
|
|
std::vector<Budget> budgets_;
|
|
std::vector<Claim> claims_;
|
|
void check_face(int pipe, bool left) const {
|
|
if (pipe < 0 || static_cast<std::size_t>(pipe) >= pipes_.size())
|
|
throw std::invalid_argument("invalid boundary pipe index");
|
|
for (auto c : claims_)
|
|
if (c.pipe == pipe && c.left == left)
|
|
throw std::invalid_argument("pipe face already connected");
|
|
}
|
|
void claim(int pipe, bool left) {
|
|
check_face(pipe, left);
|
|
claims_.push_back({pipe, left});
|
|
}
|
|
static int end(const Pipe1D &p, bool left) { return left ? 0 : static_cast<int>(p.cells.size()) - 1; }
|
|
static void ghost(Pipe1D &p, bool left, GasState s) {
|
|
if (left)
|
|
p.left_ghost = s;
|
|
else
|
|
p.right_ghost = s;
|
|
}
|
|
void resolve_orifice(Orifice &o, float dt) {
|
|
auto &p = pipes_[o.pipe];
|
|
int i = end(p, o.left);
|
|
auto s = p.state(i);
|
|
float face_area = p.cell_area(i);
|
|
if (o.area && (!std::isfinite(*o.area) || *o.area < 0))
|
|
throw std::invalid_argument("invalid dynamic orifice area");
|
|
if (o.port && (!(o.port->pressure > 0) || !(o.port->density > 0) ||
|
|
!std::isfinite(o.port->pressure) || !std::isfinite(o.port->density)))
|
|
throw std::runtime_error("nonphysical component port");
|
|
const float area = o.area ? std::min(std::max(*o.area, 0.f), face_area) : 0.f;
|
|
if (!o.port || area <= 0 || o.cd == 0) {
|
|
auto g = s;
|
|
g.velocity = -s.velocity;
|
|
ghost(p, o.left, g);
|
|
auto wall_flux =
|
|
Pipe1D::scaled(o.left ? Pipe1D::riemann(g, s) : Pipe1D::riemann(s, g), face_area);
|
|
// A stationary adiabatic wall transfers momentum, never mass or energy.
|
|
p.set_boundary_flux(o.left, {0, wall_flux.p, 0, 0});
|
|
if (o.port)
|
|
o.port->mass_flow = 0;
|
|
o.mdot = 0;
|
|
return;
|
|
}
|
|
auto &v = *o.port;
|
|
if (std::abs(v.gamma_value - gamma) > 1e-6f || std::abs(v.gas_constant_value - gas_constant) > .001f)
|
|
throw std::invalid_argument("pipe coupling requires the same calorically perfect gas");
|
|
const float dp = v.pressure - s.pressure;
|
|
const double g = gamma;
|
|
const double critical = std::pow(2.0 / (g + 1), g / (g - 1));
|
|
const double c = std::sqrt(g * s.pressure / s.rho);
|
|
const double ui = s.velocity * (o.left ? 1 : -1);
|
|
const auto face_density = [&](double pressure) {
|
|
return s.rho * std::pow(pressure / s.pressure, 1 / g);
|
|
};
|
|
const auto nozzle_flow = [&](double pressure) {
|
|
const bool from = v.pressure >= pressure;
|
|
const double up = from ? v.pressure : pressure;
|
|
const double down = from ? pressure : v.pressure;
|
|
const double density = from ? v.density : face_density(pressure);
|
|
const double ratio = std::max(down / up, critical);
|
|
const double r = std::pow(ratio, 1 / g);
|
|
const double speed2 = 2 * g / (g - 1) * up / density * (1 - std::pow(ratio, (g - 1) / g));
|
|
return (from ? 1 : -1) * density * r * std::sqrt(std::max(0.0, speed2)) * area * o.cd;
|
|
};
|
|
// Match nozzle discharge to the outgoing duct characteristic. Imposing
|
|
// sqrt(dp) against the *cell-centre* pressure is an unstable explicit
|
|
// boundary near equilibrium and ignores the duct's acoustic impedance.
|
|
const double zero_sound = c - .5 * (g - 1) * ui;
|
|
if (zero_sound <= 0)
|
|
throw std::runtime_error("orifice characteristic outside subsonic model");
|
|
const double zero_pressure = s.pressure * std::pow(zero_sound / c, 2 * g / (g - 1));
|
|
double low = std::min(zero_pressure, double(v.pressure));
|
|
double high = std::max(zero_pressure, double(v.pressure));
|
|
double duct_pressure = .5 * (low + high);
|
|
const double flow_sign = v.pressure >= zero_pressure ? 1 : -1;
|
|
const double flow_scale = area * o.cd;
|
|
const double density_ratio_scale = std::pow(double(s.pressure) / v.pressure, 1 / g);
|
|
const double critical_density_ratio = std::pow(critical, 1 / g);
|
|
// Squared-flow residual avoids the nozzle sqrt singularity at dp=0.
|
|
// Safeguarded analytic Newton; bisection remains the fallback.
|
|
if (o.length == 0 && high > low)
|
|
for (int iteration = 0; iteration < 40; ++iteration) {
|
|
++nonlinear_iterations;
|
|
const double density = face_density(duct_pressure);
|
|
const double sound = std::sqrt(g * duct_pressure / density);
|
|
const double velocity = ui + 2 * (sound - c) / (g - 1);
|
|
const double q = density * p.cell_area(i) * velocity;
|
|
const double dq = p.cell_area(i) * density / (g * duct_pressure) * (velocity + sound);
|
|
const bool from = v.pressure >= duct_pressure;
|
|
const double up = from ? v.pressure : duct_pressure;
|
|
const double down = from ? duct_pressure : v.pressure;
|
|
const double upstream_density = from ? v.density : density;
|
|
const double raw_ratio = down / up;
|
|
const double ratio = std::max(raw_ratio, critical);
|
|
const double root_ratio = raw_ratio <= critical
|
|
? critical_density_ratio
|
|
: (from ? density / s.rho * density_ratio_scale
|
|
: s.rho / (density * density_ratio_scale));
|
|
const double r1 = root_ratio * root_ratio, r2 = root_ratio * ratio;
|
|
const double coefficient = flow_scale * flow_scale * 2 * g / (g - 1) * upstream_density * up;
|
|
const double q2 = coefficient * (r1 - r2);
|
|
const double dr = raw_ratio > critical ? (from ? 1 / up : -ratio / up) : 0;
|
|
const double derivative = coefficient * ((2 / g * r1 - (g + 1) / g * r2) / ratio) * dr +
|
|
(from ? 0 : q2 * (1 + 1 / g) / duct_pressure);
|
|
const double residual = flow_sign * (q * q - q2);
|
|
if (residual > 0)
|
|
high = duct_pressure;
|
|
else
|
|
low = duct_pressure;
|
|
if (high - low < std::max(1e-9, duct_pressure * 1e-12))
|
|
break;
|
|
const double slope = flow_sign * (2 * q * dq - derivative);
|
|
const double candidate = duct_pressure - residual / slope;
|
|
if (std::abs(residual) < 4e-7 * std::max(q * q + q2, 1e-24))
|
|
break;
|
|
duct_pressure = (std::isfinite(candidate) && candidate > low && candidate < high)
|
|
? candidate
|
|
: .5 * (low + high);
|
|
}
|
|
float mdot = static_cast<float>(nozzle_flow(duct_pressure));
|
|
const bool upstream_port = mdot >= 0;
|
|
const float pup = upstream_port ? v.pressure : static_cast<float>(duct_pressure);
|
|
const float ratio =
|
|
static_cast<float>(std::max((upstream_port ? duct_pressure : v.pressure) / pup, critical));
|
|
const float throat_rho =
|
|
(upstream_port ? v.density : static_cast<float>(face_density(duct_pressure))) *
|
|
std::pow(ratio, 1 / gamma);
|
|
if (o.length > 0) {
|
|
// d(mdot)/dt = A/L dp - K mdot |mdot|/(2 rho A L).
|
|
const float donor_rho = o.mdot >= 0 ? v.density : s.rho;
|
|
// Semi-implicit resistance remains dissipative for large dt.
|
|
mdot = (o.mdot + area / o.length * dp * dt) /
|
|
(1 + dt * o.loss * std::abs(o.mdot) / (2 * donor_rho * area * o.length));
|
|
}
|
|
const bool from = mdot >= 0;
|
|
const float h = from ? gamma / (gamma - 1) * v.pressure / v.density
|
|
: gamma / (gamma - 1) * s.pressure / s.rho + .5f * s.velocity * s.velocity;
|
|
const float requested_mdot = mdot;
|
|
if (from && o.budget < budgets_.size()) {
|
|
auto &b = budgets_[o.budget];
|
|
mdot = std::min(mdot, std::max(0.f, std::min(b.mass, b.heat / std::max(h, 1.f))) / dt);
|
|
b.mass -= mdot * dt;
|
|
b.heat -= mdot * h * dt;
|
|
} else if (!from) {
|
|
const float cell_energy = p.cells[i].energy * face_area * p.cell_length(i);
|
|
mdot = std::max(mdot, -.25f * std::min(p.cell_mass(i), cell_energy / std::max(h, 1.f)) / dt);
|
|
}
|
|
if (mdot != requested_mdot)
|
|
++limited_fluxes;
|
|
o.mdot = mdot;
|
|
const float donor_rho = from ? v.density : s.rho;
|
|
const float density = o.length > 0 ? donor_rho : throat_rho;
|
|
const float velocity = mdot / (density * area);
|
|
const float face_pressure =
|
|
o.length > 0 ? .5f * (v.pressure + s.pressure) : static_cast<float>(duct_pressure);
|
|
const float direction = o.left ? 1.f : -1.f;
|
|
const float species = from ? v.air_fraction : s.air;
|
|
const float signed_mass = direction * mdot;
|
|
p.set_boundary_flux(o.left, {signed_mass,
|
|
mdot * velocity + face_pressure * area + s.pressure * (face_area - area),
|
|
signed_mass * h, signed_mass * species});
|
|
ghost(p, o.left, {density, direction * velocity, face_pressure, species});
|
|
v.mass_flow = -mdot;
|
|
v.enthalpy = h;
|
|
if (!from)
|
|
v.air_fraction = s.air;
|
|
}
|
|
void resolve_open(Open &o) {
|
|
auto &p = pipes_[o.pipe];
|
|
int i = end(p, o.left);
|
|
auto s = p.state(i);
|
|
const float c = std::sqrt(gamma * s.pressure / s.rho), gm1 = gamma - 1;
|
|
if (o.left ? s.velocity >= c : s.velocity <= -c)
|
|
throw std::runtime_error(
|
|
"supersonic inlet needs full upstream state; pressure-only boundary is insufficient");
|
|
const float entropy = s.pressure / std::pow(s.rho, gamma);
|
|
const float rho = std::pow(o.pressure / entropy, 1 / gamma);
|
|
const float cb = std::sqrt(gamma * o.pressure / rho);
|
|
const float u =
|
|
o.left ? s.velocity - 2 * c / gm1 + 2 * cb / gm1 : s.velocity + 2 * c / gm1 - 2 * cb / gm1;
|
|
const bool out = o.left ? s.velocity <= -c : s.velocity >= c;
|
|
GasState g = out ? s : GasState{rho, u, o.pressure, (o.left ? u >= 0 : u <= 0) ? 1.f : s.air};
|
|
ghost(p, o.left, g);
|
|
const float a = std::min(o.area, p.cell_area(i));
|
|
auto f = Pipe1D::scaled(o.left ? Pipe1D::riemann(g, s) : Pipe1D::riemann(s, g), a);
|
|
f.p += s.pressure * (p.cell_area(i) - a);
|
|
p.set_boundary_flux(o.left, f);
|
|
o.last_mass_flow = f.m * (o.left ? -1 : 1);
|
|
o.last_face_pressure = f.p / p.cell_area(i);
|
|
}
|
|
void resolve_connection(const Connection &o) {
|
|
auto &a = pipes_[o.a];
|
|
auto &b = pipes_[o.b];
|
|
int ia = end(a, o.left_a), ib = end(b, o.left_b);
|
|
auto sa = a.state(ia), sb = b.state(ib);
|
|
const float da = o.left_a ? -1.f : 1.f, db = o.left_b ? 1.f : -1.f;
|
|
sa.velocity *= da;
|
|
sb.velocity *= db;
|
|
const float area = std::min(a.cell_area(ia), b.cell_area(ib));
|
|
auto f = Pipe1D::scaled(Pipe1D::riemann(sa, sb), area);
|
|
a.set_boundary_flux(o.left_a,
|
|
{da * f.m, f.p + sa.pressure * (a.cell_area(ia) - area), da * f.e, da * f.y});
|
|
b.set_boundary_flux(o.left_b,
|
|
{db * f.m, f.p + sb.pressure * (b.cell_area(ib) - area), db * f.e, db * f.y});
|
|
}
|
|
};
|
|
|
|
} // namespace fluidsim::port
|