OGS
ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
Go to the documentation of this file.
1 
44 #pragma once
45 
47 
50 
52 
53 namespace ProcessLib
54 {
55 namespace ThermalTwoPhaseFlowWithPP
56 {
57 template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
58 void ThermalTwoPhaseFlowWithPPLocalAssembler<
59  ShapeFunction, IntegrationMethod,
60  GlobalDim>::assemble(double const t, double const /*dt*/,
61  std::vector<double> const& local_x,
62  std::vector<double> const& /*local_xdot*/,
63  std::vector<double>& local_M_data,
64  std::vector<double>& local_K_data,
65  std::vector<double>& local_b_data)
66 {
68  auto const& water_mol_mass =
70  auto const& air_mol_mass = MaterialLib::PhysicalConstant::MolarMass::Air;
71 
72  auto const local_matrix_size = local_x.size();
73 
74  assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
75 
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);
82 
83  auto Mgp =
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);
90 
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);
97 
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);
104 
105  NodalMatrixType laplace_operator =
106  NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
107 
108  auto Kgp =
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);
115 
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);
122 
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);
129 
130  auto Bg = local_b.template segment<nonwet_pressure_size>(
131  nonwet_pressure_matrix_index);
132  auto Bl =
133  local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
134  auto Be =
135  local_b.template segment<temperature_size>(temperature_matrix_index);
136 
137  unsigned const n_integration_points =
138  _integration_method.getNumberOfPoints();
139 
141  pos.setElementID(_element.getID());
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());
146 
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);
152 
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);
156  GlobalDimMatrixType permeability = GlobalDimMatrixType::Zero(
157  _element.getDimension(), _element.getDimension());
158  if (perm.rows() == _element.getDimension())
159  {
160  permeability = perm;
161  }
162  else if (perm.rows() == 1)
163  {
164  permeability.diagonal().setConstant(perm(0, 0));
165  }
166 
167  for (unsigned ip = 0; ip < n_integration_points; ip++)
168  {
169  double pg_int_pt = 0.;
170  double pc_int_pt = 0.;
171  double T_int_pt = 0.0;
172  NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, pg_int_pt,
173  pc_int_pt, T_int_pt);
174 
175  double const density_water =
176  two_phase_material_model.getLiquidDensity(pg_int_pt, T_int_pt);
177 
178  double const Sw = two_phase_material_model.getSaturation(
179  material_id, t, pos, pg_int_pt, T_int_pt, pc_int_pt);
180 
181  _saturation[ip] = Sw;
182 
183  double dSwdpc =
184  (pc_int_pt > two_phase_material_model.getCapillaryPressure(
185  material_id, t, pos, pg_int_pt, T_int_pt, 0.0))
186  ? 0.0
187  : two_phase_material_model.getSaturationDerivative(
188  material_id, t, pos, pg_int_pt, T_int_pt, Sw);
189 
190  double const p_vapor_nonwet =
191  _process_data.material->calculateVaporPressureNonwet(
192  pc_int_pt, T_int_pt, density_water);
193  // partial pressure of gas component
194  double const p_gas_nonwet = pg_int_pt - p_vapor_nonwet;
195  // molar fraction of gas component in nonwet phase
196  double const x_gas_nonwet = p_gas_nonwet / pg_int_pt;
197  // molar fraction of water vapor in nonwet phase
198  double const x_vapor_nonwet = p_vapor_nonwet / pg_int_pt;
199  // mass fraction of gas component in the nonwet phase
200  double const X_gas_nonwet =
201  x_gas_nonwet /
202  (x_gas_nonwet + x_vapor_nonwet * water_mol_mass / air_mol_mass);
203  double const mol_density_nonwet = pg_int_pt / IdealGasConstant / T_int_pt;
204  double const mol_density_water = density_water / water_mol_mass;
205 
206  double const d_mol_density_nonwet_d_pg = 1 / IdealGasConstant / T_int_pt;
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 =
214  -pg_int_pt / IdealGasConstant / T_int_pt / T_int_pt;
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;
219 
220  double const density_nonwet_gas =
221  p_gas_nonwet * air_mol_mass / IdealGasConstant / T_int_pt;
222  double const density_nonwet_vapor =
223  p_vapor_nonwet * water_mol_mass / IdealGasConstant / T_int_pt;
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];
227  // Derivative of nonwet phase density in terms of T
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);
231 
232  _pressure_wetting[ip] = pg_int_pt - pc_int_pt;
233  // heat capacity of nonwet phase
234  double const heat_capacity_dry_gas =
235  _process_data.material->getSpecificHeatCapacityAir(pg_int_pt,
236  T_int_pt);
237  const double heat_capacity_water_vapor =
238  _process_data.material->getSpecificHeatCapacityVapor(pg_int_pt,
239  T_int_pt);
240 
241  double const heat_capacity_water =
242  _process_data.material->getSpecificHeatCapacityWater(pg_int_pt,
243  T_int_pt);
244  double const heat_capacity_solid =
245  _process_data.material->getSpecificHeatCapacitySolid(pg_int_pt,
246  T_int_pt);
247  double const latent_heat_evaporation =
248  _process_data.latent_heat_evaporation(t, pos)[0];
249 
250  double const enthalpy_nonwet_gas =
251  _process_data.material->getAirEnthalpySimple(
252  T_int_pt, heat_capacity_dry_gas, pg_int_pt);
253 
254  double const enthalpy_wet =
255  _process_data.material->getLiquidWaterEnthalpySimple(
256  T_int_pt, heat_capacity_water, _pressure_wetting[ip]);
257 
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 =
270  heat_capacity_dry_gas + IdealGasConstant / air_mol_mass;
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;
274  // Assemble M matrix
275  // nonwetting
276  double const porosity = two_phase_material_model.getPorosity(
277  material_id, t, pos, pg_int_pt, T_int_pt, 0);
278 
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;
291 
292  Mlpc.noalias() +=
293  porosity *
294  ((1 - Sw) * d_p_vapor_nonwet_d_pc / IdealGasConstant / T_int_pt +
295  mol_density_nonwet * x_vapor_nonwet * (-dSwdpc) +
296  dSwdpc * mol_density_water) *
297  _ip_data[ip].mass_operator;
298  Mlt.noalias() +=
299  porosity *
300  ((1 - Sw) *
301  (d_p_vapor_nonwet_d_T / IdealGasConstant / T_int_pt -
302  p_vapor_nonwet / IdealGasConstant / T_int_pt / T_int_pt)) *
303  _ip_data[ip].mass_operator;
304 
305  Mep.noalias() +=
306  porosity *
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 -
311  1) *
312  (1 - Sw) * _ip_data[ip].mass_operator;
313  Mepc.noalias() +=
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 /
318  IdealGasConstant / T_int_pt) *
319  (1 - Sw) * d_p_vapor_nonwet_d_pc * _ip_data[ip].mass_operator;
320  Met.noalias() +=
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;
326 
327  // nonwet
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];
336 
337  // wet
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;
344 
345  GlobalDimVectorType const velocity_nonwet =
346  -lambda_nonwet * permeability *
347  (_ip_data[ip].dNdx * pg_nodal_values);
348  GlobalDimVectorType const velocity_wet =
349  -lambda_wet * permeability *
350  (_ip_data[ip].dNdx * (pg_nodal_values - pc_nodal_values));
351 
352  laplace_operator.noalias() = _ip_data[ip].dNdx.transpose() *
353  permeability * _ip_data[ip].dNdx *
354  _ip_data[ip].integration_weight;
355 
356  Ket.noalias() +=
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() *
363  _ip_data[ip].dNdx;
364 
365  double const heat_conductivity_dry_solid =
366  _process_data.material->getThermalConductivityDrySolid(pg_int_pt,
367  T_int_pt);
368  double const heat_conductivity_wet_solid =
369  _process_data.material->getThermalConductivityWetSolid(pg_int_pt,
370  T_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);
375  // Laplace
376  Kgp.noalias() +=
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;
387 
388  Klp.noalias() += (mol_density_nonwet * x_vapor_nonwet * lambda_nonwet) *
389  laplace_operator +
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;
401 
402  Kep.noalias() +=
403  (lambda_nonwet * density_nonwet * enthalpy_nonwet +
404  lambda_wet * density_wet * enthalpy_wet) *
405  laplace_operator +
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;
410  Kepc.noalias() +=
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;
416  Ket.noalias() +=
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;
423 
424  if (_process_data.has_gravity)
425  {
426  auto const& b = _process_data.specific_body_force;
427  NodalVectorType gravity_operator = _ip_data[ip].dNdx.transpose() *
428  permeability * b *
429  _ip_data[ip].integration_weight;
430  Bg.noalias() +=
431  (mol_density_nonwet * x_gas_nonwet * lambda_nonwet * density_nonwet) *
432  gravity_operator;
433  Bl.noalias() +=
434  (mol_density_water * lambda_wet * density_wet +
435  mol_density_nonwet * x_vapor_nonwet * lambda_nonwet * density_nonwet) *
436  gravity_operator;
437  Be.noalias() +=
438  (lambda_nonwet * density_nonwet * density_nonwet * enthalpy_nonwet +
439  lambda_wet * density_wet * density_wet * enthalpy_wet) *
440  gravity_operator;
441  } // end of has gravity
442  }
443  if (_process_data.has_mass_lumping)
444  {
445  for (unsigned row = 0; row < Mgp.cols(); row++)
446  {
447  for (unsigned column = 0; column < Mgp.cols(); column++)
448  {
449  if (row != column)
450  {
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;
469  }
470  }
471  }
472  } // end of mass-lumping
473 }
474 
475 } // namespace ThermalTwoPhaseFlowWithPP
476 } // namespace ProcessLib
Definition of the PiecewiseLinearInterpolation class.
std::optional< std::size_t > getElementID() const
void setElementID(std::size_t element_id)
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:79