diff --git a/Documentation/ProjectFile/material/fluid/density/t_type.md b/Documentation/ProjectFile/material/fluid/density/t_type.md
index e7c192685cc3b1f5aeca75dc2c3f045d53c7b3e6..94821d5db9c09283a616f2d898f95f4c9fa2c937 100644
--- a/Documentation/ProjectFile/material/fluid/density/t_type.md
+++ b/Documentation/ProjectFile/material/fluid/density/t_type.md
@@ -1,2 +1,2 @@
Type of density model. It can be Constant, LiquidDensity, TemperatureDependent,
-or IdealGasLaw.
\ No newline at end of file
+ IdealGasLaw, or WaterDensityIAPWSIF97Region1.
\ No newline at end of file
diff --git a/MaterialLib/Fluid/Density/WaterDensityIAPWSIF97Region1.h b/MaterialLib/Fluid/Density/WaterDensityIAPWSIF97Region1.h
new file mode 100644
index 0000000000000000000000000000000000000000..2d08fc68884ca0a8629ea0992d7ce6445304c364
--- /dev/null
+++ b/MaterialLib/Fluid/Density/WaterDensityIAPWSIF97Region1.h
@@ -0,0 +1,96 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file WaterDensityIAPWSIF97Region1.h
+ *
+ * Created on December 8, 2016, 4:19 PM
+ */
+
+#ifndef OGS_WATER_DENSITY_IAPWSIF97_REGION1_H
+#define OGS_WATER_DENSITY_IAPWSIF97_REGION1_H
+
+#include <string>
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/** Water density model
+ * base on the IAPWS Industrial Formulation 1997
+ * <a href="http://www.iapws.org/relguide/IF97-Rev.pdf">IF97-Rev</a>
+ */
+class WaterDensityIAPWSIF97Region1 final : public FluidProperty
+{
+public:
+ WaterDensityIAPWSIF97Region1()
+ : _gibbs_free_energy(DimensionLessGibbsFreeEnergyRegion1()){};
+
+ /// Get density model name.
+ std::string getName() const override
+ {
+ return "Water density model based on IAPWSIF97";
+ }
+
+ /// Get density value.
+ /// \param var_vals Variable values in an array of temperature and pressure.
+ double getValue(const ArrayType& var_vals) const override
+ {
+ const double T = var_vals[static_cast<int>(PropertyVariableType::T)];
+ const double p = var_vals[static_cast<int>(PropertyVariableType::p)];
+ const double tau = _ref_T / T;
+ const double pi = p / _ref_p;
+
+ return _ref_p / (_sR * T * _gibbs_free_energy.get_dgamma_dpi(tau, pi));
+ }
+
+ /**
+ * Get the partial differential of the density with respect to temperature
+ * or pressure.
+ * \param var_vals Variable values in an array of temperature and pressure.
+ * \param var_type Variable type.
+ */
+ double getdValue(const ArrayType& var_vals,
+ const PropertyVariableType var_type) const override
+ {
+ const double T = var_vals[static_cast<int>(PropertyVariableType::T)];
+ const double p = var_vals[static_cast<int>(PropertyVariableType::p)];
+
+ const double tau = _ref_T / T;
+ const double pi = p / _ref_p;
+
+ const double dgamma_dpi = _gibbs_free_energy.get_dgamma_dpi(tau, pi);
+ switch (var_type)
+ {
+ case PropertyVariableType::T:
+ return -(_ref_p -
+ tau * _ref_p *
+ _gibbs_free_energy.get_dgamma_dtau_dpi(tau, pi) /
+ dgamma_dpi) /
+ (_sR * T * T * dgamma_dpi);
+ case PropertyVariableType::p:
+ return -_gibbs_free_energy.get_dgamma_dpi_dpi(tau, pi) /
+ (_sR * T * dgamma_dpi * dgamma_dpi);
+ default:
+ return 0.0;
+ }
+ }
+
+private:
+ const DimensionLessGibbsFreeEnergyRegion1 _gibbs_free_energy;
+
+ const double _ref_T = 1386; ///< reference temperature in K.
+ const double _ref_p = 1.653e7; ///< reference pressure in Pa.
+ /// Specific water vapour gas constant in J/(kgK).
+ const double _sR = 461.526;
+};
+
+} // end of namespace
+} // end of namespace
+#endif /* OGS_WATER_DENSITY_IAPWSIF97_REGION1_H */
diff --git a/MaterialLib/Fluid/Density/createFluidDensityModel.cpp b/MaterialLib/Fluid/Density/createFluidDensityModel.cpp
index dae39c447b662c648c7486df9db342799f9fc912..89420f159a59cc9c5638fc55fc36cf044b9f68cd 100644
--- a/MaterialLib/Fluid/Density/createFluidDensityModel.cpp
+++ b/MaterialLib/Fluid/Density/createFluidDensityModel.cpp
@@ -18,6 +18,8 @@
#include "IdealGasLaw.h"
#include "LinearTemperatureDependentDensity.h"
#include "LiquidDensity.h"
+#include "WaterDensityIAPWSIF97Region1.h"
+
#include "MaterialLib/Fluid/ConstantFluidProperty.h"
namespace MaterialLib
@@ -96,12 +98,18 @@ std::unique_ptr<FluidProperty> createFluidDensityModel(
//! \ogs_file_param{material__fluid__density__IdealGasLaw__molar_mass}
new IdealGasLaw(config.getConfigParameter<double>("molar_mass")));
}
+ else if (type == "WaterDensityIAPWSIF97Region1")
+ {
+ return std::unique_ptr<FluidProperty>(
+ new WaterDensityIAPWSIF97Region1());
+ }
else
{
OGS_FATAL(
"The density type %s is unavailable.\n"
"The available types are: \n\tConstant, \n\tLiquidDensity, "
- "\n\tTemperatureDependent, \n\tIdealGasLaw.\n",
+ "\n\tTemperatureDependent, \n\tIdealGasLaw."
+ "\n\tWaterDensityIAPWSIF97Region1\n",
type.data());
}
}
diff --git a/Tests/MaterialLib/TestFluidDensity.cpp b/Tests/MaterialLib/TestFluidDensity.cpp
index a9aa4096df354733cdb11c6be2ec2dfe2b5502e1..4baac6f07be555bfddd08897ee5b2553077aa57c 100644
--- a/Tests/MaterialLib/TestFluidDensity.cpp
+++ b/Tests/MaterialLib/TestFluidDensity.cpp
@@ -132,3 +132,31 @@ TEST(Material, checkLiquidDensity)
ASSERT_NEAR(rho0 / (1. + beta * (T - T0)) / (fac_p * fac_p * K),
rho->getdValue(vars, Fluid::PropertyVariableType::p), 1.e-10);
}
+
+TEST(Material, checkWaterDensityIAPWSIF97Region1)
+{
+ const char xml[] =
+ "<density>"
+ " <type>WaterDensityIAPWSIF97Region1</type>"
+ "</density>";
+ const auto rho = createTestFluidDensityModel(xml);
+
+ ArrayType vars = {{473.15, 4.e+7}};
+ const double rho_expected = 890.943136237744;
+ ASSERT_NEAR(rho_expected, rho->getValue(vars), 1.e-10);
+
+ const double drho_dT = rho->getdValue(vars, PropertyVariableType::T);
+ const double drho_dp = rho->getdValue(vars, PropertyVariableType::p);
+
+ const double perturbation = 1.e-4;
+
+ // Test the differentiation: with respect to temperature:
+ vars[static_cast<unsigned>(PropertyVariableType::T)] += perturbation;
+ const double rho_T1 = rho->getValue(vars);
+ ASSERT_NEAR((rho_T1 - rho_expected) / perturbation, drho_dT, 1.e-6);
+
+ // Test the differentiation: with respect to pressure:
+ vars[static_cast<unsigned>(PropertyVariableType::p)] += perturbation;
+ const double rho_p1 = rho->getValue(vars);
+ ASSERT_NEAR((rho_p1 - rho_T1) / perturbation, drho_dp, 1.e-6);
+}