29 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
30 typename IntegrationMethod,
int DisplacementDim>
31 TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
32 IntegrationMethod, DisplacementDim>::
35 bool const is_axially_symmetric,
36 unsigned const integration_order,
38 : _process_data(process_data),
39 _integration_method(integration_order),
41 _is_axially_symmetric(is_axially_symmetric)
43 unsigned const n_integration_points =
46 _ip_data.reserve(n_integration_points);
49 auto const shape_matrices_u =
52 DisplacementDim>(e, is_axially_symmetric,
55 auto const shape_matrices_p =
60 auto const& solid_material =
65 for (
unsigned ip = 0; ip < n_integration_points; ip++)
67 _ip_data.emplace_back(solid_material);
69 auto const& sm_u = shape_matrices_u[ip];
70 ip_data.integration_weight =
72 sm_u.integralMeasure * sm_u.detJ;
74 ip_data.N_u_op = ShapeMatricesTypeDisplacement::template
MatrixType<
77 for (
int i = 0; i < DisplacementDim; ++i)
86 ip_data.dNdx_u = sm_u.dNdx;
88 ip_data.N_p = shape_matrices_p[ip].N;
89 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
95 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
96 typename IntegrationMethod,
int DisplacementDim>
97 std::vector<ConstitutiveVariables<DisplacementDim>>
99 IntegrationMethod, DisplacementDim>::
100 updateConstitutiveVariables(Eigen::VectorXd
const& local_x,
101 Eigen::VectorXd
const& local_x_dot,
102 double const t,
double const dt)
104 [[maybe_unused]]
auto const matrix_size =
105 gas_pressure_size + capillary_pressure_size + temperature_size +
108 assert(local_x.size() == matrix_size);
110 auto const gas_pressure =
111 local_x.template segment<gas_pressure_size>(gas_pressure_index);
113 local_x.template segment<capillary_pressure_size>(
114 capillary_pressure_index);
116 auto const temperature =
117 local_x.template segment<temperature_size>(temperature_index);
118 auto const temperature_dot =
119 local_x_dot.template segment<temperature_size>(temperature_index);
121 auto const displacement =
122 local_x.template segment<displacement_size>(displacement_index);
127 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
128 auto const& gas_phase = medium.phase(
"Gas");
129 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
130 auto const& solid_phase = medium.phase(
"Solid");
132 unsigned const n_integration_points =
133 _integration_method.getNumberOfPoints();
135 std::vector<ConstitutiveVariables<DisplacementDim>>
136 ip_constitutive_variables(n_integration_points);
138 for (
unsigned ip = 0; ip < n_integration_points; ip++)
141 auto& ip_data = _ip_data[ip];
142 auto& ip_cv = ip_constitutive_variables[ip];
144 auto const& Np = ip_data.N_p;
146 auto const& Nu = ip_data.N_u;
147 auto const& gradNu = ip_data.dNdx_u;
148 auto const& gradNp = ip_data.dNdx_p;
154 double const T = NT.dot(temperature);
155 double const pGR = Np.dot(gas_pressure);
157 double const pLR = pGR - pCap;
163 vars[
static_cast<int>(MPL::Variable::phase_pressure)] = pGR;
165 vars[
static_cast<int>(MPL::Variable::liquid_phase_pressure)] = pLR;
168 auto const K_S = ip_data.solid_material.getBulkModulus(t, pos);
171 .template value<double>(vars, pos, t, dt);
175 .template value<double>(
176 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
178 vars[
static_cast<int>(MPL::Variable::liquid_saturation)] = ip_data.s_L;
181 ip_data.k_S = MPL::formEigenTensor<DisplacementDim>(
183 .value(vars, pos, t, dt));
187 ShapeFunctionDisplacement::NPOINTS,
189 gradNu, Nu, x_coord, _is_axially_symmetric);
191 auto& eps = ip_data.eps;
192 eps.noalias() = Bu * displacement;
200 auto const sigma_total =
201 (_ip_data[ip].sigma_eff - ip_data.alpha_B *
202 (pGR - ip_data.s_L * pCap) *
203 Invariants::identity2)
206 vars[
static_cast<int>(MPL::Variable::total_stress)]
207 .emplace<SymmetricTensor>(
212 vars[
static_cast<int>(MPL::Variable::volumetric_strain)] =
213 Invariants::trace(eps);
214 vars[
static_cast<int>(MPL::Variable::equivalent_plastic_strain)] =
215 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
216 vars[
static_cast<int>(MPL::Variable::mechanical_strain)]
224 .template value<double>(vars, pos, t, dt);
226 auto const dk_rel_G_ds_L =
228 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
233 .template value<double>(vars, pos, t, dt);
235 auto const dk_rel_L_ds_L =
237 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
241 ip_data.beta_p_SR = (1. - ip_data.alpha_B) / K_S;
249 .value(vars, pos, t, dt)));
252 ip_data.beta_T_SR = Invariants::trace(ip_data.alpha_T_SR);
254 double const T_dot = NT.dot(temperature_dot);
256 dthermal_strain = ip_data.alpha_T_SR * T_dot * dt;
258 auto& eps_prev = ip_data.eps_prev;
259 auto& eps_m = ip_data.eps_m;
260 auto& eps_m_prev = ip_data.eps_m_prev;
262 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
267 auto const rho_ref_SR =
269 .template value<double>(
270 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
272 auto const lambdaSR = MPL::formEigenTensor<DisplacementDim>(
274 .value(vars, pos, t, dt));
276 double const T0 = _process_data.reference_temperature(t, pos)[0];
277 double const delta_T(T - T0);
278 ip_data.thermal_volume_strain = ip_data.beta_T_SR * delta_T;
282 .template value<double>(vars, pos, t, dt);
284 auto const phi_S_0 = 1. - phi_0;
286 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
287 auto const& m = Invariants::identity2;
288 double const div_u = m.transpose() * eps;
290 const double phi_S = phi_S_0 * (1. + ip_data.thermal_volume_strain -
291 ip_data.alpha_B * div_u);
293 const double phi_S = phi_S_0;
297 ip_data.phi = 1. - phi_S;
300 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
301 auto const rhoSR = rho_ref_SR * (1. - ip_data.thermal_volume_strain +
302 (ip_data.alpha_B - 1.) * div_u);
304 auto const rhoSR = rho_ref_SR;
307 auto const T_prev = T - T_dot * dt;
308 ip_cv.C = ip_data.updateConstitutiveRelation(vars, t, pos, dt, T_prev);
311 auto& ptm = *_process_data.phase_transition_model_;
312 ptm.computeConstitutiveVariables(&medium, vars, pos, t, dt);
314 auto const phi_L = ip_data.s_L * ip_data.phi;
315 auto const phi_G = (1. - ip_data.s_L) * ip_data.phi;
332 auto const lambdaGR = MPL::formEigenTensor<DisplacementDim>(
c.lambdaGR);
333 auto const lambdaLR = MPL::formEigenTensor<DisplacementDim>(
c.lambdaLR);
335 ip_data.lambda = phi_G * lambdaGR + phi_L * lambdaLR + phi_S * lambdaSR;
339 .template value<double>(vars, pos, t, dt);
340 ip_data.h_S = cpS * T;
341 auto const u_S = ip_data.h_S;
343 ip_data.rho_u_eff = phi_G *
c.rhoGR *
c.uG + phi_L *
c.rhoLR *
c.uL +
346 ip_data.rho_G_h_G = phi_G *
c.rhoGR *
c.hG;
347 ip_data.rho_L_h_L = phi_L *
c.rhoLR *
c.hL;
348 ip_data.rho_S_h_S = phi_S * rhoSR * ip_data.h_S;
350 ip_data.muGR =
c.muGR;
351 ip_data.muLR =
c.muLR;
353 ip_data.rhoGR =
c.rhoGR;
354 ip_data.rhoLR =
c.rhoLR;
355 ip_data.rhoSR = rhoSR;
357 ip_data.rhoCGR =
c.rhoCGR;
358 ip_data.rhoCLR =
c.rhoCLR;
359 ip_data.rhoWGR =
c.rhoWGR;
360 ip_data.rhoWLR =
c.rhoWLR;
362 ip_data.dxmCG_dpGR =
c.dxmCG_dpGR;
363 ip_data.dxmCG_dT =
c.dxmCG_dT;
364 ip_data.dxmCL_dpLR =
c.dxmCL_dpLR;
365 ip_data.dxmCL_dT =
c.dxmCL_dT;
366 ip_data.dxmWG_dpGR =
c.dxmWG_dpGR;
367 ip_data.dxmWG_dT =
c.dxmWG_dT;
368 ip_data.dxmWL_dpLR =
c.dxmWL_dpLR;
369 ip_data.dxmWL_dT =
c.dxmWL_dT;
372 ip_data.xnCG =
c.xnCG;
373 ip_data.xmCG =
c.xmCG;
374 ip_data.xmWL =
c.xmWL;
376 ip_data.diffusion_coefficient_vapour =
c.diffusion_coefficient_vapour;
377 ip_data.diffusion_coefficient_solvate =
c.diffusion_coefficient_solvate;
381 ip_data.pWGR =
c.pWGR;
386 auto const drho_LR_dT =
390 auto const drho_SR_dT =
394 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
395 * (1. - ip_data.thermal_volume_strain +
396 (ip_data.alpha_B - 1.) * div_u) -
397 rho_ref_SR * ip_data.beta_T_SR
402 auto const dphi_0_dT =
406 auto const dphi_S_0_dT = -dphi_0_dT;
407 const double dphi_S_dT = dphi_S_0_dT
408 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
409 * (1. + ip_data.thermal_volume_strain -
410 ip_data.alpha_B * div_u) +
411 phi_S_0 * ip_data.beta_T_SR
415 ip_cv.drho_u_eff_dT =
416 phi_G *
c.drho_GR_dT *
c.uG + phi_G *
c.rhoGR *
c.du_G_dT +
417 phi_L * drho_LR_dT *
c.uL + phi_L *
c.rhoLR *
c.du_L_dT +
418 phi_S * drho_SR_dT * u_S + phi_S * rhoSR * cpS +
419 dphi_S_dT * rhoSR * u_S;
427 double const ds_G_dp_cap = -ip_cv.ds_L_dp_cap;
430 double const dphi_G_dp_cap = -ip_cv.ds_L_dp_cap * ip_data.phi;
432 double const dphi_L_dp_cap = ip_cv.ds_L_dp_cap * ip_data.phi;
434 auto const dlambda_GR_dT = MPL::formEigenTensor<DisplacementDim>(
437 auto const dlambda_LR_dT = MPL::formEigenTensor<DisplacementDim>(
440 auto const dlambda_SR_dT = MPL::formEigenTensor<DisplacementDim>(
444 ip_cv.dlambda_dp_cap =
445 dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
447 ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
448 phi_S * dlambda_SR_dT + dphi_S_dT * lambdaSR;
454 double const drho_LR_dp_GR =
c.drho_LR_dp_LR;
455 double const drho_LR_dp_cap = -
c.drho_LR_dp_LR;
458 ip_cv.drho_h_eff_dp_GR =
459 phi_G *
c.drho_GR_dp_GR *
461 phi_L * drho_LR_dp_GR *
463 ip_cv.drho_h_eff_dp_cap = dphi_G_dp_cap *
c.rhoGR *
c.hG +
465 dphi_L_dp_cap *
c.rhoLR *
c.hL +
466 phi_L * drho_LR_dp_cap *
c.hL;
469 constexpr
double dphi_G_dT = 0;
470 constexpr
double dphi_L_dT = 0;
471 ip_cv.drho_h_eff_dT =
472 dphi_G_dT *
c.rhoGR *
c.hG + phi_G *
c.drho_GR_dT *
c.hG +
473 phi_G *
c.rhoGR *
c.dh_G_dT + dphi_L_dT *
c.rhoLR *
c.hL +
474 phi_L * drho_LR_dT *
c.hL + phi_L *
c.rhoLR *
c.dh_L_dT +
475 dphi_S_dT * rhoSR * ip_data.h_S + phi_S * drho_SR_dT * ip_data.h_S +
478 ip_cv.drho_u_eff_dp_GR =
480 phi_G *
c.drho_GR_dp_GR *
c.uG + phi_G *
c.rhoGR *
c.du_G_dp_GR +
482 phi_L * drho_LR_dp_GR *
c.uL + phi_L *
c.rhoLR *
c.du_L_dp_GR;
484 ip_cv.drho_u_eff_dp_cap = dphi_G_dp_cap *
c.rhoGR *
c.uG +
486 dphi_L_dp_cap *
c.rhoLR *
c.uL +
487 phi_L * drho_LR_dp_cap *
c.uL +
488 phi_L *
c.rhoLR *
c.du_L_dp_cap;
490 auto const& b = _process_data.specific_body_force;
491 auto const k_over_mu_G =
492 (ip_data.k_S * ip_data.k_rel_G / ip_data.muGR).eval();
493 auto const k_over_mu_L =
494 (ip_data.k_S * ip_data.k_rel_L / ip_data.muLR).eval();
502 ip_cv.dk_over_mu_G_dp_cap =
503 ip_data.k_S * dk_rel_G_ds_L * ip_cv.ds_L_dp_cap / ip_data.muGR;
504 ip_cv.dk_over_mu_L_dp_cap =
505 ip_data.k_S * dk_rel_L_ds_L * ip_cv.ds_L_dp_cap / ip_data.muLR;
508 k_over_mu_G *
c.rhoGR * b - k_over_mu_G * gradpGR;
510 k_over_mu_L *
c.rhoLR * b -
511 k_over_mu_L * gradpGR;
513 ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
514 c.drho_GR_dp_GR *
c.hG * w_GS +
515 c.rhoGR *
c.hG * k_over_mu_G *
c.drho_GR_dp_GR * b;
516 ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
517 -
c.rhoGR *
c.hG * k_over_mu_G -
c.rhoLR *
c.hL * k_over_mu_L;
519 ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
520 -drho_LR_dp_cap *
c.hL * w_LS -
521 c.rhoLR *
c.hL * k_over_mu_L * drho_LR_dp_cap * b;
522 ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
524 -
c.rhoLR *
c.hL * k_over_mu_L;
526 ip_cv.drho_GR_h_w_eff_dT =
527 c.drho_GR_dT *
c.hG * w_GS +
c.rhoGR *
c.dh_G_dT * w_GS +
528 drho_LR_dT *
c.hL * w_LS +
c.rhoLR *
c.dh_L_dT * w_LS;
534 double const s_L = ip_data.s_L;
535 double const s_G = 1. - ip_data.s_L;
536 double const rho_C_GR_dot = (ip_data.rhoCGR - ip_data.rhoCGR_prev) / dt;
537 double const rho_C_LR_dot = (ip_data.rhoCLR - ip_data.rhoCLR_prev) / dt;
538 double const rho_C_FR = s_G * ip_data.rhoCGR + s_L * ip_data.rhoCLR;
539 double const rho_W_FR = s_G * ip_data.rhoWGR + s_L * ip_data.rhoWLR;
541 constexpr
double drho_C_GR_dp_cap = 0;
543 s_G *
c.drho_C_GR_dp_GR / dt +
544 s_L *
c.drho_C_LR_dp_GR / dt;
545 ip_cv.dfC_3a_dp_cap =
546 ds_G_dp_cap * rho_C_GR_dot + s_G * drho_C_GR_dp_cap / dt +
547 ip_cv.ds_L_dp_cap * rho_C_LR_dot - s_L *
c.drho_C_LR_dp_LR / dt;
548 ip_cv.dfC_3a_dT = s_G *
c.drho_C_GR_dT / dt + s_L *
c.drho_C_LR_dT / dt;
550 double const drho_C_FR_dp_GR =
551 s_G *
c.drho_C_GR_dp_GR +
552 s_L *
c.drho_C_LR_dp_GR;
553 ip_cv.dfC_4_MCpG_dp_GR = drho_C_FR_dp_GR *
554 (ip_data.alpha_B - ip_data.phi) *
557 double const drho_C_FR_dT = s_G *
c.drho_C_GR_dT + s_L *
c.drho_C_LR_dT;
558 ip_cv.dfC_4_MCpG_dT =
559 drho_C_FR_dT * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR
560 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
561 - rho_C_FR * ip_data.dphi_dT * ip_data.beta_p_SR
566 drho_C_FR_dT * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_T_SR
567 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
568 + rho_C_FR * (ip_data.alpha_B - ip_data.dphi_dT) * ip_data.beta_T_SR
572 ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_data.alpha_B;
574 ip_cv.dfC_2a_dp_GR = -ip_data.phi *
c.drho_C_GR_dp_GR -
575 drho_C_FR_dp_GR * pCap *
576 (ip_data.alpha_B - ip_data.phi) *
579 double const drho_C_FR_dp_cap =
580 ds_G_dp_cap * ip_data.rhoCGR + s_G * drho_C_GR_dp_cap +
581 ip_cv.ds_L_dp_cap * ip_data.rhoCLR - s_L *
c.drho_C_LR_dp_LR;
583 ip_cv.dfC_2a_dp_cap =
584 ip_data.phi * (-
c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
585 drho_C_FR_dp_cap * pCap * (ip_data.alpha_B - ip_data.phi) *
587 rho_C_FR * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR;
590 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
591 ip_data.dphi_dT * (ip_data.rhoCLR - ip_data.rhoCGR) +
593 ip_data.phi * (
c.drho_C_LR_dT -
c.drho_C_GR_dT) -
594 drho_C_FR_dT * pCap * (ip_data.alpha_B - ip_data.phi) *
596 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
597 + rho_C_FR * pCap * ip_data.dphi_dT * ip_data.beta_p_SR
601 ip_cv.dadvection_C_dp_GR =
602 c.drho_C_GR_dp_GR * k_over_mu_G +
c.drho_C_LR_dp_GR * k_over_mu_L;
604 ip_cv.dfC_4_LCpG_dT =
605 c.drho_C_GR_dT * k_over_mu_G +
c.drho_C_LR_dT * k_over_mu_L
609 double const drho_W_FR_dp_GR =
610 s_G *
c.drho_W_GR_dp_GR +
611 s_L *
c.drho_W_LR_dp_GR;
612 double const drho_W_FR_dp_cap =
613 ds_G_dp_cap * ip_data.rhoWGR + s_G *
c.drho_W_GR_dp_cap +
614 ip_cv.ds_L_dp_cap * ip_data.rhoWLR - s_L *
c.drho_W_LR_dp_LR;
615 double const drho_W_FR_dT = s_G *
c.drho_W_GR_dT + s_L *
c.drho_W_LR_dT;
618 ip_data.phi * (
c.drho_W_LR_dp_GR -
c.drho_W_GR_dp_GR);
619 ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
620 (ip_data.alpha_B - ip_data.phi) *
622 ip_cv.dfW_2a_dp_cap =
623 ip_data.phi * (-
c.drho_W_LR_dp_LR -
c.drho_W_GR_dp_cap);
624 ip_cv.dfW_2b_dp_cap =
625 drho_W_FR_dp_cap * pCap * (ip_data.alpha_B - ip_data.phi) *
627 rho_W_FR * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR;
630 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
631 ip_data.dphi_dT * (ip_data.rhoWLR - ip_data.rhoWGR) +
633 ip_data.phi * (
c.drho_W_LR_dT -
c.drho_W_GR_dT);
635 drho_W_FR_dT * pCap * (ip_data.alpha_B - ip_data.phi) *
637 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
638 - rho_W_FR * pCap * ip_data.dphi_dT * ip_data.beta_p_SR
642 double const rho_W_GR_dot = (ip_data.rhoWGR - ip_data.rhoWGR_prev) / dt;
643 double const rho_W_LR_dot = (ip_data.rhoWLR - ip_data.rhoWLR_prev) / dt;
646 s_G *
c.drho_W_GR_dp_GR / dt +
647 s_L *
c.drho_W_LR_dp_GR / dt;
648 ip_cv.dfW_3a_dp_cap =
649 ds_G_dp_cap * rho_W_GR_dot + s_G *
c.drho_W_GR_dp_cap / dt +
650 ip_cv.ds_L_dp_cap * rho_W_LR_dot - s_L *
c.drho_W_LR_dp_LR / dt;
651 ip_cv.dfW_3a_dT = s_G *
c.drho_W_GR_dT / dt + s_L *
c.drho_W_LR_dT / dt;
653 ip_cv.dfW_4_LWpG_a_dp_GR =
c.drho_W_GR_dp_GR * k_over_mu_G
655 +
c.drho_W_LR_dp_GR * k_over_mu_L
658 ip_cv.dfW_4_LWpG_a_dp_cap = -
c.drho_W_LR_dp_LR * k_over_mu_L;
659 ip_cv.dfW_4_LWpG_a_dT =
660 c.drho_W_GR_dT * k_over_mu_G
662 +
c.drho_W_LR_dT * k_over_mu_L
667 ip_cv.dfW_4_LWpG_d_dp_GR =
668 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
669 ip_cv.dfW_4_LWpG_d_dp_cap =
670 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
671 ip_cv.dfW_4_LWpG_d_dT =
672 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
674 ip_cv.dfW_4_LWpC_a_dp_GR =
c.drho_W_LR_dp_GR * k_over_mu_L
677 ip_cv.dfW_4_LWpC_a_dp_cap = -
c.drho_W_LR_dp_LR * k_over_mu_L +
678 ip_data.rhoWLR * ip_cv.dk_over_mu_L_dp_cap;
679 ip_cv.dfW_4_LWpC_a_dT =
c.drho_W_LR_dT * k_over_mu_L
684 ip_cv.dfW_4_LWpC_d_dp_GR =
685 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
686 ip_cv.dfW_4_LWpC_d_dp_cap =
687 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
688 ip_cv.dfW_4_LWpC_d_dT =
689 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
692 return ip_constitutive_variables;
695 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
696 typename IntegrationMethod,
int DisplacementDim>
698 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
699 DisplacementDim>::setIPDataInitialConditions(std::string
const&
name,
700 double const* values,
701 int const integration_order)
703 if (integration_order !=
704 static_cast<int>(_integration_method.getIntegrationOrder()))
707 "Setting integration point initial conditions; The integration "
708 "order of the local assembler for element {:d} is different "
709 "from the integration order in the initial condition.",
713 if (
name ==
"sigma_ip")
715 if (_process_data.initial_stress !=
nullptr)
718 "Setting initial conditions for stress from integration "
719 "point data and from a parameter '{:s}' is not possible "
721 _process_data.initial_stress->name);
723 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
724 values, _ip_data, &IpData::sigma_eff);
727 if (
name ==
"saturation_ip")
732 if (
name ==
"epsilon_ip")
734 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
735 values, _ip_data, &IpData::eps);
740 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
741 typename IntegrationMethod,
int DisplacementDim>
743 IntegrationMethod, DisplacementDim>::
744 setInitialConditionsConcrete(std::vector<double>
const& local_x,
749 [[maybe_unused]]
auto const matrix_size =
750 gas_pressure_size + capillary_pressure_size + temperature_size +
753 assert(local_x.size() == matrix_size);
755 updateConstitutiveVariables(
756 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
757 Eigen::VectorXd::Zero(matrix_size), t, 0);
759 unsigned const n_integration_points =
760 _integration_method.getNumberOfPoints();
761 for (
unsigned ip = 0; ip < n_integration_points; ip++)
763 auto& ip_data = _ip_data[ip];
764 ip_data.pushBackState();
768 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
769 typename IntegrationMethod,
int DisplacementDim>
771 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
772 DisplacementDim>::assemble(
double const t,
double const dt,
773 std::vector<double>
const& local_x,
774 std::vector<double>
const& local_x_dot,
775 std::vector<double>& local_M_data,
776 std::vector<double>& local_K_data,
777 std::vector<double>& local_rhs_data)
779 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
780 temperature_size + displacement_size;
781 assert(local_x.size() == matrix_size);
783 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size>
const>(
784 local_x.data() + gas_pressure_index, gas_pressure_size);
787 Eigen::Map<VectorType<capillary_pressure_size>
const>(
788 local_x.data() + capillary_pressure_index, capillary_pressure_size);
790 auto const capillary_pressure_dot =
791 Eigen::Map<VectorType<capillary_pressure_size>
const>(
792 local_x_dot.data() + capillary_pressure_index,
793 capillary_pressure_size);
797 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
798 local_M_data, matrix_size, matrix_size);
802 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
803 local_K_data, matrix_size, matrix_size);
806 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
807 local_rhs_data, matrix_size);
813 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
814 C_index, gas_pressure_index);
815 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
816 C_index, capillary_pressure_index);
817 auto MCT = local_M.template block<C_size, temperature_size>(
818 C_index, temperature_index);
819 auto MCu = local_M.template block<C_size, displacement_size>(
820 C_index, displacement_index);
823 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
824 C_index, gas_pressure_index);
825 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
826 C_index, capillary_pressure_index);
827 auto LCT = local_K.template block<C_size, temperature_size>(
828 C_index, temperature_index);
831 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
832 W_index, gas_pressure_index);
833 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
834 W_index, capillary_pressure_index);
835 auto MWT = local_M.template block<W_size, temperature_size>(
836 W_index, temperature_index);
837 auto MWu = local_M.template block<W_size, displacement_size>(
838 W_index, displacement_index);
841 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
842 W_index, gas_pressure_index);
843 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
844 W_index, capillary_pressure_index);
845 auto LWT = local_K.template block<W_size, temperature_size>(
846 W_index, temperature_index);
849 auto MTu = local_M.template block<temperature_size, displacement_size>(
850 temperature_index, displacement_index);
853 auto KTT = local_K.template block<temperature_size, temperature_size>(
854 temperature_index, temperature_index);
857 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
858 displacement_index, gas_pressure_index);
860 local_K.template block<displacement_size, capillary_pressure_size>(
861 displacement_index, capillary_pressure_index);
864 auto fC = local_f.template segment<C_size>(C_index);
866 auto fW = local_f.template segment<W_size>(W_index);
868 auto fT = local_f.template segment<temperature_size>(temperature_index);
870 auto fU = local_f.template segment<displacement_size>(displacement_index);
875 unsigned const n_integration_points =
876 _integration_method.getNumberOfPoints();
878 updateConstitutiveVariables(
879 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
880 Eigen::Map<Eigen::VectorXd const>(local_x_dot.data(),
884 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
887 auto& ip = _ip_data[int_point];
889 auto const& Np = ip.N_p;
891 auto const& Nu = ip.N_u;
893 auto const& NpT = Np.transpose().eval();
894 auto const& NTT = NT.transpose().eval();
896 auto const& gradNp = ip.dNdx_p;
897 auto const& gradNT = gradNp;
898 auto const& gradNu = ip.dNdx_u;
900 auto const& gradNpT = gradNp.transpose().eval();
901 auto const& gradNTT = gradNT.transpose().eval();
903 auto const& Nu_op = ip.N_u_op;
904 auto const& w = ip.integration_weight;
906 auto const& m = Invariants::identity2;
908 auto const mT = m.transpose().eval();
917 ShapeFunctionDisplacement::NPOINTS,
919 gradNu, Nu, x_coord, _is_axially_symmetric);
921 auto const BuT = Bu.transpose().eval();
928 double const pCap_dot = Np.dot(capillary_pressure_dot);
929 auto& beta_T_SR = ip.beta_T_SR;
932 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
934 const double sD_G = ip.diffusion_coefficient_vapour;
935 const double sD_L = ip.diffusion_coefficient_solvate;
937 auto const D_C_G = (sD_G * I).eval();
938 auto const D_W_G = (sD_G * I).eval();
939 auto const D_C_L = (sD_L * I).eval();
940 auto const D_W_L = (sD_L * I).eval();
945 auto const s_G = 1. - s_L;
946 auto const s_L_dot = (s_L - ip.s_L_prev) / dt;
948 auto& alpha_B = ip.alpha_B;
949 auto& beta_p_SR = ip.beta_p_SR;
951 auto const& b = _process_data.specific_body_force;
957 auto const phi_G = s_G * phi;
958 auto const phi_L = s_L * phi;
959 auto const phi_S = 1. - phi;
962 auto& rho_SR = ip.rhoSR;
964 auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
967 const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
968 const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
974 auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
975 auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
976 auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
977 auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
979 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
981 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
983 auto const k_over_mu_G = (k_S * ip.k_rel_G / ip.muGR).eval();
984 auto const k_over_mu_L = (k_S * ip.k_rel_L / ip.muLR).eval();
987 k_over_mu_G * ip.rhoGR * b - k_over_mu_G * gradpGR;
990 k_over_mu_L * ip.rhoLR * b -
991 k_over_mu_L * gradpGR;
997 MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
999 NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1001 if (_process_data.apply_mass_lumping)
1007 (phi * (ip.rhoCLR - ip.rhoCGR) -
1008 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1009 s_L_dot / pCap_dot * Np * w;
1013 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1014 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1016 auto const advection_C_G = (ip.rhoCGR * k_over_mu_G).eval();
1017 auto const advection_C_L = (ip.rhoCLR * k_over_mu_L).eval();
1018 auto const diffusion_C_G_p =
1019 (phi_G * ip.rhoGR * D_C_G * ip.dxmCG_dpGR).eval();
1020 auto const diffusion_C_L_p =
1021 (phi_L * ip.rhoLR * D_C_L * ip.dxmCL_dpLR).eval();
1022 auto const diffusion_C_G_T =
1023 (phi_G * ip.rhoGR * D_C_G * ip.dxmCG_dT).eval();
1024 auto const diffusion_C_L_T =
1025 (phi_L * ip.rhoLR * D_C_L * ip.dxmCL_dT).eval();
1027 auto const advection_C = (advection_C_G + advection_C_L).eval();
1028 auto const diffusion_C_p = (diffusion_C_G_p + diffusion_C_L_p).eval();
1029 auto const diffusion_C_T = (diffusion_C_G_T + diffusion_C_L_T).eval();
1031 LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1034 gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1036 LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1038 fC.noalias() += gradNpT *
1039 (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1042 if (!_process_data.apply_mass_lumping)
1044 fC.noalias() -= NpT *
1045 (phi * (ip.rhoCLR - ip.rhoCGR) -
1046 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1051 NpT * phi * (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1057 MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1059 NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1061 if (_process_data.apply_mass_lumping)
1067 (phi * (ip.rhoWLR - ip.rhoWGR) -
1068 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1069 s_L_dot / pCap_dot * Np * w;
1073 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1075 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1077 auto const advection_W_G = (ip.rhoWGR * k_over_mu_G).eval();
1078 auto const advection_W_L = (ip.rhoWLR * k_over_mu_L).eval();
1079 auto const diffusion_W_G_p =
1080 (phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR).eval();
1081 auto const diffusion_W_L_p =
1082 (phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpLR).eval();
1083 auto const diffusion_W_G_T =
1084 (phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT).eval();
1085 auto const diffusion_W_L_T =
1086 (phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT).eval();
1088 auto const advection_W = (advection_W_G + advection_W_L).eval();
1089 auto const diffusion_W_p = (diffusion_W_G_p + diffusion_W_L_p).eval();
1090 auto const diffusion_W_T = (diffusion_W_G_T + diffusion_W_L_T).eval();
1092 LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1095 gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1097 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1099 fW.noalias() += gradNpT *
1100 (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1103 if (!_process_data.apply_mass_lumping)
1105 fW.noalias() -= NpT *
1106 (phi * (ip.rhoWLR - ip.rhoWGR) -
1107 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1112 NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1118 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1120 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1122 fT.noalias() -= NTT * rho_u_eff_dot * w;
1125 gradNTT * (ip.rhoGR * h_G * w_GS + ip.rhoLR * h_L * w_LS) * w;
1128 NTT * (ip.rhoGR * w_GS.transpose() + ip.rhoLR * w_LS.transpose()) *
1135 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1137 KUpC.noalias() += (BuT * alpha_B * s_L * m * Np) * w;
1139 fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w;
1141 if (_process_data.apply_mass_lumping)
1143 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1144 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1145 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1146 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1153 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
1154 typename IntegrationMethod,
int DisplacementDim>
1156 IntegrationMethod, DisplacementDim>::
1157 assembleWithJacobian(
double const t,
double const dt,
1158 std::vector<double>
const& local_x,
1159 std::vector<double>
const& local_xdot,
1160 const double ,
const double ,
1161 std::vector<double>& ,
1162 std::vector<double>& ,
1163 std::vector<double>& local_rhs_data,
1164 std::vector<double>& local_Jac_data)
1166 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1167 temperature_size + displacement_size;
1168 assert(local_x.size() == matrix_size);
1170 auto const temperature = Eigen::Map<VectorType<temperature_size>
const>(
1171 local_x.data() + temperature_index, temperature_size);
1173 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size>
const>(
1174 local_x.data() + gas_pressure_index, gas_pressure_size);
1177 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1178 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1180 auto const gas_pressure_dot =
1181 Eigen::Map<VectorType<gas_pressure_size>
const>(
1182 local_xdot.data() + gas_pressure_index, gas_pressure_size);
1184 auto const capillary_pressure_dot =
1185 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1186 local_xdot.data() + capillary_pressure_index,
1187 capillary_pressure_size);
1189 auto const temperature_dot = Eigen::Map<VectorType<temperature_size>
const>(
1190 local_xdot.data() + temperature_index, temperature_size);
1192 auto const displacement_dot =
1193 Eigen::Map<VectorType<displacement_size>
const>(
1194 local_xdot.data() + displacement_index, displacement_size);
1197 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1198 local_Jac_data, matrix_size, matrix_size);
1200 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1201 local_rhs_data, matrix_size);
1212 C_size, capillary_pressure_size);
1222 C_size, capillary_pressure_size);
1231 W_size, capillary_pressure_size);
1242 W_size, capillary_pressure_size);
1249 temperature_size, displacement_size);
1259 displacement_size, gas_pressure_size);
1262 displacement_size, capillary_pressure_size);
1265 auto fC = local_f.template segment<C_size>(C_index);
1267 auto fW = local_f.template segment<W_size>(W_index);
1269 auto fT = local_f.template segment<temperature_size>(temperature_index);
1271 auto fU = local_f.template segment<displacement_size>(displacement_index);
1276 unsigned const n_integration_points =
1277 _integration_method.getNumberOfPoints();
1279 auto const ip_constitutive_variables = updateConstitutiveVariables(
1280 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1281 Eigen::Map<Eigen::VectorXd const>(local_xdot.data(), local_xdot.size()),
1284 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1287 auto& ip = _ip_data[int_point];
1288 auto& ip_cv = ip_constitutive_variables[int_point];
1290 auto const& Np = ip.N_p;
1291 auto const& NT = Np;
1292 auto const& Nu = ip.N_u;
1294 auto const& NpT = Np.transpose().eval();
1295 auto const& NTT = NT.transpose().eval();
1297 auto const& gradNp = ip.dNdx_p;
1298 auto const& gradNT = gradNp;
1299 auto const& gradNu = ip.dNdx_u;
1301 auto const& gradNpT = gradNp.transpose().eval();
1302 auto const& gradNTT = gradNT.transpose().eval();
1304 auto const& Nu_op = ip.N_u_op;
1305 auto const& w = ip.integration_weight;
1307 auto const& m = Invariants::identity2;
1308 auto const mT = m.transpose().eval();
1310 auto const x_coord =
1317 ShapeFunctionDisplacement::NPOINTS,
1319 gradNu, Nu, x_coord, _is_axially_symmetric);
1321 auto const BuT = Bu.transpose().eval();
1323 double const div_u_dot = Invariants::trace(Bu * displacement_dot);
1331 double const pGR_dot = Np.dot(gas_pressure_dot);
1332 double const pCap_dot = Np.dot(capillary_pressure_dot);
1333 double const T_dot = NT.dot(temperature_dot);
1334 auto& beta_T_SR = ip.beta_T_SR;
1337 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1339 const double sD_G = ip.diffusion_coefficient_vapour;
1340 const double sD_L = ip.diffusion_coefficient_solvate;
1342 auto const D_C_G = (sD_G * I).eval();
1343 auto const D_W_G = (sD_G * I).eval();
1344 auto const D_C_L = (sD_L * I).eval();
1345 auto const D_W_L = (sD_L * I).eval();
1350 auto const s_G = 1. - s_L;
1351 auto const s_L_dot = (s_L - ip.s_L_prev) / dt;
1353 auto& alpha_B = ip.alpha_B;
1354 auto& beta_p_SR = ip.beta_p_SR;
1356 auto const& b = _process_data.specific_body_force;
1362 auto const phi_G = s_G * phi;
1363 auto const phi_L = s_L * phi;
1364 auto const phi_S = 1. - phi;
1367 auto& rho_SR = ip.rhoSR;
1369 auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
1372 const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
1373 const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
1379 auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
1380 auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
1381 auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
1382 auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
1384 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1386 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1388 auto const k_over_mu_G = (k_S * ip.k_rel_G / ip.muGR).eval();
1389 auto const k_over_mu_L = (k_S * ip.k_rel_L / ip.muLR).eval();
1392 k_over_mu_G * ip.rhoGR * b - k_over_mu_G * gradpGR;
1395 k_over_mu_L * ip.rhoLR * b -
1396 k_over_mu_L * gradpGR;
1402 MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1404 NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1406 if (_process_data.apply_mass_lumping)
1412 (phi * (ip.rhoCLR - ip.rhoCGR) -
1413 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1414 s_L_dot / pCap_dot * Np * w;
1418 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1421 .template block<C_size, temperature_size>(C_index,
1423 .noalias() += NpT * ip_cv.dfC_4_MCT_dT * T_dot * NT * w;
1425 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1428 .template block<C_size, temperature_size>(C_index,
1430 .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w;
1432 auto const advection_C_G = (ip.rhoCGR * k_over_mu_G).eval();
1433 auto const advection_C_L = (ip.rhoCLR * k_over_mu_L).eval();
1434 auto const diffusion_C_G_p =
1435 (phi_G * ip.rhoGR * D_C_G * ip.dxmCG_dpGR).eval();
1436 auto const diffusion_C_L_p =
1437 (phi_L * ip.rhoLR * D_C_L * ip.dxmCL_dpLR).eval();
1438 auto const diffusion_C_G_T =
1439 (phi_G * ip.rhoGR * D_C_G * ip.dxmCG_dT).eval();
1440 auto const diffusion_C_L_T =
1441 (phi_L * ip.rhoLR * D_C_L * ip.dxmCL_dT).eval();
1443 auto const advection_C = (advection_C_G + advection_C_L).eval();
1444 auto const diffusion_C_p = (diffusion_C_G_p + diffusion_C_L_p).eval();
1445 auto const diffusion_C_T = (diffusion_C_G_T + diffusion_C_L_T).eval();
1447 LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1450 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1452 (ip_cv.dadvection_C_dp_GR
1459 .template block<C_size, temperature_size>(C_index,
1461 .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w;
1464 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1465 NpT * ip_cv.dfC_4_MCpG_dp_GR * pGR_dot * Np * w;
1469 .template block<C_size, temperature_size>(C_index,
1471 .noalias() += NpT * ip_cv.dfC_4_MCpG_dT * pGR_dot * NT * w;
1474 gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1476 LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
1479 fC.noalias() += gradNpT *
1480 (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1483 if (!_process_data.apply_mass_lumping)
1486 auto const a = phi * (ip.rhoCLR - ip.rhoCGR) -
1487 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR;
1488 fC.noalias() -= NpT * a * s_L_dot * w;
1490 local_Jac.template block<C_size, C_size>(C_index, C_index)
1493 (ip_cv.dfC_2a_dp_GR * s_L_dot ) *
1496 local_Jac.template block<C_size, W_size>(C_index, W_index)
1499 (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.ds_L_dp_cap / dt) *
1503 .template block<C_size, temperature_size>(C_index,
1505 .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w;
1509 double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
1510 fC.noalias() -= NpT * phi * a * w;
1512 local_Jac.template block<C_size, C_size>(C_index, C_index)
1513 .noalias() += NpT * phi * ip_cv.dfC_3a_dp_GR * Np * w;
1515 local_Jac.template block<C_size, W_size>(C_index, W_index)
1516 .noalias() += NpT * phi * ip_cv.dfC_3a_dp_cap * Np * w;
1519 .template block<C_size, temperature_size>(C_index,
1523 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1526 phi * ip_cv.dfC_3a_dT) *
1533 MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1535 NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1537 if (_process_data.apply_mass_lumping)
1543 (phi * (ip.rhoWLR - ip.rhoWGR) -
1544 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1545 s_L_dot / pCap_dot * Np * w;
1549 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1551 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1553 auto const advection_W_G = (ip.rhoWGR * k_over_mu_G).eval();
1554 auto const advection_W_L = (ip.rhoWLR * k_over_mu_L).eval();
1555 auto const diffusion_W_G_p = phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR;
1556 auto const diffusion_W_L_p = phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpLR;
1557 auto const diffusion_W_G_T = phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT;
1558 auto const diffusion_W_L_T = phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT;
1560 auto const advection_W = advection_W_G + advection_W_L;
1561 auto const diffusion_W_p = diffusion_W_G_p + diffusion_W_L_p;
1562 auto const diffusion_W_T = diffusion_W_G_T + diffusion_W_L_T;
1564 LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1567 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1568 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) *
1571 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1572 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) *
1576 .template block<W_size, temperature_size>(W_index,
1578 .noalias() -= gradNpT *
1579 (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) *
1583 gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1586 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1587 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
1590 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1591 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
1595 .template block<W_size, temperature_size>(W_index,
1597 .noalias() -= gradNpT *
1598 (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
1601 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1604 fW.noalias() += gradNpT *
1605 (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1609 if (!_process_data.apply_mass_lumping)
1611 double const f = phi * (ip.rhoWLR - ip.rhoWGR);
1612 double const g = rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR;
1614 fW.noalias() -= NpT * (f - g) * s_L_dot * w;
1616 local_Jac.template block<W_size, C_size>(W_index, C_index)
1617 .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) *
1623 local_Jac.template block<W_size, W_size>(W_index, W_index)
1626 ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot +
1627 (f - g) * ip_cv.ds_L_dp_cap / dt) *
1631 .template block<W_size, temperature_size>(W_index,
1634 NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w;
1639 NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1641 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1642 NpT * phi * ip_cv.dfW_3a_dp_GR * Np * w;
1644 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1645 NpT * phi * ip_cv.dfW_3a_dp_cap * Np * w;
1648 .template block<W_size, temperature_size>(W_index,
1653 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1654 ip.dphi_dT * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
1656 phi * ip_cv.dfW_3a_dT) *
1663 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1668 .template block<temperature_size, C_size>(temperature_index,
1670 .noalias() += NTT * ip_cv.drho_h_eff_dp_GR * div_u_dot * NT * w;
1675 .template block<temperature_size, W_size>(temperature_index,
1677 .noalias() -= NTT * ip_cv.drho_h_eff_dp_cap * div_u_dot * NT * w;
1682 .template block<temperature_size, temperature_size>(
1683 temperature_index, temperature_index)
1684 .noalias() += NTT * ip_cv.drho_h_eff_dT * div_u_dot * NT * w;
1686 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1704 .template block<temperature_size, W_size>(temperature_index,
1706 .noalias() += gradNTT * ip_cv.dlambda_dp_cap * gradT * Np * w;
1710 .template block<temperature_size, temperature_size>(
1711 temperature_index, temperature_index)
1712 .noalias() += gradNTT * ip_cv.dlambda_dT * gradT * NT * w;
1715 fT.noalias() -= NTT * rho_u_eff_dot * w;
1719 .template block<temperature_size, C_size>(temperature_index,
1721 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_GR * Np * w;
1725 .template block<temperature_size, W_size>(temperature_index,
1727 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_cap * Np * w;
1732 .template block<temperature_size, temperature_size>(
1733 temperature_index, temperature_index)
1734 .noalias() += NTT * ip_cv.drho_u_eff_dT / dt * NT * w;
1738 gradNTT * (ip.rhoGR * h_G * w_GS + ip.rhoLR * h_L * w_LS) * w;
1742 .template block<temperature_size, C_size>(temperature_index,
1746 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_Npart * Np * w +
1748 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart * gradNp * w;
1752 .template block<temperature_size, W_size>(temperature_index,
1756 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_Npart) * Np * w +
1758 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart) * gradNp * w;
1762 .template block<temperature_size, temperature_size>(
1763 temperature_index, temperature_index)
1764 .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
1768 NTT * (ip.rhoGR * w_GS.transpose() + ip.rhoLR * w_LS.transpose()) *
1775 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1780 KUpC.noalias() += (BuT * alpha_B * s_L * m * Np) * w;
1785 .template block<displacement_size, W_size>(displacement_index,
1787 .noalias() += BuT * alpha_B * ip_cv.ds_L_dp_cap * pCap * m * Np * w;
1790 .template block<displacement_size, displacement_size>(
1791 displacement_index, displacement_index)
1792 .noalias() += BuT * ip_cv.C * Bu * w;
1795 fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w;
1799 .template block<displacement_size, temperature_size>(
1800 displacement_index, temperature_index)
1801 .noalias() -= BuT * (ip_cv.C * ip.alpha_T_SR) * NT * w;
1810 if (_process_data.apply_mass_lumping)
1812 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1813 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1814 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1815 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1822 LCT * temperature + MCpG * gas_pressure_dot +
1823 MCpC * capillary_pressure_dot + MCT * temperature_dot +
1824 MCu * displacement_dot;
1826 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1828 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1831 .template block<C_size, temperature_size>(C_index, temperature_index)
1832 .noalias() += LCT + MCT / dt;
1834 .template block<C_size, displacement_size>(C_index, displacement_index)
1835 .noalias() += MCu / dt;
1840 LWT * temperature + MWpG * gas_pressure_dot +
1841 MWpC * capillary_pressure_dot + MWT * temperature_dot +
1842 MWu * displacement_dot;
1844 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1846 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1849 .template block<W_size, temperature_size>(W_index, temperature_index)
1850 .noalias() += LWT + MWT / dt;
1852 .template block<W_size, displacement_size>(W_index, displacement_index)
1853 .noalias() += MWu / dt;
1857 fT.noalias() -= KTT * temperature + MTu * displacement_dot;
1860 .template block<temperature_size, temperature_size>(temperature_index,
1864 .template block<temperature_size, displacement_size>(temperature_index,
1866 .noalias() += MTu / dt;
1873 .template block<displacement_size, C_size>(displacement_index, C_index)
1876 .template block<displacement_size, W_size>(displacement_index, W_index)
1880 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
1881 typename IntegrationMethod,
int DisplacementDim>
1882 std::vector<double>
const&
1884 IntegrationMethod, DisplacementDim>::
1885 getIntPtDarcyVelocityGas(
1887 std::vector<GlobalVector*>
const& x,
1888 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
1889 std::vector<double>& cache)
const
1891 auto const num_intpts = _ip_data.size();
1893 constexpr
int process_id = 0;
1894 auto const indices =
1896 assert(!indices.empty());
1897 auto const local_x = x[process_id]->get(indices);
1901 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
1902 cache, DisplacementDim, num_intpts);
1905 Eigen::Map<
typename ShapeMatricesTypePressure::template
VectorType<
1906 gas_pressure_size>
const>(local_x.data() + gas_pressure_index,
1909 Eigen::Map<
typename ShapeMatricesTypePressure::template
VectorType<
1910 capillary_pressure_size>
const>(
1911 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1913 Eigen::Map<
typename ShapeMatricesTypePressure::template
VectorType<
1914 temperature_size>
const>(local_x.data() + temperature_index,
1917 unsigned const n_integration_points =
1918 _integration_method.getNumberOfPoints();
1923 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
1924 auto const& gas_phase = medium.phase(
"Gas");
1928 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1932 auto const& N_p = _ip_data[ip].N_p;
1936 vars[
static_cast<int>(MPL::Variable::phase_pressure)] = N_p.dot(pGR);
1942 double const dt = std::numeric_limits<double>::quiet_NaN();
1945 .template value<double>(vars, pos, t, dt);
1949 .value(vars, pos, t, dt));
1952 .template value<double>(vars, pos, t, dt);
1954 vars[
static_cast<int>(MPL::Variable::liquid_saturation)] = s_L;
1960 .template value<double>(vars, pos, t, dt);
1962 auto const k_over_mu = k_S * k_rel / mu_GR;
1966 .template value<double>(vars, pos, t, dt);
1967 auto const& b = _process_data.specific_body_force;
1970 auto const& dNdx_p = _ip_data[ip].dNdx_p;
1971 cache_matrix.col(ip).noalias() =
1972 -k_over_mu * dNdx_p * pGR + k_over_mu * rho_GR * b;
1978 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
1979 typename IntegrationMethod,
int DisplacementDim>
1980 std::vector<double>
const&
1982 IntegrationMethod, DisplacementDim>::
1983 getIntPtDarcyVelocityLiquid(
1985 std::vector<GlobalVector*>
const& x,
1986 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
1987 std::vector<double>& cache)
const
1989 auto const num_intpts = _ip_data.size();
1991 constexpr
int process_id = 0;
1992 auto const indices =
1994 assert(!indices.empty());
1995 auto const local_x = x[process_id]->get(indices);
1999 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2000 cache, DisplacementDim, num_intpts);
2003 Eigen::Map<
typename ShapeMatricesTypePressure::template
VectorType<
2004 gas_pressure_size>
const>(local_x.data() + gas_pressure_index,
2007 Eigen::Map<
typename ShapeMatricesTypePressure::template
VectorType<
2008 capillary_pressure_size>
const>(
2009 local_x.data() + capillary_pressure_index, capillary_pressure_size);
2010 auto const pLR = pGR - pCap;
2012 Eigen::Map<
typename ShapeMatricesTypePressure::template
VectorType<
2013 temperature_size>
const>(local_x.data() + temperature_index,
2016 unsigned const n_integration_points =
2017 _integration_method.getNumberOfPoints();
2022 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
2023 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
2027 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2031 auto const& N_p = _ip_data[ip].N_p;
2034 vars[
static_cast<int>(MPL::Variable::phase_pressure)] = N_p.dot(pGR);
2035 vars[
static_cast<int>(MPL::Variable::liquid_phase_pressure)] =
2042 double const dt = std::numeric_limits<double>::quiet_NaN();
2045 .template value<double>(vars, pos, t, dt);
2048 .value(vars, pos, t, dt));
2051 .template value<double>(vars, pos, t, dt);
2053 vars[
static_cast<int>(MPL::Variable::liquid_saturation)] = s_L;
2057 .template value<double>(vars, pos, t, dt);
2059 auto const k_over_mu = k_S * k_rel / mu_LR;
2061 vars[
static_cast<int>(MPL::Variable::molar_fraction)] = 1.0;
2063 auto const cCL = [&]()
2068 .template value<double>(vars, pos, t, dt);
2076 .template value<double>(vars, pos, t, dt);
2077 auto const& b = _process_data.specific_body_force;
2080 auto const& dNdx_p = _ip_data[ip].dNdx_p;
2081 cache_matrix.col(ip).noalias() =
2082 -k_over_mu * dNdx_p * pLR + k_over_mu * rho_LR * b;
2088 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
2089 typename IntegrationMethod,
int DisplacementDim>
2091 IntegrationMethod, DisplacementDim>::
2092 computeSecondaryVariableConcrete(
double const t,
double const dt,
2093 Eigen::VectorXd
const& local_x,
2094 Eigen::VectorXd
const& local_x_dot)
2096 auto const gas_pressure =
2097 local_x.template segment<gas_pressure_size>(gas_pressure_index);
2099 local_x.template segment<capillary_pressure_size>(
2100 capillary_pressure_index);
2104 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
2105 DisplacementDim>(_element, _is_axially_symmetric, gas_pressure,
2106 *_process_data.gas_pressure_interpolated);
2109 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
2111 *_process_data.capillary_pressure_interpolated);
2114 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
2115 DisplacementDim>(_element, _is_axially_symmetric, liquid_pressure,
2116 *_process_data.liquid_pressure_interpolated);
2118 auto const temperature =
2119 local_x.template segment<temperature_size>(temperature_index);
2122 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
2123 DisplacementDim>(_element, _is_axially_symmetric, temperature,
2124 *_process_data.temperature_interpolated);
2126 unsigned const n_integration_points =
2127 _integration_method.getNumberOfPoints();
2129 double saturation_avg = 0;
2131 updateConstitutiveVariables(local_x, local_x_dot, t, dt);
2133 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2135 saturation_avg += _ip_data[ip].s_L;
2137 saturation_avg /= n_integration_points;
2138 (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
virtual std::size_t getID() const final
Returns the ID of the element.
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
typename ShapeMatricesTypePressure::template MatrixType< M, N > MatrixType
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
typename ShapeMatricesTypePressure::template VectorType< N > VectorType
static const int displacement_size
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
TH2MProcessData< DisplacementDim > & _process_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
IntegrationMethod _integration_method
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
@ relative_permeability_nonwetting_phase
@ concentration
used to specify decay rate of a substance.
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
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.
std::size_t setIntegrationPointScalarData(double const *values, std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> &ip_data, MemberType member)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u