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_PARAMS_MIXTURE_HPP
00028 #define OPM_PENG_ROBINSON_PARAMS_MIXTURE_HPP
00029
00030 #include "PengRobinsonParams.hpp"
00031
00032 #include <opm/material/common/MathToolbox.hpp>
00033 #include <opm/material/Constants.hpp>
00034
00035 #include <algorithm>
00036
00037 namespace Opm
00038 {
00039
00057 template <class Scalar, class FluidSystem, unsigned phaseIdx, bool useSpe5Relations=false>
00058 class PengRobinsonParamsMixture
00059 : public PengRobinsonParams<Scalar>
00060 {
00061 enum { numComponents = FluidSystem::numComponents };
00062
00063
00064 typedef Opm::PengRobinsonParams<Scalar> PureParams;
00065
00066 typedef MathToolbox<Scalar> Toolbox;
00067
00068
00069 static const Scalar R;
00070
00071 public:
00075 template <class FluidState>
00076 void updatePure(const FluidState& fluidState)
00077 {
00078 updatePure(fluidState.temperature(phaseIdx),
00079 fluidState.pressure(phaseIdx));
00080 }
00081
00087 void updatePure(Scalar temperature, Scalar pressure)
00088 {
00089 Valgrind::CheckDefined(temperature);
00090 Valgrind::CheckDefined(pressure);
00091
00092
00093
00094
00095
00096
00097 for (unsigned i = 0; i < numComponents; ++i) {
00098 Scalar pc = FluidSystem::criticalPressure(i);
00099 Scalar omega = FluidSystem::acentricFactor(i);
00100 Scalar Tr = temperature/FluidSystem::criticalTemperature(i);
00101 Scalar RTc = R*FluidSystem::criticalTemperature(i);
00102
00103 Scalar f_omega;
00104
00105 if (useSpe5Relations) {
00106 if (omega < 0.49) f_omega = 0.37464 + omega*(1.54226 + omega*(-0.26992));
00107 else f_omega = 0.379642 + omega*(1.48503 + omega*(-0.164423 + omega*0.016666));
00108 }
00109 else
00110 f_omega = 0.37464 + omega*(1.54226 - omega*0.26992);
00111
00112 Valgrind::CheckDefined(f_omega);
00113
00114 Scalar tmp = 1 + f_omega*(1 - Opm::sqrt(Tr));
00115 tmp = tmp*tmp;
00116
00117 Scalar newA = 0.4572355*RTc*RTc/pc * tmp;
00118 Scalar newB = 0.0777961 * RTc / pc;
00119 assert(std::isfinite(Opm::scalarValue(newA)));
00120 assert(std::isfinite(Opm::scalarValue(newB)));
00121
00122 this->pureParams_[i].setA(newA);
00123 this->pureParams_[i].setB(newB);
00124 Valgrind::CheckDefined(this->pureParams_[i].a());
00125 Valgrind::CheckDefined(this->pureParams_[i].b());
00126 }
00127
00128 updateACache_();
00129 }
00130
00138 template <class FluidState>
00139 void updateMix(const FluidState& fs)
00140 {
00141 Scalar sumx = 0.0;
00142 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
00143 sumx += fs.moleFraction(phaseIdx, compIdx);
00144 sumx = std::max(Scalar(1e-10), sumx);
00145
00146
00147
00148
00149
00150 Scalar newA = 0;
00151 Scalar newB = 0;
00152 for (unsigned compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
00153 const Scalar moleFracI = fs.moleFraction(phaseIdx, compIIdx);
00154 Scalar xi = Opm::max(0.0, Opm::min(1.0, moleFracI));
00155 Valgrind::CheckDefined(xi);
00156
00157 for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
00158 const Scalar moleFracJ = fs.moleFraction(phaseIdx, compJIdx );
00159 Scalar xj = Opm::max(0.0, Opm::min(1.0, moleFracJ));
00160 Valgrind::CheckDefined(xj);
00161
00162
00163 newA += xi * xj * aCache_[compIIdx][compJIdx];
00164
00165 assert(std::isfinite(Opm::scalarValue(newA)));
00166 }
00167
00168
00169 newB += Opm::max(0.0, xi) * this->pureParams_[compIIdx].b();
00170 assert(std::isfinite(Opm::scalarValue(newB)));
00171 }
00172
00173
00174 this->setA(newA);
00175 this->setB(newB);
00176
00177 Valgrind::CheckDefined(this->a());
00178 Valgrind::CheckDefined(this->b());
00179
00180 }
00181
00190 template <class FluidState>
00191 void updateSingleMoleFraction(const FluidState& fs,
00192 unsigned )
00193 {
00194 updateMix(fs);
00195 }
00196
00200 const PureParams& pureParams(unsigned compIdx) const
00201 { return pureParams_[compIdx]; }
00202
00206 const PureParams& operator[](unsigned compIdx) const
00207 {
00208 assert(0 <= compIdx && compIdx < numComponents);
00209 return pureParams_[compIdx];
00210 }
00211
00216 void checkDefined() const
00217 {
00218 #ifndef NDEBUG
00219 for (unsigned i = 0; i < numComponents; ++i)
00220 pureParams_[i].checkDefined();
00221
00222 Valgrind::CheckDefined(this->a());
00223 Valgrind::CheckDefined(this->b());
00224 #endif
00225 }
00226
00227 protected:
00228 PureParams pureParams_[numComponents];
00229
00230 private:
00231 void updateACache_()
00232 {
00233 for (unsigned compIIdx = 0; compIIdx < numComponents; ++ compIIdx) {
00234 for (unsigned compJIdx = 0; compJIdx < numComponents; ++ compJIdx) {
00235
00236 Scalar Psi = FluidSystem::interactionCoefficient(compIIdx, compJIdx);
00237
00238 aCache_[compIIdx][compJIdx] =
00239 Opm::sqrt(this->pureParams_[compIIdx].a()
00240 * this->pureParams_[compJIdx].a())
00241 * (1 - Psi);
00242 }
00243 }
00244 }
00245
00246 Scalar aCache_[numComponents][numComponents];
00247 };
00248
00249 template <class Scalar, class FluidSystem, unsigned phaseIdx, bool useSpe5Relations>
00250 const Scalar PengRobinsonParamsMixture<Scalar, FluidSystem, phaseIdx, useSpe5Relations>::R = Opm::Constants<Scalar>::R;
00251
00252 }
00253
00254 #endif