847 lines
36 KiB
C++
847 lines
36 KiB
C++
#include "FluidSim/Atmosphere.hpp"
|
|
#include "FluidSim/Crankcase.hpp"
|
|
#include "FluidSim/Fluid.hpp"
|
|
#include "FluidSim/Junction.hpp"
|
|
#include "FluidSim/Solver.hpp"
|
|
#include <cmath>
|
|
#include <complex>
|
|
#include <functional>
|
|
#include <iostream>
|
|
#include <string>
|
|
using namespace fluidsim::port;
|
|
namespace {
|
|
PipeSolver selected_solver=PipeSolver::Rusanov;
|
|
void check(bool b, const std::string &s) {
|
|
if (!b)
|
|
throw std::runtime_error(s);
|
|
}
|
|
void near(double a, double b, double tol, const char *s) {
|
|
check(std::isfinite(a) && std::abs(a - b) <= tol, s);
|
|
}
|
|
void set(Pipe1D &p, int i, float pressure, float rho, float u = 0, float air = 1, float fuel = 0) {
|
|
p.cells[i] = {rho, rho * u, pressure / pipe_gm1 + .5f * rho * u * u, air, fuel};
|
|
}
|
|
void lossless(Pipe1D &p) {
|
|
p.solver=selected_solver;
|
|
p.damping_multiplier = 0;
|
|
p.energy_relaxation_rate = 0;
|
|
}
|
|
double mass(const Pipe1D &p) {
|
|
double v = 0;
|
|
for (int i = 0; i < int(p.cells.size()); i++)
|
|
v += p.cell_mass(i);
|
|
return v;
|
|
}
|
|
double energy(const Pipe1D &p) {
|
|
double v = 0;
|
|
for (int i = 0; i < int(p.cells.size()); i++)
|
|
v += double(p.cells[i].energy) * p.cell_area(i) * p.cell_length(i);
|
|
return v;
|
|
}
|
|
double vol_energy(const Volume0D &v) { return v.pressure() * v.volume() / (v.gamma - 1); }
|
|
} // namespace
|
|
int main(int argc,char** argv) {
|
|
if(argc>1) selected_solver=parse_pipe_solver(argv[1]);
|
|
int failed = 0, total = 0;
|
|
auto test = [&](const char *n, std::function<void()> fn) {
|
|
++total;
|
|
try {
|
|
fn();
|
|
std::cout << "PASS " << n << '\n';
|
|
} catch (const std::exception &e) {
|
|
++failed;
|
|
std::cout << "FAIL " << n << ": " << e.what() << '\n';
|
|
}
|
|
};
|
|
test("solver flux consistency and reflection symmetry", [] {
|
|
for(auto solver : {PipeSolver::Rusanov,PipeSolver::Hlle,PipeSolver::Hllc}) {
|
|
GasState a{1.2f,30,101325,.8f,.05f}, b{.5f,-10,60000,.1f,.3f};
|
|
const auto equal=Pipe1D::numerical_flux(a,a,solver);
|
|
near(equal.m,a.rho*a.velocity,1e-5,"uniform mass flux");
|
|
near(equal.p,a.rho*a.velocity*a.velocity+a.pressure,.02,"uniform momentum flux");
|
|
const auto forward=Pipe1D::numerical_flux(a,b,solver);
|
|
a.velocity=-a.velocity; b.velocity=-b.velocity;
|
|
const auto reverse=Pipe1D::numerical_flux(b,a,solver);
|
|
near(forward.m,-reverse.m,.0001,"orientation mass symmetry");
|
|
near(forward.p,reverse.p,.1,"orientation momentum symmetry");
|
|
near(forward.e,-reverse.e,1,"orientation energy symmetry");
|
|
near(forward.y,-reverse.y,.0001,"orientation air symmetry");
|
|
near(forward.fuel,-reverse.fuel,.0001,"orientation fuel symmetry");
|
|
}
|
|
});
|
|
test("HLLC preserves stationary contact and falls back at expansion vacuum", [] {
|
|
GasState a{1,0,1,.8f,.1f}, b{.125f,0,1,.1f,.5f};
|
|
auto f=Pipe1D::numerical_flux(a,b,PipeSolver::Hllc);
|
|
near(f.m,0,1e-7,"stationary contact mass diffusion");
|
|
near(f.e,0,1e-7,"stationary contact energy diffusion");
|
|
near(f.y,0,1e-7,"stationary contact tracer diffusion");
|
|
a.velocity=-10; b.velocity=10;
|
|
bool fallback=false;
|
|
f=Pipe1D::numerical_flux(a,b,PipeSolver::Hllc,std::sqrt(1.4f),std::sqrt(1.4f/.125f),&fallback);
|
|
check(fallback && std::isfinite(f.e),"strong expansion did not use safe HLLE flux");
|
|
});
|
|
test("unflanged radiation reflection magnitude and phase", [] {
|
|
for (double ka : {.1, .5, 1., 2.}) {
|
|
OpenEndRadiation reflection;
|
|
std::complex<double> input{}, output{};
|
|
double t = 0;
|
|
const double dt = 2*3.141592653589793/(ka*4096);
|
|
for (int i=0; i<4096*32; ++i) {
|
|
t += dt;
|
|
const double incident = std::cos(ka*t);
|
|
const double reflected = reflection.process(incident, dt);
|
|
if (i >= 4096*16) {
|
|
auto phase = std::exp(std::complex<double>(0, -ka*t));
|
|
input += incident*phase;
|
|
output += reflected*phase;
|
|
}
|
|
}
|
|
constexpr double n=.167, d=n+2*.6133, e=(d*d-n*n-1)/2;
|
|
const std::complex<double> z(0,ka);
|
|
const auto expected = -(1.+n*z)/(1.+d*z+e*z*z);
|
|
near(std::abs(output/input-expected), 0, .001, "radiation amplitude/phase mismatch");
|
|
check(std::abs(output/input) < 1, "radiation boundary creates acoustic energy");
|
|
}
|
|
});
|
|
test("N-port junction conserves sealed network mass and energy", [] {
|
|
std::vector<Pipe1D> pipes;
|
|
for(int i=0;i<3;++i) { pipes.emplace_back(.001f*(i+1),.3f,10);lossless(pipes.back()); }
|
|
for(int i=0;i<10;++i) set(pipes[0],i,120000,1.4f);
|
|
BoundarySystem boundaries(pipes);
|
|
Junction junction(.0003f);
|
|
const std::array<JunctionEnd,3> ends{{{0,false},{1,true},{2,true}}};
|
|
junction.connect(boundaries,pipes,ends);
|
|
boundaries.add_closed_end(0,true);boundaries.add_closed_end(1,false);boundaries.add_closed_end(2,false);
|
|
Solver solver(pipes,boundaries);solver.add_component(junction);
|
|
auto totals=[&] { std::array<double,2> t{junction.mass(),junction.internal_energy()};
|
|
for(auto& p:pipes) { t[0]+=mass(p);t[1]+=energy(p); } return t; };
|
|
const auto before=totals();
|
|
for(int i=0;i<1000;++i) solver.step(1.f/44100);
|
|
const auto after=totals();
|
|
if(std::abs(after[0]-before[0])>before[0]*2e-6) std::cerr << "mass relative drift=" << (after[0]-before[0])/before[0] << "\n";
|
|
near(after[0],before[0],before[0]*2e-6,"junction mass drift");
|
|
near(after[1],before[1],before[1]*2e-6,"junction energy drift");
|
|
check(pipes[1].state(0).pressure>101325,"junction did not transmit pressure to branch");
|
|
for(auto& p:pipes) check(p.positivity_repairs==0,"junction required positivity repairs");
|
|
});
|
|
test("Zero-storage N-port junction conserves sealed network mass and energy", [] {
|
|
std::vector<Pipe1D> pipes;
|
|
for(int i=0;i<3;++i) { pipes.emplace_back(.001f*(i+1),.3f,10);lossless(pipes.back()); }
|
|
for(int i=0;i<10;++i) set(pipes[0],i,120000,1.4f);
|
|
BoundarySystem boundaries(pipes);
|
|
boundaries.add_junction({{0,false},{1,true},{2,true}});
|
|
boundaries.add_closed_end(0,true);boundaries.add_closed_end(1,false);boundaries.add_closed_end(2,false);
|
|
Solver solver(pipes,boundaries);
|
|
auto totals=[&] { std::array<double,2> t{0,0};
|
|
for(auto& p:pipes) { t[0]+=mass(p);t[1]+=energy(p); } return t; };
|
|
const auto before=totals();
|
|
for(int i=0;i<1000;++i) solver.step(1.f/44100);
|
|
const auto after=totals();
|
|
// Float cell updates accumulate roundoff over 4000 substeps (50 ppm bound).
|
|
near(after[0],before[0],before[0]*5e-5,"junction mass drift");
|
|
near(after[1],before[1],before[1]*5e-5,"junction energy drift");
|
|
check(pipes[1].state(0).pressure>101325,"junction did not transmit pressure to branch");
|
|
for(auto& p:pipes) check(p.positivity_repairs==0,"junction required positivity repairs");
|
|
});
|
|
test("radiating acoustic volume source is independent of donor density", [] {
|
|
for (bool left : {false,true}) {
|
|
std::vector<Pipe1D> hot, cold;
|
|
hot.emplace_back(.001f,.1f,1); cold.emplace_back(.001f,.1f,1);
|
|
BoundarySystem a(hot), b(cold);
|
|
a.add_radiating_open_end(0,left,101325,.001f);
|
|
b.add_radiating_open_end(0,left,101325,.001f);
|
|
for(int i=0;i<4096;++i) {
|
|
const float velocity=2*std::sin(2*pi*150*i/44100);
|
|
set(hot[0],0,101325,.5f,velocity);
|
|
set(cold[0],0,101325,ambient_density,velocity);
|
|
a.resolve(1.f/44100); b.resolve(1.f/44100);
|
|
near(a.open_end_volume_flow(0),b.open_end_volume_flow(0),2e-7,
|
|
"temperature change contaminated acoustic volume source");
|
|
}
|
|
}
|
|
});
|
|
test("radiating mouth conserves reported mass at both orientations", [] {
|
|
for (bool left : {false,true}) {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f,.1f,10);
|
|
lossless(p[0]);
|
|
for (int i=0;i<10;++i) set(p[0],i,105000,1.2f);
|
|
BoundarySystem b(p);
|
|
b.add_radiating_open_end(0,left,101325,.001f);
|
|
b.add_closed_end(0,!left);
|
|
const double initial = mass(p[0]);
|
|
b.resolve(1e-6f);
|
|
const double exchange = -b.open_end_mass_flow(0)*1e-6;
|
|
p[0].advance(1e-6f);
|
|
near(mass(p[0])-initial,exchange,1e-11,"radiating mass flux mismatch");
|
|
}
|
|
});
|
|
test("radiating mouth reflects rarefaction at either pipe orientation", [] {
|
|
for (bool left : {false,true}) {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f,1,400);
|
|
lossless(p[0]);
|
|
const double c=std::sqrt(1.4*101325/1.2);
|
|
for(int i=0;i<400;++i) {
|
|
double x=(i+.5)/400.;
|
|
if(left) x=1-x;
|
|
float dp=100*std::exp(-std::pow((x-.3)/.04,2));
|
|
set(p[0],i,101325+dp,1.2f+dp/(c*c),(left?-1:1)*dp/(1.2*c));
|
|
}
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0,!left);
|
|
b.add_radiating_open_end(0,left,101325,.001f);
|
|
Solver solver(p,b); solver.minimum_substeps=1;
|
|
for(int i=0;i<int(1.4/c*44100);++i) solver.step(1.f/44100);
|
|
check(p[0].state(left?280:119).pressure < 101305,"open mouth did not return an inverted pressure pulse");
|
|
check(p[0].positivity_repairs==0,"radiation boundary required positivity repair");
|
|
}
|
|
});
|
|
test("radiating inflow supplies fresh ambient stagnation enthalpy", [] {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f,.1f,10);
|
|
for(int i=0;i<10;++i) set(p[0],i,99000,1.2f,0,.1f,.2f);
|
|
BoundarySystem b(p);
|
|
b.add_radiating_open_end(0,true,101325,.001f,300);
|
|
b.add_closed_end(0,false);
|
|
b.resolve(1e-6f);
|
|
const auto face=p[0].left_ghost;
|
|
check(face.velocity>0 && face.air==1 && face.fuel==0,"backflow supplied exhaust or fuel");
|
|
near(fluidsim::port::gamma/(fluidsim::port::gamma-1)*face.pressure/face.rho+.5*face.velocity*face.velocity,
|
|
fluidsim::port::gamma/(fluidsim::port::gamma-1)*gas_constant*300, .1,"ambient enthalpy mismatch");
|
|
});
|
|
test("tracer underflow preserves bulk gas and normal concentrations", [] {
|
|
for (float fraction : {1e-40f, 1e-30f, .5f}) {
|
|
Pipe1D p(.001f, 1, 8);
|
|
for (int i=0; i<8; ++i) set(p, i, 101325, 1.2f, 0, fraction);
|
|
p.left_ghost = p.state(0);
|
|
p.right_ghost = p.state(7);
|
|
const double before_mass = mass(p), before_energy = energy(p);
|
|
p.advance(1e-6f);
|
|
near(mass(p), before_mass, 0, "tracer truncation changed total mass");
|
|
near(energy(p), before_energy, 0, "tracer truncation changed energy");
|
|
for (auto cell : p.cells) {
|
|
check(std::fpclassify(cell.air) != FP_SUBNORMAL, "subnormal tracer survived");
|
|
near(cell.air, fraction < std::numeric_limits<float>::min() ? 0 : fraction,
|
|
fraction * 1e-6, "normal tracer concentration changed");
|
|
}
|
|
}
|
|
});
|
|
test("invalid uniform geometry", [] {
|
|
bool caught = false;
|
|
try {
|
|
Pipe1D p(0, .1f, 4);
|
|
} catch (const std::invalid_argument &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "zero area accepted");
|
|
});
|
|
test("closed pipe conservation", [] {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, 1, 50);
|
|
lossless(p[0]);
|
|
for (int i = 0; i < 50; i++)
|
|
set(p[0], i, i < 25 ? 120000 : 100000, 1.2f);
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0, true);
|
|
b.add_closed_end(0, false);
|
|
Solver s(p, b);
|
|
auto m = mass(p[0]), e = energy(p[0]);
|
|
for (int i = 0; i < 2000; i++)
|
|
s.step(1.f / 44100);
|
|
near(mass(p[0]), m, m * 2e-5, "closed mass drift");
|
|
near(energy(p[0]), e, e * 2e-5, "closed energy drift");
|
|
});
|
|
test("profiled pipe conservation", [] {
|
|
std::vector<float> a(40), l(40, .025f);
|
|
for (int i = 0; i < 40; i++)
|
|
a[i] = .001f * (1 + .5f * std::sin(i * .1f));
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(a, l);
|
|
lossless(p[0]);
|
|
for (int i = 0; i < 40; i++)
|
|
set(p[0], i, i < 20 ? 120000 : 100000, 1.2f);
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0, true);
|
|
b.add_closed_end(0, false);
|
|
Solver s(p, b);
|
|
auto m = mass(p[0]), e = energy(p[0]);
|
|
for (int i = 0; i < 1000; i++)
|
|
s.step(1.f / 44100);
|
|
near(mass(p[0]), m, m * 3e-5, "profile mass drift");
|
|
near(energy(p[0]), e, e * 3e-5, "profile energy drift");
|
|
});
|
|
test("volume pipe sealed conservation", [] {
|
|
Volume0D v(.001f, 180000, 400);
|
|
v.energy_relaxation_rate = 0;
|
|
auto &port = v.create_port();
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.0003f, .5f, 25);
|
|
lossless(p[0]);
|
|
BoundarySystem b(p);
|
|
float a = .0001f;
|
|
b.add_orifice(port, 0, true, a);
|
|
b.add_closed_end(0, false);
|
|
Solver s(p, b);
|
|
s.add_component(v);
|
|
auto m = mass(p[0]) + v.mass(), e = energy(p[0]) + vol_energy(v);
|
|
for (int i = 0; i < 2000; i++)
|
|
s.step(1.f / 44100);
|
|
near(mass(p[0]) + v.mass(), m, m * 3e-5, "coupling creates/deletes mass");
|
|
near(energy(p[0]) + vol_energy(v), e, e * 3e-5, "coupling creates/deletes energy");
|
|
});
|
|
test("equal pressure hot contact", [] {
|
|
Volume0D v(.001f, 101325, 600);
|
|
v.energy_relaxation_rate = 0;
|
|
auto &port = v.create_port();
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .2f, 20);
|
|
lossless(p[0]);
|
|
BoundarySystem b(p);
|
|
float a = .001f;
|
|
b.add_orifice(port, 0, true, a);
|
|
b.add_closed_end(0, false);
|
|
Solver s(p, b);
|
|
s.add_component(v);
|
|
double m = mass(p[0]) + v.mass();
|
|
for (int i = 0; i < 500; i++)
|
|
s.step(1.f / 44100);
|
|
near(mass(p[0]) + v.mass(), m, m * 2e-5, "contact produces mass");
|
|
near(v.pressure(), 101325, 5, "stationary contact pressure changes");
|
|
});
|
|
test("NASA choked nozzle mass flow", [] {
|
|
Atmosphere v;
|
|
v.pressure = 200000;
|
|
v.update_state(0);
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .2f, 10);
|
|
for (int i = 0; i < 10; i++)
|
|
set(p[0], i, 10000, 10000 / (gas_constant * 300));
|
|
BoundarySystem b(p);
|
|
float a = .0001f;
|
|
b.add_orifice(v.port, 0, true, a, .8f);
|
|
b.resolve(1e-8f);
|
|
double expected = .8 * a * 200000 / std::sqrt(300.) *
|
|
std::sqrt(fluidsim::port::gamma / gas_constant) *
|
|
std::pow(2 / (fluidsim::port::gamma + 1),
|
|
(fluidsim::port::gamma + 1) / (2 * (fluidsim::port::gamma - 1)));
|
|
near(-v.port.mass_flow, expected, expected * 2e-5,
|
|
"choked mdot differs from isentropic sonic-throat formula");
|
|
});
|
|
test("inertance units", [] {
|
|
Atmosphere v;
|
|
v.pressure = 102325;
|
|
v.update_state(0);
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .1f, 10);
|
|
BoundarySystem b(p);
|
|
float a = .0001f;
|
|
b.add_orifice_with_inertance(v.port, 0, true, a, 1, .05f);
|
|
b.resolve(1e-6f);
|
|
double expected = a / .05 * 1000 * 1e-6;
|
|
near(-v.port.mass_flow, expected, expected * 1e-4,
|
|
"mass-flow inertance must satisfy d(mdot)/dt=A*dp/L");
|
|
});
|
|
test("moving volume adiabatic law", [] {
|
|
Volume0D v(.001f);
|
|
v.energy_relaxation_rate = 0;
|
|
v.volume_change_rate = .0001f;
|
|
for (int i = 0; i < 1000; i++)
|
|
v.update(1e-3f);
|
|
near(v.volume(), .0011, 1e-8, "volume work did not change geometric volume");
|
|
near(v.pressure(), 101325 * std::pow(.001 / .0011, 1.4), 3, "p V^gamma not conserved");
|
|
});
|
|
test("profile at rest is well balanced", [] {
|
|
std::vector<float> a(80), l(80, .01f);
|
|
for (int i = 0; i < 80; i++)
|
|
a[i] = .001f * (1 + .7f * std::sin(.1f * i));
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(a, l);
|
|
lossless(p[0]);
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0, true);
|
|
b.add_closed_end(0, false);
|
|
Solver s(p, b);
|
|
for (int j = 0; j < 1000; j++)
|
|
s.step(1.f / 44100);
|
|
for (int i = 0; i < 80; i++) {
|
|
near(p[0].state(i).velocity, 0, 1e-3, "stationary nozzle accelerates");
|
|
near(p[0].state(i).pressure, 101325, 2, "stationary nozzle pressure drifts");
|
|
}
|
|
});
|
|
test("pipe connection all four orientations", [] {
|
|
for (bool la : {false, true})
|
|
for (bool lb : {false, true}) {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .3f, 20);
|
|
p.emplace_back(.001f, .3f, 20);
|
|
for (auto &q : p)
|
|
lossless(q);
|
|
for (int i = 0; i < 20; i++)
|
|
set(p[0], i, 130000, 1.4f, 0, .2f);
|
|
BoundarySystem b(p);
|
|
b.add_connection(0, la, 1, lb);
|
|
b.add_closed_end(0, !la);
|
|
b.add_closed_end(1, !lb);
|
|
Solver s(p, b);
|
|
double m = mass(p[0]) + mass(p[1]), e = energy(p[0]) + energy(p[1]);
|
|
for (int i = 0; i < 2000; i++)
|
|
s.step(1.f / 44100);
|
|
near(mass(p[0]) + mass(p[1]), m, m * 4e-5, "junction mass drift");
|
|
near(energy(p[0]) + energy(p[1]), e, e * 4e-5, "junction energy drift");
|
|
for (auto &q : p)
|
|
check(q.positivity_repairs == 0, "junction positivity repair");
|
|
}
|
|
});
|
|
test("split pipe matches unsplit pipe", [] {
|
|
std::vector<Pipe1D> a, b;
|
|
a.emplace_back(.001f, 1, 80);
|
|
b.emplace_back(.001f, .5f, 40);
|
|
b.emplace_back(.001f, .5f, 40);
|
|
lossless(a[0]);
|
|
lossless(b[0]);
|
|
lossless(b[1]);
|
|
for (int i = 0; i < 80; i++) {
|
|
float pr = 101325 + 100 * std::exp(-std::pow((i - 30) / 5.f, 2));
|
|
set(a[0], i, pr, 1.2f);
|
|
set(b[i / 40], i % 40, pr, 1.2f);
|
|
}
|
|
BoundarySystem ba(a), bb(b);
|
|
ba.add_closed_end(0, true);
|
|
ba.add_closed_end(0, false);
|
|
bb.add_closed_end(0, true);
|
|
bb.add_closed_end(1, false);
|
|
bb.add_connection(0, false, 1, true);
|
|
Solver sa(a, ba), sb(b, bb);
|
|
for (int j = 0; j < 300; j++) {
|
|
sa.step(1.f / 44100);
|
|
sb.step(1.f / 44100);
|
|
}
|
|
for (int i = 0; i < 80; i++)
|
|
near(a[0].state(i).pressure, b[i / 40].state(i % 40).pressure, .1,
|
|
"split junction changes solution");
|
|
});
|
|
test("three volumes shared donor conservation", [] {
|
|
Volume0D a(.000001f, 1e6f, 500), b(.001f), c(.001f);
|
|
a.energy_relaxation_rate = b.energy_relaxation_rate = c.energy_relaxation_rate = 0;
|
|
auto &ap = a.create_port();
|
|
auto &ap2 = a.create_port();
|
|
auto &bp = b.create_port();
|
|
auto &cp = c.create_port();
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .1f, 10);
|
|
p.emplace_back(.001f, .1f, 10);
|
|
for (auto &q : p)
|
|
lossless(q);
|
|
BoundarySystem bs(p);
|
|
float ar = .001f;
|
|
bs.add_orifice(ap, 0, true, ar);
|
|
bs.add_orifice(ap2, 1, true, ar);
|
|
bs.add_orifice(bp, 0, false, ar);
|
|
bs.add_orifice(cp, 1, false, ar);
|
|
Solver s(p, bs);
|
|
s.add_component(a);
|
|
s.add_component(b);
|
|
s.add_component(c);
|
|
double m = mass(p[0]) + mass(p[1]) + a.mass() + b.mass() + c.mass();
|
|
double e = energy(p[0]) + energy(p[1]) + vol_energy(a) + vol_energy(b) + vol_energy(c);
|
|
for (int j = 0; j < 3000; j++)
|
|
s.step(1.f / 44100);
|
|
near(mass(p[0]) + mass(p[1]) + a.mass() + b.mass() + c.mass(), m, m * 5e-5,
|
|
"shared donor mass drift");
|
|
near(energy(p[0]) + energy(p[1]) + vol_energy(a) + vol_energy(b) + vol_energy(c), e, e * 5e-5,
|
|
"shared donor energy drift");
|
|
check(a.pressure() > 0, "donor vacuum");
|
|
});
|
|
test("open boundary reports actual conserved flux both ends", [] {
|
|
for (bool left : {false, true}) {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .1f, 10);
|
|
lossless(p[0]);
|
|
for (int i = 0; i < 10; i++)
|
|
set(p[0], i, 105000, 1.2f);
|
|
BoundarySystem b(p);
|
|
b.add_open_end(0, left, 101325, .001f);
|
|
b.add_closed_end(0, !left);
|
|
double m = mass(p[0]);
|
|
b.resolve(1e-6f);
|
|
double dm = -b.open_end_mass_flow(0) * 1e-6;
|
|
p[0].advance(1e-6f);
|
|
near(mass(p[0]) - m, dm, 1e-11,
|
|
"reported radiation flow differs from finite-volume mass exchange");
|
|
}
|
|
});
|
|
test("duplicate and invalid boundary rejection", [] {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, 1, 10);
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0, true);
|
|
bool caught = false;
|
|
try {
|
|
b.add_open_end(0, true, 101325, .001f);
|
|
} catch (const std::invalid_argument &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "duplicate face accepted");
|
|
caught = false;
|
|
try {
|
|
b.add_closed_end(9, true);
|
|
} catch (const std::invalid_argument &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "invalid pipe accepted");
|
|
});
|
|
test("supersonic outlet extrapolates", [] {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, 1, 10);
|
|
for (int i = 0; i < 10; i++)
|
|
set(p[0], i, 101325, 1.2f, 700);
|
|
BoundarySystem b(p);
|
|
b.add_open_end(0, false, 50000, .001f);
|
|
b.resolve(1e-6f);
|
|
near(b.open_end_mass_flow(0), .001 * 1.2 * 700, 1e-5,
|
|
"supersonic outlet responds to downstream pressure");
|
|
});
|
|
test("underdefined supersonic inlet rejected", [] {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, 1, 10);
|
|
set(p[0], 0, 101325, 1.2f, 700);
|
|
BoundarySystem b(p);
|
|
b.add_open_end(0, true, 101325, .001f);
|
|
bool caught = false;
|
|
try {
|
|
b.resolve(1e-6f);
|
|
} catch (const std::runtime_error &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "supersonic inlet silently accepted pressure-only data");
|
|
});
|
|
test("closed and open acoustic reflection signs", [] {
|
|
for (bool closed : {false, true}) {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, 1, 400);
|
|
lossless(p[0]);
|
|
const double c = std::sqrt(1.4 * 101325 / 1.2);
|
|
for (int i = 0; i < 400; i++) {
|
|
double x = (i + .5) / 400.;
|
|
float dp = 100 * std::exp(-std::pow((x - .3) / .04, 2));
|
|
set(p[0], i, 101325 + dp, 1.2f + dp / (c * c), dp / (1.2 * c));
|
|
}
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0, true);
|
|
if (closed)
|
|
b.add_closed_end(0, false);
|
|
else
|
|
b.add_open_end(0, false, 101325, .001f);
|
|
Solver s(p, b);
|
|
s.minimum_substeps = 1;
|
|
const int steps = int((1.4 / c) * 44100);
|
|
for (int j = 0; j < steps; j++)
|
|
s.step(1.f / 44100);
|
|
float dp = p[0].state(119).pressure - 101325;
|
|
check(closed ? dp > 30 : dp < -30, "wrong reflected pressure sign/amplitude");
|
|
check(p[0].positivity_repairs == 0, "acoustic positivity repair");
|
|
}
|
|
});
|
|
test("Sod shock tube convergence", [] {
|
|
double previous = 1e9;
|
|
for (int n : {100, 200, 400}) {
|
|
Pipe1D p(1, 1, n);
|
|
lossless(p);
|
|
for (int i = 0; i < n; i++)
|
|
set(p, i, i < n / 2 ? 1 : .1f, i < n / 2 ? 1 : .125f);
|
|
double t = 0;
|
|
while (t < .2) {
|
|
p.left_ghost = p.state(0);
|
|
p.right_ghost = p.state(n - 1);
|
|
float dt = std::min(.2 - t, double(p.stable_timestep(.7f)));
|
|
p.advance(dt);
|
|
t += dt;
|
|
}
|
|
double error = 0;
|
|
for (int i = 0; i < n; i++) {
|
|
double xi = ((i + .5) / n - .5) / .2;
|
|
double rho;
|
|
if (xi < -std::sqrt(1.4))
|
|
rho = 1;
|
|
else if (xi < -.0702728)
|
|
rho = std::pow((2 / (1.4 + 1)) + (1.4 - 1) / (1.4 + 1) * (-xi / std::sqrt(1.4)),
|
|
2 / (1.4 - 1));
|
|
else if (xi < .9274526)
|
|
rho = .4263194;
|
|
else if (xi < 1.7521557)
|
|
rho = .2655737;
|
|
else
|
|
rho = .125;
|
|
error += std::abs(p.state(i).rho - rho) / n;
|
|
}
|
|
std::cout << " Sod n=" << n << " density L1=" << error << '\n';
|
|
check(error < previous * .9, "shock density error does not decrease on refinement");
|
|
previous = error;
|
|
check(p.positivity_repairs == 0, "Sod positivity repair");
|
|
}
|
|
check(previous < .025, "Sod density error excessive");
|
|
});
|
|
test("strong pulse and near-vacuum rarefaction", [] {
|
|
for (bool vacuum : {false, true}) {
|
|
Pipe1D p(1, 1, 200);
|
|
lossless(p);
|
|
for (int i = 0; i < 200; i++)
|
|
set(p, i, vacuum ? .4f : (i < 100 ? 1000 : .01f), 1, vacuum ? (i < 100 ? -2 : 2) : 0);
|
|
double t = 0;
|
|
while (t < .005) {
|
|
p.left_ghost = p.state(0);
|
|
p.right_ghost = p.state(199);
|
|
float dt = std::min(.005 - t, double(p.stable_timestep(.5f)));
|
|
p.advance(dt);
|
|
t += dt;
|
|
}
|
|
for (int i = 0; i < 200; i++)
|
|
check(p.state(i).rho > 0 && p.state(i).pressure > 0, "negative shock/vacuum state");
|
|
check(p.positivity_repairs == 0, "stress case required conservation-breaking positivity repair");
|
|
}
|
|
});
|
|
|
|
test("small-signal Helmholtz resonance", [] {
|
|
std::vector<Pipe1D> p;
|
|
float area = .000314159f;
|
|
p.emplace_back(area, .05f, 40);
|
|
Volume0D v(.001f, 101335, 300);
|
|
auto &port = v.create_port();
|
|
BoundarySystem b(p);
|
|
b.add_orifice(port, 0, true, area, 1);
|
|
b.add_open_end(0, false, 101325, area);
|
|
Solver solver(p, b);
|
|
solver.minimum_substeps = 1;
|
|
solver.add_component(v);
|
|
double prev = v.pressure() - 101325, first = 0, last = 0;
|
|
int count = 0;
|
|
for (int i = 0; i < 4410; i++) {
|
|
solver.step(1.f / 44100);
|
|
double now = v.pressure() - 101325;
|
|
check(std::abs(now) < 15, "small-signal cavity grows instead of oscillating/decaying");
|
|
if (prev > 0 && now <= 0 && i > 441) {
|
|
double t = (i + prev / (prev - now)) / 44100;
|
|
if (!count)
|
|
first = t;
|
|
last = t;
|
|
count++;
|
|
}
|
|
prev = now;
|
|
}
|
|
double frequency = (count - 1) / (last - first);
|
|
double expected =
|
|
std::sqrt(1.4 * 287 * 300) / (2 * 3.141592653589793) * std::sqrt(area / (.001 * .05));
|
|
std::cout << " Helmholtz frequency=" << frequency << " analytic lumped=" << expected << '\n';
|
|
check(count > 5, "not enough cavity oscillations");
|
|
near(frequency, expected, expected * .05, "Helmholtz mode frequency error exceeds 5%");
|
|
check(p[0].positivity_repairs == 0, "cavity requires positivity repairs");
|
|
});
|
|
test("invalid volume and timestep rejection", [] {
|
|
bool caught = false;
|
|
try {
|
|
Volume0D v(-1);
|
|
} catch (const std::invalid_argument &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "negative volume accepted");
|
|
Pipe1D p(.001f, .1f, 10);
|
|
caught = false;
|
|
try {
|
|
p.required_substeps(std::nanf(""));
|
|
} catch (const std::invalid_argument &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "NaN timestep accepted");
|
|
});
|
|
test("fresh-air species conservation across connection", [] {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .5f, 40);
|
|
p.emplace_back(.001f, .5f, 40);
|
|
for (int k = 0; k < 2; k++) {
|
|
lossless(p[k]);
|
|
for (int i = 0; i < 40; i++)
|
|
set(p[k], i, k ? 100000 : 120000, 1.2f, 0, k ? .9f : .1f,
|
|
k ? .02f : .05f);
|
|
}
|
|
auto air_mass = [&] {
|
|
double sum = 0;
|
|
for (auto &q : p)
|
|
for (int i = 0; i < 40; i++)
|
|
sum += q.cell_mass(i) * q.cells[i].air;
|
|
return sum;
|
|
};
|
|
auto fuel_mass = [&] {
|
|
double sum = 0;
|
|
for (auto &q : p)
|
|
for (int i = 0; i < 40; i++)
|
|
sum += q.cell_mass(i) * q.cells[i].fuel;
|
|
return sum;
|
|
};
|
|
double before = air_mass();
|
|
double fuel_before = fuel_mass();
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0, true);
|
|
b.add_closed_end(1, false);
|
|
b.add_connection(0, false, 1, true);
|
|
Solver solver(p, b);
|
|
for (int i = 0; i < 3000; i++)
|
|
solver.step(1.f / 44100);
|
|
near(air_mass(), before, before * 4e-5, "species created at shared face");
|
|
near(fuel_mass(), fuel_before, fuel_before * 4e-5, "fuel species created at shared face");
|
|
});
|
|
test("unsupported unequal junction rejected", [] {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .5f, 20);
|
|
p.emplace_back(.002f, .5f, 20);
|
|
BoundarySystem b(p);
|
|
bool caught = false;
|
|
try {
|
|
b.add_connection(0, false, 1, true);
|
|
} catch (const std::invalid_argument &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "unvalidated sudden area junction silently accepted");
|
|
});
|
|
test("unconnected endpoint rejected by solver", [] {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .5f, 20);
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0, true);
|
|
Solver solver(p, b);
|
|
bool caught = false;
|
|
try {
|
|
solver.step(1.f / 44100);
|
|
} catch (const std::runtime_error &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "unconnected face silently used default ghost");
|
|
});
|
|
test("closed duct acoustic eigenfrequency", [] {
|
|
for (int n : {50, 100, 200}) {
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, 1, n);
|
|
double c = std::sqrt(1.4 * 101325 / 1.2);
|
|
for (int i = 0; i < n; i++) {
|
|
double dp = 100 * std::cos(3.141592653589793 * (i + .5) / n);
|
|
set(p[0], i, 101325 + dp, 1.2f + dp / (c * c));
|
|
}
|
|
BoundarySystem b(p);
|
|
b.add_closed_end(0, true);
|
|
b.add_closed_end(0, false);
|
|
Solver solver(p, b);
|
|
solver.minimum_substeps = 1;
|
|
double prev = 1, first = 0, last = 0;
|
|
int count = 0;
|
|
for (int j = 0; j < 1764; j++) {
|
|
solver.step(1.f / 44100);
|
|
double projection = 0;
|
|
for (int i = 0; i < n; i++)
|
|
projection +=
|
|
(p[0].state(i).pressure - 101325) * std::cos(3.141592653589793 * (i + .5) / n);
|
|
if (prev > 0 && projection <= 0) {
|
|
double t = (j + prev / (prev - projection)) / 44100;
|
|
if (!count)
|
|
first = t;
|
|
last = t;
|
|
count++;
|
|
}
|
|
prev = projection;
|
|
}
|
|
double f = (count - 1) / (last - first);
|
|
std::cout << " closed mode n=" << n << " f=" << f << " analytic=" << c / 2 << '\n';
|
|
check(count > 4, "insufficient acoustic cycles");
|
|
near(f, c / 2, c / 2 * .01, "closed duct frequency exceeds 1% error");
|
|
}
|
|
});
|
|
test("inertance impulse and dissipative resistance", [] {
|
|
Atmosphere v;
|
|
v.pressure = 102325;
|
|
v.update_state(0);
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .1f, 10);
|
|
BoundarySystem b(p);
|
|
float area = .0001f;
|
|
b.add_orifice_with_inertance(v.port, 0, true, area, 1, .05f, 1);
|
|
for (int j = 0; j < 10000; j++)
|
|
b.resolve(.0001f);
|
|
double equilibrium = std::sqrt(2 * v.density() * area * area * 1000);
|
|
near(-v.port.mass_flow, equilibrium, equilibrium * .001,
|
|
"inertance resistance steady flow has wrong units or sign");
|
|
});
|
|
test("incompatible gas coupling rejected", [] {
|
|
Volume0D v(.001f, 101325, 300, 300, 1.3f);
|
|
auto &port = v.create_port();
|
|
std::vector<Pipe1D> p;
|
|
p.emplace_back(.001f, .1f, 10);
|
|
BoundarySystem b(p);
|
|
float area = .0001f;
|
|
bool caught = false;
|
|
try {
|
|
b.add_orifice(port, 0, true, area);
|
|
} catch (const std::invalid_argument &) {
|
|
caught = true;
|
|
}
|
|
check(caught, "different gamma gases silently coupled");
|
|
});
|
|
test("thermal relaxation matches analytic energy decay", [] {
|
|
Volume0D v(.001f, 101325, 600);
|
|
v.energy_relaxation_rate = 2;
|
|
double initial = vol_energy(v), target = v.mass() * 287 * 300 / (1.4 - 1);
|
|
for (int i = 0; i < 1000; i++)
|
|
v.update(.001f);
|
|
near(vol_energy(v), target + (initial - target) * std::exp(-2), initial * 3e-5,
|
|
"thermal relaxation differs from exact exponential");
|
|
});
|
|
|
|
test("pipe solver scratch follows runtime solver changes", [] {
|
|
Pipe1D pipe(.001f,.2f,16);
|
|
pipe.left_ghost=pipe.state(0);
|
|
pipe.right_ghost=pipe.state(15);
|
|
pipe.solver=PipeSolver::Rusanov;
|
|
pipe.prepare_step();
|
|
const auto rusanov_bytes=pipe.scratch_bytes();
|
|
pipe.solver=PipeSolver::Hlle;
|
|
pipe.prepare_step();
|
|
const auto hlle_bytes=pipe.scratch_bytes();
|
|
pipe.solver=PipeSolver::Muscl;
|
|
pipe.prepare_step();
|
|
const auto muscl_bytes=pipe.scratch_bytes();
|
|
check(rusanov_bytes<hlle_bytes && hlle_bytes<muscl_bytes,
|
|
"solver-specific work buffers were not allocated on demand");
|
|
for(auto solver : {PipeSolver::Rusanov,PipeSolver::Hlle,PipeSolver::Hllc,
|
|
PipeSolver::Muscl,PipeSolver::Rusanov}) {
|
|
pipe.solver=solver;
|
|
for(int step=0;step<8;++step) {
|
|
pipe.left_ghost=pipe.state(0);
|
|
pipe.right_ghost=pipe.state(15);
|
|
const float dt=.5f*pipe.prepare_step();
|
|
pipe.advance(dt);
|
|
}
|
|
}
|
|
for(int i=0;i<16;++i) {
|
|
const auto state=pipe.state(i);
|
|
check(std::isfinite(state.rho)&&std::isfinite(state.pressure)&&
|
|
state.rho>0&&state.pressure>0,"solver switch produced a nonphysical pipe state");
|
|
}
|
|
});
|
|
|
|
test("sealed crankcase compression conserves mass above five bar", [] {
|
|
Crankshaft crank;
|
|
crank.angle = 0;
|
|
Crankcase v(crank, .02725f, .11f, .054f, 1e-5f);
|
|
double m = v.mass(), invariant = v.pressure() * std::pow(v.volume(), 1.4);
|
|
for (int i = 1; i <= 1000; i++) {
|
|
crank.angle = pi * i / 1000;
|
|
v.pre_step(1e-5f);
|
|
v.update_state(1e-5f);
|
|
}
|
|
near(v.mass(), m, m * 1e-6, "crankcase pressure cap deletes mass");
|
|
near(v.pressure() * std::pow(v.volume(), 1.4), invariant, invariant * 1e-4,
|
|
"crankcase compression violates adiabatic invariant");
|
|
check(v.pressure() > 5e5, "compression did not exercise old hard pressure cap");
|
|
});
|
|
std::cout << total - failed << "/" << total << " passed\n";
|
|
return failed ? 1 : 0;
|
|
}
|