README.md 12.8 KB
Newer Older
1
# Exercise Fluidsystem
2

3
The aim of this exercise is to get familiar with the _DuMuX_ way of implementing new components (fluids) and fluid systems (mixtures). In the scope of this exercise, a new fictitious component is implemented
4
(exercise-fluidsystem a) as well as its mixture with water (exercise-fluidsystem b).
5
6
7
8
9

## Problem set-up

The domain has a size of 60 x 60 m and contains two low-permeable lenses. Initially, the domain is fully water saturated and the fictitious component is injected through the middle portion of the upper boundary 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 two lateral sides Dirichlet boundary conditions are applied (hydrostatic conditions for the pressure and zero saturation).

10
![](../extradoc/exercise-fluidsystem_setup.png)
11
12
13
14


## Preparing the exercise

15
* Navigate to the directory `dumux-course/exercises/exercise-fluidsystem`
16
17
18
19

### 1. Getting familiar with the code

Locate all the files you will need for this exercise
20
21
* The shared __main file__ : `exercise-fluidsystem.cc`
* The __input file__ for part a: `exercise-fluidsystem_a.input`
22
* The __problem file__ for part a: `2pproblem.hh`
23
* The __input file__ for part b: `exercise-fluidsystem_b.input`
24
25
26
27
28
29
30
31
32
33
34
35
36
37
* The __problem file__ for part b: `2p2cproblem.hh`
* The __spatial parameters file__: `spatialparams.hh`

Furthermore you will find the following folders:
* `binarycoefficients`: Stores headers containing data/methods on binary mixtures
* `components`: Stores headers containing data/methods on pure components
* `fluidsystems`: Stores headers containing data/methods on mixtures of pure components. Uses methods from `binarycoefficients`.

To see more components, fluidsystems and binarycoefficients implementations, have a look at the folder `dumux/material`.

### 2. Implement a new component

In the following, the basic steps required to set the desired fluid system are outlined. Here, this is done in the __problem file__, i.e. for this part of the exercise the code shown below is taken from the `2pproblem.hh` file.

38
In this part of the exercise we will consider a system consisting of two immiscible phases. Therefore, the _TypeTag_ for this problem (`ExerciseFluidsystemTwoPTypeTag`) derives from
39
the `TwoP` _TypeTag_ (immiscible two-phase model properties) and the `BoxModel` _TypeTag_ (specifies properties of the discretization scheme).
40
41

```c++
42
43
44
45
// Create new type tags
namespace TTag {
struct ExerciseFluidsystemTwoPTypeTag { using InheritsFrom = std::tuple<BoxModel, TwoP>; };
} // end namespace TTag
46
47
```

48
49
In order to be able to derive from these _TypeTags_, the declarations of the `TwoP` _TypeTag_ and `BoxModel` _TypeTag_ have to be included.
The `TwoP` _TypeTag_ can be found in the `2p/model.hh` header:
50
51
52
53
54
55

```c++
// The numerical model
#include <dumux/porousmediumflow/2p/model.hh>
```

56
while the `BoxModel` _TypeTag_ can be found in the `box/properties.hh` header:
57
58

```c++
59
// The box discretization
60
61
62
#include <dumux/discretization/box/properties.hh>
```

63
64
For a cell-centered scheme, you could derive from `CCTpfaModel` or `CCMpfaModel` instead (and, of course, include the right headers).

65
As one of the two phases we want to use water and we want to precompute tables on which the properties are then interpolated in order to save computational time. Thus, in a first step we have to include the following headers:
66
67
68
69
70
71

```c++
// The water component
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/h2o.hh>
```
72
73
The other phase will be created from our new component, where we want to implement an incompressible and a compressible variant.
The respective headers are prepared, but still incomplete. The compressible variant is still commented so that compilation does not fail when finishing the incompressible variant.
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

```c++
// The components that will be created in this exercise
#include "components/myincompressiblecomponent.hh"
// #include "components/mycompressiblecomponent.hh"
```
As mentioned above, we want to simulate two non-mixing components. The respective fluid system is found in:

```c++
// The two-phase immiscible fluid system
#include <dumux/material/fluidsystems/2pimmiscible.hh>
```

This fluid system expects __phases__ as input and so far we have only included the components, which contain data on the pure component for all physical states. Thus, we need to include

```c++
// We will only have liquid phases here
#include <dumux/material/fluidsystems/1pliquid.hh>
```

94
which creates a _liquid phase_ from a given component. Finally, using all of the included classes we set the fluid system property by choosing that
95
the water phase is liquid (OnePLiquid) and consists of the tabulated water component, and
96
the other phase is liquid as well and consists of the incompressible fictitious component. Both will make up the immiscible fluid system (TwoPImmiscible):
97
98
99
100


```c++
// we use the immiscible fluid system here
101
102
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::ExerciseFluidsystemTwoPTypeTag>
103
104
{
private:
105
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
106
    using TabulatedH2O = Components::TabulatedComponent<Components::H2O<Scalar>>;
107
    using LiquidWaterPhase = typename FluidSystems::OnePLiquid<Scalar, TabulatedH2O>;
108
109
110
111
    /*!
     * Uncomment first line and comment second line for using the incompressible component
     * Uncomment second line and comment first line for using the compressible component
     */
112
    using LiquidMyComponentPhase = typename FluidSystems::OnePLiquid<Scalar, MyIncompressibleComponent<Scalar> >;
113
    // using LiquidMyComponentPhase = typename FluidSystems::OnePLiquid<Scalar, MyCompressibleComponent<Scalar> >;
114
115

public:
116
    using type = typename FluidSystems::TwoPImmiscible<Scalar, LiquidWaterPhase, LiquidMyComponentPhase>;
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
};
```

### 2.1. Incompressible component

Open the file `myincompressiblecomponent.hh`. You can see in line 42 that a component should always derive from the _Base_ class (see `dumux/material/components/base.hh`), which defines the interface of a _DuMuX_ component with possibly required functions to be overloaded by the actual implementation. Additionally it is required for liquids to derive from the _Liquid_ class (see `dumux/material/components/liquid.hh`), for gases to derive from the _Gas_ class (see `dumux/material/components/gas.hh`) and for solids to derive from the _Solid_ class (see `dumux/material/components/solid.hh`), with functions specific to liquid, gas or solid.

```c++
/*!
 * \ingroup Components
 * \brief A ficitious component to be implemented in exercise 3.
 *
 * \tparam Scalar The type used for scalar values
 */
template <class Scalar>
class MyIncompressibleComponent
: public Components::Base<Scalar, MyIncompressibleComponent<Scalar> >
, public Components::Liquid<Scalar, MyIncompressibleComponent<Scalar> >
```

__Task__:

Implement an incompressible component into the file `myincompressiblecomponent.hh`, which has the following specifications:

| Parameter | unit | value |
| -----| --------| -------- |
143
144
| $`M`$ | $`kg/mol`$   | $`131.39 \cdot 10^{-3}`$ |
| $`\rho_{liquid}`$ | $`kg/m^3`$   | $`1460`$   |
145
146
147
148
149
150
151
| $`\mu_{liquid}`$ | $`Pa \cdot s`$   | $`5.7 \cdot 10^{-4}`$   |

In order to do so, have a look at the files `dumux/material/components/base.hh` and `dumux/material/components/liquid.hh` to see how the interfaces are defined and overload them accordingly.

In order to execute the program, change to the build directory and compile and execute the program by typing

```bash
152
cd build-cmake/exercises/exercise-fluidsystem
153
154
make exercise-fluidsystem_a
./exercise-fluidsystem_a exercise-fluidsystem_a.input
155
156
```

157
The saturation distribution of the nonwetting phase S$`_n`$ (the phase consisting of our fictitious incompressible component) at the final simulation time should look like this:
158

159
![](../extradoc/exercise-fluidsystem_a_solution.png)
160
161
162
163
164
165
166

### 2.2. Compressible component

We now want to implement a pressure-dependent density for our component. Open the file `mycompressiblecomponent.hh` and copy in the functions you implemented for the incompressible variant. Now substitute the method that returns the density by the following expression:

$`\displaystyle \rho_{MyComp} = \rho_{min} + \frac{ \rho_{max} - \rho_{min} }{ 1 + \rho_{min}*e^{-1.0*k*(\rho_{max} - \rho_{min})*p} } `$

167
where $`p`$ is the pressure and $`\rho_{min} = 1440 `$, $`\rho_{max} = 1480 `$ and $`k = 5 \cdot 10^{-7} `$. Also, make sure the header is included in the `2pproblem.hh` file by uncommenting line 54. Furthermore, the new component has to be set as a liquid phase in the fluid system, i.e. comment line 109 and uncomment line 110. The density distribution of this phase (rhoN) at the final simulation time should look like this:
168

169
![](../extradoc/exercise-fluidsystem_a_solution2.png)
170

171
172
173
You can plot the density of the phase consisting of your compressible component by setting `PlotDensity` in `exercise-fluidsystem_a.input` to `true` and starting the simulation again.
Compare the gnuplot output to the following plot of the density function from above:

174
![](../extradoc/exercise-fluidsystem_a_densityfunction.png)
175

176
177
### 3. Implement a new fluid system

178
The problem file for this part of the exercise is `2p2cproblem.hh`. We now want to implement a new fluid system consisting of two liquid phases, which are water and the previously implemented compressible component. We will consider compositional effects, which is why we now have to derive our _TypeTag_ (`ExerciseFluidsystemTwoPTwoCTypeTag`) from a _TypeTag_ (`TwoPTwoC`) that holds the miscible two-phase
179
two-component model properties:
180
181
182
183
184
185
186
187

```c++
// The numerical model
#include <dumux/porousmediumflow/2p2c/model.hh>
```

```c++
// Create a new type tag for the problem
188
189
struct ExerciseFluidsystemTwoPTwoCTypeTag { using InheritsFrom = std::tuple<BoxModel, TwoPTwoC>; };
} // end namespace TTag
190
191
192
193
194
195
196
197
198
199
200
```

The new fluid system is to be implemented in the file `fluidsystems/h2omycompressiblecomponent.hh`. This is already included in the problem and the fluid system property is set accordingly.

```c++
// The fluid system that is created in this exercise
#include "fluidsystems/h2omycompressiblecomponent.hh"
```

```c++
// The fluid system property
201
202
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::ExerciseFluidsystemTwoPTwoCTypeTag>
203
204
{
private:
205
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
206
207
208
209
210
public:
    using type = FluidSystems::H2OMyCompressibleComponent<Scalar>;
};
```

211
In the `fluidsystems/h2omycompressiblecomponent.hh` file, your implemented compressible component and the binary coefficient files are already included.
212
213

```c++
214
// the ficitious component that was created in exercise-fluidsystem a
215
#include <exercises/exercise-fluidsystem/components/mycompressiblecomponent.hh>
216
217

// the binary coefficients corresponding to this fluid system
218
#include <exercises/exercise-fluidsystem/binarycoefficients/h2omycompressiblecomponent.hh>
219
220
221
222
223
224
225
226
227
228
229
```

__Task__:

Under the assumption that one molecule of `MyCompressibleComponent` displaces exactly one molecule of water, the water phase density can be expressed as follows:

$` \rho_{w} = \frac{ \rho_{w, pure} }{ M_{H_2O} }*(M_{H_2O}*x_{H_2O} + M_{MyComponent}*x_{MyComponent}) `$

Implement this dependency in the `density()` method in the fluid system. In order to compile and execute the program run

```bash
230
cd build-cmake/exercises/exercise-fluidsystem
231
232
make exercise-fluidsystem_b
./exercise-fluidsystem_b exercise-fluidsystem_b.input
233
234
```

235
You will observe an error message and an abortion of the program. This is due to the fact that in order for the constraint solver and other mechanisms in the two-phase two-component model to work, an additional functionality in the component has to be implemented: the model has to know the vapour pressure. As in the previous exercise, check the `dumux/material/components/base.hh` file for this function and implement it into `mycompressiblecomponent.hh`. For the vapour pressure, use a value of $`3900`$  Pa.
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260

### 4. Change wettability of the porous medium

In the `spatialparams.hh` file, we can find the following function, with which we can specify which phase of the fluid system is to be considered as the wetting phase at a given position within the domain:

```cpp
/*!
 * \brief Function for defining which phase is to be considered as the wetting phase.
 *
 * \return the wetting phase index
 * \param globalPos The position of the center of the element
 */
template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{
    // Our fluid system is H2OMyCompressibleComponent
    // We want to define water as the wetting phase in
    // the entire domain (see fluid system for the phase indices)
    return FluidSystem::phase0Idx;
}
```

Change this function such that the the phase of our new component is the wetting phase __only__ within the lenses. Execute the program of task 3 again:

```bash
261
cd build-cmake/exercises/exercise-fluidsystem
262
263
make exercise-fluidsystem_b
./exercise-fluidsystem_b exercise-fluidsystem_b.input
264
```