Skip to content
Snippets Groups Projects
Commit e1c26bf5 authored by Timo Koch's avatar Timo Koch
Browse files

[electrochem] Generalize electrochemistry 

* Make current density publicly available
* Use standard grid creator
* Implement reaction source term for box and cell-centered
parent 8fb19f03
No related branches found
No related tags found
1 merge request!80Feature/generalize 2pnc
Hide whitespace changes
Inline Side-by-side
dumux/material/chemistry/electrochemistry/electrochemistry.hh
+ 34
− 23
View file @ e1c26bf5
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
#include <cmath> #include <cmath>
#include <dune/common/deprecated.hh>
#include <dumux/common/basicproperties.hh> #include <dumux/common/basicproperties.hh>
#include <dumux/common/exceptions.hh> #include <dumux/common/exceptions.hh>
#include <dumux/material/constants.hh> #include <dumux/material/constants.hh>
...@@ -97,6 +98,12 @@ class ElectroChemistry ...@@ -97,6 +98,12 @@ class ElectroChemistry
energyEqIdx = FluidSystem::numComponents //energy equation 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: public:
/*! /*!
* \brief Calculates reaction sources with an electrochemical model approach. * \brief Calculates reaction sources with an electrochemical model approach.
...@@ -107,31 +114,22 @@ public: ...@@ -107,31 +114,22 @@ public:
* For this method, the \a values parameter stores source values * For this method, the \a values parameter stores source values
*/ */
static void reactionSource(PrimaryVariables &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 //correction to account for actually relevant reaction area
//current density has to be devided by the half length of the box //current density has to be devided by the half length of the box
Scalar lengthBox = (gridYMax - gridYMin)/nCellsY; //\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];
if(electroChemistryModel == ElectroChemistryModel::Acosta) static Scalar nCellsY = GET_RUNTIME_PARAM(TypeTag, CellVector, Grid.Cells)[1];
currentDensity = currentDensity/lengthBox;
// 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 else
currentDensity = currentDensity*2/lengthBox; currentDensity *= 1.0/lengthBox;
//conversion from [A/cm^2] to [A/m^2] static Scalar transportNumberH2O = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.TransportNumberH20);
currentDensity = currentDensity*10000;
//calculation of flux terms with faraday equation //calculation of flux terms with faraday equation
values[contiH2OEqIdx] = currentDensity/(2*Constant::F); //reaction term in reaction layer values[contiH2OEqIdx] = currentDensity/(2*Constant::F); //reaction term in reaction layer
...@@ -139,13 +137,20 @@ public: ...@@ -139,13 +137,20 @@ public:
values[contiO2EqIdx] = -currentDensity/(4*Constant::F); //O2-equation 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. * \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 specificResistance = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.SpecificResistance);
static Scalar reversibleVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.ReversibleVoltage); static Scalar reversibleVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.ReversibleVoltage);
static Scalar cellVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.CellVoltage); static Scalar cellVoltage = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.CellVoltage);
...@@ -201,9 +206,15 @@ protected: ...@@ -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: private:
/*! /*!
......
dumux/material/chemistry/electrochemistry/electrochemistryni.hh
+ 30
− 25
View file @ e1c26bf5
...@@ -52,6 +52,7 @@ class ElectroChemistryNI : public ElectroChemistry<TypeTag, electroChemistryMode ...@@ -52,6 +52,7 @@ class ElectroChemistryNI : public ElectroChemistry<TypeTag, electroChemistryMode
typedef ElectroChemistry<TypeTag, electroChemistryModel> ParentType; typedef ElectroChemistry<TypeTag, electroChemistryModel> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; 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, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
...@@ -72,6 +73,12 @@ class ElectroChemistryNI : public ElectroChemistry<TypeTag, electroChemistryMode ...@@ -72,6 +73,12 @@ class ElectroChemistryNI : public ElectroChemistry<TypeTag, electroChemistryMode
energyEqIdx = FluidSystem::numComponents, //energy equation 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: public:
/*! /*!
* \brief Calculates reaction sources with an electrochemical model approach. * \brief Calculates reaction sources with an electrochemical model approach.
...@@ -82,42 +89,40 @@ public: ...@@ -82,42 +89,40 @@ public:
* For this method, the \a values parameter stores source values * For this method, the \a values parameter stores source values
*/ */
static void reactionSource(PrimaryVariables &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 //correction to account for actually relevant reaction area
//current density has to be devided by the half length of the box //current density has to be devided by the half length of the box
Scalar lengthBox = (gridYMax - gridYMin)/nCellsY; //\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];
if(electroChemistryModel == ElectroChemistryModel::Acosta) static Scalar nCellsY = GET_RUNTIME_PARAM(TypeTag, CellVector, Grid.Cells)[1];
currentDensity = currentDensity/lengthBox;
// 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 else
currentDensity = currentDensity*2/lengthBox; currentDensity *= 1.0/lengthBox;
//conversion from [A/cm^2] to [A/m^2] static Scalar transportNumberH2O = GET_RUNTIME_PARAM(TypeTag, Scalar, ElectroChemistry.TransportNumberH20);
currentDensity = currentDensity*10000; 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 //calculation of flux terms with faraday equation
values[contiH2OEqIdx] = currentDensity/(2*Constant::F); //reaction term in reaction layer 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[contiO2EqIdx] = -currentDensity/(4*Constant::F); //O2-equation
values[energyEqIdx] = (thermoneutralVoltage - cellVoltage)*currentDensity; //energy 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
test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
+ 1
− 6
View file @ e1c26bf5
...@@ -24,7 +24,6 @@ ...@@ -24,7 +24,6 @@
#ifndef DUMUX_FUELCELL_PROBLEM_HH #ifndef DUMUX_FUELCELL_PROBLEM_HH
#define DUMUX_FUELCELL_PROBLEM_HH #define DUMUX_FUELCELL_PROBLEM_HH
#include <dumux/io/cubegridcreator.hh>
#include <dumux/porousmediumflow/2pnc/implicit/model.hh> #include <dumux/porousmediumflow/2pnc/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh> #include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/material/fluidsystems/h2on2o2.hh> #include <dumux/material/fluidsystems/h2on2o2.hh>
...@@ -46,11 +45,7 @@ NEW_TYPE_TAG(FuelCellBoxProblem, INHERITS_FROM(BoxModel, FuelCellProblem)); ...@@ -46,11 +45,7 @@ NEW_TYPE_TAG(FuelCellBoxProblem, INHERITS_FROM(BoxModel, FuelCellProblem));
NEW_TYPE_TAG(FuelCellCCProblem, INHERITS_FROM(CCModel, FuelCellProblem)); NEW_TYPE_TAG(FuelCellCCProblem, INHERITS_FROM(CCModel, FuelCellProblem));
// Set the grid type // Set the grid type
#if HAVE_UG SET_TYPE_PROP(FuelCellProblem, Grid, Dune::YaspGrid<2>);
SET_TYPE_PROP(FuelCellProblem, Grid, Dune::UGGrid<2>);
#endif
// Set the grid creator
SET_TYPE_PROP(FuelCellProblem, GridCreator, Dumux::CubeGridCreator<TypeTag>);
// Set the problem property // Set the problem property
SET_TYPE_PROP(FuelCellProblem, Problem, Dumux::FuelCellProblem<TypeTag>); SET_TYPE_PROP(FuelCellProblem, Problem, Dumux::FuelCellProblem<TypeTag>);
// Set the primary variable combination for the 2pnc model // Set the primary variable combination for the 2pnc model
......
test/porousmediumflow/2pnc/implicit/test_2pnc.input
+ 20
− 36
View file @ e1c26bf5
# Parameter file for test case 2pncmin. # Parameter file for test case 2pnc.
# Everything behind a '#' is a comment. # Everything behind a '#' is a comment.
# Type "./test_2pncmin --help" for more information. # Type "./test_2pnc --help" for more information.
[TimeManager] [TimeManager]
DtInitial = 5e-1 # initial time step size [s] DtInitial = 5e-1 # [s] initial time step size
MaxTimeStepSize = 1000 # maximum time step size MaxTimeStepSize = 1000 # [s] maximum time step size
TEnd= 1e3 # duration of the simulation [s] TEnd= 1e3 # [s] duration of the simulation
[Grid] [Grid]
#Number of cells in X and Y directions UpperRight = 2.0e-3 5.5e-4 # [m] upper right corner coordinates
NumberOfCellsX = 21 Cells = 21 6 # [-] number of cells in x,y-direction
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]
[Problem] [Problem]
Name = fuelcell Name = fuelcell
...@@ -29,29 +25,17 @@ TemperatureHigh = 314.15 # [K] upper temperature limit for tabularization ...@@ -29,29 +25,17 @@ TemperatureHigh = 314.15 # [K] upper temperature limit for tabularization
InitialTemperature = 310 # [K] initial temperature for tabularization InitialTemperature = 310 # [K] initial temperature for tabularization
[ElectroChemistry] [ElectroChemistry]
SpecificResistance = 0.25 #[Ohm*cm^2] SpecificResistance = 0.25 # [Ohm*cm^2]
ReversibleVoltage = 1.191 #[V] ReversibleVoltage = 1.191 # [V]
CellVoltage = 0.5 #[V] CellVoltage = 0.5 # [V]
ThermoneutralVoltage = 1.4836 #[V] ThermoneutralVoltage = 1.4836 # [V]
RefCurrentDensity = 1.87e-8 #[A/cm] for calculating exchange current density at reference conditions 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 RefO2PartialPressure = 5.0e5 # [Pa] for calculating exchange current density at reference conditions
RefTemperature = 353.25 #[K] 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] ActivationBarrier = 73.0e3 # [J/mol]
TransferCoefficient = 0.5 #[-] (alpha) TransferCoefficient = 0.5 # [-] (alpha)
NumElectrons = 2 #[-] number of electrons transferred in reaction NumElectrons = 2 # [-] number of electrons transferred in reaction
SurfaceIncreasingFactor = 60 SurfaceIncreasingFactor = 60
TransportNumberH20 = 0.2327 #[-] account for osmotic H20 transport term from membrane, value Lena Walter 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 pO2Inlet = 0.3e5 # [Pa] partial pressure of oxygen at gas channel inlet
MaxIterations = 300 #[-] Maximum number for iterations for solving electrochemical system 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
test/porousmediumflow/2pnc/implicit/test_box2pnc.cc
+ 1
− 19
View file @ e1c26bf5
...@@ -41,19 +41,7 @@ void usage(const char *progName, const std::string &errorMsg) ...@@ -41,19 +41,7 @@ void usage(const char *progName, const std::string &errorMsg)
errorMessageOut += " [options]\n"; errorMessageOut += " [options]\n";
errorMessageOut += errorMsg; errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
"\t-TimeManager.TEnd End of the simulation [s] \n" "\t-ParameterFile Parameter file (Input file) \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";
std::cout << errorMessageOut std::cout << errorMessageOut
<< "\n"; << "\n";
...@@ -62,12 +50,6 @@ void usage(const char *progName, const std::string &errorMsg) ...@@ -62,12 +50,6 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv) int main(int argc, char** argv)
{ {
#if HAVE_UG
typedef TTAG(FuelCellBoxProblem) ProblemTypeTag; typedef TTAG(FuelCellBoxProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage); 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
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment