From 532c6c5652d48daf2e3e1b15dc149a0b22004306 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 15 Nov 2010 11:36:26 +0100 Subject: [PATCH] Work in progress on ensuring units are done correctly. --- dune/porsol/blackoil/fluid/BlackoilPVT.cpp | 21 +- .../blackoil/fluid/FluidStateBlackoil.hpp | 10 +- .../blackoil/fluid/FluidSystemBlackoil.hpp | 345 +++--------------- .../blackoil/fluid/MiscibilityWater.hpp | 26 +- 4 files changed, 87 insertions(+), 315 deletions(-) diff --git a/dune/porsol/blackoil/fluid/BlackoilPVT.cpp b/dune/porsol/blackoil/fluid/BlackoilPVT.cpp index f7c32005..4df2ced7 100644 --- a/dune/porsol/blackoil/fluid/BlackoilPVT.cpp +++ b/dune/porsol/blackoil/fluid/BlackoilPVT.cpp @@ -38,30 +38,25 @@ namespace Opm typedef std::vector > > table_t; region_number_ = 0; - - - - - - // Surface densities. Accounting for different orders in eclipse and our code. if (parser.hasField("DENSITY")) { const int region_number = 0; enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; - const std::vector > d_tmp = - parser.getDENSITY().densities_; - densities_[Aqua] = d_tmp[region_number][ECL_water]; - densities_[Liquid] = d_tmp[region_number][ECL_oil]; - densities_[Vapour] = d_tmp[region_number][ECL_gas]; + const std::vector& d = parser.getDENSITY().densities_[region_number]; + const double du = parser.units().density; + using namespace Dune::unit; + densities_[Aqua] = convert::from(d[ECL_water], du); + densities_[Vapour] = convert::from(d[ECL_gas], du); + densities_[Liquid] = convert::from(d[ECL_oil], du); } else { THROW("Input is missing DENSITY\n"); } // Water PVT if (parser.hasField("PVTW")) { - water_props_.reset(new MiscibilityWater(parser.getPVTW().pvtw_)); + water_props_.reset(new MiscibilityWater(parser.getPVTW().pvtw_, parser.units())); } else { - water_props_.reset(new MiscibilityWater(3e-4)); // Default is 0.3 cP. + water_props_.reset(new MiscibilityWater(0.5*Dune::prefix::centi*Dune::unit::Poise)); // Eclipse 100 default } // Oil PVT diff --git a/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp b/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp index 5c36e69a..bc7eb53d 100644 --- a/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp +++ b/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp @@ -233,13 +233,11 @@ public: { return temperature_; }; public: - Scalar fugacity_[numComponents]; - Scalar moleFrac_[numPhases][numComponents]; - Scalar phaseConcentration_[numPhases]; - Scalar meanMolarMass_[numPhases]; - Scalar phasePressure_[numPhases]; - Scalar saturation_[numPhases]; Scalar temperature_; + Scalar surface_volume_[numComponents]; + Scalar phase_pressure_[numPhases]; + Scalar phase_volume_[numPhases]; + Scalar saturation_[numPhases]; }; } // end namespace Opm diff --git a/dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp b/dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp index 58d3331f..6077bc1c 100644 --- a/dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp +++ b/dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp @@ -1,3 +1,17 @@ +/***************************************************************************** + * Copyright (C) 2009 by Andreas Lauser + * Institute of Hydraulic Engineering * + * University of Stuttgart, Germany * + * email: .@iws.uni-stuttgart.de * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version, as long as this copyright notice * + * is included in its original form. * + * * + * This program is distributed WITHOUT ANY WARRANTY. * + *****************************************************************************/ /* Copyright 2010 SINTEF ICT, Applied Mathematics. @@ -28,11 +42,11 @@ namespace Opm { // Forward declaration needed by associated parameters classes. -template -class FluidMatrixInteractionBlackoil; +template +class FluidSystemBlackoil; - -class FluidSystemBlackoilParametersNonmiscible +// Parameters for the black oil fluid system. +class FluidSystemBlackoilParameters { public: void init(const Dune::EclipseGridParser& ep) @@ -40,6 +54,8 @@ public: pvt_.init(ep); } private: + template + friend class FluidSystemBlackoil; BlackoilPVT pvt_; }; @@ -47,7 +63,7 @@ private: /*! * \brief A black oil fluid system. */ -template +template class FluidSystemBlackoil { public: @@ -60,6 +76,7 @@ public: enum ComponentIndex { Water = 0, Gas = 1, Oil = 2 }; enum PhaseIndex { Aqua = 0, Vapour = 1, Liquid = 2 }; + typedef BlackoilPVT::surfvol_t FluidVec; /*! * \brief Initialize system from input. @@ -96,15 +113,15 @@ public: } -#if 0 /*! * \brief Return the molar mass of a component [kg/mol]. */ static Scalar molarMass(int compIdx) { - + throw std::logic_error("molarMass() not implemented."); } + /*! * \brief Given a phase's composition, temperature, pressure, and * the partial pressures of all components, return its @@ -116,31 +133,7 @@ public: Scalar pressure, const FluidState &fluidState) { - if (phaseIdx == lPhaseIdx) { - // See: Ochs 2008 - // \todo: proper citation - Scalar rholH2O = H2O::liquidDensity(temperature, pressure); - Scalar clH2O = rholH2O/H2O::molarMass(); - - // this assumes each nitrogen molecule displaces exactly one - // water molecule in the liquid - return - clH2O*(H2O::molarMass()*fluidState.moleFrac(lPhaseIdx, H2OIdx) - + - N2::molarMass()*fluidState.moleFrac(lPhaseIdx, N2Idx)); - } - else if (phaseIdx == gPhaseIdx) { - Scalar fugH2O = - fluidState.moleFrac(gPhaseIdx, H2OIdx) * - fluidState.phasePressure(gPhaseIdx); - Scalar fugN2 = - fluidState.moleFrac(gPhaseIdx, N2Idx) * - fluidState.phasePressure(gPhaseIdx); - return - H2O::gasDensity(temperature, fugH2O) + - N2::gasDensity(temperature, fugN2); - } - DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); + throw std::logic_error("phaseDensity() not implemented."); } /*! @@ -153,151 +146,41 @@ public: Scalar pressure, const FluidState &fluidState) { - if (phaseIdx == lPhaseIdx) { - // assume pure water for the liquid phase - // TODO: viscosity of mixture - return H2O::liquidViscosity(temperature, - pressure); - } - else { - return N2::gasViscosity(temperature, - pressure); - /* Wilke method. See: - * - * S.O.Ochs: "Development of a multiphase multicomponent - * model for PEMFC - Technical report: IRTG-NUPUS", - * University of Stuttgart, 2008 - * - * and: - * - * See: R. Reid, et al.: The Properties of Gases and Liquids, 4th - * edition, McGraw-Hill, 1987, 407-410 - */ - Scalar muResult = 0; - const Scalar mu[numComponents] = { - H2O::gasViscosity(temperature, - H2O::vaporPressure(temperature)), - N2::gasViscosity(temperature, - pressure) - }; - // molar masses - const Scalar M[numComponents] = { - H2O::molarMass(), - N2::molarMass() - }; - - for (int i = 0; i < numComponents; ++i) { - Scalar divisor = 0; - for (int j = 0; j < numComponents; ++j) { - Scalar phiIJ = 1 + sqrt(mu[i]/mu[j] * - pow(M[i]/M[j], 1/4.0)); - phiIJ *= phiIJ; - phiIJ /= sqrt(8*(1 + M[i]/M[j])); - divisor += fluidState.moleFrac(phaseIdx, j)*phiIJ; - } - muResult += fluidState.moleFrac(phaseIdx, i)*mu[i] / divisor; - } - - return muResult; - } + throw std::logic_error("phaseViscosity() not implemented."); } /*! - * \brief Assuming the composition of a single phase and the - * pressure of all phases is known or that all phases are - * present, compute the thermodynamic equilibrium from the - * temperature and phase pressures. If the known phase - * index - * + * \brief Assuming the surface volumes and the pressure of all + * phases are known, compute the phase volumes and + * saturations. */ template - static void computeEquilibrium(FluidState &fluidState, - int knownPhaseIdx = -1) + static void computeEquilibrium(FluidState& fluid_state) { - const Scalar T = fluidState.temperature(); - const Scalar pg = fluidState.phasePressure(gPhaseIdx); - const Scalar pl = fluidState.phasePressure(lPhaseIdx); + // Get B and R factors. + const double* p = fluid_state.phase_pressure_; + const double* z = fluid_state.surface_volume_; + double B[3]; + B[Aqua] = params().pvt_.B(p[Aqua], z, Aqua); + B[Vapour] = params().pvt_.B(p[Vapour], z, Vapour); + B[Liquid] = params().pvt_.B(p[Liquid], z, Liquid); + double R[3]; // Only using 2 of them, though. + R[Vapour] = params().pvt_.B(p[Vapour], z, Vapour); + R[Liquid] = params().pvt_.B(p[Liquid], z, Liquid); - const Scalar betaH2O = H2O::vaporPressure(T); - const Scalar betaN2 = BinaryCoeff::H2O_N2::henry(T); + // Update phase volumes. + double* u = fluid_state.phase_volume_; + double detR = 1.0 - R[Vapour]*R[Liquid]; + u[Aqua] = B[Aqua]*z[Water]; + u[Vapour] = B[Vapour]*(z[Gas] - R[Liquid]*z[Oil])/detR; + u[Liquid] = B[Aqua]*(z[Oil] - R[Vapour]*z[Gas])/detR; - if (knownPhaseIdx < 0) - { - // we only have all phase pressures and temperature and - // know that all phases are present - Scalar xlH2O = (pg - betaN2)/(betaH2O - betaN2); - Scalar xlN2 = 1 - xlH2O; - Scalar rhol = liquidPhaseDensity_(T, pl, xlH2O, xlN2); - - Scalar cgH2O = H2O::gasDensity(T, betaH2O*xlH2O)/H2O::molarMass(); - Scalar cgN2 = N2::gasDensity(T, betaN2*xlN2)/N2::molarMass(); - - Scalar xgH2O = cgH2O/(cgH2O + cgN2); - Scalar xgN2 = cgN2/(cgH2O + cgN2); - - // set the liquid phase composition - SettablePhase liquid; - liquid.moleFrac_[H2OIdx] = xlH2O; - liquid.moleFrac_[N2Idx] = xlN2; - liquid.pressure_ = pl; - liquid.density_ = rhol; - liquid.xToX(); // compute mass fractions from mole fractions - fluidState.assignPhase(lPhaseIdx, liquid); - - // set the gas phase composition - SettablePhase gas; - gas.moleFrac_[H2OIdx] = xgH2O; - gas.moleFrac_[N2Idx] = xgN2; - gas.pressure_ = pg; - gas.density_ = cgH2O*H2O::molarMass() + cgN2*N2::molarMass(); - gas.xToX(); // compute mass fractions from mole fractions - fluidState.assignPhase(gPhaseIdx, gas); - } - else if (knownPhaseIdx == lPhaseIdx) { - // the composition of the liquid phase is given - - // retrieve the known mole fractions from the fluid state - Scalar xlH2O = fluidState.moleFrac(lPhaseIdx, H2OIdx); - Scalar xlN2 = fluidState.moleFrac(lPhaseIdx, N2Idx); - - // calculate the component contentrations in the gas phase - Scalar pH2O = betaH2O*xlH2O; // fugacity of water - Scalar pN2 = betaN2*xlN2; // fugacity of nitrogen - Scalar cgH2O = H2O::gasDensity(T, pH2O)/H2O::molarMass(); - Scalar cgN2 = N2::gasDensity(T, pN2)/N2::molarMass(); - - // convert concentrations to mole fractions - Scalar xgH2O = cgH2O/(cgH2O + cgN2) * (pH2O + pN2)/pg; - Scalar xgN2 = cgN2/(cgH2O + cgN2) * (pH2O + pN2)/pg; - - // set gas phase composition - SettablePhase gas; - gas.moleFrac_[H2OIdx] = xgH2O; - gas.moleFrac_[N2Idx] = xgN2; - gas.pressure_ = pg; - gas.density_ = cgH2O*H2O::molarMass() + cgN2*N2::molarMass(); - gas.xToX(); // update mass fractions from mole fractions - fluidState.assignPhase(gPhaseIdx, gas); - } - else if (knownPhaseIdx == gPhaseIdx) { - // the composition of the gas phase is given - - Scalar xgH2O = fluidState.moleFrac(gPhaseIdx, H2OIdx); - Scalar xgN2 = fluidState.moleFrac(gPhaseIdx, N2Idx); - Scalar pgH2O = pg*xgH2O; - Scalar pgN2 = pg*xgN2; - - Scalar xlH2O = pgH2O/betaH2O; - Scalar xlN2 = pgN2/betaN2; - - SettablePhase liquid; - liquid.moleFrac_[H2OIdx] = xlH2O; - liquid.moleFrac_[N2Idx] = xlN2; - liquid.pressure_ = pl; - liquid.density_ = liquidPhaseDensity_(T, pl, xlH2O, xlN2); - liquid.xToX(); // update mass fractions from mole fractions - fluidState.assignPhase(lPhaseIdx, liquid); + // Update saturations. + double sumu = u[Aqua] + u[Vapour] + u[Liquid]; + double* s = fluid_state.saturation_; + for (int i = 0; i < 3; ++i) { + s[i] = u[i]/sumu; } } @@ -326,28 +209,8 @@ public: Scalar pressure, const FluidState &state) { - if (phaseIdx == gPhaseIdx) { - return pressure; - Scalar fugH2O = std::max(1e-3, state.fugacity(H2OIdx)); - Scalar fugN2 = std::max(1e-3, state.fugacity(N2Idx)); - Scalar cH2O = H2O::gasDensity(temperature, fugH2O) / H2O::molarMass(); - Scalar cN2 = N2::gasDensity(temperature, fugN2) / N2::molarMass(); - - Scalar alpha = (fugH2O + fugN2)/pressure; - - if (compIdx == H2OIdx) - return fugH2O/(alpha*cH2O/(cH2O + cN2)); - else if (compIdx == N2Idx) - return fugN2/(alpha*cN2/(cH2O + cN2)); - - DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); - } - - switch (compIdx) { - case H2OIdx: return H2O::vaporPressure(temperature); - case N2Idx: return BinaryCoeff::H2O_N2::henry(temperature); - }; - DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); + throw std::logic_error("activityCoeff() not implemented."); + return 0.0; } /*! @@ -363,49 +226,8 @@ public: Scalar pressure, const FluidState &fluidState) { - if (compIIdx > compJIdx) - std::swap(compIIdx, compJIdx); - -#ifndef NDEBUG - if (compIIdx == compJIdx || - phaseIdx > numPhases - 1 || - compJIdx > numComponents - 1) - { - DUNE_THROW(Dune::InvalidStateException, - "Binary diffusion coefficient of components " - << compIIdx << " and " << compJIdx - << " in phase " << phaseIdx << " is undefined!\n"); - } -#endif - - - switch (phaseIdx) { - case lPhaseIdx: - switch (compIIdx) { - case H2OIdx: - switch (compJIdx) { - case N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(temperature, - pressure); - } - default: - DUNE_THROW(Dune::InvalidStateException, - "Binary diffusion coefficients of trace " - "substances in liquid phase is undefined!\n"); - } - case gPhaseIdx: - switch (compIIdx) { - case H2OIdx: - switch (compJIdx) { - case N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(temperature, - pressure); - } - } - } - - DUNE_THROW(Dune::InvalidStateException, - "Binary diffusion coefficient of components " - << compIIdx << " and " << compJIdx - << " in phase " << phaseIdx << " is undefined!\n"); + throw std::logic_error("diffCoeff() not implemented."); + return 0.0; }; /*! @@ -418,35 +240,8 @@ public: Scalar pressure, const FluidState &fluidState) { - if (phaseIdx == lPhaseIdx) { - Scalar cN2 = fluidState.concentration(lPhaseIdx, N2Idx); - Scalar pN2 = N2::gasPressure(temperature, cN2*N2::molarMass()); - - // TODO: correct way to deal with the solutes??? - return - fluidState.massFrac(lPhaseIdx, H2OIdx)* - H2O::liquidEnthalpy(temperature, pressure) - + - fluidState.massFrac(lPhaseIdx, N2Idx)* - N2::gasEnthalpy(temperature, pN2); - } - else { - Scalar cH2O = fluidState.concentration(gPhaseIdx, H2OIdx); - Scalar cN2 = fluidState.concentration(gPhaseIdx, N2Idx); - - Scalar pH2O = H2O::gasPressure(temperature, cH2O*H2O::molarMass()); - Scalar pN2 = N2::gasPressure(temperature, cN2*N2::molarMass()); - - Scalar result = 0; - result += - H2O::gasEnthalpy(temperature, pH2O) * - fluidState.massFrac(gPhaseIdx, H2OIdx); - result += - N2::gasEnthalpy(temperature, pN2) * - fluidState.massFrac(gPhaseIdx, N2Idx); - - return result; - } + throw std::logic_error("phaseEnthalpy() not implemented."); + return 0.0; } /*! @@ -459,33 +254,11 @@ public: Scalar pressure, const FluidState &fluidState) { - return - phaseEnthalpy(phaseIdx, temperature, pressure, fluidState) - - pressure/phaseDensity(phaseIdx, temperature, pressure, fluidState); + throw std::logic_error("phaseInternalEnergy() not implemented."); + return 0.0; } private: - static Scalar liquidPhaseDensity_(Scalar T, Scalar pl, Scalar xlH2O, Scalar xlN2) - { - // See: Ochs 2008 - // \todo: proper citation - Scalar rholH2O = H2O::liquidDensity(T, pl); - Scalar clH2O = rholH2O/H2O::molarMass(); - - // this assumes each nitrogen molecule displaces exactly one - // water molecule in the liquid - return - clH2O*(xlH2O*H2O::molarMass() - + - xlN2*N2::molarMass()); - } - - static Scalar gasPhaseDensity_(Scalar T, Scalar pg, Scalar xgH2O, Scalar xgN2) - { - return H2O::gasDensity(T, pg*xgH2O) + N2::gasDensity(T, pg*xgN2); - }; - -#endif static Params& params() { diff --git a/dune/porsol/blackoil/fluid/MiscibilityWater.hpp b/dune/porsol/blackoil/fluid/MiscibilityWater.hpp index e532e211..66708868 100644 --- a/dune/porsol/blackoil/fluid/MiscibilityWater.hpp +++ b/dune/porsol/blackoil/fluid/MiscibilityWater.hpp @@ -36,6 +36,7 @@ #include "MiscibilityProps.hpp" #include +#include // Forward declaration. class PVTW; @@ -46,24 +47,26 @@ namespace Opm { public: typedef std::vector > table_t; - MiscibilityWater(const table_t& pvtw) + MiscibilityWater(const table_t& pvtw, const Dune::EclipseUnits& units) { const int region_number = 0; if (pvtw.size() != 1) { THROW("More than one PVD-region"); } - double b = pvtw[region_number][1]; - double comp = pvtw[region_number][2]; - double visc = pvtw[region_number][3]; - const double VISCOSITY_UNIT = 1e-3; - if (b == 1.0 && comp == 0) { - viscosity_ = visc*VISCOSITY_UNIT; - } else { - THROW("Compressible water not implemented."); + using namespace Dune::unit; + ref_press_ = convert::from(pvtw[region_number][0], units.pressure); + ref_B_ = convert::from(pvtw[region_number][1], units.liqvol_r/units.liqvol_s); + comp_ = convert::from(pvtw[region_number][2], units.compressibility); + viscosity_ = convert::from(pvtw[region_number][3], units.viscosity); + if (pvtw[region_number].size() > 4 && pvtw[region_number][4] != 0.0) { + THROW("MiscibilityWater does not support 'viscosibility'."); } } MiscibilityWater(double visc) - : viscosity_(visc) + : ref_press_(0.0), + ref_B_(1.0), + comp_(0.0), + viscosity_(visc) { } virtual ~MiscibilityWater() @@ -91,6 +94,9 @@ namespace Opm return 0.0; } private: + double ref_press_; + double ref_B_; + double comp_; double viscosity_; };