OGS
ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim > Class Template Reference

Detailed Description

template<typename ShapeFunction, int GlobalDim>
class ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >

Definition at line 87 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

#include <ThermalTwoPhaseFlowWithPPLocalAssembler.h>

Inheritance diagram for ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >:
[legend]
Collaboration diagram for ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >:
[legend]

Public Member Functions

 ThermalTwoPhaseFlowWithPPLocalAssembler (MeshLib::Element const &element, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, ThermalTwoPhaseFlowWithPPProcessData const &process_data)
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
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
std::vector< double > const & getIntPtSaturation (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
std::vector< double > const & getIntPtWettingPressure (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
std::vector< double > const & getIntPtLiquidMolFracContaminant (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
std::vector< double > const & getIntPtGasMolFracWater (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
std::vector< double > const & getIntPtGasMolFracContaminant (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
virtual void setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void preAssemble (double const, double const, std::vector< double > const &)
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
 Fits to staggered scheme.
virtual int getNumberOfVectorElementsForDeformation () const
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual ~ExtrapolatableElement ()=default

Private Types

using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices
using LocalAssemblerTraits
using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType
using GlobalDimNodalMatrixType
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType
using NodalVectorType = typename ShapeMatricesType::NodalVectorType
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
using LocalMatrixType = typename LocalAssemblerTraits::LocalMatrix
using LocalVectorType = typename LocalAssemblerTraits::LocalVector

Private Attributes

MeshLib::Element const & _element
NumLib::GenericIntegrationMethod const & _integration_method
ThermalTwoPhaseFlowWithPPProcessData const & _process_data
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType > > > _ip_data
std::vector< double > _saturation
std::vector< double > _pressure_wetting
std::vector< double > _liquid_molar_fraction_contaminant
std::vector< double > _gas_molar_fraction_water
std::vector< double > _gas_molar_fraction_contaminant

Static Private Attributes

static const int nonwet_pressure_matrix_index = 0
static const int cap_pressure_matrix_index = ShapeFunction::NPOINTS
static const int contaminant_matrix_index = 2 * ShapeFunction::NPOINTS
static const int temperature_matrix_index = 3 * ShapeFunction::NPOINTS
static const int nonwet_pressure_size = ShapeFunction::NPOINTS
static const int cap_pressure_size = ShapeFunction::NPOINTS
static const int contaminant_size = ShapeFunction::NPOINTS
static const int temperature_size = ShapeFunction::NPOINTS

Member Typedef Documentation

◆ GlobalDimMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType
private

Definition at line 101 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ GlobalDimNodalMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::GlobalDimNodalMatrixType
private
Initial value:
typename ShapeMatricesType::GlobalDimNodalMatrixType
EigenFixedShapeMatrixPolicy::GlobalDimNodalMatrixType
MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType
Definition ShapeMatrixPolicy.h:75

Definition at line 97 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ GlobalDimVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
private

Definition at line 102 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ LocalAssemblerTraits

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::LocalAssemblerTraits
private
Initial value:
ProcessLib::LocalAssemblerTraits<
ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>
ProcessLib::VolumetricSourceTermLocalAssembler::ShapeMatricesType
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
Definition VolumetricSourceTermFEM.h:35
ProcessLib::LocalAssemblerTraits
detail::LocalAssemblerTraitsFixed< ShpPol, NNodes, NodalDOF, Dim > LocalAssemblerTraits
Definition LocalAssemblerTraits.h:178

Definition at line 93 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ LocalMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::LocalMatrixType = typename LocalAssemblerTraits::LocalMatrix
private

Definition at line 103 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ LocalVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::LocalVectorType = typename LocalAssemblerTraits::LocalVector
private

Definition at line 104 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ NodalMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::NodalMatrixType = typename ShapeMatricesType::NodalMatrixType
private

Definition at line 99 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ NodalRowVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType
private

Definition at line 95 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ NodalVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::NodalVectorType = typename ShapeMatricesType::NodalVectorType
private

Definition at line 100 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ ShapeMatrices

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::ShapeMatrices = typename ShapeMatricesType::ShapeMatrices
private

Definition at line 91 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ ShapeMatricesType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>
private

Definition at line 90 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Constructor & Destructor Documentation

◆ ThermalTwoPhaseFlowWithPPLocalAssembler()

template<typename ShapeFunction, int GlobalDim>
ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::ThermalTwoPhaseFlowWithPPLocalAssembler ( MeshLib::Element const & element,
std::size_t const ,
NumLib::GenericIntegrationMethod const & integration_method,
bool const is_axially_symmetric,
ThermalTwoPhaseFlowWithPPProcessData const & process_data )
inline

Definition at line 107 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

113 : _element(element),
114 _integration_method(integration_method),
115 _process_data(process_data),
116 _saturation(
117 std::vector<double>(_integration_method.getNumberOfPoints())),
118 _pressure_wetting(
119 std::vector<double>(_integration_method.getNumberOfPoints())),
120 _liquid_molar_fraction_contaminant(
121 std::vector<double>(_integration_method.getNumberOfPoints())),
122 _gas_molar_fraction_water(
123 std::vector<double>(_integration_method.getNumberOfPoints())),
124 _gas_molar_fraction_contaminant(
125 std::vector<double>(_integration_method.getNumberOfPoints()))
126 {
127 unsigned const n_integration_points =
128 _integration_method.getNumberOfPoints();
129 _ip_data.reserve(n_integration_points);
130 auto const shape_matrices =
131 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
132 GlobalDim>(element, is_axially_symmetric,
133 _integration_method);
134 for (unsigned ip = 0; ip < n_integration_points; ip++)
135 {
136 auto const& sm = shape_matrices[ip];
137 const double integration_factor = sm.integralMeasure * sm.detJ;
138 _ip_data.emplace_back(
139 sm.N, sm.dNdx,
140 sm.integralMeasure * sm.detJ *
141 _integration_method.getWeightedPoint(ip).getWeight(),
142 sm.N.transpose() * sm.N * integration_factor *
143 _integration_method.getWeightedPoint(ip).getWeight(),
144 sm.dNdx.transpose() * sm.dNdx * integration_factor *
145 _integration_method.getWeightedPoint(ip).getWeight());
146 }
147 }
ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler::_ip_data
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType > > > _ip_data
Definition ThermalTwoPhaseFlowWithPPLocalAssembler.h:226

References _element, _gas_molar_fraction_contaminant, _gas_molar_fraction_water, _integration_method, _ip_data, _liquid_molar_fraction_contaminant, _pressure_wetting, _process_data, _saturation, and NumLib::initShapeMatrices().

Member Function Documentation

◆ assemble()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::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 )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 20 of file ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h.

27{
28 using MaterialLib::PhysicalConstant::CelsiusZeroInKelvin;
29 using MaterialLib::PhysicalConstant::IdealGasConstant;
30
31 auto const local_matrix_size = local_x.size();
32
33 assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
34
35 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
36 local_M_data, local_matrix_size, local_matrix_size);
37 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
38 local_K_data, local_matrix_size, local_matrix_size);
39 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
40 local_b_data, local_matrix_size);
41
42 auto Map =
43 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
44 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
45 auto Mapc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
46 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
47 auto Max = local_M.template block<nonwet_pressure_size, contaminant_size>(
48 nonwet_pressure_matrix_index, contaminant_matrix_index);
49 auto Mat = local_M.template block<nonwet_pressure_size, temperature_size>(
50 nonwet_pressure_matrix_index, temperature_matrix_index);
51
52 auto Mwp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
53 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
54 auto Mwpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
55 cap_pressure_matrix_index, cap_pressure_matrix_index);
56 auto Mwx = local_M.template block<cap_pressure_size, contaminant_size>(
57 cap_pressure_matrix_index, contaminant_matrix_index);
58 auto Mwt = local_M.template block<cap_pressure_size, temperature_size>(
59 cap_pressure_matrix_index, temperature_matrix_index);
60
61 auto Mcp = local_M.template block<contaminant_size, nonwet_pressure_size>(
62 contaminant_matrix_index, nonwet_pressure_matrix_index);
63 auto Mcpc = local_M.template block<contaminant_size, cap_pressure_size>(
64 contaminant_matrix_index, cap_pressure_matrix_index);
65 auto Mcx = local_M.template block<contaminant_size, contaminant_size>(
66 contaminant_matrix_index, contaminant_matrix_index);
67 auto Mct = local_M.template block<contaminant_size, temperature_size>(
68 contaminant_matrix_index, temperature_matrix_index);
69
70 auto Mep = local_M.template block<temperature_size, nonwet_pressure_size>(
71 temperature_matrix_index, nonwet_pressure_matrix_index);
72 auto Mepc = local_M.template block<temperature_size, cap_pressure_size>(
73 temperature_matrix_index, cap_pressure_matrix_index);
74 auto Mex = local_M.template block<temperature_size, contaminant_size>(
75 temperature_matrix_index, contaminant_matrix_index);
76 auto Met = local_M.template block<temperature_size, temperature_size>(
77 temperature_matrix_index, temperature_matrix_index);
78
79 NodalMatrixType laplace_operator =
80 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
81
82 auto Kap =
83 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
84 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
85 auto Kapc = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
86 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
87 auto Kax = local_K.template block<nonwet_pressure_size, contaminant_size>(
88 nonwet_pressure_matrix_index, contaminant_matrix_index);
89 auto Kat = local_K.template block<nonwet_pressure_size, temperature_size>(
90 nonwet_pressure_matrix_index, temperature_matrix_index);
91
92 auto Kwp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
93 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
94 auto Kwpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
95 cap_pressure_matrix_index, cap_pressure_matrix_index);
96 auto Kwx = local_K.template block<cap_pressure_size, contaminant_size>(
97 cap_pressure_matrix_index, contaminant_matrix_index);
98 auto Kwt = local_K.template block<cap_pressure_size, temperature_size>(
99 cap_pressure_matrix_index, temperature_matrix_index);
100
101 auto Kcp = local_K.template block<contaminant_size, nonwet_pressure_size>(
102 contaminant_matrix_index, nonwet_pressure_matrix_index);
103 auto Kcpc = local_K.template block<contaminant_size, cap_pressure_size>(
104 contaminant_matrix_index, cap_pressure_matrix_index);
105 auto Kcx = local_K.template block<contaminant_size, contaminant_size>(
106 contaminant_matrix_index, contaminant_matrix_index);
107 auto Kct = local_K.template block<contaminant_size, temperature_size>(
108 contaminant_matrix_index, temperature_matrix_index);
109
110 auto Kep = local_K.template block<temperature_size, nonwet_pressure_size>(
111 temperature_matrix_index, nonwet_pressure_matrix_index);
112 auto Kepc = local_K.template block<temperature_size, cap_pressure_size>(
113 temperature_matrix_index, cap_pressure_matrix_index);
114 auto Ket = local_K.template block<temperature_size, temperature_size>(
115 temperature_matrix_index, temperature_matrix_index);
116
117 auto Ba = local_b.template segment<nonwet_pressure_size>(
118 nonwet_pressure_matrix_index);
119 auto Bw =
120 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
121 auto Bc =
122 local_b.template segment<contaminant_size>(contaminant_matrix_index);
123 auto Be =
124 local_b.template segment<temperature_size>(temperature_matrix_index);
125
126 unsigned const n_integration_points =
127 _integration_method.getNumberOfPoints();
128
129 auto const num_nodes = ShapeFunction::NPOINTS;
130 auto const pg_nodal_values =
131 Eigen::Map<const NodalVectorType>(&local_x[0], num_nodes);
132 auto const pc_nodal_values =
133 Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
134
135 MaterialPropertyLib::VariableArray vars;
136
137 GlobalDimMatrixType const& I(
138 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
139
140 for (unsigned ip = 0; ip < n_integration_points; ip++)
141 {
142 auto const& ip_data = _ip_data[ip];
143 auto const& N = ip_data.N;
144 auto const& dNdx = ip_data.dNdx;
145 auto const& w = ip_data.integration_weight;
146 auto const& mass_operator = ip_data.mass_operator;
147 auto const& diffusion_operator = ip_data.diffusion_operator;
148 double pg_int_pt = 0.;
149 double pc_int_pt = 0.;
150 double Xc_int_pt = 0.;
151 double T_int_pt = 0.;
152 NumLib::shapeFunctionInterpolate(local_x, N, pg_int_pt, pc_int_pt,
153 Xc_int_pt, T_int_pt);
154
155 ParameterLib::SpatialPosition const pos{
156 std::nullopt, _element.getID(),
157 MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
158 ShapeMatricesType>(
159 _element, N))};
160
161 _pressure_wetting[ip] = pg_int_pt - pc_int_pt;
162 double const ideal_gas_constant_times_T_int_pt =
163 IdealGasConstant * T_int_pt;
164 vars.temperature = T_int_pt;
165 vars.capillary_pressure = pc_int_pt;
166 vars.gas_phase_pressure = pg_int_pt;
167
168 auto const& medium =
169 *_process_data.media_map.getMedium(this->_element.getID());
170 auto const& liquid_phase = medium.phase("AqueousLiquid");
171 auto const& solid_phase = medium.phase("Solid");
172 auto const& gas_phase = medium.phase("Gas");
173
174 auto const& water_vapour_component = gas_phase.component("w");
175 auto const& dry_air_component = gas_phase.component("a");
176 auto const& contaminant_vapour_component = gas_phase.component("c");
177 auto const& dissolved_contaminant_component =
178 liquid_phase.component("c");
179
180 auto const porosity =
181 medium.property(MaterialPropertyLib::PropertyType::porosity)
182 .template value<double>(vars, pos, t, dt);
183
184 auto const water_mol_mass =
185 water_vapour_component
186 .property(MaterialPropertyLib::PropertyType::molar_mass)
187 .template value<double>(vars, pos, t, dt);
188 auto const air_mol_mass =
189 dry_air_component
190 .property(MaterialPropertyLib::PropertyType::molar_mass)
191 .template value<double>(vars, pos, t, dt);
192 auto const contaminant_mol_mass =
193 contaminant_vapour_component
194 .property(MaterialPropertyLib::PropertyType::molar_mass)
195 .template value<double>(vars, pos, t, dt);
196
197 double const Sw =
198 medium.property(MaterialPropertyLib::PropertyType::saturation)
199 .template value<double>(vars, pos, t, dt);
200
201 _saturation[ip] = Sw;
202 vars.liquid_saturation = Sw;
203
204 double const dSw_dpc =
205 medium.property(MaterialPropertyLib::PropertyType::saturation)
206 .template dValue<double>(
207 vars, MaterialPropertyLib::Variable::capillary_pressure,
208 pos, t, dt);
209
210 auto const density_water =
211 liquid_phase.property(MaterialPropertyLib::PropertyType::density)
212 .template value<double>(vars, pos, t, dt);
213
214 // molar densities and derivatives
215 double const mol_density_water = density_water / water_mol_mass;
216 double const mol_density_wet = mol_density_water;
217
218 double const mol_density_nonwet =
219 pg_int_pt / ideal_gas_constant_times_T_int_pt;
220 double const mol_density_tot =
221 Sw * mol_density_wet + (1 - Sw) * mol_density_nonwet;
222
223 double const d_mol_density_nonwet_dpg =
224 1 / ideal_gas_constant_times_T_int_pt;
225 double const d_mol_density_nonwet_dT = -mol_density_nonwet / T_int_pt;
226 double const d_mol_density_tot_dpc =
227 (mol_density_wet - mol_density_nonwet) * dSw_dpc;
228 double const d_mol_density_tot_dpg =
229 (1 - Sw) * d_mol_density_nonwet_dpg;
230 double const d_mol_density_tot_dT = (1 - Sw) * d_mol_density_nonwet_dT;
231
232 // specific latent heat of evaporation
233 double const latent_heat_evaporation =
234 water_vapour_component
235 .property(
236 MaterialPropertyLib::PropertyType::specific_latent_heat)
237 .template value<double>(vars, pos, t, dt);
238
239 vars.enthalpy_of_evaporation = latent_heat_evaporation;
240
241 // saturated vapour pressure
242 vars.molar_mass = air_mol_mass;
243 double const p_sat =
244 water_vapour_component
245 .property(MaterialPropertyLib::PropertyType::vapour_pressure)
246 .template value<double>(vars, pos, t, dt);
247 double const dp_sat_dT =
248 water_vapour_component
249 .property(MaterialPropertyLib::PropertyType::vapour_pressure)
250 .template dValue<double>(
251 vars, MaterialPropertyLib::Variable::temperature, pos, t,
252 dt);
253
254 // Kelvin-Laplace correction for menisci
255 double const K = std::exp(-pc_int_pt / mol_density_water /
256 ideal_gas_constant_times_T_int_pt);
257 double const dK_dT = pc_int_pt / mol_density_water /
258 ideal_gas_constant_times_T_int_pt / T_int_pt * K;
259
260 // vapour pressure inside pore space (water partial pressure in gas
261 // phase)
262 double const p_vapour_nonwet = p_sat * K;
263 double const d_p_vapour_nonwet_dT = dp_sat_dT * K + p_sat * dK_dT;
264 double const d_p_vapour_nonwet_dpc =
265 p_vapour_nonwet *
266 (-1 / mol_density_water / ideal_gas_constant_times_T_int_pt);
267
268 // Henry constant of organic contaminant
269 double const henry_contam =
270 contaminant_vapour_component
271 .property(MaterialPropertyLib::PropertyType::henry_coefficient)
272 .template value<double>(vars, pos, t, dt);
273 double d_henry_contaminant_dT =
274 contaminant_vapour_component
275 .property(MaterialPropertyLib::PropertyType::henry_coefficient)
276 .template dValue<double>(
277 vars, MaterialPropertyLib::Variable::temperature, pos, t,
278 dt);
279
280 // mass distribution coefficients of contam. and water
281 double const k_c = pg_int_pt * henry_contam / mol_density_wet;
282 double const k_w = pg_int_pt / p_vapour_nonwet;
283
284 // intermediate parameter
285 double const Ntot_c =
286 Sw * mol_density_wet * k_c + (1 - Sw) * mol_density_nonwet;
287 // phase-wise component molar fractions
288 double const x_contaminant_nonwet =
289 Xc_int_pt * mol_density_tot / Ntot_c;
290 double const x_contaminant_wet = k_c * x_contaminant_nonwet;
291 double const x_water_wet = 1 - x_contaminant_wet;
292 double const x_water_nonwet = x_water_wet / k_w;
293 double const x_air_nonwet = 1 - x_water_nonwet - x_contaminant_nonwet;
294
295 _gas_molar_fraction_water[ip] = x_water_nonwet;
296 _liquid_molar_fraction_contaminant[ip] = x_contaminant_wet;
297 _gas_molar_fraction_contaminant[ip] = x_contaminant_nonwet;
298
299 double const d_kc_dpg = henry_contam / mol_density_wet;
300 double const d_kc_dT =
301 pg_int_pt * d_henry_contaminant_dT / mol_density_wet;
302
303 // derivatives of component molar fractions w.r.t. PVs
304 double const d_x_contaminant_nonwet_dpc =
305 Xc_int_pt * (d_mol_density_tot_dpc / Ntot_c -
306 (mol_density_wet * k_c - mol_density_nonwet) *
307 dSw_dpc * mol_density_tot / Ntot_c / Ntot_c);
308 double const d_x_contaminant_nonwet_dpg =
309 Xc_int_pt * (d_mol_density_tot_dpg / Ntot_c -
310 (Sw * mol_density_wet * d_kc_dpg +
311 (1 - Sw) * d_mol_density_nonwet_dpg) *
312 mol_density_tot / Ntot_c / Ntot_c);
313 double const d_x_contaminant_nonwet_dXc = mol_density_tot / Ntot_c;
314 double const d_x_contaminant_nonwet_dT =
315 Xc_int_pt * (d_mol_density_tot_dT / Ntot_c -
316 (Sw * mol_density_wet * d_kc_dT +
317 (1 - Sw) * d_mol_density_nonwet_dT) *
318 mol_density_tot / Ntot_c / Ntot_c);
319
320 double const d_x_contaminant_wet_dpc = k_c * d_x_contaminant_nonwet_dpc;
321 double const d_x_contaminant_wet_dpg =
322 k_c * d_x_contaminant_nonwet_dpg + d_kc_dpg * x_contaminant_nonwet;
323 double const d_x_contaminant_wet_dXc = k_c * d_x_contaminant_nonwet_dXc;
324 double const d_x_contaminant_wet_dT =
325 k_c * d_x_contaminant_nonwet_dT + d_kc_dT * x_contaminant_nonwet;
326
327 double const d_x_water_wet_dpc = -d_x_contaminant_wet_dpc;
328 double const d_x_water_wet_dpg = -d_x_contaminant_wet_dpg;
329 double const d_x_water_wet_dXc = -d_x_contaminant_wet_dXc;
330 double const d_x_water_wet_dT = -d_x_contaminant_wet_dT;
331
332 double const d_x_water_nonwet_dpc =
333 (d_p_vapour_nonwet_dpc * x_water_wet +
334 p_vapour_nonwet * d_x_water_wet_dpc) /
335 pg_int_pt;
336 double const d_x_water_nonwet_dpg =
337 p_vapour_nonwet * (d_x_water_wet_dpg / pg_int_pt -
338 x_water_wet / pg_int_pt / pg_int_pt);
339 double const d_x_water_nonwet_dXc = d_x_water_wet_dXc / k_w;
340 double const d_x_water_nonwet_dT =
341 (d_p_vapour_nonwet_dT * x_water_wet +
342 p_vapour_nonwet * d_x_water_wet_dT) /
343 pg_int_pt;
344
345 double const d_x_air_nonwet_dpc =
346 -d_x_water_nonwet_dpc - d_x_contaminant_nonwet_dpc;
347 double const d_x_air_nonwet_dpg =
348 -d_x_water_nonwet_dpg - d_x_contaminant_nonwet_dpg;
349 double const d_x_air_nonwet_dXc =
350 -d_x_water_nonwet_dXc - d_x_contaminant_nonwet_dXc;
351 double const d_x_air_nonwet_dT =
352 -d_x_water_nonwet_dT - d_x_contaminant_nonwet_dT;
353
354 // mass densities
355 double const density_contaminant_wet =
356 mol_density_wet * contaminant_mol_mass * x_contaminant_wet;
357 double const density_wet = density_water + density_contaminant_wet;
358
359 double const density_air_nonwet =
360 mol_density_nonwet * air_mol_mass * x_air_nonwet;
361 double const density_water_nonwet =
362 mol_density_nonwet * water_mol_mass * x_water_nonwet;
363 double const density_contaminant_nonwet =
364 mol_density_nonwet * contaminant_mol_mass * x_contaminant_nonwet;
365 double const density_nonwet = density_air_nonwet +
366 density_water_nonwet +
367 density_contaminant_nonwet;
368
369 auto const density_solid =
370 solid_phase.property(MaterialPropertyLib::PropertyType::density)
371 .template value<double>(vars, pos, t, dt);
372
373 // derivatives of nonwet phase densities w.r.t. PVs
374 double const d_density_nonwet_dpg =
375 d_mol_density_nonwet_dpg *
376 (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
377 contaminant_mol_mass * x_contaminant_nonwet) +
378 mol_density_nonwet *
379 (air_mol_mass * d_x_air_nonwet_dpg +
380 water_mol_mass * d_x_water_nonwet_dpg +
381 contaminant_mol_mass * d_x_contaminant_nonwet_dpg);
382 double const d_density_nonwet_dpc =
383 mol_density_nonwet *
384 (air_mol_mass * d_x_air_nonwet_dpc +
385 water_mol_mass * d_x_water_nonwet_dpc +
386 contaminant_mol_mass * d_x_contaminant_nonwet_dpc);
387 double const d_density_nonwet_dXc =
388 mol_density_nonwet *
389 (air_mol_mass * d_x_air_nonwet_dXc +
390 water_mol_mass * d_x_water_nonwet_dXc +
391 contaminant_mol_mass * d_x_contaminant_nonwet_dXc);
392 double const d_density_nonwet_dT =
393 d_mol_density_nonwet_dT *
394 (air_mol_mass * x_air_nonwet + water_mol_mass * x_water_nonwet +
395 contaminant_mol_mass * x_contaminant_nonwet) +
396 mol_density_nonwet *
397 (air_mol_mass * d_x_air_nonwet_dT +
398 water_mol_mass * d_x_water_nonwet_dT +
399 contaminant_mol_mass * d_x_contaminant_nonwet_dT);
400
401 double const mol_mass_nonwet =
402 x_water_nonwet * water_mol_mass + x_air_nonwet * air_mol_mass +
403 x_contaminant_nonwet * contaminant_mol_mass;
404 // phase-wise component mass fractions
405 double const X_water_nonwet =
406 x_water_nonwet * water_mol_mass / mol_mass_nonwet;
407 double const X_air_nonwet =
408 x_air_nonwet * air_mol_mass / mol_mass_nonwet;
409 double const X_contaminant_nonwet = 1 - X_water_nonwet - X_air_nonwet;
410
411 // spec. heat capacities
412 double const heat_capacity_dry_air =
413 dry_air_component
414 .property(
415 MaterialPropertyLib::PropertyType::specific_heat_capacity)
416 .template value<double>(vars, pos, t, dt);
417 double const heat_capacity_water_vapour =
418 water_vapour_component
419 .property(
420 MaterialPropertyLib::PropertyType::specific_heat_capacity)
421 .template value<double>(vars, pos, t, dt);
422 double const heat_capacity_contaminant_vapour =
423 contaminant_vapour_component
424 .property(
425 MaterialPropertyLib::PropertyType::specific_heat_capacity)
426 .template value<double>(vars, pos, t, dt);
427
428 double const heat_capacity_water =
429 liquid_phase
430 .property(
431 MaterialPropertyLib::PropertyType::specific_heat_capacity)
432 .template value<double>(vars, pos, t, dt);
433 double const heat_capacity_solid =
434 solid_phase
435 .property(
436 MaterialPropertyLib::PropertyType::specific_heat_capacity)
437 .template value<double>(vars, pos, t, dt);
438
439 // enthalpies of gaseous components
440 // Note: h^a_G = C^a_P * (T - T0) = C^a_V * (T - T0) + RT/M^a, thus
441 // "C^a_V" should be used in the following. Same for contaminant.
442 double const enthalpy_air_nonwet =
443 heat_capacity_dry_air * (T_int_pt - CelsiusZeroInKelvin) +
444 IdealGasConstant * T_int_pt / air_mol_mass;
445 double const enthalpy_water_nonwet =
446 heat_capacity_water_vapour * (T_int_pt - CelsiusZeroInKelvin) +
447 latent_heat_evaporation;
448 double const enthalpy_contaminant_nonwet =
449 heat_capacity_contaminant_vapour *
450 (T_int_pt - CelsiusZeroInKelvin) +
451 IdealGasConstant * T_int_pt / contaminant_mol_mass;
452
453 // gas and liquid phase enthalpies
454 double const enthalpy_nonwet =
455 enthalpy_air_nonwet * X_air_nonwet +
456 enthalpy_water_nonwet * X_water_nonwet +
457 enthalpy_contaminant_nonwet * X_contaminant_nonwet;
458 double const enthalpy_wet =
459 heat_capacity_water * (T_int_pt - CelsiusZeroInKelvin);
460
461 // gas and liquid phase internal energies
462 double const internal_energy_nonwet =
463 enthalpy_nonwet - pg_int_pt / density_nonwet;
464 double const internal_energy_wet = enthalpy_wet;
465
466 // derivatives of enthalpies w.r.t. temperature
467 double const d_enthalpy_air_nonwet_dT =
468 heat_capacity_dry_air + IdealGasConstant / air_mol_mass;
469 double const d_enthalpy_contaminant_nonwet_dT =
470 heat_capacity_contaminant_vapour +
471 IdealGasConstant / contaminant_mol_mass;
472
473 double const d_enthalpy_nonwet_dT =
474 heat_capacity_water * X_water_nonwet +
475 d_enthalpy_air_nonwet_dT * X_air_nonwet +
476 d_enthalpy_contaminant_nonwet_dT * X_contaminant_nonwet;
477
478 // Assemble M matrix
479 Map.noalias() += porosity * (1 - Sw) *
480 (mol_density_nonwet * d_x_air_nonwet_dpg +
481 x_air_nonwet * d_mol_density_nonwet_dpg) *
482 mass_operator;
483 Mapc.noalias() +=
484 porosity * mol_density_nonwet *
485 ((1 - Sw) * d_x_air_nonwet_dpc - x_air_nonwet * dSw_dpc) *
486 mass_operator;
487 Max.noalias() += porosity * mol_density_nonwet * (1 - Sw) *
488 d_x_air_nonwet_dXc * mass_operator;
489 Mat.noalias() += porosity * (1 - Sw) *
490 (mol_density_nonwet * d_x_air_nonwet_dT +
491 x_air_nonwet * d_mol_density_nonwet_dT) *
492 mass_operator;
493
494 Mwp.noalias() +=
495 porosity *
496 (mol_density_wet * Sw * d_x_water_wet_dpg +
497 (1 - Sw) * x_water_nonwet * d_mol_density_nonwet_dpg +
498 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dpg) *
499 mass_operator;
500 Mwpc.noalias() +=
501 porosity *
502 (mol_density_wet *
503 (x_water_wet * dSw_dpc + Sw * d_x_water_wet_dpc) +
504 mol_density_nonwet *
505 ((1 - Sw) * d_x_water_nonwet_dpc - x_water_nonwet * dSw_dpc)) *
506 mass_operator;
507 Mwx.noalias() +=
508 porosity *
509 (mol_density_wet * Sw * d_x_water_wet_dXc +
510 mol_density_nonwet * (1 - Sw) * d_x_water_nonwet_dXc) *
511 mass_operator;
512 Mwt.noalias() += porosity *
513 ((1 - Sw) / ideal_gas_constant_times_T_int_pt *
514 (d_p_vapour_nonwet_dT - p_vapour_nonwet / T_int_pt)) *
515 mass_operator;
516
517 Mcp.noalias() +=
518 porosity *
519 (mol_density_wet * Sw * d_x_contaminant_wet_dpg +
520 (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dpg +
521 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dpg) *
522 mass_operator;
523 Mcpc.noalias() +=
524 porosity *
525 (mol_density_wet *
526 (x_contaminant_wet * dSw_dpc + Sw * d_x_contaminant_wet_dpc) +
527 mol_density_nonwet * ((1 - Sw) * d_x_contaminant_nonwet_dpc -
528 x_contaminant_nonwet * dSw_dpc)) *
529 mass_operator;
530 Mcx.noalias() +=
531 porosity *
532 (mol_density_wet * Sw * d_x_contaminant_wet_dXc +
533 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dXc) *
534 mass_operator;
535 Mct.noalias() +=
536 porosity *
537 (mol_density_wet * Sw * d_x_contaminant_wet_dT +
538 (1 - Sw) * x_contaminant_nonwet * d_mol_density_nonwet_dT +
539 mol_density_nonwet * (1 - Sw) * d_x_contaminant_nonwet_dT) *
540 mass_operator;
541
542 Mep.noalias() += porosity *
543 (d_density_nonwet_dpg * enthalpy_nonwet - 1) *
544 (1 - Sw) * mass_operator;
545 Mepc.noalias() += porosity *
546 (density_wet * internal_energy_wet -
547 density_nonwet * internal_energy_nonwet) *
548 dSw_dpc * mass_operator +
549 porosity * d_density_nonwet_dpc * enthalpy_nonwet *
550 (1 - Sw) * mass_operator;
551 Mex.noalias() += porosity * d_density_nonwet_dXc * enthalpy_nonwet *
552 (1 - Sw) * mass_operator;
553 Met.noalias() +=
554 ((1 - porosity) * density_solid * heat_capacity_solid +
555 porosity * ((1 - Sw) * (d_density_nonwet_dT * enthalpy_nonwet +
556 density_nonwet * d_enthalpy_nonwet_dT) +
557 Sw * density_wet * heat_capacity_water)) *
558 mass_operator;
559
560 // pore diffusion coefficients
561 double const diffusion_coeff_water_nonwet =
562 water_vapour_component
563 .property(MaterialPropertyLib::PropertyType::pore_diffusion)
564 .template value<double>(vars, pos, t, dt);
565 double const diffusion_coeff_contaminant_nonwet =
566 contaminant_vapour_component
567 .property(MaterialPropertyLib::PropertyType::pore_diffusion)
568 .template value<double>(vars, pos, t, dt);
569 double const diffusion_coeff_contaminant_wet =
570 dissolved_contaminant_component
571 .property(MaterialPropertyLib::PropertyType::pore_diffusion)
572 .template value<double>(vars, pos, t, dt);
573
574 // gas phase relative permeability, viscosity and mobility
575 double const k_rel_nonwet =
576 medium
577 .property(MaterialPropertyLib::PropertyType::
578 relative_permeability_nonwetting_phase)
579 .template value<double>(vars, pos, t, dt);
580 auto const mu_nonwet =
581 gas_phase.property(MaterialPropertyLib::PropertyType::viscosity)
582 .template value<double>(vars, pos, t, dt);
583 double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
584
585 // liquid phase relative permeability, viscosity and mobility
586 double const k_rel_wet =
587 medium
588 .property(
589 MaterialPropertyLib::PropertyType::relative_permeability)
590 .template value<double>(vars, pos, t, dt);
591 auto const mu_wet =
592 liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
593 .template value<double>(vars, pos, t, dt);
594 double const lambda_wet = k_rel_wet / mu_wet;
595
596 // intrinsic permeability
597 auto const permeability =
598 MaterialPropertyLib::formEigenTensor<GlobalDim>(
599 medium.property(MaterialPropertyLib::PropertyType::permeability)
600 .value(vars, pos, t, dt));
601
602 // gravity
603 auto const& b = _process_data.specific_body_force;
604
605 // gas and liquid phase velocities
606 GlobalDimVectorType const velocity_wet =
607 _process_data.has_gravity
608 ? GlobalDimVectorType(
609 -lambda_wet * permeability *
610 (dNdx * (pg_nodal_values - pc_nodal_values) -
611 density_wet * b))
612 : GlobalDimVectorType(-lambda_wet * permeability * dNdx *
613 (pg_nodal_values - pc_nodal_values));
614
615 laplace_operator.noalias() = dNdx.transpose() * permeability * dNdx * w;
616
617 // mechanical dispersivities
618 auto const solute_dispersivity_transverse =
619 medium.template value<double>(
620 MaterialPropertyLib::PropertyType::transversal_dispersivity);
621 auto const solute_dispersivity_longitudinal =
622 medium.template value<double>(
623 MaterialPropertyLib::PropertyType::longitudinal_dispersivity);
624
625 double const velocity_wet_magnitude = velocity_wet.norm();
626 GlobalDimMatrixType const hydrodynamic_dispersion =
627 velocity_wet_magnitude != 0.0
628 ? GlobalDimMatrixType(
629 (porosity * Sw * diffusion_coeff_contaminant_wet +
630 solute_dispersivity_transverse *
631 velocity_wet_magnitude) *
632 I +
633 (solute_dispersivity_longitudinal -
634 solute_dispersivity_transverse) /
635 velocity_wet_magnitude * velocity_wet *
636 velocity_wet.transpose())
637 : GlobalDimMatrixType(
638 (porosity * Sw * diffusion_coeff_contaminant_wet +
639 solute_dispersivity_transverse *
640 velocity_wet_magnitude) *
641 I);
642
643 auto const dispersion_operator =
644 dNdx.transpose() * hydrodynamic_dispersion * dNdx * w;
645
646 // Assemble K matrix
647 // The sum of all diffusive fluxes in either phase must equal to zero
648 Kap.noalias() +=
649 (mol_density_nonwet * x_air_nonwet * lambda_nonwet) *
650 laplace_operator -
651 (1 - Sw) * porosity * mol_density_nonwet *
652 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpg +
653 diffusion_coeff_contaminant_nonwet *
654 d_x_contaminant_nonwet_dpg) *
655 diffusion_operator;
656 Kapc.noalias() +=
657 -(1 - Sw) * porosity * mol_density_nonwet *
658 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dpc +
659 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dpc) *
660 diffusion_operator;
661 Kax.noalias() +=
662 -(1 - Sw) * porosity * mol_density_nonwet *
663 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dXc +
664 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dXc) *
665 diffusion_operator;
666 Kat.noalias() +=
667 -(1 - Sw) * porosity * mol_density_nonwet *
668 (diffusion_coeff_water_nonwet * d_x_water_nonwet_dT +
669 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT) *
670 diffusion_operator;
671
672 Kwp.noalias() +=
673 (mol_density_wet * x_water_wet * lambda_wet +
674 mol_density_nonwet * x_water_nonwet * lambda_nonwet) *
675 laplace_operator +
676 porosity *
677 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
678 d_x_water_wet_dpg +
679 (1 - Sw) * diffusion_coeff_water_nonwet * mol_density_nonwet *
680 d_x_water_nonwet_dpg) *
681 diffusion_operator;
682 Kwpc.noalias() +=
683 -mol_density_wet * x_water_wet * lambda_wet * laplace_operator +
684 porosity *
685 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
686 d_x_water_wet_dpc +
687 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
688 d_x_water_nonwet_dpc) *
689 diffusion_operator;
690 Kwx.noalias() +=
691 porosity *
692 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
693 d_x_water_wet_dXc +
694 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
695 d_x_water_nonwet_dXc) *
696 diffusion_operator;
697 Kwt.noalias() +=
698 porosity *
699 (Sw * mol_density_wet * diffusion_coeff_contaminant_wet *
700 d_x_water_wet_dT +
701 (1 - Sw) * mol_density_nonwet * diffusion_coeff_water_nonwet *
702 d_x_water_nonwet_dT) *
703 diffusion_operator;
704
705 Kcp.noalias() +=
706 (mol_density_wet * x_contaminant_wet * lambda_wet +
707 mol_density_nonwet * x_contaminant_nonwet * lambda_nonwet) *
708 laplace_operator +
709 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpg +
710 porosity * (1 - Sw) * mol_density_nonwet *
711 diffusion_coeff_contaminant_nonwet *
712 d_x_contaminant_nonwet_dpg * diffusion_operator;
713 Kcpc.noalias() +=
714 -mol_density_wet * x_contaminant_wet * lambda_wet *
715 laplace_operator +
716 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dpc +
717 porosity * (1 - Sw) * mol_density_nonwet *
718 diffusion_coeff_contaminant_nonwet *
719 d_x_contaminant_nonwet_dpc * diffusion_operator;
720 Kcx.noalias() +=
721 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dXc +
722 porosity * (1 - Sw) * mol_density_nonwet *
723 diffusion_coeff_contaminant_nonwet *
724 d_x_contaminant_nonwet_dXc * diffusion_operator;
725 Kct.noalias() +=
726 mol_density_wet * dispersion_operator * d_x_contaminant_wet_dT +
727 porosity * (1 - Sw) * mol_density_nonwet *
728 diffusion_coeff_contaminant_nonwet * d_x_contaminant_nonwet_dT *
729 diffusion_operator;
730
731 Kep.noalias() += (lambda_nonwet * density_nonwet * enthalpy_nonwet +
732 lambda_wet * density_wet * enthalpy_wet) *
733 laplace_operator;
734 Kepc.noalias() +=
735 -lambda_wet * enthalpy_wet * density_wet * laplace_operator;
736
737 if (medium.hasProperty(
738 MaterialPropertyLib::PropertyType::thermal_conductivity))
739 {
740 auto const lambda =
741 medium
742 .property(
743 MaterialPropertyLib::PropertyType::thermal_conductivity)
744 .value(vars, pos, t, dt);
745
746 GlobalDimMatrixType const heat_conductivity_unsaturated =
747 MaterialPropertyLib::formEigenTensor<GlobalDim>(lambda);
748
749 Ket.noalias() +=
750 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
751 }
752 else
753 {
754 auto const thermal_conductivity_solid =
755 solid_phase
756 .property(
757 MaterialPropertyLib::PropertyType::thermal_conductivity)
758 .value(vars, pos, t, dt);
759
760 auto const thermal_conductivity_fluid =
761 liquid_phase
762 .property(
763 MaterialPropertyLib::PropertyType::thermal_conductivity)
764 .template value<double>(vars, pos, t, dt) *
765 Sw;
766
767 GlobalDimMatrixType const heat_conductivity_unsaturated =
768 MaterialPropertyLib::formEffectiveThermalConductivity<
769 GlobalDim>(thermal_conductivity_solid,
770 thermal_conductivity_fluid, porosity);
771
772 Ket.noalias() +=
773 dNdx.transpose() * heat_conductivity_unsaturated * dNdx * w;
774 }
775
776 if (_process_data.has_gravity)
777 {
778 NodalVectorType gravity_operator =
779 dNdx.transpose() * permeability * b * w;
780 Ba.noalias() += (mol_density_nonwet * x_air_nonwet * lambda_nonwet *
781 density_nonwet) *
782 gravity_operator;
783 Bw.noalias() +=
784 (mol_density_wet * x_water_wet * lambda_wet * density_wet +
785 mol_density_nonwet * x_water_nonwet * lambda_nonwet *
786 density_nonwet) *
787 gravity_operator;
788 Bc.noalias() += (mol_density_wet * x_contaminant_wet * lambda_wet *
789 density_wet +
790 mol_density_nonwet * x_contaminant_nonwet *
791 lambda_nonwet * density_nonwet) *
792 gravity_operator;
793 Be.noalias() +=
794 (lambda_nonwet * density_nonwet * density_nonwet *
795 enthalpy_nonwet +
796 lambda_wet * density_wet * density_wet * enthalpy_wet) *
797 gravity_operator;
798 } // end of has gravity
799 }
800 if (_process_data.has_mass_lumping)
801 {
802 Map = Map.colwise().sum().eval().asDiagonal();
803 Mapc = Mapc.colwise().sum().eval().asDiagonal();
804 Max = Max.colwise().sum().eval().asDiagonal();
805 Mat = Mat.colwise().sum().eval().asDiagonal();
806 Mwp = Mwp.colwise().sum().eval().asDiagonal();
807 Mwpc = Mwpc.colwise().sum().eval().asDiagonal();
808 Mwx = Mwx.colwise().sum().eval().asDiagonal();
809 Mwt = Mwt.colwise().sum().eval().asDiagonal();
810 Mcp = Mcp.colwise().sum().eval().asDiagonal();
811 Mcpc = Mcpc.colwise().sum().eval().asDiagonal();
812 Mcx = Mcx.colwise().sum().eval().asDiagonal();
813 Mct = Mct.colwise().sum().eval().asDiagonal();
814 Mep = Mep.colwise().sum().eval().asDiagonal();
815 Mepc = Mepc.colwise().sum().eval().asDiagonal();
816 Mex = Mex.colwise().sum().eval().asDiagonal();
817 Met = Met.colwise().sum().eval().asDiagonal();
818 } // end of mass-lumping
819}
MaterialPropertyLib::formEigenTensor
constexpr Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.h:163
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:146
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:25
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:20

References _element, _gas_molar_fraction_contaminant, _gas_molar_fraction_water, _integration_method, _ip_data, _liquid_molar_fraction_contaminant, _pressure_wetting, _process_data, _saturation, cap_pressure_matrix_index, MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::VariableArray::capillary_pressure, MaterialLib::PhysicalConstant::CelsiusZeroInKelvin, contaminant_matrix_index, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::enthalpy_of_evaporation, MaterialPropertyLib::formEffectiveThermalConductivity(), MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::henry_coefficient, MaterialLib::PhysicalConstant::IdealGasConstant, NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::longitudinal_dispersivity, MaterialPropertyLib::molar_mass, MaterialPropertyLib::VariableArray::molar_mass, nonwet_pressure_matrix_index, ProcessLib::ThermalTwoPhaseFlowWithPP::NUM_NODAL_DOF, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::porosity, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::specific_heat_capacity, MaterialPropertyLib::specific_latent_heat, MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, temperature_matrix_index, MaterialPropertyLib::thermal_conductivity, MaterialPropertyLib::transversal_dispersivity, MaterialPropertyLib::vapour_pressure, and MaterialPropertyLib::viscosity.

◆ getIntPtGasMolFracContaminant()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::getIntPtGasMolFracContaminant ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > &  ) const
inlineoverridevirtual

Implements ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssemblerInterface.

Definition at line 205 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

210 {
211 assert(!_gas_molar_fraction_contaminant.empty());
212 return _gas_molar_fraction_contaminant;
213 }

References _gas_molar_fraction_contaminant.

◆ getIntPtGasMolFracWater()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::getIntPtGasMolFracWater ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > &  ) const
inlineoverridevirtual

Implements ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssemblerInterface.

Definition at line 195 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

200 {
201 assert(!_gas_molar_fraction_water.empty());
202 return _gas_molar_fraction_water;
203 }

References _gas_molar_fraction_water.

◆ getIntPtLiquidMolFracContaminant()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::getIntPtLiquidMolFracContaminant ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > &  ) const
inlineoverridevirtual

Implements ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssemblerInterface.

Definition at line 185 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

190 {
191 assert(!_liquid_molar_fraction_contaminant.empty());
192 return _liquid_molar_fraction_contaminant;
193 }

References _liquid_molar_fraction_contaminant.

◆ getIntPtSaturation()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::getIntPtSaturation ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > &  ) const
inlineoverridevirtual

Implements ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssemblerInterface.

Definition at line 165 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

170 {
171 assert(!_saturation.empty());
172 return _saturation;
173 }

References _saturation.

◆ getIntPtWettingPressure()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::getIntPtWettingPressure ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > &  ) const
inlineoverridevirtual

Implements ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssemblerInterface.

Definition at line 175 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

180 {
181 assert(!_pressure_wetting.empty());
182 return _pressure_wetting;
183 }

References _pressure_wetting.

◆ getShapeMatrix()

template<typename ShapeFunction, int GlobalDim>
Eigen::Map< const Eigen::RowVectorXd > ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::getShapeMatrix ( const unsigned integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 156 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

158 {
159 auto const& N = _ip_data[integration_point].N;
160
161 // assumes N is stored contiguously in memory
162 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
163 }

References _ip_data.

Member Data Documentation

◆ _element

template<typename ShapeFunction, int GlobalDim>
MeshLib::Element const& ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_element
private

Definition at line 216 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), and assemble().

◆ _gas_molar_fraction_contaminant

template<typename ShapeFunction, int GlobalDim>
std::vector<double> ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_gas_molar_fraction_contaminant
private

Definition at line 232 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), assemble(), and getIntPtGasMolFracContaminant().

◆ _gas_molar_fraction_water

template<typename ShapeFunction, int GlobalDim>
std::vector<double> ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_gas_molar_fraction_water
private

Definition at line 231 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), assemble(), and getIntPtGasMolFracWater().

◆ _integration_method

template<typename ShapeFunction, int GlobalDim>
NumLib::GenericIntegrationMethod const& ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_integration_method
private

Definition at line 218 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), and assemble().

◆ _ip_data

template<typename ShapeFunction, int GlobalDim>
std::vector< IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>, Eigen::aligned_allocator<IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType> > > ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_ip_data
private

Definition at line 226 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), assemble(), and getShapeMatrix().

◆ _liquid_molar_fraction_contaminant

template<typename ShapeFunction, int GlobalDim>
std::vector<double> ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_liquid_molar_fraction_contaminant
private

Definition at line 230 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), assemble(), and getIntPtLiquidMolFracContaminant().

◆ _pressure_wetting

template<typename ShapeFunction, int GlobalDim>
std::vector<double> ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_pressure_wetting
private

Definition at line 229 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), assemble(), and getIntPtWettingPressure().

◆ _process_data

template<typename ShapeFunction, int GlobalDim>
ThermalTwoPhaseFlowWithPPProcessData const& ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_process_data
private

Definition at line 220 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), and assemble().

◆ _saturation

template<typename ShapeFunction, int GlobalDim>
std::vector<double> ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::_saturation
private

Definition at line 228 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by ThermalTwoPhaseFlowWithPPLocalAssembler(), assemble(), and getIntPtSaturation().

◆ cap_pressure_matrix_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::cap_pressure_matrix_index = ShapeFunction::NPOINTS
staticprivate

Definition at line 235 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by assemble().

◆ cap_pressure_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::cap_pressure_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 240 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ contaminant_matrix_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::contaminant_matrix_index = 2 * ShapeFunction::NPOINTS
staticprivate

Definition at line 236 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by assemble().

◆ contaminant_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::contaminant_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 241 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ nonwet_pressure_matrix_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::nonwet_pressure_matrix_index = 0
staticprivate

Definition at line 234 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by assemble().

◆ nonwet_pressure_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::nonwet_pressure_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 239 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ temperature_matrix_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::temperature_matrix_index = 3 * ShapeFunction::NPOINTS
staticprivate

Definition at line 237 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

Referenced by assemble().

◆ temperature_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermalTwoPhaseFlowWithPP::ThermalTwoPhaseFlowWithPPLocalAssembler< ShapeFunction, GlobalDim >::temperature_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 242 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

The documentation for this class was generated from the following files: