/* bwrm7.c */ #include "bwrm.h" #include #include #define CHANGED 0 #define NOT_CHANGED 1 #define CRITICAL 2 #define NOPROBREM ((*c).tmp.errorcode == BWRMERR_NOPROBLEM) double RhoForbwrmPs( BWRM *c, double *p, double *pmax, double *pmin, double T, double rho0, int *rcode); double bwrmPs(/* calculates saturation pressure in [Pa] from T[K] */ BWRM *c, double T /* temperature [K] */ ){ #define TMP (*c).tmp int i, N; double A[2][2], b[2], x[2]; double F1, F2, dF1dv, dF1dl, dF2dv, dF2dl; double rhov, rhol, p, pmax, pmin; double D, EPS; int flag; EPS=1e-6; D=1e-8; N=1000; /* EPS: convergence limit, D: used in calc. dF/dho, N: max. of iteration */ /* --- T >= Tc ?? --- */ TMP.n_iteration = 0; TMP.errorcode = BWRMERR_NOPROBLEM; if (T > (*c).Tc ) TMP.errorcode = BWRMERR_ARGUMENT_OUTOFRANGE; if (T >= (*c).Tc) return((*c).Pc); /* ------ Out of range ?? ------------ */ if ( bwrmdPdrho(c, T, (*c).RHOc) > 0.0){ /* Hmmm! Tr is out of range where the BWR equation guarantees the accuracy. */ (*c).rhov = (*c).rhol = (*c).RHOc; (*c).P = bwrmP(c, T, (*c).RHOc); TMP.n_iteration = 1; TMP.errorcode = BWRMERR_INACCURATE_RESULT; return((*c).P); } /* --- determination initial values on rhov and rhol --- */ p = (*c).P = (*c).Pc * pow(10.0, -7.0/3.0 * (1.0+(*c).omega) * ((*c).Tc/T-1.0)); /* see SAITO SHOZABURO, HEIKOUBUSSEI-SUISAN-NO-KISO, pp.100-101 */ pmax = (*c).Pc; pmin = 0.0; rhov = 0.0001; rhol = 30.0; (*c).RHOc = (*c).RHOc; do{ flag = NOT_CHANGED; rhov = RhoForbwrmPs(c, &p, &pmax, &pmin, T, rhov, &flag); rhol = RhoForbwrmPs(c, &p, &pmax, &pmin, T, rhol, &flag); }while(flag == CHANGED); rhov -= 1.1 * D; rhol += 1.1 * D; TMP.n_iteration = 0; TMP.errorcode = BWRMERR_NOPROBLEM; /* --- begin iteration --- */ for (i = 0; i < N ; i++ ){ double F1pv, F2pv, rhov_prev, rhol_prev; rhov_prev = rhov; rhol_prev = rhol; F1 = bwrmP(c, T, rhov) - bwrmP(c, T, rhol); /* = Pv - Pl */ F2 = bwrmRTlnF(c, T, rhov) - bwrmRTlnF(c, T, rhol); /* = fv - fl */ dF1dv = bwrmdPdrho(c, T, rhov); dF1dl = -bwrmdPdrho(c, T, rhol); #define dF2(d) ((bwrmRTlnF(c, T, (d)+D) - bwrmRTlnF(c, T, (d)-D))/2.0/D) dF2dv = dF2(rhov); dF2dl = -dF2(rhol); #undef dF2 A[0][0]=dF1dv; A[0][1]=dF1dl; b[0]=F1; A[1][0]=dF2dv; A[1][1]=dF2dl; b[1]=F2; if ( bwrmSolveA22x_b(A, b, x) != 0 ){ /* printf("Fatal error : bwrmPs at T=%lfK\n", T); /**/ TMP.errorcode = BWRMERR_CANT_RESOLVE_MATRIX; goto RETURN; } rhov -= x[0]; rhol -= x[1]; if (rhov <= D) /* too small */ rhov = rhov_prev * 0.8; if (x[1] <= (*c).RHOc * 0.2) /* too large */ rhol = rhol_prev + (*c).RHOc * 0.2; {/* BEGIN: check rhov and rhol ---------------------------- */ double coef, dv, dl; coef = 0.5; while(rhov+D >= (*c).RHOc || bwrmdPdrho(c, T, rhov+D) < 0.0){ rhov = rhov_prev - x[0] * coef; coef *= 0.5; } coef = 0.5; while(rhol-D <= (*c).RHOc || bwrmdPdrho(c, T, rhol-D) < 0.0){ rhol = rhol_prev - x[1] * coef; coef *= 0.5; } dv = fabs((rhov - rhov_prev)/x[0]); dl = fabs((rhol - rhol_prev)/x[1]); if ( dv !=1 || dl !=1 ) if (dv < dl) rhol = rhol_prev - x[1] * dv; else rhov = rhov_prev - x[0] * dl; }/* END: check rhov and rhol ---------------------------- */ if ((fabs(x[0]/rhov) <= EPS) && (fabs(x[1]/rhol) <= EPS)) goto RETURN; /* printf("%E %E %E\n", rhov, rhol, bwrmP(c, T, rhov)); /**/ } /* printf("bwrmPs warning: not convergenced(T=%lfK)\n", T);/**/ TMP.errorcode = BWRMERR_CANT_CONVERGE ; RETURN: (*c).rhov = rhov; (*c).rhol = rhol; (*c).P = bwrmP(c, T, rhov); TMP.n_iteration = i; return((*c).P); #undef TMP } double RhoForbwrmPs( BWRM *c, double *p, double *pmax, double *pmin, double T, double rho0, int *rcode) { #define VAPOR 1 #define LIQUID 2 int mode; double rho, rfactor, factor; double rhomin, rhomax, Padd; int i,n; double p0; p0 = *p; if ( rho0 < (*c).RHOc ) mode = VAPOR; else mode = LIQUID; if (mode == VAPOR ){ rhomin = 0.0; rhomax = (*c).RHOc; rfactor = 0.9; }else{ rhomin = (*c).RHOc; rhomax = 100.0; rfactor = 2.0; } rho = bwrmRHO(c, (*p), T, rfactor * rho0); Padd = 0.0; factor = 0.1; i = 0; while( !((rhomin < rho) && (rho < rhomax) && NOPROBREM && (bwrmdPdrho(c, T, rho) > 0.0) ) ){ if (i == 9){i = 0; factor *= 0.1;} i++; Padd += factor; if (mode == VAPOR ){ *pmax = *p; *p = p0 - Padd * (p0 - *pmin); }else{ *pmin = *p; *p = p0 + Padd * (*pmax - p0); } rho = bwrmRHO(c, (*p), T, rho0 * rfactor); /* fprintf(stdout, "- p = %E T %.3f RHO %9.2E %9.2E (%10.3E-%10.3E) -- %d\n", *p, T, rho, (*c).RHOc, rhomin, rhomax, *rcode); /**/ *rcode = CHANGED; } return(rho); } double bwrmTs(/* calculates saturation temperature in [K] from P. This use Newton method with bwrmPs(). */ BWRM *c, double p /* pressure[Pa] */ ){ #define EPS 1.0e-2 #define D 1.0e-3 #define TMP (*c).tmp double T, Tp, Tm, F; double buf; double y, dx, df; int n = 0; int error; TMP.n_iteration = 0; error = BWRMERR_NOPROBLEM; /* --- P >= Pc ?? --- */ if (p >= (*c).Pc) { TMP.errorcode = BWRMERR_ARGUMENT_OUTOFRANGE; return((*c).Tc); } T = (*c).Tc / ( 1 + log10(p/(*c).Pc)/(-7.0/3.0*(1.0+(*c).omega))); /* see SAITO SHOZABURO, HEIKOUBUSSEI-SUISAN-NO-KISO, pp.100-101 */ do{ Tp = bwrmPs(c, T+D); Tm = bwrmPs(c, T-D); df = (Tp - Tm) / (2.0 * D); if (df <= 0.0) df = 1e+8; F = bwrmPs(c, T) - p; dx = (F) / df; T -= dx; if (T >= (*c).Tc - D){ T += dx; T = (T + (*c).Tc)/2.0; } if (n++ > 1000){ /* printf("bwrmTs warning: not convergenced(P=%lfMPa)\n", p/1.0e+6);/**/ error = BWRMERR_CANT_CONVERGE ; goto RETURN; } }while(fabs(F) > EPS ); error = BWRMERR_NOPROBLEM; goto RETURN; RETURN: (*c).rhov = bwrmRHOv(c, p, T); (*c).rhol = bwrmRHOl(c, p, T); (*c).T = bwrmT(c, p, (*c).rhov, T); TMP.errorcode = error; TMP.n_iteration = n; return((*c).T); #undef TMP #undef EPS #undef D }