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());
120 unsigned const n_integration_points =
121 this->integration_method_.getNumberOfPoints();
123 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>
124 ip_constitutive_data(n_integration_points);
125 std::vector<ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>>
126 ip_constitutive_variables(n_integration_points);
128 for (
unsigned ip = 0; ip < n_integration_points; ip++)
130 auto& ip_data = _ip_data[ip];
131 auto& ip_cv = ip_constitutive_variables[ip];
132 auto& ip_cd = ip_constitutive_data[ip];
133 auto& ip_out = this->output_data_[ip];
134 auto& current_state = this->current_states_[ip];
135 auto& prev_state = this->prev_states_[ip];
137 auto const& Np = ip_data.N_p;
139 auto const& Nu = ip_data.N_u;
140 auto const& gradNu = ip_data.dNdx_u;
141 auto const& gradNp = ip_data.dNdx_p;
143 std::nullopt, this->element_.getID(), ip,
147 this->element_, Nu))};
153 double const T = NT.dot(temperature);
154 double const T_prev = NT.dot(temperature_prev);
155 double const pGR = Np.dot(gas_pressure);
156 double const pCap = Np.dot(capillary_pressure);
161 this->process_data_.reference_temperature(t, pos)[0]};
170 models.
biot_model.
eval({pos, t, dt}, media_data, ip_cv.biot_data);
174 ShapeFunctionDisplacement::NPOINTS,
176 gradNu, Nu, x_coord, this->is_axially_symmetric_);
178 auto& eps = ip_out.eps_data.eps;
179 eps.noalias() = Bu * displacement;
181 current_state.S_L_data, ip_cv.dS_L_dp_cap);
184 current_state.S_L_data, ip_cv.chi_S_L);
188 ip_cv.C_el_data, ip_cv.beta_p_SR);
192 {pos, t, dt}, media_data, ip_cv.C_el_data, current_state.S_L_data,
193 prev_state.S_L_data, prev_state.swelling_data,
194 current_state.swelling_data, ip_cv.swelling_data);
198 ip_cv.s_therm_exp_data);
201 T_data, ip_cv.s_therm_exp_data, ip_out.eps_data,
202 Bu * displacement_prev, prev_state.mechanical_strain_data,
203 ip_cv.swelling_data, current_state.mechanical_strain_data);
206 {pos, t, dt}, T_data, current_state.mechanical_strain_data,
207 prev_state.mechanical_strain_data, prev_state.eff_stress_data,
208 current_state.eff_stress_data, this->material_states_[ip],
209 ip_cd.s_mech_data, ip_cv.equivalent_plastic_strain_data);
212 ip_cv.biot_data, ip_cv.chi_S_L, pGR_data,
213 pCap_data, ip_cv.total_stress_data);
216 {pos, t, dt}, media_data, current_state.S_L_data, pCap_data, T_data,
217 ip_cv.total_stress_data, ip_out.eps_data,
218 ip_cv.equivalent_plastic_strain_data, ip_out.permeability_data);
221 pGR_data, pCap_data, T_data,
222 current_state.rho_W_LR);
225 {pos, t, dt}, media_data, pGR_data, pCap_data, T_data,
226 current_state.rho_W_LR, ip_out.enthalpy_data,
227 ip_out.mass_mole_fractions_data, ip_out.fluid_density_data,
228 ip_out.vapour_pressure_data, current_state.constituent_density_data,
229 ip_cv.phase_transition_data);
232 ip_out.mass_mole_fractions_data,
233 ip_cv.viscosity_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);
250 ip_cv.solid_heat_capacity_data);
253 {pos, t, dt}, media_data, T_data, ip_out.porosity_data,
254 ip_cv.porosity_d_data, current_state.S_L_data, ip_cv.dS_L_dp_cap,
255 ip_cv.thermal_conductivity_data);
257 auto const& c = ip_cv.phase_transition_data;
260 current_state.S_L_data.S_L * ip_out.porosity_data.phi;
262 (1. - current_state.S_L_data.S_L) * ip_out.porosity_data.phi;
263 double const phi_S = 1. - ip_out.porosity_data.phi;
265 ip_out.enthalpy_data.h_S = ip_cv.solid_heat_capacity_data() * T;
266 auto const u_S = ip_out.enthalpy_data.h_S;
268 current_state.internal_energy_data() =
269 phi_G * ip_out.fluid_density_data.rho_GR * c.uG +
270 phi_L * ip_out.fluid_density_data.rho_LR * c.uL +
271 phi_S * ip_out.solid_density_data.rho_SR * u_S;
273 ip_cv.effective_volumetric_enthalpy_data.rho_h_eff =
274 phi_G * ip_out.fluid_density_data.rho_GR *
275 ip_out.enthalpy_data.h_G +
276 phi_L * ip_out.fluid_density_data.rho_LR *
277 ip_out.enthalpy_data.h_L +
278 phi_S * ip_out.solid_density_data.rho_SR * ip_out.enthalpy_data.h_S;
281 auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
284 c.dxmWG_dpCap * gradpCap +
289 c.dxmWL_dpCap * gradpCap +
297 ip_out.diffusion_velocity_data.d_CG =
298 ip_out.mass_mole_fractions_data.xmCG == 0.
301 : -phi_G / ip_out.mass_mole_fractions_data.xmCG *
302 c.diffusion_coefficient_vapour * gradxmCG;
304 ip_out.diffusion_velocity_data.d_WG =
305 ip_out.mass_mole_fractions_data.xmCG == 1.
308 : -phi_G / (1 - ip_out.mass_mole_fractions_data.xmCG) *
309 c.diffusion_coefficient_vapour * gradxmWG;
311 ip_out.diffusion_velocity_data.d_CL =
315 : -phi_L / xmCL * c.diffusion_coefficient_solute * gradxmCL;
317 ip_out.diffusion_velocity_data.d_WL =
318 ip_out.mass_mole_fractions_data.xmWL == 0.
321 : -phi_L / ip_out.mass_mole_fractions_data.xmWL *
322 c.diffusion_coefficient_solute * gradxmWL;
327 ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dT =
328 phi_G * c.drho_GR_dT * c.uG +
329 phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dT +
330 phi_L * c.drho_LR_dT * c.uL +
331 phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dT +
332 phi_S * ip_cv.solid_density_d_data.drho_SR_dT * u_S +
333 phi_S * ip_out.solid_density_data.rho_SR *
334 ip_cv.solid_heat_capacity_data() -
335 ip_cv.porosity_d_data.dphi_dT * ip_out.solid_density_data.rho_SR *
340 double const ds_G_dp_cap = -ip_cv.dS_L_dp_cap();
343 double const dphi_G_dp_cap =
344 -ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
346 double const dphi_L_dp_cap =
347 ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
353 double const drho_LR_dp_GR = c.drho_LR_dp_LR;
354 double const drho_LR_dp_cap = -c.drho_LR_dp_LR;
357 ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR =
360 phi_G * c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G +
363 phi_L * drho_LR_dp_GR * ip_out.enthalpy_data.h_L;
364 ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap =
365 dphi_G_dp_cap * ip_out.fluid_density_data.rho_GR *
366 ip_out.enthalpy_data.h_G +
368 dphi_L_dp_cap * ip_out.fluid_density_data.rho_LR *
369 ip_out.enthalpy_data.h_L +
370 phi_L * drho_LR_dp_cap * ip_out.enthalpy_data.h_L;
373 constexpr double dphi_G_dT = 0;
374 constexpr double dphi_L_dT = 0;
375 ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dT =
376 dphi_G_dT * ip_out.fluid_density_data.rho_GR *
377 ip_out.enthalpy_data.h_G +
378 phi_G * c.drho_GR_dT * ip_out.enthalpy_data.h_G +
379 phi_G * ip_out.fluid_density_data.rho_GR * c.dh_G_dT +
380 dphi_L_dT * ip_out.fluid_density_data.rho_LR *
381 ip_out.enthalpy_data.h_L +
382 phi_L * c.drho_LR_dT * ip_out.enthalpy_data.h_L +
383 phi_L * ip_out.fluid_density_data.rho_LR * c.dh_L_dT -
384 ip_cv.porosity_d_data.dphi_dT * ip_out.solid_density_data.rho_SR *
385 ip_out.enthalpy_data.h_S +
386 phi_S * ip_cv.solid_density_d_data.drho_SR_dT *
387 ip_out.enthalpy_data.h_S +
388 phi_S * ip_out.solid_density_data.rho_SR *
389 ip_cv.solid_heat_capacity_data();
391 ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dp_GR =
393 phi_G * c.drho_GR_dp_GR * c.uG +
394 phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dp_GR +
396 phi_L * drho_LR_dp_GR * c.uL +
397 phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dp_GR;
399 ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dp_cap =
400 dphi_G_dp_cap * ip_out.fluid_density_data.rho_GR * c.uG +
402 dphi_L_dp_cap * ip_out.fluid_density_data.rho_LR * c.uL +
403 phi_L * drho_LR_dp_cap * c.uL +
404 phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dp_cap;
406 auto const& b = this->process_data_.specific_body_force;
408 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
409 ip_cv.viscosity_data.mu_GR;
411 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
412 ip_cv.viscosity_data.mu_LR;
422 ip_cv.dk_over_mu_G_dp_cap = ip_out.permeability_data.Ki *
423 ip_out.permeability_data.dk_rel_G_dS_L *
424 ip_cv.dS_L_dp_cap() /
425 ip_cv.viscosity_data.mu_GR;
426 ip_cv.dk_over_mu_L_dp_cap = ip_out.permeability_data.Ki *
427 ip_out.permeability_data.dk_rel_L_dS_L *
428 ip_cv.dS_L_dp_cap() /
429 ip_cv.viscosity_data.mu_LR;
431 ip_out.darcy_velocity_data.w_GS =
432 k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
433 k_over_mu_G * gradpGR;
434 ip_out.darcy_velocity_data.w_LS =
435 k_over_mu_L * gradpCap +
436 k_over_mu_L * ip_out.fluid_density_data.rho_LR * b -
437 k_over_mu_L * gradpGR;
439 ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
440 c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G *
441 ip_out.darcy_velocity_data.w_GS +
442 ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G *
443 k_over_mu_G * c.drho_GR_dp_GR * b;
444 ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
445 ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G *
447 ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
450 ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
451 -drho_LR_dp_cap * ip_out.enthalpy_data.h_L *
452 ip_out.darcy_velocity_data.w_LS -
453 ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
454 k_over_mu_L * drho_LR_dp_cap * b;
455 ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
457 ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
460 ip_cv.drho_GR_h_w_eff_dT =
461 c.drho_GR_dT * ip_out.enthalpy_data.h_G *
462 ip_out.darcy_velocity_data.w_GS +
463 ip_out.fluid_density_data.rho_GR * c.dh_G_dT *
464 ip_out.darcy_velocity_data.w_GS +
465 c.drho_LR_dT * ip_out.enthalpy_data.h_L *
466 ip_out.darcy_velocity_data.w_LS +
467 ip_out.fluid_density_data.rho_LR * c.dh_L_dT *
468 ip_out.darcy_velocity_data.w_LS;
474 double const s_L = current_state.S_L_data.S_L;
475 double const s_G = 1. - s_L;
476 double const rho_C_FR =
477 s_G * current_state.constituent_density_data.rho_C_GR +
478 s_L * current_state.constituent_density_data.rho_C_LR;
479 double const rho_W_FR =
480 s_G * current_state.constituent_density_data.rho_W_GR +
481 s_L * current_state.rho_W_LR();
483 constexpr double drho_C_GR_dp_cap = 0;
486 ip_cv.dfC_3a_dp_GR = 0.;
487 ip_cv.dfC_3a_dp_cap = 0.;
488 ip_cv.dfC_3a_dT = 0.;
492 double const rho_C_GR_dot =
493 (current_state.constituent_density_data.rho_C_GR -
494 prev_state.constituent_density_data->rho_C_GR) /
496 double const rho_C_LR_dot =
497 (current_state.constituent_density_data.rho_C_LR -
498 prev_state.constituent_density_data->rho_C_LR) /
501 s_G * c.drho_C_GR_dp_GR /
503 s_L * c.drho_C_LR_dp_GR /
505 ip_cv.dfC_3a_dp_cap = ds_G_dp_cap * rho_C_GR_dot +
506 s_G * drho_C_GR_dp_cap / dt +
507 ip_cv.dS_L_dp_cap() * rho_C_LR_dot -
508 s_L * c.drho_C_LR_dp_LR / dt;
510 s_G * c.drho_C_GR_dT / dt + s_L * c.drho_C_LR_dT / dt;
513 double const drho_C_FR_dp_GR =
516 s_G * c.drho_C_GR_dp_GR +
519 s_L * c.drho_C_LR_dp_GR;
520 ip_cv.dfC_4_MCpG_dp_GR =
521 drho_C_FR_dp_GR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
524 double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT;
525 ip_cv.dfC_4_MCpG_dT =
526 drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
528 rho_C_FR * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
531 drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
532 ip_cv.s_therm_exp_data.beta_T_SR
533#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
534 + rho_C_FR * (ip_cv.biot_data() - ip_cv.porosity_d_data.dphi_dT) *
535 ip_cv.s_therm_exp_data.beta_T_SR
539 ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_cv.biot_data();
542 -ip_out.porosity_data.phi * c.drho_C_GR_dp_GR -
543 drho_C_FR_dp_GR * pCap *
544 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
547 double const drho_C_FR_dp_cap =
548 ds_G_dp_cap * current_state.constituent_density_data.rho_C_GR +
549 s_G * drho_C_GR_dp_cap +
550 ip_cv.dS_L_dp_cap() *
551 current_state.constituent_density_data.rho_C_LR -
552 s_L * c.drho_C_LR_dp_LR;
554 ip_cv.dfC_2a_dp_cap =
555 ip_out.porosity_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
556 drho_C_FR_dp_cap * pCap *
557 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
559 rho_C_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
563 ip_cv.porosity_d_data.dphi_dT *
564 (current_state.constituent_density_data.rho_C_LR -
565 current_state.constituent_density_data.rho_C_GR) +
566 ip_out.porosity_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
567 drho_C_FR_dT * pCap *
568 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
570 rho_C_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
572 ip_cv.dadvection_C_dp_GR = c.drho_C_GR_dp_GR * k_over_mu_G
575 + c.drho_C_LR_dp_GR * k_over_mu_L;
577 ip_cv.dadvection_C_dp_cap =
579 current_state.constituent_density_data.rho_C_GR *
580 ip_cv.dk_over_mu_G_dp_cap +
581 (-c.drho_C_LR_dp_LR) * k_over_mu_L +
582 current_state.constituent_density_data.rho_C_LR *
583 ip_cv.dk_over_mu_L_dp_cap;
585 ip_cv.dfC_4_LCpG_dT =
586 c.drho_C_GR_dT * k_over_mu_G + c.drho_C_LR_dT * k_over_mu_L
590 double const drho_W_FR_dp_GR =
593 s_G * c.drho_W_GR_dp_GR +
595 s_L * c.drho_W_LR_dp_GR;
596 double const drho_W_FR_dp_cap =
597 ds_G_dp_cap * current_state.constituent_density_data.rho_W_GR +
598 s_G * c.drho_W_GR_dp_cap +
599 ip_cv.dS_L_dp_cap() * current_state.rho_W_LR() -
600 s_L * c.drho_W_LR_dp_LR;
601 double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT;
604 ip_out.porosity_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
605 ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
606 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
608 ip_cv.dfW_2a_dp_cap = ip_out.porosity_data.phi *
609 (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
610 ip_cv.dfW_2b_dp_cap =
611 drho_W_FR_dp_cap * pCap *
612 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
614 rho_W_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
618 ip_cv.porosity_d_data.dphi_dT *
619 (current_state.rho_W_LR() -
620 current_state.constituent_density_data.rho_W_GR) +
621 ip_out.porosity_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
623 drho_W_FR_dT * pCap *
624 (ip_cv.biot_data() - ip_out.porosity_data.phi) *
626 rho_W_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
630 ip_cv.dfW_3a_dp_GR = 0.;
631 ip_cv.dfW_3a_dp_cap = 0.;
632 ip_cv.dfW_3a_dT = 0.;
636 double const rho_W_GR_dot =
637 (current_state.constituent_density_data.rho_W_GR -
638 prev_state.constituent_density_data->rho_W_GR) /
640 double const rho_W_LR_dot =
641 (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
644 s_G * c.drho_W_GR_dp_GR /
646 s_L * c.drho_W_LR_dp_GR /
648 ip_cv.dfW_3a_dp_cap = ds_G_dp_cap * rho_W_GR_dot +
649 s_G * c.drho_W_GR_dp_cap / dt +
650 ip_cv.dS_L_dp_cap() * rho_W_LR_dot -
651 s_L * c.drho_W_LR_dp_LR / dt;
653 s_G * c.drho_W_GR_dT / dt + s_L * c.drho_W_LR_dT / dt;
656 ip_cv.dfW_4_LWpG_a_dp_GR = c.drho_W_GR_dp_GR * k_over_mu_G
658 + c.drho_W_LR_dp_GR * k_over_mu_L
661 ip_cv.dfW_4_LWpG_a_dp_cap =
662 c.drho_W_GR_dp_cap * k_over_mu_G +
663 current_state.constituent_density_data.rho_W_GR *
664 ip_cv.dk_over_mu_G_dp_cap +
665 -c.drho_W_LR_dp_LR * k_over_mu_L +
666 current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap;
668 ip_cv.dfW_4_LWpG_a_dT =
669 c.drho_W_GR_dT * k_over_mu_G
671 + c.drho_W_LR_dT * k_over_mu_L
676 ip_cv.dfW_4_LWpG_d_dp_GR =
677 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
678 ip_cv.dfW_4_LWpG_d_dp_cap =
679 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
680 ip_cv.dfW_4_LWpG_d_dT =
681 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
683 ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L
686 ip_cv.dfW_4_LWpC_a_dp_cap =
687 -c.drho_W_LR_dp_LR * k_over_mu_L +
688 current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap;
689 ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
694 ip_cv.dfW_4_LWpC_d_dp_GR =
695 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
696 ip_cv.dfW_4_LWpC_d_dp_cap =
697 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
698 ip_cv.dfW_4_LWpC_d_dT =
699 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
701 ip_cv.dfC_4_LCpC_a_dp_GR = c.drho_C_LR_dp_GR * k_over_mu_L
704 ip_cv.dfC_4_LCpC_a_dp_cap =
705 -c.drho_C_LR_dp_LR * k_over_mu_L +
706 current_state.constituent_density_data.rho_C_LR *
707 ip_cv.dk_over_mu_L_dp_cap;
708 ip_cv.dfC_4_LCpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
713 ip_cv.dfC_4_LCpC_d_dp_GR =
714 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
715 ip_cv.dfC_4_LCpC_d_dp_cap =
716 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
717 ip_cv.dfC_4_LCpC_d_dT =
718 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
721 return {ip_constitutive_data, ip_constitutive_variables};
921 DisplacementDim>::assemble(
double const t,
double const dt,
922 std::vector<double>
const& local_x,
923 std::vector<double>
const& local_x_prev,
924 std::vector<double>& local_M_data,
925 std::vector<double>& local_K_data,
926 std::vector<double>& local_rhs_data)
928 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
929 temperature_size + displacement_size;
930 assert(local_x.size() == matrix_size);
932 auto const capillary_pressure =
933 Eigen::Map<VectorType<capillary_pressure_size>
const>(
934 local_x.data() + capillary_pressure_index, capillary_pressure_size);
936 auto const capillary_pressure_prev =
937 Eigen::Map<VectorType<capillary_pressure_size>
const>(
938 local_x_prev.data() + capillary_pressure_index,
939 capillary_pressure_size);
943 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
944 local_M_data, matrix_size, matrix_size);
948 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
949 local_K_data, matrix_size, matrix_size);
952 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
953 local_rhs_data, matrix_size);
959 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
960 C_index, gas_pressure_index);
961 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
962 C_index, capillary_pressure_index);
963 auto MCT = local_M.template block<C_size, temperature_size>(
964 C_index, temperature_index);
965 auto MCu = local_M.template block<C_size, displacement_size>(
966 C_index, displacement_index);
969 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
970 C_index, gas_pressure_index);
971 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
972 C_index, capillary_pressure_index);
973 auto LCT = local_K.template block<C_size, temperature_size>(
974 C_index, temperature_index);
977 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
978 W_index, gas_pressure_index);
979 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
980 W_index, capillary_pressure_index);
981 auto MWT = local_M.template block<W_size, temperature_size>(
982 W_index, temperature_index);
983 auto MWu = local_M.template block<W_size, displacement_size>(
984 W_index, displacement_index);
987 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
988 W_index, gas_pressure_index);
989 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
990 W_index, capillary_pressure_index);
991 auto LWT = local_K.template block<W_size, temperature_size>(
992 W_index, temperature_index);
995 auto MTu = local_M.template block<temperature_size, displacement_size>(
996 temperature_index, displacement_index);
999 auto KTT = local_K.template block<temperature_size, temperature_size>(
1000 temperature_index, temperature_index);
1003 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1004 displacement_index, gas_pressure_index);
1006 local_K.template block<displacement_size, capillary_pressure_size>(
1007 displacement_index, capillary_pressure_index);
1010 auto fC = local_f.template segment<C_size>(C_index);
1012 auto fW = local_f.template segment<W_size>(W_index);
1014 auto fT = local_f.template segment<temperature_size>(temperature_index);
1016 auto fU = local_f.template segment<displacement_size>(displacement_index);
1018 unsigned const n_integration_points =
1019 this->integration_method_.getNumberOfPoints();
1022 this->solid_material_, *this->process_data_.phase_transition_model_};
1024 auto const [ip_constitutive_data, ip_constitutive_variables] =
1025 updateConstitutiveVariables(
1026 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1027 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1028 local_x_prev.size()),
1031 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1033 auto& ip = _ip_data[int_point];
1034 auto& ip_cv = ip_constitutive_variables[int_point];
1035 auto& ip_out = this->output_data_[int_point];
1036 auto& current_state = this->current_states_[int_point];
1037 auto const& prev_state = this->prev_states_[int_point];
1039 auto const& Np = ip.N_p;
1040 auto const& NT = Np;
1041 auto const& Nu = ip.N_u;
1043 std::nullopt, this->element_.getID(), int_point,
1047 this->element_, Nu))};
1049 auto const& NpT = Np.transpose().eval();
1050 auto const& NTT = NT.transpose().eval();
1052 auto const& gradNp = ip.dNdx_p;
1053 auto const& gradNT = gradNp;
1054 auto const& gradNu = ip.dNdx_u;
1056 auto const& gradNpT = gradNp.transpose().eval();
1057 auto const& gradNTT = gradNT.transpose().eval();
1059 auto const& w = ip.integration_weight;
1061 auto const& m = Invariants::identity2;
1063 auto const mT = m.transpose().eval();
1065 auto const x_coord =
1068 this->element_, Nu);
1072 ShapeFunctionDisplacement::NPOINTS,
1074 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1076 auto const BuT = Bu.transpose().eval();
1078 double const pCap = Np.dot(capillary_pressure);
1079 double const pCap_prev = Np.dot(capillary_pressure_prev);
1081 double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
1084 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1087 ip_cv.phase_transition_data.diffusion_coefficient_vapour;
1089 ip_cv.phase_transition_data.diffusion_coefficient_solute;
1096 auto const s_L = current_state.S_L_data.S_L;
1097 auto const s_G = 1. - s_L;
1098 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1100 auto& alpha_B = ip_cv.biot_data();
1101 auto& beta_p_SR = ip_cv.beta_p_SR();
1103 auto const& b = this->process_data_.specific_body_force;
1106 auto const phi_G = s_G * ip_out.porosity_data.phi;
1107 auto const phi_L = s_L * ip_out.porosity_data.phi;
1108 auto const phi_S = 1. - ip_out.porosity_data.phi;
1110 auto const rhoGR = ip_out.fluid_density_data.rho_GR;
1111 auto const rhoLR = ip_out.fluid_density_data.rho_LR;
1114 auto const rho = phi_G * rhoGR + phi_L * rhoLR +
1115 phi_S * ip_out.solid_density_data.rho_SR;
1118 const double rho_C_FR =
1119 s_G * current_state.constituent_density_data.rho_C_GR +
1120 s_L * current_state.constituent_density_data.rho_C_LR;
1121 const double rho_W_FR =
1122 s_G * current_state.constituent_density_data.rho_W_GR +
1123 s_L * current_state.rho_W_LR();
1125 auto const rho_C_GR_dot =
1126 (current_state.constituent_density_data.rho_C_GR -
1127 prev_state.constituent_density_data->rho_C_GR) /
1129 auto const rho_C_LR_dot =
1130 (current_state.constituent_density_data.rho_C_LR -
1131 prev_state.constituent_density_data->rho_C_LR) /
1133 auto const rho_W_GR_dot =
1134 (current_state.constituent_density_data.rho_W_GR -
1135 prev_state.constituent_density_data->rho_W_GR) /
1137 auto const rho_W_LR_dot =
1138 (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
1141 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
1142 ip_cv.viscosity_data.mu_GR;
1144 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
1145 ip_cv.viscosity_data.mu_LR;
1151 MCpG.noalias() += NpT * rho_C_FR *
1152 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1154 MCpC.noalias() -= NpT * rho_C_FR *
1155 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1158 if (this->process_data_.apply_mass_lumping)
1160 if (pCap - pCap_prev != 0.)
1164 (ip_out.porosity_data.phi *
1165 (current_state.constituent_density_data.rho_C_LR -
1166 current_state.constituent_density_data.rho_C_GR) -
1167 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1169 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1173 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
1175 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1177 using DisplacementDimMatrix =
1178 Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
1180 DisplacementDimMatrix
const advection_C_G =
1181 current_state.constituent_density_data.rho_C_GR * k_over_mu_G;
1182 DisplacementDimMatrix
const advection_C_L =
1183 current_state.constituent_density_data.rho_C_LR * k_over_mu_L;
1185 DisplacementDimMatrix
const diffusion_CGpGR =
1186 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1187 DisplacementDimMatrix
const diffusion_CLpGR =
1188 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpGR;
1190 DisplacementDimMatrix
const diffusion_CGpCap =
1191 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpCap;
1192 DisplacementDimMatrix
const diffusion_CLpCap =
1193 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpCap;
1195 DisplacementDimMatrix
const diffusion_CGT =
1196 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT;
1197 DisplacementDimMatrix
const diffusion_CLT =
1198 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT;
1200 DisplacementDimMatrix
const advection_C = advection_C_G + advection_C_L;
1201 DisplacementDimMatrix
const diffusion_C_pGR =
1202 diffusion_CGpGR + diffusion_CLpGR;
1203 DisplacementDimMatrix
const diffusion_C_pCap =
1204 diffusion_CGpCap + diffusion_CLpCap;
1206 DisplacementDimMatrix
const diffusion_C_T =
1207 diffusion_CGT + diffusion_CLT;
1210 gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w;
1213 gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w;
1215 LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1218 gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w;
1220 if (!this->process_data_.apply_mass_lumping)
1224 (ip_out.porosity_data.phi *
1225 (current_state.constituent_density_data.rho_C_LR -
1226 current_state.constituent_density_data.rho_C_GR) -
1227 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1232 fC.noalias() -= NpT * ip_out.porosity_data.phi *
1233 (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1239 MWpG.noalias() += NpT * rho_W_FR *
1240 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1242 MWpC.noalias() -= NpT * rho_W_FR *
1243 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1246 if (this->process_data_.apply_mass_lumping)
1248 if (pCap - pCap_prev != 0.)
1252 (ip_out.porosity_data.phi *
1253 (current_state.rho_W_LR() -
1254 current_state.constituent_density_data.rho_W_GR) -
1255 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1257 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1261 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
1264 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1266 DisplacementDimMatrix
const advection_W_G =
1267 current_state.constituent_density_data.rho_W_GR * k_over_mu_G;
1268 DisplacementDimMatrix
const advection_W_L =
1269 current_state.rho_W_LR() * k_over_mu_L;
1271 DisplacementDimMatrix
const diffusion_WGpGR =
1272 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1273 DisplacementDimMatrix
const diffusion_WLpGR =
1274 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpGR;
1276 DisplacementDimMatrix
const diffusion_WGpCap =
1277 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpCap;
1278 DisplacementDimMatrix
const diffusion_WLpCap =
1279 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpCap;
1281 DisplacementDimMatrix
const diffusion_WGT =
1282 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
1283 DisplacementDimMatrix
const diffusion_WLT =
1284 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
1286 DisplacementDimMatrix
const advection_W = advection_W_G + advection_W_L;
1287 DisplacementDimMatrix
const diffusion_W_pGR =
1288 diffusion_WGpGR + diffusion_WLpGR;
1289 DisplacementDimMatrix
const diffusion_W_pCap =
1290 diffusion_WGpCap + diffusion_WLpCap;
1292 DisplacementDimMatrix
const diffusion_W_T =
1293 diffusion_WGT + diffusion_WLT;
1296 gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w;
1299 gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w;
1301 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1304 gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w;
1306 if (!this->process_data_.apply_mass_lumping)
1310 (ip_out.porosity_data.phi *
1311 (current_state.rho_W_LR() -
1312 current_state.constituent_density_data.rho_W_GR) -
1313 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1318 fW.noalias() -= NpT * ip_out.porosity_data.phi *
1319 (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1325 MTu.noalias() += NTT *
1326 ip_cv.effective_volumetric_enthalpy_data.rho_h_eff *
1330 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
1332 auto const rho_u_eff_dot = (current_state.internal_energy_data() -
1333 **prev_state.internal_energy_data) /
1335 fT.noalias() -= NTT * rho_u_eff_dot * w;
1337 fT.noalias() += gradNTT *
1338 (rhoGR * ip_out.enthalpy_data.h_G *
1339 ip_out.darcy_velocity_data.w_GS +
1340 rhoLR * ip_out.enthalpy_data.h_L *
1341 ip_out.darcy_velocity_data.w_LS) *
1344 fT.noalias() += gradNTT *
1345 (current_state.constituent_density_data.rho_C_GR *
1346 ip_cv.phase_transition_data.hCG *
1347 ip_out.diffusion_velocity_data.d_CG +
1348 current_state.constituent_density_data.rho_W_GR *
1349 ip_cv.phase_transition_data.hWG *
1350 ip_out.diffusion_velocity_data.d_WG) *
1353 fT.noalias() += NTT *
1354 (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() +
1355 rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) *
1362 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1364 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
1366 fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
1367 N_u_op(Nu).transpose() * rho * b) *
1370 if (this->process_data_.apply_mass_lumping)
1372 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1373 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1374 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1375 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1386 assembleWithJacobian(
double const t,
double const dt,
1387 std::vector<double>
const& local_x,
1388 std::vector<double>
const& local_x_prev,
1389 std::vector<double>& ,
1390 std::vector<double>& ,
1391 std::vector<double>& local_rhs_data,
1392 std::vector<double>& local_Jac_data)
1394 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1395 temperature_size + displacement_size;
1396 assert(local_x.size() == matrix_size);
1398 auto const temperature = Eigen::Map<VectorType<temperature_size>
const>(
1399 local_x.data() + temperature_index, temperature_size);
1401 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size>
const>(
1402 local_x.data() + gas_pressure_index, gas_pressure_size);
1404 auto const capillary_pressure =
1405 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1406 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1408 auto const displacement = Eigen::Map<VectorType<displacement_size>
const>(
1409 local_x.data() + displacement_index, displacement_size);
1411 auto const gas_pressure_prev =
1412 Eigen::Map<VectorType<gas_pressure_size>
const>(
1413 local_x_prev.data() + gas_pressure_index, gas_pressure_size);
1415 auto const capillary_pressure_prev =
1416 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1417 local_x_prev.data() + capillary_pressure_index,
1418 capillary_pressure_size);
1420 auto const temperature_prev =
1421 Eigen::Map<VectorType<temperature_size>
const>(
1422 local_x_prev.data() + temperature_index, temperature_size);
1424 auto const displacement_prev =
1425 Eigen::Map<VectorType<displacement_size>
const>(
1426 local_x_prev.data() + displacement_index, displacement_size);
1429 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1430 local_Jac_data, matrix_size, matrix_size);
1432 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1433 local_rhs_data, matrix_size);
1444 C_size, capillary_pressure_size);
1454 C_size, capillary_pressure_size);
1463 W_size, capillary_pressure_size);
1474 W_size, capillary_pressure_size);
1481 temperature_size, displacement_size);
1491 displacement_size, gas_pressure_size);
1494 displacement_size, capillary_pressure_size);
1497 auto fC = local_f.template segment<C_size>(C_index);
1499 auto fW = local_f.template segment<W_size>(W_index);
1501 auto fT = local_f.template segment<temperature_size>(temperature_index);
1503 auto fU = local_f.template segment<displacement_size>(displacement_index);
1505 unsigned const n_integration_points =
1506 this->integration_method_.getNumberOfPoints();
1509 this->solid_material_, *this->process_data_.phase_transition_model_};
1511 auto const [ip_constitutive_data, ip_constitutive_variables] =
1512 updateConstitutiveVariables(
1513 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1514 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1515 local_x_prev.size()),
1518 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1520 auto& ip = _ip_data[int_point];
1521 auto& ip_cd = ip_constitutive_data[int_point];
1522 auto& ip_cv = ip_constitutive_variables[int_point];
1523 auto& ip_out = this->output_data_[int_point];
1524 auto& current_state = this->current_states_[int_point];
1525 auto& prev_state = this->prev_states_[int_point];
1527 auto const& Np = ip.N_p;
1528 auto const& NT = Np;
1529 auto const& Nu = ip.N_u;
1531 std::nullopt, this->element_.getID(), int_point,
1535 this->element_, Nu))};
1537 auto const& NpT = Np.transpose().eval();
1538 auto const& NTT = NT.transpose().eval();
1540 auto const& gradNp = ip.dNdx_p;
1541 auto const& gradNT = gradNp;
1542 auto const& gradNu = ip.dNdx_u;
1544 auto const& gradNpT = gradNp.transpose().eval();
1545 auto const& gradNTT = gradNT.transpose().eval();
1547 auto const& w = ip.integration_weight;
1549 auto const& m = Invariants::identity2;
1550 auto const mT = m.transpose().eval();
1552 auto const x_coord =
1555 this->element_, Nu);
1559 ShapeFunctionDisplacement::NPOINTS,
1561 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1563 auto const BuT = Bu.transpose().eval();
1565 double const div_u_dot =
1566 Invariants::trace(Bu * (displacement - displacement_prev) / dt);
1568 double const pGR = Np.dot(gas_pressure);
1569 double const pCap = Np.dot(capillary_pressure);
1570 double const T = NT.dot(temperature);
1576 double const pGR_prev = Np.dot(gas_pressure_prev);
1577 double const pCap_prev = Np.dot(capillary_pressure_prev);
1578 double const T_prev = NT.dot(temperature_prev);
1579 double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
1582 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1585 ip_cv.phase_transition_data.diffusion_coefficient_vapour;
1587 ip_cv.phase_transition_data.diffusion_coefficient_solute;
1594 auto& s_L = current_state.S_L_data.S_L;
1595 auto const s_G = 1. - s_L;
1596 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1598 auto const alpha_B = ip_cv.biot_data();
1599 auto const beta_p_SR = ip_cv.beta_p_SR();
1601 auto const& b = this->process_data_.specific_body_force;
1604 auto const phi_G = s_G * ip_out.porosity_data.phi;
1605 auto const phi_L = s_L * ip_out.porosity_data.phi;
1606 auto const phi_S = 1. - ip_out.porosity_data.phi;
1608 auto const rhoGR = ip_out.fluid_density_data.rho_GR;
1609 auto const rhoLR = ip_out.fluid_density_data.rho_LR;
1612 auto const rho = phi_G * rhoGR + phi_L * rhoLR +
1613 phi_S * ip_out.solid_density_data.rho_SR;
1616 const double rho_C_FR =
1617 s_G * current_state.constituent_density_data.rho_C_GR +
1618 s_L * current_state.constituent_density_data.rho_C_LR;
1619 const double rho_W_FR =
1620 s_G * current_state.constituent_density_data.rho_W_GR +
1621 s_L * current_state.rho_W_LR();
1623 auto const rho_C_GR_dot =
1624 (current_state.constituent_density_data.rho_C_GR -
1625 prev_state.constituent_density_data->rho_C_GR) /
1627 auto const rho_C_LR_dot =
1628 (current_state.constituent_density_data.rho_C_LR -
1629 prev_state.constituent_density_data->rho_C_LR) /
1631 auto const rho_W_GR_dot =
1632 (current_state.constituent_density_data.rho_W_GR -
1633 prev_state.constituent_density_data->rho_W_GR) /
1635 auto const rho_W_LR_dot =
1636 (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
1639 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
1640 ip_cv.viscosity_data.mu_GR;
1642 ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
1643 ip_cv.viscosity_data.mu_LR;
1649 MCpG.noalias() += NpT * rho_C_FR *
1650 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1652 MCpC.noalias() -= NpT * rho_C_FR *
1653 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1656 if (this->process_data_.apply_mass_lumping)
1658 if (pCap - pCap_prev != 0.)
1662 (ip_out.porosity_data.phi *
1663 (current_state.constituent_density_data.rho_C_LR -
1664 current_state.constituent_density_data.rho_C_GR) -
1665 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1667 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1671 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
1675 .template block<C_size, temperature_size>(C_index,
1677 .noalias() += NpT * ip_cv.dfC_4_MCT_dT * (T - T_prev) / dt * NT * w;
1679 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1682 .template block<C_size, temperature_size>(C_index,
1684 .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w;
1687 current_state.constituent_density_data.rho_C_GR * k_over_mu_G;
1689 current_state.constituent_density_data.rho_C_LR * k_over_mu_L;
1691 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1693 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpLR;
1695 -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT;
1697 -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT;
1701 diffusion_C_G_p + diffusion_C_L_p;
1703 diffusion_C_G_T + diffusion_C_L_T;
1705 LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1708 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1710 (ip_cv.dadvection_C_dp_GR
1716 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1718 (ip_cv.dadvection_C_dp_cap
1725 .template block<C_size, temperature_size>(C_index,
1727 .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w;
1730 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1731 NpT * ip_cv.dfC_4_MCpG_dp_GR * (pGR - pGR_prev) / dt * Np * w;
1735 .template block<C_size, temperature_size>(C_index,
1738 NpT * ip_cv.dfC_4_MCpG_dT * (pGR - pGR_prev) / dt * NT * w;
1741 gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1769 LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
1773 gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w;
1775 if (!this->process_data_.apply_mass_lumping)
1779 ip_out.porosity_data.phi *
1780 (current_state.constituent_density_data.rho_C_LR -
1781 current_state.constituent_density_data.rho_C_GR) -
1782 rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1784 fC.noalias() -= NpT * a * s_L_dot * w;
1786 local_Jac.template block<C_size, C_size>(C_index, C_index)
1789 (ip_cv.dfC_2a_dp_GR * s_L_dot ) *
1792 local_Jac.template block<C_size, W_size>(C_index, W_index)
1795 (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.dS_L_dp_cap() / dt) *
1799 .template block<C_size, temperature_size>(C_index,
1801 .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w;
1805 double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
1806 fC.noalias() -= NpT * ip_out.porosity_data.phi * a * w;
1808 local_Jac.template block<C_size, C_size>(C_index, C_index)
1810 NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_GR * Np * w;
1812 local_Jac.template block<C_size, W_size>(C_index, W_index)
1814 NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_cap * Np * w;
1817 .template block<C_size, temperature_size>(C_index,
1820 (ip_cv.porosity_d_data.dphi_dT * a +
1821 ip_out.porosity_data.phi * ip_cv.dfC_3a_dT) *
1828 MWpG.noalias() += NpT * rho_W_FR *
1829 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1831 MWpC.noalias() -= NpT * rho_W_FR *
1832 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
1835 if (this->process_data_.apply_mass_lumping)
1837 if (pCap - pCap_prev != 0.)
1841 (ip_out.porosity_data.phi *
1842 (current_state.rho_W_LR() -
1843 current_state.constituent_density_data.rho_W_GR) -
1844 rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
1846 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1850 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
1853 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1856 current_state.constituent_density_data.rho_W_GR * k_over_mu_G;
1858 current_state.rho_W_LR() * k_over_mu_L;
1860 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR;
1862 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpLR;
1864 phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
1866 phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
1870 diffusion_W_G_p + diffusion_W_L_p;
1872 diffusion_W_G_T + diffusion_W_L_T;
1874 LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1877 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1878 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) *
1881 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1882 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) *
1886 .template block<W_size, temperature_size>(W_index,
1888 .noalias() += gradNpT *
1889 (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) *
1893 gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1896 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1897 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
1900 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1901 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
1905 .template block<W_size, temperature_size>(W_index,
1907 .noalias() -= gradNpT *
1908 (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
1911 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1915 gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w;
1918 if (!this->process_data_.apply_mass_lumping)
1920 double const f = ip_out.porosity_data.phi *
1921 (current_state.rho_W_LR() -
1922 current_state.constituent_density_data.rho_W_GR);
1923 double const g = rho_W_FR * pCap *
1924 (alpha_B - ip_out.porosity_data.phi) * beta_p_SR;
1926 fW.noalias() -= NpT * (f - g) * s_L_dot * w;
1928 local_Jac.template block<W_size, C_size>(W_index, C_index)
1929 .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) *
1935 local_Jac.template block<W_size, W_size>(W_index, W_index)
1938 ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot +
1939 (f - g) * ip_cv.dS_L_dp_cap() / dt) *
1943 .template block<W_size, temperature_size>(W_index,
1946 NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w;
1950 fW.noalias() -= NpT * ip_out.porosity_data.phi *
1951 (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1953 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1954 NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_GR * Np * w;
1956 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1957 NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_cap * Np * w;
1960 .template block<W_size, temperature_size>(W_index,
1963 (ip_cv.porosity_d_data.dphi_dT *
1964 (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
1965 ip_out.porosity_data.phi * ip_cv.dfW_3a_dT) *
1972 MTu.noalias() += NTT *
1973 ip_cv.effective_volumetric_enthalpy_data.rho_h_eff *
1979 .template block<temperature_size, C_size>(temperature_index,
1982 NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR *
1988 .template block<temperature_size, W_size>(temperature_index,
1991 NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap *
1997 .template block<temperature_size, temperature_size>(
1998 temperature_index, temperature_index)
2000 NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dT *
2004 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
2022 .template block<temperature_size, W_size>(temperature_index,
2024 .noalias() += gradNTT *
2025 ip_cv.thermal_conductivity_data.dlambda_dp_cap *
2030 .template block<temperature_size, temperature_size>(
2031 temperature_index, temperature_index)
2032 .noalias() += gradNTT * ip_cv.thermal_conductivity_data.dlambda_dT *
2036 auto const rho_u_eff_dot = (current_state.internal_energy_data() -
2037 **prev_state.internal_energy_data) /
2039 fT.noalias() -= NTT * rho_u_eff_dot * w;
2043 .template block<temperature_size, C_size>(temperature_index,
2045 .noalias() += NpT * Np *
2046 (ip_cv.effective_volumetric_internal_energy_d_data
2052 .template block<temperature_size, W_size>(temperature_index,
2054 .noalias() += NpT * Np *
2055 (ip_cv.effective_volumetric_internal_energy_d_data
2056 .drho_u_eff_dp_cap /
2062 .template block<temperature_size, temperature_size>(
2063 temperature_index, temperature_index)
2066 (ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dT /
2070 fT.noalias() += gradNTT *
2071 (rhoGR * ip_out.enthalpy_data.h_G *
2072 ip_out.darcy_velocity_data.w_GS +
2073 rhoLR * ip_out.enthalpy_data.h_L *
2074 ip_out.darcy_velocity_data.w_LS) *
2079 .template block<temperature_size, C_size>(temperature_index,
2083 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_Npart * Np * w +
2085 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart * gradNp * w;
2089 .template block<temperature_size, W_size>(temperature_index,
2093 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_Npart) * Np * w +
2095 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart) * gradNp * w;
2099 .template block<temperature_size, temperature_size>(
2100 temperature_index, temperature_index)
2101 .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
2104 fT.noalias() += NTT *
2105 (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() +
2106 rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) *
2109 fT.noalias() += gradNTT *
2110 (current_state.constituent_density_data.rho_C_GR *
2111 ip_cv.phase_transition_data.hCG *
2112 ip_out.diffusion_velocity_data.d_CG +
2113 current_state.constituent_density_data.rho_W_GR *
2114 ip_cv.phase_transition_data.hWG *
2115 ip_out.diffusion_velocity_data.d_WG) *
2122 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
2127 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
2132 .template block<displacement_size, W_size>(displacement_index,
2134 .noalias() += BuT * alpha_B * ip_cv.chi_S_L.dchi_dS_L *
2135 ip_cv.dS_L_dp_cap() * pCap * m * Np * w;
2138 .template block<displacement_size, displacement_size>(
2139 displacement_index, displacement_index)
2140 .noalias() += BuT * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
2143 fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
2144 N_u_op(Nu).transpose() * rho * b) *
2149 .template block<displacement_size, temperature_size>(
2150 displacement_index, temperature_index)
2153 (ip_cd.s_mech_data.stiffness_tensor *
2154 ip_cv.s_therm_exp_data.solid_linear_thermal_expansivity_vector) *
2165 if (this->process_data_.apply_mass_lumping)
2167 MCpG = MCpG.colwise().sum().eval().asDiagonal();
2168 MCpC = MCpC.colwise().sum().eval().asDiagonal();
2169 MWpG = MWpG.colwise().sum().eval().asDiagonal();
2170 MWpC = MWpC.colwise().sum().eval().asDiagonal();
2176 fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
2178 MCpG * (gas_pressure - gas_pressure_prev) / dt +
2179 MCpC * (capillary_pressure - capillary_pressure_prev) / dt +
2180 MCT * (temperature - temperature_prev) / dt +
2181 MCu * (displacement - displacement_prev) / dt;
2183 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
2185 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
2188 .template block<C_size, temperature_size>(C_index, temperature_index)
2189 .noalias() += LCT + MCT / dt;
2191 .template block<C_size, displacement_size>(C_index, displacement_index)
2192 .noalias() += MCu / dt;
2196 fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
2198 MWpG * (gas_pressure - gas_pressure_prev) / dt +
2199 MWpC * (capillary_pressure - capillary_pressure_prev) / dt +
2200 MWT * (temperature - temperature_prev) / dt +
2201 MWu * (displacement - displacement_prev) / dt;
2203 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
2205 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
2208 .template block<W_size, temperature_size>(W_index, temperature_index)
2209 .noalias() += LWT + MWT / dt;
2211 .template block<W_size, displacement_size>(W_index, displacement_index)
2212 .noalias() += MWu / dt;
2217 KTT * temperature + MTu * (displacement - displacement_prev) / dt;
2220 .template block<temperature_size, temperature_size>(temperature_index,
2224 .template block<temperature_size, displacement_size>(temperature_index,
2226 .noalias() += MTu / dt;
2230 fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
2233 .template block<displacement_size, C_size>(displacement_index, C_index)
2236 .template block<displacement_size, W_size>(displacement_index, W_index)