35 namespace TwoPhaseFlowWithPrho
39 std::unique_ptr<MaterialLib::Fluid::FluidProperty>
41 std::unique_ptr<MaterialLib::Fluid::FluidProperty>
43 std::unique_ptr<MaterialLib::Fluid::FluidProperty>
45 std::unique_ptr<MaterialLib::Fluid::FluidProperty>
47 std::vector<std::unique_ptr<MaterialLib::PorousMedium::Permeability>>&&
48 intrinsic_permeability_models,
49 std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
51 std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
53 std::vector<std::unique_ptr<
55 capillary_pressure_models,
57 std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
58 relative_permeability_models)
59 : _liquid_density(std::move(liquid_density)),
61 _gas_density(std::move(gas_density)),
62 _gas_viscosity(std::move(gas_viscosity)),
63 _material_ids(material_ids),
64 _intrinsic_permeability_models(std::move(intrinsic_permeability_models)),
65 _porosity_models(std::move(porosity_models)),
66 _storage_models(std::move(storage_models)),
67 _capillary_pressure_models(std::move(capillary_pressure_models)),
68 _relative_permeability_models(std::move(relative_permeability_models))
70 DBUG(
"Create material properties for Two-Phase flow with P-RHO model.");
74 const std::size_t element_id)
81 assert(element_id < _material_ids->size());
86 const double p,
const double T)
const
95 const double p,
const double T)
const
104 const double p,
const double T)
const
113 const double p,
const double T)
const
122 const int material_id,
const double t,
130 const int material_id,
const double t,
132 const double T,
const double porosity_variable)
const
140 const double ,
const double ,
const double saturation)
const
147 const double ,
const double ,
const double saturation)
const
153 const int material_id,
const double ,
162 const int material_id,
const double ,
172 int const material_id,
184 using LocalJacobianMatrix =
185 Eigen::Matrix<double, 2, 2, Eigen::RowMajor>;
186 using LocalResidualVector = Eigen::Matrix<double, 2, 1>;
187 using LocalUnknownVector = Eigen::Matrix<double, 2, 1>;
188 LocalJacobianMatrix J_loc;
190 Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(2);
191 auto const update_residual = [&](LocalResidualVector& residual)
194 auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
200 auto const update_solution = [&](LocalUnknownVector
const& increment)
210 const int maximum_iterations(20);
211 const double residuum_tolerance(1.e-14);
212 const double increment_tolerance(0);
215 decltype(linear_solver), LocalJacobianMatrix,
216 decltype(update_jacobian), LocalResidualVector,
217 decltype(update_residual), decltype(update_solution)>(
218 linear_solver, update_jacobian, update_residual, update_solution,
219 {maximum_iterations, residuum_tolerance, increment_tolerance});
221 auto const success_iterations = newton_solver.
solve(J_loc);
223 if (!success_iterations)
235 const int material_id,
double const pl,
double const X,
double const T,
248 const int material_id,
double const ,
250 double const ,
double const T,
JacobianMatrix& Jac,
double const Sw,
251 double const rho_h2_wet)
257 double const dPC_dSw =
261 if ((1 - Sw) < (rho_equili_h2_wet - rho_h2_wet))
267 Jac(0, 0) = drhoh2wet_dpg * dPC_dSw;
271 Jac(1, 0) = rho_h2_nonwet - rho_h2_wet;
279 double const pg,
double const Sw,
double const rho_wet_h2)
282 return std::min(1 - Sw, rho_equilibrium_wet_h2 - rho_wet_h2);
289 double ,
double X,
double Sw,
double rho_wet_h2,
double rho_nonwet_h2,
292 return X - (Sw * rho_wet_h2 + (1 - Sw) * rho_nonwet_h2);
299 double pl,
double S,
double rho_wet_h2,
double const T,
300 int current_material_id)
const
307 if ((1 - S) < (rho_equilibrium_wet_h2 - rho_wet_h2))
314 ((drhoh2nonwet_dpg - drhoh2wet_dpg) * (1 - S) + drhoh2wet_dpg);
315 double const beta = (drhoh2nonwet_dpg - drhoh2wet_dpg) *
317 double const dPC_dSw =
325 double const pl,
const double ,
const double S,
326 const double rho_wet_h2,
double const T,
int current_material_id)
const
333 if ((1 - S) < (rho_equilibrium_wet_h2 - rho_wet_h2))
340 ((drhoh2nonwet_dpg - drhoh2wet_dpg) * (1 - S) + drhoh2wet_dpg);
341 double const beta = (drhoh2nonwet_dpg - drhoh2wet_dpg) *
343 double const dPC_dSw =
351 double pl,
double Sw,
double rho_wet_h2,
double dSwdX,
352 int current_material_id)
const
359 double const dPC_dSw =
361 if ((1 - Sw) < (rho_equilibrium_wet_h2 - rho_wet_h2))
371 double pl,
double Sw,
double rho_wet_h2,
double dSwdP,
372 int current_material_id)
const
379 double const dPC_dSw =
381 if ((1 - Sw) < (rho_equilibrium_wet_h2 - rho_wet_h2))
void DBUG(char const *fmt, Args const &... args)
Definition of the Mesh class.
Definition of the PiecewiseLinearInterpolation class.
Base class of capillary pressure models.
std::optional< int > solve(JacobianMatrix &jacobian) const
static double calculateSaturation(double, double X, double Sw, double rho_wet_h2, double rho_nonwet_h2, double)
double calculatedXmdP(double pl, double Sw, double rho_wet_h2, double dSwdP, int current_material_id) const
TwoPhaseFlowWithPrhoMaterialProperties(MeshLib::PropertyVector< int > const *material_ids, std::unique_ptr< MaterialLib::Fluid::FluidProperty > liquid_density, std::unique_ptr< MaterialLib::Fluid::FluidProperty > viscosity, std::unique_ptr< MaterialLib::Fluid::FluidProperty > gas_density, std::unique_ptr< MaterialLib::Fluid::FluidProperty > gas_viscosity, std::vector< std::unique_ptr< MaterialLib::PorousMedium::Permeability >> &&intrinsic_permeability_models, std::vector< std::unique_ptr< MaterialLib::PorousMedium::Porosity >> &&porosity_models, std::vector< std::unique_ptr< MaterialLib::PorousMedium::Storage >> &&storage_models, std::vector< std::unique_ptr< MaterialLib::PorousMedium::CapillaryPressureSaturation >> &&capillary_pressure_models, std::vector< std::unique_ptr< MaterialLib::PorousMedium::RelativePermeability >> &&relative_permeability_models)
MaterialLib::Fluid::FluidProperty::ArrayType ArrayType
std::vector< std::unique_ptr< MaterialLib::PorousMedium::CapillaryPressureSaturation > > _capillary_pressure_models
double calculatedSwdX(double const pl, const double, const double S, const double rho_wet_h2, double const T, int current_material_id) const
std::vector< std::unique_ptr< MaterialLib::PorousMedium::Permeability > > _intrinsic_permeability_models
double getNonwetRelativePermeability(const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
void calculateResidual(const int material_id, double const pl, double const X, double const T, double Sw, double rho_h2_wet, ResidualVector &res)
Eigen::Matrix< double, jacobian_residual_size, 1 > ResidualVector
double getPorosity(const int material_id, const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double porosity_variable) const
bool computeConstitutiveRelation(double const t, ParameterLib::SpatialPosition const &x_position, 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)
double getWetRelativePermeability(const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
std::unique_ptr< MaterialLib::Fluid::FluidProperty > _gas_density
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)
std::unique_ptr< MaterialLib::Fluid::FluidProperty > _liquid_density
double getGasViscosity(const double p, const double T) const
std::unique_ptr< MaterialLib::Fluid::FluidProperty > _viscosity
std::unique_ptr< MaterialLib::Fluid::FluidProperty > _gas_viscosity
static double calculateEquilibiumRhoWetLight(double const pg, double const Sw, double const rho_wet_h2)
std::vector< std::unique_ptr< MaterialLib::PorousMedium::Porosity > > _porosity_models
double getCapillaryPressureDerivative(const int material_id, const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
double calculatedSwdP(double pl, double S, double rho_wet_h2, double const T, int current_material_id) const
int getMaterialID(const std::size_t element_id)
Eigen::Matrix< double, jacobian_residual_size, jacobian_residual_size, Eigen::RowMajor > JacobianMatrix
double getCapillaryPressure(const int material_id, const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
double getLiquidDensity(const double p, const double T) const
Eigen::MatrixXd getPermeability(const int material_id, const double t, const ParameterLib::SpatialPosition &pos, const int dim) const
double getLiquidViscosity(const double p, const double T) const
double getGasDensity(const double p, const double T) const
MeshLib::PropertyVector< int > const *const _material_ids
double calculatedXmdX(double pl, double Sw, double rho_wet_h2, double dSwdX, int current_material_id) const
std::vector< std::unique_ptr< MaterialLib::PorousMedium::RelativePermeability > > _relative_permeability_models
constexpr double HenryConstantH2
constexpr double H2
kg mol-1
constexpr double IdealGasConstant