User Guide PDF
User Guide PDF
User Guide PDF
Version 3.0
Acknowledgments
The work leading to the preparation of this document has received funding from
the European Research Council under the European Union’s Seventh Framework
Programme (FP7/2007-2013)/ERC Grant agreement no 307499. The
collaboration with Professor Fernando T. Pinho (University of Porto, Portugal),
Professor Paulo J. Oliveira (University of Beira Interior, Portugal) and Dr
Alexandre Afonso (University of Porto, Portugal) in the development of
numerical methods for computational rheology is also acknowledged.
Disclaimer
This offering is not approved nor endorsed by ESI-Group, the producer of the
OpenFOAM
R
software and owner of the OpenFOAM
R
trademark.
The recommendations expressed in this document are those of the authors and
are not necessarily the views of, or endorsement by, third parties named in this
document.
RheoTool, where this guide is included, is distributed in the hope that it will be
useful, but WITHOUT ANY WARRANTY. See the GNU General Public
License (http://www.gnu.org/licenses/) for more details.
Trademarks
Linux is a registered trademark of Linus Torvalds.
OpenFOAM is a registered trademark of ESI Group.
Paraview is a registered trademark of Kitware.
Typeset in LATEX.
c 2016-2018 Francisco Pimenta, Manuel A. Alves
Contents
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Guide organization . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Changelog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Citing rheoTool . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.5 Contacts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.6 Contributing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Installation 7
2.1 Folder organization . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Compatibility with OpenFOAM
R
and foam-extend versions . . . . 8
2.3 System requirements . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4 Downloading Eigen library . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Installing rheoTool . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.6 Differences between versions . . . . . . . . . . . . . . . . . . . . . . 10
3 Theoretical background 12
3.1 Governing equations of complex fluid flows . . . . . . . . . . . . . . 12
3.2 Stabilization of viscoelastic fluid flow simulations . . . . . . . . . . 13
3.2.1 The both-sides-diffusion (BSD) technique . . . . . . . . . . . 13
3.2.2 The log-conformation tensor approach . . . . . . . . . . . . 14
3.3 Coupling algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3.1 Pressure-velocity coupling . . . . . . . . . . . . . . . . . . . 15
3.3.2 Stress-velocity coupling . . . . . . . . . . . . . . . . . . . . . 16
3.4 High-resolution schemes . . . . . . . . . . . . . . . . . . . . . . . . 17
3.5 Moving grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.6 Electrically-driven flow models . . . . . . . . . . . . . . . . . . . . . 19
3.6.1 Poisson-Nernst-Planck model . . . . . . . . . . . . . . . . . 19
3.6.2 Splitting the electric potential . . . . . . . . . . . . . . . . . 20
3.6.3 Poisson-Boltzmann model . . . . . . . . . . . . . . . . . . . 21
3.6.4 Debye-Hückel model . . . . . . . . . . . . . . . . . . . . . . 21
3.6.5 Slip model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.6.6 Ohmic (leaky dielectric) model . . . . . . . . . . . . . . . . 23
3.7 Brownian dynamics simulations . . . . . . . . . . . . . . . . . . . . 24
3.7.1 The bead-spring model . . . . . . . . . . . . . . . . . . . . . 24
3.7.2 Governing equations of beads motion . . . . . . . . . . . . . 27
ii
CONTENTS iii
4 Overview of rheoTool 30
4.1 The constitutiveEquations library . . . . . . . . . . . . . . . . . . . 30
4.1.1 Available GNF and viscoelastic models . . . . . . . . . . . . 30
4.1.2 A note on FENE-type models . . . . . . . . . . . . . . . . . 38
4.1.3 Multi-mode modeling . . . . . . . . . . . . . . . . . . . . . . 40
4.1.4 Analysis of a code sample . . . . . . . . . . . . . . . . . . . 40
4.1.5 Advanced settings . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1.6 Adding new viscoelastic or GNF models . . . . . . . . . . . 46
4.2 The EDFModels library . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.1 Available EDF models . . . . . . . . . . . . . . . . . . . . . 47
4.2.2 The potentials splitting approach and multi-species model-
ing in the PNP, PB and DH models . . . . . . . . . . . . . . 49
4.2.3 Electrokinetic coupling loop in the PNP model . . . . . . . . 49
4.2.4 Analysis of a code sample . . . . . . . . . . . . . . . . . . . 49
4.2.5 Adding new EDF models . . . . . . . . . . . . . . . . . . . . 56
4.3 The BDmolecule library . . . . . . . . . . . . . . . . . . . . . . . . 57
4.3.1 Organization of variables . . . . . . . . . . . . . . . . . . . . 57
4.3.2 Solution sequence . . . . . . . . . . . . . . . . . . . . . . . . 57
4.3.3 External forcing type . . . . . . . . . . . . . . . . . . . . . . 61
4.3.4 External forcing interpolation . . . . . . . . . . . . . . . . . 61
4.3.5 Spring force and time-integration schemes . . . . . . . . . . 64
4.3.6 Tethering and fixing the molecules center of mass . . . . . . 66
4.3.7 Beads tracking . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.3.8 Data output for post-processing . . . . . . . . . . . . . . . . 67
4.3.9 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.4 Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.4.1 rheoFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.4.2 rheoTestFoam . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.4.3 rheoInterFoam . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.4.4 rheoEFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.4.5 rheoBDFoam . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.5 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.5.1 linearExtrapolation . . . . . . . . . . . . . . . . . . . . . . . 83
4.5.2 navierSlip . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.5.3 zeroIonicFlux . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.5.4 boltzmannEquilibrium . . . . . . . . . . . . . . . . . . . . . . 85
4.5.5 inducedPotential . . . . . . . . . . . . . . . . . . . . . . . . 85
4.5.6 slipSmoluchowski . . . . . . . . . . . . . . . . . . . . . . . . 85
4.5.7 slipSigmaDependent . . . . . . . . . . . . . . . . . . . . . . 86
4.5.8 A note on wall boundary conditions for pressure . . . . . . . 86
4.6 Utilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.6.1 GaussDefCmpw schemes for convective terms . . . . . . . . 87
4.6.2 Generic post-processing: ppUtil . . . . . . . . . . . . . . . . 89
CONTENTS iv
4.6.3 writeEfield . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.6.4 initMolecules . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.6.5 averageMolcN . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.6.6 averageMolcX . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5 Tutorials 99
5.1 rheoFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5.1.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 99
5.1.2 A note on coded FunctionObjects . . . . . . . . . . . . . . . 105
5.1.3 Case 1: flow between parallel plates . . . . . . . . . . . . . . 106
5.1.4 Case 2: lid-driven cavity flow . . . . . . . . . . . . . . . . . 108
5.1.5 Case 3: flow in a 4:1 planar contraction . . . . . . . . . . . . 109
5.1.6 Case 4: flow around a confined cylinder . . . . . . . . . . . . 112
5.1.7 Case 5: bifurcation in a 2D cross-slot flow . . . . . . . . . . 115
5.1.8 Case 6: blood flow simulation in a real-model aneurysm . . . 117
5.1.9 Case 7: viscous fluid damper (moving mesh) . . . . . . . . . 121
5.2 rheoTestFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.2.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 124
5.2.2 Case I: Herschel-Bulkley model . . . . . . . . . . . . . . . . 126
5.2.3 Case II: FENE-CR model . . . . . . . . . . . . . . . . . . . 127
5.3 rheoInterFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
5.3.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 130
5.3.2 Case 1: impacting drop . . . . . . . . . . . . . . . . . . . . . 131
5.3.3 Case 2: planar die swell . . . . . . . . . . . . . . . . . . . . 133
5.4 rheoEFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
5.4.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 135
5.4.2 Case I: EDF of power-law and PTT fluids in a microchannel 138
5.4.3 Case II: induced-charge electroosmosis around a cylinder . . 142
5.4.4 Case III: charge transport across an ion-selective membrane 144
5.4.5 Case IV: electrokinetic instabilities in a flow-focusing device 146
5.4.6 Case V: electrokinetic mixer . . . . . . . . . . . . . . . . . . 150
5.4.7 Case VI: electro-elastic instabilities in cross-shaped geometries152
5.5 rheoBDFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
5.5.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 154
5.5.2 Molecules visualization with Paraview . . . . . . . . . . . . 159
5.5.3 Case 1: λ-DNA extension in a planar extensional flow . . . . 160
5.5.4 Case 2: 7λ-DNA extension in a flow-focusing device . . . . . 163
5.5.5 Case 3: λ-DNA dynamics in LAOE . . . . . . . . . . . . . . 166
Bibliography 174
Chapter 1
Introduction
1.1 Motivation
The open-source OpenFOAM
R
toolbox can be used as a versatile finite-volume
solver for CFD simulations in general polyhedral grids. A number of constitutive
equations for Generalized Newtonian Fluids (GNF) are already available in the
toolbox for a long time. More recently, Favero et al. [1] created a library containing
a wide range of constitutive equations to model viscoelastic fluids, along with a
solver named viscoelasticFluidFoam which makes use of this library. However,
viscoelasticFluidFoam presents stability issues in certain conditions, such as, for
example, in the simulation of high Weissenberg number (Wi) flows or when there
is no solvent viscosity contribution (e.g. in the upper-convected Maxwell model).
In Ref. [2], we attempted to minimize those issues by modifying critical points
in the viscoelasticFluidFoam solver and in the handling of viscoelastic models. The
modified solver was tested in benchmark flows and second-order accuracy, both in
space and time, was observed, in addition to an enhanced stability [2]. The package
that we present in this document – rheoTool – implements the method described
in [2].
Afterwards, the capability to simulate electrically-driven flows was added to
rheoTool [3] and is available since version 2.0.
Recognizing the importance of modeling polymeric flows at different scales, a
Brownian dynamics solver has been implemented in rheoTool [4], which is available
since version 3.0.
rheoTool is more than a collection of solvers and libraries. In addition to robust
solvers for the simulation of pressure- and electrically-driven flows of both GNF
and viscoelastic fluids, we provide also tutorials and utilities that can be useful for
the users starting to apply the OpenFOAM
R
toolbox in the simulation of complex
fluid flows. In particular, some of the distinguishing features of rheoTool are:
• both GNF and viscoelastic models can be selected on run time and applied to
single-phase laminar flows. A solver for two-phase flows is also being devel-
oped and an experimental (but fully functional) version is already available.
1
CHAPTER 1. Introduction 2
• transient flow solvers use the SIMPLEC algorithm for pressure-velocity cou-
pling, instead of the PISO implementation. Large time-steps can be used
without decoupling problems, and the use of under-relaxation is not required
(except for pressure in some problems using non-orthogonal grids).
• the tool is provided with a user-guide (this document) and a selected set of
tutorials reproducing relevant benchmark or real-life flow problems.
• Chapter 5 contains several tutorials, guiding the reader into the use of
rheoTool .
1
http://openfoam.org/
2
http://www.extend-project.de/
CHAPTER 1. Introduction 3
The language and the content used in this guide assumes that the reader has a
basic knowledge on the use of the OpenFOAM
R
toolbox and is familiar with the
finite-volume method applied to CFD problems. Thus, it is out the scope of this
document to serve as an introduction on those subjects.
Although rheoTool is available for different OpenFOAM
R
and foam-extend
versions, for historical reasons Chapters 4 and 5 are still mainly based on
OpenFOAM
R
version 2.2.2 to describe the contents (except the content related
with Brownian dynamics simulations, described for OpenFOAM
R
version 5.x).
However, the small differences among different versions should not be an obstacle
to the readers using any other version.
The readers interested in the theory behind rheoTool are strongly encouraged
to first read Refs. [2], [3] and [4] before this guide.
1.3 Changelog
Version 3.0
Released on 18/09/2018.
• Add: solvers, libraries, utilities and tutorials for Brownian dynamics simu-
lations of polymer molecules.
Generic
• Add: all solvers are now compatible with dynamic meshes. Due to this
change, and for convenience, momentum equation is the first to be solved,
followed by pressure equation and then the equations for the remaining vari-
ables (extra-stresses, passive scalar, etc.).
• Add: tutorial fluidDamper showing the use of rheoFoam with dynamic
meshes.
• Add: added an explicit Navier slip boundary condition for velocity.
• Change: Namespace encapsulation of several derived classes.
• Change: rheoTool version compatible with OpenFOAM
R
v2.2.2 is discon-
tinued.
• Add: added rheoTool patch for OpenFOAM
R
v6.0.
• Add: added note in the user-guide (Section 2.5) about parallel compilation
of rheoTool .
Constitutive equations
• Add: Papanastasiou regularization is now available for yield-stress GNF
models (Hershel-Bulkley/Bingham and Casson models).
CHAPTER 1. Introduction 4
• Add: Casson model has been added to the library of constitutive equations.
• Add: the Vasquez-Cook-Mckinley (VCM) model has been added to the li-
brary of constitutive equations.
• Add: the Reactive Rod Model (RRM) model has been added to the library
of constitutive equations.
• Add: the Bautista-Manero-Puig (BMP) model has been added to the library
of constitutive equations. Both stress and log-conformation versions are
available.
Version 2.0
Released on 09/02/2018.
Electrically-driven flows
Constitutive equations
• Add: the Rolie-Poly viscoelastic model has been added to the library of
constitutive equations. Both the stress and log-conformation versions are
available.
• Change: sPTT models have been generalized to their full form by replacing
the upper-convected derivative by the Gordon-Schowalter derivative. It is
now possible to simulate PTT models with non-affine deformation, in both
the stress and log-conformation versions.
Post-Processing
• Add: class ppUtil for post-processing purposes has been added to the versions
for OpenFOAM
R
and the one existing for foam-extend has been modified.
Enable the use of multiple ppUtil in simultaneous.
• Fix: sampling error was fixed for the tutorials of versions of40 and fe40.
Multiphase flows
Generic
• Change/Fix: code cleanup and bug fix (BC evaluation of the explicit fvc::
div(phi,X) operator) in class GaussDefCmpw .
Version 1.0
Released on 6/12/2016.
Initial version.
@misc{rheoTool,
author = "F. Pimenta and M.A. Alves",
title = "rheoTool",
howpublished = "\url{https://github.com/fppimenta/rheoTool}",
year = "2016"}
Since the underlying theory of rheoTool has been mainly presented in technical
papers (Refs. [2], [3] and [4]), these can also be used for citation purposes.
1.5 Contacts
rheoTool is under continuous development and new features and improvements
will be added in the future. If you have any suggestions, comments or doubts
regarding the tool, or if you found a bug or error, feel free to contact us:
1.6 Contributing
In the open-source spirit, rheoTool is open to contributions from the community.
If you believe that your piece of code is worth to be incorporated in rheoTool ’s
next version, feel free to contact us.
Chapter 2
Installation
doc
(…)
(…)
(…) (…)
feyy solvers
tutorials
The top-level directory of rheoTool contains the versions available for different
OpenFOAM
R
(of ) and foam-extend (fe) versions (see next section for compati-
bility issues). The folder doc/, containing the user guide, is also in the top-level
directory. Inside the folder for each version, there are two directories: src/, where
the source code can be found, and tutorials/, containing several tutorial cases
showing the use of rheoTool . The src/ directory is further subdivided in a di-
rectory with the applications (solvers/) and another one containing libraries
(libs/).
After cloning/downloading rheoTool , the user is free to remove from the top-
level directory all the versions not needed and keep only the one(s) of interest.
7
CHAPTER 2. Installation 8
version 4.0. However, an effort has been made to release rheoTool also running
under other (more recent) versions of OpenFOAM
R
and foam-extend. Thus, com-
patible versions of rheoTool are provided for:
• OpenFOAM
R
v4.0, v4.1 (of40/).
• OpenFOAM
R
v6.0 (of60/).
Note that the list above includes the versions which were effectively tested.
This means that a given version of rheoTool may be compatible with other
OpenFOAM
R
or foam-extend versions not included in this list. The versions
above were tested in a Ubuntu 16.04 environment, but other operating systems
running OpenFOAM
R
can eventually support some version of rheoTool . However,
the installation is only described here for a Linux OS.
After ensuring that these conditions are fulfilled, the user is ready to start the
installation, which includes two major steps: downloading (no install) the open-
source Eigen library [7] and installing rheoTool .
or foam-extend. If this is the case, be sure that the alias pointing to the desired
version has been typed. Shortly, you should only advance to the next step if a
command like ∼$ icoFoam -help is recognized in the terminal. Note that in
this document we use the prepending ∼$ for any instruction to be typed in the
command line (thus, ∼$ icoFoam -help means that you only type icoFoam
-help). If this check is successful, run the script downloadEigen in that terminal:
CHAPTER 2. Installation 9
∼$ ./downloadEigen
This script downloads Eigen version 3.2.9 (other versions close to that would
also work adequately) from the Internet (using wget), extracts it and moves it to
directory:
$WM PROJECT USER DIR/ThirdParty/Eigen3.2.9
Eigen is used in rheoTool for computation of eigenvalues and eigenvectors and
there is no need to install the library, since the inclusion of the required headers
is enough for our purposes.
However, its location in the system must defined and exported. This is achieved
by attributing to variable EIGEN_RHEO – the one used and recognized by rheoTool –
the actual path of Eigen. The command to do so has been displayed to the terminal
after running script downloadEigen (if everything was ok) and looks like:
∼$ echo "export EIGEN_RHEO=/home/user/OpenFOAM/user-4.0/ThirdParty
/Eigen3.2.9">>/home/user/.bashrc
Do not copy this command, it is just an example of what is displayed to the
screen. Instead, copy-paste and run the command appearing in your terminal.
If, for some reason, the user wants to move Eigen to another directory (or
already has an Eigen version in another directory), then move Eigen to its final
location (if already not) and define variable EIGEN_RHEO accordingly. Note that
Eigen only needs to be installed once per system. Even if the user has
installed multiple versions of rheoTool in the same system, the above procedure
only needs to be run once (for the first version being installed), as long as the
directory containing Eigen since the first installation is not deleted, or renamed.
∼$ ./Allwmake -j N
for parallel compilation with N processors. For example, ./Allwmake -j 3
will compile in parallel using 3 processors. If the number of processors is not
specified, all available processors are used. If option −j is not passed to the script,
the compilation will use WM NCOMPPROCS processors, where this variable is
usually defined in the etc/bashrc file of your OpenFOAM
R
installation. For
foam-extend users, option −j is not recognized by the script, therefore simply run
∼$ ./Allwmake
The compilation in foam-extend will typically make use of all processors avail-
able in your system, since variable WM NCOMPPROCS is set by default in such
way.
Both the libraries and applications installed with rheoTool can be cleaned by
running the script Allwclean.
Since the user will probably not need the remaining versions of rheoTool that
remain in the top-level directory, they can simply be deleted, if already not.
To check if the installation succeeded, the user should try to run one of the
tutorials of Chapter 5.
Theoretical background
∇· u = 0 (3.1)
∂u 0
ρ + u· ∇u = −∇p + ∇· τ + f (3.2)
∂t
0
where u is the velocity vector, t is the time, p is the pressure, τ is the extra-stress
tensor and f is any external body-force, such as the electric force discussed in
Section 3.6. To simulate viscoelastic fluid flows, it is a common approach to split
the total extra-stress tensor in a solvent contribution (τs ) and a polymeric contri-
0
bution (τ), τ = τ + τs . In order to have a closed set of equations, a constitutive
equation is required for each tensor contribution, which can be generally written
as in Eqs. (3.3) and (3.4), for a wide range of models,
∇
f (τ)τ + λ(γ̇) τ + h(τ) = ηp (γ̇)(∇u + ∇uT ) (3.4)
12
CHAPTER 3. Theoretical background 13
In Eqs. (3.3) and (3.4), ηs is the solvent viscosity, ηp is the polymeric viscosity
coefficient, λ is the relaxation time, γ̇ is the shear-rate, f (τ) is a general scalar
function depending on an invariant of τ, h(τ) is a tensor-valued function depend-
∇
ing on τ and τ= ∂τ ∂t
+ u· ∇τ − τ· ∇u − ∇uT · τ represents the upper-convected
time derivative, which renders the models frame-invariant. Some models use the
∇
Gordon-Schowalter derivative (τ =τ +ζ (τ· D + D· τ), with D = 21 (∇u + ∇uT ))
instead of the upper-convected derivative, in order to take non-affine deformation
into account (controlled by parameter ζ). In rheoTool , this is the case of PTT-
type models. Other constitutive models exist, which can also make use of the
lower-convected time derivative, but those are not explored here. The constitu-
tive equation for a GNF is limited to Eq. (3.3), since elasticity is not considered
0
(τ = τs ). In Table 4.1 presented in the next Chapter, Eqs. (3.3) and (3.4) are
specified for several GNF and viscoelastic models.
Eqs. (3.1)–(3.4) represent the standard system of equations to be solved. How-
ever, due to numerical stability issues in viscoelastic fluid flow simulations, the
system is rarely solved in that form. Indeed, several techniques are available for
stabilization purposes (see, for instance, Ref. [8] for a comparison between the
most popular techniques) and the ones used in rheoTool are addressed next.
∂u
ρ + u· ∇u − ∇· (ηs + ηp )∇u = −∇p − ∇· (ηp ∇u) + ∇· τ + f (3.5)
∂t
Note that the added diffusive terms are scaled by the polymeric viscosity (ηp ),
which is a common choice in the literature (e.g. Ref. [8]), although not mandatory.
In order to simplify the reading, the possible dependence of the viscosity and
relaxation time on the shear-rate will be dropped in the respective symbols, as
already done in Eq. (3.5), although this relation still holds to keep generality.
CHAPTER 3. Theoretical background 14
Θ = RΛΘ RT (3.13)
and the conformation tensor is recovered by the inverse relation of Eq. (3.7)
Finally, the polymeric extra-stress tensor can be computed from A (Eq. 3.6)
and used in the momentum equation.
Note that for PTT-type models, which may include non-affine deformation
through the Gordon-Schowalter derivative, the tensor M (Eq. 3.11) is computed
differently: M = R ∇uT − ζD RT .
It is worth to mention that the log-conformation approach can be considered a
particular case of the kernel-conformation method [10]. However, from our expe-
rience, the log kernel is frequently the optimal kernel (in terms of robustness and
accuracy) for generic problems, so that only this one is widely used in rheoTool .
Nevertheless, for the Oldroyd-B model, the rootk kernel [10] and the square-root
transformation [11] are also included in rheoTool for demonstration purposes.
1 H 1 1 ∗
∇· (∇p)P = ∇· + − (∇p )P (3.15)
aP − H1 aP aP − H1 aP
P
where aP are the diagonal coefficients from the momentum equation, H1 = − anb
nb
is an operator representing the negative sum of the off-diagonal coefficients from
momentum equation, H = − anb u∗nb + b is an operator containing the off-
P
nb
diagonal contributions, plus source terms (except the pressure gradient) of the
momentum equation and p∗ is the pressure field known from the previous time-
step or iteration. Accordingly, the equation to correct the velocity after obtaining
the continuity-compliant pressure field from Eq. (3.15) is
H 1 1 1
u= + − (∇p∗ )P − (∇p)P (3.16)
aP aP − H1 aP aP − H1
Importantly, in order to avoid the onset of checkerboard fields, the pressure
gradient terms involved in the computation of face velocities, i.e., in Eqs. (3.15)
CHAPTER 3. Theoretical background 16
and (3.16), are directly evaluated using the pressure on the cells straddling the face,
in a Rhie-Chow-like procedure (more details in Ref. [2]). Nonetheless, when Eq.
(3.16) is used to correct the cell-centered velocity field, the pressure gradient terms
are computed ”in the usual way”, for example using Green-Gauss integration.
Rhie-Chow methods used to avoid checkerboard fields, as the one described
in the previous paragraph, are known to be affected by the use of small time-
steps and they also present time-step dependency on steady-state results [12].
In OpenFOAM
R
solvers, a common strategy to avoid such effects is to add a
corrective term to face-interpolated velocities, through functions ddtPhiCorr() or
ddtCorr(). Recently, in foam-extend the time-step dependency was solved in a
different way, by removing the transient term contribution from the aP coefficients
of the momentum equation [13]. However, this approach may be problematic when
used with the SIMPLEC algorithm, since a division by zero is prone to happen.
In rheoTool , we keep using the added corrective term, although, as mentioned in
Ref. [2], this term can be improved in order to more efficiently avoid the small
time-step dependency of steady-state solutions.
∂u
ρ + u· ∇u − ∇· (ηs + ηp )∇u = −∇p − ∇· ηp ∇u + ∇· τ + f (3.18)
∂t
in the momentum equation (resulting in Eq. 3.18), we drop the transpose velocity
gradients for simplicity, since continuity imposes ∇· ∇uT = 0.
φf = [φC ]implicit + [(α − 1)φC + βφD + (1 − α − β)(φD − 2(∇φ)C · dCD )]explicit (3.23)
HRS. Additionally, in Ref. [2] it was observed that the usual methodology of
OpenFOAM
R
to apply HRSs to non-scalar variables (tensors and vectors) can
locally introduce numerical instabilities in some viscoelastic flow problems. This
methodology consists in using a frame-invariant quantity for non-scalar variables,
such as the squared magnitude for vectors, or the trace (or double-dot product)
for tensors, to compute the α and β parameters in Eq. (3.23). It was observed
that such artificial instabilities can be significantly damped with a component-wise
handling of non-scalar variables [2], at the cost of losing frame-invariance, which
however is very weak and vanishes with grid refinement. Accordingly, non-scalar
variables are split into its components and Eq. (3.23) is applied independently to
each one of them. Note that this approach still generates one single matrix of co-
efficients for such variables, since the upwind differencing scheme coefficients are
common to all the components (they only depend on the flux). The differentiation
between components is only introduced in the explicit part of Eq. (3.23), generat-
ing a different source term for each individual tensor/vector component. This is
possible due to the use of a deferred correction approach.
where φ is any generic variable being advected and ub is the velocity at which
surface S is moving. Moreover, the space conservation law (SCL) needs to be
satisfied to ensure mass conservation [16],
Z Z
d
dV − ub · ndS = 0 (3.25)
dt V S
CHAPTER 3. Theoretical background 19
If the SCL is ensured, the continuity equation remains unchanged and R so does
the pressure equation. In practice, the SCL is imposed while computing S ub ·ndS
in Eq. (3.24), which is the flux due to mesh motion. According to Eq. (3.25), the
form taken by this term involving the volume swept by the moving faces at different
times depends on the discretization scheme of time-derivatives. More details can
be found in [16].
In addition to changing the position of its control volumes, the mesh can also
change its topology if cells are removed or added. This is at the basis of auto-
matic mesh refinement (AMR), frequently used to locally (un)refine the mesh at
particular regions of interest (e.g. zones where the gradient of a given variable is
high). The introduction/removal of cells in the mesh requires defining the fields
and their fluxes in the newly generated cells/faces, which is based on a interpola-
tion procedure that uses the values in the neighboring cells.
kEk2 kEk2
f = fE = ∇· ε EE − I = ρE E − ∇ε (3.26)
2 2
where E is the electric field, ε = ε0 εR is the electric permittivity and ρE is the
charge density (per unit volume). In order to close the system of equations for
electrically-driven flows (EDFs), additional relations must be provided to compute
the terms in Eq. (3.26). Some options, the ones available in rheoTool , are presented
next. Note that when referring generically to EDFs, we do not exclude the possi-
bility of having any other external forcing (for example due to an imposed pressure
difference), in addition to the electric forcing. When only an electric forcing exists,
we call this flow as pure EDF.
The second term of Eq. (3.26) is only non-zero for a system of two fluids, each
having a different electric permittivity.
which closes the system of equations for an EDF. In Eq. (3.29), D is the diffu-
sion coefficient, e is the elementary charge, k is Boltzmann’s constant and T is
the absolute temperature. The last term of Eq. (3.29), representing the transport
of charged species due to an electric field, can be though as a standard convec-
tive term driven by an electromigration velocity (uM,i ). However, it may also be
considered as the Laplacian operator applied to field Ψ , with a space and time
varying diffusion coefficient, Di ezi
kT i
c (this last approach is used in rheoTool for
discretization purposes).
The so-called Poisson-Nernst-Planck model (henceforth PNP model) is con-
stituted by Eqs. (3.27)-(3.29) and, coupled with the continuity and momentum
equations, is applicable to a wide range of EDFs. However, the coexistence of dif-
ferent scales of time and length in EDFs may originate a stiff system of equations
when the PNP model is used. As such, several simplified models can be derived to
mitigate these numerical issues, as described next. Note that the PNP model does
not take into account molecular crowding effects (e.g., the number of ions near a
surface may grow unbounded), so care must be taken when using it to simulate
electrolytes of mild to high ionic strength.
In the PNP model, the electric-related unknowns are ci and Ψ . Due to the
convective term in Eq. (3.29), there is a two-way coupling between the PNP and
the momentum equations.
∇· (ε∇φExt ) = 0 (3.30a)
∇· (ε∇ψ) = −ρE (3.30b)
CHAPTER 3. Theoretical background 21
with bi = − ez
kT
i
and ai = zi ci,0 exp (bi ψ). All the terms of Eq. (3.32) with a star are
evaluated explicitly.
As for the PB model, the only electric-related unknowns are the two electric
potentials, ψ and φExt , computed from Eqs. (3.30a) and (3.33). Also, there is no
influence of flow variables in the DH model (one-way coupling).
uSch = µE (3.35)
where µ = − ηεζ0 is the electroosmotic mobility (ζ is usually the surface zeta-
potential). Thus, when Eq. (3.35) is used as a boundary condition for velocity
in the momentum equation, both the electroosmotic mobility and the electric field
at the surface must be known. The electroosmotic mobility is assumed to be known
a priori – it can be a fixed value over all the surface or have a known distribution.
On the other hand, the electric field on the surface must be computed, making use
of the initial assumption that no free charge exists in the bulk electrolyte, thus
Ψ = φExt + ψ = φExt , and
∇· (ε∇Ψ ) = 0 (3.36)
When the slip model is used, the electric body-force is not included in the
momentum equation – electric effects contribute uniquely via the slip boundary
condition on the wall.
Note that slip models do not resolve any phenomena occurring in the EDL.
Thus, this approach is highly inaccurate for some flows, even though the condition
λD
W
1 is satisfied. For example, this kind of model is unable to predict the high
values of shear-rate typically found in EDLs, which can trigger elastic instabilities
CHAPTER 3. Theoretical background 23
for complex fluid flows [17] – using a slip model would simply retrieve a smooth
flow in such cases.
F 2z2
σ= (D+ c+ + D− c− ) (3.37)
RT
ρE = F z(c+ − c− ) (3.38)
where R is the universal gas constant. Imposing the conservation of each variable
leads to (after the assumptions mentioned above; more details in Ref. [18])
∂σ
+ u· ∇σ = Deff ∇2 σ (3.39)
∂t
∇· (σ∇Ψ ) = 0 (3.40)
2D− D+
where the effective diffusivity is Deff = D− +D+
.The conductivity is transported
through Eq. (3.39), while Eq. (3.40), derived from the conservation of charge-
density (then simplified on the basis of electroneutrallity), is actually used to
compute the distribution of electric potential. The electric force entering the mo-
mentum equation assumes its standard form, taking into account that the charge
density can be expressed as ρE = −∇· (ε∇Ψ ) from Gauss’ law, then
Spring
Bead Molecule
Group of Molecules
Branch
Group of Molecules
Group of Molecules
RheoTool simulation
kT
I = DI, i=j
6πηa h
i
3D a 2 2 x x
Dij = 4 |xij |
1 + 3|x2aij |2 I + 1 − |x2aij |2 |xijij |ij2 , i =6 j ∧ |xij | ≥ 2a (3.44)
h i
D 1 − 9|xij | I + 3 xij xij ,
i 6= j ∧ |xij | < 2a
32a 32a |xij |
where |xij | = |rj − ri |, η is the fluid viscosity, a is the bead radius and I is the unit
tensor (3 × 3). Note that in rheoTool , a and D are defined independently by the
∂D
user and need not to be related. For the RPY tensor, ∂rjij = 0 in Eq. (3.43). The
decomposition of D, a symmetric tensor, to obtain σ is currently performed by a
Cholesky decomposition, whereby σ results in a lower triangular tensor (matrix).
The free-draining approach is sometimes assumed in bead-spring models, which
results from ignoring the beads disturbance in the continuum velocity field. In
such situations, hydrodynamic interactions are neglected and all the off-diagonal
tensor elements of tensor D become zero. The diffusion becomes isotropic and
can simply be defined by coefficient D. The summations√ in Eq. (3.43) reduce to a
single element contribution (j = i), and σii = DI.
The exclusion volume forces impose a repulsive potential between beads, which,
however, does not avoid any possible crossover between beads or springs (there is
also no collision between beads). The following exclusion volume force is used [21],
3 N
9 kT ν EV |xij |2 xij
3 X 9
FEV
i =− √ (2Nk,s ) 9/2
exp − Nk,s 2 (3.45)
2 l l3 4 π j=1
2 l l
and the models that do not impose any restriction on the spring length. Among
the mostly used models, the Marko-Siggia, Cohen Padé and FENE models fall into
the first category, while the Hookean model falls in the second one.
The Hookean model is arguably the simplest model representing a spring [20],
N
X
FSi = Hxij (3.46)
j∈gi
3kT N k,s
where H = l2
. As can be seen, the force is linearly proportional to the
spring extension and both are unlimited (here l is simply a parameter, and not
the effective limit of maximum spring extension). Although unphysical for high
deformations, the Hookean model is at the basis of the closed-form UCM and
Oldroyd-B constitutive equations used in continuum mechanics simulations, avail-
able in rheoTool . Some of the difficulties felt in the continuum simulations with
these two models arise precisely due to the unlimited stretch of Hookean springs,
which usually translate in unbounded stresses.
For the extension-limited models, we have [20]:
• Marko-Siggia model
N
" #
X 2 |xij | 1 1 xij
FSi = Hl − + 2 (3.47)
j∈gi
3 l 4 4 (1 − |xij |/l) |xij |
N
H 3 − (|xij |/l)2
X
FSi = 2
xij (3.48)
j∈g
3 1 − (|x ij |/l)
i
N
X 1
FSi = H xij (3.49)
j∈gi
1 − (|xij |/l)2
The relation between the spring force and the spring extension is non-linear in
these three models, and all are singular for |xij | = l, which represents an asymp-
tote for the springs extension. For low spring extension, the three models closely
approach the Hookean model. For high spring extension (close to the asymptote),
both FENE and Cohen Padé models present sharper gradients of force than the
Marko-Siggia model, which directly impacts the numerical stability of the time
integration algorithm.
computational cost. Most of the methods suggested in the literature are first-order
accurate in time, using Euler schemes to discretize the time derivative (higher-order
methods are not effective due to the random nature of the Brownian term [22]).
The explicit first-order Euler method evolves the beads positions from the
previous time-step (t) to the new time-step (t + ∆t) as,
" N i
#
Dij FSj
0.5 X
∗ t
X Dii FEV
i 6
ri = ri + ∆t uf + + + σij nj (3.51)
j=1,j6=i
kT kT ∆t j=1 t
Dii Fsi
rt+∆t
i = r∗i + ∆t (3.52)
kT t+∆t
In Eq. (3.51), the intermediate beads positions (r∗ ) are obtained explicitly from
the contribution of drag, Brownian and exclusion volume forces, and also from the
off-diagonal spring force terms. This results in a non-linear system of equations
that can be solved, for example, with the iterative Newton-Raphson method. As
discussed in [4], the Newton-Raphson method requires solving a linear system of
equations in each iteration (k ),
Overview of rheoTool
In the previous Chapter, the main theoretical points behind rheoTool were briefly
discussed. This Chapter focus on the numerical implementation of the governing
equations in the OpenFOAM
R
environment, providing an overview of the func-
tionalities available in rheoTool .
30
Table 4.1: Available constitutive models in the constitutiveEquations library.
GNF models
1
Model TypeName ηs (γ̇)
Newtonian Newtonian η
2 (Bounded) max ηmin , min ηmax , k γ̇ n−1
Power-Law PowerLaw
n−1
Carreau-Yasuda CarreauYasuda η∞ + (η0 − η∞ )[1 + (k γ̇)a ] a
Bounded: min η0 , τ0 γ̇ −1
+ k γ̇ n−1
2 Herschel-Bulkley HerschelBulkley 3 Papanastasiou reg.: min η , τ γ̇ −1 [1 − exp(−mγ̇)] + k γ̇ n−1
0 0
√ q 2
Bounded: max ηmin , min ηmax , η∞ + τγ̇0
2 Casson Casson n√ q √ o2
Papanastasiou reg.: η∞ + τγ̇0 1 − exp − mγ̇
1
Corresponds to the name entry identifying the model in the source code.
2
Special care is taken in these models to avoid division by zero when γ̇ is zero or very small and n − 1 < 0. For γ̇ < VSMALL, the
value γ̇ =VSMALL is used in the computation of the shear viscosity (VSMALL = 10−300 for versions using double precision).
3
The original Papanastasiou regularization does not include the artificial upper-bounding by η0 . However, this bounding is needed
to avoid an infinite viscosity for γ̇ → 0 (e.g. startup of flow) and n < 1. The original Papanastasiou regularization is recovered
for η0 → ∞. In practice, η0 should be low enough to avoid an infinite viscosity in quiescent conditions and high enough to allow
Papanastasiou regularization to take control in the remaining situations.
Notes:
q
• γ̇ = γ̇:2γ̇ , with γ̇ = ∇u + ∇uT .
D ∂φ
• I is the identity tensor and Dt (φ) = ∂t + u· ∇φ represents the material derivative of the generic variable φ.
∇ ∂τ
• τ= ∂t + u· ∇τ − τ· ∇u − ∇uT · τ is the upper-convected derivative of τ.
∇
• τ = τ +ζ (τ· D + D· τ) is the Gordon-Schowalter derivative of τ, with D = 12 (∇u + ∇uT ).
31
Continuation of Table 4.1
Viscoelastic models solved in the standard extra-stress or conformation tensor variables
λ ∇
h i
D 1
1 + λ Dt f τ+ f τ= ηp (∇u + ∇uT )
FENE-CR FENE-CR ηs ηp λ L2 + ηλ tr(τ)
p
where f = L2 −3
λ ∇
aηp T D 1
τ+ f τ= f (∇u + ∇u ) − Dt f [λτ + aηp I]
FENE-P FENE-P ηs ηp λ L2 + aηλ tr(τ)
p L2
where f = L2 −3
and a = L2 −3
∇ δ
−(A − I) − 2k λλDR 1 − 3/tr(A) A + β tr(A)
p
λD A= 3 (A − I)
3Rolie-Poly
Rolie-Poly ηs ηp λD 3− 2χ
2
1− 21 q
and χ = tr(A)
χmax χmax
where k =
3
χ2 1
1− 2 3− 2
χmax χmax
∇ η
f τ + λB τ +α ληBp (τ· τ) + λBp (f − 1) I = ηp (∇u + ∇uT )
eXtended Pom-Pom XPomPom ηs ηp λB 2
h i q
(Λ−1)
where f = 2 λλBS e q 1 − Λ1 + Λ12 1 − α3 (ηtr(τ·τ) tr(τ)
P /λB )
2 and Λ = 1+ 3ηP /λB
3 ηp
See Ref. [23]. This model is exclusively solved in the conformation tensor variable, which is then converted to τ using, τ = λD k(A − I).
32
Continuation of Table 4.1
‡Viscoelastic models solved with the log-conformation approach
4,5
Model TypeName Θ¸τ Constitutive Equation
ηp Θ
6Oldroyd-B 1
e−Θ − I
Oldroyd-BLog τ= λ (e − I) Υ= λ
7 WhiteMetzner
ηp Θ 1
e−Θ − I
WhiteMetznerCYLog τ= λ (e − I) Υ= λ(γ̇)
(Carreau-Yasuda)
h 2 i
ηp Θ 1
e−Θ − I − αeΘ e−Θ − I
Giesekus GiesekusLog τ= λ (e − I) Υ= λ
n o
ηp Θ 1 ε
tr(eΘ ) − 3 (e−Θ − I)
PTT linear PTTlinearLog τ= λ(1−ζ) (e − I) Υ= λ 1+ 1−ζ
ε
ηp (tr(eΘ )−3)
PTT exponential PTTexpLog τ= λ(1−ζ) (e
Θ − I) Υ = λ1 e 1−ζ (e−Θ − I)
ηp f Θ f L2
e−Θ − I , where f =
FENE-CR FENE-CRLog τ= λ (e − I) Υ= λ L2 −tr(eΘ )
ηp L2 L2
Θ 1
ae−Θ − f I , where a =
FENE-P FENE-PLog τ= λ (f e − aI) Υ= λ L2 −3
and f = L2 −tr(eΘ )
Θ δ
ηp tr(e )
2k λλDR
p
8Rolie-Poly Rolie-PolyLog τ= λD k(e
Θ − I) Υ= − λ1D e−Θ (eΘ − I) + Θ Θ
1 − 3/tr(e ) e + β 3
Θ
(e − I)
‡
The solvent viscosity, the polymeric viscosity coefficient and the relaxation time for the models solved in variable Θ are the same as those for the models solved in
variable τ or A, in the previous page.
4
For the shortness of notation, we have introduced the operator: Υ = ∂Θ ∂t + u· ∇Θ − (ΩΘ − ΘΩ) − 2B.
5
The following equivalences hold true: eΘ = A = RΛRT and e−Θ = A−1 = RΛ−1 RT .
6
For this model, we also included the square-root conformation approach [11] (TypeName: Oldroyd-BSqrt) and the rootk kernel approach [10] (TypeName: Oldroyd-
BRootk ), for demonstration purposes.
7 ηp (γ̇) ηp
This log-conformation tensor approach of the White-Metzner model is only applicable when λ( γ̇) = λ is constant, i.e., for K = L, a = b and n = m. The
version based on the extra-stress tensor variable is more general and does not have this restriction.
8
The expression for k is the same as for the model solved in variable A, in the previous page, considering that A = eΘ .
33
CHAPTER 4. Overview of rheoTool 34
γ̇ : γ̇ √
r
1
γ̇ = = 2D : D, with γ̇ = ∇u + ∇uT and D = γ̇ (4.1)
2 2
In the code, the shear-rate is returned by function strainRate() as
strainRate() = sqrt(2.0)*mag(symm(fvc::grad(U())))
and it is equivalent to Eq. (4.1). Indeed,
1
∇u + ∇uT = 21 γ̇ = D
symm(fvc::grad(U())) = 2
thus,
√ q1 q √
sqrt(2.0)*mag(symm(fvc::grad(U()))) = 2 2 γ̇ : 12 γ̇ = γ̇:2γ̇ = 2D : D
which is equal to Eq. (4.1) – the definitions of operators symm(), mag() and :
(double contraction) can be found in the OpenFOAM
R
programmers’ guide. Note
that the invariant computed in Eq. (4.1) is actually the magnitude of the rate-of-
strain tensor, which is usually called shear rate or strain rate for shear-dominated
or extensional-dominated flows, respectively.
All the viscoelastic models presented in Table 4.1 can be solved in the standard
extra-stress tensor τ (Eq. 3.4) or using the log-conformation approach (Eq. 3.8).
The selection is made in dictionary constitutiveProperties, which should be
located inside the folder constant/ of the case (see more details in section 5.1.1).
For the Oldroyd-B model, we provide two additional methods for demonstration
purposes. One of them (TypeName: Oldroyd-BSqrt) consists in solving the con-
stitutive equation using the square-root of the conformation tensor, according to
Ref. [11]. The second approach (TypeName: Oldroyd-BRootk ) allows to apply a
general rootk kernel, as described in Ref. [10]. Both can be used in 2D or 3D
simulations, as any other model in the library. Since both models are only illus-
trative, their implementation and theory are not described in this guide, although
both can be easily understood after a close inspection of the source code and tak-
ing as reference the literature cited for each one. Furthermore, tutorials for both
methodologies are included in rheoTool (see the tutorial of Section 5.1.7).
! Other models:
+ VCM model (TypeName: VCM )
The Vasquez-Cook-McKinley (VCM) model [24] can be used to simulate worm-
like micellar solutions, being able to predict the shear-banding behavior typically
observed in these fluids. The model represents such fluids as a combination of
large (subscript A) and small chain (subscript B) species that can convert into
each other. A transport equation is solved for each species [24],
∂nA 1 cA n A
+ u · ∇nA = 2DA ∇2 nA + cB n2B − (4.2)
∂t 2λA λA
∂nB cB n2B cA n A
+ u · ∇nB = 2DB ∇2 nB − +2 (4.3)
∂t λA λA
CHAPTER 4. Overview of rheoTool 35
where n is the dimensionless number density of the specie, λ is the relaxation time,
D is the diffusivity coefficient and cA and cB are, respectively, the dimensionless
breakage and reformation rates, expressed as
χ A
cA = cAEq + γ̇ : (4.4)
3 nA
cB = cBEq (4.5)
In Eqs. (4.2) and (4.3), the double contraction term originally presented in [24]
is not included, in order to simplify the definition of no-flux boundary conditions
for nA and nB at impermeable walls (which reduce to a zero-gradient condition).
The contribution of these omitted terms is typically negligible.
A constitutive equation is also solved for each species [24],
∇
λA A +A − nA I − λA DA ∇2 A = cB nB B − cA A (4.6)
∇ nB I
λA B +B − − λA DB ∇2 B = −2cB nB B + 2cA A (4.7)
2
where A and B represent the conformation tensor of each species, and = λλAB . The
contribution of each species to the polymeric extra-stress tensor is given by [24],
ln L∗ + m
Dr,0
Dr = ∗3 (4.10)
L m
is a time-varying diffusion coefficient (units are s−1 ), L∗ = L/L0 is the time-varying
rod length normalized by its initial value and m is the initial aspect ratio of the
rods (see [25] for more details). In Eq. (4.9), the last term is approximated by [25]
1
∇uT :< uuuu >≈ [S · D + D · S − S · S · D − D · S · S + 2S · D · S + 3(S : D)S]
5
(4.11)
CHAPTER 4. Overview of rheoTool 36
∂Θ g(τ) −Θ
+ u· ∇Θ − (ΩΘ − ΘΩ) − 2B = e −I (4.18)
∂t λ
with
1
σ−τ0 n
ηp max 0, kσ n , Oldroyd-B–Herschel-Bulkley
max 0, σ−τ0
, Oldroyd-B–Bingham
σ
g(τ) = σ−τ0
n ε
Θ
o
max 0, σ 1 + 1−ζ tr(e ) − 3 , lin. PTT–Bingham
σ−τ0
ε (tr(eΘ )−3)
max 0, σ e 1−ζ , exp. PTT–Bingham
(4.19)
ηp
and the polymeric extra-stress tensor is recovered using τ = λ(1−ζ)
(eΘ − I).
+ ML-IKH model (TypeName: ML-IKH )
Elastoviscoplastic models can be rendered more complex (and realistic) once
thixotropy is added to them. This is the case of the Multi-Lambda Isotropic Kine-
matic Hardening (MLK-IKH) model [29]. According to this model, the equation
governing the polymeric extra-stress tensor is given by [29]
λky ∇
max 0, 1 − τeff + λλE τ= ληp (∇u + ∇uT ) (4.20)
σ
CHAPTER 4. Overview of rheoTool 38
∂A
+ u · ∇A + Ω · A − A · Ω = Dp · A + A · Dp + Dp − qdp A (4.21)
∂t
q
where Ω = ∇u − ∇uT /2 is the vorticity tensor, dp = Dp2:Dp and
(
0 , if σ < λky
Dp = (σ−λky ) τeff (4.22)
2ληp
· σ , if σ ≥ λky
is the plastic component of the rate of deformation tensor. In the previous equa-
tions, λ is a structure parameter regulating the thixotropic behavior and is ob-
tained from the multi-lambda model [29],
M
X
λ= C i λi (4.23)
i=1
transforming the conformation tensor (A) in the extra-stress tensor (τ), using the
relations in Table 4.1 (for the models expressed in the log-conformation approach,
considering that eΘ = A). The same applies for the Roly-Polie model.
In order to write the constitutive equation for FENE-type models as a function
of τ, some terms arise, which may compromise the numerical stability. Further-
more, the computational cost to evaluate the resulting expression is higher than
for the original model. As such, some authors simplify the constitutive equation
by neglecting certain terms [30]. For the FENE-CR and FENE-P models, the
complete constitutive equation written as a function of A and τ and the modified
formulation in τ are:
• FENE-CR
– Complete in A:
∇
L2
λ A= −f (tr(A))(A − I), where f (tr(A)) = L2 −tr(A)
• FENE-P
– Complete in A:
∇
L2 L2
λ A= − [f (tr(A))A − aI], where f (tr(A)) = L2 −tr(A)
and a = L2 −3
– Modified in τ:
L2 + aηλ tr(τ)
λ ∇ L2
τ+ f
τ= aηf p (∇u + ∇uT ), where f = L2 −3
p
and a = L2 −3
In rheoTool , all the formulations are available and can be used (see Section
5.1.1 to know how to select each one). The steady material functions evaluated
for canonical flows are the same for all the formulations. However, this is not true
when evaluating the transient material functions: the modified formulations have
a different behavior comparing with the complete ones, which are themselves simi-
lar. For a generic flow, the complete formulations, either in A or τ, should provide
similar results, since they are mathematically equivalent. Due to discretization
CHAPTER 4. Overview of rheoTool 40
errors and stability issues, this may not be true. Regarding the modified formula-
tions, they are not expected to behave exactly as the complete ones, even in the
limit of highly refined grids.
From our experience, we strongly recommend using the formulations written
and solved as a function of A for FENE-type models. Those are the most sta-
ble, the most accurate regarding the original theory presented in [6] and the ones
for which there is direct correspondence with the models solved with the log-
conformation approach, since those were derived from the constitutive equations
written as a function of the conformation tensor. Note that the FENE-CR and
FENE-P models available in the viscoelasticTransportModels library of viscoelas-
ticFluidFoam [1] are expressed in the modified form presented above.
• lines 1-91: this section initializes the variables used in the constitutive model.
In terms of field variables, we have (lines 23-84): tau (τ), theta (Θ), eigVals
(Λ) and eigVecs (R). All those fields must be defined by the user when
CHAPTER 4. Overview of rheoTool 41
• lines 94-151: this section implements the member function correct(), whose
purpose is to update the polymeric extra-stress field, by evolving Θ accord-
ing to the constitutive equation. From line 96 to 105, variables M, Ω and
B, defined in Eqs. (3.9)–(3.11), are computed. The function decomposeG-
radU() is a member function of the base class constitutiveEq (find it in the
file constitutiveEq.C), since it is used by all the models based on the
log-conformation tensor approach. Then, in lines 107-140, the constitutive
equation (Eq. 3.8) is built and solved, after which Θ is diagonalized to com-
pute its eigenvectors/eigenvalues (line 144). The function doing this task
(calcEig()) is also a member function of the class constitutiveEq and the al-
gorithm being used by default for that purpose is the QR method provided
by the Eigen library [7]. Another method is also available, as discussed in
Section 4.1.5. Note that the eigenvalues retrieved by function (calcEig())
are already exponentiated, so that they correspond to Λ = exp(ΛΘ ). Fi-
nally, with the currently computed eigenvectors/eigenvalues, the polymeric
extra-stress tensor (τ) is recovered from the conformation tensor (line 148),
according to the relation established in Eqs. (3.6) and (3.14) (check Table
4.1 for other models), and will be used in the divTau() function described
below.
1 #include "Oldroyd_BLog.H"
#include "addToRunTimeSelectionTable.H"
3
// * * * * * * * * * * * * * * Static Data Members * * * * * * * *
* * * * * //
5
namespace Foam{
7 namespace constitutiveEqs{
defineTypeNameAndDebug(Oldroyd_BLog, 0);
9 addToRunTimeSelectionTable(constitutiveEq, Oldroyd_BLog,
dictionary);
}
11 }
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * *
* * * * * //
13
Foam::constitutiveEqs::Oldroyd_BLog::Oldroyd_BLog
15 (
const word& name,
17 const volVectorField& U,
const surfaceScalarField& phi,
19 const dictionary& dict
)
21 :
constitutiveEq(name, U, phi),
23 tau_
(
25 IOobject
(
27 "tau" + name,
U.time().timeName(),
29 U.mesh(),
IOobject::MUST_READ,
31 IOobject::AUTO_WRITE
),
33 U.mesh()
),
35 theta_
(
37 IOobject
(
39 "theta" + name,
U.time().timeName(),
CHAPTER 4. Overview of rheoTool 43
41 U.mesh(),
IOobject::MUST_READ,
43 IOobject::AUTO_WRITE
),
45 U.mesh()
),
47 eigVals_
(
49 IOobject
(
51 "eigVals" + name,
U.time().timeName(),
53 U.mesh(),
IOobject::READ_IF_PRESENT,
55 IOobject::AUTO_WRITE
),
57 U.mesh(),
dimensionedTensor
59 (
"I",
61 dimless,
pTraits<tensor>::I
63 ),
zeroGradientFvPatchField<tensor>::typeName
65 ),
eigVecs_
67 (
IOobject
69 (
"eigVecs" + name,
71 U.time().timeName(),
U.mesh(),
73 IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
75 ),
U.mesh(),
77 dimensionedTensor
(
79 "I",
dimless,
81 pTraits<tensor>::I
),
83 zeroGradientFvPatchField<tensor>::typeName
),
85 rho_(dict.lookup("rho")),
etaS_(dict.lookup("etaS")),
87 etaP_(dict.lookup("etaP")),
lambda_(dict.lookup("lambda")),
89 {
checkForStab(dict);
91 }
93 // * * * * * * * * * * * * * * * Member Functions * * * * * * * *
* * * * * //
void Foam::constitutiveEqs::Oldroyd_BLog::correct()
95 {
CHAPTER 4. Overview of rheoTool 44
// Decompose grad(U).T()
97
volTensorField L = fvc::grad(U());
99
dimensionedScalar c1( "zero", dimensionSet(0, 0, -1, 0, 0, 0,
0), 0.);
101 volTensorField B = c1 * eigVecs_;
volTensorField omega = B;
103 volTensorField M = (eigVecs_.T() & L.T() & eigVecs_);
137 );
139 thetaEqn.relax();
thetaEqn.solve();
141
// Diagonalization of theta
143
calcEig(theta_, eigVals_, eigVecs_);
145
// Convert from theta to tau
147
tau_ = (etaP_/lambda_) * symm( (eigVecs_ & eigVals_ & eigVecs_
.T()) - Itensor);
149
CHAPTER 4. Overview of rheoTool 45
tau_.correctBoundaryConditions();
151 }
Listing 4.1: Source code for the Oldroyd-BLog constitutive model
(Oldroyd BLog.C)
1 tmp<fvVectorMatrix> constitutiveEq::divTau
(
3 const volVectorField& U
) const
5 {
7 if (isGNF())
{
9 return
(
11 fvm::laplacian( eta()/rho(), U, "laplacian(eta,U)")
+ (fvc::grad(U) & fvc::grad(eta()/rho()))
13 );
}
15 else
{
17 if (stabMeth_ == 0) // none
{
19
return
21 (
fvc::div(tau()/rho(), "div(tau)")
23 + fvm::laplacian(etaS()/rho(), U, "laplacian(eta,U
)")
);
25
}
27 else if (stabMeth_ == 1) // BSD
{
29
return
31 (
fvc::div(tau()/rho(), "div(tau)")
33 - fvc::laplacian(etaP()/rho(), U, "laplacian(eta,U
)")
+ fvm::laplacian( (etaP()+ etaS())/rho(), U, "
laplacian(eta,U)")
35 );
37 }
else // coupling
39 {
41 return
(
43 fvc::div(tau()/rho(), "div(tau)")
- (etaP()/rho()) * fvc::div(fvc::grad(U))
45 + fvm::laplacian( (etaP() + etaS())/rho(), U, "
laplacian(eta,U)")
);
CHAPTER 4. Overview of rheoTool 46
47
}
49 }
}
Listing 4.2: Source code of the virtual function divTau() defined in file
constitutiveEq.C
• inside the source .C and header .H files, find-replace the old model’s name
by the new one (e.g. Ctrl + H in gedit). This is to change the name of the
class, which is usually equal to the name given to the respective source file.
However, it is a good idea to always check first the name given to the class
of the template model.
CHAPTER 4. Overview of rheoTool 47
• add the source code of the new model to the compilation list of the library.
For that, edit file src/libs/constitutiveEquations/Make/files
by adding the path for the source code (see the entries for the other models
already there).
• make a first compilation of the new model, by running the script Allwmake
in directory src/. Note that, until this point, the source code of the new
model is the one from the original template model, where only the name of
the class has been changed. Thus, the model should compile without errors.
If not, something wrong occurred in the previous steps.
• the last step is to change the source code of the model in order to implement
the desired constitutive equation. Typically, the changes will be in three
main places: (i) the header file, where the new variables and parameters
of the model have to be declared (delete the ones from the template model
that are not needed); (ii) the constructor in the source file, where those new
entries should be added and initialized (delete the ones not needed); (iii)
function correct(), which is aimed to either update the viscosity (GNF) or
the polymeric extra-stresses (viscoelastic model). Note that you may also
need to define functions divTau() and divTauS() for your new model, if the
ones defined and used by default in the base class (see constitutiveEq.C)
are not adequate. After all the changes on the code had been completed,
compile again by running script Allwmake.
If all the steps listed above were successfully executed, then the new model
is now available to all the solvers of rheoTool . In order to use library constitu-
tiveEquations in any solver other than the ones provided with rheoTool , the user
should:
1 2,3,4 5
Model TypeName fE Governing Equations Section
∇· (ε∇φExt ) = 0
N
Poisson-Nernst-Planck P PN
NernstPlanck − F zi ci (∇Ψ − Ea ) ∇· (ε∇ψ) = −F zi ci Section 3.6.1
(PNP) i=1
i=1
∂ci + u· ∇c = ∇· (D ∇c ) + ∇· D ezi ∇Ψ c
∂t i i i i kT i
∇· (ε∇φExt ) = 0
N
Poisson-Boltzmann P ezi
N N N
PoissonBoltzmann −F zi ci,0 exp − kT ψ (∇Ψ − Ea ) ∇· (ε∇ψ) + ψF
P
(ai bi )∗ = −F
P
(ai )∗ + ψ ∗ F
P
(ai bi )∗ Section 3.6.3
(PB)
i=1 i=1 i=1 i=1
with bi = − ez
kT and
i
ai = zi ci,0 exp (bi ψ)
N
∇· (ε∇φExt ) = 0
Debye-Hückel P ezi
DebyeHuckel −F zi ci,0 1 − kT ψ (∇Ψ − Ea ) N Section 3.6.4
(DH) ∇· (ε∇ψ) = −F
P
zi ci,0 1 − ezi
kT ψ
i=1
i=1
48
CHAPTER 4. Overview of rheoTool 49
• lines 35-107: this is the constructor of the main class, where the several
fields and variables are initialized. The call to function checkForPhiE() in
line 44, which is implemented in the base class (see EDFEquation.C), is
checking for the existence of field phiE in the folder corresponding to the
starting time. If it is found, the code will interpret that the electric potential
should be split into 2 variables (phiE ⇔ φExt and psi ⇔ ψ), otherwise, the
code will consider that only a single potential should be used (psi ⇔ Ψ ).
This is how we identify if the splitting of potentials approach should be used.
We also highlight lines 89-107, where each specie of the PNP model is being
constructed and saved in a P trList <>, named species .
• lines 110-156: as suggested by the function’s name (Fe), these lines im-
plement the function returning the electric body-force for the momentum
equation. The charge density (rhoE) is computed in lines 114-120 (compare
with Eq. 3.28), and multiplied by the electric field in lines 122-155 (compare
with Eq. 3.26 and Table 4.2). The computation of the electric field may
include, or not, the contribution from the intrinsic potential, as discussed in
Section 3.6.2 – this is a user selection. Furthermore, the vector extraE is
an extra, uniform electric field, which can be optionally defined by the user
(see Table 4.2 and its footnotes).
• lines 158-252: it is inside this function, named correct(), that the electric-
related equations are solved for. The function is generally prepared to use
the splitting approach, in which case three equations are solved: the Laplace
equation for the external potential (lines 172-186), the Poisson equation for
the intrinsic (or full, unique) potential (lines 188-218) and the Nernst-Planck
transport equation for each ionic specie (220-250) – all these equations can
be seen in Table 4.2. Each equation is inserted in a while loop controlled by
the number of cycles and by the initial residual of the equation solved for.
This is to optionally converge the explicit terms inside each equation, for
each inner-iteration. In addition, and as discussed in Section 4.2.3, all the
equations are solved inside an electrokinetic coupling loop (lines 161-251),
whose number of iterations is controlled by variable nIterPNP , that is read
from dictionary fvSolution (line 86). If the variable is not specified by
the user, 2 iterations are carried out by default.
All the other EDF models also have the functions Fe() and correct() in their
structure, which are defined according to the given model. We believe that the
readers/users will easily understand those functions by reading the comments in-
cluded in the code, and by comparing the code with Table 4.2.
CHAPTER 4. Overview of rheoTool 51
#include "NernstPlanck.H"
2 #include "addToRunTimeSelectionTable.H"
58 (
IOobject
60 (
"phiE" + name,
62 phi.time().timeName(),
phi.mesh(),
64 IOobject::READ_IF_PRESENT,
solvePhiE_ == false ? (IOobject::NO_WRITE) : (IOobject::
AUTO_WRITE)
66 ),
phi.mesh(),
68 dimensionedScalar
(
70 "zero",
psi_.dimensions(),
72 pTraits<scalar>::zero
),
74 zeroGradientFvPatchField<scalar>::typeName
),
76 relPerm_(dict.lookup("relPerm")),
T_(dict.lookup("T")),
78 extraE_(dict.lookupOrDefault<dimensionedVector>("extraEField",
dimensionedVector("0", dimensionSet(1, 1, -3, 0, 0, -1, 0),
vector::zero))),
psiContrib_(dict.lookupOrDefault<bool>("psiContrib", true)),
80 phiEEqRes_(phi.mesh().solutionDict().subDict("electricControls").
subDict("phiEEqn").lookupOrDefault<scalar>("residuals", 1e-7)),
psiEqRes_(phi.mesh().solutionDict().subDict("electricControls").
subDict("psiEqn").lookupOrDefault<scalar>("residuals", 1e-7)),
82 ciEqRes_(phi.mesh().solutionDict().subDict("electricControls").
subDict("ciEqn").lookupOrDefault<scalar>("residuals", 1e-7)),
maxIterPhiE_(phi.mesh().solutionDict().subDict("electricControls").
subDict("phiEEqn").lookupOrDefault<scalar>("maxIter", 50)),
84 maxIterPsi_(phi.mesh().solutionDict().subDict("electricControls").
subDict("psiEqn").lookupOrDefault<scalar>("maxIter", 50)),
maxIterCi_(phi.mesh().solutionDict().subDict("electricControls").
subDict("ciEqn").lookupOrDefault<scalar>("maxIter", 50)),
86 nIterPNP_(phi.mesh().solutionDict().subDict("electricControls").
lookupOrDefault<int>("nIterPNP", 2)),
species_(),
88 nSpecies_(0)
{
90 PtrList<entry> specEntries(dict.lookup("species"));
nSpecies_ = specEntries.size();
92 species_.setSize(nSpecies_);
);
106 }
}
108
// * * * * * * * * Member Functions * * * * * * * * * * //
110 Foam::tmp<Foam::volVectorField> Foam::EDFEquations::NernstPlanck::Fe()
const
{
112 volScalarField rhoE( psi_ * dimensionedScalar("norm", epsilonK_.
dimensions()/dimArea, 0.) );
122 if (solvePhiE_)
{
124 if (psiContrib_)
{
126 return
(
128 -rhoE * ( fvc::grad(phiE_+psi_) - extraE_)
);
130 }
else
132 {
return
134 (
-rhoE * ( fvc::grad(phiE_) - extraE_)
136 );
}
138 }
else
140 {
if (psiContrib_)
142 {
return
144 (
-rhoE * ( fvc::grad(psi_) - extraE_)
146 );
}
148 else
{
150 return
(
152 -rhoE * (-extraE_)
);
154 }
}
156 }
160
// Electrokinetic coupling loop
162 for (int j=0; j<nIterPNP_; j++)
{
164
Info << "PNP Coupling iteration: " << j << endl;
166
scalar res=GREAT;
168 scalar iter=0;
170 //- Equation for the external potential (loop for the case
// of non-orthogonal grids)
172 if (solvePhiE_)
{
174 while (res > phiEEqRes_ && iter < maxIterPhiE_)
{
176 fvScalarMatrix phiEEqn
(
178 fvm::laplacian(phiE_)
);
180
phiEEqn.relax();
182 res=phiEEqn.solve().initialResidual();
184 iter++;
}
186 }
190 res=GREAT;
iter=0;
192
volScalarField souE(psi_ * dimensionedScalar("norm1",dimless/dimArea
,0.));
194
forAll (species_, i)
196 {
souE +=
198 (
-species_[i].zi()*species_[i].ci()*FK_
200 /(relPerm_*epsilonK_)
);
202 }
214 psiEqn.relax();
res=psiEqn.solve().initialResidual();
CHAPTER 4. Overview of rheoTool 55
216
iter++;
218 }
• copy-paste the folder of an existing model (that we will call template model)
in directory src/libs/EDFModels/models/. Rename the folder and
the files inside it (.C and .H files) with the new model’s name. Remove the
.dep file inside the folder. We recommend to use model slipSmoluchowski as
template, since it is the simplest one and does not contain other sub-classes,
which would also need to be renamed in the next step.
• inside the source .C and header .H files, find-replace the old model’s name
by the new one (e.g. Ctrl + H in gedit). This is to change the name of the
class, which is usually equal to the name given to the respective source file.
However, it is a good idea to always check first the name given to the class
of the template model.
• add the source code of the new model to the compilation list of the library.
For that, edit file src/libs/EDFModels/Make/files by adding the
path for the source code (see the entries for the other models already there).
• make a first compilation of the new model, by running the script Allwmake
in directory src/. Note that, until this point, the source code of the new
model is the one from the original template model, where only the name of
the class has been changed. Thus, the model should compile without errors.
If not, something wrong occurred in the previous steps.
• the last step is to change the source code of the model in order to implement
the desired EDF model. Depending on the characteristics of the new model,
several parts of the source and header files might need to be changed. Our
recommendation is to look to the closest model among the ones available.
After all the changes on the code had been completed, compile again by
running script Allwmake.
If all the steps listed above were successfully executed, then the new model is
now available to solver rheoEFoam (only). In order to use library EDFModels in
any solver other than rheoEFoam, the user should:
• create an EDFModel object by calling the constructor with the correct ar-
guments, for example: EDFModel elecM(phi).
• add the library EDFModels to the Make/options, and specify its path
(check the Make/options of rheoEFoam for an example).
CHAPTER 4. Overview of rheoTool 57
local fields of sPCloudInterface and the particles fields inside spc , from solidPar-
ticle class. We believe that the comments left in the code will help in this task. It
would be unfeasible (and probably useless) to explain all the code details in this
guide, thus we decided to only explore the function controlling the main Brownian
dynamics loop. This function, named update(), a public member of sPCloudInter-
face, is charged of evolving the position and configuration of the molecules each
time-step. In what follows, we briefly discuss the structure of this function (Listing
4.4):
• line 3: the beads positions at the beginning of a time-step (mx ) are copied
to mx0 . The algorithm structure ensures that the beads positions in mx0
are such that all the springs of the active molecules are shorter than l, for
bounded models (Section 3.7.3).
• lines 5-6: if the user selects to account for hydrodynamic interactions, tensor
D is computed from the RPY model (Eq. 3.44), and its Cholesky decompo-
sition results in tensor σ (Section 3.7.3). If hydrodynamic interactions are
suppressed, the diffusion is isotropic and can be represented by a scalar (D;
Eq. 3.44), which is used to compute each force acting on the beads (tensors
D and σ are not even computed in this case).
• line 8: this function computes the Brownian force contribution to the beads
velocity (last term of Eq. 3.43). The function is also defined inside sPClou
dInterface.C.
• lines 10-11: if exclusion volume forces are activated by the user, then their
contribution to the beads velocity is added through a call to function fEV(),
defined inside sPCloudInterface.C.
• line 13: function sendU() copies the beads velocity from the local mU
P trList <> to the particles field U. This function is defined inside sPClou
dInterface.C.
• line 16: function moveAndReceive() comprises two main steps. In the first
step, there is a call to the move() function of the solidParticleCloud object
spc . This executes the movement and tracking of all the beads, after adding
the drag force contribution to the particles velocity, which is done outside
class sPCloudInterface (see file solidParticle.C). In the second step,
the function updates mx with the final positions resulting from the particle
tracking. If for some reason a given particle (bead) is lost during the tracking
(e.g. if it exits the mesh through a patch), then the corresponding molecule
is labeled as non-active, and it is no more tracked. Up to this point, the
beads experienced all the forces, except the spring force.
• line 19: the current beads positions are saved in mxStar . Since the spring
force still did not contributed to these positions, the computation of the
spring vectors from mxStar would result eventually in overstretched springs
(Ri > l).
CHAPTER 4. Overview of rheoTool 59
• line 20: this call to a non-member function computes explicitly (Euler ex-
plicit) the spring force contribution to the beads velocity (mU ), using the
current beads positions (mxStar ). The local beads positions (mx ) are also
updated.
• line 23: based on the current beads positions (mx ), the spring vectors are
computed and a check is carried out for Ri < αl (see Section 3.7.4). For
each molecule, if any spring violates this condition mxStar is taken again
and the spring force contribution is now added implicitly, using the Newton-
Raphson method (Section 3.7.4). Of course, this is only done if the semi-
implicit method is selected (Euler-explicit is also available), and for any of
the bounded spring models. After the implicit call, the function checks if any
spring is overstretched. This can happen if the time-step is too large. Any
molecule having at least one spring overstretched is automatically deleted by
the algorithm and a warning message is printed to the terminal.
• lines 26-27: if the molecules are not tethered and if any of the components
of the push-back vector defined by the user is non-zero, the molecules center
of mass is pushed to its original position. This is equivalent to the transla-
tion of all the molecule’s beads by a fixed vector. In this case, we add the
corresponding velocity to mU , such that the translation can be effective in
the next particle tracking stage.
• line 29: see comment above for line 13.
• line 32: this is the second call to function moveAndReceive(), thus a second
call to the particle tracking engine. The operations carried out are the same
as described above for line 16, with the exception that, in this call, the beads
velocity does not receive the drag force contribution, since this was already
done in the first call (the distinction between calls is in the boolean argu-
ment passed to the function). At this point of the algorithm, the particles
positions and mx are synchronized and it can be ensured that none of the
active molecules has an overstretched spring. We do not care about the syn-
chronization of the beads’ velocity, because this field is not used in the next
time-step, contrarily to the beads’ positions.
• line 35: function writeM() ensures that the molecules’ data is written in the
case directory at each outputTime. It will create directories outputTime
/lagrangian/molecules/ and constant/runTimeInfo/outpuTi
me/, and write the molecules’ data therein. Both directories are needed on
restart of rheoBDFoam.
• lines 38-39: function writeStatistics() can be optionally activated by the
user, and will retrieve the molecules index/position/stretch over time.
• lines 41-48: these lines enclose the conditional return value of function up-
date(). A value of true is returned as long as at least one molecule is still
active and under tracking. Otherwise, the function returns false and, as ex-
plained in Section 4.4.5, this will force rheoBDFoam to abort the run, even if
CHAPTER 4. Overview of rheoTool 60
the endTime was still not reached (there is no interest in keeping the solver
running without any valid molecule).
While inspecting the source code, note that several intermediate operations
are carried out in a dimensionless form to reduce round-off errors, but final results
are always dimensional. Round-off errors can be an important source of numer-
ical error in Brownian dynamics simulations (the situation is worse in atomistic
simulations) if care is not taken with the different scales involved.
bool sPCloudInterface::update()
2 {
mx0_ = mx_;
4
if (isHI_)
6 computeDSigma();
8 fBrownian();
10 if (isExclusionVolumeF_)
fEV();
12
sendU();
14
// Move particles: Drag + Brownian + Exclusion Volume
16 moveAndReceive(true);
44 }
else
46 {
return false;
48 }
}
Listing 4.4: Member function update() of class sPCloudInterface. The source
code can be found in file sPCloudInterface.C.
• Analytical : this scheme is the only available for analytical forcing and cannot
be used with a numerical forcing. The velocity is interpolated using uf =
∇uT · ri , where ∇u is user-defined.
CHAPTER 4. Overview of rheoTool 62
• Gradient: in this approach, both the cell-centered field and its gradient are
used for the interpolation. Assuming that uC is the forcing known (com-
puted) at the cell center (located at xC ) and that ∇uC is its gradient (com-
puted numerically), then the velocity at any point xP inside the cell can be
approximated by: uf = uC + (xP − xC )· ∇uC . The numerical scheme used
to compute the cell-centered gradient (∇uC ) can be defined by the user un-
der keyword gradExternalForcing in dictionary fvSchemes. In general, we
recommend Gauss linear for hex-dominated grids and leastSquares in the
remaining cases.
Gradient
y =3 uf
3
2 Barycentric
1 1 1 Weights
-1 1 x
y =1
2 2 Profile at
y=0
3 3 3
uf
3 Gradient
y =-1 2 2
2
Barycentric
Weights
1 1 1 -1 1 y
Profile at
y =-3 x=0
x =-3 x =-1 x =1 x =3
The previous case has shown that the numerical interpolation methods available
are prone to errors in certain conditions and care should be taken in their use. In
the specific case presented above, the issues with the two methods could have been
reduced by refining the grid in the y-direction, such that the molecules spanned
CHAPTER 4. Overview of rheoTool 64
more than one cell in that direction, and by allowing the molecules to travel at least
one entire cell in the x -direction. Alternatively, the use of an unstructured grid
could possibly solve the issues in that particular case. The message that we want
to convey to the reader/user is to use refined meshes in BDS cases, even though
the interpolation methods have sub-cell resolution. Moreover, grid-dependency
studies are also advisable. In any generic situation, if Analytical interpolation can
be used, then this should always be the preferred method, since it is not prone to
such errors. If not possible, then method BarycentricWeights should be preferred
over Gradient.
(3) k
Model (1)
TypeName (2)
FSi fi (4,5) k
Jij, a (j6=i, j∈gi )
N
P
Hookean Hookean H xij – –
j∈gi
N h i h i
|xij | xij D S,k 2 H0 δ2 δ2
2 1 1
rki − r∗i − ∆t kT 1 1 1
P
Marko-Siggia MarkoSiggia 3 Hl l − 4 + 4(1−|xij |/l)2 |xij | Fi 3 d 4(d−1)2
+d− 4 d2
−1 − d 1− 2(d−1)3
j∈gi
N h
3−(|xij |/l)2
i h i
D S,k H 0 2δ 2 d2 −3 d2 −3
H
rki − r∗i − ∆t kT
P
Cohen Padé CohenPade 3 1−(|xij |/l)2
xij Fi 3 (d2 −1) d2 −1
−1− 2δ 2
j∈gi
N h i
D S,k 2δ 2
1
rki − r∗i − ∆t kT H0 1
P
FENE FENE H x
1−(|xij |/l)2 ij
Fi d2 −1
− (d2 −1)2
j∈gi
(1)
Corresponds to the name entry identifying the model in the source code.
(2) 3kT N
k,s
H= l2 .
(3)
Vectorial function used in the Newton-Raphson method (see Section 3.7.4). Superscript k denotes the iteration index.
(4)
The expressions in this column are for the off-diagonal ij (j 6= i, j ∈ gi ) elements of the Jacobian matrix of the Cartesian component a = x, y, z of fk . All the
N
k k
P
off-diagonal elements for which j ∈/ gi are zero. The diagonal elements are obtained from the sum of the off-diagonal elements, Jii,a = 1− Jim,a . This
m=1,m6=i
symmetric matrix is used in the Newton-Raphson method (see Section 3.7.4).
(5)
H 0 = H∆t kT
D
, d = |xkij |/l = |rkj − rki |/l, δ = (rj,a
k k
− ri,a )/l
65
CHAPTER 4. Overview of rheoTool 66
In general, the FENE and Cohen Padé models are stiffer to solve than the
Marko-Siggia model, thus requiring smaller time-steps.
An explicit time integration scheme (TypeName = explicit) is available for all
spring models. For the bounded spring models, the semi-implicit scheme (Type-
Name = semiImplicit) described in Section 3.7.4 is also available. The auxiliary
functions used in the semi-implicit scheme are specified in Table 4.3 for all bounded
models.
As mentioned in Section 3.7.4, several methods can be used to solve the lin-
ear system of equations (3.53). Four (direct) methods are available in library
BDmolecule:
TypeName – Description
LLT – performs a standard Cholesky decomposition of matrix Jka , using the llt()
function of Eigen [7]. It should not be used for tethered molecules.
TDMA – implements Thomas algorithm for a tridiagonal matrix. It can be used for
open, linear (no branches), not tethered molecules. Better performance than
QR is achieved for a large number of beads per molecule (it only visits
elements i -1, i, i +1).
The reader may be questioning the need and utility of so many methods. The-
oretically, some methods should be faster than others, but possibly more stringent
in their application (only tridiagonal matrices, diagonal-dominant matrices, etc.).
In practice, we observe that the QR method is usually the best compromise. The
TDMA method presents a CPU time typically equal or smaller than QR, but has
restrictions in its application (see above). All the methods are implemented us-
ing full-matrices, notwithstanding the fact that most of the off-diagonal matrix
elements are zero. Sparse matrices and sparse matrix (iterative/direct) solvers
might be considered in a future version if memory overhead starts to be an issue
of concern (say > 500 beads per molecule).
v5.0, where barycentric coordinates for the beads’ position have been general-
ized. It is out the scope of this user guide to explain the particle-tracking algo-
rithm. However, it is worth mentioning that the drag force (and only that one)
is re-interpolated inside the same time-step each time a particle (bead) crosses an
internal face of the mesh.
Once a bead hits a wall, it is repositioned at the location where the wall is
crossed or, optionally, at a given distance from that point, in the wall-normal di-
rection. This artificial repulsion aims to reproduce the finite size of the beads, that
would never let them to approach infinitely close to a surface in a real situation.
If a bead hits a patch, it will simply leave the mesh through such patch and
the corresponding molecule is deleted. If a 2D flow is simulated, we recommend to
make the empty direction long enough in order to avoid the unnecessary contact
of the empty boundaries with the beads. The contact of a bead with any other
boundary type will result in undefined behavior, although we expect the molecule
to be deleted in most of the cases. The corollary is to create meshes having only
patch, wall and empty boundary types for use in Brownian dynamics simulations.
Other types can be used as long as the beads do not enter the cells owning those
boundary faces.
It is well-known that the hydrodynamic interactions near a wall are different
than in the bulk. However, the current version of BDmolecule does not make any
correction at the walls. Therefore, care should be taken when simulating confined
problems, or when the molecule-wall contact is recurrent.
IDs.txt: contains three columns, where the first one is for the molecule name/ID,
the second is for the number of beads in the molecule and the third is for
the molecule’s group ID. This file includes such information for all molecules
starting the simulation, regardless of whether the molecule is deleted at some
time during the simulation (e.g. if one of the springs becomes overstretched).
Each row represents a molecule.
stretch.txt: the first column of this file is for the physical time and each of the remaining
columns represents the stretch of a molecule. There is a direct and sequential
one to one relation between the rows of IDs.txt and the columns of st
retch.txt (excluding the first column for time). We define stretch as
the maximum inter-bead distance in a molecule: stretch = max(|rj − ri |),
i, j = {1..N }. If a given molecule is deleted at some point of the simulation,
then the column corresponding to that molecule will be filled with zeros
starting from the row corresponding to that time. Since the length of a
molecule can never be zero, these entries can be used as flags to detect the
molecules that have been deleted during the simulation.
X.txt: the first column of this file is for the physical time and each of the following
triplets of columns represents the x-, y- and z-coordinates of the molecule
N
center of mass, xcm = N1
P
ri . Again, there is a direct and sequential one
i=1
to one correspondence between the rows of IDs.txt and each triplet of
columns of X.txt (excluding the first column, for time).
4.3.9 Limitations
Some of the main limitations of the Brownian dynamics library have been pre-
viously described. However, to make them clear and to avoid any misuse of the
library, they are summarized below:
• only patch, wall and empty boundary types can be present in a mesh.
• hydrodynamic interactions near surfaces are not corrected and the algorithm
for beads reflection at the walls is simplistic (accuracy in confined-flow con-
ditions can be compromised).
Of course, all these limitations can be seen as points to improve in future work.
Although the non-parallelization of the code is pointed out as an issue, we believe it
is not the most important one. Indeed, we expect the solver to be used essentially
with a steady external forcing. In such cases, if the code was parallelized we would
distribute the ensemble of M molecules by P available processors (M/P molecules
per processor), in a single run (single mesh, single fields, etc.). Since this is not
CHAPTER 4. Overview of rheoTool 69
possible, we can prepare P cases, with M/P molecules per case and run each one in
a processor. This is valid because all the molecules are similar objects and need not
to communicate between them. The two methods should perform closely in terms
of CPU time, although the non-parallel one consumes more memory (allocation
of the mesh P times vs 1 time). For the cases with transient external forcing,
parallelization would be indisputably advantageous regarding CPU time.
4.4 Solvers
The solvers available in rheoTool are summarized in Table 4.4. Each one is dis-
cussed in detail in the following sections.
Name Description
4.4.1 rheoFoam
The solver rheoFoam implements the transient, incompressible Navier-Stokes equa-
tions for single-phase flows of GNF or viscoelastic fluids. Figure 4.2 displays the
solving sequence.
Start
createField.H
Time loop: t = t + Δt L1
L2
Inner loop: i = i + 1
UEqn.H
Non-orthogonality L3
corrector
loop
Solve the pressure S3
(continuity) equation
pEqn.H
pEqn.H
rheoFoam.C
CEqn.H
End
Figure 4.2: Solving sequence of rheoFoam (steps related with mesh motion are
omitted for simplicity).
CHAPTER 4. Overview of rheoTool 71
The solver has three main loops: L1, which is advancing the time; L2, which is
an inner-loop used to converge the solution at each time-step; and L3, a loop which
can be enabled for non-orthogonal grids, in order to update (inside each time-step
and each inner-iteration) the explicit correction of the pressure Laplacian, avoiding
stability problems and reducing the error in transient computations. More than
understanding each step identified in Fig. 4.2, we want to help the reader to identify
them in the source code and relate them with the theory presented in the previous
Chapter. With this purpose in mind, we will do a tour to the source code of
the solver and the most important points will be discussed, skipping the lines not
essential to understand the algorithm.
The solver rheoFoam is composed of one main file (rheoFoam.C) and other
associated header files, among which (createFields.H,createPPutil.H,
UEqn.H,pEqn.H,CEqn.H), that can be found in directory src/solvers/rh
eoFoam/. All the header files are included from the main .C file. We will start
by digging into rheoFoam.C, whose source code is displayed in Listing 4.5.
• line 12: this # include allows the solver to access the models defined in the
constitutiveEquations library.
• lines 17-29: fields, variables, controls and the mesh are created (step S1 of
Fig. 4.2). Lines 23 and 27 point to the header files createFields.H and
createPPutil.H, respectively, located in the same directory as rheoFo
am.C. The later is responsible for the creation of post-processing utilities.
• lines 36-93: represents the time loop, i.e., L1 of Fig. 4.2. The solver will
keep running until the final specified time is reached or once the residuals of
the solved variables drop below some tolerance (this dual criterion is specific
of class simpleControl ).
• lines 66-85: this is the inner loop, L2, of Fig. 4.2. Inside this loop, all
the conservation equations are solved nInIter times inside the same time-
step. This reduces the explicitness of the method, which exists, for example,
in the non-linear convective term of the momentum equation, in the both-
sides-diffusion technique and in several terms of the constitutive equation
(for a given equation, only the terms introduced through a fvm:: operator
are implicit). Furthermore, these iterations also strengthen the coupling
between velocity and pressure.
CHAPTER 4. Overview of rheoTool 72
• lines 44: this function executes the mesh motion if the mesh is dynamic. For
a static mesh, this function has no practical effects.
• lines 46-63: these lines are only executed for dynamic meshes and include
an optional correction of fluxes (line 53). Note that the term ub of Eq. (3.24)
is subtracted from u in line 57.
• line 73: the momentum equation is solved. The header file UEqn.H (Listing
4.6) will be explored later.
• line 74: the pressure equation is solved. The header file pEqn.H (Listing
4.7) will be explored later.
• line 78: function correct() of the constitutive model is called. As seen before,
this function updates variable τ by solving the constitutive equation(s).
• line 81-84: the equation for a passive scalar is optionally solved, depending
on a user-defined selection (more details in Section 5.1.1). The header file
CEqn.H (Listing 4.8) will be explored later.
• lines 87: this is where the post-processing utilities are evaluated, if any has
been selected by the user.
1 #include "fvCFD.H"
#include "IFstream.H"
3 #include "OFstream.H"
#include "simpleControl.H"
5 #include "fvOptions.H"
#include "extrapolatedCalculatedFvPatchField.H"
7 #include "dynamicFvMesh.H"
#include "CorrectPhi.H"
9 #include "adjustCorrPhi.H"
11 #include "ppUtilInterface.H"
#include "constitutiveModel.H"
13 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * //
31
Info<< "\nStarting time loop\n" << endl;
33
// --- Time loop ---
35
while (simple.loop())
37 {
#include "readControls.H"
39 #include "CourantNo.H"
#include "setDeltaT.H"
41
Info<< "Time = " << runTime.timeName() << nl << endl;
43
mesh.update();
45
if (mesh.changing())
47 {
// Calculate absolute flux from the mapped surface
velocity
49 phi = mesh.Sf() & Uf();
51 if (correctPhi)
{
53 #include "correctPhi.H"
}
55
// Make the flux relative to the mesh motion
57 fvc::makeRelative(phi, U);
59 if (checkMeshCourantNo)
{
61 #include "meshCourantNo.H"
}
63 }
}
85 }
87 postProc.update();
runTime.write();
89
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << "
s"
91 << " ClockTime = " << runTime.elapsedClockTime() << "
s"
<< nl << endl;
93 }
97 return 0;
}
Listing 4.5: Source code of rheoFoam.C.
The source code in file createFields.H will not be discussed, since it mainly
contains standard declaration/initialization of fields and variables. However, it is
worth mentioning the creation of a constitutiveModel object, named constEq, which
stores the data of the constitutive equation selected.
According to the SIMPLEC algorithm, the momentum equation is first solved
to obtain an estimated velocity field, using the pressure field of the previous inner-
iteration or time-step. Then, we solve for the pressure field enforcing continuity.
Finally, using the correct pressure field, the previously estimated velocity is cor-
rected, both on faces and cell centers. This is what is being executed in UEqn.H
and pEqn.H. We follow the discussion by analyzing the content of UEqn.H, whose
code is displayed in Listing 4.6.
• lines 5-13: this is where the momentum equation is built. We can iden-
tify the transient term in line 7, an eventual acceleration term arising from
a non-inertial reference frame in line 8, the convective term in line 9, an
extra momentum source term in line 11 (term f in Eqs. 3.5 and 3.18) and
0
the extra-stress divergence in line 12 (∇· τ ), containing both the solvent
and the polymeric contributions (remember the output of divTau() function,
analyzed in Section 4.1.4).
• line 21: the momentum equation is solved considering the pressure gradi-
ent contribution, where the pressure field from the last time-step or inner-
iteration is used. This term was not added before in order for operator H
obtained from the momentum equation to be free of such contribution, as
discussed in Section 3.3.1. This is required to avoid the onset of checkerboard
fields, since a traditional Rhie-Chow interpolation is not used (cf. Ref. [2]).
// Momentum predictor
2
MRF.correctBoundaryVelocity(U);
4
CHAPTER 4. Overview of rheoTool 75
tmp<fvVectorMatrix> tUEqn
6 (
fvm::ddt(U)
8 + MRF.DDt(U)
+ fvm::div(phi, U)
10 ==
fvOptions(U)
12 + constEq.divTau(U)
);
14
fvVectorMatrix& UEqn = tUEqn.ref();
16
UEqn.relax();
18
fvOptions.constrain(UEqn);
20
solve(UEqn == -fvc::grad(p));
22
fvOptions.correct(U);
Listing 4.6: Source code of UEqn.H.
After having a guessed (non-conservative) velocity field, we will see how it is
used inside pEqn.H (Listing 4.7):
• lines 1-32: variables required to solve the pressure equation (Eq. 3.15) are
assembled. The sequence of steps can be easily understood, keeping in mind
that UEqn().A() retrieves diagonal coefficients (aP ) and that UEqn().H()
and UEqn().H1() stand for operators H and H1 , respectively. As previously
discussed in Section 3.3.1, pressure gradient terms entering the definition of
face fluxes (line 26) are directly evaluated on cell faces to avoid checkerboard
fields. Also, in line 8/11 there is the addition of the corrective term for
time-step dependency, described in Section 3.3.1.
For a number of cases, doing 2-3 non-orthogonal iterations keeps the solver
stable, without the need of under-relaxing the pressure. Even if only one
non-orthogonal correction is performed, the Laplacian term should still be
discretized with the corrective term to keep the accuracy in non-orthogonal
meshes.
• lines 39,44: the pressure equation (Eq. 3.15) is assembled (line 39) and
solved (line 44).
• line 48: this is the equation which corrects the face fluxes (Eq. 3.16 interpo-
lated to the faces). Again, pressure gradient terms are directly evaluated on
cell faces: the snGrad() operator in line 26, when building phiHbyA, and the
one coming from the laplacian() operator in line 39, from which the flux()
operator is derived.
• line 57: this is the equation which corrects the cell-centered velocity field
(Eq. 3.16).
• line 69: the term ub of Eq. (3.24) is subtracted from u (the operation is on
face fluxes, variable phi ).
1 volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
3 surfaceScalarField phiHbyA
(
5 "phiHbyA",
mesh.changing() == true ?
7 fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, Uf())
9 :
fvc::flux(HbyA)
11 + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
13
MRF.makeRelative(phiHbyA);
15
if (p.needReference())
17 {
fvc::makeRelative(phiHbyA, U);
19 adjustCorrPhi(phiHbyA, U, p);
fvc::makeAbsolute(phiHbyA, U);
21 }
23 tmp<volScalarField> rAtU(rAU);
29 tUEqn.clear();
57 U = HbyA - rAtU*fvc::grad(p);
U.correctBoundaryConditions();
59 fvOptions.correct(U);
61 if (mesh.changing())
{
63 Uf() = fvc::interpolate(U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());
65 Uf() += n*(phi/mesh.magSf() - (n & Uf()));
}
67
// Make the fluxes relative to the mesh motion
69 fvc::makeRelative(phi, U);
Listing 4.7: Source code of pEqn.H.
After pEqn.H is executed, both the pressure and the face fluxes are continuity-
compliant, but not the cell-centered velocity field. The conservative fluxes can now
be used to solve any transport equation. In rheoFoam, we offer the possibility to
solve a transport equation for a passive scalar. The governing equation, included
in file CEqn.H, is simply a convection-diffusion transport equation, as can be seen
in lines 5-11 of Listing 4.8.
1 // Transport of passive scalar
3 dimensionedScalar D_ = cttProperties.subDict("
passiveScalarProperties").lookup("D");
5 fvScalarMatrix CEqn
(
CHAPTER 4. Overview of rheoTool 78
7 fvm::ddt(C)
+ fvm::div(phi, C)
9 ==
fvc::laplacian(D_,C)
11 );
13 CEqn.relax();
CEqn.solve();
15
if (U.time().outputTime())
17 {
C.write();
19 }
Listing 4.8: Source code of CEqn.H.
4.4.2 rheoTestFoam
The main purpose of solver rheoTestFoam is to evaluate the behavior of the con-
stitutive models for a user-defined ∇u tensor. At the same time, it can also be
envisaged as a basic debugging tool to check for the correct implementation of the
constitutive models, since an analytical or semi-analytical solution usually exists,
which can be used for comparison.
Shortly, rheoTestFoam solves for the solvent and polymeric constitutive equa-
tions, Eqs. (3.3) and (3.4), respectively, for a prescribed ∇u tensor, assuming
homogeneous flow conditions (∇· τ = 0; u· ∇τ = 0). Since there are no approx-
imations related with spatial discretization, the resulting steady-state solution
0
τ = τ + τs is exact, and OpenFOAM
R
is simply acting as a nonlinear matrix
solver. To obtain unsteady solutions, the temporal discretization introduces nu-
merical errors during the transient period, which can be reduced using a small
time-step.
The computational domain used with this solver is composed of a single cell: a
cube with unitary edge length (1 m). The boundary conditions for u are internally
manipulated inside the code, in order to get the tensor ∇u defined by the user,
Fig. 4.3. Thus, the default mesh and boundary conditions should not
be changed by the user when working with rheoTestFoam. We note that the
following definition holds in this guide (and in OpenFOAM
R
, in general):
∂u ∂v ∂w
∂x ∂x ∂x
∂uj ∂u ∂v ∂w
(∇u)ij = = ∂y (4.27)
∂xi ∂y ∂y
∂u ∂v ∂w
∂z ∂z ∂z
• ramp mode: the user defines a list of ∇u tensors and the solver will retrieve
the steady solution for each entry. In this mode, the solver automatically
selects the ideal time-step value to be used and the steady-state is also auto-
matically detected (either the relative variation of the extra-stress magnitude
CHAPTER 4. Overview of rheoTool 79
z+
z x-
y- x y+
x+
z-
u x u y x u z x u x u y x u x
u x 1 x , 2 u x 1 x , 2 , 3 z
2 , 3 x 2 2 x 2
x 2 x x 2 x
u y u y y u z y u y u y y u z y
u y 1 x , 2
2 , 3 y 2 u y 1 x , 2
2 , 3 y
y 2 y y 2 y 2
u z u y z u z u z u y z u z
u z 1 x , 2 , 3 z u z 1 x , 2 , 3 z
2 z 2 2 z 2
z 2 z z 2 z
Figure 4.3: Boundary conditions manipulation in the single-cell mesh used with
rheoTestFoam. The constants currently used to represent the cell-centered velocity
are κ1 = κ2 = κ3 = 0, although any other values could be used. The edge length
of the cubic cell is set to δx = δy = δz = 1m.
drops below 10-8 , or after a predefined number of time-steps has been ex-
ceeded – this last condition is used to avoid infinite loops).
• transient mode: the user defines one single ∇u tensor and the solver will
return the evolution over time of the monitored variables. In this case,
both the time-step and the end time are controlled by the user. This mode
allows to determine the transient material functions of the constitutive model
selected.
4.4.3 rheoInterFoam
4.4.4 rheoEFoam
The solver rheoEFoam is the extension of rheoFoam to electrically-driven flows –
note the additional E in the solver name, which points to its E lectric component.
Thus, with no surprise, both solvers share the same basic structure and only
minor differences exist between both at the code level. For this reason, we will
not discuss again the common parts. Instead, we refer the reader to Section
4.4.1 to eventually recall the structure of rheoFoam and in this Section we only
highlight the main differences introduced in rheoEFoam. Note that rheoEFoam
allows the combination of both electrically- and pressure-driven flows with no
restrictions. Pure pressure-driven flows can be also simulated with rheoEFoam, for
which the solver becomes functionally equivalent to rheoFoam, although directly
using rheoFoam is the more efficient option in these situations.
CHAPTER 4. Overview of rheoTool 81
Starting by the file createFields.H, there is an extra line for the creation
of the electric model, as shown in Listing 4.9.
1 // Create the electric model
EHDEKModel elecM(phi);
Listing 4.9: Line in the file createFields.H of rheoEFoam, where the
electric model is created.
In file rheoEFoam.C, the header # include ”EDFModel.H” has been added in
order to enable the use of the EDFModels library, whose path had also to be added
to file Make/options. The other change is in the inner-iteration loop, which now
includes solving the electric-related equations (line 21, Listing 4.10), calling the
function correct() of the electric model (recall an example of that function in lines
188-281 of Listing 4.3). The user may also choose not to solve the equations
governing the fluid flow (see the if condition in line 7 of Listing 4.10). This is
useful when only the electric component of the problem is of interest.
// --- Inner loop iterations ---
2 for (int i=0; i<nInIter; i++)
{
4
Info << "Inner iteration: " << i << nl << endl;
6
if (solveFluid)
8 {
// --- Pressure-velocity SIMPLEC corrector
10 {
// ---- Solve U and p ----
12 #include "UEqn.H"
#include "pEqn.H"
14 }
1 tmp<fvVectorMatrix> tUEqn
(
3 fvm::ddt(U)
+ MRF.DDt(U)
5 + fvm::div(phi, U)
==
7 fvOptions(U)
+ constEq.divTau(U)
9 + elecM.Fe()/constEq.rho()
);
Listing 4.11: Momentum-balance equation in rheoEFoam (source code:
UEqn.H).
4.4.5 rheoBDFoam
Solver rheoBDFoam is intended for Brownian dynamics simulations in generic
meshes. In a glance, the source code of rheoBDFoam is simply the one of
rheoEFoam, to which we added an object of class sPCloudInterface, that is up-
dated over time. Therefore, we recommend the reader to first take a look at
Sections 4.3 and 4.4.4. The pieces of code specifically devoted to Brownian
dynamics are mainly contained in header files createLagrangianFields.H
and moleculesEqns.H. In createLagrangianFields.H, a sPCloudInter-
face object is created and named molecules. This header file is called once in
rheoBDSFoam.C, at the beginning of the simulation. On the other hand, header
file moleculesEqns.H is included in rheoBDSFoam.C each iteration of the
main time loop. Let’s take a look on its source code (Listing 4.12):
• lines 8-14: these lines enclose a special loop, which first divides the Eulerian
time-step by nSubCycles, and then performs the operations inside it. The
old time is incremented by this new time-step in each iteration of the loop,
such that an entire original (Eulerian) time-step is elapsed upon completion
of the loop. This procedure allows to use the original Eulerian time-step
everywhere outside the loop (e.g. in solving continuum equations), while a
lower time-step can be used for the operations enclosed in the loop, as long
as nSubCycles > 1 (this variable should be a positive integer different from
0; line 1). The operation executed inside the loop is defined at line 10, which
basically consists in a call to function update() of object molecules. The tasks
carried out by this function have been discussed in Section 4.3.2 and result
essentially in the update of the molecules’ state. The return value of the
function is true if at least one valid molecule remains in the computational
mesh, otherwise it is false and forces an early exit from the loop (lines 12-13).
• lines 20-24: the main time loop of the solver is exited if no valid molecules
remain inside the mesh.
face and the cell owning that face, respectively, and dPf is the vector connecting
their geometrical centers.
Description: linear extrapolation of each field component from boundary cells
to boundary faces. Shortly, this boundary condition starts by computing the gradi-
ent of each component at the center of the cell owning the boundary face (using the
previous iteration/time-step known values on the boundary face). Then, with both
the value and the gradient of each component at those locations, the components
at the boundary faces are estimated by linear extrapolation. The discretization
scheme to compute the gradients enrolled in the process is run time selectable
and can be adjusted in dictionary fvSchemes, through the entry linExtrapGrad
followed by the selected scheme, in the gradSchemes subDict. In general, using
linear extrapolation for the polymeric extra-stress tensor on walls should be pre-
ferred in relation to a zero-gradient boundary condition, which has a lower order
of accuracy [2].
Since version 2.0, an optional second-order accurate linear regression (see [32])
can be selected by adding the entry useRegression true to the dictionary defining
the BC, i.e., below keyword type, for example. If not present, the solver will execute
by default the linear extrapolation defined above.
4.5.2 navierSlip
TypeName: navierSlip
Type: fixed-value (vector).
Formula:
0
−k τ0 m τw,f ,model =nonLinearN avierSlip
nl w,f 0
ut+∆t
f = (1−URF)utf +URF |τw,f |
−α 1 − β|τ0 | τ0 ,model =slipT T
w,f w,f
where URF is the under-relaxation factor (0 < U RF ≤ 1), utf is the boundary
velocity at the previous time-step, and knl , m, α and β are input model-dependent
parameters (the model name should be also defined). 0 The wall
stress vector is
0 0 0
tangent to the wall and is given by τw,f = τf · nf − τf · nf · nf nf with τf corre-
sponding to the total (solvent and polymeric) extra-stress tensor at the wall.
Description: fully-explicit implementation of slip models according to [33].
Supports both two-phase flows and moving meshes. In the latter case, the moving
wall velocity is added to the slip component computed as above. Lower values
of URF increase the numerical stability, but are not adequate for time-dependent
flows.
4.5.3 zeroIonicFlux
TypeName: zeroIonicFlux
Type: fixed-gradient (scalar).
CHAPTER 4. Overview of rheoTool 85
4.5.4 boltzmannEquilibrium
TypeName: boltzmannEquilibrium
Type: fixed-value (scalar).
Formula: ci,f = ci,0 exp ez
kT
i
(ψf − ψ0 ) , where ci,0 is the reference concentration
(user-defined) at which the intrinsic electric potential is ψ0 (user-defined; see the
definition of the other variables in Section 3.6.1). A common choice is to set ci,0
as the bulk concentration of ions and ψ0 = 0.
Description: ionic concentration derived from the assumption of Boltzmann
equilibrium near the patch. This boundary condition is intended to be used when
the electric potential is split, in which case only the intrinsic potential is used
in the formula. This boundary condition does not guarantee the zero-flux of a
given specie in all the situations, and, in general, it is not accurate for transient
simulations.
4.5.5 inducedPotential
TypeName: inducedPotential
Type: fixed-value (scalar).
Nf
Formula: ψf = −φExt,f + ψFix + NPf 1
P
φExt,f |Sf |, where ψFix is the bias
|Sf | f=1
f=1
voltage (user-defined) of the patch, i.e., the electric potential of the surface in the
absence of an external electric field.
Description: intrinsic potential induced in a conducting surface placed over
an electric field and having a bias voltage [34]. The last term in the formula
represents the area-averaged external electric potential over the surface. This
boundary condition can be used with the Poisson-Boltzmann and Debye-Hückel
models, under the potentials splitting approach.
4.5.6 slipSmoluchowski
TypeName: slipSmoluchowski
Type: fixed-value (vector).
CHAPTER 4. Overview of rheoTool 86
4.5.7 slipSigmaDependent
TypeName: slipSigmaDependent
Type: fixed-value (vector).
m
Formula: uSch,f = µ0 σσ0f (−∇φExt,f ), as defined in Section 3.6.6. Parame-
ters µ0 , σ0 and m are user-defined.
Description: slip velocity derived from the Helmholtz-Smoluchowski theory
for a space-variable conductivity field (Section 3.6.6). This boundary condition
can be used with the Ohmic model.
Hf 1 1 ∗
uf · n = 0 ⇒ (∇p)f · n = (aP − H1 ) + − (∇p )f · n (4.28)
aP aP − H 1 aP
Eq. (4.28) is generic since it does not depend on a specific form of the mo-
mentum equation. It can also be seen that if both sides of Eq. (4.28) are forced
to be zero – the zero-gradient approach –, the equality is still satisfied and the
no-penetration condition is still valid. The advantage in using Eq. (4.28) is that
the resulting normal pressure gradient is going to effectively balance the remaining
local forces in the momentum equation. In practice, the difference between using
Eq. (4.28) or the zero-gradient approximation is only noticeable when the stresses
at the wall are significant (high HaP
f
), as for example in EDFs, where a strong elec-
tric force normal to the wall may exists. In the tutorials provided with rheoTool ,
the zero-gradient approximation is frequently used.
Since OpenFOAM
R
version 4.0, the boundary equation embodied by Eq. (4.28)
is available by default under the name fixedFluxExtrapolatedPressure. However,
CHAPTER 4. Overview of rheoTool 87
this boundary condition may not be usable in some cases, as for example in cases
with a zero velocity assigned to all boundaries. In such cases, there is a conflict
created by function adjustPhi() when it attempts to artificially adjust the fluxes.
This issue can be avoided by commenting the line where function adjustPhi() is
called in the code (typically inside file pEqn.H), provided the user is sure that
the set of boundary conditions under use verifies continuity. This limitation in the
boundary condition might be eventually solved in future releases of OpenFOAM
R
.
Importantly, in multiphase flows with non-zero surface tension and/or gravity
effects, the zero-gradient condition for pressure is generally incorrect. Indeed, for
such flows it is usual to compute the fluxes in two steps in OpenFOAM
R
, i.e., the
Hf
term aP in Eq. (4.28) is built in two steps. The last step adds the flux contribution
from surface-tension and gravity without enforcing a null contribution at the wall.
Thus, in general the right hand side of Eq. (4.28) is non-null in these cases and sim-
ply setting (∇p)f · n = 0 creates a local mass imbalance. The boundary conditions
fixedFluxPressure (available in OpenFOAM
R
and foam-extend) and fixedFluxEx-
trapolatedPressure (only available in OpenFOAM
R
starting from version 4.0) solve
this issue by equating the normal pressure gradient to this extra flux, and they
should be employed in such situations. In addition, for rheoTool versions run-
ning foam-extend, when the SIMPLEC algorithm is selected for pressure-velocity
coupling, the fixedFluxPressure BC requires that the entry Dp is defined and set
to ”rAtU ” (check the tutorials provided). This is not needed for the PIMPLE
coupling algorithm, or when using OpenFOAM
R
v4.x.
4.6 Utilities
4.6.1 GaussDefCmpw schemes for convective terms
The component-wise and deferred correction handling of HRSs, described in Sec-
tion 3.4, is included as a library in rheoTool . If the installation procedure pre-
sented in Chapter 2 has been followed, this new class of schemes will only be
available when using the family of solvers provided with rheoTool . However, there
are several ways to make the schemes available to any solver of OpenFOAM
R
.
One option (not requiring compilation) is to include this library (libgauss
DefCmpwConvectionSchemes.so) as a lib entry of controlDict in the
case directory. Another option is to compile the class inside library finiteVol-
ume, which is included by most OpenFOAM
R
solvers. To access this class from
a specific solver, lgaussDefCmpwConvectionSchemes should be added to
the Make/options file of that solver, along with its path. We note that the
component-wise and deferred correction handling of HRSs improved significantly
the stability of viscoelastic fluid flow simulations [2], but its performance and ad-
vantage when used in other type of flows need to be tested (by no way we argue
that this is a magic bullet for all purposes).
The new group of HRSs is accessible from class GaussDefCmpw and its use
is similar to the standard HRSs of OpenFOAM
R
. For example, the CUBISTA
scheme can be used by simply defining in dictionary fvSchemes: GaussDefCmpw
cubista; in front of the divergence term being discretized (remember that keywords
CHAPTER 4. Overview of rheoTool 88
in OpenFOAM
R
are case-sensitive). To obtain a list of all the schemes available,
simply type GaussDefCmpw ; without any additional argument and you will ob-
tain all the possibilities, listed in Table 4.5. There is a scheme named none, which
corresponds to removing the convective term from the equation being discretized.
Note that all the limiters implemented in class GaussDefCmpw are totally inde-
pendent from the already existing limiters of OpenFOAM
R
and all are defined
in file limiters.H. For example, you will have now GaussDefCmpw minmod
and Gauss Minmod, which are two different schemes, or, actually, two different
implementations of the same high-resolution scheme.
Table 4.5: Available High-Resolution schemes for convective terms in class gauss-
DefCmpw. The schemes are defined using the NWF approach (Eq. 3.20).
1 2
Scheme TypeName Equation
[1, 0] φ̃C ≤ 0 ∨ φ̃C ≥ 1
[3, 0] 0 < φ̃C < 1/6
SMART smart [α, β] =
[3/4, 3/8] 1/6 ≤ φ̃C ≤ 5/6
[0, 1] 5/6 < φ̃C < 1
[1, 0] φ̃C ≤ 0 ∨ φ̃C ≥ 1
[2, 0] 0 < φ̃C < 3/10
WACEB waceb [α, β] =
[3/4, 3/8] 3/10 ≤ φ̃C ≤ 5/6
[0, 1] 5/6 < φ̃C < 1
[1, 0] φ̃C ≤ 0 ∨ φ̃C ≥ 1
[1/2, 1/2] 0 < φ̃C < 1/2
SUPERBEE superbee [α, β] =
[3/2, 0] 1/2 ≤ φ̃C ≤ 2/3
[0, 1] 2/3 < φ̃C < 1
3
no convection none –
1
Corresponds to the name entry identifying the scheme in the source code.
2
See Eq. (3.20).
3
When this option is used, the convective term is deleted.
ogy between the equations presented in Section 3.4 and the source code in files
gaussDefCmpwConvectionScheme.C and limiters.H. Nevertheless, for
documentation purposes, we summarize next the operations being executed by
each member function of class GaussDefCmpw :
• lims(): this function retrieves three variables: alpha, is a list containing α for
each interval of the function defined in Eq. (3.20); beta is a list containing β
for each interval of the function defined in Eq. (3.20); bounds is a list of φeC
values, for which there is a change of branch in the function defined in Eq.
(3.20). Thus, function lims() defines Eq. (3.20) for the selected scheme.
• fvmDiv(): this function also exists for Gauss schemes and returns the matrix
of coefficients and the source term resulting from the discretization of the
implicit convective operator fvm::div(). Function phifDefC() is called from
here, whenever the selected scheme is different from upwind or none. Note
that both Gauss and GaussDefCmpw classes implement the upwind scheme
in the same way – it is the only scheme for which this happens.
• flux(): returns the field interpolated on face centers multiplied by the flux
on each face (phi ).
• copy and paste the folder of an already existing ppUtil and give it a new
name of your choice. Delete the .dep file in that folder.
• find & replace the old name of the ppUtil by the name that you gave to the
folder. This should be done for both the .C and .H files (there are several
ways to do it automatically, for example, Ctrl + H in gedit).
• add the source file that was just created to the list of files for compilation in
the Make/file of the library.
• run the Allwmake script in src/ to compile, and it should be ready to use.
4 functions
(
6
ciMonitor
8 {
funcType calcBalance;
10 enabled true;
evaluateInterval 100;
12 }
14 jMonitor
{
16 funcType calcJpatch;
ListOfPatches
18 (
"cylinder"
20 );
enabled true;
22 evaluateInterval 100;
}
24
CHAPTER 4. Overview of rheoTool 91
);
26
}
Listing 4.13: Example of a PostProcessing subDict.
For each ppUtil selected, at least three keywords must be defined:
• funcType: should specify the TypeName of the given ppUtil. To obtain a full
list of all the available utilities, simply insert any random letter.
When a given ppUtil is active and programmed to write some quantity (or
several) of interest, a folder named rheoToolPP/startTimeName/ppUtilN
ame/ is created in the case directory and the output is forwarded to there.
The class ppUtil not only allows to create case-specific post-processing tools,
as it also offers the possibility to build generic post-processing applications. In
Table 4.6 we present some generic ppUtil which are included in rheoTool and that
can be useful in a number of cases.
Table 4.6: General-purpose ppUtil available in rheoTool .
TypeName Description
4.6.3 writeEfield
By default, the solver rheoEFoam does not write the electric field to the time
directories. The purpose of utility writeEfield is to read the electric potential
variable(s) and write the electric field, for each time directory. This means that
this utility can only be called after the simulation has been run. The utility sums
up all the electric potential variables available in the directories: Ψ (psi ), φExt
(phiE ) and/or ψ (psi ).
The utility can be used by typing writeEfield in the terminal, without any
other requirements, except that at least one electric potential variable must exist
in the time directories. In addition, the time directories can not be decomposed
among processors (reconstruct the case if it is decomposed).
Note that the electric field can be also computed from the electric potential
in most of the visualization software, as for example Paraview, since a simple
differentiation operation is required.
4.6.4 initMolecules
The purpose of utility initMolecules is to generate a set of files representing the
molecules that can be read by solver rheoBDFoam. The users can interpret this
utility as a sort of blockMesh application that generates molecules instead of a
computational mesh.
This utility reads its controls from dictionary initMoleculesDict, which
should be located inside folder constant/. An example of such dictionary is
presented in Listing 4.14. We will now analyze its entries. However, we should first
remember that solver rheoBDFoam can handle simultaneously multiple groups of
molecules with different physical properties. As such, both the physical properties
and the beads positions should be defined for each group. In the example of Listing
4.14, there are two groups of molecules defined inside dictionary groups, which were
arbitrarily named G1 (line 3) and G2 (line 27). The user can introduce and define
as many groups as needed (at least one should exist). For each group, the following
data must be present (we will only analyze group G1 in Listing 4.14):
• number of molecules (line 5): this is specified in entry nMolecules and should
be an integer ≥ 1. All the molecules generated inside a group will share the
same physical properties.
CHAPTER 4. Overview of rheoTool 93
• physical properties (lines 7-12): these are the input parameters needed by
the physical model. We can find the diffusion coefficient (D), the beads
radius (a; only used if HI are active), the number of Kuhn steps per spring
(Nks), the exclusion volume parameter (nuEV ; only used if EV forces are
active) and the maximum length of a fully-stretched spring (Ls). See Section
3.7 for details about these parameters.
• spatial distribution and topology of the molecules (lines 15-24): these set-
tings should all be defined inside dictionary spatialDistributionCoeffs (lines
15-24). Currently, the molecules can be only distributed uniformly over a
straight line defined by its extreme points, p0 and p1. For a group with M
molecules, this means that the first bead of the first molecule will be at po-
sition p1 and the first bead of the remaining molecules will be spaced apart
consecutively by vector ∆s = (p1 − p0)/(M − 1) (the first bead of the last
molecule will be at position p2 ). Note that p0 and p1 can be assigned the
same vector positions, in which case all the molecules have their first bead at
the same position. We remember that solver rheoBDFoam does not account
for inter-molecular interactions (for all the effects considered, each molecule
behaves as if it was alone in the surrounding fluid). The last entry of dictio-
nary spatialDistributionCoeffs, named branches() (lines 20-23), controls the
position of the remaining beads in each molecule. The molecules are built
branch-by-branch, thus each line of branches() should fully define a branch
according to a fixed syntax:
(A B C) (d0 d1 d2 ) (E F) (G H)
a given branch can only connect to a branch which has been already entered
in the list. All the examples provided have E = false, which means that
the current bead position is obtained from the previous bead position by
adding vector (d0 d1 d2 )Ls. In these cases, this vector is simply the spring
vector. However, when E = true, each component of the vector provided
is multiplied by a random number between -1 and 1. In either case, it is
user’s responsibility to ensure that the vector provided does not violate the
maximum spring length: d02 + d12 + d22 < 1
1 groups
(
3 G1
{
5 nMolecules 1000;
7 // Physical properties
D 0.065e-12;
9 a .077e-6;
Nks 20;
11 nuEV 1.2e-21;
Ls 2.1e-6;
13
// Spatial distribution
15 spatialDistibutionCoeffs
{
17 p0 (0 0 0);
p1 (0 0.001 0);
19
branches
21 (
(0 0 11) ( 0.2 0 0) (true true)
23 );
}
25 }
27 G2
{
29 nMolecules 700;
31 // Physical properties
D 0.15e-12;
33 a .053e-6;
Nks 15;
35 nuEV 4e-22;
Ls 2.1e-6;
37
// Spatial distribution
39 spatialDistibutionCoeffs
{
41 p0 (0 0 0);
p1 (0.005 0 0);
43
branches
45 (
CHAPTER 4. Overview of rheoTool 95
indices- field of triplets where the first value is the global bead index, the second is
the local bead index inside the molecule to which it belongs, and the third
is the index of the group to which the molecule belongs.
molcID- field of labels corresponding to the index of the molecule which contains the
given bead.
origID- field of labels corresponding to the global bead index. In most of the cases,
this field is equal to the first element of field indices. This field is kept to en-
sure compatibility with the generic particle-tracking engine of OpenFOAM
R
.
origProcID- field of labels corresponding to the index of the processor containing the
bead. Since rheoBDFoam is still not able to run in parallel, this will be
always a field of zeros. Again, this field is kept to ensure compatibility with
the generic particle-tracking engine of OpenFOAM
R
.
positions- list of bead’s positions. This list has a special format defined by
R
OpenFOAM . Each entry is composed by the barycentric coordinates of
the bead (inside parentheses), the mesh cell index containing the bead,
the index of the face owning the tetrahedron containing the bead and the
point index defining that face. For more information, the user is advised to
check the source code of src/lagrangian/basic/particle/partic
le.H and src/lagrangian/basic/particle/particleIO.C pro-
vided with OpenFOAM
R
(version 5.0 or newer).
MoleculesInfo- holds information on the active molecules, as the physical properties of each
group of molecules. It also lists the index, number of beads and the group
that each active molecule belongs to.
springs- a list of triplets defining the springs of the active molecules (similarly to
faces in a mesh). The first element is the first-bead global index, the second
element is the second-bead global index and the third element is the molecule
index. A spring can be unequivocally defined with these elements.
CHAPTER 4. Overview of rheoTool 96
X Spring
Z
Y Z – bead index inside the molecule
Y – branch index
0 1 2 3 4 5 6 branches
0 1 2 3 4 5 6 (
0 y (0 0 7) ( 0.1 0 0 ) (false true)
);
x
Open linear chain
2 9
8
1
7 branches
1 0 (
(0 0 7) ( 0.1 0 0 ) (false true)
0 1 2 3 4 5 6 (0 2 3) ( 0.1 0.1 0 ) (false true)
0 1 2 3 4 5 6 (0 4 4) ( 0.05 -0.05 0 ) (false true)
0 0 2 (2 1 2) ( -0.05 -0.05 0 ) (false true)
10 );
1
14 0 11
y 15 1 3 2
12
3
x 13
0 0
3 0 branches
0 1 (
1 (0 0 3) ( 0.1 -0.1 0 ) (false true)
7
(0 2 2) ( -0.1 -0.1 0 ) (false true)
10 2 9 1 8 0 1 6 2 2 (1 1 2) ( -0.1 0.1 0 ) (false true)
4 (2 1 1) ( 0.1 0.1 0 ) (false false) (0 0)
y (2 1 3) ( -0.1 0 0 ) (false true)
0 5 3 0 1
);
2 4
1
x
Closed chain with branches
The seven files described above are generated upon execution of initMolecules,
and they are automatically written by rheoBDFoam during a simulation, since they
are needed for an eventual restart. Note that all the information in dictionary in
itMoleculesDict is only read by utility initMolecules; changing this dictionary
after having run the utility or just before running solver rheoBDFoam will change
nothing. This is not a valid way to change the molecules’ physical properties
after having created the molecules. Such changes can be done manually in file
MoleculesInfo, although it is generally not recommended (it is preferable to
CHAPTER 4. Overview of rheoTool 97
re-run initMolecules).
4.6.5 averageMolcN
The results obtained with solver rheoBDFoam usually require averaging over sev-
eral molecules due to the random nature of the Brownian term. Utility aver-
ageMolcN averages the molecular length over all the valid molecules registered
at a given time. This utility is mostly useful for homogeneous flows, where the
molecular extension does not depend on the spatial position.
The argument to this utility should be provided in the command line, following
this syntax:
∼$ averageMolcN time
where time is the startTime.
The results for each group of molecules are written to file rheoToolPP/start
Time/moleculesStats/groupName/Stretch_Naverage.txt, where the
first column corresponds to time and the second column is the ensemble averaged
molecular length. Note that untracked/overstretched molecules are not accounted
in the average.
This utility requires (uses) the registry of the molecules’ position and stretch
over time, optionally written and saved by solver rheoBDFoam (see Section 4.3.8).
Therefore, this post-processing utility can only be called after a call to the solver.
4.6.6 averageMolcX
In some cases, it is useful to compute the average evolution of the molecules’ length
over space. Utility averageMolcX is intended for this purpose, since it retrieves the
average molecular length over a line defined by the user. This line can have any
orientation and is divided in a finite number of bins. For example, consider line
connecting point (0,0,0) to point (10,0,0), that is divided in 10 bins. Then, all the
molecules whose center of mass, at a given instant, had its x -coordinate between
x = 0 and x = 1 will contribute to the average in the first bin. The second bin is
between x = 1 and x = 2, and so on for the remaining bins.
The arguments to this utility should be all provided in the command line,
following this syntax:
∼$ averageMolcX time -biased bool -startPoint ”sP” -endPoint
”eP” -nBins nB
where time is the startTime, sP is the position of the first point of the line, eP
is the position of the second point of the line and nB is the number of bins. The
option -biased should receive a boolean and stands for the averaging method
used. The biased average (bool = true) takes all the hits in a given bin and
averages among all with equal weights per hit. This results, in general, in a biased
average, because the molecules with more hits in a given bin (for example the
ones with slower velocity) will have a higher representation in the results. The
unbiased average (bool = f alse; default option if none is specified) first computes
the average length of each molecule inside the bin, collecting all the hits of that
CHAPTER 4. Overview of rheoTool 98
molecule in the bin, and then averages between all the molecules, attributing equal
weights to each one.
The results for each group of molecules are written to file rheoToolPP/start
Time/moleculesStats/groupName/Stretch_Xaverage.txt, where the
first three columns correspond to the x-,y- and z-coordinates of each bin’s center,
the fourth column contains the average molecular length in that bin and the fifth
column contains the number of values used in the average (number of molecules
for unbiased average, or total number of hits for biased average).
This utility requires (uses) the registry of the molecules’ position and stretch
over time, optionally written and saved by solver rheoBDFoam (see Section 4.3.8).
Therefore, this post-processing utility can only be called after a call to the solver.
Chapter 5
Tutorials
The tutorials in this Chapter are mainly intended for learning purposes.
It is not our primary goal to obtain highly accurate results with such
examples, but solely to show how to run the solvers, preferably using
fast-running cases. Higher accuracy can be obtained in all the cases by
increasing the resolution in space and time.
5.1 rheoFoam
5.1.1 General guidelines
Before proceeding, we note that the sequence of operations required to prepare a
case in OpenFOAM
R
does not need to be ordered as presented next (constant/
¸ 0/ ¸ system/). This sequence was organized in such a way to be (hopefully)
logic and easy to follow and execute.
| constant/
Inside folder constant/ there are two main components of the simulation: the
mesh, in folder polyMesh/, and the dictionary constitutiveProperties,
which is a dictionary specific of rheoTool . Since the mesh is an element required
by almost all OpenFOAM
R
solvers, it will not be discussed here and we assume
that a valid mesh already exists in folder polyMesh/.
99
CHAPTER 5. Tutorials 100
10 stabilization coupling;
}
12
passiveScalarProperties
14 {
solvePassiveScalar off;
16 D D [ 0 2 -1 0 0 0 0 ] 1e-9;
}
Listing 5.1: Example of a constitutiveProperties dictionary used
with rheoFoam.
The dictionary constitutiveProperties has two different sub-
dictionaries (subDict), which must necessarily exist: parameters, with information
on the constitutive model, and passiveScalarProperties, related with the scalar-
transport equation.
Regarding subDict parameters, in line 3 we define the TypeName of the consti-
tutive model to be used, which can be found in Table 4.1. In the example displayed
in Listing 5.1, we are using the Oldroyd-B model, solved with the log-conformation
approach (Oldroyd-B + Log). If we would like to use the same model without solv-
ing it with the log-conformation approach (solving the constitutive equation for
the extra-stress tensor), then the type would be simply Oldroyd-B – this naming
rule is valid for all viscoelastic models. Lines 5 to 8 specify the fluid properties
required by the constitutive model being solved. The density is a property com-
mon to all models (it is not related with the constitutive equation) and should
always be present, while the model-dependent properties can be checked in Table
4.1. Anyway, if some required parameter is not specified, the solver will retrieve
an error complaining for its absence.
The reader might be surprised with the unphysical parameters displayed in
Listing 5.1 (even just being an example), particularly the density of the fluid,
which is not realistic for any known viscoelastic liquid. The use of a unitary
density (ρ = 1 kg/m3 ) and many other unitary variables is simply to facilitate the
calculations. In the tutorials presented next, we frequently make use of this kind
of approach, since the computation of dimensionless parameters does not require
physically realistic quantities.
At line 10, the stabilization method is selected (recall that C++ is case-
sensitive): none for no stabilization; BSD to use the both-sides-diffusion technique;
coupling to use the stress-velocity coupling discussed in Section (3.3.2). Note that
CHAPTER 5. Tutorials 101
• GNF fluid: those cases only require defining pressure (divided by the den-
sity), p (in this guide represented by ρp ), and velocity, U (in this guide rep-
CHAPTER 5. Tutorials 102
resented by u) fields. For all the GNF models, except for the Newtonian
case, the solver will automatically write the shear-rate dependent viscosity
at subsequent times.
boundary conditions need to be defined to solve such equations (see the tutorial
in Section 5.1.9). On the other hand, rigid-body like motions, for example, do
not need such procedure, since the whole mesh is simply transformed by some
geometric operation.
| system/
The last steps before starting the simulation are related with the dictionaries
located in folder system/, which mainly control the numerical method. In par-
ticular, we will focus our attention on the following dictionaries: controlDict,
fvSchemes and fvSolution. All the three dictionaries must be present for
the simulation to run, as required by most of the OpenFOAM
R
solvers. Since
most of the entries in those dictionaries are transversal to both rheoFoam and any
OpenFOAM
R
solver, we will limit our description to the new features introduced
by rheoFoam.
In controlDict dictionary, the options allowing to automatically control
the time-step by imposing a Courant number limit are available in rheoFoam and
can be used (following the same principles of other OpenFOAM
R
solvers). Those
options are adjustTimeStep (on/off ), maxCo (the value of the limiting Courant
number) and maxDeltaT (the maximum admissible time-step). Furthermore, and
although not being a feature exclusive of rheoFoam, coded functionObjects can
be defined in controlDict and used with rheoFoam to extract and monitor
quantities of interest (this is not possible in foam-extend). This kind of functions
are frequently used in the tutorials of this Chapter.
Regarding dictionary fvSchemes, we remember that GaussDefCmpw schemes
(Section 4.6.1) are available for selection and can be used to discretize any con-
vective term with the generic form div(phi,variable), where variable is either
U,tau,theta or C. Still in the divSchemes subDict, the term div(grad(U))
is part of the stress-velocity coupling algorithm (see line 44 of Listing 4.2) and
should (always) be discretized using a central differencing scheme (Gauss linear ),
if used. In the gradSchemes subDict, the entry linExtrapGrad is for the gradient of
the tensor components when using linear extrapolation of polymeric extra-stress
at a given boundary, as discussed in Section 4.5.1. Apart from this, the remaining
entries in fvSchemes should be familiar to the user and the selection of appropri-
ate discretization schemes for each one is essential to keep the numerical method
accurate and stable.
The dictionary fvSolution is the only remaining to be adjusted before run-
ning the simulation. In subDict solvers, the matrix solver for each equation being
solved should be specified (remember that there will be N equations to solve for
theta/tau in a model using N modes; wildcard characters are useful in those
cases). If the user forgets to specify any, the solver will retrieve an error message
asking for it. Only the pressure equation results in a symmetric matrix of coeffi-
cients, while all the others generate non-symmetric matrices. The only exception
is the momentum equation without the convective term included, which also re-
sults in a symmetric matrix. This should be taken into account when selecting
the type of matrix solver, since some are specific for some type of matrices. In the
SIMPLE subDict, there is a new entry specific of rheoFoam, which is nInIter. This
variable was defined in Section 4.4.1 and controls the number of inner-iterations
CHAPTER 5. Tutorials 104
(see Fig. 4.2). If not defined, the solver will execute 1 inner-iteration as the de-
fault behavior. Still in the SIMPLE subDict, the entry residualControl allows the
solver to automatically stop the simulation once the residuals for all the speci-
fied variables drop below the prescribed value. This is mainly a characteristic of
steady-state solvers of OpenFOAM
R
and it can be also used in rheoFoam. How-
ever, if the goal is to run rheoFoam until the endTime specified in controlDict,
simply leave this entry empty. The last subDict in fvSolution is the relaxation-
Factors, that determines the amount of explicit (fields) and implicit (equations)
under-relaxation for each field or equation. When using the SIMPLEC algorithm
for pressure-velocity coupling, as a rule of thumb, the pressure does not need to be
explicitly under-relaxed to correct the velocity (see Ref. [2]). The only exception
occurs for non-orthogonal grids, where a small amount of under-relaxation may
eventually be needed for pressure. The under-relaxation factor, ranging between
1 (no under-relaxation) and 0 (total under-relaxation, pressure does not evolve in
time), is case-specific and should be as high as possible (it can be conditioned by
stability issues). Regarding implicit under-relaxation (equations), it only makes
sense to be used with steady-state solvers, where the absence of a time-derivative
term requires under-relaxation for stability reasons. Since rheoFoam is by default
a transient solver, where time-derivatives are present in all the transport equa-
tions, implicit under-relaxation is not needed. In practice, it is possible to run
rheoFoam without those time-derivatives by selecting a default steady-state dis-
cretization scheme in the ddtSchemes subDict of fvSchemes and, by this way,
rheoFoam will run as a typical steady-state solver of OpenFOAM
R
, requiring im-
plicit under-relaxation. However, in some situations (viscoelastic models solved
with the log-conformation approach) the user will probably face stability issues,
due to the poor diagonal dominance of the base matrix of coefficients. For this
reason, unless the user is experienced and knows what is doing, we strongly recom-
mend to use rheoFoam in transient mode. Note that, if it exists, the steady-state
will be reached after some time of a transient simulation (typically after several
relaxation times, of the order of 10 or above for higher Wi, for viscoelastic models).
If the transient analysis is not important, a high Courant number ( 1) can be
defined to reach this state faster, as long as the computations remain stable.
In all the tutorials presented next, a steady-state is reached for the range of
parameters used in the examples. However, in none of them we use the residuals
as the termination criteria. Indeed, we prefer to set a very long endTime and
monitor a relevant variable at sensitive points over time. It can be the extra-stress
near a singular point, the drag coefficient over a surface, a vortex length, or any
other relevant quantity for the problem at hand. This is usually achieved either
through a probe (when the variable exists within the solver) or a ppUtil or coded
FunctionObject (when the variable does not exist within the solver and needs to be
computed), in dictionary controlDict. The termination criteria is then based
on this variable. Of course, we have set the endTime in the tutorials based on
that analysis. The residuals displayed on the screen should not be used alone as
the termination criteria.
As mentioned in Section 4.6.2, the PostProcessing subDict enabling the use of
the ppUtil class is also defined in dictionary fvSolution. A detailed discussion
CHAPTER 5. Tutorials 105
walls
y
inlet x 2w outlet
walls
40w
! Boundary conditions
This flow is 2D, being solved in the xy-plane. A uniform velocity profile (U )
is set at the inlet, along with zero-gradient for pressure and polymeric extra-
stress components. At the outlet, fully-developed flow conditions are assumed
CHAPTER 5. Tutorials 107
(zero-gradient for all variables, except pressure, which is fixed to a constant value,
p = 0). A no-slip boundary condition is assigned at the walls (velocity is null,
polymeric extra-stress components are linearly extrapolated and a zero-gradient is
assumed for pressure in the normal direction to the wall).
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Run the solver:
∼$ rheoFoam
3–Extract profiles for u and τ along line x = 35:
∼$ sample
! Results
Figure 5.2 presents the fully-developed profiles at line x = 35. The variables
were normalized as follows: length is normalized with w, time with w /U, velocity
with U and polymeric components of the extra-stress with η0wU , as in Ref. [36].
A good agreement is observed between the numerical results and the analytical
solution written in dimensionless form [36]:
3
1 − y2
ux =
2
τxy = −3(1 − β)y
τxx = 18W i(1 − β)y 2
15 2 1.2
10 0 0.8
τxy
u
τxx
5 -2 0.4
0 -4 0.0
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
y y
(a) (b)
Figure 5.2: (a) Polymeric extra-stress components and (b) velocity, at x = 35,
for Re = 0, W i = 0.99 and β = 0.01.
The user can test the solver with a UCM fluid (β = 0) and confirm that an
accurate solution is still achieved, without facing any numerical issue. However,
running the same cases without the log-conformation approach leads to numerical
divergence (try it!).
CHAPTER 5. Tutorials 108
This tutorial probes a point in the flow over time (to check for convergence),
which has been specified in controlDict dictionary. The data is written to a
directory named probes/, whose location in the case directory depends on the
OpenFOAM
R
version.
movingLid
L fixedWalls fixedWalls
x
fixedWalls
! Boundary conditions
CHAPTER 5. Tutorials 109
The flow is assumed to be 2D, being solved in the xy-plane. For the three
stationary walls, a no-slip boundary condition is assigned, with null velocity, lin-
early extrapolated polymeric extra-stresses and zero normal gradient for pressure.
At the moving lid wall, the same boundary conditions are used for pressure and
polymeric extra-stresses. Regarding the velocity, a time-space dependent condi-
tion is employed in order to impose a smooth start of the flow, and to avoid a local
singularity with infinite acceleration at the top-right and top-left corners [9]:
Fattal & Kupferman (2005) rheoFoam Fattal & Kupferman (2005) rheoFoam
1 1.6
0.8 1.2
0.4 0.4
0.2 0
0 -0.4
-0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
u (0.5, y) x
(a) (b)
Fattal & Kupferman (2005) rheoFoam
0.02
0.016
0.012
Ek
0.008
0.004
0
0 2 4 6 8
t
(c)
Figure 5.4: (a) Velocity profile along line x = 0.5 (at t = 8), (b) Θxy profile along
line y = 0.5 (at t = 8) and (c) evolution of the average ”kinetic energy” over time.
All the results are for Re = 0.01, De = 1 and β = 0.5.
The 4:1 planar contraction is another traditional benchmark flow problem for
viscoelastic fluid flow solvers. The existence of singular points at the re-entrant
corners, where stresses grow exponentially as the corner is approached, make this
problem challenging from a numerical perspective.
This tutorial reproduces the work that we developed in Ref. [2] using an early
version of rheoFoam, where an Oldroyd-B fluid (β = 91 ) was studied for De = 0−12.
The constitutiveProperties dictionary is adjusted to reproduce the case
for De = 1 and Re = 0.01.
! Geometry & Mesh
The geometry for this case is reproduced in Fig. 5.6. The mesh corresponds to
mesh M1 of Ref. [2].
! Boundary conditions
The boundary conditions used are described in Ref. [2]. It is worth mentioning
that the time-varying inlet velocity is implemented as a codedFixedValue boundary
condition in file 0/U. The function is implemented as
( 1−cos(πt)
f ac
dirN t ≤ tlim
u(t) = (5.3)
Uav t > tlim
CHAPTER 5. Tutorials 111
(a) (b)
(c) (d)
Figure 5.5: Contours of (a) Θxx , (b) Θyy and (c) Θxy . In (d), the streamlines are
plotted. All the results are for Re = 0.01, De = 1, β = 0.5 and t = 8.
wall_vorttop
wall_liptop
y walls
inlet 8H x 2H outlet
walls
wall_lipdown
wall_vortdown
100H 100H
where Uav and dirN are vectors, and fac and tlim are scalar parameters. For
Uav = (0.25, 0, 0), dirN = (1, 0, 0), f ac = 8 and tlim = 1, this generates an inlet
velocity profile aligned with the x -axis and whose magnitude increases from 0, at
t = 0, to 0.25, at t = 1.
! Command-line
1–Build the mesh:
CHAPTER 5. Tutorials 112
∼$ blockMesh
2–Run the solver:
∼$ rheoFoam
3–Extract u and τ at the cell centers immediately upstream and downstream of
vertical line x = 0:
∼$ sample
! Results
The results obtained with mesh M1 can be found in Ref. [2].
A coded FunctionObject returns the points where the wall-parallel velocity com-
ponent changes of sign (for both the walls near to the upper lip and corner vortices).
Those points are delimiting the lip and corner vortices (if present). The user can
easily add two extra coded FunctionObjects for the vortices in the lower-half of the
contraction.
walls
inlet 2R 4R outlet
x
cylinder
walls
20R 60R
The flow is assumed to be 2D, being solved in the xy-plane. At the inlet, a
uniform velocity profile with magnitude U is imposed, the polymeric extra-stresses
are null and zero-gradient is assigned to pressure. Channel and cylinder walls are
static (velocity is null, polymeric extra-stresses are linearly extrapolated to the
walls and a zero-gradient is imposed for pressure). At the outlet, fully-developed
conditions are assumed: zero-gradient for all variables, except pressure, which is
fixed to a constant value, p = 0.
! Command-line
1–Create half of the mesh:
∼$ blockMesh
2–Reflect the half-mesh using plane xz as mirror, to obtain the full mesh:
∼$ mirrorMesh -noFunctionObjects
3–Run the solver:
∼$ rheoFoam
! Results
Figure 5.8 presents the contour plots for the first normal stress difference and
for the velocity magnitude (with superimposed streamlines), at Re = 0, W i = 0.7
and β = 0.59, using the Oldroyd-B model with the log-conformation approach.
Note that the velocity is normalized with U, time with R/U and polymeric extra-
stresses with η0RU . The drag coefficient obtained in such conditions is Cd = 117.357,
which is in reasonable agreement with Cd = 117.323 [38] and Cd = 117.315 [37].
Refining the mesh would further increase the accuracy of the numerical solution.
The drag coefficient was computed as
Z Nf
1 0
1 X 0
Cd = pI + τ · î· dS = Sf · pf I + τf · î
η0 U h η0 U h k=1
S
CHAPTER 5. Tutorials 114
where Sf is a vector normal to each face of the cylinder boundary, whose magnitude
is equal to the face’s area, h is the depth of the cylinder in the neutral (empty)
direction and î is a unitary vector aligned with the streamwise direction. The drag
coefficient is retrieved on run time by a coded FunctionObject, which can be found
in controlDict dictionary.
(a)
(b)
Figure 5.8: (a) First-normal stress difference and (b) velocity magnitude contours
with superimposed streamlines, at t = 15, for W i = 0.7, Re = 0 and β = 0.59.
It should be noted that, although the mesh displays some amount of non-
orthogonality (run checkMesh to confirm it), the non-orthogonality corrector loop
is not used (this is to say that only one iteration is performed) and the explicit
CHAPTER 5. Tutorials 115
pressure under-relaxation factor is kept relatively high (0.9). These features show
the enhanced stability of the solver. However, to study for example the time
evolution of the drag coefficient, the use of the non-orthogonality corrector loop,
combined with inner iterations, is suggested in order to increase the time accuracy
(by decreasing the explicitness of the method). For a high-enough number of
iterations in the non-orthogonal correction loop, under-relaxation of pressure may
eventually be avoided.
inlet_north
walls walls
outlet_west W outlet_east
x
9.5W
inlet_south
Figure 5.9: Cross-slot geometry composed of 4 arms with two balanced inlets
and two outlets.
! Boundary conditions
The flow is assumed to be 2D, being solved in the xy-plane. At both inlets, a
uniform velocity profile is specified, with magnitude U and pointing to the origin,
so that those two streams are flowing in opposite directions. The polymeric extra-
stresses are null and for pressure a zero-gradient is used. The walls are stationary,
thus velocity is null, polymeric extra-stresses are linearly extrapolated and the
pressure is assumed to not change in the normal direction. Fully-developed flow
conditions are assigned at the outlets: null normal gradient for all variables, except
pressure, which is fixed at p = 0.
A passive scalar (tracer field C) is added to the problem, which requires the
assignment of initial and boundary conditions. We impose a continuous injection
of C at inlet north (C = 1), while no tracer is injected at inlet south (C = 0).
In the remaining boundaries, we impose null normal gradient (meaning no flux of
CHAPTER 5. Tutorials 117
C across the walls and fully-developed flow conditions at both outlets). At time
t = 0, when the simulation is started, the y-positive portion of the cross-slot is
filled with the tracer (C = 1, y > 0 ∧ t = 0).
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Create field C by copying one already present, which is not initialized in the
interior domain:
∼$ cp 0/C.org 0/C
3–Initialize field C in the interior domain (C = 1, y > 0 ∧ t = 0):
∼$ setFields
4–Run the solver:
∼$ rheoFoam
! Results
The contours for some important variables are displayed in Fig. 5.10, for Re =
0, W i = 0.33, β = 0 (UCM model) and P e = 500. The variables are normalized
with U (velocity), W/U (time) and ηW 0U
(stresses).
The local Weissenberg number at the origin, defined in Ref. [39] as the product
of the relaxation time by the velocity gradient magnitude at the stagnation point
streamlines, is W i0 = 0.523, which is close to the benchmark value obtained in
a similar mesh, W i0 = 0.509 (Ref. [39], for mesh M1). The local Weissenberg
number at the origin is retrieved by the solver to the case directory, through a
coded FunctionObject that can be found in dictionary controlDict.
The importance of stress-velocity coupling in this case can be evaluated by re-
running the tutorial with either BSD only or no stabilization method (stabilization
= none or BSD, in dictionary constitutiveProperties).
Checkerboard fields easily develop in such conditions (this is a critical case due to
the use of a UCM fluid).
Note that this tutorial makes use of a variable time-step, controlled by a maxi-
mum Courant number fixed at 0.4. This strategy is used because only steady-state
results are of interest. This is also the reason to use a high tolerance in the sparse
matrix solvers, in fvSolution dictionary – the steady asymmetry in the flow
develops faster with these conditions, since the transient numerical error is higher.
(a) (b)
(c) (d)
Figure 5.10: (a) Velocity magnitude contours with superimposed streamlines, (b)
contours of C, (c) first-normal stress difference and (d) τxy contours, for Re = 0,
W i = 0.33, β = 0 and P e = 500.
The tutorial presented in this Section addresses the simulation of blood flow in
a real-model aneurysm. Contrarily to the previous tutorials, this case is based on
a 3D polyhedral (non-orthogonal) mesh and uses a GNF model. Furthermore, the
flow is simulated for a moderate Reynolds number, which will test the robustness
of the solver for such conditions, where inertia already plays an important role.
The Herschel-Bulkley model was selected to simulate the blood rheology, fol-
lowing Ref. [40]. The generalized Reynolds number for this model, assuming τ0 = 0
– the power-law limit –, is [40]:
2−n
ρ(2Rin1 )n U
ReGN = n
k 3n+1
4n
8n−1
This tutorial simulates the flow at ReGN = 420.
! Geometry & Mesh
The STL file of the aneurysm surface was downloaded from a repository with
real-model aneurysms [41], extracted from 3D rotational angiographies of diseased
patients. From the list of available models, case ID C0005 was selected, which
refers to an aneurysm in the internal carotid artery (ICA).
The surface is composed of one main entry vessel (in1, the ICA), which bi-
furcates into two smaller vessels (out1 and out2 ), Fig. 5.11a. The aneurysm is
located near the bifurcation point. Due to the long extension of vessels upstream
and downstream of the aneurysm in the original STL file, all the vessels were
shortened and a cylindrical extension was connected to each one, in order to min-
imize entry/exit effects in the region near the aneurysm. The locations where
CHAPTER 5. Tutorials 119
those connections were established are highlighted in Fig. 5.11a using red arrows.
The transition length between the tube connectors and the vessels is typically 10
% of the vessel radius and the radius of those tubes is equal to the equivalent
radius of the vessels at the connection point. The radius of each inlet/outlet is:
Rin1 = 1.66 mm, Rout1 = 1.17 mm and Rout2 = 0.95 mm.
in1 walls
out2
out1
(a)
(b)
Figure 5.11: (a) Geometry of the aneurysm considered in the tutorial. Red
arrows point to the transition regions between the aneurysm and the cylindrical
extensions. The radius of each inlet/outlet is: Rin1 = 1.66 mm, Rout1 = 1.17 mm
and Rout2 = 0.95 mm. (b) Detailed view of the mesh on patch in1, zooming the
cell layers near the wall. The reference axis for the geometry is centered on patch
in1, with the normal vector of the patch pointing in the negative direction of the
y-axis.
The mesh was built using 1 cfMesh, a meshing tool available in foam-extend
since version 3.2. The maximum cell size was limited to 5 mm and boundary
cell layers were generated near to the vessel walls in order to accurately solve the
gradients developed there, Fig. 5.11b. The mesh provided with this tutorial has
1
http://cfmesh.com/
CHAPTER 5. Tutorials 120
3n + 1
r n+1
n
U =U 1− (5.4)
n+1 R
where U is the mean velocity. In our case, U is constant over time, thus steady
conditions are simulated, instead of the cardiac cycle (this is to shorten the simu-
lation time, since transient solutions would require a lower time-step, leading to a
higher computational time). Regarding the walls, a no-slip boundary condition is
imposed. At the outlets, the pressure is fixed to zero and the velocity is assumed
to be fully-developed (zero gradient for velocity). Since the Reynolds number used
in the tutorial is well below the critical value for transition to the turbulent regime,
no special conditions need to be defined regarding turbulence modeling.
! Command-line
The mesh is already built and can be found in folder polyMesh/.
1–Decompose the case among 2 processors to speed-up the computations (with
2 processors, it takes around 1h to reach convergence in a laptop with an Intel
i5-3210M processor, 2.5 GHz):
∼$ decomposePar
2–Run the solver in parallel, using 2 processors:
∼$ mpirun -np 2 rheoFoam -parallel
3–Reconstruct the last time-step of the case for post-processing:
∼$ reconstructPar -latestTime
! Results
The results obtained at ReGN = 420 are displayed in Fig. 5.12, where both the
streamlines and the wall shear-stress magnitude (WSSmag) contours are shown.
The wall shear-stress magnitude is computed through a ppUtil (described in Sec-
tion 4.6.2), whose settings can be found in dictionary fvSolution, under subDict
PostProcessing.
This tutorial is defined to run with an adjustable time-step, controlled by a
maximum Courant number fixed at 50. The high Courant number is to quickly
CHAPTER 5. Tutorials 121
achieve the steady-state, without the need to cancel the time-derivatives. The
non-orthogonality corrector loop was turned on, with 1 iteration per time-step
(the pressure is also under-relaxed with a factor of 0.9), in order to avoid possible
numerical issues due to non-orthogonality. The convergence can be monitored by
a probe located at one of the exit vessels, downstream of the aneurysm.
(a)
(b)
Figure 5.12: (a) Streamlines (colored with the velocity magnitude) and (b) wall
shear-stress magnitude contours, at steady-state and for ReGN = 420.
Viscous fluid dampers react to applied loads with the generation of a back
pressure due to the high resistance experienced by a viscous fluid flowing through
narrow orifices. These mechanical elements are typically composed of a moving
piston enclosed inside a cage filled with a viscous fluid, with a narrow gap left
between the piston and the cage wall. The motion of the piston and the shaft
to which it is connected drive the fluid flow inside the cage. This tutorial aims
to simulate the fluid dynamics inside a damper subjected to an oscillatory load,
whereby the piston and shaft motion are represented by a moving mesh.
The setup adopted in this tutorial aims to closely reproduce the conditions
presented by Syrakos et al. [42]. The case represented in this tutorial is for the
Careau-Yasuda fluid, named CY-100 in [42], at 32 Hz. Other conditions can be
easily tested by simple adjustment of the input parameters.
! Geometry & Mesh
The fluid damper geometry is depicted in Fig. 5.13. The geometry is equal
to that presented in [42], and the piston edges are rounded by the same radius
of curvature (Rcurv = 0.75 mm) as therein. The damper is axisymmetric, thus
only a slice of the total domain is simulated (Fig. 5.13). A Cartesian system of
coordinates is shown in Fig. 5.13 because that is the one used by OpenFOAM
R
to
solve the problem, although keep in mind that symmetry around axis Ox holds.
The mesh is built using blockMesh and has a spatial resolution in the gap region
between patches top and piston, which is approximately half of mesh M1 in [42].
45 mm 20 mm 45 mm
top
top
3 mm
y y shaft
shaft
5 mm
x z z x
Figure 5.13: Fluid damper geometry. The geometry is symmetric around axis Ox,
thus the computational domain corresponds to only a slice of the whole geometry.
! Boundary conditions
The Navier slip boundary condition is applied to the velocity field in all the
boundaries (except on wedges), in agreement with [42]. In addition, the velocity of
patches shaft and piston receive an extra contribution due to their own oscillatory
motion,
where α = (α, 0, 0) is the amplitude vector, α = 0.012 m and ω = 64π rad/s. Note
that Syrakos et al. [42] used a sine function for the motion, whereas a cosine is used
here for convenience (the two functions are offset by π/2 rad). A zero-gradient
condition is imposed for pressure.
Although physically both the shaft and piston move, in practice we only let
the piston patch to follow exactly the motion described by Eq. 5.5, while the
faces of patch shaft stretch or shrink according to the motion imposed by patch
piston (slip boundary condition). The same happens for patch top, while the
remaining patches remain static. The mesh motion is obtained from the solution
of a Laplace equation with the previously described boundary conditions and a
variable diffusivity coefficient (see dynamicMeshDict dictionary). Note that
the boundary conditions presented in this paragraph are not related with the fluid
velocity field, whose boundary conditions were reported in the previous paragraph.
Still, such coupling can be obtained with boundary condition movingWallVelocity,
whereby the given patch velocity is used as boundary condition to the fluid velocity
(this can be used, for example, for patch piston, but not for patch shaft).
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Run the solver:
∼$ rheoFoam
3–Check file Ffl generated by the solver.
! Results
The oscillatory piston motion can be directly observed in Paraview. In addition,
the solver writes in runtime the reaction force exerted by the fluid on the piston
and shaft (file Ffl),
Z Nf
0 0
X
Ffl = −pI + τ · î· dS = Sf · −pf I + τf · î
S k=1
2000
1000
Ffl (N)
-1000
-2000
-3000
-15 -10 -5 0 5 10 15
Displacement (mm)
Figure 5.14: Reaction force exerted by the Carreau-Yasuda fluid on the piston
and shaft over time, for f = 32 Hz. The results of Syrakos et al. [42] are also
plotted for reference.
5.2 rheoTestFoam
5.2.1 General guidelines
In Section 4.4.2, rheoTestFoam was presented as a testing application for the con-
stitutive models implemented in rheoTool , being not a general-purpose solver. For
this reason, some of the steps usually required to setup a generic simulation in
OpenFOAM
R
are not necessary with rheoTestFoam, while, on the other hand,
extra-inputs need to be specified.
| constant/
One main difference of rheoTestFoam cases regarding, for example, rheoFoam
cases is in the mesh: the user should always use the same single-cell unitary mesh
when working with rheoTestFoam. Thus, changing the mesh from case to case is
unnecessary and not recommended.
The dictionary constitutiveProperties is composed of two subDict: pa-
rameters and rheoTestFoamParameters, as displayed in Listing 5.2. The entries in
parameters have exactly the same meaning as previously discussed for rheoFoam.
However, there are two entries which remain inactive in rheoTestFoam: rho and
stabilization. Remember from Section 4.4.2 that rheoTestFoam is only solving the
constitutive equations for a given ∇u tensor, thus those two parameters related
with the momentum equation are useless. Nevertheless, they should be present
(with any assigned value) to avoid a run time error. rheoTestFoamParameters is
a subDict specific of rheoTestFoam, in the same way as passiveScalarProperties is
CHAPTER 5. Tutorials 125
a particular subDict of rheoFoam. The keyword ramp stands for the operation
mode (see Section 4.4.2): ramp (true) or transient (false). The other two entries
define tensor ∇u, since we consider ∇u = gammaEpsilonDotL[i]· gradU, where
i is the index representing each entry of list gammaEpsilonDotL. If ramp = f alse,
the mode is transient and only one entry is expected in gammaEpsilonDotL – the
solver is testing the transient behavior of the constitutive model, for a (single)
given ∇u. On the other hand, in ramp mode (ramp = true), gammaEpsilonDotL
may have as many entries as defined by the user and steady-state variables will be
returned by the solver for each entry. Any combination of components is admissi-
ble for gradU, although only some correspond to canonical rheometric flows. The
one displayed in Listing 5.2 is for a pure-shear flow: u = (γ̇y, 0, 0), where γ̇ is the
gammaEpsilonDotL value.
1 parameters
{
3 type Oldroyd-B;
When running rheoTestFoam in ramp mode, the user does not have control
on deltaT (time-step), nor on the endTime. The time step is automatically set
based on the relaxation time and strain-rate values for viscoelastic fluids or is sim-
ply set to 1 s for GNF models (in this case, the value is not important since no
equation is solved implicitly). The endTime in ramp mode is not important, since
the stopping criteria is based on an hard-coded threshold for the residuals and for
the number of iterations. On the other hand, in transient mode, both variables
should be specified by the user in controlDict. Regarding the discretization
schemes (fvSchemes dictionary), only time-derivatives and grad(U) are used by
the solver. The discretization of grad(U) should be kept as Gauss linear, while
any valid time-scheme can be selected (except steady-state), although in ramp
mode this should not make any difference, since we are looking for steady-state
solutions. In dictionary fvSolution, the matrix solvers required by the consti-
tutive equations must be defined and the number of inner-iterations may also be
controlled if running in transient mode. Note that since the mesh has only one
cell, a good time accuracy can be achieved by selecting a small time-step, without
compromising the CPU time (in general, simulations will be always fast). The use
of under-relaxation is not needed, as long as time-derivatives are not disabled in
fvSchemes (this is our recommendation).
For the range of shear-rates between 0.01 s-1 and 10000 s-1 , the Herschel-Bulkley
model behavior is displayed in Fig. 5.15. The model predicts a shear-thinning
behavior for γ̇ > γ̇0 , where γ̇0 is the critical strain-rate at which η = η0 . For the
parameters defined in this example, γ̇0 = 0.13 s−1 .
1 100
10
0.1
1
τ'xy (Pa)
η (Pa.s)
0.1
0.01
0.01
0.001 0.001
0.01 0.1 1 10 100 1000 10000
-1
Shear-rate (s )
0
Figure 5.15: Shear viscosity and τxy (the only non-zero component of the sym-
metric extra-stress tensor) as a function of the shear-rate, in a steady shear
flow, for the Herschel-Bulkley model with parameters: τ0 = 0.0175 Pa, k =
8.9721 × 10−3 Pa.sn , n = 0.8601 and η0 = 0.15 Pa.s.
Pa.s, λ = 1 s and different values of L2 were tested (10, 100 and 1000). In such
0 0
τ −τ
conditions, the extensional viscosity, defined as ηE = xx ε̇ yy , is given by [30]
2 1
ηE = 3ηs + ηp + (5.6)
1 − 2λε̇/f 1 + λε̇/f
where f is the solution of the cubic equation
L2 = 10 L2 = 100 L2 = 1000
10000
1000
ηE/η0
100
10
1
0.01 0.1 1 10 100
λε ̇
(a)
Wi = 2 Wi = 5 Wi = 10
1000
100
ηE+/η0
10
0.1
0 1 2 3 4 5
t/λ
(b)
Figure 5.16: (a) Steady extensional viscosity (ηE ) as a function of W i = λε̇, for
different values of L2 (points represent the numerical results of rheoTestFoam and
the lines correspond to the analytical solution of Eq. 5.6); (b) transient extensional
viscosity ηE+ for different Wi, at fixed L2 = 100. The remaining parameters of the
FENE-CR model are ηs = 0.1 Pa.s and ηp = 0.9 Pa.s.
CHAPTER 5. Tutorials 130
5.3 rheoInterFoam
| system/
The main novelty comparing to rheoFoam is that when using rheoInterFoam
the user has the possibility to choose between PIMPLE and SIMPLEC for pressure-
velocity coupling. This is controled in dictionary fvSolution, in subDict PIM-
PLE, where the keyword SIMPLEC can be assigned to true or false. Independently
of the choice, the momentum equation is always solved. The variable nCorrectors
works in its usual way (looping the pressure equation) and the variable nInIter
assumes the same function as nOuterCorrectors. In a future release of rheoIn-
terFoam, this workflow will most likely change. There are currently these two
options because it is still not clear which one is more advantageous. Still in dictio-
nary fvSolution, other keywords must be assigned in subDict PIMPLE, which
are related with the VOF method and that the reader can found in the tutorials.
In dictionary fvSchemes, the discretization schemes for the two phases should
be defined, as well as the discretization schemes related with the VOF method,
which can also be found in the tutorials.
Besides the Courant number, the solvers using VOF, as rheoInterFoam, can
also restrict the time-step based on an interface Courant number, which should be
defined in dictionary controlDict.
The dictionary setFields, used to initialize parts of the domain with speci-
fied values, should also be present in folder system/ whenever used.
The domain is meshed uniformly with 120 cells in both the radial and axial
directions.
atmosphere
phase 2
phase 1
D
g
3D
axis U0 atmosphere
2D
plate
3D
! Boundary conditions
At the plate, no-slip boundary conditions are imposed with a null velocity,
linearly extrapolated polymeric extra-stress components and zero normal gradient
for the indicator field. We assign a fixedFluxPressure BC to the pressure, as
discussed in Section 4.5.8, for multiphase flows. The patches representing the
open boundaries (atmosphere) are assumed to not interfere with the dynamics of
the drop, so that zero-gradient is assumed for all variables, except the pressure,
which is fixed at p = 0.
Following Ref. [43], the drop has an initial velocity U0 , in the vertical direction
(pointing downwards), and its center of mass is at a distance 2D from the plate.
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Initialize the indicator and velocity fields in the drop region:
∼$ cp 0/alpha1.org 0/alpha1
∼$ cp 0/U.org 0/U
CHAPTER 5. Tutorials 133
∼$ setFields
3–Run the solver:
∼$ rheoInterFoam
! Results
The evolution of the drop width over time is plotted in Fig. 5.18. The width
is normalized with D (its initial diameter) and time with D/U0 . The evolution of
the drop width over time is written to a file by a coded FunctionObject.
1.6
1.4
W
1.2
0.8
0 1 2 3 4 5
t
Figure 5.18: Evolution of the drop width over time for F r = 2.26, Re = 5,
W i = 1 and β = 0.1. The profile obtained with rheoInterFoam is compared with
the data of Figueiredo et al. [43].
atmosphere
atmosphere
phase 2
4H
wallIn wallOut
y
inlet H phase 1 symmetry
outlet
x
! Boundary conditions
The boundary conditions used are described in Ref. [45]. In order to avoid the
definition of an analytical solution at the inlet, a long entrance channel is used.
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Initialize the color function field:
∼$ cp 0/alpha1.org 0/alpha1
∼$ setFields
3–Run the solver:
∼$ rheoInterFoam
! Results
CHAPTER 5. Tutorials 135
The results obtained with mesh M1 can be found in Ref. [45]. The users of
the foam-extend version (fe40 ) will notice an instability in the extrudate surface
close to the die exit, which eventually vanishes with time. Such instability is not
present in the OpenFOAM
R
versions.
Several methods can be used to compute the amount of swell far from the die
exit. Assuming that the free-surface of the extrudate corresponds to α = 0.5, for
α ∈ [0, 1], this isoline can be directly extracted with Paraview. Another option
is to extract the profile of the color function (α) over a vertical line far from
the die exit, using the sample utility of OpenFOAM
R
. Then, the value α = 0.5
can be simply interpolated from the nearest values. A more complex, still more
versatile option is to code a ppUtil function (Section 4.6.2) retrieving in run time
the maximum free-surface position in a predefined region. This would also allow to
check for convergence of the free-surface position. Slight differences can be found
among the results from the different methods.
5.4 rheoEFoam
5.4.1 General guidelines
The first steps to set up a case for rheoEFoam are the same as for rheoFoam
(Section 5.1.1). After that, the hydrodynamic component of the problem will be
ready, and only the electric component will remain to be defined – this is the
subject of this section.
| constant/
The electricProperties dictionary should be added to folder constant
/. It contains most of the information about the EDF model to be used, as shown
in Listing 5.3, illustrating an example for the PNP model.
parameters
2 {
type NernstPlanck;
4
T T [ 0 0 0 1 0 0 0 ] 298;
6 relPerm relPerm [ 0 0 0 0 0 0 0 ] 80.1;
8 psiContrib true;
extraEField extraEField [ 1 1 -3 0 0 -1 0 ] (5000 0
0);
10
species
12 (
cCation
14 {
z z [ 0 0 0 0 0 0 0 ] 1;
16 D D [ 0 2 -1 0 0 0 0 ] 1e-9;
}
18
cAnion
20 {
CHAPTER 5. Tutorials 136
z z [ 0 0 0 0 0 0 0 ] -1;
22 D D [ 0 2 -1 0 0 0 0 ] 1e-9;
}
24 );
}
Listing 5.3: Example of a electricProperties dictionary used with
rheoEFoam – the settings displayed are for the PNP model.
The electric properties in dictionary electricProperties are defined in-
side a subDict named parameters (line 1, Listing 5.3). The EDF model is selected
through keyword type, where the TypeName of any model in Table 4.2 can be used
(the user may type any random word to get the list of all available EDF models).
Apart from the type keyword, all the remaining entries are model-specific. In the
case of the PNP model in Listing 5.3, T is the absolute temperature and relPerm
corresponds to the relative permittivity (dielectric constant) of the electrolyte, εR ,
such that ε = εR ε0 , where ε0 is the vacuum permittivity. In line 8, the entry
psiContrib set to true indicates that variable psi (either Ψ or ψ) should con-
tribute to the electric field used in the definition of the electric body-force. The
default behavior, i.e., if the entry is not defined, is psiContrib = true. The next
line (line 9) is also optional and allows the user to define an additional uniform
electric field, extraEField (units are in SI, as always). This electric field will only
enter in the computation of the electric body-force in the momentum equation
(fE = ρE (E + extraEField ); extraEField corresponds to vector Ea in Table
4.2). Then, from line 11 to 24 each specie of the electrolyte is defined. In the
example, only two species are modeled: cCation and cAnion, each having its own
charge valence (z ) and diffusivity (D). Note that the charge valence is a signed
integer: positive for cations and negative for anions. The user can add as many
species as desired to the list, for the model under analysis. The name given to
each specie is user-defined, but consistence must be kept when further defining the
respective fields, as explained next.
The example in Listing 5.3 should not be generalized to all the EDF models, as
stated before. The best way for the user to know how the electricProperties
dictionary should look like for a given EDF model is to analyze a tutorial provided
for that model (at least one tutorial is provided for each model).
| 0/
In addition to the fields related with the hydrodynamics (pressure, velocity and
eventually extra-stress), when using rheoEFoam for EDFs, the fields specific to the
given EDF model should be specified in folder 0/ (or the equivalent starting time
folder when different from 0).
For the PNP model illustrated in Listing 5.3, we need to define 3 or 4 fields,
depending if we use one single electric potential or two, respectively. In the first
case, the fields would be psi (in this guide represented by Ψ ), cCation and cAni
on (in this guide represented by ci ). In the second case, field psi would be replaced
by psi (in this guide represented by ψ) and phiE (in this guide represented by
φExt ). Note that although having the same name in both cases, field psi has
different meanings for each one: it is either the total, unique electric potential (Ψ ),
CHAPTER 5. Tutorials 137
or the intrinsic electric potential (ψ) – in practice, the Poisson equation to be solved
is different in each case. The selection between both cases is made through variable
phiE: when present, psi is considered the intrinsic electric potential, ψ. It is
important to highlight that the fields representing the concentration of each specie
(in mol/m3 ) should keep the name defined in dictionary electricProperties,
in a one-to-one correspondence for each specie. These names are user-defined, in
opposition to the names for the electric potential variables, which are fixed: psi
and phiE.
The PNP model is the only to require the definition of fields for the concentra-
tion of each ionic specie. The other EDF models only require the electric potential
(one or two variables, depending on the model and on the user’s choice) to be
defined. The exception is the Ohmic model, which also requires a field for the
conductivity (sigma, in this guide represented by σ). Still for this model, only
one electric potential variable may exist and its name should be phiE. In case of
doubt, checking the tutorials is always a good starting point.
| system/
Since additional equations are solved for EDFs (comparing to pressure-driven
flows), the discretization schemes for the terms entering these equations need to
be defined, as well as the sparse matrix solvers and respective settings.
The discretization schemes are defined in dictionary fvSchemes. The new
entries to add to the dictionary are model-dependent and can be found in the
tutorials. If the discretization scheme for any term is missing, an error will be
retrieved complaining for it.
Regarding dictionary fvSolution, keep in mind when defining the matrix
solvers for the new fields that Poisson-type equations require, in general, a sym-
metric matrix solver, while generic transport equations (including advection) are
usually handled with an asymmetric matrix solver. Regarding under-relaxation,
our recommendations are the same as the ones expressed for rheoFoam: by default,
do not use under-relaxation, except, eventually for pressure in non-orthogonal
grids, if needed. Note that the Nernst-Planck equations for each specie in the
PNP model are collectively solved under the name ci, instead of the name given
to the specie (this should be taken into account when defining the matrix solver
and the under-relaxation factors).
Still in dictionary fvSolution, a new subDict needs to be defined for EDFs,
named electricControls, Listing 5.4. In the code sample analyzed in Section 4.2.4,
we have seen that each equation of a given EDF model is solved inside a while loop,
controlled by a maximum allowable number of iterations and the initial residual
of the equation being solved (the loop is exited when the first of the two criteria
is met). These loops are intended to converge explicit terms inside each equation,
since this can be critical for some EDFs. Thus, the controlling parameters of
these cycles – the maximum number of iterations and the threshold residual – are
defined in subDict electricControls. This needs to be done for each equation, or
default values are assumed otherwise (maxIter : 50; residuals: 10−7 ). For non-
stiff problems and when the mesh non-orthogonality is kept low, 1 iteration can
be enough. When the PNP model is used, the number of electrokinetic coupling
iterations (see Section 4.2.3) can also be defined in this subDict, as shown in line 3
CHAPTER 5. Tutorials 138
5 phiEEqn
{
7 residuals 1e-7;
maxIter 1;
9 }
11 psiEqn
{
13 residuals 1e-7;
maxIter 1;
15 }
17 ciEqn
{
19 residuals 1e-7;
maxIter 1;
21 }
}
Listing 5.4: Example of an electricControls subDict in dictionary
fvSolution – the settings displayed are for the PNP model.
I tutorials/rheoEFoam/channelEDF/PowerLaw/PoissonBoltzmann
! Overview
The analytical solution for the EDF of a power-law fluid in a slit microchannel
can be found in Ref. [46] for a generic flow behavior index (n) of the power-law
model, under the Debye-Hückel approximation. For fully-developedq flow condi-
2
2(zH) ec0 F
tions, the velocity profiles depend on n and on κ̃ = κH = λHD = εkT
,
where λD is the Debye length for a binary, symmetric electrolyte, as defined in
Eq. (3.34), and H is the channel half-width. Although we will only present results
for the PB model, we also provide the corresponding cases for the PNP and DH
CHAPTER 5. Tutorials 139
models – as you will see, the results are indistinguishable between the models, for
the conditions simulated.
! Geometry & Mesh
The geometry is a 2D (slit) microchannel with half-width H, Fig. 5.20. Due
to the periodic boundary conditions assumed on the inlet/outlet, the mesh has
one single cell in the x -direction and 300 non-uniformly distributed cells in the
y-direction, normal to the applied electric field.
walls
y E H
inlet outlet
x
symmetry
! Boundary conditions
The flow is 2D in the xy-plane. At both the inlet and outlet, periodicity is
imposed. Thus, we assume from the beginning that this condition will retrieve the
fully developed profiles for an infinitely long microchannel. A symmetry condition
is imposed at y = 0. The walls are impermeable (u = 0 and zero-gradient for
pressure) and have a fixed intrinsic electric potential. Due to the simple geometry
of the channel, the flow is generated by imposing a uniform electric field through-
out the channel, parallel to the walls (E). Therefore, there is no need to solve
for the external electric potential variable (phiE ): the electric field is directly im-
posed through the extraEfield entry of dictionary electricProperties. These
boundary conditions hold for both the PB and DH models. In order to use the PNP
model, a boundary condition must be defined at the wall for the ionic species. Since
Boltzmann equilibrium holds and only the steady-solution is sought, the easiest
way is to simply compute the ionic concentration from the potential distribution
(check the boltzmannEquilibrium boundary condition described in Section 4.5.4).
The flow is initially at rest and the intrinsic electric potential is zero.
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Run the solver:
∼$ rheoEFoam
3–Extract the profile of u along the vertical direction:
∼$ sample -latestTime
CHAPTER 5. Tutorials 140
! Results
The velocity profiles over the y-direction are depicted in Fig. 5.21 for varying n
= 0.25, 0.5, 0.75, 1 and 1.5, at fixed κ̃ = 15. The velocity profiles are normalized
by the velocity at the centerline of the channel, while the spatial coordinate is
normalized by the channel half-width. The numerical results reproduce accurately
the analytical solution [46] and show that shear-thinning fluids (n < 1) display an
apparently compressed EDL, while the opposite is observed for shear-thickening
fluids (n > 1).
0.8
0.6
|u|
0.4
0.2
0.0
0.5 0.6 0.7 0.8 0.9 1
y
Figure 5.21: Velocity magnitude over the direction transverse to the applied
electric field for a power-law fluid with different flow behavior index, at fixed
κ̃ = 15. The points represent numerical values, while the lines are the analytical
solution [46]. Note that the x -scale has been truncated and only represents one-
quarter of the channel width (this is to have a zoomed view near the wall).
As aforementioned, the same results would be obtained with the PNP and DH
models and the user may confirm that from the cases provided. All the cases
are prepared for n = 0.75 (this can be easily changed in dictionary constitu
tiveProperties) and will converge in the total time of simulation defined in
controlDict. However, we should note that for low n, where n = 0.25 can be
included, the simulation will take a much longer time to reach the steady-state,
because of the low viscosity that develops near the wall, where the shear-rate
attains very high values (it becomes even worse for higher κ̃). Furthermore, care
should be taken in defining the upper and lower bounds for the viscosity (check
the power-law model implementation in Table 4.1), since stringent bounds may
influence the numerical solution.
I tutorials/rheoEFoam/channelEDF/PTTlinear/DebyeHuckel
! Overview
Afonso et al. [47] derived an analytical expression for the mixed pressure-
/electrically-driven flow of simplified (ζ = 0) linear PTT fluids in slit channels
with an homogeneous zeta-potential at the walls, under the Debye-Hückel approx-
imation. For creeping flow conditions (Re = 0), the velocity profiles depend on
H 2 |∇p|
Γ = − ψ 0 |E|
(ratio between pressure and electric forcing), κ̃ = κH = λHD =
√ √
q
2(zH)2 ec0 F
εkT
and εDe = ελκU , where ε is the extensibility parameter of the
PTT model, U = − ψη00|E| is the Helmholtz–Smoluchowski velocity and ψ0 is the
zeta-potential at the wall. Both ∇p and E only have a single non-zero component,
in our case, the x -component. Note that in this tutorial we use exceptionally to
represent the electric permittivity (instead of ε), in order to distinguish it from ε
that we use to represent the extensibility parameter of the PTT model.
In this tutorial, we analyze the effect of varying Γ , while keeping the remaining
dimensionless parameters fixed.
! Geometry & Mesh
The geometry is similar to the one used in Part A for the power-law fluid ex-
ample. However, the physical dimensions are different and, importantly, we do not
consider periodicity between the inlet and outlet. Instead, the flow is also solved
in the x -direction. We note that this is not mandatory and that cyclic conditions
would also provide the right solution. However, considering the full channel in the
x -direction eases the definition of the pressure gradient, at the expanse of having a
higher number of cells in the mesh. Several options would allow to keep the cyclic
patches and to impose directly the pressure gradient in the momentum equation,
as the fvOptions tool, but we consider them less straightforward than our choice
(at least for less experienced users).
! Boundary conditions
The boundary conditions are similar to the ones used in Part A, with some
modifications required by the use of a viscoelastic model and due to the different
conditions assigned to patches inlet and outlet. Since pressure gradients are al-
lowed, we fix the pressure at the outlet and adjust the inlet pressure as required
to get the desired Γ . Remember that in rheoEFoam field p represents the pressure
divided by the density. A zero-gradient condition is imposed for the remaining
variables at those two patches. Regarding the wall, the polymeric extra-stresses
are linearly extrapolated.
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Run the solver:
∼$ rheoEFoam
CHAPTER 5. Tutorials 142
3–Extract the profile of u along the vertical direction for the latest time (the
x -position of the sampling line should not be important):
∼$ sample -latestTime
! Results
Fig. 5.22 shows the velocity profiles under different forcing ratios. The velocity
is normalized by the Helmholtz–Smoluchowski velocity (U, defined above), while
the spatial coordinate is normalized by the channel half-width. Note that for a
Newtonian case and Γ = 0 (pure EDF), the (normalized) velocity profile at the
centerline would have a value very close to 1 (the higher κ̃, the closer it is) and
we can see that this value is significantly higher for a linear PTT fluid due to
shear-thinning.
The parameters provided in the tutorial are for Γ = 4.
25
20
|u|
15
10
0
0 0.2 0.4 0.6 0.8 1
y
Figure 5.22: Velocity magnitude along the direction transverse to the applied
electric field and pressure gradient
√ for a simplified linear PTT fluid, at different
forcing ratios Γ , for κ̃ = 20 and εDe = 4. The points represent numerical values,
while the lines are the analytical solution [47].
elecNorth
50R θ
r
y
wallWest wallEast
2R x
cylinder
50R
elecSouth
Figure 5.23: Metallic cylinder placed over an electric field. The surrounding
domain is square (edge size: 100R) and the cylinder lays on its center. The cylin-
drical coordinate system (r, θ) is plotted in black, while the Cartesian coordinate
system (x, y) is represented in red (remember that OpenFOAM
R
uses the Carte-
sian system for computations).
! Boundary conditions
The boundary conditions are described in detail in Ref. [3]. It is worth to
note that the boundary conditions of this tutorial, for the PNP model, are dif-
ferent between the versions of rheoTool . Indeed, we set a zeroGradient BC for
the pressure at the cylinder surface in version fe40. However, in version of40 we
use the more accurate fixedFluxExtrapolatedPressure BC, which is not available in
the other versions (see Section 4.5.8). In order to keep the code of rheoEFoam
unchanged, this BC applied in the specific context of this problem – the pressure
field needs an internal reference value/point, because no Dirichlet BC is defined for
CHAPTER 5. Tutorials 144
it – requires that at least one of the remaining boundaries does not fix the velocity
(u = 0 is used in fe40 ). A simple workaround to this situation, which would
allow to keep the no-slip BC in all the boundaries (without violating continuity),
would be to comment the line where function adjustPhi() is called in rheoEFoam
(file pEqn.H). In practice, both the methods retrieve similar solutions, since the
bounding domain is significantly larger than the cylinder radius. In addition, the
major difference between using zeroGradient (fe40 ) or fixedFluxExtrapolatedPres-
sure (of40 ) for the pressure on the cylinder surface is essentially observed in the
velocity values in the first 2–3 cells adjacent to the cylinder (more evident at
θ = 0 and 180◦ ). Farther away from the surface, the effect is barely detected.
The differences are more noticeable at higher voltages, but they vanish with mesh
refinement, since this compensates the low-order accuracy of the zeroGradient BC.
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Run the solver:
∼$ rheoEFoam
3–Extract profiles of u along the line θ = 45◦ :
∼$ sample -latestTime
! Results
You may find the following relation useful in converting the velocity from the
Cartesian base where it is computed, to the cylindrical base used in Ref. [3] to
display the results:
Ur sin(θ) cos(θ) Ux
= (5.8)
Uθ cos(θ) − sin(θ) Uy
Note that two ppUtil (cf. Section 4.6.2) are used in this tutorial, returning
both the global balance of ions (calcBalance) and the current density through the
cylinder surface (calcJpatch). The latter allows to verify, in this specific case, that
the no-flux boundary conditions for the two ionic species is working as expected,
since a quasi -null current density is retrieved at the boundary where it is assigned.
reservoir
E
cyc0 H cyc1
y
membrane
6H
! Boundary conditions
The boundary conditions were described in detail in Ref. [3]. For the users
of OpenFOAM
R
v4.x, if the fixedFluxExtrapolatedPressure is intended to be used
at the membrane boundary – this is the most correct option –, then the function
adjustPhi() must be disabled in solver rheoEFoam, as discussed in Section 4.5.8
and in the previous tutorial (again, because no Dirichlet BC is used for pressure,
which automatically enables function adjustPhi()). Since this would require a
permanent change in the solver, we have chosen to use a zeroGradient BC in
replacement. The degree of approximation, as well as the differences between the
two approaches, are the same as for the tutorial in the previous Section (see the
comments therein).
! Command-line
This tutorial is slightly different from the previous ones regarding the
command-line sequence to run. This is because the problem is first solved in
a 1D configuration and the resulting solution is then disturbed and used as the
starting solution of the 2D configuration. Thus, we recommend to use directly the
Allrun script to run this tutorial, although we also explain next the main steps
accomplished by that script.
Firstly, in the directory of the main case you will find a folder named so
lution1D/. This is where the 1D problem is solved – script Allrun inside
this folder is the first call of the Allrun in the main directory. Computing the
1D solution is relatively straightforward in what respects the commands to be
executed, since it only requires building the mesh (blockMesh) and running the
solver (rheoEFoam). You may note that the hydrodynamic component of the
solver is switched off through the solveFluid keyword in dictionary fvSolution,
CHAPTER 5. Tutorials 146
which is set to false. Thus, only the Poisson and Nernst-Planck equations are being
solved. Furthermore, we note that underelaxation is being used to compute the 1D
solution in order to enable the use of high time-steps. Under the current settings,
the lower the voltage on patch reservoir, the higher the time for the 1D solution to
converge (the endTime should be increased accordingly in these situations). Our
criteria for convergence relies in the monitor for the current density.
After the 1D solution is computed, the resulting fields are mapped to the
2D domain (also created with blockMesh beforehand), using the mapFields util-
ity available by default in OpenFOAM
R
. Then, the fields for the cationic
and anionic concentration are locally disturbed by a 1 % random perturbation.
This is accomplished by a pre-processing utility named rndPerturbation, which
has been specifically created for this task and that can be found in directory
src/libs/preProcessing/rndPerturbation/. The case is ready to be
run with rheoEFoam, noticing that now solveFluid = true.
While the case is running, a ppUtil (cf. Section 4.6.2) is simultaneously being
executed (calcJpatch), which retrieves the surface-averaged current density over
time on both the reservoir and membrane patches (Fig. 5.24).
! Results
The current density can be plotted over time and the contours of charge density
can be computed and visualized in Paraview (these are just some suggestions).
The electric field intensity can be controlled by changing the voltage in file
solution1D/0/psi, under the entry for patch reservoir. Note that any change
in the initial/boundary conditions of any field must be done in the files inside folder
solution1D due to the mapping procedure. By running the Allrun script in
the main folder, the 1D solution will be always computed first.
inletNorth
walls walls
15H
y
intletWest 2H outlet
x
7H
walls walls
7H frontAndBack
0.8H
inletSouth 2H
The geometry is divided into 12 blocks for meshing purposes. The dimen-
sions of the cells at the junction corners are approximately: (∆x, ∆y, ∆z) ≈
(0.1, 0.1, 0.07)H.
! Boundary conditions
Since the Ohmic model is used in this simulation, we have to define boundary
conditions for pressure, velocity, electric potential and conductivity. We consider
that all the arms are open to the atmosphere, thus p = 0 and a zero-gradient is
assumed for the velocity. At inletNorth and inletSouth, the conductivity of the
fluid entering those arms is the same, but ten times lower than the conductivity of
the fluid entering inletWest (γ = 10). A passive scalar (neutral dye) is also entering
the geometry through inletWest, only. A zero-gradient condition is imposed for
both the conductivity and the passive-scalar concentration at the outlet. At the
walls we assume zero-gradient for all variables, except the velocity. As discussed
in Section 3.6.6, the Ohmic model is usually complemented with a conductivity-
dependent slip velocity (Eq. 3.42). In this tutorial, we use m = −0.3 for the
CHAPTER 5. Tutorials 148
power-law index (see Section 4.5.7). Regarding the electric potential, the outlet
is grounded, while the potential at inletNorth and inletSouth are the same, being
adjusted in order to impose the desired RaE . The potential at inletWest is set
according to β.
Note that the passive-scalar (dye) has a diffusivity which is one order of magni-
tude lower than the diffusivity of the ions (conductivity), in agreement with typical
real conditions.
! Command-line
1–Build the mesh:
∼$ blockMesh
2–Create fields C and sigma by copying the ones already present, but which are
not initialized in the interior domain:
∼$ cp 0/C.org 0/C
∼$ cp 0/sigma.org 0/sigma
3–Initialize fields C and sigma in the interior domain (C = 1 ∧ σ = σW for
x < −H and t = 0):
∼$ setFields
4–Run the solver:
∼$ rheoFoam
5–Post-process in Paraview (hint: slice the channel at the midplane in the z -
direction and plot the contours of C):
∼$ paraFoam
! Results
The contours for the passive-scalar are displayed in Fig. 5.26 at the mid-plane
in z, for different RaE . The reader will probably note similarities between these
instabilities and von-Kármán vortex streets.
The patterns are qualitatively similar to those obtained in Ref. [49] (see Fig. 4
therein, although our RaE is defined differently from the one used in that work).
The periodicity or the chaotic behavior of the instabilities can be further evaluated
by looking to the probes of C and U.
CHAPTER 5. Tutorials 149
RaE = 481
RaE = 579
RaE = 629
RaE = 839
RaE = 1424
RaE = 2164
RaE = 6239
Note that we do not used exactly the same parameters as in Posner et al. [49],
since our purpose was solely to illustrate the qualitative behavior of electrokinetic
instabilities, through a fast-running case. The large time-step used (Courant-
controlled, with Co = 1.5) is also inadequate to capture accurately the transient
behavior of this case.
inletNorth
y
walls
inletWest outlet
x
symmetry
! Boundary conditions
We consider a purely EDF, thus p = 0 at both inlets and outlet boundaries, and
a zero-gradient condition is simultaneously assigned to the velocity. Regarding the
CHAPTER 5. Tutorials 151
(a)
(b)
Figure 5.28: Snapshots of the passive scalar concentration field at different in-
stants: (a) injecting the phase devoid of C (inletNorth) and (b) injecting the phase
rich in C (inletWest). Note that the geometry has been reflected relative to the x -
axis (in Paraview) for presentation purposes – only the upper-half of the geometry
is simulated.
The user can play with the several degrees of freedom of the electric potential
boundary conditions in order to achieve different degrees of mixing. Note that the
time-step used in the tutorial has been adjusted to obtain results in a reasonable
amount time, but it should be lowered to achieve higher accuracy in time. Such
a high time-step, as well as the ”coarse” mesh used, would not be possible to use
if, for instance, the full PNP model was used instead – the simulation would most
likely diverge due to stability issues.
! Overview
The electrically-driven flow of high-molecular weight polyacrylamide solutions
in cross-slot and flow-focusing devices was seen to become unstable after a thresh-
old electric potential is exceeded [51]. This tutorial addresses the numerical simula-
tion of such phenomena, as presented in Ref. [51]. This tutorial merges electrically-
driven flows with viscoelastic fluid models, where the electric charge distribution
is computed by the Poisson-Boltzmann model and the extra-stress tensor of the
viscoelastic fluid is evolved using the Oldroyd-B model (β = 0.4). The settings
of the tutorial reproduce the case with ∆V = 160 V, W iB = 2.06, W iκ = 103
(CrossSlot/) and ∆V = 160 V, W iB = 1.03, W iκ = 154 (FlowFocusing/)
of Ref. [51]. The physical time of the simulations is tf = 0.5 s = 10λ.
! Geometry & Mesh
The geometry for this case is displayed in Fig. 5.29. Both the geometry and
the mesh correspond to the ones used in Ref. [51]. The flow is assumed to be 2D,
being solved in the xy-plane. Note that the patch names for the cross-slot and
flow-focusing geometries differ in the west arm, which is either an outlet or an
inlet (Fig. 5.29).
inlet_north
walls
walls
0.2H
y
outlet_west
2H outlet_east
inlet_west x
49H walls
walls
49H
inlet_south
Figure 5.29: Cross-slot and flow-focusing geometries. In the west arm, the black
arrow is for the cross-slot configuration (outlet west) and the purple arrow is for
the flow-focusing configuration (inlet west). In the tutorial, H = 5 × 10−5 m.
! Boundary conditions
CHAPTER 5. Tutorials 154
The boundary conditions used are described in Ref. [51]. The tutorials also
include the transport of a passive tracer. Note that the pressure extrapolation
boundary condition at the wall is replaced by a zero-gradient condition in version
fe40, since it is not available therein.
! Command-line
Before presenting the command line sequence, it is worth to note that the
mesh for this case is built in a slightly different way. Indeed, the blockMesh
application only builds one-quarter of the geometry/mesh, the one in the first
quadrant (+,+). The remaining of the geometry/mesh is built by 2 sequential
mirroring operations (mirrorMesh): first using the Oyz plane and then using
the Oxz plane. Afterwards, patches should be renamed accordingly, which requires
selecting, spliting and removing faces of already existing patches (topoSet) and
creation of new ones (createPatch).
1–Build the mesh:
∼$ blockMesh
∼$ cp system/mirrorMeshDict0 system/mirrorMeshDict
∼$ mirrorMesh
∼$ cp system/mirrorMeshDict1 system/mirrorMeshDict
∼$ mirrorMesh
∼$ topoSet
∼$ createPatch -overwrite
2–Initialize field C in the interior domain:
∼$ cp 0/C.org 0/C
∼$ setFields
3–Run the solver:
∼$ rheoEFoam
! Results
The dye patterns (field C ) of the unstable flows can be visualized in Paraview.
5.5 rheoBDFoam
5.5.1 General guidelines
| constant/
CHAPTER 5. Tutorials 155
• subDict outputOptions
writeStats – this option enables writing the index, position and stretch of all the molecules
over time (see Section 4.3.8). (boolean)
outputStats
– number of time-steps between consecutive calls to the statistics writer. It
Interval
is independent from the writing of time directories. Only effective if the
previous option is enabled. (integer)
writeVTK – this option enables writing the molecules’ structure and properties in VTK
format, which should be readable by any version of Paraview. (boolean)
• subDict exclusionVolumeProperties
activeExclusion
– this option allows to include or suppress exclusion volume forces between
Volume
beads (see Section 3.7.2). (boolean)
activeWall
– this option enables to impose a minimum approximation distance of the
Repulsion
beads to the walls, upon collision. If not enabled, the beads are repositioned
at the collision point, inside the computational domain (see Section 4.3.7).
(boolean)
repulsiveDistance – if activeW allRepulsion = true, this entry defines the repositioning distance
of the beads in the wall-normal direction, starting from the collision point
on the wall (see Section 4.3.7). (double)
• subDict HIProperties
• subDict electrophoresis
mobility – if active = true, this entry should correspond to the electrophoretic mobility
of each individual bead. (double)
• subDict springModelProperties
springModel – this entry should correspond to one of the available spring models (see Section
4.3.5). (string)
timeScheme – this entry should correspond to one of the available time integration schemes
(see Section 4.3.5). (string)
relTol – minimum normalized spring length variation below which the Newton-
Raphson method stops. (double)
tresholdF – fraction of the maximum spring length which triggers the algorithm to add
implicitly the diagonal spring force terms contribution. This is the α variable
defined in Section 3.7.4, which should take values in the range ]0, 1]. (double)
cutOffStretch – variable used to stabilize the Newton-Raphson method. In one given itera-
tion of the Newton-Raphson method, it may happen that a spring exceeds
its maximum admissible length (l ). As such, when the springs length needs
to be used to compute fk and Jk (see Section 3.7.4), we replace its value
by min (Ri , cutOf f Stretch × l). However, the resulting beads positions are
never corrected. In order to ensure that this stabilization method will not
affect the simulation results, this variable should always be higher than the
maximum expected fractional extension and ≤ 1. In most of the cases, a
value between 0.990 and 0.999 will be a good choice (FENE and Cohen Padé
models will typically require higher values, as 0.999 or even 1). (double)
solver – this entry should correspond to the name of one of the available matrix solvers
to be used in the Newton-Raphson method (see Section 4.3.5). (string)
externalFlow
2 {
writeFields false;
4 frozenFlow true;
6 tethered false;
pushBackCmp (0 0 0);
8 pushBackFreq 1;
10 interpolation BarycentricWeights;
gradU
12 (
0 0 0
14 0 0 0
0 0 0
16 );
}
18
outputOptions
20 {
writeStats true;
22 outputStatsInterval 10;
24 writeVTK true;
}
26
exclusionVolumeProperties
28 {
// bead-bead
30 activeExclusionVolume true;
32 // wall-bead
CHAPTER 5. Tutorials 158
activeWallRepulsion true;
34 repulsiveDistance 1e-7;
}
36
HIProperties
38 {
activeHI true;
40 }
42 electrophoresis
{
44 active false;
mobility 5.95767e-10;
46 }
48 springModelProperties
{
50 springModel MarkoSiggia;
52 timeScheme semiImplicit;
maxIter 20;
54 relTol 1e-6;
tresholdF 0.85;
56 cutOffStretch .99;
solver QR;
58 }
Listing 5.5: Example of a moleculesControls dictionary used with
rheoBDFoam.
| 0/
As we have seen for directory constant/, a simulation with rheoBDFoam
requires all the elements that would be needed for a simulation with rheoEFoam
(Section 5.4.1). Therefore, continuum fields p and U must always be present, and
phiE and/or psi should also be defined for an electrically-driven flow, or if molec-
ular electrophoresis is accounted for. Note that in case molecular electrophoresis
is considered, only the electric field represented by phiE (the externally applied
electric field) is used. Consider now three typical situations where rheoBDFoam is
likely to be used:
I- steady numerical external forcing. In this case, the continuum fields defined
in the startTime directory will be used and kept constant during the simu-
lation. The continuum is frozen in time. The usual procedure in this scenario
is to simulate the continuum flow elsewhere, in a separate simulation/case,
and then copy the steady-state continuum fields to the startTime directory
of the Brownian dynamics case.
II- transient numerical external forcing. In this situation, the continuum evolves
in time simultaneously with the molecules. Therefore the continuum fields
in the startTime directory should include initial and boundary conditions.
III- steady analytical external forcing. Although in this case the external forcing
is defined by an analytical expression, the continuum fields that would be
CHAPTER 5. Tutorials 159
needed for a numerical forcing should also be present. However, they are just
placeholders and are not used, such that their content is not important, but
needs to respect the expected syntax (have defined internal and boundary
fields compliant with the boundary types).
In addition to files for the continuum fields, the startTime directory still
needs to have a subdirectory lagrangian/molecules/ containing the fields
for the molecules. The easiest and recommended way to generate this subdirectory
and the files needed is running utility initMolecules, as described in Section 4.6.4.
Note that if any viscoelastic fluid model is used in the continuum flow, which
should probably be unusual, then also the corresponding variables should be de-
fined in the startTime directory (see Section 5.1.1). No additional files are
needed for a (non-elastic) generalized Newtonian fluid model.
| system/
Again, directory system/ for a rheoBDFoam case should have the same el-
ements as any rheoEFoam case. Moreover, most of the dictionaries do not re-
quire any change. The only modifications (additions) needed are in dictionary
fvSolution, in subDict SIMPLE :
solveFluid – this entry allows to switch on/off the resolution of the constitutive, momen-
tum and pressure equations. This option is typically set as false when the
external forcing is analytical, or when the external forcing is numerical but
frozen in time. (Switch)
solveElecM – this entry allows to switch on/off the resolution of the electric-related equa-
tions, independently of the option selected for solveFluid. (Switch)
The ideal scenario would be to have both the mesh and the beads represented
together. Therefore, open (F ile → Open) again the same filename.OpenFOAM file
and select only the internalMesh under Mesh Parts. Depending on the geometry
of the case, you may be able to slice the mesh behind the position of the molecules,
ending with the representation of both. An alternative to the slice is to extract the
feature edges of the geometry, whereby you will obtain the external edges of the
geometry and be able to see the interior. The options are numerous. Nevertheless,
this method is not very fruitful in information, since the only ”interesting” field
associated to each bead is the molecule index (molcID) to which they belong.
The most one can do is to filter the beads by the molecule index, and isolate all
the beads of a given molecule. However, there is no information on the beads
connectivity, neither on the molecules’ group, which may be relevant variables.
There is still an additional issue: the position file associated to the Lagrangian
particles will not be read by old Paraview versions, thus the beads cannot be
visualized therein.
In order to solve these issues, rheoBDFoam can optionally output in runtime
the molecules data in VTK format. This option can be enabled in dictionary
moleculesControls, as seen in the previous Section, and it will create a di-
rectory named VTKMolecules, to which the molecules’ data is written for each
outputTime. When the molecules are visualized in VTK format, we recommend
also to convert the mesh and continuum fields to VTK format by running the
default OpenFOAM
R
utility foamToVTK. Then, simply open the two stacks of
VTK files in Paraview. If using the paraFoam application, delete the object auto-
matically created in the Pipeline Browser and F ile → Open the two VTK stacks
(a stack contains the files generated for all the time-steps, under the same prefix
name). Note that the stack directly sent to directory VTK by OpenFOAM
R
is the
one for the internal field data. If you open that one, then the slicing procedure
described above is also needed. The fields associated to each bead are globalID,
groupID, localID and molcID. Following a top-down approach, filters can be used
(Treshold in particular) to select a group of molecules (e.g. groupID = 0), a par-
ticular molecule of the group (e.g. molcID = 0) and finally the k th bead in that
molecule (localID = k). You can also color the beads of a given molecule by their
localID to see the connections. The Glyph filter may help to obtain nice images
(Fig. 5.30). In general, the VTK files will be read by any Paraview version.
Summarizing, we recommend using the fist method to routinely check a case
being prepared for simulation and to monitor the simulation. The second method
should be used to debug, to obtain good-quality images, or if an old Paraview
version is to be used.
stretching of λ−DNA molecules in a planar extensional flow was selected for such
purpose, since this is a well-studied case. Before imposing the planar extensional
flow, the molecules are initialized in a stretched state and left to relax in no-flow
conditions. This allows one to estimate the relaxation time of the molecules and
also provides the starting molecular configurations for the planar extensional flow.
The λ−DNA molecules are modeled with the set of parameters provided in [21]:
N = 11 beads, D = 0.065 µm2 /s, a = 0.077 µm, Nk,s = 19.81, ν EV = 0.0012 µm3
and Lmax = 21 µm. The planar extensional flow is simulated for W i = λ˙ = 2,
where ˙ is the extensional-rate and λ is the relaxation time of λ−DNA molecules
(λ ≈ 4.1 s, as determined in the first part of the tutorial). The simulations are
run for an ensemble of 1000 molecules using the Marko-Siggia spring model.
! Geometry & Mesh
Since the forcing is analytical, a simple single-cell mesh can be used. The
dimensions are set large enough so that the molecules can be fully contained inside
the computational domain, even if fully-stretched. Furthermore, the molecules’
center of mass is held at a fixed position. The flow is unconfined.
! Boundary conditions
Due to the analytical nature of the forcing, any arbitrary, but consistent,
boundary condition can be attributed to the continuum fields (pressure and veloc-
ity); however, they will not be used. Consistent means here that if we assigned a
patch type to a given boundary while building the mesh, then we cannot assign,
for example, a symmetryPlane boundary condition to that boundary.
! Command-line
As aforementioned, in the first part of the tutorial the molecules are initialized
in a stretched configuration (L = 0.7Lmax ) in order to study their relaxation to
equilibrium in no-flow conditions. This is carried out inside directory relaxati
on/:
CHAPTER 5. Tutorials 162
! Results
Firstly, we will estimate the relaxation time of the molecules. It can be ob-
tained from the ensemble average relaxation curve output to file relaxation
/rheoToolPP/0/moleculesStats/G1/Stretch_Naverage.txt by util-
ity averageMolcN. The molecular length should decay over time as [20],
The evolution of the fractional molecular length over time, obtained in the
planar extensional flow at W i = 2 (second part of the tutorial), is plotted in Fig.
5.31b, along with the numerical results of [21].
Experimenting new Wi is as easy as changing the gradU entry in constant
/moleculesControls (the time-step may need adjustments). There is no need
to run again the initial molecules’ relaxation step. Also, the components of gradU
can easily be set to represent a shear flow.
0.8
-24
y = -0.2429x - 22.945
R² = 0.999
ln( L2 - L02 )
0.6
L/Lmax
-26
0.4
-28
0.2
-30 0
0 10 20 30 0 10 20 30
t (s) t (s)
(a) (b)
Figure 5.31: (a) Decay of the molecular length over time during relaxation in
no-flow conditions. The red line represents a linear fit to the data for 5 ≤ t ≤ 13
s. (b) Fractional molecular length evolution over time in a planar extensional flow
at W i = 2. The symbols represent the numerical results from [21].
to the ratio between the centerline velocity and the average velocity in the fully-
developed region of the north arm (valid for a pressure-driven laminar flow in
a straight channel with a square cross-section). Subscripts N and W point to
the North and West arms of the device. For the electrophoresis case, ˙ = UHN and
U = µE is the relation between the electric field magnitude and the electrophoretic
velocity magnitude. Both µ (the electrophoretic mobility) and E can be adjusted
interchangeably to impose a given U. Note that a given electrophoretic VR can
not be obtained by an equal ratio of electric potentials (exception is for V R = 1).
The simulations are run for an ensemble of 1000 molecules, divided over two
groups of 500 molecules. The Marko-Siggia spring model is used.
! Geometry & Mesh
The 3D flow-focusing geometry used is similar to the one in the tutorial of
Section 5.4.5 (see Fig. 5.25). However, the height considered in this example is
equal to the width, h = 2H, and H = 100 µm.
! Boundary conditions
Pressure-driven flow
The velocity magnitude in the north and south inlets is equal, and 20× higher
than the velocity magnitude in the west arm (V R = 20). The velocities are selected
in order to have W i = 20. The pressure is fixed at zero at the outlet and no-slip
boundary conditions hold at the walls.
Electrophoresis
For the electrophoresis case, the sole purpose of the continuum simulation is to
obtain the electric potential distribution. Thus, we need to define field phiE and
select any model for electrically-driven flows that allows one to compute Laplace’s
equation. Note that we are not interested in solving for the electroosmotic flow; we
only need the electric potential distribution, but this requires selecting a complete
model. The slipSmoluchowski model (see Section 4.2.1) is the one displaying the
simplest electricProperties dictionary to define, thus this one has been
selected. The electric field is assumed to be tangential to the walls, and the
potential at the inlets is fixed in order to have an electric field (velocity) ratio
≈ 20 between the north (south) and west arms. Then, we select µ in order to
match the U velocity corresponding to W i = 20 (we are not concerned with the
physical admissibility of µ in the tutorial).
! Command-line
The sequence of commands to run these cases is somewhat long. Thus, we
recommend to use directly the Allrun scripts inside directories pressureDr
ivenFlow/ and electrophoresis/. Each of these scripts will first run the
solver to obtain the steady-state continuum fields (either rheoFoam or rheoEFoam),
and then they run the Brownian dynamics solver (rheoBDFoam, inside directory
BDS/). Regarding the Brownian dynamics simulations, the molecules are first left
to equilibrate in the local fully-developed flow of the west arm (x = −4H) for
around 20λ ≈ 400s. This requires fixing the molecules center of mass (see Section
4.3.6). Then, this constraint is removed and the molecules flow freely until they
exit the geometry through the outlet patch. Both steps are carried out in the
CHAPTER 5. Tutorials 165
same directory (BDS/); the second simulation starts from the last time-step of the
first simulation. The second simulation stops when no molecule remains inside the
geometry. At the end of the simulation, we use utility averageMolcX (see Section
4.6.6) to extract the average molecular length over a line starting at x = −3H and
ending at x = 15H, partitioned in 100 bins. You can see that the application is
called for time = 400 because this is the startTime of the second simulation.
The comments in the Allrun scripts should help to understand all the steps.
Nevertheless, we would like to highlight two points that may seem less clear.
Firstly, the commands
∼$ foamDictionary -entry ...
simply use utility foamDictionary of OpenFOAM
R
to automatically change a given
entry in a given dictionary from the command-line (or via scripting). More infor-
mation can be found in the OpenFOAM
R
website. The second point is related
with the copy operations performed. For example, script pressureDrivenFl
ow/BDS/Allrun includes these lines,
∼$ cp -rf ../flowSimulation/0.01/p* 0/
∼$ cp -rf ../flowSimulation/0.01/U* 0/
which are simply intended to copy the steady-state (t = 0.01 s in this case) fields
U and p from the continuum simulation directory, to directory BDS/. Some lines
below, you will find,
∼$ cp -rf 0/p* 400/
∼$ cp -rf 0/U* 400/
which have a related function. These lines of code are executed after the first
Brownian dynamics simulation (equilibrating the molecules in the local flow) has
been run, and we call them to copy again the continuum fields U and p to the
startTime directory (in this case startTime = 400) of the next Brownian dy-
namics simulation, where they still do not exist, but are needed. The user may
be wondering why rheoBDFoam do not automatically wrote these fields during
the first run. This was because it was our option, by setting writeF ields = f alse
inside dictionary moleculesControls (see Section 5.5.1). This option was con-
sidered simply to save disk space (the advantage would be clearer if we would have
to save 1000 time-steps or if the mesh was much denser).
! Results
The fractional molecular length obtained for the pressure-driven flow and for
electrophoresis (V R = W i = 20) are plotted in Fig. 5.32. This data is the result
outputted by the call to averageMolcX. You can find it in rheoToolPP/400/
moleculesStats/G1/Stretch.txt and rheoToolPP/400/moleculesS
tats/G2/Stretch.txt for each of the two groups, and compute the average
between both groups. The pressure-driven flow seems to be slightly more effective
in stretching the molecules, eventually due to its higher flow heterogeneity over
the channel width. You can also visualize the molecules in Paraview using one of
the methods described in Section 5.5.1.
CHAPTER 5. Tutorials 166
It would be easy to run a case with the superposition of the two external
forcings, pressure-driven flow and electrophoresis. Moreover, we could even add
an electroosmotic flow. This exercise is left as a challenge.
Pressure-driven Electrophoresis
1
0.8
0.6
L/Lmax
0.4
0.2
0
-3 0 3 6 9 12 15
x/H
Figure 5.33: Fractional molecular length of λ-DNA in LAOE, for W i = 6.5 and
De = 0.45. The lines are for the experimental (dashed) and numerical (solid)
results in [53], while the symbols represent the numerical results obtained with
rheoTool .
The reader can easily try this, and check if the results are comparable, as
expected.
Appendix A
Table A.1: List of some relevant parameters and variables used by rheoTool
and correspondence with the nomenclature used in this guide. The list is not
exhaustive.
169
APPENDIX A. Parameters and variables in rheoTool 170
174
BIBLIOGRAPHY 175
[35] S. Xue and G. W. Barton, “An unstructured finite volume method for vis-
coelastic flow simulations with highly truncated domains,” Journal of Non-
Newtonian Fluid Mechanics, vol. 233, pp. 48 – 60, 2016. Papers presented at
BIBLIOGRAPHY 177
[38] M. Alves, F. Pinho, and P. Oliveira, “The flow of viscoelastic fluids past a
cylinder: finite-volume high-resolution methods,” Journal of Non-Newtonian
Fluid Mechanics, vol. 97, no. 2–3, pp. 207 – 232, 2001.
[46] C. Zhao and C. Yang, “An exact solution for electroosmosis of non-
newtonian fluids in microchannels,” Journal of Non-Newtonian Fluid Me-
chanics, vol. 166, no. 17–18, pp. 1076 – 1079, 2011.
BIBLIOGRAPHY 178
[53] Y. Zhou and C. M. Schroeder, “Single polymer dynamics under large ampli-
tude oscillatory extension,” Physical Review Fluids, vol. 1, p. 053301, 2016.