OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

- added formulation for p,T dependence of thermal conductivity according to IAPWS - can be used tabulated or directly - nitrogen is still a constant value (with proper reference though) - thermal conductivity of the mixture is obtained by mass fraction weighting in case of the gas, liquid uses thermal conductivity of water

Browse Source
This commit is contained in:
Philipp Nuske
2012-01-23 16:51:25 +00:00
committed by Andreas Lauser
parent 4ebe684372
commit ca400de662
7 changed files with 244 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

52
dumux/material/components/h2o.hh
View File

@@ -773,6 +773,58 @@ public:
return Common::viscosity(temperature, rho);
};
/*!
* \brief Thermal conductivity \f$\mathrm{[[W/(m K)]}\f$ of water (IAPWS) .
*
* Implementation taken from:
* freesteam - IAPWS-IF97 steam tables library
* Copyright (C) 2004-2009 John Pye
*
* Appendix B: Recommended Interpolating equation for Industrial Use
* see http://www.iapws.org/relguide/thcond.pdf
*
* \param temperature absolute temperature in K
* \param pressure of the phase in Pa
*/
static Scalar liquidThermalConductivity( Scalar temperature, Scalar pressure)
{
// Thermal conductivity of water is empirically fit.
// Evaluating that fitting-function outside the area of validity does not make sense.
assert( (pressure <= 400e6 and ((273.15<=temperature) and (temperature<=398.15)) )
or (pressure <= 200e6 and ((398.15<temperature) and (temperature<=523.15)) )
or (pressure <= 150e6 and ((523.15<temperature) and (temperature<=673.15)) )
or (pressure <= 100e6 and ((673.15<temperature) and (temperature<=1073.15)) ) );
Scalar rho = liquidDensity(temperature, pressure);
return Common::thermalConductivityIAPWS(temperature, rho);
}
/*!
* \brief Thermal conductivity \f$\mathrm{[[W/(m K)]}\f$ of water (IAPWS) .
*
* Implementation taken from:
* freesteam - IAPWS-IF97 steam tables library
* Copyright (C) 2004-2009 John Pye
*
* Appendix B: Recommended Interpolating equation for Industrial Use
* see http://www.iapws.org/relguide/thcond.pdf
*
* \param temperature absolute temperature in K
* \param pressure of the phase in Pa
*/
static Scalar gasThermalConductivity(const Scalar temperature, const Scalar pressure)
{
// Thermal conductivity of water is empirically fit.
// Evaluating that fitting-function outside the area of validity does not make sense.
assert( (pressure <= 400e6 and ((273.15<=temperature) and (temperature<=398.15)) )
or (pressure <= 200e6 and ((398.15<temperature) and (temperature<=523.15)) )
or (pressure <= 150e6 and ((523.15<temperature) and (temperature<=673.15)) )
or (pressure <= 100e6 and ((673.15<temperature) and (temperature<=1073.15)) ) );
Scalar rho = gasDensity(temperature, pressure);
return Common::thermalConductivityIAPWS(temperature, rho);
}
private:
// the unregularized specific enthalpy for liquid water
static Scalar enthalpyRegion1_(Scalar temperature, Scalar pressure)

85
dumux/material/components/iapws/common.hh
View File

@@ -151,6 +151,91 @@ public:
return 1e-6*muBar;
}
/*!
* \brief Thermal conductivity \f$\mathrm{[[W/(m^2 K/m)]}\f$ water (IAPWS) .
*
* Implementation taken from:
* freesteam - IAPWS-IF97 steam tables library
* Copyright (C) 2004-2009 John Pye
*
* Appendix B: Recommended Interpolating equation for Industrial Use
* see http://www.iapws.org/relguide/thcond.pdf
*
* \param T absolute temperature in K
* \param rho density of water in kg/m^3
*/
static Scalar thermalConductivityIAPWS(const Scalar T, const Scalar rho)
{
static constexpr Scalar thcond_tstar = 647.26 ;
static constexpr Scalar thcond_rhostar = 317.7 ;
/*static constexpr Scalar thcond_kstar = 1.0 ;*/
static constexpr Scalar thcond_b0 = -0.397070 ;
static constexpr Scalar thcond_b1 = 0.400302 ;
static constexpr Scalar thcond_b2 = 1.060000 ;
static constexpr Scalar thcond_B1 = -0.171587 ;
static constexpr Scalar thcond_B2 = 2.392190 ;
static constexpr Scalar thcond_c1 = 0.642857 ;
static constexpr Scalar thcond_c2 = -4.11717 ;
static constexpr Scalar thcond_c3 = -6.17937 ;
static constexpr Scalar thcond_c4 = 0.00308976 ;
static constexpr Scalar thcond_c5 = 0.0822994 ;
static constexpr Scalar thcond_c6 = 10.0932 ;
static constexpr Scalar thcond_d1 = 0.0701309 ;
static constexpr Scalar thcond_d2 = 0.0118520 ;
static constexpr Scalar thcond_d3 = 0.00169937 ;
static constexpr Scalar thcond_d4 = -1.0200 ;
static constexpr int thcond_a_count = 4;
static constexpr Scalar thcond_a[thcond_a_count] = {
0.0102811
,0.0299621
,0.0156146
,-0.00422464
};
Scalar Tbar = T / thcond_tstar;
Scalar rhobar = rho / thcond_rhostar;
/* fast implementation... minimised calls to 'pow' routine... */
Scalar Troot = sqrt(Tbar);
Scalar Tpow = Troot;
Scalar lam = 0;
for(int k = 0; k < thcond_a_count; ++k) {
lam += thcond_a[k] * Tpow;
Tpow *= Tbar;
}
lam += thcond_b0 + thcond_b1
* rhobar + thcond_b2
* exp(thcond_B1 * ((rhobar + thcond_B2)*(rhobar + thcond_B2)));
Scalar DTbar = abs(Tbar - 1) + thcond_c4;
Scalar DTbarpow = pow(DTbar, 3./5);
Scalar Q = 2. + thcond_c5 / DTbarpow;
Scalar S;
if(Tbar >= 1){
S = 1. / DTbar;
}else{
S = thcond_c6 / DTbarpow;
}
Scalar rhobar18 = pow(rhobar, 1.8);
Scalar rhobarQ = pow(rhobar, Q);
lam +=
(thcond_d1 / pow(Tbar,10) + thcond_d2) * rhobar18 *
exp(thcond_c1 * (1 - rhobar * rhobar18))
+ thcond_d3 * S * rhobarQ *
exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
+ thcond_d4 *
exp(thcond_c2 * pow(Troot,3) + thcond_c3 / pow(rhobar,5));
return /*thcond_kstar * */ lam;
}
};
template <class Scalar>

8
dumux/material/components/n2.hh
View File

@@ -181,6 +181,14 @@ public:
(cpVapB/2 + T*
(cpVapC/3 + T*
(cpVapD/4))));
//#warning NIST DATA STUPID INTERPOLATION
// Scalar T2 = 300.;
// Scalar T1 = 285.;
// Scalar h2 = 311200.;
// Scalar h1 = 295580.;
// Scalar h = h1+ (h2-h1) / (T2-T1) * (T-T1);
// return h ;
}
/*!

54
dumux/material/components/tabulatedcomponent.hh
View File

@@ -103,6 +103,8 @@ public:
liquidDensity_ = new Scalar[nTemp_*nPress_];
gasViscosity_ = new Scalar[nTemp_*nPress_];
liquidViscosity_ = new Scalar[nTemp_*nPress_];
gasThermalConductivity_ = new Scalar[nTemp_*nPress_];
liquidThermalConductivity_ = new Scalar[nTemp_*nPress_];
gasPressure_ = new Scalar[nTemp_*nDensity_];
liquidPressure_ = new Scalar[nTemp_*nDensity_];
@@ -141,6 +143,10 @@ public:
try { gasViscosity_[i] = RawComponent::gasViscosity(temperature, pressure); }
catch (Dune::NotImplemented) { gasViscosity_[i] = NaN; }
catch (NumericalProblem) { gasViscosity_[i] = NaN; };
try { gasThermalConductivity_[i] = RawComponent::gasThermalConductivity(temperature, pressure); }
catch (Dune::NotImplemented) { gasThermalConductivity_[i] = NaN; }
catch (NumericalProblem) { gasThermalConductivity_[i] = NaN; };
};
Scalar plMin = minLiquidPressure_(iT);
@@ -165,6 +171,10 @@ public:
try { liquidViscosity_[i] = RawComponent::liquidViscosity(temperature, pressure); }
catch (Dune::NotImplemented) { liquidViscosity_[i] = NaN; }
catch (NumericalProblem) { liquidViscosity_[i] = NaN; };
try { liquidThermalConductivity_[i] = RawComponent::liquidThermalConductivity(temperature, pressure); }
catch (Dune::NotImplemented) { liquidThermalConductivity_[i] = NaN; }
catch (NumericalProblem) { liquidThermalConductivity_[i] = NaN; };
}
}
@@ -491,6 +501,43 @@ public:
return result;
};
/*!
* \brief The thermal conductivity of gaseous water \f$\mathrm{[W / (m K)]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
{
Scalar result = interpolateGasTP_(gasThermalConductivity_,
temperature,
pressure);
if (std::isnan(result)) {
printWarning_("gasThermalConductivity", temperature, pressure);
return RawComponent::gasThermalConductivity(temperature, pressure);
}
return result;
};
/*!
* \brief The thermal conductivity of liquid water \f$\mathrm{[W / (m K)]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
{
Scalar result = interpolateLiquidTP_(liquidThermalConductivity_,
temperature,
pressure);
if (std::isnan(result)) {
printWarning_("liquidThermalConductivity", temperature, pressure);
return RawComponent::liquidThermalConductivity(temperature, pressure);
}
return result;
};
private:
// prints a warning if the result is not in range or the table has
// not been initialized
@@ -781,6 +828,9 @@ private:
static Scalar *gasViscosity_;
static Scalar *liquidViscosity_;
static Scalar *gasThermalConductivity_;
static Scalar *liquidThermalConductivity_;
// 2D fields with the temperature and density as degrees of
// freedom
static Scalar *gasPressure_;
@@ -835,6 +885,10 @@ Scalar* TabulatedComponent<Scalar, RawComponent, useVaporPressure>::gasViscosity
template <class Scalar, class RawComponent, bool useVaporPressure>
Scalar* TabulatedComponent<Scalar, RawComponent, useVaporPressure>::liquidViscosity_;
template <class Scalar, class RawComponent, bool useVaporPressure>
Scalar* TabulatedComponent<Scalar, RawComponent, useVaporPressure>::gasThermalConductivity_;
template <class Scalar, class RawComponent, bool useVaporPressure>
Scalar* TabulatedComponent<Scalar, RawComponent, useVaporPressure>::liquidThermalConductivity_;
template <class Scalar, class RawComponent, bool useVaporPressure>
Scalar* TabulatedComponent<Scalar, RawComponent, useVaporPressure>::gasPressure_;
template <class Scalar, class RawComponent, bool useVaporPressure>
Scalar* TabulatedComponent<Scalar, RawComponent, useVaporPressure>::liquidPressure_;

7
dumux/material/fluidstates/nonequilibriumfluidstate.hh
View File

@@ -167,7 +167,10 @@ public:
* \brief The specific internal energy of a fluid phase [J/kg]
*/
Scalar internalEnergy(int phaseIdx) const
{ return enthalpy_[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); }
{
return enthalpy_[phaseIdx]
- pressure(phaseIdx)/density(phaseIdx);
}
/*!
* \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -300,7 +303,7 @@ public:
Valgrind::CheckDefined(saturation_[i]);
Valgrind::CheckDefined(density_[i]);
Valgrind::CheckDefined(temperature_[i]);
//Valgrind::CheckDefined(internalEnergy_[i]);
Valgrind::CheckDefined(enthalpy_[i]);
Valgrind::CheckDefined(viscosity_[i]);
}
#endif // HAVE_VALGRIND

2
dumux/material/fluidsystems/basefluidsystem.hh
View File

@@ -144,7 +144,7 @@ public:
}
/*!
* \brief Thermal conductivity of a fluid phase [W/(m^2 K/m)].
* \brief Thermal conductivity of a fluid phase [W/(m K)].
*
* Use the conductivity of air and water as a first approximation.
* Source:

53
dumux/material/fluidsystems/h2on2fluidsystem.hh
View File

@@ -176,7 +176,7 @@ public:
//! The components for pure water
typedef TabulatedH2O H2O;
//typedef SimpleH2O H2O;
//typedef IapwsH2O H2O;
// typedef IapwsH2O H2O;
//! The components for pure nitrogen
typedef SimpleN2 N2;
@@ -629,11 +629,9 @@ public:
}
/*!
* \brief Thermal conductivity of a fluid phase [W/(m^2 K/m)].
* \brief Thermal conductivity of a fluid phase [W/(m K)].
*
* Use the conductivity of air and water as a first approximation.
* Source:
* http://en.wikipedia.org/wiki/List_of_thermal_conductivities
*
* \param fluidState An abitrary fluid state
* \param phaseIdx The index of the fluid phase to consider
@@ -641,18 +639,45 @@ public:
using Base::thermalConductivity;
template <class FluidState>
static Scalar thermalConductivity(const FluidState &fluidState,
int phaseIdx)
const int phaseIdx)
{
// TODO thermal conductivity is a function of:
// Scalar p = fluidState.pressure(phaseIdx);
// Scalar T = fluidState.temperature(phaseIdx);
// Scalar x = fluidState.moleFrac(phaseIdx,compIdx);
//#warning: so far rough estimates from wikipedia
if (phaseIdx == lPhaseIdx)
return 0.6; // conductivity of water[W / (m K ) ]
assert(0 <= phaseIdx && phaseIdx < numPhases);
// gas phase
return 0.025; // conductivity of air [W / (m K ) ]
if (phaseIdx == lPhaseIdx){// liquid phase
if(useComplexRelations){
const Scalar & temperature = fluidState.temperature(phaseIdx) ;
const Scalar & pressure = fluidState.pressure(phaseIdx);
return H2O::liquidThermalConductivity(temperature, pressure);
}
else
return 0.578078; // conductivity of water[W / (m K ) ] IAPWS evaluated at p=.1 MPa, T=8°C
}
else{// gas phase
// Isobaric Properties for Nitrogen in: NIST Standard Reference Database Number 69, Eds. P.J. Linstrom
// and W.G. Mallard
// evaluated at p=.1 MPa, T=8°C, does not change dramatically with p,T
const Scalar & lambdaPureNitrogen = 0.024572;
if (useComplexRelations){
const Scalar & xNitrogen = fluidState.moleFraction(phaseIdx, N2Idx);
const Scalar & xWater = fluidState.moleFraction(phaseIdx, H2OIdx);
const Scalar & lambdaNitrogen = xNitrogen * lambdaPureNitrogen;
// Assuming Raoult's, Daltons law and ideal gas
// in order to obtain the partial density of water in the air phase
const Scalar & temperature = fluidState.temperature(phaseIdx) ;
const Scalar & pressure = fluidState.pressure(phaseIdx);
const Scalar & averageMolarMass = fluidState.averageMolarMass(gPhaseIdx);
const Scalar & partialPressure = pressure * xWater;
const Scalar & lambdaWater = xWater * H2O::gasThermalConductivity(temperature,
partialPressure);
return lambdaNitrogen + lambdaWater;
}
else
return lambdaPureNitrogen; // conductivity of Nitrogen [W / (m K ) ]
}
}
/*!