User Guide

User Manual:

Open the PDF directly: View PDF PDF.
Page Count: 203 [warning: Documents this large are best viewed by clicking the View PDF Link!]

User Guide
Version 4.0
April 26, 2019
This document is licensed under
Creative Commons Attribution-NonCommercial-NoDerivs 3.0 Unported License
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 no307499. 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.
This offering is not approved or endorsed by OpenCFD Limited, producer and
distributor of the OpenFOAM software via, and owner of the
and OpenCFD R
trade marks.
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
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 ( for more details.
Linux is a registered trademark of Linus Torvalds.
OpenFOAM is a registered trademark of of OpenCFD Limited.
Paraview is a registered trademark of Kitware.
Typeset in L
2016-2019 Francisco Pimenta, Manuel A. Alves
1 Introduction 1
1.1 Motivation................................ 1
1.2 Guideorganization ........................... 2
1.3 Changelog................................ 3
1.4 Citing rheoTool ............................. 7
1.5 Contacts................................. 7
1.6 Contributing............................... 7
2 Installation 8
2.1 Compatibility with OpenFOAM R
and foam-extend versions . . . . 8
2.2 Differences between versions . . . . . . . . . . . . . . . . . . . . . . 8
2.3 System requirements . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4 Step-by-step instructions . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4.1 Download/clone rheoTool ................... 9
2.4.2 Download Eigen library . . . . . . . . . . . . . . . . . . . . 10
2.4.3 Install Petsc library . . . . . . . . . . . . . . . . . . . . . . . 11
2.4.4 Compile rheoTool ........................ 13
3 Theoretical background 15
3.1 Governing equations of complex fluid flows . . . . . . . . . . . . . . 15
3.2 Stabilization of viscoelastic fluid flow simulations . . . . . . . . . . 16
3.2.1 The both-sides-diffusion (BSD) technique . . . . . . . . . . . 16
3.2.2 The log-conformation tensor approach . . . . . . . . . . . . 17
3.3 Coupling algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3.1 Pressure-velocity coupling . . . . . . . . . . . . . . . . . . . 18
3.3.2 Stress-velocity coupling . . . . . . . . . . . . . . . . . . . . . 19
3.4 High-resolution schemes . . . . . . . . . . . . . . . . . . . . . . . . 20
3.5 Movinggrids .............................. 21
3.6 Segregated vs coupled solvers . . . . . . . . . . . . . . . . . . . . . 22
3.7 Electrically-driven flow models . . . . . . . . . . . . . . . . . . . . . 22
3.7.1 Poisson-Nernst-Planck model . . . . . . . . . . . . . . . . . 23
3.7.2 Splitting the electric potential . . . . . . . . . . . . . . . . . 24
3.7.3 Poisson-Boltzmann model . . . . . . . . . . . . . . . . . . . 24
3.7.4 Debye-H¨uckel model . . . . . . . . . . . . . . . . . . . . . . 25
3.7.5 Slipmodel............................ 25
3.7.6 Ohmic (leaky dielectric) model . . . . . . . . . . . . . . . . 26
3.8 Brownian dynamics simulations . . . . . . . . . . . . . . . . . . . . 27
3.8.1 The bead-spring model . . . . . . . . . . . . . . . . . . . . . 27
3.8.2 Governing equations of beads motion . . . . . . . . . . . . . 30
3.8.3 Spring force models . . . . . . . . . . . . . . . . . . . . . . . 30
3.8.4 Time integration algorithm . . . . . . . . . . . . . . . . . . . 31
4 Overview of rheoTool 33
4.1 The constitutiveEquations library ................... 33
4.1.1 Available GNF and viscoelastic models . . . . . . . . . . . . 33
4.1.2 A note on FENE-type models . . . . . . . . . . . . . . . . . 41
4.1.3 Multi-mode modeling . . . . . . . . . . . . . . . . . . . . . . 43
4.1.4 Analysis of a code sample . . . . . . . . . . . . . . . . . . . 43
4.1.5 Advanced settings . . . . . . . . . . . . . . . . . . . . . . . . 49
4.1.6 Adding new viscoelastic or GNF models . . . . . . . . . . . 49
4.2 The EDFModels library ........................ 50
4.2.1 Available EDF models . . . . . . . . . . . . . . . . . . . . . 50
4.2.2 The potentials splitting approach and multi-species model-
ing in the PNP, PB and DH models . . . . . . . . . . . . . . 52
4.2.3 Electrokinetic coupling loop in the PNP model . . . . . . . . 52
4.2.4 Coupled PNP model . . . . . . . . . . . . . . . . . . . . . . 52
4.2.5 Analysis of a code sample . . . . . . . . . . . . . . . . . . . 53
4.2.6 Adding new EDF models . . . . . . . . . . . . . . . . . . . . 60
4.3 The BDmolecule library ........................ 61
4.3.1 Organization of variables . . . . . . . . . . . . . . . . . . . . 61
4.3.2 Solution sequence . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3.3 External forcing type . . . . . . . . . . . . . . . . . . . . . . 65
4.3.4 External forcing interpolation . . . . . . . . . . . . . . . . . 65
4.3.5 Spring force and time-integration schemes . . . . . . . . . . 68
4.3.6 Tethering and fixing the molecules center of mass . . . . . . 70
4.3.7 Beadstracking ......................... 71
4.3.8 Data output for post-processing . . . . . . . . . . . . . . . . 71
4.3.9 Limitations ........................... 72
4.4 The sparseMatrixSolvers library.................... 73
4.4.1 Conditions to reuse the preconditioner/factorization . . . . . 73
4.4.2 Residuals and tolerances . . . . . . . . . . . . . . . . . . . . 74
4.4.3 Generic parameters . . . . . . . . . . . . . . . . . . . . . . . 74
4.4.4 OpenFOAM interface . . . . . . . . . . . . . . . . . . . . . . 76
4.4.5 Eigeninterface ......................... 76
4.4.6 Hypreinterface ......................... 77
4.4.7 Petscinterface.......................... 80
4.4.8 Coupledsolvers......................... 82
4.4.9 How to use sparseMatrixSolvers library in my own application? 84
4.4.10 Limitations ........................... 85
4.5 Solvers.................................. 86
4.5.1 rheoFoam ............................ 87
4.5.2 rheoTestFoam .......................... 95
4.5.3 rheoInterFoam ......................... 97
4.5.4 rheoEFoam ........................... 98
4.5.5 rheoBDFoam .......................... 99
4.6 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.6.1 linearExtrapolation .......................101
4.6.2 navierSlip ............................101
4.6.3 zeroIonicFlux ..........................102
4.6.4 boltzmannEquilibrium ......................102
4.6.5 inducedPotential ........................103
4.6.6 slipSmoluchowski ........................103
4.6.7 slipSigmaDependent ......................103
4.6.8 A note on wall boundary conditions for pressure . . . . . . . 103
4.7 Utilities .................................105
4.7.1 GaussDefCmpw schemes for convective terms . . . . . . . . 105
4.7.2 Generic post-processing: ppUtil ................107
4.7.3 writeEfield ............................109
4.7.4 initMolecules ..........................110
4.7.5 averageMolcN ..........................114
4.7.6 averageMolcX ..........................114
5 Tutorials 116
5.1 rheoFoam ................................116
5.1.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 116
5.1.2 A note on coded FunctionObjects ...............122
5.1.3 Case 1: flow between parallel plates . . . . . . . . . . . . . . 123
5.1.4 Case 2: lid-driven cavity flow . . . . . . . . . . . . . . . . . 125
5.1.5 Case 3: flow in a 4:1 planar contraction . . . . . . . . . . . . 126
5.1.6 Case 4: flow around a confined cylinder . . . . . . . . . . . . 129
5.1.7 Case 5: bifurcation in a 2D cross-slot flow . . . . . . . . . . 132
5.1.8 Case 6: blood flow simulation in a real-model aneurysm . . . 134
5.1.9 Case 7: viscous fluid damper (moving mesh) . . . . . . . . . 138
5.2 rheoTestFoam ..............................140
5.2.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 140
5.2.2 Case I: Herschel-Bulkley model . . . . . . . . . . . . . . . . 143
5.2.3 Case II: FENE-CR model . . . . . . . . . . . . . . . . . . . 143
5.3 rheoInterFoam .............................147
5.3.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 147
5.3.2 Case 1: impacting drop . . . . . . . . . . . . . . . . . . . . . 148
5.3.3 Case 2: planar die swell . . . . . . . . . . . . . . . . . . . . 150
5.4 rheoEFoam ...............................152
5.4.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 152
5.4.2 Case I: EDF of power-law and PTT fluids in a microchannel 155
5.4.3 Case II: induced-charge electroosmosis around a cylinder . . 159
5.4.4 Case III: charge transport across an ion-selective membrane 161
5.4.5 Case IV: electrokinetic instabilities in a flow-focusing device 163
5.4.6 Case V: electrokinetic mixer . . . . . . . . . . . . . . . . . . 167
5.4.7 Case VI: electro-elastic instabilities in cross-shaped geometries169
5.5 rheoBDFoam ..............................171
5.5.1 General guidelines . . . . . . . . . . . . . . . . . . . . . . . 171
5.5.2 Molecules visualization with Paraview . . . . . . . . . . . . 176
5.5.3 Case 1: λ-DNA extension in a planar extensional flow . . . . 177
5.5.4 Case 2: 7λ-DNA extension in a flow-focusing device . . . . . 180
5.5.5 Case 3: λ-DNA dynamics in LAOE . . . . . . . . . . . . . . 183
6 FAQs 186
Appendix A Parameters and variables in rheoTool 188
Bibliography 193
Chapter 1
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.
In [5] we implemented coupled solvers for electrically-driven flows, which can
be also used for pressure-driven flows. Moreover, rheoTool was interfaced to ex-
ternal libraries (Petsc, Hypre, Eigen) that widen the range of available (direct and
iterative) sparse matrix solvers.
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
CHAPTER 1. Introduction 2
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.
the log-conformation tensor methodology [6] is available for a wide range
of viscoelastic models. This minimizes the numerical instabilities frequently
observed for high Weissenberg number flows.
a stress-velocity coupling term can be selected on run time in order to avoid
checkerboard fields under specific conditions, such as in the simulation of the
Upper-Convected Maxwell (UCM) model in strong extensional flows.
high-resolution schemes for convective terms are available in a component-
wise and deferred correction approach, avoiding numerical instabilities (see
Ref. [2] for details). Additional schemes were added to the newly created
library, which are not available by default in the OpenFOAM R
a solver (rheoTestFoam) is provided to compute the relevant material func-
tions of each GNF/viscoelastic model included in the library. The user can
select any canonical flow to be tested (shear flow, extensional flow, etc.).
a number of models for electrically-driven flows is available and can be cou-
pled with any rheological model. Mixed pressure- and electrically-driven
flows are also allowed.
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).
a solver is provided for Brownian dynamics simulations of polymer molecules
in generic meshes. Molecules can be linear or circular and they can also have
branches. The external forcing can be steady or transient.
coupled and segregated solvers are available.
rheoTool can use sparse matrix solvers from Petsc, Hypre and Eigen.
the tool is provided with a user-guide (this document) and a selected set of
tutorials reproducing relevant benchmark or real-life flow problems.
rheoTool is available for both 1OpenFOAM R
and 2foam-extend versions.
1.2 Guide organization
The remainder of this guide is organized as follows:
Chapter 2describes the basic steps to install rheoTool.
CHAPTER 1. Introduction 3
Chapter 3provides a succinct overview of the theory behind the governing
equations being solved. More details can be found in Refs. [24,7].
Chapter 4presents an overview of the functionalities available in rheoTool,
and discusses technical details about the code implementation.
Chapter 5contains several tutorials, guiding the reader into the use of
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 4and 5are still mainly based on
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 4.0
Released on 04/04/2019.
Add: added interfaces to Petsc, Hypre and Eigen libraries, allowing the use
of their (direct and iterative) sparse matrix solvers. All the interfaces can
be used in parallel, except the one for Eigen. Only for the OpenFOAM R
Add: rheoTool needs Petsc as an extra dependency. Added a script to install
this package. The instructions to install rheoTool have been updated (see
Chapter 2). Only for the OpenFOAM R
Add: added coupled solvers for both pressure- and electrically-driven flows.
Most viscoelastic models can be solved within a coupled solution method.
Only for the OpenFOAM R
Change: rheoTool version compatible with OpenFOAM R
v4.0/4.1 is discon-
Electrically-driven flows
CHAPTER 1. Introduction 4
Add: added coupled Poisson-Nernst-Planck model to EDF models (Type-
Name = NernstPlanckCoupled). Only for the OpenFOAM R
Constitutive equations
Change: the eXtended Pom-Pom model implementation has been changed in
order to allow using its thermodynamically consistent version [8]. Parameter
nhas been added (see Table 4.1) and should be adjusted by the user (n= 1
for thermodynamic consistency and n= 0 otherwise).
Change: the code for the constitutive equations of viscoelastic models has
been modified to allow integration with coupled solvers. Only for the
Fix: fixed bug in rheoInterFoam for OpenFOAM R
version 6.0, which was
preventing post-processing (missing call to update()).
Change: tutorial rheoFoam/Cavity now uses sparse matrix solvers from
Eigen library and is 1.6 times faster. Only for the OpenFOAMR
Change: tutorial rheoFoam/Cylinder is now solved coupled, being 30 times
faster. Only for the OpenFOAMR
Change: tutorial rheoEFoam/ICEO/NernstPlanck is now solved with the
coupled implementation of the PNP model and is 30 times faster. The name
of the tutorial was changed to rheoEFoam/ICEO/NernstPlanckCoupled.
Only for the OpenFOAM R
Change: tutorial rheoEFoam/selecMembrane/NernstPlanck/solution1D is
now solved with the coupled implementation of the PNP model and is twice
faster. Moreover, under-relaxation is not needed anymore to avoid numerical
divergence. Only for the OpenFOAMR
Fix: fixed bug related to the old flag for stabilization methods in several
tutorials, which was aborting the runs.
Version 3.0
Released on 18/09/2018.
Brownian dynamics solver
Add: solvers, libraries, utilities and tutorials for Brownian dynamics simu-
lations of polymer molecules.
CHAPTER 1. Introduction 5
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
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-
Add: added rheoTool patch for OpenFOAM R
Add: added note in the user-guide (Section 2.4.4) about parallel compilation
of rheoTool.
Constitutive equations
Add: Papanastasiou regularization is now available for yield-stress GNF
models (Hershel-Bulkley/Bingham and Casson models).
Add: Casson model has been added to the library of constitutive equations.
Add: the Multi-Lambda Isotropic Kinematic Hardening (MLK-IKH) 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: Saramito’s elastoviscoplastic model has been added to the library of
constitutive equations. Both stress and log-conformation versions are avail-
Add: the Bautista-Manero-Puig (BMP) model has been added to the library
of constitutive equations. Both stress and log-conformation versions are
Change: implemented functions tauTotal() and tauTotalMF() in base classes.
Solver rheoTestFoam and the utility retrieving the wall shear-stresses were
modified accordingly.
Version 2.0
CHAPTER 1. Introduction 6
Released on 09/02/2018.
Electrically-driven flows
Add: solvers, libraries, utilities and tutorials for 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
Add: the (single-equation) eXtended Pom-Pom 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.
Change: the stabilization method in viscoelastic simulations has been made
general and run time selectable: none,BSD or coupling.
Change: a verification step has been added to the WhiteMetznerLog model
in order to prevent its incorrect use (see the note in the table displaying the
constitutive equations).
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
Change: (fvc::grad(U)&fvc::grad(etaS()*alpha)) has been replaced by
fvc::div(etaS()*alpha*dev2(T(fvc::grad(U)))) for the use in multi-
phase flows (constitutiveEq.C).
Fix: call to constrainPressure() in rheoInterFoam, version of40, has been
corrected for the SIMPLEC algorithm (pEqn.H). Added a section in the user-
guide on how to use properly the fixedFluxPressure BC with rheoInterFoam
in versions of222 and fe40.
Add: tutorials on the die swell problem.
CHAPTER 1. Introduction 7
Change/Fix: code cleanup and bug fix (BC evaluation of the explicit fvc::
div(phi,X) operator) in class GaussDefCmpw.
Change/Add: replace boundary condition extST by the Type-independent
linearExtrapolation boundary condition (no backward compatibility). Added
optional second-order regression.
Change: major update of the user guide to include electrically-driven flows.
Other changes were made in its content and organization, and some typos
were corrected.
Change: ensure compatibility with foam-extend 4.0 and OpenFOAMR
Version 1.0
Released on 6/12/2016.
Initial version.
1.4 Citing rheoTool
If you found rheoTool useful and want to cite it in your work, the following BibTex
entry can be used for that purpose:
author = "F. Pimenta and M.A. Alves",
title = "rheoTool",
howpublished = "\url{}",
year = "2016"}
Since the underlying theory of rheoTool has been mainly presented in technical
papers [25], 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:
RFrancisco Pimenta:
RManuel A. Alves:
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
2.1 Compatibility with OpenFOAM R
and foam-
extend versions
The development and testing of rheoTool is usually performed in the most re-
cent stable release of OpenFOAM R
. However, an effort has been made to keep
rheoTool also compatible with foam-extend, although not all features available for
the OpenFOAM R
version can be found there (see Section 2.2). Currently, we
provide versions of rheoTool for:
v6.0 (of60/).
foam-extend 4.0 (fe40/).
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
foam-extend versions not included in this list. The versions above were tested in
Ubuntu 16.04, 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.
2.2 Differences between versions
The complete version of rheoTool is the one available for OpenFOAMR
, since this
is the one used for development. There are two main components of this version
that are not available for the foam-extend version: the Brownian dynamics module
and the interfaces to external libraries (Petsc, Hypre and Eigen) providing access
to a wider range of sparse matrix solvers (including coupled solvers).
Besides those major differences, in order to make rheoTool compatible with
each OpenFOAM R
/foam-extend version, several modifications were required at
the programming level for each case. Still, the user-interface remained almost
unchanged among the different versions. The main exception is on the codedStream
FunctionObjects and coded boundary conditions, which are used in the tutorials
of Chapter 5. Indeed, while these functionalities are available in OpenFOAM R
CHAPTER 2. Installation 9
it is not the case for foam-extend. Thus, the coded boundary conditions and the
utilities implemented as codedStream FunctionObjects in OpenFOAM R
had to be assembled and compiled in a library for the foam-extend version.
A second point to be taken into account is that rheoTool may perform differ-
ently in each OpenFOAM R
/foam-extend version, as it may happen with any other
default solver of OpenFOAM R
/foam-extend. This is naturally a consequence of
the evolution of the core machinery of OpenFOAM R
/foam-extend, transversal to
many solvers and libraries. Fortunately, in most of the cases the differences will
be small. The following issues were detected in the tests that we performed:
in general, a tutorial of rheoTool for OpenFOAM R
versions may be run either
in serial or parallel while keeping the same numerical settings. However,
in the tests using the foam-extend version, it was observed that parallel
simulations are less stable than serial ones, usually requiring a smaller time-
step or some under-relaxation of the velocity (sometimes as low as 0.97).
2.3 System requirements
Only standard requirements are needed to install rheoTool:
a compatible and functional version of OpenFOAM R
or foam-extend should
be already installed.
the machine should be connected to the Internet.
users might additionally need sudo privilege to install Petsc.
This is discussed in more detail in Section 2.4.3.
2.4 Step-by-step instructions
After ensuring that the prerequisites are fulfilled, the user is ready to start the
installation, which includes four (three) major steps:
1. Download/clone rheoTool (Section 2.4.1);
2. Download Eigen library (Section 2.4.2);
3. (Only for OpenFOAM R
users) Install Petsc library (Section 2.4.3);
4. Compile rheoTool (Section 2.4.4).
2.4.1 Download/clone rheoTool
This step is common to both OpenFOAM R
and foam-extend users.
CHAPTER 2. Installation 10
rheoTool is publicly available as a GitHub repository. There is a single branch
(master) in the repository, which always contains the most recent, stable version of
the code (there is no public dev branch or similar). Once a new version is pushed
to the master branch, the previous, old version is tagged as a release. Therefore,
releases are only checkpoints for older versions and should not be used to get the
most recent version of the code.
As explained in the Readme file of the repository, rheoTool can be either down-
loaded or cloned from GitHub. The structure of rheoTool (
/fppimenta/rheoTool) is depicted in Fig. 2.1.
rheoTool ofxx
(…) (…)
Figure 2.1: Directory organization of rheoTool.
The top-level directory of rheoTool contains the versions available for different
(of ) and foam-extend (fe) versions. 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 directory 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.
2.4.2 Download Eigen library
This step is common to both OpenFOAM R
and foam-extend users.
In the directory corresponding to a given rheoTool version (e.g. of60, where
you will find file downloadEigen), open a terminal and check that file etc/bashrc
of your installed OpenFOAM R
or foam-extend version has been sourced. This is
particularly relevant if you have defined alias for different versions of OpenFOAMR
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
CHAPTER 2. Installation 11
command line (thus, $icoFoam -help means that you only type icoFoam
-help). If this check is successful, run the script downloadEigen in that terminal:
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
$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
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.
2.4.3 Install Petsc library
This step is only for OpenFOAM R
Petsc [911] has been added to the dependencies used by rheoTool starting
from version 4.0. This library provides a number of efficient and scalable sparse
matrix solvers, that can be used in rheoTool with both segregated and coupled
The script installPetsc is responsible for downloading Petsc and the de-
pendencies it needs, configuring Petsc, compiling the libraries and exporting the
variables needed to call and use Petsc from any location. While the script is run-
ning, the user will be prompted to insert its password (and Yes or No) to allow
the installation of some dependencies needed by Petsc:
Some of these dependencies might already exist in the system, and in those
cases nothing is done. Only the inexistent dependencies will be downloaded and
CHAPTER 2. Installation 12
installed. This is the only part of the script requiring sudo mode (password con-
firmation). If the user is sure that all these packages are already installed in
the system, then commenting section ## Install dependencies inside the script is
enough to avoid the need for sudo, which can be problematic for non-administrator
The script has been tested on fresh installs of Ubuntu 18.04 and 16.04, after
installing OpenFOAM R
binaries. The system OpenMPI is being used by both
and Petsc in such conditions. This is a point worth to emphasize:
both OpenFOAM R
and Petsc should be compiled with the same MPI software,
which, however, can be different from OpenMPI. At the beginning of script inst
allPetsc, the user can change the paths to the MPI library and wrappers upon
Petsc library is configured with standard options, that the user can found and
modify under section ## Configure petsc. To check the full range of options
available, please consult Petsc references [9,10]. The configuration that we set as
default will download and install the following additional packages from within
Note that Petsc could be also installed via apt-get (there are binaries available),
but this would not allow configuring and using the most recent versions of the
For most of the users, the script to install Petsc can be run without any mod-
ification. In such cases, go the directory corresponding to one of the rheoTool
versions (where you will find file installPetsc), open a terminal and ensure that
file etc/bashrc of your installed OpenFOAM R
version was sourced (as for Eigen).
Then run the script installPetsc in that terminal:
Petsc library is saved and compiled to directory:
which can be changed by the user inside the script (before running it).
Installing Petsc only needs to be performed once per system. For example,
if you have multiple versions of rheoTool in the same machine, only a single in-
stall of Petsc is needed (recommended). Note that script installPetsc modifies
your ~/.bashrc file by appending Petsc variables to it (see the code at the end
of the script). If you have multiple versions of Petsc installed, or if you change
the location or version of Petsc, do not forget to manage these variables. Run-
ning (until completion) installPetsc multiple times will duplicate these variables
and possibly cause conflicts of paths. After installing Petsc with success, it is
good idea to check your ~/.bashrc and verify if the three variables PETSC DIR,
PETSC ARCH and LD LIBRARY PATH appended at the end are correctly de-
fined and not duplicated.
CHAPTER 2. Installation 13
2.4.4 Compile rheoTool
This step is common to both OpenFOAM R
and foam-extend users.
It is recommended to save rheoTool in a location with write permission, other-
wise you will need to use sudo mode to run all the commands. A good location for
rheoTool is, for example, directory $WM PROJECT USER DIR, which is defined
by default when OpenFOAM R
or foam-extend is installed.
Note (only for OpenFOAM R
users): in some combinations of OpenFOAM R
patches/OS versions, it was observed that the linker is unable to locate
PETSc, whose path is absent in variable $LD LIBRARY PATH. Indeed, echo
$LD LIBRARY PATH will not include the path to PETSc libs in such situa-
tions. This is most frequently observed when there are aliases to different
versions, since in these situations sourcing the etc/bashrc of
cleans $LD LIBRARY PATH. When this happens, running the
Allwmake script as described below results in linking errors similar to:
/usr/bin/ld: warning:, needed by /home/user/OpenFOAM/user-6/
platforms/linux64GccDPInt32Opt/lib/, not found (
try using -rpath or -rpath-link)
/home/user/OpenFOAM/user-6/platforms/linux64GccDPInt32Opt/lib/ undefined reference to ‘VecSet’
To avoid this issue, the user needs to ensure that the path to PETSc library is
added to $LD LIBRARY PATH after sourcing the etc/bashrc of OpenFOAM R
such that it survives. Among the several possibilities, a simple one is to modify the
alias. For that, open your system ~/.bashrc in edit mode, for example running
gedit /.bashrc, find your alias (typically at the bottom of the file) and
modify it appending command
at its end, with the two commands separated by a semicolon. Thus, an alias of
the form (this is just an example)
alias of60=’source /opt/openfoam6/etc/bashrc’
before the modification should look like this (single line)
alias of60=’source /opt/openfoam6/etc/bashrc; export LD_LIBRARY_PATH=$PETSC_DIR
after the modification. Note that there are other ways that would lead to the
same result, as for example modifying the $LD LIBRARY PATH variable at the
end of the etc/bashrc file of OpenFOAM R
, creating symbolic links of PETSc
libraries inside one of the locations included by default by OpenFOAM R
How do you know if the above procedure applies to your system? If you do not
have aliases for OpenFOAM R
versions, then you should not have this problem, as
CHAPTER 2. Installation 14
long as the command exporting $LD LIBRARY PATH is executed after sourcing
’s bashrc. If you do have aliases, then this problem might happen
or not. It will happen if command echo $LD LIBRARY PATH does not include
the path to PETSc libs, which you can test before compiling rheoTool. Otherwise,
the compilation error similar to the one displayed above should be suggestive of
this issue.
After you move rheoTool to its final location, open a new terminal (to ensure
that your system ~/.bashrc is sourced and contains the path of Eigen and Petsc)
in the top-level directory of rheoTool (ensuring that the OpenFOAM R
or foam-
extend environment has been sourced, as previously) and enter the directory with
the version of rheoTool that is compatible with your OpenFOAM R
or foam-extend
version, and then go to directory src/. For example, for OpenFOAMR
v6.0, it
would be:
$cd of60/src
Now, run the script Allwmake to build the libraries and applications of
rheoTool. In order to speed up the compilation, several processors can be used, if
available. For OpenFOAM R
users, run
$./Allwmake -j N
for parallel compilation with Nprocessors. 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 jis 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 jis not recognized by the script, therefore simply run
The compilation in foam-extend will typically make use of all processors avail-
able in the system, since variable WM NCOMPPROCS is set by default in such
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 running one of the
tutorials in Chapter 5.
Chapter 3
Theoretical background
The equations governing pressure- and electrically-driven flows of incompressible,
complex fluids are discussed in this Chapter, along with some important aspects
related with their discretization in the finite-volume framework. Since a thorough
discussion on this subject can be found in Refs. [2,3,5], some intermediate steps
are skipped and only the more relevant equations are presented.
The last Section of this Chapter is concerned with the theory of Brownian
dynamics simulations using coarse-grained models. Again, we will only present
the more relevant aspects for the scope of this user guide and more details can be
found in Ref. [4].
3.1 Governing equations of complex fluid flows
The basic equations governing isothermal, single-phase, transient flows, under lam-
inar conditions, for incompressible fluids, in static grids, establish mass conserva-
tion (Eq. 3.1) and momentum balance (Eq. 3.2),
∇·u= 0 (3.1)
t +u·u=−∇p+∇·τ0+f(3.2)
where uis the velocity vector, tis the time, pis the pressure, τ0is the extra-stress
tensor and fis any external body-force, such as the electric force discussed in
Section 3.7. 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-
bution (τ), τ0=τ+τ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,
τs=ηs( ˙γ)(u+uT) (3.3)
f(τ)τ+λ( ˙γ)
τ+h(τ) = ηp( ˙γ)(u+uT) (3.4)
CHAPTER 3. Theoretical background 16
In Eqs. (3.3) and (3.4), ηsis the solvent viscosity, ηpis 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=1
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. [12] for a comparison between the
most popular techniques) and the ones used in rheoTool are addressed next.
3.2 Stabilization of viscoelastic fluid flow simu-
3.2.1 The both-sides-diffusion (BSD) technique
The both-sides-diffusion (BSD) is a technique already incorporated in the vis-
coelasticFluidFoam solver [1]. It consists in adding a diffusive term on both sides
of momentum equation (Eq. 3.2), with the difference that one of them (left-hand
side) is added implicitly, while the other one (right-hand side) is added explicitly.
Once steady-state is reached, both terms cancel each other exactly. Such method
increases the ellipticity of the momentum equation and, as such, has a stabiliz-
ing effect, mostly when there is no solvent contribution in the extra-stress tensor.
Incorporating the terms arising from the both-sides-diffusion in the momentum
equation, and making use of Eq. (3.3), then
t +u·u− ∇·(ηs+ηp)u=−∇p− ∇·(ηpu) + ∇·τ+f(3.5)
Note that the added diffusive terms are scaled by the polymeric viscosity (ηp),
which is a common choice in the literature (e.g. Ref. [12]), 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 17
3.2.2 The log-conformation tensor approach
The log-conformation tensor approach consists in a change of variable when evolv-
ing in time the polymeric extra-stress and it was devised to tackle the numerical
instability faced at high Weissenberg number flows [6,13].
The polymeric extra-stress tensor is related with the conformation tensor (A).
For the Oldroyd-B model, for example, this relation is expressed as (see Table 4.1
for several viscoelastic models)
λ(AI) (3.6)
In the log-conformation tensor methodology, a new tensor (Θ) is defined as
the natural logarithm of the conformation tensor
Θ= ln(A) = Rln(Λ)RT(3.7)
In Eq. (3.7), the conformation tensor was diagonalized (A=RΛRT) because it
is positive definite, where Ris a matrix containing in its columns the eigenvectors
of Aand Λis a matrix whose diagonal elements are the respective eigenvalues
resulting from the decomposition of A. Eq. (3.4) written in terms of (Θ) becomes
t +u·Θ=ΩΘ ΘΩ + 2B+1
λg(Θ) (3.8)
where g(Θ) is a model-specific tensorial function depending on Θ(see Table 4.1
for other viscoelastic models) and
mxx 0 0
0myy 0
0 0 mzz
0ωxy ωxz
ωxy 0ωyz
ωxz ωyz 0
mxx mxy mxz
myx myy myz
mzx mzy mzz
ωij =Λjmij +Λimji
After solving Eq. (3.8), Θis diagonalized in the form
and the conformation tensor is recovered by the inverse relation of Eq. (3.7)
A= exp(Θ) = Rexp(ΛΘ)RT(3.14)
CHAPTER 3. Theoretical background 18
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=RuTζDRT.
It is worth to mention that the log-conformation approach can be considered a
particular case of the kernel-conformation method [14]. 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 rootkkernel [14] and the square-root
transformation [15] are also included in rheoTool for demonstration purposes.
3.3 Coupling algorithms
3.3.1 Pressure-velocity coupling
Although the OpenFOAM R
toolbox is already able to solve linear systems of
equations in a coupled way, most of the solvers still rely on segregated solutions
(this is a rule for transient solvers). In segregated solvers, the equations for each
variable are solved sequentially. Even for a fully-implicit method, if the coupling
between variables is weak, then numerical divergence is prone to occur.
In the OpenFOAM R
toolbox, common algorithms for pressure-velocity cou-
pling are SIMPLE and SIMPLEC for steady-state solvers and either PISO or
PIMPLE (a combination of SIMPLE(C) and PISO) for transient solvers. From
the benchmark cases performed in Ref. [2], it was observed that SIMPLEC was
particularly suitable for transient viscoelastic fluid flows at low Reynolds numbers,
regarding stability and accuracy.
The continuity equation, implicit in the pressure variable, derived for SIM-
PLEC (a more detailed derivation is presented in Ref. [2]) leads to
where aPare the diagonal coefficients from the momentum equation, H1=P
is an operator representing the negative sum of the off-diagonal coefficients from
momentum equation, H=P
nb +bis an operator containing the off-
diagonal contributions, plus source terms (except the pressure gradient) of the
momentum equation and pis 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
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 19
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 [16].
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 aPcoefficients
of the momentum equation [17]. 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.
3.3.2 Stress-velocity coupling
Stress-velocity decoupling problems can arise for similar reasons as those described
for pressure-velocity: the cell-centered velocity loses the influence of the forces
(either polymeric extra-stress or pressure gradient) of its direct neighborhood (cells
sharing a face in common). This usually happens in the interpolation from cell-
centered to face-centered fields. In the case of polymeric extra-stresses, it is the
divergence term (∇·τ) in the momentum equation, when τis linearly interpolated
from cell centers to face centers, which can be responsible for the decoupling.
In Ref. [2], we described a new stress-velocity coupling method, where the
polymeric extra-stresses at face centers are computed as
τf=τf+ηphu|f+(u)T|fu|f+(u)T|fi (3.17)
where terms with an overbar are linearly interpolated from cell-centered values,
while the remaining velocity gradients are directly evaluated from the cell-centered
velocities straddling the face. When the definition of τfin Eq. (3.17) is inserted in
the momentum equation with the both-sides-diffusion terms already present (Eq.
3.5), then we obtain
t +u·u− ∇·(ηs+ηp)u=−∇p− ∇·ηpu+∇·τ+f(3.18)
where the term ∇·ηpuis a ”special second-order derivative” (different from
the laplacian operator of OpenFOAM R
), defined as the divergence of the velocity
gradient, where the velocity gradient at the faces is obtained by linear interpolation
of the velocity gradient evaluated on the cell centers. More details are presented
in Ref. [2], where it is shown that with mesh refinement Eq. (3.17) approaches
τf=τfand the additional terms cancel out. Note that when inserting Eq. (3.17)
CHAPTER 3. Theoretical background 20
in the momentum equation (resulting in Eq. 3.18), we drop the transpose velocity
gradients for simplicity, since continuity imposes ∇·uT=0.
3.4 High-resolution schemes
The discretization of convective terms within the finite-volume framework leads to
(u·φ) dV=X
φf(uf·Sf) = X
where φis a generic variable being advected, Sfis the face-area vector and Ffis
the volumetric flux crossing face f. While fluxes are known at the faces from the
Rhie-Chow-like interpolation (Eq. 3.16), φat face centers need to be interpolated
from known values at cell centers. OpenFOAM R
offers a wide range of schemes to
perform such interpolation, from upwind – an unconditionally stable scheme, but
only first-order accurate –, to central differences – a conditionally stable, second-
order accurate scheme. A good compromise between both extremes is provided
by High-Resolution Schemes (HRSs). When represented in a Normalized Variable
Diagram (NVD), several HRSs are piecewise-linear functions and can be defined
using the Normalized Weighting Factor (NWF) approach [18]:
where the following definitions hold
In Eq. (3.20), αand βare scalars specific to each HRS and they can be functions
of e
φC. Subscripts in Eqs. (3.21a,b) have the following meaning: for a given face,
cell C is the cell from which the flux comes (upstream), cell D (downstream) is the
cell to which the flux goes and cell U (far-upstream) is the cell upstream to cell C.
In a general unstructured mesh, cell U cannot be identified unequivocally, and φU
in Eqs. (3.21a,b) can be evaluated as [19]
φU=φD2(φ)C·dCD (3.22)
where dCD is the vector connecting the center of cells C and D. For a deferred
correction implementation of HRSs, the upwind part of the HRS is discretized
implicitly, while the remaining (difference between the HRS and the upwind differ-
encing scheme) is discretized explicitly (cf. Ref. [2]), which, using Eqs. (3.20-3.22),
results in
φf= [φC]implicit + [(α1)φC+βφD+ (1 αβ)(φD2(φ)C·dCD)]explicit (3.23)
Handling the HRSs in a deferred correction approach avoids, in some cases,
numerical instabilities introduced by the central-differencing component of the
CHAPTER 3. Theoretical background 21
HRS. Additionally, in Ref. [2] it was observed that the usual methodology of
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.
3.5 Moving grids
Some CFD problems require the simulation of a moving entity interacting with a
fluid. There are several approaches than can be used to tackle such problems, and
the choice is usually made based on a case-by-case analysis. Consider for example
the flow induced inside a sphere due to its time-dependent rotation. Such case can
be easily handled by defining adequate boundary conditions for the flow variables
on the sphere surface, without further modifications of the usual solver setup. On
the other hand, if we consider the time-dependent simulation of the flow inside an
axisymmetric stirred tank reactor, such approach is no longer adequate. Instead,
we can use, for example, a (rotating) non-inertial reference frame, which allows
to keep the mesh steady and introduces some acceleration terms in the momen-
tum equation (OpenFOAM R
allows the use of such non-inertial reference frames).
However, if the tank is not axisymmetric, the non-inertial reference frame becomes
useless and a different approach is needed. An immersed boundary method can
be used for that purpose, avoiding the use of moving grids. However, we will turn
our attention to moving meshes, i.e. a computational mesh whose control volumes
move in space over time.
For moving control volumes, the equations governing the flow need to be
changed regarding convective terms, which should account for the grid motion [20],
φf(ufub,f )·Sf(3.24)
where φis any generic variable being advected and ubis the velocity at which
surface Sis moving. Moreover, the space conservation law (SCL) needs to be
satisfied to ensure mass conservation [20],
ub·ndS= 0 (3.25)
CHAPTER 3. Theoretical background 22
If the SCL is ensured, the continuity equation remains unchanged and so does
the pressure equation. In practice, the SCL is imposed while computing RSub·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 [20].
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.
3.6 Segregated vs coupled solvers
The governing equations in an implicit CFD code can be solved either segregated
or coupled. In a segregated solution method, the equations are solved sequentially,
one at a time (equations for multidimensional variables are further split into com-
ponents). This is the standard method used in OpenFOAM R
and in most CFD
codes based on finite-volumes. In a coupled solution method, all the governing
equations are solved simultaneously. There are also semi-coupled solvers, which
lie somewhere between segregated and coupled solvers. In a semi-coupled solver,
part of the equations are solved coupled and part are solved segregated.
The segregated solution method has been and continues being a popular strat-
egy for its low computational cost per time-step and low memory usage, compared
to the coupled solution method. Nonetheless, they are less stable than coupled
solvers, usually requiring lower time-steps and/or more under-relaxation in order
to avoid numerical divergence. Thus, the higher usage of resources by coupled
solvers is sometimes compensated by its enhanced stability, which translates in a
lower total time of computation. Moreover, due to its higher implicitness, coupled
solvers can be also more accurate in transient flow simulations [5].
In [5] we discussed the implementation of coupled and semi-coupled solvers in
rheoTool, in the context of electrically-driven flows. Semi-coupled solvers proved
to be faster and more accurate (time accuracy) than segregated solvers in a number
of situations. Similar advantages could be also observed in pressure-driven flows.
3.7 Electrically-driven flow models
Consider now that the fluid under analysis is a weak electrolyte subjected to an
electric field. In such conditions, the momentum equation (Eq. 3.2) should include
the contribution from an electric body-force,
f=fE=∇·εEE kEk2
CHAPTER 3. Theoretical background 23
where Eis the electric field, ε=ε0εRis the electric permittivity and ρEis 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.
3.7.1 Poisson-Nernst-Planck model
In the absence of magnetic effects, the electric potential (Ψ) can be computed by
Gauss’ law
∇·(εΨ) = ρE(3.27)
where the electric field is E=−∇Ψin electrostatics. By definition, the charge
density is
where Fis Faraday’s constant, ziis the charge valence of specie iand ciis the
concentration of specie i(mol/m3). The sum is over the Ncharged species in
the electrolyte. The standard law governing the transport of charged species in a
weak electrolyte, under the action of an electric field and neglecting any reaction,
is embodied by the Nernst-Planck equation,
dt +u·ci=∇·(Dici) + ∇·
kT Ψ
| {z }
which closes the system of equations for an EDF. In Eq. (3.29), Dis the diffu-
sion coefficient, eis the elementary charge, kis Boltzmann’s constant and Tis
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, Diezi
kT ci(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
CHAPTER 3. Theoretical background 24
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 ciand Ψ. Due to the
convective term in Eq. (3.29), there is a two-way coupling between the PNP and
the momentum equations.
3.7.2 Splitting the electric potential
Before proceeding to the derivation of other EDF models, we introduce here a
useful approach to simulate EDF problems. In the PNP model, a single electric
potential variable has been used, Ψ. However, in certain situations this can pose
some difficulties when defining the boundary conditions to solve the Poisson equa-
tion. A common approach to avoid such issues is the decomposition of the electric
potential in two variables: the externally imposed electric potential, φExt, and the
intrinsic electric potential, ψ, such that Ψ=φExt +ψ[3]. Following this approach,
Gauss’ law is also decomposed in two equations,
∇·(εφExt) = 0 (3.30a)
∇·(εψ) = ρE(3.30b)
An additional simplification which can be used simultaneously with the split-
ting approach is to consider fE=ρEφExt in the momentum equation, i.e., the
intrinsic electric potential contribution is ignored in the electric field definition.
This can be justified by stating that this extra force not accounted for directly is
balanced by a pressure gradient, which mutually cancel each other in the momen-
tum equation [3], under the assumption that it would not affect the flow.
The splitting approach will be used in the derivation of the next two models.
3.7.3 Poisson-Boltzmann model
If we assume that the ions follow a Boltzmann equilibrium, then the PNP
model can be simplified to the so-called Poisson-Boltzmann model (henceforth
PB model), for which Gauss’ law reads
∇·(εψ) = F
zici,0exp ezi
kT (ψψ0)(3.31)
with ci,0being a reference concentration of specie i, where the intrinsic poten-
tial is ψ0. Without loss of generality, we will assume that ci,0is the bulk ionic
concentration, where the intrinsic potential is ψ0= 0.
Note that the right hand side of Eq. (3.31) represents (minus) the charge density
for the PB model. Thus, Eq. (3.31) provides the definition of Eq. (3.30b) for the
PB model, under the splitting approach.
CHAPTER 3. Theoretical background 25
For this model, the only electric-related unknowns are the two electric poten-
tials, ψand φExt, computed from Eqs. (3.30a) and (3.31). Furthermore, as can
be seen from Eq. (3.31), there is no influence of flow variables in the PB model
(one-way coupling).
In order to increase the implicitness of Eq. (3.31), its source term can be lin-
earized by expansion in Taylor series up to the first-derivative, transforming the
equation into
∇·(εψ) + ψF
with bi=ezi
kT and ai=zici,0exp (biψ). All the terms of Eq. (3.32) with a star are
evaluated explicitly.
3.7.4 Debye-H¨uckel model
Considering the PB model, if we further simplify Eq. (3.31) assuming low electric
potentials, ezi
kT ψ1 , then
∇·(εψ) = F
kT ψ(3.33)
which is the equation governing the electric potential distribution in the so-called
Debye-H¨uckel model (henceforth DH model).
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).
3.7.5 Slip model
A common characteristic of electrokinetic problems is the spontaneous formation
of an electric double layer (EDL) near a charged surface, upon contact with an
electrolyte. The thickness of the EDL can be approximated by the Debye length
(λD), a physical parameter appearing when solving the Poisson equation for the
electric potential,
F e
In several practical applications, the charge density is mainly located in the
EDL region, while the bulk electrolyte is neutral. If the Debye length is much
smaller than the characteristic dimension of the system (λD
W1) and assuming a
smooth, laminar flow inside the EDL, then it is possible to approximate the EDL
effect by a slip velocity at the surface, avoiding the need to solve the flow inside the
EDL. Such a case would be, for example, the pumping of a Newtonian electrolyte
(λD∼ O(109m)) in a microchannel of arbitrary shape (W∼ O(106m)), by
CHAPTER 3. Theoretical background 26
electroosmosis, at low voltage ( ez
kT ψ1) – the last conditions is usually relaxed.
The Helmholtz-Smoluchowski theory is frequently used to approximate the slip
velocity in such conditions,
uSch =µE(3.35)
where µ=εζ
η0is 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
W1 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
for complex fluid flows [21] – using a slip model would simply retrieve a smooth
flow in such cases.
3.7.6 Ohmic (leaky dielectric) model
The so-called Ohmic model [22] is particularly useful to simulate fluids of different
conductivities, although a generalized Ohmic model has been recently proposed
for different types of problems [23]. The model can be derived from the PNP
equations, rewritten in terms of the conductivity and free-charge density, and
assuming additionally instantaneous charge relaxation and electroneutrality [22].
The interested reader is directed to Ref. [22] for the full derivation of the Ohmic
model. Here, only the final equations are presented. Furthermore, and contrarily
to what was done for the previous models, we will restrict our analysis to a binary
electrolyte, i.e., an electrolyte composed of only one positive and one negative
species, with z+=z=z, but no restrictions in the relation between D+and
First, let’s start defining the conductivity (σ) and free-charge density (ρE) for
a binary electrolyte,
RT (D+c++Dc) (3.37)
ρE=F z(c+c) (3.38)
CHAPTER 3. Theoretical background 27
where Ris the universal gas constant. Imposing the conservation of each variable
leads to (after the assumptions mentioned above; more details in Ref. [22])
t +u·σ=Deff2σ(3.39)
∇·(σΨ) = 0 (3.40)
where the effective diffusivity is Deff =2DD+
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
In order to close the Ohmic model, the EDL effect is commonly represented
by a slip velocity, which avoids detailing the flow inside the EDL using a very
fine mesh. Since the zeta-potential of a surface depends generally on the ionic
conductivity, a σ-dependent slip velocity is typically used [22], such as
uSch(σ) = µ0σ
where µ0=εζ0
η0is a reference electroosmotic mobility, at a reference conduc-
tivity (σ0), and mis an exponent governing the power-law dependence of the
zeta-potential on the conductivity (m[0.5,0.3] is in agreement with several
works, e.g. [22]). Note that Ein Eq. (3.42) is the electric potential at the surface
where the slip velocity is computed.
3.8 Brownian dynamics simulations
In the previous Sections, polymeric fluid flows were addressed from a continuum
mechanics perspective. In this Section, we zoom-in the scale of analysis, such that
each polymer molecule is now modeled individually. We enter the kinetic the-
ory domain, which lays somewhere between continuum mechanics and atomistic
modeling. This means that even though each polymer molecule is simulated indi-
vidually, we ignore atomic-size events. Moreover, the models discussed here and
implemented in rheoTool neglect inter-molecular interactions, which is represen-
tative of dilute solutions. The interested reader is referred to [7] for a thorough
discussion on the kinetic theory of polymers.
3.8.1 The bead-spring model
The most commonly used coarse-grained models to simulate polymer molecules
are bead-spring and bead-rod models. Currently, only the bead-spring model
is available in rheoTool (Fig. 3.1). According to this model, the molecules are
CHAPTER 3. Theoretical background 28
represented by a set of Nbeads connected by NS= (N1) springs, for open
chains, or NS=Nsprings for closed chains. Each ibead owns a group denoted
as githat contains the index of all the beads to which it is directly connected to
by springs. For chains without branches, gihas one (beads at the edges) or two
elements, but more elements can be present for branched polymers, such as the
ones depicted in Fig. 3.1. Note that each bead and spring in the chain aims to
represent a group of atoms, and not an individual atom. In this guide, we use
either |xij|=|rjri|, j gi, or simply Rito denote the length of all the springs
associated with bead i.
Polymer molecules usually display a maximum contour length upon full-
extension (Lmax), which should be ideally reproduced by the numerical model.
Therefore, each spring also has a maximum length, l=Lmax/NS.
The minimum characteristic length featured in a spring is the so-called persis-
tence length (λP), below which the spring segments behave as rigid elements, with
fixed orientation. The Kuhn step size, defined as bk= 2λPis a measure commonly
used in coarse-grained models, especially in bead-rod models, where it represents
the fixed size of a single rod. For the physical representation of the springs to
remain valid, a minimum number of Kuhn steps should be used to represent a
(flexible) spring. The number of Kuhn steps per spring is denoted by Nk,s =l/bk
and is usually controlled by the choice of N.
In rheoTool, an individual molecule is represented by a set of beads and springs,
and a group of molecules is composed by an ensemble of molecules sharing the
same physical properties. Each simulation in rheoTool can handle simultaneously
several groups of molecules (Fig. 3.1).
CHAPTER 3. Theoretical background 29
Group of Molecules
Group of Molecules
Group of Molecules
RheoTool simulation
Figure 3.1: Polymer molecules representation by a bead-spring model. The orga-
nization levels used in rheoTool are also represented: beads and springs, molecules
and groups of molecules.
CHAPTER 3. Theoretical background 30
3.8.2 Governing equations of beads motion
Consider a chain with an arbitrary topology, composed by Nbeads. The time
evolution of the position vector corresponding to each bead (ri) is governed by
t =uf+
kT +6
where ufis a velocity imposed by an external forcing, Dij is the diffusion tensor,
kis Boltzmann’s constant, Tis the absolute temperature, Fjis the sum of spring
and exclusion volume (EV) forces (Fj=FS
j), ∆tis the discrete time-step, σ
is a tensor satisfying D=σσTand njis a vector whose 3 components are random
numbers uniformly distributed in the range [1; 1].
Dand σare (N×N) symmetric tensors, whose ij elements are themselves
(3 ×3) tensors. The single model available in rheoTool to represent the diffusion
tensor is the Rotne–Prager–Yamakawa (RPY) model [4],
Dij =
6πηaI=DI, i =j
|xij |h1 + 2a2
3|xij |2I+12a2
|xij |2xij xij
|xij |2i, i 6=j∧ |xij | ≥ 2a
Dh19|xij |
xij xij
|xij |i, i 6=j∧ |xij|<2a
where |xij|=|rjri|,ηis the fluid viscosity, ais the bead radius and Iis the unit
tensor (3 ×3). Note that in rheoTool,aand Dare defined independently by the
user and need not to be related. For the RPY tensor, Dij
rj=0in 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 Dbecome 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 [25],
exp 9
2Nk,s |xij|2
where νEV is the exclusion volume parameter.
3.8.3 Spring force models
Several models can be used to express the spring force acting on each bead. Firstly,
one should distinguish between the models that limit the maximum spring length
CHAPTER 3. Theoretical background 31
and the models that do not impose any restriction on the spring length. Among
the mostly used models, the Marko-Siggia, Cohen Pad´e 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 [24],
Hxij (3.46)
where H=3kT Nk,s
l2. As can be seen, the force is linearly proportional to the
spring extension and both are unlimited (here lis 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 [24]:
Marko-Siggia model
3Hl "|xij|
4 (1 − |xij |/l)2#xij
Cohen Pad´e model
1(|xij|/l)2xij (3.48)
FENE model (Warner spring law)
1(|xij|/l)2xij (3.49)
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´e models present sharper gradients of force than the
Marko-Siggia model, which directly impacts the numerical stability of the time
integration algorithm.
3.8.4 Time integration algorithm
The integration over time of Eq. (3.43) can be performed with explicit, semi-
implicit or implicit methods, which differ essentially in numerical stability and
CHAPTER 3. Theoretical background 32
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 [26]).
The explicit first-order Euler method evolves the beads positions from the
previous time-step (t) to the new time-step (t+ ∆t) as,
i+ ∆tBt
where Bt
irepresents the whole right hand side of Eq. (3.43) evaluated at the
previous time-step. This integration scheme only requires the explicit evaluation
of mathematical expressions, presenting a low computational cost per time-step.
However, the numerical stability of the scheme is highly dependent on ∆t, which
should be kept sufficiently small. The numerical stability is evaluated by the
capability of the method in respecting the constraint Ril, when bounded spring
force models are used. In practice, the time-step that satisfies such constraint is
relatively small, leading to the need of a very high number of iterations (time-
steps) to simulate a given period of physical time. Therefore, the explicit time
integration is seldom used.
In the semi-implicit scheme described in [4], the explicit Euler scheme (Eq.
3.50) is used while Ri< αl, where 0 < α 1 is defined by the user, and is usually
close to 1. Once this condition is violated, a two-steps computation is used:
i+ ∆t"uf+
kT +DiiFEV
kT +6
i+ ∆tDiiFs
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),
where fkis a vector function whose expression depends on the spring model, Jk
is the Jacobian of fkand ∆rkis a vector representing the difference in the beads
positions between the previous and the current iteration. Further details are given
in Ref. [4]. The linear system of equations (3.53) can be solved using different
Chapter 4
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.
4.1 The constitutiveEquations library
4.1.1 Available GNF and viscoelastic models
The constitutiveEquations library is a main component of rheoTool, since it con-
tains all the viscoelastic and GNF constitutive equations, which can be called from
the solvers. It was derived from the viscoelasticTransportModels library [1]. How-
ever, instead of restricting the library to viscoelastic models, we also extend it to
include GNF models, most of them already present in OpenFOAM R
. This was
done in order to allow accessing both classes of models from a single library, hence
from a single solver.
Most of the models available in the constitutiveEquations library are displayed
in Table 4.1, along with the respective expressions to be used in Eqs. (3.3), (3.4),
(3.6) and (3.8). However, some models falling in special categories as elastovis-
coplasticity and multispecies modeling, are presented in the text following the
table, in order to provide a more detailed discussion about them.
CHAPTER 4. Overview of rheoTool 34
Table 4.1: Available constitutive models in the constitutiveEquations library.
GNF models
Model 1TypeName ηs( ˙γ)
Newtonian Newtonian η
2(Bounded) Power-Law PowerLaw max ηmin,min ηmax, k ˙γn1
Carreau-Yasuda CarreauYasuda η+ (η0η)[1 + (k˙γ)a]n1
2Herschel-Bulkley HerschelBulkley Bounded: min η0, τ0˙γ1+k˙γn1
3Papanastasiou reg.: min η0, τ0˙γ1[1 exp(m˙γ)] + k˙γn1
2Casson Casson
Bounded: max ηmin,min ηmax,η+qτ0
Papanastasiou reg.: nη+qτ0
˙γ1exp m˙γo2
1Corresponds to the name entry identifying the model in the source code.
2Special care is taken in these models to avoid division by zero when ˙γis zero or very small and n1<0. For ˙γ < VSMALL, the
value ˙γ=VSMALL is used in the computation of the shear viscosity (VSMALL = 10300 for versions using double precision).
3The 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, η0should 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.
2, with ˙
Iis the identity tensor and D
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=1
CHAPTER 4. Overview of rheoTool 35
Continuation of Table 4.1
Viscoelastic models solved in the standard extra-stress or conformation tensor variables
Model TypeName ηs( ˙γ)ηp( ˙γ)λ( ˙γ) Constitutive Equation
Oldroyd-B Oldroyd-B ηsηpλτ+λ
(Carreau-Yasuda) WhiteMetznerCY ηsηp[1 + (K˙γ)a]n1
aλ[1 + (L˙γ)b]m1
bτ+λ( ˙γ)
τ=ηp( ˙γ)(u+uT)
Giesekus Giesekus ηsηpλτ+λ
ηp(τ·τ) = ηp(u+uT)
PTT linear PTTlinear ηsηpλh1 + ελ
PTT exponential PTTexp ηsηpλe
FENE-CR FENE-CR ηsηpλh1 + λD
where f=L2+λ
where f=L2+λ
L23and a=L2
3Rolie-Poly Rolie-Poly ηsηpλD
where k=3χ2
max 11
max 31
max and χ=qtr(A)
eXtended Pom-Pom XPomPom ηsηpλB
ηp(τ·τ) + ηp
λB(f1) I=ηp(u+uT)
where f= 2λB
q(Λ1) 11
Λn+1 +1
(ηPB)2iand Λ=q1 + tr(τ)
3See Ref. [27]. This model is exclusively solved in the conformation tensor variable, which is then converted to τusing, τ=ηp
CHAPTER 4. Overview of rheoTool 36
Continuation of Table 4.1
Viscoelastic models solved with the log-conformation approach
Model TypeName Θ¸τ4,5Constitutive Equation
6Oldroyd-B Oldroyd-BLog τ=ηp
(Carreau-Yasuda) WhiteMetznerCYLog τ=ηp
λ( ˙γ)eΘI
Giesekus GiesekusLog τ=ηp
PTT linear PTTlinearLog τ=ηp
λn1 + ε
PTT exponential PTTexpLog τ=ηp
λeΘI, where f=L2
λaeΘfI, where a=L2
L23and f=L2
8Rolie-Poly Rolie-PolyLog τ=ηp
eXtended Pom-Pom XPomPomLog τ=ηp
λBeΘ(f2α)eΘ+αeΘeΘ+ (α1)I
where f= 2λB
q(Λ1) 11
Λn+1 +1
3tr(eΘ(eΘ2I))and Λ=qtr(eΘ)
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.
4For the shortness of notation, we have introduced the operator: Υ=Θ
t +u·Θ(ΩΘ ΘΩ)2B.
5The following equivalences hold true: eΘ=A=RΛRTand eΘ=A1=1RT.
6For this model, we also included the square-root conformation approach [15] (TypeName: Oldroyd-BSqrt) and the rootkkernel approach [14] (TypeName: Oldroyd-
BRootk), for demonstration purposes.
7This log-conformation tensor approach of the White-Metzner model is only applicable when ηp( ˙γ)
λ( ˙γ)=ηp
λis constant, i.e., for K=L,a=band n=m. The version
based on the extra-stress tensor variable is more general and does not have this restriction.
8The expression for kis the same as for the model solved in variable A, in the previous page, considering that A=eΘ.
CHAPTER 4. Overview of rheoTool 37
In a footnote of Table 4.1, the (invariant) shear-rate used to compute shear-rate
dependent variables was defined as
2=2D:D,with ˙
γ=u+uTand D=1
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,
symm(fvc::grad(U())) =1
sqrt(2.0)*mag(symm(fvc::grad(U()))) =2q1
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. [15]. The second approach (TypeName: Oldroyd-BRootk) allows to apply a
general rootkkernel, as described in Ref. [14]. 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 [28] 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 [28],
t +u· ∇nA= 2DA2nA+1
t +u· ∇nB= 2DB2nBcBn2
+ 2cAnA
CHAPTER 4. Overview of rheoTool 38
where nis the dimensionless number density of the specie, λis the relaxation time,
Dis the diffusivity coefficient and cAand cBare, respectively, the dimensionless
breakage and reformation rates, expressed as
cA=cAEq +χ
cB=cBEq (4.5)
In Eqs. (4.2) and (4.3), the double contraction term originally presented in [28]
is not included, in order to simplify the definition of no-flux boundary conditions
for nAand nBat 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 [28],
2λADB2B=2cBnBB+ 2cAA(4.7)
where Aand Brepresent the conformation tensor of each species, and =λB
λA. The
contribution of each species to the polymeric extra-stress tensor is given by [28],
τ=G0[(A+ 2B)(nA+nB)I] (4.8)
where G0is the elastic modulus. In the absence of flow, nA= 1, nB=q2cAEq
2Iand τ=0.
+RRM (TypeName: RRM )
Also in the context of modeling the flow of wormlike micellar solutions, Dutta
and Graham proposed recently the Reactive Rod Model (RRM), accounting for
the formation/destruction of flow-induced structures [29]. In this model, micelles
are approached by rods, whose orientation tensor, S, evolves according to [29],
3+uT·S+S· ∇u2uT:<uuuu >(4.9)
L3ln L+m
is a time-varying diffusion coefficient (units are s1), L=L/L0is the time-varying
rod length normalized by its initial value and mis the initial aspect ratio of the
rods (see [29] for more details). In Eq. (4.9), the last term is approximated by [29]
uT:<uuuu >1
5[S·D+D·SS·S·DD·S·S+ 2S·D·S+ 3(S:D)S]
CHAPTER 4. Overview of rheoTool 39
and D=1
2u+uTis the rate of deformation tensor. The variation of the
normalized rod length is computed from [29]
α+β/P e 2(1 L) + kDr,0r3
where λS,α,βand kare parameters of the model, ˆ
3and P e is a P´eclet
number computed locally, P e =˙γ
Dr,0( ˙γis the strain-rate defined in Eq. 4.1). The
extra-stress due to the rods is accounted for as [29]
2DruT:<uuuu >(4.13)
where G0is the elastic modulus. Note that due to the generic strain-rate defi-
nition used in the P´eclet number, the strain-rate retrieved for a pure shear-flow
corresponds effectively to the local shear-rate, but for a pure extensional flow, the
strain-rate computed does not correspond to the extension rate (it differs by a con-
stant, that can be incorporated in parameter βof Eq. 4.12). Those considerations
are important when using rheoTestFoam. Other measures of the hydrodynamic
stresses can be used in the P´eclet number definition in order to generalize it for
any flow.
+Bautista-Manero-Puig model (TypeName: BMP and BMPLog)
The Bautista-Manero-Puig (BMP) model is still another option to simulate
worm-like micellar solutions [30]. This model can also predict thixotropy and it
results from the combination between the UCM model and Fredrickson’s kinetic
equation [30].
The polymeric extra-stress tensor is evolved according to
τ=G0(u+uT) (4.14)
and Fredrickson’s kinetic equation governing the structure parameter (ϕ) is repre-
sented by
t +u· ∇ϕ=ϕ0ϕ
In Eqs. (4.14) and (4.15), D=1
2(u+uT), G0is the instantaneous relaxation
modulus, ϕis the fluidity (η1
P), ϕ0is the zero shear-rate fluidity, ϕis the
infinite shear-rate fluidity, λis the structural relaxation time and kis a kinetic
constant for structure breaking down [30].
The model is also implemented within the log-conformation approach, taking
a form similar to the log-transformed Oldroyd-B equation presented in Table 4.1,
with λ1replaced by ϕG0and ηPreplaced by ϕ1.
+Saramito’s model (TypeName: Saramito and SaramitoLog)
Elastoviscoplastic fluids exhibit a solid-like behavior below the yield stress and
they flow as viscoelastic fluids when the yield stress is exceeded. The stress-strain
relation in each regime can assume several forms.
CHAPTER 4. Overview of rheoTool 40
The model proposed by Saramito [31] attempts to merge the Herschel–Bulkley
model for yield-stress fluids with the Oldroyd-B/PTT model for viscoelastic fluids.
Accordingly, the resulting constitutive relation is given by
f(τ)ηpmax 0,στ0
τ=ηp(u+uT) (4.16)
where τ0is the yield stress, f(τ) = 1 and σ=IIτD=pτD:τD
2is the second
invariant of the deviatoric stress tensor, τD=τtr(τ)
NI. An elastic modulus can
be defined, G=ηp
λ, which is often found in the literature in the description of
this model. We remember that
τrepresents the Gordon-Schowalter derivative
defined as,
τ+ζ(τ·D+D·τ), with D=1
2(u+uT), which reduces to the
upper-convected derivative for ζ= 0.
For n= 1 and k=ηp, Saramito’s model degenerates into a previous model pro-
posed by the same author [32], that merges the Bingham model with the Oldroyd-B
or PTT models, depending on the form taken by f(τ):
f(τ) =
1 , Oldroyd-B
1 + ελ
ηptr(τ) , linear PTT
ηptr(τ), exponential PTT
The four variants of the model are available in rheoTool, where both the Her-
schel–Bulkley and PTT variants can avoid the possible infinite Oldroyd-B elonga-
tional viscosity for W i 0.5 in extensional flows of Oldroyd-B fluids. In addition,
the four variants can also be solved with the log-conformation approach and the
governing equation becomes
t +u·Θ(ΩΘ ΘΩ)2B=g(τ)
g(τ) =
ηpmax 0,στ0
n, Oldroyd-B–Herschel-Bulkley
max 0,στ0
σ, Oldroyd-B–Bingham
max 0,στ0
σn1 + ε
1ζtr(eΘ)3o, lin. PTT–Bingham
max 0,στ0
1ζ(tr(eΘ)3) , exp. PTT–Bingham
and the polymeric extra-stress tensor is recovered using τ=ηp
+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 [33]. According to this model, the equation
governing the polymeric extra-stress tensor is given by [33]
max 0,1λky
στeff +λλE
τ=ληp(u+uT) (4.20)
CHAPTER 4. Overview of rheoTool 41
where λEis the viscoelastic relaxation time, τeff =τκback, with κback =
khA+ 2A2, and σ=IIτD
eff =qτD
2, where τD
eff =τeff tr(τeff)
NI(Nis the
number of dimensions of the problem). Eq. (4.20) closely resembles the governing
equation for τin Saramito’s model (Eq. 4.16). Tensor Ais related with material
hardening and its evolution follows [33]
t +u· ∇A+·AA·=Dp·A+A·Dp+DpqdpA(4.21)
where =u− ∇uT/2 is the vorticity tensor, dp=qDp:Dp
Dp=(0, if σ < λky
σ, 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 [33],
where each of the Mmodes of λobeys
dt =Dik1φaλn
i+k2φb(1 λi) + k3(1 λi)(4.24)
and φcan be computed in two ways [33],
φ=(max (0, σ λky) , Stress-controlled form
2dp, Rate-controlled form (4.25)
The model implementation in rheoTool allows for arbitrary Mmodes, where
lists Cand D(size M) should be provided as input by the user. Thus, the complete
list of input parameters for this model is: ηp, C, D, a, n, b, k1, k2, k3, ky, kh, q and λE,
to which we should add the solvent viscosity (ηs) and the fluid density (ρ). If no
initial conditions are provided for each λi, the solvers assume λi,0= 1 by default.
As discussed in [33], the transformation of the ML-IKH model from its scalar
version to the above tensorial version assumes some simplifications, which impose
restrictions on the model applicability (see [33] for more details).
4.1.2 A note on FENE-type models
The Finitely Extensible Non-linear Elastic (FENE) models were originally devel-
oped based on the representation of polymer molecules by elastic dumbbells [7].
In such analysis, the end-to-end vector for each molecule is naturally related with
the conformation tensor, such that the constitutive equations for this family of
models is frequently written and handled as a function of the conformation tensor.
The polymeric contribution to the momentum equation is then accounted for by
CHAPTER 4. Overview of rheoTool 42
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 [34]. For the FENE-CR and FENE-P models, the
complete constitutive equation written as a function of Aand τand the modified
formulation in τare:
Complete in A:
A=f(tr(A))(AI), where f(tr(A)) = L2
Complete in τ(see Table 4.1):
h1 + λD
τ=ηp(u+uT), where f=L2+λ
Modified in τ(usually known as FENE-MCR):
τ=ηp(u+uT), where f=L2+λ
Complete in A:
A=[f(tr(A))AaI], where f(tr(A)) = L2
L2tr(A)and a=L2
Complete in τ(see Table 4.1):
f[λτ+pI], where f=L2+λ
and a=L2
Modified in τ:
f(u+uT), where f=L2+λ
L23and a=L2
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 Aor τ, should provide
similar results, since they are mathematically equivalent. Due to discretization
CHAPTER 4. Overview of rheoTool 43
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 Afor FENE-type models. Those are the most sta-
ble, the most accurate regarding the original theory presented in [7] 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.
4.1.3 Multi-mode modeling
Similarly to the viscoelasticTransportModels library [1], the constitutiveEquations
library also supports multi-mode modeling for viscoelastic models. In such cases,
the total extra-stress tensor is the sum of the extra-stress tensor resulting from
each kth mode
k=1 τk+τk
In practice, this is achieved by assembling and solving one constitutive equa-
tion for each kth mode, that is, Eq. (3.3) – solvent contribution – and Eq. (3.4)
or (3.8) – polymer contribution – are built and solved Ntimes each time-step. A
warning should be made at this point, since this approach is probably not the most
conventional. Indeed, τsin Eq. (4.26) is commonly placed outside the summation
symbol, since multiple modes are only assigned to the polymeric contribution. To
achieve this in rheoTool, and considering the expression for τsin Eq. (3.3), the
user must split the ”single-solvent viscosity” by the Nmodes consid-
ered, in any way, such that this ”single-solvent viscosity” is recovered summing
all these Nvalues in Eq. (4.26).
4.1.4 Analysis of a code sample
For the readers still initiating their journey in OpenFOAM R
, we will explore in this
section the implementation of the Oldroyd-B constitutive model, solved with the
log-conformation tensor approach. This example will establish the link between
part of the theory described in Chapter 3and its implementation in the source
The source code displayed in Listing 4.1 is taken from file src/libs/constitut
Let’s analyze the most important lines:
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 44
starting a simulation, except eigVals and eigVecs , which can be defined
or not. If defined (typical of a restart from a previous simulation), they are
used in the first time-step; otherwise, they are both initialized as the identity
tensor/matrix, corresponding to a null extra-stress tensor (τ). Afterwards,
the fluid properties are read from a dictionary (lines 85-88), along with the
stabilization method selected by the user (line 90): none, BSD or the stress-
velocity coupling described in Section 3.3.2.
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 [35]. 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
in the viscoelasticTransportModels library [1] each model was in charge to
define its own contribution to the momentum equation, i.e., the term (∇·τ0).
In the constitutiveEquations library there is a default definition of this term
in the base class. In fact, the function divTau() is now defined in class
constitutiveEq and can be found in file constitutiveEq.C, Listing 4.2. This
function starts by distinguishing between GNF and viscoelastic models in
line 7. For a GNF model (lines 8-14), the extra-stress contribution is ∇·τ0=
∇·η( ˙γ)u+u·η( ˙γ), divided by the density to be compliant with the
usual strategy of OpenFOAM R
for single-phase, incompressible fluid flows.
Note that the second term is included to account for a shear-rate dependent
viscosity coefficient. By definition, a GNF fluid has no elasticity, thus τ=0.
For a viscoelastic fluid (lines 17-49), the output depends on the stabilization
method selected (see function checkForStab() in constitutiveEq.C for the
correspondence between indexes and the method): if none, there is no added
stabilization and ∇·τ=∇·τ−∇·ηsu; if BSD, then the both-sides-diffusion
technique is used and ∇·τ=∇·τ− ∇·ηpu+∇·(ηs+ηp)u(Eq. 3.5);
otherwise (if coupling), the stress-velocity coupling technique of Eq. (3.17)
is used and ∇·τ=∇·τ− ∇·ηpu+∇·(ηs+ηp)u. Note that using the
coupling stabilization is the method recommended for most of the cases,
CHAPTER 4. Overview of rheoTool 45
being the one used by default if no information is provided by the user.
However, some cases may require the use of no stabilization, as for example
the simulation of multimode models with ηPηS– the amount of artificial
diffusion may mask the real phenomena in transient simulations. For the
cases using stabilization, the explicit behavior effects on transient results
can be minimized by performing inner iterations at each time-step, a subject
discussed later in this guide (see Section 4.5.1). In file consitutiveEq.C, a
function divTauS() is also included, which retrieves part of the extra-stress
contribution to the momentum equation, when solving two-phase flows (this
topic will be discussed later).
1#include "Oldroyd_BLog.H"
#include "addToRunTimeSelectionTable.H"
// **************Static Data Members ********
namespace Foam{
7namespace constitutiveEqs{
defineTypeNameAndDebug(Oldroyd_BLog, 0);
9addToRunTimeSelectionTable(constitutiveEq, Oldroyd_BLog,
11 }
// ****************Constructors *********
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,
29 U.mesh(),
31 IOobject::AUTO_WRITE
33 U.mesh()
35 theta_
37 IOobject
39 "theta" + name,
CHAPTER 4. Overview of rheoTool 46
41 U.mesh(),
43 IOobject::AUTO_WRITE
45 U.mesh()
47 eigVals_
49 IOobject
51 "eigVals" + name,
53 U.mesh(),
55 IOobject::AUTO_WRITE
57 U.mesh(),
59 (
61 dimless,
63 ),
65 ),
67 (
69 (
"eigVecs" + name,
71 U.time().timeName(),
75 ),
77 dimensionedTensor
79 "I",
81 pTraits<tensor>::I
83 zeroGradientFvPatchField<tensor>::typeName
85 rho_(dict.lookup("rho")),
87 etaP_(dict.lookup("etaP")),
89 {
91 }
93 // ***************Member Functions ********
void Foam::constitutiveEqs::Oldroyd_BLog::correct()
95 {
CHAPTER 4. Overview of rheoTool 47
// Decompose grad(U).T()
volTensorField L = fvc::grad(U());
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_);
105 decomposeGradU(M, eigVals_, eigVecs_, omega, B);
107 // Solve the constitutive Eq in theta = log(c)
109 dimensionedTensor Itensor
111 "Identity",
dimensionSet(0, 0, 0, 0, 0, 0, 0),
113 tensor::I
fvSymmTensorMatrix thetaEqn
117 (
119 + fvm::div(phi(), theta_)
121 symm
123 (omega&theta_)
- (theta_&omega)
125 + 2.0 *B
+ (1.0/lambda_)
127 *(
eigVecs_ &
129 (
131 - Itensor
133 & eigVecs_.T()
135 )
137 );
139 thetaEqn.relax();
// Diagonalization of theta
calcEig(theta_, eigVals_, eigVecs_);
// Convert from theta to tau
tau_ = (etaP_/lambda_) *symm( (eigVecs_ & eigVals_ & eigVecs_
.T()) - Itensor);
CHAPTER 4. Overview of rheoTool 48
151 }
Listing 4.1: Source code for the Oldroyd-BLog constitutive model
(Oldroyd BLog.C)
1tmp<fvVectorMatrix> constitutiveEq::divTau
3const volVectorField& U
7if (isGNF())
11 fvm::laplacian( eta()/rho(), U, "laplacian(eta,U)")
+ (fvc::grad(U) & fvc::grad(eta()/rho()))
13 );
15 else
17 if (stabMeth_ == 0) // none
21 (
fvc::div(tau()/rho(), "div(tau)")
23 + fvm::laplacian(etaS()/rho(), U, "laplacian(eta,U
27 else if (stabMeth_ == 1) // BSD
31 (
fvc::div(tau()/rho(), "div(tau)")
33 - fvc::laplacian(etaP()/rho(), U, "laplacian(eta,U
+ fvm::laplacian( (etaP()+ etaS())/rho(), 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, "
CHAPTER 4. Overview of rheoTool 49
49 }
Listing 4.2: Source code of the virtual function divTau() defined in file
4.1.5 Advanced settings
As aforementioned, the eigenvectors/eigenvalues used in the models based on the
log-conformation approach are computed, by default, using the QR algorithm pro-
vided by the Eigen library [35]. However, there is also the possibility to use an
iterative Jacobi method [36]. While both options offer good accuracy and stabil-
ity, the QR algorithm was seen to be slightly faster and this is the reason of being
the default option. Switching between either methods is not run time selectable,
but hard-coded, instead. This can be controlled in member function calcEig()
of class constitutiveEq, located in file constitutiveEq.C. The Jacobi method
can be selected by uncommenting the currently commented block (// Eigen de-
composition using the iterative jacobi algorithm) and commenting the block (//
Eigen decomposition using a QR algorithm of Eigen library), i.e., all the remain-
ing lines inside the function. The source code of jacobi() function is located in
utils/jacobi.H. Note that, independently of which method is used in func-
tion calcEig(), this function will return the eigenvectors of the input tensor, and
the exponential of the eigenvalues of the same tensor. After re-compiling the
library with those modifications, all the log-conformation-based models will be
affected by those changes.
4.1.6 Adding new viscoelastic or GNF models
For the users with minimal programming skills in OpenFOAM R
, adding new vis-
coelastic or GNF models should be a straightforward task. The main steps are:
copy-paste the folder of an existing model (that we will call template
model) in directory src/libs/constitutiveEquations/constitut
iveEqs/ for viscoelastic models, or src/libs/constitutiveEquatio
ns/constitutiveEqs/GNF/ for GNF 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, as well as the folder for the log-implementation (if
not needed), in case of viscoelastic models.
inside the source .C and header .H files, find-replace the old model’s name
by the new one (e.g. Ctrl +Hin 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 50
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
add the header #include ”constitutiveModel.H” to the main solver.
create a constitutiveModel object by calling the constructor, with the correct
arguments, for example: constitutiveModel constEq(U, phi).
add the library constitutiveEquations to the Make/options, and specify its
path (check the Make/options of the solvers in rheoTool for an example).
4.2 The EDFModels library
4.2.1 Available EDF models
The list of runtime selectable EDF models is presented in Table 4.2.
Table 4.2: Available models for electrically-driven flows in the EDFModels library. The last column indicates the section in the user guide where the model
has been described.
Model 1TypeName 2,3,4fE5Governing Equations Section
6NernstPlanckCoupled F
∇·(εφExt) = 0
∇·(εψ) = F
t +u·ci=∇·(Dici) + ∇·Diezi
kT Ψci
Section 3.7.1
(PB) PoissonBoltzmann F
zici,0exp ezi
kT ψ(ΨEa)
∇·(εφExt) = 0
∇·(εψ) + ψF
with bi=ezi
kT and ai=zici,0exp (biψ)
Section 3.7.3
(DH) DebyeHuckel F
kT ψ(ΨEa)
∇·(εφExt) = 0
∇·(εψ) = F
kT ψSection 3.7.4
Slip velocity slipSmoluchowski 0∇·(εφExt) = 0 Section 3.7.5
Ohmic Ohmic [∇·(εφExt)] (φExt Ea)
t +u·σ=Deff2σ
∇·(σφExt) = 0
Section 3.7.6
1Corresponds to the name entry identifying the model in the source code.
2fEis the electric body-force entering the momentum equation.
3Eais an optional argument - a single vector - representing a uniform electric field.
4When the splitting of potentials approach is selected for the PNP, PB and DH models, the user may choose to use (φExt Ea), instead of (ΨEa). Recall that the splitting of
potentials is given by Ψ=φExt +ψ(cf. Section 3.7.2).
5For the PNP, PB and DH models, the equations are presented according to the splitting of potentials approach, which is optional. When a single electric potential is intended to be
used, then replace ψby Ψand ignore the equations in terms of φExt.
6Both models solve the PNP equations, but NernstPlanck uses a segregated solver, whereas NernstPlanckCoupled solves the PNP equations in a fully-coupled way, which can be
optionally coupled to the Navier-Stokes equations.
CHAPTER 4. Overview of rheoTool 52
4.2.2 The potentials splitting approach and multi-species
modeling in the PNP, PB and DH models
The possibility of splitting the electric potential into two variables, as described
in Section 3.7.2, is available for the PNP, PB and DH models. When used, one
Poisson equation for ψand one Laplace equation for φExt are solved, as shown in
Table 4.2. In practice, the choice between using one or two potentials is achieved
in the following way: if only one potential ( psi Ψ) is present in the starting-
time folder, then it is assumed than a unique potential is to be used, while if two
potentials (phiE φExt and psi ψ) are defined, then the splitting approach is
assumed. Under the splitting approach, the user still has the option to include or
not the contribution of the intrinsic potential in the electric field definition of the
body-force entering the momentum equation. The choice is through the variable
psiContrib, which should be defined in a dictionary, as explained later in Section
All the PNP, PB and DH models support multi-species modeling, with an ar-
bitrary number of species, each having different properties (charge valence, and
diffusivity, when it applies). On the other hand, the Ohmic model is only imple-
mented for a binary, symmetric electrolyte, although the two species may have
different diffusion coefficients.
4.2.3 Electrokinetic coupling loop in the PNP model
The PNP model has a loop for the coupling between the Nernst-Planck equations
(ionic concentration) and the Poisson equation (electric potential). This loop was
seen to be required to keep the second-order accuracy in time of the PNP model
[3]. Furthermore, it also allows the use of higher time-steps, while ensuring the
conservation of ions. Although it is allowed to select one single iteration for this
loop, we recommend the use of at least two iterations in any generic case. Moreover,
two coupling iterations is the default behavior for this model if no information is
provided by the user.
4.2.4 Coupled PNP model
This feature is only available for OpenFOAM R
The PNP system of equations can be solved coupled instead of segregated. The
segregated solution method is implemented in the NernstPlanck model, whereas
the coupled solution method is implemented in the NernstPlanckCoupled model.
The coupled solver, which has been presented in [5], displays a higher numerical
stability, but consumes more computational resources (memory and CPU time) per
time-step. The PNP system of equations can not only be solved coupled alone, as
it can also be coupled with momentum and continuity equations.
The technical aspects related with the use of the NernstPlanckCoupled model
are discussed in Section 4.4.8.
CHAPTER 4. Overview of rheoTool 53
4.2.5 Analysis of a code sample
As previously done in Section 4.1.4 for the constitutiveEquations library, in this
Section we analyze a piece of code from the EDFModels library in order to illustrate
the connection between the code and the theory previously discussed. However,
we should note that the differences, at the code level, between the different models
of the EDFModels library are bigger than in the constitutiveEquations library, as
can be deduced from Table 4.2.
The piece of code selected for that purpose, Listing 4.3, represents the imple-
mentation of the PNP model and can be found in src/libs/EDFModels/mod
els/NernstPlanck/NernstPlanck.C. Let’s start the analysis to the most
important parts:
lines 11-33: this is the constructor of a subclass (NPSpecie), nested in the
main class (NernstPlanck). Remember that the PNP, PB and DH models
were all presented in a multi-species formulation, which is the form avail-
able in the code. The NPSpecie subclass is exactly each specie of the PNP
model. For each new specie, we can see the initialization of the following
attributes (members): the concentration field (ci), the charge valence (zi)
and the diffusivity (Di). The NernstPlanck class may have Ninstances of
the NPSpecie subclass, as much as defined by the user.
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.7.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
CHAPTER 4. Overview of rheoTool 54
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 55
#include "NernstPlanck.H"
2#include "addToRunTimeSelectionTable.H"
4// *********Static Data Members ********//
namespace Foam{
6namespace EDFEquations{
defineTypeNameAndDebug(NernstPlanck, 0);
8addToRunTimeSelectionTable(EDFEquation, NernstPlanck, dictionary);
10 }
// **********Constructors **********//
12 Foam::EDFEquations::NernstPlanck::NPSpecie::NPSpecie
14 const word& name,
const surfaceScalarField& phi,
16 const dictionary& dict
18 :
20 (
22 (
24 phi.time().timeName(),
26 IOobject::MUST_READ,
28 ),
30 ),
32 Di_(dict.lookup("D"))
// **********Constructors **********//
36 Foam::EDFEquations::NernstPlanck::NernstPlanck
38 const word& name,
const surfaceScalarField& phi,
40 const dictionary& dict
42 :
EDFEquation(name, phi),
44 solvePhiE_(checkForPhiE(name, phi)),
46 (
48 (
"psi" + name,
50 phi.time().timeName(),
52 IOobject::MUST_READ,
54 ),
56 ),
CHAPTER 4. Overview of rheoTool 56
58 (
60 (
"phiE" + name,
62 phi.time().timeName(),
solvePhiE_ == false ? (IOobject::NO_WRITE) : (IOobject::
66 ),
68 dimensionedScalar
70 "zero",
72 pTraits<scalar>::zero
74 zeroGradientFvPatchField<scalar>::typeName
76 relPerm_(dict.lookup("relPerm")),
78 extraE_(dict.lookupOrDefault<dimensionedVector>("extraEField",
dimensionedVector("0", dimensionSet(1, 1, -3, 0, 0, -1, 0),
psiContrib_(dict.lookupOrDefault<bool>("psiContrib", true)),
80 phiEEqRes_(phi.mesh().solutionDict().subDict("electricControls").
subDict("phiEEqn").lookupOrDefault<scalar>("residuals", 1e-7)),
subDict("psiEqn").lookupOrDefault<scalar>("residuals", 1e-7)),
82 ciEqRes_(phi.mesh().solutionDict().subDict("electricControls").
subDict("ciEqn").lookupOrDefault<scalar>("residuals", 1e-7)),
subDict("phiEEqn").lookupOrDefault<scalar>("maxIter", 50)),
84 maxIterPsi_(phi.mesh().solutionDict().subDict("electricControls").
subDict("psiEqn").lookupOrDefault<scalar>("maxIter", 50)),
subDict("ciEqn").lookupOrDefault<scalar>("maxIter", 50)),
86 nIterPNP_(phi.mesh().solutionDict().subDict("electricControls").
lookupOrDefault<int>("nIterPNP", 2)),
88 nSpecies_(0)
90 PtrList<entry> specEntries(dict.lookup("species"));
nSpecies_ = specEntries.size();
92 species_.setSize(nSpecies_);
94 forAll (species_, specI)
96 species_.set
98 specI,
new NPSpecie
100 (
102 phi,
104 )
CHAPTER 4. Overview of rheoTool 57
106 }
// ******** Member Functions * *********//
110 Foam::tmp<Foam::volVectorField> Foam::EDFEquations::NernstPlanck::Fe()
112 volScalarField rhoE( psi_ *dimensionedScalar("norm", epsilonK_.
dimensions()/dimArea, 0.) );
114 forAll (species_, i)
116 rhoE
+= (
118 species_[i].zi()*species_[i].ci()*FK_
120 }
122 if (solvePhiE_)
124 if (psiContrib_)
126 return
128 -rhoE *( fvc::grad(phiE_+psi_) - extraE_)
130 }
132 {
134 (
-rhoE *( fvc::grad(phiE_) - extraE_)
136 );
138 }
140 {
if (psiContrib_)
142 {
144 (
-rhoE *( fvc::grad(psi_) - extraE_)
146 );
148 else
150 return
152 -rhoE *(-extraE_)
154 }
156 }
158 void Foam::EDFEquations::NernstPlanck::correct()
CHAPTER 4. Overview of rheoTool 58
// Electrokinetic coupling loop
162 for (int j=0; j<nIterPNP_; j++)
Info << "PNP Coupling iteration: " << j << endl;
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_)
182 res=phiEEqn.solve().initialResidual();
184 iter++;
186 }
188 //- Equation for the intrinsic potential
190 res=GREAT;
volScalarField souE(psi_ *dimensionedScalar("norm1",dimless/dimArea
forAll (species_, i)
196 {
souE +=
198 (
200 /(relPerm_*epsilonK_)
202 }
204 while (res > psiEqRes_ && iter < maxIterPsi_)
fvScalarMatrix psiEqn
208 (
210 ==
212 );
214 psiEqn.relax();
CHAPTER 4. Overview of rheoTool 59
218 }
220 //- Nernst-Planck equation for each ionic specie
222 forAll (species_, i)
224 res=GREAT;
volScalarField& ci = species_[i].ci();
dimensionedScalar cf(species_[i].Di() *eK_ *species_[i].zi() / (kbK_
*T_) );
while (res > ciEqRes_ && iter < maxIterCi_)
232 {
234 fvScalarMatrix ciEqn
236 fvm::ddt(ci)
+fvm::div(phi(), ci, "div(phi,ci)")
238 ==
fvm::laplacian(species_[i].Di(), ci, "laplacian(D,ci)")
240 +fvc::laplacian(ci*cf, phiE_+psi_, "laplacian(elecM)")
244 res=ciEqn.solve(phi().mesh().solver("ci")).initialResidual();
246 ci = Foam::max( dimensionedScalar("lowerLimit",ci.dimensions(), 0.),
ci );
250 }
252 }
Listing 4.3: Source code of the Poisson-Nernst-Planck model in file
CHAPTER 4. Overview of rheoTool 60
4.2.6 Adding new EDF models
The steps required to add new EDF models are similar to the ones described
previously, in Section 4.1.6, for GNF and viscoelastic models. The main steps are:
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 +Hin 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:
add the header #include ”EDFModel.H” to the main solver.
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 61
4.3 The BDmolecule library
This library is only available for OpenFOAM R
Library BDmolecule contains the classes and routines implementing the Brownian
dynamics algorithm. The library is based on the solidParticle library provided
by OpenFOAM R
to perform the tracking of rigid particles. This library has been
copied and modified in order to allow the creation of molecules, which are simply
organized ensembles of beads. The interface of library BDmolecule is embodied
by class sPCloudInterface. The source code of the library can be found in src/
4.3.1 Organization of variables
The base class particle is commonly used in OpenFOAM R
to create a single La-
grangian entity that can be tracked. Several other classes are derived from this
one, adding new features to it. A Cloud is a template class representing a set
of particles, inheriting the properties of C++ doubly-linked lists. This type of
lists is not the most versatile one when we need to perform operations between
non-consecutive elements, as it happens when computing intra-molecular forces.
Thus, directly creating a molecule from a Cloud class seemed not to be the best
option, even more because this option would create 4Mfiles (Mbeing the number
of molecules) for each time-step saved, which would overburden any file trans-
fer operation, notwithstanding the small disk space used by each file. Therefore,
we decided to separate the particle tracking from the molecule-related tasks. As
such, we created a single (blind) Cloud object to contain all the particles (beads)
from all the molecules of the simulation and perform the tracking, and an addi-
tional structure establishing the link between each particle and the corresponding
molecule, that takes care of all the molecule-related tasks. We classified object
Cloud as blind, because it does not differentiate among the particles, ignoring the
molecules and groups to which they belong. This solves the two previous issues
related with indexing and the number of files generated, but requires the exchange
of information between structures and, consequently, some internal duplication of
data, which is not problematic in terms of performance, but hinders the efficient
parallelization of the code.
In the interface class (sPCloudInterface), object spc (a Cloud) is the one ded-
icated to the particle tracking. In addition, we can find a number of P trList <>
whose name starts with mand ends with an underscore (e.g. mx ,mSigma ,
mU , ...) which are used in the remaining operations. For example, mx is a
P trList < F ield < vector >> holding the beads positions (Cartesian coordi-
nates) for each molecule. The beads positions are also a member of spc , but the
difference comparing to mx is in the data organization. While we can use mx to
easily find the position of bead i, in molecule mof group g, this cannot be done
from inside spc .
CHAPTER 4. Overview of rheoTool 62
4.3.2 Solution sequence
For a proper understanding of the source code behind the Brownian dynamics
module of rheoTool, the user first needs to understand the interplay between the
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
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.8.3).
lines 5-6: if the user selects to account for hydrodynamic interactions, tensor
Dis computed from the RPY model (Eq. 3.44), and its Cholesky decompo-
sition results in tensor σ(Section 3.8.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
Dand σ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
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
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.
CHAPTER 4. Overview of rheoTool 63
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).
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
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.8.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.8.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.
CHAPTER 4. Overview of rheoTool 64
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.5.5, this will force rheoBDFoam to abort the run, even if
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()
mx0_ = mx_;
if (isHI_)
10 if (isExclusionVolumeF_)
// Move particles: Drag + Brownian + Exclusion Volume
16 moveAndReceive(true);
18 // Compute Spring force term explicitly and update local mU and
mxStar_ = mx_;
20 spModel_->fSpring();
22 // Check for violations in spring lengths and add spring force
implicitly if needed
spModel_->checkSpringsLength(mxStar_, mx0_);
// Push back the molecules if not tethered and if pushback
vector is not negligible
26 if (!isTethered_ && mag(pBackV_)>SMALL)
// Move particles: spring force only
32 moveAndReceive(false);
34 // Write data (controlled by output time)
// Write statistics (has its own control for output)
38 if (writeStats_)
CHAPTER 4. Overview of rheoTool 65
if (nMolc_>0)
42 {
return true;
44 }
46 {
return false;
48 }
Listing 4.4: Member function update() of class sPCloudInterface. The source
code can be found in file sPCloudInterface.C.
4.3.3 External forcing type
The external forcing (drag) acting on the beads is embodied by term ufin Eq.
(3.43). The forcing can be defined analytically or computed numerically.
In case it is defined analytically, tensor u, provided by the user, defines the
gradient of the forcing (spatially homogeneous and constant over time, by default).
The computational mesh in those cases can (should) be simply a single cell. The
computational domain can be made large in all the directions if an unconfined flow
is intended, but a confined flow can also be imposed by shortening one or multiple
directions. Importantly, the boundaries used to impose confined flow conditions
must be of wall type (patch types will be crossed by the molecules, and remaining
types give undefined behavior).
In case it is computed numerically, then two types are still available in rheoTool:
hydrodynamic forcing, in which case ufis the fluid velocity; electric forcing for
polyelectrolytes immersed in an electric field, in which case uf=µEis the elec-
trophoretic velocity (µis the electrophoretic mobility and Eis the electric field).
Both types of forcing can coexist, and each one requires that the respective con-
tinuum field is available in the case directory.
Note that 2D meshes, and the corresponding continuum fields, are allowed in
Brownian dynamics simulations, but the motion equation of the beads (Eq. 3.43)
is always solved for the 3 Cartesian coordinates. This means, for example, that
the molecules will not feel the z-component effect of a planar velocity field solved
in the Oxy plane, but they still feel the Brownian force (and the remaining, except
drag) in that direction.
4.3.4 External forcing interpolation
Independently of the forcing nature, ufmust be interpolated to the current position
of the bead. Therefore, the question addressed in this section is how to compute
ufgiven a numerical field defined at cell centers (numerical forcing), or a forcing
gradient valid over all the domain (analytical forcing)? The following methods are
available in rheoTool:
CHAPTER 4. Overview of rheoTool 66
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 uis user-defined.
BarycentricWeights: this is the method that OpenFOAMR
uses by default
to compute the velocity of Lagrangian particles. Consider a generic cell
and its decomposition in tetrahedrons, where one vertex of the tetrahedron
is always the center of the cell, and the remaining vertices correspond to
vertices of one of the cell’s faces. For a rectangular cell, for example, one
would end-up with 12 non-unique tetrahedrons, two per face of the cell. It
is then possible to define a barycentric coordinate system for each tetrahe-
dron, such that any point inside it can be uniquely identified by a quadruplet
(b1b2b3b4), with
bi= 1. In addition to allow the spatial location of par-
ticles, these coordinates can also be used to weight the data at the vertexes
of the tetrahedron, acting as interpolants of the data at the tetrahedron’s
vertices, which is the method adopted when selecting BarycentricWeights.
Remember that the vertices of the tetrahedrons are either cell centers or
vertices of cells. Therefore, the numerical field computed at the cell-centers
uses co-located grids) also needs to be interpolated to ver-
tices, which is accomplished by simple inverse distance weighting.
Gradient: in this approach, both the cell-centered field and its gradient are
used for the interpolation. Assuming that uCis the forcing known (com-
puted) at the cell center (located at xC) and that uCis its gradient (com-
puted numerically), then the velocity at any point xPinside the cell can be
approximated by: uf=uC+ (xPxC)·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.
Although the BarycentricWeights and Gradient methods display sub-cell reso-
lution, they both have some pitfalls, as illustrated in Fig. 4.1. The data/grid in this
figure corresponds to a fully-developed field in the x-direction, which displays a
maximum in the y-direction, for y= 0, and is symmetric about y= 0. The profile
in the y-direction can be a parabola, a piecewise linear function or any other func-
tion satisfying the aforementioned conditions of maximum and symmetry. This
could be the case, for example, of the velocity in a fully-developed Poiseuille flow
sampled at the centerline. The behavior of the two interpolation methods for this
case are plotted in the right of Fig. 4.1. For the case depicted, the interpolated
ufby the BarycentricWeights method is simply a field of quadrangular pyramids
(apex at cell centers), with the height of each vertex corresponding to the local
forcing value.
For a profile taken over y= 0, method Gradient approaches correctly the
theoretical function (constant over the x-direction), whereas method Barycen-
tricWeights retrieves a sawtooth profile. Consider now that the molecules’ center
CHAPTER 4. Overview of rheoTool 67
of mass is artificially fixed (see Section 4.3.6) somewhere between x=1 and
x= 1 (keeping y= 0). For the Gradient method, the results will be independent
of the specific abscissa selected, as one would expected for a fully-developed field
in the x-direction. On the other hand, if the strain-rate in the sawtooth profile
retrieved by method BarycentricWeights is such that W i =λ˙γ >> 1, then the
results will strongly depend on the specific abscissa selected.
For a profile taken over x= 0, method Gradient is unable to capture any varia-
tion of the field in that direction for the central cell, since it predicts uC= (0 0 0)
for all the cells at the centerline, which is only true exactly at the cell’s center. The
issue arises from the gradient computation method, and is akin to the checkerboard
problem that can happen for the U-pcoupling in the momentum equation. On
the other hand, method BarycentricWeights performs better this time, retrieving
a linear variation of ufover the y-direction (the interpolation is exact if the origi-
nal field displays a linear variation). If we consider a group of molecules traveling
over the x-axis, and only spanning the central layer of cells, then we can easily
conclude that they would not feel any forcing under the Gradient method, whereas
method BarycentricWeights would be able to approach the forcing gradient in the
y-direction (in addition to the artificial one in the x-direction).
2 2
2 2
y =1
y =3
y =-1
y =-3
x =-3 x =-1 x =1x =3
Profile at
y = 0
Profile at
x = 0
Figure 4.1: Hypothetical 2D external forcing field (a scalar field to simplify) in a
co-located grid (values at cell centers), aiming to represent a fully-developed field
in the x-direction. The green values defined at vertices are obtained by inverse
distance weighting of the cell-centered data. On the right, the profiles of uftaken
at x= 0 and y= 0 are represented for two different interpolation methods.
CHAPTER 4. Overview of rheoTool 68
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
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.
4.3.5 Spring force and time-integration schemes
The spring force models available in library BDmolecule are presented in Table 4.3
and were briefly discussed in Section 3.8.3.
Table 4.3: Available spring force models for Brownian dynamics simulations.
Model (1)
TypeName (2)FS
ij, a (j6=i, jgi)
Hookean Hookean H
xij – –
Marko-Siggia MarkoSiggia 2
jgih|xij |
4(1−|xij |/l)2ixij
|xij |rk
kT FS,k
dh 1
Cohen Pad´e CohenPade H
jgih3(|xij |/l)2
1(|xij |/l)2ixij rk
kT FS,k
(d21) hd23
1(|xij |/l)2xij rk
kT FS,k
(1) Corresponds to the name entry identifying the model in the source code.
(2) H=3kT Nk,s
(3) Vectorial function used in the Newton-Raphson method (see Section 3.8.4). Superscript kdenotes the iteration index.
(4) The expressions in this column are for the off-diagonal ij (j6=i, j gi) elements of the Jacobian matrix of the Cartesian component a=x, y, z of fk. All the
off-diagonal elements for which j /giare zero. The diagonal elements are obtained from the sum of the off-diagonal elements, Jk
ii,a = 1
im,a. This
symmetric matrix is used in the Newton-Raphson method (see Section 3.8.4).
(5) H0=HtD
kT ,d=|xk
ij |/l =|rk
i|/l,δ= (rk
j,a rk
CHAPTER 4. Overview of rheoTool 70
In general, the FENE and Cohen Pad´e 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.8.4 is also available. The auxiliary
functions used in the semi-implicit scheme are specified in Table 4.3 for all bounded
As mentioned in Section 3.8.4, several methods can be used to solve the lin-
ear system of equations (3.53). Four (direct) methods are available in library
TypeName Description
QR – uses the Householder rank-revealing QR decomposition (function colPiv-
HouseholderQr()) of Eigen [35] to decompose matrix Jk
a. It can be used
in any situation, and it should be the default choice in case of doubt.
LLT performs a standard Cholesky decomposition of matrix Jk
a, using the llt()
function of Eigen [35]. 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).
Gaussian uses the Gaussian elimination method provided by OpenFOAMR
. Restric-
tions are the same as for TDMA method.
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).
4.3.6 Tethering and fixing the molecules center of mass
rheoTool allows the simulation of tethered molecules, i.e., molecules with some
specific part of the chain fixed in space. The only restriction is that solely the
first bead of a given molecule (the one having index = 0 inside the chain) can be
fixed. The numerical procedure to simulate tethered molecules is similar to that
for generic molecules, with the exception that we impose a null resulting force
acting on the tethered bead.
CHAPTER 4. Overview of rheoTool 71
Fixing the center of mass of a molecule is different from tethering a molecule,
and the use of both in simultaneous is not allowed (would it make sense?). When
the option of fixing the molecules center of mass is selected, the molecules are
simulated in the usual way, but at the end of ntime-steps all the beads are
translated by a same vector, which restores the original center of mass (ntime-steps
before). This option is useful, for instance, when simulating molecules under an
analytical external forcing (homogeneous in space), or to equilibrate the molecules
in a local flow. The value of n, the number of consecutive time-steps in which
the molecules movement is unconstrained, is defined by the user. The existence of
this tunable parameter is mainly to avoid some of the inconsistencies previously
described in Section 4.3.4 for the BarycentricWeights interpolation method.
The reader should keep in mind that the molecules will always move in space
(over the mesh) according to the external forcing if none of these options is selected,
independently of the external forcing nature.
4.3.7 Beads tracking
The default particle tracking engine of OpenFOAM R
is used to move the beads
over the mesh. The algorithm has been greatly improved since OpenFOAMR
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.
4.3.8 Data output for post-processing
The sPCloudInterface class can optionally save some information about the
molecules in runtime. When this option is active, a directory named rheoTo
CHAPTER 4. Overview of rheoTool 72
olPP/startTime/moleculesStates/groupName is created for each group
of molecules, and three files are written to this directory:
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(|rjri|),
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 b