From 3afbdf61bb0efea5e0acbb3d83cb5136847c90d1 Mon Sep 17 00:00:00 2001
From: Gabriele Seitz <gabriele.seitz@iws.uni-stuttgart.de>
Date: Mon, 3 Jul 2017 16:24:07 +0200
Subject: [PATCH] [1pncmin][problem] minor improvements
---
.../1pncmin/implicit/test_1pncmin.input | 12 +--
.../1pncmin/implicit/thermochemproblem.hh | 76 ++++++-------------
.../implicit/thermochemspatialparams.hh | 14 +++-
3 files changed, 44 insertions(+), 58 deletions(-)
diff --git a/test/porousmediumflow/1pncmin/implicit/test_1pncmin.input b/test/porousmediumflow/1pncmin/implicit/test_1pncmin.input
index 51c2ebb610..67aee3552c 100644
--- a/test/porousmediumflow/1pncmin/implicit/test_1pncmin.input
+++ b/test/porousmediumflow/1pncmin/implicit/test_1pncmin.input
@@ -1,6 +1,6 @@
-# Parameter file for test case 2pnc.
+# Parameter file for test case 1pncmin.
# Everything behind a '#' is a comment.
-# Type "./test_2pnc --help" for more information.
+# Type "./test_1pncmin --help" for more information.
[TimeManager]
DtInitial = 1 # [s] initial time step size
@@ -23,7 +23,7 @@ TemperatureHigh = 873.15 # [Pa]high end for tabularization of fluid p
[Problem]
Name = thermochem
-IsCharge = 0 # Bool: 1: charge; 0: discharge
+IsCharge = 1 # Bool: 1: charge; 0: discharge
[Charge]
PressureInitial = 1E5 # [Pa] Initial reservoir pressure
@@ -34,7 +34,7 @@ CaO2H2Initial = 0.3960 # [-] molefraction in the solid phase; 0 de
PressureIn = 1.02e5 # [Pa] Inlet pressure; charge
PressureOut = 1e5 # [Pa] outlet pressure
TemperatureIn = 773.15 # [K] inlet temperature: // Shao 500 °C
-TemperatureOut = 773.15 # [K] outlet temperature: // Shao noflow
+TemperatureOut = 673.15 # [K] outlet temperature: // Shao noflow
InFlow = 5 # [mol/s] Inflow of HTF // Shao 0.309 g/s (Area (0.015)^2pi m^2) --> here 1.55 mol/s
VaporIn = 0.01 # [] molefraction // Shao 0.01 [g/g]
@@ -52,7 +52,7 @@ InFlow = 5 # 0.277[mol/s] Inflow of HTF
VaporIn = 0.36 # [] molefraction
[Vtk]
-AddVelocity = 1 # Add extra information
+AddVelocity = 0 # Add extra information
VtuWritingFreq = 1 # 1: write a vtu file at every timestep, 2: write a vtu file every second timestep ...
[LinearSolver]
@@ -68,4 +68,4 @@ FreqRestart = 1000 # how often restart files are written out
FreqOutput = 50 # frequency of VTK output
FreqMassOutput = 2 # frequency of mass and evaporation rate output (Darcy)
FreqFluxOutput = 1000 # frequency of detailed flux output
-FreqVaporFluxOutput = 2 # frequency of summarized flux output
\ No newline at end of file
+FreqVaporFluxOutput = 2 # frequency of summarized flux output
diff --git a/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh b/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
index 6865561e5b..366843be44 100644
--- a/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
+++ b/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
@@ -65,6 +65,9 @@ SET_PROP(ThermoChemProblem, FluidSystem)
// Set the transport equation that is replaced by the total mass balance
// SET_INT_PROP(ThermoChemProblem, ReplaceCompEqIdx, 1 /*3*/ /*1*/);
+// Enable velocity output
+SET_BOOL_PROP(ThermoChemProblem, VtkAddVelocity, false);
+
SET_TYPE_PROP(ThermoChemProblem, LinearSolver, UMFPackBackend<TypeTag>);
@@ -137,9 +140,6 @@ class ThermoChemProblem : public ImplicitPorousMediaProblem<TypeTag>
using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
using DimVector = Dune::FieldVector<Scalar, dim> ;
-
- // Select the electrochemistry method
-// typedef typename Dumux::ElectroChemistry<TypeTag, Dumux::ElectroChemistryModel::Ochs> ElectroChemistry;
// typedef Dumux::Constants<Scalar> Constant;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
@@ -219,11 +219,8 @@ public:
if (globalPos[0] > this->bBoxMax()[0] - eps_){
values.setDirichlet(pressureIdx);
-// values.setDirichlet(firstMoleFracIdx);
-// values.setDirichlet(temperatureIdx);
- values.setOutflow(temperatureIdx);
- values.setOutflow(firstMoleFracIdx);
-// values.setDirichlet(firstMoleFracIdx);
+ values.setDirichlet(firstMoleFracIdx);
+ values.setDirichlet(temperatureIdx);
}
return values;
}
@@ -241,42 +238,42 @@ public:
PrimaryVariables priVars(0.0);
//input parameters
- Scalar pIn;
+// Scalar pIn;
Scalar pOut;
Scalar tIn;
- Scalar tOut;
- Scalar vapor;
+ Scalar tOut;
+ Scalar vaporIn;
// read input parameters
if (isCharge_ == true){
- pIn = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, PressureIn);
+// pIn = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, PressureIn);
pOut = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, PressureOut);
tIn = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, TemperatureIn);
tOut = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, TemperatureOut);
- vapor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, VaporIn);
+ vaporIn = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, VaporIn);
}
else{
- pIn = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, PressureIn);
+// pIn = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, PressureIn);
pOut = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, PressureOut);
tIn = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, TemperatureIn);
tOut = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, TemperatureOut);
- vapor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, VaporIn);
+ vaporIn = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, VaporIn);
}
if(globalPos[0] < eps_)
{
-// priVars[pressureIdx] = pIn;
- priVars[firstMoleFracIdx] = vapor; // Saturation outer boundary
+ priVars[firstMoleFracIdx] = vaporIn; // Saturation outer boundary
priVars[temperatureIdx] = tIn;
}
if(globalPos[0] > this->bBoxMax()[0] - eps_)
{
priVars[pressureIdx] = pOut;
-// priVars[firstMoleFracIdx] = 0.01; // Saturation inner boundary
-// priVars[temperatureIdx] = tOut;
-// priVars[firstMoleFracIdx] = 0;
+ priVars[firstMoleFracIdx] = 0.01; // Saturation inner boundary
+ priVars[temperatureIdx] = tOut;
+
}
+
return priVars;
}
@@ -312,7 +309,7 @@ public:
else
priVars[pressureIdx] = -GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, InFlow); //[mol/s] gas inflow
}
-
+// std::cout << " test neumann " << "\n";
return priVars;
}
@@ -338,8 +335,6 @@ public:
Scalar CaO2H2Init;
if (isCharge_ == true){
-
- std::cout << "true " << "\n";
pInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, PressureInitial);
tInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, TemperatureInitial);
h2oInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, VaporInitial);
@@ -355,7 +350,7 @@ public:
CaOInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, CaOInitial);
CaO2H2Init = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, CaO2H2Initial);
}
-
+ std::cout << "CaO2H2Init = " << CaO2H2Init << "\n";
priVars[pressureIdx] = pInit;
priVars[firstMoleFracIdx] = h2oInit;
#if NONISOTHERMAL
@@ -374,7 +369,7 @@ public:
* \param fvGeometry The finite volume geometry
* \param scvIdx The sub control volume index
*/
- PrimaryVariables solDependentSource(const Element &element,
+ PrimaryVariables source(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
@@ -383,29 +378,15 @@ public:
PrimaryVariables source(0.0);
const auto& volVars = elemVolVars[scv];
- Scalar Initial_Precipitate_Volume;
- Scalar maxPrecipitateVolumeCaO;
- Scalar maxPrecipitateVolumeCaO2H2;
-
- if (isCharge_ == true){
- Initial_Precipitate_Volume = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, CaO2H2Initial);
- maxPrecipitateVolumeCaO = 0.3960;
- maxPrecipitateVolumeCaO2H2 = Initial_Precipitate_Volume;
- }
- else {
- Initial_Precipitate_Volume = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, CaOInitial);
- maxPrecipitateVolumeCaO = Initial_Precipitate_Volume;
- maxPrecipitateVolumeCaO2H2 = 0.3960;
- }
Scalar T= volVars.temperature();
Scalar Teq = 0;
Scalar moleFractionVapor = 1e-3;
- if(volVars.moleFraction(firstMoleFracIdx) > 1e-3)
- moleFractionVapor = volVars.moleFraction(firstMoleFracIdx) ;
- if(volVars.moleFraction(firstMoleFracIdx) >= 1.0){
+ if(volVars.moleFraction(phaseIdx, firstMoleFracIdx) > 1e-3)
+ moleFractionVapor = volVars.moleFraction(phaseIdx, firstMoleFracIdx) ;
+ if(volVars.moleFraction(phaseIdx, firstMoleFracIdx) >= 1.0){
moleFractionVapor = 1;
// std::cout << " test vapor = " << "\n";
}
@@ -416,14 +397,7 @@ public:
Teq = -12845;
Teq /= (pFactor - 16.508); //the equilibrium temperature
-// if(isCharge == false) {
-// if (Teq < 573.15) {
-// std::cout << "Teq = " << Teq<< "\n";
-// // Teq=573.15;
-// }
-// }
-
- Scalar moleFracH2O_fPhase = volVars.moleFraction(firstMoleFracIdx);
+ Scalar moleFracH2O_fPhase = volVars.moleFraction(phaseIdx, firstMoleFracIdx);
Scalar moleFracCaO_sPhase = volVars.precipitateVolumeFraction(cPhaseIdx)*volVars.molarDensity(cPhaseIdx)
/(volVars.precipitateVolumeFraction(hPhaseIdx)*volVars.molarDensity(hPhaseIdx)
@@ -478,7 +452,7 @@ public:
if (- q*this->timeManager().timeStepSize() + moleFracCaO2H2_sPhase*volVars.molarDensity(hPhaseIdx) < 0){
q = moleFracCaO2H2_sPhase/this->timeManager().timeStepSize();
- // std::cout << "q_charge " << q << "\n";
+
}
source[conti0EqIdx+CaO2H2Idx] = -q;
diff --git a/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh b/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
index c9098f281d..b11543f55a 100644
--- a/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
+++ b/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
@@ -191,6 +191,18 @@ public:
return phi;
}
+ /*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
+ * on the position in the domain
+ *
+ * \param element The finite volume element
+ * \param scv The sub-control volume
+ *
+ * Solution dependent permeability function
+ */
+ Scalar permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return permLaw_.evaluatePermeability(element, scv, elemSol); }
/*!
* \brief Define the minimum porosity \f$[-]\f$ distribution
@@ -273,7 +285,7 @@ private:
Scalar eps_;
// MaterialLawParams materialParams_;
Scalar lambdaSolid_;
- bool isCharge_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Problem, isCharge);
+ bool isCharge_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Problem, IsCharge);
PorosityLaw poroLaw_;
PermeabilityLaw permLaw_;
};
--
GitLab