Commit 6f17c15a authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

Merge branch 'ThermoRichardsFlow' into 'master'

New Process: Non-isothermal Richards  flow

See merge request ogs/ogs!3419
parents 32f04f9c faa9307c
......@@ -118,6 +118,9 @@
#ifdef OGS_BUILD_PROCESS_THERMOMECHANICS
#include "ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h"
#endif
#ifdef OGS_BUILD_PROCESS_THERMORICHARDSFLOW
#include "ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.h"
#endif
#ifdef OGS_BUILD_PROCESS_TWOPHASEFLOWWITHPP
#include "ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h"
#endif
......@@ -1047,7 +1050,17 @@ void ProjectData::parseProcesses(
}
else
#endif
#ifdef OGS_BUILD_PROCESS_THERMORICHARDSFLOW
if (type == "THERMO_RICHARDS_FLOW")
{
process =
ProcessLib::ThermoRichardsFlow::createThermoRichardsFlowProcess(
name, *_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, integration_order,
process_config, _media);
}
else
#endif
#ifdef OGS_BUILD_PROCESS_THERMORICHARDSMECHANICS
if (type == "THERMO_RICHARDS_MECHANICS")
{
......
\copydoc ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowProcess
Coupling scheme. So far, only the full monolithic scheme is available.
A tag for enabling diagonal lumping of the mass matrix of the Richards equation.
Considertion of mechanics in the mass balance equation under specific stress conditions
to approximate the \f$ \nabla \cdot \dot {\mathbf u}\f$-term. See \cite Buchwald2021 for details.
\copydoc ProcessLib::RichardsMechanics::RichardsMechanicsProcessData::specific_body_force
......@@ -10,7 +10,17 @@
Year={1964},
Publisher={colorado State University}
}
@Article{Buchwald2021,
title = {Improved predictions of thermal fluid pressurization in hydro-thermal models based on consistent incorporation of thermo-mechanical effects in anisotropic porous media},
journal = {International Journal of Heat and Mass Transfer},
volume = {172},
pages = {121127},
year = {2021},
issn = {0017-9310},
doi = {https://doi.org/10.1016/j.ijheatmasstransfer.2021.121127},
url = {https://www.sciencedirect.com/science/article/pii/S0017931021002301},
author = {J. Buchwald and S. Kaiser and O. Kolditz and T. Nagel}
}
@Article{Ehlers1995,
Title={A single-surface yield function for geomaterials},
Author={Ehlers, W},
......
......@@ -68,6 +68,7 @@ enum PropertyType : int
/// ion diffusivity in the porous medium with account of the effect of
/// tortuosity and connectivity.
pore_diffusion,
poissons_ratio,
porosity,
reference_density,
reference_temperature,
......@@ -84,9 +85,13 @@ enum PropertyType : int
specific_heat_capacity,
specific_latent_heat,
storage,
storage_contribution,
swelling_stress_rate,
thermal_conductivity,
/// Thermal diffusion enhancement factor for water vapor flow
thermal_diffusion_enhancement_factor,
thermal_expansivity,
thermal_expansivity_contribution,
thermal_longitudinal_dispersivity,
thermal_osmosis_coefficient,
thermal_transversal_dispersivity,
......@@ -98,6 +103,7 @@ enum PropertyType : int
vapour_diffusion,
viscosity,
volume_fraction,
youngs_modulus,
number_of_properties
};
......@@ -134,6 +140,7 @@ static const std::array<std::string, PropertyType::number_of_properties>
"permeability",
"phase_velocity",
"pore_diffusion",
"poissons_ratio",
"porosity",
"reference_density",
"reference_temperature",
......@@ -148,9 +155,12 @@ static const std::array<std::string, PropertyType::number_of_properties>
"specific_heat_capacity",
"specific_latent_heat",
"storage",
"storage_contribution",
"swelling_stress_rate",
"thermal_conductivity",
"thermal_diffusion_enhancement_factor",
"thermal_expansivity",
"thermal_expansivity_contribution",
"thermal_longitudinal_dispersivity",
"thermal_osmosis_coefficient",
"thermal_transversal_dispersivity",
......@@ -160,7 +170,8 @@ static const std::array<std::string, PropertyType::number_of_properties>
"vapour_density",
"vapour_diffusion",
"viscosity",
"volume_fraction"}};
"volume_fraction",
"youngs_modulus"}};
/// This function converts a string (e.g. a string from the configuration-tree)
/// into one of the entries of the PropertyType enumerator.
......
get_source_files(SOURCES)
ogs_add_library(ThermoRichardsFlow ${SOURCES})
target_link_libraries(ThermoRichardsFlow PUBLIC ProcessLib PRIVATE ParameterLib)
if(OGS_BUILD_TESTING)
include(Tests.cmake)
endif()
/**
* \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
/**
* \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
/**
* \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 "CreateThermoRichardsFlowProcess.h"
#include <cassert>
#include "CreateSimplifiedElasticityModel.h"
#include "LocalAssemblerInterface.h"
#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "MaterialLib/SolidModels/MechanicsBase.h"
#include "ParameterLib/Utils.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
#include "ProcessLib/Utils/ProcessUtils.h"
#include "SimplifiedElasticityModel.h"
#include "ThermoRichardsFlowProcess.h"
#include "ThermoRichardsFlowProcessData.h"
namespace ProcessLib
{
namespace ThermoRichardsFlow
{
void checkMPLProperties(
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
{
std::array const required_medium_properties = {
MaterialPropertyLib::permeability,
MaterialPropertyLib::porosity, MaterialPropertyLib::biot_coefficient,
MaterialPropertyLib::relative_permeability,
MaterialPropertyLib::saturation};
std::array const required_liquid_properties = {
MaterialPropertyLib::viscosity,
MaterialPropertyLib::density,
};
std::array const required_solid_properties = {
MaterialPropertyLib::density};
// Thermal properties are not checked because they can be phase property or
// meduim property (will be enabled later).
for (auto const& m : media)
{
checkRequiredProperties(*m.second, required_medium_properties);
checkRequiredProperties(m.second->phase("AqueousLiquid"),
required_liquid_properties);
checkRequiredProperties(m.second->phase("Solid"),
required_solid_properties);
}
}
void checkProcessVariableComponents(ProcessVariable const& variable)
{
if (variable.getNumberOfGlobalComponents() != 1)
{
OGS_FATAL(
"Number of components of the process variable '{:s}' is different "
"from one: got {:d}.",
variable.getName(),
variable.getNumberOfGlobalComponents());
}
}
std::unique_ptr<Process> createThermoRichardsFlowProcess(
std::string name,
MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
unsigned const integration_order,
BaseLib::ConfigTree const& config,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
{
//! \ogs_file_param{prj__processes__process__type}
config.checkConfigParameter("type", "THERMO_RICHARDS_FLOW");
DBUG("Create ThermoRichardsFlowProcess.");
auto const staggered_scheme =
//! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__coupling_scheme}
config.getConfigParameterOptional<std::string>("coupling_scheme");
const bool use_monolithic_scheme =
!(staggered_scheme && (*staggered_scheme == "staggered"));
// Process variable.
//! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__process_variables}
auto const pv_config = config.getConfigSubtree("process_variables");
ProcessVariable* variable_T;
ProcessVariable* variable_p;
std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
process_variables;
if (use_monolithic_scheme) // monolithic scheme.
{
auto per_process_variables = findProcessVariables(
variables, pv_config,
{//! \ogs_file_param_special{prj__processes__process__THERMO_RICHARDS_FLOW__process_variables__temperature}
"temperature",
//! \ogs_file_param_special{prj__processes__process__THERMO_RICHARDS_FLOW__process_variables__pressure}
"pressure"});
variable_T = &per_process_variables[0].get();
variable_p = &per_process_variables[1].get();
process_variables.push_back(std::move(per_process_variables));
}
else // staggered scheme.
{
OGS_FATAL(
"So far, only the monolithic scheme is implemented for "
"THERMO_RICHARDS_FLOW");
}
checkProcessVariableComponents(*variable_T);
checkProcessVariableComponents(*variable_p);
// Specific body force parameter.
Eigen::VectorXd specific_body_force;
{
std::vector<double> const b =
//! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__specific_body_force}
config.getConfigParameter<std::vector<double>>(
"specific_body_force");
if (b.size() != mesh.getDimension())
{
OGS_FATAL(
"specific body force (gravity vector) has {:d} components, "
"but mesh dimension is {:d}",
b.size(), mesh.getDimension());
}
specific_body_force.resize(b.size());
std::copy_n(b.data(), b.size(), specific_body_force.data());
}
auto media_map =
MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
DBUG(
"Check the media properties of ThermoRichardsFlow process "
"...");
checkMPLProperties(media);
DBUG("Media properties verified.");
//! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__mass_lumping}
bool const mass_lumping = config.getConfigParameter<bool>("mass_lumping", false);
std::unique_ptr<SimplifiedElasticityModel> simplified_elasticity =
createElasticityModel(config);
ThermoRichardsFlowProcessData process_data{
std::move(media_map), std::move(specific_body_force),
mass_lumping, std::move(simplified_elasticity)};
SecondaryVariableCollection secondary_variables;
ProcessLib::createSecondaryVariables(config, secondary_variables);
return std::make_unique<ThermoRichardsFlowProcess>(
std::move(name), mesh, std::move(jacobian_assembler), parameters,
integration_order, std::move(process_variables),
std::move(process_data), std::move(secondary_variables),
use_monolithic_scheme);
}
} // namespace ThermoRichardsFlow
} // namespace ProcessLib
/**
* \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 <map>
#include <memory>
#include <vector>
namespace BaseLib
{
class ConfigTree;
}
namespace MeshLib
{
class Mesh;
}
namespace MaterialPropertyLib
{
class Medium;
}
namespace ParameterLib
{
struct ParameterBase;
} // namespace ParameterLib
namespace ProcessLib
{
class AbstractJacobianAssembler;
class Process;
class ProcessVariable;
} // namespace ProcessLib
namespace ProcessLib
{
namespace ThermoRichardsFlow
{
std::unique_ptr<Process> createThermoRichardsFlowProcess(
std::string name,
MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
unsigned const integration_order,
BaseLib::ConfigTree const& config,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
} // namespace ThermoRichardsFlow
} // namespace ProcessLib
/**
* \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
/**
* \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
{
template <typename ShapeMatricesType>
struct IntegrationPointData final
{
typename ShapeMatricesType::NodalRowVectorType N;
typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
typename ShapeMatricesType::GlobalDimVectorType v_darcy;
double saturation = std::numeric_limits<double>::quiet_NaN();
double saturation_prev = std::numeric_limits<double>::quiet_NaN();
double porosity = std::numeric_limits<double>::quiet_NaN();
double porosity_prev = std::numeric_limits<double>::quiet_NaN();
double transport_porosity = std::numeric_limits<double>::quiet_NaN();
double transport_porosity_prev = std::numeric_limits<double>::quiet_NaN();
double dry_density_solid = std::numeric_limits<double>::quiet_NaN();
double dry_density_pellet_saturated =
std::numeric_limits<double>::quiet_NaN();
double dry_density_pellet_unsaturated =
std::numeric_limits<double>::quiet_NaN();
double integration_weight;
void pushBackState()
{
saturation_prev = saturation;
porosity_prev = porosity;
transport_porosity_prev = transport_porosity;
}
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
} // namespace ThermoRichardsFlow
} // namespace ProcessLib
/**
* \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 <vector>
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "ProcessLib/LocalAssemblerInterface.h"
namespace ProcessLib
{
namespace ThermoRichardsFlow
{
struct LocalAssemblerInterface : public ProcessLib