README.md 13.1 KB
Newer Older
1
# Exercise Fluidsystem (DuMuX Course)
2

3
The aim of this exercise is to get familiar with the _DuMu<sup>x</sup>_ way of implementing new components (fluids) and fluid systems (mixtures). In the scope of this exercise, a new fictitious component is implemented (Section 2) as well as its mixture with water (Section 3).
4
5
6
7
8

## 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).

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


## Preparing the exercise

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

### 1. Getting familiar with the code

Locate all the files you will need for this exercise
19
20
* The shared __main file__ : `main.cc`
* The __input file__ for part a: `aparams.input`
21
* The __problem file__ for part a: `2pproblem.hh`
22
* The __properties file__ for part a: `2pproperties.hh`
23
* The __input file__ for part b: `bparams.input`
24
* The __problem file__ for part b: `2p2cproblem.hh`
25
* The __properties file__ for part b: `2p2cproperties.hh`
26
27
28
29
30
31
32
33
34
35
36
* 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

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

Felix Weinhardt's avatar
Felix Weinhardt committed
39
In this part of the exercise we will consider a system consisting of two immiscible phases. Therefore, the _TypeTag_ for this problem (`ExerciseFluidsystemTwoP`) derives from
40
the `TwoP` _TypeTag_ (immiscible two-phase model properties) and the `BoxModel` _TypeTag_ (specifies properties of the discretization scheme).
41
42

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

49
50
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:
51
52
53
54
55
56

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

Katharina Heck's avatar
Katharina Heck committed
57
while the `BoxModel` _TypeTag_ can be found in the `discretization/box.hh` header:
58
59

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

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

66
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. Since we need this part in both problem file and properties file, we have to include the following two headers in our problem file, i.e. `2pproblem.hh` file, and the properties file has access to them through problem file.
67
68
69
70
71
72

```c++
// The water component
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/h2o.hh>
```
73
74
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.
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

```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>
```

95
96
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 the water phase is liquid (`OnePLiquid`) and consists of the tabulated water component, and
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
template<class TypeTag>
Felix Weinhardt's avatar
Felix Weinhardt committed
102
struct FluidSystem<TypeTag, TTag::ExerciseFluidsystemTwoP>
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

    // TODO: dumux-course-task 2.2
    //****** Select the corresponding component for the task *****//
    // Uncomment first line and comment second 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
};
```

### 2.1. Incompressible component

122
Open the file `myincompressiblecomponent.hh`.  A component should always derive from the _Base_ class - thus we include `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.
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

```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
make exercise_fluidsystem_a
154
./exercise_fluidsystem_a aparams.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 `2pproperties.hh` file by uncommenting the corresponding line. Furthermore, the new component has to be set as a liquid phase in the fluid system. To do so, search for `TODO: dumux-course-task 2.2`. Comment out the corresponding line and uncomment the other. 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
You can plot the density of the phase consisting of your compressible component by setting `PlotDensity` in `aparams.input` to `true` and starting the simulation again.
172
173
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
179
180
The problem file and properties file for this part of the exercise are `2p2cproblem.hh` and `2p2cproperties.hh`, respectively. 
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_ (`ExerciseFluidsystemTwoPTwoC`) from a _TypeTag_ (`TwoPTwoC`) that holds the miscible two-phase
    two-component model properties:
181
182
183
184
185
186
187
188

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

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

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
202
template<class TypeTag>
Felix Weinhardt's avatar
Felix Weinhardt committed
203
struct FluidSystem<TypeTag, TTag::ExerciseFluidsystemTwoPTwoC>
204
205
{
private:
206
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
207
208
209
210
211
public:
    using type = FluidSystems::H2OMyCompressibleComponent<Scalar>;
};
```

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

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

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

__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}) `$

228
Implement this dependency in the `density()` method in the fluid system. In order to compile and execute the program, run the following commands:
229
230

```bash
231
cd build-cmake/exercises/exercise-fluidsystem
232
make exercise_fluidsystem_b
233
./exercise_fluidsystem_b bparams.input
234
235
```

236
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.
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258

### 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;
}
```

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

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