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 94 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 std::optional< VectorSegment > getVectorDeformationSegment () 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 108 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:82

Definition at line 104 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ GlobalDimVectorType

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

Definition at line 109 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:42
ProcessLib::LocalAssemblerTraits
detail::LocalAssemblerTraitsFixed< ShpPol, NNodes, NodalDOF, Dim > LocalAssemblerTraits
Definition LocalAssemblerTraits.h:185

Definition at line 100 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ LocalMatrixType

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

Definition at line 110 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ LocalVectorType

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

Definition at line 111 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ NodalMatrixType

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

Definition at line 106 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ NodalRowVectorType

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

Definition at line 102 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ NodalVectorType

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

Definition at line 107 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ ShapeMatrices

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

Definition at line 98 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

◆ ShapeMatricesType

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

Definition at line 97 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 114 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

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

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

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

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

217 {
218 assert(!_gas_molar_fraction_contaminant.empty());
219 return _gas_molar_fraction_contaminant;
220 }

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

207 {
208 assert(!_gas_molar_fraction_water.empty());
209 return _gas_molar_fraction_water;
210 }

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

197 {
198 assert(!_liquid_molar_fraction_contaminant.empty());
199 return _liquid_molar_fraction_contaminant;
200 }

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

177 {
178 assert(!_saturation.empty());
179 return _saturation;
180 }

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

187 {
188 assert(!_pressure_wetting.empty());
189 return _pressure_wetting;
190 }

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

165 {
166 auto const& N = _ip_data[integration_point].N;
167
168 // assumes N is stored contiguously in memory
169 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
170 }

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 223 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 239 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 238 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 225 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 233 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 237 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 236 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 227 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 235 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 242 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 247 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 243 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 248 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 241 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 246 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 244 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 249 of file ThermalTwoPhaseFlowWithPPLocalAssembler.h.

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