35template <
int DisplacementDim>
37 std::string
const& phasefield_model)
39 if (phasefield_model ==
"AT1")
43 if (phasefield_model ==
"AT2")
47 if (phasefield_model ==
"COHESIVE")
52 "phasefield_model must be 'AT1', 'AT2' or 'COHESIVE' but '{}' "
54 phasefield_model.c_str());
63template <
int DisplacementDim>
65 std::optional<std::string>
const& softening_curve)
69 if (*softening_curve ==
"Linear")
73 if (*softening_curve ==
"Exponential")
78 "softening_curve must be 'Linear' or 'Exponential' but '{}' "
80 softening_curve->c_str());
94template <
int DisplacementDim>
96 std::string
const& energy_split_model)
98 if (energy_split_model ==
"Isotropic")
102 if (energy_split_model ==
"VolumetricDeviatoric")
106 if (energy_split_model ==
"EffectiveStress")
110 if (energy_split_model ==
"OrthoVolDev")
114 if (energy_split_model ==
"OrthoMasonry")
119 "energy_split_model must be 'Isotropic', 'VolumetricDeviatoric', "
120 "'EffectiveStress', 'OrthoVolDev' or 'OrthoMasonry' but '{}' was given",
121 energy_split_model.c_str());
130 double const ls) = 0;
132 double const ls) = 0;
134 double const ls) = 0;
142 double const )
override
144 return d_ip * d_ip * (1. - k) + k;
147 double const )
override
149 return 2. * (1. - k) * d_ip;
152 double const )
override
154 return 2. * (1. - k);
170 double const ls)
override
172 double const m1 = 4.0 *
lch / acos(-1.) / ls;
177 double const m2 = std::pow(2., 5. / 3.) - 3.;
178 double const n = 2.5;
179 return std::pow(d_ip, n) /
181 m1 * (1. - d_ip) * (1. + m2 * (1. - d_ip))) *
187 double const m2 = -0.5;
189 return std::pow(d_ip, n) /
191 m1 * (1. - d_ip) * (1. + m2 * (1. - d_ip))) *
198 double const ls)
override
200 double const m1 = 4.0 *
lch / acos(-1.) / ls;
205 double const m2 = std::pow(2., 5. / 3.) - 3.;
206 double const n = 2.5;
207 double const a1 = std::pow(d_ip, n) +
208 m1 * (1. - d_ip) * (1. + m2 * (1. - d_ip));
209 double const a2 = n * std::pow(d_ip, n - 1.) -
210 2. * m1 * m2 * (1. - d_ip) - m1;
212 (n * std::pow(d_ip, n - 1.) * a1 -
213 std::pow(d_ip, n) * a2) /
218 double const m2 = -0.5;
220 double const a1 = std::pow(d_ip, n) +
221 m1 * (1. - d_ip) * (1. + m2 * (1. - d_ip));
222 double const a2 = n * std::pow(d_ip, n - 1.) -
223 2. * m1 * m2 * (1. - d_ip) - m1;
225 (n * std::pow(d_ip, n - 1.) * a1 -
226 std::pow(d_ip, n) * a2) /
232 double const ls)
override
234 double const m1 = 4.0 *
lch / acos(-1.) / ls;
239 double const m2 = std::pow(2., 5. / 3.) - 3.;
240 double const n = 2.5;
241 double a1 = std::pow(d_ip, n) +
242 m1 * (1. - d_ip) * (1. + m2 * (1. - d_ip));
243 double a2 = n * std::pow(d_ip, n - 1.) -
244 2. * m1 * m2 * (1. - d_ip) - m1;
246 (2. * a2 * a2 * std::pow(d_ip, n) -
247 a1 * std::pow(d_ip, n) *
249 n * std::pow(d_ip, n - 2.) * (n - 1.)) -
250 2. * a1 * a2 * n * std::pow(d_ip, n - 1.) +
251 a1 * a1 * n * std::pow(d_ip, n - 2.) * (n - 1.)) /
256 double const m2 = -0.5;
258 double const a1 = std::pow(d_ip, n) +
259 m1 * (1. - d_ip) * (1. + m2 * (1. - d_ip));
260 double const a2 = n * std::pow(d_ip, n - 1.) -
261 2. * m1 * m2 * (1. - d_ip) - m1;
263 (2. * a2 * a2 * std::pow(d_ip, n) -
264 a1 * std::pow(d_ip, n) *
266 n * std::pow(d_ip, n - 2.) * (n - 1.)) -
267 2. * a1 * a2 * n * std::pow(d_ip, n - 1.) +
268 a1 * a1 * n * std::pow(d_ip, n - 2.) * (n - 1.)) /
275template <
int DisplacementDim>
280 switch (phasefield_model)
283 return std::make_unique<COHESIVE_DegradationDerivative>(
284 lch, softening_curve);
286 return std::make_unique<AT_DegradationDerivative>();
290template <
typename T_DNDX,
typename T_N,
typename T_W,
typename T_D,
291 typename T_LOCAL_JAC,
typename T_LOCAL_RHS>
293 T_D& d, T_LOCAL_JAC& local_Jac,
294 T_LOCAL_RHS local_rhs,
295 double const gc,
double const ls,
298 switch (phasefield_model)
302 auto const local_Jac_AT1 =
303 (gc * 0.75 * dNdx.transpose() * dNdx * ls * w).eval();
304 local_Jac.noalias() += local_Jac_AT1;
306 local_rhs.noalias() -=
307 gc * (-0.375 * N.transpose() / ls) * w + local_Jac_AT1 * d;
312 auto const local_Jac_AT2 =
313 (gc * (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
316 local_Jac.noalias() += local_Jac_AT2;
318 local_rhs.noalias() -=
319 local_Jac_AT2 * d - gc * (N.transpose() / ls) * w;
324 auto const local_Jac_COHESIVE =
325 (2.0 / std::numbers::pi * gc *
326 (-N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) * w)
329 local_Jac.noalias() += local_Jac_COHESIVE;
331 local_rhs.noalias() -= local_Jac_COHESIVE * d;
336 OGS_FATAL(
"Invalid phase field model specified.");
341template <
typename T_VECTOR,
typename T_MATRIX,
int DisplacementDim>
343 T_VECTOR& sigma, T_VECTOR& sigma_tensile, T_VECTOR& sigma_compressive,
344 T_VECTOR& eps_tensile, T_MATRIX& D, T_MATRIX& C_tensile,
345 T_MATRIX& C_compressive,
double& strain_energy_tensile,
346 double& elastic_energy,
double const degradation, T_VECTOR
const& eps,
351 auto linear_elastic_mp =
353 DisplacementDim
> const&>(solid_material)
354 .getMaterialProperties();
356 auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
357 auto const mu = linear_elastic_mp.mu(t, x);
358 switch (energy_split_model)
362 std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
363 elastic_energy, C_tensile, C_compressive) =
366 degradation, bulk_modulus, mu, eps);
371 std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
372 elastic_energy, C_tensile, C_compressive) =
374 DisplacementDim>(degradation, bulk_modulus, mu, eps);
379 std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
380 elastic_energy, C_tensile, C_compressive) =
381 MaterialLib::Solids::Phasefield::
382 calculateIsotropicDegradedStressWithRankineEnergy<
383 DisplacementDim>(degradation, bulk_modulus, mu, eps);
388 double temperature = 0.;
391 DisplacementDim
> const&>(solid_material)
392 .getElasticTensor(t, x, temperature);
394 std::tie(eps_tensile, sigma, sigma_tensile, sigma_compressive, D,
395 strain_energy_tensile, elastic_energy, C_tensile,
398 degradation, eps, C_ortho);
403 double temperature = 0.;
406 DisplacementDim
> const&>(solid_material)
407 .getElasticTensor(t, x, temperature);
409 std::tie(eps_tensile, sigma, sigma_tensile, D,
410 strain_energy_tensile, elastic_energy, C_tensile,
413 degradation, eps, C_ortho);
417 OGS_FATAL(
"Invalid energy split model specified.");
double degradation(double const d_ip, double const k, double const) override
double degradationDf1(double const d_ip, double const k, double const) override
double degradationDf2(double const, double const k, double const) override
SofteningCurve softening_curve
double degradationDf1(double const d_ip, double const k, double const ls) override
double degradationDf2(double const d_ip, double const k, double const ls) override
double degradation(double const d_ip, double const k, double const ls) override
COHESIVE_DegradationDerivative(double const Parameter_lch, SofteningCurve Parameter_softening_curve)
virtual double degradationDf2(double const d_ip, double const k, double const ls)=0
virtual double degradation(double const d_ip, double const k, double const ls)=0
virtual ~DegradationDerivative()=default
virtual double degradationDf1(double const d_ip, double const k, double const ls)=0
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > calculateOrthoVolDevDegradedStress(double const degradation, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const &C_ortho)
void calculateStress(T_VECTOR &sigma, T_VECTOR &sigma_tensile, T_VECTOR &sigma_compressive, T_VECTOR &eps_tensile, T_MATRIX &D, T_MATRIX &C_tensile, T_MATRIX &C_compressive, double &strain_energy_tensile, double &elastic_energy, double const degradation, T_VECTOR const &eps, EnergySplitModel const &energy_split_model, double const t, ParameterLib::SpatialPosition const &x, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > calculateIsotropicDegradedStress(double const degradation, double const bulk_modulus, double const mu, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps)
void calculateCrackLocalJacobianAndResidual(T_DNDX &dNdx, T_N &N, T_W &w, T_D &d, T_LOCAL_JAC &local_Jac, T_LOCAL_RHS local_rhs, double const gc, double const ls, PhaseFieldModel &phasefield_model)
std::unique_ptr< DegradationDerivative > creatDegradationDerivative(PhaseFieldModel const &phasefield_model, double const lch, SofteningCurve const &softening_curve)
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > calculateVolDevDegradedStress(double const degradation, double const bulk_modulus, double const mu, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps)
PhaseFieldModel convertStringToPhaseFieldModel(std::string const &phasefield_model)
SofteningCurve convertStringToSofteningCurve(std::optional< std::string > const &softening_curve)
EnergySplitModel convertStringToEnergySplitModel(std::string const &energy_split_model)
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > calculateOrthoMasonryDegradedStress(double const degradation, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const &C_ortho)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType