Code:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define W_DENS 0.998 //Water density at 20 degrees C
#define E_DENS 0.789 //Ethanol density at 20 degrees C
#define MM_W 18.0153 //Molar mass of water - grams per mol. Source NIST
#define MM_E 46.0684 //Molar mass of ethanot grams per mol. Source NIST
#define MAX_NUM 81
#define GAL2L 3.785411784 //Gallons to liters factor
/*********FUNCTION PROTOYPES*************************/
double get_input(void);
double abv_to_abw(double abv); /**convert alcohol by volume to alcohol by weight*/
double abw_to_mf_e(double abw); /**convert alcohol by weight to mol fraction*/
double mf_to_abw(double mf_e); /**convert mol fraction to alcohol by weight*/
double abw_to_abv(double abw); /**convert alcohol by weight to alcohol by volume*/
double calc_total_mol(double wash_l, double abv);
/**calc mol fraction converts from abv to abw then from abw to ethanol molar fraction*/
double calc_mol_fract(double abv, double (*abv_to_abw)(double abv), double (*abw_to_mf_e)(double abw));
/**van Laar activity coefficient for ethanol and water to model for non-ideal mix*/
double calc_gamma_e(double x1, double x2);
double calc_gamma_w(double x1, double x2);
/**Find saturated vapor pressure*/
double calc_psat(double A, double B, double C, double bp, double gamma, double x);
double calc_boil_pt(double abv, double atm_p, double mol_fract,
double (*calc_gamma_e)(double x1, double x2),
double (*calc_gamma_w)(double x1, double x2),
double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x));
/**calc xw finds boiler mol fraction given boil point and Antoines Eq and Raoults law.*/
double calc_xw(double atm_p, double mol_fract, double boil_pt,
double (*calc_gamma_e)(double x1, double x2),
double (*calc_gamma_w)(double x1, double x2),
double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x));
/**calc yd finds distllate purity given boil point and boiler mol fraction */
double calc_yd(double atm_p, double mol_fract, double boil_pt,
double (*calc_gamma_e)(double x1, double x2),
double (*calc_gamma_w)(double x1, double x2),
double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x));
void write_vle(double *boil_pt, double *xw, double *yd, double *W, double *D, double *dist_avg);
//Antoine Coefficients (C)
double eA = 8.13484;
double eB = 1662.48;
double eC = 238.131;
double wA = 8.05573;
double wB = 1723.64;
double wC = 233.076;
int main()
{
//printf("\nEnter Wash Gallons: ");
double wash_l = 100; //get_input() * GAL2L;
//printf("\nEnter ABV (between 0 and 1): ");
double abv = 0.40; //get_input();
//printf("\nEnter Atmospheric pressure (mmHg): ");
double atm_p = 760; //get_input();
double boil_pt[MAX_NUM] = {0};
double xw[MAX_NUM] = {0};/**Boiler mol fraction (Purity)*/
double W[MAX_NUM] = {0};/**boiler molar volume*/
double D[MAX_NUM] = {0};/**Distillate molar volume*/
double yd[MAX_NUM] = {0};/**Molar fraction of ethanol vapor (Purity)*/
double recip[MAX_NUM] = {0};
double integral[MAX_NUM] = {0};
double dist_avg[MAX_NUM] = {0};
/**Determine initial conditions for VLE calculations*/
xw[0] = calc_mol_fract(abv, abv_to_abw, abw_to_mf_e); /**Purity to start*/
boil_pt[0] = calc_boil_pt(abv, atm_p, xw[0], calc_gamma_e, calc_gamma_w, calc_psat); /**Boil point to start*/
W[0] = calc_total_mol(wash_l, abv); /**Molar volume to start*/
printf("\nInitial boil point = %.3lf", boil_pt[0]);
printf("\nInitial wash purity = %.3lf", xw[0]);
printf("\nInitial wash molar vol = %.3lf", W[0]);
yd[0] = calc_yd(atm_p, xw[0], boil_pt[0], calc_gamma_e, calc_gamma_w, calc_psat);
/**initial values of integral[0] and D[0] distillate volume in mols are zero*/
/**for loop is Rayliegh Eq.*/
for ( int index = 1; index < MAX_NUM; index++)
{
boil_pt[index] = index * 0.2 + boil_pt[0];
/**xw[0] = Initial boiler wash purity*/
xw[index] = calc_xw(atm_p, xw[0], boil_pt[index], calc_gamma_e, calc_gamma_w, calc_psat);/**Boiler liquid purity at instant*/
/**distillate purity = vapor %ethanol in mols- condensed vapor = distillate*/
yd[index] = calc_yd(atm_p, xw[0], boil_pt[index], calc_gamma_e, calc_gamma_w, calc_psat); /**distillate purity at instant*/
recip[index] = 1 / (yd[index] - xw[index]);
integral[index] = integral[index-1]+ (xw[index]- xw[index-1])* ((recip[index] + recip[index-1])/2);
W[index] = W[0] * pow(10, integral[index]); /**Boiler liquid mols at instant*/
D[index] = W[0] - W[index];/**Distillate molar volume*/
dist_avg[index] = (W[0] *xw[0] -xw[index] *yd[index])/D[index];
write_vle(boil_pt, xw, yd, W, D, dist_avg);
}
system("PAUSE");
return 0;
}
/*******************************/
/*********FUNCTION DEFINITIONS*************************/
double get_input(void)
{
double input = 0;
scanf("%lf", &input);
return input;
}
/**************************************************/
//alcohol conversion functions abv to abw, and abw to abv results are between 0 and 1 and only valid for input between 0 and 1
//From "On the Conversion of Alcohol by Edwin Croissant - R1: 2014*02-16
double abw_to_abv(double abw)
{
const double a = -0.000039705486746795932;
const double b = 1.2709666849144778;
const double c = -0.40926819348115739;
const double d = 2.0463351302912738;
const double f = -7.8964816507513707;
const double g = 15.009692673927390;
const double h = -15.765836469736477;
const double i = 8.8142267038252680;
const double j = -2.0695760421183493;
double temp = j;
temp = temp * abw +i;
temp = temp * abw +h;
temp = temp * abw +g;
temp = temp * abw +f;
temp = temp * abw +d;
temp = temp * abw +c;
temp = temp * abw +b;
double abv = temp * abw +a;
return abv;
}
/*******************************/
double abv_to_abw(double abv)
{
const double a = 0.00018684999875047631;
const double b = 0.77602465132552556;
const double c = 0.41803095099103116;
const double d = -2.5221614925275091;
const double f = 9.5827123045656251;
const double g = -19.928886159385002;
const double h = 24.165120890385651;
const double i = -15.830262207383321;
const double j = 4.3390473620304988;
double temp = j;
temp = temp * abv +i;
temp = temp * abv +h;
temp = temp * abv +g;
temp = temp * abv +f;
temp = temp * abv +d;
temp = temp * abv +c;
temp = temp * abv +b;
double abw = temp * abv +a;
return abw;
}
/*******************************/
double abw_to_mf_e(double abw)
{
double mf_e = abw / (abw + MM_E / MM_W * (1-abw));
return mf_e;
}
/*******************************/
double calc_total_mol(double wash_l, double abv)
{
double water_mol = wash_l *(1-abv) *1000 *W_DENS /MM_W;
double eth_mol = wash_l *1000 *abv *E_DENS /MM_E;
return water_mol + eth_mol;
}
/***********************************/
double calc_mol_fract(double abv, double (*abv_to_abw)(double abv), double (*abw_to_mf_e)(double abw))
{
double abw = abv_to_abw(abv);
double eth_mol_fract = abw_to_mf_e(abw);
return eth_mol_fract;
}
/***********************************/
double calc_boil_pt(double abv, double atm_p, double mol_fract,
double (*calc_gamma_e)(double x1, double x2),
double (*calc_gamma_w)(double x1, double x2),
double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x))
{
double bp = 70; //Seed temperature at loop beginning
double x1 = mol_fract; //ethanol component of mixture in mols
double x2 = 1-x1; //water component of mixture in mols
double psat1 = 0;
double psat2 = 0;
double gamma_e = 0;
double gamma_w = 0;
gamma_e = calc_gamma_e(x1,x2);
gamma_w = calc_gamma_w(x1,x2);
do
{
psat1 = calc_psat(eA, eB, eC, bp, gamma_e, x1);
psat2 = calc_psat(wA, wB, wC, bp, gamma_w, x2);
bp = bp + 0.0001;
}
while (psat1 + psat2 < atm_p);
return bp;
}
/************************************************/
/*Calculate Van Laar activity coefficient for ethanol*/
double calc_gamma_e(double x1, double x2)
{
double const a12 = 1.68811;/*a12 and a21 are Van Laar activity model constants*/
double const a21 = 0.95268;
double result = exp((a12 * pow(x2, 2)) / (pow(((a12 * x1 / a21) + x2), 2)));
return (result);
}
/***************************************************/
/*Calculate Van Laar activity coefficient for water*/
double calc_gamma_w(double x1, double x2)
{
double const a12 = 1.68811;
double const a21 = 0.95268;
double result = exp((a21 * pow(x1, 2)) / (pow(((a21 * x2 / a12) + x1), 2)));
return (result);
}
/***************************************************/
/* Antoine Eq. to calculate saturated vapor pressure for water and ethanol */
double calc_psat(double A, double B, double C, double bp, double gamma, double x)
{
double psat = x*gamma*pow(10, (A-((B)/(C + bp))));
return psat;
}
/***************************************************/
/***************************************************/
double calc_xw(double atm_p, double mol_fract, double bp,
double (*calc_gamma_e)(double x1, double x2),
double (*calc_gamma_w)(double x1, double x2),
double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x))
{
double x1 = mol_fract; //ethanol component of mixture in mols
double x2 = 1-x1; //water component of mixture in mols
double psat1 = 0;
double psat2 = 0;
double gamma_e = 0;
double gamma_w = 0;
do
{
gamma_e = calc_gamma_e(x1,x2);
gamma_w = calc_gamma_w(x1,x2);
psat1 = calc_psat(eA, eB, eC, bp, gamma_e, x1);
psat2 = calc_psat(wA, wB, wC, bp, gamma_w, x2);
x1 -= 0.000001;
x2 = 1-x1;
}
while (psat1 + psat2 > atm_p);
return x1;
}
/**********************************************/
double calc_yd(double atm_p, double mol_fract, double bp,
double (*calc_gamma_e)(double x1, double x2),
double (*calc_gamma_w)(double x1, double x2),
double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x))
{
double x1 = mol_fract; //ethanol component of mixture in mols
double x2 = 1-x1; //water component of mixture in mols
double psat1 = 0;
double psat2 = 0;
double gamma_e = 0;
double gamma_w = 0;
do
{
gamma_e = calc_gamma_e(x1,x2);
gamma_w = calc_gamma_w(x1,x2);
psat1 = calc_psat(eA, eB, eC, bp, gamma_e, x1);
psat2 = calc_psat(wA, wB, wC, bp, gamma_w, x2);
x1 -= 0.000001;
x2 = 1-x1;
}
while (psat1 + psat2 > atm_p);
double y_d = psat1 / atm_p;
return y_d;
}
/*****************************************************************/
void write_vle(double *boil_pt, double *xw, double *yd, double *W, double *D, double *dist_avg)
{
//Write to screen
for(int i = 0; i < MAX_NUM; i++)
{
printf("\nindex\tboil_pt\t\tboiler mol\tdist purity\twash mols\tdist mols\tdist avg\n");
printf("%d\t%.2lf\t\t%.3lf\t\t%.3lf\t\t%.3lf\t\t%.3lf\t\t%.3lf\n",
i, boil_pt[i], xw[i], yd[i], W[i], D[i], dist_avg[i]);
}
//Write to file
/*FILE *pWrite;
pWrite = fopen("vle.dat", "w");
if (pWrite == NULL)
printf("\nFile not opened\n");
else
fprintf(pWrite, "\nboil_pt\t\txw\t\tyd\tW\tD\tdist_avg\n");
for(int i = 0; i < MAX_NUM; i++)
{
fprintf(pWrite,"%.2lf\t\t%.3lf\t%.3lf\t%.3lf\t\t%.3lf\t%.3lf\t%.3lf\t%.3lf\n",
boil_pt[i], xw[i], yd[i], W[i], D[i], D[i]);
}
*/
}