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 21 of file ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h.

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