diff --git a/Documentation/ProjectFile/material/fluid/viscosity/t_type.md b/Documentation/ProjectFile/material/fluid/viscosity/t_type.md
index 713d4d928b9ca1c320bc43cb5ad114d0b58834e0..26670a401a5c104a644801133f67709093f69f42 100644
--- a/Documentation/ProjectFile/material/fluid/viscosity/t_type.md
+++ b/Documentation/ProjectFile/material/fluid/viscosity/t_type.md
@@ -1,2 +1,2 @@
Type of viscosity model. It can be Constant, LinearPressure, TemperatureDependent,
- or Vogels.
+ Vogels, or WaterViscosityIAPWS
diff --git a/MaterialLib/Fluid/FluidProperty.h b/MaterialLib/Fluid/FluidProperty.h
index fce5a89dc3e900b64e456cd28cc042e618a3bb90..cc3d30876847b6c16756cbacee6e5c86ad27a797 100644
--- a/MaterialLib/Fluid/FluidProperty.h
+++ b/MaterialLib/Fluid/FluidProperty.h
@@ -25,6 +25,7 @@ enum class PropertyVariableType
{
T = 0, ///< temperature.
p = 1, ///< pressure.
+ rho = p, ///< density. For some models, rho substitutes p
number_of_variables = 2 ///< Number of property variables.
};
diff --git a/MaterialLib/Fluid/Viscosity/WaterViscosityIAPWS.cpp b/MaterialLib/Fluid/Viscosity/WaterViscosityIAPWS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6e1e885e039ca9236cf5ac84424aaeb27b0577c
--- /dev/null
+++ b/MaterialLib/Fluid/Viscosity/WaterViscosityIAPWS.cpp
@@ -0,0 +1,203 @@
+/**
+ * \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 WaterViscosityIAPWS.cpp
+ *
+ * Created on December 1, 2016, 1:41 PM
+ */
+
+#include "WaterViscosityIAPWS.h"
+
+#include <array>
+#include <cmath>
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+static const double Hi[4] = {1.67752, 2.20462, 0.6366564, -0.241605};
+static const double Hij[6][7] = {
+ {0.520094, 0.222531, -0.281378, 0.161913, -0.0325372, 0, 0},
+ {0.0850895, 0.999115, -0.906851, 0.257399, 0, 0, 0},
+ {-1.08374, 1.88797, -0.772479, 0, 0, 0, 0},
+ {-0.289555, 1.26613, -0.489837, 0, 0.0698452, 0, -0.00435673},
+ {0, 0, -0.25704, 0, 0, 0.00872102, 0},
+ {0, 0.120573, 0, 0, 0, 0, -0.000593264}};
+
+static double computeBarMu0Factor(const double barT);
+
+static std::array<double, 6> computeSeriesFactorTForMu1(const double barT,
+ const double bar_rho);
+static std::array<double, 7> computeSeriesFactorRhoForMu1(const double barT,
+ const double bar_rho);
+static double computeBarMu1Factor(
+ const std::array<double, 6>& series_factorT,
+ const std::array<double, 7>& series_factorRho);
+
+static double computedBarMu_dbarT(const double barT, double bar_rho);
+static double computedBarMu_dbarRho(const double barT, double bar_rho);
+
+double WaterViscosityIAPWS::getValue(const ArrayType& var_vals) const
+{
+ const double bar_T =
+ var_vals[static_cast<unsigned>(PropertyVariableType::T)] / _ref_T;
+ const double bar_rho =
+ var_vals[static_cast<unsigned>(PropertyVariableType::rho)] / _ref_rho;
+
+ const double mu0 = 100. * std::sqrt(bar_T) / computeBarMu0Factor(bar_T);
+
+ const auto& series_factorT = computeSeriesFactorTForMu1(bar_T, bar_rho);
+ const auto& series_factorRho = computeSeriesFactorRhoForMu1(bar_T, bar_rho);
+ const double mu1 = std::exp(
+ bar_rho * computeBarMu1Factor(series_factorT, series_factorRho));
+
+ return mu0 * mu1 * _ref_mu;
+}
+
+double WaterViscosityIAPWS::getdValue(const ArrayType& var_vals,
+ const PropertyVariableType var_type) const
+{
+ const double bar_T =
+ var_vals[static_cast<unsigned>(PropertyVariableType::T)] / _ref_T;
+ const double bar_rho =
+ var_vals[static_cast<unsigned>(PropertyVariableType::rho)] / _ref_rho;
+
+ switch (var_type)
+ {
+ case PropertyVariableType::T:
+ return _ref_mu * computedBarMu_dbarT(bar_T, bar_rho) / _ref_T;
+ case PropertyVariableType::rho:
+ return _ref_mu * computedBarMu_dbarRho(bar_T, bar_rho) / _ref_rho;
+ default:
+ return 0.;
+ }
+}
+
+double computeBarMu0Factor(const double barT)
+{
+ double sum_val = 0.;
+ double barT_i = 1.;
+ for (int i = 0; i < 4; i++)
+ {
+ sum_val += (Hi[i] / barT_i);
+ barT_i *= barT;
+ }
+ return sum_val;
+}
+
+std::array<double, 6> computeSeriesFactorTForMu1(const double barT,
+ const double bar_rho)
+{
+ std::array<double, 6> series_factorT;
+ series_factorT[0] = 1.;
+ const double barT_fac = 1 / barT - 1.0;
+ for (int i = 1; i < 6; i++)
+ {
+ series_factorT[i] = series_factorT[i - 1] * barT_fac;
+ }
+
+ return series_factorT;
+}
+
+std::array<double, 7> computeSeriesFactorRhoForMu1(const double barT,
+ const double bar_rho)
+{
+ std::array<double, 7> series_factorRho;
+ series_factorRho[0] = 1.;
+ for (int i = 1; i < 7; i++)
+ {
+ series_factorRho[i] = series_factorRho[i - 1] * (bar_rho - 1.0);
+ }
+
+ return series_factorRho;
+}
+
+double computeBarMu1Factor(const std::array<double, 6>& series_factorT,
+ const std::array<double, 7>& series_factorRho)
+{
+ double sum_val = 0.;
+ for (int i = 0; i < 6; i++)
+ {
+ double sum_val_j = 0;
+ for (int j = 0; j < 7; j++)
+ {
+ sum_val_j += Hij[i][j] * series_factorRho[j];
+ }
+ sum_val += series_factorT[i] * sum_val_j;
+ }
+
+ return sum_val;
+}
+
+double computedBarMu_dbarT(const double barT, double bar_rho)
+{
+ const double mu0_factor = computeBarMu0Factor(barT);
+ const double sqrt_barT = std::sqrt(barT);
+
+ double dmu0_factor_dbarT = 0.0;
+ double barT_i = barT * barT;
+ for (int i = 1; i < 4; i++)
+ {
+ dmu0_factor_dbarT -= static_cast<double>(i) * (Hi[i] / barT_i);
+ barT_i *= barT;
+ }
+
+ const double dbar_mu0_dbarT =
+ 50. / (mu0_factor * sqrt_barT) -
+ 100. * sqrt_barT * dmu0_factor_dbarT / (mu0_factor * mu0_factor);
+
+ const auto& series_factorT = computeSeriesFactorTForMu1(barT, bar_rho);
+ const auto& series_factorRho = computeSeriesFactorRhoForMu1(barT, bar_rho);
+
+ double dmu1_factor_dbarT = 0.0;
+ for (int i = 1; i < 6; i++)
+ {
+ double sum_val_j = 0;
+ for (int j = 0; j < 7; j++)
+ {
+ sum_val_j += Hij[i][j] * series_factorRho[j];
+ }
+ dmu1_factor_dbarT -= static_cast<double>(i) * series_factorT[i - 1] *
+ sum_val_j / (barT * barT);
+ }
+
+ const double mu1_factor =
+ computeBarMu1Factor(series_factorT, series_factorRho);
+ const double dbar_mu1_dbarT =
+ bar_rho * std::exp(bar_rho * mu1_factor) * dmu1_factor_dbarT;
+
+ return dbar_mu0_dbarT * std::exp(bar_rho * mu1_factor) +
+ dbar_mu1_dbarT * 100. * sqrt_barT / mu0_factor;
+}
+
+double computedBarMu_dbarRho(const double barT, double bar_rho)
+{
+ const auto& series_factorT = computeSeriesFactorTForMu1(barT, bar_rho);
+ const auto& series_factorRho = computeSeriesFactorRhoForMu1(barT, bar_rho);
+
+ double dmu1_factor_dbar_rho = 0.0;
+ for (int i = 0; i < 6; i++)
+ {
+ double sum_val_j = 0;
+ for (int j = 1; j < 7; j++)
+ {
+ sum_val_j +=
+ static_cast<double>(j) * Hij[i][j] * series_factorRho[j - 1];
+ }
+ dmu1_factor_dbar_rho += series_factorT[i] * sum_val_j;
+ }
+
+ const double mu0 = 100. * std::sqrt(barT) / computeBarMu0Factor(barT);
+
+ const double mu1_factor =
+ computeBarMu1Factor(series_factorT, series_factorRho);
+ return mu0 * std::exp(bar_rho * mu1_factor) *
+ (mu1_factor + bar_rho * dmu1_factor_dbar_rho);
+}
+
+} // end namespace
+} // end namespace
diff --git a/MaterialLib/Fluid/Viscosity/WaterViscosityIAPWS.h b/MaterialLib/Fluid/Viscosity/WaterViscosityIAPWS.h
new file mode 100644
index 0000000000000000000000000000000000000000..2411f9f8762d9868dc844776e258d37a4da11983
--- /dev/null
+++ b/MaterialLib/Fluid/Viscosity/WaterViscosityIAPWS.h
@@ -0,0 +1,89 @@
+/**
+ * \brief Viscosity model according to
+ * Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary
+ * Water Substance
+ *
+ * \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 WaterViscosityIAPWS.h
+ *
+ * Created on December 1, 2016, 1:41 PM
+ */
+
+#ifndef OGS_WATER_VISCOSITY_IAPWS_H
+#define OGS_WATER_VISCOSITY_IAPWS_H
+
+#include <string>
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/**
+ * \brief A class for viscosity model that is defined by
+ * The International Association for the Properties of Water and Steam
+ * <a href="http://www.iapws.org/relguide/visc.pdf">IAPWS</a>
+ *
+ * With the definition, the viscosity is a function of temperature and
+ * water density
+ *
+ * \attention The critical enhancement, \f$\bar{\mu}_2\f$, which is significant
+ * only within the boundaries specified by
+ * \f[ T (\mbox{in K}) \in (645.91, 650.77) \f]
+ * and
+ * \f[ \rho (\mbox{in kg m}^{-3}) \in (245.8, 405.3)\f],
+ * is not considered.
+ */
+class WaterViscosityIAPWS final : public FluidProperty
+{
+public:
+ WaterViscosityIAPWS() = default;
+
+ /// Get model name.
+ std::string getName() const override
+ {
+ return "IAPWS temperature-density dependent viscosity model";
+ }
+
+ /**
+ * Get density value.
+ * \param var_vals Variable values of temperature and water density in
+ * an array. The order of its elements
+ * is given in enum class PropertyVariableType.
+ *
+ * \attention The variable values must follow the SI standard,
+ * e.g temperature unit is K.
+ */
+ double getValue(const ArrayType& var_vals) const override;
+
+ /** Get the partial differential of the density with respect to temperature
+ * or liquid pressure.
+ * \param var_vals Variable values of temperature and water density in
+ * an array. The order of its elements
+ * is given in enum class PropertyVariableType.
+ *
+ * \attention The variable values must follow the SI standard,
+ * e.g temperature unit is K.
+ * \param var_type Variable type.
+ */
+ double getdValue(const ArrayType& var_vals,
+ const PropertyVariableType var_type) const override;
+
+private:
+ const double _ref_T = 647.096; ///< reference temperature in K
+ const double _ref_rho = 322.0; ///< reference density in `kg/m^3`
+ const double _ref_mu = 1.0e-6; ///< reference viscosity in Pa.s
+
+ // Coefficients Hi and Hij are given in two static arrays in the cpp file.
+};
+
+} // end namespace
+} // end namespace
+
+#endif /* OGS_WATER_VISCOSITY_IAPWS_H */
diff --git a/MaterialLib/Fluid/Viscosity/createViscosityModel.cpp b/MaterialLib/Fluid/Viscosity/createViscosityModel.cpp
index 7f70b56b7c4d910536a4cb6da330f7e4b67932f4..740912ef2c48b7b177eea4a1d4033e4db7831159 100644
--- a/MaterialLib/Fluid/Viscosity/createViscosityModel.cpp
+++ b/MaterialLib/Fluid/Viscosity/createViscosityModel.cpp
@@ -19,6 +19,7 @@
#include "LinearPressureDependentViscosity.h"
#include "TemperatureDependentViscosity.h"
#include "VogelsLiquidDynamicViscosity.h"
+#include "WaterViscosityIAPWS.h"
namespace MaterialLib
{
@@ -135,6 +136,12 @@ std::unique_ptr<FluidProperty> createViscosityModel(
fluid_type.data());
}
}
+ else if (type == "WaterViscosityIAPWS")
+ {
+ //! \ogs_file_param{material__fluid__viscosity__type}
+ config.checkConfigParameter("type", "WaterViscosityIAPWS");
+ return std::unique_ptr<FluidProperty>(new WaterViscosityIAPWS());
+ }
else
{
OGS_FATAL(
diff --git a/Tests/MaterialLib/TestFluidViscosity.cpp b/Tests/MaterialLib/TestFluidViscosity.cpp
index 6be91542b0ed56b54fc33d9285085c56b94eb04b..8cb6686322df087665120196c70343384ab6bbd0 100644
--- a/Tests/MaterialLib/TestFluidViscosity.cpp
+++ b/Tests/MaterialLib/TestFluidViscosity.cpp
@@ -45,8 +45,8 @@ TEST(Material, checkConstantViscosity)
ArrayType dummy;
ASSERT_EQ(1.e-4, mu->getValue(dummy));
- ASSERT_EQ(0.0,
- mu->getdValue(dummy, MaterialLib::Fluid::PropertyVariableType::T));
+ ASSERT_EQ(
+ 0.0, mu->getdValue(dummy, MaterialLib::Fluid::PropertyVariableType::T));
}
TEST(Material, checkTemperatureDependentViscosity)
@@ -64,9 +64,10 @@ TEST(Material, checkTemperatureDependentViscosity)
vars[0] = 350.0;
const double mu_expected = 1.e-3 * std::exp(-(vars[0] - 293) / 368);
ASSERT_NEAR(mu_expected, mu->getValue(vars), 1.e-10);
- ASSERT_NEAR(-mu_expected,
- mu->getdValue(vars, MaterialLib::Fluid::PropertyVariableType::T),
- 1.e-10);
+ ASSERT_NEAR(
+ -mu_expected,
+ mu->getdValue(vars, MaterialLib::Fluid::PropertyVariableType::T),
+ 1.e-10);
}
TEST(Material, checkLinearPressureDependentViscosity)
@@ -85,9 +86,10 @@ TEST(Material, checkLinearPressureDependentViscosity)
vars[1] = 2.e+6;
ASSERT_NEAR(1.e-3 * (1. + 1.e-6 * (vars[1] - 1.e+5)), mu->getValue(vars),
1.e-10);
- ASSERT_NEAR(1.e-9,
- mu->getdValue(vars, MaterialLib::Fluid::PropertyVariableType::p),
- 1.e-10);
+ ASSERT_NEAR(
+ 1.e-9,
+ mu->getdValue(vars, MaterialLib::Fluid::PropertyVariableType::p),
+ 1.e-10);
}
TEST(Material, checkVogelViscosity)
@@ -124,3 +126,50 @@ TEST(Material, checkVogelViscosity)
ASSERT_NEAR(0.352072e-4, mu_ch4->getValue(vars), 1.e-5);
ASSERT_NEAR(-2.35664e-6, mu_ch4->getdValue(vars, var_type), 1.e-5);
}
+
+TEST(Material, checkWaterViscosityIAPWS)
+{
+ const char xml_w[] =
+ "<viscosity>"
+ " <type>WaterViscosityIAPWS</type>"
+ "</viscosity>";
+
+ // Test data provided on http://www.iapws.org/relguide/visc.pdf
+ const auto mu_w = createTestViscosityModel(xml_w);
+ const double T[] = {298.15, 298.15, 373.15, 433.15, 433.15, 873.15,
+ 873.15, 873.15, 1173.15, 1173.15, 1173.15};
+ const double rho[] = {998, 1200, 1000, 1, 1000, 1, 100, 600, 1, 100, 400};
+
+ const double expected_mumu[] = {
+ 889.735100, 1437.649467, 307.883622, 14.538324, 217.685358, 32.619287,
+ 35.802262, 77.430195, 44.217245, 47.640433, 64.154608};
+
+ const double perturbation = 1.e-9;
+ ArrayType vars;
+ for (int i = 0; i < 11; i++)
+ {
+ // Test mu
+ vars[static_cast<unsigned>(PropertyVariableType::T)] = T[i];
+ vars[static_cast<unsigned>(PropertyVariableType::rho)] = rho[i];
+ const double mu = mu_w->getValue(vars);
+ ASSERT_NEAR(expected_mumu[i] * 1.e-6, mu, 1.e-9);
+
+ const double dmu_dT = mu_w->getdValue(vars, PropertyVariableType::T);
+ const double dmu_drho =
+ mu_w->getdValue(vars, PropertyVariableType::rho);
+
+ // Test dmu/dT
+ vars[static_cast<unsigned>(PropertyVariableType::T)] =
+ T[i] + perturbation;
+ double mu1 = mu_w->getValue(vars);
+ ASSERT_NEAR((mu1 - mu) / perturbation, dmu_dT, 1.e-7);
+
+ // Test dmu/drho
+ vars[static_cast<unsigned>(PropertyVariableType::T)] = T[i];
+ vars[static_cast<unsigned>(PropertyVariableType::rho)] =
+ rho[i] + perturbation;
+ mu1 = mu_w->getValue(vars);
+
+ ASSERT_NEAR((mu1 - mu) / perturbation, dmu_drho, 1.e-7);
+ }
+}