/[ascend]/trunk/models/johnpye/fprops/sat.c
ViewVC logotype

Annotation of /trunk/models/johnpye/fprops/sat.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2262 - (hide annotations) (download) (as text)
Thu Aug 5 08:34:43 2010 UTC (9 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 6306 byte(s)
Deleted sat2.c and sat2.h, that code is not being used anywhere in ASCEND/FPROPS now. REFPROP algorithm replaced by Akasaka approach.
Moving sat3.c code into sat.h.
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 + P2*Zc*log(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,delf*d->rho_c, delg*d->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

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22