tutorial-coupled.tex 30.9 KB
 Karin Erbertseder committed Dec 13, 2010 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 \section{Box method - A short introduction}\label{box} For the spatial discretization the so called BOX-method is used which unites the advantages of the finite-volume (FV) and finite-element (FE) methods. First, the model domain $G$ is discretized with a FE mesh consisting of nodes i and corresponding elements $E_k$. Then, a secondary FV mesh is constructed by connecting the midpoints and barycenters of the elements surrounding node i creating a box $B_i$ around node i (see Figure \ref{pc:box}a). \begin{figure} [h] \includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/box_disc} \caption{\label{pc:box} Discretization of the BOX-method} \end{figure} The FE mesh divides the box $B_i$ into subcontrolvolumes (scv's) $b^k_i$ (see Figure \ref{pc:box}b). Figure \ref{pc:box}c shows the finite element $E_k$ and the scv's $b^k_i$ inside $E_k$, which belong to four different boxes $B_i$. Also necessary for the discretization are the faces of the subcontrolvolumes (scvf's) $e^k_{ij}$ between the scv's $b^k_i$ and $b^k_j$, where $|e^k_{ij}|$ is the length of the scvf. The integration points $x^k_{ij}$ on $e^k_{ij}$ and the outer normal vector $\mathbf n^k_{ij}$ are also to be defined (see Figure \ref{pc:box}c). The advantage of the FE method is that unstructured grids can be used, while the FV method is mass conservative. The idea is to apply the FV method (balance of fluxes across the interfaces) to each FV box $B_i$ and to get the fluxes across the interfaces $e^k_{ij}$ at the integration points $x^k_{ij}$ from the FE approach. Consequently, at each scvf the following expression results: f(\tilde u(x^k_{ij})) \cdot \mathbf n^k_{ij} \: |e^k_{ij}| \qquad \textrm{with} \qquad \tilde u(x^k_{ij}) = \sum_i N_i(x^k_{ij}) \cdot \hat u_i . In the following the discretization of the balance equation is going to be derived. From the \textsc{reynolds}s transport theorem follows the general balance equation: \underbrace{\int_G \frac{\partial}{\partial t} \: u \: dG}_{1} + \underbrace{\int_{\partial G} (\mathbf{v} u + \mathbf w) \cdot \textbf n \: d\varGamma}_{2} = \underbrace{\int_G q \: dG}_{3} f(u) = \int_G \frac{\partial u}{\partial t} \: dG + \int_{G} \nabla \cdot \underbrace{\left[ \mathbf{v} u + \mathbf w(u)\right] }_{F(u)} \: dG - \int_G q \: dG = 0 where term 1 describes the changes of entity $u$ within a control volume over time, term 2 the advective, diffusive and dispersive fluxes over the interfaces of the control volume and term 3 is the source and sink term. $G$ denotes the model domain and $F(u) = F(\mathbf v, p) = F(\mathbf v(x,t), p(x,t))$.\\ Like the FE method, the BOX-method follows the principle of weighted residuals. In the function $f(u)$ the unknown $u$ is approximated by discrete values at the nodes of the FE mesh $\hat u_i$ and linear basis functions $N_i$ yielding an approximate function $f(\tilde u)$. For $u\in \lbrace \mathbf v, p, x^\kappa \rbrace$ this means \begin{minipage}[b]{0.5\textwidth} \label{eq:p} \tilde p = \sum_i N_i \hat{p_i} \label{eq:v} \tilde{\mathbf v} = \sum_i N_i \hat{\mathbf v} \label{eq:x} \tilde x^\kappa = \sum_i N_i \hat x^\kappa \end{minipage} \hfill \begin{minipage}[b]{0.5\textwidth} \label{eq:dp} \nabla \tilde p = \sum_i \nabla N_i \hat{p_i} \label{eq:dv} \nabla \tilde{\mathbf v} = \sum_i \nabla N_i \hat{\mathbf v} \label{eq:dx} \nabla \tilde x^\kappa = \sum_i \nabla N_i \hat x^\kappa . \end{minipage} Due to the approximation with node values and basis functions the differential equations are not exactly fulfilled anymore but a residual $\varepsilon$ is produced. f(u) = 0 \qquad \Rightarrow \qquad f(\tilde u) = \varepsilon Application of the principle of weighted residuals, meaning the multiplication of the residual $\varepsilon$ with a weighting function $W_j$ and claiming that this product has to vanish within the whole domain, \int_G W_j \cdot \varepsilon \: \overset {!}{=} \: 0 \qquad \textrm{with} \qquad \sum_j W_j =1 yields the following equation: \int_G W_j \frac{\partial \tilde u}{\partial t} \: dG + \int_G W_j \cdot \left[ \nabla \cdot F(\tilde u) \right] \: dG - \int_G W_j \cdot q \: dG = \int_G W_j \cdot \varepsilon \: dG \: \overset {!}{=} \: 0 . Then, the chain rule and the \textsc{green-gaussian} integral theorem are applied. \int_G W_j \frac{\partial \sum_i N_i \hat u_i}{\partial t} \: dG + \int_{\partial G} \left[ W_j \cdot F(\tilde u)\right] \cdot \mathbf n \: d\varGamma_G + \int_G \nabla W_j \cdot F(\tilde u) \: dG - \int_G W_j \cdot q \: dG = 0 A mass lumping technique is applied by assuming that the storage capacity is reduced to the nodes. This means that the integrals $M_{i,j} = \int_G W_j \: N_i \: dG$ are replaced by the mass lumping term $M^{lump}_{i,j}$ which is defined as: M^{lump}_{i,j} =\begin{cases} \int_G W_j \: dG = \int_G N_i \: dG = V_i &i = j\\ 0 &i \neq j\\ \end{cases} where $V_i$ is the volume of the FV box $B_i$ associated with node i. The application of this assumption in combination with $\int_G W_j \:q \: dG = V_i \: q$ yields V_i \frac{\partial \hat u_i}{\partial t} + \int_{\partial G} \left[ W_j \cdot F(\tilde u)\right] \cdot \mathbf n \: d\varGamma_G + \int_G \nabla W_j \cdot F(\tilde u) \: dG- V_i \cdot q = 0 Defining the weighting function $W_j$ to be piecewise constant over a control volume box $B_i$ W_j(x) = \begin{cases} 1 &x \in B_i \\ 0 &x \notin B_i\\ \end{cases} causes $\nabla W_j = 0$: \label{eq:disc1} V_i \frac{\partial \hat u_i}{\partial t} + \int_{\partial B_i} \left[ W_j \cdot F(\tilde u)\right] \cdot \mathbf n \; d{\varGamma}_{B_i} - V_i \cdot q = 0 . The consideration of the time discretization and inserting $W_j = 1$ finally leads to the discretized form which will be applied to the mathematical flow and transport equations: \label{eq:discfin} V_i \frac{\hat u_i^{n+1} - \hat u_i^{n}}{\Delta t} + \int_{\partial B_i} F(\tilde u^{n+1}) \cdot \mathbf n \; d{\varGamma}_{B_i} - V_i \: q^{n+1} \: = 0  Karin Erbertseder committed Dec 13, 2010 125 \section[Fully-coupled model]{Solving a problem using a Fully-Coupled Model}\label{tutorial-coupled}  Bernd Flemisch committed Jul 16, 2010 126 127 128 129  The process of solving a problem using \Dumux can be roughly divided into four parts: \begin{enumerate} \item The geometry of the problem and correspondingly a grid have to be defined.  Klaus Mosthaf committed Oct 19, 2010 130  \item Material properties and constitutive relationships have to be selected.  Bernd Flemisch committed Jul 16, 2010 131 132 133 134  \item Boundary conditions as well as initial conditions have to be defined. \item A suitable model has to be chosen. \end{enumerate}  135 The problem that is solved in this tutorial is illustrated in Figure \ref{tutorial-coupled:problemfigure}. A rectangular domain with no flow boundaries on the top and on the bottom, which is initially saturated with oil, is considered. Water infiltrates from the left side into the domain. Gravity effects are neglected.  Bernd Flemisch committed Jul 16, 2010 136 137 138 139 140 141 142 143 144 145 146 147  \begin{figure}[h] \psfrag{x}{x} \psfrag{y}{y} \psfrag{no flow}{no flow} \psfrag{water}{\textbf{water}} \psfrag{oil}{\textcolor{white}{\textbf{oil}}} \psfrag{p_w = 2 x 10^5 [Pa]}{$p_w = 2 \times 10^5$ [Pa]} \psfrag{p_w_initial = 2 x 10^5 [Pa]}{\textcolor{white}{\textbf{$\mathbf{p_{w_{initial}} = 2 \times 10^5}$ [Pa]}}} \psfrag{S_n = 0}{$S_n = 0$} \psfrag{S_n_initial = 0}{\textcolor{white}{$\mathbf{S_{n_{initial}} = 1}$}} \psfrag{q_w = 0 [kg/m^2s]}{$q_w = 0$ $\left[\frac{\textnormal{kg}}{\textnormal{m}^2 \textnormal{s}}\right]$}  Melanie Darcis committed Dec 07, 2010 148 \psfrag{q_n = -3 x 10^-4 [kg/m^2s]}{$q_n = -3 \times 10^{-2}$ $\left[\frac{\textnormal{kg}}{\textnormal{m}^2 \textnormal{s}}\right]$}  Bernd Flemisch committed Jul 16, 2010 149 150 151 152 153 154 155 156 157 158 \centering \includegraphics[width=0.9\linewidth,keepaspectratio]{EPS/tutorial-problemconfiguration} \caption{Geometry of the tutorial problem with initial and boundary conditions.}\label{tutorial-coupled:problemfigure} \end{figure} The equations that are solved here are the mass balances of oil and water: \begin{align} \label{massbalancewater} \frac {\partial (\phi \, S_{w}\, \varrho_{w})}{\partial t}  Philipp Nuske committed Dec 10, 2010 159  -  Melanie Darcis committed Nov 11, 2010 160  \nabla \cdot \left( \varrho_{w} \, \frac{k_{rw}}{\mu_{w}} \, \mathbf{K}\;\nabla p_w \right)  Bernd Flemisch committed Jul 16, 2010 161  -  Melanie Darcis committed Nov 11, 2010 162  q_w  Bernd Flemisch committed Jul 16, 2010 163 164 165 166  & = 0 \\ \label{massbalanceoil} \frac {\partial (\phi \, S_{o}\, \varrho_{o})}{\partial t}  Philipp Nuske committed Dec 10, 2010 167  -  Melanie Darcis committed Nov 11, 2010 168  \nabla \cdot \left( \varrho_{o} \, \frac{k_{ro}}{\mu_{o}} \, \mathbf{K}\;\nabla p_o \right)  Bernd Flemisch committed Jul 16, 2010 169  -  Melanie Darcis committed Nov 11, 2010 170  q_o  Bernd Flemisch committed Jul 16, 2010 171 172 173 174 175 176 177  & = 0 \end{align} \subsection{The main file} Listing \ref{tutorial-coupled:mainfile} shows the main file  Klaus Mosthaf committed Oct 19, 2010 178 \texttt{tutorial/tutorial\_coupled.cc} for the coupled two-phase  Bernd Flemisch committed Jul 16, 2010 179 180 181 182 183 model. This file needs to be executed to solve the problem described above. \begin{lst}[File tutorial/tutorial\_coupled.cc]\label{tutorial-coupled:mainfile} \mbox{} \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,  Andreas Lauser committed Feb 04, 2011 184  numberstyle=\tiny, numbersep=5pt, firstline=28]{../../tutorial/tutorial_coupled.cc}  Bernd Flemisch committed Jul 16, 2010 185 186 187 \end{lst} From line \ref{tutorial-coupled:include-begin} to line  Philipp Nuske committed Dec 10, 2010 188 \ref{tutorial-coupled:include-end} the \Dune and \Dumux files which  Bernd Flemisch committed Jul 16, 2010 189 190 191 192 193 194 195 196 contain the needed functions and classes are included. At line \ref{tutorial-coupled:set-type-tag} the type tag of the problem which is going to be simulated is set. All other data types can be retrieved by the \Dumux property system and only depend on this single type tag. Retrieving them is done between line \ref{tutorial-coupled:retrieve-types-begin} and \ref{tutorial-coupled:retrieve-types-end}. For an introduction to the  197 property system, see section \ref{sec:propertysystem}.  Bernd Flemisch committed Jul 16, 2010 198 199  The first thing which should be done at run time is to initialize the  Philipp Nuske committed Dec 10, 2010 200 message passing interface using \Dune's \texttt{MPIHelper} class. Line  Klaus Mosthaf committed Oct 19, 2010 201 \ref{tutorial-coupled:init-mpi} is essential if the simulation is  Bernd Flemisch committed Jul 16, 2010 202 203 204 intended to be run on more than one processor at the same time. Next, the command line arguments are parsed starting at line \ref{tutorial-coupled:parse-args-begin} until line  Klaus Mosthaf committed Oct 19, 2010 205 206 207 \ref{tutorial-coupled:parse-args-end}. In this example, it is checked if and from which time on a previous run of the simulation should be restarted. Furthermore, we parse the time when the simulation ends and the initial time step size.  Bernd Flemisch committed Jul 16, 2010 208   Klaus Mosthaf committed Oct 19, 2010 209 After this, a grid is created in line  Bernd Flemisch committed Jul 16, 2010 210 \ref{tutorial-coupled:create-grid} and the problem is instantiated for  Klaus Mosthaf committed Oct 19, 2010 211 its leaf grid view in line \ref{tutorial-coupled:instantiate-problem}.  Philipp Nuske committed Dec 10, 2010 212 Finally, on line \ref{tutorial-coupled:begin-restart} a state written to  Klaus Mosthaf committed Oct 19, 2010 213 214 disk by a previous simulation run is restored if requested by the user. The simulation procedure is started at line  Bernd Flemisch committed Jul 16, 2010 215 216 217 218 219 220 221 222 223 224 225 \ref{tutorial-coupled:execute}. \subsection{The problem class} When solving a problem using \Dumux, the most important file is the so-called \textit{problem file} as shown in listing \ref{tutorial-coupled:problemfile} of \texttt{tutorialproblem\_coupled.hh}. \begin{lst}[File tutorial/tutorialproblem\_coupled.hh]\label{tutorial-coupled:problemfile} \mbox{} \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,  Andreas Lauser committed Feb 04, 2011 226 numberstyle=\tiny, numbersep=5pt, firstline=27]{../../tutorial/tutorialproblem_coupled.hh}  Bernd Flemisch committed Jul 16, 2010 227 228 \end{lst}  Klaus Mosthaf committed Oct 18, 2010 229 First, a new type tag is created for the problem in line  Bernd Flemisch committed Jul 16, 2010 230 231 232 233 234 \ref{tutorial-coupled:create-type-tag}. In this case, the new type tag inherits all properties defined for the \texttt{BoxTwoP} type tag, which means that for this problem the two-phase box model is chosen as discretization scheme. On line \ref{tutorial-coupled:set-problem}, a problem class is attached to the new type tag, while the grid which  Klaus Mosthaf committed Oct 19, 2010 235 is going to be used is defined in line \ref{tutorial-coupled:set-grid} --  Philipp Nuske committed Dec 10, 2010 236 in this case that is \texttt{SGrid}. Since there's no uniform  Klaus Mosthaf committed Oct 19, 2010 237 mechanism to allocate grids in \Dune, the \texttt{Grid} property also contains  Benjamin Faigle committed Feb 04, 2011 238 239 240 241 a static \texttt{create()} method which provides just that. Therein, the three variables of Type \texttt{Dune::FieldVector} define the lower left corner of the domain (\texttt{L}), the upper right corner of the domain (\texttt{H}) and the number of cells in $x$ and $y$ direction (\texttt{N}). Next,  Bernd Flemisch committed Jul 16, 2010 242 the appropriate fluid system, that specifies both information about  Andreas Lauser committed Nov 10, 2010 243 244 245 246 247 248 249 250 251 252 the fluid mixture as well as about the pure substances, has to be chosen. The two-phase model defaults to the \texttt{FluidSystem2P} which assumes immiscibility of the phases, but requires that the wetting and non-wetting phases are explicitly set. In this case, liquid water which uses the relations from IAPWS'97~\cite{IAPWS1997} is chosen as the wetting phase on line \ref{tutorial-coupled:wettingPhase} and a liquid model oil is chosen as the non-wetting phase on line \ref{tutorial-coupled:nonwettingPhase}. For all parameters that depend on space, such as the properties of the soil, the specific spatial  Bernd Flemisch committed Jul 16, 2010 253 parameters for the problem of interest are specified in line  Andreas Lauser committed Nov 10, 2010 254 255 \ref{tutorial-coupled:set-spatialparameters}. The final property, which is set on line \ref{tutorial-coupled:gravity}, is optional and tells the model not to  Bernd Flemisch committed Jul 16, 2010 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 use gravity. Parameters which are specific to a set-up -- like boundary and initial conditions, source terms or temperature within the domain -- but are required to solve the differential equations of the models are specified via a \textit{problem} class. If the two-phase box model is used, this class must be derived from \texttt{TwoPBoxProblem} as done on line \ref{tutorial-coupled:def-problem}. The problem class always has at least five methods: \begin{itemize} \item A method \texttt{boundaryTypes()} specifying the kind of boundary conditions to be used for a boundary segment \item A method \texttt{dirichlet()} specifying the actual values for the Dirichlet conditions on a boundary segment \item A method \texttt{neumann()} specifying the actual values for the Neumann conditions on a boundary segment  Klaus Mosthaf committed Oct 18, 2010 273 \item A method for source or sink terms called \texttt{source()}  Bernd Flemisch committed Jul 16, 2010 274 275 276 277 278 279 280 281 282 283 \item A method called \texttt{initial()} for specifying the initial condition. \end{itemize} Methods which make statements about boundary segments of the grid (i.e. \texttt{boundaryTypes()}, \texttt{dirichlet()} and \texttt{neumann()}) get six parameters. The first parameter differs if the type of the boundary condition is defined \texttt{boundaryTypes()}: \begin{description} \item[BCtypes:] A container which stores the type of the boundary condition  Melanie Darcis committed Oct 13, 2010 284 for each equation. For the typical case where all equations have the same boundary  Bernd Flemisch committed Jul 16, 2010 285 286 287 288 289 290 291 292 293 294 295 296 297 condition at a certain position, there are two methods that set the appropriate conditions for all primary variables / equations: Either \texttt{setAllDirichlet()} or \texttt{setAllNeumann()}. \item[element:] The element of the grid where the boundary segment is located. \item[fvElemGeometry:] The finite-volume geometry induced on the finite element by the box scheme. \item[isIt:] The \texttt{IntersectionIterator} of the boundary segement as given by the grid \item[scvIdx:] The index of the sub-control volume in \texttt{fvElementGeometry} adjacent to the boundary segment. \item[boundaryFaceIdx:] The index of the boundary face in \texttt{fvElementGeometry} which represents the boundary segment. \end{description}  Benjamin Faigle committed Feb 04, 2011 298 299 300 301 302 To ensure that no boundaries are undefined, a small safeguard value \texttt{eps\_} is usually added when compariing spatial coordinates. The left boundary is hence not assigned where the first coordinate is equal to zero, but where it is smaller than a very small value \texttt{eps\_}.  Bernd Flemisch committed Jul 16, 2010 303 304 305 306 307 After the type of the boundary condition is defined, their values have to be assigned with the methods \texttt{dirichlet()} and \texttt{neumann()} which only differ by the first function parameter: \begin{description} \item[values:] A vector which stores the result of the method. What  Philipp Nuske committed Dec 10, 2010 308  the values in this vector mean is dependent on the method: For  Bernd Flemisch committed Jul 16, 2010 309 310 311 312 313  \texttt{dirichlet()} it contains the values of the primary variables, for \texttt{neumann()} the mass fluxes per area unit over the boundary segment. \end{description}  Melanie Darcis committed Oct 13, 2010 314 Similarly, the \texttt{initial()} and \texttt{source()} methods  Bernd Flemisch committed Jul 16, 2010 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 specify properties of sub-control volumes and thus only get \texttt{values}, \texttt{element}, \texttt{fvElemGeom} and \texttt{scvIdx} as parameters. In addition to these five methods, there might be some model-specific methods. If the isothermal two-phase model is used, this includes for example a \texttt{temperature()} method which returns the temperature in Kelvin of the fluids and the rock matrix in the domain. This temperature is then used by the model to calculate fluid properties which possibly depend on it, e.g. density. \subsection{Defining fluid properties}\label{tutorial-coupled:description-fluid-class} The \Dumux distribution includes some common substances which can be used out of the box. The properties of the pure substances (such as the component  Melanie Darcis committed Oct 13, 2010 331 332 nitrogen, water, or pseudo-component air) are stored in header files in the folder \verb+dumux/material/components+. Each of these files  Bernd Flemisch committed Jul 16, 2010 333 334 335 defines a class with the same name as the component but starting with a capital letter, e.g. \texttt{Water}, and are derived from \texttt{Component}.  Melanie Darcis committed Oct 13, 2010 336 Most often, when two or more components are considered, fluid interactions  Bernd Flemisch committed Jul 16, 2010 337 338 such as solubility effects come into play and properties of mixtures such as the density are of interest. These interactions are defined in  Melanie Darcis committed Oct 13, 2010 339 a specific \verb+fluidsystem+ in the folder \verb+dumux/material/fluidsystems+.  Bernd Flemisch committed Jul 16, 2010 340 341 It features methods returning fluid properties like density, enthalpy, viscosity, etc. by accessing the pure components as well as binary coefficients such as  Klaus Mosthaf committed Oct 18, 2010 342 Henry or diffusion coefficients, which are stored in  Melanie Darcis committed Oct 13, 2010 343 344 \verb+dumux/material/binarycoefficients+. New fluids which are not yet available in the \Dumux distribution can be defined analogously.  Bernd Flemisch committed Jul 16, 2010 345   Klaus Mosthaf committed Oct 19, 2010 346 347 348 % In this example, a class for the definition of a two-phase system is used. This allows the choice % of the two components oil and water and to access the parameters that are relevant for the two-phase model.  Bernd Flemisch committed Jul 16, 2010 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 \subsection{The definition of the parameters that are dependent on space}\label{tutorial-coupled:description-spatialParameters} In \Dumux, the properties of the porous medium such as \textit{intrinsic permeability}, the \textit{porosity}, the \textit{heat capacity} as well as the \textit{heat conductivity} can be defined in space using a so-called \texttt{spatial parameters} class. However, because the soil also has an effect on the material laws of the fluids (e.g. \textit{capillarity}), their selection and definition of their attributes (e.g. \textit{residual saturations}) are also accomplished in the spatial parameters. The base class \texttt{Dumux::BoxSpatialParameters} holds a general averaging procedure for vertex-centered box-methods. Listing \ref{tutorial-coupled:spatialparametersfile} shows the file \verb+tutorialspatialparameters_coupled.hh+: \begin{lst}[File tutorial/tutorialspatialparameters\_coupled.hh]\label{tutorial-coupled:spatialparametersfile} \mbox{} \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,  Andreas Lauser committed Feb 04, 2011 367 numberstyle=\tiny, numbersep=5pt, firstline=27]{../../tutorial/tutorialspatialparameters_coupled.hh}  Bernd Flemisch committed Jul 16, 2010 368 369 370 371 372 \end{lst} First, a certain material law that best describes the problem at hand has to be selected in line \ref{tutorial-coupled:rawlaw}\label{tutorial-coupled:materialLaw}. \Dumux provides several material laws in the folder  Melanie Darcis committed Oct 13, 2010 373 \verb+dumux/material/fluidmatrixinteractions+.  Philipp Nuske committed Dec 10, 2010 374 The selected one -- here it is a relation according to a regularized version of Brooks \& Corey -- is included  Bernd Flemisch committed Jul 16, 2010 375 in line \ref{tutorial-coupled:rawLawInclude}. After the selection,  Philipp Nuske committed Dec 10, 2010 376 377 378 an adapter in line \ref{tutorial-coupled:eff2abs} translates between the law for effective values (the Brooks \& Corey model) and the saturations generated as simulations results, i.e. residual saturations are considered. As the applied raw law knows best which kind of parameters are necessary,  Melanie Darcis committed Nov 11, 2010 379 it provides a parameter class \texttt{RegularizedBrooksCoreyParams} that is  Bernd Flemisch committed Jul 16, 2010 380 381 accessible via the member \texttt{Params} and defined in line \ref{tutorial-coupled:matLawObjectType}. The material law object  Philipp Nuske committed Dec 10, 2010 382 is now instantiated correctly as a private object  Bernd Flemisch committed Jul 16, 2010 383 384 385 386 387 388 389 390 391 392 393 in line \ref{tutorial-coupled:matParamsObject}. In line \ref{tutorial-coupled:permeability} the function returning the intrinsic permeability can be found. As can be seen, the function has to be called with three different arguments. (\texttt{Element}) is again the current element, which also holds information about its geometry and position, the second argument (\texttt{fvElemGeom}) holds information about the finite-volume geometry induced by the box-method, and the third defines the index of the current sub-control volume. The intrinsic permeability is a tensor and is thus returned in form of a $\texttt{dim} \times \texttt{dim}$-matrix where \texttt{dim} is the dimension  Philipp Nuske committed Dec 10, 2010 394 395 of the problem.  Bernd Flemisch committed Jul 16, 2010 396 397 398 The function \texttt{porosity()} defined in line \ref{tutorial-coupled:porosity} is called with the same arguments as the permeability function described before and returns the porosity  Philipp Nuske committed Dec 10, 2010 399 400 dependent on the position in the domain.  Bernd Flemisch committed Jul 16, 2010 401 Next, the method \texttt{materialLawParams()} defines in line  Philipp Nuske committed Dec 10, 2010 402 403 404 405 406 \ref{tutorial-coupled:matLawParams} which \verb+materialLawParams+ object should be applied at this specific position. Although in this case only one objects is returned, in general the problem may be heterogeneous, demanding for different objects at different positions in space. While the selection of the type of this object was already explained (line \ref{tutorial-coupled:rawLawInclude}), some specific parameter  Bernd Flemisch committed Jul 16, 2010 407 values of the applied material law are still needed. This is  Melanie Darcis committed Oct 13, 2010 408 409 done in the constructor body (line \ref{tutorial-coupled:setLawParams}). Depending on the type of the \texttt{materialLaw} object, the adequate \texttt{set}-methods  Bernd Flemisch committed Jul 16, 2010 410 are provided by the object to access all necessary parameters  Philipp Nuske committed Dec 10, 2010 411 412 for the applied material law. The name of the access / set functions as well as the rest of the implementation of the material description can be found in  Benjamin Faigle committed Feb 04, 2011 413 \verb+dumux/material/fluidmatrixinteractions/2p+.  Bernd Flemisch committed Jul 16, 2010 414 415 416 417 418  \subsection{Exercises} \label{tutorial-coupled:exercises} The following exercises will give you the opportunity to learn how you can change soil parameters, boundary conditions and fluid properties  419 in \Dumux.  Bernd Flemisch committed Jul 16, 2010 420 421 422  \subsubsection{Exercise 1} \renewcommand{\labelenumi}{\alph{enumi})} For Exercise 1 you only have  Melanie Darcis committed Nov 11, 2010 423 424 to make some small changes in the tutorial files.  Bernd Flemisch committed Jul 16, 2010 425 426 \begin{enumerate}  Melanie Darcis committed Nov 11, 2010 427 \item \textbf{Run the Model} \\  Melanie Darcis committed Dec 07, 2010 428 To get an impression what the results should look like you can first run the original version of the coupled tutorial model by typing \texttt{./tutorial\_coupled 5e5 10}. The first number behind the simulation name defines the timespan of the simulation run in seconds, the second number defines the initial time step size. Note that the time step size is automatically optimized during the simulation. For the visualisation with paraview please refer to \ref{quick-start-guide}.\\  Melanie Darcis committed Nov 11, 2010 429 430  \item \textbf{Changing the Model Domain and the Boundary Conditions} \\  Bernd Flemisch committed Jul 16, 2010 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445  Change the size of the model domain so that you get a rectangle with edge lengths of $\text{x} = 400 m$ and $\text{y} = 500 m$ and with discretization lengths of $\Delta \text{x} = 20$ m and $\Delta \text{y} = 20$ m. Change the boundary conditions in the file \texttt{tutorialproblem\_coupled.hh} so that water enters from the bottom and oil is extracted from the top boundary. The right and the left boundary should be closed for water and oil fluxes. Compile the main file by typing \texttt{make tutorial\_coupled} and run the model. \item \textbf{Changing Fluids} \\  Benjamin Faigle committed Feb 04, 2011 446 447 448 449 Now you can change the fluids. Use DNAPL instead of Oil and Brine instead of Water. To do that you have to select different components via the property system in the problem file: \begin{enumerate} \item Brine: The class \texttt{Dumux::Brine} acts as a adapter to the fluid system that alters a pure water class by adding some salt. Hence, the class \texttt{Dumux::Brine} uses a pure water class, such as \texttt{Dumux::H2O}, as a second template argument after the data type \texttt{} as a template argument (be aware to use the complete water class with its own template parameter). \item DNAPL: A standard set of chemical substances, such as Oil and Brinde, is already included (via a list of \texttt{\#include ..} commandos) and hence easy accessible by default. This is not the case for the class \texttt{Dumux::SimpleDNAPL}, however, which is located in the folder \texttt{dumux/material/components/}. Try to include the file as well as select the component via the property system.  Melanie Darcis committed Nov 11, 2010 450 \end{enumerate}  Benjamin Faigle committed Feb 04, 2011 451 452 If you want to take a closer look how the fluid classes are defined and which substances are already available please browse through the files in the directory \texttt{/dumux/material/components}.  Bernd Flemisch committed Jul 16, 2010 453   Melanie Darcis committed Nov 11, 2010 454 \item \textbf{Use the \Dumux fluid system} \\  Philipp Nuske committed Dec 10, 2010 455 \Dumux usually organises fluid mixtures via a \texttt{fluidsystem}. In order to include a fluidsystem you first have to comment the lines \ref{tutorial-coupled:2p-system-start} to \ref{tutorial-coupled:2p-system-end} in the problem file. If you use eclipse, this can easily be done by pressing \textit{str + shift + 7} -- the same as to cancel the comment later on.\\  456 Now include the file \texttt{fluidsystems/h2o\_n2\_system.hh} in the material folder, and set a property \texttt{FluidSystem} with the appropriate type, \texttt{Dumux::H2O\_N2\_System}. However, the complicated fluidsystem uses tabularized fluid data, which need to be initialized in the constructor body of the current problem by adding \texttt{GET\_PROP\_TYPE(TypeTag, PTAG(FluidSystem))::init();}, hence using the initialization function of the applied fluidsystem. As water flow replacing a gas is much faster, test your simulation only until 2e3 seconds and start with a time step of 1 second.\\  Melanie Darcis committed Nov 11, 2010 457 Please reverse the changes of this example, as we still use bulk phases and hence do not need such an extensive fluid system.  Bernd Flemisch committed Jul 16, 2010 458 459  \item \textbf{Changing Constitutive Relations} \\  Philipp Nuske committed Dec 10, 2010 460 461  Use an unregularized linear law with an entry pressure of $p_e = 0.0$ and maximal capillary pressure of e.g. $p_{c_{max}} = 2000.0$ instead of using a regularized Brooks-Corey law for the  Melanie Darcis committed Nov 11, 2010 462  relative permeability and for the capillary pressure saturation relationship. To do that you have  Bernd Flemisch committed Jul 16, 2010 463 464  to change the file \texttt{tutorialspatialparameters\_coupled.hh}. You can find the material laws in the folder  Melanie Darcis committed Oct 13, 2010 465  \verb+dumux/material/fluidmatrixinteractions+. The necessary parameters  Melanie Darcis committed Nov 11, 2010 466 467 468 of the linear law and the respective \texttt{set}-functions can be found in the file \\ \verb+dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh+.  Bernd Flemisch committed Jul 16, 2010 469 470 471 472  \item \textbf{Heterogeneities} \\ Set up a model domain with the soil properties given in Figure \ref{tutorial-coupled:exercise1_d}. Adjust the boundary conditions  Melanie Darcis committed Nov 11, 2010 473  so that water is again flowing from the left to the right of the  Bernd Flemisch committed Jul 16, 2010 474 475 476 477 478 479 480 481 482 \begin{figure}[h] \psfrag{K1 =}{K $= 10^{-8}\text{ m}^2$} \psfrag{phi1 =}{$\phi = 0.15$} \psfrag{K2 =}{\textcolor{white}{K $= 10^{-9}\text{ m}^2$}} \psfrag{phi2 =}{\textcolor{white}{$\phi = 0.3$}} \psfrag{600 m}{600 m} \psfrag{300 m}{300 m} \centering \includegraphics[width=0.5\linewidth,keepaspectratio]{EPS/exercise1_c.eps}  Melanie Darcis committed Nov 11, 2010 483 \caption{Exercise 1f: Set-up of a model domain with a heterogeneity. $\Delta \text{x} = 20$ m $\Delta \text{y} = 20$ m.}\label{tutorial-coupled:exercise1_d}  Bernd Flemisch committed Jul 16, 2010 484 \end{figure}  Benjamin Faigle committed Feb 04, 2011 485 486 domain. You can use the fluids of exercise 1c).\\ Hint: The current position of the element can be obtained via \texttt{element.geometry().center();}.\\  Melanie Darcis committed Nov 11, 2010 487 When does the front cross the material border? In paraview, the option \textit{View} $\rightarrow$ \textit{Animation View} is nice to get a rough feeling of the timestep sizes.  Bernd Flemisch committed Jul 16, 2010 488 489 490 491 \end{enumerate} \subsubsection{Exercise 2} For this exercise you should create a new proplem file analogous to  Benjamin Faigle committed Feb 04, 2011 492 493 the file \texttt{tutorialproblem\_coupled.hh} (e.g. with the name \texttt{ex2\_tutorialproblem\_coupled.hh} and new spatial parameters  Melanie Darcis committed Oct 13, 2010 494 just like \texttt{tutorialspatialparameters\_coupled.hh}. The new problem file needs to  Benjamin Faigle committed Feb 04, 2011 495 496 497 498 499 500 501 502 be included in the file \texttt{tutorial\_coupled.cc}.\\ The new files should contain the definition of a new classes with names that relate to the file name, such as \texttt{Ex2TutorialProblemCoupled}. Make sure that you also adjust the guardian macros in lines \ref{tutorial-coupled:guardian1} and \ref{tutorial-coupled:guardian1} in the header files (e.g. change \\ \texttt{DUMUX\_TUTORIALPROBLEM\_COUPLED\_HH} to \texttt{DUMUX\_EX2\_TUTORIALPROBLEM\_COUPLED\_HH}). Besides also adjusting the guardian macros,  Bernd Flemisch committed Jul 16, 2010 503 the new problem file should define and use a new type tag for the problem as well as a new problem class  Benjamin Faigle committed Feb 04, 2011 504 e.g. \texttt{Ex2TutorialProblemCoupled}. Make sure you assign your newly defined spatial  Bernd Flemisch committed Jul 16, 2010 505 parameter class to the \texttt{SpatialParameters} property for the new  Melanie Darcis committed Nov 11, 2010 506 type tag. \\  Bernd Flemisch committed Jul 16, 2010 507 508 509 510 511 After this, change the \texttt{create()} method of the \texttt{Grid} property so that it matches the domain described by figure \ref{tutorial-coupled:ex2_Domain}. Adapt the problem class so that the boundary conditions are consistent with figure \ref{tutorial-coupled:ex2_BC}. Initially the domain is fully saturated  Philipp Nuske committed Dec 10, 2010 512 with water and the pressure is $p_w = 5 \times 10^5 \text{Pa}$. Oil  Bernd Flemisch committed Jul 16, 2010 513 514 infiltrates from the left side. Create a grid with $20$ cells in $x$-direction and $10$ cells in $y$-direction. The simulation time  Melanie Darcis committed Nov 11, 2010 515 516 should be set to $1\times 10^6 \text{ s}$ with an initial time step size of $100 \text{ s}$.  Bernd Flemisch committed Jul 16, 2010 517 518 519 520 521 522 523 524 525  Now include your new problem file in the main file and replace the \texttt{TutorialProblemCoupled} type tag by the one you've created and compile the program. \begin{figure}[h] \psfrag{K1}{K $= 10^{-7}\text{ m}^2$} \psfrag{phi1}{$\phi = 0.2$}  Klaus Mosthaf committed Nov 04, 2010 526 \psfrag{Lin}{Brooks-Corey Law}  Philipp Nuske committed Dec 10, 2010 527 \psfrag{Lin2}{$\lambda = 1.8$, $p_e = 1000$}  Bernd Flemisch committed Jul 16, 2010 528 529 \psfrag{K2}{K $= 10^{-9}\text{ m}^2$} \psfrag{phi2}{$\phi = 0.15$}  Klaus Mosthaf committed Nov 04, 2010 530 \psfrag{BC1}{Brooks-Corey Law}  Philipp Nuske committed Dec 10, 2010 531 \psfrag{BC2}{$\lambda = 2$, $p_e = 1500$}  Bernd Flemisch committed Jul 16, 2010 532 533 534 535 536 537 538 539 540 541 542 543 \psfrag{H1y}{50 m} \psfrag{H2y}{15 m} \psfrag{H3y}{20 m} \psfrag{L1x}{100 m} \psfrag{L2x}{50 m} \psfrag{L3x}{25 m} \centering \includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/Ex2_Domain.eps} \caption{Set-up of the model domain and the soil parameters}\label{tutorial-coupled:ex2_Domain} \end{figure} \begin{figure}[h]  Philipp Nuske committed Dec 10, 2010 544 \psfrag{pw}{$p_w = 5 \times 10^5$ [\text{Pa}]}  Bernd Flemisch committed Jul 16, 2010 545 546 547 548 549 550 551 552 553 \psfrag{S}{$S_n = 1.0$} \psfrag{qw}{$q_w = 2 \times 10^{-4}$ [kg/$\text{m}^2$s]} \psfrag{qo}{$q_n = 0.0$ [kg/$\text{m}^2$s]} \psfrag{no flow}{no flow} \centering \includegraphics[width=0.8\linewidth,keepaspectratio]{EPS/Ex2_Boundary.eps} \caption{Boundary Conditions}\label{tutorial-coupled:ex2_BC} \end{figure}  Melanie Darcis committed Nov 11, 2010 554 555 556 557 558 \begin{itemize} \item Increase the simulation time to e.g. $4\times 10^7 \text{ s}$. Investigate the saturation: Is the value range reasonable? \item What happens if you increase the resolution of the grid? \end{itemize}  Bernd Flemisch committed Jul 16, 2010 559 560 561 562 \subsubsection{Exercise 3} Create a new file for benzene called \texttt{benzene.hh} and implement a new fluid system. (You may get a hint by looking at existing fluid  Melanie Darcis committed Nov 11, 2010 563 systems in the directory \verb+/dumux/material/fluidsystems+.) \\  Bernd Flemisch committed Jul 16, 2010 564 565 566 567 Use benzene as a new fluid and run the model of Exercise 2 with water and benzene. Benzene has a density of $889.51 \, \text{kg} / \text{m}^3$ and a viscosity of $0.00112 \, \text{Pa} \; \text{s}$.  568 \clearpage \newpage  569 570 571 572 %%% Local Variables: %%% mode: latex %%% TeX-master: "dumux-handbook" %%% End: