55 namespace ThermalTwoPhaseFlowWithPP
57 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
58 void ThermalTwoPhaseFlowWithPPLocalAssembler<
59 ShapeFunction, IntegrationMethod,
60 GlobalDim>::assemble(
double const t,
double const ,
61 std::vector<double>
const& local_x,
62 std::vector<double>
const& ,
63 std::vector<double>& local_M_data,
64 std::vector<double>& local_K_data,
65 std::vector<double>& local_b_data)
68 auto const& water_mol_mass =
72 auto const local_matrix_size = local_x.size();
74 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
76 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
77 local_M_data, local_matrix_size, local_matrix_size);
78 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
79 local_K_data, local_matrix_size, local_matrix_size);
80 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
81 local_b_data, local_matrix_size);
84 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
85 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
86 auto Mgpc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
87 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
88 auto Mgt = local_M.template block<nonwet_pressure_size, temperature_size>(
89 nonwet_pressure_matrix_index, temperature_matrix_index);
91 auto Mlp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
92 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
93 auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
94 cap_pressure_matrix_index, cap_pressure_matrix_index);
95 auto Mlt = local_M.template block<cap_pressure_size, temperature_size>(
96 cap_pressure_matrix_index, temperature_matrix_index);
98 auto Mep = local_M.template block<temperature_size, nonwet_pressure_size>(
99 temperature_matrix_index, nonwet_pressure_matrix_index);
100 auto Mepc = local_M.template block<temperature_size, cap_pressure_size>(
101 temperature_matrix_index, cap_pressure_matrix_index);
102 auto Met = local_M.template block<temperature_size, temperature_size>(
103 temperature_matrix_index, temperature_matrix_index);
106 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
109 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
110 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
111 auto Kgpc = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
112 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
113 auto Kgt = local_K.template block<nonwet_pressure_size, temperature_size>(
114 nonwet_pressure_matrix_index, temperature_matrix_index);
116 auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
117 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
118 auto Klpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
119 cap_pressure_matrix_index, cap_pressure_matrix_index);
120 auto Klt = local_K.template block<cap_pressure_size, temperature_size>(
121 cap_pressure_matrix_index, temperature_matrix_index);
123 auto Kep = local_K.template block<temperature_size, nonwet_pressure_size>(
124 temperature_matrix_index, nonwet_pressure_matrix_index);
125 auto Kepc = local_K.template block<temperature_size, cap_pressure_size>(
126 temperature_matrix_index, cap_pressure_matrix_index);
127 auto Ket = local_K.template block<temperature_size, temperature_size>(
128 temperature_matrix_index, temperature_matrix_index);
130 auto Bg = local_b.template segment<nonwet_pressure_size>(
131 nonwet_pressure_matrix_index);
133 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
135 local_b.template segment<temperature_size>(temperature_matrix_index);
137 unsigned const n_integration_points =
138 _integration_method.getNumberOfPoints();
142 auto const& two_phase_material_model =
143 _process_data.material->getTwoPhaseMaterialModel();
144 const int material_id =
145 two_phase_material_model.getMaterialID(pos.
getElementID().value());
147 auto const num_nodes = ShapeFunction::NPOINTS;
148 auto const pg_nodal_values =
149 Eigen::Map<const NodalVectorType>(&local_x[0], num_nodes);
150 auto const pc_nodal_values =
151 Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
153 const Eigen::MatrixXd& perm = two_phase_material_model.getPermeability(
154 material_id, t, pos, _element.getDimension());
155 assert(perm.rows() == _element.getDimension() || perm.rows() == 1);
157 _element.getDimension(), _element.getDimension());
158 if (perm.rows() == _element.getDimension())
162 else if (perm.rows() == 1)
167 for (
unsigned ip = 0; ip < n_integration_points; ip++)
169 double pg_int_pt = 0.;
170 double pc_int_pt = 0.;
171 double T_int_pt = 0.0;
173 pc_int_pt, T_int_pt);
175 double const density_water =
176 two_phase_material_model.getLiquidDensity(pg_int_pt, T_int_pt);
178 double const Sw = two_phase_material_model.getSaturation(
179 material_id, t, pos, pg_int_pt, T_int_pt, pc_int_pt);
181 _saturation[ip] = Sw;
184 (pc_int_pt > two_phase_material_model.getCapillaryPressure(
185 material_id, t, pos, pg_int_pt, T_int_pt, 0.0))
187 : two_phase_material_model.getSaturationDerivative(
188 material_id, t, pos, pg_int_pt, T_int_pt, Sw);
190 double const p_vapor_nonwet =
191 _process_data.material->calculateVaporPressureNonwet(
192 pc_int_pt, T_int_pt, density_water);
194 double const p_gas_nonwet = pg_int_pt - p_vapor_nonwet;
196 double const x_gas_nonwet = p_gas_nonwet / pg_int_pt;
198 double const x_vapor_nonwet = p_vapor_nonwet / pg_int_pt;
200 double const X_gas_nonwet =
202 (x_gas_nonwet + x_vapor_nonwet * water_mol_mass / air_mol_mass);
204 double const mol_density_water = density_water / water_mol_mass;
207 double const d_p_vapor_nonwet_d_T =
208 _process_data.material->calculateDerivativedPgwdT(
209 pc_int_pt, T_int_pt, density_water);
210 double const d_p_vapor_nonwet_d_pc =
211 _process_data.material->calculateDerivativedPgwdPC(
212 pc_int_pt, T_int_pt, density_water);
213 double const d_mol_density_nonwet_d_T =
215 double const d_x_gas_nonwet_d_pg =
216 p_vapor_nonwet / pg_int_pt / pg_int_pt;
217 double const d_x_gas_nonwet_d_pc = -d_p_vapor_nonwet_d_pc / pg_int_pt;
218 double const d_x_gas_nonwet_d_T = -d_p_vapor_nonwet_d_T / pg_int_pt;
220 double const density_nonwet_gas =
222 double const density_nonwet_vapor =
224 double const density_nonwet = density_nonwet_gas + density_nonwet_vapor;
225 double const density_wet = density_water;
226 double const density_solid = _process_data.density_solid(t, pos)[0];
228 double const d_density_nonwet_d_T =
229 _process_data.material->calculatedDensityNonwetdT (
230 p_gas_nonwet, p_vapor_nonwet, pc_int_pt, T_int_pt, density_water);
232 _pressure_wetting[ip] = pg_int_pt - pc_int_pt;
234 double const heat_capacity_dry_gas =
235 _process_data.material->getSpecificHeatCapacityAir(pg_int_pt,
237 const double heat_capacity_water_vapor =
238 _process_data.material->getSpecificHeatCapacityVapor(pg_int_pt,
241 double const heat_capacity_water =
242 _process_data.material->getSpecificHeatCapacityWater(pg_int_pt,
244 double const heat_capacity_solid =
245 _process_data.material->getSpecificHeatCapacitySolid(pg_int_pt,
247 double const latent_heat_evaporation =
248 _process_data.latent_heat_evaporation(t, pos)[0];
250 double const enthalpy_nonwet_gas =
251 _process_data.material->getAirEnthalpySimple(
252 T_int_pt, heat_capacity_dry_gas, pg_int_pt);
254 double const enthalpy_wet =
255 _process_data.material->getLiquidWaterEnthalpySimple(
256 T_int_pt, heat_capacity_water, _pressure_wetting[ip]);
258 double const enthalpy_nonwet_vapor =
259 _process_data.material->getWaterVaporEnthalpySimple(
260 T_int_pt, heat_capacity_water_vapor, pg_int_pt,
261 latent_heat_evaporation);
262 double const enthalpy_nonwet =
263 enthalpy_nonwet_gas * X_gas_nonwet +
264 enthalpy_nonwet_vapor * (1 - X_gas_nonwet);
265 double const internal_energy_nonwet =
266 enthalpy_nonwet - pg_int_pt / density_nonwet;
267 double const internal_energy_wet = enthalpy_wet;
269 double const d_enthalpy_gas_nonwet_d_T =
271 double const d_enthalpy_nonwet_d_T =
272 heat_capacity_water * (1 - X_gas_nonwet) +
273 d_enthalpy_gas_nonwet_d_T * X_gas_nonwet;
276 double const porosity = two_phase_material_model.getPorosity(
277 material_id, t, pos, pg_int_pt, T_int_pt, 0);
279 Mgp.noalias() += porosity *
280 ((1 - Sw) * (mol_density_nonwet * d_x_gas_nonwet_d_pg +
281 x_gas_nonwet * d_mol_density_nonwet_d_pg)) *
282 _ip_data[ip].mass_operator;
283 Mgpc.noalias() += porosity *
284 ((1 - Sw) * mol_density_nonwet * d_x_gas_nonwet_d_pc -
285 mol_density_nonwet * x_gas_nonwet * dSwdpc) *
286 _ip_data[ip].mass_operator;
287 Mgt.noalias() += porosity *
288 ((1 - Sw) * (mol_density_nonwet * d_x_gas_nonwet_d_T +
289 x_gas_nonwet * d_mol_density_nonwet_d_T)) *
290 _ip_data[ip].mass_operator;
295 mol_density_nonwet * x_vapor_nonwet * (-dSwdpc) +
296 dSwdpc * mol_density_water) *
297 _ip_data[ip].mass_operator;
303 _ip_data[ip].mass_operator;
307 ((x_gas_nonwet * air_mol_mass + x_vapor_nonwet * water_mol_mass) *
308 d_mol_density_nonwet_d_pg * enthalpy_nonwet -
309 mol_density_nonwet * (water_mol_mass - air_mol_mass) *
310 d_x_gas_nonwet_d_pg * enthalpy_nonwet -
312 (1 - Sw) * _ip_data[ip].mass_operator;
314 porosity * (density_wet * internal_energy_wet -
315 density_nonwet * internal_energy_nonwet) *
316 dSwdpc * _ip_data[ip].mass_operator +
317 porosity * ((water_mol_mass - air_mol_mass) * enthalpy_nonwet /
319 (1 - Sw) * d_p_vapor_nonwet_d_pc * _ip_data[ip].mass_operator;
321 ((1 - porosity) * density_solid * heat_capacity_solid +
322 porosity * ((1 - Sw) * (d_density_nonwet_d_T * enthalpy_nonwet +
323 density_nonwet * d_enthalpy_nonwet_d_T) +
324 Sw * density_wet * heat_capacity_water)) *
325 _ip_data[ip].mass_operator;
328 double const k_rel_nonwet =
329 two_phase_material_model.getNonwetRelativePermeability(
330 t, pos, _pressure_wetting[ip], T_int_pt, Sw);
331 double const mu_nonwet = two_phase_material_model.getGasViscosity(
332 _pressure_wetting[ip], T_int_pt);
333 double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
334 double const diffusion_coeff_component_gas =
335 _process_data.diffusion_coeff_component_b(t, pos)[0];
338 double const k_rel_wet =
339 two_phase_material_model.getWetRelativePermeability(
340 t, pos, pg_int_pt, T_int_pt, Sw);
341 double const mu_wet =
342 two_phase_material_model.getLiquidViscosity(pg_int_pt, T_int_pt);
343 double const lambda_wet = k_rel_wet / mu_wet;
347 (_ip_data[ip].dNdx * pg_nodal_values);
350 (_ip_data[ip].dNdx * (pg_nodal_values - pc_nodal_values));
352 laplace_operator.noalias() = _ip_data[ip].dNdx.transpose() *
354 _ip_data[ip].integration_weight;
357 _ip_data[ip].integration_weight * _ip_data[ip].N.transpose() *
358 (d_density_nonwet_d_T * enthalpy_nonwet +
359 density_nonwet * d_enthalpy_nonwet_d_T) *
360 velocity_nonwet.transpose() * _ip_data[ip].dNdx +
361 _ip_data[ip].integration_weight * _ip_data[ip].N.transpose() *
362 heat_capacity_water * density_water * velocity_wet.transpose() *
365 double const heat_conductivity_dry_solid =
366 _process_data.material->getThermalConductivityDrySolid(pg_int_pt,
368 double const heat_conductivity_wet_solid =
369 _process_data.material->getThermalConductivityWetSolid(pg_int_pt,
371 double const heat_conductivity_unsaturated =
372 _process_data.material->calculateUnsatHeatConductivity(
373 t, pos, Sw, heat_conductivity_dry_solid,
374 heat_conductivity_wet_solid);
377 (mol_density_nonwet * x_gas_nonwet * lambda_nonwet) * laplace_operator +
378 ((1 - Sw) * porosity * diffusion_coeff_component_gas *
379 mol_density_nonwet * d_x_gas_nonwet_d_pg) *
380 _ip_data[ip].diffusion_operator;
381 Kgpc.noalias() += ((1 - Sw) * porosity * diffusion_coeff_component_gas *
382 mol_density_nonwet * d_x_gas_nonwet_d_pc) *
383 _ip_data[ip].diffusion_operator;
384 Kgt.noalias() += ((1 - Sw) * porosity * diffusion_coeff_component_gas *
385 mol_density_nonwet * d_x_gas_nonwet_d_T) *
386 _ip_data[ip].diffusion_operator;
388 Klp.noalias() += (mol_density_nonwet * x_vapor_nonwet * lambda_nonwet) *
390 mol_density_water * lambda_wet * laplace_operator -
391 ((1 - Sw) * porosity * diffusion_coeff_component_gas *
392 mol_density_nonwet * d_x_gas_nonwet_d_pg) *
393 _ip_data[ip].diffusion_operator;
394 Klpc.noalias() += (-mol_density_water * lambda_wet * laplace_operator) -
395 ((1 - Sw) * porosity * diffusion_coeff_component_gas *
396 mol_density_nonwet * d_x_gas_nonwet_d_pc) *
397 _ip_data[ip].diffusion_operator;
398 Klt.noalias() += -((1 - Sw) * porosity * diffusion_coeff_component_gas *
399 mol_density_nonwet * d_x_gas_nonwet_d_T) *
400 _ip_data[ip].diffusion_operator;
403 (lambda_nonwet * density_nonwet * enthalpy_nonwet +
404 lambda_wet * density_wet * enthalpy_wet) *
406 (1 - Sw) * porosity * diffusion_coeff_component_gas *
407 mol_density_nonwet * (air_mol_mass * enthalpy_nonwet_gas -
408 water_mol_mass * enthalpy_nonwet_vapor) *
409 d_x_gas_nonwet_d_pg * _ip_data[ip].diffusion_operator;
411 -lambda_wet * enthalpy_wet * density_wet * laplace_operator +
412 (1 - Sw) * porosity * diffusion_coeff_component_gas *
413 mol_density_nonwet * (air_mol_mass * enthalpy_nonwet_gas -
414 water_mol_mass * enthalpy_nonwet_vapor) *
415 d_x_gas_nonwet_d_pc * _ip_data[ip].diffusion_operator;
417 _ip_data[ip].dNdx.transpose() * heat_conductivity_unsaturated *
418 _ip_data[ip].dNdx * _ip_data[ip].integration_weight +
419 (1 - Sw) * porosity * diffusion_coeff_component_gas *
420 mol_density_nonwet * (air_mol_mass * enthalpy_nonwet_gas -
421 water_mol_mass * enthalpy_nonwet_vapor) *
422 d_x_gas_nonwet_d_T * _ip_data[ip].diffusion_operator;
424 if (_process_data.has_gravity)
426 auto const& b = _process_data.specific_body_force;
429 _ip_data[ip].integration_weight;
431 (mol_density_nonwet * x_gas_nonwet * lambda_nonwet * density_nonwet) *
434 (mol_density_water * lambda_wet * density_wet +
435 mol_density_nonwet * x_vapor_nonwet * lambda_nonwet * density_nonwet) *
438 (lambda_nonwet * density_nonwet * density_nonwet * enthalpy_nonwet +
439 lambda_wet * density_wet * density_wet * enthalpy_wet) *
443 if (_process_data.has_mass_lumping)
445 for (
unsigned row = 0; row < Mgp.cols(); row++)
447 for (
unsigned column = 0; column < Mgp.cols(); column++)
451 Mgp(row, row) += Mgp(row, column);
452 Mgp(row, column) = 0.0;
453 Mgpc(row, row) += Mgpc(row, column);
454 Mgpc(row, column) = 0.0;
455 Mgt(row, row) += Mgt(row, column);
456 Mgt(row, column) = 0.0;
457 Mlp(row, row) += Mlp(row, column);
458 Mlp(row, column) = 0.0;
459 Mlpc(row, row) += Mlpc(row, column);
460 Mlpc(row, column) = 0.0;
461 Mlt(row, row) += Mlt(row, column);
462 Mlt(row, column) = 0.0;
463 Mep(row, row) += Mep(row, column);
464 Mep(row, column) = 0.0;
465 Mepc(row, row) += Mepc(row, column);
466 Mepc(row, column) = 0.0;
467 Met(row, row) += Met(row, column);
468 Met(row, column) = 0.0;
Definition of the PiecewiseLinearInterpolation class.
std::optional< std::size_t > getElementID() const
void setElementID(std::size_t element_id)
typename ShapeMatricesType::NodalVectorType NodalVectorType
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
constexpr double Air
kg mol-1
constexpr double Water
kg mol-1
constexpr double IdealGasConstant
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
const unsigned NUM_NODAL_DOF