/********************************************************************* * bwrm97.c * * Interface program between my bwrm-libraly and PROPATH. #1 * * May, 1997 S.MOMOKI * *********************************************************************/ #include "propath.h" #include "bwrm.h" #include #include static BWRM Pu[3]; static BWRM2 Mix; struct { int mode; double A[10]; double Tmin; double Tmax; }CPID[3]; static char FNAME[30]; /* name of function */ static int ERR_DETECT; static char messagestr[80]; /* macros for the case of error */ #define ERROR_MESSAGE Mix.b.tmp.errorcode #define UNLESS_CONVERGED if (ERROR_MESSAGE == BWRMERR_CANT_CONVERGE) #define IF_BWRM_OUTOFRANGE if (ERROR_MESSAGE == BWRMERR_ARGUMENT_OUTOFRANGE) #define ERROR_MESSAGE_P(iii) Pu[iii].tmp.errorcode #define UNLESS_CONVERGED_P(iii) \ if (ERROR_MESSAGE_P(iii) == BWRMERR_CANT_CONVERGE) #define IF_BWRM_OUTOFRANGE_P(iii) \ if (ERROR_MESSAGE_P(iii) == BWRMERR_ARGUMENT_OUTOFRANGE) /* macros for unit-conversion */ #define IF_KG_UNIT if(PROPATH.KAS == 1) #define IF_C_UNIT if(PROPATH.KPA == 1 || PROPATH.KPA == 3) #define IF_BAR_UNIT if(PROPATH.KPA == 1 || PROPATH.KPA == 2) #define IF_OUTOFRANGE if(ERR_DETECT == 1) #define CHECK_Z(z) {ERR_DETECT = 0;\ IF_KG_UNIT {\ if (z < 0.0 || 1.0 < z){\ OUTOFRANGE_SUB(FNAME, "Z", z);\ ERR_DETECT = 1;}\ else{\ z = kg2molZ(z);}\ }else{\ if (z < 0.0 || 1.0 < z){\ OUTOFRANGE_SUB(FNAME, "Z", z);\ ERR_DETECT = 1;\ }\ }} #define CHECK_T(t) {ERR_DETECT = 0;\ IF_C_UNIT {t += 273.15;}\ if (t < 200.0 ){\ OUTOFRANGE_SUB(FNAME, "T[K]", t);\ ERR_DETECT = 1;}} #define CHECK_P(p) {ERR_DETECT = 0;\ if(PROPATH.KPA == 1 || PROPATH.KPA == 2){p *= 1.0e+5;}\ if (p <= 0.0) {\ OUTOFRANGE_SUB(FNAME, "P[Pa]", p);\ ERR_DETECT = 1;}} #define CHECK_V(v) {ERR_DETECT = 0;\ IF_KG_UNIT{v *= Mix.b.M;}\ if (v < 0.0) {\ OUTOFRANGE_SUB(FNAME, "V", v);\ ERR_DETECT = 1;}} #define CHECK_H(h) {ERR_DETECT = 0; IF_KG_UNIT{h *= Mix.b.M;}} #define CHECK_S(s) {ERR_DETECT = 0; IF_KG_UNIT{s *= Mix.b.M;}} #define RECHECK_Z(z) {IF_KG_UNIT {z = mol2kgZ(z);}} #define RECHECK_T(t) {IF_C_UNIT {t -= 273.15;}} #define RECHECK_P(p) {IF_BAR_UNIT {p /= 1.0e+5;}} #define RECHECK_V(v) {IF_KG_UNIT {v /= Mix.b.M;}} #define RECHECK_H(h) {IF_KG_UNIT {h /= Mix.b.M;}} #define RECHECK_S(s) {IF_KG_UNIT {s /= Mix.b.M;}} #define SUPER_CRITICAL_PRESSURE_M ( #define OVER_PMAX_M (P > Mix.c[1].Pc || P > Mix.c[2].Pc) /* probably occur seg. fault */ #define OVER_PMAX_P(zzz) (P >= Pu[zzz].Pc) /******************************************************************** OUTOFRANGE_SUB and NOCONVERGENCE_SUB These subroutines output error message onto standard output device. Why do they call F77PRT0x those are written in FORTRAN ? Because it is impossible to use printf() in C-program under the MSC-7 and MS-Fortran5.1 environment. ********************************************************************/ #ifdef __MS_C_ONLY subroutine F77PRT01(char *z1, char *z2, float *v); subroutine F77PRT02(char *z); #endif void OUTOFRANGE_SUB(char *fname, char *valname, double val) { if (PROPATH.MESS == 1){ char z1[30],z2[30]; float v; #ifdef __MS_C_ONLY strcpy(z1, fname); strcat(z1, " "); strcpy(z2, valname); strcat(z2, " "); /* do not use over 10 characters for *name */ v = val; F77PRT01(z1, z2, &v); #else printf("*** OUT OF RANGE AT %10s, (%5s = %10.2e) ***\n", fname, valname, val); #endif } } #define OUTOFRANGE_REAL(func, name, val) \ {OUTOFRANGE_SUB(func, name, val);return -1.0E+20;} #define OUTOFRANGE_INT(func, name, val) \ {OUTOFRANGE_SUB(func, name, val);return -2;} void NOCONVERGENCE_SUB(char *fname) { char zzz[30]; if (PROPATH.MESS == 1){ strcpy(zzz, fname);strcat(zzz, " "); /* do not use over 10 characters for fname */ #ifdef __MS_C_ONLY F77PRT02(zzz); #else printf("*** NO CONVERGENCE AT %10s ***\n", fname); #endif } } #define NOCONVERGENCE_REAL(func) \ {NOCONVERGENCE_SUB(func);return -1.0E+10;} #define NOCONVERGENCE_INT(func) \ {NOCONVERGENCE_SUB(func);return -1;} /********************************************************************* KPAMES and STANKAS *********************************************************************/ subroutine KPAMES(integer KPA, integer MESS) { /* KPA (system unit) 0: [K], [Pa] 1: [C], [bar] 2: [K], [bar] 3: [C], [Pa] Others : same to 0 */ PROPATH.KPA = *KPA; /* MESS (mode) 0: quiet mode Others: verbose mode */ PROPATH.MESS = *MESS; } subroutine STNKAS(integer KSTAN, integer KAS) { /* KSTAN (reference state) 0: ideal gas state of each pure components at 298.15K, 1bar 1: satulated liquid of each pure components at 0[C], 200[kJ/kg], 1[kJ/kg] */ if (*KSTAN == 0){ OUTOFRANGE_SUB("STANKAS", "KSTAN", *KSTAN); } PROPATH.KSTAN = *KSTAN; /* KAS (unit of ratio of quantity ) 0: mol/mol 1: kg/kg Others: same to 0 */ PROPATH.KAS = *KAS; } /********************************************************************* Cpid, Hid, Sid; *********************************************************************/ /* -------------------- Cpid ----------------------- */ #define A0 CPID[i].A[0] #define A1 CPID[i].A[1] #define A2 CPID[i].A[2] #define A3 CPID[i].A[3] #define A4 CPID[i].A[4] #define A5 CPID[i].A[5] double bwrmCPID(int i, double T){ double cp, T2, T3, T4; /* if(T < CPID[i].Tmin || T > CPID[i].Tmax) OUTOFRANGE_REAL("CPID", "T", T);/**/ T2 = T*T; T3 = T2*T; T4 = T3*T; switch(CPID[i].mode){ case 2: cp = A0 + A1*T + A2*T2 + A3*T3 + A4*T4; break; case 3: cp = A0/T2 + A1/T + A2 + A3*T + A4*T2 + A5*T3; break; deault: ; } return(cp); } double bwrmCPID1(double T){return(bwrmCPID(1, T));} double bwrmCPID2(double T){return(bwrmCPID(2, T));} /* -------------------- Hid ----------------------- */ double bwrmHID(int i, double T){ double h, T2, T3, T4, T5; /* if (T < CPID[i].Tmin || T > CPID[i].Tmax) OUTOFRANGE_REAL("HID", "T", T);/**/ T2 = T*T; T3 = T2*T; T4 = T3*T, T5 = T4*T; switch(CPID[i].mode){ case 2: h = A0*T+ A1*T2/2.0 + A2*T3/3.0 + A3*T4/4.0 + A4*T5/5.0; break; case 3: h = -A0/T + A1*log(T) + A2*T + A3*T2/2.0 + A4*T3/3.0 + A5*T4/4.0; break; deault: ; } return(h); } double bwrmHID1(double T){return(bwrmHID(1, T));} double bwrmHID2(double T){return(bwrmHID(2, T));} /* -------------------- Sid ----------------------- */ double bwrmSID(int i, double T){ double s, T2, T3, T4, T5; /* if (T < CPID[i].Tmin || T > CPID[i].Tmax)OUTOFRANGE_REAL("SID", "T", T);/**/ T2 = T*T; T3 = T2*T; T4 = T3*T, T5 = T4*T; switch(CPID[i].mode){ case 2: s = A0*log(T) + A1*T + A2*T2/2.0 + A3*T3/3.0 + A4*T4/4.0; break; case 3: s = -2.0*A0/T2 - A1/T + A2*log(T) + A3*T + A4*T2/2.0 + A5*T3/3.0; break; deaults: ; } return(s); } double bwrmSID1(double T){return(bwrmSID(1, T));} double bwrmSID2(double T){return(bwrmSID(2, T));} /* -------------------- mij constant ----------------------- */ double MIJ_CONST; double bwrmMij_const(double y, double T){ return MIJ_CONST; } /********************************************************************* START1 (Not provided yet) *********************************************************************/ subroutine START1( integer error, integer kombi ){ /* fprintf(stderr, "Sorry, START1 for NS-BWR equation is not provided yet.\n"); /**/ *error = -2; } /********************************************************************* START2 *********************************************************************/ subroutine START2( integer errcode, float PR1[], /* PR1(1):M, PR1(2):Tc, PR1(3):Pc, PR1(4):Vc, PR1(5):omega -- f77 */ real PR2, /* PR2[0]:M, PR2[1]:Tc, PR2[2]:Pc, PR2[3]:Vc, PR2[4]:omega -- C */ real CP1, /* CP1(1):TYPE, CP1(2):Tmin, CP1(3):Tmax, CP1(4)--CP1(10):A[0-6], */ real CP2, real MIJ ){ int i; double Mw, Tc, Pc, Vc, om; strcpy(FNAME,"START2"); *errcode = -2; i = (int)(CP1[0]); if(i < 1 || 3 < i) {OUTOFRANGE_SUB("START2", "CP1(1)", (double)i);return;} i = (int)(CP2[0]); if(i < 1 || 3 < i) {OUTOFRANGE_SUB("START2", "CP2(1)", (double)i);return;} *errcode = 0; /* initialize CPID */ CPID[1].mode = (int)(CP1[0]); CPID[1].Tmin = CP1[1]; CPID[1].Tmax = CP1[2]; CPID[2].mode = (int)(CP2[0]); CPID[2].Tmin = CP2[1]; CPID[2].Tmax = CP2[2]; for(i = 0; i < 6; i++){ CPID[1].A[i] = CP1[i+3]; CPID[2].A[i] = CP2[i+3]; } /* initialize P1, P2 */ Mw = PR1[0]; Tc = PR1[1]; Pc = PR1[2]; Vc = PR1[3]; om = PR1[4]; InitBWRM(&Pu[1], Tc, 1.0/Vc, Mw, bwrmOMEGA_CONST, om, 0.0, bwrmCPID1, bwrmHID1, bwrmSID1); Mw = PR2[0]; Tc = PR2[1]; Pc = PR2[2]; Vc = PR2[3]; om = PR2[4]; InitBWRM(&Pu[2], Tc, 1.0/Vc, Mw, bwrmOMEGA_CONST, om, 0.0, bwrmCPID2, bwrmHID2, bwrmSID2); if (*MIJ <= 0.0 || CPID[1].Tmin >= CPID[1].Tmax || CPID[2].Tmin >= CPID[2].Tmax){ /* printf ("**** Funny parameters in START2. Check your program again.\n"); /**/ *errcode = -2; return; } /* initialize Mix */ MIJ_CONST = (double)(*MIJ); InitBWRM2(&Mix, &(Pu[1]), &(Pu[2]), 0.5, bwrmMij_const); } /********************************************************************* AKG & AKMOL *********************************************************************/ double mol2kgZ(double z) { double y1, y2; if (z < 0.0 || z > 1.0) return z; y1 = z * Pu[1].M; y2 = (1.0 - z) * Pu[2].M; return(y1/(y1+y2)); } function_real AKG(real ZZ) { double z, y1, y2; strcpy(FNAME,"AKG"); z = (double)*ZZ; if (z < 0.0 || 1.0 < z) OUTOFRANGE_REAL("AKG", "Z", z); return(mol2kgZ(z)); } double kg2molZ(double z) { double y1, y2; if (z < 0.0 || z > 1.0) return z; y1 = z / Pu[1].M; y2 = (1.0 - z) / Pu[2].M; return(y1/(y1+y2)); } function_real AKMOL(real ZZ) { double z; strcpy(FNAME,"AKMOL"); z = (double)*ZZ; if (z < 0.0 || 1.0 < z) OUTOFRANGE_REAL("AKMOL", "Z", z); return(kg2molZ(z)); } /********************************************************************* CRPM999 and FCM999 It is too difficult for me to call VC++ function with charater- array-arguments from MS-Fortran PWB. CRPM and FCM are written in FORTRAN, and they call CRPM999 and FCM999. *********************************************************************/ function_real CRPM999(integer II, integer A) { int i; double r; i = *II; strcpy(FNAME,"CRPM"); if ( i < 1 || 2 < i) OUTOFRANGE_INT("CRPM","I", i); switch(*A){ case 1: /* T */ r = Pu[i].Tc; RECHECK_T(r); return r; case 2: /* 'P': */ r = Pu[i].Pc; RECHECK_P(r); return r; case 3: /* V */ r = 1.0 / Pu[i].RHOc; IF_KG_UNIT r /= Pu[i].M; return r; case 4: /* H */ r = bwrmH(&(Pu[i]), Pu[i].Tc, Pu[i].RHOc); IF_KG_UNIT r /= Pu[i].M; return r; case 5: /* S */ r = bwrmS(&(Pu[i]), Pu[i].Tc, Pu[i].RHOc); IF_KG_UNIT r /= Pu[i].M; return r; } } function_real FCM999(integer II, integer A) { int i; double r; i = *II; strcpy(FNAME,"FCM"); if ( i < 1 || 2 < i) OUTOFRANGE_INT("FCM","I", i); switch(*A){ case 1: /* 'T': */ r = Pu[i].Tc; RECHECK_T(r); return r; case 2: /* 'P': */ r = Pu[i].Pc; RECHECK_P(r); return r; case 3: /* 'V': */ r = 1.0 / Pu[i].RHOc; IF_KG_UNIT r /= Pu[i].M; return r; case 4: /* 'R': */ /* Universal gas constant is 8.31451e+3 /* [J/(kg-mole.K)] */ r = 8.31451e+3 / Pu[i].M; return r; case 5: /* 'M': */ r = Pu[i].M; return r; case 6: /* 'E': */ r = Pu[i].omega; return r; } } /********************************************************************* ZPHASE( IPHASE) The Combination of MSC-7 and MS-Fortran5.1 dose not work properly for integer function written in C. ZPHASE works as same as IPHASE but return 'real' value instead of integer value. IPHASE itself is written in FORTRAN(bwrm98.for) and it calls ZPHASE. *********************************************************************/ function_real ZPHASE(real TT, real PP, real ZZ) #define CHECK_ERR \ IF_BWRM_OUTOFRANGE {OUTOFRANGE_REAL(FNAME,"?", 0.0);}\ UNLESS_CONVERGED {NOCONVERGENCE_REAL(FNAME);} { int err; double T, P, Z, Tb, Td, yl; strcpy(FNAME,"IPHASE"); /* range check */ Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{ return -1e+20;} T = *TT; CHECK_T(T); IF_OUTOFRANGE{ return -1e+20;} P = *PP; CHECK_P(P); IF_OUTOFRANGE{ return -1e+20;} if (P > Mix.c[1].Pc && P > Mix.c[2].Pc ){/* Super critical pressure */ return(4.0); } if (OVER_PMAX_M) NOCONVERGENCE_INT("IPHASE"); Tb = bwrmTyv_Pyl(&Mix, P, Z); CHECK_ERR; if ( T < Tb ) return(1.0); /* Td = bwrmTyl_Pyv(&Mix, P, Z); CHECK_ERR; /**/ if (Mix.v.y[1] < Z) yl = (1+Z)/2.0; else yl = Z / 2.0; Td = bwrmTyl_Pyv2(&Mix, P, Z, Tb+5.0, yl); CHECK_ERR; /**/ if ( T > Td ) return(3.0); return(2.0); } #undef CHECK_ERR /********************************************************************* QMIX ********************************************************************/ #define CHECK_ERR \ IF_BWRM_OUTOFRANGE {OUTOFRANGE_SUB(FNAME,"P", P);}\ UNLESS_CONVERGED {NOCONVERGENCE_SUB(FNAME);} function_real QMIX(real TT, real PP, real ZZ) { double T, P, Z, yv, yl, yb, Tb, Td; strcpy(FNAME,"QMIX"); /* range check */ Z = *ZZ; CHECK_Z(Z);IF_OUTOFRANGE{ return -1.0e-20;} T = *TT; CHECK_T(T);IF_OUTOFRANGE{ return -1.0e-20;} P = *PP; CHECK_P(P);IF_OUTOFRANGE{ return -1.0e-20;} if (OVER_PMAX_M) OUTOFRANGE_REAL(FNAME, "P[Pa]", P); /* for pure substance */ if (Z == 0.0 || Z == 1.0){ double ts; int i; i = 2 - (int)Z; /* Z=0 -> i = 2, Z=1 -> i = 1; /**/ ts = bwrmTs(&(Mix.c[i]), P); if (T > ts) return(1.0); if (T < ts) return(0.0); NOCONVERGENCE_REAL(FNAME); } /* define initial value */ if (Pu[1].Tc < Pu[2].Tc) { yl = 0.8 * Z; yv = Z + (1.0 - Z) * 0.2; }else{ yv = 0.8 * Z; yl = Z + (1.0 - Z) * 0.8; } /* range check */ Tb = bwrmTyv_Pyl(&Mix, P, Z); CHECK_ERR; if (T <= Tb) return 0.0; /* Td = bwrmTyl_Pyv(&Mix, P, Z); CHECK_ERR; /**/ if (Mix.v.y[1] < Z) yl = (1+Z)/2.0; else yl = Z / 2.0; Td = bwrmTyl_Pyv2(&Mix, P, Z, Tb+5.0, yl); CHECK_ERR; if (T >= Td) return 1.0; /* calc. properties at saturated condition */ bwrmYvYl_PT(&Mix, P, T, yv, yl); UNLESS_CONVERGED{ NOCONVERGENCE_REAL("QMIX"); } if ( Z < Mix.l.y[1]) return 0.0; if ( Z > Mix.v.y[1]) return 1.0; yb = Z; yl = Mix.l.y[1]; yv = Mix.v.y[1]; IF_KG_UNIT{ yb = mol2kgZ(yb); yl = mol2kgZ(yl); yv = mol2kgZ(yv); } return((yb - yl)/(yv - yl)); } #undef CHECK_ERR /********************************************************************* SUBMIX ********************************************************************/ subroutine SUBMIX( integer II, integer JJ, real TT, real PP, real ZZ, real VV, real HH, real SS ) #define CHECK_ERR \ IF_BWRM_OUTOFRANGE {OUTOFRANGE_SUB(FNAME,"?", 0.0); *JJ = -2;return;}\ UNLESS_CONVERGED {NOCONVERGENCE_SUB(FNAME); *JJ = -1;return;} { double T, P, Z, V, H, S; double Tb, Td, rho, rhov, rhol, hv, hl, h, sv, sl, s, x, v, yv, yl; strcpy(FNAME,"SUBMIX"); T = *TT; P = *PP; Z = *ZZ; H = *HH; S = *SS; V = *VV; *JJ = 0; switch(*II){ /* --------------------------------------------------------------- II = 1 (P,T,Z -> H,V,S) --------------------------------------------------------------- */ case 1: CHECK_P(P); IF_OUTOFRANGE{*JJ = -2;return;} CHECK_T(T); IF_OUTOFRANGE{*JJ = -2;return;} CHECK_Z(Z); IF_OUTOFRANGE{*JJ = -2;return;} if (OVER_PMAX_M){ OUTOFRANGE_SUB(FNAME,"P", P); *JJ = -2;return; } /* ---------------- Liquid ? --------------- */ Tb = bwrmTyv_Pyl(&Mix, P, Z); CHECK_ERR; yv = Mix.v.y[1]; /* use to determin intial value */ if ( T < Tb ) { rhol = bwrmRHOl(&(Mix.l), P, T); h = hl = bwrmH2(&(Mix), 'L', T, rhol, Z); s = sl = bwrmS2(&(Mix), 'L', T, rhol, Z); v = 1.0 /rhol; Mix.b = Mix.l; goto Return_1; } /* ---------------- Gas ? --------------- */ /* Td = bwrmTyl_Pyv(&Mix, P, Z); CHECK_ERR; /**/ if (Mix.v.y[1] < Z) yl = (1+Z) / 2.0; else yl = Z / 2.0; Td = bwrmTyl_Pyv2(&Mix, P, Z, Tb+5.0, yl); CHECK_ERR;/**/ yl = Mix.l.y[1]; /* use to determine intial value */ if ( T > Td ){ rhov = bwrmRHOv(&(Mix.v), P, T); h = hv = bwrmH2(&(Mix), 'V', T, rhov, Z); s = sv = bwrmS2(&(Mix), 'V', T, rhov, Z); v = 1.0 /rhov; Mix.b = Mix.v; goto Return_1; } /* ---------------- VLE --------------- */ if (Z == 0.0 || Z == 1.0){ OUTOFRANGE_SUB(FNAME,"Z", Z); *JJ = -2;return; } /* calc. equivalent point */ x = (T - Tb) / (Td - Tb); yl = Z + x * (yl - Z); yv = yv - x * (yv - Z); if (Pu[1].Tc > Pu[2].Tc){x = yl; yl = yv; yv = x;} bwrmYvYl_PT(&Mix, P, T, yv, yl); IF_BWRM_OUTOFRANGE {OUTOFRANGE_SUB(FNAME,"?", 0.0); *JJ = -2;return;} UNLESS_CONVERGED {NOCONVERGENCE_SUB(FNAME); *JJ = -1;return;} /* calc. properties in each phases */ yl = Mix.l.y[1]; yv = Mix.v.y[1]; rhol = Mix.l.rho; rhov = Mix.v.rho; hl = bwrmH2(&(Mix), 'L', T, rhol, Mix.l.y[1]); hv = bwrmH2(&(Mix), 'V', T, rhov, Mix.v.y[1]); sl = bwrmS2(&(Mix), 'L', T, rhol, Mix.l.y[1]); sv = bwrmS2(&(Mix), 'V', T, rhov, Mix.v.y[1]); /* calc. mixed-properties */ x = (Z - yl )/(yv - yl); rho = x * rhov + (1.0-x) * rhol; h = x * hv + (1.0-x) * hl; s = x * sv + (1.0-x) * sl; v = 1.0 /rho; Return_1: RECHECK_V(v); *VV = v; RECHECK_H(h); *HH = h; RECHECK_S(s); *SS = s; return; /* --------------------------------------------------------------- II = 1 (P,H,Z -> T,V,S) --------------------------------------------------------------- */ case 2: CHECK_P(P); IF_OUTOFRANGE{*JJ = -2;return;} CHECK_H(H); IF_OUTOFRANGE{*JJ = -2;return;} CHECK_Z(Z); IF_OUTOFRANGE{*JJ = -2;return;} IF_OUTOFRANGE{*JJ = -2;return;} if (OVER_PMAX_M){ OUTOFRANGE_SUB(FNAME,"P", P); *JJ = -2;return; } /* Tb */ Tb = bwrmTyv_Pyl(&Mix, P, Z); CHECK_ERR; /* Td */ if (Mix.v.y[1] < Z) yl = (1+Z)/2.0; else yl = Z / 2.0; Td = bwrmTyl_Pyv2(&Mix, P, Z, Tb+5.0, yl); CHECK_ERR; /* T */ T = bwrmT_Phy(&Mix, P, H, Z); CHECK_ERR; yl = Mix.l.y[1]; yv = Mix.v.y[1]; rhol = Mix.l.rho; rhov = Mix.v.rho; if ( T < Tb ) { rhol = bwrmRHOl(&(Mix.l), P, T); sl = bwrmS2(&Mix, 'L', T, rhol, Z); v = 1.0 /rhol; s = sl; Mix.b = Mix.l; }else if ( T > Td ){ rhov = bwrmRHOv(&(Mix.v), P, T); sv = bwrmS2(&Mix, 'V', T, rhov, Z); v = 1.0 /rhov; s = sv; Mix.b = Mix.v; }else{ if (Z == 0.0 || Z == 1.0){ OUTOFRANGE_SUB(FNAME,"Z", Z); *JJ = -2;return; } sl = bwrmS2(&Mix, 'L', T, rhol, yl); sv = bwrmS2(&Mix, 'V', T, rhov, yv); x = (Z - yl)/(yv - yl); rho = x * rhov + (1.0-x) * rhol; v = 1.0 /rho; s = x * sv + (1.0-x) * sl; } RECHECK_V(v); RECHECK_S(s); RECHECK_T(T); *VV = v; *SS = s; *TT = T; return; /* ---------------------------- II = 3 --------------------------- */ default: ; } OUTOFRANGE_SUB(FNAME, "I", *II); *JJ = -2; return; #undef CHECK_ERR } /********************************************************************* SUBXY ********************************************************************/ subroutine SUBXY( integer JJ, real TT, real PP, real ZL, real ZV, real VL, real VV, real HL, real HV, real SL, real SV ){ double T, P, yv, yl, rhov, rhol, vv, vl, hv, hl, sv, sl, Tb, Td; strcpy(FNAME,"SUBXY"); *JJ = 0; T = *TT; CHECK_T(T); IF_OUTOFRANGE{*JJ = -2;return;} P = *PP; CHECK_P(P); IF_OUTOFRANGE{*JJ = -2;return;} if (OVER_PMAX_M){ OUTOFRANGE_SUB(FNAME,"P[Pa]", P); *JJ = -2;return; } if (T >= Pu[1].Tc && T >= Pu[2].Tc){ OUTOFRANGE_SUB(FNAME,"T[K]", T); *JJ = -2;return; } if (T >= Pu[1].Tc && T >= Pu[2].Tc){ OUTOFRANGE_SUB(FNAME,"T[K]", T); *JJ = -2;return; } if (Pu[1].Tc < Pu[2].Tc) { yl = 0.3; yv = 0.7; }else{ yl = 0.7; yv = 0.3; } bwrmYvYl_PT(&Mix, P, T, yv, yl); IF_BWRM_OUTOFRANGE {OUTOFRANGE_SUB(FNAME,"?", 0.0); *JJ = -2;return;} UNLESS_CONVERGED {NOCONVERGENCE_SUB(FNAME); *JJ = -1;return;} yv = Mix.v.y[1]; yl = Mix.l.y[1]; rhov = Mix.v.rho; rhol = Mix.l.rho; vv = 1.0 / rhov; vl = 1.0 / rhol; hl = bwrmH2(&(Mix), 'L', T, rhol, yl); hv = bwrmH2(&(Mix), 'V', T, rhov, yv); sl = bwrmS2(&(Mix), 'L', T, rhol, yl); sv = bwrmS2(&(Mix), 'V', T, rhov, yv); IF_KG_UNIT{ yv = mol2kgZ(yv); yl = mol2kgZ(yl); vv /= Mix.v.M; vl /= Mix.l.M; hv /= Mix.v.M; hl /= Mix.l.M; sv /= Mix.v.M; sl /= Mix.l.M; /* printf("b:%f, v:%f, l:%f\n", Mix.b.M, Mix.v.M, Mix.l.M); /**/ } /* printf("%f %f %f %f %f %f %f %f %f %f \n", T, P, yl, yv, vl, vv, hv,hl, sl, sv);/**/ *ZL = (float)yl; *ZV = (float)yv; *VL = (float)vl; *VV = (float)vv; *HL = (float)hl; *HV = (float)hv; *SL = (float)sl; *SV = (float)sv; return; } /********************************************************************* SUBTB & SUBTD ********************************************************************/ #define CHECK_ERR \ IF_BWRM_OUTOFRANGE {OUTOFRANGE_SUB(FNAME,"?", 0.0); *JJ = -2;return;}\ UNLESS_CONVERGED {NOCONVERGENCE_SUB(FNAME); *JJ = -1;return;} subroutine SUBTB( integer JJ, real TT, real PP, real ZZ, real VV, real HH, real SS ){ double T, P, Z, V, H, S, rho, v; strcpy(FNAME,"SUBTB"); *JJ = 0; Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{ *JJ = -2;return ;} P = *PP; CHECK_P(P); IF_OUTOFRANGE{ *JJ = -2;return ;} if (OVER_PMAX_M){ OUTOFRANGE_SUB(FNAME,"P", P); *JJ = -2;return; } T = bwrmTyv_Pyl(&Mix, P, Z); CHECK_ERR; rho = bwrmRHOl(&(Mix.l), P, T); /* bubble point -> liquid */ v = 1.0 / rho; H = bwrmH2(&(Mix), 'L', T, rho, Z); S = bwrmS2(&(Mix), 'L', T, rho, Z); Mix.b = Mix.l; RECHECK_T(T); *TT = T; RECHECK_V(v); *VV = v; RECHECK_H(H); *HH = H; RECHECK_S(S); *SS = S; return; } subroutine SUBTD( integer JJ, real TT, real PP, real ZZ, real VV, real HH, real SS ){ double T, P, Z, V, H, S, rho, v; strcpy(FNAME,"SUBTD"); *JJ = 0; Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{ *JJ = -2;return ;} P = *PP; CHECK_P(P); IF_OUTOFRANGE{ *JJ = -2;return ;} if (OVER_PMAX_M){ OUTOFRANGE_SUB(FNAME,"P", P); *JJ = -2;return; } T = bwrmTyl_Pyv(&Mix, P, Z); CHECK_ERR; rho = bwrmRHOv(&(Mix.v), P, T); /* dew point -> vapor */ v = 1/rho; H = bwrmH2(&(Mix), 'V', T, rho, Z); S = bwrmS2(&(Mix), 'V', T, rho, Z); Mix.b = Mix.v; RECHECK_T(T); *TT = T; RECHECK_V(v); *VV = v; RECHECK_H(H); *HH = H; RECHECK_S(S); *SS = S; return; } #undef CHECK_ERR /********************************************************************* TBP & TDP ********************************************************************/ #define CHECK_ERR \ IF_BWRM_OUTOFRANGE {OUTOFRANGE_REAL(FNAME,"P, X or Y)", 0.0)}\ UNLESS_CONVERGED {NOCONVERGENCE_REAL(FNAME);} function_real TBP(real PP, real ZZ) { double P, Z, T; strcpy(FNAME,"TBP"); Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{ return -1.0e-20;} P = *PP; CHECK_P(P); IF_OUTOFRANGE{ return -1.0e-20;} if (OVER_PMAX_M) OUTOFRANGE_REAL(FNAME,"P", P); T = bwrmTyv_Pyl(&Mix, P, Z); CHECK_ERR; RECHECK_T(T); return T; } function_real TDP(real PP, real ZZ) { double P, Z, T; strcpy(FNAME,"TDP"); Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{ return -1.0e-20;} P = *PP; CHECK_P(P); IF_OUTOFRANGE{ return -1.0e-20;} if (OVER_PMAX_M) OUTOFRANGE_REAL(FNAME,"P", P); T = bwrmTyl_Pyv(&Mix, P, Z); CHECK_ERR; RECHECK_T(T); return T; } #undef CHECK_ERR /********************************************************************* SUBPB & SUBPD ********************************************************************/ #define CHECK_ERR \ IF_BWRM_OUTOFRANGE {OUTOFRANGE_SUB(FNAME,"?", 0.0); *JJ = -2;return;}\ UNLESS_CONVERGED {NOCONVERGENCE_SUB(FNAME); *JJ = -1;return;} subroutine SUBPB( integer JJ, real TT, real PP, real ZZ, real VV, real HH, real SS ){ double T, P, Z, V, H, S, rho, v; strcpy(FNAME,"SUBPB"); *JJ = 0; Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{*JJ = -2;return;} T = *TT; CHECK_T(T); IF_OUTOFRANGE{*JJ = -2;return;} P = bwrmPyv_Tyl(&Mix, T, Z); CHECK_ERR; rho = bwrmRHOl(&(Mix.l), P, T); /* bubble point -> liquid */ v = 1.0 / rho; H = bwrmH2(&(Mix), 'L', T, rho, Z); S = bwrmS2(&(Mix), 'L', T, rho, Z); Mix.b = Mix.l; RECHECK_P(P); *PP = P; RECHECK_V(v); *VV = v; RECHECK_H(H); *HH = H; RECHECK_S(S); *SS = S; return; } subroutine SUBPD( integer JJ, real TT, real PP, real ZZ, real VV, real HH, real SS ){ double T, P, Z, V, H, S, rho, v; strcpy(FNAME,"SUBPD"); *JJ = 0; Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{ *JJ = -2;return ;} T = *TT; CHECK_T(T); IF_OUTOFRANGE{ *JJ = -2;return ;} P = bwrmPyl_Tyv(&Mix, T, Z); CHECK_ERR; rho = bwrmRHOv(&(Mix.v), P, T); /* dew point -> vapor */ v = 1/rho; H = bwrmH2(&(Mix), 'V', T, rho, Z); S = bwrmS2(&(Mix), 'V', T, rho, Z); Mix.b = Mix.v; RECHECK_P(P); *PP = P; RECHECK_V(v); *VV = v; RECHECK_H(H); *HH = H; RECHECK_S(S); *SS = S; return; } #undef CHECK_ERR /********************************************************************* PBT & PDT ********************************************************************/ #define CHECK_ERR \ IF_BWRM_OUTOFRANGE {OUTOFRANGE_REAL(FNAME,"T, X or Y)", 0.0)}\ UNLESS_CONVERGED {NOCONVERGENCE_REAL(FNAME);} function_real PBT(real TT, real ZZ) { double P, Z, T; strcpy(FNAME,"PBT"); T = *TT; CHECK_T(T); IF_OUTOFRANGE{ return -1.0e-20;} Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{ return -1.0e-20;} P = bwrmPyv_Tyl(&Mix, T, Z); CHECK_ERR; RECHECK_P(P); return P; } function_real PDT(real TT, real ZZ) { double P, Z, T; strcpy(FNAME,"PDT"); T = *TT; CHECK_T(T); IF_OUTOFRANGE{ return -1.0e-20;} Z = *ZZ; CHECK_Z(Z); IF_OUTOFRANGE{ return -1.0e-20;} P = bwrmPyl_Tyv(&Mix, T, Z); CHECK_ERR; RECHECK_P(P); return P; } #undef CHECK_ERR /******************************************************************* FUNCTIONs/SUBROUTINEs for PURE COMPONENT ********************************************************************/ /********************************************************************* SUBPUR ********************************************************************/ subroutine SUBPUR( integer II, integer JJ, real TT, real PP, real VV, real HH, real SS) { int i; double T, P, Ps, rho, v, h, s; strcpy(FNAME,"SUBPUR"); *JJ = 0; i = *II; if ( (i != 1) && (i != 2) ){ OUTOFRANGE_SUB(FNAME, "I", i); *JJ = -2; return; } T = *TT; CHECK_T(T); IF_OUTOFRANGE{ *JJ = -2;return ;} P = *PP; CHECK_P(P); IF_OUTOFRANGE{ *JJ = -2;return ;} if (OVER_PMAX_P(i)){ OUTOFRANGE_SUB(FNAME,"P", P); *JJ = -2;return; } Ps = bwrmPs(&(Pu[i]), T); IF_BWRM_OUTOFRANGE_P(i) {OUTOFRANGE_SUB(FNAME,"T[K]", T); *JJ = -2;return;} UNLESS_CONVERGED_P(i) {NOCONVERGENCE_SUB(FNAME); *JJ = -1;return;} if ( P < Ps ){ rho = bwrmRHOv(&(Pu[i]), P, T); }else{ rho = bwrmRHOl(&(Pu[i]), P, T); } v = 1.0 / rho; h = bwrmH(&(Pu[i]), T, rho); s = bwrmS(&(Pu[i]), T, rho); IF_KG_UNIT{ v /= Pu[i].M; h /= Pu[i].M; s /= Pu[i].M; } *VV = v; *HH = h; *SS = s; } /********************************************************************* SUBTSP, TSPM ********************************************************************/ subroutine SUBTSP( integer II, integer JJ, real TT, real PP, real VL, real VV, real HL, real HV, real SL, real SV ){ int i; double Ts, P, vv, vl, hv, hl, sv, sl; double Mw, rhov, rhol; strcpy(FNAME,"SUBTSP"); *JJ = 0; i = *II; if ( (i != 1) && (i != 2) ){ OUTOFRANGE_SUB(FNAME, "I", i); *JJ = -2; return; } P = *PP; CHECK_P(P); IF_OUTOFRANGE{ *JJ = -2;return ;} if (OVER_PMAX_P(i)) {OUTOFRANGE_SUB(FNAME,"P", P); *JJ = -2;return;} Ts = bwrmTs(&(Pu[i]), P); IF_BWRM_OUTOFRANGE_P(i) {OUTOFRANGE_SUB(FNAME,"P", P); *JJ = -2;return;} UNLESS_CONVERGED_P(i) {NOCONVERGENCE_SUB(FNAME); *JJ = -1;return;} rhov = Pu[i].rhov; rhol = Pu[i].rhol; vv = 1.0 / rhov; vl = 1.0 / rhol; vv = 1.0 / rhov; vl = 1.0 / rhol; hv = bwrmH(&(Pu[i]), Ts, rhov); hl = bwrmH(&(Pu[i]), Ts, rhol); sv = bwrmS(&(Pu[i]), Ts, rhov); sl = bwrmS(&(Pu[i]), Ts, rhol); IF_KG_UNIT{ Mw = Pu[i].M; vv /= Mw; vl /= Mw; hv /= Mw; hl /= Mw; sv /= Mw; sl /= Mw; } RECHECK_T(Ts); *TT = Ts; *VV = vv; *VL = vl; *HV = hv; *HL = hl; *SV = sv; *SL = sl; } function_real TSPM(integer II, real PP) { int i; double P, T; strcpy(FNAME,"TSPM"); i = *II; if ( (i != 1) && (i != 2) ){ OUTOFRANGE_SUB(FNAME, "I", i); return -1.0e+20; } P = *PP; CHECK_P(P); IF_OUTOFRANGE{ return -1.0e+20;} if (OVER_PMAX_P(i)) OUTOFRANGE_REAL(FNAME,"P", P); Pu[i].tmp.errorcode = BWRMERR_NOPROBLEM; T = bwrmTs(&(Pu[i]), P); IF_BWRM_OUTOFRANGE_P(i) {OUTOFRANGE_REAL(FNAME,"P", P);} UNLESS_CONVERGED_P(i) {NOCONVERGENCE_REAL(FNAME);} RECHECK_T(T); return T; } /********************************************************************* SUBPST, PSTM ********************************************************************/ subroutine SUBPST( integer II, integer JJ, real TT, real PP, real VL, real VV, real HL, real HV, real SL, real SV ){ int i; double Ps, T, vv, vl, hv, hl, sv, sl; strcpy(FNAME,"SUBPST"); *JJ = 0; i = *II; if ( (i != 1) && (i != 2) ){ OUTOFRANGE_SUB(FNAME, "I", i); *JJ = -2; return; } T = *TT; CHECK_T(T); IF_OUTOFRANGE{ *JJ = -2;return ;} Ps = bwrmPs(&(Pu[i]), T); IF_BWRM_OUTOFRANGE_P(i) {OUTOFRANGE_SUB(FNAME,"T", T); *JJ = -2;return;} UNLESS_CONVERGED_P(i) {NOCONVERGENCE_SUB(FNAME); *JJ = -1;return;} vv = 1.0 / Pu[i].rhov; vl = 1.0 / Pu[i].rhol; hv = bwrmH(&(Pu[i]), T, Pu[i].rhov); hl = bwrmH(&(Pu[i]), T, Pu[i].rhol); sv = bwrmS(&(Pu[i]), T, Pu[i].rhov); sl = bwrmS(&(Pu[i]), T, Pu[i].rhol); IF_KG_UNIT{ vv /= Pu[i].M; vl /= Pu[i].M; hv /= Pu[i].M; hl /= Pu[i].M; sv /= Pu[i].M; sl /= Pu[i].M; } RECHECK_P(Ps); *PP = Ps; *VV = vv; *VL = vl; *HV = hv; *HL = hl; *SV = sv; *SL = sl; } function_real PSTM(integer II, real TT) { int i; double P, T; strcpy(FNAME,"TSPM"); i = *II; if ( (i != 1) && (i != 2) ){ OUTOFRANGE_SUB(FNAME, "I", i); return -1.0e+20; } T = *TT; CHECK_T(T); IF_OUTOFRANGE{ return -1.0e+20;} P = bwrmPs(&(Pu[i]), T); IF_BWRM_OUTOFRANGE_P(i) {OUTOFRANGE_REAL(FNAME,"T", T);} UNLESS_CONVERGED_P(i) {NOCONVERGENCE_REAL(FNAME);} RECHECK_P(P); return P; }