28 assemble(
double const t,
double const dt,
29 std::vector<double>
const& local_x,
30 std::vector<double>
const& ,
31 std::vector<double>& local_M_data,
32 std::vector<double>& local_K_data,
33 std::vector<double>& local_b_data)
38 auto const local_matrix_size = local_x.size();
40 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
42 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
43 local_M_data, local_matrix_size, local_matrix_size);
44 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
45 local_K_data, local_matrix_size, local_matrix_size);
46 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
47 local_b_data, local_matrix_size);
50 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
51 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
52 auto Mapc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
53 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
54 auto Max = local_M.template block<nonwet_pressure_size, contaminant_size>(
55 nonwet_pressure_matrix_index, contaminant_matrix_index);
56 auto Mat = local_M.template block<nonwet_pressure_size, temperature_size>(
57 nonwet_pressure_matrix_index, temperature_matrix_index);
59 auto Mwp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
60 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
61 auto Mwpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
62 cap_pressure_matrix_index, cap_pressure_matrix_index);
63 auto Mwx = local_M.template block<cap_pressure_size, contaminant_size>(
64 cap_pressure_matrix_index, contaminant_matrix_index);
65 auto Mwt = local_M.template block<cap_pressure_size, temperature_size>(
66 cap_pressure_matrix_index, temperature_matrix_index);
68 auto Mcp = local_M.template block<contaminant_size, nonwet_pressure_size>(
69 contaminant_matrix_index, nonwet_pressure_matrix_index);
70 auto Mcpc = local_M.template block<contaminant_size, cap_pressure_size>(
71 contaminant_matrix_index, cap_pressure_matrix_index);
72 auto Mcx = local_M.template block<contaminant_size, contaminant_size>(
73 contaminant_matrix_index, contaminant_matrix_index);
74 auto Mct = local_M.template block<contaminant_size, temperature_size>(
75 contaminant_matrix_index, temperature_matrix_index);
77 auto Mep = local_M.template block<temperature_size, nonwet_pressure_size>(
78 temperature_matrix_index, nonwet_pressure_matrix_index);
79 auto Mepc = local_M.template block<temperature_size, cap_pressure_size>(
80 temperature_matrix_index, cap_pressure_matrix_index);
81 auto Mex = local_M.template block<temperature_size, contaminant_size>(
82 temperature_matrix_index, contaminant_matrix_index);
83 auto Met = local_M.template block<temperature_size, temperature_size>(
84 temperature_matrix_index, temperature_matrix_index);
87 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
90 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
91 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
92 auto Kapc = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
93 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
94 auto Kax = local_K.template block<nonwet_pressure_size, contaminant_size>(
95 nonwet_pressure_matrix_index, contaminant_matrix_index);
96 auto Kat = local_K.template block<nonwet_pressure_size, temperature_size>(
97 nonwet_pressure_matrix_index, temperature_matrix_index);
99 auto Kwp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
100 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
101 auto Kwpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
102 cap_pressure_matrix_index, cap_pressure_matrix_index);
103 auto Kwx = local_K.template block<cap_pressure_size, contaminant_size>(
104 cap_pressure_matrix_index, contaminant_matrix_index);
105 auto Kwt = local_K.template block<cap_pressure_size, temperature_size>(
106 cap_pressure_matrix_index, temperature_matrix_index);
108 auto Kcp = local_K.template block<contaminant_size, nonwet_pressure_size>(
109 contaminant_matrix_index, nonwet_pressure_matrix_index);
110 auto Kcpc = local_K.template block<contaminant_size, cap_pressure_size>(
111 contaminant_matrix_index, cap_pressure_matrix_index);
112 auto Kcx = local_K.template block<contaminant_size, contaminant_size>(
113 contaminant_matrix_index, contaminant_matrix_index);
114 auto Kct = local_K.template block<contaminant_size, temperature_size>(
115 contaminant_matrix_index, temperature_matrix_index);
117 auto Kep = local_K.template block<temperature_size, nonwet_pressure_size>(
118 temperature_matrix_index, nonwet_pressure_matrix_index);
119 auto Kepc = local_K.template block<temperature_size, cap_pressure_size>(
120 temperature_matrix_index, cap_pressure_matrix_index);
121 auto Ket = local_K.template block<temperature_size, temperature_size>(
122 temperature_matrix_index, temperature_matrix_index);
124 auto Ba = local_b.template segment<nonwet_pressure_size>(
125 nonwet_pressure_matrix_index);
127 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
129 local_b.template segment<contaminant_size>(contaminant_matrix_index);
131 local_b.template segment<temperature_size>(temperature_matrix_index);
133 unsigned const n_integration_points =
134 _integration_method.getNumberOfPoints();
139 auto const num_nodes = ShapeFunction::NPOINTS;
140 auto const pg_nodal_values =
141 Eigen::Map<const NodalVectorType>(&local_x[0], num_nodes);
142 auto const pc_nodal_values =
143 Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
148 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
150 for (
unsigned ip = 0; ip < n_integration_points; ip++)
152 auto const& ip_data = _ip_data[ip];
153 auto const& N = ip_data.N;
154 auto const& dNdx = ip_data.dNdx;
155 auto const& w = ip_data.integration_weight;
156 auto const& mass_operator = ip_data.mass_operator;
157 auto const& diffusion_operator = ip_data.diffusion_operator;
158 double pg_int_pt = 0.;
159 double pc_int_pt = 0.;
160 double Xc_int_pt = 0.;
161 double T_int_pt = 0.;
163 Xc_int_pt, T_int_pt);
165 _pressure_wetting[ip] = pg_int_pt - pc_int_pt;
166 double const ideal_gas_constant_times_T_int_pt =
167 IdealGasConstant * T_int_pt;
173 *_process_data.media_map.getMedium(this->_element.getID());
174 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
175 auto const& solid_phase = medium.phase(
"Solid");
176 auto const& gas_phase = medium.phase(
"Gas");
178 auto const& water_vapour_component = gas_phase.component(
"w");
179 auto const& dry_air_component = gas_phase.component(
"a");
180 auto const& contaminant_vapour_component = gas_phase.component(
"c");
181 auto const& dissolved_contaminant_component =
182 liquid_phase.component(
"c");
184 auto const porosity =
186 .template value<double>(vars, pos, t, dt);
188 auto const water_mol_mass =
189 water_vapour_component
191 .template value<double>(vars, pos, t, dt);
192 auto const air_mol_mass =
195 .template value<double>(vars, pos, t, dt);
196 auto const contaminant_mol_mass =
197 contaminant_vapour_component
199 .template value<double>(vars, pos, t, dt);
203 .template value<double>(vars, pos, t, dt);
205 _saturation[ip] = Sw;
208 double const dSw_dpc =
210 .template dValue<double>(
214 auto const density_water =
216 .template value<double>(vars, pos, t, dt);
219 double const mol_density_water = density_water / water_mol_mass;
220 double const mol_density_wet = mol_density_water;
222 double const mol_density_nonwet =
223 pg_int_pt / ideal_gas_constant_times_T_int_pt;
224 double const mol_density_tot =
225 Sw * mol_density_wet + (1 - Sw) * mol_density_nonwet;
227 double const d_mol_density_nonwet_dpg =
228 1 / ideal_gas_constant_times_T_int_pt;
229 double const d_mol_density_nonwet_dT = -mol_density_nonwet / T_int_pt;
230 double const d_mol_density_tot_dpc =
231 (mol_density_wet - mol_density_nonwet) * dSw_dpc;
232 double const d_mol_density_tot_dpg =
233 (1 - Sw) * d_mol_density_nonwet_dpg;
234 double const d_mol_density_tot_dT = (1 - Sw) * d_mol_density_nonwet_dT;
237 double const latent_heat_evaporation =
238 water_vapour_component
241 .template value<double>(vars, pos, t, dt);
247 water_vapour_component
249 .template value<double>(vars, pos, t, dt);
250 double const dp_sat_dT =
251 water_vapour_component
253 .template dValue<double>(
258 double const K = std::exp(-pc_int_pt / mol_density_water /
259 ideal_gas_constant_times_T_int_pt);
260 double const dK_dT = pc_int_pt / mol_density_water /
261 ideal_gas_constant_times_T_int_pt / T_int_pt * K;
265 double const p_vapour_nonwet = p_sat * K;
266 double const d_p_vapour_nonwet_dT = dp_sat_dT * K + p_sat * dK_dT;
267 double const d_p_vapour_nonwet_dpc =
269 (-1 / mol_density_water / ideal_gas_constant_times_T_int_pt);
272 double const henry_contam =
273 contaminant_vapour_component
275 .template value<double>(vars, pos, t, dt);
276 double d_henry_contaminant_dT =
277 contaminant_vapour_component
279 .template dValue<double>(
284 double const k_c = pg_int_pt * henry_contam / mol_density_wet;
285 double const k_w = pg_int_pt / p_vapour_nonwet;
288 double const Ntot_c =
289 Sw * mol_density_wet * k_c + (1 - Sw) * mol_density_nonwet;
291 double const x_contaminant_nonwet =
292 Xc_int_pt * mol_density_tot / Ntot_c;
293 double const x_contaminant_wet = k_c * x_contaminant_nonwet;
294 double const x_water_wet = 1 - x_contaminant_wet;
295 double const x_water_nonwet = x_water_wet / k_w;
296 double const x_air_nonwet = 1 - x_water_nonwet - x_contaminant_nonwet;
298 _gas_molar_fraction_water[ip] = x_water_nonwet;
299 _liquid_molar_fraction_contaminant[ip] = x_contaminant_wet;
300 _gas_molar_fraction_contaminant[ip] = x_contaminant_nonwet;
302 double const d_kc_dpg = henry_contam / mol_density_wet;
303 double const d_kc_dT =
304 pg_int_pt * d_henry_contaminant_dT / mol_density_wet;
307 double const d_x_contaminant_nonwet_dpc =
308 Xc_int_pt * (d_mol_density_tot_dpc / Ntot_c -
309 (mol_density_wet * k_c - mol_density_nonwet) *
310 dSw_dpc * mol_density_tot / Ntot_c / Ntot_c);
311 double const d_x_contaminant_nonwet_dpg =
312 Xc_int_pt * (d_mol_density_tot_dpg / Ntot_c -
313 (Sw * mol_density_wet * d_kc_dpg +
314 (1 - Sw) * d_mol_density_nonwet_dpg) *
315 mol_density_tot / Ntot_c / Ntot_c);
316 double const d_x_contaminant_nonwet_dXc = mol_density_tot / Ntot_c;
317 double const d_x_contaminant_nonwet_dT =
318 Xc_int_pt * (d_mol_density_tot_dT / Ntot_c -
319 (Sw * mol_density_wet * d_kc_dT +
320 (1 - Sw) * d_mol_density_nonwet_dT) *
321 mol_density_tot / Ntot_c / Ntot_c);
323 double const d_x_contaminant_wet_dpc = k_c * d_x_contaminant_nonwet_dpc;
324 double const d_x_contaminant_wet_dpg =
325 k_c * d_x_contaminant_nonwet_dpg + d_kc_dpg * x_contaminant_nonwet;
326 double const d_x_contaminant_wet_dXc = k_c * d_x_contaminant_nonwet_dXc;
327 double const d_x_contaminant_wet_dT =
328 k_c * d_x_contaminant_nonwet_dT + d_kc_dT * x_contaminant_nonwet;
330 double const d_x_water_wet_dpc = -d_x_contaminant_wet_dpc;
331 double const d_x_water_wet_dpg = -d_x_contaminant_wet_dpg;
332 double const d_x_water_wet_dXc = -d_x_contaminant_wet_dXc;
333 double const d_x_water_wet_dT = -d_x_contaminant_wet_dT;
335 double const d_x_water_nonwet_dpc =
336 (d_p_vapour_nonwet_dpc * x_water_wet +
337 p_vapour_nonwet * d_x_water_wet_dpc) /
339 double const d_x_water_nonwet_dpg =
340 p_vapour_nonwet * (d_x_water_wet_dpg / pg_int_pt -
341 x_water_wet / pg_int_pt / pg_int_pt);
342 double const d_x_water_nonwet_dXc = d_x_water_wet_dXc / k_w;
343 double const d_x_water_nonwet_dT =
344 (d_p_vapour_nonwet_dT * x_water_wet +
345 p_vapour_nonwet * d_x_water_wet_dT) /
348 double const d_x_air_nonwet_dpc =
349 -d_x_water_nonwet_dpc - d_x_contaminant_nonwet_dpc;
350 double const d_x_air_nonwet_dpg =
351 -d_x_water_nonwet_dpg - d_x_contaminant_nonwet_dpg;
352 double const d_x_air_nonwet_dXc =
353 -d_x_water_nonwet_dXc - d_x_contaminant_nonwet_dXc;
354 double const d_x_air_nonwet_dT =
355 -d_x_water_nonwet_dT - d_x_contaminant_nonwet_dT;
358 double const density_contaminant_wet =
359 mol_density_wet * contaminant_mol_mass * x_contaminant_wet;
360 double const density_wet = density_water + density_contaminant_wet;
362 double const density_air_nonwet =
363 mol_density_nonwet * air_mol_mass * x_air_nonwet;
364 double const density_water_nonwet =
365 mol_density_nonwet * water_mol_mass * x_water_nonwet;
366 double const density_contaminant_nonwet =
367 mol_density_nonwet * contaminant_mol_mass * x_contaminant_nonwet;
368 double const density_nonwet = density_air_nonwet +
369 density_water_nonwet +
370 density_contaminant_nonwet;
372 auto const density_solid =
374 .template value<double>(vars, pos, t, dt);
377 double const d_density_nonwet_dpg =
378 d_mol_density_nonwet_dpg *
379 (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
380 contaminant_mol_mass * x_contaminant_nonwet) +
382 (air_mol_mass * d_x_air_nonwet_dpg +
383 water_mol_mass * d_x_water_nonwet_dpg +
384 contaminant_mol_mass * d_x_contaminant_nonwet_dpg);
385 double const d_density_nonwet_dpc =
387 (air_mol_mass * d_x_air_nonwet_dpc +
388 water_mol_mass * d_x_water_nonwet_dpc +
389 contaminant_mol_mass * d_x_contaminant_nonwet_dpc);
390 double const d_density_nonwet_dXc =
392 (air_mol_mass * d_x_air_nonwet_dXc +
393 water_mol_mass * d_x_water_nonwet_dXc +
394 contaminant_mol_mass * d_x_contaminant_nonwet_dXc);
395 double const d_density_nonwet_dT =
396 d_mol_density_nonwet_dT *
397 (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
398 contaminant_mol_mass * x_contaminant_nonwet) +
400 (air_mol_mass * d_x_air_nonwet_dT +
401 water_mol_mass * d_x_water_nonwet_dT +
402 contaminant_mol_mass * d_x_contaminant_nonwet_dT);
404 double const mol_mass_nonwet =
405 x_water_nonwet * water_mol_mass + x_air_nonwet * air_mol_mass +
406 x_contaminant_nonwet * contaminant_mol_mass;
408 double const X_water_nonwet =
409 x_water_nonwet * water_mol_mass / mol_mass_nonwet;
410 double const X_air_nonwet =
411 x_air_nonwet * air_mol_mass / mol_mass_nonwet;
412 double const X_contaminant_nonwet = 1 - X_water_nonwet - X_air_nonwet;
415 double const heat_capacity_dry_air =
419 .template value<double>(vars, pos, t, dt);
420 double const heat_capacity_water_vapour =
421 water_vapour_component
424 .template value<double>(vars, pos, t, dt);
425 double const heat_capacity_contaminant_vapour =
426 contaminant_vapour_component
429 .template value<double>(vars, pos, t, dt);
431 double const heat_capacity_water =
435 .template value<double>(vars, pos, t, dt);
436 double const heat_capacity_solid =
440 .template value<double>(vars, pos, t, dt);
445 double const enthalpy_air_nonwet =
446 heat_capacity_dry_air * (T_int_pt - CelsiusZeroInKelvin) +
447 IdealGasConstant * T_int_pt / air_mol_mass;
448 double const enthalpy_water_nonwet =
449 heat_capacity_water_vapour * (T_int_pt - CelsiusZeroInKelvin) +
450 latent_heat_evaporation;
451 double const enthalpy_contaminant_nonwet =
452 heat_capacity_contaminant_vapour *
453 (T_int_pt - CelsiusZeroInKelvin) +
454 IdealGasConstant * T_int_pt / contaminant_mol_mass;
457 double const enthalpy_nonwet =
458 enthalpy_air_nonwet * X_air_nonwet +
459 enthalpy_water_nonwet * X_water_nonwet +
460 enthalpy_contaminant_nonwet * X_contaminant_nonwet;
461 double const enthalpy_wet =
462 heat_capacity_water * (T_int_pt - CelsiusZeroInKelvin);
465 double const internal_energy_nonwet =
466 enthalpy_nonwet - pg_int_pt / density_nonwet;
467 double const internal_energy_wet = enthalpy_wet;
470 double const d_enthalpy_air_nonwet_dT =
471 heat_capacity_dry_air + IdealGasConstant / air_mol_mass;
472 double const d_enthalpy_contaminant_nonwet_dT =
473 heat_capacity_contaminant_vapour +
474 IdealGasConstant / contaminant_mol_mass;
476 double const d_enthalpy_nonwet_dT =
477 heat_capacity_water * X_water_nonwet +
478 d_enthalpy_air_nonwet_dT * X_air_nonwet +
479 d_enthalpy_contaminant_nonwet_dT * X_contaminant_nonwet;
482 Map.noalias() += porosity * (1 - Sw) *
483 (mol_density_nonwet * d_x_air_nonwet_dpg +
484 x_air_nonwet * d_mol_density_nonwet_dpg) *
487 porosity * mol_density_nonwet *
488 ((1 - Sw) * d_x_air_nonwet_dpc - x_air_nonwet * dSw_dpc) *
490 Max.noalias() += porosity * mol_density_nonwet * (1 - Sw) *
491 d_x_air_nonwet_dXc * mass_operator;
492 Mat.noalias() += porosity * (1 - Sw) *
493 (mol_density_nonwet * d_x_air_nonwet_dT +
494 x_air_nonwet * d_mol_density_nonwet_dT) *
499 (mol_density_wet * Sw * d_x_water_wet_dpg +
500 (1 - Sw) * x_water_nonwet * d_mol_density_nonwet_dpg +
501 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dpg) *
506 (x_water_wet * dSw_dpc + Sw * d_x_water_wet_dpc) +
508 ((1 - Sw) * d_x_water_nonwet_dpc - x_water_nonwet * dSw_dpc)) *
512 (mol_density_wet * Sw * d_x_water_wet_dXc +
513 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dXc) *
515 Mwt.noalias() += porosity *
516 ((1 - Sw) / ideal_gas_constant_times_T_int_pt *
517 (d_p_vapour_nonwet_dT - p_vapour_nonwet / T_int_pt)) *
522 (mol_density_wet * Sw * d_x_contaminant_wet_dpg +
523 (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dpg +
524 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dpg) *
529 (x_contaminant_wet * dSw_dpc + Sw * d_x_contaminant_wet_dpc) +
530 mol_density_nonwet * ((1 - Sw) * d_x_contaminant_nonwet_dpc -
531 x_contaminant_nonwet * dSw_dpc)) *
535 (mol_density_wet * Sw * d_x_contaminant_wet_dXc +
536 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dXc) *
540 (mol_density_wet * Sw * d_x_contaminant_wet_dT +
541 (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dT +
542 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dT) *
545 Mep.noalias() += porosity *
546 (d_density_nonwet_dpg * enthalpy_nonwet - 1) *
547 (1 - Sw) * mass_operator;
548 Mepc.noalias() += porosity *
549 (density_wet * internal_energy_wet -
550 density_nonwet * internal_energy_nonwet) *
551 dSw_dpc * mass_operator +
552 porosity * d_density_nonwet_dpc * enthalpy_nonwet *
553 (1 - Sw) * mass_operator;
554 Mex.noalias() += porosity * d_density_nonwet_dXc * enthalpy_nonwet *
555 (1 - Sw) * mass_operator;
557 ((1 - porosity) * density_solid * heat_capacity_solid +
558 porosity * ((1 - Sw) * (d_density_nonwet_dT * enthalpy_nonwet +
559 density_nonwet * d_enthalpy_nonwet_dT) +
560 Sw * density_wet * heat_capacity_water)) *
564 double const diffusion_coeff_water_nonwet =
565 water_vapour_component
567 .template value<double>(vars, pos, t, dt);
568 double const diffusion_coeff_contaminant_nonwet =
569 contaminant_vapour_component
571 .template value<double>(vars, pos, t, dt);
572 double const diffusion_coeff_contaminant_wet =
573 dissolved_contaminant_component
575 .template value<double>(vars, pos, t, dt);
578 double const k_rel_nonwet =
581 relative_permeability_nonwetting_phase)
582 .
template value<double>(vars, pos, t, dt);
583 auto const mu_nonwet =
585 .template value<double>(vars, pos, t, dt);
586 double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
589 double const k_rel_wet =
593 .template value<double>(vars, pos, t, dt);
596 .template value<double>(vars, pos, t, dt);
597 double const lambda_wet = k_rel_wet / mu_wet;
600 auto const permeability =
601 MaterialPropertyLib::formEigenTensor<GlobalDim>(
603 .value(vars, pos, t, dt));
606 auto const& b = _process_data.specific_body_force;
610 _process_data.has_gravity
612 -lambda_wet * permeability *
613 (dNdx * (pg_nodal_values - pc_nodal_values) -
616 (pg_nodal_values - pc_nodal_values));
618 laplace_operator.noalias() = dNdx.transpose() * permeability * dNdx * w;
621 auto const solute_dispersivity_transverse =
622 medium.template value<double>(
624 auto const solute_dispersivity_longitudinal =
625 medium.template value<double>(
628 double const velocity_wet_magnitude = velocity_wet.norm();
630 velocity_wet_magnitude != 0.0
632 (porosity * Sw * diffusion_coeff_contaminant_wet +
633 solute_dispersivity_transverse *
634 velocity_wet_magnitude) *
636 (solute_dispersivity_longitudinal -
637 solute_dispersivity_transverse) /
638 velocity_wet_magnitude * velocity_wet *
639 velocity_wet.transpose())
641 (porosity * Sw * diffusion_coeff_contaminant_wet +
642 solute_dispersivity_transverse *
643 velocity_wet_magnitude) *
646 auto const dispersion_operator =
647 dNdx.transpose() * hydrodynamic_dispersion * dNdx * w;
652 (mol_density_nonwet * x_air_nonwet * lambda_nonwet) *
654 (1 - Sw) * porosity * mol_density_nonwet *
655 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpg +
656 diffusion_coeff_contaminant_nonwet *
657 d_x_contaminant_nonwet_dpg) *
660 -(1 - Sw) * porosity * mol_density_nonwet *
661 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpc +
662 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dpc) *
665 -(1 - Sw) * porosity * mol_density_nonwet *
666 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dXc +
667 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dXc) *
670 -(1 - Sw) * porosity * mol_density_nonwet *
671 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dT +
672 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT) *
676 (mol_density_wet * x_water_wet * lambda_wet +
677 mol_density_nonwet * x_water_nonwet * lambda_nonwet) *
680 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
682 (1 - Sw) * diffusion_coeff_water_nonwet * mol_density_nonwet *
683 d_x_water_nonwet_dpg) *
686 -mol_density_wet * x_water_wet * lambda_wet * laplace_operator +
688 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
690 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
691 d_x_water_nonwet_dpc) *
695 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
697 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
698 d_x_water_nonwet_dXc) *
702 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
704 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
705 d_x_water_nonwet_dT) *
709 (mol_density_wet * x_contaminant_wet * lambda_wet +
710 mol_density_nonwet * x_contaminant_nonwet * lambda_nonwet) *
712 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpg +
713 porosity * (1 - Sw) * mol_density_nonwet *
714 diffusion_coeff_contaminant_nonwet *
715 d_x_contaminant_nonwet_dpg * diffusion_operator;
717 -mol_density_wet * x_contaminant_wet * lambda_wet *
719 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpc +
720 porosity * (1 - Sw) * mol_density_nonwet *
721 diffusion_coeff_contaminant_nonwet *
722 d_x_contaminant_nonwet_dpc * diffusion_operator;
724 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dXc +
725 porosity * (1 - Sw) * mol_density_nonwet *
726 diffusion_coeff_contaminant_nonwet *
727 d_x_contaminant_nonwet_dXc * diffusion_operator;
729 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dT +
730 porosity * (1 - Sw) * mol_density_nonwet *
731 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT *
734 Kep.noalias() += (lambda_nonwet * density_nonwet * enthalpy_nonwet +
735 lambda_wet * density_wet * enthalpy_wet) *
738 -lambda_wet * enthalpy_wet * density_wet * laplace_operator;
740 if (medium.hasProperty(
747 .value(vars, pos, t, dt);
750 MaterialPropertyLib::formEigenTensor<GlobalDim>(lambda);
753 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
757 auto const thermal_conductivity_solid =
761 .value(vars, pos, t, dt);
763 auto const thermal_conductivity_fluid =
767 .template value<double>(vars, pos, t, dt) *
772 GlobalDim>(thermal_conductivity_solid,
773 thermal_conductivity_fluid, porosity);
776 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
779 if (_process_data.has_gravity)
782 dNdx.transpose() * permeability * b * w;
783 Ba.noalias() += (mol_density_nonwet * x_air_nonwet * lambda_nonwet *
787 (mol_density_wet * x_water_wet * lambda_wet * density_wet +
788 mol_density_nonwet * x_water_nonwet * lambda_nonwet *
791 Bc.noalias() += (mol_density_wet * x_contaminant_wet * lambda_wet *
793 mol_density_nonwet * x_contaminant_nonwet *
794 lambda_nonwet * density_nonwet) *
797 (lambda_nonwet * density_nonwet * density_nonwet *
799 lambda_wet * density_wet * density_wet * enthalpy_wet) *
803 if (_process_data.has_mass_lumping)
805 Map = Map.colwise().sum().eval().asDiagonal();
806 Mapc = Mapc.colwise().sum().eval().asDiagonal();
807 Max = Max.colwise().sum().eval().asDiagonal();
808 Mat = Mat.colwise().sum().eval().asDiagonal();
809 Mwp = Mwp.colwise().sum().eval().asDiagonal();
810 Mwpc = Mwpc.colwise().sum().eval().asDiagonal();
811 Mwx = Mwx.colwise().sum().eval().asDiagonal();
812 Mwt = Mwt.colwise().sum().eval().asDiagonal();
813 Mcp = Mcp.colwise().sum().eval().asDiagonal();
814 Mcpc = Mcpc.colwise().sum().eval().asDiagonal();
815 Mcx = Mcx.colwise().sum().eval().asDiagonal();
816 Mct = Mct.colwise().sum().eval().asDiagonal();
817 Mep = Mep.colwise().sum().eval().asDiagonal();
818 Mepc = Mepc.colwise().sum().eval().asDiagonal();
819 Mex = Mex.colwise().sum().eval().asDiagonal();
820 Met = Met.colwise().sum().eval().asDiagonal();