88 updateConstitutiveVariables(
89 Eigen::VectorXd
const& local_x, Eigen::VectorXd
const& local_x_prev,
90 double const t,
double const dt,
94 [[maybe_unused]]
auto const matrix_size =
95 gas_pressure_size + capillary_pressure_size + temperature_size +
98 assert(local_x.size() == matrix_size);
100 auto const gas_pressure =
101 local_x.template segment<gas_pressure_size>(gas_pressure_index);
102 auto const capillary_pressure =
103 local_x.template segment<capillary_pressure_size>(
104 capillary_pressure_index);
106 auto const temperature =
107 local_x.template segment<temperature_size>(temperature_index);
108 auto const temperature_prev =
109 local_x_prev.template segment<temperature_size>(temperature_index);
111 auto const displacement =
112 local_x.template segment<displacement_size>(displacement_index);
113 auto const displacement_prev =
114 local_x_prev.template segment<displacement_size>(displacement_index);
117 *this->process_data_.media_map.getMedium(this->element_.getID());
118 auto const& gas_phase = medium.phase(
"Gas");
119 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
120 auto const& solid_phase = medium.phase(
"Solid");
123 unsigned const n_integration_points =
124 this->integration_method_.getNumberOfPoints();
126 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>
127 ip_constitutive_data(n_integration_points);
128 std::vector<ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>>
129 ip_constitutive_variables(n_integration_points);
131 for (
unsigned ip = 0; ip < n_integration_points; ip++)
133 auto& ip_data = _ip_data[ip];
134 auto& ip_cv = ip_constitutive_variables[ip];
135 auto& ip_cd = ip_constitutive_data[ip];
136 auto& ip_out = this->output_data_[ip];
137 auto& current_state = this->current_states_[ip];
138 auto& prev_state = this->prev_states_[ip];
140 auto const& Np = ip_data.N_p;
142 auto const& Nu = ip_data.N_u;
143 auto const& gradNu = ip_data.dNdx_u;
144 auto const& gradNp = ip_data.dNdx_p;
146 std::nullopt, this->element_.getID(), ip,
150 this->element_, Nu))};
156 double const T = NT.dot(temperature);
157 double const T_prev = NT.dot(temperature_prev);
158 double const pGR = Np.dot(gas_pressure);
159 double const pCap = Np.dot(capillary_pressure);
164 this->process_data_.reference_temperature(t, pos)[0]};
165 double const pLR = pGR - pCap;
174 models.
biot_model.
eval({pos, t, dt}, media_data, ip_cv.biot_data);
178 ShapeFunctionDisplacement::NPOINTS,
180 gradNu, Nu, x_coord, this->is_axially_symmetric_);
182 auto& eps = ip_out.eps_data.eps;
183 eps.noalias() = Bu * displacement;
185 current_state.S_L_data, ip_cv.dS_L_dp_cap);
188 current_state.S_L_data, ip_cv.chi_S_L);
192 ip_cv.C_el_data, ip_cv.beta_p_SR);
196 {pos, t, dt}, media_data, ip_cv.C_el_data, current_state.S_L_data,
197 prev_state.S_L_data, prev_state.swelling_data,
198 current_state.swelling_data, ip_cv.swelling_data);
202 ip_cv.s_therm_exp_data);
205 T_data, ip_cv.s_therm_exp_data, ip_out.eps_data,
206 Bu * displacement_prev, prev_state.mechanical_strain_data,
207 ip_cv.swelling_data, current_state.mechanical_strain_data);
210 {pos, t, dt}, T_data, current_state.mechanical_strain_data,
211 prev_state.mechanical_strain_data, prev_state.eff_stress_data,
212 current_state.eff_stress_data, this->material_states_[ip],
213 ip_cd.s_mech_data, ip_cv.equivalent_plastic_strain_data);
216 ip_cv.biot_data, ip_cv.chi_S_L, pGR_data,
217 pCap_data, ip_cv.total_stress_data);
220 {pos, t, dt}, media_data, current_state.S_L_data, pCap_data, T_data,
221 ip_cv.total_stress_data, ip_out.eps_data,
222 ip_cv.equivalent_plastic_strain_data, ip_out.permeability_data);
225 pGR_data, pCap_data, T_data,
226 current_state.rho_W_LR);
229 {pos, t, dt}, media_data, pGR_data, pCap_data, T_data,
230 current_state.rho_W_LR, ip_cv.viscosity_data, ip_out.enthalpy_data,
231 ip_out.mass_mole_fractions_data, ip_out.fluid_density_data,
232 ip_out.vapour_pressure_data, current_state.constituent_density_data,
233 ip_cv.phase_transition_data);
236#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
237 ip_cv.biot_data, ip_out.eps_data,
238 ip_cv.s_therm_exp_data,
240 ip_out.porosity_data, ip_cv.porosity_d_data);
243 {pos, t, dt}, media_data, T_data,
244#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
245 ip_cv.biot_data, ip_out.eps_data, ip_cv.s_therm_exp_data,
247 ip_out.solid_density_data, ip_cv.solid_density_d_data);
262 vars.
porosity = ip_out.porosity_data.phi;
264 auto const& c = ip_cv.phase_transition_data;
267 current_state.S_L_data.S_L * ip_out.porosity_data.phi;
269 (1. - current_state.S_L_data.S_L) * ip_out.porosity_data.phi;
270 double const phi_S = 1. - ip_out.porosity_data.phi;
273 ip_data.lambda = MaterialPropertyLib::formEigenTensor<DisplacementDim>(
277 .value(vars, pos, t, dt));
280 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
281 .template value<double>(vars, pos, t, dt);
282 ip_data.h_S = cpS * T;
283 auto const u_S = ip_data.h_S;
285 ip_data.rho_u_eff = phi_G * ip_out.fluid_density_data.rho_GR * c.uG +
286 phi_L * ip_out.fluid_density_data.rho_LR * c.uL +
287 phi_S * ip_out.solid_density_data.rho_SR * u_S;
290 phi_G * ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G;
292 phi_L * ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L;
294 phi_S * ip_out.solid_density_data.rho_SR * ip_data.h_S;
297 auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
300 c.dxmWG_dpCap * gradpCap +
305 c.dxmWL_dpCap * gradpCap +
313 ip_data.d_CG = ip_out.mass_mole_fractions_data.xmCG == 0.
316 : -phi_G / ip_out.mass_mole_fractions_data.xmCG *
317 c.diffusion_coefficient_vapour * gradxmCG;
323 : -phi_G / c.xmWG * c.diffusion_coefficient_vapour * gradxmWG;
329 : -phi_L / xmCL * c.diffusion_coefficient_solute * gradxmCL;
331 ip_data.d_WL = ip_out.mass_mole_fractions_data.xmWL == 0.
334 : -phi_L / ip_out.mass_mole_fractions_data.xmWL *
335 c.diffusion_coefficient_solute * gradxmWL;
340 auto const drho_LR_dT =
341 liquid_phase.property(MPL::PropertyType::density)
342 .template dValue<double>(vars, MPL::Variable::temperature, pos,
345 ip_cv.drho_u_eff_dT =
346 phi_G * c.drho_GR_dT * c.uG +
347 phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dT +
348 phi_L * drho_LR_dT * c.uL +
349 phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dT +
350 phi_S * ip_cv.solid_density_d_data.drho_SR_dT * u_S +
351 phi_S * ip_out.solid_density_data.rho_SR * cpS -
352 ip_cv.porosity_d_data.dphi_dT * ip_out.solid_density_data.rho_SR *
357 double const ds_G_dp_cap = -ip_cv.dS_L_dp_cap();
360 double const dphi_G_dp_cap =
361 -ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
363 double const dphi_L_dp_cap =
364 ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
366 auto const lambdaGR =
367 gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
368 ? MPL::formEigenTensor<DisplacementDim>(
370 .property(MPL::PropertyType::thermal_conductivity)
371 .value(vars, pos, t, dt))
372 : MPL::formEigenTensor<DisplacementDim>(0.);
374 auto const dlambda_GR_dT =
375 gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
376 ? MPL::formEigenTensor<DisplacementDim>(
377 gas_phase[MPL::PropertyType::thermal_conductivity].dValue(
378 vars, MPL::Variable::temperature, pos, t, dt))
379 : MPL::formEigenTensor<DisplacementDim>(0.);
381 auto const lambdaLR =
382 liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
383 ? MPL::formEigenTensor<DisplacementDim>(
385 .property(MPL::PropertyType::thermal_conductivity)
386 .value(vars, pos, t, dt))
387 : MPL::formEigenTensor<DisplacementDim>(0.);
389 auto const dlambda_LR_dT =
390 liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
391 ? MPL::formEigenTensor<DisplacementDim>(
392 liquid_phase[MPL::PropertyType::thermal_conductivity]
393 .dValue(vars, MPL::Variable::temperature, pos, t, dt))
394 : MPL::formEigenTensor<DisplacementDim>(0.);
396 auto const lambdaSR =
397 solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
398 ? MPL::formEigenTensor<DisplacementDim>(
400 .property(MPL::PropertyType::thermal_conductivity)
401 .value(vars, pos, t, dt))
402 : MPL::formEigenTensor<DisplacementDim>(0.);
404 auto const dlambda_SR_dT =
405 solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
406 ? MPL::formEigenTensor<DisplacementDim>(
407 solid_phase[MPL::PropertyType::thermal_conductivity]
408 .dValue(vars, MPL::Variable::temperature, pos, t, dt))
409 : MPL::formEigenTensor<DisplacementDim>(0.);
411 ip_cv.dlambda_dp_cap =
412 dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
414 ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
415 phi_S * dlambda_SR_dT -
416 ip_cv.porosity_d_data.dphi_dT * lambdaSR;
422 double const drho_LR_dp_GR = c.drho_LR_dp_LR;
423 double const drho_LR_dp_cap = -c.drho_LR_dp_LR;
426 ip_cv.drho_h_eff_dp_GR =
429 phi_G * c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G +
432 phi_L * drho_LR_dp_GR * ip_out.enthalpy_data.h_L;
433 ip_cv.drho_h_eff_dp_cap =
434 dphi_G_dp_cap * ip_out.fluid_density_data.rho_GR *
435 ip_out.enthalpy_data.h_G +
437 dphi_L_dp_cap * ip_out.fluid_density_data.rho_LR *
438 ip_out.enthalpy_data.h_L +
439 phi_L * drho_LR_dp_cap * ip_out.enthalpy_data.h_L;
442 constexpr double dphi_G_dT = 0;
443 constexpr double dphi_L_dT = 0;
444 ip_cv.drho_h_eff_dT =
445 dphi_G_dT * ip_out.fluid_density_data.rho_GR *
446 ip_out.enthalpy_data.h_G +
447 phi_G * c.drho_GR_dT * ip_out.enthalpy_data.h_G +
448 phi_G * ip_out.fluid_density_data.rho_GR * c.dh_G_dT +
449 dphi_L_dT * ip_out.fluid_density_data.rho_LR *
450 ip_out.enthalpy_data.h_L +
451 phi_L * drho_LR_dT * ip_out.enthalpy_data.h_L +
452 phi_L * ip_out.fluid_density_data.rho_LR * c.dh_L_dT -
453 ip_cv.porosity_d_data.dphi_dT * ip_out.solid_density_data.rho_SR *
455 phi_S * ip_cv.solid_density_d_data.drho_SR_dT * ip_data.h_S +
456 phi_S * ip_out.solid_density_data.rho_SR * cpS;
458 ip_cv.drho_u_eff_dp_GR =
460 phi_G * c.drho_GR_dp_GR * c.uG +
461 phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dp_GR +
463 phi_L * drho_LR_dp_GR * c.uL +
464 phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dp_GR;
466 ip_cv.drho_u_eff_dp_cap =
467 dphi_G_dp_cap * ip_out.fluid_density_data.rho_GR * c.uG +
469 dphi_L_dp_cap * ip_out.fluid_density_data.rho_LR * c.uL +
470 phi_L * drho_LR_dp_cap * c.uL +
471 phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dp_cap;
473 auto const& b = this->process_data_.specific_body_force;
475 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
476 ip_cv.viscosity_data.mu_GR;
478 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
479 ip_cv.viscosity_data.mu_LR;
489 ip_cv.dk_over_mu_G_dp_cap = ip_out.permeability_data.Ki *
490 ip_out.permeability_data.dk_rel_G_dS_L *
491 ip_cv.dS_L_dp_cap() /
492 ip_cv.viscosity_data.mu_GR;
493 ip_cv.dk_over_mu_L_dp_cap = ip_out.permeability_data.Ki *
494 ip_out.permeability_data.dk_rel_L_dS_L *
495 ip_cv.dS_L_dp_cap() /
496 ip_cv.viscosity_data.mu_LR;
498 ip_data.w_GS = k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
499 k_over_mu_G * gradpGR;
500 ip_data.w_LS = k_over_mu_L * gradpCap +
501 k_over_mu_L * ip_out.fluid_density_data.rho_LR * b -
502 k_over_mu_L * gradpGR;
504 ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
505 c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G * ip_data.w_GS +
506 ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G *
507 k_over_mu_G * c.drho_GR_dp_GR * b;
508 ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
509 ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G *
511 ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
514 ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
515 -drho_LR_dp_cap * ip_out.enthalpy_data.h_L * ip_data.w_LS -
516 ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
517 k_over_mu_L * drho_LR_dp_cap * b;
518 ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
520 ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
523 ip_cv.drho_GR_h_w_eff_dT =
524 c.drho_GR_dT * ip_out.enthalpy_data.h_G * ip_data.w_GS +
525 ip_out.fluid_density_data.rho_GR * c.dh_G_dT * ip_data.w_GS +
526 drho_LR_dT * ip_out.enthalpy_data.h_L * ip_data.w_LS +
527 ip_out.fluid_density_data.rho_LR * c.dh_L_dT * ip_data.w_LS;
533 double const s_L = current_state.S_L_data.S_L;
534 double const s_G = 1. - s_L;
535 double const rho_C_FR =
536 s_G * current_state.constituent_density_data.rho_C_GR +
537 s_L * current_state.constituent_density_data.rho_C_LR;
538 double const rho_W_FR =
539 s_G * current_state.constituent_density_data.rho_W_GR +
540 s_L * current_state.rho_W_LR();
542 constexpr double drho_C_GR_dp_cap = 0;
545 ip_cv.dfC_3a_dp_GR = 0.;
546 ip_cv.dfC_3a_dp_cap = 0.;
547 ip_cv.dfC_3a_dT = 0.;
551 double const rho_C_GR_dot =
552 (current_state.constituent_density_data.rho_C_GR -
553 prev_state.constituent_density_data->rho_C_GR) /
555 double const rho_C_LR_dot =
556 (current_state.constituent_density_data.rho_C_LR -
557 prev_state.constituent_density_data->rho_C_LR) /
560 s_G * c.drho_C_GR_dp_GR /
562 s_L * c.drho_C_LR_dp_GR /
564 ip_cv.dfC_3a_dp_cap = ds_G_dp_cap * rho_C_GR_dot +
565 s_G * drho_C_GR_dp_cap / dt +
566 ip_cv.dS_L_dp_cap() * rho_C_LR_dot -
567 s_L * c.drho_C_LR_dp_LR / dt;
569 s_G * c.drho_C_GR_dT / dt + s_L * c.drho_C_LR_dT / dt;
572 double const drho_C_FR_dp_GR =
575 s_G * c.drho_C_GR_dp_GR +
578 s_L * c.drho_C_LR_dp_GR;
579 ip_cv.dfC_4_MCpG_dp_GR =
580 drho_C_FR_dp_GR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
583 double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT;
584 ip_cv.dfC_4_MCpG_dT =
585 drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
587 rho_C_FR * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
590 drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
591 ip_cv.s_therm_exp_data.beta_T_SR
592#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
593 + rho_C_FR * (ip_cv.biot_data() - ip_cv.porosity_d_data.dphi_dT) *
594 ip_cv.s_therm_exp_data.beta_T_SR
598 ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_cv.biot_data();
601 -ip_out.porosity_data.phi * c.drho_C_GR_dp_GR -
602 drho_C_FR_dp_GR * pCap *
603 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
606 double const drho_C_FR_dp_cap =
607 ds_G_dp_cap * current_state.constituent_density_data.rho_C_GR +
608 s_G * drho_C_GR_dp_cap +
609 ip_cv.dS_L_dp_cap() *
610 current_state.constituent_density_data.rho_C_LR -
611 s_L * c.drho_C_LR_dp_LR;
613 ip_cv.dfC_2a_dp_cap =
614 ip_out.porosity_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
615 drho_C_FR_dp_cap * pCap *
616 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
618 rho_C_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
622 ip_cv.porosity_d_data.dphi_dT *
623 (current_state.constituent_density_data.rho_C_LR -
624 current_state.constituent_density_data.rho_C_GR) +
625 ip_out.porosity_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
626 drho_C_FR_dT * pCap *
627 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
629 rho_C_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
631 ip_cv.dadvection_C_dp_GR = c.drho_C_GR_dp_GR * k_over_mu_G
634 + c.drho_C_LR_dp_GR * k_over_mu_L;
636 ip_cv.dadvection_C_dp_cap =
638 current_state.constituent_density_data.rho_C_GR *
639 ip_cv.dk_over_mu_G_dp_cap +
640 (-c.drho_C_LR_dp_LR) * k_over_mu_L +
641 current_state.constituent_density_data.rho_C_LR *
642 ip_cv.dk_over_mu_L_dp_cap;
644 ip_cv.dfC_4_LCpG_dT =
645 c.drho_C_GR_dT * k_over_mu_G + c.drho_C_LR_dT * k_over_mu_L
649 double const drho_W_FR_dp_GR =
652 s_G * c.drho_W_GR_dp_GR +
654 s_L * c.drho_W_LR_dp_GR;
655 double const drho_W_FR_dp_cap =
656 ds_G_dp_cap * current_state.constituent_density_data.rho_W_GR +
657 s_G * c.drho_W_GR_dp_cap +
658 ip_cv.dS_L_dp_cap() * current_state.rho_W_LR() -
659 s_L * c.drho_W_LR_dp_LR;
660 double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT;
663 ip_out.porosity_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
664 ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
665 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
667 ip_cv.dfW_2a_dp_cap = ip_out.porosity_data.phi *
668 (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
669 ip_cv.dfW_2b_dp_cap =
670 drho_W_FR_dp_cap * pCap *
671 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
673 rho_W_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
677 ip_cv.porosity_d_data.dphi_dT *
678 (current_state.rho_W_LR() -
679 current_state.constituent_density_data.rho_W_GR) +
680 ip_out.porosity_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
682 drho_W_FR_dT * pCap *
683 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
685 rho_W_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
689 ip_cv.dfW_3a_dp_GR = 0.;
690 ip_cv.dfW_3a_dp_cap = 0.;
691 ip_cv.dfW_3a_dT = 0.;
695 double const rho_W_GR_dot =
696 (current_state.constituent_density_data.rho_W_GR -
697 prev_state.constituent_density_data->rho_W_GR) /
699 double const rho_W_LR_dot =
700 (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
703 s_G * c.drho_W_GR_dp_GR /
705 s_L * c.drho_W_LR_dp_GR /
707 ip_cv.dfW_3a_dp_cap = ds_G_dp_cap * rho_W_GR_dot +
708 s_G * c.drho_W_GR_dp_cap / dt +
709 ip_cv.dS_L_dp_cap() * rho_W_LR_dot -
710 s_L * c.drho_W_LR_dp_LR / dt;
712 s_G * c.drho_W_GR_dT / dt + s_L * c.drho_W_LR_dT / dt;
715 ip_cv.dfW_4_LWpG_a_dp_GR = c.drho_W_GR_dp_GR * k_over_mu_G
717 + c.drho_W_LR_dp_GR * k_over_mu_L
720 ip_cv.dfW_4_LWpG_a_dp_cap =
721 c.drho_W_GR_dp_cap * k_over_mu_G +
722 current_state.constituent_density_data.rho_W_GR *
723 ip_cv.dk_over_mu_G_dp_cap +
724 -c.drho_W_LR_dp_LR * k_over_mu_L +
725 current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap;
727 ip_cv.dfW_4_LWpG_a_dT =
728 c.drho_W_GR_dT * k_over_mu_G
730 + c.drho_W_LR_dT * k_over_mu_L
735 ip_cv.dfW_4_LWpG_d_dp_GR =
736 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
737 ip_cv.dfW_4_LWpG_d_dp_cap =
738 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
739 ip_cv.dfW_4_LWpG_d_dT =
740 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
742 ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L
745 ip_cv.dfW_4_LWpC_a_dp_cap =
746 -c.drho_W_LR_dp_LR * k_over_mu_L +
747 current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap;
748 ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
753 ip_cv.dfW_4_LWpC_d_dp_GR =
754 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
755 ip_cv.dfW_4_LWpC_d_dp_cap =
756 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
757 ip_cv.dfW_4_LWpC_d_dT =
758 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
760 ip_cv.dfC_4_LCpC_a_dp_GR = c.drho_C_LR_dp_GR * k_over_mu_L
763 ip_cv.dfC_4_LCpC_a_dp_cap =
764 -c.drho_C_LR_dp_LR * k_over_mu_L +
765 current_state.constituent_density_data.rho_C_LR *
766 ip_cv.dk_over_mu_L_dp_cap;
767 ip_cv.dfC_4_LCpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
772 ip_cv.dfC_4_LCpC_d_dp_GR =
773 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
774 ip_cv.dfC_4_LCpC_d_dp_cap =
775 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
776 ip_cv.dfC_4_LCpC_d_dT =
777 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
780 return {ip_constitutive_data, ip_constitutive_variables};
983 DisplacementDim>::assemble(
double const t,
double const dt,
984 std::vector<double>
const& local_x,
985 std::vector<double>
const& local_x_prev,
986 std::vector<double>& local_M_data,
987 std::vector<double>& local_K_data,
988 std::vector<double>& local_rhs_data)
990 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
991 temperature_size + displacement_size;
992 assert(local_x.size() == matrix_size);
994 auto const capillary_pressure =
995 Eigen::Map<VectorType<capillary_pressure_size>
const>(
996 local_x.data() + capillary_pressure_index, capillary_pressure_size);
998 auto const capillary_pressure_prev =
999 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1000 local_x_prev.data() + capillary_pressure_index,
1001 capillary_pressure_size);
1005 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1006 local_M_data, matrix_size, matrix_size);
1010 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1011 local_K_data, matrix_size, matrix_size);
1014 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1015 local_rhs_data, matrix_size);
1021 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
1022 C_index, gas_pressure_index);
1023 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
1024 C_index, capillary_pressure_index);
1025 auto MCT = local_M.template block<C_size, temperature_size>(
1026 C_index, temperature_index);
1027 auto MCu = local_M.template block<C_size, displacement_size>(
1028 C_index, displacement_index);
1031 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
1032 C_index, gas_pressure_index);
1033 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
1034 C_index, capillary_pressure_index);
1035 auto LCT = local_K.template block<C_size, temperature_size>(
1036 C_index, temperature_index);
1039 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
1040 W_index, gas_pressure_index);
1041 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
1042 W_index, capillary_pressure_index);
1043 auto MWT = local_M.template block<W_size, temperature_size>(
1044 W_index, temperature_index);
1045 auto MWu = local_M.template block<W_size, displacement_size>(
1046 W_index, displacement_index);
1049 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
1050 W_index, gas_pressure_index);
1051 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
1052 W_index, capillary_pressure_index);
1053 auto LWT = local_K.template block<W_size, temperature_size>(
1054 W_index, temperature_index);
1057 auto MTu = local_M.template block<temperature_size, displacement_size>(
1058 temperature_index, displacement_index);
1061 auto KTT = local_K.template block<temperature_size, temperature_size>(
1062 temperature_index, temperature_index);
1065 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1066 displacement_index, gas_pressure_index);
1068 local_K.template block<displacement_size, capillary_pressure_size>(
1069 displacement_index, capillary_pressure_index);
1072 auto fC = local_f.template segment<C_size>(C_index);
1074 auto fW = local_f.template segment<W_size>(W_index);
1076 auto fT = local_f.template segment<temperature_size>(temperature_index);
1078 auto fU = local_f.template segment<displacement_size>(displacement_index);
1080 unsigned const n_integration_points =
1081 this->integration_method_.getNumberOfPoints();
1084 this->solid_material_, *this->process_data_.phase_transition_model_};
1086 auto const [ip_constitutive_data, ip_constitutive_variables] =
1087 updateConstitutiveVariables(
1088 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1089 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1090 local_x_prev.size()),
1093 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1095 auto& ip = _ip_data[int_point];
1096 auto& ip_cv = ip_constitutive_variables[int_point];
1097 auto& ip_out = this->output_data_[int_point];
1098 auto& current_state = this->current_states_[int_point];
1099 auto const& prev_state = this->prev_states_[int_point];
1101 auto const& Np = ip.N_p;
1102 auto const& NT = Np;
1103 auto const& Nu = ip.N_u;
1105 std::nullopt, this->element_.getID(), int_point,
1109 this->element_, Nu))};
1111 auto const& NpT = Np.transpose().eval();
1112 auto const& NTT = NT.transpose().eval();
1114 auto const& gradNp = ip.dNdx_p;
1115 auto const& gradNT = gradNp;
1116 auto const& gradNu = ip.dNdx_u;
1118 auto const& gradNpT = gradNp.transpose().eval();
1119 auto const& gradNTT = gradNT.transpose().eval();
1121 auto const& w = ip.integration_weight;
1123 auto const& m = Invariants::identity2;
1125 auto const mT = m.transpose().eval();
1127 auto const x_coord =
1130 this->element_, Nu);
1134 ShapeFunctionDisplacement::NPOINTS,
1136 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1138 auto const BuT = Bu.transpose().eval();
1140 double const pCap = Np.dot(capillary_pressure);
1141 double const pCap_prev = Np.dot(capillary_pressure_prev);
1143 double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
1146 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1149 ip_cv.phase_transition_data.diffusion_coefficient_vapour;
1151 ip_cv.phase_transition_data.diffusion_coefficient_solute;
1158 auto const s_L = current_state.S_L_data.S_L;
1159 auto const s_G = 1. - s_L;
1160 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1162 auto& alpha_B = ip_cv.biot_data();
1163 auto& beta_p_SR = ip_cv.beta_p_SR();
1165 auto const& b = this->process_data_.specific_body_force;
1168 auto const phi_G = s_G * ip_out.porosity_data.phi;
1169 auto const phi_L = s_L * ip_out.porosity_data.phi;
1170 auto const phi_S = 1. - ip_out.porosity_data.phi;
1172 auto const rhoGR = ip_out.fluid_density_data.rho_GR;
1173 auto const rhoLR = ip_out.fluid_density_data.rho_LR;
1176 auto const rho = phi_G * rhoGR + phi_L * rhoLR +
1177 phi_S * ip_out.solid_density_data.rho_SR;
1180 const double rho_C_FR =
1181 s_G * current_state.constituent_density_data.rho_C_GR +
1182 s_L * current_state.constituent_density_data.rho_C_LR;
1183 const double rho_W_FR =
1184 s_G * current_state.constituent_density_data.rho_W_GR +
1185 s_L * current_state.rho_W_LR();
1187 auto const rho_C_GR_dot =
1188 (current_state.constituent_density_data.rho_C_GR -
1189 prev_state.constituent_density_data->rho_C_GR) /
1191 auto const rho_C_LR_dot =
1192 (current_state.constituent_density_data.rho_C_LR -
1193 prev_state.constituent_density_data->rho_C_LR) /
1195 auto const rho_W_GR_dot =
1196 (current_state.constituent_density_data.rho_W_GR -
1197 prev_state.constituent_density_data->rho_W_GR) /
1199 auto const rho_W_LR_dot =
1200 (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
1202 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1204 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1207 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
1208 ip_cv.viscosity_data.mu_GR;
1210 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
1211 ip_cv.viscosity_data.mu_LR;
1217 MCpG.noalias() += NpT * rho_C_FR *
1218 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1220 MCpC.noalias() -= NpT * rho_C_FR *
1221 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1224 if (this->process_data_.apply_mass_lumping)
1226 if (pCap - pCap_prev != 0.)
1230 (ip_out.porosity_data.phi *
1231 (current_state.constituent_density_data.rho_C_LR -
1232 current_state.constituent_density_data.rho_C_GR) -
1233 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1235 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1239 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
1241 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1243 using DisplacementDimMatrix =
1244 Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
1246 DisplacementDimMatrix
const advection_C_G =
1247 current_state.constituent_density_data.rho_C_GR * k_over_mu_G;
1248 DisplacementDimMatrix
const advection_C_L =
1249 current_state.constituent_density_data.rho_C_LR * k_over_mu_L;
1251 DisplacementDimMatrix
const diffusion_CGpGR =
1252 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1253 DisplacementDimMatrix
const diffusion_CLpGR =
1254 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpGR;
1256 DisplacementDimMatrix
const diffusion_CGpCap =
1257 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpCap;
1258 DisplacementDimMatrix
const diffusion_CLpCap =
1259 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpCap;
1261 DisplacementDimMatrix
const diffusion_CGT =
1262 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT;
1263 DisplacementDimMatrix
const diffusion_CLT =
1264 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT;
1266 DisplacementDimMatrix
const advection_C = advection_C_G + advection_C_L;
1267 DisplacementDimMatrix
const diffusion_C_pGR =
1268 diffusion_CGpGR + diffusion_CLpGR;
1269 DisplacementDimMatrix
const diffusion_C_pCap =
1270 diffusion_CGpCap + diffusion_CLpCap;
1272 DisplacementDimMatrix
const diffusion_C_T =
1273 diffusion_CGT + diffusion_CLT;
1276 gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w;
1279 gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w;
1281 LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1284 gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w;
1286 if (!this->process_data_.apply_mass_lumping)
1290 (ip_out.porosity_data.phi *
1291 (current_state.constituent_density_data.rho_C_LR -
1292 current_state.constituent_density_data.rho_C_GR) -
1293 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1298 fC.noalias() -= NpT * ip_out.porosity_data.phi *
1299 (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1305 MWpG.noalias() += NpT * rho_W_FR *
1306 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1308 MWpC.noalias() -= NpT * rho_W_FR *
1309 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1312 if (this->process_data_.apply_mass_lumping)
1314 if (pCap - pCap_prev != 0.)
1318 (ip_out.porosity_data.phi *
1319 (current_state.rho_W_LR() -
1320 current_state.constituent_density_data.rho_W_GR) -
1321 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1323 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1327 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
1330 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1332 DisplacementDimMatrix
const advection_W_G =
1333 current_state.constituent_density_data.rho_W_GR * k_over_mu_G;
1334 DisplacementDimMatrix
const advection_W_L =
1335 current_state.rho_W_LR() * k_over_mu_L;
1337 DisplacementDimMatrix
const diffusion_WGpGR =
1338 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1339 DisplacementDimMatrix
const diffusion_WLpGR =
1340 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpGR;
1342 DisplacementDimMatrix
const diffusion_WGpCap =
1343 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpCap;
1344 DisplacementDimMatrix
const diffusion_WLpCap =
1345 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpCap;
1347 DisplacementDimMatrix
const diffusion_WGT =
1348 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
1349 DisplacementDimMatrix
const diffusion_WLT =
1350 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
1352 DisplacementDimMatrix
const advection_W = advection_W_G + advection_W_L;
1353 DisplacementDimMatrix
const diffusion_W_pGR =
1354 diffusion_WGpGR + diffusion_WLpGR;
1355 DisplacementDimMatrix
const diffusion_W_pCap =
1356 diffusion_WGpCap + diffusion_WLpCap;
1358 DisplacementDimMatrix
const diffusion_W_T =
1359 diffusion_WGT + diffusion_WLT;
1362 gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w;
1365 gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w;
1367 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1370 gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w;
1372 if (!this->process_data_.apply_mass_lumping)
1376 (ip_out.porosity_data.phi *
1377 (current_state.rho_W_LR() -
1378 current_state.constituent_density_data.rho_W_GR) -
1379 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1384 fW.noalias() -= NpT * ip_out.porosity_data.phi *
1385 (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1391 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1393 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1395 fT.noalias() -= NTT * rho_u_eff_dot * w;
1397 fT.noalias() += gradNTT *
1398 (rhoGR * ip_out.enthalpy_data.h_G * ip.w_GS +
1399 rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) *
1402 fT.noalias() += gradNTT *
1403 (current_state.constituent_density_data.rho_C_GR *
1404 ip_cv.phase_transition_data.hCG * ip.d_CG +
1405 current_state.constituent_density_data.rho_W_GR *
1406 ip_cv.phase_transition_data.hWG * ip.d_WG) *
1410 NTT * (rhoGR * ip.w_GS.transpose() + rhoLR * ip.w_LS.transpose()) *
1417 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1419 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
1421 fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
1422 N_u_op(Nu).transpose() * rho * b) *
1425 if (this->process_data_.apply_mass_lumping)
1427 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1428 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1429 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1430 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1441 assembleWithJacobian(
double const t,
double const dt,
1442 std::vector<double>
const& local_x,
1443 std::vector<double>
const& local_x_prev,
1444 std::vector<double>& ,
1445 std::vector<double>& ,
1446 std::vector<double>& local_rhs_data,
1447 std::vector<double>& local_Jac_data)
1449 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1450 temperature_size + displacement_size;
1451 assert(local_x.size() == matrix_size);
1453 auto const temperature = Eigen::Map<VectorType<temperature_size>
const>(
1454 local_x.data() + temperature_index, temperature_size);
1456 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size>
const>(
1457 local_x.data() + gas_pressure_index, gas_pressure_size);
1459 auto const capillary_pressure =
1460 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1461 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1463 auto const displacement = Eigen::Map<VectorType<displacement_size>
const>(
1464 local_x.data() + displacement_index, displacement_size);
1466 auto const gas_pressure_prev =
1467 Eigen::Map<VectorType<gas_pressure_size>
const>(
1468 local_x_prev.data() + gas_pressure_index, gas_pressure_size);
1470 auto const capillary_pressure_prev =
1471 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1472 local_x_prev.data() + capillary_pressure_index,
1473 capillary_pressure_size);
1475 auto const temperature_prev =
1476 Eigen::Map<VectorType<temperature_size>
const>(
1477 local_x_prev.data() + temperature_index, temperature_size);
1479 auto const displacement_prev =
1480 Eigen::Map<VectorType<displacement_size>
const>(
1481 local_x_prev.data() + displacement_index, displacement_size);
1484 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1485 local_Jac_data, matrix_size, matrix_size);
1487 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1488 local_rhs_data, matrix_size);
1499 C_size, capillary_pressure_size);
1509 C_size, capillary_pressure_size);
1518 W_size, capillary_pressure_size);
1529 W_size, capillary_pressure_size);
1536 temperature_size, displacement_size);
1546 displacement_size, gas_pressure_size);
1549 displacement_size, capillary_pressure_size);
1552 auto fC = local_f.template segment<C_size>(C_index);
1554 auto fW = local_f.template segment<W_size>(W_index);
1556 auto fT = local_f.template segment<temperature_size>(temperature_index);
1558 auto fU = local_f.template segment<displacement_size>(displacement_index);
1560 unsigned const n_integration_points =
1561 this->integration_method_.getNumberOfPoints();
1564 this->solid_material_, *this->process_data_.phase_transition_model_};
1566 auto const [ip_constitutive_data, ip_constitutive_variables] =
1567 updateConstitutiveVariables(
1568 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1569 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1570 local_x_prev.size()),
1573 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1575 auto& ip = _ip_data[int_point];
1576 auto& ip_cd = ip_constitutive_data[int_point];
1577 auto& ip_cv = ip_constitutive_variables[int_point];
1578 auto& ip_out = this->output_data_[int_point];
1579 auto& current_state = this->current_states_[int_point];
1580 auto& prev_state = this->prev_states_[int_point];
1582 auto const& Np = ip.N_p;
1583 auto const& NT = Np;
1584 auto const& Nu = ip.N_u;
1586 std::nullopt, this->element_.getID(), int_point,
1590 this->element_, Nu))};
1592 auto const& NpT = Np.transpose().eval();
1593 auto const& NTT = NT.transpose().eval();
1595 auto const& gradNp = ip.dNdx_p;
1596 auto const& gradNT = gradNp;
1597 auto const& gradNu = ip.dNdx_u;
1599 auto const& gradNpT = gradNp.transpose().eval();
1600 auto const& gradNTT = gradNT.transpose().eval();
1602 auto const& w = ip.integration_weight;
1604 auto const& m = Invariants::identity2;
1605 auto const mT = m.transpose().eval();
1607 auto const x_coord =
1610 this->element_, Nu);
1614 ShapeFunctionDisplacement::NPOINTS,
1616 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1618 auto const BuT = Bu.transpose().eval();
1620 double const div_u_dot =
1621 Invariants::trace(Bu * (displacement - displacement_prev) / dt);
1623 double const pGR = Np.dot(gas_pressure);
1624 double const pCap = Np.dot(capillary_pressure);
1625 double const T = NT.dot(temperature);
1631 double const pGR_prev = Np.dot(gas_pressure_prev);
1632 double const pCap_prev = Np.dot(capillary_pressure_prev);
1633 double const T_prev = NT.dot(temperature_prev);
1634 double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
1637 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1640 ip_cv.phase_transition_data.diffusion_coefficient_vapour;
1642 ip_cv.phase_transition_data.diffusion_coefficient_solute;
1649 auto& s_L = current_state.S_L_data.S_L;
1650 auto const s_G = 1. - s_L;
1651 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1653 auto const alpha_B = ip_cv.biot_data();
1654 auto const beta_p_SR = ip_cv.beta_p_SR();
1656 auto const& b = this->process_data_.specific_body_force;
1659 auto const phi_G = s_G * ip_out.porosity_data.phi;
1660 auto const phi_L = s_L * ip_out.porosity_data.phi;
1661 auto const phi_S = 1. - ip_out.porosity_data.phi;
1663 auto const rhoGR = ip_out.fluid_density_data.rho_GR;
1664 auto const rhoLR = ip_out.fluid_density_data.rho_LR;
1667 auto const rho = phi_G * rhoGR + phi_L * rhoLR +
1668 phi_S * ip_out.solid_density_data.rho_SR;
1671 const double rho_C_FR =
1672 s_G * current_state.constituent_density_data.rho_C_GR +
1673 s_L * current_state.constituent_density_data.rho_C_LR;
1674 const double rho_W_FR =
1675 s_G * current_state.constituent_density_data.rho_W_GR +
1676 s_L * current_state.rho_W_LR();
1678 auto const rho_C_GR_dot =
1679 (current_state.constituent_density_data.rho_C_GR -
1680 prev_state.constituent_density_data->rho_C_GR) /
1682 auto const rho_C_LR_dot =
1683 (current_state.constituent_density_data.rho_C_LR -
1684 prev_state.constituent_density_data->rho_C_LR) /
1686 auto const rho_W_GR_dot =
1687 (current_state.constituent_density_data.rho_W_GR -
1688 prev_state.constituent_density_data->rho_W_GR) /
1690 auto const rho_W_LR_dot =
1691 (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
1693 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1695 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1698 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
1699 ip_cv.viscosity_data.mu_GR;
1701 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
1702 ip_cv.viscosity_data.mu_LR;
1708 MCpG.noalias() += NpT * rho_C_FR *
1709 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1711 MCpC.noalias() -= NpT * rho_C_FR *
1712 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1715 if (this->process_data_.apply_mass_lumping)
1717 if (pCap - pCap_prev != 0.)
1721 (ip_out.porosity_data.phi *
1722 (current_state.constituent_density_data.rho_C_LR -
1723 current_state.constituent_density_data.rho_C_GR) -
1724 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1726 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1730 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
1734 .template block<C_size, temperature_size>(C_index,
1736 .noalias() += NpT * ip_cv.dfC_4_MCT_dT * (T - T_prev) / dt * NT * w;
1738 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1741 .template block<C_size, temperature_size>(C_index,
1743 .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w;
1746 current_state.constituent_density_data.rho_C_GR * k_over_mu_G;
1748 current_state.constituent_density_data.rho_C_LR * k_over_mu_L;
1750 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1752 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpLR;
1754 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT;
1756 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT;
1760 diffusion_C_G_p + diffusion_C_L_p;
1762 diffusion_C_G_T + diffusion_C_L_T;
1764 LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1767 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1769 (ip_cv.dadvection_C_dp_GR
1775 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1777 (ip_cv.dadvection_C_dp_cap
1784 .template block<C_size, temperature_size>(C_index,
1786 .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w;
1789 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1790 NpT * ip_cv.dfC_4_MCpG_dp_GR * (pGR - pGR_prev) / dt * Np * w;
1794 .template block<C_size, temperature_size>(C_index,
1797 NpT * ip_cv.dfC_4_MCpG_dT * (pGR - pGR_prev) / dt * NT * w;
1800 gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1828 LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
1832 gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w;
1834 if (!this->process_data_.apply_mass_lumping)
1838 ip_out.porosity_data.phi *
1839 (current_state.constituent_density_data.rho_C_LR -
1840 current_state.constituent_density_data.rho_C_GR) -
1841 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1843 fC.noalias() -= NpT * a * s_L_dot * w;
1845 local_Jac.template block<C_size, C_size>(C_index, C_index)
1848 (ip_cv.dfC_2a_dp_GR * s_L_dot ) *
1851 local_Jac.template block<C_size, W_size>(C_index, W_index)
1854 (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.dS_L_dp_cap() / dt) *
1858 .template block<C_size, temperature_size>(C_index,
1860 .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w;
1864 double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
1865 fC.noalias() -= NpT * ip_out.porosity_data.phi * a * w;
1867 local_Jac.template block<C_size, C_size>(C_index, C_index)
1869 NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_GR * Np * w;
1871 local_Jac.template block<C_size, W_size>(C_index, W_index)
1873 NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_cap * Np * w;
1876 .template block<C_size, temperature_size>(C_index,
1879 (ip_cv.porosity_d_data.dphi_dT * a +
1880 ip_out.porosity_data.phi * ip_cv.dfC_3a_dT) *
1887 MWpG.noalias() += NpT * rho_W_FR *
1888 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1890 MWpC.noalias() -= NpT * rho_W_FR *
1891 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1894 if (this->process_data_.apply_mass_lumping)
1896 if (pCap - pCap_prev != 0.)
1900 (ip_out.porosity_data.phi *
1901 (current_state.rho_W_LR() -
1902 current_state.constituent_density_data.rho_W_GR) -
1903 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1905 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1909 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
1912 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1915 current_state.constituent_density_data.rho_W_GR * k_over_mu_G;
1917 current_state.rho_W_LR() * k_over_mu_L;
1919 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1921 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpLR;
1923 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
1925 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
1929 diffusion_W_G_p + diffusion_W_L_p;
1931 diffusion_W_G_T + diffusion_W_L_T;
1933 LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1936 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1937 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) *
1940 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1941 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) *
1945 .template block<W_size, temperature_size>(W_index,
1947 .noalias() += gradNpT *
1948 (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) *
1952 gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1955 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1956 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
1959 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1960 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
1964 .template block<W_size, temperature_size>(W_index,
1966 .noalias() -= gradNpT *
1967 (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
1970 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1974 gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w;
1977 if (!this->process_data_.apply_mass_lumping)
1979 double const f = ip_out.porosity_data.phi *
1980 (current_state.rho_W_LR() -
1981 current_state.constituent_density_data.rho_W_GR);
1982 double const g = rho_W_FR * pCap *
1983 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR;
1985 fW.noalias() -= NpT * (f - g) * s_L_dot * w;
1987 local_Jac.template block<W_size, C_size>(W_index, C_index)
1988 .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) *
1994 local_Jac.template block<W_size, W_size>(W_index, W_index)
1997 ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot +
1998 (f - g) * ip_cv.dS_L_dp_cap() / dt) *
2002 .template block<W_size, temperature_size>(W_index,
2005 NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w;
2009 fW.noalias() -= NpT * ip_out.porosity_data.phi *
2010 (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
2012 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
2013 NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_GR * Np * w;
2015 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
2016 NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_cap * Np * w;
2019 .template block<W_size, temperature_size>(W_index,
2022 (ip_cv.porosity_d_data.dphi_dT *
2023 (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
2024 ip_out.porosity_data.phi * ip_cv.dfW_3a_dT) *
2031 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
2036 .template block<temperature_size, C_size>(temperature_index,
2038 .noalias() += NTT * ip_cv.drho_h_eff_dp_GR * div_u_dot * NT * w;
2043 .template block<temperature_size, W_size>(temperature_index,
2045 .noalias() -= NTT * ip_cv.drho_h_eff_dp_cap * div_u_dot * NT * w;
2050 .template block<temperature_size, temperature_size>(
2051 temperature_index, temperature_index)
2052 .noalias() += NTT * ip_cv.drho_h_eff_dT * div_u_dot * NT * w;
2054 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
2072 .template block<temperature_size, W_size>(temperature_index,
2074 .noalias() += gradNTT * ip_cv.dlambda_dp_cap * gradT * Np * w;
2078 .template block<temperature_size, temperature_size>(
2079 temperature_index, temperature_index)
2080 .noalias() += gradNTT * ip_cv.dlambda_dT * gradT * NT * w;
2083 fT.noalias() -= NTT * rho_u_eff_dot * w;
2087 .template block<temperature_size, C_size>(temperature_index,
2089 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_GR * Np * w;
2093 .template block<temperature_size, W_size>(temperature_index,
2095 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_cap * Np * w;
2100 .template block<temperature_size, temperature_size>(
2101 temperature_index, temperature_index)
2102 .noalias() += NTT * ip_cv.drho_u_eff_dT / dt * NT * w;
2105 fT.noalias() += gradNTT *
2106 (rhoGR * ip_out.enthalpy_data.h_G * ip.w_GS +
2107 rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) *
2112 .template block<temperature_size, C_size>(temperature_index,
2116 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_Npart * Np * w +
2118 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart * gradNp * w;
2122 .template block<temperature_size, W_size>(temperature_index,
2126 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_Npart) * Np * w +
2128 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart) * gradNp * w;
2132 .template block<temperature_size, temperature_size>(
2133 temperature_index, temperature_index)
2134 .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
2138 NTT * (rhoGR * ip.w_GS.transpose() + rhoLR * ip.w_LS.transpose()) *
2141 fT.noalias() += gradNTT *
2142 (current_state.constituent_density_data.rho_C_GR *
2143 ip_cv.phase_transition_data.hCG * ip.d_CG +
2144 current_state.constituent_density_data.rho_W_GR *
2145 ip_cv.phase_transition_data.hWG * ip.d_WG) *
2152 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
2157 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
2162 .template block<displacement_size, W_size>(displacement_index,
2164 .noalias() += BuT * alpha_B * ip_cv.chi_S_L.dchi_dS_L *
2165 ip_cv.dS_L_dp_cap() * pCap * m * Np * w;
2168 .template block<displacement_size, displacement_size>(
2169 displacement_index, displacement_index)
2170 .noalias() += BuT * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
2173 fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
2174 N_u_op(Nu).transpose() * rho * b) *
2179 .template block<displacement_size, temperature_size>(
2180 displacement_index, temperature_index)
2183 (ip_cd.s_mech_data.stiffness_tensor *
2184 ip_cv.s_therm_exp_data.solid_linear_thermal_expansivity_vector) *
2195 if (this->process_data_.apply_mass_lumping)
2197 MCpG = MCpG.colwise().sum().eval().asDiagonal();
2198 MCpC = MCpC.colwise().sum().eval().asDiagonal();
2199 MWpG = MWpG.colwise().sum().eval().asDiagonal();
2200 MWpC = MWpC.colwise().sum().eval().asDiagonal();
2206 fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
2208 MCpG * (gas_pressure - gas_pressure_prev) / dt +
2209 MCpC * (capillary_pressure - capillary_pressure_prev) / dt +
2210 MCT * (temperature - temperature_prev) / dt +
2211 MCu * (displacement - displacement_prev) / dt;
2213 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
2215 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
2218 .template block<C_size, temperature_size>(C_index, temperature_index)
2219 .noalias() += LCT + MCT / dt;
2221 .template block<C_size, displacement_size>(C_index, displacement_index)
2222 .noalias() += MCu / dt;
2226 fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
2228 MWpG * (gas_pressure - gas_pressure_prev) / dt +
2229 MWpC * (capillary_pressure - capillary_pressure_prev) / dt +
2230 MWT * (temperature - temperature_prev) / dt +
2231 MWu * (displacement - displacement_prev) / dt;
2233 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
2235 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
2238 .template block<W_size, temperature_size>(W_index, temperature_index)
2239 .noalias() += LWT + MWT / dt;
2241 .template block<W_size, displacement_size>(W_index, displacement_index)
2242 .noalias() += MWu / dt;
2247 KTT * temperature + MTu * (displacement - displacement_prev) / dt;
2250 .template block<temperature_size, temperature_size>(temperature_index,
2254 .template block<temperature_size, displacement_size>(temperature_index,
2256 .noalias() += MTu / dt;
2260 fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
2263 .template block<displacement_size, C_size>(displacement_index, C_index)
2266 .template block<displacement_size, W_size>(displacement_index, W_index)