From e1c26bf5fa45ee429e5a082a71d15934c8899f9e Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Sun, 21 Feb 2016 23:53:51 +0100
Subject: [PATCH] [electrochem] Generalize electrochemistry
* Make current density publicly available
* Use standard grid creator
* Implement reaction source term for box and cell-centered
---
.../electrochemistry/electrochemistry.hh | 57 +++++++++++--------
.../electrochemistry/electrochemistryni.hh | 55 ++++++++++--------
.../2pnc/implicit/fuelcellproblem.hh | 7 +--
.../2pnc/implicit/test_2pnc.input | 56 +++++++-----------
.../2pnc/implicit/test_box2pnc.cc | 20 +------
5 files changed, 86 insertions(+), 109 deletions(-)
diff --git a/dumux/material/chemistry/electrochemistry/electrochemistry.hh b/dumux/material/chemistry/electrochemistry/electrochemistry.hh
index 86b91d100d..c67e8b5559 100644
--- a/dumux/material/chemistry/electrochemistry/electrochemistry.hh
+++ b/dumux/material/chemistry/electrochemistry/electrochemistry.hh
@@ -26,6 +26,7 @@
#include <cmath>
+#include <dune/common/deprecated.hh>
#include <dumux/common/basicproperties.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/material/constants.hh>
@@ -97,6 +98,12 @@ class ElectroChemistry
energyEqIdx = FluidSystem::numComponents //energy equation
};
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? GridView::dimension : 0 };
+
+ using GlobalPosition = typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
+ using CellVector = typename Dune::FieldVector<Scalar, GridView::dimension>;
+
public:
/*!
* \brief Calculates reaction sources with an electrochemical model approach.
@@ -107,31 +114,22 @@ public:
* For this method, the \a values parameter stores source values
*/
static void reactionSource(PrimaryVariables &values,
- const VolumeVariables &volVars)
+ Scalar currentDensity)
{
- static Scalar transportNumberH2O = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.TransportNumberH20);
- static Scalar maxIter = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.MaxIterations);
- static Scalar gridYMin = 0.0;
- static Scalar gridYMax = GET_RUNTIME_PARAM(TypeTag, Scalar, Grid.UpperRightY);
- static Scalar nCellsY = GET_RUNTIME_PARAM(TypeTag, Scalar, Grid.NumberOfCellsY);
-
- //initialise current density
- Scalar currentDensity = 0.0;
-
- //call internal method to calculate the current density
- currentDensity = calculateCurrentDensity_(volVars, maxIter);
-
//correction to account for actually relevant reaction area
//current density has to be devided by the half length of the box
- Scalar lengthBox = (gridYMax - gridYMin)/nCellsY;
-
- if(electroChemistryModel == ElectroChemistryModel::Acosta)
- currentDensity = currentDensity/lengthBox;
+ //\todo Do we have multiply with the electrochemically active surface area (ECSA) here instead?
+ static Scalar gridYMax = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, Grid.UpperRight)[1];
+ static Scalar nCellsY = GET_RUNTIME_PARAM(TypeTag, CellVector, Grid.Cells)[1];
+
+ // Warning: This assumes the reaction layer is always just one cell (cell-centered) or half a box (box) thick
+ const auto lengthBox = gridYMax/nCellsY;
+ if (isBox)
+ currentDensity *= 2.0/lengthBox;
else
- currentDensity = currentDensity*2/lengthBox;
+ currentDensity *= 1.0/lengthBox;
- //conversion from [A/cm^2] to [A/m^2]
- currentDensity = currentDensity*10000;
+ static Scalar transportNumberH2O = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.TransportNumberH20);
//calculation of flux terms with faraday equation
values[contiH2OEqIdx] = currentDensity/(2*Constant::F); //reaction term in reaction layer
@@ -139,13 +137,20 @@ public:
values[contiO2EqIdx] = -currentDensity/(4*Constant::F); //O2-equation
}
-protected:
+ DUNE_DEPRECATED_MSG("First compute the currentDensity with calculateCurrentDensity(const VolumeVariables&) and then use the method reactionSource(PrimaryVariables&, Scalar) instead")
+ static void reactionSource(PrimaryVariables &values,
+ const VolumeVariables &volVars)
+ {
+ reactionSource(values, calculateCurrentDensity(volVars));
+ }
/*!
* \brief Newton solver for calculation of the current density.
+ * \returns The current density in A/m^2
*/
- static Scalar calculateCurrentDensity_(const VolumeVariables &volVars, Scalar maxIter)
+ static Scalar calculateCurrentDensity(const VolumeVariables &volVars)
{
+ static Scalar maxIter = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.MaxIterations);
static Scalar specificResistance = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.SpecificResistance);
static Scalar reversibleVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.ReversibleVoltage);
static Scalar cellVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.CellVoltage);
@@ -201,9 +206,15 @@ protected:
}
}
- return currentDensity;
+ //conversion from [A/cm^2] to [A/m^2]
+ return currentDensity*10000;
}
+protected:
+ DUNE_DEPRECATED_MSG("This is now a public member function (the name lost the underscore postfix.)")
+ static Scalar calculateCurrentDensity_(const VolumeVariables &volVars)
+ { calculateCurrentDensity(volVars); }
+
private:
/*!
diff --git a/dumux/material/chemistry/electrochemistry/electrochemistryni.hh b/dumux/material/chemistry/electrochemistry/electrochemistryni.hh
index b55df00a82..6db0f53fa9 100644
--- a/dumux/material/chemistry/electrochemistry/electrochemistryni.hh
+++ b/dumux/material/chemistry/electrochemistry/electrochemistryni.hh
@@ -52,6 +52,7 @@ class ElectroChemistryNI : public ElectroChemistry<TypeTag, electroChemistryMode
typedef ElectroChemistry<TypeTag, electroChemistryModel> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
@@ -72,6 +73,12 @@ class ElectroChemistryNI : public ElectroChemistry<TypeTag, electroChemistryMode
energyEqIdx = FluidSystem::numComponents, //energy equation
};
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? GridView::dimension : 0 };
+
+ using GlobalPosition = typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
+ using CellVector = typename Dune::FieldVector<Scalar, GridView::dimension>;
+
public:
/*!
* \brief Calculates reaction sources with an electrochemical model approach.
@@ -82,42 +89,40 @@ public:
* For this method, the \a values parameter stores source values
*/
static void reactionSource(PrimaryVariables &values,
- const VolumeVariables &volVars)
+ Scalar currentDensity)
{
- static Scalar transportNumberH2O = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.TransportNumberH20);
- static Scalar maxIter = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.MaxIterations);
- static Scalar gridYMin = 0.0;
- static Scalar gridYMax = GET_RUNTIME_PARAM(TypeTag, Scalar, Grid.UpperRightY);
- static Scalar nCellsY = GET_RUNTIME_PARAM(TypeTag, Scalar, Grid.NumberOfCellsY);
- static Scalar thermoneutralVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.ThermoneutralVoltage);
- static Scalar cellVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.CellVoltage);
-
- //initialise current density
- Scalar currentDensity = 0.0;
-
- //call internal method to calculate the current density
- currentDensity = ParentType::calculateCurrentDensity_(volVars, maxIter);
-
//correction to account for actually relevant reaction area
//current density has to be devided by the half length of the box
- Scalar lengthBox = (gridYMax - gridYMin)/nCellsY;
-
- if(electroChemistryModel == ElectroChemistryModel::Acosta)
- currentDensity = currentDensity/lengthBox;
+ //\todo Do we have multiply with the electrochemically active surface area (ECSA) here instead?
+ static Scalar gridYMax = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, Grid.UpperRight)[1];
+ static Scalar nCellsY = GET_RUNTIME_PARAM(TypeTag, CellVector, Grid.Cells)[1];
+
+ // Warning: This assumes the reaction layer is always just one cell (cell-centered) or half a box (box) thick
+ const auto lengthBox = gridYMax/nCellsY;
+ if (isBox)
+ currentDensity *= 2.0/lengthBox;
else
- currentDensity = currentDensity*2/lengthBox;
+ currentDensity *= 1.0/lengthBox;
- //conversion from [A/cm^2] to [A/m^2]
- currentDensity = currentDensity*10000;
+ static Scalar transportNumberH2O = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.TransportNumberH20);
+ static Scalar thermoneutralVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.ThermoneutralVoltage);
+ static Scalar cellVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.CellVoltage);
//calculation of flux terms with faraday equation
values[contiH2OEqIdx] = currentDensity/(2*Constant::F); //reaction term in reaction layer
- values[contiH2OEqIdx] += currentDensity/Constant::F*transportNumberH2O; //osmotic term in membrane
+ values[contiH2OEqIdx] += currentDensity/Constant::F*transportNumberH2O; //osmotic term in membrane
values[contiO2EqIdx] = -currentDensity/(4*Constant::F); //O2-equation
values[energyEqIdx] = (thermoneutralVoltage - cellVoltage)*currentDensity; //energy equation
}
+
+ DUNE_DEPRECATED_MSG("First compute the currentDensity with calculateCurrentDensity(const VolumeVariables&) and then use the method reactionSource(PrimaryVariables&, Scalar) instead")
+ static void reactionSource(PrimaryVariables &values,
+ const VolumeVariables &volVars)
+ {
+ reactionSource(values, this->calculateCurrentDensity(volVars));
+ }
};
-}// end namespace
-#endif
+}// end namespace
+#endif
diff --git a/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh b/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
index a06a097eb0..1f8501e8ed 100644
--- a/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
+++ b/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
@@ -24,7 +24,6 @@
#ifndef DUMUX_FUELCELL_PROBLEM_HH
#define DUMUX_FUELCELL_PROBLEM_HH
-#include <dumux/io/cubegridcreator.hh>
#include <dumux/porousmediumflow/2pnc/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/material/fluidsystems/h2on2o2.hh>
@@ -46,11 +45,7 @@ NEW_TYPE_TAG(FuelCellBoxProblem, INHERITS_FROM(BoxModel, FuelCellProblem));
NEW_TYPE_TAG(FuelCellCCProblem, INHERITS_FROM(CCModel, FuelCellProblem));
// Set the grid type
-#if HAVE_UG
-SET_TYPE_PROP(FuelCellProblem, Grid, Dune::UGGrid<2>);
-#endif
-// Set the grid creator
-SET_TYPE_PROP(FuelCellProblem, GridCreator, Dumux::CubeGridCreator<TypeTag>);
+SET_TYPE_PROP(FuelCellProblem, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(FuelCellProblem, Problem, Dumux::FuelCellProblem<TypeTag>);
// Set the primary variable combination for the 2pnc model
diff --git a/test/porousmediumflow/2pnc/implicit/test_2pnc.input b/test/porousmediumflow/2pnc/implicit/test_2pnc.input
index 2bc445b56c..387623cb1d 100644
--- a/test/porousmediumflow/2pnc/implicit/test_2pnc.input
+++ b/test/porousmediumflow/2pnc/implicit/test_2pnc.input
@@ -1,19 +1,15 @@
-# Parameter file for test case 2pncmin.
+# Parameter file for test case 2pnc.
# Everything behind a '#' is a comment.
-# Type "./test_2pncmin --help" for more information.
+# Type "./test_2pnc --help" for more information.
[TimeManager]
-DtInitial = 5e-1 # initial time step size [s]
-MaxTimeStepSize = 1000 # maximum time step size
-TEnd= 1e3 # duration of the simulation [s]
+DtInitial = 5e-1 # [s] initial time step size
+MaxTimeStepSize = 1000 # [s] maximum time step size
+TEnd= 1e3 # [s] duration of the simulation
[Grid]
-#Number of cells in X and Y directions
-NumberOfCellsX = 21
-NumberOfCellsY = 6
-#Extend of entire domain
-UpperRightX = 2.0e-3 # x-coordinate of the upper-right corner of the grid [m]
-UpperRightY = 5.5e-4 # y-coordinate of the upper-right corner of the grid [m]
+UpperRight = 2.0e-3 5.5e-4 # [m] upper right corner coordinates
+Cells = 21 6 # [-] number of cells in x,y-direction
[Problem]
Name = fuelcell
@@ -29,29 +25,17 @@ TemperatureHigh = 314.15 # [K] upper temperature limit for tabularization
InitialTemperature = 310 # [K] initial temperature for tabularization
[ElectroChemistry]
-SpecificResistance = 0.25 #[Ohm*cm^2]
-ReversibleVoltage = 1.191 #[V]
-CellVoltage = 0.5 #[V]
-ThermoneutralVoltage = 1.4836 #[V]
-RefCurrentDensity = 1.87e-8 #[A/cm] for calculating exchange current density at reference conditions
-RefO2PartialPressure = 5.0e5 #[Pa] for calculating exchange current density at reference conditions
-RefTemperature = 353.25 #[K] for calculating exchange current density at reference conditions
-ActivationBarrier = 73.0e3 #[J/mol]
-TransferCoefficient = 0.5 #[-] (alpha)
-NumElectrons = 2 #[-] number of electrons transferred in reaction
+SpecificResistance = 0.25 # [Ohm*cm^2]
+ReversibleVoltage = 1.191 # [V]
+CellVoltage = 0.5 # [V]
+ThermoneutralVoltage = 1.4836 # [V]
+RefCurrentDensity = 1.87e-8 # [A/cm^2] for calculating exchange current density at reference conditions
+RefO2PartialPressure = 5.0e5 # [Pa] for calculating exchange current density at reference conditions
+RefTemperature = 353.25 # [K] for calculating exchange current density at reference conditions
+ActivationBarrier = 73.0e3 # [J/mol]
+TransferCoefficient = 0.5 # [-] (alpha)
+NumElectrons = 2 # [-] number of electrons transferred in reaction
SurfaceIncreasingFactor = 60
-TransportNumberH20 = 0.2327 #[-] account for osmotic H20 transport term from membrane, value Lena Walter
-pO2Inlet = 0.3e5 #[Pa] partial pressure of oxygen at gas channel inlet
-MaxIterations = 300 #[-] Maximum number for iterations for solving electrochemical system
-
-[Newton]
-RelTolerance = 1e-8
-TargetSteps = 10
-MaxSteps = 15
-WriteConvergence = false
-MaxTimeStepDivisions = 20
-
-[LinearSolver]
-ResidualReduction = 1e-8
-Verbosity = 0
-MaxIterations = 100
\ No newline at end of file
+TransportNumberH20 = 0.2327 # [-] account for osmotic H20 transport term from membrane, value Lena Walter
+pO2Inlet = 0.3e5 # [Pa] partial pressure of oxygen at gas channel inlet
+MaxIterations = 300 # [-] Maximum number for iterations for solving electrochemical system
diff --git a/test/porousmediumflow/2pnc/implicit/test_box2pnc.cc b/test/porousmediumflow/2pnc/implicit/test_box2pnc.cc
index 6c99ad2656..826b60d121 100644
--- a/test/porousmediumflow/2pnc/implicit/test_box2pnc.cc
+++ b/test/porousmediumflow/2pnc/implicit/test_box2pnc.cc
@@ -41,19 +41,7 @@ void usage(const char *progName, const std::string &errorMsg)
errorMessageOut += " [options]\n";
errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
- "\t-TimeManager.TEnd End of the simulation [s] \n"
- "\t-TimeManager.DtInitial Initial timestep size [s] \n"
- "\t-Grid.File Name of the file containing the grid \n"
- "\t definition in DGF format\n"
- "\t-FluidSystem.NTemperature Number of tabularization entries [-] \n"
- "\t-FluidSystem.NPressure Number of tabularization entries [-] \n"
- "\t-FluidSystem.PressureLow Low end for tabularization of fluid properties [Pa] \n"
- "\t-FluidSystem.PressureHigh High end for tabularization of fluid properties [Pa] \n"
- "\t-FluidSystem.TemperatureLow Low end for tabularization of fluid properties [Pa] \n"
- "\t-FluidSystem.TemperatureHigh High end for tabularization of fluid properties [Pa] \n"
- "\t-SimulationControl.Name The name of the output files [-] \n"
- "\t-InitialConditions.Temperature Initial temperature in the reservoir [K] \n"
- "\t-InitialConditions.DepthBOR Depth below ground surface [m] \n";
+ "\t-ParameterFile Parameter file (Input file) \n";
std::cout << errorMessageOut
<< "\n";
@@ -62,12 +50,6 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv)
{
-#if HAVE_UG
typedef TTAG(FuelCellBoxProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
-#else
-#warning You need UGGrid to run this test.
- std::cerr << "You need UGGrid to run this test." << std::endl;
- return 77;
-#endif
}
--
GitLab