added evaporative heat loss

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

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;

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;
+	double s = steam + n;
-	level -= fluid.g_to_l(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);
+		s += l;
-		level = 0;
+		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 m1 = 1000000;
-	double m2 = level + fluid.g_to_l(steam);
+	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);

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