diff --git a/CMakeLists.txt b/CMakeLists.txt index b09541e..a9bfd3a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.25) project(FastNuclearSim VERSION 1.0) set(CMAKE_CXX_STANDARD 26) -set(CMAKE_CXX_FLAGS "-g -O3 -I/usr/include/freetype2") +set(CMAKE_CXX_FLAGS "-g -I/usr/include/freetype2") file(GLOB_RECURSE SOURCES src/*.cpp) diff --git a/nasty_function.c++ b/nasty_function.c++ new file mode 100644 index 0000000..a798464 --- /dev/null +++ b/nasty_function.c++ @@ -0,0 +1,51 @@ + +#include + +double calc(double& val1, double& val2) +{ + double J_m2 = J_m*J_m; + double L_m2 = L_m*L_m; + double V_t2 = V_t*V_t; + double n_g2 = n_g*n_g; + double n_l2 = n_l*n_l; + double K2 = K*K; + double R2 = R*R; + + double part1 = ( + -J_m*R*n_g + + 2*K*R*n_g + + K*R*n_l - + L_m*n_g - + V_t + ); + + double part2 = std::sqrt( + J_m2*R2*n_g2 + + 2*J_m*K*R2*n_g*n_l + + 2*J_m*L_m*R*n_g2 + + 4*J_m*L_m*R*n_g*n_l + + 4*J_m*L_m*R*n_l2 - + 2*J_m*R*V_t*n_g - + 4*J_m*R*V_t*n_l + + K2*R2*n_l2 - + 2*K*L_m*R*n_g*n_l - + 4*K*L_m*R*n_l2 + + 2*K*R*V_t*n_l + + L_m2*n_g2 + + 4*L_m2*n_g*n_l + + 4*L_m2*n_l2 - + 2*L_m*V_t*n_g - + 4*L_m*V_t*n_l + + V_t2 + ); + + double part3 = (2*( + J_m*R - + K*R + + L_m + )); + + val1 = (part1 - part2) / part3; + val2 = (part1 + part2) / part3; +} + diff --git a/src/conversions/pressure.hpp b/src/conversions/pressure.hpp index 99bf124..b79bdc8 100644 --- a/src/conversions/pressure.hpp +++ b/src/conversions/pressure.hpp @@ -6,6 +6,8 @@ namespace sim::conversions::pressure constexpr double torr_to_pa(double t) { return t * 133.322; } constexpr double pa_to_torr(double p) { return p / 133.322; } +constexpr double mmhg_to_pa(double t) { return t * 133.3224; } +constexpr double pa_to_mmhg(double p) { return p / 133.3224; } }; diff --git a/src/coolant/condenser.cpp b/src/coolant/condenser.cpp index b80ac56..c6d1cac 100644 --- a/src/coolant/condenser.cpp +++ b/src/coolant/condenser.cpp @@ -21,7 +21,7 @@ condenser::condenser(fluid_t type, double height, double diameter, double mass, void condenser::update(double secs) { - ((sim::coolant::fluid_holder*)this)->update(secs); + ((fluid_holder*)this)->update(secs); } diff --git a/src/coolant/fluid_holder.cpp b/src/coolant/fluid_holder.cpp index c58d1d8..fd7c376 100644 --- a/src/coolant/fluid_holder.cpp +++ b/src/coolant/fluid_holder.cpp @@ -94,6 +94,13 @@ double fluid_holder::calc_pressure_mol(double heat, double volume, double pressu return (V * pressure) / (T * constants::R); } +double fluid_holder::calc_pressure_vol(double heat, double pressure, double mol) +{ + double T = conversions::temperature::c_to_k(heat); + + return 1000 * (mol * T * constants::R) / pressure; +} + double fluid_holder::get_pressure() const { return calc_pressure(heat, get_steam_volume(), fluid.g_to_mol(get_steam())); @@ -111,32 +118,28 @@ void fluid_holder::update(double secs) if(mass > 0) { - double V = get_steam_volume() * 0.001; - double P = fluid.vapor_pressure.calc_p(heat); - double T = conversions::temperature::c_to_k(heat); - double n = fluid.mol_to_g((V * P) / (T * constants::R)) - steam; + // use ideal gas law to get target steam density in mol/L + double heat_k = conversions::temperature::c_to_k(heat); + double target_pressure = fluid.vapor_pressure.calc_p(heat); + double density = target_pressure / (constants::R * heat_k) / 1000; + + double m_c = fluid.l_to_mol(1); + double n_t = fluid.l_to_mol(level) + fluid.g_to_mol(steam); + double v_l = (n_t - density * volume) / (m_c - density); + double n_l = fluid.l_to_mol(v_l); - double s = steam + n; - double l = fluid.l_to_g(level) - n; - double v = fluid.l_to_g(volume); - - if(l < 0) + if(n_l < 0) { - s += l; - l = 0; + v_l = 0; + n_l = 0; } - if(l > v) - { - l = v; - s = 0; - } + double n_diff = n_l - fluid.l_to_mol(level); + double steam_add = -fluid.mol_to_g(n_diff); - double diff = s - steam; - - steam = s; - level = fluid.g_to_l(l); - heat -= diff * fluid.jPg / mass; + level += fluid.mol_to_l(n_diff); + steam += steam_add; + heat -= steam_add * fluid.jPg / mass; } } diff --git a/src/coolant/fluid_holder.hpp b/src/coolant/fluid_holder.hpp index e2ca96c..a132544 100644 --- a/src/coolant/fluid_holder.hpp +++ b/src/coolant/fluid_holder.hpp @@ -8,7 +8,7 @@ namespace sim::coolant class fluid_holder { -protected: + public://protected: double level = 0; // litres double steam = 0; // grams @@ -40,6 +40,7 @@ public: static double calc_pressure(double temp, double volume, double mass); static double calc_pressure_mol(double temp, double volume, double pressure); + static double calc_pressure_vol(double heat, double pressure, double mol); void update(double dt); }; diff --git a/src/coolant/fluid_t.hpp b/src/coolant/fluid_t.hpp index 1e7ce8d..ed3bdc8 100644 --- a/src/coolant/fluid_t.hpp +++ b/src/coolant/fluid_t.hpp @@ -31,7 +31,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, 2257, 4.1816, {8.07131, 1730.63, 233.426}); +constexpr const fluid_t WATER = fluid_t(1000, 18, 2257, 4.1816, {8.07131, 1730.63, 233.426, 108.266, 8.14019, 1810.94, 244.485}); } diff --git a/src/coolant/pump.cpp b/src/coolant/pump.cpp index 0a43748..c951d90 100644 --- a/src/coolant/pump.cpp +++ b/src/coolant/pump.cpp @@ -100,7 +100,7 @@ void pump::update(double dt) double p_diff = (p_diff_1 + p_diff_2) / 2; double work = p_diff * dst_volume * 0.001 + get_rpm() * 60 * dt * friction; - velocity -= calc_work(work, mass); + velocity = std::max(velocity - calc_work(work, mass), 0.0); flow = dst_volume / dt; } diff --git a/src/coolant/valve.cpp b/src/coolant/valve.cpp index 1fa80ba..2a4f1b8 100644 --- a/src/coolant/valve.cpp +++ b/src/coolant/valve.cpp @@ -42,10 +42,21 @@ void valve::update(double dt) double m = max * state; double diff = (pressure1 - pressure2); double remove = diff - diff * std::pow(1 - m, dt); - - double mol = fluid_holder::calc_pressure_mol(src->get_heat(), src->get_steam_volume(), pressure1 - remove); - double mass = src->get_steam() - src->fluid.mol_to_g(mol); + double mol, mass; + + if(remove < 0) + { + mol = fluid_holder::calc_pressure_mol(src->get_heat(), src->get_steam_volume(), pressure1 - remove); + mass = src->get_steam() - src->fluid.mol_to_g(mol); + } + + else + { + mol = fluid_holder::calc_pressure_mol(dst->get_heat(), dst->get_steam_volume(), pressure2 - remove); + mass = dst->get_steam() - dst->fluid.mol_to_g(mol); + } + double heat1 = src->get_heat(); // C double heat2 = dst->get_heat(); diff --git a/src/coolant/vapor_pressure.cpp b/src/coolant/vapor_pressure.cpp index b025f46..4089ac2 100644 --- a/src/coolant/vapor_pressure.cpp +++ b/src/coolant/vapor_pressure.cpp @@ -10,11 +10,33 @@ using namespace sim::conversions; double vapor_pressure::calc_p(double t) const { - return pressure::torr_to_pa(std::pow(10, A - B / ( C + t ) )); + double p; + + if(t < T) + { + p = std::pow(10, A1 - B1 / ( C1 + t ) ); + } + + else + { + p = std::pow(10, A2 - B2 / ( C2 + t ) ); + } + + return pressure::mmhg_to_pa(p); } double vapor_pressure::calc_t(double p) const { - return B / ( A - std::log(pressure::pa_to_torr(p)) / std::log(10) ) - C; + double P = pressure::pa_to_mmhg(p); + + if(p < calc_p(T)) + { + return B1 / ( A1 - std::log(P) / std::log(10) ) - C1; + } + + else + { + return B2 / ( A2 - std::log(P) / std::log(10) ) - C2; + } } diff --git a/src/coolant/vapor_pressure.hpp b/src/coolant/vapor_pressure.hpp index 46ceb86..464e3ae 100644 --- a/src/coolant/vapor_pressure.hpp +++ b/src/coolant/vapor_pressure.hpp @@ -6,9 +6,12 @@ namespace sim::coolant struct vapor_pressure { - const double A, B, C; + const double A1, B1, C1; + const double A2, B2, C2; + const double T; - constexpr vapor_pressure(double A, double B, double C) : A(A), B(B), C(C) { } + constexpr vapor_pressure(double A1, double B1, double C1, double T, double A2, double B2, double C2) : + A1(A1), B1(B1), C1(C1), T(T), A2(A2), B2(B2), C2(C2) { } double calc_p(double t) const; double calc_t(double p) const; diff --git a/src/graphics/monitor/primary_loop.cpp b/src/graphics/monitor/primary_loop.cpp index 798a92e..4c1f567 100644 --- a/src/graphics/monitor/primary_loop.cpp +++ b/src/graphics/monitor/primary_loop.cpp @@ -105,29 +105,35 @@ void primary_loop::init() void primary_loop::update(double dt) { - std::stringstream ss; - sim::graphics::mesh rmesh; system& sys = sim::system::active; + clock_now += dt; - ss << "\n\n"; - ss << show( sys.turbine_bypass_valve->get_state() * 100 ) << " %\n"; - ss << show( sys.turbine_bypass_valve->get_flow() / 1000 ) << " kg/s\n"; - ss << "\n\n\n"; - ss << show( sys.turbine_inlet_valve->get_state() * 100 ) << " %\n"; - ss << show( sys.turbine_inlet_valve->get_flow() / 1000 ) << " kg/s\n"; - ss << "\n\n\n"; - ss << show( sys.primary_pump->get_power() * 100 ) << " %\n"; - ss << show( sys.primary_pump->get_rpm() ) << " r/min\n"; - ss << show( sys.primary_pump->get_flow_mass() / 1000 ) << " kg/s\n"; - ss << "\n\n\n"; - ss << show( sys.condenser->get_heat() ) << " C\n"; - ss << show( sys.condenser->get_steam() ) << " g\n"; - ss << show( sys.condenser->get_pressure() / 1000 ) << " kPa\n"; - ss << show( sys.condenser->get_level() / 1000 ) << " / " << show( sys.condenser->get_volume() / 1000 ) << " kL\n"; + if(clock_at + 1.0/30.0 < clock_now) + { + std::stringstream ss; + sim::graphics::mesh rmesh; + clock_at += 1.0/30.0; - rmesh.load_text(ss.str().c_str(), 0.04); - mesh2.bind(); - mesh2.set(rmesh, GL_DYNAMIC_DRAW); + ss << "\n\n"; + ss << show( sys.turbine_bypass_valve->get_state() * 100 ) << " %\n"; + ss << show( sys.turbine_bypass_valve->get_flow() / 1000 ) << " kg/s\n"; + ss << "\n\n\n"; + ss << show( sys.turbine_inlet_valve->get_state() * 100 ) << " %\n"; + ss << show( sys.turbine_inlet_valve->get_flow() / 1000 ) << " kg/s\n"; + ss << "\n\n\n"; + ss << show( sys.primary_pump->get_power() * 100 ) << " %\n"; + ss << show( sys.primary_pump->get_rpm() ) << " r/min\n"; + ss << show( sys.primary_pump->get_flow_mass() / 1000 ) << " kg/s\n"; + ss << "\n\n\n"; + ss << show( sys.condenser->get_heat() ) << " C\n"; + ss << show( sys.condenser->get_steam() ) << " g\n"; + ss << show( sys.condenser->get_pressure() / 1000 ) << " kPa\n"; + ss << show( sys.condenser->get_level() / 1000 ) << " / " << show( sys.condenser->get_volume() / 1000 ) << " kL\n"; + + rmesh.load_text(ss.str().c_str(), 0.04); + mesh2.bind(); + mesh2.set(rmesh, GL_DYNAMIC_DRAW); + } if(m_joystick_turbine_bypass.check_focus()) focus::set(std::make_unique(sys.turbine_bypass_valve.get())); diff --git a/src/graphics/monitor/primary_loop.hpp b/src/graphics/monitor/primary_loop.hpp index 20d5fbe..1440714 100644 --- a/src/graphics/monitor/primary_loop.hpp +++ b/src/graphics/monitor/primary_loop.hpp @@ -9,6 +9,7 @@ namespace sim::graphics::monitor class primary_loop { sim::graphics::glmesh mesh1, mesh2; + double clock_at = 0, clock_now = 0; sim::graphics::glmesh gm_switch_1; diff --git a/src/graphics/monitor/secondary_loop.cpp b/src/graphics/monitor/secondary_loop.cpp index 0c89ec3..5cdb065 100644 --- a/src/graphics/monitor/secondary_loop.cpp +++ b/src/graphics/monitor/secondary_loop.cpp @@ -78,27 +78,34 @@ void secondary_loop::init() void secondary_loop::update(double dt) { - std::stringstream ss; - sim::graphics::mesh rmesh; system& sys = sim::system::active; + clock_now += dt; - ss << "\n\n\n"; - ss << show( sys.evaporator->get_heat() ) << " C\n"; - ss << show( sys.evaporator->get_steam_output() ) << " g/s\n"; - ss << show( sys.evaporator->get_pressure() ) << " Pa\n"; - ss << show( sys.evaporator->get_level() / 1000 ) << " / " << show( sys.evaporator->get_volume() / 1000 ) << " kL\n"; - ss << "\n\n\n"; - ss << show( sys.secondary_pump->get_power() * 100 ) << " %\n"; - ss << show( sys.secondary_pump->get_rpm() ) << " r/min\n"; - ss << show( sys.secondary_pump->get_flow_mass() / 1000 ) << " kg/s\n"; - ss << "\n\n\n"; - ss << show( sys.freight_pump->get_power() * 100 ) << " %\n"; - ss << show( sys.freight_pump->get_rpm() ) << " r/min\n"; - ss << show( sys.freight_pump->get_flow_mass() / 1000 ) << " kg/s\n"; + if(clock_at + 1.0/30.0 < clock_now) + { + std::stringstream ss; + sim::graphics::mesh rmesh; + clock_at += 1.0/30.0; + + ss << "\n\n\n"; + ss << show( sys.evaporator->get_heat() ) << " C\n"; + ss << show( sys.evaporator->get_steam_output() ) << " g/s\n"; + ss << show( sys.evaporator->get_pressure() ) << " Pa\n"; + ss << show( sys.evaporator->get_level() / 1000 ) << " / " << show( sys.evaporator->get_volume() / 1000 ) << " kL\n"; + ss << "\n\n\n"; + ss << show( sys.secondary_pump->get_power() * 100 ) << " %\n"; + ss << show( sys.secondary_pump->get_rpm() ) << " r/min\n"; + ss << show( sys.secondary_pump->get_flow_mass() / 1000 ) << " kg/s\n"; + ss << "\n\n\n"; + ss << show( sys.freight_pump->get_power() * 100 ) << " %\n"; + ss << show( sys.freight_pump->get_rpm() ) << " r/min\n"; + ss << show( sys.freight_pump->get_flow_mass() / 1000 ) << " kg/s\n"; + + rmesh.load_text(ss.str().c_str(), 0.04); + mesh2.bind(); + mesh2.set(rmesh, GL_DYNAMIC_DRAW); + } - rmesh.load_text(ss.str().c_str(), 0.04); - mesh2.bind(); - mesh2.set(rmesh, GL_DYNAMIC_DRAW); if(m_switch_2.check_focus()) toggle_secondary_pump(); diff --git a/src/graphics/monitor/secondary_loop.hpp b/src/graphics/monitor/secondary_loop.hpp index 647f467..1e08335 100644 --- a/src/graphics/monitor/secondary_loop.hpp +++ b/src/graphics/monitor/secondary_loop.hpp @@ -9,6 +9,7 @@ namespace sim::graphics::monitor class secondary_loop { sim::graphics::glmesh mesh1, mesh2; + double clock_at = 0, clock_now = 0; sim::graphics::glmesh gm_switch_2; sim::graphics::glmesh gm_switch_3; diff --git a/src/graphics/monitor/vessel.cpp b/src/graphics/monitor/vessel.cpp index 491c960..87377a2 100644 --- a/src/graphics/monitor/vessel.cpp +++ b/src/graphics/monitor/vessel.cpp @@ -56,10 +56,17 @@ void vessel::update(double dt) std::stringstream ss; sim::graphics::mesh rmesh; sim::system& sys = sim::system::active; + clock_now += dt; + + if(clock_at + 1.0/30.0 > clock_now) + { + return; + } double temp_min, temp_max; double crod_min = INFINITY, crod_max = -INFINITY; + clock_at += 1.0/30.0; sys.reactor->get_stats(sim::reactor::rod::val_t::HEAT, temp_min, temp_max); for(int i = 0; i < sys.reactor->size; i++) diff --git a/src/graphics/monitor/vessel.hpp b/src/graphics/monitor/vessel.hpp index fce8e4d..99e3012 100644 --- a/src/graphics/monitor/vessel.hpp +++ b/src/graphics/monitor/vessel.hpp @@ -9,6 +9,7 @@ namespace sim::graphics::monitor class vessel { sim::graphics::glmesh mesh1, mesh2; + double clock_at = 0, clock_now = 0; public: diff --git a/src/main.cpp b/src/main.cpp index 4600df8..5ca443c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -36,9 +37,13 @@ int main() // tests::run(); // return 0; + std::ofstream log("log.csv"); graphics::window::create(); long clock = get_now(); + double at = 0; + + log << R"("clock","level","steam","heat","pressure")" << "\n"; while(!graphics::window::should_close()) { @@ -46,6 +51,7 @@ int main() long passed = now - clock; double dt = (double)passed / 1e6; clock += passed; + at += dt * sim::system::active.speed; sim::system::active.update(dt); @@ -53,6 +59,12 @@ int main() graphics::window::update(dt); graphics::focus::update(dt); graphics::window::render(); + + log << at << ","; + log << sim::system::active.condenser->get_level() << ","; + log << sim::system::active.condenser->get_steam() << ","; + log << sim::system::active.condenser->get_heat() << ","; + log << sim::system::active.condenser->get_pressure() << "\n"; } graphics::window::destroy(); diff --git a/src/system.cpp b/src/system.cpp index 8d2e757..7417f9e 100644 --- a/src/system.cpp +++ b/src/system.cpp @@ -84,11 +84,11 @@ void system::update(double dt) turbine->update(dt); condenser->update(dt); -// evaporator->update(dt); -// condenser_secondary->update(dt); + evaporator->update(dt); + condenser_secondary->update(dt); -// primary_pump->update(dt); -// secondary_pump->update(dt); -// freight_pump->update(dt); + primary_pump->update(dt); + secondary_pump->update(dt); + freight_pump->update(dt); } diff --git a/src/tests.cpp b/src/tests.cpp index 03cf3ac..7208c3b 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -22,23 +22,55 @@ std::ostream& operator<<(std::ostream& o, const fluid_holder& fh) void tests::run() { - fluid_holder src(WATER, 1e6, 0); - fluid_holder dst(WATER, 1e6, 0); + fluid_holder fhs[] = { + fluid_holder(WATER, 75398, 0), + fluid_holder(WATER, 75398, 0), + }; - src.add_fluid(1e5, 11); + valve vs[] = { + valve(&fhs[0], &fhs[1], 1, 0.1), + }; - valve v(&src, &dst, 1, 0.5); - double dt = 0.001; + fhs[0].level = 100; + fhs[0].steam = 0; + fhs[0].heat = 100; + + double dt = 0.1; double at = 0; + std::cout << "time"; + + for(fluid_holder& fh : fhs) + { + std::cout << "\t\tlevel (L)\tsteam (g)\theat (C)\tpressure (Pa)"; + } + + std::cout << "\n"; + for(int i = 0; i < 10000; i++) { - src.update(dt); - dst.update(dt); - v.update(dt); + for(fluid_holder& fh : fhs) + { + fh.update(dt); + } + + for(valve& v : vs) + { + v.update(dt); + } - std::cout << at << "\t" << src.get_pressure() << "\t" << dst.get_pressure() << "\n"; + std::cout << at; + + for(const fluid_holder& fh : fhs) + { + std::cout << "\t\t" << fh.get_level() << "\t" << fh.get_steam() << "\t" << fh.get_heat() << "\t" << fh.get_pressure(); + } + + std::cout << "\n"; at += dt; } + + std::cout << "\n" << fhs[0] << "\n"; + }