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_PENG_ROBINSON_MIXTURE_HPP
00028 #define OPM_PENG_ROBINSON_MIXTURE_HPP
00029
00030 #include "PengRobinson.hpp"
00031
00032 #include <opm/material/Constants.hpp>
00033
00034 #include <iostream>
00035
00036 namespace Opm {
00041 template <class Scalar, class StaticParameters>
00042 class PengRobinsonMixture
00043 {
00044 enum { numComponents = StaticParameters::numComponents };
00045 typedef Opm::PengRobinson<Scalar> PengRobinson;
00046
00047
00048 PengRobinsonMixture() {}
00049
00050
00051 static const Scalar R;
00052
00053
00054 static const Scalar u;
00055 static const Scalar w;
00056
00057 public:
00064 template <class MutableParams, class FluidState>
00065 static int computeMolarVolumes(Scalar* Vm,
00066 const MutableParams& params,
00067 unsigned phaseIdx,
00068 const FluidState& fs)
00069 {
00070 return PengRobinson::computeMolarVolumes(Vm, params, phaseIdx, fs);
00071 }
00072
00090 template <class FluidState, class Params, class LhsEval = typename FluidState::Scalar>
00091 static LhsEval computeFugacityCoefficient(const FluidState& fs,
00092 const Params& params,
00093 unsigned phaseIdx,
00094 unsigned compIdx)
00095 {
00096
00097
00098
00099 LhsEval Vm = params.molarVolume(phaseIdx);
00100
00101
00102 LhsEval bi_b = params.bPure(phaseIdx, compIdx) / params.b(phaseIdx);
00103
00104
00105 LhsEval RT = R*fs.temperature(phaseIdx);
00106 LhsEval p = fs.pressure(phaseIdx);
00107 LhsEval Z = p*Vm/RT;
00108
00109
00110 LhsEval Astar = params.a(phaseIdx)*p/(RT*RT);
00111 LhsEval Bstar = params.b(phaseIdx)*p/(RT);
00112
00113
00114 LhsEval sumMoleFractions = 0.0;
00115 for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx)
00116 sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx);
00117 LhsEval deltai = 2*Opm::sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
00118 LhsEval tmp = 0;
00119 for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
00120 tmp +=
00121 fs.moleFraction(phaseIdx, compJIdx)
00122 / sumMoleFractions
00123 * Opm::sqrt(params.aPure(phaseIdx, compJIdx))
00124 * (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx));
00125 };
00126 deltai *= tmp;
00127
00128 LhsEval base =
00129 (2*Z + Bstar*(u + std::sqrt(u*u - 4*w))) /
00130 (2*Z + Bstar*(u - std::sqrt(u*u - 4*w)));
00131 LhsEval expo = Astar/(Bstar*std::sqrt(u*u - 4*w))*(bi_b - deltai);
00132
00133 LhsEval fugCoeff =
00134 Opm::exp(bi_b*(Z - 1))/Opm::max(1e-9, Z - Bstar) *
00135 Opm::pow(base, expo);
00136
00138
00139
00140
00141
00142
00143 fugCoeff = Opm::min(1e10, fugCoeff);
00144
00145
00146
00147
00148 fugCoeff = Opm::max(1e-10, fugCoeff);
00150
00151 return fugCoeff;
00152 }
00153
00154 };
00155
00156 template <class Scalar, class StaticParameters>
00157 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::R = Opm::Constants<Scalar>::R;
00158 template<class Scalar, class StaticParameters>
00159 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::u = 2.0;
00160 template<class Scalar, class StaticParameters>
00161 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::w = -1.0;
00162
00163 }
00164
00165 #endif