diff --git a/tutorial/ex1/README.md b/tutorial/ex1/README.md
index 0548a993fcd5eb501eca1adde885759a42d68f59..a5757ce6834e0802f8f1018e84df39e433704606 100644
--- a/tutorial/ex1/README.md
+++ b/tutorial/ex1/README.md
@@ -66,10 +66,13 @@ EntryPressureCoarse = 1e4
### Task 4: Runtime parameters
<hr>
-The injection rate is currently hard-coded in `injection2pproblem.hh`.
+The injection rate is currently hard-coded in `injection2p2cproblem.hh` to $`1e-4 kg/(s m^2)`$.
```c++
-code snippet
+ // set the Neumann values for the Nitrogen component balance
+ // convert from units kg/(s*m^2) to mole/(s*m^2)
+values[Indices::contiNEqIdx] = -1e-4/FluidSystem::molarMass(FluidSystem::nCompIdx);
+values[Indices::contiWEqIdx] = 0.0;
```
We want to be able to set it at runtime. To this end,
@@ -77,7 +80,7 @@ We want to be able to set it at runtime. To this end,
```c++
// read the injection rate from the input file at run time
-const auto injectionRate = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, <TYPE>, <GROUPNAME>, <PARAMNAME>);
+const auto totalAreaSpecificInflow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, <TYPE>, <GROUPNAME>, <PARAMNAME>);
```
* Replace
@@ -86,9 +89,14 @@ const auto injectionRate = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, <TYPE>, <GROUPN
* `<GROUPNAME>` is the group in the input file
* `<PARAMNAME>` is the name of the parameter in the input file
-Note that due to the way the macro works, the names are specified as plain text without "quotation marks".
+Note that due to the way the macro works, the names are specified as plain text without "quotation marks". Follow the instructions given as a
-* Check the influence of that parameter on the simulation result by rerunning the simulation with different injection rates.
+```c++
+// TODO: dumux-course-task
+```
+in the `injection2p2cproblem.hh` file and also remember to also set the parameter totalAreaSpecificInflow in the input file.
+
+* Check the influence of that parameter on the simulation result by rerunning the simulation with different injection rates. Remember to also set the parameter totalAreaSpecificInflow in the input file.
Again, you don't need to recompile the program.
diff --git a/tutorial/ex1/exercise1.input b/tutorial/ex1/exercise1.input
index a4e055be8c3a8a96c0e5f29370deb9ee11dd1444..4f9a57e69a8ccd12d4831c321115806d7c7a92f5 100755
--- a/tutorial/ex1/exercise1.input
+++ b/tutorial/ex1/exercise1.input
@@ -1,6 +1,6 @@
[TimeManager]
-DtInitial = 250 # [s]
-TEnd = 1e7 # [s]
+DtInitial = 3600 # in seconds
+TEnd = 3.154e9 # in seconds, i.e ten years
[Grid]
LowerLeft = 0 0
@@ -8,11 +8,17 @@ UpperRight = 60 40
Cells = 24 16
[Problem]
-MaxDepth = 2700.0
+Name = injection
+OnlyPlotMaterialLaws = true
+AquiferDepth = 2700.0 # m
+InjectionDuration = 2.628e6 # in seconds, i.e. one month
-[Newton]
-WriteConvergence = 1 # write convergence behaviour to disk?
+#TODO: dumux-course-task:
+#set totalAreaSpecificInflow
[SpatialParams]
-EntryPressureFine = 4.5e4
-EntryPressureCoarse = 1e4
+PermeabilityAquitard = 1e-15 # m^2
+EntryPressureAquitard = 4.5e4 # Pa
+PermeabilityAquifer = 1e-12 # m^2
+EntryPressureAquifer = 1e4 # Pa
+
diff --git a/tutorial/ex1/injection2p2cproblem.hh b/tutorial/ex1/injection2p2cproblem.hh
index a63bc73d1c3a9d857b47ec2eb599050b332c8722..aada803f77f5bcfbcb06a4e61122de6a7a2c560a 100644
--- a/tutorial/ex1/injection2p2cproblem.hh
+++ b/tutorial/ex1/injection2p2cproblem.hh
@@ -64,21 +64,19 @@ SET_BOOL_PROP(Injection2p2cProblem, UseMoles, true);
* \brief Problem where air is injected under a low permeable layer in a depth of 2700m.
*
* The domain is sized 60m times 40m and consists of two layers, a moderately
- * permeable one (\f$ K=10e-12\f$) for \f$ y<22m\f$ and one with a lower permeablility (\f$ K=10e-13\f$)
+ * permeable one (\f$ K=10e-12\f$) and one with a lower permeablility (\f$ K=10e-13\f$)
* in the rest of the domain.
*
- * A mixture of Nitrogen and Water vapor, which is composed according to the prevailing conditions (temperature, pressure)
- * enters a water-filled aquifer. This is realized with a solution-dependent Neumann boundary condition at the right boundary
+ * Nitrogen is injected into a water-filled aquifer through a well. First, we inject for one month.
+ * Then, we continue simulating the development of the nitrogen plume for 10 years.This is realized with a solution-dependent Neumann boundary condition at the right boundary
* (\f$ 7m<y<15m\f$). The aquifer is situated 2700m below sea level. The injected fluid phase migrates upwards due to buoyancy.
* It accumulates and partially enters the lower permeable aquitard.
*
- * The model is able to use either mole or mass fractions. The property useMoles can be set to either true or false in the
- * problem file. Make sure that the according units are used in the problem setup. The default setting for useMoles is true.
- *
+ * The default setting for useMoles is true, i.e. each component is balaced in units of mole.
* This problem uses the \ref TwoPTwoCModel.
*
* To run the simulation execute the following line in shell:
- * <tt>./exercise1_2p2c -parameterFile exercise1.input>
+ * <tt>./exercise1_2p2c -ParameterFile exercise1.input>
*/
template <class TypeTag>
class Injection2p2cProblem : public ImplicitPorousMediaProblem<TypeTag>
@@ -89,45 +87,16 @@ class Injection2p2cProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- enum {
- // Grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
+ // grid world dimension
+ static constexpr auto dimWorld = GridView::dimensionworld;
- // copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- enum {
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
-
-
- wCompIdx = FluidSystem::wCompIdx,
- nCompIdx = FluidSystem::nCompIdx,
-
- contiH2OEqIdx = Indices::contiWEqIdx,
- contiN2EqIdx = Indices::contiNEqIdx
- };
-
-
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<dim>::Entity Vertex;
- typedef typename GridView::Intersection Intersection;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
-
- //! property that defines whether mole or mass fractions are used
- static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-
public:
/*!
* \brief The constructor
@@ -139,8 +108,6 @@ public:
const GridView &gridView)
: ParentType(timeManager, gridView)
{
- maxDepth_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.MaxDepth);
-
// initialize the tables of the fluid system
FluidSystem::init(/*tempMin=*/273.15,
/*tempMax=*/423.15,
@@ -149,15 +116,16 @@ public:
/*pMax=*/30e6,
/*numP=*/300);
- //stateing in the console whether mole or mass fractions are used
- if(useMoles)
- {
- std::cout<<"problem uses mole-fractions"<<std::endl;
- }
- else
- {
- std::cout<<"problem uses mass-fractions"<<std::endl;
- }
+ // name of the problem and output file
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+ // depth of the aquifer, units: m
+ aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth);
+ // the duration of the injection, units: second
+ injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration);
+
+ // TODO: dumux-course-task
+ // Get the specific inflow of 1e-4 kg/(s m^2) from the input file (totalAreaSpecificInflow_) here as it is done for the injectionDuration_.
+
}
/*!
@@ -169,8 +137,8 @@ public:
{
// Calculate storage terms
PrimaryVariables storageW, storageN;
- this->model().globalPhaseStorage(storageW, wPhaseIdx);
- this->model().globalPhaseStorage(storageN, nPhaseIdx);
+ this->model().globalPhaseStorage(storageW, Indices::wPhaseIdx);
+ this->model().globalPhaseStorage(storageN, Indices::nPhaseIdx);
// Write mass balance information for rank 0
if (this->gridView().comm().rank() == 0) {
@@ -190,8 +158,8 @@ public:
*
* This is used as a prefix for files generated by the simulation.
*/
- const std::string name() const
- { return "injection-2p2c"; }
+ std::string name() const
+ { return name_+"-2p2c"; }
/*!
* \brief Returns the temperature \f$ K \f$
@@ -231,8 +199,11 @@ public:
void boundaryTypesAtPos(BoundaryTypes &values,
const GlobalPosition &globalPos) const
{
- if (globalPos[0] < eps_)
+ if (globalPos[dimWorld-1] < eps_)
values.setAllDirichlet();
+
+ // and Neuman boundary conditions everywhere else
+ // note that we don't differentiate between Neumann and Robin boundary types
else
values.setAllNeumann();
}
@@ -247,49 +218,34 @@ public:
*/
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ initialAtPos(values, globalPos);
}
/*!
- * \brief Evaluates the boundary conditions for a Neumann
- * boundary segment in dependency on the current solution.
+ * \brief Evaluates 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 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.
+ * \param globalPos The globalPosition of the boundary interface
*/
- void solDependentNeumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &intersection,
- const int scvIdx,
- const int boundaryFaceIdx,
- const ElementVolumeVariables &elemVolVars) const
+ void neumannAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
{
+ // initialize values to zero, i.e. no-flow Neumann boundary conditions
values = 0;
- GlobalPosition globalPos;
- if (isBox)
- globalPos = element.geometry().corner(scvIdx);
- else
- globalPos = intersection.geometry().center();
-
- Scalar injectedPhaseMass = 1e-3;
- Scalar moleFracW = elemVolVars[scvIdx].moleFraction(nPhaseIdx, wCompIdx);
- if (globalPos[1] < 15 + eps_ && globalPos[1] > 7 -eps_) {
- values[contiN2EqIdx] = -(1-moleFracW)*injectedPhaseMass/FluidSystem::molarMass(nCompIdx); //mole/(m^2*s) -> kg/(s*m^2)
- values[contiH2OEqIdx] = -moleFracW*injectedPhaseMass/FluidSystem::molarMass(wCompIdx); //mole/(m^2*s) -> kg/(s*m^2)
- }
+ //if we are inside the injection zone set inflow Neumann boundary conditions
+ if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_
+ && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0])
+ {
+ // TODO: dumux-course-task
+ //instead of setting -1e-4 here directly use totalAreaSpecificInflow_ in the computation
+
+ // set the Neumann values for the Nitrogen component balance
+ // convert from units kg/(s*m^2) to mole/(s*m^2)
+ values[Indices::contiNEqIdx] = -1e-4/FluidSystem::molarMass(FluidSystem::nCompIdx);
+ values[Indices::contiWEqIdx] = 0.0;
+ }
}
// \}
@@ -308,7 +264,23 @@ public:
*/
void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ // get the water density at atmospheric conditions
+ const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5);
+
+ // assume an intially hydrostatic liquid pressure profile
+ // note: we subtract rho_w*g*h because g is defined negative
+ const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
+
+ // initially we have some nitrogen dissolved
+ // saturation mole fraction would be
+ // moleFracLiquidN2 = (pw + pc + p_vap^sat)/henry;
+ const Scalar moleFracLiquidN2 = pw*0.95/BinaryCoeff::H2O_N2::henry(temperature());
+
+ // note that because we start with a single phase system the primary variables
+ // are pl and x^w_N2. This will switch as soon after we start injecting to a two
+ // phase system so the primary variables will be pl and Sn (non-wetting saturation).
+ values[Indices::switchIdx] = moleFracLiquidN2;
+ values[Indices::pressureIdx] = pw;
}
/*!
@@ -322,44 +294,12 @@ public:
// \}
private:
- /*!
- * \brief Evaluates the initial values for a control volume
- *
- * The internal method for the initial condition
- *
- * \param values Stores the initial values for the conservation equations in
- * \f$ [ \textnormal{unit of primary variables} ] \f$
- * \param globalPos The global position
- */
- void initial_(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1e5);
-
- Scalar pl = 1e5 - densityW*this->gravity()[1]*(maxDepth_ - globalPos[1]);
- Scalar moleFracLiquidN2 = pl*0.95/BinaryCoeff::H2O_N2::henry(temperature());
- Scalar moleFracLiquidH2O = 1.0 - moleFracLiquidN2;
-
- Scalar meanM =
- FluidSystem::molarMass(wCompIdx)*moleFracLiquidH2O +
- FluidSystem::molarMass(nCompIdx)*moleFracLiquidN2;
- if(useMoles)
- {
- //mole-fraction formulation
- values[Indices::switchIdx] = moleFracLiquidN2;
- }
- else
- {
- //mass fraction formulation
- Scalar massFracLiquidN2 = moleFracLiquidN2*FluidSystem::molarMass(nCompIdx)/meanM;
- values[Indices::switchIdx] = massFracLiquidN2;
- }
- values[Indices::pressureIdx] = pl;
- }
-
-
static constexpr Scalar eps_ = 1e-6;
- Scalar maxDepth_;
+ std::string name_; //! Problem name
+ Scalar aquiferDepth_; //! Depth of the aquifer in m
+ Scalar injectionDuration_; //! Duration of the injection in seconds
+ //TODO: dumux-course-task
+ //define the Scalar totalAreaSpecificInflow_ here
};
} //end namespace
diff --git a/tutorial/ex1/injection2pniproblem.hh b/tutorial/ex1/injection2pniproblem.hh
index 6a964f537f371d9a465f75327b8cb3ea2c62b815..f7453851848625ae1af19337cac0e930887efd4d 100644
--- a/tutorial/ex1/injection2pniproblem.hh
+++ b/tutorial/ex1/injection2pniproblem.hh
@@ -137,8 +137,6 @@ public:
InjectionProblem2PNI(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
- maxDepth_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.MaxDepth);
-
// initialize the tables of the fluid system
FluidSystem::init(/*tempMin=*/273.15,
/*tempMax=*/423.15,
@@ -146,6 +144,14 @@ public:
/*pMin=*/0.0,
/*pMax=*/30e6,
/*numP=*/300);
+
+ // name of the problem and output file
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+ // depth of the aquifer, units: m
+ aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth);
+ // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2)
+ // the duration of the injection, units: second
+ injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration);
}
/*!
@@ -212,8 +218,14 @@ public:
*/
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- Scalar densityW = 1000.0;
- values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
+ // assume a constant density
+ const Scalar densityW = 1000;
+
+ // assume an intially hydrostatic liquid pressure profile
+ // note: we subtract rho_w*g*h because g is defined negative
+ const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
+
+ values[Indices::pressureIdx] = pw;
values[saturationIdx] = 0.0;
/*!
* TODO:dumux-course-task:
@@ -223,38 +235,25 @@ public:
}
/*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
+ * \brief Evaluates 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 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
- *
- * The \a values store the mass flux of each phase normal to the boundary.
- * Negative values indicate an inflow.
+ * \param globalPos The globalPosition of the boundary interface
*/
- void neumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &intersection,
- int scvIdx,
- int boundaryFaceIdx) const
+ void neumannAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
{
+ // initialize values to zero, i.e. no-flow Neumann boundary conditions
+ values = 0;
values = 0;
- GlobalPosition globalPos;
- if (isBox)
- globalPos = element.geometry().corner(scvIdx);
- else
- globalPos = intersection.geometry().center();
-
- if (globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_) {
+ if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_
+ && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0])
+ {
// inject air. negative values mean injection
- values[contiNEqIdx] = -1e-3; // kg/(s*m^2)
+ values[Indices::contiNEqIdx] = -1e-4; // kg/(s*m^2)
+ values[Indices::contiWEqIdx] = 0.0;
}
}
@@ -275,8 +274,14 @@ public:
*/
void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- Scalar densityW = 1000.0;
- values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
+ // assume a constant density
+ const Scalar densityW = 1000;
+
+ // assume an intially hydrostatic liquid pressure profile
+ // note: we subtract rho_w*g*h because g is defined negative
+ const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
+
+ values[Indices::pressureIdx] = pw;
values[saturationIdx] = 0.0;
/*!
@@ -289,9 +294,10 @@ public:
private:
- Scalar maxDepth_;
- static constexpr Scalar eps_ = 1.5e-7;
- std::string name_;
+ static constexpr Scalar eps_ = 1e-6;
+ std::string name_; //! Problem name
+ Scalar aquiferDepth_; //! Depth of the aquifer in m
+ Scalar injectionDuration_; //! Duration of the injection in seconds
};
} //end namespace
diff --git a/tutorial/ex1/injection2pproblem.hh b/tutorial/ex1/injection2pproblem.hh
index fbf085c328505ec389181ae7980d3680d4a5c295..c990ffb09b4c2b41aa55834bda0c9a2c7180f0ef 100644
--- a/tutorial/ex1/injection2pproblem.hh
+++ b/tutorial/ex1/injection2pproblem.hh
@@ -123,8 +123,6 @@ public:
InjectionProblem2P(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
- maxDepth_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.MaxDepth);
-
// initialize the tables of the fluid system
FluidSystem::init(/*tempMin=*/273.15,
/*tempMax=*/423.15,
@@ -132,6 +130,14 @@ public:
/*pMin=*/0.0,
/*pMax=*/30e6,
/*numP=*/300);
+
+ // name of the problem and output file
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+ // depth of the aquifer, units: m
+ aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth);
+ // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2)
+ // the duration of the injection, units: second
+ injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration);
}
/*!
@@ -204,50 +210,34 @@ public:
*/
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1e5);
- values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
- values[saturationIdx] = 0.0;
+ initialAtPos(values, globalPos);
}
/*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
+ * \brief Evaluates 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 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
- *
- * The \a values store the mass flux of each phase normal to the boundary.
- * Negative values indicate an inflow.
+ * \param globalPos The globalPosition of the boundary interface
*/
- void neumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &intersection,
- int scvIdx,
- int boundaryFaceIdx) const
+ void neumannAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
{
+ // initialize values to zero, i.e. no-flow Neumann boundary conditions
+ values = 0;
values = 0;
- GlobalPosition globalPos;
- if (isBox)
- globalPos = element.geometry().corner(scvIdx);
- else
- globalPos = intersection.geometry().center();
-
- if (globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_) {
+ if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_
+ && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0])
+ {
// inject air. negative values mean injection
- values[contiNEqIdx] = -1e-3; // kg/(s*m^2)
+ values[Indices::contiNEqIdx] = -1e-4; // kg/(s*m^2)
+ values[Indices::contiWEqIdx] = 0.0;
}
}
// \}
-
/*!
* \name Volume terms
*/
@@ -262,15 +252,24 @@ public:
*/
void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1e5);
- values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
+
+ // get the water density at atmospheric conditions
+ const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5);
+
+ // assume an intially hydrostatic liquid pressure profile
+ // note: we subtract rho_w*g*h because g is defined negative
+ const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
+
+ values[Indices::pressureIdx] = pw;
values[saturationIdx] = 0.0;
}
// \}
private:
- Scalar maxDepth_;
- static constexpr Scalar eps_ = 1.5e-7;
+ static constexpr Scalar eps_ = 1e-6;
+ std::string name_; //! Problem name
+ Scalar aquiferDepth_; //! Depth of the aquifer in m
+ Scalar injectionDuration_; //! Duration of the injection in seconds
};
} //end namespace
diff --git a/tutorial/ex1/injection2pspatialparams.hh b/tutorial/ex1/injection2pspatialparams.hh
index fc85aa4c21390abccc302f0b7d22148e83dc7caa..95071e887b77f6d1f5775204ff7322d69c82b224 100644
--- a/tutorial/ex1/injection2pspatialparams.hh
+++ b/tutorial/ex1/injection2pspatialparams.hh
@@ -31,6 +31,9 @@
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/io/plotmateriallaw.hh>
+
#include <dumux/porousmediumflow/2p2c/implicit/properties.hh>
namespace Dumux
@@ -93,34 +96,34 @@ public:
InjectionSpatialParams(const GridView &gridView)
: ParentType(gridView)
{
- layerBottom_ = 25.0;
+ aquiferHeightFromBottom_ = 30.0;
// intrinsic permeabilities
- fineK_ = 1e-13;
- coarseK_ = 1e-12;
+ aquitardK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityAquitard);
+ aquiferK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityAquifer);
// porosities
- finePorosity_ = 0.2;
- coarsePorosity_ = 0.4;
-
- //materialLawParams
- fineEntryPressure_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureFine);
- coarseEntryPressure_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureCoarse);
-
- // heat conductivity of granite
- lambdaSolid_ = 2.8;
+ aquitardPorosity_ = 0.2;
+ aquiferPorosity_ = 0.4;
// residual saturations
- fineMaterialParams_.setSwr(0.2);
- fineMaterialParams_.setSnr(0.0);
- coarseMaterialParams_.setSwr(0.2);
- coarseMaterialParams_.setSnr(0.0);
+ aquitardMaterialParams_.setSwr(0.2);
+ aquitardMaterialParams_.setSnr(0.0);
+ aquiferMaterialParams_.setSwr(0.2);
+ aquiferMaterialParams_.setSnr(0.0);
// parameters for the Brooks-Corey law
- fineMaterialParams_.setPe(fineEntryPressure_);
- coarseMaterialParams_.setPe(coarseEntryPressure_);
- fineMaterialParams_.setLambda(2.0);
- coarseMaterialParams_.setLambda(2.0);
+ aquitardMaterialParams_.setPe(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureAquitard));
+ aquiferMaterialParams_.setPe(GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureAquifer));
+ aquitardMaterialParams_.setLambda(2.0);
+ aquiferMaterialParams_.setLambda(2.0);
+
+ // plot the material laws using gnuplot and exit
+ if (GET_RUNTIME_PARAM(TypeTag, bool, Problem.OnlyPlotMaterialLaws))
+ {
+ plotMaterialLaws();
+ exit(0);
+ }
}
/*!
@@ -135,9 +138,9 @@ public:
const int scvIdx) const
{
const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
- if (isFineMaterial_(globalPos))
- return fineK_;
- return coarseK_;
+ if (isInAquitard_(globalPos))
+ return aquitardK_;
+ return aquiferK_;
}
/*!
@@ -152,9 +155,9 @@ public:
const int scvIdx) const
{
const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
- if (isFineMaterial_(globalPos))
- return finePorosity_;
- return coarsePorosity_;
+ if (isInAquitard_(globalPos))
+ return aquitardPorosity_;
+ return aquiferPorosity_;
}
@@ -171,9 +174,9 @@ public:
const int scvIdx) const
{
const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
- if (isFineMaterial_(globalPos))
- return fineMaterialParams_;
- return coarseMaterialParams_;
+ if (isInAquitard_(globalPos))
+ return aquitardMaterialParams_;
+ return aquiferMaterialParams_;
}
/*!
@@ -228,24 +231,37 @@ public:
// return lambdaSolid_;
// }
+ /*!
+ * \brief Creates a gnuplot output of the pc-Sw curve
+ */
+ void plotMaterialLaws()
+ {
+ PlotMaterialLaw<TypeTag> plotMaterialLaw;
+ GnuplotInterface<Scalar> gnuplot;
+ plotMaterialLaw.addpcswcurve(gnuplot, aquitardMaterialParams_, 0.2, 1.0, "upper layer (fine, aquitard)", "w lp");
+ plotMaterialLaw.addpcswcurve(gnuplot, aquiferMaterialParams_, 0.2, 1.0, "lower layer (coarse, aquifer)", "w l");
+ gnuplot.setOption("set xrange [0:1]");
+ gnuplot.setOption("set label \"residual\\nsaturation\" at 0.1,100000 center");
+ gnuplot.plot("pc-Sw");
+ }
+
private:
- bool isFineMaterial_(const GlobalPosition &globalPos) const
- { return globalPos[dimWorld-1] > layerBottom_; }
- Scalar fineK_;
- Scalar coarseK_;
- Scalar layerBottom_;
+ static constexpr Scalar eps_ = 1e-6;
+
+ bool isInAquitard_(const GlobalPosition &globalPos) const
+ { return globalPos[dimWorld-1] > aquiferHeightFromBottom_ + eps_; }
- Scalar finePorosity_;
- Scalar coarsePorosity_;
+ Scalar aquitardK_;
+ Scalar aquiferK_;
+ Scalar aquiferHeightFromBottom_;
- Scalar lambdaSolid_;
- Scalar fineEntryPressure_;
- Scalar coarseEntryPressure_;
+ Scalar aquitardPorosity_;
+ Scalar aquiferPorosity_;
- MaterialLawParams fineMaterialParams_;
- MaterialLawParams coarseMaterialParams_;
+ MaterialLawParams aquitardMaterialParams_;
+ MaterialLawParams aquiferMaterialParams_;
};
}
diff --git a/tutorial/ex2/injection2p2cproblem.hh b/tutorial/ex2/injection2p2cproblem.hh
index cf45851a5b31b92b2f3ba3afbdac0b83475d3c65..6b611593e14ed25de875bcf24d9ec1efe3e92737 100644
--- a/tutorial/ex2/injection2p2cproblem.hh
+++ b/tutorial/ex2/injection2p2cproblem.hh
@@ -74,7 +74,7 @@ SET_BOOL_PROP(Injection2p2cProblem, UseMoles, true);
* \brief Problem where air is injected under a low permeable layer in a depth of 2700m.
*
* The domain is sized 60m times 40m and consists of two layers, a moderately
- * permeable one (\f$ K=10e-12\f$) for \f$ y<22m\f$ and one with a lower permeablility (\f$ K=10e-13\f$)
+ * permeable one for \f$ y<22m\f$ and one with a lower permeablility
* in the rest of the domain.
*
* Nitrogen is injected into a water-filled aquifer through a well. First, we inject for one month.
diff --git a/tutorial/solution/ex1/injection2p2cproblem.hh b/tutorial/solution/ex1/injection2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e03f339f15072819f16f91b52a31af1f75173640
--- /dev/null
+++ b/tutorial/solution/ex1/injection2p2cproblem.hh
@@ -0,0 +1,310 @@
+// -*- 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 Problem where air is injected under a low permeable layer in a depth of 2700m.
+ */
+#ifndef DUMUX_INJECTION_2P2C_PROBLEM_HH
+#define DUMUX_INJECTION_2P2C_PROBLEM_HH
+
+#include <dumux/porousmediumflow/2p2c/implicit/model.hh>
+#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/material/fluidsystems/h2on2.hh>
+
+#include "injection2p2cspatialparams.hh"
+
+namespace Dumux
+{
+
+// foward declaration
+template <class TypeTag>
+class Injection2p2cProblem;
+
+// setup property TypeTag
+namespace Properties
+{
+
+// inherit from MyLocalResidualParams
+NEW_TYPE_TAG(Injection2p2cProblem, INHERITS_FROM(TwoPTwoC, InjectionSpatialParams));
+NEW_TYPE_TAG(Injection2p2cBoxProblem, INHERITS_FROM(BoxModel, Injection2p2cProblem));
+NEW_TYPE_TAG(Injection2p2pcCCProblem, INHERITS_FROM(CCModel, Injection2p2cProblem));
+
+// Set the grid type
+SET_TYPE_PROP(Injection2p2cProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(Injection2p2cProblem, Problem, Injection2p2cProblem<TypeTag>);
+
+// Set fluid configuration
+SET_TYPE_PROP(Injection2p2cProblem,
+ FluidSystem,
+ FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), true /*useComplexRelations*/>);
+
+// Define whether mole(true) or mass (false) fractions are used
+SET_BOOL_PROP(Injection2p2cProblem, UseMoles, true);
+
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \ingroup ImplicitTestProblems
+ * \brief Problem where air is injected under a low permeable layer in a depth of 2700m.
+ *
+ * The domain is sized 60m times 40m and consists of two layers, a moderately
+ * permeable one for \f$ y<22m\f$ and one with a lower permeablility
+ * in the rest of the domain.
+ *
+ * Nitrogen is injected into a water-filled aquifer through a well. First, we inject for one month.
+ * Then, we continue simulating the development of the nitrogen plume for 10 years.
+ * This is realized with a Neumann boundary condition at the right boundary
+ * (\f$ 7m<y<15m\f$). The aquifer is situated 2700m below sea level (the depth can be changed through the input file).
+ * The injected fluid phase migrates upwards due to buoyancy.
+ * It accumulates and partially enters the top layer lower permeable aquitard.
+ *
+ * The default setting for useMoles is true, i.e. each component is balaced in units of mole.
+ * The property useMoles can be set to either true or false in the
+ * problem file. If you change this, make sure that the according units are used in the problem setup.
+ *
+ * This problem uses the \ref TwoPTwoCModel.
+ */
+template <class TypeTag>
+class Injection2p2cProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+ typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ // grid world dimension
+ static constexpr auto dimWorld = GridView::dimensionworld;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ Injection2p2cProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ // initialize the tables of the fluid system
+ FluidSystem::init(/*tempMin=*/273.15,
+ /*tempMax=*/423.15,
+ /*numTemp=*/50,
+ /*pMin=*/0.0,
+ /*pMax=*/30e6,
+ /*numP=*/300);
+
+ // name of the problem and output file
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+ // depth of the aquifer, units: m
+ aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth);
+ // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2)
+ totalAreaSpecificInflow_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, TotalAreaSpecificInflow);
+ // the duration of the injection, units: second
+ injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration);
+ }
+
+ /*!
+ * \brief User defined output after the time integration
+ *
+ * Will be called diretly after the time integration
+ */
+ void postTimeStep()
+ {
+ // Calculate storage terms
+ PrimaryVariables storageW, storageN;
+ this->model().globalPhaseStorage(storageW, Indices::wPhaseIdx);
+ this->model().globalPhaseStorage(storageN, Indices::nPhaseIdx);
+
+ // Write mass balance information for rank 0
+ if (this->gridView().comm().rank() == 0)
+ {
+ std::cout <<"Storage: wetting=[" << storageW << "]"
+ << " nonwetting=[" << storageN << "]" << std::endl;
+ }
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns 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 in \f$ K \f$
+ */
+ Scalar temperature() const
+ { return 273.15 + 30; }
+
+ /*!
+ * \brief Returns the source term
+ *
+ * \param values Stores the source values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+ * \param globalPos The global position
+ */
+ void sourceAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ { values = 0; }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment
+ *
+ * \param values Stores the value of the boundary type
+ * \param globalPos The global position
+ */
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
+ {
+ // Set Dirichlet at the bottom of the domain
+ if (globalPos[dimWorld-1] < eps_)
+ {
+ values.setAllDirichlet();
+ }
+
+ // and Neuman boundary conditions everywhere else
+ // note that we don't differentiate between Neumann and Robin boundary types
+ else
+ {
+ values.setAllNeumann();
+ }
+ }
+
+ /*!
+ * \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
+ */
+ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ { initialAtPos(values, globalPos); }
+
+ /*!
+ * \brief Evaluates 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 globalPosition of the boundary interface
+ */
+ void neumannAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ // initialize values to zero, i.e. no-flow Neumann boundary conditions
+ values = 0;
+
+ //if we are inside the injection zone set inflow Neumann boundary conditions
+ if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_
+ && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0])
+ {
+ // set the Neumann values for the Nitrogen component balance
+ // convert from units kg/(s*m^2) to mole/(s*m^2)
+ values[Indices::contiNEqIdx] = -totalAreaSpecificInflow_/FluidSystem::molarMass(FluidSystem::nCompIdx);
+ values[Indices::contiWEqIdx] = 0.0;
+ }
+ }
+
+ // \}
+
+ /*!
+ * \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
+ */
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ // get the water density at atmospheric conditions
+ const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5);
+
+ // assume an intially hydrostatic liquid pressure profile
+ // note: we subtract rho_w*g*h because g is defined negative
+ const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
+
+ // initially we have some nitrogen dissolved
+ // saturation mole fraction would be
+ // moleFracLiquidN2 = (pw + pc + p_vap^sat)/henry;
+ const Scalar moleFracLiquidN2 = pw*0.95/BinaryCoeff::H2O_N2::henry(temperature());
+
+ // note that because we start with a single phase system the primary variables
+ // are pl and x^w_N2. This will switch as soon after we start injecting to a two
+ // phase system so the primary variables will be pl and Sn (non-wetting saturation).
+ values[Indices::switchIdx] = moleFracLiquidN2;
+ values[Indices::pressureIdx] = pw;
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ *
+ * \param globalPos The global position
+ * \note we start with a single phase system
+ */
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const
+ { return Indices::wPhaseOnly; }
+
+ // \}
+
+ //! If we should write restart files
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+private:
+ static constexpr Scalar eps_ = 1e-6;
+ std::string name_; //! Problem name
+ Scalar aquiferDepth_; //! Depth of the aquifer in m
+ Scalar totalAreaSpecificInflow_; //! Area specific inflow rate in mole/(s*m^2)
+ Scalar injectionDuration_; //! Duration of the injection in seconds
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/tutorial/solution/ex1/injection2pniproblem.hh b/tutorial/solution/ex1/injection2pniproblem.hh
index 7bba02cae53aa36be4c24d01b71f50c29bb459e9..a023716ea2ea6408e6b627722d4148eb12352f05 100644
--- a/tutorial/solution/ex1/injection2pniproblem.hh
+++ b/tutorial/solution/ex1/injection2pniproblem.hh
@@ -130,8 +130,6 @@ public:
InjectionProblem2PNI(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
- maxDepth_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.MaxDepth);
-
// initialize the tables of the fluid system
FluidSystem::init(/*tempMin=*/273.15,
/*tempMax=*/423.15,
@@ -139,6 +137,14 @@ public:
/*pMin=*/0.0,
/*pMax=*/30e6,
/*numP=*/300);
+
+ // name of the problem and output file
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+ // depth of the aquifer, units: m
+ aquiferDepth_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, AquiferDepth);
+ // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2)
+ // the duration of the injection, units: second
+ injectionDuration_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionDuration);
}
/*!
@@ -205,46 +211,39 @@ public:
*/
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- Scalar densityW = 1000.0;
- values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
+ // assume a constant density
+ const Scalar densityW = 1000;
+
+ // assume an intially hydrostatic liquid pressure profile
+ // note: we subtract rho_w*g*h because g is defined negative
+ const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
+
+ values[Indices::pressureIdx] = pw;
values[saturationIdx] = 0.0;
- values[temperatureIdx] = 283.0 + (maxDepth_ - globalPos[1])*0.03;
+ values[temperatureIdx] = 283.0 + (aquiferDepth_ - globalPos[1])*0.03;
}
/*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
+ * \brief Evaluates 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 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
- *
- * The \a values store the mass flux of each phase normal to the boundary.
- * Negative values indicate an inflow.
+ * \param globalPos The globalPosition of the boundary interface
*/
- void neumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &intersection,
- int scvIdx,
- int boundaryFaceIdx) const
+ void neumannAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
{
+ // initialize values to zero, i.e. no-flow Neumann boundary conditions
+ values = 0;
values = 0;
- GlobalPosition globalPos;
- if (isBox)
- globalPos = element.geometry().corner(scvIdx);
- else
- globalPos = intersection.geometry().center();
-
- if (globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_) {
+ if (this->timeManager().time() + this->timeManager().timeStepSize() < injectionDuration_
+ && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->bBoxMax()[0])
+ {
// inject air. negative values mean injection
- values[contiNEqIdx] = -1e-3; // kg/(s*m^2)
+ values[Indices::contiNEqIdx] = -1e-4; // kg/(s*m^2)
+ values[Indices::contiWEqIdx] = 0.0;
}
}
@@ -265,20 +264,27 @@ public:
*/
void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- Scalar densityW = 1000.0;
- values[pressureIdx] = 1e5 + (maxDepth_ - globalPos[1])*densityW*9.81;
+ // assume a constant density
+ const Scalar densityW = 1000;
+
+ // assume an intially hydrostatic liquid pressure profile
+ // note: we subtract rho_w*g*h because g is defined negative
+ const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
+
+ values[Indices::pressureIdx] = pw;
values[saturationIdx] = 0.0;
- values[temperatureIdx] = 283.0 + (maxDepth_ - globalPos[1])*0.03;
+ values[temperatureIdx] = 283.0 + (aquiferDepth_ - globalPos[1])*0.03;
if (globalPos[0] > 20 - eps_ && globalPos[0] < 30 + eps_ && globalPos[1] > 5 - eps_ && globalPos[1] < 35 + eps_)
values[temperatureIdx] = 380;
}
private:
- Scalar maxDepth_;
- static constexpr Scalar eps_ = 1.5e-7;
- std::string name_;
+ static constexpr Scalar eps_ = 1e-6;
+ std::string name_; //! Problem name
+ Scalar aquiferDepth_; //! Depth of the aquifer in m
+ Scalar injectionDuration_; //! Duration of the injection in seconds
};
} //end namespace