add simplified elasticity models for ThermoRichardsFlow process

ProcessLib/ThermoRichardsFlow/CreateSimplifiedElasticityModel.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateSimplifiedElasticityModel.h"
+
+#include <Eigen/Dense>
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Logging.h"
+#include "SimplifiedElasticityModel.h"
+#include "HydrostaticElasticityModel.h"
+#include "RigidElasticityModel.h"
+#include "UniaxialElasticityModel.h"
+#include "UserDefinedElasticityModel.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+std::unique_ptr<SimplifiedElasticityModel> createElasticityModel(
+    BaseLib::ConfigTree const& config)
+{
+    std::unique_ptr<SimplifiedElasticityModel> simplified_elasticity =
+        std::make_unique<RigidElasticityModel>();
+    if (auto const simplified_elasticity_switch =
+            //! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__simplified_elasticity}
+        config.getConfigParameterOptional<std::string>("simplified_elasticity"))
+    {
+        DBUG("Using simplified_elasticity for the Richards flow equation");
+        if (*simplified_elasticity_switch == "uniaxial")
+        {
+            DBUG("assuming local uniaxial deformation only.");
+            simplified_elasticity = std::make_unique<UniaxialElasticityModel>();
+        }
+        else if (*simplified_elasticity_switch == "hydrostatic")
+        {
+            DBUG("assuming constant hydrostatic stress locally.");
+            simplified_elasticity =
+                std::make_unique<HydrostaticElasticityModel>();
+        }
+        else if (*simplified_elasticity_switch == "user_defined")
+        {
+            DBUG("using user defined elasticity model.");
+            simplified_elasticity =
+                std::make_unique<UserDefinedElasticityModel>();
+        }
+        else if (*simplified_elasticity_switch == "rigid")
+        {
+            DBUG("using user defined elasticity model.");
+            simplified_elasticity = std::make_unique<RigidElasticityModel>();
+        }
+    }
+    return simplified_elasticity;
+}
+}  // namespace ThermoRichardsFlow
+}  // namespace ProcessLib

ProcessLib/ThermoRichardsFlow/CreateSimplifiedElasticityModel.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+#pragma once
+
+#include <memory>
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct SimplifiedElasticityModel;
+}
+}
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+std::unique_ptr<SimplifiedElasticityModel> createElasticityModel(
+    BaseLib::ConfigTree const& config);
+
+}  // namespace ThermoRichardsFlow
+}  // namespace ProcessLib

ProcessLib/ThermoRichardsFlow/HydrostaticElasticityModel.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on August 14, 2020, 10:56 AM
+ */
+
+#pragma once
+
+#include "SimplifiedElasticityModel.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct HydrostaticElasticityModel : SimplifiedElasticityModel
+{
+    HydrostaticElasticityModel()
+    {
+        DBUG("using hydrostatic simplified mechanics model");
+    }
+
+    double storageContribution(
+        MaterialPropertyLib::Phase const& solid_phase,
+        MaterialPropertyLib::VariableArray const& variable_array,
+        ParameterLib::SpatialPosition const& pos, double const t,
+        double const dt) override
+    {
+        return bulkCompressibilityFromYoungsModulus(
+            solid_phase, variable_array, pos, t, dt);
+    }
+
+    double thermalExpansivityContribution(
+        Eigen::Matrix<double, 3, 3> const& solid_linear_thermal_expansion_coefficient,
+        MaterialPropertyLib::Phase const& /*solid_phase*/,
+        MaterialPropertyLib::VariableArray const& /*variable_array*/,
+        ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+        double const /*dt*/) override
+    {
+        return -solid_linear_thermal_expansion_coefficient.trace();
+    }
+};
+
+}  // namespace ThermoRichardsFlow
+}  // namespace ProcessLib

ProcessLib/ThermoRichardsFlow/RigidElasticityModel.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on August 14, 2020, 10:56 AM
+ */
+
+#pragma once
+
+#include "SimplifiedElasticityModel.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct RigidElasticityModel : SimplifiedElasticityModel
+{
+    RigidElasticityModel() { DBUG("using rigid mechanics model"); }
+
+    double storageContribution(MaterialPropertyLib::Phase const&,
+                               MaterialPropertyLib::VariableArray const&,
+                               ParameterLib::SpatialPosition const&,
+                               double const, double const) override
+    {
+        return 0.0;
+    }
+
+    double thermalExpansivityContribution(
+        Eigen::Matrix<double, 3, 3> const&, MaterialPropertyLib::Phase const&,
+        MaterialPropertyLib::VariableArray const&,
+        ParameterLib::SpatialPosition const&, double const,
+        double const) override
+    {
+        return 0.0;
+    }
+};
+
+}  // namespace ThermoRichardsFlow
+}  // namespace ProcessLib

ProcessLib/ThermoRichardsFlow/SimplifiedElasticityModel.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on August 14, 2020, 10:56 AM
+ */
+
+#pragma once
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEigenVector.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct SimplifiedElasticityModel
+{
+    virtual ~SimplifiedElasticityModel() = default;
+    virtual double storageContribution(
+        MaterialPropertyLib::Phase const&,
+        MaterialPropertyLib::VariableArray const&,
+        ParameterLib::SpatialPosition const&, double const, double const) = 0;
+    virtual double thermalExpansivityContribution(
+        Eigen::Matrix<double, 3, 3> const&, MaterialPropertyLib::Phase const&,
+        MaterialPropertyLib::VariableArray const&,
+        ParameterLib::SpatialPosition const&, double const, double const) = 0;
+
+    static inline auto bulkCompressibilityFromYoungsModulus(
+        MaterialPropertyLib::Phase const& solid_phase,
+        MaterialPropertyLib::VariableArray const& variables,
+        ParameterLib::SpatialPosition const& x_position, double const t,
+        double const dt)
+    {
+        // assuming: nu[0]=nu(1,2), nu[1]=nu(2,3), nu[2]=nu(1,3)
+        if (!solid_phase.hasProperty(
+                MaterialPropertyLib::PropertyType::youngs_modulus))
+        {
+            return 0.0;
+        }
+        auto const E = MaterialPropertyLib::formEigenVector<3>(
+            solid_phase[MaterialPropertyLib::PropertyType::youngs_modulus]
+                .value(variables, x_position, t, dt));
+        auto const nu = MaterialPropertyLib::formEigenVector<3>(
+            solid_phase[MaterialPropertyLib::PropertyType::poissons_ratio]
+                .value(variables, x_position, t, dt));
+        return (E[0] * E[1] + E[0] * E[2] * (1 - 2 * nu[1]) +
+                E[1] * E[2] * (1 - 2 * nu[0] - 2 * nu[2])) /
+               (E[0] * E[1] * E[2]);
+    }
+};
+
+}  // namespace ThermoRichardsFlow
+}  // namespace ProcessLib

ProcessLib/ThermoRichardsFlow/UniaxialElasticityModel.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on August 14, 2020, 10:56 AM
+ */
+
+#pragma once
+
+#include "SimplifiedElasticityModel.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct UniaxialElasticityModel : SimplifiedElasticityModel
+{
+    UniaxialElasticityModel()
+    {
+        DBUG("using uniaxial simplified mechanics model");
+    }
+
+    double storageContribution(
+        MaterialPropertyLib::Phase const& solid_phase,
+        MaterialPropertyLib::VariableArray const& variables,
+        ParameterLib::SpatialPosition const& x_position, double const t,
+        double const dt) override
+    {
+        auto const E = MaterialPropertyLib::formEigenVector<3>(
+            solid_phase[MaterialPropertyLib::PropertyType::youngs_modulus]
+                .value(variables, x_position, t, dt));
+        auto const nu = MaterialPropertyLib::formEigenVector<3>(
+            solid_phase[MaterialPropertyLib::PropertyType::poissons_ratio]
+                .value(variables, x_position, t, dt));
+        auto const nu12 = nu[0];
+        auto const nu23 = nu[1];
+        auto const nu13 = nu[2];
+        auto const nu21 = nu12 * E[1] / E[0];
+        auto const nu32 = nu23 * E[2] / E[1];
+        auto const nu31 = nu13 * E[2] / E[0];
+        auto const D = 1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 -
+                       2 * nu12 * nu23 * nu31;
+        return D / (E[2] * (1 - nu12 * nu21));
+    }
+
+    double thermalExpansivityContribution(
+        Eigen::Matrix<double, 3, 3> const& solid_linear_thermal_expansion_coefficient,
+        MaterialPropertyLib::Phase const& solid_phase,
+        MaterialPropertyLib::VariableArray const& variables,
+        ParameterLib::SpatialPosition const& x_position, double const t,
+        double const dt) override
+    {
+        auto const E = MaterialPropertyLib::formEigenVector<3>(
+            solid_phase[MaterialPropertyLib::PropertyType::youngs_modulus]
+                .value(variables, x_position, t, dt));
+        auto const nu = MaterialPropertyLib::formEigenVector<3>(
+            solid_phase[MaterialPropertyLib::PropertyType::poissons_ratio]
+                .value(variables, x_position, t, dt));
+        auto const nu12 = nu[0];
+        auto const nu23 = nu[1];
+        auto const nu13 = nu[2];
+        auto const nu21 = nu12 * E[1] / E[0];
+        auto const D = (1 - nu12 * nu21);
+        return -(solid_linear_thermal_expansion_coefficient(2, 2) +
+                 solid_linear_thermal_expansion_coefficient(0, 0) *
+                     (nu13 + nu12 * nu23) / D +
+                 solid_linear_thermal_expansion_coefficient(1, 1) *
+                     (nu23 + nu13 * nu21) / D);
+    }
+};
+
+}  // namespace ThermoRichardsFlow
+}  // namespace ProcessLib

ProcessLib/ThermoRichardsFlow/UserDefinedElasticityModel.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on August 14, 2020, 10:56 AM
+ */
+
+#pragma once
+
+#include "SimplifiedElasticityModel.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct UserDefinedElasticityModel : SimplifiedElasticityModel
+{
+    UserDefinedElasticityModel()
+    {
+        DBUG("using user defined simplified elasticity model");
+    }
+
+    double storageContribution(
+        MaterialPropertyLib::Phase const& solid_phase,
+        MaterialPropertyLib::VariableArray const& variables,
+        ParameterLib::SpatialPosition const& x_position, double const t,
+        double const dt) override
+    {
+        return solid_phase[MaterialPropertyLib::PropertyType::storage_contribution]
+            .template value<double>(variables, x_position, t, dt);
+    }
+    double thermalExpansivityContribution(
+        Eigen::Matrix<double, 3,
+                      3> const& /*solid_linear_thermal_expansion_coefficient*/,
+        MaterialPropertyLib::Phase const& solid_phase,
+        MaterialPropertyLib::VariableArray const& variables,
+        ParameterLib::SpatialPosition const& x_position, double const t,
+        double const dt) override
+    {
+        return solid_phase[MaterialPropertyLib::PropertyType::
+                          thermal_expansivity_contribution]
+            .template value<double>(variables, x_position, t, dt);
+    }
+};
+
+}  // namespace ThermoRichardsFlow
+}  // namespace ProcessLib