Commit 6cd659ae authored by Vishal Jambhekar's avatar Vishal Jambhekar
Browse files

New test case for 2pncmin model added to dumux/test. For this following...

New test case for 2pncmin model added to dumux/test. For this following additional files are needed (Reviewde by Johannes)

- For variation of salinity in brine the needed constitutive relations are implemented in components brine_VerSalinity.hh and nacl.hh
- Additional constraint solvers (compositionfromfugacities2pncmin and computefromreferencephase2pncmin) are needed, specially for the case nPhaseOnly
- Additional constraint solvers (computefromreferencephase2pnc and miscible2pnccomposition) are needed for 2pnc cases in general.
- Fluidsystem brineairfluidsystem is updated
- Binerycoefficient file for brine_air case is included. 

- A test case for 2pncmin model is included for flusiong of precipitated salt in gas reserviors. The includes corresponding problem file, spatial parameters file, input file and *.cc file.
- For this the CMakeLists.txt is also updated.
- A reference solution 2pncmin test case "injectionbox2pncmin-reference" is included to test/references folder

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@15196 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 8326018c
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* 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. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \ingroup Components
*
* \brief A class for the brine fluid properties,.
*/
#ifndef DUMUX_BRINE_VARSALINITY_HH
#define DUMUX_BRINE_VARSALINITY_HH
#include <dumux/material/components/component.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/nacl.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
#include <cmath>
namespace Dumux
{
/*!
*
* \ingroup Components
*
* \brief A class for the brine fluid properties.
*
* \tparam Scalar The type used for scalar values
* \tparam H2O Static polymorphism: the Brine class can access all properties of the H2O class
*/
template <class Scalar,
class H2O_Tabulated = Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar>>>
class BrineVarSalinity : public Component<Scalar, BrineVarSalinity<Scalar, H2O_Tabulated> >
{
public:
typedef Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar>> H2O;
/*!
* \brief A human readable name for the brine.
*/
static const char *name()
{ return "Brine"; }
/*!
* \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of brine.
*
* This assumes that the salt is pure NaCl.
*/
static Scalar molarMass(Scalar salinity)
{
const Scalar M1 = H2O::molarMass();
const Scalar M2 = NaCl<Scalar>::molarMass(); // molar mass of NaCl [kg/mol]
const Scalar X2 = salinity; // mass fraction of salt in brine
return M1*M2/(M2 + X2*(M1 - M2));
};
/*!
* \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of brine. Here, it is assumed to be equal to that of pure water.
*/
static Scalar criticalTemperature()
{ return H2O::criticalTemperature(); /* [K] */ }
/*!
* \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of brine. Here, it is assumed to be equal to that of pure water.
*/
static Scalar criticalPressure()
{ return H2O::criticalPressure(); /* [N/m^2] */ }
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at brine's triple point. Here, it is assumed to be equal to that of pure water.
*/
static Scalar tripleTemperature()
{ return H2O::tripleTemperature(); /* [K] */ }
/*!
* \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at brine's triple point. Here, it is assumed to be equal to that of pure water.
*/
static Scalar triplePressure()
{ return H2O::triplePressure(); /* [N/m^2] */ }
/*!
* \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of pure brine
* at a given temperature. Here, it is assumed to be equal to that of pure water.
*
* \param T temperature of component in \f$\mathrm{[K]}\f$
*/
static Scalar vaporPressure(Scalar T)
{ return H2O::vaporPressure(T); /* [N/m^2] */ }
/*!
* \brief Specific enthalpy of gaseous brine \f$\mathrm{[J/kg]}\f$.
* Only water volatile and salt is suppose to stay in the liquid phase.
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasEnthalpy(Scalar temperature,
Scalar pressure)
{ return H2O::gasEnthalpy(temperature, pressure); /* [J/kg] */ }
/*!
* \brief Specific enthalpy of liquid brine \f$\mathrm{[J/kg]}\f$.
*
* \param T temperature of component in \f$\mathrm{[K]}\f$
* \param p pressure of component in \f$\mathrm{[Pa]}\f$
*
* Equations given in: - Palliser & McKibbin 1997
* - Michaelides 1981
* - Daubert & Danner 1989
*/
static const Scalar liquidEnthalpy(Scalar T,
Scalar p, Scalar salinity)
{
/*Numerical coefficents from PALLISER*/
static const Scalar f[] = {
2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10
};
/*Numerical coefficents from MICHAELIDES for the enthalpy of brine*/
static const Scalar a[4][3] = {
{ -9633.6, -4080.0, +286.49 },
{ +166.58, +68.577, -4.6856 },
{ -0.90963, -0.36524, +0.249667E-1 },
{ +0.17965E-2, +0.71924E-3, -0.4900E-4 }
};
Scalar theta, h_NaCl;
Scalar m, h_ls, h_ls1, d_h;
Scalar S_lSAT, delta_h;
int i, j;
Scalar hw;
theta = T - 273.15;
Scalar S = salinity;
S_lSAT = f[0] + f[1]*theta + f[2]*pow(theta,2) + f[3]*pow(theta,3);
/*Regularization*/
if (S>S_lSAT) {
S = S_lSAT;
}
hw = H2O::liquidEnthalpy(T, p)/1E3; /* kJ/kg */
/*DAUBERT and DANNER*/
/*U=*/h_NaCl = (3.6710E4*T + 0.5*(6.2770E1)*T*T - ((6.6670E-2)/3)*T*T*T
+((2.8000E-5)/4)*pow(T,4))/(58.44E3)- 2.045698e+02; /* kJ/kg */
m = (1E3/58.44)*(S/(1-S));
i = 0;
j = 0;
d_h = 0;
for (i = 0; i<=3; i++) {
for (j=0; j<=2; j++) {
d_h = d_h + a[i][j] * pow(theta, i) * pow(m, j);
}
}
delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
/* Enthalpy of brine */
h_ls1 =(1-S)*hw + S*h_NaCl + S*delta_h; /* kJ/kg */
h_ls = h_ls1*1E3; /*J/kg*/
return (h_ls);
}
/*!
* \brief Specific internal energy of steam \f$\mathrm{[J/kg]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasInternalEnergy(Scalar temperature,
Scalar pressure)
{
return H2O::gasInternalEnergy(temperature, pressure);
}
/*!
* \brief Specific internal energy of liquid brine \f$\mathrm{[J/kg]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidInternalEnergy(Scalar temperature,
Scalar pressure, Scalar salinity)
{
return
liquidEnthalpy(temperature, pressure) -
pressure/liquidDensity(temperature, pressure);
}
/*!
* \brief The density of steam at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
{ return H2O::gasDensity(temperature, pressure); }
/*!
* \brief The density of pure brine at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
* Equations given in: - Batzle & Wang (1992)
* - cited by: Adams & Bachu in Geofluids (2002) 2, 257-271
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure, Scalar salinity)
{
Scalar TempC = temperature - 273.15;
Scalar pMPa = pressure/1.0E6;
salinity = std::abs(salinity);
Scalar rhow = H2O::liquidDensity(temperature, pressure);
Scalar density = rhow +
1000*salinity*(
0.668 +
0.44*salinity +
1.0E-6*(
300*pMPa -
2400*pMPa*salinity +
TempC*(
80.0 -
3*TempC -
3300*salinity -
13*pMPa +
47*pMPa*salinity)));
assert(density > 0.0);
return density;
}
/*!
* \brief The pressure of steam in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density denstiy of component in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
{ return H2O::gasPressure(temperature, density); }
/*!
* \brief The pressure of liquid water in \f$\mathrm{[Pa]}\f$ at a given density and
* temperature.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density density of component in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar liquidPressure(Scalar temperature, Scalar density, Scalar salinity)
{
// We use the newton method for this. For the initial value we
// assume the pressure to be 10% higher than the vapor
// pressure
Scalar pressure = 1.1*vaporPressure(temperature);
Scalar eps = pressure*1e-7;
Scalar deltaP = pressure*2;
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
Scalar f = liquidDensity(temperature, pressure) - density;
Scalar df_dp;
df_dp = liquidDensity(temperature, pressure + eps);
df_dp -= liquidDensity(temperature, pressure - eps);
df_dp /= 2*eps;
deltaP = - f/df_dp;
pressure += deltaP;
}
assert(pressure > 0.0);
return pressure;
}
/*!
* \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of steam.
*
* \param temperature temperature of component
* \param pressure pressure of component
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
{ return H2O::gasViscosity(temperature, pressure); };
/*!
* \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure brine.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
* Equation given in: - Batzle & Wang (1992)
* - cited by: Bachu & Adams (2002)
* "Equations of State for basin geofluids"
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure, Scalar salinity)
{
if(temperature <= 275.) // regularisation
{ temperature = 275; }
salinity = std::abs(salinity);
Scalar T_C = temperature - 273.15;
if(salinity < 0.0)
{salinity = 0.0;}
Scalar A = (0.42*pow((pow(salinity, 0.8)-0.17), 2) + 0.045)*pow(T_C, 0.8);
Scalar mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*exp(-A);
assert(mu_brine > 0.0);
return mu_brine/1000.0; /* unit: Pa s */
}
};
} // end namespace
#endif
#ifndef NACL_HH
#define NACL_HH
/*!
* \file
*
* \brief A class for the NaCl properties
*/
#ifndef DUMUX_NACL_HH
#define DUMUX_NACL_HH
#include <dumux/common/exceptions.hh>
#include <dumux/material/components/component.hh>
#include <cmath>
#include <iostream>
namespace Dumux
{
/*!
* \brief A class for the NaCl properties
*/
template <class Scalar>
class NaCl : public Component<Scalar, NaCl<Scalar> >
{
public:
/*!
* \brief A human readable name for the NaCl.
*/
static const char *name()
{
return "NaCl";
}
/*!
* \brief The mass in [kg] for one mole of NaCl.
*/
static Scalar molarMass()
{
return 58.4428e-3 ;
} // kg/mol
/*!
* \brief The diffusion Coefficient of NaCl in water.
*/
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
{
return 2e-9;
}
/*!
* \brief The mass density of NaCl.
*/
static Scalar Density()
{
return (2165); /* 2165 kg/m³*/
}
};
} // end namespace
#endif
#endif // NACL_HH
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* 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. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Determines the fluid composition given the component
* fugacities and an arbitary equation of state.
*/
#ifndef DUMUX_COMPOSITION_FROM_FUGACITIES_2PNCMIN_HH
#define DUMUX_COMPOSITION_FROM_FUGACITIES_2PNCMIN_HH
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/common/valgrind.hh>
namespace Dumux {
/*!
* \brief Calculates the chemical equilibrium from the component
* fugacities in a phase.
*/
template <class Scalar, class FluidSystem>
class compositionFromFugacities2pncmin
{
enum {
numComponents = FluidSystem::numComponents,
numMajorComponents = FluidSystem::numMajorComponents
};
typedef typename FluidSystem::ParameterCache ParameterCache;
public:
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
/*!
* \brief Guess an initial value for the composition of the phase.
* \param fluidState Thermodynamic state of the fluids
* \param paramCache Container for cache parameters
* \param phaseIdx The phase index
* \param phasePresence The presence index of the reference phase
* \param fugVec fugacity vector of the component
*/
template <class FluidState>
static void guessInitial(FluidState &fluidState,
ParameterCache &paramCache,
int phaseIdx,
int phasePresence,
const ComponentVector &fugVec)
{
if (FluidSystem::isIdealMixture(phaseIdx))
return;
// Pure component fugacities
for (int i = 0; i < numComponents; ++ i) {
//std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
fluidState.setMoleFraction(phaseIdx,
i,
1.0/numComponents);
}
}
/*!
* \brief Calculates the chemical equilibrium from the component
* fugacities in a phase. This constraint solver is developed for drying scenarios where salt component is restricted to liquid phase and still for the sake for equilibrium calculation some residual salt must be considered in the gas phase. In such cases for existence of gas phase only, in the theoretical liquid phase, we set the mole fraction of salt to 1e-10.
* \param fluidState Thermodynamic state of the fluids
* \param paramCache Container for cache parameters
* \param phaseIdx The phase index
* \param targetFug target fugacity
*
* The phase's fugacities must already be set.
*/
template <class FluidState>
static void solve(FluidState &fluidState,
ParameterCache &paramCache,
int phaseIdx,
int phasePresence,
const ComponentVector &targetFug)
{
// use a much more efficient method in case the phase is an
// ideal mixture
if (FluidSystem::isIdealMixture(phaseIdx))
{
solveIdealMix_(fluidState, paramCache, phaseIdx, phasePresence, targetFug);
return;
}
else
DUNE_THROW(NumericalProblem, "This constraint solver is not tested for non-ideal mixtures: Please refer computefromfugacities.hh for details" );
}
protected:
// update the phase composition in case the phase is an ideal
// mixture, i.e. the component's fugacity coefficients are
// independent of the phase's composition.
template <class FluidState>
static void solveIdealMix_(FluidState &fluidState,
ParameterCache &paramCache,
int phaseIdx,
int phasePresence,
const ComponentVector &fugacities)
{
for (int i = 0; i < numComponents; ++ i) {
Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
paramCache,
phaseIdx,
i);
Scalar gamma = phi * fluidState.pressure(phaseIdx);
fluidState.setFugacityCoefficient(phaseIdx, i, phi);
fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
// Special situation for drying PM and n-phase only situation.set the mole fraction of salt to 1e-10.
if (phaseIdx == 0 && i >= numMajorComponents && phasePresence == 2 /*nPhaseOnly*/)
fluidState.setMoleFraction(phaseIdx, i, 1.0e-10);
}
paramCache.updatePhase(fluidState, phaseIdx);
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, rho);
return;
}
};
} // end namespace Dumux
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* 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. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Computes all quantities of a generic fluid state if a
* reference phase has been specified.
*
* This makes it is possible to specify just one phase and let the
* remaining ones be calculated by the constraint solver. This
* constraint solver assumes thermodynamic equilibrium
*/
#ifndef DUMUX_COMPUTE_FROM_REFERENCE_PHASE_2PNC_HH
#define DUMUX_COMPUTE_FROM_REFERENCE_PHASE_2PNC_HH
#include <dumux/material/constraintsolvers/compositionfromfugacities.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/common/valgrind.hh>
namespace Dumux {
/*!
* \brief Computes all quantities of a generic fluid state if a
* reference phase has been specified.
*
* This makes it is possible to specify just one phase and let the
* remaining ones be calculated by the constraint solver. This
* constraint solver assumes thermodynamic equilibrium. It assumes the
* following quantities to be set:
*
* - composition (mole+mass fractions) of the *reference* phase
* - temperature of the *reference* phase
* - saturations of *all* phases
* - pressures of *all* phases
*