00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef OPM_H2O_AIR_SYSTEM_HPP
00028 #define OPM_H2O_AIR_SYSTEM_HPP
00029
00030 #include "BaseFluidSystem.hpp"
00031 #include "NullParameterCache.hpp"
00032
00033 #include <opm/material/IdealGas.hpp>
00034 #include <opm/material/binarycoefficients/H2O_Air.hpp>
00035 #include <opm/material/components/Air.hpp>
00036 #include <opm/material/components/H2O.hpp>
00037 #include <opm/material/components/TabulatedComponent.hpp>
00038 #include <opm/common/Valgrind.hpp>
00039
00040 #include <opm/common/Exceptions.hpp>
00041 #include <opm/common/ErrorMacros.hpp>
00042
00043 #include <iostream>
00044 #include <cassert>
00045
00046 namespace Opm {
00047 namespace FluidSystems {
00048
00058 template <class Scalar,
00059
00060 class H2Otype = Opm::TabulatedComponent<Scalar, Opm::H2O<Scalar> >,
00061 bool useComplexRelations = true>
00062 class H2OAir
00063 : public BaseFluidSystem<Scalar, H2OAir<Scalar, H2Otype, useComplexRelations> >
00064 {
00065 typedef H2OAir<Scalar,H2Otype, useComplexRelations > ThisType;
00066 typedef BaseFluidSystem <Scalar, ThisType> Base;
00067 typedef Opm::IdealGas<Scalar> IdealGas;
00068
00069 public:
00070 template <class Evaluation>
00071 struct ParameterCache : public Opm::NullParameterCache<Evaluation>
00072 {};
00073
00075 typedef H2Otype H2O;
00077 typedef Opm::Air<Scalar> Air;
00078
00080 static const int numPhases = 2;
00081
00083 static const int liquidPhaseIdx = 0;
00085 static const int gasPhaseIdx = 1;
00086
00088 static const char* phaseName(unsigned phaseIdx)
00089 {
00090 switch (phaseIdx) {
00091 case liquidPhaseIdx: return "liquid";
00092 case gasPhaseIdx: return "gas";
00093 };
00094 OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
00095 }
00096
00098 static bool isLiquid(unsigned phaseIdx)
00099 {
00100
00101 return phaseIdx != gasPhaseIdx;
00102 }
00103
00105 static bool isCompressible(unsigned phaseIdx)
00106 {
00107
00108 return (phaseIdx == gasPhaseIdx)
00109
00110 ? true
00111 :
00112
00113 H2O::liquidIsCompressible();
00114 }
00115
00117 static bool isIdealGas(unsigned phaseIdx)
00118 {
00119 return
00120 (phaseIdx == gasPhaseIdx)
00121 ? H2O::gasIsIdeal() && Air::gasIsIdeal()
00122 : false;
00123 }
00124
00126 static bool isIdealMixture(unsigned )
00127 {
00128
00129
00130
00131
00132 return true;
00133 }
00134
00135
00136
00137
00138
00140 static const int numComponents = 2;
00141
00143 static const int H2OIdx = 0;
00145 static const int AirIdx = 1;
00146
00148 static const char* componentName(unsigned compIdx)
00149 {
00150 switch (compIdx)
00151 {
00152 case H2OIdx: return H2O::name();
00153 case AirIdx: return Air::name();
00154 };
00155 OPM_THROW(std::logic_error, "Invalid component index " << compIdx);
00156 }
00157
00159 static Scalar molarMass(unsigned compIdx)
00160 {
00161 return
00162 (compIdx == H2OIdx)
00163 ? H2O::molarMass()
00164 : (compIdx == AirIdx)
00165 ? Air::molarMass()
00166 : 1e30;
00167 }
00168
00169
00175 static Scalar criticalTemperature(unsigned compIdx)
00176 {
00177 return
00178 (compIdx == H2OIdx)
00179 ? H2O::criticalTemperature()
00180 : (compIdx == AirIdx)
00181 ? Air::criticalTemperature()
00182 : 1e30;
00183 }
00184
00190 static Scalar criticalPressure(unsigned compIdx)
00191 {
00192 return
00193 (compIdx == H2OIdx)
00194 ? H2O::criticalPressure()
00195 : (compIdx == AirIdx)
00196 ? Air::criticalPressure()
00197 : 1e30;
00198 }
00199
00205 static Scalar acentricFactor(unsigned compIdx)
00206 {
00207 return
00208 (compIdx == H2OIdx)
00209 ? H2O::acentricFactor()
00210 : (compIdx == AirIdx)
00211 ? Air::acentricFactor()
00212 : 1e30;
00213 }
00214
00215
00216
00217
00218
00225 static void init()
00226 {
00227 if (H2O::isTabulated)
00228 init(273.15,
00229 623.15,
00230 50,
00231 -10,
00232 20e6,
00233 50);
00234 }
00235
00247 static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
00248 Scalar pressMin, Scalar pressMax, unsigned nPress)
00249 {
00250 if (H2O::isTabulated) {
00251 H2O::init(tempMin, tempMax, nTemp,
00252 pressMin, pressMax, nPress);
00253 }
00254 }
00255
00257 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00258 static LhsEval density(const FluidState& fluidState,
00259 const ParameterCache<ParamCacheEval>& ,
00260 unsigned phaseIdx)
00261 {
00262 assert(0 <= phaseIdx && phaseIdx < numPhases);
00263
00264 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00265 LhsEval p;
00266 if (isCompressible(phaseIdx))
00267 p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00268 else {
00269
00270
00271 p = - 1e30;
00272 Valgrind::SetUndefined(p);
00273 }
00274
00275
00276 LhsEval sumMoleFrac = 0;
00277 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
00278 sumMoleFrac += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
00279
00280 if (phaseIdx == liquidPhaseIdx)
00281 {
00282 if (!useComplexRelations)
00283
00284 return H2O::liquidDensity(T, p);
00285 else
00286 {
00287
00288 const LhsEval& rholH2O = H2O::liquidDensity(T, p);
00289 const LhsEval& clH2O = rholH2O/H2O::molarMass();
00290
00291 const auto& xlH2O = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
00292 const auto& xlAir = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, AirIdx));
00293
00294 return clH2O*(H2O::molarMass()*xlH2O + Air::molarMass()*xlAir)/sumMoleFrac;
00295 }
00296 }
00297 else if (phaseIdx == gasPhaseIdx)
00298 {
00299 if (!useComplexRelations)
00300
00301 return
00302 IdealGas::molarDensity(T, p)
00303 * Opm::decay<LhsEval>(fluidState.averageMolarMass(gasPhaseIdx))
00304 / Opm::max(1e-5, sumMoleFrac);
00305
00306 LhsEval partialPressureH2O =
00307 Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))
00308 *Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
00309
00310 LhsEval partialPressureAir =
00311 Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, AirIdx))
00312 *Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
00313
00314 return H2O::gasDensity(T, partialPressureH2O) + Air::gasDensity(T, partialPressureAir);
00315 }
00316 OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
00317 }
00318
00320 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00321 static LhsEval viscosity(const FluidState& fluidState,
00322 const ParameterCache<ParamCacheEval>& ,
00323 unsigned phaseIdx)
00324 {
00325 assert(0 <= phaseIdx && phaseIdx < numPhases);
00326
00327 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00328 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00329
00330 if (phaseIdx == liquidPhaseIdx)
00331 {
00332
00333
00334
00335 return H2O::liquidViscosity(T, p);
00336 }
00337 else if (phaseIdx == gasPhaseIdx)
00338 {
00339 if(!useComplexRelations){
00340 return Air::gasViscosity(T, p);
00341 }
00342 else
00343 {
00344
00345
00346
00347
00348
00349
00350
00351
00352 LhsEval muResult = 0;
00353 const LhsEval mu[numComponents] = {
00354 H2O::gasViscosity(T, H2O::vaporPressure(T)),
00355 Air::gasViscosity(T, p)
00356 };
00357
00358
00359 const Scalar M[numComponents] = {
00360 H2O::molarMass(),
00361 Air::molarMass()
00362 };
00363
00364 for (unsigned i = 0; i < numComponents; ++i) {
00365 LhsEval divisor = 0;
00366 for (unsigned j = 0; j < numComponents; ++j) {
00367 LhsEval phiIJ =
00368 1 +
00369 Opm::sqrt(mu[i]/mu[j]) *
00370 std::pow(M[j]/M[i], 1./4.0);
00371
00372 phiIJ *= phiIJ;
00373 phiIJ /= std::sqrt(8*(1 + M[i]/M[j]));
00374 divisor += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, j))*phiIJ;
00375 }
00376 const auto& xAlphaI = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, i));
00377 muResult += xAlphaI*mu[i]/divisor;
00378 }
00379 return muResult;
00380 }
00381 }
00382 OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
00383 }
00384
00386 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00387 static LhsEval fugacityCoefficient(const FluidState& fluidState,
00388 const ParameterCache<ParamCacheEval>& ,
00389 unsigned phaseIdx,
00390 unsigned compIdx)
00391 {
00392 assert(0 <= phaseIdx && phaseIdx < numPhases);
00393 assert(0 <= compIdx && compIdx < numComponents);
00394
00395 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00396 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00397
00398 if (phaseIdx == liquidPhaseIdx) {
00399 if (compIdx == H2OIdx)
00400 return H2O::vaporPressure(T)/p;
00401 return Opm::BinaryCoeff::H2O_Air::henry(T)/p;
00402 }
00403
00404
00405
00406 return 1.0;
00407 }
00408
00410 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00411 static LhsEval binaryDiffusionCoefficient(const FluidState& fluidState,
00412 const ParameterCache<ParamCacheEval>& ,
00413 unsigned phaseIdx,
00414 unsigned )
00415 {
00416 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00417 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00418
00419 if (phaseIdx == liquidPhaseIdx)
00420 return BinaryCoeff::H2O_Air::liquidDiffCoeff(T, p);
00421
00422 assert(phaseIdx == gasPhaseIdx);
00423 return BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);
00424 }
00425
00427 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00428 static LhsEval enthalpy(const FluidState& fluidState,
00429 const ParameterCache<ParamCacheEval>& ,
00430 unsigned phaseIdx)
00431 {
00432 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00433 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00434 Valgrind::CheckDefined(T);
00435 Valgrind::CheckDefined(p);
00436
00437 if (phaseIdx == liquidPhaseIdx)
00438 {
00439
00440 return H2O::liquidEnthalpy(T, p);
00441 }
00442
00443 else if (phaseIdx == gasPhaseIdx)
00444 {
00445 LhsEval result = 0.0;
00446 result +=
00447 H2O::gasEnthalpy(T, p) *
00448 Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
00449
00450 result +=
00451 Air::gasEnthalpy(T, p) *
00452 Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, AirIdx));
00453 return result;
00454 }
00455 OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
00456 }
00457
00459 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00460 static LhsEval thermalConductivity(const FluidState& fluidState,
00461 const ParameterCache<ParamCacheEval>& ,
00462 unsigned phaseIdx)
00463 {
00464 assert(0 <= phaseIdx && phaseIdx < numPhases);
00465
00466 const LhsEval& temperature =
00467 Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00468 const LhsEval& pressure =
00469 Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00470
00471 if (phaseIdx == liquidPhaseIdx)
00472 return H2O::liquidThermalConductivity(temperature, pressure);
00473 else {
00474 const LhsEval& lambdaDryAir = Air::gasThermalConductivity(temperature, pressure);
00475
00476 if (useComplexRelations){
00477 const LhsEval& xAir =
00478 Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, AirIdx));
00479 const LhsEval& xH2O =
00480 Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
00481 LhsEval lambdaAir = xAir*lambdaDryAir;
00482
00483
00484
00485 LhsEval partialPressure = pressure*xH2O;
00486
00487 LhsEval lambdaH2O =
00488 xH2O*H2O::gasThermalConductivity(temperature, partialPressure);
00489 return lambdaAir + lambdaH2O;
00490 }
00491 else
00492 return lambdaDryAir;
00493 }
00494 }
00495 };
00496
00497 }
00498
00499 }
00500
00501 #endif