My Project
H2OAirFluidSystem.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_H2O_AIR_SYSTEM_HPP
28 #define OPM_H2O_AIR_SYSTEM_HPP
29 
30 #include "BaseFluidSystem.hpp"
31 #include "NullParameterCache.hpp"
32 
38 
40 
41 #include <iostream>
42 #include <cassert>
43 
44 namespace Opm {
45 
55 template <class Scalar,
56  //class H2Otype = SimpleH2O<Scalar>,
57  class H2Otype = TabulatedComponent<Scalar, H2O<Scalar> >>
59  : public BaseFluidSystem<Scalar, H2OAirFluidSystem<Scalar, H2Otype> >
60 {
63  typedef ::Opm::IdealGas<Scalar> IdealGas;
64 
65 public:
66  template <class Evaluation>
67  struct ParameterCache : public NullParameterCache<Evaluation>
68  {};
69 
71  typedef H2Otype H2O;
73  typedef ::Opm::Air<Scalar> Air;
74 
76  static const int numPhases = 2;
77 
79  static const int liquidPhaseIdx = 0;
81  static const int gasPhaseIdx = 1;
82 
84  static const char* phaseName(unsigned phaseIdx)
85  {
86  switch (phaseIdx) {
87  case liquidPhaseIdx: return "liquid";
88  case gasPhaseIdx: return "gas";
89  };
90  throw std::logic_error("Invalid phase index "+std::to_string(phaseIdx));
91  }
92 
94  static bool isLiquid(unsigned phaseIdx)
95  {
96  //assert(0 <= phaseIdx && phaseIdx < numPhases);
97  return phaseIdx != gasPhaseIdx;
98  }
99 
101  static bool isCompressible(unsigned phaseIdx)
102  {
103  //assert(0 <= phaseIdx && phaseIdx < numPhases);
104  return (phaseIdx == gasPhaseIdx)
105  // ideal gases are always compressible
106  ? true
107  :
108  // the water component decides for the liquid phase...
110  }
111 
113  static bool isIdealGas(unsigned phaseIdx)
114  {
115  return
116  (phaseIdx == gasPhaseIdx)
118  : false;
119  }
120 
122  static bool isIdealMixture(unsigned /*phaseIdx*/)
123  {
124  //assert(0 <= phaseIdx && phaseIdx < numPhases);
125  // we assume Henry's and Rault's laws for the water phase and
126  // and no interaction between gas molecules of different
127  // components, so all phases are ideal mixtures!
128  return true;
129  }
130 
131  /****************************************
132  * Component related static parameters
133  ****************************************/
134 
136  static const int numComponents = 2;
137 
139  static const int H2OIdx = 0;
141  static const int AirIdx = 1;
142 
144  static const char* componentName(unsigned compIdx)
145  {
146  switch (compIdx)
147  {
148  case H2OIdx: return H2O::name();
149  case AirIdx: return Air::name();
150  };
151  throw std::logic_error("Invalid component index "+std::to_string(compIdx));
152  }
153 
155  static Scalar molarMass(unsigned compIdx)
156  {
157  return
158  (compIdx == H2OIdx)
159  ? H2O::molarMass()
160  : (compIdx == AirIdx)
161  ? Air::molarMass()
162  : 1e30;
163  }
164 
165 
171  static Scalar criticalTemperature(unsigned compIdx)
172  {
173  return
174  (compIdx == H2OIdx)
176  : (compIdx == AirIdx)
178  : 1e30;
179  }
180 
186  static Scalar criticalPressure(unsigned compIdx)
187  {
188  return
189  (compIdx == H2OIdx)
191  : (compIdx == AirIdx)
193  : 1e30;
194  }
195 
201  static Scalar acentricFactor(unsigned compIdx)
202  {
203  return
204  (compIdx == H2OIdx)
206  : (compIdx == AirIdx)
208  : 1e30;
209  }
210 
211  /****************************************
212  * thermodynamic relations
213  ****************************************/
214 
221  static void init()
222  {
223  if (H2O::isTabulated)
224  init(/*tempMin=*/273.15,
225  /*tempMax=*/623.15,
226  /*numTemp=*/50,
227  /*pMin=*/-10,
228  /*pMax=*/20e6,
229  /*numP=*/50);
230  }
231 
243  static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
244  Scalar pressMin, Scalar pressMax, unsigned nPress)
245  {
246  if (H2O::isTabulated) {
247  H2O::init(tempMin, tempMax, nTemp,
248  pressMin, pressMax, nPress);
249  }
250  }
251 
253  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
254  static LhsEval density(const FluidState& fluidState,
255  const ParameterCache<ParamCacheEval>& /*paramCache*/,
256  unsigned phaseIdx)
257  {
258  assert(phaseIdx < numPhases);
259 
260  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
261  LhsEval p;
262  if (isCompressible(phaseIdx))
263  p = decay<LhsEval>(fluidState.pressure(phaseIdx));
264  else {
265  // random value which will hopefully cause things to blow
266  // up if it is used in a calculation!
267  p = - 1e30;
268  Valgrind::SetUndefined(p);
269  }
270 
271 
272  LhsEval sumMoleFrac = 0;
273  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
274  sumMoleFrac += decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
275 
276  if (phaseIdx == liquidPhaseIdx)
277  {
278  // assume ideal mixture: Molecules of one component don't discriminate
279  // between their own kind and molecules of the other component.
280  const LhsEval& clH2O = H2O::liquidDensity(T, p)/H2O::molarMass();
281 
282  const auto& xlH2O = decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
283  const auto& xlAir = decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, AirIdx));
284 
285  return clH2O*(H2O::molarMass()*xlH2O + Air::molarMass()*xlAir)/sumMoleFrac;
286  }
287  else if (phaseIdx == gasPhaseIdx)
288  {
289  LhsEval partialPressureH2O =
290  decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))
291  *decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
292 
293  LhsEval partialPressureAir =
294  decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, AirIdx))
295  *decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
296 
297  return H2O::gasDensity(T, partialPressureH2O) + Air::gasDensity(T, partialPressureAir);
298  }
299  throw std::logic_error("Invalid phase index "+std::to_string(phaseIdx));
300  }
301 
303  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
304  static LhsEval viscosity(const FluidState& fluidState,
305  const ParameterCache<ParamCacheEval>& /*paramCache*/,
306  unsigned phaseIdx)
307  {
308  assert(phaseIdx < numPhases);
309 
310  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
311  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
312 
313  if (phaseIdx == liquidPhaseIdx)
314  {
315  // assume pure water for the liquid phase
316  // TODO: viscosity of mixture
317  // couldn't find a way to solve the mixture problem
318  return H2O::liquidViscosity(T, p);
319  }
320  else if (phaseIdx == gasPhaseIdx)
321  {
322  /* Wilke method. See:
323  *
324  * See: R. Reid, et al.: The Properties of Gases and Liquids,
325  * 4th edition, McGraw-Hill, 1987, 407-410 or
326  * 5th edition, McGraw-Hill, 2000, p. 9.21/22
327  */
328  LhsEval muResult = 0;
329  const LhsEval mu[numComponents] = {
331  Air::gasViscosity(T, p)
332  };
333 
334  for (unsigned i = 0; i < numComponents; ++i) {
335  LhsEval divisor = 0;
336  for (unsigned j = 0; j < numComponents; ++j) {
337  LhsEval phiIJ =
338  1 +
339  sqrt(mu[i]/mu[j]) * // 1 + (mu[i]/mu[j]^1/2
340  std::pow(molarMass(j)/molarMass(i), 1./4.0); // (M[i]/M[j])^1/4
341 
342  phiIJ *= phiIJ;
343  phiIJ /= std::sqrt(8*(1 + molarMass(i)/molarMass(j)));
344  divisor += decay<LhsEval>(fluidState.moleFraction(phaseIdx, j))*phiIJ;
345  }
346  const auto& xAlphaI = decay<LhsEval>(fluidState.moleFraction(phaseIdx, i));
347  muResult += xAlphaI*mu[i]/divisor;
348  }
349  return muResult;
350  }
351  throw std::logic_error("Invalid phase index "+std::to_string(phaseIdx));
352  }
353 
355  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
356  static LhsEval fugacityCoefficient(const FluidState& fluidState,
357  const ParameterCache<ParamCacheEval>& /*paramCache*/,
358  unsigned phaseIdx,
359  unsigned compIdx)
360  {
361  assert(phaseIdx < numPhases);
362  assert(compIdx < numComponents);
363 
364  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
365  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
366 
367  if (phaseIdx == liquidPhaseIdx) {
368  if (compIdx == H2OIdx)
369  return H2O::vaporPressure(T)/p;
370  return BinaryCoeff::H2O_Air::henry(T)/p;
371  }
372 
373  // for the gas phase, assume an ideal gas when it comes to
374  // fugacity (-> fugacity == partial pressure)
375  return 1.0;
376  }
377 
379  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
380  static LhsEval binaryDiffusionCoefficient(const FluidState& fluidState,
381  const ParameterCache<ParamCacheEval>& /*paramCache*/,
382  unsigned phaseIdx,
383  unsigned /*compIdx*/)
384  {
385  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
386  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
387 
388  if (phaseIdx == liquidPhaseIdx)
390 
391  assert(phaseIdx == gasPhaseIdx);
393  }
394 
396  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
397  static LhsEval enthalpy(const FluidState& fluidState,
398  const ParameterCache<ParamCacheEval>& /*paramCache*/,
399  unsigned phaseIdx)
400  {
401  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
402  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
403  Valgrind::CheckDefined(T);
404  Valgrind::CheckDefined(p);
405 
406  if (phaseIdx == liquidPhaseIdx)
407  {
408  // TODO: correct way to deal with the solutes???
409  return H2O::liquidEnthalpy(T, p);
410  }
411 
412  else if (phaseIdx == gasPhaseIdx)
413  {
414  LhsEval result = 0.0;
415  result +=
416  H2O::gasEnthalpy(T, p) *
417  decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
418 
419  result +=
420  Air::gasEnthalpy(T, p) *
421  decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, AirIdx));
422  return result;
423  }
424  throw std::logic_error("Invalid phase index "+std::to_string(phaseIdx));
425  }
426 
428  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
429  static LhsEval thermalConductivity(const FluidState& fluidState,
430  const ParameterCache<ParamCacheEval>& /*paramCache*/,
431  unsigned phaseIdx)
432  {
433  assert(phaseIdx < numPhases);
434 
435  const LhsEval& temperature =
436  decay<LhsEval>(fluidState.temperature(phaseIdx));
437  const LhsEval& pressure =
438  decay<LhsEval>(fluidState.pressure(phaseIdx));
439 
440  if (phaseIdx == liquidPhaseIdx)
441  return H2O::liquidThermalConductivity(temperature, pressure);
442  else { // gas phase
443  const LhsEval& lambdaDryAir = Air::gasThermalConductivity(temperature, pressure);
444 
445  const LhsEval& xAir =
446  decay<LhsEval>(fluidState.moleFraction(phaseIdx, AirIdx));
447  const LhsEval& xH2O =
448  decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
449  LhsEval lambdaAir = xAir*lambdaDryAir;
450 
451  // Assuming Raoult's, Daltons law and ideal gas
452  // in order to obtain the partial density of water in the air phase
453  LhsEval partialPressure = pressure*xH2O;
454 
455  LhsEval lambdaH2O =
456  xH2O*H2O::gasThermalConductivity(temperature, partialPressure);
457  return lambdaAir + lambdaH2O;
458  }
459  }
460 };
461 
462 } // namespace Opm
463 
464 #endif
A simple class implementing the fluid properties of air.
The base class for all fluid systems.
Material properties of pure water .
Binary coefficients for water and nitrogen.
Relations valid for an ideal gas.
A parameter cache which does nothing.
A generic class which tabulates all thermodynamic properties of a given component.
Some templates to wrap the valgrind client request macros.
static Evaluation gasThermalConductivity(const Evaluation &, const Evaluation &)
Specific heat conductivity of steam .
Definition: Air.hpp:217
static Scalar criticalPressure()
Returns the critical pressure of .
Definition: Air.hpp:92
static Evaluation gasDensity(const Evaluation &temperature, const Evaluation &pressure)
The density of at a given pressure and temperature [kg/m^3].
Definition: Air.hpp:102
static Scalar molarMass()
The molar mass in of .
Definition: Air.hpp:80
static bool gasIsIdeal()
Returns true iff the gas phase is assumed to be ideal.
Definition: Air.hpp:72
static const char * name()
A human readable name for the .
Definition: Air.hpp:60
static Evaluation gasEnthalpy(const Evaluation &temperature, const Evaluation &)
Specific enthalpy of liquid water with 273.15 K as basis.
Definition: Air.hpp:180
static Scalar criticalTemperature()
Returns the critical temperature of .
Definition: Air.hpp:86
static Evaluation gasViscosity(const Evaluation &temperature, const Evaluation &)
The dynamic viscosity of at a given pressure and temperature.
Definition: Air.hpp:137
The base class for all fluid systems.
Definition: BaseFluidSystem.hpp:44
Scalar Scalar
The type used for scalar quantities.
Definition: BaseFluidSystem.hpp:49
static Evaluation henry(const Evaluation &temperature)
Henry coefficent for air in liquid water.
Definition: H2O_Air.hpp:55
static Evaluation liquidDiffCoeff(const Evaluation &temperature, const Evaluation &)
Lacking better data on water-air diffusion in liquids, we use at the moment the diffusion coefficient...
Definition: H2O_Air.hpp:101
static Evaluation gasDiffCoeff(const Evaluation &temperature, const Evaluation &pressure)
Binary diffusion coefficent for molecular water and air.
Definition: H2O_Air.hpp:70
static Scalar acentricFactor()
Returns the acentric factor of the component.
Definition: Component.hpp:110
static void init(Scalar, Scalar, unsigned, Scalar, Scalar, unsigned)
A default routine for initialization, not needed for components and must not be called.
Definition: Component.hpp:60
A fluid system with a liquid and a gaseous phase and water and air as components.
Definition: H2OAirFluidSystem.hpp:60
H2Otype H2O
The type of the water component used for this fluid system.
Definition: H2OAirFluidSystem.hpp:71
static LhsEval density(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &, unsigned phaseIdx)
Calculate the density [kg/m^3] of a fluid phase.
Definition: H2OAirFluidSystem.hpp:254
static const int AirIdx
The index of the air component.
Definition: H2OAirFluidSystem.hpp:141
static const int numComponents
Number of chemical species in the fluid system.
Definition: H2OAirFluidSystem.hpp:136
static Scalar molarMass(unsigned compIdx)
Return the molar mass of a component in [kg/mol].
Definition: H2OAirFluidSystem.hpp:155
static const int gasPhaseIdx
The index of the gas phase.
Definition: H2OAirFluidSystem.hpp:81
static bool isIdealGas(unsigned phaseIdx)
Returns true if and only if a fluid phase is assumed to be an ideal gas.
Definition: H2OAirFluidSystem.hpp:113
static const int liquidPhaseIdx
The index of the liquid phase.
Definition: H2OAirFluidSystem.hpp:79
static LhsEval binaryDiffusionCoefficient(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &, unsigned phaseIdx, unsigned)
Calculate the binary molecular diffusion coefficient for a component in a fluid phase [mol^2 * s / (k...
Definition: H2OAirFluidSystem.hpp:380
static const int H2OIdx
The index of the water component.
Definition: H2OAirFluidSystem.hpp:139
static bool isCompressible(unsigned phaseIdx)
Returns true if and only if a fluid phase is assumed to be compressible.
Definition: H2OAirFluidSystem.hpp:101
static const int numPhases
Number of fluid phases in the fluid system.
Definition: H2OAirFluidSystem.hpp:76
static LhsEval viscosity(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &, unsigned phaseIdx)
Calculate the dynamic viscosity of a fluid phase [Pa*s].
Definition: H2OAirFluidSystem.hpp:304
static bool isIdealMixture(unsigned)
Returns true if and only if a fluid phase is assumed to be an ideal mixture.
Definition: H2OAirFluidSystem.hpp:122
static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp, Scalar pressMin, Scalar pressMax, unsigned nPress)
Initialize the fluid system's static parameters using problem specific temperature and pressure range...
Definition: H2OAirFluidSystem.hpp:243
::Opm::Air< Scalar > Air
The type of the air component used for this fluid system.
Definition: H2OAirFluidSystem.hpp:73
static const char * phaseName(unsigned phaseIdx)
Return the human readable name of a fluid phase.
Definition: H2OAirFluidSystem.hpp:84
static Scalar criticalTemperature(unsigned compIdx)
Critical temperature of a component [K].
Definition: H2OAirFluidSystem.hpp:171
static bool isLiquid(unsigned phaseIdx)
Return whether a phase is liquid.
Definition: H2OAirFluidSystem.hpp:94
static LhsEval enthalpy(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &, unsigned phaseIdx)
Given a phase's composition, temperature, pressure and density, calculate its specific enthalpy [J/kg...
Definition: H2OAirFluidSystem.hpp:397
static Scalar acentricFactor(unsigned compIdx)
The acentric factor of a component [].
Definition: H2OAirFluidSystem.hpp:201
static void init()
Initialize the fluid system's static parameters.
Definition: H2OAirFluidSystem.hpp:221
static Scalar criticalPressure(unsigned compIdx)
Critical pressure of a component [Pa].
Definition: H2OAirFluidSystem.hpp:186
static LhsEval thermalConductivity(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &, unsigned phaseIdx)
Thermal conductivity of a fluid phase [W/(m K)].
Definition: H2OAirFluidSystem.hpp:429
static LhsEval fugacityCoefficient(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &, unsigned phaseIdx, unsigned compIdx)
Calculate the fugacity coefficient [Pa] of an individual component in a fluid phase.
Definition: H2OAirFluidSystem.hpp:356
static const char * componentName(unsigned compIdx)
Return the human readable name of a component.
Definition: H2OAirFluidSystem.hpp:144
static Evaluation liquidDensity(const Evaluation &temperature, const Evaluation &pressure, bool extrapolate=false)
The density of pure water in at a given pressure and temperature.
Definition: H2O.hpp:698
static const Scalar criticalTemperature()
Returns the critical temperature of water.
Definition: H2O.hpp:92
static Evaluation gasDensity(const Evaluation &temperature, const Evaluation &pressure)
The density of steam in at a given pressure and temperature.
Definition: H2O.hpp:572
static const char * name()
A human readable name for the water.
Definition: H2O.hpp:74
static Evaluation vaporPressure(Evaluation temperature)
The vapor pressure in of pure water at a given temperature.
Definition: H2O.hpp:138
static Evaluation gasViscosity(const Evaluation &temperature, const Evaluation &pressure)
The dynamic viscosity of steam.
Definition: H2O.hpp:802
static Evaluation liquidEnthalpy(const Evaluation &temperature, const Evaluation &pressure)
Specific enthalpy of liquid water .
Definition: H2O.hpp:236
static Evaluation liquidThermalConductivity(const Evaluation &temperature, const Evaluation &pressure)
Thermal conductivity of water (IAPWS) .
Definition: H2O.hpp:858
static const Scalar acentricFactor()
The acentric factor of water.
Definition: H2O.hpp:86
static bool gasIsIdeal()
Returns true iff the gas phase is assumed to be ideal.
Definition: H2O.hpp:638
static const Scalar criticalPressure()
Returns the critical pressure of water.
Definition: H2O.hpp:98
static const Scalar molarMass()
The molar mass in of water.
Definition: H2O.hpp:80
static bool liquidIsCompressible()
Returns true iff the liquid phase is assumed to be compressible.
Definition: H2O.hpp:556
static Evaluation liquidViscosity(const Evaluation &temperature, const Evaluation &pressure, bool extrapolate=false)
The dynamic viscosity of pure water.
Definition: H2O.hpp:828
static Evaluation gasThermalConductivity(const Evaluation &temperature, const Evaluation &pressure)
Thermal conductivity of water (IAPWS) .
Definition: H2O.hpp:878
static Evaluation gasEnthalpy(const Evaluation &temperature, const Evaluation &pressure)
Specific enthalpy of water steam .
Definition: H2O.hpp:183
Relations valid for an ideal gas.
Definition: IdealGas.hpp:38
A parameter cache which does nothing.
Definition: NullParameterCache.hpp:40
Definition: H2OAirFluidSystem.hpp:68