Commit 271ea80f authored by Simon Scholz's avatar Simon Scholz Committed by Timo Koch
Browse files

Biomin tutorial

parent 6e423445
......@@ -4,4 +4,6 @@ add_custom_target(test_exercises)
add_subdirectory(exercise-basic)
add_subdirectory(exercise-properties)
add_subdirectory(exercise-fluidsystem)
add_subdirectory(exercise-biomineralization)
add_subdirectory(solution)
# executables for exercisebiomin
dune_add_test(NAME exercisebiomin
SOURCES exercisebiomin.cc
CMD_ARGS exercisebiomin.input)
# add tutorial to the common target
add_dependencies(test_exercises exercisebiomin)
# add symlinks for the input files
add_input_file_links(FILES "exercisebiomin.input")
# Exercise #4 (DuMuX course)
The aim of this exercise is to a first glimpse with the _DuMuX_ way of implementing mineralization and reaction processes. In the scope of this exercise, the setting of boundary conditions is revisited and a new reaction term is implemented.
## Problem set-up
The domain has a size of 20 x 15 m and contains a sealing aquitard in the middle. The aquitard is interrupted by a "fault zone" and thereby connects the upper (drinking water) aquifer and the lower CO2-storage aquifer. Initially, the domain is fully water saturated and biofilm is present in the lower CO2-storage aquifer. Calcium and urea are injected in the upper drinking water aquifer by means of a Neumann boundary condition. The remaining parts of the upper and the entire lower boundary are Neumann no-flow while on the right side a Dirichlet boundary conditions is applied (according to the initial values).
Disclaimer: Please note, that this is not a realistic scenario. One does not want to store gaseous CO2 in such a subcritical setting.
![](../extradoc/exercisebiomin_setup.png)
## Preparing the exercise
* Navigate to the directory `dumux-course/exercises/exercise-biomineralization`
### 1. Getting familiar with the code
Locate all the files you will need for this exercise
* The shared __main file__ : `exercisebiomin.cc`
* The __input file__: `exercisebiomin.input`
* The __problem file__ : `biominproblem.hh`
* The __spatial parameters file__: `biominspatialparams.hh`
Furthermore you will find the following folders:
* `chemistry`: Provides a way to formulate reactions in therms of source/sink terms.
* `components`: Provides some additional components e.g. biofilm and urea.
* `fluidsystems`: Stores headers containing data/methods on fluid mixtures of pure components.
* `solidsystems`: Stores headers containing data/methods on solid mixtures of the components.
To see more chemistry, components, fluidsystems and solidsystems implementations, have a look at the folder `dumux/material`.
__Special note on solidsystems:__
There are two types of solid components. Reactive and inert. For each reactive component one mass balances is solved. The inert components compose the "unchanging" (inert) rock matrix.
### 2. Implement a chemical equation
In the following the basic steps required to set the new chemical equation are outlined. Here, this is done in the __chemistry__ folder in the prepared file: 'simplebiominreactions.hh' starting in line 96.
Please be aware, that the chemistry file already provides some convenience functions (e.g. ``moleFracToMolality()``).
__Task__
Add a kinetic reaction by first calculating the current mass of available biofilm. Note that the density and volume fraction of biofilm are already defined for you.
$`\displaystyle mass_{biofilm} = \rho_{biofilm} * \phi_{biofilm}`$
Next, we want to implement the rate of ureolysis. This can be done with the following simplified equation:
$`\displaystyle r_{urea} = k_{urease} * Z_{urease,biofilm} * m_{urea} / (K_{urea} + m_{urea})`$,
where $`\displaystyle r_{urea}`$ is the rate of ureolysis, $`\displaystyle k_{urease}`$ the urease enzyme activity, $`\displaystyle m_{urea}`$ the molality of urea and $`\displaystyle K_{urea}`$ the half-saturation constant for urea for the ureolysis rate.
Note, that the urease concentration $`\displaystyle Z_{urease,biofilm}`$ has a kinetic term of urease production per biofilm :
$`\displaystyle Z_{urease,biofilm} = k_{urease,biofilm} * mass_{biofilm}`$
The last step is defining the source term for each component according to the chemical equation:
$`\displaystyle Ca^{2+} + CO(NH)_{2} + 2 H_{2}O --> 2 NH_{4}^{+} + CaCO_{3}`$
which is:
Calcium ion + Urea + 2 Water --> 2 Ammonium ion + Calcite
### 3. Make use of your newly created chemical equation
To enable your newly created chemical equation, the chemistry file has to be included in your problem file. This has to be done in line 35:
```c++
#include "chemistry/simplebiominreactions.hh" // chemical reactions
```
Additionally the TypeTag of your chemistry file needs to be set in the problem file (line 125):
```c++
using Chemistry = typename Dumux::SimpleBiominReactions<TypeTag>;
```
__Task__
Now the source/sink term can be updated in the problem file. You find it in line 343. You can access the newly created chemistry file and call the reactionSource()-function from it. Make sure to call the reactionSource()-function with the correct arguments. Return the updated source terms in the end.
The volume variables can be set using the element volume variables and the sub control volume.
```c++
const auto& volVars = elemVolVars[scv];
```
In order to compile and execute the program, change to the build-directory
```bash
cd build-cmake/exercises/exercise-biomineralization
```
and type
```bash
make exercisebiomin
./exercisebiomin exercisebiomin.input
```
### 4. Seal leakage pathway in the aquitard
In the input file, you will find some parameters concerning the mineralization process. We want to seal the leakage pathway in the aquitard. The leakage pathway is assumed to be sealed when the porosity is reduced to `0.07` or less.
__Task:__
Vary input parameters in the input file. Available for variation are the overall injection duration in days (`InjBioTime`), the initial biomass (`InitBiofilm`), the overall injection rate (`InjVolumeflux`) and the injected concentrations of urea and calcium (`concUrea` and `concCa`). When changing the concentrations, keep in mind that both urea and calcium are needed for the reaction and their mass ratio should be 2 calcium to 3 urea for a good molar ratio of 1 mol urea per 1 mol of calcium, see also the molar masses in the component files.
The result for the porosity should like like this:
![](../extradoc/exercisebiomin_porosityFinal.png)
### 5. CO2 injection to test aquitard integrity
Now the sealed aquitard is tested with a CO2-Injection into the lower CO2-storage aquifer.
__Task:__
Implement a new boundary condition from the left boundary, injecting CO2 from 2 m to 3 m. Make sure, that the injection time for the calcium and urea is finished. You can use the predefined value `gasFlux` directly and divide it by the molar mass of CO2.
Run two simulations and compare them side by side by creating two input files, or overwriting the input file in the command line:
```bash
./exercisebiomin -Problem.Name biominNoUrea -Injection.ConcUrea 0
```
The result for the biomineralization process after the CO2 injection should look like this:
![](../extradoc/exercisebiomin_injectionFinal.png)
### Bonus: Paraview Magic: Compare different results using Programmable Filter
In the last step, the manual comparison of the results can be quite difficult. Paraview offers the option to use programmable python filters. To use them, make sure two result files with different names are loaded. Mark both of them and click on `Filters --> Alphabetical --> Programmable Filter`. Now a new field opens on the left side. Copy the following python lines there:
```python
CO2_0 = inputs[0].CellData['x_n^CO2'];
CO2_1 = inputs[1].CellData['x_n^CO2'];
output.CellData.append(abs(CO2_0-CO2_1),'diffCO2');
```
Click `Apply` and select `diffCO2` as new output. You should now see the difference between the two result files. You can also change the output to a not absolute value by changing the last line to:
```python
output.CellData.append((CO2_0-CO2_1),'diffCO2');
```
// -*- 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 Tutorial problem for a fully coupled TwoPNCMineralization CC Tpfa model.
*/
#ifndef DUMUX_EXERCISE_FOUR_PROBLEM_HH
#define DUMUX_EXERCISE_FOUR_PROBLEM_HH
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/2pncmin/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "solidsystems/biominsolidphase.hh" // The biomineralization solid system
#include <dumux/material/components/co2tablereader.hh>
#include "fluidsystems/biomin.hh" // The biomineralization fluid system
// TODO: dumux-course-task
// include chemistry file here
#include "biominspatialparams.hh" // Spatially dependent parameters
namespace Dumux {
/*!
* \brief Tutorial problem for a fully coupled TwoPNCMineralization CC Tpfa model.
*/
//! Provides the precalculated tabulated values of CO2 density and enthalpy.
#include <dumux/material/components/co2tables.inc>
template <class TypeTag>
class ExerciseFourBioMinProblem;
namespace Properties
{
//! Create new type tag for the problem
NEW_TYPE_TAG(ExerciseFourBioMinTypeTag, INHERITS_FROM(TwoPNCMin, BioMinSpatialparams));
NEW_TYPE_TAG(ExerciseFourBioMinCCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, ExerciseFourBioMinTypeTag));
//! Set the problem property
SET_TYPE_PROP(ExerciseFourBioMinTypeTag, Problem, ExerciseFourBioMinProblem<TypeTag>);
//! Set grid and the grid creator to be used
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(ExerciseFourBioMinTypeTag, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::cube, Dune::nonconforming>);
#elif HAVE_UG
SET_TYPE_PROP(ExerciseFourBioMinTypeTag, Grid, Dune::UGGrid<2>);
#else
SET_TYPE_PROP(ExerciseFourBioMinTypeTag, Grid, Dune::YaspGrid<2>);
#endif // HAVE_DUNE_ALUGRID
//! Set the fluid system type
SET_PROP(ExerciseFourBioMinTypeTag, FluidSystem)
{
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using CO2Tables = Dumux::CO2Tables;
using H2OType = Components::TabulatedComponent<Components::H2O<Scalar>>;
public:
using type = FluidSystems::BioMin<Scalar, CO2Tables, H2OType>;
};
SET_PROP(ExerciseFourBioMinTypeTag, SolidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using type = SolidSystems::BiominSolidPhase<Scalar>;
};
SET_BOOL_PROP(ExerciseFourBioMinTypeTag, EnableFVGridGeometryCache, false);
SET_BOOL_PROP(ExerciseFourBioMinTypeTag, EnableGridVolumeVariablesCache, false);
SET_BOOL_PROP(ExerciseFourBioMinTypeTag, EnableGridFluxVariablesCache, false);
} // end namespace properties
/*!
*
* \brief Problem biomineralization (MICP) in an experimental setup.
*/
template <class TypeTag>
class ExerciseFourBioMinProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<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 SolidSystem = typename GET_PROP_TYPE(TypeTag, SolidSystem);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
// Grid dimension
enum
{
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using Element = typename GridView::template Codim<0>::Entity;
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
// TODO: dumux-course-task
// set the chemistry TypeTag
enum
{
numComponents = FluidSystem::numComponents,
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx, //Saturation
//Indices of the components
H2OIdx = FluidSystem::H2OIdx,
CO2Idx = FluidSystem::CO2Idx,
CaIdx = FluidSystem::CaIdx,
UreaIdx = FluidSystem::UreaIdx,
BiofilmIdx = SolidSystem::BiofilmIdx+numComponents,
CalciteIdx = SolidSystem::CalciteIdx+numComponents,
//Index of the primary component of wetting and nonwetting phase
conti0EqIdx = Indices::conti0EqIdx,
// Phase State
liquidPhaseOnly = Indices::firstPhaseOnly,
};
/*!
* \brief The constructor
*
* \param fvGridGeometry The finite volume grid geometry
*/
public:
ExerciseFourBioMinProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
#if !DUNE_VERSION_NEWER(DUNE_COMMON,2,7)
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-35);
#endif
name_ = getParam<std::string>("Problem.Name");
//initial values
densityW_ = getParam<Scalar>("Initial.InitDensityW");
initPressure_ = getParam<Scalar>("Initial.InitPressure");
initBiofilm_ = getParam<Scalar>("Initial.InitBiofilm");
//injection values
injVolumeflux_ = getParam<Scalar>("Injection.InjVolumeflux");
injBioTime_ = getParam<Scalar>("Injection.InjBioTime")*86400;
injCO2_ = getParam<Scalar>("Injection.InjCO2");
concCa_ = getParam<Scalar>("Injection.ConcCa");
concUrea_ = getParam<Scalar>("Injection.ConcUrea");
unsigned int codim = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod == DiscretizationMethod::box ? dim : 0;
Kxx_.resize(fvGridGeometry->gridView().size(codim));
Kyy_.resize(fvGridGeometry->gridView().size(codim));
//initialize the fluidsystem
FluidSystem::init(/*startTemp=*/288,
/*endTemp=*/ 303,
/*tempSteps=*/5,
/*startPressure=*/1e4,
/*endPressure=*/1e6,
/*pressureSteps=*/500);
}
void setTime(Scalar time)
{
time_ = time;
}
/*!
* \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 300 degrees Kelvin.
*/
Scalar temperature() const
{
return 300; // [K]
}
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
// we don't set any BCs for the solid phases
Scalar xmax = this->fvGridGeometry().bBoxMax()[0];
// set all Neumann
bcTypes.setAllNeumann();
// set right boundary to Dirichlet
if(globalPos[0] > xmax - eps_)
bcTypes.setAllDirichlet();
return bcTypes;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
// recycle initialAtPos() for Dirichlet values
return initialAtPos(globalPos);
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param values Stores the Neumann values for the conservation equations in
* \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
* \param globalPos The position of the integration point of the boundary segment.
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
NumEqVector values(0.0);
Scalar waterFlux = injVolumeflux_; // divide by area if area not 1! [m/s]
Scalar gasFlux = injCO2_; // divide by area if area not 1! [m/s]
// Set values for Ca + urea injection above aquitard.
// Use negative values for injection.
if(globalPos[0] < eps_
&& globalPos[1] > 11.0 + eps_
&& globalPos[1] < 12.0 + eps_
&& time_ < injBioTime_)
{
values[conti0EqIdx + H2OIdx] = - waterFlux * densityW_ / FluidSystem::molarMass(H2OIdx);
values[conti0EqIdx + CaIdx] = - waterFlux * concCa_ / FluidSystem::molarMass(CaIdx);
values[conti0EqIdx + UreaIdx] = - waterFlux * concUrea_ / FluidSystem::molarMass(UreaIdx);
}
// TODO: dumux-course-task
// Set CO2 injection below aquitard after the injBioTime
else
{
values = 0.0; //mol/m²/s
}
return values;
}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param globalPos The global position
*
* For this method, the \a values parameter stores primary
* variables.
*/
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables priVars(0.0);
priVars.setState(liquidPhaseOnly);
Scalar hmax = this->fvGridGeometry().bBoxMax()[1];
priVars[pressureIdx] = initPressure_ - (globalPos[1] - hmax)*densityW_*9.81;
priVars[switchIdx] = 0.0;
priVars[CaIdx] = 0.0;
priVars[UreaIdx] = 0.0;
priVars[CalciteIdx] = 0.0;
priVars[BiofilmIdx] = 0.0;
// set Biofilm presence only in aquitard zone and below
if (globalPos[1] < 10 + eps_)
{
priVars[switchIdx] = 0.0;
priVars[BiofilmIdx] = initBiofilm_; // [m^3/m^3]
}
return priVars;
}
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
* This is the method for the case where the source term is
* potentially solution dependent and requires some quantities that
* are specific to the fully-implicit method.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param scv The sub control volume
*
* For this method, the return parameter stores the conserved quantity rate
* generated or annihilate per volume unit. Positive values mean
* that the conserved quantity is created, negative ones mean that it vanishes.
* E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
*/
NumEqVector source(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
NumEqVector source(0.0);
// TODO: dumux-course-task
// set Chemistry
// set VolumeVariables
// call reactionSource() in chemistry
return source;
}
const std::vector<Scalar>& getKxx()
{
return Kxx_;
}
const std::vector<Scalar>& getKyy()
{
return Kyy_;
}
/*!
* \brief Adds additional VTK output data to the VTKWriter.
* Function is called by the output module on every write.
*/
void updateVtkOutput(const SolutionVector& curSol)
{
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
const auto elemSol = elementSolution(element, curSol, this->fvGridGeometry());
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
VolumeVariables volVars;
volVars.update(elemSol, *this, element, scv);
const auto dofIdxGlobal = scv.dofIndex();
Kxx_[dofIdxGlobal] = volVars.permeability()[0][0];
Kyy_[dofIdxGlobal] = volVars.permeability()[1][1];
}
}
}
private:
static constexpr Scalar eps_ = 1e-6;
Scalar densityW_;
Scalar initPressure_;
Scalar initBiofilm_;
Scalar injVolumeflux_;
Scalar injBioTime_;
Scalar injCO2_;
Scalar concCa_;
Scalar concUrea_;
std::vector<double> Kxx_;
std::vector<double> Kyy_;
std::string name_;
Scalar time_;
};
} //end namespace
#endif
// -*- 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/>. *
*****************************************************************************/
/*!