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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2264 - (show annotations) (download) (as text)
Thu Aug 5 09:25:55 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 7235 byte(s)
A basic fprops_sat_p function working now.
Need to add p_t, T_t data to all materials.
1 /* ASCEND modelling environment
2 Copyright (C) 2008-2009 Carnegie Mellon University
3
4 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 */
27
28 #include "sat.h"
29 #include "helmholtz_impl.h"
30
31 # include <assert.h>
32 #include <math.h>
33 #include <stdio.h>
34
35 #define SQ(X) ((X)*(X))
36
37 #if 0
38 #define _GNU_SOURCE
39 #include <fenv.h>
40 #endif
41
42 /**
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 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 #ifdef TEST
53 fprintf(stderr,"Zc = %f\n",Zc);
54 #endif
55
56 double theta = SQ(Zc - 0.29);
57
58 #ifdef TEST
59 fprintf(stderr,"theta = %f\n",theta);
60 #endif
61
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 /**
93 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 Saturated liquid density correlation of Rackett, Spencer & Danner (1972)
105 see http://dx.doi.org/10.1002/aic.690250412
106 */
107 double fprops_rhof_T_rackett(double T, const HelmholtzData *D){
108
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 double fprops_rhog_T_chouaieb(double T, const HelmholtzData *D){
122 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 //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 double NNN = PPP + 1./(N1*D->omega + N2);
136 #else
137 /* exact values from Chouaieb for CO2 */
138 # 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 return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
145 }
146
147 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
153 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
163 #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
176 double DELTA = Jg_del * Kf_del - Jf_del * Kg_del;
177
178 #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
182 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 }
190 *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
196 }
197
198 int fprops_sat_p(double p, double *Tsat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
199 /*
200 Estimate of saturation temperature using definition of acentric factor and
201 the assumed p(T) relationship:
202 log10(p)=A + B/T
203 See Reid, Prausnitz and Poling, 4th Ed., section 2.3.
204 */
205 double T1 = d->T_c / (1. - 3./7. / (1.+d->omega) * log10(p / d->p_c));
206 double p1, rhof, rhog;
207 int i = 0;
208 while(i++ < 20){
209 int res = fprops_sat_T(T1, &p1, &rhof, &rhog, d);
210 if(res)return 1;
211 //fprintf(stderr,"%s: T1 = %f ������> p = %f bar\trhof = %f\trhog = %f\n",__func__, T1, p1/1e5, rhof, rhog);
212 if(fabs(p1 - p) < 1e-5){
213 *Tsat_out = T1;
214 *rhof_out = rhof;
215 *rhog_out = rhog;
216 return 0;
217 }
218 double hf = helmholtz_h_raw(T1, rhof, d);
219 double hg = helmholtz_h_raw(T1, rhog, d);
220 double dpdT_sat = (hg - hf) / T1 / (1./rhog - 1./rhof);
221 //fprintf(stderr,"\t\tdpdT_sat = %f bar/K\n",dpdT_sat/1e5);
222 double delta_T = -(p1 - p)/dpdT_sat;
223 if(T1 + delta_T < d->T_t)T1 = 0.5 * (d->T_t + T1);
224 else T1 += delta_T;
225 }
226 *Tsat_out = T1;
227 *rhof_out = rhof;
228 *rhog_out = rhog;
229 return 1;
230 }
231
232
233
234
235
236

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