Cantera  2.0
WaterProps.cpp
Go to the documentation of this file.
1 /**
2  * @file WaterProps.cpp
3  */
4 /*
5  * Copyright (2006) Sandia Corporation. Under the terms of
6  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
7  * U.S. Government retains certain rights in this software.
8  */
9 
11 #include "cantera/base/ctml.h"
15 
16 #include <cmath>
17 
18 namespace Cantera
19 {
20 
21 
22 /*
23  * default constructor -> object owns its own water evaluator
24  */
26  m_waterIAPWS(0),
27  m_own_sub(false)
28 {
30  m_own_sub = true;
31 }
32 
33 /*
34  * constructor -> object in slave mode, It doesn't own its
35  * own water evaluator.
36  */
38  m_waterIAPWS(0),
39  m_own_sub(false)
40 {
41  if (wptr) {
42  m_waterIAPWS = wptr->getWater();
43  m_own_sub = false;
44  } else {
46  m_own_sub = true;
47  }
48 }
49 
51  m_waterIAPWS(0),
52  m_own_sub(false)
53 {
54  if (waterIAPWS) {
55  m_waterIAPWS = waterIAPWS;
56  m_own_sub = false;
57  } else {
59  m_own_sub = true;
60  }
61 }
62 
63 /**
64  * Copy constructor
65  */
67  m_waterIAPWS(0),
68  m_own_sub(false)
69 {
70  *this = b;
71 }
72 
73 /**
74  * Destructor
75  */
77 {
78  if (m_own_sub) {
79  delete m_waterIAPWS;
80  }
81 }
82 
83 /**
84  * Assignment operator
85  */
87 {
88  if (&b == this) {
89  return *this;
90  }
91 
92  if (m_own_sub) {
93  if (m_waterIAPWS) {
94  delete m_waterIAPWS;
95  m_waterIAPWS = 0;
96  }
97  }
98  if (b.m_own_sub) {
100  m_own_sub = true;
101  } else {
103  m_own_sub = false;
104  }
105 
106  return *this;
107 }
108 
109 // Simple calculation of water density at atmospheric pressure.
110 // Valid up to boiling point.
111 /*
112  * This formulation has no dependence on the pressure and shouldn't
113  * be used where accuracy is needed.
114  *
115  * @param T temperature in kelvin
116  * @param P Pressure in pascal
117  * @param ifunc changes what's returned
118  *
119  * @return value returned depends on ifunc value:
120  * ifunc = 0 Returns the density in kg/m^3
121  * ifunc = 1 returns the derivative of the density wrt T.
122  * ifunc = 2 returns the 2nd derivative of the density wrt T
123  * ifunc = 3 returns the derivative of the density wrt P.
124  *
125  * Verification:
126  * Agrees with the CRC values (6-10) for up to 4 sig digits.
127  *
128  * units = returns density in kg m-3.
129  */
130 doublereal WaterProps::density_T(doublereal T, doublereal P, int ifunc)
131 {
132  doublereal Tc = T - 273.15;
133  const doublereal U1 = 288.9414;
134  const doublereal U2 = 508929.2;
135  const doublereal U3 = 68.12963;
136  const doublereal U4 = -3.9863;
137 
138  doublereal tmp1 = Tc + U1;
139  doublereal tmp4 = Tc + U4;
140  doublereal t4t4 = tmp4 * tmp4;
141  doublereal tmp3 = Tc + U3;
142  doublereal rho = 1000. * (1.0 - tmp1*t4t4/(U2 * tmp3));
143 
144  /*
145  * Impose an ideal gas lower bound on rho. We need this
146  * to ensure positivity of rho, even though it is
147  * grossly unrepresentative.
148  */
149  doublereal rhomin = P / (GasConstant * T);
150  if (rho < rhomin) {
151  rho = rhomin;
152  if (ifunc == 1) {
153  doublereal drhodT = - rhomin / T;
154  return drhodT;
155  } else if (ifunc == 3) {
156  doublereal drhodP = rhomin / P;
157  return drhodP;
158  } else if (ifunc == 2) {
159  doublereal d2rhodT2 = 2.0 * rhomin / (T * T);
160  return d2rhodT2;
161  }
162  }
163 
164  if (ifunc == 1) {
165  doublereal drhodT = 1000./U2 * (
166  - tmp4 * tmp4 / (tmp3)
167  - tmp1 * 2 * tmp4 / (tmp3)
168  + tmp1 * t4t4 / (tmp3*tmp3)
169  );
170  return drhodT;
171  } else if (ifunc == 3) {
172  return 0.0;
173  } else if (ifunc == 2) {
174  doublereal t3t3 = tmp3 * tmp3;
175  doublereal d2rhodT2 = 1000./U2 *
176  ((-4.0*tmp4-2.0*tmp1)/tmp3 +
177  (2.0*t4t4 + 4.0*tmp1*tmp4)/t3t3
178  - 2.0*tmp1 * t4t4/(t3t3*tmp3));
179  return d2rhodT2;
180  }
181 
182  return rho;
183 }
184 
185 // Bradley-Pitzer equation for the dielectric constant
186 // of water as a function of temperature and pressure.
187 /*!
188  * Returns the dimensionless relative dielectric constant
189  * and its derivatives.
190  *
191  * ifunc = 0 value
192  * ifunc = 1 Temperature derivative
193  * ifunc = 2 second temperature derivative
194  * ifunc = 3 return pressure first derivative
195  *
196  * Range of validity 0 to 350C, 0 to 1 kbar pressure
197  *
198  * @param T temperature (kelvin)
199  * @param P_pascal pressure in pascal
200  * @param ifunc changes what's returned from the function
201  *
202  * @return Depends on the value of ifunc:
203  * ifunc = 0 return value
204  * ifunc = 1 return temperature derivative
205  * ifunc = 2 return second temperature derivative
206  * ifunc = 3 return pressure first derivative
207  *
208  * Validation:
209  * Numerical experiments indicate that this function agrees with
210  * the Archer and Wang data in the CRC p. 6-10 to all 4 significant
211  * digits shown (0 to 100C).
212  *
213  * value at 25C, relEps = 78.38
214  *
215  */
216 doublereal WaterProps::relEpsilon(doublereal T, doublereal P_pascal,
217  int ifunc)
218 {
219  const doublereal U1 = 3.4279E2;
220  const doublereal U2 = -5.0866E-3;
221  const doublereal U3 = 9.4690E-7;
222  const doublereal U4 = -2.0525;
223  const doublereal U5 = 3.1159E3;
224  const doublereal U6 = -1.8289E2;
225  const doublereal U7 = -8.0325E3;
226  const doublereal U8 = 4.2142E6;
227  const doublereal U9 = 2.1417;
228  doublereal T2 = T * T;
229 
230  doublereal eps1000 = U1 * exp(U2 * T + U3 * T2);
231  doublereal C = U4 + U5/(U6 + T);
232  doublereal B = U7 + U8/T + U9 * T;
233 
234  doublereal Pbar = P_pascal * 1.0E-5;
235  doublereal tmpBpar = B + Pbar;
236  doublereal tmpB1000 = B + 1000.0;
237  doublereal ltmp = log(tmpBpar/tmpB1000);
238  doublereal epsRel = eps1000 + C * ltmp;
239 
240  if (ifunc == 1 || ifunc == 2) {
241  doublereal tmpC = U6 + T;
242  doublereal dCdT = - U5/(tmpC * tmpC);
243 
244  doublereal dBdT = - U8/(T * T) + U9;
245 
246  doublereal deps1000dT = eps1000 * (U2 + 2.0 * U3 * T);
247 
248  doublereal dltmpdT = (dBdT/tmpBpar - dBdT/tmpB1000);
249  if (ifunc == 1) {
250  doublereal depsReldT = deps1000dT + dCdT * ltmp + C * dltmpdT;
251  return depsReldT;
252  }
253  doublereal T3 = T2 * T;
254  doublereal d2CdT2 = - 2.0 * dCdT / tmpC;
255  doublereal d2BdT2 = 2.0 * U8 / (T3);
256 
257  doublereal d2ltmpdT2 = (d2BdT2*(1.0/tmpBpar - 1.0/tmpB1000) +
258  dBdT*dBdT*(1.0/(tmpB1000*tmpB1000) - 1.0/(tmpBpar*tmpBpar)));
259 
260  doublereal d2eps1000dT2 = (deps1000dT * (U2 + 2.0 * U3 * T) + eps1000 * (2.0 * U3));
261 
262  if (ifunc == 2) {
263  doublereal d2epsReldT2 = (d2eps1000dT2 + d2CdT2 * ltmp + 2.0 * dCdT * dltmpdT
264  + C * d2ltmpdT2);
265  return d2epsReldT2;
266  }
267  }
268  if (ifunc == 3) {
269  doublereal dltmpdP = 1.0E-5 / tmpBpar;
270  doublereal depsReldP = C * dltmpdP;
271  return depsReldP;
272  }
273 
274  return epsRel;
275 }
276 
277 
278 doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc)
279 {
280  doublereal psat = satPressure(T);
281  doublereal P;
282  if (psat > P_input) {
283  //printf("ADebye WARNING: p_input < psat: %g %g\n",
284  // P_input, psat);
285  P = psat;
286  } else {
287  P = P_input;
288  }
289  doublereal epsRelWater = relEpsilon(T, P, 0);
290  //printf("releps calc = %g, compare to 78.38\n", epsRelWater);
291  //doublereal B_Debye = 3.28640E9;
292 
293  doublereal epsilon = epsilon_0 * epsRelWater;
294  doublereal dw = density_IAPWS(T, P);
295  doublereal tmp = sqrt(2.0 * Avogadro * dw / 1000.);
296  doublereal tmp2 = ElectronCharge * ElectronCharge * Avogadro /
297  (epsilon * GasConstant * T);
298  doublereal tmp3 = tmp2 * sqrt(tmp2);
299  doublereal A_Debye = tmp * tmp3 / (8.0 * Pi);
300 
301 
302  /*
303  * dAdT = - 3/2 Ad/T + 1/2 Ad/dw d(dw)/dT - 3/2 Ad/eps d(eps)/dT
304  *
305  * dAdT = - 3/2 Ad/T - 1/2 Ad/Vw d(Vw)/dT - 3/2 Ad/eps d(eps)/dT
306  */
307  if (ifunc == 1 || ifunc == 2) {
308  doublereal dAdT = - 1.5 * A_Debye / T;
309 
310  doublereal depsRelWaterdT = relEpsilon(T, P, 1);
311  dAdT -= A_Debye * (1.5 * depsRelWaterdT / epsRelWater);
312 
313  //int methodD = 1;
314  //doublereal ddwdT = density_T_new(T, P, 1);
315  // doublereal contrib1 = A_Debye * (0.5 * ddwdT / dw);
316 
317  /*
318  * calculate d(lnV)/dT _constantP, i.e., the cte
319  */
320  doublereal cte = coeffThermalExp_IAPWS(T, P);
321  doublereal contrib2 = - A_Debye * (0.5 * cte);
322 
323  //dAdT += A_Debye * (0.5 * ddwdT / dw);
324  dAdT += contrib2;
325 
326 #ifdef DEBUG_HKM
327  //printf("dAdT = %g, contrib1 = %g, contrib2 = %g\n",
328  // dAdT, contrib1, contrib2);
329 #endif
330 
331  if (ifunc == 1) {
332  return dAdT;
333  }
334 
335  if (ifunc == 2) {
336  /*
337  * Get the second derivative of the dielectric constant wrt T
338  * -> we will take each of the terms in dAdT and differentiate
339  * it again.
340  */
341  doublereal d2AdT2 = 1.5 / T * (A_Debye/T - dAdT);
342 
343  doublereal d2epsRelWaterdT2 = relEpsilon(T, P, 2);
344 
345  //doublereal dT = -0.01;
346  //doublereal TT = T + dT;
347  //doublereal depsRelWaterdTdel = relEpsilon(TT, P, 1);
348  //doublereal d2alt = (depsRelWaterdTdel- depsRelWaterdT ) / dT;
349  //printf("diff %g %g\n",d2epsRelWaterdT2, d2alt);
350  // HKM -> checks out, i.e., they are the same.
351 
352  d2AdT2 += 1.5 * (- dAdT * depsRelWaterdT / epsRelWater
353  - A_Debye / epsRelWater *
354  (d2epsRelWaterdT2 - depsRelWaterdT * depsRelWaterdT / epsRelWater));
355 
356  doublereal deltaT = -0.1;
357  doublereal Tdel = T + deltaT;
358  doublereal cte_del = coeffThermalExp_IAPWS(Tdel, P);
359  doublereal dctedT = (cte_del - cte) / Tdel;
360 
361 
362  //doublereal d2dwdT2 = density_T_new(T, P, 2);
363 
364  doublereal contrib3 = 0.5 * (-(dAdT * cte) -(A_Debye * dctedT));
365  d2AdT2 += contrib3;
366 
367  return d2AdT2;
368  }
369  }
370  /*
371  * A_Debye = (1/(8 Pi)) sqrt(2 Na dw / 1000)
372  * (e e/(epsilon R T))^3/2
373  *
374  * dAdP = + 1/2 Ad/dw d(dw)/dP - 3/2 Ad/eps d(eps)/dP
375  *
376  * dAdP = - 1/2 Ad/Vw d(Vw)/dP - 3/2 Ad/eps d(eps)/dP
377  *
378  * dAdP = + 1/2 Ad * kappa - 3/2 Ad/eps d(eps)/dP
379  *
380  * where kappa = - 1/Vw d(Vw)/dP_T (isothermal compressibility)
381  */
382  if (ifunc == 3) {
383 
384  doublereal dAdP = 0.0;
385 
386  doublereal depsRelWaterdP = relEpsilon(T, P, 3);
387  dAdP -= A_Debye * (1.5 * depsRelWaterdP / epsRelWater);
388 
389  doublereal kappa = isothermalCompressibility_IAPWS(T,P);
390 
391  //doublereal ddwdP = density_T_new(T, P, 3);
392  dAdP += A_Debye * (0.5 * kappa);
393 
394  return dAdP;
395  }
396 
397  return A_Debye;
398 }
399 
400 doublereal WaterProps::satPressure(doublereal T)
401 {
402  doublereal pres = m_waterIAPWS->psat(T);
403  return pres;
404 }
405 
406 // Returns the density of water
407 /*
408  * This function sets the internal temperature and pressure
409  * of the underlying object at the same time.
410  *
411  * @param T Temperature (kelvin)
412  * @param P pressure (pascal)
413  */
414 doublereal WaterProps::density_IAPWS(doublereal temp, doublereal press)
415 {
416  doublereal dens = m_waterIAPWS->density(temp, press, WATER_LIQUID);
417  return dens;
418 }
419 
420 // Returns the density of water
421 /*
422  * This function uses the internal state of the
423  * underlying water object
424  */
425 doublereal WaterProps::density_IAPWS() const
426 {
427  doublereal dens = m_waterIAPWS->density();
428  return dens;
429 }
430 
431 doublereal WaterProps::coeffThermalExp_IAPWS(doublereal temp, doublereal press)
432 {
433  doublereal dens = m_waterIAPWS->density(temp, press, WATER_LIQUID);
434  if (dens < 0.0) {
435  throw CanteraError("WaterProps::coeffThermalExp_IAPWS",
436  "Unable to solve for density at T = " + fp2str(temp) + " and P = " + fp2str(press));
437  }
438  doublereal cte = m_waterIAPWS->coeffThermExp();
439  return cte;
440 }
441 
442 doublereal WaterProps::isothermalCompressibility_IAPWS(doublereal temp, doublereal press)
443 {
444  doublereal dens = m_waterIAPWS->density(temp, press, WATER_LIQUID);
445  if (dens < 0.0) {
446  throw CanteraError("WaterProps::isothermalCompressibility_IAPWS",
447  "Unable to solve for density at T = " + fp2str(temp) + " and P = " + fp2str(press));
448  }
449  doublereal kappa = m_waterIAPWS->isothermalCompressibility();
450  return kappa;
451 }
452 
453 
454 
455 
456 
457 // Parameters for the viscosityWater() function
458 
459 // \cond
460 const doublereal H[4] = {1.,
461  0.978197,
462  0.579829,
463  -0.202354
464  };
465 
466 const doublereal Hij[6][7] = {
467  { 0.5132047, 0.2151778, -0.2818107, 0.1778064, -0.04176610, 0., 0.},
468  { 0.3205656, 0.7317883, -1.070786 , 0.4605040, 0., -0.01578386, 0.},
469  { 0., 1.241044 , -1.263184 , 0.2340379, 0., 0., 0.},
470  { 0., 1.476783 , 0., -0.4924179, 0.1600435, 0., -0.003629481},
471  {-0.7782567, 0.0 , 0., 0. , 0., 0., 0.},
472  { 0.1885447, 0.0 , 0., 0. , 0., 0., 0.},
473 };
474 const doublereal TStar = 647.27; // Kelvin
475 const doublereal rhoStar = 317.763; // kg / m3
476 const doublereal presStar = 22.115E6; // Pa
477 const doublereal muStar = 55.071E-6; //Pa s
478 // \endcond
479 
480 // Returns the viscosity of water at the current conditions
481 // (kg/m/s)
482 /*
483  * This function calculates the value of the viscosity of pure
484  * water at the current T and P.
485  *
486  * The formulas used are from the paper
487  *
488  * J. V. Sengers, J. T. R. Watson, "Improved International
489  * Formulations for the Viscosity and Thermal Conductivity of
490  * Water Substance", J. Phys. Chem. Ref. Data, 15, 1291 (1986).
491  *
492  * The formulation is accurate for all temperatures and pressures,
493  * for steam and for water, even near the critical point.
494  * Pressures above 500 MPa and temperature above 900 C are suspect.
495  */
496 doublereal WaterProps::viscosityWater() const
497 {
498 
499  doublereal temp = m_waterIAPWS->temperature();
500  doublereal dens = m_waterIAPWS->density();
501 
502  //WaterPropsIAPWS *waterP = new WaterPropsIAPWS();
503  //m_waterIAPWS->setState_TR(temp, dens);
504  //doublereal pressure = m_waterIAPWS->pressure();
505  //printf("pressure = %g\n", pressure);
506  //dens = 18.02 * pressure / (GasConstant * temp);
507  //printf ("mod dens = %g\n", dens);
508 
509  doublereal rhobar = dens/rhoStar;
510  doublereal tbar = temp / TStar;
511  // doublereal pbar = pressure / presStar;
512 
513  doublereal tbar2 = tbar * tbar;
514  doublereal tbar3 = tbar2 * tbar;
515 
516  doublereal mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3);
517 
518  //printf("mu0bar = %g\n", mu0bar);
519  //printf("mu0 = %g\n", mu0bar * muStar);
520 
521  doublereal tfac1 = 1.0 / tbar - 1.0;
522  doublereal tfac2 = tfac1 * tfac1;
523  doublereal tfac3 = tfac2 * tfac1;
524  doublereal tfac4 = tfac3 * tfac1;
525  doublereal tfac5 = tfac4 * tfac1;
526 
527  doublereal rfac1 = rhobar - 1.0;
528  doublereal rfac2 = rfac1 * rfac1;
529  doublereal rfac3 = rfac2 * rfac1;
530  doublereal rfac4 = rfac3 * rfac1;
531  doublereal rfac5 = rfac4 * rfac1;
532  doublereal rfac6 = rfac5 * rfac1;
533 
534  doublereal sum = (Hij[0][0] + Hij[1][0]*tfac1 + Hij[4][0]*tfac4 + Hij[5][0]*tfac5 +
535  Hij[0][1]*rfac1 + Hij[1][1]*tfac1*rfac1 + Hij[2][1]*tfac2*rfac1 + Hij[3][1]*tfac3*rfac1 +
536  Hij[0][2]*rfac2 + Hij[1][2]*tfac1*rfac2 + Hij[2][2]*tfac2*rfac2 +
537  Hij[0][3]*rfac3 + Hij[1][3]*tfac1*rfac3 + Hij[2][3]*tfac2*rfac3 + Hij[3][3]*tfac3*rfac3 +
538  Hij[0][4]*rfac4 + Hij[3][4]*tfac3*rfac4 +
539  Hij[1][5]*tfac1*rfac5 + Hij[3][6]*tfac3*rfac6
540  );
541  doublereal mu1bar = std::exp(rhobar * sum);
542 
543  // Apply the near-critical point corrections if necessary
544  doublereal mu2bar = 1.0;
545  if ((tbar >= 0.9970) && tbar <= 1.0082) {
546  if ((rhobar >= 0.755) && (rhobar <= 1.290)) {
547  doublereal drhodp = 1.0 / m_waterIAPWS->dpdrho();
548  drhodp *= presStar / rhoStar;
549  doublereal xsi = rhobar * drhodp;
550  if (xsi >= 21.93) {
551  mu2bar = 0.922 * std::pow(xsi, 0.0263);
552  }
553  }
554  }
555 
556  doublereal mubar = mu0bar * mu1bar * mu2bar;
557 
558  return mubar * muStar;
559 }
560 
561 //! Returns the thermal conductivity of water at the current conditions
562 //! (W/m/K)
563 /*!
564  * This function calculates the value of the thermal conductivity of
565  * water at the current T and P.
566  *
567  * The formulas used are from the paper
568  * J. V. Sengers, J. T. R. Watson, "Improved International
569  * Formulations for the Viscosity and Thermal Conductivity of
570  * Water Substance", J. Phys. Chem. Ref. Data, 15, 1291 (1986).
571  *
572  * The formulation is accurate for all temperatures and pressures,
573  * for steam and for water, even near the critical point.
574  * Pressures above 500 MPa and temperature above 900 C are suspect.
575  */
577 {
578  static const doublereal Tstar = 647.27;
579  static const doublereal rhostar = 317.763;
580  static const doublereal lambdastar = 0.4945;
581  static const doublereal presstar = 22.115E6;
582  static const doublereal L[4] = {
583  1.0000,
584  6.978267,
585  2.599096,
586  -0.998254
587  };
588  static const doublereal Lji[6][5] = {
589  { 1.3293046, 1.7018363, 5.2246158, 8.7127675, -1.8525999},
590  {-0.40452437, -2.2156845, -10.124111, -9.5000611, 0.93404690},
591  { 0.24409490, 1.6511057, 4.9874687, 4.3786606, 0.0},
592  { 0.018660751, -0.76736002, -0.27297694, -0.91783782, 0.0},
593  {-0.12961068, 0.37283344, -0.43083393, 0.0, 0.0},
594  { 0.044809953, -0.11203160, 0.13333849, 0.0, 0.0},
595  };
596 
597  doublereal temp = m_waterIAPWS->temperature();
598  doublereal dens = m_waterIAPWS->density();
599 
600  doublereal rhobar = dens/rhostar;
601  doublereal tbar = temp / Tstar;
602  doublereal tbar2 = tbar * tbar;
603  doublereal tbar3 = tbar2 * tbar;
604  doublereal lambda0bar = sqrt(tbar) / (L[0] + L[1]/tbar + L[2]/tbar2 + L[3]/tbar3);
605 
606  //doublereal lambdagas = lambda0bar * lambdastar * 1.0E3;
607 
608  doublereal tfac1 = 1.0 / tbar - 1.0;
609  doublereal tfac2 = tfac1 * tfac1;
610  doublereal tfac3 = tfac2 * tfac1;
611  doublereal tfac4 = tfac3 * tfac1;
612 
613  doublereal rfac1 = rhobar - 1.0;
614  doublereal rfac2 = rfac1 * rfac1;
615  doublereal rfac3 = rfac2 * rfac1;
616  doublereal rfac4 = rfac3 * rfac1;
617  doublereal rfac5 = rfac4 * rfac1;
618 
619  doublereal sum = (Lji[0][0] + Lji[0][1]*tfac1 + Lji[0][2]*tfac2 + Lji[0][3]*tfac3 + Lji[0][4]*tfac4 +
620  Lji[1][0]*rfac1 + Lji[1][1]*tfac1*rfac1 + Lji[1][2]*tfac2*rfac1 + Lji[1][3]*tfac3*rfac1 + Lji[1][4]*tfac4*rfac1 +
621  Lji[2][0]*rfac2 + Lji[2][1]*tfac1*rfac2 + Lji[2][2]*tfac2*rfac2 + Lji[2][3]*tfac3*rfac2 +
622  Lji[3][0]*rfac3 + Lji[3][1]*tfac1*rfac3 + Lji[3][2]*tfac2*rfac3 + Lji[3][3]*tfac3*rfac3 +
623  Lji[4][0]*rfac4 + Lji[4][1]*tfac1*rfac4 + Lji[4][2]*tfac2*rfac4 +
624  Lji[5][0]*rfac5 + Lji[5][1]*tfac1*rfac5 + Lji[5][2]*tfac2*rfac5
625  );
626  doublereal lambda1bar = exp(rhobar * sum);
627 
628  doublereal mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3);
629 
630  doublereal tfac5 = tfac4 * tfac1;
631  doublereal rfac6 = rfac5 * rfac1;
632 
633  sum = (Hij[0][0] + Hij[1][0]*tfac1 + Hij[4][0]*tfac4 + Hij[5][0]*tfac5 +
634  Hij[0][1]*rfac1 + Hij[1][1]*tfac1*rfac1 + Hij[2][1]*tfac2*rfac1 + Hij[3][1]*tfac3*rfac1 +
635  Hij[0][2]*rfac2 + Hij[1][2]*tfac1*rfac2 + Hij[2][2]*tfac2*rfac2 +
636  Hij[0][3]*rfac3 + Hij[1][3]*tfac1*rfac3 + Hij[2][3]*tfac2*rfac3 + Hij[3][3]*tfac3*rfac3 +
637  Hij[0][4]*rfac4 + Hij[3][4]*tfac3*rfac4 +
638  Hij[1][5]*tfac1*rfac5 + Hij[3][6]*tfac3*rfac6
639  );
640  doublereal mu1bar = std::exp(rhobar * sum);
641 
642  doublereal t2r2 = tbar * tbar / (rhobar * rhobar);
643  doublereal drhodp = 1.0 / m_waterIAPWS->dpdrho();
644  drhodp *= presStar / rhoStar;
645  doublereal xsi = rhobar * drhodp;
646  doublereal xsipow = std::pow(xsi, 0.4678);
647  doublereal rho1 = rhobar - 1.;
648  doublereal rho2 = rho1 * rho1;
649  doublereal rho4 = rho2 * rho2;
650  doublereal temp2 = (tbar - 1.0) * (tbar - 1.0);
651 
652  /*
653  * beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho
654  *
655  * Note for ideal gases this is equal to one.
656  *
657  * beta = delta (phi0_d() + phiR_d())
658  * - tau delta (phi0_dt() + phiR_dt())
659  */
660  doublereal beta = m_waterIAPWS->coeffPresExp();
661 
662  doublereal dpdT_const_rho = beta * GasConstant * dens / 18.015268;
663  dpdT_const_rho *= Tstar / presstar;
664 
665  doublereal lambda2bar = 0.0013848 / (mu0bar * mu1bar) * t2r2 * dpdT_const_rho * dpdT_const_rho *
666  xsipow * sqrt(rhobar) * exp(-18.66*temp2 - rho4);
667 
668  doublereal lambda = (lambda0bar * lambda1bar + lambda2bar) * lambdastar;
669  return lambda;
670 }
671 
672 
673 }