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);
250 water_vapour_component
252 .template value<double>(vars, pos, t, dt);
253 double const dp_sat_dT =
254 water_vapour_component
256 .template dValue<double>(
261 double const K = std::exp(-pc_int_pt / mol_density_water /
262 ideal_gas_constant_times_T_int_pt);
263 double const dK_dT = pc_int_pt / mol_density_water /
264 ideal_gas_constant_times_T_int_pt / T_int_pt * K;
268 double const p_vapour_nonwet = p_sat * K;
269 double const d_p_vapour_nonwet_dT = dp_sat_dT * K + p_sat * dK_dT;
270 double const d_p_vapour_nonwet_dpc =
272 (-1 / mol_density_water / ideal_gas_constant_times_T_int_pt);
275 double const henry_contam =
276 contaminant_vapour_component
278 .template value<double>(vars, pos, t, dt);
279 double d_henry_contaminant_dT =
280 contaminant_vapour_component
282 .template dValue<double>(
287 double const k_c = pg_int_pt * henry_contam / mol_density_wet;
288 double const k_w = pg_int_pt / p_vapour_nonwet;
291 double const Ntot_c =
292 Sw * mol_density_wet * k_c + (1 - Sw) * mol_density_nonwet;
294 double const x_contaminant_nonwet =
295 Xc_int_pt * mol_density_tot / Ntot_c;
296 double const x_contaminant_wet = k_c * x_contaminant_nonwet;
297 double const x_water_wet = 1 - x_contaminant_wet;
298 double const x_water_nonwet = x_water_wet / k_w;
299 double const x_air_nonwet = 1 - x_water_nonwet - x_contaminant_nonwet;
301 _gas_molar_fraction_water[ip] = x_water_nonwet;
302 _liquid_molar_fraction_contaminant[ip] = x_contaminant_wet;
303 _gas_molar_fraction_contaminant[ip] = x_contaminant_nonwet;
305 double const d_kc_dpg = henry_contam / mol_density_wet;
306 double const d_kc_dT =
307 pg_int_pt * d_henry_contaminant_dT / mol_density_wet;
310 double const d_x_contaminant_nonwet_dpc =
311 Xc_int_pt * (d_mol_density_tot_dpc / Ntot_c -
312 (mol_density_wet * k_c - mol_density_nonwet) *
313 dSw_dpc * mol_density_tot / Ntot_c / Ntot_c);
314 double const d_x_contaminant_nonwet_dpg =
315 Xc_int_pt * (d_mol_density_tot_dpg / Ntot_c -
316 (Sw * mol_density_wet * d_kc_dpg +
317 (1 - Sw) * d_mol_density_nonwet_dpg) *
318 mol_density_tot / Ntot_c / Ntot_c);
319 double const d_x_contaminant_nonwet_dXc = mol_density_tot / Ntot_c;
320 double const d_x_contaminant_nonwet_dT =
321 Xc_int_pt * (d_mol_density_tot_dT / Ntot_c -
322 (Sw * mol_density_wet * d_kc_dT +
323 (1 - Sw) * d_mol_density_nonwet_dT) *
324 mol_density_tot / Ntot_c / Ntot_c);
326 double const d_x_contaminant_wet_dpc = k_c * d_x_contaminant_nonwet_dpc;
327 double const d_x_contaminant_wet_dpg =
328 k_c * d_x_contaminant_nonwet_dpg + d_kc_dpg * x_contaminant_nonwet;
329 double const d_x_contaminant_wet_dXc = k_c * d_x_contaminant_nonwet_dXc;
330 double const d_x_contaminant_wet_dT =
331 k_c * d_x_contaminant_nonwet_dT + d_kc_dT * x_contaminant_nonwet;
333 double const d_x_water_wet_dpc = -d_x_contaminant_wet_dpc;
334 double const d_x_water_wet_dpg = -d_x_contaminant_wet_dpg;
335 double const d_x_water_wet_dXc = -d_x_contaminant_wet_dXc;
336 double const d_x_water_wet_dT = -d_x_contaminant_wet_dT;
338 double const d_x_water_nonwet_dpc =
339 (d_p_vapour_nonwet_dpc * x_water_wet +
340 p_vapour_nonwet * d_x_water_wet_dpc) /
342 double const d_x_water_nonwet_dpg =
343 p_vapour_nonwet * (d_x_water_wet_dpg / pg_int_pt -
344 x_water_wet / pg_int_pt / pg_int_pt);
345 double const d_x_water_nonwet_dXc = d_x_water_wet_dXc / k_w;
346 double const d_x_water_nonwet_dT =
347 (d_p_vapour_nonwet_dT * x_water_wet +
348 p_vapour_nonwet * d_x_water_wet_dT) /
351 double const d_x_air_nonwet_dpc =
352 -d_x_water_nonwet_dpc - d_x_contaminant_nonwet_dpc;
353 double const d_x_air_nonwet_dpg =
354 -d_x_water_nonwet_dpg - d_x_contaminant_nonwet_dpg;
355 double const d_x_air_nonwet_dXc =
356 -d_x_water_nonwet_dXc - d_x_contaminant_nonwet_dXc;
357 double const d_x_air_nonwet_dT =
358 -d_x_water_nonwet_dT - d_x_contaminant_nonwet_dT;
361 double const density_contaminant_wet =
362 mol_density_wet * contaminant_mol_mass * x_contaminant_wet;
363 double const density_wet = density_water + density_contaminant_wet;
365 double const density_air_nonwet =
366 mol_density_nonwet * air_mol_mass * x_air_nonwet;
367 double const density_water_nonwet =
368 mol_density_nonwet * water_mol_mass * x_water_nonwet;
369 double const density_contaminant_nonwet =
370 mol_density_nonwet * contaminant_mol_mass * x_contaminant_nonwet;
371 double const density_nonwet = density_air_nonwet +
372 density_water_nonwet +
373 density_contaminant_nonwet;
375 auto const density_solid =
377 .template value<double>(vars, pos, t, dt);
380 double const d_density_nonwet_dpg =
381 d_mol_density_nonwet_dpg *
382 (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
383 contaminant_mol_mass * x_contaminant_nonwet) +
385 (air_mol_mass * d_x_air_nonwet_dpg +
386 water_mol_mass * d_x_water_nonwet_dpg +
387 contaminant_mol_mass * d_x_contaminant_nonwet_dpg);
388 double const d_density_nonwet_dpc =
390 (air_mol_mass * d_x_air_nonwet_dpc +
391 water_mol_mass * d_x_water_nonwet_dpc +
392 contaminant_mol_mass * d_x_contaminant_nonwet_dpc);
393 double const d_density_nonwet_dXc =
395 (air_mol_mass * d_x_air_nonwet_dXc +
396 water_mol_mass * d_x_water_nonwet_dXc +
397 contaminant_mol_mass * d_x_contaminant_nonwet_dXc);
398 double const d_density_nonwet_dT =
399 d_mol_density_nonwet_dT *
400 (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
401 contaminant_mol_mass * x_contaminant_nonwet) +
403 (air_mol_mass * d_x_air_nonwet_dT +
404 water_mol_mass * d_x_water_nonwet_dT +
405 contaminant_mol_mass * d_x_contaminant_nonwet_dT);
407 double const mol_mass_nonwet =
408 x_water_nonwet * water_mol_mass + x_air_nonwet * air_mol_mass +
409 x_contaminant_nonwet * contaminant_mol_mass;
411 double const X_water_nonwet =
412 x_water_nonwet * water_mol_mass / mol_mass_nonwet;
413 double const X_air_nonwet =
414 x_air_nonwet * air_mol_mass / mol_mass_nonwet;
415 double const X_contaminant_nonwet = 1 - X_water_nonwet - X_air_nonwet;
418 double const heat_capacity_dry_air =
422 .template value<double>(vars, pos, t, dt);
423 double const heat_capacity_water_vapour =
424 water_vapour_component
427 .template value<double>(vars, pos, t, dt);
428 double const heat_capacity_contaminant_vapour =
429 contaminant_vapour_component
432 .template value<double>(vars, pos, t, dt);
434 double const heat_capacity_water =
438 .template value<double>(vars, pos, t, dt);
439 double const heat_capacity_solid =
443 .template value<double>(vars, pos, t, dt);
448 double const enthalpy_air_nonwet =
449 heat_capacity_dry_air * (T_int_pt - CelsiusZeroInKelvin) +
450 IdealGasConstant * T_int_pt / air_mol_mass;
451 double const enthalpy_water_nonwet =
452 heat_capacity_water_vapour * (T_int_pt - CelsiusZeroInKelvin) +
453 latent_heat_evaporation;
454 double const enthalpy_contaminant_nonwet =
455 heat_capacity_contaminant_vapour *
456 (T_int_pt - CelsiusZeroInKelvin) +
457 IdealGasConstant * T_int_pt / contaminant_mol_mass;
460 double const enthalpy_nonwet =
461 enthalpy_air_nonwet * X_air_nonwet +
462 enthalpy_water_nonwet * X_water_nonwet +
463 enthalpy_contaminant_nonwet * X_contaminant_nonwet;
464 double const enthalpy_wet =
465 heat_capacity_water * (T_int_pt - CelsiusZeroInKelvin);
468 double const internal_energy_nonwet =
469 enthalpy_nonwet - pg_int_pt / density_nonwet;
470 double const internal_energy_wet = enthalpy_wet;
473 double const d_enthalpy_air_nonwet_dT =
474 heat_capacity_dry_air + IdealGasConstant / air_mol_mass;
475 double const d_enthalpy_contaminant_nonwet_dT =
476 heat_capacity_contaminant_vapour +
477 IdealGasConstant / contaminant_mol_mass;
479 double const d_enthalpy_nonwet_dT =
480 heat_capacity_water * X_water_nonwet +
481 d_enthalpy_air_nonwet_dT * X_air_nonwet +
482 d_enthalpy_contaminant_nonwet_dT * X_contaminant_nonwet;
485 Map.noalias() += porosity * (1 - Sw) *
486 (mol_density_nonwet * d_x_air_nonwet_dpg +
487 x_air_nonwet * d_mol_density_nonwet_dpg) *
490 porosity * mol_density_nonwet *
491 ((1 - Sw) * d_x_air_nonwet_dpc - x_air_nonwet * dSw_dpc) *
493 Max.noalias() += porosity * mol_density_nonwet * (1 - Sw) *
494 d_x_air_nonwet_dXc * mass_operator;
495 Mat.noalias() += porosity * (1 - Sw) *
496 (mol_density_nonwet * d_x_air_nonwet_dT +
497 x_air_nonwet * d_mol_density_nonwet_dT) *
502 (mol_density_wet * Sw * d_x_water_wet_dpg +
503 (1 - Sw) * x_water_nonwet * d_mol_density_nonwet_dpg +
504 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dpg) *
509 (x_water_wet * dSw_dpc + Sw * d_x_water_wet_dpc) +
511 ((1 - Sw) * d_x_water_nonwet_dpc - x_water_nonwet * dSw_dpc)) *
515 (mol_density_wet * Sw * d_x_water_wet_dXc +
516 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dXc) *
518 Mwt.noalias() += porosity *
519 ((1 - Sw) / ideal_gas_constant_times_T_int_pt *
520 (d_p_vapour_nonwet_dT - p_vapour_nonwet / T_int_pt)) *
525 (mol_density_wet * Sw * d_x_contaminant_wet_dpg +
526 (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dpg +
527 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dpg) *
532 (x_contaminant_wet * dSw_dpc + Sw * d_x_contaminant_wet_dpc) +
533 mol_density_nonwet * ((1 - Sw) * d_x_contaminant_nonwet_dpc -
534 x_contaminant_nonwet * dSw_dpc)) *
538 (mol_density_wet * Sw * d_x_contaminant_wet_dXc +
539 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dXc) *
543 (mol_density_wet * Sw * d_x_contaminant_wet_dT +
544 (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dT +
545 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dT) *
548 Mep.noalias() += porosity *
549 (d_density_nonwet_dpg * enthalpy_nonwet - 1) *
550 (1 - Sw) * mass_operator;
551 Mepc.noalias() += porosity *
552 (density_wet * internal_energy_wet -
553 density_nonwet * internal_energy_nonwet) *
554 dSw_dpc * mass_operator +
555 porosity * d_density_nonwet_dpc * enthalpy_nonwet *
556 (1 - Sw) * mass_operator;
557 Mex.noalias() += porosity * d_density_nonwet_dXc * enthalpy_nonwet *
558 (1 - Sw) * mass_operator;
560 ((1 - porosity) * density_solid * heat_capacity_solid +
561 porosity * ((1 - Sw) * (d_density_nonwet_dT * enthalpy_nonwet +
562 density_nonwet * d_enthalpy_nonwet_dT) +
563 Sw * density_wet * heat_capacity_water)) *
567 double const diffusion_coeff_water_nonwet =
568 water_vapour_component
570 .template value<double>(vars, pos, t, dt);
571 double const diffusion_coeff_contaminant_nonwet =
572 contaminant_vapour_component
574 .template value<double>(vars, pos, t, dt);
575 double const diffusion_coeff_contaminant_wet =
576 dissolved_contaminant_component
578 .template value<double>(vars, pos, t, dt);
581 double const k_rel_nonwet =
584 relative_permeability_nonwetting_phase)
585 .template value<double>(vars, pos, t, dt);
586 auto const mu_nonwet =
588 .template value<double>(vars, pos, t, dt);
589 double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
592 double const k_rel_wet =
596 .template value<double>(vars, pos, t, dt);
599 .template value<double>(vars, pos, t, dt);
600 double const lambda_wet = k_rel_wet / mu_wet;
603 auto const permeability =
606 .value(vars, pos, t, dt));
609 auto const& b = _process_data.specific_body_force;
613 _process_data.has_gravity
615 -lambda_wet * permeability *
616 (dNdx * (pg_nodal_values - pc_nodal_values) -
619 (pg_nodal_values - pc_nodal_values));
621 laplace_operator.noalias() = dNdx.transpose() * permeability * dNdx * w;
624 auto const solute_dispersivity_transverse =
625 medium.template value<double>(
627 auto const solute_dispersivity_longitudinal =
628 medium.template value<double>(
631 double const velocity_wet_magnitude = velocity_wet.norm();
633 velocity_wet_magnitude != 0.0
635 (porosity * Sw * diffusion_coeff_contaminant_wet +
636 solute_dispersivity_transverse *
637 velocity_wet_magnitude) *
639 (solute_dispersivity_longitudinal -
640 solute_dispersivity_transverse) /
641 velocity_wet_magnitude * velocity_wet *
642 velocity_wet.transpose())
644 (porosity * Sw * diffusion_coeff_contaminant_wet +
645 solute_dispersivity_transverse *
646 velocity_wet_magnitude) *
649 auto const dispersion_operator =
650 dNdx.transpose() * hydrodynamic_dispersion * dNdx * w;
655 (mol_density_nonwet * x_air_nonwet * lambda_nonwet) *
657 (1 - Sw) * porosity * mol_density_nonwet *
658 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpg +
659 diffusion_coeff_contaminant_nonwet *
660 d_x_contaminant_nonwet_dpg) *
663 -(1 - Sw) * porosity * mol_density_nonwet *
664 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpc +
665 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dpc) *
668 -(1 - Sw) * porosity * mol_density_nonwet *
669 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dXc +
670 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dXc) *
673 -(1 - Sw) * porosity * mol_density_nonwet *
674 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dT +
675 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT) *
679 (mol_density_wet * x_water_wet * lambda_wet +
680 mol_density_nonwet * x_water_nonwet * lambda_nonwet) *
683 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
685 (1 - Sw) * diffusion_coeff_water_nonwet * mol_density_nonwet *
686 d_x_water_nonwet_dpg) *
689 -mol_density_wet * x_water_wet * lambda_wet * laplace_operator +
691 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
693 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
694 d_x_water_nonwet_dpc) *
698 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
700 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
701 d_x_water_nonwet_dXc) *
705 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
707 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
708 d_x_water_nonwet_dT) *
712 (mol_density_wet * x_contaminant_wet * lambda_wet +
713 mol_density_nonwet * x_contaminant_nonwet * lambda_nonwet) *
715 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpg +
716 porosity * (1 - Sw) * mol_density_nonwet *
717 diffusion_coeff_contaminant_nonwet *
718 d_x_contaminant_nonwet_dpg * diffusion_operator;
720 -mol_density_wet * x_contaminant_wet * lambda_wet *
722 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpc +
723 porosity * (1 - Sw) * mol_density_nonwet *
724 diffusion_coeff_contaminant_nonwet *
725 d_x_contaminant_nonwet_dpc * diffusion_operator;
727 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dXc +
728 porosity * (1 - Sw) * mol_density_nonwet *
729 diffusion_coeff_contaminant_nonwet *
730 d_x_contaminant_nonwet_dXc * diffusion_operator;
732 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dT +
733 porosity * (1 - Sw) * mol_density_nonwet *
734 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT *
737 Kep.noalias() += (lambda_nonwet * density_nonwet * enthalpy_nonwet +
738 lambda_wet * density_wet * enthalpy_wet) *
741 -lambda_wet * enthalpy_wet * density_wet * laplace_operator;
743 if (medium.hasProperty(
750 .value(vars, pos, t, dt);
756 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
760 auto const thermal_conductivity_solid =
764 .value(vars, pos, t, dt);
766 auto const thermal_conductivity_fluid =
770 .template value<double>(vars, pos, t, dt) *
775 GlobalDim>(thermal_conductivity_solid,
776 thermal_conductivity_fluid, porosity);
779 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
782 if (_process_data.has_gravity)
785 dNdx.transpose() * permeability * b * w;
786 Ba.noalias() += (mol_density_nonwet * x_air_nonwet * lambda_nonwet *
790 (mol_density_wet * x_water_wet * lambda_wet * density_wet +
791 mol_density_nonwet * x_water_nonwet * lambda_nonwet *
794 Bc.noalias() += (mol_density_wet * x_contaminant_wet * lambda_wet *
796 mol_density_nonwet * x_contaminant_nonwet *
797 lambda_nonwet * density_nonwet) *
800 (lambda_nonwet * density_nonwet * density_nonwet *
802 lambda_wet * density_wet * density_wet * enthalpy_wet) *
806 if (_process_data.has_mass_lumping)
808 Map = Map.colwise().sum().eval().asDiagonal();
809 Mapc = Mapc.colwise().sum().eval().asDiagonal();
810 Max = Max.colwise().sum().eval().asDiagonal();
811 Mat = Mat.colwise().sum().eval().asDiagonal();
812 Mwp = Mwp.colwise().sum().eval().asDiagonal();
813 Mwpc = Mwpc.colwise().sum().eval().asDiagonal();
814 Mwx = Mwx.colwise().sum().eval().asDiagonal();
815 Mwt = Mwt.colwise().sum().eval().asDiagonal();
816 Mcp = Mcp.colwise().sum().eval().asDiagonal();
817 Mcpc = Mcpc.colwise().sum().eval().asDiagonal();
818 Mcx = Mcx.colwise().sum().eval().asDiagonal();
819 Mct = Mct.colwise().sum().eval().asDiagonal();
820 Mep = Mep.colwise().sum().eval().asDiagonal();
821 Mepc = Mepc.colwise().sum().eval().asDiagonal();
822 Mex = Mex.colwise().sum().eval().asDiagonal();
823 Met = Met.colwise().sum().eval().asDiagonal();