adr.tex 20 KB
Newer Older
1
\chapter{Advection-Diffusion-Reaction Solver}
Daniele de Grazia's avatar
Daniele de Grazia committed
2 3 4 5 6 7

%3.4/UserGuide/Tutorial/ADRSolver
%3.4/UserGuide/Examples/ADRSolver/1DAdvection
%3.4/UserGuide/Examples/ADRSolver/3DAdvectionMassTransport
%3.4/UserGuide/Examples/ADRSolver/Helmholtz2D

8
\section{Synopsis}
Daniele de Grazia's avatar
Daniele de Grazia committed
9 10 11 12 13 14 15 16 17 18 19 20 21

The ADRSolver is designed to solve partial differential equations of the form:
\begin{equation}
\alpha \dfrac{\partial u}{\partial t} + \lambda u + \nu \nabla u + \epsilon \nabla \cdot (D \nabla u) = f
\end{equation}
in either discontinuous or continuous projections of the solution field. 
For a full list of the equations which are supported, and the capabilities of each equation, 
see the table below.

\begin{table}[h!]
\begin{center}
\tiny
\renewcommand\arraystretch{2.2} 
22 23 24 25 26
\begin{tabular}{llll}
\toprule
\textbf{Equation to solve} & \textbf{EquationType} & \textbf{Dimensions}   &
\textbf{Projections} \\
\midrule
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
$\nabla^2 u = 0$ & 
    \inltt{Laplace} & All &  Continuous/Discontinuous	\\
$\nabla^2 u  =  f$ & 
    \inltt{Poisson} & All &  Continuous/Discontinuous	\\
$\nabla^2 u  + \lambda u =  f$ & 
    \inltt{Helmholtz} & All & Continuous/Discontinuous \\
$\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ & 
    \inltt{SteadyAdvectionDiffusion} & 2D only & Continuous/Discontinuous \\
$\epsilon \nabla^2 u +  \lambda u = f$ & 
    \inltt{SteadyDiffusionReaction} & 2D only &  Continuous/Discontinuous \\
$\epsilon \nabla^2 u  \mathbf{V}\nabla u + \lambda u = f$ & 
    \inltt{SteadyAdvectionDiffusionReaction} & 2D only & 
    Continuous/Discontinuous \\
$ \dfrac{\partial u}{\partial t} + \mathbf{V}\nabla u = f$ &
    \inltt{UnsteadyAdvection} & All & Continuous/Discontinuous \\
$\dfrac{\partial u}{\partial t}  = \epsilon \nabla^2 u$ & 
    \inltt{UnsteadyDiffusion} & All & Continuous/Discontinuous \\
$\dfrac{\partial u}{\partial t}  + \mathbf{V}\nabla u = \epsilon \nabla^2 u$ &
    \inltt{UnsteadyAdvectionDiffusion} & All & Continuous/Discontinuous \\
$\dfrac{\partial u}{\partial t}  + u\nabla u =  0$ &
    \inltt{UnsteadyInviscidBurger} & 1D only & Continuous/Discontinuous \\
48
\bottomrule
Daniele de Grazia's avatar
Daniele de Grazia committed
49 50 51 52 53 54
\end{tabular}
\end{center}
\caption{Equations supported by the ADRSolver with their capabilities.}
\label{t:ADR1}
\end{table}

55
\section{Usage}
Daniele de Grazia's avatar
Daniele de Grazia committed
56

57
\begin{lstlisting}[style=BashInputStyle]
Daniele de Grazia's avatar
Daniele de Grazia committed
58
ADRSolver session.xml
59
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
60

61
\section{Session file configuration}
Daniele de Grazia's avatar
Daniele de Grazia committed
62 63 64 65 66

The type of equation which is to be solved is specified through the EquationType 
SOLVERINFO option in the session file. This can be set as in table \ref{t:ADR1}.
At present, the Steady non-symmetric solvers cannot be used in parallel. \\

67
\subsection{Solver Info}
Daniele de Grazia's avatar
Daniele de Grazia committed
68 69 70 71 72 73 74 75

The solver info are listed below:
\begin{itemize}
\item \textbf{Eqtype}: This sets the type of equation to solve, according to the table above.
\item \textbf{TimeIntegrationMethod}: The following types of time integration methods have been tested with each solver:
\begin{center}
\footnotesize
\renewcommand\arraystretch{1.2} 
76 77
\begin{tabular}{lcccc}
\toprule
78 79
\textbf{EqType} & \textbf{Explicit} & \textbf{Diagonally Implicit} &
    \textbf{ IMEX} & \textbf{Implicit}   \\
80
\midrule
81 82 83 84
\inltt{UnsteadyAdvection} &  \checkmark 	 &	 & 	&\\
\inltt{UnsteadyDiffusion} &  \checkmark 	 & \checkmark 	 & 	&\\
\inltt{UnsteadyAdvectionDiffusion} &  	 &   	& \checkmark	&\\
\inltt{UnsteadyInviscidBurger} &  \checkmark 	 &	 & 	&\\
85
\bottomrule
Daniele de Grazia's avatar
Daniele de Grazia committed
86 87
\end{tabular}
\end{center}
88

Daniele de Grazia's avatar
Daniele de Grazia committed
89
\item \textbf{Projection}: The Galerkin projection used may be either:
Daniele de Grazia's avatar
Daniele de Grazia committed
90
\begin{itemize}
91 92
	\item \inltt{Continuous} for a C0-continuous Galerkin (CG) projection.
	\item \inltt{Discontinuous} for a discontinous Galerkin (DG) projection.
Daniele de Grazia's avatar
Daniele de Grazia committed
93
\end{itemize}
Daniele de Grazia's avatar
Daniele de Grazia committed
94
\item \textbf{DiffusionAdvancement}: This specifies how to treat the diffusion term. This will be restricted by the choice of time integration scheme:
Daniele de Grazia's avatar
Daniele de Grazia committed
95
\begin{itemize}
96 97
	\item \inltt{Explicit} Requires the use of an explicit time integration
	scheme.
98
	\item \inltt{Implicit} Requires the use of a diagonally implicit, IMEX or
99
	Implicit scheme.
Daniele de Grazia's avatar
Daniele de Grazia committed
100
\end{itemize}
Daniele de Grazia's avatar
Daniele de Grazia committed
101
\item \textbf{AdvectionAdvancement}: This specifies how to treat the advection term. This will be restricted by the choice of time integration scheme:
Daniele de Grazia's avatar
Daniele de Grazia committed
102
\begin{itemize}
103 104 105
	\item \inltt{Explicit} Requires the use of an explicit or IMEX time integration
	scheme.
	\item \inltt{Implicit} Not supported at present.
Daniele de Grazia's avatar
Daniele de Grazia committed
106
\end{itemize}
Daniele de Grazia's avatar
Daniele de Grazia committed
107
\item \textbf{AdvectionType}: Specifies the type of advection:
Daniele de Grazia's avatar
Daniele de Grazia committed
108
\begin{itemize}
109 110
	\item \inltt{NonConservative} (for CG only).
	\item \inltt{WeakDG} (for DG only).
Daniele de Grazia's avatar
Daniele de Grazia committed
111 112 113
\end{itemize}
\item \textbf{DiffusionType}:
\begin{itemize}
114
	\item \inltt{LDG}.
Daniele de Grazia's avatar
Daniele de Grazia committed
115 116 117
\end{itemize}
\item \textbf{UpwindType}:
\begin{itemize}
118
	\item \inltt{Upwind}.
Daniele de Grazia's avatar
Daniele de Grazia committed
119 120 121
\end{itemize}
\end{itemize}

122
\subsection{Parameters}
Daniele de Grazia's avatar
Daniele de Grazia committed
123

124 125
The following parameters can be specified in the \inltt{PARAMETERS} section of
the session file:
Daniele de Grazia's avatar
Daniele de Grazia committed
126
\begin{itemize}
127
\item \inltt{epsilon}: sets the diffusion coefficient $\epsilon$.\\ 
Daniele de Grazia's avatar
Daniele de Grazia committed
128 129
\textit{Can be used} in: SteadyDiffusionReaction, SteadyAdvectionDiffusionReaction, UnsteadyDiffusion, UnsteadyAdvectionDiffusion. \\
\textit{Default value}: 0.
130 131
\item  \inltt{d00}, \inltt{d11}, \inltt{d22}: sets the diagonal entries of the
diffusion tensor $D$. \\
Daniele de Grazia's avatar
Daniele de Grazia committed
132 133
\textit{Can be used in}: UnsteadyDiffusion \\
\textit{Default value}: All set to 1 (i.e. identity matrix). 
134
\item  \inltt{lambda}: sets the reaction coefficient  $\lambda$. \\
Daniele de Grazia's avatar
Daniele de Grazia committed
135 136 137 138
\textit{Can be used in}: SteadyDiffusionReaction, Helmholtz, SteadyAdvectionDiffusionReaction\\
\textit{Default value}: 0.
\end{itemize}

139
\subsection{Functions}
Daniele de Grazia's avatar
Daniele de Grazia committed
140

141 142
The following functions can be specified inside the \inltt{CONDITIONS} section
of the session file:
Daniele de Grazia's avatar
Daniele de Grazia committed
143 144

\begin{itemize}
145 146 147
\item \inltt{AdvectionVelocity}: specifies the advection velocity $\mathbf{V}$.
\item \inltt{InitialConditions}: specifies the initial condition for unsteady problems.
\item \inltt{Forcing}: specifies the forcing function f.
Daniele de Grazia's avatar
Daniele de Grazia committed
148 149
\end{itemize}

150
\section{Examples}
151 152
Example files for the ADRSolver are provided in
\inltt{solvers/ADRSolver/Examples}
153

Daniele de Grazia's avatar
Daniele de Grazia committed
154
\subsection{1D Advection equation}
Daniele de Grazia's avatar
Daniele de Grazia committed
155

156 157
In this example, it will be demonstrated how the Advection equation can be
solved on a one-dimensional domain.
Daniele de Grazia's avatar
Daniele de Grazia committed
158

159
\subsubsection{Advection equation}
Daniele de Grazia's avatar
Daniele de Grazia committed
160 161 162 163 164 165
We consider the hyperbolic partial differential equation:
\begin{equation}
\dfrac{\partial u}{\partial t} + \dfrac{\partial f}{\partial x} = 0,
\end{equation}
where $f =  a u$ i s the advection flux.

166 167
\subsubsection{Input file}
The input for this example is given in the example file \inlsh{Advection1D.xml}
Daniele de Grazia's avatar
Daniele de Grazia committed
168

169 170 171
The geometry section defines a 1D domain consisting of $10$ segments. On each
segment an expansion consisting of $4$ Lagrange polynomials on the
Gauss-Lobotto-Legendre points is used as specified by
Daniele de Grazia's avatar
Daniele de Grazia committed
172
\begin{lstlisting}[style=XMLStyle]
173 174 175
<EXPANSIONS>
    <E COMPOSITE="C[0]" FIELDS="u" TYPE="GLL_LAGRANGE_SEM" NUMMODES="4"/>
</EXPANSIONS>
Daniele de Grazia's avatar
Daniele de Grazia committed
176 177
\end{lstlisting}

178 179 180
Since we are solving the unsteady advection problem, we must specify this in the
solver information. We also choose to use a discontinuous flux-reconstruction
projection and use a Runge-Kutta order 4 time-integration scheme.
Daniele de Grazia's avatar
Daniele de Grazia committed
181
\begin{lstlisting}[style=XMLStyle]
182 183 184 185 186
<I PROPERTY="EQTYPE"                VALUE="UnsteadyAdvection"   />
<I PROPERTY="Projection"            VALUE="DisContinuous"       />
<I PROPERTY="AdvectionType"         VALUE="FRDG"                />
<I PROPERTY="UpwindType"            VALUE="Upwind"              />
<I PROPERTY="TimeIntegrationMethod" VALUE="ClassicalRungeKutta4"/>
Daniele de Grazia's avatar
Daniele de Grazia committed
187 188
\end{lstlisting}

189 190 191 192 193 194 195
We choose to advect our solution for $20$ time units with a time-step of $0.01$
and so provide the following parameters
\begin{lstlisting}[style=XMLStyle]
<P> FinTime         = 20                   </P>
<P> TimeStep        = 0.01                 </P>
<P> NumSteps        = FinTime/TimeStep     </P>
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
196

197 198 199 200 201 202
We also specify the advection velocity. We first define dummy parameters
\begin{lstlisting}[style=XMLStyle]
<P> advx            = 1                     </P>
<P> advy            = 0                     </P>
\end{lstlisting}
and then define the actual advection function as
Daniele de Grazia's avatar
Daniele de Grazia committed
203
\begin{lstlisting}[style=XMLStyle]
204 205 206
<FUNCTION NAME="AdvectionVelocity">
    <E VAR="Vx" VALUE="advx" />
</FUNCTION>
Daniele de Grazia's avatar
Daniele de Grazia committed
207 208
\end{lstlisting}

209 210 211 212 213 214 215
Two boundary regions are defined, one at each end of the domain, and periodicity
is enforced
\begin{lstlisting}[style=XMLStyle] 
<BOUNDARYREGIONS>
    <B ID="0"> C[1] </B>
    <B ID="1"> C[2] </B>
</BOUNDARYREGIONS>
Daniele de Grazia's avatar
Daniele de Grazia committed
216

217 218 219 220 221 222 223 224 225
<BOUNDARYCONDITIONS>
    <REGION REF="0">
        <P VAR="u" VALUE="[1]" />
    </REGION>
    <REGION REF="1">
        <P VAR="u" VALUE="[0]" />
    </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
226

227 228 229 230 231 232 233 234 235 236
Finally, we specify the initial value of the solution on the domain
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="InitialConditions">
    <E VAR="u" VALUE="exp(-20.0*x*x)" />
</FUNCTION>

<FUNCTION NAME="ExactSolution">
    <E VAR="u" VALUE="exp(-20.0*x*x)" />
</FUNCTION>
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
237

238 239 240 241
\subsubsection{Running the code}
\begin{lstlisting}[style=BashInputStyle]
ADRSolver Advection1D.xml
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
242

243 244
To visualise the output, we can convert it into either TecPlot or VTK formats
\begin{lstlisting}[style=BashInputStyle]
245 246
FieldConvert Advection1D.xml Advection1D.fld Advection1D.dat
FieldConvert Advection1D.xml Advection1D.fld Advection1D.vtu
247 248
\end{lstlisting}

Daniele de Grazia's avatar
Daniele de Grazia committed
249

Daniele de Grazia's avatar
Daniele de Grazia committed
250
\subsection{2D Helmholtz Problem}
Daniele de Grazia's avatar
Daniele de Grazia committed
251

252
In this example, it will be demonstrated how the Helmholtz equation can be solved on a two-dimensional domain.
Daniele de Grazia's avatar
Daniele de Grazia committed
253

254
\subsubsection{Helmholtz equation}
Daniele de Grazia's avatar
Daniele de Grazia committed
255 256 257 258 259 260 261 262 263

We consider the elliptic partial differential equation:

\begin{equation}
\nabla^2 u  + \lambda u =  f
\end{equation}

where $\nabla^2$ is the Laplacian and $\lambda$ is a real positive constant.

264
\subsubsection{Input file} 
Daniele de Grazia's avatar
Daniele de Grazia committed
265

266 267
The input for this example is given in the example file
\inlsh{Helmholtz2D\_modal.xml}
Daniele de Grazia's avatar
Daniele de Grazia committed
268

269 270 271 272
The geometry for this problem is a two-dimensional octagonal plane containing
both triangles and quadrilaterals. Note that a mesh composite may only contain
one type of element. Therefore, we define two composites for the domain, while
the rest are used for enforcing boundary conditions.
Daniele de Grazia's avatar
Daniele de Grazia committed
273
\begin{lstlisting}[style=XMLStyle]
274 275 276 277 278 279 280 281 282 283
<COMPOSITE>
     <C ID="0"> Q[22-47] </C>
     <C ID="1"> T[0-21] </C>
     <C ID="2"> E[0-1] </C>
     .
     .
     <C ID="10"> E[84,75,69,62,51,40,30,20,6] </C>
</COMPOSITE>

<DOMAIN> C[0-1] </DOMAIN>
Daniele de Grazia's avatar
Daniele de Grazia committed
284 285
\end{lstlisting}

286 287
For both the triangular and quadrilateral elements, we use the modified Legendre
basis with $7$ modes (maximum polynomial order is $6$).
Daniele de Grazia's avatar
Daniele de Grazia committed
288 289 290 291 292 293 294
\begin{lstlisting}[style=XMLStyle]
<EXPANSIONS>
    <E COMPOSITE="C[0]" NUMMODES="7" FIELDS="u" TYPE="MODIFIED" />
    <E COMPOSITE="C[1]" NUMMODES="7" FIELDS="u" TYPE="MODIFIED" />
</EXPANSIONS>
\end{lstlisting}

295 296 297 298
Only one parameter is needed for this problem. In this example $\lambda = 1$ and
the Continuous Galerkin Method is used as projection scheme to solve the 
Helmholtz equation, so we need to specify the following parameters and solver
information.
Daniele de Grazia's avatar
Daniele de Grazia committed
299
\begin{lstlisting}[style=XMLStyle]
300 301 302
<PARAMETERS>
    <P> Lambda = 1 </P>
</PARAMETERS>
Daniele de Grazia's avatar
Daniele de Grazia committed
303

304 305 306 307
<SOLVERINFO>
    <I PROPERTY="EQTYPE" VALUE="Helmholtz" />
    <I PROPERTY="Projection" VALUE="Continuous" />
</SOLVERINFO>
Daniele de Grazia's avatar
Daniele de Grazia committed
308 309
\end{lstlisting}

310 311 312 313
All three basic boundary condition types have been used in this example:
Dirichlet, Neumann and Robin boundary. The boundary regions are defined, each of
which corresponds to one of the edge composites defined earlier. Each boundary
region is then assigned an appropriate boundary condition.
Daniele de Grazia's avatar
Daniele de Grazia committed
314 315
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYREGIONS>
316 317 318 319 320 321 322 323 324 325 326
    <B ID="0"> C[2] </B>
    .
    .
    <B ID="8"> C[10] </B>
</BOUNDARYREGIONS>

<BOUNDARYCONDITIONS>
    <REGION REF="0">
        <D VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
    </REGION>
    <REGION REF="1">
327 328
        <R VAR="u" VALUE="sin(PI*x)*sin(PI*y)-PI*sin(PI*x)*cos(PI*y)" 
           PRIMCOEFF="1" />
329 330
    </REGION>
    <REGION REF="2">
331 332
        <N VAR="u" VALUE="(5/sqrt(61))*PI*cos(PI*x)*sin(PI*y)-
                          (6/sqrt(61))*PI*sin(PI*x)*cos(PI*y)" />
333
    </REGION>
Daniele de Grazia's avatar
Daniele de Grazia committed
334 335
        .
        .
336
</BOUNDARYCONDITIONS>
Daniele de Grazia's avatar
Daniele de Grazia committed
337 338
\end{lstlisting}

339 340 341 342
We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact 
solution of the two-dimensional Helmholtz equation is $u = sin(\pi x)cos(\pi
y)$. These functions are defined specified to initialise the problem and verify
the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors.
Daniele de Grazia's avatar
Daniele de Grazia committed
343 344
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="Forcing">
345 346
    <E VAR="u" VALUE="-(Lambda + 2*PI*PI)*sin(PI*x)*sin(PI*y)" />
</FUNCTION>
Daniele de Grazia's avatar
Daniele de Grazia committed
347

348 349 350
<FUNCTION NAME="ExactSolution">
    <E VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
</FUNCTION>
Daniele de Grazia's avatar
Daniele de Grazia committed
351 352 353
\end{lstlisting}


354 355
\subsubsection{Running the code}
\begin{lstlisting}[style=BashInputStyle]
356
ADRSolver Test_Helmholtz2D_modal.xml
357
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
358

359 360
This execution should print out a summary of input file, the $L_2$ and 
$L_{inf}$ errors and the time spent on the calculation.
Daniele de Grazia's avatar
Daniele de Grazia committed
361

362 363 364 365
\subsubsection{Post-processing}
Simulation results are written in the file Helmholtz2D\_modal.fld. We can choose
to visualise the output in Gmsh
\begin{lstlisting}[style=BashInputStyle]
366
FldToGmsh Helmholtz2D_modal.xml Helmholtz2D_modal.fld
367 368 369
\end{lstlisting}
which generates the file Helmholtz2D\_modal\_u.pos as shown in
Fig.~\ref{f:adrsolver:helmholtz2D}
Daniele de Grazia's avatar
Daniele de Grazia committed
370

371
\begin{figure}
Daniele de Grazia's avatar
Daniele de Grazia committed
372
\begin{center}
373
\includegraphics[width=6cm]{img/Helmholtz2D}
Daniele de Grazia's avatar
Daniele de Grazia committed
374
\caption{Solution of the 2D Helmholtz Problem.}
375
\label{f:adrsolver:helmholtz2D}
Daniele de Grazia's avatar
Daniele de Grazia committed
376 377 378 379
\end{center}
\end{figure}


Daniele de Grazia's avatar
Daniele de Grazia committed
380 381
\subsection{Advection dominated mass transport in a pipe}

382 383 384 385 386 387 388 389 390 391 392
The following example demonstrates the application of the ADRsolver for
modelling advection dominated mass transport in a straight pipe.
Such a transport regime is encountered frequently when modelling mass transport
in arteries. This is because the diffusion coefficient of small blood borne
molecules, for example oxygen or adenosine triphosphate, is very small
$O(10^{-10})$.

\subsubsection{Background}
The governing equation for modelling mass transport is the unsteady advection
diffusion equation:
\begin{align*}
Daniele de Grazia's avatar
Daniele de Grazia committed
393
\dfrac{\partial u}{\partial t}  + v\nabla u +  \epsilon \nabla^2 u = 0
394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412
\end{align*}

For small diffusion coefficient, $\epsilon$, the transport is dominated by
advection and this leads to a very fine boundary layer adjacent to the surface
which must be captured in order to get a realistic representation of the wall
mass transfer processes. This creates problems not only from a meshing
perspective, but also numerically where classical oscillations are observed in
the solution due to under-resolution of the boundary layer.

The Graetz-Nusselt solution is an analytical solution of a developing mass (or
heat) transfer boundary layer in a pipe. Previously this solution has been used
as a benchmark for the accuracy of numerical methods to capture the fine
boundary layer which develops for high Peclet number transport (the ratio of
advection to diffusion). The solution is derived based on the assumption that
the velocity field within the mass transfer boundary layer is linear i.e. the
Schmidt number (the relative thickness of the momentum to mass transfer boundary
layer) is sufficiently large. The analytical solution for the non-dimensional
mass transfer at the wall is given by:
\begin{align*}
Daniele de Grazia's avatar
Daniele de Grazia committed
413
S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} , 
414
\end{align*}
Daniele de Grazia's avatar
Daniele de Grazia committed
415 416
where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete 
Gamma function and $Pe$ the Peclet number given by:
417
\begin{align*}
Daniele de Grazia's avatar
Daniele de Grazia committed
418
Pe = \dfrac{2 U R}{\epsilon}
419
\end{align*}
Daniele de Grazia's avatar
Daniele de Grazia committed
420

421 422 423 424
In the following we will numerically solver mass transport in a pipe and compare
the calculated mass transfer at the wall with the Graetz-Nusselt solution. The
Peclet number of the transport regime under consideration is $750000$, which is
physiologically relevant.
Daniele de Grazia's avatar
Daniele de Grazia committed
425

426 427 428
\subsubsection{Input file}
The geometry under consideration is a pipe of radius, $R = 0.5$ and length $l =
0.5$
Daniele de Grazia's avatar
Daniele de Grazia committed
429 430 431

\begin{figure}[h!]
\begin{center}
432
\includegraphics[width=6cm]{img/pipe}
Daniele de Grazia's avatar
Daniele de Grazia committed
433 434 435 436
\caption{Pipe.}
\end{center}
\end{figure}

437 438 439 440
Since the mass transport boundary layer will be confined to a very small layer
adjacent to the wall we do not need to mesh the interior region, hence the mesh
consists of a layer of ten prismatic elements over a thickness of 0.036R. The
elements progressively grow over the thickness of domain.
Daniele de Grazia's avatar
Daniele de Grazia committed
441

442 443 444 445 446
In this example we utilise heterogeneous polynomial order, in which the
polynomial order normal to the wall is higher so that we avoid unphysical
oscillations, and hence the incorrect solution, in the mass transport boundary
layer. To do this we specify explicitly the expansion type, points type and
distribution in each direction as follows:
Daniele de Grazia's avatar
Daniele de Grazia committed
447 448 449 450 451 452 453 454 455 456 457
\begin{lstlisting}[style=XMLStyle]
<EXPANSIONS>
  <E COMPOSITE="C[0]"
     NUMMODES="3,5,3"
     BASISTYPE="Modified_A,Modified_A,Modified_B"
     NUMPOINTS="4,6,3"
     POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha1Beta0"
     FIELDS="u" />
</EXPANSIONS>
\end{lstlisting}

458 459 460
The above represents a quadratic polynomial order in the azimuthal and
streamwise direction and 4th order polynomial normal to the wall for a prismatic
element.
Daniele de Grazia's avatar
Daniele de Grazia committed
461

462
We choose to use a continuous projection and an first-order implicit-explicit
463 464
time-integration scheme. The \inltt{DiffusionAdvancement} and
\inltt{AdvectionAdvancement} parameters specify how these terms are treated.
Daniele de Grazia's avatar
Daniele de Grazia committed
465
\begin{lstlisting}[style=XMLStyle]
466 467 468 469 470 471
<I PROPERTY="EQTYPE"                VALUE="UnsteadyAdvectionDiffusion" />
<I PROPERTY="Projection"            VALUE="Continuous" />
<I PROPERTY="DiffusionAdvancement"  VALUE="Implicit" />
<I PROPERTY="AdvectionAdvancement"  VALUE="Explicit" />
<I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
<I PROPERTY="GlobalSysSoln"         VALUE="IterativeStaticCond" />
Daniele de Grazia's avatar
Daniele de Grazia committed
472 473
\end{lstlisting}

474
We integrate for a total of $30$ time units with a time-step of $0.0005$,
475
necessary to keep the simulation numerically stable.
Daniele de Grazia's avatar
Daniele de Grazia committed
476
\begin{lstlisting}[style=XMLStyle]
477 478 479
<P> TimeStep = 0.0005              </P>
<P> FinalTime = 30                 </P>
<P> NumSteps = FinalTime/TimeStep  </P>
Daniele de Grazia's avatar
Daniele de Grazia committed
480 481
\end{lstlisting}

482
The value of the $\epsilon$ parameter is $\epsilon = 1/Pe$
483 484 485
\begin{lstlisting}[style=XMLStyle]
<P> epsilon = 1.33333e-6           </P>
\end{lstlisting}
Daniele de Grazia's avatar
Daniele de Grazia committed
486

487 488 489 490 491 492 493 494 495 496
The analytical solution represents a developing mass transfer boundary layer in
a pipe. In order to reproduce this numerically we assume that the inlet
concentration is a uniform value and the outer wall concentration is zero; this
will lead to the development of the mass transport boundary layer along the
length of the pipe. Since we do not model explicitly the mass transfer in the
interior region of the pipe we assume that the inner wall surface concentration
is the same as the inlet concentration; this assumption is valid based on the
large Peclet number meaning the concentration boundary layer is confined to the
region in the immediate vicinity of the wall. The boundary conditions are
specified as follows in the input file:
Daniele de Grazia's avatar
Daniele de Grazia committed
497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYREGIONS>
   <B ID="0"> C[3] </B>  <!-- inlet -->
   <B ID="1"> C[4] </B>  <!-- outlet -->
   <B ID="2"> C[2] </B>  <!-- outer surface -->
   <B ID="3"> C[5] </B>  <!-- inner surface -->
</BOUNDARYREGIONS>

<BOUNDARYCONDITIONS>
  <REGION REF="0">
    <D VAR="u" VALUE="1" />
  </REGION>
  <REGION REF="1">
    <N VAR="u" VALUE="0" />
  </REGION>
  <REGION REF="2">
    <D VAR="u" VALUE="0" />
  </REGION>
  <REGION REF="3">
    <D VAR="u" VALUE="1" />
  </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}

521 522
The velocity field within the domain is fully devqeloped pipe flow (Poiseuille
flow), hence we can define this through an analytical function as follows:
Daniele de Grazia's avatar
Daniele de Grazia committed
523 524 525 526 527 528 529 530
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="AdvectionVelocity">
  <E VAR="Vx" VALUE="0" />
  <E VAR="Vy" VALUE="0" />
  <E VAR="Vz" VALUE="2.0*(1-(x*x+y*y)/0.25)" />
</FUNCTION>
\end{lstlisting}

531 532
We assume that the initial domain concentration is uniform everywhere and the
same as the inlet. This is defined by, 
Daniele de Grazia's avatar
Daniele de Grazia committed
533 534 535 536 537 538
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="InitialConditions">
  <E VAR="u" VALUE="1" />
</FUNCTION>
\end{lstlisting}

539 540 541 542 543
\subsubsection{Results}
To compare with the analytical expression we numerically calculate the
concentration gradient at the surface of the pipe. This is then plotted against
the analytical solution by extracting the solution along a line in the
streamwise direction, as shown in Fig.~\ref{f:adrsolver:masstransport}.
Daniele de Grazia's avatar
Daniele de Grazia committed
544 545 546

\begin{figure}[h!]
\begin{center}
547
\includegraphics[width=7cm]{img/graetz-nusselt}
Daniele de Grazia's avatar
Daniele de Grazia committed
548
\caption{Concentration gradient at the surface of the pipe.}
549
\label{f:adrsolver:masstransport}
Daniele de Grazia's avatar
Daniele de Grazia committed
550 551 552 553
\end{center}
\end{figure}


Daniele de Grazia's avatar
Daniele de Grazia committed
554