diff --git a/test/porousmediumflow/1pncmin/CMakeLists.txt b/test/porousmediumflow/1pncmin/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ba8341c614f1a2c797c95f5402f602025f1087b1
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory("implicit")
diff --git a/test/porousmediumflow/1pncmin/implicit/CMakeLists.txt b/test/porousmediumflow/1pncmin/implicit/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..941775ea72999bc0e5f2c5d2e4191e41ebd86216
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/implicit/CMakeLists.txt
@@ -0,0 +1,50 @@
+# isothermal tests
+
+add_input_file_links()
+
+add_executable(test_box2pnc test_box2pnc.cc)
+
+#add_executable_all(test_box1p2c test_box1p2c.cc)
+#
+# add_executable_all(test_cc1p2c test_cc1p2c.cc)
+
+# dune_symlink_to_source_files(FILES test_2pnc.input)
+# # dune_symlink_to_source_files(FILES test_box2pnc.input)
+# dune_symlink_to_source_files(FILES test_box1p2c.input)
+# dune_symlink_to_source_files(FILES test_cc1p2c.input)
+
+# set(BUILD_TYPE Debug)
+
+#install sources
+# install(FILES
+# thermochemproblem.hh
+# thermochemspatialparams.hh
+# test_box2pnc.cc
+# test_cc2pnc.cc
+# DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux-devel/appl/thermochemistry)
+
+
+# # isothermal tests
+# add_dumux_test(test_box2pnc test_box2pnc test_box2pnc.cc
+# python ${dumux_INCLUDE_DIRS}/bin/testing/runtest.py
+# --script fuzzy
+# --files ${CMAKE_SOURCE_DIR}/test/references/fuelcell2pncbox-reference.vtu
+# ${CMAKE_CURRENT_BINARY_DIR}/fuelcell-box-00016.vtu
+# --command "${CMAKE_CURRENT_BINARY_DIR}/test_box2pnc -ParameterFile test_2pnc.input -Problem.Name fuelcell-box")
+#
+# add_dumux_test(test_cc2pnc test_cc2pnc test_cc2pnc.cc
+# python ${dumux_INCLUDE_DIRS}/bin/testing/runtest.py
+# --script fuzzy
+# --files ${CMAKE_SOURCE_DIR}/test/references/fuelcell2pnccc-reference.vtu
+# ${CMAKE_CURRENT_BINARY_DIR}/fuelcell-cc-00016.vtu
+# --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc2pnc -ParameterFile test_2pnc.input -Problem.Name fuelcell-cc")
+#
+# dune_symlink_to_source_files(FILES test_2pnc.input)
+#
+# #install sources
+# install(FILES
+# thermochemproblem.hh
+# thermochemspatialparams.hh
+# test_box2pnc.cc
+# test_cc2pnc.cc
+# DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/2pnc)
diff --git a/test/porousmediumflow/1pncmin/implicit/test_2pnc.input b/test/porousmediumflow/1pncmin/implicit/test_2pnc.input
new file mode 100644
index 0000000000000000000000000000000000000000..51c2ebb61053a6d873964622e52a9806560a3d89
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/implicit/test_2pnc.input
@@ -0,0 +1,71 @@
+# Parameter file for test case 2pnc.
+# Everything behind a '#' is a comment.
+# Type "./test_2pnc --help" for more information.
+
+[TimeManager]
+DtInitial = 1 # [s] initial time step size
+MaxTimeStepSize = 50 # [s] maximum time step size
+TEnd= 5000 # [s] duration of the simulation
+FreqOutput = 10 # frequency of VTK output
+WriteRestartFile = 1 # Boolean. Should restart files be written? (1) Yes (0) No
+
+[Grid]
+UpperRight = 0.08 0.01 #20 10 # # [m] upper right corner coordinates
+Cells = 40 2 # 21 6 # [-] number of cells in x,y-direction
+
+[FluidSystem]
+NTemperature = 10 # [-] number of tabularization entries
+NPressure = 100 # [-] number of tabularization entries
+PressureLow = 1E5 # [Pa]low end for tabularization of fluid properties
+PressureHigh = 1E6 # [Pa]high end for tabularization of fluid properties
+TemperatureLow = 373.15 # [Pa]low end for tabularization of fluid properties
+TemperatureHigh = 873.15 # [Pa]high end for tabularization of fluid properties
+
+[Problem]
+Name = thermochem
+IsCharge = 0 # Bool: 1: charge; 0: discharge
+
+[Charge]
+PressureInitial = 1E5 # [Pa] Initial reservoir pressure
+TemperatureInitial = 773.15 # [K] reservoir temperature
+VaporInitial = 0.01 # [-] initial mole fraction of water
+CaOInitial = 0.0 # [-] molefraction in the solid phase; 0 dehydration/charge, 1 hydration/discharge
+CaO2H2Initial = 0.3960 # [-] molefraction in the solid phase; 0 dehydration/charge, 1 hydration/discharge
+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
+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]
+
+[Discharge]
+PressureInitial = 1E5 # [Pa] Initial reservoir pressure
+TemperatureInitial = 573.15 # [K] reservoir temperature
+VaporInitial = 0.35 # [-] initial mole fraction of water
+CaOInitial = 0.2 # [-] molefraction in the solid phase; 0 dehydration/charge, 1 hydration/discharge
+CaO2H2Initial = 0.0
+PressureIn = 1.05e5 # [Pa] Inlet pressure; charge
+PressureOut = 1e5 # [Pa] outlet pressure
+TemperatureIn = 573.15 # [K] inlet temperature: charge: 873 K ; discharge: 473K
+TemperatureOut = 573.15 # [K] outlet temperature: charge: 573; discharge: outflow
+InFlow = 5 # 0.277[mol/s] Inflow of HTF
+VaporIn = 0.36 # [] molefraction
+
+[Vtk]
+AddVelocity = 1 # Add extra information
+VtuWritingFreq = 1 # 1: write a vtu file at every timestep, 2: write a vtu file every second timestep ...
+
+[LinearSolver]
+ResidualReduction = 1e-6
+
+[Newton]
+MaxRelativeShift = 1e-6
+2WriteConvergence = 1
+
+[Output]
+#Frequency of restart file, flux and VTK output
+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
diff --git a/test/porousmediumflow/1pncmin/implicit/test_box2pnc.cc b/test/porousmediumflow/1pncmin/implicit/test_box2pnc.cc
new file mode 100644
index 0000000000000000000000000000000000000000..fc844ca2b91ec1d33791405d6bc8cc3c72b0294f
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/implicit/test_box2pnc.cc
@@ -0,0 +1,55 @@
+// -*- 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 Test for the 2pnc box model used for water management in PEM fuel cells.
+ */
+#include <config.h>
+#include "thermochemproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+ "\t-ParameterFile Parameter file (Input file) \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(ThermoChemBoxProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh b/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..585a6df9330f4a35347b9d589ee1d7267d3c6f04
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
@@ -0,0 +1,571 @@
+// -*- 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 Definition of a problem for water management in PEM fuel cells.
+ */
+#ifndef DUMUX_THERMOCHEM_PROBLEM_HH
+#define DUMUX_THERMOCHEM_PROBLEM_HH
+
+#include <dumux/porousmediumflow/1pncmin/implicit/model.hh>
+#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/material/fluidsystems/simplesteamaircao2h2.hh>
+#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
+
+#include "thermochemspatialparams.hh"
+
+#define NONISOTHERMAL 1
+
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class ThermoChemProblem;
+
+namespace Properties
+{
+#if NONISOTHERMAL
+NEW_TYPE_TAG(ThermoChemProblem, INHERITS_FROM(OnePNCMinNI, ThermoChemSpatialParams));
+NEW_TYPE_TAG(ThermoChemBoxProblem, INHERITS_FROM(BoxModel, ThermoChemProblem));
+NEW_TYPE_TAG(ThermoChemCCProblem, INHERITS_FROM(CCModel, ThermoChemProblem));
+#else
+NEW_TYPE_TAG(ThermoChemProblem, INHERITS_FROM(OnePNCMin, ThermoChemSpatialParams));
+NEW_TYPE_TAG(ThermoChemBoxProblem, INHERITS_FROM(BoxModel, ThermoChemProblem));
+NEW_TYPE_TAG(ThermoChemCCProblem, INHERITS_FROM(CCModel, ThermoChemProblem));
+#endif
+// Set the grid type
+SET_TYPE_PROP(ThermoChemProblem, Grid, Dune::YaspGrid<2>);
+// Set the problem property
+SET_TYPE_PROP(ThermoChemProblem, Problem, ThermoChemProblem<TypeTag>);
+// Set fluid configuration
+SET_PROP(ThermoChemProblem, FluidSystem)
+{ /*private:*/
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::SteamAirCaO2H2<Scalar>;
+};
+
+// Set the transport equation that is replaced by the total mass balance
+// SET_INT_PROP(ThermoChemProblem, ReplaceCompEqIdx, 1 /*3*/ /*1*/);
+
+SET_TYPE_PROP(ThermoChemProblem, LinearSolver, UMFPackBackend<TypeTag>);
+
+}
+
+// Set the spatial parameters
+SET_TYPE_PROP(ThermoChemProblem, SpatialParams, ThermoChemSpatialparams<TypeTag>);
+
+/*!
+ * \ingroup OnePNCModel
+ * \ingroup ImplicitTestProblems
+ * \brief Problem or water management in PEM fuel cells.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_box2pnc</tt>
+ */
+template <class TypeTag>
+class ThermoChemProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+ using ParentType = ImplicitPorousMediaProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+
+ enum { dim = GridView::dimension };
+ enum { dimWorld = GridView::dimensionworld };
+
+ enum
+ {
+ replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx),
+
+ numComponents = FluidSystem::numComponents,
+ numSComponents = FluidSystem::numSComponents,
+
+ // Indices of the primary variables
+ pressureIdx = Indices::pressureIdx, //gas-phase pressure
+ firstMoleFracIdx = Indices::firstMoleFracIdx, // mole fraction water
+
+ CaOIdx = FluidSystem::numComponents,
+ CaO2H2Idx = FluidSystem::numComponents+1,
+
+ //Equation Indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ firstTransportEqIdx = Indices::firstTransportEqIdx,
+ solidEqIdx = Indices::conti0EqIdx + FluidSystem::numComponents,
+
+ // Phase Indices
+ phaseIdx = FluidSystem::gPhaseIdx,
+ cPhaseIdx = FluidSystem::cPhaseIdx,
+ hPhaseIdx = FluidSystem::hPhaseIdx,
+
+#if NONISOTHERMAL
+ temperatureIdx = Indices::temperatureIdx,
+ energyEqIdx = Indices::energyEqIdx
+#endif
+ };
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Vertex = typename GridView::template Codim<dim>::Entity;
+ using Intersection = typename GridView::Intersection;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+ 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) };
+ enum { dofCodim = isBox ? dim : 0 };
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ ThermoChemProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ nTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, FluidSystem, NTemperature);
+ nPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, FluidSystem, NPressure);
+ pressureLow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, PressureLow);
+ pressureHigh_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, PressureHigh);
+ temperatureLow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, TemperatureLow);
+ temperatureHigh_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FluidSystem, TemperatureHigh);
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+
+
+ FluidSystem::init(/*Tmin=*/temperatureLow_,
+ /*Tmax=*/temperatureHigh_,
+ /*nT=*/nTemperature_,
+ /*pmin=*/pressureLow_,
+ /*pmax=*/pressureHigh_,
+ /*np=*/nPressure_);
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const std::string name() const
+ { return name_; }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return temperature_; }
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment
+ *
+ * \param globalPos The global position
+ */
+ BoundaryTypes boundaryTypesAtPos( const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes values;
+
+ values.setAllNeumann();
+
+ if(globalPos[0] < eps_ )
+ {
+// values.setDirichlet(pressureIdx);
+ values.setDirichlet(firstMoleFracIdx);
+ values.setDirichlet(temperatureIdx);
+ }
+
+ 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);
+ }
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet
+ * boundary segment
+ *
+ * \param values Stores the Dirichlet values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variable} ] \f$
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables priVars(0.0);
+
+ //input parameters
+ Scalar pIn;
+ Scalar pOut;
+ Scalar tIn;
+ Scalar tOut;
+ Scalar vapor;
+
+ // read input parameters
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, IsCharge)){
+ 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);
+
+ }
+ else{
+ 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);
+ }
+
+ if(globalPos[0] < eps_)
+ {
+// priVars[pressureIdx] = pIn;
+ priVars[firstMoleFracIdx] = vapor; // 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;
+ }
+ return priVars;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Neumann
+ * boundary segment in dependency on the current solution.
+ *
+ * \param values Stores the Neumann values for the conservation equations in
+ * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param intersection The intersection between element and boundary
+ * \param scvIdx The local index of the sub-control volume
+ * \param boundaryFaceIdx The index of the boundary face
+ * \param elemVolVars All volume variables for the element
+ *
+ * This method is used for cases, when the Neumann condition depends on the
+ * solution and requires some quantities that are specific to the fully-implicit method.
+ * The \a values store the mass flux of each phase normal to the boundary.
+ * Negative values indicate an inflow.
+ */
+
+ PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables priVars(0.0);
+
+ if(globalPos[0] < eps_)
+ {
+ //if(isCharge == true){
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, IsCharge)){
+ priVars[pressureIdx] = -GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, InFlow); //[mol/s] gas inflow; negative sign: inflow
+ }
+ else
+ priVars[pressureIdx] = -GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, InFlow); //[mol/s] gas inflow
+ }
+
+ return priVars;
+ }
+
+ /*!
+ * \name Volume terms
+ */
+
+ /*!
+ * \brief Evaluates the initial values for a control volume
+ *
+ * \param values Stores the initial values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variables} ] \f$
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables priVars(0.0);
+
+ Scalar pInit;
+ Scalar tInit;
+ Scalar h2oInit;
+ Scalar CaOInit;
+ Scalar CaO2H2Init;
+
+ //if(isCharge == true){
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, IsCharge)){
+ 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);
+ CaOInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, CaOInitial);
+ CaO2H2Init = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Charge, CaO2H2Initial);
+ }
+
+ else {
+ std::cout << "false " << "\n";
+ pInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, PressureInitial);
+ tInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, TemperatureInitial);
+ h2oInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, VaporInitial);
+ CaOInit = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, CaOInitial);
+ CaO2H2Init = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Discharge, CaO2H2Initial);
+ }
+
+ priVars[pressureIdx] = pInit;
+ priVars[firstMoleFracIdx] = h2oInit;
+#if NONISOTHERMAL
+ priVars[temperatureIdx] = tInit;
+#endif
+ priVars[CaOIdx] = CaOInit;
+ priVars[CaO2H2Idx] = CaO2H2Init;
+
+ return priVars;
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a sub control volume.
+ *
+ * \param element The element of the sub control volume
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The sub control volume index
+ */
+ PrimaryVariables solDependentSource(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+
+ PrimaryVariables source(0.0);
+ const auto& volVars = elemVolVars[scv];
+
+ Scalar Initial_Precipitate_Volume;
+ Scalar maxPrecipitateVolumeCaO;
+ Scalar maxPrecipitateVolumeCaO2H2;
+
+ //if(isCharge == true){
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, IsCharge)){
+ 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){
+ moleFractionVapor = 1;
+// std::cout << " test vapor = " << "\n";
+ }
+ Scalar vaporPressure = volVars.pressure(phaseIdx) *moleFractionVapor ;
+ vaporPressure *= 1.0e-5;
+ Scalar pFactor = log(vaporPressure);
+
+ 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 moleFracCaO_sPhase = volVars.precipitateVolumeFraction(cPhaseIdx)*volVars.molarDensity(cPhaseIdx)
+ /(volVars.precipitateVolumeFraction(hPhaseIdx)*volVars.molarDensity(hPhaseIdx)
+ + volVars.precipitateVolumeFraction(cPhaseIdx)*volVars.molarDensity(cPhaseIdx));
+
+ Scalar moleFracCaO2H2_sPhase = volVars.precipitateVolumeFraction(hPhaseIdx)*volVars.molarDensity(hPhaseIdx)
+ /(volVars.precipitateVolumeFraction(hPhaseIdx)*volVars.molarDensity(hPhaseIdx)
+ + volVars.precipitateVolumeFraction(cPhaseIdx)*volVars.molarDensity(cPhaseIdx));
+
+ Scalar deltaH = 112e3; // J/mol
+
+ Scalar solidDensityAverage = moleFracCaO_sPhase*volVars.molarDensity(cPhaseIdx)
+ + moleFracCaO2H2_sPhase* volVars.molarDensity(hPhaseIdx);
+
+// //discharge or hydration
+ if (T < Teq){
+
+ Scalar krh =0.08; //0.006
+
+ Scalar rHydration = - moleFracH2O_fPhase* (volVars.molarDensity(hPhaseIdx)- solidDensityAverage)
+ * krh * (T-Teq)/ Teq;
+
+ Scalar q = - rHydration ;
+
+ // make sure not more CaO reacts than present
+
+ if (- q*this->timeManager().timeStepSize() + moleFracCaO_sPhase* volVars.molarDensity(cPhaseIdx) < 0 + eps_){
+ q = moleFracCaO_sPhase/this->timeManager().timeStepSize();
+// std::cout << "q_discharge = " << q << "\n";
+ }
+
+ source[conti0EqIdx+CaO2H2Idx] = q;
+
+ source[conti0EqIdx+CaOIdx] = -q;
+
+ source[conti0EqIdx+firstMoleFracIdx] = -q;
+
+#if NONISOTHERMAL
+ source[energyEqIdx] = q * deltaH;
+#endif
+ }
+
+ // charge or dehydration
+ else if(T > Teq){
+
+ Scalar krd = 0.03; //0.05;
+
+ Scalar rDehydration = (volVars.molarDensity(cPhaseIdx)- solidDensityAverage)
+ * krd * (Teq-T)/ Teq;
+
+ Scalar q = -rDehydration;
+
+ 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;
+
+ source[conti0EqIdx+CaOIdx] = q;
+
+ source[conti0EqIdx+firstMoleFracIdx] = q;
+#if NONISOTHERMAL
+ source[energyEqIdx] = -q * deltaH;
+
+#endif
+ }
+
+ else
+ source = 0.0;
+
+ return source;
+ }
+
+
+ /*!
+ * \brief Add problem specific vtk output for the electrochemistry
+ */
+// void addOutputVtkFields()
+// {
+// // add the output field specific to the electrochemistry
+// typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+//
+// // get the number of degrees of freedom
+// // auto numDofs = this->model().numDofs();
+//
+// for (const auto& element : elements(this->gridView()))
+// {
+// FVElementGeometry fvGeometry;
+// fvGeometry.update(this->gridView(), element);
+//
+// ElementVolumeVariables elemVolVars;
+// elemVolVars.update(*this,
+// element,
+// fvGeometry,
+// false /* oldSol? */);
+//
+// // for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+// // {
+// // const auto& globalPos = isBox ? element.geometry().corner(scvIdx)
+// // : element.geometry().center();
+// //
+// // auto dofIdxGlobal = this->model().dofMapper().subIndex(element, scvIdx, dofCodim);
+// // }
+// }
+// }
+
+private:
+
+ static Scalar massTomoleFrac_(Scalar XlNaCl)
+ {
+ const Scalar Mw = 18.015e-3; /* molecular weight of water [kg/mol] */
+ const Scalar Ms = 58.44e-3; /* molecular weight of NaCl [kg/mol] */
+
+ const Scalar X_NaCl = XlNaCl;
+ /* XlNaCl: conversion from mass fraction to mol fraction */
+ auto xlNaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
+ return xlNaCl;
+ }
+
+ int nTemperature_;
+ int nPressure_;
+
+// PrimaryVariables storageLastTimestep_;
+
+ std::string name_;
+
+ Scalar pressureLow_, pressureHigh_;
+ Scalar temperatureLow_, temperatureHigh_;
+ Scalar reservoirPressure_;
+ Scalar innerPressure_;
+ Scalar outerPressure_;
+ Scalar temperature_;
+
+ static constexpr Scalar eps_ = 1e-6;
+ Scalar reservoirSaturation_;
+ std::ofstream outfile;
+
+// bool isCharge_ ;
+
+};
+
+} //end namespace
+
+#endif
diff --git a/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh b/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c9098f281d1b0c371e1cb1e53d3ab326b8ed2067
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
@@ -0,0 +1,283 @@
+// -*- 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 Definition of the spatial parameters for the fuel cell
+ * problem which uses the isothermal/non-insothermal 2pnc box model
+ */
+
+#ifndef DUMUX_THERMOCHEM_SPATIAL_PARAMS_HH
+#define DUMUX_THERMOCHEM_SPATIAL_PARAMS_HH
+
+#include <dumux/material/spatialparams/implicit1p.hh>
+#include <dumux/porousmediumflow/1pncmin/implicit/indices.hh>
+// #include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+// #include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+// #include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+// #include <dumux/material/fluidmatrixinteractions/2p/philtophoblaw.hh>
+#include <dumux/material/fluidmatrixinteractions/porosityreactivebed.hh>
+#include <dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh>
+
+namespace Dumux
+{
+
+//forward declaration
+template<class TypeTag>
+class ThermoChemSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(ThermoChemSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(ThermoChemSpatialParams, SpatialParams, Dumux::ThermoChemSpatialParams<TypeTag>);
+
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \ingroup BoxTestProblems
+ * \brief Definition of the spatial parameters for the FuelCell
+ * problem which uses the isothermal 2p2c box model
+ */
+template<class TypeTag>
+class ThermoChemSpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
+{
+ using ThisType = ThermoChemSpatialParams<TypeTag>;
+ using ParentType = ImplicitSpatialParamsOneP<TypeTag>;
+// typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using CoordScalar = typename GridView::ctype;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+
+// wPhaseIdx = FluidSystem::wPhaseIdx
+ };
+
+ using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
+// typedef Dune::FieldVector<CoordScalar,dim> DimVector;
+// typedef Dune::FieldMatrix<CoordScalar,dim,dim> DimMatrix;
+ using Tensor = Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld>;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+// typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using PorosityLaw = PorosityReactiveBed<TypeTag>;
+ using PermeabilityLaw = PermeabilityKozenyCarman<TypeTag>;
+
+public:
+ // type used for the permeability (i.e. tensor or scalar)
+ using PermeabilityType = Scalar;
+ /*!
+ * \brief The constructor
+ *
+ * \param gridView The grid view
+ */
+ ThermoChemSpatialParams(const Problem& problem, const GridView &gridView)
+ : ParentType(problem, gridView)
+ {
+ // intrinsic permeabilities
+// K_[0][0] = 5e-12;
+// K_[1][1] = 5e-12;
+// K_[0][0] = 2.23e-14;
+// K_[1][1] = 2.23e-14;
+
+ //thermal conductivity of CaO
+ lambdaSolid_ = 0.4; //[W/(m*K)] Nagel et al [2013b]
+
+ eps_ = 1e-6;
+ }
+
+ ~ThermoChemSpatialParams()
+ {}
+
+ /*!
+ * \brief Called by the Problem to initialize the spatial params.
+ */
+ void init()
+ {
+ //! Intitialize the parameter laws
+ poroLaw_.init(*this);
+ permLaw_.init(*this);
+ }
+
+ /*! Old
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvGeometry The current finite volume geometry of the element
+ * \param scvIdx The index of the sub-control volume
+ */
+// const DimMatrix intrinsicPermeability(const Element &element,
+// const FVElementGeometry &fvGeometry,
+// const int scvIdx) const
+// { return K_; }
+//
+ /*!
+ * \brief Define the initial permeability \f$[m^2]\f$ distribution
+ *
+ * \param element The finite element
+ * \param scv The sub-control volume
+ */
+ Scalar initialPermeability(const Element& element, const SubControlVolume &scv) const
+ { return 5e-12; }
+
+ /*!
+ * \brief Define the porosity \f$[-]\f$ of the spatial parameters
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the porosity needs to be defined
+ */
+// Scalar porosity(const Element &element,
+// const FVElementGeometry &fvGeometry,
+// const int scvIdx) const
+// {
+// const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
+//
+// if (globalPos[1]<eps_)
+// return porosity_;
+// else
+// return 0.2;
+//
+// }
+
+
+ /*!
+ * \brief Define the initial porosity \f$[-]\f$ distribution
+ *
+ * \param element The finite element
+ * \param scv The sub-control volume
+ */
+ Scalar initialPorosity(const Element& element, const SubControlVolume &scv) const
+ {
+ Scalar phi;
+
+ if(isCharge_==true)
+ phi = 0.8; //direct charging acc. to Nagel et al 2014
+// phi = 0.887; //indirect charging acc. to Schmitt 2016
+ else
+ phi = 0.604; //direct charging acc. to Nagel et al 2014
+// phi = 0.772; //indirect charging acc. to Schmitt 2016
+
+ return phi;
+ }
+
+
+ /*!
+ * \brief Define the minimum porosity \f$[-]\f$ distribution
+ *
+ * \param element The finite element
+ * \param scv The sub-control volume
+ */
+ Scalar minPorosity(const Element& element, const SubControlVolume &scv) const
+ { return 0.604; //intrinsic porosity of CaO2H2 (see Nagel et al. 2014)
+ // return 0.772; //indirect charging acc. to Schmitt 2016
+ }
+
+ /*!
+ * \brief Define the minimum porosity \f$[-]\f$ after clogging caused by mineralization
+ *
+ * \param element The finite element
+ * \param scv The sub-control volume
+ */
+ Scalar porosity(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return poroLaw_.evaluatePorosity(element, scv, elemSol); }
+
+
+
+ Scalar solidity(const SubControlVolume &scv) const
+ { return 1.0 - porosityAtPos(scv.center()); }
+
+ /*!
+ * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The element
+ * \param scv The sub control volume
+ * \param elemSol The element solution vector
+ */
+ Scalar solidHeatCapacity(const Element &element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ {
+ return 790;
+// 42 // specific heat capacity of CaO [J / (kg K)]
+// * 3370 // density of CaO [kg/m^3]
+// * (1 - porosity(element, fvGeometry, scvIdx)); // for CaO only!!
+ }
+
+ /*!
+ * \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The element
+ * \param scv The sub control volume
+ * \param elemSol The element solution vector
+ */
+ Scalar solidDensity(const Element &element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ {
+// return 3370; // density of CaO [kg/m^3]
+ return 2600;
+ }
+
+ /*!
+ * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
+ *
+ * \param element The element
+ * \param scv The sub control volume
+ * \param elemSol The element solution vector
+ */
+ Scalar solidThermalConductivity(const Element &element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return lambdaSolid_; }
+
+private:
+// DimMatrix K_;
+// Scalar porosity_;
+ Scalar eps_;
+// MaterialLawParams materialParams_;
+ Scalar lambdaSolid_;
+ bool isCharge_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Problem, isCharge);
+ PorosityLaw poroLaw_;
+ PermeabilityLaw permLaw_;
+};
+
+}//end namespace
+
+#endif
diff --git a/test/porousmediumflow/CMakeLists.txt b/test/porousmediumflow/CMakeLists.txt
index 9d36688958b1627f7d5a280198831ab5a8a63a39..3072fc0b2800f4a0ae6c58e73a8b18908e1b65fe 100644
--- a/test/porousmediumflow/CMakeLists.txt
+++ b/test/porousmediumflow/CMakeLists.txt
@@ -1,6 +1,7 @@
add_subdirectory("1p")
add_subdirectory("1p2c")
add_subdirectory("1pnc")
+add_subdirectory("1pncmin")
add_subdirectory("2p")
add_subdirectory("2p1c")
add_subdirectory("2p2c")