Exercise Fluidsystem

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 (exercise-fluidsystem a) as well as its mixture with water (exercise-fluidsystem b).

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

Preparing the exercise

• Navigate to the directory dumux-course/exercises/exercise-fluidsystem

1. Getting familiar with the code

Locate all the files you will need for this exercise

• The shared main file : exercise-fluidsystem.cc
• The input file for part a: exercise-fluidsystem_a.input
• The problem file for part a: 2pproblem.hh
• The input file for part b: exercise-fluidsystem_b.input
• 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.

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 the TwoP TypeTag (immiscible two-phase model properties) and the BoxModel TypeTag (specifies properties of the discretization scheme).

// Create new type tags
namespace TTag {
struct ExerciseFluidsystemTwoPTypeTag { using InheritsFrom = std::tuple<BoxModel, TwoP>; };
} // end namespace TTag

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:

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

while the BoxModel TypeTag can be found in the box/properties.hh header:

// The box discretization
#include <dumux/discretization/box.hh>

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

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:

// The water component
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/h2o.hh>

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.

// 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:

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

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

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):

// we use the immiscible fluid system here
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::ExerciseFluidsystemTwoPTypeTag>
{
private:
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using TabulatedH2O = Components::TabulatedComponent<Components::H2O<Scalar>>;
    using LiquidWaterPhase = typename FluidSystems::OnePLiquid<Scalar, TabulatedH2O>;
    /*!
     * Uncomment first line and comment second line for using the incompressible component
     * Uncomment second line and comment first line for using the compressible component
     */
    using LiquidMyComponentPhase = typename FluidSystems::OnePLiquid<Scalar, MyIncompressibleComponent<Scalar> >;
    // using LiquidMyComponentPhase = typename FluidSystems::OnePLiquid<Scalar, MyCompressibleComponent<Scalar> >;

public:
    using type = typename FluidSystems::TwoPImmiscible<Scalar, LiquidWaterPhase, LiquidMyComponentPhase>;
};

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.

/*!
 * \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:

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

cd build-cmake/exercises/exercise-fluidsystem
make exercise-fluidsystem_a
./exercise-fluidsystem_a exercise-fluidsystem_a.input

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:

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

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:

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:

3. Implement a new fluid system

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 two-component model properties:

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

// Create a new type tag for the problem
struct ExerciseFluidsystemTwoPTwoCTypeTag { using InheritsFrom = std::tuple<BoxModel, TwoPTwoC>; };
} // end namespace TTag

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.

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

// The fluid system property
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::ExerciseFluidsystemTwoPTwoCTypeTag>
{
private:
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
    using type = FluidSystems::H2OMyCompressibleComponent<Scalar>;
};

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

// the ficitious component that was created in exercise-fluidsystem a
#include <exercises/exercise-fluidsystem/components/mycompressiblecomponent.hh>

// the binary coefficients corresponding to this fluid system
#include <exercises/exercise-fluidsystem/binarycoefficients/h2omycompressiblecomponent.hh>

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

cd build-cmake/exercises/exercise-fluidsystem
make exercise-fluidsystem_b
./exercise-fluidsystem_b exercise-fluidsystem_b.input

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.

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:

/*!
 * \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:

cd build-cmake/exercises/exercise-fluidsystem
make exercise-fluidsystem_b
./exercise-fluidsystem_b exercise-fluidsystem_b.input