From 1264e623373460002b4902cfc44f9b23544f1de6 Mon Sep 17 00:00:00 2001
From: Jay Robson <jsrobson10@gmail.com>
Date: Tue, 13 Feb 2024 15:28:39 +1100
Subject: [PATCH] bug for pressure issues at low volume persists, but i'm doing
 a workaround :/

---
 src/coolant/condenser.cpp               |  2 +-
 src/coolant/condenser_secondary.cpp     |  9 +++-
 src/coolant/condenser_secondary.hpp     |  6 ++-
 src/coolant/evaporator.cpp              |  8 +--
 src/coolant/fluid_holder.cpp            | 69 +++++++++++++------------
 src/coolant/fluid_holder.hpp            | 11 ++--
 src/coolant/fluid_t.hpp                 |  2 +-
 src/coolant/pump.cpp                    |  2 +-
 src/coolant/pump.hpp                    |  3 +-
 src/coolant/valve.cpp                   |  4 +-
 src/coolant/vapor_pressure.cpp          | 29 +----------
 src/coolant/vapor_pressure.hpp          |  7 +--
 src/graphics/monitor/secondary_loop.cpp |  4 +-
 src/reactor/coolant/vessel.cpp          | 14 +++--
 src/system.cpp                          | 27 +++++-----
 src/tests.cpp                           | 63 +---------------------
 src/util/constants.hpp                  |  4 +-
 17 files changed, 89 insertions(+), 175 deletions(-)

diff --git a/src/coolant/condenser.cpp b/src/coolant/condenser.cpp
index c6d1cac..2ca882e 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)
 {
-	((fluid_holder*)this)->update(secs);
+	update_base(secs);
 }
 
 
diff --git a/src/coolant/condenser_secondary.cpp b/src/coolant/condenser_secondary.cpp
index 1f4042b..1908593 100644
--- a/src/coolant/condenser_secondary.cpp
+++ b/src/coolant/condenser_secondary.cpp
@@ -6,8 +6,8 @@
 
 using namespace sim::coolant;
 
-condenser_secondary::condenser_secondary(condenser* primary, evaporator* source) :
-		primary(primary), source(source), fluid_holder(primary->fluid, 0, 0)
+condenser_secondary::condenser_secondary(condenser* primary, evaporator* source, double volume) :
+		primary(primary), source(source), fluid_holder(primary->fluid, volume, 0)
 {
 	
 }
@@ -18,3 +18,8 @@ double condenser_secondary::add_fluid(double amount, double heat)
 	return source->add_fluid(amount, heat);
 }
 
+void condenser_secondary::update(double dt)
+{
+	heat = primary->add_heat(fluid.l_to_g(level), heat);
+}
+
diff --git a/src/coolant/condenser_secondary.hpp b/src/coolant/condenser_secondary.hpp
index 72fad7c..2976d04 100644
--- a/src/coolant/condenser_secondary.hpp
+++ b/src/coolant/condenser_secondary.hpp
@@ -15,14 +15,14 @@ class condenser_secondary : public fluid_holder
 	
 public:
 
-	condenser_secondary(condenser* primary, evaporator* source);
+	condenser_secondary(condenser* primary, evaporator* source, double volume);
 
 	virtual double add_heat(double m, double t) { return source->add_heat(m, t); }
 	virtual void add_steam(double amount, double t) { return source->add_steam(amount, t); }
 	virtual double extract_fluid(double amount) { return source->extract_fluid(amount); }
 
 	virtual double add_fluid(double amount, double heat);
-
+	
 	virtual double get_volume() const { return source->get_volume(); }
 	virtual double get_level() const { return source->get_level(); }
 	virtual double get_heat() const { return source->get_heat(); } // celsius
@@ -32,6 +32,8 @@ public:
 	virtual double get_thermal_mass() const { return source->get_thermal_mass(); } // grams
 	virtual double get_pressure() const { return source->get_pressure(); } // pascals
 	virtual double get_steam_density() const { return source->get_steam_density(); } // g/L
+	
+	void update(double dt);
 };
 
 };
diff --git a/src/coolant/evaporator.cpp b/src/coolant/evaporator.cpp
index de72c7a..50f9ac3 100644
--- a/src/coolant/evaporator.cpp
+++ b/src/coolant/evaporator.cpp
@@ -28,12 +28,6 @@ double evaporator::get_steam_output()
 
 void evaporator::update(double dt)
 {
-	((fluid_holder*)this)->update(dt);
-/*
-	double m = std::pow(0.5, dt / 3600);
-	double steam_out = steam * (1 - m);
-	steam *= m;
-
-	steam_output = steam_out / dt;*/
+	update_base(dt);
 }
 
diff --git a/src/coolant/fluid_holder.cpp b/src/coolant/fluid_holder.cpp
index c67d91d..00cbed5 100644
--- a/src/coolant/fluid_holder.cpp
+++ b/src/coolant/fluid_holder.cpp
@@ -34,6 +34,11 @@ double fluid_holder::add_heat(double m1, double t1)
 
 double fluid_holder::add_fluid(double v2, double t2)
 {
+	if(level + v2 <= 0)
+	{
+		return 0;
+	}
+	
 	if(level + v2 > volume - 1e-3)
 	{
 		v2 = volume - level - 1e-3;
@@ -83,30 +88,21 @@ void fluid_holder::add_steam(double m2, double t2)
 
 double fluid_holder::calc_pressure(double heat, double volume, double mol)
 {
-	double T = conversions::temperature::c_to_k(heat);
 	double V = volume * 0.001;
 
-	return V == 0 ? 0 : (mol * T * constants::R) / V;
+	return V == 0 ? 0 : (mol * heat * constants::R) / V;
 }
 
 double fluid_holder::calc_pressure_mol(double heat, double volume, double pressure)
 {
-	double T = conversions::temperature::c_to_k(heat);
 	double V = volume * 0.001;
 
-	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;
+	return (pressure * V) / (constants::R * heat);
 }
 
 double fluid_holder::get_pressure() const
 {
-	return calc_pressure(heat, get_steam_volume(), fluid.g_to_mol(get_steam()));
+	return calc_pressure(conversions::temperature::c_to_k(heat), get_steam_volume(), fluid.g_to_mol(steam));
 }
 
 double fluid_holder::get_steam_density() const
@@ -115,42 +111,47 @@ double fluid_holder::get_steam_density() const
 	return v > 0 ? steam / v : 0;
 }
 
-void fluid_holder::update(double secs)
+constexpr double calc_extra_steam(double K, double P, double L_m, double J_m, double n_g, double n_l, double V_t)
+{
+	double R = sim::constants::R * 1000;
+	double n = (P * (V_t - n_l * L_m)) / (R * K) - n_g;
+
+	return n;
+}
+
+void fluid_holder::update_base(double secs)
 {
 	double mass = get_thermal_mass();
 
 	if(mass > 0)
 	{
-		// use ideal gas law to get target steam pressure
-		double heat_k = conversions::temperature::c_to_k(heat);
-		double target_pressure = fluid.vapor_pressure.calc_p(heat);
+		double K = conversions::temperature::c_to_k(heat);	// K
+		double P = fluid.vapor_pressure.calc_p(K);			// Pa
+		double R = sim::constants::R;						// J/K/mol
 
-		double K = heat_k;
-		double R = 1000 * constants::R;
-		double P = target_pressure;
-		double J_m = fluid.jPg * fluid.gPmol;
-		double L_m = fluid.gPmol / fluid.gPl;
-		double n_g = fluid.g_to_mol(steam);
-		double n_l = fluid.l_to_mol(level);
-		double V_t = volume;
+		double J_m = fluid.jPg * fluid.gPmol;	// J/mol
+		double n_g = fluid.g_to_mol(steam);		// mol
+		double V_g = (volume - level) * 0.001;	// m^3
 
-		double n = (-K*R*n_g - L_m*P*n_l + P*V_t)/(K*R - L_m*P) * 0.5;
-
-		/*if(std::abs(n.imag()) > std::numeric_limits<double>::epsilon())
-		{
-			throw std::runtime_error("Nonzero imaginary component");
-		}*/
-
-		double l = level - fluid.mol_to_l(n);
+		double n = (P * V_g) / (R * K) - n_g;	// mol
+		double l = level - fluid.mol_to_l(n);	// L
 
 		if(l < 0)
 		{
 			n -= fluid.l_to_mol(l);
 			l = 0;
 		}
-
-		level = l;
+		
 		steam += fluid.mol_to_g(n);
+
+		if(steam < 0)
+		{
+			l -= fluid.g_to_l(steam);
+			n += fluid.g_to_mol(steam);
+			steam = 0;
+		}
+		
+		level = l;
 		heat -= n * J_m / mass;
 	}
 }
diff --git a/src/coolant/fluid_holder.hpp b/src/coolant/fluid_holder.hpp
index a132544..9dfc432 100644
--- a/src/coolant/fluid_holder.hpp
+++ b/src/coolant/fluid_holder.hpp
@@ -2,6 +2,7 @@
 #pragma once
 
 #include "fluid_t.hpp"
+#include "../conversions/temperature.hpp"
 
 namespace sim::coolant
 {
@@ -31,6 +32,7 @@ public:
 	virtual double get_volume() const { return volume; } // litres
 	virtual double get_level() const { return level; } // litres
 	virtual double get_heat() const { return heat; } // celsius
+	virtual double get_heat_k() const { return conversions::temperature::c_to_k(get_heat()); } // kelvin
 	virtual double get_steam() const { return steam; } // grams
 	virtual double get_steam_volume() const { return get_volume() - get_level(); } // litres
 	virtual double get_mass() const { return fluid.l_to_g(get_level()) + get_steam(); } // grams
@@ -38,11 +40,10 @@ public:
 	virtual double get_pressure() const; // pascals
 	virtual double get_steam_density() const; // g/L
 	
-	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);
+	static double calc_pressure(double heat, double pressure, double mol);
+	static double calc_pressure_mol(double heat, double pressure, double volume);
+	
+	void update_base(double dt);
 };
 
 };
diff --git a/src/coolant/fluid_t.hpp b/src/coolant/fluid_t.hpp
index ed3bdc8..7619532 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, 108.266, 8.14019, 1810.94, 244.485});
+constexpr const fluid_t WATER = fluid_t(1000, 18, 2257, 4.1816, {8.07131 + 2.124903, 1730.63, 233.426 - 273.15});
 
 }
 
diff --git a/src/coolant/pump.cpp b/src/coolant/pump.cpp
index c951d90..8a2644c 100644
--- a/src/coolant/pump.cpp
+++ b/src/coolant/pump.cpp
@@ -92,7 +92,7 @@ void pump::update(double dt)
 	double src_heat = src->get_heat();
 	double p_diff_1 = dst->get_pressure() - src->get_pressure();
 
-	double max_volume = std::min(src->get_level(), dst->get_level());
+	double max_volume = ignore_dst_level ? src->get_level() : std::min(src->get_level(), dst->get_volume() - dst->get_level());
 	double src_volume = src->extract_fluid(std::min(get_flow_target() * dt, max_volume));
 	double dst_volume = dst->add_fluid(src_volume, src_heat);
 
diff --git a/src/coolant/pump.hpp b/src/coolant/pump.hpp
index 68bc16d..3861a46 100644
--- a/src/coolant/pump.hpp
+++ b/src/coolant/pump.hpp
@@ -35,7 +35,8 @@ public:
 	const double friction; // J/rev
 	const double max_power; // W
 	const double target; // L
-
+	
+	bool ignore_dst_level = false;
 	bool powered = false;
 
 	pump(fluid_holder* src, fluid_holder* dst, double mass, double radius, double power, double l_per_rev, double friction, mode_t mode, double target);
diff --git a/src/coolant/valve.cpp b/src/coolant/valve.cpp
index 2a4f1b8..3e553f6 100644
--- a/src/coolant/valve.cpp
+++ b/src/coolant/valve.cpp
@@ -47,13 +47,13 @@ void valve::update(double dt)
 
 	if(remove < 0)
 	{
-		mol = fluid_holder::calc_pressure_mol(src->get_heat(), src->get_steam_volume(), pressure1 - remove);
+		mol = fluid_holder::calc_pressure_mol(src->get_heat_k(), 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);
+		mol = fluid_holder::calc_pressure_mol(dst->get_heat_k(), dst->get_steam_volume(), pressure2 - remove);
 		mass = dst->get_steam() - dst->fluid.mol_to_g(mol);
 	}
 
diff --git a/src/coolant/vapor_pressure.cpp b/src/coolant/vapor_pressure.cpp
index 4089ac2..4994ef0 100644
--- a/src/coolant/vapor_pressure.cpp
+++ b/src/coolant/vapor_pressure.cpp
@@ -1,42 +1,17 @@
 
 #include "vapor_pressure.hpp"
-#include "../conversions/temperature.hpp"
-#include "../conversions/pressure.hpp"
 
 #include <cmath>
 
 using namespace sim::coolant;
-using namespace sim::conversions;
 
 double vapor_pressure::calc_p(double t) const
 {
-	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);
+	return t > -C ? std::pow(10, A - B / (C + t)) : 0;
 }
 
 double vapor_pressure::calc_t(double p) const
 {
-	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;
-	}
+	return B / (A - std::log(p) / std::log(10)) - C;
 }
 
diff --git a/src/coolant/vapor_pressure.hpp b/src/coolant/vapor_pressure.hpp
index 464e3ae..46ceb86 100644
--- a/src/coolant/vapor_pressure.hpp
+++ b/src/coolant/vapor_pressure.hpp
@@ -6,12 +6,9 @@ namespace sim::coolant
 
 struct vapor_pressure
 {
-	const double A1, B1, C1;
-	const double A2, B2, C2;
-	const double T;
+	const double A, B, 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) { }
+	constexpr vapor_pressure(double A, double B, double C) : A(A), B(B), C(C) { }
 
 	double calc_p(double t) const;
 	double calc_t(double p) const;
diff --git a/src/graphics/monitor/secondary_loop.cpp b/src/graphics/monitor/secondary_loop.cpp
index 5cdb065..c47682d 100644
--- a/src/graphics/monitor/secondary_loop.cpp
+++ b/src/graphics/monitor/secondary_loop.cpp
@@ -87,10 +87,10 @@ void secondary_loop::update(double dt)
 		sim::graphics::mesh rmesh;
 		clock_at += 1.0/30.0;
 
-		ss << "\n\n\n";
+		ss << "\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_pressure() / 1000 ) << " kPa\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";
diff --git a/src/reactor/coolant/vessel.cpp b/src/reactor/coolant/vessel.cpp
index 78db76a..87f2649 100644
--- a/src/reactor/coolant/vessel.cpp
+++ b/src/reactor/coolant/vessel.cpp
@@ -55,17 +55,15 @@ double vessel::get_bubble_hl() const
 
 void vessel::update(double secs)
 {
-	double s = steam;
+	double steam_last = steam;
 
-	((sim::coolant::fluid_holder*)this)->update(secs);
-
-	double diff = steam - s;
+	update_base(secs);
+	
+	double diff = steam - steam_last;
 	double hl = get_bubble_hl();
 
-	if(diff > 0)
-	{
-		steam_suspended += diff;
-	}
+	steam_last = steam;
+	steam_suspended += diff;
 
 	if(hl > 0)
 	{
diff --git a/src/system.cpp b/src/system.cpp
index 7417f9e..e963766 100644
--- a/src/system.cpp
+++ b/src/system.cpp
@@ -37,18 +37,18 @@ system::system()
 	
 	vessel = std::make_unique<reactor::coolant::vessel>(sim::coolant::WATER, 8, 10, 6e6, 5e5);
 	reactor = std::make_unique<reactor::reactor>(sim::reactor::builder(19, 19, 1.0 / 4.0, 4, reactor::fuel::fuel_rod(0.5), vessel.get(), layout));
-	condenser = std::make_unique<coolant::condenser>(sim::coolant::WATER, 6, 4, 3e6, 0);
+	condenser = std::make_unique<coolant::condenser>(sim::coolant::WATER, 6, 4, 3e6, 30000);
 	turbine = std::make_unique<electric::turbine>(sim::coolant::WATER, condenser.get(), 6, 3, 2e6);
 	
 	sink = std::make_unique<coolant::sink>(sim::coolant::WATER, 11, 0, 0);
 	evaporator = std::make_unique<coolant::evaporator>(sim::coolant::WATER, 2, 30, 0, 1000);
-	condenser_secondary = std::make_unique<coolant::condenser_secondary>(condenser.get(), evaporator.get());
+	condenser_secondary = std::make_unique<coolant::condenser_secondary>(condenser.get(), evaporator.get(), 1000);
 	
 	turbine_inlet_valve = std::make_unique<coolant::valve>(vessel.get(), turbine.get(), 0, 0.5);
 	turbine_bypass_valve = std::make_unique<coolant::valve>(vessel.get(), condenser.get(), 0, 0.5);
 
-	primary_pump = std::make_unique<coolant::pump>(condenser.get(), vessel.get(), 1e5, 1, 1e4, 0.1, 10, coolant::pump::mode_t::SRC, 35000);
-	secondary_pump = std::make_unique<coolant::pump>(evaporator.get(), condenser_secondary.get(), 1e5, 1, 1e4, 0.1, 10, coolant::pump::mode_t::NONE, 0);
+	primary_pump = std::make_unique<coolant::pump>(condenser.get(), vessel.get(), 1e5, 1, 1e5, 0.1, 10, coolant::pump::mode_t::SRC, 35000);
+	secondary_pump = std::make_unique<coolant::pump>(evaporator.get(), condenser_secondary.get(), 1e5, 1, 1e4, 0.1, 1, coolant::pump::mode_t::NONE, 0);
 	freight_pump = std::make_unique<coolant::pump>(sink.get(), evaporator.get(), 1e5, 1, 1e4, 0.1, 10, coolant::pump::mode_t::DST, 1e6);
 }
 
@@ -74,21 +74,18 @@ system::system(system&& o)
 void system::update(double dt)
 {
 	dt *= speed;
-
-	reactor->update(dt);
-	vessel->update(dt);
-
+	
 	turbine_inlet_valve->update(dt);
 	turbine_bypass_valve->update(dt);
-
-	turbine->update(dt);
-	condenser->update(dt);
-
-	evaporator->update(dt);
-	condenser_secondary->update(dt);
-
 	primary_pump->update(dt);
 	secondary_pump->update(dt);
 	freight_pump->update(dt);
+	reactor->update(dt);
+
+	vessel->update(dt);
+	turbine->update(dt);
+	condenser->update(dt);
+	evaporator->update(dt);
+	condenser_secondary->update(dt);
 }
 
diff --git a/src/tests.cpp b/src/tests.cpp
index 7c53c6f..8ba9626 100644
--- a/src/tests.cpp
+++ b/src/tests.cpp
@@ -1,74 +1,15 @@
 
 #include "tests.hpp"
-#include "coolant/valve.hpp"
-#include "coolant/fluid_holder.hpp"
 
 #include <unistd.h>
 #include <iostream>
 
 using namespace sim;
-using namespace sim::coolant;
-
-std::ostream& operator<<(std::ostream& o, const fluid_holder& fh)
-{
-	o << "Fluid Holder\n";
-	o << "Heat          " << fh.get_heat() << " C\n";
-	o << "Steam         " << fh.get_steam() << " g\n";
-	o << "Pressure      " << fh.get_pressure() << " Pa\n";
-	o << "Volume        " << fh.get_level() / 1000 << " / " << fh.get_volume() / 1000 << " kL\n\n";
-
-	return o;
-}
 
 void tests::run()
 {
-	fluid_holder fhs[] = {
-		fluid_holder(WATER, 75398, 0),
-	};
-
-	valve vs[] = {
-	};
-
-	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++)
-	{
-		for(fluid_holder& fh : fhs)
-		{
-			fh.update(dt);
-		}
-		
-		for(valve& v : vs)
-		{
-			v.update(dt);
-		}
-
-		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";
+//	fluid_system fs(WATER, 1000, 10);
 
+//	std::cout << "Volume: " << fs.volume << "\n";
 }
 
diff --git a/src/util/constants.hpp b/src/util/constants.hpp
index 8539250..975c669 100644
--- a/src/util/constants.hpp
+++ b/src/util/constants.hpp
@@ -3,6 +3,8 @@
 
 namespace sim::constants
 {
-	constexpr const double R = 8.31446261815324; // molar gas constant, J/Kmol
+	constexpr double R = 8.31446261815324; // molar gas constant, J/mol/K
+	constexpr double R_air = 0.2870500676; // specific gas constant of dry air, J/g/K
+	constexpr double M_air = 28.9652; // molar mass of dry air, g/mol
 };