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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3022 by jacob, Tue Jul 7 11:15:55 2015 UTC revision 3023 by jacob, Mon Jul 20 23:00:30 2015 UTC
# Line 19  Line 19 
19      59 Temple Place - Suite 330      59 Temple Place - Suite 330
20      Boston, MA 02111-1307, USA.      Boston, MA 02111-1307, USA.
21  *//*  *//*
22      by Jacob Shealy, June 4-9, 2015      by Jacob Shealy, June 4-July 20, 2015
23            
24      Initial model of a simple mixture, to get the procedure right.  This      Initial model of a simple mixture, to get the procedure right.  This
25      is in preparation for a general algorithm to find mixing conditions      is in preparation for a general algorithm to find mixing conditions
# Line 66  extern const EosData eos_rpp_water; Line 66  extern const EosData eos_rpp_water;
66    
67  /* Function prototypes */  /* Function prototypes */
68  SecantSubjectFunction rachford_rice;  SecantSubjectFunction rachford_rice;
 SecantSubjectFunction pressure_rho_error;  
69  SecantSubjectFunction dew_p_error;  SecantSubjectFunction dew_p_error;
70  SecantSubjectFunction bubble_p_error;  SecantSubjectFunction bubble_p_error;
 void transit_non_physical(double T, double *rhos, const PureFluid *PF, FpropsError *err, int from_vap);  
 double rho_from_pressure(double T, double p, double rho_start, double tol, const PureFluid *PF, FpropsError *err);  
71  double dew_pressure(MixtureSpec *MS, double T, FpropsError *err);  double dew_pressure(MixtureSpec *MS, double T, FpropsError *err);
72  double bubble_pressure(MixtureSpec *MS, double T, FpropsError *err);  double bubble_pressure(MixtureSpec *MS, double T, FpropsError *err);
73    
# Line 79  void mixture_flash(MixturePhaseState *M, Line 76  void mixture_flash(MixturePhaseState *M,
76  void solve_mixture_conditions(unsigned n_sims, double *Ts, double *Ps, MixtureSpec *M, char **Names, FpropsError *err);  void solve_mixture_conditions(unsigned n_sims, double *Ts, double *Ps, MixtureSpec *M, char **Names, FpropsError *err);
77  void densities_Ts_to_mixture(MixtureState *MS, double *P_out, double *T_in, double tol, char **Names, FpropsError *err);  void densities_Ts_to_mixture(MixtureState *MS, double *P_out, double *T_in, double tol, char **Names, FpropsError *err);
78    
79  void pengrob_phi_all(double *phi_out, MixtureSpec *M, double T, double *rhos, FpropsError *err);  double pengrob_phi_pure(PureFluid *PF, double T, double P, PhaseName type, FpropsError *err);
80  double pengrob_phi_pure(PureFluid *PF, double T, double rho, FpropsError *err);  double poynting_factor(PureFluid *PF, double T, double P, FpropsError *err);
81  double poynting_factor(PureFluid *PF, double T, double rho, FpropsError *err);  
82    // void test_one(double *Ts, double *Ps);
83  void test_one(double *Ts, double *Ps);  // void test_two(double T, double P);
84  void test_two(double T, double P);  // void test_three(double T, double *rhos);
85  void test_three(double T, double *rhos);  // void test_four(double T, double P);
 void test_four(double T, double P);  
86  void test_five(double T, double P);  void test_five(double T, double P);
87  void test_six(void);  void test_six(void);
88    void test_seven(void);
89    
90  /*  /*
91      Establish mixing conditions (T,P), find individual densities, and use those      Establish mixing conditions (T,P), find individual densities, and use those
# Line 121  int main(){ Line 118  int main(){
118      /* test_two(T[1], P[1]); */      /* test_two(T[1], P[1]); */
119      /* test_three(T[1], rho1); */      /* test_three(T[1], rho1); */
120      /* test_four(T[5], P[6]); */      /* test_four(T[5], P[6]); */
121      test_five(270, 1.5e5);      /* test_five(270, 1.5e5); */
122      /* test_six(); */      /* test_six(); */
123        test_seven();
124    
125      return 0;      return 0;
126  } /* end of `main' */  } /* end of `main' */
# Line 244  void mixture_flash(MixturePhaseState *M, Line 242  void mixture_flash(MixturePhaseState *M,
242  #undef NPURE  #undef NPURE
243  }  }
244    
245  void pengrob_phi_all(double *phi_out, MixtureSpec *M, double T, double *rhos, FpropsError *err){  int cubic_solution(double coef[4], double *roots){
246  #define NPURE M->pures      double p, q, r,
247  #define D M->PF[i]->data             a, b,
248  #define RR D->R             cond;
249  #define PR D->corr.pengrob  
250      unsigned i;      p = coef[1]/coef[0];
251      double P,     /* pressure */      q = coef[2]/coef[0];
252             aT,    /* parameter `a' at the temperature */      r = coef[3]/coef[0];
253             alpha, /* intermediate variable in calculating aT */      a = (3*q - pow(p, 2)) / 3;
254             Tr,    /* reduced temperature */      b = (2*pow(p, 3) - (9*p*q) + (27*r)) / 27;
255             A,     /* two variables */      cond = (pow(b, 2) / 4) + (pow(a, 3) / 27);
            B,  
            Z;     /* compressibility */  
     double cpx_ln; /* rather involved part of the fugacity coefficient */  
       
     for(i=0;i<NPURE;i++){ /* find fugacity coefficient for each component */  
         P = fprops_p((FluidState){T, rhos[i], M->PF[i]}, err);  
         /* printf("\n  Substance number %u\n\tPressure is now %.0f Pa", i, P); */  
   
         Tr = T/D->T_c;  
         alpha = (1 + PR->kappa * (1 - sqrt(Tr)));  
         aT = PR->aTc * alpha * alpha;  
   
         A = aT * P / (RR * RR * T * T);  
         B = PR->b * P / (RR * T);  
         Z = P / (rhos[i] * RR * T);  
         /* printf("\n\tCompressibility is %.6g;\n\t`A' is %.6g\n\t`B' is %.6g", Z, A, B); */  
   
         cpx_ln = log( (Z + (1 + M_SQRT2)*B)/(Z + (1 - M_SQRT2)*B) );  
         phi_out[i] = exp( Z - 1 - log(Z - B) - (A / 2 / M_SQRT2) * cpx_ln );  
         /* printf("\n\tLog of fugacity coefficient is %.6g", Z - 1 - log(Z - B) - (A / (2*M_SQRT2))*cpx_ln); */  
         /* printf("\n\tFugacity coefficient is %.6g", phi_out[i]); */  
     }  
256    
257  #undef PR  #if 0
258  #undef RR      printf("\n  <cubic_solution>: variables in this cubic-formula solution are:"
259  #undef D              "\n\tp = %.6g"
260  #undef NPURE              "\n\tq = %.6g"
261                "\n\tr = %.6g"
262                "\n\ta = %.6g"
263                "\n\tb = %.6g"
264                "\n\tcond = %.6g"
265                , p, q, r, a, b, cond);
266    #endif
267    
268        if(cond>0){
269            roots[0] = cbrt(-b/2 + sqrt(cond)) + cbrt(-b/2 - sqrt(cond)) - p/3;
270            return 1;
271        }else if(cond<0){
272            double R, phi;
273    
274            R = pow(b, 2) / (4 * (-pow(a, 3)/27));
275            phi = acos(sqrt(R));
276    
277    #if 0
278            printf("\n  <cubic_solution>: there will be three roots"
279                    "\n\tR = %.6g"
280                    "\n\tphi = %.6g"
281                    , R, phi);
282    #endif
283    
284            int i;
285            for(i=0;i<3;i++){
286                if(b>0){
287                    roots[i] = -2 * sqrt(-a/3) * cos(phi/3 + (MIX_PI*2*i)/3) - p/3;
288                }else{
289                    roots[i] = 2 * sqrt(-a/3) * cos(phi/3 + (MIX_PI*2*i)/3) - p/3;
290                }
291    #if 0
292                printf("\n\troot %u is"
293                        "\n\t  %.1g * sqrt(%.9g) * cos(%.9g + %.9g) ="
294                        "\n\t  %.1g * %.9g * %.9g = %.9g"
295                        , i, (b>0) ? -2.0 : 2.0, -a/3, phi/3, (MIX_PI * 2 * i)/3
296                        , (b>0) ? -2.0 : 2.0, sqrt(-a/3), cos(phi/3 + (MIX_PI*2*i)/3)
297                        , roots[i]);
298    #endif
299            }
300            return 3;
301        }else{ /* cond equals zero */
302            double A;
303    
304            A = sqrt(a/3) * ((b>0) ? -1 : 1);
305            roots[0] = 2*A - p/3;
306            roots[1] = -A - p/3;
307            /* roots[2] = -A - p/3; */
308    
309            return 2;
310        }
311        return 0; /* error -- no roots returned */
312  }  }
313    
314  double pengrob_phi_pure(PureFluid *PF, double T, double rho, FpropsError *err){  double pengrob_phi_pure(PureFluid *PF, double T, double P, PhaseName type, FpropsError *err){
315  #define RR PF->data->R  #define RR PF->data->R
316  #define PR PF->data->corr.pengrob  #define PR PF->data->corr.pengrob
317      double P,     /* pressure */      unsigned i
318             Tr,    /* reduced temperature */          , z_num
319             aT,    /* parameter `a' at the temperature */          , z_cons1 /* tests whether the expressions for Z have been derived properly... */
320             alpha, /* intermediate variable in calculating aT */          , z_cons2 /* tests whether roots of cubic-form equations have been solved properly... */
321             A,     /* two variables */          ;
322             B,      double Tr       /* reduced temperature */
323             Z,     /* compressibility */          , aT        /* parameter `a' at the temperature */
324             cpx_ln; /* rather involved part of the fugacity coefficient */          , alpha     /* intermediate variable in calculating aT */
325            , A         /* two variables */
326            , B
327            , z_coef[4] /* coefficients of cubic expression for compressibility */
328            , Z_rt[3]   /* roots of cubic expression for compressibility */
329            , Z
330            , cpx_ln    /* rather involved part of the fugacity coefficient */
331            , tol = 1.e-6
332            ;
333            
334      P = fprops_p((FluidState){T, rho, PF}, err);      printf("\n  <pengrob_phi_pure>: Pressure is now %.0f Pa", P);
     printf("\n\tPressure is now %.0f Pa", P);  
335      Tr = T / PF->data->T_c;      Tr = T / PF->data->T_c;
336      alpha = (1 + PR->kappa * (1 - sqrt(Tr)));      alpha = (1 + PR->kappa * (1 - sqrt(Tr)));
337      aT = PR->aTc * alpha * alpha;      aT = PR->aTc * alpha * alpha;
338    
339      A = aT * P / (RR * RR * T * T);      A = aT * P / (RR * RR * T * T);
340      B = PR->b * P / (RR * T);      B = PR->b * P / (RR * T);
341      Z = P / (rho * RR * T);      /* if(type<=VAPOR){ */
342      printf("\n\tCompressibility is %.6g;\n\t`A' is %.6g\n\t`B' is %.6g", Z, A, B);          z_coef[0] = 1;
343            z_coef[1] = B - 1;
344            z_coef[2] = A - 2*B - 3*pow(B, 2);
345            z_coef[3] = pow(B, 3) + pow(B, 2) - A*B;
346        /* }else{ */
347    #if 0
348            z_coef[0] = - 1./A;
349            z_coef[1] = (1 - B)/A;
350            z_coef[2] = (2*B + 3*pow(B, 2))/A;
351            z_coef[3] = -(pow(B, 2) + pow(B, 3))/A + B;
352    /* #else */
353            z_coef[0] = 1;
354            z_coef[1] =
355            z_coef[2] =
356            z_coef[3] =
357    #endif
358        /* } */
359        z_num = cubic_solution(z_coef, Z_rt);
360        if(z_num==1){
361            Z = Z_rt[0];
362        }else{
363            if(type<=VAPOR){
364                Z = max_element(z_num, Z_rt);
365            }else{
366                if(min_positive_elem(&Z, z_num, Z_rt)){
367                    printf("\n  <pengrob_phi_pure>:" MIX_COMPR_ERROR " when calculating"
368                            "\n\tfugacity coefficient for %s phase; compressibility is %.6g"
369                            , (type<=VAPOR) ? "vapor/gas" : "liquid", Z);
370                    /* err = FPROPS_RANGE_ERROR; */
371                }
372            }
373        }
374    #if 0
375        printf("\n\t<pengrob_phi_pure>: There are %u roots of the equation"
376                "\n\t  %.6g x^3 + %.6g x^2 + %.6g x + %.6g = 0 :"
377                , z_num, z_coef[0], z_coef[1], z_coef[2], z_coef[3]);
378        for(i=0;i<z_num;i++){
379            printf("\n\t  %.6g", Z_rt[i]);
380        }
381    #endif
382    
383        if(type<=VAPOR){
384            z_cons1 = (fabs(Z - (1 + B - A * (Z - B)/(pow(Z, 2) + (2*Z*B) - pow(B, 2)))) < tol);
385        }else{
386            z_cons1 = (fabs(Z - (B + (pow(Z, 2) + (2*Z*B) - pow(B, 2))*(1 + B - Z)/A)) < tol);
387        }
388        z_cons2 = (fabs(z_coef[0]*pow(Z, 3) + z_coef[1]*pow(Z, 2) + z_coef[2]*Z + z_coef[3]) < tol);
389    
390        printf("\n\t<pengrob_phi_pure>: Compressibility for %s of %s is %.6g;"
391                "\n\t  `A' is %.6g"
392                "\n\t  `B' is %.6g"
393                /* "\n\t  Derivation of Z %s consistent with original equations." */
394                /* "\n\t  Value of Z %s consistent with cubic-form equation." */
395                , (type<=VAPOR) ? "vapor/gas" : "liquid", PF->name, Z, A, B
396                /* , (z_cons1) ? "is" : "is not", (z_cons2) ? "is" : "is not" */);
397    
398      cpx_ln = log( (Z + (1 + M_SQRT2)*B)/(Z + (1 - M_SQRT2)*B) );      cpx_ln = log( (Z + (1 + M_SQRT2)*B)/(Z + (1 - M_SQRT2)*B) );
399      return exp( Z - 1 - log(Z - B) - (A / 2 / M_SQRT2) * cpx_ln );      return exp( Z - 1 - log(Z - B) - (A / 2 / M_SQRT2) * cpx_ln );
# Line 318  double pengrob_phi_pure(PureFluid *PF, d Line 406  double pengrob_phi_pure(PureFluid *PF, d
406      deviation of a liquid in phase equilibrium that results from the pressure      deviation of a liquid in phase equilibrium that results from the pressure
407      effects due to not being at the saturation pressure.      effects due to not being at the saturation pressure.
408   */   */
409  double poynting_factor(PureFluid *PF, double T, double rho, FpropsError *err){  double poynting_factor(PureFluid *PF, double T, double P, FpropsError *err){
410      double p,      double /* p, */
411             p_sat,             p_sat,
412             rho_liq, rho_vap,             rho_liq, rho_vap,
413             ln_pf; /* natural logarithm of Poynting Factor */             ln_pf; /* natural logarithm of Poynting Factor */
414    
415      p = fprops_p((FluidState){T, rho, PF}, err);      /* p = fprops_p((FluidState){T, rho, PF}, err); */
416    
417      fprops_sat_T(T, &p_sat, &rho_liq, &rho_vap, PF, err);      fprops_sat_T(T, &p_sat, &rho_liq, &rho_vap, PF, err);
418    
419      ln_pf = (p - p_sat) / (rho_liq * PF->data->R * T);      ln_pf = (P - p_sat) / (rho_liq * PF->data->R * T);
420      return exp(ln_pf);      return exp(ln_pf);
421  }  }
422    
 /* Finding density from pressure */  
 /*  
     Structure to hold auxiliary data for function to find the error in pressure  
     at a given density  
  */  
 typedef struct PressureRhoData_Struct {  
     double T;  
     double P;  
     PureFluid *pfl;  
     FpropsError *err;  
 } PRData;  
   
 /*  
     Return error in pressure given a certain density  
  */  
 /* double pressure_rho_error(double rho, void *user_data){  
     PRData *prd = (PRData *)user_data;  
     FluidState fst = {prd->T, rho, prd->pfl};  
       
     return fabs(prd->P - fprops_p(fst, prd->err));  
 } */  
   
 /*  
     Return reasonable value of density at a given temperature and pressure, and  
     from a given starting density.  
  */  
 double rho_from_pressure(double T, double P, double rho_start, double tol, const PureFluid *PF, FpropsError *err){  
 #define MAX_ITER 20  
 #define ALPHA 2.35  
     unsigned i, j,  
              ix_low,  /* index of point (density, pressure) with lower error in P */  
              ix_high; /* index of point with higher error in P */  
     int transit=0;    /* records whether the search has passed through the region of non-physical (dP/d(rho)) behavior */  
     double ps[2],  
            /* delta_rho, */  
            rhos[] = {rho_start, 1.02*rho_start};  
     double rho_new[3],  
            p_new[3];  
     double A = 1.05,  
            B = 2.35,  
            C = 0.43;  
     PRData prd = {T, P, PF, err};  
   
     ps[1] = pressure_rho_error(rhos[1], &prd);  
   
     for(i=0;i<MAX_ITER;i++){  
         ps[0] = pressure_rho_error(rhos[0], &prd);  
         ix_low = (fabs(ps[1]) < fabs(ps[0])) ? 1 : 0; /* which position has smallest error? */  
         ix_high = 1 - ix_low;  
   
         if(fabs(ps[0])<tol){  
             printf("\n\t<rho_from_pressure>: Search for correct density SUCCEEDED after %u iterations "  
                     "at rho=%.6g kg/m^3"  
                     , i, rhos[0]);  
             return rhos[0];  
         }  
         if(fabs(rhos[0]-rhos[1])<tol){  
             /*  
                 Algorithm has converged on an extremum that does not satisfy the  
                 pressure given.  Transition through non-physical region if it  
                 has not already done so; otherwise, return an error.  
              */  
             if(!transit){  
                 transit = 1; /* mark that algorithm has transitioned */  
                 transit_non_physical(T, rhos, PF, err, (ps[0]<P) ? 1 : 0);  
                 printf("\n\t<rho_from_pressure>: Transitioning across non-physical region");  
             }else{  
                 printf("\n\t<rho_from_pressure>: Search for correct density FAILED after %u iterations,"  
                         "\n\t  densities equal at rho[0]=rho[1]=%.6g"  
                         , i, rhos[0]);  
                 return rhos[0];  
             }  
         }  
         if(rhos[0]==INFINITY || rhos[1]==INFINITY  
                 || rhos[0]!=rhos[0] || rhos[1]!=rhos[1]){  
             printf("\n\t<rho_from_pressure>: Search for correct density FAILED after %u iterations,"  
                     "\n\t  densities infinite or not-a-number at rho[0]=%.6g, rho[1]=%.6g"  
                     , i, rhos[0], rhos[1]);  
             return rhos[0];  
         }  
   
         /*  
               
          */  
         /* if((ps[0]-ps[1])/(rhos[0]-rhos[1]) > 0 || ){  
             delta_rho = -ps[0] * (rhos[0] - rhos[1]) / (ps[0] - ps[1]);  
             / *  
                 If magnitude of change in density is greater than a certain  
                 amount, reduce it to that magnitude.  This prevents the density  
                 changes from growing too large, and is necessary for a  
                 non-monotonic function like this one.  
              * /  
             if(fabs(delta_rho) > ALPHA * fabs(rhos[0]-rhos[1])){  
                 delta_rho *= ALPHA / fabs(delta_rho);  
             }  
             rhos[1] = rhos[ix_low]; / * reassign second to lowest-error position * /  
             ps[1] = ps[ix_low];  
             rho[0] += delta_rho;  
         } */  
         rho_new[0] = (1 + A)*rhos[ix_low] - A*rhos[ix_high];  
         rho_new[1] = (1 - B)*rhos[ix_low] + B*rho_new[0];  
         rho_new[2] = (1 - C)*rhos[ix_low] + C*rhos[ix_high];  
   
         for(j=0;j<3;j++){  
             p_new[j] = fabs(pressure_rho_error(rho_new[j], &prd));  
         }  
   
         rhos[1] = rhos[ix_low];  
         ps[1] = ps[ix_low];  
         rhos[0] = rho_new[index_of_min(3, p_new)];  
   
         if(*err!=FPROPS_NO_ERROR){  
             *err = FPROPS_NO_ERROR;  
         }  
     }  
   
     printf("\n\t<rho_from_pressure>: Reached maximum number of iterations (%u) without converging on density;"  
             "\n\trhos[0]=%.8g, rhos[1]=%.8g, P[0]=%.8g, P[1]=%.8g, tolerance=%g"  
             , MAX_ITER, rhos[0], rhos[1], ps[0], ps[1], tol);  
     return rhos[0];  
 #undef MAX_ITER  
 }  
   
 /*  
     Transition across the non-physical portion of the pressure-density curve,  
     where dP/d(rho) < 0, that is, where density decreases with increasing  
     pressure.  
   
     For this function, it is essential that the starting point be at or very  
     close to the peak of the P/rho curve, so that the first step takes the  
     density solidly into the non-physical region.  
  */  
 void transit_non_physical(double T, double *rhos, const PureFluid *PF, FpropsError *err, int from_vap){  
     unsigned i=0;  
     double delta_rho;  
   
     if(from_vap){  
         rhos[1] = 1.05 * rhos[0];  
         printf("\n  <transit_non_physical>: rhos[0]=%.6g, rhos[1]=%.6g, P[0]=%.1f Pa, P[1]=%.1f"  
                 , rhos[0], rhos[1]  
                 , fprops_p((FluidState){T, rhos[0], PF}, err)  
                 , fprops_p((FluidState){T, rhos[1], PF}, err));  
         while(fprops_p((FluidState){T, rhos[0], PF}, err) > fprops_p((FluidState){T, rhos[1], PF}, err)){  
             delta_rho = 1.5 * (rhos[1] - rhos[0]);  
             rhos[0] = rhos[1];    /* rho[0], the trailing edge of the line segment, is shifted first */  
             rhos[1] += delta_rho; /* then rho[1], the leading edge, is shifted */  
             i++;  
         }  
         rhos[0] = 1.05 * rhos[1];  
     }else{  
         rhos[0] = 0.95 * rhos[1];  
         printf("\n  <transit_non_physical>: rhos[0]=%.6g, rhos[1]=%.6g, P[0]=%.1f Pa, P[1]=%.1f"  
                 , rhos[0], rhos[1]  
                 , fprops_p((FluidState){T, rhos[0], PF}, err)  
                 , fprops_p((FluidState){T, rhos[1], PF}, err));  
         while(fprops_p((FluidState){T, rhos[0], PF}, err) > fprops_p((FluidState){T, rhos[1], PF}, err)){  
             rhos[1] = rhos[0]; /* rho[1] is now the trailing edge of the line segment, and is shifted first */  
             rhos[0] *= 0.7;    /* then rho[0], the leading edge, is shifted */  
             i++;  
         }  
         rhos[1] = 0.95 * rhos[0];  
     }  
     printf("\n\t<transit_non_physical>: Shifted to rho[0]=%.8g, after %u iterations", rhos[0], i);  
 }  
   
423  /* Finding phase-equilibrium conditions */  /* Finding phase-equilibrium conditions */
424  /*  /*
425      Apportion components into either critical/supercritical or saturation/VLE      Apportion components into either critical/supercritical or saturation/VLE
# Line 531  double dew_p_error(double P_D, void *use Line 454  double dew_p_error(double P_D, void *use
454  #define SATRHO dbd->sat_rhos  #define SATRHO dbd->sat_rhos
455  #define TOL dbd->tol  #define TOL dbd->tol
456  #define ERR dbd->err  #define ERR dbd->err
457      printf("\n  <dew_p_error>: entered the function");      /* printf("\n  <dew_p_error>: Entered the function"); */
458      unsigned i;      unsigned i;
459      double p_d_term,      double /* p_d_term, */
460             rp_d = 0.0;             rp_d = 0.0;
461    
462      DBData *dbd = (DBData *)user_data;      DBData *dbd = (DBData *)user_data;
463      printf("\n\t%s", "<dew_p_error>: unpacked the user_data struct");      /* printf("\n\t<dew_p_error>: Unpacked the user_data struct for mixture "
464                "at %.2f K, %.0f Pa"
465                , TT, P_D); */
466    
467      for(i=0;i<NPURE;i++){      for(i=0;i<NPURE;i++){
468          RHOS[i] = rho_from_pressure(TT, P_D, RHOS[i], TOL, PF[i], ERR);          rp_d += YS[i] * pengrob_phi_pure(PF[i], TT, P_D, VAPOR, ERR) / PSAT[i];
         /* RHOV[i] = rho_from_pressure(TT, P_D, RHOV[i], TOL, PF[i], ERR); */  
         /* RHOL[i] = rho_from_pressure(TT, P_D, RHOL[i], TOL, PF[i], ERR); */  
   
         /* p_d_term = YS[i] * pengrob_phi_pure(PF[i], TT, RHOS[i], ERR) / PSAT[i]; */  
   
         rp_d += YS[i] * pengrob_phi_pure(PF[i], TT, RHOS[i], ERR) / PSAT[i];  
469      }      }
470      return P_D - 1./rp_d;      return P_D - 1./rp_d;
471  #if 0  #if 0
# Line 567  double dew_p_error(double P_D, void *use Line 486  double dew_p_error(double P_D, void *use
486    
487  double bubble_p_error(double P_B, void *user_data){  double bubble_p_error(double P_B, void *user_data){
488  #define XS dbd->MS->Xs  #define XS dbd->MS->Xs
489        /* printf("\n  <bubble_p_error>: Entered the function"); */
490      unsigned i;      unsigned i;
491      double p_b_term,      double /* p_b_term, */
492             p_b = 0.0;             p_b = 0.0;
493    
494      DBData *dbd = (DBData *)user_data;      DBData *dbd = (DBData *)user_data;
495        /* printf("\n\t<bubble_p_error>: Unpacked the user_data struct for mixture "
496                "at %.2f K, %.0f Pa"
497                , TT, P_B); */
498    
499      for(i=0;i<NPURE;i++){      for(i=0;i<NPURE;i++){
500          RHOS[i] = rho_from_pressure(TT, P_B, RHOS[i], TOL, PF[i], ERR);          p_b += XS[i] * PSAT[i] / pengrob_phi_pure(PF[i], TT, P_B, VAPOR, ERR);
         /* RHOV[i] = rho_from_pressure(TT, P_B, RHOV[i], TOL, PF[i], ERR); */  
         /* RHOL[i] = rho_from_pressure(TT, P_B, RHOL[i], TOL, PF[i], ERR); */  
   
         /* p_b_term = XS[i] * pengrob_phi_pure(PF[i], TT, SATRHO[i], ERR)  
             * PSAT[i] * poynting_factor(PF[i], TT, RHOL[i], ERR);  
         p_b_term /= pengrob_phi_pure(PF[i], TT, RHOV[i], ERR); */  
   
         p_b += XS[i] * PSAT[i] / pengrob_phi_pure(PF[i], TT, RHOS[i], ERR);  
501      }      }
502      return P_B - p_b;      return P_B - p_b;
503  #undef ERR  #undef ERR
# Line 607  double dew_pressure(MixtureSpec *MS, dou Line 522  double dew_pressure(MixtureSpec *MS, dou
522      unsigned i;      unsigned i;
523      double p_d[2],      double p_d[2],
524             rp_d = 0.0;             rp_d = 0.0;
525      double p_sat[NPURE],      double p_sat[NPURE]
526             rho_l[NPURE],          , rho_l[NPURE]
527             rho_v[NPURE];          , rho_v[NPURE]
528            , xs[NPURE];
529      double tol = 1.e-5;      double tol = 1.e-5;
530        mole_fractions(NPURE, xs, XS, PF);
531    
532      /* find ideal-gas dew pressure as a starting point */      /* Find ideal-gas dew pressure as a starting point */
533      for(i=0;i<NPURE;i++){      for(i=0;i<NPURE;i++){
534          printf("\n  <dew_pressure>: Calculating saturation pressure for substance %s (at index %u)"          /* printf("\n  <dew_pressure>: Calculating saturation pressure for "
535                  "\n\tfrom temperature T=%.2f K...", PF[i]->name, i, T);                  "substance %s (at index %u)"
536                    "\n\tfrom temperature T=%.2f K...", PF[i]->name, i, T); */
537          fprops_sat_T(T, (p_sat+i), (rho_l+i), (rho_v+i), PF[i], err);          fprops_sat_T(T, (p_sat+i), (rho_l+i), (rho_v+i), PF[i], err);
538          rp_d += XS[i] / p_sat[i];          rp_d += xs[i] / p_sat[i];
539          printf("\n\tsaturation pressure was P_sat=%.0f Pa", p_sat[i]);          /* printf("\n\tsaturation pressure was P_sat=%.0f Pa", p_sat[i]); */
540      }      }
541      p_d[0] = 1./rp_d;      p_d[0] = 1./rp_d;
542      printf("\n  <dew_pressure>: Ideal-gas dew pressure is P_D = %.0f Pa", p_d[0]);      /* printf("\n  <dew_pressure>: Ideal-gas dew pressure is P_D = %.0f Pa", p_d[0]); */
543      p_d[1] = 1.01/rp_d;      p_d[1] = 1.01/rp_d;
544    
545      DBData dbd = {      DBData dbd = {
# Line 633  double dew_pressure(MixtureSpec *MS, dou Line 551  double dew_pressure(MixtureSpec *MS, dou
551          , .tol = tol          , .tol = tol
552          , .err = err          , .err = err
553      };      };
554      printf("\n  <dew_pressure>: Created DBData struct");      /* printf("\n  <dew_pressure>: Created DBData struct"); */
555      puts("");      /* puts(""); */
556    
557      secant_solve(&dew_p_error, &dbd, p_d, tol);      secant_solve(&dew_p_error, &dbd, p_d, tol);
558    
# Line 648  double dew_pressure(MixtureSpec *MS, dou Line 566  double dew_pressure(MixtureSpec *MS, dou
566      Find the bubble-point pressure for a mixture at a given temperature      Find the bubble-point pressure for a mixture at a given temperature
567   */   */
568  double bubble_pressure(MixtureSpec *MS, double T, FpropsError *err){  double bubble_pressure(MixtureSpec *MS, double T, FpropsError *err){
569      return 0.0;  #define NPURE MS->pures
570  }  #define PF MS->PF
571    #define XS MS->Xs
 /* Functions to test/demonstrate other functionality */  
 /*  
     Calculate and print mixture properties, given several temperatures and  
     densities  
  */  
 void test_one(double *Ts, double *Ps){  
     int i; /* counter variable */  
     enum FluidAbbrevs {N2,NH3,CO2,CH4,/* H2O, */NFLUIDS}; /* fluids that will be used */  
   
     char *fluids[]={  
         "nitrogen", "ammonia", "carbondioxide", "methane", "water"  
     };  
     char *fluid_names[]={  
         "Nitrogen", "Ammonia", "Carbon Dioxide", "Methane", "Water"  
     };  
     /* const EosData *IdealEos[]={  
         &eos_rpp_nitrogen,  
         &eos_rpp_ammonia,  
         &eos_rpp_carbon_dioxide,  
         &eos_rpp_methane,  
         &eos_rpp_water  
     }; */  
   
     PureFluid *Helms[NFLUIDS];  
     /* PureFluid *Pengs[NFLUIDS]; */  
     /* PureFluid *Ideals[NFLUIDS]; */  
     /* ReferenceState ref = {FPROPS_REF_REF0}; */  
     FpropsError err = FPROPS_NO_ERROR;  
   
     /*    
         Fill the `Helms' PureFluid array with data from the helmholtz equation  
         of state.  
      */  
     for(i=0;i<NFLUIDS;i++){  
         Helms[i] = fprops_fluid(fluids[i],"helmholtz",NULL);  
         /* Pengs[i] = fprops_fluid(fluids[i],"pengrob",NULL); */  
         /* Ideals[i] = ideal_prepare(IdealEos[i],&ref); */  
     }  
   
     double x[NFLUIDS];                    /* mass fractions */  
   
     double props[] = {1, 3, 2, 1.5, 2.5}; /* proportions of mass fractions */  
     mixture_x_props(NFLUIDS, x, props);   /* mass fractions from proportions */  
   
     /*  
         choose which model to use by commenting and uncommenting the array names  
         (only one name uncommented at a time)  
      */  
     MixtureSpec MX = {  
         .pures=NFLUIDS  
         , .Xs=x  
         , .PF=Helms  
         /* .PF=Ideals */  
         /* .PF=Pengs */  
     };  
   
     char **Names = fluid_names;  
     MixtureSpec *M = &MX;  
     unsigned n_sims = 5;  
   
     double rhos[n_sims][M->pures]; /* individual densities */  
     double rho_mix[n_sims],  
            u_mix[n_sims],  
            h_mix[n_sims],  
            cp_mix[n_sims],  
            cv_mix[n_sims],  
            s_mix[n_sims],  
            g_mix[n_sims],  
            a_mix[n_sims];  
   
     int usr_cont=1;  
     double tol = 1e-9;  
     char temp_str[100];  
     char *headers[n_sims];  
   
 #if 0  
 #define MIXTURE_CALC(PROP) mixture_##PROP(n_pure, xs, rhos[i], Ts[i], PFs, err)  
 #define MIXTURE_CALC(PROP) mixture_##PROP(M, Ts[i] err)  
 #else  
 #define MIXTURE_CALC(PROP) mixture_##PROP(&MS, &err)  
 #endif  
     for(i=0;i<n_sims;i++){  
         /*    
             For each set of conditions simulated, find initial densities to use as a  
             starting point, and then find more exact densities to use  
          */  
         MixtureState MS = {  
             .T=Ts[i],  
             .rhos=rhos[i],  
             .X=M,  
         };  
         initial_rhos(&MS, Ps[i], Names, &err);  
         pressure_rhos(&MS, Ps[i], tol, /* Names, */ &err);  
   
         /*  
             Check that user wants to continue (user can interrupt simulation if  
             bad results were encountered)  
          */  
         printf("\n  %s\n\t%s\n\t%s",  
                 "Continue to calculate and print solution properties?",  
                 "0 - No", "1 - Yes");  
         do{  
             printf("\n    %s ", "Choice?");  
             fgets(temp_str, 100, stdin);  
         }while((1 != sscanf(temp_str, "%i", &usr_cont) || (0>usr_cont || 1<usr_cont))  
                 && printf("\n  %s", "You must enter either 0 or 1"));  
         if(1!=usr_cont){  
             break;  
         }  
   
         /* Calculate solution properties */  
         rho_mix[i] = mixture_rho(&MS);  
         u_mix[i] = MIXTURE_CALC(u);  
         h_mix[i] = MIXTURE_CALC(h);  
         cp_mix[i] = MIXTURE_CALC(cp);  
         cv_mix[i] = MIXTURE_CALC(cv);  
         s_mix[i] = MIXTURE_CALC(s);  
         g_mix[i] = MIXTURE_CALC(g);  
         a_mix[i] = MIXTURE_CALC(a);  
   
         /* Fill in header */  
         headers[i] = (char *)malloc(40);  
         snprintf(headers[i], 40, "[T=%.1f K, P=%.2f bar]", Ts[i], Ps[i]/1.e5);  
     }  
 #undef MIXTURE_CALC  
   
     print_cases_properties(n_sims, headers, rho_mix, Ps, u_mix, h_mix, cp_mix, cv_mix, s_mix, g_mix, a_mix);  
 }  
   
 /*  
     Perform vapor-liquid equilibrium on a mixture of water and ammonia, to test  
     function for modeling flash processes  
  */  
 void test_two(double T, double P){  
     int i; /* counter variable */  
     enum FluidAbbrevs {/* N2, */NH3,/* CO2,CH4, */H2O,NFLUIDS}; /* fluids that will be used */  
   
     char *fluids[]={  
         /* "nitrogen", */ "ammonia", /* "carbondioxide", "methane", */ "water"  
     };  
     char *fluid_names[]={  
         /* "Nitrogen", */ "Ammonia", /* "Carbon Dioxide", "Methane", */ "Water"  
     };  
     /* const EosData *IdealEos[]={  
         &eos_rpp_nitrogen,  
         &eos_rpp_ammonia,  
         &eos_rpp_carbon_dioxide,  
         &eos_rpp_methane,  
         &eos_rpp_water  
     }; */  
   
     PureFluid *Helms[NFLUIDS];  
     /* PureFluid *Pengs[NFLUIDS]; */  
     /* ReferenceState ref = {FPROPS_REF_REF0}; */  
     FpropsError err = FPROPS_NO_ERROR;  
   
     /*    
         Fill the `Helms' PureFluid array with data from the helmholtz equation  
         of state.  
      */  
     for(i=0;i<NFLUIDS;i++){  
         Helms[i] = fprops_fluid(fluids[i],"helmholtz",NULL);  
         /* Pengs[i] = fprops_fluid(fluids[i],"pengrob",NULL); */  
     }  
   
     double x[NFLUIDS];                    /* mass fractions */  
   
     double props[] = {1, 3, 2, 1.5, 2.5}; /* proportions of mass fractions */  
     mixture_x_props(NFLUIDS, x, props);   /* mass fractions from proportions */  
   
     /*  
         choose which model to use by commenting and uncommenting the array names  
         (only one name uncommented at a time)  
      */  
     MixtureSpec MX = {  
         .pures=NFLUIDS  
         , .Xs=x  
         , .PF=Helms  
         /* .PF=Ideals */  
         /* .PF=Pengs */  
     };  
   
     double mp_ps[] = {0.0, 0.0, 0.0};  
     unsigned mp_pn[] = {0, 0, 0};  
     double *mp_xs[] = {NULL, NULL, NULL};  
     double *mp_rhos[] = {NULL, NULL, NULL};  
     MixturePhaseState MP = {  
         .T=T  
         , .X=&MX  
         , .phases = 0  
         , .ph_name = mp_pn  
         , .ph_frac = mp_ps  
         , .Xs = mp_xs  
         , .rhos = mp_rhos  
     };  
     MP.Xs[0] = (double *)malloc(2*sizeof(double));  
     MP.Xs[1] = (double *)malloc(2*sizeof(double));  
     MP.rhos[0] = (double *)malloc(2*sizeof(double));  
     MP.rhos[1] = (double *)malloc(2*sizeof(double));  
   
     mixture_flash(&MP, P, fluid_names, &err);  
     printf("\n\tAt %.1f K and %.0f Pa, the mixture is in vapor-liquid equilibrium:"  
             "\n\t  %.5f of the mass is in the vapor"  
             "\n\t    %s mass fraction in the vapor is %.5f"  
             "\n\t    %s mass fraction in the vapor is %.5f"  
             "\n\t  %.5f of the mass is in the liquid"  
             "\n\t    %s mass fraction in the liquid is %.5f"  
             "\n\t    %s mass fraction in the liquid is %.5f"  
             "\n",  
             MP.T, P, MP.ph_frac[0], fluid_names[0], MP.Xs[0][0],  
             fluid_names[1], MP.Xs[0][1], MP.ph_frac[1], fluid_names[0],  
             MP.Xs[1][0], fluid_names[1], MP.Xs[1][1]);  
     printf("\n\n\tCross-check on mass fractions:");  
     for(i=0;i<NFLUIDS;i++){  
         printf("\n\t  Total %s mass fraction is (%.4f x %.4f) + (%.4f x %.4f) = %.5f",  
                 fluid_names[i], MP.ph_frac[0], MP.Xs[0][i], MP.ph_frac[1], MP.Xs[1][i],  
                 ((MP.ph_frac[0] * MP.Xs[0][i]) + (MP.ph_frac[1] * MP.Xs[1][i])));  
     } puts("");  
     for(i=0;i<NFLUIDS;i++){  
         printf("\n\t  Total %s mass fraction should be %.5f", fluid_names[i], MP.X->Xs[i]);  
     } puts("");  
 }  
   
 void test_three(double T, double *rhos){  
572      unsigned i;      unsigned i;
573      enum FluidAbbrevs {N2,NH3,CO2,CH4,/* H2O, */NFLUIDS}; /* fluids that will be used */      double p_b[] = {0.0, 0.0}
574            , p_sat[NPURE]
575            , rho_l[NPURE]
576            , rho_v[NPURE]
577            , xs[NPURE]
578            , tol = 1.e-5;
579        mole_fractions(NPURE, xs, XS, PF);
580    
581      char *fluids[]={      /* Find ideal-gas bubble pressure as a starting point */
582          "nitrogen", "ammonia", "carbondioxide", "methane", "water"      for(i=0;i<NPURE;i++){
583      };          /* printf("\n  <bubble_pressure>: Calculating saturation pressure for substance "
584      char *fluid_names[]={                  "%s (at index %u)\n\tfrom temperature T=%.2f K...", PF[i]->name, i, T); */
585          "Nitrogen", "Ammonia", "Carbon Dioxide", "Methane", "Water"          fprops_sat_T(T, (p_sat+i), (rho_l+i), (rho_v+i), PF[i], err);
586      };          p_b[0] += xs[i] * p_sat[i];
587      /* const EosData *IdealEos[]={          /* printf("\n\tsaturation pressure was P_sat=%.0f Pa;"
588          &eos_rpp_nitrogen,                  "\n\t  since mole fraction is %.6f, bubble pressure is now %.0f Pa"
589          &eos_rpp_ammonia,                  , p_sat[i], xs[i], p_b[0]); */
         &eos_rpp_carbon_dioxide,  
         &eos_rpp_methane,  
         &eos_rpp_water  
     }; */  
   
     PureFluid *Helms[NFLUIDS];  
     /* PureFluid *Pengs[NFLUIDS]; */  
     /* PureFluid *Ideals[NFLUIDS]; */  
     /* ReferenceState ref = {FPROPS_REF_REF0}; */  
     FpropsError err = FPROPS_NO_ERROR;  
   
     /*    
         Fill the `Helms' PureFluid array with data from the helmholtz equation  
         of state.  
      */  
     for(i=0;i<NFLUIDS;i++){  
         Helms[i] = fprops_fluid(fluids[i],"helmholtz",NULL);  
         /* Pengs[i] = fprops_fluid(fluids[i],"pengrob",NULL); */  
         /* Ideals[i] = ideal_prepare(IdealEos[i],&ref); */  
590      }      }
591        /* printf("\n  <bubble_pressure>: Ideal-gas bubble pressure is P_B = %.0f Pa", p_b[0]); */
592        p_b[1] = 1.01*p_b[0];
593    
594      double x[NFLUIDS];                    /* mass fractions */      DBData dbd = {
595            .MS = MS
596      double props[] = {1, 3, 2, 1.5, 2.5}; /* proportions of mass fractions */          , .T = T
597      mixture_x_props(NFLUIDS, x, props);   /* mass fractions from proportions */          , .p_sat = p_sat
598            , .rhos = rho_v
599      /*          , .sat_rhos = rho_l
600          choose which model to use by commenting and uncommenting the array names          , .tol = tol
601          (only one name uncommented at a time)          , .err = err
      */  
     MixtureSpec MX = {  
         .pures=NFLUIDS  
         , .Xs=x  
         , .PF=Helms  
         /* , .PF=Ideals */  
         /* , .PF=Pengs */  
     };  
   
     MixtureState MS = {  
         .T=T  
         , .rhos=rhos  
         , .X = &MX  
     };  
   
     double tol=1.e-6;  
   
     densities_to_mixture(&MS, tol, fluid_names, &err);  
 }  
   
 void test_four(double T, double P){  
     unsigned i;  
     enum FluidAbbrevs {N2,NH3,CO2,CH4,/* H2O, */NFLUIDS}; /* fluids that will be used */  
   
     char *fluids[]={  
         "nitrogen", "ammonia", "carbondioxide", "methane", "water"  
     };  
     char *fluid_names[]={  
         "Nitrogen", "Ammonia", "Carbon Dioxide", "Methane", "Water"  
602      };      };
603        /* printf("\n  <bubble_pressure>: Created DBData struct\n"); */
604    
605      PureFluid *Helms[NFLUIDS];      secant_solve(&bubble_p_error, &dbd, p_b, tol);
     PureFluid *Pengs[NFLUIDS];  
     FpropsError err = FPROPS_NO_ERROR;  
   
     for(i=0;i<NFLUIDS;i++){  
         Helms[i] = fprops_fluid(fluids[i],"helmholtz",NULL);  
         Pengs[i] = fprops_fluid(fluids[i],"pengrob",NULL);  
     }  
   
     double rhos[NFLUIDS] = {0.0};  
     double x[NFLUIDS];                    /* mass fractions */  
   
     double props[] = {1, 3, 2, 1.5, 2.5}; /* proportions of mass fractions */  
     mixture_x_props(NFLUIDS, x, props);   /* mass fractions from proportions */  
   
     MixtureSpec MX = {  
         .pures=NFLUIDS  
         , .Xs=x  
         , .PF=Pengs  
     };  
     MixtureState MS = {  
         .T=T  
         , .rhos=rhos  
         , .X = &MX  
     };  
606    
607      double phi_1[NFLUIDS],      return p_b[0];
608             phi_2[NFLUIDS];  #undef XS
609      double tol = 1.e-6;  #undef PF
610    #undef NPURE
     initial_rhos(&MS, P, fluid_names, &err);  
     pressure_rhos(&MS, P, tol, /* fluid_names, */ &err);  
   
     pengrob_phi_all(phi_1, &MX, T, MS.rhos, &err);  
     /* pengrob_phi_mix(phi_2, &MS, P, &err); */  
   
     puts("\n  Pure-component fugacity coefficients:");  
     for(i=0;i<NFLUIDS;i++){  
         printf("\t%s %.6g\n", fluid_names[i], phi_1[i]);  
     }  
     puts("\n  Ideal-solution fugacity:");  
     for(i=0;i<NFLUIDS;i++){  
         printf("\t%s %.6g Pa\n", fluid_names[i], phi_1[i] * MS.X->Xs[i] * P);  
     }  
     puts("\n  Mixture fugacity coefficients:");  
     for(i=0;i<NFLUIDS;i++){  
         printf("\t%s %.6g\n", fluid_names[i], phi_2[i]);  
     }  
     puts("\n  Mixture fugacity:");  
     for(i=0;i<NFLUIDS;i++){  
         printf("\t%s %.6g\n", fluid_names[i], phi_2[i] * P);  
     }  
   
 /*  double tau, delta;  
     double B_mix, B_centi, B_milli,  
            d_helm;  
     puts("\n  Attempt to find virial coefficients from Helmholtz EOS:");  
 #define CH Helms[i]->data->corr.helm  
     for(i=0;i<NFLUIDS;i++){  
         tau = CH->T_star / T;  
         delta = rhos[i] / CH->rho_star;  
         printf("\tFor substance %s\n\t  Tau is %.6g\n\t  Delta is %.6g\n",  
                 fluid_names[i], tau, delta);  
   
         B_mix = helm_resid_del(tau, delta, CH) / Helms[i]->data->rho_c;  
         B_centi = helm_resid_del(tau, delta / 100., CH) / Helms[i]->data->rho_c;  
         B_milli = helm_resid_del(tau, delta / 1000., CH) / Helms[i]->data->rho_c;  
         d_helm = helm_resid_del(tau, 0, CH);  
   
         printf("\n\t  V.C. from derivative at mixture conditions: %.6g"  
                 "\n\t  V.C. from derivative at tau = 0.01 of mixture conditions: %.6g"  
                 "\n\t  V.C. from derivative at tau = 0.001 of mixture conditions: %.6g"  
                 "\n\t  derivative at tau = 0 : %.6g" "\n\n" ,B_mix ,B_centi ,B_milli  
                 ,d_helm);  
     } */  
 #undef CH  
611  }  }
612    
613    /* Functions to test/demonstrate other functionality */
614  void test_five(double T, double P){  void test_five(double T, double P){
615      unsigned i;      unsigned i;
616  #define NPURE 4  #define NPURE 4
617  #define D MS->PF[i]->data  #define D MS->PF[i]->data
618    #define TRIALS 5
619      char *fluids[] = {      char *fluids[] = {
620          "isohexane", "krypton", "carbonmonoxide", "ammonia", "water"          "isohexane", "krypton", "carbonmonoxide", "ammonia", "water"
621      };      };
# Line 1065  void test_five(double T, double P){ Line 652  void test_five(double T, double P){
652             rho_v[NPURE] = {0.0}, /* vapor-phase densities */             rho_v[NPURE] = {0.0}, /* vapor-phase densities */
653             rho_sc[NPURE] = {0.0}, /* supercritical densities */             rho_sc[NPURE] = {0.0}, /* supercritical densities */
654             rho_dummy = 0.0;             rho_dummy = 0.0;
655      double p_rho,      double p_b[TRIALS] = {0.0},    /* bubble pressure */
656             p_b = 0.0,    /* bubble pressure */             p_d[TRIALS] = {0.0};          /* dew pressure */
657             rp_d = 0.0,   /* reciprocal dew pressure */      /* double phi_l[NPURE],
            p_d;          /* dew pressure */  
     double phi_l[NPURE],  
658             phi_v[NPURE],             phi_v[NPURE],
659             pf_l[NPURE];             pf_l[NPURE]; */
660      double x_vle[NPURE], /* overall mole fractions in vapor-liquid equilibrium */      double x_vle[NPURE]    /* overall mole fractions in vapor-liquid equilibrium */
661             xs[4][NPURE]; /* per-phase mole fractions */             , xs[4][NPURE]; /* per-phase mole fractions */
662      double x_sum=0.0;    /* sum of mole fractions */      double x_sum=0.0;    /* sum of mole fractions */
663        double Ts[] = {
664            270
665            , 310
666            , 350
667            , 390
668            , 430
669            , 470
670        };
671    
672      MixtureSpec *MS_vle = ASC_NEW(MixtureSpec);      MixtureSpec *MS_vle = ASC_NEW(MixtureSpec);
673  #define VPURE MS_vle->pures  #define VPURE MS_vle->pures
674  #define VXS MS_vle->Xs  #define VXS MS_vle->Xs
675  #define VPF MS_vle->PF  #define VPF MS_vle->PF
676      VPURE = 0;      VPURE = 0;
677      // VXS = ASC_NEW_ARRAY(double,NPURE);      VXS = ASC_NEW_ARRAY(double,NPURE);
678      VXS = (double *)malloc(sizeof(double)*NPURE);      /* VXS = (double *)malloc(sizeof(double)*NPURE); */
679      // VPF = ASC_NEW_ARRAY(PureFluid *,NPURE);      VPF = ASC_NEW_ARRAY(PureFluid *,NPURE);
680      VPF = (PureFluid **)malloc(sizeof(PureFluid *)*NPURE);      /* VPF = (PureFluid **)malloc(sizeof(PureFluid *)*NPURE); */
681    
682      MixtureSpec *MS_sc = ASC_NEW(MixtureSpec);      MixtureSpec *MS_sc = ASC_NEW(MixtureSpec);
683  #define SPURE MS_sc->pures  #define SPURE MS_sc->pures
# Line 1108  void test_five(double T, double P){ Line 701  void test_five(double T, double P){
701              VPF[VPURE] = MS->PF[i]; /* add new component and mass fraction, in next place */              VPF[VPURE] = MS->PF[i]; /* add new component and mass fraction, in next place */
702              VXS[VPURE] = MS->Xs[i];              VXS[VPURE] = MS->Xs[i];
703              VPURE ++;               /* increment number of pures in VLE */              VPURE ++;               /* increment number of pures in VLE */
             /*  
             fprops_sat_T(T, (p_sat+i), (rho_l+i), (rho_v+i), MS->PF[i], &err);  
             fprops_sat_p(P, (T_sat+i), &rho_dummy, (rho_v+i), MS->PF[i], &err);  
   
             x_vle[i] = Xs[i] / D->M;  
             x_sum += x_vle[i];  
             */  
704          }else{          }else{
705              printf("\n  Adding substance %s to supercritical phase", fluid_names[i]);              printf("\n  Adding substance %s to supercritical phase", fluid_names[i]);
706              SPF[SPURE] = MS->PF[i]; /* add new component and mass fraction, in next place */              SPF[SPURE] = MS->PF[i]; /* add new component and mass fraction, in next place */
707              SXS[VPURE] = MS->Xs[i];              SXS[SPURE] = MS->Xs[i];
708              SPURE ++;               /* increment number of pures in supercritical condition */              SPURE ++;               /* increment number of pures in supercritical condition */
709              /* rho_sc[i] = D->rho_c; */              /* rho_sc[i] = D->rho_c; */
710          }          }
         /*  
         if(rho_v[i]!=0 && rho_l[i]!=0){  
             p_rho = fprops_p((FluidState){T,rho_v[i],MS->PF[i]}, &err);  
             if(p_rho<P){  
                 rho_v[i] = rho_from_pressure(T, P, rho_v[i]*P/p_rho, tol, MS->PF[i], &err);  
             }else{  
                 rho_v[i] = rho_from_pressure(T, P, rho_v[i], tol, MS->PF[i], &err);  
             }  
         }else{  
             p_rho = fprops_p((FluidState){T,rho_sc[i],MS->PF[i]}, &err);  
             if(p_rho>P){  
                 rho_sc[i] = rho_from_pressure(T, P, rho_sc[i]*P/p_rho, tol, MS->PF[i], &err);  
             }else{  
                 rho_sc[i] = rho_from_pressure(T, P, rho_sc[i], tol, MS->PF[i], &err);  
             }  
         }  
         */  
711          printf("\n  VPURE = %u, and SPURE = %u; i = %u", VPURE, SPURE, i);          printf("\n  VPURE = %u, and SPURE = %u; i = %u", VPURE, SPURE, i);
712      }      }
713      for(i=0;i<VPURE;i++){      for(i=0;i<VPURE;i++){
# Line 1146  void test_five(double T, double P){ Line 715  void test_five(double T, double P){
715          fprops_sat_p(P, (T_sat+i), &rho_dummy, (rho_v+i), VPF[i], &err);          fprops_sat_p(P, (T_sat+i), &rho_dummy, (rho_v+i), VPF[i], &err);
716    
717          x_vle[i] /= x_sum;          x_vle[i] /= x_sum;
   
         // phi_l[i] = pengrob_phi_pure(MS->PF[i], T, rho_l[i], &err);  
         // phi_v[i] = pengrob_phi_pure(MS->PF[i], T, rho_v[i], &err);  
718      }      }
     /*  
     for(i=0;i<NPURE;i++){  
         if(rho_v[i]!=0 && rho_l[i]!=0){  
         }  
     }  
     */  
719    
     /* p_d = 1./rp_d; */  
720      /*      /*
721      for(i=0;i<NPURE;i++){      for(i=0;i<NPURE;i++){
722          if(rho_v[i]!=0 && rho_l[i]!=0){          if(rho_v[i]!=0 && rho_l[i]!=0){
# Line 1193  void test_five(double T, double P){ Line 752  void test_five(double T, double P){
752                  , P, T_sat[i], rho_v[i]);                  , P, T_sat[i], rho_v[i]);
753      }      }
754    
755      p_d = dew_pressure(MS_vle, T, &err);      for(i=0;i<TRIALS;i++){
756      printf("\n  Dew pressure is %.0f Pa\n", p_d);          p_d[i] = dew_pressure(MS_vle, Ts[i], &err);
757            p_b[i] = bubble_pressure(MS_vle, Ts[i], &err);
758        }
759        /* p_d = dew_pressure(MS_vle, T, &err);
760        p_b = bubble_pressure(MS_vle, T, &err);
761        printf("\n  Dew pressure is %.0f Pa"
762                "\n  Bubble pressure is %.0f Pa"
763                , p_d, p_b);
764        puts(""); */
765        for(i=0;i<TRIALS;i++){
766            printf("\n  At temperature %.2f K"
767                    "\n\tDew pressure is %.0f Pa"
768                    "\n\tBubble pressure is %.0f Pa"
769                    , Ts[i], p_d[i], p_b[i]);
770        }
771        puts("");
772  #undef SPF  #undef SPF
773  #undef SXS  #undef SXS
774  #undef SPURE  #undef SPURE
# Line 1217  void test_six(void){ Line 791  void test_six(void){
791      char *fluids[] = {      char *fluids[] = {
792          "isohexane", "ammonia"          "isohexane", "ammonia"
793      };      };
794      char *fluid_names[] = {      /* char *fluid_names[] = {
795          "isohexane", "ammonia"          "isohexane", "ammonia"
796      };      }; */
797      double Xs[] = {0.55, 0.20, 0.10, 0.15, 0.10};      double Xs[] = {0.55, 0.20, 0.10, 0.15, 0.10};
798      char *src[] = {      char *src[] = {
799          NULL, NULL, NULL, NULL, NULL          NULL, NULL, NULL, NULL, NULL
# Line 1234  void test_six(void){ Line 808  void test_six(void){
808      mixture_fluid_spec(MS, NPURE, (void *)fluids, "pengrob", src, &merr);      mixture_fluid_spec(MS, NPURE, (void *)fluids, "pengrob", src, &merr);
809    
810      double P[NROWS][NCOLS];      double P[NROWS][NCOLS];
811      char *heads[NCOLS],      char *heads[NCOLS]
812           *sides[NROWS+1],          , *sides[NROWS+1]
813           *forms[NROWS] = {"%.6g"},          , *forms[NROWS] = {"%.6g"}
814           *conts[NROWS+1][NCOLS+1];          , *conts[NROWS+1][NCOLS+1];
815      unsigned widths[NCOLS],      unsigned widths[NCOLS],
816               temp_len; /* temporary string length */               temp_len; /* temporary string length */
817      sides[0] = "";      sides[0] = "";
# Line 1269  void test_six(void){ Line 843  void test_six(void){
843              err = FPROPS_NO_ERROR;              err = FPROPS_NO_ERROR;
844          }          }
845      }      }
     /* puts("");  
     for(i1=0;i1<NROWS+1;i1++){  
         printf("\n  -%s", sides[i1]);  
     }  
     puts("");  
     for(i1=0;i1<NCOLS;i1++){  
         printf("\n  %s", heads[i1]);  
     }  
     puts(""); */  
846    
847      /* PREPARE_TABLE(NROWS+1,NCOLS+1,heads,sides,P,forms,conts); */      /* PREPARE_TABLE(NROWS+1,NCOLS+1,heads,sides,P,forms,conts); */
848      for(i1=0;i1<NROWS;i1++){      for(i1=0;i1<NROWS;i1++){
# Line 1292  void test_six(void){ Line 857  void test_six(void){
857  #undef D  #undef D
858  #undef NPURE  #undef NPURE
859  }  }
860    
861    void test_seven(void){
862        unsigned j, i;
863    #define FUNC_REPT(REPORT) printf("\n  <test_seven>: %s", REPORT)
864    #define NPURE 4
865    #define D MS->PF[i]->data
866    #define TEMPS 5
867    #define PRESSURES 1
868        char *fluids[] = {
869            "isohexane", "krypton", "carbonmonoxide", "ammonia", "water"
870        };
871        char *fluid_names[] = {
872            "isohexane", "krypton", "carbon monoxide", "ammonia", "water"
873        };
874        double props[] = {11, 4, 2, 3, 2};
875        double Xs[NPURE];
876        mixture_x_props(NPURE, Xs, props);
877        char *src[] = {
878            NULL, NULL, NULL, NULL, NULL
879        };
880        FpropsError err = FPROPS_NO_ERROR;
881        MixtureError merr = MIXTURE_NO_ERROR;
882    
883        MixtureSpec *MS = ASC_NEW(MixtureSpec);
884    
885        MixtureSpec **MS_vle = ASC_NEW_ARRAY(MixtureSpec *,TEMPS);
886        MixtureSpec **MS_sc = ASC_NEW_ARRAY(MixtureSpec *,TEMPS);
887    #define VPURE(IX) MS_vle[ IX ]->pures
888    #define VXS(IX)   MS_vle[ IX ]->Xs
889    #define VPF(IX)   MS_vle[ IX ]->PF
890    #define SPURE(IX) MS_sc[ IX ]->pures
891    #define SXS(IX)   MS_sc[ IX ]->Xs
892    #define SPF(IX)   MS_sc[ IX ]->PF
893    
894        PhaseSpec **PS = ASC_NEW_ARRAY(PhaseSpec *,TEMPS);
895    #define NPHASE(IX) PS[ IX ]->phases
896    #define PNAME(IX) PS[ IX ]->ph_name
897    #define PFRAC(IX) PS[ IX ]->ph_frac
898    
899        MS->pures = NPURE;
900        MS->Xs = Xs;
901        MS->PF = ASC_NEW_ARRAY(PureFluid *,NPURE);
902        mixture_fluid_spec(MS, NPURE, (void *)fluids, "pengrob", src, &merr);
903    
904        double tol = 1.e-7
905            , p_b[TEMPS] = {0.0}
906            , p_d[TEMPS] = {0.0}
907            , xs_vle[NPURE]
908            , xs_sc[NPURE]
909            , phi_l[TEMPS][NPURE]
910            , phi_v[TEMPS][NPURE]
911            , rho_d[] = {0.0, 0.0}
912            , p_sat[TEMPS][NPURE]
913            , pf_l[TEMPS][NPURE]
914            , Flash_sum[TEMPS] = {0} /* sum of all vapor-phase minus all liquid-phase mole fractions in VLE */
915            ;
916        double Ts[] = {270, 310, 350, 390, 430, 470}
917            , P = 1.5e5;
918            ;
919    
920        FUNC_REPT("Declared all variables...");
921    
922        for(j=0;j<TEMPS;j++){
923            printf("\n  <test_seven>: %s %.2f %s",
924                    "Preparing all variables for temperture", Ts[j], "K...\n");
925    
926            MS_vle[j] = ASC_NEW(MixtureSpec);
927            VPURE(j) = 0;
928            VXS(j)   = ASC_NEW_ARRAY(double,NPURE);
929            VPF(j)   = ASC_NEW_ARRAY(PureFluid *,NPURE);
930    
931            MS_sc[j] = ASC_NEW(MixtureSpec);
932            SPURE(j) = 0;
933            SXS(j)   = ASC_NEW_ARRAY(double,NPURE);
934            SPF(j)   = ASC_NEW_ARRAY(PureFluid *,NPURE);
935    
936            /* NPHASE(j) = 0; */
937            /* PNAME(j) = ASC_NEW_ARRAY(PhaseName,3); */
938            /* PFRAC(j) = ASC_NEW_ARRAY(double,3); */
939    
940            FUNC_REPT("Prepared all variables");
941    
942            for(i=0;i<NPURE;i++){
943                if(Ts[j]<D->T_c){
944                    VPF(j)[VPURE(j)] = MS->PF[i];
945                    VXS(j)[VPURE(j)] = MS->Xs[i];
946                    VPURE(j) ++;
947                }else{
948                    /* NPHASE(j) = 1; */
949                    /* PNAME(j)[0] = SUPERCRIT; */
950                    SPF(j)[SPURE(j)] = MS->PF[i];
951                    SXS(j)[SPURE(j)] = MS->Xs[i];
952                    SPURE(j) ++;
953                }
954            }
955            mixture_x_props(VPURE(j), VXS(j), VXS(j));
956            mixture_x_props(SPURE(j), SXS(j), SXS(j));
957            mole_fractions(VPURE(j), xs_vle, VXS(j), VPF(j));
958            mole_fractions(SPURE(j), xs_sc, SXS(j), SPF(j));
959            p_b[j] = bubble_pressure(MS_vle[j], Ts[j], &err);
960            p_d[j] = dew_pressure(MS_vle[j], Ts[j], &err);
961    
962            if(P < p_b[j] && P > p_d[j]){
963                /* NPHASE(j) = 3; */
964                /* PNAME(j)[1] = VAPOR; */
965                /* PNAME(j)[2] = LIQUID; */
966            }
967    
968            for(i=0;i<VPURE(j);i++){
969                fprops_sat_T(Ts[j], (p_sat[j]+i), (rho_d), (rho_d+1), VPF(j)[i], &err);
970    
971                phi_l[j][i] = pengrob_phi_pure(VPF(j)[i], Ts[j], p_sat[j][i], LIQUID, &err);
972                /* phi_l[j][i] = pengrob_phi_pure(VPF(j)[i], Ts[j], p_b[j], LIQUID, &err); */
973                pf_l[j][i] = poynting_factor(VPF(j)[i], Ts[j], p_b[j], &err);
974    
975                phi_v[j][i] = pengrob_phi_pure(VPF(j)[i], Ts[j], p_b[j], VAPOR, &err);
976            }
977    
978            if(NPHASE(j)==3){
979                /* Flash_sum */
980            }
981            printf("\n  <test_seven>: %s %.2f K",
982                    "Found mixture conditions (mass and mole fractions, bubble and "
983                    "dew pressures) at temperature", Ts[j]);
984        }
985        for(j=0;j<TEMPS;j++){
986            printf("\n  At temperature %.2f K"
987                    "\n\tDew pressure is %.0f Pa"
988                    "\n\tBubble pressure is %.0f Pa"
989                    "\n\tAt bubble pressure --"
990                    , Ts[j], p_d[j], p_b[j]);
991            for(i=0;i<VPURE(j);i++){
992                printf("\n\t  liquid-phase fugacity coefficient for substance %s is %.6g"
993                        "\n\t  liquid-phase Poynting Factor for %s is %.6g"
994                        "\n\t  vapor-phase fugacity coefficient for %s is %.6g"
995                        "\n\t  overall vapor/liquid mass fraction for %s is %.6g"
996                        /* "\n\t  overall vapor/liquid mole fraction for %s is %.6g" */
997                        "\n"
998                        , VPF(j)[i]->name, phi_l[j][i]
999                        , VPF(j)[i]->name, pf_l[j][i]
1000                        , VPF(j)[i]->name, phi_v[j][i]
1001                        , VPF(j)[i]->name, VXS(j)[i]
1002                        /* , VPF(j)[i]->name, pf_l[j][i] */
1003                        );
1004            }
1005        } puts("");
1006    
1007    #if 0
1008        double c_coefs[] = {1, -4.635, -0.956902, 14.066841936}
1009            , c_roots[] = {0, 1, 2}
1010            ;
1011        unsigned c_num;
1012        c_num = cubic_solution(c_coefs, c_roots);
1013        printf("\n  <test_seven>: The equation x^3 %+.6g x^2 %+.6g x %+.6g has %u roots:"
1014                , c_coefs[1], c_coefs[2], c_coefs[3], c_num);
1015        for(i=0;i<c_num;i++){
1016            printf("\n\t%.6g", c_roots[i]);
1017        }
1018        puts("");
1019    #endif
1020    
1021    #undef SPF
1022    #undef SXS
1023    #undef SPURE
1024    #undef VPF
1025    #undef VXS
1026    #undef VPURE
1027    #undef PRESSURES
1028    #undef TEMPS
1029    #undef D
1030    #undef NPURE
1031    #undef FUNC_REPT
1032    }

Legend:
Removed from v.3022  
changed lines
  Added in v.3023

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