OGS
MaterialLib::Solids::Phasefield Namespace Reference

Functions

double heaviside (double const v)
 
double macaulayTensile (double const v)
 
double macaulayCompressive (double v)
 
template<int DisplacementDim>
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)
 
template<int DisplacementDim>
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)
 
template<int DisplacementDim>
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > calculateIsotropicDegradedStressWithRankineEnergy (double const degradation, double const bulk_modulus, double const mu, MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &eps)
 
template<int DisplacementDim>
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)
 
template<int DisplacementDim>
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)
 

Function Documentation

◆ calculateIsotropicDegradedStress()

template<int DisplacementDim>
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > MaterialLib::Solids::Phasefield::calculateIsotropicDegradedStress ( double const degradation,
double const bulk_modulus,
double const mu,
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const & eps )

Definition at line 138 of file LinearElasticIsotropicPhaseField.h.

143{
144 static int const KelvinVectorSize =
146 using KelvinVector =
148 using KelvinMatrix =
151 // calculation of deviatoric parts
152 auto const& P_dev = Invariants::deviatoric_projection;
153 KelvinVector const epsd_curr = P_dev * eps;
154
155 // Hydrostatic part for the stress and the tangent.
156 double const eps_curr_trace = Invariants::trace(eps);
157
158 KelvinMatrix C_tensile = KelvinMatrix::Zero();
159 KelvinMatrix C_compressive = KelvinMatrix::Zero();
160
161 double const strain_energy_tensile =
162 bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
163 mu * epsd_curr.transpose() * epsd_curr;
164 double const elastic_energy = degradation * strain_energy_tensile;
165 KelvinVector const sigma_tensile =
166 bulk_modulus * eps_curr_trace * Invariants::identity2 +
167 2 * mu * epsd_curr;
168 C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
169 C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
170 KelvinVector const sigma_real = degradation * sigma_tensile;
171
172 KelvinMatrix const D = degradation * C_tensile + C_compressive;
173
174 return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
175 elastic_energy, C_tensile, C_compressive);
176}
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ calculateIsotropicDegradedStressWithRankineEnergy()

template<int DisplacementDim>
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > MaterialLib::Solids::Phasefield::calculateIsotropicDegradedStressWithRankineEnergy ( double const degradation,
double const bulk_modulus,
double const mu,
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const & eps )

Definition at line 191 of file LinearElasticIsotropicPhaseField.h.

196{
197 static int const KelvinVectorSize =
199 using KelvinVector =
201 using KelvinMatrix =
204 // calculation of deviatoric parts
205 auto const& P_dev = Invariants::deviatoric_projection;
206 KelvinVector const epsd_curr = P_dev * eps;
207
208 // Hydrostatic part for the stress and the tangent.
209 double const eps_curr_trace = Invariants::trace(eps);
210
211 KelvinMatrix C_tensile = KelvinMatrix::Zero();
212 KelvinMatrix C_compressive = KelvinMatrix::Zero();
213
214 KelvinVector const sigma_tensile =
215 bulk_modulus * eps_curr_trace * Invariants::identity2 +
216 2 * mu * epsd_curr;
217 // The maximum principal stress
218 auto sigma_tensor =
220 Eigen::SelfAdjointEigenSolver<decltype(sigma_tensor)>
221 selfAdjoint_eigen_solver(sigma_tensor);
222 double const max_principle_stress =
223 selfAdjoint_eigen_solver
224 .eigenvalues()[2]; // selfAdjoint_eigen_solver function gives an
225 // ascending sequence.
226 //
227 double const E0 = 9.0 * mu * bulk_modulus / (3.0 * bulk_modulus + mu);
228 auto hs = [&](double const v) { return heaviside(v); };
229 //
230 double const strain_energy_tensile = 0.5 * max_principle_stress *
231 max_principle_stress *
232 hs(max_principle_stress) / E0;
233 double const elastic_energy =
234 degradation * (bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
235 mu * epsd_curr.transpose() * epsd_curr);
236 //
237 C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
238 C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
239 KelvinVector const sigma_real = degradation * sigma_tensile;
240
241 KelvinMatrix const D = degradation * C_tensile + C_compressive;
242 return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
243 elastic_energy, C_tensile, C_compressive);
244}
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)

References heaviside(), MathLib::KelvinVector::kelvin_vector_dimensions(), and MathLib::KelvinVector::kelvinVectorToTensor().

◆ calculateOrthoMasonryDegradedStress()

template<int DisplacementDim>
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 > > MaterialLib::Solids::Phasefield::calculateOrthoMasonryDegradedStress ( double const degradation,
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const & eps,
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const & C_ortho )

Definition at line 127 of file LinearElasticOrthotropicPhaseField.h.

131{
132 using KelvinVector =
134 using KelvinMatrix =
136
137 KelvinMatrix const C = C_ortho;
138
139 // calculation of square root of elasticity tensor
140 Eigen::SelfAdjointEigenSolver<KelvinMatrix> es(C);
141 KelvinMatrix const sqrt_C = es.operatorSqrt();
142 Eigen::SelfAdjointEigenSolver<KelvinMatrix> es_inverse(C.inverse());
143 KelvinMatrix const sqrt_C_inverse = es_inverse.operatorSqrt();
144
145 // strain in a transformed space
146 KelvinVector const epst = sqrt_C * eps;
147
148 // strain tensor in 3D
149 Eigen::Matrix3d const eps_3D =
151
152 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> const es_eps_3D(eps_3D);
153 Eigen::Vector3d const eigen_values_eps_3D = es_eps_3D.eigenvalues().real();
154 Eigen::Matrix3d const eigen_vectors_eps_3D = es_eps_3D.eigenvectors();
155 Eigen::Matrix3d const E1_eigenp =
156 eigen_vectors_eps_3D.col(0) * eigen_vectors_eps_3D.col(0).transpose();
157
158 Eigen::Matrix3d const E3_eigenp =
159 eigen_vectors_eps_3D.col(2) * eigen_vectors_eps_3D.col(2).transpose();
160
161 KelvinVector const E1_vec =
162 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(E1_eigenp);
163
164 KelvinVector const E3_vec =
165 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(E3_eigenp);
166
167 KelvinMatrix const E1oE1 = E1_vec * E1_vec.transpose();
168 KelvinMatrix const E3oE3 = E3_vec * E3_vec.transpose();
169
170 KelvinMatrix I_p = KelvinMatrix::Zero();
171 KelvinMatrix I_n = KelvinMatrix::Zero();
172
173 KelvinMatrix const I_S = KelvinMatrix::Identity();
174 if (DisplacementDim == 2)
175 {
176 if (std::abs(eigen_values_eps_3D(2) - eigen_values_eps_3D(0)) >
177 std::numeric_limits<double>::epsilon())
178 {
179 I_p = (macaulayTensile(eigen_values_eps_3D(0)) -
180 macaulayTensile(eigen_values_eps_3D(2))) /
181 (eigen_values_eps_3D(0) - eigen_values_eps_3D(2)) *
182 (I_S - (E1oE1 + E3oE3)) +
183 (heaviside(eigen_values_eps_3D(0)) * E1oE1 +
184 heaviside(eigen_values_eps_3D(2)) * E3oE3);
185 I_n = (macaulayCompressive(eigen_values_eps_3D(0)) -
186 macaulayCompressive(eigen_values_eps_3D(2))) /
187 (eigen_values_eps_3D(0) - eigen_values_eps_3D(2)) *
188 (I_S - (E1oE1 + E3oE3)) +
189 (heaviside(-eigen_values_eps_3D(0)) * E1oE1 +
190 heaviside(-eigen_values_eps_3D(2)) * E3oE3);
191 }
192 else
193 {
194 I_p = heaviside(eigen_values_eps_3D(0)) * I_S;
195 I_n = heaviside(-eigen_values_eps_3D(0)) * I_S;
196 }
197 }
198
199 // strain tensile / compressive
200 KelvinVector const eps_tensile = sqrt_C_inverse * (I_p * sqrt_C) * eps;
201 KelvinVector const eps_compressive = sqrt_C_inverse * (I_n * sqrt_C) * eps;
202
203 // C tensile / compressive
204 KelvinMatrix const C_tensile = sqrt_C * I_p * sqrt_C;
205 KelvinMatrix const C_compressive = sqrt_C * I_n * sqrt_C;
206
207 // stress tensile / compressive
208 KelvinVector const sigma_tensile = C * eps_tensile;
209 KelvinVector const sigma_compressive = C * eps_compressive;
210
211 // decomposition of strain energy density
212 double const strain_energy_tensile =
213 0.5 * sigma_tensile.adjoint() * eps_tensile;
214 double const strain_energy_compressive =
215 0.5 * sigma_compressive.adjoint() * eps_compressive;
216
217 double const elastic_energy =
218 degradation * strain_energy_tensile + strain_energy_compressive;
219
220 KelvinVector const sigma_real =
221 degradation * sigma_tensile + sigma_compressive;
222
223 KelvinMatrix const D = degradation * C_tensile + C_compressive;
224
225 return std::make_tuple(eps_tensile, sigma_real, sigma_tensile, D,
226 strain_energy_tensile, elastic_energy, C_tensile,
227 C_compressive);
228}

References heaviside(), MathLib::KelvinVector::kelvinVectorToTensor(), macaulayCompressive(), and macaulayTensile().

◆ calculateOrthoVolDevDegradedStress()

template<int DisplacementDim>
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 > > MaterialLib::Solids::Phasefield::calculateOrthoVolDevDegradedStress ( double const degradation,
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const & eps,
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const & C_ortho )

Definition at line 38 of file LinearElasticOrthotropicPhaseField.h.

42{
43 static constexpr int KelvinVectorSize =
45 using KelvinVector =
47 using KelvinMatrix =
50 // calculation of deviatoric parts
51 auto const& P_dev = Invariants::deviatoric_projection;
52
53 KelvinMatrix const C = C_ortho;
54 // calculation of square root of elasticity tensor
55 Eigen::SelfAdjointEigenSolver<KelvinMatrix> es(C);
56 KelvinMatrix const sqrt_C = es.operatorSqrt();
57 Eigen::SelfAdjointEigenSolver<KelvinMatrix> es_inverse(C.inverse());
58 KelvinMatrix const sqrt_C_inverse = es_inverse.operatorSqrt();
59
60 // strain in a transformed space
61 KelvinVector const epst = sqrt_C * eps;
62 double const epst_curr_trace = Invariants::trace(epst);
63
64 // projection tensors in transformed space
65 KelvinMatrix teps_p = KelvinMatrix::Zero();
66 KelvinMatrix teps_n = KelvinMatrix::Zero();
67 if (epst_curr_trace >= 0) /* QQQ */
68 {
69 teps_p.template topLeftCorner<3, 3>().setConstant(1. / 3);
70 }
71 else
72 {
73 teps_n.template topLeftCorner<3, 3>().setConstant(1. / 3);
74 }
75
76 teps_p.noalias() += P_dev * KelvinMatrix::Identity();
77
78 // strain tensile / compressive
79 KelvinVector const eps_tensile = sqrt_C_inverse * (teps_p * sqrt_C) * eps;
80 KelvinVector const eps_compressive =
81 sqrt_C_inverse * (teps_n * sqrt_C) * eps;
82
83 // projection tensors in original space
84 KelvinMatrix const der_eps_p = sqrt_C_inverse * (teps_p * sqrt_C);
85 KelvinMatrix const der_eps_n = sqrt_C_inverse * (teps_n * sqrt_C);
86 // C tensile / compressive
87 KelvinMatrix const C_tensile = der_eps_p.transpose() * C * der_eps_p;
88 KelvinMatrix const C_compressive = der_eps_n.transpose() * C * der_eps_n;
89
90 // stress tensile / compressive
91 KelvinVector const sigma_tensile = C_tensile * eps;
92 KelvinVector const sigma_compressive = C_compressive * eps;
93
94 // decomposition of strain energy density
95 double const strain_energy_tensile =
96 0.5 * sigma_tensile.adjoint() * eps_tensile;
97 double const strain_energy_compressive =
98 0.5 * sigma_compressive.adjoint() * eps_compressive;
99
100 double const elastic_energy =
101 degradation * strain_energy_tensile + strain_energy_compressive;
102
103 KelvinVector const sigma_real =
104 degradation * sigma_tensile + sigma_compressive;
105 KelvinMatrix const D = degradation * C_tensile + C_compressive;
106
107 return std::make_tuple(eps_tensile, sigma_real, sigma_tensile,
108 sigma_compressive, D, strain_energy_tensile,
109 elastic_energy, C_tensile, C_compressive);
110}

References MathLib::KelvinVector::kelvin_vector_dimensions().

◆ calculateVolDevDegradedStress()

template<int DisplacementDim>
std::tuple< MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinVectorType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, double, double, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim >, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > > MaterialLib::Solids::Phasefield::calculateVolDevDegradedStress ( double const degradation,
double const bulk_modulus,
double const mu,
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const & eps )

Definition at line 61 of file LinearElasticIsotropicPhaseField.h.

64{
65 static int const KelvinVectorSize =
67 using KelvinVector =
69 using KelvinMatrix =
72 // calculation of deviatoric parts
73 auto const& P_dev = Invariants::deviatoric_projection;
74 KelvinVector const epsd_curr = P_dev * eps;
75
76 // Hydrostatic part for the stress and the tangent.
77 double const eps_curr_trace = Invariants::trace(eps);
78
79 KelvinMatrix C_tensile = KelvinMatrix::Zero();
80 KelvinMatrix C_compressive = KelvinMatrix::Zero();
81
82 auto strain_energy_computation_vol = [&](auto&& macaulay)
83 {
84 auto macaulay_squared = [&macaulay](double x)
85 { return boost::math::pow<2>(macaulay(x)); };
86 return bulk_modulus / 2 * macaulay_squared(eps_curr_trace);
87 };
88
89 auto stress_computation_vol = [&](auto&& macaulay)
90 { return bulk_modulus * macaulay(eps_curr_trace) * Invariants::identity2; };
91
92 auto hs = [&](double const v) { return heaviside(v); };
93
94 auto mt = [&](double const v) { return macaulayTensile(v); };
95
96 auto mc = [&](double const v) { return macaulayCompressive(v); };
97
98 double const strain_energy_tensile = strain_energy_computation_vol(mt) +
99 mu * epsd_curr.transpose() * epsd_curr;
100
101 KelvinVector const sigma_tensile =
102 stress_computation_vol(mt) + 2 * mu * epsd_curr;
103
104 KelvinVector const sigma_compressive = stress_computation_vol(mc);
105
106 C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus *
107 hs(eps_curr_trace));
108 C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
109
110 C_compressive.template topLeftCorner<3, 3>().setConstant(
111 bulk_modulus * (1 - hs(eps_curr_trace)));
112
113 double const elastic_energy =
114 bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
115 mu * epsd_curr.transpose() * epsd_curr;
116 KelvinVector const sigma_real =
117 degradation * sigma_tensile + sigma_compressive;
118
119 KelvinMatrix const D = degradation * C_tensile + C_compressive;
120
121 return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
122 elastic_energy, C_tensile, C_compressive);
123}
KV::KelvinMatrixType< DisplacementDim > KelvinMatrix
Definition Base.h:30
KV::KelvinVectorType< DisplacementDim > KelvinVector
Definition Base.h:27

References heaviside(), MathLib::KelvinVector::kelvin_vector_dimensions(), macaulayCompressive(), and macaulayTensile().

Referenced by ProcessLib::ThermoMechanicalPhaseField::IntegrationPointData< BMatricesType, ShapeMatrixType, DisplacementDim >::updateConstitutiveRelation().

◆ heaviside()

double MaterialLib::Solids::Phasefield::heaviside ( double const v)
inline

Decompose the stiffness into tensile and compressive part. Judging by the physical observations, compression perpendicular to a crack does not cause crack propagation. Thus, the phase-field parameter is only involved into the tensile part to degrade the elastic strain energy. heaviside function returns 1.0 if the argument is positive and 0.0 if negative

Definition at line 32 of file LinearElasticIsotropicPhaseField.h.

33{
34 return (v < 0) ? 0.0 : 1.0;
35}

Referenced by calculateIsotropicDegradedStressWithRankineEnergy(), calculateOrthoMasonryDegradedStress(), calculateVolDevDegradedStress(), macaulayCompressive(), and macaulayTensile().

◆ macaulayCompressive()

double MaterialLib::Solids::Phasefield::macaulayCompressive ( double v)
inline

Definition at line 43 of file LinearElasticIsotropicPhaseField.h.

44{
45 return v * (1 - heaviside(v));
46}

References heaviside().

Referenced by calculateOrthoMasonryDegradedStress(), and calculateVolDevDegradedStress().

◆ macaulayTensile()

double MaterialLib::Solids::Phasefield::macaulayTensile ( double const v)
inline

Macaulay brackets: positive strain is tensile and negative strain for compressive

Definition at line 39 of file LinearElasticIsotropicPhaseField.h.

40{
41 return v * heaviside(v);
42}

References heaviside().

Referenced by calculateOrthoMasonryDegradedStress(), and calculateVolDevDegradedStress().