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_XYLENE_FLUID_SYSTEM_HPP
00028 #define OPM_H2O_AIR_XYLENE_FLUID_SYSTEM_HPP
00029
00030 #include <opm/material/IdealGas.hpp>
00031 #include <opm/material/components/Air.hpp>
00032 #include <opm/material/components/H2O.hpp>
00033 #include <opm/material/components/Xylene.hpp>
00034 #include <opm/material/components/TabulatedComponent.hpp>
00035
00036 #include <opm/material/binarycoefficients/H2O_Air.hpp>
00037 #include <opm/material/binarycoefficients/H2O_Xylene.hpp>
00038 #include <opm/material/binarycoefficients/Air_Xylene.hpp>
00039
00040 #include "BaseFluidSystem.hpp"
00041 #include "NullParameterCache.hpp"
00042
00043 namespace Opm {
00044 namespace FluidSystems {
00045
00051 template <class Scalar>
00052 class H2OAirXylene
00053 : public BaseFluidSystem<Scalar, H2OAirXylene<Scalar> >
00054 {
00055 typedef H2OAirXylene<Scalar> ThisType;
00056 typedef BaseFluidSystem<Scalar, ThisType> Base;
00057
00058 public:
00059 template <class Evaluation>
00060 struct ParameterCache : public Opm::NullParameterCache<Evaluation>
00061 {};
00062
00064 typedef Opm::H2O<Scalar> H2O;
00066 typedef Opm::Xylene<Scalar> NAPL;
00068 typedef Opm::Air<Scalar> Air;
00069
00071 static const int numPhases = 3;
00073 static const int numComponents = 3;
00074
00076 static const int waterPhaseIdx = 0;
00078 static const int naplPhaseIdx = 1;
00080 static const int gasPhaseIdx = 2;
00081
00083 static const int H2OIdx = 0;
00085 static const int NAPLIdx = 1;
00087 static const int airIdx = 2;
00088
00090 static void init()
00091 { }
00092
00094 static bool isLiquid(unsigned phaseIdx)
00095 {
00096
00097 return phaseIdx != gasPhaseIdx;
00098 }
00099
00101 static bool isIdealGas(unsigned phaseIdx)
00102 { return phaseIdx == gasPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
00103
00105 static bool isIdealMixture(unsigned )
00106 {
00107
00108
00109
00110
00111
00112 return true;
00113 }
00114
00116 static bool isCompressible(unsigned phaseIdx)
00117 {
00118 return
00119 (phaseIdx == gasPhaseIdx)
00120
00121 ? true
00122 : (phaseIdx == waterPhaseIdx)
00123
00124 ? H2O::liquidIsCompressible()
00125
00126 : NAPL::liquidIsCompressible();
00127 }
00128
00130 static const char* phaseName(unsigned phaseIdx)
00131 {
00132 switch (phaseIdx) {
00133 case waterPhaseIdx: return "water";
00134 case naplPhaseIdx: return "napl";
00135 case gasPhaseIdx: return "gas";
00136 };
00137 OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
00138 }
00139
00141 static const char* componentName(unsigned compIdx)
00142 {
00143 switch (compIdx) {
00144 case H2OIdx: return H2O::name();
00145 case airIdx: return Air::name();
00146 case NAPLIdx: return NAPL::name();
00147 };
00148 OPM_THROW(std::logic_error, "Invalid component index " << compIdx);
00149 }
00150
00152 static Scalar molarMass(unsigned compIdx)
00153 {
00154 return
00155 (compIdx == H2OIdx)
00156
00157 ? H2O::molarMass()
00158 : (compIdx == airIdx)
00159
00160 ? Air::molarMass()
00161
00162 : (compIdx == NAPLIdx)
00163 ? NAPL::molarMass()
00164 : 1e30;
00165 }
00166
00168 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00169 static LhsEval density(const FluidState& fluidState,
00170 const ParameterCache<ParamCacheEval>& ,
00171 unsigned phaseIdx)
00172 {
00173 if (phaseIdx == waterPhaseIdx) {
00174 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00175 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00176
00177
00178
00179 const LhsEval& rholH2O = H2O::liquidDensity(T, p);
00180 const LhsEval& clH2O = rholH2O/H2O::molarMass();
00181
00182 const auto& xwH2O = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
00183 const auto& xwAir = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
00184 const auto& xwNapl = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
00185
00186
00187 return clH2O*(H2O::molarMass()*xwH2O + Air::molarMass()*xwAir + NAPL::molarMass()*xwNapl);
00188 }
00189 else if (phaseIdx == naplPhaseIdx) {
00190
00191 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00192 return NAPL::liquidDensity(T, LhsEval(1e30));
00193 }
00194
00195 assert (phaseIdx == gasPhaseIdx);
00196 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00197 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00198
00199 const LhsEval& pH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*p;
00200 const LhsEval& pAir = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*p;
00201 const LhsEval& pNAPL = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*p;
00202 return
00203 H2O::gasDensity(T, pH2O) +
00204 Air::gasDensity(T, pAir) +
00205 NAPL::gasDensity(T, pNAPL);
00206 }
00207
00209 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00210 static LhsEval viscosity(const FluidState& fluidState,
00211 const ParameterCache<ParamCacheEval>& ,
00212 unsigned phaseIdx)
00213 {
00214 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00215 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00216
00217 if (phaseIdx == waterPhaseIdx) {
00218
00219 return H2O::liquidViscosity(T, p);
00220 }
00221 else if (phaseIdx == naplPhaseIdx) {
00222
00223 return NAPL::liquidViscosity(T, p);
00224 }
00225
00226 assert (phaseIdx == gasPhaseIdx);
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239 const LhsEval mu[numComponents] = {
00240 H2O::gasViscosity(T, H2O::vaporPressure(T)),
00241 Air::simpleGasViscosity(T, p),
00242 NAPL::gasViscosity(T, NAPL::vaporPressure(T))
00243 };
00244
00245 const Scalar M[numComponents] = {
00246 H2O::molarMass(),
00247 Air::molarMass(),
00248 NAPL::molarMass()
00249 };
00250
00251 const auto& xgAir = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
00252 const auto& xgH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
00253 const auto& xgNapl = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
00254
00255 const LhsEval& xgAW = xgAir + xgH2O;
00256 const LhsEval& muAW = (mu[airIdx]*xgAir + mu[H2OIdx]*xgH2O)/ xgAW;
00257
00258 const LhsEval& MAW = (xgAir*Air::molarMass() + xgH2O*H2O::molarMass())/xgAW;
00259
00260 Scalar phiCAW = 0.3;
00261
00262
00263
00264
00265 const LhsEval& phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW);
00266
00267 return (xgAW*muAW)/(xgAW+xgNapl*phiAWC) + (xgNapl*mu[NAPLIdx])/(xgNapl + xgAW*phiCAW);
00268 }
00269
00271 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00272 static LhsEval diffusionCoefficient(const FluidState& fluidState,
00273 const ParameterCache<ParamCacheEval>& ,
00274 unsigned phaseIdx,
00275 unsigned compIdx)
00276 {
00277 if (phaseIdx==gasPhaseIdx) {
00278 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00279 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00280
00281 const LhsEval& diffAC = Opm::BinaryCoeff::Air_Xylene::gasDiffCoeff(T, p);
00282 const LhsEval& diffWC = Opm::BinaryCoeff::H2O_Xylene::gasDiffCoeff(T, p);
00283 const LhsEval& diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);
00284
00285 const LhsEval& xga = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
00286 const LhsEval& xgw = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
00287 const LhsEval& xgc = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
00288
00289 if (compIdx==NAPLIdx) return (1.- xgw)/(xga/diffAW + xgc/diffWC);
00290 else if (compIdx==H2OIdx) return (1.- xgc)/(xgw/diffWC + xga/diffAC);
00291 else if (compIdx==airIdx) OPM_THROW(std::logic_error,
00292 "Diffusivity of air in the gas phase "
00293 "is constraint by sum of diffusive fluxes = 0 !\n");
00294 } else if (phaseIdx==waterPhaseIdx){
00295 Scalar diffACl = 1.e-9;
00296 Scalar diffWCl = 1.e-9;
00297 Scalar diffAWl = 1.e-9;
00298
00299 const LhsEval& xwa = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
00300 const LhsEval& xww = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
00301 const LhsEval& xwc = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
00302
00303 switch (compIdx) {
00304 case NAPLIdx:
00305 return (1.- xww)/(xwa/diffAWl + xwc/diffWCl);
00306 case airIdx:
00307 return (1.- xwc)/(xww/diffWCl + xwa/diffACl);
00308 case H2OIdx:
00309 OPM_THROW(std::logic_error,
00310 "Diffusivity of water in the water phase "
00311 "is constraint by sum of diffusive fluxes = 0 !\n");
00312 };
00313 } else if (phaseIdx==naplPhaseIdx) {
00314
00315 OPM_THROW(std::logic_error,
00316 "Diffusion coefficients of "
00317 "substances in liquid phase are undefined!\n");
00318 }
00319 return 0;
00320 }
00321
00323 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00324 static LhsEval fugacityCoefficient(const FluidState& fluidState,
00325 const ParameterCache<ParamCacheEval>& ,
00326 unsigned phaseIdx,
00327 unsigned compIdx)
00328 {
00329 assert(0 <= phaseIdx && phaseIdx < numPhases);
00330 assert(0 <= compIdx && compIdx < numComponents);
00331
00332 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00333 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00334
00335 if (phaseIdx == waterPhaseIdx) {
00336 if (compIdx == H2OIdx)
00337 return H2O::vaporPressure(T)/p;
00338 else if (compIdx == airIdx)
00339 return Opm::BinaryCoeff::H2O_Air::henry(T)/p;
00340 else if (compIdx == NAPLIdx)
00341 return Opm::BinaryCoeff::H2O_Xylene::henry(T)/p;
00342 }
00343
00344
00345
00346
00347
00348
00349 if (phaseIdx == naplPhaseIdx) {
00350 const LhsEval& phiNapl = NAPL::vaporPressure(T)/p;
00351 if (compIdx == NAPLIdx)
00352 return phiNapl;
00353 else if (compIdx == airIdx)
00354 return 1e6*phiNapl;
00355 else if (compIdx == H2OIdx)
00356 return 1e6*phiNapl;
00357 }
00358
00359
00360
00361 assert(phaseIdx == gasPhaseIdx);
00362 return 1.0;
00363 }
00364
00366 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00367 static LhsEval enthalpy(const FluidState& fluidState,
00368 const ParameterCache<ParamCacheEval>& ,
00369 unsigned phaseIdx)
00370 {
00371 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00372 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00373
00374 if (phaseIdx == waterPhaseIdx) {
00375 return H2O::liquidEnthalpy(T, p);
00376 }
00377 else if (phaseIdx == naplPhaseIdx) {
00378 return NAPL::liquidEnthalpy(T, p);
00379 }
00380 else if (phaseIdx == gasPhaseIdx) {
00381 const LhsEval& hgc = NAPL::gasEnthalpy(T, p);
00382 const LhsEval& hgw = H2O::gasEnthalpy(T, p);
00383 const LhsEval& hga = Air::gasEnthalpy(T, p);
00384
00385 LhsEval result = 0;
00386 result += hgw * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
00387 result += hga * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
00388 result += hgc * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));
00389
00390 return result;
00391 }
00392 OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
00393 }
00394
00395 private:
00396 template <class LhsEval>
00397 static LhsEval waterPhaseDensity_(const LhsEval& T,
00398 const LhsEval& pw,
00399 const LhsEval& xww,
00400 const LhsEval& xwa,
00401 const LhsEval& xwc)
00402 {
00403 const LhsEval& rholH2O = H2O::liquidDensity(T, pw);
00404 const LhsEval& clH2O = rholH2O/H2O::molarMass();
00405
00406
00407
00408 return clH2O*(xww*H2O::molarMass() + xwa*Air::molarMass() + xwc*NAPL::molarMass());
00409 }
00410
00411 template <class LhsEval>
00412 static LhsEval gasPhaseDensity_(const LhsEval& T,
00413 const LhsEval& pg,
00414 const LhsEval& xgw,
00415 const LhsEval& xga,
00416 const LhsEval& xgc)
00417 { return H2O::gasDensity(T, pg*xgw) + Air::gasDensity(T, pg*xga) + NAPL::gasDensity(T, pg*xgc); }
00418
00419 template <class LhsEval>
00420 static LhsEval NAPLPhaseDensity_(const LhsEval& T, const LhsEval& pn)
00421 { return NAPL::liquidDensity(T, pn); }
00422 };
00423 }
00424 }
00425
00426 #endif