| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup OnePNCMinTests | ||
| 10 | * \brief Definition of a problem for thermochemical heat storage using \f$ \textnormal{CaO}, | ||
| 11 | * \textnormal{Ca} \left( \textnormal{OH} \right)_2\f$. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_THERMOCHEM_PROBLEM_HH | ||
| 14 | #define DUMUX_THERMOCHEM_PROBLEM_HH | ||
| 15 | |||
| 16 | #include <dumux/common/properties.hh> | ||
| 17 | #include <dumux/common/parameters.hh> | ||
| 18 | #include <dumux/common/boundarytypes.hh> | ||
| 19 | #include <dumux/common/numeqvector.hh> | ||
| 20 | |||
| 21 | #include <dumux/discretization/elementsolution.hh> | ||
| 22 | #include <dumux/porousmediumflow/problem.hh> | ||
| 23 | |||
| 24 | #include "reaction.hh" | ||
| 25 | |||
| 26 | namespace Dumux { | ||
| 27 | |||
| 28 | /*! | ||
| 29 | * \ingroup OnePNCMinTests | ||
| 30 | * | ||
| 31 | * \brief Test for the 1pncmin model in combination with the NI model for a quasi batch | ||
| 32 | * reaction of Calciumoxyde to Calciumhydroxide. | ||
| 33 | * | ||
| 34 | * The boundary conditions of the batch test are such, that there are no gradients | ||
| 35 | * for temperature, pressure and gas water concentration within the reactor. | ||
| 36 | * The test only runs for the box discretization. | ||
| 37 | */ | ||
| 38 | template <class TypeTag> | ||
| 39 | class ThermoChemProblem : public PorousMediumFlowProblem<TypeTag> | ||
| 40 | { | ||
| 41 | using ParentType = PorousMediumFlowProblem<TypeTag>; | ||
| 42 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 43 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 44 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 45 | using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>; | ||
| 46 | |||
| 47 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 48 | using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView; | ||
| 49 | using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView; | ||
| 50 | using VolumeVariables = typename GridVariables::GridVolumeVariables::VolumeVariables; | ||
| 51 | |||
| 52 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 53 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 54 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
| 55 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 56 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 57 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 58 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 59 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 60 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
| 61 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | ||
| 62 | using ReactionRate = ThermoChemReaction; | ||
| 63 | |||
| 64 | enum { dim = GridView::dimension }; | ||
| 65 | enum { dimWorld = GridView::dimensionworld }; | ||
| 66 | |||
| 67 | enum | ||
| 68 | { | ||
| 69 | // Indices of the primary variables | ||
| 70 | pressureIdx = Indices::pressureIdx, //gas-phase pressure | ||
| 71 | H2OIdx = FluidSystem::compIdx(FluidSystem::MultiPhaseFluidSystem::H2OIdx), // mole fraction water | ||
| 72 | |||
| 73 | CaOIdx = FluidSystem::numComponents, | ||
| 74 | CaO2H2Idx = FluidSystem::numComponents+1, | ||
| 75 | |||
| 76 | // Equation Indices | ||
| 77 | conti0EqIdx = Indices::conti0EqIdx, | ||
| 78 | |||
| 79 | // Phase Indices | ||
| 80 | cPhaseIdx = SolidSystem::comp0Idx, | ||
| 81 | |||
| 82 | temperatureIdx = Indices::temperatureIdx, | ||
| 83 | energyEqIdx = Indices::energyEqIdx | ||
| 84 | }; | ||
| 85 | |||
| 86 | using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; | ||
| 87 | |||
| 88 | public: | ||
| 89 | 1 | ThermoChemProblem(std::shared_ptr<const GridGeometry> gridGeometry) | |
| 90 |
10/32✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 1 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
|
3 | : ParentType(gridGeometry) |
| 91 | { | ||
| 92 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
|
1 | name_ = getParam<std::string>("Problem.Name"); |
| 93 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | FluidSystem::init(/*tempMin=*/473.15, |
| 94 | /*tempMax=*/623.0, | ||
| 95 | /*numTemptempSteps=*/25, | ||
| 96 | /*startPressure=*/0, | ||
| 97 | /*endPressure=*/9e6, | ||
| 98 | /*pressureSteps=*/200); | ||
| 99 | |||
| 100 | // obtain BCs | ||
| 101 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | boundaryPressure_ = getParam<Scalar>("Problem.BoundaryPressure"); |
| 102 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | boundaryVaporMoleFrac_ = getParam<Scalar>("Problem.BoundaryMoleFraction"); |
| 103 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | boundaryTemperature_ = getParam<Scalar>("Problem.BoundaryTemperature"); |
| 104 | |||
| 105 | 1 | unsigned int codim = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethods::box ? dim : 0; | |
| 106 |
4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
|
3 | permeability_.resize(gridGeometry->gridView().size(codim)); |
| 107 |
4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
|
3 | porosity_.resize(gridGeometry->gridView().size(codim)); |
| 108 |
4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
|
3 | reactionRate_.resize(gridGeometry->gridView().size(codim)); |
| 109 | 1 | } | |
| 110 | |||
| 111 | /*! | ||
| 112 | * \brief The problem name. | ||
| 113 | * | ||
| 114 | * This is used as a prefix for files generated by the simulation. | ||
| 115 | */ | ||
| 116 | const std::string name() const | ||
| 117 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | { return name_; } |
| 118 | |||
| 119 | /*! | ||
| 120 | * \brief Sets the currently used time step size. | ||
| 121 | * | ||
| 122 | * This is necessary to limit the source terms to the maximum possible rate. | ||
| 123 | */ | ||
| 124 | ✗ | void setTimeStepSize( Scalar timeStepSize ) | |
| 125 | { | ||
| 126 |
1/2✓ Branch 1 taken 151 times.
✗ Branch 2 not taken.
|
151 | timeStepSize_ = timeStepSize; |
| 127 | ✗ | } | |
| 128 | |||
| 129 | /*! | ||
| 130 | * \name Boundary conditions | ||
| 131 | * | ||
| 132 | * \brief Specifies which kind of boundary condition should be | ||
| 133 | * used for which equation on a given boundary segment | ||
| 134 | * | ||
| 135 | * \param globalPos The global position | ||
| 136 | */ | ||
| 137 | ✗ | BoundaryTypes boundaryTypesAtPos( const GlobalPosition &globalPos) const | |
| 138 | { | ||
| 139 | ✗ | BoundaryTypes values; | |
| 140 | |||
| 141 | // we don't set any BCs for the solid phases | ||
| 142 | ✗ | values.setDirichlet(pressureIdx); | |
| 143 | ✗ | values.setDirichlet(H2OIdx); | |
| 144 | ✗ | values.setDirichlet(temperatureIdx); | |
| 145 | |||
| 146 | ✗ | return values; | |
| 147 | } | ||
| 148 | |||
| 149 | /*! | ||
| 150 | * \brief Evaluates the boundary conditions for a Dirichlet boundary segment | ||
| 151 | * | ||
| 152 | * \param globalPos The global position | ||
| 153 | */ | ||
| 154 | ✗ | PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const | |
| 155 | { | ||
| 156 | 3200 | PrimaryVariables priVars(0.0); | |
| 157 | |||
| 158 | 3200 | priVars[pressureIdx] = boundaryPressure_; | |
| 159 | 3200 | priVars[H2OIdx] = boundaryVaporMoleFrac_; | |
| 160 | 3200 | priVars[temperatureIdx] = boundaryTemperature_; | |
| 161 | 3200 | priVars[CaO2H2Idx] = 0.0; | |
| 162 | 6400 | priVars[CaOIdx] = 0.2; | |
| 163 | |||
| 164 | ✗ | return priVars; | |
| 165 | } | ||
| 166 | |||
| 167 | /*! | ||
| 168 | * \brief Evaluates the boundary conditions for a Neumann | ||
| 169 | * boundary segment in dependency on the current solution. | ||
| 170 | * | ||
| 171 | * \param element The element | ||
| 172 | * \param fvGeometry The finite volume geometry | ||
| 173 | * \param elemVolVars The element volume variables | ||
| 174 | * \param elemFluxVarsCache Flux variables caches for all faces in stencil | ||
| 175 | * \param scvf The subcontrolvolume face | ||
| 176 | * | ||
| 177 | * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ | ||
| 178 | * Negative values indicate an inflow. | ||
| 179 | */ | ||
| 180 | |||
| 181 | ✗ | NumEqVector neumann(const Element& element, | |
| 182 | const FVElementGeometry& fvGeometry, | ||
| 183 | const ElementVolumeVariables& elemVolVars, | ||
| 184 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 185 | const SubControlVolumeFace& scvf) const | ||
| 186 | { | ||
| 187 | ✗ | NumEqVector flux(0.0); | |
| 188 | ✗ | return flux; | |
| 189 | } | ||
| 190 | |||
| 191 | /*! | ||
| 192 | * \brief Evaluates the initial values for a control volume in | ||
| 193 | * \f$ [ \textnormal{unit of primary variables} ] \f$ | ||
| 194 | * | ||
| 195 | * \param globalPos The global position | ||
| 196 | */ | ||
| 197 | ✗ | PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const | |
| 198 | { | ||
| 199 | ✗ | PrimaryVariables priVars(0.0); | |
| 200 | |||
| 201 | Scalar pInit; | ||
| 202 | Scalar tInit; | ||
| 203 | Scalar h2oInit; | ||
| 204 | Scalar CaOInit; | ||
| 205 | Scalar CaO2H2Init; | ||
| 206 | |||
| 207 | ✗ | pInit = getParam<Scalar>("Problem.PressureInitial"); | |
| 208 | ✗ | tInit = getParam<Scalar>("Problem.TemperatureInitial"); | |
| 209 | ✗ | h2oInit = getParam<Scalar>("Problem.VaporInitial"); | |
| 210 | ✗ | CaOInit = getParam<Scalar>("Problem.CaOInitial"); | |
| 211 | ✗ | CaO2H2Init = getParam<Scalar>("Problem.CaO2H2Initial"); | |
| 212 | |||
| 213 | ✗ | priVars[pressureIdx] = pInit; | |
| 214 | ✗ | priVars[H2OIdx] = h2oInit; | |
| 215 | ✗ | priVars[temperatureIdx] = tInit; | |
| 216 | |||
| 217 | // these values are not used, as we didn't set BCs | ||
| 218 | // for the solid phases. For cell-centered models it is | ||
| 219 | // important to set the values to fully define Dirichlet BCs | ||
| 220 | ✗ | priVars[CaOIdx] = CaOInit; | |
| 221 | ✗ | priVars[CaO2H2Idx] = CaO2H2Init; | |
| 222 | |||
| 223 | ✗ | return priVars; | |
| 224 | } | ||
| 225 | |||
| 226 | /*! | ||
| 227 | * \brief Evaluates the source term for all phases within a given | ||
| 228 | * sub-control volume in units of \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$. | ||
| 229 | * | ||
| 230 | * This is the method for the case where the source term is | ||
| 231 | * potentially solution dependent and requires some quantities that | ||
| 232 | * are specific to the fully-implicit method. | ||
| 233 | * | ||
| 234 | * \param element The finite element | ||
| 235 | * \param fvGeometry The finite-volume geometry | ||
| 236 | * \param elemVolVars All volume variables for the element | ||
| 237 | * \param scv The subcontrolvolume | ||
| 238 | * | ||
| 239 | * For this method, the \a values parameter stores the conserved quantity rate | ||
| 240 | * generated or annihilated per volume unit. Positive values mean | ||
| 241 | * that the conserved quantity is created, negative ones mean that it vanishes. | ||
| 242 | * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. | ||
| 243 | */ | ||
| 244 | 67200 | NumEqVector source(const Element &element, | |
| 245 | const FVElementGeometry& fvGeometry, | ||
| 246 | const ElementVolumeVariables& elemVolVars, | ||
| 247 | const SubControlVolume &scv) const | ||
| 248 | { | ||
| 249 | |||
| 250 | 67200 | NumEqVector source(0.0); | |
| 251 | 67200 | const auto& volVars = elemVolVars[scv]; | |
| 252 | |||
| 253 | 67200 | Scalar qMass = rrate_.thermoChemReaction(volVars); | |
| 254 | 134400 | Scalar qMole = qMass/FluidSystem::molarMass(H2OIdx)*(1-volVars.porosity()); | |
| 255 | |||
| 256 | // make sure not more solid reacts than present | ||
| 257 | // In this test, we only consider discharge. Therefore, we use the cPhaseIdx for CaO. | ||
| 258 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 67200 times.
|
134400 | if (-qMole*timeStepSize_ + volVars.solidVolumeFraction(cPhaseIdx)* volVars.solidComponentMolarDensity(cPhaseIdx) < 0 + eps_) |
| 259 | { | ||
| 260 | ✗ | qMole = -volVars.solidVolumeFraction(cPhaseIdx)* volVars.solidComponentMolarDensity(cPhaseIdx)/timeStepSize_; | |
| 261 | } | ||
| 262 | 134400 | source[conti0EqIdx+CaO2H2Idx] = qMole; | |
| 263 | 134400 | source[conti0EqIdx+CaOIdx] = - qMole; | |
| 264 | 134400 | source[conti0EqIdx+H2OIdx] = - qMole; | |
| 265 | |||
| 266 | 67200 | Scalar deltaH = 108e3; // J/mol | |
| 267 | 134400 | source[energyEqIdx] = qMole * deltaH; | |
| 268 | |||
| 269 | 67200 | return source; | |
| 270 | } | ||
| 271 | |||
| 272 | |||
| 273 | /*! | ||
| 274 | * \brief Returns the permeability. | ||
| 275 | */ | ||
| 276 | const std::vector<Scalar>& getPerm() | ||
| 277 | { | ||
| 278 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return permeability_; |
| 279 | } | ||
| 280 | |||
| 281 | /*! | ||
| 282 | * \brief Returns the porosity. | ||
| 283 | */ | ||
| 284 | const std::vector<Scalar>& getPoro() | ||
| 285 | { | ||
| 286 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return porosity_; |
| 287 | } | ||
| 288 | |||
| 289 | /*! | ||
| 290 | * \brief Returns the reaction rate. | ||
| 291 | */ | ||
| 292 | const std::vector<Scalar>& getRRate() | ||
| 293 | { | ||
| 294 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return reactionRate_; |
| 295 | } | ||
| 296 | |||
| 297 | /*! | ||
| 298 | * \brief Adds additional VTK output data to the VTKWriter. | ||
| 299 | * | ||
| 300 | * Function is called by the output module on every write. | ||
| 301 | */ | ||
| 302 | 152 | void updateVtkOutput(const SolutionVector& curSol) | |
| 303 | { | ||
| 304 |
2/4✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 152 times.
✗ Branch 5 not taken.
|
304 | auto fvGeometry = localView(this->gridGeometry()); |
| 305 |
5/10✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 152 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 152 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 456 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 304 times.
✗ Branch 14 not taken.
|
1216 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 306 | { | ||
| 307 |
2/4✓ Branch 1 taken 304 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 304 times.
✗ Branch 5 not taken.
|
608 | const auto elemSol = elementSolution(element, curSol, this->gridGeometry()); |
| 308 |
1/2✓ Branch 1 taken 304 times.
✗ Branch 2 not taken.
|
304 | fvGeometry.bindElement(element); |
| 309 | |||
| 310 |
4/4✓ Branch 0 taken 1216 times.
✓ Branch 1 taken 304 times.
✓ Branch 2 taken 1216 times.
✓ Branch 3 taken 304 times.
|
3040 | for (auto&& scv : scvs(fvGeometry)) |
| 311 | { | ||
| 312 |
1/2✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
|
1216 | VolumeVariables volVars; |
| 313 |
1/2✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
|
1216 | volVars.update(elemSol, *this, element, scv); |
| 314 | 1216 | const auto dofIdxGlobal = scv.dofIndex(); | |
| 315 | 2432 | permeability_[dofIdxGlobal] = this->spatialParams().permeability(element, scv, elemSol); | |
| 316 |
1/2✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
|
2432 | porosity_[dofIdxGlobal] = volVars.porosity(); |
| 317 |
1/2✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
|
1216 | reactionRate_[dofIdxGlobal] = rrate_.thermoChemReaction(volVars); |
| 318 | } | ||
| 319 | } | ||
| 320 | 152 | } | |
| 321 | |||
| 322 | private: | ||
| 323 | std::string name_; | ||
| 324 | |||
| 325 | static constexpr Scalar eps_ = 1e-6; | ||
| 326 | |||
| 327 | // boundary conditions | ||
| 328 | Scalar boundaryPressure_; | ||
| 329 | Scalar boundaryVaporMoleFrac_; | ||
| 330 | Scalar boundaryTemperature_; | ||
| 331 | |||
| 332 | std::vector<double> permeability_; | ||
| 333 | std::vector<double> porosity_; | ||
| 334 | std::vector<double> reactionRate_; | ||
| 335 | |||
| 336 | ReactionRate rrate_; | ||
| 337 | Scalar timeStepSize_; | ||
| 338 | }; | ||
| 339 | |||
| 340 | } // end namespace Dumux | ||
| 341 | |||
| 342 | #endif | ||
| 343 |