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_MESITYLENE_FLUID_SYSTEM_HPP
00028 #define OPM_H2O_AIR_MESITYLENE_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/Air.hpp>
00036 #include <opm/material/components/H2O.hpp>
00037 #include <opm/material/components/SimpleH2O.hpp>
00038 #include <opm/material/components/Mesitylene.hpp>
00039 #include <opm/material/components/TabulatedComponent.hpp>
00040 #include <opm/material/binarycoefficients/H2O_Air.hpp>
00041 #include <opm/material/binarycoefficients/H2O_N2.hpp>
00042 #include <opm/material/binarycoefficients/H2O_Mesitylene.hpp>
00043 #include <opm/material/binarycoefficients/Air_Mesitylene.hpp>
00044
00045 #include <iostream>
00046
00047 namespace Opm {
00048 namespace FluidSystems {
00049
00055 template <class Scalar>
00056 class H2OAirMesitylene
00057 : public BaseFluidSystem<Scalar, H2OAirMesitylene<Scalar> >
00058 {
00059 typedef H2OAirMesitylene<Scalar> ThisType;
00060 typedef BaseFluidSystem<Scalar, ThisType> Base;
00061
00062 typedef Opm::H2O<Scalar> IapwsH2O;
00063 typedef Opm::TabulatedComponent<Scalar, IapwsH2O, false> TabulatedH2O;
00064
00065 public:
00066 template <class Evaluation>
00067 struct ParameterCache : public Opm::NullParameterCache<Evaluation>
00068 {};
00069
00071 typedef Opm::Mesitylene<Scalar> NAPL;
00072
00074 typedef Opm::Air<Scalar> Air;
00075
00077
00078 typedef TabulatedH2O H2O;
00079
00080
00082 static const int numPhases = 3;
00084 static const int numComponents = 3;
00085
00087 static const int waterPhaseIdx = 0;
00089 static const int naplPhaseIdx = 1;
00091 static const int gasPhaseIdx = 2;
00092
00094 static const int H2OIdx = 0;
00096 static const int NAPLIdx = 1;
00098 static const int airIdx = 2;
00099
00101 static void init()
00102 {
00103 init(273.15,
00104 623.15,
00105 50,
00106 0.0,
00107 20e6,
00108 50);
00109 }
00110
00122 static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
00123 Scalar pressMin, Scalar pressMax, unsigned nPress)
00124 {
00125 if (H2O::isTabulated) {
00126 TabulatedH2O::init(tempMin, tempMax, nTemp,
00127 pressMin, pressMax, nPress);
00128 }
00129 }
00130
00132 static bool isLiquid(unsigned phaseIdx)
00133 {
00134
00135 return phaseIdx != gasPhaseIdx;
00136 }
00137
00139 static bool isIdealGas(unsigned phaseIdx)
00140 { return phaseIdx == gasPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
00141
00143 static bool isCompressible(unsigned phaseIdx)
00144 {
00145
00146
00147 return (phaseIdx == gasPhaseIdx)
00148 ? true
00149 : (phaseIdx == waterPhaseIdx)
00150 ? H2O::liquidIsCompressible()
00151 : NAPL::liquidIsCompressible();
00152 }
00153
00155 static bool isIdealMixture(unsigned )
00156 {
00157
00158
00159
00160
00161 return true;
00162 }
00163
00165 static const char* phaseName(unsigned phaseIdx)
00166 {
00167 switch (phaseIdx) {
00168 case waterPhaseIdx: return "water";
00169 case naplPhaseIdx: return "napl";
00170 case gasPhaseIdx: return "gas";
00171 };
00172 OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
00173 }
00174
00176 static const char* componentName(unsigned compIdx)
00177 {
00178 switch (compIdx) {
00179 case H2OIdx: return H2O::name();
00180 case airIdx: return Air::name();
00181 case NAPLIdx: return NAPL::name();
00182 };
00183 OPM_THROW(std::logic_error, "Invalid component index " << compIdx);
00184 }
00185
00187 static Scalar molarMass(unsigned compIdx)
00188 {
00189 return
00190 (compIdx == H2OIdx)
00191 ? H2O::molarMass()
00192 : (compIdx == airIdx)
00193 ? Air::molarMass()
00194 : (compIdx == NAPLIdx)
00195 ? NAPL::molarMass()
00196 : -1e10;
00197
00198 }
00199
00201 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00202 static LhsEval density(const FluidState& fluidState,
00203 const ParameterCache<ParamCacheEval>& ,
00204 unsigned phaseIdx)
00205 {
00206 const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00207
00208 if (phaseIdx == waterPhaseIdx) {
00209
00210 const LhsEval& p =
00211 H2O::liquidIsCompressible()
00212 ? Opm::decay<LhsEval>(fluidState.pressure(phaseIdx))
00213 : 1e30;
00214
00215 const LhsEval& rholH2O = H2O::liquidDensity(T, p);
00216 const LhsEval& clH2O = rholH2O/H2O::molarMass();
00217
00218
00219
00220 return
00221 clH2O*(H2O::molarMass()*Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx)) +
00222 Air::molarMass()*Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx)) +
00223 NAPL::molarMass()*Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx)));
00224 }
00225 else if (phaseIdx == naplPhaseIdx) {
00226
00227 const LhsEval& p =
00228 NAPL::liquidIsCompressible()
00229 ? Opm::decay<LhsEval>(fluidState.pressure(phaseIdx))
00230 : 1e30;
00231 return NAPL::liquidDensity(T, p);
00232 }
00233
00234 assert (phaseIdx == gasPhaseIdx);
00235 const LhsEval& pg = Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
00236 const LhsEval& pH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*pg;
00237 const LhsEval& pAir = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*pg;
00238 const LhsEval& pNAPL = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*pg;
00239 return
00240 H2O::gasDensity(T, pH2O) +
00241 Air::gasDensity(T, pAir) +
00242 NAPL::gasDensity(T, pNAPL);
00243 }
00244
00246 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00247 static LhsEval viscosity(const FluidState& fluidState,
00248 const ParameterCache<ParamCacheEval>& ,
00249 unsigned phaseIdx)
00250 {
00251 const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00252 const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00253
00254 if (phaseIdx == waterPhaseIdx) {
00255
00256
00257 return H2O::liquidViscosity(T,
00258 p);
00259 }
00260 else if (phaseIdx == naplPhaseIdx) {
00261
00262 return NAPL::liquidViscosity(T, p);
00263 }
00264
00265 assert (phaseIdx == gasPhaseIdx);
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278 const LhsEval mu[numComponents] = {
00279 H2O::gasViscosity(T, H2O::vaporPressure(T)),
00280 Air::gasViscosity(T, p),
00281 NAPL::gasViscosity(T, NAPL::vaporPressure(T))
00282 };
00283
00284 const Scalar M[numComponents] = {
00285 H2O::molarMass(),
00286 Air::molarMass(),
00287 NAPL::molarMass()
00288 };
00289
00290 const LhsEval& xgAir = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
00291 const LhsEval& xgH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
00292 const LhsEval& xgNapl = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
00293 const LhsEval& xgAW = xgAir + xgH2O;
00294 const LhsEval& muAW = (mu[airIdx]*xgAir + mu[H2OIdx]*xgH2O)/xgAW;
00295 const LhsEval& MAW = (xgAir*Air::molarMass() + xgH2O*H2O::molarMass())/xgAW;
00296
00297 Scalar phiCAW = 0.3;
00298
00299
00300
00301
00302 const LhsEval& phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW);
00303
00304 return (xgAW*muAW)/(xgAW + xgNapl*phiAWC) + (xgNapl*mu[NAPLIdx])/(xgNapl + xgAW*phiCAW);
00305 }
00306
00308 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00309 static LhsEval diffusionCoefficient(const FluidState& ,
00310 const ParameterCache<ParamCacheEval>& ,
00311 unsigned ,
00312 unsigned )
00313 {
00314 return 0;
00315 #if 0
00316 typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
00317
00318 const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00319 const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00320 LhsEval diffCont;
00321
00322 if (phaseIdx==gasPhaseIdx) {
00323 const LhsEval& diffAC = Opm::BinaryCoeff::Air_Mesitylene::gasDiffCoeff(T, p);
00324 const LhsEval& diffWC = Opm::BinaryCoeff::H2O_Mesitylene::gasDiffCoeff(T, p);
00325 const LhsEval& diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);
00326
00327 const LhsEval& xga = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
00328 const LhsEval& xgw = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
00329 const LhsEval& xgc = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
00330
00331 if (compIdx==NAPLIdx) return (1 - xgw)/(xga/diffAW + xgc/diffWC);
00332 else if (compIdx==H2OIdx) return (1 - xgc)/(xgw/diffWC + xga/diffAC);
00333 else if (compIdx==airIdx) OPM_THROW(std::logic_error,
00334 "Diffusivity of air in the gas phase "
00335 "is constraint by sum of diffusive fluxes = 0 !\n");
00336 }
00337 else if (phaseIdx==waterPhaseIdx){
00338 const LhsEval& diffACl = 1.e-9;
00339 const LhsEval& diffWCl = 1.e-9;
00340 const LhsEval& diffAWl = 1.e-9;
00341
00342 const LhsEval& xwa = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
00343 const LhsEval& xww = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
00344 const LhsEval& xwc = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
00345
00346 switch (compIdx) {
00347 case NAPLIdx:
00348 diffCont = (1.- xww)/(xwa/diffAWl + xwc/diffWCl);
00349 return diffCont;
00350 case airIdx:
00351 diffCont = (1.- xwc)/(xww/diffWCl + xwa/diffACl);
00352 return diffCont;
00353 case H2OIdx:
00354 OPM_THROW(std::logic_error,
00355 "Diffusivity of water in the water phase "
00356 "is constraint by sum of diffusive fluxes = 0 !\n");
00357 };
00358 }
00359 else if (phaseIdx==naplPhaseIdx) {
00360 OPM_THROW(std::logic_error,
00361 "Diffusion coefficients of "
00362 "substances in liquid phase are undefined!\n");
00363 }
00364 return 0;
00365 #endif
00366 }
00367
00369 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00370 static LhsEval fugacityCoefficient(const FluidState& fluidState,
00371 const ParameterCache<ParamCacheEval>& ,
00372 unsigned phaseIdx,
00373 unsigned compIdx)
00374 {
00375 assert(0 <= phaseIdx && phaseIdx < numPhases);
00376 assert(0 <= compIdx && compIdx < numComponents);
00377
00378 const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00379 const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00380 Valgrind::CheckDefined(T);
00381 Valgrind::CheckDefined(p);
00382
00383 if (phaseIdx == waterPhaseIdx) {
00384 if (compIdx == H2OIdx)
00385 return H2O::vaporPressure(T)/p;
00386 else if (compIdx == airIdx)
00387 return Opm::BinaryCoeff::H2O_N2::henry(T)/p;
00388 else if (compIdx == NAPLIdx)
00389 return Opm::BinaryCoeff::H2O_Mesitylene::henry(T)/p;
00390 assert(false);
00391 }
00392
00393
00394
00395
00396
00397 else if (phaseIdx == naplPhaseIdx) {
00398 const LhsEval& phiNapl = NAPL::vaporPressure(T)/p;
00399 if (compIdx == NAPLIdx)
00400 return phiNapl;
00401 else if (compIdx == airIdx)
00402 return 1e6*phiNapl;
00403 else if (compIdx == H2OIdx)
00404 return 1e6*phiNapl;
00405 assert(false);
00406 }
00407
00408
00409
00410 assert(phaseIdx == gasPhaseIdx);
00411 return 1.0;
00412 }
00413
00414
00416 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00417 static LhsEval enthalpy(const FluidState& fluidState,
00418 const ParameterCache<ParamCacheEval>& ,
00419 unsigned phaseIdx)
00420 {
00421 const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00422 const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00423
00424 if (phaseIdx == waterPhaseIdx) {
00425 return H2O::liquidEnthalpy(T, p);
00426 }
00427 else if (phaseIdx == naplPhaseIdx) {
00428 return NAPL::liquidEnthalpy(T, p);
00429 }
00430 else if (phaseIdx == gasPhaseIdx) {
00431
00432 LhsEval result = 0;
00433 result += H2O::gasEnthalpy(T, p) * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
00434 result += NAPL::gasEnthalpy(T, p) * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
00435 result += Air::gasEnthalpy(T, p) * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));
00436
00437 return result;
00438 }
00439 OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
00440 }
00441
00443 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00444 static LhsEval thermalConductivity(const FluidState& fluidState,
00445 const ParameterCache<ParamCacheEval>& ,
00446 unsigned phaseIdx)
00447 {
00448 assert(0 <= phaseIdx && phaseIdx < numPhases);
00449
00450 if (phaseIdx == waterPhaseIdx){
00451 const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00452 const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00453
00454 return H2O::liquidThermalConductivity(T, p);
00455 }
00456 else if (phaseIdx == gasPhaseIdx) {
00457 const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00458 const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00459
00460 return Air::gasThermalConductivity(T, p);
00461 }
00462
00463 assert(phaseIdx == naplPhaseIdx);
00464
00465
00466
00467
00468
00469
00470
00471
00472 return 0.0143964;
00473 }
00474 };
00475 }
00476 }
00477
00478 #endif