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_REGULARIZED_BROOKS_COREY_PARAMS_HPP
00028 #define OPM_REGULARIZED_BROOKS_COREY_PARAMS_HPP
00029
00030 #include "BrooksCorey.hpp"
00031 #include "BrooksCoreyParams.hpp"
00032
00033 #include <cassert>
00034
00035 #include <opm/material/common/EnsureFinalized.hpp>
00036
00037 namespace Opm {
00044 template <class TraitsT>
00045 class RegularizedBrooksCoreyParams : public Opm::BrooksCoreyParams<TraitsT>
00046 {
00047 typedef Opm::BrooksCoreyParams<TraitsT> BrooksCoreyParams;
00048 typedef Opm::BrooksCorey<TraitsT, RegularizedBrooksCoreyParams> BrooksCorey;
00049 typedef typename TraitsT::Scalar Scalar;
00050
00051 public:
00052 typedef TraitsT Traits;
00053
00054 RegularizedBrooksCoreyParams()
00055 : BrooksCoreyParams()
00056 , pcnwLowSw_(0.01)
00057 {
00058 }
00059
00060 RegularizedBrooksCoreyParams(Scalar entryPressure, Scalar lambda)
00061 : BrooksCoreyParams(entryPressure, lambda)
00062 , pcnwLowSw_(0.01)
00063 { finalize(); }
00064
00069 void finalize()
00070 {
00071 BrooksCoreyParams::finalize();
00072
00073 pcnwLow_ = BrooksCorey::twoPhaseSatPcnw(*this, pcnwLowSw_);
00074 pcnwSlopeLow_ = dPcnw_dSw_(pcnwLowSw_);
00075 pcnwHigh_ = BrooksCorey::twoPhaseSatPcnw(*this, 1.0);
00076 pcnwSlopeHigh_ = dPcnw_dSw_(1.0);
00077 }
00078
00083 Scalar pcnwLowSw() const
00084 { EnsureFinalized::check(); return pcnwLowSw_; }
00085
00090 Scalar pcnwLow() const
00091 { EnsureFinalized::check(); return pcnwLow_; }
00092
00099 Scalar pcnwSlopeLow() const
00100 { EnsureFinalized::check(); return pcnwSlopeLow_; }
00101
00106 void setPcLowSw(Scalar value)
00107 { pcnwLowSw_ = value; }
00108
00113 Scalar pcnwHigh() const
00114 { EnsureFinalized::check(); return pcnwHigh_; }
00115
00122 Scalar pcnwSlopeHigh() const
00123 { EnsureFinalized::check(); return pcnwSlopeHigh_; }
00124
00125 private:
00126 Scalar dPcnw_dSw_(Scalar Sw) const
00127 {
00128
00129
00130 const Scalar eps = 1e-7;
00131 Scalar delta = 0.0;
00132 Scalar pc1;
00133 Scalar pc2;
00134 if (Sw + eps < 1.0) {
00135 pc2 = BrooksCorey::twoPhaseSatPcnw(*this, Sw + eps);
00136 delta += eps;
00137 }
00138 else
00139 pc2 = BrooksCorey::twoPhaseSatPcnw(*this, Sw);
00140
00141 if (Sw - eps > 0.0) {
00142 pc1 = BrooksCorey::twoPhaseSatPcnw(*this, Sw - eps);
00143 delta += eps;
00144 }
00145 else
00146 pc1 = BrooksCorey::twoPhaseSatPcnw(*this, Sw);
00147
00148 return (pc2 - pc1)/delta;
00149 }
00150
00151 Scalar pcnwLowSw_;
00152 Scalar pcnwLow_;
00153 Scalar pcnwSlopeLow_;
00154 Scalar pcnwHigh_;
00155 Scalar pcnwSlopeHigh_;
00156 };
00157 }
00158
00159 #endif