MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim > Class Template Reference

Detailed Description

template<int DisplacementDim>
class MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >

Definition at line 45 of file Coulomb.h.

#include <Coulomb.h>

Inheritance diagram for MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >:

Collaboration diagram for MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >:

Classes
struct	MaterialProperties
	Variables specific to the material model. More...

Public Member Functions
std::unique_ptr< typename FractureModelBase< DisplacementDim >::MaterialStateVariables >	createMaterialStateVariables () override
	Coulomb (NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters, double const penalty_aperture_cutoff, bool const tension_cutoff, MaterialProperties material_properties)
void	computeConstitutiveRelation (double const t, ParameterLib::SpatialPosition const &x, double const aperture0, Eigen::Ref< Eigen::VectorXd const > sigma0, Eigen::Ref< Eigen::VectorXd const > w_prev, Eigen::Ref< Eigen::VectorXd const > w, Eigen::Ref< Eigen::VectorXd const > sigma_prev, Eigen::Ref< Eigen::VectorXd > sigma, Eigen::Ref< Eigen::MatrixXd > Kep, typename FractureModelBase< DisplacementDim >::MaterialStateVariables &material_state_variables) override
Public Member Functions inherited from MaterialLib::Fracture::FractureModelBase< DisplacementDim >
virtual	~FractureModelBase ()=default
virtual void	computeConstitutiveRelation (double const t, ParameterLib::SpatialPosition const &x, double const aperture0, Eigen::Ref< Eigen::VectorXd const > sigma0, Eigen::Ref< Eigen::VectorXd const > w_prev, Eigen::Ref< Eigen::VectorXd const > w, Eigen::Ref< Eigen::VectorXd const > sigma_prev, Eigen::Ref< Eigen::VectorXd > sigma, Eigen::Ref< Eigen::MatrixXd > C, MaterialStateVariables &material_state_variables)=0

Private Attributes
NumLib::NewtonRaphsonSolverParameters const	_nonlinear_solver_parameters
double const	_penalty_aperture_cutoff
bool const	_tension_cutoff
MaterialProperties	_mp

Constructor & Destructor Documentation

Coulomb()

template<int DisplacementDim>

MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >::Coulomb ( NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters, double const penalty_aperture_cutoff, bool const tension_cutoff, MaterialProperties material_properties )

inlineexplicit

Definition at line 89 of file Coulomb.h.

        : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
          _penalty_aperture_cutoff(penalty_aperture_cutoff),
          _tension_cutoff(tension_cutoff),
          _mp(std::move(material_properties))
    {
    }

Member Function Documentation

computeConstitutiveRelation()

template<int DisplacementDim>

void MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >::computeConstitutiveRelation ( double const t, ParameterLib::SpatialPosition const & x, double const aperture0, Eigen::Ref< Eigen::VectorXd const > sigma0, Eigen::Ref< Eigen::VectorXd const > w_prev, Eigen::Ref< Eigen::VectorXd const > w, Eigen::Ref< Eigen::VectorXd const > sigma_prev, Eigen::Ref< Eigen::VectorXd > sigma, Eigen::Ref< Eigen::MatrixXd > Kep, typename FractureModelBase< DisplacementDim >::MaterialStateVariables & material_state_variables )

override

Computation of the constitutive relation for the Mohr-Coulomb model.

Parameters

t	current time
x	current position in space
aperture0	initial fracture's aperture
sigma0	initial stress
w_prev	fracture displacement at previous time step
w	fracture displacement at current time step
sigma_prev	stress at previous time step
sigma	stress at current time step
Kep	tangent matrix for stress and fracture displacements
material_state_variables	material state variables

Definition at line 50 of file Coulomb.cpp.

{
    assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
               &material_state_variables) != nullptr);
 
    StateVariables<DisplacementDim>& state =
        static_cast<StateVariables<DisplacementDim>&>(material_state_variables);
 
    MaterialPropertyValues const mat(_mp, t, x);
 
    const int index_ns = DisplacementDim - 1;
    double const aperture = w[index_ns] + aperture0;
 
    Eigen::MatrixXd Ke;
    {  // Elastic tangent stiffness
        Ke = Eigen::MatrixXd::Zero(DisplacementDim, DisplacementDim);
        for (int i = 0; i < index_ns; i++)
        {
            Ke(i, i) = mat.Ks;
        }
 
        Ke(index_ns, index_ns) = mat.Kn;
    }
 
    // Total plastic aperture compression
    // NOTE: Initial condition sigma0 seems to be associated with an initial
    // condition of the w0 = 0. Therefore the initial state is not associated
    // with a plastic aperture change.
    {  // Exact elastic predictor
        sigma.noalias() = Ke * (w - w_prev);
 
        sigma.coeffRef(index_ns) *=
            logPenaltyDerivative(aperture0, aperture, _penalty_aperture_cutoff);
        sigma.noalias() += sigma_prev;
    }
 
    // correction for an opening fracture
    if (_tension_cutoff && sigma[DisplacementDim - 1] >= 0)
    {
        Kep.setZero();
        sigma.setZero();
        state.w_p = w;
        material_state_variables.setTensileStress(true);
        return;
    }
 
    auto yield_function = [&mat](Eigen::VectorXd const& s)
    {
        double const sigma_n = s[DisplacementDim - 1];
        Eigen::VectorXd const sigma_s = s.head(DisplacementDim - 1);
        double const mag_tau = sigma_s.norm();  // magnitude
        return mag_tau + sigma_n * std::tan(mat.phi) - mat.c;
    };
 
    {  // Exit if still in elastic range by checking the shear yield function.
        double const Fs = yield_function(sigma);
        material_state_variables.setShearYieldFunctionValue(Fs);
        if (Fs < .0)
        {
            Kep = Ke;
            Kep(index_ns, index_ns) *= logPenaltyDerivative(
                aperture0, aperture, _penalty_aperture_cutoff);
            return;
        }
    }
 
    auto yield_function_derivative = [&mat](Eigen::VectorXd const& s)
    {
        Eigen::Matrix<double, DisplacementDim, 1> dFs_dS;
        dFs_dS.template head<DisplacementDim - 1>().noalias() =
            s.template head<DisplacementDim - 1>().normalized();
        dFs_dS.coeffRef(DisplacementDim - 1) = std::tan(mat.phi);
        return dFs_dS;
    };
 
    // plastic potential function: Qs = |tau| + Sn * tan da
    auto plastic_potential_derivative = [&mat](Eigen::VectorXd const& s)
    {
        Eigen::Matrix<double, DisplacementDim, 1> dQs_dS;
        dQs_dS.template head<DisplacementDim - 1>().noalias() =
            s.template head<DisplacementDim - 1>().normalized();
        dQs_dS.coeffRef(DisplacementDim - 1) = std::tan(mat.psi);
        return dQs_dS;
    };
 
    {  // Newton
 
        Eigen::FullPivLU<Eigen::Matrix<double, 1, 1, Eigen::RowMajor>>
            linear_solver;
        using ResidualVectorType = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
        using JacobianMatrix = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
 
        JacobianMatrix jacobian;
        ResidualVectorType solution;
        solution << 0;
 
        auto const update_residual = [&](ResidualVectorType& residual)
        { residual[0] = yield_function(sigma); };
 
        auto const update_jacobian = [&](JacobianMatrix& jacobian)
        {
            jacobian(0, 0) = -yield_function_derivative(sigma).transpose() *
                             Ke * plastic_potential_derivative(sigma);
        };
 
        auto const update_solution = [&](ResidualVectorType const& increment)
        {
            solution += increment;
            /*DBUG("analytical = {:g}",
                 Fs / (mat.Ks + mat.Kn * std::tan(mat.psi) * std::tan(mat.phi)))
                 */
            state.w_p = state.w_p_prev +
                        solution[0] * plastic_potential_derivative(sigma);
 
            sigma.noalias() = Ke * (w - w_prev - state.w_p + state.w_p_prev);
 
            sigma.coeffRef(index_ns) *= logPenaltyDerivative(
                aperture0, aperture, _penalty_aperture_cutoff);
            sigma.noalias() += sigma_prev;
        };
 
        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)
        {
            throw NumLib::AssemblyException(
                "FractureModel/Coulomb local nonlinear solver didn't "
                "converge.");
        }
 
        // Solution containing lambda is not needed; w_p and sigma already
        // up to date.
    }
 
    {  // Update material state shear yield function value.
        double const Fs = yield_function(sigma);
        material_state_variables.setShearYieldFunctionValue(Fs);
    }
 
    Ke(index_ns, index_ns) *=
        logPenaltyDerivative(aperture0, aperture, _penalty_aperture_cutoff);
    Eigen::RowVectorXd const A = yield_function_derivative(sigma).transpose() *
                                 Ke /
                                 (yield_function_derivative(sigma).transpose() *
                                  Ke * plastic_potential_derivative(sigma));
    Kep = Ke - Ke * plastic_potential_derivative(sigma) * A;
}

References MaterialLib::Fracture::logPenaltyDerivative(), MaterialLib::Fracture::FractureModelBase< DisplacementDim >::MaterialStateVariables::setShearYieldFunctionValue(), MaterialLib::Fracture::FractureModelBase< DisplacementDim >::MaterialStateVariables::setTensileStress(), NumLib::NewtonRaphson< LinearSolver, JacobianMatrix, JacobianMatrixUpdate, ResidualVector, ResidualUpdate, SolutionUpdate >::solve(), MaterialLib::Fracture::Coulomb::StateVariables< DisplacementDim >::w_p, and MaterialLib::Fracture::Coulomb::StateVariables< DisplacementDim >::w_p_prev.

createMaterialStateVariables()

template<int DisplacementDim>

std::unique_ptr< typename FractureModelBase< DisplacementDim >::MaterialStateVariables > MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >::createMaterialStateVariables ( )

inlineoverridevirtual

Polymorphic creator for MaterialStateVariables objects specific for a material model.

Implements MaterialLib::Fracture::FractureModelBase< DisplacementDim >.

Definition at line 50 of file Coulomb.h.

    {
        return std::make_unique<StateVariables<DisplacementDim>>();
    }

Member Data Documentation

_mp

template<int DisplacementDim>

MaterialProperties MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >::_mp

private

Definition at line 147 of file Coulomb.h.

_nonlinear_solver_parameters

template<int DisplacementDim>

NumLib::NewtonRaphsonSolverParameters const MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >::_nonlinear_solver_parameters

private

Definition at line 135 of file Coulomb.h.

_penalty_aperture_cutoff

template<int DisplacementDim>

double const MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >::_penalty_aperture_cutoff

private

Compressive normal displacements above this value will not enter the computation of the normal stiffness modulus of the fracture.

Note

Setting this to the initial aperture value allows negative apertures.

Definition at line 141 of file Coulomb.h.

_tension_cutoff

template<int DisplacementDim>

bool const MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >::_tension_cutoff

private

If set no resistance to open the fracture over the initial aperture is opposed.

Definition at line 145 of file Coulomb.h.

---

The documentation for this class was generated from the following files:

• MaterialLib/FractureModels/Coulomb.h
• MaterialLib/FractureModels/Coulomb.cpp