Beatrix Becker committed Jul 13, 2018 1 # Exercise Fluidsystem  Timo Koch committed Jul 09, 2018 2   Dennis Gläser committed Jul 17, 2018 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  Beatrix Becker committed Jul 13, 2018 4 (exercise-fluidsystem a) as well as its mixture with water (exercise-fluidsystem b).  Timo Koch committed Jul 09, 2018 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).  Beatrix Becker committed Jul 13, 2018 10 ![](../extradoc/exercise-fluidsystem_setup.png)  Timo Koch committed Jul 09, 2018 11 12 13 14  ## Preparing the exercise  Beatrix Becker committed Jul 12, 2018 15 * Navigate to the directory dumux-course/exercises/exercise-fluidsystem  Timo Koch committed Jul 09, 2018 16 17 18 19  ### 1. Getting familiar with the code Locate all the files you will need for this exercise  Beatrix Becker committed Jul 13, 2018 20 21 * The shared __main file__ : exercise-fluidsystem.cc * The __input file__ for part a: exercise-fluidsystem_a.input  Timo Koch committed Jul 09, 2018 22 * The __problem file__ for part a: 2pproblem.hh  Beatrix Becker committed Jul 13, 2018 23 * The __input file__ for part b: exercise-fluidsystem_b.input  Timo Koch committed Jul 09, 2018 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.  Beatrix Becker committed Jul 13, 2018 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  Beatrix Becker committed Jul 12, 2018 39 the TwoP _TypeTag_ (immiscible two-phase model properties) and the BoxModel _TypeTag_ (specifies properties of the discretization scheme).  Timo Koch committed Jul 09, 2018 40 41  c++  Felix Weinhardt committed Nov 29, 2018 42 43 44 45 // Create new type tags namespace TTag { struct ExerciseFluidsystemTwoPTypeTag { using InheritsFrom = std::tuple; }; } // end namespace TTag  Timo Koch committed Jul 09, 2018 46 47   Beatrix Becker committed Jul 12, 2018 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:  Timo Koch committed Jul 09, 2018 50 51 52 53 54 55  c++ // The numerical model #include   Beatrix Becker committed Jul 12, 2018 56 while the BoxModel _TypeTag_ can be found in the box/properties.hh header:  Timo Koch committed Jul 09, 2018 57 58  c++  Felix Weinhardt committed Nov 30, 2018 59 // The box discretization  Timo Koch committed Jul 09, 2018 60 61 62 #include   Beatrix Becker committed Jul 12, 2018 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:  Timo Koch committed Jul 09, 2018 66 67 68 69 70 71  c++ // The water component #include #include   Beatrix Becker committed Jul 12, 2018 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.  Timo Koch committed Jul 09, 2018 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  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   Beatrix Becker committed Jul 12, 2018 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  Beatrix Becker committed Jul 12, 2018 96 the other phase is liquid as well and consists of the incompressible fictitious component. Both will make up the immiscible fluid system (TwoPImmiscible):  Timo Koch committed Jul 09, 2018 97 98 99 100  c++ // we use the immiscible fluid system here  Felix Weinhardt committed Nov 29, 2018 101 102 template struct FluidSystem  Timo Koch committed Jul 09, 2018 103 104 { private:  Felix Weinhardt committed Nov 29, 2018 105  using Scalar = GetPropType;  Timo Koch committed Jul 09, 2018 106  using TabulatedH2O = Components::TabulatedComponent>;  Beatrix Becker committed Jul 12, 2018 107  using LiquidWaterPhase = typename FluidSystems::OnePLiquid;  Timo Koch committed Jul 09, 2018 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 */  Beatrix Becker committed Jul 12, 2018 112  using LiquidMyComponentPhase = typename FluidSystems::OnePLiquid >;  Ned Coltman committed Jul 17, 2018 113  // using LiquidMyComponentPhase = typename FluidSystems::OnePLiquid >;  Timo Koch committed Jul 09, 2018 114 115  public:  Beatrix Becker committed Jul 12, 2018 116  using type = typename FluidSystems::TwoPImmiscible;  Timo Koch committed Jul 09, 2018 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 MyIncompressibleComponent : public Components::Base > , public Components::Liquid >  __Task__: Implement an incompressible component into the file myincompressiblecomponent.hh, which has the following specifications: | Parameter | unit | value | | -----| --------| -------- |  Beatrix Becker committed Jul 13, 2018 143 144 | $M$ | $kg/mol$ | $131.39 \cdot 10^{-3}$ | | $\rho_{liquid}$ | $kg/m^3$ | $1460$ |  Timo Koch committed Jul 09, 2018 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  Beatrix Becker committed Jul 13, 2018 153 154 make exercise-fluidsystem_a ./exercise-fluidsystem_a exercise-fluidsystem_a.input  Timo Koch committed Jul 09, 2018 155 156   Beatrix Becker committed Jul 12, 2018 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:  Timo Koch committed Jul 09, 2018 158   Beatrix Becker committed Jul 13, 2018 159 ![](../extradoc/exercise-fluidsystem_a_solution.png)  Timo Koch committed Jul 09, 2018 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} } $  Felix Weinhardt committed Nov 30, 2018 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:  Timo Koch committed Jul 09, 2018 168   Beatrix Becker committed Jul 13, 2018 169 ![](../extradoc/exercise-fluidsystem_a_solution2.png)  Timo Koch committed Jul 09, 2018 170   Beatrix Becker committed Jul 17, 2018 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:  Dennis Gläser committed Jul 17, 2018 174 ![](../extradoc/exercise-fluidsystem_a_densityfunction.png)  Beatrix Becker committed Jul 17, 2018 175   Timo Koch committed Jul 09, 2018 176 177 ### 3. Implement a new fluid system  Dennis Gläser committed Jul 17, 2018 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:  Timo Koch committed Jul 09, 2018 180 181 182 183 184 185 186 187  c++ // The numerical model #include  c++ // Create a new type tag for the problem  Felix Weinhardt committed Nov 29, 2018 188 189 struct ExerciseFluidsystemTwoPTwoCTypeTag { using InheritsFrom = std::tuple; }; } // end namespace TTag  Timo Koch committed Jul 09, 2018 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  Felix Weinhardt committed Nov 29, 2018 201 202 template struct FluidSystem  Timo Koch committed Jul 09, 2018 203 204 { private:  Felix Weinhardt committed Nov 29, 2018 205  using Scalar = GetPropType;  Timo Koch committed Jul 09, 2018 206 207 208 209 210 public: using type = FluidSystems::H2OMyCompressibleComponent; };   211 In the fluidsystems/h2omycompressiblecomponent.hh file, your implemented compressible component and the binary coefficient files are already included.  Timo Koch committed Jul 09, 2018 212 213  c++  Beatrix Becker committed Jul 13, 2018 214 // the ficitious component that was created in exercise-fluidsystem a  Beatrix Becker committed Jul 13, 2018 215 #include  Timo Koch committed Jul 09, 2018 216 217  // the binary coefficients corresponding to this fluid system  Beatrix Becker committed Jul 13, 2018 218 #include  Timo Koch committed Jul 09, 2018 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  Beatrix Becker committed Jul 13, 2018 231 232 make exercise-fluidsystem_b ./exercise-fluidsystem_b exercise-fluidsystem_b.input  Timo Koch committed Jul 09, 2018 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.  Dennis Gläser committed Jul 10, 2018 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 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  Beatrix Becker committed Jul 13, 2018 262 263 make exercise-fluidsystem_b ./exercise-fluidsystem_b exercise-fluidsystem_b.input  Dennis Gläser committed Jul 10, 2018 264