Skip to content
Snippets Groups Projects
Commit 4b02c231 authored by Dmitri Naumov's avatar Dmitri Naumov
Browse files

Merge branch 'RM_DoublePorosityModel' into 'master' 

RM; double porosity model

See merge request ogs/ogs!3405
parents ab8d4bbd 51933b58
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MaterialLib/MPL/PropertyType.h
+ 7
− 0
View file @ 4b02c231
......@@ -77,6 +77,8 @@ enum PropertyType : int
/// specify retardation factor used in component transport process.
retardation_factor,
saturation,
/// capillary pressure saturation relationship for microstructure.
saturation_micro,
specific_heat_capacity,
storage,
swelling_stress_rate,
......@@ -257,6 +259,10 @@ inline PropertyType convertStringToProperty(std::string const& inString)
{
return PropertyType::saturation;
}
if (boost::iequals(inString, "saturation_micro"))
{
return PropertyType::saturation_micro;
}
if (boost::iequals(inString, "specific_heat_capacity"))
{
return PropertyType::specific_heat_capacity;
......@@ -358,6 +364,7 @@ static const std::array<std::string, PropertyType::number_of_properties>
"residual_liquid_saturation",
"retardation_factor",
"saturation",
"saturation_micro",
"specific_heat_capacity",
"storage",
"swelling_stress_rate",
......
ProcessLib/RichardsMechanics/ComputeMicroPorosity.h 0 → 100644
+ 211
− 0
View file @ 4b02c231
/**
* \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 <cassert>
#include <ostream>
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
#include "MathLib/KelvinVector.h"
#include "NumLib/NewtonRaphson.h"
namespace ProcessLib::RichardsMechanics
{
template <int DisplacementDim>
struct MicroPorosityStateSpace
{
double phi_m;
double e_sw;
double p_L_m;
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> sigma_sw;
MicroPorosityStateSpace& operator += (MicroPorosityStateSpace const& state)
{
phi_m += state.phi_m;
e_sw += state.e_sw;
p_L_m += state.p_L_m;
sigma_sw += state.sigma_sw;
return *this;
}
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
template <int DisplacementDim>
std::ostream& operator<<(std::ostream& os,
MicroPorosityStateSpace<DisplacementDim> const& state)
{
return os << "phi_m: " << state.phi_m << ", "
<< "e_sw: " << state.e_sw << ", "
<< "p_L_m: " << state.p_L_m << ", "
<< "sigma_sw: " << state.sigma_sw.transpose();
}
template <int DisplacementDim>
MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const&
I_2_C_el_inverse,
double const rho_LR_m, // for simplification equal to rho_LR_M
double const mu_LR, double const alpha_bar, double const alpha_B,
double const phi_M, double const p_L, double const p_L_m_prev,
MaterialPropertyLib::VariableArray const& /*variables_prev*/,
double const S_L_m_prev, double const phi_m_prev,
ParameterLib::SpatialPosition const pos, double const t, double const dt,
MaterialPropertyLib::Property const& saturation_micro,
MaterialPropertyLib::Property const& swelling_stress_rate,
NumLib::NewtonRaphsonSolverParameters const& nonlinear_solver_parameters)
{
namespace MPL = MaterialPropertyLib;
static constexpr int kelvin_vector_size =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
// phi_m, e_sw, p_L_m, sigma_sw
static constexpr int nls_size = 1 + 1 + 1 + kelvin_vector_size;
static constexpr int i_phi_m = 0;
static constexpr int i_e_sw = 1;
static constexpr int i_p_L_m = 2;
static constexpr int i_sigma_sw = 3;
using ResidualVectorType = Eigen::Matrix<double, nls_size, 1>;
using JacobianMatrix =
Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>;
Eigen::FullPivLU<Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>>
linear_solver;
JacobianMatrix jacobian;
// Agglomerated solution vector construction.
ResidualVectorType solution = ResidualVectorType::Zero();
if (p_L >= 0 && p_L_m_prev >= 0)
{
return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
solution.template segment<kelvin_vector_size>(i_sigma_sw)};
}
auto const update_residual = [&](ResidualVectorType& residual) {
double const delta_phi_m = solution[i_phi_m];
double const delta_e_sw = solution[i_e_sw];
auto const& delta_sigma_sw =
solution.template segment<kelvin_vector_size>(i_sigma_sw);
double const delta_p_L_m = solution[i_p_L_m];
double const phi_m = phi_m_prev + delta_phi_m;
double const p_L_m = p_L_m_prev + delta_p_L_m;
MPL::VariableArray variables_prev;
variables_prev[static_cast<int>(MPL::Variable::capillary_pressure)] =
-p_L_m_prev;
variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
S_L_m_prev;
MPL::VariableArray variables;
variables[static_cast<int>(MPL::Variable::capillary_pressure)] = -p_L_m;
double const S_L_m =
saturation_micro.template value<double>(variables, pos, t, dt);
variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L_m;
double const delta_S_L_m = S_L_m - S_L_m_prev;
auto const sigma_sw_dot =
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
swelling_stress_rate.template value<Eigen::Matrix3d>(
variables, variables_prev, pos, t, dt));
residual[i_phi_m] = delta_phi_m - (alpha_B - phi_M) * delta_e_sw;
residual[i_e_sw] = delta_e_sw + I_2_C_el_inverse.dot(delta_sigma_sw);
residual.template segment<kelvin_vector_size>(i_sigma_sw).noalias() =
delta_sigma_sw - sigma_sw_dot * dt;
residual[i_p_L_m] =
rho_LR_m *
(phi_m * delta_S_L_m - (alpha_B - phi_M) * S_L_m * delta_e_sw) +
phi_m * S_L_m * rho_LR_m * delta_e_sw -
alpha_bar / mu_LR * (p_L - p_L_m) * dt;
};
auto const update_jacobian = [&](JacobianMatrix& jacobian) {
jacobian = JacobianMatrix::Identity();
double const delta_phi_m = solution[i_phi_m];
double const delta_e_sw = solution[i_e_sw];
double const delta_p_L_m = solution[i_p_L_m];
double const phi_m = phi_m_prev + delta_phi_m;
double const p_L_m = p_L_m_prev + delta_p_L_m;
MPL::VariableArray variables_prev;
variables_prev[static_cast<int>(MPL::Variable::capillary_pressure)] =
-p_L_m_prev;
MPL::VariableArray variables;
variables[static_cast<int>(MPL::Variable::capillary_pressure)] = -p_L_m;
double const S_L_m =
saturation_micro.template value<double>(variables, pos, t, dt);
variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
S_L_m_prev;
variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L_m;
double const delta_S_L_m = S_L_m - S_L_m_prev;
double const dS_L_m_dp_cap_m = saturation_micro.template dValue<double>(
variables, MPL::Variable::capillary_pressure, pos, t, dt);
auto const dsigma_sw_dS_L_m =
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
swelling_stress_rate.template dValue<Eigen::Matrix3d>(
variables, variables_prev, MPL::Variable::liquid_saturation,
pos, t, dt));
jacobian(i_phi_m, i_e_sw) = -(alpha_B - phi_M);
jacobian.template block<1, kelvin_vector_size>(i_e_sw, i_sigma_sw) =
I_2_C_el_inverse.transpose();
jacobian.template block<kelvin_vector_size, 1>(
i_sigma_sw, i_p_L_m) = -dsigma_sw_dS_L_m * dS_L_m_dp_cap_m;
jacobian(i_p_L_m, i_phi_m) =
rho_LR_m * (delta_S_L_m + phi_m * delta_e_sw);
jacobian(i_p_L_m, i_e_sw) = -rho_LR_m * S_L_m * (alpha_B - phi_M - phi_m);
jacobian(i_p_L_m, i_p_L_m) =
alpha_bar / mu_LR * dt -
rho_LR_m * (phi_m - (alpha_B - phi_M - phi_m) * delta_e_sw) *
dS_L_m_dp_cap_m;
};
auto const update_solution = [&](ResidualVectorType const& increment) {
solution += increment;
};
auto newton_solver =
NumLib::NewtonRaphson<decltype(linear_solver), JacobianMatrix,
decltype(update_jacobian), ResidualVectorType,
decltype(update_residual),
decltype(update_solution)>(
linear_solver, update_jacobian, update_residual, update_solution,
nonlinear_solver_parameters);
auto const success_iterations = newton_solver.solve(jacobian);
if (!success_iterations)
{
OGS_FATAL(
"Could not find solution for local double structure nonlinear "
"problem.");
}
return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
solution.template segment<kelvin_vector_size>(i_sigma_sw)};
}
} // namespace ProcessLib::RichardsMechanics
Tests/CMakeLists.txt
+ 3
− 0
View file @ 4b02c231
......@@ -43,6 +43,9 @@ endif()
if(OGS_BUILD_PROCESS_LIE)
append_source_files(TEST_SOURCES ProcessLib/LIE)
endif()
if(OGS_BUILD_PROCESS_RichardsMechanics)
append_source_files(TEST_SOURCES ProcessLib/RichardsMechanics)
endif()
if(OGS_USE_PETSC)
list(REMOVE_ITEM TEST_SOURCES NumLib/TestSerialLinearSolver.cpp)
......
Tests/ProcessLib/RichardsMechanics/MicroporosityComputation.cpp 0 → 100644
+ 262
− 0
View file @ 4b02c231
/**
* \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 <gtest/gtest.h>
#include <Eigen/Eigen>
#include <cmath>
#include <memory>
#include "MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.h"
#include "MaterialLib/MPL/Properties/SaturationDependentSwelling.h"
#include "MaterialLib/MPL/Property.h"
#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
#include "MaterialLib/SolidModels/MechanicsBase.h"
#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
#include "ParameterLib/ConstantParameter.h"
#include "ProcessLib/RichardsMechanics/ComputeMicroPorosity.h"
using namespace ProcessLib::RichardsMechanics;
namespace MPL = MaterialPropertyLib;
namespace
{
template <int DisplacementDim>
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
computeElasticTangentStiffness(
MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material,
double const t, ParameterLib::SpatialPosition const& x_position,
double const dt, double const temperature)
{
MPL::VariableArray variable_array;
MPL::VariableArray variable_array_prev;
auto const null_state = solid_material.createMaterialStateVariables();
using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
variable_array[static_cast<int>(MPL::Variable::stress)].emplace<KV>(
KV::Zero());
variable_array[static_cast<int>(MPL::Variable::mechanical_strain)]
.emplace<KV>(KV::Zero());
variable_array[static_cast<int>(MPL::Variable::temperature)]
.emplace<double>(temperature);
variable_array_prev[static_cast<int>(MPL::Variable::stress)].emplace<KV>(
KV::Zero());
variable_array_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
.emplace<KV>(KV::Zero());
variable_array_prev[static_cast<int>(MPL::Variable::temperature)]
.emplace<double>(temperature);
auto&& solution = solid_material.integrateStress(
variable_array_prev, variable_array, t, x_position, dt, *null_state);
if (!solution)
{
OGS_FATAL("Computation of elastic tangent stiffness failed.");
}
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C =
std::move(std::get<2>(*solution));
return C;
}
} // namespace
TEST(RichardsMechanics, computeMicroPorosity)
{
static constexpr auto eps = 2e-14;
constexpr int DisplacementDim = 2;
static const NumLib::NewtonRaphsonSolverParameters
nonlinear_solver_parameters{1000, 1e-8, 1e-15};
//
// Create properties.
//
// Saturation
constexpr double residual_liquid_saturation = 0;
constexpr double residual_gas_saturation = 0;
constexpr double exponent = 0.4;
constexpr double p_b = 300e6;
MPL::Property const& saturation_micro = MPL::SaturationVanGenuchten{
"saturation_micro", residual_liquid_saturation, residual_gas_saturation,
exponent, p_b};
// Swelling
constexpr std::array swelling_pressures{1e7, 1e7, 1e7};
constexpr std::array swelling_exponents{1.2, 1.2, 1.2};
constexpr double lower_saturation_limit = 0;
constexpr double upper_saturation_limit = 1;
MPL::Property const& swelling_stress_rate =
MPL::SaturationDependentSwelling{
"swelling_stress_rate", swelling_pressures, swelling_exponents,
lower_saturation_limit, upper_saturation_limit, nullptr};
// LinearElastic
ParameterLib::ConstantParameter<double> E{"YoungsModulus", 150e6};
ParameterLib::ConstantParameter<double> nu{"PoissonsRatio", 0.25};
auto const solid_material =
MaterialLib::Solids::LinearElasticIsotropic<DisplacementDim>{{E, nu}};
//
// Populate constants
//
ParameterLib::SpatialPosition const pos;
constexpr double t0 = 0;
constexpr double t_max = 200e3;
double const dt = 100;
double const T = 293.15;
auto const& identity2 = MathLib::KelvinVector::Invariants<
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value>::
identity2;
auto const C_el =
computeElasticTangentStiffness(solid_material, t0, pos, dt, T);
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const&
I_2_C_el_inverse = identity2.transpose() * C_el.inverse();
double const rho_LR_m = 1e3;
double const mu_LR = 1e-3;
double const alpha_bar = 1e-14;
double const alpha_B = 1;
double const phi_M = 0.45;
auto saturation = [&](double const p_L) {
MPL::VariableArray v;
v[static_cast<int>(MPL::Variable::capillary_pressure)] = -p_L;
return saturation_micro.template value<double>(v, pos, t0, dt);
};
//
// Drive the computation by macro pressure. Saturation, steady, and
// desaturation.
//
// The time coordinates are simultaneously also the test points.
std::vector const t_coords{
// saturation is quick to rich maximum,
0 * t_max, 0.05 * t_max, 0.2 * t_max, 0.3 * t_max,
// desaturation takes much longer.
0.35 * t_max, 0.5 * t_max, t_max};
std::vector const pressure_M_values{// saturation
-200e6, -200e6, 0., 0.,
// desaturation
0., -200e6, -200e6};
MathLib::PiecewiseLinearInterpolation const pressure_M{t_coords,
pressure_M_values};
//
// Initial state
//
MicroPorosityStateSpace<DisplacementDim> state_prev = {
0.45 - 0.25, // phi_m
0, // e_sw
-200e6, // p_L_m
{0, 0, 0, 0} // sigma_sw
};
auto state = state_prev;
std::vector<MicroPorosityStateSpace<DisplacementDim>> results;
for (double t = t0 + dt; t <= t_max; t += dt)
{
double const p_L = pressure_M.getValue(t);
double const S_L_m_prev = saturation(state_prev.p_L_m);
//auto const [delta_phi_m, delta_e_sw, delta_sigma_sw, delta_p_L_m] =
auto const state_increment = computeMicroPorosity<DisplacementDim>(
I_2_C_el_inverse,
rho_LR_m, // for simplification equal to rho_LR_M
mu_LR, alpha_bar, alpha_B, phi_M, p_L, state_prev.p_L_m,
MaterialPropertyLib::VariableArray{}, S_L_m_prev, state_prev.phi_m,
pos, t, dt, saturation_micro, swelling_stress_rate,
nonlinear_solver_parameters);
// push back state
state_prev = state;
// update
state += state_increment;
if (std::find_if(begin(t_coords), end(t_coords), [&](auto const value) {
return std::abs(t - value) < eps;
}) != end(t_coords))
{
results.push_back(state);
/* Keep for possible result updates
std::cout << std::setprecision(17)
<< "MicroPorosityStateSpace<DisplacementDim>{ "
<< state.phi_m << ", " << state.e_sw << ", "
<< state.p_L_m << ", {" << state.sigma_sw[0] << ", "
<< state.sigma_sw[1] << ", " << state.sigma_sw[2] << ", "
<< state.sigma_sw[3] << " }}\n";
*/
}
}
//
// Expected results
//
std::vector expected_results = {
MicroPorosityStateSpace<DisplacementDim>{
0.20000000000000001, 0, -200000000, {0, 0, 0, 0}},
MicroPorosityStateSpace<DisplacementDim>{
0.20678794335787565,
0.012341715196137344,
-87023875.75987792,
{-1234171.5196137347, -1234171.5196137347, -1234171.5196137347, 0}},
MicroPorosityStateSpace<DisplacementDim>{
0.20984996711734061,
0.017909031122437243,
-4693546.2948520193,
{-1790903.1122437208, -1790903.1122437208, -1790903.1122437208, 0}},
MicroPorosityStateSpace<DisplacementDim>{
0.20987701765630837,
0.017958213920560351,
23234.984113500199,
{-1795821.3920560319, -1795821.3920560319, -1795821.3920560319, 0}},
MicroPorosityStateSpace<DisplacementDim>{
0.20475472126658667,
0.0086449477574298151,
-123250534.84079438,
{-864494.77574297995, -864494.77574297995, -864494.77574297995, 0}},
MicroPorosityStateSpace<DisplacementDim>{
0.20005534659406693,
0.0001006301710302317,
-199863745.41244847,
{-10063.017103021764, -10063.017103021764, -10063.017103021764,
0}}};
auto eps_equal = [](double const a, double const b) {
return std::abs(a - b) < eps;
};
auto const [mismatch_it_expected, mismatch_it_results] = std::mismatch(
begin(expected_results), end(expected_results), begin(results),
[&](auto const& a, auto const& b) {
EXPECT_TRUE(eps_equal(a.phi_m, b.phi_m)) << "with eps = " << eps;
EXPECT_TRUE(eps_equal(a.e_sw, b.e_sw)) << "with eps = " << eps;
EXPECT_TRUE(eps_equal(a.p_L_m / 1e9, b.p_L_m / 1e9))
<< "with eps = " << eps;
EXPECT_TRUE((a.sigma_sw - b.sigma_sw).norm() / 1e9 < eps)
<< "with eps = " << eps;
// Divide pressures and stresses by GPa to get approxmately
// O(1) numbers.
return eps_equal(a.phi_m, b.phi_m) && eps_equal(a.e_sw, b.e_sw) &&
eps_equal(a.p_L_m / 1e9, b.p_L_m / 1e9) &&
(a.sigma_sw - b.sigma_sw).norm() / 1e9 < eps;
});
ASSERT_TRUE(mismatch_it_expected == end(expected_results))
<< "Mismatch for expected result\n\t" << *mismatch_it_expected
<< "\n\tfound. Corresponding result is\n\t" << *mismatch_it_results;
ASSERT_TRUE(mismatch_it_results == end(results));
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment