alternate (but equally flawed) thermodynamic equalibrium function

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)
 

nasty_function.c++

@@ -0,0 +1,51 @@
+
+#include <cmath>
+
+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;
+}
+

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; }
 
 };
 

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);
 }
 
 

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;
 	}
 }
 

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);
 };

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});
 
 }
 

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;
 }
 

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();
 

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;
+	}
 }
 

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;

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<valve_joystick>(sys.turbine_bypass_valve.get()));

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;
 	

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();

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;

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++)

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:
 

src/main.cpp

@@ -3,6 +3,7 @@
 
 #include <random>
 #include <sstream>
+#include <fstream>
 #include <cmath>
 #include <cfenv>
 
@@ -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();

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);
 }
 

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";
+
 }