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);
43 local_M_data, local_matrix_size, local_matrix_size);
45 local_K_data, local_matrix_size, local_matrix_size);
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();
136 auto const num_nodes = ShapeFunction::NPOINTS;
137 auto const pg_nodal_values =
138 Eigen::Map<const NodalVectorType>(&local_x[0], num_nodes);
139 auto const pc_nodal_values =
140 Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
145 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
147 for (
unsigned ip = 0; ip < n_integration_points; ip++)
149 auto const& ip_data = _ip_data[ip];
150 auto const& N = ip_data.N;
151 auto const& dNdx = ip_data.dNdx;
152 auto const& w = ip_data.integration_weight;
153 auto const& mass_operator = ip_data.mass_operator;
154 auto const& diffusion_operator = ip_data.diffusion_operator;
155 double pg_int_pt = 0.;
156 double pc_int_pt = 0.;
157 double Xc_int_pt = 0.;
158 double T_int_pt = 0.;
160 Xc_int_pt, T_int_pt);
163 std::nullopt, _element.getID(),
168 _pressure_wetting[ip] = pg_int_pt - pc_int_pt;
169 double const ideal_gas_constant_times_T_int_pt =
170 IdealGasConstant * T_int_pt;
176 *_process_data.media_map.getMedium(this->_element.getID());
177 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
178 auto const& solid_phase = medium.phase(
"Solid");
179 auto const& gas_phase = medium.phase(
"Gas");
181 auto const& water_vapour_component = gas_phase.component(
"w");
182 auto const& dry_air_component = gas_phase.component(
"a");
183 auto const& contaminant_vapour_component = gas_phase.component(
"c");
184 auto const& dissolved_contaminant_component =
185 liquid_phase.component(
"c");
187 auto const porosity =
189 .template value<double>(vars, pos, t, dt);
191 auto const water_mol_mass =
192 water_vapour_component
194 .template value<double>(vars, pos, t, dt);
195 auto const air_mol_mass =
198 .template value<double>(vars, pos, t, dt);
199 auto const contaminant_mol_mass =
200 contaminant_vapour_component
202 .template value<double>(vars, pos, t, dt);
206 .template value<double>(vars, pos, t, dt);
208 _saturation[ip] = Sw;
211 double const dSw_dpc =
213 .template dValue<double>(
217 auto const density_water =
219 .template value<double>(vars, pos, t, dt);
222 double const mol_density_water = density_water / water_mol_mass;
223 double const mol_density_wet = mol_density_water;
225 double const mol_density_nonwet =
226 pg_int_pt / ideal_gas_constant_times_T_int_pt;
227 double const mol_density_tot =
228 Sw * mol_density_wet + (1 - Sw) * mol_density_nonwet;
230 double const d_mol_density_nonwet_dpg =
231 1 / ideal_gas_constant_times_T_int_pt;
232 double const d_mol_density_nonwet_dT = -mol_density_nonwet / T_int_pt;
233 double const d_mol_density_tot_dpc =
234 (mol_density_wet - mol_density_nonwet) * dSw_dpc;
235 double const d_mol_density_tot_dpg =
236 (1 - Sw) * d_mol_density_nonwet_dpg;
237 double const d_mol_density_tot_dT = (1 - Sw) * d_mol_density_nonwet_dT;
240 double const latent_heat_evaporation =
241 water_vapour_component
244 .template value<double>(vars, pos, t, dt);
251 water_vapour_component
253 .template value<double>(vars, pos, t, dt);
254 double const dp_sat_dT =
255 water_vapour_component
257 .template dValue<double>(
262 double const K = std::exp(-pc_int_pt / mol_density_water /
263 ideal_gas_constant_times_T_int_pt);
264 double const dK_dT = pc_int_pt / mol_density_water /
265 ideal_gas_constant_times_T_int_pt / T_int_pt * K;
269 double const p_vapour_nonwet = p_sat * K;
270 double const d_p_vapour_nonwet_dT = dp_sat_dT * K + p_sat * dK_dT;
271 double const d_p_vapour_nonwet_dpc =
273 (-1 / mol_density_water / ideal_gas_constant_times_T_int_pt);
276 double const henry_contam =
277 contaminant_vapour_component
279 .template value<double>(vars, pos, t, dt);
280 double d_henry_contaminant_dT =
281 contaminant_vapour_component
283 .template dValue<double>(
288 double const k_c = pg_int_pt * henry_contam / mol_density_wet;
289 double const k_w = pg_int_pt / p_vapour_nonwet;
292 double const Ntot_c =
293 Sw * mol_density_wet * k_c + (1 - Sw) * mol_density_nonwet;
295 double const x_contaminant_nonwet =
296 Xc_int_pt * mol_density_tot / Ntot_c;
297 double const x_contaminant_wet = k_c * x_contaminant_nonwet;
298 double const x_water_wet = 1 - x_contaminant_wet;
299 double const x_water_nonwet = x_water_wet / k_w;
300 double const x_air_nonwet = 1 - x_water_nonwet - x_contaminant_nonwet;
302 _gas_molar_fraction_water[ip] = x_water_nonwet;
303 _liquid_molar_fraction_contaminant[ip] = x_contaminant_wet;
304 _gas_molar_fraction_contaminant[ip] = x_contaminant_nonwet;
306 double const d_kc_dpg = henry_contam / mol_density_wet;
307 double const d_kc_dT =
308 pg_int_pt * d_henry_contaminant_dT / mol_density_wet;
311 double const d_x_contaminant_nonwet_dpc =
312 Xc_int_pt * (d_mol_density_tot_dpc / Ntot_c -
313 (mol_density_wet * k_c - mol_density_nonwet) *
314 dSw_dpc * mol_density_tot / Ntot_c / Ntot_c);
315 double const d_x_contaminant_nonwet_dpg =
316 Xc_int_pt * (d_mol_density_tot_dpg / Ntot_c -
317 (Sw * mol_density_wet * d_kc_dpg +
318 (1 - Sw) * d_mol_density_nonwet_dpg) *
319 mol_density_tot / Ntot_c / Ntot_c);
320 double const d_x_contaminant_nonwet_dXc = mol_density_tot / Ntot_c;
321 double const d_x_contaminant_nonwet_dT =
322 Xc_int_pt * (d_mol_density_tot_dT / Ntot_c -
323 (Sw * mol_density_wet * d_kc_dT +
324 (1 - Sw) * d_mol_density_nonwet_dT) *
325 mol_density_tot / Ntot_c / Ntot_c);
327 double const d_x_contaminant_wet_dpc = k_c * d_x_contaminant_nonwet_dpc;
328 double const d_x_contaminant_wet_dpg =
329 k_c * d_x_contaminant_nonwet_dpg + d_kc_dpg * x_contaminant_nonwet;
330 double const d_x_contaminant_wet_dXc = k_c * d_x_contaminant_nonwet_dXc;
331 double const d_x_contaminant_wet_dT =
332 k_c * d_x_contaminant_nonwet_dT + d_kc_dT * x_contaminant_nonwet;
334 double const d_x_water_wet_dpc = -d_x_contaminant_wet_dpc;
335 double const d_x_water_wet_dpg = -d_x_contaminant_wet_dpg;
336 double const d_x_water_wet_dXc = -d_x_contaminant_wet_dXc;
337 double const d_x_water_wet_dT = -d_x_contaminant_wet_dT;
339 double const d_x_water_nonwet_dpc =
340 (d_p_vapour_nonwet_dpc * x_water_wet +
341 p_vapour_nonwet * d_x_water_wet_dpc) /
343 double const d_x_water_nonwet_dpg =
344 p_vapour_nonwet * (d_x_water_wet_dpg / pg_int_pt -
345 x_water_wet / pg_int_pt / pg_int_pt);
346 double const d_x_water_nonwet_dXc = d_x_water_wet_dXc / k_w;
347 double const d_x_water_nonwet_dT =
348 (d_p_vapour_nonwet_dT * x_water_wet +
349 p_vapour_nonwet * d_x_water_wet_dT) /
352 double const d_x_air_nonwet_dpc =
353 -d_x_water_nonwet_dpc - d_x_contaminant_nonwet_dpc;
354 double const d_x_air_nonwet_dpg =
355 -d_x_water_nonwet_dpg - d_x_contaminant_nonwet_dpg;
356 double const d_x_air_nonwet_dXc =
357 -d_x_water_nonwet_dXc - d_x_contaminant_nonwet_dXc;
358 double const d_x_air_nonwet_dT =
359 -d_x_water_nonwet_dT - d_x_contaminant_nonwet_dT;
362 double const density_contaminant_wet =
363 mol_density_wet * contaminant_mol_mass * x_contaminant_wet;
364 double const density_wet = density_water + density_contaminant_wet;
366 double const density_air_nonwet =
367 mol_density_nonwet * air_mol_mass * x_air_nonwet;
368 double const density_water_nonwet =
369 mol_density_nonwet * water_mol_mass * x_water_nonwet;
370 double const density_contaminant_nonwet =
371 mol_density_nonwet * contaminant_mol_mass * x_contaminant_nonwet;
372 double const density_nonwet = density_air_nonwet +
373 density_water_nonwet +
374 density_contaminant_nonwet;
376 auto const density_solid =
378 .template value<double>(vars, pos, t, dt);
381 double const d_density_nonwet_dpg =
382 d_mol_density_nonwet_dpg *
383 (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
384 contaminant_mol_mass * x_contaminant_nonwet) +
386 (air_mol_mass * d_x_air_nonwet_dpg +
387 water_mol_mass * d_x_water_nonwet_dpg +
388 contaminant_mol_mass * d_x_contaminant_nonwet_dpg);
389 double const d_density_nonwet_dpc =
391 (air_mol_mass * d_x_air_nonwet_dpc +
392 water_mol_mass * d_x_water_nonwet_dpc +
393 contaminant_mol_mass * d_x_contaminant_nonwet_dpc);
394 double const d_density_nonwet_dXc =
396 (air_mol_mass * d_x_air_nonwet_dXc +
397 water_mol_mass * d_x_water_nonwet_dXc +
398 contaminant_mol_mass * d_x_contaminant_nonwet_dXc);
399 double const d_density_nonwet_dT =
400 d_mol_density_nonwet_dT *
401 (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
402 contaminant_mol_mass * x_contaminant_nonwet) +
404 (air_mol_mass * d_x_air_nonwet_dT +
405 water_mol_mass * d_x_water_nonwet_dT +
406 contaminant_mol_mass * d_x_contaminant_nonwet_dT);
408 double const mol_mass_nonwet =
409 x_water_nonwet * water_mol_mass + x_air_nonwet * air_mol_mass +
410 x_contaminant_nonwet * contaminant_mol_mass;
412 double const X_water_nonwet =
413 x_water_nonwet * water_mol_mass / mol_mass_nonwet;
414 double const X_air_nonwet =
415 x_air_nonwet * air_mol_mass / mol_mass_nonwet;
416 double const X_contaminant_nonwet = 1 - X_water_nonwet - X_air_nonwet;
419 double const heat_capacity_dry_air =
423 .template value<double>(vars, pos, t, dt);
424 double const heat_capacity_water_vapour =
425 water_vapour_component
428 .template value<double>(vars, pos, t, dt);
429 double const heat_capacity_contaminant_vapour =
430 contaminant_vapour_component
433 .template value<double>(vars, pos, t, dt);
435 double const heat_capacity_water =
439 .template value<double>(vars, pos, t, dt);
440 double const heat_capacity_solid =
444 .template value<double>(vars, pos, t, dt);
449 double const enthalpy_air_nonwet =
450 heat_capacity_dry_air * (T_int_pt - CelsiusZeroInKelvin) +
451 IdealGasConstant * T_int_pt / air_mol_mass;
452 double const enthalpy_water_nonwet =
453 heat_capacity_water_vapour * (T_int_pt - CelsiusZeroInKelvin) +
454 latent_heat_evaporation;
455 double const enthalpy_contaminant_nonwet =
456 heat_capacity_contaminant_vapour *
457 (T_int_pt - CelsiusZeroInKelvin) +
458 IdealGasConstant * T_int_pt / contaminant_mol_mass;
461 double const enthalpy_nonwet =
462 enthalpy_air_nonwet * X_air_nonwet +
463 enthalpy_water_nonwet * X_water_nonwet +
464 enthalpy_contaminant_nonwet * X_contaminant_nonwet;
465 double const enthalpy_wet =
466 heat_capacity_water * (T_int_pt - CelsiusZeroInKelvin);
469 double const internal_energy_nonwet =
470 enthalpy_nonwet - pg_int_pt / density_nonwet;
471 double const internal_energy_wet = enthalpy_wet;
474 double const d_enthalpy_air_nonwet_dT =
475 heat_capacity_dry_air + IdealGasConstant / air_mol_mass;
476 double const d_enthalpy_contaminant_nonwet_dT =
477 heat_capacity_contaminant_vapour +
478 IdealGasConstant / contaminant_mol_mass;
480 double const d_enthalpy_nonwet_dT =
481 heat_capacity_water * X_water_nonwet +
482 d_enthalpy_air_nonwet_dT * X_air_nonwet +
483 d_enthalpy_contaminant_nonwet_dT * X_contaminant_nonwet;
486 Map.noalias() += porosity * (1 - Sw) *
487 (mol_density_nonwet * d_x_air_nonwet_dpg +
488 x_air_nonwet * d_mol_density_nonwet_dpg) *
491 porosity * mol_density_nonwet *
492 ((1 - Sw) * d_x_air_nonwet_dpc - x_air_nonwet * dSw_dpc) *
494 Max.noalias() += porosity * mol_density_nonwet * (1 - Sw) *
495 d_x_air_nonwet_dXc * mass_operator;
496 Mat.noalias() += porosity * (1 - Sw) *
497 (mol_density_nonwet * d_x_air_nonwet_dT +
498 x_air_nonwet * d_mol_density_nonwet_dT) *
503 (mol_density_wet * Sw * d_x_water_wet_dpg +
504 (1 - Sw) * x_water_nonwet * d_mol_density_nonwet_dpg +
505 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dpg) *
510 (x_water_wet * dSw_dpc + Sw * d_x_water_wet_dpc) +
512 ((1 - Sw) * d_x_water_nonwet_dpc - x_water_nonwet * dSw_dpc)) *
516 (mol_density_wet * Sw * d_x_water_wet_dXc +
517 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dXc) *
519 Mwt.noalias() += porosity *
520 ((1 - Sw) / ideal_gas_constant_times_T_int_pt *
521 (d_p_vapour_nonwet_dT - p_vapour_nonwet / T_int_pt)) *
526 (mol_density_wet * Sw * d_x_contaminant_wet_dpg +
527 (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dpg +
528 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dpg) *
533 (x_contaminant_wet * dSw_dpc + Sw * d_x_contaminant_wet_dpc) +
534 mol_density_nonwet * ((1 - Sw) * d_x_contaminant_nonwet_dpc -
535 x_contaminant_nonwet * dSw_dpc)) *
539 (mol_density_wet * Sw * d_x_contaminant_wet_dXc +
540 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dXc) *
544 (mol_density_wet * Sw * d_x_contaminant_wet_dT +
545 (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dT +
546 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dT) *
549 Mep.noalias() += porosity *
550 (d_density_nonwet_dpg * enthalpy_nonwet - 1) *
551 (1 - Sw) * mass_operator;
552 Mepc.noalias() += porosity *
553 (density_wet * internal_energy_wet -
554 density_nonwet * internal_energy_nonwet) *
555 dSw_dpc * mass_operator +
556 porosity * d_density_nonwet_dpc * enthalpy_nonwet *
557 (1 - Sw) * mass_operator;
558 Mex.noalias() += porosity * d_density_nonwet_dXc * enthalpy_nonwet *
559 (1 - Sw) * mass_operator;
561 ((1 - porosity) * density_solid * heat_capacity_solid +
562 porosity * ((1 - Sw) * (d_density_nonwet_dT * enthalpy_nonwet +
563 density_nonwet * d_enthalpy_nonwet_dT) +
564 Sw * density_wet * heat_capacity_water)) *
568 double const diffusion_coeff_water_nonwet =
569 water_vapour_component
571 .template value<double>(vars, pos, t, dt);
572 double const diffusion_coeff_contaminant_nonwet =
573 contaminant_vapour_component
575 .template value<double>(vars, pos, t, dt);
576 double const diffusion_coeff_contaminant_wet =
577 dissolved_contaminant_component
579 .template value<double>(vars, pos, t, dt);
582 double const k_rel_nonwet =
585 relative_permeability_nonwetting_phase)
586 .template value<double>(vars, pos, t, dt);
587 auto const mu_nonwet =
589 .template value<double>(vars, pos, t, dt);
590 double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
593 double const k_rel_wet =
597 .template value<double>(vars, pos, t, dt);
600 .template value<double>(vars, pos, t, dt);
601 double const lambda_wet = k_rel_wet / mu_wet;
604 auto const permeability =
607 .value(vars, pos, t, dt));
610 auto const& b = _process_data.specific_body_force;
614 _process_data.has_gravity
616 -lambda_wet * permeability *
617 (dNdx * (pg_nodal_values - pc_nodal_values) -
620 (pg_nodal_values - pc_nodal_values));
622 laplace_operator.noalias() = dNdx.transpose() * permeability * dNdx * w;
625 auto const solute_dispersivity_transverse =
626 medium.template value<double>(
628 auto const solute_dispersivity_longitudinal =
629 medium.template value<double>(
632 double const velocity_wet_magnitude = velocity_wet.norm();
634 velocity_wet_magnitude != 0.0
636 (porosity * Sw * diffusion_coeff_contaminant_wet +
637 solute_dispersivity_transverse *
638 velocity_wet_magnitude) *
640 (solute_dispersivity_longitudinal -
641 solute_dispersivity_transverse) /
642 velocity_wet_magnitude * velocity_wet *
643 velocity_wet.transpose())
645 (porosity * Sw * diffusion_coeff_contaminant_wet +
646 solute_dispersivity_transverse *
647 velocity_wet_magnitude) *
650 auto const dispersion_operator =
651 dNdx.transpose() * hydrodynamic_dispersion * dNdx * w;
656 (mol_density_nonwet * x_air_nonwet * lambda_nonwet) *
658 (1 - Sw) * porosity * mol_density_nonwet *
659 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpg +
660 diffusion_coeff_contaminant_nonwet *
661 d_x_contaminant_nonwet_dpg) *
664 -(1 - Sw) * porosity * mol_density_nonwet *
665 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpc +
666 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dpc) *
669 -(1 - Sw) * porosity * mol_density_nonwet *
670 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dXc +
671 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dXc) *
674 -(1 - Sw) * porosity * mol_density_nonwet *
675 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dT +
676 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT) *
680 (mol_density_wet * x_water_wet * lambda_wet +
681 mol_density_nonwet * x_water_nonwet * lambda_nonwet) *
684 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
686 (1 - Sw) * diffusion_coeff_water_nonwet * mol_density_nonwet *
687 d_x_water_nonwet_dpg) *
690 -mol_density_wet * x_water_wet * lambda_wet * laplace_operator +
692 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
694 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
695 d_x_water_nonwet_dpc) *
699 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
701 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
702 d_x_water_nonwet_dXc) *
706 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
708 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
709 d_x_water_nonwet_dT) *
713 (mol_density_wet * x_contaminant_wet * lambda_wet +
714 mol_density_nonwet * x_contaminant_nonwet * lambda_nonwet) *
716 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpg +
717 porosity * (1 - Sw) * mol_density_nonwet *
718 diffusion_coeff_contaminant_nonwet *
719 d_x_contaminant_nonwet_dpg * diffusion_operator;
721 -mol_density_wet * x_contaminant_wet * lambda_wet *
723 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpc +
724 porosity * (1 - Sw) * mol_density_nonwet *
725 diffusion_coeff_contaminant_nonwet *
726 d_x_contaminant_nonwet_dpc * diffusion_operator;
728 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dXc +
729 porosity * (1 - Sw) * mol_density_nonwet *
730 diffusion_coeff_contaminant_nonwet *
731 d_x_contaminant_nonwet_dXc * diffusion_operator;
733 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dT +
734 porosity * (1 - Sw) * mol_density_nonwet *
735 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT *
738 Kep.noalias() += (lambda_nonwet * density_nonwet * enthalpy_nonwet +
739 lambda_wet * density_wet * enthalpy_wet) *
742 -lambda_wet * enthalpy_wet * density_wet * laplace_operator;
744 if (medium.hasProperty(
751 .value(vars, pos, t, dt);
757 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
761 auto const thermal_conductivity_solid =
765 .value(vars, pos, t, dt);
767 auto const thermal_conductivity_fluid =
771 .template value<double>(vars, pos, t, dt) *
776 GlobalDim>(thermal_conductivity_solid,
777 thermal_conductivity_fluid, porosity);
780 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
783 if (_process_data.has_gravity)
786 dNdx.transpose() * permeability * b * w;
787 Ba.noalias() += (mol_density_nonwet * x_air_nonwet * lambda_nonwet *
791 (mol_density_wet * x_water_wet * lambda_wet * density_wet +
792 mol_density_nonwet * x_water_nonwet * lambda_nonwet *
795 Bc.noalias() += (mol_density_wet * x_contaminant_wet * lambda_wet *
797 mol_density_nonwet * x_contaminant_nonwet *
798 lambda_nonwet * density_nonwet) *
801 (lambda_nonwet * density_nonwet * density_nonwet *
803 lambda_wet * density_wet * density_wet * enthalpy_wet) *
807 if (_process_data.has_mass_lumping)
809 Map = Map.colwise().sum().eval().asDiagonal();
810 Mapc = Mapc.colwise().sum().eval().asDiagonal();
811 Max = Max.colwise().sum().eval().asDiagonal();
812 Mat = Mat.colwise().sum().eval().asDiagonal();
813 Mwp = Mwp.colwise().sum().eval().asDiagonal();
814 Mwpc = Mwpc.colwise().sum().eval().asDiagonal();
815 Mwx = Mwx.colwise().sum().eval().asDiagonal();
816 Mwt = Mwt.colwise().sum().eval().asDiagonal();
817 Mcp = Mcp.colwise().sum().eval().asDiagonal();
818 Mcpc = Mcpc.colwise().sum().eval().asDiagonal();
819 Mcx = Mcx.colwise().sum().eval().asDiagonal();
820 Mct = Mct.colwise().sum().eval().asDiagonal();
821 Mep = Mep.colwise().sum().eval().asDiagonal();
822 Mepc = Mepc.colwise().sum().eval().asDiagonal();
823 Mex = Mex.colwise().sum().eval().asDiagonal();
824 Met = Met.colwise().sum().eval().asDiagonal();