# User Guide

User Manual:

Open the PDF directly: View PDF .

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

- Contents
- 1 Introduction
- 2 Installation
- 3 Theoretical background
- 4 Overview of rheoTool
- 4.1 The constitutiveEquations library
- 4.2 The EDFModels library
- 4.3 The BDmolecule library
- 4.4 The sparseMatrixSolvers library
- 4.4.1 Conditions to reuse the preconditioner/factorization
- 4.4.2 Residuals and tolerances
- 4.4.3 Generic parameters
- 4.4.4 OpenFOAM interface
- 4.4.5 Eigen interface
- 4.4.6 Hypre interface
- 4.4.7 Petsc interface
- 4.4.8 Coupled solvers
- 4.4.9 How to use sparseMatrixSolvers library in my own application?
- 4.4.10 Limitations

- 4.5 Solvers
- 4.6 Boundary conditions
- 4.7 Utilities

- 5 Tutorials
- 5.1 rheoFoam
- 5.1.1 General guidelines
- 5.1.2 A note on coded FunctionObjects
- 5.1.3 Case 1: flow between parallel plates
- 5.1.4 Case 2: lid-driven cavity flow
- 5.1.5 Case 3: flow in a 4:1 planar contraction
- 5.1.6 Case 4: flow around a confined cylinder
- 5.1.7 Case 5: bifurcation in a 2D cross-slot flow
- 5.1.8 Case 6: blood flow simulation in a real-model aneurysm
- 5.1.9 Case 7: viscous fluid damper (moving mesh)

- 5.2 rheoTestFoam
- 5.3 rheoInterFoam
- 5.4 rheoEFoam
- 5.4.1 General guidelines
- 5.4.2 Case I: EDF of power-law and PTT fluids in a microchannel
- 5.4.3 Case II: induced-charge electroosmosis around a cylinder
- 5.4.4 Case III: charge transport across an ion-selective membrane
- 5.4.5 Case IV: electrokinetic instabilities in a flow-focusing device
- 5.4.6 Case V: electrokinetic mixer
- 5.4.7 Case VI: electro-elastic instabilities in cross-shaped geometries

- 5.5 rheoBDFoam

- 5.1 rheoFoam
- 6 FAQs
- Appendix A Parameters and variables in rheoTool
- Bibliography

User Guide

Version 4.0

April 26, 2019

License

This document is licensed under

Creative Commons Attribution-NonCommercial-NoDerivs 3.0 Unported License

http://creativecommons.org/licenses/by-nc-nd/3.0/legalcode

Acknowledgments

The work leading to the preparation of this document has received funding from

the European Research Council under the European Union’s Seventh Framework

Programme (FP7/2007-2013)/ERC Grant agreement 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.

Disclaimer

This oﬀering is not approved or endorsed by OpenCFD Limited, producer and

distributor of the OpenFOAM software via www.openfoam.com, and owner of the

OPENFOAM R

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

document.

RheoTool, where this guide is included, is distributed in the hope that it will be

useful, but WITHOUT ANY WARRANTY. See the GNU General Public

License (http://www.gnu.org/licenses/) for more details.

Trademarks

Linux is a registered trademark of Linus Torvalds.

OpenFOAM is a registered trademark of of OpenCFD Limited.

Paraview is a registered trademark of Kitware.

Typeset in L

A

T

EX.

c

2016-2019 Francisco Pimenta, Manuel A. Alves

Contents

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 Diﬀerences 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 ﬂuid ﬂows . . . . . . . . . . . . . . 15

3.2 Stabilization of viscoelastic ﬂuid ﬂow simulations . . . . . . . . . . 16

3.2.1 The both-sides-diﬀusion (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 ﬂow 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

ii

CONTENTS iii

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 ﬁxing 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

CONTENTS iv

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 writeEﬁeld ............................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: ﬂow between parallel plates . . . . . . . . . . . . . . 123

5.1.4 Case 2: lid-driven cavity ﬂow . . . . . . . . . . . . . . . . . 125

5.1.5 Case 3: ﬂow in a 4:1 planar contraction . . . . . . . . . . . . 126

5.1.6 Case 4: ﬂow around a conﬁned cylinder . . . . . . . . . . . . 129

5.1.7 Case 5: bifurcation in a 2D cross-slot ﬂow . . . . . . . . . . 132

5.1.8 Case 6: blood ﬂow simulation in a real-model aneurysm . . . 134

5.1.9 Case 7: viscous ﬂuid 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 ﬂuids 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 ﬂow-focusing device 163

5.4.6 Case V: electrokinetic mixer . . . . . . . . . . . . . . . . . . 167

CONTENTS v

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 ﬂow . . . . 177

5.5.4 Case 2: 7λ-DNA extension in a ﬂow-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

Introduction

1.1 Motivation

The open-source OpenFOAM R

toolbox can be used as a versatile ﬁnite-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 ﬂuids, 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) ﬂows 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

modiﬁed solver was tested in benchmark ﬂows 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 ﬂows was added to

rheoTool [3] and is available since version 2.0.

Recognizing the importance of modeling polymeric ﬂows at diﬀerent 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 ﬂows, which can

be also used for pressure-driven ﬂows. 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 ﬂows of both GNF

and viscoelastic ﬂuids, 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

ﬂuid ﬂows. In particular, some of the distinguishing features of rheoTool are:

•both GNF and viscoelastic models can be selected on run time and applied to

1

CHAPTER 1. Introduction 2

single-phase laminar ﬂows. A solver for two-phase ﬂows 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 ﬂows.

•a stress-velocity coupling term can be selected on run time in order to avoid

checkerboard ﬁelds under speciﬁc conditions, such as in the simulation of the

Upper-Convected Maxwell (UCM) model in strong extensional ﬂows.

•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

toolbox.

•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 ﬂow to be tested (shear ﬂow, extensional ﬂow, etc.).

•a number of models for electrically-driven ﬂows is available and can be cou-

pled with any rheological model. Mixed pressure- and electrically-driven

ﬂows are also allowed.

•transient ﬂow 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 ﬂow 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.

1http://openfoam.org/

2http://www.extend-project.de/

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. [2–4,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

rheoTool.

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

ﬁnite-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 diﬀerent OpenFOAM R

and foam-extend

versions, for historical reasons Chapters 4and 5are still mainly based on

OpenFOAM R

version 2.2.2 to describe the contents (except the content related

with Brownian dynamics simulations, described for OpenFOAM R

version 5.x).

However, the small diﬀerences among diﬀerent 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 ﬁrst read Refs. [2], [3] and [4] before this guide.

1.3 Changelog

Version 4.0

Released on 04/04/2019.

Generic

•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

version.

•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

version.

•Add: added coupled solvers for both pressure- and electrically-driven ﬂows.

Most viscoelastic models can be solved within a coupled solution method.

Only for the OpenFOAM R

version.

•Change: rheoTool version compatible with OpenFOAM R

v4.0/4.1 is discon-

tinued.

Electrically-driven ﬂows

CHAPTER 1. Introduction 4

•Add: added coupled Poisson-Nernst-Planck model to EDF models (Type-

Name = NernstPlanckCoupled). Only for the OpenFOAM R

version.

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 modiﬁed to allow integration with coupled solvers. Only for the

OpenFOAM R

version.

Solvers

•Fix: ﬁxed bug in rheoInterFoam for OpenFOAM R

version 6.0, which was

preventing post-processing (missing call to update()).

Tutorials

•Change: tutorial rheoFoam/Cavity now uses sparse matrix solvers from

Eigen library and is 1.6 times faster. Only for the OpenFOAMR

version.

•Change: tutorial rheoFoam/Cylinder is now solved coupled, being 30 times

faster. Only for the OpenFOAMR

version.

•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

version.

•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

version.

•Fix: ﬁxed bug related to the old ﬂag 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.

Generic

CHAPTER 1. Introduction 5

•Add: all solvers are now compatible with dynamic meshes. Due to this

change, and for convenience, momentum equation is the ﬁrst to be solved,

followed by pressure equation and then the equations for the remaining vari-

ables (extra-stresses, passive scalar, etc.).

•Add: tutorial ﬂuidDamper showing the use of rheoFoam with dynamic

meshes.

•Add: added an explicit Navier slip boundary condition for velocity.

•Change: Namespace encapsulation of several derived classes.

•Change: rheoTool version compatible with OpenFOAM R

v2.2.2 is discon-

tinued.

•Add: added rheoTool patch for OpenFOAM R

v6.0.

•Add: added note in the user-guide (Section 2.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-

able.

•Add: the Bautista-Manero-Puig (BMP) model has been added to the library

of constitutive equations. Both stress and log-conformation versions are

available.

•Change: implemented functions tauTotal() and tauTotalMF() in base classes.

Solver rheoTestFoam and the utility retrieving the wall shear-stresses were

modiﬁed accordingly.

Version 2.0

CHAPTER 1. Introduction 6

Released on 09/02/2018.

Electrically-driven ﬂows

•Add: solvers, libraries, utilities and tutorials for electrically-driven ﬂows.

Constitutive equations

•Add: the Rolie-Poly viscoelastic model has been added to the library of

constitutive equations. Both the stress and log-conformation versions are

available.

•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-aﬃne 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 veriﬁcation 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).

Post-Processing

•Add: class ppUtil for post-processing purposes has been added to the versions

for OpenFOAM R

and the one existing for foam-extend has been modiﬁed.

Enable the use of multiple ppUtil in simultaneous.

•Fix: sampling error was ﬁxed for the tutorials of versions of40 and fe40.

Multiphase ﬂows

•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 ﬂows (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 ﬁxedFluxPressure BC with rheoInterFoam

in versions of222 and fe40.

•Add: tutorials on the die swell problem.

Generic

CHAPTER 1. Introduction 7

•Change/Fix: code cleanup and bug ﬁx (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 ﬂows.

Other changes were made in its content and organization, and some typos

were corrected.

•Change: ensure compatibility with foam-extend 4.0 and OpenFOAMR

v4.1.

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:

@misc{rheoTool,

author = "F. Pimenta and M.A. Alves",

title = "rheoTool",

howpublished = "\url{https://github.com/fppimenta/rheoTool}",

year = "2016"}

Since the underlying theory of rheoTool has been mainly presented in technical

papers [2–5], 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: fpimenta@fe.up.pt

RManuel A. Alves: mmalves@fe.up.pt

1.6 Contributing

In the open-source spirit, rheoTool is open to contributions from the community.

If you believe that your piece of code is worth to be incorporated in rheoTool’s

next version, feel free to contact us.

Chapter 2

Installation

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 eﬀort 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:

•OpenFOAM R

v6.0 (of60/).

•foam-extend 4.0 (fe40/).

The list above includes the versions which were eﬀectively tested. This means

that a given version of rheoTool may be compatible with other OpenFOAM R

or

foam-extend versions not included in this list. The versions above were tested in

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 Diﬀerences 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 diﬀerences, in order to make rheoTool compatible with

each OpenFOAM R

/foam-extend version, several modiﬁcations were required at

the programming level for each case. Still, the user-interface remained almost

unchanged among the diﬀerent 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

,

8

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

versions

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 diﬀer-

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 diﬀerences 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.

OpenFOAM R

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 fulﬁlled, 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 ﬁle of the repository, rheoTool can be either down-

loaded or cloned from GitHub. The structure of rheoTool (https://github.com

/fppimenta/rheoTool) is depicted in Fig. 2.1.

rheoTool ofxx

feyy

doc

src

tutorials

libs

solvers

(…) (…)

(…)

(…)

Figure 2.1: Directory organization of rheoTool.

The top-level directory of rheoTool contains the versions available for diﬀerent

OpenFOAM R

(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 ﬁnd ﬁle downloadEigen), open a terminal and check that ﬁle etc/bashrc

of your installed OpenFOAM R

or foam-extend version has been sourced. This is

particularly relevant if you have deﬁned alias for diﬀerent 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:

∼$./downloadEigen

This script downloads Eigen version 3.2.9 (other versions close to that would

also work adequately) from the Internet (using wget), extracts it and moves it to

directory:

$WM PROJECT USER DIR/ThirdParty/Eigen3.2.9

Eigen is used in rheoTool for computation of eigenvalues and eigenvectors and

there is no need to install the library, since the inclusion of the required headers

is enough for our purposes.

However, its location in the system must deﬁned and exported. This is achieved

by attributing to variable EIGEN_RHEO – the one used and recognized by rheoTool –

the actual path of Eigen. The command to do so has been displayed to the terminal

after running script downloadEigen (if everything was ok) and looks like:

∼$echo "export EIGEN_RHEO=/home/user/OpenFOAM/user-4.0/ThirdParty

/Eigen3.2.9">>/home/user/.bashrc

Do not copy this command, it is just an example of what is displayed to the

screen. Instead, copy-paste and run the command appearing in your terminal.

If, for some reason, the user wants to move Eigen to another directory (or

already has an Eigen version in another directory), then move Eigen to its ﬁnal

location (if already not) and deﬁne 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 ﬁrst version being installed), as long as the

directory containing Eigen since the ﬁrst installation is not deleted, or renamed.

2.4.3 Install Petsc library

This step is only for OpenFOAM R

users.

Petsc [9–11] has been added to the dependencies used by rheoTool starting

from version 4.0. This library provides a number of eﬃcient and scalable sparse

matrix solvers, that can be used in rheoTool with both segregated and coupled

solvers.

The script installPetsc is responsible for downloading Petsc and the de-

pendencies it needs, conﬁguring 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:

libatlas-dev,libatlas-base-dev,libblas-dev,liblapack-dev,flex,bi

son,git,make,cmake,gfortran

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-

ﬁrmation). 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

users.

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

OpenFOAM R

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 diﬀerent from OpenMPI. At the beginning of script inst

allPetsc, the user can change the paths to the MPI library and wrappers upon

need.

Petsc library is conﬁgured with standard options, that the user can found and

modify under section ## Conﬁgure petsc. To check the full range of options

available, please consult Petsc references [9,10]. The conﬁguration that we set as

default will download and install the following additional packages from within

Petsc:

hypre,parmetis,metis,ptscotch,mumps,scalapack

Note that Petsc could be also installed via apt-get (there are binaries available),

but this would not allow conﬁguring and using the most recent versions of the

software.

For most of the users, the script to install Petsc can be run without any mod-

iﬁcation. In such cases, go the directory corresponding to one of the rheoTool

versions (where you will ﬁnd ﬁle installPetsc), open a terminal and ensure that

ﬁle etc/bashrc of your installed OpenFOAM R

version was sourced (as for Eigen).

Then run the script installPetsc in that terminal:

∼$./installPetsc

Petsc library is saved and compiled to directory:

$WM PROJECT USER DIR/ThirdParty

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 modiﬁes

your ~/.bashrc ﬁle 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 conﬂicts 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-

ﬁned 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 deﬁned

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 diﬀerent

OpenFOAM R

versions, since in these situations sourcing the etc/bashrc of

OpenFOAM R

cleans $LD LIBRARY PATH. When this happens, running the

Allwmake script as described below results in linking errors similar to:

/usr/bin/ld: warning: libpetsc.so.3.10, needed by /home/user/OpenFOAM/user-6/

platforms/linux64GccDPInt32Opt/lib/libconstitutiveEquations.so, not found (

try using -rpath or -rpath-link)

/home/user/OpenFOAM/user-6/platforms/linux64GccDPInt32Opt/lib/

libsparseMatrixSolvers.so: 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, ﬁnd your alias (typically at the bottom of the ﬁle) and

modify it appending command

export LD_LIBRARY_PATH=$PETSC_DIR/$PETSC_ARCH/lib:$LD_LIBRARY_PATH

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 modiﬁcation should look like this (single line)

alias of60=’source /opt/openfoam6/etc/bashrc; export LD_LIBRARY_PATH=$PETSC_DIR

/$PETSC_ARCH/lib:$LD_LIBRARY_PATH’

after the modiﬁcation. 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 ﬁle of OpenFOAM R

, creating symbolic links of PETSc

libraries inside one of the locations included by default by OpenFOAM R

in

$LD LIBRARY PATH, etc.

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

OpenFOAM R

’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 ﬁnal 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

speciﬁed, 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 deﬁned in the etc/bashrc ﬁle of your OpenFOAM R

installation. For

foam-extend users, option −jis not recognized by the script, therefore simply run

∼$./Allwmake

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

way.

Both the libraries and applications installed with rheoTool can be ”cleaned”

by running the script Allwclean.

Since the user will probably not need the remaining versions of rheoTool that

remain in the top-level directory, they can simply be deleted, if already not.

To check if the installation succeeded, the user should try running one of the

tutorials in Chapter 5.

Chapter 3

Theoretical background

The equations governing pressure- and electrically-driven ﬂows of incompressible,

complex ﬂuids are discussed in this Chapter, along with some important aspects

related with their discretization in the ﬁnite-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 ﬂuid ﬂows

The basic equations governing isothermal, single-phase, transient ﬂows, under lam-

inar conditions, for incompressible ﬂuids, in static grids, establish mass conserva-

tion (Eq. 3.1) and momentum balance (Eq. 3.2),

∇·u= 0 (3.1)

ρ∂u

∂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 ﬂuid ﬂows, 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)

15

CHAPTER 3. Theoretical background 16

In Eqs. (3.3) and (3.4), ηsis the solvent viscosity, ηpis the polymeric viscosity

coeﬃcient, λ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

2(∇u+∇uT))

instead of the upper-convected derivative, in order to take non-aﬃne 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

speciﬁed 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 ﬂuid ﬂow 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 ﬂuid ﬂow simu-

lations

3.2.1 The both-sides-diﬀusion (BSD) technique

The both-sides-diﬀusion (BSD) is a technique already incorporated in the vis-

coelasticFluidFoam solver [1]. It consists in adding a diﬀusive term on both sides

of momentum equation (Eq. 3.2), with the diﬀerence 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 eﬀect, mostly when there is no solvent contribution in the extra-stress tensor.

Incorporating the terms arising from the both-sides-diﬀusion in the momentum

equation, and making use of Eq. (3.3), then

ρ∂u

∂t +u·∇u− ∇·(ηs+ηp)∇u=−∇p− ∇·(ηp∇u) + ∇·τ+f(3.5)

Note that the added diﬀusive 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 ﬂows [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)

τ=ηp

λ(A−I) (3.6)

In the log-conformation tensor methodology, a new tensor (Θ) is deﬁned 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 deﬁnite, 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

[6]

∂Θ

∂t +u·∇Θ=ΩΘ −ΘΩ + 2B+1

λg(Θ) (3.8)

where g(Θ) is a model-speciﬁc tensorial function depending on Θ(see Table 4.1

for other viscoelastic models) and

B=R

mxx 0 0

0myy 0

0 0 mzz

RT(3.9)

Ω=R

0ωxy ωxz

−ωxy 0ωyz

−ωxz −ωyz 0

RT(3.10)

M=R∇uTRT=

mxx mxy mxz

myx myy myz

mzx mzy mzz

(3.11)

ωij =Λjmij +Λimji

Λj−Λi

(3.12)

After solving Eq. (3.8), Θis diagonalized in the form

Θ=RΛΘRT(3.13)

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-aﬃne deformation

through the Gordon-Schowalter derivative, the tensor M(Eq. 3.11) is computed

diﬀerently: M=R∇uT−ζ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 ﬂuid ﬂows 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

∇·1

aP−H1

(∇p)P=∇·H

aP

+1

aP−H1−1

aP(∇p∗)P(3.15)

where aPare the diagonal coeﬃcients from the momentum equation, H1=−P

nb

anb

is an operator representing the negative sum of the oﬀ-diagonal coeﬃcients from

momentum equation, H=−P

nb

anbu∗

nb +bis an operator containing the oﬀ-

diagonal contributions, plus source terms (except the pressure gradient) of the

momentum equation and p∗is the pressure ﬁeld known from the previous time-

step or iteration. Accordingly, the equation to correct the velocity after obtaining

the continuity-compliant pressure ﬁeld from Eq. (3.15) is

u=H

aP

+1

aP−H1−1

aP(∇p∗)P−1

aP−H1

(∇p)P(3.16)

Importantly, in order to avoid the onset of checkerboard ﬁelds, 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 ﬁeld, the pressure gradient terms

are computed ”in the usual way”, for example using Green-Gauss integration.

Rhie-Chow methods used to avoid checkerboard ﬁelds, as the one described

in the previous paragraph, are known to be aﬀected by the use of small time-

steps and they also present time-step dependency on steady-state results [16].

In OpenFOAM R

solvers, a common strategy to avoid such eﬀects 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

diﬀerent way, by removing the transient term contribution from the aPcoeﬃcients

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 eﬃciently 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 inﬂuence 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 ﬁelds. 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+ηph∇u|f+(∇u)T|f−∇u|f+(∇u)T|fi (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 deﬁnition of τfin Eq. (3.17) is inserted in

the momentum equation with the both-sides-diﬀusion terms already present (Eq.

3.5), then we obtain

ρ∂u

∂t +u·∇u− ∇·(ηs+ηp)∇u=−∇p− ∇·ηp∇u+∇·τ+f(3.18)

where the term ∇·ηp∇uis a ”special second-order derivative” (diﬀerent from

the laplacian operator of OpenFOAM R

), deﬁned 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 reﬁnement 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 ﬁnite-volume framework leads to

ZV

(u·∇φ) dV=X

f

φf(uf·Sf) = X

f

φfFf(3.19)

where φis a generic variable being advected, Sfis the face-area vector and Ffis

the volumetric ﬂux crossing face f. While ﬂuxes 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

oﬀers a wide range of schemes to

perform such interpolation, from upwind – an unconditionally stable scheme, but

only ﬁrst-order accurate –, to central diﬀerences – 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 deﬁned

using the Normalized Weighting Factor (NWF) approach [18]:

e

φf=αe

φC+β(3.20)

where the following deﬁnitions hold

e

φf=φf−φU

φD−φU

(3.21a)

e

φC=φC−φU

φD−φU

(3.21b)

In Eq. (3.20), αand βare scalars speciﬁc 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 ﬂux comes (upstream), cell D (downstream) is the

cell to which the ﬂux goes and cell U (far-upstream) is the cell upstream to cell C.

In a general unstructured mesh, cell U cannot be identiﬁed unequivocally, and φU

in Eqs. (3.21a,b) can be evaluated as [19]

φU=φD−2(∇φ)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 (diﬀerence between the HRS and the upwind diﬀer-

encing scheme) is discretized explicitly (cf. Ref. [2]), which, using Eqs. (3.20-3.22),

results in

φf= [φC]implicit + [(α−1)φC+βφD+ (1 −α−β)(φD−2(∇φ)C·dCD)]explicit (3.23)

Handling the HRSs in a deferred correction approach avoids, in some cases,

numerical instabilities introduced by the central-diﬀerencing component of the

CHAPTER 3. Theoretical background 21

HRS. Additionally, in Ref. [2] it was observed that the usual methodology of

OpenFOAM R

to apply HRSs to non-scalar variables (tensors and vectors) can

locally introduce numerical instabilities in some viscoelastic ﬂow 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 artiﬁcial instabilities can be signiﬁcantly 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 reﬁnement. 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-

eﬃcients for such variables, since the upwind diﬀerencing scheme coeﬃcients are

common to all the components (they only depend on the ﬂux). The diﬀerentiation

between components is only introduced in the explicit part of Eq. (3.23), generat-

ing a diﬀerent 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

ﬂuid. 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 ﬂow induced inside a sphere due to its time-dependent rotation. Such case can

be easily handled by deﬁning adequate boundary conditions for the ﬂow variables

on the sphere surface, without further modiﬁcations of the usual solver setup. On

the other hand, if we consider the time-dependent simulation of the ﬂow 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 diﬀerent 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 ﬂow need to be

changed regarding convective terms, which should account for the grid motion [20],

ZS

φ(u−ub)·ndS=X

f

φf(uf−ub,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

satisﬁed to ensure mass conservation [20],

d

dtZV

dV−ZS

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 ﬂux due to mesh motion. According to Eq. (3.25), the

form taken by this term involving the volume swept by the moving faces at diﬀerent

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 reﬁnement (AMR), frequently used to locally (un)reﬁne 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 deﬁning the ﬁelds

and their ﬂuxes 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 ﬁnite-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 ﬂow simulations [5].

In [5] we discussed the implementation of coupled and semi-coupled solvers in

rheoTool, in the context of electrically-driven ﬂows. 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 ﬂows.

3.7 Electrically-driven ﬂow models

Consider now that the ﬂuid under analysis is a weak electrolyte subjected to an

electric ﬁeld. In such conditions, the momentum equation (Eq. 3.2) should include

the contribution from an electric body-force,

f=fE=∇·εEE −kEk2

2I=ρEE−kEk2

2∇ε(3.26)

CHAPTER 3. Theoretical background 23

where Eis the electric ﬁeld, ε=ε0εRis the electric permittivity and ρEis the

charge density (per unit volume). In order to close the system of equations for

electrically-driven ﬂows (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

diﬀerence), in addition to the electric forcing. When only an electric forcing exists,

we call this ﬂow as pure EDF.

The second term of Eq. (3.26) is only non-zero for a system of two ﬂuids, each

having a diﬀerent electric permittivity.

3.7.1 Poisson-Nernst-Planck model

In the absence of magnetic eﬀects, the electric potential (Ψ) can be computed by

Gauss’ law

∇·(ε∇Ψ) = −ρE(3.27)

where the electric ﬁeld is E=−∇Ψin electrostatics. By deﬁnition, the charge

density is

ρE=F

N

X

i=1

zici(3.28)

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 ﬁeld and neglecting any reaction,

is embodied by the Nernst-Planck equation,

dci

dt +u·∇ci=∇·(Di∇ci) + ∇·

Di

ezi

kT ∇Ψ

| {z }

uM,i

ci

(3.29)

which closes the system of equations for an EDF. In Eq. (3.29), Dis the diﬀu-

sion coeﬃcient, 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 ﬁeld, 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 ﬁeld Ψ, with a space and time

varying diﬀusion coeﬃcient, 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 stiﬀ system of equations

when the PNP model is used. As such, several simpliﬁed 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 eﬀects (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 diﬃculties when deﬁning 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 simpliﬁcation 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 ﬁeld deﬁnition.

This can be justiﬁed 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 aﬀect the ﬂow.

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 simpliﬁed to the so-called Poisson-Boltzmann model (henceforth

PB model), for which Gauss’ law reads

∇·(ε∇ψ) = −F

N

X

i=1

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 deﬁnition 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 inﬂuence of ﬂow 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 ﬁrst-derivative, transforming the

equation into

∇·(ε∇ψ) + ψF

N

X

i=1

(aibi)∗=−F

N

X

i=1

(ai)∗+ψ∗F

N

X

i=1

(aibi)∗(3.32)

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

N

X

i=1

zici,01−ezi

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

inﬂuence of ﬂow 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,

λD=v

u

u

u

t

εkT

F e

N

P

i=1

z2

ici,0

(3.34)

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 ﬂow inside the EDL, then it is possible to approximate the EDL

eﬀect by a slip velocity at the surface, avoiding the need to solve the ﬂow inside the

EDL. Such a case would be, for example, the pumping of a Newtonian electrolyte

(λD∼ O(10−9m)) in a microchannel of arbitrary shape (W∼ O(10−6m)), 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 ﬁeld

at the surface must be known. The electroosmotic mobility is assumed to be known

a priori – it can be a ﬁxed value over all the surface or have a known distribution.

On the other hand, the electric ﬁeld 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 eﬀects 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 ﬂows, even though the condition

λD

W1 is satisﬁed. 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 ﬂuid ﬂows [21] – using a slip model would simply retrieve a smooth

ﬂow in such cases.

3.7.6 Ohmic (leaky dielectric) model

The so-called Ohmic model [22] is particularly useful to simulate ﬂuids of diﬀerent

conductivities, although a generalized Ohmic model has been recently proposed

for diﬀerent 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 ﬁnal 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

D−.

First, let’s start deﬁning the conductivity (σ) and free-charge density (ρE) for

a binary electrolyte,

σ=F2z2

RT (D+c++D−c−) (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·∇σ=Deﬀ∇2σ(3.39)

∇·(σ∇Ψ) = 0 (3.40)

where the eﬀective diﬀusivity is Deﬀ =2D−D+

D−+D+. The conductivity is transported

through Eq. (3.39), while Eq. (3.40), derived from the conservation of charge-

density (then simpliﬁed 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

fE=ρEE=∇·(ε∇Ψ)∇Ψ(3.41)

In order to close the Ohmic model, the EDL eﬀect is commonly represented

by a slip velocity, which avoids detailing the ﬂow inside the EDL using a very

ﬁne 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σ

σ0m

E(3.42)

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 ﬂuid ﬂows 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= (N−1) 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|=|rj−ri|, 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

ﬁxed orientation. The Kuhn step size, deﬁned as bk= 2λPis a measure commonly

used in coarse-grained models, especially in bead-rod models, where it represents

the ﬁxed 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

(ﬂexible) 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

Molecule

Group of Molecules

Group of Molecules

Group of Molecules

Bead

Spring

Branch

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

[24,25]

∂ri

∂t =uf+

N

X

j=1

∂Dij

∂rj

+

N

X

j=1

DijFj

kT +6

∆t0.5i

X

j=1

σijnj(3.43)

where ufis a velocity imposed by an external forcing, Dij is the diﬀusion tensor,

kis Boltzmann’s constant, Tis the absolute temperature, Fjis the sum of spring

and exclusion volume (EV) forces (Fj=FS

j+FEV

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 diﬀusion

tensor is the Rotne–Prager–Yamakawa (RPY) model [4],

Dij =

kT

6πηaI=DI, i =j

3D

4

a

|xij |h1 + 2a2

3|xij |2I+1−2a2

|xij |2xij xij

|xij |2i, i 6=j∧ |xij | ≥ 2a

Dh1−9|xij |

32aI+3

32a

xij xij

|xij |i, i 6=j∧ |xij|<2a

(3.44)

where |xij|=|rj−ri|,ηis the ﬂuid viscosity, ais the bead radius and Iis the unit

tensor (3 ×3). Note that in rheoTool,aand Dare deﬁned 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 ﬁeld. In

such situations, hydrodynamic interactions are neglected and all the oﬀ-diagonal

tensor elements of tensor Dbecome zero. The diﬀusion becomes isotropic and

can simply be deﬁned by coeﬃcient 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],

FEV

i=−9

2

kT

l

νEV

l33

4√π3

(2Nk,s)9/2

N

X

j=1

exp −9

2Nk,s |xij|2

l2xij

l(3.45)

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 ﬁrst category, while the Hookean model falls in the second one.

The Hookean model is arguably the simplest model representing a spring [24],

FS

i=

N

X

j∈gi

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 eﬀective 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 diﬃculties 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

FS

i=

N

X

j∈gi

2

3Hl "|xij|

l−1

4+1

4 (1 − |xij |/l)2#xij

|xij|(3.47)

•Cohen Pad´e model

FS

i=

N

X

j∈gi

H

33−(|xij|/l)2

1−(|xij|/l)2xij (3.48)

•FENE model (Warner spring law)

FS

i=

N

X

j∈gi

H1

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 diﬀer essentially in numerical stability and

CHAPTER 3. Theoretical background 32

computational cost. Most of the methods suggested in the literature are ﬁrst-order

accurate in time, using Euler schemes to discretize the time derivative (higher-order

methods are not eﬀective due to the random nature of the Brownian term [26]).

The explicit ﬁrst-order Euler method evolves the beads positions from the

previous time-step (t) to the new time-step (t+ ∆t) as,

rt+∆t

i=rt

i+ ∆tBt

i(3.50)

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 suﬃciently small. The numerical stability is evaluated by the

capability of the method in respecting the constraint Ri≤l, when bounded spring

force models are used. In practice, the time-step that satisﬁes 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 deﬁned by the user, and is usually

close to 1. Once this condition is violated, a two-steps computation is used:

r∗

i=rt

i+ ∆t"uf+

N

X

j=1,j6=i

DijFS

j

kT +DiiFEV

i

kT +6

∆t0.5i

X

j=1

σijnj#t

(3.51)

rt+∆t

i=r∗

i+ ∆tDiiFs

i

kT t+∆t

(3.52)

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

oﬀ-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),

Jk∆rk=−fk(3.53)

where fkis a vector function whose expression depends on the spring model, Jk

is the Jacobian of fkand ∆rkis a vector representing the diﬀerence 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 diﬀerent

methods.

Chapter 4

Overview of rheoTool

In the previous Chapter, the main theoretical points behind rheoTool were brieﬂy

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.

33

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 ˙γn−1

Carreau-Yasuda CarreauYasuda η∞+ (η0−η∞)[1 + (k˙γ)a]n−1

a

2Herschel-Bulkley HerschelBulkley Bounded: min η0, τ0˙γ−1+k˙γn−1

3Papanastasiou reg.: min η0, τ0˙γ−1[1 −exp(−m˙γ)] + k˙γn−1

2Casson Casson

Bounded: max ηmin,min ηmax,√η∞+qτ0

˙γ2

Papanastasiou reg.: n√η∞+qτ0

˙γ1−exp −√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 n−1<0. For ˙γ < VSMALL, the

value ˙γ=VSMALL is used in the computation of the shear viscosity (VSMALL = 10−300 for versions using double precision).

3The original Papanastasiou regularization does not include the artiﬁcial upper-bounding by η0. However, this bounding is needed

to avoid an inﬁnite viscosity for ˙γ→0 (e.g. startup of ﬂow) and n < 1. The original Papanastasiou regularization is recovered

for η0→ ∞. In practice, η0should be low enough to avoid an inﬁnite viscosity in quiescent conditions and high enough to allow

Papanastasiou regularization to take control in the remaining situations.

Notes:

•˙γ=q˙

γ:˙

γ

2, with ˙

γ=∇u+∇uT.

•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

2(∇u+∇uT).

34

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λτ+λ∇

τ=ηp(∇u+∇uT)

WhiteMetzner

(Carreau-Yasuda) WhiteMetznerCY ηsηp[1 + (K˙γ)a]n−1

aλ[1 + (L˙γ)b]m−1

bτ+λ( ˙γ)∇

τ=ηp( ˙γ)(∇u+∇uT)

Giesekus Giesekus ηsηpλτ+λ∇

τ+αλ

ηp(τ·τ) = ηp(∇u+∇uT)

PTT linear PTTlinear ηsηpλh1 + ελ

ηptr(τ)iτ+λ

τ=ηp(∇u+∇uT)

PTT exponential PTTexp ηsηpλe

ελ

ηptr(τ)τ+λ

τ=ηp(∇u+∇uT)

FENE-CR FENE-CR ηsηpλh1 + λD

Dt1

fiτ+λ

f∇

τ=ηp(∇u+∇uT)

where f=L2+λ

ηptr(τ)

L2−3

FENE-P FENE-P ηsηpλ

τ+λ

f∇

τ=aηp

f(∇u+∇uT)−D

Dt1

f[λτ+aηpI]

where f=L2+λ

aηptr(τ)

L2−3and a=L2

L2−3

3Rolie-Poly Rolie-Poly ηsηpλD

λD∇

A=−(A−I)−2kλD

λR1−p3/tr(A)A+βtr(A)

3δ(A−I)

where k=3−χ2

χ2

max 1−1

χ2

max

1−χ2

χ2

max 3−1

χ2

max and χ=qtr(A)

3

eXtended Pom-Pom XPomPom ηsηpλB

fτ+λB∇

τ+αλB

ηp(τ·τ) + ηp

λB(f−1) I=ηp(∇u+∇uT)

where f= 2λB

λSe2

q(Λ−1) 1−1

Λn+1 +1

Λ2h1−α

3

tr(τ·τ)

(ηP/λB)2iand Λ=q1 + tr(τ)

3ηP/λB

3See Ref. [27]. This model is exclusively solved in the conformation tensor variable, which is then converted to τusing, τ=ηp

λDk(A−I).

35

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

λ(eΘ−I)Υ=1

λe−Θ−I

7WhiteMetzner

(Carreau-Yasuda) WhiteMetznerCYLog τ=ηp

λ(eΘ−I)Υ=1

λ( ˙γ)e−Θ−I

Giesekus GiesekusLog τ=ηp

λ(eΘ−I)Υ=1

λhe−Θ−I−αeΘe−Θ−I2i

PTT linear PTTlinearLog τ=ηp

λ(1−ζ)(eΘ−I)Υ=1

λn1 + ε

1−ζtr(eΘ)−3o(e−Θ−I)

PTT exponential PTTexpLog τ=ηp

λ(1−ζ)(eΘ−I)Υ=1

λeε

1−ζ(tr(eΘ)−3)(e−Θ−I)

FENE-CR FENE-CRLog τ=ηpf

λ(eΘ−I)Υ=f

λe−Θ−I, where f=L2

L2−tr(eΘ)

FENE-P FENE-PLog τ=ηp

λ(feΘ−aI)Υ=1

λae−Θ−fI, where a=L2

L2−3and f=L2

L2−tr(eΘ)

8Rolie-Poly Rolie-PolyLog τ=ηp

λDk(eΘ−I)Υ=−1

λDe−Θ(eΘ−I)+2kλD

λR1−p3/tr(eΘ)eΘ+βtr(eΘ)

3δ(eΘ−I)

eXtended Pom-Pom XPomPomLog τ=ηp

λB(eΘ−I)

Υ=−1

λBe−Θ(f−2α)eΘ+αeΘeΘ+ (α−1)I

where f= 2λB

λSe2

q(Λ−1) 1−1

Λn+1 +1

Λ21−α−α

3tr(eΘ(eΘ−2I))and Λ=qtr(eΘ)

3

‡The solvent viscosity, the polymeric viscosity coeﬃcient 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−Θ=A−1=RΛ−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Θ.

36

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 deﬁned as

˙γ=r˙

γ:˙

γ

2=√2D:D,with ˙

γ=∇u+∇uTand D=1

2˙

γ(4.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

2∇u+∇uT=1

2˙

γ=D

thus,

sqrt(2.0)*mag(symm(fvc::grad(U()))) =√2q1

2˙

γ:1

2˙

γ=q˙

γ:˙

γ

2=√2D:D

which is equal to Eq. (4.1) – the deﬁnitions 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 ﬂows, 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 ﬂuids. The model represents such ﬂuids 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],

∂nA

∂t +u· ∇nA= 2DA∇2nA+1

2λA

cBn2

B−cAnA

λA

(4.2)

∂nB

∂t +u· ∇nB= 2DB∇2nB−cBn2

B

λA

+ 2cAnA

λA

(4.3)

CHAPTER 4. Overview of rheoTool 38

where nis the dimensionless number density of the specie, λis the relaxation time,

Dis the diﬀusivity coeﬃcient and cAand cBare, respectively, the dimensionless

breakage and reformation rates, expressed as

cA=cAEq +χ

3˙

γ:A

nA(4.4)

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 deﬁnition of no-ﬂux 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],

λA∇

A+A−nAI−λADA∇2A=cBnBB−cAA(4.6)

λA∇

B+B−nBI

2−λADB∇2B=−2cBnBB+ 2cAA(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 ﬂow, nA= 1, nB=q2cAEq

cBEq

,

A=I,B=nB

2Iand τ=0.

+RRM (TypeName: RRM )

Also in the context of modeling the ﬂow of wormlike micellar solutions, Dutta

and Graham proposed recently the Reactive Rod Model (RRM), accounting for

the formation/destruction of ﬂow-induced structures [29]. In this model, micelles

are approached by rods, whose orientation tensor, S, evolves according to [29],

dS

dt=−6DrS−I

3+∇uT·S+S· ∇u−2∇uT:<uuuu >(4.9)

where

Dr=Dr,0

L∗3ln L∗+m

m(4.10)

is a time-varying diﬀusion coeﬃcient (units are s−1), 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·S−S·S·D−D·S·S+ 2S·D·S+ 3(S:D)S]

(4.11)

CHAPTER 4. Overview of rheoTool 39

and D=1

2∇u+∇uTis the rate of deformation tensor. The variation of the

normalized rod length is computed from [29]

dL∗

dt=λSDr,0

1−L∗

α+β/P e 2(1 −L∗) + kDr,0r3

2ˆ

S:ˆ

S(4.12)

where λS,α,βand kare parameters of the model, ˆ

S=S−I

3and P e is a P´eclet

number computed locally, P e =˙γ

Dr,0( ˙γis the strain-rate deﬁned in Eq. 4.1). The

extra-stress due to the rods is accounted for as [29]

τ=G0

L∗3S−I

3+1

2Dr∇uT:<uuuu >(4.13)

where G0is the elastic modulus. Note that due to the generic strain-rate deﬁ-

nition used in the P´eclet number, the strain-rate retrieved for a pure shear-ﬂow

corresponds eﬀectively to the local shear-rate, but for a pure extensional ﬂow, the

strain-rate computed does not correspond to the extension rate (it diﬀers 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 deﬁnition in order to generalize it for

any ﬂow.

+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τ+∇

τ=G0(∇u+∇uT) (4.14)

and Fredrickson’s kinetic equation governing the structure parameter (ϕ) is repre-

sented by

∂ϕ

∂t +u· ∇ϕ=ϕ0−ϕ

λ+k(ϕ∞−ϕ)τ:D(4.15)

In Eqs. (4.14) and (4.15), D=1

2(∇u+∇uT), G0is the instantaneous relaxation

modulus, ϕis the ﬂuidity (≡η−1

P), ϕ0is the zero shear-rate ﬂuidity, ϕ∞is the

inﬁnite shear-rate ﬂuidity, λ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 ﬂuids exhibit a solid-like behavior below the yield stress and

they ﬂow as viscoelastic ﬂuids 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 ﬂuids with the Oldroyd-B/PTT model for viscoelastic ﬂuids.

Accordingly, the resulting constitutive relation is given by

f(τ)ηpmax 0,σ−τ0

kσn1

n

τ+λ

τ=η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 deﬁned, G=ηp

λ, which is often found in the literature in the description of

this model. We remember that

τrepresents the Gordon-Schowalter derivative

deﬁned 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

eελ

ηptr(τ), exponential PTT

(4.17)

The four variants of the model are available in rheoTool, where both the Her-

schel–Bulkley and PTT variants can avoid the possible inﬁnite Oldroyd-B elonga-

tional viscosity for W i ≥0.5 in extensional ﬂows of Oldroyd-B ﬂuids. In addition,

the four variants can also be solved with the log-conformation approach and the

governing equation becomes

∂Θ

∂t +u·∇Θ−(ΩΘ −ΘΩ)−2B=g(τ)

λe−Θ−I(4.18)

with

g(τ) =

ηpmax 0,σ−τ0

kσn1

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

σeε

1−ζ(tr(eΘ)−3) , exp. PTT–Bingham

(4.19)

and the polymeric extra-stress tensor is recovered using τ=ηp

λ(1−ζ)(eΘ−I).

+ML-IKH model (TypeName: ML-IKH )

Elastoviscoplastic models can be rendered more complex (and realistic) once

thixotropy is added to them. This is the case of the Multi-Lambda Isotropic Kine-

matic Hardening (MLK-IKH) model [33]. According to this model, the equation

governing the polymeric extra-stress tensor is given by [33]

max 0,1−λky

στeﬀ +λλE∇

τ=ληp(∇u+∇uT) (4.20)

CHAPTER 4. Overview of rheoTool 41

where λEis the viscoelastic relaxation time, τeﬀ =τ−κback, with κback =

khA+ 2A2, and σ=IIτD

eﬀ =qτD

eﬀ:τD

eﬀ

2, where τD

eﬀ =τeﬀ −tr(τeﬀ)

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]

∂A

∂t +u· ∇A+Ω·A−A·Ω=Dp·A+A·Dp+Dp−qdpA(4.21)

where Ω=∇u− ∇uT/2 is the vorticity tensor, dp=qDp:Dp

2and

Dp=(0, if σ < λky

(σ−λky)

2ληp·τeﬀ

σ, if σ≥λky

(4.22)

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],

λ=

M

X

i=1

Ciλi(4.23)

where each of the Mmodes of λobeys

dλi

dt =Di−k1φ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 ﬂuid 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 simpliﬁcations, 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 modiﬁed

formulation in τare:

•FENE-CR

–Complete in A:

λ∇

A=−f(tr(A))(A−I), where f(tr(A)) = L2

L2−tr(A)

–Complete in τ(see Table 4.1):

h1 + λD

Dt1

fiτ+λ

f∇

τ=ηp(∇u+∇uT), where f=L2+λ

ηptr(τ)

L2−3

–Modiﬁed in τ(usually known as FENE-MCR):

τ+λ

f∇

τ=ηp(∇u+∇uT), where f=L2+λ

ηptr(τ)

L2−3

•FENE-P

–Complete in A:

λ∇

A=−[f(tr(A))A−aI], where f(tr(A)) = L2

L2−tr(A)and a=L2

L2−3

–Complete in τ(see Table 4.1):

τ+λ

f∇

τ=aηp

f(∇u+∇uT)−D

Dt1

f[λτ+aηpI], where f=L2+λ

aηptr(τ)

L2−3

and a=L2

L2−3

–Modiﬁed in τ:

τ+λ

f∇

τ=aηp

f(∇u+∇uT), where f=L2+λ

aηptr(τ)

L2−3and a=L2

L2−3

In rheoTool, all the formulations are available and can be used (see Section

5.1.1 to know how to select each one). The steady material functions evaluated

for canonical ﬂows are the same for all the formulations. However, this is not true

when evaluating the transient material functions: the modiﬁed formulations have

a diﬀerent behavior comparing with the complete ones, which are themselves simi-

lar. For a generic ﬂow, 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 modiﬁed formula-

tions, they are not expected to behave exactly as the complete ones, even in the

limit of highly reﬁned 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 modiﬁed 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

τ0=

N

X

k=1 τk+τk

s(4.26)

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

code.

The source code displayed in Listing 4.1 is taken from ﬁle src/libs/constitut

iveEquations/constitutiveEqs/Oldroyd-B/Oldroyd-BLog/Oldroyd_BLog.C.

Let’s analyze the most important lines:

•lines 1-91: this section initializes the variables used in the constitutive model.

In terms of ﬁeld variables, we have (lines 23-84): tau (τ), theta (Θ), eigVals

(Λ) and eigVecs (R). All those ﬁelds must be deﬁned by the user when

CHAPTER 4. Overview of rheoTool 44

starting a simulation, except eigVals and eigVecs , which can be deﬁned

or not. If deﬁned (typical of a restart from a previous simulation), they are

used in the ﬁrst time-step; otherwise, they are both initialized as the identity

tensor/matrix, corresponding to a null extra-stress tensor (τ). Afterwards,

the ﬂuid 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 ﬁeld, by evolving Θaccord-

ing to the constitutive equation. From line 96 to 105, variables M,Ωand

B, deﬁned in Eqs. (3.9)–(3.11), are computed. The function decomposeG-

radU() is a member function of the base class constitutiveEq (ﬁnd it in the

ﬁle 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

below.

•in the viscoelasticTransportModels library [1] each model was in charge to

deﬁne its own contribution to the momentum equation, i.e., the term (∇·τ0).

In the constitutiveEquations library there is a default deﬁnition of this term

in the base class. In fact, the function divTau() is now deﬁned in class

constitutiveEq and can be found in ﬁle 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 ﬂuid ﬂows.

Note that the second term is included to account for a shear-rate dependent

viscosity coeﬃcient. By deﬁnition, a GNF ﬂuid has no elasticity, thus τ=0.

For a viscoelastic ﬂuid (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 ∇·τ=∇·τ−∇·ηs∇u; if BSD, then the both-sides-diﬀusion

technique is used and ∇·τ=∇·τ− ∇·ηp∇u+∇·(ηs+ηp)∇u(Eq. 3.5);

otherwise (if coupling), the stress-velocity coupling technique of Eq. (3.17)

is used and ∇·τ=∇·τ− ∇·ηp∇u+∇·(η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 artiﬁcial

diﬀusion may mask the real phenomena in transient simulations. For the

cases using stabilization, the explicit behavior eﬀects 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 ﬁle consitutiveEq.C, a

function divTauS() is also included, which retrieves part of the extra-stress

contribution to the momentum equation, when solving two-phase ﬂows (this

topic will be discussed later).

1#include "Oldroyd_BLog.H"

#include "addToRunTimeSelectionTable.H"

3

// **************Static Data Members ********

*****//

5

namespace Foam{

7namespace constitutiveEqs{

defineTypeNameAndDebug(Oldroyd_BLog, 0);

9addToRunTimeSelectionTable(constitutiveEq, Oldroyd_BLog,

dictionary);

}

11 }

// ****************Constructors *********

*****//

13

Foam::constitutiveEqs::Oldroyd_BLog::Oldroyd_BLog

15 (

const word& name,

17 const volVectorField& U,

const surfaceScalarField& phi,

19 const dictionary& dict

)

21 :

constitutiveEq(name, U, phi),

23 tau_

(

25 IOobject

(

27 "tau" + name,

U.time().timeName(),

29 U.mesh(),

IOobject::MUST_READ,

31 IOobject::AUTO_WRITE

),

33 U.mesh()

),

35 theta_

(

37 IOobject

(

39 "theta" + name,

U.time().timeName(),

CHAPTER 4. Overview of rheoTool 46

41 U.mesh(),

IOobject::MUST_READ,

43 IOobject::AUTO_WRITE

),

45 U.mesh()

),

47 eigVals_

(

49 IOobject

(

51 "eigVals" + name,

U.time().timeName(),

53 U.mesh(),

IOobject::READ_IF_PRESENT,

55 IOobject::AUTO_WRITE

),

57 U.mesh(),

dimensionedTensor

59 (

"I",

61 dimless,

pTraits<tensor>::I

63 ),

zeroGradientFvPatchField<tensor>::typeName

65 ),

eigVecs_

67 (

IOobject

69 (

"eigVecs" + name,

71 U.time().timeName(),

U.mesh(),

73 IOobject::READ_IF_PRESENT,

IOobject::AUTO_WRITE

75 ),

U.mesh(),

77 dimensionedTensor

(

79 "I",

dimless,

81 pTraits<tensor>::I

),

83 zeroGradientFvPatchField<tensor>::typeName

),

85 rho_(dict.lookup("rho")),

etaS_(dict.lookup("etaS")),

87 etaP_(dict.lookup("etaP")),

lambda_(dict.lookup("lambda")),

89 {

checkForStab(dict);

91 }

93 // ***************Member Functions ********

*****//

void Foam::constitutiveEqs::Oldroyd_BLog::correct()

95 {

CHAPTER 4. Overview of rheoTool 47

// Decompose grad(U).T()

97

volTensorField L = fvc::grad(U());

99

dimensionedScalar c1( "zero", dimensionSet(0, 0, -1, 0, 0, 0,

0), 0.);

101 volTensorField B = c1 *eigVecs_;

volTensorField omega = B;

103 volTensorField M = (eigVecs_.T() & L.T() & eigVecs_);

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

);

115

fvSymmTensorMatrix thetaEqn

117 (

fvm::ddt(theta_)

119 + fvm::div(phi(), theta_)

==

121 symm

(

123 (omega&theta_)

- (theta_&omega)

125 + 2.0 *B

+ (1.0/lambda_)

127 *(

eigVecs_ &

129 (

inv(eigVals_)

131 - Itensor

)

133 & eigVecs_.T()

)

135 )

137 );

139 thetaEqn.relax();

thetaEqn.solve();

141

// Diagonalization of theta

143

calcEig(theta_, eigVals_, eigVecs_);

145

// Convert from theta to tau

147

tau_ = (etaP_/lambda_) *symm( (eigVecs_ & eigVals_ & eigVecs_

.T()) - Itensor);

149

CHAPTER 4. Overview of rheoTool 48

tau_.correctBoundaryConditions();

151 }

Listing 4.1: Source code for the Oldroyd-BLog constitutive model

(Oldroyd BLog.C)

1tmp<fvVectorMatrix> constitutiveEq::divTau

(

3const volVectorField& U

)const

5{

7if (isGNF())

{

9return

(

11 fvm::laplacian( eta()/rho(), U, "laplacian(eta,U)")

+ (fvc::grad(U) & fvc::grad(eta()/rho()))

13 );

}

15 else

{

17 if (stabMeth_ == 0) // none

{

19

return

21 (

fvc::div(tau()/rho(), "div(tau)")

23 + fvm::laplacian(etaS()/rho(), U, "laplacian(eta,U

)")

);

25

}

27 else if (stabMeth_ == 1) // BSD

{

29

return

31 (

fvc::div(tau()/rho(), "div(tau)")

33 - fvc::laplacian(etaP()/rho(), U, "laplacian(eta,U

)")

+ fvm::laplacian( (etaP()+ etaS())/rho(), U, "

laplacian(eta,U)")

35 );

37 }

else // coupling

39 {

41 return

(

43 fvc::div(tau()/rho(), "div(tau)")

- (etaP()/rho()) *fvc::div(fvc::grad(U))

45 + fvm::laplacian( (etaP() + etaS())/rho(), U, "

laplacian(eta,U)")

);

CHAPTER 4. Overview of rheoTool 49

47

}

49 }

}

Listing 4.2: Source code of the virtual function divTau() deﬁned in ﬁle

constitutiveEq.C

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 oﬀer 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 ﬁle 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 modiﬁcations, all the log-conformation-based models will be

aﬀected 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 ﬁles inside it (.C and .H ﬁles) with the new model’s name. Remove the

.dep ﬁle 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 ﬁles, ﬁnd-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 ﬁle.

However, it is a good idea to always check ﬁrst 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 ﬁle src/libs/constitutiveEquations/Make/files

by adding the path for the source code (see the entries for the other models

already there).

•make a ﬁrst 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 ﬁle, 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 ﬁle, 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 deﬁne functions divTau() and divTauS() for your new model, if the

ones deﬁned and used by default in the base class (see constitutiveEq.C)

are not adequate. After all the changes on the code had been completed,

compile again by running script Allwmake.

If all the steps listed above were successfully executed, then the new model

is now available to all the solvers of rheoTool. In order to use library constitu-

tiveEquations in any solver other than the ones provided with rheoTool, the user

should:

•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 ﬂows 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

Poisson-Nernst-Planck

(PNP)

6NernstPlanck

6NernstPlanckCoupled −F

N

P

i=1

zici(∇Ψ−Ea)

∇·(ε∇φExt) = 0

∇·(ε∇ψ) = −F

N

P

i=1

zici

∂ci

∂t +u·∇ci=∇·(Di∇ci) + ∇·Diezi

kT ∇Ψci

Section 3.7.1

Poisson-Boltzmann

(PB) PoissonBoltzmann −F

N

P

i=1

zici,0exp −ezi

kT ψ(∇Ψ−Ea)

∇·(ε∇φExt) = 0

∇·(ε∇ψ) + ψF

N

P

i=1

(aibi)∗=−F

N

P

i=1

(ai)∗+ψ∗F

N

P

i=1

(aibi)∗

with bi=−ezi

kT and ai=zici,0exp (biψ)

Section 3.7.3

Debye-H¨uckel

(DH) DebyeHuckel −F

N

P

i=1

zici,01−ezi

kT ψ(∇Ψ−Ea)

∇·(ε∇φExt) = 0

∇·(ε∇ψ) = −F

N

P

i=1

zici,01−ezi

kT ψSection 3.7.4

Slip velocity slipSmoluchowski 0∇·(ε∇φExt) = 0 Section 3.7.5

Ohmic Ohmic [∇·(ε∇φExt)] (∇φExt −Ea)

∂σ

∂t +u·∇σ=Deﬀ∇2σ

∇·(σ∇φ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 ﬁeld.

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.

51

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 deﬁned, 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 ﬁeld deﬁnition of the

body-force entering the momentum equation. The choice is through the variable

psiContrib, which should be deﬁned in a dictionary, as explained later in Section

5.4.1.

All the PNP, PB and DH models support multi-species modeling, with an ar-

bitrary number of species, each having diﬀerent properties (charge valence, and

diﬀusivity, 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

diﬀerent diﬀusion coeﬃcients.

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

versions.

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 diﬀerences, at the code level, between the diﬀerent 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 ﬁeld (ci), the charge valence (zi)

and the diﬀusivity (Di). The NernstPlanck class may have Ninstances of

the NPSpecie subclass, as much as deﬁned by the user.

•lines 35-107: this is the constructor of the main class, where the several

ﬁelds 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 ﬁeld 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 ﬁeld in lines 122-155 (compare

with Eq. 3.26 and Table 4.2). The computation of the electric ﬁeld 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 ﬁeld, which can be optionally deﬁned 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 speciﬁed 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 deﬁned 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 :

ci_

20 (

IOobject

22 (

name,

24 phi.time().timeName(),

phi.mesh(),

26 IOobject::MUST_READ,

IOobject::AUTO_WRITE

28 ),

phi.mesh()

30 ),

zi_(dict.lookup("z")),

32 Di_(dict.lookup("D"))

{}

34

// **********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)),

psi_

46 (

IOobject

48 (

"psi" + name,

50 phi.time().timeName(),

phi.mesh(),

52 IOobject::MUST_READ,

IOobject::AUTO_WRITE

54 ),

phi.mesh()

56 ),

phiE_

CHAPTER 4. Overview of rheoTool 56

58 (

IOobject

60 (

"phiE" + name,

62 phi.time().timeName(),

phi.mesh(),

64 IOobject::READ_IF_PRESENT,

solvePhiE_ == false ? (IOobject::NO_WRITE) : (IOobject::

AUTO_WRITE)

66 ),

phi.mesh(),

68 dimensionedScalar

(

70 "zero",

psi_.dimensions(),

72 pTraits<scalar>::zero

),

74 zeroGradientFvPatchField<scalar>::typeName

),

76 relPerm_(dict.lookup("relPerm")),

T_(dict.lookup("T")),

78 extraE_(dict.lookupOrDefault<dimensionedVector>("extraEField",

dimensionedVector("0", dimensionSet(1, 1, -3, 0, 0, -1, 0),

vector::zero))),

psiContrib_(dict.lookupOrDefault<bool>("psiContrib", true)),

80 phiEEqRes_(phi.mesh().solutionDict().subDict("electricControls").

subDict("phiEEqn").lookupOrDefault<scalar>("residuals", 1e-7)),

psiEqRes_(phi.mesh().solutionDict().subDict("electricControls").

subDict("psiEqn").lookupOrDefault<scalar>("residuals", 1e-7)),

82 ciEqRes_(phi.mesh().solutionDict().subDict("electricControls").

subDict("ciEqn").lookupOrDefault<scalar>("residuals", 1e-7)),

maxIterPhiE_(phi.mesh().solutionDict().subDict("electricControls").

subDict("phiEEqn").lookupOrDefault<scalar>("maxIter", 50)),

84 maxIterPsi_(phi.mesh().solutionDict().subDict("electricControls").

subDict("psiEqn").lookupOrDefault<scalar>("maxIter", 50)),

maxIterCi_(phi.mesh().solutionDict().subDict("electricControls").

subDict("ciEqn").lookupOrDefault<scalar>("maxIter", 50)),

86 nIterPNP_(phi.mesh().solutionDict().subDict("electricControls").

lookupOrDefault<int>("nIterPNP", 2)),

species_(),

88 nSpecies_(0)

{

90 PtrList<entry> specEntries(dict.lookup("species"));

nSpecies_ = specEntries.size();

92 species_.setSize(nSpecies_);

94 forAll (species_, specI)

{

96 species_.set

(

98 specI,

new NPSpecie

100 (

specEntries[specI].keyword(),

102 phi,

specEntries[specI].dict()

104 )

CHAPTER 4. Overview of rheoTool 57

);

106 }

}

108

// ******** Member Functions * *********//

110 Foam::tmp<Foam::volVectorField> Foam::EDFEquations::NernstPlanck::Fe()

const

{

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 }

else

132 {

return

134 (

-rhoE *( fvc::grad(phiE_) - extraE_)

136 );

}

138 }

else

140 {

if (psiContrib_)

142 {

return

144 (

-rhoE *( fvc::grad(psi_) - extraE_)

146 );

}

148 else

{

150 return

(

152 -rhoE *(-extraE_)

);

154 }

}

156 }

158 void Foam::EDFEquations::NernstPlanck::correct()

{

CHAPTER 4. Overview of rheoTool 58

160

// Electrokinetic coupling loop

162 for (int j=0; j<nIterPNP_; j++)

{

164

Info << "PNP Coupling iteration: " << j << endl;

166

scalar res=GREAT;

168 scalar iter=0;

170 //- Equation for the external potential (loop for the case

// of non-orthogonal grids)

172 if (solvePhiE_)

{

174 while (res > phiEEqRes_ && iter < maxIterPhiE_)

{

176 fvScalarMatrix phiEEqn

(

178 fvm::laplacian(phiE_)

);

180

phiEEqn.relax();

182 res=phiEEqn.solve().initialResidual();

184 iter++;

}

186 }

188 //- Equation for the intrinsic potential

190 res=GREAT;

iter=0;

192

volScalarField souE(psi_ *dimensionedScalar("norm1",dimless/dimArea

,0.));

194

forAll (species_, i)

196 {

souE +=

198 (

-species_[i].zi()*species_[i].ci()*FK_

200 /(relPerm_*epsilonK_)

);

202 }

204 while (res > psiEqRes_ && iter < maxIterPsi_)

{

206

fvScalarMatrix psiEqn

208 (

fvm::laplacian(psi_)

210 ==

souE

212 );

214 psiEqn.relax();

res=psiEqn.solve().initialResidual();

CHAPTER 4. Overview of rheoTool 59

216

iter++;

218 }

220 //- Nernst-Planck equation for each ionic specie

222 forAll (species_, i)

{

224 res=GREAT;

iter=0;

226

volScalarField& ci = species_[i].ci();

228

dimensionedScalar cf(species_[i].Di() *eK_ *species_[i].zi() / (kbK_

*T_) );

230

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)")

);

242

ciEqn.relax(phi().mesh().equationRelaxationFactor("ci"));

244 res=ciEqn.solve(phi().mesh().solver("ci")).initialResidual();

246 ci = Foam::max( dimensionedScalar("lowerLimit",ci.dimensions(), 0.),

ci );

iter++;

248

}

250 }

}

252 }

Listing 4.3: Source code of the Poisson-Nernst-Planck model in ﬁle

NernstPlanck.C.

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 ﬁles inside it (.C and .H ﬁles) with the new model’s name. Remove the

.dep ﬁle 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 ﬁles, ﬁnd-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 ﬁle.

However, it is a good idea to always check ﬁrst 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 ﬁle src/libs/EDFModels/Make/files by adding the

path for the source code (see the entries for the other models already there).

•make a ﬁrst 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 ﬁles 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

versions.

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 modiﬁed 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/

libs/brownianDynamics.

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 4Mﬁles (Mbeing the number

of molecules) for each time-step saved, which would overburden any ﬁle trans-

fer operation, notwithstanding the small disk space used by each ﬁle. 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 classiﬁed object

Cloud as blind, because it does not diﬀerentiate 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 ﬁles 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 eﬃcient

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 ﬁnd 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

diﬀerence comparing to mx is in the data organization. While we can use mx to

easily ﬁnd 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 ﬁrst needs to understand the interplay between the

local ﬁelds of sPCloudInterface and the particles ﬁelds 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 conﬁguration of the molecules each

time-step. In what follows, we brieﬂy discuss the structure of this function (Listing

4.4):

•line 3: the beads positions at the beginning of a time-step (mx ) are copied

to mx0 . The algorithm structure ensures that the beads positions in mx0

are such that all the springs of the active molecules are shorter than l, for

bounded models (Section 3.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 diﬀusion 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 deﬁned inside sPClou

dInterface.C.

•lines 10-11: if exclusion volume forces are activated by the user, then their

contribution to the beads velocity is added through a call to function fEV(),

deﬁned inside sPCloudInterface.C.

•line 13: function sendU() copies the beads velocity from the local mU

P trList <> to the particles ﬁeld U. This function is deﬁned inside sPClou

dInterface.C.

•line 16: function moveAndReceive() comprises two main steps. In the ﬁrst

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 ﬁle solidParticle.C). In the second step,

the function updates mx with the ﬁnal 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

updated.

•line 23: based on the current beads positions (mx ), the spring vectors are

computed and a check is carried out for Ri< αl (see Section 3.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 deﬁned 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 ﬁxed vector. In this case, we add the

corresponding velocity to mU , such that the translation can be eﬀective 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 ﬁrst 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 ﬁeld 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-oﬀ errors, but ﬁnal results

are always dimensional. Round-oﬀ 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 diﬀerent scales involved.

bool sPCloudInterface::update()

2{

mx0_ = mx_;

4

if (isHI_)

6computeDSigma();

8fBrownian();

10 if (isExclusionVolumeF_)

fEV();

12

sendU();

14

// Move particles: Drag + Brownian + Exclusion Volume

16 moveAndReceive(true);

18 // Compute Spring force term explicitly and update local mU and

mx.

mxStar_ = mx_;

20 spModel_->fSpring();

22 // Check for violations in spring lengths and add spring force

implicitly if needed

spModel_->checkSpringsLength(mxStar_, mx0_);

24

// Push back the molecules if not tethered and if pushback

vector is not negligible

26 if (!isTethered_ && mag(pBackV_)>SMALL)

pushToX0();

28

sendU();

30

// Move particles: spring force only

32 moveAndReceive(false);

34 // Write data (controlled by output time)

writeM();

36

// Write statistics (has its own control for output)

38 if (writeStats_)

CHAPTER 4. Overview of rheoTool 65

writeStatistics();

40

if (nMolc_>0)

42 {

return true;

44 }

else

46 {

return false;

48 }

}

Listing 4.4: Member function update() of class sPCloudInterface. The source

code can be found in ﬁle 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 deﬁned analytically or computed numerically.

In case it is deﬁned analytically, tensor ∇u, provided by the user, deﬁnes 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 unconﬁned ﬂow

is intended, but a conﬁned ﬂow can also be imposed by shortening one or multiple

directions. Importantly, the boundaries used to impose conﬁned ﬂow conditions

must be of wall type (patch types will be crossed by the molecules, and remaining

types give undeﬁned behavior).

In case it is computed numerically, then two types are still available in rheoTool:

hydrodynamic forcing, in which case ufis the ﬂuid velocity; electric forcing for

polyelectrolytes immersed in an electric ﬁeld, in which case uf=µEis the elec-

trophoretic velocity (µis the electrophoretic mobility and Eis the electric ﬁeld).

Both types of forcing can coexist, and each one requires that the respective con-

tinuum ﬁeld is available in the case directory.

Note that 2D meshes, and the corresponding continuum ﬁelds, 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 eﬀect of a planar velocity ﬁeld 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 ﬁeld deﬁned 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-deﬁned.

•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 deﬁne a barycentric coordinate system for each tetrahe-

dron, such that any point inside it can be uniquely identiﬁed by a quadruplet

(b1b2b3b4), with

4

P

i=1

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 ﬁeld computed at the cell-centers

(OpenFOAM R

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 ﬁeld 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+ (xP−xC)·∇uC. The numerical scheme used

to compute the cell-centered gradient (∇uC) can be deﬁned 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

ﬁgure corresponds to a fully-developed ﬁeld in the x-direction, which displays a

maximum in the y-direction, for y= 0, and is symmetric about y= 0. The proﬁle

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 ﬂow

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 ﬁeld of quadrangular pyramids

(apex at cell centers), with the height of each vertex corresponding to the local

forcing value.

For a proﬁle taken over y= 0, method Gradient approaches correctly the

theoretical function (constant over the x-direction), whereas method Barycen-

tricWeights retrieves a sawtooth proﬁle. Consider now that the molecules’ center

CHAPTER 4. Overview of rheoTool 67

of mass is artiﬁcially ﬁxed (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 speciﬁc abscissa selected, as one would expected for a fully-developed ﬁeld

in the x-direction. On the other hand, if the strain-rate in the sawtooth proﬁle

retrieved by method BarycentricWeights is such that W i =λ˙γ >> 1, then the

results will strongly depend on the speciﬁc abscissa selected.

For a proﬁle taken over x= 0, method Gradient is unable to capture any varia-

tion of the ﬁeld 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 ﬁeld 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 artiﬁcial one in the x-direction).

3

33

1

1

1

1

1

1

2 2

2 2

y =1

y =3

y =-1

y =-3

x =-3 x =-1 x =1x =3

x

1

-1

uf3

Profile at

y = 0

y

1

-1

3

2

Profile at

x = 0

uf

2

Gradient

Barycentric

Weights

Gradient

Barycentric

Weights

Figure 4.1: Hypothetical 2D external forcing ﬁeld (a scalar ﬁeld to simplify) in a

co-located grid (values at cell centers), aiming to represent a fully-developed ﬁeld

in the x-direction. The green values deﬁned at vertices are obtained by inverse

distance weighting of the cell-centered data. On the right, the proﬁles of uftaken

at x= 0 and y= 0 are represented for two diﬀerent 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 speciﬁc case presented above, the issues with the two methods could have been

reduced by reﬁning 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 reﬁned 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 brieﬂy discussed in Section 3.8.3.

Table 4.3: Available spring force models for Brownian dynamics simulations.

Model (1)

TypeName (2)FS

i(3)fk

i(4,5)Jk

ij, a (j6=i, j∈gi)

Hookean Hookean H

N

P

j∈gi

xij – –

Marko-Siggia MarkoSiggia 2

3Hl

N

P

j∈gih|xij |

l−1

4+1

4(1−|xij |/l)2ixij

|xij |rk

i−r∗

i−∆tD

kT FS,k

i2

3

H0

dh 1

4(d−1)2+d−1

4δ2

d2−1−δ2

d1−1

2(d−1)3i

Cohen Pad´e CohenPade H

3

N

P

j∈gih3−(|xij |/l)2

1−(|xij |/l)2ixij rk

i−r∗

i−∆tD

kT FS,k

iH0

3

2δ2

(d2−1) hd2−3

d2−1−1−d2−3

2δ2i

FENE FENE H

N

P

j∈gi

1

1−(|xij |/l)2xij rk

i−r∗

i−∆tD

kT FS,k

iH0h1

d2−1−2δ2

(d2−1)2i

(1) Corresponds to the name entry identifying the model in the source code.

(2) H=3kT Nk,s

l2.

(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 oﬀ-diagonal ij (j6=i, j ∈gi) elements of the Jacobian matrix of the Cartesian component a=x, y, z of fk. All the

oﬀ-diagonal elements for which j /∈giare zero. The diagonal elements are obtained from the sum of the oﬀ-diagonal elements, Jk

ii,a = 1 −

N

P

m=1,m6=i

Jk

im,a. This

symmetric matrix is used in the Newton-Raphson method (see Section 3.8.4).

(5) H0=H∆tD

kT ,d=|xk

ij |/l =|rk

j−rk

i|/l,δ= (rk

j,a −rk

i,a)/l

69

CHAPTER 4. Overview of rheoTool 70

In general, the FENE and Cohen Pad´e models are stiﬀer 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 speciﬁed in Table 4.3 for all bounded

models.

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

BDmolecule:

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 oﬀ-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 ﬁxing the molecules center of mass

rheoTool allows the simulation of tethered molecules, i.e., molecules with some

speciﬁc part of the chain ﬁxed in space. The only restriction is that solely the

ﬁrst bead of a given molecule (the one having index = 0 inside the chain) can be

ﬁxed. 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 diﬀerent from tethering a molecule,

and the use of both in simultaneous is not allowed (would it make sense?). When

the option of ﬁxing 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 ﬂow. The value of n, the number of consecutive time-steps in which

the molecules movement is unconstrained, is deﬁned 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 artiﬁcial repulsion aims to reproduce the ﬁnite size of the beads, that

would never let them to approach inﬁnitely 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 ﬂow 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 undeﬁned 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 diﬀerent

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 conﬁned

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 ﬁles are written to this directory:

IDs.txt: contains three columns, where the ﬁrst 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 ﬁle 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 ﬁrst column of this ﬁle 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 ﬁrst column for time). We deﬁne stretch as

the maximum inter-bead distance in a molecule: stretch = max(|rj−ri|),

i, j ={1..N}. If a given molecule is deleted at some point of the simulation,

then the column corresponding to that molecule will be ﬁlled with zeros

starting from the row corresponding to that time. Since the length of a

molecule can never b