diff --git a/src/coolant/fluid_t.hpp b/src/coolant/fluid_t.hpp index b914bd0..87fa23e 100644 --- a/src/coolant/fluid_t.hpp +++ b/src/coolant/fluid_t.hpp @@ -10,9 +10,11 @@ struct fluid_t { const double gPl; // g/L const double gPmol; // g/mol + const double jPg; // J/g latent heat of vaporisation + const double jPgk; // J/g/K heat capacity const coolant::vapor_pressure vapor_pressure; - constexpr fluid_t(double gPl, double gPmol, coolant::vapor_pressure vapor_pressure) : gPl(gPl), gPmol(gPmol), vapor_pressure(vapor_pressure) { } + constexpr fluid_t(double gPl, double gPmol, double jPg, double jPgk, coolant::vapor_pressure vapor_pressure) : gPl(gPl), gPmol(gPmol), jPg(jPg), jPgk(jPgk), vapor_pressure(vapor_pressure) { } constexpr double g_to_mol(double g) const { return g / gPmol; } constexpr double mol_to_g(double mol) const { return mol * gPmol; } @@ -22,7 +24,7 @@ struct fluid_t constexpr double l_to_mol(double l) const { return g_to_mol(l_to_g(l)); } }; -constexpr const fluid_t WATER = fluid_t(1000, 18, {8.07131, 1730.63, 233.426}); +constexpr const fluid_t WATER = fluid_t(1000, 18, 2257, 4.1816, {8.07131, 1730.63, 233.426}); } diff --git a/src/main.cpp b/src/main.cpp index 5542917..5a64e57 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -37,7 +37,7 @@ int main() sim::reactor::coolant::pipe(vessel), { "# ##", "#FCF ", - " C C ", + " CHC ", " FCF#", "## #" }); @@ -89,8 +89,12 @@ secs: ss << s << "s\n"; for(int i = 0; i < skip; i++) { - reactor.update(speed / framerate / skip); + double dt = speed / framerate / skip; + + reactor.update(dt); vessel.update(); + + vessel.extract_steam(dt, 0.01, 101000); } secs += speed / framerate; diff --git a/src/reactor/coolant/vessel.cpp b/src/reactor/coolant/vessel.cpp index 1521ac5..84d8a8b 100644 --- a/src/reactor/coolant/vessel.cpp +++ b/src/reactor/coolant/vessel.cpp @@ -19,27 +19,73 @@ void vessel::update() double T = conversions::temperature::c_to_k(heat); double n = fluid.mol_to_g((V * P) / (T * constants::R)) - steam; - steam += n; - level -= fluid.g_to_l(n); + double s = steam + n; + double l = fluid.l_to_g(level) - n; - if(fluid.g_to_l(steam) > volume) + if(l < 0) { - steam = fluid.l_to_g(volume); - level = 0; + s += l; + l = 0; } + + if(fluid.g_to_l(s) > volume) + { + s = fluid.l_to_g(volume); + l = 0; + } + + double diff = s - steam; + + steam = s; + level = fluid.g_to_l(l); + heat -= diff * fluid.jPg / (fluid.l_to_g(level) + steam) / fluid.jPgk; } double vessel::add_heat(double t1) { double t2 = get_heat(); double t = t1 - t2; - double m1 = 1000; - double m2 = level + fluid.g_to_l(steam); + double m1 = 1000000; + double m2 = (fluid.l_to_g(level) + steam) * fluid.jPgk; double m = m1 + m2; return heat = t1 - t * m2 / m; } +double vessel::extract_steam(double dt, double a, double p2) +{ + // calculate the mass moved + double p1 = get_pressure(); + double p = p1 - p2; + double m = 1; + + if(p == 0) + { + return 0; + } + + if(p < 0) + { + m = -1; + p = -p; + + return 0; + } + + double V = (volume - level) * 0.001; + double v = std::sqrt( V * p / steam ); + double mass = m * dt * a * p / v; + + if(mass > steam) + { + mass = steam; + } + + steam -= mass; + + return mass; +} + double vessel::get_pressure() const { double T = conversions::temperature::c_to_k(heat); diff --git a/src/reactor/coolant/vessel.hpp b/src/reactor/coolant/vessel.hpp index 25a4bec..e6d575f 100644 --- a/src/reactor/coolant/vessel.hpp +++ b/src/reactor/coolant/vessel.hpp @@ -23,6 +23,7 @@ public: void update(); double add_heat(double amount); + double extract_steam(double dt, double a, double p2); constexpr double get_volume() const { return volume; } constexpr double get_level() const { return level; }