OGS
ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties Class Reference

Detailed Description

Definition at line 33 of file TwoPhaseFlowWithPrhoMaterialProperties.h.

#include <TwoPhaseFlowWithPrhoMaterialProperties.h>

Collaboration diagram for ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties:
[legend]

Public Types

using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType
 

Public Member Functions

 TwoPhaseFlowWithPrhoMaterialProperties (MeshLib::PropertyVector< int > const *material_ids, std::vector< std::unique_ptr< MaterialLib::PorousMedium::CapillaryPressureSaturation > > &&capillary_pressure_models, std::vector< std::unique_ptr< MaterialLib::PorousMedium::RelativePermeability > > &&relative_permeability_models)
 
int getMaterialID (const std::size_t element_id)
 
double getNonwetRelativePermeability (const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
 
double getWetRelativePermeability (const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
 
double getCapillaryPressure (const int material_id, const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
 
double getCapillaryPressureDerivative (const int material_id, const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
 
bool computeConstitutiveRelation (double const t, ParameterLib::SpatialPosition const &x, const int material_id, double const pg, double const X, double const T, double &Sw, double &X_m, double &dsw_dpg, double &dsw_dX, double &dxm_dpg, double &dxm_dX)
 

Protected Attributes

MeshLib::PropertyVector< int > const *const _material_ids
 
std::vector< std::unique_ptr< MaterialLib::PorousMedium::CapillaryPressureSaturation > > _capillary_pressure_models
 
std::vector< std::unique_ptr< MaterialLib::PorousMedium::RelativePermeability > > _relative_permeability_models
 

Private Types

using ResidualVector = Eigen::Matrix<double, jacobian_residual_size, 1>
 
using JacobianMatrix
 
using UnknownVector = Eigen::Matrix<double, jacobian_residual_size, 1>
 

Private Member Functions

void calculateResidual (const int material_id, double const pl, double const X, double const T, double Sw, double rho_h2_wet, ResidualVector &res)
 
void calculateJacobian (const int material_id, double const t, ParameterLib::SpatialPosition const &x, double const pl, double const X, double const T, JacobianMatrix &Jac, double Sw, double rho_h2_wet)
 
double calculatedSwdP (double pl, double S, double rho_wet_h2, double const T, int current_material_id) const
 
double calculatedSwdX (double const pl, const double, const double S, const double rho_wet_h2, double const T, int current_material_id) const
 
double calculatedXmdX (double pl, double Sw, double rho_wet_h2, double dSwdX, int current_material_id) const
 
double calculatedXmdP (double pl, double Sw, double rho_wet_h2, double dSwdP, int current_material_id) const
 

Static Private Member Functions

static double calculateEquilibiumRhoWetLight (double const pg, double const Sw, double const rho_wet_h2)
 
static double calculateSaturation (double, double X, double Sw, double rho_wet_h2, double rho_nonwet_h2, double)
 

Static Private Attributes

static int const jacobian_residual_size = 2
 

Member Typedef Documentation

◆ ArrayType

◆ JacobianMatrix

◆ ResidualVector

◆ UnknownVector

Constructor & Destructor Documentation

◆ TwoPhaseFlowWithPrhoMaterialProperties()

ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::TwoPhaseFlowWithPrhoMaterialProperties ( MeshLib::PropertyVector< int > const * material_ids,
std::vector< std::unique_ptr< MaterialLib::PorousMedium::CapillaryPressureSaturation > > && capillary_pressure_models,
std::vector< std::unique_ptr< MaterialLib::PorousMedium::RelativePermeability > > && relative_permeability_models )

Definition at line 36 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

44 : _material_ids(material_ids),
45 _capillary_pressure_models(std::move(capillary_pressure_models)),
46 _relative_permeability_models(std::move(relative_permeability_models))
47{
48 DBUG("Create material properties for Two-Phase flow with P-RHO model.");
49}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
std::vector< std::unique_ptr< MaterialLib::PorousMedium::CapillaryPressureSaturation > > _capillary_pressure_models
std::vector< std::unique_ptr< MaterialLib::PorousMedium::RelativePermeability > > _relative_permeability_models

References DBUG().

Member Function Documentation

◆ calculatedSwdP()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::calculatedSwdP ( double pl,
double S,
double rho_wet_h2,
double const T,
int current_material_id ) const
private

Calculate the derivatives using the analytical way

Definition at line 223 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

226{
227 const double pg =
228 pl +
229 _capillary_pressure_models[current_material_id]->getCapillaryPressure(
230 S);
231 double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
232 if ((1 - S) < (rho_equilibrium_wet_h2 - rho_wet_h2))
233 {
234 return 0.0;
235 }
236 double const drhoh2wet_dpg = HenryConstantH2 * H2;
237 double const drhoh2nonwet_dpg = H2 / IdealGasConstant / T;
238 double const alpha =
239 ((drhoh2nonwet_dpg - drhoh2wet_dpg) * (1 - S) + drhoh2wet_dpg);
240 double const beta = (drhoh2nonwet_dpg - drhoh2wet_dpg) *
241 pg; // NOTE here should be PG^h, but we ignore vapor
242 double const dPC_dSw =
243 _capillary_pressure_models[current_material_id]->getdPcdS(S);
244 return alpha / (beta - alpha * dPC_dSw);
245}

References _capillary_pressure_models.

Referenced by computeConstitutiveRelation().

◆ calculatedSwdX()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::calculatedSwdX ( double const pl,
const double ,
const double S,
const double rho_wet_h2,
double const T,
int current_material_id ) const
private

Calculate the derivatives using the analytical way

Definition at line 249 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

252{
253 const double pg =
254 pl +
255 _capillary_pressure_models[current_material_id]->getCapillaryPressure(
256 S);
257 double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
258 if ((1 - S) < (rho_equilibrium_wet_h2 - rho_wet_h2))
259 {
260 return 0.0;
261 }
262 double const drhoh2wet_dpg = HenryConstantH2 * H2;
263 double const drhoh2nonwet_dpg = H2 / IdealGasConstant / T;
264 double const alpha =
265 ((drhoh2nonwet_dpg - drhoh2wet_dpg) * (1 - S) + drhoh2wet_dpg);
266 double const beta = (drhoh2nonwet_dpg - drhoh2wet_dpg) *
267 pg; // NOTE here should be PG^h, but we ignore vapor
268 double const dPC_dSw =
269 _capillary_pressure_models[current_material_id]->getdPcdS(S);
270 return -1 / (beta - alpha * dPC_dSw);
271}

References _capillary_pressure_models.

Referenced by computeConstitutiveRelation().

◆ calculatedXmdP()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::calculatedXmdP ( double pl,
double Sw,
double rho_wet_h2,
double dSwdP,
int current_material_id ) const
private

Calculate the derivatives using the analytical way

Definition at line 295 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

298{
299 const double pg =
300 pl +
301 _capillary_pressure_models[current_material_id]->getCapillaryPressure(
302 Sw);
303 double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
304 double const dPC_dSw =
305 _capillary_pressure_models[current_material_id]->getdPcdS(Sw);
306 if ((1 - Sw) < (rho_equilibrium_wet_h2 - rho_wet_h2))
307 {
308 return 0.0;
309 }
310 return HenryConstantH2 * H2 * (1 + dPC_dSw * dSwdP);
311}

References _capillary_pressure_models.

Referenced by computeConstitutiveRelation().

◆ calculatedXmdX()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::calculatedXmdX ( double pl,
double Sw,
double rho_wet_h2,
double dSwdX,
int current_material_id ) const
private

Calculate the derivatives using the analytical way

Definition at line 275 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

278{
279 const double pg =
280 pl +
281 _capillary_pressure_models[current_material_id]->getCapillaryPressure(
282 Sw);
283 double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
284 double const dPC_dSw =
285 _capillary_pressure_models[current_material_id]->getdPcdS(Sw);
286 if ((1 - Sw) < (rho_equilibrium_wet_h2 - rho_wet_h2))
287 {
288 return 1.0;
289 }
290 return HenryConstantH2 * H2 * dPC_dSw * dSwdX;
291}

References _capillary_pressure_models.

Referenced by computeConstitutiveRelation().

◆ calculateEquilibiumRhoWetLight()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::calculateEquilibiumRhoWetLight ( double const pg,
double const Sw,
double const rho_wet_h2 )
staticprivate

Complementary condition 1 for calculating molar fraction of light component in the liquid phase

Definition at line 203 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

205{
206 double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
207 return std::min(1 - Sw, rho_equilibrium_wet_h2 - rho_wet_h2);
208}

Referenced by calculateResidual().

◆ calculateJacobian()

void ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::calculateJacobian ( const int material_id,
double const t,
ParameterLib::SpatialPosition const & x,
double const pl,
double const X,
double const T,
JacobianMatrix & Jac,
double Sw,
double rho_h2_wet )
private

Calculates the Jacobian.

Definition at line 172 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

177{
178 const double pg =
179 pl + _capillary_pressure_models[material_id]->getCapillaryPressure(Sw);
180 const double rho_h2_nonwet = pg * H2 / IdealGasConstant / T;
181 double const rho_equili_h2_wet = pg * HenryConstantH2 * H2;
182 double const dPC_dSw =
183 _capillary_pressure_models[material_id]->getdPcdS(Sw);
184 double const drhoh2wet_dpg = HenryConstantH2 * H2;
185 Jac.setZero();
186 if ((1 - Sw) < (rho_equili_h2_wet - rho_h2_wet))
187 {
188 Jac(0, 0) = -1;
189 }
190 else
191 {
192 Jac(0, 0) = drhoh2wet_dpg * dPC_dSw;
193 Jac(0, 1) = -1;
194 }
195
196 Jac(1, 0) = rho_h2_nonwet - rho_h2_wet;
197 Jac(1, 1) = -Sw;
198}

References _capillary_pressure_models.

Referenced by computeConstitutiveRelation().

◆ calculateResidual()

void ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::calculateResidual ( const int material_id,
double const pl,
double const X,
double const T,
double Sw,
double rho_h2_wet,
ResidualVector & res )
private

Calculates the residual vector.

Definition at line 159 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

162{
163 const double pg =
164 pl + _capillary_pressure_models[material_id]->getCapillaryPressure(Sw);
165 const double rho_h2_nonwet = pg * H2 / IdealGasConstant / T;
166
167 // calculating residual
168 res[0] = calculateEquilibiumRhoWetLight(pg, Sw, rho_h2_wet);
169 res[1] = calculateSaturation(pl, X, Sw, rho_h2_wet, rho_h2_nonwet, T);
170}
static double calculateSaturation(double, double X, double Sw, double rho_wet_h2, double rho_nonwet_h2, double)
static double calculateEquilibiumRhoWetLight(double const pg, double const Sw, double const rho_wet_h2)

References _capillary_pressure_models, calculateEquilibiumRhoWetLight(), and calculateSaturation().

Referenced by computeConstitutiveRelation().

◆ calculateSaturation()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::calculateSaturation ( double ,
double X,
double Sw,
double rho_wet_h2,
double rho_nonwet_h2,
double  )
staticprivate

Complementary condition 2 for calculating the saturation

Definition at line 213 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

216{
217 return X - (Sw * rho_wet_h2 + (1 - Sw) * rho_nonwet_h2);
218}

Referenced by calculateResidual().

◆ computeConstitutiveRelation()

bool ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::computeConstitutiveRelation ( double const t,
ParameterLib::SpatialPosition const & x,
const int material_id,
double const pg,
double const X,
double const T,
double & Sw,
double & X_m,
double & dsw_dpg,
double & dsw_dX,
double & dxm_dpg,
double & dxm_dX )

Definition at line 94 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

107{
108 { // Local Newton solver
109 using LocalJacobianMatrix =
110 Eigen::Matrix<double, 2, 2, Eigen::RowMajor>;
111 using LocalResidualVector = Eigen::Matrix<double, 2, 1>;
112 using LocalUnknownVector = Eigen::Matrix<double, 2, 1>;
113 LocalJacobianMatrix J_loc;
114
115 Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(2);
116 auto const update_residual = [&](LocalResidualVector& residual)
117 { calculateResidual(material_id, pg, X, T, Sw, X_m, residual); };
118
119 auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
120 {
121 calculateJacobian(material_id, t, x, pg, X, T, jacobian, Sw,
122 X_m); // for solution dependent Jacobians
123 };
124
125 auto const update_solution = [&](LocalUnknownVector const& increment)
126 {
127 // increment solution vectors
128 Sw += increment[0];
129 X_m += increment[1];
130 };
131
132 // TODO Make the following choice of maximum iterations and convergence
133 // criteria available from the input file configuration. See Ehlers
134 // material model implementation for the example.
135 const int maximum_iterations(20);
136 const double residuum_tolerance(1.e-14);
137 const double increment_tolerance(0);
138
139 auto newton_solver = NumLib::NewtonRaphson<
140 decltype(linear_solver), LocalJacobianMatrix,
141 decltype(update_jacobian), LocalResidualVector,
142 decltype(update_residual), decltype(update_solution)>(
143 linear_solver, update_jacobian, update_residual, update_solution,
144 {maximum_iterations, residuum_tolerance, increment_tolerance});
145
146 auto const success_iterations = newton_solver.solve(J_loc);
147
148 if (!success_iterations)
149 {
150 return false;
151 }
152 }
153 dsw_dpg = calculatedSwdP(pg, Sw, X_m, T, material_id);
154 dsw_dX = calculatedSwdX(pg, X, Sw, X_m, T, material_id);
155 dxm_dpg = calculatedXmdP(pg, Sw, X_m, dsw_dpg, material_id);
156 dxm_dX = calculatedXmdX(pg, Sw, X_m, dsw_dX, material_id);
157 return true;
158}
std::optional< int > solve(JacobianMatrix &jacobian) const
double calculatedXmdP(double pl, double Sw, double rho_wet_h2, double dSwdP, int current_material_id) const
double calculatedSwdX(double const pl, const double, const double S, const double rho_wet_h2, double const T, int current_material_id) const
void calculateResidual(const int material_id, double const pl, double const X, double const T, double Sw, double rho_h2_wet, ResidualVector &res)
void calculateJacobian(const int material_id, double const t, ParameterLib::SpatialPosition const &x, double const pl, double const X, double const T, JacobianMatrix &Jac, double Sw, double rho_h2_wet)
double calculatedSwdP(double pl, double S, double rho_wet_h2, double const T, int current_material_id) const
double calculatedXmdX(double pl, double Sw, double rho_wet_h2, double dSwdX, int current_material_id) const

References calculatedSwdP(), calculatedSwdX(), calculatedXmdP(), calculatedXmdX(), calculateJacobian(), calculateResidual(), and NumLib::NewtonRaphson< LinearSolver, JacobianMatrix, JacobianMatrixUpdate, ResidualVector, ResidualUpdate, SolutionUpdate >::solve().

◆ getCapillaryPressure()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::getCapillaryPressure ( const int material_id,
const double t,
const ParameterLib::SpatialPosition & pos,
const double p,
const double T,
const double saturation ) const

Definition at line 77 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

81{
82 return _capillary_pressure_models[material_id]->getCapillaryPressure(
83 saturation);
84}

References _capillary_pressure_models.

◆ getCapillaryPressureDerivative()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::getCapillaryPressureDerivative ( const int material_id,
const double t,
const ParameterLib::SpatialPosition & pos,
const double p,
const double T,
const double saturation ) const

Definition at line 86 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

90{
91 return _capillary_pressure_models[material_id]->getdPcdS(saturation);
92}

References _capillary_pressure_models.

◆ getMaterialID()

int ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::getMaterialID ( const std::size_t element_id)

Definition at line 51 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

53{
54 if (!_material_ids)
55 {
56 return 0;
57 }
58
59 assert(element_id < _material_ids->size());
60 return (*_material_ids)[element_id];
61}
constexpr int size(int const displacement_dim)
Vectorized tensor size for given displacement dimension.

References _material_ids.

◆ getNonwetRelativePermeability()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::getNonwetRelativePermeability ( const double t,
const ParameterLib::SpatialPosition & pos,
const double p,
const double T,
const double saturation ) const

Definition at line 63 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

66{
67 return _relative_permeability_models[0]->getValue(saturation);
68}

References _relative_permeability_models.

◆ getWetRelativePermeability()

double ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::getWetRelativePermeability ( const double t,
const ParameterLib::SpatialPosition & pos,
const double p,
const double T,
const double saturation ) const

Definition at line 70 of file TwoPhaseFlowWithPrhoMaterialProperties.cpp.

73{
74 return _relative_permeability_models[1]->getValue(saturation);
75}

References _relative_permeability_models.

Member Data Documentation

◆ _capillary_pressure_models

std::vector< std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation> > ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::_capillary_pressure_models
protected

◆ _material_ids

MeshLib::PropertyVector<int> const* const ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::_material_ids
protected

Use two phase models for different material zones. Material IDs must be given as mesh element properties.

Definition at line 81 of file TwoPhaseFlowWithPrhoMaterialProperties.h.

Referenced by getMaterialID().

◆ _relative_permeability_models

std::vector< std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability> > ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::_relative_permeability_models
protected

◆ jacobian_residual_size

int const ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties::jacobian_residual_size = 2
staticprivate

Definition at line 91 of file TwoPhaseFlowWithPrhoMaterialProperties.h.


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