duty cycle calculations

1 files changed, 102 insertions, 12 deletions

diff --git a/wi.c b/wi.c
index 1b55c98..af3df8d 100644
--- a/wi.c
+++ b/wi.c
@@ -3,7 +3,8 @@
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_multiroots.h>
 
-#define CP_DFM CP_AIR /* specific heat of dry mixture at constant pressure at 273K [kJ/(kg*K)] --- CP_DFM=CP_AIR assumes no fuel */
+#define CP_DFM CP_AIR /* specific heat of dry mixture at constant pressure at 273K [kJ/(kg*K)]
+						 --- CP_DFM=CP_AIR assumes no fuel */
 #define CP_AIR 1.006 /* specific heat of dry air at constant pressure at 273K [kJ/(kg*K)] */
 #define CP_VAP 1.805 /* specific heat of water vapor at constant pressure at 273K [kJ/(kg*K)] */
 
@@ -18,17 +19,57 @@
 
 #define T_AMBIENT 298.0 /* ambient temperature [K] */
 
-struct theta_wb_params {
-	double h1;
-	double p2;
+#define N_P_BRPOINTS 11
+#define N_RPM_BRPOINTS 8
+#define V_DISP 2000e-3 /* displaced volume [m^3] */
+#define T_REF 323.15 /* reference air temperature [K] */
+#define ROH_W_REF  997.0 /* reference water density at T=25*C [kg/m^3] */
+#define P_W_REF 689475.7 /* reference water pressure [Pa] */
+#define V_RATE_MAX_W_REF 0.000066666666666668 /* maximum water flow rate at P_W_REF [m^3/s] */
+
+/* types */
+
+/* web bulb temperature function parameters */
+struct t_wb_params {
+	double h1; /* specific enthalpy at point 1 */
+	double p2; /* absolute pressure at point 2 */
 };
 
+/* function declarations */
+unsigned int duty_cycle(double p, double t, unsigned int s);
+double m_rate_w(double p, double t, unsigned int s);
+double m_rate_air(double p, double t, unsigned int s);
+double ve(double p, unsigned int s);
 double mixture_specific_enthalpy(double t, double w);
 double wet_bulb_temp(double h1, double p2);
 double eq_vapor_pressure(double t);
-int theta_wb(const gsl_vector *x, void *params, gsl_vector *f);
+int t_wb(const gsl_vector *x, void *params, gsl_vector *f);
 double eq_specific_water_content(double p, double t);
 
+/* global constants */
+const double p_brpoints[N_P_BRPOINTS] = {500e2, 750e2, 1000e2, 1250e2, 1500e2, 1750e2, 2000e2, 2250e2, 2500e2, 2750e2, 3000e2};
+const unsigned int rpm_brpoints[N_RPM_BRPOINTS] = {1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000};
+const unsigned int ve_tbl[N_P_BRPOINTS][N_RPM_BRPOINTS] = {
+	{75, 80, 85, 90, 95, 95, 93, 90},
+	{75, 80, 85, 90, 95, 95, 93, 90},
+	{75, 80, 85, 90, 95, 95, 93, 90},
+	{75, 80, 85, 90, 95, 95, 93, 90},
+	{75, 80, 85, 90, 95, 95, 93, 90}, /*  -  */
+	{75, 80, 85, 90, 95, 95, 93, 90}, /*  P  */
+	{75, 80, 85, 90, 95, 95, 93, 90}, /*  +  */
+	{75, 80, 85, 90, 95, 95, 93, 90},
+	{75, 80, 85, 90, 95, 95, 93, 90},
+	{75, 80, 85, 90, 95, 95, 93, 90},
+	{75, 80, 85, 90, 95, 95, 93, 90}
+	/*           - rpm +          */
+};
+
+
+/* global variables */
+unsigned int dc_tbl[N_P_BRPOINTS][N_RPM_BRPOINTS]; /* duty cycle table to be populated */
+
+/* function definitions */
+
 int
 main(int argc, char *argv[]) {
 	double t1 = 170.0 + 273.15; /* K */
@@ -41,6 +82,55 @@ main(int argc, char *argv[]) {
 	printf("h1: %f\nwet bulb temp: %f *C\n", h1, t_wb-273.15);
 	double w_eq = eq_specific_water_content(p2, t_wb);
 	printf("w_eq = %f\n", w_eq);
+
+
+	printf("\nDuty cycle table:\n");
+	double p;
+	unsigned int i, j, s;
+	for (i = 0; i < N_P_BRPOINTS; i++) {
+		printf("%4.0f ", p_brpoints[i]*1e-2);
+		for (j = 0; j < N_RPM_BRPOINTS; j++) {
+			p = p_brpoints[i];
+			s = rpm_brpoints[j];
+			printf("%4d ", duty_cycle(p, T_REF, s));
+		}
+		putchar('\n');
+	}
+	printf("%4s ", "");
+	for (j = 0; j < N_RPM_BRPOINTS; j++) {
+		printf("%4d ", rpm_brpoints[j]);
+	}
+	putchar('\n');
+}
+
+/* duty cycle (0-100) at air pressure p [Pa], air temperature t [K], and engine speed s [rpm] */
+unsigned int
+duty_cycle(double p, double t, unsigned int s) {
+	return m_rate_w(p, t, s) / ROH_W_REF / V_RATE_MAX_W_REF;
+}
+
+/* mass flow rate of water at air pressure p [Pa], air temperature t [K], and engine speed s [rpm] */
+double
+m_rate_w(double p, double t, unsigned int s) {
+	double h1, t_wb, w_eq;
+
+	h1 = mixture_specific_enthalpy(t, 0.0);
+	t_wb = wet_bulb_temp(h1, p);
+	w_eq = eq_specific_water_content(p, t_wb);
+	return w_eq * m_rate_air(p, t, s);
+}
+
+/* mass flow rate of air at pressure p [Pa], temperature t[K], and engine speed s [rpm] */
+double
+m_rate_air(double p, double t, unsigned int s) {
+	return 3.483e-3 * p * V_DISP * ve(p, s) * (double) s / t / 120.0;
+}
+
+/* volumetric efficiency at air pressure p [Pa] and engine speed s [rpm] */
+double
+ve(double p, unsigned int s) {
+	/* TODO */
+	return 1.0;
 }
 
 /* specific enthalpy of mixture h [kJ/kg] at temperature t [K] and specific water content w */
@@ -49,10 +139,11 @@ mixture_specific_enthalpy(double t, double w) {
 	return (CP_DFM + w*CP_VAP)*t + w*L_W;
 }
 
+/* wet bulb temperature [K] given specific mixture enthalpy h1 [kJ/kg] and air pressure p2 [Pa] */
 double
 wet_bulb_temp(double h1, double p2) {
-	struct theta_wb_params params = {h1, p2};
-	gsl_multiroot_function f = {&theta_wb, 1, &params};
+	struct t_wb_params params = {h1, p2};
+	gsl_multiroot_function f = {&t_wb, 1, &params};
 
 	double x_init = T_AMBIENT;
 	gsl_vector *x = gsl_vector_alloc(1);
@@ -79,16 +170,15 @@ wet_bulb_temp(double h1, double p2) {
 	return res;
 }
 
+/* wet bulb temperature function to be solved by GSL */
 int
-theta_wb(const gsl_vector *x, void *params, gsl_vector *f) {
-	double h1 = ((struct theta_wb_params *) params)->h1;
-	double p2 = ((struct theta_wb_params *) params)->p2;
+t_wb(const gsl_vector *x, void *params, gsl_vector *f) {
+	double h1 = ((struct t_wb_params *) params)->h1;
+	double p2 = ((struct t_wb_params *) params)->p2;
 	double t = gsl_vector_get(x, 0);
 
 	double w_eq = eq_specific_water_content(p2, t);
-
 	double y = CP_DFM*t + w_eq*(CP_VAP*t + L_W) - h1;
-
 	gsl_vector_set(f, 0, y);
 
 	return GSL_SUCCESS;