johnpye/fprops/sat.c

/*  ASCEND modelling environment
    Copyright (C) 2008-2009 Carnegie Mellon University

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*//** @file
    Routines to calculate saturation properties using Helmholtz correlation
    data. We first include some 'generic' saturation equations that make use
    of the acentric factor and critical point properties to predict saturation
    properties (pressure, vapour density, liquid density). These correlations
    seem to be only very rough in some cases, but it is hoped that they will
    be suitable as first-guess values that can then be fed into an iterative
    solver to converge to an accurate result.
*/

#include "sat.h"
#include "helmholtz_impl.h"

# include <assert.h>
#include <math.h>
#include <stdio.h>

#define SQ(X) ((X)*(X))

#if 0
#define _GNU_SOURCE
#include <fenv.h>
#endif

/**
    Estimate of saturation pressure using H W Xiang ''The new simple extended
    corresponding-states principle: vapor pressure and second virial
    coefficient'', Chemical Engineering Science,
    57 (2002) pp 1439-1449.
*/
double fprops_psat_T_xiang(double T, const HelmholtzData *d){

    double Zc = d->p_c / (8314. * d->rho_c * d->T_c);

#ifdef TEST
    fprintf(stderr,"Zc = %f\n",Zc);
#endif

    double theta = SQ(Zc - 0.29);

#ifdef TEST
    fprintf(stderr,"theta = %f\n",theta);
#endif

    double aa[] = {5.790206, 4.888195, 33.91196
                , 6.251894, 15.08591, -315.0248
                , 11.65859, 46.78273, -1672.179
    };

    double a0 = aa[0] + aa[1]*d->omega + aa[2]*theta;
    double a1 = aa[3] + aa[4]*d->omega + aa[5]*theta;
    double a2 = aa[6] + aa[7]*d->omega + aa[8]*theta;

    double Tr = T / d->T_c;
    double tau = 1 - Tr;

    double taupow = pow(tau, 1.89);
    double temp = a0 + a1 * taupow + a2 * taupow * taupow * taupow;

    double logpr = temp * log(Tr);
    double p_r = exp(logpr);

#ifdef TEST
    fprintf(stderr,"a0 = %f\n", a0);
    fprintf(stderr,"a1 = %f\n", a1);
    fprintf(stderr,"a2 = %f\n", a2);
    fprintf(stderr,"temp = %f\n", temp);
    fprintf(stderr,"T_r = %f\n", Tr);
    fprintf(stderr,"p_r = %f\n", p_r);
#endif

    return p_r * d->p_c;
}

/**
    Estimate saturation pressure using acentric factor. This algorithm
    is used for first estimates for later refinement in the program REFPROP.
*/
double fprops_psat_T_acentric(double T, const HelmholtzData *d){
    /* first guess using acentric factor */
    double p = d->p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));
    return p;
}


/**
    Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
    see http://dx.doi.org/10.1002/aic.690250412
*/
double fprops_rhof_T_rackett(double T, const HelmholtzData *D){

    double Zc = D->rho_c * D->R * D->T_c / D->p_c;
    double Tau = 1. - T/D->T_c;
    double vf = (D->R * D->T_c / D->p_c) * pow(Zc, -1 - pow(Tau, 2./7));

    return 1./vf;
}


/**
    Saturated vapour density correlation of Chouaieb, Ghazouani, Bellagi
    see http://dx.doi.org/10.1016/j.tca.2004.05.017
*/
double fprops_rhog_T_chouaieb(double T, const HelmholtzData *D){
    double Zc = D->rho_c * D->R * D->T_c / D->p_c;
    double Tau = 1. - T/D->T_c;
#if 0
# define N1 -0.1497547
# define N2 0.6006565 
# define P1 -19.348354
# define P2 -41.060325
# define P3 1.1878726
    double MMM = 2.6; /* guess, reading from Chouaieb, Fig 8 */
    //MMM = 2.4686277;
    double PPP = Zc / (P1 + P2*Zc*log(Zc) + P3/Zc);
    fprintf(stderr,"PPP = %f\n",PPP);
    //PPP = -0.6240188;
    double NNN = PPP + 1./(N1*D->omega + N2);   
#else
/* exact values from Chouaieb for CO2 */
# define MMM 2.4686277
# define NNN 1.1345838
# define PPP -0.6240188
#endif

    double alpha = exp(pow(Tau,1./3) + sqrt(Tau) + Tau + pow(Tau, MMM));
    return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
}

int fprops_sat_T(double T, double *psat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
    double tau = d->T_c / T;
    double delf = 1.1 * fprops_rhof_T_rackett(T,d) / d->rho_c;
    double delg = 0.9 * fprops_rhog_T_chouaieb(T,d) / d->rho_c;
    //feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);

    int i = 0;
    while(i++ < 50){
        //fprintf(stderr,"%s: iter %d: rhof = %f, rhog = %f\n",__func__,i,delf*d->rho_c, delg*d->rho_c);
        double phirf = helm_resid(tau,delf,d);
        double phirf_d = helm_resid_del(tau,delf,d);
        double phirf_dd = helm_resid_deldel(tau,delf,d);
        double phirg = helm_resid(tau,delg,d);
        double phirg_d = helm_resid_del(tau,delg,d);
        double phirg_dd = helm_resid_deldel(tau,delg,d);

#define J(FG) (del##FG * (1. + del##FG * phir##FG##_d))
#define K(FG) (del##FG * phir##FG##_d + phir##FG + log(del##FG))
#define J_del(FG) (1 + 2 * del##FG * phir##FG##_d + SQ(del##FG) * phir##FG##_dd)
#define K_del(FG) (2 * phir##FG##_d + del##FG * phir##FG##_dd + 1./del##FG)
        double Jf = J(f);
        double Jg = J(g);
        double Kf = K(f);
        double Kg = K(g);
        double Jf_del = J_del(f);
        double Jg_del = J_del(g);
        double Kf_del = K_del(f);
        double Kg_del = K_del(g);

        double DELTA = Jg_del * Kf_del - Jf_del * Kg_del;

#define gamma 1
        delf += gamma/DELTA * ((Kg - Kf) * Jg_del - (Jg - Jf) * Kg_del);
        delg += gamma/DELTA * ((Kg - Kf) * Jf_del - (Jg - Jf) * Kf_del);

        if(fabs(Kg - Kf) + fabs(Jg - Jf) < 1e-8){
            //fprintf(stderr,"%s: CONVERGED\n",__func__);
            *rhof_out = delf * d->rho_c;
            *rhog_out = delg * d->rho_c;
            *psat_out = helmholtz_p_raw(T, *rhog_out, d);
            return 0;
        }
    }
    *rhof_out = delf * d->rho_c;
    *rhog_out = delg * d->rho_c;
    *psat_out = helmholtz_p_raw(T, *rhog_out, d);
    fprintf(stderr,"%s: NOT CONVERGED for '%s' with T = %e (rhof=%f, rhog=%f)\n",__func__,d->name,T,*rhof_out,*rhog_out);
    return 1;

}

int fprops_sat_p(double p, double *Tsat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
    fprintf(stderr,"%s: NOT YET IMPLEMENTED\n",__func__);
    return 1;
}


1	jpye	2122	/* ASCEND modelling environment
2			Copyright (C) 2008-2009 Carnegie Mellon University
3	jpye	2114
4	jpye	2122	This program is free software; you can redistribute it and/or modify
5			it under the terms of the GNU General Public License as published by
6			the Free Software Foundation; either version 2, or (at your option)
7			any later version.
8
9			This program is distributed in the hope that it will be useful,
10			but WITHOUT ANY WARRANTY; without even the implied warranty of
11			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12			GNU General Public License for more details.
13
14			You should have received a copy of the GNU General Public License
15			along with this program; if not, write to the Free Software
16			Foundation, Inc., 59 Temple Place - Suite 330,
17			Boston, MA 02111-1307, USA.
18			//* @file
19			Routines to calculate saturation properties using Helmholtz correlation
20			data. We first include some 'generic' saturation equations that make use
21			of the acentric factor and critical point properties to predict saturation
22			properties (pressure, vapour density, liquid density). These correlations
23			seem to be only very rough in some cases, but it is hoped that they will
24			be suitable as first-guess values that can then be fed into an iterative
25			solver to converge to an accurate result.
26	jpye	2114	*/
27
28			#include "sat.h"
29	jpye	2262	#include "helmholtz_impl.h"
30	jpye	2114
31	jpye	2124	# include <assert.h>
32	jpye	2114	#include <math.h>
33	jpye	2262	#include <stdio.h>
34	jpye	2114
35			#define SQ(X) ((X)*(X))
36
37	jpye	2262	#if 0
38			#define _GNU_SOURCE
39			#include <fenv.h>
40			#endif
41
42	jpye	2122	/**
43			Estimate of saturation pressure using H W Xiang ''The new simple extended
44			corresponding-states principle: vapor pressure and second virial
45			coefficient'', Chemical Engineering Science,
46			57 (2002) pp 1439-1449.
47			*/
48	jpye	2114	double fprops_psat_T_xiang(double T, const HelmholtzData *d){
49
50			double Zc = d->p_c / (8314. * d->rho_c * d->T_c);
51
52	jpye	2117	#ifdef TEST
53	jpye	2114	fprintf(stderr,"Zc = %f\n",Zc);
54	jpye	2117	#endif
55	jpye	2114
56			double theta = SQ(Zc - 0.29);
57
58	jpye	2117	#ifdef TEST
59	jpye	2114	fprintf(stderr,"theta = %f\n",theta);
60	jpye	2117	#endif
61	jpye	2114
62			double aa[] = {5.790206, 4.888195, 33.91196
63			, 6.251894, 15.08591, -315.0248
64			, 11.65859, 46.78273, -1672.179
65			};
66
67			double a0 = aa[0] + aa[1]d->omega + aa[2]theta;
68			double a1 = aa[3] + aa[4]d->omega + aa[5]theta;
69			double a2 = aa[6] + aa[7]d->omega + aa[8]theta;
70
71			double Tr = T / d->T_c;
72			double tau = 1 - Tr;
73
74			double taupow = pow(tau, 1.89);
75			double temp = a0 + a1 * taupow + a2 * taupow * taupow * taupow;
76
77			double logpr = temp * log(Tr);
78			double p_r = exp(logpr);
79
80			#ifdef TEST
81			fprintf(stderr,"a0 = %f\n", a0);
82			fprintf(stderr,"a1 = %f\n", a1);
83			fprintf(stderr,"a2 = %f\n", a2);
84			fprintf(stderr,"temp = %f\n", temp);
85			fprintf(stderr,"T_r = %f\n", Tr);
86			fprintf(stderr,"p_r = %f\n", p_r);
87			#endif
88
89			return p_r * d->p_c;
90			}
91
92	jpye	2115	/**
93	jpye	2230	Estimate saturation pressure using acentric factor. This algorithm
94			is used for first estimates for later refinement in the program REFPROP.
95			*/
96			double fprops_psat_T_acentric(double T, const HelmholtzData *d){
97			/* first guess using acentric factor */
98			double p = d->p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));
99			return p;
100			}
101
102
103			/**
104	jpye	2120	Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
105			see http://dx.doi.org/10.1002/aic.690250412
106			*/
107	jpye	2121	double fprops_rhof_T_rackett(double T, const HelmholtzData *D){
108	jpye	2120
109			double Zc = D->rho_c * D->R * D->T_c / D->p_c;
110			double Tau = 1. - T/D->T_c;
111			double vf = (D->R * D->T_c / D->p_c) * pow(Zc, -1 - pow(Tau, 2./7));
112
113			return 1./vf;
114			}
115
116
117			/**
118			Saturated vapour density correlation of Chouaieb, Ghazouani, Bellagi
119			see http://dx.doi.org/10.1016/j.tca.2004.05.017
120			*/
121	jpye	2121	double fprops_rhog_T_chouaieb(double T, const HelmholtzData *D){
122	jpye	2120	double Zc = D->rho_c * D->R * D->T_c / D->p_c;
123			double Tau = 1. - T/D->T_c;
124			#if 0
125			# define N1 -0.1497547
126			# define N2 0.6006565
127			# define P1 -19.348354
128			# define P2 -41.060325
129			# define P3 1.1878726
130			double MMM = 2.6; /* guess, reading from Chouaieb, Fig 8 */
131	jpye	2122	//MMM = 2.4686277;
132			double PPP = Zc / (P1 + P2Zclog(Zc) + P3/Zc);
133			fprintf(stderr,"PPP = %f\n",PPP);
134			//PPP = -0.6240188;
135	jpye	2120	double NNN = PPP + 1./(N1*D->omega + N2);
136			#else
137	jpye	2122	/* exact values from Chouaieb for CO2 */
138	jpye	2120	# define MMM 2.4686277
139			# define NNN 1.1345838
140			# define PPP -0.6240188
141			#endif
142
143			double alpha = exp(pow(Tau,1./3) + sqrt(Tau) + Tau + pow(Tau, MMM));
144	jpye	2121	return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
145	jpye	2120	}
146
147	jpye	2262	int fprops_sat_T(double T, double psat_out, double rhof_out, double * rhog_out, const HelmholtzData *d){
148			double tau = d->T_c / T;
149			double delf = 1.1 * fprops_rhof_T_rackett(T,d) / d->rho_c;
150			double delg = 0.9 * fprops_rhog_T_chouaieb(T,d) / d->rho_c;
151			//feenableexcept(FE_DIVBYZERO \| FE_INVALID \| FE_OVERFLOW);
152	jpye	2120
153	jpye	2262	int i = 0;
154			while(i++ < 50){
155			//fprintf(stderr,"%s: iter %d: rhof = %f, rhog = %f\n",__func__,i,delfd->rho_c, delgd->rho_c);
156			double phirf = helm_resid(tau,delf,d);
157			double phirf_d = helm_resid_del(tau,delf,d);
158			double phirf_dd = helm_resid_deldel(tau,delf,d);
159			double phirg = helm_resid(tau,delg,d);
160			double phirg_d = helm_resid_del(tau,delg,d);
161			double phirg_dd = helm_resid_deldel(tau,delg,d);
162	jpye	2119
163	jpye	2262	#define J(FG) (del##FG * (1. + del##FG * phir##FG##_d))
164			#define K(FG) (del##FG * phir##FG##_d + phir##FG + log(del##FG))
165			#define J_del(FG) (1 + 2 * del##FG * phir##FG##_d + SQ(del##FG) * phir##FG##_dd)
166			#define K_del(FG) (2 * phir##FG##_d + del##FG * phir##FG##_dd + 1./del##FG)
167			double Jf = J(f);
168			double Jg = J(g);
169			double Kf = K(f);
170			double Kg = K(g);
171			double Jf_del = J_del(f);
172			double Jg_del = J_del(g);
173			double Kf_del = K_del(f);
174			double Kg_del = K_del(g);
175	jpye	2119
176	jpye	2262	double DELTA = Jg_del * Kf_del - Jf_del * Kg_del;
177	jpye	2122
178	jpye	2262	#define gamma 1
179			delf += gamma/DELTA * ((Kg - Kf) * Jg_del - (Jg - Jf) * Kg_del);
180			delg += gamma/DELTA * ((Kg - Kf) * Jf_del - (Jg - Jf) * Kf_del);
181	jpye	2119
182	jpye	2262	if(fabs(Kg - Kf) + fabs(Jg - Jf) < 1e-8){
183			//fprintf(stderr,"%s: CONVERGED\n",__func__);
184			rhof_out = delf d->rho_c;
185			rhog_out = delg d->rho_c;
186			psat_out = helmholtz_p_raw(T, rhog_out, d);
187			return 0;
188			}
189	jpye	2124	}
190	jpye	2262	rhof_out = delf d->rho_c;
191			rhog_out = delg d->rho_c;
192			psat_out = helmholtz_p_raw(T, rhog_out, d);
193			fprintf(stderr,"%s: NOT CONVERGED for '%s' with T = %e (rhof=%f, rhog=%f)\n",__func__,d->name,T,rhof_out,rhog_out);
194			return 1;
195	jpye	2124
196	jpye	2122	}
197
198	jpye	2262	int fprops_sat_p(double p, double Tsat_out, double rhof_out, double * rhog_out, const HelmholtzData *d){
199			fprintf(stderr,"%s: NOT YET IMPLEMENTED\n",__func__);
200			return 1;
201	jpye	2123	}
202
203
204	jpye	2124
205
206	jpye	2123
207	jpye	2124