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_SOMERTON_HPP
00028 #define OPM_SOMERTON_HPP
00029
00030 #include "SomertonParams.hpp"
00031
00032 #include <opm/material/common/Spline.hpp>
00033
00034 #include <opm/common/Valgrind.hpp>
00035 #include <opm/material/common/MathToolbox.hpp>
00036
00037 #include <algorithm>
00038
00039 namespace Opm
00040 {
00060 template <class FluidSystem,
00061 class ScalarT,
00062 class ParamsT = SomertonParams<FluidSystem::numPhases, ScalarT> >
00063 class Somerton
00064 {
00065 enum { numPhases = FluidSystem::numPhases };
00066
00067 public:
00068 typedef ParamsT Params;
00069 typedef typename Params::Scalar Scalar;
00070
00089 template <class FluidState, class Evaluation = Scalar>
00090 static Evaluation heatConductivity(const Params& params,
00091 const FluidState& fluidState)
00092 {
00093 Valgrind::CheckDefined(params.vacuumLambda());
00094
00095 Evaluation lambda = 0;
00096 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00097 Valgrind::CheckDefined(params.fullySaturatedLambda(phaseIdx));
00098
00099 if (FluidSystem::isLiquid(phaseIdx)) {
00100 const auto& sat = Opm::decay<Evaluation>(fluidState.saturation(phaseIdx));
00101 lambda +=
00102 regularizedSqrt_(Opm::max(0.0, Opm::min(1.0, sat)))
00103 * (params.fullySaturatedLambda(phaseIdx) - params.vacuumLambda());
00104 }
00105 else {
00106 lambda += params.fullySaturatedLambda(phaseIdx) - params.vacuumLambda();
00107 }
00108 };
00109
00110 lambda += params.vacuumLambda();
00111 assert(lambda >= 0);
00112 return lambda;
00113 }
00114
00115 protected:
00116 template <class Evaluation>
00117 static Evaluation regularizedSqrt_(const Evaluation& x)
00118 {
00119 typedef Opm::Spline<Scalar> Spline;
00120
00121 static const Scalar xMin = 1e-2;
00122 static const Scalar sqrtXMin = std::sqrt(xMin);
00123 static const Scalar fPrimeXMin = 1.0/(2*std::sqrt(xMin));
00124 static const Scalar fPrime0 = 2*fPrimeXMin;
00125 static const Spline sqrtRegSpline(0, xMin,
00126 0, sqrtXMin,
00127 fPrime0, fPrimeXMin);
00128
00129 if (x > xMin)
00130 return Opm::sqrt(x);
00131 else if (x <= 0)
00132 return fPrime0 * x;
00133 else
00134 return sqrtRegSpline.eval(x);
00135 }
00136 };
00137 }
00138
00139 #endif