34{
37
38 auto const local_matrix_size = local_x.size();
39
40 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
41
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);
48
49 auto Map =
50 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
52 auto Mapc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
54 auto Max = local_M.template block<nonwet_pressure_size, contaminant_size>(
56 auto Mat = local_M.template block<nonwet_pressure_size, temperature_size>(
58
59 auto Mwp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
61 auto Mwpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
63 auto Mwx = local_M.template block<cap_pressure_size, contaminant_size>(
65 auto Mwt = local_M.template block<cap_pressure_size, temperature_size>(
67
68 auto Mcp = local_M.template block<contaminant_size, nonwet_pressure_size>(
70 auto Mcpc = local_M.template block<contaminant_size, cap_pressure_size>(
72 auto Mcx = local_M.template block<contaminant_size, contaminant_size>(
74 auto Mct = local_M.template block<contaminant_size, temperature_size>(
76
77 auto Mep = local_M.template block<temperature_size, nonwet_pressure_size>(
79 auto Mepc = local_M.template block<temperature_size, cap_pressure_size>(
81 auto Mex = local_M.template block<temperature_size, contaminant_size>(
83 auto Met = local_M.template block<temperature_size, temperature_size>(
85
87 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
88
89 auto Kap =
90 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
92 auto Kapc = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
94 auto Kax = local_K.template block<nonwet_pressure_size, contaminant_size>(
96 auto Kat = local_K.template block<nonwet_pressure_size, temperature_size>(
98
99 auto Kwp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
101 auto Kwpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
103 auto Kwx = local_K.template block<cap_pressure_size, contaminant_size>(
105 auto Kwt = local_K.template block<cap_pressure_size, temperature_size>(
107
108 auto Kcp = local_K.template block<contaminant_size, nonwet_pressure_size>(
110 auto Kcpc = local_K.template block<contaminant_size, cap_pressure_size>(
112 auto Kcx = local_K.template block<contaminant_size, contaminant_size>(
114 auto Kct = local_K.template block<contaminant_size, temperature_size>(
116
117 auto Kep = local_K.template block<temperature_size, nonwet_pressure_size>(
119 auto Kepc = local_K.template block<temperature_size, cap_pressure_size>(
121 auto Ket = local_K.template block<temperature_size, temperature_size>(
123
124 auto Ba = local_b.template segment<nonwet_pressure_size>(
126 auto Bw =
128 auto Bc =
130 auto Be =
132
133 unsigned const n_integration_points =
135
138
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);
144
146
148 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
149
150 for (unsigned ip = 0; ip < n_integration_points; ip++)
151 {
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);
164
166 double const ideal_gas_constant_times_T_int_pt =
167 IdealGasConstant * T_int_pt;
171
172 auto const& medium =
174 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
175 auto const& solid_phase = medium.phase("Solid");
176 auto const& gas_phase = medium.phase("Gas");
177
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");
183
186 .template value<double>(vars, pos, t, dt);
187
188 auto const water_mol_mass =
189 water_vapour_component
191 .template value<double>(vars, pos, t, dt);
192 auto const air_mol_mass =
193 dry_air_component
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);
200
201 double const Sw =
203 .template value<double>(vars, pos, t, dt);
204
207
208 double const dSw_dpc =
210 .template dValue<double>(
212 pos, t, dt);
213
214 auto const density_water =
216 .template value<double>(vars, pos, t, dt);
217
218
219 double const mol_density_water = density_water / water_mol_mass;
220 double const mol_density_wet = mol_density_water;
221
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;
226
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;
235
236
237 double const latent_heat_evaporation =
238 water_vapour_component
239 .property(
241 .template value<double>(vars, pos, t, dt);
242
244
245
246 double const p_sat =
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>(
255 dt);
256
257
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;
262
263
264
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 =
268 p_vapour_nonwet *
269 (-1 / mol_density_water / ideal_gas_constant_times_T_int_pt);
270
271
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>(
281 dt);
282
283
284 double const k_c = pg_int_pt * henry_contam / mol_density_wet;
285 double const k_w = pg_int_pt / p_vapour_nonwet;
286
287
288 double const Ntot_c =
289 Sw * mol_density_wet * k_c + (1 - Sw) * mol_density_nonwet;
290
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;
297
301
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;
305
306
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);
322
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;
329
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;
334
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) /
338 pg_int_pt;
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) /
346 pg_int_pt;
347
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;
356
357
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;
361
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;
371
372 auto const density_solid =
374 .template value<double>(vars, pos, t, dt);
375
376
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) +
381 mol_density_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 =
386 mol_density_nonwet *
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 =
391 mol_density_nonwet *
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) +
399 mol_density_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);
403
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;
407
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;
413
414
415 double const heat_capacity_dry_air =
416 dry_air_component
417 .property(
419 .template value<double>(vars, pos, t, dt);
420 double const heat_capacity_water_vapour =
421 water_vapour_component
422 .property(
424 .template value<double>(vars, pos, t, dt);
425 double const heat_capacity_contaminant_vapour =
426 contaminant_vapour_component
427 .property(
429 .template value<double>(vars, pos, t, dt);
430
431 double const heat_capacity_water =
432 liquid_phase
433 .property(
435 .template value<double>(vars, pos, t, dt);
436 double const heat_capacity_solid =
437 solid_phase
438 .property(
440 .template value<double>(vars, pos, t, dt);
441
442
443
444
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;
455
456
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);
463
464
465 double const internal_energy_nonwet =
466 enthalpy_nonwet - pg_int_pt / density_nonwet;
467 double const internal_energy_wet = enthalpy_wet;
468
469
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;
475
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;
480
481
482 Map.noalias() +=
porosity * (1 - Sw) *
483 (mol_density_nonwet * d_x_air_nonwet_dpg +
484 x_air_nonwet * d_mol_density_nonwet_dpg) *
485 mass_operator;
486 Mapc.noalias() +=
488 ((1 - Sw) * d_x_air_nonwet_dpc - x_air_nonwet * dSw_dpc) *
489 mass_operator;
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) *
495 mass_operator;
496
497 Mwp.noalias() +=
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) *
502 mass_operator;
503 Mwpc.noalias() +=
505 (mol_density_wet *
506 (x_water_wet * dSw_dpc + Sw * d_x_water_wet_dpc) +
507 mol_density_nonwet *
508 ((1 - Sw) * d_x_water_nonwet_dpc - x_water_nonwet * dSw_dpc)) *
509 mass_operator;
510 Mwx.noalias() +=
512 (mol_density_wet * Sw * d_x_water_wet_dXc +
513 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dXc) *
514 mass_operator;
516 ((1 - Sw) / ideal_gas_constant_times_T_int_pt *
517 (d_p_vapour_nonwet_dT - p_vapour_nonwet / T_int_pt)) *
518 mass_operator;
519
520 Mcp.noalias() +=
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) *
525 mass_operator;
526 Mcpc.noalias() +=
528 (mol_density_wet *
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)) *
532 mass_operator;
533 Mcx.noalias() +=
535 (mol_density_wet * Sw * d_x_contaminant_wet_dXc +
536 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dXc) *
537 mass_operator;
538 Mct.noalias() +=
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) *
543 mass_operator;
544
546 (d_density_nonwet_dpg * enthalpy_nonwet - 1) *
547 (1 - Sw) * mass_operator;
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;
556 Met.noalias() +=
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)) *
561 mass_operator;
562
563
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);
576
577
578 double const k_rel_nonwet =
579 medium
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;
587
588
589 double const k_rel_wet =
590 medium
591 .property(
593 .template value<double>(vars, pos, t, dt);
594 auto const mu_wet =
596 .template value<double>(vars, pos, t, dt);
597 double const lambda_wet = k_rel_wet / mu_wet;
598
599
603 .value(vars, pos, t, dt));
604
605
607
608
612 -lambda_wet * permeability *
613 (dNdx * (pg_nodal_values - pc_nodal_values) -
614 density_wet * b))
616 (pg_nodal_values - pc_nodal_values));
617
618 laplace_operator.noalias() = dNdx.transpose() * permeability * dNdx * w;
619
620
621 auto const solute_dispersivity_transverse =
622 medium.template value<double>(
624 auto const solute_dispersivity_longitudinal =
625 medium.template value<double>(
627
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) *
635 I +
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) *
644 I);
645
646 auto const dispersion_operator =
647 dNdx.transpose() * hydrodynamic_dispersion * dNdx * w;
648
649
650
651 Kap.noalias() +=
652 (mol_density_nonwet * x_air_nonwet * lambda_nonwet) *
653 laplace_operator -
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) *
658 diffusion_operator;
659 Kapc.noalias() +=
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) *
663 diffusion_operator;
664 Kax.noalias() +=
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) *
668 diffusion_operator;
669 Kat.noalias() +=
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) *
673 diffusion_operator;
674
675 Kwp.noalias() +=
676 (mol_density_wet * x_water_wet * lambda_wet +
677 mol_density_nonwet * x_water_nonwet * lambda_nonwet) *
678 laplace_operator +
679 porosity *
680 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
681 d_x_water_wet_dpg +
682 (1 - Sw) * diffusion_coeff_water_nonwet * mol_density_nonwet *
683 d_x_water_nonwet_dpg) *
684 diffusion_operator;
685 Kwpc.noalias() +=
686 -mol_density_wet * x_water_wet * lambda_wet * laplace_operator +
687 porosity *
688 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
689 d_x_water_wet_dpc +
690 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
691 d_x_water_nonwet_dpc) *
692 diffusion_operator;
693 Kwx.noalias() +=
694 porosity *
695 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
696 d_x_water_wet_dXc +
697 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
698 d_x_water_nonwet_dXc) *
699 diffusion_operator;
700 Kwt.noalias() +=
701 porosity *
702 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
703 d_x_water_wet_dT +
704 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
705 d_x_water_nonwet_dT) *
706 diffusion_operator;
707
708 Kcp.noalias() +=
709 (mol_density_wet * x_contaminant_wet * lambda_wet +
710 mol_density_nonwet * x_contaminant_nonwet * lambda_nonwet) *
711 laplace_operator +
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;
716 Kcpc.noalias() +=
717 -mol_density_wet * x_contaminant_wet * lambda_wet *
718 laplace_operator +
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;
723 Kcx.noalias() +=
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;
728 Kct.noalias() +=
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 *
732 diffusion_operator;
733
734 Kep.noalias() += (lambda_nonwet * density_nonwet * enthalpy_nonwet +
735 lambda_wet * density_wet * enthalpy_wet) *
736 laplace_operator;
737 Kepc.noalias() +=
738 -lambda_wet * enthalpy_wet * density_wet * laplace_operator;
739
740 if (medium.hasProperty(
742 {
743 auto const lambda =
744 medium
745 .property(
746 MaterialPropertyLib::PropertyType::thermal_conductivity)
747 .value(vars, pos, t, dt);
748
749 GlobalDimMatrixType const heat_conductivity_unsaturated =
750 MaterialPropertyLib::formEigenTensor<GlobalDim>(lambda);
751
752 Ket.noalias() +=
753 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
754 }
755 else
756 {
757 auto const thermal_conductivity_solid =
758 solid_phase
759 .property(
760 MaterialPropertyLib::PropertyType::thermal_conductivity)
761 .value(vars, pos, t, dt);
762
763 auto const thermal_conductivity_fluid =
764 liquid_phase
765 .property(
766 MaterialPropertyLib::PropertyType::thermal_conductivity)
767 .template value<double>(vars, pos, t, dt) *
768 Sw;
769
770 GlobalDimMatrixType const heat_conductivity_unsaturated =
771 MaterialPropertyLib::formEffectiveThermalConductivity<
772 GlobalDim>(thermal_conductivity_solid,
773 thermal_conductivity_fluid, porosity);
774
775 Ket.noalias() +=
776 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
777 }
778
780 {
781 NodalVectorType gravity_operator =
782 dNdx.transpose() * permeability * b * w;
783 Ba.noalias() += (mol_density_nonwet * x_air_nonwet * lambda_nonwet *
784 density_nonwet) *
785 gravity_operator;
786 Bw.noalias() +=
787 (mol_density_wet * x_water_wet * lambda_wet * density_wet +
788 mol_density_nonwet * x_water_nonwet * lambda_nonwet *
789 density_nonwet) *
790 gravity_operator;
791 Bc.noalias() += (mol_density_wet * x_contaminant_wet * lambda_wet *
792 density_wet +
793 mol_density_nonwet * x_contaminant_nonwet *
794 lambda_nonwet * density_nonwet) *
795 gravity_operator;
796 Be.noalias() +=
797 (lambda_nonwet * density_nonwet * density_nonwet *
798 enthalpy_nonwet +
799 lambda_wet * density_wet * density_wet * enthalpy_wet) *
800 gravity_operator;
801 }
802 }
804 {
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();
821 }
822}
Medium * getMedium(std::size_t element_id)
Phase const & phase(std::size_t index) const
Component const & component(std::size_t const &index) const
double gas_phase_pressure
double capillary_pressure
double enthalpy_of_evaporation
std::size_t getID() const
Returns the ID of the element.
void setElementID(std::size_t element_id)
static const int temperature_matrix_index
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
static const int nonwet_pressure_matrix_index
static const int contaminant_matrix_index
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
static const int cap_pressure_matrix_index
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
constexpr double CelsiusZeroInKelvin
Zero degrees Celsius in Kelvin.
constexpr double IdealGasConstant
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
bool const has_mass_lumping
Eigen::VectorXd const specific_body_force
MaterialPropertyLib::MaterialSpatialDistributionMap media_map