OGS
MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim > Class Template Referencefinal

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 >:
[legend]
Collaboration diagram for MaterialLib::Fracture::Coulomb::Coulomb< DisplacementDim >:
[legend]

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.

94 : _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
95 _penalty_aperture_cutoff(penalty_aperture_cutoff),
96 _tension_cutoff(tension_cutoff),
97 _mp(std::move(material_properties))
98 {
99 }
NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters
Definition Coulomb.h:135

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
tcurrent time
xcurrent position in space
aperture0initial fracture's aperture
sigma0initial stress
w_prevfracture displacement at previous time step
wfracture displacement at current time step
sigma_prevstress at previous time step
sigmastress at current time step
Keptangent matrix for stress and fracture displacements
material_state_variablesmaterial state variables

Definition at line 50 of file Coulomb.cpp.

68{
69 assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
70 &material_state_variables) != nullptr);
71
72 StateVariables<DisplacementDim>& state =
73 static_cast<StateVariables<DisplacementDim>&>(material_state_variables);
74
75 MaterialPropertyValues const mat(_mp, t, x);
76
77 const int index_ns = DisplacementDim - 1;
78 double const aperture = w[index_ns] + aperture0;
79
80 Eigen::MatrixXd Ke;
81 { // Elastic tangent stiffness
82 Ke = Eigen::MatrixXd::Zero(DisplacementDim, DisplacementDim);
83 for (int i = 0; i < index_ns; i++)
84 {
85 Ke(i, i) = mat.Ks;
86 }
87
88 Ke(index_ns, index_ns) = mat.Kn;
89 }
90
91 // Total plastic aperture compression
92 // NOTE: Initial condition sigma0 seems to be associated with an initial
93 // condition of the w0 = 0. Therefore the initial state is not associated
94 // with a plastic aperture change.
95 { // Exact elastic predictor
96 sigma.noalias() = Ke * (w - w_prev);
97
98 sigma.coeffRef(index_ns) *=
100 sigma.noalias() += sigma_prev;
101 }
102
103 // correction for an opening fracture
104 if (_tension_cutoff && sigma[DisplacementDim - 1] >= 0)
105 {
106 Kep.setZero();
107 sigma.setZero();
108 state.w_p = w;
109 material_state_variables.setTensileStress(true);
110 return;
111 }
112
113 auto yield_function = [&mat](Eigen::VectorXd const& s)
114 {
115 double const sigma_n = s[DisplacementDim - 1];
116 Eigen::VectorXd const sigma_s = s.head(DisplacementDim - 1);
117 double const mag_tau = sigma_s.norm(); // magnitude
118 return mag_tau + sigma_n * std::tan(mat.phi) - mat.c;
119 };
120
121 { // Exit if still in elastic range by checking the shear yield function.
122 double const Fs = yield_function(sigma);
123 material_state_variables.setShearYieldFunctionValue(Fs);
124 if (Fs < .0)
125 {
126 Kep = Ke;
127 Kep(index_ns, index_ns) *= logPenaltyDerivative(
128 aperture0, aperture, _penalty_aperture_cutoff);
129 return;
130 }
131 }
132
133 auto yield_function_derivative = [&mat](Eigen::VectorXd const& s)
134 {
135 Eigen::Matrix<double, DisplacementDim, 1> dFs_dS;
136 dFs_dS.template head<DisplacementDim - 1>().noalias() =
137 s.template head<DisplacementDim - 1>().normalized();
138 dFs_dS.coeffRef(DisplacementDim - 1) = std::tan(mat.phi);
139 return dFs_dS;
140 };
141
142 // plastic potential function: Qs = |tau| + Sn * tan da
143 auto plastic_potential_derivative = [&mat](Eigen::VectorXd const& s)
144 {
145 Eigen::Matrix<double, DisplacementDim, 1> dQs_dS;
146 dQs_dS.template head<DisplacementDim - 1>().noalias() =
147 s.template head<DisplacementDim - 1>().normalized();
148 dQs_dS.coeffRef(DisplacementDim - 1) = std::tan(mat.psi);
149 return dQs_dS;
150 };
151
152 { // Newton
153
154 Eigen::FullPivLU<Eigen::Matrix<double, 1, 1, Eigen::RowMajor>>
155 linear_solver;
156 using ResidualVectorType = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
157 using JacobianMatrix = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
158
159 JacobianMatrix jacobian;
160 ResidualVectorType solution;
161 solution << 0;
162
163 auto const update_residual = [&](ResidualVectorType& residual)
164 { residual[0] = yield_function(sigma); };
165
166 auto const update_jacobian = [&](JacobianMatrix& jacobian)
167 {
168 jacobian(0, 0) = -yield_function_derivative(sigma).transpose() *
169 Ke * plastic_potential_derivative(sigma);
170 };
171
172 auto const update_solution = [&](ResidualVectorType const& increment)
173 {
174 solution += increment;
175 /*DBUG("analytical = {:g}",
176 Fs / (mat.Ks + mat.Kn * std::tan(mat.psi) * std::tan(mat.phi)))
177 */
178 state.w_p = state.w_p_prev +
179 solution[0] * plastic_potential_derivative(sigma);
180
181 sigma.noalias() = Ke * (w - w_prev - state.w_p + state.w_p_prev);
182
183 sigma.coeffRef(index_ns) *= logPenaltyDerivative(
184 aperture0, aperture, _penalty_aperture_cutoff);
185 sigma.noalias() += sigma_prev;
186 };
187
188 auto newton_solver = NumLib::NewtonRaphson(
189 linear_solver, update_jacobian, update_residual, update_solution,
191
192 auto const success_iterations = newton_solver.solve(jacobian);
193
194 if (!success_iterations)
195 {
197 "FractureModel/Coulomb local nonlinear solver didn't "
198 "converge.");
199 }
200
201 // Solution containing lambda is not needed; w_p and sigma already
202 // up to date.
203 }
204
205 { // Update material state shear yield function value.
206 double const Fs = yield_function(sigma);
207 material_state_variables.setShearYieldFunctionValue(Fs);
208 }
209
210 Ke(index_ns, index_ns) *=
212 Eigen::RowVectorXd const A = yield_function_derivative(sigma).transpose() *
213 Ke /
214 (yield_function_derivative(sigma).transpose() *
215 Ke * plastic_potential_derivative(sigma));
216 Kep = Ke - Ke * plastic_potential_derivative(sigma) * A;
217}
double logPenaltyDerivative(double const aperture0, double const aperture, double const aperture_cutoff)
Definition LogPenalty.h:23

References MaterialLib::Fracture::logPenaltyDerivative(), MaterialLib::Fracture::FractureModelBase< DisplacementDim >::MaterialStateVariables::setShearYieldFunctionValue(), MaterialLib::Fracture::FractureModelBase< DisplacementDim >::MaterialStateVariables::setTensileStress(), 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.

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

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: