/* bwrm.h : headder for subroutin on modified BWR equation . by MOMOKI,Satoru Jan. 1996 */ #ifndef BWRM_INCLUDED # define T0bwrm 273.15 #ifndef Rgas # define Rgas 8.31451e+003 /* Universal gas constant [1/(kg-mole.K)] */ #endif /************************************************************************ BWRM Structure for properties and BWR coefficients ************************************************************************/ #define NMAX_BWRM 3 /* ----------------------------------------------------------------------- coefficients_bwrm : BWRM COEFFICIENTS ----------------------------------------------------------------------- */ struct coefficients_bwrm{ double A0, B0, C0, D0, E0; /* coefficients */ double a, b, c, d, e, f, g, h; /* coefficients */ double alpha, gamma; /* coefficients */ /* for binary-mixtures */ double B0_12; double A0_12; double C0_12; double D0_12; double E0_12; };typedef struct coefficients_bwrm BWRM_COEF ; /* ----------------------------------------------------------------------- bwrm_buffer : temporary buffer to store calculated values such as T^2, T^4, and so on ----------------------------------------------------------------------- */ struct bwrm_buffer{ double yprev; double Tprev ; double rhoprev; double T, T2, T3, T4, T5, T6, T7, T8, T9, T16, T17, T18, T22, T23, T24; double rho, rho2, rho3, rho4, rho5, rho6, rho7; double rhoRT , B0RT, C0_T2, D0_T3, E0_T4, bRT, d_T, e_T4, f_T23, alphaa, alphad_T, alphae_T4, alphaf_T23, c_T2, g_T8, h_T17, gammarho2, EXP1; double B0_12RT, C0_12_T2, D0_12_T3, E0_12_T4;/* used in the case of binary mixtures */ /* messages from functions */ long int n_iteration; int errorcode; #define BWRMERR_NOPROBLEM 0 #define BWRMERR_INVALID_ARGUMENT 1 #define BWRMERR_ARGUMENT_OUTOFRANGE 2 #define BWRMERR_CANT_CONVERGE 10 #define BWRMERR_CANT_RESOLVE_MATRIX 11 #define BWRMERR_FUNNY_RESULT 12 #define BWRMERR_INACCURATE_RESULT 13 };typedef struct bwrm_buffer BWRM_BUF; /* ----------------------------------------------------------------------- BWRM (struct bwrm) : structure for properties and BWR coefficients. ----------------------------------------------------------------------- */ struct bwrm { /* ---------- CHARACTERISTIC -------------- */ int n; /* number of components */ BWRM_COEF cf; /* BWR coefficients */ double omega, (*omega_f)(double t); BWRM_BUF tmp; /* temporary buffer to store calculated values such as T^2, T^4, and so on */ double M; /* Molecular weight [kg/kg-mole](=[g/g-mol]) */ double Pc, Tc, RHOc; /* properties at critical point. P[Pa], T[K], rho[kg-mole/m3] */ double (*cp_id)(double t); /* molar-specific heat [J/kg-moleK] at ideal gas condition, cp = bwrmCp(t) + cp_id(t) */ double (*h_id)(double t), h_adj; /* molar-enthalpy [J/kg-mole] at ideal gas condition. h = bwrmdH(t) + h_id(t) + h_adj [J/kg-mol] */ double (*s_id)(double t), s_adj; /* molar-entropy at [J/kg-moleK] ideal gas condition */ /* ---------- PROPERTIES -------------- */ double P; /* pressure [Pa] */ double T; /* temperatutr [K] */ double rho; /* molar density */ double cp; /* molar specific heat [J/kg-moleK] */ double h; /* molar enthalpy [J/kg-mole] */ double s; /* molar entropy [J/kg-mole] */ double x; /* molar fraction of vapor[mol/mol] */ /* -- for pure substances -- */ double rhov, rhol; /* molar density of vapor and liauid */ double fv, fl; /* fugacity of vapor and liauid */ double muv, mul; /* chemical potential vapor and liauid */ /* -- for mixtures -- */ double f[NMAX_BWRM+1]; /* fugacity of componet i */ double y[NMAX_BWRM+1]; /* mole fraction of componet i */ double mu[NMAX_BWRM+1]; /* chemical potential of componet i */ };typedef struct bwrm BWRM; /* ----------------------------------------------------------------------- BWRM2 (struct bwrm2) : structure for properties and BWR coefficients. ----------------------------------------------------------------------- */ struct bwrm2{ BWRM b, v, l; /* properties and coefficients of bulk, vapor, liquid, respectively */ BWRM c[NMAX_BWRM+1]; /* properties and coefficients of component 1 and 2, respectively */ /* 1 : more volatile component */ /* 2 : less volatile component */ double (*mij)(double y, double t); double (*h_id)(double y, double t); /* yb, bulk molar fraction, -> b.y[1] Tc, RHOc, M. -> b.Tc, b.RHOc,..... cp0(t), molar-specific heat at ideal gas condition[J/kg-moleK] -> b.cp0_f(t) h0(t), molar-enthalpy at ideal gas condition[J/kg-mole] -> b.h0_f(t) s0(p,t), molar-enthalpy at ideal gas condition[J/kg-mole] -> b.s0_f(p,t) */ };typedef struct bwrm2 BWRM2; /****************************************************************** * * * FUCTIONS and MACROS * * * ******************************************************************/ /****************************************************************** bwrm1.c functions concerning to BWR coefficients ******************************************************************/ int InitBWRM(/* initializes BWRM structure for pure substance, 'c'. */ BWRM *c, double tc, /* critical temperature [K] */ double rhoc, /* molar density [kg-mol/m3] */ double M, /* molecular weight [kg/kg-mole] */ double (*w_f)(double t), /* function to calc. omega from T. If omwga is constant, w_f = bwrmOMEGA_CONST. */ double w, /* omega when w_f = bwrmOMEGA_CONST. */ double t, /* temperature used for w_f. */ double (*cp_id)(double t), /* function to calc. cp at ideal gas [J/kg-mol K]*/ double (*h_id)(double t), /* function to calc. h at ideal gas [J/kg-mol] */ double (*s_id)(double t) /* function to calc. s at ideal gas [J/kg-mol K]*/ ); double bwrmOMEGA_CONST(double t); #define RecalBWRM(C, T, RHO) \ {int FLAG=0; \ if (((T) != (*C).tmp.Tprev) && ((*C).omega_f != bwrmOMEGA_CONST)){ \ CalBWRM(C, (*C).Tc, (*C).RHOc, (*C).M, (*C).omega_f(T)); \ FLAG = 1;}\ if ((T) != (*C).tmp.Tprev) {bwrmCalTmpT(C, T); FLAG=1;}\ if ((RHO) != (*C).tmp.rhoprev) {bwrmCalTmpRHO(C, RHO); FLAG=1;}\ if ((FLAG) == 1) {bwrmCalTmpOther(C, T, RHO);}\ } #define RecalBWRMcf(C, T) \ {if (((T) != (*C).tmp.Tprev) && ((*C).omega_f != bwrmOMEGA_CONST)){ \ CalBWRM(C, (*C).Tc, (*C).RHOc, (*C).M, (*C).omega_f(T)); \ }} /* ------------------------------------------------------------- */ int InitBWRM2(/* initializes BWR coefficients for binary mixtures, 'mx'. 'pu1' and 'pu2' are required to be already caluculated using CalBWRM(). */ BWRM2 *mx, BWRM *pu1, BWRM *pu2, double y1, /* molar fraction of component 1, [mole/mole] */ double (*mij)(double y, double t) ); #define RecalBWRM2(M, MODE, T, RHO, Y) \ {BWRM *C; int FLAG=0;\ C = &((*M).b); /* bulk */ \ if (MODE == 'l' || MODE == 'L') C = &((*M).l); /* liquid */ \ if (MODE == 'v' || MODE == 'V') C = &((*M).v); /* vapor */ \ RecalBWRMcf(&((*M).c[1]), T); RecalBWRMcf(&((*M).c[2]), T); \ if (Y != (*C).tmp.yprev) {CalBWRM2(M, MODE, Y); FLAG=1;} \ if (T != (*C).tmp.Tprev) {bwrmCalTmpT(C, T); FLAG=1;} \ if (RHO != (*C).tmp.rhoprev) {bwrmCalTmpRHO(C, RHO); FLAG=1;} \ if (FLAG == 1) {bwrmCalTmpOther(C, T, RHO);} \ } /* ------------------------------------------------------------- */ int CalBWRM(/* calculates BWR coefficients for pure substance, 'coef'. */ BWRM *c, double tc, /* critical temperature [K] */ double rhoc, /* molar density [kg-mol/m3] */ double M, /* molecular weight [kg/kg-mole] */ double w /* omega */ ); int CalBWRM2(/* calculates BWR coefficients'c.b(vb)' with consideration of composition. Caluculation of 'c.c[1]' and 'c.c[2]' by CalBWRM() and definetion of c.mij are required before you call this function. */ BWRM2 *m, char mode, /* 'b' or 'B' : calulates (*c).b) (bulk), 'v' or 'V' : calulates (*c).b) (vapor), 'l' or 'L' : calulates (*c).b) (liquid) */ double y1 ); /* ------------------------------------------------------------- */ int bwrmCalTmp(/* calculates coef.tmp */ BWRM *coef, double t, /* [K] */ double rho ); /* [kg-mole/m3] */ #define CalBWRMTmp(c, t, rho) \ {int mode=0;\ if (t != (*c).tmp.Tprev){bwrmCalTmpT(c, t);mode=1;}\ if (rho != (*c).tmp.rhoprev){bwrmCalTmpRHO(c, rho);mode=1;}\ if (mode == 1){bwrmCalTmpOther(c, t, rho);}\ } /* calculates coef.tmp if T,rho or y were updated */ int bwrmCalTmpT(/* calculates T^2, T^4, and so on, and substitute them into 'coef.tmp'. */ BWRM *coef, double t ); #define CalBWRMTmpT(c, t) \ {if (t != (*c).tmp.Tprev){\ bwrmCalTmpT(c, t);bwrmCalTmpOther(c, t, rho);}} /* calculates T^2, T^4, .e.t.c if T was updated */ int bwrmCalTmpRHO(/* calculates rho^2, rho^4, and so on, and substitute them into 'coef.tmp'. */ BWRM *coef, double rho ); #define CalBWRMTmpRHO(c, rho) \ {if (rho != (*c).tmp.rhoprev){\ bwrmCalTmpRHO(c, rho);bwrmCalTmpOther(c, t, rho);}} /* calculates rho^2, rho^4, .e.t.c if rho was updated */ int bwrmCalTmpOther(/* calculates rho*R*T, B0*R*, and so on, and substitute them into 'coef.tmp'. */ BWRM *coef, double t, /* [K] */ double rho ); /* [kg-mole/m3] */ /* ------------------------------------------------------------- */ int BWRMcpy(BWRM *a, BWRM *b ); /* coef'a' <- coef'b' */ int BWRM2cpy(BWRM2 *a, BWRM2 *b ); /* coef'a' <- coef'b' */ int BWRMCTMPcpy(BWRM *a, BWRM *b ); /* copy coef.tmp of 'b' to 'a' */ /****************************************************************** bwrm2.c (pressure,P) ******************************************************************/ /* ------------------------------------------------------------- */ double bwrmP(/* calculates pressure 'P' in [Pa] using BWR eq. # include */ BWRM *coef, /* structure of BWR coefficients */ double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); double bwrmdPdrho(/* calculates dP/dRHO. */ BWRM *coef, double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); double bwrmdPdT(/* calculates dP/dT. */ BWRM *coef, double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); /****************************************************************** bwrm3.c (density, rho, and temperature) ******************************************************************/ /* ------------------------------------------------------------- */ double bwrmRHOv(/* calculates molar density 'rho' in [kg-mole/m3] for *VAPOR* using BWR eq. and Newton method */ BWRM *c, double p, /* pressure [Pa] */ double t /* temperature [K] */ ); double bwrmRHOl(/* calculates molar density 'rho' in [kg-mole/m3] for *LIQUID* using BWR eq. and Newton method */ BWRM *c, double p, /* pressure [Pa] */ double t /* temperature [K] */ ); double bwrmRHO(/* calculates molar density 'rho' in [kg-mole/m3] using BWR eq. and Newton method */ BWRM *coef, double p, /* pressure [Pa] */ double t, /* temperature [K] */ double rhostart ); /* ------------------------------------------------------------- */ double bwrmT(/* calculates temperature in [K] from P and rho */ BWRM *coef, double p, /* pressure [Pa] */ double rho, /* molaxr density [kg-mole/m3] */ double tstart /* initial value of temperature for iteration[K] */ ); /****************************************************************** bwrm4.c (molar enthalpy and molar entropy) ******************************************************************/ /* ------------------------------------------------------------- */ int bwrmSetHid200SI(/* Sets (*m).h_id, and determines (*m).h_adj so that the result of bwrmH(t) for saturated liquid at 273.15[K] is 200,000 * (*m).M [J/kg-mol] */ BWRM *m, double (*f)(double t) /* function that calculates 'h'[J/kg-mole] at ideal gas condition ~~~~~~~~~~ */ ); double bwrmH(/* calculates molar enthalpy 'h'[J/kg-mol]. You have to set (*c).h_id(t) and (*c).h_adj with bwrmSetHid200SI() before you use this function. */ BWRM *coef, /* structure of BWR coefficients */ double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); double bwrmdH(/* calculates the deflection from ideal gas condition on molar enthalpy, 'h - h_id'[J/kg-mol] */ BWRM *coef, /* structure of BWR coefficients */ double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); /* ------------------------------------------------------------- */ double bwrmdS(/* calculates molar entropy 's-s0'[J/kg-mol] using BWR eq. # include */ BWRM *coef, /* structure of BWR coefficients */ double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); double bwrmS(/* calculates molar entropy 's-s0'[J/kg-mol] using BWR eq. You have to set (*c).s_id(t) and (*c).s_adj with bwrmSetSid1SI() before you use this function. */ BWRM *coef, /* structure of BWR coefficients */ double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); /****************************************************************** bwrm5.c (fugacity and chemical potential) ******************************************************************/ /* ------------------------------------------------------------- */ double bwrmRTlnF(/* calculates RTln(f) of pure substances, where f is fugacity. */ BWRM *coef, double t, /* temperature [K] */ double rho /* molar density [kg-mole/m3] */ ); double bwrmF(/* calculates f of pure substances */ BWRM *coef, double t, /* temperature [K] */ double rho /* molar density [kg-mole/m3] */ ); /* ------------------------------------------------------------- */ double bwrmF2(/* calculates f (RTln(f)*y) of binary mixtures */ BWRM2 *c, /* BWR coefficients of mixtures */ char mode1,/* '1' or 1 : component 1 '2' or 2 : component 2 */ char mode2, /* 'l' or 'L' : calulates (*c).b) (liquid) 'v' or 'V' : calulates (*c).b) (vapor), else : calulates (*c).b) (bulk) */ double t, /* temperature [K] */ double rho, /* molar density [kg-mole/m3] */ double y /* molar fraction [mol/mol] */ ); /****************************************************************** bwrm6.c (molar specific heat ) ******************************************************************/ /* ------------------------------------------------------------- */ double bwrmdCp(/* calculates the deflection from ideal gas condition on molar specific heat with constant pressure. [J/kg-mol K] */ BWRM *coef, /* structure of BWR coefficients */ double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); double bwrmCp(/* calculates molar specific heat with constant pressure, 'cp' [J/kg-mol K]. ** This function require coef.cp0() ** */ BWRM *coef, /* structure of BWR coefficients */ double t, /* temperature [K] */ double rho /* molar density [kg-mol/m3] */ ); /****************************************************************** bwrm7.c (saturation temperature for pure substances ) ******************************************************************/ /* ------------------------------------------------------------- */ double bwrmPs(/* calculates saturation pressure in [Pa] from T[K] */ BWRM *c, double T /* temperature [K] */ ); double bwrmTs(/* calculates saturation temperature in [K] from P */ BWRM *c, double p /* pressure[Pa] */ ); /****************************************************************** bwrm8.c (temperature from enthalpy) ******************************************************************/ double bwrmT_Ph(/* Determines temperature from pressure,P[Pa] and molar enthalpy 'h'[J/kg-mol]. You have to set (*c).h_id(t) and (*c).h_adj with bwrmSetHid200SI() before you use this function. */ BWRM *c, /* structure of BWR coefficients */ char mode,/* 'l' or 'L' : calulates (*c).b) (liquid) 'v' or 'V' : calulates (*c).b) (vapor) */ double p, /* pressure [Pa] */ double h /* molar enthalpy [J/kg-mol] */ ); /****************************************************************** BINARY MIXTURES ******************************************************************/ #define F0_BWRM2(COEF) ( ( ((*COEF).v.f[1]) - ((*COEF).l.f[1]) ) ) /* f1v -f1l */ #define F1_BWRM2(COEF) ( ( ((*COEF).v.f[2]) - ((*COEF).l.f[2]) ) ) /* f2v -f2l */ int CALC_F0F1_BWRM2(BWRM2 *COEF, double P, double T, double YV1, double YL1); /* in bwrm11.c */ /****************************************************************** bwrm11.c (yv and yl at saturated condition for binary mixtures) ******************************************************************/ double bwrmYvYl_PT(/* calculates yl[mol/mol] & yv[mol/mol] from P[Pa] and T[K] of bibary bibary mixtures at saturated condition. Return value is yv[mol/mol]. And yl[mol/mol], yv, rhol[kg-mol/m3] and rhov are set on (*c).l.y[1], (*c).v.y[1], (*c).prop_l.rho, and (*c).v.rho. */ BWRM2 *m, double p, /* pressure [Pa] */ double t, /* temperature [K] */ double yvstart, /* initial value on y_v [mol/mol] */ double ylstart /* initial value on y_l [mol/mol] */ ); /****************************************************************** bwrm12.c (saturation temperature e.t.c for binary mixtures) ******************************************************************/ 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. This can be also used for computing T and yv at bubble point from P and yb. */ BWRM2 *m, double p, /* pressure [Pa] */ double yl /* molar fraction in liquid [mol/mol] */ ); double bwrmTyv_Pyl2(BWRM2 *m, double p, double yl, double t0, /* initial value of T */ double yv0/* initial value of yv */); double bwrmTyl_Pyv(/* calculates T[K] & yl[mol/mol] from P[Pa] and yv[mol/mol] of bibary bibary 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. This can be also used for computing T and yl at dew point from P and yb. */ BWRM2 *m, double p, /* pressure [Pa] */ double yv /* molar fraction in liquid [mol/mol] */ ); double bwrmTyl_Pyv2(BWRM2 *m, double p, double yv, double t0, /* initial value of T */ double yl0/* initial value of yv */); /****************************************************************** bwrm13.c (saturation pressure e.t.c for binary mixtures) ******************************************************************/ double bwrmPyv_Tyl(/* calculates P[Pa] & yv[mol/mol] from T[K] and yl[mol/mol] of bibary bibary mixtures at saturated condition. Return value is P[Pa]. 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. This can be also used for computing T and yv at bubble point from P and yb. */ BWRM2 *m, double t, /* temperature [K] */ double yl /* molar fraction in liquid [mol/mol] */ ); double bwrmPyl_Tyv(/* calculates P[Pa] & yl[mol/mol] from T[K] and yv[mol/mol] of bibary bibary mixtures at saturated condition. Return value is P[Pa]. 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. This can be also used for computing T and yl at dew point from P and yb. */ BWRM2 *m, double t, /* temperature [K] */ double yv /* molar fraction in liquid [mol/mol] */ ); /****************************************************************** bwrm14.c (molar enthalpy of binary mixtures) ******************************************************************/ double bwrmH2(/* calculates molar enthalpy 'h'[J/kg-mol] of binary mixtures. The function 'bwrmSetHid200SI' must be already done for BWRM strucure of each component,c[1] and c[2], respectively. */ BWRM2 *m, char mode, /* 'l' or 'L': BWRM structure,l for liquid phase is used. 'v' or 'V': BWRM structure,v for liquid phase is used. else : BWRM structure,b for bulk fluid is used. */ double t, /* [K] */ double rho, /* [kg-mol/m3] */ double y1 /* molar fraction of component 1 */ ); double bwrmT_Ph2(/* Determines temperature of binary mixtures from molar enthalpy, 'h'[J/kg-mol]. You have to set (*c).h_id(t) and (*c).h_adj with bwrmSetHid200SI() before you use this function. */ BWRM2 *m, /* structure of BWR coefficients */ char mode,/* 'l' or 'L' : calulates (*c).b) (liquid) 'v' or 'V' : calulates (*c).b) (vapor) */ double p, /* pressure [Pa] */ double h, /* molar enthalpy [J/kg-mol] */ double y /* molar fraction [mol/mol] */ ); double bwrmT_Phy(/* Determines temperature of binary mixtures from pressure,P[Pa] and molar enthalpy 'h'[J/kg-mol]. And the quality is set on (*m).b.x. You have to set (*m)c[12].h_id(t) and (*m)c[12].h_adj with bwrmSetHid200SI() before you use this function. */ BWRM2 *m, /* structure of BWR coefficients */ double p, /* pressure [Pa] */ double h, /* molar enthalpy [J/kg-mol] */ double yb /* bulk molar fraction [mol/mol] */ ); /****************************************************************** bwrm15.c (molar entropy of binary mixtures) ******************************************************************/ double bwrmS2(/* calculates molar entropy 's'[J/kg-molK] of binary mixtures. The function 'bwrmSetSid1SI' must be already done for BWRM strucure of each component,c[1] and c[2], respectively. */ BWRM2 *m, char mode, /* 'l' or 'L': BWRM structure,l for liquid phase is used. 'v' or 'V': BWRM structure,v for liquid phase is used. else : BWRM structure,b for bulk fluid is used. */ double t, /* [K] */ double rho, /* [kg-mol/m3] */ double y1 /* molar fraction of component 1 */ ); /****************************************************************** MISC. -- bwrm99.c ******************************************************************/ #include int bwrm2ERROR(FILE *fp, char *s, double P, double T, double yv, double yl); int bwrmSolveA22x_b(double A[2][2], double B[2], double X[2]); /* Instant subroutine to solve A(2,2)*X(2)=B(2) */ #define BWRM_INCLUDED #endif