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_HPP
00028 #define OPM_H2O_HPP
00029
00030 #include <cmath>
00031 #include <cassert>
00032
00033 #include <opm/material/IdealGas.hpp>
00034 #include <opm/common/Exceptions.hpp>
00035 #include <opm/common/ErrorMacros.hpp>
00036 #include <opm/common/Valgrind.hpp>
00037
00038 #include "Component.hpp"
00039
00040 #include "iapws/Common.hpp"
00041 #include "iapws/Region1.hpp"
00042 #include "iapws/Region2.hpp"
00043 #include "iapws/Region4.hpp"
00044
00045 namespace Opm {
00046
00060 template <class Scalar>
00061 class H2O : public Component<Scalar, H2O<Scalar> >
00062 {
00063 typedef IAPWS::Common<Scalar> Common;
00064 typedef IAPWS::Region1<Scalar> Region1;
00065 typedef IAPWS::Region2<Scalar> Region2;
00066 typedef IAPWS::Region4<Scalar> Region4;
00067
00068 static const Scalar Rs;
00069
00070 public:
00074 static const char* name()
00075 { return "H2O"; }
00076
00080 static const Scalar molarMass()
00081 { return Common::molarMass; }
00082
00086 static const Scalar acentricFactor()
00087 { return Common::acentricFactor; }
00088
00092 static const Scalar criticalTemperature()
00093 { return Common::criticalTemperature; }
00094
00098 static const Scalar criticalPressure()
00099 { return Common::criticalPressure; }
00100
00104 static const Scalar criticalMolarVolume()
00105 { return Common::criticalMolarVolume; }
00106
00110 static const Scalar tripleTemperature()
00111 { return Common::tripleTemperature; }
00112
00116 static const Scalar triplePressure()
00117 { return Common::triplePressure; }
00118
00131 template <class Evaluation>
00132 static Evaluation vaporPressure(Evaluation temperature)
00133 {
00134 if (temperature > criticalTemperature())
00135 temperature = criticalTemperature();
00136 if (temperature < tripleTemperature())
00137 temperature = tripleTemperature();
00138
00139 return Region4::saturationPressure(temperature);
00140 }
00153 template <class Evaluation>
00154 static Evaluation vaporTemperature(const Evaluation& pressure)
00155 {
00156 if (pressure > criticalPressure())
00157 pressure = criticalPressure();
00158 if (pressure < triplePressure())
00159 pressure = triplePressure();
00160
00161 return Region4::vaporTemperature(pressure);
00162 }
00163
00176 template <class Evaluation>
00177 static Evaluation gasEnthalpy(const Evaluation& temperature,
00178 const Evaluation& pressure)
00179 {
00180 if (!Region2::isValid(temperature, pressure))
00181 {
00182 OPM_THROW(NumericalProblem,
00183 "Enthalpy of steam is only implemented for temperatures below 623.15K and "
00184 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00185 }
00186
00187
00188 if (pressure < triplePressure() - 100) {
00189
00190
00191
00192
00193
00194
00195 return enthalpyRegion2_<Evaluation>(temperature, triplePressure() - 100);
00196 }
00197 Evaluation pv = vaporPressure(temperature);
00198 if (pressure > pv) {
00199
00200
00201 Evaluation dh_dp =
00202 Rs*temperature*
00203 Region2::tau(temperature)*
00204 Region2::dpi_dp(pv)*
00205 Region2::ddgamma_dtaudpi(temperature, pv);
00206
00207 return
00208 enthalpyRegion2_(temperature, pv) +
00209 (pressure - pv)*dh_dp;
00210 };
00211
00212 return enthalpyRegion2_(temperature, pressure);
00213 }
00214
00227 template <class Evaluation>
00228 static Evaluation liquidEnthalpy(const Evaluation& temperature,
00229 const Evaluation& pressure)
00230 {
00231 if (!Region1::isValid(temperature, pressure))
00232 {
00233 OPM_THROW(NumericalProblem,
00234 "Enthalpy of water is only implemented for temperatures below 623.15K and "
00235 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00236 }
00237
00238
00239 const Evaluation& pv = vaporPressure(temperature);
00240 if (pressure < pv) {
00241
00242
00243 const Evaluation& dh_dp =
00244 Rs * temperature*
00245 Region1::tau(temperature)*
00246 Region1::dpi_dp(pv)*
00247 Region1::ddgamma_dtaudpi(temperature, pv);
00248
00249 return
00250 enthalpyRegion1_(temperature, pv) +
00251 (pressure - pv)*dh_dp;
00252 };
00253
00254 return enthalpyRegion1_(temperature, pressure);
00255 }
00256
00269 template <class Evaluation>
00270 static Evaluation gasHeatCapacity(const Evaluation& temperature,
00271 const Evaluation& pressure)
00272 {
00273 if (!Region2::isValid(temperature, pressure))
00274 {
00275 OPM_THROW(NumericalProblem,
00276 "Heat capacity of steam is only implemented for temperatures below 623.15K and "
00277 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00278 }
00279
00280
00281 if (pressure < triplePressure() - 100)
00282 return heatCap_p_Region2_(temperature, Evaluation(triplePressure() - 100));
00283 const Evaluation& pv = vaporPressure(temperature);
00284 if (pressure > pv)
00285
00286
00287 return heatCap_p_Region2_(temperature, pv);
00288
00289 return heatCap_p_Region2_(temperature, pressure);
00290 }
00291
00304 template <class Evaluation>
00305 static Evaluation liquidHeatCapacity(const Evaluation& temperature,
00306 const Evaluation& pressure)
00307 {
00308 if (!Region1::isValid(temperature, pressure))
00309 {
00310 OPM_THROW(NumericalProblem,
00311 "heat Capacity of water is only implemented for temperatures below 623.15K and "
00312 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00313 }
00314
00315
00316 const Evaluation& pv = vaporPressure(temperature);
00317 if (pressure < pv) {
00318
00319
00320 return heatCap_p_Region1_(temperature, pv);
00321 };
00322
00323 return heatCap_p_Region1_(temperature, pressure);
00324 }
00325
00338 template <class Evaluation>
00339 static Evaluation liquidInternalEnergy(const Evaluation& temperature,
00340 const Evaluation& pressure)
00341 {
00342 if (!Region1::isValid(temperature, pressure))
00343 {
00344 OPM_THROW(NumericalProblem,
00345 "Internal Energy of water is only implemented for temperatures below 623.15K and "
00346 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00347 }
00348
00349
00350
00351 Scalar pv = vaporPressure<Scalar>(Opm::scalarValue(temperature));
00352 if (pressure < pv) {
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374 Scalar eps = 1e-7;
00375 const Evaluation& uv = internalEnergyRegion1_(temperature, Evaluation(pv));
00376 const Evaluation& uvPEps = internalEnergyRegion1_(temperature, Evaluation(pv + eps));
00377 const Evaluation& du_dp = (uvPEps - uv)/eps;
00378 return uv + du_dp*(pressure - pv);
00379 };
00380
00381 return internalEnergyRegion1_(temperature, pressure);
00382 }
00383
00396 template <class Evaluation>
00397 static Evaluation gasInternalEnergy(const Evaluation& temperature, const Evaluation& pressure)
00398 {
00399 if (!Region2::isValid(temperature, pressure))
00400 {
00401 OPM_THROW(NumericalProblem,
00402 "Internal Energy of steam is only implemented for temperatures below 623.15K and "
00403 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00404 }
00405
00406
00407 if (pressure < triplePressure() - 100) {
00408
00409
00410
00411
00412
00413
00414
00415
00416 return
00417 enthalpyRegion2_(temperature, Evaluation(triplePressure() - 100.0))
00418 -
00419 Rs*temperature;
00420 }
00421 Scalar pv = vaporPressure(Opm::scalarValue(temperature));
00422 if (pressure > pv) {
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448 Scalar eps = 1e-7;
00449 const Evaluation& uv = internalEnergyRegion2_(temperature, Evaluation(pv));
00450 const Evaluation& uvMEps = internalEnergyRegion2_(temperature, Evaluation(pv - eps));
00451 const Evaluation& du_dp = (uv - uvMEps)/eps;
00452 return uv + du_dp*(pressure - pv);
00453 };
00454
00455 return internalEnergyRegion2_(temperature, pressure);
00456 }
00457
00470 template <class Evaluation>
00471 static Evaluation liquidHeatCapacityConstVolume(const Evaluation& temperature,
00472 const Evaluation& pressure)
00473 {
00474 if (!Region1::isValid(temperature, pressure))
00475 {
00476 OPM_THROW(NumericalProblem,
00477 "Heat capacity of water is only implemented for temperatures below 623.15K and "
00478 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00479 }
00480
00481
00482
00483 Scalar pv = vaporPressure(temperature);
00484 if (pressure < pv) {
00485
00486
00487 return heatCap_v_Region1_(temperature, pv);
00488 }
00489
00490 return heatCap_v_Region1_(temperature, pressure);
00491 }
00492
00505 template <class Evaluation>
00506 static Evaluation gasHeatCapacityConstVolume(const Evaluation& temperature, const Evaluation& pressure)
00507 {
00508 if (!Region2::isValid(temperature, pressure))
00509 {
00510 OPM_THROW(NumericalProblem,
00511 "Heat capacity of steam is only implemented for temperatures below 623.15K and "
00512 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00513 }
00514
00515
00516 if (pressure < triplePressure() - 100) {
00517 return
00518 heatCap_v_Region2_(temperature, triplePressure() - 100);
00519 }
00520 Scalar pv = vaporPressure(temperature);
00521 if (pressure > pv) {
00522 return heatCap_v_Region2_(temperature, pv);
00523 };
00524
00525 return heatCap_v_Region2_(temperature, pressure);
00526 }
00527
00531 static bool gasIsCompressible()
00532 { return true; }
00533
00537 static bool liquidIsCompressible()
00538 { return true; }
00539
00552 template <class Evaluation>
00553 static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
00554 {
00555 if (!Region2::isValid(temperature, pressure))
00556 {
00557 OPM_THROW(NumericalProblem,
00558 "Density of steam is only implemented for temperatures below 623.15K and "
00559 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00560 }
00561
00562
00563 if (pressure < triplePressure() - 100) {
00564
00565
00566 const Evaluation& rho0IAPWS =
00567 1.0/volumeRegion2_(temperature,
00568 Evaluation(triplePressure() - 100));
00569 const Evaluation& rho0Id =
00570 IdealGas<Scalar>::density(Evaluation(molarMass()),
00571 temperature,
00572 Evaluation(triplePressure() - 100));
00573 return
00574 rho0IAPWS/rho0Id
00575 *IdealGas<Scalar>::density(Evaluation(molarMass()),
00576 temperature,
00577 pressure);
00578 }
00579 Evaluation pv = vaporPressure(temperature);
00580 if (pressure > pv) {
00581
00582
00583
00584
00585
00586
00587 Scalar eps = Opm::scalarValue(pv)*1e-8;
00588 Evaluation v0 = volumeRegion2_(temperature, pv);
00589 Evaluation v1 = volumeRegion2_(temperature, pv + eps);
00590 Evaluation dv_dp = (v1 - v0)/eps;
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606 Evaluation drho_dp = - 1/(v0*v0)*dv_dp;
00607
00608
00609 return 1.0/v0 + (pressure - pv)*drho_dp;
00610 };
00611
00612 return 1.0/volumeRegion2_(temperature, pressure);
00613 }
00614
00618 static bool gasIsIdeal()
00619 { return false; }
00620
00633 template <class Evaluation>
00634 static Evaluation gasPressure(const Evaluation& temperature, Scalar density)
00635 {
00636 Valgrind::CheckDefined(temperature);
00637 Valgrind::CheckDefined(density);
00638
00639
00640
00641 Evaluation pressure = IdealGas<Scalar>::pressure(temperature, density/molarMass());
00642 Scalar eps = pressure*1e-7;
00643
00644 Evaluation deltaP = pressure*2;
00645 Valgrind::CheckDefined(pressure);
00646 Valgrind::CheckDefined(deltaP);
00647 for (int i = 0; i < 5 && std::abs(Opm::scalarValue(pressure)*1e-9) < std::abs(Opm::scalarValue(deltaP)); ++i) {
00648 Evaluation f = gasDensity(temperature, pressure) - density;
00649
00650 Evaluation df_dp;
00651 df_dp = gasDensity(temperature, pressure + eps);
00652 df_dp -= gasDensity(temperature, pressure - eps);
00653 df_dp /= 2*eps;
00654
00655 deltaP = - f/df_dp;
00656
00657 pressure += deltaP;
00658 Valgrind::CheckDefined(pressure);
00659 Valgrind::CheckDefined(deltaP);
00660 }
00661
00662 return pressure;
00663 }
00664
00677 template <class Evaluation>
00678 static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
00679 {
00680 if (!Region1::isValid(temperature, pressure))
00681 {
00682 OPM_THROW(NumericalProblem,
00683 "Density of water is only implemented for temperatures below 623.15K and "
00684 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00685 }
00686
00687
00688 Evaluation pv = vaporPressure(temperature);
00689 if (pressure < pv) {
00690
00691
00692
00693
00694
00695 Scalar eps = Opm::scalarValue(pv)*1e-8;
00696 Evaluation v0 = volumeRegion1_(temperature, pv);
00697 Evaluation v1 = volumeRegion1_(temperature, pv + eps);
00698 Evaluation dv_dp = (v1 - v0)/eps;
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716 Evaluation drho_dp = - 1/(v0*v0)*dv_dp;
00717
00718
00719 return 1.0/v0 + (pressure - pv)*drho_dp;
00720 };
00721
00722 return 1/volumeRegion1_(temperature, pressure);
00723 }
00724
00738 template <class Evaluation>
00739 static Evaluation liquidPressure(const Evaluation& temperature, Scalar density)
00740 {
00741
00742
00743
00744 Evaluation pressure = 1.1*vaporPressure(temperature);
00745 Scalar eps = Opm::scalarValue(pressure)*1e-7;
00746
00747 Evaluation deltaP = pressure*2;
00748 for (int i = 0; i < 5 && std::abs(Opm::scalarValue(pressure)*1e-9) < std::abs(Opm::scalarValue(deltaP)); ++i) {
00749 Evaluation f = liquidDensity(temperature, pressure) - density;
00750
00751 Evaluation df_dp;
00752 df_dp = liquidDensity(temperature, pressure + eps);
00753 df_dp -= liquidDensity(temperature, pressure - eps);
00754 df_dp /= 2*eps;
00755
00756 deltaP = - f/df_dp;
00757
00758 pressure += deltaP;
00759 }
00760
00761 return pressure;
00762 }
00763
00778 template <class Evaluation>
00779 static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
00780 {
00781 if (!Region2::isValid(temperature, pressure))
00782 {
00783 OPM_THROW(NumericalProblem,
00784 "Viscosity of steam is only implemented for temperatures below 623.15K and "
00785 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00786 }
00787
00788 Evaluation rho = gasDensity(temperature, pressure);
00789 return Common::viscosity(temperature, rho);
00790 }
00791
00803 template <class Evaluation>
00804 static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
00805 {
00806 if (!Region1::isValid(temperature, pressure))
00807 {
00808 OPM_THROW(NumericalProblem,
00809 "Viscosity of water is only implemented for temperatures below 623.15K and "
00810 "pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
00811 };
00812
00813 const Evaluation& rho = liquidDensity(temperature, pressure);
00814 return Common::viscosity(temperature, rho);
00815 }
00816
00830 template <class Evaluation>
00831 static Evaluation liquidThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
00832 {
00833 const Evaluation& rho = liquidDensity(temperature, pressure);
00834 return Common::thermalConductivityIAPWS(temperature, rho);
00835 }
00836
00850 template <class Evaluation>
00851 static Evaluation gasThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
00852 {
00853 const Evaluation& rho = gasDensity(temperature, pressure);
00854 return Common::thermalConductivityIAPWS(temperature, rho);
00855 }
00856
00857 private:
00858
00859 template <class Evaluation>
00860 static Evaluation enthalpyRegion1_(const Evaluation& temperature, const Evaluation& pressure)
00861 {
00862 return
00863 Region1::tau(temperature) *
00864 Region1::dgamma_dtau(temperature, pressure) *
00865 Rs*temperature;
00866 }
00867
00868
00869 template <class Evaluation>
00870 static Evaluation heatCap_p_Region1_(const Evaluation& temperature, const Evaluation& pressure)
00871 {
00872 return
00873 - Opm::pow(Region1::tau(temperature), 2.0) *
00874 Region1::ddgamma_ddtau(temperature, pressure) *
00875 Rs;
00876 }
00877
00878
00879 template <class Evaluation>
00880 static Evaluation heatCap_v_Region1_(const Evaluation& temperature, const Evaluation& pressure)
00881 {
00882 double tau = Region1::tau(temperature);
00883 double num = Region1::dgamma_dpi(temperature, pressure) - tau * Region1::ddgamma_dtaudpi(temperature, pressure);
00884 double diff = std::pow(num, 2) / Region1::ddgamma_ddpi(temperature, pressure);
00885
00886 return
00887 - std::pow(tau, 2 ) *
00888 Region1::ddgamma_ddtau(temperature, pressure) * Rs +
00889 diff;
00890 }
00891
00892
00893 template <class Evaluation>
00894 static Evaluation internalEnergyRegion1_(const Evaluation& temperature, const Evaluation& pressure)
00895 {
00896 return
00897 Rs * temperature *
00898 ( Region1::tau(temperature)*Region1::dgamma_dtau(temperature, pressure) -
00899 Region1::pi(pressure)*Region1::dgamma_dpi(temperature, pressure));
00900 }
00901
00902
00903 template <class Evaluation>
00904 static Evaluation volumeRegion1_(const Evaluation& temperature, const Evaluation& pressure)
00905 {
00906 return
00907 Region1::pi(pressure)*
00908 Region1::dgamma_dpi(temperature, pressure) *
00909 Rs * temperature / pressure;
00910 }
00911
00912
00913 template <class Evaluation>
00914 static Evaluation enthalpyRegion2_(const Evaluation& temperature, const Evaluation& pressure)
00915 {
00916 return
00917 Region2::tau(temperature) *
00918 Region2::dgamma_dtau(temperature, pressure) *
00919 Rs*temperature;
00920 }
00921
00922
00923 template <class Evaluation>
00924 static Evaluation internalEnergyRegion2_(const Evaluation& temperature, const Evaluation& pressure)
00925 {
00926 return
00927 Rs * temperature *
00928 ( Region2::tau(temperature)*Region2::dgamma_dtau(temperature, pressure) -
00929 Region2::pi(pressure)*Region2::dgamma_dpi(temperature, pressure));
00930 }
00931
00932
00933 template <class Evaluation>
00934 static Evaluation heatCap_p_Region2_(const Evaluation& temperature, const Evaluation& pressure)
00935 {
00936 return
00937 - Opm::pow(Region2::tau(temperature), 2 ) *
00938 Region2::ddgamma_ddtau(temperature, pressure) *
00939 Rs;
00940 }
00941
00942
00943 template <class Evaluation>
00944 static Evaluation heatCap_v_Region2_(const Evaluation& temperature, const Evaluation& pressure)
00945 {
00946 const Evaluation& tau = Region2::tau(temperature);
00947 const Evaluation& pi = Region2::pi(pressure);
00948 const Evaluation& num = 1 + pi * Region2::dgamma_dpi(temperature, pressure) + tau * pi * Region2::ddgamma_dtaudpi(temperature, pressure);
00949 const Evaluation& diff = num * num / (1 - pi * pi * Region2::ddgamma_ddpi(temperature, pressure));
00950 return
00951 - std::pow(tau, 2 ) *
00952 Region2::ddgamma_ddtau(temperature, pressure) * Rs
00953 - diff;
00954 }
00955
00956
00957 template <class Evaluation>
00958 static Evaluation volumeRegion2_(const Evaluation& temperature, const Evaluation& pressure)
00959 {
00960 return
00961 Region2::pi(pressure)*
00962 Region2::dgamma_dpi(temperature, pressure) *
00963 Rs * temperature / pressure;
00964 }
00965 };
00966
00967 template <class Scalar>
00968 const Scalar H2O<Scalar>::Rs = Common::Rs;
00969 }
00970
00971 #endif