OGS
TwoPhaseFlowWithPrhoMaterialProperties.cpp
Go to the documentation of this file.
1
12
13#include <Eigen/LU>
14#include <utility>
15
16#include "BaseLib/Logging.h"
23#include "MeshLib/Mesh.h"
28
32namespace ProcessLib
33{
34namespace TwoPhaseFlowWithPrho
35{
37 MeshLib::PropertyVector<int> const* const material_ids,
38 std::vector<std::unique_ptr<
40 capillary_pressure_models,
41 std::vector<
42 std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
43 relative_permeability_models)
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}
50
52 const std::size_t element_id)
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}
62
64 const double /*t*/, const ParameterLib::SpatialPosition& /*pos*/,
65 const double /*p*/, const double /*T*/, const double saturation) const
66{
67 return _relative_permeability_models[0]->getValue(saturation);
68}
69
71 const double /*t*/, const ParameterLib::SpatialPosition& /*pos*/,
72 const double /*p*/, const double /*T*/, const double saturation) const
73{
74 return _relative_permeability_models[1]->getValue(saturation);
75}
76
78 const int material_id, const double /*t*/,
79 const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
80 const double /*T*/, const double saturation) const
81{
82 return _capillary_pressure_models[material_id]->getCapillaryPressure(
83 saturation);
84}
85
87 const int material_id, const double /*t*/,
88 const ParameterLib::SpatialPosition& /*pos*/, const double /*p*/,
89 const double /*T*/, const double saturation) const
90{
91 return _capillary_pressure_models[material_id]->getdPcdS(saturation);
92}
93
95 double const t,
97 int const material_id,
98 double const pg,
99 double const X,
100 double const T,
101 double& Sw,
102 double& X_m,
103 double& dsw_dpg,
104 double& dsw_dX,
105 double& dxm_dpg,
106 double& dxm_dX)
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}
160 const int material_id, double const pl, double const X, double const T,
161 double const Sw, double const rho_h2_wet, ResidualVector& res)
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}
171
173 const int material_id, double const /*t*/,
174 ParameterLib::SpatialPosition const& /*x*/, double const pl,
175 double const /*X*/, double const T, JacobianMatrix& Jac, double const Sw,
176 double const rho_h2_wet)
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}
199
204 double const pg, double const Sw, double const rho_wet_h2)
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}
209
214 double /*PL*/, double X, double Sw, double rho_wet_h2, double rho_nonwet_h2,
215 double /*T*/)
216{
217 return X - (Sw * rho_wet_h2 + (1 - Sw) * rho_nonwet_h2);
218}
219
224 double pl, double S, double rho_wet_h2, double const T,
225 int current_material_id) const
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}
250 double const pl, const double /*X*/, const double S,
251 const double rho_wet_h2, double const T, int current_material_id) const
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}
276 double pl, double Sw, double rho_wet_h2, double dSwdX,
277 int current_material_id) const
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}
296 double pl, double Sw, double rho_wet_h2, double dSwdP,
297 int current_material_id) const
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}
312} // namespace TwoPhaseFlowWithPrho
313} // namespace ProcessLib
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the Mesh class.
Definition of the PiecewiseLinearInterpolation class.
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
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
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)
double getWetRelativePermeability(const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
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)
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)
static double calculateEquilibiumRhoWetLight(double const pg, double const Sw, double const rho_wet_h2)
double getCapillaryPressureDerivative(const int material_id, const double t, const ParameterLib::SpatialPosition &pos, const double p, const double T, const double saturation) const
Eigen::Matrix< double, jacobian_residual_size, jacobian_residual_size, Eigen::RowMajor > JacobianMatrix
double calculatedSwdP(double pl, double S, double rho_wet_h2, double const T, int current_material_id) const
double getCapillaryPressure(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)
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