tutorial-coupled.tex 30.9 KB
Newer Older
Karin Erbertseder's avatar
Karin Erbertseder committed
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:

\begin{equation}
 	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 .
\end{equation}

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:

\begin{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}
\end{equation}

\begin{equation}
	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
\end{equation}
 
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}
\begin{equation}
\label{eq:p} 
	\tilde p = \sum_i N_i \hat{p_i}
\end{equation}
\begin{equation}
\label{eq:v} 
	\tilde{\mathbf v} = \sum_i N_i \hat{\mathbf v}
\end{equation}
\begin{equation}
\label{eq:x} 
	\tilde x^\kappa  = \sum_i N_i \hat x^\kappa 
\end{equation}
\end{minipage}
\hfill
\begin{minipage}[b]{0.5\textwidth}
\begin{equation}
\label{eq:dp} 
	\nabla \tilde p = \sum_i \nabla N_i \hat{p_i}
\end{equation}
\begin{equation}
\label{eq:dv} 
	\nabla \tilde{\mathbf v} = \sum_i \nabla N_i \hat{\mathbf v}
\end{equation}
\begin{equation}
\label{eq:dx} 
	\nabla \tilde x^\kappa  = \sum_i \nabla N_i \hat x^\kappa .
\end{equation}
\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.

\begin{equation}
	f(u) = 0  \qquad \Rightarrow \qquad f(\tilde u) = \varepsilon
\end{equation}

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,

\begin{equation}
	\int_G W_j \cdot \varepsilon \: \overset {!}{=} \: 0 \qquad \textrm{with} \qquad \sum_j W_j =1
\end{equation}

yields the following equation:

\begin{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 .
\end{equation}

Then, the chain rule and the \textsc{green-gaussian} integral theorem are applied.

\begin{equation}
	\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
\end{equation}

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:

\begin{equation}
	 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}
\end{equation}

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

\begin{equation}
	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
\end{equation}

Defining the weighting function $W_j$ to be piecewise constant over a control volume box $B_i$ 

\begin{equation}
	W_j(x) = \begin{cases}
	          1 &x \in B_i \\
		  0 &x \notin B_i\\
	         \end{cases}
\end{equation}

causes $\nabla W_j = 0$:

\begin{equation}
\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 .
\end{equation}

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:

\begin{equation}
\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 
\end{equation}

Karin Erbertseder's avatar
   
Karin Erbertseder committed
125
\section[Fully-coupled model]{Solving a problem using a Fully-Coupled Model}\label{tutorial-coupled}
Bernd Flemisch's avatar
Bernd Flemisch committed
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.
130
 \item Material properties and constitutive relationships have to be selected.
Bernd Flemisch's avatar
Bernd Flemisch committed
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's avatar
Bernd Flemisch committed
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]$}
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's avatar
Bernd Flemisch committed
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}
159
  -
160
  \nabla \cdot \left( \varrho_{w} \, \frac{k_{rw}}{\mu_{w}} \, \mathbf{K}\;\nabla p_w \right)
Bernd Flemisch's avatar
Bernd Flemisch committed
161
  -
162
  q_w
Bernd Flemisch's avatar
Bernd Flemisch committed
163
164
165
166
  & =
  0 \\
  \label{massbalanceoil}
  \frac {\partial (\phi \, S_{o}\, \varrho_{o})}{\partial t}
167
  -
168
  \nabla \cdot \left( \varrho_{o} \, \frac{k_{ro}}{\mu_{o}} \, \mathbf{K}\;\nabla p_o \right)
Bernd Flemisch's avatar
Bernd Flemisch committed
169
  -
170
  q_o 
Bernd Flemisch's avatar
Bernd Flemisch committed
171
172
173
174
175
176
177
  & =
  0
\end{align}

\subsection{The main file}

Listing \ref{tutorial-coupled:mainfile} shows the main file
178
\texttt{tutorial/tutorial\_coupled.cc} for the coupled two-phase
Bernd Flemisch's avatar
Bernd Flemisch committed
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,
184
  numberstyle=\tiny, numbersep=5pt, firstline=28]{../../tutorial/tutorial_coupled.cc}
Bernd Flemisch's avatar
Bernd Flemisch committed
185
186
187
\end{lst}

From line \ref{tutorial-coupled:include-begin} to line
188
\ref{tutorial-coupled:include-end} the \Dune and \Dumux files which
Bernd Flemisch's avatar
Bernd Flemisch committed
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's avatar
Bernd Flemisch committed
198
199

The first thing which should be done at run time is to initialize the
200
message passing interface using \Dune's \texttt{MPIHelper} class. Line
201
\ref{tutorial-coupled:init-mpi} is essential if the simulation is
Bernd Flemisch's avatar
Bernd Flemisch committed
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
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's avatar
Bernd Flemisch committed
208

209
After this, a grid is created in line
Bernd Flemisch's avatar
Bernd Flemisch committed
210
\ref{tutorial-coupled:create-grid} and the problem is instantiated for
211
its leaf grid view in line \ref{tutorial-coupled:instantiate-problem}.
212
Finally, on line \ref{tutorial-coupled:begin-restart} a state written to
213
214
disk by a previous simulation run is restored if requested by the user.
The simulation procedure is started at line
Bernd Flemisch's avatar
Bernd Flemisch committed
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,
226
numberstyle=\tiny, numbersep=5pt, firstline=27]{../../tutorial/tutorialproblem_coupled.hh}
Bernd Flemisch's avatar
Bernd Flemisch committed
227
228
\end{lst}

229
First, a new type tag is created for the problem in line
Bernd Flemisch's avatar
Bernd Flemisch committed
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
235
is going to be used is defined in line \ref{tutorial-coupled:set-grid} --
236
in this case that is  \texttt{SGrid}.  Since there's no uniform
237
mechanism to allocate grids in \Dune, the \texttt{Grid} property also contains
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's avatar
Bernd Flemisch committed
242
the appropriate fluid system, that specifies both information about
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's avatar
Bernd Flemisch committed
253
parameters for the problem of interest are specified in line
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's avatar
Bernd Flemisch committed
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
273
\item A method for source or sink terms called \texttt{source()}
Bernd Flemisch's avatar
Bernd Flemisch committed
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's avatar
Melanie Darcis committed
284
for each equation. For the typical case where all equations have the same boundary
Bernd Flemisch's avatar
Bernd Flemisch committed
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}
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's avatar
Bernd Flemisch committed
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
308
  the values in this vector mean is dependent on the method: For
Bernd Flemisch's avatar
Bernd Flemisch committed
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's avatar
Melanie Darcis committed
314
Similarly, the \texttt{initial()} and \texttt{source()} methods
Bernd Flemisch's avatar
Bernd Flemisch committed
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's avatar
Melanie Darcis committed
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's avatar
Bernd Flemisch committed
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's avatar
Melanie Darcis committed
336
Most often, when two or more components are considered, fluid interactions 
Bernd Flemisch's avatar
Bernd Flemisch committed
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's avatar
Melanie Darcis committed
339
a specific \verb+fluidsystem+ in the folder \verb+dumux/material/fluidsystems+.
Bernd Flemisch's avatar
Bernd Flemisch committed
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
342
Henry or diffusion coefficients, which are stored in 
Melanie Darcis's avatar
Melanie Darcis committed
343
344
\verb+dumux/material/binarycoefficients+. New fluids which are not yet
 available in the \Dumux distribution can be defined analogously.
Bernd Flemisch's avatar
Bernd Flemisch committed
345

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's avatar
Bernd Flemisch committed
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<TypeTag>} 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,
367
numberstyle=\tiny, numbersep=5pt, firstline=27]{../../tutorial/tutorialspatialparameters_coupled.hh}
Bernd Flemisch's avatar
Bernd Flemisch committed
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's avatar
Melanie Darcis committed
373
\verb+dumux/material/fluidmatrixinteractions+.
374
The selected one -- here it is a relation according to a regularized version of Brooks \& Corey --  is included
Bernd Flemisch's avatar
Bernd Flemisch committed
375
in line \ref{tutorial-coupled:rawLawInclude}. After the selection,
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,
379
it provides a parameter class \texttt{RegularizedBrooksCoreyParams} that is
Bernd Flemisch's avatar
Bernd Flemisch committed
380
381
accessible via the member \texttt{Params} and defined in line 
\ref{tutorial-coupled:matLawObjectType}. The material law object 
382
is now instantiated correctly as a private object
Bernd Flemisch's avatar
Bernd Flemisch committed
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 
394
395
of the problem. 

Bernd Flemisch's avatar
Bernd Flemisch committed
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
399
400
dependent on the position in the domain.

Bernd Flemisch's avatar
Bernd Flemisch committed
401
Next, the method \texttt{materialLawParams()} defines in line 
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's avatar
Bernd Flemisch committed
407
values of the applied material law are still needed. This is 
Melanie Darcis's avatar
Melanie Darcis committed
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's avatar
Bernd Flemisch committed
410
are provided by the object to access all necessary parameters 
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 
413
\verb+dumux/material/fluidmatrixinteractions/2p+.
Bernd Flemisch's avatar
Bernd Flemisch committed
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's avatar
Bernd Flemisch committed
420
421
422

\subsubsection{Exercise 1}
\renewcommand{\labelenumi}{\alph{enumi})} For Exercise 1 you only have
423
424
to make some small changes in the tutorial files.  

Bernd Flemisch's avatar
Bernd Flemisch committed
425
426
\begin{enumerate}

427
\item \textbf{Run the Model} \\
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}.\\
429
430

\item \textbf{Changing the Model Domain and the Boundary Conditions} \\
Bernd Flemisch's avatar
Bernd Flemisch committed
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} \\
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{<Scalar>} 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.
450
\end{enumerate}
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's avatar
Bernd Flemisch committed
453

454
\item \textbf{Use the \Dumux fluid system} \\
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<TypeTag>}. 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.\\
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's avatar
Bernd Flemisch committed
458
459

\item \textbf{Changing Constitutive Relations} \\
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
462
  relative permeability and for the capillary pressure saturation relationship. To do that you have
Bernd Flemisch's avatar
Bernd Flemisch committed
463
464
  to change the file \texttt{tutorialspatialparameters\_coupled.hh}. 
 You can find the material laws in the folder 
Melanie Darcis's avatar
Melanie Darcis committed
465
  \verb+dumux/material/fluidmatrixinteractions+. The necessary parameters
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's avatar
Bernd Flemisch committed
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
473
  so that water is again flowing from the left to the right of the
Bernd Flemisch's avatar
Bernd Flemisch committed
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}
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's avatar
Bernd Flemisch committed
484
\end{figure}
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();}.\\
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's avatar
Bernd Flemisch committed
488
489
490
491
\end{enumerate}

\subsubsection{Exercise 2}
For this exercise you should create a new proplem file analogous to
492
493
the file \texttt{tutorialproblem\_coupled.hh} (e.g. with the name 
\texttt{ex2\_tutorialproblem\_coupled.hh} and new spatial parameters 
Melanie Darcis's avatar
Melanie Darcis committed
494
just like \texttt{tutorialspatialparameters\_coupled.hh}. The new problem file needs to
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's avatar
Bernd Flemisch committed
503
the new problem file should define and use a new type tag for the problem as well as a new problem class
504
e.g. \texttt{Ex2TutorialProblemCoupled}. Make sure you assign your newly defined spatial 
Bernd Flemisch's avatar
Bernd Flemisch committed
505
parameter class to the \texttt{SpatialParameters} property for the new 
506
type tag. \\
Bernd Flemisch's avatar
Bernd Flemisch committed
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
512
with water and the pressure is $p_w = 5 \times 10^5 \text{Pa}$. Oil
Bernd Flemisch's avatar
Bernd Flemisch committed
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
515
516
should be set to $1\times 10^6 \text{ s}$ with an initial time step size of
$100 \text{ s}$.
Bernd Flemisch's avatar
Bernd Flemisch committed
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$}
526
\psfrag{Lin}{Brooks-Corey Law}
527
\psfrag{Lin2}{$\lambda = 1.8$, $p_e = 1000$}
Bernd Flemisch's avatar
Bernd Flemisch committed
528
529
\psfrag{K2}{K $= 10^{-9}\text{ m}^2$}
\psfrag{phi2}{$\phi = 0.15$}
530
\psfrag{BC1}{Brooks-Corey Law} 
531
\psfrag{BC2}{$\lambda = 2$, $p_e = 1500$}
Bernd Flemisch's avatar
Bernd Flemisch committed
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]
544
\psfrag{pw}{$p_w = 5 \times 10^5$ [\text{Pa}]}
Bernd Flemisch's avatar
Bernd Flemisch committed
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}

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's avatar
Bernd Flemisch committed
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 
563
systems in the directory \verb+/dumux/material/fluidsystems+.) \\
Bernd Flemisch's avatar
Bernd Flemisch committed
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: