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
42 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
43 local_M_data, local_matrix_size, local_matrix_size);
44 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
45 local_K_data, local_matrix_size, local_matrix_size);
46 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
47 local_b_data, local_matrix_size);
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
137 pos.setElementID(_element.getID());
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
147 GlobalDimMatrixType const& I(
148 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
149
150 for (unsigned ip = 0; ip < n_integration_points; ip++)
151 {
152 auto const& ip_data = _ip_data[ip];
153 auto const& N = ip_data.N;
154 auto const& dNdx = ip_data.dNdx;
155 auto const& w = ip_data.integration_weight;
156 auto const& mass_operator = ip_data.mass_operator;
157 auto const& diffusion_operator = ip_data.diffusion_operator;
158 double pg_int_pt = 0.;
159 double pc_int_pt = 0.;
160 double Xc_int_pt = 0.;
161 double T_int_pt = 0.;
162 NumLib::shapeFunctionInterpolate(local_x, N, pg_int_pt, pc_int_pt,
163 Xc_int_pt, T_int_pt);
164
165 _pressure_wetting[ip] = pg_int_pt - pc_int_pt;
166 double const ideal_gas_constant_times_T_int_pt =
167 IdealGasConstant * T_int_pt;
168 vars.temperature = T_int_pt;
169 vars.capillary_pressure = pc_int_pt;
170 vars.gas_phase_pressure = pg_int_pt;
171
172 auto const& medium =
173 *_process_data.media_map.getMedium(this->_element.getID());
174 auto const& liquid_phase = medium.phase("AqueousLiquid");
175 auto const& solid_phase = medium.phase("Solid");
176 auto const& gas_phase = medium.phase("Gas");
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
184 auto const porosity =
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
205 _saturation[ip] = Sw;
206 vars.liquid_saturation = Sw;
207
208 double const dSw_dpc =
210 .template dValue<double>(
212 pos, t, dt);
213
214 auto const density_water =
215 liquid_phase.property(MaterialPropertyLib::PropertyType::density)
216 .template value<double>(vars, pos, t, dt);
217
218 // molar densities and derivatives
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 // specific latent heat of evaporation
237 double const latent_heat_evaporation =
238 water_vapour_component
239 .property(
241 .template value<double>(vars, pos, t, dt);
242
243 vars.enthalpy_of_evaporation = latent_heat_evaporation;
244
245 // saturated vapour pressure
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 // Kelvin-Laplace correction for menisci
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 // vapour pressure inside pore space (water partial pressure in gas
264 // phase)
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 // Henry constant of organic contaminant
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 // mass distribution coefficients of contam. and water
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 // intermediate parameter
288 double const Ntot_c =
289 Sw * mol_density_wet * k_c + (1 - Sw) * mol_density_nonwet;
290 // phase-wise component molar fractions
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
298 _gas_molar_fraction_water[ip] = x_water_nonwet;
299 _liquid_molar_fraction_contaminant[ip] = x_contaminant_wet;
300 _gas_molar_fraction_contaminant[ip] = x_contaminant_nonwet;
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 // derivatives of component molar fractions w.r.t. PVs
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 // mass densities
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 // derivatives of nonwet phase densities w.r.t. PVs
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 // phase-wise component mass fractions
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 // spec. heat capacities
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 // enthalpies of gaseous components
443 // Note: h^a_G = C^a_P * (T - T0) = C^a_V * (T - T0) + RT/M^a, thus
444 // "C^a_V" should be used in the following. Same for contaminant.
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 // gas and liquid phase enthalpies
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 // gas and liquid phase internal energies
465 double const internal_energy_nonwet =
466 enthalpy_nonwet - pg_int_pt / density_nonwet;
467 double const internal_energy_wet = enthalpy_wet;
468
469 // derivatives of enthalpies w.r.t. temperature
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 // Assemble M matrix
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() +=
487 porosity * mol_density_nonwet *
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() +=
498 porosity *
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() +=
504 porosity *
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() +=
511 porosity *
512 (mol_density_wet * Sw * d_x_water_wet_dXc +
513 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dXc) *
514 mass_operator;
515 Mwt.noalias() += porosity *
516 ((1 - Sw) / ideal_gas_constant_times_T_int_pt *
517 (d_p_vapour_nonwet_dT - p_vapour_nonwet / T_int_pt)) *
518 mass_operator;
519
520 Mcp.noalias() +=
521 porosity *
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() +=
527 porosity *
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() +=
534 porosity *
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() +=
539 porosity *
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
545 Mep.noalias() += porosity *
546 (d_density_nonwet_dpg * enthalpy_nonwet - 1) *
547 (1 - Sw) * mass_operator;
548 Mepc.noalias() += porosity *
549 (density_wet * internal_energy_wet -
550 density_nonwet * internal_energy_nonwet) *
551 dSw_dpc * mass_operator +
552 porosity * d_density_nonwet_dpc * enthalpy_nonwet *
553 (1 - Sw) * mass_operator;
554 Mex.noalias() += porosity * d_density_nonwet_dXc * enthalpy_nonwet *
555 (1 - Sw) * mass_operator;
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 // pore diffusion coefficients
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 // gas phase relative permeability, viscosity and mobility
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 // liquid phase relative permeability, viscosity and mobility
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 // intrinsic permeability
600 auto const permeability =
601 MaterialPropertyLib::formEigenTensor<GlobalDim>(
603 .value(vars, pos, t, dt));
604
605 // gravity
606 auto const& b = _process_data.specific_body_force;
607
608 // gas and liquid phase velocities
609 GlobalDimVectorType const velocity_wet =
610 _process_data.has_gravity
612 -lambda_wet * permeability *
613 (dNdx * (pg_nodal_values - pc_nodal_values) -
614 density_wet * b))
615 : GlobalDimVectorType(-lambda_wet * permeability * dNdx *
616 (pg_nodal_values - pc_nodal_values));
617
618 laplace_operator.noalias() = dNdx.transpose() * permeability * dNdx * w;
619
620 // mechanical dispersivities
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();
629 GlobalDimMatrixType const hydrodynamic_dispersion =
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 // Assemble K matrix
650 // The sum of all diffusive fluxes in either phase must equal to zero
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(
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(
761 .value(vars, pos, t, dt);
762
763 auto const thermal_conductivity_fluid =
764 liquid_phase
765 .property(
767 .template value<double>(vars, pos, t, dt) *
768 Sw;
769
770 GlobalDimMatrixType const heat_conductivity_unsaturated =
772 GlobalDim>(thermal_conductivity_solid,
773 thermal_conductivity_fluid, porosity);
774
775 Ket.noalias() +=
776 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
777 }
778
779 if (_process_data.has_gravity)
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 } // end of has gravity
802 }
803 if (_process_data.has_mass_lumping)
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 } // end of mass-lumping
822}
823
824} // namespace ThermalTwoPhaseFlowWithPP
825} // namespace ProcessLib
Definition of the PiecewiseLinearInterpolation class.
void setElementID(std::size_t element_id)
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)
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)