CoreFlows/Models/src/fluid.h

   1 //////////////////////////////////////////////////////////////////////////////\r
   2 //\r
   3 // File:        fluid.h\r
   4 // Directory:   \r
   5 // Version:     \r
   6 //\r
   7 //////////////////////////////////////////////////////////////////////////////\r
   8 \r
   9 #ifndef fluid_included\r
  10 #define fluid_included\r
  11 \r
  12 #include <math.h>\r
  13 #include <Noms.h>\r
  14 #include <SFichier.h>\r
  15 #include <Process.h>\r
  16 \r
  17 class fluid\r
  18 {\r
  19 public :\r
  20     fluid() { fluide = 0; };\r
  21     void set_fluid(Nom &nom)\r
  22     {\r
  23         Noms liste(2); liste(0) = "Na"; liste(1) = "LBe"; \r
  24         fluide = liste.rang(nom);\r
  25         if (fluide == -1)\r
  26         {\r
  27             if (!Process::me()) Cerr << "fluide " << nom << " non reconnu parmi : " << liste << finl;\r
  28             Process::exit();\r
  29         }\r
  30     }\r
  31 \r
  32     //proprietes physiques\r
  33     /* Lois physiques du sodium issues de fits sur DEN/DANS/DM2S/STMF/LMES/RT/12-018/A */\r
  34     int fluide;\r
  35     #define Tk ((T) + 273.15)\r
  36     #define Ksil (fluide ? 3.022e-11 : 1e-9)\r
  37 \r
  38     /* inverse de la densite du liquide (hors pression) */\r
  39     inline double IRhoL(double T)\r
  40     {\r
  41       return fluide ? (9.03e-5 + Tk * (1.003e-8 + Tk * 2.01e-12))\r
  42                     : (9.829728e-4 + Tk * (2.641186e-7 + Tk * (-3.340743e-11 + Tk * 6.680973e-14)));\r
  43     }\r
  44 \r
  45     /* sa derivee */\r
  46     inline double DTIRhoL(double T)\r
  47     {\r
  48       return fluide ? (1.003e-8 + 2.01e-12 * 2)\r
  49                     : (2.641186e-7 + Tk * (-3.340743e-11 * 2 + Tk * 6.680973e-14 * 3));\r
  50     }\r
  51 \r
  52     /* densite du liquide */\r
  53     inline double RhoL(double T, double P)\r
  54     {\r
  55       return exp(Ksil * (P - 1e5)) / IRhoL(T);\r
  56     }\r
  57 \r
  58     /* enthalpie du liquide */\r
  59     inline double HL(double T, double P)\r
  60     {\r
  61       return (fluide ? 97.98e3 + Tk * (159 + Tk * (-2.72e-2 / 2 + Tk * 7.12e-6 / 3))\r
  62                      : 2992600 / Tk - 365487.2 + Tk * (1658.2 + Tk * (-.42395 + Tk * 1.4847e-4)))\r
  63           + (IRhoL(T) - Tk * DTIRhoL(T)) * (1 - exp(Ksil*(1e5 - P))) / Ksil;\r
  64     }\r
  65     \r
  66     /* derivees */\r
  67     inline double DTHL(double T, double P)\r
  68     {\r
  69       return (fluide ? 159 + Tk * (-2.72e-2 + Tk * 7.12e-6)\r
  70                      : -2992600 / (Tk*Tk) + 1658.2 + Tk * (-.42395 * 2 + Tk * 1.4847e-4 * 3))\r
  71              - Tk * (fluide ? 0 : -3.340743e-11 * 2 + Tk * 6.680973e-14 * 6) * (1 - exp(Ksil*(1e5 - P))) / Ksil;\r
  72     }\r
  73     \r
  74     /* inverses par methode de Newton */\r
  75     inline double TLh(double h, double T0, double P)\r
  76     {\r
  77       double T = T0, sec; int i;\r
  78       for (i = 0; i < 100 && dabs( sec = HL(T, P) - h ) > 1e-8; i++)\r
  79         T -= sec / DTHL(T, P);\r
  80       return i < 100 ? T : -1e10;\r
  81     }\r
  82     /* derivees */\r
  83     inline double DTRhoL(double T, double P)\r
  84     {\r
  85       return - exp(Ksil * (P - 1e5)) * DTIRhoL(T) / (IRhoL(T) * IRhoL(T));\r
  86     }\r
  87     inline double DPRhoL(double T, double P)\r
  88     {\r
  89         return Ksil * RhoL(T, P);\r
  90     }\r
  91     \r
  92     /* derivees */\r
  93     inline double DPHL(double T, double P)\r
  94     {\r
  95         return (IRhoL(T) - Tk * DTIRhoL(T)) * exp(Ksil * (1e5 - P));\r
  96     }\r
  97     \r
  98     /* energie volumique du liquide */\r
  99     inline double EL(double T, double P)\r
 100     {\r
 101       return HL(T, P) * RhoL(T, P);\r
 102     }\r
 103     /* derivees */\r
 104     inline double DTEL(double T, double P)\r
 105     {\r
 106       return DTHL(T, P) * RhoL(T, P) + HL(T, P) * DTRhoL(T, P);\r
 107     }\r
 108     inline double DPEL(double T, double P)\r
 109     {\r
 110       return DPHL(T, P) * RhoL(T, P) + HL(T, P) * DPRhoL(T, P);\r
 111     }\r
 112     \r
 113     inline double TLe(double e, double T0, double P)\r
 114     {\r
 115       double T = T0, sec; int i;\r
 116       for (i = 0; i < 100 && dabs( sec = EL(T, P) - e ) > 1e-3; i++)\r
 117         T -= sec / DTEL(T, P);\r
 118       return i < 100 ? T : -1e10;\r
 119     }\r
 120     \r
 121     /* viscosite du liquide*/\r
 122     inline double MuL( double T )\r
 123     {\r
 124       return fluide ? 4.94e-4 * exp(754.1 / Tk)\r
 125                     : exp( -6.4406 - 0.3958 * log(Tk) + 556.835 / Tk );\r
 126     }\r
 127     /* conductivite du liquide */\r
 128     inline double LambdaL( double T )\r
 129     {\r
 130       return fluide ? 3.61 + Tk * (1.517e-2 - Tk * 1.741e-6)\r
 131                     : max(124.67 + Tk * (-.11381 + Tk * (5.5226e-5 - Tk * 1.1848e-8)), 0.045);\r
 132     }\r
 133     /* tension superficielle */\r
 134     inline double SigmaL(double T)\r
 135     {\r
 136         return 0.2405 * pow(1 - Tk / 2503.7, 1.126);\r
 137     }\r
 138     /* Nusselt du liquide -> Skupinski */\r
 139     inline double NuL(double T, double P, double w, double Dh)\r
 140     {\r
 141       double Pe = RhoL(T, P) * dabs(w) * Dh * DTHL(T, P) / LambdaL(T);\r
 142       return 4.82 + 0.0185 * pow(Pe, 0.827);\r
 143     }\r
 144     /* temperature de saturation */\r
 145     inline double Tsat( double P )\r
 146     {\r
 147       double A = 7.8130, B = 11209, C = 5.249e5;\r
 148       return 2 * C / ( -B + sqrt(B*B + 4 * A * C - 4 * C * log(P / 1e6))) - 273.15;\r
 149     }\r
 150     /* sa derivee */\r
 151     inline double DTsat( double P )\r
 152     {\r
 153       double A = 7.8130, B = 11209, C = 5.249e5, Tsk = Tsat(P) + 273.15;\r
 154       return Tsk * Tsk / ( P * sqrt(B*B + 4 * A * C - 4 * C * log(P / 1e6)));\r
 155     }\r
 156     /* pression de vapeur saturante */\r
 157     inline double Psat( double T )\r
 158     {\r
 159       double A = 7.8130, B = 11209, C = 5.249e5;\r
 160       return 1e6 * exp(A - B / Tk - C / (Tk * Tk));\r
 161     }\r
 162     /* sa derivee */\r
 163     inline double DPsat( double T )\r
 164     {\r
 165       double B = 11209, C = 5.249e5;\r
 166       return (B / (Tk * Tk) + 2 * C / (Tk * Tk * Tk)) * Psat(T);\r
 167     }\r
 168     /* enthalpie massique de saturation */\r
 169     inline double Hsat( double P )\r
 170     {\r
 171       return HL(Tsat(P), P);\r
 172     }\r
 173     /* sa derivee */\r
 174     inline double DHsat( double P )\r
 175     {\r
 176       return DTsat(P) * DTHL(Tsat(P), P) + DPHL(Tsat(P), P);\r
 177     }\r
 178     /* chaleur latente, a prendre a Tsat */\r
 179     inline double Lvap( double P )\r
 180     {\r
 181       double Tc = 2503.7, Tsk = Tsat(P) + 273.15;\r
 182       return 3.9337e5 * ( 1 - Tsk / Tc) + 4.3986e6 * pow( 1 - Tsk / Tc, .29302);\r
 183     }\r
 184     \r
 185     /* sa derivee */\r
 186     inline double DLvap( double P )\r
 187     {\r
 188       double Tc = 2503.7, Tsk = Tsat(P) + 273.15;\r
 189       return DTsat(P) * (-3.9337e5 / Tc - 4.3986e6  * .29302 * pow( 1 - Tsk / Tc, .29302 - 1) / Tc);\r
 190     }\r
 191     \r
 192     /* densite de la vapeur : (gaz parfait) * f1(Tsat) * f2(DTsat)*/\r
 193     #define  f1(Ts) (2.49121 + Ts * (-5.53796e-3 + Ts * (7.5465e-6 + Ts * (-4.20217e-9 + Ts * 8.59212e-13))))\r
 194     #define Df1(Ts) (-5.53796e-3 + Ts * (2 * 7.5465e-6 + Ts * (-3 * 4.20217e-9 + Ts * 4 * 8.59212e-13)))\r
 195     #define  f2(dT) (1 + dT * (-5e-4 + dT * (6.25e-7 - dT * 4.1359e-25)))\r
 196     #define Df2(dT) (-5e-4 + dT * (2 * 6.25e-7 - dT * 3 * 4.1359e-25))\r
 197     inline double RhoV( double T, double P )\r
 198     {\r
 199       double Tsk = Tsat(P) + 273.15, dT = Tk - Tsk;\r
 200       return P * 2.7650313e-3 * f1(Tsk) * f2(dT) / Tk;\r
 201     }\r
 202     /* ses derivees */\r
 203     inline double DTRhoV( double T, double P )\r
 204     {\r
 205       double Tsk = Tsat(P) + 273.15, dT = Tk - Tsk;\r
 206       return P * 2.7650313e-3 * f1(Tsk) * (Df2(dT) - f2(dT) / Tk) / Tk;\r
 207     }\r
 208     inline double DPRhoV( double T, double P )\r
 209     {\r
 210       double Tsk = Tsat(P) + 273.15, dT = Tk - Tsk;\r
 211       return 2.7650313e-3 * (f1(Tsk) * f2(dT) + P * DTsat(P) * (Df1(Tsk) * f2(dT) - f1(Tsk) * Df2(dT))) / Tk;\r
 212     }\r
 213     #undef f1\r
 214     #undef Df1\r
 215     #undef f2\r
 216     #undef Df2\r
 217     /* inverse de la densite de la vapeur */\r
 218     inline double IRhoV( double T, double P )\r
 219     {\r
 220       return 1 / RhoV(T, P);\r
 221     }\r
 222     /* ses derivees */\r
 223     inline double DTIRhoV( double T, double P )\r
 224     {\r
 225       return -DTRhoV(T, P) / (RhoV(T, P) * RhoV(T, P));\r
 226     }\r
 227     inline double DPIRhoV( double T, double P )\r
 228     {\r
 229       return -DPRhoV(T, P) / (RhoV(T, P) * RhoV(T, P));\r
 230     }\r
 231     /* enthalpie de la vapeur */\r
 232     #define  CpVs(Ts) (-1223.89 + Ts * (14.1191 + Ts * (-1.62025e-2 + Ts * 5.76923e-6)))\r
 233     #define DCpVs(Ts) (14.1191 + Ts * (- 2 * 1.62025e-2 + Ts * 3 * 5.76923e-6))\r
 234     inline double HV( double T, double P)\r
 235     {\r
 236       double Ts = Tsat(P), dT = T - Ts, Cp0 = 910, k = -4.6e-3;\r
 237       return HL(Ts, P) + Lvap(P) + Cp0 * dT + (Cp0 - CpVs(Ts)) * (1 - exp(k * dT)) / k;\r
 238     }\r
 239     /* ses derivees */\r
 240     inline double DTHV( double T, double P)\r
 241     {\r
 242         double Ts = Tsat(P), dT = T - Ts, Cp0 = 910, k = -4.6e-3;\r
 243         return Cp0 - (Cp0 - CpVs(Ts)) * exp(k * dT);\r
 244     }\r
 245     inline double DPHV( double T, double P)\r
 246     {\r
 247         double Ts = Tsat(P), dT = T - Ts, Cp0 = 910, k = -4.6e-3;\r
 248         return DPHL(Ts, P) + DLvap(P) + DTsat(P) * (DTHL(Ts, P) - Cp0 - DCpVs(Ts) * (1 - exp(k * dT)) / k + (Cp0 - CpVs(Ts)) * exp(k * dT));\r
 249     }\r
 250     #undef CpVs\r
 251     #undef DCpVs\r
 252     /* energie volumique de la vapeur */\r
 253     inline double EV( double T, double P)\r
 254     {\r
 255       return RhoV(T, P) * HV(T, P);\r
 256     }\r
 257     /* ses derivees */\r
 258     inline double DTEV(double T, double P)\r
 259     {\r
 260       return DTRhoV(T, P) * HV(T, P) + RhoV(T, P) * DTHV(T, P);\r
 261     }\r
 262     inline double DPEV(double T, double P)\r
 263     {\r
 264       return DPRhoV(T, P) * HV(T, P) + RhoV(T, P) * DPHV(T, P);\r
 265     }\r
 266     \r
 267     /* conductivite de la vapeur */\r
 268     inline double LambdaV( double T )\r
 269     {\r
 270       return 0.045;\r
 271     }\r
 272     \r
 273     /* viscosite de la vapeur */\r
 274     inline double MuV( double T )\r
 275     {\r
 276       return 2.5e-6 + 1.5e-8 * Tk;\r
 277     }\r
 278     \r
 279     /* Nusselt de la vapeur -> Dittus/ Boetler */\r
 280     inline double NuV(double T, double P, double w, double Dh)\r
 281     {\r
 282       double Re = RhoV(T, P) * dabs(w) * Dh / MuV(T), Pr = MuV(T) * DTHV(T, P) / LambdaV(T);\r
 283       return 0.023 * pow(Re, 0.8) * pow(Pr, 0.4);\r
 284     }\r
 285     \r
 286     inline double Rho( double T, double x, double P )\r
 287     {\r
 288       return 1 / ( (1-x)/RhoL(T, P) + x / RhoV(T, P) );\r
 289     }\r
 290     \r
 291     #undef Tk\r
 292     #undef Ksil\r
 293 };\r
 294 #endif\r