diff --git a/dumux/boxmodels/2p2c/2p2cmodel.hh b/dumux/boxmodels/2p2c/2p2cmodel.hh
index d9bb29bd004d920f5a8f93a6db8fc9b9169f2862..e7bfdedfe12874301fdc5edf7c915b2d4082b56e 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -219,9 +219,6 @@ public:
setSwitched_(false);
resetPhasePresence_();
- /*this->localJacobian().updateStaticData(this->curSolFunction(),
- this->prevSolFunction());
- */
};
/*!
diff --git a/dumux/boxmodels/3p3c/3p3cindices.hh b/dumux/boxmodels/3p3c/3p3cindices.hh
index a4d2081ca0c3d78a89e8d2226de4d21103ba70c1..570dcb8bef04a67f7f0ea647f6c6bd6b75bd21c8 100644
--- a/dumux/boxmodels/3p3c/3p3cindices.hh
+++ b/dumux/boxmodels/3p3c/3p3cindices.hh
@@ -35,25 +35,13 @@
namespace Dumux
{
-/*!
- * \brief Enumerates the formulations which the 3p3c model accepts.
- */
-struct ThreePThreeCFormulation
-{
- enum {
- pgSwSn
- };
-};
-
/*!
* \brief The indices for the isothermal ThreePThreeC model.
*
* \tparam formulation The formulation, only pgSwSn
* \tparam PVOffset The first index in a primary variable vector.
*/
-template <class TypeTag,
- int formulation = ThreePThreeCFormulation::pgSwSn,
- int PVOffset = 0>
+template <class TypeTag, int PVOffset = 0>
class ThreePThreeCIndices
{
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
@@ -87,9 +75,10 @@ public:
static const int SOrX2Idx = switch2Idx; //!< Index of the either the saturation of the gas phase or the mass fraction secondary component if a phase is not present
// equation indices
- static const int contiWEqIdx = PVOffset + wCompIdx; //!< Index of the mass conservation equation for the water component
- static const int contiCEqIdx = PVOffset + cCompIdx; //!< Index of the mass conservation equation for the contaminant component
- static const int contiAEqIdx = PVOffset + aCompIdx; //!< Index of the mass conservation equation for the air component
+ static const int conti0EqIdx = PVOffset; //!< Index of the mass conservation equation for the first component
+ static const int contiWEqIdx = conti0EqIdx + wCompIdx; //!< Index of the mass conservation equation for the water component
+ static const int contiCEqIdx = conti0EqIdx + cCompIdx; //!< Index of the mass conservation equation for the contaminant component
+ static const int contiAEqIdx = conti0EqIdx + aCompIdx; //!< Index of the mass conservation equation for the air component
};
}
diff --git a/dumux/boxmodels/3p3c/3p3clocalresidual.hh b/dumux/boxmodels/3p3c/3p3clocalresidual.hh
index c9be2cba68ffa40a6b0b4d73dd2b05b5fd679868..a8d334f6204a6cf366e8a3a66555b8219e495da2 100644
--- a/dumux/boxmodels/3p3c/3p3clocalresidual.hh
+++ b/dumux/boxmodels/3p3c/3p3clocalresidual.hh
@@ -28,8 +28,8 @@
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the two-phase two-component box model.
*/
-#ifndef DUMUX_NEW_3P3C_LOCAL_RESIDUAL_BASE_HH
-#define DUMUX_NEW_3P3C_LOCAL_RESIDUAL_BASE_HH
+#ifndef DUMUX_3P3C_LOCAL_RESIDUAL_HH
+#define DUMUX_3P3C_LOCAL_RESIDUAL_HH
#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/common/math.hh>
@@ -59,41 +59,33 @@ class ThreePThreeCLocalResidual: public BoxLocalResidual<TypeTag>
protected:
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
-
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-
typedef typename GET_PROP_TYPE(TypeTag, ThreePThreeCIndices) Indices;
- enum
- {
-
- numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
- numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
-
-
- contiWEqIdx = Indices::contiWEqIdx,
- contiCEqIdx = Indices::contiCEqIdx,
- contiAEqIdx = Indices::contiAEqIdx,
+ enum {
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
+ conti0EqIdx = Indices::conti0EqIdx,
+ contiWEqIdx = Indices::contiWEqIdx,
+ contiCEqIdx = Indices::contiCEqIdx,
+ contiAEqIdx = Indices::contiAEqIdx,
- wCompIdx = Indices::wCompIdx,
- cCompIdx = Indices::cCompIdx,
- aCompIdx = Indices::aCompIdx
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
- };
+ wCompIdx = Indices::wCompIdx,
+ cCompIdx = Indices::cCompIdx,
+ aCompIdx = Indices::aCompIdx
+ };
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
-
public:
-
/*!
* \brief Evaluate the amount all conservation quantities
* (e.g. phase mass) within a sub-control volume.
@@ -112,7 +104,9 @@ public:
// used. The secondary variables are used accordingly. This
// is required to compute the derivative of the storage term
// using the implicit euler method.
- const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
+ const ElementVolumeVariables &elemVolVars =
+ usePrevSol
+ ? this->prevVolVars_()
: this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
@@ -120,48 +114,13 @@ public:
result = 0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
- switch (compIdx) {
- case wCompIdx:
- {
- Scalar waterContrib = 0.0;
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- waterContrib += volVars.porosity() * volVars.saturation(phaseIdx)
- * volVars.molarDensity(phaseIdx)
- * volVars.fluidState().moleFraction(phaseIdx, compIdx);
- }
- result[contiWEqIdx] += waterContrib;
- break;
- }
-
- case cCompIdx:
- {
- Scalar contContrib = 0.0;
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- contContrib += volVars.porosity() * volVars.saturation(phaseIdx)
- * volVars.molarDensity(phaseIdx)
- * volVars.fluidState().moleFraction(phaseIdx, compIdx);
- }
- contContrib += volVars.bulkDensTimesAdsorpCoeff()
- * volVars.fluidState().moleFraction(wPhaseIdx, compIdx)
- * volVars.molarDensity(wPhaseIdx);
- result[contiCEqIdx] += contContrib;
- break;
- }
-
- case aCompIdx:
- {
- Scalar airContrib = 0.0;
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- airContrib += volVars.porosity() * volVars.saturation(phaseIdx)
- * volVars.molarDensity(phaseIdx)
- * volVars.fluidState().moleFraction(phaseIdx, compIdx);
- }
- result[contiAEqIdx] += airContrib;
- break;
- }
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ result[conti0EqIdx + compIdx] +=
+ volVars.porosity()
+ * volVars.saturation(phaseIdx)
+ * volVars.molarDensity(phaseIdx)
+ * volVars.fluidState().moleFraction(phaseIdx, compIdx);
}
}
}
@@ -205,33 +164,27 @@ public:
{
// data attached to upstream and the downstream vertices
// of the current phase
- const VolumeVariables &up =
- this->curVolVars_(vars.upstreamIdx(phaseIdx));
- const VolumeVariables &dn =
- this->curVolVars_(vars.downstreamIdx(phaseIdx));
+ const VolumeVariables &up = this->curVolVars_(vars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn = this->curVolVars_(vars.downstreamIdx(phaseIdx));
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
- int eqIdx;
- if (compIdx == wCompIdx) eqIdx = contiWEqIdx;
- else if (compIdx == cCompIdx) eqIdx = contiCEqIdx;
- else eqIdx = contiAEqIdx;
-
// add advective flux of current component in current
// phase
// if alpha > 0 und alpha < 1 then both upstream and downstream
// nodes need their contribution
// if alpha == 1 (which is mostly the case) then, the downstream
// node is not evaluated
+ int eqIdx = conti0EqIdx + compIdx;
flux[eqIdx] += vars.KmvpNormal(phaseIdx)
* (massUpwindWeight
- * up.molarDensity(phaseIdx)
* up.mobility(phaseIdx)
+ * up.fluidState().molarDensity(phaseIdx)
* up.fluidState().moleFraction(phaseIdx, compIdx)
+
(1.0 - massUpwindWeight)
- * dn.molarDensity(phaseIdx)
* dn.mobility(phaseIdx)
+ * dn.fluidState().molarDensity(phaseIdx)
* dn.fluidState().moleFraction(phaseIdx, compIdx));
}
}
@@ -247,7 +200,9 @@ public:
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
{
- // TODO wo kommt das denn her? Dune::FieldMatrix<Scalar, numPhases, numComponents> averagedPorousDiffCoeffMatrix = vars.porousDiffCoeff();
+#warning "Disabled for now"
+ return;
+ // TODO: reference!? Dune::FieldMatrix<Scalar, numPhases, numComponents> averagedPorousDiffCoeffMatrix = vars.porousDiffCoeff();
// add diffusive flux of gas component in liquid phase
Scalar tmp = - vars.porousDiffCoeff()[wPhaseIdx][aCompIdx] * vars.molarDensityAtIP(wPhaseIdx);
tmp *= (vars.molarAConcGrad(wPhaseIdx) * vars.face().normal);
@@ -312,16 +267,15 @@ public:
}
protected:
-
Implementation *asImp_()
{
return static_cast<Implementation *> (this);
}
+
const Implementation *asImp_() const
{
return static_cast<const Implementation *> (this);
}
-
};
} // end namepace
diff --git a/dumux/boxmodels/3p3c/3p3cmodel.hh b/dumux/boxmodels/3p3c/3p3cmodel.hh
index 089c8450625bc2fd9b193193dfa92030924c730e..55a374e7569f911d6c8445fb47f5e6c040eeecb1 100644
--- a/dumux/boxmodels/3p3c/3p3cmodel.hh
+++ b/dumux/boxmodels/3p3c/3p3cmodel.hh
@@ -212,9 +212,6 @@ public:
setSwitched_(false);
resetPhasePresence_();
- /*this->localJacobian().updateStaticData(this->curSolFunction(),
- this->prevSolFunction());
- */
};
/*!
@@ -264,7 +261,9 @@ public:
*/
int phasePresence(int globalVertexIdx, bool oldSol) const
{
- return oldSol ? staticVertexDat_[globalVertexIdx].oldPhasePresence
+ return
+ oldSol
+ ? staticVertexDat_[globalVertexIdx].oldPhasePresence
: staticVertexDat_[globalVertexIdx].phasePresence;
}
@@ -276,7 +275,6 @@ public:
* \param sol The solution vector
* \param writer The writer for multi-file VTK datasets
*/
-
template<class MultiWriter>
void addOutputVtkFields(const SolutionVector &sol,
MultiWriter &writer)
@@ -286,12 +284,16 @@ public:
// create the required scalar fields
unsigned numVertices = this->problem_().gridView().size(dim);
- ScalarField *Sw = writer.allocateManagedBuffer (numVertices);
- ScalarField *Sn = writer.allocateManagedBuffer (numVertices);
- ScalarField *Sg = writer.allocateManagedBuffer (numVertices);
- ScalarField *pg = writer.allocateManagedBuffer (numVertices);
- ScalarField *pn = writer.allocateManagedBuffer (numVertices);
- ScalarField *pw = writer.allocateManagedBuffer (numVertices);
+ ScalarField *saturation[numPhases];
+ ScalarField *pressure[numPhases];
+ ScalarField *density[numPhases];
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ saturation[phaseIdx] = writer.allocateManagedBuffer(numVertices);
+ pressure[phaseIdx] = writer.allocateManagedBuffer(numVertices);
+ density[phaseIdx] = writer.allocateManagedBuffer(numVertices);
+ }
+
ScalarField *phasePresence = writer.allocateManagedBuffer (numVertices);
ScalarField *moleFraction[numPhases][numComponents];
for (int i = 0; i < numPhases; ++i)
@@ -326,44 +328,49 @@ public:
fvElemGeom,
i,
false);
- (*Sw)[globalIdx] = volVars.saturation(wPhaseIdx);
- (*Sn)[globalIdx] = volVars.saturation(nPhaseIdx);
- (*Sg)[globalIdx] = volVars.saturation(gPhaseIdx);
- (*pg)[globalIdx] = volVars.pressure(gPhaseIdx);
- (*pn)[globalIdx] = volVars.pressure(nPhaseIdx);
- (*pw)[globalIdx] = volVars.pressure(wPhaseIdx);
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- (*moleFraction[phaseIdx][compIdx])[globalIdx]
- = volVars.fluidState().moleFraction(phaseIdx,
- compIdx);
-
- Valgrind::CheckDefined(
- (*moleFraction[phaseIdx][compIdx])[globalIdx][0]);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ (*saturation[phaseIdx])[globalIdx] = volVars.fluidState().saturation(phaseIdx);
+ (*pressure[phaseIdx])[globalIdx] = volVars.fluidState().pressure(phaseIdx);
+ (*density[phaseIdx])[globalIdx] = volVars.fluidState().density(phaseIdx);
+ }
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ (*moleFraction[phaseIdx][compIdx])[globalIdx] =
+ volVars.fluidState().moleFraction(phaseIdx,
+ compIdx);
+
+ Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[globalIdx]);
}
+ }
+
(*poro)[globalIdx] = volVars.porosity();
(*perm)[globalIdx] = volVars.permeability();
(*temperature)[globalIdx] = volVars.temperature();
- (*phasePresence)[globalIdx]
- = staticVertexDat_[globalIdx].phasePresence;
+ (*phasePresence)[globalIdx] = staticVertexDat_[globalIdx].phasePresence;
};
}
- writer.attachVertexData(*Sw, "Sw");
- writer.attachVertexData(*Sn, "Sn");
- writer.attachVertexData(*Sg, "Sg");
- writer.attachVertexData(*pg, "pg");
- writer.attachVertexData(*pn, "pn");
- writer.attachVertexData(*pw, "pw");
+ writer.attachVertexData(*saturation[wPhaseIdx], "Sw");
+ writer.attachVertexData(*saturation[nPhaseIdx], "Sn");
+ writer.attachVertexData(*saturation[gPhaseIdx], "Sg");
+ writer.attachVertexData(*pressure[wPhaseIdx], "pw");
+ writer.attachVertexData(*pressure[nPhaseIdx], "pn");
+ writer.attachVertexData(*pressure[gPhaseIdx], "pg");
+ writer.attachVertexData(*density[wPhaseIdx], "rhow");
+ writer.attachVertexData(*density[nPhaseIdx], "rhon");
+ writer.attachVertexData(*density[gPhaseIdx], "rhog");
+
for (int i = 0; i < numPhases; ++i)
{
for (int j = 0; j < numComponents; ++j)
{
std::ostringstream oss;
- oss << "X_"
+ oss << "x^"
<< FluidSystem::phaseName(i)
+ << "_"
<< FluidSystem::componentName(j);
writer.attachVertexData(*moleFraction[i][j], oss.str().c_str());
}
diff --git a/dumux/boxmodels/3p3c/3p3cproperties.hh b/dumux/boxmodels/3p3c/3p3cproperties.hh
index 97b35e224cc393829e83b9ac3e4f4f071d8a9e2e..a330cd67754140a4f8190d1cb3c2d7d2416de41f 100644
--- a/dumux/boxmodels/3p3c/3p3cproperties.hh
+++ b/dumux/boxmodels/3p3c/3p3cproperties.hh
@@ -55,7 +55,6 @@ NEW_TYPE_TAG(BoxThreePThreeC, INHERITS_FROM(BoxModel));
NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
NEW_PROP_TAG(ThreePThreeCIndices); //!< Enumerations for the 3p3c models
-NEW_PROP_TAG(Formulation); //!< The formulation of the model
NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters
NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
diff --git a/dumux/boxmodels/3p3c/3p3cpropertydefaults.hh b/dumux/boxmodels/3p3c/3p3cpropertydefaults.hh
index 5ef51f5f383fdd8cbe5d9f227e36f60f142785f6..3d3a4a0da6d1ca6ad039ebab78f59bff31ea0a38 100644
--- a/dumux/boxmodels/3p3c/3p3cpropertydefaults.hh
+++ b/dumux/boxmodels/3p3c/3p3cpropertydefaults.hh
@@ -90,28 +90,14 @@ SET_PROP(BoxThreePThreeC, NumPhases)
SET_INT_PROP(BoxThreePThreeC, NumEq, 3); //!< set the number of equations to 2
-//! Set the default formulation to pg-Sw-Sn
-SET_INT_PROP(BoxThreePThreeC,
- Formulation,
- ThreePThreeCFormulation::pgSwSn);
-
/*!
* \brief Set the property for the material parameters by extracting
* it from the material law.
*/
-SET_PROP(BoxThreePThreeC, MaterialLawParams)
-{
- private:
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
-
- public:
- typedef typename MaterialLaw::Params type;
-};
+SET_TYPE_PROP(BoxThreePThreeC, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
-//! Use the 3p3c local jacobian operator for the 3p3c model
-SET_TYPE_PROP(BoxThreePThreeC,
- LocalResidual,
- ThreePThreeCLocalResidual<TypeTag>);
+//! The local residual function of the conservation equations
+SET_TYPE_PROP(BoxThreePThreeC, LocalResidual, ThreePThreeCLocalResidual<TypeTag>);
//! Use the 3p3c specific newton controller for the 3p3c model
SET_TYPE_PROP(BoxThreePThreeC, NewtonController, ThreePThreeCNewtonController<TypeTag>);
@@ -125,21 +111,11 @@ SET_TYPE_PROP(BoxThreePThreeC, VolumeVariables, ThreePThreeCVolumeVariables<Type
//! the FluxVariables property
SET_TYPE_PROP(BoxThreePThreeC, FluxVariables, ThreePThreeCFluxVariables<TypeTag>);
-//! the BoundaryVariables property
-// SET_TYPE_PROP(BoxThreePThreeC, BoundaryVariables, ThreePThreeCBoundaryVariables<TypeTag>);
-
//! the upwind factor for the mobility.
SET_SCALAR_PROP(BoxThreePThreeC, MassUpwindWeight, 1.0);
//! The indices required by the isothermal 3p3c model
-SET_PROP(BoxThreePThreeC,
- ThreePThreeCIndices)
-{ private:
- enum { Formulation = GET_PROP_VALUE(TypeTag, Formulation) };
- public:
- typedef ThreePThreeCIndices<TypeTag, Formulation, 0> type;
-};
-
+SET_TYPE_PROP(BoxThreePThreeC, ThreePThreeCIndices, ThreePThreeCIndices<TypeTag, /*PVOffset=*/0>);
}
}
diff --git a/dumux/material/components/mesitylene.hh b/dumux/material/components/mesitylene.hh
index 3c56b715dfa60635d98ba6edeffc5aa3e306b0f8..6201340d701e45d670d9a0b0c9e936eb0d94a6a4 100644
--- a/dumux/material/components/mesitylene.hh
+++ b/dumux/material/components/mesitylene.hh
@@ -107,8 +107,8 @@ public:
*/
static Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
{
- const Scalar DTemp = temperature - 273.0; // K -> degC
- return 0.5*DTemp*(spHeatCapLiquidPhase_(temperature) + spHeatCapLiquidPhase_(temperature));
+ Scalar DTemp = temperature-273.0; // K -> degC
+ return 0.5*DTemp*(spHeatCapLiquidPhase_(0.2113*DTemp) + spHeatCapLiquidPhase_(0.7887*DTemp));
}
/*!
@@ -144,8 +144,21 @@ public:
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{
return molarLiquidDensity_(temperature)*molarMass(); // [kg/m^3]
+
}
+ /*!
+ * \brief Returns true iff the gas phase is assumed to be compressible
+ */
+ static bool gasIsCompressible()
+ { return true; }
+
+ /*!
+ * \brief Returns true iff the gas phase is assumed to be ideal
+ */
+ static bool gasIsIdeal()
+ { return true; }
+
/*!
* \brief Returns true iff the liquid phase is assumed to be compressible
*/
@@ -161,8 +174,8 @@ public:
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure, bool regularize=true)
{
- temp = std::min(temperature, 500.0); // regularization
- temp = std::max(temperature, 250.0);
+ temperature = std::min(temperature, 500.0); // regularization
+ temperature = std::max(temperature, 250.0);
// reduced temperature
Scalar Tr = temperature/criticalTemperature();
@@ -186,13 +199,13 @@ public:
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
{
- temp = std::min(temperature, 500.0); // regularization
- temp = std::max(temperature, 250.0);
+ temperature = std::min(temperature, 500.0); // regularization
+ temperature = std::max(temperature, 250.0);
const Scalar A = -6.749;
const Scalar B = 2010.0;
- return std::exp(A + B/temp)*1e-3; // [Pa s]
+ return std::exp(A + B/temperature)*1e-3; // [Pa s]
}
protected:
@@ -242,31 +255,30 @@ protected:
*/
static Scalar spHeatCapLiquidPhase_(Scalar temperature)
{
- // according Reid et al. : Missenard group contrib. method (s. example 5-8)
- // Mesitylen: C9H12 : 3* CH3 ; 1* C6H5 (phenyl-ring) ; -2* H (this was to much!)
- // linear interpolation between table values [J/(mol K)]
-
+ /* according Reid et al. : Missenard group contrib. method (s. example 5-8) */
+ /* Mesitylen: C9H12 : 3* CH3 ; 1* C6H5 (phenyl-ring) ; -2* H (this was to much!) */
+ /* linear interpolation between table values [J/(mol K)]*/
Scalar H, CH3, C6H5;
- if(temperature < 298.0) {
- // extrapolation for Temperature<273
- H = 13.4 + 1.2*(temperature - 273.0)/25.0;
- CH3 = 40.0 + 1.6*(temperature - 273.0)/25.0;
- C6H5 = 113.0 + 4.2*(temperature - 273.0)/25.0;
+ if(temperature<298.) {
+ // extrapolation for Temperature<273 */
+ H = 13.4+1.2*(temperature-273.0)/25.;
+ CH3 = 40.0+1.6*(temperature-273.0)/25.;
+ C6H5 = 113.0+4.2*(temperature-273.0)/25.;
}
- else if(temperature < 323.0){
- H = 14.6 + 0.9*(temperature - 298.0)/25.0;
- CH3 = 41.6 + 1.9*(temperature - 298.0)/25.0;
- C6H5 = 117.2 + 6.2*(temperature - 298.0)/25.0;
+ else if((temperature>=298.0)&&(temperature<323.)){
+ H = 14.6+0.9*(temperature-298.0)/25.;
+ CH3 = 41.6+1.9*(temperature-298.0)/25.;
+ C6H5 = 117.2+6.2*(temperature-298.0)/25.;
}
- else if(temperature < 348.0){
- H = 15.5 + 1.2*(temperature - 323.0)/25.0;
- CH3 = 43.5 + 2.3*(temperature - 323.0)/25.0;
- C6H5 = 123.4 + 6.3*(temperature - 323.0)/25.0;
+ else if((temperature>=323.0)&&(temperature<348.)){
+ H = 15.5+1.2*(temperature-323.0)/25.;
+ CH3 = 43.5+2.3*(temperature-323.0)/25.;
+ C6H5 = 123.4+6.3*(temperature-323.0)/25.;
}
- else { // extrapolation for Temperature>373
- H = 16.7 + 2.1*(temperature - 348.0)/25.0; // leads probably to underestimates
- CH3 = 45.8 + 2.5*(temperature - 348.0)/25.0;
- C6H5 = 129.7 + 6.3*(temperature - 348.0)/25.0;
+ else if(temperature>=348.0){ /* take care: extrapolation for Temperature>373 */
+ H = 16.7+2.1*(temperature-348.0)/25.; /* leads probably to underestimates */
+ CH3 = 45.8+2.5*(temperature-348.0)/25.;
+ C6H5 = 129.7+6.3*(temperature-348.0)/25.;
}
return (C6H5 + 3*CH3 - 2*H)/molarMass();
diff --git a/dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh b/dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh
index e6a7e872b14c775c05c89238905c9c8a510339e2..5a55719211cd06f321efa9bb0fb71f60c50375ce 100644
--- a/dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh
+++ b/dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh
@@ -31,7 +31,7 @@
#include <dumux/material/idealgas.hh>
#include <dumux/material/components/air.hh>
-#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/mesitylene.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
@@ -58,7 +58,7 @@ class H2OAirMesitylene
typedef BaseFluidSystem<Scalar, ThisType> Base;
public:
- typedef Dumux::H2O<Scalar> H2O;
+ typedef Dumux::SimpleH2O<Scalar> H2O;
typedef Dumux::Mesitylene<Scalar> NAPL;
typedef Dumux::Air<Scalar> Air;
@@ -88,7 +88,7 @@ public:
}
static bool isIdealGas(int phaseIdx)
- { return !isLiquid(phaseIdx); }
+ { return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
/*!
* \brief Returns true if and only if a fluid phase is assumed to
@@ -125,7 +125,7 @@ public:
static bool isCompressible(int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
- // ideal gases are always compressible
+ // gases are always compressible
if (phaseIdx == gPhaseIdx)
return true;
else if (phaseIdx == wPhaseIdx)
@@ -203,23 +203,21 @@ public:
Scalar pressure = NAPL::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
return NAPL::liquidDensity(fluidState.temperature(phaseIdx), pressure);
}
- else if (phaseIdx == gPhaseIdx) {
- Scalar fugH2O =
- fluidState.moleFraction(gPhaseIdx, H2OIdx) *
- fluidState.pressure(gPhaseIdx);
- Scalar fugAir =
- fluidState.moleFraction(gPhaseIdx, airIdx) *
- fluidState.pressure(gPhaseIdx);
- Scalar fugNAPL =
- fluidState.moleFraction(gPhaseIdx, NAPLIdx) *
- fluidState.pressure(gPhaseIdx);
- return
- H2O::gasDensity(fluidState.temperature(phaseIdx), fugH2O) +
- Air::gasDensity(fluidState.temperature(phaseIdx), fugAir) +
- NAPL::gasDensity(fluidState.temperature(phaseIdx), fugNAPL);
- }
- DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ assert (phaseIdx == gPhaseIdx);
+ Scalar pH2O =
+ fluidState.moleFraction(gPhaseIdx, H2OIdx) *
+ fluidState.pressure(gPhaseIdx);
+ Scalar pAir =
+ fluidState.moleFraction(gPhaseIdx, airIdx) *
+ fluidState.pressure(gPhaseIdx);
+ Scalar pNAPL =
+ fluidState.moleFraction(gPhaseIdx, NAPLIdx) *
+ fluidState.pressure(gPhaseIdx);
+ return
+ H2O::gasDensity(fluidState.temperature(phaseIdx), pH2O) +
+ Air::gasDensity(fluidState.temperature(phaseIdx), pAir) +
+ NAPL::gasDensity(fluidState.temperature(phaseIdx), pNAPL);
}
/*!
@@ -239,55 +237,55 @@ public:
// assume pure NAPL viscosity
return NAPL::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
}
- else if (phaseIdx == gPhaseIdx) {
- /* Wilke method. See:
- *
- * See: R. Reid, et al.: The Properties of Gases and Liquids,
- * 4th edition, McGraw-Hill, 1987, 407-410
- * 5th edition, McGraw-Hill, 20001, p. 9.21/22
- *
- * in this case, we use a simplified version in order to avoid
- * computationally costly evaluation of sqrt and pow functions and
- * divisions
- * -- compare e.g. with Promo Class p. 32/33
- */
- Scalar muResult;
- const Scalar mu[numComponents] = {
- H2O::gasViscosity(fluidState.temperature(phaseIdx), H2O::vaporPressure(fluidState.temperature(phaseIdx))),
- Air::simpleGasViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)),
- NAPL::gasViscosity(fluidState.temperature(phaseIdx), NAPL::vaporPressure(fluidState.temperature(phaseIdx)))
- };
- // molar masses
- const Scalar M[numComponents] = {
- H2O::molarMass(),
- Air::molarMass(),
- NAPL::molarMass()
- };
- Scalar muAW = mu[airIdx]*fluidState.moleFraction(gPhaseIdx, airIdx)
- + mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
- / (fluidState.moleFraction(gPhaseIdx, airIdx)
- + fluidState.moleFraction(gPhaseIdx, H2OIdx));
- Scalar xAW = fluidState.moleFraction(gPhaseIdx, airIdx)
- + fluidState.moleFraction(gPhaseIdx, H2OIdx);
-
- Scalar MAW = (fluidState.moleFraction(gPhaseIdx, airIdx)*Air::molarMass()
- + fluidState.moleFraction(gPhaseIdx, H2OIdx)*H2O::molarMass())
- / xAW;
-
- Scalar phiCAW = 0.3; // simplification for this particular system
- /* actually like this
- * Scalar phiCAW = std::pow(1.+std::sqrt(mu[NAPLIdx]/muAW)*std::pow(MAW/M[NAPLIdx],0.25),2)
- * / std::sqrt(8.*(1.+M[NAPLIdx]/MAW));
- */
- Scalar phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW);
-
- muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gPhaseIdx, NAPLIdx)*phiAWC)
- + (fluidState.moleFraction(gPhaseIdx, NAPLIdx) * mu[NAPLIdx])
- / (fluidState.moleFraction(gPhaseIdx, NAPLIdx) + xAW*phiCAW);
- return muResult;
- }
- DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ assert (phaseIdx == gPhaseIdx);
+
+ /* Wilke method. See:
+ *
+ * See: R. Reid, et al.: The Properties of Gases and Liquids,
+ * 4th edition, McGraw-Hill, 1987, 407-410
+ * 5th edition, McGraw-Hill, 20001, p. 9.21/22
+ *
+ * in this case, we use a simplified version in order to avoid
+ * computationally costly evaluation of sqrt and pow functions and
+ * divisions
+ * -- compare e.g. with Promo Class p. 32/33
+ */
+ Scalar muResult;
+ const Scalar mu[numComponents] = {
+ H2O::gasViscosity(fluidState.temperature(phaseIdx), H2O::vaporPressure(fluidState.temperature(phaseIdx))),
+ Air::simpleGasViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)),
+ NAPL::gasViscosity(fluidState.temperature(phaseIdx), NAPL::vaporPressure(fluidState.temperature(phaseIdx)))
+ };
+ // molar masses
+ const Scalar M[numComponents] = {
+ H2O::molarMass(),
+ Air::molarMass(),
+ NAPL::molarMass()
+ };
+
+ Scalar muAW = mu[airIdx]*fluidState.moleFraction(gPhaseIdx, airIdx)
+ + mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
+ / (fluidState.moleFraction(gPhaseIdx, airIdx)
+ + fluidState.moleFraction(gPhaseIdx, H2OIdx));
+ Scalar xAW = fluidState.moleFraction(gPhaseIdx, airIdx)
+ + fluidState.moleFraction(gPhaseIdx, H2OIdx);
+
+ Scalar MAW = (fluidState.moleFraction(gPhaseIdx, airIdx)*Air::molarMass()
+ + fluidState.moleFraction(gPhaseIdx, H2OIdx)*H2O::molarMass())
+ / xAW;
+
+ Scalar phiCAW = 0.3; // simplification for this particular system
+ /* actually like this
+ * Scalar phiCAW = std::pow(1.+std::sqrt(mu[NAPLIdx]/muAW)*std::pow(MAW/M[NAPLIdx],0.25),2)
+ * / std::sqrt(8.*(1.+M[NAPLIdx]/MAW));
+ */
+ Scalar phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW);
+
+ muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gPhaseIdx, NAPLIdx)*phiAWC)
+ + (fluidState.moleFraction(gPhaseIdx, NAPLIdx) * mu[NAPLIdx])
+ / (fluidState.moleFraction(gPhaseIdx, NAPLIdx) + xAW*phiCAW);
+ return muResult;
}
diff --git a/test/boxmodels/3p3c/infiltrationproblem.hh b/test/boxmodels/3p3c/infiltrationproblem.hh
index 968f6c526fbd663fd75a43705afa6b7c9f78f18c..eae44df52527ef9f95e9fcd0a658990b9a96479d 100644
--- a/test/boxmodels/3p3c/infiltrationproblem.hh
+++ b/test/boxmodels/3p3c/infiltrationproblem.hh
@@ -69,13 +69,12 @@ SET_TYPE_PROP(InfiltrationProblem,
SET_BOOL_PROP(InfiltrationProblem, EnableGravity, true);
// Write newton convergence?
-SET_BOOL_PROP(InfiltrationProblem, NewtonWriteConvergence, false);
+SET_BOOL_PROP(InfiltrationProblem, NewtonWriteConvergence, true);
-//! Set the formulation
-SET_INT_PROP(InfiltrationProblem, Formulation, ThreePThreeCFormulation::pgSwSn);
+// Maximum tolerated relative error in the Newton method
+SET_SCALAR_PROP(InfiltrationProblem, NewtonRelTolerance, 1e-8);
}
-
/*!
* \ingroup ThreePThreeCBoxModel
*
@@ -83,16 +82,14 @@ SET_INT_PROP(InfiltrationProblem, Formulation, ThreePThreeCFormulation::pgSwSn);
template <class TypeTag >
class InfiltrationProblem : public ThreePThreeCProblem<TypeTag>
{
+ typedef ThreePThreeCProblem<TypeTag> ParentType;
+
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::Grid Grid;
-
- typedef ThreePThreeCProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
-
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, ThreePThreeCIndices) Indices;
enum {