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_N2_FLUID_SYSTEM_HPP
00028 #define OPM_H2O_N2_FLUID_SYSTEM_HPP
00029
00030 #include "BaseFluidSystem.hpp"
00031 #include "NullParameterCache.hpp"
00032
00033 #include <opm/material/IdealGas.hpp>
00034 #include <opm/material/components/N2.hpp>
00035 #include <opm/material/components/H2O.hpp>
00036 #include <opm/material/components/SimpleH2O.hpp>
00037 #include <opm/material/components/TabulatedComponent.hpp>
00038 #include <opm/material/binarycoefficients/H2O_N2.hpp>
00039 #include <opm/common/Valgrind.hpp>
00040
00041 #include <opm/common/Exceptions.hpp>
00042 #include <opm/common/ErrorMacros.hpp>
00043
00044 #include <iostream>
00045 #include <cassert>
00046
00047 namespace Opm {
00048 namespace FluidSystems {
00049
00055 template <class Scalar, bool useComplexRelations = true>
00056 class H2ON2
00057 : public BaseFluidSystem<Scalar, H2ON2<Scalar, useComplexRelations> >
00058 {
00059 typedef H2ON2<Scalar, useComplexRelations> ThisType;
00060 typedef BaseFluidSystem<Scalar, ThisType> Base;
00061
00062
00063 typedef Opm::IdealGas<Scalar> IdealGas;
00064 typedef Opm::H2O<Scalar> IapwsH2O;
00065 typedef Opm::TabulatedComponent<Scalar, IapwsH2O > TabulatedH2O;
00066 typedef Opm::N2<Scalar> SimpleN2;
00067
00068 public:
00070 template <class Evaluation>
00071 using ParameterCache = NullParameterCache<Evaluation>;
00072
00073
00074
00075
00076
00078 static const int numPhases = 2;
00079
00081 static const int liquidPhaseIdx = 0;
00083 static const int gasPhaseIdx = 1;
00084
00086 static const char* phaseName(unsigned phaseIdx)
00087 {
00088 static const char* name[] = {
00089 "liquid",
00090 "gas"
00091 };
00092
00093 assert(0 <= phaseIdx && phaseIdx < numPhases);
00094 return name[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
00109 return
00110 (phaseIdx == gasPhaseIdx)
00111 ? true
00112 :H2O::liquidIsCompressible();
00113 }
00114
00116 static bool isIdealGas(unsigned phaseIdx)
00117 {
00118
00119
00120 return
00121 (phaseIdx == gasPhaseIdx)
00122 ? H2O::gasIsIdeal() && N2::gasIsIdeal()
00123 : false;
00124 }
00125
00127 static bool isIdealMixture(unsigned )
00128 {
00129
00130
00131
00132
00133 return true;
00134 }
00135
00136
00137
00138
00139
00141 static const int numComponents = 2;
00142
00144 static const int H2OIdx = 0;
00146 static const int N2Idx = 1;
00147
00149 typedef TabulatedH2O H2O;
00150
00151
00152
00154 typedef SimpleN2 N2;
00155
00157 static const char* componentName(unsigned compIdx)
00158 {
00159 static const char* name[] = {
00160 H2O::name(),
00161 N2::name()
00162 };
00163
00164 assert(0 <= compIdx && compIdx < numComponents);
00165 return name[compIdx];
00166 }
00167
00169 static Scalar molarMass(unsigned compIdx)
00170 {
00171
00172 return (compIdx == H2OIdx)
00173 ? H2O::molarMass()
00174 : (compIdx == N2Idx)
00175 ? N2::molarMass()
00176 : 1e30;
00177 }
00178
00184 static Scalar criticalTemperature(unsigned compIdx)
00185 {
00186 return (compIdx == H2OIdx)
00187 ? H2O::criticalTemperature()
00188 : (compIdx == N2Idx)
00189 ? N2::criticalTemperature()
00190 : 1e30;
00191 }
00192
00198 static Scalar criticalPressure(unsigned compIdx)
00199 {
00200 return (compIdx == H2OIdx)
00201 ? H2O::criticalPressure()
00202 : (compIdx == N2Idx)
00203 ? N2::criticalPressure()
00204 : 1e30;
00205 }
00206
00212 static Scalar acentricFactor(unsigned compIdx)
00213 {
00214 return (compIdx == H2OIdx)
00215 ? H2O::acentricFactor()
00216 : (compIdx == N2Idx)
00217 ? N2::acentricFactor()
00218 : 1e30;
00219 }
00220
00221
00222
00223
00224
00231 static void init()
00232 {
00233 init(273.15,
00234 623.15,
00235 50,
00236 0.0,
00237 20e6,
00238 50);
00239 }
00240
00252 static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
00253 Scalar pressMin, Scalar pressMax, unsigned nPress)
00254 {
00255 if (H2O::isTabulated) {
00256 TabulatedH2O::init(tempMin, tempMax, nTemp,
00257 pressMin, pressMax, nPress);
00258 }
00259 }
00260
00268 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00269 static LhsEval density(const FluidState& fluidState,
00270 const ParameterCache<ParamCacheEval>& ,
00271 unsigned phaseIdx)
00272 {
00273 assert(0 <= phaseIdx && phaseIdx < numPhases);
00274
00275 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00276 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00277
00278 LhsEval sumMoleFrac = 0;
00279 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
00280 sumMoleFrac += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
00281
00282
00283 if (phaseIdx == liquidPhaseIdx) {
00284 if (!useComplexRelations)
00285
00286 return H2O::liquidDensity(T, p);
00287 else
00288 {
00289
00290 const auto& rholH2O = H2O::liquidDensity(T, p);
00291 const auto& clH2O = rholH2O/H2O::molarMass();
00292
00293 const auto& xlH2O = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
00294 const auto& xlN2 = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, N2Idx));
00295
00296
00297
00298 return clH2O*(H2O::molarMass()*xlH2O + N2::molarMass()*xlN2)/sumMoleFrac;
00299 }
00300 }
00301
00302
00303 assert(phaseIdx == gasPhaseIdx);
00304
00305 if (!useComplexRelations)
00306
00307 return
00308 IdealGas::molarDensity(T, p)
00309 * Opm::decay<LhsEval>(fluidState.averageMolarMass(gasPhaseIdx))
00310 / Opm::max(1e-5, sumMoleFrac);
00311
00312
00313
00314 const auto& xgH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
00315 const auto& xgN2 = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, N2Idx));
00316 const auto& rho_gH2O = H2O::gasDensity(T, p*xgH2O);
00317 const auto& rho_gN2 = N2::gasDensity(T, p*xgN2);
00318 return (rho_gH2O + rho_gN2)/Opm::max(1e-5, sumMoleFrac);
00319 }
00320
00322 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00323 static LhsEval viscosity(const FluidState& fluidState,
00324 const ParameterCache<ParamCacheEval>& ,
00325 unsigned phaseIdx)
00326 {
00327 assert(0 <= phaseIdx && phaseIdx < numPhases);
00328
00329 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00330 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00331
00332
00333 if (phaseIdx == liquidPhaseIdx)
00334
00335 return H2O::liquidViscosity(T, p);
00336
00337
00338 assert(phaseIdx == gasPhaseIdx);
00339
00340 if (!useComplexRelations)
00341
00342 return N2::gasViscosity(T, p);
00343 else {
00344
00345
00346
00347
00348
00349
00350 LhsEval muResult = 0;
00351 const LhsEval mu[numComponents] = {
00352 H2O::gasViscosity(T, H2O::vaporPressure(T)),
00353 N2::gasViscosity(T, p)
00354 };
00355
00356 LhsEval sumx = 0.0;
00357 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
00358 sumx += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
00359 sumx = Opm::max(1e-10, sumx);
00360
00361 for (unsigned i = 0; i < numComponents; ++i) {
00362 LhsEval divisor = 0;
00363 for (unsigned j = 0; j < numComponents; ++j) {
00364 LhsEval phiIJ = 1 + Opm::sqrt(mu[i]/mu[j]) * std::pow(molarMass(j)/molarMass(i), 1/4.0);
00365 phiIJ *= phiIJ;
00366 phiIJ /= std::sqrt(8*(1 + molarMass(i)/molarMass(j)));
00367 divisor +=
00368 Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, j))
00369 /sumx*phiIJ;
00370 }
00371 muResult +=
00372 Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, i))
00373 /sumx*mu[i]/divisor;
00374 }
00375 return muResult;
00376 }
00377 }
00378
00380 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00381 static LhsEval fugacityCoefficient(const FluidState& fluidState,
00382 const ParameterCache<ParamCacheEval>& ,
00383 unsigned phaseIdx,
00384 unsigned compIdx)
00385 {
00386 assert(0 <= phaseIdx && phaseIdx < numPhases);
00387 assert(0 <= compIdx && compIdx < numComponents);
00388
00389 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00390 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00391
00392
00393 if (phaseIdx == liquidPhaseIdx) {
00394 if (compIdx == H2OIdx)
00395 return H2O::vaporPressure(T)/p;
00396 return Opm::BinaryCoeff::H2O_N2::henry(T)/p;
00397 }
00398
00399 assert(phaseIdx == gasPhaseIdx);
00400
00401
00402
00403 return 1.0;
00404 }
00405
00407 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00408 static LhsEval diffusionCoefficient(const FluidState& fluidState,
00409 const ParameterCache<ParamCacheEval>& ,
00410 unsigned phaseIdx,
00411 unsigned )
00412
00413 {
00414 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00415 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00416
00417
00418 if (phaseIdx == liquidPhaseIdx)
00419 return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
00420
00421
00422 assert(phaseIdx == gasPhaseIdx);
00423 return BinaryCoeff::H2O_N2::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
00438 if (phaseIdx == liquidPhaseIdx) {
00439
00440 return H2O::liquidEnthalpy(T, p);
00441 }
00442
00443
00444 assert(phaseIdx == gasPhaseIdx);
00445
00446
00447
00448
00449
00450 const auto& XgH2O = Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
00451 const auto& XgN2 = Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, N2Idx));
00452
00453 LhsEval hH2O = XgH2O*H2O::gasEnthalpy(T, p);
00454 LhsEval hN2 = XgN2*N2::gasEnthalpy(T, p);
00455 return hH2O + hN2;
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 auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00467 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00468 if (phaseIdx == liquidPhaseIdx)
00469 return H2O::liquidThermalConductivity(T, p);
00470
00471
00472 assert(phaseIdx == gasPhaseIdx);
00473
00474 if (useComplexRelations){
00475
00476 const auto& xH2O = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
00477 const auto& xN2 = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, N2Idx));
00478
00479
00480
00481 const auto& lambdaN2 = N2::gasThermalConductivity(T, p*xN2);
00482 const auto& lambdaH2O = H2O::gasThermalConductivity(T, p*xH2O);
00483
00484 return lambdaN2 + lambdaH2O;
00485 }
00486 else
00487
00488 return N2::gasThermalConductivity(T, p);
00489 }
00490
00492 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00493 static LhsEval heatCapacity(const FluidState& fluidState,
00494 const ParameterCache<ParamCacheEval>& ,
00495 unsigned phaseIdx)
00496 {
00497 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00498 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00499 const auto& xAlphaH2O = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
00500 const auto& xAlphaN2 = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, N2Idx));
00501 const auto& XAlphaH2O = Opm::decay<LhsEval>(fluidState.massFraction(phaseIdx, H2OIdx));
00502 const auto& XAlphaN2 = Opm::decay<LhsEval>(fluidState.massFraction(phaseIdx, N2Idx));
00503
00504 if (phaseIdx == liquidPhaseIdx)
00505 return H2O::liquidHeatCapacity(T, p);
00506
00507 assert(phaseIdx == gasPhaseIdx);
00508
00509
00510
00511
00512 LhsEval c_pN2;
00513 LhsEval c_pH2O;
00514
00515 if (useComplexRelations) {
00516 c_pN2 = N2::gasHeatCapacity(T, p*xAlphaN2);
00517 c_pH2O = H2O::gasHeatCapacity(T, p*xAlphaH2O);
00518 }
00519 else {
00520
00521
00522
00523 Scalar c_vN2molar = Opm::Constants<Scalar>::R*2.39;
00524 Scalar c_pN2molar = Opm::Constants<Scalar>::R + c_vN2molar;
00525
00526 Scalar c_vH2Omolar = Opm::Constants<Scalar>::R*3.37;
00527 Scalar c_pH2Omolar = Opm::Constants<Scalar>::R + c_vH2Omolar;
00528
00529 c_pN2 = c_pN2molar/molarMass(N2Idx);
00530 c_pH2O = c_pH2Omolar/molarMass(H2OIdx);
00531 }
00532
00533
00534
00535 return XAlphaH2O*c_pH2O + XAlphaN2*c_pN2;
00536 }
00537 };
00538
00539 }
00540
00541 }
00542
00543 #endif