20 namespace ThermoMechanicalPhaseField
22 template <
typename ShapeFunction,
typename IntegrationMethod,
24 void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
26 assembleWithJacobianForStaggeredScheme(
27 double const t,
double const dt, Eigen::VectorXd
const& local_x,
28 Eigen::VectorXd
const& local_xdot,
const double ,
29 const double ,
int const process_id,
30 std::vector<double>& ,
31 std::vector<double>& ,
32 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
34 if (process_id == _phase_field_process_id)
36 assembleWithJacobianForPhaseFieldEquations(t, local_x, local_b_data,
41 if (process_id == _heat_conduction_process_id)
43 assembleWithJacobianForHeatConductionEquations(
44 t, dt, local_x, local_xdot, local_b_data, local_Jac_data);
49 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
53 template <
typename ShapeFunction,
typename IntegrationMethod,
57 assembleWithJacobianForDeformationEquations(
58 double const t,
double const dt, Eigen::VectorXd
const& local_x,
59 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
61 assert(local_x.size() ==
62 phasefield_size + displacement_size + temperature_size);
64 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
66 local_x.template segment<displacement_size>(displacement_index);
68 local_x.template segment<temperature_size>(temperature_index);
71 typename ShapeMatricesType::template MatrixType<displacement_size,
73 local_Jac_data, displacement_size, displacement_size);
76 typename ShapeMatricesType::template VectorType<displacement_size>>(
77 local_b_data, displacement_size);
82 int const n_integration_points = _integration_method.getNumberOfPoints();
83 for (
int ip = 0; ip < n_integration_points; ip++)
85 x_position.setIntegrationPoint(ip);
86 auto const& w = _ip_data[ip].integration_weight;
87 auto const& N = _ip_data[ip].N;
88 auto const& dNdx = _ip_data[ip].dNdx;
91 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
95 ShapeFunction::NPOINTS,
97 dNdx, N, x_coord, _is_axially_symmetric);
99 auto& eps = _ip_data[ip].eps;
101 auto const& C_tensile = _ip_data[ip].C_tensile;
102 auto const& C_compressive = _ip_data[ip].C_compressive;
104 auto const& sigma = _ip_data[ip].sigma;
106 auto const k = _process_data.residual_stiffness(t, x_position)[0];
107 auto rho_sr = _process_data.solid_density(t, x_position)[0];
108 auto const alpha = _process_data.linear_thermal_expansion_coefficient(
110 double const T0 = _process_data.reference_temperature;
111 auto const& b = _process_data.specific_body_force;
113 double const T_ip = N.dot(T);
114 double const delta_T = T_ip - T0;
116 double const rho_s = rho_sr / (1 + 3 *
alpha * delta_T);
118 double const d_ip = N.dot(d);
119 double const degradation = d_ip * d_ip * (1 - k) + k;
121 auto const C_eff = degradation * C_tensile + C_compressive;
122 eps.noalias() = B * u;
123 _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
alpha,
124 delta_T, degradation);
126 typename ShapeMatricesType::template MatrixType<DisplacementDim,
128 N_u = ShapeMatricesType::template MatrixType<
129 DisplacementDim, displacement_size>::Zero(DisplacementDim,
132 for (
int i = 0; i < DisplacementDim; ++i)
134 N_u.template block<1, displacement_size / DisplacementDim>(
135 i, i * displacement_size / DisplacementDim)
139 local_Jac.noalias() += B.transpose() * C_eff * B * w;
141 local_rhs.noalias() -=
142 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
146 template <
typename ShapeFunction,
typename IntegrationMethod,
150 assembleWithJacobianForHeatConductionEquations(
151 double const t,
double const dt, Eigen::VectorXd
const& local_x,
152 Eigen::VectorXd
const& local_xdot, std::vector<double>& local_b_data,
153 std::vector<double>& local_Jac_data)
155 assert(local_x.size() ==
156 phasefield_size + displacement_size + temperature_size);
158 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
160 local_x.template segment<temperature_size>(temperature_index);
162 local_xdot.template segment<temperature_size>(temperature_index);
164 typename ShapeMatricesType::template MatrixType<temperature_size,
166 local_Jac_data, temperature_size, temperature_size);
169 typename ShapeMatricesType::template VectorType<temperature_size>>(
170 local_b_data, temperature_size);
175 int const n_integration_points = _integration_method.getNumberOfPoints();
176 for (
int ip = 0; ip < n_integration_points; ip++)
178 x_position.setIntegrationPoint(ip);
179 auto const& w = _ip_data[ip].integration_weight;
180 auto const& N = _ip_data[ip].N;
181 auto const& dNdx = _ip_data[ip].dNdx;
183 auto& eps_m = _ip_data[ip].eps_m;
184 auto& heatflux = _ip_data[ip].heatflux;
186 auto rho_sr = _process_data.solid_density(t, x_position)[0];
187 auto const alpha = _process_data.linear_thermal_expansion_coefficient(
189 double const c = _process_data.specific_heat_capacity(t, x_position)[0];
191 _process_data.thermal_conductivity(t, x_position)[0];
192 auto const lambda_res =
193 _process_data.residual_thermal_conductivity(t, x_position)[0];
194 double const T0 = _process_data.reference_temperature;
196 double const d_ip = N.dot(d);
197 double const T_ip = N.dot(T);
198 double const T_dot_ip = N.dot(T_dot);
199 double const delta_T = T_ip - T0;
201 double const rho_s = rho_sr / (1 + 3 *
alpha * delta_T);
203 auto const lambda_eff =
204 d_ip * d_ip * lambda + (1 - d_ip) * (1 - d_ip) * lambda_res;
209 double const eps_m_trace = Invariants::trace(eps_m);
210 if (eps_m_trace >= 0)
212 local_Jac.noalias() += (dNdx.transpose() * lambda_eff * dNdx +
213 N.transpose() * rho_s *
c * N / dt) *
216 local_rhs.noalias() -= (N.transpose() * rho_s *
c * T_dot_ip +
217 dNdx.transpose() * lambda_eff * dNdx * T) *
220 heatflux.noalias() = -(lambda_eff * dNdx * T) * w;
224 local_Jac.noalias() += (dNdx.transpose() * lambda * dNdx +
225 N.transpose() * rho_s *
c * N / dt) *
228 local_rhs.noalias() -= (N.transpose() * rho_s *
c * T_dot_ip +
229 dNdx.transpose() * lambda * dNdx * T) *
232 heatflux.noalias() = -(lambda * dNdx * T) * w;
237 template <
typename ShapeFunction,
typename IntegrationMethod,
241 assembleWithJacobianForPhaseFieldEquations(
242 double const t, Eigen::VectorXd
const& local_x,
243 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
245 assert(local_x.size() ==
246 phasefield_size + displacement_size + temperature_size);
248 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
251 typename ShapeMatricesType::template MatrixType<phasefield_size,
253 local_Jac_data, phasefield_size, phasefield_size);
256 typename ShapeMatricesType::template VectorType<phasefield_size>>(
257 local_b_data, phasefield_size);
262 int const n_integration_points = _integration_method.getNumberOfPoints();
263 for (
int ip = 0; ip < n_integration_points; ip++)
265 x_position.setIntegrationPoint(ip);
266 auto const& w = _ip_data[ip].integration_weight;
267 auto const& N = _ip_data[ip].N;
268 auto const& dNdx = _ip_data[ip].dNdx;
270 auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
272 auto const gc = _process_data.crack_resistance(t, x_position)[0];
273 auto const ls = _process_data.crack_length_scale(t, x_position)[0];
275 double const d_ip = N.dot(d);
277 local_Jac.noalias() += (dNdx.transpose() * gc * ls * dNdx +
278 N.transpose() * 2 * strain_energy_tensile * N +
279 N.transpose() * gc / ls * N) *
282 local_rhs.noalias() -=
283 (dNdx.transpose() * gc * ls * dNdx * d +
284 N.transpose() * d_ip * 2 * strain_energy_tensile -
285 N.transpose() * gc / ls * (1 - d_ip)) *
void setElementID(std::size_t element_id)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.