Loading [MathJax]/extensions/MathMenu.js
OGS
PhaseFieldBase.h
Go to the documentation of this file.
1
10#pragma once
11
12#include <Eigen/Core>
13#include <numbers>
14#include <optional>
15
16#include "BaseLib/Logging.h"
21
22namespace MaterialLib
23{
24namespace Solids
25{
26namespace Phasefield
27{
29{
30 AT1,
31 AT2,
33};
34
35template <int DisplacementDim>
37 std::string const& phasefield_model)
38{
39 if (phasefield_model == "AT1")
40 {
42 }
43 if (phasefield_model == "AT2")
44 {
46 }
47 if (phasefield_model == "COHESIVE")
48 {
50 }
52 "phasefield_model must be 'AT1', 'AT2' or 'COHESIVE' but '{}' "
53 "was given",
54 phasefield_model.c_str());
55};
56
58{
59 Linear,
61};
62
63template <int DisplacementDim>
65 std::optional<std::string> const& softening_curve)
66{
67 if (softening_curve)
68 {
69 if (*softening_curve == "Linear")
70 {
72 }
73 if (*softening_curve == "Exponential")
74 {
76 }
78 "softening_curve must be 'Linear' or 'Exponential' but '{}' "
79 "was given",
80 softening_curve->c_str());
81 }
82 return SofteningCurve::Linear; // default
83};
84
93
94template <int DisplacementDim>
96 std::string const& energy_split_model)
97{
98 if (energy_split_model == "Isotropic")
99 {
101 }
102 if (energy_split_model == "VolumetricDeviatoric")
103 {
105 }
106 if (energy_split_model == "EffectiveStress")
107 {
109 }
110 if (energy_split_model == "OrthoVolDev")
111 {
113 }
114 if (energy_split_model == "OrthoMasonry")
115 {
117 }
118 OGS_FATAL(
119 "energy_split_model must be 'Isotropic', 'VolumetricDeviatoric', "
120 "'EffectiveStress', 'OrthoVolDev' or 'OrthoMasonry' but '{}' was given",
121 energy_split_model.c_str());
122};
123
125{
126public:
127 virtual ~DegradationDerivative() = default;
128 virtual double degradation(double const d_ip,
129 double const k,
130 double const ls) = 0; /* degradation_df0 */
131 virtual double degradationDf1(double const d_ip, double const k,
132 double const ls) = 0;
133 virtual double degradationDf2(double const d_ip, double const k,
134 double const ls) = 0;
135};
136
138{
139public:
140 double degradation(double const d_ip,
141 double const k,
142 double const /* ls */) override
143 {
144 return d_ip * d_ip * (1. - k) + k;
145 };
146 double degradationDf1(double const d_ip, double const k,
147 double const /* ls */) override
148 {
149 return 2. * (1. - k) * d_ip;
150 };
151 double degradationDf2(double const /* d_ip */, double const k,
152 double const /* ls */) override
153 {
154 return 2. * (1. - k);
155 };
156};
157
159{
160private:
161 double const lch;
163
164public:
165 COHESIVE_DegradationDerivative(double const Parameter_lch,
166 SofteningCurve Parameter_softening_curve)
167 : lch(Parameter_lch), softening_curve(Parameter_softening_curve){};
168 double degradation(double const d_ip,
169 double const k,
170 double const ls) override
171 {
172 double const m1 = 4.0 * lch / acos(-1.) / ls;
173 switch (softening_curve)
174 {
176 {
177 double const m2 = std::pow(2., 5. / 3.) - 3.;
178 double const n = 2.5;
179 return std::pow(d_ip, n) /
180 (std::pow(d_ip, n) +
181 m1 * (1. - d_ip) * (1. + m2 * (1. - d_ip))) *
182 (1.0 - k) +
183 k;
184 }
185 default:
186 {
187 double const m2 = -0.5;
188 double const n = 2.;
189 return std::pow(d_ip, n) /
190 (std::pow(d_ip, n) +
191 m1 * (1. - d_ip) * (1. + m2 * (1. - d_ip))) *
192 (1. - k) +
193 k;
194 }
195 }
196 };
197 double degradationDf1(double const d_ip, double const k,
198 double const ls) override
199 {
200 double const m1 = 4.0 * lch / acos(-1.) / ls;
201 switch (softening_curve)
202 {
204 {
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;
211 return (1. - k) *
212 (n * std::pow(d_ip, n - 1.) * a1 -
213 std::pow(d_ip, n) * a2) /
214 (a1 * a1);
215 }
216 default:
217 {
218 double const m2 = -0.5;
219 double const n = 2.;
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;
224 return (1. - k) *
225 (n * std::pow(d_ip, n - 1.) * a1 -
226 std::pow(d_ip, n) * a2) /
227 (a1 * a1);
228 }
229 }
230 };
231 double degradationDf2(double const d_ip, double const k,
232 double const ls) override
233 {
234 double const m1 = 4.0 * lch / acos(-1.) / ls;
235 switch (softening_curve)
236 {
238 {
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;
245 return (1. - k) *
246 (2. * a2 * a2 * std::pow(d_ip, n) -
247 a1 * std::pow(d_ip, n) *
248 (2. * m1 * m2 +
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.)) /
252 (std::pow(a1, 3));
253 }
254 default:
255 {
256 double const m2 = -0.5;
257 double const n = 2.;
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;
262 return (1. - k) *
263 (2. * a2 * a2 * std::pow(d_ip, n) -
264 a1 * std::pow(d_ip, n) *
265 (2. * m1 * m2 +
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.)) /
269 (std::pow(a1, 3));
270 }
271 }
272 };
273};
274
275template <int DisplacementDim>
276std::unique_ptr<DegradationDerivative> creatDegradationDerivative(
277 PhaseFieldModel const& phasefield_model, double const lch,
278 SofteningCurve const& softening_curve)
279{
280 switch (phasefield_model)
281 {
283 return std::make_unique<COHESIVE_DegradationDerivative>(
284 lch, softening_curve);
285 default:
286 return std::make_unique<AT_DegradationDerivative>();
287 }
288};
289
290template <typename T_DNDX, typename T_N, typename T_W, typename T_D,
291 typename T_LOCAL_JAC, typename T_LOCAL_RHS>
292void calculateCrackLocalJacobianAndResidual(T_DNDX& dNdx, T_N& N, T_W& w,
293 T_D& d, T_LOCAL_JAC& local_Jac,
294 T_LOCAL_RHS local_rhs,
295 double const gc, double const ls,
296 PhaseFieldModel& phasefield_model)
297{
298 switch (phasefield_model)
299 {
301 {
302 auto const local_Jac_AT1 =
303 (gc * 0.75 * dNdx.transpose() * dNdx * ls * w).eval();
304 local_Jac.noalias() += local_Jac_AT1;
305
306 local_rhs.noalias() -=
307 gc * (-0.375 * N.transpose() / ls) * w + local_Jac_AT1 * d;
308 break;
309 }
311 {
312 auto const local_Jac_AT2 =
313 (gc * (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
314 w)
315 .eval();
316 local_Jac.noalias() += local_Jac_AT2;
317
318 local_rhs.noalias() -=
319 local_Jac_AT2 * d - gc * (N.transpose() / ls) * w;
320 break;
321 }
323 {
324 auto const local_Jac_COHESIVE =
325 (2.0 / std::numbers::pi * gc *
326 (-N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) * w)
327 .eval();
328
329 local_Jac.noalias() += local_Jac_COHESIVE;
330
331 local_rhs.noalias() -= local_Jac_COHESIVE * d;
332 break;
333 }
334 default:
335 {
336 OGS_FATAL("Invalid phase field model specified.");
337 }
338 }
339};
340
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,
347 EnergySplitModel const& energy_split_model, double const t,
350{
351 auto linear_elastic_mp =
353 DisplacementDim> const&>(solid_material)
354 .getMaterialProperties();
355
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)
359 {
361 {
362 std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
363 elastic_energy, C_tensile, C_compressive) =
366 degradation, bulk_modulus, mu, eps);
367 break;
368 }
370 {
371 std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
372 elastic_energy, C_tensile, C_compressive) =
374 DisplacementDim>(degradation, bulk_modulus, mu, eps);
375 break;
376 }
378 {
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);
384 break;
385 }
387 {
388 double temperature = 0.;
391 DisplacementDim> const&>(solid_material)
392 .getElasticTensor(t, x, temperature);
393
394 std::tie(eps_tensile, sigma, sigma_tensile, sigma_compressive, D,
395 strain_energy_tensile, elastic_energy, C_tensile,
396 C_compressive) = MaterialLib::Solids::Phasefield::
398 degradation, eps, C_ortho);
399 break;
400 }
402 {
403 double temperature = 0.;
406 DisplacementDim> const&>(solid_material)
407 .getElasticTensor(t, x, temperature);
408
409 std::tie(eps_tensile, sigma, sigma_tensile, D,
410 strain_energy_tensile, elastic_energy, C_tensile,
411 C_compressive) = MaterialLib::Solids::Phasefield::
413 degradation, eps, C_ortho);
414 break;
415 }
416 default:
417 OGS_FATAL("Invalid energy split model specified.");
418 }
419};
420} // namespace Phasefield
421} // namespace Solids
422} // namespace MaterialLib
#define OGS_FATAL(...)
Definition Error.h:26
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
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 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