bug for pressure issues at low volume persists, but i'm doing a

workaround :/

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

src/coolant/condenser_secondary.cpp

@@ -6,8 +6,8 @@
 
 using namespace sim::coolant;
 
-condenser_secondary::condenser_secondary(condenser* primary, evaporator* source) :
+condenser_secondary::condenser_secondary(condenser* primary, evaporator* source, double volume) :
-		primary(primary), source(source), fluid_holder(primary->fluid, 0, 0)
+		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);
+}
+

src/coolant/condenser_secondary.hpp

@@ -15,7 +15,7 @@ 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); }
@@ -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);
 };
 
 };

src/coolant/evaporator.cpp

@@ -28,12 +28,6 @@ double evaporator::get_steam_output()
 
 void evaporator::update(double dt)
 {
-	((fluid_holder*)this)->update(dt);
+	update_base(dt);
-/*
-	double m = std::pow(0.5, dt / 3600);
-	double steam_out = steam * (1 - m);
-	steam *= m;
-
-	steam_output = steam_out / dt;*/
 }
 

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);
+	return (pressure * V) / (constants::R * heat);
-}
-
-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()));
+	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,33 +111,30 @@ 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 K = conversions::temperature::c_to_k(heat);	// K
-		double heat_k = conversions::temperature::c_to_k(heat);
+		double P = fluid.vapor_pressure.calc_p(K);			// Pa
-		double target_pressure = fluid.vapor_pressure.calc_p(heat);
+		double R = sim::constants::R;						// J/K/mol
 
-		double K = heat_k;
+		double J_m = fluid.jPg * fluid.gPmol;	// J/mol
-		double R = 1000 * constants::R;
+		double n_g = fluid.g_to_mol(steam);		// mol
-		double P = target_pressure;
+		double V_g = (volume - level) * 0.001;	// m^3
-		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 n = (-K*R*n_g - L_m*P*n_l + P*V_t)/(K*R - L_m*P) * 0.5;
+		double n = (P * V_g) / (R * K) - n_g;	// mol
-
+		double l = level - fluid.mol_to_l(n);	// L
-		/*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);
 
 		if(l < 0)
 		{
@@ -149,8 +142,16 @@ void fluid_holder::update(double secs)
 			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;
 	}
 }

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(double heat, double pressure, double mol);
-	static double calc_pressure_mol(double temp, double volume, double pressure);
+	static double calc_pressure_mol(double heat, double pressure, double volume);
-	static double calc_pressure_vol(double heat, double pressure, double mol);
 	
-	void update(double dt);
+	void update_base(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, 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});
 
 }
 

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

src/coolant/pump.hpp

@@ -36,6 +36,7 @@ public:
 	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);

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

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;
+	return t > -C ? std::pow(10, A - B / (C + t)) : 0;
-   
-	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
 {
-	double P = pressure::pa_to_mmhg(p);
+	return B / (A - std::log(p) / std::log(10)) - C;
-
-	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,12 +6,9 @@ namespace sim::coolant
 
 struct vapor_pressure
 {
-	const double A1, B1, C1;
+	const double A, B, C;
-	const double A2, B2, C2;
-	const double T;
 	
-	constexpr vapor_pressure(double A1, double B1, double C1, double T, double A2, double B2, double C2) :
+	constexpr vapor_pressure(double A, double B, double C) : A(A), B(B), C(C) { }
-			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/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";

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);
+	update_base(secs);
 	
-	double diff = steam - s;
+	double diff = steam - steam_last;
 	double hl = get_bubble_hl();
 
-	if(diff > 0)
+	steam_last = steam;
-	{
+	steam_suspended += diff;
-		steam_suspended += diff;
-	}
 
 	if(hl > 0)
 	{

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);
+	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, 10, coolant::pump::mode_t::NONE, 0);
+	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);
 }
 
@@ -75,20 +75,17 @@ 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);
 }
 

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_system fs(WATER, 1000, 10);
-		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";
 
+//	std::cout << "Volume: " << fs.volume << "\n";
 }
 

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