/[ascend]/branches/jacob/models/johnpye/fprops/test/init_mix.c
ViewVC logotype

Contents of /branches/jacob/models/johnpye/fprops/test/init_mix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2973 - (show annotations) (download) (as text)
Thu Jun 18 22:23:04 2015 UTC (3 years, 9 months ago) by jacob
File MIME type: text/x-csrc
File size: 10622 byte(s)
testing routine to solve vapor-liquid equilibrium condition when defining mixtures; it's giving segmentation faults right now.

1 /* ASCEND modelling environment
2 Copyright (C) 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, but
10 WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 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 --
17
18 Free Software Foundation, Inc.
19 59 Temple Place - Suite 330
20 Boston, MA 02111-1307, USA.
21 *//*
22 by Jacob Shealy, June 4-9, 2015
23
24 Initial model of a simple mixture, to get the procedure right. This
25 is in preparation for a general algorithm to find mixing conditions
26 in the ideal-mixture case.
27
28 TODO -
29 add error handling in functions:
30 [x] check that sum of mass fractions is really one
31 we will need to add different values to the FpropsError enum to represent errors
32 how do we want inconsistent mass fractions to be handled? Strategies:
33 1. for an n-component system, only specify (n-1) mass fractions;
34 remaining mass fraction is (1 - (sum over other mass fractions))
35 2. specify all mass fractions, but check that they sum to 1 (within
36 some tolerance), and return error otherwise.
37 3. specify all mass fractions, but check that they sum to 1 (within
38 some tolerance), and adjust value of last mass fraction if
39 incorrect sum occurs only because of that fraction. If incorrect
40 sum occurs before the last mass fraction, return an error. This
41 strategy is very dangerous, since it involves changing
42 user-provided data.
43 Maybe add colors from `color.h' later. Too involved for now.
44 */
45
46 #include "../mixtures/init_mixfuncs.h"
47 #include "../helmholtz.h"
48 #include "../pengrob.h"
49 #include "../fluids.h"
50 #include "../fprops.h"
51 #include "../refstate.h"
52 #include "../sat.h"
53
54 #include <stdio.h>
55 #include <math.h>
56
57 /* Macro/Preprocessor definitions; prefixed with `MIX_' */
58 #define NSIMS 5
59
60 extern const EosData eos_rpp_nitrogen;
61 extern const EosData eos_rpp_ammonia;
62 extern const EosData eos_rpp_carbon_dioxide;
63 extern const EosData eos_rpp_methane;
64 extern const EosData eos_rpp_water;
65
66 /* Function prototypes */
67 double rachford_rice(double V, void *user_data);
68 void mixture_flash(MixturePhaseState *M, double P, char **Names, FpropsError *err);
69 void solve_mixture_conditions(unsigned n_sims, double *Ts, double *Ps, MixtureSpec *M, char **Names, FpropsError *err);
70 void densities_Ts_to_mixture(MixtureState *MS, double *P_out, double *T_in, double tol, char **Names, FpropsError *err);
71
72 /*
73 Structure to pass constant parameters into the Rachford-Rice function
74 */
75 typedef struct{
76 unsigned pures;
77 double *zs;
78 double *Ks;
79 } RRData;
80
81 /*
82 Calculate the Rachford-Rice function for a given set of parameters and value
83 of the vapor fraction for a mixture in phase equilibrium
84 */
85 double rachford_rice(double V, void *user_data){
86 RRData *RR = (RRData *)user_data;
87 double F=0.0;
88 unsigned i;
89
90 /* Add each term of the Rachford-Rice function to the output F */
91 for(i=0;i<RR->pures;i++){
92 F += (RR->zs[i] * (RR->Ks[i] - 1)) / (1 + V*(RR->Ks[i] - 1));
93 }
94 return F;
95 }
96
97 /*
98 Determine whether the mixture given is in a one-phase (vapor or liquid) or
99 two-phase state, and if it is in two phases (VLE), what the mass fractions
100 within each phase are, and what fraction of the mass is in each phase.
101 */
102 void mixture_flash(MixturePhaseState *M, double P, char **Names, FpropsError *err){
103 #define NPURE M->X->pures
104 #define ZS M->X->Xs
105 #define XS M->Xs
106 #define PFL M->X->PF
107 #define D PFL[i]->data
108 #define RHOS M->rhos
109 unsigned i;
110 double p_sat, /* temporary saturation pressure */
111 p_b=0.0, /* bubble-point and dew-point pressures */
112 p_d,
113 rp_d=0.0; /* reciprocal dew-point pressure */
114 double Ks[NPURE]; /* K-factor: ratio (vapor mole fraction)/(liquid mole fraction) */
115 double V_frac[]={0.50, 0.52}; /* vapor mass fraction */
116 double xs[NPURE],
117 ys[NPURE],
118 zs[NPURE];
119 double XM_sum=0.0, /* sum of (X_i * MolarMass_i) */
120 YM_sum=0.0,
121 ZM_sum=0.0; /* sum of (Z_i / MolarMass_i) */
122
123 for(i=0;i<NPURE;i++){
124 zs[i] = ZS[i] / D->M; /* first step in calculating overall mole fractions */
125 ZM_sum += zs[i];
126 }
127 for(i=0;i<NPURE;i++){
128 printf("\n\tFinding partial pressure and overall mole fractions for substance:"
129 "\n\t %s", Names[i]);
130 zs[i] /= ZM_sum; /* divide by reciprocal average molar mass to obtain mole fractions */
131 printf("\n\t%.5f\t%.5f\t%.5f\t%.5f", RHOS[0][0], RHOS[0][1], RHOS[1][0], RHOS[1][1]);
132 fprops_sat_T(M->T, &p_sat, (RHOS[1]+i), (RHOS[0]+i), PFL[i], err);
133 p_b += zs[i] * p_sat;
134 rp_d += zs[i] / p_sat;
135 Ks[i] = p_sat / P; /* value of K-factor under Raoult's law */
136 }
137 p_d = 1./rp_d;
138
139 if(P<=p_d){ /* system will be all-gas */
140 M->ph = GAS_PHASE;
141 puts("\n\tMixture is all gas/vapor.");
142 }else if(P>=p_b){ /* system will be all-liquid */
143 M->ph = LIQ_PHASE;
144 puts("\n\tMixture is all liquid.");
145 }else{ /* system may be in vapor-liquid equilibrium */
146 RRData RR = {NPURE, zs, Ks};
147 double tol = 1.e-3;
148
149 secant_solve(&rachford_rice, &RR, V_frac, tol);
150
151 for(i=0;i<NPURE-1;i++){
152 printf("\n\n\tFor substance %s", Names[i]);
153 ys[i] = (zs[i] * Ks[i]) / ((1 - V_frac[0]) + V_frac[0]*Ks[i]);
154 printf("\n\t vapor mole fraction is %.5f", ys[i]);
155 xs[i] = zs[i] / (1 + (V_frac[0] * (Ks[i] - 1)));
156 printf("\n\t liquid mole fraction is %.5f", xs[i]);
157 /* printf("\n\tThe vapor-phase mole fraction is %.5f,"
158 " and the liquid-phase mole fraction is %.5f", ys[i], xs[i]); */
159 }
160 ys[NPURE-1] = 1 - my_sum(NPURE-1, ys); /* ys[0] starts as 0.0, so the fractions should sum correctly */
161 xs[NPURE-1] = 1 - my_sum(NPURE-1, xs);
162 // printf("\n\n\tFor substance %s", Names[NPURE-1]);
163 // printf("\n\t vapor mole fraction is %.5f", ys[NPURE-1]);
164 // printf("\n\t liquid mole fraction is %.5f", xs[NPURE-1]);
165
166 for(i=0;i<NPURE;i++){
167 printf("\n\n\tFor substance number %u", i);
168 printf("\n\tThat is, substance %s", Names[i]);
169 XS[0][i] = ys[i] * D->M; /* first step in calculating mass fractions */
170 printf("\n\t vapor mass fraction is %.5f", XS[0][i]);
171 XS[1][i] = xs[i] * D->M;
172 printf("\n\t liquid mass fraction is %.5f", XS[1][i]);
173
174 XM_sum += xs[i] * D->M;
175 YM_sum += ys[i] * D->M;
176 }
177 for(i=0;i<NPURE;i++){
178 printf("\n\tFinal mass fractions for substance %s", Names[i]);
179 XS[1][i] /= YM_sum;
180 XS[0][i] /= XM_sum;
181 printf("\n\tAssigning densities for substance %s", Names[i]);
182 }
183 // M->rhos = rhos;
184 // M->Xs = Xs;
185 M->phase_splits[0] = V_frac[0];
186 M->phase_splits[1] = 1 - V_frac[0];
187 printf("\n\n\tThe mole split between vapor and liquid is"
188 " %.5f vapor, %.5f liquid\n",
189 V_frac[0], 1 - V_frac[0]);
190 printf("\n\tThe mass fractions in the different phases are:"
191 "\n\t\tvapor\t -- %s %.5f"
192 "\n\t\t\t -- %s %.5f"
193 "\n\t\tliquid\t -- %s %.5f"
194 "\n\t\t\t -- %s %.5f\n",
195 // Names[0], M->Xs[0][0],
196 // Names[1], M->Xs[0][1],
197 // Names[0], M->Xs[1][0],
198 // Names[1], M->Xs[1][1]);
199 Names[0], XS[0][0],
200 Names[1], XS[0][1],
201 Names[0], XS[1][0],
202 Names[1], XS[1][1]);
203 }
204
205 #undef RHOS
206 #undef D
207 #undef PFL
208 #undef XS
209 #undef ZS
210 #undef NPURE
211 }
212
213 /*
214 Establish mixing conditions (T,P), find individual densities, and use those
215 along with the temperature to find first-law properties for individual
216 components and the whole mixture.
217
218 As a check, find pressure corresponding to temperature and component density
219 for each component -- this should equal the overall pressure
220
221 The fluids I use in the mixture are nitrogen, ammonia, carbon dioxide, and
222 methane.
223
224 I may add other substances later, such as methyl chloride, carbon monoxide,
225 nitrous oxide, and hydrogen sulfide.
226 */
227 int main(){
228 int i; /* counter variable */
229 enum FluidAbbrevs {/* N2, */NH3,/* CO2,CH4, */H2O,NFLUIDS}; /* fluids that will be used */
230
231 char *fluids[]={
232 /* "nitrogen", */ "ammonia", /* "carbondioxide", "methane", */ "water"
233 };
234 char *fluid_names[]={
235 /* "Nitrogen", */ "Ammonia", /* "Carbon Dioxide", "Methane", */ "Water"
236 };
237 const EosData *IdealEos[]={
238 &eos_rpp_nitrogen,
239 &eos_rpp_ammonia,
240 &eos_rpp_carbon_dioxide,
241 &eos_rpp_methane,
242 &eos_rpp_water
243 };
244
245 PureFluid *Helms[NFLUIDS];
246 PureFluid *Pengs[NFLUIDS];
247 PureFluid *Ideals[NFLUIDS];
248 ReferenceState ref = {FPROPS_REF_REF0};
249 FpropsError err = FPROPS_NO_ERROR;
250
251 /*
252 Fill the `Helms' PureFluid array with data from the helmholtz equation
253 of state.
254 */
255 for(i=0;i<NFLUIDS;i++){
256 Helms[i] = fprops_fluid(fluids[i],"helmholtz",NULL);
257 Pengs[i] = fprops_fluid(fluids[i],"pengrob",NULL);
258 Ideals[i] = ideal_prepare(IdealEos[i],&ref);
259 }
260
261 /* Mixing conditions (temperature,pressure), and mass fractions */
262 double T[]={250, 300, 300, 350, 350, 400}; /* K */
263 double P[]={1.5e5, 1.5e5, 1.9e5, 1.9e5, 2.1e5, 2.1e5}; /* Pa */
264 double x[NFLUIDS]; /* mass fractions */
265
266 double props[] = {1, 3, 2, 1.5, 2.5}; /* proportions of mass fractions */
267 mixture_x_props(NFLUIDS, x, props); /* mass fractions from proportions */
268
269 double rho1[]={
270 2, 3, 2.5, 1.7, 1.9
271 };
272 double rho2[]={
273 1.1, 2, 1.5, 1.7, 1.9
274 };
275
276 /*
277 choose which model to use by commenting and uncommenting the array names
278 (only one name uncommented at a time)
279 */
280 MixtureSpec MX = {
281 .pures=NFLUIDS
282 , .Xs=x
283 , .PF=Helms
284 /* .PF=Ideals */
285 /* .PF=Pengs */
286 };
287
288 double mp_ps[] = {0.0, 0.0};
289 double mp_xs[2][2] = {
290 {0.0, 0.0}, {0.0, 0.0}
291 };
292 double mp_rhos[2][2] = {
293 {0.0, 0.0}, {0.0, 0.0}
294 };
295 MixturePhaseState MP = {
296 .T=T[2]
297 , .X=&MX
298 , .phase_splits = mp_ps
299 , .Xs = mp_xs
300 , .rhos = mp_rhos
301 };
302
303 mixture_flash(&MP, P[2], fluid_names, &err);
304 printf("\n\tAt %.1f K and %.0f Pa, the mixture is in vapor-liquid equilibrium:"
305 "\n\t %.5f of the mass is in the vapor"
306 "\n\t %s mass fraction in the vapor is %.5f"
307 "\n\t %s mass fraction in the vapor is %.5f"
308 "\n\t %.5f of the mass is in the liquid"
309 /* "\n\t %s mass fraction in the liquid is %.5f" */
310 /* "\n\t %s mass fraction in the liquid is %.5f" */
311 "\n",
312 MP.T, P[2], MP.phase_splits[0], fluid_names[0], MP.Xs[0][0],
313 fluid_names[1], MP.Xs[0][1], MP.phase_splits[1]/* , fluid_names[0],
314 MP.Xs[1][0], fluid_names[1], MP.Xs[1][1] */);
315
316 return 0;
317 } /* end of `main' */
318

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