Loading [MathJax]/extensions/tex2jax.js
OGS
ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
Go to the documentation of this file.
1
11#pragma once
12
21
22namespace ProcessLib
23{
24namespace ThermalTwoPhaseFlowWithPP
25{
26template <typename ShapeFunction, int GlobalDim>
28 assemble(double const t, double const dt,
29 std::vector<double> const& local_x,
30 std::vector<double> const& /*local_x_prev*/,
31 std::vector<double>& local_M_data,
32 std::vector<double>& local_K_data,
33 std::vector<double>& local_b_data)
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>(
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);
58
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);
67
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);
76
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);
85
86 NodalMatrixType laplace_operator =
87 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
88
89 auto Kap =
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);
98
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);
107
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);
116
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);
123
124 auto Ba = local_b.template segment<nonwet_pressure_size>(
125 nonwet_pressure_matrix_index);
126 auto Bw =
127 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
128 auto Bc =
129 local_b.template segment<contaminant_size>(contaminant_matrix_index);
130 auto Be =
131 local_b.template segment<temperature_size>(temperature_matrix_index);
132
133 unsigned const n_integration_points =
134 _integration_method.getNumberOfPoints();
135
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);
141
143
144 GlobalDimMatrixType const& I(
145 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
146
147 for (unsigned ip = 0; ip < n_integration_points; ip++)
148 {
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.;
159 NumLib::shapeFunctionInterpolate(local_x, N, pg_int_pt, pc_int_pt,
160 Xc_int_pt, T_int_pt);
161
163 std::nullopt, _element.getID(),
166 _element, N))};
167
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;
171 vars.temperature = T_int_pt;
172 vars.capillary_pressure = pc_int_pt;
173 vars.gas_phase_pressure = pg_int_pt;
174
175 auto const& medium =
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");
180
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");
186
187 auto const porosity =
189 .template value<double>(vars, pos, t, dt);
190
191 auto const water_mol_mass =
192 water_vapour_component
194 .template value<double>(vars, pos, t, dt);
195 auto const air_mol_mass =
196 dry_air_component
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);
203
204 double const Sw =
206 .template value<double>(vars, pos, t, dt);
207
208 _saturation[ip] = Sw;
209 vars.liquid_saturation = Sw;
210
211 double const dSw_dpc =
213 .template dValue<double>(
215 pos, t, dt);
216
217 auto const density_water =
218 liquid_phase.property(MaterialPropertyLib::PropertyType::density)
219 .template value<double>(vars, pos, t, dt);
220
221 // molar densities and derivatives
222 double const mol_density_water = density_water / water_mol_mass;
223 double const mol_density_wet = mol_density_water;
224
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;
229
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;
238
239 // specific latent heat of evaporation
240 double const latent_heat_evaporation =
241 water_vapour_component
242 .property(
244 .template value<double>(vars, pos, t, dt);
245
246 vars.enthalpy_of_evaporation = latent_heat_evaporation;
247
248 // saturated vapour pressure
249 double const p_sat =
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>(
258 dt);
259
260 // Kelvin-Laplace correction for menisci
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;
265
266 // vapour pressure inside pore space (water partial pressure in gas
267 // phase)
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 =
271 p_vapour_nonwet *
272 (-1 / mol_density_water / ideal_gas_constant_times_T_int_pt);
273
274 // Henry constant of organic contaminant
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>(
284 dt);
285
286 // mass distribution coefficients of contam. and water
287 double const k_c = pg_int_pt * henry_contam / mol_density_wet;
288 double const k_w = pg_int_pt / p_vapour_nonwet;
289
290 // intermediate parameter
291 double const Ntot_c =
292 Sw * mol_density_wet * k_c + (1 - Sw) * mol_density_nonwet;
293 // phase-wise component molar fractions
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;
300
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;
304
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;
308
309 // derivatives of component molar fractions w.r.t. PVs
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);
325
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;
332
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;
337
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) /
341 pg_int_pt;
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) /
349 pg_int_pt;
350
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;
359
360 // mass densities
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;
364
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;
374
375 auto const density_solid =
377 .template value<double>(vars, pos, t, dt);
378
379 // derivatives of nonwet phase densities w.r.t. PVs
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) +
384 mol_density_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 =
389 mol_density_nonwet *
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 =
394 mol_density_nonwet *
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) +
402 mol_density_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);
406
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;
410 // phase-wise component mass fractions
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;
416
417 // spec. heat capacities
418 double const heat_capacity_dry_air =
419 dry_air_component
420 .property(
422 .template value<double>(vars, pos, t, dt);
423 double const heat_capacity_water_vapour =
424 water_vapour_component
425 .property(
427 .template value<double>(vars, pos, t, dt);
428 double const heat_capacity_contaminant_vapour =
429 contaminant_vapour_component
430 .property(
432 .template value<double>(vars, pos, t, dt);
433
434 double const heat_capacity_water =
435 liquid_phase
436 .property(
438 .template value<double>(vars, pos, t, dt);
439 double const heat_capacity_solid =
440 solid_phase
441 .property(
443 .template value<double>(vars, pos, t, dt);
444
445 // enthalpies of gaseous components
446 // Note: h^a_G = C^a_P * (T - T0) = C^a_V * (T - T0) + RT/M^a, thus
447 // "C^a_V" should be used in the following. Same for contaminant.
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;
458
459 // gas and liquid phase enthalpies
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);
466
467 // gas and liquid phase internal energies
468 double const internal_energy_nonwet =
469 enthalpy_nonwet - pg_int_pt / density_nonwet;
470 double const internal_energy_wet = enthalpy_wet;
471
472 // derivatives of enthalpies w.r.t. temperature
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;
478
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;
483
484 // Assemble M matrix
485 Map.noalias() += porosity * (1 - Sw) *
486 (mol_density_nonwet * d_x_air_nonwet_dpg +
487 x_air_nonwet * d_mol_density_nonwet_dpg) *
488 mass_operator;
489 Mapc.noalias() +=
490 porosity * mol_density_nonwet *
491 ((1 - Sw) * d_x_air_nonwet_dpc - x_air_nonwet * dSw_dpc) *
492 mass_operator;
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) *
498 mass_operator;
499
500 Mwp.noalias() +=
501 porosity *
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) *
505 mass_operator;
506 Mwpc.noalias() +=
507 porosity *
508 (mol_density_wet *
509 (x_water_wet * dSw_dpc + Sw * d_x_water_wet_dpc) +
510 mol_density_nonwet *
511 ((1 - Sw) * d_x_water_nonwet_dpc - x_water_nonwet * dSw_dpc)) *
512 mass_operator;
513 Mwx.noalias() +=
514 porosity *
515 (mol_density_wet * Sw * d_x_water_wet_dXc +
516 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dXc) *
517 mass_operator;
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)) *
521 mass_operator;
522
523 Mcp.noalias() +=
524 porosity *
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) *
528 mass_operator;
529 Mcpc.noalias() +=
530 porosity *
531 (mol_density_wet *
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)) *
535 mass_operator;
536 Mcx.noalias() +=
537 porosity *
538 (mol_density_wet * Sw * d_x_contaminant_wet_dXc +
539 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dXc) *
540 mass_operator;
541 Mct.noalias() +=
542 porosity *
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) *
546 mass_operator;
547
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;
559 Met.noalias() +=
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)) *
564 mass_operator;
565
566 // pore diffusion coefficients
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);
579
580 // gas phase relative permeability, viscosity and mobility
581 double const k_rel_nonwet =
582 medium
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;
590
591 // liquid phase relative permeability, viscosity and mobility
592 double const k_rel_wet =
593 medium
594 .property(
596 .template value<double>(vars, pos, t, dt);
597 auto const mu_wet =
599 .template value<double>(vars, pos, t, dt);
600 double const lambda_wet = k_rel_wet / mu_wet;
601
602 // intrinsic permeability
603 auto const permeability =
606 .value(vars, pos, t, dt));
607
608 // gravity
609 auto const& b = _process_data.specific_body_force;
610
611 // gas and liquid phase velocities
612 GlobalDimVectorType const velocity_wet =
613 _process_data.has_gravity
615 -lambda_wet * permeability *
616 (dNdx * (pg_nodal_values - pc_nodal_values) -
617 density_wet * b))
618 : GlobalDimVectorType(-lambda_wet * permeability * dNdx *
619 (pg_nodal_values - pc_nodal_values));
620
621 laplace_operator.noalias() = dNdx.transpose() * permeability * dNdx * w;
622
623 // mechanical dispersivities
624 auto const solute_dispersivity_transverse =
625 medium.template value<double>(
627 auto const solute_dispersivity_longitudinal =
628 medium.template value<double>(
630
631 double const velocity_wet_magnitude = velocity_wet.norm();
632 GlobalDimMatrixType const hydrodynamic_dispersion =
633 velocity_wet_magnitude != 0.0
635 (porosity * Sw * diffusion_coeff_contaminant_wet +
636 solute_dispersivity_transverse *
637 velocity_wet_magnitude) *
638 I +
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) *
647 I);
648
649 auto const dispersion_operator =
650 dNdx.transpose() * hydrodynamic_dispersion * dNdx * w;
651
652 // Assemble K matrix
653 // The sum of all diffusive fluxes in either phase must equal to zero
654 Kap.noalias() +=
655 (mol_density_nonwet * x_air_nonwet * lambda_nonwet) *
656 laplace_operator -
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) *
661 diffusion_operator;
662 Kapc.noalias() +=
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) *
666 diffusion_operator;
667 Kax.noalias() +=
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) *
671 diffusion_operator;
672 Kat.noalias() +=
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) *
676 diffusion_operator;
677
678 Kwp.noalias() +=
679 (mol_density_wet * x_water_wet * lambda_wet +
680 mol_density_nonwet * x_water_nonwet * lambda_nonwet) *
681 laplace_operator +
682 porosity *
683 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
684 d_x_water_wet_dpg +
685 (1 - Sw) * diffusion_coeff_water_nonwet * mol_density_nonwet *
686 d_x_water_nonwet_dpg) *
687 diffusion_operator;
688 Kwpc.noalias() +=
689 -mol_density_wet * x_water_wet * lambda_wet * laplace_operator +
690 porosity *
691 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
692 d_x_water_wet_dpc +
693 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
694 d_x_water_nonwet_dpc) *
695 diffusion_operator;
696 Kwx.noalias() +=
697 porosity *
698 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
699 d_x_water_wet_dXc +
700 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
701 d_x_water_nonwet_dXc) *
702 diffusion_operator;
703 Kwt.noalias() +=
704 porosity *
705 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
706 d_x_water_wet_dT +
707 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
708 d_x_water_nonwet_dT) *
709 diffusion_operator;
710
711 Kcp.noalias() +=
712 (mol_density_wet * x_contaminant_wet * lambda_wet +
713 mol_density_nonwet * x_contaminant_nonwet * lambda_nonwet) *
714 laplace_operator +
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;
719 Kcpc.noalias() +=
720 -mol_density_wet * x_contaminant_wet * lambda_wet *
721 laplace_operator +
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;
726 Kcx.noalias() +=
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;
731 Kct.noalias() +=
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 *
735 diffusion_operator;
736
737 Kep.noalias() += (lambda_nonwet * density_nonwet * enthalpy_nonwet +
738 lambda_wet * density_wet * enthalpy_wet) *
739 laplace_operator;
740 Kepc.noalias() +=
741 -lambda_wet * enthalpy_wet * density_wet * laplace_operator;
742
743 if (medium.hasProperty(
745 {
746 auto const lambda =
747 medium
748 .property(
750 .value(vars, pos, t, dt);
751
752 GlobalDimMatrixType const heat_conductivity_unsaturated =
754
755 Ket.noalias() +=
756 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
757 }
758 else
759 {
760 auto const thermal_conductivity_solid =
761 solid_phase
762 .property(
764 .value(vars, pos, t, dt);
765
766 auto const thermal_conductivity_fluid =
767 liquid_phase
768 .property(
770 .template value<double>(vars, pos, t, dt) *
771 Sw;
772
773 GlobalDimMatrixType const heat_conductivity_unsaturated =
775 GlobalDim>(thermal_conductivity_solid,
776 thermal_conductivity_fluid, porosity);
777
778 Ket.noalias() +=
779 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
780 }
781
782 if (_process_data.has_gravity)
783 {
784 NodalVectorType gravity_operator =
785 dNdx.transpose() * permeability * b * w;
786 Ba.noalias() += (mol_density_nonwet * x_air_nonwet * lambda_nonwet *
787 density_nonwet) *
788 gravity_operator;
789 Bw.noalias() +=
790 (mol_density_wet * x_water_wet * lambda_wet * density_wet +
791 mol_density_nonwet * x_water_nonwet * lambda_nonwet *
792 density_nonwet) *
793 gravity_operator;
794 Bc.noalias() += (mol_density_wet * x_contaminant_wet * lambda_wet *
795 density_wet +
796 mol_density_nonwet * x_contaminant_nonwet *
797 lambda_nonwet * density_nonwet) *
798 gravity_operator;
799 Be.noalias() +=
800 (lambda_nonwet * density_nonwet * density_nonwet *
801 enthalpy_nonwet +
802 lambda_wet * density_wet * density_wet * enthalpy_wet) *
803 gravity_operator;
804 } // end of has gravity
805 }
806 if (_process_data.has_mass_lumping)
807 {
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();
824 } // end of mass-lumping
825}
826
827} // namespace ThermalTwoPhaseFlowWithPP
828} // namespace ProcessLib
Definition of the PiecewiseLinearInterpolation class.
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
constexpr double CelsiusZeroInKelvin
Zero degrees Celsius in Kelvin.
Eigen::Matrix< double, GlobalDim, GlobalDim > formEffectiveThermalConductivity(MaterialPropertyLib::PropertyDataType const &solid_thermal_conductivity, const double fluid_thermal_conductivity, const double porosity)
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 &)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)