/* bwrm12.c by Satoru MOMOKI, Jan.1996 */ #include #include "bwrm.h" #define EPS0 1.0e-6 /* convergence limit */ #define EPS1 1.0e-6 #define DELTA0 1.0e-5 #define DELTA1 1.0e-5 #define MESSAGE (*m).b.tmp.errorcode #define N 2 double bwrmTyv_Pyl(/* Calculates T[K] & yv[mol/mol] from P[Pa] and yl[mol/mol] of bibary bibary mixtures at saturated condition. Return value is T[K]. And yv[mol/mol], rhol[kg-mol/m3] and rhov are set on (*c).v.y[1], (*c).prop_l.rho, and (*c).v.rho. Result of this function with yb as yl indicates properties at bubble points. */ BWRM2 *m, double p, /* pressure [Pa] */ double yl /* molar fraction in liquid [mol/mol] */ ) { double Pc; double ptmp; double yv, t; Pc = (*m).c[1].Pc * yl + (*m).c[2].Pc * (1.0 - yl); ptmp = p/Pc; if (ptmp > 0.98) ptmp = 0.98; t = yl * bwrmTs(&((*m).c[1]), (ptmp)*(*m).c[1].Pc ) + (1.0-yl) * bwrmTs(&((*m).c[2]), (ptmp)*(*m).c[2].Pc); if ((*m).c[1].Tc < (*m).c[2].Tc) yv = yl + (1.0 -yl) * 0.2; else yv = 0.9 * yl; return(bwrmTyv_Pyl2(m, p, yl, t, yv)); } double bwrmTyv_Pyl2( BWRM2 *m, double p, /* pressure [Pa] */ double yl, /* molar fraction in liquid [mol/mol] */ double t0, double yv0 ) #define t X[0] #define yv X[1] { BWRM2 Xp[N], Xm[N]; /* Xp = x + DELTAXi, Xm = x - DELTAXi */ double A[N][N], B[N], X[N], dX[N]; double absDX[N], Pc, f; double ymin, ymax, Tmin, Tmax; int nit=0; /* nit : iterative count */ int i; MESSAGE = BWRMERR_NOPROBLEM; (*m).v.rho = 0.0; (*m).l.rho = 30.0; if (yv0 < yl){ymin = 0.0; ymax = yl;} else {ymin = yl; ymax = 1.0;} Tmin = 230.0; /* --- Pure substance ? --- */ if (yl == 0.0){ t = bwrmTs(&(m->c[2]), p); CALC_F0F1_BWRM2(m, p, t, yl, yl); return t; } if (yl == 1.0){ t = bwrmTs(&(m->c[1]), p); CALC_F0F1_BWRM2(m, p, t, yl, yl); return t; } t = t0; if ( t < 250.0) t = 250.0; if ( t > (*m).c[1].Tc) t = (*m).c[1].Tc; if ( t > (*m).c[2].Tc) t = (*m).c[2].Tc; yv = yv0; for(i = 0; i < N; i++) {BWRM2cpy(&Xp[i], m); BWRM2cpy(&Xm[i], m);} for(i = 0; i < N; i++) {dX[i] = 1.2 * EPS0;} do{/* ----------------- BEGIN : Newton method ------------------ */ /* printf("%d %f %f %f : \n", nit, t, yv, yl); /**/ CALC_F0F1_BWRM2(m, p, t, yv, yl); CALC_F0F1_BWRM2(&Xp[0], p, t+DELTA0, yv, yl); CALC_F0F1_BWRM2(&Xm[0], p, t-DELTA0, yv, yl); CALC_F0F1_BWRM2(&Xp[1], p, t, yv+DELTA1, yl); CALC_F0F1_BWRM2(&Xm[1], p, t, yv-DELTA1, yl); B[0] = F0_BWRM2(m); B[1] = F1_BWRM2(m); for(i = 0; i < N; i++){ A[0][i] = (F0_BWRM2(&Xp[i]) - F0_BWRM2(&Xm[i])) / (2.0*DELTA0);/* dF0/dXi */ A[1][i] = (F1_BWRM2(&Xp[i]) - F1_BWRM2(&Xm[i])) / (2.0*DELTA1);/* dF1/dXi */ } MESSAGE = bwrmSolveA22x_b(A, B, dX); if ( MESSAGE != 0 ){MESSAGE = BWRMERR_CANT_RESOLVE_MATRIX; break;} if ( nit++ > 100 ) {MESSAGE = BWRMERR_CANT_CONVERGE; break;} /* cal f which is a factor to avoid y<0 or y>1 */ f = 1.0; { double z0, z1, zf; zf = 1.0; /* t */ z0 = t; z1 = t - dX[0]; if ( z1 <= Tmin) zf = (z0 - Tmin)/2.0 / (z0 - Tmin); /* yv */ z0 = yv; z1 = yv - dX[1]; if ( z1 <= ymin) zf = (z0 - ymin)/2.0 / (z0 - ymin); if ( z1 >= ymax) zf = (ymax - z0 )/2.0 / (ymax - z0); /* */ if ( zf < f) f = zf; } /* printf("%f %f %f %f %f %f \n", f, yv, yl, ymin, ymax, yv0);/**/ for(i=0; i < N; i++){ absDX[i] = fabs(dX[i]); X[i] -= f * dX[i]; } if (f < EPS0){ MESSAGE = BWRMERR_CANT_CONVERGE; break;} /* printf("--\n", t, yv, yl); /**/ }while( absDX[0] > EPS0 || absDX[1] > EPS1 ); /* --------------------- END : Newton method -------------------- */ switch(MESSAGE){ case BWRMERR_CANT_CONVERGE: bwrm2ERROR(stdout, "bwrTyv_Pyl warning: not convergenced", p, t, yv, yl); t = -1e+10; break; case BWRMERR_CANT_RESOLVE_MATRIX: bwrm2ERROR(stderr, "Fatal error : bwrmTyv_Pyl", p, t, yv, yl); bwrm2ERROR(stdout, "Fatal error : bwrmTyv_Pyl", p, t, yv, yl); break; default: MESSAGE = BWRMERR_NOPROBLEM; if (yv < 0.0) yv = 0.0; if (yv > 1.0) yv = 1.0; CALC_F0F1_BWRM2(m, p, t, yv, yl); /* sets yv, yl, rhov, rhol e.t.c. on 'c' */ } (*m).b.tmp.n_iteration = nit; /* printf("aaaa %d %f %f %f : \n", nit, t, yv, yl); /**/ return(t); } #undef yv #undef t double bwrmTyl_Pyv(/* Calculates T[K] & yl[mol/mol] from P[Pa] and yv[mol/mol] of binary mixtures at saturated condition. Return value is T[K]. And yl[mol/mol], rhol[kg-mol/m3] and rhov are set on (*c).l.y[1], (*c).prop_l.rho, and (*c).v.rho. Result of this function with yb as yl indicates properties at dew points. */ BWRM2 *m, double p, /* pressure [Pa] */ double yv /* molar fraction in liquid [mol/mol] */ ) { double Pc; double t, yl; Pc = (*m).c[1].Pc * yv + (*m).c[2].Pc * (1.0 - yv); t = yv * bwrmTs(&((*m).c[1]), (p/Pc)*(*m).c[1].Pc ) + (1.0-yv) * bwrmTs(&((*m).c[2]), (p/Pc)*(*m).c[2].Pc); if ((*m).c[1].Tc < (*m).c[2].Tc) yl = 0.9 * yv; else yl = yv + (1.0 - yv) * 0.2; return(bwrmTyl_Pyv2(m, p, yv, t, yl)); } double bwrmTyl_Pyv2( BWRM2 *m, double p, /* pressure [Pa] */ double yv, /* molar fraction in liquid [mol/mol] */ double t0, double yl0 ) #define t X[0] #define yl X[1] { BWRM2 Xp[N], Xm[N]; /* Xp = x + DELTAXi, Xm = x - DELTAXi */ double A[N][N], B[N], X[N], dX[N]; double absDX[N], Pc, f, ymin, ymax, Tmin, Tmax; int nit=0; /* nit : iterative count */ int i; MESSAGE = BWRMERR_NOPROBLEM; (*m).v.rho = 0.0; (*m).l.rho = 30.0; if (yl0 < yv){ymin = 0.0; ymax = yv;} else {ymin = yv; ymax = 1.0;} Tmin = 230.0; /* --- Pure substance ? --- */ if (yv == 0.0){ t = bwrmTs(&(m->c[2]), p); CALC_F0F1_BWRM2(m, p, t, yv, yv); return t; } if (yv == 1.0){ t = bwrmTs(&(m->c[1]), p); CALC_F0F1_BWRM2(m, p, t, yv, yv); return t; } t = t0; if ( t < 250.0) t = 250.0; if ( t > (*m).c[1].Tc) t = (*m).c[1].Tc; if ( t > (*m).c[2].Tc) t = (*m).c[2].Tc; yl = yl0; for(i = 0; i < N; i++) {BWRM2cpy(&Xp[i], m); BWRM2cpy(&Xm[i], m);} for(i = 0; i < N; i++) {dX[i] = 1.2 * EPS0;} do{/* ----------------- BEGIN : Newton method ------------------ */ /* printf("%d %f %f %f : \n", nit, t, yv, yl); /**/ CALC_F0F1_BWRM2(m, p, t, yv, yl); CALC_F0F1_BWRM2(&Xp[0], p, t+DELTA0, yv, yl); CALC_F0F1_BWRM2(&Xm[0], p, t-DELTA0, yv, yl); CALC_F0F1_BWRM2(&Xp[1], p, t, yv, yl+DELTA1); CALC_F0F1_BWRM2(&Xm[1], p, t, yv, yl-DELTA1); B[0] = F0_BWRM2(m); B[1] = F1_BWRM2(m); for(i = 0; i < N; i++){ A[0][i] = (F0_BWRM2(&Xp[i]) - F0_BWRM2(&Xm[i])) / (2.0*DELTA0);/* dF0/dXi */ A[1][i] = (F1_BWRM2(&Xp[i]) - F1_BWRM2(&Xm[i])) / (2.0*DELTA1);/* dF1/dXi */ } MESSAGE = bwrmSolveA22x_b(A, B, dX); if ( MESSAGE != 0 ){MESSAGE = BWRMERR_CANT_RESOLVE_MATRIX; break;} if ( nit++ > 100 ) {MESSAGE = BWRMERR_CANT_CONVERGE; break;} /* cal f which is a factor to avoid y<0 or y>1 */ f = 1.0; { double z0, z1, zf; zf = 1.0; /* t */ z0 = t; z1 = t - dX[0]; if ( z1 <= Tmin) zf = (z0 - Tmin)/2.0 / (z0 - Tmin); /* yv */ z0 = yl; z1 = yl - dX[1]; if ( z1 <= ymin) zf = (z0 - ymin)/2.0 / (z0 - ymin); if ( z1 >= ymax) zf = (ymax - z0 )/2.0 / (ymax - z0); /* */ if ( zf < f) f = zf; } /* printf("%f %f %f %f %f %f %f \n", f, t, yv, yl, ymin, ymax, yl0);/**/ for(i=0; i < N; i++){ X[i] -= f * dX[i]; absDX[i] = fabs(dX[i]); } if (f < EPS0){ MESSAGE = BWRMERR_CANT_CONVERGE; break;} }while( absDX[0] > EPS0 || absDX[1] > EPS1 ); /* printf("%d %f %f %f : \n", nit, t, yv, yl); /**/ /* --------------------- END : Newton method -------------------- */ switch(MESSAGE){ case BWRMERR_CANT_CONVERGE: bwrm2ERROR(stdout, "bwrTyl_Pyv warning: not convergenced", p, t, yv, yl); t = -1e+10; break; case BWRMERR_CANT_RESOLVE_MATRIX: bwrm2ERROR(stderr, "Fatal error : bwrmTyl_Pyv", p, t, yv, yl); bwrm2ERROR(stdout, "Fatal error : bwrmTyl_Pyv", p, t, yv, yl); break; default: if (yl < 0.0) yl = 0.0; if (yl > 1.0) yl = 1.0; CALC_F0F1_BWRM2(m, p, t, yv, yl); /* sets yv, yl, rhov, rhol e.t.c. on 'c' */ } (*m).b.tmp.n_iteration = nit; return(t); #undef MESSAGE } #undef t #undef yl