39 double const degradation,
43 static constexpr int KelvinVectorSize =
51 auto const& P_dev = Invariants::deviatoric_projection;
53 KelvinMatrix
const C = C_ortho;
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();
61 KelvinVector
const epst = sqrt_C * eps;
62 double const epst_curr_trace = Invariants::trace(epst);
65 KelvinMatrix teps_p = KelvinMatrix::Zero();
66 KelvinMatrix teps_n = KelvinMatrix::Zero();
67 if (epst_curr_trace >= 0)
69 teps_p.template topLeftCorner<3, 3>().setConstant(1. / 3);
73 teps_n.template topLeftCorner<3, 3>().setConstant(1. / 3);
76 teps_p.noalias() += P_dev * KelvinMatrix::Identity();
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;
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);
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;
91 KelvinVector
const sigma_tensile = C_tensile * eps;
92 KelvinVector
const sigma_compressive = C_compressive * eps;
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;
100 double const elastic_energy =
101 degradation * strain_energy_tensile + strain_energy_compressive;
103 KelvinVector
const sigma_real =
104 degradation * sigma_tensile + sigma_compressive;
105 KelvinMatrix
const D = degradation * C_tensile + C_compressive;
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);
128 double const degradation,
137 KelvinMatrix
const C = C_ortho;
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();
146 KelvinVector
const epst = sqrt_C * eps;
149 Eigen::Matrix3d
const eps_3D =
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();
158 Eigen::Matrix3d
const E3_eigenp =
159 eigen_vectors_eps_3D.col(2) * eigen_vectors_eps_3D.col(2).transpose();
161 KelvinVector
const E1_vec =
162 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(E1_eigenp);
164 KelvinVector
const E3_vec =
165 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(E3_eigenp);
167 KelvinMatrix
const E1oE1 = E1_vec * E1_vec.transpose();
168 KelvinMatrix
const E3oE3 = E3_vec * E3_vec.transpose();
170 KelvinMatrix I_p = KelvinMatrix::Zero();
171 KelvinMatrix I_n = KelvinMatrix::Zero();
173 KelvinMatrix
const I_S = KelvinMatrix::Identity();
174 if (DisplacementDim == 2)
176 if (std::abs(eigen_values_eps_3D(2) - eigen_values_eps_3D(0)) >
177 std::numeric_limits<double>::epsilon())
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);
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);
194 I_p =
heaviside(eigen_values_eps_3D(0)) * I_S;
195 I_n =
heaviside(-eigen_values_eps_3D(0)) * I_S;
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;
204 KelvinMatrix
const C_tensile = sqrt_C * I_p * sqrt_C;
205 KelvinMatrix
const C_compressive = sqrt_C * I_n * sqrt_C;
208 KelvinVector
const sigma_tensile = C * eps_tensile;
209 KelvinVector
const sigma_compressive = C * eps_compressive;
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;
217 double const elastic_energy =
218 degradation * strain_energy_tensile + strain_energy_compressive;
220 KelvinVector
const sigma_real =
221 degradation * sigma_tensile + sigma_compressive;
223 KelvinMatrix
const D = degradation * C_tensile + C_compressive;
225 return std::make_tuple(eps_tensile, sigma_real, sigma_tensile, D,
226 strain_energy_tensile, elastic_energy, C_tensile,
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)
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)