From 4fd4fc70296c0c417cbd082280d738259e426708 Mon Sep 17 00:00:00 2001 From: Trine Mykkeltvedt Date: Mon, 13 Jun 2022 14:11:08 +0200 Subject: [PATCH] wip co2brine --- opm/material/constraintsolvers/ChiFlash.hpp | 665 +++--------------- .../chifluid/co2brinefluidsystem.hh | 196 ++++++ .../fluidsystems/chifluid/components.hh | 1 + 3 files changed, 282 insertions(+), 580 deletions(-) create mode 100644 opm/material/fluidsystems/chifluid/co2brinefluidsystem.hh diff --git a/opm/material/constraintsolvers/ChiFlash.hpp b/opm/material/constraintsolvers/ChiFlash.hpp index 67f70ecad..6e18e4bd0 100644 --- a/opm/material/constraintsolvers/ChiFlash.hpp +++ b/opm/material/constraintsolvers/ChiFlash.hpp @@ -146,17 +146,12 @@ public: fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0))); - // Rachford Rice equation to get initial L for composition solver - //L_scalar = solveRachfordRice_g_(K_scalar, z_scalar, verbosity); - - // Do a stability test to check if cell is single-phase (do for all cells the first time). bool isStable = false; if ( L <= 0 || L == 1 ) { if (verbosity >= 1) { std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl; } - //phaseStabilityTestMichelsen_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity); phaseStabilityTest_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity); } @@ -207,126 +202,17 @@ public: } updateDerivatives_(fluid_state_scalar, z, fluid_state, single); - - - - // fluid_state.setLvalue(L_scalar); - - // std::cout << " ------ SUMMARY AFTER DERIVATIVES ------ " << std::endl; - // std::cout << " L " << fluid_state.L() << std::endl; - // std::cout << " K " << fluid_state.K(0) << ", " << fluid_state.K(1) << ", " << fluid_state.K(2) << std::endl; - // std::cout << " x " << fluid_state.moleFraction(oilPhaseIdx, 0) << ", " << fluid_state.moleFraction(oilPhaseIdx, 1) << ", " << fluid_state.moleFraction(oilPhaseIdx, 2) << std::endl; - // std::cout << " y " << fluid_state.moleFraction(gasPhaseIdx, 0) << ", " << fluid_state.moleFraction(gasPhaseIdx, 1) << ", " << fluid_state.moleFraction(gasPhaseIdx, 2) << std::endl; - - - // Update phases - /* typename FluidSystem::template ParameterCache paramCache; - paramCache.updatePhase(fluid_state, oilPhaseIdx); - paramCache.updatePhase(fluid_state, gasPhaseIdx); */ - - /* // Calculate compressibility factor - const Scalar R = Opm::Constants::R; - InputEval Z_L = (paramCache.molarVolume(oilPhaseIdx) * fluid_state.pressure(oilPhaseIdx) ) / - (R * fluid_state.temperature(oilPhaseIdx)); - InputEval Z_V = (paramCache.molarVolume(gasPhaseIdx) * fluid_state.pressure(gasPhaseIdx) ) / - (R * fluid_state.temperature(gasPhaseIdx)); - - std::cout << " the type of InputEval here is " << Dune::className() << std::endl; - // Update saturation - InputEval So = L*Z_L/(L*Z_L+(1-L)*Z_V); - InputEval Sg = 1-So; - fluid_state.setSaturation(oilPhaseIdx, So); - fluid_state.setSaturation(gasPhaseIdx, Sg); - //Update L and K to the problem for the next flash - for (int compIdx = 0; compIdx < numComponents; ++compIdx){ - fluid_state.setKvalue(compIdx, K[compIdx]); + + //print summary after flash + if (verbosity >= 1) { + std::cout << " ------ SUMMARY AFTER FLASH ------ " << std::endl; + std::cout << " L " << fluid_state.L() << std::endl; + std::cout << " K " << fluid_state.K(0) << ", " << fluid_state.K(1) << ", " << fluid_state.K(2) << std::endl; + std::cout << " x " << fluid_state.moleFraction(oilPhaseIdx, 0) << ", " << fluid_state.moleFraction(oilPhaseIdx, 1) << ", " << fluid_state.moleFraction(oilPhaseIdx, 2) << std::endl; + std::cout << " y " << fluid_state.moleFraction(gasPhaseIdx, 0) << ", " << fluid_state.moleFraction(gasPhaseIdx, 1) << ", " << fluid_state.moleFraction(gasPhaseIdx, 2) << std::endl; } - fluid_state.setCompressFactor(oilPhaseIdx, Z_L); - fluid_state.setCompressFactor(gasPhaseIdx, Z_V); - fluid_state.setLvalue(L); */ - - - // Print saturation - - /* std::cout << " output the molefraction derivatives" << std::endl; - std::cout << " for vapor comp 1 " << std::endl; - fluid_state.moleFraction(gasPhaseIdx, 0).print(); - std::cout << std::endl << " for vapor comp 2 " << std::endl; - fluid_state.moleFraction(gasPhaseIdx, 1).print(); - std::cout << std::endl << " for vapor comp 3 " << std::endl; - fluid_state.moleFraction(gasPhaseIdx, 2).print(); - std::cout << std::endl; - std::cout << " for oil comp 1 " << std::endl; - fluid_state.moleFraction(oilPhaseIdx, 0).print(); - std::cout << std::endl << " for oil comp 2 " << std::endl; - fluid_state.moleFraction(oilPhaseIdx, 1).print(); - std::cout << std::endl << " for oil comp 3 " << std::endl; - fluid_state.moleFraction(oilPhaseIdx, 2).print(); - std::cout << std::endl; - std::cout << " for pressure " << std::endl; - fluid_state.pressure(0).print(); - std::cout<< std::endl; - fluid_state.pressure(1).print(); - std::cout<< std::endl; */ - - // Update densities - // fluid_state.setDensity(oilPhaseIdx, FluidSystem::density(fluid_state, paramCache, oilPhaseIdx)); - // fluid_state.setDensity(gasPhaseIdx, FluidSystem::density(fluid_state, paramCache, gasPhaseIdx)); - - // check the residuals of the equations - /* using NewtonVector = Dune::FieldVector; - NewtonVector newtonX; - NewtonVector newtonB; - for (int compIdx=0; compIdx= 5) { - std::cout << " mole fractions for oil " << std::endl; - for (int i = 0; i < numComponents; ++i) { - std::cout << " i " << i << " " << fluid_state.moleFraction(oilPhaseIdx, i) << std::endl; - } - std::cout << " mole fractions for gas " << std::endl; - for (int i = 0; i < numComponents; ++i) { - std::cout << " i " << i << " " << fluid_state.moleFraction(gasPhaseIdx, i) << std::endl; - } - std::cout << " K " << std::endl; - for (int i = 0; i < numComponents; ++i) { - std::cout << " i " << K[i] << std::endl; - } - std::cout << "Z_L = " << Z_L < - static void phaseStabilityTestMichelsen_(bool& stable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& z, int verbosity) - { - // Declarations - bool isTrivialL, isTrivialV; - ComponentVector x, y; - typename FlashFluidState::Scalar S_l, S_v; - ComponentVector K0 = K; - ComponentVector K1 = K; - - // Check for vapour instable phase - if (verbosity == 3 || verbosity == 4) { - std::cout << "Stability test for vapor phase:" << std::endl; - } - bool stable_vapour = false; - michelsenTest_(fluidState, z, y, K0,stable_vapour,/*isGas*/true, verbosity); - - bool stable_liquid = false; - michelsenTest_(fluidState, z, x, K1,stable_liquid,/*isGas*/false, verbosity); - - //bool stable = false; - stable = stable_liquid && stable_vapour; - - if (!stable) { - for (int compIdx = 0; compIdx= 1) { - std::cout << "Stability test done for - vapour - liquid - sum:" << stable_vapour << " - " << stable_liquid << " - " << stable < - static void michelsenTest_(const FlashFluidState& fluidState, const ComponentVector z, ComponentVector& xy_out, - ComponentVector& K, bool& stable, bool isGas, int verbosity) - { - using FlashEval = typename FlashFluidState::Scalar; - - using PengRobinsonMixture = typename Opm::PengRobinsonMixture; - - // Declarations - FlashFluidState fluidState_xy = fluidState; - FlashFluidState fluidState_zi = fluidState; - ComponentVector xy_loc; - ComponentVector R; - FlashEval S_loc = 0.0; - FlashEval xy_c = 0.0; - bool isTrivial; - bool isConverged; - - int phaseIdx = (isGas ? static_cast(gasPhaseIdx) : static_cast(oilPhaseIdx)); - int phaseIdx2 = (isGas ? static_cast(oilPhaseIdx) : static_cast(gasPhaseIdx)); - - // Setup output - if (verbosity >= 3 || verbosity >= 4) { - std::cout << std::setw(10) << "Iteration" << std::setw(16) << "K-Norm" << std::setw(16) << "R-Norm" << std::endl; - } - //mixture fugacity - for (int compIdx=0; compIdx paramCache_zi; - paramCache_zi.updatePhase(fluidState_zi, oilPhaseIdx); - - for (int compIdx=0; compIdx paramCache_xy; - paramCache_xy.updatePhase(fluidState_xy, phaseIdx); - - for (int compIdx=0; compIdx= 3) { - std::cout << std::setw(10) << i << std::setw(16) << K_norm << std::setw(16) << R_norm << std::endl; - } - - // Check convergence - isTrivial = (K_norm < 1e-5); - isConverged = (R_norm < 1e-10); - - bool ok = isTrivial || isConverged; - bool done = ok || i == maxIter; - - if (done && !ok) { - isTrivial = true; - throw std::runtime_error(" Stability test did not converge"); - //@warn "Stability test failed to converge in $maxiter iterations. Assuming stability." cond xy K_norm R_norm K - } - if (ok) { - stable = isTrivial || S_loc <= 1 + 1e-5; - return; - } - //todo: make sure that no mole fraction is smaller than 1e-8 ? - //todo: take care of water! - } - - throw std::runtime_error(" Stability test did not converge"); - }//end michelsen - template static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& globalComposition, int verbosity) @@ -847,7 +554,6 @@ template fluidState_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake); fluidState_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal); } - ComponentVector R; for (int compIdx=0; compIdx } } } - // for (unsigned i = 0; i < num_equations; ++i) { - // for (unsigned j = 0; j < num_primary_variables; ++j) { - // std::cout << " " << jac[i][j] ; - // } - // std::cout << std::endl; - // } - std::cout << std::endl; + if (!converged) { throw std::runtime_error(" Newton composition update did not converge within maxIterations"); } @@ -1325,7 +1025,6 @@ template (secondary_fluid_state, secondary_z, sec_jac, sec_res); } - // assembly the major matrix here // primary variables are x, y and L constexpr size_t primary_num_pv = numMisciblePhases * numMiscibleComponents + 1; @@ -1339,26 +1038,20 @@ template for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { primary_z[comp_idx] = Opm::getValue(z[comp_idx]); } - // TODO: x and y are not needed here - { - std::vector x(numComponents), y(numComponents); - PrimaryEval l; + for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx); - x[comp_idx] = x_ii;//PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx); - primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x[comp_idx]); + primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii); const unsigned idx = comp_idx + numComponents; const auto y_ii = PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx); - y[comp_idx] = y_ii;//;PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx); primary_fluid_state.setMoleFraction(gasPhaseIdx, comp_idx, y_ii); primary_fluid_state.setKvalue(comp_idx, y_ii / x_ii); } + PrimaryEval l; l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1); primary_fluid_state.setLvalue(l); - } primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx)); primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx)); - primary_fluid_state.setTemperature(fluid_state_scalar.temperature(0)); // TODO: is PrimaryFlashFluidState::Scalar> PrimaryEval here? @@ -1387,176 +1080,97 @@ template (primary_fluid_state, primary_z, pri_jac, pri_res); } - //corresponds to julias J_p (we miss d/dt, and have d/dL instead of d/dV) - // for (unsigned i =0; i < num_equations; ++i) { - // for (unsigned j = 0; j < primary_num_pv; ++j) { - // std::cout << " " << pri_jac[i][j]; - // } - // std::cout << std::endl; - // } - // std::cout << std::endl; - - //corresponds to julias J_s - // for (unsigned i = 0; i < num_equations; ++i) { - // for (unsigned j = 0; j < secondary_num_pv; ++j) { - // std::cout << " " << sec_jac[i][j] ; - // } - // std::cout << std::endl; - // } - // std::cout << std::endl; - SecondaryNewtonMatrix xx; pri_jac.solve(xx,sec_jac); - // for (unsigned i = 0; i < num_equations; ++i) { - // for (unsigned j = 0; j < secondary_num_pv; ++j) { - // std::cout << " " << xx[i][j] ; - // } - // std::cout << std::endl; - // } - using InputEval = typename FluidState::Scalar; using ComponentVectorMoleFraction = Dune::FieldVector; - //std::vector x(numComponents), y(numComponents); ComponentVectorMoleFraction x(numComponents), y(numComponents); - InputEval L_eval = L; - // TODO: then beginning from that point - - { - const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx); - const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx); - std::vector K(numComponents); - - // const double L = fluid_state_scalar.L(); - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - K[compIdx] = fluid_state_scalar.K(compIdx); - x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);//;z[compIdx] * 1. / (L + (1 - L) * K[compIdx]); - y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);//;x[compIdx] * K[compIdx]; - } - + // use the chainrule (and using partial instead of total + // derivatives, DF / Dp = dF / dp + dF / ds * ds/dp. + // where p is the primary variables and s the secondary variables. We then obtain + // ds / dp = -inv(dF / ds)*(DF / Dp) - // then we try to set the derivatives for x, y and K against P and x. - // p_l and p_v are the same here, in the future, there might be slightly more complicated scenarios when capillary - // pressure joins - { - constexpr size_t num_deri = numComponents; - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - std::vector deri(num_deri, 0.); - // derivatives from P - // for (unsigned idx = 0; idx < num_deri; ++idx) { - // probably can use some DUNE operator for vectors or matrics - for (unsigned idx = 0; idx < num_deri; ++idx) { - deri[idx] = - xx[compIdx][0] * p_l.derivative(idx); - } - // } + const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx); + const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx); + std::vector K(numComponents); - for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) { - const double pz = -xx[compIdx][cIdx + 1]; - const auto& zi = z[cIdx]; - for (unsigned idx = 0; idx < num_deri; ++idx) { - deri[idx] += pz * zi.derivative(idx); - } - } - for (unsigned idx = 0; idx < num_deri; ++idx) { - x[compIdx].setDerivative(idx, deri[idx]); - } - // handling y - for (unsigned idx = 0; idx < num_deri; ++idx) { - deri[idx] = - xx[compIdx + numComponents][0]* p_v.derivative(idx); - } - for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) { - const double pz = -xx[compIdx + numComponents][cIdx + 1]; - const auto& zi = z[cIdx]; - for (unsigned idx = 0; idx < num_deri; ++idx) { - deri[idx] += pz * zi.derivative(idx); - } - } - for (unsigned idx = 0; idx < num_deri; ++idx) { - y[compIdx].setDerivative(idx, deri[idx]); - } - } - // handling derivatives of L - std::vector deri(num_deri, 0.); - for (unsigned idx = 0; idx < num_deri; ++idx) { - deri[idx] = - xx[2*numComponents][0] * p_v.derivative(idx); - } - for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) { - const double pz = -xx[2*numComponents][cIdx + 1]; - const auto& zi = z[cIdx]; - for (unsigned idx = 0; idx < num_deri; ++idx) { - deri[idx] += pz * zi.derivative(idx); - } - } - for (unsigned idx = 0; idx < num_deri; ++idx) { - L_eval.setDerivative(idx, deri[idx]); - } - } + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + K[compIdx] = fluid_state_scalar.K(compIdx); + x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);//;z[compIdx] * 1. / (L + (1 - L) * K[compIdx]); + y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);//;x[compIdx] * K[compIdx]; } - // x, y og L_eval - - + + // then we try to set the derivatives for x, y and K against P and x. + // p_l and p_v are the same here, in the future, there might be slightly more complicated scenarios when capillary + // pressure joins + + constexpr size_t num_deri = numComponents; + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + std::vector deri(num_deri, 0.); + // derivatives from P + for (unsigned idx = 0; idx < num_deri; ++idx) { + deri[idx] = - xx[compIdx][0] * p_l.derivative(idx); + } + for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) { + const double pz = -xx[compIdx][cIdx + 1]; + const auto& zi = z[cIdx]; + for (unsigned idx = 0; idx < num_deri; ++idx) { + deri[idx] += pz * zi.derivative(idx); + } + } + + for (unsigned idx = 0; idx < num_deri; ++idx) { + x[compIdx].setDerivative(idx, deri[idx]); + } + + // handling y + for (unsigned idx = 0; idx < num_deri; ++idx) { + deri[idx] = - xx[compIdx + numComponents][0]* p_v.derivative(idx); + } + + for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) { + const double pz = -xx[compIdx + numComponents][cIdx + 1]; + const auto& zi = z[cIdx]; + for (unsigned idx = 0; idx < num_deri; ++idx) { + deri[idx] += pz * zi.derivative(idx); + } + } + + for (unsigned idx = 0; idx < num_deri; ++idx) { + y[compIdx].setDerivative(idx, deri[idx]); + } + + // handling derivatives of L + std::vector deriL(num_deri, 0.); + for (unsigned idx = 0; idx < num_deri; ++idx) { + deriL[idx] = - xx[2*numComponents][0] * p_v.derivative(idx); + } + + for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) { + const double pz = -xx[2*numComponents][cIdx + 1]; + const auto& zi = z[cIdx]; + for (unsigned idx = 0; idx < num_deri; ++idx) { + deriL[idx] += pz * zi.derivative(idx); + } + } + + for (unsigned idx = 0; idx < num_deri; ++idx) { + L_eval.setDerivative(idx, deriL[idx]); + } + // set up the mole fractions for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]); fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]); } - fluid_state.setLvalue(L_eval); - - } - - /* template - static void evalJacobian(const ComponentVector& globalComposition, - const Vector& x, - const Vector& y, - const Eval& l, - Vector& b, - Matrix& m) - { - // TODO: all the things are going through the FluidState, which makes it difficult to get the AD correct. - FluidState fluidState(fluidStateIn); - ComponentVector K; - for (int compIdx=0; compIdx; - ParamCache paramCache; - for (int phaseIdx=0; phaseIdx static void evalDefect_(DefectVector& b, @@ -1571,7 +1185,6 @@ template fluidState.setMoleFraction(oilPhaseIdx, compIdx, x[compIdx]); fluidState.setMoleFraction(gasPhaseIdx, compIdx, x[compIdx + numMiscibleComponents]); } - // Compute fugacities using ValueType = typename FluidState::Scalar; @@ -1602,7 +1215,7 @@ template // sum(x) - sum(y) = 0 b[numMiscibleComponents*numMisciblePhases] += -x[compIdx] + x[compIdx + numMiscibleComponents]; } - }//end valDefect + }//end evalDefect template static void evalJacobian_(DefectMatrix& A, @@ -1748,114 +1361,6 @@ template } // throw std::runtime_error("Successive substitution composition update did not converge within maxIterations"); } - - template - static void successiveSubstitutionCompositionNew_(ComponentVector& K, typename FlashFluidState::Scalar& L, FlashFluidState& fluidState, const ComponentVector& z, - const bool newton_afterwards, const int verbosity) - { - // std::cout << " Evaluation in successiveSubstitutionComposition_ is " << Dune::className(L) << std::endl; - // Determine max. iterations based on if it will be used as a standalone flash or as a pre-process to Newton (or other) method. - const int maxIterations = newton_afterwards ? 3 : 10; - - // Store cout format before manipulation - std::ios_base::fmtflags f(std::cout.flags()); - - // Print initial guess - if (verbosity >= 1) - std::cout << "Initial guess: K = [" << K << "] and L = " << L << std::endl; - - if (verbosity == 2 || verbosity == 4) { - // Print header - int fugWidth = (numComponents * 12)/2; - int convWidth = fugWidth + 7; - std::cout << std::setw(10) << "Iteration" << std::setw(fugWidth) << "fL/fV" << std::setw(convWidth) << "norm2(fL/fv-1)" << std::endl; - } - // - // Successive substitution loop - // - for (int i=0; i < maxIterations; ++i){ - // Compute (normalized) liquid and vapor mole fractions - computeLiquidVapor_(fluidState, L, K, z); - - // Calculate fugacity coefficient - using ParamCache = typename FluidSystem::template ParameterCache; - ParamCache paramCache; - for (int phaseIdx=0; phaseIdx= 1) { - std::cout << "Solution converged to the following result :" << std::endl; - std::cout << "x = ["; - for (int compIdx=0; compIdx +#include + +#include "ChiParameterCache.hpp" +#include "LBCviscosity.hpp" + +namespace Opm { +/*! + * \ingroup FluidSystem + * + * \brief A two phase two component system, co2 brine + */ + + template + class Co2BrineFluidSystem + : public Opm::BaseFluidSystem > { + public: + // TODO: I do not think these should be constant in fluidsystem, will try to make it non-constant later + static constexpr int numPhases=2; + static constexpr int numComponents = 3; + // TODO: phase location should be more general + static constexpr int oilPhaseIdx = 0; + static constexpr int gasPhaseIdx = 1; + + static constexpr int Comp0Idx = 0; + static constexpr int Comp1Idx = 1; + static constexpr int Comp2Idx = 2; + + // TODO: needs to be more general + using Comp0 = Opm::JuliaCO2; + using Comp1 = Opm::ChiwomsBrine; + using Comp2 = Opm::JuliaC10; + + template + using ParameterCache = Opm::ChiParameterCache>; + using LBCviscosity = typename Opm::LBCviscosity>; + using PengRobinsonMixture = typename Opm::PengRobinsonMixture>; + + /*! + * \brief The acentric factor of a component []. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar acentricFactor(unsigned compIdx) + { + switch (compIdx) { + case Comp0Idx: return Comp0::acentricFactor(); + case Comp1Idx: return Comp1::acentricFactor(); + case Comp2Idx: return Comp2::acentricFactor(); + default: throw std::runtime_error("Illegal component index for acentricFactor"); + } + } + /*! + * \brief Critical temperature of a component [K]. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar criticalTemperature(unsigned compIdx) + { + switch (compIdx) { + case Comp0Idx: return Comp0::criticalTemperature(); + case Comp1Idx: return Comp1::criticalTemperature(); + case Comp2Idx: return Comp2::criticalTemperature(); + default: throw std::runtime_error("Illegal component index for criticalTemperature"); + } + } + /*! + * \brief Critical pressure of a component [Pa]. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar criticalPressure(unsigned compIdx) { + switch (compIdx) { + case Comp0Idx: return Comp0::criticalPressure(); + case Comp1Idx: return Comp1::criticalPressure(); + case Comp2Idx: return Comp2::criticalPressure(); + default: throw std::runtime_error("Illegal component index for criticalPressure"); + } + } + /*! + * \brief Critical volume of a component [m3]. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar criticalVolume(unsigned compIdx) + { + switch (compIdx) { + case Comp0Idx: return Comp0::criticalVolume(); + case Comp1Idx: return Comp1::criticalVolume(); + case Comp2Idx: return Comp2::criticalVolume(); + default: throw std::runtime_error("Illegal component index for criticalVolume"); + } + } + + //! \copydoc BaseFluidSystem::molarMass + static Scalar molarMass(unsigned compIdx) + { + switch (compIdx) { + case Comp0Idx: return Comp0::molarMass(); + case Comp1Idx: return Comp1::molarMass(); + case Comp2Idx: return Comp2::molarMass(); + default: throw std::runtime_error("Illegal component index for molarMass"); + } + } + + /*! + * \brief Returns the interaction coefficient for two components. + *. + */ + static Scalar interactionCoefficient(unsigned /*comp1Idx*/, unsigned /*comp2Idx*/) + { + return 0.0; + } + + //! \copydoc BaseFluidSystem::phaseName + static const char* phaseName(unsigned phaseIdx) + { + static const char* name[] = {"o", // oleic phase + "g"}; // gas phase + + assert(0 <= phaseIdx && phaseIdx < 2); + return name[phaseIdx]; + } + + //! \copydoc BaseFluidSystem::componentName + static const char* componentName(unsigned compIdx) + { + static const char* name[] = { + Comp0::name(), + Comp1::name(), + Comp2::name(), + }; + + assert(0 <= compIdx && compIdx < 3); + return name[compIdx]; + } + + /*! + * \copydoc BaseFluidSystem::density + */ + template + static LhsEval density(const FluidState& fluidState, + const ParameterCache& paramCache, + unsigned phaseIdx) + { + + LhsEval dens; + if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) { + // paramCache.updatePhase(fluidState, phaseIdx); + dens = fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx); + } + return dens; + + } + + //! \copydoc BaseFluidSystem::viscosity + template + static LhsEval viscosity(const FluidState& fluidState, + const ParameterCache& paramCache, + unsigned phaseIdx) + { + // Use LBC method to calculate viscosity + // LhsEval mu = LBCviscosity::LBCmod(fluidState, paramCache, phaseIdx); + // LhsEval mu = LBCviscosity::LBC(fluidState, paramCache, phaseIdx); + LhsEval mu; + mu = LBCviscosity::LBCmod(fluidState, paramCache, phaseIdx); + + // LhsEval mu = LBCviscosity::LBCJulia(fluidState, paramCache, phaseIdx); + return mu; + + } + + //! \copydoc BaseFluidSystem::fugacityCoefficient + template + static LhsEval fugacityCoefficient(const FluidState& fluidState, + const ParameterCache& paramCache, + unsigned phaseIdx, + unsigned compIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + assert(0 <= compIdx && compIdx < numComponents); + + // TODO: here the derivatives for the phi are dropped. Should we keep the derivatives against the pressure + // and temperature? + LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); + //Scalar phi = Opm::getValue( + // PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx)); + return phi; + } + + }; +} +#endif //OPM_CO2BRINEFLUIDSYSTEM_HH \ No newline at end of file diff --git a/opm/material/fluidsystems/chifluid/components.hh b/opm/material/fluidsystems/chifluid/components.hh index 848829edd..91a6f992a 100644 --- a/opm/material/fluidsystems/chifluid/components.hh +++ b/opm/material/fluidsystems/chifluid/components.hh @@ -225,6 +225,7 @@ public: // Critical volume [m3/kmol] static Scalar criticalVolume() {return 9.412e-5; } + // OLD :static Scalar criticalVolume() {return 9.4118e-2; } }; template