156 std::vector<double>
const& local_x,
157 std::vector<double>
const& local_x_prev,
158 std::vector<double>& ,
159 std::vector<double>& ,
160 std::vector<double>& local_rhs_data,
161 std::vector<double>& local_Jac_data)
163 auto const local_matrix_dim = pressure_size + temperature_size;
164 assert(local_x.size() == local_matrix_dim);
166 auto const T = Eigen::Map<
typename ShapeMatricesType::template VectorType<
167 temperature_size>
const>(local_x.data() + temperature_index,
169 auto const p_L = Eigen::Map<
170 typename ShapeMatricesType::template VectorType<pressure_size>
const>(
171 local_x.data() + pressure_index, pressure_size);
174 Eigen::Map<
typename ShapeMatricesType::template VectorType<
175 temperature_size>
const>(local_x_prev.data() + temperature_index,
177 auto const p_L_prev = Eigen::Map<
178 typename ShapeMatricesType::template VectorType<pressure_size>
const>(
179 local_x_prev.data() + pressure_index, pressure_size);
182 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
184 local_Jac_data, local_matrix_dim, local_matrix_dim);
187 typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
188 local_rhs_data, local_matrix_dim);
191 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
194 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
197 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
200 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
203 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
206 ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
209 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
211 ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
214 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
217 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
220 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
222 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
223 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
224 auto const& solid_phase = medium.phase(
"Solid");
226 medium.hasPhase(
"Gas") ? &medium.phase(
"Gas") :
nullptr;
230 unsigned const n_integration_points =
231 _integration_method.getNumberOfPoints();
232 for (
unsigned ip = 0; ip < n_integration_points; ip++)
234 auto const& w = _ip_data[ip].integration_weight;
236 auto const& N = _ip_data[ip].N;
237 auto const& dNdx = _ip_data[ip].dNdx;
240 std::nullopt, _element.getID(), ip,
251 double p_cap_prev_ip;
261 auto& S_L = _ip_data[ip].saturation;
262 auto const S_L_prev = _ip_data[ip].saturation_prev;
264 medium[MPL::PropertyType::biot_coefficient].template value<double>(
265 variables, x_position, t, dt);
267 auto& solid_elasticity = *_process_data.simplified_elasticity;
271 solid_elasticity.bulkCompressibilityFromYoungsModulus(
272 solid_phase, variables, x_position, t, dt);
273 auto const beta_SR = (1 - alpha) * beta_S;
277 liquid_phase[MPL::PropertyType::density].template value<double>(
278 variables, x_position, t, dt);
280 auto const& b = _process_data.specific_body_force;
282 double const drho_LR_dp =
283 liquid_phase[MPL::PropertyType::density].template dValue<double>(
284 variables, MPL::Variable::liquid_phase_pressure, x_position, t,
286 auto const beta_LR = drho_LR_dp / rho_LR;
288 S_L = medium[MPL::PropertyType::saturation].template value<double>(
289 variables, x_position, t, dt);
294 double const dS_L_dp_cap =
295 medium[MPL::PropertyType::saturation].template dValue<double>(
296 variables, MPL::Variable::capillary_pressure, x_position, t,
300 double const DeltaS_L_Deltap_cap =
301 (p_cap_ip == p_cap_prev_ip)
303 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
306 auto chi_S_L_prev = S_L_prev;
307 auto dchi_dS_L = 1.0;
308 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
310 auto const chi = [&medium, x_position, t, dt](
double const S_L)
314 return medium[MPL::PropertyType::bishops_effective_stress]
315 .template value<double>(variables, x_position, t, dt);
318 chi_S_L_prev = chi(S_L_prev);
320 dchi_dS_L = medium[MPL::PropertyType::bishops_effective_stress]
321 .template dValue<double>(
322 variables, MPL::Variable::liquid_saturation,
333 auto& phi = _ip_data[ip].porosity;
336 variables_prev.
porosity = _ip_data[ip].porosity_prev;
337 phi = medium[MPL::PropertyType::porosity].template value<double>(
338 variables, variables_prev, x_position, t, dt);
345 "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
346 "porosity {} in element/integration point {}/{}.",
347 alpha, phi, _element.getID(), ip);
351 medium[MPL::PropertyType::relative_permeability]
352 .template value<double>(variables, x_position, t, dt);
354 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
355 variables, x_position, t, dt);
357 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
358 medium[MPL::PropertyType::permeability].value(variables, x_position,
364 auto const K_pT_thermal_osmosis =
365 (solid_phase.hasProperty(
367 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
369 [MPL::PropertyType::thermal_osmosis_coefficient]
370 .value(variables, x_position, t, dt))
371 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
376 Eigen::Matrix<double, 3, 3>
const
377 solid_linear_thermal_expansion_coefficient =
381 .value(variables, x_position, t, dt));
384 solid_phase[MPL::PropertyType::density].template value<double>(
385 variables, x_position, t, dt);
390 laplace_p.noalias() +=
391 dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
392 laplace_T.noalias() +=
393 dNdx.transpose() * rho_LR * K_pT_thermal_osmosis * dNdx * w;
394 const double alphaB_minus_phi = alpha - phi;
395 double const a0 = alphaB_minus_phi * beta_SR;
396 double const specific_storage_a_p =
397 S_L * (phi * beta_LR + S_L * a0 +
398 chi_S_L * alpha * alpha *
399 solid_elasticity.storageContribution(
400 solid_phase, variables, x_position, t, dt));
401 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
403 double const dspecific_storage_a_p_dp_cap =
404 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0 +
406 solid_elasticity.storageContribution(
407 solid_phase, variables, x_position, t, dt) *
408 (chi_S_L + dchi_dS_L * S_L));
409 double const dspecific_storage_a_S_dp_cap =
410 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
412 storage_p_a_p.noalias() +=
413 N.transpose() * rho_LR * specific_storage_a_p * N * w;
415 storage_p_a_S.noalias() -= N.transpose() * rho_LR *
416 specific_storage_a_S * DeltaS_L_Deltap_cap *
420 .template block<pressure_size, pressure_size>(pressure_index,
422 .noalias() += N.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
423 rho_LR * dspecific_storage_a_p_dp_cap * N * w;
425 storage_p_a_S_Jpp.noalias() -=
426 N.transpose() * rho_LR *
427 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
428 specific_storage_a_S * dS_L_dp_cap) /
431 double const dk_rel_dS_L =
432 medium[MPL::PropertyType::relative_permeability]
433 .template dValue<double>(variables,
434 MPL::Variable::liquid_saturation,
438 .template block<pressure_size, pressure_size>(pressure_index,
440 .noalias() += dNdx.transpose() * rho_Ki_over_mu * grad_p_cap *
441 dk_rel_dS_L * dS_L_dp_cap * N * w;
444 .template block<pressure_size, pressure_size>(pressure_index,
446 .noalias() += dNdx.transpose() * rho_LR * rho_Ki_over_mu * b *
447 dk_rel_dS_L * dS_L_dp_cap * N * w;
449 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
450 dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
455 double const fluid_volumetric_thermal_expansion_coefficient =
458 const double eff_thermal_expansion =
459 S_L * (alphaB_minus_phi *
460 solid_linear_thermal_expansion_coefficient.trace() +
461 phi * fluid_volumetric_thermal_expansion_coefficient +
462 alpha * solid_elasticity.thermalExpansivityContribution(
463 solid_linear_thermal_expansion_coefficient,
464 solid_phase, variables, x_position, t, dt));
466 N.transpose() * rho_LR * eff_thermal_expansion * N * w;
472 auto const specific_heat_capacity_fluid =
474 .template value<double>(variables, x_position, t, dt);
476 auto const specific_heat_capacity_solid =
479 .template value<double>(variables, x_position, t, dt);
483 (rho_SR * specific_heat_capacity_solid * (1 - phi) +
484 (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
487 auto const thermal_conductivity =
488 MaterialPropertyLib::formEigenTensor<GlobalDim>(
490 thermal_conductivity]
491 .value(variables, x_position, t, dt));
494 -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b) -
495 K_pT_thermal_osmosis * dNdx * T);
497 K_TT.noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
498 N.transpose() * velocity_L.transpose() * dNdx *
499 rho_LR * specific_heat_capacity_fluid) *
506 dNdx.transpose() * T_ip * K_pT_thermal_osmosis * dNdx * w;
507 dK_TT_dp.noalias() -= rho_LR * specific_heat_capacity_fluid *
508 N.transpose() * (dNdx * T).transpose() *
509 k_rel * Ki_over_mu * dNdx * w;
511 dK_TT_dp.noalias() -= rho_LR * specific_heat_capacity_fluid *
512 N.transpose() * velocity_L.dot(dNdx * T) /
513 k_rel * dk_rel_dS_L * dS_L_dp_cap * N * w;
515 if (gas_phase && S_L < 1.0)
519 double const rho_wv =
520 gas_phase->
property(MPL::PropertyType::density)
521 .template value<double>(variables, x_position, t, dt);
523 double const drho_wv_dT =
524 gas_phase->
property(MPL::PropertyType::density)
525 .template dValue<double>(variables,
526 MPL::Variable::temperature,
528 double const drho_wv_dp =
529 gas_phase->
property(MPL::PropertyType::density)
530 .template dValue<double>(
531 variables, MPL::Variable::liquid_phase_pressure,
536 MPL::PropertyType::thermal_diffusion_enhancement_factor)
537 .template value<double>(variables, x_position, t, dt);
540 auto const tortuosity =
541 medium.property(MPL::PropertyType::tortuosity)
542 .template value<double>(variables, x_position, t, dt);
544 phi * (1.0 - S_L) * tortuosity *
545 gas_phase->
property(MPL::PropertyType::diffusion)
546 .template value<double>(variables, x_position, t, dt);
548 double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
549 double const D_pv = D_v * drho_wv_dp;
553 -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
554 double const specific_heat_capacity_vapour =
556 .template value<double>(variables, x_position, t, dt);
559 w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
562 K_TT.noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
563 specific_heat_capacity_vapour * w;
565 double const storage_coefficient_by_water_vapor =
566 phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
568 storage_p_a_p.noalias() +=
569 N.transpose() * storage_coefficient_by_water_vapor * N * w;
571 double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
572 M_pT.noalias() += N.transpose() * vapor_expansion_factor * N * w;
575 .template block<pressure_size, temperature_size>(
576 pressure_index, temperature_index)
577 .noalias() += dNdx.transpose() * f_Tv_D_Tv * dNdx * w;
579 local_rhs.template segment<pressure_size>(pressure_index)
580 .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
582 laplace_p.noalias() += dNdx.transpose() * D_pv * dNdx * w;
587 if (gas_phase->
hasProperty(MPL::PropertyType::specific_latent_heat))
589 double const factor = phi * (1 - S_L) / rho_LR;
592 gas_phase->
property(MPL::PropertyType::specific_latent_heat)
593 .template value<double>(variables, x_position, t, dt) *
596 double const drho_LR_dT =
597 liquid_phase.property(MPL::PropertyType::density)
598 .template dValue<double>(variables,
599 MPL::Variable::temperature,
602 double const rho_wv_over_rho_L = rho_wv / rho_LR;
605 (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
606 N.transpose() * N * w;
610 (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
611 L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
612 N.transpose() * N * w;
616 L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
619 L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
624 if (_process_data.apply_mass_lumping)
626 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
627 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
629 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
637 .template block<temperature_size, temperature_size>(temperature_index,
639 .noalias() += M_TT / dt + K_TT;
642 .template block<temperature_size, pressure_size>(temperature_index,
644 .noalias() += K_Tp + dK_TT_dp;
648 .template block<pressure_size, pressure_size>(pressure_index,
650 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
654 .template block<pressure_size, temperature_size>(pressure_index,
656 .noalias() += M_pT / dt + laplace_T;
662 local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
663 M_TT * (T - T_prev) / dt + K_TT * T;
664 local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
668 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
669 laplace_p * p_L + laplace_T * T +
670 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
671 M_pT * (T - T_prev) / dt;
674 if (gas_phase->
hasProperty(MPL::PropertyType::specific_latent_heat))
678 .template block<temperature_size, pressure_size>(
679 temperature_index, pressure_index)
680 .noalias() += M_Tp / dt;
682 local_rhs.template segment<temperature_size>(temperature_index)
683 .noalias() -= M_Tp * (p_L - p_L_prev) / dt;
690 double const t,
double const dt, std::vector<double>
const& local_x,
691 std::vector<double>
const& local_x_prev, std::vector<double>& local_M_data,
692 std::vector<double>& local_K_data, std::vector<double>& local_rhs_data)
694 auto const local_matrix_dim = pressure_size + temperature_size;
695 assert(local_x.size() == local_matrix_dim);
697 auto const T = Eigen::Map<
typename ShapeMatricesType::template VectorType<
698 temperature_size>
const>(local_x.data() + temperature_index,
700 auto const p_L = Eigen::Map<
701 typename ShapeMatricesType::template VectorType<pressure_size>
const>(
702 local_x.data() + pressure_index, pressure_size);
704 auto const p_L_prev = Eigen::Map<
705 typename ShapeMatricesType::template VectorType<pressure_size>
const>(
706 local_x_prev.data() + pressure_index, pressure_size);
709 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
711 local_K_data, local_matrix_dim, local_matrix_dim);
714 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
716 local_M_data, local_matrix_dim, local_matrix_dim);
719 typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
720 local_rhs_data, local_matrix_dim);
722 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
723 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
724 auto const& solid_phase = medium.phase(
"Solid");
726 medium.hasPhase(
"Gas") ? &medium.phase(
"Gas") :
nullptr;
730 unsigned const n_integration_points =
731 _integration_method.getNumberOfPoints();
732 for (
unsigned ip = 0; ip < n_integration_points; ip++)
734 auto const& w = _ip_data[ip].integration_weight;
736 auto const& N = _ip_data[ip].N;
737 auto const& dNdx = _ip_data[ip].dNdx;
740 std::nullopt, _element.getID(), ip,
751 double p_cap_prev_ip;
761 auto& S_L = _ip_data[ip].saturation;
762 auto const S_L_prev = _ip_data[ip].saturation_prev;
764 medium[MPL::PropertyType::biot_coefficient].template value<double>(
765 variables, x_position, t, dt);
767 auto& solid_elasticity = *_process_data.simplified_elasticity;
771 solid_elasticity.bulkCompressibilityFromYoungsModulus(
772 solid_phase, variables, x_position, t, dt);
773 auto const beta_SR = (1 - alpha) * beta_S;
777 liquid_phase[MPL::PropertyType::density].template value<double>(
778 variables, x_position, t, dt);
779 auto const& b = _process_data.specific_body_force;
781 double const drho_LR_dp =
782 liquid_phase[MPL::PropertyType::density].template dValue<double>(
783 variables, MPL::Variable::liquid_phase_pressure, x_position, t,
785 auto const beta_LR = drho_LR_dp / rho_LR;
787 S_L = medium[MPL::PropertyType::saturation].template value<double>(
788 variables, x_position, t, dt);
793 double const dS_L_dp_cap =
794 medium[MPL::PropertyType::saturation].template dValue<double>(
795 variables, MPL::Variable::capillary_pressure, x_position, t,
799 double const DeltaS_L_Deltap_cap =
800 (p_cap_ip == p_cap_prev_ip)
802 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
805 auto chi_S_L_prev = S_L_prev;
806 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
808 auto const chi = [&medium, x_position, t, dt](
double const S_L)
812 return medium[MPL::PropertyType::bishops_effective_stress]
813 .template value<double>(variables, x_position, t, dt);
816 chi_S_L_prev = chi(S_L_prev);
826 auto& phi = _ip_data[ip].porosity;
829 variables_prev.
porosity = _ip_data[ip].porosity_prev;
830 phi = medium[MPL::PropertyType::porosity].template value<double>(
831 variables, variables_prev, x_position, t, dt);
838 "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
839 "porosity {} in element/integration point {}/{}.",
840 alpha, phi, _element.getID(), ip);
843 auto const K_pT_thermal_osmosis =
844 (solid_phase.hasProperty(
846 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
848 [MPL::PropertyType::thermal_osmosis_coefficient]
849 .value(variables, x_position, t, dt))
850 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
853 medium[MPL::PropertyType::relative_permeability]
854 .template value<double>(variables, x_position, t, dt);
856 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
857 variables, x_position, t, dt);
859 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
860 medium[MPL::PropertyType::permeability].value(variables, x_position,
869 Eigen::Matrix<double, 3, 3>
const
870 solid_linear_thermal_expansion_coefficient =
874 .value(variables, x_position, t, dt));
877 solid_phase[MPL::PropertyType::density].template value<double>(
878 variables, x_position, t, dt);
884 .template block<pressure_size, pressure_size>(pressure_index,
886 .noalias() += dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
888 const double alphaB_minus_phi = alpha - phi;
889 double const a0 = alphaB_minus_phi * beta_SR;
890 double const specific_storage_a_p =
891 S_L * (phi * beta_LR + S_L * a0 +
892 chi_S_L * alpha * alpha *
893 solid_elasticity.storageContribution(
894 solid_phase, variables, x_position, t, dt));
895 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
898 .template block<pressure_size, pressure_size>(pressure_index,
900 .noalias() += N.transpose() * rho_LR *
901 (specific_storage_a_p -
902 specific_storage_a_S * DeltaS_L_Deltap_cap) *
905 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
906 dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
911 double const fluid_volumetric_thermal_expansion_coefficient =
914 const double eff_thermal_expansion =
915 S_L * (alphaB_minus_phi *
916 solid_linear_thermal_expansion_coefficient.trace() +
917 phi * fluid_volumetric_thermal_expansion_coefficient +
918 alpha * solid_elasticity.thermalExpansivityContribution(
919 solid_linear_thermal_expansion_coefficient,
920 solid_phase, variables, x_position, t, dt));
923 .template block<pressure_size, temperature_size>(pressure_index,
926 dNdx.transpose() * rho_LR * K_pT_thermal_osmosis * dNdx * w;
929 .template block<pressure_size, temperature_size>(pressure_index,
932 N.transpose() * rho_LR * eff_thermal_expansion * N * w;
938 auto const specific_heat_capacity_fluid =
940 .template value<double>(variables, x_position, t, dt);
942 auto const specific_heat_capacity_solid =
945 .template value<double>(variables, x_position, t, dt);
948 .template block<temperature_size, temperature_size>(
949 temperature_index, temperature_index)
952 (rho_SR * specific_heat_capacity_solid * (1 - phi) +
953 (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
956 auto const thermal_conductivity =
957 MaterialPropertyLib::formEigenTensor<GlobalDim>(
959 thermal_conductivity]
960 .value(variables, x_position, t, dt));
963 -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b) -
964 K_pT_thermal_osmosis * dNdx * T);
967 .template block<temperature_size, temperature_size>(
968 temperature_index, temperature_index)
969 .noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
970 N.transpose() * velocity_L.transpose() * dNdx *
971 rho_LR * specific_heat_capacity_fluid) *
974 .template block<temperature_size, pressure_size>(
975 temperature_index, pressure_index)
977 dNdx.transpose() * T_ip * K_pT_thermal_osmosis * dNdx * w;
979 if (gas_phase && S_L < 1.0)
983 double const rho_wv =
985 .template value<double>(variables, x_position, t, dt);
987 double const drho_wv_dT =
989 .template dValue<double>(variables,
990 MPL::Variable::temperature,
992 double const drho_wv_dp =
994 .template dValue<double>(
995 variables, MPL::Variable::liquid_phase_pressure,
1000 MPL::PropertyType::thermal_diffusion_enhancement_factor)
1001 .template value<double>(variables, x_position, t, dt);
1004 auto const tortuosity =
1005 medium.property(MPL::PropertyType::tortuosity)
1006 .template value<double>(variables, x_position, t, dt);
1008 phi * (1.0 - S_L) * tortuosity *
1009 gas_phase->
property(MPL::PropertyType::diffusion)
1010 .template value<double>(variables, x_position, t, dt);
1012 double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
1013 double const D_pv = D_v * drho_wv_dp;
1018 -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
1019 double const specific_heat_capacity_vapour =
1021 .template value<double>(variables, x_position, t, dt);
1024 .template block<temperature_size, temperature_size>(
1025 temperature_index, temperature_index)
1027 w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
1031 .template block<temperature_size, temperature_size>(
1032 temperature_index, temperature_index)
1033 .noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
1034 specific_heat_capacity_vapour * w;
1036 double const storage_coefficient_by_water_vapor =
1037 phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
1039 .template block<pressure_size, pressure_size>(pressure_index,
1042 N.transpose() * storage_coefficient_by_water_vapor * N * w;
1044 double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
1046 .template block<pressure_size, temperature_size>(
1047 pressure_index, temperature_index)
1048 .noalias() += N.transpose() * vapor_expansion_factor * N * w;
1050 local_rhs.template segment<pressure_size>(pressure_index)
1051 .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
1054 .template block<pressure_size, pressure_size>(pressure_index,
1056 .noalias() += dNdx.transpose() * D_pv * dNdx * w;
1061 if (gas_phase->
hasProperty(MPL::PropertyType::specific_latent_heat))
1063 double const factor = phi * (1 - S_L) / rho_LR;
1066 gas_phase->
property(MPL::PropertyType::specific_latent_heat)
1067 .template value<double>(variables, x_position, t, dt) *
1070 double const drho_LR_dT =
1071 liquid_phase.property(MPL::PropertyType::density)
1072 .template dValue<double>(variables,
1073 MPL::Variable::temperature,
1076 double const rho_wv_over_rho_L = rho_wv / rho_LR;
1078 .template block<temperature_size, temperature_size>(
1079 temperature_index, temperature_index)
1082 (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
1083 N.transpose() * N * w;
1086 .template block<temperature_size, pressure_size>(
1087 temperature_index, pressure_index)
1090 (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
1091 L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
1092 N.transpose() * N * w;
1096 .template block<temperature_size, temperature_size>(
1097 temperature_index, temperature_index)
1099 L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
1102 .template block<temperature_size, pressure_size>(
1103 temperature_index, pressure_index)
1105 L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
1110 if (_process_data.apply_mass_lumping)
1112 auto Mpp = local_M.template block<pressure_size, pressure_size>(
1113 pressure_index, pressure_index);
1114 Mpp = Mpp.colwise().sum().eval().asDiagonal();
1201 Eigen::VectorXd
const& local_x,
1202 Eigen::VectorXd
const& local_x_prev)
1205 local_x.template segment<temperature_size>(temperature_index);
1207 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
1210 local_x_prev.template segment<pressure_size>(pressure_index);
1212 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
1213 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
1214 auto const& solid_phase = medium.phase(
"Solid");
1218 unsigned const n_integration_points =
1219 _integration_method.getNumberOfPoints();
1221 double saturation_avg = 0;
1222 double porosity_avg = 0;
1224 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1226 auto const& N = _ip_data[ip].N;
1229 std::nullopt, _element.getID(), ip,
1240 double p_cap_prev_ip;
1251 auto& S_L = _ip_data[ip].saturation;
1252 auto const S_L_prev = _ip_data[ip].saturation_prev;
1253 S_L = medium[MPL::PropertyType::saturation].template value<double>(
1254 variables, x_position, t, dt);
1259 auto chi_S_L_prev = S_L_prev;
1260 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
1262 auto const chi = [&medium, x_position, t, dt](
double const S_L)
1266 return medium[MPL::PropertyType::bishops_effective_stress]
1267 .template value<double>(variables, x_position, t, dt);
1270 chi_S_L_prev = chi(S_L_prev);
1276 medium[MPL::PropertyType::biot_coefficient].template value<double>(
1277 variables, x_position, t, dt);
1279 auto& solid_elasticity = *_process_data.simplified_elasticity;
1281 solid_elasticity.bulkCompressibilityFromYoungsModulus(
1282 solid_phase, variables, x_position, t, dt);
1283 auto const beta_SR = (1 - alpha) * beta_S;
1286 auto& phi = _ip_data[ip].porosity;
1288 variables_prev.
porosity = _ip_data[ip].porosity_prev;
1289 phi = medium[MPL::PropertyType::porosity].template value<double>(
1290 variables, variables_prev, x_position, t, dt);
1295 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
1296 variables, x_position, t, dt);
1298 liquid_phase[MPL::PropertyType::density].template value<double>(
1299 variables, x_position, t, dt);
1301 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
1302 medium[MPL::PropertyType::permeability].value(variables, x_position,
1305 double const k_rel =
1306 medium[MPL::PropertyType::relative_permeability]
1307 .template value<double>(variables, x_position, t, dt);
1312 solid_phase[MPL::PropertyType::density].template value<double>(
1313 variables, x_position, t, dt);
1314 _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
1316 auto const& b = _process_data.specific_body_force;
1318 auto const K_pT_thermal_osmosis =
1319 (solid_phase.hasProperty(
1321 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
1323 [MPL::PropertyType::thermal_osmosis_coefficient]
1324 .value(variables, x_position, t, dt))
1325 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
1328 auto const& dNdx = _ip_data[ip].dNdx;
1329 _ip_data[ip].v_darcy.noalias() = -K_over_mu * dNdx * p_L -
1330 K_pT_thermal_osmosis * dNdx * T +
1331 rho_LR * K_over_mu * b;
1333 saturation_avg += S_L;
1334 porosity_avg += phi;
1336 saturation_avg /= n_integration_points;
1337 porosity_avg /= n_integration_points;
1339 (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
1340 (*_process_data.element_porosity)[_element.getID()] = porosity_avg;