OGS
LinearElasticTransverseIsotropic.cpp
Go to the documentation of this file.
1
13
14#include <Eigen/LU>
15
17
18namespace MPL = MaterialPropertyLib;
19
20namespace MaterialLib
21{
22namespace Solids
23{
24template <int DisplacementDim>
25std::optional<std::tuple<typename MechanicsBase<DisplacementDim>::KelvinVector,
26 std::unique_ptr<typename MechanicsBase<
27 DisplacementDim>::MaterialStateVariables>,
30 MaterialPropertyLib::VariableArray const& variable_array_prev,
31 MaterialPropertyLib::VariableArray const& variable_array, double const t,
32 ParameterLib::SpatialPosition const& x, double const /*dt*/,
34 MaterialStateVariables const& /*material_state_variables*/) const
35{
36 auto const& eps_m = std::get<MPL::SymmetricTensor<DisplacementDim>>(
37 variable_array.mechanical_strain);
38 auto const& eps_m_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
39 variable_array_prev.mechanical_strain);
40 auto const& sigma_prev = std::get<MPL::SymmetricTensor<DisplacementDim>>(
41 variable_array_prev.stress);
42 auto const& T = variable_array_prev.temperature;
43
44 KelvinMatrix C = getElasticTensor(t, x, T);
45
46 KelvinVector sigma = sigma_prev + C * (eps_m - eps_m_prev);
47
48 return {std::make_tuple(
49 sigma,
50 std::make_unique<
52 C)};
53}
54
55template <int DisplacementDim>
59 ParameterLib::SpatialPosition const& x) const
60{
61 using namespace MathLib::KelvinVector;
62 double const E_i = E_i_p_(t, x)[0];
63 double const E_a = E_a_p_(t, x)[0];
64 double const nu_i = nu_ii_p_(t, x)[0];
65 double const nu_ia = nu_ia_p_(t, x)[0];
66
67 double const nu_ai = nu_ia * (E_a / E_i);
68
69 double const nu_ia_p_nu_ai = nu_ia * nu_ai;
70 double const nu_i_p2 = nu_i * nu_i;
71 double const E_i_p2 = E_i * E_i;
72 // 1/D
73 double const one_over_D =
74 (E_i_p2 * E_a) / (1.0 - nu_i_p2 - 2.0 * (1.0 + nu_i) * nu_ia_p_nu_ai);
75
76 double const fac1 = one_over_D / (E_i * E_a);
77 double const fac2 = one_over_D / E_i_p2;
78 double const a_ii = (1.0 - nu_ia_p_nu_ai) * fac1;
79 double const a_ai = (1.0 - nu_i_p2) * fac2;
80 double const b_ii = (nu_i + nu_ia_p_nu_ai) * fac1;
81 double const b_ai = (nu_ia * (1.0 + nu_i)) * fac2;
82
85
86 if (DisplacementDim == 2)
87 {
88 // The following data are set based on the condition that the direction
89 // of anisotropy is parallel to the second base (base2) of the 2D local
90 // coordinate system.
91 C_ortho.template topLeftCorner<3, 3>() << a_ii, b_ai, b_ii, b_ai, a_ai,
92 b_ai, b_ii, b_ai, a_ii;
93
94 return C_ortho;
95 }
96 C_ortho.template topLeftCorner<3, 3>() << a_ii, b_ii, b_ai, b_ii, a_ii,
97 b_ai, b_ai, b_ai, a_ai;
98
99 return C_ortho;
100}
101
102template <>
105 double const t, ParameterLib::SpatialPosition const& x,
106 double const /*T*/) const
107{
108 using namespace MathLib::KelvinVector;
109
110 double const G_a = G_ia_p_(t, x)[0];
111 double const c_ai = G_a;
112
113 KelvinMatrixType<2> C_ortho = getElasticTensorLeftTopCorner(t, x);
114
115 C_ortho.template bottomRightCorner<1, 1>().diagonal() << 2 * c_ai;
116
117 auto const Q = [this, &x]() -> KelvinMatrixType<2>
118 {
119 if (!local_coordinate_system_)
120 {
122 }
123
124 Eigen::Matrix<double, 2, 2> R = Eigen::Matrix<double, 2, 2>::Identity();
125 R.template topLeftCorner<2, 2>().noalias() =
126 local_coordinate_system_->transformation<2>(x);
127 return fourthOrderRotationMatrix(R);
128 }();
129
130 // Rotate the forth-order tenser in Kelvin mapping with Q*C_ortho*Q^T and
131 // return the top left corner block of size 4x4 for two-dimensional case.
132 return Q * C_ortho * Q.transpose();
133}
134
135template <>
138 double const t, ParameterLib::SpatialPosition const& x,
139 double const /*T*/) const
140{
141 using namespace MathLib::KelvinVector;
142
143 double const E_i = E_i_p_(t, x)[0];
144 double const nu_i = nu_ii_p_(t, x)[0];
145 double const G_a = G_ia_p_(t, x)[0];
146
147 double const c_ii = E_i / (2.0 * (1 + nu_i));
148 double const c_ai = G_a;
149
150 KelvinMatrixType<3> C_ortho = getElasticTensorLeftTopCorner(t, x);
151
152 C_ortho.template bottomRightCorner<3, 3>().diagonal() << 2.0 * c_ii,
153 2 * c_ai, 2 * c_ai;
154
155 auto const Q = [this, &x]() -> KelvinMatrixType<3>
156 {
157 if (!local_coordinate_system_)
158 {
160 }
161 Eigen::Matrix3d R = Eigen::Matrix3d::Identity();
162 R.template topLeftCorner<3, 3>().noalias() =
163 local_coordinate_system_->transformation<3>(x);
164 return fourthOrderRotationMatrix(R);
165 }();
166
167 // Rotate the forth-order tenser in Kelvin mapping with Q*C_ortho*Q^T and
168 // return the 6x6 matrix for in the three-dimensional case.
169 return Q * C_ortho * Q.transpose();
170}
171
174
175} // namespace Solids
176} // namespace MaterialLib
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix
KelvinMatrix getElasticTensorLeftTopCorner(double const t, ParameterLib::SpatialPosition const &x) const
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > KelvinVector
KelvinMatrix getElasticTensor(double const t, ParameterLib::SpatialPosition const &x, double const T) const
std::optional< std::tuple< typename MechanicsBase< DisplacementDim >::KelvinVector, std::unique_ptr< typename MechanicsBase< DisplacementDim >::MaterialStateVariables >, typename MechanicsBase< DisplacementDim >::KelvinMatrix > > integrateStress(MaterialPropertyLib::VariableArray const &variable_array_prev, MaterialPropertyLib::VariableArray const &variable_array, double const t, ParameterLib::SpatialPosition const &x, double const dt, typename MechanicsBase< DisplacementDim >::MaterialStateVariables const &material_state_variables) const override
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > KelvinMatrix