https://gitlab.opengeosys.org/ogs/ogs.git
Tip revision: 66570a258a1844aab06ab5460511e34480a5026a authored by joergbuchwald on 08 November 2021, 07:59:42 UTC
sounds better
sounds better
Tip revision: 66570a2
Lubby2.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 "MathLib/KelvinVector.h"
#include "NumLib/NewtonRaphson.h"
#include "ParameterLib/Parameter.h"
#include "MechanicsBase.h"
namespace MaterialLib
{
namespace Solids
{
namespace Lubby2
{
//
// Variables specific to the material model.
//
struct Lubby2MaterialProperties
{
using P = ParameterLib::Parameter<double>;
Lubby2MaterialProperties(P const& GK0_,
P const& GM0_,
P const& KM0_,
P const& etaK0_,
P const& etaM0_,
P const& mK_,
P const& mvK_,
P const& mvM_)
: GK0(GK0_),
GM0(GM0_),
KM0(KM0_),
etaK0(etaK0_),
etaM0(etaM0_),
mK(mK_),
mvK(mvK_),
mvM(mvM_)
{
}
// basic material parameters
P const& GK0;
P const& GM0;
P const& KM0;
P const& etaK0;
P const& etaM0;
P const& mK;
P const& mvK;
P const& mvM;
};
namespace detail
{
template <int DisplacementDim>
struct LocalLubby2Properties
{
LocalLubby2Properties(double const t,
ParameterLib::SpatialPosition const& x,
Lubby2MaterialProperties const& mp)
: GM0(mp.GM0(t, x)[0]),
KM0(mp.KM0(t, x)[0]),
GK0(mp.GK0(t, x)[0]),
etaK0(mp.etaK0(t, x)[0]),
etaM0(mp.etaM0(t, x)[0]),
mK(mp.mK(t, x)[0]),
mvK(mp.mvK(t, x)[0]),
mvM(mp.mvM(t, x)[0])
{
}
void update(double const s_eff)
{
double const GM0_s_eff = GM0 * s_eff;
GK = GK0 * std::exp(mK * GM0_s_eff);
etaK = etaK0 * std::exp(mvK * GM0_s_eff);
etaM = etaM0 * std::exp(mvM * GM0_s_eff);
}
double const GM0;
double const KM0;
double const GK0;
double const etaK0;
double const etaM0;
double const mK;
double const mvK;
double const mvM;
// Solution dependent values.
double GK = std::numeric_limits<double>::quiet_NaN();
double etaK = std::numeric_limits<double>::quiet_NaN();
double etaM = std::numeric_limits<double>::quiet_NaN();
};
} // namespace detail
template <int DisplacementDim>
class Lubby2 final : public MechanicsBase<DisplacementDim>
{
public:
struct MaterialStateVariables
: public MechanicsBase<DisplacementDim>::MaterialStateVariables
{
MaterialStateVariables()
{
// Previous time step values are not initialized and are set later.
eps_K_t.resize(KelvinVectorSize);
eps_M_t.resize(KelvinVectorSize);
// Initialize current time step values
eps_K_j.setZero(KelvinVectorSize);
eps_M_j.setZero(KelvinVectorSize);
}
void setInitialConditions()
{
eps_K_j = eps_K_t;
eps_M_j = eps_M_t;
}
void pushBackState() override
{
eps_K_t = eps_K_j;
eps_M_t = eps_M_j;
}
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
/// Deviatoric strain in the viscous kelvin element during the current
/// iteration
KelvinVector eps_K_t;
KelvinVector eps_K_j;
/// Deviatoric strain in the viscous maxwell element during the current
/// iteration
KelvinVector eps_M_t;
KelvinVector eps_M_j;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
std::unique_ptr<
typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
createMaterialStateVariables() const override
{
return std::unique_ptr<
typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
new MaterialStateVariables};
}
public:
static int const KelvinVectorSize =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
static int const JacobianResidualSize =
3 * KelvinVectorSize; // Three is the number of components in the
// jacobian/residual, not the space dimension.
using ResidualVector = Eigen::Matrix<double, JacobianResidualSize, 1>;
using JacobianMatrix = Eigen::Matrix<double,
JacobianResidualSize,
JacobianResidualSize,
Eigen::RowMajor>;
public:
explicit Lubby2(
NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
Lubby2MaterialProperties& material_properties)
: _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
_mp(material_properties)
{
}
double computeFreeEnergyDensity(
double const t,
ParameterLib::SpatialPosition const& x,
double const dt,
KelvinVector const& eps,
KelvinVector const& sigma,
typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
material_state_variables) const override
{
assert(dynamic_cast<MaterialStateVariables const*>(
&material_state_variables) != nullptr);
MaterialStateVariables state(static_cast<MaterialStateVariables const&>(
material_state_variables));
auto const& eps_K = state.eps_K_j;
auto const& eps_K_prev = state.eps_K_t;
auto const& eps_M = state.eps_M_j;
auto local_lubby2_properties =
detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
// calculation of deviatoric parts
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
auto const& P_dev = Invariants::deviatoric_projection;
KelvinVector const epsd_i = P_dev * eps;
// initial guess as elastic predictor
KelvinVector sigd_j = 2.0 * (epsd_i - state.eps_M_t - state.eps_K_t);
// Calculate effective stress and update material properties
double sig_eff = Invariants::equivalentStress(sigd_j);
local_lubby2_properties.update(sig_eff);
auto const& eta_K = local_lubby2_properties.etaK;
// This is specific to the backward Euler time scheme and needs to be
// updated if other time schemes are used.
return (eps - eps_K - eps_M).dot(sigma) / 2 +
eps_K.dot(sigma - eta_K * (eps_K - eps_K_prev) / dt) / 2;
}
double getBulkModulus(double const t,
ParameterLib::SpatialPosition const& x,
KelvinMatrix const* const /*C*/) const override
{
return _mp.KM0(t, x)[0];
}
std::optional<std::tuple<KelvinVector,
std::unique_ptr<typename MechanicsBase<
DisplacementDim>::MaterialStateVariables>,
KelvinMatrix>>
integrateStress(
MaterialPropertyLib::VariableArray const& variable_array_prev,
MaterialPropertyLib::VariableArray const& variable_array,
double const t, ParameterLib::SpatialPosition const& x, double const dt,
typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
material_state_variables) const override;
private:
/// Calculates the 18x1 residual vector.
void calculateResidualBurgers(
double const dt,
const KelvinVector& strain_curr,
const KelvinVector& strain_t,
const KelvinVector& stress_curr,
const KelvinVector& stress_t,
const KelvinVector& strain_Kel_curr,
const KelvinVector& strain_Kel_t,
const KelvinVector& strain_Max_curr,
const KelvinVector& strain_Max_t,
ResidualVector& res,
detail::LocalLubby2Properties<DisplacementDim> const& properties) const;
/// Calculates the 18x18 Jacobian.
void calculateJacobianBurgers(
double const t,
ParameterLib::SpatialPosition const& x,
double const dt,
JacobianMatrix& Jac,
double s_eff,
const KelvinVector& sig_i,
const KelvinVector& eps_K_i,
detail::LocalLubby2Properties<DisplacementDim> const& properties) const;
private:
NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters;
Lubby2MaterialProperties _mp;
};
extern template class Lubby2<2>;
extern template class Lubby2<3>;
} // namespace Lubby2
} // namespace Solids
} // namespace MaterialLib