diff options
Diffstat (limited to 'wi.c')
| -rw-r--r-- | wi.c | 117 |
1 files changed, 117 insertions, 0 deletions
@@ -0,0 +1,117 @@ +#include <stdio.h> +#include <math.h> +#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_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)] */ + +#define L_W 2501.0 /* enthalpy of vaporisation of water at T=273K [kJ/kg] */ + +#define A_W (M_W/M_AIR) +#define M_W 18.0153 /* molar mass of water [g/mol] */ +#define M_AIR 28.9645 /* molar mass of air [g/mol] */ + +#define ABS_ERR_BOUND 1e-7 /* absolute error bound for root solver */ +#define MAXITER 1000 /* maximum iterations for root solver */ + +#define T_AMBIENT 298.0 /* ambient temperature [K] */ + +struct theta_wb_params { + double h1; + double p2; +}; + +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); +double eq_specific_water_content(double p, double t); + +int +main(int argc, char *argv[]) { + double t1 = 170.0 + 273.15; /* K */ + double w1 = 0.0; /* dry air */ + double h1 = mixture_specific_enthalpy(t1, w1); + double p2 = 2e5; /* Pa */ + printf("t1 = %f *C\n", t1-273.15); + printf("p2 = %f\n", p2); + double t_wb = wet_bulb_temp(h1, p2); + 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); +} + +/* specific enthalpy of mixture h [kJ/kg] at temperature t [K] and specific water content w */ +double +mixture_specific_enthalpy(double t, double w) { + return (CP_DFM + w*CP_VAP)*t + w*L_W; +} + +double +wet_bulb_temp(double h1, double p2) { + struct theta_wb_params params = {h1, p2}; + gsl_multiroot_function f = {&theta_wb, 1, ¶ms}; + + double x_init = T_AMBIENT; + gsl_vector *x = gsl_vector_alloc(1); + gsl_vector_set(x, 0, x_init); + + const gsl_multiroot_fsolver_type *t = gsl_multiroot_fsolver_dnewton; + gsl_multiroot_fsolver *s = gsl_multiroot_fsolver_alloc(t, 1); + gsl_multiroot_fsolver_set(s, &f, x); + + int status; + size_t iter = 0; + do { + iter++; + status = gsl_multiroot_fsolver_iterate(s); + if (status) { + break; + } + status = gsl_multiroot_test_residual(s->f, ABS_ERR_BOUND); + } while (status == GSL_CONTINUE && iter < MAXITER); + + double res = gsl_vector_get(s->x, 0); + gsl_multiroot_fsolver_free(s); + gsl_vector_free(x); + return res; +} + +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; + 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; +} + +/* equilibrium specific water content at temperature t [K] and pressure p [Pa] */ +double +eq_specific_water_content(double p, double t) { + double p_eq = eq_vapor_pressure(t); + return A_W * p_eq/(p - p_eq); +} + +/* equilibrium vapor pressure [Pa] at temperature t [K] according to Wexler 1976 */ +double +eq_vapor_pressure(double t) { + return 1.0 + / exp(2.9912729e3 / pow(t, 2)) + / exp(6.0170128e3 / t) + * exp(1.887643845e1) + / exp(2.8354721e-2 * t) + * exp(1.7838301e-5 * pow(t, 2)) + / exp(8.4150417e-10 * pow(t, 3)) + * exp(4.4412543e-13 * pow(t, 4)) + * exp(2.858487 * log(t)); +} + |